From 634d503791869fb2cd8ae89ce6d37a4aee7b59eb Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Tue, 30 Nov 2021 10:43:11 +0100 Subject: [PATCH] Move implementation for ShapeInfo to templates.h file --- include/deal.II/matrix_free/fe_evaluation.h | 1 - include/deal.II/matrix_free/shape_info.h | 55 ++++---- .../matrix_free/shape_info.templates.h | 122 +++++++++++------- source/matrix_free/shape_info.inst.in | 50 +++++-- 4 files changed, 130 insertions(+), 98 deletions(-) diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index 75e7f55842..4fc5c6cc2a 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -232,7 +232,6 @@ public: 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 diff --git a/include/deal.II/matrix_free/shape_info.h b/include/deal.II/matrix_free/shape_info.h index 59fe6a9bc1..241cbe517f 100644 --- a/include/deal.II/matrix_free/shape_info.h +++ b/include/deal.II/matrix_free/shape_info.h @@ -22,14 +22,19 @@ #include #include -#include +#include +#include #include -#include - DEAL_II_NAMESPACE_OPEN + +// forward declaration +template +class FiniteElement; + + namespace internal { namespace MatrixFreeFunctions @@ -132,7 +137,7 @@ namespace internal ElementType element_type; /** - * Stores the shape values of the 1D finite element evaluated on all 1D + * Stores the shape values of the 1D finite element evaluated at all 1D * quadrature points. The length of * this array is n_dofs_1d * n_q_points_1d and quadrature * points are the index running fastest. @@ -140,7 +145,7 @@ namespace internal AlignedVector shape_values; /** - * Stores the shape gradients of the 1D finite element evaluated on all + * Stores the shape gradients of the 1D finite element evaluated at all * 1D quadrature points. The length of * this array is n_dofs_1d * n_q_points_1d and quadrature * points are the index running fastest. @@ -148,7 +153,7 @@ namespace internal AlignedVector shape_gradients; /** - * Stores the shape Hessians of the 1D finite element evaluated on all + * Stores the shape Hessians of the 1D finite element evaluated at all * 1D quadrature points. The length of * this array is n_dofs_1d * n_q_points_1d and quadrature * points are the index running fastest. @@ -300,16 +305,16 @@ namespace internal bool nodal_at_cell_boundaries; /** - * Stores the shape values of the finite element evaluated on all + * Stores the shape values of the finite element evaluated at all * quadrature points for all faces and orientations (no tensor-product * structure exploited). */ Table<3, Number> shape_values_face; /** - * Stores the shape gradients of the finite element evaluated on all + * Stores the shape gradients of the finite element evaluated at all * quadrature points for all faces, orientations, and directions - * (no tensor-product structure exploited). + * (no tensor-product structure exploited). */ Table<4, Number> shape_gradients_face; }; @@ -344,10 +349,10 @@ namespace internal /** * Constructor that initializes the data fields using the reinit method. */ - template - ShapeInfo(const Quadrature & quad, - const FiniteElement &fe, - const unsigned int base_element = 0); + template + ShapeInfo(const Quadrature & quad, + const FiniteElement &fe, + const unsigned int base_element = 0); /** * Initializes the data fields. Takes a one-dimensional quadrature @@ -357,11 +362,11 @@ namespace internal * dimensional element by a tensor product and that the zeroth shape * function in zero evaluates to one. */ - template + template void - reinit(const Quadrature & quad, - const FiniteElement &fe_dim, - const unsigned int base_element = 0); + reinit(const Quadrature & quad, + const FiniteElement &fe_dim, + const unsigned int base_element = 0); /** * Return which kinds of elements are supported by MatrixFree. @@ -556,22 +561,6 @@ namespace internal // ------------------------------------------ inline functions - template - template - inline ShapeInfo::ShapeInfo(const Quadrature & quad, - const FiniteElement &fe_in, - const unsigned int base_element_number) - : element_type(tensor_general) - , n_dimensions(0) - , n_components(0) - , n_q_points(0) - , dofs_per_component_on_cell(0) - , n_q_points_face(0) - , dofs_per_component_on_face(0) - { - reinit(quad, fe_in, base_element_number); - } - template inline const UnivariateShapeData & ShapeInfo::get_shape_data(const unsigned int dimension, diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h index c8a482d752..b982760224 100644 --- a/include/deal.II/matrix_free/shape_info.templates.h +++ b/include/deal.II/matrix_free/shape_info.templates.h @@ -26,6 +26,7 @@ #include #include +#include #include #include #include @@ -62,8 +63,6 @@ namespace internal - // ----------------- actual ShapeInfo functions -------------------- - template Number get_first_array_element(const Number a) @@ -196,6 +195,8 @@ namespace internal + // ----------------- actual ShapeInfo implementation -------------------- + template ShapeInfo::ShapeInfo() : element_type(tensor_general) @@ -210,62 +211,33 @@ namespace internal template - template - bool - ShapeInfo::is_supported(const FiniteElement &fe) + template + inline ShapeInfo::ShapeInfo( + const Quadrature & quad, + const FiniteElement &fe_in, + const unsigned int base_element_number) + : element_type(tensor_general) + , n_dimensions(0) + , n_components(0) + , n_q_points(0) + , dofs_per_component_on_cell(0) + , n_q_points_face(0) + , dofs_per_component_on_face(0) { - if (dim != spacedim) - return false; - - for (unsigned int base = 0; base < fe.n_base_elements(); ++base) - { - const FiniteElement *fe_ptr = &(fe.base_element(base)); - if (fe_ptr->n_components() != 1) - return false; - - // then check if the base element is supported or not - if (dynamic_cast *>(fe_ptr) != nullptr) - { - const FE_Poly *fe_poly_ptr = - dynamic_cast *>(fe_ptr); - // Simplices are a special case since the polynomial family is not - // indicative of their support - if (dynamic_cast *>(fe_poly_ptr) || - dynamic_cast *>(fe_poly_ptr) || - dynamic_cast *>(fe_poly_ptr) || - dynamic_cast *>(fe_poly_ptr)) - return true; - - if (dynamic_cast *>( - &fe_poly_ptr->get_poly_space()) == nullptr && - dynamic_cast> *>( - &fe_poly_ptr->get_poly_space()) == nullptr && - dynamic_cast *>(fe_ptr) == - nullptr && - dynamic_cast *>(fe_ptr) == - nullptr) - return false; - } - else - return false; - } - - // if we arrived here, all base elements were supported so we can - // support the present element - return true; + reinit(quad, fe_in, base_element_number); } template - template + template void - ShapeInfo::reinit(const Quadrature & quad_in, - const FiniteElement &fe_in, - const unsigned int base_element_number) + ShapeInfo::reinit(const Quadrature & quad_in, + const FiniteElement &fe_in, + const unsigned int base_element_number) { + static_assert(dim == spacedim, + "Currently, only the case dim=spacedim is implemented"); if (quad_in.is_tensor_product() == false || dynamic_cast *>( &fe_in.base_element(base_element_number)) || @@ -1135,6 +1107,56 @@ namespace internal + template + template + bool + ShapeInfo::is_supported(const FiniteElement &fe) + { + if (dim != spacedim) + return false; + + for (unsigned int base = 0; base < fe.n_base_elements(); ++base) + { + const FiniteElement *fe_ptr = &(fe.base_element(base)); + if (fe_ptr->n_components() != 1) + return false; + + // then check if the base element is supported or not + if (dynamic_cast *>(fe_ptr) != nullptr) + { + const FE_Poly *fe_poly_ptr = + dynamic_cast *>(fe_ptr); + // Simplices are a special case since the polynomial family is not + // indicative of their support + if (dynamic_cast *>(fe_poly_ptr) || + dynamic_cast *>(fe_poly_ptr) || + dynamic_cast *>(fe_poly_ptr) || + dynamic_cast *>(fe_poly_ptr)) + return true; + + if (dynamic_cast *>( + &fe_poly_ptr->get_poly_space()) == nullptr && + dynamic_cast> *>( + &fe_poly_ptr->get_poly_space()) == nullptr && + dynamic_cast *>(fe_ptr) == + nullptr && + dynamic_cast *>(fe_ptr) == + nullptr) + return false; + } + else + return false; + } + + // if we arrived here, all base elements were supported so we can + // support the present element + return true; + } + + + template std::size_t ShapeInfo::memory_consumption() const diff --git a/source/matrix_free/shape_info.inst.in b/source/matrix_free/shape_info.inst.in index 30be5ee7aa..fe5de6c930 100644 --- a/source/matrix_free/shape_info.inst.in +++ b/source/matrix_free/shape_info.inst.in @@ -16,18 +16,28 @@ for (deal_II_dimension : DIMENSIONS; deal_II_scalar : REAL_SCALARS) { - template void internal::MatrixFreeFunctions::ShapeInfo:: - reinit( - const Quadrature &, - const FiniteElement &, - const unsigned int); + template internal::MatrixFreeFunctions::ShapeInfo:: + ShapeInfo(const Quadrature &, + const FiniteElement &, + const unsigned int); + + template void + internal::MatrixFreeFunctions::ShapeInfo::reinit( + const Quadrature &, + const FiniteElement &, + const unsigned int); #if deal_II_dimension > 1 - template void internal::MatrixFreeFunctions::ShapeInfo:: - reinit( - const Quadrature<1> &, - const FiniteElement &, - const unsigned int); + template internal::MatrixFreeFunctions::ShapeInfo:: + ShapeInfo(const Quadrature<1> &, + const FiniteElement &, + const unsigned int); + + template void + internal::MatrixFreeFunctions::ShapeInfo::reinit( + const Quadrature<1> &, + const FiniteElement &, + const unsigned int); #endif } @@ -45,18 +55,30 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) for (deal_II_dimension : DIMENSIONS; deal_II_scalar_vectorized : REAL_SCALARS_VECTORIZED) { - template void internal::MatrixFreeFunctions:: - ShapeInfo::reinit( + template internal::MatrixFreeFunctions:: + ShapeInfo::ShapeInfo( const Quadrature &, const FiniteElement &, const unsigned int); + template void + internal::MatrixFreeFunctions::ShapeInfo::reinit( + const Quadrature &, + const FiniteElement &, + const unsigned int); + #if deal_II_dimension > 1 - template void internal::MatrixFreeFunctions:: - ShapeInfo::reinit( + template internal::MatrixFreeFunctions:: + ShapeInfo::ShapeInfo( const Quadrature<1> &, const FiniteElement &, const unsigned int); + + template void + internal::MatrixFreeFunctions::ShapeInfo::reinit( + const Quadrature<1> &, + const FiniteElement &, + const unsigned int); #endif } -- 2.39.5