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="")