feat: strengthen dof manager foundation

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 |

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
 

include/fesa/fesa.hpp

@@ -937,6 +937,11 @@ inline std::vector<Diagnostic> 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<LocalIndex>(free_full_indices_.size());
   }
 
+  LocalIndex constrainedDofCount() const {
+    return static_cast<LocalIndex>(constrained_full_indices_.size());
+  }
+
   const std::vector<GlobalId>& nodeIds() const {
     return node_ids_;
   }
@@ -984,6 +997,15 @@ class DofManager {
     return free_full_indices_;
   }
 
+  const std::vector<LocalIndex>& constrainedFullIndices() const {
+    return constrained_full_indices_;
+  }
+
+  DofAddress fullDof(LocalIndex full_index) const {
+    const auto& key = all_dofs_.at(static_cast<std::size_t>(full_index));
+    return {key.first, static_cast<Dof>(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<Real> reduceFullVector(const std::vector<Real>& full) const {
+    std::vector<Real> reduced;
+    reduced.reserve(free_full_indices_.size());
+    for (LocalIndex full_index : free_full_indices_) {
+      reduced.push_back(full.at(static_cast<std::size_t>(full_index)));
+    }
+    return reduced;
+  }
+
   std::vector<Real> reconstructFullVector(const std::vector<Real>& reduced) const {
     std::vector<Real> full(static_cast<std::size_t>(fullDofCount()), 0.0);
     for (std::size_t i = 0; i < free_full_indices_.size(); ++i) {
-      full[static_cast<std::size_t>(free_full_indices_[i])] = reduced[i];
+      full[static_cast<std::size_t>(free_full_indices_[i])] = reduced.at(i);
     }
     return full;
   }
 
+  std::array<LocalIndex, 24> elementFullDofIndices(const Element& element) const {
+    std::array<LocalIndex, 24> indices{};
+    for (LocalIndex node = 0; node < 4; ++node) {
+      for (Dof dof : allDofs()) {
+        const LocalIndex local = 6 * node + dofIndex(dof);
+        indices[static_cast<std::size_t>(local)] = fullIndex(element.node_ids[static_cast<std::size_t>(node)], dof);
+      }
+    }
+    return indices;
+  }
+
+  std::array<EquationId, 24> elementEquationIds(const Element& element) const {
+    std::array<EquationId, 24> equations{};
+    for (LocalIndex node = 0; node < 4; ++node) {
+      for (Dof dof : allDofs()) {
+        const LocalIndex local = 6 * node + dofIndex(dof);
+        equations[static_cast<std::size_t>(local)] = equation(element.node_ids[static_cast<std::size_t>(node)], dof);
+      }
+    }
+    return equations;
+  }
+
  private:
   std::vector<GlobalId> node_ids_;
   std::vector<std::pair<GlobalId, int>> all_dofs_;
@@ -1011,6 +1064,7 @@ class DofManager {
   std::map<std::pair<GlobalId, int>, LocalIndex> full_index_by_key_;
   std::map<std::pair<GlobalId, int>, EquationId> equation_by_key_;
   std::vector<LocalIndex> free_full_indices_;
+  std::vector<LocalIndex> constrained_full_indices_;
 };
 
 class DenseMatrix {
@@ -1056,6 +1110,18 @@ class DenseMatrix {
   std::vector<Real> values_;
 };
 
+inline std::vector<Real> recoverFullReaction(const DenseMatrix& k_full, const std::vector<Real>& u_full, const std::vector<Real>& f_full) {
+  if (k_full.rows() != k_full.cols() || static_cast<LocalIndex>(u_full.size()) != k_full.cols() ||
+      static_cast<LocalIndex>(f_full.size()) != k_full.rows()) {
+    throw std::runtime_error("full reaction size mismatch");
+  }
+  std::vector<Real> 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<Real> x;
   std::vector<Diagnostic> 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) {
+    const auto element_full_indices = dofs.elementFullDofIndices(element);
-      for (Dof da : allDofs()) {
+    for (LocalIndex a = 0; a < 24; ++a) {
-        const LocalIndex ia = dofs.fullIndex(element.node_ids[static_cast<std::size_t>(a)], da);
+      const LocalIndex ia = element_full_indices[static_cast<std::size_t>(a)];
-        const LocalIndex la = 6 * a + dofIndex(da);
+      for (LocalIndex b = 0; b < 24; ++b) {
-        for (LocalIndex b = 0; b < 4; ++b) {
+        const LocalIndex ib = element_full_indices[static_cast<std::size_t>(b)];
-          for (Dof db : allDofs()) {
+        result.k_full.add(ia, ib, ke(a, b));
-            const LocalIndex ib = dofs.fullIndex(element.node_ids[static_cast<std::size_t>(b)], db);
-            const LocalIndex lb = 6 * b + dofIndex(db);
-            result.k_full.add(ia, ib, ke(la, lb));
-          }
-        }
       }
     }
   }
@@ -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);
+    result.state.reaction_full = recoverFullReaction(assembly.k_full, result.state.u_full, result.state.f_external_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.converged = true;
     InMemoryResultsWriter writer;
     writer.writeLinearStatic(domain, dofs, result.state.u_full, result.state.reaction_full);

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" },

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<std::size_t>(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<fesa::GlobalId>({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<fesa::LocalIndex>(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<fesa::LocalIndex>({6, 7, 8, 9, 10, 11, 12, 13}));
+  FESA_CHECK(dofs.freeFullIndices() ==
+             std::vector<fesa::LocalIndex>({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<fesa::Real> full(static_cast<std::size_t>(dofs.fullDofCount()), 0.0);
+  for (std::size_t i = 0; i < full.size(); ++i) {
+    full[i] = static_cast<fesa::Real>(100 + i);
+  }
+
+  auto reduced = dofs.reduceFullVector(full);
+  FESA_CHECK(reduced.size() == static_cast<std::size_t>(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<std::size_t>(full_index)], 0.0, 1.0e-15);
+    } else {
+      FESA_CHECK_NEAR(reconstructed[static_cast<std::size_t>(full_index)], full[static_cast<std::size_t>(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<std::size_t>(i)] == i);
+  }
+  for (fesa::LocalIndex i = 0; i < 6; ++i) {
+    FESA_CHECK(equation_ids[static_cast<std::size_t>(i)] == -1);
+    FESA_CHECK(equation_ids[static_cast<std::size_t>(18 + i)] == -1);
+  }
+  FESA_CHECK(equation_ids[static_cast<std::size_t>(6 + fesa::dofIndex(fesa::Dof::UZ))] == 0);
+  FESA_CHECK(equation_ids[static_cast<std::size_t>(12 + fesa::dofIndex(fesa::Dof::UZ))] == 1);
+  FESA_CHECK(equation_ids[static_cast<std::size_t>(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<fesa::Real> u_full(static_cast<std::size_t>(dofs.fullDofCount()), 0.0);
+  std::vector<fesa::Real> f_full(static_cast<std::size_t>(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<std::size_t>(free_uz)] = 3.0;
+  f_full[static_cast<std::size_t>(support_ux)] = 1.0;
+  f_full[static_cast<std::size_t>(free_uz)] = 7.0;
+
+  auto reaction = fesa::recoverFullReaction(k_full, u_full, f_full);
+  FESA_CHECK_NEAR(reaction[static_cast<std::size_t>(support_ux)], 5.0, 1.0e-15);
+  FESA_CHECK_NEAR(reaction[static_cast<std::size_t>(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;