modify handoff.md

CMakeLists.txt

@@ -27,6 +27,11 @@ include(cmake/FesaDependencies.cmake)
 
 add_library(fesa_core STATIC
     src/fesa/core/version.cpp
+    src/fesa/elements/beam/beam3d2.cpp
+    src/fesa/fem/beam_frame.cpp
+    src/fesa/fem/dof_manager.cpp
+    src/fesa/fem/gauss_rule.cpp
+    src/fesa/fem/line2_shape.cpp
     src/fesa/io/abaqus/parser.cpp
     src/fesa/io/abaqus/semantic_mapper.cpp
     src/fesa/model/domain.cpp

docs/HANDOFF.md

@@ -2,8 +2,8 @@
 
 ## 1. 문서 목적
 
-이 문서는 `domain-and-input-skeleton` 완료 후 새 세션에서
-`fem-and-beam-kernel` Phase를 바로 시작하기 위한 인수인계 기록이다.
+이 문서는 `fem-and-beam-kernel` 완료 후 새 세션에서
+`equation-and-linear-solve` Phase를 바로 시작하기 위한 인수인계 기록이다.
 요구사항과 설계의 기준은 이 문서가 아니라 다음 파일이다.
 
 - `/AGENTS.md`
@@ -11,9 +11,10 @@
 - `/docs/ARCHITECTURE.md`
 - `/docs/ADR.md`
 - `/docs/HARNESS.md`
+- `/docs/formulation/timoshenko-beam-3d.md`
 - `/docs/superpowers/plans/2026-07-29-fesa-phase-1.md`
-- `/phases/fem-and-beam-kernel/index.json`
-- `/phases/fem-and-beam-kernel/step0.md`부터 `step3.md`
+- `/phases/equation-and-linear-solve/index.json`
+- `/phases/equation-and-linear-solve/step0.md`부터 `step2.md`
 
 내용이 충돌하면 위 기준 문서와 `phases/`의 현재 상태를 우선한다.
 
@@ -22,58 +23,81 @@
 2026-07-31 확인 기준:
 
 - 현재 브랜치: `dev`
-- `dev`, `origin/dev`, `origin/HEAD` 기준 커밋:
-  `e7f0d53406bfdff0fbd610c2b595a8d214341353`
+- 현재 `dev` HEAD:
+  `fa48a4d9d717224ed84ec41f06c43fbc33d42780`
+- `origin/dev`, `origin/HEAD`:
+  `4ee3895915accf69c3f988bf14119fab7c318ed5`
+- 로컬 `dev`는 `origin/dev`보다 10 commit 앞서고 0 commit 뒤처져 있다.
 - 완료 Phase:
   - `solver-bootstrap`
   - `domain-and-input-skeleton`
-- 다음 Phase: `fem-and-beam-kernel`
-- 다음 Step: `0 - quadrature-and-shape-functions`
-- `fem-and-beam-kernel`의 Step 0~3은 모두 `pending`
-- 활성 제품 코드 blocker는 없다.
-- `include/fesa/fem/`, `include/fesa/elements/beam/`,
-  `docs/formulation/timoshenko-beam-3d.md`는 아직 존재하지 않는다. 해당 Step에서
-  실제 구현과 함께 만든다.
+  - `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` 브랜치는 삭제되었다.
+- 원격 push는 수행하지 않았다.
 
-이 문서를 갱신하기 직전 작업 트리는 clean이었다. 문서 수정 자체가 아직 커밋되지 않은
-경우 Harness 실행 전에 `docs/HANDOFF.md` 변경을 먼저 커밋하거나 별도로 정리해야 한다.
+이 문서를 갱신하기 직전 작업 트리는 clean이었다. 이 문서 변경은 사용자가 별도로
+요청하지 않는 한 커밋하지 않는다. 새 세션에서 Harness를 실행하기 전에
+`docs/HANDOFF.md` 변경을 먼저 커밋하거나 별도로 정리해야 한다. 그렇지 않으면
+executor의 feature branch와 step commit에 인수인계 문서가 섞일 수 있다.
 
-## 3. 완료된 기반
+## 3. 완료된 `fem-and-beam-kernel`
 
-### `solver-bootstrap`
+Phase metadata는 `/phases/fem-and-beam-kernel/index.json`에 기록되어 있으며 Step
+0~3이 모두 `completed`다.
 
-- Visual Studio 2026 MSVC v145, Windows x64, C++20 CMake Preset
-- `fesa_core`와 `fesa` CLI 골격
-- MKL, TBB, HDF5, GoogleTest package discovery와 dependency smoke test
-- typed entity ID, `Vec3`, source location, diagnostic, status 값 타입
-- Harness Python characterization test와 CMake/CTest 실행 계약
+### FEM primitives
 
-### `domain-and-input-skeleton`
+- `/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 테스트
 
-- solver semantic 값 타입:
-  - `Node`, `BeamElement`, `IsotropicElastic`, `BeamSection`
-  - `NodeSet`, `ElementSet`, `StepDefinition`
-  - `EntityOrigin`과 typed internal ID
-- build 이후 const view만 제공하는 `Domain`
-- `DomainBuilder`의 중복 ID/origin, 참조, 유한값, 물성, 길이, orientation,
-  경계조건 검증
-- flat/orphan mesh와 Part/Assembly/단일 무변환 Instance scope를 보존하는 Abaqus
-  syntax parser
-- 활성 Part만 선택해 flat `Domain`으로 정규화하는 semantic mapper
-- 실제 Part/Instance/local label provenance 보존
-- 생략된 전단강성에
-  \(A_{sy}=A_{sz}=5A/6\)과
-  `ShearPropertySource::phase1_default` 적용
-- 여러 Instance, Instance transform, missing Part 등 최소 조직 오류의 source
-  diagnostic
-- public parser/mapper 경로를 사용하는 flat/계층형 통합 테스트
+### `DofManager`
 
-현재 parser/mapper는 다음 Phase의 Beam kernel 테스트에 필요한 최소 semantic
-`Domain`을 제공한다. 전체 Abaqus 부분집합, nested/`GENERATE` set, 완전한 scope
-resolution, explicit transverse shear, no-op directive와 단일 step의 모든 option은
-후속 `abaqus-subset-completion` Phase 책임이다. 이번 Phase에서 선행 구현하지 않는다.
+- `/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`에 저장하지 않는다.
 
-## 4. 검증된 개발환경과 baseline
+현재 `DofManager::build(const Domain&)`는 `Domain::step().prescribed_dofs`를 읽어
+구속 자유도와 prescribed 값을 이미 내부에 반영한다. 다음 Phase에서 이 상태를
+중복 소유하지 않도록 주의한다.
+
+### Beam frame와 transformation
+
+- `/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 테스트
+
+### Timoshenko stiffness kernel
+
+- `/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`에서 정정했다.
+
+## 4. 검증된 baseline과 개발환경
 
 새 PowerShell 세션에서 configure 또는 Harness 실행 전에 다음 환경 변수를 설정한다.
 절대경로를 tracked CMake 파일이나 Preset에 넣지 않는다.
@@ -90,11 +114,12 @@ Test-Path "$env:HDF5_DIR\hdf5-config.cmake"
 Test-Path "$env:GTest_DIR\GTestConfig.cmake"
 ```
 
-2026-07-31 위 네 경로가 모두 존재함을 확인했다. 같은 환경에서 다음 baseline을
+2026-07-31 현재 네 package 경로가 모두 존재한다. 현재 `dev` HEAD에서 다음 baseline을
 새로 검증했다.
 
-- MSVC v145 Debug build 성공, 새 경고 없음
-- CTest 12개 중 12개 성공
+- MSBuild 18.8.2, MSVC v145 Debug build 성공
+- build 출력에 새 warning 없음
+- CTest 22개 중 22개 성공
 - Harness pytest 20개 중 20개 성공
 - pytest가 실제로 20개를 수집했으므로 0-test 성공이 아님
 
@@ -115,99 +140,143 @@ cmake --fresh --preset windows-debug
 
 ## 5. 다음 Phase 목표와 Step 순서
 
-`fem-and-beam-kernel`의 독립 deliverable은 실제 2절점 3D Timoshenko Beam의
-local/global stiffness와 해석적 sanity test다. 네 Step을 순서대로 실행한다.
+`equation-and-linear-solve`의 독립 deliverable은 Beam contribution과 nodal load를
+full symmetric system으로 조립하고, essential BC를 소거한 뒤 MKL PARDISO로
+reduced SPD system을 푸는 것이다. 세 Step을 순서대로 실행한다.
 
-### Step 0 — `quadrature-and-shape-functions`
+### Step 0 — `symmetric-csr-assembly`
 
-- `fem`에 1점/2점 1D Gauss rule, 2절점 선형 shape function과 derivative,
-  `length/2` Jacobian을 구현한다.
-- partition of unity, endpoint interpolation, derivative sum zero, 적분 정확도,
-  invalid order/length를 실패 테스트로 먼저 고정한다.
-- runtime registry나 Beam stiffness를 만들지 않는다.
+- `/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 1 — `dof-manager`
+### Step 1 — `essential-bc-elimination`
 
-- 절점당 자유도 순서는
-  \(u_x,u_y,u_z,r_x,r_y,r_z\)의 6개다.
-- `DofManager`가 full DOF, constrained/free equation numbering,
-  12개 요소 DOF mapping과 prescribed value를 포함한 full-vector reconstruction을
-  전담한다.
-- external label이나 입력 선언 순서에 흔들리지 않는 deterministic numbering을
-  테스트로 고정한다.
-- equation ID를 `Node`나 `BeamElement`에 저장하지 않는다.
-- sparse pattern, MPC, penalty DOF는 이 Step 범위가 아니다.
-
-### Step 2 — `beam-local-frame`
-
-- 두 절점과 mandatory orientation vector로 오른손 직교
-  `BeamFrame {ex, ey, ez}`를 만든다.
-- Gram-Schmidt, 정규화, 직교성, determinant \(+1\), 회전/에너지 invariant를
+- `/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, 충돌 조건과 반력을 먼저
   테스트한다.
-- zero length, zero orientation, 요소축과 평행한 orientation은 diagnostic으로
-  실패시킨다.
-- `Matrix12`는 현재 코드에 없으므로 이 Step의 transformation 계약에 필요한 최소
-  고정 크기 타입만 만든다. 범용 동적 matrix 계층이나 registry를 만들지 않는다.
 
-### Step 3 — `timoshenko-stiffness-kernel`
+### Step 2 — `pardiso-linear-solver`
 
-- production code보다 먼저
-  `/docs/formulation/timoshenko-beam-3d.md`를 작성한다.
-- 문서에는 최소한 다음을 고정한다.
-  - 12개 local/global DOF와 component 순서
-  - 오른손 국부 좌표계 정의
-  - 변형률과 부호 규약
-  - \(G=E/[2(1+\nu)]\)와 constitutive matrix
-  - 자연좌표, shape derivative와 \(J=L/2\)
-  - 축·굽힘·비틀림 2점, 전단 1점 선택적 감차적분
-  - local/global transformation 규칙
-- 실제 `Beam3D2Input`, `Beam3D2Contribution`, `compute_beam3d2` kernel을 구현한다.
-- 대칭성, 강체운동 zero energy, 축/비틀림/단축·이축 굽힘 해석해,
-  shear-dominant 문제, 세장비 sweep과 회전 invariant를 실패 테스트로 먼저 만든다.
-- 임시 가짜/닫힌형 stiffness로 파이프라인만 통과시키지 않는다.
+- `/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 외부의 독립 계산으로 확인한다.
 
-## 6. 수치·정식화 계약
+## 6. 다음 Phase의 핵심 계약과 확인할 설계점
 
-- 요소: 2절점 직선 3D Isoparametric Timoshenko Beam
-- 절점당 6 DOF:
-  \(u_x,u_y,u_z,\theta_x,\theta_y,\theta_z\)
-- 선형 shape function, \(\xi\in[-1,1]\)
-- 축·굽힘·비틀림: 2점 Gauss 적분
-- 전단: 1점 Gauss 적분
-- 재료: 등방성 선형 탄성,
-  \(G=E/[2(1+\nu)]\)
-- 단면:
-  \(A,I_y,I_z,J,A_{sy},A_{sz}\)
-- 제한:
-  \(I_{yz}=0\), 도심/전단중심 일치, offset/warping 없음
-- orientation은 필수이며 자동 추측하거나 임의 축으로 대체하지 않는다.
-- tolerance는 테스트 이름이나 주석에 수학적/scale-aware 근거를 기록한다.
-- FESA는 단위 변환을 하지 않는다.
+### 아키텍처 경계
 
-정식화는 `docs/PRD.md` 3.2절, `docs/ARCHITECTURE.md` 8절,
-`docs/ADR.md` ADR-005와 ADR-014를 우선한다.
+- `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 범위가 아니다.
 
-## 7. 아키텍처 경계와 이번 Phase 제외 범위
+### 구현 전에 명시적으로 정렬할 점
 
-- public header는 `include/fesa/`, 구현은 `src/fesa/`, 테스트는 `tests/`에 둔다.
-- `core`, `model`, `fem`, `elements`는 Abaqus, MKL, TBB, HDF5 API에 의존하지
-  않는다.
-- `fem`은 DOF, equation mapping, quadrature, shape function, Jacobian,
-  local/global mapping만 담당한다.
-- `elements/beam`은 Beam local contribution과 kernel 계약을 담당한다.
-- `Domain`은 읽기 전용 view로 사용하고 복제하거나 mutable accessor를 추가하지 않는다.
-- 이번 Phase에서는 다음을 구현하지 않는다.
-  - sparse COO/CSR pattern과 조립
-  - essential BC elimination system
-  - MKL PARDISO solve
-  - oneTBB parallel assembly
-  - HDF5 writer와 result recovery
-  - 여러 요소 타입을 위한 registry/factory
-  - 비선형, warping, offset, \(I_{yz}\ne0\), 여러 Instance
+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에 넘기는 정책과 혼동하지 않는다.
+
+위 항목은 범위를 넓히라는 의미가 아니다. 기존 코드와 Step 초안 사이의 중복 또는
+미정 계약을 구현 전에 드러내고 가장 단순한 일관된 계약을 선택하기 위한 확인
+목록이다.
+
+## 7. Child sandbox의 MSBuild 실행 조건
+
+`fem-and-beam-kernel` 실행 전 child sandbox를 실제로 조사한 결과:
+
+- 앱 설치 경로의 `codex-cli 0.146.0`은 matching `codex-resources`를 찾지 못했다.
+- 완전한 standalone 배포는 다음 위치에 있다.
+  `C:\Users\baram\.codex\packages\standalone\releases\0.146.0-x86_64-pc-windows-msvc`
+- child에서 WindowsApps 경유 PowerShell은 실패했고 Windows PowerShell 5.1은
+  정상 동작했다.
+- `[windows] sandbox = "elevated"`에서는 MSBuild FileTracker가 access denied로
+  실패했다.
+- 완전한 standalone Codex, 정리된 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`로 실제 경로를 확인한다.
+
+```powershell
+$codexReleaseBin = "C:\Users\baram\.codex\packages\standalone\releases\0.146.0-x86_64-pc-windows-msvc\bin"
+$windowsPowerShell = "$env:SystemRoot\System32\WindowsPowerShell\v1.0"
+$filteredPath = $env:PATH -split ";" | Where-Object {
+    $_ -and
+    $_ -ne $codexReleaseBin -and
+    $_ -ne $windowsPowerShell -and
+    $_ -notmatch "\\Microsoft\\WindowsApps$" -and
+    $_ -notmatch "\\OpenAI\\Codex\\bin$"
+}
+$env:PATH = (@($codexReleaseBin, $windowsPowerShell) + $filteredPath) -join ";"
+
+(Get-Command codex).Source
+(Get-Command powershell).Source
+```
+
+4. 같은 세션에서 package 환경 변수와 baseline을 확인한 뒤 Harness를 실행한다.
+5. 성공·실패와 무관하게 global config를 즉시
+   `[windows] sandbox = "elevated"`로 복원하고 실제 값을 다시 읽어 확인한다.
+
+사용자 profile 전체나 드라이브 루트를 `--codex-add-dir`로 허용하지 않는다. 설치
+버전이나 경로가 달라졌다면 위 절대경로를 맹목적으로 사용하지 말고 실제 standalone
+release와 `codex-resources` 존재를 먼저 확인한다.
+
+이전 Harness는 모든 Step과 phase commit을 완료한 뒤 stderr reader의 CP949/UTF-8
+decode 예외를 한 번 출력했다. Phase metadata와 Git commit에는 영향이 없었고 root
+검증은 모두 통과했다. 같은 메시지가 반복되면 제품 실패로 단정하지 말고 executor
+종료 시점, phase status와 Git 이력을 함께 확인한다.
 
 ## 8. 새 세션 시작 절차
 
-먼저 이 문서와 Step 파일을 읽고 현재 상태를 재검증한다.
+먼저 이 문서와 다음 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
+```
+
+그 다음 저장소 상태를 재검증한다.
 
 ```powershell
 git switch dev
@@ -216,62 +285,51 @@ git rev-parse HEAD
 git rev-parse origin/dev
 ```
 
-작업 트리가 clean이고 `dev`가 의도한 기준인지 확인한 뒤 4절의 package 환경 변수를
-설정하고 baseline을 실행한다. 그 다음:
+로컬 `dev`의 미push 10 commit을 보존한다. `origin/dev`에 맞추기 위한 reset, checkout
+또는 rebase를 수행하지 않는다. `docs/HANDOFF.md` 변경이 남아 있으면 먼저 사용자의
+의도대로 커밋하거나 정리한다.
+
+4절의 package 환경 변수를 설정하고 baseline을 실행한다. Child MSBuild 문제가
+재현될 가능성이 있으므로 7절 조건을 적용한 뒤 다음을 실행한다.
 
 ```powershell
-python scripts/execute.py fem-and-beam-kernel
+python scripts/execute.py equation-and-linear-solve
 ```
 
-executor는 `feat-fem-and-beam-kernel` 브랜치를 생성하거나 checkout하고 Step 상태와
-output metadata를 기록한다. 사용자가 명시적으로 요청하지 않은 한 `--push`를
+executor는 `feat-equation-and-linear-solve` 브랜치를 생성하거나 checkout하고 Step
+상태와 output metadata를 기록한다. 사용자가 명시적으로 요청하지 않은 한 `--push`를
 사용하지 않는다.
 
 각 Step은 다음 순서를 지킨다.
 
 1. Step 파일의 필수 문서와 선행 구현을 모두 읽는다.
-2. 성공 기준과 수학 invariant를 명시한다.
+2. 성공 기준과 CSR/constraint/backend invariant를 명시한다.
 3. 실패 테스트를 먼저 작성하고 예상한 이유로 실패함을 확인한다.
 4. 테스트를 통과시키는 최소 production code만 구현한다.
 5. focused test, 전체 CTest, Harness pytest를 실행한다.
 6. Step summary와 output metadata가 실제 결과와 일치하는지 확인한다.
-7. Phase 종료 전 전체 diff를 아키텍처·정식화·테스트 기준으로 review한다.
+7. Phase 종료 전 전체 diff를 아키텍처·수치 계약·resource lifetime 기준으로
+   review한다.
 
-## 9. 알려진 Harness 실행 이력과 복구 주의사항
+## 9. 다음 Phase 완료 조건
 
-이전 Phase에서 production 코드와 무관한 두 executor 문제가 있었다.
-
-- Step 2 child sandbox에서 MSBuild가 사용자 profile의 Microsoft SDK/FileTracker
-  경로를 읽거나 child compiler process를 실행하지 못했다. 같은 checkout의 root
-  PowerShell에서는 정확한 build와 focused/full CTest가 통과했다.
-- Step 3 child Codex 실행은 2026-07-30 usage limit에 도달해
-  `2026-08-05 15:07` 이후 재시도하라는 응답을 냈다. 당시 Step은 root 세션에서
-  TDD와 수용 조건을 직접 수행해 완료했다.
-
-새 세션에서는 이 상태가 여전히 유효하다고 단정하지 말고 Harness를 한 번 정상
-실행해 확인한다. 같은 외부 문제가 반복되면:
-
-1. 저장소 변경 전후 상태와 정확한 실패 stage를 보존한다.
-2. 제품 코드 문제인지 child sandbox/quota 문제인지 root의 동일 명령으로 분리한다.
-3. 수동 fallback 시에도 Step의 TDD, focused/full test, review를 생략하지 않는다.
-4. 실패 output을 성공으로 위장하지 말고 이후 성공 증거와 복구 commit을 분리해
-   metadata와 Git 이력에 남긴다.
-5. 광범위한 사용자 profile 경로를 `--codex-add-dir`로 허용하지 않는다.
-
-## 10. 다음 Phase 완료 조건
-
-- `phases/fem-and-beam-kernel/index.json`의 Step 0~3이 모두 `completed`
-- `phases/index.json`에서 `fem-and-beam-kernel`이 `completed`
-- formulation 문서와 production kernel의 DOF, 좌표계, 부호, 적분 규칙 일치
+- `/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 검증 통과
 - focused test와 전체 CTest 통과
 - Harness pytest가 0개가 아닌 상태로 전체 통과
-- 새 MSVC 경고 없음
-- 강체운동, 해석해, shear-dominant, 세장비와 회전 invariant 검증 통과
-- `fem`/`elements`에 Abaqus 또는 외부 library API 유입 없음
-- 코드 review의 Critical/Important 항목 해결
+- 새 MSVC warning 없음
+- MKL 타입과 handle이 `solvers/linear` adapter 밖으로 노출되지 않음
+- TBB 병렬 조립, analysis pipeline 또는 범위 밖 solver 기능을 선행 구현하지 않음
+- review의 Critical/Important 항목 해결
 - 사용자 선택 전 원격 push나 `dev` 병합을 수행하지 않음
 
 새 세션의 권장 첫 요청:
 
-> `docs/HANDOFF.md`와 `phases/fem-and-beam-kernel/step0.md`부터 `step3.md`를 읽고
-> 현재 baseline을 확인한 뒤 `fem-and-beam-kernel` Phase를 시작해주세요.
+> `docs/HANDOFF.md`와 `phases/equation-and-linear-solve/step0.md`부터 `step2.md`를
+> 읽고 현재 baseline과 child sandbox의 MSBuild 실행 조건을 확인한 뒤
+> `equation-and-linear-solve` Phase를 시작해주세요.

docs/formulation/timoshenko-beam-3d.md

@@ -0,0 +1,237 @@
+# 2절점 3D Timoshenko Beam 정식화
+
+## 1. 범위와 가정
+
+FESA Phase 1 Beam은 소변형, 선형 탄성의 2절점 직선
+isoparametric Timoshenko 요소다. 단면은 도심과 전단중심이 일치하는 주축 단면이며
+\(I_{yz}=0\)이다. 단면 offset, warping, 기하·재료 비선형은 포함하지 않는다.
+FESA는 단위 변환을 하지 않으므로 모든 입력은 하나의 일관 단위계를 사용해야 한다.
+
+정식화의 유한요소 이산화와 수치적분은 Bathe,
+*Finite Element Procedures*, 2판과 Hughes,
+*The Finite Element Method: Linear Static and Dynamic Finite Element
+Analysis*를 따른다. 3차원 isoparametric Beam의 국부 좌표계는 Bathe와
+Bolourchi(1979)의 선형화된 부분을 사용한다. 선택적 감차적분은 Hughes, Taylor와
+Kanoknukulchai(1977)의 원칙을 Beam 전단항에 적용한다. Plate 운동학은 사용하지
+않는다.
+
+## 2. 국부 좌표계와 자유도
+
+절점 1에서 절점 2로 향하는 단위벡터를 \(\mathbf e_x\)로 둔다. 입력 orientation
+벡터에서 \(\mathbf e_x\) 성분을 Gram-Schmidt로 제거하고 정규화한 벡터를
+\(\mathbf e_y\)로 두며,
+
+\[
+\mathbf e_z=\mathbf e_x\times\mathbf e_y
+\]
+
+로 정의한다. 따라서 \((\mathbf e_x,\mathbf e_y,\mathbf e_z)\)는 오른손 직교
+기저다. 길이가 0인 요소, 길이가 0인 orientation, 요소축과 평행한 orientation은
+오류다.
+
+국부 요소 자유도 벡터의 순서는 다음과 같다.
+
+\[
+\mathbf q_l =
+\begin{bmatrix}
+u_1&v_1&w_1&\theta_{x1}&\theta_{y1}&\theta_{z1}&
+u_2&v_2&w_2&\theta_{x2}&\theta_{y2}&\theta_{z2}
+\end{bmatrix}^{T}.
+\]
+
+\(u,v,w\)는 각각 국부 \(x,y,z\) 병진이고
+\(\theta_x,\theta_y,\theta_z\)는 오른손 규칙의 회전 성분이다. 전역 자유도도 절점별
+\([u_x,u_y,u_z,r_x,r_y,r_z]\) 순서를 사용한다.
+
+기준축에서 떨어진 점 \((0,y,z)\)의 선형화된 변위장은
+
+\[
+U_x=u+z\theta_y-y\theta_z,\qquad
+U_y=v-z\theta_x,\qquad
+U_z=w+y\theta_x
+\]
+
+이다.
+
+## 3. 변형률과 부호
+
+일반화 변형률과 대응 단면력의 component 순서는
+
+\[
+\boldsymbol\varepsilon_s =
+\begin{bmatrix}
+\epsilon&\gamma_y&\gamma_z&\kappa_x&\kappa_y&\kappa_z
+\end{bmatrix}^{T},
+\qquad
+\mathbf s =
+\begin{bmatrix}
+N&V_y&V_z&T&M_y&M_z
+\end{bmatrix}^{T}
+\]
+
+로 고정한다. 위 변위장과 오른손 회전 부호로부터
+
+\[
+\begin{aligned}
+\epsilon &= \frac{du}{dx},&
+\gamma_y &= \frac{dv}{dx}-\theta_z,&
+\gamma_z &= \frac{dw}{dx}+\theta_y,\\
+\kappa_x &= \frac{d\theta_x}{dx},&
+\kappa_y &= \frac{d\theta_y}{dx},&
+\kappa_z &= \frac{d\theta_z}{dx}
+\end{aligned}
+\]
+
+를 사용한다. 양의 \(N,V_y,V_z,T,M_y,M_z\)는 위 일반화 변형률과 양의 내부
+일률 \(\delta\boldsymbol\varepsilon_s^T\mathbf s\)을 이루는 방향이다. 이 부호에서
+단면 축응력은
+
+\[
+\sigma_{xx}=E\left(\epsilon+z\kappa_y-y\kappa_z\right)
+\]
+
+이고 \(M_y=\int_A z\sigma_{xx}\,dA\),
+\(M_z=-\int_A y\sigma_{xx}\,dA\)이다.
+
+## 4. 재료와 단면 constitutive matrix
+
+등방성 선형 탄성의 전단계수는 입력 \(E,\nu\)로부터
+
+\[
+G=\frac{E}{2(1+\nu)}
+\]
+
+로 계산한다. 일반화 constitutive 관계는
+
+\[
+\mathbf s=\mathbf D\boldsymbol\varepsilon_s,\qquad
+\mathbf D=
+\operatorname{diag}
+\left(EA,\;GA_{sy},\;GA_{sz},\;GJ,\;EI_y,\;EI_z\right).
+\]
+
+\(A,I_y,I_z,J,A_{sy},A_{sz}\)는 양의 값이어야 한다. 명시적인 전단강성이 없는
+입력은 semantic mapper가 커널 호출 전에
+\(A_{sy}=A_{sz}=5A/6\), `SCF=0`으로 정규화한다. 커널은 Abaqus의 slenderness
+compensation을 추가하지 않는다.
+
+## 5. 자연좌표, 형상함수와 Jacobian
+
+자연좌표는 \(\xi\in[-1,1]\)이고 두 절점의 선형 형상함수는
+
+\[
+N_1(\xi)=\frac{1-\xi}{2},\qquad
+N_2(\xi)=\frac{1+\xi}{2}.
+\]
+
+모든 여섯 component를 같은 형상함수로 보간한다. 직선 요소 길이를 \(L\)이라 하면
+
+\[
+J=\frac{dx}{d\xi}=\frac{L}{2},\qquad
+\frac{dN_a}{dx}=\frac{1}{J}\frac{dN_a}{d\xi}
+\]
+
+이다.
+
+절점 \(a\)의 자유도 순서
+\([u_a,v_a,w_a,\theta_{xa},\theta_{ya},\theta_{za}]\)에 대한
+strain-displacement block은
+
+\[
+\mathbf B_a(\xi)=
+\begin{bmatrix}
+N_{a,x}&0&0&0&0&0\\
+0&N_{a,x}&0&0&0&-N_a\\
+0&0&N_{a,x}&0&N_a&0\\
+0&0&0&N_{a,x}&0&0\\
+0&0&0&0&N_{a,x}&0\\
+0&0&0&0&0&N_{a,x}
+\end{bmatrix},
+\qquad
+\mathbf B=\begin{bmatrix}\mathbf B_1&\mathbf B_2\end{bmatrix}.
+\]
+
+## 6. 선택적 감차적분과 요소 강성
+
+Constitutive matrix를 축·비틀림·굽힘 부분과 전단 부분으로 나눈다.
+
+\[
+\begin{aligned}
+\mathbf D_{ab}&=
+\operatorname{diag}(EA,0,0,GJ,EI_y,EI_z),\\
+\mathbf D_s&=
+\operatorname{diag}(0,GA_{sy},GA_{sz},0,0,0).
+\end{aligned}
+\]
+
+국부 강성은
+
+\[
+\mathbf K_l =
+\sum_{g=1}^{2}
+\mathbf B(\xi_g)^T\mathbf D_{ab}\mathbf B(\xi_g)Jw_g
++
+\sum_{g=1}^{1}
+\mathbf B(\xi_g)^T\mathbf D_s\mathbf B(\xi_g)Jw_g
+\]
+
+로 직접 적분한다. 축·비틀림·굽힘은
+\(\xi=\pm1/\sqrt{3},w=1\)인 2점 Gauss rule을 사용하고, 전단은
+\(\xi=0,w=2\)인 1점 Gauss rule만 사용한다. 전단항을 2점 적분하지 않는다.
+선형 회전장에서 1점 전단 적분은 요소 중앙의 평균 회전을 사용하므로 일정곡률
+굽힘 mode의 불필요한 전단에너지를 제거한다.
+
+`Matrix12`는 `matrix[row][column]`의 12행 12열 고정 크기 저장을 사용한다. 적분은
+각 Gauss point에서 \(\mathbf B^T\mathbf D\mathbf B\)의 상삼각 entry를 고정된
+부동소수점 연산 순서로 누적한 뒤 하삼각에 복사한다.
+
+## 7. 전역 변환
+
+기존 `beam_transformation`이 만드는 \(\mathbf T\)는 각 병진·회전 3-vector를
+전역 성분에서 국부 성분으로 투영한다.
+
+\[
+\mathbf q_l=\mathbf T\mathbf q_g.
+\]
+
+따라서 전역 강성과 에너지 관계는
+
+\[
+\mathbf K_g=\mathbf T^T\mathbf K_l\mathbf T,\qquad
+\frac12\mathbf q_g^T\mathbf K_g\mathbf q_g
+=\frac12\mathbf q_l^T\mathbf K_l\mathbf q_l
+\]
+
+이다.
+
+## 8. 구현과 검증 invariant
+
+커널은 위 적분식으로 local stiffness를 계산하고 변환식으로 global stiffness를
+계산한다. 다음 invariant가 정식화 검증 기준이다.
+
+- local/global stiffness는 roundoff 이내에서 대칭이다.
+- 세 병진과 세 회전 강체 mode의 변형에너지는 0이다.
+- 축 및 비틀림 submatrix는 각각 \(EA/L\), \(GJ/L\)을 재현한다.
+- 일정곡률 단축·이축 굽힘 에너지는 각각
+  \(EI_y\kappa_y^2L/2\), \(EI_z\kappa_z^2L/2\)이다.
+- 일정전단 mode의 에너지는 각각
+  \(GA_{sy}\gamma_y^2L/2\), \(GA_{sz}\gamma_z^2L/2\)이다.
+- 일정곡률 mode의 전단에너지는 slenderness 변화와 무관하게 0이며 전체 에너지는
+  해석값과 일치한다.
+- 요소와 자유도를 같은 강체 회전으로 회전하면 변형에너지가 보존된다.
+
+## 9. 참고문헌과 FESA 차이
+
+- K. J. Bathe, *Finite Element Procedures*, 2nd ed.
+- T. J. R. Hughes, *The Finite Element Method: Linear Static and Dynamic
+  Finite Element Analysis*.
+- K. J. Bathe and S. Bolourchi, “Large Displacement Analysis of
+  Three-Dimensional Beam Structures,” 1979. FESA는 이 문헌의 선형화된 국부
+  좌표와 변환만 사용한다.
+- T. J. R. Hughes, R. L. Taylor, and W. Kanoknukulchai, “A Simple and
+  Efficient Finite Element for Plate Bending,” 1977. FESA는 selective
+  integration 원칙만 사용한다.
+- Abaqus 2024, “Choosing a Beam Element” 및 “BEAM GENERAL SECTION.”
+  FESA의 \(5A/6\) 생략값과 `SCF=0`은 Phase 1 명시적 계약이며 임의 일반 단면의
+  보편적인 Abaqus 기본값이 아니다. Abaqus B31 slenderness compensation도
+  구현하지 않는다.

include/fesa/elements/beam/beam3d2.hpp

@@ -0,0 +1,35 @@
+#pragma once
+
+#include <array>
+#include <optional>
+#include <vector>
+
+#include <fesa/core/diagnostic.hpp>
+#include <fesa/core/vec3.hpp>
+#include <fesa/fem/beam_frame.hpp>
+#include <fesa/model/beam_section.hpp>
+#include <fesa/model/material.hpp>
+
+namespace fesa {
+
+struct Beam3D2Input final {
+    std::array<Vec3, 2> coordinates;
+    IsotropicElastic material;
+    BeamSection section;
+};
+
+struct Beam3D2Contribution final {
+    Matrix12 local_stiffness;
+    Matrix12 global_stiffness;
+    BeamFrame frame;
+};
+
+struct BeamKernelResult final {
+    std::optional<Beam3D2Contribution> contribution;
+    std::vector<Diagnostic> diagnostics;
+};
+
+[[nodiscard]] BeamKernelResult compute_beam3d2(
+    const Beam3D2Input& input);
+
+}  // namespace fesa

include/fesa/fem/beam_frame.hpp

@@ -0,0 +1,32 @@
+#pragma once
+
+#include <array>
+#include <optional>
+#include <vector>
+
+#include <fesa/core/diagnostic.hpp>
+#include <fesa/core/vec3.hpp>
+
+namespace fesa {
+
+using Matrix12 = std::array<std::array<double, 12>, 12>;
+
+struct BeamFrame final {
+    Vec3 ex;
+    Vec3 ey;
+    Vec3 ez;
+};
+
+struct BeamFrameResult final {
+    std::optional<BeamFrame> frame;
+    std::vector<Diagnostic> diagnostics;
+};
+
+[[nodiscard]] BeamFrameResult make_beam_frame(
+    const Vec3& first,
+    const Vec3& second,
+    const Vec3& orientation);
+
+[[nodiscard]] Matrix12 beam_transformation(const BeamFrame& frame);
+
+}  // namespace fesa

include/fesa/fem/dof_manager.hpp

@@ -0,0 +1,52 @@
+#pragma once
+
+#include <array>
+#include <cstddef>
+#include <cstdint>
+#include <map>
+#include <optional>
+#include <span>
+#include <vector>
+
+#include <fesa/model/domain.hpp>
+
+namespace fesa {
+
+enum class NodeDof : std::uint8_t {
+    ux,
+    uy,
+    uz,
+    rx,
+    ry,
+    rz,
+};
+
+struct DofAddress final {
+    NodeId node;
+    NodeDof dof;
+};
+
+class DofManager final {
+public:
+    [[nodiscard]] static DofManager build(const Domain& domain);
+
+    [[nodiscard]] std::size_t full_dof_count() const noexcept;
+    [[nodiscard]] std::size_t free_equation_count() const noexcept;
+    [[nodiscard]] std::optional<std::size_t> equation(
+        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_;
+    std::vector<std::optional<std::size_t>> equations_;
+    std::vector<double> prescribed_values_;
+    std::size_t free_equation_count_ = 0;
+};
+
+}  // namespace fesa

include/fesa/fem/gauss_rule.hpp

@@ -0,0 +1,14 @@
+#pragma once
+
+#include <span>
+
+namespace fesa {
+
+struct GaussPoint1D final {
+    double xi;
+    double weight;
+};
+
+[[nodiscard]] std::span<const GaussPoint1D> gauss_rule_1d(int order);
+
+}  // namespace fesa

include/fesa/fem/line2_shape.hpp

@@ -0,0 +1,13 @@
+#pragma once
+
+#include <array>
+
+namespace fesa {
+
+[[nodiscard]] std::array<double, 2> line2_shape(double xi);
+
+[[nodiscard]] std::array<double, 2> line2_shape_derivative(double xi);
+
+[[nodiscard]] double line2_jacobian(double length);
+
+}  // namespace fesa

phases/fem-and-beam-kernel/index.json

@@ -5,22 +5,36 @@
     {
       "step": 0,
       "name": "quadrature-and-shape-functions",
-      "status": "pending"
+      "status": "completed",
+      "summary": "Added fixed 1D Gauss rules, Line2 shape primitives, and a validated length/2 Jacobian with invariant and invalid-input tests.",
+      "started_at": "2026-07-31T01:37:23+0900",
+      "completed_at": "2026-07-31T01:48:32+0900"
     },
     {
       "step": 1,
       "name": "dof-manager",
-      "status": "pending"
+      "status": "completed",
+      "summary": "Added deterministic six-DOF equation numbering, prescribed-value reconstruction, and 12-DOF Beam element mapping with validation tests.",
+      "started_at": "2026-07-31T01:48:33+0900",
+      "completed_at": "2026-07-31T01:55:21+0900"
     },
     {
       "step": 2,
       "name": "beam-local-frame",
-      "status": "pending"
+      "status": "completed",
+      "summary": "Added scale-aware right-handed Beam frames and fixed 12x12 global-to-local transformations with degeneracy, orthonormality, rotation, and energy-invariance tests.",
+      "started_at": "2026-07-31T01:55:21+0900",
+      "completed_at": "2026-07-31T02:08:27+0900"
     },
     {
       "step": 3,
       "name": "timoshenko-stiffness-kernel",
-      "status": "pending"
+      "status": "completed",
+      "summary": "Documented and implemented selectively integrated Beam3D2 local/global stiffness with rigid-body, analytical, slenderness, and rotation-invariance tests.",
+      "started_at": "2026-07-31T02:08:27+0900",
+      "completed_at": "2026-07-31T02:24:47+0900"
     }
-  ]
-}
+  ],
+  "created_at": "2026-07-31T01:37:23+0900",
+  "completed_at": "2026-07-31T02:24:48+0900"
+}

