Added coercive volumetric [Constitutive models]

src/Exports.jl

@@ -42,6 +42,7 @@ end
 @publish PhysicalModels IncompressibleNeoHookean2D
 @publish PhysicalModels IncompressibleNeoHookean2D_CV
 @publish PhysicalModels VolumetricEnergy
+@publish PhysicalModels CoerciveVolumetric
 @publish PhysicalModels MooneyRivlin3D
 @publish PhysicalModels MooneyRivlin2D
 @publish PhysicalModels NonlinearMooneyRivlin3D

src/PhysicalModels/MechanicalModels.jl

@@ -1,14 +1,14 @@
 
 
 # ============================================
-# Coercive volumetric Mechanical models
+# Volumetric Mechanical models
 # ============================================
 
 """
-Coercive volumetric energy term of the form:
+Volumetric energy term of the form:
 
 ```math
-\\Psi = \\frac{1}{\\kappa} (J-1)^2
+\\Psi = \\frac{\\kappa}{2} (J-1)^2
 ```
 """
 struct VolumetricEnergy <: Volumetric
@@ -36,6 +36,33 @@ function (obj::VolumetricEnergy)(Λ::Float64=1.0)
 end
 
 
+"""
+Coercive volumetric energy term of the form:
+
+```math
+\\Psi = \\frac{\\kappa}{4} (J^2-1+\\log(J))
+```
+"""
+struct CoerciveVolumetric <: Volumetric
+  κ::Float64
+  function CoerciveVolumetric(; κ::Float64)
+    new(κ)
+  end
+end
+
+function (obj::CoerciveVolumetric)(::Float64=1.0)
+  κ = obj.κ
+  J(F) = det(F)
+  H(F) = det(F) * inv(F)'
+  Ψ(F) = (κ / 4) * (J(F)^2 - 1 - 2logreg(J(F)))
+  ∂Ψ_∂J(F) = (κ / 2) * (J(F) - ∂log∂J(J(F)))
+  ∂Ψ2_∂J2(F) = (κ / 2) * (1 - ∂∂log∂JJ(J(F)))
+  ∂Ψu(F) = ∂Ψ_∂J(F) * H(F)
+  ∂Ψuu(F) = ∂Ψ2_∂J2(F) * (H(F) ⊗ H(F)) + ×ᵢ⁴(∂Ψ_∂J(F) * F)
+  return (Ψ, ∂Ψu, ∂Ψuu)
+end
+
+
 # ============================================
 # Regularization of Mechanical models
 # ============================================

src/PhysicalModels/PhysicalModels.jl

@@ -26,6 +26,7 @@ export ARAP2D
 export ARAP2D_regularized
 export NonlinearARAP2D
 export VolumetricEnergy
+export CoerciveVolumetric
 export MooneyRivlin3D
 export MooneyRivlin2D
 export NonlinearMooneyRivlin3D

src/TensorAlgebra/Functions.jl

@@ -13,11 +13,35 @@ end
 """
 Jacobian regularization
 """
-function logreg(J; Threshold=0.01)
-  if J >= Threshold
+function logreg(J; threshold=0.01)
+  if J >= threshold
     return log(J)
   else
-    return log(Threshold) - (3.0 / 2.0) + (2 / Threshold) * J - (1 / (2 * Threshold^2)) * J^2
+    return log(threshold) - (3.0 / 2.0) + (2 / threshold) * J - (1 / (2 * threshold^2)) * J^2
+  end
+end
+
+
+"""
+Jacobian regularization
+"""
+function ∂log∂J(J; threshold=0.01)
+  if J >= threshold
+    1 / J
+  else
+    2 / threshold - J / (threshold^2)
+  end
+end
+
+
+"""
+Jacobian regularization
+"""
+function ∂∂log∂JJ(J; threshold=0.01)
+  if J >= threshold
+    -1 / (J^2)
+  else
+    -1 / (threshold^2)
   end
 end
 

src/TensorAlgebra/TensorAlgebra.jl

@@ -27,6 +27,8 @@ export zerotensor3
 export zerotensor9
 
 export logreg
+export ∂log∂J
+export ∂∂log∂JJ
 export Tensorize
 export δᵢⱼδₖₗ2D
 export δᵢₖδⱼₗ2D

test/TestConstitutiveModels/PhysicalModelTests.jl

@@ -532,12 +532,20 @@ end
 @testset "VolumetricEnergy" begin
   #  Memory estimate: 0 bytes, allocs estimate: 0.
   ∇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)
-  test_derivatives_3D_(model, Kinematics(Mechano, Solid))
+  model = VolumetricEnergy(λ=1.0)
+  test_derivatives_3D_(model, Kinematics(Mechano, Solid), rtol=1e-12)
   test_equilibrium_at_rest_3D(model)
 end
 
 
+@testset "CoerciveVolumetric" begin
+  #  Memory estimate: 0 bytes, allocs estimate: 0.
+  ∇u = TensorValue(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0) * 1e-3
+  model = CoerciveVolumetric(κ=1.0)
+  test_derivatives_3D_(model, Kinematics(Mechano, Solid), rtol=1e-12)
+  test_equilibrium_at_rest_3D(model)
+end
+
 
 
 @testset "ThermoElectroMech_Govindjee" begin

test/TestTensorAlgebra/TensorAlgebraTests.jl

@@ -9,7 +9,9 @@ using Test
   F = one(∇u) + ∇u
   J = det(F)
   @test J == 1.0149819999999996
-  @test logreg(J; Threshold=0.01) == 0.014870878346353422
+  @test logreg(J; threshold=0.01) == 0.014870878346353422
+  @test ForwardDiff.derivative(logreg, J) ≈ ∂log∂J(J)
+  @test ForwardDiff.derivative(∂log∂J, J) ≈ ∂∂log∂JJ(J)
 end