Un-hardcode some orientation info in FE_SimplexP{_Bubbles}.

diff --git a/source/fe/fe_simplex_p.cc b/source/fe/fe_simplex_p.cc
index 18ab94784dcf3cc926d42f0e7b948c1b64214eef..e5998c5daf921d81a75017dfbececa6a832a51d9 100644 (file)
--- a/source/fe/fe_simplex_p.cc
+++ b/source/fe/fe_simplex_p.cc
@@ -91,71 +91,28 @@ namespace
         return unit_points;
       }
 
+    const auto reference_cell = ReferenceCells::get_simplex<dim>();
     // Piecewise constants are a special case: use a support point at the
     // centroid and only the centroid
     if (degree == 0)
       {
-        unit_points.emplace_back(
-          ReferenceCells::get_simplex<dim>().template barycenter<dim>());
+        unit_points.emplace_back(reference_cell.template barycenter<dim>());
         return unit_points;
       }
 
-    if (dim == 1)
-      {
-        // We don't really have dim = 1 support for simplex elements yet, but
-        // its convenient for populating the face array
-        Assert(degree <= 2, ExcNotImplemented());
-        if (degree >= 1)
-          {
-            unit_points.emplace_back(0.0);
-            unit_points.emplace_back(1.0);
+    Assert(degree <= 2, ExcNotImplemented());
+    for (const auto &vertex_no : reference_cell.vertex_indices())
+      unit_points.emplace_back(reference_cell.template vertex<dim>(vertex_no));
 
-            if (degree == 2)
-              unit_points.emplace_back(0.5);
-          }
-      }
-    else if (dim == 2)
-      {
-        Assert(degree <= 2, ExcNotImplemented());
-        if (degree >= 1)
-          {
-            unit_points.emplace_back(0.0, 0.0);
-            unit_points.emplace_back(1.0, 0.0);
-            unit_points.emplace_back(0.0, 1.0);
-
-            if (degree == 2)
-              {
-                unit_points.emplace_back(0.5, 0.0);
-                unit_points.emplace_back(0.5, 0.5);
-                unit_points.emplace_back(0.0, 0.5);
-              }
-          }
-      }
-    else if (dim == 3)
-      {
-        Assert(degree <= 2, ExcNotImplemented());
-        if (degree >= 1)
-          {
-            unit_points.emplace_back(0.0, 0.0, 0.0);
-            unit_points.emplace_back(1.0, 0.0, 0.0);
-            unit_points.emplace_back(0.0, 1.0, 0.0);
-            unit_points.emplace_back(0.0, 0.0, 1.0);
-
-            if (degree == 2)
-              {
-                unit_points.emplace_back(0.5, 0.0, 0.0);
-                unit_points.emplace_back(0.5, 0.5, 0.0);
-                unit_points.emplace_back(0.0, 0.5, 0.0);
-                unit_points.emplace_back(0.0, 0.0, 0.5);
-                unit_points.emplace_back(0.5, 0.0, 0.5);
-                unit_points.emplace_back(0.0, 0.5, 0.5);
-              }
-          }
-      }
-    else
-      {
-        Assert(false, ExcNotImplemented());
-      }
+    if (degree == 2)
+      for (const auto &line_no : reference_cell.line_indices())
+        {
+          const auto v0 = reference_cell.template vertex<dim>(
+            reference_cell.line_to_cell_vertices(line_no, 0));
+          const auto v1 = reference_cell.template vertex<dim>(
+            reference_cell.line_to_cell_vertices(line_no, 1));
+          unit_points.emplace_back((v0 + v1) / 2.0);
+        }
 
     return unit_points;
   }

diff --git a/source/fe/fe_simplex_p_bubbles.cc b/source/fe/fe_simplex_p_bubbles.cc
index 260e5a384737f339bbc488a1f2b367be2b11c436..d52d4ed11eeb8099ead9fb351618cfd22622bf82 100644 (file)
--- a/source/fe/fe_simplex_p_bubbles.cc
+++ b/source/fe/fe_simplex_p_bubbles.cc
@@ -67,8 +67,8 @@ namespace FE_P_BubblesImplementation
     FE_SimplexP<dim>        fe_p(degree);
     std::vector<Point<dim>> points = fe_p.get_unit_support_points();
 
