Enable refreshing MatrixFree geometry data for changed mapping

diff --git a/include/deal.II/matrix_free/mapping_info.h b/include/deal.II/matrix_free/mapping_info.h
index 759b5870dfa7afb5993c76db818e975fc9ab4b41..6d436321949e44eb80de39def9a678078e0e4fa1 100644 (file)
--- a/include/deal.II/matrix_free/mapping_info.h
+++ b/include/deal.II/matrix_free/mapping_info.h
@@ -145,6 +145,12 @@ namespace internal
          */
         unsigned int n_q_points;
 
+        /**
+         * Original one-dimensional quadrature formula applied on the given
+         * cell or face.
+         */
+        Quadrature<1> quadrature_1d;
+
         /**
          * Quadrature formula applied on the given cell or face.
          */
@@ -265,6 +271,12 @@ namespace internal
        */
       AlignedVector<Point<spacedim, VectorizedArrayType>> quadrature_points;
 
+      /**
+       * Clears all data fields except the descriptor vector.
+       */
+      void
+      clear_data_fields();
+
       /**
        * Returns the quadrature index for a given number of quadrature
        * points. If not in hp mode or if the index is not found, this
@@ -281,7 +293,7 @@ namespace internal
       template <typename StreamType>
       void
       print_memory_consumption(StreamType &    out,
-                               const SizeInfo &task_info) const;
+                               const TaskInfo &task_info) const;
 
       /**
        * Returns the memory consumption in bytes.
@@ -324,6 +336,20 @@ namespace internal
         const UpdateFlags update_flags_inner_faces,
         const UpdateFlags update_flags_faces_by_cells);
 
+      /**
+       * Update the information in the given cells and faces that is the
+       * result of a change in the given `mapping` class, keeping the cells,
+       * quadrature formulas and other unknowns unchanged. This call is only
+       * valid if MappingInfo::initialize() has been called before.
+       */
+      void
+      update_mapping(
+        const dealii::Triangulation<dim> &                        tria,
+        const std::vector<std::pair<unsigned int, unsigned int>> &cells,
+        const FaceInfo<VectorizedArrayType::n_array_elements> &   faces,
+        const std::vector<unsigned int> &active_fe_index,
+        const Mapping<dim> &             mapping);
+
       /**
        * Return the type of a given cell as detected during initialization.
        */
@@ -351,6 +377,29 @@ namespace internal
       print_memory_consumption(StreamType &    out,
                                const TaskInfo &task_info) const;
 
