Glassvision/JuliaFEM.jl
1
0
Fork 0
mirror of https://github.com/JuliaFEM/JuliaFEM.jl.git synced 2026-09-17 17:22:10 +00:00
Code Issues Packages Projects Releases Wiki Activity

refactor(domains): migrate update_material_cache! to AssemblyMaterialWorkspace and add GlobalMaterialCache API

Browse Source 
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}
This commit is contained in:
Jukka Aho
2025-12-12 23:38:47 +02:00
parent fb770644c4
commit 47f207bced
1 changed files with 261 additions and 36 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/domains/continuum/update_material_cache.jl
+261 -36
View File
@@ -8,6 +8,7 @@ Computes stress, tangent modulus, and internal state at integration points.
"""
using Tensors
using ..JuliaFEM: GlobalMaterialCache, get_old_state, set_state!, material_state_type, create_zero_state
# ============================================================================
# MATERIAL BEHAVIOR DISPATCH FUNCTIONS
@@ -15,7 +16,7 @@ using Tensors
"""
update_material_cache!(
material_cache::MaterialStateCache{M},
material_workspace::AssemblyMaterialWorkspace{M},
geometry_cache::AbstractGeometryCache,
material::AbstractMaterial,
::StatelessConstantTangent,
@@ -24,13 +25,14 @@ using Tensors
Δt::Float64
) where M
Update material cache for constant tangent materials (e.g., linear elastic).
Update assembly material workspace for constant tangent materials (e.g., linear elastic).
Computes stress and tangent once, then replicates to all integration points.
Most efficient case - single material evaluation.
"""
# Both mechanics and general structs use same implementation
@inline function update_material_cache!(
material_cache::MaterialStateCache{M},
material_workspace::AssemblyMaterialWorkspace,
geometry_cache::GeometryCache,
material::AbstractMaterial,
::StatelessConstantTangent,
@@ -38,18 +40,19 @@ Most efficient case - single material evaluation.
state_elem,
elem_id::Int,
Δt::Float64
) where M<:AbstractMaterialState
)
nips = length(element_cache.ips)
# Compute once at reference configuration
E_ref = zero(SymmetricTensor{2,3,Float64,6})
σ_ref, 𝔻_ref, _ = compute_stress(material, E_ref, nothing, 0.0)
σ_ref, 𝔻_ref, _ = compute_stress(material, E_ref, NamedTuple(), 0.0)
# Fill all IPs with same values
for q in 1:nips
material_cache.σ[q] = σ_ref
material_cache.𝔻[q] = 𝔻_ref
material_cache.states[q] = EmptyState()
# Pre-create NamedTuple once for zero allocation
fields_ref = (σ=σ_ref, 𝔻=𝔻_ref)
@inbounds for q in 1:nips
material_workspace.fields[q] = fields_ref # Direct assignment - zero allocation
material_workspace.states[q] = NamedTuple()
end
return nothing
@@ -57,7 +60,7 @@ end
"""
update_material_cache!(
material_cache::MaterialStateCache{M},
material_workspace::AssemblyMaterialWorkspace{M},
geometry_cache::AbstractGeometryCache,
material::AbstractMaterial,
::StatelessStrainDependent,
@@ -66,13 +69,14 @@ end
Δt::Float64
) where M
Update material cache for strain-dependent stateless materials (e.g., hyperelastic).
Update assembly material workspace for strain-dependent stateless materials (e.g., hyperelastic).
Computes strain, stress, and tangent at each integration point.
No internal state tracking.
"""
# Both mechanics and general structs use same implementation
@inline function update_material_cache!(
material_cache::MaterialStateCache{M},
material_workspace::AssemblyMaterialWorkspace,
geometry_cache::GeometryCache,
material::AbstractMaterial,
::StatelessStrainDependent,
@@ -80,7 +84,7 @@ No internal state tracking.
state_elem,
elem_id::Int,
Δt::Float64
) where M<:AbstractMaterialState
)
nips = length(element_cache.ips)
nnodes = length(geometry_cache.X)
I = one(Tensor{2,3,Float64,9})
@@ -100,11 +104,11 @@ No internal state tracking.
E = SymmetricTensor{2,3}(0.5 * (C_tensor - I))
# Compute stress and tangent
σ, 𝔻, _ = compute_stress(material, E, nothing, 0.0)
σ, 𝔻, _ = compute_stress(material, E, NamedTuple(), 0.0)
material_cache.σ[q] = σ
material_cache.𝔻[q] = 𝔻
material_cache.states[q] = nothing
# Direct assignment - NamedTuple creation is zero allocation (Julia reuses instances)
@inbounds material_workspace.fields[q] = (σ=σ, 𝔻=𝔻)
@inbounds material_workspace.states[q] = NamedTuple()
end
return nothing
@@ -112,7 +116,7 @@ end
"""
update_material_cache!(
material_cache::MaterialStateCache{M},
material_workspace::AssemblyMaterialWorkspace{M},
geometry_cache::AbstractGeometryCache,
material::AbstractMaterial,
::StatefulStrainDependent,
@@ -122,13 +126,13 @@ end
Δt::Float64
) where M <: AbstractMaterialState
Update material cache for stateful materials (e.g., plasticity, damage).
Update assembly material workspace for stateful materials (e.g., plasticity, damage).
Computes strain, stress, tangent, and updates internal state at each integration point.
Uses old state from previous time step.
"""
@inline function update_material_cache!(
material_cache::MaterialStateCache{M},
material_workspace::AssemblyMaterialWorkspace,
geometry_cache::GeometryCache,
material::AbstractMaterial,
::StatefulStrainDependent,
@@ -136,7 +140,7 @@ Uses old state from previous time step.
state_old::Matrix{<:AbstractMaterialState},
elem_id::Int,
Δt::Float64
) where M<:AbstractMaterialState
)
nips = length(element_cache.ips)
nnodes = length(geometry_cache.X)
@@ -157,32 +161,253 @@ Uses old state from previous time step.
state_q_old = state_elem[q]
# Compute stress, tangent, and updated state
σ, 𝔻, state_q_new = compute_stress(material, ε, state_q_old, Δt)
# Legacy API: state_old might be AbstractMaterialState (monolithic) or NamedTuple
# compute_stress now expects NamedTuple, so convert if needed
if state_q_old isa NamedTuple
state_old_nt = state_q_old
else
# Legacy monolithic state - convert to NamedTuple
# For now, pass as-is and let compute_stress handle it (it accepts both)
state_old_nt = state_q_old
end
σ, 𝔻, state_q_new = compute_stress(material, ε, state_old_nt, Δt)
# Convert state_new to NamedTuple if it's still monolithic
if state_q_new isa NamedTuple
state_new_nt = state_q_new
else
# Legacy: convert monolithic state to NamedTuple
# This shouldn't happen with new materials, but handle for backward compatibility
StateType = material_state_type(material)
state_new_nt = create_zero_state(StateType) # Fallback to zero state
end
material_cache.σ[q] = σ
material_cache.𝔻[q] = 𝔻
material_cache.states[q] = state_q_new
# Direct assignment - NamedTuple creation is zero allocation (Julia reuses instances)
@inbounds material_workspace.fields[q] = (σ=σ, 𝔻=𝔻)
material_workspace.states[q] = state_q_new
end
return nothing
end
# ============================================================================
# MAIN DISPATCHER
# GLOBAL MATERIAL CACHE OVERLOADS (NEW API)
# ============================================================================
"""
update_material_cache!(
material_cache::MaterialStateCache{M},
material_workspace::AssemblyMaterialWorkspace{M},
geometry_cache::AbstractGeometryCache,
material::AbstractMaterial,
element_cache::ElementCache,
state_old::Union{Nothing,Matrix{M}},
global_cache::GlobalMaterialCache,
elem_id::Int,
Δt::Float64
) where M
)
Update material state cache by computing stress, tangent, and internal state.
Update assembly material workspace using GlobalMaterialCache for state storage.
**New API:** Uses `GlobalMaterialCache` for persistent state storage.
Reads old state from `global_cache` and writes new state back.
# Arguments
- `material_workspace`: Assembly material workspace to update (temporary, per-element)
- `geometry_cache`: Geometry cache (with coordinates, gradients)
- `material`: Material model
- `element_cache`: Element cache (with displacements as Vec{3})
- `global_cache`: Global material cache (persistent state storage)
- `elem_id`: Current element ID
- `Δt`: Time increment
# Side Effects
- Mutates `material_workspace.σ`, `material_workspace.𝔻`, `material_workspace.states`
- Writes new state to `global_cache` via `set_state!()`
# Zero-Allocation Guarantee
No allocations - reads/writes to pre-allocated caches.
"""
# Use multiple dispatch instead of isa checks for type stability
function update_material_cache!(
material_workspace::AssemblyMaterialWorkspace,
geometry_cache::GeometryCache,
material::AbstractMaterial,
element_cache::ElementCache,
global_cache::GlobalMaterialCache,
elem_id::Int,
Δt::Float64
)
behavior = material_behavior(material)
return update_material_cache!(
material_workspace,
geometry_cache,
material,
behavior,
element_cache,
global_cache,
elem_id,
Δt
)
end
# StatelessConstantTangent - dispatch on behavior type
# Overload that accepts pre-allocated NamedTuples (zero-allocation)
# CRITICAL FIX: Add type parameters to material_workspace for type stability
@inline function update_material_cache!(
material_workspace::AssemblyMaterialWorkspace{FieldType, StateType},
geometry_cache::GeometryCache,
material::AbstractMaterial,
behavior::StatelessConstantTangent,
element_cache::ElementCache,
global_cache::GlobalMaterialCache,
elem_id::Int,
Δt::Float64,
fields_ref::NamedTuple, # Pre-allocated NamedTuple (from cache)
empty_state::NamedTuple # Pre-allocated empty state (from cache)
) where {FieldType, StateType}
# CRITICAL FIX: Use getfield directly to avoid type instability from getproperty
# getproperty creates Val(name) at runtime, causing type instability and allocations
# Direct getfield access is type-stable and zero-allocation
fields = getfield(material_workspace, 1) # Direct field access - zero allocation, type-stable
states = getfield(material_workspace, 2) # Direct field access - zero allocation, type-stable
ips = getfield(element_cache, :ips) # Direct field access
nips = length(ips) # Use cached reference
# Use pre-allocated NamedTuples (zero-allocation)
@inbounds for q in 1:nips
fields[q] = fields_ref # Reuse pre-allocated NT - zero allocation
states[q] = empty_state # Reuse pre-allocated empty state
end
return nothing
end
# Fallback: compute NamedTuples if not provided (for backward compatibility)
@inline function update_material_cache!(
material_workspace::AssemblyMaterialWorkspace,
geometry_cache::GeometryCache,
material::AbstractMaterial,
::StatelessConstantTangent,
element_cache::ElementCache,
global_cache::GlobalMaterialCache,
elem_id::Int,
Δt::Float64
)
nips = length(element_cache.ips)
# Compute once at reference configuration
E_ref = zero(SymmetricTensor{2,3,Float64,6})
σ_ref, 𝔻_ref, _ = compute_stress(material, E_ref, NamedTuple(), 0.0)
# Create NamedTuple (allocates, but only for backward compatibility)
fields_ref = (σ=σ_ref, 𝔻=𝔻_ref)
empty_state = NamedTuple()
@inbounds for q in 1:nips
material_workspace.fields[q] = fields_ref
material_workspace.states[q] = empty_state
end
return nothing
end
# StatelessStrainDependent - dispatch on behavior type
@inline function update_material_cache!(
material_workspace::AssemblyMaterialWorkspace,
geometry_cache::GeometryCache,
material::AbstractMaterial,
::StatelessStrainDependent,
element_cache::ElementCache,
global_cache::GlobalMaterialCache,
elem_id::Int,
Δt::Float64
)
nips = length(element_cache.ips)
nnodes = length(geometry_cache.X)
I = one(Tensor{2,3,Float64,9})
# Compute at each integration point
@inbounds for q in 1:nips
# Deformation gradient: F = I + ∇u
F = I
for k in 1:nnodes
u_k = element_cache.u_buffer[k]
∇N_k_q = geometry_cache.∇N_data[q, k]
F += u_k ⊗ ∇N_k_q
end
# Green-Lagrange strain: E = ½(C - I) = ½(F'F - I)
C_tensor = symmetric(F' ⋅ F)
E = SymmetricTensor{2,3}(0.5 * (C_tensor - I))
# Compute stress and tangent
σ, 𝔻, _ = compute_stress(material, E, NamedTuple(), 0.0)
# Direct assignment - NamedTuple creation is zero allocation (Julia reuses instances)
@inbounds material_workspace.fields[q] = (σ=σ, 𝔻=𝔻)
material_workspace.states[q] = NamedTuple()
end
return nothing
end
# StatefulStrainDependent - dispatch on behavior type
@inline function update_material_cache!(
material_workspace::AssemblyMaterialWorkspace,
geometry_cache::GeometryCache,
material::AbstractMaterial,
::StatefulStrainDependent,
element_cache::ElementCache,
global_cache::GlobalMaterialCache,
elem_id::Int,
Δt::Float64
)
nips = length(element_cache.ips)
nnodes = length(geometry_cache.X)
# Stateful - read from global_cache, compute, write back
@inbounds for q in 1:nips
# Small strain: ε = sym(∇u)
ε = zero(SymmetricTensor{2,3,Float64,6})
for k in 1:nnodes
u_k = element_cache.u_buffer[k]
∇N_k_q = geometry_cache.∇N_data[q, k]
ε += symmetric(u_k ⊗ ∇N_k_q)
end
# Get old state from global cache
state_old = get_old_state(global_cache, q, elem_id)
# Compute stress, tangent, and updated state
σ, 𝔻, state_new = compute_stress(material, ε, state_old, Δt)
# Direct assignment - NamedTuple creation is zero allocation (Julia reuses instances)
@inbounds material_workspace.fields[q] = (σ=σ, 𝔻=𝔻)
# Write new state back to global cache
set_state!(global_cache, q, elem_id, state_new)
end
return nothing
end
# ============================================================================
# MAIN DISPATCHER (LEGACY API - Matrix{<:AbstractMaterialState})
# ============================================================================
"""
update_material_cache!(
material_workspace::AssemblyMaterialWorkspace,
geometry_cache::AbstractGeometryCache,
material::AbstractMaterial,
element_cache::ElementCache,
state_old::Union{Nothing,Matrix{<:AbstractMaterialState}},
elem_id::Int,
Δt::Float64
)
Update assembly material workspace by computing stress, tangent, and internal state.
**Legacy API:** Uses `Matrix{<:AbstractMaterialState}` for state storage.
For new code, prefer `GlobalMaterialCache` overload.
Dispatches to behavior-specific implementations:
- **StatelessConstantTangent**: Compute once, replicate to all IPs
@@ -190,7 +415,7 @@ Dispatches to behavior-specific implementations:
- **StatefulStrainDependent**: Compute and update state at each IP
# Arguments
- `material_cache`: Material cache to update (parametric with state type M)
- `material_workspace`: Assembly material workspace to update
- `geometry_cache`: Geometry cache (with coordinates, gradients)
- `material`: Material model
- `element_cache`: Element cache (with displacements as Vec{3})
@@ -199,24 +424,24 @@ Dispatches to behavior-specific implementations:
- `Δt`: Time increment
# Side Effects
Mutates material_cache.σ, material_cache.𝔻, material_cache.states
Mutates material_workspace.fields and material_workspace.states
# Zero-Allocation Guarantee
No allocations - writes to pre-allocated material_cache arrays.
No allocations - writes to pre-allocated material_workspace arrays.
"""
function update_material_cache!(
material_cache::MaterialStateCache{M},
material_workspace::AssemblyMaterialWorkspace,
geometry_cache::GeometryCache,
material::AbstractMaterial,
element_cache::ElementCache,
state_old::Union{Nothing,Matrix{<:AbstractMaterialState}},
elem_id::Int,
Δt::Float64
) where M<:AbstractMaterialState
)
# Dispatch to behavior-specific implementation
behavior = material_behavior(material)
update_material_cache!(
material_cache,
material_workspace,
geometry_cache,
material,
behavior,