diff --git a/src/PhysicalModels/ElectroMechanicalModels.jl b/src/PhysicalModels/ElectroMechanicalModels.jl
index 700ee8e..a76f281 100644
--- a/src/PhysicalModels/ElectroMechanicalModels.jl
+++ b/src/PhysicalModels/ElectroMechanicalModels.jl
@@ -33,7 +33,7 @@ struct ElectroMechModel{E<:Electro,M<:Mechano} <: ElectroMechano
     ∂Ψφφ(F, E, N) = ∂Ψem_φφ(F, E)
     return (Ψ, ∂Ψu, ∂Ψφ, ∂Ψuu, ∂Ψφu, ∂Ψφφ)
   end
-  function (obj::ElectroMechModel{<:Electro,<:ViscoElastic})(Λ::Float64=1.0)
+  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)
@@ -44,6 +44,17 @@ struct ElectroMechModel{E<:Electro,M<:Mechano} <: ElectroMechano
     ∂Ψφφ(F, E, Fn, A...) = ∂Ψem_φφ(F, E)
     return (Ψ, ∂Ψu, ∂Ψφ, ∂Ψuu, ∂Ψφu, ∂Ψφφ)
   end
+  function (obj::ElectroMechModel{<:Electro,<:ViscoElastic{<:AnisoElastic}})()
+    Ψ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
 end
 
 ViscoElectricModel = ElectroMechModel{<:Electro,<:ViscoElastic}
@@ -53,12 +64,12 @@ function update_time_step!(obj::ElectroMechModel, Δt::Float64)
   update_time_step!(obj.mechano, Δt)
 end
 
-function initialize_state(obj::ElectroMechano, points::Measure)
+function initialize_state(obj::ElectroMechModel, points::Measure)
   initialize_state(obj.mechano, points)
 end
 
-function update_state!(obj::ElectroMechano, args...)
-  update_state!(obj.mechano, args...)
+function update_state!(obj::ElectroMechModel, state, F, E, args...)
+  update_state!(obj.mechano, state, F, args...)
 end
 
 
diff --git a/src/PhysicalModels/PhysicalModels.jl b/src/PhysicalModels/PhysicalModels.jl
index 1b67e84..125cae5 100644
--- a/src/PhysicalModels/PhysicalModels.jl
+++ b/src/PhysicalModels/PhysicalModels.jl
@@ -105,7 +105,7 @@ abstract type Elasto <: Mechano end
 abstract type IsoElastic <: Elasto end
 abstract type AnisoElastic <: Elasto end
 abstract type Visco <: Mechano end
-abstract type ViscoElastic <: Mechano end
+abstract type ViscoElastic{E<:Elasto} <: Mechano end
 
 abstract type MultiPhysicalModel <: PhysicalModel end
 abstract type ElectroMechano <: MultiPhysicalModel end
diff --git a/src/PhysicalModels/ViscousModels.jl b/src/PhysicalModels/ViscousModels.jl
index a38d1d5..9f7fbd8 100644
--- a/src/PhysicalModels/ViscousModels.jl
+++ b/src/PhysicalModels/ViscousModels.jl
@@ -42,26 +42,39 @@ function Dissipation(obj::ViscousIncompressible)
   D(F, Fn, A) = ViscousDissipation(obj, Se, ∂Se∂Ce, F, Fn, A)
 end
 
-struct GeneralizedMaxwell <: ViscoElastic
-  longterm::Elasto
+struct GeneralizedMaxwell{E<:Elasto} <: ViscoElastic{E}
+  longterm::E
   branches::NTuple{N,Visco} where N
   Δt::Ref{Float64}
-  function GeneralizedMaxwell(longTerm::Elasto, branches::Visco...)
-    new(longTerm,branches,0)
-  end
-  function (obj::GeneralizedMaxwell)()
-    Ψe, ∂ΨeF, ∂ΨeFF = obj.longterm()
-    DΨv   = map(b -> b(), obj.branches)
-    Ψα    = map(x -> x[1], DΨv)
-    ∂ΨαF  = map(x -> x[2], DΨv)
-    ∂ΨαFF = map(x -> x[3], DΨv)
-    Ψ(F, Fn, A...)     = mapreduce((Ψi, Ai) -> Ψi(F, Fn, Ai), +, Ψα, A; init=Ψe(F))
-    ∂Ψ∂F(F, Fn, A...)  = mapreduce((∂ΨiF, Ai) -> ∂ΨiF(F, Fn, Ai), +, ∂ΨαF, A; init=∂ΨeF(F))
-    ∂Ψ∂FF(F, Fn, A...) = mapreduce((∂ΨiFF, Ai) -> ∂ΨiFF(F, Fn, Ai), +, ∂ΨαFF, A; init=∂ΨeFF(F))
-    (Ψ, ∂Ψ∂F, ∂Ψ∂FF)
+  function GeneralizedMaxwell(longTerm::E, branches::Vararg{Visco}) where {E<:Elasto}
+    new{E}(longTerm,branches,0)
   end
 end
 
