The following example explains how to solve the Poisson problem with Lagrange finite elements using the Python interface of dune-functions. The full example file poisson-pq2.py can be found in the examples/ subdirectory.

Be warned that writing a quadrature loop for a stiffness matrix with dune-functions and Python is currently quite slow. This approach is therefore mainly useful for small problems, for educational purposes, or if you need complete control. For faster and more efficient ways to assemble finite element stiffness matrices, have a look at the dune-fufem and dune-fem modules.

The Poisson problem

Let $ \Omega $ be a bounded open subset of $ \mathbb{R}^2 $ or $ \mathbb{R}^3 $. The Poisson problem asks for a scalar function $ u : \Omega \to \mathbb{R} $ that solves

\[ -\Delta u = f \qquad \text{on $\Omega$} \]

for a given function $f : \Omega \to \mathbb{R} $. To make the solution unique, it is additionally required that the solution $ u $ satisfies boundary conditions

\[ u = g \qquad \text{on $\partial \Omega$}. \]

The finite element method replaces the partial differential equation by a so-called weak formulation, and solves it in a space $ V_h $ of continuous, piecewise polynomial functions: The goal is now to find a function $ u_h \in V_h $ still satisfying the boundary conditions, such that

\[ \int_\Omega \langle \nabla u_h, \nabla v_h \rangle \,dx = \int_\Omega f v_h\,dx \qquad \forall v_h \in V_{h,0}, \]

where $ V_{h,0} $ is the space of finite element functions that are zero on the boundary. Writing this weak form in terms of a basis $ \{ \varphi_i \} $ of $ V_h $, one obtains a sparse linear system of equations

\[ Ax = b, \]

which can be solved using standard algorithms from numerical linear algebra.

The implementation

As every Python program, this one starts by importing required external modules. In this case, we need code from numpy and scipy for the sparse linear algebra, and from dune for grids, finite element bases, and quadrature rules.

import numpy as np
import scipy.sparse.linalg
from scipy.sparse import lil_matrix
 
import dune.geometry
import dune.grid as grid
import dune.functions as functions

Main program

The main program starts by creating a grid. For simplicity, the domain here is the unit square in two space dimensions, and the grid is a 32 by 32 quadrilateral grid:

# Number of grid elements in one direction
gridSize = 32
 
# Create a grid of the unit square
grid = grid.structuredGrid([0,0],[1,1],[gridSize,gridSize])

The next line creates the basis of the finite element space. We use second-order Lagrange finite elements here, i.e., the space $ V_h $ consists of continuous functions that are quadratic in each variable on each element. The class functions.Lagrange encodes the nodal basis of that space:

# Create a second-order Lagrange FE basis

basis = functions.defaultGlobalBasis(grid, functions.Lagrange(order=2))

Then, the Poisson problem is assembled. The source term $ f $ is specified in a lambda expression f, which here encodes the constant function $ f : x \mapsto 10 $. The actual assembly then happens in the method assembleLaplaceMatrix, which we discuss below.

# Source term
f = lambda x : 10
 
# Compute stiffness matrix and load vector
A,b = assembleLaplaceMatrix(basis, f)

The assembleLaplaceMatrix method returns two objects, the stiffness matrix A and the assembled source term vector b.

The next code block manages the handling of the Dirichlet boundary conditions. In this example, the entire boundary is the Dirichet boundary. To find all degrees of freedom on the grid boundary, the method forEachBoundaryDOF visits each degree of freedom (DOF) associated to the grid boundary and calls the callback function markDOF with the global index of that degree of freedom. This callback then sets the corresponding entry in an array called isDirichlet to True.

# Determine all coefficients that are on the boundary
isDirichlet = np.zeros(len(basis))
 
def markDOF(boundaryDOFNumber):
    isDirichlet[boundaryDOFNumber] = True
 
functions.boundarydofs.forEachBoundaryDOF(basis,markDOF)

Then, the function $ g $ from the problem formulation is defined. In this example the Dirichlet values are given by the closed-form expression

\[ g(x_0,x_1) = \sin(2\pi x_0). \]

This closed form expression is evaluated at the boundary Lagrange nodes by the method basis.interpolate.

# The function that implements the Dirichlet values
dirichletValueFunction = lambda x : np.sin(2*np.pi*x[0])
 
# Get coefficients of a Lagrange-FE approximation of the Dirichlet values
dirichletValues = np.zeros(len(basis))
basis.interpolate(dirichletValues, dirichletValueFunction)

After this, the information about the boundary conditions is integrated into the linear system of equations. The code loops over all nonzero entries of the matrix. If an entry belongs to a matrix row that corresponds to a boundary degree of freedom, it is set to 1 if it is the diagonal entry, and to 0 otherwise. This pins the degree of freedom to the corresponding boundary value. This part of the code is completely independent from Dune.

# Integrate Dirichlet conditions into the matrix and load vector
rows, cols = A.nonzero()
 
for i,j in zip(rows, cols):
    if isDirichlet[i]:
        if i==j:
            A[i,j] = 1.0
        else:
            A[i,j] = 0
        b[i] = dirichletValues[i]

At this point, the linear system of equations has been completely assembled, and can now be solved. For this, the code simply calls a sparse direct solver from the Python scipy package. This is not part of Dune, either:

# Solve linear system!

x = scipy.sparse.linalg.spsolve(A, b)

The solver leaves the algebraic solution vector in the vector variable x. The solution is then written to a file, in the VTK file format. To do so, the solution vector must be reinterpreted as a finite element functions again. The asFunction method does this:

# Write result as vtu file
uh = basis.asFunction(x)
 
