Feel++ Mathematical Keywords

One of Feel++ assets is it finite element embedded language. The language follows the C++ grammar, and provides keywords as well as operations between objects which are, mathematically, tensors of rank 0, 1 or 2.

In all following tables we use the notations: $f: \mathbb{R}^n \mapsto \mathbb{R}^{m\times p}$ with $n=1,2,3$, $m=1,2,3$, $p=1,2,3$.

We denote $\Omega^e$ current mesh element.

The genesis of the language can be found in cite:[prud2006domain].

1. Mesh

Function Description

`makeSharedMesh<T>`

Create a `boost::shared_ptr` mesh

`makeMesh<T>`

Create a `shared_ptr` mesh

`makeUniqueMesh<T>`

Create a `std::unique_ptr` mesh

`loadMesh`

`createSubmesh`

Create a submesh out of a mesh from a range of element of subentities

Create an empty Mesh of $d$-simplex
``````// Create a shared mesh ptr
// use default  constructor
auto mesh = makeMesh<Simplex<d>>();

// pass a worldComm
auto mesh = makeMesh<Simplex<d>>( Environment::worldComm() );

// Create a unique mesh ptr
// use default  constructor
auto mesh = makeUniqueMesh<Simplex<d>>();

// pass a worldComm
auto mesh = makeUniqueMesh<Simplex<d>>( Environment::worldComm() );``````
Load a mesh of $d$-simplex
``auto mesh=loadMesh( _mesh=new Mesh<Simplex<d>> );``

2. Mathematical Expressions

Following tables present tools to declare and manipulate expressions.

 Feel++ Keyword Description `Px()` Variable $x$ `Py()` Variable $y$ `Pz()` Variable $z$ `cst( c )` Constant function equal to $c$

You can of course use all current operators ( + - / * ) and the usual following functions:

 Feel++ Keyword Math Object Description `abs(expr)` $|f(\overrightarrow{x})|$ element wise absolute value of $f$ `cos(expr)` $\cos(f(\overrightarrow{x}))$ element wise cos value of $f$ `sin(expr)` $\sin(f(\overrightarrow{x}))$ element wise sin value of $f$ `tan(expr)` $\tan(f(\overrightarrow{x}))$ element wise tan value of $f$ `acos(expr)` $\mathrm{acos}(f(\overrightarrow{x}))$ element wise acos value of $f$ `asin(expr)` $\mathrm{asin}(f(\overrightarrow{x}))$ element wise asin value of $f$ `atan(expr)` $\mathrm{atan}(f(\overrightarrow{x}))$ element wise atan value of $f$ `cosh(expr)` $\cosh(f(\overrightarrow{x}))$ element wise cosh value of $f$ `sinh(expr)` $\sinh(f(\overrightarrow{x}))$ element wise sinh value of $f$ `tanh(expr)` $\tanh(f(\overrightarrow{x}))$ element wise tanh value of $f$ `exp(expr)` $\exp(f(\overrightarrow{x}))$ element wise exp value of $f$ `log(expr)` $\log(f(\overrightarrow{x}))$ element wise log value of $f$ `sqrt(expr)` $\sqrt{f(\overrightarrow{x})}$ element wise sqrt value of $f$ `ceil(expr)` $\lceil{f(\overrightarrow{x})}\rceil$ element wise ceil of $f$ `floor(expr)` $\lfloor{f(\overrightarrow{x})}\rfloor$ element wise floor of $f$ `sign(expr)` $\begin{cases} 1 & \text{if}\ f(\overrightarrow{x}) \geq 0\\-1 & \text{if}\ f(\overrightarrow{x}) < 0\end{cases}$ element wise sign value of $f$ `chi(expr)` $\chi(f(\overrightarrow{x}))=\begin{cases}0 & \text{if}\ f(\overrightarrow{x}) = 0\\1 & \text{if}\ f(\overrightarrow{x}) \neq 0\\\end{cases}$ element wise boolean test of $f$ `rand()` random number in $[0,1$] regenerated at each call `rand(lo,hi)` random number in $[lo,hi$] regenerated at each call

