From ac896f2cfe5fde906d2c4c8ba23dab9039c9db87 Mon Sep 17 00:00:00 2001
From: RemDelaporteMathurin <rdelaportemathurin@gmail.com>
Date: Tue, 24 Feb 2026 09:37:52 -0500
Subject: [PATCH 1/5] initial example

---
 examples/real_function_space_nonlinear.ipynb | 418 +++++++++++++++++++
 1 file changed, 418 insertions(+)
 create mode 100644 examples/real_function_space_nonlinear.ipynb

diff --git a/examples/real_function_space_nonlinear.ipynb b/examples/real_function_space_nonlinear.ipynb
new file mode 100644
index 0000000..3567567
--- /dev/null
+++ b/examples/real_function_space_nonlinear.ipynb
@@ -0,0 +1,418 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "2424af5e",
+   "metadata": {},
+   "source": [
+    "# Real function spaces (non linear)\n",
+    "\n",
+    "Author: Remi Delaporte-Mathurin, Jørgen S. Dokken\n",
+    "\n",
+    "License: MIT\n",
+    "\n",
+    "In this example we will show how to use the \"real\" function space to solve\n",
+    "a non linear problem.\n",
+    "\n",
+    "## Mathematical formulation\n",
+    " The problem at hand is:\n",
+    "Find $u \\in H^1(\\Omega)$ such that\n",
+    "\n",
+    "$$\n",
+    "\\begin{align}\n",
+    "\\Delta u &= 0 \\quad \\text{in } \\Omega, \\\\\n",
+    "u &= 0 \\quad \\text{on } \\Gamma_1, \\\\\n",
+    "u &= K P_b \\quad \\text{on } \\Gamma_b, \\\\\n",
+    "\\end{align}\n",
+    "$$\n",
+    "\n",
+    "where $P_b$ is the partial pressure inside the cavity region. The temporal evolution of $P_b$ is governed by:\n",
+    "\n",
+    "$$\n",
+    "\\frac{dP_b}{dt} = \\frac{e}{V} \\int_{\\Gamma_b} -D \\nabla c \\cdot \\mathbf{n} dS\n",
+    "$$\n",
+    "\n",
+    "\n",
+    "We write the weak formulation:\n",
+    "\n",
+    "$$\n",
+    "\\begin{align}\n",
+    "\\int_\\Omega  \\frac{u - u_n}{\\Delta t} \\cdot v~\\mathrm{d}x + \\int_\\Omega \\nabla u \\cdot \\nabla v~\\mathrm{d}x &= 0\\\\\n",
+    "\\int_\\Omega  \\frac{P_b - P_{b_n}}{\\Delta t} \\cdot w~\\mathrm{d}x &= \\frac{e}{V} \\int_{\\Gamma_b} -D \\nabla c \\cdot \\mathbf{n} ~dS ~w.\n",
+    "\\end{align}\n",
+    "$$\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "f658dceb",
+   "metadata": {},
+   "source": [
+    "## Implementation\n",
+    "We start by import the necessary modules\n",
+    "```{admonition} Clickable functions/classes\n",
+    "Note that for the modules imported in this example, you can click on the function/class\n",
+    "name to be redirected to the corresponding documentation page.\n",
+    "```"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "06e16e63",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from mpi4py import MPI\n",
+    "import dolfinx\n",
+    "from dolfinx import log\n",
+    "from dolfinx.fem.petsc import NonlinearProblem\n",
+    "from dolfinx.io import VTXWriter, gmsh as gmshio\n",
+    "import numpy as np\n",
+    "from scifem import create_real_functionspace\n",
+    "import ufl\n",
+    "\n",
+    "import gmsh"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "027d8dc8",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Info    : Meshing 1D...                                                                                                                        \n",
+      "Info    : [  0%] Meshing curve 5 (Ellipse)\n",
+      "Info    : [ 30%] Meshing curve 6 (Line)\n",
+      "Info    : [ 50%] Meshing curve 7 (Line)\n",
+      "Info    : [ 70%] Meshing curve 8 (Line)\n",
+      "Info    : [ 90%] Meshing curve 9 (Line)\n",
+      "Info    : Done meshing 1D (Wall 0.00303789s, CPU 0.002975s)\n",
+      "Info    : Meshing 2D...\n",
+      "Info    : Meshing surface 1 (Plane, Frontal-Delaunay for Quads)\n",
+      "Info    : Simple recombination completed (Wall 0.000968826s, CPU 0.000586s): 79 quads, 12 triangles, 0 invalid quads, 0 quads with Q < 0.1, avg Q = 0.791132, min Q = 0.503081\n",
+      "Info    : Simple recombination completed (Wall 0.00135738s, CPU 0.001233s): 352 quads, 0 triangles, 0 invalid quads, 0 quads with Q < 0.1, avg Q = 0.844698, min Q = 0.493962\n",
+      "Info    : Done meshing 2D (Wall 0.00381131s, CPU 0.003215s)\n",
+      "Info    : Refining mesh...\n",
+      "Info    : Meshing order 2 (curvilinear on)...\n",
+      "Info    : [  0%] Meshing curve 5 order 2\n",
+      "Info    : [ 20%] Meshing curve 6 order 2\n",
+      "Info    : [ 40%] Meshing curve 7 order 2\n",
+      "Info    : [ 60%] Meshing curve 8 order 2\n",
+      "Info    : [ 70%] Meshing curve 9 order 2\n",
+      "Info    : [ 90%] Meshing surface 1 order 2\n",
+      "Info    : Done meshing order 2 (Wall 0.0010111s, CPU 0.001076s)\n",
+      "Info    : Done refining mesh (Wall 0.0011461s, CPU 0.00123s)\n",
+      "Info    : 1484 nodes 1565 elements\n",
+      "Info    : Meshing order 2 (curvilinear on)...\n",
+      "Info    : [  0%] Meshing curve 5 order 2\n",
+      "Info    : [ 20%] Meshing curve 6 order 2\n",
+      "Info    : [ 40%] Meshing curve 7 order 2\n",
+      "Info    : [ 60%] Meshing curve 8 order 2\n",
+      "Info    : [ 70%] Meshing curve 9 order 2\n",
+      "Info    : [ 90%] Meshing surface 1 order 2\n",
+      "Info    : Done meshing order 2 (Wall 0.00412756s, CPU 0.003862s)\n",
+      "Info    : Optimizing mesh (Netgen)...\n",
+      "Info    : Done optimizing mesh (Wall 5.89003e-07s, CPU 1e-06s)\n"
+     ]
+    }
+   ],
+   "source": [
+    "gmsh.initialize()\n",
+    "\n",
+    "L = 0.5\n",
+    "H = 0.41\n",
+    "c_x = c_y = 0.2\n",
+    "r = 0.05\n",
+    "gdim = 2\n",
+    "mesh_comm = MPI.COMM_WORLD\n",
+    "model_rank = 0\n",
+    "if mesh_comm.rank == model_rank:\n",
+    "    rectangle = gmsh.model.occ.addRectangle(0, 0, 0, L, H, tag=1)\n",
+    "    obstacle = gmsh.model.occ.addDisk(c_x, c_y, 0, r, r)\n",
+    "\n",
+    "if mesh_comm.rank == model_rank:\n",
+    "    fluid = gmsh.model.occ.cut([(gdim, rectangle)], [(gdim, obstacle)])\n",
+    "    gmsh.model.occ.synchronize()\n",
+    "\n",
+    "solid_marker = 1\n",
+    "if mesh_comm.rank == model_rank:\n",
+    "    volumes = gmsh.model.getEntities(dim=gdim)\n",
+    "    assert len(volumes) == 1\n",
+    "    gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], solid_marker)\n",
+    "    gmsh.model.setPhysicalName(volumes[0][0], solid_marker, \"Solid\")\n",
+    "\n",
+    "wall_marker, obstacle_marker = 2, 3\n",
+    "walls, obstacle = [], []\n",
+    "if mesh_comm.rank == model_rank:\n",
+    "    boundaries = gmsh.model.getBoundary(volumes, oriented=False)\n",
+    "    for boundary in boundaries:\n",
+    "        center_of_mass = gmsh.model.occ.getCenterOfMass(boundary[0], boundary[1])\n",
+    "        if np.allclose(center_of_mass, [0, H / 2, 0]):\n",
+    "            walls.append(boundary[1])\n",
+    "        elif np.allclose(center_of_mass, [L, H / 2, 0]):\n",
+    "            walls.append(boundary[1])\n",
+    "        elif np.allclose(center_of_mass, [L / 2, H, 0]) or np.allclose(\n",
+    "            center_of_mass, [L / 2, 0, 0]\n",
+    "        ):\n",
+    "            walls.append(boundary[1])\n",
+    "        else:\n",
+    "            obstacle.append(boundary[1])\n",
+    "    gmsh.model.addPhysicalGroup(1, walls, wall_marker)\n",
+    "    gmsh.model.setPhysicalName(1, wall_marker, \"Walls\")\n",
+    "    gmsh.model.addPhysicalGroup(1, obstacle, obstacle_marker)\n",
+    "    gmsh.model.setPhysicalName(1, obstacle_marker, \"Obstacle\")\n",
+    "\n",
+    "\n",
+    "res_min = r / 3\n",
+    "if mesh_comm.rank == model_rank:\n",
+    "    distance_field = gmsh.model.mesh.field.add(\"Distance\")\n",
+    "    gmsh.model.mesh.field.setNumbers(distance_field, \"EdgesList\", obstacle)\n",
+    "    threshold_field = gmsh.model.mesh.field.add(\"Threshold\")\n",
+    "    gmsh.model.mesh.field.setNumber(threshold_field, \"IField\", distance_field)\n",
+    "    gmsh.model.mesh.field.setNumber(threshold_field, \"LcMin\", res_min)\n",
+    "    gmsh.model.mesh.field.setNumber(threshold_field, \"LcMax\", 0.25 * H)\n",
+    "    gmsh.model.mesh.field.setNumber(threshold_field, \"DistMin\", r)\n",
+    "    gmsh.model.mesh.field.setNumber(threshold_field, \"DistMax\", 2 * H)\n",
+    "    min_field = gmsh.model.mesh.field.add(\"Min\")\n",
+    "    gmsh.model.mesh.field.setNumbers(min_field, \"FieldsList\", [threshold_field])\n",
+    "    gmsh.model.mesh.field.setAsBackgroundMesh(min_field)\n",
+    "\n",
+    "\n",
+    "if mesh_comm.rank == model_rank:\n",
+    "    gmsh.option.setNumber(\"Mesh.Algorithm\", 8)\n",
+    "    gmsh.option.setNumber(\"Mesh.RecombinationAlgorithm\", 2)\n",
+    "    gmsh.option.setNumber(\"Mesh.RecombineAll\", 1)\n",
+    "    gmsh.option.setNumber(\"Mesh.SubdivisionAlgorithm\", 1)\n",
+    "    gmsh.model.mesh.generate(gdim)\n",
+    "    gmsh.model.mesh.setOrder(2)\n",
+    "    gmsh.model.mesh.optimize(\"Netgen\")\n",
+    "\n",
+    "mesh_data = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)\n",
+    "mesh = mesh_data.mesh\n",
+    "assert mesh_data.facet_tags is not None\n",
+    "ft = mesh_data.facet_tags\n",
+    "ft.name = \"Facet markers\""
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "id": "a5039e02",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "06b10b43b4034fedb3a0312c5264d6ca",
+       "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"
+    },
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "141f00aa54dc4045a61869efa01fcaf0",
+       "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": [
+    "from dolfinx import plot\n",
+    "import pyvista\n",
+    "\n",
+    "pyvista.set_jupyter_backend(\"html\")\n",
+    "\n",
+    "tdim = mesh.topology.dim\n",
+    "\n",
+    "mesh.topology.create_connectivity(tdim, tdim)\n",
+    "topology, cell_types, geometry = plot.vtk_mesh(mesh, tdim)\n",
+    "grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)\n",
+    "\n",
+    "plotter = pyvista.Plotter()\n",
+    "plotter.add_mesh(grid, show_edges=True)\n",
+    "plotter.view_xy()\n",
+    "if not pyvista.OFF_SCREEN:\n",
+    "    plotter.show()\n",
+    "else:\n",
+    "    figure = plotter.screenshot(\"mesh.png\")\n",
+    "\n",
+    "fdim = tdim - 1\n",
+    "mesh.topology.create_connectivity(fdim, tdim)\n",
+    "topology, cell_types, x = plot.vtk_mesh(mesh, fdim, ft.indices)\n",
+    "\n",
+    "p = pyvista.Plotter()\n",
+    "grid = pyvista.UnstructuredGrid(topology, cell_types, x)\n",
+    "grid.cell_data[\"Facet Marker\"] = ft.values\n",
+    "grid.set_active_scalars(\"Facet Marker\")\n",
+    "p.view_xy()\n",
+    "p.add_mesh(grid, show_edges=True)\n",
+    "\n",
+    "if not pyvista.OFF_SCREEN:\n",
+    "    p.show()\n",
+    "else:\n",
+    "    figure = p.screenshot(\"facet_markers.png\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7acbb352",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "<dolfinx.mesh.Mesh object at 0x75c8a2c5da90>\n"
+     ]
+    },
+    {
+     "ename": "TypeError",
+     "evalue": "create_real_functionspace_float64(): incompatible function arguments. The following argument types are supported:\n    1. create_real_functionspace_float64(arg0: dolfinx::mesh::Mesh<double>, arg1: collections.abc.Sequence[int], /) -> dolfinx::fem::FunctionSpace<double>\n\nInvoked with types: dolfinx.cpp.mesh.Mesh_float64, tuple",
+     "output_type": "error",
+     "traceback": [
+      "\u001b[31m---------------------------------------------------------------------------\u001b[39m",
+      "\u001b[31mTypeError\u001b[39m                                 Traceback (most recent call last)",
+      "\u001b[36mCell\u001b[39m\u001b[36m \u001b[39m\u001b[32mIn[9]\u001b[39m\u001b[32m, line 3\u001b[39m\n\u001b[32m      1\u001b[39m V = dolfinx.fem.functionspace(mesh, (\u001b[33m\"\u001b[39m\u001b[33mLagrange\u001b[39m\u001b[33m\"\u001b[39m, \u001b[32m1\u001b[39m))\n\u001b[32m      2\u001b[39m \u001b[38;5;28mprint\u001b[39m(mesh)\n\u001b[32m----> \u001b[39m\u001b[32m3\u001b[39m R = \u001b[43mcreate_real_functionspace\u001b[49m\u001b[43m(\u001b[49m\u001b[43mmesh\u001b[49m\u001b[43m)\u001b[49m\n\u001b[32m      5\u001b[39m W = ufl.MixedFunctionSpace(V, R)\n",
+      "\u001b[36mFile \u001b[39m\u001b[32m~/scifem/src/scifem/spaces.py:32\u001b[39m, in \u001b[36mcreate_real_functionspace\u001b[39m\u001b[34m(mesh, value_shape)\u001b[39m\n\u001b[32m     27\u001b[39m ufl_e = basix.ufl.element(\n\u001b[32m     28\u001b[39m     \u001b[33m\"\u001b[39m\u001b[33mP\u001b[39m\u001b[33m\"\u001b[39m, mesh.basix_cell(), \u001b[32m0\u001b[39m, dtype=dtype, discontinuous=\u001b[38;5;28;01mTrue\u001b[39;00m, shape=value_shape\n\u001b[32m     29\u001b[39m )\n\u001b[32m     31\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m (dtype := mesh.geometry.x.dtype) == np.float64:\n\u001b[32m---> \u001b[39m\u001b[32m32\u001b[39m     cppV = \u001b[43m_scifem\u001b[49m\u001b[43m.\u001b[49m\u001b[43mcreate_real_functionspace_float64\u001b[49m\u001b[43m(\u001b[49m\u001b[43mmesh\u001b[49m\u001b[43m.\u001b[49m\u001b[43m_cpp_object\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mvalue_shape\u001b[49m\u001b[43m)\u001b[49m\n\u001b[32m     33\u001b[39m \u001b[38;5;28;01melif\u001b[39;00m dtype == np.float32:\n\u001b[32m     34\u001b[39m     cppV = _scifem.create_real_functionspace_float32(mesh._cpp_object, value_shape)\n",
+      "\u001b[31mTypeError\u001b[39m: create_real_functionspace_float64(): incompatible function arguments. The following argument types are supported:\n    1. create_real_functionspace_float64(arg0: dolfinx::mesh::Mesh<double>, arg1: collections.abc.Sequence[int], /) -> dolfinx::fem::FunctionSpace<double>\n\nInvoked with types: dolfinx.cpp.mesh.Mesh_float64, tuple"
+     ]
+    }
+   ],
+   "source": [
+    "V = dolfinx.fem.functionspace(mesh, (\"Lagrange\", 1))\n",
+    "R = create_real_functionspace(mesh)\n",
+    "\n",
+    "W = ufl.MixedFunctionSpace(V, R)\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f9d9b61c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\n",
+    "\n",
+    "dofs_boundary = dolfinx.fem.locate_dofs_topological(V, fdim, ft.indices[ft.values == wall_marker])\n",
+    "dofs_bubble = dolfinx.fem.locate_dofs_topological(V, fdim, ft.indices[ft.values == obstacle_marker]\n",
+    ")\n",
+    "\n",
+    "\n",
+    "u = dolfinx.fem.Function(V)\n",
+    "u_n = dolfinx.fem.Function(V)\n",
+    "pressure = dolfinx.fem.Function(R)\n",
+    "pressure_n = dolfinx.fem.Function(R)\n",
+    "\n",
+    "dt = dolfinx.fem.Constant(mesh, 1.0)\n",
+    "K_S = dolfinx.fem.Constant(mesh, 1.0)\n",
+    "e = dolfinx.fem.Constant(mesh, 2.0)\n",
+    "volume = dolfinx.fem.Constant(mesh, 4.0)\n",
+    "\n",
+    "v, pressure_v = ufl.TestFunctions(W)\n",
+    "\n",
+    "pressure_ini_expr = dolfinx.fem.Expression(\n",
+    "    dolfinx.fem.Constant(mesh, 3.5), R.element.interpolation_points\n",
+    ")\n",
+    "pressure_n.interpolate(pressure_ini_expr)\n",
+    "pressure.x.array[:] = pressure_n.x.array[:]\n",
+    "\n",
+    "bc_bubble_expr = dolfinx.fem.Expression(K_S * pressure, V.element.interpolation_points)\n",
+    "u_bc_bubble = dolfinx.fem.Function(V)\n",
+    "u_bc_bubble.interpolate(bc_bubble_expr)\n",
+    "\n",
+    "bc_bubble = dolfinx.fem.dirichletbc(u_bc_bubble, dofs_bubble)\n",
+    "bc_boundary = dolfinx.fem.dirichletbc(dolfinx.fem.Constant(mesh, 0.0), dofs_boundary, V)\n",
+    "\n",
+    "n = ufl.FacetNormal(mesh)\n",
+    "flux = ufl.inner(ufl.grad(u), n)\n",
+    "\n",
+    "ds = ufl.Measure(\"ds\", domain=mesh, subdomain_data=ft)\n",
+    "\n",
+    "F = -(u - u_n) / dt * v * ufl.dx + ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx\n",
+    "F += (pressure - pressure_n) / dt * pressure_v * ufl.dx - e / volume * flux * pressure_v * ds(2)\n",
+    "\n",
+    "forms = ufl.extract_blocks(F)\n",
+    "\n",
+    "solver = NonlinearProblem(\n",
+    "    forms,\n",
+    "    [u, pressure],\n",
+    "    bcs=[bc_boundary, bc_bubble],\n",
+    "    petsc_options_prefix=\"bubble\",\n",
+    "    petsc_options={\n",
+    "        \"ksp_type\": \"preonly\",\n",
+    "        \"pc_type\": \"lu\",\n",
+    "        \"pc_factor_mat_solver_type\": \"mumps\",\n",
+    "        \"snes_monitor\": None,\n",
+    "    },\n",
+    ")\n",
+    "\n",
+    "log.set_log_level(log.LogLevel.INFO)\n",
+    "t = 0\n",
+    "t_final = 20\n",
+    "times = []\n",
+    "\n",
+    "all_pressures = []\n",
+    "while t < t_final:\n",
+    "    t += dt.value\n",
+    "    times.append(t)\n",
+    "    print(f\"Solving at time {t:.2f}\")\n",
+    "\n",
+    "    # Solve the problem\n",
+    "    (u, pressure) = solver.solve()\n",
+    "    u_n.x.array[:] = u.x.array[:]\n",
+    "    pressure_n.x.array[:] = pressure.x.array[:]\n",
+    "    print(f\"u = {u.x.array}\")\n",
+    "    print(f\"pressure = {pressure.x.array}\")\n",
+    "\n",
+    "    u_bc_bubble.interpolate(bc_bubble_expr)\n",
+    "    all_pressures.append(pressure.x.array[0].copy())"
+   ]
+  }
+ ],
+ "metadata": {
+  "jupytext": {
+   "cell_metadata_filter": "-all",
+   "formats": "py,ipynb,md:myst",
+   "main_language": "python"
+  },
+  "kernelspec": {
+   "display_name": "scifem-dev",
+   "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.14.3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

From 289fab750e64cf5c5a0bdc274f9dc0297bfe37c1 Mon Sep 17 00:00:00 2001
From: RemDelaporteMathurin <rdelaportemathurin@gmail.com>
Date: Tue, 24 Feb 2026 11:26:03 -0500
Subject: [PATCH 2/5] added python file

---
 _toc.yml                                     |   1 +
 examples/real_function_space_nonlinear.ipynb | 418 -------------------
 examples/real_function_space_nonlinear.py    | 319 ++++++++++++++
 3 files changed, 320 insertions(+), 418 deletions(-)
 delete mode 100644 examples/real_function_space_nonlinear.ipynb
 create mode 100644 examples/real_function_space_nonlinear.py

diff --git a/_toc.yml b/_toc.yml
index cb5368b..b18d50c 100644
--- a/_toc.yml
+++ b/_toc.yml
@@ -7,6 +7,7 @@ parts:
   - caption: Features
     chapters:
       - file: "examples/real_function_space.py"
+      - file: "examples/real_function_space_nonlinear.py"
       - file: "examples/point_source.py"
       - file: "examples/xdmf_point_cloud.py"
       - file: "examples/mark_entities.py"
diff --git a/examples/real_function_space_nonlinear.ipynb b/examples/real_function_space_nonlinear.ipynb
deleted file mode 100644
index 3567567..0000000
--- a/examples/real_function_space_nonlinear.ipynb
+++ /dev/null
@@ -1,418 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "markdown",
-   "id": "2424af5e",
-   "metadata": {},
-   "source": [
-    "# Real function spaces (non linear)\n",
-    "\n",
-    "Author: Remi Delaporte-Mathurin, Jørgen S. Dokken\n",
-    "\n",
-    "License: MIT\n",
-    "\n",
-    "In this example we will show how to use the \"real\" function space to solve\n",
-    "a non linear problem.\n",
-    "\n",
-    "## Mathematical formulation\n",
-    " The problem at hand is:\n",
-    "Find $u \\in H^1(\\Omega)$ such that\n",
-    "\n",
-    "$$\n",
-    "\\begin{align}\n",
-    "\\Delta u &= 0 \\quad \\text{in } \\Omega, \\\\\n",
-    "u &= 0 \\quad \\text{on } \\Gamma_1, \\\\\n",
-    "u &= K P_b \\quad \\text{on } \\Gamma_b, \\\\\n",
-    "\\end{align}\n",
-    "$$\n",
-    "\n",
-    "where $P_b$ is the partial pressure inside the cavity region. The temporal evolution of $P_b$ is governed by:\n",
-    "\n",
-    "$$\n",
-    "\\frac{dP_b}{dt} = \\frac{e}{V} \\int_{\\Gamma_b} -D \\nabla c \\cdot \\mathbf{n} dS\n",
-    "$$\n",
-    "\n",
-    "\n",
-    "We write the weak formulation:\n",
-    "\n",
-    "$$\n",
-    "\\begin{align}\n",
-    "\\int_\\Omega  \\frac{u - u_n}{\\Delta t} \\cdot v~\\mathrm{d}x + \\int_\\Omega \\nabla u \\cdot \\nabla v~\\mathrm{d}x &= 0\\\\\n",
-    "\\int_\\Omega  \\frac{P_b - P_{b_n}}{\\Delta t} \\cdot w~\\mathrm{d}x &= \\frac{e}{V} \\int_{\\Gamma_b} -D \\nabla c \\cdot \\mathbf{n} ~dS ~w.\n",
-    "\\end{align}\n",
-    "$$\n"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "f658dceb",
-   "metadata": {},
-   "source": [
-    "## Implementation\n",
-    "We start by import the necessary modules\n",
-    "```{admonition} Clickable functions/classes\n",
-    "Note that for the modules imported in this example, you can click on the function/class\n",
-    "name to be redirected to the corresponding documentation page.\n",
-    "```"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 1,
-   "id": "06e16e63",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "from mpi4py import MPI\n",
-    "import dolfinx\n",
-    "from dolfinx import log\n",
-    "from dolfinx.fem.petsc import NonlinearProblem\n",
-    "from dolfinx.io import VTXWriter, gmsh as gmshio\n",
-    "import numpy as np\n",
-    "from scifem import create_real_functionspace\n",
-    "import ufl\n",
-    "\n",
-    "import gmsh"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 2,
-   "id": "027d8dc8",
-   "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Info    : Meshing 1D...                                                                                                                        \n",
-      "Info    : [  0%] Meshing curve 5 (Ellipse)\n",
-      "Info    : [ 30%] Meshing curve 6 (Line)\n",
-      "Info    : [ 50%] Meshing curve 7 (Line)\n",
-      "Info    : [ 70%] Meshing curve 8 (Line)\n",
-      "Info    : [ 90%] Meshing curve 9 (Line)\n",
-      "Info    : Done meshing 1D (Wall 0.00303789s, CPU 0.002975s)\n",
-      "Info    : Meshing 2D...\n",
-      "Info    : Meshing surface 1 (Plane, Frontal-Delaunay for Quads)\n",
-      "Info    : Simple recombination completed (Wall 0.000968826s, CPU 0.000586s): 79 quads, 12 triangles, 0 invalid quads, 0 quads with Q < 0.1, avg Q = 0.791132, min Q = 0.503081\n",
-      "Info    : Simple recombination completed (Wall 0.00135738s, CPU 0.001233s): 352 quads, 0 triangles, 0 invalid quads, 0 quads with Q < 0.1, avg Q = 0.844698, min Q = 0.493962\n",
-      "Info    : Done meshing 2D (Wall 0.00381131s, CPU 0.003215s)\n",
-      "Info    : Refining mesh...\n",
-      "Info    : Meshing order 2 (curvilinear on)...\n",
-      "Info    : [  0%] Meshing curve 5 order 2\n",
-      "Info    : [ 20%] Meshing curve 6 order 2\n",
-      "Info    : [ 40%] Meshing curve 7 order 2\n",
-      "Info    : [ 60%] Meshing curve 8 order 2\n",
-      "Info    : [ 70%] Meshing curve 9 order 2\n",
-      "Info    : [ 90%] Meshing surface 1 order 2\n",
-      "Info    : Done meshing order 2 (Wall 0.0010111s, CPU 0.001076s)\n",
-      "Info    : Done refining mesh (Wall 0.0011461s, CPU 0.00123s)\n",
-      "Info    : 1484 nodes 1565 elements\n",
-      "Info    : Meshing order 2 (curvilinear on)...\n",
-      "Info    : [  0%] Meshing curve 5 order 2\n",
-      "Info    : [ 20%] Meshing curve 6 order 2\n",
-      "Info    : [ 40%] Meshing curve 7 order 2\n",
-      "Info    : [ 60%] Meshing curve 8 order 2\n",
-      "Info    : [ 70%] Meshing curve 9 order 2\n",
-      "Info    : [ 90%] Meshing surface 1 order 2\n",
-      "Info    : Done meshing order 2 (Wall 0.00412756s, CPU 0.003862s)\n",
-      "Info    : Optimizing mesh (Netgen)...\n",
-      "Info    : Done optimizing mesh (Wall 5.89003e-07s, CPU 1e-06s)\n"
-     ]
-    }
-   ],
-   "source": [
-    "gmsh.initialize()\n",
-    "\n",
-    "L = 0.5\n",
-    "H = 0.41\n",
-    "c_x = c_y = 0.2\n",
-    "r = 0.05\n",
-    "gdim = 2\n",
-    "mesh_comm = MPI.COMM_WORLD\n",
-    "model_rank = 0\n",
-    "if mesh_comm.rank == model_rank:\n",
-    "    rectangle = gmsh.model.occ.addRectangle(0, 0, 0, L, H, tag=1)\n",
-    "    obstacle = gmsh.model.occ.addDisk(c_x, c_y, 0, r, r)\n",
-    "\n",
-    "if mesh_comm.rank == model_rank:\n",
-    "    fluid = gmsh.model.occ.cut([(gdim, rectangle)], [(gdim, obstacle)])\n",
-    "    gmsh.model.occ.synchronize()\n",
-    "\n",
-    "solid_marker = 1\n",
-    "if mesh_comm.rank == model_rank:\n",
-    "    volumes = gmsh.model.getEntities(dim=gdim)\n",
-    "    assert len(volumes) == 1\n",
-    "    gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], solid_marker)\n",
-    "    gmsh.model.setPhysicalName(volumes[0][0], solid_marker, \"Solid\")\n",
-    "\n",
-    "wall_marker, obstacle_marker = 2, 3\n",
-    "walls, obstacle = [], []\n",
-    "if mesh_comm.rank == model_rank:\n",
-    "    boundaries = gmsh.model.getBoundary(volumes, oriented=False)\n",
-    "    for boundary in boundaries:\n",
-    "        center_of_mass = gmsh.model.occ.getCenterOfMass(boundary[0], boundary[1])\n",
-    "        if np.allclose(center_of_mass, [0, H / 2, 0]):\n",
-    "            walls.append(boundary[1])\n",
-    "        elif np.allclose(center_of_mass, [L, H / 2, 0]):\n",
-    "            walls.append(boundary[1])\n",
-    "        elif np.allclose(center_of_mass, [L / 2, H, 0]) or np.allclose(\n",
-    "            center_of_mass, [L / 2, 0, 0]\n",
-    "        ):\n",
-    "            walls.append(boundary[1])\n",
-    "        else:\n",
-    "            obstacle.append(boundary[1])\n",
-    "    gmsh.model.addPhysicalGroup(1, walls, wall_marker)\n",
-    "    gmsh.model.setPhysicalName(1, wall_marker, \"Walls\")\n",
-    "    gmsh.model.addPhysicalGroup(1, obstacle, obstacle_marker)\n",
-    "    gmsh.model.setPhysicalName(1, obstacle_marker, \"Obstacle\")\n",
-    "\n",
-    "\n",
-    "res_min = r / 3\n",
-    "if mesh_comm.rank == model_rank:\n",
-    "    distance_field = gmsh.model.mesh.field.add(\"Distance\")\n",
-    "    gmsh.model.mesh.field.setNumbers(distance_field, \"EdgesList\", obstacle)\n",
-    "    threshold_field = gmsh.model.mesh.field.add(\"Threshold\")\n",
-    "    gmsh.model.mesh.field.setNumber(threshold_field, \"IField\", distance_field)\n",
-    "    gmsh.model.mesh.field.setNumber(threshold_field, \"LcMin\", res_min)\n",
-    "    gmsh.model.mesh.field.setNumber(threshold_field, \"LcMax\", 0.25 * H)\n",
-    "    gmsh.model.mesh.field.setNumber(threshold_field, \"DistMin\", r)\n",
-    "    gmsh.model.mesh.field.setNumber(threshold_field, \"DistMax\", 2 * H)\n",
-    "    min_field = gmsh.model.mesh.field.add(\"Min\")\n",
-    "    gmsh.model.mesh.field.setNumbers(min_field, \"FieldsList\", [threshold_field])\n",
-    "    gmsh.model.mesh.field.setAsBackgroundMesh(min_field)\n",
-    "\n",
-    "\n",
-    "if mesh_comm.rank == model_rank:\n",
-    "    gmsh.option.setNumber(\"Mesh.Algorithm\", 8)\n",
-    "    gmsh.option.setNumber(\"Mesh.RecombinationAlgorithm\", 2)\n",
-    "    gmsh.option.setNumber(\"Mesh.RecombineAll\", 1)\n",
-    "    gmsh.option.setNumber(\"Mesh.SubdivisionAlgorithm\", 1)\n",
-    "    gmsh.model.mesh.generate(gdim)\n",
-    "    gmsh.model.mesh.setOrder(2)\n",
-    "    gmsh.model.mesh.optimize(\"Netgen\")\n",
-    "\n",
-    "mesh_data = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)\n",
-    "mesh = mesh_data.mesh\n",
-    "assert mesh_data.facet_tags is not None\n",
-    "ft = mesh_data.facet_tags\n",
-    "ft.name = \"Facet markers\""
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 6,
-   "id": "a5039e02",
-   "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "application/vnd.jupyter.widget-view+json": {
-       "model_id": "06b10b43b4034fedb3a0312c5264d6ca",
-       "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"
-    },
-    {
-     "data": {
-      "application/vnd.jupyter.widget-view+json": {
-       "model_id": "141f00aa54dc4045a61869efa01fcaf0",
-       "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": [
-    "from dolfinx import plot\n",
-    "import pyvista\n",
-    "\n",
-    "pyvista.set_jupyter_backend(\"html\")\n",
-    "\n",
-    "tdim = mesh.topology.dim\n",
-    "\n",
-    "mesh.topology.create_connectivity(tdim, tdim)\n",
-    "topology, cell_types, geometry = plot.vtk_mesh(mesh, tdim)\n",
-    "grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)\n",
-    "\n",
-    "plotter = pyvista.Plotter()\n",
-    "plotter.add_mesh(grid, show_edges=True)\n",
-    "plotter.view_xy()\n",
-    "if not pyvista.OFF_SCREEN:\n",
-    "    plotter.show()\n",
-    "else:\n",
-    "    figure = plotter.screenshot(\"mesh.png\")\n",
-    "\n",
-    "fdim = tdim - 1\n",
-    "mesh.topology.create_connectivity(fdim, tdim)\n",
-    "topology, cell_types, x = plot.vtk_mesh(mesh, fdim, ft.indices)\n",
-    "\n",
-    "p = pyvista.Plotter()\n",
-    "grid = pyvista.UnstructuredGrid(topology, cell_types, x)\n",
-    "grid.cell_data[\"Facet Marker\"] = ft.values\n",
-    "grid.set_active_scalars(\"Facet Marker\")\n",
-    "p.view_xy()\n",
-    "p.add_mesh(grid, show_edges=True)\n",
-    "\n",
-    "if not pyvista.OFF_SCREEN:\n",
-    "    p.show()\n",
-    "else:\n",
-    "    figure = p.screenshot(\"facet_markers.png\")"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "7acbb352",
-   "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "<dolfinx.mesh.Mesh object at 0x75c8a2c5da90>\n"
-     ]
-    },
-    {
-     "ename": "TypeError",
-     "evalue": "create_real_functionspace_float64(): incompatible function arguments. The following argument types are supported:\n    1. create_real_functionspace_float64(arg0: dolfinx::mesh::Mesh<double>, arg1: collections.abc.Sequence[int], /) -> dolfinx::fem::FunctionSpace<double>\n\nInvoked with types: dolfinx.cpp.mesh.Mesh_float64, tuple",
-     "output_type": "error",
-     "traceback": [
-      "\u001b[31m---------------------------------------------------------------------------\u001b[39m",
-      "\u001b[31mTypeError\u001b[39m                                 Traceback (most recent call last)",
-      "\u001b[36mCell\u001b[39m\u001b[36m \u001b[39m\u001b[32mIn[9]\u001b[39m\u001b[32m, line 3\u001b[39m\n\u001b[32m      1\u001b[39m V = dolfinx.fem.functionspace(mesh, (\u001b[33m\"\u001b[39m\u001b[33mLagrange\u001b[39m\u001b[33m\"\u001b[39m, \u001b[32m1\u001b[39m))\n\u001b[32m      2\u001b[39m \u001b[38;5;28mprint\u001b[39m(mesh)\n\u001b[32m----> \u001b[39m\u001b[32m3\u001b[39m R = \u001b[43mcreate_real_functionspace\u001b[49m\u001b[43m(\u001b[49m\u001b[43mmesh\u001b[49m\u001b[43m)\u001b[49m\n\u001b[32m      5\u001b[39m W = ufl.MixedFunctionSpace(V, R)\n",
-      "\u001b[36mFile \u001b[39m\u001b[32m~/scifem/src/scifem/spaces.py:32\u001b[39m, in \u001b[36mcreate_real_functionspace\u001b[39m\u001b[34m(mesh, value_shape)\u001b[39m\n\u001b[32m     27\u001b[39m ufl_e = basix.ufl.element(\n\u001b[32m     28\u001b[39m     \u001b[33m\"\u001b[39m\u001b[33mP\u001b[39m\u001b[33m\"\u001b[39m, mesh.basix_cell(), \u001b[32m0\u001b[39m, dtype=dtype, discontinuous=\u001b[38;5;28;01mTrue\u001b[39;00m, shape=value_shape\n\u001b[32m     29\u001b[39m )\n\u001b[32m     31\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m (dtype := mesh.geometry.x.dtype) == np.float64:\n\u001b[32m---> \u001b[39m\u001b[32m32\u001b[39m     cppV = \u001b[43m_scifem\u001b[49m\u001b[43m.\u001b[49m\u001b[43mcreate_real_functionspace_float64\u001b[49m\u001b[43m(\u001b[49m\u001b[43mmesh\u001b[49m\u001b[43m.\u001b[49m\u001b[43m_cpp_object\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mvalue_shape\u001b[49m\u001b[43m)\u001b[49m\n\u001b[32m     33\u001b[39m \u001b[38;5;28;01melif\u001b[39;00m dtype == np.float32:\n\u001b[32m     34\u001b[39m     cppV = _scifem.create_real_functionspace_float32(mesh._cpp_object, value_shape)\n",
-      "\u001b[31mTypeError\u001b[39m: create_real_functionspace_float64(): incompatible function arguments. The following argument types are supported:\n    1. create_real_functionspace_float64(arg0: dolfinx::mesh::Mesh<double>, arg1: collections.abc.Sequence[int], /) -> dolfinx::fem::FunctionSpace<double>\n\nInvoked with types: dolfinx.cpp.mesh.Mesh_float64, tuple"
-     ]
-    }
-   ],
-   "source": [
-    "V = dolfinx.fem.functionspace(mesh, (\"Lagrange\", 1))\n",
-    "R = create_real_functionspace(mesh)\n",
-    "\n",
-    "W = ufl.MixedFunctionSpace(V, R)\n"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "f9d9b61c",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "\n",
-    "\n",
-    "dofs_boundary = dolfinx.fem.locate_dofs_topological(V, fdim, ft.indices[ft.values == wall_marker])\n",
-    "dofs_bubble = dolfinx.fem.locate_dofs_topological(V, fdim, ft.indices[ft.values == obstacle_marker]\n",
-    ")\n",
-    "\n",
-    "\n",
-    "u = dolfinx.fem.Function(V)\n",
-    "u_n = dolfinx.fem.Function(V)\n",
-    "pressure = dolfinx.fem.Function(R)\n",
-    "pressure_n = dolfinx.fem.Function(R)\n",
-    "\n",
-    "dt = dolfinx.fem.Constant(mesh, 1.0)\n",
-    "K_S = dolfinx.fem.Constant(mesh, 1.0)\n",
-    "e = dolfinx.fem.Constant(mesh, 2.0)\n",
-    "volume = dolfinx.fem.Constant(mesh, 4.0)\n",
-    "\n",
-    "v, pressure_v = ufl.TestFunctions(W)\n",
-    "\n",
-    "pressure_ini_expr = dolfinx.fem.Expression(\n",
-    "    dolfinx.fem.Constant(mesh, 3.5), R.element.interpolation_points\n",
-    ")\n",
-    "pressure_n.interpolate(pressure_ini_expr)\n",
-    "pressure.x.array[:] = pressure_n.x.array[:]\n",
-    "\n",
-    "bc_bubble_expr = dolfinx.fem.Expression(K_S * pressure, V.element.interpolation_points)\n",
-    "u_bc_bubble = dolfinx.fem.Function(V)\n",
-    "u_bc_bubble.interpolate(bc_bubble_expr)\n",
-    "\n",
-    "bc_bubble = dolfinx.fem.dirichletbc(u_bc_bubble, dofs_bubble)\n",
-    "bc_boundary = dolfinx.fem.dirichletbc(dolfinx.fem.Constant(mesh, 0.0), dofs_boundary, V)\n",
-    "\n",
-    "n = ufl.FacetNormal(mesh)\n",
-    "flux = ufl.inner(ufl.grad(u), n)\n",
-    "\n",
-    "ds = ufl.Measure(\"ds\", domain=mesh, subdomain_data=ft)\n",
-    "\n",
-    "F = -(u - u_n) / dt * v * ufl.dx + ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx\n",
-    "F += (pressure - pressure_n) / dt * pressure_v * ufl.dx - e / volume * flux * pressure_v * ds(2)\n",
-    "\n",
-    "forms = ufl.extract_blocks(F)\n",
-    "\n",
-    "solver = NonlinearProblem(\n",
-    "    forms,\n",
-    "    [u, pressure],\n",
-    "    bcs=[bc_boundary, bc_bubble],\n",
-    "    petsc_options_prefix=\"bubble\",\n",
-    "    petsc_options={\n",
-    "        \"ksp_type\": \"preonly\",\n",
-    "        \"pc_type\": \"lu\",\n",
-    "        \"pc_factor_mat_solver_type\": \"mumps\",\n",
-    "        \"snes_monitor\": None,\n",
-    "    },\n",
-    ")\n",
-    "\n",
-    "log.set_log_level(log.LogLevel.INFO)\n",
-    "t = 0\n",
-    "t_final = 20\n",
-    "times = []\n",
-    "\n",
-    "all_pressures = []\n",
-    "while t < t_final:\n",
-    "    t += dt.value\n",
-    "    times.append(t)\n",
-    "    print(f\"Solving at time {t:.2f}\")\n",
-    "\n",
-    "    # Solve the problem\n",
-    "    (u, pressure) = solver.solve()\n",
-    "    u_n.x.array[:] = u.x.array[:]\n",
-    "    pressure_n.x.array[:] = pressure.x.array[:]\n",
-    "    print(f\"u = {u.x.array}\")\n",
-    "    print(f\"pressure = {pressure.x.array}\")\n",
-    "\n",
-    "    u_bc_bubble.interpolate(bc_bubble_expr)\n",
-    "    all_pressures.append(pressure.x.array[0].copy())"
-   ]
-  }
- ],
- "metadata": {
-  "jupytext": {
-   "cell_metadata_filter": "-all",
-   "formats": "py,ipynb,md:myst",
-   "main_language": "python"
-  },
-  "kernelspec": {
-   "display_name": "scifem-dev",
-   "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.14.3"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 5
-}
diff --git a/examples/real_function_space_nonlinear.py b/examples/real_function_space_nonlinear.py
new file mode 100644
index 0000000..f047dee
--- /dev/null
+++ b/examples/real_function_space_nonlinear.py
@@ -0,0 +1,319 @@
+# ---
+# jupyter:
+#   jupytext:
+#     cell_metadata_filter: -all
+#     formats: py:percent,ipynb,md:myst
+#     text_representation:
+#       extension: .py
+#       format_name: percent
+#       format_version: '1.3'
+#       jupytext_version: 1.19.1
+#   kernelspec:
+#     display_name: scifem-dev
+#     language: python
+#     name: python3
+# ---
+
+# %% [markdown]
+# # Real function spaces (non linear)
+#
+# Author: Remi Delaporte-Mathurin, Jørgen S. Dokken
+#
+# License: MIT
+#
+# In this example we will show how to use the "real" function space to solve
+# a non linear problem.
+#
+# ## Mathematical formulation
+#  The problem at hand is:
+# Find $u \in H^1(\Omega)$ such that
+#
+# $$
+# \begin{align}
+# \Delta u &= 0 \quad \text{in } \Omega, \\
+# u &= 0 \quad \text{on } \Gamma_1, \\
+# u &= K P_b \quad \text{on } \Gamma_b, \\
+# \end{align}
+# $$
+#
+# where $P_b$ is the partial pressure inside the cavity region. The temporal evolution of $P_b$ is governed by:
+#
+# $$
+# \frac{dP_b}{dt} = \frac{e}{V} \int_{\Gamma_b} -D \nabla c \cdot \mathbf{n} dS
+# $$
+#
+#
+# We write the weak formulation:
+#
+# $$
+# \begin{align}
+# \int_\Omega  \frac{u - u_n}{\Delta t} \cdot v~\mathrm{d}x + \int_\Omega \nabla u \cdot \nabla v~\mathrm{d}x &= 0\\
+# \int_\Omega  \frac{P_b - P_{b_n}}{\Delta t} \cdot w~\mathrm{d}x &= \frac{e}{V} \int_{\Gamma_b} -D \nabla c \cdot \mathbf{n} ~dS ~w.
+# \end{align}
+# $$
+#
+
+# %% [markdown]
+# ## Implementation
+# We start by import the necessary modules
+# ```{admonition} Clickable functions/classes
+# Note that for the modules imported in this example, you can click on the function/class
+# name to be redirected to the corresponding documentation page.
+# ```
+
+# %%
+from mpi4py import MPI
+import dolfinx
+from dolfinx import log
+from dolfinx.fem.petsc import NonlinearProblem
+from dolfinx.io import VTXWriter, gmsh as gmshio
+import numpy as np
+from scifem import create_real_functionspace
+import ufl
+
+import gmsh
+
+# %%
+gmsh.initialize()
+
+L = 0.5
+H = 0.41
+c_x = c_y = 0.2
+r = 0.05
+gdim = 2
+mesh_comm = MPI.COMM_WORLD
+model_rank = 0
+if mesh_comm.rank == model_rank:
+    rectangle = gmsh.model.occ.addRectangle(0, 0, 0, L, H, tag=1)
+    obstacle = gmsh.model.occ.addDisk(c_x, c_y, 0, r, r)
+
+if mesh_comm.rank == model_rank:
+    fluid = gmsh.model.occ.cut([(gdim, rectangle)], [(gdim, obstacle)])
+    gmsh.model.occ.synchronize()
+
+solid_marker = 1
+if mesh_comm.rank == model_rank:
+    volumes = gmsh.model.getEntities(dim=gdim)
+    assert len(volumes) == 1
+    gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], solid_marker)
+    gmsh.model.setPhysicalName(volumes[0][0], solid_marker, "Solid")
+
+wall_marker, obstacle_marker = 2, 3
+walls, obstacle = [], []
+if mesh_comm.rank == model_rank:
+    boundaries = gmsh.model.getBoundary(volumes, oriented=False)
+    for boundary in boundaries:
+        center_of_mass = gmsh.model.occ.getCenterOfMass(boundary[0], boundary[1])
+        if np.allclose(center_of_mass, [0, H / 2, 0]):
+            walls.append(boundary[1])
+        elif np.allclose(center_of_mass, [L, H / 2, 0]):
+            walls.append(boundary[1])
+        elif np.allclose(center_of_mass, [L / 2, H, 0]) or np.allclose(
+            center_of_mass, [L / 2, 0, 0]
+        ):
+            walls.append(boundary[1])
+        else:
+            obstacle.append(boundary[1])
+    gmsh.model.addPhysicalGroup(1, walls, wall_marker)
+    gmsh.model.setPhysicalName(1, wall_marker, "Walls")
+    gmsh.model.addPhysicalGroup(1, obstacle, obstacle_marker)
+    gmsh.model.setPhysicalName(1, obstacle_marker, "Obstacle")
+
+
+res_min = r / 3
+if mesh_comm.rank == model_rank:
+    distance_field = gmsh.model.mesh.field.add("Distance")
+    gmsh.model.mesh.field.setNumbers(distance_field, "EdgesList", obstacle)
+    threshold_field = gmsh.model.mesh.field.add("Threshold")
+    gmsh.model.mesh.field.setNumber(threshold_field, "IField", distance_field)
+    gmsh.model.mesh.field.setNumber(threshold_field, "LcMin", res_min)
+    gmsh.model.mesh.field.setNumber(threshold_field, "LcMax", 0.25 * H)
+    gmsh.model.mesh.field.setNumber(threshold_field, "DistMin", r)
+    gmsh.model.mesh.field.setNumber(threshold_field, "DistMax", 2 * H)
+    min_field = gmsh.model.mesh.field.add("Min")
+    gmsh.model.mesh.field.setNumbers(min_field, "FieldsList", [threshold_field])
+    gmsh.model.mesh.field.setAsBackgroundMesh(min_field)
+
+
+if mesh_comm.rank == model_rank:
+    gmsh.option.setNumber("Mesh.Algorithm", 8)
+    gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 2)
+    gmsh.option.setNumber("Mesh.RecombineAll", 1)
+    gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1)
+    gmsh.model.mesh.generate(gdim)
+    gmsh.model.mesh.setOrder(2)
+    gmsh.model.mesh.optimize("Netgen")
+
+mesh_data = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)
+mesh = mesh_data.mesh
+assert mesh_data.facet_tags is not None
+ft = mesh_data.facet_tags
+ft.name = "Facet markers"
+
+# %%
+from dolfinx import plot
+import pyvista
+
+pyvista.set_jupyter_backend("html")
+
+tdim = mesh.topology.dim
+
+mesh.topology.create_connectivity(tdim, tdim)
+topology, cell_types, geometry = plot.vtk_mesh(mesh, tdim)
+grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
+
+plotter = pyvista.Plotter()
+plotter.add_mesh(grid, show_edges=True)
+plotter.view_xy()
+if not pyvista.OFF_SCREEN:
+    plotter.show()
+else:
+    figure = plotter.screenshot("mesh.png")
+
+fdim = tdim - 1
+mesh.topology.create_connectivity(fdim, tdim)
+topology, cell_types, x = plot.vtk_mesh(mesh, fdim, ft.indices)
+
+p = pyvista.Plotter()
+grid = pyvista.UnstructuredGrid(topology, cell_types, x)
+grid.cell_data["Facet Marker"] = ft.values
+grid.set_active_scalars("Facet Marker")
+p.view_xy()
+p.add_mesh(grid, show_edges=True)
+
+if not pyvista.OFF_SCREEN:
+    p.show()
+else:
+    figure = p.screenshot("facet_markers.png")
+
+# %%
+V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
+R = create_real_functionspace(mesh)
+
+W = ufl.MixedFunctionSpace(V, R)
+
+# %%
+dofs_boundary = dolfinx.fem.locate_dofs_topological(V, fdim, ft.indices[ft.values == wall_marker])
+dofs_bubble = dolfinx.fem.locate_dofs_topological(V, fdim, ft.indices[ft.values == obstacle_marker]
+)
+
+
+u = dolfinx.fem.Function(V)
+u_n = dolfinx.fem.Function(V)
+pressure = dolfinx.fem.Function(R)
+pressure_n = dolfinx.fem.Function(R)
+
+dt = dolfinx.fem.Constant(mesh, 0.05)
+K_S = dolfinx.fem.Constant(mesh, 2.0)
+e = dolfinx.fem.Constant(mesh, 2.0)
+volume = dolfinx.fem.Constant(mesh, 40.0)
+
+v, pressure_v = ufl.TestFunctions(W)
+
+pressure_ini_expr = dolfinx.fem.Expression(
+    dolfinx.fem.Constant(mesh, 3.5), R.element.interpolation_points
+)
+pressure_n.interpolate(pressure_ini_expr)
+pressure.x.array[:] = pressure_n.x.array[:]
+
+bc_bubble_expr = dolfinx.fem.Expression(K_S * pressure, V.element.interpolation_points)
+u_bc_bubble = dolfinx.fem.Function(V)
+u_bc_bubble.interpolate(bc_bubble_expr)
+
+bc_bubble = dolfinx.fem.dirichletbc(u_bc_bubble, dofs_bubble)
+bc_boundary = dolfinx.fem.dirichletbc(dolfinx.fem.Constant(mesh, 0.0), dofs_boundary, V)
+
+n = ufl.FacetNormal(mesh)
+flux = ufl.inner(ufl.grad(u), n)
+
+ds = ufl.Measure("ds", domain=mesh, subdomain_data=ft)
+
+F = -(u - u_n) / dt * v * ufl.dx + ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx
+F += (pressure - pressure_n) / dt * pressure_v * ufl.dx - e / volume * flux * pressure_v * ds(2)
+
+forms = ufl.extract_blocks(F)
+
+solver = NonlinearProblem(
+    forms,
+    [u, pressure],
+    bcs=[bc_boundary, bc_bubble],
+    petsc_options_prefix="bubble",
+    petsc_options={
+        "ksp_type": "preonly",
+        "pc_type": "lu",
+        "pc_factor_mat_solver_type": "mumps",
+        "snes_monitor": None,
+    },
+)
+
+
+# %%
+import matplotlib.pyplot as plt
+import matplotlib as mpl
+grid = pyvista.UnstructuredGrid(*plot.vtk_mesh(V))
+
+plotter = pyvista.Plotter()
+plotter.open_gif("u_time.gif", fps=10)
+
+grid.point_data["u"] = u.x.array
+warped = grid.warp_by_scalar("u", factor=0.1)
+
+viridis = mpl.colormaps.get_cmap("viridis").resampled(25)
+sargs = dict(
+    title_font_size=25,
+    label_font_size=20,
+    fmt="%.2e",
+    color="black",
+    position_x=0.1,
+    position_y=0.8,
+    width=0.8,
+    height=0.1,
+)
+
+renderer = plotter.add_mesh(
+    warped,
+    show_edges=True,
+    lighting=False,
+    cmap=viridis,
+    scalar_bar_args=sargs,
+    clim=[0, max(u_bc_bubble.x.array)],
+)
+
+# %%
+
+log.set_log_level(log.LogLevel.INFO)
+t = 0
+t_final = 4
+times = []
+
+all_pressures = []
+while t < t_final:
+    t += dt.value
+    times.append(t)
+    print(f"Solving at time {t:.2f}")
+
+    # Solve the problem
+    (u, pressure) = solver.solve()
+    u_n.x.array[:] = u.x.array[:]
+    pressure_n.x.array[:] = pressure.x.array[:]
+    print(f"u = {u.x.array}")
+    print(f"pressure = {pressure.x.array}")
+
+    u_bc_bubble.interpolate(bc_bubble_expr)
+    all_pressures.append(pressure.x.array[0].copy())
+
+    # Update plot
+    new_warped = grid.warp_by_scalar("u", factor=0.1)
+    warped.points[:, :] = new_warped.points
+    warped.point_data["u"][:] = u.x.array
+    plotter.write_frame()
+plotter.close()
+
+# %%
+import matplotlib.pyplot as plt
+
+plt.scatter([0], [3.5], color="red", label="Initial Pressure")
+plt.plot(times, all_pressures)
+plt.xlabel("Time")
+plt.ylabel("Pressure")
+plt.show()

