model optimisation

Author: miguelmaso

src/PhysicalModels/ViscousPolyconvex.jl

@@ -49,19 +49,19 @@ end
 
 # --- Underlying neo-Hookean model in terms of viscous distortional invariants ---
 
-function Ψv(obj::ViscousPolyconvex, C, Cv)
+function Ψv(obj::ViscousPolyconvex, C, invCv)
   μ = obj.μ
   IIIc = det(C)
-  0.5μ * (C ⊙ inv(Cv) -3*IIIc^(1/3))
+  0.5μ * (C ⊙ invCv -3*IIIc^(1/3))
 end
 
-function Sv(obj::ViscousPolyconvex, C, Cv)
+function Sv(obj::ViscousPolyconvex, C, invCv)
   μ = obj.μ
   IIIc = det(C)
-  μ * (inv(Cv) -IIIc^(1/3) * inv(C))
+  μ * (invCv -IIIc^(1/3) * inv(C))
 end
 
-function ∂Sv∂C_Cᵥfix(obj::ViscousPolyconvex, C, Cv)
+function ∂Sv∂C_Cᵥfix(obj::ViscousPolyconvex, C, invCv)
   μ = obj.μ
   IIIc = det(C)
   G    = cof(C)
@@ -72,22 +72,22 @@ end
 
 function energy(obj::ViscousPolyconvex, F, Fn, Cvn)
   C, Cn = Cauchy.((F, Fn))
-  Cv = return_mapping(obj, C, Cn, Cvn)
-  Ψv(obj, C, Cv)
+  invCv = Cv⁻¹(obj, C, Cn, Cvn)
+  Ψv(obj, C, invCv)
 end
 
 function first_piola(obj::ViscousPolyconvex, F, Fn, Cvn)
   C, Cn = Cauchy.((F, Fn))
-  Cv = return_mapping(obj, C, Cn, Cvn)
-  F * Sv(obj, C, Cv)
+  invCv = Cv⁻¹(obj, C, Cn, Cvn)
+  F * Sv(obj, C, invCv)
 end
 
 function tangent(obj::ViscousPolyconvex, F, Fn, Cvn)
   C, Cn = Cauchy.((F, Fn))
-  Cv = return_mapping(obj, C, Cn, Cvn)
-  H1 = obj.μ * ∂invCv∂C(obj, C, Cn, Cvn)
-  H2 = ∂Sv∂C_Cᵥfix(obj, C, Cv)
-  H3 = I3 ⊗₁₃²⁴ Sv(obj, C, Cv)
+  invCv = Cv⁻¹(obj, C, Cn, Cvn)
+  H1 = obj.μ * ∂Cv⁻¹∂C(obj, C, Cn, Cvn)
+  H2 = ∂Sv∂C_Cᵥfix(obj, C, invCv)
+  H3 = I3 ⊗₁₃²⁴ Sv(obj, C, invCv)
   DCDF = F' ⊗₁₃²⁴ I3 + I3 ⊗₁₄²³ F'
   0.5 * DCDF' · (H1 + H2) · DCDF + H3
 end
@@ -96,34 +96,32 @@ function dissipation(obj::ViscousPolyconvex, F, Fn, Cvn)
   γ = obj.μ / obj.τ
   Τ = obj.τ / obj.Δt[]
   C, Cn = Cauchy.((F, Fn))
-  Cv = return_mapping(obj, C, Cn, Cvn)
-  invC = inv(C)
+  invCv = Cv⁻¹(obj, C, Cn, Cvn)
+  Cv = inv(Cv)
   λ_algo = 1 / (det(invC + Τ*inv(Cvn))^(1/3) - Τ)  # λ = 3 / (Cv ⊙ invC)
-  -0.5γ * (C -λ_algo*Cv) ⊙ (invC - (1/λ_algo)*inv(Cv))
+  -0.5γ * (C -λ_algo*Cv) ⊙ (invC - (1/λ_algo)*invCv)
 end
 
 # --- Return mapping and derivatives for the underlying neo-Hookean ---
 
-function return_mapping(obj::ViscousPolyconvex, C, Cn, Cvn)
+function Cv⁻¹(obj::ViscousPolyconvex, C, Cn, Cvn)
   Τ = obj.Δt[] / obj.τ
   B = Τ * inv(C) + inv(Cvn)
-  invCv = det(B)^(-1/3) * B
-  inv(invCv)
-  # Cv = det(B)^(1/3) * inv(B)
+  det(B)^(-1/3) * B
 end
 
