Python and Finite Elements

Initially, I designed ELLIPT2D so that users could leverage a scripting interface by specifying geometry and PDE coefficients on the fly. This provided flexibility that cannot be matched in compiled environments. Of course, this idea isn't new. FREEFEM (http://www.freefem.org/) comes with its own scripting capability targeted for FE problems, and FEMLAB (http://www.femlab.com/) is invoked from the Matlab shell. Unfortunately, it is difficult to design a program that takes into account all potential needs. In particular, I encountered cases where the PDE coefficients led to singularities, which, to be properly treated, required users to modify ELLIPT2D's internals. In a problem involving photonics, for instance, the boundary conditions were such that the value at a node was equal to the value at another node times a phase shift. Such boundary conditions cannot be applied in a general geometry, but are perfectly reasonable in the considered geometry. But again, the problem could only be resolved after having access to the source and knowing the internals of ELLIPT2D.

To make ELLIPT2D user friendly, powerful, and transparent to the user, I realized it would have to be written in a scripting language such as Python.

The FE Method

The aforementioned physics and engineering problems can be cast in the general linear second order PDE form as in Figure 1, where F, g, and s are user-defined functions of the space coordinates, and v the solution to be determined. In the following, we assume that the problem, due to symmetry, can be reduced to two dimensions so that these functions depend on the coordinates x and y only. To obtain a unique solution, boundary conditions must be applied on the boundary of the domain (), as in Figure 2 where and are also user-prescribed functions. Here, = [/x, /y] denotes the Cartesian Del operator and n the normal vector pointing outwards of (n· would be the normal derivative). F in Figure 1 is a tensor, so that F· is the vector [F_xx/x+F_xy/y, F_yx/x+F_yy/y]. Commonly occurring elliptic PDEs can be derived from Figure 1 by suitable choices of F(x,y), g(x,y), and s(x,y). For example, ²v+k²v=0 (the Helmholtz equation describing the propagation of stationary waves) can be obtained simply by setting F=I (the identity matrix), g= -k² and s=0. Likewise, the same equation can be written in cylindrical geometry as F=xI, g=-k²x after interpreting x as the radial coordinate.

Given a problem like Figure 1 with the boundary conditions in Figure 2, the FE method involves the following steps:

1. Mesh generation. The domain is discretized in small cells. In the case of ELLIPT2D, these are triangles with each cell bounded by three nodes (vertices) going counterclockwise; see, for example, the triangulated mesh in Figure 3.

2. The solution v(x)=V_je_j(x,y) is expanded in a suitable set of basis functions e_j(x,y) or finite elements, so called because they have finite support. In ELLIPT2D, I use finite elements that are linear polynomials of x and y in each cell. By finite support, I mean that the amplitude of a given element is zero on all nodes except one, where it is unity. This ensures that the amplitude V_j coincides with the numerical value of v(x,y) on node j. Figure 3 is one such tent element.

3. The equation in Figure 1 is multiplied by all basis functions e_i(x,y) in turn and integrated over . Integration by parts then yields a discretized form of Figure 1 — namely the system of linear equations in Figure 4, where Figure 5 yields the stiffness matrix A_ij. Because of the finite support of the basis functions, A_ij is sparse and efficient methods exist to solve the linear system. The right-hand side vector b_i contains the source and a term representing the contribution from the boundary conditions. This step is referred to as the "matrix assembly." In ELLIPT2D, F(x,y), g(x,y), and s(x,y) are assumed to vary linearly between nodes leading to the integrals in Figure 5, which depend solely on the node values of these functions.

