Move quadrature_points into FEEvalData

diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h
index a5c4e384ba006f6d9bb6669060aea10ba3a5c26a..7e4bcae909b6ade95083de39c48c783a71572de6 100644 (file)
--- a/include/deal.II/matrix_free/fe_evaluation.h
+++ b/include/deal.II/matrix_free/fe_evaluation.h
@@ -2331,13 +2331,6 @@ public:
                     const bool  integrate_gradients,
                     VectorType &output_vector);
 
-  /**
-   * Return the q-th quadrature point in real coordinates stored in
-   * MappingInfo.
-   */
-  Point<dim, VectorizedArrayType>
-  quadrature_point(const unsigned int q_point) const;
-
   /**
    * The number of degrees of freedom of a single component on the cell for
    * the underlying evaluation object. Usually close to
@@ -2736,13 +2729,6 @@ public:
                     const bool  integrate_gradients,
                     VectorType &output_vector);
 
-  /**
-   * Returns the q-th quadrature point on the face in real coordinates stored
-   * in MappingInfo.
-   */
-  Point<dim, VectorizedArrayType>
-  quadrature_point(const unsigned int q_point) const;
-
   /**
    * The number of degrees of freedom of a single component on the cell for
    * the underlying evaluation object. Usually close to
@@ -3032,6 +3018,10 @@ inline FEEvaluationBase<dim,
         this->mapped_geometry->get_data_storage().jacobians[0].begin();
       this->J_value =
         this->mapped_geometry->get_data_storage().JxW_values.begin();
+      this->jacobian_gradients =
+        this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
+      this->quadrature_points =
+        this->mapped_geometry->get_data_storage().quadrature_points.begin();
     }
 
   this->set_data_pointers(scratch_data_array, n_components_);
@@ -3092,6 +3082,10 @@ operator=(const FEEvaluationBase<dim,
         this->mapped_geometry->get_data_storage().jacobians[0].begin();
       this->J_value =
         this->mapped_geometry->get_data_storage().JxW_values.begin();
+      this->jacobian_gradients =
+        this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
+      this->quadrature_points =
+        this->mapped_geometry->get_data_storage().quadrature_points.begin();
     }
   else
     {
@@ -6870,6 +6864,11 @@ FEEvaluation<dim,
   for (unsigned int i = 0; i < VectorizedArrayType::size(); ++i)
     this->cell_ids[i] = cell_index * VectorizedArrayType::size() + i;
 
+  if (this->mapping_data->quadrature_points.empty() == false)
+    this->quadrature_points =
+      &this->mapping_data->quadrature_points
+         [this->mapping_data->quadrature_point_offsets[this->cell]];
+
 #  ifdef DEBUG
   this->is_reinitialized           = true;
   this->dof_values_initialized     = false;
@@ -6948,66 +6947,6 @@ FEEvaluation<dim,
 
 
 
-template <int dim,
-          int fe_degree,
-          int n_q_points_1d,
-          int n_components_,
-          typename Number,
-          typename VectorizedArrayType>
-inline Point<dim, VectorizedArrayType>
-FEEvaluation<dim,
-             fe_degree,
-             n_q_points_1d,
-             n_components_,
-             Number,
-             VectorizedArrayType>::quadrature_point(const unsigned int q) const
-{
-  if (this->matrix_free == nullptr)
-    {
-      Assert((this->mapped_geometry->get_fe_values().get_update_flags() |
-              update_quadrature_points),
-             internal::ExcMatrixFreeAccessToUninitializedMappingField(
-               "update_quadrature_points"));
-    }
-  else
-    {
-      Assert(this->mapping_data->quadrature_point_offsets.empty() == false,
-             internal::ExcMatrixFreeAccessToUninitializedMappingField(
-               "update_quadrature_points"));
-    }
-
-  AssertIndexRange(q, n_q_points);
-
-  const Point<dim, VectorizedArrayType> *quadrature_points =
-    &this->mapping_data->quadrature_points
-       [this->mapping_data->quadrature_point_offsets[this->cell]];
-
-  // Cartesian/affine mesh: only first vertex of cell is stored, we must
-  // compute it through the Jacobian (which is stored in non-inverted and
-  // non-transposed form as index '1' in the jacobian field)
-  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
-    {
-      Assert(this->jacobian != nullptr, ExcNotInitialized());
-      Point<dim, VectorizedArrayType> point = quadrature_points[0];
-
-      const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
-      if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
-        for (unsigned int d = 0; d < dim; ++d)
-          point[d] += jac[d][d] * static_cast<Number>(
-                                    this->descriptor->quadrature.point(q)[d]);
-      else
-        for (unsigned int d = 0; d < dim; ++d)
-          for (unsigned int e = 0; e < dim; ++e)
-            point[d] += jac[d][e] * static_cast<Number>(
-                                      this->descriptor->quadrature.point(q)[e]);
-      return point;
-    }
-  else
-    return quadrature_points[q];
-}
-
-
-
 template <int dim,
           int fe_degree,
           int n_q_points_1d,
@@ -7615,6 +7554,15 @@ FEFaceEvaluation<dim,
     this->mapping_data->jacobian_gradients[!this->is_interior_face].data() +
     offsets;
 
+  if (this->mapping_data->quadrature_point_offsets.empty() == false)
+    {
+      AssertIndexRange(this->cell,
+                       this->mapping_data->quadrature_point_offsets.size());
+      this->quadrature_points =
+        this->mapping_data->quadrature_points.data() +
+        this->mapping_data->quadrature_point_offsets[this->cell];
+    }
+
 #  ifdef DEBUG
   this->is_reinitialized           = true;
   this->dof_values_initialized     = false;
@@ -7755,6 +7703,24 @@ FEFaceEvaluation<dim,
     this->mapping_data->jacobian_gradients[!this->is_interior_face].data() +
     offsets;
 
+  if (this->matrix_free->get_mapping_info()
+        .face_data_by_cells[this->quad_no]
+        .quadrature_point_offsets.empty() == false)
+    {
+      const unsigned int index =
+        this->cell * GeometryInfo<dim>::faces_per_cell + this->face_numbers[0];
+      AssertIndexRange(index,
+                       this->matrix_free->get_mapping_info()
+                         .face_data_by_cells[this->quad_no]
+                         .quadrature_point_offsets.size());
+      this->quadrature_points = this->matrix_free->get_mapping_info()
+                                  .face_data_by_cells[this->quad_no]
+                                  .quadrature_points.data() +
+                                this->matrix_free->get_mapping_info()
+                                  .face_data_by_cells[this->quad_no]
+                                  .quadrature_point_offsets[index];
+    }
+
 #  ifdef DEBUG
   this->is_reinitialized           = true;
   this->dof_values_initialized     = false;
@@ -8246,56 +8212,6 @@ FEFaceEvaluation<dim,
 
 
 
-template <int dim,
-          int fe_degree,
-          int n_q_points_1d,
-          int n_components_,
-          typename Number,
-          typename VectorizedArrayType>
-inline Point<dim, VectorizedArrayType>
-FEFaceEvaluation<dim,
-                 fe_degree,
-                 n_q_points_1d,
-                 n_components_,
-                 Number,
-                 VectorizedArrayType>::quadrature_point(const unsigned int q)
-  const
-{
-  AssertIndexRange(q, n_q_points);
-  if (this->dof_access_index < 2)
-    {
-      Assert(this->mapping_data->quadrature_point_offsets.empty() == false,
-             internal::ExcMatrixFreeAccessToUninitializedMappingField(
-               "update_quadrature_points"));
-      AssertIndexRange(this->cell,
-                       this->mapping_data->quadrature_point_offsets.size());
-      return this->mapping_data->quadrature_points
-        [this->mapping_data->quadrature_point_offsets[this->cell] + q];
-    }
-  else
-    {
-      Assert(this->matrix_free->get_mapping_info()
-                 .face_data_by_cells[this->quad_no]
-                 .quadrature_point_offsets.empty() == false,
-             internal::ExcMatrixFreeAccessToUninitializedMappingField(
-               "update_quadrature_points"));
-      const unsigned int index =
-        this->cell * GeometryInfo<dim>::faces_per_cell + this->face_numbers[0];
-      AssertIndexRange(index,
-                       this->matrix_free->get_mapping_info()
-                         .face_data_by_cells[this->quad_no]
-                         .quadrature_point_offsets.size());
-      return this->matrix_free->get_mapping_info()
-        .face_data_by_cells[this->quad_no]
-        .quadrature_points[this->matrix_free->get_mapping_info()
-                             .face_data_by_cells[this->quad_no]
-                             .quadrature_point_offsets[index] +
-                           q];
-    }
-}
-
-
-
 template <int dim,
           int fe_degree,
           int n_q_points_1d,

diff --git a/include/deal.II/matrix_free/fe_evaluation_data.h b/include/deal.II/matrix_free/fe_evaluation_data.h
index 3c0cfd2a7fe79c1a6570076dd48c60688d464c86..b645f2edbb168281af8fa543051260f4bbdafd74 100644 (file)
--- a/include/deal.II/matrix_free/fe_evaluation_data.h
+++ b/include/deal.II/matrix_free/fe_evaluation_data.h
@@ -185,6 +185,13 @@ public:
   Number
   JxW(const unsigned int q_point) const;
 
+  /**
+   * Returns the q-th quadrature point on the face in real coordinates stored
+   * in MappingInfo.
+   */
+  Point<dim, Number>
+  quadrature_point(const unsigned int q) const;
+
   /**
    * Return the inverse and transposed version $J^{-\mathrm T}$ of the
    * Jacobian of the mapping between the unit to the real cell defined as
@@ -608,6 +615,12 @@ protected:
    */
   const unsigned int n_quadrature_points;
 
