diff --git a/PLAN.md b/PLAN.md index edef15c..e16d467 100644 --- a/PLAN.md +++ b/PLAN.md @@ -13,7 +13,7 @@ Every new agent session must read this file together with `PROGRESS.md` before p - If an item becomes obsolete, move it to `PROGRESS.md` with a short reason instead of silently deleting it. ## Current Objective -Continue the new Phase 1 rebaseline plan in `phases/1-linear-static-mitc4-rebaseline`, starting with P1R-12 assembly, solver-adapter boundary, constrained solve, and full-vector RF recovery revalidation. The old `phases/1-linear-static-mitc4` path is historical and superseded by the paper-based MITC4 formulation reset. +Continue the new Phase 1 rebaseline plan in `phases/1-linear-static-mitc4-rebaseline`, starting with P1R-13 linear static workflow revalidation from input to `U`/`RF` result fields. The old `phases/1-linear-static-mitc4` path is historical and superseded by the paper-based MITC4 formulation reset. ## Required Reading For New Agents 1. `AGENTS.md` @@ -36,7 +36,7 @@ Continue the new Phase 1 rebaseline plan in `phases/1-linear-static-mitc4-rebase ## Active Phase Files - Active phase directory: `phases/1-linear-static-mitc4-rebaseline` - Execute with: `python scripts/execute.py 1-linear-static-mitc4-rebaseline` -- Step numbering is zero-based. `step0.md` is complete and recorded in `phases/1-linear-static-mitc4-rebaseline/step0-audit.md`; `step1.md` is complete and created the `quad_02_phase1.inp` normalized reference path; `step2.md` is complete and revalidated core harness guardrails; `step3.md` is complete and revalidated the Phase 1 parser/domain subset; `step4.md` is complete and strengthened validation/singular diagnostics; `step5.md` is complete and revalidated the DofManager/reaction foundation; `step6.md` is complete and revalidated the minimum result model plus displacement CSV comparator; `step7.md` is complete and revalidated MITC4 natural coordinates, tying points, center directors, and integration bases; `step8.md` is complete and revalidated degenerated-continuum displacement, direct covariant strain rows, and MITC shear tying rows; `step9.md` is complete and revalidated plane-stress material, convected-to-local transform, and `2 x 2 x 2` material integration scaffolding; `step10.md` is complete and revalidated MITC4 stiffness, internal force, six-DOF transform, and drilling stabilization; `step11.md` is complete and added MITC4 membrane, bending, shear, twist, drilling-sensitivity, and thin-cantilever locking-sensitivity tests; `step15.md` is the independent evaluator closeout. +- Step numbering is zero-based. `step0.md` is complete and recorded in `phases/1-linear-static-mitc4-rebaseline/step0-audit.md`; `step1.md` is complete and created the `quad_02_phase1.inp` normalized reference path; `step2.md` is complete and revalidated core harness guardrails; `step3.md` is complete and revalidated the Phase 1 parser/domain subset; `step4.md` is complete and strengthened validation/singular diagnostics; `step5.md` is complete and revalidated the DofManager/reaction foundation; `step6.md` is complete and revalidated the minimum result model plus displacement CSV comparator; `step7.md` is complete and revalidated MITC4 natural coordinates, tying points, center directors, and integration bases; `step8.md` is complete and revalidated degenerated-continuum displacement, direct covariant strain rows, and MITC shear tying rows; `step9.md` is complete and revalidated plane-stress material, convected-to-local transform, and `2 x 2 x 2` material integration scaffolding; `step10.md` is complete and revalidated MITC4 stiffness, internal force, six-DOF transform, and drilling stabilization; `step11.md` is complete and added MITC4 membrane, bending, shear, twist, drilling-sensitivity, and thin-cantilever locking-sensitivity tests; `step12.md` is complete and revalidated full-space assembly, reduced projection, deterministic sparse-pattern scaffold, solver adapter injection, and full-vector internal/reaction force state; `step15.md` is the independent evaluator closeout. - Every step file contains a sprint contract with objective, required reading, scope, allowed files, explicit non-goals, tests to write first, reference artifacts, acceptance command, evaluator checklist, and handoff requirements. - Historical phase directory: `phases/1-linear-static-mitc4` - Historical phase status: blocked/superseded. Do not resume the old P1-15/P1-16 path unless the user explicitly requests recovery of that exact phase. @@ -75,7 +75,7 @@ Each gate should be satisfied before moving to the next implementation band unle | G0 - Planning readiness | partial | Readiness task R-011 is resolved by `quad_02_phase1.inp`; R-010 and R-013 remain open. | Updated docs, PLAN.md, PROGRESS.md | | G1 - Build and validation | satisfied | Build system, test framework, and `scripts/validate_workspace.py` run real checks. | Validation command output | | G2 - Parser and domain | satisfied | Parser subset revalidated in step 3; validation and singular diagnostics revalidated in step 4. | Parser acceptance/rejection tests, validation negative tests, and validation output | -| G3 - DOF/math/results infrastructure | partial | Core aliases, DOF mapping, validation harness, model diagnostic context, DofManager, sparse-connectivity inputs, full-vector reaction formula, result model metadata, and displacement CSV comparator were revalidated in steps 2, 5, and 6; assembly remains for step 12. | P1R-02, P1R-05, and P1R-06 validation output | +| G3 - DOF/math/results infrastructure | satisfied | Core aliases, DOF mapping, validation harness, model diagnostic context, DofManager, sparse-connectivity inputs, full-vector reaction formula, result model metadata, displacement CSV comparator, full-space assembly, reduced projection, sparse-pattern scaffold, and solver adapter boundary were revalidated in steps 2, 5, 6, and 12. | P1R-02, P1R-05, P1R-06, and P1R-12 validation output | | G4 - MITC4 element readiness | satisfied | MITC4 formulation was rewritten from local papers; Steps 7 through 11 rebuilt geometry/director/local-basis scaffolding, displacement interpolation, direct covariant strain rows, MITC shear tying rows, plane-stress material, convected-to-local transform, `2 x 2 x 2` material integration scaffolding, stiffness/internal force, six-DOF transform, drilling stabilization, and patch/locking-sensitivity tests. | P1R-07 through P1R-11 validation output | | G5 - End-to-end solver | reopened | Linear static path must be revalidated through steps 13 and 14 after the MITC4 rebuild and `quad_02` compatibility path. | Future integration/reference regression output | @@ -93,7 +93,7 @@ All milestones are intended to become one or more self-contained sprint contract | P1R-09 | completed | MITC4 generator | Implement material matrix, transform, and `2 x 2 x 2` integration scaffolding. | P1R-08 | Material/integration tests | | P1R-10 | completed | MITC4 generator | Assemble MITC4 stiffness/internal force with six-DOF transform and drilling stabilization. | P1R-09, P1R-05 | Symmetry, rigid body, drilling sensitivity tests | | P1R-11 | completed | verification generator | Add MITC4 patch, locking-sensitivity, and benchmark tests. | P1R-10 | Membrane/bending/shear/twist/locking tests | -| P1R-12 | pending | assembly generator | Rebuild assembly, solver adapter boundary, constrained solve, and full-vector RF recovery. | P1R-05, P1R-10 | Assembly and full-vector reaction tests | +| P1R-12 | completed | assembly generator | Rebuild assembly, solver adapter boundary, constrained solve, and full-vector RF recovery. | P1R-05, P1R-10 | Assembly and full-vector reaction tests | | P1R-13 | pending | analysis generator | Rebuild linear static workflow from input to U/RF result fields. | P1R-03, P1R-04, P1R-06, P1R-12 | End-to-end linear static tests | | P1R-14 | pending | reference generator | Run stored reference displacement regression using accepted Phase 1-compatible S4 cases. | P1R-13 | At least one automated CSV displacement regression | | P1R-15 | pending | evaluator | Independent Phase 1 evaluator closeout. | P1R-14 | Pass/fail report, synchronized PLAN/PROGRESS | diff --git a/PROGRESS.md b/PROGRESS.md index ddbac24..3f9b243 100644 --- a/PROGRESS.md +++ b/PROGRESS.md @@ -13,10 +13,38 @@ Every new agent session must read this file together with `PLAN.md` before plann - Do not remove history unless the user explicitly asks for archival cleanup. ## Current Status -Phase 1 has a new rebaseline phase definition in `phases/1-linear-static-mitc4-rebaseline`. Steps 0 through 11 are complete. `quad_02_phase1.inp` is now the normalized Phase 1-compatible input path for the stored `quad_02` S4 reference pair, while the original `quad_02.inp` remains preserved unsupported provenance. Core numeric aliases, DOF mapping, validation harness, model diagnostic context, the Phase 1 parser/domain subset, validation/singular diagnostics, DofManager/reaction foundation, minimum result model metadata, displacement CSV comparator foundation, MITC4 geometry/director scaffolding, MITC4 displacement/strain/tying row scaffolding, MITC4 material/transform/integration scaffolding, MITC4 stiffness/drilling/internal-force scaffolding, and MITC4 patch/locking-sensitivity tests have been revalidated. The old `phases/1-linear-static-mitc4` path is historical and superseded after the MITC4 formulation reset. +Phase 1 has a new rebaseline phase definition in `phases/1-linear-static-mitc4-rebaseline`. Steps 0 through 12 are complete. `quad_02_phase1.inp` is now the normalized Phase 1-compatible input path for the stored `quad_02` S4 reference pair, while the original `quad_02.inp` remains preserved unsupported provenance. Core numeric aliases, DOF mapping, validation harness, model diagnostic context, the Phase 1 parser/domain subset, validation/singular diagnostics, DofManager/reaction foundation, minimum result model metadata, displacement CSV comparator foundation, MITC4 geometry/director scaffolding, MITC4 displacement/strain/tying row scaffolding, MITC4 material/transform/integration scaffolding, MITC4 stiffness/drilling/internal-force scaffolding, MITC4 patch/locking-sensitivity tests, full-space assembly, reduced projection, sparse-pattern scaffold, solver adapter injection, and full-vector internal/reaction force state have been revalidated. The old `phases/1-linear-static-mitc4` path is historical and superseded after the MITC4 formulation reset. ## Completed Work +### 2026-05-04 - P1R-12 assembly sparse solver path completed +Author: Codex + +Changed files: +- `include/fesa/fesa.hpp` +- `tests/test_main.cpp` +- `phases/1-linear-static-mitc4-rebaseline/index.json` +- `PLAN.md` +- `PROGRESS.md` + +Summary: +- Added deterministic reduced sparse-pattern scaffolding from `DofManager` element equation connectivity, preserving int64 equation and nonzero counts for the future sparse/MKL path. +- Added `projectToReducedSystem()` so constrained/free projection is a named assembly boundary instead of being embedded directly inside `LinearStaticAnalysis`. +- Extended `AssemblyResult` to carry full-space stiffness, full external load, reduced sparse pattern, and diagnostics. +- Reworked assembly to call the rebuilt MITC4 stiffness result path directly and preserve full-space `K_full`/`F_full` for reaction recovery. +- Added solver injection to `LinearStaticAnalysis` so the deterministic Gaussian solver remains the default test adapter while future MKL-backed solvers stay behind the `LinearSolver` interface. +- Added `AnalysisState::f_internal_full` and verified `RF = Fint_full - Fext_full = K_full * U_full - F_full`. +- Added tests for deterministic sparse pattern ordering, reduced projection with known displacement, full-space assembly/load preservation, reduced residual satisfaction, solver adapter injection, and singular solver diagnostic propagation. + +Verification: +- First ran `python scripts/validate_workspace.py` after adding Step 12 tests; it failed as expected because reduced sparse pattern, projection, solver injection, and `f_internal_full` did not exist yet. +- After implementation, `python scripts/validate_workspace.py` configured CMake, built `fesa_core` and `fesa_tests`, and ran CTest successfully. +- CTest result: 1 test executable passed. + +Follow-up: +- Continue with P1R-13 linear static workflow revalidation from Phase 1 input through `U` and `RF` result fields. +- Step 12 intentionally did not add MKL/TBB APIs, stored Abaqus reference comparison, pressure loads, or a production sparse matrix storage backend. + ### 2026-05-04 - P1R-11 MITC4 patch and benchmark tests completed Author: Codex diff --git a/include/fesa/fesa.hpp b/include/fesa/fesa.hpp index 58df44c..4ab6150 100644 --- a/include/fesa/fesa.hpp +++ b/include/fesa/fesa.hpp @@ -1083,6 +1083,52 @@ class DofManager { std::vector constrained_full_indices_; }; +struct SparsePatternEntry { + EquationId row = 0; + EquationId col = 0; +}; + +struct SparsePattern { + EquationId equation_count = 0; + std::vector entries; + + SparseIndex nonzeroCount() const { + return static_cast(entries.size()); + } + + bool contains(EquationId row, EquationId col) const { + return std::any_of(entries.begin(), entries.end(), [&](const SparsePatternEntry& entry) { + return entry.row == row && entry.col == col; + }); + } +}; + +inline SparsePattern buildReducedSparsePattern(const Domain& domain, const DofManager& dofs) { + SparsePattern pattern; + pattern.equation_count = dofs.freeDofCount(); + std::set> ordered_entries; + for (const auto& [element_id, element] : domain.elements) { + (void)element_id; + const auto equations = dofs.elementEquationIds(element); + for (EquationId row : equations) { + if (row < 0) { + continue; + } + for (EquationId col : equations) { + if (col < 0) { + continue; + } + ordered_entries.insert({row, col}); + } + } + } + pattern.entries.reserve(ordered_entries.size()); + for (const auto& entry : ordered_entries) { + pattern.entries.push_back({entry.first, entry.second}); + } + return pattern; +} + class DenseMatrix { public: DenseMatrix() = default; @@ -2102,12 +2148,34 @@ class MITC4ElementKernel { struct AssemblyResult { DenseMatrix k_full; std::vector f_full; + SparsePattern reduced_pattern; std::vector diagnostics; + + bool ok() const { + return !hasError(diagnostics); + } +}; + +struct ReducedSystem { + DenseMatrix k; + std::vector f; + std::vector free_full_indices; + std::vector diagnostics; + + bool ok() const { + return !hasError(diagnostics); + } }; inline AssemblyResult assembleSystem(const Domain& domain, const DofManager& dofs, ElementStiffnessOptions options = {}) { - AssemblyResult result{DenseMatrix(dofs.fullDofCount(), dofs.fullDofCount()), std::vector(static_cast(dofs.fullDofCount()), 0.0), {}}; - MITC4ElementKernel kernel; + AssemblyResult result; + result.k_full = DenseMatrix(dofs.fullDofCount(), dofs.fullDofCount()); + result.f_full = std::vector(static_cast(dofs.fullDofCount()), 0.0); + result.reduced_pattern = buildReducedSparsePattern(domain, dofs); + if (dofs.freeDofCount() > 0 && result.reduced_pattern.nonzeroCount() == 0) { + result.diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-SINGULAR-SPARSE-PATTERN", + "Reduced sparse pattern has no stiffness entries", "assembly")); + } for (const auto& [element_id, element] : domain.elements) { const ShellSection* section = shellSectionForElement(domain, element_id); if (section == nullptr) { @@ -2123,19 +2191,63 @@ inline AssemblyResult assembleSystem(const Domain& domain, const DofManager& dof for (std::size_t i = 0; i < 4; ++i) { coordinates[i] = domain.nodes.at(element.node_ids[i]).coordinates; } - DenseMatrix ke = kernel.stiffness(coordinates, material_it->second.elastic_modulus, material_it->second.poisson_ratio, section->thickness, options); + const auto stiffness = mitc4ElementStiffness(coordinates, material_it->second.elastic_modulus, + material_it->second.poisson_ratio, section->thickness, options); + result.diagnostics.insert(result.diagnostics.end(), stiffness.diagnostics.begin(), stiffness.diagnostics.end()); + if (!stiffness.ok()) { + continue; + } const auto element_full_indices = dofs.elementFullDofIndices(element); for (LocalIndex a = 0; a < 24; ++a) { const LocalIndex ia = element_full_indices[static_cast(a)]; for (LocalIndex b = 0; b < 24; ++b) { const LocalIndex ib = element_full_indices[static_cast(b)]; - result.k_full.add(ia, ib, ke(a, b)); + result.k_full.add(ia, ib, stiffness.global(a, b)); } } } for (const NodalLoad& load : domain.loads) { + const auto dof = dofFromAbaqus(load.dof); + if (!dof) { + result.diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-ASSEMBLY-LOAD-DOF", + "Nodal load references an invalid DOF", "cload")); + continue; + } for (GlobalId node_id : resolveNodeTarget(domain, load.target, &result.diagnostics)) { - result.f_full[static_cast(dofs.fullIndex(node_id, *dofFromAbaqus(load.dof)))] += load.magnitude; + result.f_full[static_cast(dofs.fullIndex(node_id, *dof))] += load.magnitude; + } + } + return result; +} + +inline ReducedSystem projectToReducedSystem(const AssemblyResult& assembly, const DofManager& dofs) { + ReducedSystem result; + result.k = DenseMatrix(dofs.freeDofCount(), dofs.freeDofCount()); + result.f = std::vector(static_cast(dofs.freeDofCount()), 0.0); + result.free_full_indices = dofs.freeFullIndices(); + if (assembly.k_full.rows() != assembly.k_full.cols() || + assembly.k_full.rows() != dofs.fullDofCount() || + static_cast(assembly.f_full.size()) != dofs.fullDofCount()) { + result.diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-ASSEMBLY-SIZE", + "Full-space stiffness/load sizes do not match DofManager", "assembly")); + return result; + } + if (dofs.freeDofCount() == 0) { + result.diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-SINGULAR-NO-FREE-DOFS", + "No free DOFs exist after applying constraints", "dof")); + return result; + } + if (assembly.reduced_pattern.equation_count != dofs.freeDofCount() || assembly.reduced_pattern.nonzeroCount() == 0) { + result.diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-SINGULAR-SPARSE-PATTERN", + "Reduced sparse pattern is empty or inconsistent with free DOFs", "assembly")); + return result; + } + for (LocalIndex i = 0; i < dofs.freeDofCount(); ++i) { + const LocalIndex full_i = dofs.freeFullIndices()[static_cast(i)]; + result.f[static_cast(i)] = assembly.f_full[static_cast(full_i)]; + for (LocalIndex j = 0; j < dofs.freeDofCount(); ++j) { + const LocalIndex full_j = dofs.freeFullIndices()[static_cast(j)]; + result.k(i, j) = assembly.k_full(full_i, full_j); } } return result; @@ -2237,6 +2349,7 @@ class InMemoryResultsWriter { struct AnalysisState { std::vector u_full; std::vector f_external_full; + std::vector f_internal_full; std::vector reaction_full; bool converged = false; }; @@ -2273,6 +2386,9 @@ class Analysis { }; class LinearStaticAnalysis final : public Analysis { + public: + explicit LinearStaticAnalysis(const LinearSolver* solver = nullptr) : solver_(solver) {} + protected: void solve(const Domain& domain, AnalysisResult& result) const override { DofManager dofs(domain); @@ -2286,30 +2402,39 @@ class LinearStaticAnalysis final : public Analysis { if (hasError(result.diagnostics)) { return; } - DenseMatrix k_reduced(dofs.freeDofCount(), dofs.freeDofCount()); - std::vector f_reduced(static_cast(dofs.freeDofCount()), 0.0); - for (LocalIndex i = 0; i < dofs.freeDofCount(); ++i) { - const LocalIndex full_i = dofs.freeFullIndices()[static_cast(i)]; - f_reduced[static_cast(i)] = assembly.f_full[static_cast(full_i)]; - for (LocalIndex j = 0; j < dofs.freeDofCount(); ++j) { - const LocalIndex full_j = dofs.freeFullIndices()[static_cast(j)]; - k_reduced(i, j) = assembly.k_full(full_i, full_j); - } + const auto reduced = projectToReducedSystem(assembly, dofs); + result.diagnostics.insert(result.diagnostics.end(), reduced.diagnostics.begin(), reduced.diagnostics.end()); + if (hasError(result.diagnostics)) { + return; } - GaussianEliminationSolver solver; - SolveResult solved = solver.solve(k_reduced, f_reduced); + const LinearSolver& active_solver = solver_ == nullptr ? defaultSolver() : *solver_; + SolveResult solved = active_solver.solve(reduced.k, reduced.f); result.diagnostics.insert(result.diagnostics.end(), solved.diagnostics.begin(), solved.diagnostics.end()); if (!solved.ok()) { return; } + if (static_cast(solved.x.size()) != dofs.freeDofCount()) { + result.diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-SOLVER-SIZE", + "Linear solver returned a vector with the wrong size", "solver")); + return; + } result.state.u_full = dofs.reconstructFullVector(solved.x); result.state.f_external_full = assembly.f_full; + result.state.f_internal_full = assembly.k_full.multiply(result.state.u_full); result.state.reaction_full = recoverFullReaction(assembly.k_full, result.state.u_full, result.state.f_external_full); result.state.converged = true; InMemoryResultsWriter writer; writer.writeLinearStatic(domain, dofs, result.state.u_full, result.state.reaction_full); result.result_file = writer.result(); } + + private: + static const LinearSolver& defaultSolver() { + static const GaussianEliminationSolver solver; + return solver; + } + + const LinearSolver* solver_ = nullptr; }; struct CsvDisplacementRow { diff --git a/phases/1-linear-static-mitc4-rebaseline/index.json b/phases/1-linear-static-mitc4-rebaseline/index.json index 396fc85..38f2637 100644 --- a/phases/1-linear-static-mitc4-rebaseline/index.json +++ b/phases/1-linear-static-mitc4-rebaseline/index.json @@ -14,7 +14,7 @@ { "step": 9, "name": "mitc4-material-integration", "status": "completed" }, { "step": 10, "name": "mitc4-stiffness-drilling", "status": "completed" }, { "step": 11, "name": "mitc4-patch-benchmark-tests", "status": "completed" }, - { "step": 12, "name": "assembly-sparse-solver-path", "status": "pending" }, + { "step": 12, "name": "assembly-sparse-solver-path", "status": "completed" }, { "step": 13, "name": "linear-static-workflow", "status": "pending" }, { "step": 14, "name": "stored-reference-regression", "status": "pending" }, { "step": 15, "name": "phase1-evaluator-closeout", "status": "pending" } diff --git a/tests/test_main.cpp b/tests/test_main.cpp index 047fef2..9b26545 100644 --- a/tests/test_main.cpp +++ b/tests/test_main.cpp @@ -230,6 +230,33 @@ fesa::Real singleElementCantileverTipUz(fesa::Real thickness) { return 0.5 * (node2_uz + node3_uz); } +class RecordingSolver final : public fesa::LinearSolver { + public: + explicit RecordingSolver(std::vector solution) : solution_(std::move(solution)) {} + + fesa::SolveResult solve(fesa::DenseMatrix a, std::vector b) const override { + called = true; + captured_a = std::move(a); + captured_b = std::move(b); + return {solution_, {}}; + } + + mutable bool called = false; + mutable fesa::DenseMatrix captured_a; + mutable std::vector captured_b; + + private: + std::vector solution_; +}; + +class FailingSolver final : public fesa::LinearSolver { + public: + fesa::SolveResult solve(fesa::DenseMatrix, std::vector) const override { + return {{}, {fesa::makeDiagnostic(fesa::Severity::Error, "FESA-SINGULAR-SOLVER", + "Injected reduced system singularity", "solver")}}; + } +}; + fesa::Domain singleElementValidationDomain() { fesa::Domain domain; domain.nodes[1] = {1, {0, 0, 0}}; @@ -725,6 +752,96 @@ FESA_TEST(full_vector_reaction_recovery_uses_full_system_quantities) { FESA_CHECK_NEAR(reaction[static_cast(free_uz)], 8.0, 1.0e-15); } +FESA_TEST(reduced_sparse_pattern_is_deterministic_for_phase1_connectivity) { + auto domain = parsedPhase1Domain(); + fesa::DofManager dofs(domain); + const auto pattern = fesa::buildReducedSparsePattern(domain, dofs); + + FESA_CHECK(pattern.equation_count == dofs.freeDofCount()); + FESA_CHECK(pattern.nonzeroCount() == 4); + FESA_CHECK(pattern.contains(0, 0)); + FESA_CHECK(pattern.contains(0, 1)); + FESA_CHECK(pattern.contains(1, 0)); + FESA_CHECK(pattern.contains(1, 1)); + FESA_CHECK(pattern.entries[0].row == 0); + FESA_CHECK(pattern.entries[0].col == 0); + FESA_CHECK(pattern.entries[1].row == 0); + FESA_CHECK(pattern.entries[1].col == 1); + FESA_CHECK(pattern.entries[2].row == 1); + FESA_CHECK(pattern.entries[2].col == 0); + FESA_CHECK(pattern.entries[3].row == 1); + FESA_CHECK(pattern.entries[3].col == 1); +} + +FESA_TEST(assembly_projection_uses_dof_manager_free_indices) { + auto domain = parsedPhase1Domain(); + fesa::DofManager dofs(domain); + fesa::AssemblyResult assembly; + assembly.k_full = fesa::DenseMatrix(dofs.fullDofCount(), dofs.fullDofCount()); + assembly.f_full = std::vector(static_cast(dofs.fullDofCount()), 0.0); + assembly.reduced_pattern = fesa::buildReducedSparsePattern(domain, dofs); + + const auto node2_uz = dofs.fullIndex(2, fesa::Dof::UZ); + const auto node3_uz = dofs.fullIndex(3, fesa::Dof::UZ); + const auto support_uz = dofs.fullIndex(1, fesa::Dof::UZ); + assembly.k_full(node2_uz, node2_uz) = 3.0; + assembly.k_full(node2_uz, node3_uz) = 1.0; + assembly.k_full(node3_uz, node2_uz) = 1.0; + assembly.k_full(node3_uz, node3_uz) = 2.0; + assembly.k_full(support_uz, node2_uz) = 7.0; + assembly.f_full[static_cast(node2_uz)] = -1.0; + assembly.f_full[static_cast(node3_uz)] = -2.0; + + const auto reduced = fesa::projectToReducedSystem(assembly, dofs); + FESA_CHECK(reduced.ok()); + FESA_CHECK(reduced.k.rows() == 2); + FESA_CHECK(reduced.k.cols() == 2); + FESA_CHECK(reduced.free_full_indices == dofs.freeFullIndices()); + FESA_CHECK_NEAR(reduced.k(0, 0), 3.0, 1.0e-15); + FESA_CHECK_NEAR(reduced.k(0, 1), 1.0, 1.0e-15); + FESA_CHECK_NEAR(reduced.k(1, 0), 1.0, 1.0e-15); + FESA_CHECK_NEAR(reduced.k(1, 1), 2.0, 1.0e-15); + FESA_CHECK_NEAR(reduced.f[0], -1.0, 1.0e-15); + FESA_CHECK_NEAR(reduced.f[1], -2.0, 1.0e-15); + + fesa::GaussianEliminationSolver solver; + const auto solved = solver.solve(reduced.k, reduced.f); + FESA_CHECK(solved.ok()); + FESA_CHECK_NEAR(solved.x[0], 0.0, 1.0e-12); + FESA_CHECK_NEAR(solved.x[1], -1.0, 1.0e-12); +} + +FESA_TEST(assembly_preserves_full_space_stiffness_load_and_reduced_pattern) { + auto domain = parsedPhase1Domain(); + fesa::DofManager dofs(domain); + const auto assembly = fesa::assembleSystem(domain, dofs); + FESA_CHECK(!fesa::hasError(assembly.diagnostics)); + FESA_CHECK(assembly.k_full.rows() == dofs.fullDofCount()); + FESA_CHECK(assembly.k_full.cols() == dofs.fullDofCount()); + FESA_CHECK(assembly.f_full.size() == static_cast(dofs.fullDofCount())); + FESA_CHECK(assembly.reduced_pattern.equation_count == dofs.freeDofCount()); + FESA_CHECK(assembly.reduced_pattern.nonzeroCount() == 4); + + const auto node2_uz = dofs.fullIndex(2, fesa::Dof::UZ); + const auto node3_uz = dofs.fullIndex(3, fesa::Dof::UZ); + FESA_CHECK_NEAR(assembly.f_full[static_cast(node2_uz)], -1.0, 1.0e-15); + FESA_CHECK_NEAR(assembly.f_full[static_cast(node3_uz)], -1.0, 1.0e-15); + for (fesa::LocalIndex i = 0; i < assembly.k_full.rows(); ++i) { + for (fesa::LocalIndex j = 0; j < assembly.k_full.cols(); ++j) { + FESA_CHECK_NEAR(assembly.k_full(i, j), assembly.k_full(j, i), 1.0e-8); + } + } + + const auto reduced = fesa::projectToReducedSystem(assembly, dofs); + FESA_CHECK(reduced.ok()); + const auto solved = fesa::GaussianEliminationSolver{}.solve(reduced.k, reduced.f); + FESA_CHECK(solved.ok()); + const auto residual = reduced.k.multiply(solved.x); + for (std::size_t i = 0; i < residual.size(); ++i) { + FESA_CHECK_NEAR(residual[i], reduced.f[i], 1.0e-8); + } +} + FESA_TEST(gaussian_solver_solves_and_diagnoses_singular_systems) { fesa::DenseMatrix a(2, 2); a(0, 0) = 2.0; @@ -1573,6 +1690,9 @@ FESA_TEST(linear_static_analysis_solves_u_and_recovers_full_vector_rf) { auto result = analysis.run(domain); FESA_CHECK(result.ok()); FESA_CHECK(result.state.converged); + FESA_CHECK(result.state.f_internal_full.size() == result.state.u_full.size()); + FESA_CHECK(result.state.f_external_full.size() == result.state.u_full.size()); + FESA_CHECK(result.state.reaction_full.size() == result.state.u_full.size()); FESA_CHECK(result.result_file.steps.size() == 1); const auto& frame = result.result_file.steps[0].frames[0]; FESA_CHECK(frame.field_outputs.count("U") == 1); @@ -1585,6 +1705,39 @@ FESA_TEST(linear_static_analysis_solves_u_and_recovers_full_vector_rf) { FESA_CHECK_NEAR(total_rf_z, 2.0, 1.0e-8); } +FESA_TEST(linear_static_analysis_uses_solver_adapter_and_reconstructs_full_vectors) { + auto domain = parsedPhase1Domain(); + RecordingSolver solver({0.25, -0.50}); + fesa::LinearStaticAnalysis analysis(&solver); + const auto result = analysis.run(domain); + FESA_CHECK(result.ok()); + FESA_CHECK(solver.called); + FESA_CHECK(solver.captured_a.rows() == 2); + FESA_CHECK(solver.captured_a.cols() == 2); + FESA_CHECK(solver.captured_b.size() == 2); + FESA_CHECK_NEAR(solver.captured_b[0], -1.0, 1.0e-15); + FESA_CHECK_NEAR(solver.captured_b[1], -1.0, 1.0e-15); + + fesa::DofManager dofs(domain); + FESA_CHECK_NEAR(result.state.u_full[static_cast(dofs.fullIndex(2, fesa::Dof::UZ))], 0.25, 1.0e-15); + FESA_CHECK_NEAR(result.state.u_full[static_cast(dofs.fullIndex(3, fesa::Dof::UZ))], -0.50, 1.0e-15); + FESA_CHECK_NEAR(result.state.u_full[static_cast(dofs.fullIndex(1, fesa::Dof::UZ))], 0.0, 1.0e-15); + for (std::size_t i = 0; i < result.state.reaction_full.size(); ++i) { + FESA_CHECK_NEAR(result.state.reaction_full[i], + result.state.f_internal_full[i] - result.state.f_external_full[i], 1.0e-10); + } +} + +FESA_TEST(linear_static_analysis_propagates_solver_singular_diagnostic) { + auto domain = parsedPhase1Domain(); + FailingSolver solver; + fesa::LinearStaticAnalysis analysis(&solver); + const auto result = analysis.run(domain); + FESA_CHECK(!result.ok()); + FESA_CHECK(!result.state.converged); + FESA_CHECK(fesa::containsDiagnostic(result.diagnostics, "FESA-SINGULAR-SOLVER")); +} + int main() { int failed = 0; for (const auto& test : registry()) {