- Design principle: users never see CPU/GPU differences
- Three-layer architecture: User API / Backend Abstraction / Implementations
- Auto() backend selection based on hardware availability
- Physics{ElasticityPhysicsType} as single problem type
- Internal conversion between CPU arrays and GPU arrays
- solve!() with automatic dispatch to CPU or GPU backend
- 611 lines: Complete architecture design proposal
- get_basis_functions() and get_basis_derivatives() recommended
- Separation of concerns: topology and basis as separate parameters
- Side-by-side examples for all common element types
- Complete assembly example showing migration path
- Type-stable implementation with no performance penalty
- 302 lines: Comprehensive migration documentation
- Alternative to element-by-element assembly for GPU/matrix-free
- Node-by-node loop eliminates atomic operations on GPU
- Spider pattern: nodes couple with 10-30 neighbors not all N
- NodeToElementsMap: inverse connectivity (node → elements)
- get_node_spider() finds coupled nodes for sparse stiffness
- NodalStiffnessContribution: 3×3 blocks per node
- 307 lines: Experimental architecture with working prototype
- Compressible Neo-Hookean strain energy function
- Automatic differentiation for stress and tangent computation
- Dual constructor: Lamé (μ,λ) or engineering (E,ν)
- Total Lagrangian formulation with 2nd Piola-Kirchhoff stress
- Zero-allocation AD via Tensors.jl
- When to use: rubber, large deformation, contact mechanics
- 571 lines: Complete AD-based material model documentation
- Complete mathematical foundation of Hooke's law in tensor form
- Lamé parameters derived from Young's modulus and Poisson's ratio
- compute_stress() implementation achieving ~25 ns execution
- Fourth-order elasticity tensor with symmetries
- Zero-allocation SIMD-optimized implementation
- Physical constraints and thermodynamic admissibility
- 736 lines: Authoritative implementation documentation
- Reference implementation of element-by-element assembly
- ElementAssemblyData and ElementContribution data structures
- Sparse matrix assembly in COO then CSC format
- scatter_to_global!() adds local to global system
- Penalty method for Dirichlet BCs
- Matrix-vector product interface for GMRES
- 479 lines: Complete documentation with examples and tests
- Mathematical derivation of F = I + ∇u for finite strain
- Zero-allocation implementation achieving 34 ns median
- LLVM IR analysis confirms 0 heap allocations
- 92 SIMD vector operations detected
- Small strain vs finite strain formulations
- Comparison with old deprecated eval_dbasis!() API
- 502 lines: Complete performance analysis with benchmarks
Document decision to store integration point indices instead of data in Element struct.
Key rationale: Elements should store relationships (indices), not data, for memory
efficiency and consistency with node connectivity pattern. Aligns with nodal assembly
approach and GPU-friendly architecture.
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)
- Relocate docs/blog/immutability_performance.md to docs/src/book/blog/
- Comprehensive guide on immutable material models with Tensors.jl
- Covers LinearElastic, NeoHookean, PerfectPlasticity implementations
- Includes full benchmarks: 5× speedup for linear, 21× for plasticity
- Zero allocation performance validated
- Part of three-tier docs reorganization under standard docs/src/ structure
Design for handling both nodal and element fields in nodal assembly:
Architecture:
- Nodes have geometry (immutable)
- Elements have connectivity + fields (immutable struct)
- Nodal fields: displacement, temperature, contact pressure
- Element fields: integration point data (σ, ε_plastic, α, C)
Update pattern:
- Create new field containers (NamedTuples)
- Create new elements with updated fields
- Shallow copy element vector, replace elements
- All immutable (GPU-compatible, thread-safe)
GPU kernel:
- Loops over nodes (nodal assembly)
- Accesses nodal_fields for global quantities
- Accesses element.fields for integration point data
- Gathers from connected elements (node_to_elements)
- No atomic operations (each node owns DOFs)
Material state update:
- Process elements in parallel (Threads.@threads)
- Extract nodal displacements from solution
- Compute strains at integration points
- Run material model (plasticity, damage, etc.)
- Create new elements with updated state
- Return new problem with updated fields
Newton iteration:
- Residual uses element.fields.C (current tangent)
- GMRES with matrix-free matvec (nodal assembly)
- Material update after each iteration
- All data structures immutable throughout
Benchmarks show creating new containers ~1000× faster than deepcopy
Explains why JuliaFEM uses nodal assembly instead of element assembly:
Five major advantages:
1. No atomic operations on GPU (each node writes to own DOFs)
2. Contact mechanics is natural (forces at nodes, not elements)
3. Clean domain decomposition (explicit node ownership for MPI)
4. Better cache locality (sequential node processing)
5. Adaptive refinement easier (local node operations)
Key data structure:
- NodeSet contains nodes + elements + node_to_elements connectivity
- Inverse connectivity enables gathering from connected elements
- Fields accessed via node_set.fields (type-stable)
Algorithm:
- Loop over nodes (not elements)
- Each node gathers contributions from connected elements
- Direct write to owned DOFs (no race conditions)
- Perfect for matrix-free Krylov methods
Validated with demo:
- CPU/GPU results match exactly (0.0 relative error)
- Average 3.24 elements per node (efficient gathering)
- Natural integration with contact mechanics
Compares to traditional element assembly:
- Element: scatter to nodes (atomic ops, cache misses)
- Nodal: gather from elements (no atomics, better cache)
Complete algorithm for GPU-resident FEM solver with nodal assembly:
- Data partitioning by node ownership (domain decomposition)
- GPU-resident data structures (nodes, elements, connectivity, state)
- Three GPU kernels: residual, matvec, state update
- MPI communication patterns for interface nodes
- Full Newton-GMRES loop on GPU (data stays resident)
Architecture:
- Each GPU owns subset of nodes (exclusive ownership)
- Ghost elements copied for gathering during assembly
- node_to_elements connectivity enables nodal assembly
- No atomic operations (each GPU writes to owned DOFs only)
Key features:
- Data moves to GPU once at start, back once at end
- GMRES iterations entirely on GPU (Arnoldi steps)
- Material state updates on GPU (integration points)
- MPI exchanges only for interface DOFs between iterations
- O(N) memory per GPU (matrix-free)
Handles nonlinearity:
- Element state contains σ, ε_plastic, α, C (tangent)
- Residual kernel uses current stress/tangent
- State update kernel after convergence
- Natural for contact mechanics (nodal forces)
Status: Design document for future GPU implementation
Explains the correct pattern for matrix-vector products in Krylov methods:
- Fields accessed through element_set (not passed separately)
- GPU kernel computes y=K*x (not K itself)
- O(N) memory (vs O(N²) for stored matrix)
- Type-stable field access (compile-time types)
Key insights:
- GMRES needs matvec operation, not the matrix
- ElementSet contains elements + fields together
- Zero allocations with immutable connectivity/fields
- Natural pattern for contact mechanics (nodal updates)
- Material state separate from field parameters
Compares old vs new approach:
- Old: Dict{String,Any} in element (type-unstable)
- New: NamedTuple in ElementSet (type-stable)
- Old: O(N²) matrix storage
- New: O(N) matrix-free operator
Validated with gpu_elementset_matvec_demo.jl:
- GPU/CPU results match exactly
- Fields accessed naturally through element_set
- Returns y vector (what Krylov methods need)
Analyzes field storage patterns and recommends ElementSet approach:
- Element has NO field type parameter (simpler type)
- ElementSet groups elements + shared fields
- Fields can be NamedTuple, struct, any type-stable container
- Embraces immutability (GPU-compatible, thread-safe)
- Separates mutable state from immutable parameters
Design rationale:
- Benchmarks show NamedTuple gives 9-92× speedup vs Dict
- Immutability enables GPU execution without copying
- Creating new containers ~1000× faster than deepcopy
- Matches physical reality (material properties per set)
Compares three options:
1. Fields as type parameter (type proliferation)
2. ElementSet pattern (RECOMMENDED)
3. Hybrid approach (too complex)
Addresses common concerns:
- Time-dependent fields (use interpolation)
- Material state (separate mutable arrays)
- Custom field types (any type-stable container works)
Status: Ready for implementation
Create 299-line focused design rationale explaining why type stability is essential
for CPU/GPU/MPI performance, without mandating specific storage patterns.
Executive summary (lines 16-26):
- v0.5.1 Dict{String,Any}: 9-92× performance degradation
- Type-unstable code CANNOT run on GPUs
- Significant MPI communication overhead
- Document does NOT prescribe storage location
- Demonstrates why type stability at access points is essential
- Key: Storage pattern matters less than type inference
Problem analysis (lines 28-67):
- Type instability definition: Runtime dispatch when type unknown at compile time
- Why it matters: 10-100× slower CPU, GPU compilation fails, MPI serialization
- Measured impact table: 9-92× speedup, 0 allocations with type stability
- Critical: Zero allocations required for GPU kernels
Design requirements (lines 69-136):
1. Type stability at access points (compiler must infer types)
- Fields could be element-local, global arrays, or arguments
- Access pattern must be type-stable regardless
2. Zero allocations in hot paths (GPU/MPI requirement)
- Assembly loop must allocate nothing
3. Contiguous memory layout (GPU/MPI optimization)
- CUDA transfers contiguous arrays directly
4. Immutable where possible (safe parallelism)
- Thread-safe reads without locks
Demonstrated solutions (lines 138-201) - EXAMPLES, not mandates:
1. NamedTuple container: Simple, type-stable, immutable
2. Struct with typed fields: Explicit, self-documenting
3. Passed as arguments: Maximum type stability, explicit dependencies
- All three achieve type stability
- Choice depends on use case, not performance
GPU and MPI rationale (lines 203-237):
- GPU execution: CUDA requires all code type-stable
- Mock demonstration in benchmarks/gpu_mpi_mock.jl
- MPI communication: Typed arrays use fast memcpy vs slow serialization
- Type stability enables identical code for CPU/GPU
Recommendations (lines 239-256):
- Use type-stable access patterns (REQUIRED)
- Prefer immutable data structures (threading/GPU)
- Pre-allocate caches (zero allocations)
- Use contiguous arrays (GPU/MPI transfer)
- Profile with @btime (verify zero allocations)
- Does NOT mandate: Storage location, container type, dynamic vs static
Validation (lines 258-275):
- benchmarks/field_storage_comparison.jl: 9-92× CPU speedup
- benchmarks/gpu_mpi_mock.jl: GPU/MPI patterns
- benchmarks/VALIDATION_RESULTS.md: Summary table
Conclusion (lines 277-299):
- Type stability is fundamental requirement, not implementation detail
- Enables: High CPU performance, GPU execution, efficient MPI, safe threading
- v1.0 must ensure type stability at access points
- Storage pattern is secondary concern (memory, cache, API)
- Next steps: Review, benchmark, choose pattern, implement, validate CUDA
Key difference from v1: Shorter (299 vs 1114 lines), focused on WHY not HOW,
explicitly states storage pattern is flexible, emphasizes GPU/MPI requirements.
Platform: Julia 1.12.1, November 9, 2025
Series: The JuliaFEM Book, Chapter 5
Status: Design rationale with validated measurements
Major documentation update (380 additions, 167 deletions):
- Explained topology is pure geometry (NO hardcoded node counts)
- Clarified basis determines BOTH polynomial degree AND node count
- Distinguished node count (connectivity) vs DOF count (unknowns)
- Added examples: Nedelec (edge DOFs), Raviart-Thomas (face DOFs)
- Documented Lagrange{Topology, P} parametric architecture
- Showed why Tri3/Quad4/Tet10 names are anti-pattern
- Updated all code examples to use new architecture
- Explained Serendipity vs full Lagrange tensor products
- Added performance implications and trade-offs
- Showed how one Triangle topology works for P1/P2/P3/Nedelec/etc
**Purpose:** Comprehensive justification for all technical decisions prioritizing
performance over convenience.
**Key Principles:**
- Efficiency > Educativeness (when forced to choose)
- Type stability over everything (100× performance difference)
- No free lunch - Julia doesn't make miracles
- HPC requires discipline and trade-offs
**Core Decisions Justified:**
1. **No Dynamic Field System**
- field["foo"] = x is 100× slower (Dict{String,Any})
- Type-stable structs only
- Sacrifice: Runtime flexibility
- Gain: Performance
2. **Immutable Data Structures**
- struct over mutable struct
- Sacrifice: Convenient mutation
- Gain: 2-10× speedup, thread-safety, stack allocation
3. **NTuple Over Vector**
- Compile-time size → SIMD optimization
- Sacrifice: Dynamic sizing
- Gain: Zero allocations, type stability
4. **Monolithic Over Multi-Package**
- Learned from 2015-2019 mistake
- Sacrifice: Small dependencies
- Gain: It actually works
5. **Manual Derivatives (hot paths)**
- 30× faster than AD for Tet10
- Sacrifice: More code
- Gain: Assembly loops stay fast
6. **Matrix-Free Methods**
- Design for 1M+ DOF from day 1
- Cannot retrofit later
7. **Explicit Over Implicit**
- No magic, show the steps
- Debuggable and teachable
**Hierarchy of Values:**
1. Correctness
2. Performance
3. Maintainability
4. Educativeness
5. Convenience
**What We're Giving Up:**
- Runtime flexibility (no element["custom_field"])
- Dynamic problem definition (no runtime topology changes)
- Duck typing convenience
- Small dependencies
- Beginner-friendly magic
**What We're Getting:**
- 10× single-thread speedup target
- 1M DOF contact problems
- Thread/GPU/distributed scalability
- Real HPC capability
**The Hard Truth:**
From Issue #266: "Do like Python, be slow like Python. Know what you do
before compiling, and be fast like C. There's no free lunch."
**Success Metrics:**
- ✅ Zero allocations in assembly
- ✅ Type-stable hot paths
- 🎯 10× faster than v0.5.1
- 🎯 1M DOF in < 1 hour
- 🎯 100+ thread scaling
**Use Cases:**
- "Why can't I use Dict?" → Point here
- "Why immutable?" → Point here
- "Why manual derivatives?" → Point here
- Any "why not convenience?" → Point here
**Status:** Living document, updated as we learn
See: Issue #266, TECHNICAL_VISION.md, benchmark results
**Three Manuals for Three Audiences:**
1. **User Manual** (docs/user/) - "Just Get It Done"
- For end users, engineers, students
- Simple, practical, step-by-step
- Quick start, tutorials, examples, troubleshooting
- Philosophy: Show me how to solve my problem
2. **Contributor Manual** (docs/contributor/) - "Show Me the Code"
- For developers, contributors, advanced users
- Technical, detailed, design rationale
- Testing, architecture, performance, CI/CD
- Philosophy: Explain HOW and WHY
3. **The JuliaFEM Book** (docs/book/) - "Let Me Show You How I Think"
- For researchers, theory nerds, and Jukka
- Comprehensive, educational, opinionated, personal
- Math foundations, design philosophy, history, research
- Philosophy: Mix theory, code, and personal experience
**Reorganization:**
- Moved: TESTING_PHILOSOPHY.md → contributor/testing_philosophy.md
- Moved: STATUS.md → contributor/status.md
- Moved: TEST_FIXES_NEEDED.md → contributor/test_fixes_needed.md
- Moved: lagrange_basis_functions.md → book/lagrange_basis_functions.md
- Moved: benchmarks/ → book/benchmarks/
- Created: docs/README.md (main index explaining structure)
- Created: README.md in each section explaining audience and contents
- Updated: All references in scripts and source files
**Naming:** All docs now lowercase (testing_philosophy not TESTING_PHILOSOPHY)
**Benefits:**
- Clear separation of concerns
- Users don't get overwhelmed with implementation details
- Contributors get technical depth
- Book preserves deep theory and personal insights
- Each manual optimized for its audience
**Next:** Populate each section with appropriate content
Jukka Aho