+      /**
+       * The given update flags for computing the geometry on the cells.
+       */
+      UpdateFlags update_flags_cells;
+
+      /**
+       * The given update flags for computing the geometry on the boundary
+       * faces.
+       */
+      UpdateFlags update_flags_boundary_faces;
+
+      /**
+       * The given update flags for computing the geometry on the interior
+       * faces.
+       */
+      UpdateFlags update_flags_inner_faces;
+
+      /**
+       * The given update flags for computing the geometry on the faces for
+       * cell-centric loops.
+       */
+      UpdateFlags update_flags_faces_by_cells;
+
       /**
        * Stores whether a cell is Cartesian (cell type 0), has constant
        * transform data (Jacobians) (cell type 1), or is general (cell type
@@ -400,10 +449,8 @@ namespace internal
       initialize_cells(
         const dealii::Triangulation<dim> &                        tria,
         const std::vector<std::pair<unsigned int, unsigned int>> &cells,
-        const std::vector<unsigned int> &              active_fe_index,
-        const Mapping<dim> &                           mapping,
-        const std::vector<dealii::hp::QCollection<1>> &quad,
-        const UpdateFlags                              update_flags_cells);
+        const std::vector<unsigned int> &active_fe_index,
+        const Mapping<dim> &             mapping);
 
       /**
        * Computes the information in the given faces, called within
@@ -415,10 +462,7 @@ namespace internal
         const std::vector<std::pair<unsigned int, unsigned int>> &cells,
         const std::vector<
           FaceToCellTopology<VectorizedArrayType::n_array_elements>> &faces,
-        const Mapping<dim> &                                          mapping,
-        const std::vector<dealii::hp::QCollection<1>> &               quad,
-        const UpdateFlags update_flags_boundary_faces,
-        const UpdateFlags update_flags_inner_faces);
+        const Mapping<dim> &                                          mapping);
 
       /**
        * Computes the information in the given faces, called within
@@ -428,9 +472,7 @@ namespace internal
       initialize_faces_by_cells(
         const dealii::Triangulation<dim> &                        tria,
         const std::vector<std::pair<unsigned int, unsigned int>> &cells,
-        const Mapping<dim> &                                      mapping,
-        const std::vector<dealii::hp::QCollection<1>> &           quad,
-        const UpdateFlags update_flags_faces_by_cells);
+        const Mapping<dim> &                                      mapping);
 
       /**
        * Helper function to determine which update flags must be set in the

diff --git a/include/deal.II/matrix_free/mapping_info.templates.h b/include/deal.II/matrix_free/mapping_info.templates.h
index 4074240891f14f2e5b6511bdf6b54549356d23f1..6b192976e483e2c5bfda0226b2794515f15f82bd 100644 (file)
--- a/include/deal.II/matrix_free/mapping_info.templates.h
+++ b/include/deal.II/matrix_free/mapping_info.templates.h
@@ -63,8 +63,9 @@ namespace internal
       Assert(structdim + 1 <= spacedim ||
                update_flags_inner_faces == update_default,
              ExcMessage("Volume cells do not allow for setting inner faces"));
-      quadrature = Quadrature<structdim>(quadrature_1d);
-      n_q_points = quadrature.size();
+      this->quadrature_1d = quadrature_1d;
+      quadrature          = Quadrature<structdim>(quadrature_1d);
+      n_q_points          = quadrature.size();
       quadrature_weights.resize(n_q_points);
       for (unsigned int i = 0; i < n_q_points; ++i)
         quadrature_weights[i] = quadrature.weight(i);
@@ -125,6 +126,29 @@ namespace internal
 
 
 
+    template <int structdim,
+              int spacedim,
+              typename Number,
+              typename VectorizedArrayType>
+    void
+    MappingInfoStorage<structdim, spacedim, Number, VectorizedArrayType>::
+      clear_data_fields()
+    {
+      data_index_offsets.clear();
+      JxW_values.clear();
+      normal_vectors.clear();
+      for (unsigned int i = 0; i < 2; ++i)
+        {
+          jacobians[i].clear();
+          jacobian_gradients[i].clear();
+          normals_times_jacobians[i].clear();
+        }
+      quadrature_point_offsets.clear();
+      quadrature_points.clear();
+    }
+
+
+
     template <int structdim,
               int spacedim,
               typename Number,
@@ -156,7 +180,7 @@ namespace internal
     template <typename StreamType>
     void
     MappingInfoStorage<structdim, spacedim, Number, VectorizedArrayType>::
-      print_memory_consumption(StreamType &out, const SizeInfo &task_info) const
+      print_memory_consumption(StreamType &out, const TaskInfo &task_info) const
     {
       // print_memory_statistics involves global communication, so we can
       // disable the check here only if no processor has any such data
@@ -287,19 +311,79 @@ namespace internal
       clear();
       this->mapping = &mapping;
 
+      cell_data.resize(quad.size());
+      face_data.resize(quad.size());
+      face_data_by_cells.resize(quad.size());
+
+      // dummy FE that is used to set up an FEValues object. Do not need the
+      // actual finite element because we will only evaluate quantities for
+      // the mapping that are independent of the FE
+      this->update_flags_cells = compute_update_flags(update_flags_cells, quad);
+
+      this->update_flags_boundary_faces =
+        ((update_flags_inner_faces | update_flags_boundary_faces) &
+             update_quadrature_points ?
+           update_quadrature_points :
+           update_default) |
+        update_normal_vectors | update_JxW_values | update_jacobians;
+      this->update_flags_inner_faces    = this->update_flags_boundary_faces;
+      this->update_flags_faces_by_cells = update_flags_faces_by_cells;
+
+      for (unsigned int my_q = 0; my_q < quad.size(); ++my_q)
+        {
+          const unsigned int n_hp_quads = quad[my_q].size();
+          AssertIndexRange(0, n_hp_quads);
+          cell_data[my_q].descriptor.resize(n_hp_quads);
+          for (unsigned int q = 0; q < n_hp_quads; ++q)
+            cell_data[my_q].descriptor[q].initialize(quad[my_q][q],
+                                                     update_default);
+
+          face_data[my_q].descriptor.resize(n_hp_quads);
+          for (unsigned int hpq = 0; hpq < n_hp_quads; ++hpq)
+            face_data[my_q].descriptor[hpq].initialize(
+              quad[my_q][hpq], update_flags_boundary_faces);
+
+          face_data_by_cells[my_q].descriptor.resize(n_hp_quads);
+          for (unsigned int hpq = 0; hpq < n_hp_quads; ++hpq)
+            face_data_by_cells[my_q].descriptor[hpq].initialize(quad[my_q][hpq],
+                                                                update_default);
+        }
+
       // Could call these functions in parallel, but not useful because the
       // work inside is nicely split up already
-      initialize_cells(
-        tria, cells, active_fe_index, mapping, quad, update_flags_cells);
-      initialize_faces(tria,
-                       cells,
-                       face_info.faces,
-                       mapping,
-                       quad,
-                       update_flags_boundary_faces,
-                       update_flags_inner_faces);
-      initialize_faces_by_cells(
-        tria, cells, mapping, quad, update_flags_faces_by_cells);
+      initialize_cells(tria, cells, active_fe_index, mapping);
+      initialize_faces(tria, cells, face_info.faces, mapping);
+      initialize_faces_by_cells(tria, cells, mapping);
+    }
+
+
+
+    template <int dim, typename Number, typename VectorizedArrayType>
+    void
+    MappingInfo<dim, Number, VectorizedArrayType>::update_mapping(
+      const dealii::Triangulation<dim> &                        tria,
+      const std::vector<std::pair<unsigned int, unsigned int>> &cells,
+      const FaceInfo<VectorizedArrayType::n_array_elements> &   face_info,
+      const std::vector<unsigned int> &                         active_fe_index,
+      const Mapping<dim> &                                      mapping)
+    {
+      AssertDimension(cells.size() / VectorizedArrayType::n_array_elements,
+                      cell_type.size());
+
+      for (auto &data : cell_data)
+        data.clear_data_fields();
+      for (auto &data : face_data)
+        data.clear_data_fields();
+      for (auto &data : face_data_by_cells)
+        data.clear_data_fields();
+
+      this->mapping = &mapping;
+
+      // Could call these functions in parallel, but not useful because the
+      // work inside is nicely split up already
+      initialize_cells(tria, cells, active_fe_index, mapping);
+      initialize_faces(tria, cells, face_info.faces, mapping);
+      initialize_faces_by_cells(tria, cells, mapping);
     }
 
 
@@ -583,8 +667,6 @@ namespace internal
         const std::vector<std::pair<unsigned int, unsigned int>> &cells,
         const std::vector<unsigned int> &              active_fe_index,
         const Mapping<dim> &                           mapping,
-        const std::vector<dealii::hp::QCollection<1>> &quad,
-        const UpdateFlags                              update_flags,
         MappingInfo<dim, Number, VectorizedArrayType> &mapping_info,
         std::pair<std::vector<
                     MappingInfoStorage<dim, dim, Number, VectorizedArrayType>>,
@@ -626,7 +708,8 @@ namespace internal
           fe_values(mapping_info.cell_data.size());
         for (unsigned int i = 0; i < fe_values.size(); ++i)
           fe_values[i].resize(mapping_info.cell_data[i].descriptor.size());
-        UpdateFlags update_flags_feval =
+        const UpdateFlags update_flags = mapping_info.update_flags_cells;
+        const UpdateFlags update_flags_feval =
           (update_flags & update_jacobians ? update_jacobians :
                                              update_default) |
           (update_flags & update_jacobian_grads ? update_jacobian_grads :
@@ -634,15 +717,18 @@ namespace internal
           (update_flags & update_quadrature_points ? update_quadrature_points :
                                                      update_default);
 
-        std::vector<std::vector<unsigned int>> n_q_points_1d(quad.size()),
-          step_size_cartesian(quad.size());
-        for (unsigned int my_q = 0; my_q < quad.size(); ++my_q)
+        std::vector<std::vector<unsigned int>> n_q_points_1d(fe_values.size()),
+          step_size_cartesian(fe_values.size());
+        for (unsigned int my_q = 0; my_q < fe_values.size(); ++my_q)
           {
-            n_q_points_1d[my_q].resize(quad[my_q].size());
-            step_size_cartesian[my_q].resize(quad[my_q].size());
-            for (unsigned int hpq = 0; hpq < quad[my_q].size(); ++hpq)
+            n_q_points_1d[my_q].resize(
+              mapping_info.cell_data[my_q].descriptor.size());
+            step_size_cartesian[my_q].resize(n_q_points_1d[my_q].size());
+            for (unsigned int hpq = 0; hpq < n_q_points_1d[my_q].size(); ++hpq)
               {
-                n_q_points_1d[my_q][hpq] = quad[my_q][hpq].size();
+                n_q_points_1d[my_q][hpq] = mapping_info.cell_data[my_q]
+                                             .descriptor[hpq]
+                                             .quadrature_1d.size();
 
                 // To walk on the diagonal for lexicographic ordering, we have
                 // to jump one index ahead in each direction. For direction 0,
@@ -1024,34 +1110,15 @@ namespace internal
       const dealii::Triangulation<dim> &                        tria,
       const std::vector<std::pair<unsigned int, unsigned int>> &cells,
       const std::vector<unsigned int> &                         active_fe_index,
-      const Mapping<dim> &                                      mapping,
-      const std::vector<dealii::hp::QCollection<1>> &           quad,
-      const UpdateFlags update_flags_input)
+      const Mapping<dim> &                                      mapping)
     {
-      const unsigned int n_quads = quad.size();
       const unsigned int n_cells = cells.size();
       const unsigned int vectorization_width =
         VectorizedArrayType::n_array_elements;
       Assert(n_cells % vectorization_width == 0, ExcInternalError());
       const unsigned int n_macro_cells = n_cells / vectorization_width;
-      cell_data.resize(n_quads);
       cell_type.resize(n_macro_cells);
 
-      // dummy FE that is used to set up an FEValues object. Do not need the
-      // actual finite element because we will only evaluate quantities for
-      // the mapping that are independent of the FE
-      UpdateFlags update_flags = compute_update_flags(update_flags_input, quad);
-
-      for (unsigned int my_q = 0; my_q < n_quads; ++my_q)
-        {
-          const unsigned int n_hp_quads = quad[my_q].size();
-          AssertIndexRange(0, n_hp_quads);
-          cell_data[my_q].descriptor.resize(n_hp_quads);
-          for (unsigned int q = 0; q < n_hp_quads; ++q)
-            cell_data[my_q].descriptor[q].initialize(quad[my_q][q],
-                                                     update_default);
-        }
-
       if (n_macro_cells == 0)
         return;
 
@@ -1078,7 +1145,7 @@ namespace internal
             data_cells_local.push_back(std::make_pair(
               std::vector<
                 MappingInfoStorage<dim, dim, Number, VectorizedArrayType>>(
-                n_quads),
+                cell_data.size()),
               ExtractCellHelper::
                 CompressedCellData<dim, Number, VectorizedArrayType>(
                   ExtractCellHelper::get_jacobian_size(tria))));
@@ -1090,8 +1157,6 @@ namespace internal
               cells,
               active_fe_index,
               mapping,
-              quad,
-              update_flags,
               *this,
               data_cells_local.back());
             cell_range.first = cell_range.second;
@@ -1135,10 +1200,10 @@ namespace internal
             data_cells_local.back().first[my_q].JxW_values.size());
           cell_data[my_q].jacobians[0].resize_fast(
             cell_data[my_q].JxW_values.size());
-          if (update_flags & update_jacobian_grads)
+          if (update_flags_cells & update_jacobian_grads)
             cell_data[my_q].jacobian_gradients[0].resize_fast(
               cell_data[my_q].JxW_values.size());
-          if (update_flags & update_quadrature_points)
+          if (update_flags_cells & update_quadrature_points)
             {
               cell_data[my_q].quadrature_point_offsets.resize(cell_type.size());
               cell_data[my_q].quadrature_points.resize_fast(
@@ -1276,8 +1341,6 @@ namespace internal
         const std::vector<
           FaceToCellTopology<VectorizedArrayType::n_array_elements>> &faces,
         const Mapping<dim> &                                          mapping,
-        const UpdateFlags                              update_flags_boundary,
-        const UpdateFlags                              update_flags_inner,
         MappingInfo<dim, Number, VectorizedArrayType> &mapping_info,
         std::pair<
           std::vector<
@@ -1328,16 +1391,18 @@ namespace internal
               if (is_boundary_face &&
                   fe_boundary_face_values_container[my_q][0] == nullptr)
                 fe_boundary_face_values_container[my_q][0] =
-                  std::make_shared<FEFaceValues<dim>>(mapping,
-                                                      dummy_fe,
-                                                      quadrature,
-                                                      update_flags_boundary);
+                  std::make_shared<FEFaceValues<dim>>(
+                    mapping,
+                    dummy_fe,
+                    quadrature,
+                    mapping_info.update_flags_boundary_faces);
               else if (fe_face_values_container[my_q][0] == nullptr)
                 fe_face_values_container[my_q][0] =
-                  std::make_shared<FEFaceValues<dim>>(mapping,
-                                                      dummy_fe,
-                                                      quadrature,
-                                                      update_flags_inner);
+                  std::make_shared<FEFaceValues<dim>>(
+                    mapping,
+                    dummy_fe,
+                    quadrature,
+                    mapping_info.update_flags_inner_faces);
 
               FEFaceValues<dim> &fe_face_values =
                 is_boundary_face ? *fe_boundary_face_values_container[my_q][0] :
@@ -1489,7 +1554,7 @@ namespace internal
                                 mapping,
                                 dummy_fe,
                                 quadrature,
-                                update_flags_inner);
+                                mapping_info.update_flags_inner_faces);
                           fe_subface_values_container[my_q][0]->reinit(
                             cell_it,
                             faces[face].exterior_face_no,
@@ -1700,39 +1765,9 @@ namespace internal
       const std::vector<std::pair<unsigned int, unsigned int>> &cells,
       const std::vector<
         FaceToCellTopology<VectorizedArrayType::n_array_elements>> &faces,
-      const Mapping<dim> &                                          mapping,
-      const std::vector<dealii::hp::QCollection<1>> &               quad,
-      const UpdateFlags update_flags_boundary_faces,
-      const UpdateFlags update_flags_inner_faces)
+      const Mapping<dim> &                                          mapping)
     {
       face_type.resize(faces.size(), general);
-      face_data.resize(quad.size());
-
-      // We currently always set the same flags on both inner and boundary
-      // faces to the same value. At some point, we might want to separate the
-      // two.
-      UpdateFlags update_flags_compute_boundary =
-        ((update_flags_inner_faces | update_flags_boundary_faces) &
-             update_quadrature_points ?
-           update_quadrature_points :
-           update_default) |
-        update_normal_vectors | update_JxW_values | update_jacobians;
-      UpdateFlags update_flags_compute_inner =
-        ((update_flags_inner_faces | update_flags_boundary_faces) &
-             update_quadrature_points ?
-           update_quadrature_points :
-           update_default) |
-        update_normal_vectors | update_JxW_values | update_jacobians;
-      UpdateFlags update_flags_common =
-        update_flags_inner_faces | update_flags_boundary_faces;
-
-      for (unsigned int my_q = 0; my_q < quad.size(); ++my_q)
-        {
-          face_data[my_q].descriptor.resize(quad[my_q].size());
-          for (unsigned int hpq = 0; hpq < quad[my_q].size(); ++hpq)
-            face_data[my_q].descriptor[hpq].initialize(
-              quad[my_q][hpq], update_flags_compute_inner);
-        }
 
       if (faces.size() == 0)
         return;
@@ -1761,7 +1796,7 @@ namespace internal
             data_faces_local.push_back(std::make_pair(
               std::vector<
                 MappingInfoStorage<dim - 1, dim, Number, VectorizedArrayType>>(
-                quad.size()),
+                face_data.size()),
               ExtractFaceHelper::
                 CompressedFaceData<dim, Number, VectorizedArrayType>(
                   ExtractCellHelper::get_jacobian_size(tria))));
@@ -1773,8 +1808,6 @@ namespace internal
               cells,
               faces,
               mapping,
-              update_flags_compute_boundary,
-              update_flags_compute_inner,
               *this,
               data_faces_local.back());
             face_range.first = face_range.second;
@@ -1794,6 +1827,9 @@ namespace internal
           data_faces_local[0].second.data,
           indices_compressed[i]);
 
+      const UpdateFlags update_flags_common =
+        update_flags_boundary_faces | update_flags_inner_faces;
+
       // Collect all data in the final data fields.
       // First allocate the memory
       const unsigned int n_constant_data =
@@ -1926,15 +1962,12 @@ namespace internal
     MappingInfo<dim, Number, VectorizedArrayType>::initialize_faces_by_cells(
       const dealii::Triangulation<dim> &                        tria,
       const std::vector<std::pair<unsigned int, unsigned int>> &cells,
-      const Mapping<dim> &                                      mapping,
-      const std::vector<dealii::hp::QCollection<1>> &           quad,
-      const UpdateFlags update_flags_faces_by_cells)
+      const Mapping<dim> &                                      mapping)
     {
       if (update_flags_faces_by_cells == update_default)
         return;
 
-      face_data_by_cells.resize(quad.size());
-      const unsigned int n_quads = quad.size();
+      const unsigned int n_quads = face_data_by_cells.size();
       const unsigned int vectorization_width =
         VectorizedArrayType::n_array_elements;
       UpdateFlags update_flags =
@@ -1945,13 +1978,6 @@ namespace internal
 
       for (unsigned int my_q = 0; my_q < n_quads; ++my_q)
         {
-          const unsigned int n_hp_quads = quad[my_q].size();
-          AssertIndexRange(0, n_hp_quads);
-          face_data_by_cells[my_q].descriptor.resize(n_hp_quads);
-          for (unsigned int q = 0; q < n_hp_quads; ++q)
-            face_data_by_cells[my_q].descriptor[q].initialize(quad[my_q][q],
-                                                              update_default);
-
           // since we already know the cell type, we can pre-allocate the right
           // amount of data straight away and we just need to do some basic
           // counting
@@ -2217,7 +2243,7 @@ namespace internal
     void
     MappingInfo<dim, Number, VectorizedArrayType>::print_memory_consumption(
       StreamType &    out,
-      const SizeInfo &task_info) const
+      const TaskInfo &task_info) const
     {
       out << "    Cell types:                      ";
       task_info.print_memory_statistics(out,

diff --git a/include/deal.II/matrix_free/matrix_free.h b/include/deal.II/matrix_free/matrix_free.h
index 107f680371780adb4593d9f7db8e0cc0a853cb5a..32b8ca42b486c327709b3fcd668755a2438695bf 100644 (file)
--- a/include/deal.II/matrix_free/matrix_free.h
+++ b/include/deal.II/matrix_free/matrix_free.h
@@ -684,6 +684,18 @@ public:
   copy_from(
     const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free_base);
 
+  /**
+   * Refreshes the geometry data stored in the MappingInfo fields when the
+   * underlying geometry has changed (e.g. by a mapping that can deform
+   * through a change in the spatial configuration like MappingFEField)
+   * whereas the topology of the mesh and unknowns have remained the
+   * same. Compared to reinit(), this operation only has to re-generate the
+   * geometry arrays and can thus be significantly cheaper (depending on the
+   * cost to evaluate the geometry).
+   */
+  void
+  update_mapping(const Mapping<dim> &mapping);
+
   /**
    * Clear all data fields and brings the class into a condition similar to
    * after having called the default constructor.

diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h
index 7cd7dd18c085a77b777a5481a83dbd53378f31a9..21b77dbe1a9a0b7e99dd9a0fbe4959ca8a34976d 100644 (file)
--- a/include/deal.II/matrix_free/matrix_free.templates.h
+++ b/include/deal.II/matrix_free/matrix_free.templates.h
@@ -681,6 +681,25 @@ MatrixFree<dim, Number, VectorizedArrayType>::internal_reinit(
 }
 
 
+
+template <int dim, typename Number, typename VectorizedArrayType>
+void
+MatrixFree<dim, Number, VectorizedArrayType>::update_mapping(
+  const Mapping<dim> &mapping)
+{
+  AssertDimension(shape_info.size(1), mapping_info.cell_data.size());
+  mapping_info.update_mapping(
+    dof_handlers.active_dof_handler == DoFHandlers::hp ?
+      dof_handlers.hp_dof_handler[0]->get_triangulation() :
+      dof_handlers.dof_handler[0]->get_triangulation(),
+    cell_level_index,
+    face_info,
+    dof_info[0].cell_active_fe_index,
+    mapping);
+}
+
+
+
 template <int dim, typename Number, typename VectorizedArrayType>
 template <int spacedim>
 bool