heat_eq.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "68d7991c-6664-49d6-a8a2-71b43fc0abda",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "# Solving a time-dependent problem"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f44f5cf7-1344-474e-b121-936caa13797a",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    },
    "tags": []
   },
   "source": [
    "This notebook will show you how to solve a transient problem using DOLFINx, and highlight differences between legacy DOLFIN and DOLFINx.\n",
    "We start by looking at the structure of DOLFINx:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36ccb951-19ab-4b69-a13c-e0a3356f305b",
   "metadata": {},
   "source": [
    "Relevant DOLFINx modules:\n",
    "- `dolfinx.mesh`: Classes and functions related to the computational domain\n",
    "- `dolfinx.fem`: Finite element method functionality\n",
    "- `dolfinx.io`: Input/Output (read/write) functionality\n",
    "- `dolfinx.plot`: Convenience functions for exporting plotting data\n",
    "- `dolfinx.la`: Functions related to linear algebra structures (matrices/vectors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "a9fbb5d1-0dad-44f6-9a32-93535f74db4e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from dolfinx import fem, la, mesh, plot"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28904567-9391-44c1-9e7b-b55e7486a93b",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Creating a distributed computational domain (mesh)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "54fb77a5-7b70-4e7f-9cf3-e1c88251b796",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    },
    "tags": []
   },
   "source": [
    "To create a simple computational domain in DOLFINx, we use the mesh generation utilities in `dolfinx.mesh`. In this module, we have the tools to build rectangles of triangular or quadrilateral elements and boxes of tetrahedral or hexahedral elements. We start by creating a rectangle spanning $[0,0]\\times[10,3]$, with 100 and 20 elements in each direction respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "e473b08c-ef9f-4756-a981-47b2958c7c4e",
   "metadata": {},
   "outputs": [],
   "source": [
    "from mpi4py import MPI\n",
    "\n",
    "length, height = 10, 3\n",
    "Nx, Ny = 80, 60\n",
    "extent = [[0.0, 0.0], [length, height]]\n",
    "domain = mesh.create_rectangle(MPI.COMM_WORLD, extent, [Nx, Ny], mesh.CellType.quadrilateral)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6282cedf-36df-4209-9a54-5f274f74c9ab",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    },
    "tags": []
   },
   "source": [
    "In constrast to legacy DOLFIN, we work on simple Python structures (nested listes, numpy arrays, etc).\n",
    "We also note that we have to send in an MPI communicator. This is because we want the user to be aware of how the mesh is distributed when running in parallel. \n",
    "If we use the communicator `MPI.COMM_SELF`, each process initialised when running the script would have a version of the full mesh local to its process."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c7e163a-acdd-49f3-8873-fde09f59b8d0",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Creating a mesh on each process"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "23061a93-3b46-45ab-ba95-987e01db1cca",
   "metadata": {
    "lines_to_next_cell": 0,
    "tags": []
   },
   "outputs": [],
   "source": [
    "local_domain = mesh.create_rectangle(MPI.COMM_SELF, extent, [Nx, Ny], mesh.CellType.quadrilateral)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e96dbc06-c46d-4455-8d85-a5eb4764e532",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": []
   },
   "source": [
    "We plot the mesh."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "2233da8a-ff57-4771-88af-fa8fd51c1ba7",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "import pyvista\n",
    "\n",
    "pyvista.start_xvfb()\n",
    "grid = pyvista.UnstructuredGrid(*plot.vtk_mesh(local_domain))\n",
    "plotter = pyvista.Plotter(window_size=(800, 400))\n",
    "renderer = plotter.add_mesh(grid, show_edges=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "56b553a6-8bad-4d22-b9fb-dc6657b5d99b",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    },
    "tags": []
   },
   "source": [
    "With Pyvista, we can export the plots to many formats including pngs, interactive notebook plots, and html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "49d7c755-d8ae-4c9a-8259-9c5fe8ef47ea",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    },
    "tags": [
     "hide-output"
    ]
   },
   "outputs": [],
   "source": [
    "plotter.view_xy()\n",
    "plotter.camera.zoom(2)\n",
    "plotter.export_html(\"./beam.html\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "4eb5a932-4d6c-46be-96a6-2a515f693cfd",
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<iframe src='./beam.html' scrolling=\"no\" width=\"800px\" height=\"400px\"></iframe> <!--  # noqa, -->\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%html\n",
    "<iframe src='./beam.html' scrolling=\"no\" width=\"800px\" height=\"400px\"></iframe> <!--  # noqa, -->"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f792ea5-ce66-4f91-a2d7-10f5a361b9d6",
   "metadata": {
    "lines_to_next_cell": 2,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Setting up a variational problem"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eea38f11-101f-4fd8-a8bb-af43979cc839",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    },
    "tags": []
   },
   "source": [
    "We will solve the heat equation, with a backward Euler time stepping scheme, ie"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3578231c-f2ff-4cbe-beba-aeef1e7a3c40",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "\\frac{u_{n+1}-u_n}{\\Delta t} - \\nabla \\cdot (\\mu  \\nabla u_{n+1}) &= f(x,t_{n+1}) && \\text{in } \\Omega,\\\\\n",
    "u &= u_D(x,t_{n+1}) &&\\text{on } \\partial\\Omega_\\text{D},\\\\\n",
    "\\mu\\frac{\\partial u_{n+1}}{\\partial n} &=0 &&\\text{on } \\partial\\Omega_\\text{N},\n",
    "\\end{align*}\n",
    "$$ \n",
    "with $u_D = y\\cos(0.25t)$, $f=0$. For this example, we take $\\Omega$ to be rectangle defined above, $\\Omega_\\text{D}$ if the left-hand edge of the rectangle, and $\\Omega_\\text{N}$ is the remaining three edges of the rectangle."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7fdc0d5f-e35b-4146-9839-ddf3a791f3c0",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    },
    "tags": []
   },
   "source": [
    "We start by defining the function space, the corresponding test and trial functions, as well as material and temporal parameters. Note that we use explicit imports from UFL to create the test and trial functions, to avoid confusion as to where they originate from. DOLFINx and UFL support both real and complex valued functions. However, to be able to use the PETSc linear algebra backend, which only supports a single floating type at compilation, we need to use appropriate scalar types in our variational form. This ensures that we generate consistent matrices and vectors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "1d859674-2bb7-4454-a538-f935b0f56558",
   "metadata": {},
   "outputs": [],
   "source": [
    "from ufl import (\n",
    "    SpatialCoordinate,\n",
    "    TestFunction,\n",
    "    TrialFunction,\n",
    "    as_vector,\n",
    "    dx,\n",
    "    grad,\n",
    "    inner,\n",
    "    system,\n",
    ")\n",
    "\n",
    "V = fem.functionspace(domain, (\"Lagrange\", 1))\n",
    "u = TrialFunction(V)\n",
    "v = TestFunction(V)\n",
    "un = fem.Function(V)\n",
    "f = fem.Constant(domain, 0.0)\n",
    "mu = fem.Constant(domain, 2.3)\n",
    "dt = fem.Constant(domain, 0.05)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83e93c23-e8ea-47c0-8912-38be470bc43a",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    },
    "tags": []
   },
   "source": [
    "The variational form can be written in UFL syntax, as done in legacy DOLFIN:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "b4e8878f-5f35-427c-9318-71d4bcdfdcce",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "F = inner(u - un, v) * dx + dt * mu * inner(grad(u), grad(v)) * dx\n",
    "F -= dt * inner(f, v) * dx\n",
    "(a, L) = system(F)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a7b2f012-1b98-466b-bc6a-5634f91a431a",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Creating Dirichlet boundary conditions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "736df497-ad53-4d45-ac67-09c4b5c97d58",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    },
    "tags": []
   },
   "source": [
    "### Creating a time dependent boundary condition"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a1e26c6f-40df-40f4-adef-cf734511ae2a",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    },
    "tags": []
   },
   "source": [
    "There are many ways of creating boundary conditions. In this example, we will create function $u_\\text{D}(x,t)$ dependent on both space and time. To do this, we define a function that takes a 2-dimensional array `x`.  Each column of `x` corresponds to an input coordinate $(x,y,z)$ and this function operates directly on the columns of `x`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "5060c070-1cd2-43b5-b162-d18ca7490d96",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "\n",
    "def uD_function(t):\n",
    "    return lambda x: x[1] * np.cos(0.25 * t)\n",
    "\n",
    "\n",
    "uD = fem.Function(V)\n",
    "t = 0\n",
    "uD.interpolate(uD_function(t))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd839336-90e1-4383-9b1f-1fcd2e6e82de",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    },
    "tags": []
   },
   "source": [
    "To give the user freedom to set boundary conditions on single degrees of freedom, the function `dolfinx.fem.dirichletbc` takes in the list of degrees of freedom (DOFs) as input. The DOFs on the boundary can be obtained in many ways: DOLFINx supplies a few convenience functions, such as `dolfinx.fem.locate_dofs_topological` and `dolfinx.fem.locate_dofs_geometrical`.\n",
    "Locating dofs topologically is generally advised, as certain finite elements have DOFs that do not have a geometrical coordinates associated with them (eg Nédélec and Raviart--Thomas). DOLFINx also has convenicence functions to obtain a list of all boundary facets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "e3c7dd5e-af0b-4223-a0aa-1ca3c217b93a",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def dirichlet_facets(x):\n",
    "    return np.isclose(x[0], length)\n",
    "\n",
    "\n",
    "tdim = domain.topology.dim\n",
    "bc_facets = mesh.locate_entities_boundary(domain, tdim - 1, dirichlet_facets)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "bb57b84a-e230-45f4-9a11-466fca63f61b",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "bndry_dofs = fem.locate_dofs_topological(V, tdim - 1, bc_facets)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "e207eeca-5e43-47b0-8993-460a7341eb56",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "bcs = [fem.dirichletbc(uD, bndry_dofs)]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd747071-0c80-4c48-b018-58769e121fd9",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Setting up a time dependent solver"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb7586ac-6788-41c8-a828-d1f4618c98a2",
   "metadata": {
    "lines_to_next_cell": 0,
    "slideshow": {
     "slide_type": "notes"
    },
    "tags": []
   },
   "source": [
    "As the left hand side of our problem (the matrix) is time independent, we would like avoid re-assembling it at every time step. DOLFINx gives the user more control over assembly so that this can be done. We assemble the matrix once outside the temporal loop."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "78c9e008-6532-42b5-9369-e58cb422be2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import dolfinx.fem.petsc as petsc\n",
    "\n",
    "compiled_a = fem.form(a)\n",
    "A = petsc.assemble_matrix(compiled_a, bcs=bcs)\n",
    "A.assemble()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d4515c5-0284-4293-b527-7bdbd79bd31a",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    },
    "tags": []
   },
   "source": [
    "Next, we can generate the integration kernel for the right hand side (RHS), and create the RHS vector `b` that we will assemble into at each time step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "4fc3678a-48f6-42b7-b7f5-f1bc33566fa5",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "compiled_L = fem.form(L)\n",
    "b = fem.Function(V)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ae18fe0-81a1-42f8-ba3a-38797cfc41ac",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    },
    "tags": []
   },
   "source": [
    "We next create the PETSc KSP (Krylov subspace method) solver, and set it to solve using an [algebraic multigrid method](https://hypre.readthedocs.io/en/latest/solvers-boomeramg.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "624d6a1f-7854-4a15-bfd9-fb1b4fdddccb",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from petsc4py import PETSc\n",
    "\n",
    "solver = PETSc.KSP().create(domain.comm)\n",
    "solver.setOperators(A)\n",
    "solver.setType(PETSc.KSP.Type.CG)\n",
    "pc = solver.getPC()\n",
    "pc.setType(PETSc.PC.Type.HYPRE)\n",
    "pc.setHYPREType(\"boomeramg\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7fd83dec-1a2c-4364-b013-cee06da3577d",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Plotting a time dependent problem\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "94070aa8-0e5f-4306-a51e-feb8ad429cb0",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    },
    "tags": []
   },
   "source": [
    "As we are solving a time dependent problem, we would like to create a time dependent animation of the solution. \n",
    "We do this by using [pyvista](https://docs.pyvista.org/), which uses VTK structures for plotting.\n",
    "In DOLFINx, we have the convenience function `dolfinx.plot.create_vtk_mesh` that can create meshes compatible with VTK formatting, based on meshes of (discontinuous) Lagrange function spaces."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "e74f1c58-aa5f-4609-b467-89014857d217",
   "metadata": {
    "tags": [
     "hide-cell"
    ]
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import pyvista\n",
    "\n",
    "pyvista.start_xvfb(0.5)  # Start virtual framebuffer for plotting\n",
    "plotter = pyvista.Plotter()\n",
    "plotter.open_gif(\"u_time.gif\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "ab139ce4-5204-4756-83f3-df25ec66579d",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "topology, cells, geometry = plot.vtk_mesh(V)\n",
    "uh = fem.Function(V)\n",
    "grid = pyvista.UnstructuredGrid(topology, cells, geometry)\n",
    "grid.point_data[\"uh\"] = uh.x.array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "5791a8ab-b377-468a-9fa8-47a093658873",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "viridis = plt.cm.get_cmap(\"viridis\", 25)\n",
    "sargs = dict(\n",
    "    title_font_size=25,\n",
    "    label_font_size=20,\n",
    "    fmt=\"%.2e\",\n",
    "    color=\"black\",\n",
    "    position_x=0.1,\n",
    "    position_y=0.8,\n",
    "    width=0.8,\n",
    "    height=0.1,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "8217df02-6cb5-433a-9bc2-7c1903e69e6d",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "renderer = plotter.add_mesh(\n",
    "    grid,\n",
    "    show_edges=True,\n",
    "    lighting=False,\n",
    "    cmap=viridis,\n",
    "    scalar_bar_args=sargs,\n",
    "    clim=[0, height],\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "318eda2e-3a8b-4f72-87ff-8f31ea2f9a47",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "plotter.view_xy()\n",
    "plotter.camera.zoom(1.3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e83a638-fd11-41e6-a1de-4c5828fcccc1",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Solving a time dependent problem"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e47c0e0a-54bb-441e-a30d-7ee4d48d596a",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    },
    "tags": []
   },
   "source": [
    "We are now ready to solve the time dependent problem. At each time step, we need to:\n",
    "1. Update the time dependent boundary condition and source\n",
    "2. Reassemble the right hand side vector `b`\n",
    "3. Apply boundary conditions to `b`\n",
    "4. Solve linear problem `Au = b`\n",
    "5. Update previous time step, `un = u`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "54d60e10-5b23-4d53-b38d-a499ba7cc67a",
   "metadata": {
    "lines_to_next_cell": 0
   },
   "outputs": [],
   "source": [
    "T = np.pi\n",
    "while t < T:\n",
    "    # Update boundary condition\n",
    "    t += dt.value\n",
    "    uD.interpolate(uD_function(t))\n",
    "\n",
    "    # Assemble RHS\n",
    "    b.x.array[:] = 0\n",
    "    petsc.assemble_vector(b.x.petsc_vec, compiled_L)\n",
    "\n",
    "    # Apply boundary condition\n",
    "    petsc.apply_lifting(b.x.petsc_vec, [compiled_a], [bcs])\n",
    "    b.x.scatter_reverse(la.InsertMode.add)\n",
    "    fem.petsc.set_bc(b.x.petsc_vec, bcs)\n",
    "\n",
    "    # Solve linear problem\n",
    "    solver.solve(b.x.petsc_vec, uh.x.petsc_vec)\n",
    "    uh.x.scatter_forward()\n",
    "\n",
    "    # Update un\n",
    "    un.x.array[:] = uh.x.array\n",
    "\n",
    "    # Update plotter\n",
    "    plotter.update_scalars(uh.x.array, render=False)\n",
    "    plotter.write_frame()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "6e6334c5-12bc-4ec9-bea9-118be535c24d",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "plotter.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "728caf7a-eebc-46b4-94f8-dcab930f330c",
   "metadata": {
    "lines_to_next_cell": 0,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "<img src=\"./u_time.gif\" alt=\"gif\" class=\"bg-primary mb-1\" width=\"800px\">"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "783d0bde-13de-41db-9335-22cd4576b345",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Post-processing without projections"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6675b27-2b78-4c71-98d7-ead71188ba2b",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    },
    "tags": []
   },
   "source": [
    "In legacy dolfin, the only way of post-processing a `ufl`-expression over the domain, would be by using a projection. This would not be scalable in most cases. Therefore, we have introduced `dolfinx.fem.Expression`, which can be used to evaluate a `ufl`-expression at any given (reference) point in any cell (local to process). Let us consider"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b196e678-f94e-497e-b181-cd5198fdafed",
   "metadata": {},
   "source": [
    "$$(y,x) \\cdot (\\nabla u)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "e8d62a1f-9a35-4969-9834-46e0760e5237",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "x = SpatialCoordinate(domain)\n",
    "x_grad = inner(as_vector((x[1], x[0])), grad(uh))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "a2e76d7f-c612-4ac0-bcaf-f93299022fb0",
   "metadata": {},
   "outputs": [],
   "source": [
    "W = fem.functionspace(domain, (\"DQ\", 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "fd0ee540-9e24-4dbb-9b63-b0d4c36f7d93",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "expr = fem.Expression(x_grad, W.element.interpolation_points())\n",
    "w = fem.Function(W)\n",
    "w.interpolate(expr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "5a42c894-ef3f-4ca9-acbb-2d20c66f53f0",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": [
     "hide-cell"
    ]
   },
   "outputs": [],
   "source": [
    "w_grid = pyvista.UnstructuredGrid(*plot.vtk_mesh(W))\n",
    "w_plotter = pyvista.Plotter(window_size=(800, 400))\n",
    "w_grid.point_data[\"w\"] = w.x.array[:].real\n",
    "w_plotter.add_mesh(w_grid, show_edges=True, cmap=viridis, scalar_bar_args=sargs)\n",
    "w_plotter.view_xy()\n",
    "w_plotter.export_html(\"./w.html\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "a6245fee-553e-40ca-a69f-6f2ddd266cff",
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<iframe src='./w.html' scrolling=\"no\" width=\"800px\" height=\"400px\"></iframe> <!--  # noqa, -->\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%html\n",
    "<iframe src='./w.html' scrolling=\"no\" width=\"800px\" height=\"400px\"></iframe> <!--  # noqa, -->"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf1452dc-c607-4cf0-99ba-8a1454cb7b0c",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "jupytext": {
   "formats": "ipynb"
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}