Cleanup some of the poisson example sections

- how difficult setup is in practice,

- how to move from first compile to first PDE solve.

!!! warning "Preliminary draft"

    This tutorial is currently in a very preliminary draft state.

    Content, structure, and examples may still change significantly.

## Who this is for

- Beginners who know C++ and basic FEM and want a realistic first contact with DUNE.

## Local basis from `dune-localfunctions`

Instead of hardcoding basis functions, we use:

We use a predefined Q1 Lagrange element from `dune-localfunctions`:

```cpp

using LocalFiniteElement = Dune::LagrangeCubeLocalFiniteElement<double,double,2,1>;

const auto& basis = lfe.localBasis();

```

Useful type information from the local basis traits:

For this tutorial we only need `localBasis()` to evaluate shape functions and derivatives.

- `RangeType` is `FieldVector<T,1>` for scalar-valued basis functions.

- `JacobianType` is `FieldMatrix<T,1,m>`, where `m` is the local element dimension.

  For 2D elements this is `FieldMatrix<T,1,2>`.

Useful trait types:

At each quadrature point on each cell:

- `RangeType = FieldVector<T,1>` for scalar basis values,

- `JacobianType = FieldMatrix<T,1,m>` for derivatives in local coordinates (`m` = local dimension, so in 2D: `FieldMatrix<T,1,2>`).

- `basis.evaluateFunction(local, phi)` provides $\varphi_i$,

- `basis.evaluateJacobian(local, jacobian)` provides reference gradients,

- `geometry.jacobianInverseTransposed(local)` maps gradients to global coordinates.

At each quadrature point $\hat x$:

### Gradient transformation to physical coordinates

- `basis.evaluateFunction(local, phi)` gives $\{\hat\varphi_i(\hat x)\}_i$,

- `basis.evaluateJacobian(local, jacobian)` gives $\{\nabla_{\hat x}\hat\varphi_i(\hat x)\}_i$.

Let $F_K : \hat K \to K$ be the element map with Jacobian $J_K$.

For each local shape function, the gradient transformation is

In code this is exactly the `mv(...)` call with the inverse-transposed Jacobian:

```cpp

std::vector<LocalFiniteElement::Traits::LocalBasisType::Traits::JacobianType> jacobian;

using Traits = typename LocalFiniteElement::Traits::LocalBasisType::Traits;

using JacobianType = typename Traits::JacobianType;

std::vector<JacobianType> jacobian;

basis.evaluateJacobian(local, jacobian); // reference gradients

const auto jacInvT = geo.jacobianInverseTransposed(local);

std::array<Dune::FieldVector<double,2>,4> grad;

for (int i = 0; i < 4; ++i)

std::array<Dune::FieldVector<double,2>,numElementDofs> grad;

for (int i = 0; i < numElementDofs; ++i)

  jacInvT.mv(jacobian[i][0], grad[i]);  // grad[i] = J^{-T} * gradHat[i]

```

## Assembly loop

With `dx = qp.weight() * geo.integrationElement(local)`, we add

### Element and quadrature traversal

Assembly iterates over all elements and, inside each element, over all quadrature points:

```cpp

for (const auto& e : Dune::elements(gv)) {

  const auto v = elementDofs(gv, e); // local vertex dofs -> global dofs

  const auto geo = e.geometry();

  const auto quadrature = Dune::QuadratureRules<double,2>::rule(geo.type(), 2);

  for (const auto& qp : quadrature) {

    const auto local = qp.position();

    basis.evaluateFunction(local, phi);

    basis.evaluateJacobian(local, jacobian);

    const auto jacInvT = geo.jacobianInverseTransposed(local);

    const double dx = qp.weight() * geo.integrationElement(local);

    for (int i = 0; i < numElementDofs; ++i)

      jacInvT.mv(jacobian[i][0], grad[i]);

    // local-to-global accumulation comes here

  }

}

```

`dx` is the quadrature contribution to the integration measure on the physical element.

### Local-to-global accumulation

Using the global DoF indices `v`, local contributions are scattered into global `A` and `b`:

```cpp

for (int i = 0; i < numElementDofs; ++i) {

  for (int j = 0; j < numElementDofs; ++j)

    A[v[i]][v[j]] += dx * (grad[i] * grad[j]);

  b[v[i]] += dx * phi[i][0];

}

```

for all local DoFs `i,j`.

Here `grad[i] * grad[j]` is the Euclidean dot product

$\nabla\varphi_i \cdot \nabla\varphi_j$.

# Possible extensions

Ideas for extending this Poisson tutorial:

- unstructured grid, maybe read from a file. This relates to the grid I/O tutorial and uses other grid managers than `YaspGrid`.

- using higher-order finite elements, e.g. `Q2` or `P2`, with additional DoFs on edges (or even in the element interior). This requires a mapper, as introduced in the corresponding grid data tutorial.

