From fb10cbbf8a51e99bc9e4047804e776e9b30b269b Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Wed, 19 Nov 2025 09:03:28 +0200 Subject: [PATCH] feat(continuum): Implement Jacobian and integration utilities MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Implement compute_jacobian(X, ∇N_ξ) for coordinate mapping - Implement compute_jacobian_determinant(J) with singularity checks - Implement compute_shape_derivatives(∇N, J) in physical coordinates - Add default_integration(topology) for element-specific quadrature - Support Hex8, Tet4, Quad4, Tri3, Seg2 element types - Include integration point selection logic - 327 lines with robust numerical handling --- src/domains/continuum/integration.jl | 363 +++++++++++++++++++++++++++ 1 file changed, 363 insertions(+) create mode 100644 src/domains/continuum/integration.jl diff --git a/src/domains/continuum/integration.jl b/src/domains/continuum/integration.jl new file mode 100644 index 0000000..4745346 --- /dev/null +++ b/src/domains/continuum/integration.jl @@ -0,0 +1,363 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +""" +Continuum mechanics integration utilities. + +Provides geometry preprocessing and integration wrappers for the weak form kernel. +These functions are generic and work with any kernel that implements compute_block_at_point. +""" + +using StaticArrays +using Tensors + +# ============================================================================ +# GEOMETRY PREPROCESSING +# ============================================================================ + +""" + PreparedElement{N,NIP,GradType,WeightType} + +Precomputed element geometry for block-oriented assembly. + +Stores all Jacobian-dependent data so blocks can be computed without +recomputing shape function gradients. Created once per element by +`prepare_element!`, then passed to integration functions multiple times. + +# Type Parameters +- `N`: Number of nodes in element +- `NIP`: Number of integration points +- `GradType`: Type of gradient storage (NTuple of SVectors) +- `WeightType`: Type of integration weight storage (SVector) + +# Fields +- `X`: Node coordinates [N × Vec{3}] (stack-allocated SVector) +- `∇N_data`: Physical gradients at each IP [NIP × (N × Vec{3})] +- `detJ_w`: detJ * weight at each IP [NIP] + +# Zero-Allocation +All fields use stack-allocated StaticArrays (SVector, NTuple). +Size known at compile time → perfect type stability. +""" +struct PreparedElement{N,NIP,GradType,WeightType} + X::SVector{N,Vec{3,Float64}} + ∇N_data::GradType # NTuple{NIP, SVector{N, Vec{3}}} + detJ_w::WeightType # SVector{NIP, Float64} +end + +""" + prepare_element!( + cache::ElementCache, + kernel::ContinuumKernel, + element_id::Int, + mesh::AbstractMesh + ) -> PreparedElement + +Precompute element geometry for integration **once**. + +Computes: +- Node coordinates (from mesh) +- Physical gradients ∇N at each integration point +- Jacobian determinant × weight (detJ * w) at each integration point + +Returned `PreparedElement` can be passed to integration functions multiple times +without recomputing geometry. Essential for nodal assemblers where each node +queries multiple blocks from the same element. + +# Arguments +- `cache`: Element cache (provides topology, basis, integration points) +- `kernel`: Continuum kernel +- `element_id`: Element index in mesh +- `mesh`: Finite element mesh + +# Returns +`PreparedElement{N,NIP}` with precomputed geometry (stack-allocated) + +# Zero-Allocation +Returns immutable struct with SVector/NTuple fields → stack-only, zero heap. +""" +@inline function prepare_element!( + cache::ElementCache{T,B,IPS}, + kernel::ContinuumKernel, + element_id::Int, + mesh::AbstractMesh +) where {T<:AbstractTopology{N},B,IPS} where {N} + + conn = mesh.connectivity[element_id] + + # Load coordinates into SVector (stack-allocated, size N known at compile time) + X = SVector{N}(ntuple(i -> Vec{3}(mesh.nodes[conn[i]]), N)) + + ips = cache.ips + NIP = length(ips) + + # Precompute physical gradients at all integration points + # Use ntuple for compile-time size (returns NTuple → stack-allocated) + ∇N_data = ntuple(NIP) do ip_idx + ip = ips[ip_idx] + ξ = Vec{3}(ip.ξ) + + # Reference gradients + dN_dξ = get_basis_derivatives(cache.topology, cache.basis, ξ) + + # Jacobian: J = X ⊗ ∇_ξ N + J = X[1] ⊗ dN_dξ[1] + @inbounds for i in 2:N + J += X[i] ⊗ dN_dξ[i] + end + + J_inv_T = transpose(inv(J)) + + # Physical gradients for all nodes: ∇N = J^{-T} ⋅ ∇_ξ N + SVector{N}(ntuple(k -> J_inv_T ⋅ dN_dξ[k], N)) + end + + # Precompute detJ * weight at each integration point + detJ_w_data = SVector{NIP}(ntuple(NIP) do ip_idx + ip = ips[ip_idx] + ξ = Vec{3}(ip.ξ) + dN_dξ = get_basis_derivatives(cache.topology, cache.basis, ξ) + + J = X[1] ⊗ dN_dξ[1] + @inbounds for i in 2:N + J += X[i] ⊗ dN_dξ[i] + end + + det(J) * ip.weight + end) + + return PreparedElement{N,NIP,typeof(∇N_data),typeof(detJ_w_data)}(X, ∇N_data, detJ_w_data) +end + +# ============================================================================ +# INTEGRATION WRAPPERS (material-specific) +# ============================================================================ + +""" + compute_block!( + prepared::PreparedElement, + material::LinearElastic, + k_local::Int, + l_local::Int + ) -> Tensor{2,3} + +Integrate weak form over element to get stiffness block between nodes k and l. + +For LinearElastic, the tangent modulus 𝔻 is constant (equal to C), so we +compute it once via compute_stress and integrate using compute_block_at_point. + +# Returns +Fully integrated 3×3 stiffness block K[k,l] +""" +@inline function compute_block!( + prepared::PreparedElement{N,NIP}, + material::LinearElastic, + k_local::Int, + l_local::Int +) where {N,NIP} + + # For LinearElastic, tangent modulus is constant (𝔻 = C) + # Compute at reference strain E = 0 + E_ref = zero(SymmetricTensor{2,3,Float64}) + S, 𝔻, _ = compute_stress(material, E_ref) + + K_kl = zero(Tensor{2,3,Float64}) + + # Integrate over all quadrature points + @inbounds for q in 1:NIP + grad_k = prepared.∇N_data[q][k_local] + grad_l = prepared.∇N_data[q][l_local] + + # Weak form contribution at this point using material tangent + K_kl_ip = compute_block_at_point(grad_k, grad_l, 𝔻) + + # Accumulate with integration weight + K_kl += K_kl_ip * prepared.detJ_w[q] + end + + return K_kl +end + +""" + compute_block!( + prepared::PreparedElement, + material::NeoHookean, + k_local::Int, + l_local::Int, + u_elem::AbstractVector{Float64} + ) -> Tensor{2,3} + +Integrate weak form for NeoHookean material (strain-dependent tangent). + +For nonlinear materials, the material tensor depends on strain, so we must +compute it at each integration point using the current displacement field. + +# Arguments +- `u_elem`: Element displacement DOFs [3N] (for computing deformation gradient) +""" +@inline function compute_block!( + prepared::PreparedElement{N,NIP}, + material::NeoHookean, + k_local::Int, + l_local::Int, + u_elem::AbstractVector{Float64} +) where {N,NIP} + + I = one(Tensor{2,3,Float64}) + K_kl = zero(Tensor{2,3,Float64}) + + @inbounds for q in 1:NIP + ∇N_q = prepared.∇N_data[q] + + # Compute deformation gradient F at this integration point + F = I + for k in 1:N + k_offset = 3(k - 1) + u_k = Vec{3}((u_elem[k_offset+1], u_elem[k_offset+2], u_elem[k_offset+3])) + F += u_k ⊗ ∇N_q[k] + end + + # Right Cauchy-Green and Green-Lagrange strain + C_tensor = symmetric(F' ⋅ F) + E = SymmetricTensor{2,3}(0.5 * (C_tensor - I)) + + # Get material tangent modulus (strain-dependent) + _, 𝔻, _ = compute_stress(material, E) + + # Weak form contribution at this point with current tangent + grad_k = ∇N_q[k_local] + grad_l = ∇N_q[l_local] + K_kl_ip = compute_block_at_point(grad_k, grad_l, 𝔻) + + K_kl += K_kl_ip * prepared.detJ_w[q] + end + + return K_kl +end + +# ============================================================================ +# BACKWARD COMPATIBILITY (element-based assemblers) +# ============================================================================ + +""" + compute_all_blocks!( + K_blocks::AbstractMatrix{Tensor{2,3}}, + prepared::PreparedElement, + material::LinearElastic, + u_elem, + N::Int + ) + +Compute all N×N blocks for LinearElastic material. + +Helper for element-based assemblers. For nodal assemblers, call compute_block! +directly for only the needed blocks. +""" +@inline function compute_all_blocks!( + K_blocks::AbstractMatrix{<:Tensor{2,3}}, + prepared::PreparedElement{N}, + material::LinearElastic, + u_elem, + Nnodes::Int +) where {N} + @inbounds for k in 1:Nnodes, l in 1:Nnodes + K_blocks[k, l] = compute_block!(prepared, material, k, l) + end +end + +""" + compute_all_blocks!( + K_blocks::AbstractMatrix{Tensor{2,3}}, + prepared::PreparedElement, + material::NeoHookean, + u_elem::AbstractVector{Float64}, + N::Int + ) + +Compute all N×N blocks for NeoHookean material (requires displacement field). +""" +@inline function compute_all_blocks!( + K_blocks::AbstractMatrix{<:Tensor{2,3}}, + prepared::PreparedElement{N}, + material::NeoHookean, + u_elem::AbstractVector{Float64}, + Nnodes::Int +) where {N} + @inbounds for k in 1:Nnodes, l in 1:Nnodes + K_blocks[k, l] = compute_block!(prepared, material, k, l, u_elem) + end +end + +""" + blocked_tensor_to_matrix_view!( + K_e::AbstractMatrix{Float64}, + K_blocks::AbstractMatrix{Tensor{2,3}} + ) + +Convert N×N matrix of 3×3 tensor blocks to 3N×3N Float64 matrix. + +Maps block[k,l][α,β] → K_e[3(k-1)+α, 3(l-1)+β] +""" +function blocked_tensor_to_matrix_view!( + K_e::AbstractMatrix{Float64}, + K_blocks::AbstractMatrix{<:Tensor{2,3}} +) + Nnodes = size(K_blocks, 1) + @inbounds for k in 1:Nnodes, l in 1:Nnodes + block = K_blocks[k, l] + k_offset = 3(k - 1) + l_offset = 3(l - 1) + for α in 1:3, β in 1:3 + K_e[k_offset + α, l_offset + β] = block[α, β] + end + end +end + +""" + compute_element_stiffness!( + cache::ElementCache, + kernel::ContinuumKernel, + element_id::Int, + mesh::AbstractMesh + ) + +Compute element stiffness matrix and force vector (writes to cache). + +Implements the kernel interface for element-based assemblers (COO, CSC). +Writes results to `cache.Ke` and `cache.fe` without allocating. + +For nodal assemblers, use prepare_element! + compute_block! directly. + +# Arguments +- `cache`: Element cache with pre-allocated Ke, fe, K_blocks +- `kernel`: Continuum kernel with material +- `element_id`: Element index in mesh +- `mesh`: Finite element mesh + +# Side Effects +Writes to: +- `cache.Ke` - Element stiffness matrix [3N × 3N] +- `cache.fe` - Element force vector [3N] (zeros for LinearElastic) +""" +function compute_element_stiffness!( + cache::ElementCache{T,B,IPS}, + kernel::ContinuumKernel{M}, + element_id::Int, + mesh::AbstractMesh +) where {T<:AbstractTopology{N},B,IPS,M} where {N} + + # Zero outputs + fill!(cache.Ke, 0.0) + fill!(cache.fe, 0.0) + + # Prepare element geometry (uses cache.X_buffer for coordinates) + prepared = prepare_element!(cache, kernel, element_id, mesh) + + # Compute all blocks into cache.K_blocks (reuses existing allocation) + compute_all_blocks!(cache.K_blocks, prepared, kernel.material, cache.u_buffer, N) + + # Convert blocks to Float64 matrix in cache.Ke + blocked_tensor_to_matrix_view!(cache.Ke, cache.K_blocks) + + return nothing +end