Implement functions to compute and store mapping information for a vector of cells...

diff --git a/include/deal.II/fe/mapping.h b/include/deal.II/fe/mapping.h
index 9c2568abc6704378b0c1c213c73d911d4643d47f..a4b9578b15e9dc9787ab15d550f6c1d7611bdeb2 100644 (file)
--- a/include/deal.II/fe/mapping.h
+++ b/include/deal.II/fe/mapping.h
@@ -55,7 +55,9 @@ namespace NonMatching
 {
   template <int dim>
   class FEImmersedSurfaceValues;
-}
+  template <int dim, int spacedim>
+  class MappingInfo;
+} // namespace NonMatching
 
 
 /**
@@ -1305,6 +1307,7 @@ public:
   friend class FEFaceValues<dim, spacedim>;
   friend class FESubfaceValues<dim, spacedim>;
   friend class NonMatching::FEImmersedSurfaceValues<dim>;
+  friend class NonMatching::MappingInfo<dim, spacedim>;
 };
 
 

diff --git a/include/deal.II/matrix_free/fe_point_evaluation.h b/include/deal.II/matrix_free/fe_point_evaluation.h
index 2d3aa934643abd944404776f17037c7a04e356e0..33232dba445516f78334c385ac2a5768f777f963 100644 (file)
--- a/include/deal.II/matrix_free/fe_point_evaluation.h
+++ b/include/deal.II/matrix_free/fe_point_evaluation.h
@@ -470,6 +470,20 @@ public:
   reinit(const typename Triangulation<dim, spacedim>::cell_iterator &cell,
          const ArrayView<const Point<dim>> &unit_points);
 
+  /**
+   * Reinitialize the evaluator to point to the correct precomputed mapping of
+   * the cell in the MappingInfo object.
+   */
+  void
+  reinit(const unsigned int cell_index);
+
+  /**
+   * Reinitialize the evaluator to point to the correct precomputed mapping of
+   * the face in the MappingInfo object.
+   */
+  void
+  reinit(const unsigned int cell_index, const unsigned int face_number);
+
   /**
    * This function interpolates the finite element solution, represented by
    * `solution_values`, on the cell and `unit_points` passed to reinit().
@@ -585,6 +599,22 @@ public:
   DerivativeForm<1, spacedim, dim>
   inverse_jacobian(const unsigned int point_index) const;
 
+  /**
+   * Return the Jacobian determinant multiplied by the quadrature weight. This
+   * class or the MappingInfo object passed to this function needs to be
+   * constructed with UpdateFlags containing `update_JxW_values`.
+   */
+  Number
+  JxW(const unsigned int point_index) const;
+
+  /**
+   * Return the normal vector. This class or the MappingInfo object passed to
+   * this function needs to be constructed with UpdateFlags containing
+   * `update_normal_vectors`.
+   */
+  Tensor<1, spacedim>
+  normal_vector(const unsigned int point_index) const;
+
   /**
    * Return the position in real coordinates of the given point index among
    * the points passed to reinit().
@@ -599,6 +629,11 @@ public:
   Point<dim>
   unit_point(const unsigned int point_index) const;
 
+  /**
+   * Number of quadrature points of the current cell/face.
+   */
+  const unsigned int n_q_points;
+
 private:
   /**
    * Common setup function for both constructors. Does the setup for both fast
@@ -713,10 +748,20 @@ private:
    */
   SmartPointer<NonMatching::MappingInfo<dim, spacedim>> mapping_info;
 
+  /**
+   * The current cell index to access mapping data from mapping info.
+   */
+  unsigned int current_cell_index;
+
+  /**
+   * The current face number to access mapping data from mapping info.
+   */
+  unsigned int current_face_number;
+
   /**
    * The reference points specified at reinit().
    */
