Refactor Electro and ThermoElectro models

src/Exports.jl

@@ -24,7 +24,6 @@ end
 @publish TensorAlgebra Tensorize
 
 
-@publish PhysicalModels DerivativeStrategy
 @publish PhysicalModels LinearElasticity3D
 @publish PhysicalModels LinearElasticity2D
 @publish PhysicalModels Yeoh3D
@@ -49,6 +48,7 @@ end
 @publish PhysicalModels TransverseIsotropy2D
 @publish PhysicalModels ThermalModel
 @publish PhysicalModels ThermalVolumetric
+@publish PhysicalModels ThermalDeviatoric
 @publish PhysicalModels IdealDielectric
 @publish PhysicalModels Magnetic
 @publish PhysicalModels IdealMagnetic
@@ -58,6 +58,7 @@ end
 @publish PhysicalModels ElectroMechModel
 @publish PhysicalModels ThermoElectroMechModel
 @publish PhysicalModels ThermoMechModel
+@publish PhysicalModels ThermoElectroModel
 @publish PhysicalModels ThermoMech_Bonet
 @publish PhysicalModels ThermoMech_EntropicPolyconvex
 @publish PhysicalModels FlexoElectroModel
@@ -94,10 +95,11 @@ end
 @publish PhysicalModels getIsoInvariants
 
 @publish PhysicalModels ThermalLaw
+@publish PhysicalModels ConstantEnergyLaw
+@publish PhysicalModels ConstantCvLaw
 @publish PhysicalModels EntropicElasticityLaw
 @publish PhysicalModels NonlinearMeltingLaw
 @publish PhysicalModels NonlinearSofteningLaw
-@publish PhysicalModels TrigonometricLaw
 @publish PhysicalModels PolynomialLaw
 
 @publish PhysicalModels SecondPiola

src/PhysicalModels/ElectricalModels.jl

@@ -9,3 +9,38 @@ struct IdealDielectric <: Electro
     new(ε)
   end
 end
