[NOT TO MERGE] Refactor Electro and ThermoElectro models

src/Exports.jl

@@ -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
@@ -97,7 +98,6 @@ end
 @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::Electro)()
+  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,51 +10,62 @@ 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)
@@ -75,48 +86,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 +110,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

@@ -50,6 +50,7 @@ export ThermalVolumetric
 export ElectroMechModel
 export ThermoElectroMechModel
 export ThermoMechModel
+export ThermoElectroModel
 export ThermoMech_Bonet
 export ThermoMech_EntropicPolyconvex
 export FlexoElectroModel
@@ -119,7 +120,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 +138,8 @@ include("ThermalModels.jl")
 
 include("ThermoMechanicalModels.jl")
 
+include("ThermoElectroModels.jl")
+
 include("ElectroMechanicalModels.jl")
 
 include("MagnetoMechanicalModels.jl")

src/PhysicalModels/ThermalModels.jl

@@ -72,27 +72,12 @@ function (law::NonlinearSofteningLaw)()
   return (f, ∂f, ∂∂f)
 end
 
-struct TrigonometricLaw <: ThermalLaw
-  θr::Float64
-  θM::Float64
-end
-
-function (law::TrigonometricLaw)()
-  @unpack θr, θM = law  
-  g(θ) = θ/θr * sin(2π*θ/θM)
-  G(θ) = 1/2/π * θM/θr * (1 - cos(2π*θ/θM))
-  H(θ) = 1/2/π * θM/θr * (θ - θM/2/π * sin(2π*θ/θM))
-  f(θ) = (H(θr) - H(θ)) / (H(θM) - H(θr)) + 1.0
-  ∂f(θ) = -G(θ) / (H(θM) - H(θr))
-  ∂∂f(θ) = -g(θ) / θ / (H(θM) - H(θr))
-  return (f, ∂f, ∂∂f)
-end
-
 struct PolynomialLaw <: ThermalLaw
   θr::Float64
   a::Float64
   b::Float64
   c::Float64
+  PolynomialLaw(; θr, a, b, c) = new(θr, a, b, c)
 end
 
 function (law::PolynomialLaw)()

src/PhysicalModels/ThermoElectroMechanicalModels.jl

@@ -31,7 +31,7 @@ struct ThermoElectroMechModel{T<:Thermo,E<:Electro,M<:Mechano} <: ThermoElectroM
   function (obj::ThermoElectroMechModel)(Λ::Float64=1.0)
     Ψt, ∂Ψt_θ, ∂Ψt_θθ = obj.thermo(Λ)
     Ψm, ∂Ψm_u, ∂Ψm_uu = obj.mechano(Λ)
-    Ψem, ∂Ψem_u, ∂Ψem_φ, ∂Ψem_uu, ∂Ψem_φu, ∂Ψem_φφ = _getCoupling(obj.electro, obj.mechano, Λ)
+    Ψem, ∂Ψem_u, ∂Ψem_φ, ∂Ψem_uu, ∂Ψem_φu, ∂Ψem_φφ = obj.electro()
     Ψtm, ∂Ψtm_u, ∂Ψtm_θ, ∂Ψtm_uu, ∂Ψtm_uθ, ∂Ψtm_θθ = _getCoupling(obj.thermo, obj.mechano, Λ)
     f(δθ) = (obj.fθ(δθ)::Float64)
     df(δθ) = (obj.dfdθ(δθ)::Float64)
@@ -73,7 +73,7 @@ struct ThermoElectroMech_Govindjee{T<:Thermo,E<:Electro,M<:Mechano} <: ThermoEle
 
   function (obj::ThermoElectroMech_Govindjee)(Λ::Float64=1.0)
     Ψm, _, _ = obj.mechano(Λ)
-    Ψem, _, _, _, _, _ = _getCoupling(obj.electro, obj.mechano, Λ)
+    Ψem, _, _, _, _, _ = obj.electro()
     f(δθ) = obj.fθ(δθ)
     df(δθ) = obj.dfdθ(δθ)
     g(δθ) = obj.gθ(δθ)
