]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Removed most of the duplicated code from Raviart-Thomas polynomials
authorNatalia Nebulishvili <natalia.nebulishvili@rub.de>
Tue, 17 Jun 2025 15:45:47 +0000 (17:45 +0200)
committerNatalia Nebulishvili <natalia.nebulishvili@rub.de>
Tue, 17 Jun 2025 15:45:47 +0000 (17:45 +0200)
include/deal.II/base/polynomials_raviart_thomas.h
source/base/polynomials_raviart_thomas.cc

index c93fd827f505a691b2cf38e496b919a521c93601..b7a618a6ce8439f59cad8aea72e7a573e527848d 100644 (file)
@@ -23,7 +23,7 @@
 #include <deal.II/base/polynomial.h>
 #include <deal.II/base/polynomial_space.h>
 #include <deal.II/base/tensor.h>
-#include <deal.II/base/tensor_polynomials_base.h>
+#include <deal.II/base/polynomials_vector_anisotropic.h>
 #include <deal.II/base/tensor_product_polynomials.h>
 
 #include <mutex>
@@ -46,7 +46,7 @@ DEAL_II_NAMESPACE_OPEN
  * @ingroup Polynomials
  */
 template <int dim>
-class PolynomialsRaviartThomas : public TensorPolynomialsBase<dim>
+class PolynomialsRaviartThomas : public PolynomialsVectorAnisotropic<dim>
 {
 public:
   /**
@@ -80,19 +80,6 @@ public:
    * The size of the vectors must either be zero or equal <tt>n()</tt>. In
    * the first case, the function will not compute these values.
    */
-  void
-  evaluate(const Point<dim>            &unit_point,
-           std::vector<Tensor<1, dim>> &values,
-           std::vector<Tensor<2, dim>> &grads,
-           std::vector<Tensor<3, dim>> &grad_grads,
-           std::vector<Tensor<4, dim>> &third_derivatives,
-           std::vector<Tensor<5, dim>> &fourth_derivatives) const override;
-
-  /**
-   * Return the name of the space, which is <tt>PolynomialsRaviartThomas</tt>.
-   */
-  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<TensorPolynomialsBase<dim>>
   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<Point<dim>>
-  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<dim> polynomial_space;
-
-  /**
-   * Renumbering from lexicographic to hierarchic order.
-   */
-  std::vector<unsigned int> lexicographic_to_hierarchic;
-
-  /**
-   * Renumbering from hierarchic to lexicographic order. Inverse of
-   * lexicographic_to_hierarchic.
-   */
-  std::vector<unsigned int> hierarchic_to_lexicographic;
-
-  /**
-   * Renumbering from shifted polynomial spaces to lexicographic one.
-   */
-  std::array<std::vector<unsigned int>, dim> renumber_aniso;
 };
 
 
index f9a50d266bbb0b8154fa46c68dab811a60bfd9ea..74ffe95e71c5ad9c82c1fe3e1b718987d7f1bcbe 100644 (file)
 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<std::vector<Polynomials::Polynomial<double>>>
-  create_rt_polynomials(const unsigned int dim,
-                        const unsigned int normal_degree,
-                        const unsigned int tangential_degree)
-  {
-    std::vector<std::vector<Polynomials::Polynomial<double>>> 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 <int dim>
 PolynomialsRaviartThomas<dim>::PolynomialsRaviartThomas(
   const unsigned int normal_degree,
   const unsigned int tangential_degree)
-  : TensorPolynomialsBase<dim>(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<dim>(normal_degree, tangential_degree,
+                               get_lexicographic_numbering(normal_degree, tangential_degree))
+{}
 
 
 
@@ -128,125 +43,22 @@ PolynomialsRaviartThomas<dim>::PolynomialsRaviartThomas(const unsigned int k)
 
 
 
-template <int dim>
-void
-PolynomialsRaviartThomas<dim>::evaluate(
-  const Point<dim>            &unit_point,
-  std::vector<Tensor<1, dim>> &values,
-  std::vector<Tensor<2, dim>> &grads,
-  std::vector<Tensor<3, dim>> &grad_grads,
-  std::vector<Tensor<4, dim>> &third_derivatives,
-  std::vector<Tensor<5, dim>> &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<double>         p_values;
-  std::vector<Tensor<1, dim>> p_grads;
-  std::vector<Tensor<2, dim>> p_grad_grads;
-  std::vector<Tensor<3, dim>> p_third_derivatives;
-  std::vector<Tensor<4, dim>> 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<dim> 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 <int dim>
-std::string
-PolynomialsRaviartThomas<dim>::name() const
-{
-  return "RaviartThomas";
-}
-
-
-
 template <int dim>
 unsigned int
-PolynomialsRaviartThomas<dim>::n_polynomials(const unsigned int degree)
+PolynomialsRaviartThomas<dim>::n_polynomials(
+  const unsigned int normal_degree,
+  const unsigned int tangential_degree)
 {
-  return n_polynomials(degree + 1, degree);
+  return PolynomialsVectorAnisotropic<dim>::n_polynomials(normal_degree, tangential_degree);
 }
 
 
 
 template <int dim>
 unsigned int
-PolynomialsRaviartThomas<dim>::n_polynomials(
-  const unsigned int normal_degree,
-  const unsigned int tangential_degree)
+PolynomialsRaviartThomas<dim>::n_polynomials(const unsigned int degree)
 {
-  return dim * (normal_degree + 1) *
-         Utilities::pow(tangential_degree + 1, dim - 1);
+  return PolynomialsVectorAnisotropic<dim>::n_polynomials(degree + 1, degree);
 }
 
 
@@ -325,37 +137,6 @@ PolynomialsRaviartThomas<dim>::clone() const
 
 
 
-template <int dim>
-std::vector<Point<dim>>
-PolynomialsRaviartThomas<dim>::get_polynomial_support_points() const
-{
-  Assert(dim > 0 && dim <= 3, ExcImpossibleInDim(dim));
-  const Quadrature<1> tangential(
-    tangential_degree == 0 ?
-      static_cast<Quadrature<1>>(QMidpoint<1>()) :
-      static_cast<Quadrature<1>>(QGaussLobatto<1>(tangential_degree + 1)));
-  const Quadrature<1> normal(
-    normal_degree == 0 ?
-      static_cast<Quadrature<1>>(QMidpoint<1>()) :
-      static_cast<Quadrature<1>>(QGaussLobatto<1>(normal_degree + 1)));
-  const QAnisotropic<dim> quad =
-    (dim == 1 ? QAnisotropic<dim>(normal) :
-                (dim == 2 ? QAnisotropic<dim>(normal, tangential) :
-                            QAnisotropic<dim>(normal, tangential, tangential)));
-
-  const unsigned int      n_sub = polynomial_space.n();
-  std::vector<Point<dim>> 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>;

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.