diff --git a/book/_toc.yml b/book/_toc.yml
index 17f7119..4865117 100644
--- a/book/_toc.yml
+++ b/book/_toc.yml
@@ -30,4 +30,5 @@ parts:
     - file: content/applications/task08
     - file: content/applications/task09
     - file: content/applications/task06
-    - file: content/applications/task10
\ No newline at end of file
+    - file: content/applications/task10
+    - file: content/advection
\ No newline at end of file
diff --git a/book/content/advection.ipynb b/book/content/advection.ipynb
new file mode 100644
index 0000000..dcd6a2e
--- /dev/null
+++ b/book/content/advection.ipynb
@@ -0,0 +1,389 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "b782b819",
+   "metadata": {},
+   "source": [
+    "# Advection-diffusion problem (multi-material)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "e4509407",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from mpi4py import MPI\n",
+    "\n",
+    "import dolfinx\n",
+    "import numpy as np\n",
+    "import basix\n",
+    "import festim as F\n",
+    "\n",
+    "n =50\n",
+    "\n",
+    "def bottom_boundary(x):\n",
+    "    return np.isclose(x[1], 0.0)\n",
+    "\n",
+    "def top_boundary(x):\n",
+    "    return np.isclose(x[1], 1.0)\n",
+    "\n",
+    "def bottom_left_boundary(x):\n",
+    "    return np.logical_and(np.isclose(x[0], 0.0), x[1] <= 0.5 + 1e-14)\n",
+    "\n",
+    "def half(x):\n",
+    "    return x[1] <= 0.5 + 1e-14\n",
+    "\n",
+    "mesh = dolfinx.mesh.create_unit_square(\n",
+    "    MPI.COMM_WORLD, n, n, dolfinx.mesh.CellType.triangle\n",
+    ")\n",
+    "\n",
+    "# Split domain in half and set an interface tag of 5\n",
+    "gdim = mesh.geometry.dim\n",
+    "tdim = mesh.topology.dim\n",
+    "fdim = tdim - 1\n",
+    "top_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, top_boundary)\n",
+    "bottom_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, bottom_boundary)\n",
+    "bot_left_facets = dolfinx.mesh.locate_entities_boundary(\n",
+    "    mesh, fdim, bottom_left_boundary\n",
+    ")\n",
+    "num_facets_local = (\n",
+    "    mesh.topology.index_map(fdim).size_local\n",
+    "    + mesh.topology.index_map(fdim).num_ghosts\n",
+    ")\n",
+    "facets = np.arange(num_facets_local, dtype=np.int32)\n",
+    "values = np.full_like(facets, 0, dtype=np.int32)\n",
+    "values[top_facets] = 1\n",
+    "values[bottom_facets] = 2\n",
+    "values[bot_left_facets] = 3\n",
+    "\n",
+    "bottom_cells = dolfinx.mesh.locate_entities(mesh, tdim, half)\n",
+    "num_cells_local = (\n",
+    "    mesh.topology.index_map(tdim).size_local\n",
+    "    + mesh.topology.index_map(tdim).num_ghosts\n",
+    ")\n",
+    "cells = np.full(num_cells_local, 5, dtype=np.int32)\n",
+    "cells[bottom_cells] = 4\n",
+    "ct = dolfinx.mesh.meshtags(\n",
+    "    mesh, tdim, np.arange(num_cells_local, dtype=np.int32), cells\n",
+    ")\n",
+    "all_b_facets = dolfinx.mesh.compute_incident_entities(\n",
+    "    mesh.topology, ct.find(4), tdim, fdim\n",
+    ")\n",
+    "all_t_facets = dolfinx.mesh.compute_incident_entities(\n",
+    "    mesh.topology, ct.find(5), tdim, fdim\n",
+    ")\n",
+    "interface = np.intersect1d(all_b_facets, all_t_facets)\n",
+    "values[interface] = 6\n",
+    "\n",
+    "mt = dolfinx.mesh.meshtags(mesh, mesh.topology.dim - 1, facets, values)\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "a91742d1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def create_velocity_field():\n",
+    "    def velocity_func(x):\n",
+    "        values = np.zeros((2, x.shape[1]))  # Initialize with zeros\n",
+    "        scalar_value = -1000 * x[1] * (x[1] - 0.5)  # Compute the scalar function\n",
+    "        values[0] = scalar_value  # Assign to first component\n",
+    "        values[1] = 0  # Second component remains zero\n",
+    "        return values\n",
+    "\n",
+    "    submesh, submesh_to_mesh, v_map = dolfinx.mesh.create_submesh(\n",
+    "        mesh, ct.dim, ct.find(4)\n",
+    "    )[0:3]\n",
+    "    v_cg = basix.ufl.element(\n",
+    "        \"Lagrange\", submesh.topology.cell_name(), 2, shape=(submesh.geometry.dim,)\n",
+    "    )\n",
+    "    V_velocity = dolfinx.fem.functionspace(submesh, v_cg)\n",
+    "    u = dolfinx.fem.Function(V_velocity)\n",
+    "    u.interpolate(velocity_func)\n",
+    "    return u\n",
+    "\n",
+    "\n",
+    "u = create_velocity_field()\n",
+    "\n",
+    "writer = dolfinx.io.VTXWriter(MPI.COMM_WORLD, \"velocity.bp\", u, engine=\"BP5\")\n",
+    "writer.write(t=0)\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "id": "c2a88a52",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from dolfinx import plot\n",
+    "import pyvista\n",
+    "\n",
+    "pyvista.start_xvfb()\n",
+    "pyvista.set_jupyter_backend(\"html\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "id": "42c527ce",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "ec6a05530c26467c8338481a55dcb5f5",
+       "version_major": 2,
+       "version_minor": 0
+      },
+      "text/plain": [
+       "EmbeddableWidget(value='<iframe srcdoc=\"<!DOCTYPE html>\\n<html>\\n  <head>\\n    <meta http-equiv=&quot;Content-…"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "topology, cell_types, geometry = plot.vtk_mesh(u.function_space)\n",
+    "values = np.zeros((geometry.shape[0], 3), dtype=np.float64)\n",
+    "values[:, :len(u)] = u.x.array.real.reshape((geometry.shape[0], len(u)))\n",
+    "\n",
+    "# Create a point cloud of glyphs\n",
+    "function_grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)\n",
+    "function_grid[\"v\"] = values\n",
+    "glyphs = function_grid.glyph(orient=\"v\", factor=0.0005)\n",
+    "\n",
+    "# Create a pyvista-grid for the mesh\n",
+    "grid = pyvista.UnstructuredGrid(*plot.vtk_mesh(mesh, mesh.topology.dim))\n",
+    "\n",
+    "# Create plotter\n",
+    "plotter = pyvista.Plotter()\n",
+    "plotter.add_mesh(grid, style=\"wireframe\", color=\"k\", opacity=0.1)\n",
+    "plotter.add_mesh(glyphs)\n",
+    "plotter.view_xy()\n",
+    "if not pyvista.OFF_SCREEN:\n",
+    "    plotter.show()\n",
+    "else:\n",
+    "    fig_as_array = plotter.screenshot(\"glyphs.png\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "id": "1fe8218b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "my_model = F.HydrogenTransportProblemDiscontinuous()\n",
+    "my_model.mesh = F.Mesh(mesh)\n",
+    "my_model.volume_meshtags = ct\n",
+    "my_model.facet_meshtags = mt\n",
+    "\n",
+    "material_bottom = F.Material(D_0=0.5, E_D=0 * 0.1)\n",
+    "material_top = F.Material(D_0=0.5, E_D=0 * 0.1)\n",
+    "\n",
+    "material_bottom.K_S_0 = 2.0\n",
+    "material_bottom.E_K_S = 0\n",
+    "material_top.K_S_0 = 6.0\n",
+    "material_top.E_K_S = 0\n",
+    "\n",
+    "top_domain = F.VolumeSubdomain(5, material=material_top)\n",
+    "bottom_domain = F.VolumeSubdomain(4, material=material_bottom)\n",
+    "\n",
+    "top_surface = F.SurfaceSubdomain(id=1)\n",
+    "bottom_surface = F.SurfaceSubdomain(id=2)\n",
+    "inlet_surf = F.SurfaceSubdomain(id=3)\n",
+    "my_model.subdomains = [\n",
+    "    bottom_domain,\n",
+    "    top_domain,\n",
+    "    top_surface,\n",
+    "    bottom_surface,\n",
+    "    inlet_surf,\n",
+    "]\n",
+    "\n",
+    "my_model.interfaces = [F.Interface(6, (bottom_domain, top_domain), penalty_term=5000)]\n",
+    "\n",
+    "my_model.surface_to_volume = {\n",
+    "    top_surface: top_domain,\n",
+    "    bottom_surface: bottom_domain,\n",
+    "    inlet_surf: bottom_domain,\n",
+    "}\n",
+    "\n",
+    "H = F.Species(\"H\", mobile=True)\n",
+    "\n",
+    "my_model.species = [H]\n",
+    "\n",
+    "for species in my_model.species:\n",
+    "    species.subdomains = [bottom_domain, top_domain]\n",
+    "\n",
+    "\n",
+    "my_model.boundary_conditions = [\n",
+    "    F.FixedConcentrationBC(top_surface, value=0.0, species=H),\n",
+    "    F.FixedConcentrationBC(inlet_surf, value=1.0, species=H),\n",
+    "    F.FixedConcentrationBC(bottom_surface, value=0.0, species=H),\n",
+    "]\n",
+    "\n",
+    "my_model.advection_terms = [\n",
+    "    F.AdvectionTerm(velocity=u, species=[H], subdomain=bottom_domain)\n",
+    "]\n",
+    "\n",
+    "\n",
+    "my_model.temperature = 500.0\n",
+    "\n",
+    "my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, transient=False)\n",
+    "\n",
+    "my_model.exports = [\n",
+    "    F.VTXSpeciesExport(f\"u_{subdomain.id}.bp\", field=H, subdomain=subdomain)\n",
+    "    for subdomain in my_model.volume_subdomains\n",
+    "]\n",
+    "\n",
+    "my_model.initialise()\n",
+    "my_model.run()\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "id": "55fa771a",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "96402a37d52f413ab81b52d14c693437",
+       "version_major": 2,
+       "version_minor": 0
+      },
+      "text/plain": [
+       "EmbeddableWidget(value='<iframe srcdoc=\"<!DOCTYPE html>\\n<html>\\n  <head>\\n    <meta http-equiv=&quot;Content-…"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "def plot_in_one_subdomain(subdomain, species, label=\"c\"):\n",
+    "    u = species.subdomain_to_post_processing_solution[subdomain]\n",
+    "    topology, cell_types, geometry = plot.vtk_mesh(u.function_space)\n",
+    "    u_grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)\n",
+    "    u_grid.point_data[label] = u.x.array.real\n",
+    "    u_grid.set_active_scalars(label)\n",
+    "    return u_grid\n",
+    "\n",
+    "u_plotter = pyvista.Plotter()\n",
+    "\n",
+    "for subdomain in my_model.volume_subdomains:\n",
+    "    u_grid = plot_in_one_subdomain(subdomain, H)\n",
+    "    u_plotter.add_mesh(u_grid, cmap=\"viridis\", show_edges=False)\n",
+    "    u_plotter.add_mesh(u_grid, style=\"wireframe\", color=\"white\", opacity=0.1)\n",
+    "    u_plotter.view_xy()\n",
+    "\n",
+    "if not pyvista.OFF_SCREEN:\n",
+    "    u_plotter.show()\n",
+    "else:\n",
+    "    figure = u_plotter.screenshot(\"concentration.png\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "id": "23b30e74",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "2aa9b335a45949898de63103551d5efa",
+       "version_major": 2,
+       "version_minor": 0
+      },
+      "text/plain": [
+       "EmbeddableWidget(value='<iframe srcdoc=\"<!DOCTYPE html>\\n<html>\\n  <head>\\n    <meta http-equiv=&quot;Content-…"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "u_plotter = pyvista.Plotter()\n",
+    "\n",
+    "u_grid = plot_in_one_subdomain(my_model.volume_subdomains[0], H)\n",
+    "u_grid.set_active_scalars(\"c\")\n",
+    "u_plotter.add_mesh(u_grid, cmap=\"viridis\", show_edges=False)\n",
+    "u_plotter.add_mesh(u_grid, style=\"wireframe\", color=\"white\", opacity=0.1)\n",
+    "u_plotter.view_xy()\n",
+    "\n",
+    "if not pyvista.OFF_SCREEN:\n",
+    "    u_plotter.show()\n",
+    "else:\n",
+    "    figure = u_plotter.screenshot(\"concentration.png\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "id": "725fe826",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "7507c1e5d95c47b585108f7c0e3e878b",
+       "version_major": 2,
+       "version_minor": 0
+      },
+      "text/plain": [
+       "EmbeddableWidget(value='<iframe srcdoc=\"<!DOCTYPE html>\\n<html>\\n  <head>\\n    <meta http-equiv=&quot;Content-…"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "u_plotter = pyvista.Plotter()\n",
+    "\n",
+    "u_grid = plot_in_one_subdomain(my_model.volume_subdomains[1], H)\n",
+    "u_grid.set_active_scalars(\"c\")\n",
+    "u_plotter.add_mesh(u_grid, cmap=\"viridis\", show_edges=False)\n",
+    "u_plotter.add_mesh(u_grid, style=\"wireframe\", color=\"white\", opacity=0.1)\n",
+    "u_plotter.view_xy()\n",
+    "\n",
+    "if not pyvista.OFF_SCREEN:\n",
+    "    u_plotter.show()\n",
+    "else:\n",
+    "    figure = u_plotter.screenshot(\"concentration.png\")"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "festim-workshop",
+   "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.13.5"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}