4. Finally, the system in Figure 4 is solved for the unknown amplitudes V_j. This is by far the most demanding numerical step. Accordingly, ELLIPT2D offers the choice between a Python conjugate gradient method and a C interface to the sparse matrix solver SuperLU (http://www.nersc.gov/~xiaoye/SuperLU/) for increased performance.

Implementation Issues

Granted, Python is not the first programming language that comes to mind for heavy number crunching: Matlab, Scilab, and the Interactive Data Language offer excellent interfaces to sparse matrix solvers. However, I chose Python because many data structures used in the FE method do not neatly fit into a matrix cast, the core data type of these languages. Furthermore, Python is object oriented, has a particularly clean syntax for nested lists and exception handling, and allows operator to be overloaded.

ELLIPT2D makes extensive use of three types of objects: node, cell, and sparse. The node object is central to the FE method as it contains information on how nodes are connected, that is, how finite elements couple in the stiffness matrix. This is captured in Python by a dictionary data type whose syntax is {key: value,...}. Dictionaries, also known as "hash tables" or "associative arrays," are powerful objects allowing values to be accessed by keys. They can be thought of as generalized arrays, where values are indexed not by integers, but by other objects. (One difference is that the dictionary's layout in memory is unordered.) Indeed, in Python, the syntax to access a list and dictionary value is the same and based on the [] operator. To determine the node connectivity, a dictionary {i1:[(x,y), [ia, ib, ic,...],a1],...} is used containing i1 as the node index, (x,y) a tuple of node coordinates, [ia,ib,ic,...] a list of node indices that are connected to i1, and a1 — an additional attribute set to 1 if i1 is a boundary node. The dictionary values can be almost anything (in this case, a nested list) while the keys can only be immutable objects (a key cannot be altered). This illustrates the utilization of a dictionary as a replacement for a C structure with the major advantage that the dictionary can grow dynamically (both the elements and the number of elements) with no particular memory-management responsibilities attached. This property is particularly useful in the case of unstructured meshes where the number of connected nodes is not known a priori.

Also of importance in ELLIPT2D is the cell object, a nested list [k,[i1,i2,i3],...], where k is the cell index and [i1,i2,i3] are the three nodes ordered in counterclockwise direction. The cell object is used during the assembly time. You can see from Figure 5 how the stiffness matrix integral naturally breaks into the sum of all cell contributions.

So far, I have not mentioned how the boundary conditions are taken care of. The FE method distinguishes between natural and essential boundary conditions. The mixed (Robbins) boundary conditions in Figure 2 are natural and dealt with by replacing n·F· v in the source term b_i in Figure 5 by -v and moving the term involving v to the stiffness matrix. However, Dirichlet boundary conditions where the value of v is imposed on need special treatment. Imposing the value for node i is achieved in ELLIPT2D by setting A_ij=0 except for A_ii=1, while b_i takes the prescribed value. ELLIPT2D has a DirichletBound dictionary object that takes the node index as key. The RobbinsBound and NeumannBound (similar to RobbinsBound, but for =0) take the tuple (i,j) as keys since they act on segments, but are otherwise invoked similarly to DirichletBound after the matrix assembly step.

Finally, the sparse stiffness matrix is represented in ELLIPT2D by the object sparse, a dictionary {(i,j): aij,...} with the row/column indices (i,j) and aij as key/value pairs, respectively. Only nonzero elements need to be stored. Python's object-oriented features let the usual matrix operations (+,-,*, and so on) be overloaded. To invert the matrix system, a conjugate gradient method using overloaded matrix operations has been written in 30 lines of Python code.

The Laplace Equation

To illustrate how ELLIPT2D works, I'll present an application that solves the Laplace equation on a rectangular mesh (Figure 6) with nx1 nodes in the x direction and ny1 nodes in the y direction; see Listing One(a). The node object grid is constructed by calling one of many available methods that break a rectangular domain into a regular set of triangles. For more general structured meshes, grid can be constructed from a set of interleaving (x,y) coordinates; see Listing One(b).

Next, boundary conditions must be applied. These are a sine wave on the west edge, and zero Dirichlet conditions on the south and north edges. On the east edge, zero-flux Neumann boundary conditions are imposed ( and set to zero in Figure 2). Since zero-flux conditions are the default boundary conditions of the FE method, there is no need to create a NeumannBound object for the west edge; see Listing One(c).

The ellipt2d object can now be built; see Listing One(d), and the stiffness matrix computed, as in Listing One(e). The functions F, g, and s, which are passed as arguments to the ellipt2d constructor, are strings. These can be any expressions of x and y. To avoid any singularity, be sure that F remains positive definite and that all functions are finite on the nodes. Thus, you could just as well have chosen F="2 +cos(pi*x**2*y)", for instance. For problems where these functions cannot be expressed in closed form, ELLIPT2D offers in addition the option to supply F, g, or s as node vectors, instead of strings. The vector object vector.zeros(N) is equivalent to supplying "0." on all nodes. Also note that a scalar form of F was fed into ellipt2d — that is, "1." is formally equivalent to the identity matrix (("1.","0."), ("0.","1.")).

The stiffness matrix and right-hand side vector are updated through a call, Listing One(f), and in the final step, Figure 4 is solved by invoking the conjugate gradient (CG) method; Listing One(g). The CG method is iterative in nature and requires an initial guess v0. The next arguments after the right side vector b specify the maximum tolerance and the maximum number of iterations. Methods for saving the results in various formats — toUCD for AVSExpress, toDX for openDX; see Listing One(h) — are available for graphical postprocessing.

Resonant Frequencies of a Structure

As a second test application, consider the problem of determining the resonant frequencies of a structure consisting of an infinite set of conducting wires arranged in a regular, periodic manner. Because of the periodicity, it is possible to focus on one cell, shown in Figure 7. In such a cell, it can be shown that the problem reduces to computing in Figure 8, where dm and dn are input parameters. Notice the boundary conditions, which are such that the amplitudes on the north edge are equal, up to a complex phase factor (i is the imaginary number), to those on the south edge, and likewise for the east and west amplitudes. Here, you're mainly interested in the lowest eigenvalue .

Complications arise in several respects. First, a circular concavity representing the wire's section must be carved out. This means that an unstructured mesh must be constructed, which is done by importing the module Ireg2tri, a wrapping class around the C-library TRIANGLE (http://www.cs.cmu.edu/~quake/triangle.html). More significantly, to build the eigenvalue problem, the stiffness matrix must be split into two matrices (called "amat" and "bmat" in code available electronically; see "Resource Center," page 5). Boundary conditions are applied by manually modifying the elements that couple to a boundary node. For the phased-shifted conditions, this is done by adding the north and south, and east and west elements. In the case of corner elements, the contributions must be added from all other corners. Finally, the eigenvalue is computed using an inverse iteration scheme.

Conclusion

Of course, performance is a drawback of scripting languages. But for problems involving up to several thousand nodes, this wasn't an issue. The Laplace problem (also available electronically) runs in 1.5 CPU seconds for 100 nodes in pure Python mode, whereas the complex eigenvalue problem (800 nodes) takes less than 30 CPU seconds on a 600-MHz Pentium III laptop. Based on simple Laplace equation tests, I found Jython (http://www.jython.org/) to run at roughly half Python's speed (Sun's Java 1.3.1). If the CPU time is broken up between modules, then the conjugate gradient solver is the most expensive part by far, with the matrix assembly requiring less than 10 percent of the CPU time. Hence, I recommend the use of the C interface to the SuperLU library. A global speedup factor of 10 can then be expected. Both the matrix assembly and SuperLU CPU's time requirements scale proportionally to N^1.7, where N is the number of nodes. A typical ELLIPT2D problem involves up to a few thousand nodes. Accordingly, ELLIPT2D cannot presently be extended to three-dimensional geometry where problems routinely have hundreds of thousands of nodes, although the programming paradigm applied here would lend itself to such extensions.

All is not perfect even in the best of all worlds. It was slightly disappointing to observe that more recent versions of Python run slower: Python 2.0 is roughly 15 percent slower than Version 1.6, which is slower than 1.5.2 by a few percents. From a programmer's point of view, I would have preferred that complex and real (floating) numbers be treated on the same footing. For instance, a.real generates an exception if a is a float, although you almost certainly want it to mean "If a is complex, then take the real part". This may simply be the result of complex numbers being added quite recently to the language. Indeed, Python's biggest strength is that so many interesting features (object orientation, exception handling, namespace) were built in the language from the ground up.

Acknowledgments

The boundary condition classes (DirichletBound, NeumannBound, and RobbinsBound), the C extension to TRIANGLE, the HTML documentation, as well as the Java graphics class are all contributions by John C. Mollis. This work was supported by the Department of Energy and through a National Undergraduate Fellowship.

DDJ

Listing One

import reg2tri, vector       
from ellipt2d import *

xmin, xmax = 0., 1.4
ymin, ymax = 0., 1.

N = 101
nx1 = int(sqrt(float(N)*(xmax-xmin)/(ymax-ymin)))
ny1 = int(sqrt(float(N)*(ymax-ymin)/(xmax-xmin)))

(b)
<pre>grid = reg2tri.rect2criss((xmin, ymin, xmax, ymax), nx1, ny1)

(c)
<pre>db = DirichletBound()               
for i in range(nx1):
    db[i]             = 0.0  # south
    db[nx1*(ny1-1)+i] = 0.0   # north

for i in range(ny1):
    y = grid.y(i*nx1)
    db[i*nx1] = sin(pi*(y-ymin)/(ymax-ymin))   # west

(d)
<pre>F='1.'
g, s = '0.', '0.'
equ = ellipt2d(grid, F, g, s)

(e)
<pre>[amat, b] = equ.stiffnessMat()

(f)
<pre>equ.dirichletB(db,amat,b)

(g)
<pre>v0 = b
v = amat.CGsolve(v0, b, 1.e-6, 2*len(b))

(h)
<pre>equ.toUCD(v, 'ex1.inp')

Back to Article