# Feel++

## Function Spaces

Prerequisites

The prerequisites are

### Notations

We now turn to the next crucial mathematical ingredient: the function space, whose definition depends on \Omega_h - or more precisely its partitioning \mathcal{T}_h - and the choice of basis function. Function spaces in Feel++ follow the same definition and Feel++ provides support for continuous and discontinuous Galerkin methods and in particular approximations in L^2, H^1-conforming and H^1-nonconforming, H^2, H(\mathrm{div}) and H(\mathrm{curl})[^1].

We introduce the following spaces

\begin{aligned} \mathbb{W}_h &= \{v_h \in L^2(\Omega_h): \ \forall K \in \mathcal{T}_h, v_h|_K \in \mathbb{P}_K\},\\ \mathbb{V}_h &= \mathbb{W}_h \cap C^0(\Omega_h)= \{ v_h \in \mathbb{W}_h: \ \forall F \in \mathcal{F}^i_h\ [ v_h ]_F = 0\}\\ \mathbb{H}_h &= \mathbb{W}_h \cap C^1(\Omega_h)= \{ v_h \in \mathbb{W}_h: \ \forall F \in \mathcal{F}^i_h\ [ v_h ]_F = [ \nabla v_h ]_F = 0\}\\ \mathbb{C}\mathbb{R}_h &= \{ v_h \in L^2(\Omega_h):\ \forall K \in \mathcal{T}_h, v_h|_K \in \mathbb{P}_1; \forall F \in \mathcal{F}^i_h\ \int_F [ v_h ] = 0 \}\\ \mathbb{R}{a}\mathbb{T}{u}_h &= \{ v_h \in L^2(\Omega_h):\ \forall K \in \mathcal{T}_h, v_h|_K \in \mathrm{Span}\{1,x,y,x^2-y^2\}; \forall F \in \mathcal{F}^i_h\ \int_F [ v_h ] = 0 \}\\ \mathbb{R}\mathbb{T}_h &= \{\mathbf{v}_h \in [L^2(\Omega_h)]^d:\ \forall K \in \mathcal{T}_h, v_h|_K \in \mathbb{R}\mathbb{T}_k; \forall F \in \mathcal{F}^i_h\ [{\mathbf{v}_h \cdot \mathrm{n}}]_F = 0 \}\\ \mathbb{N}_h &= \{\mathbf{v}_h \in [L^2(\Omega_h)]^d:\ \forall K \in \mathcal{T}_h, v_h|_K \in \mathbb{N}_k; \forall F \in \mathcal{F}^i_h\ [{\mathbf{v}_h \times \mathrm{n}}]_F = 0 \} \end{aligned}

where \mathbb{R}\mathbb{T}_k and \mathbb{N}_k are respectively the Raviart-Thomas and Nédélec finite elements of degree k.

The Legrendre and Dubiner basis yield implicitely discontinuous approximations, the Legendre and Dubiner boundary adapted basis, see~\cite MR1696933, are designed to handle continuous approximations whereas the Lagrange basis can yield either discontinuous or continuous (default behavior) approximations. \mathbb{R}\mathbb{T}_h and \mathbb{N}_h are implicitely spaces of vectorial functions \mathbf{f} such that \mathbf{f}: \Omega_h \subset \mathbb{R}^d \mapsto \mathbb{R}^d. As to the other basis functions, i.e. Lagrange, Legrendre, Dubiner, etc., they are parametrized by their values namely `Scalar`, `Vectorial` or `Matricial.`

 Products of function spaces must be supported. This is very powerful to describe complex multiphysics problems when coupled with operators, functionals and forms described in the next section. Extracting subspaces or component spaces are part of the interface.

#### Function Spaces

Function spaces support is provided by the `FunctionSpace` class

The `FunctionSpace` class

• constructs the table of degrees of freedom which maps local (elementwise) degrees of freedom to the global ones with respect to the geometrical entities,

• embeds the definition of the elements of the function space allowing for a tight coupling between the elements and their function spaces,