-  std::vector<Point<dim>> unit_points;
+  ArrayView<const Point<dim>> unit_points;
 
   /**
    * Bool indicating if fast path is chosen.
@@ -733,13 +778,16 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::FEPointEvaluation(
   const FiniteElement<dim> &fe,
   const UpdateFlags         update_flags,
   const unsigned int        first_selected_component)
-  : mapping(&mapping)
+  : n_q_points(numbers::invalid_unsigned_int)
+  , mapping(&mapping)
   , fe(&fe)
   , update_flags(update_flags)
   , mapping_info_on_the_fly(
       std::make_unique<NonMatching::MappingInfo<dim, spacedim>>(mapping,
                                                                 update_flags))
   , mapping_info(mapping_info_on_the_fly.get())
+  , current_cell_index(numbers::invalid_unsigned_int)
+  , current_face_number(numbers::invalid_unsigned_int)
 {
   setup(first_selected_component);
 }
@@ -751,10 +799,13 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::FEPointEvaluation(
   NonMatching::MappingInfo<dim, spacedim> &mapping_info,
   const FiniteElement<dim> &               fe,
   const unsigned int                       first_selected_component)
-  : mapping(&mapping_info.get_mapping())
+  : n_q_points(numbers::invalid_unsigned_int)
+  , mapping(&mapping_info.get_mapping())
   , fe(&fe)
   , update_flags(mapping_info.get_update_flags())
   , mapping_info(&mapping_info)
+  , current_cell_index(numbers::invalid_unsigned_int)
+  , current_face_number(numbers::invalid_unsigned_int)
 {
   setup(first_selected_component);
 }
@@ -853,14 +904,45 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::reinit(
       fe_values->reinit(cell);
     }
 
-  this->unit_points =
-    std::vector<Point<dim>>(unit_points.begin(), unit_points.end());
+  this->unit_points = unit_points;
+
+  const_cast<unsigned int &>(n_q_points) = unit_points.size();
 
   if (update_flags & update_values)
-    values.resize(unit_points.size(), numbers::signaling_nan<value_type>());
+    values.resize(n_q_points, numbers::signaling_nan<value_type>());
   if (update_flags & update_gradients)
-    gradients.resize(unit_points.size(),
-                     numbers::signaling_nan<gradient_type>());
+    gradients.resize(n_q_points, numbers::signaling_nan<gradient_type>());
+}
+
+
+
+template <int n_components, int dim, int spacedim, typename Number>
+void
+FEPointEvaluation<n_components, dim, spacedim, Number>::reinit(
+  const unsigned int cell_index)
+{
+  current_cell_index  = cell_index;
+  current_face_number = numbers::invalid_unsigned_int;
+
+  const_cast<unsigned int &>(n_q_points) =
+    mapping_info->get_unit_points(current_cell_index, current_face_number)
+      .size();
+}
+
+
+
+template <int n_components, int dim, int spacedim, typename Number>
+void
+FEPointEvaluation<n_components, dim, spacedim, Number>::reinit(
+  const unsigned int cell_index,
+  const unsigned int face_number)
+{
+  current_cell_index  = cell_index;
+  current_face_number = face_number;
+
+  const_cast<unsigned int &>(n_q_points) =
+    mapping_info->get_unit_points(current_cell_index, current_face_number)
+      .size();
 }
 
 
@@ -874,7 +956,8 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::evaluate(
   const bool precomputed_mapping = mapping_info_on_the_fly.get() == nullptr;
   if (precomputed_mapping)
     {
-      unit_points = mapping_info->get_unit_points();
+      unit_points =
+        mapping_info->get_unit_points(current_cell_index, current_face_number);
 
       if (update_flags & update_values)
         values.resize(unit_points.size(), numbers::signaling_nan<value_type>());
@@ -883,7 +966,7 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::evaluate(
                          numbers::signaling_nan<gradient_type>());
     }
 
-  if (unit_points.empty())
+  if (unit_points.size() == 0)
     return;
 
   AssertDimension(solution_values.size(), fe->dofs_per_cell);
@@ -946,11 +1029,12 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::evaluate(
                     Number>::set_gradient(val_and_grad.second,
                                           j,
                                           unit_gradients[i + j]);
-                  gradients[i + j] =
-                    apply_transformation(mapping_info->get_mapping_data()
-                                           .inverse_jacobians[i + j]
-                                           .transpose(),
-                                         unit_gradients[i + j]);
+                  const auto &mapping_data =
+                    mapping_info->get_mapping_data(current_cell_index,
+                                                   current_face_number);
+                  gradients[i + j] = apply_transformation(
+                    mapping_data.inverse_jacobians[i + j].transpose(),
+                    unit_gradients[i + j]);
                 }
             }
         }
@@ -1022,7 +1106,8 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::integrate(
   const bool precomputed_mapping = mapping_info_on_the_fly.get() == nullptr;
   if (precomputed_mapping)
     {
-      unit_points = mapping_info->get_unit_points();
+      unit_points =
+        mapping_info->get_unit_points(current_cell_index, current_face_number);
 
       if (update_flags & update_values)
         values.resize(unit_points.size(), numbers::signaling_nan<value_type>());
@@ -1086,9 +1171,12 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::integrate(
           if (integration_flags & EvaluationFlags::gradients)
             for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
               {
-                gradients[i + j] = apply_transformation(
-                  mapping_info->get_mapping_data().inverse_jacobians[i + j],
-                  gradients[i + j]);
+                const auto &mapping_data =
+                  mapping_info->get_mapping_data(current_cell_index,
+                                                 current_face_number);
+                gradients[i + j] =
+                  apply_transformation(mapping_data.inverse_jacobians[i + j],
+                                       gradients[i + j]);
                 internal::FEPointEvaluation::
                   EvaluatorTypeTraits<dim, n_components, Number>::get_gradient(
                     gradient, j, gradients[i + j]);
@@ -1250,9 +1338,10 @@ inline DerivativeForm<1, dim, spacedim>
 FEPointEvaluation<n_components, dim, spacedim, Number>::jacobian(
   const unsigned int point_index) const
 {
-  AssertIndexRange(point_index,
-                   mapping_info->get_mapping_data().jacobians.size());
-  return mapping_info->get_mapping_data().jacobians[point_index];
+  const auto &mapping_data =
+    mapping_info->get_mapping_data(current_cell_index, current_face_number);
+  AssertIndexRange(point_index, mapping_data.jacobians.size());
+  return mapping_data.jacobians[point_index];
 }
 
 
@@ -1262,9 +1351,35 @@ inline DerivativeForm<1, spacedim, dim>
 FEPointEvaluation<n_components, dim, spacedim, Number>::inverse_jacobian(
   const unsigned int point_index) const
 {
-  AssertIndexRange(point_index,
-                   mapping_info->get_mapping_data().inverse_jacobians.size());
-  return mapping_info->get_mapping_data().inverse_jacobians[point_index];
+  const auto &mapping_data =
+    mapping_info->get_mapping_data(current_cell_index, current_face_number);
+  AssertIndexRange(point_index, mapping_data.inverse_jacobians.size());
+  return mapping_data.inverse_jacobians[point_index];
+}
+
+
+
+template <int n_components, int dim, int spacedim, typename Number>
+inline Number
+FEPointEvaluation<n_components, dim, spacedim, Number>::JxW(
+  const unsigned int point_index) const
+{
+  const auto &mapping_data =
+    mapping_info->get_mapping_data(current_cell_index, current_face_number);
+  AssertIndexRange(point_index, mapping_data.JxW_values.size());
+  return mapping_data.JxW_values[point_index];
+}
+
+
+template <int n_components, int dim, int spacedim, typename Number>
+inline Tensor<1, spacedim>
+FEPointEvaluation<n_components, dim, spacedim, Number>::normal_vector(
+  const unsigned int point_index) const
+{
+  const auto &mapping_data =
+    mapping_info->get_mapping_data(current_cell_index, current_face_number);
+  AssertIndexRange(point_index, mapping_data.normal_vectors.size());
+  return mapping_data.normal_vectors[point_index];
 }
 
 
@@ -1274,9 +1389,10 @@ inline Point<spacedim>
 FEPointEvaluation<n_components, dim, spacedim, Number>::real_point(
   const unsigned int point_index) const
 {
-  AssertIndexRange(point_index,
-                   mapping_info->get_mapping_data().quadrature_points.size());
-  return mapping_info->get_mapping_data().quadrature_points[point_index];
+  const auto &mapping_data =
+    mapping_info->get_mapping_data(current_cell_index, current_face_number);
+  AssertIndexRange(point_index, mapping_data.quadrature_points.size());
+  return mapping_data.quadrature_points[point_index];
 }
 
 

diff --git a/include/deal.II/non_matching/mapping_info.h b/include/deal.II/non_matching/mapping_info.h
index 241622d3264b741433b7bbde3de91616483ff620..21347f05303a723e89f7bfb8e181139f78c44860 100644 (file)
--- a/include/deal.II/non_matching/mapping_info.h
+++ b/include/deal.II/non_matching/mapping_info.h
@@ -26,6 +26,7 @@
 
 #include <deal.II/fe/fe_dgq.h>
 #include <deal.II/fe/fe_update_flags.h>
+#include <deal.II/fe/mapping.h>
 #include <deal.II/fe/mapping_cartesian.h>
 #include <deal.II/fe/mapping_q.h>
 #include <deal.II/fe/mapping_related_data.h>
@@ -45,6 +46,10 @@ namespace NonMatching
   class MappingInfo : public Subscriptor
   {
   public:
+    using MappingData =
+      dealii::internal::FEValuesImplementation::MappingRelatedData<dim,
+                                                                   spacedim>;
+
     /**
      * Constructor.
      *
@@ -57,27 +62,113 @@ namespace NonMatching
     MappingInfo(const Mapping<dim> &mapping, const UpdateFlags update_flags);
 
     /**
-     * Reinitialize the mapping information for the incoming cell and unit
+     * Compute the mapping information for the incoming cell and unit
+     * points. This overload is needed to resolve ambiguity.
+     */
+    void
+    reinit(const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+           const std::vector<Point<dim>> &unit_points);
+
+    /**
+     * Compute the mapping information for the incoming cell and unit
      * points.
      */
     void
     reinit(const typename Triangulation<dim, spacedim>::cell_iterator &cell,
            const ArrayView<const Point<dim>> &unit_points);
 
