A domain speciﬁc embedded language in C ++ for automatic differentiation, projection, integration and variational formulations

. In this article, we present a domain speciﬁc embedded language in C ++ that can be used in various contexts such as numerical projection onto a functional space, numerical integration, variational formulations and automatic differentiation. Albeit these tools operate in different ways, the language overcomes this difﬁculty by decoupling expression constructions from evaluation. The language is implemented using expression templates and meta-programming techniques and uses various Boost libraries. The language is exercised on a number of non-trivial examples and a benchmark presents the performance behavior on a few test problems.


Introduction
Numerical analysis tools such as differentiation, integration, polynomial approximations or finite element approximations are standard and mainstream tools in scientific computing.Many excellent libraries or programs provide a high level programming interface to these methods : (i) programs that define a specific language such as the Freefem software family [12,21], the Fenics project [17,18], Getdp [11] or Getfem++ [23], or (ii) libraries or frameworks that supply some kind of domain specific language embedded in the programming language -hereafter called DSELsuch as LifeV (C++) [1,20], Sundance (C++) [19], Analysa (Scheme, which is suited for embedding sub-languages like other Lisp based languages) [5].
These high level interfaces or languages are desirable for several reasons: teaching purposes, solving complex problems with multiple physics and scales or rapid prototyping of new methods, schemes or algorithms.The goal is always to hide (ideally all) technical details behind software layers and provide only the relevant components required by the user or programmer.
The DSEL approach has advantages over generating a specific language like in case (i) : compiler construction complexities can be ignored, other libraries can concurrently be used which is often not the case of specific languages which would have to also develop their own libraries and DSELs inherit the capabilities of the language in which they are written.However, DSELs are often defined for one particular task inside a specific domain [24] and implementation or parts of implementation are not shared between different DSELs.
This article proposes a DSEL for automatic differentiation, projection, integration or variational formulations.The language implementation uses expression templates [24] and other meta-programming techniques [2].Related works are [20] and [19], but they differ with the proposed DSEL in many aspects: the former was designed only for variational formulations and requires to write the expression object by hand which can become complicated and error prone, while the latter implements the DSEL in an object oriented way without relying on meta-programming or expression templates.
Other objectives of our DSEL implementation are that it should (i) be efficient enough to integrate high performance/parallel software, (ii) be generic enough to accommodate different numerical types -for example arbitrary precision, see [22] but we won't discuss these aspects here.-A performance benchmark is available in Section 4.1.
To illustrate further what the DSEL achieves, here is a comparison between a mathematical formulation of a bilinear form Eq. and its programming counterpart, see listing 1.We clearly identify in listing 1 the variational formulation stated in Eq. (1).We shall describe the various steps to achieve this level of expression with as little overhead as possible.In Section 2, we present some concepts concerning mainly integration and variational formulations, then in Section 3 we present the main points about the DSEL.Finally in Section 4, we present some non-trivial examples to exercise the language.
This article contains many listings written in C++ however most of them are not correct C++ in order to simplify the exposition.In particular, C++ keywords like typename or inline are often not present.Also many numerical ingredients such as polynomial approximations, numerical integration methods used in this article are not described or only very roughly, another publication will cover the mathematical kernel used by the DSEL in more details [22].

Preliminaries on variational forms
In what follows, we consider a domain Ω ⊂ R d , d = 1, 2, 3 and its associated mesh T -out of d-simplices and product of simplices.

Mesh
We present first some tools that will be used later, namely how to extract parts of a mesh and the geometric mapping that maps a convex of reference -where polynomial sets and quadratures are constructed -to any convex of the mesh.

Mesh parts extraction
While applying integration and projection methods, it is common to be able to extract parts of the mesh.Hereafter we consider only elements of the mesh and elements faces.We wish to extract easily subsets of convexes out the total set consituting T .
To do this our mesh data structure which is by all means fairly standard uses the Boost.Multi index library1 to store the elements, elements faces, edges and points.This way the mesh entities are indexed either by their ids, their markers -material properties, boundary ids . .., -their location -whether the entity is internal or lies on the boundary of the domain.-Other indices could be certainly defined, however those three allow already a wide range of applications. 2hanks to Boost.Multi index, it is trivial to retrieve pairs of iterators over the entities -elements, faces, edges, points -containers depending on the usage context.The pairs of iterators are then turned into a range, see Boost.Range, 3 to be manipulated by the integration and projection tools that will be presented later.
A number of free functions are available that hide all details about the mesh class to concentrate only on the relevant parts.

