Skip to content

Variational Smoothing for Pyramids #4232

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 6 commits into from
Aug 28, 2025
Merged

Variational Smoothing for Pyramids #4232

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
9fbf15b
Added target elements for pyramids. Ref #4082.
pbehne Aug 14, 2025
9014cfd
Added tests for pyramids.
pbehne Aug 14, 2025
72070b5
Renamed distortion_is --> node_distortion_is.
pbehne Aug 25, 2025
3b307b3
Update target pyramid to square-base pyramid with equilateral triangu…
pbehne Aug 25, 2025
d04035f
Reduced number of elements for higher order elem types to keep meshes…
pbehne Aug 27, 2025
290bfcc
Reduced tangling damping for Hex27.
pbehne Aug 27, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
71 changes: 70 additions & 1 deletion src/systems/variational_smoother_system.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -927,6 +927,75 @@ VariationalSmootherSystem::get_target_elem(const ElemType & type)
} // if Tri


// Elems deriving from Pyramid
else if (type_str.compare(0, 7, "PYRAMID") == 0)
{

// The target element is a pyramid with an square base and
// equilateral triangular sides with volume equal to the volume of the
// reference element.

// A pyramid with square base sidelength s and equilateral triangular
// sides has height h = s / sqrt(2).
// The volume is v = s^2 h / 3 = s^3 / ( 3 sqrt(2)).
// 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
const auto target_height = side_length / sqrt_2;

const auto & s = side_length;
const auto & h = target_height;

// 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(s, s, 0.), 2));
owned_nodes.emplace_back(Node::build(Point(0., s, 0.), 3));
owned_nodes.emplace_back(Node::build(Point(0.5 * s, 0.5 * s, h), 4));

if (type == PYRAMID13 || type == PYRAMID14 || type == PYRAMID18)
{
const auto & on = owned_nodes;
// Define the edge midpoint nodes of the pyramid

// Base node to base node midpoint nodes
owned_nodes.emplace_back(Node::build((*on[0] + *on[1]) / 2., 5));
owned_nodes.emplace_back(Node::build((*on[1] + *on[2]) / 2., 6));
owned_nodes.emplace_back(Node::build((*on[2] + *on[3]) / 2., 7));
owned_nodes.emplace_back(Node::build((*on[3] + *on[0]) / 2., 8));

// Base node to apex node midpoint nodes
owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[4]) / 2.), 9));
owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[4]) / 2.), 10));
owned_nodes.emplace_back(Node::build(Point((*on[2] + *on[4]) / 2.), 11));
owned_nodes.emplace_back(Node::build(Point((*on[3] + *on[4]) / 2.), 12));

if (type == PYRAMID14 || type == PYRAMID18)
{
// Define the square face midpoint node of the pyramid
owned_nodes.emplace_back(
Node::build(Point((*on[0] + *on[1] + *on[2] + *on[3]) / 4.), 13));

if (type == PYRAMID18)
{
// Define the triangular face nodes
owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[1] + *on[4]) / 3.), 14));
owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[2] + *on[4]) / 3.), 15));
owned_nodes.emplace_back(Node::build(Point((*on[2] + *on[3] + *on[4]) / 3.), 16));
owned_nodes.emplace_back(Node::build(Point((*on[3] + *on[0] + *on[4]) / 3.), 17));
}
}
}

else if (type != PYRAMID5)
libmesh_error_msg("Unsupported pyramid element: " << type_str);

} // if Pyramid

// Set the target_elem equal to the reference elem
else
for (const auto & node : target_elem->reference_elem()->node_ref_range())
Expand All @@ -936,7 +1005,7 @@ VariationalSmootherSystem::get_target_elem(const ElemType & type)
for (const auto & node_ptr : owned_nodes)
target_elem->set_node(node_ptr->id(), node_ptr.get());

