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<N>(mesh)

Pch_type<MeshType,N>

\$P^N_{c,h}\$

Pchv<N>(mesh)

Pchv_type<MeshType,N>

\$[P^N_{c,h}\$^d]

Pdh<N>(mesh)

Pdh_type<MeshType,N>

\$P^N_{d,h}\$

Pdhv<N>(mesh)

Pdhv_type<MeshType,N>

\$[P^N_{d,h}\$^d]

THch<N>(mesh)

THch_type<MeshType,N>

\$[P^{N+1}_{c,h}\$^d \times P^N_{c,h}]

Dh<N>(mesh)

Dh_type<MeshType,N>

\$\mathbb{R}\mathbb{T}_h\$

Ned1h<N>(mesh)

Ned1h_type<MeshType,N>

\$\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.

Table 1. Table of Interpolation operators

C++ object

C++ Type

C++ shared object

C++ Type

Mathematical operator

I(_domain=Xh,_image=Yh)

I_t<functionspace_type<decltype(Xh)>, functionspace_type<decltype(Xh)>>

IPtr(…​)

I_ptr_t<functionspace_type<decltype(Xh)>, functionspace_type<decltype(Xh)>>

I: X_h \rightarrow Y_h

Grad(_domain=Xh,_image=Wh)

Grad_t<functionspace_type<decltype(Xh)>, functionspace_type<decltype(Wh)>>

GradPtr(…​)

Grad_ptr_t<functionspace_type<decltype(Xh)>, functionspace_type<decltype(Wh)>>

\nabla: X_h \rightarrow W_h

Curl(_domain=Wh,_image=Vh)

Curl_t<functionspace_type<decltype(Wh)>, functionspace_type<decltype(Vh)>>

CurlPtr(…​)

Curl_ptr_t<functionspace_type<decltype(Wh)>, functionspace_type<decltype(Vh)>>

\nabla \times : W_h \rightarrow V_h

Div(_domain=Vh,_image=Zh)

Div_t<functionspace_type<decltype(Vh)>, functionspace_type<decltype(Zh)>>

DivPtr(…​)

Div_ptr_t<functionspace_type<decltype(Vh)>, functionspace_type<decltype(Zh)>>

\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 Igrad = Grad( _domainSpace = Xh, _imageSpace=Gh );
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 and loading functions on disk

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.

Loading functions from disk

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");
// load /path/to/load/v.h5
u.load( _path="/path/to/load/", _type="hdf5" );

The path and type parameters need to be the same as the one used to save the function.

Extended parallel doftable

In some cases, when we use parallel data, informations from other interfaces of partitions are need. To manage this, we can add ghost degree of freedom on ghost elements at these locations. However, we have to know if data have extended parallel doftable to load and use it.

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.