phases/index.json

@@ -12,7 +12,8 @@
     },
     {
       "dir": "fem-and-beam-kernel",
-      "status": "pending"
+      "status": "completed",
+      "completed_at": "2026-07-31T02:24:48+0900"
     },
     {
       "dir": "equation-and-linear-solve",
@@ -43,4 +44,4 @@
       "status": "pending"
     }
   ]
-}
+}

src/fesa/elements/beam/beam3d2.cpp

@@ -0,0 +1,244 @@
+#include <fesa/elements/beam/beam3d2.hpp>
+
+#include <array>
+#include <cmath>
+#include <cstddef>
+#include <optional>
+#include <string>
+#include <utility>
+
+#include <fesa/fem/gauss_rule.hpp>
+#include <fesa/fem/line2_shape.hpp>
+
+namespace fesa {
+namespace {
+
+using StrainMatrix =
+    std::array<std::array<double, 12>, 6>;
+
+BeamKernelResult error_result(std::string code, std::string message) {
+    BeamKernelResult result;
+    result.diagnostics.push_back({
+        DiagnosticStage::model,
+        Severity::error,
+        std::move(code),
+        std::move(message),
+        std::nullopt,
+    });
+    return result;
+}
+
+bool is_positive_finite(const double value) {
+    return std::isfinite(value) && value > 0.0;
+}
+
+std::optional<BeamKernelResult> validate_properties(
+    const Beam3D2Input& input) {
+    if (!std::isfinite(input.material.young) ||
+        !std::isfinite(input.material.poisson)) {
+        return error_result(
+            "model.nonfinite_value",
+            "Beam kernel requires finite elastic constants.");
+    }
+    if (input.material.young <= 0.0 ||
+        input.material.poisson <= -1.0 ||
+        input.material.poisson >= 0.5) {
+        return error_result(
+            "model.invalid_material",
+            "Beam kernel requires E > 0 and -1 < nu < 0.5.");
+    }
+
+    if (!is_positive_finite(input.section.area) ||
+        !is_positive_finite(input.section.iy) ||
+        !is_positive_finite(input.section.iz) ||
+        !is_positive_finite(input.section.torsion_j) ||
+        !is_positive_finite(input.section.shear_area_y) ||
+        !is_positive_finite(input.section.shear_area_z)) {
+        return error_result(
+            "model.invalid_section",
+            "Beam kernel requires positive finite section properties.");
+    }
+    return std::nullopt;
+}
+
+StrainMatrix strain_matrix(
+    const double xi,
+    const double jacobian) {
+    StrainMatrix strain{};
+    const auto shape = line2_shape(xi);
+    const auto natural_derivative = line2_shape_derivative(xi);
+
+    for (std::size_t node = 0; node < shape.size(); ++node) {
+        const std::size_t offset = node * 6;
+        const double derivative = natural_derivative[node] / jacobian;
+        strain[0][offset] = derivative;
+        strain[1][offset + 1] = derivative;
+        strain[1][offset + 5] = -shape[node];
+        strain[2][offset + 2] = derivative;
+        strain[2][offset + 4] = shape[node];
+        strain[3][offset + 3] = derivative;
+        strain[4][offset + 4] = derivative;
+        strain[5][offset + 5] = derivative;
+    }
+    return strain;
+}
+
+template <std::size_t ComponentCount>
+void integrate_components(
+    Matrix12& stiffness,
+    const std::span<const GaussPoint1D> rule,
+    const std::array<std::size_t, ComponentCount>& components,
+    const std::array<double, 6>& constitutive,
+    const double jacobian) {
+    for (const GaussPoint1D& point : rule) {
+        const StrainMatrix strain = strain_matrix(point.xi, jacobian);
+        const double integration_weight = jacobian * point.weight;
+        for (std::size_t row = 0; row < stiffness.size(); ++row) {
+            for (std::size_t column = row;
+                 column < stiffness[row].size();
+                 ++column) {
+                double entry = 0.0;
+                for (const std::size_t component : components) {
+                    entry +=
+                        strain[component][row] *
+                        constitutive[component] *
+                        strain[component][column];
+                }
+                stiffness[row][column] +=
+                    entry * integration_weight;
+            }
+        }
+    }
+}
+
+void mirror_upper_triangle(Matrix12& matrix) {
+    for (std::size_t row = 0; row < matrix.size(); ++row) {
+        for (std::size_t column = row + 1;
+             column < matrix[row].size();
+             ++column) {
+            matrix[column][row] = matrix[row][column];
+        }
+    }
+}
+
+Matrix12 transform_stiffness(
+    const Matrix12& local,
+    const Matrix12& transformation) {
+    Matrix12 global{};
+    for (std::size_t row = 0; row < global.size(); ++row) {
+        for (std::size_t column = row;
+             column < global[row].size();
+             ++column) {
+            double entry = 0.0;
+            for (std::size_t local_row = 0;
+                 local_row < local.size();
+                 ++local_row) {
+                for (std::size_t local_column = 0;
+                     local_column < local[local_row].size();
+                     ++local_column) {
+                    entry +=
+                        transformation[local_row][row] *
+                        local[local_row][local_column] *
+                        transformation[local_column][column];
+                }
+            }
+            global[row][column] = entry;
+        }
+    }
+    mirror_upper_triangle(global);
+    return global;
+}
+
+bool is_finite(const Matrix12& matrix) {
+    for (const auto& row : matrix) {
+        for (const double value : row) {
+            if (!std::isfinite(value)) {
+                return false;
+            }
+        }
+    }
+    return true;
+}
+
+}  // namespace
+
+BeamKernelResult compute_beam3d2(const Beam3D2Input& input) {
+    if (const auto validation = validate_properties(input);
+        validation.has_value()) {
+        return *validation;
+    }
+
+    BeamFrameResult frame_result = make_beam_frame(
+        input.coordinates[0],
+        input.coordinates[1],
+        input.section.orientation);
+    if (!frame_result.frame.has_value()) {
+        return {std::nullopt, std::move(frame_result.diagnostics)};
+    }
+
+    const Vec3 axis{
+        input.coordinates[1].x - input.coordinates[0].x,
+        input.coordinates[1].y - input.coordinates[0].y,
+        input.coordinates[1].z - input.coordinates[0].z,
+    };
+    const double length = std::hypot(axis.x, axis.y, axis.z);
+    const double jacobian = length / 2.0;
+    if (!std::isfinite(jacobian) || jacobian <= 0.0) {
+        return error_result(
+            "model.zero_length_element",
+            "Beam kernel requires a representable positive Jacobian.");
+    }
+
+    const double shear_modulus =
+        input.material.young / (2.0 * (1.0 + input.material.poisson));
+    const std::array<double, 6> constitutive{
+        input.material.young * input.section.area,
+        shear_modulus * input.section.shear_area_y,
+        shear_modulus * input.section.shear_area_z,
+        shear_modulus * input.section.torsion_j,
+        input.material.young * input.section.iy,
+        input.material.young * input.section.iz,
+    };
+    for (const double value : constitutive) {
+        if (!is_positive_finite(value)) {
+            return error_result(
+                "model.nonfinite_value",
+                "Beam constitutive stiffness is not finite and positive.");
+        }
+    }
+
+    Matrix12 local{};
+    integrate_components(
+        local,
+        gauss_rule_1d(2),
+        std::array<std::size_t, 4>{0, 3, 4, 5},
+        constitutive,
+        jacobian);
+    integrate_components(
+        local,
+        gauss_rule_1d(1),
+        std::array<std::size_t, 2>{1, 2},
+        constitutive,
+        jacobian);
+    mirror_upper_triangle(local);
+
+    const Matrix12 transformation =
+        beam_transformation(*frame_result.frame);
+    Matrix12 global = transform_stiffness(local, transformation);
+    if (!is_finite(local) || !is_finite(global)) {
+        return error_result(
+            "model.nonfinite_value",
+            "Beam stiffness contains a nonfinite value.");
+    }
+
+    return {
+        Beam3D2Contribution{
+            std::move(local),
+            std::move(global),
+            *frame_result.frame,
+        },
+        {},
+    };
+}
+
+}  // namespace fesa