From 031a2d8ad43125b6b92820de151f1d90e06ece56 Mon Sep 17 00:00:00 2001
From: RemDelaporteMathurin <rdelaportemathurin@gmail.com>
Date: Tue, 24 Feb 2026 11:31:50 -0500
Subject: [PATCH 3/5] square + break code blocks

---
 examples/real_function_space_nonlinear.py | 30 ++++++++++++++---------
 1 file changed, 18 insertions(+), 12 deletions(-)

diff --git a/examples/real_function_space_nonlinear.py b/examples/real_function_space_nonlinear.py
index f047dee..0f538cc 100644
--- a/examples/real_function_space_nonlinear.py
+++ b/examples/real_function_space_nonlinear.py
@@ -66,7 +66,7 @@
 import dolfinx
 from dolfinx import log
 from dolfinx.fem.petsc import NonlinearProblem
-from dolfinx.io import VTXWriter, gmsh as gmshio
+from dolfinx.io import gmsh as gmshio
 import numpy as np
 from scifem import create_real_functionspace
 import ufl
@@ -77,8 +77,8 @@
 gmsh.initialize()
 
 L = 0.5
-H = 0.41
-c_x = c_y = 0.2
+H = L
+c_x = c_y = L/2
 r = 0.05
 gdim = 2
 mesh_comm = MPI.COMM_WORLD
