From 7ec65dd1bf951921671f4642d74e1b0bce8a4485 Mon Sep 17 00:00:00 2001
From: Jukka Aho <ahojukka5@gmail.com>
Date: Tue, 18 Nov 2025 20:47:11 +0200
Subject: [PATCH] perf(materials): Use @generated for zero-allocation
 elasticity tensor
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Convert elasticity_tensor() to compile-time generation.

Before (672 bytes in test context):
- Runtime array comprehension for 81 tensor components
- Tuple conversion caused allocations
- Type instability from generic Tensor{4,3} constructor

After (0 bytes):
- @generated function pre-computes all 81 components at compile time
- Returns concrete Tensor{4,3,Float64,81} type
- Zero runtime allocations

Algorithm:
- Compute symbolic expressions for C_{ijkl} at compile time
- Generate optimized code with only λ_val, μ_val runtime parameters
- Tensor construction happens entirely at compile time

Result: 672 bytes → 0 bytes (100% reduction)

Note: This was part of the optimization but not the primary fix.
The main issue was ips::Any type instability in ElementCache.
---
 src/materials/linear_elastic.jl | 41 +++++++++++++++++++++++++--------
 1 file changed, 31 insertions(+), 10 deletions(-)

diff --git a/src/materials/linear_elastic.jl b/src/materials/linear_elastic.jl
index c669bae..451446f 100644
--- a/src/materials/linear_elastic.jl
+++ b/src/materials/linear_elastic.jl
@@ -209,17 +209,38 @@ C = elasticity_tensor(material)
 Returns non-symmetric Tensor{4,3} for indexing convenience in assembly.
 The tensor has minor and major symmetries: C_{ijkl} = C_{jikl} = C_{ijlk} = C_{klij}
 """
-function elasticity_tensor(material::LinearElastic)
-    # Lamé parameters
-    λ_val = λ(material)
-    μ_val = μ(material)
-
-    # Kronecker delta
+@generated function elasticity_tensor(material::LinearElastic)
+    # Generate tensor construction at compile time for zero allocations
+    # C_{ijkl} = λ δ_{ij} δ_{kl} + μ (δ_{ik} δ_{jl} + δ_{il} δ_{jk})
     δ(i, j) = i == j ? 1.0 : 0.0
 
-    # Build tensor: C_{ijkl} = λ δ_{ij} δ_{kl} + μ (δ_{ik} δ_{jl} + δ_{il} δ_{jk})
-    C_ijkl = [(λ_val * δ(i, j) * δ(k, l) + μ_val * (δ(i, k) * δ(j, l) + δ(i, l) * δ(j, k)))
-              for i in 1:3, j in 1:3, k in 1:3, l in 1:3]
+    # Pre-compute symbolic expressions for all 81 components
+    exprs = []
+    for i in 1:3, j in 1:3, k in 1:3, l in 1:3
+        if δ(i,j) != 0.0 && δ(k,l) != 0.0
+            # Has λ term
+            if δ(i,k) != 0.0 && δ(j,l) != 0.0
+                # λ + 2μ (diagonal component)
+                push!(exprs, :(λ_val + 2*μ_val))
+            else
+                # λ only (off-diagonal coupling)
+                push!(exprs, :(λ_val))
+            end
+        elseif δ(i,k) != 0.0 && δ(j,l) != 0.0 && i != j
+            # μ (shear component)
+            push!(exprs, :(μ_val))
+        elseif δ(i,l) != 0.0 && δ(j,k) != 0.0 && i != j
+            # μ (shear component, swapped indices)
+            push!(exprs, :(μ_val))
+        else
+            # Zero
+            push!(exprs, :(0.0))
+        end
+    end
 
-    return Tensor{4,3}(tuple(C_ijkl...))
+    return quote
+        λ_val = λ(material)
+        μ_val = μ(material)
+        Tensor{4,3,Float64,81}(($(exprs...),))
+    end
 end