diff --git a/src/Exports.jl b/src/Exports.jl
index 9e7d964..68d15f3 100644
--- a/src/Exports.jl
+++ b/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
diff --git a/src/PhysicalModels/ElectricalModels.jl b/src/PhysicalModels/ElectricalModels.jl
index c4b5c20..cc203de 100644
--- a/src/PhysicalModels/ElectricalModels.jl
+++ b/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
diff --git a/src/PhysicalModels/ElectroMechanicalModels.jl b/src/PhysicalModels/ElectroMechanicalModels.jl
index 930df47..d773c6b 100644
--- a/src/PhysicalModels/ElectroMechanicalModels.jl
+++ b/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)
@@ -62,8 +73,13 @@ function update_time_step!(obj::ElectroMechModel, Δt::Float64)
   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
diff --git a/src/PhysicalModels/PhysicalModels.jl b/src/PhysicalModels/PhysicalModels.jl
index b45ae7e..d23847e 100644
--- a/src/PhysicalModels/PhysicalModels.jl
+++ b/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
diff --git a/src/PhysicalModels/ThermalModels.jl b/src/PhysicalModels/ThermalModels.jl
index 0a005ac..03cd8af 100644
--- a/src/PhysicalModels/ThermalModels.jl
+++ b/src/PhysicalModels/ThermalModels.jl
@@ -25,6 +25,29 @@ end
 # Thermal laws
 # ===================
 
+struct ConstantEnergyLaw <: ThermalLaw end
+
+function (law::ConstantEnergyLaw)()
+  f(θ) = 1.0
+  ∂f(θ) = 0.0
+  ∂∂f(θ) = 0.0
+  return (f, ∂f, ∂∂f)
+end
+
+struct ConstantCvLaw <: ThermalLaw
+  θr::Float64
+  ConstantCvLaw(θr) = new(θr)
+  ConstantCvLaw(; θr) = new(θr)
+end
+
+function (law::ConstantCvLaw)()
+  θr = law.θr
+  f(θ) = (θ-θr) -θ*log(θ/θr)
+  ∂f(θ) = -log(θ/θr)
+  ∂∂f(θ) = -1/θ
+  return (f, ∂f, ∂∂f)
+end
+
 struct EntropicElasticityLaw <: ThermalLaw
   θr::Float64
   γ::Float64
@@ -32,7 +55,7 @@ struct EntropicElasticityLaw <: ThermalLaw
 end
 
 function (law::EntropicElasticityLaw)()
-  @unpack θr, γ = law
+  (; θr, γ) = law
   f(θ) = (θ/θr)^(γ+1)
   ∂f(θ) = (γ+1) * θ^γ / θr^(γ+1)
   ∂∂f(θ) = γ*(γ+1) * θ^(γ-1) / θr^(γ+1)
@@ -47,7 +70,7 @@ struct NonlinearMeltingLaw <: ThermalLaw
 end
 
 function (law::NonlinearMeltingLaw)()
-  @unpack θr, θM, γ = law
+  (; θr, θM, γ) = law
   f(θ) = (1 - (θ/θM)^(γ+1)) / (1 - (θr/θM)^(γ+1))
   ∂f(θ) = -(γ+1)*θ^γ/θM^(γ+1) / (1 - (θr/θM)^(γ+1))
   ∂∂f(θ) = -γ*(γ+1)*θ^(γ-1)/θM^(γ+1) / (1 - (θr/θM)^(γ+1))
@@ -56,35 +79,19 @@ end
 
 struct NonlinearSofteningLaw <: ThermalLaw
   θr::Float64
-  θt::Float64
+  θT::Float64
   γ::Float64
   δ::Float64
-  NonlinearSofteningLaw(; θr, θt, γ, δ=0) = new(θr, θt, γ, δ)
+  NonlinearSofteningLaw(; θr, θT, γ, δ=0) = new(θr, θT, γ, δ)
 end
 
 function (law::NonlinearSofteningLaw)()