+    /**
+     * Compute the mapping information for the given cell and
+     * quadrature formula. As opposed to the other `reinit` function, this
+     * method allows to access a `JxW` factor at the points.
+     */
+    void
+    reinit(const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+           const Quadrature<dim> &quadrature);
+
+    /**
+     * Compute the mapping information for the incoming vector of cells and
+     * corresponding vector of unit points.
+     *
+     * It is possible to give an IteratorRange<FilteredIterator> to this
+     * function and together with @p n_unfiltered_cells specified this object
+     * can compress its storage while giving you access to the underlying data
+     * with the "unfiltered" index. The default argument for
+     * @p n_unfiltered_cells disables this built-in compression.
+     */
+    template <typename Iterator>
+    void
+    reinit_cells(
+      const IteratorRange<Iterator> &             cell_iterator_range,
+      const std::vector<std::vector<Point<dim>>> &unit_points_vector,
+      const unsigned int n_unfiltered_cells = numbers::invalid_unsigned_int);
+
+    /**
+     * Compute the mapping information for the incoming vector of cells and
+     * corresponding vector of quadratures. As opposed to the other `reinit`
+     * function, this method allows to access a `JxW` factor at the points.
+     */
+    template <typename Iterator>
+    void
+    reinit_cells(
+      const IteratorRange<Iterator> &     cell_iterator_range,
+      const std::vector<Quadrature<dim>> &quadrature_vector,
+      const unsigned int n_unfiltered_cells = numbers::invalid_unsigned_int);
+
+    /**
+     * Compute the mapping information for the incoming vector of cells and
+     * corresponding vector of ImmersedSurfaceQuadrature.
+     */
+    template <typename Iterator>
+    void
+    reinit_surface(
+      const IteratorRange<Iterator> &                    cell_iterator_range,
+      const std::vector<ImmersedSurfaceQuadrature<dim>> &quadrature_vector,
+      const unsigned int n_unfiltered_cells = numbers::invalid_unsigned_int);
+
+    /**
+     * Compute the mapping information for all faces of the incoming vector
+     * of cells and corresponding vector of quadratures.
+     */
+    template <typename Iterator>
+    void
+    reinit_faces(
+      const IteratorRange<Iterator> &                      cell_iterator_range,
+      const std::vector<std::vector<Quadrature<dim - 1>>> &quadrature_vector,
+      const unsigned int n_unfiltered_cells = numbers::invalid_unsigned_int);
+
     /**
      * Getter function for current unit points.
+     *
+     * @p cell_index and @p face_number are the indices
+     * into the compressed, CRS like data storage of unit points.
+     *
+     * If you have initialized this object with reinit_cells() you can access
+     * the stored unit points of the cell with the respective @p cell_index
+     * (and the default argument for @p face_number).
+     *
+     * If you have initialized this object with reinit_faces() you can access
+     * the stored unit points of the faces on the cell with the respective @p cell_index
+     * and the respective local @p face_number.
+     *
+     * If you have initialized this object with reinit() you can access the
+     * stored unit points of a single cell with the default arguments.
+     *
+     * The correct state of this object is checked in this call (in debug mode).
      */
-    const std::vector<Point<dim>> &
-    get_unit_points() const;
+    const ArrayView<const Point<dim>>
+    get_unit_points(
+      const unsigned int cell_index  = numbers::invalid_unsigned_int,
+      const unsigned int face_number = numbers::invalid_unsigned_int) const;
 
     /**
      * Getter function for computed mapping data. This function accesses
      * internal data and is therefore not a stable interface.
      */
-    const dealii::internal::FEValuesImplementation::MappingRelatedData<dim,
-                                                                       spacedim>
-      &
-      get_mapping_data() const;
+    const MappingData &
+    get_mapping_data(
+      const unsigned int cell_index  = numbers::invalid_unsigned_int,
+      const unsigned int face_number = numbers::invalid_unsigned_int) const;
 
     /**
      * Getter function for underlying mapping.
@@ -92,6 +183,23 @@ namespace NonMatching
     get_update_flags() const;
 
   private:
+    /**
+     * Enum class for reinitialized states.
+     */
+    enum class State
+    {
+      invalid,
+      single_cell,
+      cell_vector,
+      faces_on_cells_in_vector
+    };
+
+    /**
+     * Enum class that stores the currently initialized state
+     * upon the last call of reinit().
+     */
+    State state;
+
     /**
      * Compute the mapping related data for the given @p mapping,
      * @p cell and @p unit_points that is required by the FEPointEvaluation
@@ -100,13 +208,42 @@ namespace NonMatching
     void
     compute_mapping_data_for_generic_points(
       const typename Triangulation<dim, spacedim>::cell_iterator &cell,
-      const ArrayView<const Point<dim>> &                         unit_points);
+      const ArrayView<const Point<dim>> &                         unit_points,
+      MappingData &                                               mapping_data);
+
+    /**
+     * Compute the mapping related data for the given @p mapping,
+     * @p cell and @p quadrature that is required by the FEPointEvaluation
+     * class.
+     */
+    void
+    compute_mapping_data_for_immersed_surface_quadrature(
+      const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+      const ImmersedSurfaceQuadrature<dim> &                      quadrature,
+      MappingData &                                               mapping_data);
+
+    /**
+     * Compute the mapping related data for the given @p mapping, @p cell,
+     * @p face_no and @p quadrature that is required by the FEPointEvaluation
+     * class.
+     */
+    void
+    compute_mapping_data_for_face_quadrature(
+      const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+      const unsigned int                                          face_no,
+      const Quadrature<dim - 1> &                                 quadrature,
+      MappingData &                                               mapping_data);
 
     /**
      * The reference points specified at reinit().
      */
     std::vector<Point<dim>> unit_points;
 
+    /**
+     * Offset to point to the first unit point of a cell/face
+     */
+    std::vector<unsigned int> unit_points_index;
+
     /**
      * A pointer to the underlying mapping.
      */
@@ -126,8 +263,24 @@ namespace NonMatching
      * The internal data container for mapping information. The implementation
      * is subject to future changes.
      */
-    dealii::internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
-      mapping_data;
+    std::vector<MappingData> mapping_data;
+
+    /**
+     * Offset to point to the first element of a cell in internal data
+     * containers.
+     */
+    std::vector<unsigned int> cell_index_offset;
+
+    /**
+     * A vector that converts the cell index to a compressed cell index for e.g.
+     * a filtered IteratorRange.
+     */
+    std::vector<unsigned int> cell_index_to_compressed_cell_index;
+
+    /**
+     * A bool that determines weather cell index compression should be done.
+     */
+    bool do_cell_index_compression;
   };
 
   // ----------------------- template functions ----------------------
