From ea1e6de8e915c2171b4cb5b12732f833c624bbb1 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Sat, 9 May 2026 18:12:50 +0300 Subject: [PATCH] chore(physics): delete unused assembly helper stubs Remove legacy glue that lived outside the active `src/physics/` include list. - Drop `src/physics/assembly_helpers.jl`. --- src/physics/assembly_helpers.jl | 331 -------------------------------- 1 file changed, 331 deletions(-) delete mode 100644 src/physics/assembly_helpers.jl diff --git a/src/physics/assembly_helpers.jl b/src/physics/assembly_helpers.jl deleted file mode 100644 index e3809c8..0000000 --- a/src/physics/assembly_helpers.jl +++ /dev/null @@ -1,331 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -using Tensors - -""" -Assembly helper functions for elasticity using Tensors.jl. - -All functions designed for: -- Zero allocations (stack-allocated tuples and tensors) -- Type stability (concrete types throughout) -- GPU compatibility (kernel-friendly operations) -- Compiler optimization (small loops unrolled automatically) - -See `docs/book/material_modeling.md` for design rationale. -""" - -""" - shape_function_gradients(element::Element, ip::IntegrationPoint) -> NTuple{N, Vec{3, Float64}} - -Compute shape function gradients in current configuration. - -# Returns - -Tuple of gradients (zero allocation!): -```julia -∇N = (∇N₁, ∇N₂, ..., ∇Nₙ) -where ∇Nᵢ::Vec{3, Float64} -``` - -# Implementation Note - -This is a placeholder. Full implementation requires: -1. Evaluate basis in reference config: ∇N_ref -2. Compute Jacobian: J = ∂X/∂ξ = ∑ᵢ Xᵢ ⊗ ∇Nᵢ_ref -3. Transform to current config: ∇N = J⁻ᵀ · ∇N_ref - -Current code uses existing BasisInfo infrastructure. -""" -function shape_function_gradients( - bi::BasisInfo{B,T}, - X::NTuple{N,Vec{3,T}}, - ip -) where {B,T,N} - # Evaluate basis (fills bi.grad with ∂N/∂X) - eval_basis!(bi, X, ip) - - # Convert Matrix to NTuple{N, Vec{3}} - # bi.grad is (3, N) matrix - grads = ntuple(N) do i - Vec{3}(bi.grad[1, i], bi.grad[2, i], bi.grad[3, i]) - end - - return grads -end - -""" - compute_strain_from_gradients(∇N::NTuple{N, Vec{3}}, u::Vector{Float64}) -> SymmetricTensor{2,3} - -Compute strain tensor from shape function gradients and displacement. - -# Small Strain (Linear) - -``` -ε = sym(∇u) = ½(∇u + ∇uᵀ) -``` - -where ∇u = ∑ᵢ uᵢ ⊗ ∇Nᵢ - -# Arguments - -- `∇N`: Tuple of shape function gradients (from `shape_function_gradients`) -- `u`: Nodal displacement vector [u₁ₓ, u₁ᵧ, u₁ᵤ, u₂ₓ, ...] - -# Returns - -- `ε`: Symmetric strain tensor (6 unique components, stack-allocated) - -# Performance - -- **Time:** ~10 ns (validated in benchmarks) -- **Allocations:** 0 bytes (stack-allocated) -- **Type stability:** ✅ (concrete return type) - -# Example - -```julia -∇N = (Vec(0.1, 0.0, 0.0), Vec(0.0, 0.1, 0.0), ...) -u = [0.01, 0.02, 0.00, ...] # Nodal displacements - -ε = compute_strain_from_gradients(∇N, u) -# Returns: SymmetricTensor{2,3}([ε₁₁, ε₂₂, ε₃₃, ε₁₂, ε₂₃, ε₁₃]) -``` -""" -function compute_strain_from_gradients( - ∇N::NTuple{N,Vec{3,T}}, - u::Vector{T} -) where {N,T} - # Deformation gradient: F = I + ∇u = I + ∑ᵢ uᵢ ⊗ ∇Nᵢ - F = one(Tensor{2,3,T}) - - @inbounds for (i, ∇Nᵢ) in enumerate(∇N) - i_offset = 3(i - 1) - uᵢ = Vec{3}(u[i_offset+1], u[i_offset+2], u[i_offset+3]) - F += uᵢ ⊗ ∇Nᵢ - end - - # Small strain: ε = sym(∇u) = sym(F - I) - ε = symmetric(F) - one(SymmetricTensor{2,3,T}) - - return ε -end - -""" - compute_green_lagrange_strain(∇N::NTuple{N, Vec{3}}, u::Vector) -> SymmetricTensor{2,3} - -Compute Green-Lagrange strain for finite deformation. - -# Finite Strain (Nonlinear) - -``` -E = ½(∇u + ∇uᵀ + ∇uᵀ∇u) = ½(Fᵀ·F - I) = ½(C - I) -``` - -where: -- F = I + ∇u (deformation gradient) -- C = Fᵀ·F (right Cauchy-Green tensor) - -# Use When - -- `physics.finite_strain = true` -- Large deformations (>5% strain typically) -- Geometric nonlinearity important - -# Performance - -- **Time:** ~15 ns (slightly more than small strain) -- **Allocations:** 0 bytes -- **Type stability:** ✅ -""" -function compute_green_lagrange_strain( - ∇N::NTuple{N,Vec{3,T}}, - u::Vector{T} -) where {N,T} - # F = I + ∇u - F = one(Tensor{2,3,T}) - - @inbounds for (i, ∇Nᵢ) in enumerate(∇N) - i_offset = 3(i - 1) - uᵢ = Vec{3}(u[i_offset+1], u[i_offset+2], u[i_offset+3]) - F += uᵢ ⊗ ∇Nᵢ - end - - # E = ½(Fᵀ·F - I) = ½(C - I) - C = tdot(F) # Right Cauchy-Green: C = Fᵀ·F - E = T(0.5) * (C - one(C)) - - return E -end - -""" - accumulate_stiffness!(K_e, ∇N, 𝔻, w) -> K_e - -Accumulate stiffness contribution for integration point. - -# Formula - -For each node pair (i,j), accumulates 3×3 block: - -``` -K[i,j]ₐᵦ += w · ∑ₖₗ (∂Nᵢ/∂xₖ) · 𝔻ₐₖᵦₗ · (∂Nⱼ/∂xₗ) -``` - -# Loop Structure - -Three nested loops: -1. Node pairs (i,j) - 100 iterations for Tet10 -2. Spatial dimensions (a,b) - 9 iterations -3. Contraction (k,l) - 9 iterations - -Inner loops (a,b,k,l) are unrolled by compiler with `@inbounds @simd`. - -# Performance - -Per integration point (Tet10): -- **Node pair loops:** ~100 ns (10×10 nodes) -- **Per 3×3 block:** ~1 ns (compiler unrolls inner loops) -- **Total:** ~100 ns per IP - -# Arguments - -- `K_e`: Element stiffness matrix (ndofs × ndofs), modified in-place -- `∇N`: Shape function gradients (tuple from `shape_function_gradients`) -- `𝔻`: Material tangent modulus (SymmetricTensor{4,3} from material model) -- `w`: Integration weight × Jacobian determinant - -# Returns - -- `K_e` (for chaining, though modified in-place) -""" -function accumulate_stiffness!( - K_e::Matrix{T}, - ∇N::NTuple{N,Vec{3,T}}, - 𝔻::SymmetricTensor{4,3,T}, - w::T -) where {N,T} - - @inbounds for (i, ∇Nᵢ) in enumerate(∇N) - i_offset = 3(i - 1) - - for (j, ∇Nⱼ) in enumerate(∇N) - j_offset = 3(j - 1) - - # Each (i,j): 3×3 block - @inbounds for a in 1:3, b in 1:3 - Kval = zero(T) - @simd for k in 1:3, l in 1:3 - Kval += ∇Nᵢ[k] * 𝔻[a, k, b, l] * ∇Nⱼ[l] - end - K_e[i_offset+a, j_offset+b] += w * Kval - end - end - end - - return K_e -end - -""" - accumulate_internal_forces!(f_int, ∇N, σ, w) -> f_int - -Accumulate internal force contribution for integration point. - -# Formula - -For each node i: - -``` -fᵢ = w · (σ · ∇Nᵢ) -``` - -where: -- σ is Cauchy stress tensor -- ∇Nᵢ is shape function gradient -- w is integration weight - -# Performance - -- **Time:** ~50 ns per IP (10 nodes × 5 ns per node) -- **Allocations:** 0 bytes -- **Type stability:** ✅ - -# Arguments - -- `f_int`: Internal force vector (ndofs), modified in-place -- `∇N`: Shape function gradients -- `σ`: Cauchy stress tensor (from material model) -- `w`: Integration weight × Jacobian - -# Returns - -- `f_int` (for chaining, though modified in-place) -""" -function accumulate_internal_forces!( - f_int::Vector{T}, - ∇N::NTuple{N,Vec{3,T}}, - σ::SymmetricTensor{2,3,T}, - w::T -) where {N,T} - - @inbounds for (i, ∇Nᵢ) in enumerate(∇N) - i_offset = 3(i - 1) - - # fᵢ = w · (σ · ∇Nᵢ) - # Use double contraction: σ ⊡ ∇Nᵢ - f_i = w * (σ ⊡ ∇Nᵢ) - - for a in 1:3 - f_int[i_offset+a] += f_i[a] - end - end - - return f_int -end - -""" - accumulate_external_forces!(f_ext, N, b, w) -> f_ext - -Accumulate external body force contribution. - -# Formula - -``` -fᵢ = w · Nᵢ · b -``` - -where: -- b is body force vector [bₓ, bᵧ, bᵤ] -- Nᵢ is shape function value -- w is integration weight - -# Arguments - -- `f_ext`: External force vector -- `N`: Shape function values (ntuple) -- `b`: Body force vector (Vec{3}) -- `w`: Integration weight - -# Returns - -- `f_ext` (modified in-place) -""" -function accumulate_external_forces!( - f_ext::Vector{T}, - N::NTuple{N_nodes,T}, - b::Vec{3,T}, - w::T -) where {N_nodes,T} - - @inbounds for (i, Nᵢ) in enumerate(N) - i_offset = 3(i - 1) - - # fᵢ = w · Nᵢ · b - contribution = w * Nᵢ - - for a in 1:3 - f_ext[i_offset+a] += contribution * b[a] - end - end - - return f_ext -end