@@ -134,10 +134,14 @@ function ThermoElectroMech_Bonet(thermo::ThermalVolumetric, electro::E, mechano:
   ThermoElectroMech_Bonet{E,M}(thermo,electro,mechano,el,vis,elec)
 end
 
+function ThermoElectroMech_Bonet(thermo::ThermalVolumetric, electro::ThermoElectro{E}, mechano::M; el::ThermalLaw, vis::ThermalLaw) where {E<:Electro,M<:ViscoElastic}
+  ThermoElectroMech_Bonet{E,M}(thermo,electro.electro,mechano,el,vis,electro.law)
+end
+
 function (obj::ThermoElectroMech_Bonet{<:Electro,<:Elasto})()
   Ψt, ∂Ψt∂F, ∂Ψt∂θ, ∂∂Ψt∂FF, ∂∂Ψt∂θθ, ∂∂Ψt∂Fθ = obj.thermo()
   Ψm, ∂Ψm∂F, ∂∂Ψm∂FF = obj.mechano()
-  Ψem, ∂Ψem∂F, ∂Ψem∂E, ∂Ψem∂FF, ∂Ψem∂EF, ∂∂Ψem∂EE = _getCoupling(obj.electro, obj.mechano)
+  Ψem, ∂Ψem∂F, ∂Ψem∂E, ∂Ψem∂FF, ∂Ψem∂EF, ∂∂Ψem∂EE = obj.electro()
   fe, dfe, ddfe = obj.lawel()
   felec, dfelec, ddfelec = obj.lawelec()
 
@@ -159,7 +163,7 @@ function (obj::ThermoElectroMech_Bonet{<:Electro,<:ViscoElastic})()
   Ψt, ∂Ψt∂F, ∂Ψt∂θ, ∂∂Ψt∂FF, ∂∂Ψt∂θθ, ∂∂Ψt∂Fθ = obj.thermo()
   Ψe, ∂Ψe∂F, ∂∂Ψe∂FF = obj.mechano.longterm()
   Ψv, ∂Ψv∂F, ∂∂Ψv∂FF = obj.mechano.branches()
-  Ψem, ∂Ψem∂F, ∂Ψem∂E, ∂Ψem∂FF, ∂Ψem∂EF, ∂∂Ψem∂EE = _getCoupling(obj.electro, obj.mechano)
+  Ψem, ∂Ψem∂F, ∂Ψem∂E, ∂Ψem∂FF, ∂Ψem∂EF, ∂∂Ψem∂EE = obj.electro()
   fe, dfe, ddfe = obj.lawel()
   fv, dfv, ddfv = obj.lawvis()
   felec, dfelec, ddfelec = obj.lawelec()

src/PhysicalModels/ThermoElectroModels.jl

@@ -0,0 +1,27 @@
+
+struct ThermoElectroModel{E<:Electro} <: ThermoElectro{E}
+  electro::E
+  law::ThermalLaw
+
+  function ThermoElectroModel(electro::E, law::ThermalLaw) where {E <: Electro}
+    new{E}(electro, law)
+  end
+end
+
+function (obj::ThermoElectroModel)()
+  Ψem, ∂Ψem∂F, ∂Ψem∂E, ∂∂Ψem∂FF, ∂∂Ψem∂EF, ∂∂Ψem∂EE = obj.electro()
+  f, df, ddf = obj.law()
+
+  Ψ(F, E, θ)       =  f(θ)*Ψem(F,E)
+  ∂Ψ∂F(F, E, θ)    =  f(θ)*∂Ψem∂F(F,E)
+  ∂Ψ∂E(F, E, θ)    =  f(θ)*∂Ψem∂E(F,E)
+  ∂Ψ∂θ(F, E, θ)    =  df(θ)*Ψem(F,E)
+  ∂∂Ψ∂FF(F, E, θ)  =  f(θ)*∂∂Ψem∂FF(F,E)
+  ∂∂Ψ∂EE(F, E, θ)  =  f(θ)*∂∂Ψem∂EE(F,E)
+  ∂∂Ψ∂θθ(F, E, θ)  =  ddf(θ)*Ψem(F,E)
+  ∂∂Ψ∂EF(F, E, θ)  =  f(θ)*∂∂Ψem∂EF(F,E)
+  ∂∂Ψ∂Fθ(F, E, θ)  =  df(θ)*∂Ψem∂F(F,E)
+  ∂∂Ψ∂Eθ(F, E, θ)  =  df(θ)*∂Ψem∂E(F,E)
+
+  return (Ψ, ∂Ψ∂F, ∂Ψ∂E, ∂Ψ∂θ, ∂∂Ψ∂FF, ∂∂Ψ∂EE, ∂∂Ψ∂θθ, ∂∂Ψ∂EF, ∂∂Ψ∂Fθ, ∂∂Ψ∂Eθ)
+end

test/TestConstitutiveModels/ThermalLawsTests.jl

@@ -29,8 +29,8 @@ end
   end
 end
 
-@testset "TrigonometricLaw" begin
-  law = TrigonometricLaw(273.15, 400.0)
+@testset "PolynomialLaw" begin
+  law = PolynomialLaw(θr=273.15, a=1.1, b=2.2, c=3.3)
   f, df, ddf = law()
   for θ ∈ 200.0:50:400
     @test isapprox(df(θ), ForwardDiff.derivative(f, θ), rtol=1e-10)