]> https://gitweb.dealii.org/ - dealii.git/commitdiff
adapted FE_RT to polynomials_vector_anisotropic
authorNatalia Nebulishvili <natalia.nebulishvili@rub.de>
Fri, 20 Jun 2025 10:19:20 +0000 (12:19 +0200)
committerNatalia Nebulishvili <natalia.nebulishvili@rub.de>
Sun, 22 Jun 2025 13:42:14 +0000 (15:42 +0200)
include/deal.II/fe/fe_raviart_thomas.h
include/deal.II/matrix_free/shape_info.templates.h
source/fe/fe_raviart_thomas.cc
source/fe/fe_raviart_thomas_nodal.cc

index 5f0464fb696fd0b8ee8d325d6a168497836c0986..46541cdd61e99a8bb06d4ba4adc2ed99eee7c0cb 100644 (file)
@@ -144,6 +144,14 @@ public:
   virtual std::unique_ptr<FiniteElement<dim, dim>>
   clone() const override;
 
+  /**
+   * Compute the lexicographic to hierarchic numbering underlying this class,
+   * necessary for the creation of the respective vector polynomial space.
+   */
+  std::vector<unsigned int>
+  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<FiniteElement<dim, dim>>
   clone() const override;
 
+  /**
+   * Compute the lexicographic to hierarchic numbering underlying this class,
+   * necessary for the creation of the respective vector polynomial space.
+   */
+  std::vector<unsigned int>
+  get_lexicographic_numbering(const unsigned int normal_degree,
+                              const unsigned int tangential_degree) const;
+
   virtual void
   get_face_interpolation_matrix(const FiniteElement<dim> &source,
                                 FullMatrix<double>       &matrix,
index f3fc16f505b29a1c918c8bf91808d94df5dc78bd..8c852e76ffee351bce22a6f38a2985991699bc3e 100644 (file)
@@ -21,7 +21,6 @@
 #include <deal.II/base/memory_consumption.h>
 #include <deal.II/base/polynomial.h>
 #include <deal.II/base/polynomials_piecewise.h>
-#include <deal.II/base/polynomials_raviart_thomas.h>
 #include <deal.II/base/qprojector.h>
 #include <deal.II/base/tensor_product_polynomials.h>
 #include <deal.II/base/utilities.h>
index 7308a8ea93c5e4789a96161cddb369512e0669eb..cadd78a4ac67bdd8b746af88239ef41642a08b94 100644 (file)
@@ -14,7 +14,7 @@
 
 
 #include <deal.II/base/polynomial.h>
-#include <deal.II/base/polynomials_raviart_thomas.h>
+#include <deal.II/base/polynomials_vector_anisotropic.h>
 #include <deal.II/base/qprojector.h>
 #include <deal.II/base/quadrature.h>
 #include <deal.II/base/quadrature_lib.h>
@@ -43,16 +43,20 @@ DEAL_II_NAMESPACE_OPEN
 template <int dim>
 FE_RaviartThomas<dim>::FE_RaviartThomas(const unsigned int deg)
   : FE_PolyTensor<dim>(
-      PolynomialsRaviartThomas<dim>(deg),
+      PolynomialsVectorAnisotropic<dim>(deg + 1,
+                                        deg,
+                                        get_lexicographic_numbering(deg + 1,
+                                                                    deg)),
       FiniteElementData<dim>(get_dpo_vector(deg),
                              dim,
                              deg + 1,
                              FiniteElementData<dim>::Hdiv),
-      std::vector<bool>(PolynomialsRaviartThomas<dim>::n_polynomials(deg),
-                        true),
-      std::vector<ComponentMask>(PolynomialsRaviartThomas<dim>::n_polynomials(
-                                   deg),
-                                 ComponentMask(std::vector<bool>(dim, true))))
+      std::vector<bool>(
+        PolynomialsVectorAnisotropic<dim>::n_polynomials(deg + 1, deg),
+        true),
+      std::vector<ComponentMask>(
+        PolynomialsVectorAnisotropic<dim>::n_polynomials(deg + 1, deg),
+        ComponentMask(std::vector<bool>(dim, true))))
 {
   Assert(dim >= 2, ExcImpossibleInDim(dim));
   const unsigned int n_dofs = this->n_dofs_per_cell();
@@ -694,6 +698,71 @@ FE_RaviartThomas<dim>::has_support_on_face(const unsigned int shape_index,
 
 
 
+template <int dim>
+std::vector<unsigned int>
+FE_RaviartThomas<dim>::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<unsigned int> 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 <int dim>
 void
 FE_RaviartThomas<dim>::convert_generalized_support_point_values_to_dof_values(
index 3b61ab611ae8f2571458f606df80aa4918b7f96c..d25986923628c8671e22b72ccdccd3979f217ba0 100644 (file)
@@ -14,7 +14,7 @@
 
 
 #include <deal.II/base/polynomial.h>
-#include <deal.II/base/polynomials_raviart_thomas.h>
+#include <deal.II/base/polynomials_vector_anisotropic.h>
 #include <deal.II/base/qprojector.h>
 #include <deal.II/base/quadrature_lib.h>
 
@@ -66,24 +66,32 @@ namespace
 
 template <int dim>
 FE_RaviartThomasNodal<dim>::FE_RaviartThomasNodal(const unsigned int degree)
-  : FE_PolyTensor<dim>(PolynomialsRaviartThomas<dim>(degree),
-                       FiniteElementData<dim>(get_rt_dpo_vector(dim, degree),
-                                              dim,
-                                              degree + 1,
-                                              FiniteElementData<dim>::Hdiv),
-                       std::vector<bool>(1, false),
-                       std::vector<ComponentMask>(
-                         PolynomialsRaviartThomas<dim>::n_polynomials(degree),
-                         ComponentMask(std::vector<bool>(dim, true))))
+  : FE_PolyTensor<dim>(
+      PolynomialsVectorAnisotropic<dim>(degree + 1,
+                                        degree,
+                                        get_lexicographic_numbering(degree + 1,
+                                                                    degree)),
+      FiniteElementData<dim>(get_rt_dpo_vector(dim, degree),
+                             dim,
+                             degree + 1,
+                             FiniteElementData<dim>::Hdiv),
+      std::vector<bool>(1, false),
+      std::vector<ComponentMask>(
+        PolynomialsVectorAnisotropic<dim>::n_polynomials(degree + 1, degree),
+        ComponentMask(std::vector<bool>(dim, true))))
 {
   Assert(dim >= 2, ExcImpossibleInDim(dim));
 
   this->mapping_kind = {mapping_raviart_thomas};
 
+  const std::vector<unsigned int> 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<dim>(degree).get_polynomial_support_points();
+    PolynomialsVectorAnisotropic<dim>(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<dim>::has_support_on_face(
 
 
 
+template <int dim>
+std::vector<unsigned int>
+FE_RaviartThomasNodal<dim>::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<unsigned int> 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 <int dim>
 void
 FE_RaviartThomasNodal<dim>::

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.