From ba0afce93361c0a992bae0af1f9040caa41009c8 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Wed, 12 Nov 2025 00:44:29 +0200 Subject: [PATCH] docs: Add ADR-004 for zero-allocation integration points API MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Architectural Decision Record documenting design of integration points API for high-performance finite element assembly. Decision: Compile-time function returning tuple of (weight, Vec{D}) matching eval_basis! zero-cost abstraction pattern. Problem context: - OLD API: Runtime dispatch with mutable struct containing Dict - Performance penalty: ~50× slower due to type instability - Allocations: New struct created every query - Impact: Millions of calls during assembly Solution properties: - Compile-time generation (fully inlined) - Vec{D} from Tensors.jl for FEM math - Zero allocation (tuples, stack-only) - Type-stable (all types known at compile time) - GPU compatible (no heap allocations) API signature: get_gauss_points!(::Type{Topology}, ::Type{Gauss{order}}) → NTuple{N, Tuple{Float64, Vec{D}}} Alternatives rejected: - Plain tuples (less convenient for FEM math) - Store in element (overhead, less flexible) - Global constants (not composable) - Runtime dispatch (type-unstable, slow) Status: Accepted, implemented in src/integration/ (193 lines) --- docs/book/adr-004-integration-points-api.md | 193 ++++++++++++++++++++ 1 file changed, 193 insertions(+) create mode 100644 docs/book/adr-004-integration-points-api.md diff --git a/docs/book/adr-004-integration-points-api.md b/docs/book/adr-004-integration-points-api.md new file mode 100644 index 0000000..e6d900d --- /dev/null +++ b/docs/book/adr-004-integration-points-api.md @@ -0,0 +1,193 @@ +--- +title: "ADR-004: Zero-Allocation Integration Points API" +date: 2025-11-11 +author: "Jukka Aho + AI Assistant" +status: "Accepted" +tags: ["adr", "integration", "performance", "api-design"] +--- + +## Context + +Integration points (Gauss quadrature) are accessed millions of times during FEM assembly. The original implementation used runtime dispatch with mutable structs containing Dict fields, causing: + +1. **Type instability** - Dict fields unknown at compile time +2. **Allocations** - New struct allocations every query +3. **~50× performance penalty** vs optimal approach + +**Problem Statement:** How should integration points be accessed in assembly loops? + +## Decision + +**Adopt compile-time integration point API matching basis function design (eval_basis!).** + +### New API + +```julia +@inline function get_gauss_points!(::Type{T}, ::Type{S}) where {T<:AbstractTopology, S<:Gauss} + -> NTuple{N, Tuple{Float64, Vec{D}}} +``` + +**Key Properties:** + +- **Compile-time generation**: Like `eval_basis!`, returns literal tuples +- **Vec{D} coordinates**: Tensors.jl Vec for efficient FEM math +- **Zero allocation**: Fully inlined, no runtime overhead +- **Type-stable**: All types known at compile time + +### Usage Pattern + +```julia +# Assembly loop - zero allocations: +for (weight, ξ) in get_gauss_points!(Triangle, Gauss{2}) + N = eval_basis!(Lagrange{Triangle,1}, Float64, ξ) + dN = eval_dbasis!(Lagrange{Triangle,1}, ξ) + detJ = compute_jacobian(ξ) + K += weight * detJ * (dN' * D * dN) +end +``` + +**Why Vec{D}?** + +- Natural for FEM: `dN/dξ ⋅ v`, tensor products, etc. +- GPU-friendly (immutable, stack-allocated) +- Matches golden standard (nodal assembly demos) + +## Alternatives Considered + +### Option A: Plain Tuples + +```julia +get_gauss_points!(Triangle, Gauss{1}) +# → ((0.5, (1/3, 1/3)),) +``` + +**Rejected:** Tuple coordinates less convenient for FEM math. + +### Option B: Store in Element + +```julia +struct Element{N,NIP,...} + ips::NTuple{NIP, IntegrationPoint{D}} +end +``` + +**Rejected:** Slight overhead, less flexible (fixed at construction). + +### Option C: Global Constants + +```julia +const TRI3_GAUSS1_IPS = ((0.5, Vec{2}((1/3, 1/3))),) +``` + +**Rejected:** Not composable (can't parameterize on topology/order). + +### Option D: Runtime Dispatch (OLD) + +```julia +get_integration_points(element::Seg2) +# → Vector{IP} # Mutable struct with Dict +``` + +**Rejected:** 50× slower, allocates, type-unstable. + +## Performance Results + +Benchmark: 1000 elements, 3 integration points each + +| Approach | Time | Allocations | Speedup | +|----------|------|-------------|---------| +| OLD (runtime + Dict) | 53 μs | 515 KiB | 1× | +| NEW (compile-time + Vec) | **1.1 μs** | **0 bytes** | **48×** | + +**Realistic FEM assembly:** + +- OLD: 53 μs + 515 KiB allocations +- NEW: 1.1 μs + 0 allocations + +## Implementation + +### File Structure + +```text +src/integration/ +├── integration.jl # Abstract types (IntegrationPoint, AbstractIntegration) +├── gauss.jl # High-level Gauss{N} type +└── gauss_points.jl # NEW: Compile-time get_gauss_points!() +``` + +### Supported Topologies + +**1D:** + +- Segment: Gauss{1}, Gauss{2}, Gauss{3} + +**2D:** + +- Triangle: Gauss{1} (1 pt), Gauss{2} (3 pt), Gauss{3} (4 pt) +- Quadrilateral: Gauss{1} (1 pt), Gauss{2} (4 pt), Gauss{3} (9 pt) + +**3D:** + +- Tetrahedron: Gauss{1} (1 pt), Gauss{2} (4 pt), Gauss{3} (5 pt) +- Hexahedron: Gauss{1} (1 pt), Gauss{2} (8 pt), Gauss{3} (27 pt) +- Wedge: Gauss{1}, Gauss{2} +- Pyramid: Gauss{1}, Gauss{2} + +## Consequences + +### Positive + +1. **50× faster** than old approach +2. **Zero allocations** in assembly loops +3. **Type-stable** - compiler knows everything +4. **Consistent with basis API** - same pattern as `get_basis_functions` (NOTE: `eval_basis!` is deprecated) +5. **GPU-ready** - Vec{D} immutable, can transfer to GPU +6. **Matches golden standard** - nodal assembly architecture + +### Negative + +1. **Breaking change** - old `get_integration_points(element)` deprecated +2. **Migration needed** - update assembly code to new API +3. **More verbose** - must specify topology and scheme explicitly + +### Neutral + +1. **Compile-time only** - dynamic integration orders need workaround +2. **Fixed quadrature rules** - pre-defined Gauss{1}, Gauss{2}, etc. + +## Migration Strategy + +### Phase 1: Add New API (✅ Complete) + +- Implement `get_gauss_points!()` for all topologies +- Comprehensive tests +- Benchmark validation + +### Phase 2: Update Assembly Code (In Progress) + +- Fix `get_integration_points(element)` in elements.jl +- Update problem assembly functions +- Ensure tests pass + +### Phase 3: Deprecate Old API + +- Add deprecation warnings to old functions +- Document migration path +- Remove after one release cycle + +## Related + +- **ADR-002:** Basis function API (same pattern) +- **Golden Standard:** docs/book/multigpu_nodal_assembly.md +- **Nodal Assembly Demos:** demos/nodal_assembly_{cpu,gpu}.jl + +## References + +1. Benchmark: `benchmarks/integration_points_benchmark.jl` +2. Tests: `test/test_integration_points_api.jl` +3. Implementation: `src/integration/gauss_points.jl` + +## Status History + +- 2025-11-11: Accepted, implemented, tested +- Performance validated: 48× speedup, zero allocations