From 9a916960146c2388d01dc808b90b1a35be55c1fe Mon Sep 17 00:00:00 2001 From: NINI Date: Tue, 5 May 2026 23:33:01 +0900 Subject: [PATCH] refactor: extract assembly analysis workflow --- CMakeLists.txt | 12 + PLAN.md | 8 +- PROGRESS.md | 41 ++- include/fesa/Analysis/Analysis.hpp | 1 + .../fesa/Analysis/LinearStaticAnalysis.hpp | 120 ++++++++ include/fesa/Assembly/Assembly.hpp | 1 + include/fesa/Assembly/AssemblySystem.hpp | 168 +++++++++++ include/fesa/fesa.hpp | 279 +----------------- .../1-structure-alignment-refactor/index.json | 2 +- tests/test_analysis_module_includes.cpp | 129 ++++++++ tests/test_assembly_module_includes.cpp | 83 ++++++ 11 files changed, 559 insertions(+), 285 deletions(-) create mode 100644 include/fesa/Analysis/LinearStaticAnalysis.hpp create mode 100644 include/fesa/Assembly/AssemblySystem.hpp create mode 100644 tests/test_analysis_module_includes.cpp create mode 100644 tests/test_assembly_module_includes.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index d9b5697..9f47b78 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -38,6 +38,12 @@ target_link_libraries(fesa_element_module_tests PRIVATE fesa_core) add_executable(fesa_mitc4_stiffness_module_tests tests/test_mitc4_stiffness_module_includes.cpp) target_link_libraries(fesa_mitc4_stiffness_module_tests PRIVATE fesa_core) +add_executable(fesa_assembly_module_tests tests/test_assembly_module_includes.cpp) +target_link_libraries(fesa_assembly_module_tests PRIVATE fesa_core) + +add_executable(fesa_analysis_module_tests tests/test_analysis_module_includes.cpp) +target_link_libraries(fesa_analysis_module_tests PRIVATE fesa_core) + if(MSVC) target_compile_options(fesa_core PRIVATE /W4 /permissive-) target_compile_options(fesa_tests PRIVATE /W4 /permissive-) @@ -47,6 +53,8 @@ if(MSVC) target_compile_options(fesa_results_module_tests PRIVATE /W4 /permissive-) target_compile_options(fesa_element_module_tests PRIVATE /W4 /permissive-) target_compile_options(fesa_mitc4_stiffness_module_tests PRIVATE /W4 /permissive-) + target_compile_options(fesa_assembly_module_tests PRIVATE /W4 /permissive-) + target_compile_options(fesa_analysis_module_tests PRIVATE /W4 /permissive-) else() target_compile_options(fesa_core PRIVATE -Wall -Wextra -Wpedantic) target_compile_options(fesa_tests PRIVATE -Wall -Wextra -Wpedantic) @@ -56,6 +64,8 @@ else() target_compile_options(fesa_results_module_tests PRIVATE -Wall -Wextra -Wpedantic) target_compile_options(fesa_element_module_tests PRIVATE -Wall -Wextra -Wpedantic) target_compile_options(fesa_mitc4_stiffness_module_tests PRIVATE -Wall -Wextra -Wpedantic) + target_compile_options(fesa_assembly_module_tests PRIVATE -Wall -Wextra -Wpedantic) + target_compile_options(fesa_analysis_module_tests PRIVATE -Wall -Wextra -Wpedantic) endif() add_test(NAME fesa_tests COMMAND fesa_tests) @@ -65,3 +75,5 @@ add_test(NAME fesa_io_module_tests COMMAND fesa_io_module_tests) add_test(NAME fesa_results_module_tests COMMAND fesa_results_module_tests) add_test(NAME fesa_element_module_tests COMMAND fesa_element_module_tests) add_test(NAME fesa_mitc4_stiffness_module_tests COMMAND fesa_mitc4_stiffness_module_tests) +add_test(NAME fesa_assembly_module_tests COMMAND fesa_assembly_module_tests) +add_test(NAME fesa_analysis_module_tests COMMAND fesa_analysis_module_tests) diff --git a/PLAN.md b/PLAN.md index 6a08fd1..9e38c56 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 -Execute the Phase 1 structure-alignment refactor in `phases/1-structure-alignment-refactor`, continuing with P1A-08 Assembly and Analysis workflow extraction. P1A-00 completed the architecture drift audit, P1A-01 created the module scaffold, P1A-02 extracted Core/Util plus Phase 1 Boundary/Load/Property model ownership, P1A-03 extracted Math primitives plus the solver adapter boundary, P1A-04 extracted the Abaqus Phase 1 parser into IO, P1A-05 extracted Results/reference comparison, P1A-06 extracted MITC4 geometry/strain helpers into Element, and P1A-07 extracted MITC4 material/stiffness helpers into Material and Element without changing solver behavior. This phase must align the current monolithic `include/fesa/fesa.hpp` implementation with the module ownership model in `docs/ARCHITECTURE.md` without changing solver behavior. Product-level Phase 1 reference gaps R-010 and R-013 remain open and must not be hidden by the refactor. +Execute the Phase 1 structure-alignment refactor in `phases/1-structure-alignment-refactor`, continuing with P1A-09 independent architecture evaluator closeout. P1A-00 completed the architecture drift audit, P1A-01 created the module scaffold, P1A-02 extracted Core/Util plus Phase 1 Boundary/Load/Property model ownership, P1A-03 extracted Math primitives plus the solver adapter boundary, P1A-04 extracted the Abaqus Phase 1 parser into IO, P1A-05 extracted Results/reference comparison, P1A-06 extracted MITC4 geometry/strain helpers into Element, P1A-07 extracted MITC4 material/stiffness helpers into Material and Element, and P1A-08 extracted Assembly and Analysis workflow without changing solver behavior. `include/fesa/fesa.hpp` is now an include-only facade, but R-014 remains open until P1A-09 independently accepts the final architecture alignment. Product-level Phase 1 reference gaps R-010 and R-013 remain open and must not be hidden by the refactor. ## Required Reading For New Agents 1. `AGENTS.md` @@ -37,7 +37,7 @@ Execute the Phase 1 structure-alignment refactor in `phases/1-structure-alignmen ## Phase Files - Active phase directory: `phases/1-structure-alignment-refactor` - Execute with: `python scripts/execute.py 1-structure-alignment-refactor` -- Step numbering is zero-based. `step0.md` is complete and wrote `phases/1-structure-alignment-refactor/step0-architecture-map.md`; `step1.md` is complete and created module scaffold headers, source directories, CMake source discovery, and umbrella compatibility smoke coverage; `step2.md` is complete and extracted Core/Util domain, diagnostics, DofManager ownership, AnalysisModel/AnalysisState, and Phase 1 Boundary/Load/Property model ownership; `step3.md` is complete and extracted Math primitives, sparse pattern data, dense matrix support, and solver adapter boundary; `step4.md` is complete and extracted the Abaqus parser into IO; `step5.md` is complete and extracted Results/reference comparison code; `step6.md` is complete and extracted MITC4 geometry/strain helpers; `step7.md` is complete and extracted MITC4 material/stiffness helpers; `step8.md` extracts Assembly and Analysis workflow; `step9.md` is the independent architecture evaluator closeout. +- Step numbering is zero-based. `step0.md` is complete and wrote `phases/1-structure-alignment-refactor/step0-architecture-map.md`; `step1.md` is complete and created module scaffold headers, source directories, CMake source discovery, and umbrella compatibility smoke coverage; `step2.md` is complete and extracted Core/Util domain, diagnostics, DofManager ownership, AnalysisModel/AnalysisState, and Phase 1 Boundary/Load/Property model ownership; `step3.md` is complete and extracted Math primitives, sparse pattern data, dense matrix support, and solver adapter boundary; `step4.md` is complete and extracted the Abaqus parser into IO; `step5.md` is complete and extracted Results/reference comparison code; `step6.md` is complete and extracted MITC4 geometry/strain helpers; `step7.md` is complete and extracted MITC4 material/stiffness helpers; `step8.md` is complete and extracted Assembly and Analysis workflow; `step9.md` is the independent architecture evaluator closeout. - Completed phase directory: `phases/1-linear-static-mitc4-rebaseline` - Historical execution command: `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; `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; `step13.md` is complete and revalidated active AnalysisModel construction plus input-to-AnalysisState-to-U/RF result workflow; `step14.md` is complete and added the first stored Abaqus displacement regression for `quad_02_phase1`; `step15.md` is complete and recorded the independent evaluator closeout in `phases/1-linear-static-mitc4-rebaseline/step15-evaluator-report.md`. @@ -65,7 +65,7 @@ This phase is an architecture-preserving refactor. It must not change Phase 1 so | P1A-05 | completed | generator | Extract Results model, writer boundary, CSV loader, and reference comparator. | P1A-02, P1A-04 | `U`/`RF` schema and `quad_02_phase1` regression unchanged | | P1A-06 | completed | generator | Extract MITC4 geometry, director, strain, and tying helpers into Element. | P1A-03 | Geometry/strain tests and formulation signs unchanged | | P1A-07 | completed | generator | Extract MITC4 material, integration, stiffness, drilling, and internal-force helpers. | P1A-06 | Patch, drilling, stiffness, and locking-sensitivity tests unchanged | -| P1A-08 | pending | generator | Extract Assembly and Analysis workflow. | P1A-02, P1A-03, P1A-05, P1A-07 | Full-vector RF, solver injection, and end-to-end reference regression unchanged | +| P1A-08 | completed | generator | Extract Assembly and Analysis workflow. | P1A-02, P1A-03, P1A-05, P1A-07 | Full-vector RF, solver injection, and end-to-end reference regression unchanged | | P1A-09 | pending | evaluator | Independently evaluate final architecture alignment. | P1A-08 | `src/` ownership matches `ARCHITECTURE.md`; umbrella header is facade only | ## Phase 1 Definition Of Done @@ -152,7 +152,7 @@ Current reference state: - `references/quad_02_phase1.inp` is the accepted normalized Phase 1-compatible derivative input for the `quad_02` S4 reference pair. Required reference additions or decisions: -- Add `*_reactions.csv` or explicitly use internal equilibrium tests for Phase 1 `RF` until Abaqus RF output is available. +- Onboard any provided `*_reactionforces.csv` or `*_reactions.csv` artifact with a documented schema/tolerance, or explicitly use internal equilibrium tests for Phase 1 `RF` until Abaqus RF CSV is accepted. - Add more small cases until Phase 1 can pass one single-element case, one simple multi-element plate/shell case, and one curved shell benchmark. ## Phase 1 Risk Controls diff --git a/PROGRESS.md b/PROGRESS.md index 1966034..43a8f3c 100644 --- a/PROGRESS.md +++ b/PROGRESS.md @@ -13,10 +13,45 @@ 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 completed rebaseline execution path in `phases/1-linear-static-mitc4-rebaseline`. Steps 0 through 15 are complete, and P1R-15 recorded a pass-with-documented-gaps evaluator closeout. The follow-up architecture refactor phase in `phases/1-structure-alignment-refactor` is underway because remaining Assembly and Analysis workflow code still lives in `include/fesa/fesa.hpp` instead of the module directories documented in `docs/ARCHITECTURE.md`; P1A-00 through P1A-07 are complete, so the next step is P1A-08 Assembly and Analysis workflow extraction. `quad_02_phase1.inp` is 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, full-vector internal/reaction force state, active AnalysisModel construction, input-to-AnalysisState-to-U/RF result workflow, and the first stored Abaqus displacement regression have been revalidated. Full PRD Phase 1 completion still depends on the open architecture/reference gaps R-014, R-010, and R-013. The old `phases/1-linear-static-mitc4` path is historical and superseded after the MITC4 formulation reset. +Phase 1 has a completed rebaseline execution path in `phases/1-linear-static-mitc4-rebaseline`. Steps 0 through 15 are complete, and P1R-15 recorded a pass-with-documented-gaps evaluator closeout. The follow-up architecture refactor phase in `phases/1-structure-alignment-refactor` is underway; P1A-00 through P1A-08 are complete, `include/fesa/fesa.hpp` is now an include-only facade, and the next step is P1A-09 independent architecture evaluator closeout. `quad_02_phase1.inp` is 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, full-vector internal/reaction force state, active AnalysisModel construction, input-to-AnalysisState-to-U/RF result workflow, and the first stored Abaqus displacement regression have been revalidated. Full PRD Phase 1 completion still depends on the open architecture/reference gaps R-014, R-010, and R-013. The old `phases/1-linear-static-mitc4` path is historical and superseded after the MITC4 formulation reset. ## Completed Work +### 2026-05-05 - P1A-08 Assembly Analysis extraction completed +Author: Codex + +Changed files: +- `CMakeLists.txt` +- `include/fesa/Analysis/Analysis.hpp` +- `include/fesa/Analysis/LinearStaticAnalysis.hpp` +- `include/fesa/Assembly/Assembly.hpp` +- `include/fesa/Assembly/AssemblySystem.hpp` +- `include/fesa/fesa.hpp` +- `tests/test_analysis_module_includes.cpp` +- `tests/test_assembly_module_includes.cpp` +- `phases/1-structure-alignment-refactor/index.json` +- `PLAN.md` +- `PROGRESS.md` + +Summary: +- Extracted `buildReducedSparsePattern`, `recoverFullReaction`, `AssemblyResult`, `ReducedSystem`, `assembleSystem`, and `projectToReducedSystem` into `include/fesa/Assembly/AssemblySystem.hpp`. +- Extracted `AnalysisResult`, `Analysis`, `LinearStaticAnalysis`, and `runLinearStaticInputString` into `include/fesa/Analysis/LinearStaticAnalysis.hpp`. +- Kept `AnalysisState` in `include/fesa/Core/AnalysisState.hpp` because `docs/ARCHITECTURE.md` and the P1A-00 migration map place mutable analysis state under Core ownership. +- Updated Assembly and Analysis facade headers so direct module includes expose the relocated workflow without including `fesa/fesa.hpp`. +- Reduced `include/fesa/fesa.hpp` to an include-only umbrella facade with no production implementation body. +- Preserved full-space stiffness/load preservation, constrained/free reduced projection, solver adapter injection, deterministic default Gaussian solver, `RF = K_full * U_full - F_full`, and step/frame `U`/`RF` result writing behavior. +- No parser subset, MITC4 formulation, numerical convention, result schema, reference tolerance, sparse storage, HDF5, MKL, TBB, nonlinear, dynamic, pressure-load, or RBE behavior was added. + +Verification: +- First ran `python scripts\validate_workspace.py` after adding direct Assembly and Analysis include tests; it failed as expected because the module facades did not yet expose the Assembly and Analysis symbols. +- After extraction, `python scripts\validate_workspace.py` configured CMake, built `fesa_core`, `fesa_tests`, `fesa_core_module_tests`, `fesa_math_module_tests`, `fesa_io_module_tests`, `fesa_results_module_tests`, `fesa_element_module_tests`, `fesa_mitc4_stiffness_module_tests`, `fesa_assembly_module_tests`, and `fesa_analysis_module_tests`, and ran CTest successfully. +- CTest result: 9 test executables passed. + +Follow-up: +- Continue with P1A-09 independent architecture evaluator closeout. +- Keep R-014 open until P1A-09 independently accepts the final architecture alignment. +- Keep R-010 and R-013 open; this refactor does not onboard reaction CSV artifacts or add additional stored reference cases. + ### 2026-05-05 - P1A-07 MITC4 material stiffness extraction completed Author: Codex @@ -1169,8 +1204,8 @@ Verification: - `python scripts/validate_workspace.py` ran, but reported no configured validation commands. ## Known Blockers -- Phase 1 architecture is not yet accepted: remaining Assembly and Analysis workflow code still lives in `include/fesa/fesa.hpp` instead of the module layout documented in `docs/ARCHITECTURE.md`. -- No reaction-force reference artifact exists yet under `references/`. +- Phase 1 architecture is not yet accepted until P1A-09 independently evaluates the final module alignment and records pass/fail closeout. +- A reaction-force CSV may be present as untracked local reference input, but no reaction-force artifact has been onboarded with documented schema, tolerance, and automated comparison yet. - The PRD target of three stored Phase 1 reference cases is not yet satisfied; only `quad_02_phase1` is an active stored displacement regression. - The current initial `quad_01.inp` reference contains `S4R`, `Part/Assembly/Instance`, `*Density`, and `NLGEOM=YES`, so it is not a Phase 1 parser acceptance case as-is. diff --git a/include/fesa/Analysis/Analysis.hpp b/include/fesa/Analysis/Analysis.hpp index dad2066..6cb6227 100644 --- a/include/fesa/Analysis/Analysis.hpp +++ b/include/fesa/Analysis/Analysis.hpp @@ -1,5 +1,6 @@ #pragma once +#include "fesa/Analysis/LinearStaticAnalysis.hpp" #include "fesa/ModuleInfo.hpp" namespace fesa::module { diff --git a/include/fesa/Analysis/LinearStaticAnalysis.hpp b/include/fesa/Analysis/LinearStaticAnalysis.hpp new file mode 100644 index 0000000..daf23fb --- /dev/null +++ b/include/fesa/Analysis/LinearStaticAnalysis.hpp @@ -0,0 +1,120 @@ +#pragma once + +#include "fesa/Assembly/Assembly.hpp" +#include "fesa/Core/Core.hpp" +#include "fesa/IO/IO.hpp" +#include "fesa/Math/Math.hpp" +#include "fesa/Results/Results.hpp" +#include "fesa/Util/Diagnostics.hpp" + +#include +#include + +namespace fesa { + +struct AnalysisResult { + AnalysisModel model; + AnalysisState state; + ResultFile result_file; + std::vector diagnostics; + + bool ok() const { + return !hasError(diagnostics); + } +}; + +class Analysis { + public: + virtual ~Analysis() = default; + + AnalysisResult run(const Domain& domain) const { + AnalysisResult result; + initialize(domain, result); + if (hasError(result.diagnostics)) { + return result; + } + solve(domain, result); + return result; + } + + protected: + virtual void initialize(const Domain& domain, AnalysisResult& result) const { + auto diagnostics = validateDomain(domain); + result.diagnostics.insert(result.diagnostics.end(), diagnostics.begin(), diagnostics.end()); + } + + virtual void solve(const Domain& domain, AnalysisResult& result) const = 0; +}; + +class LinearStaticAnalysis final : public Analysis { + public: + explicit LinearStaticAnalysis(const LinearSolver* solver = nullptr) : solver_(solver) {} + + protected: + void solve(const Domain& domain, AnalysisResult& result) const override { + result.model = buildLinearStaticAnalysisModel(domain); + result.diagnostics.insert(result.diagnostics.end(), result.model.diagnostics.begin(), result.model.diagnostics.end()); + if (hasError(result.diagnostics)) { + return; + } + DofManager dofs(domain); + 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; + } + AssemblyResult assembly = assembleSystem(domain, dofs); + result.diagnostics.insert(result.diagnostics.end(), assembly.diagnostics.begin(), assembly.diagnostics.end()); + if (hasError(result.diagnostics)) { + return; + } + const auto reduced = projectToReducedSystem(assembly, dofs); + result.diagnostics.insert(result.diagnostics.end(), reduced.diagnostics.begin(), reduced.diagnostics.end()); + if (hasError(result.diagnostics)) { + return; + } + 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, result.model, 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; +}; + +inline AnalysisResult runLinearStaticInputString(const std::string& text, + const std::string& source_name = "", + const LinearSolver* solver = nullptr) { + AbaqusInputParser parser; + const auto parsed = parser.parseString(text, source_name); + if (!parsed.ok()) { + AnalysisResult result; + result.diagnostics = parsed.diagnostics; + return result; + } + LinearStaticAnalysis analysis(solver); + return analysis.run(parsed.domain); +} + +} // namespace fesa diff --git a/include/fesa/Assembly/Assembly.hpp b/include/fesa/Assembly/Assembly.hpp index a0c7fff..dd8703f 100644 --- a/include/fesa/Assembly/Assembly.hpp +++ b/include/fesa/Assembly/Assembly.hpp @@ -1,5 +1,6 @@ #pragma once +#include "fesa/Assembly/AssemblySystem.hpp" #include "fesa/ModuleInfo.hpp" namespace fesa::module { diff --git a/include/fesa/Assembly/AssemblySystem.hpp b/include/fesa/Assembly/AssemblySystem.hpp new file mode 100644 index 0000000..2b63836 --- /dev/null +++ b/include/fesa/Assembly/AssemblySystem.hpp @@ -0,0 +1,168 @@ +#pragma once + +#include "fesa/Core/Core.hpp" +#include "fesa/Element/Element.hpp" +#include "fesa/Load/Load.hpp" +#include "fesa/Math/Math.hpp" +#include "fesa/Property/Property.hpp" +#include "fesa/Util/Diagnostics.hpp" + +#include +#include +#include +#include +#include + +namespace fesa { + +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; +} + +inline std::vector recoverFullReaction(const DenseMatrix& k_full, + const std::vector& u_full, + const std::vector& f_full) { + if (k_full.rows() != k_full.cols() || static_cast(u_full.size()) != k_full.cols() || + static_cast(f_full.size()) != k_full.rows()) { + throw std::runtime_error("full reaction size mismatch"); + } + std::vector reaction = k_full.multiply(u_full); + for (std::size_t i = 0; i < reaction.size(); ++i) { + reaction[i] -= f_full[i]; + } + return reaction; +} + +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; + 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) { + result.diagnostics.push_back( + {Severity::Error, "FESA-ASSEMBLY-MISSING-PROPERTY", "Element has no shell section", {}}); + continue; + } + const auto material_it = domain.materials.find(Domain::key(section->material)); + if (material_it == domain.materials.end()) { + result.diagnostics.push_back( + {Severity::Error, "FESA-ASSEMBLY-MISSING-MATERIAL", "Element material is missing", {}}); + continue; + } + std::array coordinates{}; + for (std::size_t i = 0; i < 4; ++i) { + coordinates[i] = domain.nodes.at(element.node_ids[i]).coordinates; + } + 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, 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, *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; +} + +} // namespace fesa diff --git a/include/fesa/fesa.hpp b/include/fesa/fesa.hpp index a0c0ad7..42cea9a 100644 --- a/include/fesa/fesa.hpp +++ b/include/fesa/fesa.hpp @@ -1,5 +1,7 @@ #pragma once +#include "fesa/Analysis/Analysis.hpp" +#include "fesa/Assembly/Assembly.hpp" #include "fesa/Boundary/Boundary.hpp" #include "fesa/Core/Core.hpp" #include "fesa/Element/Element.hpp" @@ -11,280 +13,3 @@ #include "fesa/Property/Property.hpp" #include "fesa/Results/Results.hpp" #include "fesa/Util/Util.hpp" - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -namespace fesa { - -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; -} - -inline std::vector recoverFullReaction(const DenseMatrix& k_full, - const std::vector& u_full, - const std::vector& f_full) { - if (k_full.rows() != k_full.cols() || static_cast(u_full.size()) != k_full.cols() || - static_cast(f_full.size()) != k_full.rows()) { - throw std::runtime_error("full reaction size mismatch"); - } - std::vector reaction = k_full.multiply(u_full); - for (std::size_t i = 0; i < reaction.size(); ++i) { - reaction[i] -= f_full[i]; - } - return reaction; -} - -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; - 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) { - result.diagnostics.push_back( - {Severity::Error, "FESA-ASSEMBLY-MISSING-PROPERTY", "Element has no shell section", {}}); - continue; - } - const auto material_it = domain.materials.find(Domain::key(section->material)); - if (material_it == domain.materials.end()) { - result.diagnostics.push_back( - {Severity::Error, "FESA-ASSEMBLY-MISSING-MATERIAL", "Element material is missing", {}}); - continue; - } - std::array coordinates{}; - for (std::size_t i = 0; i < 4; ++i) { - coordinates[i] = domain.nodes.at(element.node_ids[i]).coordinates; - } - 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, 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, *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; -} - -struct AnalysisResult { - AnalysisModel model; - AnalysisState state; - ResultFile result_file; - std::vector diagnostics; - - bool ok() const { - return !hasError(diagnostics); - } -}; - -class Analysis { - public: - virtual ~Analysis() = default; - - AnalysisResult run(const Domain& domain) const { - AnalysisResult result; - initialize(domain, result); - if (hasError(result.diagnostics)) { - return result; - } - solve(domain, result); - return result; - } - - protected: - virtual void initialize(const Domain& domain, AnalysisResult& result) const { - auto diagnostics = validateDomain(domain); - result.diagnostics.insert(result.diagnostics.end(), diagnostics.begin(), diagnostics.end()); - } - - virtual void solve(const Domain& domain, AnalysisResult& result) const = 0; -}; - -class LinearStaticAnalysis final : public Analysis { - public: - explicit LinearStaticAnalysis(const LinearSolver* solver = nullptr) : solver_(solver) {} - - protected: - void solve(const Domain& domain, AnalysisResult& result) const override { - result.model = buildLinearStaticAnalysisModel(domain); - result.diagnostics.insert(result.diagnostics.end(), result.model.diagnostics.begin(), result.model.diagnostics.end()); - if (hasError(result.diagnostics)) { - return; - } - DofManager dofs(domain); - 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; - } - AssemblyResult assembly = assembleSystem(domain, dofs); - result.diagnostics.insert(result.diagnostics.end(), assembly.diagnostics.begin(), assembly.diagnostics.end()); - if (hasError(result.diagnostics)) { - return; - } - const auto reduced = projectToReducedSystem(assembly, dofs); - result.diagnostics.insert(result.diagnostics.end(), reduced.diagnostics.begin(), reduced.diagnostics.end()); - if (hasError(result.diagnostics)) { - return; - } - 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, result.model, 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; -}; - -inline AnalysisResult runLinearStaticInputString(const std::string& text, - const std::string& source_name = "", - const LinearSolver* solver = nullptr) { - AbaqusInputParser parser; - const auto parsed = parser.parseString(text, source_name); - if (!parsed.ok()) { - AnalysisResult result; - result.diagnostics = parsed.diagnostics; - return result; - } - LinearStaticAnalysis analysis(solver); - return analysis.run(parsed.domain); -} - -} // namespace fesa diff --git a/phases/1-structure-alignment-refactor/index.json b/phases/1-structure-alignment-refactor/index.json index dbaf8c9..5040914 100644 --- a/phases/1-structure-alignment-refactor/index.json +++ b/phases/1-structure-alignment-refactor/index.json @@ -10,7 +10,7 @@ { "step": 5, "name": "results-reference-extraction", "status": "completed" }, { "step": 6, "name": "mitc4-geometry-strain-extraction", "status": "completed" }, { "step": 7, "name": "mitc4-material-stiffness-extraction", "status": "completed" }, - { "step": 8, "name": "assembly-analysis-extraction", "status": "pending" }, + { "step": 8, "name": "assembly-analysis-extraction", "status": "completed" }, { "step": 9, "name": "architecture-evaluator-closeout", "status": "pending" } ] } diff --git a/tests/test_analysis_module_includes.cpp b/tests/test_analysis_module_includes.cpp new file mode 100644 index 0000000..b4c12b5 --- /dev/null +++ b/tests/test_analysis_module_includes.cpp @@ -0,0 +1,129 @@ +#include "fesa/Analysis/Analysis.hpp" + +#include +#include +#include +#include +#include + +namespace { + +void check(bool value, const char* message) { + if (!value) { + throw std::runtime_error(message); + } +} + +void checkNear(fesa::Real actual, fesa::Real expected, fesa::Real tolerance, const char* message) { + if (std::fabs(actual - expected) > tolerance) { + throw std::runtime_error(message); + } +} + +std::string phase1Input() { + return R"inp( +*Node +1, 0, 0, 0 +2, 1, 0, 0 +3, 1, 1, 0 +4, 0, 1, 0 +*Element, type=S4, elset=EALL +1, 1, 2, 3, 4 +*Nset, nset=LEFT +1, 4 +*Nset, nset=RIGHT +2, 3 +*Elset, elset=EALL +1 +*Material, name=STEEL +*Elastic +1000.0, 0.3 +*Shell Section, elset=EALL, material=STEEL +0.1 +*Boundary +LEFT, 1, 6, 0 +RIGHT, 1, 2, 0 +RIGHT, 4, 6, 0 +*Cload +2, 3, -1 +3, 3, -1 +*Step, name=Step-1 +*Static +*End Step +)inp"; +} + +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_; +}; + +} // namespace + +int main() { + const auto workflow = fesa::runLinearStaticInputString(phase1Input(), "analysis-module.inp"); + check(workflow.ok(), "linear static input workflow should remain valid"); + check(workflow.model.ok(), "analysis model should remain valid"); + check(workflow.model.step.name == "Step-1", "analysis step name changed"); + check(workflow.state.converged, "analysis convergence flag changed"); + check(workflow.state.u_full.size() == 24, "full displacement size changed"); + check(workflow.state.f_external_full.size() == 24, "full external-force size changed"); + check(workflow.state.f_internal_full.size() == 24, "full internal-force size changed"); + check(workflow.state.reaction_full.size() == 24, "full reaction size changed"); + check(workflow.result_file.steps.size() == 1, "result step count changed"); + const auto& frame = workflow.result_file.steps.front().frames.front(); + check(frame.field_outputs.count("U") == 1, "U field output missing"); + check(frame.field_outputs.count("RF") == 1, "RF field output missing"); + check(frame.field_outputs.at("U").component_labels == fesa::displacementComponentLabels(), "U labels changed"); + check(frame.field_outputs.at("RF").component_labels == fesa::reactionComponentLabels(), "RF labels changed"); + + fesa::Real total_rf_z = 0.0; + for (const auto& values : frame.field_outputs.at("RF").values) { + total_rf_z += values[2]; + } + checkNear(total_rf_z, 2.0, 1.0e-8, "full-vector RF balance changed"); + + fesa::AbaqusInputParser parser; + const auto parsed = parser.parseString(phase1Input()); + check(parsed.ok(), "phase1 analysis input parse changed"); + RecordingSolver solver({0.25, -0.50}); + fesa::LinearStaticAnalysis analysis(&solver); + const auto injected = analysis.run(parsed.domain); + check(injected.ok(), "solver-injected analysis should remain valid"); + check(solver.called, "linear solver adapter injection changed"); + check(solver.captured_a.rows() == 2 && solver.captured_a.cols() == 2, "captured reduced stiffness size changed"); + check(solver.captured_b.size() == 2, "captured reduced load size changed"); + + const fesa::DofManager dofs(parsed.domain); + checkNear(injected.state.u_full[static_cast(dofs.fullIndex(2, fesa::Dof::UZ))], 0.25, 1.0e-15, + "node 2 reconstructed UZ changed"); + checkNear(injected.state.u_full[static_cast(dofs.fullIndex(3, fesa::Dof::UZ))], -0.50, 1.0e-15, + "node 3 reconstructed UZ changed"); + for (std::size_t i = 0; i < injected.state.reaction_full.size(); ++i) { + checkNear(injected.state.reaction_full[i], + injected.state.f_internal_full[i] - injected.state.f_external_full[i], + 1.0e-10, + "full-vector reaction formula changed"); + } + + const auto parse_error = fesa::runLinearStaticInputString("*Part, name=P\n", "unsupported.inp"); + check(!parse_error.ok(), "parse errors should still be routed through analysis result"); + check(fesa::containsDiagnostic(parse_error.diagnostics, "FESA-PARSE-UNSUPPORTED-KEYWORD"), + "parse diagnostic routing changed"); + + return 0; +} diff --git a/tests/test_assembly_module_includes.cpp b/tests/test_assembly_module_includes.cpp new file mode 100644 index 0000000..3bbf93d --- /dev/null +++ b/tests/test_assembly_module_includes.cpp @@ -0,0 +1,83 @@ +#include "fesa/Assembly/Assembly.hpp" + +#include +#include +#include + +namespace { + +void check(bool value, const char* message) { + if (!value) { + throw std::runtime_error(message); + } +} + +void checkNear(fesa::Real actual, fesa::Real expected, fesa::Real tolerance, const char* message) { + if (std::fabs(actual - expected) > tolerance) { + throw std::runtime_error(message); + } +} + +fesa::Domain assemblyDomain() { + fesa::Domain domain; + domain.nodes[1] = {1, {0.0, 0.0, 0.0}}; + domain.nodes[2] = {2, {1.0, 0.0, 0.0}}; + domain.nodes[3] = {3, {1.0, 1.0, 0.0}}; + domain.nodes[4] = {4, {0.0, 1.0, 0.0}}; + domain.elements[1] = {1, fesa::ElementType::MITC4, {1, 2, 3, 4}, "EALL"}; + domain.element_sets["eall"] = {"EALL", {1}}; + domain.node_sets["left"] = {"LEFT", {1, 4}}; + domain.node_sets["right"] = {"RIGHT", {2, 3}}; + domain.materials["steel"] = {"STEEL", 1000.0, 0.30}; + domain.shell_sections.push_back({"EALL", "STEEL", 0.1}); + domain.boundary_conditions.push_back({"LEFT", 1, 6, 0.0}); + domain.boundary_conditions.push_back({"RIGHT", 1, 2, 0.0}); + domain.boundary_conditions.push_back({"RIGHT", 4, 6, 0.0}); + domain.loads.push_back({"2", 3, -1.0}); + domain.loads.push_back({"3", 3, -1.0}); + return domain; +} + +} // namespace + +int main() { + const auto domain = assemblyDomain(); + const fesa::DofManager dofs(domain); + + const auto pattern = fesa::buildReducedSparsePattern(domain, dofs); + check(pattern.equation_count == dofs.freeDofCount(), "reduced sparse pattern equation count changed"); + check(pattern.nonzeroCount() > 0, "reduced sparse pattern should remain non-empty"); + + const auto assembly = fesa::assembleSystem(domain, dofs); + check(assembly.ok(), "assembly should remain valid"); + check(assembly.k_full.rows() == dofs.fullDofCount(), "full stiffness row count changed"); + check(assembly.k_full.cols() == dofs.fullDofCount(), "full stiffness column count changed"); + check(static_cast(assembly.f_full.size()) == dofs.fullDofCount(), "full load size changed"); + checkNear(assembly.f_full[static_cast(dofs.fullIndex(2, fesa::Dof::UZ))], -1.0, 1.0e-15, + "node 2 UZ load changed"); + checkNear(assembly.f_full[static_cast(dofs.fullIndex(3, fesa::Dof::UZ))], -1.0, 1.0e-15, + "node 3 UZ load changed"); + + const auto reduced = fesa::projectToReducedSystem(assembly, dofs); + check(reduced.ok(), "reduced system projection should remain valid"); + check(reduced.k.rows() == dofs.freeDofCount(), "reduced stiffness row count changed"); + check(reduced.k.cols() == dofs.freeDofCount(), "reduced stiffness column count changed"); + check(static_cast(reduced.f.size()) == dofs.freeDofCount(), "reduced load size changed"); + check(reduced.free_full_indices == dofs.freeFullIndices(), "free full-index map changed"); + checkNear(reduced.f[0], -1.0, 1.0e-15, "first reduced load changed"); + checkNear(reduced.f[1], -1.0, 1.0e-15, "second reduced load changed"); + + fesa::DenseMatrix k(3, 3); + k(0, 0) = 4.0; + k(1, 1) = 5.0; + k(2, 2) = 6.0; + const std::vector u = {0.5, -0.25, 2.0}; + const std::vector f = {1.0, 2.0, -3.0}; + const auto rf = fesa::recoverFullReaction(k, u, f); + check(rf.size() == 3, "full reaction size changed"); + checkNear(rf[0], 1.0, 1.0e-15, "RF component 0 changed"); + checkNear(rf[1], -3.25, 1.0e-15, "RF component 1 changed"); + checkNear(rf[2], 15.0, 1.0e-15, "RF component 2 changed"); + + return 0; +}