diff --git a/.github/workflows/book_stable.yml b/.github/workflows/book_stable.yml
index aaac9969..f398eb76 100644
--- a/.github/workflows/book_stable.yml
+++ b/.github/workflows/book_stable.yml
@@ -26,7 +26,7 @@ jobs:
       PYVISTA_JUPYTER_BACKEND: "html"
 
     steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v5
 
       - name: Install common packages
         uses: ./.github/actions/install-dependencies
diff --git a/.github/workflows/deploy.yml b/.github/workflows/deploy.yml
index 3e467708..8fd83a60 100644
--- a/.github/workflows/deploy.yml
+++ b/.github/workflows/deploy.yml
@@ -35,7 +35,7 @@ jobs:
 
     steps:
       - name: Checkout
-        uses: actions/checkout@v4
+        uses: actions/checkout@v5
 
       - name: Setup Pages
         uses: actions/configure-pages@v5
diff --git a/.github/workflows/publish_docker.yml b/.github/workflows/publish_docker.yml
index bd09b935..c2b9133f 100644
--- a/.github/workflows/publish_docker.yml
+++ b/.github/workflows/publish_docker.yml
@@ -28,7 +28,7 @@ jobs:
 
     steps:
       - name: Checkout repository
-        uses: actions/checkout@v4
+        uses: actions/checkout@v5
 
       - name: Log in to the Container registry
         uses: docker/login-action@v3
diff --git a/.github/workflows/test_nightly.yml b/.github/workflows/test_nightly.yml
index ec345863..d1267d99 100644
--- a/.github/workflows/test_nightly.yml
+++ b/.github/workflows/test_nightly.yml
@@ -27,15 +27,15 @@ jobs:
       LIBGL_ALWAYS_SOFTWARE: 1
 
     steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v5
 
       - name: Special handling of some installation
         uses: ./.github/actions/install-dependencies
 
       - name: Install requirements
         run: |
-          python3 -m pip install --break-system-packages -U pip setuptools pkgconfig
-          python3 -m pip install --no-build-isolation --break-system-packages --no-cache-dir --no-binary=h5py . --upgrade
+          python3 -m pip install --break-system-packages -U pip setuptools pkgconfig poetry-core
+          python3 -m pip install --no-build-isolation --break-system-packages --no-cache-dir --no-binary=h5py .[netgen] --upgrade
 
       - name: Test building the book
         run: PYVISTA_OFF_SCREEN=false jupyter-book build  -W .
diff --git a/.github/workflows/test_stable.yml b/.github/workflows/test_stable.yml
index 64e4adc3..eacbe4c9 100644
--- a/.github/workflows/test_stable.yml
+++ b/.github/workflows/test_stable.yml
@@ -21,7 +21,7 @@ jobs:
 
     # Steps represent a sequence of tasks that will be executed as part of the job
     steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v5
         with:
           ref: release
 
@@ -29,7 +29,8 @@ jobs:
 
       - name: Install additional deps
         run: |
-          python3 -m pip install --no-binary=h5py .
+          python3 -m pip install -U pip setuptools pkgconfig poetry-core
+          python3 -m pip install --no-binary=h5py .[netgen]
 
       - name: Test complex notebooks in parallel
         working-directory: chapter1
diff --git a/_toc.yml b/_toc.yml
index 6b7b6301..3da48c6e 100644
--- a/_toc.yml
+++ b/_toc.yml
@@ -16,6 +16,7 @@ parts:
         sections:
           - file: chapter1/membrane_code
           - file: chapter1/membrane_paraview
+
   - caption: A Gallery of finite element solvers
     chapters:
       - file: chapter2/intro
@@ -39,6 +40,7 @@ parts:
       - file: chapter2/helmholtz
         sections:
           - file: chapter2/helmholtz_code
+      - file: chapter2/amr
 
   - caption: Subdomains and boundary conditions
     chapters:
