Merge pull request #206 from MultiSimOLab/fix-visco-test

Fixed return mapping test

Author: miguelmaso

src/PhysicalModels/ViscousPolyconvex.jl

@@ -94,21 +94,22 @@ end
 
 function dissipation(obj::ViscousPolyconvex, F, Fn, Cvn)
   γ = obj.μ / obj.τ
-  Τ = obj.τ / obj.Δt
+  Τ = obj.τ / obj.Δt[]
   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 ---
+# --- Return mapping and derivatives for the underlying neo-Hookean ---
 
 function return_mapping(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)
 end
 
 function ∂invCv∂C(obj::ViscousPolyconvex, C, Cn, Cvn)
@@ -118,6 +119,11 @@ function ∂invCv∂C(obj::ViscousPolyconvex, C, Cn, Cvn)
   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
 end
 
 function return_mapping(obj::ViscousPolyconvex)

test/TestConstitutiveModels/ViscousModelsTests.jl

@@ -234,9 +234,11 @@ end
 
 @testset "ViscousPolyconvex" begin
   model = ViscousPolyconvex(τ=1.15, μ=7e4)
+  update_time_step!(model, 0.1)
   Fn = I3
   F1 = I3 + 0.1*TensorValue(rand(9)...)
   C1 = F1' · F1
+  Cn = I3
   Cv = I3
 
   # --- Test neo-Hookean model ---
@@ -254,8 +256,10 @@ end
   return_mapping = HyperFEM.PhysicalModels.return_mapping
   ∂invCv∂C    = HyperFEM.PhysicalModels.∂invCv∂C
 
-  ∂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)
+  ∂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)
 
   # --- Test full model tangent operator ---
   Ψ, ∂Ψ∂F, ∂∂Ψ∂FF = model()