AMDiS meets DUNE - Evolution of a Finite-Element Framework

<div class="frontpage">
  <small>Felix Müller and Simon Praetorius</small>
  <h1>AMDiS meets DUNE<br><small>Evolution of a Finite-Element Framework</small></h1>

<small>Dresden, 29. Oktober 2020</small>
</div>

---
layout: true
class: tud-light, typo
---

# Initial Example - Poisson Equation

```c++
#include <amdis/AMDiS.hpp>                              // Include essential headers
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

prob.addMatrixOperator(sot(1.0));                     // a(v,u) = <grad(v), grad(u)>
  prob.addVectorOperator(zot([](auto x) { ... }, 6) );  // b(v) = <v, f(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
}
```

---

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

- **Object-Oriented** C++-framework, partially generic
- Based on C-library ALBERT(A) 1.x
- **Simplicial** grids in **1d**, **2d**, and **3d**, including surface grids
- Adaptive refinement by **bisection** of fixed refinement-edge, including automatic **interpolation** of data
- Distributed grids with **load-balancing**, parallel solver library PETSc, non-overlapping domain decomposition
- **Lagrange** local-functions of degree 1-4 in all dimensions
- High-level simulation classes for problem description and time loops
- Rudimentary abstraction of time-stepping schemes and non-linear solvers

.center.font-xl[.icon-rechts[] *legacy* AMDiS]

---

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

- **modular toolbox** for solving partial differential equations (PDEs) with grid-based methods
- a generic **grid interface**, allowing to interface a range of very different grid implementations
- High-level interfaces for trial and test functions and generic discretization modules

.center.font-xl[.icon-rechts[] *(new)* AMDiS as DUNE discretization module]

---

# 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
```bash
cd BASE_DIR
git clone https://gitlab.dune-project.org/core/dune-common.git
git clone https://gitlab.dune-project.org/core/dune-grid.git
# ...
```
--
2. Download AMDiS
```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 `BASE_DIR` 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) { /* ... */ }
```

---

# My first AMDiS project

Configure, compile and link it against all dependencies, by using `dunecontrol`:
```bash
./dune-common/bin/dunecontrol --only=my-first-project all
```