libmesh_assert(relative_fuzzy_equals(target_elem->volume(), ref_vol));
libmesh_assert(relative_fuzzy_equals(target_elem->volume(), ref_vol, TOLERANCE));

return std::make_pair(std::move(target_elem), std::move(owned_nodes));
}
Expand Down
248 changes: 178 additions & 70 deletions tests/mesh/mesh_smoother_test.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -97,7 +97,7 @@ private:
}
const auto delta = (pow<3>(xi) - xi) * modulation;
// Check for delta = 0 to make sure we perturb every point
output(i) = (delta > TOLERANCE) ? p(i) + delta : 1.05 * p(i);
output(i) = (delta > TOLERANCE) ? p(i) + delta : 1.03 * p(i);
}
else
output(i) = p(i); // dimension on boundary remains unchanged
Expand Down Expand Up @@ -154,7 +154,6 @@ class ParallelogramToSquare : public FunctionBase<Real>
output(2) = 0;
}
};

}

using namespace libMesh;
Expand Down Expand Up @@ -189,6 +188,12 @@ public:
CPPUNIT_TEST(testVariationalHex27);
CPPUNIT_TEST(testVariationalHex27Tangled);
CPPUNIT_TEST(testVariationalHex27MultipleSubdomains);

CPPUNIT_TEST(testVariationalPyramid5);
CPPUNIT_TEST(testVariationalPyramid13);
CPPUNIT_TEST(testVariationalPyramid14);
CPPUNIT_TEST(testVariationalPyramid18);
CPPUNIT_TEST(testVariationalPyramid18MultipleSubdomains);
#endif // LIBMESH_ENABLE_VSMOOTHER
#endif

Expand Down Expand Up @@ -282,8 +287,11 @@ public:
const auto * ref_elem = &(ReferenceElem::get(type));
const auto dim = ref_elem->dim();
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;

unsigned int n_elems_per_side = 4;
// 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];

switch (dim)
{
Expand Down Expand Up @@ -328,7 +336,7 @@ public:

// Define p1 and p2 given the integer indices
// Undistorted elem side length
const Real dr = 1. / n_elems_per_side;
const Real dr = 1. / (n_elems_per_side * Elem::type_to_default_order_map[type]);
const Point p1 = Point(dr, dim > 1 ? dr : 0, dim > 2 ? dr : 0);
const Point p2 = Point(2. * dr, dim > 1 ? dr : 0, dim > 2 ? 2. * dr : 0);

Expand Down Expand Up @@ -410,71 +418,70 @@ public:
distorted_subdomain_volumes[elem->subdomain_id()] += elem->volume();
}


// Get the mesh order
const auto &elem_orders = mesh.elem_default_orders();
libmesh_error_msg_if(
elem_orders.size() != 1,
"The variational smoother cannot be used for mixed-order meshes!");
const auto fe_order = *elem_orders.begin();

