FEniCS Introduction

There are several FEniCS projects, the necessary projects for solving simulations (FFC, FIAT, UFC, DOLFIN) and others that are tools that are nice for development (INSTANT, FERARI, SYFI, VIPER, KO, PUFFIN). Below I will go through necessary packages, as they were in the development version May 2007, assuming basic knowledge of how Finite Element Methods work.

The pieces of the FEniCS project correspond to different parts of Finite Element Modeling, and thus provide for a nice modular interfaces. Roughly the pieces correspond as:

  • FFC: Variational Form Compiler
  • FIAT: Element Quadrature Generator
  • UFC: Interface for Form Compilers to Simulation Engines
  • DOLFIN: The simulation engine

In general, a programmer only needs to write a variational form in FFC, then create a simulation in DOLFIN. FIAT and UFC serve more as backends that are called by other programs. Because of this I will only go through the details of FFC and DOLFIN explicitly with only instructions on installation for FIAT and UFC. All code used in this example is in the attached StokesInflowExample.tar.bz2 file.

Getting Started

Obtaining FEniCS

To get the development versions of the FEniCS packages, you need to install Mercurial, a Source Managment System.

Then simply do:

~/# hg clone http://www.fenics.org/hg/fiat
~/# hg clone http://www.fenics.org/hg/ufc
~/# hg clone http://www.fenics.org/hg/ffc
~/# hg clone http://www.fenics.org/hg/dolfin

Installing the Software

FFC, FIAT, and UFC follow general Python installation methods and only depend on Numpy.

~/fiat# python setup.py install --prefix={your desired install location}
~/ufc# python setup.py install --prefix={your desired install location}
~/ffc# python setup.py install --prefix={your desired install location}

DOLFIN is written in C++ and is built with the typical autotools (configure, make, make install). And requires the following (all of which probably are in your packaging system for your system).

