From 6aff145c7095551810266fd1b7225cd0b0afccd9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3?= Date: Fri, 31 Jul 2026 12:11:46 +0200 Subject: [PATCH 1/6] Added coercive volumetric --- src/Exports.jl | 1 + src/PhysicalModels/MechanicalModels.jl | 33 +++++++++++++++++-- src/PhysicalModels/PhysicalModels.jl | 1 + .../PhysicalModelTests.jl | 10 +++++- 4 files changed, 41 insertions(+), 4 deletions(-) diff --git a/src/Exports.jl b/src/Exports.jl index 4aef2f4..6670e98 100644 --- a/src/Exports.jl +++ b/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 diff --git a/src/PhysicalModels/MechanicalModels.jl b/src/PhysicalModels/MechanicalModels.jl index f55cc7e..0681ddf 100644 --- a/src/PhysicalModels/MechanicalModels.jl +++ b/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 + 2log(J(F))) + ∂Ψ_∂J(F) = (κ / 2) * (J(F) + 1/J(F)) + ∂Ψ2_∂J2(F) = (κ / 2) * (1 - 1/J(F)^2) + ∂Ψ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 # ============================================ diff --git a/src/PhysicalModels/PhysicalModels.jl b/src/PhysicalModels/PhysicalModels.jl index 2532a65..25360f7 100644 --- a/src/PhysicalModels/PhysicalModels.jl +++ b/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 diff --git a/test/TestConstitutiveModels/PhysicalModelTests.jl b/test/TestConstitutiveModels/PhysicalModelTests.jl index 1337762..fb45314 100644 --- a/test/TestConstitutiveModels/PhysicalModelTests.jl +++ b/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) + model = VolumetricEnergy(λ=1.0) test_derivatives_3D_(model, Kinematics(Mechano, Solid)) 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)) + test_equilibrium_at_rest_3D(model) +end + @testset "ThermoElectroMech_Govindjee" begin From f9bec3f47c0b658d6af5cfb14d3935a7575344af Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3?= Date: Fri, 31 Jul 2026 12:24:03 +0200 Subject: [PATCH 2/6] fixed sign in volumetric energy --- src/PhysicalModels/MechanicalModels.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/PhysicalModels/MechanicalModels.jl b/src/PhysicalModels/MechanicalModels.jl index 0681ddf..4499837 100644 --- a/src/PhysicalModels/MechanicalModels.jl +++ b/src/PhysicalModels/MechanicalModels.jl @@ -54,9 +54,9 @@ function (obj::CoerciveVolumetric)(::Float64=1.0) κ = obj.κ J(F) = det(F) H(F) = det(F) * inv(F)' - Ψ(F) = (κ / 4) * (J(F)^2 - 1 + 2log(J(F))) - ∂Ψ_∂J(F) = (κ / 2) * (J(F) + 1/J(F)) - ∂Ψ2_∂J2(F) = (κ / 2) * (1 - 1/J(F)^2) + Ψ(F) = (κ / 4) * (J(F)^2 - 1 - 2log(J(F))) + ∂Ψ_∂J(F) = (κ / 2) * (J(F) - 1/J(F)) + ∂Ψ2_∂J2(F) = (κ / 2) * (1 + 1/J(F)^2) ∂Ψu(F) = ∂Ψ_∂J(F) * H(F) ∂Ψuu(F) = ∂Ψ2_∂J2(F) * (H(F) ⊗ H(F)) + ×ᵢ⁴(∂Ψ_∂J(F) * F) return (Ψ, ∂Ψu, ∂Ψuu) From b1d87b48bdcc3a07a44fc2de3d07f03a105713a0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3?= Date: Fri, 31 Jul 2026 12:40:20 +0200 Subject: [PATCH 3/6] logreg --- src/PhysicalModels/MechanicalModels.jl | 6 +++--- src/TensorAlgebra/Functions.jl | 30 +++++++++++++++++++++++--- src/TensorAlgebra/TensorAlgebra.jl | 2 ++ 3 files changed, 32 insertions(+), 6 deletions(-) diff --git a/src/PhysicalModels/MechanicalModels.jl b/src/PhysicalModels/MechanicalModels.jl index 4499837..76e07b9 100644 --- a/src/PhysicalModels/MechanicalModels.jl +++ b/src/PhysicalModels/MechanicalModels.jl @@ -54,9 +54,9 @@ function (obj::CoerciveVolumetric)(::Float64=1.0) κ = obj.κ J(F) = det(F) H(F) = det(F) * inv(F)' - Ψ(F) = (κ / 4) * (J(F)^2 - 1 - 2log(J(F))) - ∂Ψ_∂J(F) = (κ / 2) * (J(F) - 1/J(F)) - ∂Ψ2_∂J2(F) = (κ / 2) * (1 + 1/J(F)^2) + Ψ(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) diff --git a/src/TensorAlgebra/Functions.jl b/src/TensorAlgebra/Functions.jl index 28e262c..57e5c0f 100644 --- a/src/TensorAlgebra/Functions.jl +++ b/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 diff --git a/src/TensorAlgebra/TensorAlgebra.jl b/src/TensorAlgebra/TensorAlgebra.jl index 97c079b..abe28e6 100644 --- a/src/TensorAlgebra/TensorAlgebra.jl +++ b/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 From de8fd953e7d0c2f90d1b2acc11f70b4be28a4c85 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3?= Date: Fri, 31 Jul 2026 13:02:09 +0200 Subject: [PATCH 4/6] update tolerance --- test/TestConstitutiveModels/PhysicalModelTests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/TestConstitutiveModels/PhysicalModelTests.jl b/test/TestConstitutiveModels/PhysicalModelTests.jl index fb45314..1d53a46 100644 --- a/test/TestConstitutiveModels/PhysicalModelTests.jl +++ b/test/TestConstitutiveModels/PhysicalModelTests.jl @@ -542,7 +542,7 @@ end # 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)) + test_derivatives_3D_(model, Kinematics(Mechano, Solid), rtol=1e-12) test_equilibrium_at_rest_3D(model) end From 833e2dfdfdc47312958539e07e45fc377451459a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3?= Date: Fri, 31 Jul 2026 15:26:09 +0200 Subject: [PATCH 5/6] Fixed test, added test --- test/TestTensorAlgebra/TensorAlgebraTests.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/test/TestTensorAlgebra/TensorAlgebraTests.jl b/test/TestTensorAlgebra/TensorAlgebraTests.jl index 857a605..fd021a7 100644 --- a/test/TestTensorAlgebra/TensorAlgebraTests.jl +++ b/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 From 73844aaa2671a73eb1bc1adf05ac451cbfd11cd2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3?= Date: Fri, 31 Jul 2026 15:45:16 +0200 Subject: [PATCH 6/6] The previous tests was doing nothing --- test/TestConstitutiveModels/PhysicalModelTests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/TestConstitutiveModels/PhysicalModelTests.jl b/test/TestConstitutiveModels/PhysicalModelTests.jl index 1d53a46..778614e 100644 --- a/test/TestConstitutiveModels/PhysicalModelTests.jl +++ b/test/TestConstitutiveModels/PhysicalModelTests.jl @@ -533,7 +533,7 @@ end # 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(λ=1.0) - test_derivatives_3D_(model, Kinematics(Mechano, Solid)) + test_derivatives_3D_(model, Kinematics(Mechano, Solid), rtol=1e-12) test_equilibrium_at_rest_3D(model) end