Run this tutorial

Click here to run this tutorial 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 $\mathrm{\Omega }\subset {\mathbb{R}}^d$ for $1\le d\le 3$ be a bounded connected domain with Lipschitz-boundary $\partial \mathrm{\Omega }$. We seek the solution $u\in H_0^1(\mathrm{\Omega })$ of the linear diffusion equation (with a homogeneous Dirichlet boundary condition)

$$\begin{array}{r}\begin{array}{rlrl}\text{(1)}& -\mathrm{\nabla }\cdot (\kappa \mathrm{\nabla }u)& =f& & \text{in }\mathrm{\Omega },u& =0& & \text{on }\partial \mathrm{\Omega },\end{array}\end{array}$$

in a weak sense, where $\kappa \in [L^{\mathrm{\infty }}(\mathrm{\Omega }){]}^{d\times d}$ denotes a given diffusion function and $f\in L^2(\mathrm{\Omega })$ denotes a given source function.

The variational problem associated with $\text{(1)}$ reads: find $u\in H_0^1(\mathrm{\Omega })$, such that

$$\begin{array}{}\text{(2)}& a(u,v)& =l(v)& & \text{for all }v\in H_0^1(\mathrm{\Omega }),\end{array}$$

where the bilinear form $a:H^1(\mathrm{\Omega })\times H^1(\mathrm{\Omega })\to \mathbb{R}$ and the linear functional $l\in H^{-1}(\mathrm{\Omega })$ are given by

$$\begin{array}{}\text{(3)}& a(u,v):={\int }_{\mathrm{\Omega }}(\kappa \mathrm{\nabla }u)\cdot v\text{d}x& & \text{and}& & l(v):={\int }_{\mathrm{\Omega }}f\text{d}x,\end{array}$$

respectively.

Consider for example $\text{(1)}$ with:

• $d=2$

• $\mathrm{\Omega }=[0,1{]}^2$

• $\kappa =1$

• $f=\exp (x_0{x}_1)$

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 ${\mathcal{T}}_h$ (other types of elements work analogously) as a partition of $\mathrm{\Omega }$ into elements $K\in {\mathcal{T}}_h$ with $h:=\underset{K\in {\mathcal{T}}_h}{max}\text{diam}(K)$, we consider the discrete space of continuous piecewise polynomial functions of order $k\in \mathbb{N}$,

$$\begin{array}{r}{V}_h:=\{v\in C^0(\mathrm{\Omega })|v{|}_K\in {\mathbb{P}}^k(K)\}\end{array}$$