-  @unpack θr, θt, γ, δ = law
-  u(θ) = exp(-(θ/θt)^(γ+1))
+  (; θr, θT, γ, δ) = law
+  u(θ) = exp(-(θ/θT)^(γ+1))
   C = (1-δ) * u(θr) + δ
   f(θ) = ((1-δ) * u(θ) + δ) / C
-  ∂f(θ) = -(1-δ)/C * (γ+1)/θt * (θ/θt)^γ * u(θ)
-  ∂∂f(θ) = (1-δ)/C * (γ+1)/θ^2 * (θ/θt)^(γ+1) * ((γ+1)*(θ/θt)^(γ+1)-γ) * u(θ)
-  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))
+  ∂f(θ) = -(1-δ)/C * (γ+1)/θT * (θ/θT)^γ * u(θ)
+  ∂∂f(θ) = (1-δ)/C * (γ+1)/θ^2 * (θ/θT)^(γ+1) * ((γ+1)*(θ/θT)^(γ+1)-γ) * u(θ)
   return (f, ∂f, ∂∂f)
 end
 
@@ -93,10 +100,11 @@ struct PolynomialLaw <: ThermalLaw
   a::Float64
   b::Float64
   c::Float64
+  PolynomialLaw(; θr, a, b, c) = new(θr, a, b, c)
 end
 
 function (law::PolynomialLaw)()
-  @unpack θr, a, b, c = law
+  (; θr, a, b, c) = law
   f(θ)   = a*((θ-θr)/θr)^3  + b*((θ-θr)/θr)^2 + c*(θ-θr)/θr + 1
   ∂f(θ)  = 3a*(θ-θr)^2/θr^3 + 2b*(θ-θr)/θr^2 + c/θr
   ∂∂f(θ) = 6a*(θ-θr)/θr^3 + 2b/θr^2
diff --git a/src/PhysicalModels/ThermoElectroMechanicalModels.jl b/src/PhysicalModels/ThermoElectroMechanicalModels.jl
index 1fff381..13092fd 100644
--- a/src/PhysicalModels/ThermoElectroMechanicalModels.jl
+++ b/src/PhysicalModels/ThermoElectroMechanicalModels.jl
@@ -1,6 +1,11 @@
 
+function Gridap.CellData.CellState(obj::ThermoElectroMechano, args...)
+  CellState(obj.mechano, args...)
+end
+
 function initialize_state(obj::ThermoElectroMechano, 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::ThermoElectroMechano, state, F, E, θ, args...)
@@ -31,7 +36,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 +78,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 +139,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 +168,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()
diff --git a/src/PhysicalModels/ThermoElectroModels.jl b/src/PhysicalModels/ThermoElectroModels.jl
new file mode 100644
index 0000000..522be7b
--- /dev/null
+++ b/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
diff --git a/src/PhysicalModels/ThermoMechanicalModels.jl b/src/PhysicalModels/ThermoMechanicalModels.jl
index 5ad7af1..4eba993 100644
--- a/src/PhysicalModels/ThermoMechanicalModels.jl
+++ b/src/PhysicalModels/ThermoMechanicalModels.jl
@@ -3,8 +3,13 @@
 # Common functions
 # ===================
 
+function Gridap.CellData.CellState(obj::ThermoMechano, args...)
+  CellState(obj.mechano, args...)
+end
+
 function initialize_state(obj::TM, points::Measure) where {TM<:ThermoMechano}
-  initialize_state(obj.mechano, points)
+  @warn "The function 'initialize_state' is deprecated, use 'CellState' instead."
+  CellState(obj.mechano, points)
 end
 
 function update_state!(obj::TM, state, F, θ, args...) where {TM<:ThermoMechano}
@@ -41,7 +46,7 @@ struct ThermalVolumetric{T<:Thermo} <: ThermoMechano{T,Volumetric}
 end
 
 function (obj::ThermalVolumetric)()
-  @unpack Cv, θr, α, κ = obj.thermo
+  (; Cv, θr, α, κ) = obj.thermo
   cv0 = Cv  # FIXME
   U, ∂U∂F, ∂∂U∂FF = obj.mechano()
   κr = tangent(obj.mechano)
