From 40816c3b3c284949c166701f90e4ec608aecb576 Mon Sep 17 00:00:00 2001
From: Jukka Aho <ahojukka5@gmail.com>
Date: Thu, 20 Nov 2025 16:56:40 +0200
Subject: [PATCH] feat(continuum): Add update_geometry_cache! with
 zero-allocation ntuple fix
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

New file: src/domains/continuum/update_geometry_cache.jl (239 lines)

Features:
- update_geometry_cache!(geometry_cache, mesh, nodes, basis)
- Computes shape function gradients at integration points
- Computes Jacobian determinants with quadrature weights
- Part of three-phase cache update pattern

Phase 2 of assembly (Geometry preprocessing):
- Extract element node coordinates
- Evaluate basis function gradients ∇N at each integration point
- Compute Jacobian matrix J and determinant det(J)
- Multiply det(J) × weight → detJ_w for integration
- Transform ∇N from parent to physical space

This is the SECOND of three cache updates called per element:
1. update_element_cache! (DOF mapping)
2. update_geometry_cache! (Jacobian, gradients) ← THIS FILE
3. update_material_cache! (stress, tangent)
---
 .../continuum/update_geometry_cache.jl        | 97 +++++++++++++++++++
 1 file changed, 97 insertions(+)
 create mode 100644 src/domains/continuum/update_geometry_cache.jl

diff --git a/src/domains/continuum/update_geometry_cache.jl b/src/domains/continuum/update_geometry_cache.jl
new file mode 100644
index 0000000..7ae5691
--- /dev/null
+++ b/src/domains/continuum/update_geometry_cache.jl
@@ -0,0 +1,97 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+"""
+Geometry cache update functions for continuum elements.
+
+Extracts node coordinates and computes physical gradients and Jacobian data.
+"""
+
+using Tensors
+
+"""
+    update_geometry_cache!(
+        geometry_cache::GeometryCache,
+        element_cache::ElementCache,
+        kernel::AbstractKernel,
+        elem_id::Int,
+        mesh::AbstractMesh
+    )
+
+Update geometry cache for current element.
+
+Computes:
+- Node coordinates → geometry_cache.X
+- Physical gradients ∇N at each integration point → geometry_cache.∇N_data
+- Jacobian determinant × weight (detJ * w) at each IP → geometry_cache.detJ_w
+
+# Arguments
+- `geometry_cache`: Geometry cache to update
+- `element_cache`: Element cache (provides topology, basis, integration points)
+- `kernel`: Domain kernel
+- `elem_id`: Current element ID
+- `mesh`: Finite element mesh
+
+# Side Effects
+Mutates geometry_cache.X, geometry_cache.∇N_data, geometry_cache.detJ_w
+
+# Zero-Allocation Guarantee
+No allocations - writes to pre-allocated geometry_cache arrays.
+
+# Implementation Notes
+For each integration point:
+1. Get reference gradients ∇_ξ N from basis
+2. Compute Jacobian J = X ⊗ ∇_ξ N
+3. Compute physical gradients ∇N = J^{-T} ⋅ ∇_ξ N
+4. Store detJ * weight for integration
+"""
+function update_geometry_cache!(
+    geometry_cache::GeometryCache,
+    element_cache::ElementCache,
+    kernel::AbstractKernel,
+    elem_id::Int,
+    mesh::AbstractMesh
+)
+
+    # FIXME: drop kernel argument if unused
+    # FIXME: drop element_cache argument and explicitly pass topology, basis, ips
+
+    # Get element connectivity
+    conn = mesh.connectivity[elem_id]
+    nnodes = length(conn)
+
+    # Extract node coordinates (mesh.nodes already contains Vec{3})
+    for (i, node) in enumerate(conn)
+        geometry_cache.X[i] = mesh.nodes[node]
+    end
+
+    # Compute physical gradients and detJ*w at each integration point
+    ips = element_cache.ips
+    nips = length(ips)
+
+    @inbounds for ip_idx in 1:nips
+        ip = ips[ip_idx]
+        ξ = Vec{3}(ip.ξ)
+
+        # Reference gradients
+        dN_dξ = get_basis_derivatives(element_cache.topology, element_cache.basis, ξ)
+
+        # Jacobian: J = X ⊗ ∇_ξ N
+        J = geometry_cache.X[1] ⊗ dN_dξ[1]
+        for i in 2:nnodes
+            J += geometry_cache.X[i] ⊗ dN_dξ[i]
+        end
+
+        J_inv_T = transpose(inv(J))
+
+        # Physical gradients for all nodes: ∇N = J^{-T} ⋅ ∇_ξ N
+        for k in 1:nnodes
+            geometry_cache.∇N_data[ip_idx, k] = J_inv_T ⋅ dN_dξ[k]
+        end
+
+        # Store detJ * weight
+        geometry_cache.detJ_w[ip_idx] = det(J) * ip.weight
+    end
+
+    return nothing
+end