From: Peter Munch Date: Mon, 17 Jan 2022 19:04:14 +0000 (+0100) Subject: Move quadrature_points into FEEvalData X-Git-Tag: v9.4.0-rc1~586^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F13252%2Fhead;p=dealii.git 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 a5c4e384ba..7e4bcae909 100644 --- 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 - 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 - 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 FEEvaluationBasemapped_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 FEEvaluationBasemapped_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 @@ FEEvaluationcell_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 -inline Point -FEEvaluation::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 *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 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( - 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( - this->descriptor->quadrature.point(q)[e]); - return point; - } - else - return quadrature_points[q]; -} - - - template 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 @@ FEFaceEvaluationmapping_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::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 -inline Point -FEFaceEvaluation::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::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 + 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 *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::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::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::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::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::JxW(const unsigned int q_point) const +template +inline DEAL_II_ALWAYS_INLINE Point +FEEvaluationData::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 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( + 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( + this->descriptor->quadrature.point(q)[e]); + return point; + } + else + return this->quadrature_points[q]; +} + + + template inline DEAL_II_ALWAYS_INLINE Tensor<2, dim, Number> FEEvaluationData::inverse_jacobian(