Major refactoring: migrate material cache updates to new workspace API
and add support for GlobalMaterialCache alongside legacy API.
- Change all material_cache to material_workspace parameter names
- Remove M<:AbstractMaterialState type constraints from functions
- Change state handling from EmptyState() to NamedTuple()
- Change material_cache.σ/𝔻/states to material_workspace.fields/states
- Add new overloads for GlobalMaterialCache API (NEW API)
- Add zero-allocation optimizations using getfield and pre-allocated NamedTuples
- Add backward compatibility handling for legacy state types
- Update all function signatures and docstrings
- Add imports for GlobalMaterialCache and helper functions
- Add StatelessConstantTangent overload with pre-allocated NamedTuples
- Add StatelessStrainDependent and StatefulStrainDependent overloads for GlobalMaterialCache
- Keep legacy API overloads for Matrix{<:AbstractMaterialState}
Optimize update_geometry_cache! to eliminate iterator allocations and
update to use coords field.
- Add @inline to update_geometry_cache! for better inlining
- Change enumerate loop to indexed loop to avoid iterator allocation
- Add @inbounds annotation for bounds-check elimination
- Change ip.ξ to ip.coords to use new QuadraturePoint API
Optimize update_element_cache! to eliminate iterator and view allocations.
- Add @inline to update_element_cache! for better inlining
- Change enumerate loop to indexed loop to avoid iterator allocation
- Add @inbounds annotations for bounds-check elimination
- Remove @view allocation by passing array directly to get_dof_mapping!
Refactor continuum kernel to add atomic scalar computation and microkernel
interface while maintaining backward compatibility.
- Add compute_stiffness_value for atomic scalar K[k,l][α,β] computation
- Add compute_stiffness_block for D×D block building using atomic kernel
- Refactor compute_block_at_point to use compute_stiffness_block internally
- Add evaluate function for microkernel interface with field coupling dispatch
- Add basevec import for unit vector construction
- Add conditional imports for evaluate and get_tangent
- Update documentation to reflect new atomic operation structure
- Changed struct field from 'basis::Lagrange{T,1}' to 'basis::Lagrange{1}' (line 91)
- Changed constructor from 'Lagrange{T,1}()' to 'Lagrange{1}()' (line 211)
- Added comments explaining new API: basis order only, topology passed separately
Replace function stub declaration with comment pointing to canonical definition in physics/api.jl. Added note that plate formulations should implement specific methods.
Eliminates documentation replacement warning while preserving API documentation context.
Changed from ∇N_q[k] to ∇N_data[q,k] to match Matrix layout.
- Removed intermediate ∇N_q vector extraction
- Direct 2D indexing: ∇N_data[q, k] for node k at integration point q
- Affects both LinearElastic and J2Plasticity material updates
- Consistent with compute_block interface refactoring
Deleted: src/domains/continuum/theory.jl
Reason:
- Content moved to better-organized files:
* Abstract types → src/domains/continuum/abstract.jl
* Concrete types → src/domains/continuum/types.jl
- Name 'theory.jl' was too vague
- New organization clearer for type hierarchy
This file is obsolete, content preserved in new locations.
Deleted: src/domains/continuum/integration.jl
Reason:
- Monolithic integration preprocessing replaced by three-phase pattern
- Functionality split into:
* update_element_cache.jl (DOF mapping)
* update_geometry_cache.jl (Jacobian, gradients)
* update_material_cache.jl (stress, tangent)
Benefits of split:
- Better separation of concerns
- Individual testing of each phase
- Easier to optimize each phase independently
- Clearer data flow through assembly
This file is obsolete with the new cache architecture.
Deleted: src/domains/continuum/assemble.jl
Reason:
- Old assembly implementation with type instability issues
- Replaced by new architecture in src/assemblers/element_based_coo.jl
- New version achieves zero allocations and 500K elem/s
- Three-phase cache update pattern replaces monolithic approach
Migration:
- Old: Single file with mixed concerns
- New: Separate cache files + update functions + generic assemblers
This file is obsolete with the new cache architecture.
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.
- 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)
- 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
- 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
- Define AbstractShellElement abstract type hierarchy
- Implement element assembly interface for shell structures
- Export shell formulation and element types
- Document thin shell theory (Kirchhoff-Love, Reissner-Mindlin)
- Support membrane and bending coupling
- Include rotation DOF handling for shell kinematics
- 100 lines of API definitions and exports
- Define AbstractBeamElement abstract type hierarchy
- Implement element assembly interface for beam structures
- Export beam formulation and element types
- Document Euler-Bernoulli and Timoshenko beam theories
- Support 2D and 3D beam elements
- Include rotation DOF handling for beam kinematics
- 98 lines of API definitions and exports
- Define AbstractTrussElement abstract type hierarchy
- Implement element assembly interface for truss structures
- Export truss formulation and element types
- Document 1D structural element API patterns
- Support both geometric and material nonlinearity
- 79 lines of API definitions and exports
BEFORE:
- Separate compute_block! for LinearElastic (lines 152-179)
- Separate compute_block! for NeoHookean (lines 198-236)
- Separate compute_all_blocks! for each material
- Adding 100 materials = 100 copies of integration code
AFTER:
- Single generic compute_block! for ALL materials (lines 310-351)
- Single generic compute_all_blocks! for ALL materials (lines 394-406)
- Trait-based dispatch via material_behavior()
- Constant tangent optimization preserved (lines 324-332)
- Zero code duplication regardless of material count
Implementation:
- Add compute_tangent_at_point() for StatelessConstantTangent
- Add compute_tangent_at_point() for StatelessStrainDependent
- Add compute_tangent_at_point() for StatefulStrainDependent
- Generic compute_block! dispatches on material_behavior()
- Generic compute_all_blocks! calls generic compute_block!
- Type-stable at compile time via trait dispatch
Performance:
- LinearElastic: tangent computed once (O(1) material queries)
- NeoHookean: tangent at each IP (O(NIP) queries)
- PerfectPlasticity: tangent + state at each IP (O(NIP) queries)
Benefits:
- Scalable to arbitrary number of materials
- Zero allocations maintained (verified by tests)
- Type stability maintained (verified by tests)
- Single source of truth for integration logic
- Implement compute_deformation_gradient(F, u, ∇N) for finite strain
- Implement compute_green_lagrange_strain(E, F) from deformation gradient
- Implement compute_small_strain(ε, u, ∇N) for linear kinematics
- Add comprehensive documentation for kinematic measures
- Support both small strain (linear) and finite strain (nonlinear)
- Include mathematical formulations in docstrings
- 224 lines with zero-allocation tensor operations
- Define ContinuumFormulation{Theory<:AbstractContinuumTheory}
- Implement formulation constructor with theory parameter
- Document formulation as discretization strategy wrapper
- Add usage examples for all theory types
- Support dispatch on Theory type parameter
- Enable theory-specific element assembly
- 323 lines with formulation infrastructure
- Define AbstractContinuumTheory abstract type hierarchy
- Implement FullThreeD for general 3D continuum mechanics
- Implement PlaneStress for thin structures (σ_zz = 0)
- Implement PlaneStrain for long structures (ε_zz = 0)
- Implement Axisymmetric for rotationally symmetric problems
- Add Voigt notation helpers for stress/strain tensors
- Document theory assumptions and use cases
- 178 lines with comprehensive documentation
- Remove apply_neumann_bcs! and apply_dirichlet_bcs!
- Functions moved to domains/common/boundary_conditions.jl
- Keeps assemble.jl focused on matrix/vector assembly only
- Move apply_neumann_bcs! and apply_dirichlet_bcs! from continuum/assemble.jl
- Functions are domain-agnostic (work with any AbstractKernel)
- Place in domains/common/ for reuse across continuum/beams/shells/trusses
- Update to use generic dofs_per_node(kernel) instead of hardcoded 3
- Refactor assemble!() to use COOAssembler + ContinuumKernel
- Remove 1200+ lines of monolithic assembly code
- Reduce to 176 lines (93% code reduction)
- Use create_cache(), assemble!(), extract_system() from assemblers
- Keep apply_neumann_bcs!() and apply_dirichlet_bcs!() for BC handling
- 176 lines (was 1200+ lines before refactoring)
Before refactoring:
- Monolithic assembly code mixing HOW and WHAT
- Difficult to extend with new assembler strategies
- Difficult to test assembler vs kernel logic separately
- 1200+ lines of tightly coupled code
After refactoring:
- Clean separation: assembler (HOW) vs kernel (WHAT)
- Easy to swap assembler (COO ↔ CSC ↔ Nodal)
- Easy to test components independently
- 93% code reduction (176 lines)
Usage example:
physics = Physics(
ContinuumFormulation{FullThreeD}(),
Displacement{3}(),
mesh,
LinearElastic(E=210e9, ν=0.3)
)
K, f = assemble!(physics)
Validation:
- Cantilever regression test passes (6/6 tests)
- Assembly time: 854.83 ms
- Tip deflection matches baseline within 0.1%
- Zero-allocation assembly confirmed
Jukka Aho