]> https://gitweb.dealii.org/ - dealii.git/commitdiff
FEEval::reinit(std::array): use MappingDataOnTheFly 13327/head
authorPeter Munch <peterrmuench@gmail.com>
Thu, 3 Feb 2022 09:42:03 +0000 (10:42 +0100)
committerPeter Munch <peterrmuench@gmail.com>
Thu, 3 Feb 2022 09:42:03 +0000 (10:42 +0100)
include/deal.II/matrix_free/fe_evaluation.h
include/deal.II/matrix_free/mapping_data_on_the_fly.h

index 449743975f5bc06c49d09527a24636c5818bcc95..0803f717a2f9dbc07ea4de9fbc0c7ff8becad7cb 100644 (file)
@@ -740,33 +740,6 @@ protected:
    * MatrixFree object was given at initialization.
    */
   mutable std::vector<types::global_dof_index> local_dof_indices;
-
-  /**
-   * A temporary data structure for the Jacobian information necessary if
-   * the reinit() function is used that allows to construct user batches.
-   */
-  AlignedVector<Tensor<2, dim, VectorizedArrayType>> jacobian_data;
-
-  /**
-   * A temporary data structure for the Jacobian determinant necessary if
-   * the reinit() function is used that allows to construct user batches.
-   */
-  AlignedVector<VectorizedArrayType> J_value_data;
-
-  /**
-   * A temporary data structure for the gradients of the inverse Jacobian
-   * transformation necessary if the reinit() function is used that allows
-   * to construct user batches.
-   */
-  AlignedVector<
-    Tensor<1, dim *(dim + 1) / 2, Tensor<1, dim, VectorizedArrayType>>>
-    jacobian_gradients_data;
-
-  /**
-   * A temporary data structure for the quadrature-point locations necessary if
-   * the reinit() function is used that allows to construct user batches.
-   */
-  AlignedVector<Point<dim, VectorizedArrayType>> quadrature_points_data;
 };
 
 
@@ -6818,12 +6791,6 @@ FEEvaluation<dim,
              VectorizedArrayType>::
   reinit(const std::array<unsigned int, VectorizedArrayType::size()> &cell_ids)
 {
-  Assert(this->mapped_geometry == nullptr,
-         ExcMessage("FEEvaluation was initialized without a matrix-free object."
-                    " Integer indexing is not possible"));
-  if (this->mapped_geometry != nullptr)
-    return;
-
   Assert(this->dof_info != nullptr, ExcNotInitialized());
   Assert(this->mapping_data != nullptr, ExcNotInitialized());
 
@@ -6847,41 +6814,52 @@ FEEvaluation<dim,
     }
 
   // allocate memory for internal data storage
+  if (this->mapped_geometry == nullptr)
+    this->mapped_geometry =
+      std::make_shared<internal::MatrixFreeFunctions::
+                         MappingDataOnTheFly<dim, VectorizedArrayType>>();
+
+  auto &mapping_storage = this->mapped_geometry->get_data_storage();
+
+  auto &this_jacobian_data           = mapping_storage.jacobians[0];
+  auto &this_J_value_data            = mapping_storage.JxW_values;
+  auto &this_jacobian_gradients_data = mapping_storage.jacobian_gradients[0];
+  auto &this_quadrature_points_data  = mapping_storage.quadrature_points;
 
   if (this->cell_type <= internal::MatrixFreeFunctions::GeometryType::affine)
     {
       if (this->mapping_data->jacobians[0].size() > 0)
-        this->jacobian_data.resize_fast(2);
+        this_jacobian_data.resize_fast(2);
 
       if (this->mapping_data->JxW_values.size() > 0)
-        this->J_value_data.resize_fast(1);
+        this_J_value_data.resize_fast(1);
 
       if (this->mapping_data->jacobian_gradients[0].size() > 0)
-        this->jacobian_gradients_data.resize_fast(1);
+        this_jacobian_gradients_data.resize_fast(1);
 
       if (this->mapping_data->quadrature_points.size() > 0)
-        this->quadrature_points_data.resize_fast(1);
+        this_quadrature_points_data.resize_fast(1);
     }
   else
     {
       if (this->mapping_data->jacobians[0].size() > 0)
-        this->jacobian_data.resize_fast(this->n_quadrature_points);
+        this_jacobian_data.resize_fast(this->n_quadrature_points);
 
       if (this->mapping_data->JxW_values.size() > 0)
-        this->J_value_data.resize_fast(this->n_quadrature_points);
+        this_J_value_data.resize_fast(this->n_quadrature_points);
 
       if (this->mapping_data->jacobian_gradients[0].size() > 0)
-        this->jacobian_gradients_data.resize_fast(this->n_quadrature_points);
+        this_jacobian_gradients_data.resize_fast(this->n_quadrature_points);
 
       if (this->mapping_data->quadrature_points.size() > 0)
-        this->quadrature_points_data.resize_fast(this->n_quadrature_points);
+        this_quadrature_points_data.resize_fast(this->n_quadrature_points);
     }
 
   // set pointers to internal data storage
