Generate second-order node maps from the linear tables

include/Makefile.in

@@ -715,6 +715,7 @@ include_HEADERS = \
         fe/fe_lagrange_shape_1D.h \
         fe/fe_macro.h \
         fe/fe_map.h \
+        fe/fe_reference_element_traits.h \
         fe/fe_transformation_base.h \
         fe/fe_type.h \
         fe/fe_xyz_map.h \

include/fe/fe_reference_element_traits.h

@@ -0,0 +1,170 @@
+// The libMesh Finite Element Library.
+// Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
+
+// This library is free software; you can 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.
+
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
+
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+
+#ifndef LIBMESH_FE_REFERENCE_ELEMENT_TRAITS_H
+#define LIBMESH_FE_REFERENCE_ELEMENT_TRAITS_H
+
+#include "libmesh/libmesh.h" // invalid_uint
+
+namespace libMesh
+{
+
+/**
+ * The second-order elements' node maps are not independent facts.
+ * Their edge_nodes_map is the first-order equivalent's edge_nodes_map
+ * with the mid-edge node appended, and mid-edge nodes are numbered
+ * after the vertices in edge order.  Their side_nodes_map row is the
+ * first-order equivalent's vertex row, then the mid-edge node of each
+ * consecutive vertex pair, then the face node if the side has one.
+ *
+ * The functions here compute those tables at compile time from the
+ * first-order tables, so a second-order element class stores neither:
+ * it holds the derived tables as static constexpr members and exposes
+ * them through its usual side_nodes_map / edge_nodes_map names as
+ * references to the underlying arrays.  Only the first-order tables
+ * and the face-node rules are written by hand.
+ */
+
+/**
+ * The value the element classes pad short side_nodes_map rows with,
+ * e.g. the triangular sides of a Prism.
+ */
+static constexpr unsigned int unused_side_node = 99;
+
+/**
+ * A fixed-size table that constexpr functions can return.  Element
+ * classes bind their side_nodes_map / edge_nodes_map references to
+ * \p values.
+ */
+template <unsigned int Rows, unsigned int Cols>
+struct NodeMapTable
+{
+  unsigned int values[Rows][Cols];
+};
+
+/**
+ * The face-node rule for elements with no face nodes.
+ */
+constexpr unsigned int no_face_node (const unsigned int)
+{
+  return invalid_uint;
+}
+
+/**
+ * \returns The edge_nodes_map of the second-order element whose
+ * first-order equivalent is \p FirstOrder: each edge's two vertices, then
+ * its mid-edge node, numbered after the vertices in edge order.
+ */
+template <class FirstOrder, unsigned int Edges>
+constexpr NodeMapTable<Edges, 3>
+derived_edge_nodes ()
+{
+  NodeMapTable<Edges, 3> t {};
+  for (unsigned int e = 0; e != Edges; ++e)
+    {
+      t.values[e][0] = FirstOrder::edge_nodes_map[e][0];
+      t.values[e][1] = FirstOrder::edge_nodes_map[e][1];
+      t.values[e][2] = FirstOrder::num_nodes + e;
+    }
+  return t;
+}
+
+/**
+ * \returns The number of vertices on side \p s of a first-order
+ * element, i.e. the entries of its side_nodes_map row that aren't
+ * padding.
+ */
+template <class FirstOrder>
+constexpr unsigned int n_side_vertices (const unsigned int s)
+{
+  unsigned int n = 0;
+  for (unsigned int k = 0; k != FirstOrder::nodes_per_side; ++k)
+    if (FirstOrder::side_nodes_map[s][k] != unused_side_node)
+      ++n;
+  return n;
+}
+
+/**
+ * \returns The third entry of the row of \p edges joining vertices \p a
+ * and \p b, i.e. that edge's mid-edge node.
+ */
+template <unsigned int Edges>
+constexpr unsigned int mid_edge_node (const unsigned int (&edges)[Edges][3],
+                                      const unsigned int a,
+                                      const unsigned int b)
+{
+  for (unsigned int e = 0; e != Edges; ++e)
+    if ((edges[e][0] == a && edges[e][1] == b) ||
+        (edges[e][0] == b && edges[e][1] == a))
+      return edges[e][2];
+  return invalid_uint;
+}
+
+/**
+ * \returns The side_nodes_map of a 3D second-order element with
+ * \p Sides sides of up to \p Cols nodes each, whose first-order
+ * equivalent is \p FirstOrder, whose edge_nodes_map is \p edges, and whose
+ * \p face_node(s) is the node at the center of side \p s (or
+ * \p invalid_uint if there is none).
+ */
+template <class FirstOrder, unsigned int Sides, unsigned int Cols,
+          unsigned int Edges, class FaceNode>
+constexpr NodeMapTable<Sides, Cols>
+derived_side_nodes (const unsigned int (&edges)[Edges][3],
+                    FaceNode face_node)
+{
+  NodeMapTable<Sides, Cols> t {};
+  for (unsigned int s = 0; s != Sides; ++s)
+    {
+      const unsigned int nv = n_side_vertices<FirstOrder>(s);
+      unsigned int n = 0;
+      for (unsigned int k = 0; k != nv; ++k)
+        t.values[s][n++] = FirstOrder::side_nodes_map[s][k];
+      for (unsigned int k = 0; k != nv; ++k)
+        t.values[s][n++] = mid_edge_node(edges,
+                                         FirstOrder::side_nodes_map[s][k],
+                                         FirstOrder::side_nodes_map[s][(k+1) % nv]);
+      if (face_node(s) != invalid_uint)
+        t.values[s][n++] = face_node(s);
+      for (; n != Cols; ++n)
+        t.values[s][n] = unused_side_node;
+    }
+  return t;
+}
+
+/**
+ * \returns The side_nodes_map of a 2D second-order element whose
+ * first-order equivalent is \p FirstOrder: each side's two vertices, then
+ * its mid-side node, numbered after the vertices in side order.
+ */
+template <class FirstOrder, unsigned int Sides>
+constexpr NodeMapTable<Sides, 3>
+derived_side_nodes ()
+{
+  NodeMapTable<Sides, 3> t {};
+  for (unsigned int s = 0; s != Sides; ++s)
+    {
+      t.values[s][0] = FirstOrder::side_nodes_map[s][0];
+      t.values[s][1] = FirstOrder::side_nodes_map[s][1];
+      t.values[s][2] = FirstOrder::num_nodes + s;
+    }
+  return t;
+}
+
+} // namespace libMesh
+
+#endif // LIBMESH_FE_REFERENCE_ELEMENT_TRAITS_H