Geometric mapping
Functional spaces and quadrature methods, for example, are derived from polynomial sets or families that have to be constructed over the convexes of T .Instead of doing this, it is common to construct these polynomials over a reference convex T -segment, triangle, quadrangle, tetrahedron, hexahedron, prism or pyramid -and provide a geometric mapping or transformation from the reference convex T to any convex T ∈ T ⊂ R P , P d.We need to be able also to transform subentities, such as faces or points, of the reference element to the corresponding entities, faces and points respectively, in the real element.
From now on, we denote with a (ˆ) the quantities defined on the reference element.We define τ : R P → R N that maps T to T .We shall denote K τ its gradient, B τ its pseudo-inverse and J τ its jacobian.The geometric mapping is described by (i) a n g components polynomial vector {φ g (x)} g=1...ng and (ii) the geometric points {p g } g=1...ng of T such that Fig. 1.Geometric mapping τ from the reference tetrahedron T to a real tetrahedron T in 3D.
We denote by G = (p 1 , . . ., p ng ) the N × n g matrix of geometric nodes.Equation (2) and quantities mentioned above are computed as follows, for any x ∈ T denotes the transpose of K τ (x).Equipped with the geometric mapping concept, we compute an integral on T as an integral on T : if f is a function defined on T , and using a quadrature formula: where {x q , ŵq } q=1...Q are quadrature nodes and quadrature weights defined in the reference element.
In our framework, the geometric mapping is not used directly by the developer but rather what we call the geometric mapping context which is a subclass of the geometric mapping class.The geometric mapping context is linked to an element T of the mesh such that, given a set of points {x ∈ T }, it provides information for each point in the set {x; x ∈ T and x = τ (x)} such as the jacobian value J τ (x), the gradient K τ (x) of the mapping, the pseudo-inverse B τ (x) of the gradient; or if the point x is on a face of x, then x is on a face of T and the context provides the normal to the face at this point.A shortened interface of the Context class is presented in the listing 2.
Another subclass of the geometric mapping class is Inverse which, as its name state, does the inverse of the transformation: given a point x in T ⊂ R N compute its location in the reference element x in T ⊂ R P .Inverse is particularly useful for interpolation purposes.
We define a bilinear form a : X × Y → R and a linear form : X → R, where X and Y are suitable function spaces defined on Ω.The finite element method discretizes X and Y using polynomials spaces defined on T .We denote by N X and N Y the dimension of the discrete spaces X and Y 5 and by {ψ i } i=1...NX and {ϕ i } i=1...NY a Listing 2: Geometric mapping context class basis for X and Y respectively.For any v ∈ X, we have v = i=1...N v i ψ i , and similarly for the functions of Y .We can then write the entries A ij , i = 1 . . .N X , j = 1 . . .N Y of the matrix A associated with a and the entries L i , i = 1 . . .N X of the vector L associated with as follows: To construct A and L, we follow a standard assembly process that iterates over the elements T of T since a(v, u) can be written as T ∈T a T (v, u) and similarly with and L. We then introduce (i) the restriction of the basis functions to T , {ψ T i } i=1...NX and {ϕ T i } i=1...NY where i is a local numbering over T , and (ii) the local to global mappings ι X (•, •) and ι Y (•, •) between the local numbering of the degrees of freedom and the global one.For example ι X (T, i) is the global degree of freedom to which the i-th local degree of freedom of T contributes to.The assembly process is described in the Algorithm 1.

Algorithm 1 Standard assembly procedure for A and L
In standard finite element software, the assembly is often split into two steps : (i) the local matrix A T = (a T (ψ T i , ϕ T j )) i=1...NX ,j=1...NY and vector L T = ( T (ψ T i )) i=1...NX are first constructed and (ii) the local to global mapping is used to add the contribution of the element T to A and L. This splitting is often used to optimize the local to global mappings [7] or optimize the local matrix and vector computation [17].We will also follow this strategy in the remaining sections.