+  /**
+   * A pointer to the quadrature-point information of the present cell.
+   * Only set to a useful value if on a non-Cartesian cell.
+   */
+  const Point<dim, Number> *quadrature_points;
+
   /**
    * A pointer to the Jacobian information of the present cell. Only set to a
    * useful value if on a non-Cartesian cell.
@@ -901,7 +914,9 @@ inline FEEvaluationData<dim, Number, is_face>::FEEvaluationData(
   , n_quadrature_points(descriptor == nullptr ?
                           (is_face ? data->n_q_points_face : data->n_q_points) :
                           descriptor->n_q_points)
+  , quadrature_points(nullptr)
   , jacobian(nullptr)
+  , jacobian_gradients(nullptr)
   , J_value(nullptr)
   , normal_vectors(nullptr)
   , normal_x_jacobian(nullptr)
@@ -948,7 +963,9 @@ inline FEEvaluationData<dim, Number, is_face>::FEEvaluationData(
   , descriptor(nullptr)
   , n_quadrature_points(
       mapped_geometry->get_data_storage().descriptor[0].n_q_points)
+  , quadrature_points(nullptr)
   , jacobian(nullptr)
+  , jacobian_gradients(nullptr)
   , J_value(nullptr)
   , normal_vectors(nullptr)
   , normal_x_jacobian(nullptr)
@@ -963,6 +980,10 @@ inline FEEvaluationData<dim, Number, is_face>::FEEvaluationData(
   mapping_data = &mapped_geometry->get_data_storage();
   jacobian     = mapped_geometry->get_data_storage().jacobians[0].begin();
   J_value      = mapped_geometry->get_data_storage().JxW_values.begin();
+  jacobian_gradients =
+    mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
+  quadrature_points =
+    mapped_geometry->get_data_storage().quadrature_points.begin();
 }
 
 
@@ -986,6 +1007,8 @@ FEEvaluationData<dim, Number, is_face>::operator=(const FEEvaluationData &other)
   J_value            = nullptr;
   normal_vectors     = nullptr;
   normal_x_jacobian  = nullptr;
+  jacobian_gradients = nullptr;
+  quadrature_points  = nullptr;
   quadrature_weights = other.quadrature_weights;
 
 #  ifdef DEBUG
@@ -1124,6 +1147,43 @@ FEEvaluationData<dim, Number, is_face>::JxW(const unsigned int q_point) const
 
 
 
+template <int dim, typename Number, bool is_face>
+inline DEAL_II_ALWAYS_INLINE Point<dim, Number>
+FEEvaluationData<dim, Number, is_face>::quadrature_point(
+  const unsigned int q) const
+{
+  AssertIndexRange(q, this->n_quadrature_points);
+  Assert(this->quadrature_points != nullptr,
+         internal::ExcMatrixFreeAccessToUninitializedMappingField(
+           "update_quadrature_points"));
+
+  // Cartesian/affine mesh: only first vertex of cell is stored, we must
+  // compute it through the Jacobian (which is stored in non-inverted and
+  // non-transposed form as index '1' in the jacobian field)
+  if (is_face == false &&
+      this->cell_type <= internal::MatrixFreeFunctions::affine)
+    {
+      Assert(this->jacobian != nullptr, ExcNotInitialized());
+      Point<dim, Number> point = this->quadrature_points[0];
+
+      const Tensor<2, dim, Number> &jac = this->jacobian[1];
+      if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
+        for (unsigned int d = 0; d < dim; ++d)
+          point[d] += jac[d][d] * static_cast<typename Number::value_type>(
+                                    this->descriptor->quadrature.point(q)[d]);
+      else
+        for (unsigned int d = 0; d < dim; ++d)
+          for (unsigned int e = 0; e < dim; ++e)
+            point[d] += jac[d][e] * static_cast<typename Number::value_type>(
+                                      this->descriptor->quadrature.point(q)[e]);
+      return point;
+    }
+  else
+    return this->quadrature_points[q];
+}
+
+
+
 template <int dim, typename Number, bool is_face>
 inline DEAL_II_ALWAYS_INLINE Tensor<2, dim, Number>
 FEEvaluationData<dim, Number, is_face>::inverse_jacobian(