Poisson equation (python)

Table of Contents

• Introduction
• Setup and initialization
• Define problem data
• Create grid and finite element basis
• Create constraints data structures
• Assemble the variational problem
• Solve the algebraic problem
• Write solution to VTK

Full code of this example This example on gitlab

Introduction

The following example shows how to solve a Poisson problem with Dirichlet boundary conditions using second order finite elements with the python interface of dune-fufem.

Setup and initialization

The example first imports the required modules from numpy, scipy, and dune and the local modules required for dune-fufems python bindings.

import numpy as np
import scipy.sparse.linalg
 
import dune.common
import dune.geometry
import dune.grid
import dune.functions
 
from dune.grid import DataType
from dune.functions import defaultGlobalBasis, Power, Composite, Lagrange
 
from fufemforms import *
from utilities import *
 

Define problem data

The example solves the Poisson equation with Dirichlet boundary conditions. The problem data consisting of the right-hand-side and the Dirichlet values and an indicator for the Dirichlet boundary are given as functions.

 
 
# Scalar right hand side
rhs = lambda x: 10
 
# Dirichlet boundary values
dirichletValues = lambda x : np.sin(2*np.pi*x[0])
 
def isNear(a,b):
  return np.abs(a-b) <= 1e-10
 
# Indicator function of Dirichlet boundary
isDirichlet = lambda x : isNear(x[0], 0) or isNear(x[0], 1) or isNear(x[1], 0) or isNear(x[1], 1)
 

Create grid and finite element basis

The example discretizes the problem using a structured grid on the unit square.

Note

There is an inconsistency in the terminology of the Python bindings with C++. What is denoted as grid in Python actually refers to the leaf GridView in the C++ interface of Dune.

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

As a next step the program creates a global finite element basis on the grid. The helper function Lagrange(order=2) declares that we want to use Lagrange finite elements of order two, while the defaultGlobalBasis(...) function creates the basis of this type on the given grid.

 
 
# Create a nodal Lagrange FE basis
basis = defaultGlobalBasis(grid, Lagrange(order=2))
 
print("Dimension of FE space is "+str(len(basis)))
 

Create constraints data structures

To handle the Dirichlet boundary conditions we first mark all constrained DOFs by interpolating the Dirichlet boundary indicator function and then interpolate the Dirichlet values on the constrained DOFs.

 
 
# Define Dirichlet constraints
N = len(basis)
 
# Mark all Dirichlet DOFs
isConstrained = np.zeros( N, dtype=bool )
basis.interpolate(isConstrained,isDirichlet)
 
# Interpolate the boundary values
constraintDOFValues = np.zeros( N )
basis.interpolate(constraintDOFValues, dirichletValues)
 

Assemble the variational problem

Inspired by the UFL language of the Fenics project project, dune-fufem allows to define local assemblers for linear and bilinear forms using simple expressions. We first obtain expressions for trial and test functions associated to the basis. The only difference between trialFunction() and testFunction() is that the resulting objects later refer to column and row indices, respectively. In order to use the right-hand-side functions within a form expression we declare it as coefficient function, i.e. a function not depending on test and trial functions. Now we can define the bilinear and linear form in a straight forward way.

 
 
# Trial and test functions for variational problem
u = trialFunction(basis)
v = testFunction(basis)
 
# Scalar rhs coefficient function
f = Coefficient(rhs, basis.gridView)
 
# Bilinear form and rhs functional
a = integrate(dot(grad(u),grad(v)))
l = integrate(f*v)
 

The defined bilinear form and linear form objects satisfy the interfaces for local matrix and vector assemblers and can thus directly be passed to a global assembler for assembling the actual matrix and right-hand-side vector. Once the system is assembled, we can constrain it to the affine subspace satisfying the constraints that we defined earlier.

Note

While the Python bindings support thread-parallel assembly this is disabled by default and not used in the example. The reason for this is, that Python itself is not thread-safe. Hence, we must not use the thread-parallel assembler when passing Python functions to the C++ assembler as we do for the right-and-side in this example.

 
 
# Assemble into matrix and vector
assembler = GlobalAssembler(basis)
A = assembler.assembleOperator(a)
b = assembler.assembleFunctional(l)
 
# Build constraints into matrix
incorporateEssentialConstraints(A, b, isConstrained, constraintDOFValues)
 

Solve the algebraic problem

To solve the linear system of equations we can use the build in sparse direct solvers of the scipy.sparse module.

 
 
# Solve linear system!
x = scipy.sparse.linalg.spsolve(A, b)
 

Write solution to VTK

Finally we can write the solution to a VTK file understood by the ParaView visualization program. To this end we first create a finite element function by combining the basis with the coefficient vector using basis.asFunction(...). Then we create a writer object with two levels of subsampling and pass the function to the writer, specifying the name and data type of the resulting field in the output.

 
u = basis.asFunction(x)
 
vtk = grid.vtkWriter(2)
u.addToVTKWriter("sol", vtk, DataType.PointData)
vtk.write("poisson-pq2-py")