From 3444f00e3cb95a8b1885c3d43287141daf26a215 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jesus=20Mart=C3=ADnez?= Date: Wed, 5 Nov 2025 19:28:24 +0100 Subject: [PATCH] Refactors kinematic descriptions in physical models Introduces a more explicit kinematic description and associates it directly with the physics of the solid material, thus decoupling the kinematic model from the specific electrical or magnetic model. This change simplifies the model definitions and usage in weak forms and post-processors, improving overall code clarity and maintainability. --- README.md | 8 +- src/ComputationalModels/PostProcessors.jl | 45 ++- src/Exports.jl | 1 + src/PhysicalModels/ElectricalModels.jl | 7 +- src/PhysicalModels/ElectroMechanicalModels.jl | 4 +- src/PhysicalModels/KinematicModels.jl | 44 ++- src/PhysicalModels/MagneticModels.jl | 33 +- src/PhysicalModels/MagnetoMechanicalModels.jl | 8 +- src/PhysicalModels/MechanicalModels.jl | 284 +++++++++--------- src/PhysicalModels/PINNs.jl | 4 +- src/PhysicalModels/PhysicalModels.jl | 4 +- .../ThermoElectroMechanicalModels.jl | 6 +- src/PhysicalModels/ThermoMechanicalModels.jl | 8 +- src/PhysicalModels/ViscousModels.jl | 27 +- src/WeakForms/ElectroMechanics.jl | 88 +++--- src/WeakForms/WeakForms.jl | 258 ++++++++-------- .../ElectroMechanicalTests.jl | 35 ++- .../PhysicalModelTests.jl | 159 ++++++---- .../ViscousModelsTests.jl | 4 - test/TestWeakForms/WeakForms.jl | 20 +- test/data/ViscoElasticSimulation.jl | 11 +- 21 files changed, 550 insertions(+), 508 deletions(-) diff --git a/README.md b/README.md index f5fb20b..2fedf38 100644 --- a/README.md +++ b/README.md @@ -92,9 +92,13 @@ Uφ = TrialFESpace(Vφ, D_bc.BoundaryCondition[2], 1.0) V = MultiFieldFESpace([Vu, Vφ]) U = MultiFieldFESpace([Uu, Uφ]) +# Kinematic Description +km=Kinematics(Mechano,Solid) +ke=Kinematics(Electro,Solid) + # residual and jacobian function of load factor -res(Λ) = ((u, φ), (v, vφ)) -> residual(physmodel, (u, φ), (v, vφ), dΩ) -jac(Λ) = ((u, φ), (du, dφ), (v, vφ)) -> jacobian(physmodel, (u, φ), (du, dφ), (v, vφ), dΩ) +res(Λ) = ((u, φ), (v, vφ)) -> residual(physmodel, (km,ke),(u, φ), (v, vφ), dΩ) +jac(Λ) = ((u, φ), (du, dφ), (v, vφ)) -> jacobian(physmodel, (km,ke), (u, φ), (du, dφ), (v,vφ),dΩ) # nonlinear solver ls = LUSolver() diff --git a/src/ComputationalModels/PostProcessors.jl b/src/ComputationalModels/PostProcessors.jl index b2afdce..eb18d1d 100644 --- a/src/ComputationalModels/PostProcessors.jl +++ b/src/ComputationalModels/PostProcessors.jl @@ -71,17 +71,16 @@ function (obj::PostProcessor{<:DynamicNonlinearModel,<:Any,<:Any})(Λ) obj.cache[1](obj, obj.cache[2]...) end -function Jacobian(uh) - F, _, J = Kinematics(Mechano).metrics +function Jacobian(uh,km) + F, _, J = get_Kinematics(km) J ∘ F ∘ ∇(uh) end -function Cauchy(physmodel::ThermoElectroMechano, uh, φh, θh, Ω, dΩ, Λ=1.0) +function Cauchy(physmodel::ThermoElectroMechano,kine::NTuple{3,KinematicModel}, uh, φh, θh, Ω, dΩ, Λ=1.0) DΨ = physmodel(Λ) - Kinematic_mec = Kinematics(Mechano) - Kinematic_elec = Kinematics(Electro) - F, _, _ = get_Kinematics(Kinematic_mec) - E = get_Kinematics(Kinematic_elec) + + F, _, _ = get_Kinematics(kine[1]) + E = get_Kinematics(kine[2]) ∂Ψu = DΨ[2] refL2 = ReferenceFE(lagrangian, Float64, 0) ref = ReferenceFE(lagrangian, Float64, 1) @@ -101,38 +100,36 @@ function Cauchy(physmodel::ThermoElectroMechano, uh, φh, θh, Ω, dΩ, Λ=1.0) end -function Cauchy(model::Elasto, uh, unh, state_vars, Ω, dΩ, t, Δt) - σh = Cauchy(model, uh) +function Cauchy(model::Elasto,km::KinematicModel,uh, unh, state_vars, Ω, dΩ, t, Δt) + σh = Cauchy(model,km,uh) interpolate_L2_tensor(σh, Ω, dΩ) end -function Cauchy(model::ViscoElastic, uh, unh, state_vars, Ω, dΩ, t, Δt) - σh = Cauchy(model, uh, unh, state_vars, Δt) +function Cauchy(model::ViscoElastic, km::KinematicModel, uh, unh, state_vars, Ω, dΩ, t, Δt) + σh = Cauchy(model, km, uh, unh, state_vars, Δt) interpolate_L2_tensor(σh, Ω, dΩ) end -function Cauchy(model::Elasto, uh, vars...) +function Cauchy(model::Elasto, km::KinematicModel,uh, vars...) _, ∂Ψu, _ = model() - F, _, _ = get_Kinematics(model.Kinematic) + F, _, _ = get_Kinematics(km) ∂Ψu ∘ (F∘∇(uh)) end -function Cauchy(model::ViscoElastic, uh, unh, states, Δt) +function Cauchy(model::ViscoElastic, km::KinematicModel, uh, unh, states, Δt) _, ∂Ψu, _ = model(Δt=Δt) - F, _, _ = get_Kinematics(model.Kinematic) + F, _, _ = get_Kinematics(km) ∂Ψu ∘ (F∘∇(uh), F∘∇(unh), states...) end -function Entropy(physmodel::ThermoElectroMechano, uh, φh, θh, Ω, dΩ, Λ=1.0) +function Entropy(physmodel::ThermoElectroMechano, kine::NTuple{3,KinematicModel}, uh, φh, θh, Ω, dΩ, Λ=1.0) DΨ = physmodel(Λ) - Kinematic_mec = Kinematics(Mechano) - Kinematic_elec = Kinematics(Electro) - F, _, _ = get_Kinematics(Kinematic_mec) - E = get_Kinematics(Kinematic_elec) + F,_,_ = get_Kinematics(kine[1]; Λ=Λ) + E = get_Kinematics(kine[2]; Λ=Λ) η = DΨ[11] refL2 = ReferenceFE(lagrangian, Float64, 0) ref = ReferenceFE(lagrangian, Float64, 1) @@ -142,12 +139,10 @@ function Entropy(physmodel::ThermoElectroMechano, uh, φh, θh, Ω, dΩ, Λ=1.0) return ηh end -function D0(physmodel::ThermoElectroMechano, uh, φh, θh, Ω, dΩ, Λ=1.0) +function D0(physmodel::ThermoElectroMechano, kine::NTuple{3,KinematicModel}, uh, φh, θh, Ω, dΩ, Λ=1.0) DΨ = physmodel(Λ) - Kinematic_mec = Kinematics(Mechano) - Kinematic_elec = Kinematics(Electro) - F, _, _ = get_Kinematics(Kinematic_mec) - E = get_Kinematics(Kinematic_elec) + F,_,_ = get_Kinematics(kine[1]; Λ=Λ) + E = get_Kinematics(kine[2]; Λ=Λ) ∂ΨE = DΨ[3] refL2 = ReferenceFE(lagrangian, Float64, 0) ref = ReferenceFE(lagrangian, Float64, 1) diff --git a/src/Exports.jl b/src/Exports.jl index 272e6db..7eb5f74 100644 --- a/src/Exports.jl +++ b/src/Exports.jl @@ -79,6 +79,7 @@ end @publish PhysicalModels EnergyInterpolationScheme @publish PhysicalModels update_state! @publish PhysicalModels Kinematics +@publish PhysicalModels Solid @publish PhysicalModels KinematicModel @publish PhysicalModels EvolutiveKinematics @publish PhysicalModels get_Kinematics diff --git a/src/PhysicalModels/ElectricalModels.jl b/src/PhysicalModels/ElectricalModels.jl index 33558ec..c4b5c20 100644 --- a/src/PhysicalModels/ElectricalModels.jl +++ b/src/PhysicalModels/ElectricalModels.jl @@ -3,10 +3,9 @@ # Electrical models # =================== -struct IdealDielectric{A} <: Electro +struct IdealDielectric <: Electro ε::Float64 - Kinematic::A - function IdealDielectric(; ε::Float64, Kinematic::KinematicModel=Kinematics(Electro)) - new{typeof(Kinematic)}(ε, Kinematic) + function IdealDielectric(; ε::Float64) + new(ε) end end diff --git a/src/PhysicalModels/ElectroMechanicalModels.jl b/src/PhysicalModels/ElectroMechanicalModels.jl index c51fa9d..45344e7 100644 --- a/src/PhysicalModels/ElectroMechanicalModels.jl +++ b/src/PhysicalModels/ElectroMechanicalModels.jl @@ -59,8 +59,8 @@ end function _getCoupling(elec::Electro, mec::Mechano, Λ::Float64) - _, H, J = get_Kinematics(mec.Kinematic; Λ=Λ) - + J(F) = det(F) + H(F) = det(F) * inv(F)' # Energy # HE(F, E) = H(F) * E HEHE(F, E) = HE(F, E) ⋅ HE(F, E) diff --git a/src/PhysicalModels/KinematicModels.jl b/src/PhysicalModels/KinematicModels.jl index 1979841..783e3dc 100644 --- a/src/PhysicalModels/KinematicModels.jl +++ b/src/PhysicalModels/KinematicModels.jl @@ -1,35 +1,57 @@ abstract type KinematicModel end +abstract type KinematicDescription end + +abstract type Solid <: KinematicDescription end +abstract type SolidShell <: KinematicDescription end + -struct KinematicDescription{Kind} end get_Kinematics(::KinematicModel; Λ::Float64) = @abstractmethod -struct Kinematics{T} <: KinematicModel + +struct Kinematics{T1,T2} <: KinematicModel metrics - function Kinematics(::Type{M}; F::Function=(∇u) -> one(∇u) + ∇u) where {M <: Mechano} + function Kinematics(::Type{Mechano}, ::Type{Solid}; F::Function=(∇u) -> one(∇u) + ∇u) J(F) = det(F) H(F) = det(F) * inv(F)' metrics = (F, H, J) - new{Mechano}(metrics) + new{Mechano,Solid}(metrics) end - function Kinematics(::Type{Electro}; E::Function=(∇φ) -> -∇φ) + # function Kinematics(::Type{Mechano}, ::Type{SolidShell}) + # function F(∇u, N) + # F̂ = ∇u + one(∇u) + # Ĵ = det(F̂) + # Ĥ = (det(F̂) * inv(F̂)') + # HN = (Ĥ * N) / norm(Ĥ * N) + # F̂ + (1 / Ĵ - 1.0) * (HN ⊗ N) + # end + # J(F) = det(F) + # H(F) = det(F) * inv(F)' + # metrics = (F, H, J) + # new{Mechano,SolidShell}(metrics) + # end + + function Kinematics(::Type{Electro}, ::Type{Solid}; E::Function=(∇φ) -> -∇φ) metrics = (E) - new{Electro}(metrics) + new{Electro,Solid}(metrics) end - function Kinematics(::Type{Magneto}; H::Function=(∇φ) -> -∇φ) + function Kinematics(::Type{Magneto}, ::Type{Solid}; H::Function=(∇φ) -> -∇φ) metrics = (H) - new{Magneto}(metrics) + new{Magneto,Solid}(metrics) + end + function Kinematics(::Type{Thermo}, ::Type{Solid}) + new{Thermo,Solid}(nothing) end end get_Kinematics(obj::Kinematics; Λ=1.0) = obj.metrics -function getIsoInvariants(obj::Kinematics{Mechano}) +function getIsoInvariants(obj::Kinematics{Mechano,<:KinematicDescription}) F, H, J = obj.metrics I1(F) = tr(F' * F) I2(F) = tr(H(F)' * H(F)) @@ -38,7 +60,7 @@ function getIsoInvariants(obj::Kinematics{Mechano}) end -function getIsoInvariants(obj_m::Kinematics{Mechano},obj_e::Kinematics{Electro}) +function getIsoInvariants(obj_m::Kinematics{Mechano,<:KinematicDescription}, obj_e::Kinematics{Electro,<:KinematicDescription}) F, H, J = obj_m.metrics E = obj_e.metrics I1(F) = tr(F' * F) @@ -65,7 +87,7 @@ struct EvolutiveKinematics{T} <: KinematicModel function EvolutiveKinematics(::Type{Mechano}, δ::Float64; F::Function=(t) -> ((∇u) -> one(∇u) + ∇u)) # F_(∇u) = one(∇u) + ∇u # F(t) = (∇u)->(Fmapping(t) ∘ F_)(∇u) - J(F) = 0.5*(det(F) + sqrt(det(F)^2 + δ^2)) + J(F) = 0.5 * (det(F) + sqrt(det(F)^2 + δ^2)) H(F) = J(F) * inv(F)' metrics = (F, H, J) new{Mechano}(metrics) diff --git a/src/PhysicalModels/MagneticModels.jl b/src/PhysicalModels/MagneticModels.jl index 480c53e..c0b0943 100644 --- a/src/PhysicalModels/MagneticModels.jl +++ b/src/PhysicalModels/MagneticModels.jl @@ -8,10 +8,9 @@ struct Magnetic <: Magneto μ::Float64 αr::Ref{Float64} χe::Float64 - Kinematic::Kinematics{Magneto} -function Magnetic(; μ::Float64, αr::Ref{Float64}, χe::Float64=0.0, Kinematic::Kinematics{Magneto}=Kinematics(Magneto)) - new(μ, αr, χe, Kinematic) +function Magnetic(; μ::Float64, αr::Ref{Float64}, χe::Float64=0.0) + new(μ, αr, χe) end function (obj::Magnetic)(Λ::Float64=1.0) μ, αr, χe = obj.μ, obj.αr, obj.χe @@ -27,12 +26,11 @@ end end -struct IdealMagnetic{A} <: Magneto +struct IdealMagnetic <: Magneto μ::Float64 χe::Float64 - Kinematic::A - function IdealMagnetic(; μ::Float64, χe::Float64=0.0, Kinematic::KinematicModel=Kinematics(Magneto)) - new{typeof(Kinematic)}(μ, χe, Kinematic) + function IdealMagnetic(; μ::Float64, χe::Float64=0.0) + new(μ, χe) end function (obj::IdealMagnetic)(Λ::Float64=1.0) @@ -73,12 +71,11 @@ struct IdealMagnetic{A} <: Magneto end -struct IdealMagnetic2D{A} <: Magneto +struct IdealMagnetic2D <: Magneto μ::Float64 χe::Float64 - Kinematic::A - function IdealMagnetic2D(; μ::Float64, χe::Float64=0.0, Kinematic::KinematicModel=Kinematics(Magneto)) - new{typeof(Kinematic)}(μ, χe, Kinematic) + function IdealMagnetic2D(; μ::Float64, χe::Float64=0.0) + new(μ, χe) end function (obj::IdealMagnetic2D)(Λ::Float64=1.0) @@ -115,16 +112,15 @@ struct IdealMagnetic2D{A} <: Magneto end -struct HardMagnetic{A} <: Magneto +struct HardMagnetic <: Magneto μ::Float64 αr::Float64 χe::Float64 χr::Float64 βmok::Float64 βcoup::Float64 - Kinematic::A - function HardMagnetic(; μ::Float64, αr::Float64, χe::Float64=0.0, χr::Float64=8.0, βmok::Float64=1.0, βcoup::Float64=1.0, Kinematic::KinematicModel=Kinematics(Magneto)) - new{typeof(Kinematic)}(μ, αr, χe, χr, βmok, βcoup, Kinematic) + function HardMagnetic(; μ::Float64, αr::Float64, χe::Float64=0.0, χr::Float64=8.0, βmok::Float64=1.0, βcoup::Float64=1.0) + new(μ, αr, χe, χr, βmok, βcoup) end function (obj::HardMagnetic)(Λ::Float64=1.0) @@ -164,16 +160,15 @@ struct HardMagnetic{A} <: Magneto end -struct HardMagnetic2D{A} <: Magneto +struct HardMagnetic2D <: Magneto μ::Float64 αr::Float64 χe::Float64 χr::Float64 βmok::Float64 βcoup::Float64 - Kinematic::A - function HardMagnetic2D(; μ::Float64, αr::Float64, χe::Float64=0.0, χr::Float64=8.0, βmok::Float64=1.0, βcoup::Float64=1.0, Kinematic::KinematicModel=Kinematics(Magneto)) - new{typeof(Kinematic)}(μ, αr, χe, χr, βmok, βcoup, Kinematic) + function HardMagnetic2D(; μ::Float64, αr::Float64, χe::Float64=0.0, χr::Float64=8.0, βmok::Float64=1.0, βcoup::Float64=1.0) + new(μ, αr, χe, χr, βmok, βcoup) end function (obj::HardMagnetic2D)(Λ::Float64=1.0) diff --git a/src/PhysicalModels/MagnetoMechanicalModels.jl b/src/PhysicalModels/MagnetoMechanicalModels.jl index 3eb0016..2bab4c2 100644 --- a/src/PhysicalModels/MagnetoMechanicalModels.jl +++ b/src/PhysicalModels/MagnetoMechanicalModels.jl @@ -51,9 +51,10 @@ function _getCoupling(mag::HardMagnetic, mec::Mechano, Λ::Float64) # Hard magnetics in ultra-soft magnetorheological elastomers enhance fracture toughness and # delay crack propagation, Journal of the Mechanics and Physics of Solids, - _, H, J = get_Kinematics(mec.Kinematic; Λ=Λ) - μ, αr, χe, χr, βcoup, βmok = mag.μ, mag.αr, mag.χe, mag.χr, mag.βcoup, mag.βmok + μ, αr, χe, χr, βcoup, βmok = mag.μ, mag.αr, mag.χe, mag.χr, mag.βcoup, mag.βmok + J(F) = det(F) + H(F) = det(F) * inv(F)' αr *= Λ #------------------------------------------------------------------------------------- # FIRST TERM @@ -131,8 +132,9 @@ function _getCoupling(mag::HardMagnetic2D, mec::Mechano, Λ::Float64) # Hard magnetics in ultra-soft magnetorheological elastomers enhance fracture toughness and # delay crack propagation, Journal of the Mechanics and Physics of Solids, - _, H, J = get_Kinematics(mec.Kinematic; Λ=Λ) μ, αr, χe, χr, βcoup, βmok = mag.μ, mag.αr, mag.χe, mag.χr, mag.βcoup, mag.βmok + J(F) = det(F) + H(F) = det(F) * inv(F)' αr *= Λ # #------------------------------------------------------------------------------------- diff --git a/src/PhysicalModels/MechanicalModels.jl b/src/PhysicalModels/MechanicalModels.jl index 17a95fb..4c9228e 100644 --- a/src/PhysicalModels/MechanicalModels.jl +++ b/src/PhysicalModels/MechanicalModels.jl @@ -2,17 +2,16 @@ # Regularization of Mechanical models # ============================================ -struct HessianRegularization{M<:Mechano,B} <: Mechano - mechano::M +struct HessianRegularization <: Mechano + mechano::Mechano δ::Float64 - Kinematic::B - function HessianRegularization(mechano::M; δ::Float64=1.0e-6) where {M<:Mechano} - new{M,typeof(mechano.Kinematic)}(mechano, δ, mechano.Kinematic) + function HessianRegularization(mechano::Mechano; δ::Float64=1.0e-6) + new(mechano, δ) end - function HessianRegularization(; mechano::M, δ::Float64=1.0e-6) where {M<:Mechano} - new{M,typeof(mechano.Kinematic)}(mechano, δ, mechano.Kinematic) + function HessianRegularization(; mechano::Mechano, δ::Float64=1.0e-6) + new(mechano, δ) end function (obj::HessianRegularization)(Λ::Float64=1.0) @@ -30,25 +29,25 @@ struct HessianRegularization{M<:Mechano,B} <: Mechano end -struct Hessian∇JRegularization{M<:Mechano,B} <: Mechano - mechano::M +struct Hessian∇JRegularization <: Mechano + mechano::Mechano δ::Float64 κ::Float64 - Kinematic::B - function Hessian∇JRegularization(mechano::M; δ::Float64=1.0e-6, κ::Float64=1.0) where {M<:Mechano} - new{M,typeof(mechano.Kinematic)}(mechano, δ, κ, mechano.Kinematic) + function Hessian∇JRegularization(mechano::Mechano; δ::Float64=1.0e-6, κ::Float64=1.0) + new(mechano, δ, κ) end - function Hessian∇JRegularization(; mechano::M, δ::Float64=1.0e-6, κ::Float64=1.0) where {M<:Mechano} - new{M,typeof(mechano.Kinematic)}(mechano, δ, κ, mechano.Kinematic) + function Hessian∇JRegularization(; mechano::Mechano, δ::Float64=1.0e-6, κ::Float64=1.0) + new(mechano, δ, κ) end function (obj::Hessian∇JRegularization)(Λ::Float64=1.0) Ψs, ∂Ψs, ∂2Ψs = obj.mechano() - _, H, J = get_Kinematics(obj.mechano.Kinematic; Λ=Λ) δ, κ = obj.δ, obj.κ + J(F) = det(F) + H(F) = det(F) * inv(F)' Ψ(F, Jh) = Ψs(F) + 0.5 * κ * (J(F) - Jh)^2 ∂Ψ(F, Jh) = ∂Ψs(F) + κ * (J(F) - Jh) * H(F) ∂2Ψ_(F, Jh) = ∂2Ψs(F) + κ * (H(F) ⊗ H(F)) + κ * (J(F) - Jh) * _∂H∂F_2D() @@ -97,9 +96,8 @@ end struct ComposedIsoElastic <: IsoElastic Model1::IsoElastic Model2::IsoElastic - Kinematic function ComposedIsoElastic(model1::IsoElastic, model2::IsoElastic) - new(model1, model2, model1.Kinematic) + new(model1, model2) end function (obj::ComposedIsoElastic)(Λ::Float64=1.0) Ψ1, ∂Ψu1, ∂Ψuu1 = obj.Model1(Λ) @@ -119,9 +117,8 @@ end struct ComposedAnisoElastic <: AnisoElastic Model1::IsoElastic Model2::AnisoElastic - Kinematic function ComposedAnisoElastic(model1::IsoElastic, model2::AnisoElastic) - new(model1, model2, model1.Kinematic) + new(model1, model2) end function (obj::ComposedAnisoElastic)(Λ::Float64=1.0) DΨ1 = obj.Model1(Λ) @@ -164,20 +161,22 @@ transpose(x::NTuple{N,<:Tuple{<:Function,<:Function,<:Function}}) where N = map( # Mechanical models # =================== -struct Yeoh3D{A} <: IsoElastic +struct Yeoh3D <: IsoElastic λ::Float64 C10::Float64 C20::Float64 C30::Float64 - Kinematic::A - function Yeoh3D(; λ::Float64, C10::Float64, C20::Float64, C30::Float64, Kinematic::KinematicModel=Kinematics(Mechano)) - new{typeof(Kinematic)}(λ, C10, C20, C30, Kinematic) + function Yeoh3D(; λ::Float64, C10::Float64, C20::Float64, C30::Float64) + new(λ, C10, C20, C30) end function (obj::Yeoh3D)(Λ::Float64=1.0; Threshold=0.01) - _, H, J = get_Kinematics(obj.Kinematic; Λ=Λ) + λ, C10, C20, C30 = obj.λ, obj.C10, obj.C20, obj.C30 + J(F) = det(F) + H(F) = det(F) * inv(F)' + # Free energy I1(F) = tr((F)' * F) ψ(F) = C10 * (I1(F) - 3) + C20 * (I1(F) - 3)^2 + C30 * (I1(F) - 3)^3 - 2 * C10 * log(J(F)) + 0.5 * λ * (J(F) - 1)^2 @@ -199,13 +198,12 @@ struct Yeoh3D{A} <: IsoElastic end end -struct LinearElasticity2D{A} <: IsoElastic +struct LinearElasticity2D <: IsoElastic λ::Float64 μ::Float64 ρ::Float64 - Kinematic::A - function LinearElasticity2D(; λ::Float64, μ::Float64, ρ::Float64=0.0, Kinematic::KinematicModel=Kinematics(Mechano)) - new{typeof(Kinematic)}(λ, μ, ρ, Kinematic) + function LinearElasticity2D(; λ::Float64, μ::Float64, ρ::Float64=0.0) + new(λ, μ, ρ) end function (obj::LinearElasticity2D)(Λ::Float64=1.0) @@ -219,13 +217,12 @@ struct LinearElasticity2D{A} <: IsoElastic end -mutable struct LinearElasticity3D{A} <: IsoElastic +mutable struct LinearElasticity3D <: IsoElastic λ::Float64 μ::Float64 ρ::Float64 - Kinematic::A - function LinearElasticity3D(; λ::Float64, μ::Float64, ρ::Float64=0.0, Kinematic::KinematicModel=Kinematics(Mechano)) - new{typeof(Kinematic)}(λ, μ, ρ, Kinematic) + function LinearElasticity3D(; λ::Float64, μ::Float64, ρ::Float64=0.0) + new(λ, μ, ρ) end function (obj::LinearElasticity3D)(Λ::Float64=1.0) @@ -239,16 +236,16 @@ mutable struct LinearElasticity3D{A} <: IsoElastic end -struct VolumetricEnergy{A} <: IsoElastic +struct VolumetricEnergy <: IsoElastic λ::Float64 - Kinematic::A - function VolumetricEnergy(; λ::Float64, Kinematic::KinematicModel=Kinematics(Elasto)) - new{typeof(Kinematic)}(λ, Kinematic) + function VolumetricEnergy(; λ::Float64) + new(λ) end function (obj::VolumetricEnergy)(Λ::Float64=1.0) - _, H, J = get_Kinematics(obj.Kinematic; Λ=Λ) λ = obj.λ + J(F) = det(F) + H(F) = det(F) * inv(F)' Ψ(F) = (λ / 2.0) * (J(F) - 1)^2 ∂Ψ_∂J(F) = λ * (J(F) - 1) ∂Ψ2_∂J2(F) = λ @@ -259,18 +256,18 @@ struct VolumetricEnergy{A} <: IsoElastic end -struct NeoHookean3D{A} <: IsoElastic +struct NeoHookean3D <: IsoElastic λ::Float64 μ::Float64 ρ::Float64 - Kinematic::A - function NeoHookean3D(; λ::Float64, μ::Float64, ρ::Float64=0.0, Kinematic::KinematicModel=Kinematics(Mechano)) - new{typeof(Kinematic)}(λ, μ, ρ, Kinematic) + function NeoHookean3D(; λ::Float64, μ::Float64, ρ::Float64=0.0) + new(λ, μ, ρ) end function (obj::NeoHookean3D)(Λ::Float64=1.0; Threshold=0.01) - _, H, J = get_Kinematics(obj.Kinematic; Λ=Λ) λ, μ = obj.λ, obj.μ + J(F) = det(F) + H(F) = det(F) * inv(F)' Ψ(F) = μ / 2 * tr((F)' * F) - μ * logreg(J(F)) + (λ / 2) * (J(F) - 1)^2 - 3.0 * (μ / 2.0) ∂log∂J(J) = J >= Threshold ? 1 / J : (2 / Threshold - J / (Threshold^2)) @@ -285,19 +282,20 @@ struct NeoHookean3D{A} <: IsoElastic end -struct MooneyRivlin3D{A} <: IsoElastic +struct MooneyRivlin3D <: IsoElastic λ::Float64 μ1::Float64 μ2::Float64 ρ::Float64 - Kinematic::A - function MooneyRivlin3D(; λ::Float64, μ1::Float64, μ2::Float64, ρ::Float64=0.0, Kinematic::KinematicModel=Kinematics(Mechano)) - new{typeof(Kinematic)}(λ, μ1, μ2, ρ, Kinematic) + function MooneyRivlin3D(; λ::Float64, μ1::Float64, μ2::Float64, ρ::Float64=0.0) + new(λ, μ1, μ2, ρ) end function (obj::MooneyRivlin3D)(Λ::Float64=1.0; Threshold=0.01) - _, H, J = get_Kinematics(obj.Kinematic; Λ=Λ) λ, μ1, μ2 = obj.λ, obj.μ1, obj.μ2 + J(F) = det(F) + H(F) = det(F) * inv(F)' + Ψ(F) = μ1 / 2 * tr((F)' * F) + μ2 / 2.0 * tr((H(F))' * H(F)) - (μ1 + 2 * μ2) * logreg(J(F)) + (λ / 2.0) * (J(F) - 1)^2 - (3.0 / 2.0) * (μ1 + μ2) ∂Ψ_∂F(F) = μ1 * F @@ -315,20 +313,20 @@ struct MooneyRivlin3D{A} <: IsoElastic end -struct MooneyRivlin2D{A} <: IsoElastic +struct MooneyRivlin2D <: IsoElastic λ::Float64 μ1::Float64 μ2::Float64 ρ::Float64 - Kinematic::A - function MooneyRivlin2D(; λ::Float64, μ1::Float64, μ2::Float64, ρ::Float64=0.0, Kinematic::KinematicModel=Kinematics(Mechano)) - new{typeof(Kinematic)}(λ, μ1, μ2, ρ, Kinematic) + function MooneyRivlin2D(; λ::Float64, μ1::Float64, μ2::Float64, ρ::Float64=0.0) + new(λ, μ1, μ2, ρ) end function (obj::MooneyRivlin2D)(Λ::Float64=1.0; Threshold=0.01) - _, H, J = get_Kinematics(obj.Kinematic; Λ=Λ) λ, μ1, μ2 = obj.λ, obj.μ1, obj.μ2 + J(F) = det(F) + H(F) = det(F) * inv(F)' Ψ(F) = (μ1 / 2 + μ2 / 2) * tr((F)' * F) + μ2 / 2.0 * J(F)^2 - (μ1 + 2 * μ2) * logreg(J(F)) + (λ / 2.0) * (J(F) - 1)^2 ∂Ψ_(F) = ForwardDiff.gradient(F -> Ψ(F), get_array(F)) @@ -341,21 +339,22 @@ struct MooneyRivlin2D{A} <: IsoElastic end -struct NonlinearMooneyRivlin3D{A} <: IsoElastic +struct NonlinearMooneyRivlin3D <: IsoElastic λ::Float64 μ1::Float64 μ2::Float64 α::Float64 β::Float64 ρ::Float64 - Kinematic::A - function NonlinearMooneyRivlin3D(; λ::Float64, μ1::Float64, μ2::Float64, α::Float64, β::Float64, ρ::Float64=0.0, Kinematic::KinematicModel=Kinematics(Mechano)) - new{typeof(Kinematic)}(λ, μ1, μ2, α, β, ρ, Kinematic) + function NonlinearMooneyRivlin3D(; λ::Float64, μ1::Float64, μ2::Float64, α::Float64, β::Float64, ρ::Float64=0.0) + new(λ, μ1, μ2, α, β, ρ) end function (obj::NonlinearMooneyRivlin3D)(Λ::Float64=1.0; Threshold=0.01) - _, H, J = get_Kinematics(obj.Kinematic; Λ=Λ) λ, μ1, μ2, α, β = obj.λ, obj.μ1, obj.μ2, obj.α, obj.β + J(F) = det(F) + H(F) = det(F) * inv(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) * logreg(J(F)) + (λ / 2.0) * (J(F) - 1)^2 @@ -375,22 +374,21 @@ struct NonlinearMooneyRivlin3D{A} <: IsoElastic end -struct NonlinearMooneyRivlin2D{A} <: IsoElastic +struct NonlinearMooneyRivlin2D <: IsoElastic λ::Float64 μ1::Float64 μ2::Float64 α::Float64 β::Float64 ρ::Float64 - Kinematic::A - function NonlinearMooneyRivlin2D(; λ::Float64, μ1::Float64, μ2::Float64, α::Float64, β::Float64, ρ::Float64=0.0, Kinematic::KinematicModel=Kinematics(Mechano)) - new{typeof(Kinematic)}(λ, μ1, μ2, α, β, ρ, Kinematic) + function NonlinearMooneyRivlin2D(; λ::Float64, μ1::Float64, μ2::Float64, α::Float64, β::Float64, ρ::Float64=0.0) + new(λ, μ1, μ2, α, β, ρ) end function (obj::NonlinearMooneyRivlin2D)(Λ::Float64=1.0; Threshold=0.01) - _, H, J = get_Kinematics(obj.Kinematic; Λ=Λ) λ, μ1, μ2, α, β = obj.λ, obj.μ1, obj.μ2, obj.α, obj.β - + J(F) = det(F) + H(F) = det(F) * inv(F)' Ψ(F) = μ1 / (2.0 * α * 3.0^(α - 1)) * (tr((F)' * F) + 1.0)^α + μ2 / (2.0 * β * 3.0^(β - 1)) * (tr((F)' * F) + J(F)^2)^β - (μ1 + 2.0 * μ2) * logreg(J(F)) + (λ / 2.0) * (J(F) - 1)^2 @@ -412,7 +410,7 @@ struct NonlinearMooneyRivlin2D{A} <: IsoElastic end -struct NonlinearMooneyRivlin2D_CV{A} <: IsoElastic +struct NonlinearMooneyRivlin2D_CV <: IsoElastic λ::Float64 μ1::Float64 μ2::Float64 @@ -420,15 +418,14 @@ struct NonlinearMooneyRivlin2D_CV{A} <: IsoElastic β::Float64 γ::Float64 ρ::Float64 - Kinematic::A - function NonlinearMooneyRivlin2D_CV(; λ::Float64, μ1::Float64, μ2::Float64, α::Float64, β::Float64, γ::Float64, ρ::Float64=0.0, Kinematic::KinematicModel=Kinematics(Mechano)) - new{typeof(Kinematic)}(λ, μ1, μ2, α, β, γ, ρ, Kinematic) + function NonlinearMooneyRivlin2D_CV(; λ::Float64, μ1::Float64, μ2::Float64, α::Float64, β::Float64, γ::Float64, ρ::Float64=0.0) + new(λ, μ1, μ2, α, β, γ, ρ) end function (obj::NonlinearMooneyRivlin2D_CV)(Λ::Float64=1.0) - _, H, J = get_Kinematics(obj.Kinematic; Λ=Λ) λ, μ1, μ2, α, β, γ = obj.λ, obj.μ1, obj.μ2, obj.α, obj.β, obj.γ - + J(F) = det(F) + H(F) = det(F) * inv(F)' Ψ(F) = μ1 / (2.0 * α * 3.0^(α - 1)) * (tr((F)' * F) + 1.0)^α + μ2 / (2.0 * β * 3.0^(β - 1)) * (tr((F)' * F) + J(F)^2)^β - (μ1 + 2.0 * μ2) * log(J(F)) + (λ) * (J(F)^(γ) + J(F)^(-γ)) @@ -448,7 +445,7 @@ struct NonlinearMooneyRivlin2D_CV{A} <: IsoElastic end -struct NonlinearMooneyRivlin_CV{A} <: IsoElastic +struct NonlinearMooneyRivlin_CV <: IsoElastic λ::Float64 μ1::Float64 μ2::Float64 @@ -456,15 +453,14 @@ struct NonlinearMooneyRivlin_CV{A} <: IsoElastic β::Float64 γ::Float64 ρ::Float64 - Kinematic::A - function NonlinearMooneyRivlin_CV(; λ::Float64, μ1::Float64, μ2::Float64, α::Float64, β::Float64, γ::Float64, ρ::Float64=0.0, Kinematic::KinematicModel=Kinematics(Mechano)) - new{typeof(Kinematic)}(λ, μ1, μ2, α, β, γ, ρ, Kinematic) + function NonlinearMooneyRivlin_CV(; λ::Float64, μ1::Float64, μ2::Float64, α::Float64, β::Float64, γ::Float64, ρ::Float64=0.0) + new(λ, μ1, μ2, α, β, γ, ρ) end function (obj::NonlinearMooneyRivlin_CV)(Λ::Float64=1.0) - _, H, J = get_Kinematics(obj.Kinematic; Λ=Λ) λ, μ1, μ2, α, β, γ = obj.λ, obj.μ1, obj.μ2, obj.α, obj.β, obj.γ - + J(F) = det(F) + H(F) = det(F) * inv(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)^(-γ)) @@ -486,21 +482,20 @@ struct NonlinearMooneyRivlin_CV{A} <: IsoElastic end -struct NonlinearNeoHookean_CV{A} <: IsoElastic +struct NonlinearNeoHookean_CV <: IsoElastic λ::Float64 μ::Float64 α::Float64 γ::Float64 ρ::Float64 - Kinematic::A - function NonlinearNeoHookean_CV(; λ::Float64, μ::Float64, α::Float64, γ::Float64, ρ::Float64=0.0, Kinematic::KinematicModel=Kinematics(Mechano)) - new{typeof(Kinematic)}(λ, μ, α, γ, ρ, Kinematic) + function NonlinearNeoHookean_CV(; λ::Float64, μ::Float64, α::Float64, γ::Float64, ρ::Float64=0.0) + new(λ, μ, α, γ, ρ) end function (obj::NonlinearNeoHookean_CV)(Λ::Float64=1.0) - _, H, J = get_Kinematics(obj.Kinematic; Λ=Λ) λ, μ, α, γ = obj.λ, obj.μ, obj.α, obj.γ - + J(F) = det(F) + H(F) = det(F) * inv(F)' Ψ(F) = μ / (2.0 * α * 3.0^(α - 1)) * (tr((F)' * F))^α - μ * log(J(F)) + λ * (J(F)^(γ) + J(F)^(-γ)) ∂Ψ_∂F(F) = ((μ / (3.0^(α - 1)) * (tr((F)' * F))^(α - 1))) * F @@ -518,21 +513,20 @@ struct NonlinearNeoHookean_CV{A} <: IsoElastic end -struct NonlinearIncompressibleMooneyRivlin2D_CV{A} <: IsoElastic +struct NonlinearIncompressibleMooneyRivlin2D_CV <: IsoElastic λ::Float64 μ::Float64 α::Float64 γ::Float64 ρ::Float64 - Kinematic::A - function NonlinearIncompressibleMooneyRivlin2D_CV(; λ::Float64, μ::Float64, α::Float64, γ::Float64, ρ::Float64=0.0, Kinematic::KinematicModel=Kinematics(Mechano)) - new{typeof(Kinematic)}(λ, μ, α, γ, ρ, Kinematic) + function NonlinearIncompressibleMooneyRivlin2D_CV(; λ::Float64, μ::Float64, α::Float64, γ::Float64, ρ::Float64=0.0) + new(λ, μ, α, γ, ρ) end function (obj::NonlinearIncompressibleMooneyRivlin2D_CV)(Λ::Float64=1.0) - _, H, J = get_Kinematics(obj.Kinematic; Λ=Λ) λ, μ, α, γ = obj.λ, obj.μ, obj.α, obj.γ - + J(F) = det(F) + H(F) = det(F) * inv(F)' e(F) = (tr((F)' * F) + 1.0) * J(F)^(-2 / 3) ∂e_∂F(F) = 2 * J(F)^(-2 / 3) * F ∂e_∂J(F) = -(2 / 3) * (tr((F)' * F) + 1.0) * J(F)^(-5 / 3) @@ -566,18 +560,17 @@ struct NonlinearIncompressibleMooneyRivlin2D_CV{A} <: IsoElastic end -struct EightChain{A} <: IsoElastic +struct EightChain <: IsoElastic μ::Float64 N::Float64 - Kinematic::A - function EightChain(; μ::Float64, N::Float64, Kinematic::KinematicModel=Kinematics(Elasto)) - new{typeof(Kinematic)}(μ, N, Kinematic) + function EightChain(; μ::Float64, N::Float64) + new(μ, N) end function (obj::EightChain)(Λ::Float64=1.0) - _, H, J = get_Kinematics(obj.Kinematic; Λ=Λ) μ, N = obj.μ, obj.N - + J(F) = det(F) + H(F) = det(F) * inv(F)' Ψ(F) = begin C = F' * F C_iso = J(F)^(-2 / 3) * C @@ -629,19 +622,19 @@ struct EightChain{A} <: IsoElastic end -struct TransverseIsotropy3D{A} <: AnisoElastic +struct TransverseIsotropy3D <: AnisoElastic μ::Float64 α::Float64 β::Float64 ρ::Float64 - Kinematic::A - function TransverseIsotropy3D(; μ::Float64, α::Float64, β::Float64, ρ::Float64=0.0, Kinematic::KinematicModel=Kinematics(Mechano)) - new{typeof(Kinematic)}(μ, α, β, ρ, Kinematic) + function TransverseIsotropy3D(; μ::Float64, α::Float64, β::Float64, ρ::Float64=0.0) + new(μ, α, β, ρ) end function (obj::TransverseIsotropy3D)(Λ::Float64=1.0; Threshold=0.01) - _, H, J = get_Kinematics(obj.Kinematic; Λ=Λ) + J(F) = det(F) + H(F) = det(F) * inv(F)' I4(F, N) = (F * N) ⋅ (F * N) I5(F, N) = (H(F) * N) ⋅ (H(F) * N) μ, α, β = obj.μ, obj.α, obj.β @@ -664,18 +657,18 @@ struct TransverseIsotropy3D{A} <: AnisoElastic end -struct TransverseIsotropy2D{A} <: AnisoElastic +struct TransverseIsotropy2D <: AnisoElastic μ::Float64 α::Float64 β::Float64 ρ::Float64 - Kinematic::A - function TransverseIsotropy2D(; μ::Float64, α::Float64, β::Float64, ρ::Float64=0.0, Kinematic::KinematicModel=Kinematics(Mechano)) - new{typeof(Kinematic)}(μ, α, β, ρ, Kinematic) + function TransverseIsotropy2D(; μ::Float64, α::Float64, β::Float64, ρ::Float64=0.0) + new(μ, α, β, ρ) end function (obj::TransverseIsotropy2D)(Λ::Float64=1.0; Threshold=0.01) - _, H, J = get_Kinematics(obj.Kinematic; Λ=Λ) + J(F) = det(F) + H(F) = det(F) * inv(F)' I4(F, N) = (F * N) ⋅ (F * N) I5(F, N) = (H(F) * N) ⋅ (H(F) * N) μ, α, β = obj.μ, obj.α, obj.β @@ -709,10 +702,9 @@ 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)) + function HGO_4Fibers(; c1::Vector{Float64}, c2::Vector{Float64}) @assert length(c1) == length(c2) == 4 - new(c1, c2, Kinematic) + new(c1, c2) end function (obj::HGO_4Fibers)(Λ::Float64=1.0; Threshold=0.01) @@ -760,9 +752,8 @@ end struct HGO_1Fiber <: AnisoElastic c1::Float64 c2::Float64 - Kinematic::Kinematics{Mechano} - function HGO_1Fiber(; c1::Float64, c2::Float64, Kinematic::KinematicModel=Kinematics(Mechano)) - new(c1, c2, Kinematic) + function HGO_1Fiber(; c1::Float64, c2::Float64) + new(c1, c2) end function (obj::HGO_1Fiber)(Λ::Float64=1.0; Threshold=0.01) @@ -775,12 +766,12 @@ struct HGO_1Fiber <: AnisoElastic ∂Ψ∂F(F, N1) = begin M1 = N1 / norm(N1) - c1 * exp(c2 * ((F * M1) ⋅ (F * M1) - 1.0)^2.0) * ((F * M1) ⋅ (F * M1) - 1.0) * ((F * M1) ⊗ M1) + c1 * exp(c2 * ((F * M1) ⋅ (F * M1) - 1.0)^2.0) * ((F * M1) ⋅ (F * M1) - 1.0) * ((F * M1) ⊗ M1) end ∂Ψ2∂F∂F(F, N1) = begin M1 = N1 / norm(N1) - c1 * exp(c2 * ((F * M1) ⋅ (F * M1) - 1.0)^2.0) * ((4 * c2 * (((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 * exp(c2 * ((F * M1) ⋅ (F * M1) - 1.0)^2.0) * ((4 * c2 * (((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))) end return (Ψ, ∂Ψ∂F, ∂Ψ2∂F∂F) @@ -788,19 +779,19 @@ struct HGO_1Fiber <: AnisoElastic end -struct IncompressibleNeoHookean3D{A} <: IsoElastic +struct IncompressibleNeoHookean3D <: IsoElastic λ::Float64 μ::Float64 ρ::Float64 δ::Float64 - Kinematic::A - function IncompressibleNeoHookean3D(; λ::Float64, μ::Float64, ρ::Float64=0.0, δ::Float64=0.1, Kinematic::KinematicModel=Kinematics(Mechano)) - new{typeof(Kinematic)}(λ, μ, ρ, δ, Kinematic) + function IncompressibleNeoHookean3D(; λ::Float64, μ::Float64, ρ::Float64=0.0, δ::Float64=0.1) + new(λ, μ, ρ, δ) end function (obj::IncompressibleNeoHookean3D)(Λ::Float64=1.0) - _, H, J_ = get_Kinematics(obj.Kinematic; Λ=Λ) λ, μ, δ = obj.λ, obj.μ, obj.δ + J_(F) = det(F) + H(F) = det(F) * inv(F)' J(F) = 0.5 * (J_(F) + sqrt(J_(F)^2 + δ^2)) ∂J(F) = 0.5 * (1.0 + J_(F) / sqrt(J_(F)^2 + δ^2)) ∂2J(F) = 0.5 * δ^2 / ((J_(F)^2 + δ^2)^(3 / 2)) @@ -829,7 +820,9 @@ struct IncompressibleNeoHookean3D{A} <: IsoElastic return (Ψ, ∂Ψu, ∂Ψuu) end - function (obj::IncompressibleNeoHookean3D)(::KinematicDescription{:SecondPiola}, Λ::Float64=1.0) +end + + function SecondPiola(obj::IncompressibleNeoHookean3D,Λ::Float64=1.0) Ψ(C) = obj.μ / 2 * tr(C) * det(C)^(-1 / 3) S(C) = begin detC = det(C) @@ -845,23 +838,20 @@ struct IncompressibleNeoHookean3D{A} <: IsoElastic end return (Ψ, S, ∂S∂C) end -end - -struct IncompressibleNeoHookean2D{A} <: IsoElastic +struct IncompressibleNeoHookean2D <: IsoElastic λ::Float64 μ::Float64 ρ::Float64 δ::Float64 - Kinematic::A - function IncompressibleNeoHookean2D(; λ::Float64, μ::Float64, ρ::Float64=0.0, δ::Float64=0.1, Kinematic::KinematicModel=Kinematics(Mechano)) - new{typeof(Kinematic)}(λ, μ, ρ, δ, Kinematic) + function IncompressibleNeoHookean2D(; λ::Float64, μ::Float64, ρ::Float64=0.0, δ::Float64=0.1) + new(λ, μ, ρ, δ) end function (obj::IncompressibleNeoHookean2D)(Λ::Float64=1.0) - _, H, J_ = get_Kinematics(obj.Kinematic; Λ=Λ) λ, μ, δ = obj.λ, obj.μ, obj.δ - + J_(F) = det(F) + H(F) = det(F) * inv(F)' J(F) = 0.5 * (J_(F) + sqrt(J_(F)^2 + δ^2)) ∂J(F) = 0.5 * (1.0 + J_(F) / sqrt(J_(F)^2 + δ^2)) ∂2J(F) = 0.5 * δ^2 / ((J_(F)^2 + δ^2)^(3 / 2)) @@ -890,20 +880,19 @@ struct IncompressibleNeoHookean2D{A} <: IsoElastic end end -struct IncompressibleNeoHookean2D_CV{A} <: IsoElastic +struct IncompressibleNeoHookean2D_CV <: IsoElastic λ::Float64 μ::Float64 γ::Float64 ρ::Float64 - Kinematic::A - function IncompressibleNeoHookean2D_CV(; λ::Float64, μ::Float64, γ::Float64, ρ::Float64=0.0, Kinematic::KinematicModel=Kinematics(Mechano)) - new{typeof(Kinematic)}(λ, μ, γ, ρ, Kinematic) + function IncompressibleNeoHookean2D_CV(; λ::Float64, μ::Float64, γ::Float64, ρ::Float64=0.0) + new(λ, μ, γ, ρ) end function (obj::IncompressibleNeoHookean2D_CV)(Λ::Float64=1.0) - _, H, J = get_Kinematics(obj.Kinematic; Λ=Λ) λ, μ, γ = obj.λ, obj.μ, obj.γ - + J(F) = det(F) + H(F) = det(F) * inv(F)' Ψ1(F) = μ / 2 * (tr((F)' * F) + 1.0) * J(F)^(-2 / 3) Ψ2(F) = λ * (J(F)^(γ) + J(F)^(-γ)) Ψ(F) = Ψ1(F) + Ψ2(F) @@ -923,19 +912,19 @@ struct IncompressibleNeoHookean2D_CV{A} <: IsoElastic end -struct ARAP2D_regularized{A} <: IsoElastic +struct ARAP2D_regularized <: IsoElastic μ::Float64 ρ::Float64 δ::Float64 - Kinematic::A - function ARAP2D_regularized(; μ::Float64, ρ::Float64=0.0, δ::Float64=0.1, Kinematic::KinematicModel=Kinematics(Mechano)) - new{typeof(Kinematic)}(μ, ρ, δ, Kinematic) + function ARAP2D_regularized(; μ::Float64, ρ::Float64=0.0, δ::Float64=0.1) + new(μ, ρ, δ) end function (obj::ARAP2D_regularized)(Λ::Float64=1.0) - _, H, J_ = get_Kinematics(obj.Kinematic; Λ=Λ) - μ, δ = obj.μ, obj.δ + μ, δ = obj.μ, obj.δ + J_(F) = det(F) + H(F) = det(F) * inv(F)' J(F) = 0.5 * (J_(F) + sqrt(J_(F)^2 + δ^2)) ∂J(F) = 0.5 * (1.0 + J_(F) / sqrt(J_(F)^2 + δ^2)) ∂2J(F) = 0.5 * δ^2 / ((J_(F)^2 + δ^2)^(3 / 2)) @@ -964,18 +953,17 @@ struct ARAP2D_regularized{A} <: IsoElastic end -struct ARAP2D{A} <: IsoElastic +struct ARAP2D <: IsoElastic μ::Float64 ρ::Float64 - Kinematic::A - function ARAP2D(; μ::Float64, ρ::Float64=0.0, Kinematic::KinematicModel=Kinematics(Mechano)) - new{typeof(Kinematic)}(μ, ρ, Kinematic) + function ARAP2D(; μ::Float64, ρ::Float64=0.0) + new(μ, ρ) end function (obj::ARAP2D)(Λ::Float64=1.0) - _, H, J = get_Kinematics(obj.Kinematic; Λ=Λ) μ = obj.μ - + J(F) = det(F) + H(F) = det(F) * inv(F)' Ψ(F) = μ * 0.5 * J(F)^(-1) * (tr((F)' * F)) ∂Ψ_∂F(F) = μ * F * J(F)^(-1) ∂Ψ_∂J(F) = -μ / 2 * (tr((F)' * F)) * J(F)^(-2) @@ -992,21 +980,19 @@ struct ARAP2D{A} <: IsoElastic end -struct IncompressibleNeoHookean3D_2dP{A} <: Mechano +struct IncompressibleNeoHookean3D_2dP <: Mechano μ::Float64 τ::Float64 Δt::Float64 ρ::Float64 - Kinematic::A - function IncompressibleNeoHookean3D_2dP(; μ::Float64, τ::Float64, Δt::Float64, ρ::Float64=0.0, Kinematic::KinematicModel=Kinematics(Mechano)) - new{typeof(Kinematic)}(μ, τ, Δt, ρ, Kinematic) + function IncompressibleNeoHookean3D_2dP(; μ::Float64, τ::Float64, Δt::Float64, ρ::Float64=0.0) + new(μ, τ, Δt, ρ) end function (obj::IncompressibleNeoHookean3D_2dP)(Λ::Float64=1.0; Threshold=0.01) - _, H, J = get_Kinematics(obj.Kinematic; Λ=Λ) μ = obj.μ - + H(F) = det(F) * inv(F)' Ψ(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) diff --git a/src/PhysicalModels/PINNs.jl b/src/PhysicalModels/PINNs.jl index 8e516b4..2564d5f 100644 --- a/src/PhysicalModels/PINNs.jl +++ b/src/PhysicalModels/PINNs.jl @@ -48,8 +48,8 @@ struct ThermoElectroMech_PINNs{A,B,C,D} <: ThermoElectroMechano function (obj::ThermoElectroMech_PINNs)(Λ::Float64=1.0) function Ψ(F, E, δθ) - Kinematic_mec = Kinematics(Mechano) - Kinematic_elec = Kinematics(Electro) + Kinematic_mec = Kinematics(Mechano,Solid) + Kinematic_elec = Kinematics(Electro,Solid) I1, I2, I3, I4, I5 = getIsoInvariants(Kinematic_mec, Kinematic_elec) I_ = InvariantScaling([I1(F), I2(F), I3(F), I4(F, E), I5(E), δθ], obj.ϵ, obj.β) h = LayerComputation(obj.W[1], obj.b[1], [I_[1], I_[2], I_[3], I_[4], I_[5], I_[6]], SoftPlus) diff --git a/src/PhysicalModels/PhysicalModels.jl b/src/PhysicalModels/PhysicalModels.jl index 06be110..4550a45 100644 --- a/src/PhysicalModels/PhysicalModels.jl +++ b/src/PhysicalModels/PhysicalModels.jl @@ -73,8 +73,8 @@ export ThermoMechano export ThermoElectro export FlexoElectro export EnergyInterpolationScheme +export SecondPiola -export KinematicDescription export DerivativeStrategy export initializeStateVariables @@ -82,6 +82,8 @@ export updateStateVariables! export update_state! export Kinematics +export KinematicDescription +export Solid export KinematicModel export EvolutiveKinematics export get_Kinematics diff --git a/src/PhysicalModels/ThermoElectroMechanicalModels.jl b/src/PhysicalModels/ThermoElectroMechanicalModels.jl index 848ca58..689197d 100644 --- a/src/PhysicalModels/ThermoElectroMechanicalModels.jl +++ b/src/PhysicalModels/ThermoElectroMechanicalModels.jl @@ -65,7 +65,8 @@ struct ThermoElectroMech_Govindjee{T<:Thermo,E<:Electro,M<:Mechano} <: ThermoEle g(δθ) = obj.gθ(δθ) dg(δθ) = obj.dgdθ(δθ) - _, _, J = get_Kinematics(obj.mechano.Kinematic; Λ=Λ) + J(F) = det(F) + H(F) = det(F) * inv(F)' Ψer(F) = obj.thermo.α * (J(F) - 1.0) * obj.thermo.θr ΨL1(δθ) = obj.thermo.Cv * obj.thermo.θr * (1 - obj.β) * ((δθ + obj.thermo.θr) / obj.thermo.θr * (1.0 - log((δθ + obj.thermo.θr) / obj.thermo.θr)) - 1.0) ΨL3(δθ) = g(δθ) - g(0.0) - dg(0.0) * δθ @@ -125,7 +126,8 @@ struct ThermoElectroMech_Bonet{T<:Thermo,E<:Electro,M<:Mechano} <: ThermoElectro ∂gv(δθ) = (δθ+θr)^γv / θr^(γv+1) ∂∂gv(δθ) = γv*(δθ+θr)^(γv-1) / θr^(γv+1) - _, H, J = get_Kinematics(obj.mechano.Kinematic; Λ=Λ) + J(F) = det(F) + H(F) = det(F) * inv(F)' η(F)=α*(J(F) - 1.0)+Cv/γv ∂η∂J(F)=α diff --git a/src/PhysicalModels/ThermoMechanicalModels.jl b/src/PhysicalModels/ThermoMechanicalModels.jl index 603ce3d..b5f6a4d 100644 --- a/src/PhysicalModels/ThermoMechanicalModels.jl +++ b/src/PhysicalModels/ThermoMechanicalModels.jl @@ -62,8 +62,8 @@ struct ThermoMech_EntropicPolyconvex{T<:Thermo,M<:Mechano} <: ThermoMechano ϕ = obj.ϕ s = obj.s - _, H, J = get_Kinematics(obj.mechano.Kinematic; Λ=Λ) - + J(F) = det(F) + H(F) = det(F) * inv(F)' I1(F) = tr(F' * F) I2(F) = tr(H(F)' * H(F)) I3(F) = J(F) @@ -94,8 +94,8 @@ end function _getCoupling(term::Thermo, mec::Mechano, Λ::Float64) - _, H, J = get_Kinematics(mec.Kinematic; Λ=Λ) - + J(F) = det(F) + H(F) = det(F) * inv(F)' ∂Ψtm_∂J(F, δθ) = -6.0 * term.α * J(F) * δθ ∂Ψtm_u(F, δθ) = ∂Ψtm_∂J(F, δθ) * H(F) ∂Ψtm_θ(F, δθ) = -3.0 * term.α * (J(F)^2.0 - 1.0) diff --git a/src/PhysicalModels/ViscousModels.jl b/src/PhysicalModels/ViscousModels.jl index 133f7c2..ec1ed69 100644 --- a/src/PhysicalModels/ViscousModels.jl +++ b/src/PhysicalModels/ViscousModels.jl @@ -6,15 +6,14 @@ using ..TensorAlgebra # Visco elastic models # ===================== -struct ViscousIncompressible{T} <: Visco +struct ViscousIncompressible <: Visco elasto::Elasto τ::Float64 - Kinematic::T - function ViscousIncompressible(elasto; τ::Float64, kinematic::KinematicModel=Kinematics(Visco)) - new{typeof(kinematic)}(elasto, τ, kinematic) + function ViscousIncompressible(elasto; τ::Float64) + new(elasto, τ) end function (obj::ViscousIncompressible)(Λ::Float64=1.0; Δt::Float64) - Ψe, Se, ∂Se∂Ce = obj.elasto(KinematicDescription{:SecondPiola}()) + Ψe, Se, ∂Se∂Ce = SecondPiola(obj.elasto) Ψ(F, Fn, state) = Energy(obj, Δt, Ψe, Se, ∂Se∂Ce, F, Fn, state) ∂Ψ∂F(F, Fn, state) = Piola(obj, Δt, Se, ∂Se∂Ce, F, Fn, state) ∂Ψ∂F∂F(F, Fn, state) = Tangent(obj, Δt, Se, ∂Se∂Ce, F, Fn, state) @@ -27,19 +26,17 @@ function initializeStateVariables(::ViscousIncompressible, points::Measure) CellState(v, points) end -function updateStateVariables!(stateVar, obj::ViscousIncompressible, Δt, u, un) - F, _, _ = get_Kinematics(obj.Kinematic) - _, Se, ∂Se∂Ce = obj.elasto(KinematicDescription{:SecondPiola}()) +function updateStateVariables!(stateVar, obj::ViscousIncompressible, Δt, F, Fn) + _, Se, ∂Se∂Ce = SecondPiola(obj.elasto) return_mapping(A, F, Fn) = ReturnMapping(obj, Δt, Se, ∂Se∂Ce, F, Fn, A) - update_state!(return_mapping, stateVar, F∘∇(u)', F∘∇(un)') + update_state!(return_mapping, stateVar, F, Fn) end -struct GeneralizedMaxwell{T} <: ViscoElastic +struct GeneralizedMaxwell <: ViscoElastic longterm::Elasto branches::NTuple{N,Visco} where N - Kinematic::T - function GeneralizedMaxwell(longTerm::Elasto, branches::Visco...; kinematic::KinematicModel=Kinematics(Elasto)) - new{typeof(kinematic)}(longTerm,branches,kinematic) + function GeneralizedMaxwell(longTerm::Elasto, branches::Visco...) + new(longTerm,branches) end function (obj::GeneralizedMaxwell)(Λ::Float64=1.0; Δt::Float64) Ψe, ∂Ψeu, ∂Ψeuu = obj.longterm(Λ) @@ -56,10 +53,10 @@ function initializeStateVariables(model::GeneralizedMaxwell, points::Measure) map(b -> initializeStateVariables(b, points), model.branches) end -function updateStateVariables!(stateVars, model::GeneralizedMaxwell, Δt, u, un) +function updateStateVariables!(stateVars, model::GeneralizedMaxwell, Δt, F, Fn) @assert length(model.branches) == length(stateVars) for (branch, state) in zip(model.branches, stateVars) - updateStateVariables!(state, branch, Δt, u, un) + updateStateVariables!(state, branch, Δt, F, Fn) end end diff --git a/src/WeakForms/ElectroMechanics.jl b/src/WeakForms/ElectroMechanics.jl index e4a2559..c4fbe86 100644 --- a/src/WeakForms/ElectroMechanics.jl +++ b/src/WeakForms/ElectroMechanics.jl @@ -7,18 +7,18 @@ # Stagered residual # ----------------- -function residual(physicalmodel::ElectroMechano, ::Type{Mechano}, (u, φ), v, dΩ, Λ=1.0) +function residual(physicalmodel::ElectroMechano, ::Type{Mechano}, kine::NTuple{2,KinematicModel}, (u, φ), v, dΩ, Λ=1.0) DΨ = physicalmodel(Λ) - F,_,_ = get_Kinematics(physicalmodel.mechano.Kinematic; Λ=Λ) - E = get_Kinematics(physicalmodel.electro.Kinematic; Λ=Λ) + F,_,_ = get_Kinematics(kine[1]; Λ=Λ) + E = get_Kinematics(kine[2]; Λ=Λ) ∂Ψu = DΨ[2] ∫(∇(v)' ⊙ (∂Ψu ∘ (F∘∇(u)', E∘∇(φ))))dΩ end -function residual(physicalmodel::ElectroMechano, ::Type{Electro}, (u, φ), vφ, dΩ, Λ=1.0) +function residual(physicalmodel::ElectroMechano, ::Type{Electro}, kine::NTuple{2,KinematicModel}, (u, φ), vφ, dΩ, Λ=1.0) DΨ = physicalmodel(Λ) - F,_,_ = get_Kinematics(physicalmodel.mechano.Kinematic; Λ=Λ) - E = get_Kinematics(physicalmodel.electro.Kinematic; Λ=Λ) + F,_,_ = get_Kinematics(kine[1]; Λ=Λ) + E = get_Kinematics(kine[2]; Λ=Λ) ∂Ψφ = DΨ[3] -1.0*∫(∇(vφ) ⋅ (∂Ψφ ∘ (F∘∇(u)', E∘∇(φ))))dΩ end @@ -27,26 +27,26 @@ end # Stagered jacobian # ----------------- -function jacobian(physicalmodel::ElectroMechano, ::Type{Mechano}, (u, φ), du, v, dΩ, Λ=1.0) +function jacobian(physicalmodel::ElectroMechano, ::Type{Mechano}, kine::NTuple{2,KinematicModel}, (u, φ), du, v, dΩ, Λ=1.0) DΨ = physicalmodel(Λ) - F,_,_ = get_Kinematics(physicalmodel.mechano.Kinematic; Λ=Λ) - E = get_Kinematics(physicalmodel.electro.Kinematic; Λ=Λ) + F,_,_ = get_Kinematics(kine[1]; Λ=Λ) + E = get_Kinematics(kine[2]; Λ=Λ) ∂Ψuu = DΨ[4] ∫(∇(v)' ⊙ ((∂Ψuu ∘ (F∘∇(u)', E∘∇(φ))) ⊙ ∇(du)'))dΩ end -function jacobian(physicalmodel::ElectroMechano, ::Type{Electro}, (u, φ), dφ, vφ, dΩ, Λ=1.0) +function jacobian(physicalmodel::ElectroMechano, ::Type{Electro}, kine::NTuple{2,KinematicModel}, (u, φ), dφ, vφ, dΩ, Λ=1.0) DΨ = physicalmodel(Λ) - F,_,_ = get_Kinematics(physicalmodel.mechano.Kinematic; Λ=Λ) - E = get_Kinematics(physicalmodel.electro.Kinematic; Λ=Λ) + F,_,_ = get_Kinematics(kine[1]; Λ=Λ) + E = get_Kinematics(kine[2]; Λ=Λ) ∂Ψφφ = DΨ[6] ∫(∇(vφ)' ⋅ ((∂Ψφφ ∘ (F∘∇(u)', E∘∇(φ))) ⋅ ∇(dφ)))dΩ end -function jacobian(physicalmodel::ElectroMechano, ::Type{ElectroMechano}, (u, φ), (du, dφ), (v, vφ), dΩ, Λ=1.0) +function jacobian(physicalmodel::ElectroMechano, ::Type{ElectroMechano}, kine::NTuple{2,KinematicModel},(u, φ), (du, dφ), (v, vφ), dΩ, Λ=1.0) DΨ = physicalmodel(Λ) - F,_,_ = get_Kinematics(physicalmodel.mechano.Kinematic; Λ=Λ) - E = get_Kinematics(physicalmodel.electro.Kinematic; Λ=Λ) + F,_,_ = get_Kinematics(kine[1]; Λ=Λ) + E = get_Kinematics(kine[2]; Λ=Λ) ∂Ψφu = DΨ[5] -1.0*∫(∇(dφ) ⋅ ((∂Ψφu ∘ (F∘∇(u)', E∘∇(φ))) ⊙ ∇(v)'))dΩ - ∫(∇(vφ) ⋅ ((∂Ψφu ∘ (F∘∇(u)', E∘∇(φ))) ⊙ ∇(du)'))dΩ @@ -56,16 +56,16 @@ end # Monolithic strategy # ------------------- -function residual(physicalmodel::ElectroMechano, (u, φ), (v, vφ), dΩ, Λ=1.0) - residual(physicalmodel, Mechano, (u, φ), v, dΩ, Λ) + - residual(physicalmodel, Electro, (u, φ), vφ, dΩ, Λ) +function residual(physicalmodel::ElectroMechano, kine::NTuple{2,KinematicModel}, (u, φ), (v, vφ), dΩ, Λ=1.0) + residual(physicalmodel, Mechano, kine, (u, φ), v, dΩ, Λ) + + residual(physicalmodel, Electro, kine, (u, φ), vφ, dΩ, Λ) end -function jacobian(physicalmodel::ElectroMechano, (u, φ), (du, dφ), (v, vφ), dΩ, Λ=1.0) - jacobian(physicalmodel, Mechano, (u, φ), du, v, dΩ, Λ)+ - jacobian(physicalmodel, Electro, (u, φ), dφ, vφ, dΩ, Λ)+ - jacobian(physicalmodel, ElectroMechano, (u, φ), (du, dφ), (v, vφ), dΩ, Λ) +function jacobian(physicalmodel::ElectroMechano, kine::NTuple{2,KinematicModel}, (u, φ), (du, dφ), (v, vφ), dΩ, Λ=1.0) + jacobian(physicalmodel, Mechano, kine, (u, φ), du, v, dΩ, Λ)+ + jacobian(physicalmodel, Electro, kine, (u, φ), dφ, vφ, dΩ, Λ)+ + jacobian(physicalmodel, ElectroMechano, kine, (u, φ), (du, dφ), (v, vφ), dΩ, Λ) end @@ -77,18 +77,18 @@ end # Stagered residual # ----------------- -function residual(physicalmodel::ViscoElectricModel, ::Type{Mechano}, (u, φ), v, dΩ, Λ, Δt, un, A) +function residual(physicalmodel::ViscoElectricModel, ::Type{Mechano}, kine::NTuple{2,KinematicModel}, (u, φ), v, dΩ, Λ, Δt, un, A) DΨ = physicalmodel(Λ, Δt=Δt) - F,_,_ = get_Kinematics(physicalmodel.mechano.Kinematic; Λ=Λ) - E = get_Kinematics(physicalmodel.electro.Kinematic; Λ=Λ) + F,_,_ = get_Kinematics(kine[1]; Λ=Λ) + E = get_Kinematics(kine[2]; Λ=Λ) ∂Ψu = DΨ[2] ∫(∇(v)' ⊙ (∂Ψu ∘ (F∘∇(u)', F∘∇(un)', E∘∇(φ), A...)))dΩ end -function residual(physicalmodel::ViscoElectricModel, ::Type{Electro}, (u, φ), vφ, dΩ, Λ, Δt, un, A) +function residual(physicalmodel::ViscoElectricModel, ::Type{Electro}, kine::NTuple{2,KinematicModel}, (u, φ), vφ, dΩ, Λ, Δt, un, A) DΨ = physicalmodel(Λ, Δt=Δt) - F,_,_ = get_Kinematics(physicalmodel.mechano.Kinematic; Λ=Λ) - E = get_Kinematics(physicalmodel.electro.Kinematic; Λ=Λ) + F,_,_ = get_Kinematics(kine[1]; Λ=Λ) + E = get_Kinematics(kine[2]; Λ=Λ) ∂Ψφ = DΨ[3] -1.0*∫(∇(vφ) ⋅ (∂Ψφ ∘ (F∘∇(u)', F∘∇(un)', E∘∇(φ), A...)))dΩ end @@ -97,26 +97,26 @@ end # Stagered jacobian # ----------------- -function jacobian(physicalmodel::ViscoElectricModel, ::Type{Mechano}, (u, φ), du, v, dΩ, Λ, Δt, un, A) +function jacobian(physicalmodel::ViscoElectricModel, ::Type{Mechano}, kine::NTuple{2,KinematicModel}, (u, φ), du, v, dΩ, Λ, Δt, un, A) DΨ = physicalmodel(Λ, Δt=Δt) - F,_,_ = get_Kinematics(physicalmodel.mechano.Kinematic; Λ=Λ) - E = get_Kinematics(physicalmodel.electro.Kinematic; Λ=Λ) + F,_,_ = get_Kinematics(kine[1]; Λ=Λ) + E = get_Kinematics(kine[2]; Λ=Λ) ∂Ψuu = DΨ[4] ∫(∇(v)' ⊙ ((∂Ψuu ∘ (F∘∇(u)', F∘∇(un)', E∘∇(φ), A...)) ⊙ (∇(du)')))dΩ end -function jacobian(physicalmodel::ViscoElectricModel, ::Type{Electro}, (u, φ), dφ, vφ, dΩ, Λ, Δt, un, A) +function jacobian(physicalmodel::ViscoElectricModel, ::Type{Electro}, kine::NTuple{2,KinematicModel}, (u, φ), dφ, vφ, dΩ, Λ, Δt, un, A) DΨ = physicalmodel(Λ, Δt=Δt) - F,_,_ = get_Kinematics(physicalmodel.mechano.Kinematic; Λ=Λ) - E = get_Kinematics(physicalmodel.electro.Kinematic; Λ=Λ) + F,_,_ = get_Kinematics(kine[1]; Λ=Λ) + E = get_Kinematics(kine[2]; Λ=Λ) ∂Ψφφ = DΨ[6] ∫(∇(vφ)' ⋅ ((∂Ψφφ ∘ (F∘∇(u)', F∘∇(un)', E∘∇(φ), A...)) ⋅ ∇(dφ)))dΩ end -function jacobian(physicalmodel::ViscoElectricModel, ::Type{ElectroMechano}, (u, φ), (du, dφ), (v, vφ), dΩ, Λ, Δt, un, A) +function jacobian(physicalmodel::ViscoElectricModel, ::Type{ElectroMechano}, kine::NTuple{2,KinematicModel}, (u, φ), (du, dφ), (v, vφ), dΩ, Λ, Δt, un, A) DΨ = physicalmodel(Λ, Δt=Δt) - F,_,_ = get_Kinematics(physicalmodel.mechano.Kinematic; Λ=Λ) - E = get_Kinematics(physicalmodel.electro.Kinematic; Λ=Λ) + F,_,_ = get_Kinematics(kine[1]; Λ=Λ) + E = get_Kinematics(kine[2]; Λ=Λ) ∂Ψφu = DΨ[5] -1.0*∫(∇(dφ) ⋅ ((∂Ψφu ∘ (F∘∇(u)', F∘∇(un)', E∘∇(φ), A...)) ⊙ (∇(v)')))dΩ - ∫(∇(vφ) ⋅ ((∂Ψφu ∘ (F∘∇(u)', F∘∇(un)', E∘∇(φ), A...)) ⊙ (∇(du)')))dΩ @@ -126,13 +126,13 @@ end # Monolithic strategy # ------------------- -function residual(physicalmodel::ViscoElectricModel, (u, φ), (v, vφ), dΩ, Λ, Δt, un, A) - residual(physicalmodel, Mechano, (u, φ), v, dΩ, Λ, Δt, un, A) + - residual(physicalmodel, Electro, (u, φ), vφ, dΩ, Λ, Δt, un, A) +function residual(physicalmodel::ViscoElectricModel, kine::NTuple{2,KinematicModel}, (u, φ), (v, vφ), dΩ, Λ, Δt, un, A) + residual(physicalmodel, Mechano, kine, (u, φ), v, dΩ, Λ, Δt, un, A) + + residual(physicalmodel, Electro, kine, (u, φ), vφ, dΩ, Λ, Δt, un, A) end -function jacobian(physicalmodel::ViscoElectricModel, (u, φ), (du, dφ), (v, vφ), dΩ, Λ, Δt, un, A) - jacobian(physicalmodel, Mechano, (u, φ), du, v, dΩ, Λ, Δt, un, A) + - jacobian(physicalmodel, Electro, (u, φ), dφ, vφ, dΩ, Λ, Δt, un, A) + - jacobian(physicalmodel, ElectroMechano, (u, φ), (du, dφ), (v, vφ), dΩ, Λ, Δt, un, A) +function jacobian(physicalmodel::ViscoElectricModel, kine::NTuple{2,KinematicModel}, (u, φ), (du, dφ), (v, vφ), dΩ, Λ, Δt, un, A) + jacobian(physicalmodel, Mechano, kine,(u, φ), du, v, dΩ, Λ, Δt, un, A) + + jacobian(physicalmodel, Electro, kine,(u, φ), dφ, vφ, dΩ, Λ, Δt, un, A) + + jacobian(physicalmodel, ElectroMechano, kine,(u, φ), (du, dφ), (v, vφ), dΩ, Λ, Δt, un, A) end diff --git a/src/WeakForms/WeakForms.jl b/src/WeakForms/WeakForms.jl index 75c6e2c..6e73622 100644 --- a/src/WeakForms/WeakForms.jl +++ b/src/WeakForms/WeakForms.jl @@ -36,9 +36,9 @@ end Calculate the residual using the given constitutive model and finite element functions. """ -function residual(physicalmodel::Mechano, u, v, dΩ, Λ=1.0, Δt=0.0, vars...) +function residual(physicalmodel::Mechano, km::KinematicModel ,u, v, dΩ, Λ=1.0, Δt=0.0, vars...) _, ∂Ψu, _ = physicalmodel(Λ) - F, _, _ = get_Kinematics(physicalmodel.Kinematic; Λ=Λ) + F, _, _ = get_Kinematics(km; Λ=Λ) ∫(∇(v)' ⊙ (∂Ψu ∘ (F∘∇(u)')))dΩ end @@ -47,21 +47,21 @@ end Calculate the jacobian using the given constitutive model and finite element functions. """ -function jacobian(physicalmodel::Mechano, u, du, v, dΩ, Λ=1.0, Δt=0.0, vars...) +function jacobian(physicalmodel::Mechano, km::KinematicModel, u, du, v, dΩ, Λ=1.0, Δt=0.0, vars...) _, _, ∂Ψuu = physicalmodel(Λ) - F, _, _ = get_Kinematics(physicalmodel.Kinematic; Λ=Λ) + F, _, _ = get_Kinematics(km; Λ=Λ) ∫(∇(v)' ⊙ ((∂Ψuu ∘ (F∘∇(u)')) ⊙ ∇(du)'))dΩ end -function residual(physicalmodel::ViscoElastic, u, v, dΩ, t, Δt, un, A) +function residual(physicalmodel::ViscoElastic, km::KinematicModel,u, v, dΩ, t, Δt, un, A) _, ∂Ψu, _ = physicalmodel(t, Δt=Δt) - F, _, _ = get_Kinematics(physicalmodel.Kinematic, Λ=t) + F, _, _ = get_Kinematics(km, Λ=t) ∫(∇(v)' ⊙ (∂Ψu ∘ (F∘∇(u)', F∘∇(un)', A...)))dΩ end -function jacobian(physicalmodel::ViscoElastic, u, du, v, dΩ, t, Δt, un, A) +function jacobian(physicalmodel::ViscoElastic, km::KinematicModel, u, du, v, dΩ, t, Δt, un, A) _, _, ∂Ψuu = physicalmodel(t, Δt=Δt) - F, _, _ = get_Kinematics(physicalmodel.Kinematic, Λ=t) + F, _, _ = get_Kinematics(km, Λ=t) ∫(∇(v)' ⊙ (inner ∘ (∂Ψuu∘(F∘∇(u)', F∘∇(un)', A...), ∇(du)')))dΩ end @@ -81,94 +81,94 @@ end # Stagered strategy # ----------------- -function residual(physicalmodel::ThermoElectroMechano, ::Type{Mechano}, (u, φ, θ), v, dΩ, Λ=1.0) +function residual(physicalmodel::ThermoElectroMechano, ::Type{Mechano}, kine::NTuple{3,KinematicModel},(u, φ, θ), v, dΩ, Λ=1.0) DΨ= physicalmodel(Λ) - F,_,_ = get_Kinematics(physicalmodel.mechano.Kinematic; Λ=Λ) - E = get_Kinematics(physicalmodel.electro.Kinematic; Λ=Λ) + F,_,_ = get_Kinematics(kine[1]; Λ=Λ) + E = get_Kinematics(kine[2]; Λ=Λ) ∂Ψu=DΨ[2] return ∫(∇(v)' ⊙ (∂Ψu ∘ (F∘(∇(u)'), E∘(∇(φ)), θ)))dΩ end -function residual(physicalmodel::ThermoElectroMechano, ::Type{Electro}, (u, φ, θ), vφ, dΩ, Λ=1.0) +function residual(physicalmodel::ThermoElectroMechano, ::Type{Electro}, kine::NTuple{3,KinematicModel}, (u, φ, θ), vφ, dΩ, Λ=1.0) DΨ= physicalmodel(Λ) - F,_,_ = get_Kinematics(physicalmodel.mechano.Kinematic; Λ=Λ) - E = get_Kinematics(physicalmodel.electro.Kinematic; Λ=Λ) + F,_,_ = get_Kinematics(kine[1]; Λ=Λ) + E = get_Kinematics(kine[2]; Λ=Λ) ∂Ψφ=DΨ[3] return -1.0*∫(∇(vφ)' ⋅ (∂Ψφ ∘ (F∘(∇(u)'), E∘(∇(φ)), θ)))dΩ end -function residual(physicalmodel::ThermoElectroMechano, ::Type{Thermo}, (u, φ, θ), vθ, dΩ, Λ=1.0) +function residual(physicalmodel::ThermoElectroMechano, ::Type{Thermo}, kine::NTuple{3,KinematicModel}, (u, φ, θ), vθ, dΩ, Λ=1.0) κ=physicalmodel.thermo.κ return ∫(κ * ∇(θ) ⋅ ∇(vθ))dΩ end -function jacobian(physicalmodel::ThermoElectroMechano, ::Type{Mechano}, (u, φ, θ), du, v, dΩ, Λ=1.0) +function jacobian(physicalmodel::ThermoElectroMechano, ::Type{Mechano}, kine::NTuple{3,KinematicModel}, (u, φ, θ), du, v, dΩ, Λ=1.0) DΨ= physicalmodel(Λ) - F,_,_ = get_Kinematics(physicalmodel.mechano.Kinematic; Λ=Λ) - E = get_Kinematics(physicalmodel.electro.Kinematic; Λ=Λ) + F,_,_ = get_Kinematics(kine[1]; Λ=Λ) + E = get_Kinematics(kine[2]; Λ=Λ) ∂Ψuu=DΨ[5] ∫(∇(v)' ⊙ ((∂Ψuu ∘ (F∘(∇(u)'), E∘(∇(φ)), θ)) ⊙ (∇(du)')))dΩ end -function jacobian(physicalmodel::ThermoElectroMechano, ::Type{Electro}, (u, φ, θ), dφ, vφ, dΩ, Λ=1.0) +function jacobian(physicalmodel::ThermoElectroMechano, ::Type{Electro}, kine::NTuple{3,KinematicModel}, (u, φ, θ), dφ, vφ, dΩ, Λ=1.0) DΨ= physicalmodel(Λ) - F,_,_ = get_Kinematics(physicalmodel.mechano.Kinematic; Λ=Λ) - E = get_Kinematics(physicalmodel.electro.Kinematic; Λ=Λ) + F,_,_ = get_Kinematics(kine[1]; Λ=Λ) + E = get_Kinematics(kine[2]; Λ=Λ) ∂Ψφφ=DΨ[6] ∫(∇(vφ) ⋅ ((∂Ψφφ ∘ (F∘(∇(u)'), E∘(∇(φ)), θ)) ⋅ ∇(dφ)))dΩ end -function jacobian(physicalmodel::ThermoElectroMechano, ::Type{Thermo}, dθ, vθ, dΩ, Λ=1.0) +function jacobian(physicalmodel::ThermoElectroMechano, ::Type{Thermo}, kine::NTuple{3,KinematicModel}, dθ, vθ, dΩ, Λ=1.0) κ=physicalmodel.thermo.κ ∫(κ * ∇(dθ) ⋅ ∇(vθ))dΩ end -function jacobian(physicalmodel::ThermoElectroMechano, ::Type{Thermo}, (u, φ, θ)::Tuple, dθ, vθ, dΩ, Λ=1.0) +function jacobian(physicalmodel::ThermoElectroMechano, ::Type{Thermo}, kine::NTuple{3,KinematicModel}, (u, φ, θ)::Tuple, dθ, vθ, dΩ, Λ=1.0) κ=physicalmodel.thermo.κ ∫((κ ∘ (u, φ, θ)) * ∇(dθ) ⋅ ∇(vθ))dΩ end -function jacobian(physicalmodel::ThermoElectroMechano, ::Type{ElectroMechano}, (u, φ, θ), (du, dφ), (v,vφ), dΩ, Λ=1.0) +function jacobian(physicalmodel::ThermoElectroMechano, ::Type{ElectroMechano}, kine::NTuple{3,KinematicModel},(u, φ, θ), (du, dφ), (v,vφ), dΩ, Λ=1.0) DΨ= physicalmodel(Λ) - F,_,_ = get_Kinematics(physicalmodel.mechano.Kinematic; Λ=Λ) - E = get_Kinematics(physicalmodel.electro.Kinematic; Λ=Λ) + F,_,_ = get_Kinematics(kine[1]; Λ=Λ) + E = get_Kinematics(kine[2]; Λ=Λ) ∂Ψφu=DΨ[8] -1.0*∫(∇(dφ) ⋅ ((∂Ψφu ∘ (F∘(∇(u)'), E∘(∇(φ)), θ)) ⊙ (∇(v)')))dΩ - ∫(∇(vφ) ⋅ ((∂Ψφu ∘ (F∘(∇(u)'), E∘(∇(φ)), θ)) ⊙ (∇(du)')))dΩ end -function jacobian(physicalmodel::ThermoElectroMechano, ::Type{ThermoMechano}, (u, φ, θ), dθ, v, dΩ, Λ=1.0) +function jacobian(physicalmodel::ThermoElectroMechano, ::Type{ThermoMechano},kine::NTuple{3,KinematicModel}, (u, φ, θ), dθ, v, dΩ, Λ=1.0) DΨ= physicalmodel(Λ) - F,_,_ = get_Kinematics(physicalmodel.mechano.Kinematic; Λ=Λ) - E = get_Kinematics(physicalmodel.electro.Kinematic; Λ=Λ) + F,_,_ = get_Kinematics(kine[1]; Λ=Λ) + E = get_Kinematics(kine[2]; Λ=Λ) ∂Ψuθ=DΨ[9] ∫(∇(v)' ⊙ (∂Ψuθ ∘ (F∘(∇(u)'), E∘(∇(φ)), θ)) * dθ)dΩ end -function jacobian(physicalmodel::ThermoElectroMechano, ::Type{ThermoElectro}, (u, φ, θ), dθ, vφ, dΩ, Λ=1.0) +function jacobian(physicalmodel::ThermoElectroMechano, ::Type{ThermoElectro}, kine::NTuple{3,KinematicModel},(u, φ, θ), dθ, vφ, dΩ, Λ=1.0) DΨ= physicalmodel(Λ) - F,_,_ = get_Kinematics(physicalmodel.mechano.Kinematic; Λ=Λ) - E = get_Kinematics(physicalmodel.electro.Kinematic; Λ=Λ) + F,_,_ = get_Kinematics(kine[1]; Λ=Λ) + E = get_Kinematics(kine[2]; Λ=Λ) ∂Ψφθ=DΨ[10] -1.0*∫(∇(vφ) ⋅ ((∂Ψφθ ∘ (F∘(∇(u)'), E∘(∇(φ)), θ)) * dθ))dΩ end # Monolithic strategy # ------------------- -function residual(physicalmodel::ThermoElectroMechano, (u, φ, θ), (v, vφ, vθ), dΩ, Λ=1.0) - residual(physicalmodel, Mechano, (u, φ, θ), v, dΩ, Λ) + - residual(physicalmodel, Electro, (u, φ, θ), vφ, dΩ, Λ) + - residual(physicalmodel, Thermo, (u, φ, θ), vθ, dΩ, Λ) +function residual(physicalmodel::ThermoElectroMechano, kine::NTuple{3,KinematicModel}, (u, φ, θ), (v, vφ, vθ), dΩ, Λ=1.0) + residual(physicalmodel, Mechano, kine, (u, φ, θ), v, dΩ, Λ) + + residual(physicalmodel, Electro, kine, (u, φ, θ), vφ, dΩ, Λ) + + residual(physicalmodel, Thermo, kine,(u, φ, θ), vθ, dΩ, Λ) end -function jacobian(physicalmodel::ThermoElectroMechano, (u, φ, θ), (du, dφ, dθ), (v, vφ, vθ), dΩ, Λ=1.0) - jacobian(physicalmodel, Mechano, (u, φ, θ), du, v, dΩ, Λ)+ - jacobian(physicalmodel, Electro, (u, φ, θ), dφ, vφ, dΩ, Λ)+ - jacobian(physicalmodel, Thermo, dθ, vθ, dΩ, Λ)+ - jacobian(physicalmodel, ElectroMechano, (u, φ, θ), (du, dφ), (v,vφ), dΩ, Λ)+ - jacobian(physicalmodel, ThermoMechano, (u, φ, θ), dθ, v, dΩ, Λ)+ - jacobian(physicalmodel, ThermoElectro, (u, φ, θ), dθ, vφ, dΩ, Λ) +function jacobian(physicalmodel::ThermoElectroMechano, kine::NTuple{3,KinematicModel}, (u, φ, θ), (du, dφ, dθ), (v, vφ, vθ), dΩ, Λ=1.0) + jacobian(physicalmodel, Mechano, kine, (u, φ, θ), du, v, dΩ, Λ)+ + jacobian(physicalmodel, Electro, kine, (u, φ, θ), dφ, vφ, dΩ, Λ)+ + jacobian(physicalmodel, Thermo, kine, dθ, vθ, dΩ, Λ)+ + jacobian(physicalmodel, ElectroMechano, kine, (u, φ, θ), (du, dφ), (v,vφ), dΩ, Λ)+ + jacobian(physicalmodel, ThermoMechano, kine, (u, φ, θ), dθ, v, dΩ, Λ)+ + jacobian(physicalmodel, ThermoElectro, kine, (u, φ, θ), dθ, vφ, dΩ, Λ) end @@ -178,53 +178,53 @@ end # Stagered strategy # ----------------- -function residual(physicalmodel::ThermoMechano, ::Type{Mechano}, (u, θ), v, dΩ, Λ=1.0) +function residual(physicalmodel::ThermoMechano, ::Type{Mechano}, kine::NTuple{2,KinematicModel}, (u, θ), v, dΩ, Λ=1.0) DΨ= physicalmodel(Λ) - F,_,_ = get_Kinematics(physicalmodel.mechano.Kinematic; Λ=Λ) + F,_,_ = get_Kinematics(kine[1]; Λ=Λ) ∂Ψu=DΨ[2] ∫(∇(v)' ⊙ (∂Ψu ∘ (F∘(∇(u)'), θ)))dΩ end -function residual(physicalmodel::ThermoMechano, ::Type{Thermo}, (u, θ), vθ, dΩ, Λ=1.0) +function residual(physicalmodel::ThermoMechano, ::Type{Thermo}, kine::NTuple{2,KinematicModel}, (u, θ), vθ, dΩ, Λ=1.0) κ=physicalmodel.thermo.κ ∫(κ * ∇(θ) ⋅ ∇(vθ))dΩ end -function jacobian(physicalmodel::ThermoMechano, ::Type{Mechano}, (u, θ), du, v, dΩ, Λ=1.0) +function jacobian(physicalmodel::ThermoMechano, ::Type{Mechano}, kine::NTuple{2,KinematicModel}, (u, θ), du, v, dΩ, Λ=1.0) DΨ= physicalmodel(Λ) - F,_,_ = get_Kinematics(physicalmodel.mechano.Kinematic; Λ=Λ) + F,_,_ = get_Kinematics(kine[1]; Λ=Λ) ∂Ψuu=DΨ[4] ∫(∇(v)' ⊙ ((∂Ψuu ∘ (F∘(∇(u)'), θ)) ⊙ (∇(du)')))dΩ end -function jacobian(physicalmodel::ThermoMechano, ::Type{Thermo}, dθ, vθ, dΩ, Λ=1.0) +function jacobian(physicalmodel::ThermoMechano, ::Type{Thermo}, kine::NTuple{2,KinematicModel},dθ, vθ, dΩ, Λ=1.0) κ=physicalmodel.thermo.κ ∫(κ * ∇(dθ) ⋅ ∇(vθ))dΩ end -function jacobian(physicalmodel::ThermoMechano, ::Type{Thermo}, (u, θ)::Tuple, dθ, vθ, dΩ, Λ=1.0) +function jacobian(physicalmodel::ThermoMechano, ::Type{Thermo}, kine::NTuple{2,KinematicModel}, (u, θ)::Tuple, dθ, vθ, dΩ, Λ=1.0) κ=physicalmodel.thermo.κ ∫((κ ∘ (u, θ)) * ∇(dθ) ⋅ ∇(vθ))dΩ end -function jacobian(physicalmodel::ThermoMechano, ::Type{ThermoMechano}, (u, θ), (du, dθ), v, dΩ, Λ=1.0) +function jacobian(physicalmodel::ThermoMechano, ::Type{ThermoMechano}, kine::NTuple{2,KinematicModel}, (u, θ), (du, dθ), v, dΩ, Λ=1.0) DΨ= physicalmodel(Λ) - F,_,_ = get_Kinematics(physicalmodel.mechano.Kinematic; Λ=Λ) + F,_,_ = get_Kinematics(kine[1]; Λ=Λ) ∂Ψuθ=DΨ[6] ∫(∇(v)' ⊙ (∂Ψuθ ∘ (F∘(∇(u)'), θ)) * dθ)dΩ end # Monolithic strategy # ------------------- -function residual(physicalmodel::ThermoMechano, (u, θ), (v, vθ), dΩ, Λ=1.0) - residual(physicalmodel, Mechano, (u, θ), v, dΩ, Λ) + - residual(physicalmodel, Thermo, (u, θ), vθ, dΩ, Λ) +function residual(physicalmodel::ThermoMechano, kine::NTuple{2,KinematicModel},(u, θ), (v, vθ), dΩ, Λ=1.0) + residual(physicalmodel, Mechano, kine, (u, θ), v, dΩ, Λ) + + residual(physicalmodel, Thermo, kine, (u, θ), vθ, dΩ, Λ) end -function jacobian(physicalmodel::ThermoMechano, (u, θ), (du, dθ), (v, vθ), dΩ, Λ=1.0) - jacobian(physicalmodel, Mechano, (u, θ), du, v, dΩ, Λ)+ - jacobian(physicalmodel, Thermo, dθ, vθ, dΩ, Λ)+ - jacobian(physicalmodel, ThermoMechano, (u, θ), (du, dθ), v, dΩ, Λ) +function jacobian(physicalmodel::ThermoMechano, kine::NTuple{2,KinematicModel}, (u, θ), (du, dθ), (v, vθ), dΩ, Λ=1.0) + jacobian(physicalmodel, Mechano, kine, (u, θ), du, v, dΩ, Λ)+ + jacobian(physicalmodel, Thermo, kine,dθ, vθ, dΩ, Λ)+ + jacobian(physicalmodel, ThermoMechano, kine,(u, θ), (du, dθ), v, dΩ, Λ) end # =================== @@ -242,60 +242,60 @@ include("ElectroMechanics.jl") # ----------------- -function residual(physicalmodel::FlexoElectro, ::Type{Mechano}, (u, φ, ϕ₁, ϕ₂, ϕ₃), v, dΩ, X, Λ=1.0) +function residual(physicalmodel::FlexoElectro, ::Type{Mechano}, kine::NTuple{2,KinematicModel}, (u, φ, ϕ₁, ϕ₂, ϕ₃), v, dΩ, X, Λ=1.0) DΨ= physicalmodel(Λ) κ=physicalmodel.κ - F,_,_ = get_Kinematics(physicalmodel.electromechano.mechano.Kinematic; Λ=Λ) - E = get_Kinematics(physicalmodel.electromechano.electro.Kinematic; Λ=Λ) + F,_,_ = get_Kinematics(kine[1]; Λ=Λ) + E = get_Kinematics(kine[2]; Λ=Λ) ∂Ψu=DΨ[2] Φ=DΨ[7] ∫((∇(v)' ⊙ (∂Ψu ∘ (F∘(∇(u)',X), E∘(∇(φ)))+κ*(∇(u)'-(Φ(ϕ₁,ϕ₂,ϕ₃))))))dΩ end -function residual(physicalmodel::FlexoElectro, ::Type{Electro}, (u, φ), vφ, dΩ, X, Λ=1.0) +function residual(physicalmodel::FlexoElectro, ::Type{Electro}, kine::NTuple{2,KinematicModel}, (u, φ), vφ, dΩ, X, Λ=1.0) DΨ= physicalmodel(Λ) - F,_,_ = get_Kinematics(physicalmodel.electromechano.mechano.Kinematic; Λ=Λ) - E = get_Kinematics(physicalmodel.electromechano.Electro.Kinematic; Λ=Λ) + F,_,_ = get_Kinematics(kine[1]; Λ=Λ) + E = get_Kinematics(kine[2]; Λ=Λ) ∂Ψφ=DΨ[3] -1.0*∫((∇(vφ) ⋅ (∂Ψφ ∘ (F∘(∇(u)',X), E∘(∇(φ))))))dΩ end -function residual(physicalmodel::FlexoElectro, ::Type{FlexoElectro}, (u, ϕ₁, ϕ₂, ϕ₃), (δϕ₁, δϕ₂, δϕ₃), dΩ, X, Λ=1.0) +function residual(physicalmodel::FlexoElectro, ::Type{FlexoElectro}, kine::NTuple{2,KinematicModel},(u, ϕ₁, ϕ₂, ϕ₃), (δϕ₁, δϕ₂, δϕ₃), dΩ, X, Λ=1.0) DΨ= physicalmodel(Λ) κ=physicalmodel.κ Φ=DΨ[7] ∫((κ⋅(-∇(u)'+(Φ(ϕ₁,ϕ₂,ϕ₃)))) ⊙ (Φ(δϕ₁,δϕ₂,δϕ₃)))dΩ end -function jacobian(physicalmodel::FlexoElectro, ::Type{Mechano}, (u, φ), du, v, dΩ, X, Λ=1.0) +function jacobian(physicalmodel::FlexoElectro, ::Type{Mechano}, kine::NTuple{2,KinematicModel},(u, φ), du, v, dΩ, X, Λ=1.0) DΨ= physicalmodel(Λ) κ=physicalmodel.κ - F,_,_ = get_Kinematics(physicalmodel.electroMechano.mechano.Kinematic; Λ=Λ) - E = get_Kinematics(physicalmodel.electroMechano.electro.Kinematic; Λ=Λ) + F,_,_ = get_Kinematics(kine[1]; Λ=Λ) + E = get_Kinematics(kine[2]; Λ=Λ) ∂Ψuu=DΨ[4] ∫(∇(v)' ⊙ ((∂Ψuu ∘ (F∘(∇(u)',X), E∘(∇(φ)))) ⊙ (∇(du)')))dΩ+ ∫(κ⋅(∇(v)'⊙ ∇(du)'))dΩ end -function jacobian(physicalmodel::FlexoElectro, ::Type{Electro}, (u, φ), dφ, vφ, dΩ, X, Λ=1.0) +function jacobian(physicalmodel::FlexoElectro, ::Type{Electro}, kine::NTuple{2,KinematicModel},(u, φ), dφ, vφ, dΩ, X, Λ=1.0) DΨ= physicalmodel(Λ) - F,_,_ = get_Kinematics(physicalmodel.electromechano.mechano.Kinematic; Λ=Λ) - E = get_Kinematics(physicalmodel.electromechano.electro.Kinematic; Λ=Λ) + F,_,_ = get_Kinematics(kine[1]; Λ=Λ) + E = get_Kinematics(kine[2]; Λ=Λ) ∂Ψφφ=DΨ[6] ∫(∇(vφ)' ⋅ ((∂Ψφφ ∘ (F∘(∇(u)',X), E∘(∇(φ)))) ⋅ ∇(dφ)))dΩ end -function jacobian(physicalmodel::FlexoElectro, ::Type{ElectroMechano}, (u, φ), (du, dφ), (v, vφ), dΩ, X, Λ=1.0) +function jacobian(physicalmodel::FlexoElectro, ::Type{ElectroMechano}, kine::NTuple{2,KinematicModel},(u, φ), (du, dφ), (v, vφ), dΩ, X, Λ=1.0) DΨ= physicalmodel(Λ) - F,_,_ = get_Kinematics(physicalmodel.electromechano.mechano.Kinematic; Λ=Λ) - E = get_Kinematics(physicalmodel.electromechano.electro.Kinematic; Λ=Λ) + F,_,_ = get_Kinematics(kine[1]; Λ=Λ) + E = get_Kinematics(kine[2]; Λ=Λ) ∂Ψφu=DΨ[5] -1.0*∫(∇(dφ) ⋅ ((∂Ψφu ∘ (F∘(∇(u)',X), E∘(∇(φ)))) ⊙ (∇(v)')))dΩ - ∫(∇(vφ) ⋅ ((∂Ψφu ∘ (F∘(∇(u)',X), E∘(∇(φ)))) ⊙ (∇(du)')))dΩ end -function jacobian(physicalmodel::FlexoElectro, ::Type{FlexoElectro}, (ϕ₁,ϕ₂,ϕ₃), (du, dϕ₁,dϕ₂,dϕ₃), (v, δϕ₁,δϕ₂,δϕ₃), dΩ, X, Λ=1.0) +function jacobian(physicalmodel::FlexoElectro, ::Type{FlexoElectro}, kine::NTuple{2,KinematicModel},(ϕ₁,ϕ₂,ϕ₃), (du, dϕ₁,dϕ₂,dϕ₃), (v, δϕ₁,δϕ₂,δϕ₃), dΩ, X, Λ=1.0) DΨ= physicalmodel(Λ) κ=physicalmodel.κ Φ=DΨ[7] @@ -308,17 +308,17 @@ end # Monolithic strategy # ------------------- -function residual(physicalmodel::FlexoElectro, (u, φ, ϕ₁, ϕ₂, ϕ₃), (v, vφ, δϕ₁, δϕ₂, δϕ₃), dΩ, X, Λ=1.0) - residual(physicalmodel, Mechano, (u, φ, ϕ₁, ϕ₂, ϕ₃), v, dΩ, X, Λ) + - residual(physicalmodel, Electro, (u, φ), vφ, dΩ, X, Λ)+ - residual(physicalmodel, FlexoElectro, (u, ϕ₁, ϕ₂, ϕ₃), (δϕ₁, δϕ₂, δϕ₃), dΩ, X, Λ) +function residual(physicalmodel::FlexoElectro, kine::NTuple{2,KinematicModel},(u, φ, ϕ₁, ϕ₂, ϕ₃), (v, vφ, δϕ₁, δϕ₂, δϕ₃), dΩ, X, Λ=1.0) + residual(physicalmodel, Mechano, kine,(u, φ, ϕ₁, ϕ₂, ϕ₃), v, dΩ, X, Λ) + + residual(physicalmodel, Electro, kine,(u, φ), vφ, dΩ, X, Λ)+ + residual(physicalmodel, FlexoElectro, kine,(u, ϕ₁, ϕ₂, ϕ₃), (δϕ₁, δϕ₂, δϕ₃), dΩ, X, Λ) end -function jacobian(physicalmodel::FlexoElectro, (u, φ, ϕ₁,ϕ₂,ϕ₃), (du, dφ, dϕ₁,dϕ₂,dϕ₃), (v, vφ, δϕ₁,δϕ₂,δϕ₃), dΩ, X, Λ=1.0) - jacobian(physicalmodel, Mechano, (u, φ), du, v, dΩ, X, Λ)+ - jacobian(physicalmodel, Electro, (u, φ), dφ, vφ, dΩ, X, Λ)+ - jacobian(physicalmodel, ElectroMechano, (u, φ), (du, dφ), (v, vφ), dΩ, X, Λ)+ - jacobian(physicalmodel, FlexoElectro, (ϕ₁,ϕ₂,ϕ₃), (du, dϕ₁,dϕ₂,dϕ₃), (v, δϕ₁,δϕ₂,δϕ₃), dΩ, X, Λ) +function jacobian(physicalmodel::FlexoElectro, kine::NTuple{2,KinematicModel},(u, φ, ϕ₁,ϕ₂,ϕ₃), (du, dφ, dϕ₁,dϕ₂,dϕ₃), (v, vφ, δϕ₁,δϕ₂,δϕ₃), dΩ, X, Λ=1.0) + jacobian(physicalmodel, Mechano, kine,(u, φ), du, v, dΩ, X, Λ)+ + jacobian(physicalmodel, Electro, kine,(u, φ), dφ, vφ, dΩ, X, Λ)+ + jacobian(physicalmodel, ElectroMechano, kine,(u, φ), (du, dφ), (v, vφ), dΩ, X, Λ)+ + jacobian(physicalmodel, FlexoElectro, kine,(ϕ₁,ϕ₂,ϕ₃), (du, dϕ₁,dϕ₂,dϕ₃), (v, δϕ₁,δϕ₂,δϕ₃), dΩ, X, Λ) end @@ -331,108 +331,94 @@ end # Stagered strategy # ----------------- -function residual(physicalmodel::ThermoElectroMech_PINNs, ::Type{Mechano}, (u, φ, θ), v, dΩ, Λ=1.0) +function residual(physicalmodel::ThermoElectroMech_PINNs, ::Type{Mechano}, kine::NTuple{2,KinematicModel}, (u, φ, θ), v, dΩ, Λ=1.0) DΨ= physicalmodel(Λ) - Kinematic_mec = Kinematics(Mechano) - Kinematic_elec = Kinematics(Electro) - F, _, _ = get_Kinematics(Kinematic_mec) - E = get_Kinematics(Kinematic_elec) + F, _, _ = get_Kinematics(kine[1]) + E = get_Kinematics(kine[2]) ∂Ψu=DΨ[2] return ∫(∇(v)' ⊙ (∂Ψu ∘ (F∘(∇(u)'), E∘(∇(φ)), θ)))dΩ end -function residual(physicalmodel::ThermoElectroMech_PINNs, ::Type{Electro}, (u, φ, θ), vφ, dΩ, Λ=1.0) +function residual(physicalmodel::ThermoElectroMech_PINNs, ::Type{Electro}, kine::NTuple{2,KinematicModel},(u, φ, θ), vφ, dΩ, Λ=1.0) DΨ= physicalmodel(Λ) - Kinematic_mec = Kinematics(Mechano) - Kinematic_elec = Kinematics(Electro) - F, _, _ = get_Kinematics(Kinematic_mec) - E = get_Kinematics(Kinematic_elec) + F, _, _ = get_Kinematics(kine[1]) + E = get_Kinematics(kine[2]) ∂Ψφ=DΨ[3] return -1.0*∫(∇(vφ)' ⋅ (∂Ψφ ∘ (F∘(∇(u)'), E∘(∇(φ)), θ)))dΩ end -function residual(physicalmodel::ThermoElectroMech_PINNs, ::Type{Thermo}, (u, φ, θ), vθ, dΩ, Λ=1.0) +function residual(physicalmodel::ThermoElectroMech_PINNs, ::Type{Thermo}, kine::NTuple{2,KinematicModel},(u, φ, θ), vθ, dΩ, Λ=1.0) κ=physicalmodel.κ return ∫(κ * ∇(θ) ⋅ ∇(vθ))dΩ end -function jacobian(physicalmodel::ThermoElectroMech_PINNs, ::Type{Mechano}, (u, φ, θ), du, v, dΩ, Λ=1.0) +function jacobian(physicalmodel::ThermoElectroMech_PINNs, ::Type{Mechano}, kine::NTuple{2,KinematicModel},(u, φ, θ), du, v, dΩ, Λ=1.0) DΨ= physicalmodel(Λ) - Kinematic_mec = Kinematics(Mechano) - Kinematic_elec = Kinematics(Electro) - F, _, _ = get_Kinematics(Kinematic_mec) - E = get_Kinematics(Kinematic_elec) + F, _, _ = get_Kinematics(kine[1]) + E = get_Kinematics(kine[2]) ∂Ψuu=DΨ[5] ∫(∇(v)' ⊙ ((∂Ψuu ∘ (F∘(∇(u)'), E∘(∇(φ)), θ)) ⊙ (∇(du)')))dΩ end -function jacobian(physicalmodel::ThermoElectroMech_PINNs, ::Type{Electro}, (u, φ, θ), dφ, vφ, dΩ, Λ=1.0) +function jacobian(physicalmodel::ThermoElectroMech_PINNs, ::Type{Electro}, kine::NTuple{2,KinematicModel},(u, φ, θ), dφ, vφ, dΩ, Λ=1.0) DΨ= physicalmodel(Λ) - Kinematic_mec = Kinematics(Mechano) - Kinematic_elec = Kinematics(Electro) - F, _, _ = get_Kinematics(Kinematic_mec) - E = get_Kinematics(Kinematic_elec) + F, _, _ = get_Kinematics(kine[1]) + E = get_Kinematics(kine[2]) ∂Ψφφ=DΨ[6] ∫(∇(vφ) ⋅ ((∂Ψφφ ∘ (F∘(∇(u)'), E∘(∇(φ)), θ)) ⋅ ∇(dφ)))dΩ end -function jacobian(physicalmodel::ThermoElectroMech_PINNs, ::Type{Thermo}, dθ, vθ, dΩ, Λ=1.0) +function jacobian(physicalmodel::ThermoElectroMech_PINNs, ::Type{Thermo}, kine::NTuple{2,KinematicModel},dθ, vθ, dΩ, Λ=1.0) κ=physicalmodel.κ ∫(κ * ∇(dθ) ⋅ ∇(vθ))dΩ end -function jacobian(physicalmodel::ThermoElectroMech_PINNs, ::Type{Thermo}, (u, φ, θ)::Tuple, dθ, vθ, dΩ, Λ=1.0) +function jacobian(physicalmodel::ThermoElectroMech_PINNs, ::Type{Thermo}, kine::NTuple{2,KinematicModel},(u, φ, θ)::Tuple, dθ, vθ, dΩ, Λ=1.0) κ=physicalmodel.κ ∫((κ ∘ (u, φ, θ)) * ∇(dθ) ⋅ ∇(vθ))dΩ end -function jacobian(physicalmodel::ThermoElectroMech_PINNs, ::Type{ElectroMechano}, (u, φ, θ), (du, dφ), (v,vφ), dΩ, Λ=1.0) +function jacobian(physicalmodel::ThermoElectroMech_PINNs, ::Type{ElectroMechano}, kine::NTuple{2,KinematicModel},(u, φ, θ), (du, dφ), (v,vφ), dΩ, Λ=1.0) DΨ= physicalmodel(Λ) - Kinematic_mec = Kinematics(Mechano) - Kinematic_elec = Kinematics(Electro) - F, _, _ = get_Kinematics(Kinematic_mec) - E = get_Kinematics(Kinematic_elec) + F, _, _ = get_Kinematics(kine[1]) + E = get_Kinematics(kine[2]) ∂Ψφu=DΨ[8] -1.0*∫(∇(dφ) ⋅ ((∂Ψφu ∘ (F∘(∇(u)'), E∘(∇(φ)), θ)) ⊙ (∇(v)')))dΩ - ∫(∇(vφ) ⋅ ((∂Ψφu ∘ (F∘(∇(u)'), E∘(∇(φ)), θ)) ⊙ (∇(du)')))dΩ end -function jacobian(physicalmodel::ThermoElectroMech_PINNs, ::Type{ThermoMechano}, (u, φ, θ), dθ, v, dΩ, Λ=1.0) +function jacobian(physicalmodel::ThermoElectroMech_PINNs, ::Type{ThermoMechano}, kine::NTuple{2,KinematicModel}, (u, φ, θ), dθ, v, dΩ, Λ=1.0) DΨ= physicalmodel(Λ) - Kinematic_mec = Kinematics(Mechano) - Kinematic_elec = Kinematics(Electro) - F, _, _ = get_Kinematics(Kinematic_mec) - E = get_Kinematics(Kinematic_elec) + F, _, _ = get_Kinematics(kine[1]) + E = get_Kinematics(kine[2]) ∂Ψuθ=DΨ[9] ∫(∇(v)' ⊙ (∂Ψuθ ∘ (F∘(∇(u)'), E∘(∇(φ)), θ)) * dθ)dΩ end -function jacobian(physicalmodel::ThermoElectroMech_PINNs, ::Type{ThermoElectro}, (u, φ, θ), dθ, vφ, dΩ, Λ=1.0) +function jacobian(physicalmodel::ThermoElectroMech_PINNs, ::Type{ThermoElectro}, kine::NTuple{2,KinematicModel}, (u, φ, θ), dθ, vφ, dΩ, Λ=1.0) DΨ= physicalmodel(Λ) - Kinematic_mec = Kinematics(Mechano) - Kinematic_elec = Kinematics(Electro) - F, _, _ = get_Kinematics(Kinematic_mec) - E = get_Kinematics(Kinematic_elec) + F, _, _ = get_Kinematics(kine[1]) + E = get_Kinematics(kine[2]) ∂Ψφθ=DΨ[10] -1.0*∫(∇(vφ) ⋅ ((∂Ψφθ ∘ (F∘(∇(u)'), E∘(∇(φ)), θ)) * dθ))dΩ end # Monolithic strategy # ------------------- -function residual(physicalmodel::ThermoElectroMech_PINNs, (u, φ, θ), (v, vφ, vθ), dΩ, Λ=1.0) - residual(physicalmodel, Mechano, (u, φ, θ), v, dΩ, Λ) + - residual(physicalmodel, Electro, (u, φ, θ), vφ, dΩ, Λ) + - residual(physicalmodel, Thermo, (u, φ, θ), vθ, dΩ, Λ) -end - -function jacobian(physicalmodel::ThermoElectroMech_PINNs, (u, φ, θ), (du, dφ, dθ), (v, vφ, vθ), dΩ, Λ=1.0) - jacobian(physicalmodel, Mechano, (u, φ, θ), du, v, dΩ, Λ)+ - jacobian(physicalmodel, Electro, (u, φ, θ), dφ, vφ, dΩ, Λ)+ - jacobian(physicalmodel, Thermo, dθ, vθ, dΩ, Λ)+ - jacobian(physicalmodel, ElectroMechano, (u, φ, θ), (du, dφ), (v,vφ), dΩ, Λ)+ - jacobian(physicalmodel, ThermoMechano, (u, φ, θ), dθ, v, dΩ, Λ)+ - jacobian(physicalmodel, ThermoElectro, (u, φ, θ), dθ, vφ, dΩ, Λ) +function residual(physicalmodel::ThermoElectroMech_PINNs, kine::NTuple{2,KinematicModel},(u, φ, θ), (v, vφ, vθ), dΩ, Λ=1.0) + residual(physicalmodel, Mechano, kine,(u, φ, θ), v, dΩ, Λ) + + residual(physicalmodel, Electro, kine,(u, φ, θ), vφ, dΩ, Λ) + + residual(physicalmodel, Thermo, kine,(u, φ, θ), vθ, dΩ, Λ) +end + +function jacobian(physicalmodel::ThermoElectroMech_PINNs, kine::NTuple{2,KinematicModel}, (u, φ, θ), (du, dφ, dθ), (v, vφ, vθ), dΩ, Λ=1.0) + jacobian(physicalmodel, Mechano, kine,(u, φ, θ), du, v, dΩ, Λ)+ + jacobian(physicalmodel, Electro, kine,(u, φ, θ), dφ, vφ, dΩ, Λ)+ + jacobian(physicalmodel, Thermo, kine,dθ, vθ, dΩ, Λ)+ + jacobian(physicalmodel, ElectroMechano, kine,(u, φ, θ), (du, dφ), (v,vφ), dΩ, Λ)+ + jacobian(physicalmodel, ThermoMechano, kine,(u, φ, θ), dθ, v, dΩ, Λ)+ + jacobian(physicalmodel, ThermoElectro, kine,(u, φ, θ), dθ, vφ, dΩ, Λ) end diff --git a/test/TestConstitutiveModels/ElectroMechanicalTests.jl b/test/TestConstitutiveModels/ElectroMechanicalTests.jl index 0826a85..8dcf6be 100644 --- a/test/TestConstitutiveModels/ElectroMechanicalTests.jl +++ b/test/TestConstitutiveModels/ElectroMechanicalTests.jl @@ -1,5 +1,6 @@ using Gridap.TensorValues using HyperFEM.PhysicalModels +using HyperFEM.TensorAlgebra const ∇φ = VectorValue(1.0:3.0...) @@ -30,8 +31,12 @@ const ∇un = TensorValue(1.0:9.0...) * 5e-4 Ψ, ∂Ψ∂F, ∂Ψ∂F∂F = modelelectro() - F, _, _ = get_Kinematics(model1.Kinematic) - E = get_Kinematics(modelID.Kinematic) + + Ke=Kinematics(Electro,Solid) + E = get_Kinematics(Ke) + K=Kinematics(Mechano,Solid) + F, _, _ = get_Kinematics(K) + @test Ψ(F(∇u),E(∇φ) ,(M1,M2,M3,M4)) == -27.513663654827152 @test isapprox(norm(∂Ψ∂F(F(∇u),E(∇φ) ,(M1,M2,M3,M4))) , 47.45420724735932,rtol=1e-14) @test isapprox(norm(∂Ψ∂F∂F(F(∇u),E(∇φ) ,(M1,M2,M3,M4))), 14.707913034885005,rtol=1e-14) @@ -44,8 +49,10 @@ end modelID = IdealDielectric(ε=4.0) modelelectro =modelMR+modelID Ψ, ∂Ψu, ∂Ψφ, ∂Ψuu, ∂Ψφu, ∂Ψφφ = modelelectro() - F, _, _ = get_Kinematics(modelMR.Kinematic) - E = get_Kinematics(modelID.Kinematic) + Ke=Kinematics(Electro,Solid) + E = get_Kinematics(Ke) + K=Kinematics(Mechano,Solid) + F, _, _ = get_Kinematics(K) @test Ψ(F(∇u), E(∇φ)) == -27.514219755428428 @test norm(∂Ψu(F(∇u), E(∇φ))) == 47.42294370458073 @@ -73,14 +80,14 @@ end Kin_mec = EvolutiveKinematics(Mechano; F=(t) -> ((∇u1, x) -> ∇u1 + one(∇u1) + t * ∇umacro + t * (A ⊙ x))) Kin_elec = EvolutiveKinematics(Electro; E=(t) -> ((∇φ) -> -∇φ + t * Emacro)) - physmec = MooneyRivlin3D(λ=10.0, μ1=1.0, μ2=1.0, Kinematic=Kin_mec) - physelec = IdealDielectric(ε=1.0, Kinematic=Kin_elec) + physmec = MooneyRivlin3D(λ=10.0, μ1=1.0, μ2=1.0) + physelec = IdealDielectric(ε=1.0) physmodel = FlexoElectroModel(mechano=physmec, electro=physelec, κ=1000.0) Ψ, ∂Ψu, ∂Ψφ, ∂Ψuu, ∂Ψφu, ∂Ψφφ, Φ = physmodel(1.0) - F, _, _ = get_Kinematics(physmec.Kinematic; Λ=1.0) - E = get_Kinematics(physelec.Kinematic; Λ=1.0) + F, _, _ = get_Kinematics(Kin_mec; Λ=1.0) + E = get_Kinematics(Kin_elec; Λ=1.0) X = VectorValue(2.4, 1.9, 3.3) @test (Ψ(F(∇u1, X), E(∇φ))) == 13.408299698687056 @@ -99,8 +106,10 @@ end visco_elastic = GeneralizedMaxwell(hyper_elastic, viscous_branch1) dielectric = IdealDielectric(ε=1.0) model =dielectric+visco_elastic - F, _, _ = get_Kinematics(model.mechano.Kinematic) - E = get_Kinematics(model.electro.Kinematic) + Ke=Kinematics(Electro,Solid) + E = get_Kinematics(Ke) + K=Kinematics(Mechano,Solid) + F, _, _ = get_Kinematics(K) Uvn = TensorValue(1.,2.,3.,2.,4.,5.,3.,5.,6.) * 2e-4 + I3 Uvn *= det(Uvn)^(-1/3) λvn = 1e-3 @@ -119,8 +128,10 @@ end visco_elastic = GeneralizedMaxwell(hyper_elastic, viscous_branch1, viscous_branch2) dielectric = IdealDielectric(ε=1.0) model =dielectric+visco_elastic - F, _, _ = get_Kinematics(model.mechano.Kinematic) - E = get_Kinematics(model.electro.Kinematic) + Ke=Kinematics(Electro,Solid) + E = get_Kinematics(Ke) + K=Kinematics(Mechano,Solid) + F, _, _ = get_Kinematics(K) Uvn = TensorValue(1.,2.,3.,2.,4.,5.,3.,5.,6.) * 2e-4 + I3 Uvn *= det(Uvn)^(-1/3) λvn = 1e-3 diff --git a/test/TestConstitutiveModels/PhysicalModelTests.jl b/test/TestConstitutiveModels/PhysicalModelTests.jl index ee6ebc1..b3ebde1 100644 --- a/test/TestConstitutiveModels/PhysicalModelTests.jl +++ b/test/TestConstitutiveModels/PhysicalModelTests.jl @@ -20,22 +20,24 @@ const μParams = [6456.9137547089595, 896.4633794151492, 0.7841068624959612, 1.5386288924587603] -function test_derivatives__(model::PhysicalModel, ∇u; rtol, kwargs...) +function test_derivatives__(model::PhysicalModel, K::KinematicModel ,∇u; rtol, kwargs...) Ψ, ∂Ψu, ∂Ψuu = model() ∂Ψu_(F) = TensorValue(ForwardDiff.gradient(Ψ, get_array(F))) ∂Ψuu_(F) = TensorValue(ForwardDiff.hessian(Ψ, get_array(F))) - F, _, _ = get_Kinematics(model.Kinematic) + K=Kinematics(Mechano,Solid) + F, _, _ = get_Kinematics(K) + @test isapprox(∂Ψu(F(∇u)), ∂Ψu_(F(∇u)), rtol=rtol, kwargs...) @test isapprox(∂Ψuu(F(∇u)), ∂Ψuu_(F(∇u)), rtol=rtol, kwargs...) end -function test_derivatives_2D_(model::PhysicalModel; rtol=1e-14, kwargs...) - test_derivatives__(model, ∇u2, rtol=rtol, kwargs...) +function test_derivatives_2D_(model::PhysicalModel, K::KinematicModel; rtol=1e-14, kwargs...) + test_derivatives__(model, K, ∇u2, rtol=rtol, kwargs...) end -function test_derivatives_3D_(model::PhysicalModel; rtol=1e-14, kwargs...) - test_derivatives__(model, ∇u3, rtol=rtol, kwargs...) +function test_derivatives_3D_(model::PhysicalModel, K::KinematicModel; rtol=1e-14, kwargs...) + test_derivatives__(model, K, ∇u3, rtol=rtol, kwargs...) end @@ -64,8 +66,8 @@ end ∂Ψ∂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(model1.Kinematic) + K=Kinematics(Mechano,Solid) + F, _, _ = get_Kinematics(K) @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) @@ -89,7 +91,9 @@ end ∂Ψ∂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) + K=Kinematics(Mechano,Solid) + F, _, _ = get_Kinematics(K) + @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) @@ -119,7 +123,9 @@ end ∂Ψ∂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) + K=Kinematics(Mechano,Solid) + F, _, _ = get_Kinematics(K) + @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) @@ -137,7 +143,8 @@ end Ψ1, ∂Ψu1, ∂Ψuu1 = model1() Ψ2, ∂Ψu2, ∂Ψuu2 = model2() - F, _, _ = get_Kinematics(model.Kinematic) + K=Kinematics(Mechano,Solid) + F, _, _ = get_Kinematics(K) @test Ψ(F(∇u3), N) == Ψ1(F(∇u3)) + Ψ2(F(∇u3), N) @test norm(∂Ψu(F(∇u3), N)) == norm(∂Ψu1(F(∇u3))+∂Ψu2(F(∇u3), N) ) @test norm(∂Ψuu(F(∇u3), N)) == norm(∂Ψuu1(F(∇u3))+∂Ψuu2(F(∇u3), N)) @@ -159,7 +166,8 @@ end Ψ2, ∂Ψu2, ∂Ψuu2 = model2() Ψ3, ∂Ψu3, ∂Ψuu3 = model3() - F, _, _ = get_Kinematics(model.Kinematic) + K=Kinematics(Mechano,Solid) + F, _, _ = get_Kinematics(K) @test Ψ(F(∇u3), [N1,N2]) == Ψ1(F(∇u3)) + Ψ2(F(∇u3), [N1]) + Ψ3(F(∇u3), [N2]) @test isapprox(norm(∂Ψu(F(∇u3), [N1,N2])), norm(∂Ψu1(F(∇u3))+∂Ψu2(F(∇u3), N1)+∂Ψu3(F(∇u3), N2)), rtol=1e-14) @test isapprox(norm(∂Ψuu(F(∇u3), [N1,N2])), norm(∂Ψuu1(F(∇u3))+∂Ψuu2(F(∇u3), N1)+∂Ψuu3(F(∇u3), N2)), rtol=1e-14) @@ -169,13 +177,13 @@ end @testset "NonlinearMooneyRivlin_CV" begin model = NonlinearMooneyRivlin_CV(λ=3.0, μ1=1.0, μ2=1.0, α=2.0, β=1.0, γ=6.0) - test_derivatives_3D_(model, rtol=1e-13) + test_derivatives_3D_(model, Kinematics(Mechano,Solid),rtol=1e-13) end @testset "NonlinearNeoHookean_CV" begin model = NonlinearNeoHookean_CV(λ=3.0, μ=1.0, α=2.0, γ=6.0) - test_derivatives_3D_(model, rtol=1e-13) + test_derivatives_3D_(model, Kinematics(Mechano,Solid),rtol=1e-13) end @@ -183,8 +191,6 @@ end 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) model = IncompressibleNeoHookean3D_2dP(μ=1.0, τ=1.0, Δt=1.0) Ψ, Se, ∂Se = model() - F, H, _ = get_Kinematics(model.Kinematic) - # Se_(Ce) =2*TensorValue(ForwardDiff.gradient(x -> Ψ(x), get_array(Ce))) @@ -201,77 +207,77 @@ end @testset "LinearElasticity2D" begin model = LinearElasticity2D(λ=3.0, μ=1.0) - test_derivatives_2D_(model) + test_derivatives_2D_(model, Kinematics(Mechano,Solid)) end @testset "LinearElasticity3D" begin model = LinearElasticity3D(λ=3.0, μ=1.0) - test_derivatives_3D_(model) + test_derivatives_3D_(model, Kinematics(Mechano,Solid)) end @testset "NeoHookean3D" begin model = NeoHookean3D(λ=3.0, μ=1.0) - test_derivatives_3D_(model, rtol=1e-13) + test_derivatives_3D_(model, Kinematics(Mechano,Solid), rtol=1e-13) end @testset "MooneyRivlin2D" begin model = MooneyRivlin2D(λ=3.0, μ1=1.0, μ2=2.0) - test_derivatives_2D_(model) + test_derivatives_2D_(model, Kinematics(Mechano,Solid)) end @testset "MooneyRivlin3D" begin model = MooneyRivlin3D(λ=3.0, μ1=1.0, μ2=2.0) - test_derivatives_3D_(model, rtol=1e-13) + test_derivatives_3D_(model,Kinematics(Mechano,Solid), rtol=1e-13) end @testset "NonlinearMooneyRivlin2D" begin model = NonlinearMooneyRivlin2D(λ=(μParams[1] + μParams[2]) * 1e2, μ1=μParams[1], μ2=μParams[2], α=μParams[3], β=μParams[4]) - test_derivatives_2D_(model) + test_derivatives_2D_(model, Kinematics(Mechano,Solid)) end @testset "Yeoh3D" begin model = Yeoh3D(λ=3.0, C10=1.0, C20=1.0, C30=1.0) - test_derivatives_3D_(model) + test_derivatives_3D_(model, Kinematics(Mechano,Solid)) end @testset "NonlinearMooneyRivlin2D_CV" begin model = NonlinearMooneyRivlin2D_CV(λ=(μParams[1] + μParams[2]) * 1e2, μ1=μParams[1], μ2=μParams[2], α=μParams[3], β=μParams[4], γ=μParams[4]) - test_derivatives_2D_(model) + test_derivatives_2D_(model, Kinematics(Mechano,Solid)) end @testset "NonlinearMooneyRivlin3D" begin model = NonlinearMooneyRivlin3D(λ=(μParams[1] + μParams[2]) * 1e2, μ1=μParams[1], μ2=μParams[2], α=μParams[3], β=μParams[4]) - test_derivatives_3D_(model, rtol=1e-13) + test_derivatives_3D_(model,Kinematics(Mechano,Solid), rtol=1e-13) end @testset "IncompressibleNeoHookean2D" begin model = IncompressibleNeoHookean2D(λ=(μParams[1] + μParams[2]) * 1e2, μ=μParams[1]) - test_derivatives_2D_(model) + test_derivatives_2D_(model, Kinematics(Mechano,Solid)) end @testset "IncompressibleNeoHookean2D_CV" begin model = IncompressibleNeoHookean2D_CV(λ=(μParams[1] + μParams[2]) * 1e2, μ=μParams[1], γ=3.0) - test_derivatives_2D_(model) + test_derivatives_2D_(model, Kinematics(Mechano,Solid)) end @testset "NonlinearIncompressibleMooneyRivlin2D_CV" begin model = NonlinearIncompressibleMooneyRivlin2D_CV(λ=(μParams[1] + μParams[2]) * 1e2, μ=μParams[1], α=μParams[3], γ=3.0) - test_derivatives_2D_(model) + test_derivatives_2D_(model, Kinematics(Mechano,Solid)) end @testset "EightChain" begin model = EightChain(μ=μParams[1], N=μParams[2]) - test_derivatives_3D_(model, rtol=1e-13) + test_derivatives_3D_(model, Kinematics(Mechano,Solid),rtol=1e-13) end @@ -283,7 +289,8 @@ end model = TransverseIsotropy2D(μ=μParams[5], α=μParams[6], β=μParams[7]) Ψ, ∂Ψu, ∂Ψuu = model() - F, _, J_ = get_Kinematics(model.Kinematic) + K=Kinematics(Mechano,Solid) + F, _, _ = get_Kinematics(K) # ∂Ψu_(F,N) =TensorValue(ForwardDiff.gradient(x -> Ψ(x,get_array(N)), get_array(F))) # ∂Ψuu_(F,N) =TensorValue(ForwardDiff.hessian(x -> Ψ(x,get_array(N)), get_array(F))) @@ -307,7 +314,8 @@ end model = TransverseIsotropy3D(μ=μParams[5], α=μParams[6], β=μParams[7]) Ψ, ∂Ψu, ∂Ψuu = model() - F, _, _ = get_Kinematics(model.Kinematic) + K=Kinematics(Mechano,Solid) + F, _, _ = get_Kinematics(K) @test Ψ(F(∇u), N) == 269927.3350807581 @test norm(∂Ψu(F(∇u), N)) == 947447.8711645481 @test norm(∂Ψuu(F(∇u), N)) == 3.8258646319087776e6 @@ -330,8 +338,10 @@ end df(δθ::Float64)::Float64 = 1.0 modelTEM = ThermoElectroMechModel(modelT, modelID, modelMR, fθ=f, dfdθ=df) Ψ, ∂Ψu, ∂Ψφ, ∂Ψθ, ∂Ψuu, ∂Ψφφ, ∂Ψθθ, ∂Ψφu, ∂Ψuθ, ∂Ψφθ = modelTEM() - F, _, _ = get_Kinematics(modelMR.Kinematic) - E = get_Kinematics(modelID.Kinematic) + K=Kinematics(Mechano,Solid) + F, _, _ = get_Kinematics(K) + Ke=Kinematics(Electro,Solid) + E = get_Kinematics(Ke) @test (Ψ(F(∇u), E(∇φ), θt)) == -95.74389746463744 @test norm(∂Ψu(F(∇u), E(∇φ), θt)) == 185.1315441384458 @@ -355,7 +365,8 @@ end df(δθ::Float64)::Float64 = 1.0 modelTM = ThermoMechModel(modelT, modelMR, fθ=f, dfdθ=df) Ψ, ∂Ψu, ∂Ψθ, ∂Ψuu, ∂Ψθθ, ∂Ψuθ = modelTM() - F, _, _ = get_Kinematics(modelMR.Kinematic) + K=Kinematics(Mechano,Solid) + F, _, _ = get_Kinematics(K) @test (Ψ(F(∇u), θt)) == -2.190116215314799 @test norm(∂Ψu(F(∇u), θt)) == 50.34457217400186 @@ -381,7 +392,8 @@ end consmodel = ThermoMech_EntropicPolyconvex(modterm, modmec, β=β, G=G, ϕ=ϕ, s=s) Ψ, ∂Ψu, ∂Ψθ, ∂Ψuu, ∂Ψθθ, ∂Ψuθ = consmodel() - F, _, _ = get_Kinematics(modmec.Kinematic) + K=Kinematics(Mechano,Solid) + F, _, _ = get_Kinematics(K) @test (Ψ(F(∇u), θt)) == -129.4022076861008 @test norm(∂Ψu(F(∇u), θt)) == 437.9269386687991 @@ -404,8 +416,10 @@ end modelTEM = ThermoElectroMech_Bonet(modelT, modelID, modelMR) Ψ, ∂Ψu, ∂ΨE, ∂Ψθ, ∂ΨFF, ∂ΨEE, ∂2Ψθθ, ∂ΨEF, ∂ΨFθ, ∂ΨEθ, η = modelTEM() - F, _, _ = get_Kinematics(modelMR.Kinematic) - E = get_Kinematics(modelID.Kinematic) + K=Kinematics(Mechano,Solid) + F, _, _ = get_Kinematics(K) + Ke=Kinematics(Electro,Solid) + E = get_Kinematics(Ke) ∂Ψ_∂F(F, E, θ) = TensorValue(ForwardDiff.gradient(F -> Ψ(F, get_array(E), θ), get_array(F))) ∂Ψ_∂E(F, E, θ) = VectorValue(ForwardDiff.gradient(E -> Ψ(get_array(F), E, θ), get_array(E))) @@ -438,8 +452,9 @@ end model = VolumetricEnergy(λ=0.0 ) Ψ, ∂Ψu, ∂Ψuu= model() - F, _, _ = get_Kinematics(model.Kinematic) - + K=Kinematics(Mechano,Solid) + F, _, _ = get_Kinematics(K) + ∂Ψu_(F) =TensorValue(ForwardDiff.gradient(x -> Ψ(x), get_array(F))) ∂Ψuu_(F) =TensorValue(ForwardDiff.hessian(x -> Ψ(x), get_array(F))) @@ -466,8 +481,10 @@ end modelTEM = ThermoElectroMech_Govindjee(modelT, modelID, modelMR, fθ=f, dfdθ=df, gθ=g, dgdθ=dg, β=0.0) Ψ, ∂Ψu, ∂ΨE, ∂Ψθ, ∂ΨFF, ∂ΨEE, ∂2Ψθθ, ∂ΨEF, ∂ΨFθ, ∂ΨEθ, η = modelTEM() - F, _, _ = get_Kinematics(modelMR.Kinematic) - E = get_Kinematics(modelID.Kinematic) + K=Kinematics(Mechano,Solid) + F, _, _ = get_Kinematics(K) + Ke=Kinematics(Electro,Solid) + E = get_Kinematics(Ke) @test Ψ(F(∇u), E(∇φ), θt) == -7.104365408674424 @test norm(∂Ψu(F(∇u), E(∇φ), θt)) == 11.921289845756304 @@ -521,10 +538,11 @@ end ∇u = TensorValue(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0) * 1e-3 ∇φ = VectorValue(1.0, 2.0, 3.0) θt = 3.4 - 1.0 - Kinematic_mec = Kinematics(Mechano) - Kinematic_elec = Kinematics(Electro) - F, _, _ = get_Kinematics(Kinematic_mec) - E = get_Kinematics(Kinematic_elec) + + K=Kinematics(Mechano,Solid) + F, _, _ = get_Kinematics(K) + Ke=Kinematics(Electro,Solid) + E = get_Kinematics(Ke) @test isapprox(Ψ(F(∇u), E(∇φ), θt), 34.24573625846419, atol=1e-12) @test norm(∂ΨF(F(∇u), E(∇φ), θt)) == 12.190784442767743 @@ -550,8 +568,10 @@ end modelID = IdealMagnetic2D(μ=1.2566e-6, χe=0.0) Ψ, ∂Ψu, ∂Ψφ, ∂Ψuu, ∂Ψφu, ∂Ψφφ = modelID() - F, _, _ = get_Kinematics(modelMR.Kinematic) - H0 = get_Kinematics(modelID.Kinematic) + K=Kinematics(Mechano,Solid) + F, _, _ = get_Kinematics(K) + Km=Kinematics(Magneto,Solid) + H0 = get_Kinematics(Km) # ∂Ψφ_(H) =VectorValue(ForwardDiff.gradient(x -> Ψ(get_array( F(∇u)), x,get_array(N) ), get_array(H))) @@ -578,8 +598,10 @@ end modelID = IdealMagnetic(μ=1.2566e-6, χe=0.0) Ψ, ∂Ψu, ∂Ψφ, ∂Ψuu, ∂Ψφu, ∂Ψφφ = modelID() - F, _, _ = get_Kinematics(modelMR.Kinematic) - H0 = get_Kinematics(modelID.Kinematic) + K=Kinematics(Mechano,Solid) + F, _, _ = get_Kinematics(K) + Km=Kinematics(Magneto,Solid) + H0 = get_Kinematics(Km) # ∂Ψu_(F) =TensorValue(ForwardDiff.gradient(x -> Ψ(x, get_array( H0(∇φ)),get_array(N) ), get_array(F))) @@ -614,8 +636,10 @@ end modelID = HardMagnetic(μ=1.2566e-6, αr=40e-3, χe=0.0, χr=8.0) modelmagneto=modelMR+modelID Ψ, ∂Ψu, ∂Ψφ, ∂Ψuu, ∂Ψφu, ∂Ψφφ = modelmagneto() - F, _, _ = get_Kinematics(modelMR.Kinematic) - H0 = get_Kinematics(modelID.Kinematic) + K=Kinematics(Mechano,Solid) + F, _, _ = get_Kinematics(K) + Km=Kinematics(Magneto,Solid) + H0 = get_Kinematics(Km) # ∂Ψu_(F) =TensorValue(ForwardDiff.gradient(x -> Ψ(x, get_array( H0(∇φ)),get_array(N) ), get_array(F))) @@ -646,7 +670,9 @@ end ra=Ref(a) modelID = Magnetic(μ=1.2566e-6, αr=ra ,χe=0.0) Ψ, ∂Ψφ, ∂Ψφφ = modelID() - H0 = get_Kinematics(modelID.Kinematic) + + Km=Kinematics(Magneto,Solid) + H0 = get_Kinematics(Km) N = VectorValue(0.0, 0.0, 1.0) # ∂Ψφ_(H) =VectorValue(ForwardDiff.gradient(x -> Ψ( x,get_array(N) ), get_array(H))) @@ -668,7 +694,9 @@ end ra=Ref(a) modelID = Magnetic(μ=1.2566e-6, αr=ra ,χe=0.0) Ψ, ∂Ψφ, ∂Ψφφ = modelID() - H0 = get_Kinematics(modelID.Kinematic) + + Km=Kinematics(Magneto,Solid) + H0 = get_Kinematics(Km) N = VectorValue(0.0, 0.0) # ∂Ψφ_(H) =VectorValue(ForwardDiff.gradient(x -> Ψ( x,get_array(N) ), get_array(H))) @@ -695,9 +723,10 @@ end modelID = IdealMagnetic2D(μ=1.2566e-6, χe=0.0) modelmagneto=modelMR+modelID Ψ, ∂Ψu, ∂Ψφ, ∂Ψuu, ∂Ψφu, ∂Ψφφ = modelmagneto() - F, _, _ = get_Kinematics(modelMR.Kinematic) - H0 = get_Kinematics(modelID.Kinematic) - + K=Kinematics(Mechano,Solid) + F, _, _ = get_Kinematics(K) + Km=Kinematics(Magneto,Solid) + H0 = get_Kinematics(Km) # ∂Ψu_(F) =TensorValue(ForwardDiff.gradient(x -> Ψ(x, get_array( H0(∇φ)),get_array(N) ), get_array(F))) # ∂Ψφ_(H) =VectorValue(ForwardDiff.gradient(x -> Ψ(get_array( F(∇u)), x,get_array(N) ), get_array(H))) @@ -737,8 +766,10 @@ end modelID = HardMagnetic2D(μ=1.2566e-6, αr=40e-3, χe=0.0, χr=8.0) modelmagneto=modelMR+modelID Ψ, ∂Ψu, ∂Ψφ, ∂Ψuu, ∂Ψφu, ∂Ψφφ = modelmagneto() - F, _, _ = get_Kinematics(modelMR.Kinematic) - H0 = get_Kinematics(modelID.Kinematic) + K=Kinematics(Mechano,Solid) + F, _, _ = get_Kinematics(K) + Km=Kinematics(Magneto,Solid) + H0 = get_Kinematics(Km) # ∂Ψu_(F) =TensorValue(ForwardDiff.gradient(x -> Ψ(x, get_array( H0(∇φ)),get_array(N) ), get_array(F))) @@ -772,7 +803,7 @@ end @testset "ARAP2D" begin model = ARAP2D(μ=μParams[1]) - test_derivatives_2D_(model, rtol=1e-13) + test_derivatives_2D_(model, Kinematics(Mechano,Solid), rtol=1e-13) end @@ -780,7 +811,7 @@ end @testset "ARAP2D_regularized" begin model = ARAP2D_regularized(μ=μParams[1]) - test_derivatives_2D_(model, rtol=1e-13) + test_derivatives_2D_(model, Kinematics(Mechano,Solid),rtol=1e-13) end @@ -792,7 +823,9 @@ end modelreg = HessianRegularization(model) Ψ, ∂Ψu, ∂Ψuu = modelreg() - F, _, J_ = get_Kinematics(modelreg.Kinematic) + K=Kinematics(Mechano,Solid) + F, _, _ = get_Kinematics(K) + # ∂Ψu_(F) =TensorValue(ForwardDiff.gradient(x -> Ψ(x), get_array(F))) # ∂Ψuu_(F) =TensorValue(ForwardDiff.hessian(x -> Ψ(x), get_array(F))) @@ -816,7 +849,9 @@ end modelreg = Hessian∇JRegularization(model) Ψ, ∂Ψu, ∂Ψuu = modelreg() - F, _, J_ = get_Kinematics(modelreg.Kinematic) + K=Kinematics(Mechano,Solid) + F, _, J_ = get_Kinematics(K) + # ∂Ψu_(F) =TensorValue(ForwardDiff.gradient(x -> Ψ(x), get_array(F))) # ∂Ψuu_(F) =TensorValue(ForwardDiff.hessian(x -> Ψ(x), get_array(F))) diff --git a/test/TestConstitutiveModels/ViscousModelsTests.jl b/test/TestConstitutiveModels/ViscousModelsTests.jl index 8ec6263..5c4e578 100644 --- a/test/TestConstitutiveModels/ViscousModelsTests.jl +++ b/test/TestConstitutiveModels/ViscousModelsTests.jl @@ -94,10 +94,6 @@ end struct EmptyElastic <: Elasto - Kinematic::KinematicModel - function EmptyElastic() - new(Kinematics(Elasto)) - end function (::EmptyElastic)(Λ=0.0) Ψ(F) = 0.0 ∂Ψu(F) = TensorValue(zeros(9)...) diff --git a/test/TestWeakForms/WeakForms.jl b/test/TestWeakForms/WeakForms.jl index aa67b2e..d000845 100644 --- a/test/TestWeakForms/WeakForms.jl +++ b/test/TestWeakForms/WeakForms.jl @@ -1,5 +1,4 @@ using Gridap -using HyperFEM: jacobian @testset "Assembly Jacobian ThermoMechanics" begin @@ -35,14 +34,17 @@ using HyperFEM: jacobian θ(x) = x[1] - 1.0 θh = interpolate_everywhere(θ, Vθ) + kt=Kinematics(Thermo, Solid) + km=Kinematics(Mechano,Solid) + k=(km,kt) function jach(uh, θh) - jac((du, dθ), (v, vθ)) = jacobian(modelTM, (uh, θh), (du, dθ), (v, vθ), dΩ) + jac((du, dθ), (v, vθ)) = jacobian(modelTM, k,(uh, θh), (du, dθ), (v, vθ), dΩ) end function jac_mech(uh, θh) - jac(du, v) = jacobian(modelTM, Mechano, (uh, θh), du, v, dΩ) + jac(du, v) = jacobian(modelTM, Mechano, k,(uh, θh), du, v, dΩ) end function jac_termoh(uh, θh) - jac(dθ, vθ) = jacobian(modelTM, Thermo, dθ, vθ, dΩ) + jac(dθ, vθ) = jacobian(modelTM, Thermo,k, dθ, vθ, dΩ) end jac_ = assemble_matrix(jach(uh, θh), V, V) @@ -95,14 +97,18 @@ end φ(x) = x[1] φh = interpolate_everywhere(φ, Vφ) + + ke=Kinematics(Electro, Solid) + km=Kinematics(Mechano,Solid) + k=(km,ke) function jach(uh, φh) - jac((du, dφ), (v, vφ)) = jacobian(modelelectro, (uh, φh), (du, dφ), (v, vφ), dΩ) + jac((du, dφ), (v, vφ)) = jacobian(modelelectro,k, (uh, φh), (du, dφ), (v, vφ), dΩ) end function jac_mech(uh, φh) - jac(du, v) = jacobian(modelelectro, Mechano, (uh, φh), du, v, dΩ) + jac(du, v) = jacobian(modelelectro, Mechano, k,(uh, φh), du, v, dΩ) end function jac_elech(uh, φh) - jac(dφ, vφ) = jacobian(modelelectro, Electro, (uh, φh), dφ, vφ, dΩ) + jac(dφ, vφ) = jacobian(modelelectro, Electro, k,(uh, φh), dφ, vφ, dΩ) end jac_ = assemble_matrix(jach(uh, φh), V, V) diff --git a/test/data/ViscoElasticSimulation.jl b/test/data/ViscoElasticSimulation.jl index 66fdd28..6b2d340 100644 --- a/test/data/ViscoElasticSimulation.jl +++ b/test/data/ViscoElasticSimulation.jl @@ -61,8 +61,9 @@ function visco_elastic_simulation(;t_end=15, write_vtk=true, verbose=true) unh = FEFunction(Uun, zero_free_values(Uun)) state_vars = initializeStateVariables(cons_model, dΩ) - res(Λ) = (u,v)->residual(cons_model, u, v, dΩ, t_end * Λ, Δt, unh, state_vars) - jac(Λ) = (u,du,v)->jacobian(cons_model, u, du, v, dΩ, t_end * Λ, Δt, unh, state_vars) + k=Kinematics(Mechano,Solid) + res(Λ) = (u,v)->residual(cons_model, k, u, v, dΩ, t_end * Λ, Δt, unh, state_vars) + jac(Λ) = (u,du,v)->jacobian(cons_model, k, u, du, v, dΩ, t_end * Λ, Δt, unh, state_vars) ls = LUSolver() nls = NewtonSolver(ls; maxiter=20, atol=1.e-6, rtol=1.e-6, verbose=verbose) @@ -70,12 +71,14 @@ function visco_elastic_simulation(;t_end=15, write_vtk=true, verbose=true) λx = Float64[] σΓ = Float64[] + F,_,_ = get_Kinematics(k) + function driverpost(post) - σh11, _... = Cauchy(cons_model, uh, unh, state_vars, Ω, dΩ, 0.0, Δt) + σh11, _... = Cauchy(cons_model, Kinematics(Mechano,Solid),uh, unh, state_vars, Ω, dΩ, 0.0, Δt) σΓ1 = sum(∫(σh11)dΓ1) / sum(∫(1.0)dΓ1) push!(σΓ, σΓ1) push!(λx, 1.0 + component_LInf(uh, :x, Ω) / long) - updateStateVariables!(state_vars, cons_model, Δt, uh, unh) + updateStateVariables!(state_vars, cons_model, Δt, F∘(∇(uh)'), F∘(∇(unh)')) end post_model = PostProcessor(comp_model, driverpost; is_vtk=false, filepath="")