- [weakdeps] CUDA package for optional GPU support
- [extensions] JuliaFEMCUDAExt loaded when 'using CUDA'
- Zero overhead when CUDA not loaded
- Enables GPU acceleration without mandatory dependency
- Complete cantilever beam example from mesh to visualization
- Gmsh mesh creation with physical groups for BCs
- Material definition (Young's modulus and Poisson's ratio)
- Dirichlet (fixed) and Neumann (pressure) boundary conditions
- ElasticityPhysics problem setup and solve!() call
- Results visualization with stress and displacement
- 361 lines: Step-by-step user tutorial for beginners
- 4 design documents for GPU implementation (2588 lines total)
- gpu_physics_architecture.md: Physics{Elasticity} GPU-first design
- gpu_cpu_backend_architecture.md: Backend selection and dispatch
- gpu_cpu_migration_guide.md: Migration from old API to new
- gpu_elasticity_refactoring.md: Complete refactoring strategy
- Zero CPU-GPU transfer during solve, matrix-free CG
- Elements store geometry, no mesh dependency
- Breaking changes allowed for GPU performance
- 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
- JuliaFEMCUDAExt package extension (loaded with 'using CUDA')
- ElasticityDataGPU: All data on device (nodes, elements, BCs, node-to-elem map)
- initialize_gpu_data!() transfers Physics to GPU with renumbering
- nodal_assembly_kernel!() computes stiffness 3×3 blocks per node
- compute_element_stresses_kernel!() element-level stress computation
- apply_surface_traction_kernel!() Neumann BC on surfaces
- apply_dirichlet_kernel!() penalty method for essential BCs
- cg_solve_matfree_gpu!() matrix-free CG solver
- solve_newton_krylov_gpu!() inexact Newton with GMRES-style restart
- initialize_backend(::GPU) and solve_backend!() dispatch methods
- 934 lines: Pure GPU implementation with zero CPU-GPU transfers during solve
- NodeToElementsMap: Inverse connectivity (node → elements touching it)
- ElementNodeInfo: Tracks element ID and local node index
- get_node_spider() finds all nodes coupling with given node
- NodalStiffnessContribution: Storage for 3×3 blocks per node
- matrix_vector_product_nodal() computes K_i*u at single node
- print_spider_info() debugging diagnostics
- 234 lines: Infrastructure for node-by-node assembly
- ElementAssemblyData: Global sparse matrix and force vectors
- ElementContribution: Local element contributions before scatter
- scatter_to_global!() adds element matrices to global system
- compute_residual!() calculates r = f_int - f_ext
- apply_dirichlet_bc!() penalty method for essential BCs
- get_dof_indices() node connectivity to global DOF mapping
- matrix_vector_product() sparse K*v multiplication
- 341 lines: Traditional element-by-element assembly infrastructure
- ElasticityDataCPU struct wraps ElementAssemblyData
- initialize_backend() assembles global system from immutable Elements
- compute_element_stiffness() uses Tensors.jl (blocked by get_basis_derivatives)
- cg_solve() implements Conjugate Gradient iterative solver
- Supports Dirichlet boundary conditions from Physics API
- 228 lines: Traditional element assembly approach for CPU
New file src/backend/abstract.jl defining backend abstraction:
- AbstractBackend base type for computation backend
- Auto() automatic backend selection (GPU if available, else CPU)
- GPU() force GPU backend (errors if CUDA unavailable)
- CPU(nthreads) force CPU backend with thread count
- select_backend() chooses concrete backend based on hardware
- AbstractElasticityData for backend-specific data structures
- ElasticitySolution struct for solve results
- solve!() dispatch point with backend parameter
- 241 lines with comprehensive API documentation
New file src/materials/finite_strain_plasticity.jl implementing J2 plasticity for large deformations:
- FiniteStrainPlasticityState storing F_p (plastic deformation gradient), α_bar (backstress), κ
- FiniteStrainPlasticity struct with E, ν, σ_y, H parameters
- Hyperelastic stress response using Neo-Hookean
- Exponential map integration for plastic flow
- Pull-back/push-forward operations for intermediate configuration
- Consistent algorithmic tangent for Newton convergence
- 293 lines with comprehensive finite deformation theory
New file src/materials/neo_hookean.jl implementing simplest hyperelasticity:
- NeoHookean struct with shear modulus μ and Lamé parameter λ
- Convenience constructor from E and ν engineering constants
- strain_energy() computes ψ = μ/2·(I₁-3) - μ·ln(J) + λ/2·ln²(J)
- Stress S = 2·∂ψ/∂C via automatic differentiation
- Tangent 𝔻 = 4·∂²ψ/∂C² via automatic differentiation
- Uses Tensors.jl built-in AD (no ForwardDiff dependency)
- Total Lagrangian formulation with 2nd PK stress
- 253 lines with comprehensive theory documentation
New file src/materials/abstract_material.jl defining material model architecture:
- AbstractMaterial base type for all materials
- AbstractElasticMaterial for stateless materials (no history)
- AbstractPlasticMaterial for stateful materials (plastic strain, etc.)
- compute_stress() interface: (material, ε, state_old, Δt) → (σ, 𝔻, state_new)
- State management convention for Newton iterations
- Thread-safe and GPU-compatible design principles
- 229 lines with comprehensive documentation and examples
New file src/physics_api.jl defining user-facing elasticity API:
- ElasticityPhysicsType (alias Elasticity) for problem configuration
- DirichletBC struct for prescribed displacements
- NeumannBC struct for surface tractions/pressures
- Physics{P} container for problem with elements and BCs
- Works with both CPU and GPU backends
- 183 lines with comprehensive examples
New file src/physics/deformation_gradient.jl:
- compute_deformation_gradient() computes F = I + ∇u at integration points
- StrainFormulation types: FiniteStrain() and SmallStrain()
- Uses Tensors.jl for all tensor operations (Vec, Tensor)
- Zero-allocation design with @inline functions
- GPU-ready immutable operations
- Comprehensive mathematical documentation with references
- 243 lines including commented high-level API for future integration
New file src/physics/assembly_helpers.jl with FEM assembly utilities:
- shape_function_gradients() computes ∇N in current configuration
- compute_strain_from_gradients() small strain ε = sym(∇u)
- compute_green_lagrange_strain() finite strain E = ½(C-I)
- accumulate_stiffness!() adds element stiffness contributions
- accumulate_internal_forces!() computes f_int = ∫σ·∇N dV
- accumulate_external_forces!() computes f_ext = ∫N·b dV
- Zero-allocation design with Tensors.jl Vec and SymmetricTensor
- 331 lines with comprehensive performance documentation
New file src/physics/abstract.jl defining physics system architecture:
- AbstractPhysics base type for all physics implementations
- get_unknown_field_name() returns primary field (displacement, temperature, etc.)
- get_formulation_type() returns :incremental, :total, or :rate
- get_unknown_field_dimension() returns DOFs per node
- assemble!() dispatch point for physics-specific assembly
- Comprehensive docstrings covering multi-physics coupling and GPU compatibility
- 138 lines documenting design philosophy and future extension
New file src/geometry/jacobian.jl implementing geometric transformations:
- compute_jacobian(X, dN_dξ) computes J = ∂x/∂ξ using tensor products
- physical_derivatives(J, dN_dξ) transforms derivatives to physical space
- Full Tensors.jl integration with Vec and Tensor types
- Zero-allocation tuple-based API for performance
- AbstractVector overloads for compatibility
- Comprehensive docstrings with 2D/3D examples
- 169 lines with mathematical definitions and usage patterns
Modified src/integration/gauss.jl to fix IntegrationPoint creation:
- Changed from generator expression to ntuple for proper type inference
- Collect quad_data first (was zip iterator, cannot be indexed)
- Remove explicit type parameter {D} - let Julia infer from arguments
- Fixes type stability issue in integration point generation
- Maintains zero-allocation design with tuple return
New file implementing Gauss quadrature point generation:
- get_gauss_points!(topology, scheme) returns tuple of (weight, Vec{D}) pairs
- Supports all 7 topologies: Segment, Triangle, Quadrilateral, Tetrahedron, Hexahedron, Wedge, Pyramid
- Orders 1-3 for each topology (exact integration up to quintic/cubic)
- Uses Tensors.jl Vec types for coordinates (GPU-friendly, zero-allocation)
- Fully inlined (@inline) for compile-time optimization
- 300 lines of quadrature rules from standard FEM references
Modified src/topology/topology.jl to reflect new architecture:
- Clarify topology defines geometric shape only, not node count
- Document that node count comes from basis functions
- Add examples showing same topology with different bases (Quad4/8/9)
- Update docstring to reference new topology types (Segment, Triangle, etc.)
- Emphasize corner nodes only in topology API
- Remove references to old node-count-baked types (Tri3, Quad4, etc.)
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.
Jukka Aho