The following example explains how to solve a mixed formulation of the Poisson problem with Lagrange elements for the solution and Raviart-Thomas elements for the gradient. The full example file poisson-mfem.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 convenient ways to assemble finite element stiffness matrices, have a look at the dune-fufem and dune-fem modules.

Mixed formulation of the Poisson problem

Let $\Omega$ be a bounded open subset of $\mathbb{R}^d$ . 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}$ . In this example, we consider boundary conditions

$\begin{alignat*}{2} u & = u_0 \qquad & \text{on $\Gamma_D$} \\ \langle \nabla u, n \rangle & = g & \text{on $\Gamma_N$}, \end{alignat*}$

where $\Gamma_D$ and $\Gamma_N$ form a partition of the domain boundary $\partial \Omega$ .

The mixed formulation of this problem introduces the gradient of as a new variable $\sigma$ . One obtains a system of two first-order equations for two unknowns

$\begin{align*} \sigma - \nabla u & = 0 \\ -\operatorname{div} \sigma & = f, \end{align*}$

and the Neumann boundary condition turns into

$\‍[ \langle \sigma, n \rangle = g \qquad \text{on $\Gamma_N$.} \‍]$

Multiplying with test functions $\tau$ and (where the normal component of $\tau$ is supposed to be zero $\Gamma_n$ ) yields

$\begin{align*} \int_\Omega \langle \sigma, \tau \rangle\,dx - \int_\Omega \langle \nabla u, \tau \rangle \,dx & = 0 \\ \int_\Omega \operatorname{div} \sigma \cdot v\,dx &= -\int_\Omega fv\,dx. \end{align*}$

To obtain a symmetric bilinear form one applies integration by parts to the first equation. One gets: Find $\sigma \in H_D(\text{div})$ and $u \in L^2(\Omega)$ such that

$\begin{alignat*}{2} \int_\Omega \langle \sigma, \tau \rangle\,dx + \int_\Omega \langle \operatorname{div} \tau, u \rangle\,dx &= \int_{\Gamma_D} \langle \tau, n \rangle u\,ds & \qquad & \forall \tau \in H_0(\text{div})\\ \int_\Omega \operatorname{div} \sigma \cdot v\,dx \qquad \qquad & = - \int_\Omega fv\,dx && \forall v \in L^2(\Omega). \end{alignat*}$

Here, the space $H(\text{div})$ is defined as

$\‍[ H(\text{div}) := \Big\{ \tau \in L^2(\Omega)^d \mid \operatorname{div} \tau \in L^2(\Omega) \Big\}. \‍]$

It is a vector space of vector-valued functions, and a strict superset of $H^1(\Omega)^d$ . The space $H_D(\text{div})$ denotes the subspace of $H(\text{div})$ of vector fields complying with the boundary condition on $\Gamma_N$ . The space $H_D(\text{div})$ denotes the subspace of $H(\text{div})$ of vector fields having zero normal trace on $\Gamma_N$ .

A natural approximation space for $H(\text{div})$ is the space of Raviart-Thomas elements. It consists of particular piecewise polynomial vector fields whose normal components are continuous across element boundaries, but whose tangential components may jump. The solution lives in , and therefore piecewise constant finite elements can be used. Even though the boundary condition for is not explicitly enforced, it can be shown to hold in the limit of vanishing element diameter.

The implementation

The Python implementation starts by importing required external modules. It needs code from numpy and scipy for the sparse linear algebra, and from dune for grids, finite element bases, and quadrature rules.

import numpy
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 first creates the grid. For simplicity, the domain here is the unit square in two space dimensions, and the grid is a 50 by 50 quadrilateral grid:

upperRight = [1, 1]   # Upper right corner of the domain
elements = [50, 50]   # Number of grid elements per direction
 
# Create a grid of the unit square
gridView = grid.structuredGrid([0,0],upperRight,elements)

The next line creates the bases of the two finite element spaces: The Raviart-Thomas basis for $\sigma$ , and a zero-order Lagrange basis (i.e., piecewise constant functions) for . In dune-functions, such pairs of bases are represented in a tree data structure with one root and two leaves: Each of the leaves represents one of the two bases, and the root combines them into a single object. This is the purpose of the functions.composite method in the following code snippet. Setting up the basis tree also allows to select the way how global degrees of freedom are numbered. Here, we use the default settings, which lead to plain integers (actually, arrays of length 1) being used.

basis = functions.defaultGlobalBasis(gridView,
                                     functions.Composite(functions.RaviartThomas(order=0),
                                                         functions.Lagrange(order=0)))
 
fluxBasis     = functions.subspaceBasis(basis, 0);
pressureBasis = functions.subspaceBasis(basis, 1);