src/fesa/fem/beam_frame.cpp

@@ -0,0 +1,146 @@
+#include <fesa/fem/beam_frame.hpp>
+
+#include <algorithm>
+#include <array>
+#include <cmath>
+#include <cstddef>
+#include <limits>
+#include <string>
+#include <utility>
+
+namespace fesa {
+namespace {
+
+constexpr double kParallelTolerance =
+    64.0 * std::numeric_limits<double>::epsilon();
+
+double length(const Vec3 value) {
+    return std::hypot(value.x, value.y, value.z);
+}
+
+double max_abs_component(const Vec3 value) {
+    return std::max(
+        std::abs(value.x),
+        std::max(std::abs(value.y), std::abs(value.z)));
+}
+
+double dot(const Vec3 first, const Vec3 second) {
+    return first.x * second.x + first.y * second.y + first.z * second.z;
+}
+
+Vec3 cross(const Vec3 first, const Vec3 second) {
+    return {
+        first.y * second.z - first.z * second.y,
+        first.z * second.x - first.x * second.z,
+        first.x * second.y - first.y * second.x,
+    };
+}
+
+Vec3 normalized(const Vec3 value) {
+    const double scale = max_abs_component(value);
+    const Vec3 scaled{
+        value.x / scale,
+        value.y / scale,
+        value.z / scale,
+    };
+    const double scaled_length = length(scaled);
+    return {
+        scaled.x / scaled_length,
+        scaled.y / scaled_length,
+        scaled.z / scaled_length,
+    };
+}
+
+BeamFrameResult error_result(std::string code, std::string message) {
+    BeamFrameResult result;
+    result.diagnostics.push_back({
+        DiagnosticStage::model,
+        Severity::error,
+        std::move(code),
+        std::move(message),
+        std::nullopt,
+    });
+    return result;
+}
+
+}  // namespace
+
+BeamFrameResult make_beam_frame(
+    const Vec3& first,
+    const Vec3& second,
+    const Vec3& orientation) {
+    if (!is_finite(first) || !is_finite(second)) {
+        return error_result(
+            "model.nonfinite_value",
+            "Beam frame requires finite node coordinates.");
+    }
+
+    const Vec3 axis{
+        second.x - first.x,
+        second.y - first.y,
+        second.z - first.z,
+    };
+    if (!is_finite(axis)) {
+        return error_result(
+            "model.nonfinite_value",
+            "Beam frame axis is not representable as a finite vector.");
+    }
+
+    if (max_abs_component(axis) == 0.0) {
+        return error_result(
+            "model.zero_length_element",
+            "Beam frame requires distinct node coordinates.");
+    }
+
+    if (!is_finite(orientation)) {
+        return error_result(
+            "model.nonfinite_value",
+            "Beam frame requires a finite orientation vector.");
+    }
+
+    if (max_abs_component(orientation) == 0.0) {
+        return error_result(
+            "model.invalid_orientation",
+            "Beam frame requires a nonzero orientation vector.");
+    }
+
+    const Vec3 ex = normalized(axis);
+    const Vec3 orientation_unit = normalized(orientation);
+    const double axial_projection = dot(orientation_unit, ex);
+    const Vec3 transverse{
+        orientation_unit.x - axial_projection * ex.x,
+        orientation_unit.y - axial_projection * ex.y,
+        orientation_unit.z - axial_projection * ex.z,
+    };
+    const double transverse_length = length(transverse);
+    if (transverse_length <= kParallelTolerance) {
+        return error_result(
+            "model.invalid_orientation",
+            "Beam orientation is parallel to the element axis.");
+    }
+
+    const Vec3 ey = normalized(transverse);
+    const Vec3 ez = cross(ex, ey);
+    return {
+        BeamFrame{ex, ey, ez},
+        {},
+    };
+}
+
+Matrix12 beam_transformation(const BeamFrame& frame) {
+    Matrix12 transformation{};
+    const std::array<Vec3, 3> basis{frame.ex, frame.ey, frame.ez};
+    constexpr std::array<std::size_t, 4> block_offsets{0, 3, 6, 9};
+
+    for (const std::size_t offset : block_offsets) {
+        for (std::size_t row = 0; row < basis.size(); ++row) {
+            transformation[offset + row][offset] = basis[row].x;
+            transformation[offset + row][offset + 1] = basis[row].y;
+            transformation[offset + row][offset + 2] = basis[row].z;
+        }
+    }
+
+    return transformation;
+}
+
+}  // namespace fesa

