Fixes on viscous models

src/PhysicalModels/KinematicModels.jl

@@ -18,7 +18,7 @@ struct Kinematics{A,B} <: KinematicModel
 
     function Kinematics(::Type{Visco}; F::Function=(∇u) -> one(∇u) + ∇u)
         C(F) = F' * F
-        Ce(C,Uvα⁻¹) = Uvα⁻¹ * C * Uvα⁻¹
+        Ce(C,Uvα⁻¹) = Uvα⁻¹' * C * Uvα⁻¹
         metrics = (F, C, Ce)
         A = typeof(metrics)
         new{A,Visco}(metrics)

src/PhysicalModels/PhysicalModels.jl

@@ -443,8 +443,10 @@ end
 struct ComposedElasticModel <: Elasto
   Model1::Elasto
   Model2::Elasto
+  Kinematic
   function ComposedElasticModel(model1::Elasto,model2::Elasto)
-    new(model1,model2)
+    @assert model1.Kinematic === model2.Kinematic
+    new(model1,model2,model1.Kinematic)
   end
   function (obj::ComposedElasticModel)(Λ::Float64=1.0)
     DΨ1 = obj.Model1(Λ)
@@ -547,8 +549,8 @@ struct VolumetricEnergy{A} <: Elasto
     Ψ(F) = (λ / 2.0) * (J(F) - 1)^2 
     ∂Ψ_∂J(F) =  λ * (J(F) - 1)
     ∂Ψ2_∂J2(F) = λ
-    ∂Ψu(F) =  ∂Ψ_∂J(F) * H(F)
-    ∂Ψuu(F) =∂Ψ2_∂J2(F) * (H(F) ⊗ H(F)) + ×ᵢ⁴(∂Ψ_∂J(F) * F)
+    ∂Ψu(F) = ∂Ψ_∂J(F) * H(F)
+    ∂Ψuu(F) = ∂Ψ2_∂J2(F) * (H(F) ⊗ H(F)) + ×ᵢ⁴(∂Ψ_∂J(F) * F)
     return (Ψ, ∂Ψu, ∂Ψuu)
   end
 end
@@ -1102,16 +1104,16 @@ struct IncompressibleNeoHookean3D{A} <: Elasto
   function (obj::IncompressibleNeoHookean3D)(::KinematicDescription{:SecondPiola}, Λ::Float64=1.0)
     Ψ(C) = obj.μ / 2 * tr(C) * det(C)^(-1/3)
     S(C) = begin
-      J = det(C)
+      detC = det(C)
       invC = inv(C)
-      obj.μ * J^(-1/3) * I3 - obj.μ / 3 * tr(C) * J^(-1/3) * invC
+      obj.μ * detC^(-1/3) * I3 - obj.μ / 3 * tr(C) * detC^(-1/3) * invC
     end
     ∂S∂C(C) = begin
-      J = det(C)
+      detC = det(C)
       trC = tr(C)
       invC = inv(C)
       IinvC = I3 ⊗ invC
-      1/3 * obj.μ * J^(-1/3) * (4/3*trC*invC⊗invC -(IinvC+IinvC') -trC/J*×ᵢ⁴(C))
+      1/3 * obj.μ * detC^(-1/3) * (4/3*trC*invC⊗invC -(IinvC+IinvC') -trC/detC*×ᵢ⁴(C))
     end
     return (Ψ, S, ∂S∂C)
   end

src/PhysicalModels/ViscousModels.jl

@@ -15,9 +15,10 @@ struct ViscousIncompressible{T} <: Visco
   end
   function (obj::ViscousIncompressible)(Λ::Float64=1.0; Δt::Float64)
     Ψe, Se, ∂Se∂Ce       = obj.ShortTerm(KinematicDescription{:SecondPiola}())
+    Ψ(F, Fn, state)      = Energy(obj, Δt, Ψe, Se, ∂Se∂Ce, F, Fn, state)
     ∂Ψ∂F(F, Fn, state)   = Piola(obj, Δt, Se, ∂Se∂Ce, F, Fn, state)
     ∂Ψ∂F∂F(F, Fn, state) = Tangent(obj, Δt, Se, ∂Se∂Ce, F, Fn, state)
-    return Ψe, ∂Ψ∂F, ∂Ψ∂F∂F
+    return Ψ, ∂Ψ∂F, ∂Ψ∂F∂F
   end
 end
 
@@ -117,7 +118,6 @@ function return_mapping_algorithm!(obj::ViscousIncompressible, Δt::Float64,
   for _ in 1:maxiter
     #---------- Update -----------#
     Δu = -∂res \ res[:]
-    # Ce += reshape(Δu[1:end-1], 3, 3)
     Ce += TensorValue{3,3}(Tuple(Δu[1:end-1]))  # TODO: Check reconstruction of TensorValue. ERROR: MethodError: no method matching (TensorValue{3, 3})(::Vector{Float64})
     λα += Δu[end]
     #---- Residual and jacobian ---------#
@@ -156,7 +156,7 @@ function JacobianReturnMapping(γα, Ce, Se, Se_trial, ∂Se∂Ce, F, λα)
     ∂res1_∂Ce = ∂Se∂Ce - (1-γα) * λα * ×ᵢ⁴(Ce)
     ∂res1_∂λα = -(1-γα) * Ge
     ∂res2_∂Ce = Ge
-    res = [get_array(res1)[:]; res2]  #TODO: Check the TensorValue interface
+    res = [get_array(res1)[:]; res2]
     ∂res = MMatrix{10,10}(zeros(10, 10))
     ∂res[1:9, 1:9] = get_array(∂res1_∂Ce)
     ∂res[1:9, 10] = get_array(∂res1_∂λα)[:]
@@ -180,7 +180,7 @@ Viscous 1st Piola-Kirchhoff stress
 - `Pα::SMatrix`
 """
 function ViscousPiola(Se::Function, Ce::TensorValue, invUv::TensorValue, F::TensorValue)
-    Sα = invUv * Se(Ce) * invUv
+    Sα = invUv' * Se(Ce) * invUv
     F * Sα
 end
 
@@ -232,7 +232,7 @@ end
 
 
 """
-Tangent operator of Ce for at fixed Uv
+Tangent operator of Ce at fixed Uv
 """
 function ∂Ce_∂C_Uvfixed(invUv)
   invUv ⊗₁₃²⁴ invUv
@@ -298,7 +298,7 @@ function ViscousTangentOperator(obj::ViscousIncompressible, Δt::Float64,
   # Sv:(D(δC_{Uvfixed})[ΔC])
   #------------------------------------------
   Sv = invUv_Se * invUv
-  C3 = Sv ⊗₁₃²⁴ I3
+  C3 = I3 ⊗₁₃²⁴ Sv
   #------------------------------------------
   # Total Contribution
   #------------------------------------------
@@ -307,6 +307,32 @@ function ViscousTangentOperator(obj::ViscousIncompressible, Δt::Float64,
 end
 
 
+function Energy(obj::ViscousIncompressible, Δt::Float64,
+                Ψe::Function, Se_::Function, ∂Se∂Ce_::Function,
+                F::TensorValue, Fn::TensorValue, stateVars::VectorValue)
+  state_vars = get_array(stateVars)
+  Uvn = TensorValue{3,3}(Tuple(state_vars[1:9]))  # TODO: Update tensor slicing until next Gridap version has been released
+  λαn = state_vars[10]
+  #------------------------------------------
+  # Get kinematics
+  #------------------------------------------
+  invUvn  = inv(Uvn)
+  _, C_, Ce_ = get_Kinematics(obj.Kinematic)
+  C = C_(F)
+  Cn = C_(Fn)
+  Ceᵗʳ = Ce_(C, invUvn)
+  Cen  = Ce_(Cn, invUvn)
+  #------------------------------------------
+  # Return mapping algorithm
+  #------------------------------------------
+  Ce, _ = return_mapping_algorithm!(obj, Δt, Se_, ∂Se∂Ce_, F, Ceᵗʳ, Cen, λαn)
+  #------------------------------------------
+  # Elastic energy
+  #------------------------------------------
+  Ψe(Ce)
+end
+
+
 """
   First Piola-Kirchhoff for the incompressible case
 

test/TestConstitutiveModels/ViscousModelsTests.jl

@@ -0,0 +1,167 @@
+
+using Gridap.TensorValues
+using Gridap.Arrays
+using HyperFEM.TensorAlgebra
+using HyperFEM.PhysicalModels
+using StaticArrays
+using Test
+import Base:≈
+≈(a::MultiValue, b::MultiValue; kwargs...) = ≈(get_array(a), get_array(b); kwargs...)
+
+
+μ      = 1.367e4  # Pa
+N      = 7.860e5  # -
+λ      = μ*100    # Pa
+μ1     = 3.153e5  # Pa
+τ1     = 10.72    # s
+μ2     = 5.639e5  # Pa
+τ2     = 0.82     # s
+μ3     = 1.981e5  # Pa
+τ3     = 498.8    # s
+
+
+function isochoric_F(F)
+  J = det(F)
+  @assert J > 0 "Non-physical deformation, got det(F) < 0 (det(F) = $J)"
+  J^(-1/3) * F
+end
+
+
+function numerical_piola(Ψ, F, ϵ=1e-6)
+  P = MMatrix{3,3}(zeros(Float64,9))
+  for i in 1:9
+    Fp = mutable(F)
+    Fm = mutable(F)
+    Fp[i] += ϵ
+    Fm[i] -= ϵ
+    Ψp = Ψ(TensorValue(Fp))
+    Ψm = Ψ(TensorValue(Fm))
+    P[i] = (Ψp - Ψm) / 2ϵ
+  end
+  return TensorValue(P)
+end
+
+
+function numerical_tangent(Ψ, F, ϵ=1e-5)
+  H = MMatrix{9,9}(zeros(Float64,81))
+  for j in 1:9
+    ej = zeros(9); ej[j] = ϵ
+    for i in 1:9
+      ei = zeros(9); ei[i] = ϵ
+      if i == j
+        Ψp = Ψ(F + TensorValue(ei...))
+        Ψ0 = Ψ(F)
+        Ψm = Ψ(F - TensorValue(ei...))
+        H[i,j] = (Ψp - 2Ψ0 + Ψm) / ϵ^2
+      else
+        Ψpp = Ψ(F + TensorValue(( ei+ej)...))
+        Ψpm = Ψ(F + TensorValue(( ei-ej)...))
+        Ψmp = Ψ(F + TensorValue((-ei+ej)...))
+        Ψmm = Ψ(F + TensorValue((-ei-ej)...))
+        H[i,j] = (Ψpp - Ψpm - Ψmp + Ψmm) / 4ϵ^2
+      end
+    end
+  end
+  return TensorValue(H)
+end
+
+
+function richardson_expansion(func, x, ϵ)
+  (4.0func(x,ϵ) - func(x,2ϵ)) / 3.0
+end
+
+
+function test_viscous_derivatives_numerical(model; rtolP=1e-12, rtolH=1e-12)
+  Ψ, ∂Ψu, ∂Ψuu = model(Δt = 1e-2)
+  F = TensorValue(1.:9...) * 1e-3 + I3
+  Fn = TensorValue(1.:9...) * 5e-4 + I3
+  Uvn = isochoric_F(TensorValue(1.,2.,3.,2.,4.,5.,3.,5.,6.) * 2e-4 + I3)
+  λvn = 1e-3
+  Avn = VectorValue(Uvn.data..., λvn)
+  piola = richardson_expansion((F, ϵ) -> numerical_piola(F -> Ψ(F, Fn, Avn), F, ϵ), F, 1e-5)
+  tangent = richardson_expansion((F, ϵ) -> numerical_tangent(F -> Ψ(F, Fn, Avn), F, ϵ), F, 1e-4)
+  @test isapprox(∂Ψu(F, Fn, Avn), piola, rtol=rtolP)
+  @test isapprox(∂Ψuu(F, Fn, Avn), tangent, rtol=rtolH)
+end
+
+
+function test_elastic_derivatives_numerical(model; rtolP=1e-12, rtolH=1e-12)
+  Ψ, ∂Ψu, ∂Ψuu = model()
+  F = TensorValue(1.:9...) * 1e-3 + I3
+  piola = richardson_expansion((F,ϵ) -> numerical_piola(Ψ,F,ϵ), F, 1e-5)
+  tangent = richardson_expansion((F,ϵ) -> numerical_tangent(Ψ,F,ϵ), F, 1e-4)
+  @test isapprox(∂Ψu(F), piola, rtol=rtolP)
+  @test isapprox(∂Ψuu(F), tangent, rtol=rtolH)
+end
+
+
+struct EmptyElastic <: Elasto
+  Kinematic::KinematicModel
+  function EmptyElastic()
+    new(Kinematics(Elasto))
+  end
+  function (::EmptyElastic)(Λ=0.0)
+    Ψ(F) = 0.0
+    ∂Ψu(F) = TensorValue(zeros(9)...)
+    ∂Ψuu(F) = TensorValue(zeros(81)...)
+    return (Ψ, ∂Ψu, ∂Ψuu)
+  end
+end
+
+
+@testset "VolumetricEnergy" begin
+  hyper_elastic_model = VolumetricEnergy(λ=λ)
+  test_elastic_derivatives_numerical(hyper_elastic_model, rtolP=1e-10, rtolH=1e-9)
+end;
+
+@testset "EightChain" begin
+  hyper_elastic_model = EightChain(μ=μ, N=N)
+  test_elastic_derivatives_numerical(hyper_elastic_model, rtolP=1e-3, rtolH=1e-2)
+end;
+
+@testset "EightChain+VolumetricEnergy" begin
+  hyper_elastic_model = EightChain(μ=μ, N=N) + VolumetricEnergy(λ=λ)
+  test_elastic_derivatives_numerical(hyper_elastic_model, rtolP=1e-5, rtolH=1e-4)
+end;
+
+@testset "NeoHookean3D" begin
+  hyper_elastic_model = NeoHookean3D(λ=λ, μ=μ)
+  test_elastic_derivatives_numerical(hyper_elastic_model, rtolP=1e-10, rtolH=1e-9)
+end;
+
+@testset "ViscousIncompressible" begin
+  visco = ViscousIncompressible(IncompressibleNeoHookean3D(λ=0., μ=μ1), τ1)
+  test_viscous_derivatives_numerical(visco, rtolP=1e-3, rtolH=1e-3)
+end
+
+@testset "ViscousIncompressible2" begin
+  visco = ViscousIncompressible(IncompressibleNeoHookean3D(λ=0., μ=1.0), 10.0)
+  Ψ, ∂Ψu, ∂Ψuu = visco(Δt = 0.1)
+  F    =  1e-2*TensorValue(1,2,3,4,5,6,7,8,9) + I3
+  Fn   =  5e-3*TensorValue(1,2,3,4,5,6,7,8,9) + I3
+  Fvn  =  2e-2*TensorValue(1.0,2.0,3.0,4.0,5.0,8.7,6.5,4.3,6.5) + I3
+  Cvn  =  Fvn'*Fvn
+  Uvn  =  sqrt(Cvn)
+  Avn  =  VectorValue(Uvn.data...,0.0)
+  @test isapprox(norm(∂Ψu(F, Fn, Avn)), 0.20303772905627682, rtol=1e-10)
+  @test isapprox(norm(∂Ψuu(F, Fn, Avn)), 4.847586088299776, rtol=1e-10)
+end
+
+@testset "GeneralizedMaxwell EightChain 0-branch" begin
+  hyper_elastic_model = EightChain(μ=μ, N=N) + VolumetricEnergy(λ=λ)
+  cons_model = GeneralizedMaxwell(hyper_elastic_model)
+  test_viscous_derivatives_numerical(cons_model, rtolP=1e-5, rtolH=1e-4)
+end;
+
+@testset "GeneralizedMaxwell NeoHookean 0-branch" begin
+  hyper_elastic_model = NeoHookean3D(λ=λ, μ=μ)
+  cons_model = GeneralizedMaxwell(hyper_elastic_model)
+  test_viscous_derivatives_numerical(cons_model, rtolP=1e-10, rtolH=1e-9)
+end;
+
+@testset "GeneralizedMaxwell NeoHookean 1-branch" begin
+  hyper_elastic_model = NeoHookean3D(λ=λ, μ=μ)
+  branch1 = ViscousIncompressible(IncompressibleNeoHookean3D(λ=0., μ=μ1), τ1)
+  cons_model = GeneralizedMaxwell(hyper_elastic_model, branch1)
+  test_viscous_derivatives_numerical(cons_model, rtolP=1e-3, rtolH=1e-2)
+end;

test/TestConstitutiveModels/runtests.jl

@@ -1,7 +1,13 @@
+using HyperFEM
+
 @testset "ConstitutiveModels" verbose = true begin
 
     @time begin
         include("PhysicalModels.jl")
     end
 
+    @time begin
+        include("ViscousModelsTests.jl")
+    end
+
 end