8. Dual-mixed formulation for Poisson equation

This demo is implemented in a single Python file, demo_mixed-poisson-dual.py, which contains both the variational forms and the solver.

This demo illustrates how to solve Poisson equation using an alternative mixed formulation. In particular, it illustrates how to

  • Use mixed and non-continuous finite element spaces
  • Set essential boundary conditions for subspaces

8.1. Equation and problem definition

A formulation of Poisson equation involves 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, where thery are referred to as Darcy flow.

After multiplying by test functions \(\tau\) and \(v\), integrating over the domain, and integrating term \(\nabla \cdot \sigma \ v\) by parts, one obtains the following variational formulation: find \(\sigma \in \Sigma\) and \(v \in V\) satisfying

\[ \begin{align}\begin{aligned}\begin{split}\int_{\Omega} (\sigma \cdot \tau + \nabla u \cdot \ \tau) \ {\rm d} x &= 0 \quad \forall \ \tau \in \Sigma, \\\end{split}\\\int_{\Omega} \sigma \cdot \nabla v \ {\rm d} x &= - \int_{\Omega} f \ v \ {\rm d} x - \int_{\Gamma_{N}} g \ v \ {\rm d} x \quad \forall \ v \in V.\end{aligned}\end{align} \]

Compared to classical mixed formulation used in demo Mixed formulation for Poisson equation, the Dirichlet condition is here essential one and Neumann condition is natural.

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 discontinuous Raviart-Thomas elements of polynomial order \(k\) and let \(V_h\) be Lagrange elements of polynomial order \(k+1\).

We will use the same definitions of functions and boundaries as in the demo for Poisson’s 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)

With the above input the solution for \(u\) and \(\sigma\) will look as follows:

../../../../_images/mixed-poisson-dual_u.png ../../../../_images/mixed-poisson-dual_sigma.png

8.2. Implementation

This demo is implemented in the demo_mixed-poisson-dual.py file.

First, the dolfin module is imported:

from dolfin import *

Then, we need to create a Mesh covering the unit square. In this example, we will let the mesh consist of 32 x 32 squares with each square divided into two triangles:

# Create mesh
mesh = UnitSquareMesh(32, 32)

Next, we need to define the function spaces. We define the two function spaces \(\Sigma_h = DRT\) and \(V_h = CG\) separately, before combining these into a mixed function space:

# Define function spaces and mixed (product) space
DRT = FunctionSpace(mesh, "DRT", 2)
CG = FunctionSpace(mesh, "CG", 3)
W = DRT * CG

The second argument to FunctionSpace specifies the type of finite element family, while the third argument specifies the polynomial degree. The UFL user manual contains a list of all available finite element families and more details. The * operator creates a mixed (product) space W from the two separate spaces DRT and CG. Hence,

\[W = \{ (\tau, v) \ \text{such that} \ \tau \in DRT, v \in CG \}.\]

Next, we need to specify the trial functions (the unknowns) and the test functions on this space. This can be done as follows

# Define trial and test functions
(sigma, u) = TrialFunctions(W)
(tau, v) = TestFunctions(W)

In order to define the variational form, it only remains to define the source functions \(f\) and \(g\). This is done just as for the mixed Poisson demo:

# Define source functions
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)")
g = Expression("sin(5.0*x[0])")

We are now ready to define the variational forms a and L.

# Define variational form
a = (dot(sigma, tau) + dot(grad(u), tau) + dot(sigma, grad(v)))*dx
L = - f*v*dx - g*v*ds

It only remains to prescribe the Dirichlet boundary condition for \(u\). Essential boundary conditions are specified through the class DirichletBC which takes three arguments: the function space the boundary condition is supposed to be applied to, the data for the boundary condition, and the relevant part of the boundary.

We want to apply the boundary condition to the second subspace of the mixed space. Subspaces of a MixedFunctionSpace can be accessed by the method sub. In our case, this reads W.sub(1). (Do not use the separate space CG as this would mess up the numbering.)

Specifying the relevant part of the boundary can be done as for the Poisson demo:

# Define Dirichlet BC
def boundary(x):
    return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS

Now, all the pieces are in place for the construction of the essential boundary condition:

bc = DirichletBC(W.sub(1), 0.0, boundary)

To compute the solution we use the bilinear and linear forms, and the boundary condition, but we also need to create a Function to store the solution(s). The (full) solution will be stored in the w, which we initialise using the FunctionSpace W. The actual computation is performed by calling solve. The separate components sigma and u of the solution can be extracted by calling the split function. Finally, we plot the solutions to examine the result.

# Compute solution
w = Function(W)
solve(a == L, w, bc)
(sigma, u) = w.split()

# Plot sigma and u
plot(sigma)
plot(u)
interactive()

8.3. Complete code


from dolfin import *

# Create mesh
mesh = UnitSquareMesh(32, 32)

# Define function spaces and mixed (product) space
DRT = FunctionSpace(mesh, "DRT", 2)
CG  = FunctionSpace(mesh, "CG", 3)
W = DRT * CG

# Define trial and test functions
(sigma, u) = TrialFunctions(W)
(tau, v) = TestFunctions(W)

# Define source functions
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)")
g = Expression("sin(5.0*x[0])")

# Define variational form
a = (dot(sigma, tau) + dot(grad(u), tau) + dot(sigma, grad(v)))*dx
L = - f*v*dx - g*v*ds

# Define Dirichlet BC
def boundary(x):
    return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS
bc = DirichletBC(W.sub(1), 0.0, boundary)

# Compute solution
w = Function(W)
solve(a == L, w, bc)
(sigma, u) = w.split()

# Plot sigma and u
plot(sigma)
plot(u)
interactive()