src/fesa/fem/dof_manager.cpp

@@ -0,0 +1,111 @@
+#include <fesa/fem/dof_manager.hpp>
+
+#include <algorithm>
+#include <stdexcept>
+#include <vector>
+
+namespace fesa {
+namespace {
+
+constexpr std::size_t dofs_per_node = 6;
+
+}  // namespace
+
+DofManager DofManager::build(const Domain& domain) {
+    DofManager manager;
+
+    std::vector<NodeId> node_ids;
+    node_ids.reserve(domain.nodes().size());
+    for (const Node& node : domain.nodes()) {
+        node_ids.push_back(node.id);
+    }
+    std::ranges::sort(node_ids);
+
+    for (std::size_t index = 0; index < node_ids.size(); ++index) {
+        manager.node_full_dof_bases_.emplace(
+            node_ids[index].value(), index * dofs_per_node);
+    }
+
+    const std::size_t full_dof_count = node_ids.size() * dofs_per_node;
+    manager.equations_.resize(full_dof_count);
+    manager.prescribed_values_.resize(full_dof_count);
+    std::vector<bool> constrained(full_dof_count, false);
+
+    for (const PrescribedDof& prescribed : domain.step().prescribed_dofs) {
+        const NodeDof dof =
+            static_cast<NodeDof>(prescribed.dof - std::uint8_t{1});
+        const std::size_t full =
+            manager.full_dof({prescribed.node, dof});
+        constrained[full] = true;
+        manager.prescribed_values_[full] = prescribed.value;
+    }
+
+    for (std::size_t full = 0; full < full_dof_count; ++full) {
+        if (!constrained[full]) {
+            manager.equations_[full] = manager.free_equation_count_;
+            ++manager.free_equation_count_;
+        }
+    }
+
+    return manager;
+}
+
+std::size_t DofManager::full_dof_count() const noexcept {
+    return equations_.size();
+}
+
+std::size_t DofManager::free_equation_count() const noexcept {
+    return free_equation_count_;
+}
+
+std::optional<std::size_t> DofManager::equation(
+    const DofAddress address) const {
+    return equations_[full_dof(address)];
+}
+
+std::array<std::size_t, 12> DofManager::element_full_dofs(
+    const BeamElement& element) const {
+    std::array<std::size_t, 12> full_dofs{};
+    for (std::size_t node = 0; node < element.nodes.size(); ++node) {
+        const std::size_t base = node_full_dof_base(element.nodes[node]);
+        for (std::size_t component = 0; component < dofs_per_node;
+             ++component) {
+            full_dofs[node * dofs_per_node + component] = base + component;
+        }
+    }
+    return full_dofs;
+}
+
+std::vector<double> DofManager::reconstruct_full(
+    const std::span<const double> reduced) const {
+    if (reduced.size() != free_equation_count_) {
+        throw std::invalid_argument{
+            "Reduced vector size must equal the free equation count."};
+    }
+
+    std::vector<double> full = prescribed_values_;
+    for (std::size_t index = 0; index < equations_.size(); ++index) {
+        if (equations_[index].has_value()) {
+            full[index] = reduced[*equations_[index]];
+        }
+    }
+    return full;
+}
+
+std::size_t DofManager::full_dof(const DofAddress address) const {
+    const auto component = static_cast<std::uint8_t>(address.dof);
+    if (component >= dofs_per_node) {
+        throw std::invalid_argument{"Node DOF must be one of ux through rz."};
+    }
+    return node_full_dof_base(address.node) + component;
+}
+
+std::size_t DofManager::node_full_dof_base(const NodeId node) const {
+    const auto found = node_full_dof_bases_.find(node.value());
+    if (found == node_full_dof_bases_.end()) {
+        throw std::out_of_range{"Node ID is not present in the DofManager."};
+    }
+    return found->second;
+}
+
+}  // namespace fesa