@@ -193,23 +193,21 @@
 W = ufl.MixedFunctionSpace(V, R)
 
 # %%
-dofs_boundary = dolfinx.fem.locate_dofs_topological(V, fdim, ft.indices[ft.values == wall_marker])
-dofs_bubble = dolfinx.fem.locate_dofs_topological(V, fdim, ft.indices[ft.values == obstacle_marker]
-)
-
-
 u = dolfinx.fem.Function(V)
 u_n = dolfinx.fem.Function(V)
 pressure = dolfinx.fem.Function(R)
 pressure_n = dolfinx.fem.Function(R)
 
+v, pressure_v = ufl.TestFunctions(W)
+
+
+# %%
 dt = dolfinx.fem.Constant(mesh, 0.05)
 K_S = dolfinx.fem.Constant(mesh, 2.0)
 e = dolfinx.fem.Constant(mesh, 2.0)
 volume = dolfinx.fem.Constant(mesh, 40.0)
 
-v, pressure_v = ufl.TestFunctions(W)
-
+# %%
 pressure_ini_expr = dolfinx.fem.Expression(
     dolfinx.fem.Constant(mesh, 3.5), R.element.interpolation_points
 )
@@ -220,9 +218,14 @@
 u_bc_bubble = dolfinx.fem.Function(V)
 u_bc_bubble.interpolate(bc_bubble_expr)
 