- non-homogeneous boundary conditions, e.g. a subset of the boundary with Dirichlet conditions and the rest Neumann, possibly using boundary intersection normal vectors.

- build a local element matrix and local element vector before scattering into the global matrix.

- thread-parallel assembly.

!!! abstract "Learning goals"

    After this chapter, you should be able to:

    - state the Poisson model and weak form for homogeneous Dirichlet data,

    - set up grid, DoF numbering, sparse matrix, and rhs vector,

    - assemble the system using `dune-localfunctions` local basis interfaces,

    - solve with `dune-istl` and export a `.vtu` file for ParaView.

## Model problem

We solve

We will solve the Poisson equation in the unit square domain with homogeneous Dirichlet boundary conditions,

$$

-\Delta u = 1 \quad\text{in }\Omega=(0,1)^2,

\qquad

u = 0 \quad\text{on }\partial\Omega.

u = 0 \quad\text{on }\partial\Omega,

$$

Weak form: find $u\in H_0^1(\Omega)$ such that

using the finite element method. Therefore, we transform the PDE into it's weak form: <br>find $u\in H_0^1(\Omega)$ such that

$$

a(u,v)=\ell(v) \quad \forall v\in H_0^1(\Omega),

\ell(v)=\int_\Omega v\,dx.

$$

We use first-order conforming Lagrange elements (`Q1` on a structured square grid, one DoF per vertex).

We use first-order conforming Lagrange elements (`Q1` on a structured square grid, one DoF per vertex) for the discretization of the function space $H^1$ and incorporate the Dirichlet boundary condition by a modification of the algebraic system.

## Three-step path

1. [Grid and linear algebra setup](../poisson/setup.md)

2. [Assembly and boundary conditions](../poisson/assembly-bc.md)

3. [Solve and VTK output](../poisson/solve-vtk.md)

## Summary and next steps

You now have the Poisson overview and implementation roadmap.

Continue with [Step 1: Grid and linear algebra setup](../poisson/setup.md).

# Step 1: Grid and linear algebra setup

This step creates the grid, global DoF numbering, and sparse matrix layout.

This step creates the computational grid, global DoF numbering, and sparse matrix layout.

!!! abstract "Learning goals"

    After this page, you should be able to:

    - create a unit-square `YaspGrid` with concise constructor syntax,

    - define one global DoF per vertex using `indexSet`,

    - initialize `BCRSMatrix<double>` and `BlockVector<double>`.

    - create a unit-square `YaspGrid`,

    - define one global DoF per vertex using an index set,

    - initialize `BCRSMatrix` and `BlockVector`.

## Grid setup

## Vertex DoF numbering

For first-order conforming elements, one DoF is attached to each vertex. Thus, we get as total number of degrees of freedom the number of vertices in the grid:

First-order conforming elements have one DoF attached to each vertex. Thus, we get as total number of degrees of freedom the number of vertices in the grid:

```cpp

const auto numDofs = gv.indexSet().size(2);

```

In each element of the grid, the local corner vertices are associated to a global vertex index. This can be retrieved from the grid's index set, by collecting the vertex sub indices for each element. For performance reasons, it is good to collect this index mapping of local to global vertex indices:

In each element of the grid, the local corner vertices are associated to a global vertex index. This can be retrieved from the grid's index set, by collecting the vertex sub-indices for each element. For performance reasons, it is good to collect this index mapping of local to global vertex indices into a container,

```cpp

template<class GridView, class Element>

!!! note

    If the number of corners or the number of DoFs per element is a statically known constant, one could store this mapping in a `std::array`, instead of a dynamic `std::vector`. This requires to know the element type of the grid and its dimension.

    For example, a 2d `YaspGrid` allows only square elements, and a first order Lagrange local finite-element has one degree of freedom per vertex. This makes $4=2^{\text{dim}}$ DoFs per Element.

## Sparse matrix pattern

    ```c++

    constexpr int numElementDofs = 1<<dim;

    std::array<std::size_t,numElementDofs> dofs;

    for (int localVertex = 0; localVertex < numElementDofs; ++localVertex)

      dofs[localVertex] = indexSet.subIndex(e, localVertex, dim);

    ```

Before assembly, we collect all nonzero couplings with `MatrixIndexSet`:

```cpp

Dune::MatrixIndexSet pattern(numDofs, numDofs);

Dune::BCRSMatrix<double> A;

pattern.exportIdx(A);

A = 0.0;

```

### Algorithm sketch: setup matrix pattern

## Sparse matrix pattern

Before assembly, we collect all nonzero couplings with `MatrixIndexSet`:

```cpp

Dune::MatrixIndexSet pattern(numDofs, numDofs);

This pattern contains all row/column pairs that can receive element contributions.

For a vertex-based first-order method, all local vertex DoFs on one element couple with each other, so we insert the full local `nCorners x nCorners` block.

The rhs uses scalar blocks too:

The rhs will be assembled into a vector with `numDofs` entries,

```cpp

Dune::BlockVector<double> b(numDofs);