From 7fadd48e1b296327c89436ca7df24aec78d7b015 Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Fri, 9 Oct 2026 01:18:20 +0900 Subject: [PATCH] fix(elements): reject rounded loss of Hex8 bulk stiffness --- src/fesa/elements/hex8_element.cpp | 10 ++ tests/unit/elements/hex8_element_test.cpp | 176 ++++++++++++++++++++ tests/unit/results/result_recovery_test.cpp | 48 ++++++ 3 files changed, 234 insertions(+) diff --git a/src/fesa/elements/hex8_element.cpp b/src/fesa/elements/hex8_element.cpp index 6a4dfe9..a42c897 100644 --- a/src/fesa/elements/hex8_element.cpp +++ b/src/fesa/elements/hex8_element.cpp @@ -92,6 +92,16 @@ Result Hex8Element::Create(EntityIndex index, return Result::Failure( candidate.Failure("Hex8 Lame moduli must remain finite with " "representable positive shear modulus.")); + // Validate the rounded moduli actually used by the kernel, not merely the + // positive exact-real E/nu formula: cancellation near nu=-1 can erase bulk. + const double bulk = candidate.lambda_ + (2.0 * candidate.mu_) / 3.0; + const double volumetric_stress = + 3.0 * candidate.lambda_ + 2.0 * candidate.mu_; + if (!std::isfinite(bulk) || bulk <= 0.0 || + !std::isfinite(volumetric_stress) || volumetric_stress <= 0.0) + return Result::Failure( + candidate.Failure("Hex8 rounded Lame moduli require finite " + "representable positive bulk stiffness.")); if (density && (!std::isfinite(*density) || *density <= 0.0)) return Result::Failure(candidate.Failure( "Hex8 reference density must be finite and positive.")); diff --git a/tests/unit/elements/hex8_element_test.cpp b/tests/unit/elements/hex8_element_test.cpp index 5206693..3aa7f90 100644 --- a/tests/unit/elements/hex8_element_test.cpp +++ b/tests/unit/elements/hex8_element_test.cpp @@ -9,6 +9,7 @@ #include #include #include +#include namespace { using Coordinates = fesa::Hex8Element::Coordinates; @@ -298,6 +299,7 @@ TEST(Hex8Element, GeometryMaterialDensityScalingAndCoordinateCovariance) { const auto k = element.Value().ComputeStiffness().Value().values; const auto m = element.Value().ComputeMass().Value().values; const auto kr = rotated_element.Value().ComputeStiffness().Value().values; + const auto mr = rotated_element.Value().ComputeMass().Value().values; const auto ks = scaled.Value().ComputeStiffness().Value().values; const auto ms = scaled.Value().ComputeMass().Value().values; Tensor r{}; @@ -316,6 +318,9 @@ TEST(Hex8Element, GeometryMaterialDensityScalingAndCoordinateCovariance) { for (std::size_t n = 0; n < 3; ++n) oracle += r[i][l] * k(3 * a + l, 3 * b + n) * r[j][n]; EXPECT_NEAR(kr(3 * a + i, 3 * b + j), oracle, 3e-12); + // Isotropic translation mass is unchanged by a rigid coordinate + // rotation; this equals blockwise R*M*R^T in global components. + EXPECT_NEAR(mr(3 * a + i, 3 * b + j), m(3 * a + i, 3 * b + j), 3e-14); EXPECT_NEAR(ks(3 * a + i, 3 * b + j), 7.5 * k(3 * a + i, 3 * b + j), 3e-11); EXPECT_NEAR(ms(3 * a + i, 3 * b + j), @@ -614,4 +619,175 @@ TEST(Hex8Element, RejectsUnrepresentableShearAndMassInsteadOfReturningZeros) { ASSERT_TRUE(element.HasValue()); EXPECT_FALSE(element.Value().ComputeMass().HasValue()); } +TEST(Hex8Element, RejectsRoundedLossOfPositiveVolumetricModulus) { + // Binary64 Lame arithmetic at this valid nu yields lambda+2mu/3 == 0 + // even though exact E/[3(1-2nu)] is approximately1/9, which is positive. + const double nu = std::nextafter(-1.0, 0.0); + auto element = Make(Brick(), 1.0, nu); + ASSERT_FALSE(element.HasValue()); + EXPECT_NE(element.GetStatus().Diagnostics()[0].message.find("bulk"), + std::string::npos); +} +TEST(Hex8Element, OrdinaryNegativePoissonRatioRetainsPositiveVolumetricEnergy) { + auto element = Make(Brick(), 120.0, -.5); + ASSERT_TRUE(element.HasValue()); + const auto u = Affine( + Brick(), Tensor{{{1.01, 0.0, 0.0}, {0.0, 1.01, 0.0}, {0.0, 0.0, 1.01}}}); + const auto response = element.Value().EvaluateLinear(u); + ASSERT_TRUE(response.HasValue()); + const double bulk = 120.0 / (3.0 * (1.0 - 2.0 * (-.5))); + const double trace = .03; + EXPECT_NEAR(response.Value().strain_energy, 48.0 * .5 * bulk * trace * trace, + 2e-14); + for (const auto& point : response.Value().points) + for (std::size_t i = 0; i < 3; ++i) + EXPECT_NEAR(point.cauchy_stress[i], bulk * trace, 2e-14); +} + +TEST(Hex8Element, FiniteInputsThatOverflowDerivedArithmeticFailClosed) { + const double huge = std::numeric_limits::max(); + auto bad_material = Make(Brick(), huge, .49); + ASSERT_FALSE(bad_material.HasValue()); + EXPECT_NE(bad_material.GetStatus().Diagnostics()[0].message.find("finite"), + std::string::npos); + auto x = Brick(); + for (auto& a : x) + for (auto& value : a) value *= 1e200; + auto bad_geometry = Make(x); + ASSERT_FALSE(bad_geometry.HasValue()); + EXPECT_NE(bad_geometry.GetStatus().Diagnostics()[0].message.find("nonfinite"), + std::string::npos); + auto element = Make(Brick(), 120.0, .25, huge); + ASSERT_TRUE(element.HasValue()); + auto bad_mass = element.Value().ComputeMass(); + ASSERT_FALSE(bad_mass.HasValue()); + EXPECT_NE(bad_mass.GetStatus().Diagnostics()[0].message.find("nonfinite"), + std::string::npos); + // Current detJ is finite (~6e300), but the StVK energy and pushed-forward + // stress arithmetic overflow despite every input displacement being finite. + const auto u = + Affine(Brick(), + Tensor{{{1e100, 0.0, 0.0}, {0.0, 1e100, 0.0}, {0.0, 0.0, 1e100}}}); + auto bad_trial = element.Value().EvaluateNonlinear(u); + ASSERT_FALSE(bad_trial.HasValue()); + EXPECT_NE(bad_trial.GetStatus().Diagnostics()[0].message.find("trial"), + std::string::npos); + EXPECT_NE(bad_trial.GetStatus().Diagnostics()[0].message.find("nonfinite"), + std::string::npos); +} + +// Independent Gauss-Legendre abscissae from polynomial root iteration. This +// test does not call the kernel's shapes, Jacobian, inverse or quadrature. +std::vector> LegendreRule(unsigned count) { + std::vector> points; + for (unsigned root = 0; root < count; ++root) { + double x = std::cos(std::acos(-1.0) * (static_cast(root) + .75) / + (count + .5)); + double derivative = 0.0; + for (unsigned iteration = 0; iteration < 30; ++iteration) { + double previous = 1.0, current = x; + for (unsigned n = 2; n <= count; ++n) { + const double next = + ((2.0 * n - 1.0) * x * current - (n - 1.0) * previous) / n; + previous = current; + current = next; + } + derivative = count * (x * current - previous) / (x * x - 1.0); + const double correction = current / derivative; + x -= correction; + if (std::abs(correction) < 2e-16) break; + } + points.push_back({x, 2.0 / ((1.0 - x * x) * derivative * derivative)}); + } + return points; +} +Triple Cross(const Triple& a, const Triple& b) { + return {a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], + a[0] * b[1] - a[1] * b[0]}; +} +std::array WarpedQuadratureEnergies(const fesa::Vector& u, + unsigned count) { + const auto rule = LegendreRule(count); + std::array result{}; + for (const auto& xi : rule) + for (const auto& eta : rule) + for (const auto& zeta : rule) { + const Triple natural{xi[0], eta[0], zeta[0]}; + // Derivatives of the analytic warped map, rather than recovered nodal + // J. + const Triple c0{1.0, .09 * zeta[0], .15 * eta[0]}; + const Triple c1{.12 * zeta[0], 2.0, .15 * xi[0]}; + const Triple c2{.12 * eta[0], .09 * xi[0], 3.0}; + Tensor reciprocal{Cross(c1, c2), Cross(c2, c0), Cross(c0, c1)}; + const double det = c0[0] * reciprocal[0][0] + c0[1] * reciprocal[0][1] + + c0[2] * reciprocal[0][2]; + for (auto& row : reciprocal) + for (auto& value : row) value /= det; + Tensor natural_gradient{}; + Triple velocity{}; + for (std::size_t a = 0; a < 8; ++a) { + double shape = .125; + for (std::size_t j = 0; j < 3; ++j) + shape *= 1.0 + kSigns[a][j] * natural[j]; + for (std::size_t i = 0; i < 3; ++i) + velocity[i] += shape * u[3 * a + i]; + for (std::size_t j = 0; j < 3; ++j) { + double derivative = .125 * kSigns[a][j]; + for (std::size_t other = 0; other < 3; ++other) + if (other != j) + derivative *= 1.0 + kSigns[a][other] * natural[other]; + for (std::size_t i = 0; i < 3; ++i) + natural_gradient[i][j] += derivative * u[3 * a + i]; + } + } + Tensor gradient{}; + for (std::size_t i = 0; i < 3; ++i) + for (std::size_t j = 0; j < 3; ++j) + for (std::size_t d = 0; d < 3; ++d) + gradient[i][j] += natural_gradient[i][d] * reciprocal[d][j]; + double trace = 0.0, strain_squared = 0.0, velocity_squared = 0.0; + for (std::size_t i = 0; i < 3; ++i) { + trace += gradient[i][i]; + velocity_squared += velocity[i] * velocity[i]; + for (std::size_t j = 0; j < 3; ++j) { + const double strain = .5 * (gradient[i][j] + gradient[j][i]); + strain_squared += strain * strain; + } + } + const double weight = xi[1] * eta[1] * zeta[1] * det; + result[0] += weight * (24.0 * trace * trace + 48.0 * strain_squared); + result[1] += weight * velocity_squared; // .5*rho*v^2 with rho=2. + } + return result; +} +TEST(Hex8Element, WarpedNonaffineQuadratureDifferenceIsQuantifiedAndRefines) { + auto element = Make(Warped()); + ASSERT_TRUE(element.HasValue()); + fesa::Vector u{24U}; + for (std::size_t i = 0; i < 24; ++i) + u[i] = .03 * std::sin(.7 * static_cast(i)); + const auto k = element.Value().ComputeStiffness().Value().values; + const auto m = element.Value().ComputeMass().Value().values; + const auto q2 = WarpedQuadratureEnergies(u, 2), + q3 = WarpedQuadratureEnergies(u, 3); + const auto q5 = WarpedQuadratureEnergies(u, 5), + q7 = WarpedQuadratureEnergies(u, 7); + const double budget = 512.0 * std::numeric_limits::epsilon(); + EXPECT_NEAR(.5 * u.Dot(k.Multiply(u)), q2[0], budget * q2[0]); + EXPECT_NEAR(.5 * u.Dot(m.Multiply(u)), q2[1], budget * q2[1]); + const double stiffness_discrepancy = std::abs(q2[0] - q7[0]) / q7[0]; + const double mass_discrepancy = std::abs(q2[1] - q7[1]) / q7[1]; + // Mandated2x2 is retained: higher-order integration gives a distinct warped + // discretization. Rational stiffness integrands converge with refinement. + EXPECT_GT(stiffness_discrepancy, 1e-8); + EXPECT_GT(mass_discrepancy, 1e-8); + EXPECT_LT(std::abs(q5[0] - q7[0]), .05 * std::abs(q3[0] - q7[0])); + // N_a*N_b*detJ has per-axis degree<=4, so order3 integrates mass exactly. + EXPECT_NEAR(q3[1], q7[1], budget * q7[1]); + EXPECT_NEAR(q5[1], q7[1], budget * q7[1]); + std::cout << "warped K2-vs-K7 energy relative=" << stiffness_discrepancy + << " M2-vs-M7 kinetic relative=" << mass_discrepancy + << " K3-vs-K7=" << std::abs(q3[0] - q7[0]) / q7[0] + << " K5-vs-K7=" << std::abs(q5[0] - q7[0]) / q7[0] << '\n'; +} } // namespace diff --git a/tests/unit/results/result_recovery_test.cpp b/tests/unit/results/result_recovery_test.cpp index 67ea101..c80bf81 100644 --- a/tests/unit/results/result_recovery_test.cpp +++ b/tests/unit/results/result_recovery_test.cpp @@ -21,6 +21,7 @@ #include "fesa/assembly/parallel_for.h" #include "fesa/assembly/sparse_assembler.h" #include "fesa/elements/element.h" +#include "fesa/elements/hex8_element.h" #include "fesa/fem/dof_manager.h" #include "fesa/math/vector3.h" #include "fesa/model/domain.h" @@ -1149,3 +1150,50 @@ TEST(ResultRecovery, InvalidLaterShellLeavesEntirePriorStateUnchanged) { EXPECT_EQ(state.Equilibrium(), prior_equilibrium); EXPECT_EQ(state.VerificationMetrics(), prior_metrics); } + +TEST(ResultRecovery, RealSolidPayloadFailsClosedAndPreservesLegacyState) { + const auto fixture = MakeFixture(); + const auto input = MakeAxialEquilibriumState(fixture); + auto state = fesa::AnalysisState::Create(*fixture.dofs, {"prior-step", 9U}); + state.Displacement()[0] = 37.0; + state.InternalForce()[0] = 19.0; + const auto prior_bundle = MakeFakeBeamBundle(*fixture.domain, 0U, 11.0); + state.EndpointResults() = + std::get(prior_bundle.payload).endpoint_rows; + fesa::Hex8Element::Coordinates x{{{-1, -1, -1}, + {1, -1, -1}, + {1, 1, -1}, + {-1, 1, -1}, + {-1, -1, 1}, + {1, -1, 1}, + {1, 1, 1}, + {-1, 1, 1}}}; + const auto identity = fixture.domain->Elements()[0U].SourceId(); + fesa::ElementDofLayout layout{ + identity, + {0, 1, 2, 3, 4, 5, 6, 7}, + {fesa::DofComponent::kUx, fesa::DofComponent::kUy, + fesa::DofComponent::kUz}}; + const fesa::LinearElasticMaterial material{"solid", 120.0, .25, {}}; + auto solid = fesa::Hex8Element::Create(0U, layout, x, material); + ASSERT_TRUE(solid.HasValue()); + auto recovered = solid.Value().Recover(fesa::Vector{24U}); + ASSERT_TRUE(recovered.HasValue()); + ASSERT_EQ(std::get(recovered.Value().payload) + .rows.size(), + 8U); + // A legacy six-DOF beam scatter forwards the real solid payload to the + // generic visitor; it must reject the incompatible payload before commit. + FakeRecoveryElement runtime{MakeRuntimeLayout(*fixture.domain, 0U), + std::move(recovered.Value())}; + const fesa::Status status = fesa::ResultRecovery::Recover( + *fixture.model, {std::cref(runtime)}, *fixture.dofs, *fixture.stiffness, + input.Displacement(), input.ExternalForce(), state); + ExpectStatusCode(status, "unsupported-solid-result-recovery"); + EXPECT_DOUBLE_EQ(state.Displacement()[0], 37.0); + EXPECT_DOUBLE_EQ(state.InternalForce()[0], 19.0); + ASSERT_EQ(state.EndpointResults().size(), 2U); + EXPECT_DOUBLE_EQ(state.EndpointResults()[0].section_resultant[0], 11.0); + EXPECT_EQ(state.Identity().step_name, "prior-step"); + EXPECT_EQ(state.Identity().frame_index, 9U); +}