]> https://gitweb.dealii.org/ - dealii.git/commitdiff
unit_support_points and unit_face_support_points for wedges and pyramids 13693/head
authorPeter Munch <peterrmuench@gmail.com>
Sun, 8 May 2022 15:13:39 +0000 (17:13 +0200)
committerPeter Munch <peterrmuench@gmail.com>
Mon, 9 May 2022 21:33:06 +0000 (23:33 +0200)
include/deal.II/base/polynomials_wedge.h
source/base/polynomials_wedge.cc
source/fe/fe_pyramid_p.cc
source/fe/fe_wedge_p.cc
tests/bits/unit_support_points_01.cc [moved from tests/bits/unit_support_points.cc with 100% similarity]
tests/bits/unit_support_points_01.output [moved from tests/bits/unit_support_points.output with 100% similarity]
tests/bits/unit_support_points_03.cc [new file with mode: 0644]
tests/bits/unit_support_points_03.output [new file with mode: 0644]

index c29410fc0a67114450324c629bb3819fd4790188..dfd86334f767647a5f963f8d06b913f02a5d342b 100644 (file)
 
 #include <deal.II/base/config.h>
 
+#include <deal.II/base/ndarray.h>
 #include <deal.II/base/polynomials_barycentric.h>
 #include <deal.II/base/scalar_polynomials_base.h>
 
 DEAL_II_NAMESPACE_OPEN
 
