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 [1].

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

load a 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
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<n>(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.

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<g\)

element wise less

\( f⇐g\)

element wise less or equal

>

\( 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

1.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\} \) 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 \)

dn(f)

\( \nabla f \cdot \overrightarrow \)

normal derivative of test function

0

\( 1 \times 1 \) \( m=p=1 \)

dn(f)

\( \nabla f \ \overrightarrow \)

normal derivative of test function

1

\( m \times 1 \) \(p=1 \)

dnt(f)

\( \nabla f \cdot \overrightarrow \)

normal derivative of trial function

0

\( 1 \times1 \) \(m=p=1\)

dnt(f)

\( \nabla f \ \overrightarrow \)

normal derivative of trial function

1

\( m \times 1 \) \(p=1\)

dnv(f)

\( \nabla f \cdot \ \overrightarrow \)

evaluation of normal derivative

0

\( 1 \times 1 \) \(m=p=1\)

dnv(f)

\( \nabla f \ \overrightarrow \)

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 \)

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 \)

Unresolved directive in ref/Keywords/README.adoc - include::image.adoc[]

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->step(k)->add("T",T);
        e->step(k)->add("K",K);
        e->step(k)->add("Kd",Kd);
        e->step(k)->add("T_K",T_K);
        e->step(k)->add("T_Kd",T_Kd);
        e->step(k)->add("D_K",D_K);
        e->step(k)->add("D_Kd",D_Kd);
        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. C. Prud’homme, “A domain specific embedded language in C++ for automatic differentiation, projection, integration and variational formulations,” Scientific Programming, vol. 14, no. 2, pp. 81–110, 2006.