New generated Stokes operators

This issue shall document and discuss the introduction of the Stokes composite operators that are built from the new generated operators. Also point out if I made a mistake and feel free to correct this description.

List of Stokes operators

It seems(?) that we are particularly interested in three types of Stokes operators, particularly differing in the momentum terms.

Let

K = \begin{bmatrix} 
A & B^T \\
B & 0
\end{bmatrix}

be the discrete Stokes operator for a \mathbb{P}_2-\mathbb{P}_1 Taylor-Hood discretization.

The operators differ in the way the A block is defined. Let's summarize.

Proposed name (prefix) Strong form of momentum term ("velocity part") Weak form that defines A
P2P1StokesConstant - \Delta u \displaystyle\int_\Omega \nabla u : \nabla v
P2P1StokesEpsilon \displaystyle- \nabla \cdot \Big(2 \mu \varepsilon(u)\Big) \displaystyle\int_\Omega 2 \mu \Big(\varepsilon(u) : \varepsilon(v)\Big)
P2P1StokesFull \displaystyle - \nabla \cdot \Big( 2 \mu \varepsilon(u) \Big) - \frac{2}{\mathrm{dim}} \nabla (\nabla \cdot u) \displaystyle\int_\Omega 2 \mu \Big(\varepsilon(u) : \varepsilon(v)\Big) - \frac{2}{\mathrm{dim}} (\nabla \cdot u) \cdot (\nabla \cdot v)

where \varepsilon(w) = \frac{1}{2} (\nabla w + (\nabla w)^T).

\mu = \mu(x) will eventually be a scalar finite element function (in \mathbb{P}_2) that is supplied by the user.

Possible file structure

I propose the following structure - open to suggestions.

hyteg/src/
- hyteg/
- terraneo/
- ...
- hyteg_operators/                 # submodule with generated operators
- hyteg_operators_composites/      # new module
  - stokes/
    - viscousblock/
      - P2ViscousBlockLaplaceOperator.*pp
      - P2ViscousBlockEpsilonOperator.*pp
      - P2ViscousBlockFullOperator.*pp
    - divergence/
      - P2VectortoP1DivergenceOperator.*pp
    - gradient/
      - P1toP2VectorGradientOperator.*pp
    - P2P1StokesConstantOperator.*pp       
    - P2P1StokesEpsilonOperator.*pp
    - P2P1StokesFullOperator.*pp       
  - ...

Blending

For each relevant blending map, additional operators will be added to the respective files. For instance

// P2P1StokesEpsilonOperator.hpp

class P2P1StokesEpsilonOperator { ... };
class P2P1StokesEpsilonIcosahedralShellMapOperator { ... };

Some notes

  • I am not sure yet which quadrature rules to choose for the blending cases. For now I chose rules that integrate degree 3 polynomials exactly for blending, and degree 2 polynomials without blending (the latter should be sufficient to maintain convergence order(?)). But we may want to discuss this.
  • The new composites shall only use/depend on the newly generated operators and never introduce any of the existing constant stencil ops, elementwise ops etc. If something is missing, we should adapt the HFG.
  • Lumped and diagonal operators are already in the making. Those are needed for all sorts of preconditioners.
  • I currently only consider Taylor-Hood-style operators. Stabilized \mathbb{P}_1-\mathbb{P}_1 operators should be straightforward to generate, too.
  • In the future, we want to enable the HFG to generate the vector-operators directly such that they do not have to be composed manually. This will likely remove the need for all operators in the subdirectories viscousblock, gradient, divergence, etc. But for starters, we will probably go with the composed versions.

Checklist

Edited by Nils Kohl