+dofs_boundary = dolfinx.fem.locate_dofs_topological(V, fdim, ft.indices[ft.values == wall_marker])
+dofs_bubble = dolfinx.fem.locate_dofs_topological(V, fdim, ft.indices[ft.values == obstacle_marker]
+)
+
 bc_bubble = dolfinx.fem.dirichletbc(u_bc_bubble, dofs_bubble)
 bc_boundary = dolfinx.fem.dirichletbc(dolfinx.fem.Constant(mesh, 0.0), dofs_boundary, V)
 
+# %%
 n = ufl.FacetNormal(mesh)
 flux = ufl.inner(ufl.grad(u), n)
 
@@ -233,6 +236,8 @@
 
 forms = ufl.extract_blocks(F)
 
+# %%
+
 solver = NonlinearProblem(
     forms,
     [u, pressure],
@@ -246,9 +251,7 @@
     },
 )
 
-
 # %%
-import matplotlib.pyplot as plt
 import matplotlib as mpl
 grid = pyvista.UnstructuredGrid(*plot.vtk_mesh(V))
 
@@ -309,6 +312,9 @@
     plotter.write_frame()
 plotter.close()
 
+# %% [markdown]
+# ![title](u_time.gif)
+
 # %%
 import matplotlib.pyplot as plt
 

