Merge pull request #50 from MultiSimOLab/membrane

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