A finite element solver for Laplace's equation

We use finite elements when the size of the data is so large that its results prohibit dealing with finite differences. To illustrate this case, we would like to explore the numerical solution of the Laplace equation, subject to certain boundary conditions.

We will start by defining the computational domain and produce a mesh dividing this domain using triangles as local finite elements. This will be our starting point to solve this problem using finite elements, as we will be placing on the computational domain a piecewise continuous function, whose pieces are linear and supported on each of the triangles.

We start by calling the necessary modules to build the mesh (other modules will be called as they are required):

>>> import numpy
>>> from numpy import linspace
>>> import scipy
>>> import matplotlib.pyplot as plt
>>> from scipy.spatial import Delaunay

First we define the region:

>>> xmin = 0 ; xmax = 1 ; nXpoints = 10
>>> ymin = 0 ; ymax = 1 ; nYpoints = 10
>>> horizontal = linspace(xmin,xmax,nXpoints)
>>> vertical = linspace(ymin,ymax,nYpoints)
>>> y, x = numpy.meshgrid(horizontal, vertical)
>>> vertices = numpy.array([x.flatten(),y.flatten()])

We may now create the triangulation:

>>> triangulation = Delaunay(vertices.T)
>>> index2point = lambda index: triangulation.points[index]
>>> all_centers = index2point(triangulation.vertices).mean(axis=1)
>>> trngl_set=triangulation.vertices

We then have the following triangulation:

>>> plt.triplot(vertices[0],vertices[1],triangles=trngl_set)
>>> plt.show()

This produces the following graph:

A finite element solver for Laplace's equation

In this case, the problem we have chosen is a standard one in mathematical methods in Physics and Engineering, consisting of solving the two-dimensional Laplace's equation on the unit square region, with zero Dirichlet boundary conditions on three sides and, on the fourth side, a constant. Physically, this problem could represent diffusion of temperature on a two-dimensional plate. Mathematically, the problem is formulated in the following form:

A finite element solver for Laplace's equation

The solution of this form can be given in terms of Fourier series as follows:

A finite element solver for Laplace's equation

This is important as you can check the correctness of the obtained numerical solution before attempting to use your numerical scheme to tackle more complex problems in complex computational domains. It should be mentioned, however, that there are alternatives in Python that implement the finite element method to solve partial differential equations. In this regard, the reader could consult the Fenics project (http://fenicsproject.org/book/) and the SfePy project (http://sfepy.org/doc-devel/index.html).

We code the solution in the usual fashion. We first compute the stiff matrix A (which for obvious reasons is sparse). Then, the construction of the vector, R, holding global boundary conditions is defined (the way we have constructed our mesh makes defining this vector straightforward). With them, the solution to the system comes from the solution X obtained from solving a matrix equation of the form AX=R using a subset of the matrices A and R corresponding to the nodes different from those on the boundaries. This should be no trouble for SciPy. Let us start with the stiff matrix:

>>> from numpy import  cross 
>>> from scipy.sparse import dok_matrix 
>>> points=triangulation.points.shape[0]
>>> stiff_matrix=dok_matrix((points,points))
>>> for triangle in triangulation.vertices:
        helper_matrix=dok_matrix((points,points))
        pt1,pt2,pt3=index2point(triangle)
        area=abs(0.5*cross(pt2-pt1,pt3-pt1))
        coeffs=0.5*numpy.vstack((pt2-pt3,pt3-pt1,pt1-pt2))/area
        #helper_matrix[triangle,triangle] = \ 
        array(mat(coeffs)*mat(coeffs).T)
        u=None 
        u=numpy.array(numpy.mat(coeffs)*numpy.mat(coeffs).T) 
        for i in range(len(triangle)):
            for j in range(len(triangle)):
                helper_matrix[triangle[i],triangle[j]] = u[i,j] 
        stiff_matrix=stiff_matrix+helper_matrix

Note that this is the cumbersome way to update the matrix stiff_matrix. This is due to the fact that the matrix is sparse, and the current choice of representation does not behave well with indexing.

To compute the global boundary vector we need to collect all edges on the boundary first and then assign to the nodes with x=1 that the function is one and to the others that the function is zero. Because of the way we set up the mesh this is easy as the nodes on which the function will take the value of one are always the last entries in the global boundary vector. This is accomplished by the following lines of code:

>>> allNodes = numpy.unique(trngl_set) 
>>> boundaryNodes = numpy.unique(triangulation.convex_hull) 
>>> NonBoundaryNodes = numpy.array([]) 
>>> for x in allNodes: 
        if x not in boundaryNodes: 
            NonBoundaryNodes = numpy.append(NonBoundaryNodes,x) 
    NonBoundaryNodes = NonBoundaryNodes.astype(int) 
    nbnodes = len(boundaryNodes) # number of boundary nodes 
    FbVals=numpy.zeros([nbnodes,1]) # Values on the boundary 
    FbVals[(nbnodes-nXpoints+1):-1]=numpy.ones([nXpoints-2, 1])

We are ready to find the numerical solution to the problem with the values obtained in our previous step:

>>> totalNodes = len(allNodes) 
>>> stiff_matrixDense = stiff_matrix.todense() 
>>> stiffNonb = \ 
    stiff_matrixDense[numpy.ix_(NonBoundaryNodes,NonBoundaryNodes)] 
>>> stiffAtb = \ 
    stiff_matrixDense[numpy.ix_(NonBoundaryNodes,boundaryNodes)] 
>>> U=numpy.zeros([totalNodes, 1]) 
>>> U[NonBoundaryNodes] = numpy.linalg.solve( - stiffNonb , \
    stiffAtb * FbVals ) 
>>> U[boundaryNodes] = FbVals

This produces the following image depicting the diffusion of temperature inside the square:

A finite element solver for Laplace's equation

This graph was obtained in the following way:

>>> X = vertices[0] 
>>> Y = vertices[1] 
>>> Z = U.T.flatten() 
>>> from mpl_toolkits.mplot3d import axes3d
>>> fig = plt.figure() 
>>> ax = fig.add_subplot(111, projection='3d') 
>>> surf = ax.plot_trisurf(X, Y, Z, cmap=cm.jet, linewidth=0) 
>>> fig.colorbar(surf) 
>>> fig.tight_layout() 
>>> ax.set_xlabel('X',fontsize=16)
>>> ax.set_ylabel('Y',fontsize=16)
>>> ax.set_zlabel(r"$\phi$",fontsize=36)
>>> plt.show()

An important point in numerical analysis is to evaluate the quality of the numerical solution obtained to any problem. In this case, we have chosen a problem whose analytical solution is available (see the preceding code), so one could check (not prove) the validity of the numerical algorithm implemented to solve our problem. In this case the analytical solution can be coded in the following manner:

>>> from numpy import pi, sinh, sin, cos, sum
>>> def f(x,y): 
        return sum( 2*(1.0/(n*pi) - \
        cos(n*pi)/(n*pi))*(sinh(n*pi*x)/ \
        sinh(n*pi))*sin(n*pi*y) 
                for n in range(1,200)) 
>>> Ze = f(X,Y) 
>>> ZdiffZe = Ze - Z 
>>> plt.plot(ZdiffZe) 
>>> plt.show()

This produces the following graph showing the difference between the exact solution (evaluated up to 200 terms) and the numerical solution of the problem (via the corresponding IPython notebook you could perform some further analysis on the numerical solution just to become more confident on the rightness of the obtained result):

A finite element solver for Laplace's equation