src/fesa/fem/gauss_rule.cpp

@@ -0,0 +1,33 @@
+#include <fesa/fem/gauss_rule.hpp>
+
+#include <array>
+#include <stdexcept>
+
+namespace fesa {
+
+namespace {
+
+constexpr std::array<GaussPoint1D, 1> order_one{{
+    {0.0, 2.0},
+}};
+
+constexpr double inverse_sqrt_three = 0.57735026918962576451;
+constexpr std::array<GaussPoint1D, 2> order_two{{
+    {-inverse_sqrt_three, 1.0},
+    {inverse_sqrt_three, 1.0},
+}};
+
+}  // namespace
+
+std::span<const GaussPoint1D> gauss_rule_1d(const int order) {
+    switch (order) {
+        case 1:
+            return order_one;
+        case 2:
+            return order_two;
+        default:
+            throw std::invalid_argument{"Gauss rule order must be 1 or 2."};
+    }
+}
+
+}  // namespace fesa

src/fesa/fem/line2_shape.cpp

@@ -0,0 +1,31 @@
+#include <fesa/fem/line2_shape.hpp>
+
+#include <cmath>
+#include <stdexcept>
+
+namespace fesa {
+
+std::array<double, 2> line2_shape(const double xi) {
+    return {0.5 * (1.0 - xi), 0.5 * (1.0 + xi)};
+}
+
+std::array<double, 2> line2_shape_derivative(const double) {
+    return {-0.5, 0.5};
+}
+
+double line2_jacobian(const double length) {
+    if (!std::isfinite(length) || length <= 0.0) {
+        throw std::invalid_argument{
+            "Line2 Jacobian requires a positive finite length."};
+    }
+
+    const double jacobian = length / 2.0;
+    if (jacobian == 0.0) {
+        throw std::invalid_argument{
+            "Line2 Jacobian must be representable as a positive double."};
+    }
+
+    return jacobian;
+}
+
+}  // namespace fesa

tests/CMakeLists.txt

@@ -170,3 +170,104 @@ add_test(
     COMMAND "$<TARGET_FILE:fesa_deck_to_domain_tests>"
         --gtest_filter=ActiveInstance.*
 )
+
+add_executable(fesa_fem_primitives_tests
+    unit/fem/beam_frame_test.cpp
+    unit/fem/dof_manager_test.cpp
+    unit/fem/gauss_rule_test.cpp
+    unit/fem/line2_shape_test.cpp
+)
+
+target_compile_features(fesa_fem_primitives_tests PRIVATE cxx_std_20)
+target_compile_options(
+    fesa_fem_primitives_tests
+    PRIVATE
+        /W4
+        /permissive-
+        /EHsc
+)
+
+target_link_libraries(fesa_fem_primitives_tests
+    PRIVATE
+        fesa_core
+        GTest::gtest_main
+)
+
+add_test(
+    NAME Quadrature
+    COMMAND "$<TARGET_FILE:fesa_fem_primitives_tests>"
+        --gtest_filter=Quadrature.*
+)
+
+add_test(
+    NAME ShapeFunction
+    COMMAND "$<TARGET_FILE:fesa_fem_primitives_tests>"
+        --gtest_filter=ShapeFunction.*
+)
+
+add_test(
+    NAME Jacobian
+    COMMAND "$<TARGET_FILE:fesa_fem_primitives_tests>"
+        --gtest_filter=Jacobian.*
+)
+
+add_test(
+    NAME DofManager
+    COMMAND "$<TARGET_FILE:fesa_fem_primitives_tests>"
+        --gtest_filter=DofManager.*
+)
+
+add_test(
+    NAME EquationNumbering
+    COMMAND "$<TARGET_FILE:fesa_fem_primitives_tests>"
+        --gtest_filter=EquationNumbering.*
+)
+
+add_test(
+    NAME BeamFrame
+    COMMAND "$<TARGET_FILE:fesa_fem_primitives_tests>"
+        --gtest_filter=BeamFrame.*
+)
+
+add_test(
+    NAME BeamTransformation
+    COMMAND "$<TARGET_FILE:fesa_fem_primitives_tests>"
+        --gtest_filter=BeamTransformation.*
+)
+
+add_executable(fesa_beam3d2_tests
+    unit/elements/beam3d2_test.cpp
+)
+
+target_compile_features(fesa_beam3d2_tests PRIVATE cxx_std_20)
+target_compile_options(
+    fesa_beam3d2_tests
+    PRIVATE
+        /W4
+        /permissive-
+        /EHsc
+)
+
+target_link_libraries(fesa_beam3d2_tests
+    PRIVATE
+        fesa_core
+        GTest::gtest_main
+)
+
+add_test(
+    NAME Beam3D2
+    COMMAND "$<TARGET_FILE:fesa_beam3d2_tests>"
+        --gtest_filter=Beam3D2.*
+)
+
+add_test(
+    NAME Timoshenko
+    COMMAND "$<TARGET_FILE:fesa_beam3d2_tests>"
+        --gtest_filter=Timoshenko.*
+)
+
+add_test(
+    NAME RigidBody
+    COMMAND "$<TARGET_FILE:fesa_beam3d2_tests>"
+        --gtest_filter=RigidBody.*
+)

tests/unit/elements/beam3d2_test.cpp

