Run this tutorial

Click here to run this tutorial on mybinder.org: try on mybinder.org
Please note that starting the notebook server may take a couple of minutes.

Tutorial: continuous FEM for the stationary heat equation#

This tutorial shows how to solve the stationary heat equation using continuous Finite Elements with dune-gdt.

This is work in progress [WIP], still missing:

  • Neumann boundary values

  • Robin boundary values

%matplotlib notebook

import numpy as np
np.warnings.filterwarnings('ignore') # silence numpys warnings

diffusion with homogeneous Dirichlet boundary condition#

analytical problem#

Let ΩRd for 1d3 be a bounded connected domain with Lipschitz-boundary Ω. We seek the solution uH01(Ω) of the linear diffusion equation (with a homogeneous Dirichlet boundary condition)

(1)(κu)=fin Ω,u=0on Ω,

in a weak sense, where κ[L(Ω)]d×d denotes a given diffusion function and fL2(Ω) denotes a given source function.

The variational problem associated with (1) reads: find uH01(Ω), such that

(2)a(u,v)=l(v)for all vH01(Ω),

where the bilinear form a:H1(Ω)×H1(Ω)R and the linear functional lH1(Ω) are given by

(3)a(u,v):=Ω(κu)vdxandl(v):=Ωfdx,

respectively.

Consider for example (1) with:

  • d=2

  • Ω=[0,1]2

  • κ=1

  • f=exp(x0x1)

from dune.xt.grid import Dim
from dune.xt.functions import ConstantFunction, ExpressionFunction

d = 2
omega = ([0, 0], [1, 1])

kappa = ConstantFunction(dim_domain=Dim(d), dim_range=Dim(1), value=[1.], name='kappa')
# note that we need to prescribe the approximation order,
# which determines the quadrature on each element
f = ExpressionFunction(
    dim_domain=Dim(d), variable='x', expression='exp(x[0]*x[1])', order=3, name='f')

continuous Finite Elements#

Let us for simplicity consider a simplicial grid Th (other types of elements work analogously) as a partition of Ω into elements KTh with h:=maxKThdiam(K), we consider the discrete space of continuous piecewise polynomial functions of order kN,

Vh:={vC0(Ω)|v|KPk(K)}

