class FEFaceValues;
template <int dim, int spacedim>
class FESubfaceValues;
+namespace NonMatching
+{
+ template <int dim>
+ class FEImmersedSurfaceValues;
+}
template <int dim, int spacedim>
class FESystem;
friend class FEValues<dim, spacedim>;
friend class FEFaceValues<dim, spacedim>;
friend class FESubfaceValues<dim, spacedim>;
+ friend class NonMatching::FEImmersedSurfaceValues<dim>;
friend class FESystem<dim, spacedim>;
// explicitly check for sensible template arguments, but not on windows
#include <deal.II/hp/q_collection.h>
+#include <deal.II/non_matching/immersed_surface_quadrature.h>
+
#include <array>
#include <cmath>
#include <memory>
class FEFaceValues;
template <int dim, int spacedim>
class FESubfaceValues;
+namespace NonMatching
+{
+ template <int dim>
+ class FEImmersedSurfaceValues;
+}
/**
dealii::internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
&output_data) const = 0;
+ /**
+ * The equivalent of Mapping::fill_fe_values(), but for the case that the
+ * quadrature is an ImmersedSurfaceQuadrature. See there for a comprehensive
+ * description of the input parameters. This function is called by
+ * FEImmersedSurfaceValues::reinit().
+ */
+ virtual void
+ fill_fe_immersed_surface_values(
+ const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+ const NonMatching::ImmersedSurfaceQuadrature<dim> & quadrature,
+ const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
+ dealii::internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
+ &output_data) const;
+
/**
* @}
*/
friend class FEValues<dim, spacedim>;
friend class FEFaceValues<dim, spacedim>;
friend class FESubfaceValues<dim, spacedim>;
+ friend class NonMatching::FEImmersedSurfaceValues<dim>;
};
internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
&output_data) const override;
+ // documentation can be found in Mapping::fill_fe_immersed_surface_values()
+ virtual void
+ fill_fe_immersed_surface_values(
+ const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+ const NonMatching::ImmersedSurfaceQuadrature<dim> & quadrature,
+ const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
+ dealii::internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
+ &output_data) const override;
+
/**
* @}
*/
dealii::internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
&output_data) const override;
+ // documentation can be found in Mapping::fill_fe_immersed_surface_values()
+ virtual void
+ fill_fe_immersed_surface_values(
+ const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+ const NonMatching::ImmersedSurfaceQuadrature<dim> & quadrature,
+ const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
+ dealii::internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
+ &output_data) const override;
+
/**
* @}
*/
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2021 - 2021 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+#ifndef dealii_non_matching_fe_immersed_values_h
+#define dealii_non_matching_fe_immersed_values_h
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+
+#include <deal.II/fe/fe.h>
+#include <deal.II/fe/fe_update_flags.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping.h>
+
+#include <deal.II/grid/tria_iterator.h>
+
+#include <deal.II/non_matching/immersed_surface_quadrature.h>
+
+DEAL_II_NAMESPACE_OPEN
+namespace NonMatching
+{
+ /**
+ * Finite element evaluated in the quadrature points of an
+ * ImmersedSurfaceQuadrature of a cell.
+ *
+ * The shape functions values and their derivatives are the same as for an
+ * FEValues-object, but the JxW-values are computed with the transformation
+ * described in the documentation of ImmersedSurfaceQuadrature. Further, the
+ * normal_vector-function returns the normal to the immersed surface.
+ *
+ * The reinit-function of this class exist mostly to be consistent with the
+ * other FEValues-like classes. The immersed quadrature rule will typically
+ * vary between each cell of the triangulation. Thus, an
+ * FEImmersedSurfaceValues object can, typically, not be reused for different
+ * cells.
+ *
+ * See also documentation in FEValuesBase.
+ *
+ * @ingroup feaccess
+ */
+ template <int dim>
+ class FEImmersedSurfaceValues : public FEValuesBase<dim, dim>
+ {
+ public:
+ /**
+ * Constructor. Gets cell-independent data from mapping and finite element
+ * objects, matching the quadrature rule and update flags.
+ *
+ * @note Currently this class is only implemented for MappingCartesian.
+ */
+ FEImmersedSurfaceValues(const Mapping<dim> & mapping,
+ const FiniteElement<dim> & element,
+ const ImmersedSurfaceQuadrature<dim> &quadrature,
+ const UpdateFlags update_flags);
+
+ /**
+ * Reinitialize quantities (normals, JxW-values, etc) for the given cell of
+ * type "iterator into a Triangulation object".
+ */
+ void
+ reinit(const typename Triangulation<dim>::cell_iterator &cell);
+
+ /**
+ * Reinitialize quantitites (shape function values, gradients, etc) for the
+ * given cell of type "iterator into a DoFHandler object", and the finite
+ * element associated with this object.
+ */
+ template <bool level_dof_access>
+ void
+ reinit(
+ const TriaIterator<DoFCellAccessor<dim, dim, level_dof_access>> &cell);
+
+ /**
+ * Returns the surface gradient of the shape function with index
+ * @p function_no at the quadrature point with index @p quadrature_point.
+ *
+ * The surface gradient is defined as the projection of the gradient to the
+ * tangent plane of the surface:
+ * $ \nabla u - (n \cdot \nabla u) n $,
+ * where $n$ is the unit normal to the surface.
+ *
+ * @dealiiRequiresUpdateFlags{update_gradients | update_normal_vectors}
+ */
+ Tensor<1, dim>
+ shape_surface_grad(const unsigned int function_no,
+ const unsigned int quadrature_point) const;
+
+ /**
+ * Return one vector component of the surface gradient of the shape function
+ * at a quadrature point. See the definition of the surface gradient in the
+ * shape_surface_grad function.
+ *
+ * @p function_no Index of the shape function to be evaluated.
+ *
+ * @p point_no Index of the quadrature point at which the function is to be
+ * evaluated.
+ *
+ * @p component Vector component to be evaluated.
+ *
+ * @dealiiRequiresUpdateFlags{update_gradients | update_normal_vectors}
+ */
+ Tensor<1, dim>
+ shape_surface_grad_component(const unsigned int function_no,
+ const unsigned int quadrature_point,
+ const unsigned int component) const;
+
+ /**
+ * Return a reference to the copy of the quadrature formula stored by this
+ * object.
+ */
+ const NonMatching::ImmersedSurfaceQuadrature<dim> &
+ get_quadrature() const;
+
+ protected:
+ /**
+ * Do work common to the constructors.
+ */
+ void
+ initialize(const UpdateFlags update_flags);
+
+ /**
+ * The reinit() functions do only that part of the work that requires
+ * knowledge of the type of iterator. After setting present_cell(), they
+ * pass on to this function, which does the real work, and which is
+ * independent of the actual type of the cell iterator.
+ */
+ void
+ do_reinit();
+
+ /**
+ * Copy of the quadrature formula that was passed to the constructor.
+ */
+ const ImmersedSurfaceQuadrature<dim> quadrature;
+ };
+
+} // namespace NonMatching
+DEAL_II_NAMESPACE_CLOSE
+
+#endif /* dealii_non_matching_fe_immersed_values_h */
+template <int dim, int spacedim>
+void
+Mapping<dim, spacedim>::fill_fe_immersed_surface_values(
+ const typename Triangulation<dim, spacedim>::cell_iterator &,
+ const NonMatching::ImmersedSurfaceQuadrature<dim> &,
+ const typename Mapping<dim, spacedim>::InternalDataBase &,
+ dealii::internal::FEValuesImplementation::MappingRelatedData<dim, spacedim> &)
+ const
+{
+ AssertThrow(false, ExcNotImplemented());
+}
+
+
+
template <int dim, int spacedim>
void
Mapping<dim, spacedim>::transform_points_real_to_unit_cell(
+template <int dim, int spacedim>
+void
+MappingCartesian<dim, spacedim>::fill_fe_immersed_surface_values(
+ const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+ const NonMatching::ImmersedSurfaceQuadrature<dim> & quadrature,
+ const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
+ dealii::internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
+ &output_data) const
+{
+ AssertDimension(dim, spacedim);
+
+ // Convert data object to internal data for this class. Fails with an
+ // exception if that is not possible.
+ Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
+ ExcInternalError());
+ const InternalData &data = static_cast<const InternalData &>(internal_data);
+
+
+ update_cell_extents(cell, CellSimilarity::none, data);
+
+ maybe_update_cell_quadrature_points(cell,
+ data,
+ output_data.quadrature_points);
+
+ if (data.update_each & update_normal_vectors)
+ for (unsigned int i = 0; i < output_data.normal_vectors.size(); ++i)
+ {
+ // The normals are n = J^{-T} * \hat{n} before normalizing.
+ Tensor<1, dim> normal;
+ const Tensor<1, dim> &ref_space_normal = quadrature.normal_vector(i);
+ for (unsigned int d = 0; d < dim; ++d)
+ {
+ normal[d] = ref_space_normal[d] / data.cell_extents[d];
+ }
+ normal /= normal.norm();
+ output_data.normal_vectors[i] = normal;
+ }
+
+ if (data.update_each & update_JxW_values)
+ for (unsigned int i = 0; i < output_data.JxW_values.size(); ++i)
+ {
+ const Tensor<1, dim> &ref_space_normal = quadrature.normal_vector(i);
+
+ // J^{-T} \times \hat{n}
+ Tensor<1, dim> invJTxNormal;
+ double det_jacobian = 1.;
+ for (unsigned int d = 0; d < dim; ++d)
+ {
+ det_jacobian *= data.cell_extents[d];
+ invJTxNormal[d] = ref_space_normal[d] / data.cell_extents[d];
+ }
+ output_data.JxW_values[i] =
+ det_jacobian * invJTxNormal.norm() * quadrature.weight(i);
+ }
+
+ maybe_update_volume_elements(data);
+ maybe_update_jacobians(data, CellSimilarity::none, output_data);
+ maybe_update_jacobian_derivatives(data, CellSimilarity::none, output_data);
+ maybe_update_inverse_jacobians(data, CellSimilarity::none, output_data);
+}
+
+
+
template <int dim, int spacedim>
void
MappingCartesian<dim, spacedim>::transform(
+template <int dim, int spacedim>
+void
+MappingQ<dim, spacedim>::fill_fe_immersed_surface_values(
+ const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+ const NonMatching::ImmersedSurfaceQuadrature<dim> & quadrature,
+ const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
+ dealii::internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
+ &output_data) const
+{
+ Assert(dim == spacedim, ExcNotImplemented());
+
+ // ensure that the following static_cast is really correct:
+ Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
+ ExcInternalError());
+ const InternalData &data = static_cast<const InternalData &>(internal_data);
+
+ const unsigned int n_q_points = quadrature.size();
+
+ data.mapping_support_points = this->compute_mapping_support_points(cell);
+ data.cell_of_current_support_points = cell;
+
+ internal::MappingQImplementation::maybe_compute_q_points<dim, spacedim>(
+ QProjector<dim>::DataSetDescriptor::cell(),
+ data,
+ output_data.quadrature_points);
+
+ internal::MappingQImplementation::maybe_update_Jacobians<dim, spacedim>(
+ CellSimilarity::none, QProjector<dim>::DataSetDescriptor::cell(), data);
+
+ internal::MappingQImplementation::maybe_update_jacobian_grads<dim, spacedim>(
+ CellSimilarity::none,
+ QProjector<dim>::DataSetDescriptor::cell(),
+ data,
+ output_data.jacobian_grads);
+
+ internal::MappingQImplementation::maybe_update_jacobian_pushed_forward_grads<
+ dim,
+ spacedim>(CellSimilarity::none,
+ QProjector<dim>::DataSetDescriptor::cell(),
+ data,
+ output_data.jacobian_pushed_forward_grads);
+
+ internal::MappingQImplementation::maybe_update_jacobian_2nd_derivatives<
+ dim,
+ spacedim>(CellSimilarity::none,
+ QProjector<dim>::DataSetDescriptor::cell(),
+ data,
+ output_data.jacobian_2nd_derivatives);
+
+ internal::MappingQImplementation::
+ maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
+ CellSimilarity::none,
+ QProjector<dim>::DataSetDescriptor::cell(),
+ data,
+ output_data.jacobian_pushed_forward_2nd_derivatives);
+
+ internal::MappingQImplementation::maybe_update_jacobian_3rd_derivatives<
+ dim,
+ spacedim>(CellSimilarity::none,
+ QProjector<dim>::DataSetDescriptor::cell(),
+ data,
+ output_data.jacobian_3rd_derivatives);
+
+ internal::MappingQImplementation::
+ maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
+ CellSimilarity::none,
+ QProjector<dim>::DataSetDescriptor::cell(),
+ data,
+ output_data.jacobian_pushed_forward_3rd_derivatives);
+
+ const UpdateFlags update_flags = data.update_each;
+ const std::vector<double> &weights = quadrature.get_weights();
+
+ if (update_flags & (update_normal_vectors | update_JxW_values))
+ {
+ AssertDimension(output_data.JxW_values.size(), n_q_points);
+
+ Assert(!(update_flags & update_normal_vectors) ||
+ (output_data.normal_vectors.size() == n_q_points),
+ ExcDimensionMismatch(output_data.normal_vectors.size(),
+ n_q_points));
+
+
+ for (unsigned int point = 0; point < n_q_points; ++point)
+ {
+ const double det = data.contravariant[point].determinant();
+
+ // check for distorted cells.
+
+ // TODO: this allows for anisotropies of up to 1e6 in 3D and
+ // 1e12 in 2D. might want to find a finer
+ // (dimension-independent) criterion
+ Assert(det > 1e-12 * Utilities::fixed_power<dim>(
+ cell->diameter() / std::sqrt(double(dim))),
+ (typename Mapping<dim, spacedim>::ExcDistortedMappedCell(
+ cell->center(), det, point)));
+
+ // The normals are n = J^{-T} * \hat{n} before normalizing.
+ Tensor<1, spacedim> normal;
+ for (unsigned int d = 0; d < spacedim; d++)
+ normal[d] =
+ data.covariant[point][d] * quadrature.normal_vector(point);
+
+ output_data.JxW_values[point] = weights[point] * det * normal.norm();
+
+ if (update_flags & update_normal_vectors)
+ {
+ normal /= normal.norm();
+ output_data.normal_vectors[point] = normal;
+ }
+ }
+ }
+
+ // copy values from InternalData to vector given by reference
+ if (update_flags & update_jacobians)
+ {
+ AssertDimension(output_data.jacobians.size(), n_q_points);
+ for (unsigned int point = 0; point < n_q_points; ++point)
+ output_data.jacobians[point] = data.contravariant[point];
+ }
+
+ // copy values from InternalData to vector given by reference
+ if (update_flags & update_inverse_jacobians)
+ {
+ AssertDimension(output_data.inverse_jacobians.size(), n_q_points);
+ for (unsigned int point = 0; point < n_q_points; ++point)
+ output_data.inverse_jacobians[point] =
+ data.covariant[point].transpose();
+ }
+}
+
+
+
template <int dim, int spacedim>
void
MappingQ<dim, spacedim>::fill_mapping_data_for_generic_points(
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR})
SET(_src
+ fe_immersed_values.cc
mesh_classifier.cc
quadrature_generator.cc
coupling.cc
immersed_surface_quadrature.cc
)
-SET(_inst
+ SET(_inst
+ fe_immersed_values.inst.in
mesh_classifier.inst.in
quadrature_generator.inst.in
coupling.inst.in
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2021 - 2021 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/grid/tria_iterator.h>
+
+#include <deal.II/non_matching/fe_immersed_values.h>
+
+
+DEAL_II_NAMESPACE_OPEN
+namespace NonMatching
+{
+ template <int dim>
+ FEImmersedSurfaceValues<dim>::FEImmersedSurfaceValues(
+ const Mapping<dim> & mapping,
+ const FiniteElement<dim> & element,
+ const ImmersedSurfaceQuadrature<dim> &quadrature,
+ const UpdateFlags update_flags)
+ : FEValuesBase<dim, dim>(quadrature.size(),
+ element.dofs_per_cell,
+ update_default,
+ mapping,
+ element)
+ , quadrature(quadrature)
+ {
+ initialize(update_flags);
+ }
+
+
+
+ template <int dim>
+ void
+ FEImmersedSurfaceValues<dim>::reinit(
+ const typename Triangulation<dim>::cell_iterator &cell)
+ {
+ // Check that mapping and reference cell type are compatible:
+ Assert(this->get_mapping().is_compatible_with(cell->reference_cell()),
+ ExcMessage(
+ "You are trying to call FEValues::reinit() with a cell of type " +
+ cell->reference_cell().to_string() +
+ " with a Mapping that is not compatible with it."));
+
+ // No FE in this cell, so no assertion necessary here.
+ this->maybe_invalidate_previous_present_cell(cell);
+ this->present_cell = {cell};
+
+ // This was the part of the work that is dependent on the actual data type
+ // of the iterator. Now pass on to the function doing the real work.
+ do_reinit();
+ }
+
+
+
+ template <int dim>
+ template <bool level_dof_access>
+ void
+ FEImmersedSurfaceValues<dim>::reinit(
+ const TriaIterator<DoFCellAccessor<dim, dim, level_dof_access>> &cell)
+ {
+ // Check that mapping and reference cell type are compatible:
+ Assert(this->get_mapping().is_compatible_with(cell->reference_cell()),
+ ExcMessage(
+ "You are trying to call FEValues::reinit() with a cell of type " +
+ cell->reference_cell().to_string() +
+ " with a Mapping that is not compatible with it."));
+
+ // Assert that the finite elements passed to the constructor and used by the
+ // DoFHandler used by this cell, are the same
+ Assert(static_cast<const FiniteElementData<dim> &>(*this->fe) ==
+ static_cast<const FiniteElementData<dim> &>(cell->get_fe()),
+ (typename FEValuesBase<dim>::ExcFEDontMatch()));
+
+ this->maybe_invalidate_previous_present_cell(cell);
+ this->present_cell = {cell};
+
+ // This was the part of the work that is dependent on the actual data type
+ // of the iterator. Now pass on to the function doing the real work.
+ do_reinit();
+ }
+
+
+
+ template <int dim>
+ void
+ FEImmersedSurfaceValues<dim>::do_reinit()
+ {
+ // First call the mapping and let it generate the data specific to the
+ // mapping.
+ if (this->update_flags & update_mapping)
+ {
+ this->get_mapping().fill_fe_immersed_surface_values(
+ this->present_cell,
+ quadrature,
+ *this->mapping_data,
+ this->mapping_output);
+ }
+
+ // Call the finite element and, with the data already filled by the mapping,
+ // let it compute the data for the mapped shape function values, gradients
+ // etc.
+ this->get_fe().fill_fe_values(this->present_cell,
+ CellSimilarity::none,
+ this->quadrature,
+ this->get_mapping(),
+ *this->mapping_data,
+ this->mapping_output,
+ *this->fe_data,
+ this->finite_element_output);
+ }
+
+
+
+ template <int dim>
+ Tensor<1, dim>
+ FEImmersedSurfaceValues<dim>::shape_surface_grad(
+ const unsigned int function_no,
+ const unsigned int quadrature_point) const
+ {
+ const unsigned int component = 0;
+ return shape_surface_grad_component(function_no,
+ quadrature_point,
+ component);
+ }
+
+
+
+ template <int dim>
+ Tensor<1, dim>
+ FEImmersedSurfaceValues<dim>::shape_surface_grad_component(
+ const unsigned int function_no,
+ const unsigned int quadrature_point,
+ const unsigned int component) const
+ {
+ const Tensor<1, dim> gradient =
+ this->shape_grad_component(function_no, quadrature_point, component);
+ const Tensor<1, dim> &normal = this->normal_vector(quadrature_point);
+
+ return gradient - (normal * gradient) * normal;
+ }
+
+
+
+ template <int dim>
+ const ImmersedSurfaceQuadrature<dim> &
+ FEImmersedSurfaceValues<dim>::get_quadrature() const
+ {
+ return quadrature;
+ }
+
+
+
+ template <int dim>
+ inline void
+ FEImmersedSurfaceValues<dim>::initialize(const UpdateFlags update_flags)
+ {
+ UpdateFlags flags = this->compute_update_flags(update_flags);
+
+ if (flags & (update_JxW_values | update_normal_vectors))
+ flags |= update_covariant_transformation;
+
+ // Initialize the base classes.
+ if (flags & update_mapping)
+ this->mapping_output.initialize(this->n_quadrature_points, flags);
+ this->finite_element_output.initialize(this->n_quadrature_points,
+ *this->fe,
+ flags);
+
+ // Then get objects into which the FE and the Mapping can store
+ // intermediate data used across calls to reinit. We can do this in
+ // parallel.
+ Threads::Task<
+ std::unique_ptr<typename FiniteElement<dim, dim>::InternalDataBase>>
+ fe_get_data = Threads::new_task(&FiniteElement<dim, dim>::get_data,
+ *this->fe,
+ flags,
+ *this->mapping,
+ this->quadrature,
+ this->finite_element_output);
+
+ Threads::Task<std::unique_ptr<typename Mapping<dim>::InternalDataBase>>
+ mapping_get_data;
+ if (flags & update_mapping)
+ mapping_get_data = Threads::new_task(&Mapping<dim>::get_data,
+ *this->mapping,
+ flags,
+ this->quadrature);
+
+ this->update_flags = flags;
+
+ // Then collect answers from the two task above.
+ this->fe_data = std::move(fe_get_data.return_value());
+ if (flags & update_mapping)
+ this->mapping_data = std::move(mapping_get_data.return_value());
+ else
+ this->mapping_data =
+ std::make_unique<typename Mapping<dim>::InternalDataBase>();
+ }
+
+
+#ifndef DOXYGEN
+# include "fe_immersed_values.inst"
+#endif
+
+} // namespace NonMatching
+DEAL_II_NAMESPACE_CLOSE
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2021 - 2021 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+for (deal_II_dimension : DIMENSIONS)
+ {
+ template class FEImmersedSurfaceValues<deal_II_dimension>;
+ }
+
+for (deal_II_dimension : DIMENSIONS; lda : BOOL)
+ {
+ template void FEImmersedSurfaceValues<deal_II_dimension>::reinit(
+ const TriaIterator<
+ DoFCellAccessor<deal_II_dimension, deal_II_dimension, lda>> &);
+ }
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2021 - 2021 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+#include <deal.II/base/point.h>
+#include <deal.II/base/tensor.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/mapping_cartesian.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/hp/mapping_collection.h>
+
+#include <deal.II/non_matching/fe_immersed_values.h>
+#include <deal.II/non_matching/immersed_surface_quadrature.h>
+
+#include "../tests.h"
+
+
+using namespace dealii;
+using NonMatching::FEImmersedSurfaceValues;
+
+// Set up a triangulation with a single Cartesian cell and a finite element
+// space with FE_Q-elements. Test that we can construct an
+// FEImmersedSurfaceValues object and compute known values when using a
+// Cartesian Mapping.
+template <int dim>
+class Test
+{
+public:
+ Test(const Mapping<dim> &mapping);
+
+ void
+ run();
+
+private:
+ // Construct a triangulation with a single cartesian cell, with extents.
+ // [0,2] in 1D,
+ // [0,2]x[0,3] in 2D, and
+ // [0,2]x[0,3]x[0,4] in 3D.
+ // Use extends that are different in all directions so that it's easier to
+ // check that the various quantities are mapped correctly.
+ void
+ setup_single_cartesian_cell_triangulation();
+
+ // Make the ImmersedSurfaceQuadraturewith contain a single point in the center
+ // of the cell, with weight .5 and a normal different in all components.
+ void
+ setup_single_point_quadrature();
+
+ // Test that the constant n_quadrature_points corresponds to what we sent in.
+ void
+ test_n_quadrature_points();
+
+ // Test that we can use the get_quadrature function to get the quadrature
+ // we passed to the FEImmersedSurfaceValues constructor.
+ void
+ test_get_quadrature();
+
+ // Print the quadrature point in real space to make sure that it's mapped to
+ // real space correctly.
+ void
+ test_point_mapped_correctly();
+
+ // Print the normal in real space to make sure that it's mapped to real space
+ // correctly.
+ void
+ test_normal();
+
+ // Print JxW to check that the value is correct.
+ void
+ test_JxW();
+
+ // Test that we can call the shape_surface_grad function.
+ void
+ test_shape_surface_grad();
+
+ const FE_Q<dim> element;
+ Triangulation<dim> triangulation;
+ DoFHandler<dim> dof_handler;
+ const SmartPointer<const Mapping<dim>> mapping;
+ NonMatching::ImmersedSurfaceQuadrature<dim> quadrature;
+};
+
+
+
+template <int dim>
+Test<dim>::Test(const Mapping<dim> &mapping)
+ : element(1)
+ , dof_handler(triangulation)
+ , mapping(&mapping)
+{}
+
+
+
+template <int dim>
+void
+Test<dim>::run()
+{
+ setup_single_cartesian_cell_triangulation();
+ dof_handler.distribute_dofs(element);
+ setup_single_point_quadrature();
+
+ test_n_quadrature_points();
+ test_get_quadrature();
+ test_point_mapped_correctly();
+ test_normal();
+ test_JxW();
+ test_shape_surface_grad();
+}
+
+
+
+template <int dim>
+void
+Test<dim>::setup_single_cartesian_cell_triangulation()
+{
+ const Point<dim> lower_left;
+ Point<dim> upper_right;
+ for (unsigned int d = 0; d < dim; ++d)
+ {
+ upper_right[d] = 2 + d;
+ }
+
+ GridGenerator::hyper_rectangle(triangulation, lower_left, upper_right);
+}
+
+
+
+template <int dim>
+void
+Test<dim>::setup_single_point_quadrature()
+{
+ Point<dim> point;
+ const double weight = .5;
+ Tensor<1, dim> normal;
+ for (unsigned int d = 0; d < dim; ++d)
+ {
+ point[d] = .5;
+ normal[d] = d + 1;
+ }
+ normal /= normal.norm();
+ quadrature.push_back(point, weight, normal);
+}
+
+
+
+template <int dim>
+void
+Test<dim>::test_n_quadrature_points()
+{
+ FEImmersedSurfaceValues<dim> fe_values(*mapping,
+ element,
+ quadrature,
+ update_default);
+ fe_values.reinit(triangulation.begin_active());
+
+ deallog << "n_quadrature_points = " << fe_values.n_quadrature_points
+ << std::endl;
+}
+
+
+
+template <int dim>
+void
+Test<dim>::test_get_quadrature()
+{
+ FEImmersedSurfaceValues<dim> fe_values(*mapping,
+ element,
+ quadrature,
+ update_default);
+ fe_values.reinit(triangulation.begin_active());
+
+ const NonMatching::ImmersedSurfaceQuadrature<dim> &stored_quadrature =
+ fe_values.get_quadrature();
+
+ for (unsigned int q = 0; q < stored_quadrature.size(); q++)
+ {
+ deallog << "(point, weight, normal) = ([" << stored_quadrature.point(q)
+ << "], " << stored_quadrature.weight(q) << ", ["
+ << stored_quadrature.normal_vector(q) << "])" << std::endl;
+ }
+}
+
+
+
+template <int dim>
+void
+Test<dim>::test_point_mapped_correctly()
+{
+ FEImmersedSurfaceValues<dim> fe_values(*mapping,
+ element,
+ quadrature,
+ update_quadrature_points);
+ fe_values.reinit(triangulation.begin_active());
+
+ deallog << "point = " << fe_values.quadrature_point(0) << std::endl;
+}
+
+
+
+template <int dim>
+void
+Test<dim>::test_normal()
+{
+ FEImmersedSurfaceValues<dim> fe_values(*mapping,
+ element,
+ quadrature,
+ update_normal_vectors);
+ fe_values.reinit(triangulation.begin_active());
+
+ deallog << "normal = " << fe_values.normal_vector(0) << std::endl;
+}
+
+
+
+template <int dim>
+void
+Test<dim>::test_JxW()
+{
+ FEImmersedSurfaceValues<dim> fe_values(*mapping,
+ element,
+ quadrature,
+ update_JxW_values);
+ fe_values.reinit(triangulation.begin_active());
+
+ deallog << "JxW = " << fe_values.JxW(0) << std::endl;
+}
+
+
+
+template <int dim>
+void
+Test<dim>::test_shape_surface_grad()
+{
+ FEImmersedSurfaceValues<dim> fe_values(
+ *mapping, element, quadrature, update_gradients | update_normal_vectors);
+ fe_values.reinit(dof_handler.begin_active());
+
+ const unsigned int function_index = 0;
+ const unsigned int q_index = 0;
+ const unsigned int component = 0;
+
+ deallog << "shape_surface_grad = "
+ << fe_values.shape_surface_grad(function_index, q_index) << std::endl;
+
+ deallog << "shape_surface_grad_component = "
+ << fe_values.shape_surface_grad_component(function_index,
+ q_index,
+ component)
+ << std::endl;
+}
+
+
+
+template <int dim>
+void
+run_test()
+{
+ deallog << "dim = " << dim << std::endl;
+
+ const MappingCartesian<dim> mapping_cartesian;
+ const unsigned int polynomial_degree = 1;
+ const MappingQ<dim> mapping_q(polynomial_degree);
+
+ const std::vector<std::string> mapping_names = {"MappingCartesian",
+ "MappingQ"};
+ const hp::MappingCollection<dim> mappings(mapping_cartesian, mapping_q);
+
+ for (unsigned int i = 0; i < mappings.size(); i++)
+ {
+ deallog << mapping_names.at(i) << std::endl;
+ Test<dim> test(mappings[i]);
+ test.run();
+ deallog << std::endl;
+ }
+}
+
+
+
+int
+main()
+{
+ initlog();
+ run_test<1>();
+ run_test<2>();
+ run_test<3>();
+}
--- /dev/null
+
+DEAL::dim = 1
+DEAL::MappingCartesian
+DEAL::n_quadrature_points = 1
+DEAL::(point, weight, normal) = ([0.500000], 0.500000, [1.00000])
+DEAL::point = 1.00000
+DEAL::normal = 1.00000
+DEAL::JxW = 0.500000
+DEAL::shape_surface_grad = 0.00000
+DEAL::shape_surface_grad_component = 0.00000
+DEAL::
+DEAL::MappingQ
+DEAL::n_quadrature_points = 1
+DEAL::(point, weight, normal) = ([0.500000], 0.500000, [1.00000])
+DEAL::point = 1.00000
+DEAL::normal = 1.00000
+DEAL::JxW = 0.500000
+DEAL::shape_surface_grad = 0.00000
+DEAL::shape_surface_grad_component = 0.00000
+DEAL::
+DEAL::dim = 2
+DEAL::MappingCartesian
+DEAL::n_quadrature_points = 1
+DEAL::(point, weight, normal) = ([0.500000 0.500000], 0.500000, [0.447214 0.894427])
+DEAL::point = 1.00000 1.50000
+DEAL::normal = 0.600000 0.800000
+DEAL::JxW = 1.11803
+DEAL::shape_surface_grad = -0.0800000 0.0600000
+DEAL::shape_surface_grad_component = -0.0800000 0.0600000
+DEAL::
+DEAL::MappingQ
+DEAL::n_quadrature_points = 1
+DEAL::(point, weight, normal) = ([0.500000 0.500000], 0.500000, [0.447214 0.894427])
+DEAL::point = 1.00000 1.50000
+DEAL::normal = 0.600000 0.800000
+DEAL::JxW = 1.11803
+DEAL::shape_surface_grad = -0.0800000 0.0600000
+DEAL::shape_surface_grad_component = -0.0800000 0.0600000
+DEAL::
+DEAL::dim = 3
+DEAL::MappingCartesian
+DEAL::n_quadrature_points = 1
+DEAL::(point, weight, normal) = ([0.500000 0.500000 0.500000], 0.500000, [0.267261 0.534522 0.801784])
+DEAL::point = 1.00000 1.50000 2.00000
+DEAL::normal = 0.445976 0.594635 0.668965
+DEAL::JxW = 3.59563
+DEAL::shape_surface_grad = -0.0593923 0.00414365 0.0359116
+DEAL::shape_surface_grad_component = -0.0593923 0.00414365 0.0359116
+DEAL::
+DEAL::MappingQ
+DEAL::n_quadrature_points = 1
+DEAL::(point, weight, normal) = ([0.500000 0.500000 0.500000], 0.500000, [0.267261 0.534522 0.801784])
+DEAL::point = 1.00000 1.50000 2.00000
+DEAL::normal = 0.445976 0.594635 0.668965
+DEAL::JxW = 3.59563
+DEAL::shape_surface_grad = -0.0593923 0.00414365 0.0359116
+DEAL::shape_surface_grad_component = -0.0593923 0.00414365 0.0359116
+DEAL::
projects / dealii.git / commitdiff