]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Optimize computations in MappingQ 15035/head
authorMaximilian Bergbauer <maximilian.bergbauer@tum.de>
Wed, 5 Apr 2023 15:26:44 +0000 (17:26 +0200)
committerMaximilian Bergbauer <maximilian.bergbauer@tum.de>
Tue, 18 Apr 2023 13:00:32 +0000 (15:00 +0200)
include/deal.II/base/qprojector.h
include/deal.II/fe/mapping_q.h
include/deal.II/fe/mapping_q_internal.h
include/deal.II/non_matching/mapping_info.h
source/base/qprojector.cc
source/fe/mapping_q.cc

index 381f8aa0fac423b9c06696aaf87a7f70e52e1648..4e4749752b14edd06ade19aa6273af37290a423d 100644 (file)
@@ -109,6 +109,20 @@ public:
                   const SubQuadrature &quadrature,
                   const unsigned int   face_no);
 
+  /**
+   * Compute the cell quadrature formula corresponding to using
+   * <tt>quadrature</tt> on face <tt>face_no</tt> taking into account the
+   * orientation of the face. For further details, see the general doc for this
+   * class.
+   */
+  static Quadrature<dim>
+  project_to_oriented_face(const ReferenceCell &reference_cell,
+                           const SubQuadrature &quadrature,
+                           const unsigned int   face_no,
+                           const bool           face_orientation,
+                           const bool           face_flip,
+                           const bool           face_rotation);
+
   /**
    * Compute the quadrature points on the cell if the given quadrature formula
    * is used on face <tt>face_no</tt>, subface number <tt>subface_no</tt>
@@ -148,6 +162,25 @@ public:
                      const RefinementCase<dim - 1> &ref_case =
                        RefinementCase<dim - 1>::isotropic_refinement);
 
+  /**
+   * Compute the cell quadrature formula corresponding to using
+   * <tt>quadrature</tt> on subface <tt>subface_no</tt> of face
+   * <tt>face_no</tt> with SubfaceCase<dim> <tt>ref_case</tt>. The last
+   * argument is only used in 3d.
+   *
+   * @note Only the points are transformed. The quadrature weights are the
+   * same as those of the original rule.
+   */
+  static Quadrature<dim>
+  project_to_oriented_subface(const ReferenceCell &            reference_cell,
+                              const SubQuadrature &            quadrature,
+                              const unsigned int               face_no,
+                              const unsigned int               subface_no,
+                              const bool                       face_orientation,
+                              const bool                       face_flip,
+                              const bool                       face_rotation,
+                              const internal::SubfaceCase<dim> ref_case);
+
   /**
    * Take a collection of face quadrature formulas and generate a cell
    * quadrature formula from it where the quadrature points of the given
index 98366d47d4d020ee29aaaaf53dd9b9053c27b9eb..28f34de1bf96fd46865421311a0fa6bf48534b05 100644 (file)
@@ -255,6 +255,37 @@ public:
     internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
       &output_data) const;
 
+  /**
+   * As opposed to the fill_fe_face_values()
+   * function that relies on pre-computed information of InternalDataBase, this
+   * function chooses the flexible evaluation path on the cell and points
+   * passed in to the current function.
+   *
+   * @param[in] cell The cell where to evaluate the mapping.
+   *
+   * @param[in] face_number The face number where to evaluate the mapping.
+   *
+   * @param[in] face_quadrature The quadrature points where the
+   * transformation (Jacobians, positions) should be computed.
+   *
+   * @param[in] internal_data A reference to an object previously created
+   * that may be used to store information the mapping can compute once on the
+   * reference cell. See the documentation of the Mapping::InternalDataBase
+   * class for an extensive description of the purpose of these objects.
+   *
+   * @param[out] output_data A struct containing the evaluated quantities such
+   * as the Jacobian resulting from application of the mapping on the given
+   * cell with its underlying manifolds.
+   */
+  void
+  fill_mapping_data_for_face_quadrature(
+    const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+    const unsigned int                                          face_number,
+    const Quadrature<dim - 1> &                                 face_quadrature,
+    const typename Mapping<dim, spacedim>::InternalDataBase &   internal_data,
+    internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
+      &output_data) const;
+
   /**
    * @name Interface with FEValues and friends
    * @{
@@ -391,6 +422,13 @@ public:
     virtual std::size_t
     memory_consumption() const override;
 
+    /**
+     * Location of quadrature points of faces or subfaces in 3d with all
+     * possible orientations. Can be accessed with the correct offset provided
+     * via QProjector::DataSetDescriptor. Not needed/used for cells.
+     */
+    AlignedVector<Point<dim>> quadrature_points;
+
     /**
      * Values of shape functions. Access by function @p shape.
      *
index ddd4c145b2bc15ac100d9e8b83e90e46d1014130..e144e330966f93a7221b6ae751cebab709e9ec44 100644 (file)
@@ -1116,6 +1116,8 @@ namespace internal
       const CellSimilarity::Similarity cell_similarity,
       const typename dealii::MappingQ<dim, spacedim>::InternalData &data,
       std::vector<Point<spacedim>> &                 quadrature_points,
+      std::vector<DerivativeForm<1, dim, spacedim>> &jacobians,
+      std::vector<DerivativeForm<1, spacedim, dim>> &inverse_jacobians,
       std::vector<DerivativeForm<2, dim, spacedim>> &jacobian_grads)
     {
       const UpdateFlags update_flags = data.update_each;
@@ -1246,6 +1248,26 @@ namespace internal
             data.volume_elements[point] =
               data.contravariant[point].determinant();
 
+      // copy values from InternalData to vector given by reference
+      if (update_flags & update_jacobians)
+        {
+          const unsigned int n_q_points = data.contravariant.size();
+          AssertDimension(jacobians.size(), n_q_points);
+          if (cell_similarity != CellSimilarity::translation)
+            for (unsigned int point = 0; point < n_q_points; ++point)
+              jacobians[point] = data.contravariant[point];
+        }
+
+      // copy values from InternalData to vector given by reference
+      if (update_flags & update_inverse_jacobians)
+        {
+          const unsigned int n_q_points = data.contravariant.size();
+          AssertDimension(inverse_jacobians.size(), n_q_points);
+          if (cell_similarity != CellSimilarity::translation)
+            for (unsigned int point = 0; point < n_q_points; ++point)
+              inverse_jacobians[point] = data.covariant[point].transpose();
+        }
+
       if (evaluation_flag & EvaluationFlags::hessians)
         {
           constexpr int desymmetrize_3d[6][2] = {
@@ -1286,112 +1308,114 @@ namespace internal
     }
 
 
-    /**
-     * Compute the locations of quadrature points on the object described by
-     * the first argument (and the cell for which the mapping support points
-     * have already been set), but only if the update_flags of the @p data
-     * argument indicate so.
-     */
+
     template <int dim, int spacedim>
     inline void
-    maybe_compute_q_points(
-      const typename QProjector<dim>::DataSetDescriptor             data_set,
+    maybe_update_q_points_Jacobians_generic(
+      const CellSimilarity::Similarity cell_similarity,
       const typename dealii::MappingQ<dim, spacedim>::InternalData &data,
-      std::vector<Point<spacedim>> &quadrature_points)
+      const ArrayView<const Point<dim>> &                           unit_points,
+      const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
+      const unsigned int                                  polynomial_degree,
+      const std::vector<unsigned int> &renumber_lexicographic_to_hierarchic,
+      std::vector<Point<spacedim>> &   quadrature_points,
+      std::vector<DerivativeForm<1, dim, spacedim>> &jacobians,
+      std::vector<DerivativeForm<1, spacedim, dim>> &inverse_jacobians)
     {
-      const UpdateFlags update_flags = data.update_each;
+      const UpdateFlags                   update_flags = data.update_each;
+      const std::vector<Point<spacedim>> &support_points =
+        data.mapping_support_points;
 
-      if (update_flags & update_quadrature_points)
-        for (unsigned int point = 0; point < quadrature_points.size(); ++point)
+      const unsigned int n_points = unit_points.size();
+      const unsigned int n_lanes  = VectorizedArray<double>::size();
+
+      // Use the more heavy VectorizedArray code path if there is more than
+      // one point left to compute
+      for (unsigned int i = 0; i < n_points; i += n_lanes)
+        if (n_points - i > 1)
           {
-            const double *  shape = &data.shape(point + data_set, 0);
-            Point<spacedim> result =
-              (shape[0] * data.mapping_support_points[0]);
-            for (unsigned int k = 1; k < data.n_shape_functions; ++k)
-              for (unsigned int i = 0; i < spacedim; ++i)
-                result[i] += shape[k] * data.mapping_support_points[k][i];
-            quadrature_points[point] = result;
-          }
-    }
+            Point<dim, VectorizedArray<double>> p_vec;
+            for (unsigned int j = 0; j < n_lanes; ++j)
+              if (i + j < n_points)
+                for (unsigned int d = 0; d < dim; ++d)
+                  p_vec[d][j] = unit_points[i + j][d];
+              else
+                for (unsigned int d = 0; d < dim; ++d)
+                  p_vec[d][j] = unit_points[i][d];
 