To run things sequentially one can use UMFPack instead of PETSc and thus is a bit easier install procedure. Python bindings are not necessary but more of a convenience as well (—disable-pydolfin if you don't want them). For example a typical install for my system looks like:

~/dolfin # ./configure CXX="g++" CC="gcc" --disable-mpi --disable-petsc --disable-pydolfin --with-umfpack-include=~/UMFPACKv4.4/UMFPACK/Include --with-umfpack-lib=~/UMFPACKv4.4/UMFPACK/Lib --with-amd-lib=~/UMFPACKv4.4/AMD/Lib 
...
~/dolfin # make
...
~/doflin # make install

For a more detailed discussion on installing these packages consult Installing Dolfin or Detailed Guide to Installing DOLFIN. I personally find these guides hard to read but if you are use to compiling with autotools then things should be pretty easy.

Writing the Variational Form with FFC

To input your variational form, FEniCS can use either FFC (which calls FIAT on the backend) or SyFi (which uses symbolics as its backend). To keep things simple, I will explain the older FFC package since at the time of writing it is the more stable package.

Poisson Example

FFC uses its own very simple language for defining variational forms. Let us look at a simple example for a Poisson Equation (which can be found at ffc/src/demo/Poisson.form).

(1)
\begin{eqnarray} a(u,v) &=& \int_{\Omega}(\nabla u\cdot \nabla v) d\Omega \\ L(v) &=& (f, v) = \int_{\Omega}(f v)d\Omega \\ a(u,v) &=& L(v) \quad \forall v \in V \end{eqnarray}

with $u = 0$ on $\partial \Omega$ and $V = \{ v : v\in L^2(\Omega)^n \text{ and } v = 0 \text{ on } \partial \Omega\}$

First we need to define our element type, as such:

element = FiniteElement("Lagrange", "triangle", 1)

The FiniteElement typically consists of an element family, a geometry, and an order, but can vary depending on the element family.

Next we define the functions we are going to use in the form.

v = TestFunction(element)
u = TrialFunction(element)
f = Function(element)

The TrialFunction will be the function we wish to compute and the TestFunction is the function that varies within the function space defined by our element. The Function is a user defined function that will be inputted by the simulation engine. Finally we define the variational form $a(v,u) = L(v)$ with our defined functions.

a = dot(grad(v), grad(u))*dx 
L = v*f*dx

As you can see the there are some typical function that can be applied to the form the variational form. One nice feature here is that we can completely change spaces by changing the element and other fancy things. Here are all the parts put together:

# The bilinear form a(v, u) and linear form L(v) for
# Poisson's equation.
#
# Compile this form with FFC: ffc Poisson.form

element = FiniteElement("Lagrange", "triangle", 1)

v = TestFunction(element)
u = TrialFunction(element)
f = Function(element)

a = dot(grad(v), grad(u))*dx
L = v*f*dx

Stokes Example

The Stokes equations are an example of a mixed system, so I wanted to give you a feel for what that looks like as well, this example can be found at ffc/src/demo/Stokes.form. The variational form looks like:

(2)
\begin{eqnarray} a(u,v) - b(p,v) &=& (f,v) \quad \forall v \in V\\ b(q,u) &=& 0 \quad \forall q \in \Pi\\ \end{eqnarray}

where $u = g$ on $\partial\Omega$, $V = \{v : v\in H^1(\Omega)^n\text{, and }v=0 \text{ on }\partial\Omega\}, \Pi = L^2(\Omega)$ and

(3)
\begin{eqnarray} a(u,v) &=& \int_{\Omega}(\nabla u\cdot \nabla v) d\Omega \\ b(q,u) &=& \int_{\Omega}q (\nabla\cdot v )d\Omega \end{eqnarray}

The ffc file will look like:

# The bilinear form a(v, u) and Linear form L(v) for the Stokes
# equations using a mixed formulation (Taylor-Hood elements).

# Compile this form with FFC: ffc Stokes.form

P2 = VectorElement("Lagrange", "triangle", 2)
P1 = FiniteElement("Lagrange", "triangle", 1)
TH = P2 + P1

(v, q) = TestFunctions(TH)
(u, p) = TrialFunctions(TH)

f = Function(P2)

a = (dot(grad(v), grad(u)) - div(v)*p + q*div(u))*dx
L = dot(v, f)*dx

Notice that we can just add the two spaces together to get a mixed space, and since Stokes is a vector valued equation we use VectorElement.

Compiling and Using the Variational Form

For any form, one only needs to call:

~ # ffc -l dolfin Stokes.form

This will produce C++ code according the the UFC specification in Stokes.h with DOLFIN wrappers. When writing the simulation you will then import this into your simulation.

Writing the Simulation with Dolfin

For a very simple example, see the file dolfin/src/demo/pde/poisson/main.cpp. I'm going to skip straight to the Stokes example, since if you can write it you will see immediately what is going on in the poisson example. For our example the base file will be dolfin/src/demo/pde/stokes/taylor-hood/main.cpp, but I will change a few things to some thing a bit more simple. The example will model an inflow on a square mesh.

There are a few basic steps to writing your simulation with Dolfin:

  • Define your mesh, and boundary conditions domains,
  • Set up the PDE
  • Hand the PDE to the linear solver and get the solution

Defining the mesh

First we define a mesh, and subdomains on the boundary for applying boundary conditions. For this problem we will use a unit square and mark the top and bottom for no slip conditions, the left edge for inflow, and the right edge for outflow.

 // Read mesh and sub domain markers
  UnitSquare mesh(64,64);
  MeshFunction<unsigned int> sub_domains(mesh, mesh.topology().dim() -1);
 
  // Mark all facets as sub domain 3
  sub_domains = 3;
 
  // Mark no-slip facets as sub domain 0
  NoslipBoundary noslip_boundary;
  noslip_boundary.mark(sub_domains, 0);
 
  // Mark inflow as sub domain 1
  InflowBoundary inflow_boundary;
  inflow_boundary.mark(sub_domains, 1);
 
  // Mark outflow as sub domain 2
  OutflowBoundary outflow_boundary;
  outflow_boundary.mark(sub_domains, 2);

The UnitSquare is an easy way to get a mesh to see other methods for finding the mesh look at dolfin/src/demo/mesh. A MeshFunction is a discrete function on the mesh, discrete not as in a solution to a PDE but rather that it has values based upon mesh entities (such as vertices, edges, faces, cells …). The mesh entity is set by its dimension, since we will not have whole cells on the boundary we mark entities of dimension 1 less than a cell. So for example in 3d a tetrahedron is the cell with a dimension of 3, but faces are the highest dimension on the boundary which will be dimension 2; triangles will have cells of dimension 2 and boundaries with of dimension 1. For more on how the Dolfin mesh library, see [4].

Next we use our Subdomain definitions to mark the mesh function, thus to tell where to apply the correct boundary condition we only need to look at the value of the mesh function for each mesh entity.

// -----------------------------------------------------------------------------
// Defining the Boundary subdomains
 
// Sub domain for no-slip (everything except inflow and outflow)
class NoslipBoundary : public SubDomain
{
  bool inside(const real* x, bool on_boundary) 
  {
    return on_boundary && (x[1] < DOLFIN_EPS || x[1] > 1.0 - DOLFIN_EPS);
  }
};
 
// Sub domain for inflow (right)
class InflowBoundary : public SubDomain
{
  bool inside(const real* x, bool on_boundary)
  {
    return x[0] > 1.0 - DOLFIN_EPS && on_boundary;
  }
};
 
// Sub domain for outflow (left)
class OutflowBoundary : public SubDomain
{
  bool inside(const real* x, bool on_boundary)
  {
    return x[0] < DOLFIN_EPS && on_boundary;
  }
};

Setting up the PDE

Next to set up the PDE. Here we just call into the form that we already built. Notice we pass a function in for the RHS linear form L.

  // Set up PDE
  Function f(mesh, 0.0);
  StokesBilinearForm a;
  StokesLinearForm L(f);
  PDE pde(a, L, mesh, bcs);

But as you might see I haven't told you how to set up the boundary conditions … well just like all PDE Code this can be messy.

  // Create functions for boundary conditions
  NoslipFunction noslip_func(mesh);
  InflowFunction inflow_func(mesh);
  Function zero(mesh, 0.0);
 
  // Define sub systems for boundary conditions
  SubSystem velocity(0);
  SubSystem pressure(1);
 
  // No-slip boundary condition for velocity
  BoundaryCondition bc0(noslip_func, sub_domains, 0, velocity);
 
  // Inflow boundary condition for velocity
  BoundaryCondition bc1(inflow_func, sub_domains, 1, velocity);
 
  // Boundary condition for pressure at outflow
  BoundaryCondition bc2(zero, sub_domains, 2, pressure);
 
  // Collect boundary conditions
  Array <BoundaryCondition*> bcs(&bc0, &bc1, &bc2);

Here we define the functions for the different SubDomains on the boundary. We use the SubSystems to tell which forms the BoundaryCondition applies to, since we are using a mixed system. Then we tell the BoundaryCondition what its function,mesh function, marker, and subsystem are, gather the bcs into an array and hand them to the PDE. You might note that we do not need to set the pressure BoundaryCondition, but I do it here to demonstrate how it would be done.

To define the functions we use:

// -----------------------------------------------------------------------------
// Defining the functions on the boundaries
 
// Function for no-slip boundary condition for velocity
class NoslipFunction : public Function
{
public:
 
  NoslipFunction(Mesh& mesh) : Function(mesh) {}
 
  void eval(real* values, const real* x)
  {
    values[0] = 0.0;
    values[1] = 0.0;
  }
 
};
 
// Function for inflow boundary condition for velocity
class InflowFunction : public Function
{
public:
 
  InflowFunction(Mesh& mesh) : Function(mesh) {}
 
  void eval(real* values, const real* x)
  {
    values[0] = -1.0;
    values[1] = 0.0;
  }
 
};

Solving the PDE

Finally we give the PDE to a linear solver and then get the answer to play with.

  // Solve PDE
  Function u;
  Function p;
  pde.set("PDE linear solver", "direct");
  pde.solve(u, p);
 
  // Plot solution
  //  plot(u);
  //  plot(p);
 
  // Save solution
  File ufile("velocity.pvd");
  ufile << u;
  File pfile("pressure.pvd");
  pfile << p;

Here I use Dolfin's parameter system to set the linear solve to be "direct", which will default to UMFPack if you installed it. Otherwise it could be Petsc or uBlas, depending on your settings when you build Dolfin. The plot function is useful to see the answers in Viper if you have it installed, but I usually like to save it to VTK format and view it with Paraview

So that's all there is to it. Here is a picture of the computed velocity:

StokesInflowVelocity.png

The Full Dolfin Code

// -----------------------------------------------------------------------------
// File: main.cpp
//
// This demo solves the Stokes equations for a simple inflow/outflow problem,
// Stokes.h should define the variational form for quadratic element for the
// velocity and first degree elements for the pressure (Taylor-Hood elements).
//
// -----------------------------------------------------------------------------
 
#include <dolfin.h>
#include "Stokes.h"
 
using namespace dolfin;
 
// -----------------------------------------------------------------------------
// Defining the Boundary subdomains
 
// Sub domain for no-slip (everything except inflow and outflow)
class NoslipBoundary : public SubDomain
{
  bool inside(const real* x, bool on_boundary)
  {
    return on_boundary && (x[1] < DOLFIN_EPS || x[1] > 1.0 - DOLFIN_EPS);
  }
};
 
// Sub domain for inflow (right)
class InflowBoundary : public SubDomain
{
  bool inside(const real* x, bool on_boundary)
  {
    return x[0] > 1.0 - DOLFIN_EPS && on_boundary;
  }
};
 
// Sub domain for outflow (left)
class OutflowBoundary : public SubDomain
{
  bool inside(const real* x, bool on_boundary)
  {
    return x[0] < DOLFIN_EPS && on_boundary;
  }
};
 
// -----------------------------------------------------------------------------
// Defining the functions on the boundaries
 
// Function for no-slip boundary condition for velocity
class NoslipFunction : public Function
{
public:
 
  NoslipFunction(Mesh& mesh) : Function(mesh) {}
 
  void eval(real* values, const real* x)
  {
    values[0] = 0.0;
    values[1] = 0.0;
  }
 
};
 
// Function for inflow boundary condition for velocity
class InflowFunction : public Function
{
public:
 
  InflowFunction(Mesh& mesh) : Function(mesh) {}
 
  void eval(real* values, const real* x)
  {
    values[0] = -1.0;
    values[1] = 0.0;
  }
 
};
 
// -----------------------------------------------------------------------------
// The main simulation routine
 
int main()
{
 
  // Read mesh and sub domain markers
  UnitSquare mesh(64,64);
  MeshFunction<unsigned int> sub_domains(mesh, mesh.topology().dim() -1);
 
  // Mark all facets as sub domain 3
  sub_domains = 3;
 
  // Mark no-slip facets as sub domain 0
  NoslipBoundary noslip_boundary;
  noslip_boundary.mark(sub_domains, 0);
 
  // Mark inflow as sub domain 1
  InflowBoundary inflow_boundary;
  inflow_boundary.mark(sub_domains, 1);
 
  // Mark outflow as sub domain 2
  OutflowBoundary outflow_boundary;
  outflow_boundary.mark(sub_domains, 2);
 
  // Create functions for boundary conditions
  NoslipFunction noslip_func(mesh);
  InflowFunction inflow_func(mesh);
  Function zero(mesh, 0.0);
 
  // Define sub systems for boundary conditions
  SubSystem velocity(0);
  SubSystem pressure(1);
 
  // No-slip boundary condition for velocity
  BoundaryCondition bc0(noslip_func, sub_domains, 0, velocity);
 
  // Inflow boundary condition for velocity
  BoundaryCondition bc1(inflow_func, sub_domains, 1, velocity);
 
  // Boundary condition for pressure at outflow
  BoundaryCondition bc2(zero, sub_domains, 2, pressure);
 
  // Collect boundary conditions
  Array <BoundaryCondition*> bcs(&bc0, &bc1, &bc2);
 
  // Set up PDE
  Function f(mesh, 0.0);
  StokesBilinearForm a;
  StokesLinearForm L(f);
  PDE pde(a, L, mesh, bcs);
 
  // Solve PDE
  Function u;
  Function p;
  pde.set("PDE linear solver", "direct");
  pde.solve(u, p);
 
  // Plot solution
  //  plot(u);
  //  plot(p);
 
  // Save solution
  File ufile("velocity.pvd");
  ufile << u;
  File pfile("pressure.pvd");
  pfile << p;
}
// -----------------------------------------------------------------------------
// End of main.cpp

References

1. FEniCS Website, http://www.fenics.org
4. Logg, Ander, "Efficient representation of computational meshes", Finite Element Center preprint series.
codefemfenicsnumericaltutorial
Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-Share Alike 2.5 License.