From 49340e3643b8506390f64c358c58964cce3af7ac Mon Sep 17 00:00:00 2001
From: RemDelaporteMathurin <rdelaportemathurin@gmail.com>
Date: Tue, 24 Feb 2026 12:28:22 -0500
Subject: [PATCH 4/5] added some write up

---
 examples/real_function_space_nonlinear.py | 93 ++++++++++++++++++-----
 1 file changed, 74 insertions(+), 19 deletions(-)

diff --git a/examples/real_function_space_nonlinear.py b/examples/real_function_space_nonlinear.py
index 0f538cc..7e44938 100644
--- a/examples/real_function_space_nonlinear.py
+++ b/examples/real_function_space_nonlinear.py
@@ -24,6 +24,11 @@
 # In this example we will show how to use the "real" function space to solve
 # a non linear problem.
 #
+# A circular gas cavity inside a solid domain has an initial partial pressure $P_b$. The concentration $u$ on the cavity surface follows Henry's solubility law: the concentration is proportional to the partial pressure $P_b$.
+#
+# As particles leave the cavity by solution/diffusion, the partial pressure starts to decrease - affecting in turn the concentration on the cavity surface.
+#
+#
 # ## Mathematical formulation
 #  The problem at hand is:
 # Find $u \in H^1(\Omega)$ such that
