General Concepts (including how to solve a Laplace problem)

In the following we introduce how to construct and use the basic components of a finite element method ending the first part with solving a simple Laplace problem. The required steps are:

  • constructing a tessellation of the computational domain

  • working with functions defined over the grid

  • setting up discrete functions

  • defining the mathematical model to solve

  • solving the (non linear) system arising from the discretization of the model by the Galerkin method

After that we study the more complex problem of solving a time dependent non linear problem.

Setting up the Mesh

After some general import statements we start our tutorial by setting up a simple Cartesian grid over the domain \([0,1]^2\) subdivided into four intervals in each dimension. We can then show the grid using the plot method - note that direct plotting with MatPlotLib is only supported in 2D, other options for visualization, e.g., using ParaView will be discussed later.

Using the dune.grid.structuredGrid function, is the simplest way to construct a grid, more complicated (unstructured) grids can be constructed using a dictionary containing the vertices and the element connectivity or by for example providing a gmsh file. This will be discussed in later section, see for example here for a 3d gmsh file or here for a simple example using a dictonary.

import io, numpy, matplotlib
matplotlib.rc( 'image', cmap='jet' )
from matplotlib import pyplot

from dune.grid import structuredGrid as leafGridView
gridView = leafGridView([0, 0], [1, 1], [4, 4])

Grid Functions

These are function that are defined over a given grid and are evaluated by using an element of the grid and local coordinated within that element:

value = gridFunction(element,localCoordinate)

Alternatively one can obtain a LocalFunction from a grid function which can be bound to an element and then evaluate via local coordinate:

localFunction = gridFunction.localFunction()
for e in grid.elements:
    value = localFunction(x)

There are multiple ways to construct grid functions. The easiest way it to use UFL expression. Many methods expecting grid functions are argument can also directly handle UFL expression. We can for example integrate a UFL expression over the grid:

from ufl import SpatialCoordinate, triangle
x = SpatialCoordinate(triangle)

exact = 1/2*(x[0]**2+x[1]**2) - 1/3*(x[0]**3 - x[1]**3) + 1

from dune.fem.function import integrate
print( integrate(gridView, exact, order=5) )

and plot them using matplotlib or write a vtk file for postprocessing

from dune.fem.plotting import plotPointData as plot
plot(exact, grid=gridView)
gridView.writeVTK('exact', pointdata={'exact': exact})

In some cases it can be necessary to convert a UFL expression into a grid function explicitly - for example to be able to evaluate it over each element in the grid (in a later section) we provide more detail on how to access the grid interface).

from dune.fem.function import uflFunction
exact_gf = uflFunction(gridView, name="ufl", order=1, ufl=exact)
mass = 0
for element in gridView.elements:
    mass += exact_gf(element,[0.5,0.5]) * element.geometry.volume

Another way to obtain a grid function is to use the gridFunction decorator. This can be obtained from dune.grid but then without UFL support. Using the decorator from dune.fem.function the resulting grid function can be used seamlessly within UFL expressions:

from dune.fem.function import gridFunction
def exactLocal(element,xLocal):
    x = element.geometry.toGlobal(xLocal)
    return 1/2.*(x[0]**2+x[1]**2) - 1/3*(x[0]**3 - x[1]**3) + 1

we can use the same approach but with a function using global coordinates but can then be used like any other grid function:

def exactGlobal(x):
    return 1/2.*(x[0]**2+x[1]**2) - 1/3*(x[0]**3 - x[1]**3) + 1

lf = exactGlobal.localFunction()

mass = 0
for element in gridView.elements:
    mass += lf([0.5,0.5]) * element.geometry.volume

print( integrate(gridView, [exact,exactLocal,exactGlobal], order=5) )
(1.333333, 1.333333, 1.333333)

As pointed out the dune.fem grid function can be used like any other UFL coefficient to form UFL expressions:

