diff --git a/src/Exports.jl b/src/Exports.jl
index 22dea9c..0caca57 100644
--- a/src/Exports.jl
+++ b/src/Exports.jl
@@ -88,8 +88,6 @@ end
 
 @publish PhysicalModels SecondPiola
 @publish PhysicalModels Dissipation
-@publish PhysicalModels initializeStateVariables
-@publish PhysicalModels updateStateVariables!
 @publish PhysicalModels initialize_state
 @publish PhysicalModels update_time_step!
 
diff --git a/src/PhysicalModels/ElectroMechanicalModels.jl b/src/PhysicalModels/ElectroMechanicalModels.jl
index 15cd126..700ee8e 100644
--- a/src/PhysicalModels/ElectroMechanicalModels.jl
+++ b/src/PhysicalModels/ElectroMechanicalModels.jl
@@ -61,14 +61,6 @@ function update_state!(obj::ElectroMechano, args...)
   update_state!(obj.mechano, args...)
 end
 
-function initializeStateVariables(obj::ElectroMechano, points::Measure)
-  initializeStateVariables(obj.mechano, points)
-end
-
-function updateStateVariables!(state, obj::ElectroMechano, args...)
-  updateStateVariables!(state, obj.mechano, args)
-end
-
 
 function _getCoupling(elec::Electro, mec::Mechano, Λ::Float64)
   J(F) = det(F)
diff --git a/src/PhysicalModels/MechanicalModels.jl b/src/PhysicalModels/MechanicalModels.jl
index e7777b2..6a80cba 100644
--- a/src/PhysicalModels/MechanicalModels.jl
+++ b/src/PhysicalModels/MechanicalModels.jl
@@ -144,18 +144,16 @@ Base.hcat(a::AnisoElastic...) = MultiAnisoElastic(a)
 
 
 function (obj::MultiAnisoElastic)(args...)
-  DΨ = map(a -> a(args...), obj.Models)
-  Ψα, ∂Ψα∂F, ∂Ψα∂FF = transpose(DΨ)
-  Ψ(F, N) = mapreduce((Ψi, Ni) -> Ψi(F, Ni), +, Ψα, N)
-  ∂Ψ∂F(F, N) = mapreduce((∂Ψi∂F, Ni) -> ∂Ψi∂F(F, Ni), +, ∂Ψα∂F, N)
+  DΨ     = map(a -> a(args...), obj.Models)
+  Ψα     = map(x -> x[1], DΨ)
+  ∂Ψα∂F  = map(x -> x[2], DΨ)
+  ∂Ψα∂FF = map(x -> x[3], DΨ)
+  Ψ(F, N)     = mapreduce((Ψi, Ni) -> Ψi(F, Ni), +, Ψα, N)
+  ∂Ψ∂F(F, N)  = mapreduce((∂Ψi∂F, Ni) -> ∂Ψi∂F(F, Ni), +, ∂Ψα∂F, N)
   ∂Ψ∂FF(F, N) = mapreduce((∂Ψi∂FF, Ni) -> ∂Ψi∂FF(F, Ni), +, ∂Ψα∂FF, N)
   (Ψ, ∂Ψ∂F, ∂Ψ∂FF)
 end
 
-transpose(x::NTuple{N,<:Tuple{<:Function,<:Function,<:Function}}) where N = map(i -> getindex.(x, i), 1:3)
-
-
-
 
 # ===================
 # Mechanical models
diff --git a/src/PhysicalModels/PhysicalModels.jl b/src/PhysicalModels/PhysicalModels.jl
index 5b7c5a5..1b67e84 100644
--- a/src/PhysicalModels/PhysicalModels.jl
+++ b/src/PhysicalModels/PhysicalModels.jl
@@ -79,8 +79,6 @@ export Dissipation
 
 export DerivativeStrategy
 
-export initializeStateVariables
-export updateStateVariables!
 export initialize_state
 export update_time_step!
 
