diff --git a/src/PhysicalModels/ViscousModels.jl b/src/PhysicalModels/ViscousModels.jl
index fad5515..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
   
@@ -228,7 +205,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 +235,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 +244,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 +283,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 +298,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
   #------------------------------------------
@@ -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
   #------------------------------------------
diff --git a/src/TensorAlgebra/TensorAlgebra.jl b/src/TensorAlgebra/TensorAlgebra.jl
index be3977a..c495533 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
@@ -31,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
  
@@ -49,20 +48,20 @@ export Ellipsoid
 
 
 """
-  msqrt(A::TensorValue)::TensorValue
+  sqrt(A::TensorValue{3})::TensorValue{3}
 
   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
 
 
@@ -118,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
@@ -128,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,6 +151,54 @@ end
 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}, 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{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}, 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
+
+
+"""
+  **`⊗₁²³(A::VectorValue{D}, B::TensorValue{D})::TensorValue{D,D*D}`**
+
+  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}
   str = ""
   for iA in 1:D*D
@@ -150,6 +210,12 @@ end
 end
 
 
+"""
+  **`⊗₁²³(A::TensorValue{D}, B::VectorValue{D})::TensorValue{D,D*D}`**
+
+  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}
   str = ""
   for iV in 1:D
@@ -161,155 +227,25 @@ end
 end
 
 
-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])
-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{D}, B::TensorValue{D})::TensorValue{D,D*D}`**
+
+  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}, 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
 
+
 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],
@@ -646,124 +582,58 @@ 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}
-
-  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.
+  **`contraction_IP_PJKL(A::TensorValue{D}, H::TensorValue{D*D})::TensorValue{D*D}`**
 
-  # 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}
-
-  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.
+  **`contraction_IP_JPKL(A::TensorValue{D}, H::TensorValue{D*D})::TensorValue{D*D}`**
 
-  # 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 08331f6..d557fec 100644
--- a/test/TestTensorAlgebra/TensorAlgebra.jl
+++ b/test/TestTensorAlgebra/TensorAlgebra.jl
@@ -1,5 +1,7 @@
 using Gridap.TensorValues
+using Gridap.Arrays
 using HyperFEM.TensorAlgebra
+using Test
 using BenchmarkTools
 
  
@@ -12,24 +14,23 @@ 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
 
+@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)
+  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
- 
- 
+
 
 # @benchmark (A ⊗₁₃²⁴ B)
 # @benchmark (A ⊗₁₂³ V1)
@@ -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