diff --git a/chapter2/amr.ipynb b/chapter2/amr.ipynb
new file mode 100644
index 00000000..28bcf9d8
--- /dev/null
+++ b/chapter2/amr.ipynb
@@ -0,0 +1,510 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "4535be91",
+   "metadata": {},
+   "source": [
+    "# Adaptive mesh refinement with NetGen and DOLFINx\n",
+    "\n",
+    "Author: Jørgen S. Dokken\n",
+    "\n",
+    "```{admonition} NetGen and linux/arm64\n",
+    "NetGen is not available on linux/arm64, so this tutorial will not work on that platform.\n",
+    "If you are on such a platform, one can use the AMD docker images, with a performance penalty due to emulation.\n",
+    "```"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "b0282fe2",
+   "metadata": {},
+   "source": [
+    "In this tutorial, we will consider an adaptive mesh refinement method, applied to\n",
+    "the Laplace eigenvalue problem.\n",
+    "This demo is an adaptation of [Firedrake - Adaptive Mesh Refinement](https://www.firedrakeproject.org/firedrake/demos/netgen_mesh.py.html).\n",
+    "In this tutorial we will use the mesh generator [NetGen](https://ngsolve.org/) from NGSolve.\n",
+    "First, we import the packages needed for this demo:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "32fe6a58",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from mpi4py import MPI\n",
+    "from petsc4py import PETSc\n",
+    "from slepc4py import SLEPc\n",
+    "from packaging.version import Version\n",
+    "import dolfinx.fem.petsc\n",
+    "import numpy as np\n",
+    "import ufl\n",
+    "import pyvista\n",
+    "import ngsPETSc.utils.fenicsx as ngfx\n",
+    "\n",
+    "from netgen.geom2d import SplineGeometry"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "4782ba3b",
+   "metadata": {},
+   "source": [
+    "## Generating a higher-order mesh with NetGen\n",
+    "Next, we generate a PacMan-like geometry using NetGen."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d6046434",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "geo = SplineGeometry()\n",
+    "pnts = [(0, 0), (1, 0), (1, 1), (0, 1), (-1, 1), (-1, 0), (-1, -1), (0, -1)]\n",
+    "p1, p2, p3, p4, p5, p6, p7, p8 = [geo.AppendPoint(*pnt) for pnt in pnts]\n",
+    "curves = [\n",
+    "    [[\"line\", p1, p2], \"line\"],\n",
+    "    [[\"spline3\", p2, p3, p4], \"curve\"],\n",
+    "    [[\"spline3\", p4, p5, p6], \"curve\"],\n",
+    "    [[\"spline3\", p6, p7, p8], \"curve\"],\n",
+    "    [[\"line\", p8, p1], \"line\"],\n",
+    "]\n",
+    "for c, bc in curves:\n",
+    "    geo.Append(c, bc=bc)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "e041abef",
+   "metadata": {},
+   "source": [
+    "## Loading a mesh into DOLFINx\n",
+    "The ngsPETSc package provides a communication layer between NetGen and DOLFINx.\n",
+    "We initialize this layer by passing in a NetGen-model, as well as an MPI communicator,\n",
+    "which will be used to distribute the mesh."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "fc64bf64",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "geoModel = ngfx.GeometricModel(geo, MPI.COMM_WORLD)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "67cded8b",
+   "metadata": {},
+   "source": [
+    "Next, we generate the mesh with the function :py:func:`ngsPETSc.utils.fenicsx.GeometricModel.model_to_mesh`.\n",
+    "Which takes in the target geometric dimension of the mesh (2 for triangular meshes, 3 for tetrahedral), the\n",
+    "maximum mesh size (`hmax`) and a few optional parameters."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "dbae564a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mesh, (ct, ft), region_map = geoModel.model_to_mesh(gdim=2, hmax=0.5)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "f90ca932",
+   "metadata": {},
+   "source": [
+    "We use pyvista to visualize the mesh."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "34407248",
+   "metadata": {
+    "tags": [
+     "hide-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "pyvista.start_xvfb(1.0)\n",
+    "\n",
+    "grid = pyvista.UnstructuredGrid(*dolfinx.plot.vtk_mesh(mesh))\n",
+    "grid.cell_data[\"ct\"] = ct.values\n",
+    "\n",
+    "plotter = pyvista.Plotter()\n",
+    "plotter.add_mesh(\n",
+    "    grid, show_edges=True, scalars=\"ct\", cmap=\"blues\", show_scalar_bar=False\n",
+    ")\n",
+    "plotter.view_xy()\n",
+    "if not pyvista.OFF_SCREEN:\n",
+    "    plotter.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "6a18b03f",
+   "metadata": {},
+   "source": [
+    "We have read in any cell and facet markers that have been defined in the NetGen model,\n",
+    "as well as a map from their names to their integer ids in `ct`, `ft` and `region_map` respectively.\n",
+    "We can curve the grids with the command `curveField`.\n",
+    "In this example, we use third order Lagrange elements to represent the geometry."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "adc50da0",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "order = 3\n",
+    "curved_mesh = geoModel.curveField(order)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "0a3449b0",
+   "metadata": {},
+   "source": [
+    "Again, we visualize the curved mesh with pyvista."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1eb97f0f",
+   "metadata": {
+    "tags": [
+     "hide-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "curved_grid = pyvista.UnstructuredGrid(*dolfinx.plot.vtk_mesh(curved_mesh))\n",
+    "curved_grid.cell_data[\"ct\"] = ct.values\n",
+    "plotter.add_mesh(\n",
+    "    curved_grid, show_edges=False, scalars=\"ct\", cmap=\"blues\", show_scalar_bar=False\n",
+    ")\n",
+    "plotter.add_mesh(grid, style=\"wireframe\", color=\"black\")\n",
+    "plotter.view_xy()\n",
+    "if not pyvista.OFF_SCREEN:\n",
+    "    plotter.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "7613cc7f",
+   "metadata": {},
+   "source": [
+    "## Solving the eigenvalue problem\n",
+    "In this section we will solve the eigenvalue problem:\n",
+    "\n",
+    "Find $u_h\\in H_0^1(\\Omega)$ and $\\lambda\\in\\mathbb{R}$ such that\n",
+    "\n",
+    "$$\n",
+    "\\begin{align}\n",
+    "\\int_\\Omega \\nabla u \\cdot \\nabla v~\\mathrm{d} x &= \\lambda \\int_\\Omega u v~\\mathrm{d} x \\qquad\n",
+    "\\forall v \\in H_0^1(\\Omega).\n",
+    "\\end{align}\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "717f48a2",
+   "metadata": {
+    "lines_to_next_cell": 2
+   },
+   "source": [
+    "Next, we define a convenience function to solve the eigenvalue problem using [SLEPc](https://slepc.upv.es/)\n",
+    "given a discretized domain, its facet markers and the region map."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d59bc925",
+   "metadata": {
+    "lines_to_next_cell": 2
+   },
+   "outputs": [],
+   "source": [
+    "def solve(\n",
+    "    mesh: dolfinx.mesh.Mesh,\n",
+    "    facet_tags: dolfinx.mesh.MeshTags,\n",
+    "    region_map: dict[tuple[int, str], tuple[int, ...]],\n",
+    ") -> tuple[float, dolfinx.fem.Function, dolfinx.fem.Function]:\n",
+    "    # We define the lhs and rhs bilinear forms\n",
+    "    V = dolfinx.fem.functionspace(mesh, (\"Lagrange\", 3))\n",
+    "    u = ufl.TrialFunction(V)\n",
+    "    v = ufl.TestFunction(V)\n",
+    "    a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx\n",
+    "    m = ufl.inner(u, v) * ufl.dx\n",
+    "\n",
+    "    # We identify the boundary facets and their corresponding dofs\n",
+    "    straight_facets = facet_tags.indices[\n",
+    "        np.isin(facet_tags.values, region_map[(1, \"line\")])\n",
+    "    ]\n",
+    "    curved_facets = facet_tags.indices[\n",
+    "        np.isin(facet_tags.values, region_map[(1, \"curve\")])\n",
+    "    ]\n",
+    "    boundary_facets = np.concatenate([straight_facets, curved_facets])\n",
+    "    mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)\n",
+    "    boundary_dofs = dolfinx.fem.locate_dofs_topological(\n",
+    "        V, mesh.topology.dim - 1, boundary_facets\n",
+    "    )\n",
+    "\n",
+    "    # We create a zero boundary condition for these dofs to be in the suitable space, and\n",
+    "    # set up the discrete matrices `A` and `M`\n",
+    "    bc = dolfinx.fem.dirichletbc(0.0, boundary_dofs, V)\n",
+    "    A = dolfinx.fem.petsc.assemble_matrix(dolfinx.fem.form(a), bcs=[bc])\n",
+    "    A.assemble()\n",
+    "    if Version(dolfinx.__version__) < Version(\"0.10.0\"):\n",
+    "        diag_kwargs = {\"diagonal\": 0.0}\n",
+    "    else:\n",
+    "        diag_kwargs = {\"diag\": 0.0}\n",
+    "\n",
+    "    M = dolfinx.fem.petsc.assemble_matrix(dolfinx.fem.form(m), bcs=[bc], **diag_kwargs)\n",
+    "    M.assemble()\n",
+    "\n",
+    "    # Next, we define the SLEPc Eigenvalue Problem Solver (EPS), and set up to use a shift\n",
+    "    # and invert (SINVERT) spectral transformation where the preconditioner factorisation\n",
+    "    # is computed using [MUMPS](https://mumps-solver.org/index.php).\n",
+    "\n",
+    "    E = SLEPc.EPS().create(mesh.comm)\n",
+    "    E.setType(SLEPc.EPS.Type.ARNOLDI)\n",
+    "    E.setProblemType(SLEPc.EPS.ProblemType.GHEP)\n",
+    "    E.setDimensions(1, SLEPc.DECIDE)\n",
+    "    E.setOperators(A, M)\n",
+    "    ST = E.getST()\n",
+    "    ST.setType(SLEPc.ST.Type.SINVERT)\n",
+    "    PC = ST.getKSP().getPC()\n",
+    "    PC.setType(\"lu\")\n",
+    "    PC.setFactorSolverType(\"mumps\")\n",
+    "    E.setST(ST)\n",
+    "    E.solve()\n",
+    "    assert E.getConvergedReason() >= 0, \"Eigenvalue solver did not converge\"\n",
+    "\n",
+    "    # We get the real and imaginary parts of the first eigenvector along with the eigenvalue.\n",
+    "    uh_r = dolfinx.fem.Function(V)\n",
+    "    uh_i = dolfinx.fem.Function(V)\n",
+    "    lam = E.getEigenpair(0, uh_r.x.petsc_vec, uh_i.x.petsc_vec)\n",
+    "    E.destroy()\n",
+    "    uh_r.x.scatter_forward()\n",
+    "    uh_i.x.scatter_forward()\n",
+    "    return (lam, uh_r, uh_i)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "1c251b79",
+   "metadata": {
+    "lines_to_next_cell": 2
+   },
+   "source": [
+    "## Error-indicator\n",
+    "In this example, we will use an error-indicator $\\eta$ to decide what cells should be refined.\n",
+    "Specifically, the estimator is:\n",
+    "\n",
+    "\\begin{align*}\n",
+    " \\eta = \\sum_{K\\in \\mathcal{T}_h(\\Omega)}\\left(h^2\\int_K \\vert \\lambda u_h + \\Delta u_h\\vert^2~\\mathrm{d}x\\right)\n",
+    "+ \\sum_{E\\in\\mathcal{F}_i}\\frac{h}{2} \\vert [\\nabla \\cdot \\mathbf{n}_E ]\\vert^2~\\mathrm{d}s\n",
+    "\\end{align*}\n",
+    "\n",
+    "where $\\mathcal{T}_h$ is the collection of cells in the mesh, $\\mathcal{F}_i$ the collection"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "95f75bd8",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def mark_cells(uh_r: dolfinx.fem.Function, lam: float):\n",
+    "    mesh = uh_r.function_space.mesh\n",
+    "    W = dolfinx.fem.functionspace(mesh, (\"DG\", 0))\n",
+    "    w = ufl.TestFunction(W)\n",
+    "    eta = dolfinx.fem.Function(W)\n",
+    "    f = dolfinx.fem.Constant(mesh, 1.0)\n",
+    "    h = dolfinx.fem.Function(W)\n",
+    "    h.x.array[:] = mesh.h(mesh.topology.dim, np.arange(len(h.x.array), dtype=np.int32))\n",
+    "    n = ufl.FacetNormal(mesh)\n",
+    "\n",
+    "    G = (  # compute cellwise error estimator\n",
+    "        ufl.inner(h**2 * (f + ufl.div(ufl.grad(uh_r))) ** 2, w) * ufl.dx\n",
+    "        + ufl.inner(h(\"+\") / 2 * ufl.jump(ufl.grad(uh_r), n) ** 2, w(\"+\")) * ufl.dS\n",
+    "        + ufl.inner(h(\"-\") / 2 * ufl.jump(ufl.grad(uh_r), n) ** 2, w(\"-\")) * ufl.dS\n",
+    "    )\n",
+    "    dolfinx.fem.petsc.assemble_vector(eta.x.petsc_vec, dolfinx.fem.form(G))\n",
+    "    sqrt_eta = dolfinx.fem.Function(W)\n",
+    "    sqrt_eta.x.array[:] = np.sqrt(eta.x.array[:])\n",
+    "\n",
+    "    sqrt_eta_max = sqrt_eta.x.petsc_vec.max()[1]\n",
+    "\n",
+    "    theta = 0.5\n",
+    "    should_refine = ufl.conditional(ufl.gt(sqrt_eta, theta * sqrt_eta_max), 1, 0)\n",
+    "    markers = dolfinx.fem.Function(W)\n",
+    "    ip = W.element.interpolation_points\n",
+    "    if Version(dolfinx.__version__) < Version(\"0.10.0\"):\n",
+    "        ip = ip()\n",
+    "    markers.interpolate(dolfinx.fem.Expression(should_refine, ip))\n",
+    "    return np.flatnonzero(np.isclose(markers.x.array.astype(np.int32), 1))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "d9243490",
+   "metadata": {},
+   "source": [
+    "## Running the adaptive refinement algorithm\n",
+    "Next, we will run the adaptive mesh refinement algorithm."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "ab8f6701",
+   "metadata": {},
+   "source": [
+    "We will track the progress of the adaptive mesh refinement as a GIF."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "293b959a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plotter = pyvista.Plotter()\n",
+    "plotter.open_gif(\"amr.gif\", fps=1)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "5fe0ac8b",
+   "metadata": {
+    "lines_to_next_cell": 2
+   },
+   "source": [
+    "We make a convenience function to attach the relevant data to the plotter at a given\n",
+    "refinement step."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "969aaaa5",
+   "metadata": {
+    "tags": [
+     "hide-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "def write_frame(plotter: pyvista.Plotter, uh_r: dolfinx.fem.Function):\n",
+    "    # Scale uh_r to be consistent between refinement steps, as it can be multiplied by -1\n",
+    "    uh_r_min = curved_mesh.comm.allreduce(uh_r.x.array.min(), op=MPI.MIN)\n",
+    "    uh_r_max = curved_mesh.comm.allreduce(uh_r.x.array.max(), op=MPI.MAX)\n",
+    "    uh_sign = np.sign(uh_r_min)\n",
+    "    if np.isclose(uh_sign, 0):\n",
+    "        uh_sign = np.sign(uh_r_max)\n",
+    "    assert not np.isclose(uh_sign, 0), \"uh_r has zero values, cannot determine sign.\"\n",
+    "    uh_r.x.array[:] *= uh_sign\n",
+    "\n",
+    "    # Update plot with refined mesh\n",
+    "    grid = pyvista.UnstructuredGrid(*dolfinx.plot.vtk_mesh(mesh))\n",
+    "    curved_grid = pyvista.UnstructuredGrid(*dolfinx.plot.vtk_mesh(uh_r.function_space))\n",
+    "    curved_grid.point_data[\"u\"] = uh_r.x.array\n",
+    "    curved_grid = curved_grid.tessellate()\n",
+    "    curved_actor = plotter.add_mesh(\n",
+    "        curved_grid,\n",
+    "        show_edges=False,\n",
+    "    )\n",
+    "\n",
+    "    actor = plotter.add_mesh(grid, style=\"wireframe\", color=\"black\")\n",
+    "    plotter.view_xy()\n",
+    "    plotter.write_frame()\n",
+    "    plotter.remove_actor(actor)\n",
+    "    plotter.remove_actor(curved_actor)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "d798464c",
+   "metadata": {},
+   "source": [
+    "We set some parameters for checking convergence of the algorithm, and provide the exact eigenvalue\n",
+    "for comparison.\n",
+    "```{admonition} Using ngsPETSc for mesh refinement\n",
+    "In `ngsPETSc`, we provide the function `GeometricModel.refineMarkedElements` which we\n",
+    "pass the entities we would like to refine, and the topological dimensions of those entities.\n",
+    "The function returns a refined mesh, with corresponding cell and facet markers extracted from\n",
+    "the NetGen model.\n",
+    "```"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ebd3fef7",
+   "metadata": {
+    "tags": [
+     "scroll-output"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "max_iterations = 15\n",
+    "exact = 3.375610652693620492628**2\n",
+    "termination_criteria = 1e-5\n",
+    "for i in range(max_iterations):\n",
+    "    lam, uh_r, _ = solve(curved_mesh, ft, region_map)\n",
+    "\n",
+    "    relative_error = (lam - exact) / abs(exact)\n",
+    "    PETSc.Sys.Print(\n",
+    "        f\"Iteration {i + 1}/{max_iterations}, {lam=:.5e}, {exact=:.5e}, {relative_error=:.2e}\"\n",
+    "    )\n",
+    "\n",
+    "    cells_to_mark = mark_cells(uh_r, lam)\n",
+    "    mesh, (_, ft) = geoModel.refineMarkedElements(mesh.topology.dim, cells_to_mark)\n",
+    "    curved_mesh = geoModel.curveField(order)\n",
+    "    write_frame(plotter, uh_r)\n",
+    "\n",
+    "    if relative_error < termination_criteria:\n",
+    "        PETSc.Sys.Print(f\"Converged in {i + 1} iterations.\")\n",
+    "        break\n",
+    "plotter.close()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "5bac3843",
+   "metadata": {},
+   "source": [
+    "<img src=\"./amr.gif\" alt=\"gif\" class=\"bg-primary mb-1\" width=\"800px\">"
+   ]
+  }
+ ],
+ "metadata": {
+  "jupytext": {
+   "cell_metadata_filter": "tags,-all",
+   "formats": "ipynb,py:light",
+   "main_language": "python"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/chapter2/amr.py b/chapter2/amr.py
new file mode 100644
index 00000000..bd940f54
--- /dev/null
+++ b/chapter2/amr.py
@@ -0,0 +1,309 @@
+# ---
+# jupyter:
+#   jupytext:
+#     cell_metadata_filter: tags,-all
+#     formats: ipynb,py:light
+#     text_representation:
+#       extension: .py
+#       format_name: light
+#       format_version: '1.5'
+#       jupytext_version: 1.17.2
+# ---
+
+# # Adaptive mesh refinement with NetGen and DOLFINx
+#
+# Author: Jørgen S. Dokken
+#
+# ```{admonition} NetGen and linux/arm64
+# NetGen is not available on linux/arm64, so this tutorial will not work on that platform.
+# If you are on such a platform, one can use the AMD docker images, with a performance penalty due to emulation.
+# ```
+
+# In this tutorial, we will consider an adaptive mesh refinement method, applied to
+# the Laplace eigenvalue problem.
+# This demo is an adaptation of [Firedrake - Adaptive Mesh Refinement](https://www.firedrakeproject.org/firedrake/demos/netgen_mesh.py.html).
+# In this tutorial we will use the mesh generator [NetGen](https://ngsolve.org/) from NGSolve.
+# First, we import the packages needed for this demo:
+
+# +
+from mpi4py import MPI
+from petsc4py import PETSc
+from slepc4py import SLEPc
+from packaging.version import Version
+import dolfinx.fem.petsc
+import numpy as np
+import ufl
+import pyvista
+import ngsPETSc.utils.fenicsx as ngfx
+
+from netgen.geom2d import SplineGeometry
+# -
+
+# ## Generating a higher-order mesh with NetGen
+# Next, we generate a PacMan-like geometry using NetGen.
+
+geo = SplineGeometry()
+pnts = [(0, 0), (1, 0), (1, 1), (0, 1), (-1, 1), (-1, 0), (-1, -1), (0, -1)]
+p1, p2, p3, p4, p5, p6, p7, p8 = [geo.AppendPoint(*pnt) for pnt in pnts]
+curves = [
+    [["line", p1, p2], "line"],
+    [["spline3", p2, p3, p4], "curve"],
+    [["spline3", p4, p5, p6], "curve"],
+    [["spline3", p6, p7, p8], "curve"],
+    [["line", p8, p1], "line"],
+]
+for c, bc in curves:
+    geo.Append(c, bc=bc)
+
+# ## Loading a mesh into DOLFINx
+# The ngsPETSc package provides a communication layer between NetGen and DOLFINx.
+# We initialize this layer by passing in a NetGen-model, as well as an MPI communicator,
+# which will be used to distribute the mesh.
+
+geoModel = ngfx.GeometricModel(geo, MPI.COMM_WORLD)
+
+# Next, we generate the mesh with the function :py:func:`ngsPETSc.utils.fenicsx.GeometricModel.model_to_mesh`.
+# Which takes in the target geometric dimension of the mesh (2 for triangular meshes, 3 for tetrahedral), the
+# maximum mesh size (`hmax`) and a few optional parameters.
+
+mesh, (ct, ft), region_map = geoModel.model_to_mesh(gdim=2, hmax=0.5)
+
+# We use pyvista to visualize the mesh.
+
+# + tags=["hide-input"]
+pyvista.start_xvfb(1.0)
+
+grid = pyvista.UnstructuredGrid(*dolfinx.plot.vtk_mesh(mesh))
+grid.cell_data["ct"] = ct.values
+
+plotter = pyvista.Plotter()
+plotter.add_mesh(
+    grid, show_edges=True, scalars="ct", cmap="blues", show_scalar_bar=False
+)
+plotter.view_xy()
+if not pyvista.OFF_SCREEN:
+    plotter.show()
+# -
+
+# We have read in any cell and facet markers that have been defined in the NetGen model,
+# as well as a map from their names to their integer ids in `ct`, `ft` and `region_map` respectively.
+# We can curve the grids with the command `curveField`.
+# In this example, we use third order Lagrange elements to represent the geometry.
+
+order = 3
+curved_mesh = geoModel.curveField(order)
+
+# Again, we visualize the curved mesh with pyvista.
+
+# + tags=["hide-input"]
+curved_grid = pyvista.UnstructuredGrid(*dolfinx.plot.vtk_mesh(curved_mesh))
+curved_grid.cell_data["ct"] = ct.values
+plotter.add_mesh(
+    curved_grid, show_edges=False, scalars="ct", cmap="blues", show_scalar_bar=False
+)
+plotter.add_mesh(grid, style="wireframe", color="black")
+plotter.view_xy()
+if not pyvista.OFF_SCREEN:
+    plotter.show()
+# -
+
+# ## Solving the eigenvalue problem
+# In this section we will solve the eigenvalue problem:
+#
+# Find $u_h\in H_0^1(\Omega)$ and $\lambda\in\mathbb{R}$ such that
+#
+# $$
+# \begin{align}
+# \int_\Omega \nabla u \cdot \nabla v~\mathrm{d} x &= \lambda \int_\Omega u v~\mathrm{d} x \qquad
+# \forall v \in H_0^1(\Omega).
+# \end{align}
+# $$
+
+# Next, we define a convenience function to solve the eigenvalue problem using [SLEPc](https://slepc.upv.es/)
+# given a discretized domain, its facet markers and the region map.
+
+
+def solve(
+    mesh: dolfinx.mesh.Mesh,
+    facet_tags: dolfinx.mesh.MeshTags,
+    region_map: dict[tuple[int, str], tuple[int, ...]],
+) -> tuple[float, dolfinx.fem.Function, dolfinx.fem.Function]:
+    # We define the lhs and rhs bilinear forms
+    V = dolfinx.fem.functionspace(mesh, ("Lagrange", 3))
+    u = ufl.TrialFunction(V)
+    v = ufl.TestFunction(V)
+    a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
+    m = ufl.inner(u, v) * ufl.dx
+
+    # We identify the boundary facets and their corresponding dofs
+    straight_facets = facet_tags.indices[
+        np.isin(facet_tags.values, region_map[(1, "line")])
+    ]
+    curved_facets = facet_tags.indices[
+        np.isin(facet_tags.values, region_map[(1, "curve")])
+    ]
+    boundary_facets = np.concatenate([straight_facets, curved_facets])
+    mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
+    boundary_dofs = dolfinx.fem.locate_dofs_topological(
+        V, mesh.topology.dim - 1, boundary_facets
+    )
+
+    # We create a zero boundary condition for these dofs to be in the suitable space, and
+    # set up the discrete matrices `A` and `M`
+    bc = dolfinx.fem.dirichletbc(0.0, boundary_dofs, V)
+    A = dolfinx.fem.petsc.assemble_matrix(dolfinx.fem.form(a), bcs=[bc])
+    A.assemble()
+    if Version(dolfinx.__version__) < Version("0.10.0"):
+        diag_kwargs = {"diagonal": 0.0}
+    else:
+        diag_kwargs = {"diag": 0.0}
+
+    M = dolfinx.fem.petsc.assemble_matrix(dolfinx.fem.form(m), bcs=[bc], **diag_kwargs)
+    M.assemble()
+
+    # Next, we define the SLEPc Eigenvalue Problem Solver (EPS), and set up to use a shift
+    # and invert (SINVERT) spectral transformation where the preconditioner factorisation
+    # is computed using [MUMPS](https://mumps-solver.org/index.php).
+
+    E = SLEPc.EPS().create(mesh.comm)
+    E.setType(SLEPc.EPS.Type.ARNOLDI)
+    E.setProblemType(SLEPc.EPS.ProblemType.GHEP)
+    E.setDimensions(1, SLEPc.DECIDE)
+    E.setOperators(A, M)
+    ST = E.getST()
+    ST.setType(SLEPc.ST.Type.SINVERT)
+    PC = ST.getKSP().getPC()
+    PC.setType("lu")
+    PC.setFactorSolverType("mumps")
+    E.setST(ST)
+    E.solve()
+    assert E.getConvergedReason() >= 0, "Eigenvalue solver did not converge"
+
+    # We get the real and imaginary parts of the first eigenvector along with the eigenvalue.
+    uh_r = dolfinx.fem.Function(V)
+    uh_i = dolfinx.fem.Function(V)
+    lam = E.getEigenpair(0, uh_r.x.petsc_vec, uh_i.x.petsc_vec)
+    E.destroy()
+    uh_r.x.scatter_forward()
+    uh_i.x.scatter_forward()
+    return (lam, uh_r, uh_i)
+
+
+# ## Error-indicator
+# In this example, we will use an error-indicator $\eta$ to decide what cells should be refined.
+# Specifically, the estimator is:
+#
+# \begin{align*}
+#  \eta = \sum_{K\in \mathcal{T}_h(\Omega)}\left(h^2\int_K \vert \lambda u_h + \Delta u_h\vert^2~\mathrm{d}x\right)
+# + \sum_{E\in\mathcal{F}_i}\frac{h}{2} \vert [\nabla \cdot \mathbf{n}_E ]\vert^2~\mathrm{d}s
+# \end{align*}
+#
+# where $\mathcal{T}_h$ is the collection of cells in the mesh, $\mathcal{F}_i$ the collection
+
+
+def mark_cells(uh_r: dolfinx.fem.Function, lam: float):
+    mesh = uh_r.function_space.mesh
+    W = dolfinx.fem.functionspace(mesh, ("DG", 0))
+    w = ufl.TestFunction(W)
+    eta = dolfinx.fem.Function(W)
+    f = dolfinx.fem.Constant(mesh, 1.0)
+    h = dolfinx.fem.Function(W)
+    h.x.array[:] = mesh.h(mesh.topology.dim, np.arange(len(h.x.array), dtype=np.int32))
+    n = ufl.FacetNormal(mesh)
+
+    G = (  # compute cellwise error estimator
+        ufl.inner(h**2 * (f + ufl.div(ufl.grad(uh_r))) ** 2, w) * ufl.dx
+        + ufl.inner(h("+") / 2 * ufl.jump(ufl.grad(uh_r), n) ** 2, w("+")) * ufl.dS
+        + ufl.inner(h("-") / 2 * ufl.jump(ufl.grad(uh_r), n) ** 2, w("-")) * ufl.dS
+    )
+    dolfinx.fem.petsc.assemble_vector(eta.x.petsc_vec, dolfinx.fem.form(G))
+    sqrt_eta = dolfinx.fem.Function(W)
+    sqrt_eta.x.array[:] = np.sqrt(eta.x.array[:])
+
+    sqrt_eta_max = sqrt_eta.x.petsc_vec.max()[1]
+
+    theta = 0.5
+    should_refine = ufl.conditional(ufl.gt(sqrt_eta, theta * sqrt_eta_max), 1, 0)
+    markers = dolfinx.fem.Function(W)
+    ip = W.element.interpolation_points
+    if Version(dolfinx.__version__) < Version("0.10.0"):
+        ip = ip()
+    markers.interpolate(dolfinx.fem.Expression(should_refine, ip))
+    return np.flatnonzero(np.isclose(markers.x.array.astype(np.int32), 1))
+
+
+# ## Running the adaptive refinement algorithm
+# Next, we will run the adaptive mesh refinement algorithm.
+
+# We will track the progress of the adaptive mesh refinement as a GIF.
+
+plotter = pyvista.Plotter()
+plotter.open_gif("amr.gif", fps=1)
+
+# We make a convenience function to attach the relevant data to the plotter at a given
+# refinement step.
+
+
+# + tags=["hide-input"]
+def write_frame(plotter: pyvista.Plotter, uh_r: dolfinx.fem.Function):
+    # Scale uh_r to be consistent between refinement steps, as it can be multiplied by -1
+    uh_r_min = curved_mesh.comm.allreduce(uh_r.x.array.min(), op=MPI.MIN)
+    uh_r_max = curved_mesh.comm.allreduce(uh_r.x.array.max(), op=MPI.MAX)
+    uh_sign = np.sign(uh_r_min)
+    if np.isclose(uh_sign, 0):
+        uh_sign = np.sign(uh_r_max)
+    assert not np.isclose(uh_sign, 0), "uh_r has zero values, cannot determine sign."
+    uh_r.x.array[:] *= uh_sign
+
+    # Update plot with refined mesh
+    grid = pyvista.UnstructuredGrid(*dolfinx.plot.vtk_mesh(mesh))
+    curved_grid = pyvista.UnstructuredGrid(*dolfinx.plot.vtk_mesh(uh_r.function_space))
+    curved_grid.point_data["u"] = uh_r.x.array
+    curved_grid = curved_grid.tessellate()
+    curved_actor = plotter.add_mesh(
+        curved_grid,
+        show_edges=False,
+    )
+
+    actor = plotter.add_mesh(grid, style="wireframe", color="black")
+    plotter.view_xy()
+    plotter.write_frame()
+    plotter.remove_actor(actor)
+    plotter.remove_actor(curved_actor)
+
+
+# -
+
+# We set some parameters for checking convergence of the algorithm, and provide the exact eigenvalue
+# for comparison.
+# ```{admonition} Using ngsPETSc for mesh refinement
+# In `ngsPETSc`, we provide the function `GeometricModel.refineMarkedElements` which we
+# pass the entities we would like to refine, and the topological dimensions of those entities.
+# The function returns a refined mesh, with corresponding cell and facet markers extracted from
+# the NetGen model.
+# ```
+
+# + tags=["scroll-output"]
+max_iterations = 15
+exact = 3.375610652693620492628**2
+termination_criteria = 1e-5
+for i in range(max_iterations):
+    lam, uh_r, _ = solve(curved_mesh, ft, region_map)
+
+    relative_error = (lam - exact) / abs(exact)
+    PETSc.Sys.Print(
+        f"Iteration {i + 1}/{max_iterations}, {lam=:.5e}, {exact=:.5e}, {relative_error=:.2e}"
+    )
+
+    cells_to_mark = mark_cells(uh_r, lam)
+    mesh, (_, ft) = geoModel.refineMarkedElements(mesh.topology.dim, cells_to_mark)
+    curved_mesh = geoModel.curveField(order)
+    write_frame(plotter, uh_r)
+
+    if relative_error < termination_criteria:
+        PETSc.Sys.Print(f"Converged in {i + 1} iterations.")
+        break
+plotter.close()
+# -
+
+# <img src="./amr.gif" alt="gif" class="bg-primary mb-1" width="800px">
diff --git a/docker/Dockerfile b/docker/Dockerfile
index b6ff1bd0..eb728854 100644
--- a/docker/Dockerfile
+++ b/docker/Dockerfile
@@ -16,22 +16,32 @@ RUN curl -sL https://deb.nodesource.com/setup_18.x -o nodesource_setup.sh && \
     bash nodesource_setup.sh && \
     apt -y install nodejs
 
-# Upgrade setuptools and pip
-RUN python3 -m pip install -U setuptools pip pkgconfig
-
-# BUILD VTK
-# ENV VTK_VERSION="v9.2.6"
-# RUN git clone --recursive --branch ${VTK_VERSION} --single-branch https://gitlab.kitware.com/vtk/vtk.git vtk && \
-#     cmake -G Ninja -DVTK_WHEEL_BUILD=ON -DVTK_WRAP_PYTHON=ON vtk/ && \
-#     ninja && \
-#     python3 setup.py bdist_wheela
-
-RUN echo ${TARGETPLATFORM}
-RUN if [ "$TARGETPLATFORM" = "linux/arm64" ]; then python3 -m pip install "https://github.com/finsberg/vtk-aarch64/releases/download/vtk-9.3.0-cp312/vtk-9.3.0.dev0-cp312-cp312-linux_aarch64.whl"; fi
-RUN if [ "$TARGETPLATFORM" = "linux/amd64" ]; then python3 -m pip install vtk; fi
+# Install netgen from source
+RUN  apt-get update && \
+ apt-get -y install python3 python3-tk libpython3-dev libxmu-dev tk-dev tcl-dev cmake git g++ libglu1-mesa-dev liblapacke-dev libocct-data-exchange-dev libocct-draw-dev occt-misc libtbb-dev libxi-dev
+WORKDIR /ngsuite
+RUN git clone --branch=v6.2.2505 --single-branch https://github.com/NGSolve/netgen.git netgen-src
+WORKDIR /ngsuite/netgen-src
+RUN git submodule update --init --recursive
+WORKDIR /ngsuite
+RUN mkdir netgen-build
+RUN mkdir netgen-install
+RUN cmake -DCMAKE="-cxx-flags=-flax-vector-conversions" -DCMAKE_INSTALL_PREFIX=/ngsuite/netgen-install /ngsuite/netgen-src
+RUN make
+RUN make install
+ENV NETGENDIR=/ngsuite/netgen-install/bin
+ENV PATH=$NETGENDIR:$PATH
+ENV PYTHONPATH=$NETGENDIR/../lib/python3.12/site-packages:${PYTHONPATH}:${PYTHONPATH}:$PATH
+
+# Install vtk
+RUN python3 -m pip install vtk
 
 ADD pyproject.toml /tmp/pyproject.toml
-RUN python3 -m pip install --no-cache-dir --no-binary=h5py -v .
+WORKDIR /tmp
+RUN if [ "$TARGETPLATFORM" = "linux/arm64" ]; then python3 -m pip install --no-cache-dir --no-binary=h5py -v .; fi
+RUN if [ "$TARGETPLATFORM" = "linux/amd64" ]; then python3 -m pip install --no-cache-dir --no-binary=h5py -v .[netgen]; fi
+
+RUN python3 -m pip install -v ngsPETSc@git+https://github.com/NGSolve/ngsPETSc.git scipy --no-deps
 RUN python3 -m pip cache purge
 
 ENTRYPOINT ["jupyter", "lab", "--ip", "0.0.0.0", "--no-browser", "--allow-root"]
diff --git a/pyproject.toml b/pyproject.toml
index 9b898e16..5a894e06 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,5 +1,5 @@
 [build-system]
-requires = ["setuptools>=64.4.0", "wheel", "pip>=22.3"]
+requires = ["setuptools>=64.4.0", "wheel", "pip>=22.3", "poetry-core"]
 build-backend = "setuptools.build_meta"
 
 [project]
@@ -18,6 +18,9 @@ dependencies = [
 
 [project.optional-dependencies]
 dev = ["pdbpp", "ipython", "jupytext", "ruff", "pre-commit"]
+netgen = [
+    "ngsPETSc@git+https://github.com/NGSolve/ngsPETSc.git",
+] # Has to be optional as we cant get netgen on linux/arm64
 
 [tool.setuptools]
 packages = []