@@ -64,6 +69,28 @@ function (obj::ThermalVolumetric)()
 end
 
 
+struct ThermalDeviatoric{M<:Mechano} <: ThermoMechano{Nothing,M}
+  mechano::M
+  law::ThermalLaw
+
+  function ThermalDeviatoric(mechano::M, law::ThermalLaw) where {M<:Mechano}
+    new{M}(mechano, law)
+  end
+end
+
+function (obj::ThermalDeviatoric{<:IsoElastic})()
+  Ψm, ∂Ψm∂F, ∂∂Ψm∂FF = obj.mechano()
+  f, df, ddf = obj.law()
+  Ψ(F,θ)      = Ψm(F) * f(θ)
+  ∂Ψ∂F(F,θ)   = ∂Ψm∂F(F) * f(θ)
+  ∂∂Ψ∂FF(F,θ) = ∂∂Ψm∂FF(F) * f(θ)
+  ∂Ψ∂θ(F,θ)   = Ψm(F) * df(θ)
+  ∂∂Ψ∂θθ(F,θ) = Ψm(F) * ddf(θ)
+  ∂∂Ψ∂Fθ(F,θ) = ∂Ψm∂F(F) * df(θ)
+  return (Ψ, ∂Ψ∂F, ∂Ψ∂θ, ∂∂Ψ∂FF, ∂∂Ψ∂θθ, ∂∂Ψ∂Fθ)
+end
+
+
 struct ThermoMechModel{T<:Thermo,M<:Mechano} <: ThermoMechano{T,M}
   thermo::T
   mechano::M
diff --git a/src/PhysicalModels/ViscousModels.jl b/src/PhysicalModels/ViscousModels.jl
index 4113aa9..d247c0c 100644
--- a/src/PhysicalModels/ViscousModels.jl
+++ b/src/PhysicalModels/ViscousModels.jl
@@ -26,11 +26,20 @@ function update_time_step!(obj::ViscousIncompressible, Δt::Float64)
   obj.Δt[] = Δt
 end
 
-function initialize_state(::ViscousIncompressible, points::Measure)
-  v = VectorValue(1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0,0.0)
+function Gridap.CellData.CellState(::ViscousIncompressible, F0::TensorValue, points::Measure)
+  v = VectorValue(F0..., 0.0)
   CellState(v, points)
 end
 
+function Gridap.CellData.CellState(obj::ViscousIncompressible, points::Measure)
+  CellState(obj, I3, points)
+end
+
+function initialize_state(obj::ViscousIncompressible, points::Measure)
+  @warn "The function 'initialize_state' is deprecated, use 'CellState' instead."
+  CellState(obj, points)
+end
+
 function update_state!(obj::ViscousIncompressible, state, F, Fn)
   _, Se, ∂Se∂Ce = SecondPiola(obj.elasto)
   return_mapping(A, F, Fn) = ReturnMapping(obj, Se, ∂Se∂Ce, F, Fn, A)
@@ -79,8 +88,13 @@ function update_time_step!(obj::NVisco, Δt::Float64)
   Δt
 end
 
+function Gridap.CellData.CellState(obj::NVisco, args...)
+  map(b -> CellState(b, args...), obj)
+end
+
 function initialize_state(obj::NVisco, points::Measure)
-  map(b -> initialize_state(b, points), obj)
+  @warn "The function 'initialize_state' is deprecated, use 'CellState' instead."
+  map(b -> CellState(b, points), obj)
 end
 
 function update_state!(obj::NVisco, states, F, Fn)
@@ -125,8 +139,13 @@ function update_time_step!(obj::GeneralizedMaxwell, Δt::Float64)
   update_time_step!(obj.branches, Δt)
 end
 
+function Gridap.CellData.CellState(obj::GeneralizedMaxwell, args...)
+  CellState(obj.branches, args...)
+end
+
 function initialize_state(obj::GeneralizedMaxwell, points::Measure)