@@ -141,8 +294,12 @@ namespace NonMatching
   {
     update_flags_mapping = update_default;
     // translate update flags
-    if (update_flags & update_jacobians)
+    if (update_flags & update_jacobians || update_flags & update_JxW_values)
       update_flags_mapping |= update_jacobians;
+    if (update_flags & update_JxW_values)
+      update_flags_mapping |= update_JxW_values;
+    if (update_flags & update_normal_vectors)
+      update_flags_mapping |= update_normal_vectors;
     if (update_flags & update_gradients ||
         update_flags & update_inverse_jacobians)
       update_flags_mapping |= update_inverse_jacobians;
@@ -157,30 +314,451 @@ namespace NonMatching
   void
   MappingInfo<dim, spacedim>::reinit(
     const typename Triangulation<dim, spacedim>::cell_iterator &cell,
-    const ArrayView<const Point<dim>> &                         unit_points)
+    const std::vector<Point<dim>> &                             unit_points_in)
+  {
+    reinit(cell, make_array_view(unit_points_in.begin(), unit_points_in.end()));
+  }
+
+
+
+  template <int dim, int spacedim>
+  void
+  MappingInfo<dim, spacedim>::reinit(
+    const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+    const ArrayView<const Point<dim>> &                         unit_points_in)
+  {
+    unit_points =
+      std::vector<Point<dim>>(unit_points_in.begin(), unit_points_in.end());
+
+    mapping_data.resize(1);
+    compute_mapping_data_for_generic_points(cell,
+                                            unit_points_in,
+                                            mapping_data[0]);
+
+    state = State::single_cell;
+  }
+
+
+
+  template <int dim, int spacedim>
+  void
+  MappingInfo<dim, spacedim>::reinit(
+    const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+    const Quadrature<dim> &                                     quadrature)
+  {
+    const auto &points  = quadrature.get_points();
+    const auto &weights = quadrature.get_weights();
+
+    reinit(cell, points);
+
+    if (update_flags_mapping & update_JxW_values)
+      for (unsigned int q = 0; q < points.size(); ++q)
+        mapping_data[0].JxW_values[q] =
+          determinant(Tensor<2, dim>(mapping_data[0].jacobians[q])) *
+          weights[q];
+  }
+
+
+
+  template <int dim, int spacedim>
+  template <typename Iterator>
+  void
+  MappingInfo<dim, spacedim>::reinit_cells(
+    const IteratorRange<Iterator> &             cell_iterator_range,
+    const std::vector<std::vector<Point<dim>>> &unit_points_vector,
+    const unsigned int                          n_unfiltered_cells)
   {
-    this->unit_points =
-      std::vector<Point<dim>>(unit_points.begin(), unit_points.end());
-    compute_mapping_data_for_generic_points(cell, unit_points);
+    do_cell_index_compression =
+      n_unfiltered_cells != numbers::invalid_unsigned_int;
+
+    const unsigned int n_cells = unit_points_vector.size();
+    AssertDimension(n_cells,
+                    std::distance(cell_iterator_range.begin(),
+                                  cell_iterator_range.end()));
+
+    // fill unit points index offset vector
+    unit_points_index.reserve(n_cells + 1);
+    unit_points_index.push_back(0);
+    for (const auto &unit_points : unit_points_vector)
+      unit_points_index.push_back(unit_points_index.back() +
+                                  unit_points.size());
+
+    const unsigned int n_unit_points = unit_points_index.back();
+
+    unit_points.resize(n_unit_points);
+    mapping_data.resize(n_cells);
+
+    if (do_cell_index_compression)
+      cell_index_to_compressed_cell_index.resize(n_unfiltered_cells,
+                                                 numbers::invalid_unsigned_int);
+    unsigned int cell_index = 0;
+    for (const auto &cell : cell_iterator_range)
+      {
+        auto it = unit_points.begin() + unit_points_index[cell_index];
+        for (const auto &unit_point : unit_points_vector[cell_index])
+          {
+            *it = unit_point;
+            ++it;
+          }
+
+        compute_mapping_data_for_generic_points(cell,
+                                                unit_points_vector[cell_index],
+                                                mapping_data[cell_index]);
+
+        if (do_cell_index_compression)
+          cell_index_to_compressed_cell_index[cell->active_cell_index()] =
+            cell_index;
+
+        ++cell_index;
+      }
+
+    state = State::cell_vector;
   }
 
 
 
   template <int dim, int spacedim>
-  const std::vector<Point<dim>> &
-  MappingInfo<dim, spacedim>::get_unit_points() const
+  template <typename Iterator>
+  void
+  MappingInfo<dim, spacedim>::reinit_cells(
+    const IteratorRange<Iterator> &     cell_iterator_range,
+    const std::vector<Quadrature<dim>> &quadrature_vector,
+    const unsigned int                  n_unfiltered_cells)
   {
-    return unit_points;
+    const unsigned int n_cells = quadrature_vector.size();
+    AssertDimension(n_cells,
+                    std::distance(cell_iterator_range.begin(),
+                                  cell_iterator_range.end()));
+
+    std::vector<std::vector<Point<dim>>> unit_points_vector(n_cells);
+    for (unsigned int cell_index = 0; cell_index < n_cells; ++cell_index)
+      unit_points_vector[cell_index] = std::vector<Point<dim>>(
+        quadrature_vector[cell_index].get_points().begin(),
+        quadrature_vector[cell_index].get_points().end());
+
+    reinit_cells(cell_iterator_range, unit_points_vector, n_unfiltered_cells);
+
+    if (update_flags_mapping & update_JxW_values)
+      for (unsigned int cell_index = 0; cell_index < n_cells; ++cell_index)
+        {
+          const auto &weights = quadrature_vector[cell_index].get_weights();
+          for (unsigned int q = 0; q < weights.size(); ++q)
+            mapping_data[cell_index].JxW_values[q] =
+              determinant(
+                Tensor<2, dim>(mapping_data[cell_index].jacobians[q])) *
+              weights[q];
+        }
   }
 
 
 
   template <int dim, int spacedim>