3. Geometry

3.1. Points

3.1.1. Current Point:

 Feel++ Keyword Math Object Description Dimension `P()` $\overrightarrow{P}$ $(P_x, P_y, P_z)^T$ $d \times 1$ `Px()` $P_x$ $x$ coordinate of $\overrightarrow{P}$ $1 \times 1$ `Py()` $P_y$ $y$ coordinate of $\overrightarrow{P}$ (value is 0 in 1D) $1 \times 1$ `Pz()` $P_z$ $z$ coordinate of $\overrightarrow{P}$ (value is 0 in 1D and 2D) $1 \times 1$

3.1.2. Element Barycenter Point:

Feel++ Keyword Math Object Description Dimension

`C()`

$\overrightarrow{C}$

$(C_x, C_y, C_z)^T$

$d \times 1$

`Cx()`

$C_x$

$x$ coordinate of $\overrightarrow{C}$

$1 \times 1$

`Cy()`

$C_y$

$y$ coordinate of $\overrightarrow{C}$ (value is 0 in 1D)

$1 \times 1$

`Cz()`

$C_z$

$z$ coordinate of $\overrightarrow{C}$ (value is 0 in 1D and 2D)

$1 \times 1$

3.1.3. Normal at Current Point:

 Feel++ Keyword Math Object Description Dimension `N()` $\overrightarrow{N}$ $(N_x, N_y, N_z)^T$ $d \times 1$ `Nx()` $N_x$ $x$ coordinate of $\overrightarrow{N}$ $1 \times 1$ `Ny()` $N_y$ $y$ coordinate of $\overrightarrow{N}$ (value is 0 in 1D) $1 \times 1$ `Nz()` $N_z$ $z$ coordinate of $\overrightarrow{N}$ (value is 0 in 1D and 2D) $1 \times 1$

3.1.4. Extra informations on elements and face of element:

Feel++ Keyword Math Object Description Dimension

`h()`

$h$

maximum edge size in element

1

`hMin()`

$h_{\mathrm{min}}$

minimum edge size in element

1

`hFace()`

$h_{\mathrm{f}}$

maximum edge size of current face

1

`meas()`

stem:[

\cdot

]

measure of the current element

1

`measFace()`

stem:[

\cdot

]

measure of the current face

1

`eid()`

id of the current element

1

`emarker()`

marker of the current element

1

`emarker2()`

marker2 of the current element

1

`fmarker()`

marker of the face in the current element when iterating over faces

1

`epid()`

pid of the current element

1

3.2. Geometric Transformations

3.2.1. Jacobian Matrix

You can access to the jacobian matrix, $J$, of the geometric transformation, using the keyword: `J()` There are some tools to manipulate this jacobian.

 Feel++ Keyword Math Object Description `detJ()` $\det(J)$ Determinant of jacobian matrix `invJT()` $(J^{-1})^T$ Transposed inverse of jacobian matrix

4. Vectors and Matrices

4.1. Building Vectors

Usual syntax to create vectors:

 Feel++ Keyword Math Object Description Dimension `vec(v_1,v_2,…​,v_n)` $\begin{pmatrix} v_1\\v_2\\ \vdots \\v_n \end{pmatrix}$ Column Vector with $n$ rows entries being expressions $n \times 1$

You can also use expressions and the unit base vectors:

 Feel++ Keyword Math Object Description `oneX()` $\begin{pmatrix} 1\\0\\0 \end{pmatrix}$ Unit vector $\overrightarrow{i}$ `oneY()` $\begin{pmatrix} 0\\1\\0 \end{pmatrix}$ Unit vector $\overrightarrow{j}$ `oneZ()` $\begin{pmatrix} 0\\0\\1 \end{pmatrix}$ Unit vector $\overrightarrow{k}$

4.2. Building Matrices

Table 1. Matrix and vectors creation
Feel++ Keyword Math Object Description Dimension

`mat<m,n>(m_11,m_12,…​,m_mn)`

