From b9b0963d50aab4c987ac59515f205e34f62c6184 Mon Sep 17 00:00:00 2001 From: NINI Date: Mon, 4 May 2026 13:18:28 +0900 Subject: [PATCH] feat: strengthen dof manager foundation --- PLAN.md | 8 +- PROGRESS.md | 27 +++- include/fesa/fesa.hpp | 90 ++++++++++--- .../index.json | 2 +- tests/test_main.cpp | 118 ++++++++++++++++++ 5 files changed, 223 insertions(+), 22 deletions(-) diff --git a/PLAN.md b/PLAN.md index 72d29d0..23809d2 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-05 six-DOF DofManager and reaction-foundation coverage. 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-06 results model and displacement CSV comparator foundation. 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; `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; `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, and model diagnostic context were revalidated in step 2; DofManager/results/assembly remain for steps 5, 6, and 12. | P1R-02 validation output | +| G3 - DOF/math/results infrastructure | partial | Core aliases, DOF mapping, validation harness, model diagnostic context, DofManager, sparse-connectivity inputs, and full-vector reaction formula were revalidated in steps 2 and 5; results and assembly remain for steps 6 and 12. | P1R-02 and P1R-05 validation output | | G4 - MITC4 element readiness | reopened | MITC4 formulation was rewritten from local papers; element implementation must be rebuilt or revalidated through steps 7 through 11. | Revised `docs/MITC4_FORMULATION.md`, future element tests | | 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 | @@ -86,7 +86,7 @@ All milestones are intended to become one or more self-contained sprint contract |---|---|---|---|---|---| | P1R-03 | completed | parser generator | Revalidate Phase 1 parser and immutable Domain subset. | none | Supported keywords accepted; unsupported features rejected | | P1R-04 | completed | validation generator | Rebuild validation and singular diagnostic coverage. | P1R-03 | Missing-reference and singular-prone negative tests | -| P1R-05 | pending | DOF generator | Rebuild six-DOF DofManager, constrained/free mapping, equation numbering, and full-vector reconstruction. | none | DOF mapping and reaction foundation tests | +| P1R-05 | completed | DOF generator | Rebuild six-DOF DofManager, constrained/free mapping, equation numbering, and full-vector reconstruction. | none | DOF mapping and reaction foundation tests | | P1R-06 | pending | results generator | Rebuild minimum results model and displacement CSV comparator. | none | U/RF schema tests and CSV comparator tests | | P1R-07 | pending | MITC4 generator | Implement MITC4 geometry, node order, tying points, directors, and local bases. | none | Shape/basis/diagnostic tests | | P1R-08 | pending | MITC4 generator | Implement degenerated-continuum displacement, covariant strain rows, and MITC shear tying. | P1R-07 | Finite-difference and tying interpolation tests | diff --git a/PROGRESS.md b/PROGRESS.md index a4950a7..72c0094 100644 --- a/PROGRESS.md +++ b/PROGRESS.md @@ -13,10 +13,35 @@ 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 4 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, and validation/singular diagnostics 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 5 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, and DofManager/reaction foundation 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-05 DofManager and reaction foundation 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 DofManager tests for six-DOF global ordering, noncontiguous node-id stability, constrained/free partitioning, equation numbering, full-vector reduction/reconstruction, element sparse-connectivity inputs, and a full-system reaction recovery formula. +- Extended `DofManager` with full DOF addresses, constrained full-index lists, full-vector reduction, element full-index connectivity, and element equation-id connectivity while keeping equation ids outside `Node` and `Element`. +- Added `recoverFullReaction(K_full, U_full, F_full)` and routed `LinearStaticAnalysis` through it so reaction recovery remains a full-vector operation rather than a reduced-vector shortcut. +- Updated assembly to consume element full DOF indices from `DofManager`, preserving DofManager ownership of sparse-pattern inputs. + +Verification: +- First ran `python scripts/validate_workspace.py` after adding tests; it failed as expected because the new DofManager/reaction APIs did not exist yet. +- After implementing the DofManager and reaction-foundation APIs, `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-06 results model and displacement CSV comparator foundation. +- Keep RF reference CSV availability open; current RF verification remains internal full-vector equilibrium until a stored `*_reactions.csv` artifact is provided. + ### 2026-05-04 - P1R-04 validation and singular diagnostics completed Author: Codex diff --git a/include/fesa/fesa.hpp b/include/fesa/fesa.hpp index febdc03..3357c93 100644 --- a/include/fesa/fesa.hpp +++ b/include/fesa/fesa.hpp @@ -937,6 +937,11 @@ inline std::vector validateDomain(const Domain& domain) { return diagnostics; } +struct DofAddress { + GlobalId node_id = 0; + Dof dof = Dof::UX; +}; + class DofManager { public: explicit DofManager(const Domain& domain) { @@ -951,6 +956,9 @@ class DofManager { } } for (const BoundaryCondition& boundary : domain.boundary_conditions) { + if (!validAbaqusDofRange(boundary.first_dof, boundary.last_dof)) { + continue; + } for (GlobalId node_id : resolveNodeTarget(domain, boundary.target)) { for (int dof = boundary.first_dof; dof <= boundary.last_dof; ++dof) { constrained_.insert(std::make_pair(node_id, dof - 1)); @@ -964,6 +972,7 @@ class DofManager { free_full_indices_.push_back(full_index); } else { equation_by_key_[key] = -1; + constrained_full_indices_.push_back(full_index); } } } @@ -976,6 +985,10 @@ class DofManager { return static_cast(free_full_indices_.size()); } + LocalIndex constrainedDofCount() const { + return static_cast(constrained_full_indices_.size()); + } + const std::vector& nodeIds() const { return node_ids_; } @@ -984,6 +997,15 @@ class DofManager { return free_full_indices_; } + const std::vector& constrainedFullIndices() const { + return constrained_full_indices_; + } + + DofAddress fullDof(LocalIndex full_index) const { + const auto& key = all_dofs_.at(static_cast(full_index)); + return {key.first, static_cast(key.second)}; + } + LocalIndex fullIndex(GlobalId node_id, Dof dof) const { return full_index_by_key_.at(std::make_pair(node_id, dofIndex(dof))); } @@ -996,14 +1018,45 @@ class DofManager { return constrained_.count(std::make_pair(node_id, dofIndex(dof))) != 0; } + std::vector reduceFullVector(const std::vector& full) const { + std::vector reduced; + reduced.reserve(free_full_indices_.size()); + for (LocalIndex full_index : free_full_indices_) { + reduced.push_back(full.at(static_cast(full_index))); + } + return reduced; + } + std::vector reconstructFullVector(const std::vector& reduced) const { std::vector full(static_cast(fullDofCount()), 0.0); for (std::size_t i = 0; i < free_full_indices_.size(); ++i) { - full[static_cast(free_full_indices_[i])] = reduced[i]; + full[static_cast(free_full_indices_[i])] = reduced.at(i); } return full; } + std::array elementFullDofIndices(const Element& element) const { + std::array indices{}; + for (LocalIndex node = 0; node < 4; ++node) { + for (Dof dof : allDofs()) { + const LocalIndex local = 6 * node + dofIndex(dof); + indices[static_cast(local)] = fullIndex(element.node_ids[static_cast(node)], dof); + } + } + return indices; + } + + std::array elementEquationIds(const Element& element) const { + std::array equations{}; + for (LocalIndex node = 0; node < 4; ++node) { + for (Dof dof : allDofs()) { + const LocalIndex local = 6 * node + dofIndex(dof); + equations[static_cast(local)] = equation(element.node_ids[static_cast(node)], dof); + } + } + return equations; + } + private: std::vector node_ids_; std::vector> all_dofs_; @@ -1011,6 +1064,7 @@ class DofManager { std::map, LocalIndex> full_index_by_key_; std::map, EquationId> equation_by_key_; std::vector free_full_indices_; + std::vector constrained_full_indices_; }; class DenseMatrix { @@ -1056,6 +1110,18 @@ class DenseMatrix { 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()) { + 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 SolveResult { std::vector x; std::vector diagnostics; @@ -1378,17 +1444,12 @@ inline AssemblyResult assembleSystem(const Domain& domain, const DofManager& dof 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); - for (LocalIndex a = 0; a < 4; ++a) { - for (Dof da : allDofs()) { - const LocalIndex ia = dofs.fullIndex(element.node_ids[static_cast(a)], da); - const LocalIndex la = 6 * a + dofIndex(da); - for (LocalIndex b = 0; b < 4; ++b) { - for (Dof db : allDofs()) { - const LocalIndex ib = dofs.fullIndex(element.node_ids[static_cast(b)], db); - const LocalIndex lb = 6 * b + dofIndex(db); - result.k_full.add(ia, ib, ke(la, lb)); - } - } + 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)); } } } @@ -1546,10 +1607,7 @@ class LinearStaticAnalysis final : public Analysis { } result.state.u_full = dofs.reconstructFullVector(solved.x); result.state.f_external_full = assembly.f_full; - result.state.reaction_full = assembly.k_full.multiply(result.state.u_full); - for (std::size_t i = 0; i < result.state.reaction_full.size(); ++i) { - result.state.reaction_full[i] -= result.state.f_external_full[i]; - } + 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); diff --git a/phases/1-linear-static-mitc4-rebaseline/index.json b/phases/1-linear-static-mitc4-rebaseline/index.json index 7e2ce13..ff7ed83 100644 --- a/phases/1-linear-static-mitc4-rebaseline/index.json +++ b/phases/1-linear-static-mitc4-rebaseline/index.json @@ -7,7 +7,7 @@ { "step": 2, "name": "core-harness-guardrails", "status": "completed" }, { "step": 3, "name": "parser-domain-subset", "status": "completed" }, { "step": 4, "name": "validation-singular-diagnostics", "status": "completed" }, - { "step": 5, "name": "dof-manager-reaction-foundation", "status": "pending" }, + { "step": 5, "name": "dof-manager-reaction-foundation", "status": "completed" }, { "step": 6, "name": "results-comparator-foundation", "status": "pending" }, { "step": 7, "name": "mitc4-geometry-directors", "status": "pending" }, { "step": 8, "name": "mitc4-covariant-strain-tying", "status": "pending" }, diff --git a/tests/test_main.cpp b/tests/test_main.cpp index e587e5f..dc32a24 100644 --- a/tests/test_main.cpp +++ b/tests/test_main.cpp @@ -139,6 +139,22 @@ fesa::Domain singleElementValidationDomain() { return domain; } +fesa::Domain noncontiguousDofDomain() { + fesa::Domain domain; + domain.nodes[30] = {30, {1, 1, 0}}; + domain.nodes[10] = {10, {0, 0, 0}}; + domain.nodes[40] = {40, {0, 1, 0}}; + domain.nodes[20] = {20, {1, 0, 0}}; + domain.elements[5] = {5, fesa::ElementType::MITC4, {10, 20, 30, 40}, "EALL"}; + domain.element_sets["eall"] = {"EALL", {5}}; + domain.node_sets["clamped"] = {"CLAMPED", {20}}; + domain.materials["mat"] = {"MAT", 1000.0, 0.3}; + domain.shell_sections.push_back({"EALL", "MAT", 0.1}); + domain.boundary_conditions.push_back({"CLAMPED", 1, 6, 0.0}); + domain.boundary_conditions.push_back({"30", 1, 2, 0.0}); + return domain; +} + } // namespace FESA_TEST(core_types_and_dof_mapping_are_stable) { @@ -501,6 +517,108 @@ FESA_TEST(dof_manager_owns_equation_numbering_and_reconstruction) { FESA_CHECK_NEAR(full[static_cast(dofs.fullIndex(1, fesa::Dof::UX))], 0.0, 1.0e-15); } +FESA_TEST(dof_manager_preserves_global_order_for_noncontiguous_nodes) { + auto domain = noncontiguousDofDomain(); + fesa::DofManager dofs(domain); + + FESA_CHECK(dofs.nodeIds() == std::vector({10, 20, 30, 40})); + FESA_CHECK(dofs.fullDofCount() == 24); + for (std::size_t node_offset = 0; node_offset < dofs.nodeIds().size(); ++node_offset) { + const fesa::GlobalId node_id = dofs.nodeIds()[node_offset]; + for (fesa::Dof dof : fesa::allDofs()) { + const fesa::LocalIndex expected = static_cast(6 * node_offset + fesa::dofIndex(dof)); + FESA_CHECK(dofs.fullIndex(node_id, dof) == expected); + const auto address = dofs.fullDof(expected); + FESA_CHECK(address.node_id == node_id); + FESA_CHECK(address.dof == dof); + } + } +} + +FESA_TEST(dof_manager_partitions_constraints_and_equations_deterministically) { + auto domain = noncontiguousDofDomain(); + fesa::DofManager dofs(domain); + + FESA_CHECK(dofs.constrainedDofCount() == 8); + FESA_CHECK(dofs.freeDofCount() == 16); + FESA_CHECK(dofs.constrainedFullIndices() == std::vector({6, 7, 8, 9, 10, 11, 12, 13})); + FESA_CHECK(dofs.freeFullIndices() == + std::vector({0, 1, 2, 3, 4, 5, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23})); + FESA_CHECK(dofs.equation(10, fesa::Dof::UX) == 0); + FESA_CHECK(dofs.equation(20, fesa::Dof::RZ) == -1); + FESA_CHECK(dofs.equation(30, fesa::Dof::UX) == -1); + FESA_CHECK(dofs.equation(30, fesa::Dof::UZ) == 6); + FESA_CHECK(dofs.equation(40, fesa::Dof::RZ) == 15); +} + +FESA_TEST(dof_manager_reduces_and_reconstructs_full_vectors) { + auto domain = noncontiguousDofDomain(); + fesa::DofManager dofs(domain); + + std::vector full(static_cast(dofs.fullDofCount()), 0.0); + for (std::size_t i = 0; i < full.size(); ++i) { + full[i] = static_cast(100 + i); + } + + auto reduced = dofs.reduceFullVector(full); + FESA_CHECK(reduced.size() == static_cast(dofs.freeDofCount())); + FESA_CHECK_NEAR(reduced[0], full[0], 1.0e-15); + FESA_CHECK_NEAR(reduced[6], full[14], 1.0e-15); + FESA_CHECK_NEAR(reduced.back(), full[23], 1.0e-15); + + auto reconstructed = dofs.reconstructFullVector(reduced); + for (fesa::LocalIndex full_index = 0; full_index < dofs.fullDofCount(); ++full_index) { + const auto address = dofs.fullDof(full_index); + if (dofs.isConstrained(address.node_id, address.dof)) { + FESA_CHECK_NEAR(reconstructed[static_cast(full_index)], 0.0, 1.0e-15); + } else { + FESA_CHECK_NEAR(reconstructed[static_cast(full_index)], full[static_cast(full_index)], 1.0e-15); + } + } +} + +FESA_TEST(dof_manager_provides_element_sparse_connectivity_inputs) { + auto domain = parsedPhase1Domain(); + fesa::DofManager dofs(domain); + const auto& element = domain.elements.at(1); + + auto full_indices = dofs.elementFullDofIndices(element); + auto equation_ids = dofs.elementEquationIds(element); + + for (fesa::LocalIndex i = 0; i < 24; ++i) { + FESA_CHECK(full_indices[static_cast(i)] == i); + } + for (fesa::LocalIndex i = 0; i < 6; ++i) { + FESA_CHECK(equation_ids[static_cast(i)] == -1); + FESA_CHECK(equation_ids[static_cast(18 + i)] == -1); + } + FESA_CHECK(equation_ids[static_cast(6 + fesa::dofIndex(fesa::Dof::UZ))] == 0); + FESA_CHECK(equation_ids[static_cast(12 + fesa::dofIndex(fesa::Dof::UZ))] == 1); + FESA_CHECK(equation_ids[static_cast(6 + fesa::dofIndex(fesa::Dof::UX))] == -1); +} + +FESA_TEST(full_vector_reaction_recovery_uses_full_system_quantities) { + auto domain = noncontiguousDofDomain(); + fesa::DofManager dofs(domain); + fesa::DenseMatrix k_full(dofs.fullDofCount(), dofs.fullDofCount()); + std::vector u_full(static_cast(dofs.fullDofCount()), 0.0); + std::vector f_full(static_cast(dofs.fullDofCount()), 0.0); + + const fesa::LocalIndex support_ux = dofs.fullIndex(20, fesa::Dof::UX); + const fesa::LocalIndex free_uz = dofs.fullIndex(30, fesa::Dof::UZ); + k_full(support_ux, support_ux) = 10.0; + k_full(support_ux, free_uz) = 2.0; + k_full(free_uz, support_ux) = 2.0; + k_full(free_uz, free_uz) = 5.0; + u_full[static_cast(free_uz)] = 3.0; + f_full[static_cast(support_ux)] = 1.0; + f_full[static_cast(free_uz)] = 7.0; + + auto reaction = fesa::recoverFullReaction(k_full, u_full, f_full); + FESA_CHECK_NEAR(reaction[static_cast(support_ux)], 5.0, 1.0e-15); + FESA_CHECK_NEAR(reaction[static_cast(free_uz)], 8.0, 1.0e-15); +} + FESA_TEST(gaussian_solver_solves_and_diagnoses_singular_systems) { fesa::DenseMatrix a(2, 2); a(0, 0) = 2.0;