diff --git a/CMakeLists.txt b/CMakeLists.txt index 44d0ea3..50dc7a3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -21,15 +21,21 @@ target_compile_definitions(fesa_tests PRIVATE FESA_SOURCE_DIR="${CMAKE_CURRENT_S add_executable(fesa_core_module_tests tests/test_core_module_includes.cpp) target_link_libraries(fesa_core_module_tests PRIVATE fesa_core) +add_executable(fesa_math_module_tests tests/test_math_module_includes.cpp) +target_link_libraries(fesa_math_module_tests PRIVATE fesa_core) + if(MSVC) target_compile_options(fesa_core PRIVATE /W4 /permissive-) target_compile_options(fesa_tests PRIVATE /W4 /permissive-) target_compile_options(fesa_core_module_tests PRIVATE /W4 /permissive-) + target_compile_options(fesa_math_module_tests PRIVATE /W4 /permissive-) else() target_compile_options(fesa_core PRIVATE -Wall -Wextra -Wpedantic) target_compile_options(fesa_tests PRIVATE -Wall -Wextra -Wpedantic) target_compile_options(fesa_core_module_tests PRIVATE -Wall -Wextra -Wpedantic) + target_compile_options(fesa_math_module_tests PRIVATE -Wall -Wextra -Wpedantic) endif() add_test(NAME fesa_tests COMMAND fesa_tests) add_test(NAME fesa_core_module_tests COMMAND fesa_core_module_tests) +add_test(NAME fesa_math_module_tests COMMAND fesa_math_module_tests) diff --git a/PLAN.md b/PLAN.md index 49e7c02..d72d168 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-03 Math and solver adapter extraction. P1A-00 completed the architecture drift audit, P1A-01 created the module scaffold, and P1A-02 extracted Core/Util plus Phase 1 Boundary/Load/Property model ownership 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-04 IO parser 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, and P1A-03 extracted Math primitives plus the solver adapter boundary 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. ## 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` extracts Math and solver adapter boundaries; `step4.md` extracts the Abaqus parser into IO; `step5.md` extracts Results and reference comparison code; `step6.md` extracts MITC4 geometry/strain helpers; `step7.md` extracts 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` extracts the Abaqus parser into IO; `step5.md` extracts Results and reference comparison code; `step6.md` extracts MITC4 geometry/strain helpers; `step7.md` extracts MITC4 material/stiffness helpers; `step8.md` extracts 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`. @@ -60,7 +60,7 @@ This phase is an architecture-preserving refactor. It must not change Phase 1 so | P1A-00 | completed | planner/evaluator | Audit `fesa.hpp` architecture drift and create a symbol-to-module migration map. | P1R-15 | Complete migration map and validation baseline | | P1A-01 | completed | generator | Create module directory scaffold, CMake source boundaries, and umbrella facade policy. | P1A-00 | Module include smoke tests and build stability | | P1A-02 | completed | generator | Extract Core/Util domain, diagnostics, aliases, DOF mapping, `AnalysisModel`, `DofManager`, and Phase 1 Boundary/Load/Property model ownership. | P1A-01 | Core has no dependency on higher layers; Boundary/Load/Property types are no longer hidden in the umbrella header; DOF tests unchanged | -| P1A-03 | pending | generator | Extract Math and solver adapter boundaries. | P1A-02 | Linear solver interface remains adapter-ready; int64 paths unchanged | +| P1A-03 | completed | generator | Extract Math and solver adapter boundaries. | P1A-02 | Linear solver interface remains adapter-ready; int64 paths unchanged | | P1A-04 | pending | generator | Extract Abaqus parser into IO. | P1A-02 | Parser subset and unsupported-feature diagnostics unchanged | | P1A-05 | pending | 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 | pending | generator | Extract MITC4 geometry, director, strain, and tying helpers into Element. | P1A-03 | Geometry/strain tests and formulation signs unchanged | diff --git a/PROGRESS.md b/PROGRESS.md index 4cbc594..20908dc 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 the current production implementation is concentrated in `include/fesa/fesa.hpp` instead of the module directories documented in `docs/ARCHITECTURE.md`; P1A-00, P1A-01, and P1A-02 are complete, so the next step is P1A-03 Math and solver adapter 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 because the current production implementation is concentrated in `include/fesa/fesa.hpp` instead of the module directories documented in `docs/ARCHITECTURE.md`; P1A-00 through P1A-03 are complete, so the next step is P1A-04 IO parser 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. ## Completed Work +### 2026-05-05 - P1A-03 Math solver extraction completed +Author: Codex + +Changed files: +- `CMakeLists.txt` +- `include/fesa/Core/Domain.hpp` +- `include/fesa/Math/DenseMatrix.hpp` +- `include/fesa/Math/LinearSolver.hpp` +- `include/fesa/Math/Math.hpp` +- `include/fesa/Math/SparsePattern.hpp` +- `include/fesa/Math/Vector.hpp` +- `include/fesa/fesa.hpp` +- `tests/test_math_module_includes.cpp` +- `phases/1-structure-alignment-refactor/index.json` +- `PLAN.md` +- `PROGRESS.md` + +Summary: +- Extracted `Vec3` and vector helper functions into `include/fesa/Math/Vector.hpp`; `Core/Domain.hpp` now consumes `Vec3` from Math instead of defining it locally. +- Extracted sparse pattern data structures into `include/fesa/Math/SparsePattern.hpp`. +- Extracted `DenseMatrix` into `include/fesa/Math/DenseMatrix.hpp`. +- Extracted `SolveResult`, `LinearSolver`, and `GaussianEliminationSolver` into `include/fesa/Math/LinearSolver.hpp`. +- Added `fesa_math_module_tests`, a direct module include smoke test that does not include `fesa/fesa.hpp` and checks vector math, int64 sparse index/equation boundaries, dense matrix multiply, solver success, and singular-solver diagnostics. +- Preserved `LinearSolver` as the adapter boundary and introduced no MKL, TBB, HDF5, or production sparse storage dependency. +- Remaining large groups in `fesa.hpp` are IO parser, Assembly helpers (`buildReducedSparsePattern`, `recoverFullReaction`), MITC4 Element/Material helpers, Results/reference comparison, and Analysis workflow. + +Verification: +- First ran `python scripts/validate_workspace.py` after adding the direct Math include test; it failed as expected because `fesa/Math/Math.hpp` did not yet expose Math primitives or solver types. +- After extraction, `python scripts/validate_workspace.py` configured CMake, built `fesa_core`, `fesa_tests`, `fesa_core_module_tests`, and `fesa_math_module_tests`, and ran CTest successfully. +- CTest result: 3 test executables passed. + +Follow-up: +- Continue with P1A-04 IO parser extraction. +- Keep R-014 open until P1A-09 independently accepts the final architecture alignment. + ### 2026-05-05 - P1A-02 Core domain DOF extraction completed Author: Codex diff --git a/include/fesa/Core/Domain.hpp b/include/fesa/Core/Domain.hpp index af1f17b..62e6ea0 100644 --- a/include/fesa/Core/Domain.hpp +++ b/include/fesa/Core/Domain.hpp @@ -3,6 +3,7 @@ #include "fesa/Boundary/Boundary.hpp" #include "fesa/Core/Types.hpp" #include "fesa/Load/Load.hpp" +#include "fesa/Math/Vector.hpp" #include "fesa/Property/Property.hpp" #include "fesa/Util/String.hpp" @@ -13,12 +14,6 @@ namespace fesa { -struct Vec3 { - Real x = 0.0; - Real y = 0.0; - Real z = 0.0; -}; - struct Node { GlobalId id = 0; Vec3 coordinates; diff --git a/include/fesa/Math/DenseMatrix.hpp b/include/fesa/Math/DenseMatrix.hpp new file mode 100644 index 0000000..56da9f2 --- /dev/null +++ b/include/fesa/Math/DenseMatrix.hpp @@ -0,0 +1,53 @@ +#pragma once + +#include "fesa/Core/Types.hpp" + +#include +#include + +namespace fesa { + +class DenseMatrix { + public: + DenseMatrix() = default; + DenseMatrix(LocalIndex rows, LocalIndex cols) : rows_(rows), cols_(cols), values_(static_cast(rows * cols), 0.0) {} + + LocalIndex rows() const { + return rows_; + } + + LocalIndex cols() const { + return cols_; + } + + Real& operator()(LocalIndex row, LocalIndex col) { + return values_[static_cast(row * cols_ + col)]; + } + + Real operator()(LocalIndex row, LocalIndex col) const { + return values_[static_cast(row * cols_ + col)]; + } + + void add(LocalIndex row, LocalIndex col, Real value) { + (*this)(row, col) += value; + } + + std::vector multiply(const std::vector& x) const { + std::vector y(static_cast(rows_), 0.0); + for (LocalIndex i = 0; i < rows_; ++i) { + Real sum = 0.0; + for (LocalIndex j = 0; j < cols_; ++j) { + sum += (*this)(i, j) * x[static_cast(j)]; + } + y[static_cast(i)] = sum; + } + return y; + } + + private: + LocalIndex rows_ = 0; + LocalIndex cols_ = 0; + std::vector values_; +}; + +} // namespace fesa diff --git a/include/fesa/Math/LinearSolver.hpp b/include/fesa/Math/LinearSolver.hpp new file mode 100644 index 0000000..63b0204 --- /dev/null +++ b/include/fesa/Math/LinearSolver.hpp @@ -0,0 +1,81 @@ +#pragma once + +#include "fesa/Math/DenseMatrix.hpp" +#include "fesa/Util/Diagnostics.hpp" + +#include +#include +#include +#include +#include + +namespace fesa { + +struct SolveResult { + std::vector x; + std::vector diagnostics; + + bool ok() const { + return !hasError(diagnostics); + } +}; + +class LinearSolver { + public: + virtual ~LinearSolver() = default; + virtual SolveResult solve(DenseMatrix a, std::vector b) const = 0; +}; + +class GaussianEliminationSolver final : public LinearSolver { + public: + SolveResult solve(DenseMatrix a, std::vector b) const override { + const LocalIndex n = a.rows(); + SolveResult result; + if (a.rows() != a.cols() || static_cast(b.size()) != n) { + result.diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-SOLVER-SIZE", "Linear system size mismatch", "solver")); + return result; + } + for (LocalIndex col = 0; col < n; ++col) { + LocalIndex pivot = col; + Real pivot_abs = std::fabs(a(col, col)); + for (LocalIndex row = col + 1; row < n; ++row) { + const Real candidate = std::fabs(a(row, col)); + if (candidate > pivot_abs) { + pivot_abs = candidate; + pivot = row; + } + } + if (pivot_abs < 1.0e-12) { + result.diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-SINGULAR-SOLVER", + "Reduced system is singular or ill-conditioned", "solver")); + return result; + } + if (pivot != col) { + for (LocalIndex j = col; j < n; ++j) { + std::swap(a(col, j), a(pivot, j)); + } + std::swap(b[static_cast(col)], b[static_cast(pivot)]); + } + const Real diag = a(col, col); + for (LocalIndex row = col + 1; row < n; ++row) { + const Real factor = a(row, col) / diag; + a(row, col) = 0.0; + for (LocalIndex j = col + 1; j < n; ++j) { + a(row, j) -= factor * a(col, j); + } + b[static_cast(row)] -= factor * b[static_cast(col)]; + } + } + result.x.assign(static_cast(n), 0.0); + for (LocalIndex i = n; i-- > 0;) { + Real sum = b[static_cast(i)]; + for (LocalIndex j = i + 1; j < n; ++j) { + sum -= a(i, j) * result.x[static_cast(j)]; + } + result.x[static_cast(i)] = sum / a(i, i); + } + return result; + } +}; + +} // namespace fesa diff --git a/include/fesa/Math/Math.hpp b/include/fesa/Math/Math.hpp index 6c6c6f0..5e2c25e 100644 --- a/include/fesa/Math/Math.hpp +++ b/include/fesa/Math/Math.hpp @@ -1,5 +1,9 @@ #pragma once +#include "fesa/Math/DenseMatrix.hpp" +#include "fesa/Math/LinearSolver.hpp" +#include "fesa/Math/SparsePattern.hpp" +#include "fesa/Math/Vector.hpp" #include "fesa/ModuleInfo.hpp" namespace fesa::module { diff --git a/include/fesa/Math/SparsePattern.hpp b/include/fesa/Math/SparsePattern.hpp new file mode 100644 index 0000000..c2c49f1 --- /dev/null +++ b/include/fesa/Math/SparsePattern.hpp @@ -0,0 +1,30 @@ +#pragma once + +#include "fesa/Core/Types.hpp" + +#include +#include + +namespace fesa { + +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; + }); + } +}; + +} // namespace fesa diff --git a/include/fesa/Math/Vector.hpp b/include/fesa/Math/Vector.hpp new file mode 100644 index 0000000..b4be202 --- /dev/null +++ b/include/fesa/Math/Vector.hpp @@ -0,0 +1,66 @@ +#pragma once + +#include "fesa/Core/Types.hpp" + +#include +#include +#include +#include + +namespace fesa { + +struct Vec3 { + Real x = 0.0; + Real y = 0.0; + Real z = 0.0; +}; + +inline Vec3 operator+(const Vec3& a, const Vec3& b) { + return {a.x + b.x, a.y + b.y, a.z + b.z}; +} + +inline Vec3 operator-(const Vec3& a, const Vec3& b) { + return {a.x - b.x, a.y - b.y, a.z - b.z}; +} + +inline Vec3 operator*(Real scalar, const Vec3& value) { + return {scalar * value.x, scalar * value.y, scalar * value.z}; +} + +inline Real dot(const Vec3& a, const Vec3& b) { + return a.x * b.x + a.y * b.y + a.z * b.z; +} + +inline Vec3 cross(const Vec3& a, const Vec3& b) { + return {a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x}; +} + +inline Real norm(const Vec3& value) { + return std::sqrt(dot(value, value)); +} + +inline bool isFinite(Real value) { + return std::isfinite(value); +} + +inline bool isFinite(const Vec3& value) { + return isFinite(value.x) && isFinite(value.y) && isFinite(value.z); +} + +inline std::optional normalizedIfValid(const Vec3& value, Real tolerance = 1.0e-12) { + const Real length = norm(value); + if (!isFinite(length) || length <= tolerance) { + return std::nullopt; + } + return (1.0 / length) * value; +} + +inline Vec3 normalized(const Vec3& value) { + const Real length = norm(value); + if (length <= std::numeric_limits::epsilon()) { + throw std::runtime_error("zero-length vector"); + } + return (1.0 / length) * value; +} + +} // namespace fesa diff --git a/include/fesa/fesa.hpp b/include/fesa/fesa.hpp index 6dff244..caabaf7 100644 --- a/include/fesa/fesa.hpp +++ b/include/fesa/fesa.hpp @@ -3,6 +3,7 @@ #include "fesa/Boundary/Boundary.hpp" #include "fesa/Core/Core.hpp" #include "fesa/Load/Load.hpp" +#include "fesa/Math/Math.hpp" #include "fesa/ModuleInfo.hpp" #include "fesa/Property/Property.hpp" #include "fesa/Util/Util.hpp" @@ -27,54 +28,6 @@ namespace fesa { -inline Vec3 operator+(const Vec3& a, const Vec3& b) { - return {a.x + b.x, a.y + b.y, a.z + b.z}; -} - -inline Vec3 operator-(const Vec3& a, const Vec3& b) { - return {a.x - b.x, a.y - b.y, a.z - b.z}; -} - -inline Vec3 operator*(Real scalar, const Vec3& value) { - return {scalar * value.x, scalar * value.y, scalar * value.z}; -} - -inline Real dot(const Vec3& a, const Vec3& b) { - return a.x * b.x + a.y * b.y + a.z * b.z; -} - -inline Vec3 cross(const Vec3& a, const Vec3& b) { - return {a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x}; -} - -inline Real norm(const Vec3& value) { - return std::sqrt(dot(value, value)); -} - -inline bool isFinite(Real value) { - return std::isfinite(value); -} - -inline bool isFinite(const Vec3& value) { - return isFinite(value.x) && isFinite(value.y) && isFinite(value.z); -} - -inline std::optional normalizedIfValid(const Vec3& value, Real tolerance = 1.0e-12) { - const Real length = norm(value); - if (!isFinite(length) || length <= tolerance) { - return std::nullopt; - } - return (1.0 / length) * value; -} - -inline Vec3 normalized(const Vec3& value) { - const Real length = norm(value); - if (length <= std::numeric_limits::epsilon()) { - throw std::runtime_error("zero-length vector"); - } - return (1.0 / length) * value; -} - struct KeywordLine { std::string name; std::map parameters; @@ -518,26 +471,6 @@ class AbaqusInputParser { } }; -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(); @@ -564,49 +497,6 @@ inline SparsePattern buildReducedSparsePattern(const Domain& domain, const DofMa return pattern; } -class DenseMatrix { - public: - DenseMatrix() = default; - DenseMatrix(LocalIndex rows, LocalIndex cols) : rows_(rows), cols_(cols), values_(static_cast(rows * cols), 0.0) {} - - LocalIndex rows() const { - return rows_; - } - - LocalIndex cols() const { - return cols_; - } - - Real& operator()(LocalIndex row, LocalIndex col) { - return values_[static_cast(row * cols_ + col)]; - } - - Real operator()(LocalIndex row, LocalIndex col) const { - return values_[static_cast(row * cols_ + col)]; - } - - void add(LocalIndex row, LocalIndex col, Real value) { - (*this)(row, col) += value; - } - - std::vector multiply(const std::vector& x) const { - std::vector y(static_cast(rows_), 0.0); - for (LocalIndex i = 0; i < rows_; ++i) { - Real sum = 0.0; - for (LocalIndex j = 0; j < cols_; ++j) { - sum += (*this)(i, j) * x[static_cast(j)]; - } - y[static_cast(i)] = sum; - } - return y; - } - - private: - LocalIndex rows_ = 0; - LocalIndex cols_ = 0; - std::vector values_; -}; - 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()) { @@ -619,73 +509,6 @@ inline std::vector recoverFullReaction(const DenseMatrix& k_full, const st return reaction; } -struct SolveResult { - std::vector x; - std::vector diagnostics; - - bool ok() const { - return !hasError(diagnostics); - } -}; - -class LinearSolver { - public: - virtual ~LinearSolver() = default; - virtual SolveResult solve(DenseMatrix a, std::vector b) const = 0; -}; - -class GaussianEliminationSolver final : public LinearSolver { - public: - SolveResult solve(DenseMatrix a, std::vector b) const override { - const LocalIndex n = a.rows(); - SolveResult result; - if (a.rows() != a.cols() || static_cast(b.size()) != n) { - result.diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-SOLVER-SIZE", "Linear system size mismatch", "solver")); - return result; - } - for (LocalIndex col = 0; col < n; ++col) { - LocalIndex pivot = col; - Real pivot_abs = std::fabs(a(col, col)); - for (LocalIndex row = col + 1; row < n; ++row) { - const Real candidate = std::fabs(a(row, col)); - if (candidate > pivot_abs) { - pivot_abs = candidate; - pivot = row; - } - } - if (pivot_abs < 1.0e-12) { - result.diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-SINGULAR-SOLVER", - "Reduced system is singular or ill-conditioned", "solver")); - return result; - } - if (pivot != col) { - for (LocalIndex j = col; j < n; ++j) { - std::swap(a(col, j), a(pivot, j)); - } - std::swap(b[static_cast(col)], b[static_cast(pivot)]); - } - const Real diag = a(col, col); - for (LocalIndex row = col + 1; row < n; ++row) { - const Real factor = a(row, col) / diag; - a(row, col) = 0.0; - for (LocalIndex j = col + 1; j < n; ++j) { - a(row, j) -= factor * a(col, j); - } - b[static_cast(row)] -= factor * b[static_cast(col)]; - } - } - result.x.assign(static_cast(n), 0.0); - for (LocalIndex i = n; i-- > 0;) { - Real sum = b[static_cast(i)]; - for (LocalIndex j = i + 1; j < n; ++j) { - sum -= a(i, j) * result.x[static_cast(j)]; - } - result.x[static_cast(i)] = sum / a(i, i); - } - return result; - } -}; - struct ShapeData { std::array n{}; std::array dr{}; diff --git a/phases/1-structure-alignment-refactor/index.json b/phases/1-structure-alignment-refactor/index.json index 2bbafad..4632d30 100644 --- a/phases/1-structure-alignment-refactor/index.json +++ b/phases/1-structure-alignment-refactor/index.json @@ -5,7 +5,7 @@ { "step": 0, "name": "architecture-drift-audit", "status": "completed", "artifact": "step0-architecture-map.md" }, { "step": 1, "name": "module-scaffold-and-facade", "status": "completed" }, { "step": 2, "name": "core-domain-dof-extraction", "status": "completed" }, - { "step": 3, "name": "math-solver-extraction", "status": "pending" }, + { "step": 3, "name": "math-solver-extraction", "status": "completed" }, { "step": 4, "name": "io-parser-extraction", "status": "pending" }, { "step": 5, "name": "results-reference-extraction", "status": "pending" }, { "step": 6, "name": "mitc4-geometry-strain-extraction", "status": "pending" }, diff --git a/tests/test_math_module_includes.cpp b/tests/test_math_module_includes.cpp new file mode 100644 index 0000000..307e80c --- /dev/null +++ b/tests/test_math_module_includes.cpp @@ -0,0 +1,71 @@ +#include "fesa/Math/Math.hpp" + +#include +#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); + } +} + +} // namespace + +int main() { + static_assert(std::is_same_v, "SparseIndex must remain int64"); + static_assert(std::is_same_v, "EquationId must remain int64"); + + const fesa::Vec3 x{1.0, 0.0, 0.0}; + const fesa::Vec3 y{0.0, 1.0, 0.0}; + const auto z = fesa::cross(x, y); + checkNear(fesa::dot(x, y), 0.0, 1.0e-12, "Vector dot product changed"); + checkNear(z.z, 1.0, 1.0e-12, "Vector cross product changed"); + checkNear(fesa::norm(fesa::Vec3{3.0, 4.0, 0.0}), 5.0, 1.0e-12, "Vector norm changed"); + check(fesa::normalizedIfValid(fesa::Vec3{0.0, 0.0, 0.0}) == std::nullopt, "Invalid normalization changed"); + + fesa::SparsePattern pattern; + pattern.equation_count = 3; + pattern.entries.push_back({0, 0}); + pattern.entries.push_back({0, 2}); + check(pattern.nonzeroCount() == 2, "Sparse nonzero count changed"); + check(pattern.contains(0, 2), "Sparse pattern lookup changed"); + check(!pattern.contains(2, 0), "Sparse pattern lookup became nondirectional"); + + fesa::DenseMatrix matrix(2, 2); + matrix(0, 0) = 2.0; + matrix(0, 1) = 1.0; + matrix(1, 0) = 1.0; + matrix(1, 1) = 3.0; + const auto product = matrix.multiply({1.0, 2.0}); + checkNear(product[0], 4.0, 1.0e-12, "DenseMatrix multiply row 0 changed"); + checkNear(product[1], 7.0, 1.0e-12, "DenseMatrix multiply row 1 changed"); + + const fesa::LinearSolver& solver = fesa::GaussianEliminationSolver{}; + const auto solved = solver.solve(matrix, {1.0, 2.0}); + check(solved.ok(), "Gaussian solver unexpectedly failed"); + checkNear(solved.x[0], 0.2, 1.0e-12, "Gaussian solver x0 changed"); + checkNear(solved.x[1], 0.6, 1.0e-12, "Gaussian solver x1 changed"); + + fesa::DenseMatrix singular(2, 2); + singular(0, 0) = 1.0; + singular(0, 1) = 2.0; + singular(1, 0) = 2.0; + singular(1, 1) = 4.0; + const auto failed = solver.solve(singular, {1.0, 2.0}); + check(!failed.ok(), "Singular solve unexpectedly passed"); + check(fesa::containsDiagnostic(failed.diagnostics, "FESA-SINGULAR-SOLVER"), "Singular diagnostic code changed"); + + return 0; +}