From df1bf94369e19a3bb1db936bfbca55046afe3e0a Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Tue, 20 Oct 2026 09:27:58 +0900 Subject: [PATCH] Preserve empty MITC4 hidden-constraint boundary --- src/fesa/analysis/hht_dynamic_algebra.cpp | 20 +++--- .../analysis/hht_dynamic_algebra_test.cpp | 62 ++++++++++++++++++- .../unit/analysis/hht_reduced_shell_test.cpp | 20 ++++++ 3 files changed, 90 insertions(+), 12 deletions(-) diff --git a/src/fesa/analysis/hht_dynamic_algebra.cpp b/src/fesa/analysis/hht_dynamic_algebra.cpp index 7b4c8cd..9311989 100644 --- a/src/fesa/analysis/hht_dynamic_algebra.cpp +++ b/src/fesa/analysis/hht_dynamic_algebra.cpp @@ -64,16 +64,6 @@ bool Within(const double value, const double scale, const double tolerance) { Status Validate(const Domain& domain, const DofManager& dofs, const DynamicAssembly& assembly, const double length) { const std::size_t count = dofs.FullDofCount(); - if (assembly.frame.mitc4_drilling && - assembly.frame.mitc4_drilling->mode == Mitc4DrillingMode::kDirectSpring) { - if (assembly.hidden_constraint.Size() != assembly.rotations.size()) - return Failure(domain, "dynamic-invalid-algebra-state", - "Canonical hidden values must cover every node."); - for (std::size_t i = 0U; i < assembly.hidden_constraint.Size(); ++i) - if (!std::isfinite(assembly.hidden_constraint[i])) - return Failure(domain, "dynamic-invalid-algebra-state", - "Canonical hidden values must remain finite."); - } if (!std::isfinite(length) || !(length > 0.0) || count % 6U != 0U || assembly.rotations.size() != count / 6U) { return Failure(domain, "dynamic-invalid-algebra-state", @@ -559,6 +549,16 @@ Status HiddenEvidence(const Domain& domain, const DynamicAssembly& assembly, const bool require_acceleration) { const bool direct = frame->mitc4_drilling && frame->mitc4_drilling->mode == Mitc4DrillingMode::kDirectSpring; + if (direct) { + for (const std::size_t node : nodes) { + if (node >= assembly.hidden_constraint.Size()) + return Failure(domain, "dynamic-invalid-algebra-state", + "Canonical hidden values must cover every admissible null node."); + if (!std::isfinite(assembly.hidden_constraint[node])) + return Failure(domain, "dynamic-invalid-algebra-state", + "Canonical hidden values must remain finite."); + } + } Vector psi{nodes.size()}; Vector psi_scale{nodes.size()}; Vector velocity_defect{nodes.size()}; diff --git a/tests/unit/analysis/hht_dynamic_algebra_test.cpp b/tests/unit/analysis/hht_dynamic_algebra_test.cpp index a6da082..06e1c87 100644 --- a/tests/unit/analysis/hht_dynamic_algebra_test.cpp +++ b/tests/unit/analysis/hht_dynamic_algebra_test.cpp @@ -28,7 +28,8 @@ struct Fixture { Fixture MakeFixture(const std::size_t free_nodes = 1U, const bool partial_rotation = false, const bool translations_free = true, - const double youngs_modulus = 1000.0) { + const double youngs_modulus = 1000.0, + const bool all_rotations_constrained = false) { fesa::ModelDefinition definition{}; definition.source_path = "hht-algebra.inp"; definition.nodes = { @@ -66,7 +67,10 @@ Fixture MakeFixture(const std::size_t free_nodes = 1U, step.boundaries.push_back( {name, 1, 6, 0.0, {definition.source_path, 10U}}); } else { - if (partial_rotation) + if (all_rotations_constrained) + step.boundaries.push_back( + {name, 4, 6, 0.0, {definition.source_path, 10U}}); + else if (partial_rotation) step.boundaries.push_back( {name, 4, 5, 0.0, {definition.source_path, 10U}}); if (!translations_free) @@ -705,3 +709,57 @@ TEST(HhtEndpointFields, MeasuresActualFullChartAcrossNonlinearRootAndPartialBoun } } // namespace + +TEST(HhtDynamicAlgebra, NonDirectionalAssemblyRetainsValidEmptyNullBoundary) { + for (const std::size_t free_nodes : {4U, 0U}) { + SCOPED_TRACE(free_nodes); + const auto fixture = MakeFixture(free_nodes, false, true, 1000.0, true); + const auto shells = fesa::hht_internal::BuildRuntimeShells(*fixture.domain, + fixture.dofs); + ASSERT_TRUE(shells.HasValue()); + auto trial = EmptyTrial(fixture); + trial.external_force[0U] = 0.4; + fesa::Vector load_scale(24U); + load_scale[0U] = 0.4; + fesa::SerialParallelFor serial; + const auto parameters = fesa::MakeHhtParameters(-0.05); + ASSERT_TRUE(parameters.HasValue()); + const auto assembly = fesa::hht_internal::EvaluateAssembly( + *fixture.domain, fixture.dofs, shells.Value(), trial, load_scale, + 0.1, parameters.Value(), serial, false); + ASSERT_TRUE(assembly.HasValue()); + ASSERT_TRUE(assembly.Value().frame.mitc4_drilling); + EXPECT_EQ(assembly.Value().frame.mitc4_drilling->mode, + fesa::Mitc4DrillingMode::kDirectSpring); + EXPECT_EQ(assembly.Value().hidden_constraint.Size(), 0U); + EXPECT_FALSE(assembly.Value().hidden_gradient.has_value()); + const auto nodes = fesa::hht_internal::NullNodes(fixture.dofs, + assembly.Value().rotations); + ASSERT_TRUE(nodes.empty()); + fesa::MklPardisoSolver positive; + fesa::MklPardisoTangentSolver tangent; + const auto acceleration = fesa::hht_internal::InitializeAcceleration( + *fixture.domain, fixture.dofs, assembly.Value(), 2.0, positive, tangent); + EXPECT_TRUE(acceleration.HasValue()) + << (acceleration.GetStatus().Diagnostics().empty() ? "" : + acceleration.GetStatus().Diagnostics()[0U].message); + auto frame = assembly.Value().frame; + if (acceleration.HasValue()) { + frame.acceleration = acceleration.Value(); + // rho*t*area/4 = .2 per node; the .4 force gives acceleration 2. + EXPECT_NEAR(acceleration.Value()[0U], free_nodes == 0U ? 0.0 : 2.0, 1e-12); + for (const auto full : fixture.dofs.ConstrainedDofs()) + EXPECT_DOUBLE_EQ(acceleration.Value()[full], 0.0); + } + const auto complete = fesa::hht_internal::CompleteConstraintEvidence( + *fixture.domain, fixture.dofs, assembly.Value(), 2.0, nodes, + positive, tangent, &frame); + EXPECT_TRUE(complete.IsOk()) + << (complete.Diagnostics().empty() ? "" : complete.Diagnostics()[0U].message); + if (complete.IsOk()) { + EXPECT_EQ(frame.hidden_constraint.Size(), 0U); + EXPECT_EQ(frame.hidden_velocity_defect.Size(), 0U); + EXPECT_EQ(frame.hidden_acceleration_defect.Size(), 0U); + } + } +} diff --git a/tests/unit/analysis/hht_reduced_shell_test.cpp b/tests/unit/analysis/hht_reduced_shell_test.cpp index 4e0d477..10577c7 100644 --- a/tests/unit/analysis/hht_reduced_shell_test.cpp +++ b/tests/unit/analysis/hht_reduced_shell_test.cpp @@ -1938,4 +1938,24 @@ TEST(Mitc4DrillingEvidence, DirectHiddenCertificateUsesCanonicalPrimitive) { EXPECT_DOUBLE_EQ(frame.hidden_velocity_defect[row], 0.0); EXPECT_DOUBLE_EQ(frame.hidden_acceleration_defect[row], 0.0); } + const auto accepted = frame; + for (const bool missing : {true, false}) { + auto invalid = assembly.Value(); + if (missing) + invalid.hidden_constraint = fesa::Vector(0U); + else + invalid.hidden_constraint[0U] = std::numeric_limits::infinity(); + const auto rejected = fesa::hht_internal::CompleteConstraintEvidence( + *fixture.domain, fixture.dofs, invalid, length.Value(), null_nodes, + positive, tangent, &frame); + EXPECT_FALSE(rejected.IsOk()); + for (const auto& pair : { + std::make_pair(&frame.hidden_constraint, &accepted.hidden_constraint), + std::make_pair(&frame.hidden_velocity_defect, &accepted.hidden_velocity_defect), + std::make_pair(&frame.hidden_acceleration_defect, &accepted.hidden_acceleration_defect)}) { + ASSERT_EQ(pair.first->Size(), pair.second->Size()); + for (std::size_t row = 0U; row < pair.first->Size(); ++row) + ExpectSameBits((*pair.first)[row], (*pair.second)[row]); + } + } }