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

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

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

element wise less

⇐

\( f\leq g\)

element wise less or equal

>

\( 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->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 );
    }