The method subspaceBasis gives access to the leaves of a composite basis. The numbers 0 and 1 passed as arguments single out the desired leaves. Note that accessing the two bases as subspace bases of a tree is not the same as simply having two separate basis objects: The actual basis functions are not affected by this, but the numbers differ.

Then, the problem stiffness matrix and right-hand side vector are assembled. The code uses two separate methods for that, which are discussed below. The source term is specified in a lambda expression f, which here encodes the constant function $f : x \mapsto 2$ .

# Compute the stiffness matrix
stiffnessMatrix = assembleMixedPoissonMatrix(basis)
 
# Compute the load vector
volumeTerm = lambda x : 2
 
b = assembleMixedPoissonRhs(basis, volumeTerm)

The next code block manages the handling of the Dirichlet boundary conditions. In the dual formulation used here, Dirichlet values are only applied to the gradient field $\sigma$ , and only to its normal components. This is easy to implement with Raviart-Thomas elements, whose degrees of freedom are exactly the normal components of the vector field.

In this example, the gradient of is prescribed on the upper and lower horizontal edge of the domain boundary. To find all degrees of freedom on these grid boundary parts, the method markBoundaryDOFsByIndicator fills the array isDirichlet with the corresponding information: The array contains a boolean variable per degree of freedom, which is set to True if a degree of freedom is attached to a grid intersection whose center is on the relevant boundary.

# This marks the top and bottom boundary of the domain
fluxDirichletIndicator = lambda x : 1.*((x[1] > upperRight[1] - 1e-8) or (x[1] < 1e-8))
 
isDirichlet = numpy.zeros(len(basis))
 
# Mark all DOFs located in a boundary intersection marked
# by the fluxDirichletIndicator function.
markBoundaryDOFsByIndicator(fluxBasis, isDirichlet, fluxDirichletIndicator);

Then, the function from the Neumann boundary condition $\langle \sigma, n \rangle = g$ is defined. In this example, it is given by the closed-form expression

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

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

normal = lambda x : numpy.array([0.,1.]) if numpy.abs(x[1]-1) < 1e-8 else numpy.array([0., -1.])
fluxDirichletValues = lambda x : numpy.sin(2.*numpy.pi*x[0]) * normal(x)
 
# TODO: This should be constrained to boundary DOFs
fluxDirichletCoefficients = numpy.zeros(len(basis))
fluxBasis.interpolate(fluxDirichletCoefficients, fluxDirichletValues);

After this, the information about the boundary condition is integrated into the linear system of equations. This happens by calling a method that appears earlier in the file:

incorporateEssentialConstraints(stiffnessMatrix, b, isDirichlet, fluxDirichletCoefficients)

The implementation of this method is

incorporateEssentialConstraints(stiffnessMatrix, b, isDirichlet, fluxDirichletCoefficients)

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.

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.

