From: David Wells <drwells@email.unc.edu>
Date: Wed, 12 Oct 2022 19:09:33 +0000 (-0400)
Subject: Implement ReferenceCell::n_face_orientations().
X-Git-Tag: v9.5.0-rc1~898^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=34a04440076a298bfdabfec9ee70182f7646b41a;p=dealii.git

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 5fa19e26fa..0866e39ab8 100644
--- 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 857c3fac99..bf93b102d1 100644
--- 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 029c2ae128..1978e07a06 100644
--- 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 84417923ed..168613170d 100644
--- 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
index 0000000000..08e5e513c5
--- /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
index 0000000000..fc028bede4
--- /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]::