+            const auto result =
+              internal::evaluate_tensor_product_value_and_gradient(
+                polynomials_1d,
+                support_points,
+                p_vec,
+                polynomial_degree == 1,
+                renumber_lexicographic_to_hierarchic);
 
+            if (update_flags & update_quadrature_points)
+              for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
+                for (unsigned int d = 0; d < spacedim; ++d)
+                  quadrature_points[i + j][d] = result.first[d][j];
 
-    /**
-     * Update the co- and contravariant matrices as well as their determinant,
-     * for the cell
-     * described stored in the data object, but only if the update_flags of the @p data
-     * argument indicate so.
-     *
-     * Skip the computation if possible as indicated by the first argument.
-     */
-    template <int dim, int spacedim>
-    inline void
-    maybe_update_Jacobians(
-      const CellSimilarity::Similarity                          cell_similarity,
-      const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
-      const typename dealii::MappingQ<dim, spacedim>::InternalData &data)
-    {
-      const UpdateFlags update_flags = data.update_each;
+            if (cell_similarity == CellSimilarity::translation)
+              continue;
 
-      if (update_flags & update_contravariant_transformation)
-        // if the current cell is just a
-        // translation of the previous one, no
-        // need to recompute jacobians...
-        if (cell_similarity != CellSimilarity::translation)
+            if (update_flags & update_contravariant_transformation)
+              for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
+                for (unsigned int d = 0; d < spacedim; ++d)
+                  for (unsigned int e = 0; e < dim; ++e)
+                    data.contravariant[i + j][d][e] = result.second[e][d][j];
+
+            if (update_flags & update_volume_elements)
+              for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
+                data.volume_elements[i + j] =
+                  data.contravariant[i + j].determinant();
+
+            if (update_flags & update_jacobians)
+              for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
+                jacobians[i + j] = data.contravariant[i + j];
+
+            if (update_flags & update_covariant_transformation)
+              for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
+                data.covariant[i + j] =
+                  data.contravariant[i + j].covariant_form();
+
+            if (update_flags & update_inverse_jacobians)
+              for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
+                inverse_jacobians[i + j] = data.covariant[i + j].transpose();
+          }
+        else
           {
-            const unsigned int n_q_points = data.contravariant.size();
+            const auto result =
+              internal::evaluate_tensor_product_value_and_gradient(
+                polynomials_1d,
+                support_points,
+                unit_points[i],
+                polynomial_degree == 1,
+                renumber_lexicographic_to_hierarchic);
 
-            std::fill(data.contravariant.begin(),
-                      data.contravariant.end(),
-                      DerivativeForm<1, dim, spacedim>());
+            if (update_flags & update_quadrature_points)
+              quadrature_points[i] = result.first;
 
-            Assert(data.n_shape_functions > 0, ExcInternalError());
+            if (cell_similarity == CellSimilarity::translation)
+              continue;
 
-            for (unsigned int point = 0; point < n_q_points; ++point)
+            if (update_flags & update_contravariant_transformation)
               {
-                double result[spacedim][dim];
-
-                // peel away part of sum to avoid zeroing the
-                // entries and adding for the first time
-                for (unsigned int i = 0; i < spacedim; ++i)
-                  for (unsigned int j = 0; j < dim; ++j)
-                    result[i][j] = data.derivative(point + data_set, 0)[j] *
-                                   data.mapping_support_points[0][i];
-                for (unsigned int k = 1; k < data.n_shape_functions; ++k)
-                  for (unsigned int i = 0; i < spacedim; ++i)
-                    for (unsigned int j = 0; j < dim; ++j)
-                      result[i][j] += data.derivative(point + data_set, k)[j] *
-                                      data.mapping_support_points[k][i];
-
-                // write result into contravariant data. for
-                // j=dim in the case dim<spacedim, there will
-                // never be any nonzero data that arrives in
-                // here, so it is ok anyway because it was
-                // initialized to zero at the initialization
-                for (unsigned int i = 0; i < spacedim; ++i)
-                  for (unsigned int j = 0; j < dim; ++j)
-                    data.contravariant[point][i][j] = result[i][j];
+                DerivativeForm<1, spacedim, dim> jac_transposed = result.second;
+                data.contravariant[i] = jac_transposed.transpose();
               }
-          }
 
-      if (update_flags & update_covariant_transformation)
-        if (cell_similarity != CellSimilarity::translation)
-          {
-            const unsigned int n_q_points = data.contravariant.size();
-            for (unsigned int point = 0; point < n_q_points; ++point)
-              {
-                data.covariant[point] =
-                  (data.contravariant[point]).covariant_form();
-              }
-          }
+            if (update_flags & update_volume_elements)
+              data.volume_elements[i] = data.contravariant[i].determinant();
 
-      if (update_flags & update_volume_elements)
-        if (cell_similarity != CellSimilarity::translation)
-          {
-            const unsigned int n_q_points = data.contravariant.size();
-            for (unsigned int point = 0; point < n_q_points; ++point)
-              data.volume_elements[point] =
-                data.contravariant[point].determinant();
+            if (update_flags & update_jacobians)
+              jacobians[i] = data.contravariant[i];
+
+            if (update_flags & update_covariant_transformation)
+              data.covariant[i] = data.contravariant[i].covariant_form();
+
+            if (update_flags & update_inverse_jacobians)
+              inverse_jacobians[i] = data.covariant[i].transpose();
           }
     }
 
@@ -2012,19 +2036,11 @@ namespace internal
               output_data.normal_vectors[i] =
                 Point<spacedim>(output_data.boundary_forms[i] /
                                 output_data.boundary_forms[i].norm());
-
-          if (update_flags & update_jacobians)
-            for (unsigned int point = 0; point < n_q_points; ++point)
-              output_data.jacobians[point] = data.contravariant[point];
-
-          if (update_flags & update_inverse_jacobians)
-            for (unsigned int point = 0; point < n_q_points; ++point)
-              output_data.inverse_jacobians[point] =
-                data.covariant[point].transpose();
         }
     }
 
 
+
     /**
      * Do the work of MappingQ::fill_fe_face_values() and
      * MappingQ::fill_fe_subface_values() in a generic way,
@@ -2041,6 +2057,9 @@ namespace internal
       const typename QProjector<dim>::DataSetDescriptor             data_set,
       const Quadrature<dim - 1> &                                   quadrature,
       const typename dealii::MappingQ<dim, spacedim>::InternalData &data,
+      const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
+      const unsigned int                                  polynomial_degree,
+      const std::vector<unsigned int> &renumber_lexicographic_to_hierarchic,
       internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
         &output_data)
     {
@@ -2050,16 +2069,25 @@ namespace internal
             CellSimilarity::none,
             data,
             output_data.quadrature_points,
+            output_data.jacobians,
+            output_data.inverse_jacobians,
             output_data.jacobian_grads);
         }
       else
         {
-          maybe_compute_q_points<dim, spacedim>(data_set,
-                                                data,
-                                                output_data.quadrature_points);
-          maybe_update_Jacobians<dim, spacedim>(CellSimilarity::none,
-                                                data_set,
-                                                data);
+          internal::MappingQImplementation::
+            maybe_update_q_points_Jacobians_generic(
+              CellSimilarity::none,
+              data,
+              make_array_view(&data.quadrature_points[data_set],
+                              &data.quadrature_points[data_set] +
+                                quadrature.size()),
+              polynomials_1d,
+              polynomial_degree,
+              renumber_lexicographic_to_hierarchic,
+              output_data.quadrature_points,
+              output_data.jacobians,
+              output_data.inverse_jacobians);
           maybe_update_jacobian_grads<dim, spacedim>(
             CellSimilarity::none, data_set, data, output_data.jacobian_grads);
         }
index 5dc9d8f6656a622b23050a44691d2cc30e000675..9dc630e15fa34748a6a1e83f2aabf5ef7c92277c 100644 (file)
@@ -212,10 +212,11 @@ namespace NonMatching
      * class.
      */
     void
