From 5aba74f407c92e0ced0237467afe53458de2a0ee Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Tue, 1 Sep 2026 13:38:13 -0500 Subject: [PATCH 1/7] Log partition-of-unity test I don't think this is a slow one but it's worth logging. --- tests/fe/fe_test.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tests/fe/fe_test.h b/tests/fe/fe_test.h index 468d982d5f..b61f19ebf7 100644 --- a/tests/fe/fe_test.h +++ b/tests/fe/fe_test.h @@ -670,6 +670,8 @@ class FETest : public FETestBase { void testPartitionOfUnity() { + LOG_UNIT_TEST; + if (!this->_elem) return; From 251806e26c8835f7e06ed9fe90cf03d02177fc50 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Tue, 1 Sep 2026 13:46:38 -0500 Subject: [PATCH 2/7] Fix broken Pyramid master points I'm not sure how we had the wrong data here. --- src/geom/cell_pyramid.C | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/geom/cell_pyramid.C b/src/geom/cell_pyramid.C index 04bdf5298a..62d686acc5 100644 --- a/src/geom/cell_pyramid.C +++ b/src/geom/cell_pyramid.C @@ -44,10 +44,10 @@ const Real Pyramid::_master_points[14][3] = {1, 0, 0}, {0, 1, 0}, {-1, 0, 0}, - {0, -0.5, 0.5}, - {0.5, 0, 0.5}, - {0, 0.5, 0.5}, - {-0.5, 0, 0.5}, + {-0.5, -0.5, 0.5}, + {0.5, -0.5, 0.5}, + {0.5, 0.5, 0.5}, + {-0.5, 0.5, 0.5}, {0, 0, 0} }; From 9bb050888e8d491f0dbd4d983c412c9018f6bc7c Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Tue, 1 Sep 2026 15:50:25 -0500 Subject: [PATCH 3/7] Extend Pyramid::_master_points to Pyramid18 --- include/geom/cell_pyramid.h | 2 +- src/geom/cell_pyramid.C | 8 ++++++-- 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/include/geom/cell_pyramid.h b/include/geom/cell_pyramid.h index fe347d1953..7f8a9c4f39 100644 --- a/include/geom/cell_pyramid.h +++ b/include/geom/cell_pyramid.h @@ -211,7 +211,7 @@ class Pyramid : public Cell /** * Master element node locations */ - static const Real _master_points[14][3]; + static const Real _master_points[18][3]; #ifdef LIBMESH_ENABLE_AMR diff --git a/src/geom/cell_pyramid.C b/src/geom/cell_pyramid.C index 62d686acc5..1169693a8d 100644 --- a/src/geom/cell_pyramid.C +++ b/src/geom/cell_pyramid.C @@ -33,7 +33,7 @@ const int Pyramid::num_sides; const int Pyramid::num_edges; const int Pyramid::num_children; -const Real Pyramid::_master_points[14][3] = +const Real Pyramid::_master_points[18][3] = { {-1, -1, 0}, {1, -1, 0}, @@ -48,7 +48,11 @@ const Real Pyramid::_master_points[14][3] = {0.5, -0.5, 0.5}, {0.5, 0.5, 0.5}, {-0.5, 0.5, 0.5}, - {0, 0, 0} + {0, 0, 0}, + {0, -Real(2)/3, Real(1)/3}, + { Real(2)/3, 0, Real(1)/3}, + {0, Real(2)/3, Real(1)/3}, + {-Real(2)/3, 0, Real(1)/3} }; const unsigned int Pyramid::edge_sides_map[8][2] = From 07d0e879e4a699853a3b097506ce618668eee4b3 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Tue, 1 Sep 2026 16:00:58 -0500 Subject: [PATCH 4/7] Fix get_refspace_nodes(Pyramid18) --- src/fe/fe_abstract.C | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/fe/fe_abstract.C b/src/fe/fe_abstract.C index 20d435ae44..74c4ffe05e 100644 --- a/src/fe/fe_abstract.C +++ b/src/fe/fe_abstract.C @@ -584,10 +584,10 @@ void FEAbstract::get_refspace_nodes(const ElemType itemType, std::vector case PYRAMID18: { // triangle centers - nodes[14] = Point (-2/Real(3),0.,1/Real(3)); - nodes[15] = Point (0.,2/Real(3),1/Real(3)); - nodes[16] = Point (2/Real(3),0.,1/Real(3)); - nodes[17] = Point (0.,-2/Real(3),1/Real(3)); + nodes[14] = Point (0., -2/Real(3),1/Real(3)); + nodes[15] = Point (2/Real(3),0.,1/Real(3)); + nodes[16] = Point (0.,2/Real(3),1/Real(3)); + nodes[17] = Point (-2/Real(3),0.,1/Real(3)); libmesh_fallthrough(); } From 70021d5a1c7d5a5fe607d5a2d4fd3ef1ddb68382 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Tue, 1 Sep 2026 21:48:53 -0500 Subject: [PATCH 5/7] Fix master_point() on Tri7/Prism20/Prism21 too --- include/geom/cell_prism.h | 2 +- include/geom/face_tri.h | 2 +- src/geom/cell_prism.C | 7 +++++-- src/geom/face_tri.C | 5 +++-- 4 files changed, 10 insertions(+), 6 deletions(-) diff --git a/include/geom/cell_prism.h b/include/geom/cell_prism.h index fbe7688280..20bfa80560 100644 --- a/include/geom/cell_prism.h +++ b/include/geom/cell_prism.h @@ -211,7 +211,7 @@ class Prism : public Cell /** * Master element node locations */ - static const Real _master_points[18][3]; + static const Real _master_points[21][3]; /** * This maps the \f$ j^{th} \f$ Prism vertex to the (in this case) 3 diff --git a/include/geom/face_tri.h b/include/geom/face_tri.h index bfc249a0cb..595ee04014 100644 --- a/include/geom/face_tri.h +++ b/include/geom/face_tri.h @@ -205,7 +205,7 @@ class Tri : public Face /** * Master element node locations */ - static const Real _master_points[6][3]; + static const Real _master_points[7][3]; /** * This maps the \f$ j^{th} \f$ node to the (in this case) 2 side diff --git a/src/geom/cell_prism.C b/src/geom/cell_prism.C index 618d681051..180720f34c 100644 --- a/src/geom/cell_prism.C +++ b/src/geom/cell_prism.C @@ -34,7 +34,7 @@ const int Prism::num_sides; const int Prism::num_edges; const int Prism::num_children; -const Real Prism::_master_points[18][3] = +const Real Prism::_master_points[21][3] = { {0, 0, -1}, {1, 0, -1}, @@ -53,7 +53,10 @@ const Real Prism::_master_points[18][3] = {0, 0.5, 1}, {0.5, 0, 0}, {0.5, 0.5, 0}, - {0, 0.5, 0} + {0, 0.5, 0}, + {Real(1)/3, Real(1)/3, -1}, + {Real(1)/3, Real(1)/3, 1}, + {Real(1)/3, Real(1)/3, 0} }; const unsigned int Prism::edge_sides_map[9][2] = diff --git a/src/geom/face_tri.C b/src/geom/face_tri.C index f63a640954..b6ebd6a879 100644 --- a/src/geom/face_tri.C +++ b/src/geom/face_tri.C @@ -35,14 +35,15 @@ const int Tri::num_children; // Note: we can omit initialization of the third entry of each row because // static variables are automatically zero-initialized. -const Real Tri::_master_points[6][3] = +const Real Tri::_master_points[7][3] = { {0, 0}, {1, 0}, {0, 1}, {0.5, 0}, {0.5, 0.5}, - {0, 0.5} + {0, 0.5}, + {Real(1)/3, Real(1)/3} }; const unsigned int Tri::adjacent_sides_map[/*num_vertices*/3][/*n_adjacent_sides*/2] = From 3bddc3557443b0893944f6670ea5a440c7ecd318 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Tue, 1 Sep 2026 13:37:27 -0500 Subject: [PATCH 6/7] Add reference-point unit tests Rochi found bugs in our higher-order Pyramids reference elements, and these tests catch them. --- tests/fe/fe_test.h | 16 ++++++++++++++++ tests/geom/elem_test.C | 18 ++++++++++++++++++ 2 files changed, 34 insertions(+) diff --git a/tests/fe/fe_test.h b/tests/fe/fe_test.h index b61f19ebf7..7f5db6861d 100644 --- a/tests/fe/fe_test.h +++ b/tests/fe/fe_test.h @@ -25,6 +25,7 @@ #define FETEST \ CPPUNIT_TEST( testFEInterface ); \ + CPPUNIT_TEST( testRefspaceNodes ); \ CPPUNIT_TEST( testU ); \ CPPUNIT_TEST( testPartitionOfUnity ); \ CPPUNIT_TEST( testGradU ); \ @@ -668,6 +669,21 @@ class FETest : public FETestBase { #endif // LIBMESH_ENABLE_EXCEPTIONS } + void testRefspaceNodes() + { + LOG_UNIT_TEST; + + if (!this->_elem) + return; + + std::vector nodes; + this->_fe->get_refspace_nodes(this->_elem->type(), nodes); + for (auto n : index_range(nodes)) + LIBMESH_ASSERT_REALVEC_EQUAL(nodes[n], + this->_elem->master_point(n), + TOLERANCE*TOLERANCE); + } + void testPartitionOfUnity() { LOG_UNIT_TEST; diff --git a/tests/geom/elem_test.C b/tests/geom/elem_test.C index 41ced6f6b7..7a52880f50 100644 --- a/tests/geom/elem_test.C +++ b/tests/geom/elem_test.C @@ -50,6 +50,23 @@ public: } } + void test_ref_elem() + { + LOG_UNIT_TEST; + + for (const auto & elem : this->_mesh->active_local_element_ptr_range()) + { + const Elem & ref_elem = *elem->reference_elem(); + + for (auto i : elem->node_index_range()) + { + LIBMESH_ASSERT_REALVEC_EQUAL(elem->master_point(i), + ref_elem.point(i), + TOLERANCE*TOLERANCE); + } + } + } + void test_quality() { LOG_UNIT_TEST; @@ -944,6 +961,7 @@ public: #define ELEMTEST \ CPPUNIT_TEST( test_bounding_box ); \ + CPPUNIT_TEST( test_ref_elem ); \ CPPUNIT_TEST( test_quality ); \ CPPUNIT_TEST( test_node_edge_map_consistency ); \ CPPUNIT_TEST( test_maps ); \ From c3bcaaa0f53f30d58a291773f32db52829178f9f Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Tue, 1 Sep 2026 16:17:05 -0500 Subject: [PATCH 7/7] Consolidate redundant definition of reference pts --- src/fe/fe_abstract.C | 228 ++--------------------------------------- tests/fe/fe_test.h | 8 ++ tests/geom/elem_test.C | 8 ++ 3 files changed, 22 insertions(+), 222 deletions(-) diff --git a/src/fe/fe_abstract.C b/src/fe/fe_abstract.C index 74c4ffe05e..741e0cc089 100644 --- a/src/fe/fe_abstract.C +++ b/src/fe/fe_abstract.C @@ -31,6 +31,7 @@ #include "libmesh/periodic_boundary.h" #include "libmesh/quadrature.h" #include "libmesh/quadrature_gauss.h" +#include "libmesh/reference_elem.h" #include "libmesh/remote_elem.h" #include "libmesh/tensor_value.h" #include "libmesh/threads.h" @@ -405,230 +406,13 @@ void FEAbstract::get_refspace_nodes(const ElemType itemType, std::vector Utility::enum_to_string(itemType)); nodes.resize(n_nodes); - switch(itemType) - { - case NODEELEM: - { - nodes[0] = Point (0.,0.,0.); - return; - } - case EDGE3: - { - nodes[2] = Point (0.,0.,0.); - libmesh_fallthrough(); - } - case EDGE2: - { - nodes[0] = Point (-1.,0.,0.); - nodes[1] = Point (1.,0.,0.); - return; - } - case EDGE4: // not nested with EDGE3 - { - nodes[0] = Point (-1.,0.,0.); - nodes[1] = Point (1.,0.,0.); - nodes[2] = Point (-1./3.,0.,0.); - nodes[3] - Point (1./3.,0.,0.); - return; - } - case TRI7: - { - nodes[6] = Point (1./3.,1./3.,0.); - libmesh_fallthrough(); - } - case TRI6: - { - nodes[3] = Point (.5,0.,0.); - nodes[4] = Point (.5,.5,0.); - nodes[5] = Point (0.,.5,0.); - libmesh_fallthrough(); - } - case TRI3: - case TRISHELL3: - { - nodes[0] = Point (0.,0.,0.); - nodes[1] = Point (1.,0.,0.); - nodes[2] = Point (0.,1.,0.); - return; - } - case QUAD9: - case QUADSHELL9: - { - nodes[8] = Point (0.,0.,0.); - libmesh_fallthrough(); - } - case QUAD8: - case QUADSHELL8: - { - nodes[4] = Point (0.,-1.,0.); - nodes[5] = Point (1.,0.,0.); - nodes[6] = Point (0.,1.,0.); - nodes[7] = Point (-1.,0.,0.); - libmesh_fallthrough(); - } - case QUAD4: - case QUADSHELL4: - { - nodes[0] = Point (-1.,-1.,0.); - nodes[1] = Point (1.,-1.,0.); - nodes[2] = Point (1.,1.,0.); - nodes[3] = Point (-1.,1.,0.); - return; - } - case TET14: - { - nodes[10] = Point (1/Real(3),1/Real(3),0.); - nodes[11] = Point (1/Real(3),0.,1/Real(3)); - nodes[12] = Point (1/Real(3),1/Real(3),1/Real(3)); - nodes[13] = Point (0.,1/Real(3),1/Real(3)); - libmesh_fallthrough(); - } - case TET10: - { - nodes[4] = Point (.5,0.,0.); - nodes[5] = Point (.5,.5,0.); - nodes[6] = Point (0.,.5,0.); - nodes[7] = Point (0.,0.,.5); - nodes[8] = Point (.5,0.,.5); - nodes[9] = Point (0.,.5,.5); - libmesh_fallthrough(); - } - case TET4: - { - nodes[0] = Point (0.,0.,0.); - nodes[1] = Point (1.,0.,0.); - nodes[2] = Point (0.,1.,0.); - nodes[3] = Point (0.,0.,1.); - return; - } - case HEX27: - { - nodes[20] = Point (0.,0.,-1.); - nodes[21] = Point (0.,-1.,0.); - nodes[22] = Point (1.,0.,0.); - nodes[23] = Point (0.,1.,0.); - nodes[24] = Point (-1.,0.,0.); - nodes[25] = Point (0.,0.,1.); - nodes[26] = Point (0.,0.,0.); - libmesh_fallthrough(); - } - case HEX20: - { - nodes[8] = Point (0.,-1.,-1.); - nodes[9] = Point (1.,0.,-1.); - nodes[10] = Point (0.,1.,-1.); - nodes[11] = Point (-1.,0.,-1.); - nodes[12] = Point (-1.,-1.,0.); - nodes[13] = Point (1.,-1.,0.); - nodes[14] = Point (1.,1.,0.); - nodes[15] = Point (-1.,1.,0.); - nodes[16] = Point (0.,-1.,1.); - nodes[17] = Point (1.,0.,1.); - nodes[18] = Point (0.,1.,1.); - nodes[19] = Point (-1.,0.,1.); - libmesh_fallthrough(); - } - case HEX8: - { - nodes[0] = Point (-1.,-1.,-1.); - nodes[1] = Point (1.,-1.,-1.); - nodes[2] = Point (1.,1.,-1.); - nodes[3] = Point (-1.,1.,-1.); - nodes[4] = Point (-1.,-1.,1.); - nodes[5] = Point (1.,-1.,1.); - nodes[6] = Point (1.,1.,1.); - nodes[7] = Point (-1.,1.,1.); - return; - } - case PRISM21: - { - nodes[20] = Point (1/Real(3),1/Real(3),0); - libmesh_fallthrough(); - } - case PRISM20: - { - nodes[18] = Point (1/Real(3),1/Real(3),-1); - nodes[19] = Point (1/Real(3),1/Real(3),1); - libmesh_fallthrough(); - } - case PRISM18: - { - nodes[15] = Point (.5,0.,0.); - nodes[16] = Point (.5,.5,0.); - nodes[17] = Point (0.,.5,0.); - libmesh_fallthrough(); - } - case PRISM15: - { - nodes[6] = Point (.5,0.,-1.); - nodes[7] = Point (.5,.5,-1.); - nodes[8] = Point (0.,.5,-1.); - nodes[9] = Point (0.,0.,0.); - nodes[10] = Point (1.,0.,0.); - nodes[11] = Point (0.,1.,0.); - nodes[12] = Point (.5,0.,1.); - nodes[13] = Point (.5,.5,1.); - nodes[14] = Point (0.,.5,1.); - libmesh_fallthrough(); - } - case PRISM6: - { - nodes[0] = Point (0.,0.,-1.); - nodes[1] = Point (1.,0.,-1.); - nodes[2] = Point (0.,1.,-1.); - nodes[3] = Point (0.,0.,1.); - nodes[4] = Point (1.,0.,1.); - nodes[5] = Point (0.,1.,1.); - return; - } - case PYRAMID18: - { - // triangle centers - nodes[14] = Point (0., -2/Real(3),1/Real(3)); - nodes[15] = Point (2/Real(3),0.,1/Real(3)); - nodes[16] = Point (0.,2/Real(3),1/Real(3)); - nodes[17] = Point (-2/Real(3),0.,1/Real(3)); - - libmesh_fallthrough(); - } - case PYRAMID14: - { - // base center - nodes[13] = Point (0.,0.,0.); - libmesh_fallthrough(); - } - case PYRAMID13: - { - // base midedge - nodes[5] = Point (0.,-1.,0.); - nodes[6] = Point (1.,0.,0.); - nodes[7] = Point (0.,1.,0.); - nodes[8] = Point (-1,0.,0.); - - // lateral midedge - nodes[9] = Point (-.5,-.5,.5); - nodes[10] = Point (.5,-.5,.5); - nodes[11] = Point (.5,.5,.5); - nodes[12] = Point (-.5,.5,.5); - - libmesh_fallthrough(); - } - case PYRAMID5: - { - // base corners - nodes[0] = Point (-1.,-1.,0.); - nodes[1] = Point (1.,-1.,0.); - nodes[2] = Point (1.,1.,0.); - nodes[3] = Point (-1.,1.,0.); - // apex - nodes[4] = Point (0.,0.,1.); - return; - } + // For any element with a well-defined n_nodes we should also have a + // well-defined reference element. + const Elem & ref_elem = ReferenceElem::get(itemType); - default: - libmesh_error_msg("ERROR: Unknown element type " << Utility::enum_to_string(itemType)); - } + for (auto i : make_range(n_nodes)) + nodes[i] = ref_elem.point(i); } diff --git a/tests/fe/fe_test.h b/tests/fe/fe_test.h index 7f5db6861d..d2a6e49836 100644 --- a/tests/fe/fe_test.h +++ b/tests/fe/fe_test.h @@ -676,6 +676,14 @@ class FETest : public FETestBase { if (!this->_elem) return; + // We don't currently support reference infinite elements + if (this->_elem->infinite()) + return; + + // We'll probably never support reference polytopes + if (this->_elem->runtime_topology()) + return; + std::vector nodes; this->_fe->get_refspace_nodes(this->_elem->type(), nodes); for (auto n : index_range(nodes)) diff --git a/tests/geom/elem_test.C b/tests/geom/elem_test.C index 7a52880f50..943e4f8554 100644 --- a/tests/geom/elem_test.C +++ b/tests/geom/elem_test.C @@ -56,6 +56,14 @@ public: for (const auto & elem : this->_mesh->active_local_element_ptr_range()) { + // We don't currently support reference infinite elements + if (elem->infinite()) + continue; + + // We'll probably never support reference polytopes + if (elem->runtime_topology()) + continue; + const Elem & ref_elem = *elem->reference_elem(); for (auto i : elem->node_index_range())