]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add class NonMatching::FEImmersedSurfaceValues 12024/head
authorSimon Sticko <simon@sticko.se>
Wed, 17 Nov 2021 09:42:24 +0000 (10:42 +0100)
committerSimon Sticko <simon@sticko.se>
Wed, 17 Nov 2021 19:28:19 +0000 (20:28 +0100)
That has the same shape function values/derivatives as FEValues,
but has JxW-values corresponding to integration over a surface
immersed in the cell. This is achieved by passing an
ImmersedSurfaceQuadrature instead of a Quadrature to the constructor.

Add a virtual function fill_fe_immersed_surface_values on Mapping to
be able to compute the correct JxW-values and surface normals. Make it
throw ExcNotImplemented by default, but override it with an
implementation on MappingCartesian and MappingQ.

13 files changed:
include/deal.II/fe/fe.h
include/deal.II/fe/mapping.h
include/deal.II/fe/mapping_cartesian.h
include/deal.II/fe/mapping_q.h
include/deal.II/non_matching/fe_immersed_values.h [new file with mode: 0644]
source/fe/mapping.cc
source/fe/mapping_cartesian.cc
source/fe/mapping_q.cc
source/non_matching/CMakeLists.txt
source/non_matching/fe_immersed_values.cc [new file with mode: 0644]
source/non_matching/fe_immersed_values.inst.in [new file with mode: 0644]
tests/non_matching/fe_immersed_surface_values.cc [new file with mode: 0644]
tests/non_matching/fe_immersed_surface_values.output [new file with mode: 0644]

index 4536a687478ca6800ab2d68ba2b4125181d07abd..453e5e519c763263b026823917acad80a779a3ed 100644 (file)
@@ -40,6 +40,11 @@ template <int dim, int spacedim>
 class FEFaceValues;
 template <int dim, int spacedim>
 class FESubfaceValues;
+namespace NonMatching
+{
+  template <int dim>
+  class FEImmersedSurfaceValues;
+}
 template <int dim, int spacedim>
 class FESystem;
 
@@ -3068,6 +3073,7 @@ protected:
   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
index 6c3809987387690265e8dc6c764432f228bd33c7..f5510b58a23c0e14b9e2ccd09bcb9662bec80823 100644 (file)
@@ -28,6 +28,8 @@
 
 #include <deal.II/hp/q_collection.h>
 
+#include <deal.II/non_matching/immersed_surface_quadrature.h>
+
 #include <array>
 #include <cmath>
 #include <memory>
@@ -48,6 +50,11 @@ template <int dim, int spacedim>
 class FEFaceValues;
 template <int dim, int spacedim>
 class FESubfaceValues;
+namespace NonMatching
+{
+  template <int dim>
+  class FEImmersedSurfaceValues;
+}
 
 
 /**
@@ -981,6 +988,20 @@ protected:
     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;
+
   /**
    * @}
    */
@@ -1282,6 +1303,7 @@ public:
   friend class FEValues<dim, spacedim>;
   friend class FEFaceValues<dim, spacedim>;
   friend class FESubfaceValues<dim, spacedim>;
+  friend class NonMatching::FEImmersedSurfaceValues<dim>;
 };
 
 
index 127d570f2525882554641e9373fd55ebcc46cea6..e2838213d7911f1a25b767e2d96500ef4f1049c7 100644 (file)
@@ -287,6 +287,15 @@ private:
     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;
+
   /**
    * @}
    */
index 29eccf2bad004fbb97acacfbc5b2f79fbd3e6651..7cf23032516b2fc58adb8ddefb4b59a21f583255 100644 (file)
@@ -615,6 +615,15 @@ protected:
     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;
+
   /**
    * @}
    */
diff --git a/include/deal.II/non_matching/fe_immersed_values.h b/include/deal.II/non_matching/fe_immersed_values.h
new file mode 100644 (file)
index 0000000..daf0034
--- /dev/null
@@ -0,0 +1,153 @@
+// ---------------------------------------------------------------------
+//
+// 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 */
index b075dd6b657f56394f7dad266652bf422420fdd7..0aeb35b77e189de02e6130d16efed44c2e8eec2e 100644 (file)
@@ -82,6 +82,20 @@ Mapping<dim, spacedim>::get_bounding_box(
 
 
 
+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(
index 2d29cbc65bd7e635b62381687a41df6038239627..84cc8a3ddc5e4d121b6cd41d91ccc83365cfaad0 100644 (file)
@@ -637,6 +637,69 @@ MappingCartesian<dim, spacedim>::fill_fe_subface_values(
 
 
 
+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(
index 2f4c30aeb45def22039e1cb1229a820479be7f16..ed4bfe9d0e030ad7fad71737b4197bff7d03ff7c 100644 (file)
@@ -1259,6 +1259,139 @@ MappingQ<dim, spacedim>::fill_fe_subface_values(
 
 
 
+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(
index aec821f02a52a20ad423cf9b2028727e646a26bf..aa3d565b8bfc8fe5c046c17c6f9bfb88de9eecb1 100644 (file)
 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
diff --git a/source/non_matching/fe_immersed_values.cc b/source/non_matching/fe_immersed_values.cc
new file mode 100644 (file)
index 0000000..a844126
--- /dev/null
@@ -0,0 +1,216 @@
+// ---------------------------------------------------------------------
+//
+// 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
diff --git a/source/non_matching/fe_immersed_values.inst.in b/source/non_matching/fe_immersed_values.inst.in
new file mode 100644 (file)
index 0000000..639964c
--- /dev/null
@@ -0,0 +1,26 @@
+// ---------------------------------------------------------------------
+//
+// 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>> &);
+  }
diff --git a/tests/non_matching/fe_immersed_surface_values.cc b/tests/non_matching/fe_immersed_surface_values.cc
new file mode 100644 (file)
index 0000000..9e67f21
--- /dev/null
@@ -0,0 +1,305 @@
+// ---------------------------------------------------------------------
+//
+// 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>();
+}
diff --git a/tests/non_matching/fe_immersed_surface_values.output b/tests/non_matching/fe_immersed_surface_values.output
new file mode 100644 (file)
index 0000000..797c582
--- /dev/null
@@ -0,0 +1,58 @@
+
+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::

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.