--
## Recommendations
- Always create a **git repository** for each project, e.g., https://gitlab.mn.tu-dresden.de
- Frequently **push** changes to the repository, by [git add, git commit, git push](https://training.github.com/downloads/github-git-cheat-sheet.pdf)
- After a project is finished (or paper published), create a [git tag](https://git-scm.com/book/en/v2/Git-Basics-Tagging)
  of the current state

```bash
git init                                git tag -a v1.4 -m "description of version 1.4"
git remote add origin https://URL       git push origin v1.4
git add .
git commit -m "Initial commit"
git push -u origin master
```

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

# `#include "live example"`

---
layout: true
class: tud-light, typo
---

# Grids

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

---

# Grids

- describe the domain \( \Omega \), a \( d \)-manifold embedded in \( \mathbb{R}^m \)
- provide a triangulation \(  \mathcal{T}_h \) of the domain
- AMDiS uses implementations provided by Dune
- different strengths: unstructured, parallel, conforming, ...
- common interface, change implementation by changing a few lines of code
- topology and geometry separated in dune-grid and dune-geometry modules

---

# Grid Components

<img align="right" width="40%" src="images/reftree.png">
- provide access to grid elements using a `GridView`
  - topological information
  - `LevelGridView`: contains elements at a certain depth in the refinement tree
  - `LeafGridView`: contains the leaf elements in the refinement tree
  <br> .icon-rechts[] usually the view we are interested it
- can also access vertices, edges, faces from an element
- called entities in Dune

---

# Available Grid Implementations

see https://dune-project.org/doc/grids/ and ["The Dune framework" chapter 3.1.6](https://www.sciencedirect.com/science/article/pii/S089812212030256X#d1e5072)
- `YASPGrid`: parallel, structured cube grid of any dimension
- `AlbertaGrid`: simplicial grids in two and three dimensions with bisection refinement based on [ALBERTA](https://www.alberta-fem.de)
- `ALUGrid` (dune-alugrid): parallel, unstructured, simplices or cubes, conforming refinement based on bisection supported for simplices only
<br> .icon-rechts[] closest to legacy AMDiS mesh in functionality
- `UGGrid` (dune-uggrid): parallel, unstructured grid with mixed element types
- `FoamGrid` (dune-foamgrid): one and two-dimensional grids embedded in three-dimensional space including non-manifold grids with branches

---
# Meta-Grids
- wrapper around grid implementation
- change certain aspects and add functionality

- `GeometryGrid`: maps element geometries using a user-defined functor, can be used to implement moving domains
- `CurvedGrid` (dune-curvedgrid): maps flat geometries to curved geometries

.center[
<img src="images/sphere_p1_ref0.png" width="20%" style="margin-right:20px">
<img src="images/sphere_p4_ref0.png" width="20%" style="margin-left:20px">
]

---

# Finite Element Spaces

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

---

# Finite Element Spaces

- weak formulation \( a(v,u) = b(v) \) is given
  - bilinear form \( a \)
  - linear form \( b \)
  - \( v \in V_0 \) test space .icon-rechts[] row basis
  - \( u \in V_1 \) trial space .icon-rechts[] col(umn) basis

- example Poisson
\[ \Delta u = f \text{ on } \Omega, \qquad u = g \text{ on } \partial \Omega \]

\[ \Rightarrow \text{find } u \in V_1 = H^1_g \text{ s.t. } \underbrace{\int_\Omega \nabla v \cdot \nabla u}_{=: a(v,u)}
= \underbrace{\int_\Omega f v}_{:= b(v)} \qquad \forall v \in V_0 = H^1_0 \]

---

# Structure of an AMDiS Basis

- implementation done by dune-functions and dune-localfunctions

--
- dune-functions

- FE space defined by a set of global basis functions
  - structure of the basis can be modeled, e.g. vector-valued basis can be viewed as a tensor product of scalars
  - indexed by multiindices

--
- dune-localfunctions

- FE space locally defined by restriction of global basis to grid elements
  - grid elements can be characterized by their reference element
  - indexed by local indices

---

# The Basis Tree

- AMDiS does most of the work of using the indices: interpolation, writing files, assembly of global stiffness matrix can be defined on a node in the basis tree
  <br> .icon-rechts[] user needs to define the basis tree

--
- 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 = \underbrace{(V_x \otimes V_y \otimes V_z)}_{=: V_v} \otimes V_p \)
  - composite (different components) \( V = 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

--
- AMDiS requires integer constants as indices for composite nodes!
  <br>.icon-rechts[] shorthand `Dune::Indices::_i`

---

# 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>()`

```c++
using namespace Dune::Functions::BasisFactory;
auto factory = composite(
                 power<2>(
                   lagrange<2>()),
                 lagrange<1>());
ProblemStat prob{"example", grid, factory};

using namespace Dune::Indices;                // use _0 and _1
prob.solution(_0, 1).interpolate(Expression); // interpolate onto y-component of velocity
prob.solution(_1).interpolate(Expression);    // interpolate onto pressure
```

---

# The `ProblemStat` class

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

---

# The `ProblemStat` class

- main interface class
--

- setup functions:
  - constructor call: name, grid, basis
  - `initialize(Flag)`: creates the components, must always be called
  - `addMatrixOperator`, `addVectorOperator`: add operators to define the PDE
  - `addDirichletBC`, `addPeriodicBC`: add constraints
  - `addMarker`: add a marker for automatic adaptation
--

- adaptation and solving the problem
  - `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

---

# The `ProblemStat` class

- default implementations handle most cases
  - interface allows custom operators, markers
  - flexible basis structure
  - many solvers and filewriters supported out of the box
--

- user can inherit from `ProblemStat` and replace the implementations freely, e.g. override `assemble`
- 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 members can be accessed to allow maximum flexibility

---

# Operators

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

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

---

# Operators

- define a bilinear or linear form
- contain the information how to assemble an element matrix or vector
--

- AMDiS contains a lot of predefined operators for most situations
- `zot(Expr, QuadOrder)`: creates a zero-order term \( \langle v, c\, u\rangle \)
- `fot(Expr, Tag, QuadOrder)`: creates a first-order term
  - `tag::partial_trial{i}`:  \( \langle v, c\,\partial_i u\rangle \)
  - `tag::partial_test{i}`:  \( \langle\partial_i v, c\, u\rangle \)
  - `tag::grad_test`: \( \langle\nabla v, \mathbf{b}\, u\rangle \)
  - `tag::grad_trial`: \( \langle v, \mathbf{b}\cdot\nabla u\rangle \)
- `sot(Expr, QuadOrder)`: creates a second-order term \( \langle\nabla v, c\,\nabla u\rangle \), or \( \langle\nabla v, A\,\nabla u\rangle \)
- `sot_ij(Expr, i, j, QuadOrder)`: creates a second-order term of partial derivatives \( \langle\partial_i v, c\,\partial_j u\rangle \)

---

# Using Operators

- more operators available using `makeOperator`
- list of tags in the [online documentation](https://amdis.readthedocs.io/en/latest/reference/Operators/#tags)

```c++
template <class Tag, class Expr, class QuadOrder>
auto makeOperator(Tag tag, Expr&& expr, QuadOrder&& order);
```
--

- add operators to a problem

```c++
template <class Operator, class RowPath, class ColPath>
void addMatrixOperator(Operator const& op, RowPath row, ColPath col);

template <class Operator, class Path>
void addVectorOperator(Operator const& op, Path path);
```
- `Path` is the tree index of the node the operator relates to

---

# Constraints

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

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

---

# Constraints

- AMDiS implements dirichlet and periodic constraints (essential boundary conditions) by changing the stiffness matrix and right hand side
- can be applied to a specific basis node using row- and column indices

```c++
template <class Predicate, class RowPath, class ColPath, class Values>
void addDirichletBC(Predicate const& p, RowPath row, ColPath col, Values const& values);

template <class RowPath, class ColPath, class Values>
void addDirichletBC(BoundaryType id, RowPath row, ColPath col, Values const& values);

void addPeriodicBC(BoundaryType id, WorldMatrix const& A, WorldVector const& b);
```

- `Predicate` is a function from world coordinates to `bool`: constraint applies on all boundary parts where it evaluates to true
- `BoundaryType` uses a boundary manager to set ids to certain parts

---

# Other Constraints

- natural boundary conditions can be applied by providing an operator
- add operator on a part of the boundary

```c++
template <class Operator, class RowPath, class ColPath>
void addMatrixOperator(BoundaryType b, Operator const& op, RowPath row, ColPath col)

template <class Operator, class Path>
void addVectorOperator(BoundaryType b, Operator const& op, TreePath path);
```

--
- examples usage:
```c++
auto opNeumann = makeOperator(tag::test{}, 1.0);
prob.addVectorOperator(BoundaryType{2}, opNeumann);
```

---

# Adaptation

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

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

---

# Adaptation

- controlled by an `AdaptInfo` object
- loads parameters from a parameter file
<br>.icon-rechts[] program behaviour can be changed at runtime, no recompilation necessary
--

- marking step
  - AMDiS has built-in markers
  - `prob.markElements(adaptInfo);`
  - manual marking with `grid.mark`
- `prob.adaptGrid(adaptInfo);`
- `prob.assemble(adaptInfo);`
- `prob.solve(adaptInfo);`
- `prob.writeFiles(adaptInfo);`

---

# The `AdaptStationary` Class
- default adaptation procedure implemented by `AdaptStationary` class
- must `#include <amdis/AdaptStationary.hpp>`

```c++
#include <amdis/AMDiS.hpp>                              // Include essential headers
#include <amdis/AdaptStationary.hpp>
int main(int argc, char** argv)
{
  //  make problem and add operators/constraints as above

AdaptInfo adaptInfo{"adapt"};
  AdaptStationary adaptStat{"adapt", prob, adaptInfo};

adaptStat.adapt();                                    // Replaces the part below
  //  prob.assemble(adaptInfo);
  //  prob.solve(adaptInfo);
  //  prob.writeFiles(adaptInfo);
}
```
---

# Solvers

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

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

---

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

- template code allows different matrix/vector implementations
- solvers can be selected at runtime using a parameter file
  <br>.icon-rechts[] many solver types available just by changing one line
--

- backend chosen when compiling AMDiS (default ISTL)

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

- list of options in the [online manual](https://amdis.readthedocs.io/en/latest/reference/Initfile/#solvers-and-preconditioners)

---

# Writing Files

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

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

---

# Writing Files

- write files automatically after each time step or given intervals
- choose options for each basis node individually ([see Initfile](https://amdis.readthedocs.io/en/latest/reference/Initfile/#filewriters))
- view `VTK` files using [ParaView](https://www.paraview.org/)

```bash
example->output[0]->format: vtk
example->output[0]->filename: example
example->output[0]->name: u
example->output[0]->output directory: ./output
```

--
### Backup and Restore
- can also generate backups to pause and restart computations

```bash
example->output->format: backup
example->restore->grid: FILENAME
example->restore->solution: FILENAME

void ProblemStat::restore(Flag); // instead of initialize(Flag)
```

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

# See the example in action!

---

---

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

int main(int argc, char** argv)
{
  using namespace AMDiS;
  using namespace Dune::Functions::BasisFactory;

Environment env{argc, argv};

* auto f = [](auto x) { return -1.0; };
* auto b = [](auto x) { return x[0] < 1e-8 || x[1] < 1e-8; };
* auto g = [](auto x) { return 0.0; };

Dune::YaspGrid<2> grid{ {1.0,1.0}, {8,8} };
* ProblemStat prob{"my-first-project", grid, lagrange<2>()};
  prob.initialize(INIT_ALL);

prob.addMatrixOperator(sot(1.0));
* prob.addVectorOperator(zot([&](auto x) { return f(x); }, 6));
  prob.addDirichletBC([&](auto x) { return b(x); },
                      [&](auto x) { return g(x); });

AdaptInfo adaptInfo{"adapt"};
  prob.assemble(adaptInfo);
  prob.solve(adaptInfo);
  prob.writeFiles(adaptInfo);
}
```

---

# 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\otimes 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 and DiscreteFunction

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

---

# DOFVector and DiscreteFunction

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

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

---

# DOFVector and DiscreteFunction

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

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

---

# DOFVector and DiscreteFunction

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

// getting the composite discrete-function u
auto u = valueOf(U);
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};
```

- Extract the type of a DOFVector:

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

---

# DOFVector and DiscreteFunction

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

// getting the composite discrete-function u
auto u = valueOf(U);
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};
```

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

---

# 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:

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

auto opL = makeOperator(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};
```

---

# Grid Construction

- explicit using the grid class constructor

- no common interface
  - pass grid to `ProblemStat` instance via constructor
  - only use the grid wrapper stored in `ProblemStat` using `prob.grid()`

```c++
  Dune::YaspGrid<2> baseGrid({1.0, 1.0}, {4, 4});
  ProblemStat prob("example", baseGrid, basisFactory);
  auto& grid = *prob.grid();
  // do something with grid
  ```

--
  - advanced: use the grid wrapper without `ProblemStat`

```c++
  Dune::YaspGrid<2> baseGrid({1.0, 1.0}, {4, 4});
  AdaptiveGrid grid(baseGrid);
  // do something with grid
  ```

---

# Grid Construction

- using the `MeshCreator` helper class

```c++
  using Grid = Dune::YaspGrid<2>;
  auto gridPtr = MeshCreator<Grid>{"meshName"}.create();
  auto& grid = *gridPtr;
  // do something with grid
  ```

--
  - creates grid from macro file using `meshName->macro file name`
  - if not set creates structured grid

```bash
  meshName->macro file name: example_macro.2d

meshName->structured: cube        % cube, simplex
  meshName->min corner: 0.0 0.0     % cube domain with corner
  meshName->max corner: 1.0 1.0     %   coordinates given here
  meshName->num cells: 4 4          % elements in each direction
  meshName->global refinements: 2   % initial refinement steps
  meshName->load balance: 0         % perform initial load balancing
  ```

---

# Grid Construction

- create automatically by `ProblemStat<Traits>`
- `Traits` declares grid and basis: see `amdis/ProblemStatTraits.hpp`

```c++
// basis composed of lagrange bases of different degrees
template <class Grid, int... degrees>
struct LagrangeBasis;

// LagrangeBasis using Dune::YaspGrid<dim>
template <int dim, int... degrees>
struct YaspGridBasis;

// taylor-hood basis with pressure degree k
template <class Grid, int k = 1>
struct TaylorHoodBasis;
```

- Example:
```c++
// 2d YaspGrid with power<3>(lagrange<1>()) basis
ProblemStat<YaspGridBasis<2, 1,1,1>> prob("prob");
```

---

# Traversing Grid Elements

- see [Oliver Sander's book on Dune](https://cloudstore.zih.tu-dresden.de/index.php/s/NBsqPs7J9cEdmwA) for in-depth explanation

```c++
auto const& gridView = grid.leafGridView(); // or
// auto const& gridView = grid.levelGridView(2);

for (auto const& element : elements(gridView))
{
  // do something with element
}
```

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

---

# Element Geometries

- `geometry()`: obtain mapping \( F: T_{ref} \to T \) from reference to real element
- exports coordinate types in \( T_{ref} \text{ and } T \)

```c++
auto geometry = element.geometry();
using LocalCoordinate = typename decltype(geometry)::LocalCoordinate;
using GlobalCoordinate = typename decltype(geometry)::GlobalCoordinate;
```

---

# Element Geometries

- `geometry()`: obtain mapping \( F: T_{ref} \to T \) from reference to real element
- exports coordinate types in \( T_{ref} \text{ and } T \)
- `father()`: get father element (element one level above in the hierarchy)
- `geometryInFather()`: obtain mapping from reference element to father's reference element

---

# Element Geometries

- `global(x)`: evaluates the mapping \( F \)
- `local(x)`: evaluates the inverse \( F^{-1} \)
- `integrationElement(x)`: returns \( \sqrt{\det(\nabla F^T \nabla F)} \)
- `jacobianInverseTransposed(x)`: returns \( (\nabla F)^{-T} \)
- `jacobianTransposed(x)`: returns \( (\nabla F)^{T} \)
- `volume()`: volume of the image of \( F \)
- `center()`: coordinates of the centerpoint of the image of \( F \)
- `corner(i)`: coordinates of the i-th corner of the image of \( F \)

<br>.icon-rechts[] can compute integrals over elements

---

# Intersections

- sometimes integrals on the boundary or the interface between two elements are needed .icon-rechts[] `Intersection` class

```c++
for (auto const& is : intersections(gridView, element))
{
  // do something with is
}
```

- can access elements on either side
- provide normal to the intersection for computations
--

- example: iterate over all boundary intersections:

---

# Boundary Intersections in AMDiS

- available via `BoundaryManager` class

```c++
auto& boundaryManager = *prob.boundaryManager();
```

- set boundary ids [left,right, front,back, bottom,top] for cube domains

```c++
void setBoxBoundary(std::array<BoundaryType, 2*dow> const& ids);
```

- use predicate-id-pair or indicator function

```c++
template <class Indicator>
void setIndicator(Indicator const& indicator);

template <class Predicate>
void setPredicate(Predicate const& pred, BoundaryType id);
```

---

# Boundary Intersections in AMDiS

- retrieve id using intersections

```c++
template <class Intersection>
BoundaryType boundaryId(Intersection const& intersection) const;
```

```c++
for (auto const& element : elements(gridView))
  if (element.hasBoundaryIntersections())
    for (auto const& intersection : intersections(gridView, element))
      if (intersection.boundary())
      {
        auto id = prob.boundaryManager()->boundaryId(intersection);
        // do something specific to id
      }
```

---

# Global Refinement

- already known: initial global refine
- `grid.globalRefine(count)` or `prob.globalRefine(1)`
- grid elements change and with them the local basis functions
--

<br>.icon-rechts[] basis must be updated and user data must be interpolated onto new elements
--

- in DUNE: manually update the basis and everything depending on it
- in AMDiS: update happens automatically when using the grid and basis wrappers
  - global basis is *notified* by the grid when `globalRefine` is called
  - `DOFVector` get interpolated after the basis is updated
  - user classes can *observe* changes in other classes by deriving from the `Observer<Event>` class
  - **Observer - Notifier** pattern

---

# Local Adaptivity

- refinement can be controlled by marks on the grid elements

```c++
auto marker = GridFunctionMarker{"marker", Grid, Expression};
prob.addMarker(marker);
```

- `Expression` maps center point of elements to the requested refinement level
- if requested level is higher than element level: mark for refinement
- if requested level is lower than element level: mark for coarsening
--

- start adaptation process with `prob.adaptGrid(adaptInfo)`
<br>.icon-rechts[] coarsen and refine in one step, maximum one level difference
<br>.icon-rechts[] automatic update of associated data just like when using `globalRefine`

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

# Example: Moving Stripe

---

---

```c++
#include <amdis/AMDiS.hpp>
#include <amdis/LocalOperators.hpp>
#include <dune/alugrid/grid.hh>
int main(int argc, char** argv)
{
  using namespace AMDiS;
  Environment env{argc, argv};
  using Grid = Dune::ALUGrid<2, 2, Dune::simplex, Dune::conforming>;
  ProblemStat<LagrangeBasis<Grid, 2>> prob{"my-first-project"};
  prob.initialize(INIT_ALL);
  auto f = [](auto x) { return x[0] * x[1]; };
  double poi = 0.0;
  // Get DOFVector and interpolate f
  // Refine up to level 10 in the stripe [poi-0.25, poi+0.25] x [0.0, 1.0]
  auto g = [&](auto x) { return ?; };
  auto marker = GridFunctionMarker{"marker", prob.grid(), g};
  prob.addMarker(marker);
  AdaptInfo adaptInfo("adapt");
  for (; poi <= 2.0; poi += 0.25) {
    for (std::size_t i = 0; i < 8; ++i) {
      // mark the grid and call adapt
    }
    prob.writeFiles(adaptInfo);
  }
}
```

---

```bash
mesh->global refinements: 3
mesh->structured: simplex
mesh->max corner: 2.0 1.0
mesh->num cells: 4 2

my-first-project->mesh: mesh
my-first-project->solver: default
my-first-project->solver->relative tolerance: 1.e-6
my-first-project->solver->info: 1

my-first-project->output->format: vtk
my-first-project->output->directory: .
my-first-project->output->filename: my_first_project
my-first-project->output->name: solution
my-first-project->output->ParaView mode: 1
my-first-project->output->animation: 1
```

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

# see the example in action

---

---

```c++
#include <amdis/AMDiS.hpp>
#include <amdis/LocalOperators.hpp>
#include <dune/alugrid/grid.hh>
int main(int argc, char** argv)
{
  using namespace AMDiS;
  Environment env{argc, argv};
  using Grid = Dune::ALUGrid<2, 2, Dune::simplex, Dune::conforming>;
  using Traits = LagrangeBasis<Grid, 2>;
  ProblemStat<Traits> prob{"my-first-project"};
  prob.initialize(INIT_ALL);
  auto f = [](auto x) { return x[0] * x[1]; };
  double poi = 0.0;
  prob.solution() << f;
  auto g = [&](auto x) { return (x[0] <= poi+0.25 && x[0] >= poi-0.25) ? 10 : 3; };
  auto marker = GridFunctionMarker{"marker", prob.grid(), g};
  prob.addMarker(marker);
  AdaptInfo adaptInfo("adapt");
  for (; poi <= 2.0; poi += 0.25)
  {
    for (std::size_t i = 0; i < 8; ++i)
    {
      prob.markElements(adaptInfo);
      prob.adaptGrid(adaptInfo);
    }
    prob.writeFiles(adaptInfo);
  }
}
```

---

# Linear Algebra Backends and Solvers

- backend chosen when compiling AMDiS (default ISTL)

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

- list of options in the [online manual](https://amdis.readthedocs.io/en/latest/reference/Initfile/#solvers-and-preconditioners)

---
# Linear Algebra Backends and Solvers

- Backends define their own types for matrices and vectors
- Each backend has its own set of solvers and preconditioners
- Solvers are named like in the backend libraries
- 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

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

```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
  ...
```

---

# Parallel Grids and Parallel Solvers
- AMDiS is by default aware of parallelization
- Parallel grid __and__ linear algebra backend required
- Grids with *ghost* elements

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

---

# **Parallel Grids** and Parallel Solvers
- AMDiS is by default aware of parallelization
- Parallel grid __and__ linear algebra backend required
- Grids with *overlapping* elements

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

---

# **Parallel Grids** and Parallel Solvers
- 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).

---

# Operators and LocalOperators

- Assemble element-matrix / element-vector on each element
- Attached to node in basis tree
- `Operator` provides a `LocalOperator` by `localOperator(Operator)`

`Operator` interface to implement:
```c++
void update(GridView const&);  // Update the operator data

// Transform an operator into a local-operator
friend LocalOperator localOperator(Operator const&);
```

---

# Operators and LocalOperators

- Assemble element-matrix / element-vector on each element
- Attached to node in basis tree
- `Operator` provides a `LocalOperator` by `localOperator(Operator)`

`LocalOperator` interface to implement:
```c++
void bind(Element const&);  // Bind to the currently visited element in grid traversal
void unbind();              // Might free memory and unbinds other data from the element

// Assemble a local element matrix on the element that is bound
void assemble(ContextGeo const&, RowNode const&, ColNode const&, ElementMatrix&) const;

// or:
// Assemble a local element vector on the element that is bound
void assemble(ContextGeo const&, Node const&, ElementVector&) const;
```

---

# Operators and LocalOperators

```c++
class ZeroOrderTest
{
public:
  template <class Element>
  void bind(Element const& element)
  {
    lf.bind(element);
  }

void unbind()
  {
    lf.unbind();
  }

// ...

private:
  LocalFunction lf;
};
```

---

# Operators and LocalOperators

```c++
template <class EG, class Node, class Vec>
void assemble(EG const& elementGeo, Node const& node, Vec& elementVector) const
{
  auto const& quad = Dune::QuadratureRules<double,EG::dimension>::rule(10);
  for (auto const& qp : quad)
  {
    // The multiplicative factor in the integral transformation formula
    auto dx = elementGeo.integrationElement(qp.position()) * qp.weight();

// Evaluation of local function at quadrature point
    auto fAtQP = lf(qp.position());

auto const& shapeValues = node.localBasisValuesAt(qp.position());
    for (std::size_t i = 0; i < node.size(); ++i) {
      auto local_i = node.localIndex(i);
      elementVector[local_i] += fAtQP * shapeValues[i] * dx;
    }
  }
}

```

---

# Summary and Outlook

</div>
</div>

---

# Summary and Outlook

<div class="pure-g">
  <div class="pure-u-1-2"><h2>Summary</h2>
  <ul>
    <li>DUNE base - AMDiS frontend</li>
    <li>Flexible choice of grids</li>
    <li>Flexible choice of finite-element space</li>
    <li>Various linear-algebra backends</li>
    <li>Expression templates for grid-functions</li>
    <li>Operators applicable to basis node</li>
  </ul>
  </div>
  <div class="pure-u-1-2"><h2>ToDo</h2>
  <ul>
    <li>Better handling of Constraints</li>
    <li>Load balancing with automatic data communication</li>
    <li>Assembler/LinearForm/BiLinearForm get restructured</li>
    <li>Support for blocked data-structures</li>
  </ul>
  </div>
</div>

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

# Find the code on

.font-xl[
[gitlab.mn.tu-dresden.de/amdis/amdis-core](https://gitlab.mn.tu-dresden.de/amdis/amdis-core)
]

# Documentation