diff --git a/src/domains/continuum/update_material_cache.jl b/src/domains/continuum/update_material_cache.jl index 069ea94..2d20a09 100644 --- a/src/domains/continuum/update_material_cache.jl +++ b/src/domains/continuum/update_material_cache.jl @@ -8,190 +8,11 @@ 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 +using ..JuliaFEM: GlobalMaterialCache, get_old_state, set_state! +using ..JuliaFEM: continuum_kinematics, SmallStrainKinematics, GreenLagrangeKinematics # ============================================================================ -# MATERIAL BEHAVIOR DISPATCH FUNCTIONS -# ============================================================================ - -""" - update_material_cache!( - material_workspace::AssemblyMaterialWorkspace{M}, - geometry_cache::AbstractGeometryCache, - material::AbstractMaterial, - ::StatelessConstantTangent, - element_cache::ElementCache, - state_elem, - Δt::Float64 - ) where M - -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_workspace::AssemblyMaterialWorkspace, - geometry_cache::GeometryCache, - material::AbstractMaterial, - ::StatelessConstantTangent, - element_cache::ElementCache, - state_elem, - 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) - - # Fill all IPs with same values - # 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 -end - -""" - update_material_cache!( - material_workspace::AssemblyMaterialWorkspace{M}, - geometry_cache::AbstractGeometryCache, - material::AbstractMaterial, - ::StatelessStrainDependent, - element_cache::ElementCache, - state_elem, - Δt::Float64 - ) where M - -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_workspace::AssemblyMaterialWorkspace, - geometry_cache::GeometryCache, - material::AbstractMaterial, - ::StatelessStrainDependent, - element_cache::ElementCache, - state_elem, - 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] = (σ=σ, 𝔻=𝔻) - @inbounds material_workspace.states[q] = NamedTuple() - end - - return nothing -end - -""" - update_material_cache!( - material_workspace::AssemblyMaterialWorkspace{M}, - geometry_cache::AbstractGeometryCache, - material::AbstractMaterial, - ::StatefulStrainDependent, - element_cache::ElementCache, - state_old::Union{Nothing,Matrix{<:AbstractMaterialState}}, - elem_id::Int, - Δt::Float64 - ) where M <: AbstractMaterialState - -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_workspace::AssemblyMaterialWorkspace, - geometry_cache::GeometryCache, - material::AbstractMaterial, - ::StatefulStrainDependent, - element_cache::ElementCache, - state_old::Matrix{<:AbstractMaterialState}, - elem_id::Int, - Δt::Float64 -) - nips = length(element_cache.ips) - nnodes = length(geometry_cache.X) - - # Extract state_old for this element - state_elem = @view state_old[:, elem_id] - - # Compute and update state at each integration point - @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 at this IP - state_q_old = state_elem[q] - - # Compute stress, tangent, and updated state - # 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 - - # 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 - -# ============================================================================ -# GLOBAL MATERIAL CACHE OVERLOADS (NEW API) +# GLOBAL MATERIAL CACHE — behavior-dispatched implementations # ============================================================================ """ @@ -205,10 +26,11 @@ end Δt::Float64 ) -Update assembly material workspace using GlobalMaterialCache for state storage. +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. +Reads the old state from `global_cache`, computes stress and tangent, and +writes the new state back. Behavior dispatches on `material_behavior(material)` +so stateless / strain-dependent / stateful materials all share this entry. # Arguments - `material_workspace`: Assembly material workspace to update (temporary, per-element) @@ -226,7 +48,6 @@ Reads old state from `global_cache` and writes new state back. # 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, @@ -234,7 +55,7 @@ function update_material_cache!( element_cache::ElementCache, global_cache::GlobalMaterialCache, elem_id::Int, - Δt::Float64 + Δt::Float64, ) behavior = material_behavior(material) return update_material_cache!( @@ -245,13 +66,13 @@ function update_material_cache!( element_cache, global_cache, elem_id, - Δt + Δ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 +# StatelessConstantTangent — compute stress and tangent once at the reference +# configuration and replicate to every IP. Both NamedTuples are pre-allocated +# in the cache, so the per-IP loop is zero-allocation. @inline function update_material_cache!( material_workspace::AssemblyMaterialWorkspace{FieldType, StateType}, geometry_cache::GeometryCache, @@ -261,27 +82,25 @@ end 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) + fields_ref::NamedTuple, + empty_state::NamedTuple, ) 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) + # `getfield` is type-stable; `getproperty` would allocate here. + fields = getfield(material_workspace, 1) + states = getfield(material_workspace, 2) + ips = getfield(element_cache, :ips) + nips = length(ips) + @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 + fields[q] = fields_ref + states[q] = empty_state end return nothing end -# Fallback: compute NamedTuples if not provided (for backward compatibility) +# Convenience overload that builds the NamedTuples on demand. Used by the +# behavior dispatcher when no pre-allocated tuples are passed in. @inline function update_material_cache!( material_workspace::AssemblyMaterialWorkspace, geometry_cache::GeometryCache, @@ -290,15 +109,13 @@ end element_cache::ElementCache, global_cache::GlobalMaterialCache, elem_id::Int, - Δt::Float64 + Δ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 @@ -309,7 +126,7 @@ end return nothing end -# StatelessStrainDependent - dispatch on behavior type +# StatelessStrainDependent — strain at each IP, no persistent state. @inline function update_material_cache!( material_workspace::AssemblyMaterialWorkspace, geometry_cache::GeometryCache, @@ -318,15 +135,14 @@ end element_cache::ElementCache, global_cache::GlobalMaterialCache, elem_id::Int, - Δt::Float64 + Δ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 + # Deformation gradient F = I + ∇u F = I for k in 1:nnodes u_k = element_cache.u_buffer[k] @@ -334,14 +150,12 @@ end F += u_k ⊗ ∇N_k_q end - # Green-Lagrange strain: E = ½(C - I) = ½(F'F - I) + # Green–Lagrange strain E = ½(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 @@ -349,7 +163,8 @@ end return nothing end -# StatefulStrainDependent - dispatch on behavior type +# StatefulStrainDependent — read old state from `global_cache`, compute the +# stress / tangent / new state at each IP, and write the new state back. @inline function update_material_cache!( material_workspace::AssemblyMaterialWorkspace, geometry_cache::GeometryCache, @@ -358,98 +173,42 @@ end element_cache::ElementCache, global_cache::GlobalMaterialCache, elem_id::Int, - Δt::Float64 + Δt::Float64, ) nips = length(element_cache.ips) nnodes = length(geometry_cache.X) + kin = continuum_kinematics(material) + I = one(Tensor{2,3,Float64,9}) - # 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) + ε_gl_or_small = if kin isa SmallStrainKinematics + ε = 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 + ε + elseif kin isa GreenLagrangeKinematics + 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 + C_tensor = symmetric(F' ⋅ F) + SymmetricTensor{2,3}(0.5 * (C_tensor - I)) + else + error("unknown continuum kinematics $kin") 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) + σ, 𝔻, state_new = compute_stress(material, ε_gl_or_small, state_old, Δt) + @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 -- **StatelessStrainDependent**: Compute at each IP (no state) -- **StatefulStrainDependent**: Compute and update state at each IP - -# Arguments -- `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}) -- `state_old`: Global material state [nips, nelems] (nothing for stateless) -- `elem_id`: Current element ID -- `Δt`: Time increment - -# Side Effects -Mutates material_workspace.fields and material_workspace.states - -# Zero-Allocation Guarantee -No allocations - writes to pre-allocated material_workspace arrays. -""" -function update_material_cache!( - material_workspace::AssemblyMaterialWorkspace, - geometry_cache::GeometryCache, - material::AbstractMaterial, - element_cache::ElementCache, - state_old::Union{Nothing,Matrix{<:AbstractMaterialState}}, - elem_id::Int, - Δt::Float64 -) - # Dispatch to behavior-specific implementation - behavior = material_behavior(material) - update_material_cache!( - material_workspace, - geometry_cache, - material, - behavior, - element_cache, - state_old, - elem_id, - Δt - ) - - return nothing -end