Skip to content

Updated tests for Mechanical models #2

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
merged 9 commits into from
Sep 8, 2025
Merged

Updated tests for Mechanical models #2

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
d8537ec
initial commit
miguelmaso Sep 5, 2025
77b9ef0
some mech tests
miguelmaso Sep 5, 2025
eeb64c3
more tests
miguelmaso Sep 5, 2025
5b21a7f
refactor linear elasticity
miguelmaso Sep 8, 2025
a07a41b
remove physics from tensor algebra
miguelmaso Sep 8, 2025
332e597
testing linear elasticity
miguelmaso Sep 8, 2025
373b55c
shorter isapprox
miguelmaso Sep 8, 2025
f7d2314
simplified mechanical models
miguelmaso Sep 8, 2025
c5ed163
removed nonsense test
miguelmaso Sep 8, 2025
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
35 changes: 13 additions & 22 deletions src/PhysicalModels/PhysicalModels.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -56,6 +56,7 @@ export MagnetoMechModel
export GeneralizedMaxwell
export ViscousIncompressible

export PhysicalModel
export Mechano
export Elasto
export Visco
Expand Down Expand Up @@ -135,7 +136,6 @@ struct HessianRegularization{A,B} <: Mechano
new{typeof(Mechano),typeof(Mechano.Kinematic)}(Mechano, δ, Mechano.Kinematic)
end