// Function to assert the distortion is as expected
auto distortion_is = [
&n_elems_per_side, &dim, &boundary_info, &fe_order
](const Node &node, bool distortion, Real distortion_tol = TOLERANCE)
{
// Get boundary ids associated with the node
std::vector<boundary_id_type> boundary_ids;
boundary_info.boundary_ids(&node, boundary_ids);

// This tells us what type of node we are: internal, sliding, or fixed
const auto num_dofs = dim - boundary_ids.size();
/*
* The following cases of num_dofs are possible, ASSUMING all boundaries
* are non-overlapping
* 3D: 3-0, 3-1, 3-2, 3-3
* = 3 2 1 0
* internal sliding, sliding, fixed
* 2D: 2-0, 2-1, 2-2
* = 2 1 0
* internal sliding, fixed
* 1D: 1-0, 1-1
* = 1 0
* internal fixed
*
* We expect that R is an integer in [0, n_elems_per_side] for
* num_dofs of the node's cooridinantes, while the remaining coordinates
* are fixed to the boundary with value 0 or 1. In other words, at LEAST
* dim - num_dofs coordinantes should be 0 or 1.
*/

std::size_t num_zero_or_one = 0;

bool distorted = false;
for (const auto d : make_range(dim))
{
const Real r = node(d);
const Real R = r * n_elems_per_side * fe_order;

const bool d_distorted = std::abs(R - std::round(R)) > distortion_tol;
distorted |= d_distorted;
num_zero_or_one +=
(absolute_fuzzy_equals(r, 0.) || absolute_fuzzy_equals(r, 1.));
}

CPPUNIT_ASSERT_GREATEREQUAL(dim - num_dofs, num_zero_or_one);

// We can never expect a fixed node to be distorted
if (num_dofs == 0)
// if (num_dofs < dim)
return true;
return distorted == distortion;

// Get the mesh order
const auto & elem_orders = mesh.elem_default_orders();
libmesh_error_msg_if(elem_orders.size() != 1,
"The variational smoother cannot be used for mixed-order meshes!");
// 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);

// Function to assert the node distortion is as expected
auto node_distortion_is = [&n_elems_per_side, &dim, &boundary_info, &scale_factor](
const Node & node, bool distortion, Real distortion_tol = TOLERANCE) {
// Get boundary ids associated with the node
std::vector<boundary_id_type> boundary_ids;
boundary_info.boundary_ids(&node, boundary_ids);

// This tells us what type of node we are: internal, sliding, or fixed
const auto num_dofs = dim - boundary_ids.size();
/*
* The following cases of num_dofs are possible, ASSUMING all boundaries
* are non-overlapping
* 3D: 3-0, 3-1, 3-2, 3-3
* = 3 2 1 0
* internal sliding, sliding, fixed
* 2D: 2-0, 2-1, 2-2
* = 2 1 0
* internal sliding, fixed
* 1D: 1-0, 1-1
* = 1 0
* internal fixed
*
* We expect that R is an integer in [0, n_elems_per_side] for
* num_dofs of the node's cooridinantes, while the remaining coordinates
* are fixed to the boundary with value 0 or 1. In other words, at LEAST
* dim - num_dofs coordinantes should be 0 or 1.
*/

std::size_t num_zero_or_one = 0;

bool distorted = false;
for (const auto d : make_range(dim))
{
const Real r = node(d);
const Real R = r * n_elems_per_side * scale_factor;
CPPUNIT_ASSERT_GREATER(-distortion_tol * distortion_tol, r);
CPPUNIT_ASSERT_GREATER(-distortion_tol * distortion_tol, 1 - r);

const bool d_distorted = std::abs(R - std::round(R)) > distortion_tol;
distorted |= d_distorted;
num_zero_or_one += (absolute_fuzzy_equals(r, 0.) || absolute_fuzzy_equals(r, 1.));
}

CPPUNIT_ASSERT_GREATEREQUAL(dim - num_dofs, num_zero_or_one);

// We can never expect a fixed node to be distorted
if (num_dofs == 0)
// if (num_dofs < dim)
return true;
return distorted == distortion;
};

// Make sure our DistortHyperCube transformation has distorted the mesh
for (auto node : mesh.node_ptr_range())
CPPUNIT_ASSERT(distortion_is(*node, true));
CPPUNIT_ASSERT(node_distortion_is(*node, true));

// Transform the square mesh of triangles to a parallelogram mesh of
// triangles. This will allow the Variational Smoother to smooth the mesh
Expand Down Expand Up @@ -511,13 +518,74 @@ public:
const auto & subdomain_id = pair.first;
CPPUNIT_ASSERT(
relative_fuzzy_equals(libmesh_map_find(distorted_subdomain_volumes, subdomain_id),
libmesh_map_find(smoothed_subdomain_volumes, subdomain_id)));
libmesh_map_find(smoothed_subdomain_volumes, subdomain_id),
TOLERANCE));
}
}
else
// Make sure we're not too distorted anymore
for (auto node : mesh.node_ptr_range())
CPPUNIT_ASSERT(distortion_is(*node, false, 1e-2));
{
// Make sure we're not too distorted anymore
std::set<dof_id_type> nodes_checked;
for (const auto * elem : mesh.active_element_ptr_range())
{
for (const auto local_node_id : make_range(elem->n_nodes()))
{
const auto & node = elem->node_ref(local_node_id);
if (nodes_checked.find(node.id()) != nodes_checked.end())
continue;

nodes_checked.insert(node.id());

// Check for special case where pyramidal base-to-apex midpoint
// nodes are not smoothed to the actual midpoints.
if (type_is_pyramid)
{
if (local_node_id > 8 && local_node_id < 13)
{
// Base-Apex midpoint nodes
//
// Due to the nature of the pyramidal target element and the
// cubic nature of the test mesh, for a dilation weight o
// 0.5, smoothed base-apex midpoint nodes are smoothed to
// the point base + x * (apex - base) instead of
// base + 0.5 (apex - base), where x is in (0,1) and
// depends on the number of nodes in the element.
// We hard-code this check below.

const auto & base = elem->node_ref(local_node_id - 9);
const auto & apex = elem->node_ref(4);
const Real x = (type == PYRAMID18) ? 0.569332 : 0.549876;

CPPUNIT_ASSERT(node.relative_fuzzy_equals(base + x * (apex - base), 1e-3));
continue;
}
else if (local_node_id > 13)
{
// Triangular face midpoint nodes
//
// Due to the nature of the pyramidal target element and the
// cubic nature of the test mesh, for a dilation weight o
// 0.5, smoothed triangular face midpoint nodes are
// smoothed a weighted average of the constituent
// vertices instead of an unweighted average.
// We hard-code this check below.
const auto & base1 = elem->node_ref(local_node_id - 14);
const auto & base2 = elem->node_ref((local_node_id - 13) % 4);
const auto & apex = elem->node_ref(4);

const auto node_approx = (0.3141064847 * base1 +
0.3141064847 * base2 +
0.3717870306 * apex);
CPPUNIT_ASSERT(node.relative_fuzzy_equals(node_approx, 1e-3));
continue;
}
}

CPPUNIT_ASSERT(node_distortion_is(node, false, 1e-2));

}
}
}
}

void testLaplaceQuad()
Expand Down Expand Up @@ -646,7 +714,47 @@ public:
ReplicatedMesh mesh(*TestCommWorld);
VariationalMeshSmoother variational(mesh);

testVariationalSmoother(mesh, variational, HEX27, false, true, 0.75);
testVariationalSmoother(mesh, variational, HEX27, false, true, 0.55);
}

void testVariationalPyramid5()
{
ReplicatedMesh mesh(*TestCommWorld);
VariationalMeshSmoother variational(mesh);

testVariationalSmoother(mesh, variational, PYRAMID5);
}

void testVariationalPyramid13()
{
ReplicatedMesh mesh(*TestCommWorld);
VariationalMeshSmoother variational(mesh);

testVariationalSmoother(mesh, variational, PYRAMID13);
}

void testVariationalPyramid14()
{
ReplicatedMesh mesh(*TestCommWorld);
VariationalMeshSmoother variational(mesh);

testVariationalSmoother(mesh, variational, PYRAMID14);
}

void testVariationalPyramid18()
{
ReplicatedMesh mesh(*TestCommWorld);
VariationalMeshSmoother variational(mesh);

testVariationalSmoother(mesh, variational, PYRAMID18);
}

void testVariationalPyramid18MultipleSubdomains()
{
ReplicatedMesh mesh(*TestCommWorld);
VariationalMeshSmoother variational(mesh);

testVariationalSmoother(mesh, variational, PYRAMID18, true);
}
#endif // LIBMESH_ENABLE_VSMOOTHER
};
Expand Down