@@ -43,7 +48,7 @@
 # $$
 #
 #
-# We write the weak formulation:
+# We write the weak formulation with $v$ and $w$ suitable test functions:
 #
 # $$
 # \begin{align}
@@ -64,19 +69,21 @@
 # %%
 from mpi4py import MPI
 import dolfinx
-from dolfinx import log
 from dolfinx.fem.petsc import NonlinearProblem
 from dolfinx.io import gmsh as gmshio
 import numpy as np
-from scifem import create_real_functionspace
+from scifem import create_real_functionspace, assemble_scalar
 import ufl
 
 import gmsh
 
+# %% [markdown]
+# We generate the mesh using GMSH:
+
 # %%
 gmsh.initialize()
 
-L = 0.5
+L = 2
 H = L
 c_x = c_y = L/2
 r = 0.05
@@ -186,12 +193,18 @@
 else:
     figure = p.screenshot("facet_markers.png")
 
+# %% [markdown]
+# We create two functionspaces `V` and `R` for $u$ and $P_b$, respectively, as well as a `ufl.MixedFunctionSpace` for test functions.
+
 # %%
 V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
 R = create_real_functionspace(mesh)
 
 W = ufl.MixedFunctionSpace(V, R)
 
+# %% [markdown]
+# We then create appropriate functions and test functions:
+
 # %%
 u = dolfinx.fem.Function(V)
 u_n = dolfinx.fem.Function(V)
