Skip to content

Adds HGO_4Fibers anisotropic model #42

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
jmartfrut merged 1 commit into from
Nov 2, 2025
Merged

Adds HGO_4Fibers anisotropic model #42

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion src/Exports.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -64,6 +64,7 @@ end
@publish PhysicalModels Hessian∇JRegularization
@publish PhysicalModels ViscousIncompressible
@publish PhysicalModels GeneralizedMaxwell
@publish PhysicalModels HGO_4Fibers

@publish PhysicalModels Mechano
@publish PhysicalModels Thermo
Expand All @@ -81,7 +82,7 @@ end
@publish PhysicalModels EvolutiveKinematics
@publish PhysicalModels get_Kinematics
@publish PhysicalModels getIsoInvariants

@publish PhysicalModels initializeStateVariables
@publish PhysicalModels updateStateVariables!

Expand Down
7 changes: 7 additions & 0 deletions src/PhysicalModels/ElectroMechanicalModels.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -81,6 +81,13 @@ function _getCoupling(elec::Electro, mec::Mechano, Λ::Float64)
end


function (+)(Model1::Electro, Model2::Mechano)
ElectroMechModel(Model1, Model2)
end
function (+)(Model1::Mechano, Model2::Electro)
ElectroMechModel(Model2, Model1)
end

struct FlexoElectroModel{EM<:ElectroMechano} <: FlexoElectro
electromechano::EM
κ::Float64
Expand Down
8 changes: 8 additions & 0 deletions src/PhysicalModels/MagnetoMechanicalModels.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -202,3 +202,11 @@ function _getCoupling(mag::HardMagnetic2D, mec::Mechano, Λ::Float64)

return (Ψ, ∂Ψ_u, ∂Ψ_φ, ∂Ψ_uu, ∂Ψ_φu, ∂Ψ_φφ)
end


function (+)(Model1::Magneto, Model2::Mechano)
MagnetoMechModel(Model1, Model2)
end
function (+)(Model1::Mechano, Model2::Magneto)
MagnetoMechModel(Model2, Model1)
end
76 changes: 56 additions & 20 deletions src/PhysicalModels/MechanicalModels.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,3 @@

# ============================================
# Regularization of Mechanical models
# ============================================
Expand Down Expand Up @@ -99,45 +98,45 @@ struct ComposedIsoElastic <: IsoElastic
Model1::IsoElastic
Model2::IsoElastic
Kinematic
function ComposedIsoElastic(model1::IsoElastic,model2::IsoElastic)
new(model1,model2,model1.Kinematic)
function ComposedIsoElastic(model1::IsoElastic, model2::IsoElastic)
new(model1, model2, model1.Kinematic)
end
function (obj::ComposedIsoElastic)(Λ::Float64=1.0)
Ψ1, ∂Ψu1, ∂Ψuu1 = obj.Model1(Λ)
Ψ2, ∂Ψu2, ∂Ψuu2 = obj.Model2(Λ)
Ψ(x)=Ψ1(x)+Ψ2(x)
∂Ψ(x)=∂Ψu1(x)+∂Ψu2(x)
∂∂Ψ(x)=∂Ψuu1(x)+∂Ψuu2(x)
Ψ1, ∂Ψu1, ∂Ψuu1 = obj.Model1(Λ)
Ψ2, ∂Ψu2, ∂Ψuu2 = obj.Model2(Λ)
Ψ(x) = Ψ1(x) + Ψ2(x)
∂Ψ(x) = ∂Ψu1(x) + ∂Ψu2(x)
∂∂Ψ(x) = ∂Ψuu1(x) + ∂Ψuu2(x)
return (Ψ, ∂Ψ, ∂∂Ψ)
end
end


function (+)(Model1::IsoElastic, Model2::IsoElastic)
ComposedIsoElastic(Model1,Model2)
ComposedIsoElastic(Model1, Model2)
end

