modify handoff.md

CMakeLists.txt

@@ -26,6 +26,8 @@ include(CTest)
 include(cmake/FesaDependencies.cmake)
 
 add_library(fesa_core STATIC
+    src/fesa/assembly/serial_assembler.cpp
+    src/fesa/constraints/essential_bc.cpp
     src/fesa/core/version.cpp
     src/fesa/elements/beam/beam3d2.cpp
     src/fesa/fem/beam_frame.cpp
@@ -36,6 +38,7 @@ add_library(fesa_core STATIC
     src/fesa/io/abaqus/semantic_mapper.cpp
     src/fesa/model/domain.cpp
     src/fesa/model/domain_builder.cpp
+    src/fesa/solvers/linear/pardiso_linear_solver.cpp
 )
 
 target_include_directories(fesa_core
@@ -45,6 +48,7 @@ target_include_directories(fesa_core
 
 target_compile_features(fesa_core PUBLIC cxx_std_20)
 target_compile_options(fesa_core PRIVATE /W4 /permissive- /EHsc)
+target_link_libraries(fesa_core PRIVATE MKL::MKL)
 
 add_executable(fesa
     src/fesa/cli/main.cpp

docs/HANDOFF.md

@@ -2,8 +2,8 @@
 
 ## 1. 문서 목적
 
-이 문서는 `fem-and-beam-kernel` 완료 후 새 세션에서
-`equation-and-linear-solve` Phase를 바로 시작하기 위한 인수인계 기록이다.
+이 문서는 `equation-and-linear-solve` 완료 후 새 세션에서
+`results-and-pipeline` Phase를 바로 시작하기 위한 인수인계 기록이다.
 요구사항과 설계의 기준은 이 문서가 아니라 다음 파일이다.
 
 - `/AGENTS.md`
@@ -11,12 +11,13 @@
 - `/docs/ARCHITECTURE.md`
 - `/docs/ADR.md`
 - `/docs/HARNESS.md`
-- `/docs/formulation/timoshenko-beam-3d.md`
+- `/docs/HDF5_SCHEMA.md` — 다음 Phase Step 1에서 writer보다 먼저 생성할 문서
 - `/docs/superpowers/plans/2026-07-29-fesa-phase-1.md`
-- `/phases/equation-and-linear-solve/index.json`
-- `/phases/equation-and-linear-solve/step0.md`부터 `step2.md`
+- `/phases/results-and-pipeline/index.json`
+- `/phases/results-and-pipeline/step0.md`부터 `step3.md`
 
-내용이 충돌하면 위 기준 문서와 `phases/`의 현재 상태를 우선한다.
+내용이 충돌하면 `AGENTS.md`, 제품·아키텍처 문서와 `phases/`의 현재 상태를
+우선한다. 이 문서는 현재 구현과 실행환경에서 특히 놓치기 쉬운 계약을 보충한다.
 
 ## 2. 현재 저장소 상태
 
@@ -24,78 +25,126 @@
 
 - 현재 브랜치: `dev`
 - 현재 `dev` HEAD:
-  `fa48a4d9d717224ed84ec41f06c43fbc33d42780`
+  `a02024929ca62c714335dd9811eea6eab816af38`
 - `origin/dev`, `origin/HEAD`:
-  `4ee3895915accf69c3f988bf14119fab7c318ed5`
-- 로컬 `dev`는 `origin/dev`보다 10 commit 앞서고 0 commit 뒤처져 있다.
+  `d39340217421006c33bb7c61a0111d320cbaab2e`
+- 로컬 `dev`는 `origin/dev`보다 9 commit 앞서고 0 commit 뒤처져 있다.
 - 완료 Phase:
   - `solver-bootstrap`
   - `domain-and-input-skeleton`
   - `fem-and-beam-kernel`
-- 다음 Phase: `equation-and-linear-solve`
-- 다음 Step: `0 - symmetric-csr-assembly`
-- `equation-and-linear-solve`의 Step 0~2는 모두 `pending`이다.
-- `fem-and-beam-kernel`은 `dev`에 fast-forward 병합되었고 로컬
-  `feat-fem-and-beam-kernel` 브랜치는 삭제되었다.
+  - `equation-and-linear-solve`
+- 다음 Phase: `results-and-pipeline`
+- 다음 Step: `0 - result-database`
+- `results-and-pipeline`의 Step 0~3은 모두 `pending`이다.
+- `equation-and-linear-solve`는 `dev`에 fast-forward 병합되었고 로컬
+  `feat-equation-and-linear-solve` 브랜치는 삭제되었다.
 - 원격 push는 수행하지 않았다.
 
-이 문서를 갱신하기 직전 작업 트리는 clean이었다. 이 문서 변경은 사용자가 별도로
-요청하지 않는 한 커밋하지 않는다. 새 세션에서 Harness를 실행하기 전에
-`docs/HANDOFF.md` 변경을 먼저 커밋하거나 별도로 정리해야 한다. 그렇지 않으면
-executor의 feature branch와 step commit에 인수인계 문서가 섞일 수 있다.
+이 문서를 갱신하기 직전 작업 트리는 clean이었다. 현재
+`docs/HANDOFF.md` 변경은 사용자가 별도로 요청하지 않는 한 커밋하지 않는다.
+새 세션에서 Harness를 실행하기 전에 이 변경을 먼저 커밋하거나 별도로 정리해야
+한다. 그렇지 않으면 executor의 feature branch와 Step commit에 인수인계 문서가
+섞일 수 있다.
 
-## 3. 완료된 `fem-and-beam-kernel`
+로컬 `dev`의 미push 9 commit을 보존한다. 원격에 맞추기 위한 reset, 강제 checkout,
+rebase 또는 force push를 수행하지 않는다.
 
-Phase metadata는 `/phases/fem-and-beam-kernel/index.json`에 기록되어 있으며 Step
-0~3이 모두 `completed`다.
+## 3. 완료된 `equation-and-linear-solve`
 
-### FEM primitives
+Phase metadata는 `/phases/equation-and-linear-solve/index.json`에 기록되어 있으며
+Step 0~2가 모두 `completed`다.
 
-- `/include/fesa/fem/gauss_rule.hpp`
-- `/include/fesa/fem/line2_shape.hpp`
-- 1점/2점 1D Gauss rule
-- 2절점 선형 shape function, 자연좌표 derivative와 \(J=L/2\)
-- partition of unity, endpoint interpolation, 적분 정확도와 invalid-input 테스트
+주요 commit:
 
-### `DofManager`
+- `8342774` — deterministic serial symmetric CSR assembly
+- `1bf277c` — essential BC elimination과 reaction recovery
+- `e216660` — MKL PARDISO linear solver adapter
+- `2dd17be` — Phase 완료 metadata
+- `741fc9e` — 독립 review 지적과 HANDOFF 계약 수정
+- `a020249` — 완료된 Step 계약 문서 정합화
+
+독립 재검토 결과 남은 Critical, Important, Minor 항목 없이 merge-ready 판정을
+받았으며, 아래 계약이 테스트로 고정되어 있다.
+
+### 3.1 Symmetric CSR과 serial assembly
+
+관련 파일:
+
+- `/include/fesa/assembly/symmetric_csr.hpp`
+- `/include/fesa/assembly/equation_system.hpp`
+- `/include/fesa/assembly/serial_assembler.hpp`
+- `/src/fesa/assembly/serial_assembler.cpp`
+- `/tests/unit/assembly/serial_assembler_test.cpp`
+
+핵심 계약:
+
+- `SymmetricCsr`는 0-based upper triangle만 저장한다.
+- 각 row의 column index는 strictly increasing이다.
+- 연결되지 않은 절점의 자유도를 포함해 모든 row가 diagonal entry를 갖는다.
+  contribution이 없으면 값은 0이다. 따라서 이후 solver가 구조 오류가 아니라
+  singular equation system으로 진단할 수 있다.
+- sparsity pattern 생성과 numeric merge는 분리되어 있다.
+- numeric contribution은
+  `(row,column,element-origin,local-order)` 순서로 정렬·합산한다.
+- 비결합적인 `1 + 1 + 1e16` 규모의 테스트가 element origin 순서를 bit pattern으로
+  고정한다. Domain storage order로 합산하도록 바꾸면 안 된다.
+- Beam kernel failure를 0 stiffness로 대체하지 않고 assembly 오류로 전달한다.
+- `EquationSystem`은 원래 full `stiffness`와 `force`를 소유하며 constraint 처리 전
+  상태로 보존한다.
+
+이 serial assembler는 후속 `deterministic-parallel-assembly` Phase의 oracle이다.
+이번 다음 Phase에서 TBB assembly를 선행 구현하지 않는다.
+
+### 3.2 `DofManager`와 essential BC
+
+관련 파일:
 
 - `/include/fesa/fem/dof_manager.hpp`
-- internal `NodeId` 정렬에 기반한 deterministic full DOF numbering
-- 절점당 순서:
-  \(u_x,u_y,u_z,r_x,r_y,r_z\)
-- free equation mapping, prescribed value, 12-DOF Beam mapping과 full-vector 복원
-- equation ID는 `Node`나 `BeamElement`에 저장하지 않는다.
+- `/include/fesa/constraints/essential_bc.hpp`
+- `/src/fesa/constraints/essential_bc.cpp`
+- `/tests/unit/constraints/essential_bc_test.cpp`
 
-현재 `DofManager::build(const Domain&)`는 `Domain::step().prescribed_dofs`를 읽어
-구속 자유도와 prescribed 값을 이미 내부에 반영한다. 다음 Phase에서 이 상태를
-중복 소유하지 않도록 주의한다.
+핵심 계약:
 
-### Beam frame와 transformation
+- `DofManager`가 full DOF, free equation mapping과 prescribed value를 단독 소유한다.
+- full index 기반 `equation(std::size_t)`와
+  `prescribed_value(std::size_t)` query가 추가되어 있다.
+- `eliminate_essential_bcs(const EquationSystem&, const DofManager&)`에는 별도
+  prescribed 목록을 전달하지 않는다.
+- `ReducedSystem`은 reduced `stiffness`와 `force`만 소유한다. free/full mapping이나
+  prescribed full vector를 중복 저장하지 않는다.
+- full solution은 `DofManager::reconstruct_full()`로 복원한다.
+- 0과 비영 prescribed value의 RHS shift, all-constrained order 0 system, 원본 system
+  불변성이 검증되어 있다.
+- 반력은 reduced system이 아니라 원래 full system의 `r=Ku-f`에서 계산한다.
+- 입력뿐 아니라 RHS shift와 reaction 산술 결과가 NaN/Inf가 되는 경우도 성공
+  결과로 반환하지 않는다.
 
-- `/include/fesa/fem/beam_frame.hpp`
-- scale-aware normalization과 Gram-Schmidt
-- 오른손 기저 \(e_z=e_x\times e_y\)
-- 요소축과 orientation의 평행 판정에 scale-aware tolerance 적용
-- 12x12 global-to-local transformation과 회전/에너지 invariant 테스트
+`Domain::step().prescribed_dofs`를 다시 순회해 별도 constraint 상태를 만들지 않는다.
 
-### Timoshenko stiffness kernel
+### 3.3 MKL PARDISO adapter
 
-- `/include/fesa/elements/beam/beam3d2.hpp`
-- `/docs/formulation/timoshenko-beam-3d.md`
-- local DOF 순서:
-  `[u,v,w,theta_x,theta_y,theta_z]` per node
-- strain component 순서:
-  \([\epsilon,\gamma_y,\gamma_z,\kappa_x,\kappa_y,\kappa_z]\)
-- \(\gamma_y=v'-\theta_z\), \(\gamma_z=w'+\theta_y\)
-- constitutive diagonal:
-  \([EA,GA_{sy},GA_{sz},GJ,EI_y,EI_z]\)
-- 축·비틀림·굽힘은 2점, 전단은 1점 Gauss 적분
-- global stiffness는 \(T^\mathsf{T}K_\mathrm{local}T\)
-- 대칭성, 여섯 강체모드, 축, 비틀림, 양 축 굽힘과 전단, 세장비 sweep,
-  좌표 회전 invariant 테스트
+관련 파일:
 
-통합 설명은 구현처럼 upper triangle을 적분한 뒤 대칭 위치로 복사한다고
-`fa48a4d`에서 정정했다.
+- `/include/fesa/solvers/linear/linear_solver.hpp`
+- `/include/fesa/solvers/linear/pardiso_linear_solver.hpp`
+- `/src/fesa/solvers/linear/pardiso_linear_solver.cpp`
+- `/tests/unit/solvers/linear/pardiso_linear_solver_test.cpp`
+
+핵심 계약:
+
+- public `LinearSolver` 계약에는 MKL 타입이 노출되지 않는다.
+- `PardisoLinearSolver`는 noncopyable이다.
+- MKL LP64 `MKL_INT == std::int32_t`, `mtype=2`, `iparm[26]=1` matrix checker,
+  `iparm[34]=1` 0-based indexing을 사용한다.
+- analysis, factorization, solve와 release phase를 adapter 내부에서 관리한다.
+- release는 destructor fallback 외에 명시적으로 호출되어 release error도
+  `DiagnosticStage::solver`의 `solver.release_failed`로 변환된다.
+- matrix/RHS validation, singular diagnostic, repeated solve와 adapter 밖의 독립
+  상대잔차 계산이 검증되어 있다.
+- order 0 reduced system은 constraint 계층에서는 유효하지만 PARDISO 입력으로는
+  거부된다. analysis orchestration이 all-constrained case를 별도로 처리해야 한다.
 
 ## 4. 검증된 baseline과 개발환경
 
@@ -114,12 +163,12 @@ Test-Path "$env:HDF5_DIR\hdf5-config.cmake"
 Test-Path "$env:GTest_DIR\GTestConfig.cmake"
 ```
 
-2026-07-31 현재 네 package 경로가 모두 존재한다. 현재 `dev` HEAD에서 다음 baseline을
-새로 검증했다.
+2026-07-31 현재 네 package 경로가 모두 존재한다. 현재 `dev` HEAD에서 다음
+baseline을 검증했다.
 
 - MSBuild 18.8.2, MSVC v145 Debug build 성공
 - build 출력에 새 warning 없음
-- CTest 22개 중 22개 성공
+- CTest 29개 중 29개 성공
 - Harness pytest 20개 중 20개 성공
 - pytest가 실제로 20개를 수집했으므로 0-test 성공이 아님
 
@@ -138,103 +187,184 @@ CMake cache가 없거나 package 경로가 바뀐 경우에만 같은 환경 변
 cmake --fresh --preset windows-debug
 ```
 
+HDF5 command-line tool은 현재 PATH에 없지만 다음 파일은 존재한다.
+
+```text
+C:\Program Files\HDF_Group\HDF5\2.1.1\bin\h5ls.exe
+```
+
+`h5ls.exe`는 HDF5 DLL 외에 Intel `libmmd.dll`도 요구한다. HDF5 `bin`만 PATH에
+추가하면 Windows exit `0xC0000135`로 실패하므로 Step 1과 Step 3 Acceptance
+Criteria에서는 HDF5와 oneAPI `bin`을 모두 현재 세션 PATH에 추가한다. 시스템
+PATH를 영구 변경하지 않는다.
+
+```powershell
+$env:PATH = @(
+    "C:\Program Files\HDF_Group\HDF5\2.1.1\bin",
+    "C:\Program Files (x86)\Intel\oneAPI\2026.1\bin",
+    $env:PATH
+) -join ";"
+h5ls --version
+```
+
+2026-07-31 현재 이 설정에서 `h5ls: Version 2.1.1`과 exit code 0을 확인했다.
+
 ## 5. 다음 Phase 목표와 Step 순서
 
-`equation-and-linear-solve`의 독립 deliverable은 Beam contribution과 nodal load를
-full symmetric system으로 조립하고, essential BC를 소거한 뒤 MKL PARDISO로
-reduced SPD system을 푸는 것이다. 세 Step을 순서대로 실행한다.
+`results-and-pipeline`의 독립 deliverable은 nodal result semantic model과 최소 HDF5
+schema를 만들고, 이미 구현된 parser부터 linear solve까지 production 경로를
+조율해 `fesa solve ... --output ...` 수직 슬라이스를 완성하는 것이다.
 
-### Step 0 — `symmetric-csr-assembly`
+이 milestone은 입력-해석-출력 경로가 연결되었다는 뜻일 뿐 Beam의 수치 자격 완료나
+Phase 1 내부 배포 완료를 의미하지 않는다. 요소 결과 회복, 완전한 자기완결 HDF5,
+Abaqus reference qualification은 뒤의 별도 Phase에 남아 있다.
 
-- `/include/fesa/assembly/`와 `/src/fesa/assembly/`에 필요한 최소 계약을 만든다.
-- `SymmetricCsr {order,row_offsets,column_indices,values}`와
-  `EquationSystem {stiffness,force}`를 구현한다.
-- `Domain`의 Beam, material, section, node와 `DofManager::element_full_dofs()`를
-  사용해 `compute_beam3d2()`의 global contribution을 조립한다.
-- sparsity pattern 생성과 numeric merge를 분리한다.
-- `(row,column,element-origin,local-order)`의 안정된 순서로 합산한다.
-- hand-calculated 2-element system, duplicate contribution, external ID 순서 변화와
-  CSR invariant를 실패 테스트로 먼저 고정한다.
-- 이 Step의 serial assembler는 후속 병렬 조립의 oracle이다. TBB 병렬화는
-  `deterministic-parallel-assembly` Phase까지 미룬다.
+### Step 0 — `result-database`
 
-### Step 1 — `essential-bc-elimination`
+- `/include/fesa/results/`와 `/src/fesa/results/`에 HDF5와 독립적인 최소 semantic
+  result model을 만든다.
+- 현재 Step 계약의 실체는 `ResultDatabase -> ResultStep -> ResultFrame ->
+  NodalFrame`이다.
+- 이 Phase에는 전역 좌표계 nodal displacement/rotation과 reaction/moment만 담는다.
+- node ID와 6-component field의 size 일치, duplicate node, duplicate step/frame,
+  nonfinite 값을 실패 테스트로 먼저 고정한다.
+- aggregate를 그대로 공개하면서 invalid state를 나중에 검사할지, 검증된 factory를
+  둘지 구현 전에 하나로 정한다. 현재 Step의 “invalid result model을 거부” 조건을
+  실제 호출 가능한 API로 표현해야 한다.
+- element result, field/history 범용 hierarchy, velocity, acceleration, temperature를
+  미리 만들지 않는다.
 
-- `/include/fesa/constraints/`와 `/src/fesa/constraints/`에 essential-BC elimination을
-  구현한다.
-- 원래 full `EquationSystem`은 변경하지 않고 reduced symmetric system을 만든다.
-- 0과 비영 prescribed value의 RHS shift를 모두 처리한다.
-- free solution과 prescribed 값을 full vector로 복원한다.
-- 반력은 reduced matrix가 아니라 원래 시스템의 \(r=Ku-f\)에서 계산한다.
-- hand calculation으로 0/비영 지정값, all-constrained, 충돌 조건과 반력을 먼저
-  테스트한다.
+### Step 1 — `minimal-hdf5-schema`
 
-### Step 2 — `pardiso-linear-solver`
+- writer code보다 먼저 `/docs/HDF5_SCHEMA.md`에 schema `1.0.0`의 정확한 계약을
+  작성한다.
+- 최소 writer/reader adapter는 `/include/fesa/io/hdf5/`와
+  `/src/fesa/io/hdf5/`에 둔다.
+- schema/version, node origin `(part,instance,local label)`, dense ID map, 좌표,
+  connectivity, 적용 전단강성과 input/default source, nodal displacement/reaction을
+  round trip한다.
+- public reader와 `h5ls`로 writer 산출물을 다시 연다. production writer 내부 상태나
+  test-only parser로 검증하지 않는다.
+- 모든 `hid_t`와 HDF5 resource는 move-only RAII wrapper 내부에 둔다.
+- 모든 HDF5 실패를 `DiagnosticStage::results`의 diagnostic 또는 adapter 경계에서
+  포착되는 오류로 변환한다.
+- reference CSV와 아직 존재하지 않는 결과를 위한 빈 group을 추가하지 않는다.
 
-- `/include/fesa/solvers/linear/`와 `/src/fesa/solvers/linear/`에 backend 경계를
-  구현한다.
-- public `LinearSolver` 계약에는 MKL 타입을 노출하지 않는다.
-- `PardisoLinearSolver`가 analysis, factorization, solve, release 수명을 RAII로
-  관리한다.
-- `mtype=2`, LP64 index, `iparm[34]=1`, matrix checker를 사용한다.
-- 3x3 SPD, repeated solve, invalid CSR, dimension mismatch와 singular matrix를
-  실패 테스트로 먼저 고정한다.
-- 해와 상대잔차를 adapter 외부의 독립 계산으로 확인한다.
+### Step 2 — `linear-static-analysis`
 
-## 6. 다음 Phase의 핵심 계약과 확인할 설계점
+- `/include/fesa/analysis/`와 `/src/fesa/analysis/`에 orchestration만 구현한다.
+- 이미 검증된 `Domain`을 입력받으며 parser, CLI, HDF5를 호출하지 않는다.
+- 실행 순서는 `DofManager -> serial assembly -> essential BC -> PARDISO -> full
+  reconstruction -> reaction -> nodal ResultDatabase`다.
+- free equation 수가 0이면 PARDISO를 호출하지 않고 prescribed full vector와 원래
+  평형식으로 결과를 만든다.
+- nodal result의 ID 순서와 full vector의 6-DOF block 순서가 반드시 일치해야 한다.
+  `DofManager` numbering이 internal `NodeId` 정렬 기반이므로 Domain storage order를
+  묵시적으로 사용하지 않는다.
+- hand-check 가능한 한 요소 Domain으로 displacement, reaction, solver residual과
+  `Ku-f-r` 평형을 실패 테스트로 먼저 고정한다.
+- 수치 kernel, constraint 또는 solver 로직을 analysis에 복제하지 않는다.
 
-### 아키텍처 경계
+### Step 3 — `cli-pipeline-integration`
 
-- `assembly`는 sparse pattern, local-to-global mapping, contribution 정렬·병합과
-  CSR 변환만 담당한다. Beam formulation이나 PARDISO handle을 소유하지 않는다.
-- `constraints`는 essential BC와 full/reduced vector 변환만 담당한다.
-- `solvers/linear`만 MKL PARDISO API에 의존한다.
-- `core`, `model`, `fem`, `elements`에 MKL 또는 TBB API를 유입하지 않는다.
-- penalty, MPC, Lagrange multiplier, iterative solver, 병렬 조립과 analysis
-  pipeline은 이번 Phase 범위가 아니다.
+- `run_solver(const AnalysisRequest&)`가 parser, semantic mapper, Domain, analysis와
+  HDF5 writer를 application 경계에서 조율한다.
+- 현재 parser API는 `parse_deck(path)`, semantic mapping API는
+  `map_deck_to_domain(deck)`이다.
+- CLI 계약은 다음 두 가지다.
 
-### 구현 전에 명시적으로 정렬할 점
+```text
+fesa solve <model.inp> --output <results.h5>
+fesa --version
+```
 
-1. `SymmetricCsr`의 upper/lower triangle 저장 계약이 Step 파일에 아직 하나로
-   고정되어 있지 않다. Step 0에서 테스트와 주석으로 하나를 명시하고 Step 1~2와
-   일관되게 사용한다.
-2. `DofManager`는 이미 free equation mapping과 prescribed 값을 내부 소유한다.
-   Step 1 초안의 별도 `prescribed` 인수와 `ReducedSystem` 상태를 그대로 중복하지
-   말고, 기존 계약을 읽은 뒤 필요한 최소 query만 추가한다.
-3. `Domain`은 material/section span을 제공하지만 현재 ID lookup public accessor는
-   없다. 조립을 위해 mutable `Domain`이나 범용 registry를 만들지 말고, Phase 1에
-   필요한 deterministic lookup을 최소 범위로 구현한다.
-4. `BeamKernelResult`는 contribution 또는 diagnostic을 반환한다. 조립 계층이
-   실패를 무시하거나 zero matrix로 대체하지 않도록 오류 전달 계약을 테스트한다.
-5. all-constrained 문제는 reduced order가 0일 수 있다. constraint 계층에서 유효한
-   결과로 다루되, 빈 system을 PARDISO에 넘기는 정책과 혼동하지 않는다.
+- `/tests/fixtures/abaqus/minimal_cantilever.inp`는 이미 존재한다.
+- `MinimalCantileverPipeline`은 `run_solver` 또는 CLI와 public HDF5 reader만 사용한다.
+- 성공 시 node ID, finite displacement, reaction/equilibrium diagnostic과 schema path를
+  검증한다.
+- parse, semantic, equation, solver 또는 results 실패는 원래 diagnostic stage와
+  source 정보를 보존하고 CLI nonzero exit로 전달한다.
+- hierarchical fixture의 아직 미지원 keyword를 Step 3에서 우회 처리하지 않는다.
+  전체 Abaqus subset은 다음 `abaqus-subset-completion` Phase의 책임이다.
 
-위 항목은 범위를 넓히라는 의미가 아니다. 기존 코드와 Step 초안 사이의 중복 또는
-미정 계약을 구현 전에 드러내고 가장 단순한 일관된 계약을 선택하기 위한 확인
-목록이다.
+## 6. 구현 전에 명시적으로 정렬할 설계점
 
-## 7. Child sandbox의 MSBuild 실행 조건
+아래 항목은 범위를 늘리라는 의미가 아니다. 현재 Step 문서와 최종 PRD 사이의
+모호함을 구현 전에 드러내고 가장 단순한 일관된 계약을 선택하기 위한 확인 목록이다.
 
-`fem-and-beam-kernel` 실행 전 child sandbox를 실제로 조사한 결과:
+1. **Result model 유효성 API**
+   - Step 0은 invalid model을 거부하라고 하지만 제시된 타입은 public aggregate다.
+   - raw aggregate + 별도 validation과 validated factory 중 하나를 선택하고 테스트가
+     production validation 경로를 통과하게 한다.
+
+2. **HDF5 reader의 model metadata 반환 범위**
+   - Step 1의 `Hdf5ReadResult` 초안은 `ResultDatabase`만 반환하지만 round-trip 조건은
+     node origin, 좌표, connectivity, shear source까지 재검증하라고 한다.
+   - 이 metadata를 `ResultDatabase`에 억지로 넣지 말고, reader inspection model을
+     최소로 추가하거나 read result 계약을 정렬한다. `Domain`을 HDF5 API 타입으로
+     오염시키지 않는다.
+
+3. **schema `1.0.0`의 최소/최종 범위**
+   - 이번 Phase는 nodal vertical slice만 구현하고 `result-contract-completion`이
+     element 결과와 완전한 자기완결 계약을 뒤에서 채운다.
+   - 빈 미래 hierarchy는 만들지 않되, `/docs/HDF5_SCHEMA.md`에 이번 최소 required
+     dataset과 이후 additive compatibility 규칙을 분명히 구분한다.
+
+4. **deterministic nodal ordering**
+   - `DofManager` full vector는 sorted internal `NodeId` 순서다.
+   - `NodalFrame::node_ids`와 6-component displacement/reaction 배열을 같은 순서로
+     만드는 최소 query 또는 정렬 로직을 한 곳에서만 소유한다.
+
+5. **all-constrained analysis**
+   - reduced order 0은 constraint 성공이고 PARDISO invalid input이다.
+   - Step 2가 solver 호출을 생략하는 명시적 branch를 갖고 prescribed displacement와
+     full reaction을 계산한다.
+
+6. **analysis와 application 경계**
+   - Step 2 `LinearStaticAnalysis`는 Domain-to-ResultDatabase만 담당한다.
+   - Step 3 `run_solver`가 parser와 HDF5를 담당한다. ARCHITECTURE의 포괄적 lifecycle
+     설명을 이유로 HDF5 API를 analysis에 직접 넣지 않는다.
+
+## 7. 아키텍처와 범위 경계
+
+- `results`는 HDF5 API에 의존하지 않는 semantic model이다.
+- `io/hdf5`만 HDF5 C API, schema version과 resource lifetime을 안다.
+- `analysis`는 기존 production 모듈을 조율하지만 수치 kernel과 외부 API를
+  재구현하지 않는다.
+- `run_solver`와 CLI는 application 경계다. CLI parsing을 `fesa_core`의 analysis
+  객체에 넣지 않는다.
+- `core`, `model`, `fem`, `elements`는 HDF5 API에 의존하지 않는다.
+- 외부 ID와 internal dense index mapping을 혼동하지 않는다.
+- HDF5 파일에 단위 변환을 추가하지 않는다. FESA 입력과 결과는 일관 단위계를
+  전제로 한다.
+- element section result, point stress, reference CSV, TBB parallel assembly,
+  다중 Step/Instance와 미지원 Abaqus keyword를 선행 구현하지 않는다.
+- 성공 경로에 fake stiffness, fake result 또는 test-only solver를 사용하지 않는다.
+
+## 8. Child sandbox의 MSBuild 실행 조건
+
+이전 두 Harness Phase에서 child sandbox를 조사한 결과:
 
 - 앱 설치 경로의 `codex-cli 0.146.0`은 matching `codex-resources`를 찾지 못했다.
-- 완전한 standalone 배포는 다음 위치에 있다.
+- 완전한 standalone 배포는 다음 위치에 있으며 현재 `codex.exe`와
+  `codex-resources`가 모두 존재한다.
   `C:\Users\baram\.codex\packages\standalone\releases\0.146.0-x86_64-pc-windows-msvc`
-- child에서 WindowsApps 경유 PowerShell은 실패했고 Windows PowerShell 5.1은
-  정상 동작했다.
+- WindowsApps 경유 PowerShell에서는 child command가 실패했고 Windows PowerShell
+  5.1에서는 정상 동작했다.
 - `[windows] sandbox = "elevated"`에서는 MSBuild FileTracker가 access denied로
   실패했다.
-- 완전한 standalone Codex, 정리된 PATH와
+- standalone Codex, WindowsApps가 제거된 PATH와
   `[windows] sandbox = "unelevated"` 조합에서는 child MSBuild가 성공했다.
 
 현재 `C:\Users\baram\.codex\config.toml`은 원래 값인
 `[windows] sandbox = "elevated"`로 복원되어 있다. 다음 Harness 실행에서도 같은
-문제가 재현되면 사용자가 이미 선택한 1번 방식에 따라 다음 순서를 사용한다.
+문제가 재현되면 사용자가 이전에 승인한 임시 전환 방식에 따라 다음 순서를 사용한다.
 
 1. 다른 Codex 작업에 미칠 영향을 확인하고 global config 원래 값을 기록한다.
 2. Harness 실행 동안만 `[windows] sandbox = "unelevated"`로 바꾼다.
 3. 현재 PowerShell 세션의 PATH 앞에 standalone `bin`과 Windows PowerShell 5.1을
-   둔다. `Get-Command codex`와 `Get-Command powershell`로 실제 경로를 확인한다.
+   둔다. 모든 `WindowsApps` entry를 제거해야 하며 alias 디렉터리 하나만 제거하면
+   direct package 경로가 남을 수 있다.
 
 ```powershell
 $codexReleaseBin = "C:\Users\baram\.codex\packages\standalone\releases\0.146.0-x86_64-pc-windows-msvc\bin"
@@ -243,7 +373,7 @@ $filteredPath = $env:PATH -split ";" | Where-Object {
     $_ -and
     $_ -ne $codexReleaseBin -and
     $_ -ne $windowsPowerShell -and
-    $_ -notmatch "\\Microsoft\\WindowsApps$" -and
+    $_ -notmatch "WindowsApps" -and
     $_ -notmatch "\\OpenAI\\Codex\\bin$"
 }
 $env:PATH = (@($codexReleaseBin, $windowsPowerShell) + $filteredPath) -join ";"
@@ -252,8 +382,8 @@ $env:PATH = (@($codexReleaseBin, $windowsPowerShell) + $filteredPath) -join ";"
 (Get-Command powershell).Source
 ```
 
-4. 같은 세션에서 package 환경 변수와 baseline을 확인한 뒤 Harness를 실행한다.
-5. 성공·실패와 무관하게 global config를 즉시
+4. package 환경 변수와 baseline을 확인한 뒤 같은 세션에서 Harness를 실행한다.
+5. 성공·실패와 무관하게 `finally`에 해당하는 정리 단계에서 global config를 즉시
    `[windows] sandbox = "elevated"`로 복원하고 실제 값을 다시 읽어 확인한다.
 
 사용자 profile 전체나 드라이브 루트를 `--codex-add-dir`로 허용하지 않는다. 설치
@@ -261,19 +391,19 @@ $env:PATH = (@($codexReleaseBin, $windowsPowerShell) + $filteredPath) -join ";"
 release와 `codex-resources` 존재를 먼저 확인한다.
 
 이전 Harness는 모든 Step과 phase commit을 완료한 뒤 stderr reader의 CP949/UTF-8
-decode 예외를 한 번 출력했다. Phase metadata와 Git commit에는 영향이 없었고 root
-검증은 모두 통과했다. 같은 메시지가 반복되면 제품 실패로 단정하지 말고 executor
-종료 시점, phase status와 Git 이력을 함께 확인한다.
+decode 예외를 출력한 적이 있다. exit code, phase metadata와 Git commit이 정상이면
+그 메시지만으로 제품 실패로 단정하지 않는다.
 
-## 8. 새 세션 시작 절차
+## 9. 새 세션 시작 절차
 
 먼저 이 문서와 다음 Phase 파일을 모두 읽는다.
 
 ```text
-/phases/equation-and-linear-solve/index.json
-/phases/equation-and-linear-solve/step0.md
-/phases/equation-and-linear-solve/step1.md
-/phases/equation-and-linear-solve/step2.md
+/phases/results-and-pipeline/index.json
+/phases/results-and-pipeline/step0.md
+/phases/results-and-pipeline/step1.md
+/phases/results-and-pipeline/step2.md
+/phases/results-and-pipeline/step3.md
 ```
 
 그 다음 저장소 상태를 재검증한다.
@@ -283,53 +413,56 @@ git switch dev
 git status --short --branch
 git rev-parse HEAD
 git rev-parse origin/dev
+git rev-list --left-right --count origin/dev...dev
 ```
 
-로컬 `dev`의 미push 10 commit을 보존한다. `origin/dev`에 맞추기 위한 reset, checkout
-또는 rebase를 수행하지 않는다. `docs/HANDOFF.md` 변경이 남아 있으면 먼저 사용자의
-의도대로 커밋하거나 정리한다.
-
-4절의 package 환경 변수를 설정하고 baseline을 실행한다. Child MSBuild 문제가
-재현될 가능성이 있으므로 7절 조건을 적용한 뒤 다음을 실행한다.
+`docs/HANDOFF.md` 변경이 남아 있으면 먼저 사용자의 의도대로 커밋하거나 정리한다.
+4절의 package 환경 변수를 설정하고 baseline을 실행한다. 필요하면 8절의 child
+sandbox 조건을 적용한 뒤 다음을 실행한다.
 
 ```powershell
-python scripts/execute.py equation-and-linear-solve
+python scripts/execute.py results-and-pipeline
 ```
 
-executor는 `feat-equation-and-linear-solve` 브랜치를 생성하거나 checkout하고 Step
-상태와 output metadata를 기록한다. 사용자가 명시적으로 요청하지 않은 한 `--push`를
+executor는 `feat-results-and-pipeline` 브랜치를 생성하거나 checkout하고 Step 상태와
+output metadata를 기록한다. 사용자가 명시적으로 요청하지 않은 한 `--push`를
 사용하지 않는다.
 
 각 Step은 다음 순서를 지킨다.
 
 1. Step 파일의 필수 문서와 선행 구현을 모두 읽는다.
-2. 성공 기준과 CSR/constraint/backend invariant를 명시한다.
-3. 실패 테스트를 먼저 작성하고 예상한 이유로 실패함을 확인한다.
-4. 테스트를 통과시키는 최소 production code만 구현한다.
-5. focused test, 전체 CTest, Harness pytest를 실행한다.
-6. Step summary와 output metadata가 실제 결과와 일치하는지 확인한다.
-7. Phase 종료 전 전체 diff를 아키텍처·수치 계약·resource lifetime 기준으로
-   review한다.
+2. 성공 기준과 semantic/schema/lifecycle invariant를 명시한다.
+3. 모호한 계약은 6절을 기준으로 구현 전에 정렬한다.
+4. 실패 테스트를 먼저 작성하고 예상한 이유로 실패함을 확인한다.
+5. 테스트를 통과시키는 최소 production code만 구현한다.
+6. focused test, 전체 CTest와 Harness pytest를 실행한다.
+7. Step summary와 output metadata가 실제 결과와 일치하는지 확인한다.
+8. Phase 종료 전 전체 diff를 아키텍처, schema, diagnostic과 resource lifetime
+   기준으로 review한다.
 
-## 9. 다음 Phase 완료 조건
+## 10. 다음 Phase 완료 조건
 
-- `/phases/equation-and-linear-solve/index.json`의 Step 0~2가 모두 `completed`
-- `/phases/index.json`에서 `equation-and-linear-solve`가 `completed`
-- deterministic serial symmetric CSR의 구조·값·stable merge 테스트 통과
-- full load vector와 Beam global contribution의 hand calculation 일치
-- 0/비영 essential BC의 RHS shift와 full-vector 복원 검증 통과
-- 반력이 원래 \(r=Ku-f\)에서 정확히 복원됨
-- PARDISO SPD solve, repeated solve, invalid input와 singular diagnostic 검증 통과
+- `/phases/results-and-pipeline/index.json`의 Step 0~3이 모두 `completed`
+- `/phases/index.json`에서 `results-and-pipeline`이 `completed`
+- invalid nodal result model과 finite-data invariant 테스트 통과
+- `/docs/HDF5_SCHEMA.md`가 writer보다 먼저 작성되고 실제 산출물과 일치
+- public HDF5 reader round trip과 `h5ls` schema inspection 통과
+- 모든 HDF5 handle이 adapter 내부 RAII wrapper에서 해제됨
+- one-element linear static displacement, reaction, residual과 equilibrium 검증 통과
+- all-constrained analysis가 PARDISO 없이 성공
+- `fesa solve tests\fixtures\abaqus\minimal_cantilever.inp --output ...` 성공
+- CLI 실패 입력이 nonzero exit와 원래 stage/source diagnostic을 반환
 - focused test와 전체 CTest 통과
 - Harness pytest가 0개가 아닌 상태로 전체 통과
 - 새 MSVC warning 없음
-- MKL 타입과 handle이 `solvers/linear` adapter 밖으로 노출되지 않음
-- TBB 병렬 조립, analysis pipeline 또는 범위 밖 solver 기능을 선행 구현하지 않음
+- HDF5 API가 `io/hdf5` 밖의 public semantic contract로 노출되지 않음
+- element result, full Abaqus subset, TBB assembly, reference 비교를 선행 구현하지 않음
+- 이 결과를 Beam 수치 자격 완료나 Phase 1 내부 배포 완료로 표시하지 않음
 - review의 Critical/Important 항목 해결
 - 사용자 선택 전 원격 push나 `dev` 병합을 수행하지 않음
 
 새 세션의 권장 첫 요청:
 
-> `docs/HANDOFF.md`와 `phases/equation-and-linear-solve/step0.md`부터 `step2.md`를
-> 읽고 현재 baseline과 child sandbox의 MSBuild 실행 조건을 확인한 뒤
-> `equation-and-linear-solve` Phase를 시작해주세요.
+> `docs/HANDOFF.md`와 `phases/results-and-pipeline/step0.md`부터 `step3.md`를 읽고
+> 현재 baseline, HDF5 도구 경로와 child sandbox의 MSBuild 실행 조건을 확인한 뒤
+> `results-and-pipeline` Phase를 시작해주세요.

include/fesa/assembly/equation_system.hpp

@@ -0,0 +1,14 @@
+#pragma once
+
+#include <vector>
+
+#include <fesa/assembly/symmetric_csr.hpp>
+
+namespace fesa {
+
+struct EquationSystem final {
+    SymmetricCsr stiffness;
+    std::vector<double> force;
+};
+
+}  // namespace fesa

include/fesa/assembly/serial_assembler.hpp

@@ -0,0 +1,13 @@
+#pragma once
+
+#include <fesa/assembly/equation_system.hpp>
+#include <fesa/fem/dof_manager.hpp>
+#include <fesa/model/domain.hpp>
+
+namespace fesa {
+
+[[nodiscard]] EquationSystem assemble_serial(
+    const Domain& domain,
+    const DofManager& dofs);
+
+}  // namespace fesa

include/fesa/assembly/symmetric_csr.hpp

@@ -0,0 +1,18 @@
+#pragma once
+
+#include <cstddef>
+#include <cstdint>
+#include <vector>
+
+namespace fesa {
+
+// Zero-based CSR containing only the upper triangle. Column indices are
+// strictly increasing within each row.
+struct SymmetricCsr final {
+    std::size_t order;
+    std::vector<std::int32_t> row_offsets;
+    std::vector<std::int32_t> column_indices;
+    std::vector<double> values;
+};
+
+}  // namespace fesa

include/fesa/constraints/essential_bc.hpp

@@ -0,0 +1,31 @@
+#pragma once
+
+#include <optional>
+#include <span>
+#include <vector>
+
+#include <fesa/assembly/equation_system.hpp>
+#include <fesa/core/diagnostic.hpp>
+#include <fesa/fem/dof_manager.hpp>
+
+namespace fesa {
+
+struct ReducedSystem final {
+    SymmetricCsr stiffness;
+    std::vector<double> force;
+};
+
+struct ConstraintResult final {
+    std::optional<ReducedSystem> reduced_system;
+    std::vector<Diagnostic> diagnostics;
+};
+
+[[nodiscard]] ConstraintResult eliminate_essential_bcs(
+    const EquationSystem& original,
+    const DofManager& dofs);
+
+[[nodiscard]] std::vector<double> recover_reaction(
+    const EquationSystem& original,
+    std::span<const double> full_displacement);
+
+}  // namespace fesa

include/fesa/fem/dof_manager.hpp

@@ -34,13 +34,19 @@ public:
     [[nodiscard]] std::size_t free_equation_count() const noexcept;
     [[nodiscard]] std::optional<std::size_t> equation(
         DofAddress address) const;
+    [[nodiscard]] std::optional<std::size_t> equation(
+        std::size_t full_dof) const;
+    [[nodiscard]] std::optional<double> prescribed_value(
+        DofAddress address) const;
+    [[nodiscard]] std::optional<double> prescribed_value(
+        std::size_t full_dof) const;
+    [[nodiscard]] std::size_t full_dof(DofAddress address) const;
     [[nodiscard]] std::array<std::size_t, 12> element_full_dofs(
         const BeamElement& element) const;
     [[nodiscard]] std::vector<double> reconstruct_full(
         std::span<const double> reduced) const;
 
 private:
-    [[nodiscard]] std::size_t full_dof(DofAddress address) const;
     [[nodiscard]] std::size_t node_full_dof_base(NodeId node) const;
 
     std::map<std::int64_t, std::size_t> node_full_dof_bases_;

include/fesa/model/domain.hpp

@@ -57,6 +57,8 @@ public:
 
     [[nodiscard]] const Node& node(NodeId id) const;
     [[nodiscard]] const Node& node(const EntityOrigin& origin) const;
+    [[nodiscard]] const IsotropicElastic& material(MaterialId id) const;
+    [[nodiscard]] const BeamSection& section(SectionId id) const;
 
 private:
     friend class DomainBuilder;

include/fesa/solvers/linear/linear_solver.hpp

@@ -0,0 +1,26 @@
+#pragma once
+
+#include <span>
+#include <vector>
+
+#include <fesa/assembly/symmetric_csr.hpp>
+#include <fesa/core/diagnostic.hpp>
+
+namespace fesa {
+
+struct LinearSolveResult final {
+    std::vector<double> solution;
+    double relative_residual{};
+    std::vector<Diagnostic> diagnostics;
+};
+
+class LinearSolver {
+public:
+    virtual ~LinearSolver() = default;
+
+    [[nodiscard]] virtual LinearSolveResult solve(
+        const SymmetricCsr& matrix,
+        std::span<const double> rhs) = 0;
+};
+
+}  // namespace fesa

include/fesa/solvers/linear/pardiso_linear_solver.hpp

@@ -0,0 +1,22 @@
+#pragma once
+
+#include <span>
+
+#include <fesa/solvers/linear/linear_solver.hpp>
+
+namespace fesa {
+
+class PardisoLinearSolver final : public LinearSolver {
+public:
+    PardisoLinearSolver();
+    ~PardisoLinearSolver() override;
+    PardisoLinearSolver(const PardisoLinearSolver&) = delete;
+    PardisoLinearSolver& operator=(
+        const PardisoLinearSolver&) = delete;
+
+    [[nodiscard]] LinearSolveResult solve(
+        const SymmetricCsr& matrix,
+        std::span<const double> rhs) override;
+};
+
+}  // namespace fesa

phases/equation-and-linear-solve/index.json

@@ -5,17 +5,28 @@
     {
       "step": 0,
       "name": "symmetric-csr-assembly",
-      "status": "pending"
+      "status": "completed",
+      "started_at": "2026-07-31T15:19:35+0900",
+      "summary": "Added upper-triangle symmetric CSR with every-row diagonals and deterministic serial Beam stiffness/load assembly with stable floating-point merge tests.",
+      "completed_at": "2026-07-31T15:45:02+0900"
     },
     {
       "step": 1,
       "name": "essential-bc-elimination",
-      "status": "pending"
+      "status": "completed",
+      "started_at": "2026-07-31T15:45:03+0900",
+      "summary": "Added DofManager-owned essential BC elimination with zero/nonzero RHS shifts, all-constrained handling, full reconstruction, and finite original-system reaction recovery.",
+      "completed_at": "2026-07-31T16:00:44+0900"
     },
     {
       "step": 2,
       "name": "pardiso-linear-solver",
-      "status": "pending"
+      "status": "completed",
+      "started_at": "2026-07-31T16:00:45+0900",
+      "summary": "Added a noncopyable LP64 MKL PARDISO SPD solver with explicit diagnostic-aware RAII release, CSR validation, residual diagnostics, and regression tests.",
+      "completed_at": "2026-07-31T16:16:28+0900"
     }
-  ]
+  ],
+  "created_at": "2026-07-31T15:19:35+0900",
+  "completed_at": "2026-07-31T16:16:28+0900"
 }

phases/equation-and-linear-solve/step0.md

@@ -33,8 +33,9 @@ struct EquationSystem final {
 
 - sparsity pattern builder와 numeric contribution merge를 분리한다.
 - `(row,column,element-origin,local-order)`의 안정된 순서로 합산한다.
+- 연결되지 않은 자유도를 포함해 모든 CSR row에 diagonal entry를 보존한다.
 - hand-calculated 2-element system, duplicate contribution, external ID 순서 변화,
-  CSR invariant를 실패 테스트로 먼저 작성한다.
+  비결합적 부동소수점 합산 순서와 CSR invariant를 실패 테스트로 먼저 작성한다.
 
 ## Acceptance Criteria
 

phases/equation-and-linear-solve/step1.md

@@ -9,30 +9,29 @@
 - `/include/fesa/assembly/symmetric_csr.hpp`
 - `/include/fesa/assembly/equation_system.hpp`
 - `/include/fesa/fem/dof_manager.hpp`
-- `/include/fesa/model/step_definition.hpp`
 
 ## 작업
 
-0과 비영 지정변위를 지원하는 essential-BC elimination과 full-vector 복원을 구현한다.
+`DofManager`가 소유한 equation mapping과 지정변위로 essential-BC elimination을
+수행하고, 기존 `DofManager::reconstruct_full()`로 full vector를 복원한다.
 
 ```cpp
 struct ReducedSystem final {
     SymmetricCsr stiffness;
     std::vector<double> force;
-    std::vector<std::size_t> free_to_full;
-    std::vector<double> prescribed_full;
 };
 [[nodiscard]] ConstraintResult eliminate_essential_bcs(
     const EquationSystem& original,
-    const DofManager& dofs,
-    std::span<const PrescribedDof> prescribed);
+    const DofManager& dofs);
 [[nodiscard]] std::vector<double> recover_reaction(
     const EquationSystem& original,
     std::span<const double> full_displacement);
 ```
 
 - 작은 hand calculation으로 RHS shift, 0/비영 prescribed value, all constrained,
-  충돌 조건, \(r=Ku-f\) 반력 복원을 먼저 테스트한다.
+  \(r=Ku-f\) 반력 복원을 먼저 테스트한다.
+- 별도 prescribed 인수나 full/free mapping 상태를 중복하지 않는다.
+- 입력과 산술 결과의 NaN/Inf를 성공 결과로 반환하지 않는다.
 
 ## Acceptance Criteria
 
@@ -46,7 +45,7 @@ ctest --preset windows-debug --output-on-failure
 
 1. 비영 지정값 테스트의 실패를 먼저 확인한다.
 2. 원래 EquationSystem을 보존한 채 reduced system을 생성한다.
-3. 복원 변위와 원래 평형식 반력을 assertion한다.
+3. `DofManager` 복원 변위와 원래 평형식 반력을 assertion한다.
 4. 전체 테스트와 index를 갱신한다.
 
 ## 금지사항

phases/equation-and-linear-solve/step2.md

@@ -32,6 +32,8 @@ class PardisoLinearSolver final : public LinearSolver {
 public:
     PardisoLinearSolver();
     ~PardisoLinearSolver() override;
+    PardisoLinearSolver(const PardisoLinearSolver&) = delete;
+    PardisoLinearSolver& operator=(const PardisoLinearSolver&) = delete;
     [[nodiscard]] LinearSolveResult solve(
         const SymmetricCsr&,
         std::span<const double>) override;
@@ -42,6 +44,7 @@ public:
   테스트한다.
 - `mtype=2`, LP64 index, `iparm[34]=1`, matrix checker, analysis/factor/solve/release
   phase를 사용한다.
+- analysis/factor/solve/release의 모든 MKL 오류를 solver diagnostic으로 변환한다.
 
 ## Acceptance Criteria
 

phases/index.json

@@ -17,7 +17,8 @@
     },
     {
       "dir": "equation-and-linear-solve",
-      "status": "pending"
+      "status": "completed",
+      "completed_at": "2026-07-31T16:16:28+0900"
     },
     {
       "dir": "results-and-pipeline",

src/fesa/assembly/serial_assembler.cpp

@@ -0,0 +1,244 @@
+#include <fesa/assembly/serial_assembler.hpp>
+
+#include <algorithm>
+#include <array>
+#include <cstddef>
+#include <cstdint>
+#include <iterator>
+#include <limits>
+#include <stdexcept>
+#include <string>
+#include <tuple>
+#include <utility>
+#include <vector>
+
+#include <fesa/elements/beam/beam3d2.hpp>
+
+namespace fesa {
+namespace {
+
+struct NumericContribution final {
+    std::size_t row;
+    std::size_t column;
+    const EntityOrigin* element_origin;
+    std::size_t local_order;
+    double value;
+};
+
+std::int32_t csr_index(const std::size_t value) {
+    if (value >
+        static_cast<std::size_t>(
+            std::numeric_limits<std::int32_t>::max())) {
+        throw std::overflow_error{
+            "Symmetric CSR exceeds the 32-bit index range."};
+    }
+    return static_cast<std::int32_t>(value);
+}
+
+auto origin_key(const EntityOrigin& origin) {
+    return std::tie(
+        origin.instance_name,
+        origin.local_label,
+        origin.part_name);
+}
+
+SymmetricCsr build_sparsity_pattern(
+    const Domain& domain,
+    const DofManager& dofs) {
+    using Coordinate = std::pair<std::size_t, std::size_t>;
+    std::vector<Coordinate> coordinates;
+    coordinates.reserve(
+        dofs.full_dof_count() +
+        domain.beam_elements().size() * 78);
+
+    for (std::size_t full = 0; full < dofs.full_dof_count(); ++full) {
+        coordinates.emplace_back(full, full);
+    }
+
+    for (const BeamElement& element : domain.beam_elements()) {
+        const std::array<std::size_t, 12> full_dofs =
+            dofs.element_full_dofs(element);
+        for (std::size_t local_row = 0; local_row < full_dofs.size();
+             ++local_row) {
+            for (std::size_t local_column = local_row;
+                 local_column < full_dofs.size();
+                 ++local_column) {
+                coordinates.emplace_back(
+                    std::min(
+                        full_dofs[local_row],
+                        full_dofs[local_column]),
+                    std::max(
+                        full_dofs[local_row],
+                        full_dofs[local_column]));
+            }
+        }
+    }
+
+    std::ranges::sort(coordinates);
+    coordinates.erase(
+        std::ranges::unique(coordinates).begin(),
+        coordinates.end());
+
+    SymmetricCsr pattern{
+        dofs.full_dof_count(),
+        std::vector<std::int32_t>(dofs.full_dof_count() + 1, 0),
+        {},
+        {},
+    };
+    pattern.column_indices.reserve(coordinates.size());
+    for (const auto [row, column] : coordinates) {
+        if (row >= pattern.order || column >= pattern.order) {
+            throw std::invalid_argument{
+                "DofManager element mapping exceeds the full system order."};
+        }
+        ++pattern.row_offsets[row + 1];
+        pattern.column_indices.push_back(csr_index(column));
+    }
+    for (std::size_t row = 0; row < pattern.order; ++row) {
+        const std::size_t offset =
+            static_cast<std::size_t>(pattern.row_offsets[row]) +
+            static_cast<std::size_t>(pattern.row_offsets[row + 1]);
+        pattern.row_offsets[row + 1] = csr_index(offset);
+    }
+    pattern.values.resize(pattern.column_indices.size(), 0.0);
+    return pattern;
+}
+
+std::runtime_error kernel_error(
+    const BeamElement& element,
+    const BeamKernelResult& result) {
+    std::string message =
+        "Beam element " + std::to_string(element.origin.local_label) +
+        " kernel failed";
+    for (const Diagnostic& diagnostic : result.diagnostics) {
+        message += ": " + diagnostic.code + " - " + diagnostic.message;
+    }
+    return std::runtime_error{std::move(message)};
+}
+
+std::vector<NumericContribution> collect_numeric_contributions(
+    const Domain& domain,
+    const DofManager& dofs) {
+    std::vector<NumericContribution> contributions;
+    contributions.reserve(domain.beam_elements().size() * 78);
+
+    for (const BeamElement& element : domain.beam_elements()) {
+        const BeamKernelResult result = compute_beam3d2({
+            {
+                domain.node(element.nodes[0]).position,
+                domain.node(element.nodes[1]).position,
+            },
+            domain.material(element.material),
+            domain.section(element.section),
+        });
+        if (!result.contribution.has_value()) {
+            throw kernel_error(element, result);
+        }
+
+        const std::array<std::size_t, 12> full_dofs =
+            dofs.element_full_dofs(element);
+        for (std::size_t local_row = 0; local_row < full_dofs.size();
+             ++local_row) {
+            for (std::size_t local_column = local_row;
+                 local_column < full_dofs.size();
+                 ++local_column) {
+                contributions.push_back({
+                    std::min(
+                        full_dofs[local_row],
+                        full_dofs[local_column]),
+                    std::max(
+                        full_dofs[local_row],
+                        full_dofs[local_column]),
+                    &element.origin,
+                    local_row * full_dofs.size() + local_column,
+                    result.contribution
+                        ->global_stiffness[local_row][local_column],
+                });
+            }
+        }
+    }
+
+    std::ranges::sort(
+        contributions,
+        [](const NumericContribution& left,
+           const NumericContribution& right) {
+            return std::tuple{
+                       left.row,
+                       left.column,
+                       origin_key(*left.element_origin),
+                       left.local_order} <
+                   std::tuple{
+                       right.row,
+                       right.column,
+                       origin_key(*right.element_origin),
+                       right.local_order};
+        });
+    return contributions;
+}
+
+void merge_numeric_contributions(
+    SymmetricCsr& matrix,
+    const std::vector<NumericContribution>& contributions) {
+    std::size_t contribution_index = 0;
+    while (contribution_index < contributions.size()) {
+        const NumericContribution& first =
+            contributions[contribution_index];
+        double value = 0.0;
+        do {
+            value += contributions[contribution_index].value;
+            ++contribution_index;
+        } while (
+            contribution_index < contributions.size() &&
+            contributions[contribution_index].row == first.row &&
+            contributions[contribution_index].column == first.column);
+
+        const auto row_begin =
+            matrix.column_indices.begin() +
+            matrix.row_offsets[first.row];
+        const auto row_end =
+            matrix.column_indices.begin() +
+            matrix.row_offsets[first.row + 1];
+        const auto entry = std::lower_bound(
+            row_begin,
+            row_end,
+            csr_index(first.column));
+        if (entry == row_end || *entry != csr_index(first.column)) {
+            throw std::logic_error{
+                "Numeric contribution is absent from the CSR pattern."};
+        }
+        matrix.values[static_cast<std::size_t>(
+            std::distance(matrix.column_indices.begin(), entry))] = value;
+    }
+}
+
+std::vector<double> assemble_force(
+    const Domain& domain,
+    const DofManager& dofs) {
+    std::vector<double> force(dofs.full_dof_count(), 0.0);
+    for (const NodalLoad& load : domain.step().nodal_loads) {
+        for (std::size_t component = 0; component < load.values.size();
+             ++component) {
+            const auto dof = static_cast<NodeDof>(component);
+            force[dofs.full_dof({load.node, dof})] +=
+                load.values[component];
+        }
+    }
+    return force;
+}
+
+}  // namespace
+
+EquationSystem assemble_serial(
+    const Domain& domain,
+    const DofManager& dofs) {
+    SymmetricCsr stiffness = build_sparsity_pattern(domain, dofs);
+    const std::vector<NumericContribution> contributions =
+        collect_numeric_contributions(domain, dofs);
+    merge_numeric_contributions(stiffness, contributions);
+    return {
+        std::move(stiffness),
+        assemble_force(domain, dofs),
+    };
+}
+
+}  // namespace fesa

src/fesa/constraints/essential_bc.cpp

@@ -0,0 +1,254 @@
+#include <fesa/constraints/essential_bc.hpp>
+
+#include <algorithm>
+#include <cmath>
+#include <cstddef>
+#include <cstdint>
+#include <limits>
+#include <optional>
+#include <stdexcept>
+#include <string>
+#include <utility>
+#include <vector>
+
+namespace fesa {
+namespace {
+
+std::optional<std::string> validate_system(
+    const EquationSystem& system) {
+    const SymmetricCsr& matrix = system.stiffness;
+    if (system.force.size() != matrix.order) {
+        return "Force size must equal the matrix order.";
+    }
+    if (matrix.order >
+        static_cast<std::size_t>(
+            std::numeric_limits<std::int32_t>::max())) {
+        return "Matrix order exceeds the 32-bit CSR index range.";
+    }
+    if (matrix.row_offsets.size() != matrix.order + 1 ||
+        matrix.row_offsets.empty() ||
+        matrix.row_offsets.front() != 0) {
+        return "CSR row offsets must contain order + 1 entries "
+               "starting at zero.";
+    }
+    if (matrix.column_indices.size() != matrix.values.size()) {
+        return "CSR column and value counts must match.";
+    }
+    if (!std::ranges::all_of(
+            matrix.values, [](const double value) {
+                return std::isfinite(value);
+            }) ||
+        !std::ranges::all_of(
+            system.force, [](const double value) {
+                return std::isfinite(value);
+            })) {
+        return "Matrix and force values must be finite.";
+    }
+
+    std::int32_t previous_offset = 0;
+    for (const std::int32_t offset : matrix.row_offsets) {
+        if (offset < previous_offset || offset < 0 ||
+            static_cast<std::size_t>(offset) >
+                matrix.column_indices.size()) {
+            return "CSR row offsets must be nondecreasing and in range.";
+        }
+        previous_offset = offset;
+    }
+    if (static_cast<std::size_t>(matrix.row_offsets.back()) !=
+        matrix.column_indices.size()) {
+        return "The final CSR row offset must equal the entry count.";
+    }
+
+    for (std::size_t row = 0; row < matrix.order; ++row) {
+        std::int32_t previous_column = -1;
+        const std::size_t begin =
+            static_cast<std::size_t>(matrix.row_offsets[row]);
+        const std::size_t end =
+            static_cast<std::size_t>(matrix.row_offsets[row + 1]);
+        for (std::size_t entry = begin; entry < end; ++entry) {
+            const std::int32_t column = matrix.column_indices[entry];
+            if (column < static_cast<std::int32_t>(row) ||
+                column >= static_cast<std::int32_t>(matrix.order) ||
+                column <= previous_column) {
+                return "CSR rows must contain sorted unique "
+                       "upper-triangle columns.";
+            }
+            previous_column = column;
+        }
+    }
+    return std::nullopt;
+}
+
+Diagnostic equation_error(std::string code, std::string message) {
+    return {
+        DiagnosticStage::equation,
+        Severity::error,
+        std::move(code),
+        std::move(message),
+        std::nullopt,
+    };
+}
+
+std::int32_t csr_index(const std::size_t value) {
+    if (value >
+        static_cast<std::size_t>(
+            std::numeric_limits<std::int32_t>::max())) {
+        throw std::overflow_error{
+            "Reduced system exceeds the 32-bit CSR index range."};
+    }
+    return static_cast<std::int32_t>(value);
+}
+
+}  // namespace
+
+ConstraintResult eliminate_essential_bcs(
+    const EquationSystem& original,
+    const DofManager& dofs) {
+    ConstraintResult result;
+    if (const auto error = validate_system(original);
+        error.has_value()) {
+        result.diagnostics.push_back(equation_error(
+            "equation.invalid_system", *error));
+        return result;
+    }
+    if (original.stiffness.order != dofs.full_dof_count()) {
+        result.diagnostics.push_back(equation_error(
+            "equation.dof_count_mismatch",
+            "Matrix order must equal the DofManager full DOF count."));
+        return result;
+    }
+
+    const std::size_t full_count = original.stiffness.order;
+    std::vector<std::size_t> full_to_free(full_count, full_count);
+    std::vector<double> prescribed_full(full_count, 0.0);
+    for (std::size_t full = 0; full < full_count; ++full) {
+        const std::optional<std::size_t> equation =
+            dofs.equation(full);
+        if (equation.has_value()) {
+            full_to_free[full] = *equation;
+        } else {
+            prescribed_full[full] =
+                dofs.prescribed_value(full).value();
+        }
+    }
+
+    std::vector<double> reduced_force(
+        dofs.free_equation_count(), 0.0);
+    for (std::size_t full = 0; full < full_count; ++full) {
+        if (full_to_free[full] != full_count) {
+            reduced_force[full_to_free[full]] = original.force[full];
+        }
+    }
+
+    SymmetricCsr reduced{
+        dofs.free_equation_count(),
+        std::vector<std::int32_t>(
+            dofs.free_equation_count() + 1, 0),
+        {},
+        {},
+    };
+    reduced.column_indices.reserve(
+        original.stiffness.column_indices.size());
+    reduced.values.reserve(original.stiffness.values.size());
+
+    for (std::size_t full_row = 0; full_row < full_count; ++full_row) {
+        const std::size_t begin = static_cast<std::size_t>(
+            original.stiffness.row_offsets[full_row]);
+        const std::size_t end = static_cast<std::size_t>(
+            original.stiffness.row_offsets[full_row + 1]);
+        for (std::size_t entry = begin; entry < end; ++entry) {
+            const std::size_t full_column =
+                static_cast<std::size_t>(
+                    original.stiffness.column_indices[entry]);
+            const double stiffness = original.stiffness.values[entry];
+
+            if (full_to_free[full_row] == full_count) {
+                if (full_to_free[full_column] != full_count) {
+                    reduced_force[full_to_free[full_column]] -=
+                        stiffness * prescribed_full[full_row];
+                }
+            } else if (full_to_free[full_column] == full_count) {
+                reduced_force[full_to_free[full_row]] -=
+                    stiffness * prescribed_full[full_column];
+            } else {
+                reduced.column_indices.push_back(
+                    csr_index(full_to_free[full_column]));
+                reduced.values.push_back(stiffness);
+            }
+        }
+        if (full_to_free[full_row] != full_count) {
+            reduced.row_offsets[full_to_free[full_row] + 1] =
+                csr_index(reduced.column_indices.size());
+        }
+    }
+    if (!std::ranges::all_of(
+            reduced_force, [](const double value) {
+                return std::isfinite(value);
+            })) {
+        result.diagnostics.push_back(equation_error(
+            "equation.nonfinite_result",
+            "Essential BC elimination produced a nonfinite force."));
+        return result;
+    }
+
+    result.reduced_system = ReducedSystem{
+        std::move(reduced),
+        std::move(reduced_force),
+    };
+    return result;
+}
+
+std::vector<double> recover_reaction(
+    const EquationSystem& original,
+    const std::span<const double> full_displacement) {
+    if (const auto error = validate_system(original);
+        error.has_value()) {
+        throw std::invalid_argument{*error};
+    }
+    if (full_displacement.size() != original.stiffness.order) {
+        throw std::invalid_argument{
+            "Full displacement size must equal the matrix order."};
+    }
+    if (!std::ranges::all_of(
+            full_displacement, [](const double value) {
+                return std::isfinite(value);
+            })) {
+        throw std::invalid_argument{
+            "Full displacement values must be finite."};
+    }
+
+    std::vector<double> reaction = original.force;
+    for (double& value : reaction) {
+        value = -value;
+    }
+    for (std::size_t row = 0;
+         row < original.stiffness.order;
+         ++row) {
+        const std::size_t begin = static_cast<std::size_t>(
+            original.stiffness.row_offsets[row]);
+        const std::size_t end = static_cast<std::size_t>(
+            original.stiffness.row_offsets[row + 1]);
+        for (std::size_t entry = begin; entry < end; ++entry) {
+            const std::size_t column =
+                static_cast<std::size_t>(
+                    original.stiffness.column_indices[entry]);
+            const double stiffness = original.stiffness.values[entry];
+            reaction[row] +=
+                stiffness * full_displacement[column];
+            if (column != row) {
+                reaction[column] +=
+                    stiffness * full_displacement[row];
+            }
+        }
+    }
+    if (!std::ranges::all_of(
+            reaction, [](const double value) {
+                return std::isfinite(value);
+            })) {
+        throw std::invalid_argument{
+            "Reaction recovery produced a nonfinite value."};
+    }
+    return reaction;
+}
+
+}  // namespace fesa

src/fesa/fem/dof_manager.cpp

@@ -60,7 +60,25 @@ std::size_t DofManager::free_equation_count() const noexcept {
 
 std::optional<std::size_t> DofManager::equation(
     const DofAddress address) const {
-    return equations_[full_dof(address)];
+    return equation(full_dof(address));
+}
+
+std::optional<std::size_t> DofManager::equation(
+    const std::size_t full_dof) const {
+    return equations_.at(full_dof);
+}
+
+std::optional<double> DofManager::prescribed_value(
+    const DofAddress address) const {
+    return prescribed_value(full_dof(address));
+}
+
+std::optional<double> DofManager::prescribed_value(
+    const std::size_t full_dof) const {
+    if (equations_.at(full_dof).has_value()) {
+        return std::nullopt;
+    }
+    return prescribed_values_.at(full_dof);
 }
 
 std::array<std::size_t, 12> DofManager::element_full_dofs(

src/fesa/model/domain.cpp

@@ -55,6 +55,24 @@ const Node& Domain::node(const EntityOrigin& origin) const {
     return nodes_[found->second];
 }
 
+const IsotropicElastic& Domain::material(const MaterialId id) const {
+    const auto found = material_indices_.find(id.value());
+    if (found == material_indices_.end()) {
+        throw std::out_of_range{
+            "Material ID is not present in the Domain."};
+    }
+    return materials_[found->second];
+}
+
+const BeamSection& Domain::section(const SectionId id) const {
+    const auto found = section_indices_.find(id.value());
+    if (found == section_indices_.end()) {
+        throw std::out_of_range{
+            "Section ID is not present in the Domain."};
+    }
+    return sections_[found->second];
+}
+
 Domain::OriginKey Domain::origin_key(const EntityOrigin& origin) {
     return {origin.instance_name, origin.local_label};
 }

src/fesa/solvers/linear/pardiso_linear_solver.cpp

@@ -0,0 +1,330 @@
+#include <fesa/solvers/linear/pardiso_linear_solver.hpp>
+
+#include <algorithm>
+#include <array>
+#include <cmath>
+#include <cstddef>
+#include <cstdint>
+#include <limits>
+#include <optional>
+#include <string>
+#include <string_view>
+#include <type_traits>
+#include <utility>
+#include <vector>
+
+#include <mkl.h>
+
+namespace fesa {
+namespace {
+
+static_assert(
+    std::is_same_v<MKL_INT, std::int32_t>,
+    "FESA requires the oneMKL LP64 interface.");
+
+Diagnostic solver_error(std::string code, std::string message) {
+    return {
+        DiagnosticStage::solver,
+        Severity::error,
+        std::move(code),
+        std::move(message),
+        std::nullopt,
+    };
+}
+
+std::optional<std::string> validate_matrix(
+    const SymmetricCsr& matrix) {
+    if (matrix.order == 0) {
+        return "PARDISO requires a nonempty reduced system.";
+    }
+    if (matrix.order >
+        static_cast<std::size_t>(
+            std::numeric_limits<MKL_INT>::max())) {
+        return "Matrix order exceeds the oneMKL LP64 index range.";
+    }
+    if (matrix.row_offsets.size() != matrix.order + 1 ||
+        matrix.row_offsets.front() != 0) {
+        return "CSR row offsets must contain order + 1 entries "
+               "starting at zero.";
+    }
+    if (matrix.column_indices.size() != matrix.values.size()) {
+        return "CSR column and value counts must match.";
+    }
+
+    MKL_INT previous_offset = 0;
+    for (const MKL_INT offset : matrix.row_offsets) {
+        if (offset < previous_offset || offset < 0 ||
+            static_cast<std::size_t>(offset) >
+                matrix.column_indices.size()) {
+            return "CSR row offsets must be nondecreasing and in range.";
+        }
+        previous_offset = offset;
+    }
+    if (static_cast<std::size_t>(matrix.row_offsets.back()) !=
+        matrix.column_indices.size()) {
+        return "The final CSR row offset must equal the entry count.";
+    }
+
+    for (std::size_t row = 0; row < matrix.order; ++row) {
+        MKL_INT previous_column = -1;
+        bool has_diagonal = false;
+        const std::size_t begin =
+            static_cast<std::size_t>(matrix.row_offsets[row]);
+        const std::size_t end =
+            static_cast<std::size_t>(matrix.row_offsets[row + 1]);
+        for (std::size_t entry = begin; entry < end; ++entry) {
+            const MKL_INT column = matrix.column_indices[entry];
+            if (column < static_cast<MKL_INT>(row) ||
+                column >= static_cast<MKL_INT>(matrix.order) ||
+                column <= previous_column) {
+                return "CSR rows must contain sorted unique "
+                       "upper-triangle columns.";
+            }
+            if (!std::isfinite(matrix.values[entry])) {
+                return "CSR values must be finite.";
+            }
+            has_diagonal =
+                has_diagonal || column == static_cast<MKL_INT>(row);
+            previous_column = column;
+        }
+        if (!has_diagonal) {
+            return "Every CSR row must contain its diagonal entry.";
+        }
+    }
+    return std::nullopt;
+}
+
+class PardisoSession final {
+public:
+    PardisoSession() {
+        pardisoinit(handles_.data(), &matrix_type_, parameters_.data());
+        parameters_[26] = 1;
+        parameters_[34] = 1;
+    }
+
+    PardisoSession(const PardisoSession&) = delete;
+    PardisoSession& operator=(const PardisoSession&) = delete;
+
+    ~PardisoSession() noexcept {
+        static_cast<void>(release());
+    }
+
+    MKL_INT release() noexcept {
+        if (!active_) {
+            return 0;
+        }
+        active_ = false;
+
+        constexpr MKL_INT release_all = -1;
+        MKL_INT error = 0;
+        pardiso(
+            handles_.data(),
+            &max_factorizations_,
+            &matrix_number_,
+            &matrix_type_,
+            &release_all,
+            &order_,
+            matrix_->values.data(),
+            matrix_->row_offsets.data(),
+            matrix_->column_indices.data(),
+            permutation_.data(),
+            &right_hand_side_count_,
+            parameters_.data(),
+            &message_level_,
+            right_hand_side_->data(),
+            solution_->data(),
+            &error);
+        return error;
+    }
+
+    MKL_INT execute(
+        const MKL_INT phase,
+        const SymmetricCsr& matrix,
+        std::vector<double>& right_hand_side,
+        std::vector<double>& solution) {
+        order_ = static_cast<MKL_INT>(matrix.order);
+        matrix_ = &matrix;
+        right_hand_side_ = &right_hand_side;
+        solution_ = &solution;
+        permutation_.resize(matrix.order);
+        active_ = true;
+
+        MKL_INT error = 0;
+        pardiso(
+            handles_.data(),
+            &max_factorizations_,
+            &matrix_number_,
+            &matrix_type_,
+            &phase,
+            &order_,
+            matrix.values.data(),
+            matrix.row_offsets.data(),
+            matrix.column_indices.data(),
+            permutation_.data(),
+            &right_hand_side_count_,
+            parameters_.data(),
+            &message_level_,
+            right_hand_side.data(),
+            solution.data(),
+            &error);
+        return error;
+    }
+
+private:
+    std::array<void*, 64> handles_{};
+    std::array<MKL_INT, 64> parameters_{};
+    std::vector<MKL_INT> permutation_;
+    const SymmetricCsr* matrix_{};
+    std::vector<double>* right_hand_side_{};
+    std::vector<double>* solution_{};
+    MKL_INT order_{};
+    MKL_INT max_factorizations_{1};
+    MKL_INT matrix_number_{1};
+    MKL_INT matrix_type_{2};
+    MKL_INT right_hand_side_count_{1};
+    MKL_INT message_level_{};
+    bool active_{};
+};
+
+double relative_residual(
+    const SymmetricCsr& matrix,
+    const std::span<const double> rhs,
+    const std::span<const double> solution) {
+    std::vector<double> residual(rhs.begin(), rhs.end());
+    for (double& value : residual) {
+        value = -value;
+    }
+
+    for (std::size_t row = 0; row < matrix.order; ++row) {
+        const std::size_t begin =
+            static_cast<std::size_t>(matrix.row_offsets[row]);
+        const std::size_t end =
+            static_cast<std::size_t>(matrix.row_offsets[row + 1]);
+        for (std::size_t entry = begin; entry < end; ++entry) {
+            const std::size_t column =
+                static_cast<std::size_t>(matrix.column_indices[entry]);
+            const double value = matrix.values[entry];
+            residual[row] += value * solution[column];
+            if (column != row) {
+                residual[column] += value * solution[row];
+            }
+        }
+    }
+
+    double residual_squared = 0.0;
+    double rhs_squared = 0.0;
+    for (std::size_t index = 0; index < rhs.size(); ++index) {
+        residual_squared += residual[index] * residual[index];
+        rhs_squared += rhs[index] * rhs[index];
+    }
+    const double residual_norm = std::sqrt(residual_squared);
+    const double rhs_norm = std::sqrt(rhs_squared);
+    return rhs_norm == 0.0 ? residual_norm : residual_norm / rhs_norm;
+}
+
+std::string pardiso_failure(
+    const std::string_view phase,
+    const MKL_INT error) {
+    return "PARDISO " + std::string{phase} +
+           " failed with error " + std::to_string(error) + ".";
+}
+
+}  // namespace
+
+PardisoLinearSolver::PardisoLinearSolver() = default;
+
+PardisoLinearSolver::~PardisoLinearSolver() = default;
+
+LinearSolveResult PardisoLinearSolver::solve(
+    const SymmetricCsr& matrix,
+    const std::span<const double> rhs) {
+    LinearSolveResult result;
+    if (const auto error = validate_matrix(matrix);
+        error.has_value()) {
+        result.diagnostics.push_back(
+            solver_error("solver.invalid_csr", *error));
+        return result;
+    }
+    if (rhs.size() != matrix.order) {
+        result.diagnostics.push_back(solver_error(
+            "solver.dimension_mismatch",
+            "Right-hand side size must equal the matrix order."));
+        return result;
+    }
+    if (!std::ranges::all_of(rhs, [](const double value) {
+            return std::isfinite(value);
+        })) {
+        result.diagnostics.push_back(solver_error(
+            "solver.invalid_rhs",
+            "Right-hand side values must be finite."));
+        return result;
+    }
+
+    std::vector<double> right_hand_side(rhs.begin(), rhs.end());
+    std::vector<double> solution(matrix.order, 0.0);
+    {
+        PardisoSession session;
+
+        constexpr MKL_INT analyze = 11;
+        MKL_INT error = session.execute(
+            analyze, matrix, right_hand_side, solution);
+        if (error != 0) {
+            result.diagnostics.push_back(solver_error(
+                "solver.analysis_failed",
+                pardiso_failure("analysis", error)));
+        }
+
+        if (error == 0) {
+            constexpr MKL_INT factorize = 22;
+            error = session.execute(
+                factorize, matrix, right_hand_side, solution);
+            if (error != 0) {
+                result.diagnostics.push_back(solver_error(
+                    "solver.factorization_failed",
+                    pardiso_failure("factorization", error)));
+            }
+        }
+
+        if (error == 0) {
+            constexpr MKL_INT solve_system = 33;
+            error = session.execute(
+                solve_system, matrix, right_hand_side, solution);
+            if (error != 0) {
+                result.diagnostics.push_back(solver_error(
+                    "solver.solve_failed",
+                    pardiso_failure("solve", error)));
+            }
+        }
+
+        if (const MKL_INT release_error = session.release();
+            release_error != 0) {
+            result.diagnostics.push_back(solver_error(
+                "solver.release_failed",
+                pardiso_failure("release", release_error)));
+        }
+    }
+    if (!result.diagnostics.empty()) {
+        return result;
+    }
+    if (!std::ranges::all_of(solution, [](const double value) {
+            return std::isfinite(value);
+        })) {
+        result.diagnostics.push_back(solver_error(
+            "solver.nonfinite_solution",
+            "PARDISO produced a nonfinite solution."));
+        return result;
+    }
+
+    result.relative_residual =
+        relative_residual(matrix, rhs, solution);
+    if (!std::isfinite(result.relative_residual)) {
+        result.diagnostics.push_back(solver_error(
+            "solver.nonfinite_residual",
+            "The independently computed relative residual is nonfinite."));
+        return result;
+    }
+    result.solution = std::move(solution);
+    return result;
+}
+
+}  // namespace fesa

tests/CMakeLists.txt

@@ -271,3 +271,107 @@ add_test(
     COMMAND "$<TARGET_FILE:fesa_beam3d2_tests>"
         --gtest_filter=RigidBody.*
 )
+
+add_executable(fesa_serial_assembly_tests
+    unit/assembly/serial_assembler_test.cpp
+)
+
+target_compile_features(fesa_serial_assembly_tests PRIVATE cxx_std_20)
+target_compile_options(
+    fesa_serial_assembly_tests
+    PRIVATE
+        /W4
+        /permissive-
+        /EHsc
+)
+
+target_link_libraries(fesa_serial_assembly_tests
+    PRIVATE
+        fesa_core
+        GTest::gtest_main
+)
+
+add_test(
+    NAME SparsePattern
+    COMMAND "$<TARGET_FILE:fesa_serial_assembly_tests>"
+        --gtest_filter=SparsePattern.*
+)
+
+add_test(
+    NAME SerialAssembly
+    COMMAND "$<TARGET_FILE:fesa_serial_assembly_tests>"
+        --gtest_filter=SerialAssembly.*
+)
+
+add_test(
+    NAME SymmetricCsr
+    COMMAND "$<TARGET_FILE:fesa_serial_assembly_tests>"
+        --gtest_filter=SymmetricCsr.*
+)
+
+add_executable(fesa_constraint_tests
+    unit/constraints/essential_bc_test.cpp
+)
+
+target_compile_features(fesa_constraint_tests PRIVATE cxx_std_20)
+target_compile_options(
+    fesa_constraint_tests
+    PRIVATE
+        /W4
+        /permissive-
+        /EHsc
+)
+
+target_link_libraries(fesa_constraint_tests
+    PRIVATE
+        fesa_core
+        GTest::gtest_main
+)
+
+add_test(
+    NAME EssentialBc
+    COMMAND "$<TARGET_FILE:fesa_constraint_tests>"
+        --gtest_filter=EssentialBc.*
+)
+
+add_test(
+    NAME ConstraintElimination
+    COMMAND "$<TARGET_FILE:fesa_constraint_tests>"
+        --gtest_filter=ConstraintElimination.*
+)
+
+add_test(
+    NAME Reaction
+    COMMAND "$<TARGET_FILE:fesa_constraint_tests>"
+        --gtest_filter=Reaction.*
+)
+
+add_executable(fesa_linear_solver_tests
+    unit/solvers/linear/pardiso_linear_solver_test.cpp
+)
+
+target_compile_features(fesa_linear_solver_tests PRIVATE cxx_std_20)
+target_compile_options(
+    fesa_linear_solver_tests
+    PRIVATE
+        /W4
+        /permissive-
+        /EHsc
+)
+
+target_link_libraries(fesa_linear_solver_tests
+    PRIVATE
+        fesa_core
+        GTest::gtest_main
+)
+
+add_test(
+    NAME PardisoLinearSolver
+    COMMAND "$<TARGET_FILE:fesa_linear_solver_tests>"
+)
+
+set_property(
+    TEST PardisoLinearSolver
+    PROPERTY ENVIRONMENT_MODIFICATION
+        ${FESA_DEPENDENCY_RUNTIME_MODIFICATIONS}
+)

tests/unit/assembly/serial_assembler_test.cpp

@@ -0,0 +1,494 @@
+#include <fesa/assembly/serial_assembler.hpp>
+
+#include <algorithm>
+#include <array>
+#include <bit>
+#include <cstddef>
+#include <cstdint>
+#include <limits>
+#include <stdexcept>
+#include <string>
+#include <utility>
+#include <vector>
+
+#include <fesa/model/domain_builder.hpp>
+
+#include <gtest/gtest.h>
+
+namespace {
+
+fesa::IsotropicElastic unit_material() {
+    return {
+        fesa::MaterialId{0},
+        "Unit",
+        1.0,
+        0.25,
+    };
+}
+
+fesa::BeamSection unit_section() {
+    return {
+        fesa::SectionId{0},
+        "Unit",
+        1.0,
+        1.0,
+        1.0,
+        1.0,
+        1.0,
+        1.0,
+        fesa::ShearPropertySource::input,
+        fesa::Vec3{0.0, 1.0, 0.0},
+        {},
+    };
+}
+
+fesa::Domain finish_domain(
+    fesa::DomainBuilder builder,
+    fesa::StepDefinition step) {
+    builder.set_step(std::move(step));
+    auto result = std::move(builder).build();
+    if (!result.domain.has_value()) {
+        throw std::runtime_error{"Test Domain failed validation."};
+    }
+    return std::move(*result.domain);
+}
+
+fesa::Domain build_chain_domain(const bool reverse_storage_order) {
+    fesa::DomainBuilder builder;
+    std::array<fesa::Node, 3> nodes{{
+        {
+            fesa::NodeId{20},
+            fesa::EntityOrigin{"BeamPart", "Beam-1", 2},
+            fesa::Vec3{2.0, 0.0, 0.0},
+        },
+        {
+            fesa::NodeId{4},
+            fesa::EntityOrigin{"BeamPart", "Beam-1", 100},
+            fesa::Vec3{0.0, 0.0, 0.0},
+        },
+        {
+            fesa::NodeId{10},
+            fesa::EntityOrigin{"BeamPart", "Beam-1", 50},
+            fesa::Vec3{1.0, 0.0, 0.0},
+        },
+    }};
+    if (reverse_storage_order) {
+        std::ranges::reverse(nodes);
+    }
+    for (fesa::Node& node : nodes) {
+        builder.add_node(std::move(node));
+    }
+
+    builder.add_material(unit_material());
+    builder.add_section(unit_section());
+
+    std::array<fesa::BeamElement, 2> elements{{
+        {
+            fesa::ElementId{1},
+            fesa::EntityOrigin{"BeamPart", "Beam-1", 200},
+            {fesa::NodeId{10}, fesa::NodeId{20}},
+            fesa::MaterialId{0},
+            fesa::SectionId{0},
+        },
+        {
+            fesa::ElementId{0},
+            fesa::EntityOrigin{"BeamPart", "Beam-1", 10},
+            {fesa::NodeId{4}, fesa::NodeId{10}},
+            fesa::MaterialId{0},
+            fesa::SectionId{0},
+        },
+    }};
+    if (reverse_storage_order) {
+        std::ranges::reverse(elements);
+    }
+    for (fesa::BeamElement& element : elements) {
+        builder.add_beam_element(std::move(element));
+    }
+
+    return finish_domain(
+        std::move(builder),
+        {
+            "Load",
+            {},
+            {
+                {
+                    fesa::NodeId{20},
+                    {3.0, 4.0, 5.0, 6.0, 7.0, 8.0},
+                },
+                {
+                    fesa::NodeId{10},
+                    {1.0, 2.0, 3.0, 4.0, 5.0, 6.0},
+                },
+                {
+                    fesa::NodeId{10},
+                    {-1.0, 10.0, -3.0, 0.0, 0.0, 0.0},
+                },
+            },
+        });
+}
+
+fesa::Domain build_parallel_domain(
+    const std::vector<std::int64_t>& element_labels) {
+    fesa::DomainBuilder builder;
+    builder.add_node({
+        fesa::NodeId{4},
+        fesa::EntityOrigin{"BeamPart", "Beam-1", 1},
+        fesa::Vec3{0.0, 0.0, 0.0},
+    });
+    builder.add_node({
+        fesa::NodeId{10},
+        fesa::EntityOrigin{"BeamPart", "Beam-1", 2},
+        fesa::Vec3{1.0, 0.0, 0.0},
+    });
+    builder.add_material(unit_material());
+    builder.add_section(unit_section());
+
+    for (std::size_t index = 0; index < element_labels.size(); ++index) {
+        builder.add_beam_element({
+            fesa::ElementId{static_cast<std::int64_t>(index)},
+            fesa::EntityOrigin{
+                "BeamPart",
+                "Beam-1",
+                element_labels[index],
+            },
+            {fesa::NodeId{4}, fesa::NodeId{10}},
+            fesa::MaterialId{0},
+            fesa::SectionId{0},
+        });
+    }
+    return finish_domain(std::move(builder), {"Load", {}, {}});
+}
+
+fesa::Domain build_rounding_domain(const bool large_element_first) {
+    fesa::DomainBuilder builder;
+    builder.add_node({
+        fesa::NodeId{4},
+        fesa::EntityOrigin{"BeamPart", "Beam-1", 1},
+        fesa::Vec3{0.0, 0.0, 0.0},
+    });
+    builder.add_node({
+        fesa::NodeId{10},
+        fesa::EntityOrigin{"BeamPart", "Beam-1", 2},
+        fesa::Vec3{1.0, 0.0, 0.0},
+    });
+
+    for (std::size_t index = 0; index < 3; ++index) {
+        fesa::IsotropicElastic material = unit_material();
+        material.id =
+            fesa::MaterialId{static_cast<std::int64_t>(index)};
+        material.name = "Material-" + std::to_string(index);
+        material.young = index == 2 ? 1.0e16 : 1.0;
+        builder.add_material(std::move(material));
+    }
+    builder.add_section(unit_section());
+
+    std::array<std::size_t, 3> storage_order{0, 1, 2};
+    if (large_element_first) {
+        storage_order = {2, 0, 1};
+    }
+    for (const std::size_t index : storage_order) {
+        builder.add_beam_element({
+            fesa::ElementId{static_cast<std::int64_t>(index)},
+            fesa::EntityOrigin{
+                "BeamPart",
+                "Beam-1",
+                static_cast<std::int64_t>((index + 1) * 10),
+            },
+            {fesa::NodeId{4}, fesa::NodeId{10}},
+            fesa::MaterialId{static_cast<std::int64_t>(index)},
+            fesa::SectionId{0},
+        });
+    }
+    return finish_domain(std::move(builder), {"Load", {}, {}});
+}
+
+fesa::Domain build_orphan_node_domain() {
+    fesa::DomainBuilder builder;
+    builder.add_node({
+        fesa::NodeId{0},
+        fesa::EntityOrigin{"BeamPart", "Beam-1", 1},
+        fesa::Vec3{0.0, 0.0, 0.0},
+    });
+    builder.add_node({
+        fesa::NodeId{1},
+        fesa::EntityOrigin{"BeamPart", "Beam-1", 2},
+        fesa::Vec3{1.0, 0.0, 0.0},
+    });
+    builder.add_node({
+        fesa::NodeId{2},
+        fesa::EntityOrigin{"BeamPart", "Beam-1", 3},
+        fesa::Vec3{2.0, 0.0, 0.0},
+    });
+    builder.add_material(unit_material());
+    builder.add_section(unit_section());
+    builder.add_beam_element({
+        fesa::ElementId{0},
+        fesa::EntityOrigin{"BeamPart", "Beam-1", 1},
+        {fesa::NodeId{0}, fesa::NodeId{1}},
+        fesa::MaterialId{0},
+        fesa::SectionId{0},
+    });
+    return finish_domain(std::move(builder), {"Load", {}, {}});
+}
+
+double csr_value(
+    const fesa::SymmetricCsr& matrix,
+    std::size_t row,
+    std::size_t column) {
+    if (column < row) {
+        std::swap(row, column);
+    }
+    const auto begin = matrix.column_indices.begin() + matrix.row_offsets[row];
+    const auto end =
+        matrix.column_indices.begin() + matrix.row_offsets[row + 1];
+    const auto found = std::lower_bound(
+        begin, end, static_cast<std::int32_t>(column));
+    if (found == end || *found != static_cast<std::int32_t>(column)) {
+        throw std::out_of_range{"CSR entry is not present."};
+    }
+    return matrix.values[static_cast<std::size_t>(
+        std::distance(matrix.column_indices.begin(), found))];
+}
+
+TEST(SparsePattern, BuildsExpectedTwoElementChainStructure) {
+    const fesa::Domain domain = build_chain_domain(false);
+    const fesa::DofManager dofs = fesa::DofManager::build(domain);
+
+    const fesa::EquationSystem system =
+        fesa::assemble_serial(domain, dofs);
+
+    EXPECT_EQ(
+        system.stiffness.row_offsets,
+        (std::vector<std::int32_t>{
+            0,
+            12,
+            23,
+            33,
+            42,
+            50,
+            57,
+            69,
+            80,
+            90,
+            99,
+            107,
+            114,
+            120,
+            125,
+            129,
+            132,
+            134,
+            135,
+        }));
+
+    std::vector<std::int32_t> expected_columns;
+    for (std::int32_t row = 0; row < 18; ++row) {
+        const std::int32_t last_column = row < 6 ? 11 : 17;
+        for (std::int32_t column = row; column <= last_column; ++column) {
+            expected_columns.push_back(column);
+        }
+    }
+    EXPECT_EQ(system.stiffness.column_indices, expected_columns);
+}
+
+TEST(SymmetricCsr, StoresSortedUpperTriangleWithValidOffsets) {
+    const fesa::Domain domain = build_chain_domain(false);
+    const fesa::DofManager dofs = fesa::DofManager::build(domain);
+    const fesa::SymmetricCsr matrix =
+        fesa::assemble_serial(domain, dofs).stiffness;
+
+    ASSERT_EQ(matrix.order, 18);
+    ASSERT_EQ(matrix.row_offsets.size(), matrix.order + 1);
+    ASSERT_EQ(matrix.row_offsets.front(), 0);
+    ASSERT_EQ(
+        static_cast<std::size_t>(matrix.row_offsets.back()),
+        matrix.column_indices.size());
+    ASSERT_EQ(matrix.column_indices.size(), matrix.values.size());
+
+    for (std::size_t row = 0; row < matrix.order; ++row) {
+        const auto begin =
+            matrix.column_indices.begin() + matrix.row_offsets[row];
+        const auto end =
+            matrix.column_indices.begin() + matrix.row_offsets[row + 1];
+        EXPECT_TRUE(std::ranges::is_sorted(begin, end));
+        EXPECT_EQ(std::adjacent_find(begin, end), end);
+        for (auto entry = begin; entry != end; ++entry) {
+            EXPECT_GE(*entry, static_cast<std::int32_t>(row));
+            EXPECT_LT(*entry, static_cast<std::int32_t>(matrix.order));
+        }
+    }
+}
+
+TEST(SparsePattern, IncludesZeroDiagonalForEveryOrphanNodeDof) {
+    const fesa::Domain domain = build_orphan_node_domain();
+    const fesa::SymmetricCsr matrix =
+        fesa::assemble_serial(
+            domain, fesa::DofManager::build(domain))
+            .stiffness;
+
+    ASSERT_EQ(matrix.order, 18);
+    for (std::size_t row = 0; row < matrix.order; ++row) {
+        const auto begin =
+            matrix.column_indices.begin() + matrix.row_offsets[row];
+        const auto end =
+            matrix.column_indices.begin() + matrix.row_offsets[row + 1];
+        const auto diagonal = std::lower_bound(
+            begin, end, static_cast<std::int32_t>(row));
+        ASSERT_NE(diagonal, end);
+        EXPECT_EQ(*diagonal, static_cast<std::int32_t>(row));
+    }
+    for (std::size_t row = 12; row < 18; ++row) {
+        EXPECT_DOUBLE_EQ(csr_value(matrix, row, row), 0.0);
+    }
+}
+
+TEST(SerialAssembly, AssemblesHandCalculatedAxialChainAndFullLoad) {
+    const fesa::Domain domain = build_chain_domain(false);
+    const fesa::DofManager dofs = fesa::DofManager::build(domain);
+
+    const fesa::EquationSystem system =
+        fesa::assemble_serial(domain, dofs);
+
+    EXPECT_DOUBLE_EQ(csr_value(system.stiffness, 0, 0), 1.0);
+    EXPECT_DOUBLE_EQ(csr_value(system.stiffness, 0, 6), -1.0);
+    EXPECT_DOUBLE_EQ(csr_value(system.stiffness, 6, 6), 2.0);
+    EXPECT_DOUBLE_EQ(csr_value(system.stiffness, 6, 12), -1.0);
+    EXPECT_DOUBLE_EQ(csr_value(system.stiffness, 12, 12), 1.0);
+
+    EXPECT_EQ(
+        system.force,
+        (std::vector<double>{
+            0.0,
+            0.0,
+            0.0,
+            0.0,
+            0.0,
+            0.0,
+            0.0,
+            12.0,
+            0.0,
+            4.0,
+            5.0,
+            6.0,
+            3.0,
+            4.0,
+            5.0,
+            6.0,
+            7.0,
+            8.0,
+        }));
+}
+
+TEST(SerialAssembly, MergesDuplicateElementContributions) {
+    const fesa::Domain single_domain = build_parallel_domain({10});
+    const fesa::Domain duplicate_domain =
+        build_parallel_domain({30, 10, 20});
+    const fesa::EquationSystem single = fesa::assemble_serial(
+        single_domain, fesa::DofManager::build(single_domain));
+    const fesa::EquationSystem duplicate = fesa::assemble_serial(
+        duplicate_domain, fesa::DofManager::build(duplicate_domain));
+
+    EXPECT_EQ(
+        duplicate.stiffness.row_offsets,
+        single.stiffness.row_offsets);
+    EXPECT_EQ(
+        duplicate.stiffness.column_indices,
+        single.stiffness.column_indices);
+    ASSERT_EQ(duplicate.stiffness.values.size(), single.stiffness.values.size());
+    for (std::size_t index = 0; index < single.stiffness.values.size();
+         ++index) {
+        const double value = single.stiffness.values[index];
+        EXPECT_DOUBLE_EQ(
+            duplicate.stiffness.values[index], (value + value) + value);
+    }
+}
+
+TEST(SerialAssembly, ReducesByElementOriginNotDomainStorageOrder) {
+    const fesa::Domain origin_order_domain =
+        build_rounding_domain(false);
+    const fesa::Domain storage_order_domain =
+        build_rounding_domain(true);
+
+    const double origin_order_value = csr_value(
+        fesa::assemble_serial(
+            origin_order_domain,
+            fesa::DofManager::build(origin_order_domain))
+            .stiffness,
+        0,
+        0);
+    const double storage_order_value = csr_value(
+        fesa::assemble_serial(
+            storage_order_domain,
+            fesa::DofManager::build(storage_order_domain))
+            .stiffness,
+        0,
+        0);
+    const double expected = (1.0 + 1.0) + 1.0e16;
+
+    EXPECT_EQ(
+        std::bit_cast<std::uint64_t>(origin_order_value),
+        std::bit_cast<std::uint64_t>(expected));
+    EXPECT_EQ(
+        std::bit_cast<std::uint64_t>(storage_order_value),
+        std::bit_cast<std::uint64_t>(expected));
+}
+
+TEST(SerialAssembly, IsIndependentOfDomainStorageAndExternalLabelOrder) {
+    const fesa::Domain first_domain = build_chain_domain(false);
+    const fesa::Domain second_domain = build_chain_domain(true);
+
+    const fesa::EquationSystem first = fesa::assemble_serial(
+        first_domain, fesa::DofManager::build(first_domain));
+    const fesa::EquationSystem second = fesa::assemble_serial(
+        second_domain, fesa::DofManager::build(second_domain));
+
+    EXPECT_EQ(first.stiffness.order, second.stiffness.order);
+    EXPECT_EQ(
+        first.stiffness.row_offsets, second.stiffness.row_offsets);
+    EXPECT_EQ(
+        first.stiffness.column_indices,
+        second.stiffness.column_indices);
+    EXPECT_EQ(first.stiffness.values, second.stiffness.values);
+    EXPECT_EQ(first.force, second.force);
+}
+
+TEST(SerialAssembly, PropagatesBeamKernelFailure) {
+    fesa::DomainBuilder builder;
+    builder.add_node({
+        fesa::NodeId{0},
+        fesa::EntityOrigin{"BeamPart", "Beam-1", 1},
+        fesa::Vec3{0.0, 0.0, 0.0},
+    });
+    builder.add_node({
+        fesa::NodeId{1},
+        fesa::EntityOrigin{"BeamPart", "Beam-1", 2},
+        fesa::Vec3{1.0, 0.0, 0.0},
+    });
+    auto material = unit_material();
+    material.young = std::numeric_limits<double>::max();
+    builder.add_material(std::move(material));
+    auto section = unit_section();
+    section.area = std::numeric_limits<double>::max();
+    builder.add_section(std::move(section));
+    builder.add_beam_element({
+        fesa::ElementId{0},
+        fesa::EntityOrigin{"BeamPart", "Beam-1", 1},
+        {fesa::NodeId{0}, fesa::NodeId{1}},
+        fesa::MaterialId{0},
+        fesa::SectionId{0},
+    });
+    const fesa::Domain domain =
+        finish_domain(std::move(builder), {"Load", {}, {}});
+
+    try {
+        static_cast<void>(fesa::assemble_serial(
+            domain, fesa::DofManager::build(domain)));
+        FAIL() << "Expected a Beam kernel failure.";
+    } catch (const std::runtime_error& error) {
+        EXPECT_NE(
+            std::string{error.what()}.find("model.nonfinite_value"),
+            std::string::npos);
+    }
+}
+
+}  // namespace

tests/unit/constraints/essential_bc_test.cpp

@@ -0,0 +1,210 @@
+#include <fesa/constraints/essential_bc.hpp>
+
+#include <algorithm>
+#include <array>
+#include <cstddef>
+#include <cstdint>
+#include <limits>
+#include <stdexcept>
+#include <utility>
+#include <vector>
+
+#include <fesa/model/domain_builder.hpp>
+
+#include <gtest/gtest.h>
+
+namespace {
+
+fesa::DofManager build_dofs(
+    std::vector<fesa::PrescribedDof> prescribed) {
+    fesa::DomainBuilder builder;
+    builder.add_node({
+        fesa::NodeId{0},
+        fesa::EntityOrigin{"", "", 1},
+        fesa::Vec3{0.0, 0.0, 0.0},
+    });
+    builder.set_step({"Load", std::move(prescribed), {}});
+    auto result = std::move(builder).build();
+    if (!result.domain.has_value()) {
+        throw std::runtime_error{"Test Domain failed validation."};
+    }
+    return fesa::DofManager::build(*result.domain);
+}
+
+fesa::EquationSystem hand_system() {
+    return {
+        {
+            6,
+            {0, 3, 5, 6, 7, 8, 9},
+            {0, 1, 2, 1, 2, 2, 3, 4, 5},
+            {4.0, 1.0, 2.0, 3.0, 6.0, 5.0, 1.0, 1.0, 1.0},
+        },
+        {7.0, 8.0, 9.0, 0.0, 0.0, 0.0},
+    };
+}
+
+TEST(EssentialBc, ZeroPrescribedValueDoesNotShiftReducedForce) {
+    const std::vector<fesa::PrescribedDof> prescribed{
+        {fesa::NodeId{0}, 2, 0.0},
+    };
+    const fesa::DofManager dofs = build_dofs(prescribed);
+
+    const fesa::ConstraintResult result =
+        fesa::eliminate_essential_bcs(hand_system(), dofs);
+
+    ASSERT_TRUE(result.reduced_system.has_value());
+    EXPECT_TRUE(result.diagnostics.empty());
+    EXPECT_EQ(
+        result.reduced_system->force,
+        (std::vector<double>{7.0, 9.0, 0.0, 0.0, 0.0}));
+}
+
+TEST(
+    ConstraintElimination,
+    ShiftsNonzeroPrescribedValueAndPreservesOriginalSystem) {
+    const std::vector<fesa::PrescribedDof> prescribed{
+        {fesa::NodeId{0}, 2, 2.0},
+    };
+    const fesa::DofManager dofs = build_dofs(prescribed);
+    const fesa::EquationSystem original = hand_system();
+    const fesa::EquationSystem before = original;
+
+    const fesa::ConstraintResult result =
+        fesa::eliminate_essential_bcs(original, dofs);
+
+    ASSERT_TRUE(result.reduced_system.has_value());
+    EXPECT_TRUE(result.diagnostics.empty());
+    const fesa::ReducedSystem& reduced = *result.reduced_system;
+    EXPECT_EQ(reduced.stiffness.order, 5);
+    EXPECT_EQ(
+        reduced.stiffness.row_offsets,
+        (std::vector<std::int32_t>{0, 2, 3, 4, 5, 6}));
+    EXPECT_EQ(
+        reduced.stiffness.column_indices,
+        (std::vector<std::int32_t>{0, 1, 1, 2, 3, 4}));
+    EXPECT_EQ(
+        reduced.stiffness.values,
+        (std::vector<double>{4.0, 2.0, 5.0, 1.0, 1.0, 1.0}));
+    EXPECT_EQ(
+        reduced.force,
+        (std::vector<double>{5.0, -3.0, 0.0, 0.0, 0.0}));
+
+    EXPECT_EQ(original.stiffness.order, before.stiffness.order);
+    EXPECT_EQ(
+        original.stiffness.row_offsets,
+        before.stiffness.row_offsets);
+    EXPECT_EQ(
+        original.stiffness.column_indices,
+        before.stiffness.column_indices);
+    EXPECT_EQ(original.stiffness.values, before.stiffness.values);
+    EXPECT_EQ(original.force, before.force);
+
+    EXPECT_EQ(
+        dofs.reconstruct_full(
+            std::vector<double>{1.0, 3.0, 0.0, 0.0, 0.0}),
+        (std::vector<double>{1.0, 2.0, 3.0, 0.0, 0.0, 0.0}));
+}
+
+TEST(ConstraintElimination, AcceptsAllDofsConstrained) {
+    const std::vector<fesa::PrescribedDof> prescribed{
+        {fesa::NodeId{0}, 1, 1.0},
+        {fesa::NodeId{0}, 2, 2.0},
+        {fesa::NodeId{0}, 3, 3.0},
+        {fesa::NodeId{0}, 4, 4.0},
+        {fesa::NodeId{0}, 5, 5.0},
+        {fesa::NodeId{0}, 6, 6.0},
+    };
+    const fesa::DofManager dofs = build_dofs(prescribed);
+
+    const fesa::ConstraintResult result =
+        fesa::eliminate_essential_bcs(hand_system(), dofs);
+
+    ASSERT_TRUE(result.reduced_system.has_value());
+    const fesa::ReducedSystem& reduced = *result.reduced_system;
+    EXPECT_EQ(reduced.stiffness.order, 0);
+    EXPECT_EQ(
+        reduced.stiffness.row_offsets,
+        (std::vector<std::int32_t>{0}));
+    EXPECT_TRUE(reduced.stiffness.column_indices.empty());
+    EXPECT_TRUE(reduced.stiffness.values.empty());
+    EXPECT_TRUE(reduced.force.empty());
+    EXPECT_EQ(
+        dofs.reconstruct_full({}),
+        (std::vector<double>{1.0, 2.0, 3.0, 4.0, 5.0, 6.0}));
+}
+
+TEST(ConstraintElimination, RejectsNonfiniteSystemData) {
+    const fesa::DofManager dofs = build_dofs({});
+    fesa::EquationSystem nonfinite_matrix = hand_system();
+    nonfinite_matrix.stiffness.values[0] =
+        std::numeric_limits<double>::quiet_NaN();
+    fesa::EquationSystem nonfinite_force = hand_system();
+    nonfinite_force.force[0] =
+        std::numeric_limits<double>::infinity();
+
+    const fesa::ConstraintResult matrix_result =
+        fesa::eliminate_essential_bcs(nonfinite_matrix, dofs);
+    const fesa::ConstraintResult force_result =
+        fesa::eliminate_essential_bcs(nonfinite_force, dofs);
+
+    ASSERT_FALSE(matrix_result.diagnostics.empty());
+    EXPECT_EQ(
+        matrix_result.diagnostics.front().code,
+        "equation.invalid_system");
+    ASSERT_FALSE(force_result.diagnostics.empty());
+    EXPECT_EQ(
+        force_result.diagnostics.front().code,
+        "equation.invalid_system");
+}
+
+TEST(ConstraintElimination, RejectsNonfiniteReducedForce) {
+    const fesa::DofManager dofs = build_dofs({
+        {fesa::NodeId{0}, 2, 2.0},
+    });
+    fesa::EquationSystem overflowing = hand_system();
+    overflowing.stiffness.values[1] =
+        std::numeric_limits<double>::max();
+
+    const fesa::ConstraintResult result =
+        fesa::eliminate_essential_bcs(overflowing, dofs);
+
+    EXPECT_FALSE(result.reduced_system.has_value());
+    ASSERT_FALSE(result.diagnostics.empty());
+    EXPECT_EQ(
+        result.diagnostics.front().code,
+        "equation.nonfinite_result");
+}
+
+TEST(Reaction, UsesOriginalFullEquilibriumEquation) {
+    const std::vector<double> reaction = fesa::recover_reaction(
+        hand_system(),
+        std::vector<double>{1.0, 2.0, 3.0, 0.0, 0.0, 0.0});
+
+    EXPECT_EQ(
+        reaction,
+        (std::vector<double>{5.0, 17.0, 20.0, 0.0, 0.0, 0.0}));
+}
+
+TEST(Reaction, RejectsNonfiniteFullDisplacement) {
+    std::vector<double> displacement(6, 0.0);
+    displacement[0] = std::numeric_limits<double>::quiet_NaN();
+
+    EXPECT_THROW(
+        static_cast<void>(
+            fesa::recover_reaction(hand_system(), displacement)),
+        std::invalid_argument);
+}
+
+TEST(Reaction, RejectsNonfiniteRecoveredValue) {
+    fesa::EquationSystem overflowing = hand_system();
+    overflowing.stiffness.values[0] =
+        std::numeric_limits<double>::max();
+
+    EXPECT_THROW(
+        static_cast<void>(fesa::recover_reaction(
+            overflowing,
+            std::vector<double>{2.0, 0.0, 0.0, 0.0, 0.0, 0.0})),
+        std::invalid_argument);
+}
+
+}  // namespace

tests/unit/fem/dof_manager_test.cpp

@@ -130,6 +130,16 @@ TEST(DofManager, NumbersOnlyFreeDofsAndReconstructsPrescribedValues) {
     EXPECT_EQ(
         dofs.equation({fesa::NodeId{20}, fesa::NodeDof::ux}),
         std::optional<std::size_t>{10});
+    EXPECT_EQ(dofs.equation(std::size_t{0}), std::nullopt);
+    EXPECT_EQ(
+        dofs.prescribed_value(std::size_t{0}),
+        std::optional<double>{1.25});
+    EXPECT_EQ(
+        dofs.equation(std::size_t{1}),
+        std::optional<std::size_t>{0});
+    EXPECT_EQ(
+        dofs.prescribed_value(std::size_t{1}),
+        std::nullopt);
 
     std::vector<double> reduced(dofs.free_equation_count());
     for (std::size_t equation = 0; equation < reduced.size(); ++equation) {
@@ -190,6 +200,13 @@ TEST(DofManager, RejectsInvalidAddressesAndReducedVectorSize) {
         static_cast<void>(
             dofs.equation({fesa::NodeId{99}, fesa::NodeDof::ux})),
         std::out_of_range);
+    EXPECT_THROW(
+        static_cast<void>(dofs.equation(dofs.full_dof_count())),
+        std::out_of_range);
+    EXPECT_THROW(
+        static_cast<void>(
+            dofs.prescribed_value(dofs.full_dof_count())),
+        std::out_of_range);
     EXPECT_THROW(
         static_cast<void>(
             dofs.reconstruct_full(std::vector<double>(17, 0.0))),

tests/unit/solvers/linear/pardiso_linear_solver_test.cpp

@@ -0,0 +1,154 @@
+#include <fesa/solvers/linear/pardiso_linear_solver.hpp>
+
+#include <algorithm>
+#include <cmath>
+#include <cstddef>
+#include <cstdint>
+#include <span>
+#include <string_view>
+#include <type_traits>
+#include <vector>
+
+#include <gtest/gtest.h>
+
+namespace {
+
+static_assert(
+    !std::is_copy_constructible_v<fesa::PardisoLinearSolver>);
+static_assert(
+    !std::is_copy_assignable_v<fesa::PardisoLinearSolver>);
+
+fesa::SymmetricCsr spd_matrix() {
+    return {
+        3,
+        {0, 2, 4, 5},
+        {0, 1, 1, 2, 2},
+        {4.0, 1.0, 3.0, 1.0, 2.0},
+    };
+}
+
+double independent_relative_residual(
+    const fesa::SymmetricCsr& matrix,
+    const std::span<const double> rhs,
+    const std::span<const double> solution) {
+    std::vector<double> residual(rhs.begin(), rhs.end());
+    for (double& value : residual) {
+        value = -value;
+    }
+
+    for (std::size_t row = 0; row < matrix.order; ++row) {
+        const std::size_t begin =
+            static_cast<std::size_t>(matrix.row_offsets[row]);
+        const std::size_t end =
+            static_cast<std::size_t>(matrix.row_offsets[row + 1]);
+        for (std::size_t entry = begin; entry < end; ++entry) {
+            const std::size_t column =
+                static_cast<std::size_t>(matrix.column_indices[entry]);
+            const double value = matrix.values[entry];
+            residual[row] += value * solution[column];
+            if (column != row) {
+                residual[column] += value * solution[row];
+            }
+        }
+    }
+
+    double residual_squared = 0.0;
+    double rhs_squared = 0.0;
+    for (std::size_t index = 0; index < rhs.size(); ++index) {
+        residual_squared += residual[index] * residual[index];
+        rhs_squared += rhs[index] * rhs[index];
+    }
+    const double residual_norm = std::sqrt(residual_squared);
+    const double rhs_norm = std::sqrt(rhs_squared);
+    return rhs_norm == 0.0 ? residual_norm : residual_norm / rhs_norm;
+}
+
+bool has_diagnostic(
+    const fesa::LinearSolveResult& result,
+    const std::string_view code) {
+    return std::ranges::find(
+               result.diagnostics, code, &fesa::Diagnostic::code) !=
+           result.diagnostics.end();
+}
+
+TEST(PardisoLinearSolver, SolvesThreeByThreeSpdSystem) {
+    fesa::PardisoLinearSolver pardiso;
+    fesa::LinearSolver& solver = pardiso;
+    const fesa::SymmetricCsr matrix = spd_matrix();
+    const std::vector<double> rhs{6.0, 10.0, 8.0};
+
+    const fesa::LinearSolveResult result = solver.solve(matrix, rhs);
+
+    ASSERT_TRUE(result.diagnostics.empty());
+    ASSERT_EQ(result.solution.size(), 3);
+    EXPECT_NEAR(result.solution[0], 1.0, 1.0e-12);
+    EXPECT_NEAR(result.solution[1], 2.0, 1.0e-12);
+    EXPECT_NEAR(result.solution[2], 3.0, 1.0e-12);
+    const double independent =
+        independent_relative_residual(matrix, rhs, result.solution);
+    EXPECT_NEAR(result.relative_residual, independent, 1.0e-15);
+    EXPECT_LT(result.relative_residual, 1.0e-12);
+}
+
+TEST(PardisoLinearSolver, SupportsRepeatedSolve) {
+    fesa::PardisoLinearSolver solver;
+    const fesa::SymmetricCsr matrix = spd_matrix();
+
+    const fesa::LinearSolveResult first =
+        solver.solve(matrix, std::vector<double>{6.0, 10.0, 8.0});
+    const fesa::LinearSolveResult second =
+        solver.solve(matrix, std::vector<double>{-3.5, 2.5, 4.5});
+
+    ASSERT_TRUE(first.diagnostics.empty());
+    ASSERT_TRUE(second.diagnostics.empty());
+    ASSERT_EQ(second.solution.size(), 3);
+    EXPECT_NEAR(second.solution[0], -1.0, 1.0e-12);
+    EXPECT_NEAR(second.solution[1], 0.5, 1.0e-12);
+    EXPECT_NEAR(second.solution[2], 2.0, 1.0e-12);
+}
+
+TEST(LinearSolver, RejectsInvalidUpperTriangleCsr) {
+    fesa::PardisoLinearSolver solver;
+    fesa::SymmetricCsr invalid = spd_matrix();
+    invalid.column_indices[2] = 0;
+
+    const fesa::LinearSolveResult result =
+        solver.solve(invalid, std::vector<double>{6.0, 10.0, 8.0});
+
+    EXPECT_TRUE(result.solution.empty());
+    EXPECT_TRUE(has_diagnostic(result, "solver.invalid_csr"));
+}
+
+TEST(LinearSolver, RejectsRhsDimensionMismatch) {
+    fesa::PardisoLinearSolver solver;
+
+    const fesa::LinearSolveResult result =
+        solver.solve(spd_matrix(), std::vector<double>{1.0, 2.0});
+
+    EXPECT_TRUE(result.solution.empty());
+    EXPECT_TRUE(has_diagnostic(result, "solver.dimension_mismatch"));
+}
+
+TEST(PardisoLinearSolver, ReportsSingularMatrix) {
+    fesa::PardisoLinearSolver solver;
+    const fesa::SymmetricCsr singular{
+        2,
+        {0, 2, 3},
+        {0, 1, 1},
+        {1.0, 1.0, 1.0},
+    };
+
+    const fesa::LinearSolveResult result =
+        solver.solve(singular, std::vector<double>{2.0, 2.0});
+
+    EXPECT_TRUE(result.solution.empty());
+    EXPECT_FALSE(result.diagnostics.empty());
+    EXPECT_TRUE(std::ranges::all_of(
+        result.diagnostics,
+        [](const fesa::Diagnostic& diagnostic) {
+            return diagnostic.stage == fesa::DiagnosticStage::solver &&
+                   diagnostic.severity == fesa::Severity::error;
+        }));
+}
+
+}  // namespace