From 9037778e53806855ee067db2ed8a2a38444a4eb2 Mon Sep 17 00:00:00 2001 From: miguelmaso Date: Sun, 5 Apr 2026 12:29:41 +0200 Subject: [PATCH 01/16] Added softening law --- src/Exports.jl | 1 + src/PhysicalModels/ThermoMechanicalModels.jl | 14 ++++++++++++++ .../TestConstitutiveModels/ThermalLawsTests.jl | 18 ++++++++++++++++++ 3 files changed, 33 insertions(+) diff --git a/src/Exports.jl b/src/Exports.jl index 3818ab2..2272ad5 100644 --- a/src/Exports.jl +++ b/src/Exports.jl @@ -96,6 +96,7 @@ end @publish PhysicalModels ThermalLaw @publish PhysicalModels VolumetricLaw @publish PhysicalModels DeviatoricLaw +@publish PhysicalModels SofteningLaw @publish PhysicalModels EntropicMeltingLaw @publish PhysicalModels LogisticLaw @publish PhysicalModels InterceptLaw diff --git a/src/PhysicalModels/ThermoMechanicalModels.jl b/src/PhysicalModels/ThermoMechanicalModels.jl index 9c07e22..9116847 100644 --- a/src/PhysicalModels/ThermoMechanicalModels.jl +++ b/src/PhysicalModels/ThermoMechanicalModels.jl @@ -92,6 +92,20 @@ function derivatives(law::DeviatoricLaw) return (f, ∂f, ∂∂f) end +struct SofteningLaw <: ThermalLaw + θr::Float64 + θt::Float64 + γ::Float64 +end + +function derivatives(law::SofteningLaw) + @unpack θr, θt, γ = law + f(θ) = exp((θr/θt)^γ-(θ/θt)^γ) + ∂f(θ) = -γ/θt * (θ/θt)^(γ-1) * f(θ) + ∂∂f(θ) = 1/θ * (γ -1 -γ*(θ/θt)^γ) * ∂f(θ) + return (f, ∂f, ∂∂f) +end + struct InterceptLaw <: ThermalLaw θr::Float64 γ::Float64 diff --git a/test/TestConstitutiveModels/ThermalLawsTests.jl b/test/TestConstitutiveModels/ThermalLawsTests.jl index cfaba98..85add8a 100644 --- a/test/TestConstitutiveModels/ThermalLawsTests.jl +++ b/test/TestConstitutiveModels/ThermalLawsTests.jl @@ -10,3 +10,21 @@ using Test @test isapprox(ddf(θ), ForwardDiff.derivative(df, θ), rtol=1e-10) end end + +@testset "TrigonometricLaw" begin + law = TrigonometricLaw(273.15, 400.0) + f, df, ddf = derivatives(law) + for θ ∈ 200.0:50:400 + @test isapprox(df(θ), ForwardDiff.derivative(f, θ), rtol=1e-10) + @test isapprox(ddf(θ), ForwardDiff.derivative(df, θ), rtol=1e-10) + end +end + +@testset "SofteningLaw" begin + law = SofteningLaw(273.15, 300.0, 2.0) + f, df, ddf = derivatives(law) + for θ ∈ 200.0:50:400 + @test isapprox(df(θ), ForwardDiff.derivative(f, θ), rtol=1e-10) + @test isapprox(ddf(θ), ForwardDiff.derivative(df, θ), rtol=1e-10) + end +end From f34e9443ce3547262d2e1bf60e9c475208190086 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3?= Date: Tue, 7 Apr 2026 14:20:49 +0200 Subject: [PATCH 02/16] added parameter --- src/PhysicalModels/ThermoMechanicalModels.jl | 8 +++++--- test/TestConstitutiveModels/ThermalLawsTests.jl | 2 +- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/src/PhysicalModels/ThermoMechanicalModels.jl b/src/PhysicalModels/ThermoMechanicalModels.jl index 9116847..ca45008 100644 --- a/src/PhysicalModels/ThermoMechanicalModels.jl +++ b/src/PhysicalModels/ThermoMechanicalModels.jl @@ -96,12 +96,14 @@ struct SofteningLaw <: ThermalLaw θr::Float64 θt::Float64 γ::Float64 + δ::Float64 end function derivatives(law::SofteningLaw) - @unpack θr, θt, γ = law - f(θ) = exp((θr/θt)^γ-(θ/θt)^γ) - ∂f(θ) = -γ/θt * (θ/θt)^(γ-1) * f(θ) + @unpack θr, θt, γ, δ = law + h(θ) = exp((θr/θt)^γ-(θ/θt)^γ) * δ + f(θ) = h(θ) + 1 - δ + ∂f(θ) = -γ/θt * (θ/θt)^(γ-1) * h(θ) ∂∂f(θ) = 1/θ * (γ -1 -γ*(θ/θt)^γ) * ∂f(θ) return (f, ∂f, ∂∂f) end diff --git a/test/TestConstitutiveModels/ThermalLawsTests.jl b/test/TestConstitutiveModels/ThermalLawsTests.jl index 85add8a..4b8f244 100644 --- a/test/TestConstitutiveModels/ThermalLawsTests.jl +++ b/test/TestConstitutiveModels/ThermalLawsTests.jl @@ -21,7 +21,7 @@ end end @testset "SofteningLaw" begin - law = SofteningLaw(273.15, 300.0, 2.0) + law = SofteningLaw(273.15, 300.0, 2.0, 0.5) f, df, ddf = derivatives(law) for θ ∈ 200.0:50:400 @test isapprox(df(θ), ForwardDiff.derivative(f, θ), rtol=1e-10) From 14b32ff4ec674bf605725c1d3a6857a1ad6dc801 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3?= Date: Fri, 10 Apr 2026 11:11:02 +0200 Subject: [PATCH 03/16] thermal laws according to paper --- src/Exports.jl | 5 ++- src/PhysicalModels/ThermoMechanicalModels.jl | 37 ++++++++++++++++++- .../ThermalLawsTests.jl | 23 ++++++++---- 3 files changed, 54 insertions(+), 11 deletions(-) diff --git a/src/Exports.jl b/src/Exports.jl index 2272ad5..244da4d 100644 --- a/src/Exports.jl +++ b/src/Exports.jl @@ -96,8 +96,9 @@ end @publish PhysicalModels ThermalLaw @publish PhysicalModels VolumetricLaw @publish PhysicalModels DeviatoricLaw -@publish PhysicalModels SofteningLaw -@publish PhysicalModels EntropicMeltingLaw +@publish PhysicalModels EntropicElasticityLaw +@publish PhysicalModels NonlinearMeltingLaw +@publish PhysicalModels NonlinearSofteningLaw @publish PhysicalModels LogisticLaw @publish PhysicalModels InterceptLaw @publish PhysicalModels TrigonometricLaw diff --git a/src/PhysicalModels/ThermoMechanicalModels.jl b/src/PhysicalModels/ThermoMechanicalModels.jl index ca45008..7140196 100644 --- a/src/PhysicalModels/ThermoMechanicalModels.jl +++ b/src/PhysicalModels/ThermoMechanicalModels.jl @@ -65,13 +65,26 @@ function derivatives(law::VolumetricLaw) return (f, ∂f, ∂∂f) end -struct EntropicMeltingLaw <: ThermalLaw +struct EntropicElasticityLaw <: ThermalLaw + θr::Float64 + γ::Float64 +end + +function derivatives(law::EntropicElasticityLaw) + @unpack θr, γ = law + f(θ) = (θ/θr)^(γ+1) + ∂f(θ) = (γ+1) * θ^γ / θr^(γ+1) + ∂∂f(θ) = γ*(γ+1) * θ^(γ-1) / θr^(γ+1) + return (f, ∂f, ∂∂f) +end + +struct NonlinearMeltingLaw <: ThermalLaw θr::Float64 θM::Float64 γ::Float64 end -function derivatives(law::EntropicMeltingLaw) +function derivatives(law::NonlinearMeltingLaw) @unpack θr, θM, γ = law f(θ) = (1 - (θ/θM)^(γ+1)) / (1 - (θr/θM)^(γ+1)) ∂f(θ) = -(γ+1)*θ^γ/θM^(γ+1) / (1 - (θr/θM)^(γ+1)) @@ -79,6 +92,26 @@ function derivatives(law::EntropicMeltingLaw) return (f, ∂f, ∂∂f) end +struct NonlinearSofteningLaw <: ThermalLaw + θr::Float64 + θt::Float64 + γ::Float64 + δ::Float64 +end + +function derivatives(law::NonlinearSofteningLaw) + @unpack θr, θt, γ, δ = law + u(θ) = exp(-(θ/θt)^(γ+1)) + C = (1-δ) * u(θr) + δ + f(θ) = ((1-δ) * u(θ) + δ) / C + ∂f(θ) = -(1-δ)/C * (γ+1)/θt * (θ/θt)^γ * u(θ) + ∂∂f(θ) = (1-δ)/C * (γ+1)/θ^2 * (θ/θt)^(γ+1) * ((γ+1)*(θ/θt)^(γ+1)-γ) * u(θ) + return (f, ∂f, ∂∂f) +end + +@deprecate EntropicMeltingLaw NonlinearMeltingLaw +@deprecate SofteningLaw NonlinearSofteningLaw + struct DeviatoricLaw <: ThermalLaw θr::Float64 γ::Float64 diff --git a/test/TestConstitutiveModels/ThermalLawsTests.jl b/test/TestConstitutiveModels/ThermalLawsTests.jl index 4b8f244..35b6a2e 100644 --- a/test/TestConstitutiveModels/ThermalLawsTests.jl +++ b/test/TestConstitutiveModels/ThermalLawsTests.jl @@ -2,17 +2,26 @@ using ForwardDiff using HyperFEM using Test -@testset "LogisticLaw" begin - law = LogisticLaw(273.15, log(300.1), 0.11) +@testset "EntropicElasticityLaw" begin + law = EntropicElasticityLaw(273.15, 0.55) f, df, ddf = derivatives(law) - for θ ∈ 200.0:50:400 # NOTE: The numerical derivative of erf is a bad approximation, while the analyitical function uses the exact value. + for θ ∈ 200.0:50:400 @test isapprox(df(θ), ForwardDiff.derivative(f, θ), rtol=1e-3) @test isapprox(ddf(θ), ForwardDiff.derivative(df, θ), rtol=1e-10) end end -@testset "TrigonometricLaw" begin - law = TrigonometricLaw(273.15, 400.0) +@testset "NonlinearMeltingLaw" begin + law = NonlinearMeltingLaw(273.15, 400.0, 0.55) + f, df, ddf = derivatives(law) + for θ ∈ 200.0:50:400 + @test isapprox(df(θ), ForwardDiff.derivative(f, θ), rtol=1e-10) + @test isapprox(ddf(θ), ForwardDiff.derivative(df, θ), rtol=1e-10) + end +end + +@testset "NonlinearSofteningLaw" begin + law = NonlinearSofteningLaw(273.15, 300.0, 2.0, 0.5) f, df, ddf = derivatives(law) for θ ∈ 200.0:50:400 @test isapprox(df(θ), ForwardDiff.derivative(f, θ), rtol=1e-10) @@ -20,8 +29,8 @@ end end end -@testset "SofteningLaw" begin - law = SofteningLaw(273.15, 300.0, 2.0, 0.5) +@testset "TrigonometricLaw" begin + law = TrigonometricLaw(273.15, 400.0) f, df, ddf = derivatives(law) for θ ∈ 200.0:50:400 @test isapprox(df(θ), ForwardDiff.derivative(f, θ), rtol=1e-10) From 04a898c03d83d4ae78e6c7763743cd5356c0197b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3?= Date: Fri, 10 Apr 2026 11:45:17 +0200 Subject: [PATCH 04/16] removed unused thermal laws --- src/PhysicalModels/ThermoMechanicalModels.jl | 48 -------------------- 1 file changed, 48 deletions(-) diff --git a/src/PhysicalModels/ThermoMechanicalModels.jl b/src/PhysicalModels/ThermoMechanicalModels.jl index 7140196..0c1792f 100644 --- a/src/PhysicalModels/ThermoMechanicalModels.jl +++ b/src/PhysicalModels/ThermoMechanicalModels.jl @@ -125,36 +125,6 @@ function derivatives(law::DeviatoricLaw) return (f, ∂f, ∂∂f) end -struct SofteningLaw <: ThermalLaw - θr::Float64 - θt::Float64 - γ::Float64 - δ::Float64 -end - -function derivatives(law::SofteningLaw) - @unpack θr, θt, γ, δ = law - h(θ) = exp((θr/θt)^γ-(θ/θt)^γ) * δ - f(θ) = h(θ) + 1 - δ - ∂f(θ) = -γ/θt * (θ/θt)^(γ-1) * h(θ) - ∂∂f(θ) = 1/θ * (γ -1 -γ*(θ/θt)^γ) * ∂f(θ) - return (f, ∂f, ∂∂f) -end - -struct InterceptLaw <: ThermalLaw - θr::Float64 - γ::Float64 - δ::Float64 -end - -function derivatives(law::InterceptLaw) - @unpack θr, γ, δ = law - f(θ) = (θ/θr)^(-γ) * (1-δ) + δ - ∂f(θ) = -γ*θ^(-γ-1) * θr^γ * (1-δ) - ∂∂f(θ) = γ*(γ+1)*θ^(-γ-2) * θr^γ * (1-δ) - return (f, ∂f, ∂∂f) -end - struct TrigonometricLaw <: ThermalLaw θr::Float64 θM::Float64 @@ -186,24 +156,6 @@ function derivatives(law::PolynomialLaw) return (f, ∂f, ∂∂f) end -struct LogisticLaw <: ThermalLaw - θr::Float64 - μ::Float64 - σ::Float64 -end - -function derivatives(law::LogisticLaw) - @unpack θr, μ, σ = law - z(x) = (log(x) - μ) / σ - std_pdf(x) = 1/(σ*sqrt(2 * π)) * exp(-z(x)^2 / 2) - std_cdf(x) = 0.5 * (1 + erf(z(x) / sqrt(2))) - ξR = 1 / (1-std_cdf(θr)) - f(θ) = ξR * (1-std_cdf(θ)) - ∂f(θ) = -ξR / θ * std_pdf(θ) - ∂∂f(θ) = ξR / θ^2 * std_pdf(θ) * (1 + z(θ)/σ) - return (f, ∂f, ∂∂f) -end - struct ThermoMech_Bonet{T<:Thermo,M<:Mechano} <: ThermoMechano{T,M} thermo::T mechano::M From 1cd70fd732936ae9bdec3548cbe7df4f5c5a457c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3?= Date: Fri, 10 Apr 2026 11:46:58 +0200 Subject: [PATCH 05/16] removed erf and SpecialFunctions.jl --- src/PhysicalModels/PhysicalModels.jl | 1 - src/TensorAlgebra/Functions.jl | 22 -------------------- test/Project.toml | 1 - test/TestTensorAlgebra/TensorAlgebraTests.jl | 7 ------- 4 files changed, 31 deletions(-) diff --git a/src/PhysicalModels/PhysicalModels.jl b/src/PhysicalModels/PhysicalModels.jl index 843fc9b..9d1f562 100644 --- a/src/PhysicalModels/PhysicalModels.jl +++ b/src/PhysicalModels/PhysicalModels.jl @@ -11,7 +11,6 @@ using StaticArrays using ..TensorAlgebra using ..TensorAlgebra: _∂H∂F_2D using ..TensorAlgebra: trAA -using ..TensorAlgebra: erf import Base: + import Gridap: update_state! diff --git a/src/TensorAlgebra/Functions.jl b/src/TensorAlgebra/Functions.jl index 013b29c..36e3f5b 100644 --- a/src/TensorAlgebra/Functions.jl +++ b/src/TensorAlgebra/Functions.jl @@ -22,28 +22,6 @@ function logreg(J; Threshold=0.01) end -""" -Fast and dependency-free implementation of erf function, up to 1e-6 precision. -""" -@inline function erf(x::Real) - p = 0.3275911 - a1 = 0.254829592 - a2 = -0.284496736 - a3 = 1.421413741 - a4 = -1.453152027 - a5 = 1.061405429 - - ax = abs(x) - t = 1.0 / (1.0 + p*ax) - - y = (((((a5*t + a4)*t + a3)*t + a2)*t + a1)*t) - - r = 1.0 - y*exp(-ax*ax) - - return copysign(r, x) -end - - function _∂H∂F_2D() TensorValue(0.0, 0.0, 0.0, 1.0, 0.0, 0.0, -1.0, 0.0, 0.0, -1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0) end diff --git a/test/Project.toml b/test/Project.toml index 1e56ebb..d898967 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -3,7 +3,6 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" Gridap = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e" GridapSolvers = "6d3209ee-5e3c-4db7-a716-942eb12ed534" JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" -SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" diff --git a/test/TestTensorAlgebra/TensorAlgebraTests.jl b/test/TestTensorAlgebra/TensorAlgebraTests.jl index b0b191d..d38df12 100644 --- a/test/TestTensorAlgebra/TensorAlgebraTests.jl +++ b/test/TestTensorAlgebra/TensorAlgebraTests.jl @@ -115,10 +115,3 @@ end cofA = det(A) * inv(A') @test isapprox(cof(A), cofA) end - - -@testset "erf" begin - for x ∈ 0:.17234:2 - @test isapprox(TensorAlgebra.erf(x), SpecialFunctions.erf(x), atol=1e-6) - end -end From 262ae457c42db4ff9e792ae8280729718fe30af2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3?= Date: Mon, 13 Apr 2026 09:30:10 +0200 Subject: [PATCH 06/16] hotfix: removed pkg not found --- test/TestTensorAlgebra/TensorAlgebraTests.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/test/TestTensorAlgebra/TensorAlgebraTests.jl b/test/TestTensorAlgebra/TensorAlgebraTests.jl index d38df12..c84ecff 100644 --- a/test/TestTensorAlgebra/TensorAlgebraTests.jl +++ b/test/TestTensorAlgebra/TensorAlgebraTests.jl @@ -1,7 +1,6 @@ using Gridap.TensorValues using Gridap.Arrays using HyperFEM.TensorAlgebra -using SpecialFunctions using Test From 4f3830d597bcd8dc41c46934e12f727194153d7c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3?= Date: Mon, 13 Apr 2026 09:46:34 +0200 Subject: [PATCH 07/16] Hotfix: removed undef var --- src/Exports.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/Exports.jl b/src/Exports.jl index 244da4d..42cd037 100644 --- a/src/Exports.jl +++ b/src/Exports.jl @@ -99,7 +99,6 @@ end @publish PhysicalModels EntropicElasticityLaw @publish PhysicalModels NonlinearMeltingLaw @publish PhysicalModels NonlinearSofteningLaw -@publish PhysicalModels LogisticLaw @publish PhysicalModels InterceptLaw @publish PhysicalModels TrigonometricLaw @publish PhysicalModels PolynomialLaw From a1073b4a2e9fc9a2f49fbea075f7c07f59f96375 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3?= Date: Mon, 13 Apr 2026 10:06:34 +0200 Subject: [PATCH 08/16] Hotfix: removed undef var --- src/Exports.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/Exports.jl b/src/Exports.jl index 42cd037..e9f8dea 100644 --- a/src/Exports.jl +++ b/src/Exports.jl @@ -99,7 +99,6 @@ end @publish PhysicalModels EntropicElasticityLaw @publish PhysicalModels NonlinearMeltingLaw @publish PhysicalModels NonlinearSofteningLaw -@publish PhysicalModels InterceptLaw @publish PhysicalModels TrigonometricLaw @publish PhysicalModels PolynomialLaw From 52d0d3b85837ff5eb14c99c09eaeea278d7d9938 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3?= Date: Wed, 15 Apr 2026 18:47:56 +0200 Subject: [PATCH 09/16] addde isochoric neo-hookean --- src/Exports.jl | 1 + src/PhysicalModels/MechanicalModels.jl | 44 +++++++++++++++++++ src/PhysicalModels/PhysicalModels.jl | 1 + .../PhysicalModelTests.jl | 16 +++++++ 4 files changed, 62 insertions(+) diff --git a/src/Exports.jl b/src/Exports.jl index e9f8dea..0a22da1 100644 --- a/src/Exports.jl +++ b/src/Exports.jl @@ -30,6 +30,7 @@ end @publish PhysicalModels Yeoh3D @publish PhysicalModels Gent2D @publish PhysicalModels NeoHookean3D +@publish PhysicalModels IsochoricNeoHookean3D @publish PhysicalModels IncompressibleNeoHookean3D @publish PhysicalModels IncompressibleNeoHookean2D @publish PhysicalModels IncompressibleNeoHookean2D_CV diff --git a/src/PhysicalModels/MechanicalModels.jl b/src/PhysicalModels/MechanicalModels.jl index 72489cf..ef0ced7 100644 --- a/src/PhysicalModels/MechanicalModels.jl +++ b/src/PhysicalModels/MechanicalModels.jl @@ -1026,6 +1026,50 @@ struct ARAP2D <: IsoElastic end +struct IsochoricNeoHookean3D <: IsoElastic + μ::Float64 +end + +function IsochoricNeoHookean3D(; μ::Real) + IsochoricNeoHookean3D(float(μ)) +end + +function (obj::IsochoricNeoHookean3D)() + Ψ(F) = obj.μ / 2 * (F⊙F * det(F)^(-2/3) - 3) + ∂Ψ∂F(F) = begin + μ = obj.μ + J = det(F) + Ic = F⊙F + obj.μ * J^(-2/3) * (F - 1/3*Ic*inv(F)') + end + ∂Ψ∂FF(F) = begin + μ = obj.μ + J = det(F) + Ic = F⊙F + invF = inv(F) + H = cof(F) + TensorValue(ForwardDiff.jacobian(∂Ψ∂F, get_array(F))) + end +end + +function SecondPiola(obj::IsochoricNeoHookean3D) + Ψ(C) = obj.μ / 2 * (tr(C) * det(C)^(-1/3) - 3) + S(C) = begin + J2 = det(C) + invC = inv(C) + obj.μ * J2^(-1 / 3) * I3 - obj.μ / 3 * tr(C) * J2^(-1 / 3) * invC + end + ∂S∂C(C) = begin + J2 = det(C) + trC = tr(C) + invC = inv(C) + IinvC = I3 ⊗ invC + 1 / 3 * obj.μ * J2^(-1 / 3) * (4 / 3 * trC * invC ⊗ invC - (IinvC + IinvC') - trC / J2 * ×ᵢ⁴(C)) + end + return (Ψ, S, ∂S∂C) +end + + struct IncompressibleNeoHookean3D_2dP <: Mechano μ::Float64 τ::Float64 diff --git a/src/PhysicalModels/PhysicalModels.jl b/src/PhysicalModels/PhysicalModels.jl index 9d1f562..815e993 100644 --- a/src/PhysicalModels/PhysicalModels.jl +++ b/src/PhysicalModels/PhysicalModels.jl @@ -18,6 +18,7 @@ import Gridap: update_state! export Yeoh3D export Gent2D export NeoHookean3D +export IsochoricNeoHookean3D export IncompressibleNeoHookean3D export IncompressibleNeoHookean2D export IncompressibleNeoHookean2D_CV diff --git a/test/TestConstitutiveModels/PhysicalModelTests.jl b/test/TestConstitutiveModels/PhysicalModelTests.jl index d0e61b1..b8666e6 100644 --- a/test/TestConstitutiveModels/PhysicalModelTests.jl +++ b/test/TestConstitutiveModels/PhysicalModelTests.jl @@ -51,6 +51,14 @@ function test_equilibrium_at_rest_3D(obj::Mechano, atol=1e-10) @test isapprox(Ψ(I3), 0.0, atol=atol) end +function test_second_piola_3D_(model::PhysicalModel; rtol=1e-12, kwargs...) + F = I3 + ∇u3 + C = F'·F + Ψ, S, ∂S∂C = SecondPiola(model) + @test isapprox(S, 2*TensorValue(ForwardDiff.gradient(Ψ, get_array(C))), rtol=rtol, kwargs...) + @test isapprox(∂S∂C, TensorValue(ForwardDiff.hessian(S, get_array(C))), rtol=rtol, kwargs...) +end + @@ -204,6 +212,14 @@ end end +@testset "IsochoricNeoHookean3D" begin + model = IsochoricNeoHookean3D(λ=3) + test_derivatives_3D_(model, Kinematics(Mechano,Solid)) + test_second_piola_3D_(model) + test_equilibrium_at_rest_3D(model) +end + + @testset "IncompressibleNeoHookean3D_2dP" begin # Memory estimate: 0 bytes, allocs estimate: 0. Ce = TensorValue(0.01 + 1.0, 0.02, 0.03, 0.04, 0.05 + 1.0, 0.06, 0.07, 0.08, 0.09 + 1.0) From e2d19206887165b5e6409cee22f1e822f9436c8f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3?= Date: Wed, 15 Apr 2026 18:49:07 +0200 Subject: [PATCH 10/16] fixed model initialization in test --- 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 b8666e6..6f13c4f 100644 --- a/test/TestConstitutiveModels/PhysicalModelTests.jl +++ b/test/TestConstitutiveModels/PhysicalModelTests.jl @@ -213,7 +213,7 @@ end @testset "IsochoricNeoHookean3D" begin - model = IsochoricNeoHookean3D(λ=3) + model = IsochoricNeoHookean3D(μ=3) test_derivatives_3D_(model, Kinematics(Mechano,Solid)) test_second_piola_3D_(model) test_equilibrium_at_rest_3D(model) From 9de1f343c8070483284d0898377d1e4f1cc65da9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3?= Date: Wed, 15 Apr 2026 19:03:53 +0200 Subject: [PATCH 11/16] added missing return --- src/PhysicalModels/MechanicalModels.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/PhysicalModels/MechanicalModels.jl b/src/PhysicalModels/MechanicalModels.jl index ef0ced7..2473eea 100644 --- a/src/PhysicalModels/MechanicalModels.jl +++ b/src/PhysicalModels/MechanicalModels.jl @@ -1050,6 +1050,7 @@ function (obj::IsochoricNeoHookean3D)() H = cof(F) TensorValue(ForwardDiff.jacobian(∂Ψ∂F, get_array(F))) end + return (Ψ, ∂Ψ∂F, ∂Ψ∂FF) end function SecondPiola(obj::IsochoricNeoHookean3D) From 18e0ab1ce59f47a33c22490ca1194bc8f8a870d8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3?= Date: Wed, 15 Apr 2026 19:13:47 +0200 Subject: [PATCH 12/16] hotfix test --- src/TensorAlgebra/TensorsDefinitions.jl | 6 +++--- test/TestConstitutiveModels/PhysicalModelTests.jl | 4 ++-- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/TensorAlgebra/TensorsDefinitions.jl b/src/TensorAlgebra/TensorsDefinitions.jl index 22bf40b..0438d6b 100644 --- a/src/TensorAlgebra/TensorsDefinitions.jl +++ b/src/TensorAlgebra/TensorsDefinitions.jl @@ -39,15 +39,15 @@ const I9 = I_(9) -function Id(A::VectorValue{2, Float64}) +function Id(::VectorValue{2, Float64}) return I2 end -function Id(A::VectorValue{3, Float64}) +function Id(::VectorValue{3, Float64}) return I3 end -function Id(A::VectorValue{4, Float64}) +function Id(::VectorValue{4, Float64}) return I4 end diff --git a/test/TestConstitutiveModels/PhysicalModelTests.jl b/test/TestConstitutiveModels/PhysicalModelTests.jl index 6f13c4f..8b69ce4 100644 --- a/test/TestConstitutiveModels/PhysicalModelTests.jl +++ b/test/TestConstitutiveModels/PhysicalModelTests.jl @@ -55,8 +55,8 @@ function test_second_piola_3D_(model::PhysicalModel; rtol=1e-12, kwargs...) F = I3 + ∇u3 C = F'·F Ψ, S, ∂S∂C = SecondPiola(model) - @test isapprox(S, 2*TensorValue(ForwardDiff.gradient(Ψ, get_array(C))), rtol=rtol, kwargs...) - @test isapprox(∂S∂C, TensorValue(ForwardDiff.hessian(S, get_array(C))), rtol=rtol, kwargs...) + @test isapprox(S(C), 2*TensorValue(ForwardDiff.gradient(Ψ, get_array(C))), rtol=rtol, kwargs...) + @test isapprox(∂S∂C(C), TensorValue(ForwardDiff.hessian(S, get_array(C))), rtol=rtol, kwargs...) end From c9024f4c22c0aa7a7075b3c7b4aa255dd14464db Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3?= Date: Wed, 15 Apr 2026 19:33:47 +0200 Subject: [PATCH 13/16] partial fix --- test/TestConstitutiveModels/PhysicalModelTests.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/test/TestConstitutiveModels/PhysicalModelTests.jl b/test/TestConstitutiveModels/PhysicalModelTests.jl index 8b69ce4..c8f563e 100644 --- a/test/TestConstitutiveModels/PhysicalModelTests.jl +++ b/test/TestConstitutiveModels/PhysicalModelTests.jl @@ -14,6 +14,10 @@ import Base: - (-)(A::TensorValue, B::SMatrix) = get_array(A) - B +import Gridap.inner + +inner(a::SMatrix, b::SMatrix) = sum(get_array(a).data .* get_array(b).data) + const ∇u2 = TensorValue(1.0, 2.0, 3.0, 4.0) * 1e-3 const ∇u3 = TensorValue(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0) * 1e-3 const μParams = [6456.9137547089595, 896.4633794151492, @@ -214,7 +218,7 @@ end @testset "IsochoricNeoHookean3D" begin model = IsochoricNeoHookean3D(μ=3) - test_derivatives_3D_(model, Kinematics(Mechano,Solid)) + test_derivatives_3D_(model, Kinematics(Mechano,Solid), rtol=1e-12) test_second_piola_3D_(model) test_equilibrium_at_rest_3D(model) end From ee431959a9884f9f10ef4e36276e7d5df8dc4d94 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3?= Date: Wed, 15 Apr 2026 19:39:10 +0200 Subject: [PATCH 14/16] fixed method call --- test/TestConstitutiveModels/PhysicalModelTests.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/TestConstitutiveModels/PhysicalModelTests.jl b/test/TestConstitutiveModels/PhysicalModelTests.jl index c8f563e..badce86 100644 --- a/test/TestConstitutiveModels/PhysicalModelTests.jl +++ b/test/TestConstitutiveModels/PhysicalModelTests.jl @@ -59,8 +59,8 @@ function test_second_piola_3D_(model::PhysicalModel; rtol=1e-12, kwargs...) F = I3 + ∇u3 C = F'·F Ψ, S, ∂S∂C = SecondPiola(model) - @test isapprox(S(C), 2*TensorValue(ForwardDiff.gradient(Ψ, get_array(C))), rtol=rtol, kwargs...) - @test isapprox(∂S∂C(C), TensorValue(ForwardDiff.hessian(S, get_array(C))), rtol=rtol, kwargs...) + @test isapprox(S(C), 2*TensorValue(ForwardDiff.gradient(Ψ, get_array(C))), rtol=rtol, kwargs...) + @test isapprox(∂S∂C(C), TensorValue(ForwardDiff.jacobian(S, get_array(C))), rtol=rtol, kwargs...) end From 49f9468705e53a9a39b5480394148e88bc09f1ff Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3?= Date: Wed, 15 Apr 2026 19:53:31 +0200 Subject: [PATCH 15/16] try fix tensor algebra --- src/PhysicalModels/MechanicalModels.jl | 38 +++++++++---------- .../PhysicalModelTests.jl | 8 ++-- 2 files changed, 22 insertions(+), 24 deletions(-) diff --git a/src/PhysicalModels/MechanicalModels.jl b/src/PhysicalModels/MechanicalModels.jl index 2473eea..50b6d20 100644 --- a/src/PhysicalModels/MechanicalModels.jl +++ b/src/PhysicalModels/MechanicalModels.jl @@ -1054,18 +1054,18 @@ function (obj::IsochoricNeoHookean3D)() end function SecondPiola(obj::IsochoricNeoHookean3D) - Ψ(C) = obj.μ / 2 * (tr(C) * det(C)^(-1/3) - 3) - S(C) = begin - J2 = det(C) - invC = inv(C) - obj.μ * J2^(-1 / 3) * I3 - obj.μ / 3 * tr(C) * J2^(-1 / 3) * invC - end - ∂S∂C(C) = begin - J2 = det(C) - trC = tr(C) - invC = inv(C) - IinvC = I3 ⊗ invC - 1 / 3 * obj.μ * J2^(-1 / 3) * (4 / 3 * trC * invC ⊗ invC - (IinvC + IinvC') - trC / J2 * ×ᵢ⁴(C)) + μ = obj.μ + H(F) = cof(F) + Ψ(C) = μ / 2 * tr(C) * (det(C))^(-1 / 3) + ∂Ψ∂C(C) = μ / 2 * I3 * (det(C))^(-1 / 3) + ∂Ψ∂dC(C) = -μ / 6 * tr(C) * (det(C))^(-4 / 3) + S(C) = let HC = H(C) + 2 * (∂Ψ∂C(C) + ∂Ψ∂dC(C) * HC) + end + ∂2Ψ∂CdC(C) = -μ / 6 * I3 * (det(C))^(-4 / 3) + ∂2Ψ∂2dC(C) = 2 * μ / 9 * tr(C) * (det(C))^(-7 / 3) + ∂S∂C(C) = let HC = H(C) + 2 * (∂2Ψ∂2dC(C) * (HC ⊗ HC) + ∂2Ψ∂CdC(C) ⊗ HC + HC ⊗ ∂2Ψ∂CdC(C) + ∂Ψ∂dC(C) * ×ᵢ⁴(C)) end return (Ψ, S, ∂S∂C) end @@ -1087,16 +1087,14 @@ struct IncompressibleNeoHookean3D_2dP <: Mechano Ψ(Ce) = μ / 2 * tr(Ce) * (det(Ce))^(-1 / 3) ∂Ψ∂Ce(Ce) = μ / 2 * I3 * (det(Ce))^(-1 / 3) ∂Ψ∂dCe(Ce) = -μ / 6 * tr(Ce) * (det(Ce))^(-4 / 3) - Se(Ce) = - let HCe = H(Ce) - 2 * (∂Ψ∂Ce(Ce) + ∂Ψ∂dCe(Ce) * HCe) - end + Se(Ce) = let HCe = H(Ce) + 2 * (∂Ψ∂Ce(Ce) + ∂Ψ∂dCe(Ce) * HCe) + end ∂2Ψ∂CedCe(Ce) = -μ / 6 * I3 * (det(Ce))^(-4 / 3) ∂2Ψ∂2dCe(Ce) = 2 * μ / 9 * tr(Ce) * (det(Ce))^(-7 / 3) - ∂Se∂Ce(Ce) = - let HCe = H(Ce) - 2 * (∂2Ψ∂2dCe(Ce) * (HCe ⊗ HCe) + ∂2Ψ∂CedCe(Ce) ⊗ HCe + HCe ⊗ ∂2Ψ∂CedCe(Ce) + ∂Ψ∂dCe(Ce) * ×ᵢ⁴(Ce)) - end + ∂Se∂Ce(Ce) = let HCe = H(Ce) + 2 * (∂2Ψ∂2dCe(Ce) * (HCe ⊗ HCe) + ∂2Ψ∂CedCe(Ce) ⊗ HCe + HCe ⊗ ∂2Ψ∂CedCe(Ce) + ∂Ψ∂dCe(Ce) * ×ᵢ⁴(Ce)) + end return (Ψ, Se, ∂Se∂Ce) end diff --git a/test/TestConstitutiveModels/PhysicalModelTests.jl b/test/TestConstitutiveModels/PhysicalModelTests.jl index badce86..5cd966f 100644 --- a/test/TestConstitutiveModels/PhysicalModelTests.jl +++ b/test/TestConstitutiveModels/PhysicalModelTests.jl @@ -9,14 +9,14 @@ using HyperFEM.IO import Base: - - (-)(A::SMatrix, B::TensorValue) = A - get_array(B) # NOTE: These functions are required for LinearElasticity to work with ForwardDiff (-)(A::TensorValue, B::SMatrix) = get_array(A) - B +import Gridap:inner +inner(a::SMatrix, b::SMatrix) = sum(get_array(a).data .* get_array(b).data) # This function is required to work with ForwardDiff -import Gridap.inner - -inner(a::SMatrix, b::SMatrix) = sum(get_array(a).data .* get_array(b).data) +import HyperFEM.TensorAlgebra:cof +cof(a::SMatrix) = det(a) * inv(a)' # This function is required to work with ForwardDiff const ∇u2 = TensorValue(1.0, 2.0, 3.0, 4.0) * 1e-3 const ∇u3 = TensorValue(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0) * 1e-3 From 670a7d18f7e0dd0a7252cc1555c8cd64f7362639 Mon Sep 17 00:00:00 2001 From: miguelmaso Date: Thu, 16 Apr 2026 08:58:17 +0200 Subject: [PATCH 16/16] missing + function --- test/TestConstitutiveModels/PhysicalModelTests.jl | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/test/TestConstitutiveModels/PhysicalModelTests.jl b/test/TestConstitutiveModels/PhysicalModelTests.jl index 5cd966f..c22bae9 100644 --- a/test/TestConstitutiveModels/PhysicalModelTests.jl +++ b/test/TestConstitutiveModels/PhysicalModelTests.jl @@ -8,15 +8,18 @@ using HyperFEM.TensorAlgebra using HyperFEM.IO -import Base: - -(-)(A::SMatrix, B::TensorValue) = A - get_array(B) # NOTE: These functions are required for LinearElasticity to work with ForwardDiff +import Base: +,- +(+)(A::SMatrix, B::TensorValue) = A + get_array(B) # + is required by SecondPiola to work with ForwardDiff +(+)(A::TensorValue, B::SMatrix) = get_array(A) + B +(-)(A::SMatrix, B::TensorValue) = A - get_array(B) # - is required by LinearElasticity to work with ForwardDiff (-)(A::TensorValue, B::SMatrix) = get_array(A) - B -import Gridap:inner -inner(a::SMatrix, b::SMatrix) = sum(get_array(a).data .* get_array(b).data) # This function is required to work with ForwardDiff +import Gridap: inner +inner(a::SMatrix, b::SMatrix) = sum(a.data .* b.data) # inner function is required by SecondPiola to work with ForwardDiff + +import HyperFEM.TensorAlgebra: cof +cof(a::SMatrix) = det(a) * inv(a)' # cof is required by SecondPiola to work with ForwardDiff -import HyperFEM.TensorAlgebra:cof -cof(a::SMatrix) = det(a) * inv(a)' # This function is required to work with ForwardDiff const ∇u2 = TensorValue(1.0, 2.0, 3.0, 4.0) * 1e-3 const ∇u3 = TensorValue(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0) * 1e-3