struct ComposedAnisoElastic <: AnisoElastic
Model1::IsoElastic
Model2::AnisoElastic
Kinematic
function ComposedAnisoElastic(model1::IsoElastic,model2::AnisoElastic)
new(model1,model2,model1.Kinematic)
function ComposedAnisoElastic(model1::IsoElastic, model2::AnisoElastic)
new(model1, model2, model1.Kinematic)
end
function (obj::ComposedAnisoElastic)(Λ::Float64=1.0)
DΨ1 = obj.Model1(Λ)
DΨ2 = obj.Model2(Λ)
Ψ, ∂Ψ, ∂∂Ψ = map((ψ1,ψ2) -> (x,N) -> ψ1(x) + ψ2(x,N), DΨ1, DΨ2)
Ψ, ∂Ψ, ∂∂Ψ = map((ψ1, ψ2) -> (x, N...) -> ψ1(x) + ψ2(x, N...), DΨ1, DΨ2)
return (Ψ, ∂Ψ, ∂∂Ψ)
end
end
end


function (+)(Model1::IsoElastic, Model2::AnisoElastic)
ComposedAnisoElastic(Model1,Model2)
ComposedAnisoElastic(Model1, Model2)
end
function (+)(Model1::AnisoElastic, Model2::IsoElastic)
ComposedAnisoElastic(Model2,Model1)
ComposedAnisoElastic(Model2, Model1)
end

struct MultiAnisoElastic <: AnisoElastic
Expand All @@ -155,8 +154,8 @@ function (obj::MultiAnisoElastic)(args...)
∂Ψ∂FF(F, N) = mapreduce((∂Ψi∂FF, Ni) -> ∂Ψi∂FF(F, Ni), +, ∂Ψα∂FF, N)
(Ψ, ∂Ψ∂F, ∂Ψ∂FF)
end
transpose(x::NTuple{N, <:Tuple{<:Function, <:Function, <:Function}}) where N = map(i -> getindex.(x, i), 1:3)

transpose(x::NTuple{N,<:Tuple{<:Function,<:Function,<:Function}}) where N = map(i -> getindex.(x, i), 1:3)



Expand Down Expand Up @@ -707,6 +706,43 @@ struct TransverseIsotropy2D{A} <: AnisoElastic
end


struct HGO_4Fibers <: AnisoElastic
c1::Vector{Float64}
c2::Vector{Float64}
Kinematic::Kinematics{Mechano}
function HGO_4Fibers(; c1::Vector{Float64}, c2::Vector{Float64}, Kinematic::KinematicModel=Kinematics(Mechano))
@assert length(c1) == length(c2) == 4
new(c1, c2, Kinematic)
end

function (obj::HGO_4Fibers)(Λ::Float64=1.0; Threshold=0.01)
_, H, J = get_Kinematics(obj.Kinematic; Λ=Λ)
c1, c2 = obj.c1, obj.c2

Ψ(F, M1, M2, M3, M4) = c1[1] / (4 * c2[1]) * (exp(c2[1] * ((F * M1) ⋅ (F * M1) - 1.0)^2.0) - 1.0) +
c1[2] / (4 * c2[2]) * (exp(c2[2] * ((F * M2) ⋅ (F * M2) - 1.0)^2.0) - 1.0) +
c1[3] / (4 * c2[3]) * (exp(c2[3] * ((F * M3) ⋅ (F * M3) - 1.0)^2.0) - 1.0) +
c1[4] / (4 * c2[4]) * (exp(c2[4] * ((F * M4) ⋅ (F * M4) - 1.0)^2.0) - 1.0)

∂Ψ∂F(F, M1, M2, M3, M4) = c1[1] * exp(c2[1] * ((F * M1) ⋅ (F * M1) - 1.0)^2.0) * ((F * M1) ⋅ (F * M1) - 1.0) * ((F * M1) ⊗ M1) +
c1[2] * exp(c2[2] * ((F * M2) ⋅ (F * M2) - 1.0)^2.0) * ((F * M2) ⋅ (F * M2) - 1.0) * ((F * M2) ⊗ M2) +
c1[3] * exp(c2[3] * ((F * M3) ⋅ (F * M3) - 1.0)^2.0) * ((F * M3) ⋅ (F * M3) - 1.0) * ((F * M3) ⊗ M3) +
c1[4] * exp(c2[4] * ((F * M4) ⋅ (F * M4) - 1.0)^2.0) * ((F * M4) ⋅ (F * M4) - 1.0) * ((F * M4) ⊗ M4)