where ${\mathbb{P}}^k(K)$ denotes the space of polynomials of (total) degree up to $k$ (note that $V_h\subset H^1(\mathrm{\Omega })$ and $V_h$ does not include the Dirichlet boundary condition, thus $V_h\not\subset H_0^1(\mathrm{\Omega }$.). We obtain a finite-dimensional variational problem by Galerkin-projection of $\text{(2)}$ onto $V_h$, thas is: we seek the approximate solution $u_h\in V_h\cap H_0^1(\mathrm{\Omega })$, such that

$$\begin{array}{}\text{(4)}& a(u_h,v_h)& =l(v_h)& & \text{for all }{v}_h\in V_h\cap H_0^1(\mathrm{\Omega }).\end{array}$$

A basis of $V_h$ is given by the Lagrangian basis functions

$$\begin{array}{r}\Phi :=\{{\phi }_1,\dots ,{\phi }_N\}\end{array}$$

of order $k$ (e.g., the usual hat-functions for $k=1$, which we consider from here on), with $N:=\text{dim}(V_h)$. 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 $F_K:\hat{K}\to K$ for all $K\in {\mathcal{T}}_h$, such that $\hat{x}\mapsto x:=F_K(\hat{x})$, where $\hat{K}$ is the reference element associated with $K$,

• a set of Lagrangian shape functions $\{{\hat{\phi }}_1,\dots ,{\hat{\phi }}_{d+1}\}\in {\mathbb{P}}^1(K)$, each associated with a vertex ${\hat{a}}_1,\dots ,{\hat{a}}_{d+1}$ of the reference element $\hat{K}$ and

• a DoF mapping ${\sigma }_K:\{1,\dots ,d+1\}\to \{1,\dots ,N\}$ for all $K\in {\mathcal{T}}_h$, which associates the local index $\hat{i}\in \{1,\dots ,d+1\}$ of a vertex ${\hat{a}}_{\hat{i}}$ of the reference element $\hat{K}$ with the global index $i\in \{1,\dots ,N\}$ of the vertex $a_i:=F_K({\hat{a}}_{\hat{i}})$ in the grid ${\mathcal{T}}_h$.

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

$$\begin{array}{r}\begin{array}{}\text{(5)}& {\phi }_i{|}_K& ={\hat{\phi }}_{\hat{i}}\circ F_K^{-1}& & \text{and}\text{(6)}& (\mathrm{\nabla }{\phi }_i){|}_K& =\mathrm{\nabla }({\hat{\phi }}_{\hat{i}}\circ F_K^{-1})=\mathrm{\nabla }{F}_K^{-1}\cdot (\mathrm{\nabla }{\hat{\phi }}_{\hat{i}}\circ F_K^{-1}),\end{array}\end{array}$$

owing to the chain rule, with $i:={\sigma }_K(\hat{i})$ for all $1\le \hat{i}\le d+1$ and all $K\in {\mathcal{T}}_h$.

To obtain the algebraic analogue to $\text{(4)}$, we

• replace the bilinear form and functional by discrete counterparts acting on $V_h$, namely $a_h:V_h\times V_h\to \mathbb{R}$ and $l_h\in V_h^{\prime }$ (the construction of which is detailed further below) and

• assemble the respective basis representations of $a_h$ and $l_h$ w.r.t. the basis of $V_h$ into a matrix $\underset{―}{a_h}\in {\mathbb{R}}^{N\times N}$ and vector $\underset{―}{l_h}\in {\mathbb{R}}^N$, given by $(\underset{―}{a_h}{)}_{i,j}:=a_h({\phi }_j,{\phi }_i)$ and $(\underset{―}{l_h}{)}_i:=l_h({\phi }_i)$, respectively, for $1\le i,j\le N$;

• and obtain the restrictions of $\underset{―}{a_h}$ and $\underset{―}{l_h}$ to $V_h\cap H_0^1(\mathrm{\Omega })$ by modifying all entries associated with basis functions defined on the Dirichlet boundary: for each index $i\in \{1,\dots N\}$, where the Lagrange-point defining the basis function ${\phi }_i$ lies on the Dirichlet boundary $\partial \mathrm{\Omega }$, we set the $i$th row of $\underset{―}{a_h}$ to a unit row and clear the $i$th entry of $\underset{―}{l_h}$.

The algebraic version of $\text{(4)}$ then reads: find the vector of degrees of freedom (DoF) $\underset{―}{u_h}\in {\mathbb{R}}^N$, such that

$$\begin{array}{}\text{(7)}& \underset{―}{a_h}\underset{―}{u_h}=\underset{―}{l_h}.\end{array}$$

After solving the above linear system, we recover the solution of $\text{(4)}$ from its basis representation $u_h=\sum _{i=1}^N{\underset{―}{u_h}}_i{\phi }_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 ${\psi }_i$ is localizable w.r.t. the grid, e.g.

$$\begin{array}{}\text{(8)}& l({\psi }_i)=\sum _{K\in {\mathcal{T}}_h}\underset{=:l^K({\psi }_i)}{\underbrace{{\int }_{K}f{\psi }_i\text{d}x}},\end{array}$$

we first consider local functionals (such as $l^K\in L^2(K{)}^{\prime }$), where local means: with respect to a grid element $K$. Using the reference map $F_K$ and $\text{(5)}$ from above, we transform the evaluation of $l^K({\psi }_i)$ to the reference element,

$$\begin{array}{r}\begin{array}{rl}{l}^K({\psi }_i)& ={\int }_{K}f{\psi }_i\text{d}x={\int }_{\hat{K}}|\text{det}\mathrm{\nabla }{F}_K|\underset{=:f^K}{\underbrace{(f\circ F_K)}}({\hat{\psi }}_{\hat{i}}\circ F_K^{-1}\circ F_K)\text{d}\hat{x}\\ \text{(9)}& & ={\int }_{\hat{K}}|\text{det}\mathrm{\nabla }{F}_K|f^K{\hat{\psi }}_{\hat{i}}\text{d}\hat{x},\end{array}\end{array}$$

where $f^K:\hat{K}\to \mathbb{R}$ is the local function associated with $f$, $i={\sigma }_K(\hat{i})$ and ${\hat{\psi }}_{\hat{i}}$ is the corresponding shape function.

Note that, apart from the integration domain ($\hat{K}$ instead of $K$) and the transformation factor ($|\text{det}\mathrm{\nabla }{F}^K|$), the structure of the local functional from $\text{(8)}$ is reflected in $\text{(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 $\hat{\psi }$ and a point in the reference element $\hat{x}$, the integrand is determined by ${\mathrm{\Xi }}_{\text{prod}}^{1,K}:{\mathbb{P}}^k(\hat{K})\times \hat{K}\to \mathbb{R}$, $\hat{\psi },\hat{x}\mapsto f^K(\hat{x})\hat{\psi }(\hat{x})$, which is modelled by LocalElementProductIntegrand in dune-gdt (see below); and

• an approximation of the integral in $\text{(9)}$ by a numerical quadrature: given any unary element integrand ${\mathrm{\Xi }}^{1,K}$, and $Q\in \mathbb{N}$ quadrature points ${\hat{x}}_1,\dots ,{\hat{x}}_Q$ and weights ${\omega }_1,\dots ,{\omega }_Q\in \mathbb{R}$, we approximate $l_h^K({\psi }_i):=\sum _{q=1}^Q|\text{det}\mathrm{\nabla }{F}_K({\hat{x}}_q)|{\omega }_q{\mathrm{\Xi }}^{1,K}({\hat{\psi }}_{\hat{i}},{\hat{x}}_q)\approx {\int }_{\hat{K}}{\mathrm{\Xi }}^{1,K}({\hat{\psi }}_{\hat{i}},\hat{x})\text{d}\hat{x}=l^K({\psi }_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 $f^K$ is 3 by our construction and the polynomial order of $\hat{\psi }$ is 1, since we are using piecewise linear shape functions, yielding a polynomial order of 4 for ${\mathrm{\Xi }}_{\text{prod}}^{1,K}$).

Given local functionals, the purpose of the VectorFunctional in dune-gdt is to assemble $\underset{―}{l_h}$ from $\text{(7)}$ by

• creating an appropriate vector of length $N$

• iterating over all grid elements $K\in {\mathcal{T}}_h$

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

• evaluating the local functionals $l_h^K$ for each localized basis function

• adding the results to the respective entry of $\underset{―}{l_h}$, determined by the DoF-mapping of the discrete function space ContinuousLagrangeSpace

In our example, we define $l_h$ 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, $\hat{x}$, evaluates to $f(F_K(\hat{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:\mathrm{\Omega }\to \mathbb{R}$, given an ansatz function $\hat{\phi }$, a test function $\hat{\psi }$ and a point $\hat{x}\in \hat{K}$, this integrand is determined by ${\mathrm{\Xi }}_{\text{prod}}^{2,K}:{\mathbb{P}}^k(\hat{K})\times {\mathbb{P}}^k(\hat{K})\times \hat{K}\to \mathbb{R}$, $\hat{\phi },\hat{\psi },\hat{x}\mapsto w^K(\hat{x})\hat{\phi }(\hat{x})\hat{\psi }(\hat{x})$. Thus, LocalElementProductIntegrand is often used in bilinear forms to assemble $L^2$ products (with weight $w=1$), which we achieve by LocalElementProductIntegrand(GF(grid, 1)). However, even with $w=1$, the integrand ${\mathrm{\Xi }}_{\text{prod}}^{2,K}$ still depends on the test and ansatz function. Using with_ansatz(GF(grid, f)), we fix $f^K$ as the ansatz function to obtain exactly the unary integrand we require, which only depends on the test function, ${\mathrm{\Xi }}_{\text{prod}}^{1,K}:{\mathbb{P}}^k(\hat{K})\times \hat{K}\to \mathbb{R}$, $\hat{\psi },\hat{x}\mapsto {\mathrm{\Xi }}_{\text{prod}}^{2,K}(f^K,\hat{\psi },\hat{x})=f^K(\hat{x})\hat{\psi }(\hat{x})$, which is exactly what we need to approximate $l_{\text{src}}^K({\psi }_i)={\int }_{K}f{\psi }_i\text{d}x$.

• the above code creates the vector $\underset{―}{l_h}$ (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 ${\psi }_i,{\phi }_j$ follows in a similar manner. We obtain by localization

$$\begin{array}{rl}a({\psi }_i,{\phi }_j)& ={\int }_{\mathrm{\Omega }}(\kappa \mathrm{\nabla }{\phi }_j)\cdot \mathrm{\nabla }{\psi }_i\text{d}x=\sum _{K\in {\mathcal{T}}_h}\underset{=:a^K({\psi }_i,{\phi }_j)}{\underbrace{{\int }_K(\kappa \mathrm{\nabla }{\phi }_j)\cdot \mathrm{\nabla }{\psi }_i\text{d}x}}\end{array}$$

and by transformation and the chain rule, using $\text{(6)}$ and $F_K^{-1}\circ F_K=\text{id}$

$$\begin{array}{r}\begin{array}{rl}{a}^K({\psi }_i,{\phi }_j)& ={\int }_{\hat{K}}|\text{det}\mathrm{\nabla }{F}_K|(\underset{=:{\kappa }^K}{\underbrace{(\kappa \circ F_K)}}\underset{=:{\mathrm{\nabla }}_K{\hat{\phi }}_{\hat{j}}}{\underbrace{(\mathrm{\nabla }{F}_K^{-1}\cdot \mathrm{\nabla }{\hat{\phi }}_{\hat{j}})}})\cdot \underset{=:{\mathrm{\nabla }}_K{\hat{\psi }}_{\hat{i}}}{\underbrace{(\mathrm{\nabla }{F}_K^{-1}\cdot \mathrm{\nabla }{\hat{\psi }}_{\hat{i}})}}\text{d}\hat{x}\\ & ={\int }_{\hat{K}}|\text{det}\mathrm{\nabla }{F}_K|({\kappa }^K{\mathrm{\nabla }}_K{\hat{\phi }}_{\hat{j}})\cdot {\mathrm{\nabla }}_K{\hat{\psi }}_{\hat{i}}\text{d}\hat{x},\end{array}\end{array}$$

where ${\kappa }^K$ denote the local function of $\kappa$ as above, and where ${\mathrm{\nabla }}_K{\hat{\phi }}_{\hat{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) $$\begin{array}{rl}{\mathrm{\Xi }}_{\text{laplace}}^{2,K}:{\mathbb{P}}^k(\hat{K})\times {\mathbb{P}}^k(\hat{K})\times \hat{K}& \to \mathbb{R}\\ \hat{\phi },\hat{\xi },\hat{x}& \mapsto ({\kappa }^K(\hat{x}){\mathrm{\nabla }}_K\hat{\phi }(\hat{x}))\cdot {\mathrm{\nabla }}_K\hat{\psi }(\hat{x})\end{array}$$ and

• an approximation of the integral by a numerical quadrature: given any binary element integrand ${\mathrm{\Xi }}^{2,K}$, and $Q\in \mathbb{N}$ quadrature points ${\hat{x}}_1,\dots ,{\hat{x}}_Q$ and weights ${\omega }_1,\dots ,{\omega }_Q\in \mathbb{R}$, we approximate $\begin{array}{r}{a}_h^K({\psi }_i,{\phi }_j):=\sum _{q=1}^Q|\text{det}\mathrm{\nabla }{F}_K({\hat{x}}_q)|{\omega }_q{\mathrm{\Xi }}^{2,K}({\hat{\psi }}_{\hat{i}},{\hat{\phi }}_{\hat{j}},{\hat{x}}_q)\approx {\int }_{\hat{K}}{\mathrm{\Xi }}^{2,K}({\hat{\psi }}_{\hat{i}},{\hat{\phi }}_{\hat{j}},\hat{x})\text{d}\hat{x}=a^K({\psi }_i,{\phi }_i),\end{array}$ 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 $\underset{―}{a_h}$ from $\text{(7)}$ by

• creating an appropriate (sparse) matrix of size $N\times N$

• iterating over all grid elements $K\in {\mathcal{T}}_h$

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

• evaluating the local bilinear form $a_h^K$ for each combination localized ansatz and test basis functions

• adding the results to the respective entry of $\underset{―}{a_h}$, 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 $\underset{―}{a_h}$ (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 $V_h$ 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. $V_h$ 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 $\text{(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, $\underset{―}{u_h}\in {\mathbb{R}}^N$, it is convenient to interpert it as a discrete function again, $u_h\in V_h$ 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 $\text{(1)}$ from above, but with non-homogeneous Dirichlet boundary values. That is:

Let $\mathrm{\Omega }\subset {\mathbb{R}}^d$ for $1\le d\le 3$ be a bounded connected domain with Lipschitz-boundary $\partial \mathrm{\Omega }$. We seek the solution $u\in H^1(\mathrm{\Omega })$ of the linear diffusion equation (with a non-homogeneous Dirichlet boundary condition)

$$\begin{array}{r}\begin{array}{rlrl}\text{(10)}& -\mathrm{\nabla }\cdot (\kappa \mathrm{\nabla }u)& =f& & \text{in }\mathrm{\Omega },u& =g_{\text{D}}& & \text{on }\partial \mathrm{\Omega },\end{array}\end{array}$$

in a weak sense, where $\kappa \in [L^{\mathrm{\infty }}(\mathrm{\Omega }){]}^{d\times d}$ denotes a given diffusion function, $f\in L^2(\mathrm{\Omega })$ denotes a given source function and $g_{\text{D}}\in L^2(\partial \mathrm{\Omega })$ denotes given Dirichlet boundary values.

The variational problem associated with $\text{(10)}$ reads: find $u\in H^1(\mathrm{\Omega })$, such that

$$\begin{array}{}\text{(11)}& a(u,v)& =l(v)& & \text{for all }v\in H_0^1(\mathrm{\Omega }),\end{array}$$

with the same bilinear form $a$ and linear functional $l$ as above in $\text{(3)}$.

Dirichlet shift

Suppose ${\hat{g}}_{\text{D}}\in H^1(\mathrm{\Omega })$ is an extension of the Dirichlet boundary values to $\mathrm{\Omega }$, such that

$$\begin{array}{r}{\hat{g}}_{\text{D}}{|}_{\partial \mathrm{\Omega }}=g_{\text{D}}\end{array}$$

in the sense of traces. We then have for the solution $u\in H^1(\mathrm{\Omega })$ of $\text{(11)}$, that

$$\begin{array}{r}u=u_0+{\hat{g}}_{\text{D}}\end{array}$$

for some $u_0\in H_0^1(\mathrm{\Omega })$. Thus, we may reformulate $\text{(11)}$ as follows: Find $u_0\in H_0^1(\mathrm{\Omega })$, such that

$$\begin{array}{rlrl}a(u_0+{\hat{g}}_{\text{D}},v)& =l(v)& & \text{for all }v\in H_0^1(\mathrm{\Omega }),\end{array}$$

or equivalently

$$\begin{array}{rlrl}a(u_0,v)& =l(v)-a({\hat{g}}_{\text{D}},v)& & \text{for all }v\in H_0^1(\mathrm{\Omega }).\end{array}$$

We have thus shifted problem $\text{(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 $\text{(1)}$ with:

• $d=2$

• $\mathrm{\Omega }=[0,1{]}^2$

• $\kappa =1$

• $f=\exp (x_0{x}_1)$

• $g_{\text{D}}=x_0{x}_1$

Except for $g_{\text{D}}$, we may thus use the same grid, boundary info and data functions as above.

Dirichlet interpolation

To obtain ${\hat{g}}_{\text{D}}\in H^1(\mathrm{\Omega })$, we use the Lagrange interpolation of $g_{\text{D}}$ 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 $x_0{x}_1$ and on all inner DoFs $0$.

assembling the shifted system

• We assemble the unconstrained matrix- and vector representation of $a_h$ and $l_h$ w.r.t. $V_h$ 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 $V_h\cap H_0^1(\mathrm{\Omega })$.

dirichlet_constraints.apply(a_h.matrix, rhs_vector)

• Then, we solve the shifted and constrained system for the DoF vector of $u_{0,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