diff --git a/src/PhysicalModels/PhysicalModels.jl b/src/PhysicalModels/PhysicalModels.jl
index b5c320d..015e39c 100644
--- a/src/PhysicalModels/PhysicalModels.jl
+++ b/src/PhysicalModels/PhysicalModels.jl
@@ -56,6 +56,7 @@ export MagnetoMechModel
 export GeneralizedMaxwell
 export ViscousIncompressible
 
+export PhysicalModel
 export Mechano
 export Elasto
 export Visco
@@ -135,7 +136,6 @@ struct HessianRegularization{A,B} <: Mechano
     new{typeof(Mechano),typeof(Mechano.Kinematic)}(Mechano, δ, Mechano.Kinematic)
   end
 
-
   function (obj::HessianRegularization)(Λ::Float64=1.0)
     Ψs, ∂Ψs, ∂2Ψs = obj.Mechano()
     δ = obj.δ
@@ -146,14 +146,11 @@ struct HessianRegularization{A,B} <: Mechano
       val = real(vecval.values)
       TensorValue(vec * diagm(max.(δ, val)) * vec')
     end
-
     return (Ψs, ∂Ψs, ∂2Ψ)
   end
-
-
-
 end
 
+
 struct Hessian∇JRegularization{A,B} <: Mechano
   Mechano::A
   δ::Float64
@@ -163,7 +160,6 @@ struct Hessian∇JRegularization{A,B} <: Mechano
     new{typeof(Mechano),typeof(Mechano.Kinematic)}(Mechano, δ, κ, Mechano.Kinematic)
   end
 
-
   function (obj::Hessian∇JRegularization)(Λ::Float64=1.0)
     Ψs, ∂Ψs, ∂2Ψs = obj.Mechano()
     _, H, J = get_Kinematics(obj.Mechano.Kinematic; Λ=Λ)
@@ -179,8 +175,6 @@ struct Hessian∇JRegularization{A,B} <: Mechano
       val = real(vecval.values)
       TensorValue(vec * diagm(max.(δ, val)) * vec')
     end
-
-    # ∂2Ψ(F, Jh) = TensorValue(real(λ(F, Jh).vectors) * diagm(max.(δ, real(λ(F, Jh).values))) * real(λ(F, Jh).vectors)')
     return (Ψ, ∂Ψ, ∂2Ψ)
   end
 end
@@ -502,16 +496,15 @@ struct LinearElasticity2D{A} <: Mechano
 
   function (obj::LinearElasticity2D)(Λ::Float64=1.0)
     λ, μ = obj.λ, obj.μ
-    I22 = I2()
-    ∂Ψuu(F) = _δδ_μ_2D(μ) + _δδ_λ_2D(λ)
-    ∂Ψu(F) = ∂Ψuu(F) ⊙ (F - I22)
-    Ψ(F) = 0.5 * (F - I22) ⊙ (∂Ψuu(F) ⊙ (F - I22))
+    ε(F) = 0.5(F + F') -I2()
+    ∂Ψuu(F) = μ * (δᵢₖδⱼₗ2D + δᵢₗδⱼₖ2D) + λ * δᵢⱼδₖₗ2D
+    ∂Ψu(F) = ∂Ψuu(F) ⊙ (F - I2())
+    Ψ(F) = μ * sum(ε(F).*ε(F)) + 0.5 * λ * tr(ε(F))^2
     return (Ψ, ∂Ψu, ∂Ψuu)
   end
-
-
 end
 
+
 mutable struct LinearElasticity3D{A} <: Mechano
   λ::Float64
   μ::Float64
@@ -523,17 +516,15 @@ mutable struct LinearElasticity3D{A} <: Mechano
 
   function (obj::LinearElasticity3D)(Λ::Float64=1.0)
     λ, μ = obj.λ, obj.μ
-    I33 = I3()
-    # ∂Ψuu(F) = _δδ_μ_3D(μ) + _δδ_λ_3D(λ)
+    ε(F) = 0.5(F + F') -I3()
     ∂Ψuu(F) = μ * (δᵢₖδⱼₗ3D + δᵢₗδⱼₖ3D) + λ * δᵢⱼδₖₗ3D
-    ∂Ψu(F) = ∂Ψuu(F) ⊙ (F - I33)
-    Ψ(F) = 0.5 * (F - I33) ⊙ (∂Ψuu(F) ⊙ (F - I33))
+    ∂Ψu(F) = ∂Ψuu(F) ⊙ (F - I3())
+    Ψ(F) = μ * sum(ε(F).*ε(F)) + 0.5 * λ * tr(ε(F))^2
     return (Ψ, ∂Ψu, ∂Ψuu)
   end
-
-
 end
 
+
 struct VolumetricEnergy{A} <: Mechano
   λ::Float64
   Kinematic::A
@@ -794,11 +785,11 @@ struct NonlinearMooneyRivlin_CV{A} <: Mechano
 
 
 
-    Ψ(F::TensorValue{3, 3, Float64, 9}) = μ1 / (2.0 * α * 3.0^(α - 1)) * (tr((F)' * F))^α + 
+    Ψ(F) = μ1 / (2.0 * α * 3.0^(α - 1)) * (tr((F)' * F))^α + 
            μ2 / (2.0 * β * 3.0^(β - 1)) * (tr((H(F))' * H(F)))^β - 
            (μ1 + 2*μ2) * log(J(F)) + λ * (J(F)^(γ) + J(F)^(-γ))
 
-    ∂Ψ_∂F(F::TensorValue{3, 3, Float64, 9}) = ((μ1 / (3.0^(α - 1)) * (trAA(F))^(α - 1))) * F
+    ∂Ψ_∂F(F) = ((μ1 / (3.0^(α - 1)) * (trAA(F))^(α - 1))) * F
     ∂Ψ_∂H(F) = ((μ2 / (3.0^(β - 1)) * (tr((H(F))' * H(F)))^(β - 1))) * H(F)
     ∂Ψ_∂J(F) = -(μ1 + 2*μ2) * (1.0 / J(F)) + λ * γ * (J(F)^(γ - 1) - J(F)^(-γ - 1))
     ∂Ψu(F) = ∂Ψ_∂F(F) + ∂Ψ_∂H(F) × F + ∂Ψ_∂J(F) * H(F)
diff --git a/src/TensorAlgebra/TensorAlgebra.jl b/src/TensorAlgebra/TensorAlgebra.jl
index a244890..d38e060 100644
--- a/src/TensorAlgebra/TensorAlgebra.jl
+++ b/src/TensorAlgebra/TensorAlgebra.jl
@@ -23,6 +23,9 @@ export I2
 export I4
 export logreg
 export Tensorize
+export δᵢⱼδₖₗ2D
+export δᵢₖδⱼₗ2D
+export δᵢₗδⱼₖ2D
 export δᵢⱼδₖₗ3D
 export δᵢₖδⱼₗ3D
 export δᵢₗδⱼₖ3D
@@ -90,9 +93,13 @@ function _Kroneckerδδ(δδ::Function, N::Int)
     i, j, k, l = _full_idx4(α,N)
     δδ(i,j,k,l) ? 1.0 : 0.0
   end,
-  9*9))
+  N*N*N*N))
 end
 
+const δᵢⱼδₖₗ2D = _Kroneckerδδ((i,j,k,l) -> i==j && k==l, 2)
+const δᵢₖδⱼₗ2D = _Kroneckerδδ((i,j,k,l) -> i==k && j==l, 2)
+const δᵢₗδⱼₖ2D = _Kroneckerδδ((i,j,k,l) -> i==l && j==k, 2)
+
 const δᵢⱼδₖₗ3D = _Kroneckerδδ((i,j,k,l) -> i==j && k==l, 3)
 const δᵢₖδⱼₗ3D = _Kroneckerδδ((i,j,k,l) -> i==k && j==l, 3)
 const δᵢₗδⱼₖ3D = _Kroneckerδδ((i,j,k,l) -> i==l && j==k, 3)
@@ -102,217 +109,6 @@ function _∂H∂F_2D()
   TensorValue(0.0, 0.0, 0.0, 1.0, 0.0, 0.0, -1.0, 0.0, 0.0, -1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0)
 end
 
-function _δδ_μ_2D(μ::Float64)
-  TensorValue{4,4,Float64,16}(
-    2 * μ,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    μ,
-    μ,
-    0.0,
-    0.0,
-    μ,
-    μ,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    2.0 * μ)
-end
-
-function _δδ_λ_2D(λ::Float64)
-  TensorValue{4,4,Float64,16}(
-    λ,
-    0.0,
-    0.0,
-    λ,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    λ,
-    0.0,
-    0.0,
-    λ)
-end
-
-function _δδ_μ_3D(μ::Float64)
-  TensorValue{9,9,Float64,81}(
-    2.0 * μ,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    μ,
-    0.0,
-    μ,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    μ,
-    0.0,
-    0.0,
-    0.0,
-    μ,
-    0.0,
-    0.0,
-    0.0,
-    μ,
-    0.0,
-    μ,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    2.0 * μ,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    μ,
-    0.0,
-    μ,
-    0.0,
-    0.0,
-    0.0,
-    μ,
-    0.0,
-    0.0,
-    0.0,
-    μ,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    μ,
-    0.0,
-    μ,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    2.0 * μ)
-end
-
-function _δδ_λ_3D(λ::Float64)
-  TensorValue{9,9,Float64,81}(
-    λ,
-    0.0,
-    0.0,
-    0.0,
-    λ,
-    0.0,
-    0.0,
-    0.0,
-    λ,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    λ,
-    0.0,
-    0.0,
-    0.0,
-    λ,
-    0.0,
-    0.0,
-    0.0,
-    λ,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    0.0,
-    λ,
-    0.0,
-    0.0,
-    0.0,
-    λ,
-    0.0,
-    0.0,
-    0.0,
-    λ)
-end
-
-
 
 function Gridap.TensorValues.outer(A::TensorValue{D,D,Float64}, B::TensorValue{D,D,Float64}) where {D}
   return (A ⊗₁₂³⁴ B)
diff --git a/test/TestConstitutiveModels/PhysicalModels.jl b/test/TestConstitutiveModels/PhysicalModels.jl
index f439b82..ef89ef4 100644
--- a/test/TestConstitutiveModels/PhysicalModels.jl
+++ b/test/TestConstitutiveModels/PhysicalModels.jl
@@ -1,54 +1,52 @@
 using Gridap
+using ForwardDiff
+using StaticArrays
+using Test
+using HyperFEM.PhysicalModels
 
 
+import Base: -
 
+(-)(A::SMatrix, B::TensorValue) = A - get_array(B)  # NOTE: These functions are required for LinearElasticity to work with ForwardDiff
+(-)(A::TensorValue, B::SMatrix) = get_array(A) - B
 
 
-@testset "NonlinearMooneyRivlin_CV" begin
-  ∇u = TensorValue(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0) * 1e-3
+const ∇u2 = TensorValue(1.0, 2.0, 3.0, 4.0) * 1e-3
+const ∇u3 = TensorValue(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0) * 1e-3
+const μParams = [6456.9137547089595, 896.4633794151492,
+    1.999999451256222, 1.9999960497608036, 11747.646562400318,
+    0.7841068624959612, 1.5386288924587603]
 
-  model = NonlinearMooneyRivlin_CV(λ=3.0, μ1=1.0, μ2=1.0, α=2.0, β=1.0, γ=6.0)
 
+function test_derivatives__(model::PhysicalModel, ∇u; rtol, kwargs...)
   Ψ, ∂Ψu, ∂Ψuu = model()
-  F, _, _ = get_Kinematics(model.Kinematic)
-
-  # @benchmark norm(∂Ψuu(F(∇u)))
-  # @code_warntype  norm(∂Ψuu(F(∇u)))
+  ∂Ψu_(F) = TensorValue(ForwardDiff.gradient(Ψ, get_array(F)))
+  ∂Ψuu_(F) = TensorValue(ForwardDiff.hessian(Ψ, get_array(F)))
 
-  #  ∂Ψu_(F) =TensorValue(ForwardDiff.gradient(x -> Ψ(x), get_array(F)))
-  #  ∂Ψuu_(F) =TensorValue(ForwardDiff.hessian(x -> Ψ(x), get_array(F)))
+  F, _, _ = get_Kinematics(model.Kinematic)
+  @test isapprox(∂Ψu(F(∇u)), ∂Ψu_(F(∇u)), rtol=rtol, kwargs...)
+  @test isapprox(∂Ψuu(F(∇u)), ∂Ψuu_(F(∇u)), rtol=rtol, kwargs...)
+end
 
-  #  norm(∂Ψu_(F(∇u))) -norm(∂Ψu(F(∇u)))
-  #  norm(∂Ψuu_(F(∇u))) -norm(∂Ψuu(F(∇u)))
+function test_derivatives_2D_(model::PhysicalModel; rtol=1e-14, kwargs...)
+  test_derivatives__(model, ∇u2, rtol=rtol, kwargs...)
+end
 
+function test_derivatives_3D_(model::PhysicalModel; rtol=1e-14, kwargs...)
+  test_derivatives__(model, ∇u3, rtol=rtol, kwargs...)
+end
 
 
-  @test Ψ(F(∇u)) == 8.274742322531269
-  @test norm(∂Ψu(F(∇u))) == 5.647570016731348
-  @test norm(∂Ψuu(F(∇u))) == 653.1484437383998
 
+@testset "NonlinearMooneyRivlin_CV" begin
+  model = NonlinearMooneyRivlin_CV(λ=3.0, μ1=1.0, μ2=1.0, α=2.0, β=1.0, γ=6.0)
+  test_derivatives_3D_(model, rtol=1e-13)
 end
 
 
-
 @testset "NonlinearNeoHookean_CV" begin
-  ∇u = TensorValue(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0) * 1e-3
-
   model = NonlinearNeoHookean_CV(λ=3.0, μ=1.0, α=2.0, γ=6.0)
-
-  Ψ, ∂Ψu, ∂Ψuu = model()
-  F, _, _ = get_Kinematics(model.Kinematic)
-
-  #  ∂Ψu_(F) =TensorValue(ForwardDiff.gradient(x -> Ψ(x), get_array(F)))
-  #  ∂Ψuu_(F) =TensorValue(ForwardDiff.hessian(x -> Ψ(x), get_array(F)))
-
-  #  norm(∂Ψu_(F(∇u))) -norm(∂Ψu(F(∇u)))
-  #  norm(∂Ψuu_(F(∇u))) -norm(∂Ψuu(F(∇u)))
-
-  @test Ψ(F(∇u)) == 6.774247349061977
-  @test norm(∂Ψu(F(∇u))) == 5.578324662235092
-  @test norm(∂Ψuu(F(∇u))) == 645.1103360183206
-
+  test_derivatives_3D_(model, rtol=1e-13)
 end
 
 
@@ -71,260 +69,82 @@ end
 
 end
 
+
 @testset "LinearElasticity2D" begin
-  ∇u = TensorValue(1.0, 2.0, 3.0, 4.0) * 1e-3
   model = LinearElasticity2D(λ=3.0, μ=1.0)
-  Ψ, ∂Ψu, ∂Ψuu = model()
-  F, _, _ = get_Kinematics(model.Kinematic)
-  @test (Ψ(F(∇u))) == 6.699999999999821e-5
-  @test norm(∂Ψu(F(∇u))) == 0.029461839725311915
-  @test norm(∂Ψuu(F(∇u))) == 8.48528137423857
+  test_derivatives_2D_(model)
 end
 
+
 @testset "LinearElasticity3D" begin
-  ∇u = TensorValue(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0) * 1e-3
   model = LinearElasticity3D(λ=3.0, μ=1.0)
-  Ψ, ∂Ψu, ∂Ψuu = model()
-  F, _, _ = get_Kinematics(model.Kinematic)
-
-  @test (Ψ(F(∇u))) == 0.0006104999999999824
-  @test norm(∂Ψu(F(∇u))) == 0.09933277404764056
-  @test norm(∂Ψuu(F(∇u))) == 11.874342087037917
+  test_derivatives_3D_(model)
 end
 
 
 @testset "NeoHookean3D" begin
-  ∇u = TensorValue(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0) * 1e-3
   model = NeoHookean3D(λ=3.0, μ=1.0)
-  Ψ, ∂Ψu, ∂Ψuu = model()
-  F, _, _ = get_Kinematics(model.Kinematic)
-  @test Ψ(F(∇u)) == 0.0006083121396460722
-  @test norm(∂Ψu(F(∇u))) == 0.099612127449168118
-  @test norm(∂Ψuu(F(∇u))) == 12.073268944343628
+  test_derivatives_3D_(model, rtol=1e-13)
 end
 
 
 @testset "MooneyRivlin2D" begin
-  ∇u = TensorValue(1.0, 2.0, 3.0, 4.0) * 1e-3
   model = MooneyRivlin2D(λ=3.0, μ1=1.0, μ2=2.0)
-  Ψ, ∂Ψu, ∂Ψuu = model()
-  F, _, _ = get_Kinematics(model.Kinematic)
-  @test Ψ(F(∇u)) == 4.000175692713462
-  @test norm(∂Ψu(F(∇u))) == 0.07481942119475013
-  @test norm(∂Ψuu(F(∇u))) == 21.74472726344396
+  test_derivatives_2D_(model)
 end
 
 @testset "MooneyRivlin3D" begin
-  ∇u = TensorValue(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0) * 1e-3
   model = MooneyRivlin3D(λ=3.0, μ1=1.0, μ2=2.0)
-  Ψ, ∂Ψu, ∂Ψuu = model()
-  F, _, _ = get_Kinematics(model.Kinematic)
-  @test Ψ(F(∇u)) == 0.001598259078230413
-  @test norm(∂Ψu(F(∇u))) == 0.24833325775972206
-  @test norm(∂Ψuu(F(∇u))) == 30.36786840739546
+  test_derivatives_3D_(model, rtol=1e-13)
 end
 
 
 @testset "NonlinearMooneyRivlin2D" begin
-  ∇u = TensorValue(1.0, 2.0, 3.0, 4.0) * 1e-3
-  ∇u0 = TensorValue(0.0, 0.0, 0.0, 0.0) * 1e-3
-
-  μParams = [6456.9137547089595, 896.4633794151492,
-    1.999999451256222,
-    1.9999960497608036,
-    11747.646562400318,
-    0.7841068624959612, 1.5386288924587603]
   model = NonlinearMooneyRivlin2D(λ=(μParams[1] + μParams[2]) * 1e2, μ1=μParams[1], μ2=μParams[2], α=μParams[3], β=μParams[4])
-
-  Ψ, ∂Ψu, ∂Ψuu = model()
-  F, _, _ = get_Kinematics(model.Kinematic)
-
-  # ∂Ψu_(F) =TensorValue(ForwardDiff.gradient(x -> Ψ(x), get_array(F)))
-  # ∂Ψuu_(F) =TensorValue(ForwardDiff.hessian(x -> Ψ(x), get_array(F)))
-
-  # norm(∂Ψu_(F(∇u))) -norm(∂Ψu(F(∇u)))
-  # norm(∂Ψuu_(F(∇u))) -norm(∂Ψuu(F(∇u)))
-
-  @test Ψ(F(∇u)) == 5524.542944928555
-  @test norm(∂Ψu(F(∇u))) == 5322.887691298287
-  @test norm(∂Ψuu(F(∇u))) == 1.5136029879040532e6
-
-
+  test_derivatives_2D_(model)
 end
 
 
-
-
 @testset "Yeoh3D" begin
-  ∇u = TensorValue(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0) * 1e-3
-
   model = Yeoh3D(λ=3.0, C10=1.0, C20=1.0, C30=1.0)
-  Ψ, ∂Ψu, ∂Ψuu = model()
-
-  F, _, _ = get_Kinematics(model.Kinematic)
-
-  #    ∂Ψu_(F) =TensorValue(ForwardDiff.gradient(x -> Ψ(x), get_array(F)))
-  #  ∂Ψuu_(F) =TensorValue(ForwardDiff.hessian(x -> Ψ(x), get_array(F)))
-
-  #  norm(∂Ψu_(F(∇u))) -norm(∂Ψu(F(∇u)))
-  #  norm(∂Ψuu_(F(∇u))) -norm(∂Ψuu(F(∇u)))
-
-
-
-  @test Ψ(F(∇u)) == 0.0018248918516909718
-  @test norm(∂Ψu(F(∇u))) == 0.33825920882848515
-  @test norm(∂Ψuu(F(∇u))) == 40.845986774511886
+  test_derivatives_3D_(model)
 end
 
 
-
-
 @testset "NonlinearMooneyRivlin2D_CV" begin
-  ∇u = TensorValue(1.0, 2.0, 3.0, 4.0) * 1e-3
-  ∇u0 = TensorValue(0.0, 0.0, 0.0, 0.0) * 1e-3
-
-  μParams = [6456.9137547089595, 896.4633794151492,
-    1.999999451256222,
-    1.9999960497608036,
-    11747.646562400318,
-    0.7841068624959612, 1.5386288924587603]
   model = NonlinearMooneyRivlin2D_CV(λ=(μParams[1] + μParams[2]) * 1e2, μ1=μParams[1], μ2=μParams[2], α=μParams[3], β=μParams[4], γ=μParams[4])
-
-  Ψ, ∂Ψu, ∂Ψuu = model()
-  F, _, _ = get_Kinematics(model.Kinematic)
-
-  # ∂Ψu_(F) =TensorValue(ForwardDiff.gradient(x -> Ψ(x), get_array(F)))
-  # ∂Ψuu_(F) =TensorValue(ForwardDiff.hessian(x -> Ψ(x), get_array(F)))
-
-  # norm(∂Ψu_(F(∇u))) -norm(∂Ψu(F(∇u)))
-  # norm(∂Ψuu_(F(∇u))) -norm(∂Ψuu(F(∇u)))
-
-  @test Ψ(F(∇u)) == 1.47626389512027e6
-  @test norm(∂Ψu(F(∇u))) == 41486.412892304914
-  @test norm(∂Ψuu(F(∇u))) == 1.171028839080193e7
-
-
+  test_derivatives_2D_(model)
 end
 
 
-
 @testset "NonlinearMooneyRivlin3D" begin
-  ∇u = TensorValue(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0) * 1e-3
-
-  μParams = [6456.9137547089595, 896.4633794151492,
-    1.999999451256222,
-    1.9999960497608036,
-    11747.646562400318,
-    0.7841068624959612, 1.5386288924587603]
   model = NonlinearMooneyRivlin3D(λ=(μParams[1] + μParams[2]) * 1e2, μ1=μParams[1], μ2=μParams[2], α=μParams[3], β=μParams[4])
-
-  Ψ, ∂Ψu, ∂Ψuu = model()
-  F, _, _ = get_Kinematics(model.Kinematic)
-
-  # ∂Ψu_(F) =TensorValue(ForwardDiff.gradient(x -> Ψ(x), get_array(F)))
-  # ∂Ψuu_(F) =TensorValue(ForwardDiff.hessian(x -> Ψ(x), get_array(F)))
-
-  # norm(∂Ψu_(F(∇u))) -norm(∂Ψu(F(∇u)))
-  # norm(∂Ψuu_(F(∇u))) -norm(∂Ψuu(F(∇u)))
-
-  @test Ψ(F(∇u)) == 5600.526853280392
-  @test norm(∂Ψu(F(∇u))) == 19622.49002309361
-  @test norm(∂Ψuu(F(∇u))) == 2.313386876448522e6
-
+  test_derivatives_3D_(model, rtol=1e-13)
 end
 
 
 @testset "IncompressibleNeoHookean2D" begin
-  ∇u = TensorValue(1.0, 2.0, 3.0, 4.0) * 1e-3
-  ∇u0 = TensorValue(1.0, 2.0, 3.0, 4.0) * 0.0
-
-  μParams = [6456.9137547089595, 896.4633794151492,
-    1.999999451256222,
-    1.9999960497608036,
-    11747.646562400318,
-    0.7841068624959612, 1.5386288924587603]
   model = IncompressibleNeoHookean2D(λ=(μParams[1] + μParams[2]) * 1e2, μ=μParams[1])
-
-  Ψ, ∂Ψu, ∂Ψuu = model()
-  F, _, J_ = get_Kinematics(model.Kinematic)
-
-  # ∂Ψu_(F) =TensorValue(ForwardDiff.gradient(x -> Ψ(x), get_array(F)))
-  # ∂Ψuu_(F) =TensorValue(ForwardDiff.hessian(x -> Ψ(x), get_array(F)))
-
-  # norm(∂Ψu_(F(∇u))) - norm(∂Ψu(F(∇u)))
-  #  norm(∂Ψuu_(F(∇u))) - norm(∂Ψuu(F(∇u)))
-
-  @test Ψ(F(∇u)) == 9678.62138906247
-  @test norm(∂Ψu(F(∇u))) == 5225.228283839419
-  @test norm(∂Ψuu(F(∇u))) == 1.4859216839324834e6
+  test_derivatives_2D_(model)
 end
 
 @testset "IncompressibleNeoHookean2D_CV" begin
-  ∇u = TensorValue(1.0, 2.0, 3.0, 4.0) * 1e-3
-  ∇u0 = TensorValue(1.0, 2.0, 3.0, 4.0) * 0.0
-
-  μParams = [6456.9137547089595, 896.4633794151492,
-    1.999999451256222,
-    1.9999960497608036,
-    11747.646562400318,
-    0.7841068624959612, 1.5386288924587603]
   model = IncompressibleNeoHookean2D_CV(λ=(μParams[1] + μParams[2]) * 1e2, μ=μParams[1], γ=3.0)
-
-  Ψ, ∂Ψu, ∂Ψuu = model()
-  F, _, J_ = get_Kinematics(model.Kinematic)
-
-  ∂Ψu_(F) = TensorValue(ForwardDiff.gradient(x -> Ψ(x), get_array(F)))
-  ∂Ψuu_(F) = TensorValue(ForwardDiff.hessian(x -> Ψ(x), get_array(F)))
-
-  norm(∂Ψu_(F(∇u))) - norm(∂Ψu(F(∇u)))
-  norm(∂Ψuu_(F(∇u))) - norm(∂Ψuu(F(∇u)))
-
-  @test Ψ(F(∇u)) == 1.4805254328174442e6
-  @test norm(∂Ψu(F(∇u))) == 93109.55194565043
-  @test norm(∂Ψuu(F(∇u))) == 2.6282687148741454e7
+  test_derivatives_2D_(model)
 end
 
 
 @testset "NonlinearIncompressibleMooneyRivlin2D_CV" begin
-  ∇u = TensorValue(1.0, 2.0, 3.0, 4.0) * 1e-3
-  ∇u0 = TensorValue(1.0, 2.0, 3.0, 4.0) * 0.0
-
-  μParams = [6456.9137547089595, 896.4633794151492,
-    1.999999451256222,
-    1.9999960497608036,
-    11747.646562400318,
-    0.7841068624959612, 1.5386288924587603]
   model = NonlinearIncompressibleMooneyRivlin2D_CV(λ=(μParams[1] + μParams[2]) * 1e2, μ=μParams[1], α=μParams[3], γ=3.0)
-
-  Ψ, ∂Ψu, ∂Ψuu = model()
-  F, _, J_ = get_Kinematics(model.Kinematic)
-
-  # ∂Ψu_(F) =TensorValue(ForwardDiff.gradient(x -> Ψ(x), get_array(F)))
-  # ∂Ψuu_(F) =TensorValue(ForwardDiff.hessian(x -> Ψ(x), get_array(F)))
-
-  # norm(∂Ψu_(F(∇u))) - norm(∂Ψu(F(∇u)))
-  #  norm(∂Ψuu_(F(∇u))) - norm(∂Ψuu(F(∇u)))
-
-
-  @test Ψ(F(∇u)) == 1.4853683856379557e6
-  @test norm(∂Ψu(F(∇u))) == 93139.40884969762
-  @test norm(∂Ψuu(F(∇u))) == 2.629114137261709e7
+  test_derivatives_2D_(model)
 end
 
 
-
-
-
 @testset "TransverseIsotropy2D" begin
   ∇u = TensorValue(1.0, 2.0, 3.0, 4.0) * 1e-3
   ∇u0 = TensorValue(1.0, 2.0, 3.0, 4.0) * 0.0
 
   N = VectorValue(1.0, 2.0) / sqrt(5.0)
-  μParams = [6456.9137547089595, 896.4633794151492,
-    1.999999451256222,
-    1.9999960497608036,
-    11747.646562400318,
-    0.7841068624959612, 1.5386288924587603]
   model = TransverseIsotropy2D(μ=μParams[5], α=μParams[6], β=μParams[7])
   Ψ, ∂Ψu, ∂Ψuu = model()
 
@@ -349,11 +169,6 @@ end
 @testset "TransverseIsotropy3D" begin
   ∇u = TensorValue(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0) * 1e-3
   N = VectorValue(1.0, 2.0, 3.0)
-  μParams = [6456.9137547089595, 896.4633794151492,
-    1.999999451256222,
-    1.9999960497608036,
-    11747.646562400318,
-    0.7841068624959612, 1.5386288924587603]
   model = TransverseIsotropy3D(μ=μParams[5], α=μParams[6], β=μParams[7])
 
   Ψ, ∂Ψu, ∂Ψuu = model()
@@ -843,58 +658,16 @@ end
 
 
 @testset "ARAP2D" begin
-  ∇u = TensorValue(1.0, 2.0, 3.0, 4.0) * 1e-3
-  ∇u0 = TensorValue(1.0, 2.0, 3.0, 4.0) * 0.0
-
-  μParams = [6456.9137547089595, 896.4633794151492,
-    1.999999451256222,
-    1.9999960497608036,
-    11747.646562400318,
-    0.7841068624959612, 1.5386288924587603]
   model = ARAP2D(μ=μParams[1])
-
-  Ψ, ∂Ψu, ∂Ψuu = model()
-  F, _, J_ = get_Kinematics(model.Kinematic)
-
-  # ∂Ψu_(F) =TensorValue(ForwardDiff.gradient(x -> Ψ(x), get_array(F)))
-  # ∂Ψuu_(F) =TensorValue(ForwardDiff.hessian(x -> Ψ(x), get_array(F)))
-
-  # norm(∂Ψu_(F(∇u))) - norm(∂Ψu(F(∇u)))
-  #  norm(∂Ψuu_(F(∇u))) - norm(∂Ψuu(F(∇u)))
-
-  @test Ψ(F(∇u)) == 6457.022976353012
-  @test norm(∂Ψu(F(∇u))) == 52.980951554554586
-  @test norm(∂Ψuu(F(∇u))) == 18172.854404408303
+  test_derivatives_2D_(model, rtol=1e-13)
 end
 
 
 
 
 @testset "ARAP2D_regularized" begin
-  ∇u = TensorValue(1.0, 2.0, 3.0, 4.0) * 1e-3
-  ∇u0 = TensorValue(1.0, 2.0, 3.0, 4.0) * 0.0
-
-  μParams = [6456.9137547089595, 896.4633794151492,
-    1.999999451256222,
-    1.9999960497608036,
-    11747.646562400318,
-    0.7841068624959612, 1.5386288924587603]
   model = ARAP2D_regularized(μ=μParams[1])
-
-  Ψ, ∂Ψu, ∂Ψuu = model()
-  F, _, J_ = get_Kinematics(model.Kinematic)
-
-  # ∂Ψu_(F) =TensorValue(ForwardDiff.gradient(x -> Ψ(x), get_array(F)))
-  # ∂Ψuu_(F) =TensorValue(ForwardDiff.hessian(x -> Ψ(x), get_array(F)))
-
-  # norm(∂Ψu_(F(∇u))) - norm(∂Ψu(F(∇u)))
-  #  norm(∂Ψuu_(F(∇u))) - norm(∂Ψuu(F(∇u)))
-  #  norm(∂Ψu(F(∇u0)))
-
-
-  @test Ψ(F(∇u)) == 6440.959849358168
-  @test norm(∂Ψu(F(∇u))) == 52.8548808805944
-  @test norm(∂Ψuu(F(∇u))) == 18128.952058660318
+  test_derivatives_2D_(model, rtol=1e-13)
 end
 
 
@@ -902,11 +675,6 @@ end
   ∇u = TensorValue(1.0, 2.0, 3.0, 4.0) * 1e-3
   ∇u0 = TensorValue(1.0, 2.0, 3.0, 4.0) * 0.0
 
-  μParams = [6456.9137547089595, 896.4633794151492,
-    1.999999451256222,
-    1.9999960497608036,
-    11747.646562400318,
-    0.7841068624959612, 1.5386288924587603]
   model = ARAP2D_regularized(μ=μParams[1])
   modelreg = HessianRegularization(Mechano=model)
 
diff --git a/test/TestTensorAlgebra/TensorAlgebraBenchmarks.jl b/test/TestTensorAlgebra/TensorAlgebraBenchmarks.jl
new file mode 100644
index 0000000..749098c
--- /dev/null
+++ b/test/TestTensorAlgebra/TensorAlgebraBenchmarks.jl
@@ -0,0 +1,49 @@
+using HyperFEM.TensorAlgebra
+using BenchmarkTools
+
+
+function _δδ_μ_2D(μ::Float64)
+  TensorValue{4,4,Float64,16}(
+    2 * μ,
+    0.0,
+    0.0,
+    0.0,
+    0.0,
+    μ,
+    μ,
+    0.0,
+    0.0,
+    μ,
+    μ,
+    0.0,
+    0.0,
+    0.0,
+    0.0,
+    2.0 * μ)
+end
+
+function _δδ_λ_2D(λ::Float64)
+  TensorValue{4,4,Float64,16}(
+    λ,
+    0.0,
+    0.0,
+    λ,
+    0.0,
+    0.0,
+    0.0,
+    0.0,
+    0.0,
+    0.0,
+    0.0,
+    0.0,
+    λ,
+    0.0,
+    0.0,
+    λ)
+end
+
+@benchmark _δδ_μ_2D(1.0)
+@benchmark 1.0 * ((δᵢₖδⱼₗ2D + δᵢₗδⱼₖ2D))
+
+@benchmark _δδ_λ_2D(1.0)
+@benchmark 1.0 * δᵢⱼδₖₗ2D
diff --git a/test/runtests.jl b/test/runtests.jl
index 4d26fae..8e0e807 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -5,6 +5,7 @@ using Gridap, GridapGmsh, GridapSolvers, DrWatson, TimerOutputs
 using GridapSolvers.NonlinearSolvers
 using Gridap.FESpaces
 using Gridap.CellData
+using Gridap.TensorValues
 using HyperFEM: jacobian, IterativeSolver, solve!
 using WriteVTK
 using Revise
@@ -15,15 +16,9 @@ using BenchmarkTools
 
 import Base.isapprox
 
-function isapprox(a::TensorValue, b::TensorValue; kwargs...)
-    return isapprox(get_array(a), get_array(b); kwargs...)
-end
+isapprox(A::MultiValue, B::MultiValue; kwargs...) = isapprox(get_array(A), get_array(B); kwargs...)
 
 
-function isapprox(a::VectorValue, b::VectorValue; kwargs...)
-    return isapprox(get_array(a), get_array(b); kwargs...)
-end
-
 @testset "HyperFEMTests" verbose = true begin
 
     @time begin