-  const dealii::internal::FEValuesImplementation::MappingRelatedData<dim,
-                                                                     spacedim> &
-  MappingInfo<dim, spacedim>::get_mapping_data() const
+  template <typename Iterator>
+  void
+  MappingInfo<dim, spacedim>::reinit_surface(
+    const IteratorRange<Iterator> &                    cell_iterator_range,
+    const std::vector<ImmersedSurfaceQuadrature<dim>> &quadrature_vector,
+    const unsigned int                                 n_unfiltered_cells)
   {
-    return mapping_data;
+    do_cell_index_compression =
+      n_unfiltered_cells != numbers::invalid_unsigned_int;
+
+    if (update_flags_mapping & (update_JxW_values | update_normal_vectors))
+      update_flags_mapping |= update_covariant_transformation;
+
+    const unsigned int n_cells = quadrature_vector.size();
+    AssertDimension(n_cells,
+                    std::distance(cell_iterator_range.begin(),
+                                  cell_iterator_range.end()));
+
+    // fill unit points index offset vector
+    unit_points_index.reserve(n_cells + 1);
+    unit_points_index.push_back(0);
+    for (const auto &quadrature : quadrature_vector)
+      unit_points_index.push_back(unit_points_index.back() +
+                                  quadrature.get_points().size());
+
+    const unsigned int n_unit_points = unit_points_index.back();
+
+    unit_points.resize(n_unit_points);
+    mapping_data.resize(n_cells);
+
+    if (do_cell_index_compression)
+      cell_index_to_compressed_cell_index.resize(n_unfiltered_cells,
+                                                 numbers::invalid_unsigned_int);
+    unsigned int cell_index = 0;
+    for (const auto &cell : cell_iterator_range)
+      {
+        const auto &quadrature = quadrature_vector[cell_index];
+
+        auto it = unit_points.begin() + unit_points_index[cell_index];
+        for (const auto &unit_point : quadrature.get_points())
+          {
+            *it = unit_point;
+            ++it;
+          }
+
+        compute_mapping_data_for_immersed_surface_quadrature(
+          cell, quadrature, mapping_data[cell_index]);
+
+        if (do_cell_index_compression)
+          cell_index_to_compressed_cell_index[cell->active_cell_index()] =
+            cell_index;
+
+        ++cell_index;
+      }
+
+    state = State::cell_vector;
+  }
+
+
+
+  template <int dim, int spacedim>
+  template <typename Iterator>
+  void
+  MappingInfo<dim, spacedim>::reinit_faces(
+    const IteratorRange<Iterator> &                      cell_iterator_range,
+    const std::vector<std::vector<Quadrature<dim - 1>>> &quadrature_vector,
+    const unsigned int                                   n_unfiltered_cells)
+  {
+    do_cell_index_compression =
+      n_unfiltered_cells != numbers::invalid_unsigned_int;
+
+    const unsigned int n_cells = quadrature_vector.size();
+    AssertDimension(n_cells,
+                    std::distance(cell_iterator_range.begin(),
+                                  cell_iterator_range.end()));
+
+    // fill cell index offset vector
+    cell_index_offset.resize(n_cells);
+    unsigned int n_faces    = 0;
+    unsigned int cell_index = 0;
+    for (const auto &cell : cell_iterator_range)
+      {
+        cell_index_offset[cell_index] = n_faces;
+        n_faces += cell->n_faces();
+        ++cell_index;
+      }
+
+    // fill unit points index offset vector
+    unit_points_index.resize(n_faces + 1);
+    cell_index                 = 0;
+    unsigned int n_unit_points = 0;
+    for (const auto &cell : cell_iterator_range)
+      {
+        for (const auto &f : cell->face_indices())
+          {
+            const unsigned int current_face_index =
+              cell_index_offset[cell_index] + f;
+
+            unit_points_index[current_face_index] = n_unit_points;
+            n_unit_points +=
+              quadrature_vector[cell_index][f].get_points().size();
+          }
+
+        ++cell_index;
+      }
+    unit_points_index[n_faces] = n_unit_points;
+
+    // compress indices
+    if (do_cell_index_compression)
+      cell_index_to_compressed_cell_index.resize(n_unfiltered_cells,
+                                                 numbers::invalid_unsigned_int);
+
+    // fill unit points and mapping data for every face of all cells
+    unit_points.resize(n_unit_points);
+    mapping_data.resize(n_faces);
+    cell_index = 0;
+    QProjector<dim> q_projector;
+    for (const auto &cell : cell_iterator_range)
+      {
+        const auto &quadratures_on_faces = quadrature_vector[cell_index];
+
+        Assert(quadratures_on_faces.size() == cell->n_faces(),
+               ExcDimensionMismatch(quadratures_on_faces.size(),
+                                    cell->n_faces()));
+
+        for (const auto &f : cell->face_indices())
+          {
+            const auto &quadrature_on_face = quadratures_on_faces[f];
+
+            const auto quadrature_on_cell =
+              q_projector.project_to_face(cell->reference_cell(),
+                                          quadrature_on_face,
+                                          f);
+
+            const auto &unit_points_on_cell = quadrature_on_cell.get_points();
+
+            const unsigned int current_face_index =
+              cell_index_offset[cell_index] + f;
+
+            auto it =
+              unit_points.begin() + unit_points_index[current_face_index];
+            for (const auto &unit_point : unit_points_on_cell)
+              {
+                *it = unit_point;
+                ++it;
+              }
+
+            compute_mapping_data_for_face_quadrature(
+              cell, f, quadrature_on_face, mapping_data[current_face_index]);
+          }
+        if (do_cell_index_compression)
+          cell_index_to_compressed_cell_index[cell->active_cell_index()] =
+            cell_index;
+
+        ++cell_index;
+      }
+
+    state = State::faces_on_cells_in_vector;
+  }
+
+
+
+  template <int dim, int spacedim>
+  const ArrayView<const Point<dim>>
+  MappingInfo<dim, spacedim>::get_unit_points(
+    const unsigned int cell_index,
+    const unsigned int face_number) const
+  {
+    if (cell_index == numbers::invalid_unsigned_int &&
+        face_number == numbers::invalid_unsigned_int)
+      {
+        Assert(state == State::single_cell,
+               ExcMessage(
+                 "This mapping info is not reinitialized for a single cell!"));
+        return unit_points;
+      }
+    else if (face_number == numbers::invalid_unsigned_int)
+      {
+        Assert(state == State::cell_vector,
+               ExcMessage(
+                 "This mapping info is not reinitialized for a cell vector!"));
+        if (do_cell_index_compression)
+          {
+            Assert(cell_index_to_compressed_cell_index[cell_index] !=
+                     numbers::invalid_unsigned_int,
+                   ExcMessage("Mapping info object was not initialized for this"
+                              " active cell index!"));
+            const auto it_begin =
+              unit_points.begin() +
+              unit_points_index
+                [cell_index_to_compressed_cell_index[cell_index]];
+            const auto it_end =
+              unit_points.begin() +
+              unit_points_index
+                [cell_index_to_compressed_cell_index[cell_index] + 1];
+            return make_array_view(it_begin, it_end);
+          }
+        else
+          {
+            const auto it_begin =
+              unit_points.begin() + unit_points_index[cell_index];
+            const auto it_end =
+              unit_points.begin() + unit_points_index[cell_index + 1];
+            return make_array_view(it_begin, it_end);
+          }
+      }
+    else if (cell_index != numbers::invalid_unsigned_int)
+      {
+        Assert(state == State::faces_on_cells_in_vector,
+               ExcMessage("This mapping info is not reinitialized for faces"
+                          " on cells in a vector!"));
+        if (do_cell_index_compression)
+          {
+            Assert(
+              cell_index_to_compressed_cell_index[cell_index] !=
+                numbers::invalid_unsigned_int,
+              ExcMessage(
+                "Mapping info object was not initialized for this active cell index"
+                " and corresponding face numbers!"));
+            const unsigned int current_face_index =
+              cell_index_offset
+                [cell_index_to_compressed_cell_index[cell_index]] +
+              face_number;
+            const auto it_begin =
+              unit_points.begin() + unit_points_index[current_face_index];
+            const auto it_end =
+              unit_points.begin() + unit_points_index[current_face_index + 1];
+            return make_array_view(it_begin, it_end);
+          }
+        else
+          {
+            const unsigned int current_face_index =
+              cell_index_offset[cell_index] + face_number;
+            const auto it_begin =
+              unit_points.begin() + unit_points_index[current_face_index];
+            const auto it_end =
+              unit_points.begin() + unit_points_index[current_face_index + 1];
+            return make_array_view(it_begin, it_end);
+          }
+      }
+    else
+      AssertThrow(
+        false,
+        ExcMessage(
+          "cell_index has to be specified if face_number is specified!"));
+  }
+
+
+
+  template <int dim, int spacedim>
+  const typename MappingInfo<dim, spacedim>::MappingData &
+  MappingInfo<dim, spacedim>::get_mapping_data(
+    const unsigned int cell_index,
+    const unsigned int face_number) const
+  {
+    if (cell_index == numbers::invalid_unsigned_int &&
+        face_number == numbers::invalid_unsigned_int)
+      {
+        Assert(state == State::single_cell,
+               ExcMessage(
+                 "This mapping info is not reinitialized for a single cell!"));
+        return mapping_data[0];
+      }
+    else if (face_number == numbers::invalid_unsigned_int)
+      {
+        Assert(state == State::cell_vector,
+               ExcMessage(
+                 "This mapping info is not reinitialized for a cell vector!"));
+        if (do_cell_index_compression)
+          {
+            Assert(cell_index_to_compressed_cell_index[cell_index] !=
+                     numbers::invalid_unsigned_int,
+                   ExcMessage("Mapping info object was not initialized for this"
+                              " active cell index!"));
+            return mapping_data
+              [cell_index_to_compressed_cell_index[cell_index]];
+          }
+        else
+          return mapping_data[cell_index];
+      }
+    else if (cell_index != numbers::invalid_unsigned_int)
+      {
+        Assert(state == State::faces_on_cells_in_vector,
+               ExcMessage("This mapping info is not reinitialized for faces"
+                          " on cells in a vector!"));
+        if (do_cell_index_compression)
+          {
+            Assert(
+              cell_index_to_compressed_cell_index[cell_index] !=
+                numbers::invalid_unsigned_int,
+              ExcMessage(
+                "Mapping info object was not initialized for this active cell index"
+                " and corresponding face numbers!"));
+            return mapping_data
+              [cell_index_offset
+                 [cell_index_to_compressed_cell_index[cell_index]] +
+               face_number];
+          }
+        else
+          return mapping_data[cell_index_offset[cell_index] + face_number];
+      }
+    else
+      AssertThrow(
+        false,
+        ExcMessage(
+          "cell_index has to be specified if face_number is specified!"));
   }
 
 
