From 470387db15bd84343e50784e96fa12672e78da12 Mon Sep 17 00:00:00 2001 From: Patrick Behne Date: Tue, 21 Oct 2025 16:13:30 -0600 Subject: [PATCH 1/8] Added reference elements for tets. Ref #4082. --- src/systems/variational_smoother_system.C | 63 +++++++++++++++++++++++ 1 file changed, 63 insertions(+) diff --git a/src/systems/variational_smoother_system.C b/src/systems/variational_smoother_system.C index 1f182b3fb3..33b1255a0a 100644 --- a/src/systems/variational_smoother_system.C +++ b/src/systems/variational_smoother_system.C @@ -1083,6 +1083,69 @@ VariationalSmootherSystem::get_target_elem(const ElemType & type) } // if Pyramid + // Elems deriving from Tet + else if (type_str.compare(0, 3, "TET") == 0) + { + + // The ideal target element is a a regular tet with equilateral + // triangles for all faces, with volume equal to the volume of the + // reference element. + + // The volume of a tet is given by v = b * h / 3, where b is the area of + // the base face and h is the height of the apex node. The area of an + // equilateral triangle with side length s is b = sqrt(3) s^2 / 4. + // For all faces to have side length s, the height of the apex node is + // h = sqrt(2/3) * s. Then the volume is v = sqrt(2) * s^3 / 12. + // Solving for s, the side length that will preserve the volume of the + // reference element is s = (12 * v / sqrt(2))^(1/3), where v is the volume + // of the non-optimal reference element (i.e., a right tet). + + const Real sqrt_2 = std::sqrt(Real(2)); + const Real sqrt_3 = std::sqrt(Real(3)); + // Side length that preserves the volume of the reference element + const auto side_length = std::pow(12. * ref_vol / sqrt_2, 1. / 3.); + // tet height with the property that all faces are equilateral triangles + const auto target_height = sqrt_2 / sqrt_3 * side_length; + + const auto & s = side_length; + const auto & h = target_height; + + // For regular tet + // x y z node_id + owned_nodes.emplace_back(Node::build(Point(0., 0., 0.), 0)); + owned_nodes.emplace_back(Node::build(Point(s, 0., 0.), 1)); + owned_nodes.emplace_back(Node::build(Point(0.5 * s, 0.5 * sqrt_3 * s, 0.), 2)); + owned_nodes.emplace_back(Node::build(Point(0.5 * s, sqrt_3 / 6. * s, h), 3)); + + if (type == TET10 || type == TET14) + { + const auto & on = owned_nodes; + // Define the edge midpoint nodes of the tet + + // Base node to base node midpoint nodes + owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[1]) / 2.), 4)); + owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[2]) / 2.), 5)); + owned_nodes.emplace_back(Node::build(Point((*on[2] + *on[0]) / 2.), 6)); + // Base node to apex node midpoint nodes + owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[3]) / 2.), 7)); + owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[3]) / 2.), 8)); + owned_nodes.emplace_back(Node::build(Point((*on[2] + *on[3]) / 2.), 9)); + + if (type == TET14) + { + // Define the face midpoint nodes of the tet + owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[1] + *on[2]) / 3.), 10)); + owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[1] + *on[3]) / 3.), 11)); + owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[2] + *on[3]) / 3.), 12)); + owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[2] + *on[3]) / 3.), 13)); + } + } + + else if (type != TET4) + libmesh_error_msg("Unsupported tet element: " << type_str); + + } // if Tet + // Set the target_elem equal to the reference elem else for (const auto & node : target_elem->reference_elem()->node_ref_range()) From e2608fbf66ed9747cee5a40f281020c6f9352f33 Mon Sep 17 00:00:00 2001 From: Patrick Behne Date: Tue, 21 Oct 2025 16:15:45 -0600 Subject: [PATCH 2/8] Added unit tests for smoothing tets. --- .../systems/variational_smoother_constraint.h | 26 +- tests/mesh/mesh_smoother_test.C | 296 +++++++++++++++++- 2 files changed, 308 insertions(+), 14 deletions(-) diff --git a/include/systems/variational_smoother_constraint.h b/include/systems/variational_smoother_constraint.h index 27f5f5d503..3222ead25b 100644 --- a/include/systems/variational_smoother_constraint.h +++ b/include/systems/variational_smoother_constraint.h @@ -421,19 +421,6 @@ class VariationalSmootherConstraint : public System::Constraint */ void constrain_node_to_line(const Node & node, const Point & line_vec); - /** - * Given a mesh and a node in the mesh, the vector will be filled with - * every node directly attached to the given one. IF NO neighbors are found, - * all nodes on the containing side are treated as neighbors. This is useful - * when the node does not lie on an edge, such as the central face node in - * HEX27 elements. - */ - static void find_nodal_or_face_neighbors( - const MeshBase & mesh, - const Node & node, - const std::unordered_map> & nodes_to_elem_map, - std::vector & neighbors); - /** * Determines whether two neighboring nodes share a common boundary id. * @param boundary_node The first of the two prospective nodes. @@ -517,6 +504,19 @@ class VariationalSmootherConstraint : public System::Constraint virtual ~VariationalSmootherConstraint() override; virtual void constrain() override; + + /** + * Given a mesh and a node in the mesh, the vector will be filled with + * every node directly attached to the given one. IF NO neighbors are found, + * all nodes on the containing side are treated as neighbors. This is useful + * when the node does not lie on an edge, such as the central face node in + * HEX27 elements. + */ + static void find_nodal_or_face_neighbors( + const MeshBase & mesh, + const Node & node, + const std::unordered_map> & nodes_to_elem_map, + std::vector & neighbors); }; diff --git a/tests/mesh/mesh_smoother_test.C b/tests/mesh/mesh_smoother_test.C index eac1b3c6c6..41334acdfe 100644 --- a/tests/mesh/mesh_smoother_test.C +++ b/tests/mesh/mesh_smoother_test.C @@ -7,6 +7,7 @@ #include #include #include +#include #include #include #include @@ -264,6 +265,11 @@ public: CPPUNIT_TEST(testVariationalPrism21); CPPUNIT_TEST(testVariationalPrism21MultipleSubdomains); + CPPUNIT_TEST(testVariationalTet4); + CPPUNIT_TEST(testVariationalTet10); + CPPUNIT_TEST(testVariationalTet14); + CPPUNIT_TEST(testVariationalTet14MultipleSubdomains); + #if defined(LIBMESH_HAVE_GZSTREAM) CPPUNIT_TEST(testVariationalMixed2D); CPPUNIT_TEST(testVariationalMixed3D); @@ -345,6 +351,54 @@ public: } } + // Helper function to determine how many dimensions of a given point lie at + // the center and faces of a sub-cube + std::pair + numCenteredAndFacedDimensions(const Point & point, const Real & side_length, const Real & tol) + { + const auto half_length = 0.5 * side_length; + unsigned int num_centered = 0; + unsigned int num_face = 0; + for (const auto d : make_range(3)) + { + const auto num_half_lengths = point(d) / half_length; + if (std::abs(num_half_lengths - std::round(num_half_lengths)) > tol) + // This coordinante does not occur at a sub-cube face or center + continue; + + if ((unsigned int)std::round(num_half_lengths) % 2) + // An odd number of half lengths means point(d) is at a sub-cube center + ++num_centered; + else + // An odd number of half lengths means point(d) is on a sub-cube face + ++num_face; + } + + return std::make_pair(num_centered, num_face); + } + + // Helper function to determine whether a given point lies at the center of + // a sub-cube + bool pointIsCubeCenter(const Point & point, const Real & side_length, const Real & tol = 1e-3) + { + return numCenteredAndFacedDimensions(point, side_length, tol).first == 3; + } + + // Helper function to determine whether a given point lies at the vertex of + // a sub-cube + bool pointIsCubeVertex(const Point & point, const Real & side_length, const Real & tol = 1e-3) + { + return numCenteredAndFacedDimensions(point, side_length, tol).second == 3; + } + + // Helper function to determine whether a given point lies at the center of + // a sub-cube face + bool pointIsCubeFaceCenter(const Point & point, const Real & side_length, const Real & tol = 1e-3) + { + const auto result = numCenteredAndFacedDimensions(point, side_length, tol); + return (result.first == 2) && (result.second == 1); + } + void testVariationalSmoother(ReplicatedMesh & mesh, VariationalMeshSmoother & smoother, const ElemType type, @@ -364,10 +418,12 @@ public: const bool type_is_tri = Utility::enum_to_string(type).compare(0, 3, "TRI") == 0; const bool type_is_pyramid = Utility::enum_to_string(type).compare(0, 7, "PYRAMID") == 0; const bool type_is_prism = Utility::enum_to_string(type).compare(0, 5, "PRISM") == 0; + const bool type_is_tet = Utility::enum_to_string(type).compare(0, 3, "TET") == 0; // Used fewer elems for higher order types, as extra midpoint nodes will add // enough complexity unsigned int n_elems_per_side = 4 / Elem::type_to_default_order_map[type]; + const auto side_length = 1.0 / n_elems_per_side; switch (dim) { @@ -501,7 +557,8 @@ public: // For pyramids, the factor 2 accounts for the account that cubes of side // length s are divided into pyramids of base side length s and height s/2. // The factor 4 lets us catch triangular face midpoint nodes for PYRAMID18 elements. - const auto scale_factor = *elem_orders.begin() * (type_is_pyramid ? 2 * 4 : 1); + // Similar reasoning for tets. + const auto scale_factor = *elem_orders.begin() * ((type_is_pyramid || type_is_tet) ? 2 * 4 : 1); // Function to assert the node distortion is as expected auto node_distortion_is = [&n_elems_per_side, &dim, &boundary_info, &scale_factor, &type_is_prism]( @@ -628,6 +685,9 @@ public: } else { + std::unordered_map> nodes_to_elem_map; + MeshTools::build_nodes_to_elem_map(mesh, nodes_to_elem_map); + // Make sure we're not too distorted anymore std::set nodes_checked; for (const auto * elem : mesh.active_element_ptr_range()) @@ -686,6 +746,207 @@ public: } } + // Check for special case where some tet midpoint nodes are not + // smoothed to the actual midpoints. + else if (type_is_tet && !elem->is_vertex(local_node_id)) + { + // We have a non-vertex node. Determine what "type" of + // midpoint node with respect to the mesh geometry. + // First, get the nodes that neighbor this node + std::vector neighbors; + VariationalSmootherConstraint::find_nodal_or_face_neighbors( + mesh, node, nodes_to_elem_map, neighbors); + + switch (neighbors.size()) + { + case 2: { + // Midpoint node + // Determine what "type" of midpoint node + + // Is one of the neighbors at the center of a sub-cube? + const auto is_0_cube_center = + pointIsCubeCenter(*neighbors[0], side_length); + const auto is_1_cube_center = + pointIsCubeCenter(*neighbors[1], side_length); + libmesh_assert(!(is_0_cube_center && is_1_cube_center)); + + if (is_0_cube_center || is_1_cube_center) + { + const auto & cube_center = + is_0_cube_center ? *neighbors[0] : *neighbors[1]; + const auto & other = + is_0_cube_center ? *neighbors[1] : *neighbors[0]; + + // Since one neighbor is a sub-cube center, the + // other will either be a sub-cube face center or + // a sub-cube vertex. Determine which one. + if (pointIsCubeFaceCenter(other, side_length)) + { + const Real x = (type == TET10) ? 0.42895 : 0.41486; + CPPUNIT_ASSERT(node.relative_fuzzy_equals( + other + x * (cube_center - other), 1e-3)); + } + + else if (pointIsCubeVertex(other, side_length)) + { + const Real x = (type == TET10) ? 0.55389 : 0.58094; + CPPUNIT_ASSERT(node.relative_fuzzy_equals( + other + x * (cube_center - other), 1e-3)); + } + } + + else + { + // The remaining possibilities are + // cube-vertex-to-cube-vertex midpoints and + // cube-vertex-to-cube-face-center midpoints. + const auto is_0_cube_vertex = + pointIsCubeVertex(*neighbors[0], side_length); + const auto is_1_cube_vertex = + pointIsCubeVertex(*neighbors[1], side_length); + const auto is_0_cube_face_center = + pointIsCubeFaceCenter(*neighbors[0], side_length); + const auto is_1_cube_face_center = + pointIsCubeFaceCenter(*neighbors[1], side_length); + + if (is_0_cube_vertex && is_1_cube_vertex) + // Due to symmetry, this type of midpoint is + // the geometric midpoint. Let the + // node_distortion_is function check it + CPPUNIT_ASSERT(node_distortion_is(node, false, 1e-2)); + + else + { + libmesh_error_msg_if( + (is_0_cube_center || is_0_cube_face_center) && + (is_1_cube_center || is_1_cube_face_center), + "We should never get here!"); + const auto & cube_vertex = + is_0_cube_vertex ? *neighbors[0] : *neighbors[1]; + const auto & cube_face_center = + is_0_cube_face_center ? *neighbors[0] : *neighbors[1]; + const Real x = (type == TET10) ? 0.61299 : 0.65126; + CPPUNIT_ASSERT(node.relative_fuzzy_equals( + cube_vertex + x * (cube_face_center - cube_vertex), 1e-3)); + } + } + + continue; + break; + } + + case 6: { + // Face node + // Only keep the face vertex nodes + neighbors.erase(std::remove_if(neighbors.begin(), + neighbors.end(), + [&elem](const Node * n) { + return !elem->is_vertex( + elem->local_node(n->id())); + }), + neighbors.end()); + + // There are three types of faces: + // + // 1) Isoscelese triangle with vertices at two + // sub-cube vertices and the sub-cube center, + // 2) Isosceles triangle with vertices at two + // sub-cube vertices and a sub-cube face center, + // 3) Scalene triangle with vertices at a sub-cube + // center, a sub-cube vertex, and a sub-cube face + // center. + // + // Determine which kind of face this node belongs to based on the vertex + // neighbors + + unsigned int num_vertices_at_cube_center = 0; + unsigned int num_vertices_at_cube_vertices = 0; + unsigned int num_vertices_at_cube_face_centers = 0; + for (const auto * neighbor : neighbors) + { + if (pointIsCubeCenter(*neighbor, side_length)) + num_vertices_at_cube_center += 1; + else if (pointIsCubeVertex(*neighbor, side_length)) + num_vertices_at_cube_vertices += 1; + else if (pointIsCubeFaceCenter(*neighbor, side_length)) + num_vertices_at_cube_face_centers += 1; + } + + // We will express the face node at a linear + // combination of the face vertices + Node node_approx = Node(0, 0, 0); + + // Isosceles triangular face + if (num_vertices_at_cube_vertices == 2) + { + // Isosceles triangular face (type 1) + if (num_vertices_at_cube_center) + for (const auto * neighbor : neighbors) + { + Real weight; + if (pointIsCubeVertex(*neighbor, side_length)) + weight = 0.30601; + else if (pointIsCubeCenter(*neighbor, side_length)) + weight = 0.38799; + else + libmesh_error_msg("We should never get here!"); + + node_approx += weight * (*neighbor); + } + + // Isosceles triangular face (type 2) + else if (num_vertices_at_cube_face_centers) + for (const auto * neighbor : neighbors) + { + Real weight; + if (pointIsCubeVertex(*neighbor, side_length)) + weight = 0.28078; + else if (pointIsCubeFaceCenter(*neighbor, side_length)) + weight = 0.43843; + else + libmesh_error_msg("We should never get here!"); + + node_approx += weight * (*neighbor); + } + + else + libmesh_error_msg("We should never get here!"); + } + + // Scalene triangular face (type 3) + else if (num_vertices_at_cube_center && num_vertices_at_cube_vertices && + num_vertices_at_cube_face_centers) + for (const auto * neighbor : neighbors) + { + Real weight; + if (pointIsCubeCenter(*neighbor, side_length)) + weight = 0.33102; + else if (pointIsCubeVertex(*neighbor, side_length)) + weight = 0.27147; + else if (pointIsCubeFaceCenter(*neighbor, side_length)) + weight = 0.39751; + else + libmesh_error_msg("We should never get here!"); + + node_approx += weight * (*neighbor); + } + + else + libmesh_error_msg("We should never get here!"); + + CPPUNIT_ASSERT(node.relative_fuzzy_equals(node_approx, 1e-3)); + + continue; + break; + } + default: { + libmesh_error_msg(node << " has unexpected number of neighbors (" + << neighbors.size() << ")"); + break; + } + } + } + CPPUNIT_ASSERT(node_distortion_is(node, false, 1e-2)); } @@ -988,6 +1249,38 @@ public: testVariationalSmoother(mesh, variational, PRISM21, true); } + void testVariationalTet4() + { + ReplicatedMesh mesh(*TestCommWorld); + VariationalMeshSmoother variational(mesh); + + testVariationalSmoother(mesh, variational, TET4); + } + + void testVariationalTet10() + { + ReplicatedMesh mesh(*TestCommWorld); + VariationalMeshSmoother variational(mesh); + + testVariationalSmoother(mesh, variational, TET10); + } + + void testVariationalTet14() + { + ReplicatedMesh mesh(*TestCommWorld); + VariationalMeshSmoother variational(mesh); + + testVariationalSmoother(mesh, variational, TET14); + } + + void testVariationalTet14MultipleSubdomains() + { + ReplicatedMesh mesh(*TestCommWorld); + VariationalMeshSmoother variational(mesh); + + testVariationalSmoother(mesh, variational, TET14, true); + } + void testVariationalMixed2D() { ReplicatedMesh mesh(*TestCommWorld); @@ -1003,6 +1296,7 @@ public: testVariationalSmootherRegression(mesh); } + #endif // LIBMESH_ENABLE_VSMOOTHER }; From 43a8c42bce27b9624610602f9199dde4f853bb26 Mon Sep 17 00:00:00 2001 From: Patrick Behne Date: Wed, 22 Oct 2025 09:17:04 -0600 Subject: [PATCH 3/8] Removed duplicate variables. --- src/systems/variational_smoother_system.C | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/systems/variational_smoother_system.C b/src/systems/variational_smoother_system.C index 33b1255a0a..8668230b50 100644 --- a/src/systems/variational_smoother_system.C +++ b/src/systems/variational_smoother_system.C @@ -896,6 +896,7 @@ VariationalSmootherSystem::get_target_elem(const ElemType & type) const auto ref_vol = target_elem->reference_elem()->volume(); // Update the nodes of the target element, depending on type + const Real sqrt_2 = std::sqrt(Real(2)); const Real sqrt_3 = std::sqrt(Real(3)); std::vector> owned_nodes; @@ -1028,7 +1029,6 @@ VariationalSmootherSystem::get_target_elem(const ElemType & type) // Solving for s: s = (3 sqrt(2) v)^(1/3), where v is the volume of the // non-optimal reference element. - const Real sqrt_2 = std::sqrt(Real(2)); // Side length that preserves the volume of the reference element const auto side_length = std::pow(3. * sqrt_2 * ref_vol, 1. / 3.); // Pyramid height with the property that all faces are equilateral triangles @@ -1100,8 +1100,6 @@ VariationalSmootherSystem::get_target_elem(const ElemType & type) // reference element is s = (12 * v / sqrt(2))^(1/3), where v is the volume // of the non-optimal reference element (i.e., a right tet). - const Real sqrt_2 = std::sqrt(Real(2)); - const Real sqrt_3 = std::sqrt(Real(3)); // Side length that preserves the volume of the reference element const auto side_length = std::pow(12. * ref_vol / sqrt_2, 1. / 3.); // tet height with the property that all faces are equilateral triangles From 5e574f7c52669d41cba61ecf302664f7783e17b1 Mon Sep 17 00:00:00 2001 From: Patrick Behne Date: Fri, 24 Oct 2025 08:57:48 -0600 Subject: [PATCH 4/8] Addressed Roy's review. --- include/mesh/mesh_tools.h | 13 ++++++ .../systems/variational_smoother_constraint.h | 13 ------ src/mesh/mesh_tools.C | 39 +++++++++++++++++ src/systems/variational_smoother_constraint.C | 43 +------------------ src/systems/variational_smoother_system.C | 4 +- tests/mesh/mesh_smoother_test.C | 22 +++++----- 6 files changed, 67 insertions(+), 67 deletions(-) diff --git a/include/mesh/mesh_tools.h b/include/mesh/mesh_tools.h index 9564378221..cf702fbbc7 100644 --- a/include/mesh/mesh_tools.h +++ b/include/mesh/mesh_tools.h @@ -347,6 +347,19 @@ void find_nodal_neighbors(const MeshBase & mesh, const std::unordered_map> & nodes_to_elem_map, std::vector & neighbors); +/** + * Given a mesh and a node in the mesh, the vector will be filled with + * every node directly attached to the given one. IF NO nodal neighbors are found, + * all nodes on the containing side are treated as neighbors. This is useful + * when the node does not lie on an edge, such as the central face node in + * HEX27, TET14, etc. elements. + */ +void find_nodal_or_face_neighbors( + const MeshBase & mesh, + const Node & node, + const std::unordered_map> & nodes_to_elem_map, + std::vector & neighbors); + /** * Given a mesh hanging_nodes will be filled with an associative array keyed off the * global id of all the hanging nodes in the mesh. It will hold an array of the diff --git a/include/systems/variational_smoother_constraint.h b/include/systems/variational_smoother_constraint.h index 3222ead25b..190318b5c5 100644 --- a/include/systems/variational_smoother_constraint.h +++ b/include/systems/variational_smoother_constraint.h @@ -504,19 +504,6 @@ class VariationalSmootherConstraint : public System::Constraint virtual ~VariationalSmootherConstraint() override; virtual void constrain() override; - - /** - * Given a mesh and a node in the mesh, the vector will be filled with - * every node directly attached to the given one. IF NO neighbors are found, - * all nodes on the containing side are treated as neighbors. This is useful - * when the node does not lie on an edge, such as the central face node in - * HEX27 elements. - */ - static void find_nodal_or_face_neighbors( - const MeshBase & mesh, - const Node & node, - const std::unordered_map> & nodes_to_elem_map, - std::vector & neighbors); }; diff --git a/src/mesh/mesh_tools.C b/src/mesh/mesh_tools.C index 405b5c4506..a2f91c539f 100644 --- a/src/mesh/mesh_tools.C +++ b/src/mesh/mesh_tools.C @@ -1061,6 +1061,45 @@ void find_nodal_neighbors(const MeshBase &, find_nodal_neighbors_helper(node.id(), node_to_elem_vec, neighbors); } +void find_nodal_or_face_neighbors( + const MeshBase & mesh, + const Node & node, + const std::unordered_map> & nodes_to_elem_map, + std::vector & neighbors) +{ + // Find all the nodal neighbors... that is the nodes directly connected + // to this node through one edge. + find_nodal_neighbors(mesh, node, nodes_to_elem_map, neighbors); + + // If no neighbors are found, use all nodes on the containing side as + // neighbors. + if (!neighbors.size()) + { + // Grab the element containing node + const auto * elem = libmesh_map_find(nodes_to_elem_map, node.id()).front(); + // Find the element side containing node + for (const auto &side : elem->side_index_range()) + { + const auto &nodes_on_side = elem->nodes_on_side(side); + const auto it = + std::find_if(nodes_on_side.begin(), nodes_on_side.end(), [&](auto local_node_id) { + return elem->node_id(local_node_id) == node.id(); + }); + + if (it != nodes_on_side.end()) + { + for (const auto &local_node_id : nodes_on_side) + // No need to add node itself as a neighbor + if (const auto *node_ptr = elem->node_ptr(local_node_id); + *node_ptr != node) + neighbors.push_back(node_ptr); + break; + } + } + } + libmesh_assert(neighbors.size()); +} + void find_hanging_nodes_and_parents(const MeshBase & mesh, diff --git a/src/systems/variational_smoother_constraint.C b/src/systems/variational_smoother_constraint.C index 44c942fd8c..f1b9b6d709 100644 --- a/src/systems/variational_smoother_constraint.C +++ b/src/systems/variational_smoother_constraint.C @@ -521,45 +521,6 @@ void VariationalSmootherConstraint::constrain_node_to_line(const Node & node, } } -void VariationalSmootherConstraint::find_nodal_or_face_neighbors( - const MeshBase & mesh, - const Node & node, - const std::unordered_map> & nodes_to_elem_map, - std::vector & neighbors) -{ - // Find all the nodal neighbors... that is the nodes directly connected - // to this node through one edge. - MeshTools::find_nodal_neighbors(mesh, node, nodes_to_elem_map, neighbors); - - // If no neighbors are found, use all nodes on the containing side as - // neighbors. - if (!neighbors.size()) - { - // Grab the element containing node - const auto * elem = libmesh_map_find(nodes_to_elem_map, node.id()).front(); - // Find the element side containing node - for (const auto &side : elem->side_index_range()) - { - const auto &nodes_on_side = elem->nodes_on_side(side); - const auto it = - std::find_if(nodes_on_side.begin(), nodes_on_side.end(), [&](auto local_node_id) { - return elem->node_id(local_node_id) == node.id(); - }); - - if (it != nodes_on_side.end()) - { - for (const auto &local_node_id : nodes_on_side) - // No need to add node itself as a neighbor - if (const auto *node_ptr = elem->node_ptr(local_node_id); - *node_ptr != node) - neighbors.push_back(node_ptr); - break; - } - } - } - libmesh_assert(neighbors.size()); -} - // Utility function to determine whether two nodes share a boundary ID. // The motivation for this is that a sliding boundary node on a triangular // element can have a neighbor boundary node in the same element that is not @@ -630,7 +591,7 @@ VariationalSmootherConstraint::get_neighbors_for_subdomain_constraint( // to this node through one edge, or if none exists, use other nodes on the // containing face std::vector neighbors; - VariationalSmootherConstraint::find_nodal_or_face_neighbors( + MeshTools::find_nodal_or_face_neighbors( mesh, node, nodes_to_elem_map, neighbors); // Each constituent set corresponds to neighbors sharing a face on the @@ -716,7 +677,7 @@ VariationalSmootherConstraint::get_neighbors_for_boundary_constraint( // to this node through one edge, or if none exists, use other nodes on the // containing face std::vector neighbors; - VariationalSmootherConstraint::find_nodal_or_face_neighbors( + MeshTools::find_nodal_or_face_neighbors( mesh, node, nodes_to_elem_map, neighbors); // Each constituent set corresponds to neighbors sharing a face on the diff --git a/src/systems/variational_smoother_system.C b/src/systems/variational_smoother_system.C index 8668230b50..5aff4ef695 100644 --- a/src/systems/variational_smoother_system.C +++ b/src/systems/variational_smoother_system.C @@ -1097,11 +1097,11 @@ VariationalSmootherSystem::get_target_elem(const ElemType & type) // For all faces to have side length s, the height of the apex node is // h = sqrt(2/3) * s. Then the volume is v = sqrt(2) * s^3 / 12. // Solving for s, the side length that will preserve the volume of the - // reference element is s = (12 * v / sqrt(2))^(1/3), where v is the volume + // reference element is s = (6 * sqrt(2) * v)^(1/3), where v is the volume // of the non-optimal reference element (i.e., a right tet). // Side length that preserves the volume of the reference element - const auto side_length = std::pow(12. * ref_vol / sqrt_2, 1. / 3.); + const auto side_length = std::pow(6. * sqrt_2 * ref_vol, 1 / Real(3)); // tet height with the property that all faces are equilateral triangles const auto target_height = sqrt_2 / sqrt_3 * side_length; diff --git a/tests/mesh/mesh_smoother_test.C b/tests/mesh/mesh_smoother_test.C index 41334acdfe..835944950f 100644 --- a/tests/mesh/mesh_smoother_test.C +++ b/tests/mesh/mesh_smoother_test.C @@ -7,7 +7,6 @@ #include #include #include -#include #include #include #include @@ -366,7 +365,7 @@ public: // This coordinante does not occur at a sub-cube face or center continue; - if ((unsigned int)std::round(num_half_lengths) % 2) + if (std::lround(num_half_lengths) % 2) // An odd number of half lengths means point(d) is at a sub-cube center ++num_centered; else @@ -750,11 +749,12 @@ public: // smoothed to the actual midpoints. else if (type_is_tet && !elem->is_vertex(local_node_id)) { + const Real tol = 1e-5; // We have a non-vertex node. Determine what "type" of // midpoint node with respect to the mesh geometry. // First, get the nodes that neighbor this node std::vector neighbors; - VariationalSmootherConstraint::find_nodal_or_face_neighbors( + MeshTools::find_nodal_or_face_neighbors( mesh, node, nodes_to_elem_map, neighbors); switch (neighbors.size()) @@ -784,14 +784,14 @@ public: { const Real x = (type == TET10) ? 0.42895 : 0.41486; CPPUNIT_ASSERT(node.relative_fuzzy_equals( - other + x * (cube_center - other), 1e-3)); + other + x * (cube_center - other), tol)); } else if (pointIsCubeVertex(other, side_length)) { const Real x = (type == TET10) ? 0.55389 : 0.58094; CPPUNIT_ASSERT(node.relative_fuzzy_equals( - other + x * (cube_center - other), 1e-3)); + other + x * (cube_center - other), tol)); } } @@ -813,7 +813,7 @@ public: // Due to symmetry, this type of midpoint is // the geometric midpoint. Let the // node_distortion_is function check it - CPPUNIT_ASSERT(node_distortion_is(node, false, 1e-2)); + CPPUNIT_ASSERT(node_distortion_is(node, false)); else { @@ -827,7 +827,7 @@ public: is_0_cube_face_center ? *neighbors[0] : *neighbors[1]; const Real x = (type == TET10) ? 0.61299 : 0.65126; CPPUNIT_ASSERT(node.relative_fuzzy_equals( - cube_vertex + x * (cube_face_center - cube_vertex), 1e-3)); + cube_vertex + x * (cube_face_center - cube_vertex), tol)); } } @@ -934,7 +934,7 @@ public: else libmesh_error_msg("We should never get here!"); - CPPUNIT_ASSERT(node.relative_fuzzy_equals(node_approx, 1e-3)); + CPPUNIT_ASSERT(node.relative_fuzzy_equals(node_approx, tol)); continue; break; @@ -944,10 +944,10 @@ public: << neighbors.size() << ")"); break; } - } - } + } // switch (neighbors.size()) + } // if (type_is_tet) - CPPUNIT_ASSERT(node_distortion_is(node, false, 1e-2)); + CPPUNIT_ASSERT(node_distortion_is(node, false)); } } From 6211c747db1d66ba68f2151a9211aa222168bc1d Mon Sep 17 00:00:00 2001 From: Patrick Behne Date: Fri, 24 Oct 2025 09:13:35 -0600 Subject: [PATCH 5/8] Addressed John's review. --- src/systems/variational_smoother_system.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/systems/variational_smoother_system.C b/src/systems/variational_smoother_system.C index 5aff4ef695..3d362d3351 100644 --- a/src/systems/variational_smoother_system.C +++ b/src/systems/variational_smoother_system.C @@ -1101,7 +1101,7 @@ VariationalSmootherSystem::get_target_elem(const ElemType & type) // of the non-optimal reference element (i.e., a right tet). // Side length that preserves the volume of the reference element - const auto side_length = std::pow(6. * sqrt_2 * ref_vol, 1 / Real(3)); + const auto side_length = std::cbrt(6. * sqrt_2 * ref_vol); // tet height with the property that all faces are equilateral triangles const auto target_height = sqrt_2 / sqrt_3 * side_length; From f56be277097ca139b255bfecbb2ab5301e738ff3 Mon Sep 17 00:00:00 2001 From: Patrick Behne Date: Fri, 24 Oct 2025 09:21:41 -0600 Subject: [PATCH 6/8] Tightened tolerances in pointIsCube* helper functions. --- tests/mesh/mesh_smoother_test.C | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/mesh/mesh_smoother_test.C b/tests/mesh/mesh_smoother_test.C index 835944950f..8dcf5615fc 100644 --- a/tests/mesh/mesh_smoother_test.C +++ b/tests/mesh/mesh_smoother_test.C @@ -378,21 +378,21 @@ public: // Helper function to determine whether a given point lies at the center of // a sub-cube - bool pointIsCubeCenter(const Point & point, const Real & side_length, const Real & tol = 1e-3) + bool pointIsCubeCenter(const Point & point, const Real & side_length, const Real & tol = TOLERANCE) { return numCenteredAndFacedDimensions(point, side_length, tol).first == 3; } // Helper function to determine whether a given point lies at the vertex of // a sub-cube - bool pointIsCubeVertex(const Point & point, const Real & side_length, const Real & tol = 1e-3) + bool pointIsCubeVertex(const Point & point, const Real & side_length, const Real & tol = TOLERANCE) { return numCenteredAndFacedDimensions(point, side_length, tol).second == 3; } // Helper function to determine whether a given point lies at the center of // a sub-cube face - bool pointIsCubeFaceCenter(const Point & point, const Real & side_length, const Real & tol = 1e-3) + bool pointIsCubeFaceCenter(const Point & point, const Real & side_length, const Real & tol = TOLERANCE) { const auto result = numCenteredAndFacedDimensions(point, side_length, tol); return (result.first == 2) && (result.second == 1); From bd4defd8ec46019f84a37596eee19cde2e0d6d65 Mon Sep 17 00:00:00 2001 From: Patrick Behne Date: Fri, 24 Oct 2025 11:43:33 -0600 Subject: [PATCH 7/8] Tightened test tolerances. --- tests/mesh/mesh_smoother_test.C | 24 +++++++++++++----------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/tests/mesh/mesh_smoother_test.C b/tests/mesh/mesh_smoother_test.C index 8dcf5615fc..156937b134 100644 --- a/tests/mesh/mesh_smoother_test.C +++ b/tests/mesh/mesh_smoother_test.C @@ -20,6 +20,7 @@ #include "test_comm.h" #include "libmesh_cppunit.h" +#include #include @@ -646,6 +647,7 @@ public: } smoother.smooth(); + mesh.write("smoothed.e"); if (type_is_tri) { @@ -749,7 +751,7 @@ public: // smoothed to the actual midpoints. else if (type_is_tet && !elem->is_vertex(local_node_id)) { - const Real tol = 1e-5; + const Real tol = TOLERANCE; // We have a non-vertex node. Determine what "type" of // midpoint node with respect to the mesh geometry. // First, get the nodes that neighbor this node @@ -782,14 +784,14 @@ public: // a sub-cube vertex. Determine which one. if (pointIsCubeFaceCenter(other, side_length)) { - const Real x = (type == TET10) ? 0.42895 : 0.41486; + const Real x = (type == TET10) ? 0.42895041 : 0.41486385; CPPUNIT_ASSERT(node.relative_fuzzy_equals( other + x * (cube_center - other), tol)); } else if (pointIsCubeVertex(other, side_length)) { - const Real x = (type == TET10) ? 0.55389 : 0.58094; + const Real x = (type == TET10) ? 0.55388920 : 0.58093516; CPPUNIT_ASSERT(node.relative_fuzzy_equals( other + x * (cube_center - other), tol)); } @@ -825,7 +827,7 @@ public: is_0_cube_vertex ? *neighbors[0] : *neighbors[1]; const auto & cube_face_center = is_0_cube_face_center ? *neighbors[0] : *neighbors[1]; - const Real x = (type == TET10) ? 0.61299 : 0.65126; + const Real x = (type == TET10) ? 0.61299101 : 0.65125580; CPPUNIT_ASSERT(node.relative_fuzzy_equals( cube_vertex + x * (cube_face_center - cube_vertex), tol)); } @@ -885,9 +887,9 @@ public: { Real weight; if (pointIsCubeVertex(*neighbor, side_length)) - weight = 0.30601; + weight = 0.30600747; else if (pointIsCubeCenter(*neighbor, side_length)) - weight = 0.38799; + weight = 0.38798506; else libmesh_error_msg("We should never get here!"); @@ -900,9 +902,9 @@ public: { Real weight; if (pointIsCubeVertex(*neighbor, side_length)) - weight = 0.28078; + weight = 0.28078090; else if (pointIsCubeFaceCenter(*neighbor, side_length)) - weight = 0.43843; + weight = 0.43843820; else libmesh_error_msg("We should never get here!"); @@ -920,11 +922,11 @@ public: { Real weight; if (pointIsCubeCenter(*neighbor, side_length)) - weight = 0.33102; + weight = 0.33102438; else if (pointIsCubeVertex(*neighbor, side_length)) - weight = 0.27147; + weight = 0.27147230; else if (pointIsCubeFaceCenter(*neighbor, side_length)) - weight = 0.39751; + weight = 0.39750332; else libmesh_error_msg("We should never get here!"); From 5376a24af0bd4166d0d5cda764dd63d14ec4f6c4 Mon Sep 17 00:00:00 2001 From: Patrick Behne Date: Fri, 24 Oct 2025 11:45:19 -0600 Subject: [PATCH 8/8] Removed debug statement. --- tests/mesh/mesh_smoother_test.C | 2 -- 1 file changed, 2 deletions(-) diff --git a/tests/mesh/mesh_smoother_test.C b/tests/mesh/mesh_smoother_test.C index 156937b134..310a0c0ada 100644 --- a/tests/mesh/mesh_smoother_test.C +++ b/tests/mesh/mesh_smoother_test.C @@ -20,7 +20,6 @@ #include "test_comm.h" #include "libmesh_cppunit.h" -#include #include @@ -647,7 +646,6 @@ public: } smoother.smooth(); - mesh.write("smoothed.e"); if (type_is_tri) {