From 35bdac5b3777000924dd219440c5e51f411d9ca7 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Sat, 9 May 2026 16:30:34 +0300 Subject: [PATCH] refactor(src): remove material_cache.jl src/assemblers/material_cache.jl | 782 --------------------------------------- 1 file changed, 782 deletions(-) --- src/assemblers/material_cache.jl | 782 ------------------------------- 1 file changed, 782 deletions(-) delete mode 100644 src/assemblers/material_cache.jl diff --git a/src/assemblers/material_cache.jl b/src/assemblers/material_cache.jl deleted file mode 100644 index 0b8635a..0000000 --- a/src/assemblers/material_cache.jl +++ /dev/null @@ -1,782 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -""" -Assembly material workspace implementations for zero-allocation assembly. - -Contains mutable (AssemblyMaterialWorkspace) and immutable (ImmutableMaterialStateCache) variants. - -**Purpose:** Per-element temporary workspace during assembly. Stores stress (σ), tangent (đ”»), -and temporary state for ONE element at a time. Reset between elements. - -**Note:** This is different from `GlobalMaterialCache`, which stores persistent state -across all elements and time steps. -""" - -using Tensors -using ..JuliaFEM: material_field_type, material_state_type, create_zero_field, create_zero_state - -""" - AssemblyMaterialWorkspace{FieldType, StateType} - -Per-element temporary workspace for material fields and state during assembly. - -**Array of Structs (AoS) Pattern**: Matches prototype implementation for zero-allocation access. - -**Compositional Design**: Field structure inferred from material's `supported_physics()` trait. - -Contains pre-allocated arrays for material fields and temporary state. -Mutated per element during assembly, then reset for next element. - -**Purpose:** Temporary workspace during stiffness matrix assembly. -**Scope:** ONE element at a time (reset between elements). -**Lifetime:** Assembly loop only (not persistent). - -**Zero-Allocation Design**: Uses Array of Structs (AoS) pattern - Vector of NamedTuples. -Each integration point has its own NamedTuple of fields, enabling cache-friendly access -when looping through IPs. - -# Type Parameters -- `FieldType`: NamedTuple type for material fields (e.g., `(σ=..., đ”»=...)` for mechanics) -- `StateType`: NamedTuple type for state (e.g., `(Δ_p=..., α=..., Îș=...)` for plasticity) - -# Fields - -- `fields::Vector{FieldType}`: Material fields at each IP [NIP] - ONE element (AoS pattern) -- `states::Vector{StateType}`: Temporary state at each IP [NIP] - ONE element - -# Zero-Allocation Usage - -**Direct field access** (zero allocation): -```julia -workspace.fields[q].σ # → Stress at IP q (0 bytes!) -workspace.fields[q].đ”» # → Tangent at IP q (0 bytes!) -``` - -**Update pattern** (pre-create NamedTuple outside hot loop): -```julia -# Pre-create NamedTuple ONCE (outside hot loop) -fields_ref = (σ=σ_val, đ”»=đ”»_val) # ~896 bytes, but only once - -# In hot loop - reuse same NamedTuple (zero allocation) -for q in 1:nips - workspace.fields[q] = fields_ref # ~0-36 bytes (just assignment) -end -``` - -# Examples -```julia -# Mechanics only -workspace = create_material_cache(LinearElastic(...), 8) - -# Access fields -σ = workspace.fields[1].σ # → Stress at IP 1 (0 bytes!) -đ”» = workspace.fields[1].đ”» # → Tangent at IP 1 (0 bytes!) - -# Multiphysics (future) -workspace = create_material_cache(ThermoElastic(...), 8) -σ = workspace.fields[1].σ # → Stress (0 bytes!) -đ”» = workspace.fields[1].đ”» # → Tangent (0 bytes!) -q = workspace.fields[1].q # → Heat flux (0 bytes!) -k = workspace.fields[1].k # → Thermal conductivity (0 bytes!) -``` - -**Implementation:** Uses mutable struct with Vector of NamedTuples. Access via compile-time -known struct field indices for zero-allocation reads. Updates reuse pre-created NamedTuples -for zero-allocation writes. - -# See Also -- `GlobalMaterialCache`: Persistent state storage for time-stepping (all elements) -- `material_field_type()`: Trait function to infer field structure -""" -# AssemblyMaterialWorkspace uses Array of Structs (AoS) pattern matching prototype -# Each integration point has its own field container (better cache locality) -# Uses mutable struct wrapper to enable zero-allocation in-place updates -@generated function _create_field_container_type(::Type{FieldType}) where {FieldType<:NamedTuple} - field_names = fieldnames(FieldType) - field_types = [fieldtype(FieldType, name) for name in field_names] - - # Create mutable struct with same fields as FieldType - struct_fields = Expr[] - for (name, T) in zip(field_names, field_types) - push!(struct_fields, Expr(:(::), name, T)) - end - - struct_name = Symbol("FieldContainer_$(hash(FieldType))") - struct_def = Expr(:struct, true, :($struct_name), Expr(:block, struct_fields...)) - - return struct_def -end - -# AssemblyMaterialWorkspace uses Array of Structs (AoS) pattern matching prototype -# Each integration point has its own NamedTuple of fields (better cache locality) -struct AssemblyMaterialWorkspace{FieldType<:NamedTuple, StateType<:NamedTuple} <: AbstractMaterialStateCache{FieldType, StateType} - fields::Vector{FieldType} # Vector of NamedTuples - one per integration point (AoS pattern) - states::Vector{StateType} # Temporary state at each IP -end - -# ============================================================================ -# ZERO-ALLOCATION FIELD EXTRACTION HELPERS -# ============================================================================ - -# ============================================================================ -# Macro-based zero-allocation field access -# ============================================================================ - -""" - @field_vector(workspace, field_name) - -Macro to extract field vector with ZERO allocations using compile-time field index lookup. - -This macro generates code that uses `getfield` with compile-time constant indices, -completely bypassing NamedTuple property access overhead. - -# Examples -```julia -workspace = create_material_cache(LinearElastic(...), 8) - -# Zero-allocation vector extraction -đ”»_vec = @field_vector(workspace, :đ”») # → Vector{SymmetricTensor{4,3,Float64,36}} -σ_vec = @field_vector(workspace, :σ) # → Vector{SymmetricTensor{2,3,Float64,6}} - -# Then use in hot loops -for q in 1:8 - C = đ”»_vec[q] # Zero allocation! -end -``` -""" -# Helper @generated function that generates zero-allocation field access code -# This is called by the macro to generate compile-time constant getfield calls -# CRITICAL: The generated code uses nested getfield with compile-time constant indices -# This should be zero-allocation if the compiler can infer types properly -@generated function _get_field_vector_impl(workspace::AssemblyMaterialWorkspace{FieldType}, ::Val{FieldName}) where {FieldType, FieldName} - # Check if FieldType has this field - if hasfield(FieldType, FieldName) - # Find field index in FieldType NamedTuple (compile-time!) - field_names = fieldnames(FieldType) - field_idx = findfirst(==(FieldName), field_names) - - if field_idx === nothing - error("FieldType $FieldType does not have field :$FieldName") - end - - # Get the field type for type stability - field_type = fieldtype(FieldType, FieldName) - vec_type = Vector{field_type} - - # Generate code that extracts field from each element in workspace.fields - # workspace.fields is Vector{FieldType}, where FieldType is a NamedTuple - # We need to extract field FieldName from each NamedTuple in the vector - # NOTE: This still allocates a new Vector, but it's the same as get_tangent_vector - # The benefit is compile-time field index lookup (type stability) - # For true zero-allocation, we'd need to pre-allocate a buffer in the cache - return quote - # Extract vector by accessing field at compile-time known index - # This allocates a new Vector (same as get_tangent_vector), but with type stability - n = length(workspace.fields) - result = Vector{$field_type}(undef, n) - @inbounds for i in 1:n - result[i] = getfield(workspace.fields[i], $field_idx) - end - return result::$vec_type - end - else - error("FieldType $FieldType does not have field :$FieldName") - end -end - -export @field_vector - -macro field_vector(workspace, field_name) - # Extract the Symbol from the field_name argument - # Handle :field_name (QuoteNode), field_name (Symbol), and quoted expressions - field_sym = if field_name isa QuoteNode - field_name.value - elseif field_name isa Symbol - field_name - elseif field_name isa Expr && field_name.head == :quote && length(field_name.args) == 1 - field_name.args[1] - elseif field_name isa Expr && field_name.head == :macrocall - # Handle @doc macro expansion - skip it - return nothing - else - error("Expected Symbol, QuoteNode, or quoted Symbol, got $(typeof(field_name)): $field_name") - end - - # Generate code that calls the @generated function - # The @generated function will specialize on the workspace type and field name - return :(_get_field_vector_impl($(esc(workspace)), Val($(QuoteNode(field_sym))))) -end - -""" - get_tangent_vector(workspace::AssemblyMaterialWorkspace, buffer::Vector) -> Vector - -Extract tangent vector from AoS structure using pre-allocated buffer (zero-allocation). - -Updates buffer in-place and returns reference to buffer. -This eliminates allocations from Vector() constructor in list comprehension. - -# Arguments -- `workspace`: Assembly material workspace -- `buffer`: Pre-allocated buffer (must have length >= length(workspace.fields)) - -# Returns -- Reference to buffer (updated in-place) - -# Zero-Allocation -Direct assignment to buffer elements is zero-allocation (no Vector() constructor). -""" -@inline function get_tangent_vector( - workspace::AssemblyMaterialWorkspace{FieldType}, - buffer::Vector{T} -) where {FieldType, T} - if !hasfield(FieldType, :đ”») - error("FieldType $FieldType does not have :đ”» field") - end - # Update buffer in-place (zero allocation - direct assignment) - n = length(workspace.fields) - @inbounds for i in 1:n - buffer[i] = workspace.fields[i].đ”» - end - return buffer -end - -""" - get_tangent_vector(workspace::AssemblyMaterialWorkspace) -> Vector - -Extract tangent vector from AoS structure (allocates new Vector). - -**DEPRECATED**: Use `get_tangent_vector(workspace, buffer)` with pre-allocated buffer -for zero-allocation access. - -This version allocates a new Vector via list comprehension. -""" -function get_tangent_vector(workspace::AssemblyMaterialWorkspace{FieldType}) where {FieldType} - if !hasfield(FieldType, :đ”») - error("FieldType $FieldType does not have :đ”» field") - end - # Extract vector by accessing each field's đ”» component - # This allocates once when called, but is outside the hot loop - return [workspace.fields[i].đ”» for i in 1:length(workspace.fields)] -end - -""" - get_stress_vector(workspace::AssemblyMaterialWorkspace) -> Vector - -Extract stress vector from AoS structure. -Creates vector by extracting σ from each field - called once outside hot loop. -""" -function get_stress_vector(workspace::AssemblyMaterialWorkspace{FieldType}) where {FieldType} - if !hasfield(FieldType, :σ) - error("FieldType $FieldType does not have :σ field") - end - # Extract vector by accessing each field's σ component - # This allocates once when called, but is outside the hot loop - return [workspace.fields[i].σ for i in 1:length(workspace.fields)] -end - -# ============================================================================ -# ZERO-COST FIELD ACCESS VIA getproperty (COMPILE-TIME MAGIC!) -# ============================================================================ - -""" - workspace.σ[ip] # Zero-cost field access! - -Enable natural field access syntax with zero-allocation using `@generated` functions. - -# Examples -```julia -workspace = create_material_cache(LinearElastic(...), 8) - -# Natural syntax - zero allocation! -σ = workspace.σ[1] # → Stress at IP 1 (0 bytes!) -đ”» = workspace.đ”»[1] # → Tangent at IP 1 (0 bytes!) - -# Works for multiphysics too -q = workspace.q[1] # → Heat flux (0 bytes!) -k = workspace.k[1] # → Thermal conductivity (0 bytes!) -``` - -# Implementation -Uses `@generated` functions with `Base.getproperty` to enable compile-time field lookup. -The field name is known at compile time, so we generate direct field access code. -""" -# Use @generated for compile-time field lookup -# Generate specialized methods for each field name at compile time -@generated function Base.getproperty(workspace::AssemblyMaterialWorkspace{FieldType}, name::Val{Name}) where {FieldType, Name} - # Check if FieldType has this field - if hasfield(FieldType, Name) - # Find field index in FieldType NamedTuple - field_names = fieldnames(FieldType) - field_idx = findfirst(==(Name), field_names) - - if field_idx === nothing - error("FieldType $FieldType does not have field :$Name") - end - - # Generate code that extracts vector by accessing each field's component - # This creates a vector - called once outside hot loop - return :([getfield(workspace.fields[i], $field_idx) for i in 1:length(workspace.fields)]) - elseif Name === :fields - return :(getfield(workspace, 1)) - elseif Name === :states - return :(getfield(workspace, 2)) - else - # Field doesn't exist - generate error at compile time - return :(error("AssemblyMaterialWorkspace{$(FieldType)} has no field :$Name. Available material fields: $(fieldnames(FieldType)), struct fields: (:fields, :states)")) - end -end - -# Runtime fallback for Symbol (less efficient but works) -function Base.getproperty(workspace::AssemblyMaterialWorkspace{FieldType}, name::Symbol) where {FieldType} - # Convert to Val for compile-time dispatch - return getproperty(workspace, Val(name)) -end - -# ============================================================================ -# CONVENIENCE ACCESSORS -# ============================================================================ - -""" - get_stress(workspace::AssemblyMaterialWorkspace, ip::Int) - -Get stress tensor from workspace (mechanics only). - -# Backward Compatibility -Replaces `workspace.σ[ip]` with `get_stress(workspace, ip)`. - -# Examples -```julia -workspace = create_material_cache(LinearElastic(...), 8) -σ = get_stress(workspace, 1) # → Stress at IP 1 -``` -""" -# ============================================================================ -# ZERO-ALLOCATION FIELD ACCESSORS -# ============================================================================ - -# Accessor using @generated for compile-time field lookup -@generated function get_stress(workspace::AssemblyMaterialWorkspace{FieldType}, ip::Int) where {FieldType} - # Check if FieldType has :σ field - if hasfield(FieldType, :σ) - # Get field type for type stability - field_type = fieldtype(FieldType, :σ) - - # Find field index in FieldType NamedTuple - field_names = fieldnames(FieldType) - σ_idx = findfirst(==(:σ), field_names) - if σ_idx === nothing - error("FieldType $FieldType does not have :σ field") - end - - # Access via fields[ip].σ - zero allocation (compile-time known indices) - return :(@inbounds return getfield(workspace.fields[ip], $σ_idx)::$field_type) - else - error("FieldType $FieldType does not have :σ field") - end -end - -""" - get_tangent(workspace::AssemblyMaterialWorkspace, ip::Int) - -Get tangent modulus from workspace (mechanics only). - -# Backward Compatibility -Replaces `workspace.đ”»[ip]` with `get_tangent(workspace, ip)`. - -# Examples -```julia -workspace = create_material_cache(LinearElastic(...), 8) -đ”» = get_tangent(workspace, 1) # → Tangent at IP 1 -``` -""" -# Accessor using @generated for compile-time field lookup -# ============================================================================ -# TYPE-STABLE FIELD ACCESS HELPERS (Zero-Allocation) -# ============================================================================ - -""" - @generated function _get_tangent_field_index(::Type{FieldType}) where {FieldType<:NamedTuple} - -Get compile-time field index for `:đ”»` field in FieldType. - -Returns the field index as a compile-time constant, enabling type-stable `getfield` access. -""" -@generated function _get_tangent_field_index(::Type{FieldType}) where {FieldType<:NamedTuple} - field_names = fieldnames(FieldType) - field_idx = findfirst(==(:đ”»), field_names) - - if field_idx === nothing - error("FieldType $FieldType does not have field :đ”»") - end - - # Return the compile-time constant index - return field_idx -end - -""" - extract_tangent!(buffer::Vector{SymmetricTensor{4,3,Float64,36}}, - fields::Vector{FieldType}, - ::Type{FieldType}) where {FieldType<:NamedTuple} - -Extract tangent field `:đ”»` from fields vector into buffer (type-stable, zero-allocation). - -Uses compile-time field index lookup to avoid Symbol-based getfield which causes type instability. -""" -@inline function extract_tangent!( - buffer::Vector{SymmetricTensor{4,3,Float64,36}}, - fields::Vector{FieldType}, - ::Type{FieldType} -) where {FieldType<:NamedTuple} - # Get compile-time field index for :đ”» - field_idx = _get_tangent_field_index(FieldType) # Compile-time constant! - n = length(fields) - @inbounds for i in 1:n - # Use compile-time known index - type-stable and zero-allocation - buffer[i] = getfield(fields[i], field_idx)::SymmetricTensor{4,3,Float64,36} - end - return nothing -end - -@generated function get_tangent(workspace::AssemblyMaterialWorkspace{FieldType}, ip::Int) where {FieldType} - # Check if FieldType has :đ”» field - if hasfield(FieldType, :đ”») - # Get field type for type stability - field_type = fieldtype(FieldType, :đ”») - - # Find field index in FieldType NamedTuple - field_names = fieldnames(FieldType) - đ”»_idx = findfirst(==(:đ”»), field_names) - if đ”»_idx === nothing - error("FieldType $FieldType does not have :đ”» field") - end - - # Access via fields[ip].đ”» - zero allocation (compile-time known indices) - return :(@inbounds return getfield(workspace.fields[ip], $đ”»_idx)::$field_type) - else - error("FieldType $FieldType does not have :đ”» field") - end -end - -""" - get_field(workspace::AssemblyMaterialWorkspace, field_name::Symbol, ip::Int) - -Get any field from workspace by name. - -# Examples -```julia -workspace = create_material_cache(LinearElastic(...), 8) -get_field(workspace, :σ, 1) # → Stress -get_field(workspace, :đ”», 1) # → Tangent - -# Multiphysics -workspace = create_material_cache(ThermoElastic(...), 8) -get_field(workspace, :q, 1) # → Heat flux -get_field(workspace, :k, 1) # → Thermal conductivity -``` -""" -@generated function get_field(workspace::AssemblyMaterialWorkspace{FieldType}, field_name::Val{Name}, ip::Int) where {FieldType, Name} - # Check if FieldType has this field - if hasfield(FieldType, Name) - # Find field index in FieldType NamedTuple - field_names = fieldnames(FieldType) - field_idx = findfirst(==(Name), field_names) - if field_idx === nothing - error("FieldType $FieldType does not have field :$Name") - end - - # Get field type for type stability - field_type = fieldtype(FieldType, Name) - - # Access via fields[ip].Name - zero allocation (compile-time known indices) - return :(@inbounds return getfield(workspace.fields[ip], $field_idx)::$field_type) - else - error("FieldType $FieldType does not have field :$Name") - end -end - -# Non-generated fallback for runtime Symbol (less efficient but works) -function get_field(workspace::AssemblyMaterialWorkspace{FieldType}, field_name::Symbol, ip::Int) where {FieldType} - return get_field(workspace, Val(field_name), ip) -end - -# set_fields! - CRITICAL: This function MUST be zero-allocation -# Uses @generated function to generate code that constructs NamedTuple at compile time -# The generated code uses getfield with compile-time indices to extract values -# Compiler should optimize NamedTuple construction to zero allocation -@generated function set_fields!(workspace::AssemblyMaterialWorkspace{FieldType}, ip::Int, field_values::NamedTuple) where {FieldType} - field_names = fieldnames(FieldType) - n_fields = length(field_names) - - # Generate code that extracts values using getfield with compile-time indices - # This avoids runtime property access overhead - field_accesses = [:(getfield(field_values, $i)) for i in 1:n_fields] - - # Construct NamedTuple using compile-time known structure - # The compiler should optimize this to zero allocation if: - # 1. FieldType is known at compile time (it is, via @generated) - # 2. Field values are already allocated (they are, from compute_stress) - # 3. NamedTuple wrapper can be optimized away (compiler optimization) - names_tuple = Expr(:tuple, [QuoteNode(n) for n in field_names]...) - values_tuple = Expr(:tuple, field_accesses...) - - # Generate: workspace.fields[ip] = NamedTuple{(:σ, :đ”»)}((σ_val, đ”»_val)) - # This should be zero-allocation after compiler optimization - return :(@inbounds workspace.fields[ip] = NamedTuple{$names_tuple}($values_tuple); return nothing) -end - -""" - ImmutableMaterialStateCache{M,NIP} - -Immutable assembly material workspace using NTuple for zero-allocation access. - -Unlike `AssemblyMaterialWorkspace`, this version: -- Uses `NTuple` instead of `Vector` (stack-allocated, no heap access) -- Is immutable (must create new instance per element) -- Has **zero allocations** during cache access -- Enables full compiler optimization (sizes known at compile time) - -# Type Parameters -- `M`: Material state type (EmptyState for stateless) -- `NIP`: Number of integration points (compile-time constant) - -# Fields -- `σ::NTuple{NIP, SymmetricTensor{2,3,Float64,6}}`: Stress at each IP -- `đ”»::NTuple{NIP, SymmetricTensor{4,3,Float64,36}}`: Tangent modulus at each IP -- `states::NTuple{NIP, M}`: Internal state at each IP - -# Zero-Allocation Access - -```julia -# Indexing is zero-allocation: -tangent = cache.đ”»[q] # 0 bytes! -stress = cache.σ[q] # 0 bytes! -``` - -# Performance - -**Eliminates type instability** from `Vector` indexing: -- Before: `đ”»::SYMMETRICTENSOR{4, 3, FLOAT64}` (UPPERCASE = unstable) -- After: `đ”»::SymmetricTensor{4, 3, Float64}` (lowercase = concrete) - -**Pros:** -- Zero allocations during access -- Full compile-time type inference -- Stack-allocated (no GC pressure) - -**Cons:** -- Immutable (must create new instance per element) -- Cannot be reused across elements - -# Usage - -```julia -# Create new cache per element: -material_cache = create_material_cache( - ImmutableMaterialStateCache, - geometry_cache, material, element_cache -) - -# Then use normally in compute_block!: -K_kl = compute_block!(geometry_cache, material_cache, k, l) -``` -""" -# Legacy type - not part of new compositional design -# Use AssemblyMaterialWorkspace{FieldType, StateType} instead -struct ImmutableMaterialStateCache{M<:AbstractMaterialState,NIP} - σ::NTuple{NIP,SymmetricTensor{2,3,Float64,6}} # 6 independent components for 2nd order symmetric - đ”»::NTuple{NIP,SymmetricTensor{4,3,Float64,36}} # 36 independent components for 4th order symmetric - states::NTuple{NIP,M} -end - -""" - reset!(workspace::AssemblyMaterialWorkspace) - -Reset assembly material workspace to zero values. - -# Side Effects -Mutates all arrays in workspace to zero. -""" -function reset!(workspace::AssemblyMaterialWorkspace{FieldType, StateType}) where {FieldType, StateType} - # Reset all fields to zero - # CRITICAL FIX: Use getfield directly to avoid type instability from getproperty - fields = getfield(workspace, 1) # Direct field access - zero allocation, type-stable - states = getfield(workspace, 2) # Direct field access - zero allocation, type-stable - # CRITICAL FIX: Pre-compute zero_field and zero_state ONCE (they're constants for stateless materials) - # For StatelessConstantTangent, these are the same every time, so we can reuse them - # But we need to compute them here since FieldType and StateType are type parameters - zero_field = create_zero_field(FieldType) - zero_state = create_zero_state(StateType) - n = length(fields) # Direct length call - zero allocation - @inbounds for i in 1:n - fields[i] = zero_field - states[i] = zero_state - end - return nothing -end - -# Zero-allocation overload: Accept pre-allocated zero values to avoid create_zero_field allocation -function reset!( - workspace::AssemblyMaterialWorkspace{FieldType, StateType}, - zero_field::FieldType, - zero_state::StateType -) where {FieldType, StateType} - # Reset all fields to zero using pre-allocated values (zero-allocation) - fields = getfield(workspace, 1) # Direct field access - zero allocation, type-stable - states = getfield(workspace, 2) # Direct field access - zero allocation, type-stable - n = length(fields) # Direct length call - zero allocation - @inbounds for i in 1:n - fields[i] = zero_field - states[i] = zero_state - end - return nothing -end - -# ============================================================================ -# CONSTRUCTORS -# ============================================================================ - -""" - create_material_cache(material::M, max_nips::Int) -> AssemblyMaterialWorkspace{FieldType, StateType} - where {M <: AbstractMaterial} - -Create pre-allocated assembly material workspace with field and state types inferred from material traits. - -Uses trait system to determine: -- `FieldType` from `material_field_type(material)` (inferred from `supported_physics()`) -- `StateType` from `material_state_type(material)` (inferred from `required_state_variables()`) - -**Purpose:** Create temporary workspace for ONE element during assembly. -**Note:** For persistent state storage, use `create_global_material_cache()` instead. - -# Arguments -- `material`: Material model (type M determines field and state types) -- `max_nips`: Maximum integration points per element - -# Returns -- `AssemblyMaterialWorkspace{FieldType, StateType}` with field structure inferred from material - -# Type Stability -Return type is fully inferrable: -- `M` is concrete material type (known at compile time) -- `FieldType = material_field_type(material)` is concrete NamedTuple type (trait dispatch) -- `StateType = material_state_type(material)` is concrete NamedTuple type (trait dispatch) -- **Zero allocations** in hot loops! - -# Examples -```julia -# Stateless material (mechanics) -mat = LinearElastic(E=210e9, Îœ=0.3) -workspace = create_material_cache(mat, 8) -# → AssemblyMaterialWorkspace{(:σ, :đ”»), ()} -workspace.fields[1].σ # → Stress -workspace.fields[1].đ”» # → Tangent - -# Stateful material (mechanics with plasticity) -mat = PerfectPlasticity(E=210e9, Îœ=0.3, σ_y=250e6) -workspace = create_material_cache(mat, 8) -# → AssemblyMaterialWorkspace{(:σ, :đ”»), (:Δ_p, :α, :Îș)} -workspace.fields[1].σ # → Stress -workspace.states[1] # → (Δ_p=..., α=..., Îș=...) -``` - -# See Also -- `create_global_material_cache()`: For persistent state storage (all elements, time-stepping) -- `material_field_type()`: Trait function to infer field structure -""" -function create_material_cache(material::M, max_nips::Int) where M<:AbstractMaterial - # Infer field type from material traits - FieldType = material_field_type(material) - StateType = material_state_type(material) - - # Create zero-initialized field NamedTuple - zero_field = create_zero_field(FieldType) - - # Create Vector of NamedTuples - one per integration point (AoS pattern) - # This matches the prototype's MaterialContext pattern - fields = [zero_field for _ in 1:max_nips] - - # Create zero-initialized states - zero_state = create_zero_state(StateType) - states = [zero_state for _ in 1:max_nips] - - return AssemblyMaterialWorkspace{FieldType, StateType}(fields, states) -end - -# Backward compatibility alias -const create_assembly_workspace = create_material_cache - -""" - create_material_cache( - ::Type{ImmutableMaterialStateCache}, - geometry_cache::ImmutableGeometryCache{N,NIP}, - material::AbstractMaterial, - element_cache::ElementCache - ) -> ImmutableMaterialStateCache{M,NIP} - -Create immutable material state cache with NTuple fields (zero allocations). - -# Process -1. Compute stress/tangent at all integration points -2. Convert Vectors to NTuples (compile-time sizes) -3. Return immutable cache - -# Zero-Allocation Benefits - -Unlike mutable `AssemblyMaterialWorkspace`, this version: -- Uses NTuple (stack-allocated, no heap access) -- Enables full compiler optimization (sizes known at compile time) -- Eliminates type instability from Vector indexing - -# Example - -```julia -geometry_cache = create_geometry_cache( - ImmutableGeometryCache, element_cache, kernel, elem_id, mesh -) -material_cache = create_material_cache( - ImmutableMaterialStateCache, geometry_cache, material, element_cache -) -# Now both caches are zero-allocation! -``` -""" -function create_material_cache( - ::Type{ImmutableMaterialStateCache}, - geometry_cache::ImmutableGeometryCache{N,NIP}, - material::AbstractMaterial, - element_cache::ElementCache -) where {N,NIP} - # Compute stress and tangent at all integration points - σ_vec = Vector{SymmetricTensor{2,3,Float64,6}}(undef, NIP) - đ”»_vec = Vector{SymmetricTensor{4,3,Float64,36}}(undef, NIP) - - # Get strain field (if needed for material evaluation) - # For now, assume zero strain (elastic initialization) - # This will be updated in actual assembly loop - - if needs_state(material) - # Stateful material - states_vec = Vector{PlasticityState}(undef, NIP) - for q in 1:NIP - Δ = zero(SymmetricTensor{2,3,Float64,6}) # Zero strain - state = PlasticityState() # Initial state - σ_vec[q], đ”»_vec[q], states_vec[q] = update_material!(material, Δ, state) - end - # Convert to NTuple - σ_tuple = ntuple(i -> σ_vec[i], Val(NIP)) - đ”»_tuple = ntuple(i -> đ”»_vec[i], Val(NIP)) - states_tuple = ntuple(i -> states_vec[i], Val(NIP)) - return ImmutableMaterialStateCache{PlasticityState,NIP}(σ_tuple, đ”»_tuple, states_tuple) - else - # Stateless material - for q in 1:NIP - Δ = zero(SymmetricTensor{2,3,Float64,6}) - σ_vec[q], đ”»_vec[q] = evaluate_material(material, Δ) - end - # Convert to NTuple - σ_tuple = ntuple(i -> σ_vec[i], Val(NIP)) - đ”»_tuple = ntuple(i -> đ”»_vec[i], Val(NIP)) - states_tuple = ntuple(i -> EmptyState(), Val(NIP)) - return ImmutableMaterialStateCache{EmptyState,NIP}(σ_tuple, đ”»_tuple, states_tuple) - end -end