Feel++ Book Contributors https://github.com/feelpp/book.feelpp.org/graphs/contributors :sources: ../../../../codes/ :sourcedir: ../../../../codes/

## 1. Programming Syntax Quick Reference

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 (prud2006domain).

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

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

### 1.3. Geometry

#### 1.3.1. Points

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

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

#### 1.3.2. Geometric Transformations

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

### 1.4. Vectors and Matrix

#### 1.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}$

#### 1.4.2. Building Matrix

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$

#### 1.4.3. Manipulating Vectors and Matrix

Let $A$ and $B$ be two matrix (or two vectors) of same dimension $m \times n$.

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

`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 matrix $A$ Generalized on non-squared Matrix Generalized on Vectors

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

### 1.5. Operators

#### 1.5.1. Operations

You can use the usual operations and logical operators.

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

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

#### 1.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.txt").string();
if ( Environment::worldComm().isMasterRank() )
{
// Generates the datafile
// we assume an unitsquare as mesh
std::ofstream datafile( datafilename );
for(double t = -1; t < 2; t+=0.32)
datafile << t << "\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, interpType ) );
Kd.on(_range=elements(mesh), _expr=fitDiff( idv(T), datafilename, 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 );
}``````

## Appendix A: Bibliography

1. prud2006domain