From: Peter Munch Date: Thu, 3 Feb 2022 09:42:03 +0000 (+0100) Subject: FEEval::reinit(std::array): use MappingDataOnTheFly X-Git-Tag: v9.4.0-rc1~527^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=1b51f4a1d764b4f48397c4d4e934b64eafbe6ae1;p=dealii.git FEEval::reinit(std::array): use MappingDataOnTheFly --- diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index 449743975f..0803f717a2 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -740,33 +740,6 @@ protected: * MatrixFree object was given at initialization. */ mutable std::vector 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> jacobian_data; - - /** - * A temporary data structure for the Jacobian determinant necessary if - * the reinit() function is used that allows to construct user batches. - */ - AlignedVector 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> quadrature_points_data; }; @@ -6818,12 +6791,6 @@ FEEvaluation:: reinit(const std::array &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 @@ FEEvaluationmapped_geometry == nullptr) + this->mapped_geometry = + std::make_shared>(); + + 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 @@ FEEvaluationmapping_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 @@ FEEvaluationmapping_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 @@ FEEvaluationdescriptor->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] + diff --git a/include/deal.II/matrix_free/mapping_data_on_the_fly.h b/include/deal.II/matrix_free/mapping_data_on_the_fly.h index ebefb06f1f..a5d3763310 100644 --- a/include/deal.II/matrix_free/mapping_data_on_the_fly.h +++ b/include/deal.II/matrix_free/mapping_data_on_the_fly.h @@ -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 & 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 & + 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 fe_values; + std::unique_ptr> fe_values; /** * Get 1D quadrature formula to be used for reinitializing shape info. @@ -159,32 +174,33 @@ namespace internal const Mapping & mapping, const Quadrature<1> &quadrature, const UpdateFlags update_flags) - : fe_values(mapping, - fe_dummy, - Quadrature(quadrature), - MappingInfoStorage::compute_update_flags( - update_flags)) + : fe_values(std::make_unique>( + mapping, + fe_dummy, + Quadrature(quadrature), + MappingInfoStorage::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::cell_iterator MappingDataOnTheFly::get_cell() const { - return fe_values.get_cell(); + return fe_values->get_cell(); } @@ -263,7 +279,7 @@ namespace internal inline const dealii::FEValues & MappingDataOnTheFly::get_fe_values() const { - return fe_values; + return *fe_values; } @@ -277,6 +293,15 @@ namespace internal + template + inline MappingInfoStorage & + MappingDataOnTheFly::get_data_storage() + { + return mapping_info_storage; + } + + + template inline const Quadrature<1> & MappingDataOnTheFly::get_quadrature() const