From ca25489d9bfebec583f6f6c8df355127d098a93e Mon Sep 17 00:00:00 2001 From: Natalia Nebulishvili Date: Fri, 20 Jun 2025 12:19:20 +0200 Subject: [PATCH] adapted FE_RT to polynomials_vector_anisotropic --- include/deal.II/fe/fe_raviart_thomas.h | 16 ++++ .../matrix_free/shape_info.templates.h | 1 - source/fe/fe_raviart_thomas.cc | 83 ++++++++++++++-- source/fe/fe_raviart_thomas_nodal.cc | 95 ++++++++++++++++--- 4 files changed, 176 insertions(+), 19 deletions(-) diff --git a/include/deal.II/fe/fe_raviart_thomas.h b/include/deal.II/fe/fe_raviart_thomas.h index 5f0464fb69..46541cdd61 100644 --- a/include/deal.II/fe/fe_raviart_thomas.h +++ b/include/deal.II/fe/fe_raviart_thomas.h @@ -144,6 +144,14 @@ public: virtual std::unique_ptr> clone() const override; + /** + * Compute the lexicographic to hierarchic numbering underlying this class, + * necessary for the creation of the respective vector polynomial space. + */ + std::vector + get_lexicographic_numbering(const unsigned int normal_degree, + const unsigned int tangential_degree) const; + /** * This function returns @p true, if the shape function @p shape_index has * non-zero function values somewhere on the face @p face_index. @@ -352,6 +360,14 @@ public: virtual std::unique_ptr> clone() const override; + /** + * Compute the lexicographic to hierarchic numbering underlying this class, + * necessary for the creation of the respective vector polynomial space. + */ + std::vector + get_lexicographic_numbering(const unsigned int normal_degree, + const unsigned int tangential_degree) const; + virtual void get_face_interpolation_matrix(const FiniteElement &source, FullMatrix &matrix, diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h index f3fc16f505..8c852e76ff 100644 --- a/include/deal.II/matrix_free/shape_info.templates.h +++ b/include/deal.II/matrix_free/shape_info.templates.h @@ -21,7 +21,6 @@ #include #include #include -#include #include #include #include diff --git a/source/fe/fe_raviart_thomas.cc b/source/fe/fe_raviart_thomas.cc index 7308a8ea93..cadd78a4ac 100644 --- a/source/fe/fe_raviart_thomas.cc +++ b/source/fe/fe_raviart_thomas.cc @@ -14,7 +14,7 @@ #include -#include +#include #include #include #include @@ -43,16 +43,20 @@ DEAL_II_NAMESPACE_OPEN template FE_RaviartThomas::FE_RaviartThomas(const unsigned int deg) : FE_PolyTensor( - PolynomialsRaviartThomas(deg), + PolynomialsVectorAnisotropic(deg + 1, + deg, + get_lexicographic_numbering(deg + 1, + deg)), FiniteElementData(get_dpo_vector(deg), dim, deg + 1, FiniteElementData::Hdiv), - std::vector(PolynomialsRaviartThomas::n_polynomials(deg), - true), - std::vector(PolynomialsRaviartThomas::n_polynomials( - deg), - ComponentMask(std::vector(dim, true)))) + std::vector( + PolynomialsVectorAnisotropic::n_polynomials(deg + 1, deg), + true), + std::vector( + PolynomialsVectorAnisotropic::n_polynomials(deg + 1, deg), + ComponentMask(std::vector(dim, true)))) { Assert(dim >= 2, ExcImpossibleInDim(dim)); const unsigned int n_dofs = this->n_dofs_per_cell(); @@ -694,6 +698,71 @@ FE_RaviartThomas::has_support_on_face(const unsigned int shape_index, +template +std::vector +FE_RaviartThomas::get_lexicographic_numbering( + const unsigned int normal_degree, + const unsigned int tangential_degree) const +{ + const unsigned int n_dofs_face = + Utilities::pow(tangential_degree + 1, dim - 1); + std::vector lexicographic_numbering; + // component 1 + for (unsigned int j = 0; j < n_dofs_face; ++j) + { + lexicographic_numbering.push_back(j); + if (normal_degree > 1) + for (unsigned int i = n_dofs_face * 2 * dim; + i < n_dofs_face * 2 * dim + normal_degree - 1; + ++i) + lexicographic_numbering.push_back(i + j * (normal_degree - 1)); + lexicographic_numbering.push_back(n_dofs_face + j); + } + + // component 2 + unsigned int layers = (dim == 3) ? tangential_degree + 1 : 1; + for (unsigned int k = 0; k < layers; ++k) + { + unsigned int k_add = k * (tangential_degree + 1); + for (unsigned int j = n_dofs_face * 2; + j < n_dofs_face * 2 + tangential_degree + 1; + ++j) + lexicographic_numbering.push_back(j + k_add); + + if (normal_degree > 1) + for (unsigned int i = n_dofs_face * (2 * dim + (normal_degree - 1)); + i < n_dofs_face * (2 * dim + (normal_degree - 1)) + + (normal_degree - 1) * (tangential_degree + 1); + ++i) + { + lexicographic_numbering.push_back(i + k_add * tangential_degree); + } + for (unsigned int j = n_dofs_face * 3; + j < n_dofs_face * 3 + tangential_degree + 1; + ++j) + lexicographic_numbering.push_back(j + k_add); + } + + // component 3 + if (dim == 3) + { + for (unsigned int i = 4 * n_dofs_face; i < 5 * n_dofs_face; ++i) + lexicographic_numbering.push_back(i); + if (normal_degree > 1) + for (unsigned int i = + 6 * n_dofs_face + n_dofs_face * 2 * (normal_degree - 1); + i < 6 * n_dofs_face + n_dofs_face * 3 * (normal_degree - 1); + ++i) + lexicographic_numbering.push_back(i); + for (unsigned int i = 5 * n_dofs_face; i < 6 * n_dofs_face; ++i) + lexicographic_numbering.push_back(i); + } + + return lexicographic_numbering; +} + + + template void FE_RaviartThomas::convert_generalized_support_point_values_to_dof_values( diff --git a/source/fe/fe_raviart_thomas_nodal.cc b/source/fe/fe_raviart_thomas_nodal.cc index 3b61ab611a..d259869236 100644 --- a/source/fe/fe_raviart_thomas_nodal.cc +++ b/source/fe/fe_raviart_thomas_nodal.cc @@ -14,7 +14,7 @@ #include -#include +#include #include #include @@ -66,24 +66,32 @@ namespace template FE_RaviartThomasNodal::FE_RaviartThomasNodal(const unsigned int degree) - : FE_PolyTensor(PolynomialsRaviartThomas(degree), - FiniteElementData(get_rt_dpo_vector(dim, degree), - dim, - degree + 1, - FiniteElementData::Hdiv), - std::vector(1, false), - std::vector( - PolynomialsRaviartThomas::n_polynomials(degree), - ComponentMask(std::vector(dim, true)))) + : FE_PolyTensor( + PolynomialsVectorAnisotropic(degree + 1, + degree, + get_lexicographic_numbering(degree + 1, + degree)), + FiniteElementData(get_rt_dpo_vector(dim, degree), + dim, + degree + 1, + FiniteElementData::Hdiv), + std::vector(1, false), + std::vector( + PolynomialsVectorAnisotropic::n_polynomials(degree + 1, degree), + ComponentMask(std::vector(dim, true)))) { Assert(dim >= 2, ExcImpossibleInDim(dim)); this->mapping_kind = {mapping_raviart_thomas}; + const std::vector numbering = + get_lexicographic_numbering(degree + 1, degree); + // First, initialize the generalized support points and quadrature weights, // since they are required for interpolation. this->generalized_support_points = - PolynomialsRaviartThomas(degree).get_polynomial_support_points(); + PolynomialsVectorAnisotropic(degree + 1, degree, numbering) + .get_polynomial_support_points(); AssertDimension(this->generalized_support_points.size(), this->n_dofs_per_cell()); @@ -234,6 +242,71 @@ FE_RaviartThomasNodal::has_support_on_face( +template +std::vector +FE_RaviartThomasNodal::get_lexicographic_numbering( + const unsigned int normal_degree, + const unsigned int tangential_degree) const +{ + const unsigned int n_dofs_face = + Utilities::pow(tangential_degree + 1, dim - 1); + std::vector lexicographic_numbering; + // component 1 + for (unsigned int j = 0; j < n_dofs_face; ++j) + { + lexicographic_numbering.push_back(j); + if (normal_degree > 1) + for (unsigned int i = n_dofs_face * 2 * dim; + i < n_dofs_face * 2 * dim + normal_degree - 1; + ++i) + lexicographic_numbering.push_back(i + j * (normal_degree - 1)); + lexicographic_numbering.push_back(n_dofs_face + j); + } + + // component 2 + unsigned int layers = (dim == 3) ? tangential_degree + 1 : 1; + for (unsigned int k = 0; k < layers; ++k) + { + unsigned int k_add = k * (tangential_degree + 1); + for (unsigned int j = n_dofs_face * 2; + j < n_dofs_face * 2 + tangential_degree + 1; + ++j) + lexicographic_numbering.push_back(j + k_add); + + if (normal_degree > 1) + for (unsigned int i = n_dofs_face * (2 * dim + (normal_degree - 1)); + i < n_dofs_face * (2 * dim + (normal_degree - 1)) + + (normal_degree - 1) * (tangential_degree + 1); + ++i) + { + lexicographic_numbering.push_back(i + k_add * tangential_degree); + } + for (unsigned int j = n_dofs_face * 3; + j < n_dofs_face * 3 + tangential_degree + 1; + ++j) + lexicographic_numbering.push_back(j + k_add); + } + + // component 3 + if (dim == 3) + { + for (unsigned int i = 4 * n_dofs_face; i < 5 * n_dofs_face; ++i) + lexicographic_numbering.push_back(i); + if (normal_degree > 1) + for (unsigned int i = + 6 * n_dofs_face + n_dofs_face * 2 * (normal_degree - 1); + i < 6 * n_dofs_face + n_dofs_face * 3 * (normal_degree - 1); + ++i) + lexicographic_numbering.push_back(i); + for (unsigned int i = 5 * n_dofs_face; i < 6 * n_dofs_face; ++i) + lexicographic_numbering.push_back(i); + } + + return lexicographic_numbering; +} + + + template void FE_RaviartThomasNodal:: -- 2.39.5