diff --git a/src/ComputationalModels/BoundaryConditions.jl b/src/ComputationalModels/BoundaryConditions.jl
index 59f397d..89c52a4 100644
--- a/src/ComputationalModels/BoundaryConditions.jl
+++ b/src/ComputationalModels/BoundaryConditions.jl
@@ -10,6 +10,9 @@ struct MultiFieldBC <: BoundaryCondition
     BoundaryCondition::Vector{BoundaryCondition}
 end
 
+include("EvolutionFunctions.jl")
+include("CartesianTags.jl")
+
 
 struct MultiFieldTC{A} <: TimedependentCondition
     vh::A # could be a multifield or single field
diff --git a/src/ComputationalModels/CartesianTags.jl b/src/ComputationalModels/CartesianTags.jl
new file mode 100644
index 0000000..27badca
--- /dev/null
+++ b/src/ComputationalModels/CartesianTags.jl
@@ -0,0 +1,46 @@
+
+"Shortcuts for the tags of cartesian discrete models."
+module CartesianTags
+
+  "Tags indicating points, edges and faces at X0."
+  const faceX0 = [1,3,5,7,13,15,17,19,25]
+
+  "Tags indicating points, edges and faces at X1."
+  const faceX1 = [2,4,6,8,14,16,18,20,26]
+
+  "Tags indicating points, edges and faces at Y0."
+  const faceY0 = [1,2,5,6,9,11,17,18,23]
+
+  "Tags indicating points, edges and faces at Y1."
+  const faceY1 = [3,4,7,8,10,12,19,20,24]
+
+  "Tags indicating points, edges and faces at Z0."
+  const faceZ0 = [1,2,3,4,9,10,13,14,21]
+
+  "Tags indicating points, edges and faces at Z1."
+  const faceZ1 = [5,6,7,8,11,12,15,16,22]
+
+  "Tag indicating the point at corner X0, Y0, Z0."
+  const corner000 = [1]
+
+  "Tag indicating the point at corner X1, Y0, Z0."
+  const corner100 = [2]
+
+  "Tag indicating the point at corner X0, Y1, Z0."
+  const corner010 = [3]
+
+  "Tag indicating the point at corner X1, Y1, Z0."
+  const corner110 = [4]
+
+  "Tag indicating the point at corner X0, Y0, Z1."
+  const corner001 = [5]
+
+  "Tag indicating the point at corner X1, Y0, Z1."
+  const corner101 = [6]
+
+  "Tag indicating the point at corner X0, Y1, Z1."
+  const corner011 = [7]
+
+  "Tag indicating the point at corner X1, Y1, Z1."
+  const corner111 = [8]
+end
diff --git a/src/ComputationalModels/EvolutionFunctions.jl b/src/ComputationalModels/EvolutionFunctions.jl
new file mode 100644
index 0000000..e77e2e6
--- /dev/null
+++ b/src/ComputationalModels/EvolutionFunctions.jl
@@ -0,0 +1,42 @@
+
+
+"Return a triangular evolution function ranging from 0 to 1, centered at T, having edges at 0 and 2T."
+function triangular(T::Float64)
+  triangular(0.0, T)
+end
+
+"Return a triangular evolution function ranging from 0 to 1, centered at Tmax, having edges at T0 and 2Tmax-T0."
+function triangular(T0::Float64, Tmax::Float64)
+  t::Float64 -> begin
+    Δ = Tmax - T0
+    u = (t - T0) / Δ
+    v = (t - Tmax) / Δ
+    max(min(u, 1.0-v), .0)
+  end
+end
+
+"Return the Heaviside function."
+function step(T::Float64)
+  t::Float64 -> t > T ? 1.0 : 0.0
+end
+
+"Return a sigmoid-like function with edges at 0 and 2ϵ."
+function smoothstep(ϵ::Float64)
+  smoothstep(ϵ, ϵ)    
+end
+
+"Return a sigmoid-like function centered at T and edges at ±ϵ."
+function smoothstep(T::Float64, ϵ::Float64)
+  t::Float64 -> begin
+    u::Float64 = (t - T + ϵ) / (2 * ϵ)
+    if u < 0.0 return 0.0
+    elseif u < 1.0 return 3*u^2 - 2*u^3
+    else return 1.0
+    end
+  end
+end
+
+"Return a constant function which is always evaluated to 1."
+function constant()
+  t::Float64 -> 1.0
+end
diff --git a/src/ComputationalModels/PostProcessors.jl b/src/ComputationalModels/PostProcessors.jl
index e5f6584..6cc52c4 100644
--- a/src/ComputationalModels/PostProcessors.jl
+++ b/src/ComputationalModels/PostProcessors.jl
@@ -97,8 +97,16 @@ function Cauchy(physmodel::ThermoElectroMechano, uh, φh, θh, Ω, dΩ, Λ=1.0)
 end
 
 
