docs: close euler beam kernel phase

CMakeLists.txt

@@ -0,0 +1,12 @@
+cmake_minimum_required(VERSION 3.20)
+
+project(FESA LANGUAGES CXX)
+
+file(GLOB_RECURSE FESA_CORE_SOURCES CONFIGURE_DEPENDS src/fesa/*.cpp)
+
+add_library(fesa_core STATIC ${FESA_CORE_SOURCES})
+target_compile_features(fesa_core PUBLIC cxx_std_17)
+target_include_directories(fesa_core PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/src)
+
+enable_testing()
+add_subdirectory(tests)

docs/build-test-reports/euler-beam-3d-build-test.md

@@ -0,0 +1,75 @@
+# 3D Euler-Bernoulli Beam Kernel Build/Test Report
+
+## Metadata
+- feature_id: euler-beam-3d
+- source_implementation_plan: docs/implementation-plans/euler-beam-3d-implementation-plan.md
+- status: pass-for-reference-verification
+- owner_agent: build-test-executor-agent
+- date: 2026-06-12
+
+## Changed Files Observed
+- `docs/research/euler-beam-3d-research.md`
+- `docs/formulations/euler-beam-3d-formulation.md`
+- `docs/numerical-reviews/euler-beam-3d-review.md`
+- `docs/io-definitions/euler-beam-3d-io.md`
+- `docs/reference-models/euler-beam-3d-reference-models.md`
+- `docs/implementation-plans/euler-beam-3d-implementation-plan.md`
+- `src/fesa/model/element.hpp`
+- `tests/unit/model_element_test.cpp`
+- `src/fesa/elements/euler_beam_3d.hpp`
+- `src/fesa/elements/euler_beam_3d.cpp`
+- `tests/unit/euler_beam_3d_local_stiffness_test.cpp`
+- `tests/unit/euler_beam_3d_transform_recovery_test.cpp`
+- `phases/euler-beam-3d/index.json`
+
+## Execution Environment
+- os: Windows, MSVC toolchain through Visual Studio 17 2022
+- generator: Visual Studio 17 2022
+- platform: x64
+- config: Debug
+- build_dir: build/msvc-debug
+- active_override_env_vars: `HARNESS_VALIDATION_COMMANDS`
+- command_discovery_path: `HARNESS_VALIDATION_COMMANDS` with absolute CMake and CTest executable paths
+
+## Command Log Summary
+
+| order | command | exit_code | evidence_tail |
+| --- | --- | --- | --- |
+| 1 | `ctest --test-dir build/msvc-debug --output-on-failure -C Debug -R "model_element_test\|euler_beam_3d_(local_stiffness\|transform_recovery)_test"` | 0 | 3/3 tests passed: `euler_beam_3d_local_stiffness_test`, `euler_beam_3d_transform_recovery_test`, `model_element_test`. |
+| 2 | `python -m unittest discover -s scripts -p "test_*.py"` | 0 | Ran 104 tests, OK. |
+| 3 | `python scripts/validate_workspace.py` with absolute CMake/CTest commands in `HARNESS_VALIDATION_COMMANDS` | 0 | Configure/build succeeded; CTest passed 23/23 tests; validation succeeded. |
+
+## Validation Results
+
+| validation_stage | result | evidence |
+| --- | --- | --- |
+| harness self-test | pass | `Ran 104 tests ... OK`. |
+| configure | pass | CMake generated `build/msvc-debug` with Visual Studio 17 2022 x64. |
+| build | pass | `fesa_core` and all unit/integration test executables built in Debug. |
+| CTest | pass | Full validation passed 23/23 tests. |
+| feature-specific tests | pass | 3/3 feature tests passed with the step 9 CTest filter. |
+
+## Failure Classification
+- classification: N/A
+- primary_failure: N/A
+- first_failed_command: N/A
+- evidence_tail: all required commands passed.
+
+## Failed Test Inventory
+
+| test_name | label | command | failure_summary |
+| --- | --- | --- | --- |
+| N/A | N/A | N/A | No failing tests. |
+
+## Scope Boundary
+This report only records build/test evidence for the kernel increment. It does not approve reference verification, physics sanity, or release readiness. No Abaqus reference solver was run and no reference CSV artifacts were generated or modified.
+
+## No-Change Assertion
+- source_files_modified: false in this report step.
+- test_files_modified: false in this report step.
+- cmake_files_modified: false.
+- reference_artifacts_modified: false.
+- notes: the report documents observed implementation changes from earlier steps but does not add new code.
+
+## Open Issues
+- Full solver release still requires parser integration, section/property mapping, assembly/static solver integration, HDF5 output, stored Abaqus reference artifacts, reference comparison, and physics sanity evidence.

docs/build-test-reports/solver-core-skeleton.md

@@ -0,0 +1,80 @@
+# Solver Core Skeleton Build/Test Report
+
+## Metadata
+- phase: solver-core-skeleton
+- scope: C++ skeleton classes only
+- status: passed
+
+## Commands Run
+
+```powershell
+python -m unittest scripts.test_header_declaration_only
+```
+
+- exit_code: 0
+- summary: Solver headers contain declarations only; function bodies are implemented in `.cpp` files.
+
+```powershell
+python -m unittest discover -s scripts -p "test_*.py"
+```
+
+- exit_code: 0
+- summary: 99 Python Harness tests passed.
+
+```powershell
+python scripts/validate_workspace.py
+```
+
+- exit_code: 0
+- summary: CMake configure, MSVC Debug build, and full CTest suite passed.
+
+```powershell
+ctest --test-dir build/msvc-debug --output-on-failure -C Debug -R solver_core_skeleton_integration_test
+```
+
+- exit_code: 0
+- summary: `solver_core_skeleton_integration_test` passed.
+
+## CTest Tests Added
+
+- `harness_smoke_test`
+- `core_diagnostic_test`
+- `core_ids_test`
+- `core_primitives_test`
+- `core_status_test`
+- `model_analysis_step_test`
+- `model_boundary_condition_test`
+- `model_domain_test`
+- `model_element_test`
+- `model_load_test`
+- `model_material_test`
+- `model_node_test`
+- `model_property_test`
+- `analysis_model_view_test`
+- `dof_manager_dof_key_test`
+- `dof_manager_numbering_test`
+- `analysis_state_vectors_test`
+- `analysis_flow_linear_static_analysis_test`
+- `analysis_flow_template_test`
+- `results_containers_test`
+- `solver_core_skeleton_integration_test`
+
+## Structural Tests Added
+
+- `scripts.test_header_declaration_only`
+
+## Build Structure
+
+- `fesa_core` now builds as a static library from `src/fesa/**/*.cpp`.
+- Solver headers under `src/fesa/**/*.hpp` declare functions only; method bodies live in matching `.cpp` translation units.
+
+## Known Limitations
+
+- No element stiffness, residual, tangent, or stress recovery calculation is implemented.
+- No material law evaluation is implemented.
+- No sparse assembly implementation is implemented beyond `DofManager` sparse pattern ownership.
+- No linear solver backend is implemented.
+- No HDF5 writer or reader is implemented.
+- No Abaqus `.inp` parser is implemented.
+- No reference comparison against Abaqus CSV artifacts is implemented.
+- `LinearStaticAnalysis` currently prepares the analysis model, DOF map, and zero-valued state only.

docs/formulations/euler-beam-3d-formulation.md

@@ -0,0 +1,210 @@
+# 3D Euler-Bernoulli Beam Kernel Formulation
+
+## Metadata
+- feature_id: euler-beam-3d
+- source_requirement: docs/requirements/euler-beam-3d.md
+- source_research: docs/research/euler-beam-3d-research.md
+- status: ready-for-numerical-review
+- owner_agent: formulation-agent
+- date: 2026-06-12
+
+## Scope and Assumptions
+- analysis_type: linear static kernel only.
+- element_type: two-node straight prismatic 3D Euler-Bernoulli beam.
+- deformation: small displacement and small rotation.
+- linearity: linear elastic.
+- material_model_boundary: constant `E` and `G`.
+- section_model_boundary: constant `A`, `J`, `Iy`, `Iz`.
+- coordinate_system: local Cartesian beam basis and global Cartesian basis.
+- units: user-consistent.
+
+## Primary Variables and DOFs
+- nodal_variables: displacement vector and rotation vector at each node.
+- dof_ordering: `[u1, v1, w1, rx1, ry1, rz1, u2, v2, w2, rx2, ry2, rz2]`.
+- local_axis_convention:
+  - local `x` is the normalized vector from node 1 to node 2.
+  - local `y` is the normalized projection of the user orientation vector onto the plane normal to local `x`.
+  - local `z = x cross y`.
+- sign_convention: internal local end force vector is `f_local = K_local * u_local` in the same DOF order.
+
+## Strong Form and Boundary Conditions
+The kernel is an element-level direct-stiffness implementation of the linear Euler-Bernoulli frame equations. For a prismatic member in local coordinate `x`:
+
+- axial equilibrium uses constant axial stiffness `EA`.
+- torsional equilibrium uses constant torsional stiffness `GJ`.
+- bending in local `x-y` uses flexural rigidity `EIz`.
+- bending in local `x-z` uses flexural rigidity `EIy`.
+
+Essential and natural boundary conditions are not applied by this kernel. They are downstream solver assembly concerns.
+
+## Weak or Variational Form
+The element stiffness represents the second variation of strain energy:
+
+```text
+U = 1/2 integral_0^L [
+      EA * (du/dx)^2
+    + GJ * (drx/dx)^2
+    + EIz * (d2v/dx2)^2
+    + EIy * (d2w/dx2)^2
+] dx
+```
+
+The element internal force is the derivative of this energy with respect to the nodal DOFs:
+
+```text
+f_local = K_local * u_local
+```
+
+External load vectors are out of scope for this increment.
+
+## Discretization
+- axial interpolation: linear two-node interpolation for `u`.
+- torsion interpolation: linear two-node interpolation for `rx`.
+- bending interpolation: cubic Hermite interpolation for each transverse displacement and corresponding nodal rotation.
+- nodal_layout: node 1 at local coordinate `x=0`, node 2 at `x=L`.
+- partition_of_unity_check: axial and torsion linear shape functions sum to one.
+- kronecker_delta_check: axial, torsion, and Hermite bending shape functions recover the nodal displacement/rotation DOFs.
+
+## Kinematics
+- axial strain: `epsilon_x = du/dx`.
+- torsion rate: `kappa_x = drx/dx`.
+- bending curvature for local `x-y` bending: curvature associated with `v` and rotations `rz`.
+- bending curvature for local `x-z` bending: curvature associated with `w` and rotations `ry`.
+- strain_measure: infinitesimal strain and small rotation curvature measures only.
+
+## Constitutive Contract
+- axial force result: `N = EA * du/dx`.
+- torsional moment result: `Mx = GJ * drx/dx`.
+- bending moment about local `z`: uses `EIz`.
+- bending moment about local `y`: uses `EIy`.
+- material_state_variables: none.
+
+## Element Equations
+Let:
+
+```text
+a  = E*A/L
+t  = G*J/L
+by = E*Iy
+bz = E*Iz
+cy1 = 12*by/L^3
+cy2 =  6*by/L^2
+cy3 =  4*by/L
+cy4 =  2*by/L
+cz1 = 12*bz/L^3
+cz2 =  6*bz/L^2
+cz3 =  4*bz/L
+cz4 =  2*bz/L
+```
+
+`K_local` is a 12x12 symmetric matrix initialized to zero, with these nonzero entries before symmetric mirroring:
+
+```text
+K(0,0)= a       K(0,6)=-a       K(6,6)= a
+
+K(3,3)= t       K(3,9)=-t       K(9,9)= t
+
+K(1,1)= cz1     K(1,5)= cz2     K(1,7)=-cz1    K(1,11)= cz2
+K(5,5)= cz3     K(5,7)=-cz2     K(5,11)= cz4
+K(7,7)= cz1     K(7,11)=-cz2    K(11,11)= cz3
+
+K(2,2)= cy1     K(2,4)=-cy2     K(2,8)=-cy1    K(2,10)=-cy2
+K(4,4)= cy3     K(4,8)= cy2     K(4,10)= cy4
+K(8,8)= cy1     K(8,10)= cy2    K(10,10)= cy3
+```
+
+The implementation must mirror upper-triangular entries to the lower triangle.
+
+## Mapping and Numerical Integration
+- reference_coordinates: local beam coordinate `x in [0,L]`.
+- isoparametric_mapping: straight member mapping only.
+- jacobian: `dx/dxi = L/2` if a parent coordinate is introduced; the direct closed-form stiffness does not require runtime quadrature.
+- determinant_checks: `L > 0` and finite.
+- gauss_points_and_weights: not used by the closed-form kernel.
+- integration_policy: analytical closed-form.
+
+## Transformation
+Let `R` be the 3x3 matrix whose rows are local basis vectors expressed in global coordinates:
+
+```text
+R = [ x_local^T
+      y_local^T
+      z_local^T ]
+```
+
+For each node, the same `R` maps global translational and rotational components to local components. The 12x12 transform `T` is block diagonal:
+
+```text
+u_local = T * u_global
+K_global = T^T * K_local * T
+```
+
+Global end-force recovery uses the same convention:
+
+```text
+u_local = T * u_global
+f_local = K_local * u_local
+f_global = T^T * f_local
+```
+
+## Output Recovery
+- displacement: input only for this kernel.
+- reaction: not recovered by this kernel; downstream assembly/constraints own reactions.
+- element_force: local and global element end forces recover from stiffness times displacement.
+- strain: out of scope.
+- stress: out of scope.
+- nodal_extrapolation: not applicable.
+
+## Invalid Input Handling
+- `length <= 0`, nonfinite length, or near-zero geometry length must throw `std::invalid_argument`.
+- nonpositive or nonfinite `E`, `G`, `A`, `J`, `Iy`, or `Iz` must throw `std::invalid_argument`.
+- zero orientation vector or orientation parallel to local `x` must throw `std::invalid_argument`.
+
+## Unit Test Tolerance
+- representative coefficient checks: absolute tolerance `1.0e-10` for the planned deterministic values.
+- symmetry checks: absolute tolerance `1.0e-10`.
+- rigid mode force norm checks: absolute tolerance `1.0e-9` for simple test magnitudes.
+
+## Algorithm Pseudocode
+
+```text
+validate length and section constants
+compute named stiffness coefficients
+initialize 12x12 local matrix to zero
+fill upper triangular axial, torsion, and bending terms
+mirror upper triangular terms
+
+for global operations:
+  compute normalized local x from node2 - node1
+  project orientation onto plane normal to local x and normalize as local y
+  compute local z = x cross y
+  build block diagonal T from R
+  compute K_global = T^T * K_local * T
+  recover forces with u_local = T*u_global, f_local = K_local*u_local, f_global = T^T*f_local
+```
+
+## Numerical Risks
+- rigid_body_modes: the unconstrained local/global stiffness should produce near-zero force for six rigid body modes; tests should include at least rigid translation in this increment.
+- patch_test: future solver-level patch tests require parser/assembly integration.
+- symmetry: local and global matrices must remain symmetric.
+- positive_definiteness: constrained systems may become positive definite; the isolated element is positive semidefinite.
+- hourglass: not applicable to this closed-form beam kernel.
+- shear_locking: excluded by Euler-Bernoulli assumptions and no shear deformation terms.
+- volumetric_locking: not applicable.
+- distortion: curved or offset geometry is out of scope.
+- singular_jacobian: represented by zero or near-zero length.
+
+## Open Issues and Downstream Handoff
+
+### Numerical Review Agent
+- Confirm signs for the local `w` and `ry` bending block.
+- Confirm the transformation convention and invalid orientation handling are implementation-ready.
+
+### I/O Definition Agent
+- Define future orientation input mapping and unsupported cases without claiming parser support in this kernel increment.
+
+### Reference Model Agent
+- Define future axial, torsion, bending, and skew-oriented reference models without generating artifacts.
+
+### Implementation Planning Agent
+- Transcribe the matrix entries and transformation convention into C++ tests before production implementation.

docs/implementation-plans/euler-beam-3d-implementation-plan.md

@@ -0,0 +1,127 @@
+# 3D Euler-Bernoulli Beam Kernel Implementation Plan
+
+## Metadata
+- feature_id: euler-beam-3d
+- source_requirement: docs/requirements/euler-beam-3d.md
+- source_research: docs/research/euler-beam-3d-research.md
+- source_formulation: docs/formulations/euler-beam-3d-formulation.md
+- source_numerical_review: docs/numerical-reviews/euler-beam-3d-review.md
+- source_io_definition: docs/io-definitions/euler-beam-3d-io.md
+- source_reference_models: docs/reference-models/euler-beam-3d-reference-models.md
+- status: ready-for-implementation
+- owner_agent: implementation-planning-agent
+- date: 2026-06-12
+
+## Readiness Check
+
+| input | required_status | observed_status | decision |
+| --- | --- | --- | --- |
+| requirement | sufficient draft | draft, kernel scope explicit | proceed |
+| research | ready-for-formulation | ready-for-formulation | proceed |
+| formulation | ready-for-numerical-review | ready-for-numerical-review | proceed |
+| numerical_review | pass-for-implementation-planning | pass-for-implementation-planning | proceed |
+| io_definition | ready-for-implementation-planning or explicit future scope | ready-for-implementation-planning | proceed |
+| reference_models | planned artifacts or future gate | needs-reference-artifacts | proceed for kernel only |
+
+## Implementation Scope
+- included_behavior:
+  - semantic `beam2` model topology.
+  - local 12x12 stiffness for a two-node 3D Euler-Bernoulli beam.
+  - local end-force recovery as `K_local * u_local`.
+  - local/global basis construction, global stiffness transform, and global end-force recovery.
+- excluded_behavior:
+  - parser integration, assembly, solver integration, HDF5 output, reference CSV generation, reference comparison, stress recovery, distributed loads, mass, geometric stiffness, releases, offsets, warping, shear deformation, nonlinear and dynamic behavior.
+- non_goals:
+  - no external linear algebra dependency.
+  - no CMake edits unless the existing source/test globs fail to pick up new files.
+
+## Work Breakdown
+
+| task_id | order | purpose | upstream_trace | depends_on | expected_test_first |
+| --- | --- | --- | --- | --- | --- |
+| EB3D-TASK-001 | 1 | Add `ElementTopology::beam2` semantic topology. | EB3D-REQ-001, EB3D-REQ-002 | none | EB3D-TEST-001 |
+| EB3D-TASK-002 | 2 | Add local stiffness and local end-force kernel. | EB3D-REQ-004, EB3D-REQ-005, EB3D-REQ-007, EB3D-REQ-009 | EB3D-TASK-001 | EB3D-TEST-002 |
+| EB3D-TASK-003 | 3 | Add global transform and global end-force recovery. | EB3D-REQ-006, EB3D-REQ-008, EB3D-REQ-010, EB3D-REQ-011 | EB3D-TASK-002 | EB3D-TEST-003 |
+| EB3D-TASK-004 | 4 | Record build/test evidence. | project validation policy | EB3D-TASK-003 | N/A |
+| EB3D-TASK-005 | 5 | Record release-readiness note with limitations. | EB3D-REQ-012 through EB3D-REQ-015 | EB3D-TASK-004 | N/A |
+
+## TDD Test Plan
+
+| test_id | order | test_type | red_condition | green_condition | linked_task | command |
+| --- | --- | --- | --- | --- | --- | --- |
+| EB3D-TEST-001 | 1 | unit | compile fails because `ElementTopology::beam2` is missing | model element test passes with `beam2` and existing topologies preserved | EB3D-TASK-001 | `ctest --test-dir build/msvc-debug --output-on-failure -C Debug -R model_element_test` |
+| EB3D-TEST-002 | 2 | unit | compile fails because `euler_beam_3d.hpp` API is missing | local stiffness entries, symmetry, `K*u`, and invalid inputs pass | EB3D-TASK-002 | `ctest --test-dir build/msvc-debug --output-on-failure -C Debug -R euler_beam_3d_local_stiffness_test` |
+| EB3D-TEST-003 | 3 | unit | compile fails because global transform API is missing | identity transform, rotated symmetry, rigid translation, axial force, and invalid orientation pass | EB3D-TASK-003 | `ctest --test-dir build/msvc-debug --output-on-failure -C Debug -R euler_beam_3d_transform_recovery_test` |
+
+## CMake/CTest Plan
+- target_candidates: existing FESA library target and unit-test executable pattern.
+- add_test_needs: only if existing test globs do not register new tests.
+- labels: unit.
+- msvc_config: Debug.
+- expected_feature_commands:
+  - `ctest --test-dir build/msvc-debug --output-on-failure -C Debug -R model_element_test`
+  - `ctest --test-dir build/msvc-debug --output-on-failure -C Debug -R euler_beam_3d_local_stiffness_test`
+  - `ctest --test-dir build/msvc-debug --output-on-failure -C Debug -R euler_beam_3d_transform_recovery_test`
+- workspace_validation: `python scripts/validate_workspace.py`.
+
+## Candidate Files and Ownership
+
+| file_candidate | purpose | owner_boundary | notes |
+| --- | --- | --- | --- |
+| `src/fesa/model/element.hpp` | add semantic topology enum value | model semantic layer | No equation IDs or section constants. |
+| `src/fesa/model/element.cpp` | only if topology string/conversion exists | model semantic layer | Preserve existing behavior. |
+| `tests/unit/model_element_test.cpp` | topology unit test | test-first for model change | Must be edited before production enum change. |
+| `src/fesa/elements/euler_beam_3d.hpp` | kernel declarations and small value types | element kernel layer | C++17 standard library only. |
+| `src/fesa/elements/euler_beam_3d.cpp` | local/global stiffness and recovery implementation | element kernel layer | No parser, assembly, or HDF5 code. |
+| `tests/unit/euler_beam_3d_local_stiffness_test.cpp` | local kernel tests | test-first for local kernel | Covers representative entries and invalid inputs. |
+| `tests/unit/euler_beam_3d_transform_recovery_test.cpp` | transform/recovery tests | test-first for global transform | Covers identity, rotation, rigid translation, axial force, invalid orientation. |
+
+## Acceptance Traceability Matrix
+
+| requirement_id | task_id | test_id | reference_model_id | acceptance_criterion | status |
+| --- | --- | --- | --- | --- | --- |
+| EB3D-REQ-001 | EB3D-TASK-001 | EB3D-TEST-001 | N/A | semantic topology exists for two-node beam | planned |
+| EB3D-REQ-002 | EB3D-TASK-002 | EB3D-TEST-002 | N/A | local DOF order matches formulation entries | planned |
+| EB3D-REQ-004 | EB3D-TASK-002 | EB3D-TEST-002 | N/A | section constants are used and invalid constants rejected | planned |
+| EB3D-REQ-005 | EB3D-TASK-002 | EB3D-TEST-002 | N/A | local 12x12 stiffness entries match formulation | planned |
+| EB3D-REQ-006 | EB3D-TASK-003 | EB3D-TEST-003 | `eb3d-skew-transform` future | global stiffness uses documented transform | planned |
+| EB3D-REQ-007 | EB3D-TASK-002 | EB3D-TEST-002 | N/A | local end forces equal `K*u` | planned |
+| EB3D-REQ-008 | EB3D-TASK-003 | EB3D-TEST-003 | `eb3d-skew-transform` future | global end forces use consistent transform | planned |
+| EB3D-REQ-009 | EB3D-TASK-002/003 | EB3D-TEST-002/003 | N/A | local and global stiffness matrices are symmetric | planned |
+| EB3D-REQ-010 | EB3D-TASK-003 | EB3D-TEST-003 | N/A | rigid global translation has near-zero end forces | planned |
+| EB3D-REQ-011 | EB3D-TASK-003 | EB3D-TEST-003 | `eb3d-skew-transform` future | transform and recovery conventions are consistent | planned |
+| EB3D-REQ-012 | EB3D-TASK-002/003 | EB3D-TEST-002/003 | N/A | excluded behaviors are not implemented | planned |
+| EB3D-REQ-013 | EB3D-TASK-004 | N/A | all future models | no `reference/` artifacts are modified | planned |
+| EB3D-REQ-014 | EB3D-TASK-005 | N/A | N/A | no full Abaqus compatibility claim | planned |
+| EB3D-REQ-015 | EB3D-TASK-005 | N/A | N/A | kernel completion remains separate from full solver release | planned |
+
+## Validation Commands
+
+```powershell
+ctest --test-dir build/msvc-debug --output-on-failure -C Debug -R model_element_test
+ctest --test-dir build/msvc-debug --output-on-failure -C Debug -R euler_beam_3d_local_stiffness_test
+ctest --test-dir build/msvc-debug --output-on-failure -C Debug -R euler_beam_3d_transform_recovery_test
+python -m unittest discover -s scripts -p "test_*.py"
+python scripts/validate_workspace.py
+```
+
+## Risks and Downstream Handoff
+
+### Implementation Agent
+- Follow `RED -> GREEN -> VERIFY` for each production change.
+- Keep matrix storage row-major with index `row * 12 + col`.
+- Throw `std::invalid_argument` for invalid length, section constants, and orientation.
+
+### Build/Test Executor Agent
+- Record CTest and workspace validation evidence after the C++ steps.
+
+### Correction Agent
+- Compile failures are likely missing CMake glob registration or signature mismatch.
+- Test failures are likely sign convention, transform convention, or tolerance mistakes.
+
+### Reference Verification Agent
+- No reference comparison is expected until future HDF5 and reference artifact work is approved.
+
+## Open Issues
+- Exact parser keyword subset and HDF5 schema are future integration work.
+- Reference artifacts remain unavailable by design in this kernel phase.

docs/io-definitions/euler-beam-3d-io.md

@@ -0,0 +1,94 @@
+# 3D Euler-Bernoulli Beam I/O Definition
+
+## Metadata
+- feature_id: euler-beam-3d
+- source_requirement: docs/requirements/euler-beam-3d.md
+- source_formulation: docs/formulations/euler-beam-3d-formulation.md
+- source_numerical_review: docs/numerical-reviews/euler-beam-3d-review.md
+- source_research: docs/research/euler-beam-3d-research.md
+- status: ready-for-implementation-planning
+- owner_agent: io-definition-agent
+- date: 2026-06-12
+
+## Abaqus Input Scope
+- input_format: Abaqus input file (`.inp`) in a future parser integration phase.
+- compatibility_disclaimer: FESA will support only the keyword subset defined for this feature; full Abaqus compatibility is not claimed.
+- kernel_increment_status: no parser, HDF5 writer, or reference comparison implementation is included in the current kernel phase.
+
+| keyword | support_status | level | required_parameters | mapped_internal_concept | notes |
+| --- | --- | --- | --- | --- | --- |
+| `*NODE` | future-supported | model | none | node id and global coordinates | Required for future beam geometry. |
+| `*ELEMENT` | future-supported | model | `TYPE` | two-node beam topology | Candidate Abaqus-equivalent type: `B31`; exact subset remains future parser work. |
+| `*MATERIAL` | future-supported | model | `NAME` | material identity | Required before solver integration. |
+| `*ELASTIC` | future-supported | model | none | `E`, Poisson ratio or `E,G` derivation policy | Future contract must define `G` source. |
+| beam section keyword | future-supported | model | section assignment | `A`, `J`, `Iy`, `Iz` | Exact Abaqus keyword subset is not implemented in this kernel phase. |
+| orientation data | future-supported | model | orientation vector | local axis construction | Must define nonparallel orientation vector. |
+| `*BOUNDARY` | future-supported | history | DOF range | constrained DOFs | DOF numbers `1..6` map to `ux,uy,uz,rx,ry,rz`. |
+| `*CLOAD` | future-supported | history | node, DOF, magnitude | concentrated load | Distributed loads are out of this kernel increment. |
+| `*STEP` / procedure | future-supported | history | procedure data | linear static step | Solver integration gate. |
+| output keywords | future-supported | history | output variables | HDF5 field/history requests | Authoritative output remains `results.h5`. |
+
+## Syntax Policy
+- case_insensitivity: future parser should treat Abaqus keyword names case-insensitively.
+- comma_separated_fields: future parser should follow the existing parser subset policy when this feature is integrated.
+- comment_lines: lines beginning with `**`.
+- unsupported_keywords: unsupported with a clear diagnostic unless a future I/O contract explicitly marks them ignored-with-warning.
+- ascii_assumption: input text policy follows the parser project contract.
+
+## Model Data Mapping
+- nodes: node label and three global Cartesian coordinates.
+- elements: two node labels and semantic topology `beam2`.
+- material: linear elastic constants sufficient to obtain `E` and `G`.
+- section: constants `A`, `J`, `Iy`, `Iz` assigned to the beam element set.
+- coordinates: global Cartesian coordinates plus a local orientation vector.
+- units: user-consistent and recorded in future metadata.
+
+## History Data Mapping
+- boundary_conditions: Abaqus DOF numbers map as `1=ux`, `2=uy`, `3=uz`, `4=rx`, `5=ry`, `6=rz`.
+- loads: future `*CLOAD` concentrated nodal loads may use the same DOF numbering.
+- output_requests: future nodal displacement, reaction, and element internal force requests must map to HDF5 datasets.
+
+## Internal Model Contract
+- element_type: two-node 3D Euler-Bernoulli beam semantic model, not an Abaqus compatibility claim.
+- connectivity: exactly two distinct nodes.
+- orientation: nonzero vector not parallel to the beam axis.
+- section_constants: positive finite `A`, `J`, `Iy`, `Iz`.
+- material_constants: positive finite `E` and `G`.
+- unsupported_cases: shear deformation, warping, releases, offsets, distributed loads, mass, geometric stiffness, nonlinear kinematics, dynamics, and thermal coupling.
+
+## Output HDF5 Schema
+This kernel phase does not write HDF5. Future solver integration should extend authoritative `results.h5` output with:
+
+| quantity | dataset_path | shape | dtype | required_attributes | location | notes |
+| --- | --- | --- | --- | --- | --- | --- |
+| displacement | `/steps/<step>/frames/<frame>/nodal/displacement` | `<nnode, 6>` | float64 | component names, units, coordinate system | nodal | Components `ux,uy,uz,rx,ry,rz`. |
+| reaction | `/steps/<step>/frames/<frame>/nodal/reaction` | `<nnode, 6>` | float64 | component names, units, coordinate system | nodal | Requires constraints/assembly. |
+| internal force | `/steps/<step>/frames/<frame>/element/internal_force` | `<nelem, 12>` | float64 | component names, units, coordinate system, element ids | element end | Components match beam local/global recovery contract. |
+| stress | `/steps/<step>/frames/<frame>/element/stress` | TBD | float64 | component names, units, location | TBD | Placeholder only; stress recovery is not approved in this kernel increment. |
+
+## FESA HDF5 to Reference CSV Comparison Schema
+Future reference comparison must read FESA HDF5 rows and compare them to Abaqus-generated CSV files under `reference/<model-id>/`.
+
+Required future CSV mappings:
+- `<model-id>_displacements.csv`: step, frame, node id, `ux`, `uy`, `uz`, and future rotational components if approved for comparison.
+- `<model-id>_reactions.csv`: step, frame, node id, force and moment reaction components.
+- `<model-id>_internalforces.csv`: step, frame, element id, end location, component, value.
+- `<model-id>_stresses.csv`: not required until stress recovery is approved.
+
+## Validation Rules
+- Reject elements whose connectivity does not contain exactly two distinct nodes.
+- Reject missing or nonpositive `E`, `G`, `A`, `J`, `Iy`, or `Iz`.
+- Reject missing, zero, or parallel orientation vectors when global transform construction is required.
+- Reject unsupported beam keywords with a clear diagnostic in future parser work.
+- Do not generate or modify Abaqus reference CSV files in this kernel phase.
+
+## Open Issues and Downstream Handoff
+
+### Reference Model Agent
+- Future `model.inp` files must stay within the subset documented here or record open parser issues.
+
+### Implementation Planning Agent
+- Current implementation may add only kernel-level C++ tests and source files. Parser/HDF5 work remains out of scope.
+
+### Reference Verification Agent
+- Use HDF5 as authoritative FESA output when future reference comparison is implemented.

docs/numerical-reviews/euler-beam-3d-review.md

@@ -0,0 +1,76 @@
+# 3D Euler-Bernoulli Beam Kernel Numerical Review
+
+## Metadata
+- feature_id: euler-beam-3d
+- source_formulation: docs/formulations/euler-beam-3d-formulation.md
+- status: pass-for-implementation-planning
+- owner_agent: numerical-review-agent
+- date: 2026-06-12
+
+## Review Verdict
+- verdict: pass-for-implementation-planning
+- reason: the formulation has a bounded kernel scope, explicit DOF order, closed-form local stiffness, unambiguous transform convention, and concrete invalid-input behavior.
+
+## Critical Findings
+- No blocking numerical defects were found for kernel implementation planning.
+- The formulation is not a full solver release contract. Parser, assembly, HDF5 output, reference comparison, and physics sanity remain downstream gates.
+
+## Numerical Risk Assessment
+- rigid_body_modes: expected six zero-energy modes for the unconstrained element. Implementation tests must include at least rigid global translation and local/global force consistency.
+- patch_test: solver-level patch testing is blocked until assembly and parser integration exist.
+- symmetry: local stiffness is symmetric by construction; global stiffness remains symmetric if `K_global = T^T K_local T` is used consistently.
+- positive_definiteness: isolated element stiffness is positive semidefinite, not positive definite. Positive definiteness requires enough constraints in a solver model.
+- hourglass: not applicable to closed-form Euler-Bernoulli beam stiffness.
+- shear_locking: Timoshenko shear terms are excluded, so shear locking is not introduced by this kernel increment.
+- volumetric_locking: not applicable.
+- distortion: curved geometry and offsets are out of scope; zero or near-zero length must be rejected.
+- singular_jacobian: zero length is the relevant singular mapping case.
+- conditioning: very slender elements or large stiffness ratios can be ill-conditioned; tests should use moderate deterministic values.
+- convergence: not applicable to this linear element kernel.
+
+## Consistency Checks
+- units: pass. `EA/L`, `GJ/L`, `EI/L^3`, `EI/L^2`, and `EI/L` have the expected force/displacement or moment/rotation dimensions for the matching DOF pairs.
+- dimensions: pass. The element vector is 12x1 and stiffness is 12x12.
+- signs: pass for the documented DOF order; local `v-rz` and `w-ry` coupling signs are explicitly stated.
+- dof_ordering: pass. The formulation keeps node 1 DOFs followed by node 2 DOFs.
+- coordinate_transforms: pass. `u_local = T*u_global`, `K_global = T^T*K_local*T`, and `f_global = T^T*f_local` are mutually consistent.
+- matrix_vector_dimensions: pass.
+- integration_weights: not applicable because the kernel uses closed-form stiffness.
+- output_locations: pass for element end forces only; stress/strain recovery is explicitly out of scope.
+
+## Verification Readiness
+- unit_tests:
+  - model topology accepts `ElementTopology::beam2`.
+  - local stiffness representative entries match named coefficients.
+  - local stiffness is symmetric within `1.0e-10`.
+  - local force recovery equals `K*u`.
+  - invalid length and nonpositive section constants throw `std::invalid_argument`.
+  - axis-aligned transform produces global stiffness equal to local stiffness.
+  - rotated global stiffness remains symmetric.
+  - rigid global translation produces near-zero global end forces.
+  - simple axial extension produces equal and opposite axial end forces.
+  - parallel orientation vector throws `std::invalid_argument`.
+- patch_tests: future solver integration gate.
+- mms_or_mes: not required for this closed-form kernel increment.
+- benchmark_reference_comparison: future gate after reference artifacts exist.
+- missing_evidence: no missing evidence blocks kernel implementation.
+
+## Required Revisions
+
+### Formulation Agent
+- None before implementation planning.
+
+### Research Agent
+- None before implementation planning.
+
+### Reference Model Agent
+- Future reference model artifacts remain required before release readiness, but they are not required for this kernel implementation.
+
+## Downstream Handoff
+
+### Implementation Planning Agent
+- Use the verification readiness list as the minimum C++ TDD checklist.
+- Keep tests deterministic and avoid ill-conditioned constants.
+
+### Reference Model Agent
+- Define future solution-verification models separately from these kernel-only unit tests.

docs/reference-models/euler-beam-3d-reference-models.md

@@ -0,0 +1,136 @@
+# 3D Euler-Bernoulli Beam Reference Models
+
+## Metadata
+- feature_id: euler-beam-3d
+- source_requirement: docs/requirements/euler-beam-3d.md
+- source_research: docs/research/euler-beam-3d-research.md
+- source_formulation: docs/formulations/euler-beam-3d-formulation.md
+- source_numerical_review: docs/numerical-reviews/euler-beam-3d-review.md
+- source_io_definition: docs/io-definitions/euler-beam-3d-io.md
+- status: needs-reference-artifacts
+- owner_agent: reference-model-agent
+- date: 2026-06-12
+
+## Reference Strategy
+- verification_scope: future solver-level verification after parser, assembly, solve, HDF5 output, and comparison tooling exist.
+- code_verification: current kernel unit tests cover local stiffness, transformation, and force recovery without external reference artifacts.
+- solution_verification: future cantilever and skew-orientation models will compare displacement, reaction, and internal force quantities.
+- benchmark_reference_comparison: future Abaqus-generated CSV files under `reference/<model-id>/` are required before release readiness.
+- excluded_validation_scope: physical experiment validation is not in scope.
+
+## Model Inventory
+
+| model_id | category | purpose | status | required_artifacts |
+| --- | --- | --- | --- | --- |
+| `eb3d-axial-cantilever` | analytical | axial bar-as-beam displacement, reaction, and axial end force | needs-reference-artifacts | `model.inp`, `metadata.json`, displacement, reaction, internal force CSV, README |
+| `eb3d-torsion-cantilever` | analytical | torsional rotation, reaction moment, and torsional end moment | needs-reference-artifacts | `model.inp`, `metadata.json`, displacement/rotation, reaction, internal force CSV, README |
+| `eb3d-bend-y-cantilever` | analytical | bending response using `EIy` | needs-reference-artifacts | `model.inp`, `metadata.json`, displacement, reaction, internal force CSV, README |
+| `eb3d-bend-z-cantilever` | analytical | bending response using `EIz` | needs-reference-artifacts | `model.inp`, `metadata.json`, displacement, reaction, internal force CSV, README |
+| `eb3d-skew-transform` | regression | skew-oriented beam transform and component mapping | needs-reference-artifacts | `model.inp`, `metadata.json`, displacement, reaction, internal force CSV, README |
+
+## Model Records
+
+### `eb3d-axial-cantilever`
+- category: analytical.
+- verified_requirements: EB3D-REQ-005, EB3D-REQ-007, future solver I/O requirements.
+- analysis_type: linear static.
+- element_type: two-node 3D Euler-Bernoulli beam, Abaqus-equivalent candidate `B31`.
+- expected_physical_quantities: axial displacement, support reaction, element axial end force.
+- artifact_status: needs-reference-artifacts.
+
+### `eb3d-torsion-cantilever`
+- category: analytical.
+- verified_requirements: EB3D-REQ-005, EB3D-REQ-007.
+- analysis_type: linear static.
+- expected_physical_quantities: twist, reaction moment about beam axis, torsional end moment.
+- artifact_status: needs-reference-artifacts.
+
+### `eb3d-bend-y-cantilever`
+- category: analytical.
+- verified_requirements: EB3D-REQ-005, EB3D-REQ-007.
+- analysis_type: linear static.
+- expected_physical_quantities: transverse displacement, support reaction, bending end forces using `EIy`.
+- artifact_status: needs-reference-artifacts.
+
+### `eb3d-bend-z-cantilever`
+- category: analytical.
+- verified_requirements: EB3D-REQ-005, EB3D-REQ-007.
+- analysis_type: linear static.
+- expected_physical_quantities: transverse displacement, support reaction, bending end forces using `EIz`.
+- artifact_status: needs-reference-artifacts.
+
+### `eb3d-skew-transform`
+- category: regression.
+- verified_requirements: EB3D-REQ-006, EB3D-REQ-008, EB3D-REQ-011.
+- analysis_type: linear static.
+- expected_physical_quantities: transformed displacement, reaction, and internal force components.
+- artifact_status: needs-reference-artifacts.
+
+## Abaqus Input Requirements
+- input_file: `reference/<model-id>/model.inp`.
+- supported_keyword_subset: must follow docs/io-definitions/euler-beam-3d-io.md when future parser work is approved.
+- model_data: nodes, one or more beam elements, material, section, orientation, and units.
+- history_data: linear static step, boundary conditions, concentrated loads or moments, and output requests.
+- unsupported_keyword_policy: unsupported unless explicitly added by a future I/O contract.
+
+## Artifact Bundle Contract
+
+```text
+reference/
+  <model-id>/
+    model.inp
+    metadata.json
+    <model-id>_displacements.csv
+    <model-id>_reactions.csv
+    <model-id>_internalforces.csv
+    README.md
+```
+
+`<model-id>_stresses.csv` is not required until stress recovery is approved. No artifact in this bundle is created or modified in the current phase.
+
+## Metadata JSON Contract
+Future `metadata.json` files must include:
+
+- `feature_id`
+- `model_id`
+- `artifact_status`
+- `input_file`
+- Abaqus version/source and generation owner
+- generation date
+- source documents
+- units and coordinate system
+- analysis type and element types
+- material and section values
+- boundary and load summaries
+- output requests
+- reference CSV schema version
+- reference CSV files
+- tolerance policy
+- limitations
+
+## Coverage Matrix
+
+| requirement_id | model_id | compared_quantity | fesa_hdf5_dataset | reference_csv | tolerance | verification_method | status |
+| --- | --- | --- | --- | --- | --- | --- | --- |
+| EB3D-REQ-005 | `eb3d-axial-cantilever` | displacement/internal force | `/steps/<step>/frames/<frame>/nodal/displacement`, `/element/internal_force` | displacement/internal force CSV | TBD | hdf5-to-reference-csv | needs-reference-artifacts |
+| EB3D-REQ-006 | `eb3d-skew-transform` | displacement/reaction/internal force | displacement, reaction, internal force HDF5 datasets | displacement/reaction/internal force CSV | TBD | hdf5-to-reference-csv | needs-reference-artifacts |
+| EB3D-REQ-007 | `eb3d-torsion-cantilever` | rotation/reaction/internal force | displacement, reaction, internal force HDF5 datasets | displacement/reaction/internal force CSV | TBD | hdf5-to-reference-csv | needs-reference-artifacts |
+| EB3D-REQ-008 | `eb3d-skew-transform` | global end force | `/steps/<step>/frames/<frame>/element/internal_force` | internal force CSV | TBD | hdf5-to-reference-csv | needs-reference-artifacts |
+| EB3D-REQ-011 | `eb3d-skew-transform` | transform consistency | displacement and internal force datasets | displacement/internal force CSV | TBD | hdf5-to-reference-csv | needs-reference-artifacts |
+
+## Artifact Acceptance Checklist
+- No files under `reference/` are created or modified in this phase.
+- Every future model must document provenance before use as verification evidence.
+- Required CSV files keep the model status at `needs-reference-artifacts` until generated through an approved reference phase.
+- Stress CSV files remain out of scope until stress recovery is approved.
+
+## Open Issues and Downstream Handoff
+
+### I/O Definition Agent
+- Finalize the exact beam section and orientation keyword subset before creating `model.inp` files.
+
+### Implementation Planning Agent
+- Kernel TDD may proceed without reference artifacts because this phase only implements element-level matrix and recovery routines.
+
+### Reference Verification Agent
+- Do not compare results until authoritative HDF5 output and approved CSV artifacts exist.

docs/releases/euler-beam-3d-release.md

@@ -0,0 +1,131 @@
+# 3D Euler-Bernoulli Beam Kernel Release Readiness Note
+
+## Metadata
+- feature_id: euler-beam-3d
+- source_requirement: docs/requirements/euler-beam-3d.md
+- source_formulation: docs/formulations/euler-beam-3d-formulation.md
+- source_numerical_review: docs/numerical-reviews/euler-beam-3d-review.md
+- source_io_definition: docs/io-definitions/euler-beam-3d-io.md
+- source_reference_model: docs/reference-models/euler-beam-3d-reference-models.md
+- source_implementation_plan: docs/implementation-plans/euler-beam-3d-implementation-plan.md
+- source_build_test_report: docs/build-test-reports/euler-beam-3d-build-test.md
+- source_reference_verification_report: N/A
+- source_physics_evaluation_report: N/A
+- status: not-release-ready-kernel-increment-complete
+- owner_agent: release-agent
+- date: 2026-06-12
+
+## Release Scope
+
+| item | included | excluded | notes |
+| --- | --- | --- | --- |
+| analysis_type | element kernel calculations | full linear static solve | Solver integration is not complete. |
+| element_type | two-node straight prismatic 3D Euler-Bernoulli beam kernel | Timoshenko beams, curved beams, offsets, releases, warping | Kernel scope only. |
+| material_model | positive finite `E` and `G` constants supplied to the kernel | material parser or library integration | No material database integration. |
+| section_model | positive finite `A`, `J`, `Iy`, `Iz` constants supplied to the kernel | beam section parser/property integration | Future I/O phase required. |
+| output_quantities | local/global stiffness and local/global end forces returned in memory | HDF5 results, reference CSV views, stress recovery | Authoritative output path is future work. |
+
+## Completed Kernel Increment
+- Added semantic `ElementTopology::beam2`.
+- Added local 12x12 Euler-Bernoulli beam stiffness for axial, torsion, and two bending planes.
+- Added local end-force recovery as `K_local * u_local`.
+- Added local axis construction from two nodes plus an orientation vector.
+- Added global stiffness transformation using `K_global = T^T K_local T`.
+- Added global end-force recovery using the same transformation convention.
+- Added unit tests for topology, local stiffness entries, local symmetry, local force recovery, invalid local inputs, identity transform, rotated symmetry, rigid translation, axial extension, parallel orientation, and zero orientation.
+
+## Missing For Full Feature Release
+- Parser implementation for the approved Abaqus subset.
+- Beam section/property semantic model integration.
+- Material-to-section mapping for `E`, `G`, `A`, `J`, `Iy`, and `Iz`.
+- Assembler and static solver integration.
+- HDF5 result emission for beam displacement, reaction, and internal force quantities.
+- Stored Abaqus reference artifacts under `reference/<model-id>/`.
+- Reference comparison report.
+- Physics sanity report.
+- Final release report with all upstream gates passing.
+
+## Gate Evidence Inventory
+
+| gate | source | expected_status | observed_status | verdict |
+| --- | --- | --- | --- | --- |
+| requirements | docs/requirements/euler-beam-3d.md | sufficient for kernel | draft with explicit kernel scope | pass for kernel |
+| research | docs/research/euler-beam-3d-research.md | ready-for-formulation | ready-for-formulation | pass |
+| formulation | docs/formulations/euler-beam-3d-formulation.md | ready-for-numerical-review | ready-for-numerical-review | pass |
+| numerical_review | docs/numerical-reviews/euler-beam-3d-review.md | pass-for-implementation-planning | pass-for-implementation-planning | pass |
+| io_definition | docs/io-definitions/euler-beam-3d-io.md | future integration contract | ready-for-implementation-planning | pass for planning |
+| reference_model | docs/reference-models/euler-beam-3d-reference-models.md | artifacts present for release | needs-reference-artifacts | fail for release |
+| implementation | source and tests | implemented kernel | implemented kernel | pass for kernel |
+| build_test | docs/build-test-reports/euler-beam-3d-build-test.md | pass | pass-for-reference-verification | pass |
+| reference_verification | docs/reference-verifications/euler-beam-3d-reference-verification.md | pass-for-physics-evaluation | missing | fail for release |
+| physics_evaluation | docs/physics-evaluations/euler-beam-3d-physics-evaluation.md | pass-for-release-agent | missing | fail for release |
+
+## Acceptance Traceability
+
+| requirement_id | release_disposition |
+| --- | --- |
+| EB3D-REQ-001 through EB3D-REQ-011 | Implemented and unit-tested at kernel level. Full solver evidence deferred. |
+| EB3D-REQ-012 | Excluded behaviors remain out of scope. |
+| EB3D-REQ-013 | No reference artifacts were created or modified. |
+| EB3D-REQ-014 | No full Abaqus compatibility claim is made. |
+| EB3D-REQ-015 | Kernel completion is explicitly separate from full solver release readiness. |
+
+## Validation Evidence
+
+| command_or_report | expected | observed | notes |
+| --- | --- | --- | --- |
+| `ctest --test-dir build/msvc-debug --output-on-failure -C Debug -R "model_element_test\|euler_beam_3d_(local_stiffness\|transform_recovery)_test"` | pass | pass, 3/3 tests | Feature-specific CTest filter passed. |
+| `python -m unittest discover -s scripts -p "test_*.py"` | pass | pass, 104 tests | Harness self-test passed. |
+| `python scripts/validate_workspace.py` | pass | pass, 23/23 CTests | Used absolute CMake/CTest commands through `HARNESS_VALIDATION_COMMANDS`. |
+
+## Known Limitations
+
+| limitation | category | user_impact | disposition |
+| --- | --- | --- | --- |
+| No parser support for beam input yet. | input | Users cannot run beam `.inp` models end-to-end. | deferred |
+| No section/property model integration. | model | Constants must be supplied directly to the kernel API. | deferred |
+| No assembly/static solver integration. | solver | Beam elements do not contribute to a global system yet. | deferred |
+| No HDF5 result output for beam quantities. | output | No official solver output exists for beam results. | deferred |
+| No reference artifacts or comparison report. | verification | Full solver correctness against Abaqus is not established. | blocker for release |
+| No stress or strain recovery. | physics/output | Stress CSV and stress HDF5 datasets are not available. | deferred |
+| No shear deformation, warping, releases, offsets, mass, geometric stiffness, nonlinear, dynamic, or thermal behavior. | physics | Feature is limited to the approved Euler-Bernoulli kernel. | documented |
+
+## Release Notes Draft
+
+### Feature Summary
+- Kernel increment for a two-node 3D Euler-Bernoulli beam element is complete.
+
+### Verification Scope
+- Unit-tested local stiffness, local end forces, global transform, global end forces, and invalid geometry/section handling.
+
+### Main Limitations
+- This is not an end-to-end solver feature. Parser, assembly, HDF5 output, reference comparison, and physics sanity remain future gates.
+
+### Artifacts
+- docs/build-test-reports/euler-beam-3d-build-test.md
+- docs/requirements/euler-beam-3d.md
+- docs/formulations/euler-beam-3d-formulation.md
+- docs/reference-models/euler-beam-3d-reference-models.md
+
+## Release Verdict
+- verdict: not-release-ready-kernel-increment-complete
+- reason: kernel implementation and build/test validation are complete, but full feature release lacks parser integration, HDF5 output, reference artifacts, reference verification, and physics sanity evidence.
+
+## Handoff Recommendation
+
+| target_agent | reason | required_input |
+| --- | --- | --- |
+| Coordinator Agent | Decide next phase sequencing. | This release note and phase index. |
+| I/O Definition Agent | Finalize exact Abaqus beam subset and section/orientation keywords. | Current I/O contract and kernel API behavior. |
+| Reference Model Agent | Generate approved future reference artifact contracts and then artifacts in an explicit phase. | Reference model contract. |
+| Implementation Agent | Integrate kernel into parser/model/assembly only after upstream contracts are approved. | Implementation plan and build/test report. |
+
+## No-Change Assertion
+- source_files_modified: false in this release note step.
+- test_files_modified: false in this release note step.
+- cmake_files_modified: false.
+- reference_artifacts_modified: false.
+- tolerance_policies_modified: false.
+
+## Open Issues
+- Full solver release remains blocked until the missing gates listed above are completed.

docs/requirements/euler-beam-3d.md

@@ -0,0 +1,201 @@
+# 3D Euler-Bernoulli Beam Kernel Requirements
+
+## Metadata
+- feature_id: euler-beam-3d
+- title: Kernel-first 3D Euler-Bernoulli beam element
+- status: draft
+- owner_agent: requirement-agent
+- date: 2026-06-12
+
+## Purpose
+Define the first implementable requirements baseline for a two-node, straight, prismatic, small-displacement 3D Euler-Bernoulli beam element kernel. This increment is limited to element-level stiffness and end-force recovery so downstream formulation, numerical review, and C++ implementation can proceed without claiming full solver readiness.
+
+## Assumptions
+- The element is a two-node straight beam with nonzero length.
+- The element uses small-displacement, small-rotation, linear elastic kinematics.
+- Each node has six mechanical DOFs in this order: `ux, uy, uz, rx, ry, rz`.
+- The element displacement vector order is node 1 DOFs followed by node 2 DOFs, each in `ux, uy, uz, rx, ry, rz` order.
+- Section and material constants are constant over the element: `E`, `G`, `A`, `J`, `Iy`, `Iz`.
+- `Iy` and `Iz` are second moments of area about the element local principal bending axes, and `J` is the torsional constant about the local beam axis.
+- Units are user-consistent; the kernel does not perform unit conversion in this increment.
+- The local element basis is right-handed and orthonormal. The local `x` axis follows node 1 to node 2; the source of local `y` and `z` orientation for parser or solver integration remains open.
+- Boundary conditions, global equation numbering, sparse assembly, linear solve, HDF5 writing, and reference comparison are downstream solver concerns, not kernel responsibilities in this increment.
+
+## In Scope
+- Two-node, straight, prismatic 3D Euler-Bernoulli beam element kernel.
+- Six mechanical DOFs per node in `ux, uy, uz, rx, ry, rz` order.
+- Linear elastic section constants: `E`, `G`, `A`, `J`, `Iy`, `Iz`.
+- Local 12x12 stiffness matrix support.
+- Global 12x12 stiffness matrix support when a valid local-to-global transformation is supplied.
+- Local and global element end-force recovery from compatible 12-entry nodal displacement vectors.
+- Kernel-level verification hooks for stiffness symmetry, rigid body modes, and local/global transformation consistency.
+
+## Out Of Scope
+- Shear deformation, Timoshenko beam behavior, shear correction factors, and shear locking treatment.
+- Warping DOFs, warping stiffness, restrained warping, or bimoment recovery.
+- End releases, offsets, rigid links, eccentric sections, or joint constraints.
+- Distributed loads, equivalent nodal load vectors, body forces, or pressure loads.
+- Mass matrix, damping matrix, geometric stiffness, buckling, nonlinear kinematics, dynamics, or thermal coupling.
+- Abaqus reference CSV generation or modification.
+- Full Abaqus compatibility claims.
+- Parser integration, end-to-end assembly, solver execution, HDF5 result output, reference comparison, physics sanity, or release readiness.
+
+## Analysis Definition
+- analysis_type: linear static kernel only
+- element_type: two-node 3D Euler-Bernoulli beam
+- nodes_per_element: 2
+- dofs_per_node: 6
+- dof_order: `ux, uy, uz, rx, ry, rz`
+- element_vector_size: 12
+- material_model: linear elastic constants `E` and `G`
+- section_model: constant section constants `A`, `J`, `Iy`, `Iz`
+- coordinate_system: local Cartesian beam frame and global Cartesian frame
+- loads: not in this increment
+- boundary_conditions: not in this increment
+- official_solver_output: not in this increment
+
+## Input Requirements
+- The kernel must receive or derive two distinct node positions sufficient to determine element length and local `x` direction.
+- The kernel must receive positive finite constants for `E`, `G`, `A`, `J`, `Iy`, and `Iz`.
+- The kernel must receive a valid right-handed orthonormal local basis or equivalent transformation data before producing a global stiffness matrix or global end-force vector.
+- The kernel must receive 12-entry nodal displacement vectors ordered as node 1 then node 2, with each node using `ux, uy, uz, rx, ry, rz`.
+
+## Output Requirements
+- The kernel must produce a local 12x12 stiffness matrix in the documented element DOF order.
+- The kernel must produce a global 12x12 stiffness matrix in the documented global DOF order when supplied with a valid local-to-global transformation.
+- The kernel must recover local element end forces from a compatible local nodal displacement vector.
+- The kernel must recover global element end forces from a compatible global nodal displacement vector when supplied with the same transformation used for the global stiffness matrix.
+- This increment must not write `results.h5` and must not create deterministic CSV views or Abaqus reference CSV files.
+
+## Verification Quantities
+| quantity | required_for_kernel | verification intent | downstream owner |
+| --- | --- | --- | --- |
+| displacement DOFs | yes | Confirm 12-entry vectors use the required `ux, uy, uz, rx, ry, rz` per-node ordering. | Implementation Planning Agent |
+| reactions / end forces | yes | Confirm local and global end-force recovery is consistent with the corresponding stiffness matrix and displacement vector. | Formulation Agent; Implementation Planning Agent |
+| stiffness symmetry | yes | Confirm local and global 12x12 stiffness matrices are symmetric within the formulation-defined numeric tolerance. | Numerical Review Agent; Implementation Planning Agent |
+| rigid body modes | yes | Confirm an unconstrained isolated element has six zero-energy rigid body modes within the formulation-defined numeric tolerance. | Formulation Agent; Numerical Review Agent |
+| local/global transformation consistency | yes | Confirm transformed stiffness and force recovery are mutually consistent for the same basis and displacement vector. | Formulation Agent; Implementation Planning Agent |
+| stress | no | Stress recovery is not part of this kernel increment. | Future requirements |
+| strain | no | Strain recovery is not part of this kernel increment. | Future requirements |
+| energy or residual | downstream hook | Strain energy consistency may be added as a kernel verification hook after formulation review. | Numerical Review Agent |
+
+## Tolerance Policy
+- absolute_tolerance: needs-downstream-decision by Numerical Review Agent and Implementation Planning Agent
+- relative_tolerance: needs-downstream-decision by Numerical Review Agent and Implementation Planning Agent
+- norm_based_tolerance: needs-downstream-decision by Numerical Review Agent and Implementation Planning Agent
+- Units: user-consistent structural units; no unit conversion in this increment.
+- Coordinate system: local and global Cartesian frames with a right-handed orthonormal element basis.
+- Rule: any numerical kernel test added later must state the tolerance source and whether it checks component-wise values, matrix norms, eigenvalue/residual norms, or force vector norms.
+
+## Reference Artifact Requirements
+Reference artifacts are not required for kernel completion in this increment.
+
+Future end-to-end solver verification may use a bundle under `reference/<model-id>/`, but this requirements baseline does not authorize creating, modifying, or regenerating Abaqus reference CSV files.
+
+| artifact | kernel increment status | future solver-release status |
+| --- | --- | --- |
+| `model.inp` | not required | required after parser and assembly scope is approved |
+| `metadata.json` | not required | required with provenance for approved reference artifacts |
+| `<model-id>_displacements.csv` | not required | required for reference comparison |
+| `<model-id>_reactions.csv` | not required | required for reaction/end-force comparison |
+| `<model-id>_internalforces.csv` | not required | required when internal force comparison scope is approved |
+| `<model-id>_stresses.csv` | not required | not required until stress recovery is in scope |
+
+## Must Requirements
+| id | requirement |
+| --- | --- |
+| EB3D-REQ-001 | The beam kernel must model only a two-node, straight, prismatic, small-displacement 3D Euler-Bernoulli beam element in this increment. |
+| EB3D-REQ-002 | The beam kernel must use six mechanical DOFs per node ordered as `ux, uy, uz, rx, ry, rz`. |
+| EB3D-REQ-003 | The beam kernel must use a 12-entry element vector ordered as all node 1 DOFs followed by all node 2 DOFs. |
+| EB3D-REQ-004 | The beam kernel must accept linear elastic constants `E` and `G` and section constants `A`, `J`, `Iy`, and `Iz`. |
+| EB3D-REQ-005 | The beam kernel must produce a local 12x12 stiffness matrix in the documented local DOF order. |
+| EB3D-REQ-006 | The beam kernel must produce a global 12x12 stiffness matrix when supplied with a valid right-handed orthonormal local-to-global transformation. |
+| EB3D-REQ-007 | The beam kernel must recover local element end forces from a compatible local nodal displacement vector. |
+| EB3D-REQ-008 | The beam kernel must recover global element end forces from a compatible global nodal displacement vector and the same transformation convention used for global stiffness. |
+| EB3D-REQ-009 | The local and global stiffness matrices must be symmetric within the formulation-defined tolerance. |
+| EB3D-REQ-010 | The unconstrained isolated element stiffness must preserve six rigid body modes within the formulation-defined tolerance. |
+| EB3D-REQ-011 | The local and global stiffness and end-force recovery paths must be mutually consistent under the same transformation. |
+| EB3D-REQ-012 | The kernel increment must not include shear deformation, warping, end releases, offsets, distributed loads, mass matrix, geometric stiffness, nonlinear kinematics, dynamics, or thermal coupling. |
+| EB3D-REQ-013 | The kernel increment must not create, modify, or regenerate Abaqus reference CSV files. |
+| EB3D-REQ-014 | The requirements and downstream implementation must not claim full Abaqus compatibility for this feature. |
+| EB3D-REQ-015 | The kernel completion criteria must be separate from full solver release readiness. |
+
+## Requirement Verification Matrix
+| id | source | priority | verification_method | acceptance_criteria | tolerance | downstream_agents | status |
+| --- | --- | --- | --- | --- | --- | --- | --- |
+| EB3D-REQ-001 | user step scope | must | formulation review; C++ unit tests in later implementation steps | Formulation and implementation address only the scoped Euler-Bernoulli beam kernel and do not add excluded behaviors. | not numeric | Formulation Agent; Implementation Planning Agent | draft |
+| EB3D-REQ-002 | user step scope | must | unit test; code review | Tests or implementation checks demonstrate the exact per-node DOF order `ux, uy, uz, rx, ry, rz`. | not numeric | Implementation Planning Agent | draft |
+| EB3D-REQ-003 | user step scope | must | unit test; code review | Tests or implementation checks demonstrate the 12-entry node 1 then node 2 vector order. | not numeric | Implementation Planning Agent | draft |
+| EB3D-REQ-004 | user step scope | must | input validation test in later implementation steps | Nonpositive or nonfinite constants are rejected or diagnosed according to the downstream implementation contract. | not numeric for this baseline | I/O Definition Agent; Implementation Planning Agent | draft |
+| EB3D-REQ-005 | user step scope | must | formulation review; local stiffness unit test | Local stiffness is 12x12, uses the documented local DOF order, and matches formulation-specified benchmark values. | formulation-defined | Formulation Agent; Numerical Review Agent; Implementation Planning Agent | draft |
+| EB3D-REQ-006 | user step scope | must | transformation unit test | Global stiffness is 12x12 and matches the formulation-specified transformation of the local matrix. | formulation-defined | Formulation Agent; Implementation Planning Agent | draft |
+| EB3D-REQ-007 | user step scope | must | local force recovery unit test | Local end-force output is consistent with the local stiffness matrix and local displacement vector. | formulation-defined | Formulation Agent; Implementation Planning Agent | draft |
+| EB3D-REQ-008 | user step scope | must | global force recovery unit test | Global end-force output is consistent with global stiffness, global displacement, and the local/global transformation convention. | formulation-defined | Formulation Agent; Implementation Planning Agent | draft |
+| EB3D-REQ-009 | user verification notes | must | matrix symmetry unit test | Local and global stiffness matrices satisfy symmetry checks within the selected matrix tolerance. | needs-downstream-decision | Numerical Review Agent; Implementation Planning Agent | draft |
+| EB3D-REQ-010 | user verification notes | must | rigid body mode unit test; numerical review | Six rigid body modes produce zero or tolerance-bounded internal force/energy for an unconstrained isolated element. | needs-downstream-decision | Formulation Agent; Numerical Review Agent | draft |
+| EB3D-REQ-011 | user verification notes | must | transformation consistency unit test | Stiffness transformation and end-force transformation use one consistent convention for the same local basis. | formulation-defined | Formulation Agent; Implementation Planning Agent | draft |
+| EB3D-REQ-012 | user exclusions | must | scope review; code review | No excluded behavior appears in requirements, formulation, or kernel implementation for this increment. | not numeric | Coordinator Agent; Implementation Planning Agent | draft |
+| EB3D-REQ-013 | user forbidden list | must | git diff review | No file under `reference/` and no Abaqus reference CSV is created or modified in this increment. | not numeric | Reference Model Agent; Coordinator Agent | draft |
+| EB3D-REQ-014 | ADR-003; user verification notes | must | documentation review; I/O contract review | Documentation states this is a limited feature subset and does not claim full Abaqus compatibility. | not numeric | I/O Definition Agent; Coordinator Agent | draft |
+| EB3D-REQ-015 | project acceptance gates | must | release gate review | Kernel completion can be marked separately, while full solver release remains blocked until parser, assembly, HDF5 output, reference comparison, physics sanity, and release readiness gates pass. | not numeric | Release Agent; Coordinator Agent | draft |
+
+## Acceptance Criteria
+### This Documentation Step
+- `docs/requirements/euler-beam-3d.md` exists and contains metadata, assumptions, non-goals, must requirements, verification quantities, acceptance criteria, and open issues.
+- The document does not create source, test, or reference artifact files.
+- The document does not claim full Abaqus compatibility.
+- Every `must` requirement has an acceptance criterion or downstream verification hook in the Requirement Verification Matrix.
+- The following validation commands pass:
+
+```powershell
+python -m unittest discover -s scripts -p "test_*.py"
+python scripts/validate_workspace.py
+```
+
+### Kernel Completion
+Kernel completion may be claimed only after downstream implementation steps provide C++ tests and implementation evidence for:
+- local 12x12 stiffness matrix generation;
+- global 12x12 stiffness matrix transformation;
+- local and global end-force recovery;
+- DOF ordering;
+- stiffness symmetry;
+- six rigid body modes;
+- local/global transformation consistency.
+
+Kernel completion does not require Abaqus reference CSV files, `results.h5`, parser integration, global sparse assembly, or release readiness.
+
+### Full Solver Release Readiness
+Full solver release readiness is not in scope for this increment. It requires later evidence for:
+- parser and I/O contract approval for the beam feature;
+- Domain, AnalysisModel, DofManager, assembly, constraints, solve, and ResultsWriter integration;
+- HDF5 result schema coverage;
+- approved reference artifacts under `reference/<model-id>/`;
+- FESA `results.h5` to Abaqus reference CSV comparison;
+- physics sanity review;
+- release readiness review with known limitations.
+
+## Open Issues
+- Parser integration: define the approved Abaqus keyword subset, beam section mapping, orientation input, and unsupported keyword diagnostics before parser work starts.
+- Reference artifact availability: no Abaqus reference CSV files are required or generated in this kernel increment; approved reference artifacts are still needed for later solver-level verification.
+- Full end-to-end assembly: global sparse assembly, boundary condition application, reaction recovery, solve flow, and HDF5 output are not covered by this kernel baseline.
+- Local orientation convention: downstream formulation and I/O contract must define how local `y` and `z` axes are supplied, derived, validated, and reported.
+- Tolerance values: downstream numerical review must select tolerances for matrix symmetry, rigid body mode residuals, transformation consistency, and force recovery.
+
+## Downstream Handoff
+### Research Agent
+- Confirm authoritative Euler-Bernoulli 3D beam references and benchmark cases for stiffness, transformation, and rigid body checks.
+
+### Formulation Agent
+- Define local stiffness, transformation convention, end-force recovery convention, rigid body mode checks, and edge cases for invalid geometry or constants.
+
+### Numerical Review Agent
+- Review stiffness symmetry, rank/rigid body modes, sign conventions, transformation convention, conditioning risks, and numerical tolerances.
+
+### I/O Definition Agent
+- Define future parser and semantic model mapping for beam geometry, section constants, material constants, and local orientation without claiming full Abaqus compatibility.
+
+### Reference Model Agent
+- Define future reference model bundles only after parser and end-to-end solver scope is approved; do not generate reference CSV files in this increment.
+
+### Implementation Planning Agent
+- Plan C++ tests before production code for local stiffness, global transformation, end-force recovery, DOF order, symmetry, rigid body modes, and invalid inputs.

docs/research/euler-beam-3d-research.md

@@ -0,0 +1,81 @@
+# 3D Euler-Bernoulli Beam Kernel Research Brief
+
+## Metadata
+- feature_id: euler-beam-3d
+- source_requirement: docs/requirements/euler-beam-3d.md
+- status: ready-for-formulation
+- owner_agent: research-agent
+- date: 2026-06-12
+
+## Research Questions
+- What theory is sufficient for a straight, prismatic, two-node 3D Euler-Bernoulli beam kernel?
+- Which terms must appear in the local stiffness matrix for axial, torsion, and two uncoupled bending planes?
+- Which implementation checks can verify the kernel without running Abaqus or creating reference CSV artifacts?
+- Which input cases must be rejected before matrix construction?
+
+## Source Inventory
+
+| source_type | title | author_or_org | URL_or_identifier | access_date | reliability_tier | notes |
+| --- | --- | --- | --- | --- | --- | --- |
+| textbook | Concepts and Applications of Finite Element Analysis | Cook, Malkus, Plesha, Witt | ISBN 978-0-471-35605-9 | 2026-06-12 | Tier 2 | Beam element stiffness and coordinate transformation background. |
+| textbook | A First Course in the Finite Element Method | Daryl L. Logan | ISBN 978-1-305-63511-1 | 2026-06-12 | Tier 2 | Direct-stiffness beam and frame examples. |
+| textbook | Matrix Structural Analysis | William McGuire, Richard H. Gallagher, Ronald D. Ziemian | ISBN 978-0-471-12379-7 | 2026-06-12 | Tier 2 | Space-frame element local stiffness and transformation conventions. |
+| project requirement | 3D Euler-Bernoulli Beam Kernel Requirements | FESA | docs/requirements/euler-beam-3d.md | 2026-06-12 | Tier 1 project contract | Defines approved kernel scope and exclusions. |
+
+## Extracted Facts
+- Verified fact: a 3D straight prismatic Euler-Bernoulli beam frame element can be assembled from one axial response, one torsional response, and two uncoupled cubic-Hermite bending responses in the element local frame.
+- Verified fact: for a local element axis `x`, bending displacement in local `y` couples with rotation about local `z` and uses `EIz`; bending displacement in local `z` couples with rotation about local `y` and uses `EIy`.
+- Verified fact: a two-node space-frame element with six DOFs per node has a 12-entry vector. The FESA requirement fixes per-node order as `ux, uy, uz, rx, ry, rz`.
+- Verified fact: the unconstrained element local stiffness is symmetric positive semidefinite and has six rigid body modes before constraints are applied.
+- Project contract: no shear deformation, warping, end releases, offsets, mass, geometric stiffness, nonlinear kinematics, parser integration, HDF5 output, or reference artifact generation is allowed in this kernel increment.
+
+## Candidate Benchmarks
+
+| benchmark_id | source | benchmark_type | physics | target_quantities | artifact_needs | applicability |
+| --- | --- | --- | --- | --- | --- | --- |
+| EB3D-BENCH-001 | project/formulation | code verification | local matrix algebra | symmetry of all `K(i,j)` | C++ unit test only | Does not verify solver assembly. |
+| EB3D-BENCH-002 | direct stiffness theory | analytical | axial extension | `EA/L` terms and equal/opposite axial end forces | C++ unit test only | Axis-aligned local response. |
+| EB3D-BENCH-003 | direct stiffness theory | analytical | torsion | `GJ/L` terms and equal/opposite torsional end moments | C++ unit test only | Saint-Venant torsion constant assumed supplied. |
+| EB3D-BENCH-004 | beam theory | analytical | bending in local `x-y` | `12EIz/L^3`, `6EIz/L^2`, `4EIz/L`, `2EIz/L` terms | C++ unit test only | Euler-Bernoulli small-displacement bending. |
+| EB3D-BENCH-005 | beam theory | analytical | bending in local `x-z` | `12EIy/L^3`, `6EIy/L^2`, `4EIy/L`, `2EIy/L` terms | C++ unit test only | Euler-Bernoulli small-displacement bending. |
+| EB3D-BENCH-006 | matrix structural analysis | invariant | rigid body motion | local end forces are near zero for rigid translations/rotations | C++ unit test only | Numerical tolerance needed. |
+| EB3D-BENCH-007 | matrix structural analysis | invariant | coordinate transform | axis-aligned global stiffness equals local stiffness when basis is identity | C++ unit test only | Does not validate parser orientation input. |
+
+## Verification Relevance
+- code_verification: local stiffness entries, symmetry, end-force recovery as `K*u`, invalid input handling, and transformation invariants can be checked by deterministic C++ unit tests.
+- solution_verification: future cantilever axial, torsion, and bending reference models can compare displacements, reactions, and internal forces after parser/assembly/HDF5 integration exists.
+- validation: no physical experiment validation is in scope for this kernel increment.
+- reference_comparison: future reference comparison requires Abaqus-generated artifacts under `reference/<model-id>/`; this phase only defines the contract and must not create those files.
+
+## Applicability Limits
+- linear_or_nonlinear: linear only.
+- deformation: small displacement and small rotation.
+- element_type: straight two-node prismatic Euler-Bernoulli beam.
+- material_model: constant linear elastic `E` and `G`.
+- geometry: nonzero length, no offsets, no curved beams.
+- boundary_conditions: not applied inside the kernel.
+- loads: no distributed or equivalent nodal load vector in this increment.
+- coordinate_system: right-handed local Cartesian basis and global Cartesian basis.
+- units: user-consistent; no conversion.
+
+## Risks and Open Issues
+- Orientation vector parallel or nearly parallel to the beam axis must be rejected because the local `y` direction is undefined.
+- Zero or near-zero length must be rejected before stiffness coefficient computation.
+- Nonpositive or nonfinite `E`, `G`, `A`, `J`, `Iy`, or `Iz` must be rejected.
+- Very slender beams can produce ill-conditioned local stiffness matrices; kernel tests should avoid relying on a condition-number estimate as a pass/fail criterion.
+- Exact parser keyword mapping for future beam input remains open and belongs to the I/O contract, not the kernel implementation.
+
+## Downstream Handoff
+
+### Formulation Agent
+- Define the local stiffness matrix with named axial, torsion, `EIz`, and `EIy` coefficients.
+- Define the local-to-global transform convention so global stiffness and force recovery are unambiguous.
+
+### Numerical Review Agent
+- Review symmetry, rigid body modes, sign convention, coefficient dimensions, and invalid geometry handling.
+
+### Reference Model Agent
+- Future reference models should include axial cantilever, torsion cantilever, two bending cantilevers, and one skew transform case.
+
+### Implementation Planning Agent
+- Unit tests should cover representative local entries, symmetry, `K*u` force recovery, invalid constants, axis-aligned transform identity, rotated symmetry, rigid translation, and axial end forces.

phases/euler-beam-3d/index.json

@@ -0,0 +1,113 @@
+{
+  "project": "FESA Structural Solver",
+  "phase": "euler-beam-3d",
+  "steps": [
+    {
+      "step": 0,
+      "name": "requirements-baseline",
+      "status": "completed",
+      "summary": "3D Euler beam kernel requirements baseline added",
+      "allowed_paths": [
+        "docs/requirements/euler-beam-3d.md"
+      ],
+      "started_at": "2026-06-12T17:44:12+0900",
+      "completed_at": "2026-06-12T17:50:31+0900"
+    },
+    {
+      "step": 1,
+      "name": "research-evidence",
+      "status": "completed",
+      "summary": "3D Euler beam research evidence added",
+      "allowed_paths": [
+        "docs/research/euler-beam-3d-research.md"
+      ]
+    },
+    {
+      "step": 2,
+      "name": "formulation-spec",
+      "status": "completed",
+      "summary": "3D Euler beam formulation contract added",
+      "allowed_paths": [
+        "docs/formulations/euler-beam-3d-formulation.md"
+      ]
+    },
+    {
+      "step": 3,
+      "name": "numerical-review",
+      "status": "completed",
+      "summary": "3D Euler beam numerical review added",
+      "allowed_paths": [
+        "docs/numerical-reviews/euler-beam-3d-review.md"
+      ]
+    },
+    {
+      "step": 4,
+      "name": "io-reference-contract",
+      "status": "completed",
+      "summary": "3D Euler beam I/O and reference model contracts added",
+      "allowed_paths": [
+        "docs/io-definitions/euler-beam-3d-io.md",
+        "docs/reference-models/euler-beam-3d-reference-models.md"
+      ]
+    },
+    {
+      "step": 5,
+      "name": "implementation-plan",
+      "status": "completed",
+      "summary": "3D Euler beam implementation plan added",
+      "allowed_paths": [
+        "docs/implementation-plans/euler-beam-3d-implementation-plan.md"
+      ]
+    },
+    {
+      "step": 6,
+      "name": "model-beam-topology",
+      "status": "completed",
+      "summary": "beam2 model topology added with unit test",
+      "allowed_paths": [
+        "src/fesa/model/element.hpp",
+        "src/fesa/model/element.cpp",
+        "tests/unit/model_element_test.cpp"
+      ]
+    },
+    {
+      "step": 7,
+      "name": "local-stiffness-kernel",
+      "status": "completed",
+      "summary": "local 3D Euler beam stiffness and local end-force kernel added",
+      "allowed_paths": [
+        "src/fesa/elements/",
+        "tests/unit/euler_beam_3d_*_test.cpp"
+      ]
+    },
+    {
+      "step": 8,
+      "name": "global-transform-recovery",
+      "status": "completed",
+      "summary": "global transform and global end-force recovery added for 3D Euler beam",
+      "allowed_paths": [
+        "src/fesa/elements/",
+        "tests/unit/euler_beam_3d_*_test.cpp"
+      ]
+    },
+    {
+      "step": 9,
+      "name": "build-test-report",
+      "status": "completed",
+      "summary": "3D Euler beam build/test report added",
+      "allowed_paths": [
+        "docs/build-test-reports/euler-beam-3d-build-test.md"
+      ]
+    },
+    {
+      "step": 10,
+      "name": "release-readiness-note",
+      "status": "completed",
+      "summary": "3D Euler beam release readiness note added",
+      "allowed_paths": [
+        "docs/releases/euler-beam-3d-release.md"
+      ]
+    }
+  ],
+  "created_at": "2026-06-12T17:44:12+0900"
+}

phases/euler-beam-3d/step0.md

@@ -0,0 +1,63 @@
+# Step 0: requirements-baseline
+
+## Read These Files First
+
+Read the following files before editing:
+
+- `/AGENTS.md`
+- `/docs/PRD.md`
+- `/docs/ARCHITECTURE.md`
+- `/docs/ADR.md`
+- `/docs/requirements/README.md`
+
+## Task
+
+Create `/docs/requirements/euler-beam-3d.md` for a kernel-first 3D Euler-Bernoulli beam feature.
+
+Scope the first implementable increment narrowly:
+
+- two-node, straight, prismatic, small-displacement 3D Euler-Bernoulli beam element
+- six mechanical DOFs per node in this order: `ux, uy, uz, rx, ry, rz`
+- linear elastic section constants: `E`, `G`, `A`, `J`, `Iy`, `Iz`
+- local and global 12x12 stiffness matrix support
+- local and global element end-force recovery from nodal displacement vectors
+- no shear deformation, warping, end releases, offsets, distributed loads, mass matrix, geometric stiffness, nonlinear kinematics, dynamics, or thermal coupling in this increment
+- no Abaqus reference CSV generation in this increment
+- no full Abaqus compatibility claim
+
+The document must contain:
+
+- metadata: `feature_id: euler-beam-3d`, owner agent, status
+- explicit assumptions and non-goals
+- `must` requirements with stable IDs such as `EB3D-REQ-001`
+- verification quantities: displacement DOFs, reactions/end forces, stiffness symmetry, rigid body modes, local/global transformation consistency
+- acceptance criteria that distinguish kernel completion from full solver release readiness
+- open issues for parser integration, reference artifact availability, and full end-to-end assembly
+
+Do not create C++ files in this step.
+
+## Tests To Write First
+
+- No C++ test is required in this documentation-only step.
+- Still run the harness validation commands in the acceptance criteria.
+
+## Acceptance Criteria
+
+```powershell
+python -m unittest discover -s scripts -p "test_*.py"
+python scripts/validate_workspace.py
+```
+
+## Verification Notes
+
+1. Confirm the requirements do not claim full Abaqus compatibility.
+2. Confirm each `must` requirement has an acceptance criterion or a downstream verification hook.
+3. Update `phases/euler-beam-3d/index.json` step 0:
+   - success: `"status": "completed"`, `"summary": "3D Euler beam kernel requirements baseline added"`
+   - failure after retries: `"status": "error"`, `"error_message": "<specific error>"`
+   - blocked: `"status": "blocked"`, `"blocked_reason": "<specific reason>"`
+
+## Forbidden
+
+- Do not create or modify reference CSV files.
+- Do not modify source or test files.

phases/euler-beam-3d/step1.md

@@ -0,0 +1,60 @@
+# Step 1: research-evidence
+
+## Read These Files First
+
+Read the following files before editing:
+
+- `/AGENTS.md`
+- `/docs/PRD.md`
+- `/docs/ARCHITECTURE.md`
+- `/docs/ADR.md`
+- `/docs/requirements/euler-beam-3d.md`
+- `/docs/research/README.md`
+
+## Task
+
+Create `/docs/research/euler-beam-3d-research.md`.
+
+Summarize the evidence needed to implement the approved kernel-first 3D Euler-Bernoulli beam increment. The document must be implementation-oriented and must include:
+
+- supported theory: straight prismatic Euler-Bernoulli beam, axial, torsion, and two uncoupled bending planes
+- assumptions and applicability limits
+- local DOF ordering and sign convention used by the planned matrix
+- source reliability classification
+- benchmark-style checks that do not require external reference solver execution:
+  - local stiffness symmetry
+  - axial-only response
+  - torsion-only response
+  - bending about local `y`
+  - bending about local `z`
+  - rigid body translation/rotation zero internal forces in local coordinates
+  - global transform identity for an axis-aligned beam
+- risks: orientation vector parallel to element axis, near-zero length, nonpositive section constants, ill-conditioning for very slender elements
+
+If internet access or a FEM wiki is used, cite sources briefly. Do not include long copyrighted excerpts.
+
+Do not create C++ files in this step.
+
+## Tests To Write First
+
+- No C++ test is required in this documentation-only step.
+
+## Acceptance Criteria
+
+```powershell
+python -m unittest discover -s scripts -p "test_*.py"
+python scripts/validate_workspace.py
+```
+
+## Verification Notes
+
+1. Confirm the research document ties each source or theory point to a planned implementation check.
+2. Confirm unresolved items are listed as risks or open issues instead of silently assumed.
+3. Update `phases/euler-beam-3d/index.json` step 1:
+   - success: `"status": "completed"`, `"summary": "3D Euler beam research evidence added"`
+   - failure after retries: `"status": "error"`, `"error_message": "<specific error>"`
+   - blocked: `"status": "blocked"`, `"blocked_reason": "<specific reason>"`
+
+## Forbidden
+
+- Do not modify source, test, reference, or I/O contract files.

phases/euler-beam-3d/step10.md

@@ -0,0 +1,60 @@
+# Step 10: release-readiness-note
+
+## Read These Files First
+
+Read the following files before editing:
+
+- `/AGENTS.md`
+- `/docs/PRD.md`
+- `/docs/ADR.md`
+- `/docs/releases/README.md`
+- `/docs/requirements/euler-beam-3d.md`
+- `/docs/io-definitions/euler-beam-3d-io.md`
+- `/docs/reference-models/euler-beam-3d-reference-models.md`
+- `/docs/build-test-reports/euler-beam-3d-build-test.md`
+
+## Task
+
+Create `/docs/releases/euler-beam-3d-release.md`.
+
+This is a readiness note, not a release approval. It must state:
+
+- status: `not-release-ready-kernel-increment-complete` unless all upstream gates and reference artifacts somehow exist
+- completed scope: local/global stiffness and end-force kernel for two-node 3D Euler-Bernoulli beam
+- missing for full feature release:
+  - parser implementation for the approved Abaqus subset
+  - section/property semantic model integration
+  - assembler/static solver integration
+  - HDF5 result emission for beam quantities
+  - stored Abaqus reference artifacts
+  - reference comparison report
+  - physics sanity report
+- known limitations from requirements
+- next recommended phase or phase dependencies
+
+Do not change source or tests in this step.
+
+## Tests To Write First
+
+- No C++ test is required in this documentation-only step.
+
+## Acceptance Criteria
+
+```powershell
+python -m unittest discover -s scripts -p "test_*.py"
+python scripts/validate_workspace.py
+```
+
+## Verification Notes
+
+1. Confirm the note does not claim full release readiness.
+2. Confirm all missing gates are explicit.
+3. Update `phases/euler-beam-3d/index.json` step 10:
+   - success: `"status": "completed"`, `"summary": "3D Euler beam release readiness note added"`
+   - failure after retries: `"status": "error"`, `"error_message": "<specific error>"`
+   - blocked: `"status": "blocked"`, `"blocked_reason": "<specific reason>"`
+
+## Forbidden
+
+- Do not create or modify reference artifacts.
+- Do not modify source, tests, requirements, formulations, or I/O contracts.

phases/euler-beam-3d/step2.md

@@ -0,0 +1,63 @@
+# Step 2: formulation-spec
+
+## Read These Files First
+
+Read the following files before editing:
+
+- `/AGENTS.md`
+- `/docs/PRD.md`
+- `/docs/ARCHITECTURE.md`
+- `/docs/ADR.md`
+- `/docs/requirements/euler-beam-3d.md`
+- `/docs/research/euler-beam-3d-research.md`
+- `/docs/formulations/README.md`
+
+## Task
+
+Create `/docs/formulations/euler-beam-3d-formulation.md`.
+
+The formulation must define the implementation contract for the C++ kernel. Include:
+
+- local coordinate system: local `x` from node 1 to node 2; local `y` from the user orientation vector projected normal to local `x`; local `z = x cross y`
+- local DOF order: `[u1, v1, w1, rx1, ry1, rz1, u2, v2, w2, rx2, ry2, rz2]`
+- section constants: `E`, `G`, `A`, `J`, `Iy`, `Iz`
+- local 12x12 stiffness matrix terms for:
+  - axial: `EA/L`
+  - torsion: `GJ/L`
+  - bending in local `x-y` using `EIz`
+  - bending in local `x-z` using `EIy`
+- transformation matrix convention and global stiffness equation `K_global = T^T K_local T`
+- end-force recovery equation `f_local = K_local u_local`
+- validation tolerances for unit tests, using deterministic double comparisons
+- singular and invalid input handling:
+  - zero or near-zero length throws `std::invalid_argument`
+  - nonpositive `E`, `G`, `A`, `J`, `Iy`, or `Iz` throws `std::invalid_argument`
+  - orientation vector parallel to beam axis throws `std::invalid_argument`
+
+Use compact matrices and named scalar coefficients so the implementation step can transcribe directly.
+
+Do not create C++ files in this step.
+
+## Tests To Write First
+
+- No C++ test is required in this documentation-only step.
+
+## Acceptance Criteria
+
+```powershell
+python -m unittest discover -s scripts -p "test_*.py"
+python scripts/validate_workspace.py
+```
+
+## Verification Notes
+
+1. Confirm the matrix is symmetric by construction.
+2. Confirm local/global transformation convention is unambiguous.
+3. Update `phases/euler-beam-3d/index.json` step 2:
+   - success: `"status": "completed"`, `"summary": "3D Euler beam formulation contract added"`
+   - failure after retries: `"status": "error"`, `"error_message": "<specific error>"`
+   - blocked: `"status": "blocked"`, `"blocked_reason": "<specific reason>"`
+
+## Forbidden
+
+- Do not modify requirements, source, tests, or reference artifacts.

phases/euler-beam-3d/step3.md

@@ -0,0 +1,56 @@
+# Step 3: numerical-review
+
+## Read These Files First
+
+Read the following files before editing:
+
+- `/AGENTS.md`
+- `/docs/PRD.md`
+- `/docs/ARCHITECTURE.md`
+- `/docs/ADR.md`
+- `/docs/requirements/euler-beam-3d.md`
+- `/docs/research/euler-beam-3d-research.md`
+- `/docs/formulations/euler-beam-3d-formulation.md`
+- `/docs/numerical-reviews/README.md`
+
+## Task
+
+Create `/docs/numerical-reviews/euler-beam-3d-review.md`.
+
+Review the formulation for implementation readiness. The review must include:
+
+- dimension and units check for each stiffness coefficient
+- symmetry and rigid body mode expectations
+- local axis construction risks
+- positive semi-definite local stiffness expectation before constraints
+- test obligations before production code
+- known numerical limits for slender beams and very small section constants
+- pass/fail verdict for kernel implementation planning
+
+If the formulation is not ready, mark the document status as `needs-upstream-decision` and update this phase step as blocked with a concrete reason.
+
+Do not create C++ files in this step.
+
+## Tests To Write First
+
+- No C++ test is required in this documentation-only step.
+
+## Acceptance Criteria
+
+```powershell
+python -m unittest discover -s scripts -p "test_*.py"
+python scripts/validate_workspace.py
+```
+
+## Verification Notes
+
+1. Confirm the review does not change requirements or formulation.
+2. Confirm implementation-owned risks are turned into concrete tests.
+3. Update `phases/euler-beam-3d/index.json` step 3:
+   - success: `"status": "completed"`, `"summary": "3D Euler beam numerical review added"`
+   - failure after retries: `"status": "error"`, `"error_message": "<specific error>"`
+   - blocked: `"status": "blocked"`, `"blocked_reason": "<specific reason>"`
+
+## Forbidden
+
+- Do not modify requirements, formulation, source, tests, or reference artifacts.

phases/euler-beam-3d/step4.md

@@ -0,0 +1,68 @@
+# Step 4: io-reference-contract
+
+## Read These Files First
+
+Read the following files before editing:
+
+- `/AGENTS.md`
+- `/docs/PRD.md`
+- `/docs/ARCHITECTURE.md`
+- `/docs/ADR.md`
+- `/docs/requirements/euler-beam-3d.md`
+- `/docs/formulations/euler-beam-3d-formulation.md`
+- `/docs/numerical-reviews/euler-beam-3d-review.md`
+- `/docs/io-definitions/README.md`
+- `/docs/reference-models/README.md`
+
+## Task
+
+Create two documents:
+
+- `/docs/io-definitions/euler-beam-3d-io.md`
+- `/docs/reference-models/euler-beam-3d-reference-models.md`
+
+The I/O contract must describe the planned Abaqus subset without claiming it is implemented in this kernel increment:
+
+- element keyword mapping candidate: two-node beam topology equivalent to Abaqus `B31`
+- section keyword candidate: beam section constants sufficient for `A`, `J`, `Iy`, `Iz`, `E`, and `G`
+- orientation data requirement for constructing local axes
+- boundary/load DOFs: `1..6` map to `ux, uy, uz, rx, ry, rz`
+- HDF5 output quantities expected after solver integration: nodal displacement, reaction, element internal force, stress placeholders if stress recovery is later approved
+- unsupported input cases and diagnostics
+
+The reference model contract must list required future models without generating artifacts:
+
+- axial cantilever bar-as-beam
+- torsion cantilever
+- bending cantilever about local `y`
+- bending cantilever about local `z`
+- skew-oriented beam transform check
+
+For each future model, specify expected files under `reference/<model-id>/` and which CSV quantities are required. State that Abaqus reference CSVs must not be generated or modified in this phase.
+
+Do not create C++ files in this step.
+
+## Tests To Write First
+
+- No C++ test is required in this documentation-only step.
+
+## Acceptance Criteria
+
+```powershell
+python -m unittest discover -s scripts -p "test_*.py"
+python scripts/validate_workspace.py
+```
+
+## Verification Notes
+
+1. Confirm reference artifacts are specified, not created.
+2. Confirm unsupported parser or solver paths are explicit open issues.
+3. Update `phases/euler-beam-3d/index.json` step 4:
+   - success: `"status": "completed"`, `"summary": "3D Euler beam I/O and reference model contracts added"`
+   - failure after retries: `"status": "error"`, `"error_message": "<specific error>"`
+   - blocked: `"status": "blocked"`, `"blocked_reason": "<specific reason>"`
+
+## Forbidden
+
+- Do not create or modify files under `/reference/`.
+- Do not modify source or tests.

phases/euler-beam-3d/step5.md

@@ -0,0 +1,75 @@
+# Step 5: implementation-plan
+
+## Read These Files First
+
+Read the following files before editing:
+
+- `/AGENTS.md`
+- `/docs/PRD.md`
+- `/docs/ARCHITECTURE.md`
+- `/docs/ADR.md`
+- `/docs/SOLVER_AGENT_DESIGN.md`
+- `/docs/implementation-plans/README.md`
+- `/docs/requirements/euler-beam-3d.md`
+- `/docs/research/euler-beam-3d-research.md`
+- `/docs/formulations/euler-beam-3d-formulation.md`
+- `/docs/numerical-reviews/euler-beam-3d-review.md`
+- `/docs/io-definitions/euler-beam-3d-io.md`
+- `/docs/reference-models/euler-beam-3d-reference-models.md`
+- `/src/fesa/model/element.hpp`
+- `/tests/unit/model_element_test.cpp`
+
+## Task
+
+Create `/docs/implementation-plans/euler-beam-3d-implementation-plan.md`.
+
+The implementation plan must be ready for C++ TDD work and must include:
+
+- readiness check for requirements, research, formulation, numerical review, I/O, reference model contract
+- explicit kernel-first implementation scope
+- tasks matching the remaining phase steps:
+  - model beam topology
+  - local stiffness kernel
+  - global transform and end-force recovery
+  - build/test report
+  - release readiness note
+- test IDs and RED/GREEN conditions for each production change
+- candidate files:
+  - `src/fesa/model/element.hpp`
+  - `src/fesa/model/element.cpp`
+  - `tests/unit/model_element_test.cpp`
+  - `src/fesa/elements/euler_beam_3d.hpp`
+  - `src/fesa/elements/euler_beam_3d.cpp`
+  - `tests/unit/euler_beam_3d_local_stiffness_test.cpp`
+  - `tests/unit/euler_beam_3d_transform_recovery_test.cpp`
+- CTest commands:
+  - `ctest --test-dir build/msvc-debug --output-on-failure -C Debug -R model_element_test`
+  - `ctest --test-dir build/msvc-debug --output-on-failure -C Debug -R euler_beam_3d_local_stiffness_test`
+  - `ctest --test-dir build/msvc-debug --output-on-failure -C Debug -R euler_beam_3d_transform_recovery_test`
+- acceptance traceability from requirements to tests
+
+Do not create C++ files in this step.
+
+## Tests To Write First
+
+- No C++ test is required in this documentation-only step.
+
+## Acceptance Criteria
+
+```powershell
+python -m unittest discover -s scripts -p "test_*.py"
+python scripts/validate_workspace.py
+```
+
+## Verification Notes
+
+1. Confirm every production file named in the plan has a related test.
+2. Confirm the plan does not ask Implementation Agent to change requirements, formulation, I/O contracts, reference artifacts, or tolerance policies.
+3. Update `phases/euler-beam-3d/index.json` step 5:
+   - success: `"status": "completed"`, `"summary": "3D Euler beam implementation plan added"`
+   - failure after retries: `"status": "error"`, `"error_message": "<specific error>"`
+   - blocked: `"status": "blocked"`, `"blocked_reason": "<specific reason>"`
+
+## Forbidden
+
+- Do not modify source, tests, reference artifacts, or CMake files.

phases/euler-beam-3d/step6.md

@@ -0,0 +1,79 @@
+# Step 6: model-beam-topology
+
+## Read These Files First
+
+Read the following files before editing:
+
+- `/AGENTS.md`
+- `/docs/ARCHITECTURE.md`
+- `/docs/ADR.md`
+- `/docs/implementation-plans/euler-beam-3d-implementation-plan.md`
+- `/src/fesa/model/element.hpp`
+- `/src/fesa/model/element.cpp`
+- `/tests/unit/model_element_test.cpp`
+
+## Task
+
+Use TDD to add a semantic model topology for a two-node 3D beam.
+
+Required production behavior:
+
+- Add `beam2` to `fesa::model::ElementTopology`.
+- Preserve existing `truss2`, `bar2`, and `unknown` behavior.
+- Do not store equation IDs on `Element`.
+- Do not add section constants to `Element` in this step.
+
+## Tests To Write First
+
+Modify `/tests/unit/model_element_test.cpp` first.
+
+Add a failing assertion that constructs an element with:
+
+- id `ElementId{10}`
+- topology `ElementTopology::beam2`
+- node ids `{NodeId{1}, NodeId{2}}`
+- property id `PropertyId{7}`
+
+The test must verify:
+
+- `element.topology() == ElementTopology::beam2`
+- `element.node_ids().size() == 2`
+- existing `bar2` behavior still works or the test still covers it
+
+RED command:
+
+```powershell
+ctest --test-dir build/msvc-debug --output-on-failure -C Debug -R model_element_test
+```
+
+Expected RED: compile failure because `ElementTopology::beam2` is not defined.
+
+Then implement the minimal enum addition.
+
+GREEN command:
+
+```powershell
+ctest --test-dir build/msvc-debug --output-on-failure -C Debug -R model_element_test
+```
+
+## Acceptance Criteria
+
+```powershell
+ctest --test-dir build/msvc-debug --output-on-failure -C Debug -R model_element_test
+python -m unittest discover -s scripts -p "test_*.py"
+python scripts/validate_workspace.py
+```
+
+## Verification Notes
+
+1. Confirm the test failed before editing production code.
+2. Confirm no unrelated model refactor was made.
+3. Update `phases/euler-beam-3d/index.json` step 6:
+   - success: `"status": "completed"`, `"summary": "beam2 model topology added with unit test"`
+   - failure after retries: `"status": "error"`, `"error_message": "<specific error>"`
+   - blocked: `"status": "blocked"`, `"blocked_reason": "<specific reason>"`
+
+## Forbidden
+
+- Do not create parser code.
+- Do not add beam stiffness code in this step.

phases/euler-beam-3d/step7.md

@@ -0,0 +1,110 @@
+# Step 7: local-stiffness-kernel
+
+## Read These Files First
+
+Read the following files before editing:
+
+- `/AGENTS.md`
+- `/docs/ARCHITECTURE.md`
+- `/docs/ADR.md`
+- `/docs/formulations/euler-beam-3d-formulation.md`
+- `/docs/numerical-reviews/euler-beam-3d-review.md`
+- `/docs/implementation-plans/euler-beam-3d-implementation-plan.md`
+
+Also read any files created by previous steps in this phase.
+
+## Task
+
+Use TDD to create the local 12x12 stiffness kernel for the 3D Euler-Bernoulli beam.
+
+Create:
+
+- `/src/fesa/elements/euler_beam_3d.hpp`
+- `/src/fesa/elements/euler_beam_3d.cpp`
+- `/tests/unit/euler_beam_3d_local_stiffness_test.cpp`
+
+Required API:
+
+```cpp
+namespace fesa::elements {
+
+using Vector12 = std::array<double, 12>;
+using Matrix12 = std::array<double, 144>;
+
+struct EulerBeam3DSection {
+    double young_modulus;
+    double shear_modulus;
+    double area;
+    double torsion_constant;
+    double second_moment_y;
+    double second_moment_z;
+};
+
+Matrix12 euler_beam_3d_local_stiffness(double length, const EulerBeam3DSection& section);
+Vector12 euler_beam_3d_local_end_forces(double length,
+                                        const EulerBeam3DSection& section,
+                                        const Vector12& local_displacements);
+
+} // namespace fesa::elements
+```
+
+Implementation rules:
+
+- Matrix storage is row-major: index `(row, col)` is `row * 12 + col`.
+- Validate `length > 0` and all section constants are positive; throw `std::invalid_argument` otherwise.
+- Use only C++17 standard library.
+- Do not introduce an external linear algebra dependency.
+- Do not edit CMake files; source and test globs should pick up new files.
+
+## Tests To Write First
+
+Create `/tests/unit/euler_beam_3d_local_stiffness_test.cpp` before production code.
+
+Test behavior:
+
+- For `L = 2.0`, `E = 210.0`, `G = 80.0`, `A = 3.0`, `J = 4.0`, `Iy = 5.0`, `Iz = 6.0`, verify representative matrix entries:
+  - axial: `K(0,0) = EA/L`, `K(0,6) = -EA/L`, `K(6,6) = EA/L`
+  - torsion: `K(3,3) = GJ/L`, `K(3,9) = -GJ/L`, `K(9,9) = GJ/L`
+  - local `x-y` bending uses `EIz`: `K(1,1) = 12*E*Iz/L^3`, `K(1,5) = 6*E*Iz/L^2`, `K(5,5) = 4*E*Iz/L`
+  - local `x-z` bending uses `EIy`: `K(2,2) = 12*E*Iy/L^3`, `K(2,4) = -6*E*Iy/L^2`, `K(4,4) = 4*E*Iy/L`
+- Verify all entries are symmetric within `1.0e-10`.
+- Verify `euler_beam_3d_local_end_forces` returns `K * u` for a displacement vector with at least three nonzero components.
+- Verify invalid length and nonpositive section constants throw `std::invalid_argument`.
+
+RED command:
+
+```powershell
+ctest --test-dir build/msvc-debug --output-on-failure -C Debug -R euler_beam_3d_local_stiffness_test
+```
+
+Expected RED: test executable missing or compile failure because the header/API does not exist.
+
+Then implement the minimal API.
+
+GREEN command:
+
+```powershell
+ctest --test-dir build/msvc-debug --output-on-failure -C Debug -R euler_beam_3d_local_stiffness_test
+```
+
+## Acceptance Criteria
+
+```powershell
+ctest --test-dir build/msvc-debug --output-on-failure -C Debug -R euler_beam_3d_local_stiffness_test
+python -m unittest discover -s scripts -p "test_*.py"
+python scripts/validate_workspace.py
+```
+
+## Verification Notes
+
+1. Confirm RED failed before production code was added.
+2. Confirm no parser, assembly, or results writer code was changed.
+3. Update `phases/euler-beam-3d/index.json` step 7:
+   - success: `"status": "completed"`, `"summary": "local 3D Euler beam stiffness and local end-force kernel added"`
+   - failure after retries: `"status": "error"`, `"error_message": "<specific error>"`
+   - blocked: `"status": "blocked"`, `"blocked_reason": "<specific reason>"`
+
+## Forbidden
+
+- Do not add shear deformation or Timoshenko terms.
+- Do not modify reference artifacts.

phases/euler-beam-3d/step8.md

@@ -0,0 +1,97 @@
+# Step 8: global-transform-recovery
+
+## Read These Files First
+
+Read the following files before editing:
+
+- `/AGENTS.md`
+- `/docs/ARCHITECTURE.md`
+- `/docs/ADR.md`
+- `/docs/formulations/euler-beam-3d-formulation.md`
+- `/docs/numerical-reviews/euler-beam-3d-review.md`
+- `/docs/implementation-plans/euler-beam-3d-implementation-plan.md`
+- `/src/fesa/elements/euler_beam_3d.hpp`
+- `/src/fesa/elements/euler_beam_3d.cpp`
+- `/tests/unit/euler_beam_3d_local_stiffness_test.cpp`
+
+## Task
+
+Use TDD to add local/global transformation and global end-force recovery to the 3D Euler-Bernoulli beam kernel.
+
+Extend the API in `/src/fesa/elements/euler_beam_3d.hpp`:
+
+```cpp
+using Vector3 = std::array<double, 3>;
+
+struct EulerBeam3DGeometry {
+    Vector3 node1;
+    Vector3 node2;
+    Vector3 orientation;
+};
+
+Matrix12 euler_beam_3d_global_stiffness(const EulerBeam3DGeometry& geometry,
+                                        const EulerBeam3DSection& section);
+Vector12 euler_beam_3d_global_end_forces(const EulerBeam3DGeometry& geometry,
+                                         const EulerBeam3DSection& section,
+                                         const Vector12& global_displacements);
+```
+
+Implementation rules:
+
+- local `x` is normalized `node2 - node1`
+- local `y` is the normalized projection of `orientation` onto the plane normal to local `x`
+- local `z = x cross y`
+- use the same 3x3 rotation block for translational and rotational DOFs
+- compute `K_global = T^T K_local T`
+- compute global end forces by transforming global displacements to local, recovering local forces, then transforming forces back to global
+- throw `std::invalid_argument` for zero-length element, zero orientation vector, or orientation parallel to the beam axis
+
+## Tests To Write First
+
+Create `/tests/unit/euler_beam_3d_transform_recovery_test.cpp` before production code.
+
+Test behavior:
+
+- Axis-aligned beam from `(0,0,0)` to `(2,0,0)` with orientation `(0,1,0)` gives global stiffness equal to local stiffness.
+- A rotated beam preserves stiffness symmetry.
+- A rigid global translation vector produces near-zero global end forces.
+- A simple global axial extension on the axis-aligned beam produces equal and opposite axial end forces.
+- Parallel orientation vector throws `std::invalid_argument`.
+
+RED command:
+
+```powershell
+ctest --test-dir build/msvc-debug --output-on-failure -C Debug -R euler_beam_3d_transform_recovery_test
+```
+
+Expected RED: test executable missing or compile failure because the transform API does not exist.
+
+Then implement the minimal API.
+
+GREEN command:
+
+```powershell
+ctest --test-dir build/msvc-debug --output-on-failure -C Debug -R euler_beam_3d_transform_recovery_test
+```
+
+## Acceptance Criteria
+
+```powershell
+ctest --test-dir build/msvc-debug --output-on-failure -C Debug -R "euler_beam_3d_(local_stiffness|transform_recovery)_test"
+python -m unittest discover -s scripts -p "test_*.py"
+python scripts/validate_workspace.py
+```
+
+## Verification Notes
+
+1. Confirm RED failed before production code was added.
+2. Confirm local stiffness tests still pass.
+3. Update `phases/euler-beam-3d/index.json` step 8:
+   - success: `"status": "completed"`, `"summary": "global transform and global end-force recovery added for 3D Euler beam"`
+   - failure after retries: `"status": "error"`, `"error_message": "<specific error>"`
+   - blocked: `"status": "blocked"`, `"blocked_reason": "<specific reason>"`
+
+## Forbidden
+
+- Do not add assembly or solver integration in this step.
+- Do not modify reference artifacts.

phases/euler-beam-3d/step9.md

@@ -0,0 +1,57 @@
+# Step 9: build-test-report
+
+## Read These Files First
+
+Read the following files before editing:
+
+- `/AGENTS.md`
+- `/docs/build-test-reports/README.md`
+- `/docs/implementation-plans/euler-beam-3d-implementation-plan.md`
+- `/src/fesa/elements/euler_beam_3d.hpp`
+- `/src/fesa/elements/euler_beam_3d.cpp`
+- `/tests/unit/euler_beam_3d_local_stiffness_test.cpp`
+- `/tests/unit/euler_beam_3d_transform_recovery_test.cpp`
+
+## Task
+
+Create `/docs/build-test-reports/euler-beam-3d-build-test.md`.
+
+The report must record:
+
+- feature id
+- changed files observed for this phase
+- command log summary with exit codes and short evidence tails
+- validation results for:
+  - harness self-test
+  - CMake configure/build
+  - CTest
+  - feature-specific tests
+- failure classification if anything failed
+- explicit statement that this report does not approve reference verification or release readiness
+
+Do not change source or tests in this step.
+
+## Tests To Write First
+
+- No C++ test is required in this documentation-only step.
+
+## Acceptance Criteria
+
+```powershell
+ctest --test-dir build/msvc-debug --output-on-failure -C Debug -R "model_element_test|euler_beam_3d_(local_stiffness|transform_recovery)_test"
+python -m unittest discover -s scripts -p "test_*.py"
+python scripts/validate_workspace.py
+```
+
+## Verification Notes
+
+1. Record actual command exit codes and concise output evidence.
+2. If validation fails for an environment reason, mark the report `needs-environment-fix` and update this phase step as blocked.
+3. Update `phases/euler-beam-3d/index.json` step 9:
+   - success: `"status": "completed"`, `"summary": "3D Euler beam build/test report added"`
+   - failure after retries: `"status": "error"`, `"error_message": "<specific error>"`
+   - blocked: `"status": "blocked"`, `"blocked_reason": "<specific reason>"`
+
+## Forbidden
+
+- Do not modify source, tests, requirements, formulations, I/O contracts, or reference artifacts.

phases/index.json

@@ -0,0 +1,12 @@
+{
+  "phases": [
+    {
+      "dir": "solver-core-skeleton",
+      "status": "completed"
+    },
+    {
+      "dir": "euler-beam-3d",
+      "status": "pending"
+    }
+  ]
+}

phases/solver-core-skeleton/index.json

@@ -0,0 +1,100 @@
+{
+  "project": "FESA Structural Solver",
+  "phase": "solver-core-skeleton",
+  "steps": [
+    {
+      "step": 0,
+      "name": "cmake-ctest-bootstrap",
+      "status": "pending",
+      "allowed_paths": [
+        "CMakeLists.txt",
+        "tests/"
+      ],
+      "started_at": "2026-06-12T02:09:10+0900",
+      "summary": "CMake/CTest bootstrap with fesa_core interface target and smoke test",
+      "status": "completed"
+    },
+    {
+      "step": 1,
+      "name": "core-primitives",
+      "status": "completed",
+      "allowed_paths": [
+        "src/fesa/core/",
+        "tests/unit/core_*_test.cpp"
+      ],
+      "summary": "Core ID, diagnostic, and status primitives added with tests"
+    },
+    {
+      "step": 2,
+      "name": "domain-model-entities",
+      "status": "completed",
+      "allowed_paths": [
+        "src/fesa/model/",
+        "tests/unit/model_*_test.cpp"
+      ],
+      "summary": "Model entities and Domain ownership API added with tests"
+    },
+    {
+      "step": 3,
+      "name": "analysis-model-view",
+      "status": "completed",
+      "allowed_paths": [
+        "src/fesa/analysis/",
+        "tests/unit/analysis_model_*_test.cpp"
+      ],
+      "summary": "AnalysisModel step view added without Domain copies"
+    },
+    {
+      "step": 4,
+      "name": "dof-manager",
+      "status": "completed",
+      "allowed_paths": [
+        "src/fesa/fem/",
+        "tests/unit/dof_manager_*_test.cpp"
+      ],
+      "summary": "DofManager deterministic numbering and constrained/free mapping added"
+    },
+    {
+      "step": 5,
+      "name": "analysis-state",
+      "status": "completed",
+      "allowed_paths": [
+        "src/fesa/analysis/",
+        "tests/unit/analysis_state_*_test.cpp"
+      ],
+      "summary": "AnalysisState vector ownership and residual update added"
+    },
+    {
+      "step": 6,
+      "name": "analysis-template-flow",
+      "status": "completed",
+      "allowed_paths": [
+        "src/fesa/analysis/",
+        "tests/unit/analysis_flow_*_test.cpp"
+      ],
+      "summary": "Analysis template method and LinearStaticAnalysis skeleton added"
+    },
+    {
+      "step": 7,
+      "name": "results-containers",
+      "status": "completed",
+      "allowed_paths": [
+        "src/fesa/results/",
+        "tests/unit/results_*_test.cpp"
+      ],
+      "summary": "ResultStep, ResultFrame, FieldOutput, and HistoryOutput containers added"
+    },
+    {
+      "step": 8,
+      "name": "solver-skeleton-integration-report",
+      "status": "completed",
+      "allowed_paths": [
+        "tests/integration/",
+        "docs/build-test-reports/"
+      ],
+      "summary": "Solver skeleton integration test and build/test report added"
+    }
+  ],
+  "created_at": "2026-06-12T02:09:10+0900",
+  "completed_at": "2026-06-12T02:42:00+0900"
+}

phases/solver-core-skeleton/step0.md

@@ -0,0 +1,83 @@
+# Step 0: cmake-ctest-bootstrap
+
+## 읽어야 할 파일
+
+먼저 아래 파일들을 읽고 프로젝트의 아키텍처와 검증 규칙을 파악하라:
+
+- `/AGENTS.md`
+- `/docs/PRD.md`
+- `/docs/ARCHITECTURE.md`
+- `/docs/ADR.md`
+- `/scripts/validate_workspace.py`
+
+## 작업
+
+C++ solver skeleton을 구현할 수 있도록 최소 CMake/CTest 부트스트랩을 만든다.
+
+요구사항:
+
+- 루트 `/CMakeLists.txt`를 생성한다.
+- C++ 표준은 C++17 이상으로 고정한다.
+- MSVC에서 warning-as-error를 강제하지 않는다.
+- 아직 production source가 없으므로 solver target은 `INTERFACE` library `fesa_core`로 시작한다.
+- `fesa_core`는 repository root의 `src`를 include directory로 노출한다.
+- `enable_testing()`과 `tests/` 하위 CMake 구성을 연결한다.
+- `/tests/CMakeLists.txt`, `/tests/unit/CMakeLists.txt`, `/tests/integration/CMakeLists.txt`를 생성한다.
+- `tests/unit/*_test.cpp`와 `tests/integration/*_test.cpp` 파일을 각각 독립 test executable로 등록한다.
+- 각 test executable은 `fesa_core`에 link한다.
+- `/tests/unit/harness_smoke_test.cpp`를 생성하고, 표준 라이브러리만 사용해 CTest가 동작함을 확인하는 최소 `main()`을 둔다.
+- npm, JavaScript, TypeScript fallback은 추가하지 않는다.
+
+권장 CMake 구조:
+
+```cmake
+cmake_minimum_required(VERSION 3.20)
+project(FESA LANGUAGES CXX)
+
+add_library(fesa_core INTERFACE)
+target_compile_features(fesa_core INTERFACE cxx_std_17)
+target_include_directories(fesa_core INTERFACE ${CMAKE_CURRENT_SOURCE_DIR}/src)
+
+enable_testing()
+add_subdirectory(tests)
+```
+
+`tests/unit/CMakeLists.txt`와 `tests/integration/CMakeLists.txt`는 `file(GLOB CONFIGURE_DEPENDS ...)`와 `foreach`를 사용해 새 `*_test.cpp`가 자동으로 CTest에 등록되도록 만든다.
+
+## Tests To Write First
+
+- `/tests/unit/harness_smoke_test.cpp`
+  - `main()`이 `std::string{"fesa"}.size() == 4` 같은 deterministic smoke assertion을 검증한다.
+  - 실패 시 non-zero를 반환한다.
+
+RED 확인:
+
+1. `tests/unit/harness_smoke_test.cpp`를 먼저 만든다.
+2. `ctest --test-dir build/msvc-debug --output-on-failure -C Debug -R harness_smoke_test`를 실행해 아직 CMake 구성이 없어 실패함을 확인한다.
+3. 그 뒤 CMake 파일을 작성한다.
+
+## Acceptance Criteria
+
+```powershell
+python -m unittest discover -s scripts -p "test_*.py"
+python scripts/validate_workspace.py
+ctest --test-dir build/msvc-debug --output-on-failure -C Debug -R harness_smoke_test
+```
+
+## 검증 절차
+
+1. 위 AC 커맨드를 실행한다.
+2. 아키텍처 체크리스트를 확인한다:
+   - ARCHITECTURE.md의 `src/`, `tests/unit/` 구조를 따르는가?
+   - ADR-002의 C++17/MSVC/CMake/CTest 기본값을 벗어나지 않았는가?
+   - AGENTS.md의 `python scripts/validate_workspace.py` 기본 검증 경로를 유지하는가?
+3. 결과에 따라 `phases/solver-core-skeleton/index.json`의 step 0을 업데이트한다:
+   - 성공: `"status": "completed"`, `"summary": "CMake/CTest bootstrap with fesa_core interface target and smoke test"`
+   - 3회 수정 시도 후 실패: `"status": "error"`, `"error_message": "구체적 에러 내용"`
+   - 사용자 개입 필요: `"status": "blocked"`, `"blocked_reason": "구체적 사유"` 후 중단
+
+## 금지사항
+
+- C++ production solver class를 이 step에서 만들지 마라.
+- 외부 test framework를 추가하지 마라.
+- JavaScript/TypeScript/npm fallback을 추가하지 마라.

phases/solver-core-skeleton/step1.md

@@ -0,0 +1,105 @@
+# Step 1: core-primitives
+
+## 읽어야 할 파일
+
+먼저 아래 파일들을 읽고 프로젝트의 아키텍처와 이전 step 산출물을 파악하라:
+
+- `/AGENTS.md`
+- `/docs/PRD.md`
+- `/docs/ARCHITECTURE.md`
+- `/docs/ADR.md`
+- `/CMakeLists.txt`
+- `/tests/unit/CMakeLists.txt`
+
+이전 step에서 만들어진 CMake/CTest 구성을 꼼꼼히 읽고, 새 unit test가 자동 등록되는 방식과 일관성을 유지하라.
+
+## 작업
+
+solver skeleton의 공통 primitive를 `/src/fesa/core/` 아래 header-only로 구현한다.
+
+필수 파일:
+
+- `/src/fesa/core/ids.hpp`
+- `/src/fesa/core/diagnostic.hpp`
+- `/src/fesa/core/status.hpp`
+- `/tests/unit/core_primitives_test.cpp`
+
+필수 interface:
+
+```cpp
+namespace fesa::core {
+
+struct NodeId { int value; };
+struct ElementId { int value; };
+struct MaterialId { int value; };
+struct PropertyId { int value; };
+struct StepId { int value; };
+
+enum class Severity { info, warning, error };
+
+struct Diagnostic {
+    Severity severity;
+    std::string code;
+    std::string message;
+};
+
+class Status {
+public:
+    static Status ok();
+    static Status failure(Diagnostic diagnostic);
+
+    bool is_ok() const;
+    const std::vector<Diagnostic>& diagnostics() const;
+    void add(Diagnostic diagnostic);
+};
+
+} // namespace fesa::core
+```
+
+구현 규칙:
+
+- 모든 ID type은 strong typedef 역할을 하며 서로 암시적으로 대체되지 않아야 한다.
+- ID에는 equation id 또는 DOF numbering 정보를 넣지 않는다.
+- `Status::ok()`는 diagnostics가 비어 있고 `is_ok() == true`여야 한다.
+- `Status::failure(...)`는 최소 하나의 diagnostic을 포함하고 `is_ok() == false`여야 한다.
+- 외부 라이브러리에 의존하지 않는다.
+
+## Tests To Write First
+
+- `/tests/unit/core_primitives_test.cpp`
+  - `NodeId{1}`과 `ElementId{1}`이 서로 다른 타입임을 compile-time으로 확인한다.
+  - `Status::ok().is_ok()`가 true임을 확인한다.
+  - `Status::failure(...)`가 false이고 diagnostic code/message를 보존함을 확인한다.
+  - `Status::add(...)`로 warning을 추가해도 diagnostic 순서가 보존됨을 확인한다.
+
+RED 확인:
+
+1. 테스트 파일을 먼저 작성한다.
+2. `python scripts/validate_workspace.py` 또는 targeted CTest를 실행해 `fesa/core/ids.hpp` 등 missing include로 실패함을 확인한다.
+3. 그 뒤 production header를 작성한다.
+
+## Acceptance Criteria
+
+```powershell
+python -m unittest discover -s scripts -p "test_*.py"
+python scripts/validate_workspace.py
+ctest --test-dir build/msvc-debug --output-on-failure -C Debug -R core_primitives_test
+```
+
+## 검증 절차
+
+1. 위 AC 커맨드를 실행한다.
+2. 아키텍처 체크리스트를 확인한다:
+   - `core`가 외부 라이브러리에 의존하지 않는가?
+   - ID type에 equation numbering을 저장하지 않는가?
+   - C++ production header 변경에 대응하는 test file이 있는가?
+3. 결과에 따라 `phases/solver-core-skeleton/index.json`의 step 1을 업데이트한다:
+   - 성공: `"status": "completed"`, `"summary": "Core ID, diagnostic, and status primitives added with tests"`
+   - 3회 수정 시도 후 실패: `"status": "error"`, `"error_message": "구체적 에러 내용"`
+   - 사용자 개입 필요: `"status": "blocked"`, `"blocked_reason": "구체적 사유"` 후 중단
+
+## 금지사항
+
+- Node, Element, Domain 같은 model class를 이 step에서 만들지 마라.
+- MKL, TBB, HDF5 API를 include하지 마라.
+- JavaScript/TypeScript/npm fallback을 추가하지 마라.

phases/solver-core-skeleton/step2.md

@@ -0,0 +1,175 @@
+# Step 2: domain-model-entities
+
+## 읽어야 할 파일
+
+먼저 아래 파일들을 읽고 프로젝트의 아키텍처와 이전 step 산출물을 파악하라:
+
+- `/AGENTS.md`
+- `/docs/PRD.md`
+- `/docs/ARCHITECTURE.md`
+- `/docs/ADR.md`
+- `/src/fesa/core/ids.hpp`
+- `/src/fesa/core/status.hpp`
+- `/tests/unit/core_primitives_test.cpp`
+
+이전 step에서 만들어진 core primitive를 꼼꼼히 읽고, ID ownership과 diagnostic convention을 유지하라.
+
+## 작업
+
+입력 파일에서 생성된 전체 semantic model을 소유하는 최소 model layer를 `/src/fesa/model/`에 구현한다.
+
+필수 파일:
+
+- `/src/fesa/model/node.hpp`
+- `/src/fesa/model/element.hpp`
+- `/src/fesa/model/material.hpp`
+- `/src/fesa/model/property.hpp`
+- `/src/fesa/model/boundary_condition.hpp`
+- `/src/fesa/model/load.hpp`
+- `/src/fesa/model/analysis_step.hpp`
+- `/src/fesa/model/domain.hpp`
+- `/tests/unit/model_domain_test.cpp`
+
+필수 interface:
+
+```cpp
+namespace fesa::model {
+
+class Node {
+public:
+    Node(core::NodeId id, std::array<double, 3> coordinates);
+    core::NodeId id() const;
+    const std::array<double, 3>& coordinates() const;
+};
+
+enum class ElementTopology { truss2, bar2, unknown };
+
+class Element {
+public:
+    Element(core::ElementId id, ElementTopology topology,
+            std::vector<core::NodeId> node_ids, core::PropertyId property_id);
+    core::ElementId id() const;
+    ElementTopology topology() const;
+    const std::vector<core::NodeId>& node_ids() const;
+    core::PropertyId property_id() const;
+};
+
+class Material {
+public:
+    Material(core::MaterialId id, std::string name);
+    core::MaterialId id() const;
+    const std::string& name() const;
+};
+
+class Property {
+public:
+    Property(core::PropertyId id, std::string name, core::MaterialId material_id);
+    core::PropertyId id() const;
+    const std::string& name() const;
+    core::MaterialId material_id() const;
+};
+
+enum class DofComponent { ux, uy, uz, rx, ry, rz, temperature };
+
+class BoundaryCondition {
+public:
+    BoundaryCondition(core::NodeId node_id, DofComponent component, double value);
+    core::NodeId node_id() const;
+    DofComponent component() const;
+    double value() const;
+};
+
+class Load {
+public:
+    Load(core::NodeId node_id, DofComponent component, double value);
+    core::NodeId node_id() const;
+    DofComponent component() const;
+    double value() const;
+};
+
+class AnalysisStep {
+public:
+    AnalysisStep(core::StepId id, std::string name);
+    core::StepId id() const;
+    const std::string& name() const;
+    void add_boundary_condition(BoundaryCondition bc);
+    void add_load(Load load);
+    const std::vector<BoundaryCondition>& boundary_conditions() const;
+    const std::vector<Load>& loads() const;
+};
+
+class Domain {
+public:
+    void add_node(Node node);
+    void add_element(Element element);
+    void add_material(Material material);
+    void add_property(Property property);
+    void add_step(AnalysisStep step);
+
+    const std::vector<Node>& nodes() const;
+    const std::vector<Element>& elements() const;
+    const std::vector<Material>& materials() const;
+    const std::vector<Property>& properties() const;
+    const std::vector<AnalysisStep>& steps() const;
+
+    const Node* find_node(core::NodeId id) const;
+    const Element* find_element(core::ElementId id) const;
+    const Material* find_material(core::MaterialId id) const;
+    const Property* find_property(core::PropertyId id) const;
+    const AnalysisStep* find_step(core::StepId id) const;
+};
+
+} // namespace fesa::model
+```
+
+구현 규칙:
+
+- Domain은 semantic model 객체를 소유한다.
+- `nodes()`, `elements()`, `materials()`, `properties()`, `steps()`는 const reference를 반환한다.
+- `find_*`는 없으면 `nullptr`을 반환한다.
+- Node와 Element 내부에 equation id, constrained/free equation mapping, sparse pattern 정보를 저장하지 않는다.
+- `DofComponent`는 아직 equation number가 아니라 물리 DOF component만 표현한다.
+- Abaqus keyword 문자열이나 parser detail을 model object에 저장하지 않는다.
+
+## Tests To Write First
+
+- `/tests/unit/model_domain_test.cpp`
+  - Node가 id와 3D coordinates를 보존한다.
+  - Element가 topology, connectivity, property id를 보존한다.
+  - Material과 Property가 id/name/material link를 보존한다.
+  - AnalysisStep이 boundary condition과 load를 저장한다.
+  - Domain이 add/find를 통해 각 객체를 조회한다.
+  - 없는 id는 `nullptr`을 반환한다.
+  - Node/Element public interface에 equation id를 노출하지 않는다는 점을 테스트 코드 사용 방식으로 확인한다.
+
+RED 확인:
+
+1. 테스트 파일을 먼저 작성한다.
+2. targeted CTest를 실행해 missing model headers로 실패함을 확인한다.
+3. 그 뒤 production headers를 작성한다.
+
+## Acceptance Criteria
+
+```powershell
+python -m unittest discover -s scripts -p "test_*.py"
+python scripts/validate_workspace.py
+ctest --test-dir build/msvc-debug --output-on-failure -C Debug -R model_domain_test
+```
+
+## 검증 절차
+
+1. 위 AC 커맨드를 실행한다.
+2. 아키텍처 체크리스트를 확인한다:
+   - Domain이 전체 모델 정의를 소유하는가?
+   - model layer가 parser keyword 문자열이나 analysis state를 소유하지 않는가?
+   - Node/Element에 equation id가 분산 저장되지 않는가?
+3. 결과에 따라 `phases/solver-core-skeleton/index.json`의 step 2를 업데이트한다:
+   - 성공: `"status": "completed"`, `"summary": "Model entities and Domain ownership API added with tests"`
+   - 3회 수정 시도 후 실패: `"status": "error"`, `"error_message": "구체적 에러 내용"`
+   - 사용자 개입 필요: `"status": "blocked"`, `"blocked_reason": "구체적 사유"` 후 중단
+
+## 금지사항
+
+- Parser, assembler, solver backend를 만들지 마라.
+- `Domain`을 실행 중 state container로 사용하지 마라.
+- JavaScript/TypeScript/npm fallback을 추가하지 마라.

phases/solver-core-skeleton/step3.md

@@ -0,0 +1,97 @@
+# Step 3: analysis-model-view
+
+## 읽어야 할 파일
+
+먼저 아래 파일들을 읽고 프로젝트의 아키텍처와 이전 step 산출물을 파악하라:
+
+- `/AGENTS.md`
+- `/docs/PRD.md`
+- `/docs/ARCHITECTURE.md`
+- `/docs/ADR.md`
+- `/src/fesa/model/domain.hpp`
+- `/src/fesa/model/analysis_step.hpp`
+- `/tests/unit/model_domain_test.cpp`
+
+이전 step에서 만들어진 Domain과 AnalysisStep을 꼼꼼히 읽고, Domain 복사 없는 step view를 유지하라.
+
+## 작업
+
+현재 step에서 활성화되는 해석 객체 view를 제공하는 `AnalysisModel`을 `/src/fesa/analysis/`에 구현한다.
+
+필수 파일:
+
+- `/src/fesa/analysis/analysis_model.hpp`
+- `/tests/unit/analysis_model_view_test.cpp`
+
+필수 interface:
+
+```cpp
+namespace fesa::analysis {
+
+class AnalysisModel {
+public:
+    AnalysisModel(const model::Domain& domain, core::StepId step_id);
+
+    const model::Domain& domain() const;
+    const model::AnalysisStep& step() const;
+    const std::vector<const model::Element*>& active_elements() const;
+    const std::vector<const model::BoundaryCondition*>& active_boundary_conditions() const;
+    const std::vector<const model::Load*>& active_loads() const;
+
+    const model::Property* property_for(const model::Element& element) const;
+    const model::Material* material_for(const model::Property& property) const;
+};
+
+} // namespace fesa::analysis
+```
+
+구현 규칙:
+
+- `AnalysisModel`은 `Domain`을 복사하지 않고 const reference 또는 pointer view만 보유한다.
+- Phase skeleton에서는 모든 Domain elements가 active라고 간주한다.
+- Boundary condition과 load는 선택된 `AnalysisStep`에서 가져온다.
+- `property_for`와 `material_for`는 Domain lookup을 사용한다.
+- 없는 step id는 구조화된 exception 대신 `std::invalid_argument`로 실패해도 된다. 이 skeleton 단계에서는 별도 error hierarchy를 만들지 않는다.
+- `AnalysisModel`은 displacement, residual, equation number를 소유하지 않는다.
+
+## Tests To Write First
+
+- `/tests/unit/analysis_model_view_test.cpp`
+  - Domain에 node/element/material/property/step을 만든다.
+  - `AnalysisModel(domain, step_id)`가 원본 Domain reference를 유지함을 pointer identity로 확인한다.
+  - 모든 Domain element가 `active_elements()`에 const pointer로 나타난다.
+  - step의 boundary condition과 load가 active view에 const pointer로 나타난다.
+  - `property_for(element)`와 `material_for(property)`가 Domain 소유 객체를 반환한다.
+  - 없는 step id는 `std::invalid_argument`를 던진다.
+
+RED 확인:
+
+1. 테스트 파일을 먼저 작성한다.
+2. targeted CTest를 실행해 missing `analysis_model.hpp`로 실패함을 확인한다.
+3. 그 뒤 production header를 작성한다.
+
+## Acceptance Criteria
+
+```powershell
+python -m unittest discover -s scripts -p "test_*.py"
+python scripts/validate_workspace.py
+ctest --test-dir build/msvc-debug --output-on-failure -C Debug -R analysis_model_view_test
+```
+
+## 검증 절차
+
+1. 위 AC 커맨드를 실행한다.
+2. 아키텍처 체크리스트를 확인한다:
+   - AnalysisModel이 Domain을 복사하지 않는가?
+   - AnalysisModel이 현재 step view만 제공하고 mutable state를 소유하지 않는가?
+   - active BC/load가 AnalysisStep에서 온 것임이 테스트되는가?
+3. 결과에 따라 `phases/solver-core-skeleton/index.json`의 step 3을 업데이트한다:
+   - 성공: `"status": "completed"`, `"summary": "AnalysisModel step view added without Domain copies"`
+   - 3회 수정 시도 후 실패: `"status": "error"`, `"error_message": "구체적 에러 내용"`
+   - 사용자 개입 필요: `"status": "blocked"`, `"blocked_reason": "구체적 사유"` 후 중단
+
+## 금지사항
+
+- DofManager나 equation numbering을 이 step에서 구현하지 마라.
+- AnalysisState를 이 step에서 구현하지 마라.
+- JavaScript/TypeScript/npm fallback을 추가하지 마라.

phases/solver-core-skeleton/step4.md

@@ -0,0 +1,112 @@
+# Step 4: dof-manager
+
+## 읽어야 할 파일
+
+먼저 아래 파일들을 읽고 프로젝트의 아키텍처와 이전 step 산출물을 파악하라:
+
+- `/AGENTS.md`
+- `/docs/PRD.md`
+- `/docs/ARCHITECTURE.md`
+- `/docs/ADR.md`
+- `/src/fesa/model/boundary_condition.hpp`
+- `/src/fesa/analysis/analysis_model.hpp`
+- `/tests/unit/analysis_model_view_test.cpp`
+
+이전 step에서 만들어진 model object와 AnalysisModel을 꼼꼼히 읽고, equation numbering이 Node/Element로 새지 않도록 유지하라.
+
+## 작업
+
+node별 DOF 정의, constrained/free mapping, equation numbering, sparse pattern ownership의 최소 골격을 `/src/fesa/fem/`에 구현한다.
+
+필수 파일:
+
+- `/src/fesa/fem/dof_key.hpp`
+- `/src/fesa/fem/dof_manager.hpp`
+- `/tests/unit/dof_manager_numbering_test.cpp`
+
+필수 interface:
+
+```cpp
+namespace fesa::fem {
+
+struct DofKey {
+    core::NodeId node_id;
+    model::DofComponent component;
+};
+
+bool operator==(const DofKey& lhs, const DofKey& rhs);
+
+class DofManager {
+public:
+    void define_node_dofs(core::NodeId node_id, std::vector<model::DofComponent> components);
+    void apply_boundary_condition(const model::BoundaryCondition& bc);
+    void number_equations();
+
+    int total_dof_count() const;
+    int free_dof_count() const;
+    int constrained_dof_count() const;
+
+    bool is_constrained(DofKey key) const;
+    int equation_id(DofKey key) const;
+    std::optional<int> free_equation_id(DofKey key) const;
+
+    std::vector<double> expand_free_vector(const std::vector<double>& free_values) const;
+    const std::vector<std::pair<int, int>>& sparse_pattern() const;
+};
+
+} // namespace fesa::fem
+```
+
+구현 규칙:
+
+- Equation id는 `DofManager` 내부에만 저장한다.
+- `equation_id`는 전체 DOF ordering의 dense id를 반환한다.
+- `free_equation_id`는 constrained DOF면 `std::nullopt`를 반환한다.
+- `number_equations()`는 deterministic ordering을 보장한다:
+  - node id value 오름차순
+  - component enum 순서 오름차순
+- `expand_free_vector`는 constrained DOF 값을 0.0으로 채우고 free DOF 값만 배치한다.
+- `sparse_pattern()`은 skeleton 단계에서 free equation ids의 dense full matrix pattern을 deterministic pair list로 보유해도 된다.
+- missing DOF 조회는 `std::invalid_argument`를 던진다.
+
+## Tests To Write First
+
+- `/tests/unit/dof_manager_numbering_test.cpp`
+  - 두 node에 `ux`, `uy`를 정의하고 deterministic equation id ordering을 확인한다.
+  - boundary condition 적용 후 constrained/free count가 맞는지 확인한다.
+  - constrained key의 `free_equation_id`가 `std::nullopt`임을 확인한다.
+  - free vector가 full vector로 재구성될 때 constrained DOF가 0.0으로 채워지는지 확인한다.
+  - sparse pattern pair list가 free equation id 기반 deterministic dense pattern을 가진다.
+  - model::Node나 model::Element를 수정하지 않고도 equation numbering이 가능함을 확인한다.
+
+RED 확인:
+
+1. 테스트 파일을 먼저 작성한다.
+2. targeted CTest를 실행해 missing `dof_manager.hpp`로 실패함을 확인한다.
+3. 그 뒤 production headers를 작성한다.
+
+## Acceptance Criteria
+
+```powershell
+python -m unittest discover -s scripts -p "test_*.py"
+python scripts/validate_workspace.py
+ctest --test-dir build/msvc-debug --output-on-failure -C Debug -R dof_manager_numbering_test
+```
+
+## 검증 절차
+
+1. 위 AC 커맨드를 실행한다.
+2. 아키텍처 체크리스트를 확인한다:
+   - DofManager가 equation numbering을 전담하는가?
+   - Node/Element public interface가 equation id로 오염되지 않았는가?
+   - sparse pattern ownership이 DofManager 내부에 있는가?
+3. 결과에 따라 `phases/solver-core-skeleton/index.json`의 step 4를 업데이트한다:
+   - 성공: `"status": "completed"`, `"summary": "DofManager deterministic numbering and constrained/free mapping added"`
+   - 3회 수정 시도 후 실패: `"status": "error"`, `"error_message": "구체적 에러 내용"`
+   - 사용자 개입 필요: `"status": "blocked"`, `"blocked_reason": "구체적 사유"` 후 중단
+
+## 금지사항
+
+- Solver backend, assembly matrix, MKL adapter를 구현하지 마라.
+- Node 또는 Element에 equation id field를 추가하지 마라.
+- JavaScript/TypeScript/npm fallback을 추가하지 마라.

phases/solver-core-skeleton/step5.md

@@ -0,0 +1,112 @@
+# Step 5: analysis-state
+
+## 읽어야 할 파일
+
+먼저 아래 파일들을 읽고 프로젝트의 아키텍처와 이전 step 산출물을 파악하라:
+
+- `/AGENTS.md`
+- `/docs/PRD.md`
+- `/docs/ARCHITECTURE.md`
+- `/docs/ADR.md`
+- `/src/fesa/fem/dof_manager.hpp`
+- `/tests/unit/dof_manager_numbering_test.cpp`
+
+이전 step에서 만들어진 DofManager를 꼼꼼히 읽고, 해석 중 변하는 물리량은 AnalysisState가 소유하도록 유지하라.
+
+## 작업
+
+displacement 중심의 최소 `AnalysisState`를 `/src/fesa/analysis/`에 구현한다.
+
+필수 파일:
+
+- `/src/fesa/analysis/analysis_state.hpp`
+- `/tests/unit/analysis_state_vectors_test.cpp`
+
+필수 interface:
+
+```cpp
+namespace fesa::analysis {
+
+struct IterationState {
+    double time = 0.0;
+    int increment = 0;
+    int iteration = 0;
+};
+
+class AnalysisState {
+public:
+    explicit AnalysisState(int total_dof_count);
+
+    const std::vector<double>& displacement() const;
+    const std::vector<double>& velocity() const;
+    const std::vector<double>& acceleration() const;
+    const std::vector<double>& temperature() const;
+    const std::vector<double>& external_force() const;
+    const std::vector<double>& internal_force() const;
+    const std::vector<double>& residual() const;
+
+    void set_displacement(std::vector<double> values);
+    void set_external_force(std::vector<double> values);
+    void set_internal_force(std::vector<double> values);
+    void update_residual();
+
+    IterationState& iteration_state();
+    const IterationState& iteration_state() const;
+
+    void set_element_state(core::ElementId element_id, std::vector<double> state);
+    const std::vector<double>* element_state(core::ElementId element_id) const;
+};
+
+} // namespace fesa::analysis
+```
+
+구현 규칙:
+
+- 모든 vector는 `total_dof_count` 크기로 초기화한다.
+- `update_residual()`은 `external_force - internal_force`를 component-wise로 계산한다.
+- setter는 입력 vector 크기가 맞지 않으면 `std::invalid_argument`를 던진다.
+- temperature는 Phase skeleton에서 0.0 vector로 둔다.
+- element state는 향후 integration point state 확장을 위한 최소 map으로 둔다.
+- AnalysisState는 Domain, AnalysisModel, DofManager를 소유하지 않는다.
+
+## Tests To Write First
+
+- `/tests/unit/analysis_state_vectors_test.cpp`
+  - 생성 시 displacement/velocity/acceleration/temperature/force/residual vector 크기와 0.0 초기값을 확인한다.
+  - displacement setter가 값을 보존한다.
+  - external/internal force 설정 후 residual이 `Fext - Fint`가 되는지 확인한다.
+  - 크기가 맞지 않는 vector setter가 `std::invalid_argument`를 던진다.
+  - time/increment/iteration 값이 저장된다.
+  - element state를 element id로 저장/조회한다.
+
+RED 확인:
+
+1. 테스트 파일을 먼저 작성한다.
+2. targeted CTest를 실행해 missing `analysis_state.hpp`로 실패함을 확인한다.
+3. 그 뒤 production header를 작성한다.
+
+## Acceptance Criteria
+
+```powershell
+python -m unittest discover -s scripts -p "test_*.py"
+python scripts/validate_workspace.py
+ctest --test-dir build/msvc-debug --output-on-failure -C Debug -R analysis_state_vectors_test
+```
+
+## 검증 절차
+
+1. 위 AC 커맨드를 실행한다.
+2. 아키텍처 체크리스트를 확인한다:
+   - AnalysisState가 해석 중 변하는 물리량을 소유하는가?
+   - Domain이나 AnalysisModel을 복사/소유하지 않는가?
+   - displacement 중심이되 velocity/acceleration/temperature 확장 지점을 유지하는가?
+3. 결과에 따라 `phases/solver-core-skeleton/index.json`의 step 5를 업데이트한다:
+   - 성공: `"status": "completed"`, `"summary": "AnalysisState vector ownership and residual update added"`
+   - 3회 수정 시도 후 실패: `"status": "error"`, `"error_message": "구체적 에러 내용"`
+   - 사용자 개입 필요: `"status": "blocked"`, `"blocked_reason": "구체적 사유"` 후 중단
+
+## 금지사항
+
+- Solver backend나 numerical integration loop를 구현하지 마라.
+- HDF5 writer를 이 step에서 구현하지 마라.
+- JavaScript/TypeScript/npm fallback을 추가하지 마라.

phases/solver-core-skeleton/step6.md

@@ -0,0 +1,114 @@
+# Step 6: analysis-template-flow
+
+## 읽어야 할 파일
+
+먼저 아래 파일들을 읽고 프로젝트의 아키텍처와 이전 step 산출물을 파악하라:
+
+- `/AGENTS.md`
+- `/docs/PRD.md`
+- `/docs/ARCHITECTURE.md`
+- `/docs/ADR.md`
+- `/src/fesa/analysis/analysis_model.hpp`
+- `/src/fesa/analysis/analysis_state.hpp`
+- `/src/fesa/fem/dof_manager.hpp`
+- `/tests/unit/analysis_state_vectors_test.cpp`
+
+이전 step에서 만들어진 AnalysisModel, DofManager, AnalysisState를 꼼꼼히 읽고, ARCHITECTURE.md의 Template Method 실행 흐름과 일관성을 유지하라.
+
+## 작업
+
+공통 해석 실행 흐름을 고정하는 `Analysis` base class와 선형 정적 해석 skeleton을 `/src/fesa/analysis/`에 구현한다.
+
+필수 파일:
+
+- `/src/fesa/analysis/analysis.hpp`
+- `/src/fesa/analysis/linear_static_analysis.hpp`
+- `/tests/unit/analysis_flow_template_test.cpp`
+
+필수 interface:
+
+```cpp
+namespace fesa::analysis {
+
+class Analysis {
+public:
+    virtual ~Analysis() = default;
+    void run();
+
+protected:
+    virtual void initialize() {}
+    virtual void build_analysis_model() {}
+    virtual void build_dof_map() {}
+    virtual void build_sparse_pattern() {}
+    virtual void assemble() {}
+    virtual void apply_boundary_conditions() {}
+    virtual void solve() {}
+    virtual void update_state() {}
+    virtual void write_results() {}
+};
+
+class LinearStaticAnalysis : public Analysis {
+public:
+    LinearStaticAnalysis(const model::Domain& domain, core::StepId step_id);
+
+    const AnalysisModel* analysis_model() const;
+    const AnalysisState* state() const;
+
+protected:
+    void build_analysis_model() override;
+    void build_dof_map() override;
+    void update_state() override;
+};
+
+} // namespace fesa::analysis
+```
+
+구현 규칙:
+
+- `Analysis::run()`은 반드시 다음 순서로 hook을 호출한다:
+  `initialize -> build_analysis_model -> build_dof_map -> build_sparse_pattern -> assemble -> apply_boundary_conditions -> solve -> update_state -> write_results`
+- `LinearStaticAnalysis`는 skeleton 단계에서 실제 stiffness assembly나 solve를 하지 않는다.
+- `LinearStaticAnalysis::build_analysis_model()`은 `AnalysisModel`을 생성한다.
+- `LinearStaticAnalysis::build_dof_map()`은 active element connectivity에서 등장한 node에 `ux`, `uy`, `uz`를 정의하고 active BC를 적용한 뒤 equation numbering을 수행한다.
+- `LinearStaticAnalysis::update_state()`는 total DOF count 크기의 `AnalysisState`를 준비한다.
+- MKL, TBB, HDF5 adapter는 만들지 않는다.
+
+## Tests To Write First
+
+- `/tests/unit/analysis_flow_template_test.cpp`
+  - test-only derived `RecordingAnalysis`가 hook 호출 순서를 vector에 기록한다.
+  - `run()` 호출 후 ARCHITECTURE.md 순서와 정확히 일치하는지 확인한다.
+  - 최소 Domain과 Step으로 `LinearStaticAnalysis`를 실행하면 `analysis_model()`과 `state()`가 null이 아니게 되는지 확인한다.
+  - `LinearStaticAnalysis`가 실제 solver 결과를 계산한다고 주장하지 않음을 테스트 이름과 assertion 범위에 반영한다.
+
+RED 확인:
+
+1. 테스트 파일을 먼저 작성한다.
+2. targeted CTest를 실행해 missing `analysis.hpp` 또는 `linear_static_analysis.hpp`로 실패함을 확인한다.
+3. 그 뒤 production headers를 작성한다.
+
+## Acceptance Criteria
+
+```powershell
+python -m unittest discover -s scripts -p "test_*.py"
+python scripts/validate_workspace.py
+ctest --test-dir build/msvc-debug --output-on-failure -C Debug -R analysis_flow_template_test
+```
+
+## 검증 절차
+
+1. 위 AC 커맨드를 실행한다.
+2. 아키텍처 체크리스트를 확인한다:
+   - Analysis::run()이 Template Method 흐름을 고정하는가?
+   - LinearStaticAnalysis가 skeleton 범위를 넘어 solver backend를 구현하지 않는가?
+   - 외부 라이브러리 API가 solver core에 노출되지 않는가?
+3. 결과에 따라 `phases/solver-core-skeleton/index.json`의 step 6을 업데이트한다:
+   - 성공: `"status": "completed"`, `"summary": "Analysis template method and LinearStaticAnalysis skeleton added"`
+   - 3회 수정 시도 후 실패: `"status": "error"`, `"error_message": "구체적 에러 내용"`
+   - 사용자 개입 필요: `"status": "blocked"`, `"blocked_reason": "구체적 사유"` 후 중단
+
+## 금지사항
+
+- 실제 stiffness matrix assembly, linear solve, MKL PARDISO adapter를 구현하지 마라.
+- HDF5 writer를 구현하지 마라.
+- JavaScript/TypeScript/npm fallback을 추가하지 마라.

phases/solver-core-skeleton/step7.md

@@ -0,0 +1,118 @@
+# Step 7: results-containers
+
+## 읽어야 할 파일
+
+먼저 아래 파일들을 읽고 프로젝트의 아키텍처와 이전 step 산출물을 파악하라:
+
+- `/AGENTS.md`
+- `/docs/PRD.md`
+- `/docs/ARCHITECTURE.md`
+- `/docs/ADR.md`
+- `/src/fesa/analysis/analysis_state.hpp`
+- `/src/fesa/analysis/analysis.hpp`
+- `/tests/unit/analysis_flow_template_test.cpp`
+
+이전 step에서 만들어진 AnalysisState와 Analysis flow를 꼼꼼히 읽고, 결과 container가 HDF5 API에 직접 의존하지 않도록 유지하라.
+
+## 작업
+
+HDF5 writer 구현 전 단계의 results data model skeleton을 `/src/fesa/results/`에 구현한다.
+
+필수 파일:
+
+- `/src/fesa/results/results.hpp`
+- `/tests/unit/results_containers_test.cpp`
+
+필수 interface:
+
+```cpp
+namespace fesa::results {
+
+enum class FieldLocation { nodal, element, integration_point };
+
+struct FieldOutput {
+    std::string name;
+    FieldLocation location;
+    std::vector<std::string> components;
+    std::vector<int> entity_ids;
+    std::vector<double> values;
+};
+
+struct HistoryOutput {
+    std::string name;
+    std::vector<double> time;
+    std::vector<double> values;
+};
+
+class ResultFrame {
+public:
+    ResultFrame(int frame_id, double time);
+    int frame_id() const;
+    double time() const;
+    void add_field_output(FieldOutput output);
+    void add_history_output(HistoryOutput output);
+    const std::vector<FieldOutput>& field_outputs() const;
+    const std::vector<HistoryOutput>& history_outputs() const;
+};
+
+class ResultStep {
+public:
+    explicit ResultStep(std::string name);
+    const std::string& name() const;
+    ResultFrame& add_frame(int frame_id, double time);
+    const std::vector<ResultFrame>& frames() const;
+};
+
+} // namespace fesa::results
+```
+
+구현 규칙:
+
+- Result hierarchy는 `ResultStep -> ResultFrame -> FieldOutput/HistoryOutput` 구조를 따른다.
+- Field output은 row identity를 위해 `entity_ids`를 보존한다.
+- Field output values layout은 skeleton 단계에서 row-major flat vector로 둔다.
+- HDF5 file/dataset, schema writer, CSV export는 구현하지 않는다.
+- validation은 크기 consistency만 최소로 확인한다:
+  - `components`가 비어 있으면 `std::invalid_argument`
+  - `entity_ids.size() * components.size() == values.size()`가 아니면 `std::invalid_argument`
+
+## Tests To Write First
+
+- `/tests/unit/results_containers_test.cpp`
+  - ResultStep이 이름과 frame 목록을 보존한다.
+  - ResultFrame이 frame id/time을 보존한다.
+  - nodal displacement FieldOutput을 추가하면 components/entity ids/values가 보존된다.
+  - invalid FieldOutput shape가 `std::invalid_argument`를 던진다.
+  - HistoryOutput이 time/value series를 보존한다.
+
+RED 확인:
+
+1. 테스트 파일을 먼저 작성한다.
+2. targeted CTest를 실행해 missing `results.hpp`로 실패함을 확인한다.
+3. 그 뒤 production header를 작성한다.
+
+## Acceptance Criteria
+
+```powershell
+python -m unittest discover -s scripts -p "test_*.py"
+python scripts/validate_workspace.py
+ctest --test-dir build/msvc-debug --output-on-failure -C Debug -R results_containers_test
+```
+
+## 검증 절차
+
+1. 위 AC 커맨드를 실행한다.
+2. 아키텍처 체크리스트를 확인한다:
+   - Results hierarchy가 ARCHITECTURE.md와 일치하는가?
+   - HDF5 API가 results container에 직접 노출되지 않는가?
+   - reference comparison을 위한 entity row identity가 보존되는가?
+3. 결과에 따라 `phases/solver-core-skeleton/index.json`의 step 7을 업데이트한다:
+   - 성공: `"status": "completed"`, `"summary": "ResultStep, ResultFrame, FieldOutput, and HistoryOutput containers added"`
+   - 3회 수정 시도 후 실패: `"status": "error"`, `"error_message": "구체적 에러 내용"`
+   - 사용자 개입 필요: `"status": "blocked"`, `"blocked_reason": "구체적 사유"` 후 중단
+
+## 금지사항
+
+- HDF5 writer/reader를 구현하지 마라.
+- deterministic CSV export를 구현하지 마라.
+- JavaScript/TypeScript/npm fallback을 추가하지 마라.

phases/solver-core-skeleton/step8.md

@@ -0,0 +1,89 @@
+# Step 8: solver-skeleton-integration-report
+
+## 읽어야 할 파일
+
+먼저 아래 파일들을 읽고 프로젝트의 아키텍처와 이전 step 산출물을 파악하라:
+
+- `/AGENTS.md`
+- `/docs/PRD.md`
+- `/docs/ARCHITECTURE.md`
+- `/docs/ADR.md`
+- `/src/fesa/model/domain.hpp`
+- `/src/fesa/analysis/linear_static_analysis.hpp`
+- `/src/fesa/results/results.hpp`
+- `/tests/unit/results_containers_test.cpp`
+- `/docs/build-test-reports/README.md`
+
+이전 step에서 만들어진 전체 skeleton API를 꼼꼼히 읽고, integration smoke test가 실제 solver 수치해석 완료를 주장하지 않도록 범위를 제한하라.
+
+## 작업
+
+solver skeleton의 주요 class가 함께 컴파일되고 기본 data flow를 구성할 수 있음을 통합 테스트와 build/test report로 남긴다.
+
+필수 파일:
+
+- `/tests/integration/solver_core_skeleton_integration_test.cpp`
+- `/docs/build-test-reports/solver-core-skeleton.md`
+
+통합 테스트 요구사항:
+
+- Domain에 두 개 Node, 하나 Element, Material, Property, AnalysisStep을 구성한다.
+- AnalysisStep에는 최소 하나의 BoundaryCondition과 Load를 추가한다.
+- `LinearStaticAnalysis`를 생성하고 `run()`을 호출한다.
+- `analysis_model()`과 `state()`가 생성되는지 확인한다.
+- `ResultStep`과 `ResultFrame`을 만들고 nodal displacement `FieldOutput`을 추가한다.
+- 이 테스트는 stiffness assembly, linear solve, HDF5 write, reference comparison을 검증하지 않는다.
+
+보고서 요구사항:
+
+- `docs/build-test-reports/solver-core-skeleton.md`에 아래 항목을 기록한다.
+  - phase: solver-core-skeleton
+  - scope: C++ skeleton classes only
+  - commands run
+  - exit code summary
+  - CTest tests added
+  - known limitations
+- known limitations에는 최소한 다음을 명시한다:
+  - 실제 element stiffness/residual/tangent 계산 없음
+  - 실제 linear solver backend 없음
+  - HDF5 writer 없음
+  - Abaqus parser 및 reference comparison 없음
+
+## Tests To Write First
+
+- `/tests/integration/solver_core_skeleton_integration_test.cpp`
+  - 전체 skeleton header를 include한다.
+  - 위 통합 테스트 요구사항을 `main()` assertion으로 검증한다.
+
+RED 확인:
+
+1. 통합 테스트 파일을 먼저 작성한다.
+2. targeted CTest를 실행해 아직 integration test registration 또는 required API 문제로 실패하면 실패 내용을 확인한다.
+3. Step 0의 CMake glob이 `tests/integration/*_test.cpp`를 자동 등록해야 한다. 등록되지 않는다면 현재 step에서 CMake 파일을 수정하지 말고 `blocked`로 표시하고 사용자에게 allowed_paths 확장을 요청한다.
+
+## Acceptance Criteria
+
+```powershell
+python -m unittest discover -s scripts -p "test_*.py"
+python scripts/validate_workspace.py
+ctest --test-dir build/msvc-debug --output-on-failure -C Debug -R solver_core_skeleton_integration_test
+```
+
+## 검증 절차
+
+1. 위 AC 커맨드를 실행한다.
+2. 아키텍처 체크리스트를 확인한다:
+   - Domain -> AnalysisModel -> DofManager/AnalysisState -> Results data flow가 컴파일되는가?
+   - 통합 테스트가 skeleton 범위를 넘어 수치 정확성을 주장하지 않는가?
+   - build/test report가 실제 실행 evidence와 known limitations를 기록하는가?
+3. 결과에 따라 `phases/solver-core-skeleton/index.json`의 step 8을 업데이트한다:
+   - 성공: `"status": "completed"`, `"summary": "Solver skeleton integration test and build/test report added"`
+   - 3회 수정 시도 후 실패: `"status": "error"`, `"error_message": "구체적 에러 내용"`
+   - 사용자 개입 필요: `"status": "blocked"`, `"blocked_reason": "구체적 사유"` 후 중단
+
+## 금지사항
+
+- 실제 FEM stiffness, residual, tangent, material law 계산을 구현하지 마라.
+- Abaqus reference artifact를 생성, 수정, 복원하지 마라.
+- HDF5 writer를 구현하지 마라.
+- JavaScript/TypeScript/npm fallback을 추가하지 마라.

scripts/execute.py

@@ -12,6 +12,7 @@ import fnmatch
 import json
 import os
 import re
+import shutil
 import subprocess
 import sys
 import threading
@@ -24,6 +25,12 @@ from typing import Optional
 ROOT = Path(__file__).resolve().parent.parent
 
 
+def configure_output_encoding():
+    for stream in (sys.stdout, sys.stderr):
+        if hasattr(stream, "reconfigure"):
+            stream.reconfigure(encoding="utf-8", errors="replace")
+
+
 @contextlib.contextmanager
 def progress_indicator(label: str):
     """터미널 진행 표시기. with 문으로 사용하며 .elapsed 로 경과 시간을 읽는다."""
@@ -373,8 +380,9 @@ class StepExecutor:
 
         prompt = preamble + step_file.read_text(encoding="utf-8")
         result = subprocess.run(
-            ["codex", "exec", "--dangerously-bypass-approvals-and-sandbox", "--json", prompt],
-            cwd=self._root, capture_output=True, text=True, timeout=1800,
+            self._codex_exec_command(),
+            cwd=self._root, capture_output=True, text=True, input=prompt,
+            encoding="utf-8", errors="replace", timeout=1800,
         )
 
         if result.returncode != 0:
@@ -395,6 +403,26 @@ class StepExecutor:
 
     # --- 헤더 & 검증 ---
 
+    @staticmethod
+    def _codex_command() -> str:
+        override = os.environ.get("HARNESS_CODEX_COMMAND", "").strip()
+        if override:
+            return override
+        return (
+            shutil.which("codex.cmd")
+            or shutil.which("codex.exe")
+            or shutil.which("codex")
+            or "codex"
+        )
+
+    def _codex_exec_command(self) -> list[str]:
+        cmd = [self._codex_command(), "exec"]
+        model = os.environ.get("HARNESS_CODEX_MODEL", "").strip()
+        if model:
+            cmd.extend(["-m", model])
+        cmd.extend(["--dangerously-bypass-approvals-and-sandbox", "--json", "-"])
+        return cmd
+
     def _print_header(self):
         print(f"\n{'='*60}")
         print(f"  Harness Step Executor")
@@ -550,6 +578,7 @@ class StepExecutor:
 
 
 def main():
+    configure_output_encoding()
     parser = argparse.ArgumentParser(description="Harness Step Executor")
     parser.add_argument("phase_dir", help="Phase directory name (e.g. 0-mvp)")
     parser.add_argument("--push", action="store_true", help="Push branch after completion")

scripts/test_execute.py

@@ -1,5 +1,6 @@
 import importlib.util
 import json
+import os
 import subprocess
 import sys
 import tempfile
@@ -281,6 +282,94 @@ class ExecuteRunnerSafetyTests(unittest.TestCase):
 
             self.assertEqual(cm.exception.code, 1)
 
+    def test_invoke_codex_passes_prompt_through_stdin(self):
+        execute = load_execute()
+        with tempfile.TemporaryDirectory() as tmp:
+            root = Path(tmp)
+            write_phase(root)
+            executor = make_executor(execute, root)
+            step = {"step": 1, "name": "Docs"}
+            long_preamble = "x" * 40000
+
+            def fake_run(cmd, **kwargs):
+                return subprocess.CompletedProcess(cmd, 0, '{"event":"done"}\n', "")
+
+            with patch.object(execute.subprocess, "run", side_effect=fake_run) as run_mock:
+                with patch.object(executor, "_codex_command", return_value="codex.cmd"):
+                    with patch.dict(os.environ, {"HARNESS_CODEX_MODEL": ""}):
+                        executor._invoke_codex(step, long_preamble)
+
+            cmd = run_mock.call_args.args[0]
+            kwargs = run_mock.call_args.kwargs
+            self.assertEqual(
+                cmd,
+                ["codex.cmd", "exec", "--dangerously-bypass-approvals-and-sandbox", "--json", "-"],
+            )
+            self.assertEqual(kwargs["input"], long_preamble + "# Step 1\n")
+            self.assertEqual(kwargs["cwd"], str(root))
+            self.assertEqual(kwargs["encoding"], "utf-8")
+            self.assertEqual(kwargs["errors"], "replace")
+
+    def test_codex_command_prefers_windows_cmd_or_exe_shim(self):
+        execute = load_execute()
+
+        def fake_which(name):
+            return {
+                "codex.cmd": "C:/tools/codex.cmd",
+                "codex.exe": "C:/tools/codex.exe",
+                "codex": "C:/tools/codex",
+            }.get(name)
+
+        with patch.object(execute.shutil, "which", side_effect=fake_which):
+            self.assertEqual(execute.StepExecutor._codex_command(), "C:/tools/codex.cmd")
+
+    def test_codex_command_uses_env_override_first(self):
+        execute = load_execute()
+
+        with patch.dict(os.environ, {"HARNESS_CODEX_COMMAND": "C:/app/codex.exe"}):
+            self.assertEqual(execute.StepExecutor._codex_command(), "C:/app/codex.exe")
+
+    def test_codex_exec_command_uses_model_env_override(self):
+        execute = load_execute()
+        with tempfile.TemporaryDirectory() as tmp:
+            root = Path(tmp)
+            write_phase(root)
+            executor = make_executor(execute, root)
+
+            with patch.object(executor, "_codex_command", return_value="codex.cmd"):
+                with patch.dict(os.environ, {"HARNESS_CODEX_MODEL": "gpt-5"}):
+                    self.assertEqual(
+                        executor._codex_exec_command(),
+                        [
+                            "codex.cmd",
+                            "exec",
+                            "-m",
+                            "gpt-5",
+                            "--dangerously-bypass-approvals-and-sandbox",
+                            "--json",
+                            "-",
+                        ],
+                    )
+
+    def test_configure_output_encoding_sets_utf8_replace(self):
+        execute = load_execute()
+
+        class Stream:
+            def __init__(self):
+                self.calls = []
+
+            def reconfigure(self, **kwargs):
+                self.calls.append(kwargs)
+
+        stdout = Stream()
+        stderr = Stream()
+        with patch.object(execute.sys, "stdout", stdout):
+            with patch.object(execute.sys, "stderr", stderr):
+                execute.configure_output_encoding()
+
+        self.assertEqual(stdout.calls, [{"encoding": "utf-8", "errors": "replace"}])
+        self.assertEqual(stderr.calls, [{"encoding": "utf-8", "errors": "replace"}])
+
 
 if __name__ == "__main__":
     unittest.main()

scripts/test_header_declaration_only.py

@@ -0,0 +1,45 @@
+import re
+import unittest
+from pathlib import Path
+
+
+ROOT = Path(__file__).resolve().parents[1]
+HEADER_ROOT = ROOT / "src" / "fesa"
+
+
+def _looks_like_function_definition(prefix):
+    if "(" not in prefix or ")" not in prefix:
+        return False
+
+    stripped = prefix.strip()
+    control_prefixes = ("if ", "for ", "while ", "switch ", "catch ")
+    declaration_prefixes = ("namespace ", "class ", "struct ", "enum ")
+    return not stripped.startswith(control_prefixes + declaration_prefixes)
+
+
+class HeaderDeclarationOnlyTests(unittest.TestCase):
+    def test_solver_headers_do_not_contain_function_bodies(self):
+        violations = []
+
+        for header in sorted(HEADER_ROOT.rglob("*.hpp")):
+            candidate = ""
+            for line_number, line in enumerate(header.read_text(encoding="utf-8").splitlines(), start=1):
+                stripped = line.strip()
+                if not stripped:
+                    continue
+
+                candidate = f"{candidate} {stripped}".strip()
+                if "{" in stripped and _looks_like_function_definition(candidate.split("{", 1)[0]):
+                    violations.append(f"{header.relative_to(ROOT)}:{line_number}: function body in header")
+
+                if re.search(r"\([^;{}]*\)\s*=\s*(default|delete)\s*;", stripped):
+                    violations.append(f"{header.relative_to(ROOT)}:{line_number}: function definition in header")
+
+                if stripped.endswith(";") or stripped.endswith("}") or stripped.endswith(":"):
+                    candidate = ""
+
+        self.assertEqual([], violations)
+
+
+if __name__ == "__main__":
+    unittest.main()

src/fesa/analysis/analysis.cpp

@@ -0,0 +1,30 @@
+#include <fesa/analysis/analysis.hpp>
+
+namespace fesa::analysis {
+
+Analysis::~Analysis() = default;
+
+void Analysis::run()
+{
+    initialize();
+    build_analysis_model();
+    build_dof_map();
+    build_sparse_pattern();
+    assemble();
+    apply_boundary_conditions();
+    solve();
+    update_state();
+    write_results();
+}
+
+void Analysis::initialize() {}
+void Analysis::build_analysis_model() {}
+void Analysis::build_dof_map() {}
+void Analysis::build_sparse_pattern() {}
+void Analysis::assemble() {}
+void Analysis::apply_boundary_conditions() {}
+void Analysis::solve() {}
+void Analysis::update_state() {}
+void Analysis::write_results() {}
+
+} // namespace fesa::analysis

src/fesa/analysis/analysis.hpp

@@ -0,0 +1,23 @@
+#pragma once
+
+namespace fesa::analysis {
+
+class Analysis {
+public:
+    virtual ~Analysis();
+
+    void run();
+
+protected:
+    virtual void initialize();
+    virtual void build_analysis_model();
+    virtual void build_dof_map();
+    virtual void build_sparse_pattern();
+    virtual void assemble();
+    virtual void apply_boundary_conditions();
+    virtual void solve();
+    virtual void update_state();
+    virtual void write_results();
+};
+
+} // namespace fesa::analysis

src/fesa/analysis/analysis_model.cpp

@@ -0,0 +1,60 @@
+#include <fesa/analysis/analysis_model.hpp>
+
+#include <stdexcept>
+
+namespace fesa::analysis {
+
+AnalysisModel::AnalysisModel(const model::Domain& domain, core::StepId step_id)
+    : domain_(domain), step_(domain.find_step(step_id))
+{
+    if (step_ == nullptr) {
+        throw std::invalid_argument("analysis step not found");
+    }
+
+    for (const auto& element : domain_.elements()) {
+        active_elements_.push_back(&element);
+    }
+    for (const auto& boundary_condition : step_->boundary_conditions()) {
+        active_boundary_conditions_.push_back(&boundary_condition);
+    }
+    for (const auto& load : step_->loads()) {
+        active_loads_.push_back(&load);
+    }
+}
+
+const model::Domain& AnalysisModel::domain() const
+{
+    return domain_;
+}
+
+const model::AnalysisStep& AnalysisModel::step() const
+{
+    return *step_;
+}
+
+const std::vector<const model::Element*>& AnalysisModel::active_elements() const
+{
+    return active_elements_;
+}
+
+const std::vector<const model::BoundaryCondition*>& AnalysisModel::active_boundary_conditions() const
+{
+    return active_boundary_conditions_;
+}
+
+const std::vector<const model::Load*>& AnalysisModel::active_loads() const
+{
+    return active_loads_;
+}
+
+const model::Property* AnalysisModel::property_for(const model::Element& element) const
+{
+    return domain_.find_property(element.property_id());
+}
+
+const model::Material* AnalysisModel::material_for(const model::Property& property) const
+{
+    return domain_.find_material(property.material_id());
+}
+
+} // namespace fesa::analysis

src/fesa/analysis/analysis_model.hpp

@@ -0,0 +1,29 @@
+#pragma once
+
+#include <fesa/model/domain.hpp>
+
+#include <vector>
+
+namespace fesa::analysis {
+
+class AnalysisModel {
+public:
+    AnalysisModel(const model::Domain& domain, core::StepId step_id);
+
+    const model::Domain& domain() const;
+    const model::AnalysisStep& step() const;
+    const std::vector<const model::Element*>& active_elements() const;
+    const std::vector<const model::BoundaryCondition*>& active_boundary_conditions() const;
+    const std::vector<const model::Load*>& active_loads() const;
+    const model::Property* property_for(const model::Element& element) const;
+    const model::Material* material_for(const model::Property& property) const;
+
+private:
+    const model::Domain& domain_;
+    const model::AnalysisStep* step_;
+    std::vector<const model::Element*> active_elements_;
+    std::vector<const model::BoundaryCondition*> active_boundary_conditions_;
+    std::vector<const model::Load*> active_loads_;
+};
+
+} // namespace fesa::analysis

src/fesa/analysis/analysis_state.cpp

@@ -0,0 +1,124 @@
+#include <fesa/analysis/analysis_state.hpp>
+
+#include <cstddef>
+#include <stdexcept>
+#include <utility>
+
+namespace fesa::analysis {
+
+AnalysisState::AnalysisState(int total_dof_count)
+    : displacement_(vector_of(total_dof_count)),
+      velocity_(vector_of(total_dof_count)),
+      acceleration_(vector_of(total_dof_count)),
+      temperature_(vector_of(total_dof_count)),
+      external_force_(vector_of(total_dof_count)),
+      internal_force_(vector_of(total_dof_count)),
+      residual_(vector_of(total_dof_count))
+{
+}
+
+const std::vector<double>& AnalysisState::displacement() const
+{
+    return displacement_;
+}
+
+const std::vector<double>& AnalysisState::velocity() const
+{
+    return velocity_;
+}
+
+const std::vector<double>& AnalysisState::acceleration() const
+{
+    return acceleration_;
+}
+
+const std::vector<double>& AnalysisState::temperature() const
+{
+    return temperature_;
+}
+
+const std::vector<double>& AnalysisState::external_force() const
+{
+    return external_force_;
+}
+
+const std::vector<double>& AnalysisState::internal_force() const
+{
+    return internal_force_;
+}
+
+const std::vector<double>& AnalysisState::residual() const
+{
+    return residual_;
+}
+
+void AnalysisState::set_displacement(std::vector<double> values)
+{
+    assign_same_size(displacement_, std::move(values));
+}
+
+void AnalysisState::set_external_force(std::vector<double> values)
+{
+    assign_same_size(external_force_, std::move(values));
+}
+
+void AnalysisState::set_internal_force(std::vector<double> values)
+{
+    assign_same_size(internal_force_, std::move(values));
+}
+
+void AnalysisState::update_residual()
+{
+    for (std::size_t index = 0; index < residual_.size(); ++index) {
+        residual_[index] = external_force_[index] - internal_force_[index];
+    }
+}
+
+IterationState& AnalysisState::iteration_state()
+{
+    return iteration_state_;
+}
+
+const IterationState& AnalysisState::iteration_state() const
+{
+    return iteration_state_;
+}
+
+void AnalysisState::set_element_state(core::ElementId element_id, std::vector<double> state)
+{
+    for (auto& entry : element_states_) {
+        if (entry.first.value == element_id.value) {
+            entry.second = std::move(state);
+            return;
+        }
+    }
+    element_states_.push_back({element_id, std::move(state)});
+}
+
+const std::vector<double>* AnalysisState::element_state(core::ElementId element_id) const
+{
+    for (const auto& entry : element_states_) {
+        if (entry.first.value == element_id.value) {
+            return &entry.second;
+        }
+    }
+    return nullptr;
+}
+
+std::vector<double> AnalysisState::vector_of(int size)
+{
+    if (size < 0) {
+        throw std::invalid_argument("negative dof count");
+    }
+    return std::vector<double>(static_cast<std::size_t>(size), 0.0);
+}
+
+void AnalysisState::assign_same_size(std::vector<double>& target, std::vector<double> values)
+{
+    if (target.size() != values.size()) {
+        throw std::invalid_argument("vector size mismatch");
+    }
+    target = std::move(values);
+}
+
+} // namespace fesa::analysis

src/fesa/analysis/analysis_state.hpp

@@ -0,0 +1,53 @@
+#pragma once
+
+#include <fesa/core/ids.hpp>
+
+#include <vector>
+
+namespace fesa::analysis {
+
+struct IterationState {
+    double time = 0.0;
+    int increment = 0;
+    int iteration = 0;
+};
+
+class AnalysisState {
+public:
+    explicit AnalysisState(int total_dof_count);
+
+    const std::vector<double>& displacement() const;
+    const std::vector<double>& velocity() const;
+    const std::vector<double>& acceleration() const;
+    const std::vector<double>& temperature() const;
+    const std::vector<double>& external_force() const;
+    const std::vector<double>& internal_force() const;
+    const std::vector<double>& residual() const;
+
+    void set_displacement(std::vector<double> values);
+    void set_external_force(std::vector<double> values);
+    void set_internal_force(std::vector<double> values);
+    void update_residual();
+
+    IterationState& iteration_state();
+    const IterationState& iteration_state() const;
+
+    void set_element_state(core::ElementId element_id, std::vector<double> state);
+    const std::vector<double>* element_state(core::ElementId element_id) const;
+
+private:
+    static std::vector<double> vector_of(int size);
+    static void assign_same_size(std::vector<double>& target, std::vector<double> values);
+
+    std::vector<double> displacement_;
+    std::vector<double> velocity_;
+    std::vector<double> acceleration_;
+    std::vector<double> temperature_;
+    std::vector<double> external_force_;
+    std::vector<double> internal_force_;
+    std::vector<double> residual_;
+    IterationState iteration_state_;
+    std::vector<std::pair<core::ElementId, std::vector<double>>> element_states_;
+};
+
+} // namespace fesa::analysis

src/fesa/analysis/linear_static_analysis.cpp

@@ -0,0 +1,48 @@
+#include <fesa/analysis/linear_static_analysis.hpp>
+
+namespace fesa::analysis {
+
+LinearStaticAnalysis::LinearStaticAnalysis(const model::Domain& domain, core::StepId step_id)
+    : domain_(domain), step_id_(step_id)
+{
+}
+
+const AnalysisModel* LinearStaticAnalysis::analysis_model() const
+{
+    return analysis_model_.get();
+}
+
+const AnalysisState* LinearStaticAnalysis::state() const
+{
+    return state_.get();
+}
+
+void LinearStaticAnalysis::build_analysis_model()
+{
+    analysis_model_ = std::make_unique<AnalysisModel>(domain_, step_id_);
+}
+
+void LinearStaticAnalysis::build_dof_map()
+{
+    dof_manager_ = std::make_unique<fem::DofManager>();
+    for (const auto* element : analysis_model_->active_elements()) {
+        for (const auto node_id : element->node_ids()) {
+            dof_manager_->define_node_dofs(node_id, {
+                model::DofComponent::ux,
+                model::DofComponent::uy,
+                model::DofComponent::uz
+            });
+        }
+    }
+    for (const auto* boundary_condition : analysis_model_->active_boundary_conditions()) {
+        dof_manager_->apply_boundary_condition(*boundary_condition);
+    }
+    dof_manager_->number_equations();
+}
+
+void LinearStaticAnalysis::update_state()
+{
+    state_ = std::make_unique<AnalysisState>(dof_manager_->total_dof_count());
+}
+
+} // namespace fesa::analysis

src/fesa/analysis/linear_static_analysis.hpp

@@ -0,0 +1,32 @@
+#pragma once
+
+#include <fesa/analysis/analysis.hpp>
+#include <fesa/analysis/analysis_model.hpp>
+#include <fesa/analysis/analysis_state.hpp>
+#include <fesa/fem/dof_manager.hpp>
+
+#include <memory>
+
+namespace fesa::analysis {
+
+class LinearStaticAnalysis : public Analysis {
+public:
+    LinearStaticAnalysis(const model::Domain& domain, core::StepId step_id);
+
+    const AnalysisModel* analysis_model() const;
+    const AnalysisState* state() const;
+
+protected:
+    void build_analysis_model() override;
+    void build_dof_map() override;
+    void update_state() override;
+
+private:
+    const model::Domain& domain_;
+    core::StepId step_id_;
+    std::unique_ptr<AnalysisModel> analysis_model_;
+    std::unique_ptr<fem::DofManager> dof_manager_;
+    std::unique_ptr<AnalysisState> state_;
+};
+
+} // namespace fesa::analysis

src/fesa/core/diagnostic.hpp

@@ -0,0 +1,19 @@
+#pragma once
+
+#include <string>
+
+namespace fesa::core {
+
+enum class Severity {
+    info,
+    warning,
+    error
+};
+
+struct Diagnostic {
+    Severity severity;
+    std::string code;
+    std::string message;
+};
+
+} // namespace fesa::core

src/fesa/core/ids.hpp

@@ -0,0 +1,25 @@
+#pragma once
+
+namespace fesa::core {
+
+struct NodeId {
+    int value;
+};
+
+struct ElementId {
+    int value;
+};
+
+struct MaterialId {
+    int value;
+};
+
+struct PropertyId {
+    int value;
+};
+
+struct StepId {
+    int value;
+};
+
+} // namespace fesa::core

src/fesa/core/status.cpp

@@ -0,0 +1,34 @@
+#include <fesa/core/status.hpp>
+
+#include <utility>
+
+namespace fesa::core {
+
+Status Status::ok()
+{
+    return Status{};
+}
+
+Status Status::failure(Diagnostic diagnostic)
+{
+    Status status;
+    status.add(std::move(diagnostic));
+    return status;
+}
+
+bool Status::is_ok() const
+{
+    return diagnostics_.empty();
+}
+
+const std::vector<Diagnostic>& Status::diagnostics() const
+{
+    return diagnostics_;
+}
+
+void Status::add(Diagnostic diagnostic)
+{
+    diagnostics_.push_back(std::move(diagnostic));
+}
+
+} // namespace fesa::core

src/fesa/core/status.hpp

@@ -0,0 +1,22 @@
+#pragma once
+
+#include <fesa/core/diagnostic.hpp>
+
+#include <vector>
+
+namespace fesa::core {
+
+class Status {
+public:
+    static Status ok();
+    static Status failure(Diagnostic diagnostic);
+
+    bool is_ok() const;
+    const std::vector<Diagnostic>& diagnostics() const;
+    void add(Diagnostic diagnostic);
+
+private:
+    std::vector<Diagnostic> diagnostics_;
+};
+
+} // namespace fesa::core

src/fesa/elements/euler_beam_3d.cpp

@@ -0,0 +1,261 @@
+#include <fesa/elements/euler_beam_3d.hpp>
+
+#include <cmath>
+#include <cstddef>
+#include <stdexcept>
+
+namespace fesa::elements {
+namespace {
+
+constexpr int matrix_size = 12;
+constexpr double geometry_tolerance = 1.0e-12;
+
+std::size_t index(int row, int column)
+{
+    return static_cast<std::size_t>(row * matrix_size + column);
+}
+
+void require_positive_finite(double value, const char* name)
+{
+    if (!std::isfinite(value) || value <= 0.0) {
+        throw std::invalid_argument(name);
+    }
+}
+
+void validate_section(double length, const EulerBeam3DSection& section)
+{
+    require_positive_finite(length, "Euler beam length must be positive and finite");
+    require_positive_finite(section.young_modulus, "Euler beam Young modulus must be positive and finite");
+    require_positive_finite(section.shear_modulus, "Euler beam shear modulus must be positive and finite");
+    require_positive_finite(section.area, "Euler beam area must be positive and finite");
+    require_positive_finite(section.torsion_constant, "Euler beam torsion constant must be positive and finite");
+    require_positive_finite(section.second_moment_y, "Euler beam Iy must be positive and finite");
+    require_positive_finite(section.second_moment_z, "Euler beam Iz must be positive and finite");
+}
+
+void set_symmetric(Matrix12& matrix, int row, int column, double value)
+{
+    matrix[index(row, column)] = value;
+    matrix[index(column, row)] = value;
+}
+
+Vector3 subtract(const Vector3& lhs, const Vector3& rhs)
+{
+    return {
+        lhs[0] - rhs[0],
+        lhs[1] - rhs[1],
+        lhs[2] - rhs[2]
+    };
+}
+
+double dot(const Vector3& lhs, const Vector3& rhs)
+{
+    return lhs[0] * rhs[0] + lhs[1] * rhs[1] + lhs[2] * rhs[2];
+}
+
+double norm(const Vector3& vector)
+{
+    return std::sqrt(dot(vector, vector));
+}
+
+Vector3 scale(const Vector3& vector, double scalar)
+{
+    return {
+        vector[0] * scalar,
+        vector[1] * scalar,
+        vector[2] * scalar
+    };
+}
+
+Vector3 normalize(const Vector3& vector, const char* message)
+{
+    const double vector_norm = norm(vector);
+    if (!std::isfinite(vector_norm) || vector_norm <= geometry_tolerance) {
+        throw std::invalid_argument(message);
+    }
+    return scale(vector, 1.0 / vector_norm);
+}
+
+Vector3 cross(const Vector3& lhs, const Vector3& rhs)
+{
+    return {
+        lhs[1] * rhs[2] - lhs[2] * rhs[1],
+        lhs[2] * rhs[0] - lhs[0] * rhs[2],
+        lhs[0] * rhs[1] - lhs[1] * rhs[0]
+    };
+}
+
+double geometry_length(const EulerBeam3DGeometry& geometry)
+{
+    return norm(subtract(geometry.node2, geometry.node1));
+}
+
+std::array<Vector3, 3> local_basis(const EulerBeam3DGeometry& geometry)
+{
+    const auto axis = subtract(geometry.node2, geometry.node1);
+    const auto local_x = normalize(axis, "Euler beam geometry length must be positive and finite");
+
+    const double orientation_along_x = dot(geometry.orientation, local_x);
+    const auto projected_orientation = subtract(
+        geometry.orientation,
+        scale(local_x, orientation_along_x));
+    const auto local_y = normalize(
+        projected_orientation,
+        "Euler beam orientation must be nonzero and not parallel to the element axis");
+    const auto local_z = cross(local_x, local_y);
+
+    return {local_x, local_y, local_z};
+}
+
+Matrix12 transform_matrix(const EulerBeam3DGeometry& geometry)
+{
+    const auto basis = local_basis(geometry);
+    Matrix12 transform{};
+    for (int block_offset : {0, 3, 6, 9}) {
+        for (int row = 0; row < 3; ++row) {
+            for (int column = 0; column < 3; ++column) {
+                transform[index(block_offset + row, block_offset + column)] =
+                    basis[static_cast<std::size_t>(row)][static_cast<std::size_t>(column)];
+            }
+        }
+    }
+    return transform;
+}
+
+Vector12 multiply(const Matrix12& matrix, const Vector12& vector)
+{
+    Vector12 result{};
+    for (int row = 0; row < matrix_size; ++row) {
+        for (int column = 0; column < matrix_size; ++column) {
+            result[static_cast<std::size_t>(row)] +=
+                matrix[index(row, column)] * vector[static_cast<std::size_t>(column)];
+        }
+    }
+    return result;
+}
+
+Vector12 multiply_transpose(const Matrix12& matrix, const Vector12& vector)
+{
+    Vector12 result{};
+    for (int row = 0; row < matrix_size; ++row) {
+        for (int column = 0; column < matrix_size; ++column) {
+            result[static_cast<std::size_t>(row)] +=
+                matrix[index(column, row)] * vector[static_cast<std::size_t>(column)];
+        }
+    }
+    return result;
+}
+
+Matrix12 multiply(const Matrix12& lhs, const Matrix12& rhs)
+{
+    Matrix12 result{};
+    for (int row = 0; row < matrix_size; ++row) {
+        for (int column = 0; column < matrix_size; ++column) {
+            for (int inner = 0; inner < matrix_size; ++inner) {
+                result[index(row, column)] += lhs[index(row, inner)] * rhs[index(inner, column)];
+            }
+        }
+    }
+    return result;
+}
+
+Matrix12 multiply_transpose_left(const Matrix12& lhs, const Matrix12& rhs)
+{
+    Matrix12 result{};
+    for (int row = 0; row < matrix_size; ++row) {
+        for (int column = 0; column < matrix_size; ++column) {
+            for (int inner = 0; inner < matrix_size; ++inner) {
+                result[index(row, column)] += lhs[index(inner, row)] * rhs[index(inner, column)];
+            }
+        }
+    }
+    return result;
+}
+
+} // namespace
+
+Matrix12 euler_beam_3d_local_stiffness(double length, const EulerBeam3DSection& section)
+{
+    validate_section(length, section);
+
+    const double length2 = length * length;
+    const double length3 = length2 * length;
+    const double axial = section.young_modulus * section.area / length;
+    const double torsion = section.shear_modulus * section.torsion_constant / length;
+
+    const double eiy = section.young_modulus * section.second_moment_y;
+    const double cy1 = 12.0 * eiy / length3;
+    const double cy2 = 6.0 * eiy / length2;
+    const double cy3 = 4.0 * eiy / length;
+    const double cy4 = 2.0 * eiy / length;
+
+    const double eiz = section.young_modulus * section.second_moment_z;
+    const double cz1 = 12.0 * eiz / length3;
+    const double cz2 = 6.0 * eiz / length2;
+    const double cz3 = 4.0 * eiz / length;
+    const double cz4 = 2.0 * eiz / length;
+
+    Matrix12 stiffness{};
+
+    set_symmetric(stiffness, 0, 0, axial);
+    set_symmetric(stiffness, 0, 6, -axial);
+    set_symmetric(stiffness, 6, 6, axial);
+
+    set_symmetric(stiffness, 3, 3, torsion);
+    set_symmetric(stiffness, 3, 9, -torsion);
+    set_symmetric(stiffness, 9, 9, torsion);
+
+    set_symmetric(stiffness, 1, 1, cz1);
+    set_symmetric(stiffness, 1, 5, cz2);
+    set_symmetric(stiffness, 1, 7, -cz1);
+    set_symmetric(stiffness, 1, 11, cz2);
+    set_symmetric(stiffness, 5, 5, cz3);
+    set_symmetric(stiffness, 5, 7, -cz2);
+    set_symmetric(stiffness, 5, 11, cz4);
+    set_symmetric(stiffness, 7, 7, cz1);
+    set_symmetric(stiffness, 7, 11, -cz2);
+    set_symmetric(stiffness, 11, 11, cz3);
+
+    set_symmetric(stiffness, 2, 2, cy1);
+    set_symmetric(stiffness, 2, 4, -cy2);
+    set_symmetric(stiffness, 2, 8, -cy1);
+    set_symmetric(stiffness, 2, 10, -cy2);
+    set_symmetric(stiffness, 4, 4, cy3);
+    set_symmetric(stiffness, 4, 8, cy2);
+    set_symmetric(stiffness, 4, 10, cy4);
+    set_symmetric(stiffness, 8, 8, cy1);
+    set_symmetric(stiffness, 8, 10, cy2);
+    set_symmetric(stiffness, 10, 10, cy3);
+
+    return stiffness;
+}
+
+Vector12 euler_beam_3d_local_end_forces(double length,
+                                        const EulerBeam3DSection& section,
+                                        const Vector12& local_displacements)
+{
+    const auto stiffness = euler_beam_3d_local_stiffness(length, section);
+    return multiply(stiffness, local_displacements);
+}
+
+Matrix12 euler_beam_3d_global_stiffness(const EulerBeam3DGeometry& geometry,
+                                        const EulerBeam3DSection& section)
+{
+    const double length = geometry_length(geometry);
+    const auto local_stiffness = euler_beam_3d_local_stiffness(length, section);
+    const auto transform = transform_matrix(geometry);
+    return multiply_transpose_left(transform, multiply(local_stiffness, transform));
+}
+
+Vector12 euler_beam_3d_global_end_forces(const EulerBeam3DGeometry& geometry,
+                                         const EulerBeam3DSection& section,
+                                         const Vector12& global_displacements)
+{
+    const double length = geometry_length(geometry);
+    const auto transform = transform_matrix(geometry);
+    const auto local_displacements = multiply(transform, global_displacements);
+    const auto local_forces = euler_beam_3d_local_end_forces(length, section, local_displacements);
+    return multiply_transpose(transform, local_forces);
+}
+
+} // namespace fesa::elements

src/fesa/elements/euler_beam_3d.hpp

@@ -0,0 +1,39 @@
+#pragma once
+
+#include <array>
+
+namespace fesa::elements {
+
+using Vector12 = std::array<double, 12>;
+using Matrix12 = std::array<double, 144>;
+using Vector3 = std::array<double, 3>;
+
+struct EulerBeam3DSection {
+    double young_modulus;
+    double shear_modulus;
+    double area;
+    double torsion_constant;
+    double second_moment_y;
+    double second_moment_z;
+};
+
+struct EulerBeam3DGeometry {
+    Vector3 node1;
+    Vector3 node2;
+    Vector3 orientation;
+};
+
+Matrix12 euler_beam_3d_local_stiffness(double length, const EulerBeam3DSection& section);
+
+Vector12 euler_beam_3d_local_end_forces(double length,
+                                        const EulerBeam3DSection& section,
+                                        const Vector12& local_displacements);
+
+Matrix12 euler_beam_3d_global_stiffness(const EulerBeam3DGeometry& geometry,
+                                        const EulerBeam3DSection& section);
+
+Vector12 euler_beam_3d_global_end_forces(const EulerBeam3DGeometry& geometry,
+                                         const EulerBeam3DSection& section,
+                                         const Vector12& global_displacements);
+
+} // namespace fesa::elements

src/fesa/fem/dof_key.cpp

@@ -0,0 +1,10 @@
+#include <fesa/fem/dof_key.hpp>
+
+namespace fesa::fem {
+
+bool operator==(const DofKey& lhs, const DofKey& rhs)
+{
+    return lhs.node_id.value == rhs.node_id.value && lhs.component == rhs.component;
+}
+
+} // namespace fesa::fem

src/fesa/fem/dof_key.hpp

@@ -0,0 +1,15 @@
+#pragma once
+
+#include <fesa/core/ids.hpp>
+#include <fesa/model/boundary_condition.hpp>
+
+namespace fesa::fem {
+
+struct DofKey {
+    core::NodeId node_id;
+    model::DofComponent component;
+};
+
+bool operator==(const DofKey& lhs, const DofKey& rhs);
+
+} // namespace fesa::fem

src/fesa/fem/dof_manager.cpp

@@ -0,0 +1,133 @@
+#include <fesa/fem/dof_manager.hpp>
+
+#include <algorithm>
+#include <cstddef>
+#include <stdexcept>
+
+namespace fesa::fem {
+
+void DofManager::define_node_dofs(core::NodeId node_id, std::vector<model::DofComponent> components)
+{
+    for (const auto component : components) {
+        DofKey key{node_id, component};
+        if (find_record(key) == records_.end()) {
+            records_.push_back({key});
+        }
+    }
+}
+
+void DofManager::apply_boundary_condition(const model::BoundaryCondition& bc)
+{
+    auto record = find_record({bc.node_id(), bc.component()});
+    if (record == records_.end()) {
+        throw std::invalid_argument("boundary condition references undefined dof");
+    }
+    record->constrained = true;
+}
+
+void DofManager::number_equations()
+{
+    std::sort(records_.begin(), records_.end(), [](const Record& lhs, const Record& rhs) {
+        if (lhs.key.node_id.value != rhs.key.node_id.value) {
+            return lhs.key.node_id.value < rhs.key.node_id.value;
+        }
+        return static_cast<int>(lhs.key.component) < static_cast<int>(rhs.key.component);
+    });
+
+    int free_id = 0;
+    for (int equation_id = 0; equation_id < static_cast<int>(records_.size()); ++equation_id) {
+        records_[equation_id].equation_id = equation_id;
+        if (records_[equation_id].constrained) {
+            records_[equation_id].free_equation_id = std::nullopt;
+        } else {
+            records_[equation_id].free_equation_id = free_id++;
+        }
+    }
+
+    sparse_pattern_.clear();
+    for (int row = 0; row < free_id; ++row) {
+        for (int column = 0; column < free_id; ++column) {
+            sparse_pattern_.push_back({row, column});
+        }
+    }
+}
+
+int DofManager::total_dof_count() const
+{
+    return static_cast<int>(records_.size());
+}
+
+int DofManager::free_dof_count() const
+{
+    return static_cast<int>(
+        std::count_if(records_.begin(), records_.end(), [](const Record& record) {
+            return !record.constrained;
+        })
+    );
+}
+
+int DofManager::constrained_dof_count() const
+{
+    return total_dof_count() - free_dof_count();
+}
+
+bool DofManager::is_constrained(DofKey key) const
+{
+    return require_record(key).constrained;
+}
+
+int DofManager::equation_id(DofKey key) const
+{
+    return require_record(key).equation_id;
+}
+
+std::optional<int> DofManager::free_equation_id(DofKey key) const
+{
+    return require_record(key).free_equation_id;
+}
+
+std::vector<double> DofManager::expand_free_vector(const std::vector<double>& free_values) const
+{
+    if (free_values.size() != static_cast<std::size_t>(free_dof_count())) {
+        throw std::invalid_argument("free vector size does not match dof manager");
+    }
+
+    std::vector<double> full(records_.size(), 0.0);
+    for (const auto& record : records_) {
+        if (record.free_equation_id.has_value()) {
+            full[static_cast<std::size_t>(record.equation_id)] =
+                free_values[static_cast<std::size_t>(*record.free_equation_id)];
+        }
+    }
+    return full;
+}
+
+const std::vector<std::pair<int, int>>& DofManager::sparse_pattern() const
+{
+    return sparse_pattern_;
+}
+
+std::vector<DofManager::Record>::iterator DofManager::find_record(DofKey key)
+{
+    return std::find_if(records_.begin(), records_.end(), [key](const Record& record) {
+        return record.key == key;
+    });
+}
+
+std::vector<DofManager::Record>::const_iterator DofManager::find_record(DofKey key) const
+{
+    return std::find_if(records_.begin(), records_.end(), [key](const Record& record) {
+        return record.key == key;
+    });
+}
+
+const DofManager::Record& DofManager::require_record(DofKey key) const
+{
+    const auto record = find_record(key);
+    if (record == records_.end()) {
+        throw std::invalid_argument("dof is not defined");
+    }
+    return *record;
+}
+
+} // namespace fesa::fem

src/fesa/fem/dof_manager.hpp

@@ -0,0 +1,42 @@
+#pragma once
+
+#include <fesa/fem/dof_key.hpp>
+
+#include <optional>
+#include <utility>
+#include <vector>
+
+namespace fesa::fem {
+
+class DofManager {
+public:
+    void define_node_dofs(core::NodeId node_id, std::vector<model::DofComponent> components);
+    void apply_boundary_condition(const model::BoundaryCondition& bc);
+    void number_equations();
+
+    int total_dof_count() const;
+    int free_dof_count() const;
+    int constrained_dof_count() const;
+    bool is_constrained(DofKey key) const;
+    int equation_id(DofKey key) const;
+    std::optional<int> free_equation_id(DofKey key) const;
+    std::vector<double> expand_free_vector(const std::vector<double>& free_values) const;
+    const std::vector<std::pair<int, int>>& sparse_pattern() const;
+
+private:
+    struct Record {
+        DofKey key;
+        bool constrained = false;
+        int equation_id = -1;
+        std::optional<int> free_equation_id;
+    };
+
+    std::vector<Record>::iterator find_record(DofKey key);
+    std::vector<Record>::const_iterator find_record(DofKey key) const;
+    const Record& require_record(DofKey key) const;
+
+    std::vector<Record> records_;
+    std::vector<std::pair<int, int>> sparse_pattern_;
+};
+
+} // namespace fesa::fem

src/fesa/model/analysis_step.cpp

@@ -0,0 +1,42 @@
+#include <fesa/model/analysis_step.hpp>
+
+#include <utility>
+
+namespace fesa::model {
+
+AnalysisStep::AnalysisStep(core::StepId id, std::string name)
+    : id_(id), name_(std::move(name))
+{
+}
+
+core::StepId AnalysisStep::id() const
+{
+    return id_;
+}
+
+const std::string& AnalysisStep::name() const
+{
+    return name_;
+}
+
+void AnalysisStep::add_boundary_condition(BoundaryCondition bc)
+{
+    boundary_conditions_.push_back(std::move(bc));
+}
+
+void AnalysisStep::add_load(Load load)
+{
+    loads_.push_back(std::move(load));
+}
+
+const std::vector<BoundaryCondition>& AnalysisStep::boundary_conditions() const
+{
+    return boundary_conditions_;
+}
+
+const std::vector<Load>& AnalysisStep::loads() const
+{
+    return loads_;
+}
+
+} // namespace fesa::model

src/fesa/model/analysis_step.hpp

@@ -0,0 +1,29 @@
+#pragma once
+
+#include <fesa/model/boundary_condition.hpp>
+#include <fesa/model/load.hpp>
+
+#include <string>
+#include <vector>
+
+namespace fesa::model {
+
+class AnalysisStep {
+public:
+    AnalysisStep(core::StepId id, std::string name);
+
+    core::StepId id() const;
+    const std::string& name() const;
+    void add_boundary_condition(BoundaryCondition bc);
+    void add_load(Load load);
+    const std::vector<BoundaryCondition>& boundary_conditions() const;
+    const std::vector<Load>& loads() const;
+
+private:
+    core::StepId id_;
+    std::string name_;
+    std::vector<BoundaryCondition> boundary_conditions_;
+    std::vector<Load> loads_;
+};
+
+} // namespace fesa::model

src/fesa/model/boundary_condition.cpp

@@ -0,0 +1,25 @@
+#include <fesa/model/boundary_condition.hpp>
+
+namespace fesa::model {
+
+BoundaryCondition::BoundaryCondition(core::NodeId node_id, DofComponent component, double value)
+    : node_id_(node_id), component_(component), value_(value)
+{
+}
+
+core::NodeId BoundaryCondition::node_id() const
+{
+    return node_id_;
+}
+
+DofComponent BoundaryCondition::component() const
+{
+    return component_;
+}
+
+double BoundaryCondition::value() const
+{
+    return value_;
+}
+
+} // namespace fesa::model

src/fesa/model/boundary_condition.hpp

@@ -0,0 +1,31 @@
+#pragma once
+
+#include <fesa/core/ids.hpp>
+
+namespace fesa::model {
+
+enum class DofComponent {
+    ux,
+    uy,
+    uz,
+    rx,
+    ry,
+    rz,
+    temperature
+};
+
+class BoundaryCondition {
+public:
+    BoundaryCondition(core::NodeId node_id, DofComponent component, double value);
+
+    core::NodeId node_id() const;
+    DofComponent component() const;
+    double value() const;
+
+private:
+    core::NodeId node_id_;
+    DofComponent component_;
+    double value_;
+};
+
+} // namespace fesa::model

src/fesa/model/domain.cpp

@@ -0,0 +1,107 @@
+#include <fesa/model/domain.hpp>
+
+#include <utility>
+
+namespace fesa::model {
+
+void Domain::add_node(Node node)
+{
+    nodes_.push_back(std::move(node));
+}
+
+void Domain::add_element(Element element)
+{
+    elements_.push_back(std::move(element));
+}
+
+void Domain::add_material(Material material)
+{
+    materials_.push_back(std::move(material));
+}
+
+void Domain::add_property(Property property)
+{
+    properties_.push_back(std::move(property));
+}
+
+void Domain::add_step(AnalysisStep step)
+{
+    steps_.push_back(std::move(step));
+}
+
+const std::vector<Node>& Domain::nodes() const
+{
+    return nodes_;
+}
+
+const std::vector<Element>& Domain::elements() const
+{
+    return elements_;
+}
+
+const std::vector<Material>& Domain::materials() const
+{
+    return materials_;
+}
+
+const std::vector<Property>& Domain::properties() const
+{
+    return properties_;
+}
+
+const std::vector<AnalysisStep>& Domain::steps() const
+{
+    return steps_;
+}
+
+const Node* Domain::find_node(core::NodeId id) const
+{
+    for (const auto& node : nodes_) {
+        if (node.id().value == id.value) {
+            return &node;
+        }
+    }
+    return nullptr;
+}
+
+const Element* Domain::find_element(core::ElementId id) const
+{
+    for (const auto& element : elements_) {
+        if (element.id().value == id.value) {
+            return &element;
+        }
+    }
+    return nullptr;
+}
+
+const Material* Domain::find_material(core::MaterialId id) const
+{
+    for (const auto& material : materials_) {
+        if (material.id().value == id.value) {
+            return &material;
+        }
+    }
+    return nullptr;
+}
+
+const Property* Domain::find_property(core::PropertyId id) const
+{
+    for (const auto& property : properties_) {
+        if (property.id().value == id.value) {
+            return &property;
+        }
+    }
+    return nullptr;
+}
+
+const AnalysisStep* Domain::find_step(core::StepId id) const
+{
+    for (const auto& step : steps_) {
+        if (step.id().value == id.value) {
+            return &step;
+        }
+    }
+    return nullptr;
+}
+
+} // namespace fesa::model

src/fesa/model/domain.hpp

@@ -0,0 +1,41 @@
+#pragma once
+
+#include <fesa/model/analysis_step.hpp>
+#include <fesa/model/element.hpp>
+#include <fesa/model/material.hpp>
+#include <fesa/model/node.hpp>
+#include <fesa/model/property.hpp>
+
+#include <vector>
+
+namespace fesa::model {
+
+class Domain {
+public:
+    void add_node(Node node);
+    void add_element(Element element);
+    void add_material(Material material);
+    void add_property(Property property);
+    void add_step(AnalysisStep step);
+
+    const std::vector<Node>& nodes() const;
+    const std::vector<Element>& elements() const;
+    const std::vector<Material>& materials() const;
+    const std::vector<Property>& properties() const;
+    const std::vector<AnalysisStep>& steps() const;
+
+    const Node* find_node(core::NodeId id) const;
+    const Element* find_element(core::ElementId id) const;
+    const Material* find_material(core::MaterialId id) const;
+    const Property* find_property(core::PropertyId id) const;
+    const AnalysisStep* find_step(core::StepId id) const;
+
+private:
+    std::vector<Node> nodes_;
+    std::vector<Element> elements_;
+    std::vector<Material> materials_;
+    std::vector<Property> properties_;
+    std::vector<AnalysisStep> steps_;
+};
+
+} // namespace fesa::model

src/fesa/model/element.cpp

@@ -0,0 +1,38 @@
+#include <fesa/model/element.hpp>
+
+#include <utility>
+
+namespace fesa::model {
+
+Element::Element(core::ElementId id,
+                 ElementTopology topology,
+                 std::vector<core::NodeId> node_ids,
+                 core::PropertyId property_id)
+    : id_(id),
+      topology_(topology),
+      node_ids_(std::move(node_ids)),
+      property_id_(property_id)
+{
+}
+
+core::ElementId Element::id() const
+{
+    return id_;
+}
+
+ElementTopology Element::topology() const
+{
+    return topology_;
+}
+
+const std::vector<core::NodeId>& Element::node_ids() const
+{
+    return node_ids_;
+}
+
+core::PropertyId Element::property_id() const
+{
+    return property_id_;
+}
+
+} // namespace fesa::model

src/fesa/model/element.hpp

@@ -0,0 +1,35 @@
+#pragma once
+
+#include <fesa/core/ids.hpp>
+
+#include <vector>
+
+namespace fesa::model {
+
+enum class ElementTopology {
+    truss2,
+    bar2,
+    beam2,
+    unknown
+};
+
+class Element {
+public:
+    Element(core::ElementId id,
+            ElementTopology topology,
+            std::vector<core::NodeId> node_ids,
+            core::PropertyId property_id);
+
+    core::ElementId id() const;
+    ElementTopology topology() const;
+    const std::vector<core::NodeId>& node_ids() const;
+    core::PropertyId property_id() const;
+
+private:
+    core::ElementId id_;
+    ElementTopology topology_;
+    std::vector<core::NodeId> node_ids_;
+    core::PropertyId property_id_;
+};
+
+} // namespace fesa::model

src/fesa/model/load.cpp

@@ -0,0 +1,25 @@
+#include <fesa/model/load.hpp>
+
+namespace fesa::model {
+
+Load::Load(core::NodeId node_id, DofComponent component, double value)
+    : node_id_(node_id), component_(component), value_(value)
+{
+}
+
+core::NodeId Load::node_id() const
+{
+    return node_id_;
+}
+
+DofComponent Load::component() const
+{
+    return component_;
+}
+
+double Load::value() const
+{
+    return value_;
+}
+
+} // namespace fesa::model

src/fesa/model/load.hpp

@@ -0,0 +1,21 @@
+#pragma once
+
+#include <fesa/model/boundary_condition.hpp>
+
+namespace fesa::model {
+
+class Load {
+public:
+    Load(core::NodeId node_id, DofComponent component, double value);
+
+    core::NodeId node_id() const;
+    DofComponent component() const;
+    double value() const;
+
+private:
+    core::NodeId node_id_;
+    DofComponent component_;
+    double value_;
+};
+
+} // namespace fesa::model

src/fesa/model/material.cpp

@@ -0,0 +1,22 @@
+#include <fesa/model/material.hpp>
+
+#include <utility>
+
+namespace fesa::model {
+
+Material::Material(core::MaterialId id, std::string name)
+    : id_(id), name_(std::move(name))
+{
+}
+
+core::MaterialId Material::id() const
+{
+    return id_;
+}
+
+const std::string& Material::name() const
+{
+    return name_;
+}
+
+} // namespace fesa::model

src/fesa/model/material.hpp

@@ -0,0 +1,21 @@
+#pragma once
+
+#include <fesa/core/ids.hpp>
+
+#include <string>
+
+namespace fesa::model {
+
+class Material {
+public:
+    Material(core::MaterialId id, std::string name);
+
+    core::MaterialId id() const;
+    const std::string& name() const;
+
+private:
+    core::MaterialId id_;
+    std::string name_;
+};
+
+} // namespace fesa::model

src/fesa/model/node.cpp

@@ -0,0 +1,20 @@
+#include <fesa/model/node.hpp>
+
+namespace fesa::model {
+
+Node::Node(core::NodeId id, std::array<double, 3> coordinates)
+    : id_(id), coordinates_(coordinates)
+{
+}
+
+core::NodeId Node::id() const
+{
+    return id_;
+}
+
+const std::array<double, 3>& Node::coordinates() const
+{
+    return coordinates_;
+}
+
+} // namespace fesa::model

src/fesa/model/node.hpp

@@ -0,0 +1,21 @@
+#pragma once
+
+#include <fesa/core/ids.hpp>
+
+#include <array>
+
+namespace fesa::model {
+
+class Node {
+public:
+    Node(core::NodeId id, std::array<double, 3> coordinates);
+
+    core::NodeId id() const;
+    const std::array<double, 3>& coordinates() const;
+
+private:
+    core::NodeId id_;
+    std::array<double, 3> coordinates_;
+};
+
+} // namespace fesa::model

src/fesa/model/property.cpp

@@ -0,0 +1,27 @@
+#include <fesa/model/property.hpp>
+
+#include <utility>
+
+namespace fesa::model {
+
+Property::Property(core::PropertyId id, std::string name, core::MaterialId material_id)
+    : id_(id), name_(std::move(name)), material_id_(material_id)
+{
+}
+
+core::PropertyId Property::id() const
+{
+    return id_;
+}
+
+const std::string& Property::name() const
+{
+    return name_;
+}
+
+core::MaterialId Property::material_id() const
+{
+    return material_id_;
+}
+
+} // namespace fesa::model

src/fesa/model/property.hpp

@@ -0,0 +1,23 @@
+#pragma once
+
+#include <fesa/core/ids.hpp>
+
+#include <string>
+
+namespace fesa::model {
+
+class Property {
+public:
+    Property(core::PropertyId id, std::string name, core::MaterialId material_id);
+
+    core::PropertyId id() const;
+    const std::string& name() const;
+    core::MaterialId material_id() const;
+
+private:
+    core::PropertyId id_;
+    std::string name_;
+    core::MaterialId material_id_;
+};
+
+} // namespace fesa::model

src/fesa/results/results.cpp

@@ -0,0 +1,70 @@
+#include <fesa/results/results.hpp>
+
+#include <stdexcept>
+#include <utility>
+
+namespace fesa::results {
+
+ResultFrame::ResultFrame(int frame_id, double time)
+    : frame_id_(frame_id), time_(time)
+{
+}
+
+int ResultFrame::frame_id() const
+{
+    return frame_id_;
+}
+
+double ResultFrame::time() const
+{
+    return time_;
+}
+
+void ResultFrame::add_field_output(FieldOutput output)
+{
+    if (output.components.empty()) {
+        throw std::invalid_argument("field output must have components");
+    }
+    if (output.entity_ids.size() * output.components.size() != output.values.size()) {
+        throw std::invalid_argument("field output values do not match row shape");
+    }
+    field_outputs_.push_back(std::move(output));
+}
+
+void ResultFrame::add_history_output(HistoryOutput output)
+{
+    history_outputs_.push_back(std::move(output));
+}
+
+const std::vector<FieldOutput>& ResultFrame::field_outputs() const
+{
+    return field_outputs_;
+}
+
+const std::vector<HistoryOutput>& ResultFrame::history_outputs() const
+{
+    return history_outputs_;
+}
+
+ResultStep::ResultStep(std::string name)
+    : name_(std::move(name))
+{
+}
+
+const std::string& ResultStep::name() const
+{
+    return name_;
+}
+
+ResultFrame& ResultStep::add_frame(int frame_id, double time)
+{
+    frames_.push_back(ResultFrame{frame_id, time});
+    return frames_.back();
+}
+
+const std::vector<ResultFrame>& ResultStep::frames() const
+{
+    return frames_;
+}
+
+} // namespace fesa::results

src/fesa/results/results.hpp

@@ -0,0 +1,59 @@
+#pragma once
+
+#include <string>
+#include <vector>
+
+namespace fesa::results {
+
+enum class FieldLocation {
+    nodal,
+    element,
+    integration_point
+};
+
+struct FieldOutput {
+    std::string name;
+    FieldLocation location;
+    std::vector<std::string> components;
+    std::vector<int> entity_ids;
+    std::vector<double> values;
+};
+
+struct HistoryOutput {
+    std::string name;
+    std::vector<double> time;
+    std::vector<double> values;
+};
+
+class ResultFrame {
+public:
+    ResultFrame(int frame_id, double time);
+
+    int frame_id() const;
+    double time() const;
+    void add_field_output(FieldOutput output);
+    void add_history_output(HistoryOutput output);
+    const std::vector<FieldOutput>& field_outputs() const;
+    const std::vector<HistoryOutput>& history_outputs() const;
+
+private:
+    int frame_id_;
+    double time_;
+    std::vector<FieldOutput> field_outputs_;
+    std::vector<HistoryOutput> history_outputs_;
+};
+
+class ResultStep {
+public:
+    explicit ResultStep(std::string name);
+
+    const std::string& name() const;
+    ResultFrame& add_frame(int frame_id, double time);
+    const std::vector<ResultFrame>& frames() const;
+
+private:
+    std::string name_;
+    std::vector<ResultFrame> frames_;
+};
+
+} // namespace fesa::results

tests/CMakeLists.txt

@@ -0,0 +1,2 @@
+add_subdirectory(unit)
+add_subdirectory(integration)

tests/integration/CMakeLists.txt

@@ -0,0 +1,8 @@
+file(GLOB FESA_INTEGRATION_TEST_SOURCES CONFIGURE_DEPENDS "${CMAKE_CURRENT_SOURCE_DIR}/*_test.cpp")
+
+foreach(test_source IN LISTS FESA_INTEGRATION_TEST_SOURCES)
+    get_filename_component(test_name "${test_source}" NAME_WE)
+    add_executable("${test_name}" "${test_source}")
+    target_link_libraries("${test_name}" PRIVATE fesa_core)
+    add_test(NAME "${test_name}" COMMAND "${test_name}")
+endforeach()

tests/integration/solver_core_skeleton_integration_test.cpp

@@ -0,0 +1,64 @@
+#include <fesa/analysis/linear_static_analysis.hpp>
+#include <fesa/results/results.hpp>
+
+namespace {
+
+int fail()
+{
+    return 1;
+}
+
+fesa::model::Domain make_domain()
+{
+    fesa::model::Domain domain;
+    domain.add_node({fesa::core::NodeId{1}, {0.0, 0.0, 0.0}});
+    domain.add_node({fesa::core::NodeId{2}, {1.0, 0.0, 0.0}});
+    domain.add_material({fesa::core::MaterialId{3}, "steel"});
+    domain.add_property({fesa::core::PropertyId{4}, "bar", fesa::core::MaterialId{3}});
+    domain.add_element({
+        fesa::core::ElementId{5},
+        fesa::model::ElementTopology::truss2,
+        {fesa::core::NodeId{1}, fesa::core::NodeId{2}},
+        fesa::core::PropertyId{4}
+    });
+    fesa::model::AnalysisStep step{fesa::core::StepId{6}, "static-step"};
+    step.add_boundary_condition({fesa::core::NodeId{1}, fesa::model::DofComponent::ux, 0.0});
+    step.add_load({fesa::core::NodeId{2}, fesa::model::DofComponent::ux, 10.0});
+    domain.add_step(step);
+    return domain;
+}
+
+} // namespace
+
+int main()
+{
+    const auto domain = make_domain();
+    fesa::analysis::LinearStaticAnalysis analysis{domain, fesa::core::StepId{6}};
+    analysis.run();
+
+    if (analysis.analysis_model() == nullptr || analysis.state() == nullptr) {
+        return fail();
+    }
+    if (analysis.analysis_model()->active_elements().size() != 1 ||
+        analysis.analysis_model()->active_boundary_conditions().size() != 1 ||
+        analysis.analysis_model()->active_loads().size() != 1) {
+        return fail();
+    }
+
+    fesa::results::ResultStep result_step{"static-step"};
+    auto& frame = result_step.add_frame(0, 0.0);
+    frame.add_field_output({
+        "U",
+        fesa::results::FieldLocation::nodal,
+        {"ux", "uy", "uz"},
+        {1, 2},
+        {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}
+    });
+
+    if (result_step.frames().size() != 1 ||
+        result_step.frames()[0].field_outputs().size() != 1) {
+        return fail();
+    }
+
+    return 0;
+}

tests/unit/CMakeLists.txt

@@ -0,0 +1,8 @@
+file(GLOB FESA_UNIT_TEST_SOURCES CONFIGURE_DEPENDS "${CMAKE_CURRENT_SOURCE_DIR}/*_test.cpp")
+
+foreach(test_source IN LISTS FESA_UNIT_TEST_SOURCES)
+    get_filename_component(test_name "${test_source}" NAME_WE)
+    add_executable("${test_name}" "${test_source}")
+    target_link_libraries("${test_name}" PRIVATE fesa_core)
+    add_test(NAME "${test_name}" COMMAND "${test_name}")
+endforeach()

tests/unit/analysis_flow_linear_static_analysis_test.cpp

@@ -0,0 +1,37 @@
+#include <fesa/analysis/linear_static_analysis.hpp>
+
+namespace {
+
+fesa::model::Domain make_domain()
+{
+    fesa::model::Domain domain;
+    domain.add_node({fesa::core::NodeId{1}, {0.0, 0.0, 0.0}});
+    domain.add_node({fesa::core::NodeId{2}, {1.0, 0.0, 0.0}});
+    domain.add_material({fesa::core::MaterialId{3}, "steel"});
+    domain.add_property({fesa::core::PropertyId{4}, "bar", fesa::core::MaterialId{3}});
+    domain.add_element({
+        fesa::core::ElementId{5},
+        fesa::model::ElementTopology::truss2,
+        {fesa::core::NodeId{1}, fesa::core::NodeId{2}},
+        fesa::core::PropertyId{4}
+    });
+    fesa::model::AnalysisStep step{fesa::core::StepId{6}, "static"};
+    step.add_boundary_condition({fesa::core::NodeId{1}, fesa::model::DofComponent::ux, 0.0});
+    step.add_load({fesa::core::NodeId{2}, fesa::model::DofComponent::ux, 10.0});
+    domain.add_step(step);
+    return domain;
+}
+
+} // namespace
+
+int main()
+{
+    const auto domain = make_domain();
+    fesa::analysis::LinearStaticAnalysis analysis{domain, fesa::core::StepId{6}};
+    analysis.run();
+
+    if (analysis.analysis_model() == nullptr || analysis.state() == nullptr) {
+        return 1;
+    }
+    return analysis.state()->displacement().size() == 6 ? 0 : 1;
+}

tests/unit/analysis_flow_template_test.cpp

@@ -0,0 +1,44 @@
+#include <fesa/analysis/analysis.hpp>
+
+#include <string>
+#include <vector>
+
+class RecordingAnalysis : public fesa::analysis::Analysis {
+public:
+    const std::vector<std::string>& calls() const
+    {
+        return calls_;
+    }
+
+protected:
+    void initialize() override { calls_.push_back("initialize"); }
+    void build_analysis_model() override { calls_.push_back("build_analysis_model"); }
+    void build_dof_map() override { calls_.push_back("build_dof_map"); }
+    void build_sparse_pattern() override { calls_.push_back("build_sparse_pattern"); }
+    void assemble() override { calls_.push_back("assemble"); }
+    void apply_boundary_conditions() override { calls_.push_back("apply_boundary_conditions"); }
+    void solve() override { calls_.push_back("solve"); }
+    void update_state() override { calls_.push_back("update_state"); }
+    void write_results() override { calls_.push_back("write_results"); }
+
+private:
+    std::vector<std::string> calls_;
+};
+
+int main()
+{
+    RecordingAnalysis analysis;
+    analysis.run();
+    const std::vector<std::string> expected{
+        "initialize",
+        "build_analysis_model",
+        "build_dof_map",
+        "build_sparse_pattern",
+        "assemble",
+        "apply_boundary_conditions",
+        "solve",
+        "update_state",
+        "write_results"
+    };
+    return analysis.calls() == expected ? 0 : 1;
+}

tests/unit/analysis_model_view_test.cpp

@@ -0,0 +1,73 @@
+#include <fesa/analysis/analysis_model.hpp>
+
+#include <stdexcept>
+
+namespace {
+
+fesa::model::Domain make_domain()
+{
+    fesa::model::Domain domain;
+    domain.add_node({fesa::core::NodeId{1}, {0.0, 0.0, 0.0}});
+    domain.add_node({fesa::core::NodeId{2}, {1.0, 0.0, 0.0}});
+    domain.add_material({fesa::core::MaterialId{3}, "steel"});
+    domain.add_property({fesa::core::PropertyId{4}, "bar", fesa::core::MaterialId{3}});
+    domain.add_element({
+        fesa::core::ElementId{5},
+        fesa::model::ElementTopology::truss2,
+        {fesa::core::NodeId{1}, fesa::core::NodeId{2}},
+        fesa::core::PropertyId{4}
+    });
+
+    fesa::model::AnalysisStep step{fesa::core::StepId{6}, "static"};
+    step.add_boundary_condition({fesa::core::NodeId{1}, fesa::model::DofComponent::ux, 0.0});
+    step.add_load({fesa::core::NodeId{2}, fesa::model::DofComponent::ux, 10.0});
+    domain.add_step(step);
+    return domain;
+}
+
+int fail()
+{
+    return 1;
+}
+
+} // namespace
+
+int main()
+{
+    const auto domain = make_domain();
+    const fesa::analysis::AnalysisModel model{domain, fesa::core::StepId{6}};
+
+    if (&model.domain() != &domain) {
+        return fail();
+    }
+    if (&model.step() != domain.find_step(fesa::core::StepId{6})) {
+        return fail();
+    }
+    if (model.active_elements().size() != 1 ||
+        model.active_elements()[0] != domain.find_element(fesa::core::ElementId{5})) {
+        return fail();
+    }
+    if (model.active_boundary_conditions().size() != 1 ||
+        model.active_loads().size() != 1) {
+        return fail();
+    }
+
+    const auto* property = model.property_for(*model.active_elements()[0]);
+    if (property == nullptr || property != domain.find_property(fesa::core::PropertyId{4})) {
+        return fail();
+    }
+
+    const auto* material = model.material_for(*property);
+    if (material == nullptr || material != domain.find_material(fesa::core::MaterialId{3})) {
+        return fail();
+    }
+
+    try {
+        const fesa::analysis::AnalysisModel missing{domain, fesa::core::StepId{99}};
+        (void)missing;
+        return fail();
+    } catch (const std::invalid_argument&) {
+    }
+
+    return 0;
+}

tests/unit/analysis_state_vectors_test.cpp

@@ -0,0 +1,80 @@
+#include <fesa/analysis/analysis_state.hpp>
+
+#include <stdexcept>
+#include <vector>
+
+namespace {
+
+int fail()
+{
+    return 1;
+}
+
+bool all_zero(const std::vector<double>& values)
+{
+    for (const double value : values) {
+        if (value != 0.0) {
+            return false;
+        }
+    }
+    return true;
+}
+
+} // namespace
+
+int main()
+{
+    fesa::analysis::AnalysisState state{3};
+    if (state.displacement().size() != 3 ||
+        state.velocity().size() != 3 ||
+        state.acceleration().size() != 3 ||
+        state.temperature().size() != 3 ||
+        state.external_force().size() != 3 ||
+        state.internal_force().size() != 3 ||
+        state.residual().size() != 3) {
+        return fail();
+    }
+    if (!all_zero(state.displacement()) || !all_zero(state.residual())) {
+        return fail();
+    }
+
+    state.set_displacement({1.0, 2.0, 3.0});
+    if (state.displacement()[2] != 3.0) {
+        return fail();
+    }
+
+    state.set_external_force({10.0, 20.0, 30.0});
+    state.set_internal_force({1.0, 2.0, 3.0});
+    state.update_residual();
+    if (state.residual()[0] != 9.0 ||
+        state.residual()[1] != 18.0 ||
+        state.residual()[2] != 27.0) {
+        return fail();
+    }
+
+    try {
+        state.set_displacement({1.0});
+        return fail();
+    } catch (const std::invalid_argument&) {
+    }
+
+    state.iteration_state().time = 1.25;
+    state.iteration_state().increment = 2;
+    state.iteration_state().iteration = 3;
+    if (state.iteration_state().time != 1.25 ||
+        state.iteration_state().increment != 2 ||
+        state.iteration_state().iteration != 3) {
+        return fail();
+    }
+
+    state.set_element_state(fesa::core::ElementId{7}, {4.0, 5.0});
+    const auto* element_state = state.element_state(fesa::core::ElementId{7});
+    if (element_state == nullptr || element_state->size() != 2 || (*element_state)[1] != 5.0) {
+        return fail();
+    }
+    if (state.element_state(fesa::core::ElementId{8}) != nullptr) {
+        return fail();
+    }
+
+    return 0;
+}

tests/unit/core_diagnostic_test.cpp

@@ -0,0 +1,11 @@
+#include <fesa/core/diagnostic.hpp>
+
+int main()
+{
+    const fesa::core::Diagnostic diagnostic{
+        fesa::core::Severity::info,
+        "core.info",
+        "diagnostic"
+    };
+    return diagnostic.code == "core.info" ? 0 : 1;
+}

tests/unit/core_ids_test.cpp

@@ -0,0 +1,9 @@
+#include <fesa/core/ids.hpp>
+
+#include <type_traits>
+
+int main()
+{
+    static_assert(!std::is_same_v<fesa::core::NodeId, fesa::core::ElementId>);
+    return fesa::core::NodeId{7}.value == 7 ? 0 : 1;
+}

tests/unit/core_primitives_test.cpp

@@ -0,0 +1,50 @@
+#include <fesa/core/diagnostic.hpp>
+#include <fesa/core/ids.hpp>
+#include <fesa/core/status.hpp>
+
+#include <string>
+#include <type_traits>
+
+namespace {
+
+int fail()
+{
+    return 1;
+}
+
+} // namespace
+
+int main()
+{
+    static_assert(!std::is_same_v<fesa::core::NodeId, fesa::core::ElementId>);
+
+    const auto ok = fesa::core::Status::ok();
+    if (!ok.is_ok() || !ok.diagnostics().empty()) {
+        return fail();
+    }
+
+    fesa::core::Diagnostic error{
+        fesa::core::Severity::error,
+        "core.error",
+        "core failure"
+    };
+    auto failed = fesa::core::Status::failure(error);
+    if (failed.is_ok() || failed.diagnostics().size() != 1) {
+        return fail();
+    }
+    if (failed.diagnostics()[0].code != "core.error" ||
+        failed.diagnostics()[0].message != "core failure") {
+        return fail();
+    }
+
+    failed.add({fesa::core::Severity::warning, "core.warning", "check warning"});
+    if (failed.diagnostics().size() != 2) {
+        return fail();
+    }
+    if (failed.diagnostics()[0].code != "core.error" ||
+        failed.diagnostics()[1].code != "core.warning") {
+        return fail();
+    }
+
+    return 0;
+}

tests/unit/core_status_test.cpp

@@ -0,0 +1,7 @@
+#include <fesa/core/status.hpp>
+
+int main()
+{
+    const auto status = fesa::core::Status::ok();
+    return status.is_ok() ? 0 : 1;
+}

tests/unit/dof_manager_dof_key_test.cpp

@@ -0,0 +1,8 @@
+#include <fesa/fem/dof_key.hpp>
+
+int main()
+{
+    const fesa::fem::DofKey lhs{fesa::core::NodeId{1}, fesa::model::DofComponent::ux};
+    const fesa::fem::DofKey rhs{fesa::core::NodeId{1}, fesa::model::DofComponent::ux};
+    return lhs == rhs ? 0 : 1;
+}

tests/unit/dof_manager_numbering_test.cpp

@@ -0,0 +1,83 @@
+#include <fesa/fem/dof_manager.hpp>
+
+#include <stdexcept>
+#include <vector>
+
+namespace {
+
+int fail()
+{
+    return 1;
+}
+
+} // namespace
+
+int main()
+{
+    fesa::fem::DofManager dofs;
+    dofs.define_node_dofs(fesa::core::NodeId{2}, {
+        fesa::model::DofComponent::ux,
+        fesa::model::DofComponent::uy
+    });
+    dofs.define_node_dofs(fesa::core::NodeId{1}, {
+        fesa::model::DofComponent::uy,
+        fesa::model::DofComponent::ux
+    });
+    dofs.apply_boundary_condition({
+        fesa::core::NodeId{1},
+        fesa::model::DofComponent::ux,
+        0.0
+    });
+    dofs.number_equations();
+
+    const fesa::fem::DofKey n1ux{fesa::core::NodeId{1}, fesa::model::DofComponent::ux};
+    const fesa::fem::DofKey n1uy{fesa::core::NodeId{1}, fesa::model::DofComponent::uy};
+    const fesa::fem::DofKey n2ux{fesa::core::NodeId{2}, fesa::model::DofComponent::ux};
+    const fesa::fem::DofKey n2uy{fesa::core::NodeId{2}, fesa::model::DofComponent::uy};
+
+    if (dofs.total_dof_count() != 4 ||
+        dofs.constrained_dof_count() != 1 ||
+        dofs.free_dof_count() != 3) {
+        return fail();
+    }
+    if (dofs.equation_id(n1ux) != 0 ||
+        dofs.equation_id(n1uy) != 1 ||
+        dofs.equation_id(n2ux) != 2 ||
+        dofs.equation_id(n2uy) != 3) {
+        return fail();
+    }
+    if (!dofs.is_constrained(n1ux) ||
+        dofs.free_equation_id(n1ux).has_value()) {
+        return fail();
+    }
+    if (!dofs.free_equation_id(n1uy).has_value() ||
+        *dofs.free_equation_id(n1uy) != 0 ||
+        *dofs.free_equation_id(n2ux) != 1 ||
+        *dofs.free_equation_id(n2uy) != 2) {
+        return fail();
+    }
+
+    const auto full = dofs.expand_free_vector({11.0, 22.0, 33.0});
+    if (full.size() != 4 ||
+        full[0] != 0.0 ||
+        full[1] != 11.0 ||
+        full[2] != 22.0 ||
+        full[3] != 33.0) {
+        return fail();
+    }
+
+    const auto& pattern = dofs.sparse_pattern();
+    if (pattern.size() != 9 ||
+        pattern.front() != std::pair<int, int>{0, 0} ||
+        pattern.back() != std::pair<int, int>{2, 2}) {
+        return fail();
+    }
+
+    try {
+        (void)dofs.equation_id({fesa::core::NodeId{99}, fesa::model::DofComponent::ux});
+        return fail();
+    } catch (const std::invalid_argument&) {
+    }
+
+    return 0;
+}

tests/unit/euler_beam_3d_local_stiffness_test.cpp

@@ -0,0 +1,125 @@
+#include <fesa/elements/euler_beam_3d.hpp>
+
+#include <cmath>
+#include <cstddef>
+#include <stdexcept>
+
+namespace {
+
+constexpr double tolerance = 1.0e-10;
+
+bool close(double actual, double expected, double tol = tolerance)
+{
+    return std::abs(actual - expected) <= tol;
+}
+
+double entry(const fesa::elements::Matrix12& matrix, int row, int column)
+{
+    return matrix[static_cast<std::size_t>(row * 12 + column)];
+}
+
+fesa::elements::Vector12 multiply(const fesa::elements::Matrix12& matrix,
+                                  const fesa::elements::Vector12& vector)
+{
+    fesa::elements::Vector12 result{};
+    for (int row = 0; row < 12; ++row) {
+        for (int column = 0; column < 12; ++column) {
+            result[static_cast<std::size_t>(row)] +=
+                entry(matrix, row, column) * vector[static_cast<std::size_t>(column)];
+        }
+    }
+    return result;
+}
+
+bool throws_invalid_length(const fesa::elements::EulerBeam3DSection& section)
+{
+    try {
+        (void)fesa::elements::euler_beam_3d_local_stiffness(0.0, section);
+    } catch (const std::invalid_argument&) {
+        return true;
+    }
+    return false;
+}
+
+bool throws_invalid_section(fesa::elements::EulerBeam3DSection section)
+{
+    section.second_moment_z = 0.0;
+    try {
+        (void)fesa::elements::euler_beam_3d_local_stiffness(2.0, section);
+    } catch (const std::invalid_argument&) {
+        return true;
+    }
+    return false;
+}
+
+} // namespace
+
+int main()
+{
+    const double length = 2.0;
+    const fesa::elements::EulerBeam3DSection section{
+        210.0,
+        80.0,
+        3.0,
+        4.0,
+        5.0,
+        6.0
+    };
+
+    const auto stiffness = fesa::elements::euler_beam_3d_local_stiffness(length, section);
+
+    if (!close(entry(stiffness, 0, 0), 315.0) ||
+        !close(entry(stiffness, 0, 6), -315.0) ||
+        !close(entry(stiffness, 6, 6), 315.0)) {
+        return 1;
+    }
+
+    if (!close(entry(stiffness, 3, 3), 160.0) ||
+        !close(entry(stiffness, 3, 9), -160.0) ||
+        !close(entry(stiffness, 9, 9), 160.0)) {
+        return 2;
+    }
+
+    if (!close(entry(stiffness, 1, 1), 1890.0) ||
+        !close(entry(stiffness, 1, 5), 1890.0) ||
+        !close(entry(stiffness, 5, 5), 2520.0)) {
+        return 3;
+    }
+
+    if (!close(entry(stiffness, 2, 2), 1575.0) ||
+        !close(entry(stiffness, 2, 4), -1575.0) ||
+        !close(entry(stiffness, 4, 4), 2100.0)) {
+        return 4;
+    }
+
+    for (int row = 0; row < 12; ++row) {
+        for (int column = 0; column < 12; ++column) {
+            if (!close(entry(stiffness, row, column), entry(stiffness, column, row))) {
+                return 5;
+            }
+        }
+    }
+
+    fesa::elements::Vector12 displacement{};
+    displacement[0] = 0.1;
+    displacement[5] = 0.2;
+    displacement[8] = -0.05;
+
+    const auto expected_forces = multiply(stiffness, displacement);
+    const auto recovered_forces =
+        fesa::elements::euler_beam_3d_local_end_forces(length, section, displacement);
+    for (std::size_t index = 0; index < expected_forces.size(); ++index) {
+        if (!close(recovered_forces[index], expected_forces[index])) {
+            return 6;
+        }
+    }
+
+    if (!throws_invalid_length(section)) {
+        return 7;
+    }
+    if (!throws_invalid_section(section)) {
+        return 8;
+    }
+
+    return 0;
+}

tests/unit/euler_beam_3d_transform_recovery_test.cpp

@@ -0,0 +1,151 @@
+#include <fesa/elements/euler_beam_3d.hpp>
+
+#include <cmath>
+#include <cstddef>
+#include <stdexcept>
+
+namespace {
+
+constexpr double tolerance = 1.0e-9;
+
+bool close(double actual, double expected, double tol = tolerance)
+{
+    return std::abs(actual - expected) <= tol;
+}
+
+double entry(const fesa::elements::Matrix12& matrix, int row, int column)
+{
+    return matrix[static_cast<std::size_t>(row * 12 + column)];
+}
+
+fesa::elements::EulerBeam3DSection section()
+{
+    return fesa::elements::EulerBeam3DSection{
+        210.0,
+        80.0,
+        3.0,
+        4.0,
+        5.0,
+        6.0
+    };
+}
+
+bool matrices_close(const fesa::elements::Matrix12& lhs, const fesa::elements::Matrix12& rhs)
+{
+    for (std::size_t index = 0; index < lhs.size(); ++index) {
+        if (!close(lhs[index], rhs[index])) {
+            return false;
+        }
+    }
+    return true;
+}
+
+bool symmetric(const fesa::elements::Matrix12& matrix)
+{
+    for (int row = 0; row < 12; ++row) {
+        for (int column = 0; column < 12; ++column) {
+            if (!close(entry(matrix, row, column), entry(matrix, column, row))) {
+                return false;
+            }
+        }
+    }
+    return true;
+}
+
+bool throws_parallel_orientation()
+{
+    const fesa::elements::EulerBeam3DGeometry geometry{
+        {0.0, 0.0, 0.0},
+        {2.0, 0.0, 0.0},
+        {1.0, 0.0, 0.0}
+    };
+
+    try {
+        (void)fesa::elements::euler_beam_3d_global_stiffness(geometry, section());
+    } catch (const std::invalid_argument&) {
+        return true;
+    }
+    return false;
+}
+
+bool throws_zero_orientation()
+{
+    const fesa::elements::EulerBeam3DGeometry geometry{
+        {0.0, 0.0, 0.0},
+        {2.0, 0.0, 0.0},
+        {0.0, 0.0, 0.0}
+    };
+
+    try {
+        (void)fesa::elements::euler_beam_3d_global_stiffness(geometry, section());
+    } catch (const std::invalid_argument&) {
+        return true;
+    }
+    return false;
+}
+
+} // namespace
+
+int main()
+{
+    const auto beam_section = section();
+    const fesa::elements::EulerBeam3DGeometry axis_aligned{
+        {0.0, 0.0, 0.0},
+        {2.0, 0.0, 0.0},
+        {0.0, 1.0, 0.0}
+    };
+
+    const auto local = fesa::elements::euler_beam_3d_local_stiffness(2.0, beam_section);
+    const auto global_identity =
+        fesa::elements::euler_beam_3d_global_stiffness(axis_aligned, beam_section);
+    if (!matrices_close(global_identity, local)) {
+        return 1;
+    }
+
+    const fesa::elements::EulerBeam3DGeometry rotated{
+        {0.0, 0.0, 0.0},
+        {0.0, 2.0, 0.0},
+        {1.0, 0.0, 0.0}
+    };
+    const auto rotated_stiffness =
+        fesa::elements::euler_beam_3d_global_stiffness(rotated, beam_section);
+    if (!symmetric(rotated_stiffness)) {
+        return 2;
+    }
+
+    fesa::elements::Vector12 rigid_translation{};
+    rigid_translation[0] = 0.2;
+    rigid_translation[1] = -0.1;
+    rigid_translation[2] = 0.3;
+    rigid_translation[6] = 0.2;
+    rigid_translation[7] = -0.1;
+    rigid_translation[8] = 0.3;
+    const auto rigid_forces =
+        fesa::elements::euler_beam_3d_global_end_forces(rotated, beam_section, rigid_translation);
+    for (double force : rigid_forces) {
+        if (!close(force, 0.0)) {
+            return 3;
+        }
+    }
+
+    fesa::elements::Vector12 axial_extension{};
+    axial_extension[6] = 0.1;
+    const auto axial_forces = fesa::elements::euler_beam_3d_global_end_forces(
+        axis_aligned,
+        beam_section,
+        axial_extension);
+    if (!close(axial_forces[0], -31.5) ||
+        !close(axial_forces[6], 31.5) ||
+        !close(axial_forces[0] + axial_forces[6], 0.0)) {
+        return 4;
+    }
+
+    if (!throws_parallel_orientation()) {
+        return 5;
+    }
+    if (!throws_zero_orientation()) {
+        return 6;
+    }
+
+    return 0;
+}

tests/unit/harness_smoke_test.cpp

@@ -0,0 +1,7 @@
+#include <string>
+
+int main()
+{
+    const std::string project = "fesa";
+    return project.size() == 4 ? 0 : 1;
+}

tests/unit/model_analysis_step_test.cpp

@@ -0,0 +1,7 @@
+#include <fesa/model/analysis_step.hpp>
+
+int main()
+{
+    const fesa::model::AnalysisStep step{fesa::core::StepId{1}, "static"};
+    return step.name() == "static" ? 0 : 1;
+}

tests/unit/model_boundary_condition_test.cpp

@@ -0,0 +1,11 @@
+#include <fesa/model/boundary_condition.hpp>
+
+int main()
+{
+    const fesa::model::BoundaryCondition bc{
+        fesa::core::NodeId{1},
+        fesa::model::DofComponent::ux,
+        0.0
+    };
+    return bc.node_id().value == 1 ? 0 : 1;
+}

tests/unit/model_domain_test.cpp

@@ -0,0 +1,79 @@
+#include <fesa/model/domain.hpp>
+
+#include <array>
+#include <vector>
+
+namespace {
+
+int fail()
+{
+    return 1;
+}
+
+} // namespace
+
+int main()
+{
+    fesa::model::Domain domain;
+
+    domain.add_node(fesa::model::Node{fesa::core::NodeId{1}, {1.0, 2.0, 3.0}});
+    domain.add_element(fesa::model::Element{
+        fesa::core::ElementId{10},
+        fesa::model::ElementTopology::truss2,
+        {fesa::core::NodeId{1}, fesa::core::NodeId{2}},
+        fesa::core::PropertyId{20}
+    });
+    domain.add_material(fesa::model::Material{fesa::core::MaterialId{30}, "steel"});
+    domain.add_property(fesa::model::Property{
+        fesa::core::PropertyId{20},
+        "bar",
+        fesa::core::MaterialId{30}
+    });
+
+    fesa::model::AnalysisStep step{fesa::core::StepId{40}, "load"};
+    step.add_boundary_condition({fesa::core::NodeId{1}, fesa::model::DofComponent::ux, 0.0});
+    step.add_load({fesa::core::NodeId{2}, fesa::model::DofComponent::ux, 10.0});
+    domain.add_step(step);
+
+    const auto* node = domain.find_node(fesa::core::NodeId{1});
+    if (node == nullptr || node->coordinates()[2] != 3.0) {
+        return fail();
+    }
+
+    const auto* element = domain.find_element(fesa::core::ElementId{10});
+    if (element == nullptr ||
+        element->topology() != fesa::model::ElementTopology::truss2 ||
+        element->node_ids().size() != 2 ||
+        element->property_id().value != 20) {
+        return fail();
+    }
+
+    const auto* material = domain.find_material(fesa::core::MaterialId{30});
+    if (material == nullptr || material->name() != "steel") {
+        return fail();
+    }
+
+    const auto* property = domain.find_property(fesa::core::PropertyId{20});
+    if (property == nullptr ||
+        property->name() != "bar" ||
+        property->material_id().value != 30) {
+        return fail();
+    }
+
+    const auto* analysis_step = domain.find_step(fesa::core::StepId{40});
+    if (analysis_step == nullptr ||
+        analysis_step->boundary_conditions().size() != 1 ||
+        analysis_step->loads().size() != 1) {
+        return fail();
+    }
+
+    if (domain.find_node(fesa::core::NodeId{999}) != nullptr ||
+        domain.find_element(fesa::core::ElementId{999}) != nullptr ||
+        domain.find_material(fesa::core::MaterialId{999}) != nullptr ||
+        domain.find_property(fesa::core::PropertyId{999}) != nullptr ||
+        domain.find_step(fesa::core::StepId{999}) != nullptr) {
+        return fail();
+    }
+
+    return 0;
+}

tests/unit/model_element_test.cpp

@@ -0,0 +1,26 @@
+#include <fesa/model/element.hpp>
+
+int main()
+{
+    const fesa::model::Element element{
+        fesa::core::ElementId{1},
+        fesa::model::ElementTopology::bar2,
+        {fesa::core::NodeId{1}, fesa::core::NodeId{2}},
+        fesa::core::PropertyId{3}
+    };
+
+    const fesa::model::Element beam{
+        fesa::core::ElementId{10},
+        fesa::model::ElementTopology::beam2,
+        {fesa::core::NodeId{1}, fesa::core::NodeId{2}},
+        fesa::core::PropertyId{7}
+    };
+
+    if (element.topology() != fesa::model::ElementTopology::bar2) {
+        return 1;
+    }
+    if (beam.topology() != fesa::model::ElementTopology::beam2) {
+        return 1;
+    }
+    return beam.node_ids().size() == 2 ? 0 : 1;
+}

tests/unit/model_load_test.cpp

@@ -0,0 +1,11 @@
+#include <fesa/model/load.hpp>
+
+int main()
+{
+    const fesa::model::Load load{
+        fesa::core::NodeId{2},
+        fesa::model::DofComponent::ux,
+        5.0
+    };
+    return load.value() == 5.0 ? 0 : 1;
+}

tests/unit/model_material_test.cpp

@@ -0,0 +1,7 @@
+#include <fesa/model/material.hpp>
+
+int main()
+{
+    const fesa::model::Material material{fesa::core::MaterialId{1}, "steel"};
+    return material.name() == "steel" ? 0 : 1;
+}

tests/unit/model_node_test.cpp

@@ -0,0 +1,7 @@
+#include <fesa/model/node.hpp>
+
+int main()
+{
+    const fesa::model::Node node{fesa::core::NodeId{1}, {0.0, 1.0, 2.0}};
+    return node.coordinates()[1] == 1.0 ? 0 : 1;
+}

tests/unit/model_property_test.cpp

@@ -0,0 +1,11 @@
+#include <fesa/model/property.hpp>
+
+int main()
+{
+    const fesa::model::Property property{
+        fesa::core::PropertyId{2},
+        "section",
+        fesa::core::MaterialId{3}
+    };
+    return property.material_id().value == 3 ? 0 : 1;
+}

tests/unit/results_containers_test.cpp

@@ -0,0 +1,69 @@
+#include <fesa/results/results.hpp>
+
+#include <stdexcept>
+
+namespace {
+
+int fail()
+{
+    return 1;
+}
+
+} // namespace
+
+int main()
+{
+    fesa::results::ResultStep step{"static"};
+    auto& frame = step.add_frame(1, 0.0);
+    if (step.name() != "static" ||
+        step.frames().size() != 1 ||
+        frame.frame_id() != 1 ||
+        frame.time() != 0.0) {
+        return fail();
+    }
+
+    frame.add_field_output({
+        "U",
+        fesa::results::FieldLocation::nodal,
+        {"ux", "uy"},
+        {1, 2},
+        {0.0, 0.1, 1.0, 1.1}
+    });
+    if (frame.field_outputs().size() != 1 ||
+        frame.field_outputs()[0].entity_ids[1] != 2 ||
+        frame.field_outputs()[0].values[3] != 1.1) {
+        return fail();
+    }
+
+    frame.add_history_output({"load-factor", {0.0, 1.0}, {0.0, 10.0}});
+    if (frame.history_outputs().size() != 1 ||
+        frame.history_outputs()[0].values[1] != 10.0) {
+        return fail();
+    }
+
+    try {
+        frame.add_field_output({
+            "bad",
+            fesa::results::FieldLocation::nodal,
+            {},
+            {1},
+            {0.0}
+        });
+        return fail();
+    } catch (const std::invalid_argument&) {
+    }
+
+    try {
+        frame.add_field_output({
+            "bad-shape",
+            fesa::results::FieldLocation::nodal,
+            {"ux", "uy"},
+            {1},
+            {0.0}
+        });
+        return fail();
+    } catch (const std::invalid_argument&) {
+    }
+
+    return 0;
+}