∂Ψ2∂F∂F(F, M1, M2, M3, M4) = c1[1] * exp(c2[1] * ((F * M1) ⋅ (F * M1) - 1.0)^2.0) * ((4 * c2[1] * (((F * M1) ⋅ (F * M1) - 1.0)^2.0) + 2.0) * (((F * M1) ⊗ M1) ⊗ ((F * M1) ⊗ M1)) + ((F * M1) ⋅ (F * M1) - 1.0) * (I3 ⊗₁₃²⁴ (M1 ⊗ M1))) +
c1[2] * exp(c2[2] * ((F * M2) ⋅ (F * M2) - 1.0)^2.0) * ((4 * c2[2] * (((F * M2) ⋅ (F * M2) - 1.0)^2.0) + 2.0) * (((F * M2) ⊗ M2) ⊗ ((F * M2) ⊗ M2)) + ((F * M2) ⋅ (F * M2) - 1.0) * (I3 ⊗₁₃²⁴ (M2 ⊗ M2))) +
c1[3] * exp(c2[3] * ((F * M3) ⋅ (F * M3) - 1.0)^2.0) * ((4 * c2[3] * (((F * M3) ⋅ (F * M3) - 1.0)^2.0) + 2.0) * (((F * M3) ⊗ M3) ⊗ ((F * M3) ⊗ M3)) + ((F * M3) ⋅ (F * M3) - 1.0) * (I3 ⊗₁₃²⁴ (M3 ⊗ M3))) +
c1[4] * exp(c2[4] * ((F * M4) ⋅ (F * M4) - 1.0)^2.0) * ((4 * c2[4] * (((F * M4) ⋅ (F * M4) - 1.0)^2.0) + 2.0) * (((F * M4) ⊗ M4) ⊗ ((F * M4) ⊗ M4)) + ((F * M4) ⋅ (F * M4) - 1.0) * (I3 ⊗₁₃²⁴ (M4 ⊗ M4)))


return (Ψ, ∂Ψ∂F, ∂Ψ2∂F∂F)
end
end




struct IncompressibleNeoHookean3D{A} <: IsoElastic
λ::Float64
μ::Float64
Expand Down
3 changes: 3 additions & 0 deletions src/PhysicalModels/PhysicalModels.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
module PhysicalModels

using Gridap
using Gridap.CellData
using Gridap.Helpers
using DrWatson

Expand Down Expand Up @@ -53,6 +54,7 @@ export ThermoElectroMech_Bonet
export MagnetoMechModel
export GeneralizedMaxwell
export ViscousIncompressible
export HGO_4Fibers

export PhysicalModel
export Mechano
Expand Down Expand Up @@ -108,6 +110,7 @@ abstract type ThermoMechano <: MultiPhysicalModel end
abstract type ThermoElectro <: MultiPhysicalModel end
abstract type FlexoElectro <: MultiPhysicalModel end
abstract type MagnetoMechano <: MultiPhysicalModel end
abstract type InternalFibers end

include("KinematicModels.jl")

Expand Down
1 change: 1 addition & 0 deletions src/PhysicalModels/ThermoMechanicalModels.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -107,3 +107,4 @@ function _getCoupling(term::Thermo, mec::Mechano, Λ::Float64)

return (Ψtm, ∂Ψtm_u, ∂Ψtm_θ, ∂Ψtm_uu, ∂Ψtm_uθ, ∂Ψtm_θθ)
end

