From f7a63b9346e87c018e9fe058bffcfb70f819a47a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Miguel=20Mas=C3=B3?= <miguel.maso@upc.edu>
Date: Fri, 19 Sep 2025 14:08:07 +0200
Subject: [PATCH 1/5] added evolution functions

---
 src/ComputationalModels/BoundaryConditions.jl |  2 +
 src/ComputationalModels/EvolutionFunctions.jl | 42 +++++++++++++++++++
 2 files changed, 44 insertions(+)
 create mode 100644 src/ComputationalModels/EvolutionFunctions.jl

diff --git a/src/ComputationalModels/BoundaryConditions.jl b/src/ComputationalModels/BoundaryConditions.jl
index 59f397d..43518c8 100644
--- a/src/ComputationalModels/BoundaryConditions.jl
+++ b/src/ComputationalModels/BoundaryConditions.jl
@@ -10,6 +10,8 @@ struct MultiFieldBC <: BoundaryCondition
     BoundaryCondition::Vector{BoundaryCondition}
 end
 
+include("EvolutionFunctions.jl")
+
 
 struct MultiFieldTC{A} <: TimedependentCondition
     vh::A # could be a multifield or single field
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

From c19d63e315d766269f9e01fd503419d21ad3169f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Miguel=20Mas=C3=B3?= <miguel.maso@upc.edu>
Date: Tue, 30 Sep 2025 17:48:36 +0200
Subject: [PATCH 2/5] added cartesian tags

---
 src/ComputationalModels/BoundaryConditions.jl |  1 +
 src/ComputationalModels/CartesianTags.jl      | 46 +++++++++++++++++++
 2 files changed, 47 insertions(+)
 create mode 100644 src/ComputationalModels/CartesianTags.jl

diff --git a/src/ComputationalModels/BoundaryConditions.jl b/src/ComputationalModels/BoundaryConditions.jl
index 43518c8..89c52a4 100644
--- a/src/ComputationalModels/BoundaryConditions.jl
+++ b/src/ComputationalModels/BoundaryConditions.jl
@@ -11,6 +11,7 @@ struct MultiFieldBC <: BoundaryCondition
 end
 
 include("EvolutionFunctions.jl")
+include("CartesianTags.jl")
 
 
 struct MultiFieldTC{A} <: TimedependentCondition
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

From 66ea51fd8d8e714039dfeadc3e51c0df7051bec9 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Miguel=20Mas=C3=B3?= <miguel.maso@upc.edu>
Date: Tue, 30 Sep 2025 17:56:25 +0200
Subject: [PATCH 3/5] postprocess stress

---
 src/ComputationalModels/PostProcessors.jl | 10 +++++++++-
 1 file changed, 9 insertions(+), 1 deletion(-)

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Ω)

From 99f634d009655469947b67cbaf3c732e9960d698 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Miguel=20Mas=C3=B3?= <miguel.maso@upc.edu>
Date: Tue, 30 Sep 2025 17:57:18 +0200
Subject: [PATCH 4/5] added elastic res and jac in weak forms

---
 src/WeakForms/WeakForms.jl | 26 ++++++++++++++++++++++++--
 1 file changed, 24 insertions(+), 2 deletions(-)

diff --git a/src/WeakForms/WeakForms.jl b/src/WeakForms/WeakForms.jl
index 216bf8f..12ad578 100644
--- a/src/WeakForms/WeakForms.jl
+++ b/src/WeakForms/WeakForms.jl
@@ -47,7 +47,29 @@ end
 """
     residual(...)::Function
 
-Return the residual for a visco-elastic model as a FUNCTION.
+Return the residual for an elastic model as a Function.
+"""
+function residual(model::Elasto, un, vars, dΩ; t, Δt)
+    _, ∂Ψu, _ = model(t)
+    F, _, _   = get_Kinematics(model.Kinematic, Λ=t)
+    (u, v) -> ∫(∇(v)' ⊙ (∂Ψu ∘ (F∘∇(u)')))dΩ
+end
+
+"""
+    jacobian(...)::Function
+
+Return the jacobian for an elastic model as a Function.
+"""
+function jacobian(model::Elasto, un, vars, dΩ; t, Δt)
+    _, _, ∂Ψuu = model(t)
+    F, _, _   = get_Kinematics(model.Kinematic, Λ=t)
+    (u, du, v) -> ∫(∇(v)' ⊙ (inner ∘ (∂Ψ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)
@@ -58,7 +80,7 @@ end
 """
     jacobian(...)::Function
 
-Return the jacobian for a visco-elastic model as a 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)

From bb58e0f7359e924b59d6a1cd6aee513583959d5f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Miguel=20Mas=C3=B3?= <miguel.maso@upc.edu>
Date: Wed, 1 Oct 2025 11:34:19 +0200
Subject: [PATCH 5/5] unified interface

---
 src/WeakForms/WeakForms.jl | 66 ++++++++++++--------------------------
 1 file changed, 20 insertions(+), 46 deletions(-)

diff --git a/src/WeakForms/WeakForms.jl b/src/WeakForms/WeakForms.jl
index 12ad578..55a00cf 100644
--- a/src/WeakForms/WeakForms.jl
+++ b/src/WeakForms/WeakForms.jl
@@ -27,65 +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
-
-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Ω
-end
 
 """
-    residual(...)::Function
+    residual(...)::Gridap.CellData.Integrand
 
-Return the residual for an elastic model as a Function.
+Calculate the residual using the given constitutive model and finite element functions.
 """
-function residual(model::Elasto, un, vars, dΩ; t, Δt)
-    _, ∂Ψu, _ = model(t)
-    F, _, _   = get_Kinematics(model.Kinematic, Λ=t)
-    (u, v) -> ∫(∇(v)' ⊙ (∂Ψu ∘ (F∘∇(u)')))dΩ
+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
 
 """
-    jacobian(...)::Function
+    jacobian(...)::Gridap.CellData.Integrand
 
-Return the jacobian for an elastic model as a Function.
+Calculate the jacobian using the given constitutive model and finite element functions.
 """
-function jacobian(model::Elasto, un, vars, dΩ; t, Δt)
-    _, _, ∂Ψuu = model(t)
-    F, _, _   = get_Kinematics(model.Kinematic, Λ=t)
-    (u, du, v) -> ∫(∇(v)' ⊙ (inner ∘ (∂Ψuu∘(F∘∇(u)'), ∇(du)')))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
 
-"""
-    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Ω
+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
 
-"""
-    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Ω
+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