-    compute_mapping_data_for_generic_points(
+    compute_mapping_data_for_quadrature(
       const typename Triangulation<dim, spacedim>::cell_iterator &cell,
-      const ArrayView<const Point<dim>> &                         unit_points,
-      MappingData &                                               mapping_data);
+      CellSimilarity::Similarity &cell_similarity,
+      const Quadrature<dim> &     quadrature,
+      MappingData &               mapping_data);
 
     /**
      * Compute the mapping related data for the given @p mapping,
@@ -250,6 +251,12 @@ namespace NonMatching
      */
     std::vector<unsigned int> unit_points_index;
 
+    /**
+     * A pointer to the internal data of the underlying mapping.
+     */
+    std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
+      internal_mapping_data;
+
     /**
      * A pointer to the underlying mapping.
      */
@@ -318,6 +325,16 @@ namespace NonMatching
 
     // always save quadrature points for now
     update_flags_mapping |= update_quadrature_points;
+
+    // construct internal_mapping_data for MappingQ to be able to reuse it in
+    // reinit() calls to avoid memory allocations
+    if (const MappingQ<dim, spacedim> *mapping_q =
+          dynamic_cast<const MappingQ<dim, spacedim> *>(&mapping))
+      {
+        internal_mapping_data =
+          std::make_unique<typename MappingQ<dim, spacedim>::InternalData>(
+            mapping_q->get_degree());
+      }
   }
 
 
@@ -328,7 +345,7 @@ namespace NonMatching
     const typename Triangulation<dim, spacedim>::cell_iterator &cell,
     const std::vector<Point<dim>> &                             unit_points_in)
   {
-    reinit(cell, make_array_view(unit_points_in.begin(), unit_points_in.end()));
+    reinit(cell, Quadrature<dim>(unit_points_in));
   }
 
 
@@ -339,16 +356,9 @@ namespace NonMatching
     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;
-    is_reinitialized();
+    reinit(cell,
+           std::vector<Point<dim>>(unit_points_in.begin(),
+                                   unit_points_in.end()));
   }
 
 
@@ -359,16 +369,18 @@ namespace NonMatching
     const typename Triangulation<dim, spacedim>::cell_iterator &cell,
     const Quadrature<dim> &                                     quadrature)
   {
-    const auto &points  = quadrature.get_points();
-    const auto &weights = quadrature.get_weights();
+    unit_points = quadrature.get_points();
 
-    reinit(cell, points);
+    mapping_data.resize(1);
+    CellSimilarity::Similarity cell_similarity =
+      CellSimilarity::Similarity::none;
+    compute_mapping_data_for_quadrature(cell,
+                                        cell_similarity,
+                                        quadrature,
+                                        mapping_data[0]);
 
-    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];
+    state = State::single_cell;
+    is_reinitialized();
   }
 
 
@@ -380,11 +392,34 @@ namespace NonMatching
     const IteratorRange<Iterator> &             cell_iterator_range,
     const std::vector<std::vector<Point<dim>>> &unit_points_vector,
     const unsigned int                          n_unfiltered_cells)
+  {
+    const unsigned int n_cells = unit_points_vector.size();
+    AssertDimension(n_cells,
+                    std::distance(cell_iterator_range.begin(),
+                                  cell_iterator_range.end()));
+
+    std::vector<Quadrature<dim>> quadrature_vector(n_cells);
+    for (unsigned int cell_index = 0; cell_index < n_cells; ++cell_index)
+      quadrature_vector[cell_index] =
+        Quadrature<dim>(quadrature_vector[cell_index].get_points());
+
+    reinit_cells(cell_iterator_range, quadrature_vector, n_unfiltered_cells);
+  }
+
+
+
+  template <int dim, int spacedim>
+  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)
   {
     do_cell_index_compression =
       n_unfiltered_cells != numbers::invalid_unsigned_int;
 
-    const unsigned int n_cells = unit_points_vector.size();
+    const unsigned int n_cells = quadrature_vector.size();
     AssertDimension(n_cells,
                     std::distance(cell_iterator_range.begin(),
                                   cell_iterator_range.end()));
@@ -392,9 +427,8 @@ namespace NonMatching
     // 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());
+    for (const auto &quadrature : quadrature_vector)
+      unit_points_index.push_back(unit_points_index.back() + quadrature.size());
 
     const unsigned int n_unit_points = unit_points_index.back();
 
@@ -404,19 +438,23 @@ namespace NonMatching
     if (do_cell_index_compression)
       cell_index_to_compressed_cell_index.resize(n_unfiltered_cells,
                                                  numbers::invalid_unsigned_int);
+    CellSimilarity::Similarity cell_similarity =
+      CellSimilarity::Similarity::none;
     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])
+        for (const auto &unit_point :
+             quadrature_vector[cell_index].get_points())
           {
             *it = unit_point;
             ++it;
           }
 
-        compute_mapping_data_for_generic_points(cell,
-                                                unit_points_vector[cell_index],
-                                                mapping_data[cell_index]);
+        compute_mapping_data_for_quadrature(cell,
+                                            cell_similarity,
+                                            quadrature_vector[cell_index],
+                                            mapping_data[cell_index]);
 
         if (do_cell_index_compression)
           cell_index_to_compressed_cell_index[cell->active_cell_index()] =
@@ -431,41 +469,6 @@ namespace NonMatching
 
 
 
