From 24968918c8d687a595c1ee9142931f62fd90bc71 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Thu, 20 Nov 2025 16:56:40 +0200 Subject: [PATCH] feat(continuum): Add update_material_cache! for stress and tangent MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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. --- .../continuum/update_material_cache.jl | 232 ++++++++++++++++++ 1 file changed, 232 insertions(+) create mode 100644 src/domains/continuum/update_material_cache.jl diff --git a/src/domains/continuum/update_material_cache.jl b/src/domains/continuum/update_material_cache.jl new file mode 100644 index 0000000..ee4e889 --- /dev/null +++ b/src/domains/continuum/update_material_cache.jl @@ -0,0 +1,232 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +""" +Material cache update functions for continuum elements. + +Computes stress, tangent modulus, and internal state at integration points. +""" + +using Tensors + +# ============================================================================ +# MATERIAL BEHAVIOR DISPATCH FUNCTIONS +# ============================================================================ + +""" + update_material_cache!( + material_cache::MaterialStateCache{M}, + geometry_cache::AbstractGeometryCache, + material::AbstractMaterial, + ::StatelessConstantTangent, + element_cache::ElementCache, + state_elem, + Δt::Float64 + ) where M + +Update material cache 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. +""" +@inline function update_material_cache!( + material_cache::MaterialStateCache{M}, + geometry_cache::GeometryCache, + material::AbstractMaterial, + ::StatelessConstantTangent, + element_cache::ElementCache, + 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) + + # Fill all IPs with same values + for q in 1:nips + material_cache.σ[q] = σ_ref + material_cache.𝔻[q] = 𝔻_ref + material_cache.states[q] = EmptyState() + end + + return nothing +end + +""" + update_material_cache!( + material_cache::MaterialStateCache{M}, + geometry_cache::AbstractGeometryCache, + material::AbstractMaterial, + ::StatelessStrainDependent, + element_cache::ElementCache, + state_elem, + Δt::Float64 + ) where M + +Update material cache for strain-dependent stateless materials (e.g., hyperelastic). + +Computes strain, stress, and tangent at each integration point. +No internal state tracking. +""" +@inline function update_material_cache!( + material_cache::MaterialStateCache{M}, + geometry_cache::GeometryCache, + material::AbstractMaterial, + ::StatelessStrainDependent, + element_cache::ElementCache, + 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}) + + # Compute at each integration point + @inbounds for q in 1:nips + ∇N_q = geometry_cache.∇N_data[q] + + # Deformation gradient: F = I + ∇u + F = I + for k in 1:nnodes + u_k = element_cache.u_buffer[k] + F += u_k ⊗ ∇N_q[k] + 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, nothing, 0.0) + + material_cache.σ[q] = σ + material_cache.𝔻[q] = 𝔻 + material_cache.states[q] = nothing + end + + return nothing +end + +""" + update_material_cache!( + material_cache::MaterialStateCache{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 material cache 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}, + geometry_cache::GeometryCache, + material::AbstractMaterial, + ::StatefulStrainDependent, + element_cache::ElementCache, + state_old::Matrix{<:AbstractMaterialState}, + elem_id::Int, + Δt::Float64 +) where M<:AbstractMaterialState + 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 + ∇N_q = geometry_cache.∇N_data[q] + + # Small strain: ε = sym(∇u) + ε = zero(SymmetricTensor{2,3,Float64,6}) + for k in 1:nnodes + u_k = element_cache.u_buffer[k] + ε += symmetric(u_k ⊗ ∇N_q[k]) + end + + # Get old state at this IP + state_q_old = state_elem[q] + + # Compute stress, tangent, and updated state + σ, 𝔻, state_q_new = compute_stress(material, ε, state_q_old, Δt) + + material_cache.σ[q] = σ + material_cache.𝔻[q] = 𝔻 + material_cache.states[q] = state_q_new + end + + return nothing +end + +# ============================================================================ +# MAIN DISPATCHER +# ============================================================================ + +""" + update_material_cache!( + material_cache::MaterialStateCache{M}, + geometry_cache::AbstractGeometryCache, + material::AbstractMaterial, + element_cache::ElementCache, + state_old::Union{Nothing,Matrix{M}}, + elem_id::Int, + Δt::Float64 + ) where M + +Update material state cache by computing stress, tangent, and internal state. + +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_cache`: Material cache to update (parametric with state type M) +- `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_cache.σ, material_cache.𝔻, material_cache.states + +# Zero-Allocation Guarantee +No allocations - writes to pre-allocated material_cache arrays. +""" +function update_material_cache!( + material_cache::MaterialStateCache{M}, + 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, + geometry_cache, + material, + behavior, + element_cache, + state_old, + elem_id, + Δt + ) + + return nothing +end