• stores an interpolation data structure (e.g. region tree) for rapid localisation of point sets (determining in which element they reside).

 C++ Function C++ Type Function Space [1] `Pch(mesh)` `Pch_type` P^N_{c,h} `Pchv(mesh)` `Pchv_type` [P^N_{c,h}^d] `Pdh(mesh)` `Pdh_type` P^N_{d,h} `Pdhv(mesh)` `Pdhv_type` [P^N_{d,h}^d] `THch(mesh)` `THch_type` [P^{N+1}_{c,h}^d \times P^N_{c,h}] `Dh(mesh)` `Dh_type` \mathbb{R}\mathbb{T}_h `Ned1h(mesh)` `Ned1h_type` \mathbb{N}_h

[[[1]]]: see Notations for the function spaces definitions.

Here are some examples how to define function spaces with Lagrange basis functions.

``````#include <feel/feeldiscr/pch.hpp>

// Mesh with triangles
using MeshType = Mesh<Simplex<2>>;
// Space spanned by P3 Lagrange finite element
FunctionSpace<MeshType,bases<Lagrange<3>>> Xh;
// is equivalent to (they are the same type)
Pch_type<MeshType,3> Xh;

// using the auto keyword
MeshType mesh = loadMesh( _mesh=new MeshType );
auto Xh = Pch<3>( mesh );
// is equivalent to
auto Xh = FunctionSpace<MeshType,bases<Lagrange<3>>>::New( mesh );
auto Xh = Pch_type<MeshType,3>::New( mesh );``````

#### Functions

 One important feature in `FunctionSpace` is that it embeds the definition of element which allows for the strict definition of an `Element` of a `FunctionSpace` and thus ensures the correctness of the code.

An element has its representation as a vector, also in the case of product of multiple spaces.

``````#include <feel/feeldiscr/pch.hpp>

// Mesh with triangles
using MeshType = Mesh<Simplex<2>>;
auto mesh = loadMesh( _mesh=new MeshType );

// define P3 Lagrange finite element space
auto P3ch = Pch<3>(mesh);

// definie an element from P3ch, initialized to 0
auto u = P3ch.element();
// definie an element from P3ch, initialized to x^2+y^2
auto v = P3ch.element(Px()*Px()+Py()*Py());``````

#### Components

``````FunctionSpace<Mesh<Simplex<2> >,
bases<Lagrange<2,Vectorial>, Lagrange<1,Scalar>,
Lagrange<1,Scalar> > > P2P1P1;
auto U = P2P1P1.element();
// Views: changing a view changes U and vice versa
// view on element associated to P2
auto u = U.element<0>();
// extract view of first component
auto ux = u.comp(X);
// view on element associated to 1st P1
auto p = U.element<1>();
// view on element associated to 2nd P1
auto q = U.element<2>();``````

### Interpolation

Feel++ has a very powerful interpolation framework which allows to:

• transfer functions from one mesh to another

• transfer functions from one space type to another.

this is done seamlessly in parallel. The framework provides a set of C++ classes and C++ free-functions enabled short, concise and expressive handling of interpolation.

#### Using interpolation operator

Building interpolation operator I_h : P^1_{c,h} \rightarrow P^0_{td.h}
``````using MeshType = Mesh<Simplex<2>>;
auto mesh loadMesh( _mesh=new MeshType );
auto P1h = Pch<1>( mesh );
auto P0h = Pdh<0>( mesh );
auto Ih = I( _domain=P1h, _image=P0h );``````

#### De Rahm Diagram

The De Rahm diagram reads as follows: the range of each of the operators coincides with the null space of the next operator in the sequence below, and the last map is a surjection.

\begin{array}{ccccccc} H^1(\Omega)& \overset{\nabla}{\longrightarrow}& H^{\mathrm{curl}}(\Omega)& \overset{\nabla \times}{\longrightarrow}& H^{\mathrm{div}}(\Omega)& \overset{\nabla \cdot}{\longrightarrow}& L^2(\Omega) \end{array}

An important result is that the diagram transfers to the discrete level

