Mixed formulation for Poisson equation (C++)

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(div) spaces
• Define a (vector-valued) expression using additional geometry information

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.

After multiplying by test functions \(\tau\) and \(v\), integrating over the domain, and integrating the gradient term 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 \cdot \tau \ u) \ {\rm d} x &= \int_{\Gamma} \tau \cdot n \ u \ {\rm d} s \quad \forall \ \tau \in \Sigma, \\\end{split}\\\int_{\Omega} \nabla \cdot \sigma v \ {\rm d} x &= - \int_{\Omega} f \ v \ {\rm d} x \quad \forall \ v \in V.\end{aligned}\end{align} \]

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}) \text{ 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-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 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:

Implementation

The implementation is split in two files, a form file containing the definition of the variational forms expressed in UFL and the solver which is implemented in a C++ file.

Running this demo requires the files: main.cpp, MixedPoisson.ufl and CMakeLists.txt.

UFL form file

The UFL file is implemented in MixedPoisson.ufl, and the explanation of the UFL file can be found at here.

C++ program

The solver is implemented in the main.cpp file.

At the top we include the DOLFIN header file and the generated header file containing the variational forms. For convenience we also include the DOLFIN namespace.

#include <dolfin.h>
#include "MixedPoisson.h"

using namespace dolfin;

Then follows the definition of the coefficient functions (for \(f\) and \(G\)), which are derived from the DOLFIN Expression class.

// Source term (right-hand side)
class Source : public Expression
{
  void eval(Array<double>& values, const Array<double>& x) const
  {
    double dx = x[0] - 0.5;
    double dy = x[1] - 0.5;
    values[0] = 10*exp(-(dx*dx + dy*dy) / 0.02);
  }
};

// Boundary source for flux boundary condition
class BoundarySource : public Expression
{
public:

  BoundarySource(const Mesh& mesh) : Expression(2), mesh(mesh) {}

  void eval(Array<double>& values, const Array<double>& x,
            const ufc::cell& ufc_cell) const
  {
    dolfin_assert(ufc_cell.local_facet >= 0);

    Cell cell(mesh, ufc_cell.index);
    Point n = cell.normal(ufc_cell.local_facet);

    const double g = sin(5*x[0]);
    values[0] = g*n[0];
    values[1] = g*n[1];
  }

private:

  const Mesh& mesh;

};

Then follows the definition of the essential boundary part of the boundary of the domain, which is derived from the SubDomain class.

// Sub domain for essential boundary condition
class EssentialBoundary : public SubDomain
{
  bool inside(const Array<double>& x, bool on_boundary) const
  {
    return x[1] < DOLFIN_EPS or x[1] > 1.0 - DOLFIN_EPS;
  }
};

Inside the main() function we first create the mesh and then we define the (mixed) function space for the variational formulation. We also define the bilinear form a and linear form L relative to this function space.

int main()
{
  // Create mesh
  auto mesh = std::make_shared<UnitSquareMesh>(32, 32);

  // Construct function space
  auto W = std::make_shared<MixedPoisson::FunctionSpace>(mesh);
  MixedPoisson::BilinearForm a(W, W);
  MixedPoisson::LinearForm L(W);

Then we create the source (\(f\)) and assign it to the linear form.

// Create source and assign to L
auto f = std::make_shared<Source>();
L.f = f;

It only remains to prescribe the boundary condition for the flux. 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 first subspace of the mixed space. This space can be accessed through the sub member function of the FunctionSpace class.

Next, we need to construct the data for the boundary condition. An essential boundary condition is handled by replacing degrees of freedom by the degrees of freedom evaluated at the given data. The \(BDM\) finite element spaces are vector-valued spaces and hence the degrees of freedom act on vector-valued objects. The effect is that the user is required to construct a \(G\) such that \(G \cdot n = g\). Such a \(G\) can be constructed by letting \(G = g n\). This is what the derived expression class BoundarySource defined above does.

// Define boundary condition
auto G = std::make_shared<BoundarySource>(*mesh);
auto boundary = std::make_shared<EssentialBoundary>();
DirichletBC bc(W->sub(0), G, 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 Function w, which we initialise using the FunctionSpace W. The actual computation is performed by calling solve.

// Compute solution
Function w(W);
solve(a == L, w, bc);

Now, the separate components sigma and u of the solution can be extracted by taking components. These can then be written to XDMF files.

  // Extract sub functions (function views)
  Function& sigma = w[0];
  Function& u = w[1];

  XDMFFile("sigma.xdmf").write(sigma);
  XDMFFile("u.xdmf").write(u);

  return 0;
}