+function Cauchy(model::Elasto, uh, unh, state_vars, Ω, dΩ, t, Δt)
+    _, ∂Ψu, _ = model(t)
+    F, _, _ = get_Kinematics(model.Kinematic)
+    σ = ∂Ψu ∘ (F∘∇(uh))
+    return interpolate_σ_everywhere(σ, Ω, dΩ)
+end
+
+
 function Cauchy(model::ViscoElastic, uh, unh, state_vars, Ω, dΩ, t, Δt)
-    _, ∂Ψu, _ = model(Δt=Δt)
+    _, ∂Ψu, _ = model(t, Δt=Δt)
     F, _, _ = get_Kinematics(model.Kinematic)
     σ = ∂Ψu ∘ (F∘∇(uh), F∘∇(unh), state_vars...)
     return interpolate_σ_everywhere(σ, Ω, dΩ)
diff --git a/src/WeakForms/WeakForms.jl b/src/WeakForms/WeakForms.jl
index 216bf8f..55a00cf 100644
--- a/src/WeakForms/WeakForms.jl
+++ b/src/WeakForms/WeakForms.jl
@@ -27,43 +27,39 @@ end
 # ===================
 # Mechanics
 # ===================
- 
-function residual(physicalmodel::Mechano, u, v, dΩ, Λ=1.0)
-    DΨ= physicalmodel(Λ)
-    F,_,_ = get_Kinematics(physicalmodel.Kinematic; Λ=Λ)
-    ∂Ψu=DΨ[2]
 
-    ∫(∇(v)' ⊙ (∂Ψu ∘ (F∘(∇(u)'))))dΩ
-
-end
+"""
+    residual(...)::Gridap.CellData.Integrand
 
-function jacobian(physicalmodel::Mechano, u, du, v, dΩ, Λ=1.0)
-    DΨ= physicalmodel(Λ)
-    F,_,_ = get_Kinematics(physicalmodel.Kinematic; Λ=Λ)
-    ∂Ψuu=DΨ[3]
-    ∫(∇(v)' ⊙ ((∂Ψuu ∘ (F∘(∇(u)'))) ⊙ (∇(du)')))dΩ
+Calculate the residual using the given constitutive model and finite element functions.
+"""
+function residual(physicalmodel::Mechano, u, v, dΩ, Λ=1.0, Δt=0.0, vars...)
+  _, ∂Ψu, _ = physicalmodel(Λ)
+  F, _, _   = get_Kinematics(physicalmodel.Kinematic; Λ=Λ)
+  ∫(∇(v)' ⊙ (∂Ψu ∘ (F∘∇(u)')))dΩ
 end
 
 """
-    residual(...)::Function
+    jacobian(...)::Gridap.CellData.Integrand
 
-Return the residual for a visco-elastic model as a FUNCTION.
+Calculate the jacobian using the given constitutive model and finite element functions.
 """
-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Ω
+function jacobian(physicalmodel::Mechano, u, du, v, dΩ, Λ=1.0, Δt=0.0, vars...)
+  _, _, ∂Ψuu = physicalmodel(Λ)
+  F, _, _    = get_Kinematics(physicalmodel.Kinematic; Λ=Λ)
+  ∫(∇(v)' ⊙ ((∂Ψuu ∘ (F∘∇(u)')) ⊙ ∇(du)'))dΩ
 end
 
-"""
-    jacobian(...)::Function
+function residual(physicalmodel::ViscoElastic, u, v, dΩ, t, Δt, un, A)
+  _, ∂Ψu, _ = physicalmodel(t, Δt=Δt)
+  F, _, _   = get_Kinematics(physicalmodel.Kinematic, Λ=t)
+  ∫(∇(v)' ⊙ (∂Ψu ∘ (F∘∇(u)', F∘∇(un)', A...)))dΩ
+end
 
-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Ω
+function jacobian(physicalmodel::ViscoElastic, u, du, v, dΩ, t, Δt, un, A)
+  _, _, ∂Ψuu = physicalmodel(t, Δt=Δt)
+  F, _, _   = get_Kinematics(physicalmodel.Kinematic, Λ=t)
+  ∫(∇(v)' ⊙ (inner ∘ (∂Ψuu∘(F∘∇(u)', F∘∇(un)', A...), ∇(du)')))dΩ
 end