x = scipy.sparse.linalg.spsolve(stiffnessMatrix, 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:

# TODO: Improve file writing.  Currently this simply projects everything
# onto a first-order Lagrange space
vtkWriter = gridView.vtkWriter(subsampling=2)
 
fluxFunction = fluxBasis.asFunction(x)
fluxFunction.addToVTKWriter("flux", vtkWriter, grid.DataType.PointVector)
 
pressureFunction = pressureBasis.asFunction(x)
pressureFunction.addToVTKWriter("pressure", vtkWriter, grid.DataType.PointData)
 
vtkWriter.write("poisson-mfem-result")

This will write a file called poisson-mfem-result.vtu, which contains the grid and the solution in a field called "value", and the approximate gradient vector field $\sigma$ in a field called "gradient". It can be opened, for example, with the ParaView visualization program.

The assembler

The assembler is implemented in four methods above the main method of the program. Two of them implement the local assembler, i.e., they compute stiffness matrix and load vector for a single element. The two global assembler methods then piece together the local matrices and vectors into a single global matrix and vector, respectively.

Global assembler

The function that assembles the global stiffness matrix is assembleMixedPoissonMatrix. 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 elementMatrix, the element stiffness matrix. Since the basis functions of both finite element spaces are zero on most elements by construction, this element stiffness matrix is much smaller than the global stiffness matrix, but it is dense. Computation of elementMatrix happens by calling the method getLocalMatrix (discussed below). Then, the entries of the element stiffness matrix are added to the appropriate places of the global stiffnessMatrix object. The basis object provides the correspondence between the entries of the element stiffness matrices and the entries of the global stiffness matrix. Note how this does not require to distinguish between value and gradient degrees of freedom at all.

def assembleMixedPoissonMatrix(basis):
 
    # Total number of DOFs (both FE spaces together)
    n = len(basis)
 
    # Make an empty stiffness matrix
    stiffnessMatrix = lil_matrix( (n,n) )
 
    # A view on the composite FE basis on a single element
    localView = basis.localView()
 
    # A loop over all elements of the grid
    for element in basis.gridView.elements:
 
        # Bind the local basis view to the current element
        localView.bind(element)
 
        # Assemble the element stiffness matrix
        elementMatrix = getLocalMatrix(localView)
 
        # Add element stiffness matrix onto the global stiffness matrix
        for i in range(len(localView)):
 
            # The global index of the i-th local degree of freedom of the element
            row = localView.index(i)[0]  # The index is an array of length 1
 
            for j in range(len(localView)):
 
                # The global index of the j-th local degree of freedom of the element
                col = localView.index(j)[0]
                stiffnessMatrix[row,col] += elementMatrix[i, j];
 
    # Transform the stiffness matrix to CSR format, and return it
    return stiffnessMatrix.tocsr()

Note that stiffnessMatrix 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.

The next method does the global assembly for the load vector. Since a nonzero volume term only appears for the second equation, only the second basis (accessed via functions.subspaceBasis(basis,1)) is used. The volume term function given by the main method takes as argument a point in world coordinates. However, the local assembler needs a function that takes arguments in element-local coordinates. The C++ interface provides the class Dune::Functions::AnalyticGridViewFunction to do this conversion, but this is not exposed as Python yet. Instead, we project the volumeTerm onto the Lagrange finite element space. This serves the same purpose.

def assembleMixedPoissonRhs(basis, volumeTerm):
 
    # Get the basis for the pressure variable
    pressureBasis = functions.subspaceBasis(basis, 1)
 
    # Represent the volume term function as a FE function in the pressure space
    volumeTermCoeff = numpy.zeros(len(basis))
    pressureBasis.interpolate(volumeTermCoeff, volumeTerm)
    volumeTermGF = pressureBasis.asFunction(volumeTermCoeff)
 
    # A view on a single element
    localVolumeTerm = volumeTermGF.localFunction()
 
    # Create empty vector of correct size
    b = numpy.zeros(len(basis))
 
    # A view on the FE basis on a single element
    localView = basis.localView()
 
    # A loop over all elements of the grid
    for element in basis.gridView.elements:
 
        # Bind the local FE basis view to the current element
        localView.bind(element)
 
        # Get the local contribution to the right-hand side vector
        localVolumeTerm.bind(element)
        localRhs = getVolumeTerm(localView, localVolumeTerm)
 
        # Add the local vector to the global one
        for i in range(len(localRhs)):
            # The global index of the i-th degree of freedom of the element
            row = localView.index(i)[0]    # The index is an array of length 1.
            b[row] += localRhs[i]
 
    return b

Local assembler

The element stiffness matrices are computed by a free function called getLocalMatrix. Let $(\theta_i)$ be the Raviart-Thomas shape functions and $(\varphi_j)$ the zero-order Lagrange shape functions. 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.

For each element , the stiffness matrix has a two-by-two block form

$\‍[ (A_T) = \begin{pmatrix} C & B \\ B^T & 0 \end{pmatrix}, \‍]$

with entries

$\‍[ C_{ij} = \int_\Omega \langle \theta_i, \theta_j\,dx \qquad B_{ij} = \int_\Omega \operatorname{div} \theta_i, \phi_j \rangle\,dx \‍]$

Both submatrices are dense objects. In the code, they are referred to as "gradient-gradient coupling" and "value-gradient coupling", respectively.

The computation of the integrals happens via numerical quadrature, with a quadrature rule that is defined on a reference element $T_\text{ref}$ . The transformation of the divergence uses the transformation of the full Jacobian $\nabla \theta_i$ as explained in the Poisson equation (Python) example.

def getLocalMatrix(localView):
 
    # Get the grid element from the local FE basis view
    element = localView.element()
 
    # Make dense element stiffness matrix
    n = len(localView)    # Number of local degrees of freedom (gradient + value)
    elementMatrix = numpy.zeros((n,n))
 
    # Get set of shape functions for this element
    fluxLocalFiniteElement     = localView.tree().child(0).finiteElement()
    pressureLocalFiniteElement = localView.tree().child(1).finiteElement()
 
    # The actual shape functions on the reference element
    localFluxBasis = fluxLocalFiniteElement.localBasis
    localPressureBasis = pressureLocalFiniteElement.localBasis
 
    nFlux = len(fluxLocalFiniteElement)
    nPressure = len(pressureLocalFiniteElement)
 
    # Get a quadrature rule
    fluxOrder     = element.mydimension * fluxLocalFiniteElement.localBasis.order
    pressureOrder = element.mydimension * pressureLocalFiniteElement.localBasis.order
    quadOrder     = numpy.max((2*fluxOrder, (fluxOrder-1)+pressureOrder))
 
    quadRule = dune.geometry.quadratureRule(element.type, quadOrder)
 
    # Loop over all quadrature points
    for quadPoint in quadRule:
 
        # Position of the current quadrature point in the reference element
        quadPos = quadPoint.position
 
        # The inverse Jacobian of the map from the reference element to the element
        geometryJacobianInverse = element.geometry.jacobianInverse(quadPos)
 
        # The multiplicative factor in the integral transformation formula
        integrationElement = element.geometry.integrationElement(quadPos)
 
        # --------------------------------------------------------------
        #  Shape functions - gradient variable
        # --------------------------------------------------------------
 
        # Values of the gradient variable shape functions on the current element
        fluxValues = localFluxBasis.evaluateFunction(quadPos)
 
        # Gradients of the gradient variable shape functions on the reference element
        fluxReferenceJacobians = localFluxBasis.evaluateJacobian(quadPos)
 
        fluxDivergence = numpy.zeros(nFlux)
 
        # Domain transformation of Jacobians and computation of div = trace(Jacobian)
        # TODO: Extend the Dune Python interface to allow to do this without casting to numpy.array
        for i in range(nFlux):
            fluxDivergence[i] = numpy.trace(numpy.array(fluxReferenceJacobians[i]) * geometryJacobianInverse)
 
        # --------------------------------------------------------------
        #  Shape functions - value variable
        # --------------------------------------------------------------
 
        # Values of the pressure shape functions
        pressureValues = localPressureBasis.evaluateFunction(quadPos)
 
        # --------------------------------------------------------------
        #  Gradient-variable--gradient-variable coupling
        # --------------------------------------------------------------
 
        for i in range(nFlux):
 
            row = localView.tree().child(0).localIndex(i)
 
            for j in range(nFlux):
 
                col = localView.tree().child(0).localIndex(j)
 
                # Compute the actual matrix entry contribution
                elementMatrix[row,col] += numpy.dot(fluxValues[i], fluxValues[j]) * quadPoint.weight * integrationElement
 
        # --------------------------------------------------------------
        #  Flux--pressure coupling
        # --------------------------------------------------------------
 
        for i in range(nFlux):
 
            fluxIndex = localView.tree().child(0).localIndex(i)
 
            for j in range(nPressure):
 
                pressureIndex = localView.tree().child(1).localIndex(j)
 
                # Matrix entry contribution
                tmp = (fluxDivergence[i] * pressureValues[j][0]) * quadPoint.weight * integrationElement
 
                # Add to both off-diagonal block matrices
                elementMatrix[fluxIndex, pressureIndex] += tmp
                elementMatrix[pressureIndex, fluxIndex] += tmp
 
    return elementMatrix

The final method computes the load vector for one element. Since only the second equation has a nontrivial linear term, the entries corresponding to Raviart-Thomas degrees of freedom are zero. The entries for the Lagrange degrees of freedom are

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

The element load vector is computed by the following code:

def getVolumeTerm(localView, localVolumeTerm):
 
    # Get the grid element from the local FE basis view
    element = localView.element()
 
    # Construct an empty vector of the correct size
    localRhs = numpy.zeros(len(localView))
 
    # Get set of shape functions for this element
    # Only the pressure part has a non-zero right-hand side
    pressureLocalFiniteElement = localView.tree().child(1).finiteElement()
 
    # A quadrature rule
    quadOrder = 2*element.mydimension*pressureLocalFiniteElement.localBasis.order
    quadRule = dune.geometry.quadratureRule(element.type, quadOrder)
 
    nPressure = len(pressureLocalFiniteElement)
 
    # Loop over all quadrature points
    for quadPoint in quadRule:
 
        # Position of the current quadrature point in the reference element
        quadPos = quadPoint.position
 
        # The multiplicative factor in the integral transformation formula
        integrationElement = element.geometry.integrationElement(quadPos)
 
        # Evaluate the load density at the quadrature point
        functionValue = localVolumeTerm(quadPos)
 
        # Evaluate all shape function values at this point
        pressureValues = pressureLocalFiniteElement.localBasis.evaluateFunction(quadPos)
 
        # Actually compute the vector entries
        for j in range(nPressure):
            pressureIndex = localView.tree().child(1).localIndex(j)
            localRhs[pressureIndex] += - pressureValues[j][0] * functionValue * quadPoint.weight * integrationElement
 
    return localRhs

Previous	Next
Poisson equation (Python)