Generate second-order node maps from the linear tables by rochi00 · Pull Request #4554 · libMesh/libmesh

rochi00 · 2026-09-16T16:44:18Z

Put the second-order element side/edge topology and reference-node
coordinates into one constexpr derivation in
fe_reference_element_traits.h: a second-order side row is the linear
corners, then the edge midpoints from the edge table, then the face
centre; higher-order reference nodes are subentity centroids.

The element classes' side_nodes_map/edge_nodes_map are now generated
from that derivation at compile time instead of being 26 hand-typed
literals. Every table is bit-identical to the one it replaces and is
constant-initialized. The maps keep their names, their [side][node]
indexing and their public visibility (MOOSE ray tracing reads them
directly); only their type changes, to a ReferenceElementTable that
indexes like the old arrays.

Put the reference-element side/edge topology into constexpr LIBMESH_DEVICE_INLINE tables and lookups in a new fe_reference_element_traits.h header, usable from host code and Kokkos kernels alike. Only the irreducible facts are stored: edge tables once per element family, side tables only for the linear elements. A second-order side row is fully determined as the linear corners, then the midpoint of each consecutive corner pair from the edge table, then the face's center node, so the side and edge topologies cannot disagree; each type's full table is materialized from that derivation at compile time, keeping runtime lookups direct indexing.

Only the vertices' reference coordinates are tabulated: every higher-order node sits at the centroid of its subentity's vertices (mid-edge nodes at edge midpoints, face nodes at face-corner centroids, interior nodes at the vertex centroid), so try_reference_node() computes them through the same side/edge tables instead of duplicating ~500 lines of coordinates. The cubic EDGE4 trisection nodes are the one true exception. This derivation is what flagged the Pyramid reference-point inconsistencies fixed in libMesh#4539.

Generate the second-order elements' static side_nodes_map and edge_nodes_map from the shared constexpr traits instead of literal initializers, so the host element classes and Kokkos FE code agree on one copy of the reference-element topology. The maps keep their names, their [side][node] indexing, and their public visibility (MOOSE's ray tracing reads Hex20::side_nodes_map and friends directly); only their type changes, from a raw array to a ReferenceElementTable whose operator[] hands back the row array, so std::begin()/std::end() on a row keep working too. Initialization is still constant: every table is materialized at compile time.

For every element type in the ElemTest suite, check the shared tables against the Elem they describe: try_reference_node() against master_point(), and the side and edge lookups against local_side_node() / nodes_on_side() and local_edge_node() / nodes_on_edge(), counts included. Coverage is all-or-nothing per type so a partially supported element can't slip through untested. This is the check that flagged the Pyramid reference-point bugs fixed in libMesh#4539.

moosebuild · 2026-09-16T18:37:29Z

Job Coverage, step Generate coverage on 5a536de wanted to post the following:

Coverage

	54e0d5	#4554 5a536d
	Total	Total	+/-	New
Rate	65.99%	65.22%	-0.77%	100.00%
Hits	79562	78629	-933	2
Misses	40996	41931	+935	0

Diff coverage report

Full coverage report

This comment will be updated on new commits.

roystgnr

Didn't quite make it through the whole thing, but this should give you enough to chew on for a first pass. ;-)

I love this in theory, but in practice? Adding 1000 lines of code to replace under 300 is never a good sign for a refactor. I assume some of this is prerequisite for Kokkos work, giving us static lookups for data we previously had to get via virtual functions? Those cases count as genuinely new features and justify some new lines, but I'll bet we can factor this way down.

roystgnr · 2026-09-16T18:13:25Z

+{
+
+template <unsigned int N>
+struct ReferenceElementVector


This looks like an atavism that never got used?

roystgnr · 2026-09-16T18:15:20Z

   * element node numbers.
   */
-  static const unsigned int side_nodes_map[num_sides][nodes_per_side];
+  static const ReferenceElementTable<num_sides, nodes_per_side> side_nodes_map;


These are public members; is there any way we can avoid changing their type and thus breaking backwards compatibility?

roystgnr · 2026-09-16T18:19:10Z

+}
+
+LIBMESH_DEVICE_INLINE constexpr ReferenceElementTable<8, 3>
+pyramid_edge_nodes()


Let's have a less-broad name for this? pyramid_quadratic_edge_nodes() maybe? It's valid for Pyramid13/14/18, but not Pyramid5 and not for many hypothetical future subclasses like a cubic Pyramid30.

Same issue with tet, hex, and prism cases; pyramid is just where I noticed it first.

roystgnr · 2026-09-16T18:22:13Z

+
+
+LIBMESH_DEVICE_INLINE constexpr ReferenceElementTable<5, 4>
+prism6_side_nodes()


Can we move the element-specific cases to element-specific headers? E.g. cell_prism6.h for this and cell_prism.h for the quadratic version?

roystgnr · 2026-09-16T18:48:21Z

+}
+
+LIBMESH_DEVICE_INLINE bool
+requires_side_specific_topology(ElemType parent)


These should probably be static member functions (or static arrays like type_to_n_sides_map in the simpler cases) in Elem::.

roystgnr · 2026-09-16T19:06:11Z

+}
+
+LIBMESH_DEVICE_INLINE constexpr bool
+try_local_edge_node(ElemType parent,


Documentation?

roystgnr · 2026-09-16T19:09:23Z

+}
+
+LIBMESH_DEVICE_INLINE constexpr ElemType
+linear_sibling_or_invalid(ElemType type)


Like "parent", AMR/C called dibs on "sibling" first.

Is this just Elem::first_order_equivalent_type(), but rewritten with a different name? We cannot keep letting AI try to slip that by us! That does not get us from an API with N functions to one with N², it gets us from N to ∞.

roystgnr · 2026-09-16T19:15:53Z

+  }
+}
+
+// Corner k of a linear element's side, from the stored linear tables.


What does "corner" mean here? vertex? But we don't seem to be doing any tests on k here; this will happily return a non-vertex.

This also seems like something that could be a composition rather than a single function.

roystgnr · 2026-09-16T19:17:39Z

+try_local_edge_node(ElemType parent,
+                    unsigned int edge,
+                    unsigned int edge_node,
+                    unsigned int & node)


We have an invalid_uint; could we just return that for the false case rather than a separate bool?

roystgnr · 2026-09-16T19:33:09Z

+}
+
+LIBMESH_DEVICE_INLINE constexpr bool
+try_local_edge_node(ElemType parent,


We can drop the try_ on all these names; this is the same thing as our local_edge_node but taking an input ElemType (which again shouldn't be called "parent") rather than Elem * this, right? We can just overload Elem::local_edge_node with a static, right? I'd prefer to throw an error (ideally there'd be some way to libmesh_assert() in runtime invocations while doing static_assert at compile-time, though I confess I'm not sure how) when we're asking for an invalid node, but if we really do need to call functions like these in cases where it might return false, we can return invalid_uint instead.

rochi00 added 5 commits September 16, 2026 09:23

Re-bootstrap

5a536de

roystgnr requested changes Sep 16, 2026

View reviewed changes



		LIBMESH_DEVICE_INLINE constexpr ReferenceElementTable<5, 4>
		prism6_side_nodes()

Uh oh!

Conversation

rochi00 commented Sep 16, 2026

Uh oh!

moosebuild commented Sep 16, 2026

Coverage

Uh oh!

roystgnr left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants