tboudreaux/MeanField
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Audit physical compensation artifacts and verify quartic affine displacement representation

Browse Source
This commit is contained in:
Emily Boudreaux
2026-09-18 13:31:14 -04:00
parent d79db9445b
commit 2c4a5e05a6
2 changed files with 916 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

117
tests/fem/reference_tables.cpp
View File

@@ -556,3 +556,120 @@ TEST_CASE("Tensor Gradient Transposes Match Dense Assembly And Adjoint Products"
mfem::IntRules.Get(mfem::Geometry::CUBE, 2*order+3));
}
}
TEST_CASE(
"Quartic Curved Hex Geometry Distinguishes Q3 And Q4 Affine Physical Displacements",
tags::unit &tags::geometry &tags::quadrature
) {
constexpr int dimension = 3;
constexpr double warpAmplitude = 0.2;
constexpr double dilation = 0.2;
constexpr double exactDeterminant = (1.0 + dilation) * (1.0 + dilation) * (1.0 + dilation);
int ranks = 1;
MPI_Comm_size(MPI_COMM_WORLD, &ranks);
// Subdivide only y so every rank can retain the full quartic x variation.
// Its coordinate Jacobian is 1 + 2*a*x*(1-x)*(1-2*x) >= 1-a/2 > 0.
auto serial = mfem::Mesh::MakeCartesian3D(1, ranks, 1, mfem::Element::HEXAHEDRON);
serial.SetCurvature(4, false, dimension, mfem::Ordering::byNODES);
mfem::VectorFunctionCoefficient warp(dimension, [](const mfem::Vector &x, mfem::Vector &y) {
y = x;
const double bubble = x(0) * (1.0 - x(0));
y(0) += warpAmplitude * bubble * bubble;
});
serial.Transform(warp);
mfem::ParMesh mesh(MPI_COMM_WORLD, serial);
REQUIRE(mesh.GetNodes() != nullptr);
CHECK(mesh.GetNodes()->FESpace()->GetMaxElementOrder() == 4);
// u(x)=s*x pulls back to the quartic coordinate polynomial. This tests the
// ambient vector space, not the production surface prescription's r^2 lift.
mfem::VectorFunctionCoefficient affineDisplacement(
dimension, [](const mfem::Vector &x, mfem::Vector &u) {
u = x;
u *= dilation;
}
);
const auto &rule = mfem::IntRules.Get(mfem::Geometry::CUBE, 24);
for (const int order : {3, 4}) {
CAPTURE(order);
mfem::H1_FECollection collection(order, dimension, mfem::BasisType::GaussLobatto);
mfem::ParFiniteElementSpace space(&mesh, &collection, dimension, mfem::Ordering::byNODES);
CHECK(collection.FiniteElementForGeometry(mfem::Geometry::CUBE)->GetDof() ==
(order + 1) * (order + 1) * (order + 1));
mfem::ParGridFunction field(&space);
field.ProjectCoefficient(affineDisplacement);
mfem::Vector state;
field.GetTrueDofs(state);
field.SetFromTrueDofs(state);
// Accumulate rather than assert per quadrature point, and reject
// non-finite samples explicitly so a max reduction cannot hide them.
enum Error { coordinates, coordinateJacobian, values, gradients, determinants, count };
std::array<double, count> localErrors{};
int localFinite = 1;
double localGeometryMinimum = std::numeric_limits<double>::infinity();
auto recordError = [&](const Error error, const double value) {
if (!std::isfinite(value)) localFinite = 0;
else localErrors[error] = std::max(localErrors[error], std::abs(value));
};
mfem::Vector position(dimension), value(dimension);
mfem::DenseMatrix gradient(dimension);
for (int element = 0; element < mesh.GetNE(); ++element) {
auto &transformation = *mesh.GetElementTransformation(element);
for (int point = 0; point < rule.GetNPoints(); ++point) {
const auto &ip = rule.IntPoint(point);
transformation.SetIntPoint(&ip);
transformation.Transform(ip, position);
const double bubble = ip.x * (1.0 - ip.x);
recordError(coordinates, position(0) - (ip.x + warpAmplitude * bubble * bubble));
const double expectedJacobian = 1.0 + 2.0 * warpAmplitude * bubble * (1.0 - 2.0 * ip.x);
const double coordinateDeterminant = ranks * transformation.Jacobian().Det();
recordError(coordinateJacobian, coordinateDeterminant - expectedJacobian);
if (std::isfinite(coordinateDeterminant))
localGeometryMinimum = std::min(localGeometryMinimum, coordinateDeterminant);
field.GetVectorValue(transformation, ip, value);
field.GetVectorGradient(transformation, gradient);
for (int component = 0; component < dimension; ++component) {
recordError(values, value(component) - dilation * position(component));
for (int axis = 0; axis < dimension; ++axis)
recordError(gradients, gradient(component, axis) -
(component == axis ? dilation : 0.0));
}
for (int component = 0; component < dimension; ++component) gradient(component, component) += 1.0;
recordError(determinants, gradient.Det() - exactDeterminant);
}
}
std::array<double, count> errors{};
int finite = 0;
double geometryMinimum = 0.0;
MPI_Allreduce(localErrors.data(), errors.data(), count, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
MPI_Allreduce(&localFinite, &finite, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD);
MPI_Allreduce(&localGeometryMinimum, &geometryMinimum, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
REQUIRE(finite == 1);
CHECK(errors[coordinates] < 2.0e-13);
CHECK(errors[coordinateJacobian] < 2.0e-12);
CHECK(geometryMinimum >= 1.0 - warpAmplitude / 2.0 - 2.0e-12);
mean_field::fem::PreparedGeometryInspector inspector;
const double fastMinimum = inspector.Inspect(space, state);
const double referenceMinimum = inspector.Inspect(space, state, true);
REQUIRE(std::isfinite(fastMinimum));
REQUIRE(std::isfinite(referenceMinimum));
CHECK(fastMinimum > 0.0);
CHECK(std::abs(fastMinimum - referenceMinimum) < 2.0e-12);
CHECK(inspector.GetTensorElementCount() == mesh.GetNE());
if (order == 4) {
CHECK(errors[values] < 2.0e-13);
CHECK(errors[gradients] < 2.0e-12);
CHECK(errors[determinants] < 5.0e-12);
CHECK(std::abs(fastMinimum - exactDeterminant) < 5.0e-12);
} else {
CHECK(errors[values] > 1.0e-5);
CHECK(errors[gradients] > 1.0e-4);
CHECK(errors[determinants] > 1.0e-4);
CHECK(std::abs(fastMinimum - exactDeterminant) > 1.0e-4);
}
}
}