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
in a weak sense, where
The variational problem associated with
where the bilinear form
respectively.
Consider for example
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
where
A basis of
of order
the invertible affine reference map
for all , such that , where is the reference element associated with ,a set of Lagrangian shape functions
, each associated with a vertex of the reference element anda DoF mapping
for all , which associates the local index of a vertex of the reference element with the global index of the vertex in the grid .
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
To obtain the algebraic analogue to
replace the bilinear form and functional by discrete counterparts acting on
, namely and (the construction of which is detailed further below) andassemble the respective basis representations of
and w.r.t. the basis of into a matrix and vector , given by and , respectively, for ;and obtain the restrictions of
and to by modifying all entries associated with basis functions defined on the Dirichlet boundary: for each index , where the Lagrange-point defining the basis function lies on the Dirichlet boundary , we set the th row of to a unit row and clear the th entry of .
The algebraic version of
After solving the above linear system, we recover the solution of
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
we first consider local functionals (such as
where
Note that, apart from the integration domain (
This leads us to the definition of a local functional in dune-gdt
: ignoring the user input (namely the data function
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 , the integrand is determined by , , which is modelled byLocalElementProductIntegrand
indune-gdt
(see below); andan approximation of the integral in
by a numerical quadrature: given any unary element integrand , and quadrature points and weights , we approximate , which is modelled byLocalElementIntegralFunctional
indune-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
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 ).
Given local functionals, the purpose of the VectorFunctional
in dune-gdt
is to assemble
creating an appropriate vector of length
iterating over all grid elements
localizing the basis of
w.r.t. each grid elementevaluating the local functionals
for each localized basis functionadding the results to the respective entry of
, determined by the DoF-mapping of the discrete function spaceContinuousLagrangeSpace
In our example, we define
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 aGridFunction
, which accepts all kind of functions, discrete functions or numbers. ThusGF(grid, 1)
creates a grid function which is localizable w.r.t. each grid elements and evaluates to 1, whenever evaluated; whereasGF(grid, f)
, when localized to a grid element and evaluated at a point on the associated reference element, , evaluates to . 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 , given an ansatz function , a test function and a point , this integrand is determined by , . Thus,LocalElementProductIntegrand
is often used in bilinear forms to assemble products (with weight ), which we achieve byLocalElementProductIntegrand(GF(grid, 1))
. However, even with , the integrand still depends on the test and ansatz function. Usingwith_ansatz(GF(grid, f))
, we fix as the ansatz function to obtain exactly the unary integrand we require, which only depends on the test function, , , which is exactly what we need to approximate .the above code creates the vector
(available as thevector
attribute ofl_h
) filled with0
, 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
and by transformation and the chain rule, using
where
Similar to local functionals, a local bilinear form is determined
by a binary element integrand, in our case the
LocalLaplaceIntegrand
(see below) $ $ andan approximation of the integral by a numerical quadrature: given any binary element integrand
, and quadrature points and weights , we approximate which is modelled byLocalElementIntegralBilinearForm
indune-gdt
(see below).
Given local bilinear forms, the purpose of the MatrixOperator
in dune-gdt
is to assemble
creating an appropriate (sparse) matrix of size
iterating over all grid elements
localizing the basis of
w.r.t. each grid elementevaluating the local bilinear form
for each combination localized ansatz and test basis functionsadding the results to the respective entry of
, determined by the DoF-mapping of the discrete function spaceContinuousLagrangeSpace
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 functionkappa
to a matrix-valuedGridFunction
the above code creates the matrix
(available as thematrix
attribute ofa_h
) with sparse0
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 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.
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
u_h_vector = a_h.apply_inverse(l_h.vector)
postprocessing the solution#
To make use of the DoF vector of the approximate solution, 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
Let
in a weak sense, where
The variational problem associated with
with the same bilinear form
Dirichlet shift#
Suppose
in the sense of traces. We then have for the solution
for some
or equivalently
We have thus shifted problem
Consider for example
Except for
Dirichlet interpolation#
To obtain 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
assembling the shifted system#
We assemble the unconstrained matrix- and vector representation of
and w.r.t. 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
.
dirichlet_constraints.apply(a_h.matrix, rhs_vector)
Then, we solve the shifted and constrained system for the DoF vector of
.
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