# Poisson equation

This demo illustrates how to:

- Create a `function space`

- Solve a linear partial differential equation

Download sources

- `Python script`

- `Jupyter notebook`

## Equation and problem definition

For a domain \(\Omega \subset \mathbb{R}^n\) with boundary \(\partial \Omega = \Gamma_{D} \cup \Gamma_{N}\), the Poisson equation with particular boundary conditions reads:

\[\begin{split} \begin{align} - \nabla^{2} u &= f \quad {\rm in} \ \Omega, \\ u &= 0 \quad {\rm on} \ \Gamma_{D}, \\ \nabla u \cdot n &= g \quad {\rm on} \ \Gamma_{N}. \\ \end{align} \end{split}\]

where \(f\) and \(g\) are input data and \(n\) denotes the outward directed boundary normal. The variational problem reads: find \(u \in V\) such that

\[ a(u, v) = L(v) \quad \forall \ v \in V, \]

where \(V\) is a suitable function space and

\[\begin{split} \begin{align} a(u, v) &:= \int_{\Omega} \nabla u \cdot \nabla v \, {\rm d} x, \\ L(v) &:= \int_{\Omega} f v \, {\rm d} x + \int_{\Gamma_{N}} g v \, {\rm d} s. \end{align} \end{split}\]

The expression \(a(u, v)\) is the bilinear form and \(L(v)\) is the linear form. It is assumed that all functions in \(V\) satisfy the Dirichlet boundary conditions (\(u = 0 \ {\rm on} \ \Gamma_{D}\)).

In this demo we consider:

- \(\Omega = [0,2] \times [0,1]\) (a rectangle)

- \(\Gamma_{D} = \{(0, y) \cup (2, y) \subset \partial \Omega\}\)

- \(\Gamma_{N} = \{(x, 0) \cup (x, 1) \subset \partial \Omega\}\)

- \(g = \sin(5x)\)

- \(f = 10\exp(-((x - 0.5)^2 + (y - 0.5)^2) / 0.02)\)

## Implementation

The modules that will be used are imported:

    from pathlib import Path

    from mpi4py import MPI
    from petsc4py.PETSc import ScalarType  # type: ignore

    import numpy as np

    import ufl
    from dolfinx import fem, io, mesh, plot
    from dolfinx.fem.petsc import LinearProblem

Note that it is important to first `from mpi4py import MPI` to ensure that MPI is correctly initialised.

We create a rectangular `Mesh` using `create_rectangle`, and create a finite element `function space` \(V\) on the mesh.

    msh = mesh.create_rectangle(
        comm=MPI.COMM_WORLD,
        points=((0.0, 0.0), (2.0, 1.0)),
        n=(32, 16),
        cell_type=mesh.CellType.triangle,
    )
    V = fem.functionspace(msh, ("Lagrange", 1))

The second argument to `functionspace` is a tuple `(family, degree)`, where `family` is the finite element family, and `degree` specifies the polynomial degree. In this case `V` is a space of continuous Lagrange finite elements of degree 1. For further details of how one can specify finite elements as tuples, see `ElementMetaData`.

To apply the Dirichlet boundary conditions, we find the mesh facets (entities of topological co-dimension 1) that lie on the boundary \(\Gamma_D\) using `locate_entities_boundary`. The function is provided with a ‘marker’ function that returns `True` for points `x` on the boundary and `False` otherwise.

    facets = mesh.locate_entities_boundary(
        msh,
        dim=(msh.topology.dim - 1),
        marker=lambda x: np.isclose(x[0], 0.0) | np.isclose(x[0], 2.0),
    )

We now find the degrees-of-freedom that are associated with the boundary facets using `locate_dofs_topological`:

    dofs = fem.locate_dofs_topological(V=V, entity_dim=1, entities=facets)

and use `dirichletbc` to create a `DirichletBC` class that represents the boundary condition:

    bc = fem.dirichletbc(value=ScalarType(0), dofs=dofs, V=V)

Next, the variational problem is defined:

    u = ufl.TrialFunction(V)
    v = ufl.TestFunction(V)
    x = ufl.SpatialCoordinate(msh)
    f = 10 * ufl.exp(-((x[0] - 0.5) ** 2 + (x[1] - 0.5) ** 2) / 0.02)
    g = ufl.sin(5 * x[0])
    a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
    L = ufl.inner(f, v) * ufl.dx + ufl.inner(g, v) * ufl.ds

A `LinearProblem` object is created that brings together the variational problem, the Dirichlet boundary condition, and which specifies the linear solver. In this case an LU solver is used, and we ask that PETSc throws an error if the solver does not converge. The `solve` computes the solution.

    problem = LinearProblem(
        a,
        L,
        bcs=[bc],
        petsc_options_prefix="demo_poisson_",
        petsc_options={"ksp_type": "preonly", "pc_type": "lu", "ksp_error_if_not_converged": True},
    )
    uh = problem.solve()
    assert isinstance(uh, fem.Function)

The solution can be written to a `XDMFFile` file visualization with ParaView or VisIt:

    out_folder = Path("out_poisson")
    out_folder.mkdir(parents=True, exist_ok=True)
    with io.XDMFFile(msh.comm, out_folder / "poisson.xdmf", "w") as file:
        file.write_mesh(msh)
        file.write_function(uh)

and displayed using pyvista.

    try:
        import pyvista

        cells, types, x = plot.vtk_mesh(V)
        grid = pyvista.UnstructuredGrid(cells, types, x)
        grid.point_data["u"] = uh.x.array.real
        grid.set_active_scalars("u")
        plotter = pyvista.Plotter()
        plotter.add_mesh(grid, show_edges=True)
        warped = grid.warp_by_scalar()
        plotter.add_mesh(warped)
        if pyvista.OFF_SCREEN:
            plotter.screenshot(out_folder / "uh_poisson.png")
        else:
            plotter.show()
    except ModuleNotFoundError:
        print("'pyvista' is required to visualise the solution.")
        print("To install pyvista with pip: 'python3 -m pip install pyvista'.")