Construction of a
We focus now on the construction of the elementary contribution a T (ψ T i , ϕ T j ) and T (ψ T i ) which is the case of our methodology.

Basis functions
We turn to the treatment of the basis functions in our framework, and in particular we describe the computation of f (τ (x)) for any x ∈ T as in Eq. ( 5).We define the finite element basis functions on the reference element.If f belongs to X, then we have for a given T ∈ T and its associated geometric mapping τ : where Similar computations, albeit more involved, can be derived for the second order derivatives.
The basis function concept we developed is similar to the geometric mapping.In our framework, the degrees of freedom are associated with the elements of the mesh.More precisely they are ordered with respect to the geometric subentities of the elements -vertices, edges, faces and volumes -for global continuous functions to ensure a continuous expansion whereas in the case of global discontinuous functions it does not matter how the degress of freedom are ordered or organized within the element.This allows for flexible construction of polynomial sets such as Lagrange, Raviart-Thomas or modal basis with global continuous expansion or not.The article [22] presents these aspects in details.Essentially the Basis base class, see listing 3, provides an interface for obtaining the value of the basis function and its derivatives at any given point in the reference element.Similar to the geometric mapping, we also define a Context subclass that provides information on the basis functions at a given set of points {x; x ∈ T }.
Equipped with these tools and concepts and if we consider a function f ∈ X, we have Listing 3: Basis functions interface Listing 4: Approximation Space Interface

Approximation space
We define the notion of approximation space in C++ that maps closely the mathematical counterpart.An approximation space is a template class parametrized by a mesh class and the basis functions type -for example the standard Lagrange finite elements.-An approximation space wraps the mesh, the table of degrees of freedom (DoF), the basis function type and provides access to all them.Note that the geometric mapping is provided by the mesh class.
A Space defines its own element type as a subclass: it ensures coherence and consistency when manipulating finite element functions.An Element derives from your preferred numerical vector type: we use uBLAS6 for our linear algebra data structures and algorithms.The interface is roughly described in the listing 5.
An extension of the Space concept, is the MixedSpace which is a product of two spaces.This can actually be extended to a product of several spaces of different types -implemented using the MPL [2].-This concept is useful for mixed formulations.MixedSpace defines also its own element type with some extra member to retrieve the underlying space elements.At the moment, MixedSpace is a function concept for a product of two functional spaces.Extending this concept to a product of N function spaces would be useful.

Linear and bilinear forms
One last concept needed to have the language expressive is the notion of forms.They follow closely their mathematical counterparts: they are template classes with arguments being the space or product of spaces they take their input from and the representation we can make out of these forms.In what follows, we consider only the case where the linear and bilinear forms are represented by vectors and matrices respectively.In a future work, we will eventually propose the possibility to have vector-free and matrix-free representations: that would require to store the definition of the forms.
Listing 7 displays the basic interface and usage of the form classes.Note that the linear and bilinear form classes are the glue between their representation and the mathematical expression given by Expr, it will -fill the matrix with non-zero entries depending on the approximation space(s) and the mathematical expression; -allow a per-component/per-space construction(blockwise); -check that the numerical types of the expression and the representation are consistent; -when operator=(Expr const&) is called, the expression is evaluated and will fill the representation entries.
The concepts of MixedLinearForm and MixedBilinearForm that would correspond to mixed linear and bilinear forms respectively -taking their values in the product of two functional spaces -exist also and follow the same ideas.
With the high level concepts described we can now focus on the language.

Language
The expression template technique won't be described as it is nowadays a mainstream technique [2][3][4]20,24].The construction of the expression template objects in the coming sections is standard.

