From 670a598ef5b49eace842dd1c96f93ea519d10222 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Tue, 30 Nov 2021 22:37:52 +0100 Subject: [PATCH] Make evaluation_kernels only include FEEvaluationData file --- include/deal.II/base/array_view.h | 3 +- include/deal.II/matrix_free/dof_info.h | 8 +- .../deal.II/matrix_free/evaluation_kernels.h | 7 +- .../evaluation_template_factory.templates.h | 1 - ...template_factory_hanging_nodes.templates.h | 2 +- include/deal.II/matrix_free/fe_evaluation.h | 1714 ++++++----------- .../matrix_free/fe_evaluation_base_data.h | 763 ++++++++ include/deal.II/matrix_free/mapping_info.h | 304 +-- .../matrix_free/mapping_info.templates.h | 226 +-- .../matrix_free/mapping_info_storage.h | 346 ++++ .../mapping_info_storage.templates.h | 267 +++ 11 files changed, 1955 insertions(+), 1686 deletions(-) create mode 100644 include/deal.II/matrix_free/fe_evaluation_base_data.h create mode 100644 include/deal.II/matrix_free/mapping_info_storage.h create mode 100644 include/deal.II/matrix_free/mapping_info_storage.templates.h diff --git a/include/deal.II/base/array_view.h b/include/deal.II/base/array_view.h index f38a7d14ae..6a4e2017dd 100644 --- a/include/deal.II/base/array_view.h +++ b/include/deal.II/base/array_view.h @@ -414,7 +414,8 @@ inline void ArrayView::reinit(value_type *starting_element, const std::size_t n_elements) { - *this = ArrayView(starting_element, n_elements); + this->starting_element = starting_element; + this->n_elements = n_elements; } diff --git a/include/deal.II/matrix_free/dof_info.h b/include/deal.II/matrix_free/dof_info.h index 08178f7eab..7b0583ea9e 100644 --- a/include/deal.II/matrix_free/dof_info.h +++ b/include/deal.II/matrix_free/dof_info.h @@ -31,7 +31,6 @@ #include #include -#include #include #include #include @@ -49,7 +48,10 @@ namespace internal { template class HangingNodes; - } + + template + struct FPArrayComparator; + } // namespace MatrixFreeFunctions } // namespace internal #endif @@ -86,7 +88,7 @@ namespace internal std::pair, types::global_dof_index> next_constraint; std::map, types::global_dof_index, - FPArrayComparator> + FPArrayComparator>> constraints; }; diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index 0532e46907..dc84931d63 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -26,6 +26,7 @@ #include #include +#include #include #include #include @@ -35,12 +36,6 @@ DEAL_II_NAMESPACE_OPEN -// forward declaration -template -class FEEvaluationBaseData; - - - namespace internal { // Select evaluator type from element shape function type diff --git a/include/deal.II/matrix_free/evaluation_template_factory.templates.h b/include/deal.II/matrix_free/evaluation_template_factory.templates.h index 5a98c3781b..3f6272b615 100644 --- a/include/deal.II/matrix_free/evaluation_template_factory.templates.h +++ b/include/deal.II/matrix_free/evaluation_template_factory.templates.h @@ -24,7 +24,6 @@ #include #include #include -#include DEAL_II_NAMESPACE_OPEN diff --git a/include/deal.II/matrix_free/evaluation_template_factory_hanging_nodes.templates.h b/include/deal.II/matrix_free/evaluation_template_factory_hanging_nodes.templates.h index cf03c7c2d5..3c37bdd339 100644 --- a/include/deal.II/matrix_free/evaluation_template_factory_hanging_nodes.templates.h +++ b/include/deal.II/matrix_free/evaluation_template_factory_hanging_nodes.templates.h @@ -23,7 +23,7 @@ #include #include #include -#include +#include DEAL_II_NAMESPACE_OPEN diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index 4fc5c6cc2a..d74ef03107 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -33,6 +33,7 @@ #include #include #include +#include #include #include #include @@ -45,461 +46,6 @@ DEAL_II_NAMESPACE_OPEN -namespace internal -{ - DeclException0(ExcAccessToUninitializedField); - - DeclException1( - ExcMatrixFreeAccessToUninitializedMappingField, - std::string, - << "You are requesting information from an FEEvaluation/FEFaceEvaluation " - << "object for which this kind of information has not been computed. What " - << "information these objects compute is determined by the update_* flags you " - << "pass to MatrixFree::reinit() via MatrixFree::AdditionalData. Here, " - << "the operation you are attempting requires the <" << arg1 - << "> flag to be set, but it was apparently not specified " - << "upon initialization."); -} // namespace internal - -template > -class FEEvaluation; - - - -/** - * This base class of the FEEvaluation and FEFaceEvaluation classes handles - * mapping-related information independent of the degrees of freedom and - * finite element in use. This class provides access functionality for user - * code but is otherwise invisible without any public constructor. The usage - * is through the class FEEvaluation instead. - * - * This class has four template arguments: - * - * @tparam dim Dimension in which this class is to be used - * - * @tparam Number Number format, usually @p double or @p float - * - * @tparam is_face Whether the class is used for a cell integrator (with - * quadrature dimension the same as the space dimension) or for a face - * integrator (with quadrature dimension one less) - * - * @tparam VectorizedArrayType Type of array to be woked on in a vectorized - * fashion, defaults to VectorizedArray - * - * @note Currently only VectorizedArray is supported as - * VectorizedArrayType. - * - * - * @ingroup matrixfree - */ -template -class FEEvaluationBaseData -{ - static_assert( - std::is_same::value, - "Type of Number and of VectorizedArrayType do not match."); - -public: - static constexpr unsigned int dimension = dim; - - /** - * Destructor. - */ - ~FEEvaluationBaseData(); - - /** - * Return the index offset within the geometry fields for the cell the @p - * reinit() function has been called for. This index can be used to access - * an index into a field that has the same compression behavior as the - * Jacobian of the geometry, e.g., to store an effective coefficient tensors - * that combines a coefficient with the geometry for lower memory transfer - * as the available data fields. - */ - unsigned int - get_mapping_data_index_offset() const; - - /** - * Return the type of the cell the @p reinit() function has been called for. - * Valid values are @p cartesian for Cartesian cells (which allows for - * considerable data compression), @p affine for cells with affine mappings, - * and @p general for general cells without any compressed storage applied. - */ - internal::MatrixFreeFunctions::GeometryType - get_cell_type() const; - - /** - * Return a reference to the ShapeInfo object currently in use. - */ - const internal::MatrixFreeFunctions::ShapeInfo & - get_shape_info() const; - - /** - * Return a reference to the DoFInfo object currently in use. - */ - const internal::MatrixFreeFunctions::DoFInfo & - get_dof_info() const; - - /** - * Return the determinant of the Jacobian from the unit to the real cell - * times the quadrature weight. - */ - VectorizedArrayType - JxW(const unsigned int q_point) 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 - * $J_{ij} = d x_i / d\hat x_j$. The $(i,j)$ entry of the returned tensor - * contains $d\hat x_j/dx_i$, i.e., columns refer to reference space - * coordinates and rows to real cell coordinates. Thus, the returned tensor - * represents a covariant transformation, which is used in the - * FEEvaluationBase::get_gradient() function to transform the unit cell - * gradients to gradients on the real cell by a multiplication $J^{-\mathrm - * T} \hat{\nabla} u_h$. - */ - Tensor<2, dim, VectorizedArrayType> - inverse_jacobian(const unsigned int q_point) const; - - /** - * Return the unit normal vector on a face. Note that both sides of a face - * use the same orientation of the normal vector: For the faces enumerated - * as `interior` in FaceToCellTopology and selected with the - * `is_interior_face=true` flag of the constructor, this corresponds to the - * outer normal vector, whereas for faces enumerated as `exterior` in - * FaceToCellTopology and selected with the `is_interior_face=false` flag of - * the constructor, the normal points into the element as a consequence of - * the single normal vector. - * - * @note Only implemented in case `is_face == true`. - */ - Tensor<1, dim, VectorizedArrayType> - get_normal_vector(const unsigned int q_point) const; - - /** - * Provides a unified interface to access data in a vector of - * VectorizedArray fields of length MatrixFree::n_cell_batches() + - * MatrixFree::n_ghost_cell_batches() for both cells (plain read) and faces - * (indirect addressing). - */ - VectorizedArrayType - read_cell_data(const AlignedVector &array) const; - - /** - * Provides a unified interface to set data in a vector of - * VectorizedArray fields of length MatrixFree::n_cell_batches() + - * MatrixFree::n_ghost_cell_batches() for both cells (plain read) and faces - * (indirect addressing). - */ - void - set_cell_data(AlignedVector &array, - const VectorizedArrayType & value) const; - - /** - * The same as above, just for std::array of length of VectorizedArrayType for - * arbitrary data type. - */ - template - std::array - read_cell_data(const AlignedVector> - &array) const; - - /** - * The same as above, just for std::array of length of VectorizedArrayType for - * arbitrary data type. - */ - template - void - set_cell_data( - AlignedVector> &array, - const std::array & value) const; - - /** - * Return the id of the cells this FEEvaluation or FEFaceEvaluation is - * associated with. - */ - std::array - get_cell_ids() const; - - /** - * Return the id of the cells/faces this FEEvaluation/FEFaceEvaluation is - * associated with. - */ - std::array - get_cell_or_face_ids() const; - - /** - * Return the numbering of local degrees of freedom within the evaluation - * routines of FEEvaluation in terms of the standard numbering on finite - * elements. - */ - const std::vector & - get_internal_dof_numbering() const; - - /** - * Return an ArrayView to internal memory for temporary use. Note that some - * of this memory is overwritten during evaluate() and integrate() calls so - * do not assume it to be stable over those calls. The maximum size you can - * write into is 3*dofs_per_cell+2*n_q_points. - */ - ArrayView - get_scratch_data() const; - - /** - * Return the number of the quadrature formula of the present cell. - */ - unsigned int - get_quadrature_index() const; - - /** - * Return index of the current cell or face. - */ - unsigned int - get_current_cell_index() const; - - /** - * Return the active FE index for this class for efficient indexing in the - * hp-case. - */ - unsigned int - get_active_fe_index() const; - - /** - * Return the active quadrature index for this class for efficient indexing in - * the hp-case. - */ - unsigned int - get_active_quadrature_index() const; - - /** - * Return the underlying MatrixFree object. - */ - const MatrixFree & - get_matrix_free() const; - -protected: - /** - * Constructor. Made protected to prevent users from directly using this - * class. Takes all data stored in MatrixFree. If applied to problems with - * more than one quadrature formula selected during construction of - * `matrix_free`, `quad_no` allows to select the appropriate formula. - */ - FEEvaluationBaseData( - const MatrixFree &matrix_free, - const unsigned int dof_no, - const unsigned int first_selected_component, - const unsigned int quad_no, - const unsigned int fe_degree, - const unsigned int n_q_points, - const bool is_interior_face, - const unsigned int active_fe_index, - const unsigned int active_quad_index, - const unsigned int face_type); - - /** - * Constructor that comes with reduced functionality and works similar as - * FEValues. - */ - FEEvaluationBaseData( - const Mapping & mapping, - const FiniteElement &fe, - const Quadrature<1> & quadrature, - const UpdateFlags update_flags, - const unsigned int first_selected_component, - const FEEvaluationBaseData - *other); - - /** - * Copy constructor. If FEEvaluationBase was constructed from a mapping, fe, - * quadrature, and update flags, the underlying geometry evaluation based on - * FEValues will be deep-copied in order to allow for using in parallel with - * threads. - */ - FEEvaluationBaseData(const FEEvaluationBaseData &other); - - /** - * Copy assignment operator. If FEEvaluationBase was constructed from a - * mapping, fe, quadrature, and update flags, the underlying geometry - * evaluation based on FEValues will be deep-copied in order to allow for - * using in parallel with threads. - */ - FEEvaluationBaseData & - operator=(const FEEvaluationBaseData &other); - - /** - * This is the general array for all data fields. - */ - AlignedVector *scratch_data_array; - - /** - * This is the user-visible part of scratch_data_array, only showing the - * last part of scratch_data_array. The first part is consumed by - * values_dofs, values_quad, etc. - */ - VectorizedArrayType *scratch_data; - - /** - * The number of the quadrature formula of the present cell. - */ - const unsigned int quad_no; - - /** - * A pointer to the underlying data. - */ - const MatrixFree *matrix_info; - - /** - * A pointer to the underlying DoF indices and constraint description - * for the component specified at construction. Also contained in - * matrix_info, but it simplifies code if we store a reference to it. - */ - const internal::MatrixFreeFunctions::DoFInfo *dof_info; - - /** - * A pointer to the underlying transformation data from unit to real cells - * for the given quadrature formula specified at construction. Also - * contained in matrix_info, but it simplifies code if we store a reference - * to it. - */ - const internal::MatrixFreeFunctions::MappingInfoStorage< - (is_face ? dim - 1 : dim), - dim, - Number, - VectorizedArrayType> *mapping_data; - - /** - * The active FE index for this class for efficient indexing in the hp-case. - */ - const unsigned int active_fe_index; - - /** - * The active quadrature index for this class for efficient indexing in the - * hp-case. - */ - const unsigned int active_quad_index; - - /** - * A pointer to the underlying quadrature formula specified at construction. - * Also contained in matrix_info, but it simplifies code if we store a - * reference to it. - */ - const typename internal::MatrixFreeFunctions::MappingInfoStorage< - (is_face ? dim - 1 : dim), - dim, - Number, - VectorizedArrayType>::QuadratureDescriptor *descriptor; - - /** - * The number of quadrature points in the current evaluation context. - */ - const unsigned int n_quadrature_points; - - /** - * A pointer to the unit cell shape data, i.e., values, gradients and - * Hessians in 1D at the quadrature points that constitute the tensor - * product. Also contained in matrix_info, but it simplifies code if we - * store a reference to it. - */ - const internal::MatrixFreeFunctions::ShapeInfo *data; - - /** - * A pointer to the Jacobian information of the present cell. Only set to a - * useful value if on a non-Cartesian cell. - */ - const Tensor<2, dim, VectorizedArrayType> *jacobian; - - /** - * A pointer to the Jacobian determinant of the present cell. If on a - * Cartesian cell or on a cell with constant Jacobian, this is just the - * Jacobian determinant, otherwise the Jacobian determinant times the - * quadrature weight. - */ - const VectorizedArrayType *J_value; - - /** - * A pointer to the normal vectors at faces. - */ - const Tensor<1, dim, VectorizedArrayType> *normal_vectors; - - /** - * A pointer to the normal vectors times the jacobian at faces. - */ - const Tensor<1, dim, VectorizedArrayType> *normal_x_jacobian; - - /** - * A pointer to the quadrature weights of the underlying quadrature formula. - */ - const Number *quadrature_weights; - - /** - * After a call to reinit(), stores the number of the cell we are currently - * working with. - */ - unsigned int cell; - - /** - * Flag holding information whether a face is an interior or exterior face - * according to the defined direction of the normal. For cells it defines if - * the dof values should be read from the actual cell corresponding to the - * interior face or the neighboring cell corresponding to the exterior face. - */ - bool is_interior_face; - - /** - * Stores the index an FEFaceEvaluation object is currently pointing into - * (interior face, exterior face, data associated with cell). - */ - internal::MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index; - - /** - * Stores the current number of a face within the given cell in case - * `is_face==true`, using values between `0` and `2*dim`. - */ - unsigned int face_no; - - /** - * Stores the orientation of the given face with respect to the standard - * orientation, 0 if in standard orientation. - */ - unsigned int face_orientation; - - /** - * Stores the subface index of the given face. Usually, this variable takes - * the value numbers::invalid_unsigned_int to indicate integration over the - * full face, but in case the current physical face has a neighbor that is - * more refined, it is a subface and must scale the entries in ShapeInfo - * appropriately. - */ - unsigned int subface_index; - - /** - * Stores the type of the cell we are currently working with after a call to - * reinit(). Valid values are @p cartesian, @p affine and @p general, which - * have different implications on how the Jacobian transformations are - * stored internally in MappingInfo. - */ - internal::MatrixFreeFunctions::GeometryType cell_type; - - /** - * Geometry data that can be generated FEValues on the fly with the - * respective constructor. - */ - std::shared_ptr> - mapped_geometry; - - // Make FEEvaluation objects friends for access to protected member - // mapped_geometry. - template - friend class FEEvaluation; -}; - - - /** * This is the base class for the FEEvaluation classes. This class needs * usually not be called in user code and does not have any public @@ -540,8 +86,8 @@ protected: template > + bool is_face, + typename VectorizedArrayType> class FEEvaluationBase : public FEEvaluationBaseData { @@ -1143,12 +689,70 @@ public: //@} + /** + * Provides a unified interface to access data in a vector of + * VectorizedArray fields of length MatrixFree::n_cell_batches() + + * MatrixFree::n_ghost_cell_batches() for both cells (plain read) and faces + * (indirect addressing). + */ + VectorizedArrayType + read_cell_data(const AlignedVector &array) const; + + /** + * Provides a unified interface to set data in a vector of + * VectorizedArray fields of length MatrixFree::n_cell_batches() + + * MatrixFree::n_ghost_cell_batches() for both cells (plain read) and faces + * (indirect addressing). + */ + void + set_cell_data(AlignedVector &array, + const VectorizedArrayType & value) const; + + /** + * The same as above, just for std::array of length of VectorizedArrayType for + * arbitrary data type. + */ + template + std::array + read_cell_data(const AlignedVector> + &array) const; + + /** + * The same as above, just for std::array of length of VectorizedArrayType for + * arbitrary data type. + */ + template + void + set_cell_data( + AlignedVector> &array, + const std::array & value) const; + + /** + * Return the id of the cells this FEEvaluation or FEFaceEvaluation is + * associated with. + */ + std::array + get_cell_ids() const; + + /** + * Return the id of the cells/faces this FEEvaluation/FEFaceEvaluation is + * associated with. + */ + std::array + get_cell_or_face_ids() const; + /** * Return the first selected component. */ unsigned int get_first_selected_component() const; + /** + * Return the underlying MatrixFree object. + */ + const MatrixFree & + get_matrix_free() const; + protected: /** * Constructor. Made protected to prevent users from directly using this @@ -1232,6 +836,11 @@ protected: FEEvaluationBase & operator=(const FEEvaluationBase &other); + /** + * Destructor. + */ + ~FEEvaluationBase(); + /** * A unified function to read from and write into vectors based on the given * template operation. It can perform the operation for @p read_dof_values, @@ -1286,6 +895,16 @@ protected: void apply_hanging_node_constraints() const; + /** + * This is the general array for all data fields. + */ + AlignedVector *scratch_data_array; + + /** + * A pointer to the underlying data. + */ + const MatrixFree *matrix_info; + /** * This field stores the values for local degrees of freedom (e.g. after * reading out from a vector but before applying unit cell transformations @@ -3470,83 +3089,174 @@ namespace internal #ifndef DOXYGEN -/*----------------------- FEEvaluationBaseData ------------------------*/ +namespace internal +{ + template + inline unsigned int + get_active_fe_index(const MatrixFreeType &data, + const unsigned int dof_no, + const unsigned int fe_degree, + const unsigned int first_selected_component, + const unsigned int active_fe_index) + { + return fe_degree != numbers::invalid_unsigned_int ? + data.get_dof_info(dof_no).fe_index_from_degree( + first_selected_component, fe_degree) : + (active_fe_index != numbers::invalid_unsigned_int ? + active_fe_index : + 0); + } -template -inline FEEvaluationBaseData:: - FEEvaluationBaseData( - const MatrixFree &data_in, - const unsigned int dof_no, - const unsigned int first_selected_component, - const unsigned int quad_no_in, - const unsigned int fe_degree, - const unsigned int n_q_points, - const bool is_interior_face, - const unsigned int active_fe_index_in, - const unsigned int active_quad_index_in, - const unsigned int face_type) - : scratch_data_array(data_in.acquire_scratch_data()) - , quad_no(quad_no_in) - , matrix_info(&data_in) - , dof_info(&data_in.get_dof_info(dof_no)) - , mapping_data( + + + template + inline unsigned int + get_active_quad_index(const MatrixFree &mf, + const unsigned int quad_no, + const unsigned int fe_degree, + const unsigned int n_q_points, + const unsigned int active_fe_index, + const unsigned int active_quad_index) + { + const auto &mapping_storage = internal::MatrixFreeFunctions:: + MappingInfoCellsOrFaces::get( + mf.get_mapping_info(), quad_no); + if (fe_degree == numbers::invalid_unsigned_int) + return active_quad_index != numbers::invalid_unsigned_int ? + active_quad_index : + std::min(active_fe_index, + mapping_storage.descriptor.size() - 1); + else + return mapping_storage.quad_index_from_n_q_points(n_q_points); + } +} // namespace internal + + +/*----------------------- FEEvaluationBase ----------------------------------*/ + +template +inline FEEvaluationBase:: + FEEvaluationBase(const MatrixFree &data_in, + const unsigned int dof_no, + const unsigned int first_selected_component, + const unsigned int quad_no, + const unsigned int fe_degree, + const unsigned int n_q_points, + const bool is_interior_face, + const unsigned int active_fe_index, + const unsigned int active_quad_index, + const unsigned int face_type) + : FEEvaluationBaseData( + data_in.get_shape_info( + dof_no, + quad_no, + data_in.get_dof_info(dof_no) + .component_to_base_index[first_selected_component], + internal::get_active_fe_index(data_in, + dof_no, + fe_degree, + first_selected_component, + active_fe_index), + internal::get_active_quad_index( + data_in, + quad_no, + fe_degree, + n_q_points, + internal::get_active_fe_index(data_in, + dof_no, + fe_degree, + first_selected_component, + active_fe_index), + active_quad_index)), + data_in.get_dof_info(dof_no), internal::MatrixFreeFunctions:: MappingInfoCellsOrFaces::get( data_in.get_mapping_info(), - quad_no)) - , active_fe_index(fe_degree != numbers::invalid_unsigned_int ? - data_in.get_dof_info(dof_no).fe_index_from_degree( - first_selected_component, - fe_degree) : - (active_fe_index_in != numbers::invalid_unsigned_int ? - active_fe_index_in : - 0)) - , active_quad_index( - fe_degree != numbers::invalid_unsigned_int ? - (mapping_data->quad_index_from_n_q_points(n_q_points)) : - (active_quad_index_in != numbers::invalid_unsigned_int ? - active_quad_index_in : - std::min(active_fe_index, - mapping_data->descriptor.size() - 1))) - , descriptor( - &mapping_data->descriptor - [is_face ? - (active_quad_index * std::max(1, dim - 1) + - (face_type == numbers::invalid_unsigned_int ? 0 : face_type)) : - active_quad_index]) - , n_quadrature_points(descriptor->n_q_points) - , data(&data_in.get_shape_info( - dof_no, - quad_no_in, - dof_info->component_to_base_index[first_selected_component], - active_fe_index, - active_quad_index)) - , jacobian(nullptr) - , J_value(nullptr) - , normal_vectors(nullptr) - , normal_x_jacobian(nullptr) - , quadrature_weights(descriptor->quadrature_weights.begin()) - , cell(numbers::invalid_unsigned_int) - , is_interior_face(is_interior_face) - , dof_access_index( - is_face ? - (is_interior_face ? - internal::MatrixFreeFunctions::DoFInfo::dof_access_face_interior : - internal::MatrixFreeFunctions::DoFInfo::dof_access_face_exterior) : - internal::MatrixFreeFunctions::DoFInfo::dof_access_cell) - , cell_type(internal::MatrixFreeFunctions::general) -{ - Assert(matrix_info->mapping_initialized() == true, ExcNotInitialized()); - AssertDimension(matrix_info->get_task_info().vectorization_length, - VectorizedArrayType::size()); - AssertDimension(n_quadrature_points, descriptor->n_q_points); + quad_no), + quad_no, + is_interior_face, + internal::get_active_fe_index(data_in, + dof_no, + fe_degree, + first_selected_component, + active_fe_index), + internal::get_active_quad_index( + data_in, + quad_no, + fe_degree, + n_q_points, + internal::get_active_fe_index(data_in, + dof_no, + fe_degree, + first_selected_component, + active_fe_index), + active_quad_index), + face_type) + , scratch_data_array(data_in.acquire_scratch_data()) + , matrix_info(&data_in) + , n_fe_components(data_in.get_dof_info(dof_no).start_components.back()) + , dof_values_initialized(false) + , values_quad_initialized(false) + , gradients_quad_initialized(false) + , hessians_quad_initialized(false) + , values_quad_submitted(false) + , gradients_quad_submitted(false) + , first_selected_component(first_selected_component) +{ + set_data_pointers(); + Assert( + this->dof_info->start_components.back() == 1 || + static_cast(n_components_) <= + static_cast( + this->dof_info->start_components + [this->dof_info->component_to_base_index[first_selected_component] + + 1]) - + first_selected_component, + ExcMessage( + "You tried to construct a vector-valued evaluator with " + + std::to_string(n_components) + + " components. However, " + "the current base element has only " + + std::to_string( + this->dof_info->start_components + [this->dof_info->component_to_base_index[first_selected_component] + + 1] - + first_selected_component) + + " components left when starting from local element index " + + std::to_string( + first_selected_component - + this->dof_info->start_components + [this->dof_info->component_to_base_index[first_selected_component]]) + + " (global index " + std::to_string(first_selected_component) + ")")); + + // do not check for correct dimensions of data fields here, should be done + // in derived classes } -template -inline FEEvaluationBaseData:: - FEEvaluationBaseData( +template +inline FEEvaluationBase:: + FEEvaluationBase( const Mapping & mapping, const FiniteElement &fe, const Quadrature<1> & quadrature, @@ -3554,115 +3264,129 @@ inline FEEvaluationBaseData:: const unsigned int first_selected_component, const FEEvaluationBaseData *other) - : scratch_data_array(new AlignedVector()) - , quad_no(numbers::invalid_unsigned_int) - , active_fe_index(numbers::invalid_unsigned_int) - , active_quad_index(numbers::invalid_unsigned_int) - , descriptor(nullptr) - , n_quadrature_points( - Utilities::fixed_power < is_face ? dim - 1 : dim > (quadrature.size())) + : FEEvaluationBaseData( + other != nullptr && + other->mapped_geometry->get_quadrature() == quadrature ? + other->mapped_geometry : + std::make_shared< + internal::MatrixFreeFunctions:: + MappingDataOnTheFly>( + mapping, + quadrature, + update_flags)) + , scratch_data_array(new AlignedVector()) , matrix_info(nullptr) - , dof_info(nullptr) - , mapping_data(nullptr) - , - // select the correct base element from the given FE component - data(new internal::MatrixFreeFunctions::ShapeInfo( - Quadrature(quadrature), - fe, - fe.component_to_base_index(first_selected_component).first)) - , jacobian(nullptr) - , J_value(nullptr) - , normal_vectors(nullptr) - , normal_x_jacobian(nullptr) - , quadrature_weights(nullptr) - , cell(0) - , cell_type(internal::MatrixFreeFunctions::general) - , is_interior_face(true) - , dof_access_index(internal::MatrixFreeFunctions::DoFInfo::dof_access_cell) -{ - Assert(other == nullptr || other->mapped_geometry.get() != nullptr, - ExcInternalError()); - if (other != nullptr && - other->mapped_geometry->get_quadrature() == quadrature) - mapped_geometry = other->mapped_geometry; - else - mapped_geometry = - std::make_shared>( - mapping, quadrature, update_flags); - cell = 0; + , n_fe_components(n_components_) + , dof_values_initialized(false) + , values_quad_initialized(false) + , gradients_quad_initialized(false) + , hessians_quad_initialized(false) + , values_quad_submitted(false) + , gradients_quad_submitted(false) + // keep the number of the selected component within the current base element + // for reading dof values + , first_selected_component(first_selected_component) +{ + const unsigned int base_element_number = + fe.component_to_base_index(first_selected_component).first; + Assert(fe.element_multiplicity(base_element_number) == 1 || + fe.element_multiplicity(base_element_number) - + first_selected_component >= + n_components_, + ExcMessage("The underlying element must at least contain as many " + "components as requested by this class")); + (void)base_element_number; + + Assert(this->data == nullptr, ExcInternalError()); + this->data = + new internal::MatrixFreeFunctions::ShapeInfo( + Quadrature(quadrature), + fe, + fe.component_to_base_index(first_selected_component).first); - 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(); + set_data_pointers(); } -template -inline FEEvaluationBaseData:: - FEEvaluationBaseData( - const FEEvaluationBaseData - &other) - : scratch_data_array(other.matrix_info == nullptr ? +template +inline FEEvaluationBase:: + FEEvaluationBase(const FEEvaluationBase &other) + : FEEvaluationBaseData(other) + , scratch_data_array(other.matrix_info == nullptr ? new AlignedVector() : other.matrix_info->acquire_scratch_data()) - , quad_no(other.quad_no) - , active_fe_index(other.active_fe_index) - , active_quad_index(other.active_quad_index) - , descriptor(other.descriptor == nullptr ? nullptr : other.descriptor) - , n_quadrature_points(other.n_quadrature_points) , matrix_info(other.matrix_info) - , dof_info(other.dof_info) - , mapping_data(other.mapping_data) - , data(other.matrix_info == nullptr ? - new internal::MatrixFreeFunctions::ShapeInfo( - *other.data) : - other.data) - , jacobian(nullptr) - , J_value(nullptr) - , normal_vectors(nullptr) - , normal_x_jacobian(nullptr) - , quadrature_weights(other.descriptor == nullptr ? - nullptr : - descriptor->quadrature_weights.begin()) - , cell(numbers::invalid_unsigned_int) - , cell_type(internal::MatrixFreeFunctions::general) - , is_interior_face(other.is_interior_face) - , dof_access_index(other.dof_access_index) -{ - // Create deep copy of mapped geometry for use in parallel... - if (other.mapped_geometry.get() != nullptr) + , n_fe_components(other.n_fe_components) + , dof_values_initialized(false) + , values_quad_initialized(false) + , gradients_quad_initialized(false) + , hessians_quad_initialized(false) + , values_quad_submitted(false) + , gradients_quad_submitted(false) + , first_selected_component(other.first_selected_component) +{ + if (other.matrix_info == nullptr) { - mapped_geometry = std::make_shared< + Assert(other.mapped_geometry.get() != nullptr, ExcInternalError()); + this->data = + new internal::MatrixFreeFunctions::ShapeInfo( + *other.data); + + // Create deep copy of mapped geometry for use in parallel + this->mapped_geometry = std::make_shared< internal::MatrixFreeFunctions:: MappingDataOnTheFly>( other.mapped_geometry->get_fe_values().get_mapping(), other.mapped_geometry->get_quadrature(), other.mapped_geometry->get_fe_values().get_update_flags()); - mapping_data = &mapped_geometry->get_data_storage(); - cell = 0; + this->mapping_data = &this->mapped_geometry->get_data_storage(); + this->cell = 0; - jacobian = mapped_geometry->get_data_storage().jacobians[0].begin(); - J_value = mapped_geometry->get_data_storage().JxW_values.begin(); + this->jacobian = + this->mapped_geometry->get_data_storage().jacobians[0].begin(); + this->J_value = + this->mapped_geometry->get_data_storage().JxW_values.begin(); } + + set_data_pointers(); } -template -inline FEEvaluationBaseData & -FEEvaluationBaseData::operator=( - const FEEvaluationBaseData &other) +template +inline FEEvaluationBase & +FEEvaluationBase:: +operator=(const FEEvaluationBase &other) { - AssertDimension(quad_no, other.quad_no); - AssertDimension(active_fe_index, other.active_fe_index); - AssertDimension(active_quad_index, other.active_quad_index); - // release old memory if (matrix_info == nullptr) { - delete data; + delete this->data; delete scratch_data_array; } else @@ -3670,182 +3394,156 @@ FEEvaluationBaseData::operator=( matrix_info->release_scratch_data(scratch_data_array); } - matrix_info = other.matrix_info; - dof_info = other.dof_info; - descriptor = other.descriptor; - mapping_data = other.mapping_data; + this->FEEvaluationBaseData:: + operator=(other); + + matrix_info = other.matrix_info; + if (other.matrix_info == nullptr) { - data = new internal::MatrixFreeFunctions::ShapeInfo( - *other.data); + Assert(other.mapped_geometry.get() != nullptr, ExcInternalError()); + this->data = + new internal::MatrixFreeFunctions::ShapeInfo( + *other.data); scratch_data_array = new AlignedVector(); - } - else - { - data = other.data; - scratch_data_array = matrix_info->acquire_scratch_data(); - } - quadrature_weights = - (descriptor != nullptr ? descriptor->quadrature_weights.begin() : nullptr); - cell = numbers::invalid_unsigned_int; - cell_type = internal::MatrixFreeFunctions::general; - is_interior_face = other.is_interior_face; - dof_access_index = other.dof_access_index; - - // Create deep copy of mapped geometry for use in parallel... - if (other.mapped_geometry.get() != nullptr) - { - mapped_geometry = std::make_shared< + // Create deep copy of mapped geometry for use in parallel + this->mapped_geometry = std::make_shared< internal::MatrixFreeFunctions:: MappingDataOnTheFly>( other.mapped_geometry->get_fe_values().get_mapping(), other.mapped_geometry->get_quadrature(), other.mapped_geometry->get_fe_values().get_update_flags()); - cell = 0; - 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(); - } - - return *this; -} - - - -template -inline FEEvaluationBaseData:: - ~FEEvaluationBaseData() -{ - if (matrix_info != nullptr) - { - try - { - matrix_info->release_scratch_data(scratch_data_array); - } - catch (...) - {} - } - else - { - delete scratch_data_array; - delete data; - data = nullptr; + this->cell = 0; + this->mapping_data = &this->mapped_geometry->get_data_storage(); + this->jacobian = + this->mapped_geometry->get_data_storage().jacobians[0].begin(); + this->J_value = + this->mapped_geometry->get_data_storage().JxW_values.begin(); } - scratch_data_array = nullptr; -} - - - -template -inline unsigned int -FEEvaluationBaseData:: - get_mapping_data_index_offset() const -{ - if (matrix_info == nullptr) - return 0; else { - AssertIndexRange(cell, this->mapping_data->data_index_offsets.size()); - return this->mapping_data->data_index_offsets[cell]; + scratch_data_array = matrix_info->acquire_scratch_data(); } -} - - - -template -inline internal::MatrixFreeFunctions::GeometryType -FEEvaluationBaseData::get_cell_type() - const -{ - Assert(cell != numbers::invalid_unsigned_int, ExcNotInitialized()); - return cell_type; -} + AssertDimension(n_fe_components, other.n_fe_components); + AssertDimension(first_selected_component, other.first_selected_component); + set_data_pointers(); -template -inline const internal::MatrixFreeFunctions::ShapeInfo & -FEEvaluationBaseData:: - get_shape_info() const -{ - Assert(data != nullptr, ExcInternalError()); - return *data; + return *this; } -template -inline const internal::MatrixFreeFunctions::DoFInfo & -FEEvaluationBaseData::get_dof_info() - const +template +inline FEEvaluationBase::~FEEvaluationBase() { - Assert(dof_info != nullptr, - ExcMessage( - "FEEvaluation was not initialized with a MatrixFree object!")); - return *dof_info; + if (matrix_info != nullptr) + { + try + { + matrix_info->release_scratch_data(scratch_data_array); + } + catch (...) + {} + } + else + { + delete scratch_data_array; + delete this->data; + } } -template -inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, VectorizedArrayType> -FEEvaluationBaseData:: - get_normal_vector(const unsigned int q_point) const +template +inline void +FEEvaluationBase:: + set_data_pointers() { - AssertIndexRange(q_point, n_quadrature_points); - Assert(normal_vectors != nullptr, - internal::ExcMatrixFreeAccessToUninitializedMappingField( - "update_normal_vectors")); - if (this->cell_type <= internal::MatrixFreeFunctions::flat_faces) - return normal_vectors[0]; - else - return normal_vectors[q_point]; -} + Assert(this->scratch_data_array != nullptr, ExcInternalError()); + const unsigned int tensor_dofs_per_component = + Utilities::fixed_power(this->data->data.front().fe_degree + 1); + const unsigned int dofs_per_component = + this->data->dofs_per_component_on_cell; + const unsigned int n_quadrature_points = this->n_quadrature_points; + const unsigned int size_scratch_data = + std::max(tensor_dofs_per_component + 1, dofs_per_component) * + n_components_ * 3 + + 2 * n_quadrature_points; + const unsigned int allocated_size = + size_scratch_data + n_components_ * dofs_per_component + + (n_components_ * ((dim * (dim + 1)) / 2 + 2 * dim + 1) * + n_quadrature_points); + this->scratch_data_array->resize_fast(allocated_size); -template -inline DEAL_II_ALWAYS_INLINE VectorizedArrayType -FEEvaluationBaseData::JxW( - const unsigned int q_point) const -{ - AssertIndexRange(q_point, n_quadrature_points); - Assert(J_value != nullptr, - internal::ExcMatrixFreeAccessToUninitializedMappingField( - "update_values|update_gradients")); - if (this->cell_type <= internal::MatrixFreeFunctions::affine) + // set the pointers to the correct position in the data array + for (unsigned int c = 0; c < n_components_; ++c) { - Assert(this->quadrature_weights != nullptr, ExcInternalError()); - return J_value[0] * this->quadrature_weights[q_point]; + values_dofs[c] = + this->scratch_data_array->begin() + c * dofs_per_component; } - else - return J_value[q_point]; + values_quad = + this->scratch_data_array->begin() + n_components * dofs_per_component; + gradients_quad = this->scratch_data_array->begin() + + n_components * (dofs_per_component + n_quadrature_points); + gradients_from_hessians_quad = + this->scratch_data_array->begin() + + n_components * (dofs_per_component + (dim + 1) * n_quadrature_points); + hessians_quad = + this->scratch_data_array->begin() + + n_components * (dofs_per_component + (2 * dim + 1) * n_quadrature_points); + + this->scratch_data.reinit(this->scratch_data_array->begin() + + n_components_ * dofs_per_component + + (n_components_ * + ((dim * (dim + 1)) / 2 + 2 * dim + 1) * + n_quadrature_points), + size_scratch_data); } -template -inline DEAL_II_ALWAYS_INLINE Tensor<2, dim, VectorizedArrayType> -FEEvaluationBaseData:: - inverse_jacobian(const unsigned int q_point) const +template +inline const MatrixFree & +FEEvaluationBase:: + get_matrix_free() const { - AssertIndexRange(q_point, n_quadrature_points); - Assert(this->jacobian != nullptr, - internal::ExcMatrixFreeAccessToUninitializedMappingField( - "update_gradients")); - if (this->cell_type <= internal::MatrixFreeFunctions::affine) - return jacobian[0]; - else - return jacobian[q_point]; + Assert(matrix_info != nullptr, + ExcMessage( + "FEEvaluation was not initialized with a MatrixFree object!")); + return *matrix_info; } -template +template inline std::array -FEEvaluationBaseData::get_cell_ids() - const +FEEvaluationBase:: + get_cell_ids() const { Assert(this->matrix_info != nullptr, ExcNotInitialized()); @@ -3864,7 +3562,7 @@ FEEvaluationBaseData::get_cell_ids() { // cell or interior face face (element-centric loop) for (unsigned int i = 0; i < n_lanes; ++i) - cells[i] = cell * n_lanes + i; + cells[i] = this->cell * n_lanes + i; } else if (is_face && this->dof_access_index == @@ -3905,9 +3603,9 @@ FEEvaluationBaseData::get_cell_ids() { // face-centric faces const unsigned int *cells_ = - is_interior_face ? - &this->matrix_info->get_face_info(cell).cells_interior[0] : - &this->matrix_info->get_face_info(cell).cells_exterior[0]; + this->is_interior_face ? + &this->matrix_info->get_face_info(this->cell).cells_interior[0] : + &this->matrix_info->get_face_info(this->cell).cells_exterior[0]; for (unsigned int i = 0; i < VectorizedArrayType::size(); ++i) if (cells_[i] != numbers::invalid_unsigned_int) cells[i] = cells_[i]; @@ -3920,6 +3618,7 @@ FEEvaluationBaseData::get_cell_ids() namespace internal { template inline void process_cell_data( - const FEEvaluationBaseData &phi, + const FEEvaluationBase & phi, const MatrixFree *matrix_info, GlobalVectorType & array, VectorizedArrayType2 & out, @@ -3946,277 +3649,69 @@ namespace internal for (unsigned int i = 0; i < VectorizedArrayType::size(); ++i) if (cells[i] != numbers::invalid_unsigned_int) fu(out[i], - array[cells[i] / VectorizedArrayType::size()] - [cells[i] % VectorizedArrayType::size()]); - } -} // namespace internal - - - -template -std::array -FEEvaluationBaseData:: - get_cell_or_face_ids() const -{ - const unsigned int v_len = VectorizedArrayType::size(); - std::array cells; - - // initialize array - for (unsigned int i = 0; i < v_len; ++i) - cells[i] = numbers::invalid_unsigned_int; - - if ((is_face && - this->dof_access_index == - internal::MatrixFreeFunctions::DoFInfo::dof_access_cell && - this->is_interior_face == false) || - (!is_face && !this->is_interior_face)) - { - // cell-based face-loop: plus face - for (unsigned int i = 0; i < v_len; ++i) - { - // compute actual (non vectorized) cell ID - const unsigned int cell_this = this->cell * v_len + i; - // compute face ID - unsigned int fn = - this->matrix_info->get_cell_and_face_to_plain_faces()(this->cell, - this->face_no, - i); - - if (fn == numbers::invalid_unsigned_int) - continue; // invalid face ID: no neighbor on boundary - - // get cell ID on both sides of face - auto cell_m = this->matrix_info->get_face_info(fn / v_len) - .cells_interior[fn % v_len]; - auto cell_p = this->matrix_info->get_face_info(fn / v_len) - .cells_exterior[fn % v_len]; - - // compare the IDs with the given cell ID - if (cell_m == cell_this) - cells[i] = cell_p; // neighbor has the other ID - else if (cell_p == cell_this) - cells[i] = cell_m; - } - } - else - { - for (unsigned int i = 0; i < v_len; ++i) - cells[i] = cell * v_len + i; - } - - return cells; -} - - - -template -inline VectorizedArrayType -FEEvaluationBaseData::read_cell_data( - const AlignedVector &array) const -{ - VectorizedArrayType out = Number(1.); - internal::process_cell_data( - *this, this->matrix_info, array, out, [](auto &local, const auto &global) { - local = global; - }); - return out; -} - - - -template -inline void -FEEvaluationBaseData::set_cell_data( - AlignedVector &array, - const VectorizedArrayType & in) const -{ - internal::process_cell_data( - *this, this->matrix_info, array, in, [](const auto &local, auto &global) { - global = local; - }); -} - - - -template -template -inline std::array -FEEvaluationBaseData::read_cell_data( - const AlignedVector> &array) const -{ - std::array out; - internal::process_cell_data( - *this, this->matrix_info, array, out, [](auto &local, const auto &global) { - local = global; - }); - return out; -} - - - -template -template -inline void -FEEvaluationBaseData::set_cell_data( - AlignedVector> &array, - const std::array & in) const -{ - internal::process_cell_data( - *this, this->matrix_info, array, in, [](const auto &local, auto &global) { - global = local; - }); -} - - - -template -inline const std::vector & -FEEvaluationBaseData:: - get_internal_dof_numbering() const -{ - return data->lexicographic_numbering; -} - - - -template -inline ArrayView -FEEvaluationBaseData:: - get_scratch_data() const -{ - return ArrayView( - const_cast(scratch_data), - scratch_data_array->end() - scratch_data); -} - - - -template -inline unsigned int -FEEvaluationBaseData:: - get_quadrature_index() const -{ - return this->quad_no; -} - - - -template -inline unsigned int -FEEvaluationBaseData:: - get_current_cell_index() const -{ - if (is_face && this->dof_access_index == - internal::MatrixFreeFunctions::DoFInfo::dof_access_cell) - return this->cell * GeometryInfo::faces_per_cell + this->face_no; - else - return this->cell; -} - - - -template -inline unsigned int -FEEvaluationBaseData:: - get_active_fe_index() const -{ - return active_fe_index; -} - - - -template -inline unsigned int -FEEvaluationBaseData:: - get_active_quadrature_index() const -{ - return active_quad_index; -} - - - -template -inline const MatrixFree & -FEEvaluationBaseData:: - get_matrix_free() const -{ - Assert(matrix_info != nullptr, - ExcMessage( - "FEEvaluation was not initialized with a MatrixFree object!")); - return *matrix_info; -} + array[cells[i] / VectorizedArrayType::size()] + [cells[i] % VectorizedArrayType::size()]); + } +} // namespace internal -/*----------------------- FEEvaluationBase ----------------------------------*/ template -inline FEEvaluationBase:: - FEEvaluationBase(const MatrixFree &data_in, - const unsigned int dof_no, - const unsigned int first_selected_component, - const unsigned int quad_no_in, - const unsigned int fe_degree, - const unsigned int n_q_points, - const bool is_interior_face, - const unsigned int active_fe_index, - const unsigned int active_quad_index, - const unsigned int face_type) - : FEEvaluationBaseData( - data_in, - dof_no, - first_selected_component, - quad_no_in, - fe_degree, - n_q_points, - is_interior_face, - active_fe_index, - active_quad_index, - face_type) - , n_fe_components(data_in.get_dof_info(dof_no).start_components.back()) - , dof_values_initialized(false) - , values_quad_initialized(false) - , gradients_quad_initialized(false) - , hessians_quad_initialized(false) - , values_quad_submitted(false) - , gradients_quad_submitted(false) - , first_selected_component(first_selected_component) +std::array +FEEvaluationBase:: + get_cell_or_face_ids() const { - set_data_pointers(); - Assert( - this->dof_info->start_components.back() == 1 || - static_cast(n_components_) <= - static_cast( - this->dof_info->start_components - [this->dof_info->component_to_base_index[first_selected_component] + - 1]) - - first_selected_component, - ExcMessage( - "You tried to construct a vector-valued evaluator with " + - std::to_string(n_components) + - " components. However, " - "the current base element has only " + - std::to_string( - this->dof_info->start_components - [this->dof_info->component_to_base_index[first_selected_component] + - 1] - - first_selected_component) + - " components left when starting from local element index " + - std::to_string( - first_selected_component - - this->dof_info->start_components - [this->dof_info->component_to_base_index[first_selected_component]]) + - " (global index " + std::to_string(first_selected_component) + ")")); + const unsigned int v_len = VectorizedArrayType::size(); + std::array cells; - // do not check for correct dimensions of data fields here, should be done - // in derived classes + // initialize array + for (unsigned int i = 0; i < v_len; ++i) + cells[i] = numbers::invalid_unsigned_int; + + if ((is_face && + this->dof_access_index == + internal::MatrixFreeFunctions::DoFInfo::dof_access_cell && + this->is_interior_face == false) || + (!is_face && !this->is_interior_face)) + { + // cell-based face-loop: plus face + for (unsigned int i = 0; i < v_len; ++i) + { + // compute actual (non vectorized) cell ID + const unsigned int cell_this = this->cell * v_len + i; + // compute face ID + unsigned int fn = + this->matrix_info->get_cell_and_face_to_plain_faces()(this->cell, + this->face_no, + i); + + if (fn == numbers::invalid_unsigned_int) + continue; // invalid face ID: no neighbor on boundary + + // get cell ID on both sides of face + auto cell_m = this->matrix_info->get_face_info(fn / v_len) + .cells_interior[fn % v_len]; + auto cell_p = this->matrix_info->get_face_info(fn / v_len) + .cells_exterior[fn % v_len]; + + // compare the IDs with the given cell ID + if (cell_m == cell_this) + cells[i] = cell_p; // neighbor has the other ID + else if (cell_p == cell_this) + cells[i] = cell_m; + } + } + else + { + for (unsigned int i = 0; i < v_len; ++i) + cells[i] = this->cell * v_len + i; + } + + return cells; } @@ -4226,48 +3721,16 @@ template -inline FEEvaluationBase:: - FEEvaluationBase( - const Mapping & mapping, - const FiniteElement &fe, - const Quadrature<1> & quadrature, - const UpdateFlags update_flags, - const unsigned int first_selected_component, - const FEEvaluationBaseData - *other) - : FEEvaluationBaseData( - mapping, - fe, - quadrature, - update_flags, - first_selected_component, - other) - , n_fe_components(n_components_) - , dof_values_initialized(false) - , values_quad_initialized(false) - , gradients_quad_initialized(false) - , hessians_quad_initialized(false) - , values_quad_submitted(false) - , gradients_quad_submitted(false) - // keep the number of the selected component within the current base element - // for reading dof values - , first_selected_component(first_selected_component) +inline VectorizedArrayType +FEEvaluationBase:: + read_cell_data(const AlignedVector &array) const { - set_data_pointers(); - - const unsigned int base_element_number = - fe.component_to_base_index(first_selected_component).first; - Assert(fe.element_multiplicity(base_element_number) == 1 || - fe.element_multiplicity(base_element_number) - - first_selected_component >= - n_components_, - ExcMessage("The underlying element must at least contain as many " - "components as requested by this class")); - (void)base_element_number; + VectorizedArrayType out = Number(1.); + internal::process_cell_data( + *this, this->matrix_info, array, out, [](auto &local, const auto &global) { + local = global; + }); + return out; } @@ -4277,27 +3740,15 @@ template -inline FEEvaluationBase:: - FEEvaluationBase(const FEEvaluationBase &other) - : FEEvaluationBaseData(other) - , n_fe_components(other.n_fe_components) - , dof_values_initialized(false) - , values_quad_initialized(false) - , gradients_quad_initialized(false) - , hessians_quad_initialized(false) - , values_quad_submitted(false) - , gradients_quad_submitted(false) - , first_selected_component(other.first_selected_component) +inline void +FEEvaluationBase:: + set_cell_data(AlignedVector &array, + const VectorizedArrayType & in) const { - set_data_pointers(); + internal::process_cell_data( + *this, this->matrix_info, array, in, [](const auto &local, auto &global) { + global = local; + }); } @@ -4307,24 +3758,18 @@ template -inline FEEvaluationBase & +template +inline std::array FEEvaluationBase:: -operator=(const FEEvaluationBase &other) + read_cell_data(const AlignedVector> + &array) const { - this->FEEvaluationBaseData:: - operator=(other); - AssertDimension(n_fe_components, other.n_fe_components); - AssertDimension(first_selected_component, other.first_selected_component); - - return *this; + std::array out; + internal::process_cell_data( + *this, this->matrix_info, array, out, [](auto &local, const auto &global) { + local = global; + }); + return out; } @@ -4334,48 +3779,17 @@ template +template inline void FEEvaluationBase:: - set_data_pointers() + set_cell_data( + AlignedVector> &array, + const std::array & in) const { - Assert(this->scratch_data_array != nullptr, ExcInternalError()); - - const unsigned int tensor_dofs_per_component = - Utilities::fixed_power(this->data->data.front().fe_degree + 1); - const unsigned int dofs_per_component = - this->data->dofs_per_component_on_cell; - const unsigned int n_quadrature_points = this->n_quadrature_points; - - const unsigned int shift = - std::max(tensor_dofs_per_component + 1, dofs_per_component) * - n_components_ * 3 + - 2 * n_quadrature_points; - const unsigned int allocated_size = - shift + n_components_ * dofs_per_component + - (n_components_ * ((dim * (dim + 1)) / 2 + 2 * dim + 1) * - n_quadrature_points); - this->scratch_data_array->resize_fast(allocated_size); - - // set the pointers to the correct position in the data array - for (unsigned int c = 0; c < n_components_; ++c) - { - values_dofs[c] = - this->scratch_data_array->begin() + c * dofs_per_component; - } - values_quad = - this->scratch_data_array->begin() + n_components * dofs_per_component; - gradients_quad = this->scratch_data_array->begin() + - n_components * (dofs_per_component + n_quadrature_points); - gradients_from_hessians_quad = - this->scratch_data_array->begin() + - n_components * (dofs_per_component + (dim + 1) * n_quadrature_points); - hessians_quad = - this->scratch_data_array->begin() + - n_components * (dofs_per_component + (2 * dim + 1) * n_quadrature_points); - this->scratch_data = this->scratch_data_array->begin() + - n_components_ * dofs_per_component + - (n_components_ * ((dim * (dim + 1)) / 2 + 2 * dim + 1) * - n_quadrature_points); + internal::process_cell_data( + *this, this->matrix_info, array, in, [](const auto &local, auto &global) { + global = local; + }); } @@ -6571,7 +5985,7 @@ inline FEEvaluationAccess &data_in, const unsigned int dof_no, const unsigned int first_selected_component, - const unsigned int quad_no_in, + const unsigned int quad_no, const unsigned int fe_degree, const unsigned int n_q_points, const bool is_interior_face, @@ -6582,7 +5996,7 @@ inline FEEvaluationAccess:: const MatrixFree &data_in, const unsigned int dof_no, const unsigned int first_selected_component, - const unsigned int quad_no_in, + const unsigned int quad_no, const unsigned int fe_degree, const unsigned int n_q_points, const bool is_interior_face, @@ -6690,7 +6104,7 @@ inline FEEvaluationAccess:: data_in, dof_no, first_selected_component, - quad_no_in, + quad_no, fe_degree, n_q_points, is_interior_face, @@ -7006,7 +6420,7 @@ inline FEEvaluationAccess:: const MatrixFree &data_in, const unsigned int dof_no, const unsigned int first_selected_component, - const unsigned int quad_no_in, + const unsigned int quad_no, const unsigned int fe_degree, const unsigned int n_q_points, const bool is_interior_face, @@ -7017,7 +6431,7 @@ inline FEEvaluationAccess:: data_in, dof_no, first_selected_component, - quad_no_in, + quad_no, fe_degree, n_q_points, is_interior_face, @@ -7417,7 +6831,7 @@ inline FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>:: FEEvaluationAccess(const MatrixFree<1, Number, VectorizedArrayType> &data_in, const unsigned int dof_no, const unsigned int first_selected_component, - const unsigned int quad_no_in, + const unsigned int quad_no, const unsigned int fe_degree, const unsigned int n_q_points, const bool is_interior_face, @@ -7428,7 +6842,7 @@ inline FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>:: data_in, dof_no, first_selected_component, - quad_no_in, + quad_no, fe_degree, n_q_points, is_interior_face, @@ -8377,7 +7791,7 @@ FEEvaluationvalues_quad, this->gradients_quad, this->hessians_quad, - this->scratch_data); + this->scratch_data.begin()); else internal::FEEvaluationFactory::evaluate( n_components, @@ -8387,7 +7801,7 @@ FEEvaluationvalues_quad, this->gradients_quad, this->hessians_quad, - this->scratch_data); + this->scratch_data.begin()); if (hessians_on_general_cells) for (unsigned int i = 0; i < n_components * dim * n_q_points; ++i) @@ -8808,7 +8222,7 @@ FEEvaluationvalues_quad, this->gradients_quad, this->hessians_quad, - this->scratch_data, + this->scratch_data.begin(), false); else internal::FEEvaluationFactory::integrate( @@ -8819,7 +8233,7 @@ FEEvaluationvalues_quad, this->gradients_quad, this->hessians_quad, - this->scratch_data, + this->scratch_data.begin(), false); # ifdef DEBUG @@ -8888,7 +8302,7 @@ FEEvaluationvalues_quad, this->gradients_quad, this->hessians_quad, - this->scratch_data, + this->scratch_data.begin(), this->data, integration_flag) == false) { @@ -9256,7 +8670,7 @@ FEFaceEvaluationbegin_values(), this->begin_gradients(), this->begin_hessians(), - this->scratch_data, + this->scratch_data.begin(), evaluation_flag_actual & EvaluationFlags::values, evaluation_flag_actual & EvaluationFlags::gradients, evaluation_flag_actual & EvaluationFlags::hessians, @@ -9277,7 +8691,7 @@ FEFaceEvaluationbegin_values(), this->begin_gradients(), this->begin_hessians(), - this->scratch_data, + this->scratch_data.begin(), evaluation_flag_actual & EvaluationFlags::values, evaluation_flag_actual & EvaluationFlags::gradients, evaluation_flag_actual & EvaluationFlags::hessians, @@ -9293,7 +8707,7 @@ FEFaceEvaluationbegin_values(), this->begin_gradients(), this->begin_hessians(), - this->scratch_data, + this->scratch_data.begin(), evaluation_flag_actual & EvaluationFlags::values, evaluation_flag_actual & EvaluationFlags::gradients, evaluation_flag_actual & EvaluationFlags::hessians, @@ -9451,7 +8865,7 @@ FEFaceEvaluationbegin_values(), this->begin_gradients(), this->begin_hessians(), - this->scratch_data, + this->scratch_data.begin(), evaluation_flag_actual & EvaluationFlags::values, evaluation_flag_actual & EvaluationFlags::gradients, evaluation_flag_actual & EvaluationFlags::hessians, @@ -9467,7 +8881,7 @@ FEFaceEvaluationbegin_values(), this->begin_gradients(), this->begin_hessians(), - this->scratch_data, + this->scratch_data.begin(), evaluation_flag_actual & EvaluationFlags::values, evaluation_flag_actual & EvaluationFlags::gradients, evaluation_flag_actual & EvaluationFlags::hessians, @@ -9560,7 +8974,7 @@ FEFaceEvaluationbegin_values(), this->begin_gradients(), this->begin_hessians(), - this->scratch_data, + this->scratch_data.begin(), evaluation_flag & EvaluationFlags::values, evaluation_flag & EvaluationFlags::gradients, evaluation_flag & EvaluationFlags::hessians, @@ -9602,7 +9016,7 @@ FEFaceEvaluationbegin_values(), this->begin_gradients(), this->begin_hessians(), - this->scratch_data, + this->scratch_data.begin(), evaluation_flag & EvaluationFlags::values, evaluation_flag & EvaluationFlags::gradients, evaluation_flag & EvaluationFlags::hessians, @@ -9627,7 +9041,7 @@ FEFaceEvaluationbegin_values(), this->begin_gradients(), this->begin_hessians(), - this->scratch_data, + this->scratch_data.begin(), evaluation_flag & EvaluationFlags::values, evaluation_flag & EvaluationFlags::gradients, evaluation_flag & EvaluationFlags::hessians, @@ -9744,7 +9158,7 @@ FEFaceEvaluationbegin_values(), this->begin_gradients(), this->begin_hessians(), - this->scratch_data, + this->scratch_data.begin(), evaluation_flag & EvaluationFlags::values, evaluation_flag & EvaluationFlags::gradients, evaluation_flag & EvaluationFlags::hessians, @@ -9777,7 +9191,7 @@ FEFaceEvaluationbegin_values(), this->begin_gradients(), this->begin_hessians(), - this->scratch_data, + this->scratch_data.begin(), evaluation_flag & EvaluationFlags::values, evaluation_flag & EvaluationFlags::gradients, evaluation_flag & EvaluationFlags::hessians, diff --git a/include/deal.II/matrix_free/fe_evaluation_base_data.h b/include/deal.II/matrix_free/fe_evaluation_base_data.h new file mode 100644 index 0000000000..dcce525ab9 --- /dev/null +++ b/include/deal.II/matrix_free/fe_evaluation_base_data.h @@ -0,0 +1,763 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 - 2021 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +#ifndef dealii_matrix_free_fe_evaluation_base_data_h +#define dealii_matrix_free_fe_evaluation_base_data_h + + +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + + +DEAL_II_NAMESPACE_OPEN + + + +namespace internal +{ + DeclException0(ExcAccessToUninitializedField); + + DeclException1( + ExcMatrixFreeAccessToUninitializedMappingField, + std::string, + << "You are requesting information from an FEEvaluation/FEFaceEvaluation " + << "object for which this kind of information has not been computed. What " + << "information these objects compute is determined by the update_* flags you " + << "pass to MatrixFree::reinit() via MatrixFree::AdditionalData. Here, " + << "the operation you are attempting requires the <" << arg1 + << "> flag to be set, but it was apparently not specified " + << "upon initialization."); +} // namespace internal + +// forward declarations +template > +class FEEvaluation; + +template > +class FEEvaluationBase; + +namespace internal +{ + namespace MatrixFreeFunctions + { + template + class MappingDataOnTheFly; + } +} // namespace internal + + +/** + * This base class of the FEEvaluation and FEFaceEvaluation classes handles + * mapping-related information independent of the degrees of freedom and + * finite element in use. This class provides access functionality for user + * code but is otherwise invisible without any public constructor. The usage + * is through the class FEEvaluation instead. + * + * This class has four template arguments: + * + * @tparam dim Dimension in which this class is to be used + * + * @tparam Number Number format, usually @p double or @p float + * + * @tparam is_face Whether the class is used for a cell integrator (with + * quadrature dimension the same as the space dimension) or for a face + * integrator (with quadrature dimension one less) + * + * @tparam VectorizedArrayType Type of array to be woked on in a vectorized + * fashion, defaults to VectorizedArray + * + * @note Currently only VectorizedArray is supported as + * VectorizedArrayType. + * + * + * @ingroup matrixfree + */ +template +class FEEvaluationBaseData +{ + static_assert( + std::is_same::value, + "Type of Number and of VectorizedArrayType do not match."); + + using ShapeInfoType = + internal::MatrixFreeFunctions::ShapeInfo; + using MappingInfoStorageType = + internal::MatrixFreeFunctions::MappingInfoStorage<(is_face ? dim - 1 : dim), + dim, + Number, + VectorizedArrayType>; + using DoFInfo = internal::MatrixFreeFunctions::DoFInfo; + +public: + static constexpr unsigned int dimension = dim; + + /** + * Return the index offset within the geometry fields for the cell the @p + * reinit() function has been called for. This index can be used to access + * an index into a field that has the same compression behavior as the + * Jacobian of the geometry, e.g., to store an effective coefficient tensors + * that combines a coefficient with the geometry for lower memory transfer + * as the available data fields. + */ + unsigned int + get_mapping_data_index_offset() const; + + /** + * Return the type of the cell the @p reinit() function has been called for. + * Valid values are @p cartesian for Cartesian cells (which allows for + * considerable data compression), @p affine for cells with affine mappings, + * and @p general for general cells without any compressed storage applied. + */ + internal::MatrixFreeFunctions::GeometryType + get_cell_type() const; + + /** + * Return a reference to the ShapeInfo object currently in use. + */ + const internal::MatrixFreeFunctions::ShapeInfo & + get_shape_info() const; + + /** + * Return a reference to the DoFInfo object currently in use. + */ + const internal::MatrixFreeFunctions::DoFInfo & + get_dof_info() const; + + /** + * Return an ArrayView to internal memory for temporary use. Note that some + * of this memory is overwritten during evaluate() and integrate() calls so + * do not assume it to be stable over those calls. The maximum size you can + * write into is 3*dofs_per_cell+2*n_q_points. + */ + ArrayView + get_scratch_data() const; + + /** + * Return the determinant of the Jacobian from the unit to the real cell + * times the quadrature weight. + */ + VectorizedArrayType + JxW(const unsigned int q_point) 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 + * $J_{ij} = d x_i / d\hat x_j$. The $(i,j)$ entry of the returned tensor + * contains $d\hat x_j/dx_i$, i.e., columns refer to reference space + * coordinates and rows to real cell coordinates. Thus, the returned tensor + * represents a covariant transformation, which is used in the + * FEEvaluationBase::get_gradient() function to transform the unit cell + * gradients to gradients on the real cell by a multiplication $J^{-\mathrm + * T} \hat{\nabla} u_h$. + */ + Tensor<2, dim, VectorizedArrayType> + inverse_jacobian(const unsigned int q_point) const; + + /** + * Return the unit normal vector on a face. Note that both sides of a face + * use the same orientation of the normal vector: For the faces enumerated + * as `interior` in FaceToCellTopology and selected with the + * `is_interior_face=true` flag of the constructor, this corresponds to the + * outer normal vector, whereas for faces enumerated as `exterior` in + * FaceToCellTopology and selected with the `is_interior_face=false` flag of + * the constructor, the normal points into the element as a consequence of + * the single normal vector. + * + * @note Only implemented in case `is_face == true`. + */ + Tensor<1, dim, VectorizedArrayType> + get_normal_vector(const unsigned int q_point) const; + + /** + * Return the numbering of local degrees of freedom within the evaluation + * routines of FEEvaluation in terms of the standard numbering on finite + * elements. + */ + const std::vector & + get_internal_dof_numbering() const; + + /** + * Return the number of the quadrature formula of the present cell. + */ + unsigned int + get_quadrature_index() const; + + /** + * Return index of the current cell or face. + */ + unsigned int + get_current_cell_index() const; + + /** + * Return the active FE index for this class for efficient indexing in the + * hp-case. + */ + unsigned int + get_active_fe_index() const; + + /** + * Return the active quadrature index for this class for efficient indexing in + * the hp-case. + */ + unsigned int + get_active_quadrature_index() const; + +protected: + /** + * Constructor. Made protected to prevent users from directly using this + * class. Takes all data stored in MatrixFree. If applied to problems with + * more than one quadrature formula selected during construction of + * `matrix_free`, `quad_no` allows to select the appropriate formula. + */ + FEEvaluationBaseData(const ShapeInfoType & shape_info, + const DoFInfo & dof_info, + const MappingInfoStorageType &mapping_data, + const unsigned int quad_no, + const bool is_interior_face, + const unsigned int active_fe_index, + const unsigned int active_quad_index, + const unsigned int face_type); + + /** + * Constructor that comes with reduced functionality and works similar as + * FEValues. + */ + FEEvaluationBaseData( + const std::shared_ptr< + internal::MatrixFreeFunctions:: + MappingDataOnTheFly> &mapping_data); + + /** + * Copy constructor. + */ + FEEvaluationBaseData(const FEEvaluationBaseData &other); + + /** + * Copy assignment operator. If FEEvaluationBase was constructed from a + * mapping, fe, quadrature, and update flags, the underlying geometry + * evaluation based on FEValues will be deep-copied in order to allow for + * using in parallel with threads. + */ + FEEvaluationBaseData & + operator=(const FEEvaluationBaseData &other); + + /** + * A pointer to the unit cell shape data, i.e., values, gradients and + * Hessians in 1D at the quadrature points that constitute the tensor + * product. Also contained in matrix_info, but it simplifies code if we + * store a reference to it. + */ + const ShapeInfoType *data; + + /** + * A pointer to the underlying DoF indices and constraint description + * for the component specified at construction. Also contained in + * matrix_info, but it simplifies code if we store a reference to it. + */ + const DoFInfo *dof_info; + + /** + * A pointer to the underlying transformation data from unit to real cells + * for the given quadrature formula specified at construction. Also + * contained in matrix_info, but it simplifies code if we store a reference + * to it. + */ + const MappingInfoStorageType *mapping_data; + + /** + * The number of the quadrature formula of the present cell among all + * quadrature formulas available in the MatrixFree objects pass to derived + * classes. + */ + const unsigned int quad_no; + + /** + * The active FE index for this class for efficient indexing in the hp-case. + */ + const unsigned int active_fe_index; + + /** + * The active quadrature index for this class for efficient indexing in the + * hp-case. + */ + const unsigned int active_quad_index; + + /** + * A pointer to the underlying quadrature formula specified at construction. + * Also contained in mapping_data, but it simplifies code if we store a + * reference to it. + */ + const typename MappingInfoStorageType::QuadratureDescriptor *descriptor; + + /** + * The number of quadrature points in the current evaluation context. + */ + const unsigned int n_quadrature_points; + + /** + * A pointer to the Jacobian information of the present cell. Only set to a + * useful value if on a non-Cartesian cell. + */ + const Tensor<2, dim, VectorizedArrayType> *jacobian; + + /** + * A pointer to the Jacobian determinant of the present cell. If on a + * Cartesian cell or on a cell with constant Jacobian, this is just the + * Jacobian determinant, otherwise the Jacobian determinant times the + * quadrature weight. + */ + const VectorizedArrayType *J_value; + + /** + * A pointer to the normal vectors at faces. + */ + const Tensor<1, dim, VectorizedArrayType> *normal_vectors; + + /** + * A pointer to the normal vectors times the jacobian at faces. + */ + const Tensor<1, dim, VectorizedArrayType> *normal_x_jacobian; + + /** + * A pointer to the quadrature weights of the underlying quadrature formula. + */ + const Number *quadrature_weights; + + /** + * This is the user-visible part of FEEvaluationBase::scratch_data_array, + * only showing the part that can be consumed by various users. It is set + * during the allocation of the internal data structures in + * FEEvaluationBase. + */ + mutable ArrayView scratch_data; + + /** + * After a call to reinit(), stores the number of the cell we are currently + * working with. + */ + unsigned int cell; + + /** + * Flag holding information whether a face is an interior or exterior face + * according to the defined direction of the normal. For cells it defines if + * the dof values should be read from the actual cell corresponding to the + * interior face or the neighboring cell corresponding to the exterior face. + */ + bool is_interior_face; + + /** + * Stores the index an FEFaceEvaluation object is currently pointing into + * (interior face, exterior face, data associated with cell). + */ + internal::MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index; + + /** + * Stores the current number of a face within the given cell in case + * `is_face==true`, using values between `0` and `2*dim`. + */ + unsigned int face_no; + + /** + * Stores the orientation of the given face with respect to the standard + * orientation, 0 if in standard orientation. + */ + unsigned int face_orientation; + + /** + * Stores the subface index of the given face. Usually, this variable takes + * the value numbers::invalid_unsigned_int to indicate integration over the + * full face, but in case the current physical face has a neighbor that is + * more refined, it is a subface and must scale the entries in ShapeInfo + * appropriately. + */ + unsigned int subface_index; + + /** + * Stores the type of the cell we are currently working with after a call to + * reinit(). Valid values are @p cartesian, @p affine and @p general, which + * have different implications on how the Jacobian transformations are + * stored internally in MappingInfo. + */ + internal::MatrixFreeFunctions::GeometryType cell_type; + + /** + * Geometry data that can be generated FEValues on the fly with the + * respective constructor, as an alternative to the entry point with + * MatrixFree. + */ + std::shared_ptr> + mapped_geometry; + + // Make FEEvaluation and FEEvaluationBase objects friends for access to + // protected member mapped_geometry. + template + friend class FEEvaluationBase; + + template + friend class FEEvaluation; +}; + + + +/*----------------------- Inline functions ----------------------------------*/ + +#ifndef DOXYGEN + +template +inline FEEvaluationBaseData:: + FEEvaluationBaseData(const ShapeInfoType & shape_info, + const DoFInfo & dof_info, + const MappingInfoStorageType &mapping_data, + const unsigned int quad_no, + const bool is_interior_face, + const unsigned int active_fe_index, + const unsigned int active_quad_index, + const unsigned int face_type) + : data(&shape_info) + , dof_info(&dof_info) + , mapping_data(&mapping_data) + , quad_no(quad_no) + , active_fe_index(active_fe_index) + , active_quad_index(active_quad_index) + , descriptor( + &mapping_data.descriptor + [is_face ? + (active_quad_index * std::max(1, dim - 1) + + (face_type == numbers::invalid_unsigned_int ? 0 : face_type)) : + active_quad_index]) + , n_quadrature_points(descriptor->n_q_points) + , jacobian(nullptr) + , J_value(nullptr) + , normal_vectors(nullptr) + , normal_x_jacobian(nullptr) + , quadrature_weights(descriptor->quadrature_weights.begin()) + , cell(numbers::invalid_unsigned_int) + , is_interior_face(is_interior_face) + , dof_access_index( + is_face ? + (is_interior_face ? + internal::MatrixFreeFunctions::DoFInfo::dof_access_face_interior : + internal::MatrixFreeFunctions::DoFInfo::dof_access_face_exterior) : + internal::MatrixFreeFunctions::DoFInfo::dof_access_cell) + , face_no(0) + , face_orientation(0) + , subface_index(0) + , cell_type(internal::MatrixFreeFunctions::general) +{} + + + +template +inline FEEvaluationBaseData:: + FEEvaluationBaseData( + const std::shared_ptr< + internal::MatrixFreeFunctions:: + MappingDataOnTheFly> &mapped_geometry) + : data(nullptr) + , dof_info(nullptr) + , mapping_data(nullptr) + , quad_no(numbers::invalid_unsigned_int) + , active_fe_index(numbers::invalid_unsigned_int) + , active_quad_index(numbers::invalid_unsigned_int) + , descriptor(nullptr) + , n_quadrature_points( + mapped_geometry->get_data_storage().descriptor[0].n_q_points) + , jacobian(nullptr) + , J_value(nullptr) + , normal_vectors(nullptr) + , normal_x_jacobian(nullptr) + , quadrature_weights(nullptr) + , cell(0) + , cell_type(internal::MatrixFreeFunctions::general) + , is_interior_face(true) + , dof_access_index(internal::MatrixFreeFunctions::DoFInfo::dof_access_cell) + , mapped_geometry(mapped_geometry) +{ + 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(); +} + + + +template +inline FEEvaluationBaseData:: + FEEvaluationBaseData(const FEEvaluationBaseData &other) + : data(other.data) + , dof_info(other.dof_info) + , mapping_data(other.mapping_data) + , quad_no(other.quad_no) + , active_fe_index(other.active_fe_index) + , active_quad_index(other.active_quad_index) + , descriptor(other.descriptor) + , n_quadrature_points(other.n_quadrature_points) + , jacobian(nullptr) + , J_value(nullptr) + , normal_vectors(nullptr) + , normal_x_jacobian(nullptr) + , quadrature_weights(other.quadrature_weights) + , cell(numbers::invalid_unsigned_int) + , is_interior_face(other.is_interior_face) + , dof_access_index( + is_face ? + (is_interior_face ? + internal::MatrixFreeFunctions::DoFInfo::dof_access_face_interior : + internal::MatrixFreeFunctions::DoFInfo::dof_access_face_exterior) : + internal::MatrixFreeFunctions::DoFInfo::dof_access_cell) + , face_no(0) + , face_orientation(0) + , subface_index(0) + , cell_type(internal::MatrixFreeFunctions::general) +{} + + + +template +inline FEEvaluationBaseData & +FEEvaluationBaseData::operator=( + const FEEvaluationBaseData &other) +{ + AssertDimension(quad_no, other.quad_no); + AssertDimension(active_fe_index, other.active_fe_index); + AssertDimension(active_quad_index, other.active_quad_index); + AssertDimension(n_quadrature_points, descriptor->n_q_points); + + data = other.data; + dof_info = other.dof_info; + mapping_data = other.mapping_data; + descriptor = other.descriptor; + jacobian = nullptr; + J_value = nullptr; + normal_vectors = nullptr; + normal_x_jacobian = nullptr; + quadrature_weights = other.quadrature_weights; + cell = numbers::invalid_unsigned_int; + is_interior_face = other.is_interior_face; + dof_access_index = + is_face ? + (is_interior_face ? + internal::MatrixFreeFunctions::DoFInfo::dof_access_face_interior : + internal::MatrixFreeFunctions::DoFInfo::dof_access_face_exterior) : + internal::MatrixFreeFunctions::DoFInfo::dof_access_cell; + face_no = 0; + face_orientation = 0; + subface_index = 0; + cell_type = internal::MatrixFreeFunctions::general; + + return *this; +} + + + +template +inline unsigned int +FEEvaluationBaseData:: + get_mapping_data_index_offset() const +{ + Assert(mapping_data != nullptr, ExcInternalError()); + + if (dof_info == nullptr) + return 0; + else + { + AssertIndexRange(cell, mapping_data->data_index_offsets.size()); + return mapping_data->data_index_offsets[cell]; + } +} + + + +template +inline internal::MatrixFreeFunctions::GeometryType +FEEvaluationBaseData::get_cell_type() + const +{ + Assert(cell != numbers::invalid_unsigned_int, ExcNotInitialized()); + return cell_type; +} + + + +template +inline const internal::MatrixFreeFunctions::ShapeInfo & +FEEvaluationBaseData:: + get_shape_info() const +{ + Assert(data != nullptr, ExcInternalError()); + return *data; +} + + + +template +inline const internal::MatrixFreeFunctions::DoFInfo & +FEEvaluationBaseData::get_dof_info() + const +{ + Assert(dof_info != nullptr, + ExcMessage( + "FEEvaluation was not initialized with a MatrixFree object!")); + return *dof_info; +} + + + +template +inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, VectorizedArrayType> +FEEvaluationBaseData:: + get_normal_vector(const unsigned int q_point) const +{ + AssertIndexRange(q_point, n_quadrature_points); + Assert(normal_vectors != nullptr, + internal::ExcMatrixFreeAccessToUninitializedMappingField( + "update_normal_vectors")); + if (this->cell_type <= internal::MatrixFreeFunctions::flat_faces) + return normal_vectors[0]; + else + return normal_vectors[q_point]; +} + + + +template +inline DEAL_II_ALWAYS_INLINE VectorizedArrayType +FEEvaluationBaseData::JxW( + const unsigned int q_point) const +{ + AssertIndexRange(q_point, n_quadrature_points); + Assert(J_value != nullptr, + internal::ExcMatrixFreeAccessToUninitializedMappingField( + "update_values|update_gradients")); + if (this->cell_type <= internal::MatrixFreeFunctions::affine) + { + Assert(this->quadrature_weights != nullptr, ExcInternalError()); + return J_value[0] * this->quadrature_weights[q_point]; + } + else + return J_value[q_point]; +} + + + +template +inline DEAL_II_ALWAYS_INLINE Tensor<2, dim, VectorizedArrayType> +FEEvaluationBaseData:: + inverse_jacobian(const unsigned int q_point) const +{ + AssertIndexRange(q_point, n_quadrature_points); + Assert(this->jacobian != nullptr, + internal::ExcMatrixFreeAccessToUninitializedMappingField( + "update_gradients")); + if (this->cell_type <= internal::MatrixFreeFunctions::affine) + return jacobian[0]; + else + return jacobian[q_point]; +} + + + +template +inline const std::vector & +FEEvaluationBaseData:: + get_internal_dof_numbering() const +{ + return data->lexicographic_numbering; +} + + + +template +inline unsigned int +FEEvaluationBaseData:: + get_quadrature_index() const +{ + return this->quad_no; +} + + + +template +inline unsigned int +FEEvaluationBaseData:: + get_current_cell_index() const +{ + if (is_face && this->dof_access_index == + internal::MatrixFreeFunctions::DoFInfo::dof_access_cell) + return this->cell * GeometryInfo::faces_per_cell + this->face_no; + else + return this->cell; +} + + + +template +inline unsigned int +FEEvaluationBaseData:: + get_active_fe_index() const +{ + return active_fe_index; +} + + + +template +inline unsigned int +FEEvaluationBaseData:: + get_active_quadrature_index() const +{ + return active_quad_index; +} + + + +template +inline ArrayView +FEEvaluationBaseData:: + get_scratch_data() const +{ + return scratch_data; +} + + + +#endif // ifndef DOXYGEN + + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/include/deal.II/matrix_free/mapping_info.h b/include/deal.II/matrix_free/mapping_info.h index a97a97252b..4fc62e48be 100644 --- a/include/deal.II/matrix_free/mapping_info.h +++ b/include/deal.II/matrix_free/mapping_info.h @@ -34,6 +34,7 @@ #include #include +#include #include @@ -45,285 +46,6 @@ namespace internal { namespace MatrixFreeFunctions { - /** - * An enum to identify various types of cells and faces. The most general - * type is what we typically compute in the FEValues context but for many - * geometries we can save significant storage. - * - * @ingroup matrixfree - */ - enum GeometryType : unsigned char - { - /** - * The cell or face is Cartesian. - */ - cartesian = 0, - - /** - * The cell or face can be described with an affine mapping. - */ - affine = 1, - - /** - * The face is flat, i.e., the normal factor on a face is the same on - * all quadrature points. This type is not assigned for cells. - */ - flat_faces = 2, - - /** - * There is no special information available for compressing the - * representation of the object under consideration. - */ - general = 3 - }; - - - - /** - * Definition of a structure that stores all cached data related to the - * evaluated geometry from the mapping. In order to support hp-adaptivity - * and compressed storage (in particular for Jacobians, JxW values, and - * normals), storage length can be different for different rows. Thus, it - * allows to jump at the data of individual rows similar to compressed row - * storage in sparse matrices. We have two different start indices for - * fields with different sizes. The first category of offsets are the - * indices for Jacobians of the transformation from unit to real cell (we - * store the inverse Jacobian), second derivatives, JxW values, and normal - * vectors. We keep separate arrays for all these data structures because - * a user code might access only some of them. In such a case, one array - * will be gone through in a contiguous order with access to all entries, - * which makes it easy for the processor to prefetch data. Having all data - * in a single array would require some strides in the access pattern, - * which is much more complicated for the processor to predict (and indeed - * leads to prefetching of data that does not get used on Intel processors - * such as BroadwellEP). - * - * The second category of indices are the offsets for the quadrature - * points. Quadrature points can be compressed less than the other fields - * and thus need longer fields. Quadrature point indices are often used in - * other contexts such as evaluation of right hand sides. - * - * The third component is a descriptor of data from the unit cells, called - * QuadratureDescriptor, which contains the quadrature weights and - * permutations of how to go through quadrature points in case of face - * data. The latter comes in a vector for the support of hp-adaptivity, - * with several data fields for the individual quadrature formulas. - * - * @ingroup matrixfree - */ - template - struct MappingInfoStorage - { - static_assert( - std::is_same::value, - "Type of Number and of VectorizedArrayType do not match."); - - struct QuadratureDescriptor - { - /** - * Constructor. Does nothing. - */ - QuadratureDescriptor(); - - /** - * Set up the lengths in the various members of this struct. - */ - template - void - initialize(const Quadrature &quadrature, - const UpdateFlags update_flags_inner_faces = update_default); - - /** - * Set up the lengths in the various members of this struct. - */ - void - initialize(const Quadrature<1> &quadrature_1d, - const UpdateFlags update_flags_inner_faces = update_default); - - /** - * Returns the memory consumption in bytes. - */ - std::size_t - memory_consumption() const; - - /** - * Number of quadrature points applied on the given cell or face. - */ - unsigned int n_q_points; - - /** - * Original one-dimensional quadrature formula applied on the given - * cell or face. - */ - Quadrature<1> quadrature_1d; - - /** - * Quadrature formula applied on the given cell or face. - */ - Quadrature quadrature; - - /** - * Quadrature weights separated by dimension for use in specific - * situations. - */ - std::array, structdim> tensor_quadrature_weights; - - /** - * A cached vector of quadrature weights in the given number format - * (non-vectorized, as it is cheap to broadcast the value to all lanes - * when it is used in a vectorized context). - */ - AlignedVector quadrature_weights; - - /** - * For quadrature on faces, the evaluation of basis functions is not - * in the correct order if a face is not in the standard orientation - * to a given element. This data structure is used to re-order the - * data evaluated on quadrature points to represent the correct order. - */ - dealii::Table<2, unsigned int> face_orientations; - }; - - /** - * A class describing the layout of the sections in the @p data_storage - * field and also includes some data that depends on the number of - * quadrature points in the hp-context such as the inner quadrature - * formula and re-indexing for faces that are not in the standard - * orientation. - */ - std::vector descriptor; - - /** - * Collection of quadrature formulae applied on the given face. - * - * @note Only filled for faces, since faces might be quadrilateral or - * triangle shaped. - */ - std::vector> q_collection; - - /** - * Stores the index offset into the arrays @p jxw_values, @p jacobians, - * @p normal_vectors and the second derivatives. Note that affine cells - * have shorter fields of length 1, where the others have lengths equal - * to the number of quadrature points of the given cell. - */ - AlignedVector data_index_offsets; - - /** - * The storage of the Jacobian determinant (times the quadrature weight - * in case the transformation is non-affine) on quadrature - * points. - * - * Indexed by @p data_index_offsets. - */ - AlignedVector JxW_values; - - /** - * Stores the normal vectors. - * - * Indexed by @p data_index_offsets. - */ - AlignedVector> normal_vectors; - - /** - * The storage of covariant transformation on quadrature points, i.e., - * the inverse and transposed Jacobians of the transformation from the - * unit to the real cell. - * - * Indexed by @p data_index_offsets. - * - * Contains two fields for access from both sides for interior faces, - * but the default case (cell integrals or boundary integrals) only - * fills the zeroth component and ignores the first one. - */ - std::array>, 2> - jacobians; - - /** - * The storage of the gradients of the inverse Jacobian - * transformation. Because of symmetry, only the upper diagonal and - * diagonal part are needed. The first index runs through the - * derivatives, starting with the diagonal and then continuing row-wise, - * i.e., $\partial^2/\partial x_1 \partial x_2$ first, then - * $\partial^2/\partial x_1 \partial x_3$, and so on. The second index - * is the spatial coordinate. - * - * Indexed by @p data_index_offsets. - * - * Contains two fields for access from both sides for interior faces, - * but the default case (cell integrals or boundary integrals) only - * fills the zeroth component and ignores the first one. - */ - std::array< - AlignedVector>>, - 2> - jacobian_gradients; - - /** - * Stores the Jacobian transformations times the normal vector (this - * represents a shortcut that is accessed often and can thus get higher - * performance). - * - * Indexed by @p data_index_offsets. - */ - std::array>, 2> - normals_times_jacobians; - - /** - * Stores the index offset of a particular cell into the quadrature - * points array in real coordinates. Note that Cartesian cells have - * shorter fields (length is @p n_q_points_1d) than non-Cartesian cells - * (length is @p n_q_points) or faces. - */ - AlignedVector quadrature_point_offsets; - - /** - * Stores the quadrature points in real coordinates, including a - * compression scheme for Cartesian cells where we do not need to store - * the full data on all points. - * - * Indexed by @p quadrature_point_offsets. - */ - AlignedVector> quadrature_points; - - /** - * Clears all data fields except the descriptor vector. - */ - void - clear_data_fields(); - - /** - * Returns the quadrature index for a given number of quadrature - * points. If not in hp-mode or if the index is not found, this - * function always returns index 0. Hence, this function does not - * check whether the given degree is actually present. - */ - unsigned int - quad_index_from_n_q_points(const unsigned int n_q_points) const; - - /** - * Prints a detailed summary of memory consumption in the different - * structures of this class to the given output stream. - */ - template - void - print_memory_consumption(StreamType & out, - const TaskInfo &task_info) const; - - /** - * Returns the memory consumption in bytes. - */ - std::size_t - memory_consumption() const; - }; - - - /** * The class that stores all geometry-dependent data related with cell * interiors for use in the matrix-free class. @@ -555,12 +277,12 @@ namespace internal template struct MappingInfoCellsOrFaces { - static const MappingInfoStorage * + static const MappingInfoStorage & get(const MappingInfo &mapping_info, const unsigned int quad_no) { AssertIndexRange(quad_no, mapping_info.cell_data.size()); - return &mapping_info.cell_data[quad_no]; + return mapping_info.cell_data[quad_no]; } }; @@ -568,12 +290,12 @@ namespace internal struct MappingInfoCellsOrFaces { static const MappingInfoStorage - * + & get(const MappingInfo &mapping_info, const unsigned int quad_no) { AssertIndexRange(quad_no, mapping_info.face_data.size()); - return &mapping_info.face_data[quad_no]; + return mapping_info.face_data[quad_no]; } }; @@ -651,22 +373,6 @@ namespace internal /* ------------------- inline functions ----------------------------- */ - template - inline unsigned int - MappingInfoStorage:: - quad_index_from_n_q_points(const unsigned int n_q_points) const - { - for (unsigned int i = 0; i < descriptor.size(); ++i) - if (n_q_points == descriptor[i].n_q_points) - return i; - return 0; - } - - - template inline GeometryType MappingInfo::get_cell_type( diff --git a/include/deal.II/matrix_free/mapping_info.templates.h b/include/deal.II/matrix_free/mapping_info.templates.h index 2c3072008f..d5c35a4972 100644 --- a/include/deal.II/matrix_free/mapping_info.templates.h +++ b/include/deal.II/matrix_free/mapping_info.templates.h @@ -30,6 +30,7 @@ #include #include +#include #include DEAL_II_NAMESPACE_OPEN @@ -39,231 +40,6 @@ namespace internal { namespace MatrixFreeFunctions { - /* ------------------------ MappingInfoStorage implementation ---------- */ - - template - MappingInfoStorage:: - QuadratureDescriptor::QuadratureDescriptor() - : n_q_points(numbers::invalid_unsigned_int) - {} - - - - template - template - void - MappingInfoStorage:: - QuadratureDescriptor::initialize( - const Quadrature &quadrature, - const UpdateFlags update_flags_inner_faces) - { - Assert(structdim + 1 <= spacedim || - update_flags_inner_faces == update_default, - ExcMessage("Volume cells do not allow for setting inner faces")); - this->quadrature = quadrature; - n_q_points = quadrature.size(); - quadrature_weights.resize(n_q_points); - for (unsigned int i = 0; i < n_q_points; ++i) - quadrature_weights[i] = quadrature.weight(i); - - // note: quadrature_1d and tensor_quadrature_weights are not set up - - // TODO: set up face_orientations - (void)update_flags_inner_faces; - } - - - - template - void - MappingInfoStorage:: - QuadratureDescriptor::initialize( - const Quadrature<1> &quadrature_1d, - const UpdateFlags update_flags_inner_faces) - { - Assert(structdim + 1 <= spacedim || - update_flags_inner_faces == update_default, - ExcMessage("Volume cells do not allow for setting inner faces")); - this->quadrature_1d = quadrature_1d; - quadrature = Quadrature(quadrature_1d); - n_q_points = quadrature.size(); - quadrature_weights.resize(n_q_points); - for (unsigned int i = 0; i < n_q_points; ++i) - quadrature_weights[i] = quadrature.weight(i); - - for (int d = 0; d < structdim; ++d) - { - tensor_quadrature_weights[d].resize(quadrature_1d.size()); - for (unsigned int i = 0; i < quadrature_1d.size(); ++i) - tensor_quadrature_weights[d][i] = quadrature_1d.weight(i); - } - - // face orientation for faces in 3D - if (structdim == spacedim - 1 && spacedim == 3 && - update_flags_inner_faces != update_default) - { - const unsigned int n = quadrature_1d.size(); - face_orientations.reinit(8, n * n); - for (unsigned int j = 0, i = 0; j < n; ++j) - for (unsigned int k = 0; k < n; ++k, ++i) - { - // face_orientation=true, face_flip=false, face_rotation=false - face_orientations[0][i] = i; - // face_orientation=false, face_flip=false, face_rotation=false - face_orientations[1][i] = j + k * n; - // face_orientation=true, face_flip=true, face_rotation=false - face_orientations[2][i] = (n - 1 - k) + (n - 1 - j) * n; - // face_orientation=false, face_flip=true, face_rotation=false - face_orientations[3][i] = (n - 1 - j) + (n - 1 - k) * n; - // face_orientation=true, face_flip=false, face_rotation=true - face_orientations[4][i] = j + (n - 1 - k) * n; - // face_orientation=false, face_flip=false, face_rotation=true - face_orientations[5][i] = k + (n - 1 - j) * n; - // face_orientation=true, face_flip=true, face_rotation=true - face_orientations[6][i] = (n - 1 - j) + k * n; - // face_orientation=false, face_flip=true, face_rotation=true - face_orientations[7][i] = (n - 1 - k) + j * n; - } - } - } - - - - template - std::size_t - MappingInfoStorage:: - QuadratureDescriptor::memory_consumption() const - { - std::size_t memory = sizeof(this) + quadrature.memory_consumption() + - quadrature_weights.memory_consumption() + - face_orientations.memory_consumption(); - for (int d = 0; d < structdim; ++d) - memory += tensor_quadrature_weights[d].memory_consumption(); - return memory; - } - - - - template - void - MappingInfoStorage:: - clear_data_fields() - { - data_index_offsets.clear(); - JxW_values.clear(); - normal_vectors.clear(); - for (unsigned int i = 0; i < 2; ++i) - { - jacobians[i].clear(); - jacobian_gradients[i].clear(); - normals_times_jacobians[i].clear(); - } - quadrature_point_offsets.clear(); - quadrature_points.clear(); - } - - - - template - std::size_t - MappingInfoStorage:: - memory_consumption() const - { - return MemoryConsumption::memory_consumption(descriptor) + - MemoryConsumption::memory_consumption(data_index_offsets) + - MemoryConsumption::memory_consumption(JxW_values) + - MemoryConsumption::memory_consumption(normal_vectors) + - MemoryConsumption::memory_consumption(jacobians[0]) + - MemoryConsumption::memory_consumption(jacobians[1]) + - MemoryConsumption::memory_consumption(jacobian_gradients[0]) + - MemoryConsumption::memory_consumption(jacobian_gradients[1]) + - MemoryConsumption::memory_consumption(normals_times_jacobians[0]) + - MemoryConsumption::memory_consumption(normals_times_jacobians[1]) + - MemoryConsumption::memory_consumption(quadrature_point_offsets) + - MemoryConsumption::memory_consumption(quadrature_points); - } - - - - template - template - void - MappingInfoStorage:: - print_memory_consumption(StreamType &out, const TaskInfo &task_info) const - { - // print_memory_statistics involves global communication, so we can - // disable the check here only if no processor has any such data - const std::size_t size = - Utilities::MPI::sum(jacobians[0].size(), task_info.communicator); - if (size > 0) - { - out << " Memory JxW data: "; - task_info.print_memory_statistics( - out, - MemoryConsumption::memory_consumption(data_index_offsets) + - MemoryConsumption::memory_consumption(JxW_values)); - out << " Memory Jacobian data: "; - task_info.print_memory_statistics( - out, - MemoryConsumption::memory_consumption(jacobians[0]) + - MemoryConsumption::memory_consumption(jacobians[1])); - out << " Memory second derivative data: "; - task_info.print_memory_statistics( - out, - MemoryConsumption::memory_consumption(jacobian_gradients[0]) + - MemoryConsumption::memory_consumption(jacobian_gradients[1])); - } - const std::size_t normal_size = - Utilities::MPI::sum(normal_vectors.size(), task_info.communicator); - if (normal_size > 0) - { - out << " Memory normal vectors data: "; - task_info.print_memory_statistics( - out, - MemoryConsumption::memory_consumption(normal_vectors) + - MemoryConsumption::memory_consumption( - normals_times_jacobians[0]) + - MemoryConsumption::memory_consumption( - normals_times_jacobians[1])); - } - - const std::size_t quad_size = - Utilities::MPI::sum(quadrature_points.size(), task_info.communicator); - if (quad_size > 0) - { - out << " Memory quadrature points: "; - task_info.print_memory_statistics( - out, - MemoryConsumption::memory_consumption(quadrature_point_offsets) + - MemoryConsumption::memory_consumption(quadrature_points)); - } - } - - - - /* ------------------------ MappingInfo implementation ----------------- */ - template void MappingInfo::clear() diff --git a/include/deal.II/matrix_free/mapping_info_storage.h b/include/deal.II/matrix_free/mapping_info_storage.h new file mode 100644 index 0000000000..d121e96a06 --- /dev/null +++ b/include/deal.II/matrix_free/mapping_info_storage.h @@ -0,0 +1,346 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2011 - 2021 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +#ifndef dealii_matrix_free_mapping_info_storage_h +#define dealii_matrix_free_mapping_info_storage_h + + +#include + +#include +#include +#include +#include +#include + +#include + +#include + +#include + + +DEAL_II_NAMESPACE_OPEN + + +namespace internal +{ + namespace MatrixFreeFunctions + { + // forward declaration + struct TaskInfo; + + /** + * An enum to identify various types of cells and faces. The most general + * type is what we typically compute in the FEValues context but for many + * geometries we can save significant storage. + * + * @ingroup matrixfree + */ + enum GeometryType : unsigned char + { + /** + * The cell or face is Cartesian. + */ + cartesian = 0, + + /** + * The cell or face can be described with an affine mapping. + */ + affine = 1, + + /** + * The face is flat, i.e., the normal factor on a face is the same on + * all quadrature points. This type is not assigned for cells. + */ + flat_faces = 2, + + /** + * There is no special information available for compressing the + * representation of the object under consideration. + */ + general = 3 + }; + + + + /** + * Definition of a structure that stores all cached data related to the + * evaluated geometry from the mapping. In order to support hp-adaptivity + * and compressed storage (in particular for Jacobians, JxW values, and + * normals), storage length can be different for different rows. Thus, it + * allows to jump at the data of individual rows similar to compressed row + * storage in sparse matrices. We have two different start indices for + * fields with different sizes. The first category of offsets are the + * indices for Jacobians of the transformation from unit to real cell (we + * store the inverse Jacobian), second derivatives, JxW values, and normal + * vectors. We keep separate arrays for all these data structures because + * a user code might access only some of them. In such a case, one array + * will be gone through in a contiguous order with access to all entries, + * which makes it easy for the processor to prefetch data. Having all data + * in a single array would require some strides in the access pattern, + * which is much more complicated for the processor to predict (and indeed + * leads to prefetching of data that does not get used on Intel processors + * such as BroadwellEP). + * + * The second category of indices are the offsets for the quadrature + * points. Quadrature points can be compressed less than the other fields + * and thus need longer fields. Quadrature point indices are often used in + * other contexts such as evaluation of right hand sides. + * + * The third component is a descriptor of data from the unit cells, called + * QuadratureDescriptor, which contains the quadrature weights and + * permutations of how to go through quadrature points in case of face + * data. The latter comes in a vector for the support of hp-adaptivity, + * with several data fields for the individual quadrature formulas. + * + * @ingroup matrixfree + */ + template + struct MappingInfoStorage + { + static_assert( + std::is_same::value, + "Type of Number and of VectorizedArrayType do not match."); + + struct QuadratureDescriptor + { + /** + * Constructor. Does nothing. + */ + QuadratureDescriptor(); + + /** + * Set up the lengths in the various members of this struct. + */ + template + void + initialize(const Quadrature &quadrature, + const UpdateFlags update_flags_inner_faces = update_default); + + /** + * Set up the lengths in the various members of this struct. + */ + void + initialize(const Quadrature<1> &quadrature_1d, + const UpdateFlags update_flags_inner_faces = update_default); + + /** + * Returns the memory consumption in bytes. + */ + std::size_t + memory_consumption() const; + + /** + * Number of quadrature points applied on the given cell or face. + */ + unsigned int n_q_points; + + /** + * Original one-dimensional quadrature formula applied on the given + * cell or face. + */ + Quadrature<1> quadrature_1d; + + /** + * Quadrature formula applied on the given cell or face. + */ + Quadrature quadrature; + + /** + * Quadrature weights separated by dimension for use in specific + * situations. + */ + std::array, structdim> tensor_quadrature_weights; + + /** + * A cached vector of quadrature weights in the given number format + * (non-vectorized, as it is cheap to broadcast the value to all lanes + * when it is used in a vectorized context). + */ + AlignedVector quadrature_weights; + + /** + * For quadrature on faces, the evaluation of basis functions is not + * in the correct order if a face is not in the standard orientation + * to a given element. This data structure is used to re-order the + * data evaluated on quadrature points to represent the correct order. + */ + dealii::Table<2, unsigned int> face_orientations; + }; + + /** + * A class describing the layout of the sections in the @p data_storage + * field and also includes some data that depends on the number of + * quadrature points in the hp-context such as the inner quadrature + * formula and re-indexing for faces that are not in the standard + * orientation. + */ + std::vector descriptor; + + /** + * Collection of quadrature formulae applied on the given face. + * + * @note Only filled for faces, since faces might be quadrilateral or + * triangle shaped. + */ + std::vector> q_collection; + + /** + * Stores the index offset into the arrays @p jxw_values, @p jacobians, + * @p normal_vectors and the second derivatives. Note that affine cells + * have shorter fields of length 1, where the others have lengths equal + * to the number of quadrature points of the given cell. + */ + AlignedVector data_index_offsets; + + /** + * The storage of the Jacobian determinant (times the quadrature weight + * in case the transformation is non-affine) on quadrature + * points. + * + * Indexed by @p data_index_offsets. + */ + AlignedVector JxW_values; + + /** + * Stores the normal vectors. + * + * Indexed by @p data_index_offsets. + */ + AlignedVector> normal_vectors; + + /** + * The storage of covariant transformation on quadrature points, i.e., + * the inverse and transposed Jacobians of the transformation from the + * unit to the real cell. + * + * Indexed by @p data_index_offsets. + * + * Contains two fields for access from both sides for interior faces, + * but the default case (cell integrals or boundary integrals) only + * fills the zeroth component and ignores the first one. + */ + std::array>, 2> + jacobians; + + /** + * The storage of the gradients of the inverse Jacobian + * transformation. Because of symmetry, only the upper diagonal and + * diagonal part are needed. The first index runs through the + * derivatives, starting with the diagonal and then continuing row-wise, + * i.e., $\partial^2/\partial x_1 \partial x_2$ first, then + * $\partial^2/\partial x_1 \partial x_3$, and so on. The second index + * is the spatial coordinate. + * + * Indexed by @p data_index_offsets. + * + * Contains two fields for access from both sides for interior faces, + * but the default case (cell integrals or boundary integrals) only + * fills the zeroth component and ignores the first one. + */ + std::array< + AlignedVector>>, + 2> + jacobian_gradients; + + /** + * Stores the Jacobian transformations times the normal vector (this + * represents a shortcut that is accessed often and can thus get higher + * performance). + * + * Indexed by @p data_index_offsets. + */ + std::array>, 2> + normals_times_jacobians; + + /** + * Stores the index offset of a particular cell into the quadrature + * points array in real coordinates. Note that Cartesian cells have + * shorter fields (length is @p n_q_points_1d) than non-Cartesian cells + * (length is @p n_q_points) or faces. + */ + AlignedVector quadrature_point_offsets; + + /** + * Stores the quadrature points in real coordinates, including a + * compression scheme for Cartesian cells where we do not need to store + * the full data on all points. + * + * Indexed by @p quadrature_point_offsets. + */ + AlignedVector> quadrature_points; + + /** + * Clears all data fields except the descriptor vector. + */ + void + clear_data_fields(); + + /** + * Returns the quadrature index for a given number of quadrature + * points. If not in hp-mode or if the index is not found, this + * function always returns index 0. Hence, this function does not + * check whether the given degree is actually present. + */ + unsigned int + quad_index_from_n_q_points(const unsigned int n_q_points) const; + + /** + * Prints a detailed summary of memory consumption in the different + * structures of this class to the given output stream. + */ + template + void + print_memory_consumption(StreamType & out, + const TaskInfo &task_info) const; + + /** + * Returns the memory consumption in bytes. + */ + std::size_t + memory_consumption() const; + }; + + + + /* ------------------- inline functions ----------------------------- */ + + template + inline unsigned int + MappingInfoStorage:: + quad_index_from_n_q_points(const unsigned int n_q_points) const + { + for (unsigned int i = 0; i < descriptor.size(); ++i) + if (n_q_points == descriptor[i].n_q_points) + return i; + return 0; + } + + } // end of namespace MatrixFreeFunctions +} // end of namespace internal + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/include/deal.II/matrix_free/mapping_info_storage.templates.h b/include/deal.II/matrix_free/mapping_info_storage.templates.h new file mode 100644 index 0000000000..69509d1024 --- /dev/null +++ b/include/deal.II/matrix_free/mapping_info_storage.templates.h @@ -0,0 +1,267 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2011 - 2021 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#ifndef dealii_matrix_free_mapping_info_storage_templates_h +#define dealii_matrix_free_mapping_info_storage_templates_h + +#include + +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + + +namespace internal +{ + namespace MatrixFreeFunctions + { + template + MappingInfoStorage:: + QuadratureDescriptor::QuadratureDescriptor() + : n_q_points(numbers::invalid_unsigned_int) + {} + + + + template + template + void + MappingInfoStorage:: + QuadratureDescriptor::initialize( + const Quadrature &quadrature, + const UpdateFlags update_flags_inner_faces) + { + Assert(structdim + 1 <= spacedim || + update_flags_inner_faces == update_default, + ExcMessage("Volume cells do not allow for setting inner faces")); + this->quadrature = quadrature; + n_q_points = quadrature.size(); + quadrature_weights.resize(n_q_points); + for (unsigned int i = 0; i < n_q_points; ++i) + quadrature_weights[i] = quadrature.weight(i); + + // note: quadrature_1d and tensor_quadrature_weights are not set up + + // TODO: set up face_orientations + (void)update_flags_inner_faces; + } + + + + template + void + MappingInfoStorage:: + QuadratureDescriptor::initialize( + const Quadrature<1> &quadrature_1d, + const UpdateFlags update_flags_inner_faces) + { + Assert(structdim + 1 <= spacedim || + update_flags_inner_faces == update_default, + ExcMessage("Volume cells do not allow for setting inner faces")); + this->quadrature_1d = quadrature_1d; + quadrature = Quadrature(quadrature_1d); + n_q_points = quadrature.size(); + quadrature_weights.resize(n_q_points); + for (unsigned int i = 0; i < n_q_points; ++i) + quadrature_weights[i] = quadrature.weight(i); + + for (int d = 0; d < structdim; ++d) + { + tensor_quadrature_weights[d].resize(quadrature_1d.size()); + for (unsigned int i = 0; i < quadrature_1d.size(); ++i) + tensor_quadrature_weights[d][i] = quadrature_1d.weight(i); + } + + // face orientation for faces in 3D + if (structdim == spacedim - 1 && spacedim == 3 && + update_flags_inner_faces != update_default) + { + const unsigned int n = quadrature_1d.size(); + face_orientations.reinit(8, n * n); + for (unsigned int j = 0, i = 0; j < n; ++j) + for (unsigned int k = 0; k < n; ++k, ++i) + { + // face_orientation=true, face_flip=false, face_rotation=false + face_orientations[0][i] = i; + // face_orientation=false, face_flip=false, face_rotation=false + face_orientations[1][i] = j + k * n; + // face_orientation=true, face_flip=true, face_rotation=false + face_orientations[2][i] = (n - 1 - k) + (n - 1 - j) * n; + // face_orientation=false, face_flip=true, face_rotation=false + face_orientations[3][i] = (n - 1 - j) + (n - 1 - k) * n; + // face_orientation=true, face_flip=false, face_rotation=true + face_orientations[4][i] = j + (n - 1 - k) * n; + // face_orientation=false, face_flip=false, face_rotation=true + face_orientations[5][i] = k + (n - 1 - j) * n; + // face_orientation=true, face_flip=true, face_rotation=true + face_orientations[6][i] = (n - 1 - j) + k * n; + // face_orientation=false, face_flip=true, face_rotation=true + face_orientations[7][i] = (n - 1 - k) + j * n; + } + } + } + + + + template + std::size_t + MappingInfoStorage:: + QuadratureDescriptor::memory_consumption() const + { + std::size_t memory = sizeof(this) + quadrature.memory_consumption() + + quadrature_weights.memory_consumption() + + face_orientations.memory_consumption(); + for (int d = 0; d < structdim; ++d) + memory += tensor_quadrature_weights[d].memory_consumption(); + return memory; + } + + + + template + void + MappingInfoStorage:: + clear_data_fields() + { + data_index_offsets.clear(); + JxW_values.clear(); + normal_vectors.clear(); + for (unsigned int i = 0; i < 2; ++i) + { + jacobians[i].clear(); + jacobian_gradients[i].clear(); + normals_times_jacobians[i].clear(); + } + quadrature_point_offsets.clear(); + quadrature_points.clear(); + } + + + + template + std::size_t + MappingInfoStorage:: + memory_consumption() const + { + return MemoryConsumption::memory_consumption(descriptor) + + MemoryConsumption::memory_consumption(data_index_offsets) + + MemoryConsumption::memory_consumption(JxW_values) + + MemoryConsumption::memory_consumption(normal_vectors) + + MemoryConsumption::memory_consumption(jacobians[0]) + + MemoryConsumption::memory_consumption(jacobians[1]) + + MemoryConsumption::memory_consumption(jacobian_gradients[0]) + + MemoryConsumption::memory_consumption(jacobian_gradients[1]) + + MemoryConsumption::memory_consumption(normals_times_jacobians[0]) + + MemoryConsumption::memory_consumption(normals_times_jacobians[1]) + + MemoryConsumption::memory_consumption(quadrature_point_offsets) + + MemoryConsumption::memory_consumption(quadrature_points); + } + + + + template + template + void + MappingInfoStorage:: + print_memory_consumption(StreamType &out, const TaskInfo &task_info) const + { + // print_memory_statistics involves global communication, so we can + // disable the check here only if no processor has any such data + const std::size_t size = + Utilities::MPI::sum(jacobians[0].size(), task_info.communicator); + if (size > 0) + { + out << " Memory JxW data: "; + task_info.print_memory_statistics( + out, + MemoryConsumption::memory_consumption(data_index_offsets) + + MemoryConsumption::memory_consumption(JxW_values)); + out << " Memory Jacobian data: "; + task_info.print_memory_statistics( + out, + MemoryConsumption::memory_consumption(jacobians[0]) + + MemoryConsumption::memory_consumption(jacobians[1])); + out << " Memory second derivative data: "; + task_info.print_memory_statistics( + out, + MemoryConsumption::memory_consumption(jacobian_gradients[0]) + + MemoryConsumption::memory_consumption(jacobian_gradients[1])); + } + const std::size_t normal_size = + Utilities::MPI::sum(normal_vectors.size(), task_info.communicator); + if (normal_size > 0) + { + out << " Memory normal vectors data: "; + task_info.print_memory_statistics( + out, + MemoryConsumption::memory_consumption(normal_vectors) + + MemoryConsumption::memory_consumption( + normals_times_jacobians[0]) + + MemoryConsumption::memory_consumption( + normals_times_jacobians[1])); + } + + const std::size_t quad_size = + Utilities::MPI::sum(quadrature_points.size(), task_info.communicator); + if (quad_size > 0) + { + out << " Memory quadrature points: "; + task_info.print_memory_statistics( + out, + MemoryConsumption::memory_consumption(quadrature_point_offsets) + + MemoryConsumption::memory_consumption(quadrature_points)); + } + } + + } // namespace MatrixFreeFunctions +} // end of namespace internal + +DEAL_II_NAMESPACE_CLOSE + +#endif -- 2.39.5