\begin{array}{ccccccc} H^1(\Omega)& \overset{\nabla}{\longrightarrow}& H^{\mathrm{curl}}(\Omega)& \overset{\nabla \times}{\longrightarrow}& H^{\mathrm{div}}(\Omega)& \overset{\nabla \cdot}{\longrightarrow}& L^2(\Omega) \\ \left\downarrow\right.\pi_{c,h}& ~ & \left\downarrow\right.\pi_{\mathrm{curl},h}& ~ & \left\downarrow\right.\pi_{\mathrm{div},_h}& ~ & \left\downarrow\right.\pi_{d,h}& ~ \\ U_h& \overset{\nabla}{\longrightarrow}& V_h& \overset{\nabla \times}{\longrightarrow}& W_h& \overset{\nabla \cdot}{\longrightarrow}& Z_h\\ \end{array}

The diagram above is commutative which means that we have the following properties:

\begin{aligned} \nabla(\pi_{c,h} u) &= \pi_{\mathrm{curl},h}( \nabla u ),\\ \nabla\times(\pi_{\mathrm{curl},h} u) &= \pi_{\mathrm{div},h}( \nabla\times u ),\\ \nabla\cdot(\pi_{\mathrm{div},h} u) &= \pi_{d,h}( \nabla\cdot u ) \end{aligned}
 The diagram can be restricted to functions satisfying the homogeneous Dirichlet boundary conditions
\begin{array}{ccccccc} H^1_0(\Omega)& \overset{\nabla}{\longrightarrow}& H_0^{\mathrm{curl}}(\Omega)& \overset{\nabla \times}{\longrightarrow}& H_0^{\mathrm{div}}(\Omega)& \overset{\nabla \cdot}{\longrightarrow}& L^2_0(\Omega) \end{array}

Interpolation operators are provided as is or as shared pointers. The table below presents the alternatives.

 C++ object C++ Type C++ shared object C++ Type Mathematical operator `I(_domain=Xh,_image=Yh)` ```I_t, functionspace_type>``` `IPtr(…​)` ```I_ptr_t, functionspace_type>``` I: X_h \rightarrow Y_h `Grad(_domain=Xh,_image=Wh)` ```Grad_t, functionspace_type>``` `GradPtr(…​)` ```Grad_ptr_t, functionspace_type>``` \nabla: X_h \rightarrow W_h `Curl(_domain=Wh,_image=Vh)` ```Curl_t, functionspace_type>``` `CurlPtr(…​)` ```Curl_ptr_t, functionspace_type>``` \nabla \times : W_h \rightarrow V_h `Div(_domain=Vh,_image=Zh)` ```Div_t, functionspace_type>``` `DivPtr(…​)` ```Div_ptr_t, functionspace_type>``` \nabla \cdot: V_h \rightarrow Z_h
Building the discrete operators associated to the De Rahm diagram in Feel++
``````auto mesh = loadMesh( _mesh=new Mesh<Simplex<Dim>>());
auto Xh = Pch<1>(mesh);
auto Gh = Ned1h<0>(mesh);
auto Ch = Dh<0>(mesh);
auto P0h = Pdh<0>(mesh);
auto Icurl = Curl( _domainSpace = Gh, _imageSpace=Ch );
auto Idiv = Div( _domainSpace = Ch, _imageSpace=P0h );

auto u = Xh->element(<expr>);
auto w = Igrad(u); // w in Gh
auto x = Icurl(w); // z in Ch
auto y = Idiv(x); // y in P0h``````

#### Saving functions on disk

To save a function on disk to use it later, for example in another application, you can use the `save` function.

The saved file will be named after the name registered for the variable in the constructor (default : `u`).

``````auto Vh = Pch<1>( mesh );
auto u = Vh->element("v");
// do something with u
...
// save /path/to/save/v.h5
u.save( _path="/path/to/save", _type="hdf5" );``````

The `path` parameter creates a directory at this path to store all the degrees of liberty of this function.

The `type` parameter can be `binary`, `text` or `hdf5` . The first two will create one file per processor, whereas "hdf5" will creates only one file.

 To load a function, the mesh need to be exactly the same as the one used when saving it.
``````auto Vh = Pch<1>( mesh );
auto u = Vh->element("v");
The `path` and `type` parameters need to be the same as the one used to save the function.
In order to pass above this restriction, the two function `load` and `save` has been updated to use hdf5 format. With this format, extended parallel doftable or not, the function will work without any issues. More than that, we can load elements with extended parallel doftable and resave it without it, and vice versa. This last feature isn’t available with other formats than hdf5.