From 775602860f88a77f9168980a58bece33ca96020f Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Wed, 12 Aug 2026 21:54:14 +0900 Subject: [PATCH] feat(linear-static-mitc4-shell): step 12 - shell-linear-static-flow --- src/fesa/results/result_recovery.cpp | 42 +++- .../analysis/linear_static_analysis_test.cpp | 232 ++++++++++++++++++ .../integration/app/fesa_application_test.cpp | 75 ++++++ 3 files changed, 347 insertions(+), 2 deletions(-) diff --git a/src/fesa/results/result_recovery.cpp b/src/fesa/results/result_recovery.cpp index 6623d8b..722cd4a 100644 --- a/src/fesa/results/result_recovery.cpp +++ b/src/fesa/results/result_recovery.cpp @@ -85,6 +85,39 @@ double indexedNorm(const Vector& values, return result; } +std::array freeEquationInternalTermNorms( + const SparseMatrix& stiffness, + const Vector& displacement, + const DofManager& dofs) { + std::vector freeColumn(dofs.fullDofCount(), 0U); + for (const std::size_t fullDof : dofs.freeDofs()) { + freeColumn[fullDof] = 1U; + } + + double freeTermNorm = 0.0; + double constrainedTermNorm = 0.0; + for (const std::size_t row : dofs.freeDofs()) { + double freeTerm = 0.0; + double constrainedTerm = 0.0; + for (std::size_t position = stiffness.rowOffsets()[row]; + position < stiffness.rowOffsets()[row + 1U]; + ++position) { + const std::size_t column = stiffness.columnIndices()[position]; + const double contribution = + stiffness.values()[position] * displacement[column]; + if (freeColumn[column] != 0U) { + freeTerm += contribution; + } else { + constrainedTerm += contribution; + } + } + freeTermNorm = std::hypot(freeTermNorm, freeTerm); + constrainedTermNorm = + std::hypot(constrainedTermNorm, constrainedTerm); + } + return {freeTermNorm, constrainedTermNorm}; +} + bool strictlyIncreasing(const std::vector& values) { return std::adjacent_find( values.begin(), values.end(), @@ -639,10 +672,15 @@ Status ResultRecovery::recover(const AnalysisModel& model, } const double residualNorm = indexedNorm(residual, dofs.freeDofs()); - const double internalNorm = indexedNorm(internalForce, dofs.freeDofs()); + const auto internalTermNorms = freeEquationInternalTermNorms( + fullStiffness, state.displacement(), dofs); const double externalNorm = indexedNorm(state.externalForce(), dofs.freeDofs()); - const double denominator = (std::max)(internalNorm, externalNorm); + // Normalize against the three terms of + // Kff*df + Kfc*dc - Ff. Using only the already-cancelled K*d term would + // classify prescribed-only equilibrium roundoff as a unit residual. + const double denominator = (std::max)({ + internalTermNorms[0U], internalTermNorms[1U], externalNorm}); if (!std::isfinite(residualNorm) || !std::isfinite(denominator)) { return recoveryFailure( "nonfinite-recovery-value", diff --git a/tests/integration/analysis/linear_static_analysis_test.cpp b/tests/integration/analysis/linear_static_analysis_test.cpp index 867c053..5f59f01 100644 --- a/tests/integration/analysis/linear_static_analysis_test.cpp +++ b/tests/integration/analysis/linear_static_analysis_test.cpp @@ -11,6 +11,7 @@ #include #include #include +#include #include #include #include @@ -95,6 +96,60 @@ Tip, 1, )inp" + std::to_string(tipForce) + R"inp( )inp"; } +std::string shellDeck( + const std::string& boundaryBlock, + const std::string& loadBlock = {}) { + return std::string{R"inp(*Part, name=ShellPart +*Node +1, 0., 0., 0. +2, 1., 0., 0. +3, 1., 1., 0. +4, 0., 1., 0. +*Element, type=S4 +1, 1, 2, 3, 4 +*Elset, elset=ShellSet +1 +*Shell Section, elset=ShellSet, material=Steel +0.1 +*End Part +*Assembly, name=Assembly +*Instance, name=Shell-1, part=ShellPart +*End Instance +*Nset, nset=N1, instance=Shell-1 +1 +*Nset, nset=N2, instance=Shell-1 +2 +*Nset, nset=N3, instance=Shell-1 +3 +*Nset, nset=N4, instance=Shell-1 +4 +*Nset, nset=All, instance=Shell-1 +1, 2, 3, 4 +*End Assembly +*Material, name=Steel +*Elastic +1000., 0.25 +)inp"} + boundaryBlock + R"inp(*Step, name=Load, nlgeom=NO +*Static +0.1, 1., 0.01, 1. +)inp" + loadBlock + R"inp(*End Step +)inp"; +} + +std::string allConstrainedShellDeck() { + return shellDeck("*Boundary\nAll, 1, 6\n"); +} + +std::string prescribedShellDeck() { + return shellDeck(R"inp(*Boundary +N1, 1, 6 +N2, 1, 1, 0.1 +N2, 2, 6 +N3, 2, 6 +N4, 1, 6 +)inp"); +} + // The pure Template Method spy makes the eight public lifecycle hooks observable // without coupling the ordering assertion to any solver backend. class SpyAnalysis final : public fesa::Analysis { @@ -167,6 +222,61 @@ private: mutable int solveCalls_{0}; }; +class RecordingMklSolver final : public fesa::LinearSolver { +public: + fesa::Status factorize(const fesa::SparseMatrix& matrix) override { + ++factorizeCalls_; + factorizedDimension_ = matrix.rows(); + if (matrix.rows() == 1U && matrix.columns() == 1U && + matrix.values().size() == 1U) { + scalarStiffness_ = matrix.values()[0U]; + } + return backend_.factorize(matrix); + } + + fesa::Status solve( + const fesa::Vector& rhs, fesa::Vector& solution) const override { + ++solveCalls_; + if (rhs.size() == 0U) { + rhs_.clear(); + } else { + rhs_.assign(rhs.data(), rhs.data() + rhs.size()); + } + return backend_.solve(rhs, solution); + } + + int factorizeCalls() const noexcept { return factorizeCalls_; } + int solveCalls() const noexcept { return solveCalls_; } + std::size_t factorizedDimension() const noexcept { + return factorizedDimension_; + } + double scalarStiffness() const noexcept { return scalarStiffness_; } + const std::vector& rhs() const noexcept { return rhs_; } + +private: + fesa::MklPardisoSolver backend_; + int factorizeCalls_{0}; + mutable int solveCalls_{0}; + std::size_t factorizedDimension_{0U}; + double scalarStiffness_{0.0}; + mutable std::vector rhs_; +}; + +class NonfiniteLinearSolver final : public fesa::LinearSolver { +public: + fesa::Status factorize(const fesa::SparseMatrix&) override { + return fesa::Status::ok(); + } + + fesa::Status solve( + const fesa::Vector&, fesa::Vector& solution) const override { + for (std::size_t index = 0U; index < solution.size(); ++index) { + solution[index] = (std::numeric_limits::quiet_NaN)(); + } + return fesa::Status::ok(); + } +}; + class SpyResultsWriter final : public fesa::ResultsWriter { public: explicit SpyResultsWriter(std::vector& events) @@ -198,6 +308,7 @@ public: const std::vector& diagnostics) override { outputPath_ = outputPath; nodeCount_ = domain.nodes().size(); + shellElementCount_ = domain.shellElements().size(); state_ = std::make_unique(state); diagnostics_ = diagnostics; return fesa::Status::ok(); @@ -214,6 +325,9 @@ public: return outputPath_; } std::size_t nodeCount() const noexcept { return nodeCount_; } + std::size_t shellElementCount() const noexcept { + return shellElementCount_; + } const std::vector& diagnostics() const noexcept { return diagnostics_; } @@ -221,6 +335,7 @@ public: private: std::filesystem::path outputPath_; std::size_t nodeCount_{0U}; + std::size_t shellElementCount_{0U}; std::unique_ptr state_; std::vector diagnostics_; }; @@ -298,3 +413,120 @@ TEST(LinearStaticCli, RealPipelineHandlesAnalyticalAndNonzeroPrescription) { EXPECT_EQ(state.gaussResults().size(), 2U); EXPECT_EQ(state.stressResults().size(), 2U); } + +// MITC4-FLOW-001 +TEST(Mitc4ShellCli, UsesExistingLifecycleAndExactlyOneFactorization) { + TempDirectory directory{"shell-order"}; + const auto input = directory.path() / "all-constrained-shell.inp"; + const auto output = directory.path() / "results.h5"; + writeText(input, allConstrainedShellDeck()); + + std::vector adapterEvents; + fesa::SerialParallelFor serial; + SpyLinearSolver solver{adapterEvents}; + SpyResultsWriter writer{adapterEvents}; + fesa::LinearStaticAnalysis analysis{serial, solver, writer}; + + ASSERT_TRUE(analysis.run({input, output}).isOk()); + EXPECT_EQ(solver.factorizeCalls(), 1); + EXPECT_EQ(solver.solveCalls(), 1); + EXPECT_EQ(writer.writeCalls(), 1); + EXPECT_EQ( + adapterEvents, + (std::vector{ + "solver-factorize", "solver-solve", "writer-write"})); +} + +// MITC4-FLOW-002 +TEST(Mitc4ShellCli, AppliesKfcForNonzeroPrescribedDisplacement) { + TempDirectory directory{"shell-prescribed"}; + const auto input = directory.path() / "prescribed-shell.inp"; + const auto output = directory.path() / "captured-results.h5"; + writeText(input, prescribedShellDeck()); + + fesa::SerialParallelFor serial; + RecordingMklSolver solver; + CapturingResultsWriter writer; + fesa::LinearStaticAnalysis analysis{serial, solver, writer}; + + const auto status = analysis.run({input, output}); + for (const auto& diagnostic : status.diagnostics()) { + EXPECT_TRUE(status.isOk()) + << diagnostic.code << ": " << diagnostic.message; + } + ASSERT_TRUE(status.isOk()); + ASSERT_EQ(solver.factorizeCalls(), 1); + ASSERT_EQ(solver.solveCalls(), 1); + ASSERT_EQ(solver.factorizedDimension(), 1U); + ASSERT_EQ(solver.rhs().size(), 1U); + EXPECT_NEAR(solver.scalarStiffness(), 440.0 / 9.0, 1.0e-12); + // Ff is exactly zero, so this nonzero RHS is solely -Kfc*dc. + EXPECT_NEAR(solver.rhs()[0U], -4.0 / 9.0, 1.0e-12); + + EXPECT_EQ(writer.outputPath(), output); + EXPECT_EQ(writer.nodeCount(), 4U); + EXPECT_EQ(writer.shellElementCount(), 1U); + const auto& state = writer.state(); + ASSERT_EQ(state.displacement().size(), 24U); + EXPECT_NEAR(state.displacement()[6U], 0.1, 1.0e-12); + EXPECT_NEAR(state.displacement()[12U], -1.0 / 110.0, 1.0e-12); + EXPECT_NEAR(state.verificationMetrics()[0U], 0.0, 1.0e-10); + EXPECT_EQ(state.shellResults().size(), 4U); + EXPECT_GT(state.physicalStrainEnergy(), 0.0); +} + +// MITC4-FLOW-003 +TEST(Mitc4ShellCli, RejectsSingularAndAcceptsZeroByZeroFreeSystem) { + TempDirectory directory{"shell-singular-all"}; + const auto singularInput = directory.path() / "singular-shell.inp"; + const auto constrainedInput = directory.path() / "constrained-shell.inp"; + writeText(singularInput, shellDeck("")); + writeText(constrainedInput, allConstrainedShellDeck()); + + fesa::SerialParallelFor serial; + fesa::MklPardisoSolver singularSolver; + std::vector singularEvents; + SpyResultsWriter singularWriter{singularEvents}; + fesa::LinearStaticAnalysis singularAnalysis{ + serial, singularSolver, singularWriter}; + const auto singular = singularAnalysis.run( + {singularInput, directory.path() / "singular.h5"}); + ASSERT_FALSE(singular.isOk()); + EXPECT_EQ(singular.failureCategory(), fesa::FailureCategory::solver); + EXPECT_EQ(singularWriter.writeCalls(), 0); + + RecordingMklSolver constrainedSolver; + CapturingResultsWriter constrainedWriter; + fesa::LinearStaticAnalysis constrainedAnalysis{ + serial, constrainedSolver, constrainedWriter}; + const auto constrained = constrainedAnalysis.run( + {constrainedInput, directory.path() / "constrained.h5"}); + ASSERT_TRUE(constrained.isOk()); + EXPECT_EQ(constrainedSolver.factorizeCalls(), 1); + EXPECT_EQ(constrainedSolver.factorizedDimension(), 0U); + EXPECT_EQ(constrainedSolver.solveCalls(), 1); + EXPECT_TRUE(constrainedSolver.rhs().empty()); + EXPECT_EQ(constrainedWriter.state().shellResults().size(), 4U); +} + +// MITC4-FLOW-004 +TEST(Mitc4ShellCli, DoesNotWriteAnInvalidRecoveryCandidate) { + TempDirectory directory{"shell-invalid-candidate"}; + const auto input = directory.path() / "invalid-recovery-shell.inp"; + writeText(input, prescribedShellDeck()); + + fesa::SerialParallelFor serial; + NonfiniteLinearSolver solver; + std::vector adapterEvents; + SpyResultsWriter writer{adapterEvents}; + fesa::LinearStaticAnalysis analysis{serial, solver, writer}; + const auto status = analysis.run( + {input, directory.path() / "must-not-exist.h5"}); + + ASSERT_FALSE(status.isOk()); + EXPECT_EQ(status.failureCategory(), fesa::FailureCategory::model); + ASSERT_FALSE(status.diagnostics().empty()); + EXPECT_EQ(status.diagnostics().front().code, "nonfinite-recovery-value"); + EXPECT_EQ(writer.writeCalls(), 0); + EXPECT_TRUE(adapterEvents.empty()); +} diff --git a/tests/integration/app/fesa_application_test.cpp b/tests/integration/app/fesa_application_test.cpp index 2642e87..04ad264 100644 --- a/tests/integration/app/fesa_application_test.cpp +++ b/tests/integration/app/fesa_application_test.cpp @@ -154,6 +154,38 @@ Tip, 1, 10. )inp"; } +std::string allConstrainedShellDeck() { + return R"inp(*Part, name=ShellPart +*Node +1, 0., 0., 0. +2, 1., 0., 0. +3, 1., 1., 0. +4, 0., 1., 0. +*Element, type=S4 +1, 1, 2, 3, 4 +*Elset, elset=ShellSet +1 +*Shell Section, elset=ShellSet, material=Steel +0.1 +*End Part +*Assembly, name=Assembly +*Instance, name=Shell-1, part=ShellPart +*End Instance +*Nset, nset=All, instance=Shell-1 +1, 2, 3, 4 +*End Assembly +*Material, name=Steel +*Elastic +1000., 0.25 +*Boundary +All, 1, 6 +*Step, name=Load, nlgeom=NO +*Static +0.1, 1., 0.01, 1. +*End Step +)inp"; +} + Hdf5Handle openFile(const std::filesystem::path& path) { const hid_t file = H5Fopen( path.string().c_str(), H5F_ACC_RDONLY, H5P_DEFAULT); @@ -259,6 +291,37 @@ void expectFesaHdf5Identity( 0U); } +void expectShellHdf5Identity( + const std::filesystem::path& output, + const std::filesystem::path& input) { + ASSERT_TRUE(std::filesystem::exists(output)); + ASSERT_GT(H5Fis_hdf5(output.string().c_str()), 0); + const auto file = openFile(output); + Hdf5Handle metadata{ + H5Gopen2(file.get(), "/metadata", H5P_DEFAULT), H5Gclose}; + ASSERT_GE(metadata.get(), 0); + EXPECT_EQ( + readStringAttribute(metadata.get(), "feature_id"), + "linear-static-mitc4-shell"); + EXPECT_GT( + H5Lexists( + file.get(), + "/steps/Step-1/frames/0/element/shell/generalized_strain", + H5P_DEFAULT), + 0); + EXPECT_EQ( + datasetDimensions( + file.get(), + "/steps/Step-1/frames/0/nodal/displacement"), + (std::vector{4U, 6U})); + const std::string normalizedInput = + std::filesystem::absolute(input).lexically_normal().generic_u8string(); + EXPECT_EQ( + readStringAttribute(metadata.get(), "source_input_identity").find( + "path=" + normalizedInput + ";content_identity="), + 0U); +} + struct AppRun { int exitCode; std::string standardError; @@ -407,3 +470,15 @@ TEST(LinearStaticCli, OutputRequestsDoNotFilterMandatoryResults) { ASSERT_EQ(requestedDiagnostics.size(), 1U); EXPECT_GT(requestedDiagnostics[0U], 0U); } + +// MITC4-FLOW-001: shell input uses the unchanged application route and syntax. +TEST(Mitc4ShellCli, WritesShellHdf5ThroughTheExistingApplicationRoute) { + TempDirectory directory{"shell-route"}; + const auto input = directory.path() / "shell.inp"; + const auto output = directory.path() / "shell-results.h5"; + writeText(input, allConstrainedShellDeck()); + + const auto result = runApplication(explicitOutputArguments(input, output)); + ASSERT_EQ(result.exitCode, 0) << result.standardError; + expectShellHdf5Identity(output, input); +}