$\begin{pmatrix} m_{11} & m_{12} & ...\\ m_{21} & m_{22} & ...\\ \vdots & & \end{pmatrix}$

$m\times n$ Matrix entries being expressions

$m \times n$

`ones<m,n>()`

$\begin{pmatrix} 1 & 1 & ...\\ 1 & 1 & ...\\ \vdots & & \end{pmatrix}$

$m\times n$ Matrix Filled with 1

$m \times n$

`zero<m,n>()`

$\begin{pmatrix} 0 & 0 & ...\\ 0 & 0 & ...\\ \vdots & & \end{pmatrix}$

$m\times n$ Matrix Filled with 0

$m \times n$

`constant<m,n>(c)`

$\begin{pmatrix} c & c & ...\\ c & c & ...\\ \vdots & & \end{pmatrix}$

$m\times n$ Matrix Filled with a constant c

$m \times n$

`eye<n>()`

$\begin{pmatrix} 1 & 0 & ...\\ 0 & 1 & ...\\ \vdots & & \end{pmatrix}$

Unit diagonal Matrix of size $n\times n$

$n \times n$

`Id<n>()`

$\begin{pmatrix} 1 & 0 & ...\\ 0 & 1 & ...\\ \vdots & & \end{pmatrix}$

Unit diagonal Matrix of size $n\times n$

$n \times n$

4.3. Manipulating Vectors and Matrices

Let $A$ and $B$ be two matrices (or two vectors of the operation is defined) of dimension $m \times n$. $A$ and $B$ may depend on space variables

Table 2. Matrix operations
Feel++ Keyword Math Object Description Dimension of the output

`inv(A)`

$A^{-1}$

Inverse of matrix $A$

$n \times n$

`det(A)`

$\det (A)$

Determinant of matrix $A$

$1 \times 1$

`sym(A)`

$\text{Sym}(A)$

Symmetric part of matrix $A$: $\frac{1}{2}(A+A^T)$

$n \times n$

`antisym(A)`

$\text{Asym}(A)$

Antisymmetric part of $A$: $\frac{1}{2}(A-A^T)$

$n \times n$

`trace(A)`

$\text{tr}(A)$

Trace of square matrix $A$

$1 \times 1$

`trans(B)`

$B^T$

Transpose of matrix $B$ Can be used on non-squared Matrix Can be used on Vectors

$n \times m$

`inner(A,B)`

$A \cdot B \\ A:B = \text{tr}(A*B^T)$

Scalar product of two vectors Generalized scalar product of two matrix

$1 \times 1$

`cross(A,B)`

$A\times B$

Cross product of two vectors

$n \times 1$

`eig(A)`

$\begin{pmatrix}\lambda_1 \\ \lambda_2 \\ \lambda_3\end{pmatrix}$

vector of real $n$ eigenvalues $\lambda_1,\lambda_2,\lambda_3)$ of $A$. $A$ must be symmetric.

$n \times 1$

`vonmises(A)`

$\begin{matrix} \text{1D:}& A_v = & A_{11}\\ \text{2D:}& A_v = & (A_{11}^2- A_{11}A_{22}+ A_{22}^2+3A_{12}^2)^{1/2} \\ \text{3D:}& A_v = & ((A_{11} - A_{22})^2 + (A_{22} - A_{33})^2 + \\ & & (A_{33} - A_{11})^2 + 6(A_{12}^2 + A_{23}^2 + A_{31}^2))^{1/2} \end{matrix}$

if $A$ is the Cauchy stress tensor, it computes the Von Mises yield criterion $A_v$ (scalar value). A material starts yielding when the von Mises stress $A_v$ reaches a value known as yield strength

$1 \times 1$

`sum(A)`

$rr$

`S = sum(A)` returns the sum of the elements of expression `A` along the first array dimension whose size does not equal 1.

* If `A` is a vectorial expression, then `sum(A)` returns the sum of the elements.

* If `A` is a matricial expression, then `sum(A)` returns a row vector containing the sum of each column

stem:[1 \times n]

 we have $eig(A)$