-  initialize_state(obj.branches, points)
+  @warn "The function 'initialize_state' is deprecated, use 'CellState' instead."
+  CellState(obj.branches, points)
 end
 
 function update_state!(obj::GeneralizedMaxwell{<:IsoElastic}, states, F, Fn)
diff --git a/test/TestConstitutiveModels/ElectroMechanicalTests.jl b/test/TestConstitutiveModels/ElectroMechanicalTests.jl
index ed703f3..815e467 100644
--- a/test/TestConstitutiveModels/ElectroMechanicalTests.jl
+++ b/test/TestConstitutiveModels/ElectroMechanicalTests.jl
@@ -1,6 +1,7 @@
 using Gridap.TensorValues
 using HyperFEM.PhysicalModels
 using HyperFEM.TensorAlgebra
+using ForwardDiff
 
 
 const ∇φ = VectorValue(1.0:3.0...)
@@ -8,6 +9,20 @@ const ∇u = TensorValue(1.0:9.0...) * 1e-3
 const ∇un = TensorValue(1.0:9.0...) * 5e-4
 
 
+@testset "IdealDielectric" begin
+  model = IdealDielectric(ε=4.0*8.85e-12)
+  Ψ, ∂Ψ∂F, ∂Ψ∂E, ∂Ψ∂FF, ∂Ψ∂EF, ∂Ψ∂EE = model()
+  E0 = VectorValue(rand(3)) * 6000 / 0.001
+  F1 = I3 + 0.1*TensorValue(rand(9)...)
+  @test get_array(∂Ψ∂F(F1,E0))  ≈ ForwardDiff.gradient(Fi -> Ψ(Fi, get_array(E0)), get_array(F1))
+  @test get_array(∂Ψ∂E(F1,E0))  ≈ ForwardDiff.gradient(Ei -> Ψ(get_array(F1), Ei), get_array(E0))
+  @test get_array(∂Ψ∂FF(F1,E0)) ≈ ForwardDiff.jacobian(Fi -> ∂Ψ∂F(Fi, get_array(E0)), get_array(F1))
+  @test get_array(∂Ψ∂EF(F1,E0)) ≈ ForwardDiff.jacobian(Fi -> ∂Ψ∂E(Fi, get_array(E0)), get_array(F1))
+  @test get_array(∂Ψ∂EF(F1,E0)) ≈ ForwardDiff.jacobian(Ei -> ∂Ψ∂F(get_array(F1), Ei), get_array(E0))'
+  @test get_array(∂Ψ∂EE(F1,E0)) ≈ ForwardDiff.jacobian(Ei -> ∂Ψ∂E(get_array(F1), Ei), get_array(E0))
+end
+
+
 @testset "Electro+4*HGO_1Fiber" begin
   c1  =  [0.6639232500447778, 0.5532987701062146, 0.9912576142028674, 0.4951942011240962]
   c2  =  [0.800583033264982, 0.3141082734275339, 0.8063905248474006, 0.5850486948450955]
diff --git a/test/TestConstitutiveModels/ThermalLawsTests.jl b/test/TestConstitutiveModels/ThermalLawsTests.jl
index 5fc74c9..e5877c3 100644
--- a/test/TestConstitutiveModels/ThermalLawsTests.jl
+++ b/test/TestConstitutiveModels/ThermalLawsTests.jl
@@ -2,12 +2,30 @@ using ForwardDiff
 using HyperFEM
 using Test
 
+@testset "ConstantEnergyLaw" begin
+  law = ConstantEnergyLaw()
+  f, df, ddf = law()
+  for θ ∈ 200.0:50:400
+    @test df(θ) ≈ ForwardDiff.derivative(f, θ)
+    @test ddf(θ) ≈ ForwardDiff.derivative(df, θ)
+  end
+end
+
+@testset "ConstantCvLaw" begin
+  law = ConstantCvLaw(θr=273.15)
+  f, df, ddf = law()
+  for θ ∈ 200.0:50:400
+    @test df(θ) ≈ ForwardDiff.derivative(f, θ)
+    @test ddf(θ) ≈ ForwardDiff.derivative(df, θ)
+  end
+end
+
 @testset "EntropicElasticityLaw" begin
   law = EntropicElasticityLaw(θr=273.15, γ=0.55)
   f, df, ddf = law()
   for θ ∈ 200.0:50:400
