AMDiS Workshop 2021

layout: true
class: tud-dark, typo, splash
---

<div class="frontpage">
  <small>Dr. Simon Praetorius</small>
  <h1>AMDiS Workshop 2021<br>
  <small>Institute of Scientific Computing</small></h1>

<small>Apr 2021</small>
</div>
---
layout: true
class: tud-light, typo

---
class: center, middle
# Introduction
---

# What is this course about?

### 1. DUNE: The **D**istributed and **U**nified **N**umerics **E**nvironment

> [...] a modular toolbox for solving partial differential
> equations (PDEs) with grid-based methods [...] - http://www.dune-project.org

### 2. AMDiS: **A**daptive **M**ulti-**Di**mensional **S**imulations

> A Finite-Element discretization module allowing for fast prototyping of PDEs on adaptively
> refined grids - https://gitlab.mn.tu-dresden.de/amdis/amdis-core

.footnote[Slides are based on [DUNE-PDELab Course 2021](https://dune-pdelab-course.readthedocs.io/en/latest/index.html)]

---

# The DUNE Framework

.pure-g[.pure-u-1-2[
- Modules
  * Code is split into separate modules.
  * Applications use only the modules they need.
  * Modules are sorted according to level of maturaty.
  * Everybody can provide her/his own modules.
- Portability
- Open Development Process
- Free Software Licence
].pure-u-1-2[
<img src="images/dune-modules.png">
]]

.footnote[P. Bastian, M. Blatt, A. Dedner, C. Engwer, R. Klöfkorn, M. Ohlberger, O. Sander.
A Generic Grid Interface for Parallel and Adaptive Scientific Computing. (2008)]

---

# A brief history of DUNE

.center[<img src="images/dune-history.png" width="80%">]

---

# DUNE Release 2.7.1

**dune-common**
: foundation classes, infrastructure

**dune-geometry**
: geometric mappings, quadrature rules, reference element

**dune-grid**
: grid interface, grid construction, visualization

**dune-istl**
: (Iterative Solver Template Library) generic sparse matrix/vector classes, solvers (Krylov methods, AMG, etc.)

**dune-localfunctions**
: generic interface for local finite element functions. Abstract definition following Ciarlet.
  Collection of different finite elements.

.center[<img src="images/dune-logo.png" width="80%">]

---

# DUN(E)iverse

- modular structure
- write your own DUNE modules
- available under different licenses
- **Discretization Modules**
  * `dune-pdelab`: Discretization module based on `dune-localfunctions`.
  * `dune-fem`: Discretization module based on `dune-localfunctions`.
  * `dune-functions`: A new initiative to provide unified interfaces for functions and function spaces.
  * **AMDiS**: based on `dune-functions`
- **Additional Grid Implementations**...
- **Extension Modules**...

---

# DUN(E)iverse

- modular structure
- write your own DUNE modules
- available under different licenses
- **Discretization Modules**
- **Additional Grid Implementations**
  * `dune-grid-glue`: allows to compute overlapping and nonoverlapping couplings of Dune grids, as
    required for most domain decomposition algorithms.
  * `dune-subgrid`: allows you to work on a subset of a given DUNE grid.
  * `dune-foamgrid`: non-manifold grids of 1d or 2d entities in higher-dim. world.
  * `dune-prismgrid`: is a tensorgrid of a 2D simplex grid and a 1D grid.
  * `dune-curvedgrid`: Arbitrary parametrization of your element geometries.
- **Extension Modules**...

---

# DUN(E)iverse

- modular structure
- write your own DUNE modules
- available under different licenses
- **Discretization Modules**
- **Additional Grid Implementations**
- **Extension Modules**
  * `dune-python`: python bindings for central DUNE components
  * `dune-typetree`: classes to organize types in trees
  * `dune-dpg`: construct optimal Discontinuous-Petrov-Galerkin test spaces
  * `dune-tpmc`: cut-cell construction using level-sets

---

# DUN(E)iverse

- modular structure
- write your own DUNE modules
- available under different licenses
- **Discretization Modules**
- **Additional Grid Implementations**
- **Extension Modules**

## Goals: allow people to...
- get credit for their innovations
- experiment without breaking the core
- develop at different speeds

---

# Getting the DUNE libraries and sources

## 1. Install Linux Package

- DUNE is available as system package
- Use your package manager
- Either download `.tar.gz` files with packed sources from https://dune-project.org
- Or install directly in Debian systems

```bash
sudo apt-get install libdune-common-dev libdune-geometry-dev libdune-grid-dev ...
```

> **Note:** Check the DUNE version your distribution provides. Should be >= 2.7!

> **Better:** For more control about DUNE versions and to use any existing module, download sources directly

---

# Getting the DUNE libraries and sources

## 2. Download the Sources from Git
- Central source repository: `<gitlab>:=`https://gitlab.dune-project.org
- Modules grouped into `core`, `staging`, `extensions`, discretization/user modules...
- **Suggestion:** Clone all dune modules into same base directory, called `<dune-base>` in the following.

```bash
$ mkdir <dune-base> && cd <dune-base>
$ git clone <gitlab>/core/dune-[common,geometry,grid,istl,localfunctions]
$ git clone <gitlab>/staging/dune-[functions,typetree,uggrid]
$ git clone <gitlab>/extensions/dune-[alugrid,foamgrid,vtk]
```

- Also DUNE modules in other repositories, e.g.,

```bash
$ git clone https://gitlab.mn.tu-dresden.de/iwr/dune-[curvedgrid,curvedgeometry,gmsh4]
```

> **Remark:** Clone a specific version of dune, e.g., `--branch releases/2.7`

---

# Dependency Management

## dunecontrol
- control of module-interplay
- suggestions & dependencies
- integrates with cmake & git
- works with Linux, Mac and Mingw

> **Note:** Dependencies should form a DAG

`dunecontrol cmake`
: configure packages via `cmake`, include necessary paths

` dunecontrol make`
: build packages in correct order

... works without `make install`

.footnote[Image source: gnome]

---

# Module Management

## duneproject
- Create a new module
- Specify dependencies
- Additional information in `dune.module` file:
  * Maintainer
  * Module Version
  * URL, Description
  * python package requirements
- Interactive tool

`duneproject [<project_name>]`
: Create a new directory `<project_name>` containing the module information

.footnote[Image source: gnome]

---

# Building the Source

- Run `dunecontrol` from within `<dune-base>`
- It searches all subdirectories recursively for a file `dune.module` describing each module and
  its dependencies
- Script is provided by `dune-common` module $\rightarrow$ create alias/link!

```bash
$ export PATH=<dune-base>/dune-common/bin:$PATH   # maybe put into .bashrc

$ cd <dune-base>
$ dunecontrol cmake   # configure all found modules
$ dunecontrol make    # build all found modules
```

> **Note:** Both steps can be combined to `dunecontrol all`

---

# Building the Source

```bash
$ cd <dune-base>
$ dunecontrol --opts=config.opts cmake   # configure all found modules
$ dunecontrol --opts=config.opts make    # build all found modules
```

```bash
$ cat config.opts
CMAKE_FLAGS="-DCMAKE_BUILD_TYPE=Release"
MAKE_FLAGS="-j4"
```

---

# Building the Source

```bash
$ cd <dune-base>
$ dunecontrol --opts=config.opts --module=dune-grid cmake   # configure dune-grid + deps
$ dunecontrol --opts=config.opts --module=dune-grid make    # build dune-grid + deps
```

---

# Building the Source

```bash
$ cd <dune-base>
$ dunecontrol --opts=config.opts --only=dune-grid cmake   # configure dune-grid
$ dunecontrol --opts=config.opts --only=dune-grid make    # build dune-grid
```
or
```bash
$ cd <dune-base>/dune-grid
$ cmake build-cmake/                                # configure dune-grid
$ cmake --build build-cmake/ [--target <target>]    # build dune-grid
```

---

# Resources for Help

- DUNE - The Distributed and Unified Numerics Environment, O. Sander, (2020)
- The DUNE website https://www.dune-project.org
- DUNE-PDELab Course:<br>https://dune-pdelab-course.readthedocs.io/en/latest/index.html
- Mailing list for user questions and discussions:<br>
  [dune@lists.dune-project.org](dune@lists.dune-project.org),
  [Un-/Subscribe](https://lists.dune-project.org/mailman/listinfo/dune)
- The DUNE Core-developer group:
.pure-g[.pure-u-1-2[.font-m[<ul>
  <li>Peter Bastian (Universität Heidelberg)</li>
  <li>Markus Blatt (HPC-Simulation-Software & Services)</li>
  <li>Ansgar Burchardt (TU Dresden)</li>
  <li>Andreas Dedner (University of Warwick, UK)</li>
  <li>Christian Engwer (Universität Münster)</li>
  <li>Jorrit Fahlke (Jö)</li>
  </ul>]].pure-u-1-2[.font-m[<ul>
  <li>Carsten Gräser (Freie Universität Berlin)</li>
  <li>Christoph Grüninger</li>
  <li>Dominic Kempf (Universität Heidelberg)</li>
  <li>Robert Klöfkorn (Lund University, Sweden)</li>
  <li>Simon Praetorius (TU Dresden)</li>
  <li>Oliver Sander (TU Dresden)</li>
  </ul>]]]

---
layout: true
class: tud-dark, typo
---
class: center, middle, tud-dark
# Exercise 1
---

# Exercise 1

1. Find out the requirements and suggestions of the dune module `dune-grid`.
2. Download and build `dune-grid` and all its dependencies.
3. Create an options-file with cmake and make options:
  * build type `Release`
  * Enable all warnings
  * Compile with 4 processors
  * Change the buildsystem from *Makefile* to *Ninja Build*
  * Set a build directory
4. Recompile with the options passed to `dunecontrol`

See also https://dune-project.org/doc/installation/ for some more information
---
layout: true
class: tud-light, typo
---
class: center, middle
# The DUNE Grid Interface

> "People think that computer science is the art of geniuses but the actual reality
> is the opposite, just many people doing things that build on each other, like a
> wall of mini stones." - *Donald E. Knuth*

---

# The DUNE Grid Interface

Weak formulation of boundary value problem:

\[
\text{Find }u\in U\;\text{ s.t. }\;a(u, v) = l(v)\; \forall v\in V
\]
with $a(u, v)$ and $l(v)$ are (bi)linear forms, e.g.,

\[
a(u,v) = \int_\Omega \nabla u \cdot \nabla v\, dx,\quad l(v) = \int_\Omega f(x)\,v\,dx,
\]
with spatial domain $\Omega\in\mathbb{R}^d$

Grids are necessary for at least three reasons:
1. Piecewise description of the complicated domain $\Omega$
2. Piecewise approximation of functions (by polynomials)
3. Piecewise computation of integrals (by numerical quadrature)

---

# The DUNE Grid Interface

- The DUNE Grid Module is one of the five DUNE Core Modules.
- DUNE wants to provide an interfaces for grid-based methods. Therefore the concept of a *Grid* is the central part of DUNE.
- `dune-grid` provides the interfaces, following the concept of a Grid.
- Its implementation follows the three *design principles* of DUNE:
  * **Flexibility:** Separation of data structures and algorithms.
  * **Efficiency:** Generic programming techniques.
  * **Legacy Code:** Reuse existing finite element software.

---

# Designed to support a wide range of Grids

.center[<img src="images/dune-grids.png" width="80%">]

---

# DUNE Grid Interface Features

- Provide abstract interface to grids with:
  * Arbitrary dimension embedded in a world dimension,
  * multiple element types,
  * conforming or nonconforming,
  * hierarchical, local refinement,
  * arbitrary refinement rules (conforming or nonconforming),
  * parallel data distribution and communication,
  * dynamic load balancing.
- Reuse existing implementations (ALU, UG, Alberta) + special implementations (YaspGrid, FoamGrid).
- Meta-Grids built on-top of the interface (GeometryGrid, SubGrid, MultiDomainGrid, SubdomainGrid, CurvedGrid)

.footnote[P. Bastian, M. Blatt, A. Dedner, C. Engwer, R. Klöfkorn, M. Ohlberger, O. Sander.
A Generic Grid Interface for Parallel and Adaptive Scientific Computing. (2008)]

---

# The Grid

A formal specification of grids is required to enable an accurate description of the grid interface.

.pure-g[.pure-u-1-2[
  .center[<img src="images/hierarchic-grid.png">]
].pure-u-1-2[
In DUNE a Grid is always a hierarchic grid of
dimension `d`, existing in a `w` dimensional space.
The Grid is parametrized by
- the dimension `d`,
- the world dimension `w`
- and the maximum level `J`.
]]

---

# The Grid... A Container of Entities...

.pure-g[.pure-u-1-2[
  .center[<img src="images/grid-cell.png" width="80%">]
].pure-u-1-2[
- vertices
- edges
- faces
- ...
- cells
]]

In order to do dimension independent programming, we need a dimension independent
naming for different entities.

---

# The Grid... A Container of Entities...

.pure-g[.pure-u-1-2[
  .center[<img src="images/grid-cell-codim.png" width="80%">]
].pure-u-1-2[
- vertices (Entity `codim = d`)
- edges (Entity `codim = d - 1`)
- faces (Entity `codim = d - 2`)
- ...
- cells (Entity `codim = 0`)
]]

In order to do dimension independent programming, we need a dimension independent
naming for different entities.

We distinguish entities according to their codimension.
Entities of `codim = c` contain subentities of `codim = c + 1`. This gives a recursive
construction down to `codim = d`.

---

# The DUNE Grid Interface

The DUNE Grid Interface is a collection of classes and methods
```c++
#include <dune/grid/yaspgrid.hh>

// ...

// Create a 2d structured grid of [0,1] x [0,1]
using Grid = Dune::YaspGrid<2>;
Grid grid{ {1.0, 1.0}, {4, 4} };

// traverse the grid
auto gv = grid.leafGridView();
for (auto const& cell : elements(gv)) {
  // do something
}
```

*We will now get to know the most important classes and see how they interact.*

---

# Modifying a Grid

The DUNE Grid interface follows the View-only Concept.

## View-Only Concept
- Views offer (read-only) access to the data
  * Read-only access to grid entities allow the consequent use of `const`.
  * Access to entities is only through iterators for a certain view.
    - This allows on-the-fly implementations.
- Data can only be modified in the primary container (*the Grid*)

## Modification Methods:
- Global Refinement
- Local Refinement & Adaption
- Load Balancing

---

# Views to the Grid
A Grid offers two major views:

.pure-g[.pure-u-1-2[
<img src="images/grid-levels.png" width="90%"><br><br><img src="images/gridviews.png" width="90%">
].pure-u-1-2[
**levelwise:** all entities associated with the same level.
<div style="visibility:hidden;margin-bottom:10px;"><i>Note: not all levels must cover the whole
domain.</i></div>

**leafwise:** all leaf entities (entities which are not refined).
<div style="visibility:hidden">The leaf view can be seen as the projection of a
levels onto a flat grid. It again covers the whole
domain.</div>
]]

---

# Views to the Grid
A Grid offers two major views:

.pure-g[.pure-u-1-2[
<img src="images/grid-levels2.png" width="90%"><br><br><img src="images/gridviews.png" width="90%">
].pure-u-1-2[
**levelwise:** all entities associated with the same level.
<div style="margin-bottom:10px;"><i>Note: not all levels must cover the whole
domain.</i></div>

**leafwise:** all leaf entities (entities which are not refined).
<div>The leaf view can be seen as the projection of a
levels onto a flat grid. It again covers the whole
domain.</div>
]]

---

# Views to the Grid
## `Dune::GridView`

- The `Dune::GridView` class consolidates all information depending on the current View.
--

- Every Grid must provide
  * `Grid::LeafGridView` and
  * `Grid::LevelGridView`.
--

- The Grid creates a new view every time you ask it for one, so you need to store a copy of it.
- Accessing the Views:
  * `Grid::leafGridView()` and
  * `Grid::levelGridView(int level)`.

---

# Iterating over grid entities

Typically, most code uses the grid to iterate over some of its entities (e.g. cells) and perform
some calculations with each of those entities.
- GridView supports iteration over all entities of one codimension.
- Iteration uses C++11 range-based for loops:

```c++
for (auto const& cell : elements(gv)) {
  // do some work with cell
}
```

---

# Iteration functions

You can do similar calls for other entity types:
```c++
// Iterates over cells (codim = 0)
for (const auto& c : elements(gv))
// Iterates over vertices (dim = 0)
for (const auto& v : vertices(gv))
// Iterates over facets (codim = 1)
for (const auto& f : facets(gv))
// Iterates over edges (dim = 1)
for (const auto& e : edges(gv))

// Iterates over entities with a given codimension (here : 2)
for (const auto& e : entities(gv, Dune::Codim<2>{}))
// Iterates over entities with a given dimension (here : 2)
for (const auto& e : entities(gv, Dune::Dim<2>{}))
```

---
layout: true
class: tud-light, typo
---
class: center, middle
# Entities
---

# Entities

Iterating over a grid view, we get access to the entities.
```c++
for (const auto& cell : elements(gv)) {
  cell.?????(); // what can we do here ?
}
```

- Entities cannot be modified.
- Entities can be copied and stored (but copies might be expensive!).
- Entities provide topological and geometrical information.

---

# Entities
### An Entity `T` provides both topological information
- Type of the entity (triangle, quadrilateral, etc.).
- Relations to other entities.

### and geometrical information
- Position of the entity in the grid.

.pure-g[.pure-u-2-3[
  <img src="images/refmap.png" width="90%">
].pure-u-1-3[
### Entity `T` is defined by...
  - Reference Element $T_\text{ref}$
  - Transformation $F_T$

`GridView::Codim<c>::Entity` implements the entity concept.
]]

---

# Storing Entities

`GridView::Codim<c>::Entity`
- Entities can be copied and stored like any normal object.
- Important: There can be multiple entity objects for a single logical grid entity (because they can be copied)
- Memory expensive, but fast.

`GridView::Codim<c>::EntitySeed`
- Store minimal information to find an entity again.
- Create like this:
```c++
auto entity_seed = entity.seed();
```
- The grid can create a new Entity object from an EntitySeed:
```c++
auto entity = grid.entity(entity_seed);
```
- Memory efficient, but run-time overhead to recreate entity.

---

# Reference Elements
`Dune::GeometryType` identifies the type of the entity's reference element.
`cell.type()` returns the GeometryType of an entity.

---

# Geometry Types
`GeometryType` is a simple identifier for a reference element
- Obtain from entity or geometry object using `.type()`
- GeometryType for specific reference elements in namespace `Dune::GeometryTypes`:

```c++
#include <dune/geometry/type.hh>
...
namespace GeometryTypes = Dune::GeometryTypes;
Dune::GeometryType gt;

gt = GeometryTypes::vertex;
gt = GeometryTypes::line;
gt = GeometryTypes::triangle;
gt = GeometryTypes::square;
gt = GeometryTypes::hexahedron;
gt = GeometryTypes::cube(dim);
gt = GeometryTypes::simplex(dim);
```

- GeometryTypes are cheap, always store and pass around copies (don't use references)

---

# ReferenceElement (I)
A reference element provides topological and geometrical information about the embedding of subentities:

- Numbering of subentities within the reference element
- Geometrical mappings from reference elements of subentities to the current reference element

---

# ReferenceElement (II)
- Reference elements are templated on the dimension and the coordinate field type
```c++
#include <dune/geometry/referenceelement.hh>
...
Dune::ReferenceElement<double, dim> ref_el = ...;
```
- The function `Dune::referenceElement()` will extract the reference element from most objects that have one:
```c++
auto ref_el = Dune::referenceElement(entity.geometry()); // or
auto ref_el = Dune::referenceElement(entity);
```
When using this function, you don’t have to figure out the template parameters.
- ReferenceElements are cheap, always store and pass around copies (don't use references)

---

# Geometry
Transformation $F_T$
- Maps from one space to an other.
- Main purpose is to map from the reference element to global coordinates.
- Provides transposed inverse of the Jacobian $J^{-T}_F$.
- Gramian determinant of transformation Jacobian for quadrature.

.center[<img src="images/refmap.png" width="70%">]

---

# Geometry Interface (I)
- Obtain Geometry from entity
```c++
auto geo = entity.geometry();
```
- Convert local coordinate to global coordinate
```c++
auto x_global = geo.global(x_local);
```
- Convert global coordinate to local coordinate
```c++
auto x_local = geo.local(x_global);
```

- Get center of geometry in global coordinates
```c++
auto center = geo.center();
```
- Get number of corners of the geometry (e.g. 3 for a triangle)
```c++
auto num_corners = geo.corners();
```
- Get global coordinates of i-th geometry corner ($0 \leq i < $`geo.corners()`)
```c++
auto corner_global = geo.corner(i);
```

---

# Geometry Interface (II)

- Get type of reference element
```c++
auto geometry_type = geo.type(); // square , triangle , ...
```
- Find out whether geometry is affine
```c++
if (geo.affine()) {
  // do something optimized
}
```
- Get volume of geometry in global coordinate system
```c++
auto volume = geo.volume();
```
- Get integration element for a local coordinate (required for numerical integration)
```c++
auto mu = geo.integrationElement(x_local);
```

---

# Gradient Transformation
Assume $f : \Omega \to \mathbb{R}$

evaluated on a cell $T$, i.e., $f(F_T(\hat{x}))$.

The gradient of $f$ is then given by \[J_F^{-T}(\hat{x})\hat{\nabla} f(F_T(\hat{x}))\]
```c++
auto x_global = geo.global(x_local);
auto J_inv = geo.jacobianInverseTransposed(x_local);
auto tmp = grad(f)(x_global); // grad(f) supplied by user
<grad_type> gradient;         // something like FieldVector<double,dimworld>
J_inv.mv(tmp, gradient);
```

---

# Obtaining Quadrature Rules
Numerical quadrature rules given by
\[
  \int_{T_\text{ref}} \hat{f}(\hat{x})\,\text{d}\hat{x} \approx \sum_{i=0}^{N} w_i \hat{f}(\hat{x}_i)
\]

- `dune-geometry` provides pre-defined quadrature rules for common geometry types:
```c++
#include <dune/geometry/quadraturerules.hh>
...
Dune::GeometryType gt = ...;
auto const& rule = Dune::QuadratureRules<double, dim>::rule(gt, order);
```
- The rule factory is parameterized by the number type (typically use `Grid::ctype`)
  and the dimension of the integration domain, e.g. `Entity::mydimension`
- The rule is exact for polynomials up to the given `order`
- Use `auto const&` for the type of the rule to avoid expensive copies
- Optional third parameter to select type of rule (Jacobi, Legendre, Lobatto)

---

# Using Quadrature Rules
- A `QuadratureRule` is a range of `QuadraturePoint`s.
- `QuadraturePoint` provides weight and position:
  * `QuadraturePoint::weight()`
  * `QuadraturePoint::position()`

### Example
```c++
auto fLocal = some_function_to_integrate(...);
double integral = 0.0;
auto geo = cell.geometry();
for (const auto& qp : rule)
{
  integral += fLocal(qp.position()) ⁎ qp.weight() * geo.integrationElement(qp.position());
}
```

> **Attention:** When integrating over cells in a grid, keep in mind that the quadrature
> point coordinates are local to the reference element.
---
layout: true
class: tud-dark, typo
---
class: center, middle, tud-dark
# Exercise 2
---

# Exercise 2

1. Create a new dune module `dune-grid-exercise` using the tool `duneproject`
2. Create a structured grid for the domain $\Omega=[0,2]\times[0,2]$ with 8 elements in each
   direction. Therefore, use either `YaspGrid` or `SPGrid`.
3. Implement a function that computes the numerical integration of a given function $f:\Omega\to\mathbb{R}$
   over the domain $\Omega$. Therefore, traverse the grid and create a `QuadratureRule` on each grid
   entity. For each quadrature point, evaluate the function $f$ in global coordinates by mapping element-local
   coordinates using the element geometry.

---

# Exercise 2
Concrete setup:
- Compute the numerical quadrature of the function

\[
   f(x) = \sin(x_0) \cdot \cos(x_1)
\]
---
layout: true
class: tud-light, typo
---
class: center, middle
# Intersections
---

# Intersections

- Grids may be non conforming.
- Entities can intersect with neighbours and boundary.
- Represented by Intersection objects.
- Intersections hold topological and geometrical information.
- Intersections depend on the view.
- **Note:** Intersections are always of codimension 1!

---

# Intersection Interface
- Is this an intersection with the domain boundary?
```c++
bool b = intersection.boundary();
```
- Is there an entity on the outside of the intersection?
```c++
bool b = intersection.neighbor();
```
- Get the cell on the inside
```c++
auto inside_cell = intersection.inside();
```
- Get the cell on the outside
```c++
// Do this only if intersection.neighbor() == true
auto outside_cell = intersection.outside();
```

---

# Intersection: Geometries

- Get mapping from intersection reference element to global coordinates
```c++
auto world_geo = intersection.geometry();
```
- Get mapping from intersection reference element to reference element of inside cell
```c++
auto inside_geo = intersection.geometryInInside();
```
- Get mapping from intersection reference element to reference element of outside cell
```c++
auto outside_geo = intersection.geometryInOutside();
```

---

# Intersection: Normals

- Get unit outer normal for local coordinate.
```c++
auto unit_outer_normal = intersection.unitOuterNormal(x_local);
```
- Get unit outer normal for center of intersection (good for affine geometries).
```c++
auto unit_outer_normal = intersection.centerUnitOuterNormal();
```
- Get unit outer normal scaled with integration element (convenient for numerical quadrature).
```c++
auto integration_outer_normal
  = intersection.integrationOuterNormal(x_local);
```

---

# Iterating over intersections
## Example
In order to iterate over the intersections of a given grid cell with respect to some
`GridView`, use a range-based for loop with the argument `intersections(gv,cell)`.

The following code iterates over all cells in a GridView and over all intersections of
each cell:
```c++
for (const auto& cell : elements(gv)) {
  for (const auto& is : intersections(gv, cell)) {
    if (is.boundary()) {
      // handle potential Neumann boundary
    }
    if (is.neighbor()) {
      // code for Discontinuous Galerkin or Finite Volume
    }
  }
}
```

---

# Sequential finite volume solver
*Elementwise divergence* of a vector field:
\[
  \int_T \nabla\cdot f(x)\,\text{d}x = \int_{\partial T} f(x)\cdot\mathbf{n}_{\partial T}\,\text{d}s
\]

Consider the *first-order linear PDE*
$
  \partial_t u + \nabla\cdot(v u) = 0
$
with given vector field $v(x)$ and unknown solution $u(x, t)$.
The structure-explicit cell-centered finite volume method reads
\[
  \bar{u}_T^{k+1} = \bar{u}_T^k - \frac{\Delta t}{|T|} \sum_{(T,T')\in I(T)} \phi(v\cdot\mathbf{n}_T, \bar{u}_T^k, \bar{u}_{T'}^k) |I(T,T')|
\]
with the numerical flux function $\phi$ chosen as upwind flux here.

> keep this in mind for exercise 3

---
layout: true
class: tud-light, typo
---
class: center, middle
# Attaching Data to the Grid
---

# Attaching Data to the Grid

For computations we need to associate data with grid entities:
- spatially varying parameters
- entries in the solution vector or the stiffness matrix
- polynomial degree for p-adaptivity
- status information during assembling
- ...

... meaning we want to
- associate data with subsets of entities
- subsets could be "vertices of level l", "faces of leaf elements" ...
- data should be stored in arrays for efficiency
- associate index/id with each entity

---

# Indices and Ids
**Index Set:** Provides a map $m : E \to \mathbb{N}_0$, where $E$ is a subset of the entities of a grid
view.

We define the subsets $E_g^c$ of a grid view
\[
  E_g^c = \{e \in E \vert e\text{ has codimension }c\text{ and geometry type }g\}.
\]
- unique within the subsets $E_g^c$.
- consecutive and zero-starting within the subsets $E_g^c$.
- distinct leaf and a level index.

**Id Set:** provides a map $m : E \to \mathbb{I}$, where $\mathbb{I}$ is a discrete set of ids.
- unique within $E$.
- ids need not to be consecutive nor positive.
- persistent with respect to grid modifications.

---

# Example: Store the lengths of all edges
The following example demonstrates how to
- query an index set for the number of contained entities of a certain codimension
  (so that we can allocate a vector of correct size).
- obtain the index of a grid entity from an index set and use it to store associated
  data.

```c++
// Get the IndexSet associated to a GridView gv
const auto& index_set = gv.indexSet();

// Create a vector with one entry for each edge
auto edge_lengths = std::vector<double>(index_set.size(1));

// Loop over all edges and store their length
for (const auto& edge : edges(gv))
  edge_lengths[ index_set.index(edge) ] = edge.geometry().volume();
```

---

# Example: 2D Multi-Element Grid – Indices
Locally refined grid:

---

# Example: 2D Multi-Element Grid – Indices
Locally refined grid:

---

# Example: 2D Multi-Element Grid – Indices
Locally refined grid: **Indices**

.center[Consecutive index for vertices]

---

# Example: 2D Multi-Element Grid – Indices
Locally refined grid: **Indices**

.center[... and cells]

---

# Example: 2D Multi-Element Grid – Indices
Locally refined grid: **Indices**

.center[Old cell indices on coarse grid level]

---

# Example: 2D Multi-Element Grid – Indices
Locally refined grid: **Indices**

.center[Consecutive cell indices on coarse and refined grid]

---

# Example: 2D Multi-Element Grid – Indices
Locally refined grid: **Ids**

.center[Persistent Ids on coarse and refined grid]

---

# Mapper
**Mappers extend** the functionality of **Index Sets**.
- associate data with an arbitrary subsets $E' \subseteq E$
  of the entities $E$ of a grid.
- the data $D(E')$ associated with $E'$ is stored in an array.
- map from the consecutive, zero-starting index $I_{E'} = \{0,\ldots , |E'|-1\}$ to the data
  set $D(E')$.

Mappers can be easily implemented upon the Index Sets and Id Sets.

## Example
You will be using the
```c++
Dune::MultipleCodimMultipleGeomTypeMapper<GridView>
```

---

```c++
#include <dune/grid/common/mcmgmapper.hh>
...
using GridView = SomeGrid::LeafGridView;

// Layout description
Dune::MCMGLayout layout = [](Dune::GeometryType gt, int griddim) {
  return gt.dim() == griddim;
};

// mapper for elements (codim==0) on leaf
using Mapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
Mapper mapper{gridview, layout};

// iterate over the leaf
for (const auto& entity : elements(gridview))
{
  int index = mapper.index(entity);
  // iterate over all intersections of this cell
  for (const auto& i : intersections(gridview, entity))
  {
    // neighbor intersection
    if (i.neighbor()) {
      int nindex = mapper.index(i.outside());
      sparsityPattern[index].insert(nindex);
    }
  }
}
```
---
layout: true
class: tud-light, typo
---
class: center, middle
# Input and Output
---

# Input and Output

- We need to provide a grid either by geometric description or loaded from file
- Powerful grid generators exist, exporting grids in multiple formats
- For the visualization of grid we also need to write the geometry and maybe attached data
  in a common file format

## In this section:
1. Grid factories
2. Structured grid generation
3. Reading grids from file
4. Writing output to a file

---

# Grid Factories

- A grid is a combination of elements, described by their vertex coordinates + a connectivity
  of the element nodes.
- All Dune grids provide a **Grid Factory** that allows to create a grid by these two information
```c++
Dune::GridFactor<SomeGrid> factory;
```

- The `factory` then provides two main methods:
```c++
void insertVertex (const FieldVector<ct,dimworld>& pos);
void insertElement (const GeometryType& type, const std::vector<unsigned int>& vertices)
```

- Construct the grid after inserting vertices and element connectivity:
```c++
std::unique_ptr<SomeGrid> createGrid ();
```

---

# Grid Factories
## Example
```c++
Dune::GridFactory<SomeGrid> factory;

factory.insertVertex({0, 0, 0});
factory.insertVertex({1, 0, 0});
factory.insertVertex({0, 1, 0});
factory.insertVertex({1, 1, 0});
factory.insertVertex({0, 0, 1});
factory.insertVertex({1, 0, 1});
factory.insertVertex({0, 1, 1});
factory.insertVertex({1, 1, 1});

namespace GT = Dune::GeometryTypes;
factory.insertElement(GT::tetrahedron, {0, 1, 3, 7});  // insertion-index 0
factory.insertElement(GT::tetrahedron, {0, 5, 1, 7});  // insertion-index 1
factory.insertElement(GT::tetrahedron, {0, 4, 5, 7});  // insertion-index 2
factory.insertElement(GT::tetrahedron, {0, 6, 4, 7});  // insertion-index 3
factory.insertElement(GT::tetrahedron, {0, 2, 6, 7});  // insertion-index 4
factory.insertElement(GT::tetrahedron, {0, 3, 2, 7});  // insertion-index 5

std::unique_ptr<SomeGrid> gridPtr = factory.createGrid();
```

---

# Grid Factories
- Identify elements in the constructed grid, by **insertion-index**
```c++
unsigned int insertionIndex (const typename Codim<0>::Entity& entity);
```

- **Note:** The insertion-index might be different from the entity index in the Index Set

### Example
```c++
auto const& indexSet = gridPtr->leafIndexSet();
for (auto const& cell : elements(*gridPtr))
{
  auto index = indexSet.index(cell);
  auto insertion_index = factory.insertionIndex(cell);

// in general: index != insertion_index
}
```

---

# Structured Grid Generation
- Some standard Grid Factories exist for rectangular domains
```c++
#include <dune/grid/utility/structuredgridfactory.hh>
...
Dune::StructuredGridFactory<SomeGrid>
```

- These structured factories provide static methods to construct cube/simplex grids of a prescribed
  domain:

```c++
std::unique_ptr<SomeGrid> createCubeGrid (const FieldVector<ct,dimworld>& lowerLeft,
                                          const FieldVector<ct,dimworld>& upperRight,
                                          const std::array<unsigned int,dim>& elements);

std::unique_ptr<SomeGrid> createSimplexGrid (const FieldVector<ct,dimworld>& lowerLeft,
                                             const FieldVector<ct,dimworld>& upperRight,
                                             const std::array<unsigned int,dim>& elements);
```

---
# Structured Grid Generation
## Example
A domain $\Omega=[0,2]\times[0,4]$ with 100 elements in x direction and 200 elements in y direction

```c++
using GridType = Dune::UGGrid<2>;
using Factory = Dune::StructuredGridFactory<GridType>;

auto gridPtr = Factory::createCubeGrid({0.0,0.0}, {2.0,4.0}, {100u,200u});
```

But also simplex elements can be used:

```c++
using GridType = Dune::UGGrid<2>;
using Factory = Dune::StructuredGridFactory<GridType>;

auto gridPtr = Factory::createSimplexGrid({0.0,0.0}, {2.0,4.0}, {100u,200u});
```

---

# Reading Grids from File

External tools like **GMsh** or **ParaView** or CAD tools export their grids in a specific file format.
In order to create a Dune Grid from these file, the content needs to be parsed and a `GridFactory` can
be filled with the vertices and element connectivity.

For some file formats this is already implemented in `dune-grid`:
- GMsh (2.x)
- StarCD
- DGF
- AmiraMesh
- AlbertaGrid

With external modules also other formats are supported, e.g., GMsh (4.x) with `dune-gmsh4` and VTK with `dune-vtk`.

---

# Reading Grids from File
## Example
Reading a grid from a GMsh 2.x file
```c++
#include <dune/grid/io/file/gmshreader>
...
using Reader = Dune::GmshReader<SomeGrid>;
std::unique_ptr<SomeGrid> gridPtr = Reader::read(<filename>);
```

Some readers provide interfaces to also read data, e.g. with the `GMshReader`
```c++
Dune::GridFactory<SomeGrid> factory;
std::vector<int> boundaryTags;
std::vector<int> elementTags;
Reader::read(factory, <filename>, boundaryTags, elementTags);

std::unique_ptr<SomeGrid> gridPtr = factory.createGrid();
```

---

# Writing Output to a File

For the visualization of simulation results, you need to write the grid and attached data to a file.
A common output format is VTK.

Dune provides a `VTKWriter` and also sequence writers for animation files:
```c++
#include <dune/grid/io/file/vtk.hh>
...
Dune::VTKWriter<GridView>;
Dune::VTKSequenceWriter<GridView>;
```

A writer is constructed from a `GridView` and allows to attach data on the cells or on the vertices:
```c++
void addCellData (...);
void addVertexData (...);
```
where the argument is something that can be evaluated on the corresponding grid entities:
- A container of values with contiguous indexing of the entity values
- A localizable function (see later)

---
# Writing Output to a File
## Example
Writing a vector of values attached to the grid cells to a file
```c++
Dune::VTKWriter writer{gv, Dune::VTK::nonconforming};  // write a non-conforming grid

auto u = std::vector<double>(gv.size(0));
writer.addCellData(u, "c", 1);                         // vector contains 1 value per cell

writer.write("concentration", "output", "");           // write into the directory output/
```

A `VTKSequenceWriter` is constructed with the output location:
```c++
Dune::VTKSequenceWriter writer{gv, "concentration", "output", "", Dune::VTK::nonconforming};
writer.addCellData(u, "c", 1);
writer.write(time);
```
---
layout: true
class: tud-dark, typo
---
class: center, middle, tud-dark
# Exercise 3
---

# Exercise 3

The partial differential equation considered in this exercise is the linear transport
equation,
\[
  \partial_t u(x,t) + \nabla\cdot(v(x) u(x,t)) = 0\quad\text{in }\Omega \\
  u(x,t) = u_\text{in}(x,t)\quad\text{on }\Gamma_\text{in}
\]
The unknown solution is denoted by $u(x, t)$ and the velocity field by $v(x)$. The
domain $\Omega$ is some open subset of $\mathbb{R}^d$. For this exercise, we choose `d = 2` where
intersections are 1-D edges. The inflow boundary $\Gamma_\text{in}$ is the set of points $x$ on
the boundary of $\Omega$ for which the velocity vector $v(x)$ points inwards.

We want to numerically solve this equation by a cell-centered finite volume
scheme. We discretize the domain $\Omega$ by a triangulation $\mathcal{T}_h$ and approximate the
solution $u$ by a function $u_h$ that is constant on each cell $T\in\mathcal{T}_h$. We denote the
value of $u_h$ on a cell $T$ by $u_T$.

---

# Exercise 3

Using the explicit Euler time discretization, the scheme can be written as
\[
  u_T(t_{k+1}) = u_T(t_k) - \frac{t_{k+1} - t_k}{|T|} \sum_{e\in\partial T} |e| u^e(t_k) \mathbb{n}_T^e\cdot v^e\,.
\]

The notation is as follows: $|T|$ is the area of the cell $T$, the sum runs over all
intersections $e$ of $T$ with either the boundary or a neighboring cell, $|e|$ is the length
of edge $e$, $\mathbf{n}^e_T$ is the unit outer normal of edge $e$ and $v^e$ is the velocity at the center
of edge $e$. Finally, $u^e$ denotes the upwind concentration. If $\mathbf{n}^e_T \cdot v^e > 0$, this is
$u_T$. Otherwise it is either the concentration in the neighboring cell or given by the boundary
condition $u_\text{in}$, depending on the location of $e$.

---

# Exercise 3

For the concrete setup, consider:
- Grid $\Omega=[0,1]\times[0,1]$ with 100 elements in each direction
- $u_T(t_0) = u_0(center(T))$ with

\[
  u_0(x) = \left\{\begin{array}{ll} 1 & \text{if }\|x-(0.15,0.15)\|_\infty\leq 0.05 \\ 0 & \text{otherwise}\end{array}\right.
\]

- Velocity $v=(1,1)$
- Boundary condition $u_\text{in}=0$
- Time interval $t\in[0,1]$ with timestep size $t_{k+1}-t_k = 0.5\cdot h$ with $h$ the grid size.
---
layout: true
class: tud-light, typo
---
class: center, middle
# DUNE LocalFunctions
---

# DUNE LocalFunctions

Finite-Elements are defined as triple $(T,P,\Sigma)$:
1. An element of the grid $T$
2. The space of functions $P$ over the domain $T$
3. A finite set $\Sigma$ of *linear forms* $\sigma_i:P\to\mathbb{R}$, called *degrees of freedom*

### Example
- Let $T$ be a d-dim simplex
- $P$ the set of polynomials of order $\leq p$
- $\sigma_i:P\to\mathbb{R}, i=0,\ldots,N-1$ are point evaluations at $N$ *Lagrange* points in $T$
- The dual to $\sigma_i$ corresponds to the Lagrange polynomials

$\to$ local Lagrange Finite-Element of order $p$

.footnote[P. Ciarlet. The finite element method for elliptic problems. SIAM, 2002.]

---

# DUNE LocalFunctions

.center[<img src="images/local-finite-element.png" width="60%">]

In Dune, a local finite element consists of three things:
1. A set of shape functions on a fixed reference element $T_\text{ref}$,
2. a mechanism to evaluate the degrees of freedom $\{\sigma_i\}_{i=0}^{N-1}$
   for any sufficiently smooth function $f$ on $T_\text{ref}$,
3. a map that assigns each degree of freedom to a subentity of $T_\text{ref}$,
   and gives it a local index on that same subentity.

---

# DUNE LocalFunctions
## LocalFiniteElement

The interface `LocalFiniteElement` consists of three functions to access the three components:

```c++
const LocalBasisType& localBasis() const;
const LocalInterpolationType& localInterpolation() const;
const LocalCoefficientsType& localCoefficients() const;
```

Additionally, some more information is exported:

```c++
GeometryType type() const;    // The shape of the T_ref element
unsigned int size() const;    // The number N of shape functions
```

### Associated types
The types `LocalBasisType`, `LocalInterpolationType`, and `LocalCoefficientsType` are available in a
`Traits` type provided by each implementation of the `LocalFiniteElement` interface.

---

# DUNE LocalFunctions
## LocalFiniteElement - Example

An order $p$ Lagrange local finite-element on simplices can be constructed directly:

```c++
#include <dune/localfunctions/lagrange/lagrangesimplex.hh>

template <class D, class R, int d, int p>
class LagrangeSimplexLocalFiniteElement;
```

with `D,R` field-types of the domain/range (e.g. `double`), dimension `d` of the domain:

```c++
using P1_2d = LagrangeSimplexLocalFiniteElement<double,double,2,1>;
```

### Access the members

```c++
P1_2d localFE;
auto const& localBasis = localFE.localBasis();
auto const& localInter = localFE.localInterpolation();
auto const& localCoeff = localFE.localCoefficients();
```

---

# DUNE LocalFunctions
## LocalFiniteElement - Example

An order $p$ Lagrange local finite-element on simplices can be constructed directly:

```c++
#include <dune/localfunctions/lagrange/lagrangesimplex.hh>

template <class D, class R, int d, int p>
class LagrangeSimplexLocalFiniteElement;
```

with `D,R` field-types of the domain/range (e.g. `double`), dimension `d` of the domain:

```c++
using P1_2d = LagrangeSimplexLocalFiniteElement<double,double,2,1>;
```

### Associated types

```c++
using P1LocaBasis = typename P1_2d::Traits::LocalBasisType;
using P1LocaInter = typename P1_2d::Traits::LocalInterpolationType;
using P1LocaCoeff = typename P1_2d::Traits::LocalCoefficientsType;
```

---

# DUNE LocalFunctions
## LocalBasis

The shape functions $\phi_i: T_\text{ref}\to\mathbb{R}^m$ of the finite-element are implemented in classes
that follow the `LocalBasis` interface.

- Shape functions can be evaluated in coordinates of $T_\text{ref}$ (all at once)
- Derivatives w.r.t. to local coordinates can be evaluated (all at once)

### Interface:
```c++
void evaluateFunction (const Traits::DomainType& local,
                       std::vector<Traits::RangeType>& out) const;
void evaluateJacobian (const Traits::DomainType& local,
                       std::vector<Traits::JacobianType>& out) const;
```

Corresponding types are defined in `Traits` type of the local basis implementation class.

---

# DUNE LocalFunctions
## LocalBasis - Example

Evaluation of a function given as tuple of local coefficients: $f(\hat{x}) = \sum_i f_i \phi_i(\hat{x})$ for $\hat{x}\in T_\text{ref}$

```c++
auto const& localBasis = localFE.localBasis();
using LocalBasis = std::decay_t<decltype(localBasis)>;

using RangeType = typename LocalBasis::Traits::RangeType;
std::vector<RangeType> phi;              // {phi_i(local)}_i
localBasis.evaluateFunction(local, phi); // local coordinate in T_ref

RangeType f_x(0);
for (std::size_t i = 0; i < phi.size(); ++i)
  f_x.axpy(f_i[i], phi[i]);  // f_x += f_i[i] * phi[i]
```

---

# DUNE LocalFunctions
## LocalInterpolation

For sufficiently smooth functions $f$ defined on the reference
element, compute the values $f_0:=\sigma_0(f), \ldots, f_{N-1} :=\sigma_{N-1}(f)$ of Ciarlet's definition

```c++
template <class F, class R>
void interpolate (const F& f, std::vector<R>& out) const;
```
with $f:T_\text{ref}\to R$

## Example

Given the function $f(\hat{x})$ we want to compute its Lagrange coefficients:

```c++
auto const& localInter = localFE.localInterpolation();

std::vector<double> f_i;
localInter.interpolate(f, f_i);
```

---

# DUNE LocalFunctions
## LocalCoefficients

.center[<img src="images/P1.png" width="50%">]

- Local basis functions should be combined/associated on multiple elements

---

# DUNE LocalFunctions
## LocalCoefficients

.center[<img src="images/local-coefficients.png" width="60%">]

- Degrees of freedom are associated to element subentities + a local index
- The association is encoded in a `LocalKey`
- Triples of numbers `(s, c, i)` denote the number `s` of the associated subentity, the codimension `c` of that subentity, and the number `i` of the degree of freedom within that subentity.

---

# DUNE LocalFunctions
## LocalCoefficients

The interface functions:

```c++
const LocalKey& localKey (std::size_t i) const;
std::size_t size () const;
```

where `LocalKey` is a class with the following interface

```c++
unsigned int subEntity () const;  // s
unsigned int codim () const;      // c
unsigned int index () const;      // i
```

---

# DUNE LocalFunctions
## LocalCoefficients - Example
Store the coefficients in a global vector of a grid, using grid entities

```c++
std::vector<double> coefficients(gridView.size(dim)); // global coeff vector
// fill coefficients -> see exercise

auto const& indexSet = gridView.indexSet();
for (auto const& cell : elements(gridView))
{
  // ...
  auto const& localCoeff = localFE.localCoefficients();
  // ...

RangeType f_x(0);
  for (std::size_t i = 0; i < phi.size(); ++i) {
    auto const& key = localCoeff.localKey(i);
    auto idx = indexSet.subIndex(cell, key.subEntity(), key.codim()); // global vertex index
    f_x.axpy(coefficients[idx], phi[i]);  // f_x += coefficients[i] * phi[i]
  }
}
```

---

# DUNE LocalFunctions
## Global numbering of degrees of freedom

- If we have an association of each local DOF to a subentity of the corresponding element,
this information can be combined using the `MCMGMapper` to define global unique indices for each local degree of freedom / shape function.

- **Problem 1:** What if two or more DOFs are associated to a single grid entity?
- **Problem 2:** How to identify the association from local-key information to the corresponding local-key on the neighbouring element?

### Possible solution
- Use *twist-free* grids, e.g. YaspGrid, ALUGrid: view on a subentity from inner and outer element is the same (e.g. same orientation of an edge)
- Transform/reorder the local indices on a subentity to match the indexinv on the neighbouring element $\to$ Twist-Utilities (not standardized)

---

# DUNE LocalFunctions
## Global numbering of degrees of freedom

---
layout: true
class: tud-dark, typo
---
class: center, middle, tud-dark
# Exercise 4
---

# Exercise 4

1. Given a global function $f:\Omega\to\mathbb{R}$
2. On a simplex grid, define local finite-element spaces
3. Define a global finite-element space by combinding the local finite-elements and associating to
   each degree of freedom a global index.
3. Interpolate the global function $f$ into this global finite-element space $\to f_h$.
4. Test your approximation by computing the integral of $f$ and $f_h$ over the domain $\Omega$

---

# Exercise 4
## Concrete Setup

- You domain $\Omega=[0,2]\times[0,2]$ is triangulated with subdivided boxes, 8 in each direction
- The local finite-element space should represent P2 elements
- The local degrees of freedom are associated to the vertices and edges of your grid (one for each entity)
- Interpolate the function

\[
   f(x) = \sin(x_0) \cdot \cos(x_1)
\]

---
layout: true
class: tud-light, typo
---
class: center, middle
# DUNE Functions
---

# DUNE Functions

- Association of local degrees-of-freedom to global indices
- The `LocalBasis` is a "local-view" onto the `GlobalBasis`
- The same local basis-functions be combined in multiple ways to define global basis-functions

## Examples:
- `LagrangeBasis`
- `LagrangeDGBasis`
- `RaviartThomasBasis`

---

# DUNE Functions
## DefaultGlobalBasis

Interface for global basis functions:

```c++
template <class Impl>
class DefaultGlobalBasis;
```

with interface functions

```c++
const GridView& gridView () const;  // get the associated GridView
void update (const GridView & gv);  // update the indices after grid has changed
size_type dimension () const;       // total number of basis functions
LocalView localView () const;       // get the local-view onto the global basis
```

---

# DUNE Functions
## DefaultGlobalBasis - Example

Lagrange basis functions:

```c++
#include <dune/functions/functionspacebases/lagrangebasis.hh>

namespace Dune::Functions
{
  template <class GridView, int p, class R = double>
  class LagrangeBasis;
}
```

can be constructed from `GridView`.

**Alternative:** use factory functions
```c++
using namespace Dune::Functions::BasisFactory;
auto basis = makeBasis(gridView, lagrange<p>());
```

---

# DUNE Functions
## DefaultLocalView

- View onto a global basis on a single element
- Must be bound to an element before it can be used
- Map local indices of DOFs to global indices

.center[<img src="images/P1.png" width="40%">]

---

# DUNE Functions
## DefaultLocalView - Interface

Binding the `LocalView` to an element:

```c++
void bind (const Element& e);     // Bind the view to a grid element
void unbind ();                   // Release bound element
bool isBound () const;            // Check whether LocalView is bound
const Element& element() const;   // Return the bound element.
```

Mapping indices:

```c++
size_type size () const;     // Number of DOFs on the current element
size_type maxSize () const;  // Maximum local size for any element on the GridView

MultiIndex index (size_type i) const; // map local index i to global Multi-Index
```

Local finite-element:
```c++
auto const& localFE = basis.localView().tree().finiteElement(); // later more...
```

---

# DUNE Functions
## Example - Same as in exercise 4
```c++
auto basis = makeBasis(gridView, lagrange<2>());
std::vector<double> coefficients(basis.dimension());

auto localView = basis.localView();
for (auto const& cell : elements(gridView)) {
  auto geo = cell.geometry();

localView.bind(cell);
  auto const& localFE = localView.tree().finiteElement();

// local-function that can be evaluated in local coordinates
  auto localF = [&](auto const& local) { return f(geo.global(local)); };

std::vector<double> localCoefficients;
  localFE.localInterpolation().interpolate(localF, localCoefficients);

for (std::size_t i = 0; i < localFE.size(); ++i)
    coefficients[localView.index(i)] = localCoefficients[i];
}
```

---

# DUNE Functions
## Interpolation

- Interpolate global function into global function-space
- Utility function:

```c++
auto basis = makeBasis(gridView, lagrange<2>());
std::vector<double> coefficients(basis.dimension());

Dune::Functions::interpolate(basis, coefficients, f);
```

> **AMDiS:** `interpolate(basis, coefficients, f)`

---

# DUNE Functions
## Basis Trees

- More complicated bases than Lagrange require composition, e.g., multiple components of the same or
  mixed types of function spaces.
- Example: Taylor-Hood basis consisting of velocity and pressure components
<br><br> <img align="right" width="340" height="200" src="images/basistree.png">
$ \qquad \qquad V_{\text{TH}} = \underbrace{(V_x \otimes V_y \otimes V_z)}_{=: V_v} \otimes V_p $
  - composite (different components) $ V_{\text{TH}} = V_v \otimes V_p $
  - power (identical components) $ V_v = V_x^3 $
  - leaf (implementation)

- Tree nodes are indexed starting from the root node with tuple of integers

---

# DUNE Functions
# Example: Taylor-Hood-Basis in 2d

- composition of velocity space and pressure space
  <br>.icon-rechts[] `composite(<velocity>, <pressure>)`

--
- velocity space itself is a composition of 2 scalar components
  <br>.icon-rechts[] `power<2>(<velocity_component>)`

--
- we choose nodal basis of order 2 for velocity and order 1 for pressure
  <br>.icon-rechts[] `lagrange<order>()` or `lagrange(order)`

```c++
using namespace Dune::Functions::BasisFactory;
auto factory = composite(
                 power<2>(
                   lagrange<2>()),
                 lagrange<1>());

auto V_TH = makeBasis(gridView, factory);
```

---

# DUNE Functions
## Basis Trees

- Accessing the tree in the `DefaultLocalView`:

```c++
// Return the local ansatz tree associated to the bound entity
const Tree& tree() const;
```

- All nodes in the tree allow to access its childs:

```c++
auto& child (std::size_t i) [const];            // for power nodes
auto& child (Dune::index_constant<i>) [const];  // for composite nodes
// or
Dune::TypeTree::child (<tree>, indices...);
```

- composite Nodes require integer constants as indices!
  <br>.icon-rechts[] shorthand `Dune::Indices::_i`

---

## DUNE Functions
## Basis trees

- Interface for all basis trees/nodes:

```c++
// map index of a local basis of node into the index of all local bases on the tree
size_type localIndex(size_type i) const;
// Number of all local basis functions in this node
size_type size() const;
// Get a unique index of the node in the tree
size_type treeIndex() const;
// The element the node/local view is bound to
const Element& element() const;
```

- If *tree* is a leaf-node, it might provide an additional interface function:

```c++
// The local finite-element
const FiniteElement& finiteElement() const
```

---

## DUNE Functions
## Basis trees

The `localIndex(i)` function can be used to arrange all the leaf node local-basis indices into all the
basis functions on the element:

```c++
localView.bind(cell);                               // assume cell == (simplex,2d)
auto const& tree = localView.tree();                // Taylor-Hood
auto const& node_v = tree.child(Dune::Indices::_0); // node_v.size() == 2x6 == 12
auto const& node_p = tree.child(Dune::Indices::_1); // node_p.size() == 3
...
for (std::size_t i = 0; i < node_p.finiteElement().size(); ++i)
  coefficients[localView.index(node_p.localIndex(i))] = localCoefficients_p[i];
```

---

# DUNE Functions
## Basis Trees

- Sub-bases can be constructed using the tree accessor by node indexing.
- **Example:** view onto subspace basis

```c++
auto V_v = subspaceBasis(V_TH, Dune::Indices::_0);
auto V_x = subspaceBasis(V_TH, Dune::Indices::_0,0);
auto V_y = subspaceBasis(V_TH, Dune::Indices::_0,1);
auto V_z = subspaceBasis(V_TH, Dune::Indices::_0,2);

auto V_p = subspaceBasis(V_TH, Dune::Indices::_1);
```

---

# DUNE Functions
## Basis Trees

- Sub-bases can be constructed using the tree accessor by node indexing.
- **Example:** view onto subspace basis

```c++
auto V_v = subspaceBasis(V_TH, Dune::Indices::_0);
auto V_x = subspaceBasis(V_v, 0); // alternative
auto V_y = subspaceBasis(V_v, 1);
auto V_z = subspaceBasis(V_v, 2);

auto V_p = subspaceBasis(V_TH, Dune::Indices::_1);
```

Access the corresponding node directly as tree
```c++
auto localView = V_p.localView();
...
localView.bind(cell);                               // assume cell == (simplex,2d)
auto const& tree = localView.tree();                // pressure node, tree.size() == 3
...
for (std::size_t i = 0; i < tree.finiteElement().size(); ++i)
  coefficients[localView.index(tree.localIndex(i))] = localCoefficients_p[i];
```

---

# DUNE Functions
## Index merging strategies

- `PowerNode` and `CompositeNode` combines indices of sub-bases into global contiguous indexing
  scheme
- Let $\{\phi_i\}$ be the set of basis functions of $V_h$. A function $u_h\in V_h$ can be
represented as

\[
  u_h(x) = \sum_i u_i \phi_i(x)
\]

with coefficients $\{u_i\}$.

---

# DUNE Functions
## Index merging strategies

- Assume we have a product space $ V = U_1\otimes U_2 $
- How to number all the basis functions in this product space?

\[
  u_1(x) = \sum_{i\in I_1} u_i \phi^1_i(x) \\
  u_2(x) = \sum_{i\in I_2} u_i \phi^2_i(x) \\
\]

with $u = (u_1, u_2) $ and $I_1, I_2$ subsets of the global index set.

---

# DUNE Functions
## Index merging strategies

```c++
auto factory = power<2>(lagrange<1>(), flatLexicographic());
// with lexicographic index-merging strategy
// I1 = (0,1,2,3), I2 = (4,5,6,7)
```

.center[
<img src="images/index-merging-1.png" height="350">
]

.footnote1[(1) Image from: S. Müthing, *A Flexible Framework for Multi Physics and Multi Domain PDE Simulations*, 2015]

---

# DUNE Functions
## Index merging strategies

```c++
auto factory = power<2>(lagrange<1>(), flatInterleaved());
// with interleaved index-merging strategy
// I1 = (0,2,4,6), I2 = (1,3,5,7)
```

.center[
<img src="images/index-merging-1.png" height="350">
]

.footnote1[(1) Image from: S. Müthing, *A Flexible Framework for Multi Physics and Multi Domain PDE Simulations*, 2015]

---

# DUNE Functions
## Blocked Indexing

- Instead of creating one contiguous global index set, one can create multiple indexsets, e.g., for each component a
  separate indexing
- This allows to use blocked/hierarchic data structures
- The `LocalView` returns not a single index but a `MultiIndex` that addresses entries in blocked containers

.center[
<img src="images/blocking.png" height="250">
]

---

# DUNE Functions
## Blocked Indexing

In the example from above, if we use `blockedLexicgraphic()` index-merging strategy

```c++
auto factory = power<2>(lagrange<1>(), blockedLexicgraphic());
```

we get the numbering:

\[
  u_1(x) = \sum_{i\in I} u_{0,i} \phi_i(x) \\
  u_2(x) = \sum_{i\in I} u_{1,i} \phi_i(x) \\
\]

with coefficients indexed by *MultiIndices* $(0,i)$ and $(1,i)$ for $u_1$ and $u_2$ respectively.

```c++
std::array<std::vector<double>, 2> U; // data structure
```

---

# DUNE Functions
## Blocked Indexing

In the example from above, if we use `blockedInterleaved()` index-merging strategy

```c++
auto factory = power<2>(lagrange<1>(), blockedInterleaved());
```

we get the numbering:

\[
  u_1(x) = \sum_{i\in I} u_{i,0} \phi_i(x) \\
  u_2(x) = \sum_{i\in I} u_{i,1} \phi_i(x) \\
\]

with coefficients indexed by *MultiIndices* $(i,0)$ and $(i,1)$ for $u_1$ and $u_2$ respectively.

```c++
std::vector<FieldVector<double, 2>> U; // data structure
```

---

# DUNE Functions
## DiscreteGlobalBasisFunction

- Combination of coefficient vector and basis
- Provides evaluation function

```c++
auto fh = Dune::Functions::makeDiscreteGlobalBasisFunction<Range>(basis, coefficients);

auto localFh = localFunction(fh);
for (auto const& cell : elements(gridView)) {
  localFh.bind(cell);
  auto fhAtLocal = localFh(local);
}
```

> **AMDiS:** `DiscreteFunction{coefficients, basis}`

---

# DUNE Functions
## DiscreteGlobalBasisFunction

- Use `subspaceBasis` to represent a component of your discrete function:

```c++
using namespace Dune::Indices;
auto V_v = subspaceBasis(V_TH, _0);
auto V_p = subspaceBasis(V_TH, _1);

using VelocityRange = FieldVector<double,2>;
using PressureRange = double;

// discrete function to evaluate on a grid
auto vh = makeDiscreteGlobalBasisFunction<VelocityRange>(V_v, coefficients);
auto ph = makeDiscreteGlobalBasisFunction<PressureRange>(V_p, coefficients);

// function interpolation
interpolate(V_v, coefficients, [](auto const& x) { return VelocityRange{1.,1.}; });
interpolate(V_p, coefficients, [](auto const& x) { return 42; });
```

---
layout: true
class: tud-dark, typo
---
class: center, middle, tud-dark
# Exercise 5
---

# Exercise 5

1. Let $\Omega=[0,2]^2$ be a square domain triangulated with simplices.
2. Construct a global basis $V=P_2\otimes P_2$ with blocked-interleaved index merging strategy
3. Define a vector container for the coefficients suitable for this basis
3. Interpolate the global function $f$ into this global finite-element space $V \to f_h$.
4. Compute the integral $\int_{\Omega}\operatorname{div}(f_h)\,\text{d}x$ over the domain $\Omega$

\[
  f(x) = \big(sin(x_0) cos(x_1) + 3x_0x_1^2, -cos(x_1) sin(x_1) - x_1^3\big)^T
\]

---
layout: true
class: tud-light, typo
---
class: center, middle
# AMDiS
## Adaptive Multi-Dimensional Simulations
---

# AMDiS - Poisson equation

```c++
#include <amdis/AMDiS.hpp>                              // Include essential headers
#include <amdis/LocalOperators.hpp>
int main(int argc, char** argv)
{
  using namespace AMDiS;                                // Namespace of all AMDiS functions
  using namespace Dune::Functions::BasisFactory;        // Use DUNE facilities directly

Environment env{argc, argv};                          // Initialize Initfile, MPI...

Dune::YaspGrid<2> grid{ {1.0,1.0}, {8,8} };           // Grid: Unit square domain
  ProblemStat prob{"example", grid, lagrange<2>()};     // Lagrange elements of deg. 2
  prob.initialize(INIT_ALL);                            // Prepare data structures

auto c = [](auto x) { ... };
  prob.addMatrixOperator(sot(1.0));                     // a(v,u) = <grad(v), grad(u)>
  prob.addVectorOperator(zot(c, 6));                    // b(v) = <v, c(x)>
  prob.addDirichletBC([&](auto x) { return b(x); },     // boundary
                      [&](auto x) { return g(x); });    // u = g | on boundary

AdaptInfo adaptInfo{"adapt"};
  prob.assemble(adaptInfo);                             // Assemble the linear system
  prob.solve(adaptInfo);                                // Solve the linear system
  prob.writeFiles(adaptInfo);                           // Write solution to file
}
```

---

# Installing AMDiS and its dependencies

Find install instructions in [README.md](https://gitlab.mn.tu-dresden.de/amdis/amdis-core/-/blob/master/README.md)
file or at https://amdis.readthedocs.io

1. Download DUNE modules: common, geometry, grid, localfunctions, istl, typetree, functions
2. Download AMDiS (into the same `<dune-base>` directory as all DUNE modules)
```bash
git clone https://gitlab.mn.tu-dresden.de/amdis/amdis-core.git
```
--
3. Configure and build everything
```bash
./dune-common/bin/dunecontrol cmake -DCMAKE_BUILD_TYPE=Release
./dune-common/bin/dunecontrol make -j10
```

---

# My first AMDiS project

Similar to DUNE, we provide a script `amdisproject` that helps you create a simple initial empty project:

```bash
./amdis-core/bin/amdisproject my-first-project
```

Some question are asked about name, version, maintainer, and dependencies. After answering everything,
a directory `my-first-project` is created in `<dune-base>` with a `CMakeLists.txt` file and some
example code.
--

```
- README.md                                     - doc/doxygen/CMakeLists.txt
- CMakeLists.txt                                - doc/doxygen/Doxylocal
- config.h.cmake                                - cmake/modules/CMakeLists.txt
- dune.module                                   - cmake/modules/MyFirstProjectMacros.cmake
- src/CMakeLists.txt                            - init/CMakeLists.txt
- src/my-first-project.cpp                      - init/my-first-project.dat
- amdis/my-first-project/CMakeLists.txt         - macro/CMakeLists.txt
- amdis/my-first-project/MyFirstProject.hpp     - macro/my-first-project.Xd.amc
```

---

# My first AMDiS project

The source file `src/my-first-project.cpp` contains the following code snippet:

```c++
#include <amdis/AMDiS.hpp>

using namespace AMDiS;
int main(int argc, char** argv)
{
  Environment env{argc, argv};

// your code comes here
}
```

--
- Main include header `<amdis/AMDiS.hpp>` includes *some* parts of the library.

--
- All classes and functions can be found in the namespace `AMDiS`

--
- `Environment` initializes some global states, like MPI and PETSc

--
- The `CMakeLists.txt` file in the `src/` folder contains the executable.

---

# My first AMDiS project

The directory `amdis/my-first-project/` can be used to place header files shared my several source
files in this project. An example is given in `MyFirstProject.hpp`:

```c++
#pragma once

namespace AMDiS::MyFirstProject
{
  // add your classes here
}
```
that can be included in the source file:

```c++
#include <amdis/AMDiS.hpp>
#include <amdis/my-first-project/MyFirstProject.hpp>

using namespace AMDiS;
int main(int argc, char** argv) { /* ... */ }
```

Configure, compile and link it against all dependencies, by using `dunecontrol`:
```bash
./dune-common/bin/dunecontrol --only=my-first-project all
```
---
layout: true
class: tud-light, typo
---
class: center, middle
# AMDiS - The `ProblemStat` class
---

# AMDiS - Variational problem

Find $(u_0, u_1)\in U_0\times U_1$ s.t.
\[
  A((v_0,v_1),(u_0,u_1)) = B((v_0,v_1))\quad\forall (v_0,v_1)\in U_0'\times U_1'
\]
with (bi-)linear form

\[
  A((v_0,v_1),(u_0,u_1)) = a_{00}(v_0,u_0) + a_{01}(v_0,u_1) + a_{10}(v_1,u_0) + a_{11}(v_1,u_1) + \\
  \bar{a}((v_0,v_1),(u_0,u_1)) \\
  B((v_0,v_1)) = b_0(v_0) + b_1(v_1) + \bar{b}((v_0,v_1))
\]

Here $a_{ij}, \bar{a}, b_i$, and $\bar{b}$ are "simple" components of the (bi-)linear form and are called *Operators* in the following.

---

# AMDiS - Variational problem

Often, the *Operators* have a common structure
\[
  a(v,u) = \int_{\Omega} \mathfrak{a}(c,v,u,\nabla v,\nabla u)\,\text{dx} = \sum_{T\in\mathcal{T}_h} \int_T \mathfrak{a}(c,v,u,\nabla v,\nabla u)\,\text{dx}
\]
with $c:\Omega\to\mathbb{R}$ a coefficient function. Sometimes we have vector-valued or matrix-valued coefficients. All these are defined with the generic $c$ expression.

We call the integrand $\mathfrak{a}$ the *OperatorTerm* associated to the *Operator* $a$. The operator is an operator term bound to the grid (view) it is assembled on.

## Examples:

- $\mathfrak{a}(c,v,u,\nabla v,\nabla u)=c(x) v(x) u(x)$ called `test_trial` or `zot`.
- $\mathfrak{a}(c,v,u,\nabla v,\nabla u)=c(x) \nabla v(x) \nabla u(x)$ called `gradtest_gradtrial` or `sot`.

---

# AMDiS - Variational problem

The *Problem* is described by
- A *grid* it is defined on
- The (composite) *basis* $U_0\times U_1$
- *Constraints* on the trial and test basis, e.g., Dirichlet constraints or periodic constraints
- An operator *matrix* $A$ and an operator *rhs*-vector $B$
  * Composed of individual operators, defined by operator terms
- A *solution* vector $U=(u_0,u_1)$
  * Might be decomposed into components

All these problem components + some additional management/visualization tools are combined into the
class `ProblemStat`.

---

# AMDiS - The `ProblemStat` class

**Main interface class:** The class template `ProblemStat` is parametrized with `BasisTraits`, i.e.,
a description of the grid and global basis.
```c++
template <class BasisTraits>
class ProblemStat;
```

Can be constructed in multiple ways
```c++
ProblemStat (std::string name);               // requires `BasisTraits`

template <class Grid>
ProblemStat (std::string name, Grid& grid);   // requires `BasisTraits`

template <class Grid, class Basis>
  requires Concepts::GlobalBasis<Basis>
ProblemStat (std::string name, Grid& grid, Basis& globalBasis);         // CTAD

template <class Grid, class BF>
  requires Concepts::PreBasisFactory<BF>
ProblemStat (std::string name, Grid& grid, BF const& basisFactory);     // CTAD
```

---

# AMDiS - The `ProblemStat` class
## Main Setup functions

Create the basis and the grid, initialize internal data structures for the solution vector and
the linear system, create adaption markers and file writers.

```c++
// Initialize the problem.
//
// Parameters read in initialize()
//   - `[GRID_NAME]->global refinements`:  nr of initial global refinements
//
void initialize (Flag initFlag, ProblemStat* adoptProblem = [...], Flag adoptFlag = [...]);

// Read the grid and solution from backup files and initialize the problem
//
// Parameters read in restore() for problem with name 'PROB'
//   - `[PROB]->restore->grid`:      name of the grid backup file
//   - `[PROB]->restore->solution`:  name of the solution backup file
//
void restore (Flag initFlag);
```

---

# AMDiS - The `ProblemStat` class
## Access to grid and basis

After initialization you can access the internally stored grid and global basis:

```c++
// Return the stored grid
std::shared_ptr<Grid> grid ();

// Return the stored basis
std::shared_ptr<GlobalBasis> globalBasis ();

// Return the gridView of the basis
GridView gridView () const;
```

---

# AMDiS - The `ProblemStat` class
## Access to data structures

The `ProblemStat` internally stores the data structures for the linear system to assemble and for the
solution vector. All these data structure can be accessed directly:

```c++
// Return the system-matrix -> BiLinearForm<GlobalBasis, GlobalBasis, [...]>
std::shared_ptr<SystemMatrix> systemMatrix();

// Return the system rhs vector -> LinearForm<GlobalBasis, [...]>
std::shared_ptr<SystemVector> rhsVector();

// Return the solution vector -> DOFVector<GlobalBasis, [...]>
std::shared_ptr<SolutionVector> solutionVector();
```

The solution can also be extracted as *DiscreteFunction* $\to$ see later
```c++
// Return a discrete function wrapper representing a (sub-)view onto the solution
template <class Range = void, class... Indices>
auto solution (Indices... ii);
```

---

# AMDiS - The `ProblemStat` class
## Describe your PDE

(Bi-)linear form is composed of smaller components, called *OperatorTerms*.

```c++
// Add an operator term to the bilinear form
template <class OperatorTerm, class RowTreePath = [...], class ColTreePath = [...]>
void addMatrixOperator (OperatorTerm const& op, RowTreePath row = {}, ColTreePath col = {})

// Add an operator term to the linear form
template <class OperatorTerm, class TreePath = [...]>
void addVectorOperator (OperatorTerm const& op, TreePath path = {})
```

If terms are to be evaluated on (part of) the boundary of the grid only:

```c++
// Add an operator term to the bilinear form
template <class OperatorTerm, class RTP = [...], class CTP = [...]>
void addMatrixOperator (BoundaryType b, OperatorTerm const& op, RTP row = {}, CTP col = {})

// Add an operator term to the linear form
template <class OperatorTerm, class TP = [...]>
void addVectorOperator (BoundaryType b, OperatorTerm const& op, TP path = {})
```

---

# AMDiS - The `ProblemStat` class
## Describe your PDE

Constraints can be added by

```c++
// Enforce Dirichlet boundary values for the solution vector on boundary regions
template <class Boundary, [class RowTreePath, class ColTreePath,] class Values>
void addDirichletBC (Boundary const& predicate,
                     [RowTreePath row, ColTreePath col,]
                     Values const& values);

// Add a periodic boundary conditions to the system, by specifying a face transformation
// y = A*x + b of coordinates. We assume, that A is orthonormal.
void addPeriodicBC (BoundaryType id, WorldMatrix const& A, WorldVector const& b);

// General constraints specification
void addConstraint (BoundaryCondition<SystemMatrix,SolutionVector,SystemVector> constraint);
```

---

# AMDiS - The `ProblemStat` class
## Adaptation and solving the problem

The following steps can be called to solve the PDE system:
- `markElements`: use added markers to mark elements for coarsen/refine
- `adaptGrid`: coarsen/refine the grid
- `assemble`: assemble the stiffness matrix and right hand side
- `solve`: solve the linear system above
- `writeFiles`: write files to output

Default implementations handle most cases
- `ProblemStat` interface allows to customize markers, interpolation during refinement, file writers,
  linear solvers [, assemblers],...
- Many solvers and file writers supported out of the box

---

# AMDiS - The `ProblemStat` class
## Customization

User can inherit from `ProblemStat` and replace the implementations freely, e.g. override the `assemble()` method
- adaptation subroutines can be called individually, so a user can stop after a certain step, run
  some user code and then resume with the default AMDiS implementations
- almost all member functions can be accessed / overridden to allow maximum flexibility

---
layout: true
class: tud-dark, typo
---
class: center, middle, tud-dark
# Exercise 6
---

# Exercise 6

- Create a new amdis module `amdis-exercise` using the tool `amdisproject`
- Copy the initial example of the Poisson equation in the file `src/amdis-exercise.cc`
- Implement the missing functions with some meaningful content.
- Configure, compile and run the code.

---

# Exercise 6

- See the default init-file in the project folder, adapt its content to allow configuration
   of the Poisson problem. See https://amdis.readthedocs.io/en/latest/reference/Initfile/ for
   some infos about the parameter file
- Configure the automatic output of the solution of the problem
- Run your code again with

`./build-cmake/src/amdis-exercise init/amdis-exercise.dat`

- Visualize the result in ParaView
---
layout: true
class: tud-light, typo
---
class: center, middle
# DOFVector and DiscreteFunction
---

# DOFVector and DiscreteFunction

- A `DOFVector` is a container for storing the coefficients of a function $ u_h\in V_h $
- It is attached to a global basis to give its coefficients a meaning.

--
- A `DiscreteFunction` is a "GridFunction" representing $ u_h $

Let $\{\phi_i\}$ be the set of basis functions of $V_h$. A function $u_h\in V_h$ can be
represented as

\[
  u_h(x) = \sum_i u_i \phi_i(x)
\]

with coefficients $\{u_i\}$.

The pair $ (\{u_i\},\{\phi_i\}) $ is called `DOFVector` and the function $ u_h(x) $ is called `DiscreteFunction`.
(see [documentation](https://amdis.readthedocs.io/en/latest/reference/DOFVector))

---

# DOFVector and DiscreteFunction

- Assume we have a product space $ V = V_1\times V_2 $

--
- Coefficients $\{u_i\}$ in a `DOFVector` `U` might be stored in one *flat* vector container .icon-rechts[] some coefficients
  correspond to basis functions of $V_1=span(\{\phi_i^1\})$ the others to basis functions of $V_2=span(\{\phi_i^2\})$

\[
  u_1(x) = \sum_{i\in I_1} u_i \phi^1_i(x) \\
  u_2(x) = \sum_{i\in I_2} u_i \phi^2_i(x) \\
\]

with $u = (u_1, u_2) $.

---

# DOFVector

```c++
DOFVector U{gridView, power<2>(lagrange<1>())};
```

---

# DOFVector

```c++
DOFVector U{gridView, power<2>(lagrange<1>(), flatLexicographic())};
// with lexicographic index-merging strategy
// I1 = (0,1,2,3), I2 = (4,5,6,7)
```

.center[
<img src="images/index-merging-1.png" height="350">
]

.footnote1[(1) Image from: S. Müthing, *A Flexible Framework for Multi Physics and Multi Domain PDE Simulations*, 2015]

---

# DOFVector

```c++
DOFVector U{gridView, power<2>(lagrange<1>(), flatInterleaved())};
// with interleaved index-merging strategy
// I1 = (0,2,4,6), I2 = (1,3,5,7)
```

.center[
<img src="images/index-merging-1.png" height="350">
]

.footnote1[(1) Image from: S. Müthing, *A Flexible Framework for Multi Physics and Multi Domain PDE Simulations*, 2015]

---

# DOFVector

The class is parametrized:
```c++
template <class Basis, class ValueType = double, class TraitsType = BackendTraits<Basis>>
class DOFVector;
```

Can be constructed with
```c++
// Provide a full global basis
DOFVector::DOFVector(Basis& basis,
          DataTransferOperation op = DataTransferOperation::INTERPOLATE);

// construct directly from a basis-factory
template <class BF>
  requires Concept::PreBasisFactor<BF>
DOFVector::DOFVector(GridView const& gridView, BF basisFactory,
          DataTransferOperation op = DataTransferOperation::INTERPOLATE);
```

Extract the type of a DOFVector:

```c++
using DOFVectorType = decltype(DOFVector{std::declval<GridView>(), power<2>(lagrange<1>())})
```

---

# DOFVector

Interface functions:
```c++
/// Reference to the global basis
Basis const& basis () const;

/// Coefficient vector
Coefficients const& coefficients () const; // -> see VectorBackend
Coefficients&       coefficients ();

/// Write DOFVector to file
void backup (std::string const& filename);

/// Read backup data from file
void restore (std::string const& filename);

/// Resize the interface coefficient vector to basis dimension
void resize ();
void resizeZero ();
```

---

# DOFVector and DiscreteFunction

Transform the `DOFVector` into a `DiscreteFunction`
```c++
template <class Range = void, class DV, class... Indices>
auto valueOf (DV& dofVec, Indices... ii) -> DiscreteFunction<...>;
```

## Example
```c++
DOFVector U{gridView, power<2>(lagrange<1>())};

// getting the composite discrete-function u
auto u = valueOf(U);                          // automatic range type
auto V = valueOf<FieldVector<double,2>>(U);   // specify a range type
DiscreteFunction v{U};

// getting the component functions u_i
auto u_0 = valueOf(U,0);
auto u_1 = valueOf(U,1);       // NOTE: u_i == u.child(i)
DiscreteFunction v_1{U,1};
```

---

# DiscreteFunction

Interface functions for (global) DiscreteFunction:

```c++
/// Evaluate DiscreteFunction in global coordinates. NOTE: expensive
Range operator() (Domain const& x) const;

/// Return treePath associated with this discrete function
TreePath const& treePath () const;

/// Return global basis bound to the DOFVector
Basis const& basis () const;

/// Return const coefficient vector
Coefficients const& coefficients () const;

/// Get a sub-view on the discrete function related to the extended tree path
template <class Range = void, class... Indices>
auto child (Indices... ii) const
```

---

# DiscreteFunction

A DiscreteFunction can be localized to an element, to allow evaluation in element-local coordinates:
```c++
auto localFunction (DiscreteFunction<...>);
```

This element-function has the following interface:
```c++
/// Bind the discrete function to the element
void bind (Element const& element);

/// Unbind from the element (optional, cleanup)
void unbind ();

/// Evaluate LocalFunction at bound element in local coordinates
Range operator() (const Domain& local) const  [[expects: isBound() ]];

/// Polynomial order of the discrete function (if available)
int order () const [[expects: isBound() ]];

/// Return the bound element
Element const& localContext () const [[expects: isBound() ]];
```

---

# DiscreteFunction
## Derivatives

A local basis provides derivatives w.r.t. local coordinates. Those can be transformed into global coordinate derivatives using the element geometry. Thus, a discrete function can easily be differentiated:

```c++
auto derivativeOf (DiscreteFunction<...>, tag::derivative_type);

// some specific derivative types:

auto gradientOf (DiscreteFunction<...>);    // scalar valued discrete function
auto divergenceOf (DiscreteFunction<...>);  // vector valued discrete function

// partial derivative of scalar valued discrete function in direction i
auto partialDerivativeOf (DiscreteFunction<...>, std::size_t i);
```

The derivatives are "GridFunctions" (see later) similar to `DiscreteFunction`.

---

# DOFVector and DiscreteFunction
## The solution DOFVector

The `ProblemStat` returns both, the `DOFVector` storing the solution and its `DiscreteFunction`:

```c++
ProblemStat prob{"prob", grid, power<2>(lagrange<1>())};

auto& U = *prob.solutionVector();           // NOTE: shared_ptr returned
auto u = prob.solution();

DiscreteFunction u_i = prob.solution(i);    // provide treepath argument
```

---
layout: true
class: tud-light, typo
---
class: center, middle
# GridFunctions
---

# GridFunctions

```c++
auto f = valueOf(F);

FieldVector<double,2> x{0.2, 0.3};  // global coordinate
auto f_x = f(x);                    // evaluate f at x
```

---

# GridFunctions

```c++
auto f = valueOf(F);
auto lf = localFunction(f);

for (auto const& e : elements(gridView)) {
  lf.bind(e);

FieldVector<double,2> a{0.2, 0.3};  // local coordinate
  auto lf_a = lf(a);                  // evaluate lf at a
}
```

---

# GridFunctions

<div class="pure-g">
  <div class="pure-u-3-4"><ul>
    <li>A <a href="https://amdis.readthedocs.io/en/latest/reference/GridFunctions">GridFunction</a> <code class="remark-inline-code">f</code> is a function-object that can be restricted to a grid element</li>
    <li>Restriction is called <strong>LocalFunction</strong> <code class="remark-inline-code">lf=localFunction(f)</code> and must be bound to the element, i.e., <code class="remark-inline-code">lf.bind(element)</code></li>
    <li>GridFunctions are built from <a href="https://amdis.readthedocs.io/en/latest/reference/GridFunctions/#list-of-grid-functions-and-expressions">Expressions</a> <code class="remark-inline-code">expr</code> that are a composition of some elementary terms, i.e., <code class="remark-inline-code">f=gridFunction(expr, gridView)</code></li>
  </ul></div>
  <div class="pure-u-1-4"><img src="images/macro_lf.png"></div>
</div>

```c++
auto expr1 = prob.solution(0);                          // discrete-function
auto expr2 = 1.0;                                       // constant
auto expr3 = std::ref(expr2);                           // reference
auto expr4 = [](auto x) { return x[0] + x[1]; };        // functor in global coords
// composition...
auto expr5 = prob.solution(0) + 1.0;                    // NOTE: one of the arguments
auto expr6 = max(evalAtQP(expr4), expr3);               //       must be GridFunction
auto no_expr = expr4 + 1.0;                             // ERROR
```

---

# Example of a composite GridFunction

Let $c$ be a (scalar) GridFunction

\[
  u(x) = \frac{1}{2}\left(c(x)^2(1-c(x))^2 + \frac{1}{\epsilon}\|\nabla c(x)\|^2\right)
\]

This can be constructed as composite GridFunction as follows:

```c++
auto c = valueOf(C); // e.g. c is discrete function
auto u = 0.5*(pow<2>(c * (1 - c)) + (1.0/eps) * unary_dot(gradientOf(c)));
```

or, if you have a more complicated expression, just put it into a lambda function:

```c++
auto v = invokeAtQP([eps](auto c_x, auto grad_c_x) {
  return 0.5*(std::pow(c_x * (1 - c_x), 2) + (1.0/eps) * unary_dot(grad_c_x));
}, c, gradientOf(c));
```

---

# Application of GridFunctions

- Use GridFunctions to build Operators: (see later)

```c++
auto opB = operatorTerm(tag::BiLinearForm, Expression);
prob.addMatrixOperator(opB, Row, Col);

auto opL = operatorTerm(tag::LinearForm, Expression);
prob.addVectorOperator(opL, Row);
```

--
- Use GridFunctions in BoundaryConditions:

```c++
prob.addDirichletBC(Predicate, Row, Col, Expression);
```

--
- Interpolate a GridFunction to a DOFVector:

```c++
prob.solution(Path).interpolate(Expression);
prob.solution(Path) << Expression;
```

---

# Application of GridFunctions

- Integrate a GridFunction on a GridView:

```c++
auto value = integrate(Expression, gridView);
```

--
- Define a refinement marker

```c++
auto marker = GridFunctionMarker{"marker", Grid, Expression};
```

---
layout: true
class: tud-dark, typo
---
class: center, middle, tud-dark
# Exercise 7
---

# Exercise 7

This exercise is about examining the abilities of DOFVectors, DiscreteFunctions and Expressions

1) Create a structured grid of $\Omega=[-2,2]^2$

2) Create a `DOFVector` of a scalar Lagrange P2 basis.

3) Interpolate the function $c(x)$ into the `DOFVector`:

\[
  c(x) = 0.5(1 - \tanh(3\rho(x)/\epsilon)),\text{ with }\rho(x)=\|x\| - 1
\]

for $\epsilon=0.1$

4) Integrate the following expression over the domain $\Omega$

\[
  \frac{1}{2}\left(c(x)^2(1-c(x))^2 + \frac{1}{\epsilon}\|\nabla c(x)\|^2\right)
\]

---

# Exercise 7
## Advanced exercise

- Replace the scalar basis with a power-basis
- Interpolate $c(x)$ into the first component of that DOFVector
- Interpolate $(c(x),0)$ into the full DOFVector
- Compute the divergence of the expression and the partial derivative  $\partial_0$ of the first
  component of the expression. Print the corresponding integral value.

---
layout: true
class: tud-light, typo
---
class: center, middle
# Adaptivity
---

# Adaptivity - Grid Refinement

- Grids are always hierarchical grids
- Gids allow multiple different forms of refinement:
  * bisection (longest edge, first edge, ...)
  * red refinement
  * blue refinement, trisection, multiple closures

---

# Adaptivity - Grid Refinement

- Some grid manager allow only global refinement:
  * structured grids: `YaspGrid`, `SPGrid`
  * special grids: `FoamGrid`

- Other grid managers support only a specific type of elements:
  * simplices: `AlbertaGrid`, `FoamGrid`
  * cubes: `YaspGrid`, `SPGrid`

- Red refinement with multiple closures:
  * Red-Green: `UGGrid`
  * None: `ALUGrid` (cubes), `UGGrid`

- Different refinement strategies:
  * bisection: `AlbertaGrid` (restricted), `ALUGrid` (always possible)
  * red-refinement: `UGGrid`, `FoamGrid`

---

# Adaptivity - Grid Refinement

.center[<img src="images/2d-red-refinement.png">]

---

# Adaptivity - Grid Refinement

- In DUNE: Grid and data is separated
- Grid adaptivity happens in two steps:
  1. Mark elements for refinement coarsening
  2. Initiate the grid adaption by `grid.adapt()` (actually 3 steps: preAdapt, adapt, postAdapt)
- In AMDiS, a grid wrapper establishes connection between Grid and data: `AdaptiveGrid`. This can be
  obtained by `ProblemStat::grid()`.

```c++
// Refines all grid elements count times.
void globalRefine (int count);

// Marks an entity to be refined/coarsened in a subsequent adapt
bool mark (int count, Element const& e);

// Return refinement mark for entity
int getMark (Element const& e) const;
```

---

# Adaptivity - Grid Refinement

```c++
// Prepare the grid for adaptation and notify observers of the preAdapt event
// returns true if an element might be coarsened
bool preAdapt();

// Adapt the grid and notify observers of the adapt event
// returns true if the grid was refined
bool adapt();

// Perform cleanup after grid adaptation and notify observers of the postAdapt event
void postAdapt();
```

---

# Adaptivity - Marking Elements

- Before the grid can be adapted, elements must be marked for refinement and coarsening.
- Only 1 level of refinement/coarsening at a time
- In the `mark` function, set `count=1` for refinement and `count=-1` for coarsening:

```c++
auto& grid = *prob.grid(); // ProblemStat returns shared_ptr
for (const auto& element : elements(gridView))
{
  if (indicator(element))
    grid.mark( 1, element); // Mark for refinement
  else
    grid.mark(-1, element); // Mark for coarsening
}
grid.preAdapt();
grid.adapt();
grid.postAdapt();
```

---

# Adaptivity - Marking Elements
- AMDiS provides classes to perform the marking of the grid
- Either based on error estimates, or marking indicators

```c++
template <class Grid, class GridFunction>
class GridFunctionMarker;
```

- Main interface methods: `markGrid(AdaptInfo& adaptInfo)`
- The `GridFunctionMarker` expects a grid-function that returns the target refinement level.
- Grid-function is evaluated in the center of the elements
- Element mark is computed from target level and current level.

---

# Adaptivity - Marking Elements
## Example

```c++
ProblemStat prob{"name", baseGrid, basis};
prob.initialize(INIT_ALL);

// marker with analytical function as level indicator
GridFunctionMarker marker{"marker", prob.grid(), [](auto const& x) -> int {
  return ((x - FieldVector<double,2>{0.5,0.5}).two_norm() < 0.3) ? 4 : 2;
}};

// mark elements for refinement/coarsening
AdaptInfo adaptInfo{"adapt"};
for (int i = 0; i < 10 && marker.markGrid(adaptInfo); ++i) {
  prob.grid()->preAdapt();
  bool refined = prob.grid()->adapt();
  prob.grid()->postAdapt();

if (!refined) break;
}
```

---

# Adaptivity - Marking Elements
- The `ProblemStat` manages markers attached to the grid
- Markers in the problem are invoked whenever the function `prob.markElements(adaptInfo)` is called.
- Register your own marker:

```c++
// register a new marker to the problem
template <class Marker>
void addMarker (Marker& m);

// Remove a marker with the given name from the problem
void removeMarker (std::string name);

// Remove a marker from the problem
void removeMarker (Marker<Grid> const& marker);
```

---

# Adaptivity - Marking Elements
## Example

```c++
ProblemStat prob{"name", baseGrid, basis};
prob.initialize(INIT_ALL);

prob.addMarker(marker);

AdaptInfo adaptInfo{"adapt"};
for (int i = 0; i < 10; ++i) {
  prob.markElements(adaptInfo);     // mark elements for refinement/coarsening
  if (!prob.adaptGrid(adaptInfo))   // perform the grid adaption
    break;
}
prob.removeMarker("marker");
```

---

# Adaptivity - Data Interpolation
- When the grid is adapted, the number of entities changes
- Global basis depends on the number of entities, thus needs to be updated
- DOFVector depends on the basis, thus needs to be updated

## Tree/chains of updates
- `grid.preAdapt()`
  * `dofVector.update(event::preAdapt)`
- `grid.adapt()`
  * `basis.update(event::adapt)`
      - `dofVector.update(event::adapt)`
      - `bilinearForm.update(event::adapt)`
      - `linearForm.update(event::adapt)`
- `grid.postAdapt()`
  * `dofVector.update(event::postAdapt)`

---

# Adaptivity - Data Interpolation
## Basis update
- Get a new GridView from the grid.
- Update the index mapping
- Update a parallel communicator (e.g. parallel index map)

## DOFVector update
- In `preAdapt` call `dataTransfer.preAdapt()`
- In `adapt` update size corresponding to updated basis, call `dataTransfer.adapt()`
- In `postAdapt` call `dataTransfer.postAdapt()`

## (Bi)LinearForm update
- Update all attached operators
- Reset sparsity pattern

---

# Adaptivity - Data Interpolation
## Data Transfer
- Each `DOFVector` has attached a `DataTransferInterface`.
- Responsible for "interpolating" data between grid changes
- Two implementations: `NoDataTransfer`, interpolation bases `DataTransfer`

### Interface
```c++
// Collect data that is needed before grid adaption
virtual void preAdapt (Container const& container, bool mightCoarsen);

// Interpolate data to new grid after grid adaption
virtual void adapt (Container& container);

// Perform cleanup after grid adaption
virtual void postAdapt (Container& container);
```

---

# Adaptivity - Data Interpolation
## Example

```c++
ProblemStat prob{"name", baseGrid, basisFactory};
prob.initialize(INIT_ALL);

DOFVector data1{prob.globalBasis()}; // default data transfer = interpolation
DOFVector data2{prob.globalBasis(), DataTransferOperation::INTERPOLATE};
DOFVector data3{prob.globalBasis(), DataTransferOperation::NO_OPERATION};

// change the data transfer operation
data1.setDataTransfer(DataTransferOperation::NO_OPERATION);

// write your own data transfer code, must be derived from DataTransferInterface
MyFancyDataTransfer myFancyDataTransfer{prob.basis()};
data2.setDataTransfer(myFancyDataTransfer);

// use the wrapped grid - AdaptiveGrid - to construct an adaptive DOFVector
DOFVector data4{prob.gridView(), lagrange<3>(), DataTransferOperation::INTERPOLATE};

// When using a default grid, then NO GRID ADAPTION!
DOFVector data5{baseGrid.leafGridView(), lagrange<3>()};
```

---

# Adaptivity - Observer Pattern
- How to react on grid changes? How to react on changes in the basis?
- Design patter: **Observer - Notifier**
- Observer:
  * waits for any calls to its `update(event)` function
  * must be registered to a notifier
- Notifier:
  * calls the update method on all attached observers, when `notify(event)` is called.

---

# Adaptivity - Observer Pattern
## Example

```c++
template <class Basis>
class MyClass
  : public  Notifier<event::adapt>  // can sent a signal that MyClass is updated
  , private Observer<event::adapt>  // wait for signals from the basis
{
public:
  MyClass (Basis const& basis)
    : Observer<event::adapt>(basis) // here, basis is the notifier we register to
  {}

protected:
  // Observer method must be implemented
  void updateImpl (event::adapt e) override {
    /* do some update */
    Notifier<event::adapt>::notify(e);
  }
};
```
---
layout: true
class: tud-dark, typo
---
class: center, middle, tud-dark
# Exercise 8
---

# Exercise 8

This exercise is about adaptivity and data transfer

1) Create an unstructured (adaptive) triangular grid of $\Omega=[-2,2]^2$, e.g., `UGGrid` or `ALUGrid`

2) Create a `ProblemStat` with a Lagrange P2 basis

3) Let $c(x)$ be a phase-field function, defined as

\[
  c(x) = 0.5(1 - \tanh(3\rho(x)/\epsilon)),\text{ with }\rho(x)=\|x - center\| - 0.3
\]

for $\epsilon=0.1$ and $center=(0,-1.5)^T$

4) Refine the grid along the interface with $0.05 <= c(x) <= 0.95$ to level 4 and level 0 outside the interface

---

# Exercise 8
## Advanced exercise

- Interpolate the function $\sin(4 x_1)$ into the solution DOFVector
- Iteratively move the center of the phase-field function upwards, by steps $dy=0.1$
- Update the grid refinement correspondingly
- Use a `GridFunctionMarker` that gets registered to the `ProblemStat`
- Visualize in each step the solution vector
- Do the same, but this time set as data transfer operation `NO_OPERATION`.

---
layout: true
class: tud-light, typo
---
class: center, middle
# Linear Algebra Backends and Solvers
---

# Linear Algebra Backends and Solvers

- AMDiS includes support for 4 linear algebra backends
  - dune-istl: parallel, use block structure, (P)AMG
  - PETSc: parallel, many solvers
  - MTL4: sequential, HYPRE AMG support
  - EIGEN: sequential

- Backend chosen when compiling AMDiS (default ISTL)

```bash
cmake -DBACKEND=[ISTL, PETSC, MTL, EIGEN] build-cmake
```

---
# Linear Algebra Backends and Solvers

- Backends define their own types for matrices and vectors
- Each backend has its own set of linear solvers (and preconditioners)
- Linear solver and preconditioner can be chose and configured in the parameter file:
```bash
prob->solver: <solver_name>
prob->solver->parameter: <value>
```
- Solvers are named like in the backend libraries

## Backend: ISTL
- We follow dune-istl naming of parameters in the `ISTL` backend

```bash
prob->solver: pcg
prob->solver->info: 1
prob->solver->maxit: 1000
prob->solver->reduction: 1.e-10
prob->solver->precon: ilu
```

---
# Linear Algebra Backends and Solvers

## Backend: PETSC
- We follow PETSc naming of parameters in the `PETSC` backend

```bash
prob->solver: cg
prob->solver->info: 1
prob->solver->max it: 1000
prob->solver->rtol: 1.e-10
prob->solver->pc: ilu
```
---
# Linear Algebra Backends and Solvers

## Backend: MTL
- We follow MTL4 naming of parameters in the `MTL` backend

```bash
prob->solver: cg
prob->solver->info: 1
prob->solver->print cycle: 10
prob->solver->max iteration: 1000
prob->solver->relative tolerance: 1.e-10
prob->solver->precon: ilu
```

---

# ISTL Linear Algebra Backend

- List of iterative solvers:
  - `cg`, `pcg`, `fcg`, `cfcg`, `bicgstab` (default), `minres`, `gmres`, `fgmres`
- List of direct solvers (squential):
  - `umfpack` (direct), `ldl`, `spqr`, `cholmod`, `superlu`
- List of preconditioners:
  - `jacobi` (diag), `gauss_seidel` (gs), `sor`, `ssor`, `richardson` (default), `solver`, `ilu(0)`, `ildl`, `amg`, `fastamg`, `kamg`
- List of parallel preconditioners:
  - `pssor`, `bjacobi`, (mabye in the future: `paamg`)

---

# ISTL Linear Algebra Backend

- You can configure solvers recursively, using the `solver` preconditioner:

```bash
prob->solver: pcg
prob->solver->precon: solver
prob->solver->precon->solver: pcg
prob->solver->precon->solver->precon: solver
prob->solver->precon->solver->precon->solver: pcg
...
```

---

# ISTL Linear Algebra Backend
### Directly use the dune-istl backend

```c++
#include <dune/istl/solvers.hh>

auto const& M = prob.systemMatrix()->impl().matrix();
auto      & X = prob.solutionVector()->impl().vector();
auto const& Y = prob.rhsVector()->impl().vector();

Dune::MatrixAdapter<TYPEOF(M),TYPEOF(X),TYPEOF(Y)> linOp{M};
Dune::Richardson<TYPEOF(X),TYPEOF(Y)> precon{0.5};
Dune::GradientSolver solver{linOp, precon, /*reduct*/1.e-6, /*maxit*/100, /*verbose*/1};

Dune::InverseOperatorResult statistics;
solver.apply(X, Y, statistics);
```

---

# ISTL Linear Algebra Backend
### Directly use the dune-istl backend

```c++
template <class Traits>
class MyProblemStat : public ProblemStat<Traits>
{
public:
  using ProblemStat<Traits>::ProblemStat;

void solve(AdaptInfo& adaptInfo, bool=true, bool=false) override
  {
    // put the solver code here
  }
};
```

---

# ISTL Linear Algebra Backend
### Algebraic MultiGrid

- Powerful and efficient algebraic multigrid preconditioner
- Konfiguration with Initfile parameters

```bash
ellipt->solver: pcg
ellipt->solver->info: 1
ellipt->solver->maxit: 10000
ellipt->solver->reduction: 1.e-10
ellipt->solver->precon: amg  % [fastamg,amg,kamg]
ellipt->solver->precon->smoother: ssor
ellipt->solver->precon->smoother->iteration: 2
ellipt->solver->precon->smoother->relaxationFactor: 0.75
ellipt->solver->precon->preSmoothSteps: 1
ellipt->solver->precon->postSmoothSteps 1
ellipt->solver->precon->gamma: 1
```

.footnote1[(3) M. Blatt†, O. Ippisch, P. Bastian, *A Massively Parallel Algebraic MultigridPreconditioner based on Aggregation for EllipticProblems with Heterogeneous Coefficients*, 2013]

---

# ISTL Linear Algebra Backend
### Algebraic MultiGrid

- Ellipt problem example:

```bash
Level 0 has 263169 unknowns, 263169 unknowns per proc (procs=1)
aggregating finished.
Level 1 has 65153 unknowns, 65153 unknowns per proc (procs=1)
aggregating finished.
Level 2 has 16195 unknowns, 16195 unknowns per proc (procs=1)
aggregating finished.
Level 3 has 3974 unknowns, 3974 unknowns per proc (procs=1)
aggregating finished.
Level 4 has 966 unknowns, 966 unknowns per proc (procs=1)
operator complexity: 1.17411
Using a direct coarse solver (UMFPack)
Building hierarchy of 5 levels (including coarse solver) took 1.69902 seconds.
```

---

# ISTL Linear Algebra Backend
### Algebraic MultiGrid

- AMG precondition:
```bash
=== Dune::GeneralizedPCGSolver
*=== rate=0.521386, T=11.6638, TIT=0.197691, IT=59
solution of discrete system needed 13.3763 seconds
```
- FastAMG preconditioner:
```bash
=== Dune::GeneralizedPCGSolver
 Preprocessing Dirichlet took 0.0134938
*=== rate=0.535372, T=5.90178, TIT=0.10003, IT=59
solution of discrete system needed 7.64907 seconds
```
- ILU preconditioner:
```bash
=== Dune::GeneralizedPCGSolver
*=== rate=0.820344, T=17.9574, TIT=0.153482, IT=117
solution of discrete system needed 23.569 seconds
```

---
# PETSC Linear Algebra Backend

- Install PETSc: e.g. `sudo apt install petsc-dev`
- Set Backend in AMDiS: `cmake -DBACKEND=PETSC build-cmake`
--

- PETSc allows *sequential* and *parallel* solvers/preconditioners
- We follow PETSc naming of parameters in the `PETSC` backend

```bash
prob->solver: cg
prob->solver->info: 1
prob->solver->max it: 1000
prob->solver->rtol: 1.e-10
prob->solver->pc: ilu
```

```bash
iteration    0: resid 2.020775e+00
iteration   10: resid 3.742518e-02
iteration   20: resid 6.077763e-03
iteration   30: resid 1.013754e-04
iteration   40: resid 1.188473e-06
iteration   50: resid 8.657256e-09
Linear solve converged due to CONVERGED_RTOL iterations 56
```

---
# PETSC Linear Algebra Backend

- Install PETSc: e.g. `sudo apt install petsc-dev`
- Set Backend in AMDiS: `cmake -DBACKEND=PETSC build-cmake`
- PETSc allows *sequential* and *parallel* solvers/preconditioners
- We follow PETSc naming of parameters in the `PETSC` backend

```bash
prob->solver: cg
prob->solver->info: 1
prob->solver->max it: 1000
prob->solver->rtol: 1.e-10
prob->solver->pc: ilu
```

- Note: meaning of `rtol` (*relative tolerance*) is different from ISTL

---
# PETSC Linear Algebra Backend
### Direct solvers

- PETSc does not have the notion of direct solvers, but lu preconditoner with matsolver packages
- Either use `preonly` solver or `richardson` solver

```bash
prob->solver: preonly
prob->solver->pc: lu
prob->solver->pc->mat solver: umfpack
```

---
# PETSC Linear Algebra Backend
### Direct solvers

- PETSc does not have the notion of direct solvers, but lu preconditoner with matsolver packages
- Either use `preonly` solver or `richardson` solver

```bash
prob->solver: richardson
prob->solver->max it: 1
prob->solver->pc: lu
prob->solver->pc->mat solver: umfpack
```

---
# PETSC Linear Algebra Backend
### Direct solvers

- PETSc does not have the notion of direct solvers, but lu preconditoner with matsolver packages
- Either use `preonly` solver or `richardson` solver

```bash
prob->solver: direct
prob->solver->mat solver: umfpack
```

---
# PETSC Linear Algebra Backend
### KSP preconditioner

- Use a solver as preconditioner
- Recursive parametrization of the solver/preconditioner

```bash
prob->solver: cg
prob->solver->pc: ksp
prob->solver->pc->ksp: cg
prob->solver->pc->ksp->pc: ksp
prob->solver->pc->ksp->pc->ksp: cg
...
```

---
# PETSC Linear Algebra Backend

- More complicated parametrization of solvers and preconditioners by command-line parameters
- Set `prefix` to name individual solvers:

```bash
prob->solver: cg
prob->solver->prefix: xyz
...
```

```bash
<executable> <initfile> -xyz_pc_type ilu
```

---
# PETSC Linear Algebra Backend

- Set `->info` parameter to 10 to enable `KSP_View` for the solver, i.e. list all PETSc parameters
  set for the solver and preconditioner...

```bash
KSP Object: 1 MPI processes
  type: cg
  maximum iterations=10000
  tolerances:  relative=1e-10, absolute=1e-50, divergence=10000.
  left preconditioning
  using nonzero initial guess
  using DEFAULT norm type for convergence test
PC Object: 1 MPI processes
  type: ilu
  PC has not been set up so information may be incomplete
    ILU: out-of-place factorization
    0 levels of fill
    tolerance for zero pivot 2.22045e-14
    matrix ordering: natural
  ...
```

---
# ProblemStat::solve
## The process of solving the linear system

```c++
template <class Traits>
ProblemStat<Traits>::solve(AdaptInfo& adaptInfo, bool createMatrixData = true,
                                                 bool storeMatrixData = false) {
  // prepare the solution vector
  solution_->resize();

// copy matrix into the solver (if necessary)
  if (createMatrixData)
    linearSolver_->init(systemMatrix_->impl());

// solve the linear system
  Dune::InverseOperatorResult stat;
  linearSolver_->apply(solution_->impl(), rhs_->impl(), stat);

// maybe destroy matrix data
  if (!storeMatrixData)
    linearSolver_->finish();

// some output of the residuum...
}
```

---

# PETSc LinearSolver
## Example: The basic implementation of the PETSc LinearSolver
```c++
/// Implements \ref LinearSolverInterface::init()
void init (Matrix const& A) override
{
  finish();
  KSPCreate(<comm>, &ksp_);
  KSPSetOperators(ksp_, A.matrix(), A.matrix());

// read parameters and configure the ksp object
  initKSP(ksp_, prefix_);
  initialized_ = true;
}

/// Implements \ref LinearSolverInterface::finish()
void finish () override
{
  if (initialized_)
    KSPDestroy(&ksp_);
  initialized_ = false;
}
```

---

# PETSc LinearSolver
## Example: The basic implementation of the PETSc LinearSolver
```c++
/// Implements \ref LinearSolverInterface::apply()
void apply (Vector& x, Vector const& b, Dune::InverseOperatorResult& stat) override
{
  KSPSolve(ksp_, b.vector(), x.vector());

if (info_ > 0)
    KSPReasonView(ksp_, PETSC_VIEWER_STDOUT_WORLD);

KSPConvergedReason reason;
  KSPGetConvergedReason(ksp_, &reason);

stat.converged = (reason > 0);
}
```

---

# PETSc LinearSolver - Implement your own
```c++
template <class M, class X, class Y>
class MyPetscSolver : public LinearSolverInterface<M,X,Y>
{
  KSP ksp_;

void init (Matrix const& A) override {
    KSPCreate(PETSC_COMM_SELF, &ksp_);
    KSPSetOperators(ksp_, A.matrix(), A.matrix());
    KSPSetType(ksp_, KSPBCGSL);
  }

void finish () override  {
    KSPDestroy(&ksp_);
  }

void apply (Vector& x, Vector const& b, Dune::InverseOperatorResult& stat) override {
    KSPSolve(ksp_, b.vector(), x.vector());
  }
};
```

---

# PETSc LinearSolver - Implement your own
```c++
ProblemStat prob{"name", grid, basis_factory};
using Prob = decltype(prob);

// get the types of matrix and vector backend from ProblemStat
using M = typename Prob::Mat;
using X = typename Prob::Vec;
using Y = X;

// register the solver in the problem
MyPetscSolver<M,X,Y> mySolver{};
prob.setSolver(mySolver);
```

---
layout: true
class: tud-dark, typo
---
class: center, middle, tud-dark
# Exercise 9
---

# Exercise 9

This exercise is about solving the linear system that came from a finite-element discretization of
a Poisson equation on the unit square:

\[
-\Delta u = f \text{ in }\Omega=[0,1]^2,\quad u=0\text{ on }\partial\Omega
\]

with $f=1$.

- Discretize this equation using Lagrange finite-elements of order 2
- Assemble and solve the equation using the parameter file
- Try different backends and different solvers and preconditioners

---

# Exercise 9
## Advanced exercise

- Extract the matrix and the vectors from the `ProblemStat`
- Solve directly using the solver/preconditioner routines from the backend
- Put this solver strategy into a Solver class and add to the Problem
- Create your own `ProblemStat` class that derives from the default `ProblemStat` and override the
  solve method to incorporate your own solver code there.

---
layout: true
class: tud-light, typo
---
class: center, middle
# Parallel Grids and Parallel Solvers
---

# Parallel Grids and Parallel Solvers
- AMDiS is by default aware of parallelization
- Parallel grid __and__ parallel linear algebra backend required

.center[
<img height="300" src="images/grid-partitions.png">
]
Decomposition of an adaptively refined grid into eight subdomains

.footnote1[(4) Image from: O. Sander, *DUNE - The Distributed and Unified Numerics Environment*, 2020]

---
# Parallel Grids
## Distributed Grids
- Grids with *ghost* elements

.center[
<img height="300" src="images/distributed-grid-ghosts.png">
]

.footnote1[(4) Image from: O. Sander, *DUNE - The Distributed and Unified Numerics Environment*, 2020]

---

# Parallel Grids and Parallel Solvers
## Distributed Grids
- Grids with *overlapping* elements

.center[
<img height="300" src="images/distributed-grid-overlap.png">
]

.footnote1[(4) Image from: O. Sander, *DUNE - The Distributed and Unified Numerics Environment*, 2020]

---

# Parallel Grids
## Creating a Parallel Grid
- Some grids can be constructed distributed directly
  - Structured grids (`YaspGrid` or `SPGrid`) with parameters for overlap
  - Unstructured grid `ALUGrid` from special distributed grid files, e.g., with [dune-vtk](https://gitlab.mn.tu-dresden.de/iwr/dune-vtk)
- Other grids require the initial construction on rank 0 and perform a distribution afterwards, using `grid.loadBalance()`.
  - Example: `UGGrid`
- Some grids are not parallel on its own, e.g. `AlbertaGrid`, `FoamGrid`
  - Require "parallelization wrapper" in form of a Meta-Grid, see [dune-metagrid](https://gitlab.dune-project.org/extensions/dune-metagrid).

### Example 1
```c++
std::bitset<dim> periodic("00");  // Not periodic in either of the two directions
int overlapSize = 1;              // Thickness of the overlap layer
Dune::YaspGrid<2> grid({1.0,1.0}, {10,10}, periodic, overlapSize, MPI_COMM_WORLD)
```

---

### Example 2
```c++
using Factory = Dune::StructuredGridFactory< Dune::UGGrid<2> >;
auto gridPtr = Factory::createSimplexGrid({0.0,0.0}, {1.0,1.0}, {10u,10u});
gridPtr->loadBalance();
```

---

# Parallel Grids
## Partition Type of Entities

- Each entity in a process has a **partition type** assigned to it.
- There are five different possible partition types:
  `interior` (grey), `border` (blue), `overlap` (green), `front` (magenta), and `ghost` (yellow).

.center[
<img height="200" src="images/partition-type1.png">
]

Partition with `interior`, `overlap`, and `ghost` elements

.footnote1[(4) Image from: O. Sander, *DUNE - The Distributed and Unified Numerics Environment*, 2020]

---

# Parallel Grids
## Partition Type of Entities

.center[
<img height="200" src="images/partition-type2.png">
]

Partition with `interior` and `overlap` elements

.footnote1[(4) Image from: O. Sander, *DUNE - The Distributed and Unified Numerics Environment*, 2020]

---

# Parallel Grids
## Partition Type of Entities

.center[
<img height="200" src="images/partition-type3.png">
]

Partition with `interior` and `ghost` elements

.footnote1[(4) Image from: O. Sander, *DUNE - The Distributed and Unified Numerics Environment*, 2020]

---

# Parallel Grids
## Partition Type of Entities
- You can traverse all entities of a specific group of partition types
- Groups are defined as (unions of) sets: e.g. `Partititons::interior`, `Partititons::border`, or
  the union `Partititons::interior + Partititons::border`, or all (local) entities: `Partititons::all`

```c++
std::size_t counter = 0;
for (const auto& vertex : vertices(gridView, Partitions::interior + Partitions::border))
{
  if (vertex.partitionType() == InteriorEntity)
    counter++;
}
```

> **Note:** `gridView.indexSet()` enumerates the entities of the `all` partition!

---

# Parallel Grids
## Grid Partititoning
- Grids provide default partititon, when calling `loadBalance()`
- Sometimes, you want to use a more powerful partitioning algorithms
- Some grids provide extended interface (not unified)

### Example: UGGrid
```c++
bool loadBalance(const std::vector<Rank>& targetProcesses,
                 unsigned int fromLevel);
```

Create initial partitioning using *ParMETIS*
```c++
#include <dune/grid/utility/parmetisgridpartitioner.hh>
...
std::vector<unsigned int> part
  = ParMetisGridPartitioner<GridView>::partition(gridView, mpiHelper);
gridPtr->loadBalance(part, 0);
```

---

### Example: UGGrid
```c++
bool loadBalance(const std::vector<Rank>& targetProcesses,
                 unsigned int fromLevel);
```

Refine partititoning by redistribution:
```c++
#include <dune/grid/utility/parmetisgridpartitioner.hh>
...
std::vector<unsigned int> part
  = ParMetisGridPartitioner<GridView>::repartition(gridView, mpiHelper, 1000f);
gridPtr->loadBalance(part, 0);
```

---

# Parallel Grids and Parallel Solvers

- Data structures and linear algebra backend must be aware of distributed grids
- Two paradigms:
  1. **Domain decomposition** (local data structures, solvers handle parallelization)
  2. **Distributed data structures** (data structures handle parallelization, "standard" solvers)

In the first category falls dune-istl, in the second category PETSc (and in the future PMTL)

- In AMDiS, domain decomposition requires grid **overlap**:<br>
  `grid.overlapSize() > 0`
- In AMDiS, distributed data structured allows grids with **ghost** cells:<br>
  `grid.ghostSize() >= 0`

---

# Parallel Grids and Parallel Solvers
## ISTL Backend
- Create a parallel grid with overlap, e.g. `Dune::YaspGrid`, `Dune::SPGrid`.
- Assembling happens on the `Partition::all` (local) entities.
- Only iterative solvers can be used!
- Local data structures are always up-to-date, since overlapping ensures synchronization

## Examples
- Choose parallel preconditioner: ParSSOR (`pssor`), BlockJacobi (`bjacobi`), or AMG

```tex
prob->solver: pcg
prob->solver->precon: bjacobi
prob->solver->precon->sub precon: ilu   % precon to be applied on subdomains
```

---

# Parallel Grids and Parallel Solvers
## PETSc Backend

- Create a parallel grid w/ or w/o overlap/ghost, e.g. `Dune::UGGrid`, `Dune::ALUGrid`,...
- Basis provides `DistributedCommunication` that holds a parallel DOF-map, mapping local to global indices.
- Assembling happens on the `Partition::interior` entities.
- Ghost/Overlapping entities are used to collect data from neighbouring processors during synchronization.
- Classification of each DOF as owner and non-owner DOF on each processor, i.e., a DOF belongs to exactly one
  process.
- Vectors and Matrices perform communication of data on insertion and before access.
- Data is automatically synchronized. But: Do not mix read and write access to data $\to$ expensive

---

# Parallel Grids and Parallel Solvers
## PETSc Backend
### Examples

- PETSc Matrix type: *MATMPIAIJ* (parallel sparse matrix), PETSc Vector type: *VECMPI* (parallel vector with ghost padding)
- Support all solvers and preconditioners for these matrix/vector types

```tex
prob->solver: cg
prob->solver->pc: bjacobi
prob->solver->pc->sub ksp: preonly    % solver to be applied on subdomains
prob->solver->pc->sub ksp->pc: ilu    % preconditioner to use for ths sub solver
```

---
layout: true
class: tud-light, typo
---
class: center, middle
# Code organization and Object-Oriented Interfaces
---

# Code organization and Object-Oriented Interfaces

## Classes to organize solution process

`ProblemStat`
: Holds all data to assemble and solve a stationary problem

`ProblemInstat`
: Contains some procedures to handle instationary problems

`AdaptStationary`
: Describes the (adaptive) solution process of a stationary problem

`AdaptInstationary`
: Describes the (adaptive) solution of an evolutionary problem

> What is the relation between these classes?

> How to control the solution process?

> How to interact with the classes and change its behaviour?

---

# Inheritance Diagram

.pure-g[.pure-u-1-2[
## Stationary Problem
<div class="mermaid">
    classDiagram
    ProblemIterationInterface <|-- StandardProblemIteration
    ProblemIterationInterface: +beginIteration(AdaptInfo)
    ProblemIterationInterface: +oneIteration(AdaptInfo, Flag)
    ProblemIterationInterface: +endIteration(AdaptInfo)

StandardProblemIteration <|-- ProblemStat
    StandardProblemIteration: +ProblemStatBase*

ProblemStatBase <|-- ProblemStat
    ProblemStatBase: +markElements(AdaptInfo)
    ProblemStatBase: +adaptGrid(int n)
    ProblemStatBase: +globalCoarsen(int n)
    ProblemStatBase: +globalRefine(int n)
    ProblemStatBase: +buildAfterAdapt(AdaptInfo, Flag, bool, bool)
    ProblemStatBase: +solve(AdaptInfo, bool, bool)
    ProblemStatBase: +estimate(AdaptInfo)

ProblemStat: +initialize(Flag, ProblemStat, Flag)
    ProblemStat: +restore(Flag)
    ProblemStat: +solution(...)
    ProblemStat: +solutionVector(...)
</div>
].pure-u-1-2[
## Instationary Problem
<div class="mermaid">
    classDiagram
    ProblemTimeInterface <|-- ProblemInstatBase
    ProblemTimeInterface: +initTimeInterface()
    ProblemTimeInterface: +setTime(AdaptInfo)
    ProblemTimeInterface: +initTimestep(AdaptInfo)
    ProblemTimeInterface: +closeTimestep(AdaptInfo)
    ProblemTimeInterface: +solveInitialProblem(AdaptInfo)
    ProblemTimeInterface: +transferInitialSolution(AdaptInfo)

ProblemInstatBase <|-- ProblemInstat
    ProblemInstatBase: +time()
    ProblemInstatBase: +tau()
    ProblemInstatBase: +invTau()

ProblemInstat: +initialize(Flag)
    ProblemInstat: +ProblemStat*
    ProblemInstat: +oldSolution(...)
    ProblemInstat: +oldSolutionVector(...)
</div>
]]

---

# Stationary Solution Process
## `AdaptStationary::adapt()`

1. Initial iteration
```c++
problemIteration->beginIteration(adaptInfo);
problemIteration->oneIteration(adaptInfo, NO_ADAPTION);
problemIteration->endIteration(adaptInfo);
```

2. Adaption loop
```c++
while (!adaptInfo.spaceToleranceReached()) {
      problemIteration->beginIteration(adaptInfo);
      problemIteration->oneIteration(adaptInfo, FULL_ITERATION);
      problemIteration->endIteration(adaptInfo);
}
```

### Iteration Flags

`NO_ADAPTION`
: BUILD | SOLVE | ESTIMATE

`FULL_ITERATION`
: MARK | ADAPT | BUILD | SOLVE | ESTIMATE

---

# Stationary Solution Process
## `StandardProblemIteration::oneIteration()`

<div class="mermaid">
graph LR
  id0["mark elements"]
  id1["adapt grid"]
  id2["assemble"]
  id3["solve"]
  id4["estimate"]
  id0 --> id1 --> id2 --> id3 --> id4
</div>

1. Adapt the grid and assemble the system
```c++
buildAndAdapt(adaptInfo, todo) {
      Flag markFlag = problem.markElements(adaptInfo);
      if (markFlag) problem.adaptGrid(adaptInfo);
      problem.buildAfterAdapt(adaptInfo, markFlag, true, true);   // build matrix
      problem.buildAfterAdapt(adaptInfo, markFlag, false, true);  // build rhs
}
```

2. Solve the system and estimate the error
```c++
problem.solve(adaptInfo, true, false);  // createMatrixData and !storeMatrixData
problem.estimate(adaptInfo);
```

---

# Stationary Solution Process
## Example: Nonlinear Iteration

Implementation of a Newton iteration inside of `oneIteration()`:

```c++
template <class Problem>
struct NewtonIteration : public ProblemIterationInterface
{
  Problem& prob_;
  using SolutionVector = typename Problem::SolutionVector;
  std::shared_ptr<SolutionVector> stepSolution_;

NewtonIteration (Problem& prob) : prob_(prob)
  { // create a vector fo the step solution
    stepSolution_.reset(new SolutionVector(*prob_->solutionVector()));
  }

void beginIteration(AdaptInfo& adaptInfo) override { ... }
  Flag oneIteration(AdaptInfo& adaptInfo, Flat toDo) override;
  void endIteration(AdaptInfo& adaptInfo) override { ... }
  // implementation of some more abstract functions...
}
```

---

# Stationary Solution Process
## Example: Nonlinear Iteration

Implementation of a Newton iteration inside of `oneIteration()`:

```c++
template <class Problem>
Flag NewtonIteration<Problem>::oneIteration(AdaptInfo& adaptInfo, Flat toDo) override
{
  if (adaptInfo.spaceIteration() == 1) {
    prob_.buildAfterCoarsen(adaptInfo, 0, true, true);  // assemble DF(u0)
    prob_.solver()->init(prob_.systemMatrix()->impl());
  } else
    prob_.buildAfterCoarsen(adaptInfo, 0, false, true); // assemble F(u)

stepSolution_->resizeZero();                          // Initial guess is zero

Dune::InverseOperatorResult stat;                     // solve the linear system
  prob.solver()->apply(stepSolution_->impl(), prob_.rhsVector()->impl(), stat);
                                                        // u = u + d
  prob_.solutionVector()->impl().vector() += stepSolution_->impl().vector();
  adaptInfo.setEstSum(/*error estimation*/, "");
  return Flag(0);
}
```

---

# Stationary Solution Process
## Example: Nonlinear Iteration
How to use the `NewtonIteration`?

```c++
ProblemStat prob{"name", grid, basisFactory};
prob.initialize(INIT_ALL);

NewtonIteration newton{prob}; // CTAD

auto u = prob.solution();
prob.addMatrixOperator(...);  // fill jacobian DF(u)
prob.addVectorOperator(...);  // fill objective function -F(u)

AdaptInfo adaptInfo{"adapt"};
AdaptStationary adaptStat{"adapt", newton, adaptInfo};
adaptStat.adapt();
```

---

# Stationary Solution Process
## Example: Nonlinear Iteration
In parameter file we can control the Newton iterations:
```tex
adapt->max iteration: 10    % maximal newton iterations
adapt[]->tolerance: 1.e-6   % Newton break tolerance for whole system
```

---

# Intationary Solution Process
## `AdaptInstationary::adapt()`

1. Initial timestep
```c++
if (adaptInfo.timestepNumber() == 0) {
      adaptInfo.setTime(adaptInfo.startTime());
      problemTime->setTime(adaptInfo);
      problemTime->solveInitialProblem(adaptInfo); // maybe initialAdaptInfo
      problemTime->transferInitialSolution(adaptInfo);
}
```

2. Adaption loop
```c++
while (!adaptInfo.reachedEndTime()) {
      problemTime->initTimestep(adaptInfo);
      oneTimestep(); // different strategies
      problemTime->closeTimestep(adaptInfo);
}
```

With `oneTimestep()` either `explicitTimeStrategy()` (0), `implicitTimeStrategy()` (1), `simpleAdaptiveTimeStrategy()` (2)

---

# Intationary Solution Process
## `AdaptInstationary::explicitTimeStrategy()`

```c++
adaptInfo.setTimestepIteration(0);

// estimate before first adaption
if (adaptInfo.time() <= adaptInfo.startTime())
  problemIteration->oneIteration(adaptInfo, ESTIMATE);

// increment time
adaptInfo.setTime(adaptInfo.time() + adaptInfo.timestep());
problemTime->setTime(adaptInfo);

adaptInfo.setSpaceIteration(0);

// do the iteration
problemIteration->beginIteration(adaptInfo);
problemIteration->oneIteration(adaptInfo, FULL_ITERATION);
problemIteration->endIteration(adaptInfo);

adaptInfo.setLastProcessedTimestep(adaptInfo.timestep());
adaptInfo.incTimestepNumber();
```

---

# Instationary Solution Process
## Example: Monitor solution in each timestep

```c++
template <class Problem>
struct MonitorProblemInstat : public ProblemTimeIteration
{
  using SolutionVector = typename Problem::SolutionVector;
  std::shared_ptr<SolutionVector> oldSolution_;

MonitorProblemInstat (Problem& prob) : prob_(prob) { // create old solution vector
    oldSolution_.reset(new SolutionVector(*prob_->solutionVector()));
  }

void initTimestep(AdaptInfo& adaptInfo) override {
    *oldSolution_ = *prob_.solutionVector();  // copy the solution to oldSolution
  }
  void closeTimestep(AdaptInfo& adaptInfo) override {
    prob_.writeFiles(adaptInfo);              // write the solution to file
  }
  void transferInitialSolution(AdaptInfo& adaptInfo) override { ... }
  // some more interface methods
};
```

---

# Instationary Solution Process
## Example: Monitor solution in each timestep

Use that instationary Problem to run your solution process:
```c++
ProblemStat prob{"name", grid, basisFactory};
prob.initialize(INIT_ALL);

MonitorProblemInstat probInstat{prob}; // CTAD

auto u = prob.solution();
auto uOld = valueOf(*probInstat.oldSolution_); // TODO: add nice interface method
prob.addMatrixOperator(...);
prob.addVectorOperator(...);

AdaptInfo adaptInfo{"adapt"};
AdaptInstationary adaptInstat{"adapt", prob, adaptInfo, probInstat, adaptInfo};
adaptInstat.adapt();
```

---

# Instationary Solution Process
## The `ProblemInstat`
The default implementation of a `ProblemTimeInterface` is called `ProblemInstat` and provides lots of
utility functions to time-dependent processes:

```c++
ProblemStat probInstat{"name", prob};

auto invTau = std::ref(probInstat.invTau()); // a reference to 1/timestep
prob.addMatrixOperator(zot(invTau));
prob.addVectorOperator(zot(invTau * probInstat.oldSolution())); // uOld

prob.addVectorOperator(zot(invokeAtQP([&t=probInstat.time()] { // reference to time
  return std::sin(t);
})));
```

> **NOTE**: the `std::ref` cannot be multiplied with `double`. This would lead to `double` and thus the
> reference character is lost. Use a lambda for this:<br>`invokeAtQP([&a=a,&b=b]{ return a*b; })`

---

# AMDiS Extensions

- The `ProblemInstat` always needs a `ProblemStat`. Why not combine these two?
- It would be nice to have a method where we can put all the operators

### BaseProblem
- Wrapper class that implements `ProblemIterationInterface` and `ProblemTimeInterface` and also `ProblemInstatBase`
- You need to implement the methods `fillOperators()` and `fillBoundaryConditions()`

---

# AMDiS Extensions : `BaseProblem`

```c++
struct BaseProblemInterface
    : public virtual ProblemIterationInterface
    , public virtual ProblemTimeInterface
{
  virtual void initialize(Flag initFlag, BaseProblemInterface*, Flag adoptFlag) = 0;
  virtual void initData(AdaptInfo& adaptInfo) = 0;
  virtual void finalizeData(AdaptInfo& adaptInfo) = 0;

/// A method where operators are added to the problem
  virtual void fillOperators(AdaptInfo& adaptInfo) = 0;
  /// A method where boundary conditions are added to the problem
  virtual void fillBoundaryConditions(AdaptInfo& adaptInfo) = 0;
  /// Build up the operators
  virtual void initBaseProblem(AdaptInfo& adaptInfo) {
    initData(adaptInfo);
    fillOperators(adaptInfo);
    fillBoundaryConditions(adaptInfo);
    finalizeData(adaptInfo);
  }
  /// Calls writeFiles of the problem
  virtual void writeFiles(AdaptInfo& adaptInfo) = 0;
};
```

---

# AMDiS Extensions : `BaseProblem`
## Example

```c++
template <class Traits>
struct PoissonProblem : BaseProblem<Traits>
{
  template <class Grid, class BasisFactory>
  PoissonProblem(std::string name, Grid& grid, const BasisFactory& bf)
    : BaseProblem<Traits>(name, grid, bf) {}

void fillOperators(AdaptInfo& adaptInfo) override {
    this->problem().addMatrixOperator(sot(1.0));
    this->problem().addVectorOperator(zot(1.0));
  }
  void fillBoundaryConditions(AdaptInfo& adaptInfo) override {
    this->problem().addDirichletBC([]{return true;}, 0.0);
  }
};

// deduction guide
template <class Grid, class BasisFactory>
PoissonProblem(std::string name, Grid& grid, const BasisFactory& bf)
  -> PoissonProblem<Impl::DeducedProblemTraits_t<Grid,BasisFactory>>;
```

---

# AMDiS Extensions : `BaseProblem`
## Example

```c++
PoissonProblem prob{"ellipt", grid, lagrange<2>()};
prob.initialize(INIT_ALL);

AdaptInfo adaptInfo{"adapt"};
prob.initBaseProblem(adaptInfo);

// put the BaseProblem as probInstat and probTime
AdaptInstationary adaptInstat{"adapt", prob, adaptInfo, prob, adaptInfo};
adaptInstat.adapt();
```

---

# AMDiS Extensions
- What if we want to couple multiple problems

### CouplingIterationInterface
- Combine multiple `ProblemIterationInterface`s

### CouplingTimeInterface
- Combine multiple `ProblemtimeInterface`s

### CouplingBaseProblem and VariadicCouplingBaseProblem
- Combine multiple `BaseProblem`s

---

# AMDiS Extensions : `CouplingIterationInterface`

.center[<img src="images/coupling-iteration-interface.png" width="70%">]

---

# AMDiS Extensions : `CouplingIterationInterface`
## Example

```c++
AdaptiveGrid grid{baseGrid};  // create a shared grid
ProblemStat prob1{"prob1", grid, basisFactory1};
ProblemStat prob2{"prob2", grid, basisFactory2};

// register the two problems in the coupling iteration interface
CouplingIterationInterface probIter{};
probIter.addIterationInterface(prob1);
probIter.addIterationInterface(prob2);

AdaptInfo adaptInfo{"adapt"};
AdaptStationary adaptStat{"adapt", probIter, adaptInfo};
adaptStat.adapt();
```

---

# AMDiS Extensions : `CouplingTimeInterface`

.center[<img src="images/coupling-time-interface.png" width="70%">]

---

# AMDiS Extensions : `CouplingTimeInterface`
## Example

```c++
// create coupling iteration interface probInstat as before

ProblemInstat probInstat1{"prob1", prob1};
ProblemInstat probInstat2{"prob2", prob2};

// register the two instationary probs in the coupling time interface
CouplingTimeInterface probTime{};
probTime.addTimeInterface(probInstat1);
probTime.addTimeInterface(probInstat2);

AdaptInfo adaptInfo{"adapt"};
AdaptInstationary adaptInstat{"adapt", probIter, adaptInfo, probTime, adaptInfo};
adaptInstat.adapt();
```

---

# AMDiS Extensions : `CouplingBaseProblem`
## Example
```c++
AdaptiveGrid grid{baseGrid}; // create a shared grid
CahnhilliardProblem chProb{"ch", grid};
NavierStokesProblem nsProb{"ns", grid};

// couple two base problems
CouplingBaseProblem chnsProb{"chns", chProb, nsProb};
chnsProb.initialize(INIT_ALL);

AdaptInfo adaptInfo{"adapt"};
chnsProb.initBaseProblem(adaptInfo);

// put the BaseProblem as probInstat and probTime
AdaptInstationary adaptInstat{"adapt", chnsProb, adaptInfo, chnsProb, adaptInfo};
adaptInstat.adapt();
```

---

# ProblemStat `Traits`

The `Traits` type passed to the `ProblemStat` hold all the type necessary to build the problem
- The **grid**, wrapped into an `AdaptiveGrid`
- The **gridView** to traverse over
- A factory for the global basis
- The datatype for the matrix and vector coefficients
- Traits for the linear-algebra backend

---

# ProblemStat `Traits`
## Example: `NavierStokesTraits`
```c++
template <class BaseGrid>
struct NavierStokesTraits {
  using Grid = AdaptiveGrid_t<BaseGrid>;
  using GridView = typename Grid::LeafGridView;

template <class... Args>
  static auto create(const Args&... args) {
    return AMDiS::GlobalBasis{args..., composite(power<2>(lagrange<2>()), lagrange<1>())};
  }

using GlobalBasis = decltype(create(std::declval<GridView>()));
  using CoefficientType = double;
  using LinAlgTraits = BackendTraits<GlobalBasis>;
};
```

---

# ProblemStat `Traits`
## Example: `NavierStokesTraits` and `CahnHilliardTraits`

or simply
```c++
struct NavierStokesBasisFactory
{
  static auto create() { return composite(power<2>(lagrange<2>()), lagrange<1>()); }
};

struct CahnHilliardBasisFactory
{
  static auto create() { return power<2>(lagrange<2>()); }
};

template <class Grid>
using NavierStokesTraits = DefaultBasisCreator<Grid, NavierStokesBasisFactory>;

template <class Grid>
using CahnHilliardTraits = DefaultBasisCreator<Grid, CahnHilliardBasisFactory>;
```

---

# ProblemStat `Traits`
## Example: `BaseProblem`

```c++
template <class Grid>
struct NavierStokesProblem
  : BaseProblem<NavierStokesTraits<Grid>>
{
  NavierStokesProblem(Grid& grid)
    : BaseProblem<NavierStokesTraits<Grid>>("ns", grid)
  {}
  // ...
};

template <class Grid>
struct CahnHilliardProblem
  : BaseProblem<CahnHilliardTraits<Grid>>
{
  CahnHilliardProblem(Grid& grid)
    : BaseProblem<CahnHilliardTraits<Grid>>("ch", grid)
  {}
  // ...
};
```