include/geom/cell_hex20.h

@@ -22,6 +22,8 @@
 
 // Local includes
 #include "libmesh/cell_hex.h"
+#include "libmesh/cell_hex8.h"
+#include "libmesh/fe_reference_element_traits.h"
 
 namespace libMesh
 {
@@ -217,16 +219,18 @@ class Hex20 final : public Hex
   static const int nodes_per_edge = 3;
 
   /**
-   * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ side to
-   * element node numbers.
+   * These map the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ edge or
+   * side to element node numbers.  They are derived from the
+   * first-order Hex8 tables; see fe_reference_element_traits.h.
    */
-  static const unsigned int side_nodes_map[num_sides][nodes_per_side];
-
-  /**
-   * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ edge to
-   * element node numbers.
-   */
-  static const unsigned int edge_nodes_map[num_edges][nodes_per_edge];
+  static constexpr NodeMapTable<num_edges, nodes_per_edge>
+  _edge_nodes = derived_edge_nodes<Hex8, num_edges>();
+  static constexpr const unsigned int (&edge_nodes_map)[num_edges][nodes_per_edge] = _edge_nodes.values;
+
+  static constexpr NodeMapTable<num_sides, nodes_per_side>
+  _side_nodes = derived_side_nodes<Hex8, num_sides, nodes_per_side>
+    (_edge_nodes.values, no_face_node);
+  static constexpr const unsigned int (&side_nodes_map)[num_sides][nodes_per_side] = _side_nodes.values;
 
   /**
    * A specialization for computing the volume of a Hex20.

include/geom/cell_hex27.h

@@ -22,6 +22,8 @@
 
 // Local includes
 #include "libmesh/cell_hex.h"
+#include "libmesh/cell_hex8.h"
+#include "libmesh/fe_reference_element_traits.h"
 
 namespace libMesh
 {
@@ -232,16 +234,18 @@ class Hex27 final : public Hex
   static const int nodes_per_edge = 3;
 
   /**
-   * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ side to
-   * element node numbers.
+   * These map the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ edge or
+   * side to element node numbers.  They are derived from the
+   * first-order Hex8 tables; see fe_reference_element_traits.h.
    */
-  static const unsigned int side_nodes_map[num_sides][nodes_per_side];
-
-  /**
-   * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ edge to
-   * element node numbers.
-   */
-  static const unsigned int edge_nodes_map[num_edges][nodes_per_edge];
+  static constexpr NodeMapTable<num_edges, nodes_per_edge>
+  _edge_nodes = derived_edge_nodes<Hex8, num_edges>();
+  static constexpr const unsigned int (&edge_nodes_map)[num_edges][nodes_per_edge] = _edge_nodes.values;
+
+  static constexpr NodeMapTable<num_sides, nodes_per_side>
+  _side_nodes = derived_side_nodes<Hex8, num_sides, nodes_per_side>
+    (_edge_nodes.values, [](unsigned int s) { return 20 + s; });
+  static constexpr const unsigned int (&side_nodes_map)[num_sides][nodes_per_side] = _side_nodes.values;
 
   /**
    * A specialization for computing the volume of a Hex27.

include/geom/cell_hex8.h

@@ -170,13 +170,35 @@ class Hex8 final : public Hex
    * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ side to
    * element node numbers.
    */
-  static const unsigned int side_nodes_map[num_sides][nodes_per_side];
+  static constexpr unsigned int side_nodes_map[num_sides][nodes_per_side] =
+    {
+      {0, 3, 2, 1}, // Side 0
+      {0, 1, 5, 4}, // Side 1
+      {1, 2, 6, 5}, // Side 2
+      {2, 3, 7, 6}, // Side 3
+      {3, 0, 4, 7}, // Side 4
+      {4, 5, 6, 7}  // Side 5
+    };
 
   /**
    * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ edge to
    * element node numbers.
    */
-  static const unsigned int edge_nodes_map[num_edges][nodes_per_edge];
+  static constexpr unsigned int edge_nodes_map[num_edges][nodes_per_edge] =
+    {
+      {0, 1}, // Edge 0
+      {1, 2}, // Edge 1
+      {2, 3}, // Edge 2
+      {0, 3}, // Edge 3
+      {0, 4}, // Edge 4
+      {1, 5}, // Edge 5
+      {2, 6}, // Edge 6
+      {3, 7}, // Edge 7
+      {4, 5}, // Edge 8
+      {5, 6}, // Edge 9
+      {6, 7}, // Edge 10
+      {4, 7}  // Edge 11
+    };
 
   /**
    * Class static helper function that computes the centroid of a

include/geom/cell_prism15.h

@@ -22,6 +22,8 @@
 
 // Local includes
 #include "libmesh/cell_prism.h"
+#include "libmesh/cell_prism6.h"
+#include "libmesh/fe_reference_element_traits.h"
 
 namespace libMesh
 {
@@ -222,16 +224,18 @@ class Prism15 final : public Prism
   static const int nodes_per_edge = 3;
 
   /**
-   * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ side to
-   * element node numbers.
+   * These map the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ edge or
+   * side to element node numbers.  They are derived from the
+   * first-order Prism6 tables; see fe_reference_element_traits.h.
    */
-  static const unsigned int side_nodes_map[num_sides][nodes_per_side];
-
-  /**
-   * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ edge to
-   * element node numbers.
-   */
-  static const unsigned int edge_nodes_map[num_edges][nodes_per_edge];
+  static constexpr NodeMapTable<num_edges, nodes_per_edge>
+  _edge_nodes = derived_edge_nodes<Prism6, num_edges>();
+  static constexpr const unsigned int (&edge_nodes_map)[num_edges][nodes_per_edge] = _edge_nodes.values;
+
+  static constexpr NodeMapTable<num_sides, nodes_per_side>
+  _side_nodes = derived_side_nodes<Prism6, num_sides, nodes_per_side>
+    (_edge_nodes.values, no_face_node);
+  static constexpr const unsigned int (&side_nodes_map)[num_sides][nodes_per_side] = _side_nodes.values;
 
   /**
    * A specialization for computing the volume of a Prism15.

include/geom/cell_prism18.h

@@ -22,6 +22,8 @@
 
 // Local includes
 #include "libmesh/cell_prism.h"
+#include "libmesh/cell_prism6.h"
+#include "libmesh/fe_reference_element_traits.h"
 
 namespace libMesh
 {
@@ -237,16 +239,18 @@ class Prism18 final : public Prism
   static const int nodes_per_edge = 3;
 
   /**
-   * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ side to
-   * element node numbers.
+   * These map the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ edge or
+   * side to element node numbers.  They are derived from the
+   * first-order Prism6 tables; see fe_reference_element_traits.h.
    */
-  static const unsigned int side_nodes_map[num_sides][nodes_per_side];
-
-  /**
-   * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ edge to
-   * element node numbers.
-   */
-  static const unsigned int edge_nodes_map[num_edges][nodes_per_edge];
+  static constexpr NodeMapTable<num_edges, nodes_per_edge>
+  _edge_nodes = derived_edge_nodes<Prism6, num_edges>();
+  static constexpr const unsigned int (&edge_nodes_map)[num_edges][nodes_per_edge] = _edge_nodes.values;
+
+  static constexpr NodeMapTable<num_sides, nodes_per_side>
+  _side_nodes = derived_side_nodes<Prism6, num_sides, nodes_per_side>
+    (_edge_nodes.values, [](unsigned int s) { return (s >= 1 && s <= 3) ? 14 + s : invalid_uint; });
+  static constexpr const unsigned int (&side_nodes_map)[num_sides][nodes_per_side] = _side_nodes.values;
 
   /**
    * A specialization for computing the volume of a Prism18.