where Pk(K) denotes the space of polynomials of (total) degree up to k (note that VhH1(Ω) and Vh does not include the Dirichlet boundary condition, thus VhH01(Ω.). We obtain a finite-dimensional variational problem by Galerkin-projection of (2) onto Vh, thas is: we seek the approximate solution uhVhH01(Ω), such that

(4)a(uh,vh)=l(vh)for all vhVhH01(Ω).

A basis of Vh is given by the Lagrangian basis functions

Φ:={φ1,,φN}

of order k (e.g., the usual hat-functions for k=1, which we consider from here on), with N:=dim(Vh). As usual, each of these global basis functions, if restricted to a grid element, is given by the concatenation of a local shape function and the reference map: given

  • the invertible affine reference map FK:K^K for all KTh, such that x^x:=FK(x^), where K^ is the reference element associated with K,

  • a set of Lagrangian shape functions {φ^1,,φ^d+1}P1(K), each associated with a vertex a^1,,a^d+1 of the reference element K^ and

  • a DoF mapping σK:{1,,d+1}{1,,N} for all KTh, which associates the local index i^{1,,d+1} of a vertex a^i^ of the reference element K^ with the global index i{1,,N} of the vertex ai:=FK(a^i^) in the grid Th.

The DoF mapping as well as a localizable global basis is provided by a discrete function space in dune-gdt. We thus have

(5)φi|K=φ^i^FK1and(6)(φi)|K=(φ^i^FK1)=FK1(φ^i^FK1),

owing to the chain rule, with i:=σK(i^) for all 1i^d+1 and all KTh.

To obtain the algebraic analogue to (4), we

  • replace the bilinear form and functional by discrete counterparts acting on Vh, namely ah:Vh×VhR and lhVh (the construction of which is detailed further below) and

  • assemble the respective basis representations of ah and lh w.r.t. the basis of Vh into a matrix ahRN×N and vector lhRN, given by (ah)i,j:=ah(φj,φi) and (lh)i:=lh(φi), respectively, for 1i,jN;

  • and obtain the restrictions of ah and lh to VhH01(Ω) by modifying all entries associated with basis functions defined on the Dirichlet boundary: for each index i{1,N}, where the Lagrange-point defining the basis function φi lies on the Dirichlet boundary Ω, we set the ith row of ah to a unit row and clear the ith entry of lh.

The algebraic version of (4) then reads: find the vector of degrees of freedom (DoF) uhRN, such that

(7)ahuh=lh.

After solving the above linear system, we recover the solution of (4) from its basis representation uh=i=1Nuhiφi.

We consider for example a structured simplicial grid with 16 triangles.

from dune.xt.grid import Simplex, make_cube_grid, AllDirichletBoundaryInfo, visualize_grid

grid = make_cube_grid(Dim(d), Simplex(),
                      lower_left=omega[0], upper_right=omega[1], num_elements=[2, 2])
grid.global_refine(1)  # we need to refine once to obtain a symmetric grid

print(f'grid has {grid.size(0)} elements, '
      f'{grid.size(d - 1)} edges and {grid.size(d)} vertices')

boundary_info = AllDirichletBoundaryInfo(grid)

_ = visualize_grid(grid)
grid has 16 elements, 28 edges and 13 vertices
GridParameterBlock: Parameter 'refinementedge' not specified, defaulting to 'ARBITRARY'.
from dune.gdt import ContinuousLagrangeSpace

V_h = ContinuousLagrangeSpace(grid, order=1)

print(f'V_h has {V_h.num_DoFs} DoFs')

assert V_h.num_DoFs == grid.size(d)
V_h has 13 DoFs

approximate functionals#

Since the application of the functional to a global basis function ψi is localizable w.r.t. the grid, e.g.

(8)l(ψi)=KThKfψidx=:lK(ψi),

we first consider local functionals (such as lKL2(K)), where local means: with respect to a grid element K. Using the reference map FK and (5) from above, we transform the evaluation of lK(ψi) to the reference element,

lK(ψi)=Kfψidx=K^|detFK|(fFK)=:fK(ψ^i^FK1FK)dx^(9)=K^|detFK|fKψ^i^dx^,

where fK:K^R is the local function associated with f, i=σK(i^) and ψ^i^ is the corresponding shape function.

Note that, apart from the integration domain (K^ instead of K) and the transformation factor (|detFK|), the structure of the local functional from (8) is reflected in (9).

This leads us to the definition of a local functional in dune-gdt: ignoring the user input (namely the data function f for a moment), a local functional is determined by

  • an integrand, depending on a single test function, that we can evaluate at points on the reference element. We call such integrands unary element integrands. In the above example, given a test basis function ψ^ and a point in the reference element x^, the integrand is determined by Ξprod1,K:Pk(K^)×K^R, ψ^,x^fK(x^)ψ^(x^), which is modelled by LocalElementProductIntegrand in dune-gdt (see below); and

  • an approximation of the integral in (9) by a numerical quadrature: given any unary element integrand Ξ1,K, and QN quadrature points x^1,,x^Q and weights ω1,,ωQR, we approximate lhK(ψi):=q=1Q|detFK(x^q)|ωqΞ1,K(ψ^i^,x^q)K^Ξ1,K(ψ^i^,x^)dx^=lK(ψi), which is modelled by LocalElementIntegralFunctional in dune-gdt (see below).

    Note that the order of the quadrature is determined automatically, since the integrand computes its polynomial degree given all data functions and basis functions (in the above example, the polynomial order of fK is 3 by our construction and the polynomial order of ψ^ is 1, since we are using piecewise linear shape functions, yielding a polynomial order of 4 for Ξprod1,K).

Given local functionals, the purpose of the VectorFunctional in dune-gdt is to assemble lh from (7) by

  • creating an appropriate vector of length N

  • iterating over all grid elements KTh

  • localizing the basis of Vh w.r.t. each grid element K

  • evaluating the local functionals lhK for each localized basis function

  • adding the results to the respective entry of lh, determined by the DoF-mapping of the discrete function space ContinuousLagrangeSpace

In our example, we define lh as:

from dune.xt.functions import GridFunction as GF

from dune.gdt import (
    VectorFunctional,
    LocalElementProductIntegrand,
    LocalElementIntegralFunctional,
)

l_h = VectorFunctional(grid, source_space=V_h)
l_h += LocalElementIntegralFunctional(
         LocalElementProductIntegrand(GF(grid, 1)).with_ansatz(GF(grid, f)))

A few notes regarding the above code:

  • there exists a large variety of data functions, but in order to all handle them in dune-gdt we require them to be localizable w.r.t. a grid (i.e. to have local functions as above). This is achieved by wrapping them into a GridFunction, which accepts all kind of functions, discrete functions or numbers. Thus GF(grid, 1) creates a grid function which is localizable w.r.t. each grid elements and evaluates to 1, whenever evaluated; whereas GF(grid, f), when localized to a grid element K and evaluated at a point on the associated reference element, x^, evaluates to f(FK(x^)). See also the tutorial on data functions.

  • the LocalElementProductIntegrand is actually a binary element integrand modelling a weighted product, as in: with a weight function w:ΩR, given an ansatz function φ^, a test function ψ^ and a point x^K^, this integrand is determined by Ξprod2,K:Pk(K^)×Pk(K^)×K^R, φ^,ψ^,x^wK(x^)φ^(x^)ψ^(x^). Thus, LocalElementProductIntegrand is often used in bilinear forms to assemble L2 products (with weight w=1), which we achieve by LocalElementProductIntegrand(GF(grid, 1)). However, even with w=1, the integrand Ξprod2,K still depends on the test and ansatz function. Using with_ansatz(GF(grid, f)), we fix fK as the ansatz function to obtain exactly the unary integrand we require, which only depends on the test function, Ξprod1,K:Pk(K^)×K^R, ψ^,x^Ξprod2,K(fK,ψ^,x^)=fK(x^)ψ^(x^), which is exactly what we need to approximate lsrcK(ψi)=Kfψidx.

  • the above code creates the vector lh (available as the vector attribute of l_h) filled with 0, but does not yet assemble the functional into it, which we can check by:

assert len(l_h.vector) == V_h.num_DoFs

print(l_h.vector.sup_norm())
0.0

approximate bilinear forms#

The approximation of the application of the bilinear form a to two global basis function ψi,φj follows in a similar manner. We obtain by localization

a(ψi,φj)=Ω(κφj)ψidx=KThK(κφj)ψidx=:aK(ψi,φj)

and by transformation and the chain rule, using (6) and FK1FK=id

aK(ψi,φj)=K^|detFK|((κFK)=:κK(FK1φ^j^)=:Kφ^j^)(FK1ψ^i^)=:Kψ^i^dx^=K^|detFK|(κKKφ^j^)Kψ^i^dx^,

where κK denote the local function of κ as above, and where Kφ^j^ denote suitably transformed global gradients of the local shape functions, for the integrand to reflect the same structure as above.

Similar to local functionals, a local bilinear form is determined

  • by a binary element integrand, in our case the LocalLaplaceIntegrand (see below) $Ξlaplace2,K:Pk(K^)×Pk(K^)×K^Rφ^,ξ^,x^(κK(x^)Kφ^(x^))Kψ^(x^)$ and

  • an approximation of the integral by a numerical quadrature: given any binary element integrand Ξ2,K, and QN quadrature points x^1,,x^Q and weights ω1,,ωQR, we approximate ahK(ψi,φj):=q=1Q|detFK(x^q)|ωqΞ2,K(ψ^i^,φ^j^,x^q)K^Ξ2,K(ψ^i^,φ^j^,x^)dx^=aK(ψi,φi), which is modelled by LocalElementIntegralBilinearForm in dune-gdt (see below).

Given local bilinear forms, the purpose of the MatrixOperator in dune-gdt is to assemble ah from (7) by

  • creating an appropriate (sparse) matrix of size N×N

  • iterating over all grid elements KTh

  • localizing the basis of Vh w.r.t. each grid element K

  • evaluating the local bilinear form ahK for each combination localized ansatz and test basis functions

  • adding the results to the respective entry of ah, determined by the DoF-mapping of the discrete function space ContinuousLagrangeSpace

from dune.gdt import (
    MatrixOperator,
    make_element_sparsity_pattern,
    LocalLaplaceIntegrand,
    LocalElementIntegralBilinearForm,
)

a_h = MatrixOperator(grid, source_space=V_h, range_space=V_h,
                     sparsity_pattern=make_element_sparsity_pattern(V_h))
a_h += LocalElementIntegralBilinearForm(
        LocalLaplaceIntegrand(GF(grid, kappa, dim_range=(Dim(d), Dim(d)))))

A few notes regarding the above code:

  • the LocalLaplaceIntegrand expects a matrix-valued function, which we achieve by converting the scalar function kappa to a matrix-valued GridFunction

  • the above code creates the matrix ah (available as the matrix attribute of a_h) with sparse 0 entries, but does not yet assemble the bilinear form into it, which we can check by:

assert a_h.matrix.rows == a_h.matrix.cols == V_h.num_DoFs

print(a_h.matrix.sup_norm())
0.0

handling the Dirichlet boundary condition#

As noted above, we handle the Dirichlet boundary condition on the algebraic level by modifying the assembled matrices and vector. We therefore require a means to identify all DoFs of Vh associated with the Dirichlet boundary modelled by boundary_info.

from dune.gdt import DirichletConstraints

dirichlet_constraints = DirichletConstraints(boundary_info, V_h)

Similar to the bilinear forms and functionals above, the dirichlet_constraints are not yet assembled, which we can check as follows:

print(dirichlet_constraints.dirichlet_DoFs)
set()

walking the grid#

Until now, we constructed a bilinear form a_h, a linear functional l_h and Dirichlet constrinaints dirichlet_constraints, which are all localizable w.r.t. the grid, that is: i order to compute their application or to assemble them, it is sufficient to apply them to each element of the grid.

Internally, this is realized by means of the Walker from dune-xt-grid, which allows to register all kinds of grid functors which are applied locally on each element. All bilinear forms, operators and functionals (as well as other constructs such as the Dirichlet constraints) in dune-gdt are implemented as such functors.

Thus, we may assemble everything in one grid walk:

from dune.xt.grid import Walker

walker = Walker(grid)
walker.append(a_h)
walker.append(l_h)
walker.append(dirichlet_constraints)
walker.walk()

We can check that the assembled bilinear form and functional as well as the Dirichlet constraints actually contain some data:

print(f'a_h = {a_h.matrix.__repr__()}')
print()
print(f'l_h = {l_h.vector.__repr__()}')
print()
print(f'Dirichlet DoFs: {dirichlet_constraints.dirichlet_DoFs}')
a_h = IstlSparseMatrix(
  [4 0 -1 ... 0 0 -1]
  [0 2 -1 ... 0 0 0]
  [-1 -1 4 ... 0 0 0]
   ...
  [0 0 0 ... 2 0 -1]
  [0 0 0 ... 0 1 -1]
  [-1 0 0 ... -1 -1 4]
)

l_h = IstlVector([0.215393 0.088827 0.0887864 ... 0.130627 0.0897916 0.147315])

Dirichlet DoFs: {1, 3, 4, 5, 6, 8, 10, 11}

solving the linear system#

After walking the grid, the bilinra form and linear functional are assembled w.r.t. Vh and we constrain them to include the handling of the Dirichlet boundary condition.

dirichlet_constraints.apply(a_h.matrix, l_h.vector)

Since the bilinear form is implemented as a MatrixOperator, we may simply invert the operator to obtain the DoF vector of the solution of (4).

u_h_vector = a_h.apply_inverse(l_h.vector)

postprocessing the solution#

To make use of the DoF vector of the approximate solution, uhRN, it is convenient to interpert it as a discrete function again, uhVh by means of the Galerkin isomorphism. This can be achieved by the DiscreteFunction in dune-gdt.

All discrete functions are in particular grid functions and can thus be compared to analytical solutions, used as input in discretization schemes or visualized.

Note: if visualization fails for some reason, call paraview on the command line and open u_h.vtu!

from dune.gdt import DiscreteFunction, visualize_function

u_h = DiscreteFunction(V_h, u_h_vector, name='u_h')
_ = visualize_function(u_h)

everything in a single function#

For a better overview, the above discretization code is also available in a single function in the file discretize_elliptic_cg.py:

import inspect
from discretize_elliptic_cg import discretize_elliptic_cg_dirichlet_zero

print(inspect.getsource(discretize_elliptic_cg_dirichlet_zero))
def discretize_elliptic_cg_dirichlet_zero(grid, diffusion, source):

    """
    Discretizes the stationary heat equation with homogeneous Dirichlet boundary values everywhere
    with dune-gdt using continuous Lagrange finite elements.

    Parameters
    ----------
    grid
        The grid, given as a GridProvider from dune.xt.grid.
    diffusion
        Diffusion function with values in R^{d x d}, anything that dune.xt.functions.GridFunction
        can handle.
    source
        Right hand side source function with values in R, anything that
        dune.xt.functions.GridFunction can handle.

    Returns
    -------
    u_h
        The computed solution as a dune.gdt.DiscreteFunction for postprocessing.
    """

    d = grid.dimension
    diffusion = GF(grid, diffusion, dim_range=(Dim(d), Dim(d)))
    source = GF(grid, source)

    boundary_info = AllDirichletBoundaryInfo(grid)

    V_h = ContinuousLagrangeSpace(grid, order=1)

    l_h = VectorFunctional(grid, source_space=V_h)
    l_h += LocalElementIntegralFunctional(
        LocalElementProductIntegrand(GF(grid, 1)).with_ansatz(source)
    )

    a_h = MatrixOperator(
        grid,
        source_space=V_h,
        range_space=V_h,
        sparsity_pattern=make_element_sparsity_pattern(V_h),
    )
    a_h += LocalElementIntegralBilinearForm(LocalLaplaceIntegrand(diffusion))

    dirichlet_constraints = DirichletConstraints(boundary_info, V_h)

    walker = Walker(grid)
    walker.append(a_h)
    walker.append(l_h)
    walker.append(dirichlet_constraints)
    walker.walk()

    u_h = DiscreteFunction(V_h, name="u_h")
    dirichlet_constraints.apply(a_h.matrix, l_h.vector)
    a_h.apply_inverse(l_h.vector, u_h.dofs.vector)

    return u_h

Calling it gives the same solution as above:

u_h = discretize_elliptic_cg_dirichlet_zero(grid, kappa, f)
_ = visualize_function(u_h)

diffusion with non-homogeneous Dirichlet boundary condition#

analytical problem#

Consider problem (1) from above, but with non-homogeneous Dirichlet boundary values. That is:

Let ΩRd for 1d3 be a bounded connected domain with Lipschitz-boundary Ω. We seek the solution uH1(Ω) of the linear diffusion equation (with a non-homogeneous Dirichlet boundary condition)

(10)(κu)=fin Ω,u=gDon Ω,

in a weak sense, where κ[L(Ω)]d×d denotes a given diffusion function, fL2(Ω) denotes a given source function and gDL2(Ω) denotes given Dirichlet boundary values.

The variational problem associated with (10) reads: find uH1(Ω), such that

(11)a(u,v)=l(v)for all vH01(Ω),

with the same bilinear form a and linear functional l as above in (3).

Dirichlet shift#

Suppose g^DH1(Ω) is an extension of the Dirichlet boundary values to Ω, such that

g^D|Ω=gD

in the sense of traces. We then have for the solution uH1(Ω) of (11), that

u=u0+g^D

for some u0H01(Ω). Thus, we may reformulate (11) as follows: Find u0H01(Ω), such that

a(u0+g^D,v)=l(v)for all vH01(Ω),

or equivalently

a(u0,v)=l(v)a(g^D,v)for all vH01(Ω).

We have thus shifted problem (11) to be of familiar form, i.e., similarly to above we consider a linear diffusion equation with homogeneous Dirichlet boundary conditions, but a modified source term.

Consider for example (1) with:

  • d=2

  • Ω=[0,1]2

  • κ=1

  • f=exp(x0x1)

  • gD=x0x1

Except for gD, we may thus use the same grid, boundary info and data functions as above.

Dirichlet interpolation#

To obtain g^DH1(Ω), we use the Lagrange interpolation of gD and set all DoFs associated with inner Lagrange nodes to zero. This is realized by the boundary_interpolation.

from dune.xt.grid import DirichletBoundary
from dune.gdt import boundary_interpolation

g_D = ExpressionFunction(dim_domain=Dim(d),
                         variable='x', expression='x[0]*x[1]', order=2, name='g_D')

g_D_hat = boundary_interpolation(GF(grid, g_D), V_h, boundary_info, DirichletBoundary())

_ = visualize_function(g_D_hat)

As we observe, the values on all boundary DoFs are x0x1 and on all inner DoFs 0.

assembling the shifted system#

  • We assemble the unconstrained matrix- and vector representation of ah and lh w.r.t. Vh similarly as above.

l_h = VectorFunctional(grid, source_space=V_h)
l_h += LocalElementIntegralFunctional(
        LocalElementProductIntegrand(GF(grid, 1)).with_ansatz(GF(grid, f)))

a_h = MatrixOperator(grid, source_space=V_h, range_space=V_h,
                     sparsity_pattern=make_element_sparsity_pattern(V_h))
a_h += LocalElementIntegralBilinearForm(
        LocalLaplaceIntegrand(GF(grid, kappa, dim_range=(Dim(d), Dim(d)))))

walker = Walker(grid)
walker.append(a_h)
walker.append(l_h)
walker.walk()
  • We then obtain the shifted system by directly modifying the right hand side vector.

rhs_vector = l_h.vector - a_h.matrix@g_D_hat.dofs.vector
  • Afterwards, we constrain the shifted system to VhH01(Ω).

dirichlet_constraints.apply(a_h.matrix, rhs_vector)
  • Then, we solve the shifted and constrained system for the DoF vector of u0,h.

u_0_h_vector = a_h.apply_inverse(rhs_vector)
  • We may then obtain the solution by shifting it back.

u_h = DiscreteFunction(V_h, u_0_h_vector + g_D_hat.dofs.vector)

_ = visualize_function(u_h)

everything in a single function#

As above, solving with non-homogeneous Dirichlet values is also available in a single function.

from discretize_elliptic_cg import discretize_elliptic_cg_dirichlet

u_h = discretize_elliptic_cg_dirichlet(grid, kappa, f, g_D)

_ = visualize_function(u_h)

Download the code: dune_gdt_tutorial_on_cg_fem_for_the_stationary_heat_equation.md dune_gdt_tutorial_on_cg_fem_for_the_stationary_heat_equation.ipynb