6 changes: 3 additions & 3 deletions test/TestConstitutiveModels/ElectroMechanicalTests.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@ const ∇un = TensorValue(1.0:9.0...) * 5e-4
@testset "ElectroMechano" begin
modelMR = MooneyRivlin3D(λ=3.0, μ1=1.0, μ2=2.0)
modelID = IdealDielectric(ε=4.0)
modelelectro = ElectroMechModel(mechano=modelMR, electro=modelID)
modelelectro =modelMR+modelID
Ψ, ∂Ψu, ∂Ψφ, ∂Ψuu, ∂Ψφu, ∂Ψφφ = modelelectro()
F, _, _ = get_Kinematics(modelMR.Kinematic)
E = get_Kinematics(modelID.Kinematic)
Expand Down Expand Up @@ -66,7 +66,7 @@ end
viscous_branch1 = ViscousIncompressible(short_term, τ=6.)
visco_elastic = GeneralizedMaxwell(hyper_elastic, viscous_branch1)
dielectric = IdealDielectric(ε=1.0)
model = ElectroMechModel(dielectric, visco_elastic)
model =dielectric+visco_elastic
F, _, _ = get_Kinematics(model.mechano.Kinematic)
E = get_Kinematics(model.electro.Kinematic)
Uvn = TensorValue(1.,2.,3.,2.,4.,5.,3.,5.,6.) * 2e-4 + I3
Expand All @@ -86,7 +86,7 @@ end
viscous_branch2 = ViscousIncompressible(short_term, τ=60.)
visco_elastic = GeneralizedMaxwell(hyper_elastic, viscous_branch1, viscous_branch2)
dielectric = IdealDielectric(ε=1.0)
model = ElectroMechModel(dielectric, visco_elastic)
model =dielectric+visco_elastic
F, _, _ = get_Kinematics(model.mechano.Kinematic)
E = get_Kinematics(model.electro.Kinematic)
Uvn = TensorValue(1.,2.,3.,2.,4.,5.,3.,5.,6.) * 2e-4 + I3
Expand Down
61 changes: 58 additions & 3 deletions test/TestConstitutiveModels/PhysicalModelTests.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,61 @@ function test_derivatives_3D_(model::PhysicalModel; rtol=1e-14, kwargs...)
end


@testset "HGO_4Fibers" begin
c1 = [0.6639232500447778, 0.5532987701062146, 0.9912576142028674, 0.4951942011240962]
c2 = [0.800583033264982, 0.3141082734275339, 0.8063905248474006, 0.5850486948450955]
M1 = [ 0.36799630150742724, 0.8353476258002335, 0.57704047269419]
M1 = VectorValue(M1/norm(M1))
M2 = [ 0.3857610953527303, 0.024655018338846868, 0.6133770006613235]
M2 = VectorValue(M2/norm(M2))
M3 = [ 0.4516424464747618, 0.6609741557924332, 0.8441070681368911]
M3 = VectorValue(M3/norm(M3))
M4 = [0.7650638541897623, 0.8268625401770648, 0.3103412304991431]
M4 = VectorValue(M4/norm(M4))

model = HGO_4Fibers(c1=c1, c2=c2)
Ψ, ∂Ψ∂F, ∂Ψ∂F∂F = model()
∂Ψ∂F_(F,M1,M2,M3,M4) = TensorValue(ForwardDiff.gradient(F -> Ψ(F,get_array(M1),get_array(M2),get_array(M3),get_array(M4)), get_array(F)))
∂Ψ∂F∂F_(F,M1,M2,M3,M4) = TensorValue(ForwardDiff.hessian(F -> Ψ(F,get_array(M1),get_array(M2),get_array(M3),get_array(M4)), get_array(F)))


F, _, _ = get_Kinematics(model.Kinematic)
@test Ψ(F(∇u3), M1,M2,M3,M4) == 0.0005561006012767033
@test isapprox(norm(∂Ψ∂F(F(∇u3), M1,M2,M3,M4)) , norm(∂Ψ∂F_(F(∇u3), M1,M2,M3,M4)),rtol=1e-14)
@test isapprox(norm(∂Ψ∂F∂F(F(∇u3), M1,M2,M3,M4)), norm(∂Ψ∂F∂F_(F(∇u3), M1,M2,M3,M4)),rtol=1e-14)
end