@@ -207,7 +785,8 @@ namespace NonMatching
   void
   MappingInfo<dim, spacedim>::compute_mapping_data_for_generic_points(
     const typename Triangulation<dim, spacedim>::cell_iterator &cell,
-    const ArrayView<const Point<dim>> &                         unit_points)
+    const ArrayView<const Point<dim>> &                         unit_points,
+    MappingData &                                               mapping_data)
   {
     if (const MappingQ<dim, spacedim> *mapping_q =
           dynamic_cast<const MappingQ<dim, spacedim> *>(&(*mapping)))
@@ -246,6 +825,58 @@ namespace NonMatching
             mapping_data.quadrature_points[q] = fe_values.quadrature_point(q);
       }
   }
+
+
+
+  template <int dim, int spacedim>
+  void
+  MappingInfo<dim, spacedim>::
+    compute_mapping_data_for_immersed_surface_quadrature(
+      const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+      const ImmersedSurfaceQuadrature<dim> &                      quadrature,
+      MappingData &                                               mapping_data)
+  {
+    update_flags_mapping |=
+      mapping->requires_update_flags(update_flags_mapping);
+
+    mapping_data.initialize(quadrature.get_points().size(),
+                            update_flags_mapping);
+
+    auto internal_mapping_data =
+      mapping->get_data(update_flags_mapping, quadrature);
+
+    mapping->fill_fe_immersed_surface_values(cell,
+                                             quadrature,
+                                             *internal_mapping_data,
+                                             mapping_data);
+  }
+
+
+
+  template <int dim, int spacedim>
+  void
+  MappingInfo<dim, spacedim>::compute_mapping_data_for_face_quadrature(
+    const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+    const unsigned int                                          face_no,
+    const Quadrature<dim - 1> &                                 quadrature,
+    MappingData &                                               mapping_data)
+  {
+    update_flags_mapping |=
+      mapping->requires_update_flags(update_flags_mapping);
+
+    mapping_data.initialize(quadrature.get_points().size(),
+                            update_flags_mapping);
+
+    auto internal_mapping_data =
+      mapping->get_face_data(update_flags_mapping,
+                             hp::QCollection<dim - 1>(quadrature));
+
+    mapping->fill_fe_face_values(cell,
+                                 face_no,
+                                 hp::QCollection<dim - 1>(quadrature),
+                                 *internal_mapping_data,
+                                 mapping_data);
+  }
 } // namespace NonMatching
 
 DEAL_II_NAMESPACE_CLOSE