-  template <int dim, int spacedim>
-  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)
-  {
-    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>
   template <typename Iterator>
   void
@@ -809,47 +812,35 @@ namespace NonMatching
 
   template <int dim, int spacedim>
   void
-  MappingInfo<dim, spacedim>::compute_mapping_data_for_generic_points(
+  MappingInfo<dim, spacedim>::compute_mapping_data_for_quadrature(
     const typename Triangulation<dim, spacedim>::cell_iterator &cell,
-    const ArrayView<const Point<dim>> &                         unit_points,
+    CellSimilarity::Similarity &                                cell_similarity,
+    const Quadrature<dim> &                                     quadrature,
     MappingData &                                               mapping_data)
   {
+    update_flags_mapping |=
+      mapping->requires_update_flags(update_flags_mapping);
+
+    mapping_data.initialize(quadrature.size(), update_flags_mapping);
+
+    // reuse internal_mapping_data for MappingQ to avoid memory allocations
     if (const MappingQ<dim, spacedim> *mapping_q =
           dynamic_cast<const MappingQ<dim, spacedim> *>(&(*mapping)))
       {
-        mapping_q->fill_mapping_data_for_generic_points(cell,
-                                                        unit_points,
-                                                        update_flags_mapping,
-                                                        mapping_data);
-      }
-    else if (const MappingCartesian<dim, spacedim> *mapping_cartesian =
-               dynamic_cast<const MappingCartesian<dim, spacedim> *>(
-                 &(*mapping)))
-      {
-        mapping_cartesian->fill_mapping_data_for_generic_points(
-          cell, unit_points, update_flags_mapping, mapping_data);
+        (void)mapping_q;
+        auto &data =
+          dynamic_cast<typename MappingQ<dim, spacedim>::InternalData &>(
+            *internal_mapping_data);
+        data.initialize(update_flags_mapping, quadrature, quadrature.size());
       }
     else
       {
-        FE_DGQ<dim, spacedim>           dummy_fe(1);
-        dealii::FEValues<dim, spacedim> fe_values(
-          *mapping,
-          dummy_fe,
-          Quadrature<dim>(
-            std::vector<Point<dim>>(unit_points.begin(), unit_points.end())),
-          update_flags_mapping);
-        fe_values.reinit(cell);
-        mapping_data.initialize(unit_points.size(), update_flags_mapping);
-        if (update_flags_mapping & update_jacobians)
-          for (unsigned int q = 0; q < unit_points.size(); ++q)
-            mapping_data.jacobians[q] = fe_values.jacobian(q);
-        if (update_flags_mapping & update_inverse_jacobians)
-          for (unsigned int q = 0; q < unit_points.size(); ++q)
-            mapping_data.inverse_jacobians[q] = fe_values.inverse_jacobian(q);
-        if (update_flags_mapping & update_quadrature_points)
-          for (unsigned int q = 0; q < unit_points.size(); ++q)
-            mapping_data.quadrature_points[q] = fe_values.quadrature_point(q);
+        internal_mapping_data =
+          mapping->get_data(update_flags_mapping, quadrature);
       }
+
+    cell_similarity = mapping->fill_fe_values(
+      cell, cell_similarity, quadrature, *internal_mapping_data, mapping_data);
   }
 
 
@@ -865,11 +856,23 @@ namespace NonMatching
     update_flags_mapping |=
       mapping->requires_update_flags(update_flags_mapping);
 
-    mapping_data.initialize(quadrature.get_points().size(),
-                            update_flags_mapping);
+    mapping_data.initialize(quadrature.size(), update_flags_mapping);
 
-    auto internal_mapping_data =
-      mapping->get_data(update_flags_mapping, quadrature);
+    // reuse internal_mapping_data for MappingQ to avoid memory allocations
+    if (const MappingQ<dim, spacedim> *mapping_q =
+          dynamic_cast<const MappingQ<dim, spacedim> *>(&(*mapping)))
+      {
+        (void)mapping_q;
+        auto &data =
+          dynamic_cast<typename MappingQ<dim, spacedim>::InternalData &>(
+            *internal_mapping_data);
+        data.initialize(update_flags_mapping, quadrature, quadrature.size());
+      }
+    else
+      {
+        internal_mapping_data =
+          mapping->get_data(update_flags_mapping, quadrature);
+      }
 
     mapping->fill_fe_immersed_surface_values(cell,
                                              quadrature,
@@ -890,18 +893,40 @@ namespace NonMatching
     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_data.initialize(quadrature.size(), update_flags_mapping);
 
-    mapping->fill_fe_face_values(cell,
-                                 face_no,
-                                 hp::QCollection<dim - 1>(quadrature),
-                                 *internal_mapping_data,
-                                 mapping_data);
+    // reuse internal_mapping_data for MappingQ to avoid memory allocations
+    if (const MappingQ<dim, spacedim> *mapping_q =
+          dynamic_cast<const MappingQ<dim, spacedim> *>(&(*mapping)))
+      {
+        auto &data =
+          dynamic_cast<typename MappingQ<dim, spacedim>::InternalData &>(
+            *internal_mapping_data);
+        data.initialize_face(update_flags_mapping,
+                             QProjector<dim>::project_to_oriented_face(
+                               ReferenceCells::get_hypercube<dim>(),
+                               quadrature,
+                               face_no,
+                               cell->face_orientation(face_no),
+                               cell->face_flip(face_no),
+                               cell->face_rotation(face_no)),
+                             quadrature.size());
+
+        mapping_q->fill_mapping_data_for_face_quadrature(
+          cell, face_no, quadrature, *internal_mapping_data, mapping_data);
+      }
+    else
+      {
+        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
 
index 92600865dcdbaaa16100d88dcd9f66d515e2cfa3..e96aed68ecf1650b576cd4592d80e83ef4200bed 100644 (file)
@@ -143,6 +143,193 @@ namespace internal
             q_points.push_back(cell_point);
           }
       }
+
+      std::vector<Point<2>>
+      mutate_points_with_offset(const std::vector<Point<2>> &points,
+                                const unsigned int           offset)
+      {
+        switch (offset)
+          {
+            case 0:
+              return points;
+            case 1:
+            case 2:
+            case 3:
+              return rotate(points, offset);
+            case 4:
+              return reflect(points);
+            case 5:
+            case 6:
+            case 7:
+              return rotate(reflect(points), 8 - offset);
+            default:
+              Assert(false, ExcInternalError());
+          }
+        return {};
+      }
+
+      Quadrature<2>
+      mutate_quadrature(const Quadrature<2> &quadrature,
+                        const bool           face_orientation,
+                        const bool           face_flip,
+                        const bool           face_rotation)
+      {
+        static const unsigned int offset[2][2][2] = {
+          {{4, 5},   // face_orientation=false; face_flip=false;
+                     // face_rotation=false and true
+           {6, 7}},  // face_orientation=false; face_flip=true;
+                     // face_rotation=false and true
+          {{0, 1},   // face_orientation=true;  face_flip=false;
+                     // face_rotation=false and true
+           {2, 3}}}; // face_orientation=true; face_flip=true;
+                     // face_rotation=false and true
+
+        return Quadrature<2>(
+          mutate_points_with_offset(
+            quadrature.get_points(),
+            offset[face_orientation][face_flip][face_rotation]),
+          quadrature.get_weights());
+      }
+
+      std::pair<unsigned int, RefinementCase<2>>
+      select_subface_no_and_refinement_case(
+        const unsigned int             subface_no,
+        const bool                     face_orientation,
+        const bool                     face_flip,
+        const bool                     face_rotation,
+        const internal::SubfaceCase<3> ref_case)
+      {
+        constexpr int dim = 3;
+        // for each subface of a given FaceRefineCase
+        // there is a corresponding equivalent
+        // subface number of one of the "standard"
+        // RefineCases (cut_x, cut_y, cut_xy). Map
+        // the given values to those equivalent
+        // ones.
+
+        // first, define an invalid number
+        static const unsigned int e = numbers::invalid_unsigned_int;
+
+        static const RefinementCase<dim - 1>
+          equivalent_refine_case[internal::SubfaceCase<dim>::case_isotropic + 1]
+                                [GeometryInfo<3>::max_children_per_face] = {
+                                  // case_none. there should be only
+                                  // invalid values here. However, as
+                                  // this function is also called (in
+                                  // tests) for cells which have no
+                                  // refined faces, use isotropic
+                                  // refinement instead
+                                  {RefinementCase<dim - 1>::cut_xy,
+                                   RefinementCase<dim - 1>::cut_xy,
+                                   RefinementCase<dim - 1>::cut_xy,
+                                   RefinementCase<dim - 1>::cut_xy},
+                                  // case_x
+                                  {RefinementCase<dim - 1>::cut_x,
+                                   RefinementCase<dim - 1>::cut_x,
+                                   RefinementCase<dim - 1>::no_refinement,
+                                   RefinementCase<dim - 1>::no_refinement},
+                                  // case_x1y
+                                  {RefinementCase<dim - 1>::cut_xy,
+                                   RefinementCase<dim - 1>::cut_xy,
+                                   RefinementCase<dim - 1>::cut_x,
+                                   RefinementCase<dim - 1>::no_refinement},
+                                  // case_x2y
+                                  {RefinementCase<dim - 1>::cut_x,
+                                   RefinementCase<dim - 1>::cut_xy,
+                                   RefinementCase<dim - 1>::cut_xy,
+                                   RefinementCase<dim - 1>::no_refinement},
+                                  // case_x1y2y
+                                  {RefinementCase<dim - 1>::cut_xy,
+                                   RefinementCase<dim - 1>::cut_xy,
+                                   RefinementCase<dim - 1>::cut_xy,
+                                   RefinementCase<dim - 1>::cut_xy},
+                                  // case_y
+                                  {RefinementCase<dim - 1>::cut_y,
+                                   RefinementCase<dim - 1>::cut_y,
+                                   RefinementCase<dim - 1>::no_refinement,
+                                   RefinementCase<dim - 1>::no_refinement},
+                                  // case_y1x
+                                  {RefinementCase<dim - 1>::cut_xy,
+                                   RefinementCase<dim - 1>::cut_xy,
+                                   RefinementCase<dim - 1>::cut_y,
+                                   RefinementCase<dim - 1>::no_refinement},
+                                  // case_y2x
+                                  {RefinementCase<dim - 1>::cut_y,
+                                   RefinementCase<dim - 1>::cut_xy,
+                                   RefinementCase<dim - 1>::cut_xy,
+                                   RefinementCase<dim - 1>::no_refinement},
+                                  // case_y1x2x
+                                  {RefinementCase<dim - 1>::cut_xy,
+                                   RefinementCase<dim - 1>::cut_xy,
+                                   RefinementCase<dim - 1>::cut_xy,
+                                   RefinementCase<dim - 1>::cut_xy},
+                                  // case_xy (case_isotropic)
+                                  {RefinementCase<dim - 1>::cut_xy,
+                                   RefinementCase<dim - 1>::cut_xy,
+                                   RefinementCase<dim - 1>::cut_xy,
+                                   RefinementCase<dim - 1>::cut_xy}};
+
+        static const unsigned int
+          equivalent_subface_number[internal::SubfaceCase<dim>::case_isotropic +
+                                    1][GeometryInfo<3>::max_children_per_face] =
+            {// case_none, see above
+             {0, 1, 2, 3},
+             // case_x
+             {0, 1, e, e},
+             // case_x1y
+             {0, 2, 1, e},
+             // case_x2y
+             {0, 1, 3, e},
+             // case_x1y2y
+             {0, 2, 1, 3},
+             // case_y
+             {0, 1, e, e},
+             // case_y1x
+             {0, 1, 1, e},
+             // case_y2x
+             {0, 2, 3, e},
+             // case_y1x2x
+             {0, 1, 2, 3},
+             // case_xy (case_isotropic)
+             {0, 1, 2, 3}};
+
+        // If face-orientation or face_rotation are
+        // non-standard, cut_x and cut_y have to be
+        // exchanged.
+        static const RefinementCase<dim - 1> ref_case_permutation[4] = {
+          RefinementCase<dim - 1>::no_refinement,
+          RefinementCase<dim - 1>::cut_y,
+          RefinementCase<dim - 1>::cut_x,
+          RefinementCase<dim - 1>::cut_xy};
+
+        // set a corresponding (equivalent)
+        // RefineCase and subface number
+        const RefinementCase<dim - 1> equ_ref_case =
+          equivalent_refine_case[ref_case][subface_no];
+        const unsigned int equ_subface_no =
+          equivalent_subface_number[ref_case][subface_no];
+        // make sure, that we got a valid subface and RefineCase
+        Assert(equ_ref_case != RefinementCase<dim>::no_refinement,
+               ExcInternalError());
+        Assert(equ_subface_no != e, ExcInternalError());
+        // now, finally respect non-standard faces
+        const RefinementCase<dim - 1> final_ref_case =
+          (face_orientation == face_rotation ?
+             ref_case_permutation[equ_ref_case] :
+             equ_ref_case);
+
+        const unsigned int final_subface_no =
+          GeometryInfo<dim>::child_cell_on_face(RefinementCase<dim>(
+                                                  final_ref_case),
+                                                4,
+                                                equ_subface_no,
+                                                face_orientation,
+                                                face_flip,
+                                                face_rotation,
+                                                equ_ref_case);
+
+        return std::make_pair(final_subface_no, final_ref_case);
+      }
     } // namespace
   }   // namespace QProjector
 } // namespace internal
