Tutorial

From The FEniCS project

Getting started

This all looks good, but where do I start? This is a good question and the answer depends on your interest.

You may be interested in using a particular component, like calling FIAT or SyFi to tabulate finite element basis functions, or you may be interested in using FFC as a code generator. However, you are most likely reading this because you want to solve differential equations and that's what we will concentrate on here - using FEniCS as a problem-solving environment (PSE) for differential equations.

DOLFIN functions as the main user interface of FEniCS and is available both in the form of a C++ library and as a Python module. DOLFIN depends on the other components of FEniCS as backends for specific tasks. For example, FIAT provides DOLFIN with finite element basis functions, FFC generates code for working with finite element variational formulations and Viper provides built-in plotting.

Installation

The first thing you need to do is to install the software. For a minimal installation, you need to install DOLFIN and UFC, but it is highly recommended that you also install FFC, FIAT, Instant and Viper. All these packages can be found in the Download section.

If you are lucky and are running Debian or a Debian-based system like Ubuntu, you may install all required packages through the package manager:

apt-get install fenics

This will install DOLFIN and all the other packages.

If you cannot use the prebuilt packages, you need to download the source code and build the packages manually. With the exception of DOLFIN, all packages are simple Python modules that are easy to install by

sudo python setup.py install

To build DOLFIN, use the following commands to build and install:

./configure
make
sudo make install

The configure script will search your system and remind you about any missing dependencies.

For more detailed build instructions, refer to the manuals of each package and don't be afraid to ask for help on one of the mailing lists (Contact section).

Solving Poisson's equation

If installation went smoothly and you have installed DOLFIN and its dependencies, you are now ready to solve a differential equation. The Python script listed below demonstrates how to solve Poisson's equation on the unit square and plot the solution. To run this script, save it to a file and name it demo.py. Then simply run

python demo.py

If all goes well, you should see the following plot of the computed solution.

In this script, we first create a mesh of the unit square and a finite element. We use standard piecewise linears but many other families of finite elements are available, including arbitrary order continuous and discontinuous Lagrange elements, BDM elements and RT elements. We then define a function Source for the right-hand side and a function Flux which is used for the Neumann boundary condition. We also define a subset of the boundary where we specify a Dirichlet boundary condition.

Solving Poisson's equation is then just a matter of writing down the variational problem, computing and plotting the solution. We also store the solution to a VTK file for postprocessing in MayaVi or Paraview.

The script reads as follows:

from dolfin import *

# Create mesh and finite element
mesh = UnitSquare(32, 32)
element = FiniteElement("Lagrange", "triangle", 1)

# Source term
class Source(Function):
    def __init__(self, element, mesh):
        Function.__init__(self, element, mesh)
    def eval(self, values, x):
        dx = x[0] - 0.5
        dy = x[1] - 0.5
        values[0] = 500.0*exp(-(dx*dx + dy*dy)/0.02)

# Neumann boundary condition
class Flux(Function):
    def __init__(self, element, mesh):
        Function.__init__(self, element, mesh)
    def eval(self, values, x):
        if x[0] > DOLFIN_EPS:
            values[0] = 25.0*sin(5.0*DOLFIN_PI*x[1])
        else:
            values[0] = 0.0

# Sub domain for Dirichlet boundary condition
class DirichletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return bool(on_boundary and x[0] < DOLFIN_EPS)

# Define variational problem
v = TestFunction(element)
u = TrialFunction(element)
f = Source(element, mesh)
g = Flux(element, mesh)

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

# Define boundary condition
u0 = Function(mesh, 0.0)
boundary = DirichletBoundary()
bc = DirichletBC(u0, mesh, boundary)

# Solve PDE and plot solution
pde = LinearPDE(a, L, mesh, bc)
u = pde.solve()
plot(u)

# Save solution to file
file = File("poisson.pvd")
file << u

Don't worry if this looks too simple. DOLFIN provides a simple solver for linear variational problems (which just assembles a linear system and solves it) but for more complicated problems, you may want to assemble specific matrices and vectors yourself and use these as building blocks in an iterative procedure.

For a C++ implementation of the above example, take a look at the quickstart chapter in the DOLFIN manual.