Expression evaluation at a set of points in a convex
Let C be a convex in R d , d 1,2,3 -a n-simplex n d like lines, triangles or tetrahedrons or products of simplices like quadrangles, hexahedrons or prisms, -and Ĉ be a convex of reference in R d , d 1, 2, 3 associated to C where we define quadrature points for integration or points to construct polynomials for finite elements and other approximation methods, see [10,15].
We wish to evaluate In our code f is represented by an expression template -and not a standard C++ function or a functor, -see [24].For example, consider f : x ∈ C → cos(πx) sin(πy), we write it in C++ as cos(π*Px())*sin(π*Py()).The expression graph is shown on Fig. 2.Here Px() and Py() are free functions that construct objects that are evaluated as the x and y coordinates of the points x; x ∈ C.
Constructing the C++ object that represents the expression is done with standard expression template approach.However evaluating the expression is problematic as some ingredients are not known yet to the expression object such as the geometric mapping.So using a standard expression template approach certainly allow high level expressivity but cannot be applied to evaluate the expression.
To remedy this issue, we propose a very simple but very powerful solution which delegates the evaluation of the expression to another object than the expression object itself.In our case, the evaluation is delegated to a sublass of each object of the expression.
The Expression class, which is the glue between the various object types forming the expression template, is roughly sketched in listing 8.
Expression<Expr>::Eval is a template class parametrized by the geometric mapping context associated with each geometric element of the mesh.The constructor takes the expression Expr and the geometric mapping context as arguments to pass geometric data -coordinates of the current point, normals, measure of the elementdown to all objects of the expression so that they can use it as needed.As already mentioned, Px() constructs a C++ Py() is implemented in a similar way.Regarding the mathematical functors cos and sin, they also follow the same idea as shown in listing 10.

Nodal projection
We described the mechanism to evaluate an expression at a set of points in a convex, we now turn to nodal projection of a function f onto an approximation space X -for example X = {u ∈ P k (T )} where T is a triangulation of Ω and P k is the set spanned by the Lagrange polynomials of degree k. -We denote π X f the nodal projection of f onto X.
The nodal projection is an extension of the previous section at a set of convexes and ŜP being the set of coordinates of the degrees of freedom (DoF) associated with X.The nodal projection is described by Algorithm 2.

x) end for end for
We define a free function project(<space>,[elements,]<expression>) that takes two or three arguments : the approximation space onto which we project the function, the expression representing the function we want to project and optionally a range of elements that restricts the projection to this set of elements, see Section 2.1.1.project() constructs an template class parametrized by the arguments types passed to project(), see listing 11.
Listing 12 shows an example of nodal projection.Other types of projection like L 2 or H 1 projections require other ingredients presented in the coming sections.

Numerical integration
We now turn to numerical integration of Ω f (x)dx where f is the function to be integration over Ω. Numerical integration requires the evaluation of the function f at quadrature points in the convexes of the mesh associated to Ω.In our code, we used the quadrature constructions presented in [15] for n-simplices and simplices products.
The integration process is described by Algorithm 3. We introduce a new keyword to reflect the integration action, see listing 13, which is a free function instantiating an integrator parametrized by (i) the set of geometric elements, see Section 2.1.1,where the integration is done, (ii) the expression to integrate and (iii) the integration method, see listing 14.

Variational formulations
The framework, presented in the last sections, can be extended to handle variational formulation with only minor changes to the evaluation class.We consider for now a convex C ⊂ R d , d = 1, 2, 3 and its associated reference convex Ĉ, an approximation space X -for example P k (C) -a bilinear form a X × X → R defined by Listing 16 shows the C++ counterpart of Eq. ( 13).The t in idt(.)allows to distinguish trial and test functions: for example, id(.) identifies the test function values whereas idt(.)identifies the trial function values.
Given u, v ∈ X, we wish to compute the value a(u, v) which can be approximated as follows ŵq ψi (x q ) ψj (x q )J(x q ) where ( ŵq , xq ) q=1...Q are the quadrature nodes and weights defined in Ĉ, J(x q ) is the jacobian of the geometric transformation between Ĉ and C at the point xq and ( ψi ) i=1...NX is the basis of X. Recall Section 2.2.1, we have the basis context subclass that allows to store values and derivatives of basis functions at a set of points.In the case of Eq. , the basis context subclass stores and provides an interface to ( ψi (x q )) q=1...Q,i=1...NX .
Again the basis functions context is not known to the expression object.In order to accommodate the language with these concepts, it suffices to add new template parameters to the evaluation subclass of each classes allowed in an expression.However these parameters have default values that allows to handle the case of the previous section as well as linear forms and bilinear forms, see listing 17.In particular test basis functions context type are defaulted to boost::none t7 and trial basis functions context type to the test ones.If boost::none t is used, then no language keyword may be used in the expression that will need basis functions operators such as id(.) or idt(.).Let's examine for example the implementation of the operator dx(<element:u>) which provides the first component of the first derivative of the basis function associated with u, see listing 18.
The types of the test and trial basis functions, Basis test t and Basis trial t are tested with respect to the basis function type basis type of the element passed to the operator.If the types are not the same then at compile time the evaluation of this operator returns 0. At first glance for standard scalar equations -heat equation say -it allows just to ensure that the evaluation makes sense, however when dealing with mixed formulation -e.g.Stokes equations -this feature becomes very important if not crucial for correctness and performance-wise.Indeed it will disable automatically the terms associated with the basis functions which are not evaluated during local assembly of the mixed formulation: for example when filling the block corresponding to the velocity space, all the others termsvelocity-pressure and pressure terms -are evaluated to 0.
An immediate remark is that we may have lots of computation for nothing, i.e. computing 0. However a simple check with g++ -O2 --save-temps shows that g++ optimizes out expressions containing 0 at compile time and removes the corresponding code.The compile time check of the basis functions types is then crucial to get the previous g++ optimization. 8s mentionned in the introduction, DSEL inherits the capabilities of the underlying language; we have shown here a very good illustration of this statement.Now that we know how to do integrate a variational form on a convex, it is easy to extend it to a set of convexes.

