Very short introduction to FEniCS

• started in 2003, collaboration between University of Chicago and Chalmers University of Technology

• 2011 - version 1.0 released

• 2012 - the book with main contribution by 5 institutions (Simula Research Laboratory, University of Cambridge, University of Chicago, Texas Tech University, KTH Royal Institute of Technology)

  A. Logg, K.-A. Mardal, G. N. Wells et al. (2012). Automated Solution of Differential Equations by the Finite Element Method. Springer, 2012. Download from Launchpad.

• 2015 - version 1.5 released

• open source license (GNU LGPL v3), open source developement on bitbucket

• good support with mailing-list and Q&A forum

FEniCS components

Core components

• DOLFIN C++/Python interface of FEniCS, providing a consistent PSE (Problem Solving Environment).
• UFL (Unified Form Language) is a specific language for declaration of finite element discretizations of variational forms.
• FFC (FEniCS Form Compiler) from UFL code generates C++ code for assembling element tensors.
• Instant Python module that allows for instant inlining and JIT (Just-In-Time) compilation of C++ code in Python.
• FIAT (FInite element Automatic Tabulator) generates finite elements of arbitrary order on lines, triangles and tetrahedra.
• UFLACS (UFL Analyser and Compiler System) optimizing frontend for FFC.
• mshr is FEniCS mesh generator. Uses CGAL and Tetgen as backends for generating meshing geometries described by CSG (Constructive Solid Geometry).

External libraries

• MPI, OpenMP parallel programming frameworks.
• PETSc (Portable, Extensible Toolkit for Scientific Computation) parallel linear algebra backend. Provides data structeres for holding vectors and matrices and lots of linear and non-linear solvers and preconditioners.
• SLEPc (the Scalable Library for Eigenvalue Problem Computations) extension of PETSc for solution of eigen-problems.
• SCOTCH, ParMETIS mesh partitioning and graph coloring backends.
• VTK, HDF5, XDMF visualization and IO backends.

And now ...

Task 1. Start interactive Python session and type in following code.

>>> from dolfin import *

>>> mesh = UnitSquareMesh(16, 16)
>>> V = FunctionSpace(mesh, 'Lagrange', 3)
Calling FFC just-in-time (JIT) compiler, this may take some time.

We see JIT compilation of finite element code, i.e. C++ code is generated by FFC (FEniCS from compiler) and compiled by C++ compiler. This is done once and will not be done again with different mesh. The result is cached in ~/.instant.

>>> f = Expression("sin(6.0*pi*x[0])*sin(2.0*pi*x[1])")
Calling DOLFIN just-in-time (JIT) compiler, this may take some time.

This is compiled C++ expression which is evaluated very quickly when evaluated but requires JIT compilation.

>>> def boundary(x, on_boundary):
...     return on_boundary
...
>>> bc = DirichletBC(V, 0.0, boundary)

The function boundary defines boundary of the domain and bc represents Dirichlet boundary condition on space V.

>>> u = TrialFunction(V)
>>> v = TestFunction(V)
>>> a = inner(grad(u), grad(v))*dx
>>> L = f*v*dx

This code defines bilinear form a and linear form L using UFL (Unified form language).

>>> u = Function(V)
>>> solve(a == L, u, bc)
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Solving linear variational problem.

Finally we create finite-element function u on space V. (TrialFunction and TestFunction were only thought arguments of multi-linear forms - not a real function with its values in memory.) The we ask DOLFIN to assemble linear system for the respective problem and solve it by some linear algebra backend. For the former C++ code to assemble forms a and L is again generated by FFC and compiled by C++ compiler. This is again mesh independent so that it won’t be done again when refining a mesh. (Note that the solution process can be controlled in a much detailed way.)

>>> plot(u, interactive=True)

Task 2. Prepare .py file with the code above and try executing it from shell. Try also running it in parallel using mpirun command.

Task 3. Modify boundary condition on \(x=0,1\) to homogeneous Neumann.

Task 4. Modify \(-\Delta\) operator to non-linear \(-\mathrm{div} (1+k\,u^2) \nabla\) for some large \(k\). (Write linear form depending non-linearly on unknown Function and provide F == 0 instead of a == L to solve function. Use Constant class for k to avoid form recompilation when changing k.)

Reference solution

from dolfin import *

# Build mesh and function space
mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, 'Lagrange', 3)

# Right-hand side
f = Expression("sin(6.0*pi*x[0])*sin(2.0*pi*x[1])")

# Define boundary condition
def boundary(x, on_boundary):
    return on_boundary
bc = DirichletBC(V, 0.0, boundary)

# Prepare variational formulation
u = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(u), grad(v))*dx
L = f*v*dx

# Init function and solve variational problem
u = Function(V)
solve(a == L, u, bc)

# Plot solution
plot(u, interactive=True)


# Now for Neumann
def boundary(x, on_boudary):
    return near(x[1], 0.0) or near(x[1], 1.0)
bc = DirichletBC(V, 0.0, boundary)
solve(a == L, u, bc)
plot(u, interactive=True)


# And now for non-linear elliptic operator
# TODO: This implementation does not stress that linear form (without
#       TrialFunction) with non-linearity in Function is needed. But
#       attendees will probably encounter this when modifying the code.
k = 1e6
F = ( (1.0 + Constant(k)*u*u)*inner(grad(u), grad(v)) - f*v )*dx
solve(F == 0, u, bc)
plot(u, interactive=True)