From 1b7744c2a161d8bd5af3c4e5465e39b81ff3063f Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Thu, 22 Oct 2020 18:31:13 +0200 Subject: [PATCH] Improve MatrixFree::is_supported() --- .../matrix_free/matrix_free.templates.h | 33 +------------ include/deal.II/matrix_free/shape_info.h | 7 +++ .../matrix_free/shape_info.templates.h | 46 +++++++++++++++++++ .../numerics/vector_tools_project.templates.h | 3 +- source/matrix_free/shape_info.inst.in | 9 ++++ 5 files changed, 65 insertions(+), 33 deletions(-) diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h index b6152e0076..64c08905f0 100644 --- a/include/deal.II/matrix_free/matrix_free.templates.h +++ b/include/deal.II/matrix_free/matrix_free.templates.h @@ -464,38 +464,7 @@ bool MatrixFree::is_supported( const FiniteElement &fe) { - if (dim != spacedim) - return false; - - // first check for degree, number of base_elemnt and number of its components - if (fe.degree == 0 || fe.n_base_elements() != 1) - return false; - - const FiniteElement *fe_ptr = &(fe.base_element(0)); - if (fe_ptr->n_components() != 1) - return false; - - // then check of the base element is supported - if (dynamic_cast *>(fe_ptr) != nullptr) - { - const FE_Poly *fe_poly_ptr = - dynamic_cast *>(fe_ptr); - if (dynamic_cast *>( - &fe_poly_ptr->get_poly_space()) != nullptr) - return true; - if (dynamic_cast> *>( - &fe_poly_ptr->get_poly_space()) != nullptr) - return true; - } - if (dynamic_cast *>(fe_ptr) != nullptr) - return true; - if (dynamic_cast *>(fe_ptr) != nullptr) - return true; - - // if the base element is not in the above list it is not supported - return false; + return internal::MatrixFreeFunctions::ShapeInfo::is_supported(fe); } diff --git a/include/deal.II/matrix_free/shape_info.h b/include/deal.II/matrix_free/shape_info.h index 0637a5e775..64daf31106 100644 --- a/include/deal.II/matrix_free/shape_info.h +++ b/include/deal.II/matrix_free/shape_info.h @@ -331,6 +331,13 @@ namespace internal const FiniteElement &fe_dim, const unsigned int base_element = 0); + /** + * Return which kinds of elements are supported by MatrixFree. + */ + template + static bool + is_supported(const FiniteElement &fe); + /** * Return data of univariate shape functions which defines the * dimension @p dimension of tensor product shape functions diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h index 41a90beb54..f41881d869 100644 --- a/include/deal.II/matrix_free/shape_info.templates.h +++ b/include/deal.II/matrix_free/shape_info.templates.h @@ -61,6 +61,8 @@ namespace internal return a; } + + template Number get_first_array_element(const VectorizedArray a) @@ -68,6 +70,8 @@ namespace internal return a[0]; } + + template ShapeInfo::ShapeInfo() : element_type(tensor_general) @@ -81,6 +85,48 @@ 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); + 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 template void diff --git a/include/deal.II/numerics/vector_tools_project.templates.h b/include/deal.II/numerics/vector_tools_project.templates.h index 1988f3decd..69083a0bdc 100644 --- a/include/deal.II/numerics/vector_tools_project.templates.h +++ b/include/deal.II/numerics/vector_tools_project.templates.h @@ -706,7 +706,8 @@ namespace VectorTools // If we can, use the matrix-free implementation bool use_matrix_free = MatrixFree::is_supported( - dof.get_fe()); + dof.get_fe()) && + dof.get_fe().n_base_elements() == 1; // enforce_zero_boundary and project_to_boundary_first // are not yet supported. diff --git a/source/matrix_free/shape_info.inst.in b/source/matrix_free/shape_info.inst.in index f2c859da1d..293099d241 100644 --- a/source/matrix_free/shape_info.inst.in +++ b/source/matrix_free/shape_info.inst.in @@ -23,6 +23,15 @@ for (deal_II_dimension : DIMENSIONS; deal_II_scalar : REAL_SCALARS) const unsigned int); } +for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) + { +#if deal_II_dimension <= deal_II_space_dimension + template bool + internal::MatrixFreeFunctions::ShapeInfo::is_supported( + const FiniteElement &); +#endif + } + for (deal_II_dimension : DIMENSIONS; -- 2.39.5