-    const Point<dim> centroid =
-      fe_p.reference_cell().template barycenter<dim>();
+    const auto       reference_cell = fe_p.reference_cell();
+    const Point<dim> centroid       = reference_cell.template barycenter<dim>();
 
     switch (dim)
       {
@@ -85,11 +85,27 @@ namespace FE_P_BubblesImplementation
           {
             if (degree == 2)
               {
-                const double q13 = 1.0 / 3.0;
-                points.emplace_back(q13, q13, 0.0);
-                points.emplace_back(q13, 0.0, q13);
-                points.emplace_back(0.0, q13, q13);
-                points.emplace_back(q13, q13, q13);
+                for (const auto &face_no : reference_cell.face_indices())
+                  {
+                    Point<dim> midpoint;
+                    for (const auto face_vertex_no :
+                         reference_cell.face_reference_cell(0).vertex_indices())
+                      {
+                        const auto vertex_no =
+                          reference_cell.face_to_cell_vertices(
+                            face_no,
+                            face_vertex_no,
+                            ReferenceCell::default_combined_face_orientation());
+
+                        midpoint +=
+                          reference_cell.template vertex<dim>(vertex_no);
+                      }
+
+                    midpoint /=
+                      reference_cell.face_reference_cell(0).n_vertices();
+                    points.push_back(midpoint);
+                  }
+
                 points.push_back(centroid);
               }
             return points;
@@ -117,8 +133,8 @@ namespace FE_P_BubblesImplementation
   BarycentricPolynomials<dim>
   get_basis(const unsigned int degree)
   {
-    const Point<dim> centroid =
-      ReferenceCells::get_simplex<dim>().template barycenter<dim>();
+    const auto       reference_cell = ReferenceCells::get_simplex<dim>();
+    const Point<dim> centroid       = reference_cell.template barycenter<dim>();
 
     auto M = [](const unsigned int d) {
       return BarycentricPolynomial<dim, double>::monomial(d);
@@ -140,8 +156,8 @@ namespace FE_P_BubblesImplementation
 
             // in 2D and 3D we add a centroid bubble function
             auto c_bubble = BarycentricPolynomial<dim>() + 1;
-            for (unsigned int d = 0; d < dim + 1; ++d)
-              c_bubble = c_bubble * M(d);
+            for (const auto &vertex : reference_cell.vertex_indices())
+              c_bubble = c_bubble * M(vertex);
             c_bubble = c_bubble / c_bubble.value(centroid);
 
             std::vector<BarycentricPolynomial<dim>> bubble_functions;
@@ -154,14 +170,22 @@ namespace FE_P_BubblesImplementation
                 // need 'face bubble' functions in addition to the centroid.
                 // Furthermore we need to subtract them off from the other
                 // functions so that we end up with an interpolatory basis
-                auto b0 = 27 * M(0) * M(1) * M(2);
-                bubble_functions.push_back(b0 - b0.value(centroid) * c_bubble);
-                auto b1 = 27 * M(0) * M(1) * M(3);
-                bubble_functions.push_back(b1 - b1.value(centroid) * c_bubble);
-                auto b2 = 27 * M(0) * M(2) * M(3);
-                bubble_functions.push_back(b2 - b2.value(centroid) * c_bubble);
-                auto b3 = 27 * M(1) * M(2) * M(3);
-                bubble_functions.push_back(b3 - b3.value(centroid) * c_bubble);
+                for (const auto &face_no : reference_cell.face_indices())
+                  {
+                    std::vector<unsigned int> vertices;
+                    for (const auto face_vertex_no :
+                         reference_cell.face_reference_cell(0).vertex_indices())
+                      vertices.push_back(reference_cell.face_to_cell_vertices(
+                        face_no,
+                        face_vertex_no,
+                        ReferenceCell::default_combined_face_orientation()));
+
+                    Assert(vertices.size() == 3, ExcInternalError());
+                    auto b =
+                      27.0 * M(vertices[0]) * M(vertices[1]) * M(vertices[2]);
+                    bubble_functions.push_back(b -
+                                               b.value(centroid) * c_bubble);
+                  }
 
                 bubble_functions.push_back(c_bubble);
               }