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)
GPU port of nodal assembly strategy with CUDA kernels demonstrating
atomic-free assembly on GPU using node-parallel approach.
GPU kernel design:
- One thread per node (not per element)
- Each thread gathers from touching elements
- No atomic operations (node ownership)
- Coalesced memory access via node ordering
Kernel structure:
- Thread ID maps to node ID
- Loop over elements touching this node
- Loop over element nodes for block contributions
- Compute 3×3 stiffness blocks with Tensors.jl
- Accumulate locally, write once to global
Data layout:
- node_to_elements: CSR-like structure on GPU
- Element data: Array of Structs (immutable elements)
- Node displacement: Flat vector (3*n_nodes)
- Result: Flat vector (3*n_nodes)
Performance characteristics:
- Memory bandwidth bound (not compute bound)
- Benefits from coalescing (sequential node access)
- Scalable to multi-GPU (domain decomposition)
- No synchronization within kernel
Comparison to element assembly:
- Element: N_elem threads, atomic scatter
- Nodal: N_nodes threads, no atomics
Reference: CPU version in nodal_assembly_cpu.jl
Complete CPU reference implementation of Newton-Krylov solver with
Anderson acceleration for nonlinear elasticity with plasticity.
Solver components:
- Newton outer loop (nonlinear iterations)
- GMRES inner loop (linear solve, matrix-free)
- Anderson acceleration (convergence improvement)
- Adaptive GMRES tolerance (Eisenstat-Walker formula)
Matrix-free strategy:
- No tangent matrix assembly
- Jacobian-vector product via finite differences: J·v ≈ [r(u+ε·v)-r(u)]/ε
- Residual assembly: r(u) = f_int(u) - f_ext
- Each GMRES iteration = 2 residual evaluations
Plasticity handling:
- Radial return mapping at each Gauss point
- Material state tracking (ε_p, α) during iterations
- State update only on Newton convergence
- Von Mises yield criterion with perfect plasticity
Reference for GPU implementation:
- Validates numerical correctness
- Establishes performance baseline
- Documents algorithm flow for GPU port
- Shows data dependencies and kernel opportunities
Problem: 3D elasticity with J2 plasticity, Tet4 mesh
GPU implementation for 10-node tetrahedral elements demonstrating
higher-order finite elements with quadratic shape functions.
Tet10 specifics:
- 10 nodes per element (vertices + edge midpoints)
- 4-point Gauss quadrature (order 2)
- Quadratic shape functions (N_i second-order polynomials)
- Shape function derivatives via analytical formulas
Challenges vs Tet4:
- More integration points (4 vs 1)
- More DOFs per element (30 vs 12)
- More complex shape functions
- Larger local stiffness (10×10 vs 4×4 blocks)
GPU kernel modifications:
- Loop over 4 Gauss points instead of 1
- Evaluate quadratic shape functions at each IP
- Accumulate contributions from all IPs
- Scatter 30 DOFs per element (not 12)
Benefits of Tet10:
- Better stress/strain representation
- Fewer elements needed for accuracy
- Curved boundary representation
- Higher convergence rate
Same problem setup: 3D cantilever with steel properties
Test validates GPU higher-order element implementation (430 lines).
First working GPU assembly implementation (proof-of-concept stage)
demonstrating complete FEM solve staying on GPU for 2D elasticity.
Implementation:
- Element-parallel GPU kernel for Quad4 elements
- 2×2 Gauss quadrature on GPU
- Plain vector approach (before Tensors.jl integration)
- Matrix-free Jacobian-vector product
- Complete Newton-Krylov loop on GPU
- BC enforcement via masking
Test case: 10×10 Quad4 mesh (100 elements, 242 DOFs)
- Material: Steel (E=200 GPa, ν=0.3)
- BC: Fixed left edge, displacement on right edge
Architecture validation:
- GPU assembly matches CPU (error < 1e-15)
- Entire solve stays on GPU (no ping-pong)
- Only transfers: mesh (once) + u0/u_final (boundary)
Note: This is the initial version using plain vectors and manual
indexing. See gpu_assembly_poc_tensors.jl for corrected version
using proper Tensors.jl material API (606 lines).
Complete GPU-accelerated FEM solve for 3D cantilever beam using
nodal assembly strategy and matrix-free Newton-Krylov solver.
GPU implementation:
- Gmsh mesh generation (same as CPU version)
- Data transfer to GPU (nodes, connectivity, BC)
- GPU kernels for nodal assembly (element contributions)
- Matrix-free Jacobian-vector product on GPU
- GMRES solver on GPU (Krylov.jl with CuArrays)
- CPU fallback for Anderson acceleration
Problem characteristics:
- Geometry: 10m × 1m × 1m cantilever beam
- Elements: Tet4 from Gmsh
- Material: Steel (E=210 GPa, ν=0.3)
- BC: Fixed left end, tip force on right end
Architecture:
- Single GPU transfer: mesh + BC → GPU at start
- Entire Newton-Krylov loop stays on GPU
- Single result transfer: u_final ← GPU at end
- No ping-pong between CPU and GPU during solve
Demonstrates complete GPU FEM pipeline from meshing to solution
with realistic geometry and material properties.
Compares traditional element assembly vs nodal assembly on CPU for
cantilever beam example, validating assembly equivalence and measuring
performance characteristics.
Comparison:
- Element assembly: Traditional FEM (loop over elements, atomic scatter)
- Nodal assembly: Modern approach (loop over nodes, block operations)
Validation:
- Residual equivalence (element vs nodal assembly)
- Stiffness operator equivalence (matvec comparison)
- Assembly time comparison
- Memory allocation tracking
Problem: Same cantilever beam as cantilever_beam_simple.jl
- Tet4 mesh from Gmsh
- Steel properties
- Fixed left, force on right
Demonstrates CPU assembly strategies before GPU implementation,
establishing baseline for GPU performance comparison.
Demonstrates modern Physics API for 3D elasticity on realistic geometry
using Gmsh mesh generation and both direct/iterative solvers.
Features:
- Gmsh mesh generation (10m × 1m × 1m cantilever beam)
- Tet4 elements with controlled mesh size (lc=1.5)
- Physics API setup (Elasticity, continuum formulation)
- Steel material properties (E=210 GPa, ν=0.3)
- Boundary conditions: Fixed left end, force on right end
Problem setup:
- Geometry: Cantilever beam (aspect ratio 10:1:1)
- Discretization: Tet4 elements from Gmsh
- Loading: Tip force applied via Neumann BC
- Constraints: Fixed end via Dirichlet BC
Workflow demonstration:
1. Mesh generation with Gmsh API
2. Physics problem creation
3. Solver setup (direct or iterative)
4. Post-processing and visualization
Educational example showing complete FEM workflow from meshing
to solution with modern JuliaFEM API (183 lines).
Demonstrates modern Physics API for solving elasticity problems using
CPU backend with element assembly.
Features:
- Simple 2-element beam mesh (Hex8 elements, 12 nodes, 36 DOFs)
- Immutable Element API with field-based material properties
- Physics problem setup (Elasticity, continuum formulation)
- Material properties: Steel (E=210 GPa, ν=0.3)
Demonstrates workflow:
1. Create mesh (nodes dictionary + connectivity tuples)
2. Create Physics problem (Elasticity with continuum formulation)
3. Build elements with immutable API (fields tuple)
4. Add elements to physics
Educational example showing modern API usage for elasticity
problems with clean separation between geometry and physics (131 lines).
Documents architectural correction from manual Voigt indexing to proper
Tensors.jl material modeling in GPU assembly proof-of-concept.
Problem identified:
- Initial POC used plain vectors instead of SymmetricTensor
- Hardcoded constitutive matrix instead of material API
- Manual index arithmetic for stress components
- Didn't match established material_modeling.md architecture
Solution implemented:
- SymmetricTensor{2,2} for 2D strain and stress
- Material API: compute_stress(material, ε)
- LinearElastic struct with Lamé parameters
- Clean tensor operations matching theory
- GPU compatible (Tensors.jl works on CUDA)
Key architectural changes:
1. Material model struct (LinearElastic with E, ν)
2. Material API with Hooke's law (σ = λ·tr(ε)·I + 2μ·ε)
3. SymmetricTensor strain computation (εxx, εyy, γxy/2)
4. Stress-to-force conversion (Bᵀ·σ operator)
Reference: demos/gpu_assembly_poc_tensors.jl (264 lines)
Implementation roadmap for GPU-accelerated nonlinear solver pipeline derived
from CPU reference implementation (newton_krylov_anderson_cpu.jl).
Breakdown of solver pipeline:
- Outer loop: Newton iterations (residual assembly + line search)
- Middle loop: GMRES iterations (matrix-free matvec + Arnoldi)
- Inner operation: Element residual assembly with plasticity
Key GPU kernels identified:
1. Element residual assembly (workhorse kernel, nodal scatter with atomics)
2. Vector operations (standard cuBLAS: axpy, dot, norm)
Four-phase implementation strategy:
1. Single kernel test (residual assembly CPU vs GPU)
2. Matrix-free matvec test (Jacobian-vector product validation)
3. GMRES on GPU (Krylov.jl with CuArrays)
4. Complete pipeline (GPU main loop, CPU Anderson acceleration)
Includes plastic state GPU representation (NTuple vs SymmetricTensor),
kernel launch configuration, and atomic scatter pattern (296 lines).
Shell script for quick validation of GPU benchmarks without running
full suites (which can take minutes on large problem sizes).
Features:
- Julia installation check (version validation)
- GPU availability detection (nvidia-smi query)
- CUDA.jl functionality verification
- Quick state management test (10K elements, not 1M)
- Quick matrix-free test (1K DOFs, not 1M)
- Both CPU and GPU paths tested
- Error handling with informative messages
Runs small problem sizes to verify:
- Code compiles and loads correctly
- CUDA kernels launch without errors
- Basic functionality works before long benchmarks
- Development workflow (test before full run)
Executable: chmod +x benchmarks/test_gpu_benchmarks.sh (159 lines)
Comprehensive benchmarking of J2 plasticity with radial return mapping:
Tests performed:
1. Single evaluation: elastic path (below yield) vs plastic path
2. Zero-allocation verification for both branches
3. Type stability validation
4. State management overhead (fresh vs history)
5. Hardening parameter sensitivity analysis
6. Assembly loop simulation (realistic FEM usage)
7. Comparison to LinearElastic and NeoHookean
8. Strain level scalability (elastic to plastic transition)
Validates:
- Plastic path overhead (radial return vs elastic)
- State handling performance (PlasticityState vs NoState)
- Zero allocations maintained even with mutable state
- Type stability for both converged and trial states
Includes von Mises stress computation, yield surface check, and
algorithmic tangent calculation (320 lines).
Detailed performance analysis of automatic differentiation overhead in
hyperelastic stress computation for Neo-Hookean material model.
Key analyses:
1. Single stress evaluation timing (typical FEM assembly use case)
2. Allocation verification (zero-allocation requirement)
3. Component breakdown (strain energy vs stress vs tangent)
4. Scaling with problem size (assembly loop performance)
Compares NeoHookean (AD) against LinearElastic (manual derivatives) to
quantify AD overhead in production FEM assembly loops.
Results inform whether AD is suitable for hot paths vs manual derivatives
for performance-critical material models (260 lines).
Comprehensive benchmark comparing three nonlinear solver strategies:
1. Traditional Newton (full Jacobian assembly + direct solve)
2. Matrix-free Newton-Krylov (GMRES, no Jacobian matrix)
3. Matrix-free with Anderson acceleration (accelerated convergence)
Problem: 3D nonlinear elasticity with cubic nonlinearity
- r(u) = K·u + β·(K·u)³ - f
- Jacobian-vector product via finite differences: J·v ≈ [r(u+ε·v) - r(u)]/ε
Key findings validated:
- Matrix-free eliminates Jacobian assembly cost
- Anderson acceleration reduces iteration count
- GPU acceleration for large problems (memory bandwidth bound)
- GMRES with adaptive tolerance (Eisenstat-Walker formula)
Includes both CPU and GPU implementations with performance comparison
showing memory usage, iteration counts, and wall-clock times for systems
ranging from 1K to 1M DOFs (835 lines, full implementation).
Complete execution output from material_models_benchmark.jl validation:
Performance results:
- Linear Elastic: 5.1× speedup (Tensors.jl vs Voigt/Dict)
- Neo-Hookean Manual: 2.0× speedup over old approach
- Perfect Plasticity: 21.0× speedup (zero allocations vs Dict)
- Average speedup: 9.4× (validates 5-50× claim range)
Key validation:
- All new implementations: ZERO allocations (confirmed)
- Manual derivatives: 21.1× faster than automatic differentiation
- Type stability: All @code_warntype checks pass (no red flags)
- AbstractMaterialState hierarchy: State handling identical for all materials
Demonstrates Newton iteration state handling for both stateless (LinearElastic,
NoState) and stateful (PerfectPlasticity, PlasticityState) materials.
Compares 5 strategies for accessing integration points during assembly:
- OLD: Runtime dispatch + mutable struct with Dict (type-unstable)
- Option A: Compile-time function returning tuples
- Option B: Store in element as NTuple (current approach)
- Option C: Compile-time with Vec{D} from Tensors.jl
- Option D: Pre-computed global constants
- Option E: Function returning pre-computed constant
Validates that compile-time generation (Options C-E) matches the golden
standard architecture from nodal assembly demos. Includes realistic FEM
assembly comparison showing performance difference between old runtime
dispatch and new compile-time approach.
Recommendation: Option C (compile-time with Vec) or D/E (pre-computed)
for zero allocation and full inlining, matching eval_basis! pattern.
Compares two state update strategies for GPU optimization:
- Strategy 1: Array of Structs (AoS) - immutable elements with embedded state
- Strategy 2: Structure of Arrays (SoA) - separate geometry and mutable state
Validates that SoA achieves 5-10× better memory bandwidth due to coalesced
access patterns. Benchmarks both CPU and GPU implementations with detailed
performance metrics including bandwidth utilization.
Key findings:
- SoA enables coalesced memory access (consecutive threads → consecutive memory)
- AoS suffers from pointer chasing and non-coalesced access
- SoA has zero allocations (in-place updates vs element reconstruction)
- Tests with 10k, 100k, 1M elements showing scalability
- Validates zero-allocation claim for compute_deformation_gradient()
- Allocation analysis with @allocated macro
- Performance benchmarking with BenchmarkTools
- LLVM IR analysis for optimization verification
- Tests finite strain formulation F = I + ∇u
- Hex8 element with 10% stretch in x-direction
- Jacobian computation J = ∑ X_i ⊗ dN_i/dξ
- Measures median time in nanoseconds/microseconds
- Confirms type stability and inlineability
- 240 lines analyzing deformation gradient computation performance
- Implements BOTH assembly methods for direct comparison on same problem
- Traditional element assembly: builds 12×12 K_e matrices, scatters to global K
- Nodal assembly: computes 3×3 K_ij blocks directly, accumulates per node
- Test problem: linear elasticity on simple Tet4 mesh
- TestLinearElastic material with Lamé parameters (λ, μ from E, ν)
- B-matrix computation: strain-displacement operator (6×3 per node, Voigt notation)
- Element stiffness: K_e = ∫ B^T C B dV with Gauss integration
- Nodal contribution: spider pattern with 3×3 blocks for coupled nodes
- Matrix-vector product comparison: K*v computed both ways
- Validates numerical equivalence: ‖K_element - K_nodal‖ < tol
- Performance characteristics: element (matrix scatter) vs nodal (direct blocks)
- Architectural differences demonstration: gather-scatter vs direct accumulation
- 542 lines validating nodal assembly correctness and comparing approaches
- Tests compute_deformation_gradient() for both FiniteStrain and SmallStrain formulations
- Identity case: u=0 → F=I, det(F)=1
- Pure translation: constant u → ∇u=0 → F=I (rigid body motion)
- Pure stretch: uniaxial extension (10%, 20%) → diagonal F
- Simple shear: u_x = γ·y → off-diagonal F components
- Validates F = I + ∇u (finite strain) vs F = I (small strain approximation)
- Physical constraint: det(F) > 0 (orientation preservation)
- Incompressibility check: det(F) ≈ 1 for volume-preserving deformation
- Symmetry verification for Right Cauchy-Green tensor C = F^T·F
- Type stability and zero allocation checks
- Integration with new API: get_basis_derivatives(Hexahedron(), Lagrange{}, ξ)
- Tests Hex8 elements with various deformation patterns
- 393 lines validating fundamental kinematics with Tensors.jl
- Demonstrates correct new API usage: Topology + Basis + IntegrationPoint separation
- Mock BasisValues struct with shape functions N and derivatives dN_dξ (SVector)
- Tests linear tetrahedron (P1, 4 nodes) evaluation at center and corner nodes
- Tests linear triangle (P1, 3 nodes) evaluation and partition of unity
- Validates constant derivatives for linear elements
- Integration with Gauss quadrature: evaluate_basis at all integration points
- Complete FEM workflow demonstration: Topology → Integration → Basis → Assembly
- Multiple element types from same topology (P1 vs P2 with same integration points)
- Type stability and zero allocation verification with StaticArrays
- 291 lines demonstrating separation of concerns: Topology ≠ Basis ≠ Integration
- Tests NeoHookean construction with both Lamé parameters and engineering constants
- Strain energy computation: reference state (ψ=0), uniaxial extension, invalid deformations
- Stress computation: small deformation, large deformation (50% extension), pure shear
- Second Piola-Kirchhoff stress: S = 2·∂ψ/∂C computed via automatic differentiation
- Tangent modulus validation: 4th-order symmetric tensor, finite difference consistency
- Verifies stress-energy relationship: S = 2·gradient(strain_energy, C)
- Small strain limit: Neo-Hookean → linear elasticity as ε → 0
- Incompressibility check for nearly incompressible materials (ν → 0.5)
- Automatic differentiation accuracy verification
- Zero allocation and type stability checks
- 295 lines validating finite deformation hyperelasticity with Tensors.jl
- Tests compute_jacobian() for 2D triangles and 3D tetrahedra
- Validates identity, scaling, and rotation transformations
- Tests physical_derivatives() conversion from reference to physical coordinates
- Verifies constant strain condition (∑ dNᵢ/dx = 0)
- Element quality checks via determinant (positive = proper orientation)
- Detects degenerate elements (det ≈ 0)
- Type stability and zero allocation verification
- Manual calculation consistency checks for known Jacobians
- Tests both tuple and vector interfaces
- 261 lines covering fundamental isoparametric mapping operations
Unit tests for FiniteStrainPlasticity with multiplicative decomposition.
Test coverage:
- Material construction with validation (E, ν, σ_y, H parameters)
- State initialization (F_p, α_bar, κ)
- Small strain limit verification
- Identity and pure rotation deformation (frame indifference)
- Uniaxial extension (elastic and plastic regimes)
- Simple shear deformation
- Incremental loading with state persistence
- Plastic incompressibility constraint (det(F_p) ≈ 1)
- Kinematic hardening behavior (backstress evolution)
- State persistence across load steps
- Type stability verification
Validates core assembly implementation by solving single Tet10 element
under uniaxial tension and comparing to analytical solution.
Test coverage:
- Linear elastic material model validation
- Strain computation from gradients (uniaxial extension)
- Assembly helpers zero allocation verification
- Type stability verification
- Stiffness matrix properties (symmetry, positive definiteness)
- Internal forces accumulation
Validates complete assembly infrastructure works correctly.
- New API: get_basis_functions() returns tuple of functions
- New API: get_basis_derivatives() returns tuple of gradient functions
- basis_api.jl: 210 lines implementing modern functional API
- Re-generated lagrange_generated.jl with 242 new lines
- abstract.jl: Add nnodes() method for Lagrange type
- Backward compatible: old eval_basis! API unchanged
- See ADR-003 for design rationale
- Implement update() that returns new element (immutable pattern)
- Supports keyword arguments for ergonomic field updates
- Preserves backward compatibility with update!() (legacy)
- Dual-API approach: modern immutable + legacy mutable both supported
- 82 lines including documentation and examples
- See docs/book/fundamentals_element_creation.md for usage guide
- Add new basis API exports: get_basis_functions, get_basis_derivatives
- Export both update() (immutable) and update!() (legacy)
- Re-enable assemble! and postprocess! exports
- Document new basis API with ADR-003 reference
- Include basis_api.jl for dual-API support (modern + legacy)
- Complete test suite for ElasticityPhysics solver
- Cantilever beam mesh: 190 nodes, 434 Tet4 elements
- Gmsh-generated mesh file (cantilever_beam.msh)
- Material: Steel (E=200 GPa, ν=0.3)
- Boundary conditions: Fixed end, pressure load on free end
- Validates convergence and displacement field
- Test passes: 430 CG iterations, max displacement 4.1 cm
Jukka Aho