Automatic differentiation
Automatic differentiation as described in [3,4] operates differently from the previous algorithms when evaluating an expression template.We are no longer evaluating an expression at a set of points but rather dealing with expressions that manipulate a differentiation numerical type ADType<P,N> -P is the number of parameters and N the order of differentiation -holding a value and corresponding derivatives -in pratice up to order N = 2 and P 100 -To implement the automatic differentiation engine, we construct the expression object as before.However the evaluation is changed: in our case we create a new subclass for evaluating the differentiation.Each class which could be possibly used for automatic differentiation need to implement this subclass.Listing 19 shows an excerpt of the implementation of the expression templates glue class Expression.

Listing 17: Evaluation subclass modifications for Variational Formulations
The automatic differentiation evaluation engine does not need to be much further discussed as it follows ideas already published elsewhere, see [4].We just demonstrate here that decoupling expression construction from its evaluation allows to share the code constructing the expression template object.

A Remark on mixing automatic differentiation and variational formulations
We could actually push the envelop quite a lot and use the automatic differentiation type within a variational formulation to differentiate with respect to some parameters for sensitivity analysis, optimization or control of engineering components.That is to say, mix the two evaluation engines : first the integral expression is evaluated and during the integral evaluation, the automatic differentiation language takes over to finalize it.The enabler of the second stage in the evaluation would be the automatic differentation data type and its operator=( Expr const& ).This is truly what meta-programming is about -code that generates code that generates code . . .-However this has not been tested in the examples shown in Section 4 and in particular performance could be a concern depending on the quality of the generated code.This is one of the topics for future research as it may open important areas of applications for the language such as the ones mentionned previously -sensitivity analysis, optimization and control.
This leads to another advantage of sharing the expression object construction: it ensures that we have the same unit outward normal at the current node of the current face Nx(),Ny(),Nz() x, y and z component of the unit outward normal P() coordinates of the current node Px(),Py(),Pz() x, y and z component of current node set of supported mathematical functions or functors for automatic differentation on the one hand and projection, integration on the other hand which means less development work and less bugs.

Overview of the language
The implementation of the language itself, the construction of the expression object is done in a very standard way conceptually, see [24].However from a technical point, it uses state of the art tools such as the Boost.Mpl [2] and Boost.Preprocessor [2,16].
In particular the Boost.preprocessorlibrary is extremely useful when it comes to generate the objects and functions of the language from a variety of numerical types and operators -using for example the macro POOST PP LIST FOR EACH PRODUCT.
Regarding the grammar of the language, it follows the C++ one for a specific set of keywords wich are displayed in Tables 1, 2 and 3.The keywords choice follows closely the Freefem++ language, see [12].

