# Mixed formulation for the Poisson equation

This demo illustrates how to solve Poisson equation using a mixed (two-field) formulation. In particular, it illustrates how to

• Use mixed and non-continuous finite element spaces.

• Set essential boundary conditions for subspaces and $$H(\mathrm{div})$$ spaces.

## Equation and problem definition

An alternative formulation of Poisson equation can be formulated by introducing an additional (vector) variable, namely the (negative) flux: $$\sigma = \nabla u$$. The partial differential equations then read

$\begin{split} \sigma - \nabla u &= 0 \quad {\rm in} \ \Omega, \\ \nabla \cdot \sigma &= - f \quad {\rm in} \ \Omega, \end{split}$

with boundary conditions

$\begin{split} u = u_0 \quad {\rm on} \ \Gamma_{D}, \\ \sigma \cdot n = g \quad {\rm on} \ \Gamma_{N}. \end{split}$

The same equations arise in connection with flow in porous media, and are also referred to as Darcy flow. Here $$n$$ denotes the outward pointing normal vector on the boundary. Looking at the variational form, we see that the boundary condition for the flux ($$\sigma \cdot n = g$$) is now an essential boundary condition (which should be enforced in the function space), while the other boundary condition ($$u = u_0$$) is a natural boundary condition (which should be applied to the variational form). Inserting the boundary conditions, this variational problem can be phrased in the general form: find $$(\sigma, u) \in \Sigma_g \times V$$ such that

$a((\sigma, u), (\tau, v)) = L((\tau, v)) \quad \forall \ (\tau, v) \in \Sigma_0 \times V,$

where the variational forms $$a$$ and $$L$$ are defined as

$\begin{split} a((\sigma, u), (\tau, v)) &= \int_{\Omega} \sigma \cdot \tau + \nabla \cdot \tau \ u + \nabla \cdot \sigma \ v \ {\rm d} x, \\ L((\tau, v)) &= - \int_{\Omega} f v \ {\rm d} x + \int_{\Gamma_D} u_0 \tau \cdot n \ {\rm d} s, \end{split}$

and $$\Sigma_g = \{ \tau \in H({\rm div})$$ such that $$\tau \cdot n|_{\Gamma_N} = g \}$$ and $$V = L^2(\Omega)$$.

To discretize the above formulation, two discrete function spaces $$\Sigma_h \subset \Sigma$$ and $$V_h \subset V$$ are needed to form a mixed function space $$\Sigma_h \times V_h$$. A stable choice of finite element spaces is to let $$\Sigma_h$$ be the Brezzi-Douglas-Fortin-Marini elements of polynomial order $$k$$ and let $$V_h$$ be discontinuous elements of polynomial order $$k-1$$.

We will use the same definitions of functions and boundaries as in the demo for the Poisson equation. These are:

• $$\Omega = [0,1] \times [0,1]$$ (a unit square)

• $$\Gamma_{D} = \{(0, y) \cup (1, y) \in \partial \Omega\}$$

• $$\Gamma_{N} = \{(x, 0) \cup (x, 1) \in \partial \Omega\}$$

• $$u_0 = 0$$

• $$g = \sin(5x)$$ (flux)

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

## Implementation


import numpy as np

from basix.ufl import mixed_element, element
from dolfinx import fem, io, mesh
from ufl import (Measure, SpatialCoordinate, TestFunctions, TrialFunctions,
div, exp, inner)

from mpi4py import MPI
from petsc4py import PETSc

domain = mesh.create_unit_square(
MPI.COMM_WORLD,
32, 32,
)

k = 1
Q_el = element("BDMCF", domain.basix_cell(), k)
P_el = element("DG", domain.basix_cell(), k - 1)
V_el = mixed_element([Q_el, P_el])
V = fem.FunctionSpace(domain, V_el)

(sigma, u) = TrialFunctions(V)
(tau, v) = TestFunctions(V)

x = SpatialCoordinate(domain)
f = 10.0 * exp(-((x - 0.5) * (x - 0.5) + (x - 0.5) * (x - 0.5)) / 0.02)

dx = Measure("dx", domain)
a = inner(sigma, tau) * dx + inner(u, div(tau)) * dx + inner(div(sigma), v) * dx
L = -inner(f, v) * dx

def boundary_top(x):
return np.isclose(x, 1.0)

fdim = domain.topology.dim - 1
facets_top = mesh.locate_entities_boundary(domain, fdim, boundary_top)
Q, _ = V.sub(0).collapse()
dofs_top = fem.locate_dofs_topological((V.sub(0), Q), fdim, facets_top)

def f1(x):
values = np.zeros((2, x.shape))
values[1, :] = np.sin(5 * x)
return values

f_h1 = fem.Function(Q)
f_h1.interpolate(f1)
bc_top = fem.dirichletbc(f_h1, dofs_top, V.sub(0))

def boundary_bottom(x):
return np.isclose(x, 0.0)

facets_bottom = mesh.locate_entities_boundary(domain, fdim, boundary_bottom)
dofs_bottom = fem.locate_dofs_topological((V.sub(0), Q), fdim, facets_bottom)

def f2(x):
values = np.zeros((2, x.shape))
values[1, :] = -np.sin(5 * x)
return values

f_h2 = fem.Function(Q)
f_h2.interpolate(f2)
bc_bottom = fem.dirichletbc(f_h2, dofs_bottom, V.sub(0))

bcs = [bc_top, bc_bottom]

problem = fem.petsc.LinearProblem(a, L, bcs=bcs, petsc_options={
"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"})
try:
w_h = problem.solve()
except PETSc.Error as e:
if e.ierr == 92:
print("The required PETSc solver/preconditioner is not available. Exiting.")
print(e)
exit(0)
else:
raise e

sigma_h, u_h = w_h.split()

with io.XDMFFile(domain.comm, "out_mixed_poisson/u.xdmf", "w") as file:
file.write_mesh(domain)
file.write_function(u_h)