Implement ReferenceCell::n_face_orientations().

This encodes some previously hard-coded values in FESystem.

diff --git a/include/deal.II/grid/reference_cell.h b/include/deal.II/grid/reference_cell.h
index 5fa19e26fa1f7a5c1e82f0ea5e11a54fec94d49d..0866e39ab8f5a905bee1aaa12366c97ed542ac8c 100644 (file)
--- a/include/deal.II/grid/reference_cell.h
+++ b/include/deal.II/grid/reference_cell.h
@@ -554,6 +554,14 @@ public:
   Tensor<1, dim>
   unit_normal_vectors(const unsigned int face_no) const;
 
+  /**
+   * Return the number of orientations for a face in the ReferenceCell. For
+   * example, for hexahedra this is 8 for every face since quadrilaterals have
+   * 8 possible orientations.
+   */
+  unsigned int
+  n_face_orientations(const unsigned int face_no) const;
+
   /**
    * Determine the orientation of the current entity described by its
    * vertices @p var_1 relative to an entity described by @p var_0.
@@ -2235,6 +2243,25 @@ ReferenceCell::unit_normal_vectors(const unsigned int face_no) const
 
 
 
+inline unsigned int
+ReferenceCell::n_face_orientations(const unsigned int face_no) const
+{
+  AssertIndexRange(face_no, n_faces());
+  if (get_dimension() == 1)
+    return 1;
+  if (get_dimension() == 2)
+    return 2;
+  else if (face_reference_cell(face_no) == ReferenceCells::Quadrilateral)
+    return 8;
+  else if (face_reference_cell(face_no) == ReferenceCells::Triangle)
+    return 6;
+
+  Assert(false, ExcInternalError());
+  return numbers::invalid_unsigned_int;
+}
+
+
+
 inline bool
 ReferenceCell::standard_vs_true_line_orientation(
   const unsigned int  line,

diff --git a/source/fe/fe.cc b/source/fe/fe.cc
index 857c3fac9992033e78a024c6c6d2f4d1eb0bb011..bf93b102d14789754215d78860d8becdf99bbdc8 100644 (file)
--- a/source/fe/fe.cc
+++ b/source/fe/fe.cc
@@ -158,10 +158,7 @@ FiniteElement<dim, spacedim>::FiniteElement(
         {
           adjust_quad_dof_index_for_face_orientation_table[f] =
             Table<2, int>(this->n_dofs_per_quad(f),
-                          this->reference_cell().face_reference_cell(f) ==
-                              ReferenceCells::Quadrilateral ?
-                            8 :
-                            6);
+                          this->reference_cell().n_face_orientations(f));
           adjust_quad_dof_index_for_face_orientation_table[f].fill(0);
         }
     }

diff --git a/source/fe/fe_raviart_thomas.cc b/source/fe/fe_raviart_thomas.cc
index 029c2ae1284aef26ad15574fe08a714d54200b22..1978e07a06a08130e9f1a6a67c36d572729203b6 100644 (file)
--- a/source/fe/fe_raviart_thomas.cc
+++ b/source/fe/fe_raviart_thomas.cc
@@ -269,9 +269,11 @@ FE_RaviartThomas<dim>::initialize_quad_dof_index_permutation_and_sign_change()
   AssertDimension(this->n_unique_faces(), 1);
   const unsigned int face_no = 0;
 
-  Assert(this->adjust_quad_dof_index_for_face_orientation_table[0]
-             .n_elements() == 8 * this->n_dofs_per_quad(face_no),
-         ExcInternalError());
+  Assert(
+    this->adjust_quad_dof_index_for_face_orientation_table[0].n_elements() ==
+      this->reference_cell().n_face_orientations(face_no) *
+        this->n_dofs_per_quad(face_no),
+    ExcInternalError());
 
   // The 3D RaviartThomas space has tensor_degree*tensor_degree face dofs
   const unsigned int n = this->tensor_degree();

diff --git a/source/fe/fe_system.cc b/source/fe/fe_system.cc
index 84417923edbef2bfab463afd1db96d8ef9affaca..168613170d10831f7f3412a2cbdac11c15173a84 100644 (file)
--- a/source/fe/fe_system.cc
+++ b/source/fe/fe_system.cc
@@ -2008,10 +2008,19 @@ FESystem<dim, spacedim>::initialize(
       for (unsigned int face_no = 0; face_no < this->n_unique_faces();
            ++face_no)
         {
+          // This logic isn't valid for higher-order wedges because the first
+          // two faces are triangles, so no table will be set up for
+          // quadrilateral faces
+          Assert(this->n_dofs_per_quad(face_no) == 0 ||
+                   this->reference_cell() != ReferenceCells::Wedge,
+                 ExcNotImplemented());
+
           // the array into which we want to write should have the correct size
           // already.
           Assert(this->adjust_quad_dof_index_for_face_orientation_table[face_no]
-                     .n_elements() == 8 * this->n_dofs_per_quad(face_no),
+                     .n_elements() ==
+                   this->reference_cell().n_face_orientations(face_no) *
+                     this->n_dofs_per_quad(face_no),
                  ExcInternalError());
 
           // to obtain the shifts for this composed element, copy the shift
@@ -2025,7 +2034,10 @@ FESystem<dim, spacedim>::initialize(
               for (unsigned int c = 0; c < this->element_multiplicity(b); ++c)
                 {
                   for (unsigned int i = 0; i < temp.size(0); ++i)
-                    for (unsigned int j = 0; j < 8; ++j)
+                    for (unsigned int j = 0;
+                         j <
+                         this->reference_cell().n_face_orientations(face_no);
+                         ++j)
                       this->adjust_quad_dof_index_for_face_orientation_table
                         [face_no](index + i, j) = temp(i, j);
                   index += temp.size(0);

diff --git a/tests/simplex/fe_p_bubbles_03.cc b/tests/simplex/fe_p_bubbles_03.cc
new file mode 100644 (file)
index 0000000..08e5e51
--- /dev/null
+++ b/tests/simplex/fe_p_bubbles_03.cc
@@ -0,0 +1,79 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2022 - 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.
+//
+// ---------------------------------------------------------------------
+
+// Ensure that FE_SimplexP_Bubbles works with FESystem at higher-order. We
+// previously hardcoded some assumptions which prevented non-hypercube
+// elements with more than zero DoF per quad from working with FESystem.
+
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
+#include <deal.II/fe/fe_system.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <int dim, int spacedim>
+void
+test(const FiniteElement<dim, spacedim> &fe, const Quadrature<dim> &quad)
+{
+  deallog.push(fe.get_name());
+
+  for (const auto &point : quad.get_points())
+    {
+      deallog << point << " : " << std::endl;
+      for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
+        deallog << fe.shape_value(i, point) << ' ';
+      for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
+        deallog << fe.shape_grad(i, point) << ' ';
+      for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
+        deallog << fe.shape_grad_grad(i, point) << ' ';
+      deallog << std::endl;
+    }
+  deallog << std::endl;
+
+  deallog.pop();
+}
+
+
+
+template <int dim, int spacedim = dim>
+void
+test_unit_support_points()
+{
+  deallog << "Test support points for dim = " << dim
+          << " and spacedim = " << spacedim << std::endl;
+  for (unsigned int degree = 1; degree < 3; ++degree)
+    {
+      deallog << "approximation degree = " << degree << std::endl;
+      FESystem<dim, spacedim> fe(FE_SimplexP_Bubbles<dim, spacedim>(degree),
+                                 dim);
+      deallog << "element tensor degree = " << fe.tensor_degree() << std::endl;
+      Quadrature<dim> quad(
+        fe.reference_cell().template get_midpoint_quadrature<dim>());
+      test(fe, quad);
+    }
+}
+
+
+
+int
+main()
+{
+  initlog();
+
+  test_unit_support_points<3>();
+}

diff --git a/tests/simplex/fe_p_bubbles_03.output b/tests/simplex/fe_p_bubbles_03.output
new file mode 100644 (file)
index 0000000..fc028be
--- /dev/null
+++ b/tests/simplex/fe_p_bubbles_03.output
@@ -0,0 +1,12 @@
+
+DEAL::Test support points for dim = 3 and spacedim = 3
+DEAL::approximation degree = 1
+DEAL::element tensor degree = 1
+DEAL:FESystem<3>[FE_SimplexP_Bubbles<3>(1)^3]::0.250000 0.250000 0.250000 : 
+DEAL:FESystem<3>[FE_SimplexP_Bubbles<3>(1)^3]::0.250000 0.250000 0.250000 0.250000 0.250000 0.250000 0.250000 0.250000 0.250000 0.250000 0.250000 0.250000 -1.00000 -1.00000 -1.00000 -1.00000 -1.00000 -1.00000 -1.00000 -1.00000 -1.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 
+DEAL:FESystem<3>[FE_SimplexP_Bubbles<3>(1)^3]::
+DEAL::approximation degree = 2
+DEAL::element tensor degree = 3
+DEAL:FESystem<3>[FE_SimplexP_Bubbles<3>(2)^3]::0.250000 0.250000 0.250000 : 
+DEAL:FESystem<3>[FE_SimplexP_Bubbles<3>(2)^3]::0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 1.00000 1.00000 1.00000 -0.187500 -0.187500 -0.187500 -0.187500 -0.187500 -0.187500 -0.187500 -0.187500 -0.187500 0.187500 0.00000 0.00000 0.187500 0.00000 0.00000 0.187500 0.00000 0.00000 0.00000 0.187500 0.00000 0.00000 0.187500 0.00000 0.00000 0.187500 0.00000 0.00000 0.00000 0.187500 0.00000 0.00000 0.187500 0.00000 0.00000 0.187500 -4.44089e-16 -0.250000 -0.250000 -4.44089e-16 -0.250000 -0.250000 -4.44089e-16 -0.250000 -0.250000 0.250000 0.250000 0.00000 0.250000 0.250000 0.00000 0.250000 0.250000 0.00000 -0.250000 -4.44089e-16 -0.250000 -0.250000 -4.44089e-16 -0.250000 -0.250000 -4.44089e-16 -0.250000 -0.250000 -0.250000 -4.44089e-16 -0.250000 -0.250000 -4.44089e-16 -0.250000 -0.250000 -4.44089e-16 0.250000 0.00000 0.250000 0.250000 0.00000 0.250000 0.250000 0.00000 0.250000 0.00000 0.250000 0.250000 0.00000 0.250000 0.250000 0.00000 0.250000 0.250000 0.00000 0.00000 -1.68750 0.00000 0.00000 -1.68750 0.00000 0.00000 -1.68750 0.00000 -1.68750 0.00000 0.00000 -1.68750 0.00000 0.00000 -1.68750 0.00000 -1.68750 0.00000 0.00000 -1.68750 0.00000 0.00000 -1.68750 0.00000 0.00000 1.68750 1.68750 1.68750 1.68750 1.68750 1.68750 1.68750 1.68750 1.68750 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 1.50000 2.00000 2.00000 2.00000 1.50000 2.00000 2.00000 2.00000 1.50000 1.50000 2.00000 2.00000 2.00000 1.50000 2.00000 2.00000 2.00000 1.50000 1.50000 2.00000 2.00000 2.00000 1.50000 2.00000 2.00000 2.00000 1.50000 1.50000 -0.500000 -0.500000 -0.500000 -1.00000 -0.500000 -0.500000 -0.500000 -1.00000 1.50000 -0.500000 -0.500000 -0.500000 -1.00000 -0.500000 -0.500000 -0.500000 -1.00000 1.50000 -0.500000 -0.500000 -0.500000 -1.00000 -0.500000 -0.500000 -0.500000 -1.00000 -1.00000 -0.500000 -0.500000 -0.500000 1.50000 -0.500000 -0.500000 -0.500000 -1.00000 -1.00000 -0.500000 -0.500000 -0.500000 1.50000 -0.500000 -0.500000 -0.500000 -1.00000 -1.00000 -0.500000 -0.500000 -0.500000 1.50000 -0.500000 -0.500000 -0.500000 -1.00000 -1.00000 -0.500000 -0.500000 -0.500000 -1.00000 -0.500000 -0.500000 -0.500000 1.50000 -1.00000 -0.500000 -0.500000 -0.500000 -1.00000 -0.500000 -0.500000 -0.500000 1.50000 -1.00000 -0.500000 -0.500000 -0.500000 -1.00000 -0.500000 -0.500000 -0.500000 1.50000 -3.55271e-15 -3.55271e-15 -3.55271e-15 -3.55271e-15 2.00000 4.00000 -3.55271e-15 4.00000 2.00000 -3.55271e-15 -3.55271e-15 -3.55271e-15 -3.55271e-15 2.00000 4.00000 -3.55271e-15 4.00000 2.00000 -3.55271e-15 -3.55271e-15 -3.55271e-15 -3.55271e-15 2.00000 4.00000 -3.55271e-15 4.00000 2.00000 2.00000 2.00000 -2.00000 2.00000 2.00000 -2.00000 -2.00000 -2.00000 -4.00000 2.00000 2.00000 -2.00000 2.00000 2.00000 -2.00000 -2.00000 -2.00000 -4.00000 2.00000 2.00000 -2.00000 2.00000 2.00000 -2.00000 -2.00000 -2.00000 -4.00000 2.00000 -3.55271e-15 4.00000 -3.55271e-15 -3.55271e-15 -3.10862e-15 4.00000 -3.10862e-15 2.00000 2.00000 -3.55271e-15 4.00000 -3.55271e-15 -3.55271e-15 -3.10862e-15 4.00000 -3.10862e-15 2.00000 2.00000 -3.55271e-15 4.00000 -3.55271e-15 -3.55271e-15 -3.10862e-15 4.00000 -3.10862e-15 2.00000 2.00000 4.00000 -3.55271e-15 4.00000 2.00000 -3.10862e-15 -3.55271e-15 -3.10862e-15 -3.55271e-15 2.00000 4.00000 -3.55271e-15 4.00000 2.00000 -3.10862e-15 -3.55271e-15 -3.10862e-15 -3.55271e-15 2.00000 4.00000 -3.55271e-15 4.00000 2.00000 -3.10862e-15 -3.55271e-15 -3.10862e-15 -3.55271e-15 2.00000 -2.00000 2.00000 -2.00000 -4.00000 -2.00000 2.00000 -2.00000 2.00000 2.00000 -2.00000 2.00000 -2.00000 -4.00000 -2.00000 2.00000 -2.00000 2.00000 2.00000 -2.00000 2.00000 -2.00000 -4.00000 -2.00000 2.00000 -2.00000 2.00000 -4.00000 -2.00000 -2.00000 -2.00000 2.00000 2.00000 -2.00000 2.00000 2.00000 -4.00000 -2.00000 -2.00000 -2.00000 2.00000 2.00000 -2.00000 2.00000 2.00000 -4.00000 -2.00000 -2.00000 -2.00000 2.00000 2.00000 -2.00000 2.00000 2.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 13.5000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 13.5000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 13.5000 0.00000 0.00000 0.00000 0.00000 13.5000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 13.5000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 13.5000 0.00000 0.00000 0.00000 0.00000 13.5000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 13.5000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 13.5000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 13.5000 13.5000 13.5000 13.5000 13.5000 13.5000 13.5000 13.5000 13.5000 13.5000 13.5000 13.5000 13.5000 13.5000 13.5000 13.5000 13.5000 13.5000 13.5000 13.5000 13.5000 13.5000 13.5000 13.5000 13.5000 13.5000 13.5000 -32.0000 -16.0000 -16.0000 -16.0000 -32.0000 -16.0000 -16.0000 -16.0000 -32.0000 -32.0000 -16.0000 -16.0000 -16.0000 -32.0000 -16.0000 -16.0000 -16.0000 -32.0000 -32.0000 -16.0000 -16.0000 -16.0000 -32.0000 -16.0000 -16.0000 -16.0000 -32.0000 
+DEAL:FESystem<3>[FE_SimplexP_Bubbles<3>(2)^3]::