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

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

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>=g$ element wise greater or equal == $f==g$ element wise equal != $f!=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 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 );
}