@@ -151,9 +149,6 @@ Initialize the state variables for the given constitutive model and discretizati
 function initialize_state(::PhysicalModel, points::Measure)
   return nothing
 end
-function initializeStateVariables(::PhysicalModel, points::Measure)
-  return nothing
-end
 
 
 """
@@ -161,8 +156,6 @@ Update the state variables. The state variables must be initialized using the fu
 """
 function update_state!(::PhysicalModel, vars...)
 end
-function updateStateVariables!(::Any, ::PhysicalModel, vars...)
-end
 
 
 """
diff --git a/src/PhysicalModels/ViscousModels.jl b/src/PhysicalModels/ViscousModels.jl
index 17460cb..a38d1d5 100644
--- a/src/PhysicalModels/ViscousModels.jl
+++ b/src/PhysicalModels/ViscousModels.jl
@@ -37,17 +37,6 @@ function update_state!(obj::ViscousIncompressible, state, F, Fn)
   update_state!(return_mapping, state, F, Fn)
 end
 
-function initializeStateVariables(::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)
-  CellState(v, points)
-end
-
-function updateStateVariables!(state, obj::ViscousIncompressible, F, Fn)
-  _, Se, ∂Se∂Ce = SecondPiola(obj.elasto)
-  return_mapping(A, F, Fn) = ReturnMapping(obj, Se, ∂Se∂Ce, F, Fn, A)
-  update_state!(return_mapping, state, F, Fn)
-end
-
 function Dissipation(obj::ViscousIncompressible)
   _, Se, ∂Se∂Ce = SecondPiola(obj.elasto)
   D(F, Fn, A) = ViscousDissipation(obj, Se, ∂Se∂Ce, F, Fn, A)
@@ -61,14 +50,15 @@ struct GeneralizedMaxwell <: ViscoElastic
     new(longTerm,branches,0)
   end
   function (obj::GeneralizedMaxwell)()
-    Ψe, ∂Ψeu, ∂Ψeuu = obj.longterm()
-    DΨv = map(b -> b(), obj.branches)
-    Ψα, ∂Ψαu, ∂Ψαuu = map(i -> getindex.(DΨv, i), 1:3)
-    Ψα, ∂Ψαu, ∂Ψαuu = transpose(DΨv)
-    Ψ(F, Fn, A...) = mapreduce((Ψi, Ai) -> Ψi(F, Fn, Ai), +, Ψα, A; init=Ψe(F))
-    ∂Ψu(F, Fn, A...) = mapreduce((∂Ψiu, Ai) -> ∂Ψiu(F, Fn, Ai), +, ∂Ψαu, A; init=∂Ψeu(F))
-    ∂Ψuu(F, Fn, A...) = mapreduce((∂Ψiuu, Ai) -> ∂Ψiuu(F, Fn, Ai), +, ∂Ψαuu, A; init=∂Ψeuu(F))
-    return (Ψ, ∂Ψu, ∂Ψuu)
+    Ψ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
 end
 
@@ -86,17 +76,6 @@ function update_state!(obj::GeneralizedMaxwell, states, F, Fn)
   map((b, s) -> update_state!(b, s, F, Fn), obj.branches, states)
 end
 
-function initializeStateVariables(obj::GeneralizedMaxwell, points::Measure)
-  map(b -> initializeStateVariables(b, points), obj.branches)
-end
-
-function updateStateVariables!(states, obj::GeneralizedMaxwell, F, Fn)
-  @assert length(obj.branches) == length(states)
-  for (branch, state) in zip(obj.branches, states)
-    updateStateVariables!(state, branch, F, Fn)
-  end
-end
-
 function Dissipation(obj::GeneralizedMaxwell)
   Dα = map(Dissipation, obj.branches)
   D(F, Fn, A...) = mapreduce((Di, Ai) -> Di(F, Fn, Ai), +, Dα, A)