diff --git a/src/domains/continuum/kernel.jl b/src/domains/continuum/kernel.jl index f607167..1b1ef61 100644 --- a/src/domains/continuum/kernel.jl +++ b/src/domains/continuum/kernel.jl @@ -6,12 +6,20 @@ Continuum mechanics kernel - defines the weak form only. This module defines: 1. The ContinuumKernel type -2. The weak form: compute_block_at_point (atomic operation) +2. The weak form: compute_stiffness_value (most atomic operation) +3. Block builder: compute_stiffness_block (builds D×D blocks) +4. NEW: Microkernel interface via evaluate() wrapper Everything else (geometry preprocessing, integration, assembly, DOF mapping) belongs elsewhere. """ using Tensors +using Tensors: basevec # For unit vector construction + +# Import for microkernel interface and material cache accessors +if isdefined(Main, :JuliaFEM) && isdefined(Main.JuliaFEM, :evaluate) + import ..JuliaFEM: evaluate, get_tangent +end """ ContinuumKernel{Theory<:AbstractContinuumTheory, Mat<:AbstractMaterial} <: AbstractKernel @@ -73,94 +81,201 @@ get_field(kernel::ContinuumKernel) = kernel.field # ============================================================================ """ - compute_block_at_point( - grad_k::Vec{3}, - grad_l::Vec{3}, - C::SymmetricTensor{4,3} - ) -> Tensor{2,3} + compute_stiffness_value( + grad_k::Vec{D}, + grad_l::Vec{D}, + C::Tensor{4,D}, + α::Int, + β::Int + ) -> Float64 -Compute the weak form contribution at a single integration point. +Compute **single scalar** stiffness value K[k,l][α,β] at integration point. -This is the **atomic kernel operation** - pure weak form math, no geometry, no loops. -Given shape function gradients and material tensor, compute the 3×3 stiffness block. +This is the **most atomic kernel operation** - computes one DOF-pair contribution. -# Weak Form +# Theory -For displacement field u, the weak form of linear momentum is: +For displacement field u with components uₐ (α = 1,2,3 for 3D), the weak form is: -``` -∫_Ω δε : C : ε dV = ∫_Ω δu ⋅ b dV + ∫_∂Ω δu ⋅ t dS -``` + K[k,l][α,β] = ∫_Ω Bₖ,α : C : Bₗ,β dV where: -- ε = ½(∇u + (∇u)ᵀ) is the strain (symmetric part of displacement gradient) -- C is the 4th-order material stiffness tensor -- b is body force, t is surface traction +- Bₖ,α = ½(∇Nₖ ⊗ eα + eα ⊗ ∇Nₖ) = strain-displacement operator +- C = 4th-order elasticity tensor +- eα = unit vector in direction α -Discretizing u = ∑ Nᵢ uᵢ, the stiffness matrix coupling nodes k and l is: +Expanded in index notation: + + K[k,l][α,β] = ½ C[α,i,β,j] (∂Nₖ/∂xᵢ) (∂Nₗ/∂xⱼ) + + ½ C[i,α,β,j] (∂Nₖ/∂xᵢ) (∂Nₗ/∂xⱼ) + + ½ C[α,i,j,β] (∂Nₖ/∂xᵢ) (∂Nₗ/∂xⱼ) + + ½ C[i,α,j,β] (∂Nₖ/∂xᵢ) (∂Nₗ/∂xⱼ) -``` -K[k,l][α,β] = ∫_Ω Bₖ,α : C : Bₗ,β dV -``` - -where Bₖ,α = ½(∇Nₖ ⊗ eα + eα ⊗ ∇Nₖ) is the strain-displacement matrix. - -This function computes the integrand (before multiplying by detJ*w). +Using symmetry of C, this simplifies to the implementation below. # Arguments -- `grad_k`: Physical gradient ∇Nₖ at integration point -- `grad_l`: Physical gradient ∇Nₗ at integration point -- `C`: Material stiffness tensor (from elasticity_tensor(material) or compute_stress) +- `grad_k`: Physical gradient ∇Nₖ at integration point (Vec{D}) +- `grad_l`: Physical gradient ∇Nₗ at integration point (Vec{D}) +- `C`: Material stiffness tensor (Tensor{4,D}) +- `α`: DOF component at node k (1,2,3 for x,y,z) +- `β`: DOF component at node l (1,2,3 for x,y,z) # Returns -3×3 stiffness block K[k,l] at this integration point (before detJ*w scaling) +Scalar contribution to K[k,l][α,β] (before detJ*w scaling) # Performance -Zero allocations - all tensors stack-allocated. +Zero allocations, fully inlined, SIMD-friendly. # Example - ```julia -# At an integration point: +# At integration point: grad_k = Vec{3}((0.1, 0.2, 0.3)) grad_l = Vec{3}((0.4, 0.5, 0.6)) -C = elasticity_tensor(LinearElastic(E=210e9, ν=0.3)) +C = elasticity_tensor(material) -# Compute weak form contribution -K_kl_ip = compute_block_at_point(grad_k, grad_l, C) - -# Integrate: K[k,l] += K_kl_ip * detJ * weight +# Compute K[k,l][2,3] (y-component of node k, z-component of node l) +k_23 = compute_stiffness_value(grad_k, grad_l, C, 2, 3) ``` + +# Design Note +This is more atomic than 3×3 blocks. Caller can: +1. Build blocks: `K_kl[α,β] = compute_stiffness_value(...)` +2. Direct assembly: `K_global[dof_k_α, dof_l_β] += value * detJ * w` +""" +@inline function compute_stiffness_value( + grad_k::Vec{D,Float64}, + grad_l::Vec{D,Float64}, + C::Tensor{4,D,Float64}, + α::Int, + β::Int +) where D + # Build strain-displacement operators Bₖ,α and Bₗ,β + # These are 2nd-order tensors (D×D matrices) + + # eα and eβ are unit vectors + e_α = basevec(Vec{D}, α) + e_β = basevec(Vec{D}, β) + + # Bₖ,α = ½(∇Nₖ ⊗ eα + eα ⊗ ∇Nₖ) + B_k_α = 0.5 * (grad_k ⊗ e_α + e_α ⊗ grad_k) + + # Bₗ,β = ½(∇Nₗ ⊗ eβ + eβ ⊗ ∇Nₗ) + B_l_β = 0.5 * (grad_l ⊗ e_β + e_β ⊗ grad_l) + + # Compute weak form: K[α,β] = Bₖ,α : C : Bₗ,β + # Double contraction: sum over all indices + return dcontract(B_k_α, dcontract(C, B_l_β)) +end + +""" + compute_stiffness_block( + grad_k::Vec{D}, + grad_l::Vec{D}, + C::Tensor{4,D} + ) -> Tensor{2,D} + +Compute D×D stiffness block K[k,l] at integration point. + +This builds a block by calling the atomic `compute_stiffness_value()` kernel. +Use this when you want 3×3 (or 2×2) blocks. Use atomic kernel for direct assembly. + +# Arguments +- `grad_k`: Physical gradient ∇Nₖ +- `grad_l`: Physical gradient ∇Nₗ +- `C`: Material stiffness tensor + +# Returns +D×D tensor K[k,l] at this integration point (before detJ*w scaling) + +# Example +```julia +# 3×3 block for 3D +K_kl = compute_stiffness_block(grad_k, grad_l, C) + +# Direct assembly +K_global[3*(k-1)+1:3*k, 3*(l-1)+1:3*l] += K_kl * detJ * w +``` +""" +@inline function compute_stiffness_block( + grad_k::Vec{D,Float64}, + grad_l::Vec{D,Float64}, + C::Tensor{4,D,Float64} +) where D + # Build D×D block by calling atomic kernel + K_kl = zero(Tensor{2,D,Float64}) + + @inbounds for α in 1:D, β in 1:D + K_kl += compute_stiffness_value(grad_k, grad_l, C, α, β) * + (basevec(Vec{D}, α) ⊗ basevec(Vec{D}, β)) + end + + return K_kl +end + +""" + compute_block_at_point(grad_k, grad_l, C) -> Tensor{2,3} + +Computes 3×3 stiffness block from gradients and elasticity tensor. + +Maintains API compatibility with `src/assemblers/kernel_interface.jl` which +uses SymmetricTensor. Internally uses dimension-generic `compute_stiffness_block()`. """ @inline function compute_block_at_point( grad_k::Vec{3,Float64}, grad_l::Vec{3,Float64}, C::SymmetricTensor{4,3,Float64,36} ) - # Basis vectors for displacement components - e_1 = Vec{3}((1.0, 0.0, 0.0)) - e_2 = Vec{3}((0.0, 1.0, 0.0)) - e_3 = Vec{3}((0.0, 0.0, 1.0)) - e = (e_1, e_2, e_3) - - K_kl_ip = zero(Tensor{2,3,Float64,9}) - - # Loop over displacement components (α, β ∈ {x, y, z}) - @inbounds for α in 1:3, β in 1:3 - e_α, e_β = e[α], e[β] - - # Strain-displacement B-matrices (symmetric part of displacement gradient) - # Bₖ,α = ½(∇Nₖ ⊗ eα + eα ⊗ ∇Nₖ) - B_k_α = 0.5 * (grad_k ⊗ e_α + e_α ⊗ grad_k) - B_l_β = 0.5 * (grad_l ⊗ e_β + e_β ⊗ grad_l) - - # Weak form: K[α,β] = Bₖ,α : C : Bₗ,β - # (double contraction of 2nd-order tensors with 4th-order material tensor) - k_αβ = dcontract(B_k_α, dcontract(C, B_l_β)) - - # Accumulate to 3×3 block - K_kl_ip += k_αβ * (e_α ⊗ e_β) - end - - return K_kl_ip + # Convert SymmetricTensor to Tensor for computation + C_tensor = Tensor{4,3}(C) + return compute_stiffness_block(grad_k, grad_l, C_tensor) +end + +""" + evaluate( + kernel::ContinuumKernel{Theory}, + ::Displacement{3}, ::Displacement{3}, + k::Int, l::Int, α::Int, β::Int, + material_cache, geometry_cache, q::Int + ) -> Float64 + +Microkernel interface for continuum mechanics stiffness assembly. + +Wraps `compute_stiffness_value()` to provide dispatch-based field coupling. + +# Implementation + +Extracts gradients and elasticity tensor from caches, then calls +`compute_stiffness_value()` for the actual computation. + +# Example + +```julia +# Direct computation +grad_k = geometry_cache.∇N_data[q, k] +grad_l = geometry_cache.∇N_data[q, l] + @inbounds C = get_tangent(material_workspace, q) +value = compute_stiffness_value(grad_k, grad_l, C, α, β) + +# Microkernel interface +value = evaluate(kernel, Displacement{3}(), Displacement{3}(), + k, l, α, β, material_workspace, geometry_cache, q) +``` + +# Performance + +No overhead - compiler inlines to identical code. +""" +@inline function evaluate( + kernel::ContinuumKernel{Theory}, + ::Displacement{3}, ::Displacement{3}, + k::Int, l::Int, α::Int, β::Int, + material_cache, geometry_cache, q::Int +) where {Theory<:AbstractContinuumTheory} + # Extract from caches using indices + grad_k = geometry_cache.∇N_data[q, k] + grad_l = geometry_cache.∇N_data[q, l] + @inbounds C = get_tangent(material_cache, q) + + # Call existing implementation (no duplication!) + return compute_stiffness_value(grad_k, grad_l, C, α, β) end