@testset "HGO_4Fibers+MooneyRivlin3D" begin
c1 = [0.6639232500447778, 0.5532987701062146, 0.9912576142028674, 0.4951942011240962]
c2 = [0.800583033264982, 0.3141082734275339, 0.8063905248474006, 0.5850486948450955]
M1 = [ 0.36799630150742724, 0.8353476258002335, 0.57704047269419]
M1 = VectorValue(M1/norm(M1))
M2 = [ 0.3857610953527303, 0.024655018338846868, 0.6133770006613235]
M2 = VectorValue(M2/norm(M2))
M3 = [ 0.4516424464747618, 0.6609741557924332, 0.8441070681368911]
M3 = VectorValue(M3/norm(M3))
M4 = [0.7650638541897623, 0.8268625401770648, 0.3103412304991431]
M4 = VectorValue(M4/norm(M4))



model1 = MooneyRivlin3D(λ=(1e3 + 1e3) * 1e2, μ1=1e3, μ2=1e3)
model2 = HGO_4Fibers(c1=c1, c2=c2)
model=model1+model2

Ψ, ∂Ψ∂F, ∂Ψ∂F∂F = model()
∂Ψ∂F_(F,M1,M2,M3,M4) = TensorValue(ForwardDiff.gradient(F -> Ψ(F,get_array(M1),get_array(M2),get_array(M3),get_array(M4)), get_array(F)))
∂Ψ∂F∂F_(F,M1,M2,M3,M4) = TensorValue(ForwardDiff.hessian(F -> Ψ(F,get_array(M1),get_array(M2),get_array(M3),get_array(M4)), get_array(F)))


F, _, _ = get_Kinematics(model.Kinematic)
@test Ψ(F(∇u3), M1,M2,M3,M4) == 23.21318362353833
@test isapprox(norm(∂Ψ∂F(F(∇u3), M1,M2,M3,M4)) , norm(∂Ψ∂F_(F(∇u3), M1,M2,M3,M4)),rtol=1e-10)
@test isapprox(norm(∂Ψ∂F∂F(F(∇u3), M1,M2,M3,M4)), norm(∂Ψ∂F∂F_(F(∇u3), M1,M2,M3,M4)),rtol=1e-10)
end



@testset "Iso+Aniso" begin
N = VectorValue(1.0, 2.0, 3.0)
Expand Down Expand Up @@ -525,7 +580,7 @@ end

modelMR = MooneyRivlin3D(λ=3.0, μ1=1.0, μ2=2.0)
modelID = HardMagnetic(μ=1.2566e-6, αr=40e-3, χe=0.0, χr=8.0)
modelmagneto = MagnetoMechModel(modelID, modelMR)
modelmagneto=modelMR+modelID
Ψ, ∂Ψu, ∂Ψφ, ∂Ψuu, ∂Ψφu, ∂Ψφφ = modelmagneto()
F, _, _ = get_Kinematics(modelMR.Kinematic)
H0 = get_Kinematics(modelID.Kinematic)
Expand Down Expand Up @@ -606,7 +661,7 @@ end

modelMR = MooneyRivlin2D(λ=3.0, μ1=1.0, μ2=2.0)
modelID = IdealMagnetic2D(μ=1.2566e-6, χe=0.0)
modelmagneto = MagnetoMechModel(modelID, modelMR)
modelmagneto=modelMR+modelID
Ψ, ∂Ψu, ∂Ψφ, ∂Ψuu, ∂Ψφu, ∂Ψφφ = modelmagneto()
F, _, _ = get_Kinematics(modelMR.Kinematic)
H0 = get_Kinematics(modelID.Kinematic)
Expand Down Expand Up @@ -648,7 +703,7 @@ end

modelMR = MooneyRivlin2D(λ=3.0, μ1=1.0, μ2=2.0)
modelID = HardMagnetic2D(μ=1.2566e-6, αr=40e-3, χe=0.0, χr=8.0)
modelmagneto = MagnetoMechModel(modelID, modelMR)
modelmagneto=modelMR+modelID
Ψ, ∂Ψu, ∂Ψφ, ∂Ψuu, ∂Ψφu, ∂Ψφφ = modelmagneto()
F, _, _ = get_Kinematics(modelMR.Kinematic)
H0 = get_Kinematics(modelID.Kinematic)
Expand Down