function (obj::HessianRegularization)(Λ::Float64=1.0)
Ψs, ∂Ψs, ∂2Ψs = obj.Mechano()
δ = obj.δ
Expand All @@ -146,14 +146,11 @@ struct HessianRegularization{A,B} <: Mechano
val = real(vecval.values)
TensorValue(vec * diagm(max.(δ, val)) * vec')
end

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



end


struct Hessian∇JRegularization{A,B} <: Mechano
Mechano::A
δ::Float64
Expand All @@ -163,7 +160,6 @@ struct Hessian∇JRegularization{A,B} <: Mechano
new{typeof(Mechano),typeof(Mechano.Kinematic)}(Mechano, δ, κ, Mechano.Kinematic)
end


function (obj::Hessian∇JRegularization)(Λ::Float64=1.0)
Ψs, ∂Ψs, ∂2Ψs = obj.Mechano()
_, H, J = get_Kinematics(obj.Mechano.Kinematic; Λ=Λ)
Expand All @@ -179,8 +175,6 @@ struct Hessian∇JRegularization{A,B} <: Mechano
val = real(vecval.values)
TensorValue(vec * diagm(max.(δ, val)) * vec')
end

# ∂2Ψ(F, Jh) = TensorValue(real(λ(F, Jh).vectors) * diagm(max.(δ, real(λ(F, Jh).values))) * real(λ(F, Jh).vectors)')
return (Ψ, ∂Ψ, ∂2Ψ)
end
end
Expand Down Expand Up @@ -502,16 +496,15 @@ struct LinearElasticity2D{A} <: Mechano

function (obj::LinearElasticity2D)(Λ::Float64=1.0)
λ, μ = obj.λ, obj.μ
I22 = I2()
∂Ψuu(F) = _δδ_μ_2D(μ) + _δδ_λ_2D(λ)
∂Ψu(F) = ∂Ψuu(F) ⊙ (F - I22)
Ψ(F) = 0.5 * (F - I22) ⊙ (∂Ψuu(F) ⊙ (F - I22))
ε(F) = 0.5(F + F') -I2()
∂Ψuu(F) = μ * (δᵢₖδⱼₗ2D + δᵢₗδⱼₖ2D) + λ * δᵢⱼδₖₗ2D
∂Ψu(F) = ∂Ψuu(F) ⊙ (F - I2())
Ψ(F) = μ * sum(ε(F).*ε(F)) + 0.5 * λ * tr(ε(F))^2
return (Ψ, ∂Ψu, ∂Ψuu)
end


end


mutable struct LinearElasticity3D{A} <: Mechano
λ::Float64
μ::Float64
Expand All @@ -523,17 +516,15 @@ mutable struct LinearElasticity3D{A} <: Mechano

function (obj::LinearElasticity3D)(Λ::Float64=1.0)
λ, μ = obj.λ, obj.μ
I33 = I3()
# ∂Ψuu(F) = _δδ_μ_3D(μ) + _δδ_λ_3D(λ)
ε(F) = 0.5(F + F') -I3()
∂Ψuu(F) = μ * (δᵢₖδⱼₗ3D + δᵢₗδⱼₖ3D) + λ * δᵢⱼδₖₗ3D
∂Ψu(F) = ∂Ψuu(F) ⊙ (F - I33)
Ψ(F) = 0.5 * (F - I33) ⊙ (∂Ψuu(F) ⊙ (F - I33))
∂Ψu(F) = ∂Ψuu(F) ⊙ (F - I3())
Ψ(F) = μ * sum(ε(F).*ε(F)) + 0.5 * λ * tr(ε(F))^2
return (Ψ, ∂Ψu, ∂Ψuu)
end


end


struct VolumetricEnergy{A} <: Mechano
λ::Float64
Kinematic::A
Expand Down Expand Up @@ -794,11 +785,11 @@ struct NonlinearMooneyRivlin_CV{A} <: Mechano



Ψ(F::TensorValue{3, 3, Float64, 9}) = μ1 / (2.0 * α * 3.0^(α - 1)) * (tr((F)' * F))^α +
Ψ(F) = μ1 / (2.0 * α * 3.0^(α - 1)) * (tr((F)' * F))^α +
μ2 / (2.0 * β * 3.0^(β - 1)) * (tr((H(F))' * H(F)))^β -
(μ1 + 2*μ2) * log(J(F)) + λ * (J(F)^(γ) + J(F)^(-γ))

∂Ψ_∂F(F::TensorValue{3, 3, Float64, 9}) = ((μ1 / (3.0^(α - 1)) * (trAA(F))^(α - 1))) * F
∂Ψ_∂F(F) = ((μ1 / (3.0^(α - 1)) * (trAA(F))^(α - 1))) * F
∂Ψ_∂H(F) = ((μ2 / (3.0^(β - 1)) * (tr((H(F))' * H(F)))^(β - 1))) * H(F)
∂Ψ_∂J(F) = -(μ1 + 2*μ2) * (1.0 / J(F)) + λ * γ * (J(F)^(γ - 1) - J(F)^(-γ - 1))
∂Ψu(F) = ∂Ψ_∂F(F) + ∂Ψ_∂H(F) × F + ∂Ψ_∂J(F) * H(F)
Expand Down
220 changes: 8 additions & 212 deletions src/TensorAlgebra/TensorAlgebra.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,9 @@ export I2
export I4
export logreg
export Tensorize
export δᵢⱼδₖₗ2D
export δᵢₖδⱼₗ2D
export δᵢₗδⱼₖ2D
export δᵢⱼδₖₗ3D
export δᵢₖδⱼₗ3D
export δᵢₗδⱼₖ3D
Expand Down Expand Up @@ -90,9 +93,13 @@ function _Kroneckerδδ(δδ::Function, N::Int)
i, j, k, l = _full_idx4(α,N)
δδ(i,j,k,l) ? 1.0 : 0.0
end,
9*9))
N*N*N*N))
end

const δᵢⱼδₖₗ2D = _Kroneckerδδ((i,j,k,l) -> i==j && k==l, 2)
const δᵢₖδⱼₗ2D = _Kroneckerδδ((i,j,k,l) -> i==k && j==l, 2)
const δᵢₗδⱼₖ2D = _Kroneckerδδ((i,j,k,l) -> i==l && j==k, 2)

const δᵢⱼδₖₗ3D = _Kroneckerδδ((i,j,k,l) -> i==j && k==l, 3)
const δᵢₖδⱼₗ3D = _Kroneckerδδ((i,j,k,l) -> i==k && j==l, 3)
const δᵢₗδⱼₖ3D = _Kroneckerδδ((i,j,k,l) -> i==l && j==k, 3)
Expand All @@ -102,217 +109,6 @@ 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

function _δδ_μ_2D(μ::Float64)
TensorValue{4,4,Float64,16}(
2 * μ,
0.0,
0.0,
0.0,
0.0,
μ,
μ,
0.0,
0.0,
μ,
μ,
0.0,
0.0,
0.0,
0.0,
2.0 * μ)
end

function _δδ_λ_2D(λ::Float64)
TensorValue{4,4,Float64,16}(
λ,
0.0,
0.0,
λ,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
λ,
0.0,
0.0,
λ)
end

function _δδ_μ_3D(μ::Float64)
TensorValue{9,9,Float64,81}(
2.0 * μ,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
μ,
0.0,
μ,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
μ,
0.0,
0.0,
0.0,
μ,
0.0,
0.0,
0.0,
μ,
0.0,
μ,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
2.0 * μ,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
μ,
0.0,
μ,
0.0,
0.0,
0.0,
μ,
0.0,
0.0,
0.0,
μ,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
μ,
0.0,
μ,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
2.0 * μ)
end

function _δδ_λ_3D(λ::Float64)
TensorValue{9,9,Float64,81}(
λ,
0.0,
0.0,
0.0,
λ,
0.0,
0.0,
0.0,
λ,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
λ,
0.0,
0.0,
0.0,
λ,
0.0,
0.0,
0.0,
λ,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
λ,
0.0,
0.0,
0.0,
λ,
0.0,
0.0,
0.0,
λ)
end



function Gridap.TensorValues.outer(A::TensorValue{D,D,Float64}, B::TensorValue{D,D,Float64}) where {D}
return (A ⊗₁₂³⁴ B)
Expand Down
Loading