# Mixed formulation for the Poisson equation¶

We’re considering the Poisson equation \(\nabla^2 u = -f\) using a mixed formulation on two coupled fields. We start by introducing the negative flux \(\sigma = \nabla u\) as an auxiliary vector-valued variable. This leaves us with the PDE on a unit square \(\Omega = [0,1] \times [0,1]\) with boundary \(\Gamma\)

for some known function \(f\). The solution to this equation will be some functions \(u\in V\) and \(\sigma\in \Sigma\) for some suitable function space \(V\) and \(\Sigma\) that satisfy these equations. We multiply by arbitrary test functions \(\tau \in \Sigma\) and \(\nu \in V\), integrate over the domain and then integrate by parts to obtain a weak formulation of the variational problem: find \(\sigma\in \Sigma\) and \(\nu\in V\) such that:

The flux boundary condition \(\sigma \cdot n = g\) becomes an *essential*
boundary condition to be enforced on the function space, while the boundary
condition \(u = u_0\) turn into a *natural* boundary condition which
enters into the variational form, such that the variational problem can be
written as: find \((\sigma, u)\in \Sigma_g \times V\) such that

with the variational forms \(a\) and \(L\) defined as

The essential boundary condition is reflected in function spaces \(\Sigma_g = \{ \tau \in H({\rm div}) \text{ such that } \tau \cdot n|_{\Gamma_N} = g \}\) and \(V = L^2(\Omega)\).

We need to choose a stable combination of discrete function spaces \(\Sigma_h \subset \Sigma\) and \(V_h \subset V\) to form a mixed function space \(\Sigma_h \times V_h\). One such choice is Brezzi-Douglas-Marini elements of polynomial order \(k\) for \(\Sigma_h\) and discontinuous elements of polynomial order \(k-1\) for \(V_h\).

For the remaining functions and boundaries we choose:

To produce a numerical solution to this PDE in Firedrake we procede as follows:

The mesh is chosen as a \(32\times32\) element unit square.

```
from firedrake import *
mesh = UnitSquareMesh(32, 32)
```

As argued above, a stable choice of function spaces for our problem is the
combination of order \(k\) Brezzi-Douglas-Marini (BDM) elements and order
\(k - 1\) discontinuous Galerkin elements (DG). We use \(k = 1\) and
combine the BDM and DG spaces into a mixed function space `W`

.

```
BDM = FunctionSpace(mesh, "BDM", 1)
DG = FunctionSpace(mesh, "DG", 0)
W = BDM * DG
```

We obtain test and trial functions on the subspaces of the mixed function spaces as follows:

```
sigma, u = TrialFunctions(W)
tau, v = TestFunctions(W)
```

Next we declare our source function `f`

over the DG space and initialise it
with our chosen right hand side function value.

```
x, y = SpatialCoordinate(mesh)
f = Function(DG).interpolate(
10*exp(-(pow(x - 0.5, 2) + pow(y - 0.5, 2)) / 0.02))
```

After dropping the vanishing boundary term on the right hand side, the bilinear and linear forms of the variational problem are defined as:

```
a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx
L = - f*v*dx
```

The strongly enforced boundary conditions on the BDM space on the top and bottom of the domain are declared as:

```
bc0 = DirichletBC(W.sub(0), as_vector([0.0, -sin(5*x)]), 3)
bc1 = DirichletBC(W.sub(0), as_vector([0.0, sin(5*x)]), 4)
```

Note that it is necessary to apply these boundary conditions to the first
subspace of the mixed function space using `W.sub(0)`

. This way the
association with the mixed space is preserved. Declaring it on the BDM space
directly is *not* the same and would in fact cause the application of the
boundary condition during the later solve to fail.

Now we’re ready to solve the variational problem. We define w to be a function to hold the solution on the mixed space.

```
w = Function(W)
```

Then we solve the linear variational problem `a == L`

for `w`

under the
given boundary conditions `bc0`

and `bc1`

. Afterwards we extract the
components `sigma`

and `u`

on each of the subspaces with `split`

.

```
solve(a == L, w, bcs=[bc0, bc1])
sigma, u = w.split()
```

Lastly we write the component of the solution corresponding to the primal variable on the DG space to a file in VTK format for later inspection with a visualisation tool such as ParaView

```
File("poisson_mixed.pvd").write(u)
```

We could use the built in plot function of firedrake by calling
`plot`

to plot a surface graph. Before that,
matplotlib.pyplot should be installed and imported:

```
try:
import matplotlib.pyplot as plt
except:
warning("Matplotlib not imported")
try:
fig, axes = plt.subplots()
colors = tripcolor(u, axes=axes)
fig.colorbar(colors)
except Exception as e:
warning("Cannot plot figure. Error msg '%s'" % e)
```

Don’t forget to show the image:

```
try:
plt.show()
except Exception as e:
warning("Cannot show figure. Error msg '%s'" % e)
```

This demo is based on the corresponding DOLFIN mixed Poisson demo and can be found as a script in poisson_mixed.py.