Refactors kinematic descriptions in physical models

Introduces a more explicit kinematic description and associates
it directly with the physics of the solid material, thus decoupling
the kinematic model from the specific electrical or magnetic model.

This change simplifies the model definitions and usage in weak forms
and post-processors, improving overall code clarity and maintainability.

Author: jmartfrut

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()

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)

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

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

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)

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)

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)

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 *= Λ
 
   # #-------------------------------------------------------------------------------------