Introduction to Python for FEniCS

Introduction

redirect the "WIP" (work in progress) wiki page Egm6936.s11/FENICS WIP to a more permanent wiki page. Egm6321.f12 (talk) 14:15, 3 August 2012 (UTC)

find the expansion for FEniCS (FE = Finite Element, what does "niCS" stand for ?) Egm6321.f12 (talk) 14:15, 3 August 2012 (UTC)

Origin of the name FEniCS

FEniCS is an acronym with FE representing finite element, CS representing computational software, and according to Anders Logg, a Senior Research Scientist with the FEniCS Project, "ni sits nicely in the middle". Andy Terrel also notes that the FEniCS software package was originally compiled at the University of Chicago, whose mascot is a phoenix, which likely inspired the name. A discussion of the topic can be found here on the FEniCS launchpad.

Article contents

This article acts as an introduction to the Python programming language and specifically its use with FEniCS.

The Python code below, taken from an example problem found here:Egm6936.s11/FENICS/Python, is explained on a line-by-line basis.

The structure of the example code can be applied to many problems but will require varying degrees of modification.

This article is designed to be self-reliant but an official Python tutorial can be found here: http://docs.python.org/tutorial/

The FEniCS Fundamentals page covers similar material, in addition to much more, found here: http://fenicsproject.org/documentation/tutorial/fundamentals.html

The official DOLFIN library can be found here: DOLFIN Library

Creating a Python code file

A Python code file can quickly be created from the Linux terminal. Entering the command "gedit" will open a blank document. Simply write the code in this document and save it into the desired directory with a .py extension. The file can be reopened for editing by entering "gedit filename.py". The code can be run by entering "Python filename.py".

Code can also be written using iPython, an interactive Python with many benefits and additional capabilities. Begin iPython by entering "ipython" in the terminal. You can quit iPython at any time by entering "quit". The code can be saved to your current directory by entering "save filename.py line#" where "line#" represents the lines you wish to save. For example, "save example.py 1-7 9 12" would save lines 1 through 7, 9, and 12 to a new file named example.py.

Full example code

#!/usr/bin/python
from dolfin import *

#define mesh and function space
mesh=IntervalMesh(20, 0, 1)
V = FunctionSpace(mesh, "Lagrange", 1)

