and external volume force $`\mathbf{f}`$, w.r.t. to boundary conditions $`\mathbf{u}=\mathbf{g}`$ on $`\partial\Omega`$.
For an overview, see also [Wikipedia](https://en.wikipedia.org/wiki/Navier%E2%80%93Stokes_equations).
Discretization
--------------
### Weak formulation
In a weak variational formulation, we try to find $`\mathbf{u}(t,\cdot)\in \mathbf{V}_\mathbf{g}:=\{\mathbf{v}\in H^1(\Omega)^d\,:\, \operatorname{tr}_{\partial\Omega}\mathbf{v} = \mathbf{g}\}`$ and $`p(t,\cdot)\in Q:=L_0^2(\Omega)`$, such that
for $`t\in[0,T]`$, given $`\mathbf{u}(0,\cdot)\equiv\mathbf{u}^0`$.
The pair $`\big(\mathbf{u}(t,\cdot), p(t,\cdot)\big)`$ is thus in the product space $`\mathbf{V}_\mathbf{g}\times Q`$ where $`\mathbf{V}_\mathbf{g}`$ is also a product of $`H^1`$ spaces, i.e., $`\mathbf{V}_\mathbf{g}=V_{g_0}\times V_{g_1}\times V_{g_2}\simeq [V_{g_i}]^d`$. Later, this space will be approximated with a Taylor-Hood space of piecewise quadratic and linear functions.
### Time-discretization
A linearized time-discretization with a semi-implicit backward Euler scheme then reads:
For $`k=0,1,2,\ldots, N`$ find $`\mathbf{u}^{k+1}\in\mathbf{V}_\mathbf{g},\, p\in Q`$, s.t.
with $`0=t^0<t^1<\ldots<t^N=T`$, timestep $`\tau=t^{k+1} - t^k`$, and given initial velocity $`\mathbf{u}^0`$.
### Finite-Element formulation
Let $`\Omega`$ be partitioned into quasi-uniform elements by a conforming triangulation $`\mathcal{T}_h`$. The functions $`\mathbf{u}`$ and $`p`$ are approximated by functions from the discrete functionspace $`\mathbb{T}_h := [\mathbb{V}_{2,\,g_i}]^d\times \mathbb{V}_1`$ with $`\mathbb{V}_{p} := \{v\in C^0(\Omega)\,:\, v|_T\in\mathbb{P}_p(T),\, \forall T\in\mathcal{T}_h\}`$ and $`\mathbb{V}_{p,\,g} := \{v\in \mathbb{V}_p\,:\,\operatorname{tr}_{\partial\Omega} v=g\}`$, where $`\mathbb{P}_p`$ is the space of polynomials of order at most $`p`$. The space $`\mathbb{T}_h`$ will be denoted as Taylor-Hood space.
Let `Grid grid` be the Dune grid type representing the triangulation $`\mathcal{T}_h`$ of the domain $`\Omega`$, then
the basis of the Taylor-Hood space can be written as a composition of Lagrange bases:
```c++
usingnamespaceDune::Functions::BasisFactory;
autobasis=composite(
power<Grid::dimensionworld>(
lagrange<2>(),flatInterleaved()),
lagrange<1>(),flatLexicographic());
```
At this point, any Dirichlet boundary conditions are ignored in the definition of the bases. We
build the Taylor-Hood basis as a composition of a product of $`d\times\mathbb{V}_2`$ bases for
the velocity components and a $`\mathbb{V}_1`$ basis for the pressure. We have to use
`composite(...)` to combine the different types for velocity and pressure space, while we
can use `power<d>(...)` for the `d` velocity component spaces of the same type.
For efficiency, we use flat indexing strategies `flatInterleaved` and `flatLexicographic`
to build a global continuous numbering of the basis functions.
See the tutorial [Grids and Discrete Functions](../tutorials/grids-and-discretefunctions.md)
for more details.
Implementation
--------------
We split the implementation into three parts, 1. the time derivative, 2. the stokes
operator and 3. any external forces. For addressing the different components of
the Taylor-Hood basis, we introduce the treepaths `_v = Dune::Indices::_0` and
`_p = Dune::Indices::_1` for velocity and pressure, respectively.
### Problem framework
We have an in-stationary problem consisting of a sequence of stationary equations
for each timestep, thus we combine a [`ProblemInstat`](../reference/Problem.md#class-probleminstat)
with a [`ProblemStat`](../reference/Problem.md#class-problemstat) .
```c++
usingnamespaceAMDiS;
ProblemStatprob{"stokes",grid,basis};
prob.initialize(INIT_ALL);
ProblemInstatprobInstat{"stokes",prob};
probInstat.initialize(INIT_UH_OLD);
```
For convenience, we use class template argument deduction here.
### Time derivative
For simplicity of this example, we use a backward Euler time stepping and a linearization
