From ad934e37f6746dd786d2cd9a9bec6092a328234f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Miguel=20Mas=C3=B3?= <miguel.maso@upc.edu>
Date: Sat, 13 Sep 2025 18:31:34 +0200
Subject: [PATCH 01/13] added outer products generated

---
 src/TensorAlgebra/TensorAlgebra.jl | 161 ++++++++++-------------------
 1 file changed, 57 insertions(+), 104 deletions(-)

diff --git a/src/TensorAlgebra/TensorAlgebra.jl b/src/TensorAlgebra/TensorAlgebra.jl
index be3977a..ddc88e3 100644
--- a/src/TensorAlgebra/TensorAlgebra.jl
+++ b/src/TensorAlgebra/TensorAlgebra.jl
@@ -14,8 +14,9 @@ export (+)
 export (⊗₁₂³)
 export (⊗₁₃²)
 export (⊗₁²³)
-export (⊗₁₃²⁴)
 export (⊗₁₂³⁴)
+export (⊗₁₃²⁴)
+export (⊗₁₄²³)
 export (⊗₁²)
 export I3
 export I9
@@ -205,111 +206,63 @@ function (⊗₁₃²)(A::TensorValue{2,2,Float64}, V::VectorValue{2,Float64})
     A[4] * V[2])
 end
 
-function (⊗₁₃²⁴)(A::TensorValue{3,3,Float64}, B::TensorValue{3,3,Float64})
-
-  TensorValue{9,9,Float64,81}(A[1] * B[1],
-    A[2] * B[1],
-    A[3] * B[1],
-    A[1] * B[2],
-    A[2] * B[2],
-    A[3] * B[2],
-    A[1] * B[3],
-    A[2] * B[3],
-    A[3] * B[3],
-    A[4] * B[1],
-    A[5] * B[1],
-    A[6] * B[1],
-    A[4] * B[2],
-    A[5] * B[2],
-    A[6] * B[2],
-    A[4] * B[3],
-    A[5] * B[3],
-    A[6] * B[3],
-    A[7] * B[1],
-    A[8] * B[1],
-    A[9] * B[1],
-    A[7] * B[2],
-    A[8] * B[2],
-    A[9] * B[2],
-    A[7] * B[3],
-    A[8] * B[3],
-    A[9] * B[3],
-    A[1] * B[4],
-    A[2] * B[4],
-    A[3] * B[4],
-    A[1] * B[5],
-    A[2] * B[5],
-    A[3] * B[5],
-    A[1] * B[6],
-    A[2] * B[6],
-    A[3] * B[6],
-    A[4] * B[4],
-    A[5] * B[4],
-    A[6] * B[4],
-    A[4] * B[5],
-    A[5] * B[5],
-    A[6] * B[5],
-    A[4] * B[6],
-    A[5] * B[6],
-    A[6] * B[6],
-    A[7] * B[4],
-    A[8] * B[4],
-    A[9] * B[4],
-    A[7] * B[5],
-    A[8] * B[5],
-    A[9] * B[5],
-    A[7] * B[6],
-    A[8] * B[6],
-    A[9] * B[6],
-    A[1] * B[7],
-    A[2] * B[7],
-    A[3] * B[7],
-    A[1] * B[8],
-    A[2] * B[8],
-    A[3] * B[8],
-    A[1] * B[9],
-    A[2] * B[9],
-    A[3] * B[9],
-    A[4] * B[7],
-    A[5] * B[7],
-    A[6] * B[7],
-    A[4] * B[8],
-    A[5] * B[8],
-    A[6] * B[8],
-    A[4] * B[9],
-    A[5] * B[9],
-    A[6] * B[9],
-    A[7] * B[7],
-    A[8] * B[7],
-    A[9] * B[7],
-    A[7] * B[8],
-    A[8] * B[8],
-    A[9] * B[8],
-    A[7] * B[9],
-    A[8] * B[9],
-    A[9] * B[9])
-end
-
-function (⊗₁₃²⁴)(A::TensorValue{2,2,Float64}, B::TensorValue{2,2,Float64})
-
-  TensorValue{4,4,Float64,16}(A[1] * B[1],
-    A[2] * B[1],
-    A[1] * B[2],
-    A[2] * B[2],
-    A[3] * B[1],
-    A[4] * B[1],
-    A[3] * B[2],
-    A[4] * B[2],
-    A[1] * B[3],
-    A[2] * B[3],
-    A[1] * B[4],
-    A[2] * B[4],
-    A[3] * B[3],
-    A[4] * B[3],
-    A[3] * B[4],
-    A[4] * B[4])
+
+"""
+  ⊗₁₃²⁴(A::TensorValue, B::TensorValue)::TensorValue
+
+  Computes the **outer product** of two second-order tensors (matrices), returning a fourth-order tensor 
+  represented in a `D² x D²` flattened matrix using comined indices.
+
+  # Arguments
+  - `A::TensorValue{D}`: A second-order tensor.
+  - `B::TensorValue{D}`: Another second-order tensor.
+
+  # Returns
+  - `TensorValue{D*D}`: A statically-sized matrix representing the fourth-order tensor.
+"""
+@generated function (⊗₁₃²⁴)(A::TensorValue{D}, B::TensorValue{D}) where D
+  str = ""
+  for l in 1:D
+    for k in 1:D
+      for j in 1:D
+        for i in 1:D
+          str *= "A[$i,$k]*B[$j,$l],"
+        end
+      end
+    end
+  end
+  Meta.parse("TensorValue{D*D}($str)")
 end
 