+
+function (obj::IdealDielectric)()
+  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)
+  Ψem(F, E) = (-obj.ε / (2.0 * J(F))) * HEHE(F, E)
+
+  # First Derivatives #
+  ∂Ψem_∂H(F, E) = (-obj.ε / (J(F))) * (HE(F, E) ⊗ E)
+  ∂Ψem_∂J(F, E) = (+obj.ε / (2.0 * J(F)^2.0)) * HEHE(F, E)
+  ∂Ψem_∂E(F, E) = (-obj.ε / (J(F))) * (H(F)' * HE(F, E))
+  ∂Ψem∂F(F, E) = ∂Ψem_∂H(F, E) × F + ∂Ψem_∂J(F, E) * H(F)
+  ∂Ψem∂E(F, E) = ∂Ψem_∂E(F, E)
+
+  # Second Derivatives #
+  ∂Ψem_HH(F, E) = (-obj.ε / (J(F))) * (I3 ⊗₁₃²⁴ (E ⊗ E))
+  ∂Ψem_HJ(F, E) = (+obj.ε / (J(F))^2.0) * (HE(F, E) ⊗ E)
+  ∂Ψem_JJ(F, E) = (-obj.ε / (J(F))^3.0) * HEHE(F, E)
+  ∂Ψem∂FF(F, E) = (F × (∂Ψem_HH(F, E) × F)) +
+                  H(F) ⊗₁₂³⁴ (∂Ψem_HJ(F, E) × F) +
+                  (∂Ψem_HJ(F, E) × F) ⊗₁₂³⁴ H(F) +
+                  ∂Ψem_JJ(F, E) * (H(F) ⊗₁₂³⁴ H(F)) +
+                  ×ᵢ⁴(∂Ψem_∂H(F, E) + ∂Ψem_∂J(F, E) * F)
+
+  ∂Ψem_EH(F, E) = (-obj.ε / (J(F))) * ((I3 ⊗₁₃² HE(F, E)) + (H(F)' ⊗₁₂³ E))
+  ∂Ψem_EJ(F, E) = (+obj.ε / (J(F))^2.0) * (H(F)' * HE(F, E))
+  ∂Ψem∂EF(F, E) = (∂Ψem_EH(F, E) × F) + (∂Ψem_EJ(F, E) ⊗₁²³ H(F))
+
+  ∂Ψem∂EE(F, E) = (-obj.ε / (J(F))) * (H(F)' * H(F))
+
+  return (Ψem, ∂Ψem∂F, ∂Ψem∂E, ∂Ψem∂FF, ∂Ψem∂EF, ∂Ψem∂EE)
+end

src/PhysicalModels/ElectroMechanicalModels.jl

@@ -10,60 +10,76 @@ struct ElectroMechModel{E<:Electro,M<:Mechano} <: ElectroMechano{E,M}
   function ElectroMechModel(; electro::E, mechano::M) where {E<:Electro,M<:Mechano}
     new{E,M}(electro, mechano)
   end
+end
 
-  function (obj::ElectroMechModel{<:Electro,<:IsoElastic})(Λ::Float64=1.0)
-    Ψm, ∂Ψm_u, ∂Ψm_uu = obj.mechano()
-    Ψem, ∂Ψem_u, ∂Ψem_φ, ∂Ψem_uu, ∂Ψem_φu, ∂Ψem_φφ = _getCoupling(obj.electro, obj.mechano)
-    Ψ(F, E) = Ψm(F) + Ψem(F, E)
-    ∂Ψu(F, E) = ∂Ψm_u(F) + ∂Ψem_u(F, E)
-    ∂Ψφ(F, E) = ∂Ψem_φ(F, E)
-    ∂Ψuu(F, E) = ∂Ψm_uu(F) + ∂Ψem_uu(F, E)
-    ∂Ψφu(F, E) = ∂Ψem_φu(F, E)
-    ∂Ψφφ(F, E) = ∂Ψem_φφ(F, E)
-    return (Ψ, ∂Ψu, ∂Ψφ, ∂Ψuu, ∂Ψφu, ∂Ψφφ)
-  end
-  function (obj::ElectroMechModel{<:Electro,<:AnisoElastic})(Λ::Float64=1.0)
-    Ψm, ∂Ψm_u, ∂Ψm_uu = obj.mechano()
-    Ψem, ∂Ψem_u, ∂Ψem_φ, ∂Ψem_uu, ∂Ψem_φu, ∂Ψem_φφ = _getCoupling(obj.electro, obj.mechano)
-    Ψ(F, E, N) = Ψm(F, N) + Ψem(F, E)
-    ∂Ψu(F, E, N) = ∂Ψm_u(F, N) + ∂Ψem_u(F, E)
-    ∂Ψφ(F, E, N) = ∂Ψem_φ(F, E)
-    ∂Ψuu(F, E, N) = ∂Ψm_uu(F, N) + ∂Ψem_uu(F, E)
-    ∂Ψφu(F, E, N) = ∂Ψem_φu(F, E)
-    ∂Ψφφ(F, E, N) = ∂Ψem_φφ(F, E)
-    return (Ψ, ∂Ψu, ∂Ψφ, ∂Ψuu, ∂Ψφu, ∂Ψφφ)
-  end
-  function (obj::ElectroMechModel{<:Electro,<:ViscoElastic{<:IsoElastic}})(Λ::Float64=1.0)
-    Ψm, ∂Ψm_u, ∂Ψm_uu = obj.mechano()
-    Ψem, ∂Ψem_u, ∂Ψem_φ, ∂Ψem_uu, ∂Ψem_φu, ∂Ψem_φφ = _getCoupling(obj.electro, obj.mechano)
-    Ψ(F, E, Fn, A...) = Ψm(F, Fn, A...) + Ψem(F, E)
-    ∂Ψu(F, E, Fn, A...) = ∂Ψm_u(F, Fn, A...) + ∂Ψem_u(F, E)
-    ∂Ψφ(F, E, Fn, A...) = ∂Ψem_φ(F, E)
-    ∂Ψuu(F, E, Fn, A...) = ∂Ψm_uu(F, Fn, A...) + ∂Ψem_uu(F, E)
-    ∂Ψφu(F, E, Fn, A...) = ∂Ψem_φu(F, E)
-    ∂Ψφφ(F, E, Fn, A...) = ∂Ψem_φφ(F, E)
-    return (Ψ, ∂Ψu, ∂Ψφ, ∂Ψuu, ∂Ψφu, ∂Ψφφ)
-  end
-  function (obj::ElectroMechModel{<:Electro,<:ViscoElastic{<:AnisoElastic}})(Λ::Float64=1.0)
-    Ψm, ∂Ψm_u, ∂Ψm_uu = obj.mechano()
-    Ψem, ∂Ψem_u, ∂Ψem_φ, ∂Ψem_uu, ∂Ψem_φu, ∂Ψem_φφ = _getCoupling(obj.electro, obj.mechano)
-    Ψ(F, E, n, Fn, A...) = Ψm(F, n, Fn, A...) + Ψem(F, E)
-    ∂Ψu(F, E, n, Fn, A...) = ∂Ψm_u(F, n, Fn, A...) + ∂Ψem_u(F, E)
-    ∂Ψφ(F, E, n, Fn, A...) = ∂Ψem_φ(F, E)
-    ∂Ψuu(F, E, n, Fn, A...) = ∂Ψm_uu(F, n, Fn, A...) + ∂Ψem_uu(F, E)
-    ∂Ψφu(F, E, n, Fn, A...) = ∂Ψem_φu(F, E)
-    ∂Ψφφ(F, E, n, Fn, A...) = ∂Ψem_φφ(F, E)
-    return (Ψ, ∂Ψu, ∂Ψφ, ∂Ψuu, ∂Ψφu, ∂Ψφφ)
-  end
+function (+)(Model1::Electro, Model2::Mechano)
+  ElectroMechModel(Model1, Model2)
+end
+
+function (+)(Model1::Mechano, Model2::Electro)
+  ElectroMechModel(Model2, Model1)
+end
+
+function (obj::ElectroMechModel{<:Electro,<:IsoElastic})(Λ::Float64=1.0)
+  Ψm, ∂Ψm∂F, ∂Ψm∂FF = obj.mechano()
+  Ψem, ∂Ψem∂F, ∂Ψem∂E, ∂Ψem∂FF, ∂Ψem∂EF, ∂Ψem∂EE = obj.electro()
+  Ψ(F, E)     = Ψm(F)     + Ψem(F, E)
+  ∂Ψ∂F(F, E)  = ∂Ψm∂F(F)  + ∂Ψem∂F(F, E)
+  ∂Ψ∂E(F, E)  =             ∂Ψem∂E(F, E)
+  ∂Ψ∂FF(F, E) = ∂Ψm∂FF(F) + ∂Ψem∂FF(F, E)
+  ∂Ψ∂EF(F, E) =             ∂Ψem∂EF(F, E)
+  ∂Ψ∂EE(F, E) =             ∂Ψem∂EE(F, E)
+  return (Ψ, ∂Ψ∂F, ∂Ψ∂E, ∂Ψ∂FF, ∂Ψ∂EF, ∂Ψ∂EE)
+end
+
+function (obj::ElectroMechModel{<:Electro,<:AnisoElastic})(Λ::Float64=1.0)
+  Ψm, ∂Ψm∂F, ∂Ψm∂FF = obj.mechano()
+  Ψem, ∂Ψem∂F, ∂Ψem∂E, ∂Ψem∂FF, ∂Ψem∂EF, ∂Ψem∂EE = obj.electro()
+  Ψ(F, E, N)     = Ψm(F, N)     + Ψem(F, E)
+  ∂Ψ∂F(F, E, N)  = ∂Ψm∂F(F, N)  + ∂Ψem∂F(F, E)
+  ∂Ψ∂E(F, E, N)  =                ∂Ψem∂E(F, E)
+  ∂Ψ∂FF(F, E, N) = ∂Ψm∂FF(F, N) + ∂Ψem∂FF(F, E)
+  ∂Ψ∂EF(F, E, N) =                ∂Ψem∂EF(F, E)
+  ∂Ψ∂EE(F, E, N) =                ∂Ψem∂EE(F, E)
+  return (Ψ, ∂Ψ∂F, ∂Ψ∂E, ∂Ψ∂FF, ∂Ψ∂EF, ∂Ψ∂EE)
+end
+
+function (obj::ElectroMechModel{<:Electro,<:ViscoElastic{<:IsoElastic}})(Λ::Float64=1.0)
+  Ψm, ∂Ψm∂F, ∂Ψm∂FF = obj.mechano()
+  Ψem, ∂Ψem∂F, ∂Ψem∂E, ∂Ψem∂FF, ∂Ψem∂EF, ∂Ψem∂EE = obj.electro()
+  Ψ(F, E, Fn, A...)     = Ψm(F, Fn, A...)     + Ψem(F, E)
+  ∂Ψ∂F(F, E, Fn, A...)  = ∂Ψm∂F(F, Fn, A...)  + ∂Ψem∂F(F, E)
+  ∂Ψ∂E(F, E, Fn, A...)  =                       ∂Ψem∂E(F, E)
+  ∂Ψ∂FF(F, E, Fn, A...) = ∂Ψm∂FF(F, Fn, A...) + ∂Ψem∂FF(F, E)
+  ∂Ψ∂EF(F, E, Fn, A...) =                       ∂Ψem∂EF(F, E)
+  ∂Ψ∂EE(F, E, Fn, A...) =                       ∂Ψem∂EE(F, E)
+  return (Ψ, ∂Ψ∂F, ∂Ψ∂E, ∂Ψ∂FF, ∂Ψ∂EF, ∂Ψ∂EE)
+end
+
+function (obj::ElectroMechModel{<:Electro,<:ViscoElastic{<:AnisoElastic}})(Λ::Float64=1.0)
+  Ψm, ∂Ψm∂F, ∂Ψm∂FF = obj.mechano()
+  Ψem, ∂Ψem∂F, ∂Ψem∂E, ∂Ψem∂FF, ∂Ψem∂EF, ∂Ψem∂EE = obj.electro()
+  Ψ(F, E, n, Fn, A...)     = Ψm(F, n, Fn, A...)     + Ψem(F, E)
+  ∂Ψ∂F(F, E, n, Fn, A...)  = ∂Ψm∂F(F, n, Fn, A...)  + ∂Ψem∂F(F, E)
+  ∂Ψ∂E(F, E, n, Fn, A...)  =                          ∂Ψem∂E(F, E)
+  ∂Ψ∂FF(F, E, n, Fn, A...) = ∂Ψm∂FF(F, n, Fn, A...) + ∂Ψem∂FF(F, E)
+  ∂Ψ∂EF(F, E, n, Fn, A...) =                          ∂Ψem∂EF(F, E)
+  ∂Ψ∂EE(F, E, n, Fn, A...) =                          ∂Ψem∂EE(F, E)
+  return (Ψ, ∂Ψ∂F, ∂Ψ∂E, ∂Ψ∂FF, ∂Ψ∂EF, ∂Ψ∂EE)
 end
 
 function update_time_step!(obj::ElectroMechModel, Δt::Float64)
   update_time_step!(obj.electro, Δt)
   update_time_step!(obj.mechano, Δt)
 end
 
+function Gridap.CellData.CellState(obj::ElectroMechModel, args...)
+  CellState(obj.mechano, args...)
+end
+
 function initialize_state(obj::ElectroMechModel, points::Measure)
-  initialize_state(obj.mechano, points)
+  @warn "The function 'initialize_state' is deprecated, use 'CellState' instead."
+  CellState(obj.mechano, points)
 end
 
 function update_state!(obj::ElectroMechModel, state, F, E, args...)
@@ -75,48 +91,6 @@ function Dissipation(obj::ElectroMechModel)
   D(F, E, X...) = Dvis(F, X...)
 end
 
-
-function _getCoupling(elec::Electro, mec::Mechano, Λ::Float64=0.0)
-  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)
-  Ψem(F, E) = (-elec.ε / (2.0 * J(F))) * HEHE(F, E)
-  # First Derivatives #
-  ∂Ψem_∂H(F, E) = (-elec.ε / (J(F))) * (HE(F, E) ⊗ E)
-  ∂Ψem_∂J(F, E) = (+elec.ε / (2.0 * J(F)^2.0)) * HEHE(F, E)
-  ∂Ψem_∂E(F, E) = (-elec.ε / (J(F))) * (H(F)' * HE(F, E))
-  ∂Ψem_u(F, E) = ∂Ψem_∂H(F, E) × F + ∂Ψem_∂J(F, E) * H(F)
-  ∂Ψem_φ(F, E) = ∂Ψem_∂E(F, E)
-
-  # Second Derivatives #
-  ∂Ψem_HH(F, E) = (-elec.ε / (J(F))) * (I3 ⊗₁₃²⁴ (E ⊗ E))
-  ∂Ψem_HJ(F, E) = (+elec.ε / (J(F))^2.0) * (HE(F, E) ⊗ E)
-  ∂Ψem_JJ(F, E) = (-elec.ε / (J(F))^3.0) * HEHE(F, E)
-  ∂Ψem_uu(F, E) = (F × (∂Ψem_HH(F, E) × F)) +
-                  H(F) ⊗₁₂³⁴ (∂Ψem_HJ(F, E) × F) +
-                  (∂Ψem_HJ(F, E) × F) ⊗₁₂³⁴ H(F) +
-                  ∂Ψem_JJ(F, E) * (H(F) ⊗₁₂³⁴ H(F)) +
-                  ×ᵢ⁴(∂Ψem_∂H(F, E) + ∂Ψem_∂J(F, E) * F)
-
-  ∂Ψem_EH(F, E) = (-elec.ε / (J(F))) * ((I3 ⊗₁₃² HE(F, E)) + (H(F)' ⊗₁₂³ E))
-  ∂Ψem_EJ(F, E) = (+elec.ε / (J(F))^2.0) * (H(F)' * HE(F, E))
-  ∂Ψem_φu(F, E) = (∂Ψem_EH(F, E) × F) + (∂Ψem_EJ(F, E) ⊗₁²³ H(F))
-
-  ∂Ψem_φφ(F, E) = (-elec.ε / (J(F))) * (H(F)' * H(F))
-
-  return (Ψem, ∂Ψem_u, ∂Ψem_φ, ∂Ψem_uu, ∂Ψem_φu, ∂Ψem_φφ)
-end
-
-
-function (+)(Model1::Electro, Model2::Mechano)
-  ElectroMechModel(Model1, Model2)
-end
-function (+)(Model1::Mechano, Model2::Electro)
-  ElectroMechModel(Model2, Model1)
-end
-
 struct FlexoElectroModel{EM<:ElectroMechano} <: FlexoElectro{EM}
   electromechano::EM
   κ::Float64
@@ -141,7 +115,7 @@ struct FlexoElectroModel{EM<:ElectroMechano} <: FlexoElectro{EM}
     f3(δϕ) = δϕ ⊗₁² e₃
     Φ(ϕ₁, ϕ₂, ϕ₃) = (f1 ∘ (ϕ₁) + f2 ∘ (ϕ₂) + f3 ∘ (ϕ₃))
 
-    Ψ, ∂Ψu, ∂Ψφ, ∂Ψuu, ∂Ψφu, ∂Ψφφ = obj.electromechano(Λ)
-    return Ψ, ∂Ψu, ∂Ψφ, ∂Ψuu, ∂Ψφu, ∂Ψφφ, Φ
+    Ψ, ∂Ψ∂F, ∂Ψ∂E, ∂Ψ∂FF, ∂Ψ∂EF, ∂Ψ∂EE = obj.electromechano(Λ)
+    return Ψ, ∂Ψ∂F, ∂Ψ∂E, ∂Ψ∂FF, ∂Ψ∂EF, ∂Ψ∂EE, Φ
   end
 end

src/PhysicalModels/PhysicalModels.jl

@@ -3,7 +3,6 @@ module PhysicalModels
 using Gridap
 using Gridap.CellData
 using Gridap.Helpers
-using UnPack
 using ForwardDiff
 using LinearAlgebra
 using StaticArrays
@@ -47,9 +46,11 @@ export HardMagnetic
 export HardMagnetic2D
 export ThermalModel
 export ThermalVolumetric
+export ThermalDeviatoric
 export ElectroMechModel
 export ThermoElectroMechModel
 export ThermoMechModel
+export ThermoElectroModel
 export ThermoMech_Bonet
 export ThermoMech_EntropicPolyconvex
 export FlexoElectroModel
@@ -81,8 +82,6 @@ export EnergyInterpolationScheme
 export SecondPiola
 export Dissipation
 
-export DerivativeStrategy
-
 export initialize_state
 export update_time_step!
 
@@ -97,8 +96,6 @@ export getIsoInvariants
 export HessianRegularization
 export Hessian∇JRegularization
 
-struct DerivativeStrategy{Kind} end
-
 abstract type PhysicalModel end
 abstract type Mechano <: PhysicalModel end
 abstract type Electro <: PhysicalModel end
@@ -119,7 +116,7 @@ abstract type MultiPhysicalModel <: PhysicalModel end
 abstract type ElectroMechano{E,M} <: MultiPhysicalModel end
 abstract type ThermoElectroMechano{T,E,M} <: MultiPhysicalModel end
 abstract type ThermoMechano{T,M} <: MultiPhysicalModel end
-abstract type ThermoElectro{E,M} <: MultiPhysicalModel end
+abstract type ThermoElectro{E} <: MultiPhysicalModel end
 abstract type FlexoElectro{EM} <: MultiPhysicalModel end
 abstract type MagnetoMechano{G,M} <: MultiPhysicalModel end
 
@@ -137,6 +134,8 @@ include("ThermalModels.jl")
 
 include("ThermoMechanicalModels.jl")
 
+include("ThermoElectroModels.jl")
+
 include("ElectroMechanicalModels.jl")
 
 include("MagnetoMechanicalModels.jl")
@@ -155,13 +154,17 @@ Base.broadcastable(m::PhysicalModel) = Ref(m) # Allows to use the @. syntax for
 """
 Initialize the state variables for the given constitutive model and discretization.
 """
+function Gridap.CellData.CellState(::PhysicalModel, args...)
+  return nothing
+end
+
 function initialize_state(::PhysicalModel, points::Measure)
   return nothing
 end
 
 
 """
-Update the state variables. The state variables must be initialized using the function 'initialize_state'.
+Update the state variables. The state variables must be initialized using the function 'CellState' with the constitutive model.
 """
 function update_state!(::PhysicalModel, vars...)
 end