From e8feaffd7ab0d5a2f33aa200fea7ed8fda4cec46 Mon Sep 17 00:00:00 2001 From: Natalia Nebulishvili Date: Tue, 17 Jun 2025 17:45:47 +0200 Subject: [PATCH] Removed most of the duplicated code from Raviart-Thomas polynomials --- .../deal.II/base/polynomials_raviart_thomas.h | 62 +---- source/base/polynomials_raviart_thomas.cc | 237 +----------------- 2 files changed, 12 insertions(+), 287 deletions(-) diff --git a/include/deal.II/base/polynomials_raviart_thomas.h b/include/deal.II/base/polynomials_raviart_thomas.h index c93fd827f5..b7a618a6ce 100644 --- a/include/deal.II/base/polynomials_raviart_thomas.h +++ b/include/deal.II/base/polynomials_raviart_thomas.h @@ -23,7 +23,7 @@ #include #include #include -#include +#include #include #include @@ -46,7 +46,7 @@ DEAL_II_NAMESPACE_OPEN * @ingroup Polynomials */ template -class PolynomialsRaviartThomas : public TensorPolynomialsBase +class PolynomialsRaviartThomas : public PolynomialsVectorAnisotropic { public: /** @@ -80,19 +80,6 @@ public: * The size of the vectors must either be zero or equal n(). In * the first case, the function will not compute these values. */ - void - evaluate(const Point &unit_point, - std::vector> &values, - std::vector> &grads, - std::vector> &grad_grads, - std::vector> &third_derivatives, - std::vector> &fourth_derivatives) const override; - - /** - * Return the name of the space, which is PolynomialsRaviartThomas. - */ - std::string - name() const override; /** * Return the number of polynomials in the space without requiring to @@ -102,7 +89,7 @@ public: static unsigned int n_polynomials(const unsigned int normal_degree, const unsigned int tangential_degree); - + /** * Variant of the n_polynomials() function taking only a single argument * `degree`, assuming `degree + 1` in the normal direction and `degree` in @@ -124,49 +111,6 @@ public: */ virtual std::unique_ptr> clone() const override; - - /** - * Compute the generalized support points in the ordering used by the - * polynomial shape functions. Note that these points are not support points - * in the classical sense as the Lagrange polynomials of the different - * components have different points, which need to be combined in terms of - * Piola transforms. - */ - std::vector> - get_polynomial_support_points() const; - -private: - /** - * The given degree in the normal direction. - */ - const unsigned int normal_degree; - - /** - * The given degree in the tangential direction. - */ - const unsigned int tangential_degree; - - /** - * An object representing the polynomial space for a single component. We - * can re-use it by rotating the coordinates of the evaluation point. - */ - const AnisotropicPolynomials polynomial_space; - - /** - * Renumbering from lexicographic to hierarchic order. - */ - std::vector lexicographic_to_hierarchic; - - /** - * Renumbering from hierarchic to lexicographic order. Inverse of - * lexicographic_to_hierarchic. - */ - std::vector hierarchic_to_lexicographic; - - /** - * Renumbering from shifted polynomial spaces to lexicographic one. - */ - std::array, dim> renumber_aniso; }; diff --git a/source/base/polynomials_raviart_thomas.cc b/source/base/polynomials_raviart_thomas.cc index f9a50d266b..74ffe95e71 100644 --- a/source/base/polynomials_raviart_thomas.cc +++ b/source/base/polynomials_raviart_thomas.cc @@ -25,99 +25,14 @@ DEAL_II_NAMESPACE_OPEN -namespace -{ - // Create nodal Raviart-Thomas polynomials as the tensor product of Lagrange - // polynomials on Gauss-Lobatto points of the given degrees in the normal and - // tangential directions, respectively (we could also choose Lagrange - // polynomials on Gauss points but those are slightly more expensive to handle - // in classes). - std::vector>> - create_rt_polynomials(const unsigned int dim, - const unsigned int normal_degree, - const unsigned int tangential_degree) - { - std::vector>> pols(dim); - if (normal_degree > 0) - pols[0] = Polynomials::generate_complete_Lagrange_basis( - QGaussLobatto<1>(normal_degree + 1).get_points()); - else - pols[0] = Polynomials::generate_complete_Lagrange_basis( - QMidpoint<1>().get_points()); - if (tangential_degree > 0) - for (unsigned int d = 1; d < dim; ++d) - pols[d] = Polynomials::generate_complete_Lagrange_basis( - QGaussLobatto<1>(tangential_degree + 1).get_points()); - else - for (unsigned int d = 1; d < dim; ++d) - pols[d] = Polynomials::generate_complete_Lagrange_basis( - QMidpoint<1>().get_points()); - - return pols; - } -} // namespace - - template PolynomialsRaviartThomas::PolynomialsRaviartThomas( const unsigned int normal_degree, const unsigned int tangential_degree) - : TensorPolynomialsBase(std::min(normal_degree, tangential_degree), - n_polynomials(normal_degree, tangential_degree)) - , normal_degree(normal_degree) - , tangential_degree(tangential_degree) - , polynomial_space( - create_rt_polynomials(dim, normal_degree, tangential_degree)) -{ - // create renumbering of the unknowns from the lexicographic order to the - // actual order required by the finite element class with unknowns on - // faces placed first - const unsigned int n_pols = polynomial_space.n(); - lexicographic_to_hierarchic = - get_lexicographic_numbering(normal_degree, tangential_degree); - - hierarchic_to_lexicographic = - Utilities::invert_permutation(lexicographic_to_hierarchic); - - // since we only store an anisotropic polynomial for the first component, - // we set up a second numbering to switch out the actual coordinate - // direction - renumber_aniso[0].resize(n_pols); - for (unsigned int i = 0; i < n_pols; ++i) - renumber_aniso[0][i] = i; - if (dim == 2) - { - // switch x and y component (i and j loops) - renumber_aniso[1].resize(n_pols); - for (unsigned int j = 0; j < normal_degree + 1; ++j) - for (unsigned int i = 0; i < tangential_degree + 1; ++i) - renumber_aniso[1][j * (tangential_degree + 1) + i] = - j + i * (normal_degree + 1); - } - if (dim == 3) - { - // switch x, y, and z component (i, j, k) -> (j, k, i) - renumber_aniso[1].resize(n_pols); - for (unsigned int k = 0; k < tangential_degree + 1; ++k) - for (unsigned int j = 0; j < normal_degree + 1; ++j) - for (unsigned int i = 0; i < tangential_degree + 1; ++i) - renumber_aniso[1][(k * (normal_degree + 1) + j) * - (tangential_degree + 1) + - i] = - j + (normal_degree + 1) * (k + i * (tangential_degree + 1)); - - // switch x, y, and z component (i, j, k) -> (k, i, j) - renumber_aniso[2].resize(n_pols); - for (unsigned int k = 0; k < normal_degree + 1; ++k) - for (unsigned int j = 0; j < tangential_degree + 1; ++j) - for (unsigned int i = 0; i < tangential_degree + 1; ++i) - renumber_aniso[2][(k * (tangential_degree + 1) + j) * - (tangential_degree + 1) + - i] = - k + (normal_degree + 1) * (i + j * (tangential_degree + 1)); - } -} + : PolynomialsVectorAnisotropic(normal_degree, tangential_degree, + get_lexicographic_numbering(normal_degree, tangential_degree)) +{} @@ -128,125 +43,22 @@ PolynomialsRaviartThomas::PolynomialsRaviartThomas(const unsigned int k) -template -void -PolynomialsRaviartThomas::evaluate( - const Point &unit_point, - std::vector> &values, - std::vector> &grads, - std::vector> &grad_grads, - std::vector> &third_derivatives, - std::vector> &fourth_derivatives) const -{ - Assert(values.size() == this->n() || values.empty(), - ExcDimensionMismatch(values.size(), this->n())); - Assert(grads.size() == this->n() || grads.empty(), - ExcDimensionMismatch(grads.size(), this->n())); - Assert(grad_grads.size() == this->n() || grad_grads.empty(), - ExcDimensionMismatch(grad_grads.size(), this->n())); - Assert(third_derivatives.size() == this->n() || third_derivatives.empty(), - ExcDimensionMismatch(third_derivatives.size(), this->n())); - Assert(fourth_derivatives.size() == this->n() || fourth_derivatives.empty(), - ExcDimensionMismatch(fourth_derivatives.size(), this->n())); - - std::vector p_values; - std::vector> p_grads; - std::vector> p_grad_grads; - std::vector> p_third_derivatives; - std::vector> p_fourth_derivatives; - - const unsigned int n_sub = polynomial_space.n(); - p_values.resize((values.empty()) ? 0 : n_sub); - p_grads.resize((grads.empty()) ? 0 : n_sub); - p_grad_grads.resize((grad_grads.empty()) ? 0 : n_sub); - p_third_derivatives.resize((third_derivatives.empty()) ? 0 : n_sub); - p_fourth_derivatives.resize((fourth_derivatives.empty()) ? 0 : n_sub); - - for (unsigned int d = 0; d < dim; ++d) - { - // First we copy the point. The polynomial space for component d - // consists of polynomials of degree k in x_d and degree k+1 in the - // other variables. in order to simplify this, we use the same - // AnisotropicPolynomial space and simply rotate the coordinates - // through all directions. - Point p; - for (unsigned int c = 0; c < dim; ++c) - p[c] = unit_point[(c + d) % dim]; - - polynomial_space.evaluate(p, - p_values, - p_grads, - p_grad_grads, - p_third_derivatives, - p_fourth_derivatives); - - for (unsigned int i = 0; i < p_values.size(); ++i) - values[lexicographic_to_hierarchic[i + d * n_sub]][d] = - p_values[renumber_aniso[d][i]]; - - for (unsigned int i = 0; i < p_grads.size(); ++i) - for (unsigned int d1 = 0; d1 < dim; ++d1) - grads[lexicographic_to_hierarchic[i + d * n_sub]][d][(d1 + d) % dim] = - p_grads[renumber_aniso[d][i]][d1]; - - for (unsigned int i = 0; i < p_grad_grads.size(); ++i) - for (unsigned int d1 = 0; d1 < dim; ++d1) - for (unsigned int d2 = 0; d2 < dim; ++d2) - grad_grads[lexicographic_to_hierarchic[i + d * n_sub]][d] - [(d1 + d) % dim][(d2 + d) % dim] = - p_grad_grads[renumber_aniso[d][i]][d1][d2]; - - for (unsigned int i = 0; i < p_third_derivatives.size(); ++i) - for (unsigned int d1 = 0; d1 < dim; ++d1) - for (unsigned int d2 = 0; d2 < dim; ++d2) - for (unsigned int d3 = 0; d3 < dim; ++d3) - third_derivatives[lexicographic_to_hierarchic[i + d * n_sub]][d] - [(d1 + d) % dim][(d2 + d) % dim] - [(d3 + d) % dim] = - p_third_derivatives[renumber_aniso[d][i]][d1] - [d2][d3]; - - for (unsigned int i = 0; i < p_fourth_derivatives.size(); ++i) - for (unsigned int d1 = 0; d1 < dim; ++d1) - for (unsigned int d2 = 0; d2 < dim; ++d2) - for (unsigned int d3 = 0; d3 < dim; ++d3) - for (unsigned int d4 = 0; d4 < dim; ++d4) - fourth_derivatives[lexicographic_to_hierarchic[i + d * n_sub]] - [d][(d1 + d) % dim][(d2 + d) % dim] - [(d3 + d) % dim][(d4 + d) % dim] = - p_fourth_derivatives[renumber_aniso[d][i]] - [d1][d2][d3][d4]; - } -} - - - -template -std::string -PolynomialsRaviartThomas::name() const -{ - return "RaviartThomas"; -} - - - template unsigned int -PolynomialsRaviartThomas::n_polynomials(const unsigned int degree) +PolynomialsRaviartThomas::n_polynomials( + const unsigned int normal_degree, + const unsigned int tangential_degree) { - return n_polynomials(degree + 1, degree); + return PolynomialsVectorAnisotropic::n_polynomials(normal_degree, tangential_degree); } template unsigned int -PolynomialsRaviartThomas::n_polynomials( - const unsigned int normal_degree, - const unsigned int tangential_degree) +PolynomialsRaviartThomas::n_polynomials(const unsigned int degree) { - return dim * (normal_degree + 1) * - Utilities::pow(tangential_degree + 1, dim - 1); + return PolynomialsVectorAnisotropic::n_polynomials(degree + 1, degree); } @@ -325,37 +137,6 @@ PolynomialsRaviartThomas::clone() const -template -std::vector> -PolynomialsRaviartThomas::get_polynomial_support_points() const -{ - Assert(dim > 0 && dim <= 3, ExcImpossibleInDim(dim)); - const Quadrature<1> tangential( - tangential_degree == 0 ? - static_cast>(QMidpoint<1>()) : - static_cast>(QGaussLobatto<1>(tangential_degree + 1))); - const Quadrature<1> normal( - normal_degree == 0 ? - static_cast>(QMidpoint<1>()) : - static_cast>(QGaussLobatto<1>(normal_degree + 1))); - const QAnisotropic quad = - (dim == 1 ? QAnisotropic(normal) : - (dim == 2 ? QAnisotropic(normal, tangential) : - QAnisotropic(normal, tangential, tangential))); - - const unsigned int n_sub = polynomial_space.n(); - std::vector> points(dim * n_sub); - points.resize(n_polynomials(normal_degree, tangential_degree)); - for (unsigned int d = 0; d < dim; ++d) - for (unsigned int i = 0; i < n_sub; ++i) - for (unsigned int c = 0; c < dim; ++c) - points[lexicographic_to_hierarchic[i + d * n_sub]][(c + d) % dim] = - quad.point(renumber_aniso[d][i])[c]; - return points; -} - - - template class PolynomialsRaviartThomas<1>; template class PolynomialsRaviartThomas<2>; template class PolynomialsRaviartThomas<3>; -- 2.39.5