From 0f8c583c9d1259e2f15859173935eafa4350ba19 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Wed, 20 Mar 2024 16:23:40 +0100 Subject: [PATCH] FEPointEvaluation: Work towards dim < spacedim --- include/deal.II/base/derivative_form.h | 32 +- .../deal.II/base/floating_point_comparator.h | 2 +- .../deal.II/matrix_free/fe_point_evaluation.h | 327 ++++++------------ .../matrix_free/fe_remote_evaluation.h | 4 +- .../matrix_free/shape_info.templates.h | 72 ++-- .../mg_transfer_global_coarsening.templates.h | 8 +- include/deal.II/non_matching/mapping_info.h | 24 +- source/matrix_free/fe_point_evaluation.cc | 3 - source/matrix_free/shape_info.inst.in | 43 ++- 9 files changed, 207 insertions(+), 308 deletions(-) diff --git a/include/deal.II/base/derivative_form.h b/include/deal.II/base/derivative_form.h index 624f1781b1..455bc11188 100644 --- a/include/deal.II/base/derivative_form.h +++ b/include/deal.II/base/derivative_form.h @@ -609,9 +609,7 @@ template inline DerivativeForm<1, spacedim, dim, Number> transpose(const DerivativeForm<1, dim, spacedim, Number> &DF) { - DerivativeForm<1, spacedim, dim, Number> tt; - tt = DF.transpose(); - return tt; + return DF.transpose(); } @@ -691,6 +689,34 @@ apply_diagonal_transformation( } + +/** + * Similar to the previous apply_transformation(). + * Each row of the result corresponds to one of the rows of @p D_X transformed + * by @p grad_F, equivalent to $\mathrm{D\_X} \, \mathrm{grad\_F}^T$ in matrix + * notation. + * + * @relatesalso DerivativeForm + */ +// rank=2 +template +inline DerivativeForm<1, + spacedim, + dim, + typename ProductType::type> +apply_diagonal_transformation( + const DerivativeForm<1, dim, spacedim, Number1> &grad_F, + const Tensor<2, dim, Number2> &D_X) +{ + DerivativeForm<1, spacedim, dim, typename ProductType::type> + dest; + for (unsigned int i = 0; i < dim; ++i) + dest[i] = apply_diagonal_transformation(grad_F, D_X[i]); + + return dest; +} + + DEAL_II_NAMESPACE_CLOSE #endif diff --git a/include/deal.II/base/floating_point_comparator.h b/include/deal.II/base/floating_point_comparator.h index 9a1d877457..c3a468784a 100644 --- a/include/deal.II/base/floating_point_comparator.h +++ b/include/deal.II/base/floating_point_comparator.h @@ -232,7 +232,7 @@ FloatingPointComparator::compare( const DerivativeForm &t1, const DerivativeForm &t2) const { - for (unsigned int i = 0; i < dim; ++i) + for (unsigned int i = 0; i < spacedim; ++i) { const ComparisonResult result = compare(t1[i], t2[i]); if (result != ComparisonResult::equal) diff --git a/include/deal.II/matrix_free/fe_point_evaluation.h b/include/deal.II/matrix_free/fe_point_evaluation.h index ca0c20726a..357d1922b5 100644 --- a/include/deal.II/matrix_free/fe_point_evaluation.h +++ b/include/deal.II/matrix_free/fe_point_evaluation.h @@ -56,7 +56,11 @@ namespace internal * Struct to distinguish between the value and gradient types of different * numbers of components used by the FlexibleEvaluator class. */ - template + template struct EvaluatorTypeTraits { using ScalarNumber = @@ -69,6 +73,10 @@ namespace internal using vectorized_value_type = Tensor<1, n_components, VectorizedArrayType>; using gradient_type = Tensor<1, n_components, Tensor<1, dim, Number>>; + using real_gradient_type = std::conditional_t< + n_components == spacedim, + Tensor<2, spacedim, Number>, + Tensor<1, n_components, Tensor<1, spacedim, Number>>>; using scalar_gradient_type = Tensor<1, n_components, Tensor<1, dim, ScalarNumber>>; using vectorized_gradient_type = @@ -134,10 +142,22 @@ namespace internal } static void - set_zero_gradient(gradient_type &value, const unsigned int vector_lane) + get_gradient(interface_vectorized_gradient_type &value, + const unsigned int vector_lane, + const DerivativeForm<1, dim, n_components, Number> &result) { for (unsigned int i = 0; i < n_components; ++i) for (unsigned int d = 0; d < dim; ++d) + internal::VectorizedArrayTrait::get_from_vectorized( + value[d][i], vector_lane) = result[i][d]; + } + + static void + set_zero_gradient(real_gradient_type &value, + const unsigned int vector_lane) + { + for (unsigned int i = 0; i < n_components; ++i) + for (unsigned int d = 0; d < spacedim; ++d) internal::VectorizedArrayTrait::get(value[i][d], vector_lane) = 0.; } @@ -204,24 +224,24 @@ namespace internal } static void - access(gradient_type &value, - const unsigned int vector_lane, - const unsigned int component, - const Tensor<1, dim, ScalarNumber> &shape_gradient) + access(real_gradient_type &value, + const unsigned int vector_lane, + const unsigned int component, + const Tensor<1, spacedim, ScalarNumber> &shape_gradient) { - for (unsigned int d = 0; d < dim; ++d) + for (unsigned int d = 0; d < spacedim; ++d) internal::VectorizedArrayTrait::get(value[component][d], vector_lane) += shape_gradient[d]; } - static Tensor<1, dim, ScalarNumber> - access(const gradient_type &value, - const unsigned int vector_lane, - const unsigned int component) + static Tensor<1, spacedim, ScalarNumber> + access(const real_gradient_type &value, + const unsigned int vector_lane, + const unsigned int component) { - Tensor<1, dim, ScalarNumber> result; - for (unsigned int d = 0; d < dim; ++d) + Tensor<1, spacedim, ScalarNumber> result; + for (unsigned int d = 0; d < spacedim; ++d) result[d] = internal::VectorizedArrayTrait::get(value[component][d], vector_lane); @@ -229,8 +249,8 @@ namespace internal } }; - template - struct EvaluatorTypeTraits + template + struct EvaluatorTypeTraits { using ScalarNumber = typename internal::VectorizedArrayTrait::value_type; @@ -241,6 +261,7 @@ namespace internal using scalar_value_type = ScalarNumber; using vectorized_value_type = VectorizedArrayType; using gradient_type = Tensor<1, dim, Number>; + using real_gradient_type = Tensor<1, spacedim, Number>; using scalar_gradient_type = Tensor<1, dim, ScalarNumber>; using vectorized_gradient_type = Tensor<1, dim, VectorizedArrayType>; using interface_vectorized_gradient_type = vectorized_gradient_type; @@ -306,9 +327,10 @@ namespace internal } static void - set_zero_gradient(gradient_type &value, const unsigned int vector_lane) + set_zero_gradient(real_gradient_type &value, + const unsigned int vector_lane) { - for (unsigned int d = 0; d < dim; ++d) + for (unsigned int d = 0; d < spacedim; ++d) internal::VectorizedArrayTrait::get(value[d], vector_lane) = 0.; } @@ -370,23 +392,23 @@ namespace internal } static void - access(gradient_type &value, - const unsigned int vector_lane, + access(real_gradient_type &value, + const unsigned int vector_lane, const unsigned int, - const scalar_gradient_type &shape_gradient) + const Tensor<1, spacedim, ScalarNumber> &shape_gradient) { - for (unsigned int d = 0; d < dim; ++d) + for (unsigned int d = 0; d < spacedim; ++d) internal::VectorizedArrayTrait::get(value[d], vector_lane) += shape_gradient[d]; } - static scalar_gradient_type - access(const gradient_type &value, - const unsigned int vector_lane, + static Tensor<1, spacedim, ScalarNumber> + access(const real_gradient_type &value, + const unsigned int vector_lane, const unsigned int) { - scalar_gradient_type result; - for (unsigned int d = 0; d < dim; ++d) + Tensor<1, spacedim, ScalarNumber> result; + for (unsigned int d = 0; d < spacedim; ++d) result[d] = internal::VectorizedArrayTrait::get(value[d], vector_lane); return result; @@ -394,7 +416,11 @@ namespace internal }; template - struct EvaluatorTypeTraits + struct EvaluatorTypeTraits> { using ScalarNumber = typename internal::VectorizedArrayTrait::value_type; @@ -405,6 +431,7 @@ namespace internal using scalar_value_type = Tensor<1, dim, ScalarNumber>; using vectorized_value_type = Tensor<1, dim, VectorizedArrayType>; using gradient_type = Tensor<2, dim, Number>; + using real_gradient_type = gradient_type; using scalar_gradient_type = Tensor<2, dim, ScalarNumber>; using vectorized_gradient_type = Tensor<2, dim, VectorizedArrayType>; using interface_vectorized_gradient_type = @@ -538,7 +565,7 @@ namespace internal } static void - access(gradient_type &value, + access(real_gradient_type &value, const unsigned int vector_lane, const unsigned int component, const Tensor<1, dim, ScalarNumber> &shape_gradient) @@ -550,9 +577,9 @@ namespace internal } static Tensor<1, dim, ScalarNumber> - access(const gradient_type &value, - const unsigned int vector_lane, - const unsigned int component) + access(const real_gradient_type &value, + const unsigned int vector_lane, + const unsigned int component) { Tensor<1, dim, ScalarNumber> result; for (unsigned int d = 0; d < dim; ++d) @@ -563,164 +590,6 @@ namespace internal } }; - template - struct EvaluatorTypeTraits<1, 1, Number> - { - using ScalarNumber = - typename internal::VectorizedArrayTrait::value_type; - using VectorizedArrayType = - typename dealii::internal::VectorizedArrayTrait< - Number>::vectorized_value_type; - using value_type = Number; - using scalar_value_type = ScalarNumber; - using vectorized_value_type = VectorizedArrayType; - using gradient_type = Tensor<1, 1, Number>; - using scalar_gradient_type = Tensor<1, 1, ScalarNumber>; - using vectorized_gradient_type = Tensor<1, 1, VectorizedArrayType>; - using interface_vectorized_gradient_type = vectorized_gradient_type; - - static void - read_value(const ScalarNumber vector_entry, - const unsigned int, - scalar_value_type &result) - { - result = vector_entry; - } - - static scalar_value_type - sum_value(const scalar_value_type &result) - { - return result; - } - - static scalar_value_type - sum_value(const vectorized_value_type &result) - { - return result.sum(); - } - - static ScalarNumber - sum_value(const unsigned int, const vectorized_value_type &result) - { - return result.sum(); - } - - static void - set_gradient(const vectorized_gradient_type &value, - const unsigned int vector_lane, - scalar_gradient_type &result) - { - result[0] = value[0][vector_lane]; - } - - static void - set_gradient(const vectorized_gradient_type &value, - const unsigned int, - vectorized_gradient_type &result) - { - result = value; - } - - static void - get_gradient(vectorized_gradient_type &value, - const unsigned int vector_lane, - const scalar_gradient_type &result) - { - value[0][vector_lane] = result[0]; - } - - static void - get_gradient(vectorized_gradient_type &value, - const unsigned int, - const vectorized_gradient_type &result) - { - value = result; - } - - static void - set_zero_gradient(gradient_type &value, const unsigned int vector_lane) - { - internal::VectorizedArrayTrait::get(value[0], vector_lane) = 0.; - } - - static void - set_value(const vectorized_value_type &value, - const unsigned int vector_lane, - scalar_value_type &result) - { - result = value[vector_lane]; - } - - static void - set_value(const vectorized_value_type &value, - const unsigned int, - vectorized_value_type &result) - { - result = value; - } - - static void - get_value(vectorized_value_type &value, - const unsigned int vector_lane, - const scalar_value_type &result) - { - value[vector_lane] = result; - } - - static void - get_value(vectorized_value_type &value, - const unsigned int, - const vectorized_value_type &result) - { - value = result; - } - - static void - set_zero_value(value_type &value, const unsigned int vector_lane) - { - internal::VectorizedArrayTrait::get(value, vector_lane) = 0.; - } - - static void - access(value_type &value, - const unsigned int vector_lane, - const unsigned int, - const ScalarNumber &shape_value) - { - internal::VectorizedArrayTrait::get(value, vector_lane) += - shape_value; - } - - static ScalarNumber - access(const value_type &value, - const unsigned int vector_lane, - const unsigned int) - { - return internal::VectorizedArrayTrait::get(value, vector_lane); - } - - static void - access(gradient_type &value, - const unsigned int vector_lane, - const unsigned int, - const scalar_gradient_type &shape_gradient) - { - internal::VectorizedArrayTrait::get(value[0], vector_lane) += - shape_gradient[0]; - } - - static scalar_gradient_type - access(const gradient_type &value, - const unsigned int vector_lane, - const unsigned int) - { - scalar_gradient_type result; - result[0] = - internal::VectorizedArrayTrait::get(value[0], vector_lane); - return result; - } - }; - template bool is_fast_path_supported(const FiniteElement &fe, @@ -761,11 +630,11 @@ public: using VectorizedArrayType = typename dealii::internal::VectorizedArrayTrait< Number>::vectorized_value_type; using ETT = typename internal::FEPointEvaluation:: - EvaluatorTypeTraits; + EvaluatorTypeTraits; using value_type = typename ETT::value_type; using scalar_value_type = typename ETT::scalar_value_type; using vectorized_value_type = typename ETT::vectorized_value_type; - using gradient_type = typename ETT::gradient_type; + using gradient_type = typename ETT::real_gradient_type; using interface_vectorized_gradient_type = typename ETT::interface_vectorized_gradient_type; @@ -788,10 +657,10 @@ protected: * objects, this parameter allows to select a range of `n_components` * components starting from this parameter. */ - FEPointEvaluationBase(const Mapping &mapping, - const FiniteElement &fe, - const UpdateFlags update_flags, - const unsigned int first_selected_component = 0); + FEPointEvaluationBase(const Mapping &mapping, + const FiniteElement &fe, + const UpdateFlags update_flags, + const unsigned int first_selected_component = 0); /** * Constructor to make the present class able to re-use the geometry @@ -814,7 +683,7 @@ protected: */ FEPointEvaluationBase( NonMatching::MappingInfo &mapping_info, - const FiniteElement &fe, + const FiniteElement &fe, const unsigned int first_selected_component = 0, const bool is_interior = true); @@ -995,7 +864,7 @@ protected: /** * Pointer to the FiniteElement object passed to the constructor. */ - SmartPointer> fe; + SmartPointer> fe; /** * Description of the 1d polynomial basis for tensor product elements used @@ -1237,11 +1106,12 @@ public: using VectorizedArrayType = typename dealii::internal::VectorizedArrayTrait< Number>::vectorized_value_type; using ETT = typename internal::FEPointEvaluation:: - EvaluatorTypeTraits; + EvaluatorTypeTraits; using value_type = typename ETT::value_type; using scalar_value_type = typename ETT::scalar_value_type; using vectorized_value_type = typename ETT::vectorized_value_type; - using gradient_type = typename ETT::gradient_type; + using unit_gradient_type = typename ETT::gradient_type; + using gradient_type = typename ETT::real_gradient_type; using interface_vectorized_gradient_type = typename ETT::interface_vectorized_gradient_type; @@ -1263,10 +1133,10 @@ public: * objects, this parameter allows to select a range of `n_components` * components starting from this parameter. */ - FEPointEvaluation(const Mapping &mapping, - const FiniteElement &fe, - const UpdateFlags update_flags, - const unsigned int first_selected_component = 0); + FEPointEvaluation(const Mapping &mapping, + const FiniteElement &fe, + const UpdateFlags update_flags, + const unsigned int first_selected_component = 0); /** * Constructor to make the present class able to re-use the geometry @@ -1286,7 +1156,7 @@ public: */ FEPointEvaluation( NonMatching::MappingInfo &mapping_info, - const FiniteElement &fe, + const FiniteElement &fe, const unsigned int first_selected_component = 0); /** @@ -1622,11 +1492,12 @@ public: using VectorizedArrayType = typename dealii::internal::VectorizedArrayTrait< Number>::vectorized_value_type; using ETT = typename internal::FEPointEvaluation:: - EvaluatorTypeTraits; + EvaluatorTypeTraits; using value_type = typename ETT::value_type; using scalar_value_type = typename ETT::scalar_value_type; using vectorized_value_type = typename ETT::vectorized_value_type; - using gradient_type = typename ETT::gradient_type; + using unit_gradient_type = typename ETT::gradient_type; + using gradient_type = typename ETT::real_gradient_type; using interface_vectorized_gradient_type = typename ETT::interface_vectorized_gradient_type; @@ -1635,7 +1506,7 @@ public: */ FEFacePointEvaluation( NonMatching::MappingInfo &mapping_info, - const FiniteElement &fe, + const FiniteElement &fe, const bool is_interior = true, const unsigned int first_selected_component = 0); @@ -1869,15 +1740,17 @@ private: const bool sum_into_values); }; + + // ----------------------- template and inline function ---------------------- template FEPointEvaluationBase:: - FEPointEvaluationBase(const Mapping &mapping, - const FiniteElement &fe, - const UpdateFlags update_flags, - const unsigned int first_selected_component) + FEPointEvaluationBase(const Mapping &mapping, + const FiniteElement &fe, + const UpdateFlags update_flags, + const unsigned int first_selected_component) : n_q_batches(numbers::invalid_unsigned_int) , n_q_points(numbers::invalid_unsigned_int) , n_q_points_scalar(numbers::invalid_unsigned_int) @@ -1904,7 +1777,7 @@ template FEPointEvaluationBase:: FEPointEvaluationBase( NonMatching::MappingInfo &mapping_info, - const FiniteElement &fe, + const FiniteElement &fe, const unsigned int first_selected_component, const bool is_interior) : n_q_batches(numbers::invalid_unsigned_int) @@ -1972,9 +1845,7 @@ FEPointEvaluationBase:: template FEPointEvaluationBase:: - FEPointEvaluationBase( - FEPointEvaluationBase - &&other) noexcept + FEPointEvaluationBase(FEPointEvaluationBase &&other) noexcept : n_q_batches(other.n_q_batches) , n_q_points(other.n_q_points) , n_q_points_scalar(other.n_q_points_scalar) @@ -2373,7 +2244,7 @@ FEPointEvaluationBase:: template FEPointEvaluation::FEPointEvaluation( NonMatching::MappingInfo &mapping_info, - const FiniteElement &fe, + const FiniteElement &fe, const unsigned int first_selected_component) : FEPointEvaluationBase( mapping_info, @@ -2385,10 +2256,10 @@ FEPointEvaluation::FEPointEvaluation( template FEPointEvaluation::FEPointEvaluation( - const Mapping &mapping, - const FiniteElement &fe, - const UpdateFlags update_flags, - const unsigned int first_selected_component) + const Mapping &mapping, + const FiniteElement &fe, + const UpdateFlags update_flags, + const unsigned int first_selected_component) : FEPointEvaluationBase( mapping, fe, @@ -2743,13 +2614,13 @@ FEPointEvaluation::evaluate_fast( v < stride && (stride == 1 || offset < this->n_q_points_scalar); ++v, ++offset) { - gradient_type unit_gradient; + unit_gradient_type unit_gradient; ETT::set_gradient(gradient, v, unit_gradient); this->gradients[offset] = this->cell_type <= internal::MatrixFreeFunctions::GeometryType::cartesian ? - apply_diagonal_transformation(this->inverse_jacobian_ptr[0], - unit_gradient) : + apply_diagonal_transformation( + this->inverse_jacobian_ptr[0].transpose(), unit_gradient) : apply_transformation( this ->inverse_jacobian_ptr[this->cell_type <= @@ -2997,7 +2868,7 @@ FEPointEvaluation::integrate_fast( v < stride && (stride == 1 || offset < this->n_q_points_scalar); ++v, ++offset) { - const auto grad_w = + const gradient_type grad_w = do_JxW ? this->gradients[offset] * this->JxW_ptr[offset] : this->gradients[offset]; ETT::get_gradient( @@ -3178,7 +3049,7 @@ template FEFacePointEvaluation:: FEFacePointEvaluation( NonMatching::MappingInfo &mapping_info, - const FiniteElement &fe, + const FiniteElement &fe, const bool is_interior, const unsigned int first_selected_component) : FEPointEvaluationBase( @@ -3665,13 +3536,13 @@ FEFacePointEvaluation:: v < stride && (stride == 1 || offset < this->n_q_points_scalar); ++v, ++offset) { - gradient_type unit_gradient; + unit_gradient_type unit_gradient; ETT::set_gradient(gradient, v, unit_gradient); this->gradients[offset] = this->cell_type <= internal::MatrixFreeFunctions::GeometryType::cartesian ? - apply_diagonal_transformation(this->inverse_jacobian_ptr[0], - unit_gradient) : + apply_diagonal_transformation( + this->inverse_jacobian_ptr[0].transpose(), unit_gradient) : apply_transformation( this ->inverse_jacobian_ptr[this->cell_type <= diff --git a/include/deal.II/matrix_free/fe_remote_evaluation.h b/include/deal.II/matrix_free/fe_remote_evaluation.h index 930dd9ef17..c3e512070e 100644 --- a/include/deal.II/matrix_free/fe_remote_evaluation.h +++ b/include/deal.II/matrix_free/fe_remote_evaluation.h @@ -42,10 +42,10 @@ namespace internal struct PrecomputedEvaluationData { using value_type = typename internal::FEPointEvaluation:: - EvaluatorTypeTraits::value_type; + EvaluatorTypeTraits::value_type; using gradient_type = typename internal::FEPointEvaluation:: - EvaluatorTypeTraits::gradient_type; + EvaluatorTypeTraits::gradient_type; /** * Values at quadrature points. diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h index dec73fa992..2eafca6dbf 100644 --- a/include/deal.II/matrix_free/shape_info.templates.h +++ b/include/deal.II/matrix_free/shape_info.templates.h @@ -84,23 +84,23 @@ namespace internal - template + template void get_element_type_specific_information( - const FiniteElement &fe_in, - const FiniteElement &fe, - const unsigned int base_element_number, - ElementType &element_type, - std::vector &scalar_lexicographic, - std::vector &lexicographic_numbering) + const FiniteElement &fe_in, + const FiniteElement &fe, + const unsigned int base_element_number, + ElementType &element_type, + std::vector &scalar_lexicographic, + std::vector &lexicographic_numbering) { element_type = tensor_general; - const auto fe_poly = dynamic_cast *>(&fe); + const auto fe_poly = dynamic_cast *>(&fe); - if (dynamic_cast *>(&fe) != nullptr || - dynamic_cast *>(&fe) != nullptr || - dynamic_cast *>(&fe) != nullptr) + if (dynamic_cast *>(&fe) != nullptr || + dynamic_cast *>(&fe) != nullptr || + dynamic_cast *>(&fe) != nullptr) { scalar_lexicographic.resize(fe.n_dofs_per_cell()); for (unsigned int i = 0; i < scalar_lexicographic.size(); ++i) @@ -115,14 +115,16 @@ namespace internal Polynomials::PiecewisePolynomial> *>( &fe_poly->get_poly_space()) != nullptr)) scalar_lexicographic = fe_poly->get_poly_space_numbering_inverse(); - else if (const auto fe_dgp = dynamic_cast *>(&fe)) + else if (const auto fe_dgp = + dynamic_cast *>(&fe)) { scalar_lexicographic.resize(fe_dgp->n_dofs_per_cell()); for (unsigned int i = 0; i < fe_dgp->n_dofs_per_cell(); ++i) scalar_lexicographic[i] = i; element_type = truncated_tensor; } - else if (const auto fe_q_dg0 = dynamic_cast *>(&fe)) + else if (const auto fe_q_dg0 = + dynamic_cast *>(&fe)) { scalar_lexicographic = fe_q_dg0->get_poly_space_numbering_inverse(); element_type = tensor_symmetric_plus_dg0; @@ -190,7 +192,7 @@ namespace internal Assert(fe_name[template_starts + 1] == (dim == 1 ? '1' : (dim == 2 ? '2' : '3')), ExcInternalError()); - fe_name[template_starts + 1] = std::to_string(dim_to).c_str()[0]; + fe_name[template_starts + 1] = std::to_string(dim_to)[0]; } return FETools::get_fe_by_name(fe_name); } @@ -237,9 +239,6 @@ namespace internal const FiniteElement &fe_in, const unsigned int base_element_number) { - static_assert(dim == spacedim, - "Currently, only the case dim=spacedim is implemented"); - // ShapeInfo for RT elements. Here, data is of size 2 instead of 1. // data[0] is univariate_shape_data in normal direction and // data[1] is univariate_shape_data in tangential direction @@ -251,7 +250,7 @@ namespace internal const auto quad = quad_in.get_tensor_basis()[0]; - const FiniteElement &fe = + const FiniteElement &fe = fe_in.base_element(base_element_number); n_dimensions = dim; n_components = fe_in.n_components(); @@ -312,11 +311,11 @@ namespace internal return; } else if (quad_in.is_tensor_product() == false || - dynamic_cast *>( + dynamic_cast *>( &fe_in.base_element(base_element_number)) != nullptr || - dynamic_cast *>( + dynamic_cast *>( &fe_in.base_element(base_element_number)) != nullptr || - dynamic_cast *>( + dynamic_cast *>( &fe_in.base_element(base_element_number)) != nullptr) { // specialization for arbitrary finite elements and quadrature rules @@ -504,9 +503,10 @@ namespace internal const auto quad = quad_in.get_tensor_basis()[0]; - const FiniteElement &fe = fe_in.base_element(base_element_number); - n_dimensions = dim; - n_components = fe_in.n_components(); + const FiniteElement &fe = + fe_in.base_element(base_element_number); + n_dimensions = dim; + n_components = fe_in.n_components(); Assert(fe.n_components() == 1, ExcMessage("FEEvaluation only works for scalar finite elements.")); @@ -553,8 +553,8 @@ namespace internal scalar_lexicographic, 0); - if (dim > 1 && (dynamic_cast *>(&fe) || - dynamic_cast *>(&fe))) + if (dim > 1 && (dynamic_cast *>(&fe) || + dynamic_cast *>(&fe))) { auto &subface_interpolation_matrix_0 = univariate_shape_data.subface_interpolation_matrices[0]; @@ -571,7 +571,8 @@ namespace internal subface_interpolation_matrix_scalar_0.resize(nn * nn); subface_interpolation_matrix_scalar_1.resize(nn * nn); - const bool is_feq = dynamic_cast *>(&fe) != nullptr; + const bool is_feq = + dynamic_cast *>(&fe) != nullptr; std::vector> fe_q_points = is_feq ? QGaussLobatto<1>(nn).get_points() : @@ -613,7 +614,7 @@ namespace internal if (element_type == tensor_general && univariate_shape_data.check_and_set_shapes_symmetric()) { - if (dynamic_cast *>(&fe) && + if (dynamic_cast *>(&fe) && fe.tensor_degree() > 1) element_type = tensor_symmetric_no_collocation; else if (univariate_shape_data.check_shapes_collocation()) @@ -1151,9 +1152,6 @@ namespace internal bool ShapeInfo::is_supported(const FiniteElement &fe) { - if (dim != spacedim) - return false; - if (dynamic_cast *>(&fe)) return true; @@ -1170,12 +1168,12 @@ namespace internal 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) != - nullptr || - dynamic_cast *>(fe_poly_ptr) != - nullptr || - dynamic_cast *>(fe_poly_ptr) != - nullptr) + if (dynamic_cast *>( + fe_poly_ptr) != nullptr || + dynamic_cast *>( + fe_poly_ptr) != nullptr || + dynamic_cast *>( + fe_poly_ptr) != nullptr) return true; if (dynamic_cast::prolongate_and_add_internal_comp( VectorType &dst, const VectorType &src) const { - using Traits = - internal::FEPointEvaluation::EvaluatorTypeTraits; + using Traits = internal::FEPointEvaluation:: + EvaluatorTypeTraits; using value_type = typename Traits::value_type; std::vector evaluation_point_results; @@ -5023,8 +5023,8 @@ MGTwoLevelTransferNonNested::restrict_and_add_internal_comp( VectorType &dst, const VectorType &src) const { - using Traits = - internal::FEPointEvaluation::EvaluatorTypeTraits; + using Traits = internal::FEPointEvaluation:: + EvaluatorTypeTraits; using value_type = typename Traits::value_type; std::vector evaluation_point_results; diff --git a/include/deal.II/non_matching/mapping_info.h b/include/deal.II/non_matching/mapping_info.h index ea51fa8435..1742fd46d3 100644 --- a/include/deal.II/non_matching/mapping_info.h +++ b/include/deal.II/non_matching/mapping_info.h @@ -186,7 +186,7 @@ namespace NonMatching dealii::internal::MatrixFreeFunctions::GeometryType compute_geometry_type( const double diameter, - const std::vector> + const std::vector> &inverse_jacobians) { const auto jac_0 = inverse_jacobians[0]; @@ -195,7 +195,7 @@ namespace NonMatching bool jacobian_constant = true; for (unsigned int q = 1; q < inverse_jacobians.size(); ++q) { - const DerivativeForm<1, dim, spacedim> &jac = inverse_jacobians[q]; + const DerivativeForm<1, spacedim, dim> &jac = inverse_jacobians[q]; for (unsigned int d = 0; d < dim; ++d) for (unsigned int e = 0; e < spacedim; ++e) if (std::fabs(jac_0[d][e] - jac[d][e]) > zero_tolerance_double) @@ -206,7 +206,7 @@ namespace NonMatching // check whether the Jacobian is diagonal to machine // accuracy - bool cell_cartesian = jacobian_constant; + bool cell_cartesian = jacobian_constant && dim == spacedim; for (unsigned int d = 0; d < dim; ++d) for (unsigned int e = 0; e < dim; ++e) if (d != e) @@ -296,8 +296,8 @@ namespace NonMatching * @param additional_data Additional data for the class to specify the * behavior during reinitialization. */ - MappingInfo(const Mapping &mapping, - const UpdateFlags update_flags, + MappingInfo(const Mapping &mapping, + const UpdateFlags update_flags, const AdditionalData additional_data = AdditionalData()); /** @@ -448,7 +448,7 @@ namespace NonMatching * Getter function for real points. The offset can be obtained with * compute_data_index_offset(). */ - const Point * + const Point * get_real_point(const unsigned int offset) const; /** @@ -797,9 +797,9 @@ namespace NonMatching template MappingInfo::MappingInfo( - const Mapping &mapping, - const UpdateFlags update_flags, - const AdditionalData additional_data) + const Mapping &mapping, + const UpdateFlags update_flags, + const AdditionalData additional_data) : mapping(&mapping) , update_flags(update_flags) , update_flags_mapping(update_default) @@ -1960,8 +1960,8 @@ namespace NonMatching v) = mapping_data.jacobians[q * n_lanes + v][d][s]; if (update_flags_mapping & UpdateFlags::update_inverse_jacobians) - for (unsigned int d = 0; d < dim; ++d) - for (unsigned int s = 0; s < spacedim; ++s) + for (unsigned int d = 0; d < spacedim; ++d) + for (unsigned int s = 0; s < dim; ++s) dealii::internal::VectorizedArrayTrait::get( inverse_jacobians[is_interior ? 0 : 1] [compressed_offset][s][d], @@ -2076,7 +2076,7 @@ namespace NonMatching template - inline const Point * + inline const Point * MappingInfo::get_real_point( const unsigned int offset) const { diff --git a/source/matrix_free/fe_point_evaluation.cc b/source/matrix_free/fe_point_evaluation.cc index 958e1ef2e1..149569b8fa 100644 --- a/source/matrix_free/fe_point_evaluation.cc +++ b/source/matrix_free/fe_point_evaluation.cc @@ -43,9 +43,6 @@ namespace internal { // check if supported const bool flag = [&]() { - if (dim != spacedim) - return false; - const FiniteElement *fe_ptr = &(fe.base_element(base_element_number)); if (fe_ptr->n_components() != 1) diff --git a/source/matrix_free/shape_info.inst.in b/source/matrix_free/shape_info.inst.in index c07dcea0e8..d013a949e0 100644 --- a/source/matrix_free/shape_info.inst.in +++ b/source/matrix_free/shape_info.inst.in @@ -13,30 +13,35 @@ // ------------------------------------------------------------------------ -for (deal_II_dimension : DIMENSIONS; deal_II_scalar : REAL_SCALARS) +for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : DIMENSIONS; + deal_II_scalar : REAL_SCALARS) { - template internal::MatrixFreeFunctions::ShapeInfo:: - ShapeInfo(const Quadrature &, - const FiniteElement &, - const unsigned int); +#if deal_II_dimension <= deal_II_space_dimension + template internal::MatrixFreeFunctions::ShapeInfo< + deal_II_scalar>::ShapeInfo(const Quadrature &, + const FiniteElement &, + const unsigned int); template void internal::MatrixFreeFunctions::ShapeInfo::reinit( const Quadrature &, - const FiniteElement &, + const FiniteElement &, const unsigned int); -#if deal_II_dimension > 1 - template internal::MatrixFreeFunctions::ShapeInfo:: - ShapeInfo(const Quadrature<1> &, - const FiniteElement &, - const unsigned int); +# if deal_II_dimension > 1 + template internal::MatrixFreeFunctions::ShapeInfo< + deal_II_scalar>::ShapeInfo(const Quadrature<1> &, + const FiniteElement &, + const unsigned int); template void internal::MatrixFreeFunctions::ShapeInfo::reinit( const Quadrature<1> &, - const FiniteElement &, + const FiniteElement &, const unsigned int); +# endif #endif } @@ -51,33 +56,35 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) -for (deal_II_dimension : DIMENSIONS; +for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : DIMENSIONS; deal_II_scalar_vectorized : REAL_SCALARS_VECTORIZED) { +#if deal_II_dimension <= deal_II_space_dimension template internal::MatrixFreeFunctions:: ShapeInfo::ShapeInfo( const Quadrature &, - const FiniteElement &, + const FiniteElement &, const unsigned int); template void internal::MatrixFreeFunctions::ShapeInfo::reinit( const Quadrature &, - const FiniteElement &, + const FiniteElement &, const unsigned int); -#if deal_II_dimension > 1 +# if deal_II_dimension > 1 template internal::MatrixFreeFunctions:: ShapeInfo::ShapeInfo( const Quadrature<1> &, - const FiniteElement &, + const FiniteElement &, const unsigned int); template void internal::MatrixFreeFunctions::ShapeInfo::reinit( const Quadrature<1> &, - const FiniteElement &, + const FiniteElement &, const unsigned int); +# endif #endif } -- 2.39.5