5. Operators

5.1. Operations

You can use the usual operations and logical operators.

5.3. Two Valued Operators

 Feel++ Keyword Math Object Description Rank Dimension `jump(f)` $[f]=f_0\overrightarrow{N_0}+f_1\overrightarrow{N_1}$ jump of test function 0 $n \times 1$ $m=1$ `jump(f)` $[\overrightarrow{f}]=\overrightarrow{f_0}\cdot\overrightarrow{N_0}+\overrightarrow{f_1}\cdot\overrightarrow{N_1}$ jump of test function 0 $1 \times 1$ $m=2$ `jumpt(f)` $[f]=f_0\overrightarrow{N_0}+f_1\overrightarrow{N_1}$ jump of trial function 0 $n \times 1$ $m=1$ `jumpt(f)` $[\overrightarrow{f}]=\overrightarrow{f_0}\cdot\overrightarrow{N_0}+\overrightarrow{f_1}\cdot\overrightarrow{N_1}$ jump of trial function 0 $1 \times 1$ $m=2$ `jumpv(f)` $[f]=f_0\overrightarrow{N_0}+f_1\overrightarrow{N_1}$ jump of function evaluation 0 $n \times 1$ $m=1$ `jumpv(f)` $[\overrightarrow{f}]=\overrightarrow{f_0}\cdot\overrightarrow{N_0}+\overrightarrow{f_1}\cdot\overrightarrow{N_1}$ jump of function evaluation 0 $1 \times 1$ $m=2$ `average(f)` ${f}=\frac{1}{2}(f_0+f_1)$ average of test function rank$( f(\overrightarrow{x}))$ $n \times n$ $m=n$ `averaget(f)` ${f}=\frac{1}{2}(f_0+f_1)$ average of trial function rank$( f(\overrightarrow{x}))$ $n \times n$ $m=n$ `averagev(f)` ${f}=\frac{1}{2}(f_0+f_1)$ average of function evaluation rank$( f(\overrightarrow{x}))$ $n \times n$ $m=n$ `leftface(f)` $f_0$ left test function rank$( f(\overrightarrow{x}))$ $n \times n$ $m=n$ `leftfacet(f)` $f_0$ left trial function rank$( f(\overrightarrow{x}))$ $n \times n$ $m=n$ `leftfacev(f)` $f_0$ left function evaluation rank$( f(\overrightarrow{x}))$ $n \times n$ $m=n$ `rightface(f)` $f_1$ right test function rank$( f(\overrightarrow{x}))$ $n \times n$ $m=n$ `rightfacet(f)` $f_1$ right trial function rank$( f(\overrightarrow{x}))$ $n \times n$ $m=n$ `rightfacev(f)` $f_1$ right function evaluation rank$( f(\overrightarrow{x}))$ $n \times n$ $m=n$ `maxface(f)` $\max(f_0,f_1)$ maximum of right and left test function rank$( f(\overrightarrow{x}))$ $n \times p$ `maxfacet(f)` $\max(f_0,f_1)$ maximum of right and lef trial function rank$( f(\overrightarrow{x}))$ $n \times p$ `maxfacev(f)` $\max(f_0,f_1)$ maximum of right and left function evaluation rank$( f(\overrightarrow{x}))$ $n \times p$ `minface(f)` $\min(f_0,f_1)$ minimum of right and left test function rank$( f(\overrightarrow{x}))$ $n \times p$ `minfacet(f)` $\min(f_0,f_1)$ minimum of right and left trial function rank$( f(\overrightarrow{x}))$ $n \times p$ `minfacev(f)` $\min(f_0,f_1)$ minimum of right and left function evaluation rank$( f(\overrightarrow{x}))$ $n \times p$

6. Fitting

Keywords `fit` and `fitDiff` provide the interpolated data and the derivative of the interpolated data respectively. This is useful when material laws, properties or some terms in variational formulation are given as a data file.