vtkWriter = grid.vtkWriter(subsampling=2)
uh.addToVTKWriter("u", vtkWriter, dune.grid.DataType.PointData)
vtkWriter.write("poisson-pq2-result")

This will write a file called poisson-pq2-result.vtu, which contains the grid and the solution function $ u_h $. It can be opened, for example, with the ParaView visualization program.

The assembler

The assembler is implemented in two methods above the main method of the program.

Global assembler

The central function is assembleLaplaceMatrix, which computes the stiffness matrix $ A $ with entries

\[ A_{ij} = \int_\Omega \langle \nabla \varphi_i, \nabla \varphi_j \rangle dx \]

and the load vector $ b $ with entries

\[ b_i = \int_\Omega f \varphi_i\,dx. \]

As described in many text books, this assembly is implemented as a loop over the grid elements. On each element, all relevant integrals are computed and stored in a matrix localA, the element stiffness matrix. Since the basis functions $ \varphi_i $ are zero on most elements by construction, this element stiffness matrix is much smaller than the global stiffness matrix A, but it is dense. Computation of localA happens by calling yet another method localAssembler (discussed below). Afterwards, the entries of the element stiffness matrix are added to the appropriate places of the global stiffness matrix A. The basis object provides the correspondence between the entries of the element stiffness matrices and the entries of the global stiffness matrix.

def assembleLaplaceMatrix(basis, volumeTerm):
 
    # Total number of degrees of freedom
    n = len(basis)
 
    # Make empty sparse matrix
    A = lil_matrix((n,n))
 
    # Make empty vector
    b = np.zeros(n)
 
    # View on the finite element basis on a single element
    localView = basis.localView()
 
    # Loop over all grid elements
    grid = basis.gridView
    for element in grid.elements:
 
        # Bind the localView to the current element
        localView.bind(element)
 
        # Number of degrees of freedom on the current element
        localN = len(localView)
 
        # Assemble the local stiffness matrix and load vector
        localA, localb = localAssembler(localView, volumeTerm)
 
        # Copy the local entries into the global matrix and vector
        for i in range(localN):
 
            gi = localView.index(i)[0]
 
            b[gi] += localb[i]
 
            for j in range(localN):
                gj = localView.index(j)[0]
                A[gi, gj] += localA[i, j]
 
    # Convert matrix to CSR format
    return A.tocsr(), b

Note that A is sparse (i.e., it contains mainly zeros). This property has to be exploited. Therefore, a lil_matrix object from the Python scipy package is used to build it. After assembly has been completed, the matrix is converted to Compressed Sparse Row (CSR) format, because that is more efficient for the solver.

Local assembler

The element stiffness matrices are computed by a free function called localAssembler. For a given element $ T $, it computes the element stiffness matrix

\[ (A_T)_{ij} = \int_T \langle \nabla \varphi_i, \nabla \varphi_j \rangle \,dx \]

and the element load vector $ b_T $ with entries

\[ (b_T)_i = \int_T f \varphi_i\,dx. \]

Both are dense objects. In this example, where the domain is two-dimensional, the grid consists of quadrilaterals, and the finite element space consists of second-order Lagrange elements, each $ A_T $ is a $ 9 \times 9$ matrix. The element and the set of basis functions on the element are passed to the method via a localView object. See the Dune::Functions::DefaultLocalView page for documentation of the corresponding C++ interface.

The computation of the integrals happens via numerical quadrature, and the quadrature rule is defined on a reference element $ T_\text{ref} $. The integral therefore needs to be transformed to the reference element, using a mapping $ F : T_\text{ref} \to T $:

\[ (A_T)_{pq} = \int_{T_\text{ref}} \Big \langle (\nabla_\xi \hat{\varphi}_{T,p}(\xi)) \cdot (\nabla_\xi F_T(\xi))^{-1}, (\nabla_\xi \hat{\varphi}_{T,q}(\xi)) \cdot (\nabla_\xi F_T(\xi))^{-1} \Big\rangle | \det \nabla_\xi F_T(\xi)|\,d\xi. \]

The integrals over $ T_\text{ref} $ are then computed by a quadrature rule.

def localAssembler(localView, volumeTerm):
 
    # Number of degrees of freedom on this element
    n = len(localView)
 
    # The grid element
    element = localView.element()
 
    # The set of shape functions
    localBasis = localView.tree().finiteElement().localBasis
 
    # Initialize element matrix and load vector
    localA = np.zeros((n,n))
    localB = np.zeros(n)
 
    # Create a quadrature rule and integrate
    quadRule = dune.geometry.quadratureRule(element.type, order=4)
    for pt in quadRule:
 
        # The determinant that appears in the integral transformation formula
        integrationElement = element.geometry.integrationElement(pt.position)
 
        # Evaluate all shape functions (The return value is an array!)
        values = localBasis.evaluateFunction(pt.position)
 
        # Evaluate the shape function Jacobians on the reference element (array of arrays)
        referenceJacobians = localBasis.evaluateJacobian(pt.position)
 
        # Transform the reference Jacobians to the actual element
        geometryJacobianInverse = element.geometry.jacobianInverse(pt.position)
        jacobians = [ np.dot(np.array(g)[0], geometryJacobianInverse) for g in referenceJacobians ]
 
        quadPosGlobal = element.geometry.toGlobal(pt.position)
 
        for i in range( n ):
            for j in range( n ):
                localA[i,j] += pt.weight * integrationElement * np.dot(jacobians[i], jacobians[j])
 
            localB[i] += pt.weight * integrationElement * values[i] * volumeTerm(quadPosGlobal)
 
    return localA, localB

Previous	Next
Poisson equation (C++)	Mixed FE discretization of the Poisson equation (Python)