-  this->jacobian           = this->jacobian_data.data();
-  this->J_value            = this->J_value_data.data();
-  this->jacobian_gradients = this->jacobian_gradients_data.data();
-  this->quadrature_points  = this->quadrature_points_data.data();
+  this->jacobian           = this_jacobian_data.data();
+  this->J_value            = this_J_value_data.data();
+  this->jacobian_gradients = this_jacobian_gradients_data.data();
+  this->quadrature_points  = this_quadrature_points_data.data();
 
   // fill internal data storage lane by lane
   for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
@@ -6904,26 +6882,26 @@ FEEvaluation<dim,
           const unsigned int q = 0;
 
           if (this->mapping_data->JxW_values.size() > 0)
-            this->J_value_data[q][v] =
+            this_J_value_data[q][v] =
               this->mapping_data->JxW_values[offsets + q][lane];
 
           if (this->mapping_data->jacobians[0].size() > 0)
             for (unsigned int q = 0; q < 2; ++q)
               for (unsigned int i = 0; i < dim; ++i)
                 for (unsigned int j = 0; j < dim; ++j)
-                  this->jacobian_data[q][i][j][v] =
+                  this_jacobian_data[q][i][j][v] =
                     this->mapping_data->jacobians[0][offsets + q][i][j][lane];
 
           if (this->mapping_data->jacobian_gradients[0].size() > 0)
             for (unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
               for (unsigned int j = 0; j < dim; ++j)
-                this->jacobian_gradients_data[q][i][j][v] =
+                this_jacobian_gradients_data[q][i][j][v] =
                   this->mapping_data
                     ->jacobian_gradients[0][offsets + q][i][j][lane];
 
           if (this->mapping_data->quadrature_points.size() > 0)
             for (unsigned int i = 0; i < dim; ++i)
-              this->quadrature_points_data[q][i][v] =
+              this_quadrature_points_data[q][i][v] =
                 this->mapping_data->quadrature_points
                   [this->mapping_data
                      ->quadrature_point_offsets[cell_batch_index] +
@@ -6945,20 +6923,20 @@ FEEvaluation<dim,
                   q;
 
               if (this->mapping_data->JxW_values.size() > 0)
-                this->J_value_data[q][v] =
+                this_J_value_data[q][v] =
                   this->mapping_data->JxW_values[offsets + q_src][lane];
 
               if (this->mapping_data->jacobians[0].size() > 0)
                 for (unsigned int i = 0; i < dim; ++i)
                   for (unsigned int j = 0; j < dim; ++j)
-                    this->jacobian_data[q][i][j][v] =
+                    this_jacobian_data[q][i][j][v] =
                       this->mapping_data
                         ->jacobians[0][offsets + q_src][i][j][lane];
 
               if (this->mapping_data->jacobian_gradients[0].size() > 0)
                 for (unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
                   for (unsigned int j = 0; j < dim; ++j)
-                    this->jacobian_gradients_data[q][i][j][v] =
+                    this_jacobian_gradients_data[q][i][j][v] =
                       this->mapping_data
                         ->jacobian_gradients[0][offsets + q_src][i][j][lane];
 
@@ -6993,13 +6971,13 @@ FEEvaluation<dim,
                                 this->descriptor->quadrature.point(q)[e]);
 
                       for (unsigned int i = 0; i < dim; ++i)
-                        this->quadrature_points_data[q][i][v] = point[i][lane];
+                        this_quadrature_points_data[q][i][v] = point[i][lane];
                     }
                   else
                     {
                       // general case: quadrature points are available
                       for (unsigned int i = 0; i < dim; ++i)
-                        this->quadrature_points_data[q][i][v] =
+                        this_quadrature_points_data[q][i][v] =
                           this->mapping_data->quadrature_points
                             [this->mapping_data
                                ->quadrature_point_offsets[cell_batch_index] +
index ebefb06f1f70edc4ffcf7d05921e7526e1ed99e7..a5d376331049e57b97644b3ccad8ab382934af47 100644 (file)
@@ -60,6 +60,11 @@ namespace internal
     class MappingDataOnTheFly
     {
     public:
+      /**
+       * Default constructor.
+       */
+      MappingDataOnTheFly() = default;
+
       /**
        * Constructor, similar to FEValues. Since this class only evaluates the
        * geometry, no finite element has to be specified and the simplest
@@ -115,6 +120,16 @@ namespace internal
       const MappingInfoStorage<dim, dim, Number> &
       get_data_storage() const;
 
+      /**
+       * Return a non-const reference to the underlying storage field
+       * of type MappingInfoStorage of the same format as the data fields
+       * in MappingInfo. This function can be used to manually fill the
+       * data if the default constructor has been called and reinit() was
+       * not called for a given cell.
+       */
+      MappingInfoStorage<dim, dim, Number> &
+      get_data_storage();
+
       /**
        * Return a reference to 1D quadrature underlying this object.
        */
@@ -137,7 +152,7 @@ namespace internal
       /**
        * An underlying FEValues object that performs the (scalar) evaluation.
        */
-      dealii::FEValues<dim> fe_values;
+      std::unique_ptr<dealii::FEValues<dim>> fe_values;
 
       /**
        * Get 1D quadrature formula to be used for reinitializing shape info.
@@ -159,32 +174,33 @@ namespace internal
       const Mapping<dim> & mapping,
       const Quadrature<1> &quadrature,
       const UpdateFlags    update_flags)
-      : fe_values(mapping,
-                  fe_dummy,
-                  Quadrature<dim>(quadrature),
-                  MappingInfoStorage<dim, dim, Number>::compute_update_flags(
-                    update_flags))
+      : fe_values(std::make_unique<dealii::FEValues<dim>>(
+          mapping,
+          fe_dummy,
+          Quadrature<dim>(quadrature),
+          MappingInfoStorage<dim, dim, Number>::compute_update_flags(
+            update_flags)))
       , quadrature_1d(quadrature)
     {
       mapping_info_storage.descriptor.resize(1);
       mapping_info_storage.descriptor[0].initialize(quadrature);
       mapping_info_storage.data_index_offsets.resize(1);
-      mapping_info_storage.JxW_values.resize(fe_values.n_quadrature_points);
-      mapping_info_storage.jacobians[0].resize(fe_values.n_quadrature_points);
+      mapping_info_storage.JxW_values.resize(fe_values->n_quadrature_points);
+      mapping_info_storage.jacobians[0].resize(fe_values->n_quadrature_points);
       if ((update_flags & update_quadrature_points) != 0u)
         {
           mapping_info_storage.quadrature_point_offsets.resize(1, 0);
           mapping_info_storage.quadrature_points.resize(
-            fe_values.n_quadrature_points);
+            fe_values->n_quadrature_points);
         }
-      if (fe_values.get_update_flags() & update_normal_vectors)
+      if (fe_values->get_update_flags() & update_normal_vectors)
         {
           mapping_info_storage.normal_vectors.resize(
-            fe_values.n_quadrature_points);
+            fe_values->n_quadrature_points);
           mapping_info_storage.normals_times_jacobians[0].resize(
-            fe_values.n_quadrature_points);
+            fe_values->n_quadrature_points);
         }
-      Assert(!(fe_values.get_update_flags() & update_jacobian_grads),
+      Assert(!(fe_values->get_update_flags() & update_jacobian_grads),
              ExcNotImplemented());
     }
 
@@ -209,28 +225,28 @@ namespace internal
       if (present_cell == cell)
         return;
       present_cell = cell;
-      fe_values.reinit(present_cell);
-      for (unsigned int q = 0; q < fe_values.get_quadrature().size(); ++q)
+      fe_values->reinit(present_cell);
+      for (unsigned int q = 0; q < fe_values->get_quadrature().size(); ++q)
         {
-          if (fe_values.get_update_flags() & update_JxW_values)
-            mapping_info_storage.JxW_values[q] = fe_values.JxW(q);
-          if (fe_values.get_update_flags() & update_jacobians)
+          if (fe_values->get_update_flags() & update_JxW_values)
+            mapping_info_storage.JxW_values[q] = fe_values->JxW(q);
+          if (fe_values->get_update_flags() & update_jacobians)
             {
-              Tensor<2, dim> jac = fe_values.jacobian(q);
+              Tensor<2, dim> jac = fe_values->jacobian(q);
               jac                = invert(transpose(jac));
               for (unsigned int d = 0; d < dim; ++d)
                 for (unsigned int e = 0; e < dim; ++e)
                   mapping_info_storage.jacobians[0][q][d][e] = jac[d][e];
             }
-          if (fe_values.get_update_flags() & update_quadrature_points)
+          if (fe_values->get_update_flags() & update_quadrature_points)
             for (unsigned int d = 0; d < dim; ++d)
               mapping_info_storage.quadrature_points[q][d] =
-                fe_values.quadrature_point(q)[d];
-          if (fe_values.get_update_flags() & update_normal_vectors)
+                fe_values->quadrature_point(q)[d];
+          if (fe_values->get_update_flags() & update_normal_vectors)
             {
               for (unsigned int d = 0; d < dim; ++d)
                 mapping_info_storage.normal_vectors[q][d] =
-                  fe_values.normal_vector(q)[d];
+                  fe_values->normal_vector(q)[d];
               mapping_info_storage.normals_times_jacobians[0][q] =
                 mapping_info_storage.normal_vectors[q] *
                 mapping_info_storage.jacobians[0][q];
@@ -254,7 +270,7 @@ namespace internal
     inline typename dealii::Triangulation<dim>::cell_iterator
     MappingDataOnTheFly<dim, Number>::get_cell() const
     {
-      return fe_values.get_cell();
+      return fe_values->get_cell();
     }
 
 
@@ -263,7 +279,7 @@ namespace internal
     inline const dealii::FEValues<dim> &
     MappingDataOnTheFly<dim, Number>::get_fe_values() const
     {
-      return fe_values;
+      return *fe_values;
     }
 
 
@@ -277,6 +293,15 @@ namespace internal
 
 
 
+    template <int dim, typename Number>
+    inline MappingInfoStorage<dim, dim, Number> &
+    MappingDataOnTheFly<dim, Number>::get_data_storage()
+    {
+      return mapping_info_storage;
+    }
+
+
+
     template <int dim, typename Number>
     inline const Quadrature<1> &
     MappingDataOnTheFly<dim, Number>::get_quadrature() const

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.