From 8acc9549dd1389a421654dc1653bf24b76fb2b9b Mon Sep 17 00:00:00 2001
From: David Wells <drwells@email.unc.edu>
Date: Fri, 16 Dec 2022 19:19:01 -0500
Subject: [PATCH] ReferenceCell: port line_to_cell_vertices() over

---
 include/deal.II/grid/reference_cell.h | 88 +++++++++++++++++++++++++++
 source/fe/mapping_q.cc                |  6 +-
 2 files changed, 91 insertions(+), 3 deletions(-)

diff --git a/include/deal.II/grid/reference_cell.h b/include/deal.II/grid/reference_cell.h
index b149d0dc64..c44a65188c 100644
--- a/include/deal.II/grid/reference_cell.h
+++ b/include/deal.II/grid/reference_cell.h
@@ -423,6 +423,18 @@ public:
   std::array<unsigned int, 2>
   standard_line_to_face_and_line_index(const unsigned int line) const;
 
+  /**
+   * Map line vertex number to cell vertex number, i.e., return the cell vertex
+   * number of the <tt>vertex</tt>th vertex of line <tt>line</tt>.
+   *
+   * The order of the lines, as well as their direction (which in turn
+   * determines which vertices are first and second on a line) is the canonical
+   * one in deal.II, as described in the general documentation of this class.
+   */
+  unsigned int
+  line_to_cell_vertices(const unsigned int line,
+                        const unsigned int vertex) const;
+
   /**
    * Map face line number to cell line number.
    */
@@ -1491,6 +1503,82 @@ ReferenceCell::standard_line_to_face_and_line_index(
 }
 
 
+inline unsigned int
+ReferenceCell::line_to_cell_vertices(const unsigned int line,
+                                     const unsigned int vertex) const
+{
+  AssertIndexRange(vertex, 2);
+  AssertIndexRange(line, n_lines());
+
+  if (*this == ReferenceCells::Vertex)
+    return vertex;
+  else if (*this == ReferenceCells::Line)
+    return vertex;
+  else if (*this == ReferenceCells::Quadrilateral)
+    {
+      static constexpr ndarray<unsigned int, 4, 2> table = {
+        {{{0, 2}}, {{1, 3}}, {{0, 1}}, {{2, 3}}}};
+      return table[line][vertex];
+    }
+  else if (*this == ReferenceCells::Hexahedron)
+    {
+      // first four lines comprise the bottom face, next four are the top,
+      // and the last four are 'bottom to top'
+      static constexpr ndarray<unsigned int, 12, 2> table = {{{{0, 2}},
+                                                              {{1, 3}},
+                                                              {{0, 1}},
+                                                              {{2, 3}},
+                                                              {{4, 6}},
+                                                              {{5, 7}},
+                                                              {{4, 5}},
+                                                              {{6, 7}},
+                                                              {{0, 4}},
+                                                              {{1, 5}},
+                                                              {{2, 6}},
+                                                              {{3, 7}}}};
+      return table[line][vertex];
+    }
+  else if (*this == ReferenceCells::Triangle)
+    {
+      static constexpr ndarray<unsigned int, 3, 2> table = {
+        {{{0, 1}}, {{1, 2}}, {{2, 0}}}};
+      return table[line][vertex];
+    }
+  else if (*this == ReferenceCells::Tetrahedron)
+    {
+      static constexpr ndarray<unsigned int, 6, 2> table = {
+        {{{0, 1}}, {{1, 2}}, {{2, 0}}, {{0, 3}}, {{1, 3}}, {{2, 3}}}};
+      return table[line][vertex];
+    }
+  else if (*this == ReferenceCells::Pyramid)
+    {
+      static constexpr ndarray<unsigned int, 8, 2> table = {{{{0, 2}},
+                                                             {{1, 3}},
+                                                             {{0, 1}},
+                                                             {{2, 3}},
+                                                             {{4, 0}},
+                                                             {{1, 4}},
+                                                             {{2, 4}},
+                                                             {{4, 3}}}};
+      return table[line][vertex];
+    }
+  else if (*this == ReferenceCells::Wedge)
+    {
+      static constexpr ndarray<unsigned int, 9, 2> table = {{{{1, 0}},
+                                                             {{2, 1}},
+                                                             {{0, 2}},
+                                                             {{3, 4}},
+                                                             {{4, 5}},
+                                                             {{5, 3}},
+                                                             {{0, 3}},
+                                                             {{1, 4}},
+                                                             {{2, 5}}}};
+      return table[line][vertex];
+    }
+  Assert(false, ExcInternalError());
+  return numbers::invalid_unsigned_int;
+}
+
 
 inline unsigned int
 ReferenceCell::face_to_cell_lines(const unsigned int  face,
diff --git a/source/fe/mapping_q.cc b/source/fe/mapping_q.cc
index 9760fb508d..bfb4089a7a 100644
--- a/source/fe/mapping_q.cc
+++ b/source/fe/mapping_q.cc
@@ -1680,10 +1680,10 @@ MappingQ<dim, spacedim>::add_line_support_points(
                cell->get_manifold() :
                line->get_manifold());
 
+          const auto reference_cell = ReferenceCells::get_hypercube<dim>();
           const std::array<Point<spacedim>, 2> vertices{
-            {cell->vertex(GeometryInfo<dim>::line_to_cell_vertices(line_no, 0)),
-             cell->vertex(
-               GeometryInfo<dim>::line_to_cell_vertices(line_no, 1))}};
+            {cell->vertex(reference_cell.line_to_cell_vertices(line_no, 0)),
+             cell->vertex(reference_cell.line_to_cell_vertices(line_no, 1))}};
 
           const std::size_t n_rows =
             support_point_weights_perimeter_to_interior[0].size(0);
-- 
2.39.5