A note on the supported numerical types
The language supports the standard ones available in C++ -boolean, integral and floating-point types -and some additional ones std::complex<> complex data type -ADtype<.,.> the automatic differentiation type -dd_real double-double precision data type from the QD library, see [13] -qd_real quad-double precision data type from the QD library -mp_real arbitrary precision data type from the ARPREC library, see [6] Other data types can be relatively easily supported thanks to Boost.Preprocessor.There are plans for example to support an interval arithmetic data type -the one provided by Boost.Interval, -which would be interesting to measure, for example, the impact of small perturbations or uncertainty on the results.

Test cases
We present now various examples to illustrate the techniques developed previously.We shall show only the relevant part of each example.Also, for each example, we shall present different possible formulations.

Performance
In the following section, various results are presented with some timings for the construction of bilinear and linear forms.The calculations have been conducted on an AMD64 opteron9 running linux 2.6.9 and GNU/Debian/Linux.The g++-4.0.1 compiler was used to compile the code with the following options Listing 20: g++ options for optimization These options are relatively standard and give good results all the time.Other optimization options have been tried but did not yield significant performance improvements.
We consider various Advection-Diffusion-Reaction problems to evaluate the performance of our framework proposed in [20].
µ∆u + σu = 0 DR, (16) µ∆u + β∆∇u + σu = 0 DAR (17) on a unit square and cube domain.For each problem we consider both the case of constant and space-dependent coefficients.The following expressions were assumed for the space-dependent coefficients µ(x, y, z) = x 3 + y 2 (2D) σ(x, y, z) = x 3 + y 2 (2D) The corresponding C++ code for the 3D is presented in listing 21.The C++ for the 2D cases is very similar but simpler.

Results
The benchmark table, see 4, reports the timings for filling the matrix entries associated with the problems D, DR and DAR.Lagrange element of order 1-3 dof in 2D, 4 in 3D -and order 2-6 dof in 2D, 10 dof in 3D -have been tested.The integration rule changes depending on the polynomial order we wish to integrate exactly.For P1 Lagrange elements 3 quadrature points are used in 2D, 4 in 3D and for P2 Lagrange elements 4 are used in 2D and 15 in 3D.

Analysis
In order to facilitate the study of these timing results, they are displayed on Figs 3 and 4 for the 2D cases and 3D cases respectively with each time the matrix assembly time with respect to the number of elements along with the ratio between the cst timings and xyz timings for P1 and P2 cases.We can observe a few things: -matrix assembly time versus number of elements is linear (in log − log) which is no surprise, -the difference of performances between the different equations D, DR and DAR are quite small even with the DAR and the 3D cases which means that we do a good job at sharing as much computations as possible between the different terms of the equation, -the overhead due to non-constant coefficients is very small in all cases.which is not the case for example in [20].
In particular, as shown by the ratio figures, the ratios are always less than 2 and in most cases -among them the most expensive ones -they are in [1.1; 1.3].In [20] they used functions or functors to treat the non constant coefficients, they get factors between 2 and 5 (5 in the worst case 3D P 2 ).This means that the expression templates mechanism and the g++ optimizer do a very good job at optimizing out the expression evaluation.
These results illustrate very well the pertinence of using meta-programming in general and expression templates in particular in demanding scientific computing codes.Also note that there was no particular optimization based on some knowledge of the underlying equation terms with respect to the local matrix assembly, so there is room for improvements -for example the symmetry of the bilinear form or the precomputation of stiffness, mass and convection matrices on the reference element, see [17].Regarding compilation time, it is certainly true that using expressions template and meta-programming has an impact on the compile time.However if most of the library is templatized then the compile time cost is usually paid when compiling the end-user application and no more when compiling the library itself.To give an idea, the performance benchmark code was compiled with all cases: P1 and P2 in 2D and 3D for all problems D, DR and DAR.In other words, a single executable was generated to run all the cases presented earlier.On an AMD64 opteron,10 the compilation takes between 2 and 3 minutes using the compiler options from listing 20.The book by David Abrahams and Aleksey Gurtovoy [2] provides a complete discussion on meta-programming and its impact on compile time.