@@ -0,0 +1,372 @@
+#include <fesa/elements/beam/beam3d2.hpp>
+
+#include <algorithm>
+#include <array>
+#include <cmath>
+#include <cstddef>
+#include <limits>
+
+#include <gtest/gtest.h>
+
+namespace {
+
+constexpr double kYoung = 210.0;
+constexpr double kPoisson = 0.3;
+constexpr double kArea = 0.4;
+constexpr double kIy = 0.03;
+constexpr double kIz = 0.05;
+constexpr double kTorsionJ = 0.02;
+constexpr double kShearAreaY = 0.25;
+constexpr double kShearAreaZ = 0.2;
+
+fesa::Beam3D2Input make_x_axis_input(const double length) {
+    return {
+        {{{0.0, 0.0, 0.0}, {length, 0.0, 0.0}}},
+        {fesa::MaterialId{1}, "elastic", kYoung, kPoisson},
+        {
+            fesa::SectionId{1},
+            "section",
+            kArea,
+            kIy,
+            kIz,
+            kTorsionJ,
+            kShearAreaY,
+            kShearAreaZ,
+            fesa::ShearPropertySource::input,
+            {0.0, 1.0, 0.0},
+            {},
+        },
+    };
+}
+
+double quadratic_form(
+    const fesa::Matrix12& matrix,
+    const std::array<double, 12>& vector) {
+    double result = 0.0;
+    for (std::size_t row = 0; row < matrix.size(); ++row) {
+        for (std::size_t column = 0; column < matrix[row].size(); ++column) {
+            result += vector[row] * matrix[row][column] * vector[column];
+        }
+    }
+    return result;
+}
+
+double strain_energy(
+    const fesa::Matrix12& matrix,
+    const std::array<double, 12>& vector) {
+    return 0.5 * quadratic_form(matrix, vector);
+}
+
+double maximum_abs_entry(const fesa::Matrix12& matrix) {
+    double maximum = 0.0;
+    for (const auto& row : matrix) {
+        for (const double value : row) {
+            maximum = std::max(maximum, std::abs(value));
+        }
+    }
+    return maximum;
+}
+
+double squared_norm(const std::array<double, 12>& vector) {
+    double result = 0.0;
+    for (const double value : vector) {
+        result += value * value;
+    }
+    return result;
+}
+
+double dot(const fesa::Vec3 first, const fesa::Vec3 second) {
+    return first.x * second.x + first.y * second.y + first.z * second.z;
+}
+
+fesa::Vec3 cross(const fesa::Vec3 first, const fesa::Vec3 second) {
+    return {
+        first.y * second.z - first.z * second.y,
+        first.z * second.x - first.x * second.z,
+        first.x * second.y - first.y * second.x,
+    };
+}
+
+fesa::Vec3 rotate_about_axis(
+    const fesa::Vec3 value,
+    const fesa::Vec3 unit_axis,
+    const double angle) {
+    const double cosine = std::cos(angle);
+    const double sine = std::sin(angle);
+    const fesa::Vec3 axis_cross_value = cross(unit_axis, value);
+    const double axis_projection = dot(unit_axis, value);
+    return {
+        value.x * cosine + axis_cross_value.x * sine +
+            unit_axis.x * axis_projection * (1.0 - cosine),
+        value.y * cosine + axis_cross_value.y * sine +
+            unit_axis.y * axis_projection * (1.0 - cosine),
+        value.z * cosine + axis_cross_value.z * sine +
+            unit_axis.z * axis_projection * (1.0 - cosine),
+    };
+}
+
+std::array<double, 12> rotate_dofs(
+    const std::array<double, 12>& values,
+    const fesa::Vec3 unit_axis,
+    const double angle) {
+    std::array<double, 12> rotated{};
+    for (const std::size_t offset : std::array<std::size_t, 4>{0, 3, 6, 9}) {
+        const fesa::Vec3 value{
+            values[offset],
+            values[offset + 1],
+            values[offset + 2],
+        };
+        const fesa::Vec3 rotated_value =
+            rotate_about_axis(value, unit_axis, angle);
+        rotated[offset] = rotated_value.x;
+        rotated[offset + 1] = rotated_value.y;
+        rotated[offset + 2] = rotated_value.z;
+    }
+    return rotated;
+}
+
+void expect_relative_near(
+    const double actual,
+    const double expected,
+    const double relative_tolerance =
+        512.0 * std::numeric_limits<double>::epsilon()) {
+    const double scale = std::max(1.0, std::abs(expected));
+    EXPECT_NEAR(actual, expected, relative_tolerance * scale);
+}
+
+TEST(Beam3D2, ProducesFiniteSymmetricLocalAndGlobalStiffness) {
+    const auto result = fesa::compute_beam3d2(make_x_axis_input(2.5));
+
+    ASSERT_TRUE(result.contribution.has_value());
+    EXPECT_TRUE(result.diagnostics.empty());
+    const auto& contribution = *result.contribution;
+    for (std::size_t row = 0; row < 12; ++row) {
+        for (std::size_t column = 0; column < 12; ++column) {
+            EXPECT_TRUE(std::isfinite(
+                contribution.local_stiffness[row][column]));
+            EXPECT_TRUE(std::isfinite(
+                contribution.global_stiffness[row][column]));
+            EXPECT_DOUBLE_EQ(
+                contribution.local_stiffness[row][column],
+                contribution.local_stiffness[column][row]);
+            EXPECT_DOUBLE_EQ(
+                contribution.global_stiffness[row][column],
+                contribution.global_stiffness[column][row]);
+        }
+    }
+}
+
+TEST(RigidBody, SixIndependentModesHaveZeroStrainEnergy) {
+    constexpr double length = 3.0;
+    const auto result = fesa::compute_beam3d2(make_x_axis_input(length));
+    ASSERT_TRUE(result.contribution.has_value());
+    const auto& stiffness = result.contribution->local_stiffness;
+
+    const std::array<std::array<double, 12>, 6> modes{{
+        {1.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+         1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
+        {0.0, 1.0, 0.0, 0.0, 0.0, 0.0,
+         0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
+        {0.0, 0.0, 1.0, 0.0, 0.0, 0.0,
+         0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
+        {0.0, 0.0, 0.0, 1.0, 0.0, 0.0,
+         0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
+        {0.0, 0.0, 0.0, 0.0, 1.0, 0.0,
+         0.0, 0.0, -length, 0.0, 1.0, 0.0},
+        {0.0, 0.0, 0.0, 0.0, 0.0, 1.0,
+         0.0, length, 0.0, 0.0, 0.0, 1.0},
+    }};
+
+    for (std::size_t mode = 0; mode < modes.size(); ++mode) {
+        const double roundoff_bound =
+            4096.0 * std::numeric_limits<double>::epsilon() *
+            maximum_abs_entry(stiffness) * squared_norm(modes[mode]);
+        EXPECT_NEAR(
+            quadratic_form(stiffness, modes[mode]),
+            0.0,
+            roundoff_bound)
+            << "rigid mode " << mode;
+    }
+}
+
+TEST(Timoshenko, ReproducesAnalyticalAxialSubmatrix) {
+    constexpr double length = 2.5;
+    const auto result = fesa::compute_beam3d2(make_x_axis_input(length));
+    ASSERT_TRUE(result.contribution.has_value());
+    const auto& stiffness = result.contribution->local_stiffness;
+    const double expected = kYoung * kArea / length;
+
+    expect_relative_near(stiffness[0][0], expected);
+    expect_relative_near(stiffness[0][6], -expected);
+    expect_relative_near(stiffness[6][0], -expected);
+    expect_relative_near(stiffness[6][6], expected);
+}
+
+TEST(Timoshenko, ReproducesAnalyticalTorsionalSubmatrix) {
+    constexpr double length = 2.5;
+    const auto result = fesa::compute_beam3d2(make_x_axis_input(length));
+    ASSERT_TRUE(result.contribution.has_value());
+    const auto& stiffness = result.contribution->local_stiffness;
+    const double shear_modulus = kYoung / (2.0 * (1.0 + kPoisson));
+    const double expected = shear_modulus * kTorsionJ / length;
+
+    expect_relative_near(stiffness[3][3], expected);
+    expect_relative_near(stiffness[3][9], -expected);
+    expect_relative_near(stiffness[9][3], -expected);
+    expect_relative_near(stiffness[9][9], expected);
+}
+
+TEST(Timoshenko, ReproducesConstantCurvatureEnergyAboutLocalY) {
+    constexpr double length = 2.5;
+    constexpr double curvature = 0.4;
+    const auto result = fesa::compute_beam3d2(make_x_axis_input(length));
+    ASSERT_TRUE(result.contribution.has_value());
+    std::array<double, 12> displacement{};
+    displacement[4] = -0.5 * curvature * length;
+    displacement[10] = 0.5 * curvature * length;
+
+    const double expected =
+        0.5 * kYoung * kIy * curvature * curvature * length;
+    expect_relative_near(
+        strain_energy(result.contribution->local_stiffness, displacement),
+        expected);
+}
+
+TEST(Timoshenko, ReproducesConstantCurvatureEnergyAboutLocalZ) {
+    constexpr double length = 2.5;
+    constexpr double curvature = 0.4;
+    const auto result = fesa::compute_beam3d2(make_x_axis_input(length));
+    ASSERT_TRUE(result.contribution.has_value());
+    std::array<double, 12> displacement{};
+    displacement[5] = -0.5 * curvature * length;
+    displacement[11] = 0.5 * curvature * length;
+
+    const double expected =
+        0.5 * kYoung * kIz * curvature * curvature * length;
+    expect_relative_near(
+        strain_energy(result.contribution->local_stiffness, displacement),
+        expected);
+}
+
+TEST(Timoshenko, ReproducesConstantShearEnergyInLocalY) {
+    constexpr double length = 2.5;
+    constexpr double shear_strain = 0.04;
+    const auto result = fesa::compute_beam3d2(make_x_axis_input(length));
+    ASSERT_TRUE(result.contribution.has_value());
+    std::array<double, 12> displacement{};
+    displacement[1] = -0.5 * shear_strain * length;
+    displacement[7] = 0.5 * shear_strain * length;
+
+    const double shear_modulus = kYoung / (2.0 * (1.0 + kPoisson));
+    const double expected =
+        0.5 * shear_modulus * kShearAreaY *
+        shear_strain * shear_strain * length;
+    expect_relative_near(
+        strain_energy(result.contribution->local_stiffness, displacement),
+        expected);
+}
+
+TEST(Timoshenko, ReproducesConstantShearEnergyInLocalZ) {
+    constexpr double length = 2.5;
+    constexpr double shear_strain = 0.04;
+    const auto result = fesa::compute_beam3d2(make_x_axis_input(length));
+    ASSERT_TRUE(result.contribution.has_value());
+    std::array<double, 12> displacement{};
+    displacement[2] = -0.5 * shear_strain * length;
+    displacement[8] = 0.5 * shear_strain * length;
+
+    const double shear_modulus = kYoung / (2.0 * (1.0 + kPoisson));
+    const double expected =
+        0.5 * shear_modulus * kShearAreaZ *
+        shear_strain * shear_strain * length;
+    expect_relative_near(
+        strain_energy(result.contribution->local_stiffness, displacement),
+        expected);
+}
+
+TEST(Timoshenko, AvoidsShearLockingAcrossSlendernessSweep) {
+    constexpr double length = 2.0;
+    constexpr double curvature = 0.2;
+    for (const double slenderness : std::array{2.0, 10.0, 100.0, 1000.0}) {
+        fesa::Beam3D2Input input = make_x_axis_input(length);
+        const double side = length / slenderness;
+        input.section.area = side * side;
+        input.section.iy = side * side * side * side / 12.0;
+        input.section.iz = input.section.iy;
+        input.section.shear_area_y = 5.0 * input.section.area / 6.0;
+        input.section.shear_area_z = input.section.shear_area_y;
+
+        const auto result = fesa::compute_beam3d2(input);
+        ASSERT_TRUE(result.contribution.has_value());
+        std::array<double, 12> displacement{};
+        displacement[5] = -0.5 * curvature * length;
+        displacement[11] = 0.5 * curvature * length;
+
+        const double expected =
+            0.5 * kYoung * input.section.iz *
+            curvature * curvature * length;
+        const double relative_tolerance =
+            4096.0 * std::numeric_limits<double>::epsilon() *
+            slenderness * slenderness;
+        EXPECT_NEAR(
+            strain_energy(
+                result.contribution->local_stiffness,
+                displacement),
+            expected,
+            relative_tolerance * expected)
+            << "L/h=" << slenderness;
+    }
+}
+
+TEST(Beam3D2, PreservesGlobalEnergyUnderRigidCoordinateRotation) {
+    fesa::Beam3D2Input original = make_x_axis_input(2.5);
+    original.coordinates = {{
+        {1.0, -2.0, 0.5},
+        {3.0, -0.5, 1.25},
+    }};
+    original.section.orientation = {-1.0, 2.0, 3.0};
+
+    const double inverse_sqrt_fourteen = 1.0 / std::sqrt(14.0);
+    const fesa::Vec3 rotation_axis{
+        inverse_sqrt_fourteen,
+        2.0 * inverse_sqrt_fourteen,
+        3.0 * inverse_sqrt_fourteen,
+    };
+    constexpr double angle = 0.73;
+    fesa::Beam3D2Input rotated = original;
+    for (std::size_t node = 0; node < rotated.coordinates.size(); ++node) {
+        rotated.coordinates[node] = rotate_about_axis(
+            original.coordinates[node],
+            rotation_axis,
+            angle);
+    }
+    rotated.section.orientation = rotate_about_axis(
+        original.section.orientation,
+        rotation_axis,
+        angle);
+
+    const auto original_result = fesa::compute_beam3d2(original);
+    const auto rotated_result = fesa::compute_beam3d2(rotated);
+    ASSERT_TRUE(original_result.contribution.has_value());
+    ASSERT_TRUE(rotated_result.contribution.has_value());
+
+    const std::array<double, 12> original_displacement{
+        0.1, -0.3, 0.2, 0.04, -0.02, 0.03,
+        -0.2, 0.5, -0.1, -0.01, 0.06, -0.05,
+    };
+    const std::array<double, 12> rotated_displacement =
+        rotate_dofs(original_displacement, rotation_axis, angle);
+    const double original_energy = strain_energy(
+        original_result.contribution->global_stiffness,
+        original_displacement);
+    const double rotated_energy = strain_energy(
+        rotated_result.contribution->global_stiffness,
+        rotated_displacement);
+
+    expect_relative_near(
+        rotated_energy,
+        original_energy,
+        4096.0 * std::numeric_limits<double>::epsilon());
+}
+
+}  // namespace

tests/unit/fem/beam_frame_test.cpp

@@ -0,0 +1,278 @@
+#include <fesa/fem/beam_frame.hpp>
+
+#include <algorithm>
+#include <array>
+#include <cmath>
+#include <limits>
+#include <string_view>
+
+#include <gtest/gtest.h>
+
+namespace {
+
+constexpr double kUnitTolerance =
+    128.0 * std::numeric_limits<double>::epsilon();
+
+double dot(const fesa::Vec3 first, const fesa::Vec3 second) {
+    return first.x * second.x + first.y * second.y + first.z * second.z;
+}
+
+fesa::Vec3 cross(const fesa::Vec3 first, const fesa::Vec3 second) {
+    return {
+        first.y * second.z - first.z * second.y,
+        first.z * second.x - first.x * second.z,
+        first.x * second.y - first.y * second.x,
+    };
+}
+
+fesa::Vec3 rotate_about_axis(
+    const fesa::Vec3 value,
+    const fesa::Vec3 unit_axis,
+    const double angle) {
+    const double cosine = std::cos(angle);
+    const double sine = std::sin(angle);
+    const fesa::Vec3 axis_cross_value = cross(unit_axis, value);
+    const double axis_projection = dot(unit_axis, value);
+    return {
+        value.x * cosine + axis_cross_value.x * sine +
+            unit_axis.x * axis_projection * (1.0 - cosine),
+        value.y * cosine + axis_cross_value.y * sine +
+            unit_axis.y * axis_projection * (1.0 - cosine),
+        value.z * cosine + axis_cross_value.z * sine +
+            unit_axis.z * axis_projection * (1.0 - cosine),
+    };
+}
+
+void expect_vec_near(
+    const fesa::Vec3 actual,
+    const fesa::Vec3 expected,
+    const double tolerance = kUnitTolerance) {
+    EXPECT_NEAR(actual.x, expected.x, tolerance);
+    EXPECT_NEAR(actual.y, expected.y, tolerance);
+    EXPECT_NEAR(actual.z, expected.z, tolerance);
+}
+
+void expect_error(
+    const fesa::BeamFrameResult& result,
+    const std::string_view code) {
+    EXPECT_FALSE(result.frame.has_value());
+    ASSERT_EQ(result.diagnostics.size(), 1);
+    EXPECT_EQ(result.diagnostics[0].stage, fesa::DiagnosticStage::model);
+    EXPECT_EQ(result.diagnostics[0].severity, fesa::Severity::error);
+    EXPECT_EQ(result.diagnostics[0].code, code);
+}
+
+std::array<double, 12> multiply(
+    const fesa::Matrix12& matrix,
+    const std::array<double, 12>& vector) {
+    std::array<double, 12> product{};
+    for (std::size_t row = 0; row < product.size(); ++row) {
+        for (std::size_t column = 0; column < vector.size(); ++column) {
+            product[row] += matrix[row][column] * vector[column];
+        }
+    }
+    return product;
+}
+
+double squared_norm(const std::array<double, 12>& vector) {
+    double result = 0.0;
+    for (const double value : vector) {
+        result += value * value;
+    }
+    return result;
+}
+
+TEST(BeamFrame, AlignsLocalAxesForGlobalXAxisBeam) {
+    const auto result = fesa::make_beam_frame(
+        {0.0, 0.0, 0.0},
+        {2.0, 0.0, 0.0},
+        {0.0, 3.0, 0.0});
+
+    ASSERT_TRUE(result.frame.has_value());
+    EXPECT_TRUE(result.diagnostics.empty());
+    expect_vec_near(result.frame->ex, {1.0, 0.0, 0.0});
+    expect_vec_near(result.frame->ey, {0.0, 1.0, 0.0});
+    expect_vec_near(result.frame->ez, {0.0, 0.0, 1.0});
+}
+
+TEST(BeamFrame, ProducesRightHandedOrthonormalBasisForArbitraryAxis) {
+    const auto result = fesa::make_beam_frame(
+        {1.0, 2.0, 3.0},
+        {3.0, 4.0, 4.0},
+        {0.0, 1.0, 0.0});
+
+    ASSERT_TRUE(result.frame.has_value());
+    const fesa::BeamFrame& frame = *result.frame;
+    const double inverse_sqrt_five = 1.0 / std::sqrt(5.0);
+    expect_vec_near(frame.ex, {2.0 / 3.0, 2.0 / 3.0, 1.0 / 3.0});
+    expect_vec_near(
+        frame.ey,
+        {
+            -4.0 * inverse_sqrt_five / 3.0,
+            5.0 * inverse_sqrt_five / 3.0,
+            -2.0 * inverse_sqrt_five / 3.0,
+        });
+    expect_vec_near(
+        frame.ez,
+        {-inverse_sqrt_five, 0.0, 2.0 * inverse_sqrt_five});
+
+    EXPECT_NEAR(dot(frame.ex, frame.ex), 1.0, kUnitTolerance);
+    EXPECT_NEAR(dot(frame.ey, frame.ey), 1.0, kUnitTolerance);
+    EXPECT_NEAR(dot(frame.ez, frame.ez), 1.0, kUnitTolerance);
+    EXPECT_NEAR(dot(frame.ex, frame.ey), 0.0, kUnitTolerance);
+    EXPECT_NEAR(dot(frame.ey, frame.ez), 0.0, kUnitTolerance);
+    EXPECT_NEAR(dot(frame.ez, frame.ex), 0.0, kUnitTolerance);
+    EXPECT_NEAR(dot(frame.ex, cross(frame.ey, frame.ez)), 1.0, kUnitTolerance);
+}
+
+TEST(BeamFrame, RotatesCovariantlyWithBeamGeometry) {
+    constexpr fesa::Vec3 first{1.0, -2.0, 0.5};
+    constexpr fesa::Vec3 second{4.0, 1.0, 2.5};
+    constexpr fesa::Vec3 orientation{-1.0, 2.0, 3.0};
+    const double inverse_sqrt_fourteen = 1.0 / std::sqrt(14.0);
+    const fesa::Vec3 rotation_axis{
+        inverse_sqrt_fourteen,
+        2.0 * inverse_sqrt_fourteen,
+        3.0 * inverse_sqrt_fourteen,
+    };
+    constexpr double angle = 0.73;
+
+    const auto original =
+        fesa::make_beam_frame(first, second, orientation);
+    const auto rotated = fesa::make_beam_frame(
+        rotate_about_axis(first, rotation_axis, angle),
+        rotate_about_axis(second, rotation_axis, angle),
+        rotate_about_axis(orientation, rotation_axis, angle));
+
+    ASSERT_TRUE(original.frame.has_value());
+    ASSERT_TRUE(rotated.frame.has_value());
+    expect_vec_near(
+        rotated.frame->ex,
+        rotate_about_axis(original.frame->ex, rotation_axis, angle),
+        256.0 * std::numeric_limits<double>::epsilon());
+    expect_vec_near(
+        rotated.frame->ey,
+        rotate_about_axis(original.frame->ey, rotation_axis, angle),
+        256.0 * std::numeric_limits<double>::epsilon());
+    expect_vec_near(
+        rotated.frame->ez,
+        rotate_about_axis(original.frame->ez, rotation_axis, angle),
+        256.0 * std::numeric_limits<double>::epsilon());
+}
+
+TEST(BeamFrame, IsInvariantToFiniteOrientationScale) {
+    const double maximum = std::numeric_limits<double>::max();
+    const auto unit_scale = fesa::make_beam_frame(
+        {0.0, 0.0, 0.0},
+        {1.0, 1.0, 1.0},
+        {1.0, -1.0, 0.0});
+    const auto maximum_scale = fesa::make_beam_frame(
+        {0.0, 0.0, 0.0},
+        {1.0, 1.0, 1.0},
+        {maximum, -maximum, 0.0});
+
+    ASSERT_TRUE(unit_scale.frame.has_value());
+    ASSERT_TRUE(maximum_scale.frame.has_value());
+    expect_vec_near(maximum_scale.frame->ex, unit_scale.frame->ex);
+    expect_vec_near(maximum_scale.frame->ey, unit_scale.frame->ey);
+    expect_vec_near(maximum_scale.frame->ez, unit_scale.frame->ez);
+}
+
+TEST(BeamFrame, RejectsZeroLengthAxis) {
+    expect_error(
+        fesa::make_beam_frame(
+            {1.0, 2.0, 3.0},
+            {1.0, 2.0, 3.0},
+            {0.0, 1.0, 0.0}),
+        "model.zero_length_element");
+}
+
+TEST(BeamFrame, RejectsZeroOrientation) {
+    expect_error(
+        fesa::make_beam_frame(
+            {0.0, 0.0, 0.0},
+            {1.0, 0.0, 0.0},
+            {0.0, 0.0, 0.0}),
+        "model.invalid_orientation");
+}
+
+TEST(BeamFrame, RejectsOrientationParallelToElementAxis) {
+    expect_error(
+        fesa::make_beam_frame(
+            {0.0, 0.0, 0.0},
+            {2.0, 0.0, 0.0},
+            {5.0, 0.0, 0.0}),
+        "model.invalid_orientation");
+}
+
+TEST(BeamFrame, RejectsNumericallyParallelOrientationAtLargeScale) {
+    constexpr double transverse_ratio =
+        32.0 * std::numeric_limits<double>::epsilon();
+    expect_error(
+        fesa::make_beam_frame(
+            {0.0, 0.0, 0.0},
+            {2.0, 0.0, 0.0},
+            {1.0e200, transverse_ratio * 1.0e200, 0.0}),
+        "model.invalid_orientation");
+}
+
+TEST(BeamTransformation, IsIdentityForGlobalAlignedFrame) {
+    const fesa::Matrix12 transformation = fesa::beam_transformation({
+        {1.0, 0.0, 0.0},
+        {0.0, 1.0, 0.0},
+        {0.0, 0.0, 1.0},
+    });
+
+    for (std::size_t row = 0; row < transformation.size(); ++row) {
+        for (std::size_t column = 0;
+             column < transformation[row].size();
+             ++column) {
+            EXPECT_DOUBLE_EQ(
+                transformation[row][column],
+                row == column ? 1.0 : 0.0);
+        }
+    }
+}
+
+TEST(BeamTransformation, MapsGlobalComponentsIntoLocalDofOrder) {
+    const fesa::Matrix12 transformation = fesa::beam_transformation({
+        {0.0, 1.0, 0.0},
+        {0.0, 0.0, 1.0},
+        {1.0, 0.0, 0.0},
+    });
+    const std::array<double, 12> global{
+        1.0, 2.0, 3.0, 4.0, 5.0, 6.0,
+        7.0, 8.0, 9.0, 10.0, 11.0, 12.0,
+    };
+
+    EXPECT_EQ(
+        multiply(transformation, global),
+        (std::array<double, 12>{
+            2.0, 3.0, 1.0, 5.0, 6.0, 4.0,
+            8.0, 9.0, 7.0, 11.0, 12.0, 10.0,
+        }));
+}
+
+TEST(BeamTransformation, PreservesVectorEnergyForArbitraryFrame) {
+    const auto frame_result = fesa::make_beam_frame(
+        {-2.0, 1.0, 3.0},
+        {1.0, 5.0, 5.0},
+        {2.0, -1.0, 4.0});
+    ASSERT_TRUE(frame_result.frame.has_value());
+    const fesa::Matrix12 transformation =
+        fesa::beam_transformation(*frame_result.frame);
+    const std::array<double, 12> global{
+        3.0, -1.0, 2.0, 0.5, 4.0, -2.0,
+        -3.0, 6.0, 1.5, 2.5, -0.25, 5.0,
+    };
+
+    const std::array<double, 12> local =
+        multiply(transformation, global);
+    const double global_energy = squared_norm(global);
+    const double tolerance =
+        256.0 * std::numeric_limits<double>::epsilon() *
+        std::max(1.0, global_energy);
+    EXPECT_NEAR(squared_norm(local), global_energy, tolerance);
+}
+
+}  // namespace

tests/unit/fem/dof_manager_test.cpp

@@ -0,0 +1,199 @@
+#include <fesa/fem/dof_manager.hpp>
+
+#include <array>
+#include <cstdint>
+#include <stdexcept>
+#include <utility>
+#include <vector>
+
+#include <fesa/model/domain_builder.hpp>
+
+#include <gtest/gtest.h>
+
+namespace {
+
+fesa::Domain build_domain(
+    std::vector<fesa::Node> nodes,
+    fesa::StepDefinition step,
+    const std::array<fesa::NodeId, 2> element_nodes = {
+        fesa::NodeId{4},
+        fesa::NodeId{10},
+    }) {
+    fesa::DomainBuilder builder;
+    for (fesa::Node& node : nodes) {
+        builder.add_node(std::move(node));
+    }
+    builder.add_material({
+        fesa::MaterialId{0},
+        "Steel",
+        210.0e9,
+        0.3,
+    });
+    builder.add_section({
+        fesa::SectionId{0},
+        "General",
+        0.04,
+        1.2e-4,
+        1.4e-4,
+        2.0e-4,
+        0.03,
+        0.031,
+        fesa::ShearPropertySource::input,
+        fesa::Vec3{0.0, 1.0, 0.0},
+        {},
+    });
+    builder.add_beam_element({
+        fesa::ElementId{0},
+        fesa::EntityOrigin{"BeamPart", "Beam-1", 1},
+        element_nodes,
+        fesa::MaterialId{0},
+        fesa::SectionId{0},
+    });
+    builder.set_step(std::move(step));
+
+    auto result = std::move(builder).build();
+    EXPECT_TRUE(result.diagnostics.empty());
+    EXPECT_TRUE(result.domain.has_value());
+    return std::move(*result.domain);
+}
+
+std::vector<fesa::Node> nodes_in_external_label_order() {
+    return {
+        {
+            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},
+        },
+    };
+}
+
+TEST(EquationNumbering, UsesInternalNodeIdOrderNotExternalLabelOrStorageOrder) {
+    const fesa::Domain domain =
+        build_domain(nodes_in_external_label_order(), {"Load", {}, {}});
+
+    const fesa::DofManager dofs = fesa::DofManager::build(domain);
+
+    EXPECT_EQ(dofs.full_dof_count(), 18);
+    EXPECT_EQ(dofs.free_equation_count(), 18);
+    EXPECT_EQ(
+        dofs.equation({fesa::NodeId{4}, fesa::NodeDof::ux}),
+        std::optional<std::size_t>{0});
+    EXPECT_EQ(
+        dofs.equation({fesa::NodeId{4}, fesa::NodeDof::rz}),
+        std::optional<std::size_t>{5});
+    EXPECT_EQ(
+        dofs.equation({fesa::NodeId{10}, fesa::NodeDof::ux}),
+        std::optional<std::size_t>{6});
+    EXPECT_EQ(
+        dofs.equation({fesa::NodeId{20}, fesa::NodeDof::rz}),
+        std::optional<std::size_t>{17});
+}
+
+TEST(DofManager, NumbersOnlyFreeDofsAndReconstructsPrescribedValues) {
+    const fesa::Domain domain = build_domain(
+        nodes_in_external_label_order(),
+        {
+            "Load",
+            {
+                {fesa::NodeId{4}, 1, 1.25},
+                {fesa::NodeId{10}, 6, -0.5},
+            },
+            {},
+        });
+
+    const fesa::DofManager dofs = fesa::DofManager::build(domain);
+
+    EXPECT_EQ(dofs.full_dof_count(), 18);
+    EXPECT_EQ(dofs.free_equation_count(), 16);
+    EXPECT_EQ(
+        dofs.equation({fesa::NodeId{4}, fesa::NodeDof::ux}),
+        std::nullopt);
+    EXPECT_EQ(
+        dofs.equation({fesa::NodeId{4}, fesa::NodeDof::uy}),
+        std::optional<std::size_t>{0});
+    EXPECT_EQ(
+        dofs.equation({fesa::NodeId{10}, fesa::NodeDof::rx}),
+        std::optional<std::size_t>{8});
+    EXPECT_EQ(
+        dofs.equation({fesa::NodeId{10}, fesa::NodeDof::rz}),
+        std::nullopt);
+    EXPECT_EQ(
+        dofs.equation({fesa::NodeId{20}, fesa::NodeDof::ux}),
+        std::optional<std::size_t>{10});
+
+    std::vector<double> reduced(dofs.free_equation_count());
+    for (std::size_t equation = 0; equation < reduced.size(); ++equation) {
+        reduced[equation] = 100.0 + static_cast<double>(equation);
+    }
+
+    EXPECT_EQ(
+        dofs.reconstruct_full(reduced),
+        (std::vector<double>{
+            1.25,
+            100.0,
+            101.0,
+            102.0,
+            103.0,
+            104.0,
+            105.0,
+            106.0,
+            107.0,
+            108.0,
+            109.0,
+            -0.5,
+            110.0,
+            111.0,
+            112.0,
+            113.0,
+            114.0,
+            115.0,
+        }));
+}
+
+TEST(DofManager, MapsElementNodesToTwelveFullDofsInComponentOrder) {
+    const fesa::Domain domain = build_domain(
+        nodes_in_external_label_order(),
+        {"Load", {}, {}},
+        {fesa::NodeId{10}, fesa::NodeId{4}});
+    const fesa::DofManager dofs = fesa::DofManager::build(domain);
+
+    EXPECT_EQ(
+        dofs.element_full_dofs(domain.beam_elements().front()),
+        (std::array<std::size_t, 12>{
+            6, 7, 8, 9, 10, 11,
+            0, 1, 2, 3, 4, 5,
+        }));
+}
+
+TEST(DofManager, RejectsInvalidAddressesAndReducedVectorSize) {
+    const fesa::Domain domain =
+        build_domain(nodes_in_external_label_order(), {"Load", {}, {}});
+    const fesa::DofManager dofs = fesa::DofManager::build(domain);
+
+    EXPECT_THROW(
+        static_cast<void>(dofs.equation({
+            fesa::NodeId{4},
+            static_cast<fesa::NodeDof>(6),
+        })),
+        std::invalid_argument);
+    EXPECT_THROW(
+        static_cast<void>(
+            dofs.equation({fesa::NodeId{99}, fesa::NodeDof::ux})),
+        std::out_of_range);
+    EXPECT_THROW(
+        static_cast<void>(
+            dofs.reconstruct_full(std::vector<double>(17, 0.0))),
+        std::invalid_argument);
+}
+
+}  // namespace

tests/unit/fem/gauss_rule_test.cpp

@@ -0,0 +1,53 @@
+#include <fesa/fem/gauss_rule.hpp>
+
+#include <array>
+#include <limits>
+#include <stdexcept>
+
+#include <gtest/gtest.h>
+
+namespace {
+
+TEST(Quadrature, OnePointRuleIntegratesDegreeZeroAndOneExactly) {
+    const auto rule = fesa::gauss_rule_1d(1);
+
+    ASSERT_EQ(rule.size(), 1);
+    EXPECT_DOUBLE_EQ(rule[0].weight, 2.0);
+    EXPECT_DOUBLE_EQ(rule[0].weight * rule[0].xi, 0.0);
+}
+
+TEST(Quadrature, TwoPointRuleIntegratesThroughDegreeThreeWithinRoundoff) {
+    const auto rule = fesa::gauss_rule_1d(2);
+
+    ASSERT_EQ(rule.size(), 2);
+
+    double degree_zero = 0.0;
+    double degree_one = 0.0;
+    double degree_two = 0.0;
+    double degree_three = 0.0;
+    for (const auto& point : rule) {
+        const double xi_squared = point.xi * point.xi;
+        degree_zero += point.weight;
+        degree_one += point.weight * point.xi;
+        degree_two += point.weight * xi_squared;
+        degree_three += point.weight * xi_squared * point.xi;
+    }
+
+    // Sixteen ulps at unit scale cover only the rounding in these short sums.
+    constexpr double tolerance =
+        16.0 * std::numeric_limits<double>::epsilon();
+    EXPECT_NEAR(degree_zero, 2.0, tolerance);
+    EXPECT_NEAR(degree_one, 0.0, tolerance);
+    EXPECT_NEAR(degree_two, 2.0 / 3.0, tolerance);
+    EXPECT_NEAR(degree_three, 0.0, tolerance);
+}
+
+TEST(Quadrature, RejectsUnsupportedOrders) {
+    for (const int order : std::array{0, 3, -1}) {
+        EXPECT_THROW(
+            static_cast<void>(fesa::gauss_rule_1d(order)),
+            std::invalid_argument);
+    }
+}
+
+}  // namespace

tests/unit/fem/line2_shape_test.cpp

@@ -0,0 +1,71 @@
+#include <fesa/fem/line2_shape.hpp>
+
+#include <array>
+#include <limits>
+#include <stdexcept>
+
+#include <gtest/gtest.h>
+
+namespace {
+
+TEST(ShapeFunction, InterpolatesBothEndpoints) {
+    const auto first_endpoint = fesa::line2_shape(-1.0);
+    const auto second_endpoint = fesa::line2_shape(1.0);
+
+    EXPECT_DOUBLE_EQ(first_endpoint[0], 1.0);
+    EXPECT_DOUBLE_EQ(first_endpoint[1], 0.0);
+    EXPECT_DOUBLE_EQ(second_endpoint[0], 0.0);
+    EXPECT_DOUBLE_EQ(second_endpoint[1], 1.0);
+}
+
+TEST(ShapeFunction, FormsPartitionOfUnityWithinRoundoff) {
+    for (const double xi : std::array{-1.0, -0.25, 0.0, 0.4, 1.0}) {
+        const auto shape = fesa::line2_shape(xi);
+
+        // Four ulps at unit scale cover the two affine evaluations and sum.
+        constexpr double tolerance =
+            4.0 * std::numeric_limits<double>::epsilon();
+        EXPECT_NEAR(shape[0] + shape[1], 1.0, tolerance) << "xi=" << xi;
+    }
+}
+
+TEST(ShapeFunction, DerivativesSumToZeroAtEveryNaturalCoordinate) {
+    for (const double xi : std::array{-1.0, -0.25, 0.0, 0.4, 1.0}) {
+        const auto derivative = fesa::line2_shape_derivative(xi);
+
+        EXPECT_DOUBLE_EQ(derivative[0], -0.5);
+        EXPECT_DOUBLE_EQ(derivative[1], 0.5);
+        EXPECT_DOUBLE_EQ(derivative[0] + derivative[1], 0.0);
+    }
+}
+
+TEST(Jacobian, ReturnsHalfThePhysicalLength) {
+    EXPECT_DOUBLE_EQ(fesa::line2_jacobian(4.0), 2.0);
+    EXPECT_DOUBLE_EQ(fesa::line2_jacobian(0.25), 0.125);
+}
+
+TEST(Jacobian, DoesNotAdjustSmallPositiveFiniteLengths) {
+    constexpr double length = 1.0e-300;
+
+    EXPECT_DOUBLE_EQ(fesa::line2_jacobian(length), length / 2.0);
+}
+
+TEST(Jacobian, RejectsNonpositiveAndNonfiniteLengths) {
+    const double infinity = std::numeric_limits<double>::infinity();
+    const double nan = std::numeric_limits<double>::quiet_NaN();
+
+    for (const double length : std::array{0.0, -1.0, infinity, nan}) {
+        EXPECT_THROW(
+            static_cast<void>(fesa::line2_jacobian(length)),
+            std::invalid_argument);
+    }
+}
+
+TEST(Jacobian, RejectsLengthWhoseHalfUnderflowsToZero) {
+    EXPECT_THROW(
+        static_cast<void>(fesa::line2_jacobian(
+            std::numeric_limits<double>::denorm_min())),
+        std::invalid_argument);
+}
+
+}  // namespace