diff --git a/Exercises_2.ipynb b/Exercises_2.ipynb
index b17326a..6f4fc31 100644
--- a/Exercises_2.ipynb
+++ b/Exercises_2.ipynb
@@ -1,557 +1,586 @@
 {
  "cells": [
   {
    "cell_type": "markdown",
    "id": "17e4b0c9-8fc8-4f4f-b41f-0b05c5ee22e5",
-   "metadata": {},
+   "metadata": {
+    "tags": []
+   },
    "source": [
     "### Setup"
    ]
   },
   {
    "cell_type": "markdown",
    "id": "b9a9a851-94d7-454c-be68-96ad35c2d0f2",
    "metadata": {},
    "source": [
-    "We start by importing some Python libraries. `numpy` to work with arrays of numbers, `matplotlib` for plotting, and `tamaas` for generating rough surfaces and solving contact problems.\n",
+    "We start by importing some Python libraries. `numpy` to work with arrays of numbers, and `matplotlib` for plotting.\n",
     "\n",
     "Clic on a code cell and press `shift+enter` to run it."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "id": "12ef8dc7-8413-4842-8938-f76a977732d8",
    "metadata": {},
    "outputs": [],
    "source": [
     "import numpy as np\n",
     "from matplotlib import pyplot as plt\n",
-    "from mpl_toolkits.mplot3d import Axes3D # to make 3d plots\n",
+    "from mpl_toolkits.mplot3d import Axes3D # to make 3d plots"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "72986cf4-f4d9-4975-8893-655bc91e93b9",
+   "metadata": {},
+   "source": [
+    "We will use `tamaas` for generating rough surfaces and solving contact problems. It must be installed before importing."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7b9e07bf-e46b-415a-818c-c899d656f317",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "!pip install --user tamaas==2.3.0.post1"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "fa04bae9-2b4f-4cc6-823a-4169a58abb31",
+   "metadata": {},
+   "outputs": [],
+   "source": [
     "import tamaas as tm\n",
     "tm.initialize()"
    ]
   },
   {
    "cell_type": "markdown",
    "id": "4fdb80fa-43ad-43e9-8f6a-947977dba487",
    "metadata": {},
    "source": [
     "# Generation of rough surfaces"
    ]
   },
   {
    "cell_type": "markdown",
    "id": "8d4ab560-e4ad-4fda-a78b-70e2de654f3d",
    "metadata": {},
    "source": [
     "Tamaas is capable of generating fractal rough surfaces, called **self-affine** surfaces. These surfaces are characterised by the shape of their power spectral density (PSD), which is linked to the spatial Fourier transform. The PSD of a rough surface indicates how present are undulations of short or large wavelenghts.\n",
     "\n",
     "The PSD of a self-affine surface is a power-law (a line in a logarithmic plot) inside a range of spatial frequencies. We can control the values of the **smallest and largest wavelengths**. The slope of the PSD is influenced by the **Hurst exponent** $H$. Finally, the heights can be scaled to match a certain magnitude, for example quantified by the **root mean square** of the height distribution.\n",
     "\n",
     "Play with the parameters below and see how they influence the final rough surface."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "id": "1960d707-de86-468c-b867-2578d9f53288",
    "metadata": {},
    "outputs": [],
    "source": [
     "# PLAY WITH THESE PARAMETERS ###################################################\n",
     "L = 1     # side length\n",
     "n = 128   # number of points (discretization)\n",
     "\n",
     "l_smallest = L / 64   # smallest wavelength (not smaller than L/n)\n",
     "l_largest  = L / 2    # largest wavelength\n",
     "h_RMS      = 0.02     # \"average\" height\n",
     "Hurst      = 0.8      # Hurst exponent (0 ≤ H ≤ 1)\n",
     "seed       = 1        # change this value to generate other surfaces with identical properties\n",
     "################################################################################"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "id": "df672cb1-776a-4e2f-8007-611def4c0bac",
    "metadata": {},
    "outputs": [],
    "source": [
     "# Create PSD object\n",
     "spectrum = tm.Isopowerlaw2D()\n",
     "\n",
     "# Set spectrum parameters\n",
     "spectrum.q0 = round(L / l_largest) # round: must be integers\n",
     "spectrum.q1 = round(L / l_largest) # q1 is the \"roll-off\" frequency, usually equal to q0\n",
     "spectrum.q2 = round(L / l_smallest)\n",
     "spectrum.hurst = Hurst\n",
     "\n",
     "# Create surface generator\n",
     "generator = tm.SurfaceGeneratorFilter2D([n, n])\n",
     "generator.spectrum = spectrum\n",
     "generator.random_seed = seed\n",
     "\n",
     "# Generate surface\n",
     "h = generator.buildSurface()\n",
     "h *= h_RMS / tm.Statistics2D.computeRMSHeights(h)\n",
     "\n",
     "\n",
     "# Plot\n",
     "x = np.linspace(-L/2, L/2, n)\n",
     "y = x.copy()\n",
     "x, y = np.meshgrid(x, y, indexing=\"ij\")\n",
     "\n",
     "fig = plt.figure(figsize=(5, 5), dpi=150)\n",
     "ax = fig.add_subplot(111, projection=\"3d\")\n",
     "surf = ax.plot_surface(x, y, h, rstride=1, cstride=1, cmap=plt.cm.inferno)\n",
     "\n",
     "z_ratio = 0.2\n",
     "ax.set_zlim([-L/2 * z_ratio, L/2 * z_ratio])\n",
     "ax.set_box_aspect([1, 1, z_ratio])\n",
     "plt.locator_params(axis=\"z\", nbins=2) # reduce number of ticks on z axis\n",
     "\n",
     "plt.show()"
    ]
   },
   {
    "cell_type": "markdown",
    "id": "ca664fdb-ee5b-4a8c-bafa-95e6771fd2cc",
    "metadata": {},
    "source": [
     "# Solving a Hertzian contact"
    ]
   },
   {
    "cell_type": "markdown",
    "id": "580921ca-ce4b-48d0-88fc-a53c7bfa1107",
    "metadata": {},
    "source": [
     "## Definition and solving"
    ]
   },
   {
    "cell_type": "markdown",
    "id": "69e95c09-0f32-425a-86ea-07bbaa8aebef",
    "metadata": {},
    "source": [
     "One simple case of Hertz contact is the contact between a sphere and a flat surface. Tamaas is made for solving the contact between a **rigid surface** having a given shape (near to flat) and a **flat elastic surface**, knowing the total applied normal load."
    ]
   },
   {
    "cell_type": "markdown",
    "id": "e0b81529-7afc-4faf-8d8b-20c3f7c60f44",
    "metadata": {},
    "source": [
     "### Geometry"
    ]
   },
   {
    "cell_type": "markdown",
    "id": "130cc0bc-f4e5-4d0b-8e61-8c7351604df4",
    "metadata": {},
    "source": [
     "Start by defining the geometry of the rigid surface, which must be a sphere."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "id": "ccba9013-0391-480a-ae39-d618b1a34454",
    "metadata": {},
    "outputs": [],
    "source": [
     "# CHANGE THESE VALUES ##########################################################\n",
     "L = 1     # side length\n",
     "n = 128   # number of points (discretization)\n",
     "R = 10    # rigid indenter radius\n",
     "################################################################################\n",
     "\n",
     "x = np.linspace(-L/2, L/2, n)\n",
     "y = x.copy()\n",
     "x, y = np.meshgrid(x, y, indexing=\"ij\")\n",
     "\n",
     "\n",
     "def create_indenter(x, y, R):\n",
     "    return # DEFINE SPHERICAL SURFACE ##########################################\n",
     "\n",
     "h = create_indenter(x, y, R)"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "id": "10b70cf7-ec28-4f45-b8fa-27387892e38f",
    "metadata": {},
    "outputs": [],
    "source": [
     "# Plot\n",
     "fig = plt.figure(figsize=(5, 5), dpi=150)\n",
     "ax = fig.add_subplot(111, projection=\"3d\")\n",
     "surf = ax.plot_surface(x, y, h, rstride=1, cstride=1, cmap=plt.cm.inferno)\n",
     "plt.show()"
    ]
   },
   {
    "cell_type": "markdown",
    "id": "17733385-28c3-4929-aee6-0df3f5d6f279",
    "metadata": {},
    "source": [
     "### Physical parameters"
    ]
   },
   {
    "cell_type": "markdown",
    "id": "2a4b9eef-1d00-40c1-a556-c8e82ef9522b",
    "metadata": {},
    "source": [
     "Then, we define the material parameters and the normal load applied by the rigid sphere."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "id": "ce342c8c-6213-45e2-a7c2-2a3c01f53e4a",
    "metadata": {},
    "outputs": [],
    "source": [
     "E  = 1.0    # Young's modulus\n",
     "nu = 0.3    # Poisson's ratio\n",
     "F  = 1e-4   # Normal load # CHANGE THIS VALUE ##################################"
    ]
   },
   {
    "cell_type": "markdown",
    "id": "7e6da2e6-6b62-48a4-852d-7f0ea8a503fd",
    "metadata": {},
    "source": [
     "### Tamaas model and solver"
    ]
   },
   {
    "cell_type": "markdown",
    "id": "bf357417-7c81-4a3b-aa37-c168cfcd16cb",
    "metadata": {},
    "source": [
     "In Tamaas, we must define a **model** and a **solver**."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "id": "bcc86ebc-4b78-459f-9667-151990030d63",
    "metadata": {},
    "outputs": [],
    "source": [
     "model = tm.ModelFactory.createModel(tm.model_type.basic_2d, [L, L], [n, n])\n",
     "model.E = E\n",
     "model.nu = nu"
    ]
   },
   {
    "cell_type": "markdown",
    "id": "386442ba-12ca-4fab-8622-f6b8ec37a94a",
    "metadata": {},
    "source": [
     "We create variables to access the displacements and the tractions at the surface:"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "id": "c4af11d9-007c-4d99-91bd-e143319c01b2",
    "metadata": {},
    "outputs": [],
    "source": [
     "u = model.displacement\n",
     "p = model.traction"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "id": "7e8e084d-f442-42bf-8564-e21cefb99f57",
    "metadata": {},
    "outputs": [],
    "source": [
     "solver = tm.PolonskyKeerRey(model, h, 1e-12)"
    ]
   },
   {
    "cell_type": "markdown",
    "id": "b057fc64-a939-4eb0-915e-a72fa7c263cf",
    "metadata": {},
    "source": [
     "The creation of the solver takes three arguments: a model, the shape of the indenter, and a numerical accuracy for the solve."
    ]
   },
   {
    "cell_type": "markdown",
    "id": "2cf8c5f5-84a1-4072-9379-da3e1360b288",
    "metadata": {},
    "source": [
     "### Solve"
    ]
   },
   {
    "cell_type": "markdown",
    "id": "8d9e028f-37cc-4f34-8641-68dbd797cf9e",
    "metadata": {},
    "source": [
     "The solver needs an **average pressure**, obtained from the normal load and the total area of the system."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "id": "1ac2ba4b-4ca3-449a-80d4-c56e56a636b6",
    "metadata": {},
    "outputs": [],
    "source": [
     "solver.solve(F / (L**2))"
    ]
   },
   {
    "cell_type": "markdown",
    "id": "5265926a-4b28-4f03-9465-8d16de1b56f3",
    "metadata": {},
    "source": [
     "We can plot the displacements to see that the elastic surface matches the shape of the sphere."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "id": "ae0e5655-16ca-4e5f-b99a-2b34392a942e",
    "metadata": {},
    "outputs": [],
    "source": [
     "plt.figure()\n",
     "plt.plot(x[:, n//2], h[:, n//2], \"k\", label=\"indenter\")\n",
     "plt.plot(x[:, n//2], u[:, n//2], label=\"elastic surface\")\n",
     "# plt.axis(\"equal\")\n",
     "plt.xlabel(\"$x$\")\n",
     "plt.ylabel(\"$u$\")\n",
     "plt.legend()\n",
     "plt.title(\"Deformation\")\n",
     "plt.show()"
    ]
   },
   {
    "cell_type": "markdown",
    "id": "6cbfcea1-5daf-4a3c-b5a4-520f0249e48b",
    "metadata": {},
    "source": [
     "### Comparison with theory"
    ]
   },
   {
    "cell_type": "markdown",
    "id": "02d2211e-3e48-4328-bd27-020cdaa787d8",
    "metadata": {},
    "source": [
     "The surface tractions obtained by solving the contact with Tamaas can be compared to the theoretical pressure distribution. Complete the function below to compute this theoretical distribution."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "id": "555ad42b-b5d1-4956-87d6-4f881696d994",
    "metadata": {},
    "outputs": [],
    "source": [
     "def compute_p_theory(x, y, E, nu, F):\n",
     "    return # DEFINE THEORETICAL PRESSURE DISTRIBUTION ##########################\n",
     "\n",
     "p_theory = compute_p_theory(x, y, E, nu, F)"
    ]
   },
   {
    "cell_type": "markdown",
    "id": "c7c763f8-6f88-4f10-b601-c1fa7cb8cdfe",
    "metadata": {},
    "source": [
     "Both pressure distributions can be compared visually."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "id": "80cc852d-0c8c-4ba7-a612-9c2278d3ef24",
    "metadata": {},
    "outputs": [],
    "source": [
     "plt.figure()\n",
     "plt.plot(x[:, n//2], p_theory[:, n//2], \"k--\", label=\"theory\")\n",
     "plt.plot(x[:, n//2], p[:, n//2], label=\"simulation\")\n",
     "plt.xlabel(\"$x$\")\n",
     "plt.ylabel(\"$u$\")\n",
     "plt.legend()\n",
     "plt.title(\"Pressure\")\n",
     "plt.show()"
    ]
   },
   {
    "cell_type": "markdown",
    "id": "3d1fc51a-d953-466c-888b-1e7a6d7b0c46",
    "metadata": {},
    "source": [
     "...and numerically."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "id": "81b219d1-2809-4863-9f14-e37c540a117c",
    "metadata": {},
    "outputs": [],
    "source": [
     "print(\n",
     "\"\"\"Pressure at center:\n",
     "    theory:     {}\n",
     "    simulation: {}\"\"\"\n",
     ".format(p_theory[n//2, n//2], p[n//2, n//2]))"
    ]
   },
   {
    "cell_type": "markdown",
    "id": "f31fb95c-b93d-4f79-9a19-b87223c63441",
    "metadata": {},
    "source": [
     "## Convergence analysis"
    ]
   },
   {
    "cell_type": "markdown",
    "id": "4c54676e-902a-4530-b7e9-a3443f7a6235",
    "metadata": {},
    "source": [
     "The accuracy of the simulation depends on multiple factors. Write a function to numerically quantify the error between a simulated and a theoretical pressure distribution."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "id": "9867b042-10d6-4762-9ce6-93064f7fd80b",
    "metadata": {},
    "outputs": [],
    "source": [
     "def compute_error(p, p_theory):\n",
     "   return # DEFINE AN ERROR MEASUREMENT ########################################"
    ]
   },
   {
    "cell_type": "markdown",
    "id": "6d6cf67b-78b4-41e7-87d9-6cbbc59a205c",
    "metadata": {},
    "source": [
     "We define below a function to solve for the contact given some parameters: the discretization $n$, the normal load $F$ and the radius $R$."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "id": "bcc2a1c6-51f6-47d9-aed9-d0a1d4b485c8",
    "metadata": {},
    "outputs": [],
    "source": [
     "def solve(n, F, R):\n",
     "    # Goemetry\n",
     "    x = np.linspace(-L/2, L/2, n)\n",
     "    y = x.copy()\n",
     "    x, y = np.meshgrid(x, y, indexing=\"ij\")\n",
     "    h = create_indenter(x, y, R)\n",
     "    \n",
     "    # Model\n",
     "    model = tm.ModelFactory.createModel(tm.model_type.basic_2d, [L, L], [n, n])\n",
     "    model.E = E\n",
     "    model.nu = nu\n",
     "    p = model.traction\n",
     "    \n",
     "    # Solver\n",
     "    solver = tm.PolonskyKeerRey(model, h, 1e-12)\n",
     "    \n",
     "    # Solve\n",
     "    solver.solve(F / (L**2))\n",
     "    \n",
     "    return p, x, y"
    ]
   },
   {
    "cell_type": "markdown",
    "id": "6f77fafd-fd33-4126-8057-7f771c644e53",
    "metadata": {},
    "source": [
     "We can assess the effect of the discretization $n$ on the accuracy with a loop. Which values of $n$ give the best accuracy?"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "id": "c9c90509-e1c2-4861-aaf8-df47c0700b17",
    "metadata": {},
    "outputs": [],
    "source": [
     "n_list = [] # CREATE A LIST OF VALUES FOR n ####################################\n",
     "error_list = []\n",
     "\n",
     "F = 1e-4\n",
     "R = 10.0\n",
     "\n",
     "for n in n_list:\n",
     "    p, x, y = solve(n, F, R)\n",
     "    p_theory = compute_p_theory(x, y, E, nu, F)\n",
     "    error = compute_error(p, p_theory)\n",
     "    error_list.append(error)\n",
     "    print(\"n = {}, L/n = {:.6f}, a = {:.6f}, error = {:.2f}%\".format(n, L/n, (3*F*R/(4*E/(1-nu**2)))**(1/3), error * 100))\n",
     "\n",
     "plt.figure()\n",
     "plt.loglog(n_list, np.array(error_list) * 100, \".-\")\n",
     "plt.grid()\n",
     "plt.xlabel(\"$n$\")\n",
     "plt.ylabel(\"error [%]\")\n",
     "plt.title(\"Convergence\")\n",
     "plt.show()"
    ]
   },
   {
    "cell_type": "markdown",
    "id": "af2deaaa-3907-43bf-9c66-8c01e97f97f3",
    "metadata": {},
    "source": [
     "Now, perform the same kind of loop, but with other varying parameters. For each simulation, see how the contact radius $a$ compares with the system size $L$ and the grid size $L/n$, and find the conditions under which the accuracy is better."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "id": "b469a7bb-eadc-4972-9b4a-040b9a7359ef",
    "metadata": {},
    "outputs": [],
    "source": [
     "# WRITE OTHER TESTS ############################################################"
    ]
   }
  ],
  "metadata": {
   "kernelspec": {
    "display_name": "Python",
    "language": "python",
    "name": "python3"
   },
   "language_info": {
    "codemirror_mode": {
     "name": "ipython",
     "version": 3
    },
    "file_extension": ".py",
    "mimetype": "text/x-python",
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
    "version": "3.8.10"
   }
  },
  "nbformat": 4,
  "nbformat_minor": 5
 }