A variationnal inequality
This test case is described in [9].We consider a rectangular tank of length a = 0.1 m, and height b = 0.05 m.A cylindrical tube crosses it with a diameter of 0.015 m.The tank is filled with ice and there is a thin layer of solid/liquid polymer on top of it.For symmetry reason, we consider only half of the tanka = 0.05 m. -The problem is formulated as follows: The temperature θ is solution of the parabolic equation where and ε = 0.05.β(θ) is the volumetric heat capacity and µ(θ) is the thermal conductivity.Finally, χ represent the characteristic function.
We impose the following conditions: θ = 0 on boundaries 1 and 2 -∂θ ∂n = 0 on boundary 3 and 4 θ = −0.1 + 0.05t on boundary 5 where t represents the time The actual code is presented in listing 22.
Figure 6 shows the contour lines of the temperature at various time steps and the mesh colored by the thermal conductivity over the entire domain: we see that the ice is melting and transformed into water -in dark grey the ice and light grey the water.

Stokes and navier-stokes
We consider here the standard driven cavity test case with the following setting described by the Fig. 7. First we use a stable mixed approximation for the velocity and pressure spaces -say the Taylor-Hood element (P 2 -P 1 ).-The variational formulation reads as follows: Find (u, p) ∈ X h × M h such that for all (v, q) ∈ X h × M h We can also use an equal order approximation -say P1-P1 -and add a stabilization term like the jump of the pressure gradient over the internal faces as proposed in [8].The formulation reads then as follows: Find

Particule in a shear flow
We consider now a rigid particule in a shear Stokes flow treated using a penalization method, see [14].Denote Ω = [−1, 1] 2 , we consider a cirular particule positionned at (x p = 0, y p = 0) ∈ Ω of radius r p = 0.1 and a penalization parameter ε.Denote χ p the characteristic function of the particule defined as χ p = χ((x − xp) 2 + (y − yp) 2 > r p ).We seek (u, p) ∈ X h × M h such that ∀(v, q) ∈ X h × M h Figure 8 shows the velocity vector field and identifies the particule in the mesh using its characteristic function χ p .

Conclusion
We have developed a unified DSEL for different aspects of numerical analysis, namely differentiation, integration, projection and variational formulations.The main results of this article are that (i) such a DSEL is feasible: an implementation has been done and exercised with some non-trivial test cases and (ii) decoupling the expression object construction from its evaluation using a delegate subclass allows for very powerful notations in C++: one unique engine is used for the expression construction and as many engines as needed for the expression evaluation -we have seen two different engines: one for projection, integration and variational formulations and one for automatic differentiation.-This technique can certainly be successfully used in other contexts.
Future developments will include a study whether the work done in [17] can be applied to accelerate the assembly steps during integration.Also, although vectorial notations in the language can be used, they need to be formalized within the language to avoid ambiguities that would yield wrong results.

Listing 1 :
Variational Formulation in C++; the t extension of gradt and idt identifies the trial functions

Listing 5 :
Approximation Space Interface Listing 6: Mixed Space Example

Listing 8 : 9 :
Evaluation Delegation to SubclassListing Current Evaluation Point Coordinates class that returns the x coordinate of the points where the evaluation is effected.Its implementation is presented in listing 9.

Listing 10 :
Current Evaluation Point Coordinates

Listing 13: integrate prototype Listing 14 :
Integrator class and integrate free function Listing 15 shows an example of the syntactic sugar brought by the language.

Listing 21 :
Benchmark for elliptic problems

Fig. 3 .
Fig. 3. Cases P 1 and P 2 in 2D.Matrix assembly time versus number of elements and ratios between non-constant coefficients assembly and constant coefficients assembly.

Fig. 4 .
Fig. 4. Cases P 1 and P 2 in 3D.Matrix assembly time versus number of elements and ratios between non-constant coefficients assembly and constant coefficients assembly.

Fig. 6 .Fig. 7 .
Fig. 6.Temperature and thermal conductivity in the tank at various time steps.

Algorithm 2
Nodal projection on X ι X is the local/global correspondance table for T ⊂ T do pi = 1, . . ., N X points coordinates associated with the DoF in T c ← {T, G, (p i ) i=1...N X } geometric mapping context, see 2.1.2for i = 1, . . ., N X points coordinates associated with the DoF do

Table 1
Tables of geometric/element operators available at the current element, face and node

Table 2
Tables of mathematical functions that can be applied to expressions