@@ -201,15 +214,28 @@
 v, pressure_v = ufl.TestFunctions(W)
 
 
+# %% [markdown]
+# Let's define the problems constants:
+
 # %%
-dt = dolfinx.fem.Constant(mesh, 0.05)
+dt = dolfinx.fem.Constant(mesh, 0.1)
 K_S = dolfinx.fem.Constant(mesh, 2.0)
 e = dolfinx.fem.Constant(mesh, 2.0)
 volume = dolfinx.fem.Constant(mesh, 40.0)
+P_b_initial = dolfinx.fem.Constant(mesh, 3.5)
+u_out = dolfinx.fem.Constant(mesh, 0.0)
+
+# %% [markdown]
+# ```{note}
+# If `dt` is too large, the problem becomes unstable
+# ```
+#
+#
+# We can now define the initial condition and boundary conditions:
 
 # %%
 pressure_ini_expr = dolfinx.fem.Expression(
-    dolfinx.fem.Constant(mesh, 3.5), R.element.interpolation_points
+    P_b_initial, R.element.interpolation_points
 )
 pressure_n.interpolate(pressure_ini_expr)
 pressure.x.array[:] = pressure_n.x.array[:]
@@ -223,7 +249,10 @@
 )
 
 bc_bubble = dolfinx.fem.dirichletbc(u_bc_bubble, dofs_bubble)
