Merge pull request #14348 from drwells/n-face-permutations

dealii · Oct 13, 2022 · b2edb79 · b2edb79
2 parents 19e1464 + 34a0444
commit b2edb79
Show file tree

Hide file tree

Showing 6 changed files with 138 additions and 9 deletions.
diff --git a/include/deal.II/grid/reference_cell.h b/include/deal.II/grid/reference_cell.h
@@ -554,6 +554,14 @@ class ReferenceCell
   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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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]::