@@ -250,6 +437,38 @@ QProjector<3>::project_to_face(const ReferenceCell &  reference_cell,
 }
 
 
+template <int dim>
+Quadrature<dim>
+QProjector<dim>::project_to_oriented_face(const ReferenceCell &reference_cell,
+                                          const Quadrature<dim - 1> &quadrature,
+                                          const unsigned int         face_no,
+                                          const bool,
+                                          const bool,
+                                          const bool)
+{
+  return QProjector<dim>::project_to_face(reference_cell, quadrature, face_no);
+}
+
+
+
+template <>
+Quadrature<3>
+QProjector<3>::project_to_oriented_face(const ReferenceCell &reference_cell,
+                                        const Quadrature<2> &quadrature,
+                                        const unsigned int   face_no,
+                                        const bool           face_orientation,
+                                        const bool           face_flip,
+                                        const bool           face_rotation)
+{
+  Assert(reference_cell == ReferenceCells::Hexahedron, ExcNotImplemented());
+
+  const Quadrature<2> mutation = internal::QProjector::mutate_quadrature(
+    quadrature, face_orientation, face_flip, face_rotation);
+
+  return QProjector<3>::project_to_face(reference_cell, mutation, face_no);
+}
+
+
 
 template <>
 void
@@ -440,6 +659,60 @@ QProjector<3>::project_to_subface(const ReferenceCell &    reference_cell,
 
 
 
+template <int dim>
+Quadrature<dim>
+QProjector<dim>::project_to_oriented_subface(
+  const ReferenceCell &      reference_cell,
+  const Quadrature<dim - 1> &quadrature,
+  const unsigned int         face_no,
+  const unsigned int         subface_no,
+  const bool,
+  const bool,
+  const bool,
+  const internal::SubfaceCase<dim>)
+{
+  return QProjector<dim>::project_to_subface(
+    reference_cell,
+    quadrature,
+    face_no,
+    subface_no,
+    RefinementCase<dim - 1>::isotropic_refinement);
+}
+
+
+
+template <>
+Quadrature<3>
+QProjector<3>::project_to_oriented_subface(
+  const ReferenceCell &          reference_cell,
+  const Quadrature<2> &          quadrature,
+  const unsigned int             face_no,
+  const unsigned int             subface_no,
+  const bool                     face_orientation,
+  const bool                     face_flip,
+  const bool                     face_rotation,
+  const internal::SubfaceCase<3> ref_case)
+{
+  Assert(reference_cell == ReferenceCells::Hexahedron, ExcNotImplemented());
+
+  const Quadrature<2> mutation = internal::QProjector::mutate_quadrature(
+    quadrature, face_orientation, face_flip, face_rotation);
+
+  const std::pair<unsigned int, RefinementCase<2>>
+    final_subface_no_and_ref_case =
+      internal::QProjector::select_subface_no_and_refinement_case(
+        subface_no, face_orientation, face_flip, face_rotation, ref_case);
+
+  return QProjector<3>::project_to_subface(
+    reference_cell,
+    mutation,
+    face_no,
+    final_subface_no_and_ref_case.first,
+    final_subface_no_and_ref_case.second);
+}
+
+
+
 template <>
 Quadrature<1>
 QProjector<1>::project_to_all_faces(const ReferenceCell &     reference_cell,
@@ -751,29 +1024,10 @@ QProjector<3>::project_to_all_faces(const ReferenceCell &     reference_cell,
       std::vector<double> weights;
       weights.reserve(n_points_total);
 
-      Table<2, std::vector<Point<2>>> mutations(quadrature.size(), 8);
-
-      // do the following for all possible mutations of a face (mutation==0
-      // corresponds to a face with standard orientation, no flip and no
-      // rotation)
-      for (unsigned int face = 0; face < quadrature.size(); ++face)
-        {
-          const std::vector<Point<2>> &quad = quadrature[face].get_points();
-          mutations(face, 0)                = quad;
-          mutations(face, 1) = internal::QProjector::rotate(quad, 1);
-          mutations(face, 2) = internal::QProjector::rotate(quad, 2);
-          mutations(face, 3) = internal::QProjector::rotate(quad, 3);
-          mutations(face, 4) = internal::QProjector::reflect(quad);
-          mutations(face, 5) =
-            internal::QProjector::rotate(mutations(face, 4), 3);
-          mutations(face, 6) =
-            internal::QProjector::rotate(mutations(face, 4), 2);
-          mutations(face, 7) =
-            internal::QProjector::rotate(mutations(face, 4), 1);
-        }
-
-      for (unsigned int i = 0; i < 8; ++i)
+      for (unsigned int offset = 0; offset < 8; ++offset)
         {
+          const auto mutation = internal::QProjector::mutate_points_with_offset(
+            quadrature[0].get_points(), offset);
           // project to each face and append results
           for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
                ++face)
@@ -781,7 +1035,11 @@ QProjector<3>::project_to_all_faces(const ReferenceCell &     reference_cell,
               const unsigned int q_index = quadrature.size() == 1 ? 0 : face;
 
               internal::QProjector::project_to_hex_face_and_append(
-                mutations(q_index, i), face, q_points);
+                q_index > 0 ? internal::QProjector::mutate_points_with_offset(
+                                quadrature[face].get_points(), offset) :
+                              mutation,
+                face,
+                q_points);
 
               std::copy(quadrature[q_index].get_weights().begin(),
                         quadrature[q_index].get_weights().end(),
@@ -789,7 +1047,6 @@ QProjector<3>::project_to_all_faces(const ReferenceCell &     reference_cell,
             }
         }
 
-
       Assert(q_points.size() == n_points_total, ExcInternalError());
       Assert(weights.size() == n_points_total, ExcInternalError());
 
@@ -908,18 +1165,7 @@ QProjector<3>::project_to_all_subfaces(const ReferenceCell &reference_cell,
 
   Assert(reference_cell == ReferenceCells::Hexahedron, ExcNotImplemented());
 
-  const unsigned int           dim  = 3;
-  const std::vector<Point<2>> &quad = quadrature.get_points();
-  std::vector<Point<2>> q_reflected = internal::QProjector::reflect(quad);
-  std::array<std::vector<Point<2>>, 8> mutations{
-    {quad,
-     internal::QProjector::rotate(quad, 1),
-     internal::QProjector::rotate(quad, 2),
-     internal::QProjector::rotate(quad, 3),
-     q_reflected,
-     internal::QProjector::rotate(q_reflected, 3),
-     internal::QProjector::rotate(q_reflected, 2),
-     internal::QProjector::rotate(q_reflected, 1)}};
+  const unsigned int dim = 3;
 
   const unsigned int n_points = quadrature.size(),
                      n_faces  = GeometryInfo<dim>::faces_per_cell,
@@ -934,8 +1180,12 @@ QProjector<3>::project_to_all_subfaces(const ReferenceCell &reference_cell,
 
   // do the following for all possible mutations of a face (mutation==0
   // corresponds to a face with standard orientation, no flip and no rotation)
-  for (const auto &mutation : mutations)
+  for (unsigned int offset = 0; offset < 8; ++offset)
     {
+      const auto mutation =
+        internal::QProjector::mutate_points_with_offset(quadrature.get_points(),
+                                                        offset);
+
       // project to each face and copy results
       for (unsigned int face = 0; face < n_faces; ++face)
         for (unsigned int ref_case = RefinementCase<dim - 1>::cut_xy;
@@ -1435,161 +1685,14 @@ QProjector<3>::DataSetDescriptor::subface(
     0  // cut_xy
   };
 
-
-  // for each subface of a given FaceRefineCase
-  // there is a corresponding equivalent
-  // subface number of one of the "standard"
-  // RefineCases (cut_x, cut_y, cut_xy). Map
-  // the given values to those equivalent
-  // ones.
-
-  // first, define an invalid number
-  static const unsigned int e = numbers::invalid_unsigned_int;
-
-  static const RefinementCase<dim - 1>
-    equivalent_refine_case[internal::SubfaceCase<dim>::case_isotropic + 1]
-                          [GeometryInfo<3>::max_children_per_face] = {
-                            // case_none. there should be only
-                            // invalid values here. However, as
-                            // this function is also called (in
-                            // tests) for cells which have no
-                            // refined faces, use isotropic
-                            // refinement instead
-                            {RefinementCase<dim - 1>::cut_xy,
-                             RefinementCase<dim - 1>::cut_xy,
-                             RefinementCase<dim - 1>::cut_xy,
-                             RefinementCase<dim - 1>::cut_xy},
-                            // case_x
-                            {RefinementCase<dim - 1>::cut_x,
-                             RefinementCase<dim - 1>::cut_x,
-                             RefinementCase<dim - 1>::no_refinement,
-                             RefinementCase<dim - 1>::no_refinement},
-                            // case_x1y
-                            {RefinementCase<dim - 1>::cut_xy,
-                             RefinementCase<dim - 1>::cut_xy,
-                             RefinementCase<dim - 1>::cut_x,
-                             RefinementCase<dim - 1>::no_refinement},
-                            // case_x2y
-                            {RefinementCase<dim - 1>::cut_x,
-                             RefinementCase<dim - 1>::cut_xy,
-                             RefinementCase<dim - 1>::cut_xy,
-                             RefinementCase<dim - 1>::no_refinement},
-                            // case_x1y2y
-                            {RefinementCase<dim - 1>::cut_xy,
-                             RefinementCase<dim - 1>::cut_xy,
-                             RefinementCase<dim - 1>::cut_xy,
-                             RefinementCase<dim - 1>::cut_xy},
-                            // case_y
-                            {RefinementCase<dim - 1>::cut_y,
-                             RefinementCase<dim - 1>::cut_y,
-                             RefinementCase<dim - 1>::no_refinement,
-                             RefinementCase<dim - 1>::no_refinement},
-                            // case_y1x
-                            {RefinementCase<dim - 1>::cut_xy,
-                             RefinementCase<dim - 1>::cut_xy,
-                             RefinementCase<dim - 1>::cut_y,
-                             RefinementCase<dim - 1>::no_refinement},
-                            // case_y2x
-                            {RefinementCase<dim - 1>::cut_y,
-                             RefinementCase<dim - 1>::cut_xy,
-                             RefinementCase<dim - 1>::cut_xy,
-                             RefinementCase<dim - 1>::no_refinement},
-                            // case_y1x2x
-                            {RefinementCase<dim - 1>::cut_xy,
-                             RefinementCase<dim - 1>::cut_xy,
-                             RefinementCase<dim - 1>::cut_xy,
-                             RefinementCase<dim - 1>::cut_xy},
-                            // case_xy (case_isotropic)
-                            {RefinementCase<dim - 1>::cut_xy,
-                             RefinementCase<dim - 1>::cut_xy,
-                             RefinementCase<dim - 1>::cut_xy,
-                             RefinementCase<dim - 1>::cut_xy}};
-
-  static const unsigned int
-    equivalent_subface_number[internal::SubfaceCase<dim>::case_isotropic + 1]
-                             [GeometryInfo<3>::max_children_per_face] = {
-                               // case_none, see above
-                               {0, 1, 2, 3},
-                               // case_x
-                               {0, 1, e, e},
-                               // case_x1y
-                               {0, 2, 1, e},
-                               // case_x2y
-                               {0, 1, 3, e},
-                               // case_x1y2y
-                               {0, 2, 1, 3},
-                               // case_y
-                               {0, 1, e, e},
-                               // case_y1x
-                               {0, 1, 1, e},
-                               // case_y2x
-                               {0, 2, 3, e},
-                               // case_y1x2x
-                               {0, 1, 2, 3},
-                               // case_xy (case_isotropic)
-                               {0, 1, 2, 3}};
-
-  // If face-orientation or face_rotation are
-  // non-standard, cut_x and cut_y have to be
-  // exchanged.
-  static const RefinementCase<dim - 1> ref_case_permutation[4] = {
-    RefinementCase<dim - 1>::no_refinement,
-    RefinementCase<dim - 1>::cut_y,
-    RefinementCase<dim - 1>::cut_x,
-    RefinementCase<dim - 1>::cut_xy};
-
-  // set a corresponding (equivalent)
-  // RefineCase and subface number
-  const RefinementCase<dim - 1> equ_ref_case =
-    equivalent_refine_case[ref_case][subface_no];
-  const unsigned int equ_subface_no =
-    equivalent_subface_number[ref_case][subface_no];
-  // make sure, that we got a valid subface and RefineCase
-  Assert(equ_ref_case != RefinementCase<dim>::no_refinement,
-         ExcInternalError());
-  Assert(equ_subface_no != e, ExcInternalError());
-  // now, finally respect non-standard faces
-  const RefinementCase<dim - 1> final_ref_case =
-    (face_orientation == face_rotation ? ref_case_permutation[equ_ref_case] :
-                                         equ_ref_case);
-
-  // what we have now is the number of
-  // the subface in the natural
-  // orientation of the *face*. what we
-  // need to know is the number of the
-  // subface concerning the standard face
-  // orientation as seen from the *cell*.
-
-  // this mapping is not trivial, but we
-  // have done exactly this stuff in the
-  // child_cell_on_face function. in
-  // order to reduce the amount of code
-  // as well as to make maintaining the
-  // functionality easier we want to
-  // reuse that information. So we note
-  // that on the bottom face (face 4) of
-  // a hex cell the local x and y
-  // coordinates of the face and the cell
-  // coincide, thus also the refinement
-  // case of the face corresponds to the
-  // refinement case of the cell
-  // (ignoring cell refinement along the
-  // z direction). Using this knowledge
-  // we can (ab)use the
-  // child_cell_on_face function to do
-  // exactly the transformation we are in
-  // need of now
-  const unsigned int final_subface_no =
-    GeometryInfo<dim>::child_cell_on_face(RefinementCase<dim>(final_ref_case),
-                                          4,
-                                          equ_subface_no,
-                                          face_orientation,
-                                          face_flip,
-                                          face_rotation,
-                                          equ_ref_case);
+  const std::pair<unsigned int, RefinementCase<2>>
+    final_subface_no_and_ref_case =
+      internal::QProjector::select_subface_no_and_refinement_case(
+        subface_no, face_orientation, face_flip, face_rotation, ref_case);
 
   return (((face_no * total_subfaces_per_face +
-            ref_case_offset[final_ref_case - 1] + final_subface_no) +
+            ref_case_offset[final_subface_no_and_ref_case.second - 1] +
+            final_subface_no_and_ref_case.first) +
            orientation_offset[face_orientation][face_flip][face_rotation]) *
           n_quadrature_points);
 }
index d945c3a0668cd9ddf2a9cd4867bd729fc72e031e..2961f7302611fecd64508fc78c0db04c7c722c17 100644 (file)
@@ -83,21 +83,14 @@ template <int dim, int spacedim>
 void
 MappingQ<dim, spacedim>::InternalData::initialize(
   const UpdateFlags      update_flags,
-  const Quadrature<dim> &q,
+  const Quadrature<dim> &quadrature,
   const unsigned int     n_original_q_points)
 {
   // store the flags in the internal data object so we can access them
   // in fill_fe_*_values()
   this->update_each = update_flags;
 
-  const unsigned int n_q_points = q.size();
-
-  const bool needs_higher_order_terms =
-    this->update_each &
-    (update_jacobian_pushed_forward_grads | update_jacobian_2nd_derivatives |
-     update_jacobian_pushed_forward_2nd_derivatives |
-     update_jacobian_3rd_derivatives |
-     update_jacobian_pushed_forward_3rd_derivatives);
+  const unsigned int n_q_points = quadrature.size();
 
   if (this->update_each & update_covariant_transformation)
     covariant.resize(n_original_q_points);
@@ -108,7 +101,7 @@ MappingQ<dim, spacedim>::InternalData::initialize(
   if (this->update_each & update_volume_elements)
     volume_elements.resize(n_original_q_points);
 
-  tensor_product_quadrature = q.is_tensor_product();
+  tensor_product_quadrature = quadrature.is_tensor_product();
 
   // use of MatrixFree only for higher order elements and with more than one
   // point where tensor products do not make sense
@@ -122,7 +115,7 @@ MappingQ<dim, spacedim>::InternalData::initialize(
       if (tensor_product_quadrature)
         {
           const std::array<Quadrature<1>, dim> &quad_array =
-            q.get_tensor_basis();
+            quadrature.get_tensor_basis();
           for (unsigned int i = 1; i < dim && tensor_product_quadrature; ++i)
             {
               if (quad_array[i - 1].size() != quad_array[i].size())
@@ -158,37 +151,43 @@ MappingQ<dim, spacedim>::InternalData::initialize(
               // numbering manually (building an FE_Q<dim> is relatively
               // expensive due to constraints)
               const FE_DGQ<1> fe(polynomial_degree);
-              shape_info.reinit(q.get_tensor_basis()[0], fe);
+              shape_info.reinit(quadrature.get_tensor_basis()[0], fe);
               shape_info.lexicographic_numbering =
                 FETools::lexicographic_to_hierarchic_numbering<dim>(
                   polynomial_degree);
-              shape_info.n_q_points = q.size();
+              shape_info.n_q_points = n_q_points;
               shape_info.dofs_per_component_on_cell =
                 Utilities::pow(polynomial_degree + 1, dim);
             }
         }
     }
 
+  const bool needs_higher_order_terms =
+    this->update_each &
+    (update_jacobian_pushed_forward_grads | update_jacobian_2nd_derivatives |
+     update_jacobian_pushed_forward_2nd_derivatives |
+     update_jacobian_3rd_derivatives |
+     update_jacobian_pushed_forward_3rd_derivatives);
+
+  const bool needs_higher_order_terms_generic =
+    !tensor_product_quadrature &&
+    (needs_higher_order_terms || this->update_each & update_jacobian_grads);
+
   // Only fill the big arrays on demand in case we cannot use the tensor
   // product quadrature code path
-  if (dim == 1 || !tensor_product_quadrature || needs_higher_order_terms)
+  if (dim == 1 || needs_higher_order_terms_generic || needs_higher_order_terms)
     {
-      // see if we need the (transformation) shape function values
-      // and/or gradients and resize the necessary arrays
-      if (this->update_each & update_quadrature_points)
-        shape_values.resize(n_shape_functions * n_q_points);
-
-      if (this->update_each &
-          (update_covariant_transformation |
-           update_contravariant_transformation | update_JxW_values |
-           update_boundary_forms | update_normal_vectors | update_jacobians |
-           update_jacobian_grads | update_inverse_jacobians |
-           update_jacobian_pushed_forward_grads |
-           update_jacobian_2nd_derivatives |
-           update_jacobian_pushed_forward_2nd_derivatives |
-           update_jacobian_3rd_derivatives |
-           update_jacobian_pushed_forward_3rd_derivatives))
-        shape_derivatives.resize(n_shape_functions * n_q_points);
+      // compute shapes and derivatives for codim1 (for
+      // do_transform_real_to_unit_cell_internal_codim1)
+      if (dim == (spacedim - 1))
+        {
+          // see if we need the (transformation) shape function values
+          // and/or gradients and resize the necessary arrays
+          if (this->update_each & update_quadrature_points)
+            shape_values.resize(n_shape_functions * n_q_points);
+          if (this->update_each & update_jacobians)
+            shape_derivatives.resize(n_shape_functions * n_q_points);
+        }
 
       if (this->update_each &
           (update_jacobian_grads | update_jacobian_pushed_forward_grads))
@@ -203,7 +202,7 @@ MappingQ<dim, spacedim>::InternalData::initialize(
         shape_fourth_derivatives.resize(n_shape_functions * n_q_points);
 
       // now also fill the various fields with their correct values
-      compute_shape_function_values(q.get_points());
+      compute_shape_function_values(quadrature.get_points());
     }
 }
 
@@ -213,16 +212,21 @@ template <int dim, int spacedim>
 void
 MappingQ<dim, spacedim>::InternalData::initialize_face(
   const UpdateFlags      update_flags,
-  const Quadrature<dim> &q,
+  const Quadrature<dim> &quadrature,
   const unsigned int     n_original_q_points)
 {
-  initialize(update_flags, q, n_original_q_points);
+  initialize(update_flags, quadrature, n_original_q_points);
+
+  const unsigned int n_q_points = quadrature.size();
+  quadrature_points.resize(n_q_points);
+  for (unsigned int q = 0; q < n_q_points; ++q)
+    quadrature_points[q] = quadrature.get_points()[q];
 
   if (dim > 1 && tensor_product_quadrature)
     {
       constexpr unsigned int facedim = dim - 1;
       const FE_DGQ<1>        fe(polynomial_degree);
-      shape_info.reinit(q.get_tensor_basis()[0], fe);
+      shape_info.reinit(quadrature.get_tensor_basis()[0], fe);
       shape_info.lexicographic_numbering =
         FETools::lexicographic_to_hierarchic_numbering<facedim>(
           polynomial_degree);
@@ -985,19 +989,22 @@ MappingQ<dim, spacedim>::fill_fe_values(
           computed_cell_similarity,
           data,
           output_data.quadrature_points,
+          output_data.jacobians,
+          output_data.inverse_jacobians,
           output_data.jacobian_grads);
     }
   else
     {
-      internal::MappingQImplementation::maybe_compute_q_points<dim, spacedim>(
-        QProjector<dim>::DataSetDescriptor::cell(),
-        data,
-        output_data.quadrature_points);
-
-      internal::MappingQImplementation::maybe_update_Jacobians<dim, spacedim>(
+      internal::MappingQImplementation::maybe_update_q_points_Jacobians_generic(
         computed_cell_similarity,
-        QProjector<dim>::DataSetDescriptor::cell(),
-        data);
+        data,
+        make_array_view(quadrature.get_points()),
+        polynomials_1d,
+        polynomial_degree,
+        renumber_lexicographic_to_hierarchic,
+        output_data.quadrature_points,
+        output_data.jacobians,
+        output_data.inverse_jacobians);
 
       internal::MappingQImplementation::maybe_update_jacobian_grads<dim,
                                                                     spacedim>(
@@ -1135,27 +1142,6 @@ MappingQ<dim, spacedim>::fill_fe_values(
           }
     }
 
-
-
-  // copy values from InternalData to vector given by reference
-  if (update_flags & update_jacobians)
-    {
-      AssertDimension(output_data.jacobians.size(), n_q_points);
-      if (computed_cell_similarity != CellSimilarity::translation)
-        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);
-      if (computed_cell_similarity != CellSimilarity::translation)
-        for (unsigned int point = 0; point < n_q_points; ++point)
-          output_data.inverse_jacobians[point] =
-            data.covariant[point].transpose();
-    }
-
   return computed_cell_similarity;
 }
 
@@ -1205,6 +1191,9 @@ MappingQ<dim, spacedim>::fill_fe_face_values(
       quadrature[0].size()),
     quadrature[0],
     data,
+    polynomials_1d,
+    polynomial_degree,
+    renumber_lexicographic_to_hierarchic,
     output_data);
 }
 
@@ -1255,6 +1244,9 @@ MappingQ<dim, spacedim>::fill_fe_subface_values(
       cell->subface_case(face_no)),
     quadrature,
     data,
+    polynomials_1d,
+    polynomial_degree,
+    renumber_lexicographic_to_hierarchic,
     output_data);
 }
 
@@ -1281,13 +1273,16 @@ MappingQ<dim, spacedim>::fill_fe_immersed_surface_values(
   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(),
+  internal::MappingQImplementation::maybe_update_q_points_Jacobians_generic(
+    CellSimilarity::none,
     data,
-    output_data.quadrature_points);
-
-  internal::MappingQImplementation::maybe_update_Jacobians<dim, spacedim>(
-    CellSimilarity::none, QProjector<dim>::DataSetDescriptor::cell(), data);
+    make_array_view(quadrature.get_points()),
+    polynomials_1d,
+    polynomial_degree,
+    renumber_lexicographic_to_hierarchic,
+    output_data.quadrature_points,
+    output_data.jacobians,
+    output_data.inverse_jacobians);
 
   internal::MappingQImplementation::maybe_update_jacobian_grads<dim, spacedim>(
     CellSimilarity::none,
@@ -1372,23 +1367,6 @@ MappingQ<dim, spacedim>::fill_fe_immersed_surface_values(
             }
         }
     }
-
-  // copy values from InternalData to vector given by reference
-  if ((update_flags & update_jacobians) != 0u)
-    {
-      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) != 0u)
-    {
-      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();
-    }
 }
 
 
@@ -1411,85 +1389,60 @@ MappingQ<dim, spacedim>::fill_mapping_data_for_generic_points(
          ExcNotImplemented());
 
   output_data.initialize(unit_points.size(), update_flags);
-  const std::vector<Point<spacedim>> support_points =
-    this->compute_mapping_support_points(cell);
 
-  const unsigned int n_points = unit_points.size();
-  const unsigned int n_lanes  = VectorizedArray<double>::size();
+  auto internal_data =
+    this->get_data(update_flags,
+                   Quadrature<dim>(std::vector<Point<dim>>(unit_points.begin(),
+                                                           unit_points.end())));
+  const InternalData &data = static_cast<const InternalData &>(*internal_data);
+  data.mapping_support_points = this->compute_mapping_support_points(cell);
 
-  // Use the more heavy VectorizedArray code path if there is more than
-  // one point left to compute
-  for (unsigned int i = 0; i < n_points; i += n_lanes)
-    if (n_points - i > 1)
-      {
-        Point<dim, VectorizedArray<double>> p_vec;
-        for (unsigned int j = 0; j < n_lanes; ++j)
-          if (i + j < n_points)
-            for (unsigned int d = 0; d < dim; ++d)
-              p_vec[d][j] = unit_points[i + j][d];
-          else
-            for (unsigned int d = 0; d < dim; ++d)
-              p_vec[d][j] = unit_points[i][d];
-
-        const auto result =
-          internal::evaluate_tensor_product_value_and_gradient(
-            polynomials_1d,
-            support_points,
-            p_vec,
-            polynomial_degree == 1,
-            renumber_lexicographic_to_hierarchic);
-
-        if (update_flags & update_quadrature_points)
-          for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
-            for (unsigned int d = 0; d < spacedim; ++d)
-              output_data.quadrature_points[i + j][d] = result.first[d][j];
+  internal::MappingQImplementation::maybe_update_q_points_Jacobians_generic(
+    CellSimilarity::none,
+    data,
+    unit_points,
+    polynomials_1d,
+    polynomial_degree,
+    renumber_lexicographic_to_hierarchic,
+    output_data.quadrature_points,
+    output_data.jacobians,
+    output_data.inverse_jacobians);
+}
 
-        if (update_flags & update_jacobians)
-          for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
-            for (unsigned int d = 0; d < spacedim; ++d)
-              for (unsigned int e = 0; e < dim; ++e)
-                output_data.jacobians[i + j][d][e] = result.second[e][d][j];
 
-        if (update_flags & update_inverse_jacobians)
-          {
-            DerivativeForm<1, spacedim, dim, VectorizedArray<double>> jac(
-              result.second);
-            const DerivativeForm<1, spacedim, dim, VectorizedArray<double>>
-              inv_jac = jac.covariant_form();
-            for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
-              for (unsigned int d = 0; d < dim; ++d)
-                for (unsigned int e = 0; e < spacedim; ++e)
-                  output_data.inverse_jacobians[i + j][d][e] = inv_jac[d][e][j];
-          }
-      }
-    else
-      {
-        const auto result =
-          internal::evaluate_tensor_product_value_and_gradient(
-            polynomials_1d,
-            support_points,
-            unit_points[i],
-            polynomial_degree == 1,
-            renumber_lexicographic_to_hierarchic);
-
-        if (update_flags & update_quadrature_points)
-          output_data.quadrature_points[i] = result.first;
-
-        if (update_flags & update_jacobians)
-          {
-            DerivativeForm<1, spacedim, dim> jac = result.second;
-            output_data.jacobians[i]             = jac.transpose();
-          }
 
-        if (update_flags & update_inverse_jacobians)
-          {
-            DerivativeForm<1, spacedim, dim> jac(result.second);
-            DerivativeForm<1, spacedim, dim> inv_jac = jac.covariant_form();
-            for (unsigned int d = 0; d < dim; ++d)
-              for (unsigned int e = 0; e < spacedim; ++e)
-                output_data.inverse_jacobians[i][d][e] = inv_jac[d][e];
-          }
-      }
+template <int dim, int spacedim>
+void
+MappingQ<dim, spacedim>::fill_mapping_data_for_face_quadrature(
+  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+  const unsigned int                                          face_no,
+  const Quadrature<dim - 1> &                                 face_quadrature,
+  const typename Mapping<dim, spacedim>::InternalDataBase &   internal_data,
+  internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
+    &output_data) const
+{
+  if (face_quadrature.get_points().empty())
+    return;
+
+  // 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);
+
+  data.mapping_support_points = this->compute_mapping_support_points(cell);
+
+  internal::MappingQImplementation::do_fill_fe_face_values(
+    *this,
+    cell,
+    face_no,
+    numbers::invalid_unsigned_int,
+    QProjector<dim>::DataSetDescriptor::cell(),
+    face_quadrature,
+    data,
+    polynomials_1d,
+    polynomial_degree,
+    renumber_lexicographic_to_hierarchic,
+    output_data);
 }
 
 

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.