-    @test isapprox(df(θ), ForwardDiff.derivative(f, θ), rtol=1e-3)
-    @test isapprox(ddf(θ), ForwardDiff.derivative(df, θ), rtol=1e-10)
+    @test df(θ) ≈ ForwardDiff.derivative(f, θ)
+    @test ddf(θ) ≈ ForwardDiff.derivative(df, θ)
   end
 end
 
@@ -15,25 +33,25 @@ end
   law = NonlinearMeltingLaw(θr=273.15, θM=400.0, γ=0.55)
   f, df, ddf = law()
   for θ ∈ 200.0:50:400
-    @test isapprox(df(θ), ForwardDiff.derivative(f, θ), rtol=1e-10)
-    @test isapprox(ddf(θ), ForwardDiff.derivative(df, θ), rtol=1e-10)
+    @test df(θ) ≈ ForwardDiff.derivative(f, θ)
+    @test ddf(θ) ≈ ForwardDiff.derivative(df, θ)
   end
 end
 
 @testset "NonlinearSofteningLaw" begin
-  law = NonlinearSofteningLaw(θr=273.15, θt=300.0, γ=2.0, δ=0.5)
+  law = NonlinearSofteningLaw(θr=273.15, θT=300.0, γ=2.0, δ=0.5)
   f, df, ddf = law()
   for θ ∈ 200.0:50:400
-    @test isapprox(df(θ), ForwardDiff.derivative(f, θ), rtol=1e-10)
-    @test isapprox(ddf(θ), ForwardDiff.derivative(df, θ), rtol=1e-10)
+    @test df(θ) ≈ ForwardDiff.derivative(f, θ)
+    @test ddf(θ) ≈ ForwardDiff.derivative(df, θ)
   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)
-    @test isapprox(ddf(θ), ForwardDiff.derivative(df, θ), rtol=1e-10)
+    @test df(θ) ≈ ForwardDiff.derivative(f, θ)
+    @test ddf(θ) ≈ ForwardDiff.derivative(df, θ)
   end
 end
diff --git a/test/data/StaggeredViscoElectricSimulation.jl b/test/data/StaggeredViscoElectricSimulation.jl
index 7cac8c2..74ea0cd 100644
--- a/test/data/StaggeredViscoElectricSimulation.jl
+++ b/test/data/StaggeredViscoElectricSimulation.jl
@@ -80,7 +80,7 @@ function staggered_visco_electric_simulation(; t_end=2, writevtk=true, verbose=t
 
   Fh  = F∘∇(uh⁺)'
   Fh⁻ = F∘∇(uh⁻)'
-  A   = initialize_state(cons_model, dΩ)
+  A   = CellState(cons_model, dΩ)
 
   # Electrical staggered step
   res_elec(Λ) = (φ, vφ) -> residual(cons_model, Electro, (ku, ke), (uh⁺, φ), vφ, dΩ, 0.0, Fh⁻, A...)
diff --git a/test/data/ViscoElasticSimulation.jl b/test/data/ViscoElasticSimulation.jl
index 14e67e3..d14ea16 100644
--- a/test/data/ViscoElasticSimulation.jl
+++ b/test/data/ViscoElasticSimulation.jl
@@ -66,7 +66,7 @@ function visco_elastic_simulation(;t_end=15, writevtk=true, verbose=true)
 
   uh = FEFunction(Uu, zero_free_values(Uu))
   unh = FEFunction(Uun, zero_free_values(Uun))
-  state_vars = initialize_state(cons_model, dΩ)
+  state_vars = CellState(cons_model, dΩ)
   F,_,_ = get_Kinematics(k)
   Fnh = F∘∇(unh)'