From 70527936da46dc99287c53f4687989dd2fe42d86 Mon Sep 17 00:00:00 2001
From: miguelmaso <miguel.maso@upc.edu>
Date: Fri, 11 Sep 2026 13:46:57 +0200
Subject: [PATCH] Zero-allocation tangent for ViscousPolyconvex
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Add generated implementations of `-` (binary and unary) for `TensorValue`,
mirroring the existing `+`. Without them, subtraction of fourth-order tensors
fell back to the generic Gridap implementation, which exceeds Julia's inlining
threshold for `TensorValue{9,9}` and heap-allocates the result.

Those subtractions appear in `∂Sv∂C_Cᵥfix` and `∂Cv⁻¹∂C`, i.e. exactly on the
path of the viscoelastic tangent operator, which is why only `∂∂Ψ∂FF` allocated
while `Ψ` and `∂Ψ∂F` did not.

`+` is also generalized from `TensorValue{D,D}` to `TensorValue{D1,D2}` so that
non-square tensors are covered as well. No change in the mathematics.

Visco-polyconvex ∂∂Ψ∂FF: 187 allocs / 15552 B / 5367 ns -> 0 allocs / 0 B / 599 ns.
Visco-elastic ∂∂Ψ∂FF: 1021 allocs / 49600 ns -> 501 allocs / 34200 ns.

Also interpolate the arguments of the tensor algebra benchmarks. Without `$`,
the benchmarked expressions read non-const globals, so the call is dynamically
dispatched and the result is boxed: this is the source of the spurious "1 alloc
of 0.656 kB" reported for `IIsym`, `push_forward_C_to_F` and `×ᵢ⁴`, which are
already allocation-free.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
---
 .../TensorAlgebraBenchmarks.jl                | 10 ++++----
 src/TensorAlgebra/Operations.jl               | 24 ++++++++++++++++---
 src/TensorAlgebra/TensorAlgebra.jl            |  2 ++
 3 files changed, 28 insertions(+), 8 deletions(-)

diff --git a/benchmark/TensorAlgebraBenchmarks/TensorAlgebraBenchmarks.jl b/benchmark/TensorAlgebraBenchmarks/TensorAlgebraBenchmarks.jl
index d7974b2..9334a52 100644
--- a/benchmark/TensorAlgebraBenchmarks/TensorAlgebraBenchmarks.jl
+++ b/benchmark/TensorAlgebraBenchmarks/TensorAlgebraBenchmarks.jl
@@ -49,8 +49,8 @@ H = TensorValue(1.:81...)
 
 SUITE["Tensor algebra"]["δδ_μ_2d"] = @benchmarkable δᵢₖδⱼₗ2D + δᵢₗδⱼₖ2D
 SUITE["Tensor algebra"]["δδ_λ_2d"] = @benchmarkable 1.0 * δᵢⱼδₖₗ2D
-SUITE["Tensor algebra"]["Cofactor"] = @benchmarkable cof(A)
-SUITE["Tensor algebra"]["Det(A)Inv(A')"] = @benchmarkable det(A)*inv(A')
-SUITE["Tensor algebra"]["×ᵢ⁴"] = @benchmarkable ×ᵢ⁴(A)
-SUITE["Tensor algebra"]["IIsym"] = @benchmarkable IIsym(A)
-SUITE["Tensor algebra"]["push_forward_C_to_F"] = @benchmarkable push_forward_C_to_F(F,H)
+SUITE["Tensor algebra"]["Cofactor"] = @benchmarkable cof($A)
+SUITE["Tensor algebra"]["Det(A)Inv(A')"] = @benchmarkable det($A)*inv($A')
+SUITE["Tensor algebra"]["×ᵢ⁴"] = @benchmarkable ×ᵢ⁴($A)
+SUITE["Tensor algebra"]["IIsym"] = @benchmarkable IIsym($A)
+SUITE["Tensor algebra"]["push_forward_C_to_F"] = @benchmarkable push_forward_C_to_F($F,$H)
diff --git a/src/TensorAlgebra/Operations.jl b/src/TensorAlgebra/Operations.jl
index 24a974b..985f3f6 100644
--- a/src/TensorAlgebra/Operations.jl
+++ b/src/TensorAlgebra/Operations.jl
@@ -10,12 +10,30 @@ function (*)(Ten1::TensorValue, Ten2::TensorValue)
 end
 
 
-@inline @generated function (+)(A::TensorValue{D,D}, B::TensorValue{D,D}) where {D}
+@inline @generated function (+)(A::TensorValue{D1,D2}, B::TensorValue{D1,D2}) where {D1,D2}
   str = ""
-  for i in 1:D*D
+  for i in 1:D1*D2
     str *= "A.data[$i] + B.data[$i], "
   end
-  Meta.parse("TensorValue{D,D}($str)")
+  Meta.parse("TensorValue{D1,D2}($str)")
+end
+
+
+@inline @generated function (-)(A::TensorValue{D1,D2}, B::TensorValue{D1,D2}) where {D1,D2}
+  str = ""
+  for i in 1:D1*D2
+    str *= "A.data[$i] - B.data[$i], "
+  end
+  Meta.parse("TensorValue{D1,D2}($str)")
+end
+
+
+@inline @generated function (-)(A::TensorValue{D1,D2}) where {D1,D2}
+  str = ""
+  for i in 1:D1*D2
+    str *= "-A.data[$i], "
+  end
+  Meta.parse("TensorValue{D1,D2}($str)")
 end
 
 
diff --git a/src/TensorAlgebra/TensorAlgebra.jl b/src/TensorAlgebra/TensorAlgebra.jl
index dfc7839..76abf9a 100644
--- a/src/TensorAlgebra/TensorAlgebra.jl
+++ b/src/TensorAlgebra/TensorAlgebra.jl
@@ -6,9 +6,11 @@ using StaticArrays
 using LinearAlgebra
 import Base: *
 import Base: +
+import Base: -
 
 export (*)
 export (+)
+export (-)
 export (⊗₁₂³)
 export (⊗₁₃²)
 export (⊗₁²³)