-function ∂invCv∂C(obj::ViscousPolyconvex, C, Cn, Cvn)
+function ∂Cv⁻¹∂C(obj::ViscousPolyconvex, C, Cn, Cvn)
   Τ = obj.Δt[] / obj.τ
-  B = Τ * inv(C) + inv(Cvn)
-  G = cof(C)
-  IIIb = det(B)
-  IIIc = det(C)
-  Τ * IIIb^(-1/3) * IIIc^(-1) * (IIsym(I3) -1/3*B⊗inv(B)) * (-IIIc^(-1) * G⊗G + ×ᵢ⁴(C))
-  # invC = inv(C)
-  # B = Τ * invC + inv(Cvn)
-  # ∂invCv∂B = det(B)^(-1/3) * (IIsym(I3) - (1/3) * (B ⊗ inv(B)))
-  # ∂B∂C = -Τ * IIsym(invC)
-  # ∂invCv∂B · ∂B∂C
+  invC = inv(C)
+  B = Τ * invC + inv(Cvn)
+  ∂invCv∂B = det(B)^(-1/3) * (IIsym(I3) - (1/3) * (B ⊗ inv(B)))
+  ∂B∂C = -Τ * IIsym(invC)
+  ∂invCv∂B · ∂B∂C
+end
+
+function return_mapping(obj::ViscousPolyconvex, C, Cn, Cvn)
+  invCv = Cv⁻¹(obj, C, Cn, Cvn)
+  inv(invCv)
 end
 
 function return_mapping(obj::ViscousPolyconvex)

test/TestConstitutiveModels/ViscousModelsTests.jl

@@ -246,20 +246,18 @@ end
   Sv          = HyperFEM.PhysicalModels.Sv
   ∂Sv∂C_Cᵥfix = HyperFEM.PhysicalModels.∂Sv∂C_Cᵥfix
 
-  Sv_ref = 2*TensorValue(ForwardDiff.gradient(C -> Ψv(model, TensorValue(C), Cv), get_array(C1)))
-  ∂Sv_ref = 2*TensorValue(ForwardDiff.hessian(C -> Ψv(model, TensorValue(C), Cv), get_array(C1)))
+  Sv_ref = 2*TensorValue(ForwardDiff.gradient(C -> Ψv(model, TensorValue(C), inv(Cv)), get_array(C1)))
+  ∂Sv_ref = 2*TensorValue(ForwardDiff.hessian(C -> Ψv(model, TensorValue(C), inv(Cv)), get_array(C1)))
 
-  @test isapprox(Sv(model, C1, Cv), Sv_ref, rtol=1e-8)
-  @test isapprox(∂Sv∂C_Cᵥfix(model, C1, Cv), ∂Sv_ref, rtol=1e-8)
+  @test isapprox(Sv(model, C1, inv(Cv)), Sv_ref, rtol=1e-8)
+  @test isapprox(∂Sv∂C_Cᵥfix(model, C1, inv(Cv)), ∂Sv_ref, rtol=1e-8)
 
   # --- Test return mapping ---
-  return_mapping = HyperFEM.PhysicalModels.return_mapping
-  ∂invCv∂C    = HyperFEM.PhysicalModels.∂invCv∂C
+  Cv⁻¹    = HyperFEM.PhysicalModels.Cv⁻¹
+  ∂Cv⁻¹∂C = HyperFEM.PhysicalModels.∂Cv⁻¹∂C
 
-  ∂invCv_ref = TensorValue(ForwardDiff.jacobian(
-    C -> get_array(inv(return_mapping(model, 0.5*TensorValue(C+C'), Cn, Cv))),  # Enforce symmetry of C within ForwardDiff
-    get_array(C1)))
-  @test isapprox(∂invCv∂C(model, C1, Cn, Cv), ∂invCv_ref, rtol=1e-8)
+  ∂invCv_ref = TensorValue(ForwardDiff.jacobian(C -> get_array(Cv⁻¹(model, 0.5*TensorValue(C+C'), Cn, Cv)), get_array(C1)))  # Enforce symmetry of C within ForwardDiff
+  @test isapprox(∂Cv⁻¹∂C(model, C1, Cn, Cv), ∂invCv_ref, rtol=1e-8)
 
   # --- Test full model tangent operator ---
   Ψ, ∂Ψ∂F, ∂∂Ψ∂FF = model()