diff --git a/src/ComputationalModels/BoundaryConditions.jl b/src/ComputationalModels/BoundaryConditions.jl
index 663d198..59f397d 100644
--- a/src/ComputationalModels/BoundaryConditions.jl
+++ b/src/ComputationalModels/BoundaryConditions.jl
@@ -139,6 +139,15 @@ end
 
 function residual_Neumann(::NothingBC, kwargs...) end
 
+"""
+    residual_Neumann(...)::Function
+
+Return the Neumann residual as a FUNCTION.
+"""
+function residual_Neumann(bc::NeumannBC, dΓ::Vector, Λ::Float64)
+    v -> mapreduce((fi, dΓi) -> ∫(v ⋅ fi(Λ))dΓi, +, bc.values, dΓ)
+end
+
 function residual_Neumann(bc::NeumannBC, v, dΓ, Λ)
     bc_func_ = Vector{Function}(undef, length(bc.tags))
     for (i, f) in enumerate(bc.values)
@@ -155,20 +164,22 @@ function residual_Neumann(bc::NeumannBC, v, dΓ, Λ⁺, Λ⁻)
     return mapreduce(f -> f(v), +, bc_func_)
 end
 
-function get_Neumann_dΓ(model, ::NothingBC, degree)
+"""
+    get_Neumann_dΓ(...)::Vector{Gridap.CellData.GenericMeasure}
+
+Return a collection of boundary triangulations at the specified Neumann boundaries.
+"""
+function get_Neumann_dΓ(model, ::NothingBC, degree::Int)
     Vector{Gridap.CellData.GenericMeasure}(undef, 1)
 end
 
-function get_Neumann_dΓ(model, bc::NeumannBC, degree)
-    dΓ = Vector{Gridap.CellData.GenericMeasure}(undef, length(bc.tags))
-    for i in 1:length(bc.tags)
-        Γ = BoundaryTriangulation(model, tags=bc.tags[i])
-        dΓ[i] = Measure(Γ, degree)
-    end
-    return dΓ
+function get_Neumann_dΓ(model, bc::NeumannBC, degree::Int)
+    all_Γ = map(tag -> BoundaryTriangulation(model, tags=tag), bc.tags)
+    all_dΓ = map(Γi -> Measure(Γi, degree), all_Γ)
+    all_dΓ
 end
 
-function get_Neumann_dΓ(model, bc::MultiFieldBC, degree::Int64)
+function get_Neumann_dΓ(model, bc::MultiFieldBC, degree::Int)
     dΓ = Vector{Vector{Gridap.CellData.GenericMeasure}}(undef, length(bc.BoundaryCondition))
     for (i, bc_i) in enumerate(bc.BoundaryCondition)
         dΓ[i] = get_Neumann_dΓ(model, bc_i, degree)
diff --git a/src/ComputationalModels/PostProcessors.jl b/src/ComputationalModels/PostProcessors.jl
index b4d8e46..e5f6584 100644
--- a/src/ComputationalModels/PostProcessors.jl
+++ b/src/ComputationalModels/PostProcessors.jl
@@ -97,6 +97,31 @@ function Cauchy(physmodel::ThermoElectroMechano, uh, φh, θh, Ω, dΩ, Λ=1.0)
 end
 
 
