Merge branch 'main' into thermo-visco

Author: miguelmaso

benchmark/TensorAlgebraBenchmarks/TensorAlgebraBenchmarks.jl

@@ -47,5 +47,5 @@ A = A + A' + I3
 
 SUITE["Tensor algebra"]["δδ_μ_2d"] = @benchmarkable δᵢₖδⱼₗ2D + δᵢₗδⱼₖ2D
 SUITE["Tensor algebra"]["δδ_λ_2d"] = @benchmarkable 1.0 * δᵢⱼδₖₗ2D
-SUITE["Tensor algebra"]["Cofactor"] = cof(A)
-SUITE["Tensor algebra"]["Det(A)Inv(A')"] = det(A)*inv(A')
+SUITE["Tensor algebra"]["Cofactor"] = @benchmarkable cof(A)
+SUITE["Tensor algebra"]["Det(A)Inv(A')"] = @benchmarkable det(A)*inv(A')

src/PhysicalModels/MechanicalModels.jl

@@ -328,7 +328,7 @@ struct MooneyRivlin2D <: IsoElastic
     J(F) = det(F)
     H(F) = det(F) * inv(F)'
     Ψ(F) = (μ1 / 2 + μ2 / 2) * tr((F)' * F) + μ2 / 2.0 * J(F)^2 - (μ1 + 2 * μ2) * logreg(J(F)) +
-           (λ / 2.0) * (J(F) - 1)^2
+           (λ / 2.0) * (J(F) - 1)^2 -(μ1 + μ2) - μ2/2
     ∂Ψ_(F) = ForwardDiff.gradient(F -> Ψ(F), get_array(F))
     ∂2Ψ_(F) = ForwardDiff.jacobian(F -> ∂Ψ_(F), get_array(F))
 
@@ -356,7 +356,8 @@ struct NonlinearMooneyRivlin3D <: IsoElastic
     H(F) = det(F) * inv(F)'
 
     Ψ(F) = μ1 / (2.0 * α1 * 3.0^(α1 - 1)) * (tr((F)' * F))^α1 + μ2 / (2.0 * α2 * 3.0^(α2 - 1)) * (tr((H(F))' * H(F)))^α2 - (μ1 + 2 * μ2) * logreg(J(F)) +
-           (λ / 2.0) * (J(F) - 1)^2
+           (λ / 2.0) * (J(F) - 1)^2 +
+          -μ1/(2.0 * α1 * 3.0^(α1 - 1))*3^α1 -μ2/(2.0 * α2 * 3.0^(α2 - 1))*3^α2
 
     ∂Ψ_∂F(F) = (μ1 / (3.0^(α1 - 1)) * (tr((F)' * F))^(α1 - 1)) * F
     ∂Ψ_∂H(F) = (μ2 / (3.0^(α2 - 1)) * (tr((H(F))' * H(F)))^(α2 - 1)) * H(F)
@@ -389,7 +390,8 @@ struct NonlinearMooneyRivlin2D <: IsoElastic
     J(F) = det(F)
     H(F) = det(F) * inv(F)'
     Ψ(F) = μ1 / (2.0 * α1 * 3.0^(α1 - 1)) * (tr((F)' * F) + 1.0)^α1 + μ2 / (2.0 * α2 * 3.0^(α2 - 1)) * (tr((F)' * F) + J(F)^2)^α2 - (μ1 + 2.0 * μ2) * logreg(J(F)) +
-           (λ / 2.0) * (J(F) - 1)^2
+           (λ / 2.0) * (J(F) - 1)^2 +
+          -μ1/(2.0 * α1 * 3.0^(α1 - 1)) * 3^α1 -μ2/(2.0 * α2 * 3.0^(α2 - 1)) * 3^α2
 
     ∂Ψ_∂F(F) = ((μ1 / (3.0^(α1 - 1)) * (tr((F)' * F) + 1.0)^(α1 - 1)) + μ2 / (3.0^(α2 - 1)) * (tr((F)' * F) + J(F)^2)^(α2 - 1)) * F
     ∂log∂J(J) = J >= Threshold ? 1 / J : (2 / Threshold - J / (Threshold^2))
@@ -426,7 +428,8 @@ struct NonlinearMooneyRivlin2D_CV <: IsoElastic
     J(F) = det(F)
     H(F) = det(F) * inv(F)'
     Ψ(F) = μ1 / (2.0 * α1 * 3.0^(α1 - 1)) * (tr((F)' * F) + 1.0)^α1 + μ2 / (2.0 * α2 * 3.0^(α2 - 1)) * (tr((F)' * F) + J(F)^2)^α2 - (μ1 + 2.0 * μ2) * log(J(F)) +
-           (λ) * (J(F)^(γ) + J(F)^(-γ))
+           λ * (J(F)^(γ) + J(F)^(-γ)) +
+          -μ1/(2.0 * α1 * 3.0^(α1 - 1)) * 3^α1 -μ2/(2.0 * α2 * 3.0^(α2 - 1)) * 3^α2 -2λ
 
     ∂Ψ_∂F(F) = ((μ1 / (3.0^(α1 - 1)) * (tr((F)' * F) + 1.0)^(α1 - 1)) + μ2 / (3.0^(α2 - 1)) * (tr((F)' * F) + J(F)^2)^(α2 - 1)) * F
     ∂Ψ_∂J(F) = μ2 / (3.0^(α2 - 1)) * J(F) * (tr((F)' * F) + J(F)^2)^(α2 - 1) - (μ1 + 2.0 * μ2) * (1.0 / J(F)) + λ * γ * (J(F)^(γ - 1) - J(F)^(-γ - 1))
@@ -461,7 +464,10 @@ struct NonlinearMooneyRivlin_CV <: IsoElastic
     H(F) = det(F) * inv(F)'
     Ψ(F) = μ1 / (2.0 * α1 * 3.0^(α1 - 1)) * (tr((F)' * F))^α1 +
            μ2 / (2.0 * α2 * 3.0^(α2 - 1)) * (tr((H(F))' * H(F)))^α2 -
-           (μ1 + 2 * μ2) * log(J(F)) + λ * (J(F)^(γ) + J(F)^(-γ))
+           (μ1 + 2 * μ2) * log(J(F)) + λ * (J(F)^(γ) + J(F)^(-γ)) +
+          -μ1 / (2.0 * α1 * 3.0^(α1 - 1)) * 3^α1 +
+          -μ2 / (2.0 * α2 * 3.0^(α2 - 1)) * 3^α2 +
+          -2λ
 
     ∂Ψ_∂F(F) = ((μ1 / (3.0^(α1 - 1)) * (trAA(F))^(α1 - 1))) * F
     ∂Ψ_∂H(F) = ((μ2 / (3.0^(α2 - 1)) * (tr((H(F))' * H(F)))^(α2 - 1))) * H(F)
@@ -494,7 +500,8 @@ struct NonlinearNeoHookean_CV <: IsoElastic
     λ, μ, α, γ = obj.λ, obj.μ, obj.α, obj.γ
     J(F) = det(F)
     H(F) = det(F) * inv(F)'
-    Ψ(F) = μ / (2.0 * α * 3.0^(α - 1)) * (tr((F)' * F))^α - μ * log(J(F)) + λ * (J(F)^(γ) + J(F)^(-γ))
+    Ψ(F) = μ / (2.0 * α * 3.0^(α - 1)) * (tr((F)' * F))^α - μ * log(J(F)) + λ * (J(F)^(γ) + J(F)^(-γ)) +
+          -μ / (2.0 * α * 3.0^(α - 1)) * 3.0^α - 2λ
 
     ∂Ψ_∂F(F) = ((μ / (3.0^(α - 1)) * (tr((F)' * F))^(α - 1))) * F
     ∂Ψ_∂J(F) = -μ * (1.0 / J(F)) + λ * γ * (J(F)^(γ - 1) - J(F)^(-γ - 1))
@@ -532,8 +539,8 @@ struct NonlinearIncompressibleMooneyRivlin2D_CV <: IsoElastic
     ∂e2_∂J2(F) = (10 / 9) * J(F)^(-8 / 3) * (tr((F)' * F) + 1.0)
     ∂e2_∂FJ(F) = -(4 / 3) * J(F)^(-5 / 3) * F
 
-    Ψ1(F) = μ / (2 * α) * (e(F))^α
-    Ψ2(F) = (λ) * (J(F)^(γ) + J(F)^(-γ))
+    Ψ1(F) = μ / (2 * α) * (e(F))^α -μ / (2α) * 3^α
+    Ψ2(F) = (λ) * (J(F)^(γ) + J(F)^(-γ)) -2λ
     Ψ(F) = Ψ1(F) + Ψ2(F)
 
     ∂Ψ1_∂F(F) = (μ / 2) * (((e(F))^(α - 1.0)) * ∂e_∂F(F))
@@ -574,7 +581,8 @@ struct EightChain <: IsoElastic
       C_iso = J(F)^(-2 / 3) * C
       β = sqrt(tr(C_iso) / 3 / N)
       L = β * (3.0 - β^2) / (1.0 - β^2)
-      μ * N * (β * L + log(L / sinh(L)))
+      L0 = (3N-1)/(N-1)/sqrt(N)
+      μ * N * (β * L + log(L / sinh(L)) - L0/sqrt(N) - log(L0 / sinh(L0)))
     end
 
     ∂Ψ∂F(F) = begin
@@ -858,7 +866,7 @@ struct IncompressibleNeoHookean2D <: IsoElastic
     J1 = 0.5 * (1.0 + sqrt(1.0 + δ^2))
     ∂J1 = 0.5 * (1.0 + 1.0 / sqrt(1.0^2 + δ^2))
     β = μ * (J1^(-2 / 3) - J1^(-5 / 3) * ∂J1)
-    Ψ1(F) = μ / 2 * (tr((F)' * F) + 1.0) * J(F)^(-2 / 3)
+    Ψ1(F) = μ / 2 * (tr((F)' * F) + 1.0) * J(F)^(-2 / 3) - 3μ/2 * J(I2)^(-2 / 3)
     Ψ2(F) = (λ / 2) * (J_(F) - 1)^2
     Ψ(F) = Ψ1(F) + Ψ2(F) - β * log(J_(F))
 
@@ -892,8 +900,8 @@ struct IncompressibleNeoHookean2D_CV <: IsoElastic
     λ, μ, γ = obj.λ, obj.μ, obj.γ
     J(F) = det(F)
     H(F) = det(F) * inv(F)'
-    Ψ1(F) = μ / 2 * (tr((F)' * F) + 1.0) * J(F)^(-2 / 3)
-    Ψ2(F) = λ * (J(F)^(γ) + J(F)^(-γ))
+    Ψ1(F) = μ / 2 * (tr((F)' * F) + 1.0) * J(F)^(-2 / 3) -3μ/2
+    Ψ2(F) = λ * (J(F)^(γ) + J(F)^(-γ)) -2λ
     Ψ(F) = Ψ1(F) + Ψ2(F)
 
     ∂Ψ1_∂J(F) = -μ / 3 * (tr((F)' * F) + 1.0) * J(F)^(-5 / 3)
@@ -931,7 +939,7 @@ struct ARAP2D_regularized <: IsoElastic
     J1 = 0.5 * (1.0 + sqrt(1.0 + δ^2))
     ∂J1 = 0.5 * (1.0 + 1.0 / sqrt(1.0^2 + δ^2))
     β = μ * (J1^(-1) - J1^(-2) * ∂J1)
-    Ψ(F) = μ * 0.5 * J(F)^(-1) * (tr((F)' * F)) - β * log(J_(F))
+    Ψ(F) = μ * 0.5 * J(F)^(-1) * (tr((F)' * F)) - β * log(J_(F)) -μ*J(I2)^-1
 
     ∂Ψ1_∂J(F) = -μ / 2 * (tr((F)' * F)) * J(F)^(-2)
     ∂Ψ2_∂J(F) = -β / J_(F)
@@ -963,7 +971,7 @@ struct ARAP2D <: IsoElastic
     μ = obj.μ
     J(F) = det(F)
     H(F) = det(F) * inv(F)'
-    Ψ(F) = μ * 0.5 * J(F)^(-1) * (tr((F)' * F))
+    Ψ(F) = μ * 0.5 * J(F)^(-1) * (tr((F)' * F)) -μ
     ∂Ψ_∂F(F) = μ * F * J(F)^(-1)
     ∂Ψ_∂J(F) = -μ / 2 * (tr((F)' * F)) * J(F)^(-2)
 

test/TestConstitutiveModels/PhysicalModelTests.jl

@@ -4,6 +4,7 @@ using JSON
 using StaticArrays
 using Test
 using HyperFEM.PhysicalModels
+using HyperFEM.TensorAlgebra
 
 
 
@@ -40,6 +41,16 @@ function test_derivatives_3D_(model::PhysicalModel, K::KinematicModel; rtol=1e-1
   test_derivatives__(model, K, ∇u3, rtol=rtol, kwargs...)
 end
 
+function test_equilibrium_at_rest_2D(obj::Mechano; atol=1e-10)
+  Ψ, _... = obj()
+  @test isapprox(Ψ(I2), 0.0, atol=atol)
+end
+
+function test_equilibrium_at_rest_3D(obj::Mechano, atol=1e-10)
+  Ψ, _... = obj()
+  @test isapprox(Ψ(I3), 0.0, atol=atol)
+end
+
 
 
 
@@ -178,12 +189,14 @@ end
 @testset "NonlinearMooneyRivlin_CV" begin
   model = NonlinearMooneyRivlin_CV(λ=3.0, μ1=1.0, μ2=1.0, α1=2.0, α2=1.0, γ=6.0)
   test_derivatives_3D_(model, Kinematics(Mechano,Solid),rtol=1e-13)
+  test_equilibrium_at_rest_3D(model)
 end
 
 
 @testset "NonlinearNeoHookean_CV" begin
   model = NonlinearNeoHookean_CV(λ=3.0, μ=1.0, α=2.0, γ=6.0)
   test_derivatives_3D_(model, Kinematics(Mechano,Solid),rtol=1e-13)
+  test_equilibrium_at_rest_3D(model)
 end
 
 
@@ -208,76 +221,89 @@ end
 @testset "LinearElasticity2D" begin
   model = LinearElasticity2D(λ=3.0, μ=1.0)
   test_derivatives_2D_(model, Kinematics(Mechano,Solid))
+  test_equilibrium_at_rest_2D(model)
 end
 
 
 @testset "LinearElasticity3D" begin
   model = LinearElasticity3D(λ=3.0, μ=1.0)
   test_derivatives_3D_(model, Kinematics(Mechano,Solid))
+  test_equilibrium_at_rest_3D(model)
 end
 
 
 @testset "NeoHookean3D" begin
   model = NeoHookean3D(λ=3.0, μ=1.0)
   test_derivatives_3D_(model, Kinematics(Mechano,Solid), rtol=1e-13)
+  test_equilibrium_at_rest_3D(model)
 end
 
 
 @testset "MooneyRivlin2D" begin
   model = MooneyRivlin2D(λ=3.0, μ1=1.0, μ2=2.0)
   test_derivatives_2D_(model, Kinematics(Mechano,Solid))
+  test_equilibrium_at_rest_2D(model)
 end
 
 @testset "MooneyRivlin3D" begin
   model = MooneyRivlin3D(λ=3.0, μ1=1.0, μ2=2.0)
   test_derivatives_3D_(model,Kinematics(Mechano,Solid), rtol=1e-13)
+  test_equilibrium_at_rest_3D(model)
 end
 
 
 @testset "NonlinearMooneyRivlin2D" begin
   model = NonlinearMooneyRivlin2D(λ=(μParams[1] + μParams[2]) * 1e2, μ1=μParams[1], μ2=μParams[2], α1=μParams[3], α2=μParams[4])
   test_derivatives_2D_(model, Kinematics(Mechano,Solid))
+  test_equilibrium_at_rest_2D(model)
 end
 
 
 @testset "Yeoh3D" begin
   model = Yeoh3D(λ=3.0, C10=1.0, C20=1.0, C30=1.0)
   test_derivatives_3D_(model, Kinematics(Mechano,Solid))
+  test_equilibrium_at_rest_3D(model)
 end
 
 
 @testset "NonlinearMooneyRivlin2D_CV" begin
   model = NonlinearMooneyRivlin2D_CV(λ=(μParams[1] + μParams[2]) * 1e2, μ1=μParams[1], μ2=μParams[2], α1=μParams[3], α2=μParams[4], γ=μParams[4])
   test_derivatives_2D_(model, Kinematics(Mechano,Solid))
+  test_equilibrium_at_rest_2D(model, atol=1e-9)
 end
 
 
 @testset "NonlinearMooneyRivlin3D" begin
   model = NonlinearMooneyRivlin3D(λ=(μParams[1] + μParams[2]) * 1e2, μ1=μParams[1], μ2=μParams[2], α1=μParams[3], α2=μParams[4])
   test_derivatives_3D_(model,Kinematics(Mechano,Solid), rtol=1e-13)
+  test_equilibrium_at_rest_3D(model)
 end
 
 
 @testset "IncompressibleNeoHookean2D" begin
   model = IncompressibleNeoHookean2D(λ=(μParams[1] + μParams[2]) * 1e2, μ=μParams[1])
   test_derivatives_2D_(model, Kinematics(Mechano,Solid))
+  test_equilibrium_at_rest_2D(model)
 end
 
 @testset "IncompressibleNeoHookean2D_CV" begin
   model = IncompressibleNeoHookean2D_CV(λ=(μParams[1] + μParams[2]) * 1e2, μ=μParams[1], γ=3.0)
   test_derivatives_2D_(model, Kinematics(Mechano,Solid))
+  test_equilibrium_at_rest_2D(model)
 end
 
 
 @testset "NonlinearIncompressibleMooneyRivlin2D_CV" begin
   model = NonlinearIncompressibleMooneyRivlin2D_CV(λ=(μParams[1] + μParams[2]) * 1e2, μ=μParams[1], α=μParams[3], γ=3.0)
   test_derivatives_2D_(model, Kinematics(Mechano,Solid))
+  test_equilibrium_at_rest_2D(model)
 end
 
 
 @testset "EightChain" begin
   model = EightChain(μ=μParams[1], N=μParams[2])
   test_derivatives_3D_(model, Kinematics(Mechano,Solid),rtol=1e-13)
+  test_equilibrium_at_rest_3D(model)
 end
 
 
@@ -302,6 +328,7 @@ end
   @test Ψ(F(∇u), N) == 0.27292220826242186
   @test norm(∂Ψu(F(∇u), N)) == 100.64088114687468
   @test norm(∂Ψuu(F(∇u), N)) == 46792.35008576098
+  @test isapprox(Ψ(I2, N), 0.0, atol=1e-10)
 end
 
 
@@ -311,14 +338,16 @@ 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)
+  N /= norm(N)
   model = TransverseIsotropy3D(μ=μParams[5], α1=μParams[6], α2=μParams[7])
 
   Ψ, ∂Ψu, ∂Ψuu = model()
   K=Kinematics(Mechano,Solid)
   F, _, _ = get_Kinematics(K)
-  @test Ψ(F(∇u), N) == 269927.3350807581
-  @test norm(∂Ψu(F(∇u), N)) == 947447.8711645481
-  @test norm(∂Ψuu(F(∇u), N)) == 3.8258646319087776e6
+  @test Ψ(F(∇u), N) == 2.5259068330070704
+  @test norm(∂Ψu(F(∇u), N)) == 309.14297430663385
+  @test norm(∂Ψuu(F(∇u), N)) == 81316.15339475962
+  @test isapprox(Ψ(I3, N), 0.0, atol=1e-10)
 end
 
 
@@ -450,17 +479,8 @@ end
 @testset "VolumetricEnergy" begin
   ∇u = TensorValue(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0) * 1e-3
   model = VolumetricEnergy(λ=0.0 )
- 
-   Ψ, ∂Ψu, ∂Ψuu= model()
-  K=Kinematics(Mechano,Solid)
-  F, _, _ = get_Kinematics(K)
-  
-  ∂Ψu_(F) =TensorValue(ForwardDiff.gradient(x -> Ψ(x), get_array(F)))
-  ∂Ψuu_(F) =TensorValue(ForwardDiff.hessian(x -> Ψ(x), get_array(F)))
-
-  @test isapprox(∂Ψu(F(∇u)), ∂Ψu_(F(∇u)); rtol=1e-14)
-  @test isapprox(∂Ψuu(F(∇u)), ∂Ψuu_(F(∇u)); rtol=1e-14)
- 
+  test_derivatives_3D_(model, Kinematics(Mechano,Solid))
+  test_equilibrium_at_rest_3D(model) 
 end
 
 
@@ -741,7 +761,7 @@ end
 
 
 
-  @test Ψ(F(∇u), H0(∇φ), N) == 4.000172569336671
+  @test Ψ(F(∇u), H0(∇φ), N) == 0.0001725693366710852
   @test norm(∂Ψu(F(∇u), H0(∇φ), N)) == 0.07482084634773895
   @test norm(∂Ψφ(F(∇u), H0(∇φ), N)) == 2.793633631007779e-6
   @test norm(∂Ψuu(F(∇u), H0(∇φ), N)) == 21.74472389462642
@@ -782,7 +802,7 @@ end
   # norm(∂Ψφφ_(H0(∇φ))) -norm(∂Ψφφ(F(∇u), H0(∇φ), N))
 
 
-  @test Ψ(F(∇u), H0(∇φ), N) == 4.000172469501178
+  @test Ψ(F(∇u), H0(∇φ), N) == 0.0001724695011788059
   @test norm(∂Ψu(F(∇u), H0(∇φ), N)) == 0.07482089298212842
   @test norm(∂Ψφ(F(∇u), H0(∇φ), N)) == 2.8384487487963508e-6
   @test norm(∂Ψuu(F(∇u), H0(∇φ), N)) == 21.744723980670503
@@ -801,6 +821,7 @@ end
 @testset "ARAP2D" begin
   model = ARAP2D(μ=μParams[1])
   test_derivatives_2D_(model, Kinematics(Mechano,Solid), rtol=1e-13)
+  test_equilibrium_at_rest_2D(model)
 end
 
 
@@ -809,6 +830,7 @@ end
 @testset "ARAP2D_regularized" begin
   model = ARAP2D_regularized(μ=μParams[1])
   test_derivatives_2D_(model, Kinematics(Mechano,Solid),rtol=1e-13)
+  test_equilibrium_at_rest_2D(model)
 end
 
 
@@ -830,9 +852,10 @@ end
   #  norm(∂Ψuu_(F(∇u))) - norm(∂Ψuu(F(∇u)))
   #  norm(∂Ψu(F(∇u0)))
 
-  @test Ψ(F(∇u)) == 6440.959849358168
+  @test Ψ(F(∇u)) == 0.10816855558641691
   @test norm(∂Ψu(F(∇u))) == 52.8548808805944
   @test isapprox(norm(∂Ψuu(F(∇u))), 18128.524371074407, rtol=1e-14)
+  test_equilibrium_at_rest_2D(model)
 end
 
 
@@ -857,9 +880,10 @@ end
   #  norm(∂Ψu(F(∇u0)))
 
 
-  @test Ψ(F(∇u), J_(F(∇u))) == 6457.022976353012
+  @test Ψ(F(∇u), J_(F(∇u))) == 0.10922164405292278
   @test norm(∂Ψu(F(∇u), J_(F(∇u)))) == 52.980951554554586
   @test norm(∂Ψuu(F(∇u), J_(F(∇u)))) == 18172.854611409108
+  test_equilibrium_at_rest_2D(model)
 end
 
 

test/TestConstitutiveModels/runtests.jl

@@ -3,16 +3,16 @@ using Test
 
 @testset "ConstitutiveModels" begin
 
-    @time begin
-        include("PhysicalModelTests.jl")
-    end
+  @time begin
+    include("PhysicalModelTests.jl")
+  end
 
-    @time begin
-        include("ViscousModelsTests.jl")
-    end
+  @time begin
+    include("ViscousModelsTests.jl")
+  end
 
-    @time begin
-        include("ElectroMechanicalTests.jl")
-    end
+  @time begin
+    include("ElectroMechanicalTests.jl")
+  end
 
 end