print( integrate(gridView, abs(exact-exactLocal), order=5) )
gf = uflFunction(gridView, name="ufl", order=1, ufl=exact+exactLocal*exactGlobal)
fig = pyplot.figure(figsize=(20,10))

Converting UFL expressions to grid functions leads to JIT code generation and is therefore efficient when used in other C++ algorithm (like integrate). On the other hand using the gridFunction decorator leads to a callback into Python for each evaluation and is therefore much less efficient. An alternative approach is based on writing small C++ snippets implementing the grid function:

from dune.fem.function import cppFunction
#include <cmath>
#include <dune/common/fvector.hh>
template <class GridView>
auto aTimesExact(double a) {
  // the return value needs to be a lambda which either returns a `double`
  // or a `Dune::FieldVector<double,R>` with `R>=1`
  return [a](const auto& en,const auto& xLocal) -> auto {
    auto x = en.geometry().global(xLocal);
    return a*(1./2.*(std::pow(x[0],2)+std::pow(x[1],2)) - 1./3.*(std::pow(x[0],3) - std::pow(x[1],3)) + 1.);
exactCpp = cppFunction(gridView, name="exactCpp", order=2,
print( integrate(gridView, abs(2*exact-exactCpp), order=5) )

As the above example shows it is easy to pass in parameters to the C++ implementation - here a double 2. Note that it is not possible to obtain a reference to this double so to make sure a change to the constant on the Python side carries over to the C++ side an option is to use a dune.common.FieldVector instead. These parameters can also be more complex e.g. other grid # function - note that a UFL expression can not be directly passed in - it first needs to be converted into a grid function using uflFunction.

#include <cmath>
#include <dune/common/fvector.hh>
template <class GridView, class GF>
auto aTimesExact(const GF &gf,Dune::FieldVector<double,1> &a) {
  return [lgf=localFunction(gf),&a] (const auto& en,const auto& xLocal) mutable -> auto {
    lgf.bind(en); // lambda must be mutable so that the non const function can be called
    return a[0]*lgf(xLocal);
from dune.common import FieldVector
a = FieldVector([2])
exactCpp2 = cppFunction(gridView, name="exactCpp", order=2,
print( integrate(gridView, abs(exactCpp-exactCpp2), order=5) )
a[0] = 4
print( integrate(gridView, abs(2*exactCpp-exactCpp2), order=5) )

The above is just one of a few ways C++ code snippets can be used together with Python code to improve efficiency or extend the existing binding to Dune. In the above example there is no advantage of using the C++ code over the code generated based on the UFL expression. For more complicated functions e.g. with many if statements or based on more information from the given element like it’s neighbors the expressibility of UFL might not be sufficient or lead to hard to read code. In these cases directly providing C++ code (or Python code) can be a reasonable alternative.

Discrete Functions

Note that the grid functions set up so far did not involve any discretization, they are exactly evaluated at the given point. A special type of grid functions are discrete functions living in a discrete (finite dimensional) space. .. index:: Functions; discrete functions

from import lagrange as solutionSpace
space = solutionSpace(gridView, order=2)

The easiest way to construct a discrete function is to use the interpolate method on the discrete function space.

u_h = space.interpolate(exact, name='u_h')

On an existing discrete function the interpolate method can be used to reinitialize it

u_h.interpolate( cppFunction(gridView, name="exactCpp", order=2,
                 args=[2.]) )
u_h.interpolate( lambda x: 1/2.*(x[0]**2+x[1]**2) - 1/3*(x[0]**3 - x[1]**3) + 1 )

Note that in the last example the Python lambda is used as a callback automatically using the same concept used in the gridFunction decorator. As pointed out above there are some methods where these conversions are implicit and no explicit generation of a grid function has be carried out.

If a discrete function is already available it is possible to call copy to obtain further discrete functions:

u_h_n = u_h.copy(name="previous")

Finally, clear can be called on a discrete function which sets all coefficient to zero and assign can be used to copy all coefficients between two discrete function over the same space:

u_h_n.assign( u_h )

All the things we did above with grid functions can be done with discrete functions, e.g., evaluate locally

localUh = u_h.localFunction()
mass = 0
for element in gridView.elements:
    localUh.bind(element) # using u_h(element,[0.5,0.5]) also works
    mass += localUh([0.5,0.5]) * element.geometry.volume

or plot using matplotlib and write a vtk file for postprocessing (using binary data format to reduce size)

from dune.grid import OutputType
gridView.writeVTK('uh', pointdata=[u_h], outputType=OutputType.appendedraw)

Note: the discrete function u_h already has a name attribute given in the interpolate call. This is used by default in the vtk file. An alternative name can be given by using a dictionary as shown previously.

Of course a discrete function can also be used as a coefficient in a UFL expression

print( integrate(gridView, abs(exact-u_h), order=5) )

The main difference between a grid function and a discrete function is that the latter has a dofVector attached to it. We can iterate over that easily

for d in u_h.dofVector:
    print(d,end=", ")
1.0, 1.0260416666666667, 1.0833333333333333, 1.140625, 1.1666666666666667, 1.0364583333333333, 1.0625, 1.1197916666666667, 1.1770833333333333, 1.203125, 1.1666666666666667, 1.1927083333333333, 1.25, 1.3072916666666667, 1.3333333333333335, 1.421875, 1.4479166666666665, 1.5052083333333335, 1.5625, 1.5885416666666667, 1.8333333333333333, 1.859375, 1.9166666666666665, 1.9739583333333333, 2.0, 1.0071614583333333, 1.052734375, 1.1139322916666667, 1.1595052083333333, 1.0436197916666667, 1.0891927083333333, 1.150390625, 1.1959635416666667, 1.173828125, 1.2194010416666667, 1.2805989583333333, 1.326171875, 1.4290364583333333, 1.474609375, 1.5358072916666665, 1.5813802083333335, 1.8404947916666665, 1.8860677083333333, 1.947265625, 1.9928385416666665, 1.0084635416666667, 1.0345052083333333, 1.091796875, 1.1490885416666667, 1.1751302083333335, 1.087890625, 1.1139322916666667, 1.1712239583333333, 1.228515625, 1.2545572916666667, 1.2766927083333333, 1.302734375, 1.3600260416666667, 1.4173177083333333, 1.443359375, 1.6061197916666665, 1.6321614583333333, 1.689453125, 1.7467447916666665, 1.7727864583333335, 1.015625, 1.0611979166666667, 1.1223958333333333, 1.16796875, 1.0950520833333333, 1.140625, 1.2018229166666667, 1.2473958333333333, 1.2838541666666665, 1.3294270833333333, 1.390625, 1.4361979166666667, 1.61328125, 1.6588541666666665, 1.7200520833333333, 1.765625,

In a later section we will see how to extract the underlying dof vector in the form of a numpy or a petsc vector.

Models and Schemes

We consider a scalar boundary value problem \begin{align*} -\triangle u &= f & \text{in}\;\Omega:=(0,1)^2 \\ \nabla u\cdot n &= g_N & \text{on}\;\Gamma_N \\ u &= g_D & \text{on}\;\Gamma_D \end{align*} and \(f=f(x)\) is some forcing term. For the boundary conditions we set \(\Gamma_D={0}\times[0,1]\) and take \(\Gamma_N\) to be the remaining boundary of \(\Omega\).

We will solve this problem in variational form \begin{align*} \int \nabla u \cdot \nabla \varphi \ - \int_{\Omega} f(x) \varphi\ dx - \int_{\Gamma_N} g_N(x) v\ ds = 0. \end{align*} We choose \(f,g_N,g_D\) so that the exact solution is given by \begin{align*} u(x) = \left(\frac{1}{2}(x^2 + y^2) - \frac{1}{3}(x^3 - y^3)\right) + 1 \end{align*} Note: in a later section we discuss more general boundary conditions.

from ufl import TestFunction, TrialFunction
from dune.ufl import DirichletBC
u = TrialFunction(space)
v = TestFunction(space)

from ufl import dx, grad, div, grad, dot, inner, sqrt, conditional, FacetNormal, ds
a = dot(grad(u), grad(v)) * dx

f   = -div( grad(exact) )
g_N = grad(exact)
n   = FacetNormal(space)
b   = f*v*dx + dot(g_N,n)*conditional(x[0]>=1e-8,1,0)*v*ds
dbc = DirichletBC(space,exact,x[0]<=1e-8)

With the model described as a ufl form, we can construct a scheme class that provides the solve method which we can use to compute the solution:

from dune.fem.scheme import galerkin as solutionScheme
scheme = solutionScheme([a == b, dbc], solver='cg')
scheme.solve(target = u_h)
{'converged': True, 'iterations': 1, 'linear_iterations': 49}

We can compute the error between the exact and the discrete solution by using the integrate function described above:

h1error = dot(grad(u_h - exact), grad(u_h - exact))
error = sqrt(integrate(gridView, h1error, order=5))
print("Number of elements:",gridView.size(0),
      "number of dofs:",space.size,"H^1 error:", error)
Number of elements: 16 number of dofs: 81 H^1 error: 0.00658807845868413

To verify that the discrete scheme is converging to the exact solution we can compute the experimental order of convergence (EOC): \begin{align*} {\rm eoc} = \frac{\log{e_h/e_H}}{\log{h/H}} \end{align*} where \(h,H\) refer to the spacing of two grids and \(e_h,e_H\) are measures for the accuracy of the discrete solution, i.e., \(e_h=\|u-u_h\|,e_H=\|u-u_h\) using a suitable norm (here \(h_1\) as used before):

import dune.fem
loops = 2
for eocLoop in range(loops):
    error_old = error
    # gridView.hierarchicalGrid.globalRefine(1)
    scheme.solve(target = u_h)
    error = sqrt(integrate(gridView, h1error, order=5))
    eoc = round(numpy.log(error/error_old)/numpy.log(0.5),2)
          "Number of elements:", gridView.size(0),
          "number of dofs:", space.size,"H^1 error:", error)
EOC: -0.0 Number of elements: 16 number of dofs: 81 H^1 error: 0.0065880784587151
EOC: -0.0 Number of elements: 16 number of dofs: 81 H^1 error: 0.0065880784587151

We have already seen how grid function can be used within UFL expressions and forms - in the next section we will give another example for this in the context of a time dependent problems. In addition we can also use the Constant class to add constants to UFL expressions/forms which then can be changed without requiring any recompilation of the model. An example would again be in a time dependent problem a time varying coefficient. Being able to change the value of the time in the model without recompilation is crucial for an efficient code. We will demonstrate this here by adding a mass term with a constant mass \(m\). We will not change the right hand side so the exact solution we used so far will only be valid for \(m=0\).

from dune.ufl import Constant

m = Constant(1,name="mass")   # start with m=1
a = dot(grad(u), grad(v)) * dx + m*u*v * dx
scheme = solutionScheme([a == b, dbc], solver='cg')

scheme.solve(target = u_h)

h1error = dot(grad(u_h - exact), grad(u_h - exact))
error = sqrt(integrate(gridView, h1error, order=5))
print("Number of elements:",gridView.size(0),
      "number of dofs:",space.size,"H^1 error:", error)
Number of elements: 16 number of dofs: 81 H^1 error: 0.5635954310234561

We can print the value of a Constant with name foo via and change it’s value using the same attribute:

scheme.model.mass = 0      # go back to original problem

scheme.solve(target = u_h)

h1error = dot(grad(u_h - exact), grad(u_h - exact))
error = sqrt(integrate(gridView, h1error, order=5))
print("Number of elements:",gridView.size(0),
      "number of dofs:",space.size,"H^1 error:", error)
Number of elements: 16 number of dofs: 81 H^1 error: 0.0065880784586842085


Parallelization is available using either distributed memory based on MPI or multithreading using OpenMP.


It is straightforward to use MPI for parallelization. It requires a parallel grid, in fact most of the DUNE grids work in parallel except albertaGrid or polyGrid. Most iterative solvers in DUNE work for parallel runs. Some of the preconditioning methods also work in parallel. Running a parallel job can be done by using mpirun

mpirun -np 4 python

in order to use 4 MPI processes. Example scripts that run in parallel are, for example, the Re-entrant Corner Problem.


It is straightforward to enable some multithreading support. Note that this will speedup assembly and evaluation of the spatial operators but not in general the linear algebra backend so that the overall speedup of the code might not be as high as expected. Since we rely mostly on external packages for the linear algebra, speedup in this step will depend on the multithreading support available in the chosen linear algebra backend. We discuss in a later section how to switch between different implementation for solving linear systems of equations, e.g., how to use the thread parallel solvers from scipy.

By default only a single thread is used. To enable multithreading simply add

from dune.fem import threading
threading.use = 4 # use 4 threads
Using 1 threads
Using 4 threads

At startup a maximum number of threads is selected based on the hardware concurrency. This number can be changed by setting the environment variable DUNE_NUM_THREADS which sets both the maximum and set the number of threads to use. To get this number or set the number of threads to use to the maximum use

print("Maximum number of threads available:",threading.max)
Maximum number of threads available: 12
Using 12 threads

Listing Available Dune Components

The available realization of a given interface, i.e., the available grid implementations, depends on the modules found during configuration. Getting access to all available components is straightforward:

from dune.utility import components
# to get a list of all available components:
# to get for example all available grid implementations:
available categories are:
available entries for this category are:
entry         function          module
agglomerate   polyGrid          dune.vem.vem
alberta       albertaGrid       dune.grid
alu           aluGrid           dune.alugrid
aluconform    aluConformGrid    dune.alugrid
alucube       aluCubeGrid       dune.alugrid
alusimplex    aluSimplexGrid    dune.alugrid
oned          onedGrid          dune.grid
sp            spIsotropicGrid   dune.spgrid
spanisotropic spAnisotropicGrid dune.spgrid
spbisection   spBisectionGrid   dune.spgrid
spisotropic   spIsotropicGrid   dune.spgrid
ug            ugGrid            dune.grid
yasp          yaspGrid          dune.grid
[ ]:

Available discrete function spaces are:

available entries for this category are:
entry          function       module
bbdg           bbdgSpace      dune.vem.vem
bdm            bdm  
combined       combined
dglagrange     dglagrange
dglegendre     dglegendre
dglegendrehp   dglegendrehp
dgonb          dgonb
dgonbhp        dgonbhp
finitevolume   finiteVolume
lagrange       lagrange
lagrangehp     lagrangehp
p1bubble       p1Bubble
product        product
rannacherturek rannacherTurek
raviartthomas  raviartThomas
vem            vemSpace       dune.vem.vem

Available discrete functions are:

available entries for this category are:
entry      function          module
cpp        cppFunction       dune.fem.function
discrete   discreteFunction  dune.fem.function
global     globalFunction    dune.fem.function
levels     levelFunction     dune.fem.function
local      localFunction     dune.fem.function
partitions partitionFunction dune.fem.function
ufl        uflFunction       dune.fem.function

Available schemes:

available entries for this category are:
entry      function         module
bbdg       bbdgScheme       dune.vem.vem
burgers    burgers          dune.fem.scheme
dg         dg               dune.fem.scheme
dggalerkin dgGalerkin       dune.fem.scheme
galerkin   galerkin         dune.fem.scheme
h1         h1               dune.fem.scheme
h1galerkin h1Galerkin       dune.fem.scheme
linearized linearized       dune.fem.scheme
rungekutta rungeKuttaSolver dune.femdg
stokes     stokes           dune.fem.scheme
vem        vemScheme        dune.vem.vem

This page was generated from the notebook concepts_nb.ipynb and is part of the tutorial for the dune-fem python bindings DOI