diff --git a/_toc.yml b/_toc.yml
index 1f7600b..d0c1c13 100644
--- a/_toc.yml
+++ b/_toc.yml
@@ -58,7 +58,10 @@ parts:
           - file: chapter2/stokes-v2
       - file: chapter2/elliptic-curl
       - file: chapter2/elliptic-div
-      - file: chapter2/linear-elasticity
+      - file: chapter2/linear-elasticity-intro
+        sections:
+          - file: chapter2/linear-elasticity-v1
+          - file: chapter2/linear-elasticity-v2
   - caption: Nonlinear Problems 
     chapters:
       - file: chapter3/poisson
diff --git a/chapter2/linear-elasticity-intro.md b/chapter2/linear-elasticity-intro.md
new file mode 100644
index 0000000..e5efecc
--- /dev/null
+++ b/chapter2/linear-elasticity-intro.md
@@ -0,0 +1,139 @@
+# Linear Elasticity
+
+## The PDE problem
+The governing equations for small elastic deformations of a body $\Omega$ can be expressed as:
+
+$$
+\begin{aligned}
+    -\nabla \cdot \sigma(u) &= f & \text{in } \Omega \\
+    \sigma(u) &= C : \epsilon(u) & \text{in } \Omega \\
+\end{aligned}
+$$
+
+where : 
+ - $\Omega$ is the domain of interest,
+ - $u$ is the displacement vector field,
+ - $\sigma(u)$ is a second-order symmetric tensor, called stress tensor, 
+ - $f$ represents the load vector field (force per unit volume), 
+ <!-- $\kappa$ and $\mu$ are Lamé's elasticity parameters for the material, $I$ denotes the identity tensor,  -->
+ - $\epsilon(u)$ is a second-order symmetric tensor, called strain tensor and by definition $\epsilon(u) := \frac{1}{2}(\nabla u + (\nabla u)^T)$.
+ - $C$ is the elasticity tensor, which relates stress and strain. 
+
+Then, the strong formulation of linear elasticity is :
+Find $u \in V$ such that
+
+$$
+\begin{aligned}
+    -\nabla \cdot \sigma (u) &= f & \text{in } & \Omega \\
+    u &= 0 & \text{on } & \partial \Omega_D \\
+    \sigma(u) \cdot n &= g_T & \text{on } & \partial \Omega_T
+\end{aligned}
+$$
+
+where : 
+- $V = \{ v \in (H^1(\Omega))^3 : v = 0 \text{ on } \partial \Omega_D \}$,
+
+- $g_T$ is the traction vector on the part $\partial \Omega_T$ of the boundary,
+
+- $n$ is the outward normal vector on the boundary $\partial \Omega$,
+
+- $\partial \Omega_D \cap \partial \Omega_T = \emptyset$ and $\partial \Omega = \partial \Omega_D \cup \partial \Omega_T$. 
+- $\partial \Omega_D$ is the Dirichlet boundary condition where the displacement is fixed to zero, and $\partial \Omega_T$ is the Neumann boundary condition where the traction is prescribed.
+
+## The Variational Formulation
+The variational formulation of the linear elasticity equations involves forming the inner product of the PDE with a vector test function $v \in V$ and integrating over the domain $\Omega$. This yields:
+
+$$
+\int_{\Omega} - \nabla \cdot \sigma(u) \cdot v \, \mathrm{d} x = \int_{\Omega} f \cdot v \, \mathrm{d} x
+$$
+
+Integrating the term $\nabla \cdot \sigma(u) \cdot v$ by parts, considering boundary conditions, we obtain:
+
+$$
+\int_{\Omega} \sigma(u) : \nabla v \, \mathrm{d} x = \int_{\Omega} f \cdot v \, \mathrm{d} x + \int_{\partial \Omega_T} g_T \cdot v \, \mathrm{d} s
+$$
+
+By using the symmetry of the stress tensor $\sigma$ and its definition from $(2)$, we can notice that :
+
+$$
+\begin{aligned}
+\int_{\Omega} \sigma(u) : \nabla v \, \mathrm{d} x &= \int_{\Omega} \sigma(u) : \epsilon(v) \, \mathrm{d} x = \int_{\Omega} C : \epsilon(u) : \epsilon(v) \, \mathrm{d} x \\ &= \int_{\Omega} \epsilon(u) : C : \epsilon(v) \, \mathrm{d} x
+\end{aligned}
+$$
+
+This leads to the following variational formulation:
+
+$$
+\boxed{
+\begin{aligned}
+&\text{Find } u \in V \text{ such that:} \\
+&\qquad a(u, v) = L(v) \quad \forall v \in V
+\end{aligned}
+}
+$$
+
+with
+
+$$
+\begin{aligned}
+&a : 
+\begin{cases}
+V \times V \rightarrow \mathbb{R} \\
+(u, v) \longmapsto \int_{\Omega} \epsilon(u) : C : \epsilon(v) \, \mathrm{d} x
+\end{cases} \\
+&L : 
+\begin{cases}
+V \rightarrow \mathbb{R} \\
+v \longmapsto \int_{\Omega} f \cdot v \, \mathrm{d} x + \int_{\partial \Omega_T} g_T \cdot v \, \mathrm{d} s
+\end{cases}
+\end{aligned}
+$$
+
+## Isotropic Materials
+For isotropic materials, the elasticity tensor $C$ can be expressed in terms of the Lamé parameters $\lambda$ and $\mu$ as follows:
+
+$$
+C := \lambda (\nabla \cdot u) I_3 + 2\mu \epsilon(u)
+$$
+Then, the stress tensor can be expressed as:
+$$\sigma(u) = \lambda (\nabla \cdot u) I_3 + 2\mu \epsilon(u)$$
+
+
+This leads to the variational formulation:
+
+$$
+\boxed{
+\begin{aligned}
+&\text{Find } u \in V \text{ such that:}
+\\
+&\qquad a(u, v) = L(v) \quad \forall v \in V
+\end{aligned}
+}
+$$
+
+with
+
+$$
+\begin{aligned}
+&a :\begin{cases}
+V \times V \rightarrow \mathbb{R} \\
+(u, v) \longmapsto \int_{\Omega} \sigma(u) : \epsilon(v) \, \mathrm{d} x
+\end{cases} \\
+&L :\begin{cases}
+V \rightarrow \mathbb{R} \\
+v \longmapsto \int_{\Omega} f \cdot v \, \mathrm{d} x + \int_{\partial \Omega_T} g_T \cdot v \, \mathrm{d} s
+\end{cases}
+\end{aligned}
+$$
+
+With this formulation, the problem is well-posed under the assumption that the material is isotropic and the boundary conditions are properly defined. While $\frac{\lambda}{\mu}$ is not too large (typically $\frac{\lambda}{\mu} \leq 10^4$), the problem remains well-posed numerically. However, as $\frac{\lambda}{\mu}$ increases, the problem can become ill-posed, leading to numerical difficulties in finding a solution. The first notebook of this chapter illustrates the case of isotropic materials with $\frac{\lambda}{\mu} \leq 10^4$ and the second notebook is trying to illustrate the case of isotropic materials with $\frac{\lambda}{\mu} > 10^4$.
+
+<p align="center">
+
+| Material  | $\lambda$ (GPa) | $\mu$ (GPa) |
+|-----------|:--------------:|:-----------:|
+| Steel     | 120            | 80          |
+| Concrete  | 17             | 14          |
+| Rubber    | 0.16           | 0.00033     |
+
+</p>
diff --git a/chapter2/linear-elasticity-v1.ipynb b/chapter2/linear-elasticity-v1.ipynb
new file mode 100644
index 0000000..4fda0a8
--- /dev/null
+++ b/chapter2/linear-elasticity-v1.ipynb
@@ -0,0 +1,372 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "9f28f9af",
+   "metadata": {},
+   "source": [
+    "# Linear Elasticity Problem with Non-mixed Element Method\n",
+    "\n",
+    "## The Variational Formulation (Pure Displacement Formulation)\n",
+    "This section describes the variational formulation of the linear elasticity problem in a pure displacement formulation. The problem is defined on a domain $\\Omega$ with boundary $\\partial\\Omega$, where we seek to find a displacement field $u$ that satisfies the equilibrium equations under given boundary conditions. \n",
+    "\n",
+    "This leads to the following variational formulation:\n",
+    "\n",
+    "$$\n",
+    "\\boxed{\n",
+    "\\begin{aligned}\n",
+    "&\\text{Find } u \\in V \\text{ such that:} \\\\\n",
+    "&\\qquad a(u, v) = L(v) \\quad \\forall v \\in V\n",
+    "\\end{aligned}\n",
+    "}\n",
+    "$$\n",
+    "\n",
+    "with\n",
+    "\n",
+    "$$\n",
+    "\\begin{aligned}\n",
+    "&a : \n",
+    "\\begin{cases}\n",
+    "V \\times V \\rightarrow \\mathbb{R} \\\\\n",
+    "(u, v) \\longmapsto \\int_\\Omega \\sigma (u) : \\boldsymbol{\\varepsilon}(v) \\, dx\n",
+    "\\end{cases} \\\\[0.3cm]\n",
+    "&L : \n",
+    "\\begin{cases}\n",
+    "V \\rightarrow \\mathbb{R} \\\\\n",
+    "v \\longmapsto \\int_\\Omega f \\cdot v \\, dx + \\int_{\\partial\\Omega_T} g_T \\cdot v \\, ds\n",
+    "\\end{cases}\n",
+    "\\end{aligned}\n",
+    "$$\n",
+    "\n",
+    "\n",
+    "This variational formulation is essential for solving linear elasticity problems numerically using PSYDAC."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "5a958607",
+   "metadata": {},
+   "source": [
+    "## Formal Model"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "00b02838",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from sympde.expr import BilinearForm, LinearForm, integral\n",
+    "from sympde.expr     import find, EssentialBC, Norm, SemiNorm\n",
+    "from sympde.topology import VectorFunctionSpace, Cube, element_of, Union, NormalVector\n",
+    "from sympde.calculus import grad, inner, div, Transpose\n",
+    "from sympde.core import Constant\n",
+    "from sympde.core import Matrix\n",
+    "\n",
+    "from sympy import Identity, Tuple, sin, pi, cos\n",
+    "\n",
+    "from psydac.api.discretization import discretize\n",
+    "\n",
+    "\n",
+    "domain = Cube()\n",
+    "Gamma_DH = Union(domain.get_boundary(axis=2,ext=-1),domain.get_boundary(axis=2,ext=1))\n",
+    "Gamma_DNH = Union(domain.get_boundary(axis=0,ext=-1),domain.get_boundary(axis=0,ext=1))\n",
+    "Gamma_N = domain.boundary.complement(Union(Gamma_DH, Gamma_DNH))\n",
+    "\n",
+    "nn = NormalVector(domain.boundary)\n",
+    "\n",
+    "\n",
+    "lambda_ = Constant('lambda_', real=True)\n",
+    "mu   = Constant('mu',   real=True)\n",
+    "\n",
+    "epsilon = lambda w: 0.5*(grad(w) + Transpose(grad(w)))\n",
+    "sigma = lambda w: lambda_ * div(w) * Identity(3) + 2 * mu * epsilon(w)\n",
+    "\n",
+    "V = VectorFunctionSpace('V', domain)\n",
+    "\n",
+    "x,y,z = domain.coordinates\n",
+    "\n",
+    "u,v = [element_of(V, name=i) for i in ['u', 'v']]\n",
+    "\n",
+    "\n",
+    "# # bilinear form\n",
+    "a = BilinearForm((u,v), integral(domain , inner(sigma(u), epsilon(v))))\n",
+    "\n",
+    "mu_val=1.\n",
+    "lambda_val=1.25\n",
+    "\n",
+    "# Exact solution\n",
+    "ue1 = 0\n",
+    "ue2 = 0\n",
+    "ue3 = cos(pi*x)*cos(pi*y)*sin(pi*z)\n",
+    "ue = Tuple(ue1, ue2, ue3)\n",
+    "\n",
+    "#Source term\n",
+    "f1 = pi*pi * (lambda_ + mu) *   sin(pi*x) * cos(pi*y) * cos(pi*z)\n",
+    "f2 = pi*pi * (lambda_ + mu) *   cos(pi*x) * sin(pi*y) * cos(pi*z)\n",
+    "f3 = pi*pi * (4*mu + lambda_) * cos(pi*x) * cos(pi*y) * sin(pi*z)\n",
+    "f = Tuple(f1, f2, f3)\n",
+    "t = Tuple(0, -pi * lambda_ * cos(pi*x) * cos(pi*z), 0)\n",
+    "gi = ue\n",
+    "\n",
+    "# linear forms\n",
+    "l = LinearForm(v, integral(domain, inner(f, v)))\n",
+    "ln = LinearForm(v, integral(Gamma_N, inner(t, v)))\n",
+    "\n",
+    "# Dirichlet boundary conditions\n",
+    "bc = [EssentialBC(u, 0, Gamma_DH)]\n",
+    "bc += [EssentialBC(u, gi, Gamma_DNH)]  \n",
+    "\n",
+    "# Variational problem\n",
+    "equation = find(u, forall=v, lhs=a(u, v), rhs=l(v)+ln(v), bc=bc)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "4f983ece",
+   "metadata": {},
+   "source": [
+    "## Discretization"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "51095918",
+   "metadata": {},
+   "source": [
+    "We shall need the **discretize** function from **PsyDAC**."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "cedaf5b0",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "d = 2 # discretization parameters\n",
+    "degree = [d,d,d]\n",
+    "h = 6\n",
+    "ncells = [h,h,h]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "801069cc",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Create computational domain from topological domain\n",
+    "domain_h = discretize(domain, ncells=ncells, comm=None)\n",
+    "\n",
+    "# Create discrete spline space\n",
+    "Vh = discretize(V, domain_h, degree=degree)\n",
+    "\n",
+    "# Discretize equation\n",
+    "equation_h = discretize(equation, domain_h, [Vh, Vh])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "7b29fbcf",
+   "metadata": {},
+   "source": [
+    "## Solving the PDE"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3d2c99bc",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "equation_h.set_solver('gmres', info=False, tol=1e-8)\n",
+    "uh = equation_h.solve(mu=mu_val, lambda_=lambda_val)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "685a113f",
+   "metadata": {},
+   "source": [
+    "## Computing the $L^2$ norm"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ec0c3c42",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "u = element_of(V, name='u')\n",
+    "\n",
+    "error = Matrix([u[0]-ue[0], u[1]-ue[1], u[2]-ue[2]])\n",
+    "\n",
+    "# create the formal Norm object\n",
+    "l2norm = Norm(error, domain, kind='l2')\n",
+    "\n",
+    "# discretize the norm\n",
+    "l2norm_h = discretize(l2norm, domain_h, Vh)\n",
+    "\n",
+    "# assemble the norm\n",
+    "l2_error = l2norm_h.assemble(u=uh)\n",
+    "\n",
+    "print(l2_error)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "e1854531",
+   "metadata": {},
+   "source": [
+    "## Computing the $H^1$ semi-norm"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8a2181b0",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# create the formal Norm object\n",
+    "h1_seminorm = SemiNorm(error, domain, kind='h1')\n",
+    "\n",
+    "# discretize the norm\n",
+    "h1_seminorm_h = discretize(h1_seminorm, domain_h, Vh)\n",
+    "\n",
+    "# assemble the norm\n",
+    "h1_semierror = h1_seminorm_h.assemble(u=uh)\n",
+    "\n",
+    "print(h1_semierror)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "3fc201db",
+   "metadata": {},
+   "source": [
+    "## Computing the $H^1$ norm"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4f6e01df",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# create the formal Norm object\n",
+    "h1norm = Norm(error, domain, kind='h1')\n",
+    "\n",
+    "# discretize the norm\n",
+    "h1norm_h = discretize(h1norm, domain_h, Vh)\n",
+    "\n",
+    "# assemble the norm\n",
+    "h1_error = h1norm_h.assemble(u=uh)\n",
+    "\n",
+    "print(h1_error)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "38d7101f",
+   "metadata": {},
+   "source": [
+    "## Plotting the simulation results, exact solution and error"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "db432cfd",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# %% plotting exact solution, simulation solution and error\n",
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "\n",
+    "h = 100\n",
+    "# Create a grid for plotting\n",
+    "x_vals = np.linspace(0, 1, h)\n",
+    "y_vals = np.linspace(0, 1, h)\n",
+    "z_vals = np.linspace(0, 1, h)\n",
+    "X, Y, Z = np.meshgrid(x_vals, y_vals, z_vals)\n",
+    "x_lim = 1\n",
+    "y_lim = 1\n",
+    "z_lim = 1\n",
+    "\n",
+    "# Calculate the exact solution on the grid\n",
+    "exact_u1 = np.zeros_like(X)\n",
+    "exact_u2 = np.zeros_like(Y)\n",
+    "exact_u3 = np.cos(np.pi * X) * np.cos(np.pi * Y) * np.sin(np.pi * Z)\n",
+    "\n",
+    "# Calculate simulated solution on the grid\n",
+    "u_simulated = Vh.eval_field(uh, x_vals, y_vals, z_vals)\n",
+    "u1 = u_simulated[0]\n",
+    "u2 = u_simulated[1]\n",
+    "u3 = u_simulated[2]\n",
+    "\n",
+    "# Calculate the error\n",
+    "error_u1 = u1 - exact_u1\n",
+    "error_u2 = u2 - exact_u2\n",
+    "error_u3 = u3 - exact_u3\n",
+    "\n",
+    "\n",
+    "\n",
+    "# %%\n",
+    "z_plane_list = np.linspace(0, 1, 5)  # List of z-planes to plot\n",
+    "fig = plt.figure(figsize=(24, 6 * len(z_plane_list)))\n",
+    "\n",
+    "for i, z_plane in enumerate(z_plane_list):\n",
+    "    list_index = int(z_plane * (h - 1))  # Index for the z-plane\n",
+    "\n",
+    "    # Simulated u3\n",
+    "    ax1 = fig.add_subplot(len(z_plane_list), 3, 3 * i + 1, projection='3d')\n",
+    "    ax1.plot_surface(X[:, :, list_index], Y[:, :, list_index], u3[:, :, list_index], cmap='viridis')\n",
+    "    ax1.set_title('Simulated u3 on z={:.2f} plane'.format(z_plane))\n",
+    "    ax1.set_xlabel('X-axis')\n",
+    "    ax1.set_ylabel('Y-axis')\n",
+    "    ax1.set_zlabel('u3_simulated')\n",
+    "    ax1.set_xlim(0, x_lim)\n",
+    "    ax1.set_ylim(0, y_lim)\n",
+    "    ax1.set_zlim(-z_lim, z_lim)\n",
+    "\n",
+    "    # Exact u3\n",
+    "    ax2 = fig.add_subplot(len(z_plane_list), 3, 3 * i + 2, projection='3d')\n",
+    "    ax2.plot_surface(X[:, :, list_index], Y[:, :, list_index], exact_u3[:, :, list_index], cmap='viridis')\n",
+    "    ax2.set_title('Exact u3 on z={:.2f} plane'.format(z_plane))\n",
+    "    ax2.set_xlabel('X-axis')\n",
+    "    ax2.set_ylabel('Y-axis')\n",
+    "    ax2.set_zlabel('u3_exact')\n",
+    "    ax2.set_xlim(0, x_lim)\n",
+    "    ax2.set_ylim(0, y_lim)\n",
+    "    ax2.set_zlim(-z_lim, z_lim)\n",
+    "\n",
+    "    # Error in u3\n",
+    "    ax3 = fig.add_subplot(len(z_plane_list), 3, 3 * i + 3, projection='3d')\n",
+    "    ax3.plot_surface(X[:, :, list_index], Y[:, :, list_index], error_u3[:, :, list_index], cmap='viridis')\n",
+    "    ax3.set_title('Error in u3 on z={:.2f} plane'.format(z_plane))\n",
+    "    ax3.set_xlabel('X-axis')\n",
+    "    ax3.set_ylabel('Y-axis')\n",
+    "    ax3.set_zlabel('Error in u3')\n",
+    "    ax3.set_xlim(0, x_lim)\n",
+    "    ax3.set_ylim(0, y_lim)\n",
+    "    ax3.set_zlim(error_u3[:, :, list_index].min(), error_u3[:, :, list_index].max())\n",
+    "    mappable = ax3.plot_surface(X[:, :, list_index], Y[:, :, list_index], error_u3[:, :, list_index], cmap='viridis')\n",
+    "    fig.colorbar(mappable, ax=ax3, shrink=0.5, aspect=5)\n",
+    "\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  }
+ ],
+ "metadata": {},
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/chapter2/linear-elasticity-v2.ipynb b/chapter2/linear-elasticity-v2.ipynb
new file mode 100644
index 0000000..5b7e7ad
--- /dev/null
+++ b/chapter2/linear-elasticity-v2.ipynb
@@ -0,0 +1,510 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "9f28f9af",
+   "metadata": {},
+   "source": [
+    "# Linear Elasticity Problem with Mixed Element Method\n",
+    "This notebook treats the linear elasticity problem for materials with high values of the Lamé parameters ($\\frac{\\lambda}{\\mu} \\gg 1$).\n",
+    "\n",
+    "\n",
+    "## The pressure field\n",
+    "In the case of high values of the Lamé parameter $ \\lambda $, the pressure field $ p $ can be defined as:\n",
+    "\n",
+    "$$\n",
+    "p := -\\lambda \\nabla \\cdot u \\hspace{1cm} \\text{ in } \\Omega\n",
+    "$$\n",
+    "\n",
+    "As $u \\in H^1(\\Omega)$, the pressure field $p \\in L^2(\\Omega)$.\n",
+    "\n",
+    "## The strong formulation of the problem\n",
+    "The strong formulation of the problem can be expressed as follows:\n",
+    "Find $(u, p) \\in V \\times L^2(\\Omega)$ such that\n",
+    "\n",
+    "$$\n",
+    "\\begin{align}\n",
+    "    -\\nabla \\cdot \\sigma(u,p) &= f & \\text{in } & \\Omega \\\\\n",
+    "    u &= 0 & \\text{on } & \\partial \\Omega_D \\\\\n",
+    "    \\sigma(u,p) \\cdot n &= g_T & \\text{on } & \\partial \\Omega_T \\\\\n",
+    "    \\nabla \\cdot u + \\frac{1}{\\lambda} p &= 0 & \\text{in } & \\Omega\n",
+    "\\end{align}\n",
+    "$$\n",
+    "\n",
+    "## The Variational Formulation (Pressure-Displacement Formulation)\n",
+    "To derive the weak formulation of the problem, we multiply the first equation by a test function $v \\in V$ and integrate by parts over the domain $\\Omega$ and we add the fourth equation multiplied by a test function $q \\in L^2(\\Omega)$ integrated over the domain $\\Omega$:\n",
+    "\n",
+    "$$\n",
+    "\\boxed{\n",
+    "\\begin{aligned}\n",
+    "&\\text{Find } (u, p) \\in V \\times L^2(\\Omega) \\text{ such that:} \\\\\n",
+    "&\\qquad \\tilde{a}((u, p), (v, q)) = l(v, q) \\quad \\forall (v, q) \\in V \\times L^2(\\Omega)\n",
+    "\\end{aligned}\n",
+    "}\n",
+    "$$\n",
+    "\n",
+    "with\n",
+    "\n",
+    "$$\n",
+    "\\begin{aligned}\n",
+    "&\\tilde{a} : \n",
+    "\\begin{cases}\n",
+    "\\left( V \\times L^2(\\Omega) \\right)^2 \\rightarrow \\mathbb{R} \\\\\n",
+    "((u, p), (v, q)) \\longmapsto \\int_\\Omega \\left( 2\\mu\\, \\boldsymbol{\\varepsilon}(u) : \\boldsymbol{\\varepsilon}(v) - p\\, \\mathrm{div}\\, v + (\\mathrm{div}\\, u)\\, q + \\frac{1}{\\lambda} p q \\right) \\, dx\n",
+    "\\end{cases} \\\\[0.3cm]\n",
+    "&l : \n",
+    "\\begin{cases}\n",
+    "V \\rightarrow \\mathbb{R} \\\\\n",
+    "v \\longmapsto \\int_\\Omega f \\cdot v \\, dx + \\int_{\\partial\\Omega_N} t_N \\cdot v \\, ds\n",
+    "\\end{cases}\n",
+    "\\end{aligned}\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "5a958607",
+   "metadata": {},
+   "source": [
+    "## Formal Model"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "00b02838",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from sympde.expr import BilinearForm, LinearForm, integral\n",
+    "from sympde.expr     import find, EssentialBC, Norm, SemiNorm\n",
+    "from sympde.topology import VectorFunctionSpace, ScalarFunctionSpace, Cube, element_of\n",
+    "from sympde.calculus import grad, inner, div, Transpose\n",
+    "from sympde.core import Constant\n",
+    "\n",
+    "from sympy import Tuple, sin, pi, cos\n",
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "from psydac.api.discretization import discretize\n",
+    "from psydac.fem.basic          import FemField\n",
+    "\n",
+    "domain = Cube()\n",
+    "\n",
+    "\n",
+    "lambda_ = Constant('lambda_', real=True)\n",
+    "mu    = Constant('mu',    real=True)\n",
+    "\n",
+    "V1 = VectorFunctionSpace('V1', domain,kind='H1')\n",
+    "V2 = ScalarFunctionSpace('V2', domain,kind='L2')\n",
+    "X = V1*V2\n",
+    "\n",
+    "x,y,z = domain.coordinates\n",
+    "\n",
+    "u,v = [element_of(V1, name=i) for i in ['u', 'v']]\n",
+    "p,q = [element_of(V2, name=i) for i in ['p', 'q']]\n",
+    "\n",
+    "epsilon = lambda w: 0.5*(grad(w) + Transpose(grad(w)))\n",
+    "\n",
+    "#bilinear form\n",
+    "a = BilinearForm(((u,p),(v,q)), integral(domain, inner(2*mu*epsilon(u),epsilon(v)) - p*div(v) + div(u)*q + (1/lambda_)*p*q))\n",
+    "\n",
+    "mu_val =1.\n",
+    "lambda_val =1.\n",
+    "\n",
+    "# Exact solution\n",
+    "ue1 = 0\n",
+    "ue2 = 0\n",
+    "ue3 = sin(pi*x)*sin(pi*y)*sin(pi*z)\n",
+    "ue = Tuple(ue1, ue2, ue3)\n",
+    "pe = -lambda_ * pi * sin(pi*x) * sin(pi*y) * cos(pi*z)\n",
+    "\n",
+    "#Second member\n",
+    "f1 = - pi*pi * (lambda_ + mu) *  cos(pi*x) * sin(pi*y) * cos(pi*z)\n",
+    "f2 = - pi*pi * (lambda_ + mu) *  sin(pi*x) * cos(pi*y) * cos(pi*z)\n",
+    "f3 =   pi*pi * (4*mu + lambda_) * sin(pi*x) * sin(pi*y) * sin(pi*z)\n",
+    "f = Tuple(f1, f2, f3)\n",
+    "\n",
+    "# linear form\n",
+    "l = LinearForm((v,q), integral(domain, inner(f, v)))\n",
+    "\n",
+    "# Essential boundary conditions\n",
+    "bc = [EssentialBC(u, 0, domain.boundary)]\n",
+    "\n",
+    "equation = find((u,p), forall=(v,q), lhs=a((u,p),(v,q)), rhs=l(v,q), bc=bc)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "4f983ece",
+   "metadata": {},
+   "source": [
+    "## Discretization"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "51095918",
+   "metadata": {},
+   "source": [
+    "We shall need the **discretize** function from **PsyDAC**."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "cedaf5b0",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "d = 2 # discretization parameters\n",
+    "degree = [d,d,d]\n",
+    "h = 6\n",
+    "ncells = [h,h,h]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "801069cc",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Create computational domain from topological domain\n",
+    "domain_h = discretize(domain, ncells=ncells, comm=None)\n",
+    "\n",
+    "# Create discrete spline space\n",
+    "V1h = discretize(V1, domain_h, degree=degree)\n",
+    "V2h = discretize(V2, domain_h, degree=degree)\n",
+    "Xh = discretize(X,   domain_h, degree=degree)\n",
+    "\n",
+    "# Discretize equation\n",
+    "equation_h = discretize(equation, domain_h, [Xh, Xh])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "7b29fbcf",
+   "metadata": {},
+   "source": [
+    "## Solving the PDE"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3d2c99bc",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "equation_h.set_solver('gmres', info=False, tol=1e-8)\n",
+    "phi_h = equation_h.solve(mu=mu_val, lambda_=lambda_val)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "685a113f",
+   "metadata": {},
+   "source": [
+    "## Computing the $L^2$ norm on $p$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ec0c3c42",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ph = FemField(V2h)\n",
+    "ph.coeffs[:] = phi_h.coeffs[3][:]\n",
+    "\n",
+    "# create the formal Norm object\n",
+    "l2norm_p  = Norm(pe - p, domain, kind='l2')\n",
+    "\n",
+    "# discretize the Norm object\n",
+    "l2norm_ph = discretize(l2norm_p, domain_h, V2h)\n",
+    "\n",
+    "# assemble the norm\n",
+    "l2_error_p = l2norm_ph.assemble(p=ph, lambda_=lambda_val, mu=mu_val)\n",
+    "\n",
+    "print(\"L2 error in pressure:\", l2_error_p)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "b46cb158",
+   "metadata": {},
+   "source": [
+    "## Computing the $L^2$ norm on $u$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "bf324fe3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "uh = FemField(V1h)\n",
+    "uh.coeffs[0][:] = phi_h.coeffs[0][:]\n",
+    "uh.coeffs[1][:] = phi_h.coeffs[1][:]\n",
+    "uh.coeffs[2][:] = phi_h.coeffs[2][:]\n",
+    "\n",
+    "# create the formal Norm object\n",
+    "error_u   = [ue[0]-u[0], ue[1]-u[1], ue[2]-u[2]]\n",
+    "l2norm_u = Norm(error_u, domain, kind='l2')\n",
+    "\n",
+    "# discretize the Norm object\n",
+    "l2norm_uh = discretize(l2norm_u, domain_h, V1h)\n",
+    "\n",
+    "# assemble the Norm object\n",
+    "l2_error_u = l2norm_uh.assemble(u=uh, lambda_=lambda_val, mu=mu_val)\n",
+    "\n",
+    "print(\"L2 error in displacement:\", l2_error_u)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "e1854531",
+   "metadata": {},
+   "source": [
+    "## Computing the $H^1$ semi-norm on $u$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8a2181b0",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# create the formal Norm object\n",
+    "h1_seminorm = SemiNorm(error_u, domain, kind='h1')\n",
+    "\n",
+    "# discretize the norm\n",
+    "h1_seminorm_h = discretize(h1_seminorm, domain_h, V1h)\n",
+    "\n",
+    "# assemble the norm\n",
+    "h1_semierror = h1_seminorm_h.assemble(u=uh)\n",
+    "\n",
+    "print(h1_semierror)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "3fc201db",
+   "metadata": {},
+   "source": [
+    "## Computing the $H^1$ norm on $u$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4f6e01df",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# create the formal Norm object\n",
+    "h1norm = Norm(error_u, domain, kind='h1')\n",
+    "\n",
+    "# discretize the norm\n",
+    "h1norm_h = discretize(h1norm, domain_h, V1h)\n",
+    "\n",
+    "# assemble the norm\n",
+    "h1_error = h1norm_h.assemble(u=uh)\n",
+    "\n",
+    "print(h1_error)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "38d7101f",
+   "metadata": {},
+   "source": [
+    "## Plotting the simulation results, exact solution and error"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "db432cfd",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# %% plotting exact solution, simulation solution and error\n",
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "\n",
+    "h = 100\n",
+    "# Create a grid for plotting\n",
+    "x_vals = np.linspace(0, 1, h)\n",
+    "y_vals = np.linspace(0, 1, h)\n",
+    "z_vals = np.linspace(0, 1, h)\n",
+    "X, Y, Z = np.meshgrid(x_vals, y_vals, z_vals)\n",
+    "x_lim = 1\n",
+    "y_lim = 1\n",
+    "z_lim = 1\n",
+    "\n",
+    "# Calculate the exact solution on the grid\n",
+    "exact_u1 = np.zeros_like(X)\n",
+    "exact_u2 = np.zeros_like(Y)\n",
+    "exact_u3 = np.sin(np.pi * X) * np.sin(np.pi * Y) * np.sin(np.pi * Z)\n",
+    "exact_p = -lambda_val * np.pi * np.sin(np.pi * X) * np.sin(np.pi * Y) * np.cos(np.pi * Z)\n",
+    "\n",
+    "uh_sim = Xh.eval_field(phi_h, x_vals, y_vals, z_vals)\n",
+    "u1 = uh_sim[0]\n",
+    "u2 = uh_sim[1]\n",
+    "u3 = uh_sim[2]\n",
+    "\n",
+    "# Calculate the error\n",
+    "error_u1 = u1 - exact_u1\n",
+    "error_u2 = u2 - exact_u2\n",
+    "error_u3 = u3 - exact_u3\n",
+    "\n",
+    "# Displacement plots\n",
+    "z_plane_list = np.linspace(0, 1, 5)  #List of z-planes to plot\n",
+    "fig = plt.figure(figsize=(24, 6 * len(z_plane_list)))\n",
+    "\n",
+    "for i, z_plane in enumerate(z_plane_list):\n",
+    "    list_index = int(z_plane * (h - 1))  # Index for the z-plane\n",
+    "\n",
+    "    # Simulated u3\n",
+    "    ax1 = fig.add_subplot(len(z_plane_list), 3, 3 * i + 1, projection='3d')\n",
+    "    ax1.plot_surface(X[:, :, list_index], Y[:, :, list_index], u3[:, :, list_index], cmap='viridis')\n",
+    "    ax1.set_title('Simulated u3 on z={:.2f} plane'.format(z_plane))\n",
+    "    ax1.set_xlabel('X-axis')\n",
+    "    ax1.set_ylabel('Y-axis')\n",
+    "    ax1.set_zlabel('u3_simulated')\n",
+    "    ax1.set_xlim(0, x_lim)\n",
+    "    ax1.set_ylim(0, y_lim)\n",
+    "    ax1.set_zlim(-z_lim, z_lim)\n",
+    "\n",
+    "    # Exact u3\n",
+    "    ax2 = fig.add_subplot(len(z_plane_list), 3, 3 * i + 2, projection='3d')\n",
+    "    ax2.plot_surface(X[:, :, list_index], Y[:, :, list_index], exact_u3[:, :, list_index], cmap='viridis')\n",
+    "    ax2.set_title('Exact u3 on z={:.2f} plane'.format(z_plane))\n",
+    "    ax2.set_xlabel('X-axis')\n",
+    "    ax2.set_ylabel('Y-axis')\n",
+    "    ax2.set_zlabel('u3_exact')\n",
+    "    ax2.set_xlim(0, x_lim)\n",
+    "    ax2.set_ylim(0, y_lim)\n",
+    "    ax2.set_zlim(-z_lim, z_lim)\n",
+    "\n",
+    "    # Error in u3\n",
+    "    ax3 = fig.add_subplot(len(z_plane_list), 3, 3 * i + 3, projection='3d')\n",
+    "    ax3.plot_surface(X[:, :, list_index], Y[:, :, list_index], error_u3[:, :, list_index], cmap='viridis')\n",
+    "    ax3.set_title('Error in u3 on z={:.2f} plane'.format(z_plane))\n",
+    "    ax3.set_xlabel('X-axis')\n",
+    "    ax3.set_ylabel('Y-axis')\n",
+    "    ax3.set_zlabel('Error in u3')\n",
+    "    ax3.set_xlim(0, x_lim)\n",
+    "    ax3.set_ylim(0, y_lim)\n",
+    "    ax3.set_zlim(error_u3[:, :, list_index].min(), error_u3[:, :, list_index].max())\n",
+    "    mappable = ax3.plot_surface(X[:, :, list_index], Y[:, :, list_index], error_u3[:, :, list_index], cmap='viridis')\n",
+    "    fig.colorbar(mappable, ax=ax3, shrink=0.5, aspect=5)\n",
+    "\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6226718f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def compute_divergence(u1, u2, u3, x, y, z):\n",
+    "\tdu1_dx = np.gradient(u1, x, axis=0)\n",
+    "\tdu2_dy = np.gradient(u2, y, axis=1)\n",
+    "\tdu3_dz = np.gradient(u3, z, axis=2)\n",
+    "\treturn du1_dx + du2_dy + du3_dz\n",
+    "\n",
+    "p_sim = - lambda_val * compute_divergence(u1, u2, u3, x_vals, y_vals, z_vals)\n",
+    "\n",
+    "# Pressure plots\n",
+    "fig, ax = plt.subplots(1, 3, figsize=(18, 6))\n",
+    "# Simulated pressure\n",
+    "c1 = ax[0].contourf(X[:, :, list_index], Y[:, :, list_index], p_sim[:, :, list_index], cmap='viridis')\n",
+    "ax[0].set_title('Simulated Pressure on z={:.2f} plane'.format(z_plane))\n",
+    "ax[0].set_xlabel('X-axis')\n",
+    "ax[0].set_ylabel('Y-axis')\n",
+    "ax[0].set_xlim(0, x_lim)\n",
+    "ax[0].set_ylim(0, y_lim)\n",
+    "fig.colorbar(c1, ax=ax[0], shrink=0.5, aspect=5)\n",
+    "\n",
+    "# Exact pressure\n",
+    "c2 = ax[1].contourf(X[:, :, list_index], Y[:, :, list_index], exact_p[:, :, list_index], cmap='viridis')\n",
+    "ax[1].set_title('Exact Pressure on z={:.2f} plane'.format(z_plane))\n",
+    "ax[1].set_xlabel('X-axis')\n",
+    "ax[1].set_ylabel('Y-axis')\n",
+    "ax[1].set_xlim(0, x_lim)\n",
+    "ax[1].set_ylim(0, y_lim)\n",
+    "fig.colorbar(c2, ax=ax[1], shrink=0.5, aspect=5)\n",
+    "\n",
+    "# Error in pressure\n",
+    "c3 = ax[2].contourf(X[:, :, list_index], Y[:, :, list_index], p_sim[:, :, list_index] - exact_p[:, :, list_index], cmap='viridis')\n",
+    "ax[2].set_title('Error in Pressure on z={:.2f} plane'.format(z_plane))\n",
+    "ax[2].set_xlabel('X-axis')\n",
+    "ax[2].set_ylabel('Y-axis')\n",
+    "ax[2].set_xlim(0, x_lim)\n",
+    "ax[2].set_ylim(0, y_lim)\n",
+    "fig.colorbar(c3, ax=ax[2], shrink=0.5, aspect=5)\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6bdd2964",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig = plt.figure(figsize=(24, 6 * len(z_plane_list)))\n",
+    "for i, z_plane in enumerate(z_plane_list):\n",
+    "    list_index = int(z_plane * (h - 1))\n",
+    "\n",
+    "    # Simulated pressure\n",
+    "    ax1 = fig.add_subplot(len(z_plane_list), 3, 3 * i + 1, projection='3d')\n",
+    "    ax1.plot_surface(X[:, :, list_index], Y[:, :, list_index], p_sim[:, :, list_index], cmap='viridis')\n",
+    "    ax1.set_title('Simulated Pressure on z={:.2f} plane'.format(z_plane))\n",
+    "    ax1.set_xlabel('X-axis')\n",
+    "    ax1.set_ylabel('Y-axis')\n",
+    "    ax1.set_zlabel('p_sim')\n",
+    "    ax1.set_xlim(0, x_lim)\n",
+    "    ax1.set_ylim(0, y_lim)\n",
+    "    ax1.set_zlim(p_sim[:, :, list_index].min(), p_sim[:, :, list_index].max())\n",
+    "\n",
+    "    # Exact pressure\n",
+    "    ax2 = fig.add_subplot(len(z_plane_list), 3, 3 * i + 2, projection='3d')\n",
+    "    ax2.plot_surface(X[:, :, list_index], Y[:, :, list_index], exact_p[:, :, list_index], cmap='viridis')\n",
+    "    ax2.set_title('Exact Pressure on z={:.2f} plane'.format(z_plane))\n",
+    "    ax2.set_xlabel('X-axis')\n",
+    "    ax2.set_ylabel('Y-axis')\n",
+    "    ax2.set_zlabel('p_exact')\n",
+    "    ax2.set_xlim(0, x_lim)\n",
+    "    ax2.set_ylim(0, y_lim)\n",
+    "    ax2.set_zlim(exact_p[:, :, list_index].min(), exact_p[:, :, list_index].max())\n",
+    "\n",
+    "    # Error in pressure\n",
+    "    ax3 = fig.add_subplot(len(z_plane_list), 3, 3 * i + 3, projection='3d')\n",
+    "    mappable = ax3.plot_surface(X[:, :, list_index], Y[:, :, list_index], p_sim[:, :, list_index] - exact_p[:, :, list_index], cmap='viridis')\n",
+    "    ax3.set_title('Error in Pressure on z={:.2f} plane'.format(z_plane))\n",
+    "    ax3.set_xlabel('X-axis')\n",
+    "    ax3.set_ylabel('Y-axis')\n",
+    "    ax3.set_zlabel('Error in pressure')\n",
+    "    ax3.set_xlim(0, x_lim)\n",
+    "    ax3.set_ylim(0, y_lim)\n",
+    "    ax3.set_zlim((p_sim[:, :, list_index] - exact_p[:, :, list_index]).min(), (p_sim[:, :, list_index] - exact_p[:, :, list_index]).max())\n",
+    "    fig.colorbar(mappable, ax=ax3, shrink=0.5, aspect=5)\n",
+    "\n",
+    "plt.tight_layout()\n",
+    "plt.show()\n"
+   ]
+  }
+ ],
+ "metadata": {},
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/chapter2/linear-elasticity.ipynb b/chapter2/linear-elasticity.ipynb
deleted file mode 100644
index f4d1366..0000000
--- a/chapter2/linear-elasticity.ipynb
+++ /dev/null
@@ -1,225 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "markdown",
-   "id": "9f28f9af",
-   "metadata": {},
-   "source": [
-    "# Linear Elasticity Problem\n",
-    "\n",
-    "Analysis of deformable structures is essential in engineering, with the equations of linear elasticity being fundamental in this field. In this section, we present the variational formulation of linear elasticity equations using the principle of virtual work.\n",
-    "\n",
-    "## The PDE problem\n",
-    "The governing equations for small elastic deformations of a body $ \\Omega $ can be expressed as:\n",
-    "\n",
-    "$$\n",
-    "\\begin{align}\n",
-    "    -\\nabla \\cdot \\sigma(u) &= f & \\text{in } \\Omega \\\\\n",
-    "    \\sigma(u) &= \\kappa \\text{tr}(\\epsilon(u))I + 2 \\mu \\epsilon(u)\n",
-    "\\end{align}\n",
-    "$$\n",
-    "\n",
-    "where $ \\sigma $ is the stress tensor, $ f $ represents the body force per unit volume, $ \\kappa $ and $ \\mu $ are Lamé's elasticity parameters for the material, $ I $ denotes the identity tensor, and $ \\epsilon $ is the symmetric strain tensor. The displacement vector field is denoted by $ u $.\n",
-    "\n",
-    "By substituting $ \\epsilon(u) $ into $ \\sigma $, we obtain:\n",
-    "\n",
-    "$$\n",
-    "\\sigma(u) = \\kappa (\\nabla \\cdot u)I + \\mu(\\nabla u + (\\nabla u)^T)\n",
-    "$$\n",
-    "\n",
-    "## The Variational Formulation\n",
-    "The variational formulation of the linear elasticity equations involves forming the inner product of the PDE with a vector test function $ v \\in \\hat{V} $ and integrating over the domain $ \\Omega $. This yields:\n",
-    "\n",
-    "$$\n",
-    "\\int_{\\Omega} \\sigma : \\nabla v \\, \\mathrm{d} x = \\int_{\\Omega} f \\cdot v \\, \\mathrm{d} x\n",
-    "$$\n",
-    "\n",
-    "Integrating the term $ \\nabla \\cdot \\sigma \\cdot v $ by parts, considering boundary conditions, we obtain:\n",
-    "\n",
-    "$$\n",
-    "\\int_{\\Omega} \\sigma : \\nabla v \\, \\mathrm{d} x = \\int_{\\Omega} f \\cdot v \\, \\mathrm{d} x + \\int_{\\partial \\Omega_T} T \\cdot v \\, \\mathrm{d} s\n",
-    "$$\n",
-    "\n",
-    "where $ T $ represents the traction vector on the part $ \\partial \\Omega_T $ of the boundary where it's prescribed.\n",
-    "\n",
-    "This leads to the variational formulation: Find $ u \\in V $ such that\n",
-    "\n",
-    "$$\n",
-    "a(u, v) = L(v) \\quad \\forall v \\in \\hat{V}\n",
-    "$$\n",
-    "\n",
-    "where\n",
-    "\n",
-    "$$\n",
-    "\\begin{align}\n",
-    "    a(u, v) &= \\int_{\\Omega} \\sigma(u) : \\nabla v \\, \\mathrm{d} x \\\\\n",
-    "    L(v) &= \\int_{\\Omega} f \\cdot v \\, \\mathrm{d} x + \\int_{\\partial \\Omega_T} T \\cdot v \\, \\mathrm{d} s\n",
-    "\\end{align}\n",
-    "$$\n",
-    "\n",
-    "This formulation can be alternatively expressed as:\n",
-    "\n",
-    "$$\n",
-    "a(u, v) = \\int_{\\Omega} \\sigma(u) : \\epsilon(v) \\, \\mathrm{d} x\n",
-    "$$\n",
-    "\n",
-    "where $ \\epsilon(v) = \\frac{1}{2} (\\nabla v + (\\nabla v)^T) $ is the symmetric strain tensor.\n",
-    "\n",
-    "This variational formulation is essential for solving linear elasticity problems numerically using methods like the finite element method (FEM)."
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "5a958607",
-   "metadata": {},
-   "source": [
-    "## Formal Model"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "d742586c",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "from sympde.expr import BilinearForm, LinearForm, integral\n",
-    "from sympde.expr     import find, EssentialBC, Norm, SemiNorm\n",
-    "from sympde.topology import VectorFunctionSpace, Cube, element_of\n",
-    "from sympde.calculus import grad, dot, inner, outer, cross, div\n",
-    "from sympde.core import Constant\n",
-    "from sympde.core import Matrix, Vector, Transpose\n",
-    "\n",
-    "domain = Cube()\n",
-    "\n",
-    "I = Matrix([[1, 0, 0], [0, 1, 0], [0, 0, 1]], name='I')\n",
-    "\n",
-    "kappa = Constant('kappa', is_real=True)\n",
-    "mu    = Constant('mu',    is_real=True)\n",
-    "rho   = Constant('rho',   is_real=True)\n",
-    "\n",
-    "epsilon = lambda w: 0.5*(grad(w) + Transpose(grad(w)))\n",
-    "sigma   = lambda w: kappa * div(w) * I + 2 * mu * epsilon(w)\n",
-    "\n",
-    "V = VectorFunctionSpace('V', domain)\n",
-    "\n",
-    "x,y,z = domain.coordinates\n",
-    "\n",
-    "u,v = [element_of(V, name=i) for i in ['u', 'v']]"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "b16cc745",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# bilinear form\n",
-    "a = BilinearForm((u,v), integral(domain , inner(sigma(u), epsilon(v))))"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "0a0f152c",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "L = 1\n",
-    "W = 1 #0.2\n",
-    "delta = W / L\n",
-    "\n",
-    "#mu = 1\n",
-    "#rho = 1\n",
-    "#kappa = 1.25\n",
-    "g = 0.4 * delta**2\n",
-    "\n",
-    "# linear form\n",
-    "f = Vector([0, 0, -rho*g], name='f')\n",
-    "l = LinearForm(v, integral(domain, dot(f,v)))\n",
-    "\n",
-    "# Dirichlet boundary conditions\n",
-    "bc = [EssentialBC(u, 0, domain.boundary)]\n",
-    "\n",
-    "# Variational problem\n",
-    "equation   = find(u, forall=v, lhs=a(u, v), rhs=l(v), bc=bc)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "4f983ece",
-   "metadata": {},
-   "source": [
-    "## Discretization"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "51095918",
-   "metadata": {},
-   "source": [
-    "We shall need the **discretize** function from **PsyDAC**."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "a2a0a2a1",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "from psydac.api.discretization import discretize"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "00e54163",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "degree = [2,2,2]\n",
-    "ncells = [8,8,8]"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "5999c62b",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# Create computational domain from topological domain\n",
-    "domain_h = discretize(domain, ncells=ncells, comm=None)\n",
-    "\n",
-    "# Create discrete spline space\n",
-    "Vh = discretize(V, domain_h, degree=degree)\n",
-    "\n",
-    "# Discretize equation\n",
-    "equation_h = discretize(equation, domain_h, [Vh, Vh])"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "7b29fbcf",
-   "metadata": {},
-   "source": [
-    "## Solving the PDE"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "541192ee",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "uh = equation_h.solve(mu=1, rho=1, kappa=1.25)"
-   ]
-  }
- ],
- "metadata": {},
- "nbformat": 4,
- "nbformat_minor": 5
-}