feat(continuum): Implement Jacobian and integration utilities

- 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
2026-08-31 08:16:23 +00:00 · 2025-11-19 09:03:28 +02:00
parent 5ce552c962
commit fb10cbbf8a
1 changed files with 363 additions and 0 deletions
@@ -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