Run this tutorial

Click here to run this tutorial on mybinder.org:# 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 \(\Omega \subset \mathbb{R}^d\) for \(1 \leq d \leq 3\) be a bounded connected domain with Lipschitz-boundary \(\partial\Omega\). We seek the solution \(u \in H^1_0(\Omega)\) of the **linear diffusion equation** (with a homogeneous Dirichlet boundary condition)

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

The variational problem associated with \(\eqref{pde}\) reads: find \(u \in H^1_0(\Omega)\), such that

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

respectively.

Consider for example \(\eqref{pde}\) with:

\(d = 2\)

\(\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 \(\Omega\) into elements \(K \in \mathcal{T}_h\) with \(h := \max_{K \in \mathcal{T}_h} \text{diam}(K)\), we consider the discrete space of continuous piecewise polynomial functions of order \(k \in \mathbb{N}\),

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

A basis of \(V_h\) is given by the Lagrangian basis functions

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{\varphi}_1, \dots, \hat{\varphi}_{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}\) anda

**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

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

To obtain the algebraic analogue to \(\eqref{discrete_variational_problem}\), 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'\) (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 \(\underline{a_h} \in \mathbb{R}^{N \times N}\) and vector \(\underline{l_h} \in \mathbb{R}^N\), given by \((\underline{a_h})_{i, j} := a_h(\varphi_j, \varphi_i)\) and \((\underline{l_h})_i := l_h(\varphi_i)\), respectively, for \(1 \leq i, j \leq N\);

and obtain the restrictions of \(\underline{a_h}\) and \(\underline{l_h}\) to \(V_h \cap H^1_0(\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 \(\varphi_i\) lies on the Dirichlet boundary \(\partial\Omega\), we set the \(i\)th row of \(\underline{a_h}\) to a unit row and clear the \(i\)th entry of \(\underline{l_h}\).

The algebraic version of \(\eqref{discrete_variational_problem}\) then reads: find the vector of degrees of freedom (DoF) \(\underline{u_h} \in \mathbb{R}^N\), such that

After solving the above linear system, we recover the solution of \(\eqref{discrete_variational_problem}\) from its basis representation \(u_h = \sum_{i=1}^{N}\underline{u_h}_i \varphi_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.

we first consider local functionals (such as \(l^K \in L^2(K)'\)), where *local* means: *with respect to a grid element \(K\)*. Using the reference map \(F_K\) and \(\eqref{basis_transformation}\) from above, we transform the evaluation of \(l^K(\psi_i)\) to the reference element,

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}\nabla F^K|\)), the structure of the local functional from \(\eqref{localized_rhs}\) is reflected in \(\eqref{transformed_localized_rhs}\).

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 integrand**s. 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 \(\Xi^{1, K}_\text{prod}: \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); andan approximation of the integral in \(\eqref{transformed_localized_rhs}\) by a numerical

**quadrature**: given any unary element integrand \(\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}\nabla F_K(\hat{x}_q)|\,\omega_q\,\Xi^{1,K}(\hat{\psi}_\hat{i}, \hat{x}_q) \approx \int_\hat{K} \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 \(\Xi_\text{prod}^{1,K}\)).

Given local functionals, the purpose of the `VectorFunctional`

in `dune-gdt`

is to assemble \(\underline{l_h}\) from \(\eqref{algebraic_problem}\) 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 \(\underline{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: \Omega \to \mathbb{R}\), given an ansatz function \(\hat{\varphi}\), a test function \(\hat{\psi}\) and a point \(\hat{x} \in \hat{K}\), this integrand is determined by \(\Xi_\text{prod}^{2,K}: \mathbb{P}^k(\hat{K}) \times \mathbb{P}^k(\hat{K}) \times \hat{K} \to \mathbb{R}\), \(\hat{\varphi}, \hat{\psi}, \hat{x} \mapsto w^K(\hat{x})\,\hat{\varphi}(\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 \(\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, \(\Xi_\text{prod}^{1,K}: \mathbb{P}^k(\hat{K}) \times \hat{K} \to \mathbb{R}\), \(\hat{\psi}, \hat{x} \mapsto \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 \(\underline{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, \varphi_j\) follows in a similar manner. We obtain by localization

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

where \(\kappa^K\) denote the *local function* of \(\kappa\) as above, and where \(\nabla_K\hat{\varphi}_\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{align} \Xi_\text{laplace}^{2, K}: \mathbb{P}^k(\hat{K}) \times \mathbb{P}^k(\hat{K}) \times \hat{K} &\to \mathbb{R}\\ \hat{\varphi}, \hat{\xi}, \hat{x} &\mapsto \big(\kappa^K(\hat{x})\,\nabla_K\hat{\varphi}(\hat{x})\big)\cdot\nabla_K\hat{\psi}(\hat{x}) \end{align}\)$ andan approximation of the integral by a numerical

**quadrature**: given any binary element integrand \(\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{align} a_h^K(\psi_i, \varphi_j) := \sum_{q = 1}^Q |\text{det}\nabla F_K(\hat{x}_q)|\,\omega_q\,\Xi^{2,K}(\hat{\psi}_\hat{i}, \hat{\varphi}_\hat{j}, \hat{x}_q) \approx \int_\hat{K} \Xi^{2,K}(\hat{\psi}_\hat{i}, \hat{\varphi}_\hat{j}, \hat{x})\,\text{d}\hat{x} = a^K(\psi_i, \varphi_i), \end{align}\) 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 \(\underline{a_h}\) from \(\eqref{algebraic_problem}\) 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 \(\underline{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 \(\underline{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 \(\eqref{discrete_variational_problem}\).

```
u_h_vector = a_h.apply_inverse(l_h.vector)
```

### postprocessing the solution#

To make use of the DoF vector of the approximate solution, \(\underline{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 \(\eqref{pde}\) from above, but with non-homogeneous Dirichlet boundary values. That is:

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

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

The variational problem associated with \(\eqref{inhomogeneous_pde}\) reads: find \(u \in H^1(\Omega)\), such that

with the same bilinear form \(a\) and linear functional \(l\) as above in \(\eqref{a_and_l}\).

### Dirichlet shift#

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

in the sense of traces. We then have for the solution \(u \in H^1(\Omega)\) of \(\eqref{inhomogeneous_weak_formulation}\), that

for some \(u_0 \in H^1_0(\Omega)\). Thus, we may reformulate \(\eqref{inhomogeneous_weak_formulation}\) as follows: Find \(u_0 \in H^1_0(\Omega)\), such that

or equivalently

We have thus shifted problem \(\eqref{inhomogeneous_weak_formulation}\) 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 \(\eqref{pde}\) with:

\(d = 2\)

\(\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(\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^1_0(\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`