# 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.

 Feel++ Keyword Math Object Description `+` $f+g$ tensor sum `-` $f-g$ tensor substraction `*` $f*g$ tensor product `/` $f/g$ tensor tensor division ($g$ scalar field) `<` $f` \( f>g$ element wise greater `>=` $f \geq g$ element wise greater or equal `==` $f==g$ element wise equal `!=` $f \neq g$ element wise not equal `-` $-g$ element wise unary minus `&&` $f$ and $g$ element wise logical and `||` $f$ or $g$ element wise logical or `!` $!g$ element wise logical not

### 5.2. Differential Operators

Feel++ finit element language use test and trial functions. Keywords are different according to the kind of the manipulated function.
Usual operators are for test functions.
t-operators for trial functions.
v-operators to get an evaluation.

Suppose that $f \in X_h$ reads

$f=\sum_{i=0}^{\mathcal{N}} f_i \phi_i$

where $X_h = \mathrm{span}\{ \phi_i, i=1,\ldots,\mathcal{N}\}$ is a finite element space.

 Feel++ Keyword Math Object Description Rank Dimension `id(f)` $\{\phi_i\}$ test function rank$(f(\overrightarrow{x}))$ $m \times p$ `idt(f)` $\{\phi_i\}$ trial function rank$(f(\overrightarrow{x}))$ $m \times p$ `idv(f)` $f$ evaluation function rank$(f(\overrightarrow{x}))$ $m \times p$ `grad(f)` $\nabla f$ gradient of test function rank$(f(\overrightarrow{x}))+1$ $m \times n$ $p=1$ `gradt(f)` $\nabla f$ grdient of trial function rank$(f(\overrightarrow{x}))+1$ $m \times n$ $p=1$ `gradv(f)` $\nabla f$ evaluation function gradient rank$(f(\overrightarrow{x}))+1$ $m \times n$ $p=1$ `div(f)` $\nabla\cdot f$ divergence of test function rank$(f(\overrightarrow{x}))-1$ $1 \times 1$ `divt(f)` $\nabla\cdot f$ divergence of trial function rank$(f(\overrightarrow{x}))-1$ $1 \times 1$ `divv(f)` $\nabla\cdot f$ evaluation function divergence rank$(f(\overrightarrow{x}))-1$ $1 \times 1$ `curl(f)` $\nabla\times f$ curl of test function 1 $n \times 1$ $m=n$ `curlt(f)` $\nabla\times f$ curl of trial function 1 $n \times 1$ $m=n$ `curlv(f)` $\nabla\times f$ evaluation function curl 1 $n \times 1$ $m=n$ `laplacian(f)` $\Delta f$ laplacian of test function 0 $1 \times 1$ $m=p=1$ `laplaciant(f)` $\Delta f$ laplacian of trial function 0 $1 \times 1$ $m=p=1$ `laplacianv(f)` $\Delta f$ laplacian of the function $f$ 0 $1 \times 1$ $m=p=1$ `hess(f)` $\nabla^2 f$ hessian of test function 2 $n \times n$ $m=p=1$ `trace(f)` $\mathrm{trace}( f )$ trace of test matrix field 0 $1 \times 1$ $m=p=d$ `tracet(f)` $\mathrm{trace}( f )$ trace of trial matrix field 0 $1 \times 1$ $m=p=d$ `tracev(f)` $\mathrm{trace}( f )$ trace of matrix field $f$ 0 $1 \times 1$ $m=p=d$ `normal(f)` $f \cdot \overrightarrow{N}$ normal component of test function rank$(f(\overrightarrow{x}))-1$ `normalt(f)` $f \cdot \overrightarrow{N}$ normal component of trial function rank$(f(\overrightarrow{x}))-1$ `normalv(f)` $f \cdot \overrightarrow{N}$ normal component of function $f$ rank$(f(\overrightarrow{x}))-1$ `dn(f)` $\nabla f \cdot \overrightarrow{N}$ normal derivative of test function 0 $1 \times 1$ $m=p=1$ `dn(f)` $\nabla f \ \overrightarrow{N}$ normal derivative of test function 1 $m \times 1$ $p=1$ `dnt(f)` $\nabla f \cdot \overrightarrow{N}$ normal derivative of trial function 0 $1 \times1$ $m=p=1$ `dnt(f)` $\nabla f \ \overrightarrow{N}$ normal derivative of trial function 1 $m \times 1$ $p=1$ `dnv(f)` $\nabla f \cdot \ \overrightarrow{N}$ evaluation of normal derivative 0 $1 \times 1$ $m=p=1$ `dnv(f)` $\nabla f \ \overrightarrow{N}$ evaluation of normal derivative 1 $m \times 1$ $p=1$ `dx(f)` $\nabla f \cdot \overrightarrow{i}$ derivative of test function in $x$ 0 $1 \times 1$ $m=p=1$ `dy(f)` $\nabla f \cdot \overrightarrow{j}$ derivative of test function in $y$ 0 $1 \times 1$ $m=p=1$ `dz(f)` $\nabla f \cdot \overrightarrow{k}$ derivative of test function in $z$ 0 $1 \times 1$ $m=p=1$

### 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. Multiscale Image

in order to deal with large scale image, several tools have been developed to facilitate access to images.

Denote $f \in \mathbb{R}^N$ where $N$ is the dimension of the image $N=N_x N_y$ with $N_x, N_y$ are the number of pixels in direction $X$ and $Y$ respectively. $f$ contains the pixel value.

 define image format pending NOTE: define transformation pending
 Feel++ Keyword Math Object Description `msi(f)` $T(f)$ the transformation from the coarse to fine level

## 7. 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)]]));``````

### 7.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 );
}``````