#define left boundary
class BoundaryLeft(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and \
            (near(x[0], 0.))

#define right boundary
class BoundaryRight(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and \
            (near(x[0], 1.))

#define inputs
u_0 = Constant(0.)
u_1 = Constant(1.)

#initialize subdomains
boundaryLeft = BoundaryLeft()
boundaryRight = BoundaryRight()

#apply essential boundary conditions
bcs = [DirichletBC(V, u_0, boundaryLeft),
       DirichletBC(V, u_1, boundaryRight)
      ]

#define functions
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0.)
g = Constant(0.)
nu = Constant(1.)
ac = Constant(1.)
 
#define problem
a = nu*dot(grad(u), grad(v))*dx+ac*u*v*dx
L = f*v*dx + g*v*ds

#solve problem
u = Function(V)
solve(a == L, u, bcs)
 
#plot solution
plot(u)

Methodology

Importing classes via the import statement

from dolfin import *

Code explanation

The code here imports classes like Interval, Functionspace, Function, etc. from the DOLFIN library. "dolfin" is a python package allowing access to the C++ software library for finite element computing, DOLFIN. FEniCS relies heavily on these classes so this will normally be the first line in your Python code. Including "from" causes the statement to repeatedly define names from the module. The "*" tells the statement to continue until all names in the module have been defined, hence, the full statement imports and defines all classes in the dolfin module.

read the "Elements of Style" by Strunk and White (there is also an online 1918 version, an older version). never use the article "This" without a noun that follows it, since "This" by itself is very vague; it could stand for anything in the above-mentioned list of objects. Egm6321.f12 (talk) 14:42, 3 August 2012 (UTC)

Python manual reference

give a link to the online python manual regarding the python statements "from" and "import". explain also the meaning of the "*".

see The Python Language Reference, section 6.12 on the "import" statement.

from a Unix user perspective, the "*" represents "everything", so that the above code line may mean to import all existing classes from the libray "dolfin" into the present python code. you need to verify this guess.

Egm6321.f12 (talk) 14:42, 3 August 2012 (UTC)

The import statement, including "from" and "*".

Defining the mesh and functionspace classes

the wikipedia style (and the styles of orther publishers) is to capitalize only the first character in a title. you don't have to capitalize the first characters of every word in a title; see the above title. Egm6321.f12 (talk) 14:50, 3 August 2012 (UTC)

mesh=Interval(20, 0, 1)
V = FunctionSpace(mesh, "Lagrange", 1)

Code explanation

The code here first defines your mesh class as another DOLFIN class, Interval, following the syntax "Interval"(divisions, minimum, maximum). The example interval runs from 0 to 1 with 20 divisions.

explain also what "mesh" is as a python object, i.e., is it a class, a variable, etc. ? also what is "Interval" ? a function ? where can we find its code listing ? Egm6321.f12 (talk) 14:50, 3 August 2012 (UTC)

We then define a discrete function space V over the mesh. The functionspace class follows the syntax "FunctionSpace"(your mesh, type of element, degree of basis functions). Above, the standard Lagrange family of elements is used which can also be denoted by 'CG' for Continuous Galerkin. A degree of 1 is also shown. This denotes the standard linear Lagrange element, later yielding a computed u that is continuous and linearly varying over each mesh division.

Additional mesh definitions

The following lines can be used to create a unit square or rectangular mesh. Each cell created by the x and y divisions is divided in half to create two triangular elements.

mesh="Unitsquare"(x-divisions, y-divisions)

mesh="Rectangle"(x-min, y-min, x-max, y-max, x-divisions, y-divisions)

Defining classes for boundary conditions

class BoundaryLeft(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and \
            (near(x[0], 0.))

class BoundaryRight(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and \
            (near(x[0], 1.))

Code explanation

Here the code defines new classes, BoundaryLeft and BoundaryRight, for the left and right boundaries. Note that variables x and y would be denoted by x[0] and x[1], respectively. We define the function "inside" to return true for all x points both on the boundary (on_boundary) and first 0 (near(x[0], 0.), then 1 (near(x[0], 1.). The function "near"(value1,value2) checks that value1 is within machine precision of value2.

Defining values at boundaries

u_0 = Constant(0.)
u_1 = Constant(1.)

Code explanation

Here we create variables of the DOLFIN class Constant to assign as our values at the boundaries. The problem statement specifies that u(0)=0 and u(1)=1, demonstrated above.

Initializing subdomains for boundaries

boundaryLeft = BoundaryLeft()
boundaryRight = BoundaryRight()

Code explanation

The code here creates two variables of the previously created class types BoundaryLeft() and BoundaryRight(). These variables are used in the next step when applying the boundary conditions.

Applying essential boundary conditions using DirichletBC

bcs = [DirichletBC(V, u_0, boundaryLeft),
       DirichletBC(V, u_1, boundaryRight)
      ]

Code explanation

Here we apply the Dirichlet boundary conditions using the dolfin class DirichletBC following the syntax "DirichletBC"(function space, boundary value, boundary variables). Note that the left and right boundaries are both included.

Defining TrialFunction and TestFunction classes and given constants

u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0.)
g = Constant(0.)
nu = Constant(1.)
ac = Constant(1.)

Code explanation

Here we define several variables including the dolfin classes TrialFunction and TestFunction using the function space V. Constants f, g, nu, and ac are also defined to values given in the problem statement.

Defining the variational problem

a = nu*inner(grad(u), grad(v))*dx+ac*u*v*dx
L = f*v*dx + g*v*ds

Code explanation

Here we define our a(u,v) and L(v), obtained by the derivation of the weak variational form. Note the syntax for inner product inner and gradient grad.

Solving the problem

u = Function(V)
solve(a == L, u, bcs)

Code explanation

Here the first line creates a function, u, over the function space V. Having defined a, L, and information about our essential boundary conditions, we use solve to compute the solution, a finite element function u.

Plotting the solution

plot(u)

Code explanation

Here we plot the solution, u, seen below. Including interative=True allows the user to manipulate the plot and is necessary to keep it on screen.