+
+namespace internal
+{
+  /**
+   * Decompose the shape-function index of a linear wedge into an index
+   * to access the right shape function within the triangle and and within
+   * the line.
+   */
+  static const constexpr dealii::ndarray<unsigned int, 6, 2> wedge_table_1{
+    {{{0, 0}}, {{1, 0}}, {{2, 0}}, {{0, 1}}, {{1, 1}}, {{2, 1}}}};
+
+  /**
+   * Decompose the shape-function index of a quadratic wedge into an index
+   * to access the right shape function within the triangle and and within
+   * the line.
+   */
+  static const constexpr dealii::ndarray<unsigned int, 18, 2> wedge_table_2{
+    {{{0, 0}},
+     {{1, 0}},
+     {{2, 0}},
+     {{0, 1}},
+     {{1, 1}},
+     {{2, 1}},
+     {{3, 0}},
+     {{4, 0}},
+     {{5, 0}},
+     {{3, 1}},
+     {{4, 1}},
+     {{5, 1}},
+     {{0, 2}},
+     {{1, 2}},
+     {{2, 2}},
+     {{3, 2}},
+     {{4, 2}},
+     {{5, 2}}}};
+} // namespace internal
+
+
 /**
  * Polynomials defined on wedge entities. This class is basis of
  * FE_WedgeP.
index 402fbf954ea345520d464fc545cc4fbb8afae77a..186b228f1572d68943e145b5ca24e7890e3916dc 100644 (file)
@@ -14,7 +14,6 @@
 // ---------------------------------------------------------------------
 
 
-#include <deal.II/base/ndarray.h>
 #include <deal.II/base/polynomials_barycentric.h>
 #include <deal.II/base/polynomials_wedge.h>
 
@@ -50,49 +49,14 @@ ScalarLagrangePolynomialWedge<dim>::ScalarLagrangePolynomialWedge(
 {}
 
 
-namespace
-{
-  /**
-   * Decompose the shape-function index of a linear wedge into an index
-   * to access the right shape function within the triangle and and within
-   * the line.
-   */
-  static const constexpr ndarray<unsigned int, 6, 2> wedge_table_1{
-    {{{0, 0}}, {{1, 0}}, {{2, 0}}, {{0, 1}}, {{1, 1}}, {{2, 1}}}};
-
-  /**
-   * Decompose the shape-function index of a quadratic wedge into an index
-   * to access the right shape function within the triangle and and within
-   * the line.
-   */
-  static const constexpr ndarray<unsigned int, 18, 2> wedge_table_2{{{{0, 0}},
-                                                                     {{1, 0}},
-                                                                     {{2, 0}},
-                                                                     {{0, 1}},
-                                                                     {{1, 1}},
-                                                                     {{2, 1}},
-                                                                     {{3, 0}},
-                                                                     {{4, 0}},
-                                                                     {{5, 0}},
-                                                                     {{3, 1}},
-                                                                     {{4, 1}},
-                                                                     {{5, 1}},
-                                                                     {{0, 2}},
-                                                                     {{1, 2}},
-                                                                     {{2, 2}},
-                                                                     {{3, 2}},
-                                                                     {{4, 2}},
-                                                                     {{5, 2}}}};
-} // namespace
-
-
 
 template <int dim>
 double
 ScalarLagrangePolynomialWedge<dim>::compute_value(const unsigned int i,
                                                   const Point<dim> & p) const
 {
-  const auto pair = this->degree() == 1 ? wedge_table_1[i] : wedge_table_2[i];
+  const auto pair = this->degree() == 1 ? internal::wedge_table_1[i] :
+                                          internal::wedge_table_2[i];
 
   const Point<2> p_tri(p[0], p[1]);
   const auto     v_tri = poly_tri.compute_value(pair[0], p_tri);
@@ -110,7 +74,8 @@ Tensor<1, dim>
 ScalarLagrangePolynomialWedge<dim>::compute_grad(const unsigned int i,
                                                  const Point<dim> & p) const
 {
-  const auto pair = this->degree() == 1 ? wedge_table_1[i] : wedge_table_2[i];
+  const auto pair = this->degree() == 1 ? internal::wedge_table_1[i] :
+                                          internal::wedge_table_2[i];
 
   const Point<2> p_tri(p[0], p[1]);
   const auto     v_tri = poly_tri.compute_value(pair[0], p_tri);
index 60bd0d0a03df2247a6116fc53fce3a153dc35c78..2ef5822bfb47e35c5b758f14db7f55b87f030574 100644 (file)
@@ -97,9 +97,21 @@ FE_PyramidPoly<dim, spacedim>::FE_PyramidPoly(
 
   if (degree == 1)
     {
-      for (const unsigned int i : ReferenceCells::Pyramid.vertex_indices())
+      for (const auto i : this->reference_cell().vertex_indices())
         this->unit_support_points.emplace_back(
-          ReferenceCells::Pyramid.vertex<dim>(i));
+          this->reference_cell().template vertex<dim>(i));
+
+      this->unit_face_support_points.resize(this->reference_cell().n_faces());
+
+      for (const auto f : this->reference_cell().face_indices())
+        {
+          const auto face_reference_cell =
+            this->reference_cell().face_reference_cell(f);
+
+          for (const auto i : face_reference_cell.vertex_indices())
+            this->unit_face_support_points[f].emplace_back(
+              face_reference_cell.template vertex<dim - 1>(i));
+        }
     }
   else
     Assert(false, ExcNotImplemented());
index 42e6c1dcecc70f0c5d2d93708d0e37e8efa945a4..263699c606efc0b908eca6bf38333762dabaa865 100644 (file)
@@ -105,19 +105,36 @@ FE_WedgePoly<dim, spacedim>::FE_WedgePoly(
 {
   AssertDimension(dim, 3);
 
-  if (degree == 1)
-    {
-      for (const unsigned int i : ReferenceCells::Wedge.vertex_indices())
-        this->unit_support_points.emplace_back(
-          ReferenceCells::Wedge.vertex<dim>(i));
-    }
-  else
+  Assert(1 <= degree && degree <= 2, ExcNotImplemented());
+
+  FE_SimplexP<2> fe_triangle(degree);
+  FE_Q<1>        fe_line(degree);
+  FE_Q<2>        fe_quad(degree);
+
+  for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
     {
-      // TODO: Wedge elements work for higher degrees, but we don't currently
-      // fill their support points. Leaving the array empty is valid, however,
-      // and will simply result in an error when someone tries to access the
-      // array.
+      const auto pair = this->degree == 1 ? internal::wedge_table_1[i] :
+                                            internal::wedge_table_2[i];
+
+      this->unit_support_points.emplace_back(
+        fe_triangle.get_unit_support_points()[pair[0]][0],
+        fe_triangle.get_unit_support_points()[pair[0]][1],
+        fe_line.get_unit_support_points()[pair[1]][0]);
     }
+
+  this->unit_face_support_points.resize(this->reference_cell().n_faces());
+
+  for (const auto f : this->reference_cell().face_indices())
+    if (this->reference_cell().face_reference_cell(f) ==
+        ReferenceCells::Triangle)
+      for (const auto &p : fe_triangle.get_unit_support_points())
+        this->unit_face_support_points[f].emplace_back(p[0], p[1]);
+    else if (this->reference_cell().face_reference_cell(f) ==
+             ReferenceCells::Quadrilateral)
+      for (const auto &p : fe_quad.get_unit_support_points())
+        this->unit_face_support_points[f].emplace_back(p[0], p[1]);
+    else
+      Assert(false, ExcInternalError());
 }
 
 
diff --git a/tests/bits/unit_support_points_03.cc b/tests/bits/unit_support_points_03.cc
new file mode 100644 (file)
index 0000000..b630f5f
--- /dev/null
@@ -0,0 +1,96 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2022 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+// Test FiniteElement::unit_support_points() and ::unit_face_support_points()
+// for different shapes.
+
+#include <deal.II/fe/fe_pyramid_p.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_wedge_p.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <int dim, int spacedim>
+void
+test(const FiniteElement<dim, spacedim> &fe)
+{
+  deallog << fe.get_name() << std::endl;
+
+  const auto &unit_support_points = fe.get_unit_support_points();
+
+  AssertDimension(unit_support_points.size(), fe.n_dofs_per_cell());
+
+  for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
+    for (unsigned int j = 0; j < fe.n_dofs_per_cell(); ++j)
+      {
+        const auto value = fe.shape_value(i, unit_support_points[j]);
+        Assert((i == j ? (std::abs(value - 1.0) < 1e-8) :
+                         (std::abs(value) < 1e-8)),
+               ExcInternalError());
+      }
+
+  FE_Q<dim - 1>        fe_q(fe.degree);
+  FE_SimplexP<dim - 1> fe_p(fe.degree);
+
+  for (const auto f : fe.reference_cell().face_indices())
+    {
+      const auto &unit_face_support_points = fe.get_unit_face_support_points(f);
+
+      AssertDimension(unit_face_support_points.size(), fe.n_dofs_per_face(f));
+
+      const auto &fe_face =
+        fe.reference_cell().face_reference_cell(f).is_hyper_cube() ?
+          static_cast<FiniteElement<dim - 1> &>(fe_q) :
+          static_cast<FiniteElement<dim - 1> &>(fe_p);
+
+      for (unsigned int i = 0; i < fe.n_dofs_per_face(f); ++i)
+        for (unsigned int j = 0; j < fe.n_dofs_per_face(f); ++j)
+          {
+            const auto value =
+              fe_face.shape_value(i, unit_face_support_points[j]);
+            Assert((i == j ? (std::abs(value - 1.0) < 1e-8) :
+                             (std::abs(value) < 1e-8)),
+                   ExcInternalError());
+          }
+    }
+
+
+
+  deallog << "OK!" << std::endl;
+}
+
+int
+main()
+{
+  initlog();
+
+  test(FE_Q<2>(1));
+  test(FE_Q<2>(2));
+  test(FE_Q<3>(1));
+  test(FE_Q<3>(2));
+
+  test(FE_SimplexP<2>(1));
+  test(FE_SimplexP<2>(2));
+  test(FE_SimplexP<3>(1));
+  test(FE_SimplexP<3>(2));
+
+  test(FE_PyramidP<3>(1));
+
+  test(FE_WedgeP<3>(1));
+  test(FE_WedgeP<3>(2));
+}
diff --git a/tests/bits/unit_support_points_03.output b/tests/bits/unit_support_points_03.output
new file mode 100644 (file)
index 0000000..8c41550
--- /dev/null
@@ -0,0 +1,23 @@
+
+DEAL::FE_Q<2>(1)
+DEAL::OK!
+DEAL::FE_Q<2>(2)
+DEAL::OK!
+DEAL::FE_Q<3>(1)
+DEAL::OK!
+DEAL::FE_Q<3>(2)
+DEAL::OK!
+DEAL::FE_SimplexP<2>(1)
+DEAL::OK!
+DEAL::FE_SimplexP<2>(2)
+DEAL::OK!
+DEAL::FE_SimplexP<3>(1)
+DEAL::OK!
+DEAL::FE_SimplexP<3>(2)
+DEAL::OK!
+DEAL::FE_PyramidP<3>(1)
+DEAL::OK!
+DEAL::FE_WedgeP<3>(1)
+DEAL::OK!
+DEAL::FE_WedgeP<3>(2)
+DEAL::OK!

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.