``````auto Xh = Pch<2>(mesh);
auto K = Xh->element();
K.on(_range=elements(mesh), _expr=fit(idv(T)[,dataFile(string),[type(string)]]));
Kd.on(_range=elements(mesh), _expr=fitDiff(idv(T)[,dataFile(string),[type(string)]]));``````

6.1. Options

 `--fit.datafile` the data file to interpolate (two columns) `--fit.type` P0, P1, Spline, Akima `--fit.P0` left = 0, right = 1, center = 2 `--fit.P1_right` Kind of extention : zero = 0, constant = 1, extrapol = 2 `--fit.P1_left` Kind of extention : zero = 0, constant = 1, extrapol = 2 `--fit.Spline_right` natural = 0, clamped = 1 `--fit.Spline_left` natural = 0, clamped = 1
``````    auto mesh = loadMesh(_mesh=new Mesh<Simplex<2>>);
auto Xh = Pch<1>(mesh);
auto T = Xh->element(); // the temperature (say)
auto K = Xh->element(); // K(T) - the dependant of the temperature conductivity
auto Kd= Xh->element(); // K'(T)
auto T_K = Xh->element(); // T_ = true
auto T_Kd= Xh->element(); //
auto D_K = Xh->element(); // D_ = difference
auto D_Kd= Xh->element(); //
T.on(_range=elements(mesh), _expr=Px());
T_K.on(_range=elements(mesh),_expr=(5*Px()+sin(Px())));
T_Kd.on(_range=elements(mesh),_expr=(5+cos(Px())));
//double f(double t = 0.) { return 5. * t + sin(t); }
auto f = [](double x = 0.) { return 5. * x + sin(x); };
#if 1
auto e = exporter(_mesh = mesh );
#endif
std::string datafilename = (fs::current_path()/"data.csv").string();
if ( Environment::worldComm().isMasterRank() )
{
// Generates the datafile
// we assume an unitsquare as mesh
std::ofstream datafile( datafilename );
datafile << "x, y\n";
for(double t = -1; t < 2; t+=0.32)
datafile << t << " , " << f(t) << "\n";
datafile.close();
}
Environment::worldComm().barrier();

std::vector<std::string> interpTypeRange = { "P0" , "P1", "Spline", "Akima" };
for(int k = 0; k < interpTypeRange.size(); ++k )
{
std::string const& interpType = interpTypeRange[k];
BOOST_TEST_MESSAGE( boost::format("test %1%")% interpType );
// evaluate K(T) with the interpolation from the datafile
K.on(_range=elements(mesh), _expr=fit( idv(T), datafilename, "x", "y", interpType ) );
Kd.on(_range=elements(mesh), _expr=fitDiff( idv(T), datafilename, "x", "y", interpType ) );

D_K.on(_range=elements(mesh),_expr=vf::abs(idv(K)-idv(T_K)));
D_Kd.on(_range=elements(mesh),_expr=vf::abs(idv(Kd)-idv(T_Kd)));

auto max_K = D_K.max();
auto max_Kd= D_Kd.max();
#if 1
e->save();
#endif
/// Note : the error has nothing to do with the mesh size but the step on the datafile
switch( InterpolationTypeMap[interpType] )
{
case InterpolationType::P0: //P0 interpolation
{
BOOST_CHECK_SMALL(max_K, 0.95);
BOOST_CHECK_SMALL(max_Kd, 6.0001); // the derivative is null
break;
}
case InterpolationType::P1: // P1 interpolation
{
BOOST_CHECK_SMALL(max_K, 0.01);
BOOST_CHECK_SMALL(max_Kd, 0.15);
break;
}
case InterpolationType::Spline: // CSpline interpolation
{
BOOST_CHECK_SMALL(max_K, 0.01);
BOOST_CHECK_SMALL(max_Kd, 0.15);
break;
}
case InterpolationType::Akima: // Akima interpolation
{
BOOST_CHECK_SMALL(max_K, 0.016);
BOOST_CHECK_SMALL(max_Kd, 0.03);
break;
}
}

BOOST_TEST_MESSAGE( boost::format("test %1% done")% interpType );
}``````