+function Cauchy(model::ViscoElastic, uh, unh, state_vars, Ω, dΩ, t, Δt)
+    _, ∂Ψu, _ = model(Δt=Δt)
+    F, _, _ = get_Kinematics(model.Kinematic)
+    σ = ∂Ψu ∘ (F∘∇(uh), F∘∇(unh), state_vars...)
+    return interpolate_σ_everywhere(σ, Ω, dΩ)
+end
+
+
+function interpolate_σ_everywhere(σ, Ω, dΩ)
+    ref_L2 = ReferenceFE(lagrangian, Float64, 0)
+    ref_fe = ReferenceFE(lagrangian, Float64, 1)
+    VL2 = FESpace(Ω, ref_L2, conformity=:L2)
+    V = FESpace(Ω, ref_fe, conformity=:H1)
+    n1 = VectorValue(1.0, 0.0, 0.0)
+    n2 = VectorValue(0.0, 1.0, 0.0)
+    n3 = VectorValue(0.0, 0.0, 1.0)
+    σ11h = interpolate_everywhere(L2_Projection(n1 ⋅ σ ⋅ n1, dΩ, VL2), V)
+    σ12h = interpolate_everywhere(L2_Projection(n1 ⋅ σ ⋅ n2, dΩ, VL2), V)
+    σ13h = interpolate_everywhere(L2_Projection(n1 ⋅ σ ⋅ n3, dΩ, VL2), V)
+    σ22h = interpolate_everywhere(L2_Projection(n2 ⋅ σ ⋅ n2, dΩ, VL2), V)
+    σ23h = interpolate_everywhere(L2_Projection(n2 ⋅ σ ⋅ n3, dΩ, VL2), V)
+    σ33h = interpolate_everywhere(L2_Projection(n3 ⋅ σ ⋅ n3, dΩ, VL2), V)
+    ph   = interpolate_everywhere(L2_Projection(tr ∘ σ, dΩ, VL2), V)
+    return (σ11h, σ12h, σ13h, σ22h, σ23h, σ33h, ph)
+end
 
 
 function Entropy(physmodel::ThermoElectroMechano, uh, φh, θh, Ω, dΩ, Λ=1.0)
diff --git a/src/WeakForms/WeakForms.jl b/src/WeakForms/WeakForms.jl
index 2c9e20d..216bf8f 100644
--- a/src/WeakForms/WeakForms.jl
+++ b/src/WeakForms/WeakForms.jl
@@ -30,7 +30,7 @@ end
  
 function residual(physicalmodel::Mechano, u, v, dΩ, Λ=1.0)
     DΨ= physicalmodel(Λ)
-       F,_,_ = get_Kinematics(physicalmodel.Kinematic; Λ=Λ)
+    F,_,_ = get_Kinematics(physicalmodel.Kinematic; Λ=Λ)
     ∂Ψu=DΨ[2]
 
     ∫(∇(v)' ⊙ (∂Ψu ∘ (F∘(∇(u)'))))dΩ
@@ -39,11 +39,38 @@ end
 
 function jacobian(physicalmodel::Mechano, u, du, v, dΩ, Λ=1.0)
     DΨ= physicalmodel(Λ)
-      F,_,_ = get_Kinematics(physicalmodel.Kinematic; Λ=Λ)
+    F,_,_ = get_Kinematics(physicalmodel.Kinematic; Λ=Λ)
     ∂Ψuu=DΨ[3]
     ∫(∇(v)' ⊙ ((∂Ψuu ∘ (F∘(∇(u)'))) ⊙ (∇(du)')))dΩ
 end
 
+"""
+    residual(...)::Function
+
+Return the residual for a visco-elastic model as a FUNCTION.
+"""
+function residual(model::ViscoElastic, un, vars, dΩ; t, Δt)
+    _, ∂Ψu, _ = model(t, Δt=Δt)
+    F, _, _   = get_Kinematics(model.Kinematic, Λ=t)
+    (u, v) -> ∫(∇(v)' ⊙ (∂Ψu∘(F∘∇(u)', F∘∇(un)', vars...)))dΩ
+end
+
+"""
+    jacobian(...)::Function
+
+Return the jacobian for a visco-elastic model as a FUNCTION.
+"""
+function jacobian(model::ViscoElastic, un, vars, dΩ; t, Δt)
+    _, _, ∂Ψuu = model(t, Δt=Δt)
+    F, _, _   = get_Kinematics(model.Kinematic, Λ=t)
+    (u, du, v) -> ∫(∇(v)' ⊙ (inner ∘ (∂Ψuu∘(F∘∇(u)', F∘∇(un)', vars...), ∇(du)')))dΩ
+end
+
+
+# ===================
+# Mass term
+# ===================
+
 function mass_term(u, v, Coeff, dΩ)
     ∫(Coeff* (u⋅v))dΩ
 end