diff --git a/tests/non_matching/mapping_info.cc b/tests/non_matching/mapping_info.cc
new file mode 100644 (file)
index 0000000..f6c7988
--- /dev/null
+++ b/tests/non_matching/mapping_info.cc
@@ -0,0 +1,477 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2023 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.
+//
+// ---------------------------------------------------------------------
+
+/*
+ * Test the NonMatching::MappingInfo class together with FEPointEvaluation and
+ * compare to NonMatching::FEValues
+ */
+
+#include <deal.II/base/function_signed_distance.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/filtered_iterator.h>
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/matrix_free/fe_point_evaluation.h>
+
+#include <deal.II/non_matching/fe_values.h>
+#include <deal.II/non_matching/mapping_info.h>
+#include <deal.II/non_matching/mesh_classifier.h>
+#include <deal.II/non_matching/quadrature_generator.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+void
+test(const bool filtered_compression)
+{
+  constexpr unsigned int dim    = 2;
+  constexpr unsigned int degree = 1;
+
+  FE_Q<dim> fe_q(degree);
+
+  Triangulation<dim> tria;
+
+  MappingQ<dim> mapping(degree);
+
+  GridGenerator::subdivided_hyper_cube(tria, 4);
+
+  DoFHandler<dim> dof_handler(tria);
+
+  dof_handler.distribute_dofs(fe_q);
+
+  Functions::SignedDistance::Sphere<dim> level_set;
+
+  Vector<double> level_set_vec(dof_handler.n_dofs());
+
+  VectorTools::interpolate(dof_handler, level_set, level_set_vec);
+
+  NonMatching::MeshClassifier<dim> mesh_classifier(dof_handler, level_set_vec);
+  mesh_classifier.reclassify();
+
+  hp::QCollection<1> q_collection((QGauss<1>(degree)));
+
+  NonMatching::DiscreteQuadratureGenerator<dim> quadrature_generator(
+    q_collection, dof_handler, level_set_vec);
+
+  NonMatching::DiscreteFaceQuadratureGenerator<dim> face_quadrature_generator(
+    q_collection, dof_handler, level_set_vec);
+
+  // FEPointEvaluation
+  NonMatching::MappingInfo<dim> mapping_info_cell(
+    mapping, update_values | update_gradients | update_JxW_values);
+
+  NonMatching::MappingInfo<dim> mapping_info_surface(mapping,
+                                                     update_values |
+                                                       update_gradients |
+                                                       update_JxW_values |
+                                                       update_normal_vectors);
+
+  NonMatching::MappingInfo<dim> mapping_info_faces(mapping,
+                                                   update_values |
+                                                     update_gradients |
+                                                     update_JxW_values |
+                                                     update_normal_vectors);
+
+  auto physical =
+    [&](const typename Triangulation<dim>::active_cell_iterator &i) {
+      return (mesh_classifier.location_to_level_set(i) ==
+              NonMatching::LocationToLevelSet::intersected) ||
+             (mesh_classifier.location_to_level_set(i) ==
+              NonMatching::LocationToLevelSet::inside);
+    };
+
+  if (filtered_compression)
+    {
+      const auto physical_active_cell_iterators =
+        tria.active_cell_iterators() | physical;
+
+      unsigned int n_physical_cells = 0;
+      for (const auto &cell : physical_active_cell_iterators)
+        {
+          ++n_physical_cells;
+        }
+
+      std::vector<Quadrature<dim>> quad_vec_cell;
+      quad_vec_cell.reserve(n_physical_cells);
+      std::vector<NonMatching::ImmersedSurfaceQuadrature<dim>> quad_vec_surface;
+      quad_vec_surface.reserve(n_physical_cells);
+      std::vector<std::vector<Quadrature<dim - 1>>> quad_vec_faces(
+        n_physical_cells);
+
+      unsigned int counter = 0;
+      for (const auto &cell : physical_active_cell_iterators)
+        {
+          quadrature_generator.generate(cell);
+          quad_vec_cell.push_back(quadrature_generator.get_inside_quadrature());
+          quad_vec_surface.push_back(
+            quadrature_generator.get_surface_quadrature());
+
+          for (auto f : cell->face_indices())
+            {
+              face_quadrature_generator.generate(cell, f);
+              quad_vec_faces[counter].push_back(
+                face_quadrature_generator.get_inside_quadrature());
+            }
+
+          ++counter;
+        }
+
+      mapping_info_cell.reinit_cells(physical_active_cell_iterators,
+                                     quad_vec_cell,
+                                     tria.n_active_cells());
+      mapping_info_surface.reinit_surface(physical_active_cell_iterators,
+                                          quad_vec_surface,
+                                          tria.n_active_cells());
+      mapping_info_faces.reinit_faces(physical_active_cell_iterators,
+                                      quad_vec_faces,
+                                      tria.n_active_cells());
+    }
+  else
+    {
+      std::vector<Quadrature<dim>>                             quad_vec_cell;
+      std::vector<NonMatching::ImmersedSurfaceQuadrature<dim>> quad_vec_surface;
+      std::vector<std::vector<Quadrature<dim - 1>>>            quad_vec_faces(
+        tria.n_active_cells());
+      for (const auto &cell : tria.active_cell_iterators())
+        {
+          quadrature_generator.generate(cell);
+          quad_vec_cell.push_back(quadrature_generator.get_inside_quadrature());
+          quad_vec_surface.push_back(
+            quadrature_generator.get_surface_quadrature());
+
+          for (auto f : cell->face_indices())
+            {
+              face_quadrature_generator.generate(cell, f);
+              quad_vec_faces[cell->active_cell_index()].push_back(
+                face_quadrature_generator.get_inside_quadrature());
+            }
+        }
+
+      mapping_info_cell.reinit_cells(tria.active_cell_iterators(),
+                                     quad_vec_cell);
+      mapping_info_surface.reinit_surface(tria.active_cell_iterators(),
+                                          quad_vec_surface);
+      mapping_info_faces.reinit_faces(tria.active_cell_iterators(),
+                                      quad_vec_faces);
+    }
+
+  Vector<double> src(dof_handler.n_dofs()), dst_cell(dof_handler.n_dofs()),
+    dst_surface(dof_handler.n_dofs()), dst_faces(dof_handler.n_dofs());
+
+  for (auto &v : src)
+    v = random_value<double>();
+
+  FEPointEvaluation<1, dim, dim, double> fe_point_cell(mapping_info_cell, fe_q);
+  FEPointEvaluation<1, dim, dim, double> fe_point_surface(mapping_info_surface,
+                                                          fe_q);
+  FEPointEvaluation<1, dim, dim, double> fe_point_faces_m(mapping_info_faces,
+                                                          fe_q);
+  FEPointEvaluation<1, dim, dim, double> fe_point_faces_p(mapping_info_faces,
+                                                          fe_q);
+
+  std::vector<double> solution_values_in(fe_q.dofs_per_cell);
+  std::vector<double> solution_values_neighbor_in(fe_q.dofs_per_cell);
+  std::vector<double> solution_values_cell_out(fe_q.dofs_per_cell);
+  std::vector<double> solution_values_surface_out(fe_q.dofs_per_cell);
+  std::vector<double> solution_values_faces_out(fe_q.dofs_per_cell);
+  for (const auto &cell : dof_handler.active_cell_iterators())
+    {
+      if (!(physical(cell)))
+        continue;
+
+      cell->get_dof_values(src,
+                           solution_values_in.begin(),
+                           solution_values_in.end());
+
+      auto test_fe_point = [&](FEPointEvaluation<1, dim, dim, double> &fe_point,
+                               std::vector<double> &solution_values_out) {
+        fe_point.reinit(cell->active_cell_index());
+
+        fe_point.evaluate(solution_values_in,
+                          EvaluationFlags::values | EvaluationFlags::gradients);
+
+        for (unsigned int q = 0; q < fe_point.n_q_points; ++q)
+          {
+            fe_point.submit_value(fe_point.JxW(q) * fe_point.get_value(q), q);
+            fe_point.submit_gradient(fe_point.JxW(q) * fe_point.get_gradient(q),
+                                     q);
+          }
+
+        fe_point.integrate(solution_values_out,
+                           EvaluationFlags::values |
+                             EvaluationFlags::gradients);
+      };
+
+      test_fe_point(fe_point_cell, solution_values_cell_out);
+      test_fe_point(fe_point_surface, solution_values_surface_out);
+
+      for (const auto f : cell->face_indices())
+        {
+          if (cell->at_boundary(f) || !physical(cell->neighbor(f)))
+            continue;
+
+          fe_point_faces_m.reinit(cell->active_cell_index(), f);
+          fe_point_faces_p.reinit(cell->neighbor(f)->active_cell_index(),
+                                  cell->neighbor_of_neighbor(f));
+
+          cell->neighbor(f)->get_dof_values(src,
+                                            solution_values_neighbor_in.begin(),
+                                            solution_values_neighbor_in.end());
+
+          fe_point_faces_m.evaluate(solution_values_in,
+                                    EvaluationFlags::values |
+                                      EvaluationFlags::gradients);
+
+          fe_point_faces_p.evaluate(solution_values_neighbor_in,
+                                    EvaluationFlags::values |
+                                      EvaluationFlags::gradients);
+
+          for (unsigned int q = 0; q < fe_point_faces_m.n_q_points; ++q)
+            {
+              fe_point_faces_m.submit_value(fe_point_faces_m.JxW(q) *
+                                              (fe_point_faces_m.get_value(q) -
+                                               fe_point_faces_p.get_value(q)),
+                                            q);
+              fe_point_faces_m.submit_gradient(
+                fe_point_faces_m.JxW(q) * (fe_point_faces_m.get_gradient(q) -
+                                           fe_point_faces_p.get_gradient(q)),
+                q);
+            }
+
+          fe_point_faces_m.integrate(solution_values_faces_out,
+                                     EvaluationFlags::values |
+                                       EvaluationFlags::gradients);
+
+          cell->distribute_local_to_global(
+            Vector<double>(solution_values_faces_out.begin(),
+                           solution_values_faces_out.end()),
+            dst_faces);
+        }
+
+      cell->distribute_local_to_global(
+        Vector<double>(solution_values_cell_out.begin(),
+                       solution_values_cell_out.end()),
+        dst_cell);
+
+      cell->distribute_local_to_global(
+        Vector<double>(solution_values_surface_out.begin(),
+                       solution_values_surface_out.end()),
+        dst_surface);
+    }
+
+
+  // FEValues
+  const QGauss<1> quadrature_1D(degree);
+
+  NonMatching::RegionUpdateFlags region_update_flags;
+  region_update_flags.inside = update_values | update_gradients |
+                               update_JxW_values | update_quadrature_points;
+  region_update_flags.surface = update_values | update_gradients |
+                                update_JxW_values | update_quadrature_points |
+                                update_normal_vectors;
+
+  hp::FECollection<dim> fe_collection(fe_q);
+
+  NonMatching::FEValues<dim> non_matching_fe_values(fe_collection,
+                                                    quadrature_1D,
+                                                    region_update_flags,
+                                                    mesh_classifier,
+                                                    dof_handler,
+                                                    level_set_vec);
+
+  NonMatching::FEInterfaceValues<dim> non_matching_fe_interface_values(
+    fe_collection,
+    quadrature_1D,
+    region_update_flags,
+    mesh_classifier,
+    dof_handler,
+    level_set_vec);
+
+  Vector<double> dst_cell_2(dof_handler.n_dofs()),
+    dst_surface_2(dof_handler.n_dofs()), dst_faces_2(dof_handler.n_dofs());
+
+  for (const auto &cell : dof_handler.active_cell_iterators())
+    {
+      non_matching_fe_values.reinit(cell);
+
+      cell->get_dof_values(src,
+                           solution_values_in.begin(),
+                           solution_values_in.end());
+
+      const auto &inside_fe_values =
+        non_matching_fe_values.get_inside_fe_values();
+
+      const auto &surface_fe_values =
+        non_matching_fe_values.get_surface_fe_values();
+
+
+      auto test_fe_values = [&](const auto &         fe_values,
+                                std::vector<double> &solution_values_out,
+                                Vector<double> &     dst) {
+        if (fe_values)
+          {
+            std::vector<Tensor<1, dim>> solution_gradients(
+              fe_values->n_quadrature_points);
+            std::vector<double> solution_values(fe_values->n_quadrature_points);
+
+            for (const auto q : fe_values->quadrature_point_indices())
+              {
+                double         values = 0.;
+                Tensor<1, dim> gradients;
+
+                for (const auto i : fe_values->dof_indices())
+                  {
+                    gradients +=
+                      solution_values_in[i] * fe_values->shape_grad(i, q);
+                    values +=
+                      solution_values_in[i] * fe_values->shape_value(i, q);
+                  }
+                solution_gradients[q] = gradients * fe_values->JxW(q);
+                solution_values[q]    = values * fe_values->JxW(q);
+              }
+
+            for (const auto i : fe_values->dof_indices())
+              {
+                double sum_gradients = 0.;
+                double sum_values    = 0.;
+                for (const auto q : fe_values->quadrature_point_indices())
+                  {
+                    sum_gradients +=
+                      solution_gradients[q] * fe_values->shape_grad(i, q);
+                    sum_values +=
+                      solution_values[q] * fe_values->shape_value(i, q);
+                  }
+
+                solution_values_cell_out[i] = sum_gradients + sum_values;
+              }
+
+            cell->distribute_local_to_global(
+              Vector<double>(solution_values_cell_out.begin(),
+                             solution_values_cell_out.end()),
+              dst);
+          }
+      };
+
+      test_fe_values(inside_fe_values, solution_values_cell_out, dst_cell_2);
+      test_fe_values(surface_fe_values,
+                     solution_values_surface_out,
+                     dst_surface_2);
+
+      for (const auto f : cell->face_indices())
+        {
+          if (cell->at_boundary(f))
+            continue;
+
+          non_matching_fe_interface_values.reinit(
+            cell,
+            f,
+            numbers::invalid_unsigned_int,
+            cell->neighbor(f),
+            cell->neighbor_of_neighbor(f),
+            numbers::invalid_unsigned_int);
+
+          const auto &fe_interface_values =
+            non_matching_fe_interface_values.get_inside_fe_values();
+
+          if (fe_interface_values)
+            {
+              const auto &fe_face_values_m =
+                fe_interface_values->get_fe_face_values(0);
+
+              const auto &fe_face_values_p =
+                fe_interface_values->get_fe_face_values(1);
+
+              cell->neighbor(f)->get_dof_values(
+                src,
+                solution_values_neighbor_in.begin(),
+                solution_values_neighbor_in.end());
+
+              std::vector<Tensor<1, dim>> solution_gradients(
+                fe_face_values_m.n_quadrature_points);
+              std::vector<double> solution_values(
+                fe_face_values_m.n_quadrature_points);
+
+              for (const auto q : fe_face_values_m.quadrature_point_indices())
+                {
+                  double         values = 0.;
+                  Tensor<1, dim> gradients;
+
+                  for (const auto i : fe_face_values_m.dof_indices())
+                    {
+                      gradients += solution_values_in[i] *
+                                     fe_face_values_m.shape_grad(i, q) -
+                                   solution_values_neighbor_in[i] *
+                                     fe_face_values_p.shape_grad(i, q);
+                      values += solution_values_in[i] *
+                                  fe_face_values_m.shape_value(i, q) -
+                                solution_values_neighbor_in[i] *
+                                  fe_face_values_p.shape_value(i, q);
+                    }
+                  solution_gradients[q] = gradients * fe_face_values_m.JxW(q);
+                  solution_values[q]    = values * fe_face_values_m.JxW(q);
+                }
+
+              for (const auto i : fe_face_values_m.dof_indices())
+                {
+                  double sum_gradients = 0.;
+                  double sum_values    = 0.;
+                  for (const auto q :
+                       fe_face_values_m.quadrature_point_indices())
+                    {
+                      sum_gradients += solution_gradients[q] *
+                                       fe_face_values_m.shape_grad(i, q);
+                      sum_values +=
+                        solution_values[q] * fe_face_values_m.shape_value(i, q);
+                    }
+
+                  solution_values_faces_out[i] = sum_gradients + sum_values;
+                }
+
+              cell->distribute_local_to_global(
+                Vector<double>(solution_values_faces_out.begin(),
+                               solution_values_faces_out.end()),
+                dst_faces_2);
+            }
+        }
+    }
+
+  deallog << "check difference l2 norm cell: "
+          << dst_cell.l2_norm() - dst_cell_2.l2_norm() << std::endl;
+
+  deallog << "check difference l2 norm surface: "
+          << dst_surface.l2_norm() - dst_surface_2.l2_norm() << std::endl;
+
+  deallog << "check difference l2 norm faces: "
+          << dst_faces.l2_norm() - dst_faces_2.l2_norm() << std::endl;
+}
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1);
+
+  initlog();
+
+  test(true);
+  deallog << std::endl;
+  test(false);
+}

diff --git a/tests/non_matching/mapping_info.output b/tests/non_matching/mapping_info.output
new file mode 100644 (file)
index 0000000..b473752
--- /dev/null
+++ b/tests/non_matching/mapping_info.output
@@ -0,0 +1,8 @@
+
+DEAL::check difference l2 norm cell: 0.00000
+DEAL::check difference l2 norm surface: 0.00000
+DEAL::check difference l2 norm faces: 0.00000
+DEAL::
+DEAL::check difference l2 norm cell: 0.00000
+DEAL::check difference l2 norm surface: 0.00000
+DEAL::check difference l2 norm faces: 0.00000