+
+"""
+  ⊗₁₄²³(A::TensorValue, B::TensorValue)::TensorValue
+
+  Computes the **outer product** of two second-order tensors (matrices), returning a fourth-order tensor 
+  represented in a `D² x D²` flattened matrix using comined indices.
+
+  # Arguments
+  - `A::TensorValue{D}`: A second-order tensor.
+  - `B::TensorValue{D}`: Another second-order tensor.
+
+  # Returns
+  - `TensorValue{D*D}`: A statically-sized matrix representing the fourth-order tensor.
+"""
+@generated function (⊗₁₄²³)(A::TensorValue{D}, B::TensorValue{D}) where D
+  str = ""
+  for l in 1:D
+    for k in 1:D
+      for j in 1:D
+        for i in 1:D
+          str *= "A[$i,$l]*B[$j,$k],"
+        end
+      end
+    end
+  end
+  Meta.parse("TensorValue{D*D}($str)")
+end
+
+
 function (×ᵢ⁴)(A::TensorValue{3,3,Float64})
 
   TensorValue(0.0, 0.0, 0.0, 0.0, A[9], -A[8], 0.0, -A[6], A[5], 0.0, 0.0, 0.0, -A[9],

From 48eacd21c4a3b5c2ee4dfb6e557f2ea98fb9c3f0 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Miguel=20Mas=C3=B3?= <miguel.maso@upc.edu>
Date: Sat, 13 Sep 2025 18:31:51 +0200
Subject: [PATCH 02/13] testing outer products

---
 test/TestTensorAlgebra/TensorAlgebra.jl | 16 ++++++++++++++--
 1 file changed, 14 insertions(+), 2 deletions(-)

diff --git a/test/TestTensorAlgebra/TensorAlgebra.jl b/test/TestTensorAlgebra/TensorAlgebra.jl
index 08331f6..91d2237 100644
--- a/test/TestTensorAlgebra/TensorAlgebra.jl
+++ b/test/TestTensorAlgebra/TensorAlgebra.jl
@@ -1,5 +1,6 @@
 using Gridap.TensorValues
 using HyperFEM.TensorAlgebra
+using Test
 using BenchmarkTools
 
  
@@ -28,8 +29,19 @@ end
   @test norm(get_array(V1) ⊗ get_array(V2)) == 17.04406054905931
 
 end
- 
- 
+
+
+@testset "outer products" begin
+  A = TensorValue(1.0, 2.0, 3.0, 4.0)
+  B = TensorValue(5.0, 6.0, 7.0, 8.0)
+  u = VectorValue(1.0, 2.0)
+  v = VectorValue(3.0, 4.0)
+  @test A ⊗ B     == TensorValue(5.0, 10.0, 15.0, 20.0, 6.0, 12.0, 18.0, 24.0, 7.0, 14.0, 21.0, 28.0, 8.0, 16.0, 24.0, 32.0)
+  @test A ⊗₁₂³⁴ B == TensorValue(5.0, 10.0, 15.0, 20.0, 6.0, 12.0, 18.0, 24.0, 7.0, 14.0, 21.0, 28.0, 8.0, 16.0, 24.0, 32.0)
+  @test A ⊗₁₃²⁴ B == TensorValue(5.0, 10.0, 6.0, 12.0, 15.0, 20.0, 18.0, 24.0, 7.0, 14.0, 8.0, 16.0, 21.0, 28.0, 24.0, 32.0)
+  @test A ⊗₁₄²³ B == TensorValue(5.0, 10.0, 6.0, 12.0, 7.0, 14.0, 8.0, 16.0, 15.0, 20.0, 18.0, 24.0, 21.0, 28.0, 24.0, 32.0)
+end
+
 
 # @benchmark (A ⊗₁₃²⁴ B)
 # @benchmark (A ⊗₁₂³ V1)

From 26cae9296c8c85be63660378d06abfa557038a75 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Miguel=20Mas=C3=B3?= <miguel.maso@upc.edu>
Date: Sat, 13 Sep 2025 18:33:46 +0200
Subject: [PATCH 03/13] removed old outer

---
 src/TensorAlgebra/TensorAlgebra.jl | 58 ------------------------------
 1 file changed, 58 deletions(-)

diff --git a/src/TensorAlgebra/TensorAlgebra.jl b/src/TensorAlgebra/TensorAlgebra.jl
index ddc88e3..510cc74 100644
--- a/src/TensorAlgebra/TensorAlgebra.jl
+++ b/src/TensorAlgebra/TensorAlgebra.jl
@@ -32,8 +32,6 @@ export δᵢₖδⱼₗ3D
 export δᵢₗδⱼₖ3D
 export sqrt
 export cof
-export outer_13_24
-export outer_14_23
 export contraction_IP_JPKL
 export contraction_IP_PJKL
  
@@ -598,62 +596,6 @@ end
 end
 
 
-"""
-  outer_13_24(Matrix1::TensorValue, Matrix2::TensorValue)::TensorValue
-
-  Computes the **outer product** of two second-order tensors (matrices), returning a fourth-order tensor 
-  represented in a `D² x D²` matrix form (Voigt or flattened index notation) using combined indices.
-
-  # Arguments
-  - `Matrix1::TensorValue{D, D}`: A second-order tensor (e.g., a stress, strain, or deformation tensor).
-  - `Matrix2::TensorValue{D, D}`: Another second-order tensor to be used in the outer product.
-
-  # Returns
-  - `TensorValue{D^2, D^2}`: A statically-sized matrix representing the fourth-order tensor formed from the outer product.
-"""
-function outer_13_24(Matrix1::TensorValue{D}, Matrix2::TensorValue{D}) where D
-  Out = zeros(Float64,D*D,D*D)
-  for j in 1:D
-    for i in 1:D
-      for l in 1:D
-        for k in 1:D
-          Out[i + D * (j - 1), k + D * (l - 1)] += Matrix1[i,k] * Matrix2[j,l]
-        end
-      end
-    end
-  end
-  TensorValue{D*D,D*D}(Out)
-end
-
-
-"""
-  outer_14_23(Matrix1::TensorValue, Matrix2::TensorValue)::TensorValue
-
-  Computes the **outer product** of two second-order tensors (matrices), returning a fourth-order tensor 
-  represented in a `D² x D²` matrix form (flattened notation) using combined indices.
-
-  # Arguments
-  - `Matrix1::TensorValue{D, D}`: A second-order tensor (e.g., a stress, strain, or deformation tensor).
-  - `Matrix2::TensorValue{D, D}`: Another second-order tensor to be used in the outer product.
-
-  # Returns
-  - `TensorValue{D^2, D^2}`: A statically-sized matrix representing the fourth-order tensor formed from the outer product.
-"""
-function outer_14_23(Matrix1::TensorValue{D}, Matrix2::TensorValue{D}) where D
-  Out = zeros(Float64,D*D,D*D)
-  for j in 1:D
-    for i in 1:D
-      for l in 1:D
-        for k in 1:D
-          Out[i + D * (j - 1), k + D * (l - 1)] += Matrix1[i,l] * Matrix2[j,k]
-        end
-      end
-    end
-  end
-  TensorValue{D*D,D*D}(Out)
-end
-
-
 """
   contraction_IP_JPKL(Matrix1::TensorValue{D}, Matrix2::TensorValue{D*D})::TensorValue{D*D}
 

From c47545d8a009d35eb980808ba8d7cf9ff176ab5b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Miguel=20Mas=C3=B3?= <miguel.maso@upc.edu>
Date: Sat, 13 Sep 2025 18:37:19 +0200
Subject: [PATCH 04/13] docstrings

---
 src/TensorAlgebra/TensorAlgebra.jl | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/TensorAlgebra/TensorAlgebra.jl b/src/TensorAlgebra/TensorAlgebra.jl
index 510cc74..4104cd0 100644
--- a/src/TensorAlgebra/TensorAlgebra.jl
+++ b/src/TensorAlgebra/TensorAlgebra.jl
@@ -48,20 +48,20 @@ export Ellipsoid
 
 
 """
-  msqrt(A::TensorValue)::TensorValue
+  sqrt(A::TensorValue)::TensorValue
 
   Compute the square root of a 3x3 matrix by means of eigen decomposition.
 
   # Arguments
-  - `A::TensorValue`: the matrix to calculate the square root
+  - `A::TensorValue{3}`: the matrix to calculate the square root
 
   # Returns
-  - `::TensorValue`: the squared root matrix
+  - `::TensorValue{3}`: the squared root matrix
 """
-function sqrt(A::TensorValue{3,3})
+function sqrt(A::TensorValue{3})
   λ, Q = eigen(get_array(A))  # TODO: the get_array must be removed as long as it is supported after https://github.com/gridap/Gridap.jl/pull/1157
   λ = sqrt.(λ)
-  TensorValue{3,3}(λ[1]*Q[1:3]*Q[1:3]' + λ[2]*Q[4:6]*Q[4:6]' + λ[3]*Q[7:9]*Q[7:9]')
+  TensorValue{3}(λ[1]*Q[1:3]*Q[1:3]' + λ[2]*Q[4:6]*Q[4:6]' + λ[3]*Q[7:9]*Q[7:9]')
 end
 
 

From a941fc141bc5cb061cb4b8d1a4c621fea04725d6 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Miguel=20Mas=C3=B3?= <miguel.maso@upc.edu>
Date: Mon, 15 Sep 2025 12:16:11 +0200
Subject: [PATCH 05/13] move implementation

---
 src/TensorAlgebra/TensorAlgebra.jl | 112 ++++++++++++++---------------
 1 file changed, 56 insertions(+), 56 deletions(-)

diff --git a/src/TensorAlgebra/TensorAlgebra.jl b/src/TensorAlgebra/TensorAlgebra.jl
index 4104cd0..195c1cf 100644
--- a/src/TensorAlgebra/TensorAlgebra.jl
+++ b/src/TensorAlgebra/TensorAlgebra.jl
@@ -138,6 +138,62 @@ end
 end
 
 
+"""
+  ⊗₁₃²⁴(A::TensorValue, B::TensorValue)::TensorValue
+
+  Computes the **outer product** of two second-order tensors (matrices), returning a fourth-order tensor 
+  represented in a `D² x D²` flattened matrix using comined indices.
+
+  # Arguments
+  - `A::TensorValue{D}`: A second-order tensor.
+  - `B::TensorValue{D}`: Another second-order tensor.
+
+  # Returns
+  - `TensorValue{D*D}`: A statically-sized matrix representing the fourth-order tensor.
+"""
+@generated function (⊗₁₃²⁴)(A::TensorValue{D}, B::TensorValue{D}) where D
+  str = ""
+  for l in 1:D
+    for k in 1:D
+      for j in 1:D
+        for i in 1:D
+          str *= "A[$i,$k]*B[$j,$l],"
+        end
+      end
+    end
+  end
+  Meta.parse("TensorValue{D*D}($str)")
+end
+
+
+"""
+  ⊗₁₄²³(A::TensorValue, B::TensorValue)::TensorValue
+
+  Computes the **outer product** of two second-order tensors (matrices), returning a fourth-order tensor 
+  represented in a `D² x D²` flattened matrix using comined indices.
+
+  # Arguments
+  - `A::TensorValue{D}`: A second-order tensor.
+  - `B::TensorValue{D}`: Another second-order tensor.
+
+  # Returns
+  - `TensorValue{D*D}`: A statically-sized matrix representing the fourth-order tensor.
+"""
+@generated function (⊗₁₄²³)(A::TensorValue{D}, B::TensorValue{D}) where D
+  str = ""
+  for l in 1:D
+    for k in 1:D
+      for j in 1:D
+        for i in 1:D
+          str *= "A[$i,$l]*B[$j,$k],"
+        end
+      end
+    end
+  end
+  Meta.parse("TensorValue{D*D}($str)")
+end
+
+
 @generated function (⊗₁²³)(V::VectorValue{D,Float64}, A::TensorValue{D,D,Float64}) where {D}
   str = ""
   for iA in 1:D*D
@@ -205,62 +261,6 @@ function (⊗₁₃²)(A::TensorValue{2,2,Float64}, V::VectorValue{2,Float64})
 end
 
 
-"""
-  ⊗₁₃²⁴(A::TensorValue, B::TensorValue)::TensorValue
-
-  Computes the **outer product** of two second-order tensors (matrices), returning a fourth-order tensor 
-  represented in a `D² x D²` flattened matrix using comined indices.
-
-  # Arguments
-  - `A::TensorValue{D}`: A second-order tensor.
-  - `B::TensorValue{D}`: Another second-order tensor.
-
-  # Returns
-  - `TensorValue{D*D}`: A statically-sized matrix representing the fourth-order tensor.
-"""
-@generated function (⊗₁₃²⁴)(A::TensorValue{D}, B::TensorValue{D}) where D
-  str = ""
-  for l in 1:D
-    for k in 1:D
-      for j in 1:D
-        for i in 1:D
-          str *= "A[$i,$k]*B[$j,$l],"
-        end
-      end
-    end
-  end
-  Meta.parse("TensorValue{D*D}($str)")
-end
-
-
-"""
-  ⊗₁₄²³(A::TensorValue, B::TensorValue)::TensorValue
-
-  Computes the **outer product** of two second-order tensors (matrices), returning a fourth-order tensor 
-  represented in a `D² x D²` flattened matrix using comined indices.
-
-  # Arguments
-  - `A::TensorValue{D}`: A second-order tensor.
-  - `B::TensorValue{D}`: Another second-order tensor.
-
-  # Returns
-  - `TensorValue{D*D}`: A statically-sized matrix representing the fourth-order tensor.
-"""
-@generated function (⊗₁₄²³)(A::TensorValue{D}, B::TensorValue{D}) where D
-  str = ""
-  for l in 1:D
-    for k in 1:D
-      for j in 1:D
-        for i in 1:D
-          str *= "A[$i,$l]*B[$j,$k],"
-        end
-      end
-    end
-  end
-  Meta.parse("TensorValue{D*D}($str)")
-end
-
-
 function (×ᵢ⁴)(A::TensorValue{3,3,Float64})
 
   TensorValue(0.0, 0.0, 0.0, 0.0, A[9], -A[8], 0.0, -A[6], A[5], 0.0, 0.0, 0.0, -A[9],

From cfaca4a8c345331e9160cc2cdd94888aed42d293 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Miguel=20Mas=C3=B3?= <miguel.maso@upc.edu>
Date: Mon, 15 Sep 2025 12:25:12 +0200
Subject: [PATCH 06/13] docstrings

---
 src/TensorAlgebra/TensorAlgebra.jl | 41 +++++++++++++++---------------
 1 file changed, 20 insertions(+), 21 deletions(-)

diff --git a/src/TensorAlgebra/TensorAlgebra.jl b/src/TensorAlgebra/TensorAlgebra.jl
index 195c1cf..d57454a 100644
--- a/src/TensorAlgebra/TensorAlgebra.jl
+++ b/src/TensorAlgebra/TensorAlgebra.jl
@@ -48,7 +48,7 @@ export Ellipsoid
 
 
 """
-  sqrt(A::TensorValue)::TensorValue
+  sqrt(A::TensorValue{3})::TensorValue{3}
 
   Compute the square root of a 3x3 matrix by means of eigen decomposition.
 
@@ -117,6 +117,12 @@ function Gridap.TensorValues.outer(A::VectorValue{D,Float64}, B::VectorValue{D,F
   return (A ⊗₁² B)
 end
 
+
+"""
+  ⊗₁²(A::VectorValue{D}, B::VectorValue{D})::TensorValue{D,D}
+
+  Outer product of two first-order tensors (vectors), returning a second-order tensor (matrix).
+"""
 @generated function (⊗₁²)(A::VectorValue{D,Float64}, B::VectorValue{D,Float64}) where {D}
   str = ""
   for iB in 1:D
@@ -127,6 +133,13 @@ end
   Meta.parse("TensorValue{D,D, Float64}($str)")
 end
 
+
+"""
+  ⊗₁₃²⁴(A::TensorValue{D}, B::TensorValue{D})::TensorValue{D*D}
+
+  Outer product of two second-order tensors (matrices), returning a fourth-order tensor 
+  represented in a `D² x D²` flattened matrix using combined indices.
+"""
 @generated function (⊗₁₂³⁴)(A::TensorValue{D,D,Float64}, B::TensorValue{D,D,Float64}) where {D}
   str = ""
   for iB in 1:D*D
@@ -139,17 +152,10 @@ end
 
 
 """
-  ⊗₁₃²⁴(A::TensorValue, B::TensorValue)::TensorValue
+  ⊗₁₃²⁴(A::TensorValue{D}, B::TensorValue{D})::TensorValue{D*D}
 
-  Computes the **outer product** of two second-order tensors (matrices), returning a fourth-order tensor 
-  represented in a `D² x D²` flattened matrix using comined indices.
-
-  # Arguments
-  - `A::TensorValue{D}`: A second-order tensor.
-  - `B::TensorValue{D}`: Another second-order tensor.
-
-  # Returns
-  - `TensorValue{D*D}`: A statically-sized matrix representing the fourth-order tensor.
+  Outer product of two second-order tensors (matrices), returning a fourth-order tensor 
+  represented in a `D² x D²` flattened matrix using combined indices.
 """
 @generated function (⊗₁₃²⁴)(A::TensorValue{D}, B::TensorValue{D}) where D
   str = ""
@@ -167,17 +173,10 @@ end
 
 
 """
-  ⊗₁₄²³(A::TensorValue, B::TensorValue)::TensorValue
+  ⊗₁₄²³(A::TensorValue{D}, B::TensorValue{D})::TensorValue{D*D}
 
-  Computes the **outer product** of two second-order tensors (matrices), returning a fourth-order tensor 
-  represented in a `D² x D²` flattened matrix using comined indices.
-
-  # Arguments
-  - `A::TensorValue{D}`: A second-order tensor.
-  - `B::TensorValue{D}`: Another second-order tensor.
-
-  # Returns
-  - `TensorValue{D*D}`: A statically-sized matrix representing the fourth-order tensor.
+  Outer product of two second-order tensors (matrices), returning a fourth-order tensor 
+  represented in a `D² x D²` flattened matrix using combined indices.
 """
 @generated function (⊗₁₄²³)(A::TensorValue{D}, B::TensorValue{D}) where D
   str = ""

From 2fcf770fdce0427206f5b220220d1edd5fca604f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Miguel=20Mas=C3=B3?= <miguel.maso@upc.edu>
Date: Mon, 15 Sep 2025 12:40:22 +0200
Subject: [PATCH 07/13] more docstrings

---
 test/TestTensorAlgebra/TensorAlgebra.jl | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/test/TestTensorAlgebra/TensorAlgebra.jl b/test/TestTensorAlgebra/TensorAlgebra.jl
index 91d2237..3e53901 100644
--- a/test/TestTensorAlgebra/TensorAlgebra.jl
+++ b/test/TestTensorAlgebra/TensorAlgebra.jl
@@ -36,10 +36,15 @@ end
   B = TensorValue(5.0, 6.0, 7.0, 8.0)
   u = VectorValue(1.0, 2.0)
   v = VectorValue(3.0, 4.0)
+  @test u ⊗ v   == TensorValue(3.0, 6.0, 4.0, 8.0)
+  @test u ⊗₁² v == TensorValue(3.0, 6.0, 4.0, 8.0)
   @test A ⊗ B     == TensorValue(5.0, 10.0, 15.0, 20.0, 6.0, 12.0, 18.0, 24.0, 7.0, 14.0, 21.0, 28.0, 8.0, 16.0, 24.0, 32.0)
   @test A ⊗₁₂³⁴ B == TensorValue(5.0, 10.0, 15.0, 20.0, 6.0, 12.0, 18.0, 24.0, 7.0, 14.0, 21.0, 28.0, 8.0, 16.0, 24.0, 32.0)
   @test A ⊗₁₃²⁴ B == TensorValue(5.0, 10.0, 6.0, 12.0, 15.0, 20.0, 18.0, 24.0, 7.0, 14.0, 8.0, 16.0, 21.0, 28.0, 24.0, 32.0)
   @test A ⊗₁₄²³ B == TensorValue(5.0, 10.0, 6.0, 12.0, 7.0, 14.0, 8.0, 16.0, 15.0, 20.0, 18.0, 24.0, 21.0, 28.0, 24.0, 32.0)
+  @test u ⊗₁²³ A == TensorValue{2,4}(1.0, 2.0, 2.0, 4.0, 3.0, 6.0, 4.0, 8.0)
+  @test A ⊗₁₂³ u == TensorValue{2,4}(1.0, 2.0, 3.0, 4.0, 2.0, 4.0, 6.0, 8.0)
+  @test A ⊗₁₃² u == TensorValue{2,4}(1.0, 2.0, 2.0, 4.0, 3.0, 4.0, 6.0, 8.0)
 end
 
 

From b475e145d9343ba50b0ae9ad66ef9df80fcb7002 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Miguel=20Mas=C3=B3?= <miguel.maso@upc.edu>
Date: Mon, 15 Sep 2025 12:40:28 +0200
Subject: [PATCH 08/13] tests

---
 src/TensorAlgebra/TensorAlgebra.jl | 28 +++++++++++++++++++++++-----
 1 file changed, 23 insertions(+), 5 deletions(-)

diff --git a/src/TensorAlgebra/TensorAlgebra.jl b/src/TensorAlgebra/TensorAlgebra.jl
index d57454a..e3ff029 100644
--- a/src/TensorAlgebra/TensorAlgebra.jl
+++ b/src/TensorAlgebra/TensorAlgebra.jl
@@ -119,7 +119,7 @@ end
 
 
 """
-  ⊗₁²(A::VectorValue{D}, B::VectorValue{D})::TensorValue{D,D}
+  **`⊗₁²(A::VectorValue{D}, B::VectorValue{D})::TensorValue{D,D}`**
 
   Outer product of two first-order tensors (vectors), returning a second-order tensor (matrix).
 """
@@ -135,7 +135,7 @@ end
 
 
 """
-  ⊗₁₃²⁴(A::TensorValue{D}, B::TensorValue{D})::TensorValue{D*D}
+  **`⊗₁₃²⁴(A::TensorValue{D}, B::TensorValue{D})::TensorValue{D*D}`**
 
   Outer product of two second-order tensors (matrices), returning a fourth-order tensor 
   represented in a `D² x D²` flattened matrix using combined indices.
@@ -152,7 +152,7 @@ end
 
 
 """
-  ⊗₁₃²⁴(A::TensorValue{D}, B::TensorValue{D})::TensorValue{D*D}
+  **`⊗₁₃²⁴(A::TensorValue{D}, B::TensorValue{D})::TensorValue{D*D}`**
 
   Outer product of two second-order tensors (matrices), returning a fourth-order tensor 
   represented in a `D² x D²` flattened matrix using combined indices.
@@ -173,7 +173,7 @@ end
 
 
 """
-  ⊗₁₄²³(A::TensorValue{D}, B::TensorValue{D})::TensorValue{D*D}
+  **`⊗₁₄²³(A::TensorValue{D}, B::TensorValue{D})::TensorValue{D*D}`**
 
   Outer product of two second-order tensors (matrices), returning a fourth-order tensor 
   represented in a `D² x D²` flattened matrix using combined indices.
@@ -193,6 +193,12 @@ end
 end
 
 
+"""
+  **`⊗₁²³(A::TensorValue{D}, B::TensorValue{D})::TensorValue{D*D}`**
+
+  Outer product of a first-order and second-order tensors (vector by matrix),
+  returning a third-order tensor represented in a `D x D²` flattened matrix using combined indices.
+"""
 @generated function (⊗₁²³)(V::VectorValue{D,Float64}, A::TensorValue{D,D,Float64}) where {D}
   str = ""
   for iA in 1:D*D
@@ -204,6 +210,12 @@ end
 end
 
 
+"""
+  **`⊗₁²³(A::TensorValue{D}, B::TensorValue{D})::TensorValue{D*D}`**
+
+  Outer product of a second-order and first-order tensors (matrix by vector),
+  returning a third-order tensor represented in a `D x D²` flattened matrix using combined indices.
+"""
 @generated function (⊗₁₂³)(A::TensorValue{D,D,Float64}, V::VectorValue{D,Float64}) where {D}
   str = ""
   for iV in 1:D
@@ -215,7 +227,13 @@ end
 end
 
 
-function (⊗₁₃²)(A::TensorValue{3,3,Float64}, V::VectorValue{3,Float64})
+"""
+  **`⊗₁²³(A::TensorValue{D}, B::TensorValue{D})::TensorValue{D*D}`**
+
+  Outer product of a second-order and first-order tensors (matrix by vector),
+  returning a third-order tensor represented in a `D x D²` flattened matrix using combined indices.
+"""
+@generated function (⊗₁₃²)(A::TensorValue{3,3,Float64}, V::VectorValue{3,Float64})
 
   TensorValue{3,9,Float64,27}(A[1] * V[1],
     A[2] * V[1],

From 40c8f3d250b2bc5d3fb7677358e363a223622cd9 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Miguel=20Mas=C3=B3?= <miguel.maso@upc.edu>
Date: Mon, 15 Sep 2025 12:47:41 +0200
Subject: [PATCH 09/13] generated outer product

---
 src/TensorAlgebra/TensorAlgebra.jl | 64 ++++++++----------------------
 1 file changed, 16 insertions(+), 48 deletions(-)

diff --git a/src/TensorAlgebra/TensorAlgebra.jl b/src/TensorAlgebra/TensorAlgebra.jl
index e3ff029..5641ed0 100644
--- a/src/TensorAlgebra/TensorAlgebra.jl
+++ b/src/TensorAlgebra/TensorAlgebra.jl
@@ -194,9 +194,9 @@ end
 
 
 """
-  **`⊗₁²³(A::TensorValue{D}, B::TensorValue{D})::TensorValue{D*D}`**
+  **`⊗₁²³(A::VectorValue{D}, B::TensorValue{D})::TensorValue{D,D*D}`**
 
-  Outer product of a first-order and second-order tensors (vector by matrix),
+  Outer product of a first-order and second-order tensors (vector and matrix),
   returning a third-order tensor represented in a `D x D²` flattened matrix using combined indices.
 """
 @generated function (⊗₁²³)(V::VectorValue{D,Float64}, A::TensorValue{D,D,Float64}) where {D}
@@ -211,9 +211,9 @@ end
 
 
 """
-  **`⊗₁²³(A::TensorValue{D}, B::TensorValue{D})::TensorValue{D*D}`**
+  **`⊗₁²³(A::TensorValue{D}, B::VectorValue{D})::TensorValue{D,D*D}`**
 
-  Outer product of a second-order and first-order tensors (matrix by vector),
+  Outer product of a second-order and first-order tensors (matrix and vector),
   returning a third-order tensor represented in a `D x D²` flattened matrix using combined indices.
 """
 @generated function (⊗₁₂³)(A::TensorValue{D,D,Float64}, V::VectorValue{D,Float64}) where {D}
@@ -228,53 +228,21 @@ end
 
 
 """
-  **`⊗₁²³(A::TensorValue{D}, B::TensorValue{D})::TensorValue{D*D}`**
+  **`⊗₁²³(A::TensorValue{D}, B::TensorValue{D})::TensorValue{D,D*D}`**
 
-  Outer product of a second-order and first-order tensors (matrix by vector),
+  Outer product of a second-order and first-order tensors (matrix and vector),
   returning a third-order tensor represented in a `D x D²` flattened matrix using combined indices.
 """
-@generated function (⊗₁₃²)(A::TensorValue{3,3,Float64}, V::VectorValue{3,Float64})
-
-  TensorValue{3,9,Float64,27}(A[1] * V[1],
-    A[2] * V[1],
-    A[3] * V[1],
-    A[1] * V[2],
-    A[2] * V[2],
-    A[3] * V[2],
-    A[1] * V[3],
-    A[2] * V[3],
-    A[3] * V[3],
-    A[4] * V[1],
-    A[5] * V[1],
-    A[6] * V[1],
-    A[4] * V[2],
-    A[5] * V[2],
-    A[6] * V[2],
-    A[4] * V[3],
-    A[5] * V[3],
-    A[6] * V[3],
-    A[7] * V[1],
-    A[8] * V[1],
-    A[9] * V[1],
-    A[7] * V[2],
-    A[8] * V[2],
-    A[9] * V[2],
-    A[7] * V[3],
-    A[8] * V[3],
-    A[9] * V[3])
-end
-
-
-function (⊗₁₃²)(A::TensorValue{2,2,Float64}, V::VectorValue{2,Float64})
-
-  TensorValue{2,4,Float64,8}(A[1] * V[1],
-    A[2] * V[1],
-    A[1] * V[2],
-    A[2] * V[2],
-    A[3] * V[1],
-    A[4] * V[1],
-    A[3] * V[2],
-    A[4] * V[2])
+@generated function (⊗₁₃²)(A::TensorValue{D}, V::VectorValue{D}) where D
+  str = ""
+  for k in 1:D
+    for j in 1:D
+      for i in 1:D
+        str *= "A[$i,$k]*V[$j],"
+      end
+    end
+  end
+  Meta.parse("TensorValue{D,D*D}($str)")
 end
 
 

From 2105d797d4f64955289b973a1a262fc21741e1fa Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Miguel=20Mas=C3=B3?= <miguel.maso@upc.edu>
Date: Mon, 15 Sep 2025 12:48:59 +0200
Subject: [PATCH 10/13] removed old dangerous tests

---
 test/TestTensorAlgebra/TensorAlgebra.jl | 19 +------------------
 1 file changed, 1 insertion(+), 18 deletions(-)

diff --git a/test/TestTensorAlgebra/TensorAlgebra.jl b/test/TestTensorAlgebra/TensorAlgebra.jl
index 3e53901..86a08ba 100644
--- a/test/TestTensorAlgebra/TensorAlgebra.jl
+++ b/test/TestTensorAlgebra/TensorAlgebra.jl
@@ -13,25 +13,8 @@ using BenchmarkTools
   @test logreg(J; Threshold=0.01) == 0.014870878346353422
 end
 
-@testset "outer" begin
-  A = TensorValue(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0) * 1e-3
-  B = TensorValue(4.6, 2.1, 1.7, 3.2, 6.5, 1.4, 9.2, 8.0, 9.0) * 1e-3
-  V1 = VectorValue(1.0, 2.0, 3.0)
-  V2 = VectorValue(1.5, 2.5, 3.5)
-  @test norm(A ⊗ B) == 0.0002984572833756952
-  @test norm(A ⊗₁₂³⁴ B)== 0.0002984572833756952
-  @test norm(V1 ⊗ V2) == 17.04406054905931
-  @test norm(V1 ⊗₁² V2) == 17.04406054905931
-  @test norm(V1 ⊗₁²³ A)== 0.06316644678941503
-  @test norm(A ⊗₁₂³ V1)== 0.06316644678941503
-  @test norm(A ⊗₁₃² V1)== 0.06316644678941503
-  @test norm(A ⊗₁₃²⁴ B)== 0.00029845728337569516
-  @test norm(get_array(V1) ⊗ get_array(V2)) == 17.04406054905931
 
-end
-
-
-@testset "outer products" begin
+@testset "outer" begin
   A = TensorValue(1.0, 2.0, 3.0, 4.0)
   B = TensorValue(5.0, 6.0, 7.0, 8.0)
   u = VectorValue(1.0, 2.0)

From 93e04451f54201a67733cfb607ea6bf0f7c8ddc8 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Miguel=20Mas=C3=B3?= <miguel.maso@upc.edu>
Date: Mon, 15 Sep 2025 15:31:01 +0200
Subject: [PATCH 11/13] replaced outer products visco

---
 src/PhysicalModels/ViscousModels.jl | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/PhysicalModels/ViscousModels.jl b/src/PhysicalModels/ViscousModels.jl
index fad5515..914c16b 100644
--- a/src/PhysicalModels/ViscousModels.jl
+++ b/src/PhysicalModels/ViscousModels.jl
@@ -228,7 +228,7 @@ function ∂Ce_∂C(::ViscousIncompressible, γα, ∂Se∂Ce_, invUvn, Ce, Ce_t
     Ge = cof(Ce)
     ∂Se∂Ce = ∂Se∂Ce_(Ce)
     ∂Se∂Ce_trial = ∂Se∂Ce_(Ce_trial)
-    ∂Ce_trial_∂C = outer_13_24(invUvn, invUvn)
+    ∂Ce_trial_∂C = invUvn ⊗₁₃²⁴ invUvn
     #------------------------------------------
     # Derivative of return mapping with respect to Ce and λα
     #------------------------------------------   
@@ -258,7 +258,7 @@ end
 Tangent operator of Ce for at fixed Uv
 """
 function ∂Ce_∂C_Uvfixed(invUv)
-  return outer_13_24(invUv, invUv)
+  invUv ⊗₁₃²⁴ invUv
 end
 
 
@@ -267,7 +267,7 @@ end
 """
 function ∂Ce_∂invUv(C, invU)
   invU_C = invU * C
-  outer_13_24(invU_C, I3) + outer_13_24(I3, invU_C)
+  invU_C ⊗₁₃²⁴ I3 + I3 ⊗₁₃²⁴ invU_C
 end
 
 
@@ -306,7 +306,7 @@ function ViscousTangentOperator(obj::ViscousIncompressible, Δt::Float64,
   DCe_DC_Uvfixed = ∂Ce_∂C_Uvfixed(invUv)
   DCe_DinvUv = ∂Ce_∂invUv(C, invUv)
   DinvUv_DC = inv(DCe_DinvUv) * (DCe_DC - DCe_DC_Uvfixed)
-  DCDF = outer_13_24(F', I3) + outer_14_23(I3, F')
+  DCDF = F' ⊗₁₃²⁴ I3 + I3 ⊗₁₄²³ F'
   #------------------------------------------
   # 0.5*δC_{Uvfixed}:DSe[ΔC]
   #------------------------------------------
@@ -321,7 +321,7 @@ function ViscousTangentOperator(obj::ViscousIncompressible, Δt::Float64,
   # Sv:(D(δC_{Uvfixed})[ΔC])
   #------------------------------------------
   Sv = invUv_Se * invUv
-  C3 = outer_13_24(Sv, I3)
+  C3 = Sv ⊗₁₃²⁴ I3
   #------------------------------------------
   # Total Contribution
   #------------------------------------------

From 0b9de10d5cb252f5cc8e1361bb37510d3b3690cb Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Miguel=20Mas=C3=B3?= <miguel.maso@upc.edu>
Date: Mon, 15 Sep 2025 15:36:16 +0200
Subject: [PATCH 12/13] removed comments

---
 src/PhysicalModels/ViscousModels.jl | 69 +++--------------------------
 1 file changed, 5 insertions(+), 64 deletions(-)

diff --git a/src/PhysicalModels/ViscousModels.jl b/src/PhysicalModels/ViscousModels.jl
index 914c16b..93682fc 100644
--- a/src/PhysicalModels/ViscousModels.jl
+++ b/src/PhysicalModels/ViscousModels.jl
@@ -63,29 +63,6 @@ function updateStateVariables!(model::GeneralizedMaxwell, Δt, u, un, stateVars)
 end
 
 
-# """
-#   _getKinematic(::Visco, ∇u, Uv⁻¹)
-
-# Compute the kinematics of a viscous model.
-
-# # Arguments
-# - `::Visco`: A viscous model
-# - `∇u`: The deformation gradient at the considered time step
-# - `Uv⁻¹`: The inverse of the viscous strain at the considered time step
-
-# # Returns
-# - `F`
-# - `C`
-# - `Ce`
-# """
-# function _getKinematic(::Visco, ∇u, Uv⁻¹)
-#   F = one(∇u) + ∇u
-#   C = F' * F
-#   Ce = Uv⁻¹ * C * Uv⁻¹
-#   return (F, C, Ce)
-# end
-
-
 """
   ViscousStrain(Ce::TensorValue, C::TensorValue)::TensorValue
   
@@ -348,15 +325,9 @@ end
 function Piola(obj::ViscousIncompressible, Δt::Float64,
                 Se_::Function, ∂Se∂Ce_::Function,
                 F::TensorValue, Fn::TensorValue, stateVars::VectorValue)
-                # ∇u_::TensorValue, ∇un_::TensorValue, stateVars::VectorValue)
   state_vars = get_array(stateVars)
-  Uvn = TensorValue{3,3}(Tuple(state_vars[1:9]))
-  # Uvn = SMatrix{3,3}(state_vars[1:9])
+  Uvn = TensorValue{3,3}(Tuple(state_vars[1:9]))  # TODO: Update tensor slicing until next Gridap version has been released
   λαn = state_vars[10]
-#   ∇u = get_array(∇u_)
-#   ∇un = get_array(∇un_)
-  # F = get_array(F_)
-  # Fn = get_array(Fn_)
   #------------------------------------------
   # Get kinematics
   #------------------------------------------
@@ -366,24 +337,16 @@ function Piola(obj::ViscousIncompressible, Δt::Float64,
   Cn = C_(Fn)
   Ceᵗʳ = Ce_(C, invUvn)
   Cen  = Ce_(Cn, invUvn)
-#   F, C, Ceᵗʳ = _getKinematic(obj, ∇u, invUvn)
-#   _, _, Cen  = _getKinematic(obj, ∇un, invUvn)
   #------------------------------------------
   # Return mapping algorithm
   #------------------------------------------
   Ce, _ = return_mapping_algorithm!(obj, Δt, Se_, ∂Se∂Ce_, F, Ceᵗʳ, Cen, λαn)
   #------------------------------------------
-  # Get invUv and Sα
+  # Get invUv and Pα
   #------------------------------------------
   _, _, invUv = ViscousStrain(Ce, C)
   Pα = ViscousPiola(Se_, Ce, invUv, F)
-  #------------------------------------------
-  # Tangent operator
-  #------------------------------------------
-  # check_type = Pα isa TensorValue{3,3,Float64}
-  # if !check_type throw("Pα is a $(typeof(Pα))") end
   Pα
-  # return TensorValue(Pα)
 end
 
 
@@ -405,15 +368,9 @@ Visco-Elastic model for incompressible case
 function Tangent(obj::ViscousIncompressible, Δt::Float64,
                  Se_::Function, ∂Se∂Ce_::Function,
                  F::TensorValue, Fn::TensorValue, stateVars::VectorValue)
-                #  ∇u_::TensorValue, ∇un_::TensorValue, stateVars::VectorValue)
   state_vars = get_array(stateVars)
-  Uvn = TensorValue{3,3}(Tuple(state_vars[1:9]))
-  # Uvn = SMatrix{3,3}(state_vars[1:9])
+  Uvn = TensorValue{3,3}(Tuple(state_vars[1:9]))  # TODO: Update tensor slicing until next Gridap version has been released
   λαn = state_vars[10]
-#   ∇u = get_array(∇u_)
-#   ∇un = get_array(∇un_)
-#   F = get_array(F_)
-#   Fn = get_array(Fn_)
   #------------------------------------------
   # Get kinematics
   #------------------------------------------
@@ -423,8 +380,6 @@ function Tangent(obj::ViscousIncompressible, Δt::Float64,
   Cn = C_(Fn)
   Ceᵗʳ = Ce_(C, invUvn)
   Cen  = Ce_(Cn, invUvn)
-#   F, C, Ceᵗʳ = _getKinematic(obj, ∇u, invUvn)
-#   _, _, Cen  = _getKinematic(obj, ∇un, invUvn)
   #------------------------------------------
   # Return mapping algorithm
   #------------------------------------------
@@ -461,13 +416,8 @@ function ReturnMapping(obj::ViscousIncompressible, Δt::Float64,
                        Se_::Function, ∂Se∂Ce_::Function,
                        F::TensorValue, Fn::TensorValue, stateVars::VectorValue)
   state_vars = get_array(stateVars)
-  Uvn = TensorValue{3,3}(Tuple(state_vars[1:9]))
-  # Uvn = SMatrix{3,3}(state_vars[1:9])
+  Uvn = TensorValue{3,3}(Tuple(state_vars[1:9]))  # TODO: Update tensor slicing until next Gridap version has been released
   λαn = state_vars[10]
-  # ∇u = get_array(∇u_)
-  # ∇un = get_array(∇un_)
-  # F = get_array(F_)
-  # Fn = get_array(Fn_)
   #------------------------------------------
   # Get kinematics
   #------------------------------------------
@@ -477,8 +427,6 @@ function ReturnMapping(obj::ViscousIncompressible, Δt::Float64,
   Cn = C_(Fn)
   Ceᵗʳ = Ce_(C, invUvn)
   Cen  = Ce_(Cn, invUvn)
-#   F, C, Ceᵗʳ = _getKinematic(obj, ∇u, invUvn)
-#   _, _, Cen  = _getKinematic(obj, ∇un, invUvn)
   #------------------------------------------
   # Return mapping algorithm
   #------------------------------------------
@@ -496,13 +444,8 @@ function ViscousDissipation(obj::ViscousIncompressible, Δt::Float64,
                        Se_::Function, ∂Se∂Ce_::Function,
                        F::TensorValue, Fn::TensorValue, stateVars::VectorValue)
   state_vars = get_array(stateVars)
-  Uvn = TensorValue{3,3}(Tuple(state_vars[1:9]))
-  # Uvn = SMatrix{3,3}(state_vars[1:9])
+  Uvn = TensorValue{3,3}(Tuple(state_vars[1:9]))  # TODO: Update tensor slicing until next Gridap version has been released
   λαn = state_vars[10]
-  # ∇u = get_array(∇u_)
-  # ∇un = get_array(∇un_)
-  # F = get_array(F_)
-  # Fn = get_array(Fn_)
   #------------------------------------------
   # Get kinematics
   #------------------------------------------
@@ -512,8 +455,6 @@ function ViscousDissipation(obj::ViscousIncompressible, Δt::Float64,
   Cn = C_(Fn)
   Ceᵗʳ = Ce_(C, invUvn)
   Cen  = Ce_(Cn, invUvn)
-#   F, C, Ceᵗʳ = _getKinematic(obj, ∇u, invUvn)
-#   _, _, Cen  = _getKinematic(obj, ∇un, invUvn)
   #------------------------------------------
   # Return mapping algorithm
   #------------------------------------------

From 70ababec8a1f7aadc3f21079c64734d2ee1ba381 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Miguel=20Mas=C3=B3?= <miguel.maso@upc.edu>
Date: Mon, 15 Sep 2025 16:39:27 +0200
Subject: [PATCH 13/13] moved to fast implementation contraction

---
 src/TensorAlgebra/TensorAlgebra.jl      | 70 +++++++++++--------------
 test/TestTensorAlgebra/TensorAlgebra.jl | 11 +++-
 2 files changed, 40 insertions(+), 41 deletions(-)

diff --git a/src/TensorAlgebra/TensorAlgebra.jl b/src/TensorAlgebra/TensorAlgebra.jl
index 5641ed0..c495533 100644
--- a/src/TensorAlgebra/TensorAlgebra.jl
+++ b/src/TensorAlgebra/TensorAlgebra.jl
@@ -582,68 +582,58 @@ end
 
 
 """
-  contraction_IP_JPKL(Matrix1::TensorValue{D}, Matrix2::TensorValue{D*D})::TensorValue{D*D}
+  **`contraction_IP_PJKL(A::TensorValue{D}, H::TensorValue{D*D})::TensorValue{D*D}`**
 
-  Performs a specific tensor contraction between a second-order tensor `Matrix1` (of size `D × D`)
-  and a fourth-order tensor `Matrix2` (represented as a `D² x D²` matrix in Voigt or flattened index notation),
-  returning the result as a new fourth-order tensor (in the same flattened form).
-  The contraction follows the **index contraction pattern**, where addition is performed for repeated indices.
-
-  # Arguments
-  - `Matrix1::TensorValue{D, D}`: A second-order tensor (e.g., a stress or deformation tensor).
-  - `Matrix2::TensorValue{D^2, D^2}`: A fourth-order tensor written in matrix form using combined indices.
-
-  # Returns
-  - `TensorValue{D^2, D^2}`: The resulting fourth-order tensor in the same matrix representation.
+  Performs a tensor contraction between a second-order tensor (of size `D × D`)
+  and a fourth-order tensor (represented as a `D² × D²` matrix in flattened index notation).
+  The operation follows the **index contraction pattern**, where addition is performed for repeated indices.
 """
-function contraction_IP_JPKL(Matrix1::TensorValue{D}, Matrix2::TensorValue{D²}) where {D, D²}
+@generated function contraction_IP_PJKL(A::TensorValue{D}, H::TensorValue{D²}) where {D, D²}
   @assert D*D == D² "Second and Fourth-order tensors size mismatch"
-  Out  =  zeros(Float64,D*D,D*D)
-  for i in 1:D
-    for j in 1:D
-      for k in 1:D
-        for l in 1:D
+  str = ""
+  for l in 1:D
+    for k in 1:D
+      for j in 1:D
+        for i in 1:D
           for p in 1:D
-            Out[i + D * (j - 1), k + D * (l - 1)] += Matrix1[i,p] * Matrix2[j + D * (p - 1), k + D * (l - 1)]
+            a = _flat_idx(p,j,D)
+            b = _flat_idx(k,l,D)
+            str *= "+A[$i,$p]*H[$a,$b]"
           end
+          str *= ","
         end
       end
     end
   end
-  return TensorValue{D*D,D*D}(Out)
+  Meta.parse("TensorValue{D²,D², Float64}($str)")
 end
 
 
 """
-  contraction_IP_JPKL(Matrix1::TensorValue{D}, Matrix2::TensorValue{D*D})::TensorValue{D*D}
+  **`contraction_IP_JPKL(A::TensorValue{D}, H::TensorValue{D*D})::TensorValue{D*D}`**
 
-  Performs a specific tensor contraction between a second-order tensor `Matrix1` (of size `D × D`)
-  and a fourth-order tensor `Matrix2` (represented as a `D² x D²` matrix in Voigt or flattened index notation),
-  returning the result as a new fourth-order tensor (in the same flattened form).
-  The contraction follows the **index contraction pattern**, where addition is performed for repeated indices.
-
-  # Arguments
-  - `Matrix1::TensorValue{D, D}`: A second-order tensor (e.g., a stress or deformation tensor).
-  - `Matrix2::TensorValue{D^2, D^2}`: A fourth-order tensor written in matrix form using combined indices.
-
-  # Returns
-  - `TensorValue{D^2, D^2}`: The resulting fourth-order tensor in the same matrix representation.
+  Performs a tensor contraction between a second-order tensor (of size `D × D`)
+  and a fourth-order tensor (represented as a `D² × D²` matrix in flattened index notation).
+  The operation follows the **index contraction pattern**, where addition is performed for repeated indices.
 """
-function contraction_IP_PJKL(Matrix1::TensorValue{D}, Matrix2::TensorValue{D²}) where {D, D²}
+@generated function contraction_IP_JPKL(A::TensorValue{D}, H::TensorValue{D²}) where {D, D²}
   @assert D*D == D² "Second and Fourth-order tensors size mismatch"
-  Out  =  zeros(Float64,D*D,D*D)
-  for i in 1:D
-    for j in 1:D
-      for k in 1:D
-        for l in 1:D
+  str = ""
+  for l in 1:D
+    for k in 1:D
+      for j in 1:D
+        for i in 1:D
           for p in 1:D
-            Out[i + D * (j - 1), k + D * (l - 1)] += Matrix1[i,p] * Matrix2[p + D * (j - 1), k + D * (l - 1)]
+            a = _flat_idx(j,p,D)
+            b = _flat_idx(k,l,D)
+            str *= "+A[$i,$p]*H[$a,$b]"
           end
+          str *= ","
         end
       end
     end
   end
-  return TensorValue{D*D,D*D}(Out)
+  Meta.parse("TensorValue{D²,D², Float64}($str)")
 end
 
 
diff --git a/test/TestTensorAlgebra/TensorAlgebra.jl b/test/TestTensorAlgebra/TensorAlgebra.jl
index 86a08ba..d557fec 100644
--- a/test/TestTensorAlgebra/TensorAlgebra.jl
+++ b/test/TestTensorAlgebra/TensorAlgebra.jl
@@ -1,4 +1,5 @@
 using Gridap.TensorValues
+using Gridap.Arrays
 using HyperFEM.TensorAlgebra
 using Test
 using BenchmarkTools
@@ -86,7 +87,7 @@ end
 end
 
 
-@testset "Identity" begin
+@testset "identity" begin
   I2_ = TensorValue(1.0, 0.0, 0.0, 1.0)
   I3_ = TensorValue(1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0)
   I4_ = TensorValue(1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0)
@@ -109,3 +110,11 @@ end
   # @benchmark I9_
   # @benchmark I9
 end
+
+
+@testset "contraction" begin
+  A = TensorValue(1.0, 2.0, 3.0, 4.0)
+  H = TensorValue(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0)
+  @test contraction_IP_PJKL(A,H) == TensorValue(7.0, 10.0, 15.0, 22.0, 23.0, 34.0, 31.0, 46.0, 39.0, 58.0, 47.0, 70.0, 55.0, 82.0, 63.0, 94.0)
+  @test contraction_IP_JPKL(A,H) == TensorValue(10.0, 14.0, 14.0, 20.0, 26.0, 38.0, 30.0, 44.0, 42.0, 62.0, 46.0, 68.0, 58.0, 86.0, 62.0, 92.0)
+end