fixed return mapping

Author: miguelmaso

src/PhysicalModels/ViscousPolyconvex.jl

@@ -68,67 +68,61 @@ function ∂Sv∂C_Cᵥfix(obj::ViscousPolyconvex, C, Cv)
   μ * IIIc^(-2/3) * (2/3 * (1/IIIc) * G ⊗ G - ×ᵢ⁴(C))
 end
 
-function ∂invCv∂C(obj::ViscousPolyconvex, C, Cn, Cvn)
-  τ = obj.τ
-  Δt = obj.Δt[]
-  invC = inv(C)
-  invCvn = inv(Cvn)
-  B = Δt/τ * invC + invCvn
-  G = cof(C)
-  IIIb = det(B)
-  IIIc = det(C)
-  0.5*Δt * IIIb^(-1/3) * IIIc^(-1) * (IIsym(I3) -1/3*B⊗inv(B)) ⊙ (-IIIc^(-1) * G⊗G + ×ᵢ⁴(C))
-end
-
 # --- Implementation of derivatives ---
 
 function energy(obj::ViscousPolyconvex, F, Fn, Cvn)
-  C = Cauchy(F)
-  Cn = Cauchy(Fn)
+  C, Cn = Cauchy.((F, Fn))
   Cv = return_mapping(obj, C, Cn, Cvn)
   Ψv(obj, C, Cv)
 end
 
 function first_piola(obj::ViscousPolyconvex, F, Fn, Cvn)
-  C = Cauchy(F)
-  Cn = Cauchy(Fn)
+  C, Cn = Cauchy.((F, Fn))
   Cv = return_mapping(obj, C, Cn, Cvn)
   F * Sv(obj, C, Cv)
 end
 
 function tangent(obj::ViscousPolyconvex, F, Fn, Cvn)
-  C = Cauchy(F)
-  Cn = Cauchy(Fn)
+  C, Cn = Cauchy.((F, Fn))
   Cv = return_mapping(obj, C, Cn, Cvn)
-  H1 = ∂Sv∂C_Cᵥfix(obj, C, Cv)
-  H2 = obj.μ * ∂invCv∂C(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)
   DCDF = F' ⊗₁₃²⁴ I3 + I3 ⊗₁₄²³ F'
-  DCDF' · (H1 + H2) · DCDF + H3
+  0.5 * DCDF' · (H1 + H2) · DCDF + H3
 end
 
 function dissipation(obj::ViscousPolyconvex, F, Fn, Cvn)
   γ = obj.μ / obj.τ
   Τ = obj.τ / obj.Δt
-  C = Cauchy(F)
-  Cn = Cauchy(Fn)
+  C, Cn = Cauchy.((F, Fn))
   Cv = return_mapping(obj, C, Cn, Cvn)
   invC = inv(C)
   λ_algo = 1 / (det(invC + Τ*inv(Cvn))^(1/3) - Τ)  # λ = 3 / (Cv ⊙ invC)
   -0.5γ * (C -λ_algo*Cv) ⊙ (invC - (1/λ_algo)*inv(Cv))
 end
 
+# --- Return mapping and derivatives for underlying neo-Hookean ---
+
 function return_mapping(obj::ViscousPolyconvex, C, Cn, Cvn)
-  Τ = obj.τ / obj.Δt[]
+  Τ = obj.Δt[] / obj.τ
   B = Τ * inv(C) + inv(Cvn)
   invCv = det(B)^(-1/3) * B
   inv(invCv)
 end
 
+function ∂invCv∂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))
+end
+
 function return_mapping(obj::ViscousPolyconvex)
   (A, F, Fn) -> begin
-    C = Cauchy(F)
-    Cn = Cauchy(Fn)
+    C, Cn = Cauchy.((F, Fn))
     Cv = return_mapping(obj, C, Cn, A)
     (true, Cv)
   end

test/TestConstitutiveModels/ViscousModelsTests.jl

@@ -232,39 +232,33 @@ end
 end
 
 
-@testset "ViscousPolyconvex_Base" begin
+@testset "ViscousPolyconvex" begin
   model = ViscousPolyconvex(τ=1.15, μ=7e4)
+  Fn = I3
   F1 = I3 + 0.1*TensorValue(rand(9)...)
   C1 = F1' · F1
   Cv = I3
 
+  # --- Test neo-Hookean model ---
   Ψv          = HyperFEM.PhysicalModels.Ψv
   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)))
 
-  @test isapprox(Sv(model, C1, Cv), Sv_ref, rtol=1e-8)  # Pass
-  @test isapprox(∂Sv∂C_Cᵥfix(model, C1, Cv), ∂Sv_ref, rtol=1e-8) # Pass
-end
+  @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 return mapping ---
+  return_mapping = HyperFEM.PhysicalModels.return_mapping
+  ∂invCv∂C    = HyperFEM.PhysicalModels.∂invCv∂C
 
-@testset "ViscousPolyconvex" begin
-  long_term = NeoHookean3D(λ=1e7, μ=8e4)
-  branch_1 = ViscousPolyconvex(τ=1.15, μ=7e4)
-  model = GeneralizedMaxwell(long_term, branch_1)
-  update_time_step!(model, 0.01)
-  Ψ, ∂Ψ∂F, ∂∂Ψ∂FF = model()
-  Fn = I3
-  A  = I3
-  F1 = I3 + 0.1*TensorValue(rand(9)...)
+  ∂invCv_ref = TensorValue(ForwardDiff.jacobian(C -> get_array(return_mapping(model, TensorValue(C), C1, Cv)), get_array(C1)))
+  @test isapprox(∂invCv∂C(model, C1, C1, Cv), ∂invCv_ref, rtol=1e-8)
 
-  # ∂Ψ∂F_ref  = TensorValue(ForwardDiff.gradient(F -> Ψ(TensorValue(F), Fn, A), get_array(F1)))
-  # ∂∂Ψ∂FF_ref = TensorValue(ForwardDiff.hessian(F -> Ψ(TensorValue(F), Fn, A), get_array(F1)))
-  ∂∂Ψ∂FF_ref = TensorValue(ForwardDiff.jacobian(F -> get_array(∂Ψ∂F(TensorValue(F), Fn, A)), get_array(F1)))
-
-  @test isapprox(∂∂Ψ∂FF(F1, Fn, A), ∂∂Ψ∂FF_ref, rtol=1e-8)
-  # @test isapprox(∂Ψ∂F(F1, Fn, A), ∂Ψ∂F_ref, rtol=1e-8)  # Fail!!!
-  # @test isapprox(∂∂Ψ∂FF(F1, Fn, A), ∂∂Ψ∂FF_ref, rtol=1e-8)
+  # --- Test full model tangent operator ---
+  Ψ, ∂Ψ∂F, ∂∂Ψ∂FF = model()
+  ∂∂Ψ∂FF_ref = TensorValue(ForwardDiff.jacobian(F -> get_array(∂Ψ∂F(TensorValue(F), Fn, Cv)), get_array(F1)))
+  @test isapprox(∂∂Ψ∂FF(F1, Fn, Cv), ∂∂Ψ∂FF_ref, rtol=1e-8)
 end