Poisson equation with pure Neumann boundary conditions

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

Implementation

This description goes through the implementation in demo_neumann-poisson.py of a solver for the above described Poisson equation step-by-step.

First, the dolfin module is imported:

from dolfin import *

We proceed by defining a mesh of the domain. We use a built-in mesh provided by the class UnitSquareMesh. In order to create a mesh consisting of \(64 \times 64\) squares, we do as follows:

# Create mesh
mesh = UnitSquareMesh.create(64, 64, CellType.Type_quadrilateral)

Next, we need to define the function space.

# Build function space with Lagrange multiplier
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
R = FiniteElement("Real", mesh.ufl_cell(), 0)
W = FunctionSpace(mesh, P1 * R)

The second argument to FunctionSpace specifies underlying finite element, here a mixed element is obtained by * operator.

Now, we want to define the variational problem, but first we need to specify the trial functions (the unknowns) and the test functions. This can be done using TrialFunctions and TestFunctions. It only remains to define the source function \(f\), before we define the bilinear and linear forms. It is given by a simple mathematical formula, and can easily be declared using the Expression class. Note that the string defining f uses C++ syntax since, for efficiency, DOLFIN will generate and compile C++ code for these expressions at run-time. The following code shows how this is done and defines the variational problem:

# Define variational problem
(u, c) = TrialFunction(W)
(v, d) = TestFunctions(W)
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)", degree=2)
g = Expression("-sin(5*x[0])", degree=2)
a = (inner(grad(u), grad(v)) + c*v + u*d)*dx
L = f*v*dx + g*v*ds

Since we have natural (Neumann) boundary conditions in this problem, we do not have to implement boundary conditions. This is because Neumann boundary conditions are default in DOLFIN.

To compute the solution we use the bilinear form, the 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 w, which we initialize using the FunctionSpace W. The actual computation is performed by calling solve. The separate components u and c of the solution can be extracted by calling the split function. Finally, we output the solution to a VTK file to examine the result.

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

# Save solution in VTK format
file = File("neumann_poisson.pvd")
file << u