New file: test/domains/continuum/test_kernel_functions.jl
Tests for:
- ContinuumKernel construction
- compute_block_at_point (weak form at single point)
- Validates stiffness contribution at integration point
- Checks tensor dimensions and symmetry
Validates:
- Kernel properly wraps formulation and material
- compute_block_at_point produces symmetric Tensor{2,3}
- Integration point contributions are reasonable magnitude
- No NaN or Inf values
Low-level validation of the atomic weak form operation
before integration loop aggregation.
New file: test/domains/continuum/test_helpers.jl
Helper functions:
- create_test_mesh_tet4(n) - generates n-element Tet4 mesh
- create_test_kernel() - creates LinearElastic kernel
- create_test_cache() - creates COOCache
- setup_assembly_test(n) - complete setup in one call
Purpose:
- Reduce code duplication across test files
- Provide consistent test data
- Make tests more readable
- Easy to extend for other element types
Used by test_full_assembly.jl, test_compute_block.jl, etc.
New file: test/domains/continuum/test_full_assembly.jl
Tests for:
- Complete assembly workflow from mesh to sparse matrix
- COO assembly with LinearElastic material
- Validates K matrix properties (symmetric, positive definite)
- Validates force vector dimensions
- Zero allocation verification
End-to-end test:
- Create mesh (Tet4 elements)
- Create kernel (ContinuumKernel + LinearElastic)
- Create assembler and cache
- Call assemble!
- Extract K and f
- Validate results
Critical integration test ensuring all components work together:
- Cache updates (3 phases)
- Compute blocks (integration)
- Scatter operations (triplets)
- Sparse matrix construction
This is the PRIMARY validation that the entire assembly pipeline
produces correct results with zero allocations.
New file: test/domains/continuum/test_dofs_per_node.jl
Tests for:
- dofs_per_node(field) for various field types
- Displacement{3} → 3 DOFs per node
- Temperature → 1 DOF per node
- DisplacementRotation → 6 DOFs per node (future)
Validates correct DOF count for field types used in
DOF mapping and cache dimensioning.
New file: test/domains/continuum/test_dof_mapping.jl
Tests for:
- get_dof_mapping!(dofs, node_ids, field)
- Maps node IDs to global DOF indices
- Validates displacement field (3 DOF/node)
- Validates temperature field (1 DOF/node)
Test cases:
- Single node → [1, 2, 3] for displacement
- Multiple nodes → [1,2,3, 4,5,6, ...] sequential DOFs
- Validates no off-by-one errors
- Validates correct DOF ordering
Ensures DOF mapping used in update_element_cache! is correct.
New file: src/domains/continuum/update_material_cache.jl (247 lines)
Features:
- update_material_cache!(material_cache, kernel, geometry_cache, ...)
- Computes stress tensors at integration points
- Computes tangent modulus tensors
- Updates material state for history-dependent materials
- Part of three-phase cache update pattern
Phase 3 of assembly (Material evaluation):
- Loop over integration points
- Compute strain tensor from ∇N and displacements
- Call material.compute_stress(ε, state_old, Δt)
- Store σ (stress) and 𝔻 (tangent modulus)
- Update state_new for next increment
Implementation:
- Handles linear case (u_global = nothing)
- Handles nonlinear case (with displacement field)
- Calls compute_stress (not inlined, complex material law)
- Stores results in material_cache.σ and material_cache.𝔻
State management:
- state_old: material state at start of increment
- state_new: material state at end of increment
- After convergence: state_old ← state_new
This is the THIRD of three cache updates called per element:
1. update_element_cache! (DOF mapping)
2. update_geometry_cache! (Jacobian, gradients)
3. update_material_cache! (stress, tangent) ← THIS FILE
After these three updates, compute_block! uses the caches to
compute element stiffness blocks K_kl.
New file: src/domains/continuum/update_geometry_cache.jl (239 lines)
Features:
- update_geometry_cache!(geometry_cache, mesh, nodes, basis)
- Computes shape function gradients at integration points
- Computes Jacobian determinants with quadrature weights
- Part of three-phase cache update pattern
Phase 2 of assembly (Geometry preprocessing):
- Extract element node coordinates
- Evaluate basis function gradients ∇N at each integration point
- Compute Jacobian matrix J and determinant det(J)
- Multiply det(J) × weight → detJ_w for integration
- Transform ∇N from parent to physical space
This is the SECOND of three cache updates called per element:
1. update_element_cache! (DOF mapping)
2. update_geometry_cache! (Jacobian, gradients) ← THIS FILE
3. update_material_cache! (stress, tangent)
New file: src/domains/continuum/update_element_cache.jl (179 lines)
Features:
- update_element_cache!(element_cache, elem_id, mesh, N, field)
- Updates DOF mapping for element nodes
- Resets K_blocks and f_blocks to zero
- Part of three-phase cache update pattern
Phase 1 of assembly (DOF mapping):
- Extract element connectivity
- Map node IDs to global DOFs
- Store in element_cache.dofs
- Clear previous element's data
Implementation:
- Uses get_dof_mapping! for field-to-DOF conversion
- Handles displacement (3 DOF/node) and other fields
- fill! for zero initialization
- @inline for performance
This is the FIRST of three cache updates called per element:
1. update_element_cache! (DOF mapping) ← THIS FILE
2. update_geometry_cache! (Jacobian, gradients)
3. update_material_cache! (stress, tangent)
New file: src/domains/continuum/types.jl
Concrete types defined:
- EmptyState <: AbstractMaterialState (for stateless materials)
- Future: Can add J2PlasticityState, DamageState, etc.
Purpose:
- Centralize material state types
- Separate from abstract interface definitions
- Enable type-stable state management in MaterialStateCache
EmptyState used by LinearElastic material (no history variables).
Plastic materials will have custom state types with fields for
equivalent plastic strain, back stress, etc.
New file: src/domains/continuum/abstract.jl
Abstract types defined:
- AbstractContinuumTheory - Supertype for FullThreeD, PlaneStress, etc.
- AbstractKernel - Supertype for ContinuumKernel and future variants
- AbstractMaterialState - Supertype for material state (EmptyState, plasticity, etc.)
Purpose:
- Establish type hierarchy for dispatch
- Document interface expectations
- Enable future extensions (shells, beams, etc.)
These were previously defined in other files, now centralized
for clarity and maintainability.
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
New functions:
- get_dof_mapping!(dofs, node_ids, field) - maps nodes to global DOFs
- get_field(field) - returns field object from various inputs
Features:
- Handles 1D, 2D, 3D displacement fields
- Supports temperature fields
- Validates node IDs are positive integers
- Added @inline for performance (called per element)
These functions support the new cache update architecture.
- Added @inline annotation for hot path function
- Called once per integration point per element pair
- Critical for achieving 484K elem/s throughput
- Part of selective inlining strategy (97% of max performance)
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
Changes to continuum domain:
- Added includes for abstract.jl and types.jl (new files)
- Replaced integration.jl with three update_*_cache.jl files
- Removed assemble.jl include
Changes to assemblers:
- Added includes for element_cache.jl, geometry_cache.jl, material_cache.jl
- Added exports for GeometryCache, MaterialStateCache types
- Added exports for update functions: update_geometry_cache!, update_element_cache!, update_material_cache!
Changes to fields API:
- Added exports for get_dof_mapping!, get_field functions
This restructuring separates cache management into dedicated modules
and implements the three-phase assembly pattern.
- Implement full 3D cantilever beam FEM validation
- Test LinearElastic material with known analytical solution
- Verify tip displacement against reference value
- Test assembly pipeline from mesh to solution
- Include boundary conditions (fixed end, tip load)
- Validate solver convergence and accuracy
- Document expected displacement and tolerance
- Serve as integration test for complete FEM workflow
- 359 lines of end-to-end validation test
- Test compute_stiffness_block! allocations for all materials
- Verify LinearElastic stiffness assembly is allocation-free
- Verify NeoHookean stiffness assembly is allocation-free
- Verify PerfectPlasticity stiffness assembly is allocation-free
- Test all continuum theory types (3D, PlaneStress, PlaneStrain, Axisymmetric)
- Use @test @allocations macro for precise allocation tracking
- Validate material tangent computation maintains zero allocations
- 260 lines of stiffness assembly allocation tests
- Test compute_stress! allocations for all material types
- Verify LinearElastic kernel is allocation-free
- Verify NeoHookean kernel is allocation-free
- Verify PerfectPlasticity kernel is allocation-free
- Test all continuum theory types (3D, PlaneStress, PlaneStrain, Axisymmetric)
- Use @test @allocations macro for precise allocation tracking
- Ensure material trait dispatch maintains zero allocations
- 257 lines of allocation verification tests
- Implement Discrete Kirchhoff Triangle (DKT) plate bending element
- Define DKTPlate formulation type with material and thickness parameters
- Implement assemble_stiffness! for plate bending problems
- Compute element stiffness matrix using DKT basis functions
- Support transverse displacement (w) and rotation (θx, θy) DOFs
- Include numerical integration over triangular domain
- Implement element force vector assembly
- Support distributed and point loads on plate surface
- Document DKT theory and implementation details
- 645 lines of complete DKT plate element implementation
- Implement assemble_stiffness! with MaterialBehavior trait dispatch
- Support StatelessStrainDependent materials (LinearElastic, NeoHookean)
- Support StatefulStrainDependent materials (PerfectPlasticity)
- Implement zero-allocation element stiffness assembly
- Use generic material kernel integration
- Replace material-specific assembly functions with unified implementation
- Include integration point loops with Jacobian computation
- Support all continuum theory types (3D, PlaneStress, PlaneStrain, Axisymmetric)
- 522 lines of generic continuum assembly implementation
- Define NodalAssembly type for nodal force assembly
- Implement direct nodal force vector accumulation
- Support pre-allocated buffers for zero-allocation assembly
- Provide nodal-to-global DOF mapping
- Include nodal load and constraint data structures
- Document nodal assembly workflow for point loads and BCs
- 234 lines of nodal assembly infrastructure
- Define ElementAssembly type for element matrix/vector assembly
- Implement local stiffness matrix and force vector containers
- Support pre-allocated buffers for zero-allocation assembly
- Provide DOF connectivity and element-to-global mapping
- Include element-level integration point data structures
- Document element assembly workflow and memory layout
- 341 lines of element assembly infrastructure
- Move legacy Problem-based assembly to src/legacy/
- Maintain backward compatibility for existing code
- Document deprecation path to new Physics-based API
- Preserve assembly_problem!, solve_problem! functions
- Support legacy element and boundary condition patterns
- Include migration guide in deprecation warnings
- 478 lines of legacy assembly implementation
- Define common assembly patterns for all element types
- Implement element-level and global assembly helpers
- Support both sparse and dense assembly strategies
- Provide integration point loop abstractions
- Include DOF mapping and scatter operations
- Document assembly workflow for structural elements
- 201 lines of framework infrastructure
- Implement Discrete Kirchhoff Triangle shape functions
- Compute rotation field interpolation with C1 continuity
- Calculate bending strain-displacement matrix
- Support transverse displacement and rotation DOFs
- Include shape function derivatives for plate bending
- Implement discrete Kirchhoff constraints at element level
- 416 lines with comprehensive DKT formulation
- Define AbstractPlateElement abstract type hierarchy
- Implement element assembly interface for plate structures
- Export DKT (Discrete Kirchhoff Triangle) plate element
- Document thin plate theory (Kirchhoff assumptions)
- Support bending and transverse shear
- Include rotation DOF handling for plate kinematics
- 188 lines of API definitions and exports
Jukka Aho