direct array parameters instead of cache structs:
- ∇N_data::Matrix{Vec{3,Float64}} - shape function gradients
- detJ_w::Vector{Float64} - integration weights
- 𝔻::Vector{SymmetricTensor{4,3,Float64,36}} - material tangent
Eliminates field access overhead from geometry_cache and material_cache.
Renamed local variable detJ_w → w to avoid shadowing parameter.
Renamed local variable 𝔻 → D to avoid shadowing parameter.
Updated docstrings with new signatures and argument descriptions.
Maintains zero-allocation performance (verified via benchmarks).
New file: src/assemblers/scatter_blocks_to_triplets_symmetric_direct.jl (124 lines)
Features:
- scatter_blocks_to_triplets_symmetric_direct!(I, J, V, counter, capacity, K_blocks, dofs, N)
- Direct array access bypassing cache indirection
- Returns counter as Int (not Ref{Int})
- Zero dynamic dispatch (verified with @code_llvm)
BREAKTHROUGH OPTIMIZATION:
- Cache struct indirection causes dispatch (cache.I[idx])
- Direct array access enables full optimization (I[idx])
- Result: 18 → 0 dispatch sites, 500K elem/s achieved
Implementation:
- Same algorithm as symmetric version
- Takes raw arrays I, J, V instead of cache
- Counter passed by value, returned as Int
- All operations @inbounds and @inline
Performance impact:
- Eliminated ALL dynamic dispatch in assembly loop
- Key to achieving zero allocations
- Critical for 500K elem/s throughput
This is the PRODUCTION version used in element_based_coo.jl.
Other scatter functions kept for reference/alternative use cases.
New file: src/assemblers/scatter_blocks_to_triplets_symmetric.jl (92 lines)
Features:
- scatter_blocks_to_triplets_symmetric!(cache, K_blocks, dofs, N)
- Exploits matrix symmetry (only store upper triangle + diagonal)
- Reduces triplet count by ~50%
- Uses cache for I, J, V, counter
Implementation:
- Outer loop: k in 1:N
- Inner loop: l in k:N (only k ≤ l, upper triangle)
- Store both (i,j) and (j,i) entries for off-diagonal
- Store only (i,i) for diagonal
Optimization over general scatter:
- Half the triplets for symmetric matrices
- Lower memory usage
- Faster sparse matrix construction
Original cache-based version before direct scatter optimization.
New file: src/assemblers/scatter_blocks_to_triplets.jl (74 lines)
Features:
- scatter_blocks_to_triplets!(I, J, V, counter, capacity, K_blocks, dofs, N)
- General scatter for unsymmetric matrices
- Stores all N×N blocks as triplets
- Returns updated counter
Implementation:
- Double loop over node pairs (N × N)
- Triple loop over DOF pairs (3 × 3 per block)
- Direct triplet storage: I[idx], J[idx], V[idx]
- Capacity checking with bounds validation
Not currently used (symmetric version preferred), but available
for unsymmetric problems like convection or non-symmetric contact.
New file: src/assemblers/scatter_blocks_to_force.jl (54 lines)
Features:
- scatter_blocks_to_force!(f_global, f_blocks, dofs, N)
- Scatters element force blocks to global force vector
- Handles 3 DOFs per node (displacement field)
- Uses @inbounds and @inline for performance
Implementation:
- Double loop over nodes (N × 3 DOFs)
- Direct vector indexing f_global[dof] += f_i[component]
- Zero allocations, zero dispatch
Part of the direct scatter strategy that achieved zero-allocation
assembly. Critical for 500K elem/s performance.
New file: src/assemblers/node_based_coo.jl
Features:
- assemble! implementation for nodal assembly
- Assembles contributions node-by-node instead of element-by-element
- Uses node_to_elements connectivity
- Accumulates blocks for all elements touching each node
Architecture:
- Outer loop over nodes (not elements)
- Inner loop over elements containing each node
- Natural for contact mechanics (contact is nodal)
Status: Experimental, proof-of-concept implementation.
Not yet optimized like element-based assembly.
New file: src/assemblers/nodal_cache.jl
Features:
- NodalCache struct with NodeCache
- Node-based assembly pattern (alternative to element-based)
- Includes reset!, extract_system functions
- create_node_cache constructor
Use case:
- Nodal assembly (assemble node-by-node, not element-by-element)
- Experimental architecture for contact mechanics
- Potential for better parallelization and domain decomposition
Status: Framework in place, not yet optimized like COO assembly.
New file: src/assemblers/material_cache.jl (247 lines)
Features:
- Parametric MaterialStateCache{StateType}
- Stores stress tensors (σ)
- Stores tangent modulus tensors (𝔻)
- Stores material state history (state, state_new)
- update_material_cache! function
State management:
- EmptyState for stateless materials (LinearElastic)
- Custom state types for plasticity (J2PlasticityState, etc.)
- State evolution tracked across load increments
Type parameter:
- StateType: Material state type (EmptyState, J2PlasticityState, etc.)
- Enables type-stable state access
Also includes ImmutableMaterialStateCache for read-only views
with @inline accessor functions.
New file: src/assemblers/element_cache.jl
Features:
- Parametric struct ElementCache{Topo, Basis, IPs}
- Stores element stiffness blocks (K_blocks)
- Stores element force blocks (f_blocks)
- Stores DOF mapping (dofs)
- create_element_cache constructor
Type parameters:
- Topo: Element topology type (Tet4, Hex8, etc.)
- Basis: Basis function type (Lagrange{Tet4,1}, etc.)
- IPs: Integration points tuple type
This cache is reused across all elements, updated once per element
in the assembly loop. Part of three-phase cache update pattern.
New file: src/assemblers/csc_cache.jl
Features:
- CSCCache with pre-allocated sparse matrix structure
- Faster than COO for fixed sparsity patterns
- Includes reset!, extract_system functions
- Uses build_sparsity_pattern for initialization
Use case:
- Problems with known, unchanging sparsity pattern
- Faster assembly than COO (no sorting overhead)
- Lower memory usage (no duplicate entries)
New file: src/assemblers/coo_cache.jl (184 lines)
Features:
- Parametric struct COOCache{EC<:ElementCache, MC<:MaterialStateCache}
- Eliminates type instability from cache field accesses
- Stores triplets (I, J, V) for sparse matrix construction
- Includes reset! and extract_system functions
Performance impact:
- Enables zero allocations in assembly loop
- Required for achieving 500K elem/s throughput
- Critical optimization for type stability
Documentation includes:
- COO format explanation
- Performance characteristics
- Use cases and trade-offs
Major changes:
- Replaced cache-based scatter with direct array scatter
- Extract counter once before loop, write once after loop
- Use scatter_blocks_to_triplets_symmetric_direct! for zero dispatch
- Use scatter_blocks_to_force! for force vector assembly
- Removed Ref{Int} indirection in counter management
Performance improvements:
- Zero allocations in assembly loop (verified with benchmarks)
- Zero dynamic dispatch (verified with @code_llvm)
- 500K elements/second throughput (5× baseline improvement)
Three-phase cache update pattern:
- update_element_cache! for DOF mapping
- update_geometry_cache! for Jacobian and gradients
- update_material_cache! for stress and tangent modulus
- Added includes for coo_cache.jl, csc_cache.jl, nodal_cache.jl
- Removed old COOCache, CSCCache, NodalCache definitions (now in separate files)
- Removed old ElementCache, NodeCache definitions (moved to element_cache.jl)
- Kept only high-level cache coordination logic
- Implement NodalAssembler placeholder for future GPU implementation
- Add create_cache() stub for NodalCache creation
- Add assemble!() stub with planned algorithm documentation
- Add compute_node_contributions!() stub for node-level assembly
- Document GPU parallelization strategy (one thread per node)
- 178 lines of placeholder and documentation
Planned GPU algorithm:
1. Launch one thread per node
2. Each thread gets touching elements for its node
3. Compute contributions from all touching elements
4. Atomic add to global K, f (thread-safe on GPU)
Expected performance:
- 2-10x speedup on GPU for large problems (> 100k nodes)
- Better cache locality for nodal DOFs
- Natural parallelization pattern
Status:
- Not yet implemented
- Raises error directing users to COO/CSC assemblers
- Will require CUDA.jl or similar GPU framework
- Implement CSCAssembler using pre-built CSC structure
- Implement create_cache() for CSCCache with sparsity pattern
- Implement assemble!() with in-place merge to CSC arrays
- Implement merge_to_csc!() using two-pointer algorithm
- Implement scatter_to_force!() for force vector assembly
- 298 lines of optimized CSC assembly
Algorithm:
1. Pre-build sparsity pattern once (during cache creation)
2. Loop over elements
3. Compute element stiffness using kernel (in-place)
4. Get DOF mapping (in-place)
5. Merge Ke directly into CSC structure (two-pointer merge)
6. Accumulate fe to global force vector
Performance characteristics:
- 4.1x faster than COO
- 16.6x less memory than COO
- Best for production code and nonlinear problems
Two-pointer merge:
- Efficient in-place insertion into CSC arrays
- No sorting or duplicate removal needed
- Inspired by Ferrite.jl, adapted for JuliaFEM
Critical for performance:
- Structure reused across assembly calls
- Ideal for nonlinear iterations (Newton's method)
- Ideal for time stepping (same topology)
- Implement COOAssembler using coordinate (triplet) format
- Implement create_cache() for COOCache creation
- Implement assemble!() with zero-allocation element traversal
- Implement scatter_to_triplets!() for in-place triplet accumulation
- Implement scatter_to_force!() for force vector assembly
- 247 lines of COO assembly implementation
Algorithm:
1. Loop over elements
2. Compute element stiffness using kernel (in-place)
3. Get DOF mapping (in-place)
4. Scatter Ke to triplet arrays (I, J, V)
5. Scatter fe to global force vector
6. Build sparse matrix at end: sparse(I, J, V)
Performance characteristics:
- Baseline reference implementation (1.0x)
- Simple and robust
- Moderate memory usage
- Best for prototyping and debugging
Zero-allocation assembly:
- All arrays pre-allocated in cache
- Element cache reused for all elements
- No heap allocations during assembly loop
- Define AbstractKernel interface for domain-specific assembly
- Specify required methods: compute_element_stiffness!(), dofs_per_node(), get_dof_mapping!()
- Document zero-allocation requirements for all interface methods
- Provide comprehensive examples for continuum, plate, beam kernels
- Add validation helpers: validate_kernel_implementation()
- Document dispatch strategies for material models
- Changed dofs parameter to AbstractVector{Int} for view compatibility
- 329 lines of interface specification and validation
Interface contract:
- compute_element_stiffness!(): Write Ke, fe to ElementCache in-place
- dofs_per_node(): Return number of DOFs per node (pure function)
- get_dof_mapping!(): Fill global DOF indices to pre-allocated buffer
Design philosophy:
- Assemblers are generic (work with any kernel)
- Kernels are domain-specific (continuum, plate, beam, etc.)
- Interface enforces zero-allocation assembly
- Implement COOCache for coordinate format assembly
- Implement CSCCache for compressed sparse column assembly
- Implement NodalCache for node-based assembly (future GPU)
- Add reset!() methods for cache reuse in nonlinear iterations
- Add extract_system() methods to get K, f from caches
- Implement build_sparsity_pattern() for CSC structure pre-building
- Extract mesh type parameters at runtime for capacity estimation
- 407 lines of cache implementation
Zero-allocation guarantee:
- All arrays pre-allocated during cache creation
- Assembly calls reuse existing arrays
- Critical for nonlinear solvers and time stepping
Memory efficiency:
- COO: Triplet arrays sized for element connectivity
- CSC: Pre-built sparsity pattern, reused structure
- Nodal: Includes node-to-elements inverse connectivity
Jukka Aho