+function (obj::GeneralizedMaxwell{<:IsoElastic})()
+  Ψe, ∂ΨeF, ∂ΨeFF = obj.longterm()
+  DΨv   = map(b -> b(), obj.branches)
+  Ψα    = map(x -> x[1], DΨv)
+  ∂ΨαF  = map(x -> x[2], DΨv)
+  ∂ΨαFF = map(x -> x[3], DΨv)
+  Ψ(F, Fn, A...)     = mapreduce((Ψi, Ai) -> Ψi(F, Fn, Ai), +, Ψα, A; init=Ψe(F))
+  ∂Ψ∂F(F, Fn, A...)  = mapreduce((∂ΨiF, Ai) -> ∂ΨiF(F, Fn, Ai), +, ∂ΨαF, A; init=∂ΨeF(F))
+  ∂Ψ∂FF(F, Fn, A...) = mapreduce((∂ΨiFF, Ai) -> ∂ΨiFF(F, Fn, Ai), +, ∂ΨαFF, A; init=∂ΨeFF(F))
+  (Ψ, ∂Ψ∂F, ∂Ψ∂FF)
+end
+
+function (obj::GeneralizedMaxwell{<:AnisoElastic})()
+  Ψe, ∂ΨeF, ∂ΨeFF = obj.longterm()
+  DΨv   = map(b -> b(), obj.branches)
+  Ψα    = map(x -> x[1], DΨv)
+  ∂ΨαF  = map(x -> x[2], DΨv)
+  ∂ΨαFF = map(x -> x[3], DΨv)
+  Ψ(F, n, Fn, A...)     = mapreduce((Ψi, Ai) -> Ψi(F, Fn, Ai), +, Ψα, A; init=Ψe(F,n))
+  ∂Ψ∂F(F, n, Fn, A...)  = mapreduce((∂ΨiF, Ai) -> ∂ΨiF(F, Fn, Ai), +, ∂ΨαF, A; init=∂ΨeF(F,n))
+  ∂Ψ∂FF(F, n, Fn, A...) = mapreduce((∂ΨiFF, Ai) -> ∂ΨiFF(F, Fn, Ai), +, ∂ΨαFF, A; init=∂ΨeFF(F,n))
+  (Ψ, ∂Ψ∂F, ∂Ψ∂FF)
+end
+
 function update_time_step!(obj::GeneralizedMaxwell, Δt::Float64)
   foreach(b -> update_time_step!(b, Δt), obj.branches)
   obj.Δt[] = Δt
@@ -71,16 +84,26 @@ function initialize_state(obj::GeneralizedMaxwell, points::Measure)
   map(b -> initialize_state(b, points), obj.branches)
 end
 
-function update_state!(obj::GeneralizedMaxwell, states, F, Fn)
+function update_state!(obj::GeneralizedMaxwell{<:IsoElastic}, states, F, Fn)
+  @assert length(obj.branches) == length(states)
+  map((b, s) -> update_state!(b, s, F, Fn), obj.branches, states)
+end
+
+function update_state!(obj::GeneralizedMaxwell{<:AnisoElastic}, states, F, N, Fn)
   @assert length(obj.branches) == length(states)
   map((b, s) -> update_state!(b, s, F, Fn), obj.branches, states)
 end
 
-function Dissipation(obj::GeneralizedMaxwell)
+function Dissipation(obj::GeneralizedMaxwell{<:IsoElastic})
   Dα = map(Dissipation, obj.branches)
   D(F, Fn, A...) = mapreduce((Di, Ai) -> Di(F, Fn, Ai), +, Dα, A)
 end
 
+function Dissipation(obj::GeneralizedMaxwell{<:AnisoElastic})
+  Dα = map(Dissipation, obj.branches)
+  D(F, N, Fn, A...) = mapreduce((Di, Ai) -> Di(F, Fn, Ai), +, Dα, A)
+end
+
 
 # =====================
 # Internal functions