-bc_boundary = dolfinx.fem.dirichletbc(dolfinx.fem.Constant(mesh, 0.0), dofs_boundary, V)
+bc_boundary = dolfinx.fem.dirichletbc(u_out, dofs_boundary, V)
+
+# %% [markdown]
+# Next we define the variational formulation and call `ufl.extract_blocks` to form the blocked formulations:
 
 # %%
 n = ufl.FacetNormal(mesh)
@@ -231,13 +260,15 @@
 
 ds = ufl.Measure("ds", domain=mesh, subdomain_data=ft)
 
-F = -(u - u_n) / dt * v * ufl.dx + ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx
-F += (pressure - pressure_n) / dt * pressure_v * ufl.dx - e / volume * flux * pressure_v * ds(2)
+F = (u - u_n) / dt * v * ufl.dx + ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx
+F += (pressure - pressure_n) / dt * pressure_v * ufl.dx + e / volume * flux * pressure_v * ds(3)
 
 forms = ufl.extract_blocks(F)
 
-# %%
+# %% [markdown]
+# We can now create a nonlinear solver with the blocked formulations and the functions `u` and `pressure` as a list:
 
+# %%
 solver = NonlinearProblem(
     forms,
     [u, pressure],
@@ -251,6 +282,9 @@
     },
 )
 
+# %% [markdown]
+# Set up transient pyvista visualisation:
+
 # %%
 import matplotlib as mpl
 grid = pyvista.UnstructuredGrid(*plot.vtk_mesh(V))
@@ -282,14 +316,18 @@
     clim=[0, max(u_bc_bubble.x.array)],
 )
 
-# %%
+# %% [markdown]
+# Time stepping loop:
 
-log.set_log_level(log.LogLevel.INFO)
+# %%
 t = 0
-t_final = 4
-times = []
+t_final = 15
 
+times = []
 all_pressures = []
+outgassing_fluxes = []
+
+
 while t < t_final:
     t += dt.value
     times.append(t)
@@ -297,11 +335,12 @@
 
     # Solve the problem
     (u, pressure) = solver.solve()
+
+    # Update previous solution
     u_n.x.array[:] = u.x.array[:]
     pressure_n.x.array[:] = pressure.x.array[:]
-    print(f"u = {u.x.array}")
-    print(f"pressure = {pressure.x.array}")
 
+    # Update bubble BC
     u_bc_bubble.interpolate(bc_bubble_expr)
     all_pressures.append(pressure.x.array[0].copy())
 
@@ -310,16 +349,32 @@
     warped.points[:, :] = new_warped.points
     warped.point_data["u"][:] = u.x.array
     plotter.write_frame()
+
+    # compute outgassing flux
+    outgassing_flux = -assemble_scalar(flux * ds(2))
+    outgassing_fluxes.append(outgassing_flux)
 plotter.close()
 
 # %% [markdown]
 # ![title](u_time.gif)
 
+# %% [markdown]
+# We see that, as expected, the partial pressure $P_b$ decreases with time.
+
 # %%
 import matplotlib.pyplot as plt
 
-plt.scatter([0], [3.5], color="red", label="Initial Pressure")
-plt.plot(times, all_pressures)
+fig, axs = plt.subplots(2, 1, figsize=(10, 8), sharex=True)
+
+
+axs[0].scatter([0], [P_b_initial.value], color="red", label="Initial Pressure")
+axs[0].plot(times, all_pressures)
+axs[0].set_ylabel("Pressure")
+axs[0].set_ylim(bottom=0)
+
+axs[1].plot(times, outgassing_fluxes)
+axs[1].set_ylabel("Outgassing Flux")
+axs[1].set_ylim(bottom=0)
+
 plt.xlabel("Time")
-plt.ylabel("Pressure")
 plt.show()

From 96c1eb751dfc3817c3fde5e5c56242dad1a962c0 Mon Sep 17 00:00:00 2001
From: RemDelaporteMathurin <rdelaportemathurin@gmail.com>
Date: Fri, 27 Feb 2026 08:53:02 -0500
Subject: [PATCH 5/5] c/u inconsistency

---
 examples/real_function_space_nonlinear.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/examples/real_function_space_nonlinear.py b/examples/real_function_space_nonlinear.py
index 7e44938..76c5fae 100644
--- a/examples/real_function_space_nonlinear.py
+++ b/examples/real_function_space_nonlinear.py
@@ -2,7 +2,7 @@
 # jupyter:
 #   jupytext:
 #     cell_metadata_filter: -all
-#     formats: py:percent,ipynb,md:myst
+#     formats: ipynb,md:myst,py:percent
 #     text_representation:
 #       extension: .py
 #       format_name: percent
@@ -44,7 +44,7 @@
 # where $P_b$ is the partial pressure inside the cavity region. The temporal evolution of $P_b$ is governed by:
 #
 # $$
-# \frac{dP_b}{dt} = \frac{e}{V} \int_{\Gamma_b} -D \nabla c \cdot \mathbf{n} dS
+# \frac{dP_b}{dt} = \frac{e}{V} \int_{\Gamma_b} -D \nabla u \cdot \mathbf{n} dS
 # $$
 #
 #
@@ -53,7 +53,7 @@
 # $$
 # \begin{align}
 # \int_\Omega  \frac{u - u_n}{\Delta t} \cdot v~\mathrm{d}x + \int_\Omega \nabla u \cdot \nabla v~\mathrm{d}x &= 0\\
-# \int_\Omega  \frac{P_b - P_{b_n}}{\Delta t} \cdot w~\mathrm{d}x &= \frac{e}{V} \int_{\Gamma_b} -D \nabla c \cdot \mathbf{n} ~dS ~w.
+# \int_\Omega  \frac{P_b - P_{b_n}}{\Delta t} \cdot w~\mathrm{d}x &= \frac{e}{V} \int_{\Gamma_b} -D \nabla u \cdot \mathbf{n} ~dS \cdot ~w.
 # \end{align}
 # $$
 #