Object-oriented design to automate a high order non-linear solver based on asymptotic numerical method

The Manitoo library is devoted to the resolution of analytical non-linear problems using a high order method called asymptotic numerical method. We describe here the Object Oriented design of this library and especially the choices made to obtain a quite generic and ﬂexible numerical solver. Through classical examples, we present a comparison with some existing tools implemented in Matlab and Fortran 77


Introduction
For about twenty years, many works on the resolution of nonlinear problems using Asymptotic Numerical Method (ANM) have been proposed [1] and applied to a wide range of problems in fluid and solid mechanics. For instance, various applications concerning the design of marine structures [2], non-linear vibrations [3], sheet metal forming [4], biomechanics [5] or multi-scale instabilities [6] have been implemented in many research informatics tools using Fortran 77, Fortran 90 [7,8], Matlab [9,3] and also in an industrial code [10]. Within ANM, a solution branch is computed using series expansion. The step length of this branch is then defined a posteriori and it yields naturally automatic path following techniques. This automation of the non-linear computation is very important in many cases and especially for the numerical computation of physical problems involving instabilities [11][12][13][14][15][16]. The robustness of the algorithm is assessed for instance in [2], where thousand thin shell computations have been performed in order to predict thin shell buckling with random imperfections. In Computational Fluid Dynamics, the discretized problem often involves millions of degrees of freedom so parallel implementations are needed.
The robustness and efficiency of ANM have been established to solve such large scale problems [17].
ANM needs high order derivatives that can be obtained by recurrence formulae. These recurrence formulae computation is easy for algebraically simple equations, as Navier-Stokes equations, but it becomes more intricate for many other physical problems [18]. A first answer to this question is the MANLAB software [9] which permits to solve generic problems with few unknowns, providing the problem written in a quadratic form. Another approach, called DIAMANT and based on Automatic Differentiation (AD) by operator overloading, has been recently proposed [19,20] and applied to small size academic problems in the Matlab and Fortran contexts [21,22]. Due to Object Oriented limitation of these languages [23], it seems then difficult to obtain a really generic, reusable and efficient ANM library. Moreover, as for most of informatics libraries, and based on a 15 years old experience, we know that the ANM library mainly requires maintainability, a wide extensibility and portability. Indeed, as mechanical models based on finite elements (or other approximation methods) often need a large number of degrees of freedom, we have to consider the use of parallel computing.
Since the 90's, Object Oriented Programming has been commonly used to design complex scientific applications. As mentioned in the context of finite element method by [24], object oriented programming permits to develop numerical tools with portability on different computer architectures such as clusters, which was not possible with the prior sequential Fortran codes. Moreover, inheritance and polymorphism are important characteristics to obtain a quite generic solver. One could argue for the loss of performance due to the overhead of oriented object implementation. Such a loss is not so obvious and Veldhuizhen demonstrates that C++ could surperform Fortran [25,26].
Despite of the propagation of object oriented and C++ programming in the field of numerical simulation of engineering problems, the ANM community has not developed an up-to-date numerical tool.
Considering all these remarks, we have been developing a new C++ library, called MANITOO which is devoted to the resolution of non-linear problems with ANM, since 2008. The main goal of this tool is to reduce development costs without losing computational performance compared to the former library developed in Fortran 77. Here we propose to deal with the object oriented design of Manitoo which has never been published and is about to be mature.
In the second part, we make a description of the ANM and present the corresponding algorithms while the third part is devoted to the Object Oriented design of the library. The fourth part shows some applications and comparisons with existing tools. We conclude with some remarks on future developments.

General structure
ANM consists in solving an analytical non-linear problem with a path-following (or continuation) method associated with a high order perturbation technique. First, unknowns are expanded in Taylor series with respect to a scalar path parameter. Then the non-linear problem, as a problem depending on unknowns, is also expanded in Taylor series leading to a system of linear equations at each order. Expressing high order equations with respect to Taylor coefficients requires a hard programming work. Coupling ANM with AD (namely the DIAMANT approach) allows to automate the computation of Taylor series terms in a peculiar direction using well-known recurrence formulae.
As in [27] and to obtain a generic library, the model describing the non-linear problem is distinguished from the analysis. Moreover, we split the analysis in ANM solver and linear solver.  We focus here on the design of the ANM solver therefore few indications about model and linear solvers integration will be given.

General overview
First works concerning Asymptotic Numerical Method (ANM) referred it as the computation of a solution branch of a non-linear analytical problem [1]. It has been extended by [28] into a continuation method resulting in an assembly of successive solution branches. We are using the term of ANM as an alias of path following technique via an asymptotic numerical method based on power series expansion in the sense of [28].
Here is a brief overview of an ANM version summarized from [29]. More interested users may refer to this book.
Let us consider a quite generic quasi-static problem expressed as:   with U a convex domain of finite dimension and R a non-linear analytical function. Damil et al. [1] used a high order perturbation technique to solve this kind of problem by looking for solutions (u, k) as a function of a parameter denoted by t close to a given solution (u 0 , k 0 ): This can easily be related to the power series expansion and Taylor series terms of u and k expressed as functions of the parameter t: Hence the computation of a solution branch is equivalent to the computation of the expansion coefficients. The latter are obviously related to the high order derivatives of the analytical function R. Every derivative R k can be split in a linear part depending on the unknowns (u k , k k ) and a non-linear part depending on the previously computed terms. To compute a solution branch up to a given order n, one has to solve the following successive linear problems: 0 ; k 0 Þk n ¼ÀR nl n ðu 0 ; u 1 ; ...; u nÀ1 ; k 0 ; k 1 ; ...; k nÀ1 Þ with L u T ¼ @R @u the tangent operator with respect to u; L k T ¼ @R @k the tangent operator with respect to k and R nl n the non-linear part of R n .  As this set of equations is subdertermined, a path-parameter equation is added and can be expressed in a quite generic form: where hÁ, Ái u is a bilinear form defined in U and a arc is a positive coefficient. Many other parameters are possible, see for instance [30].A complete discussion on this point can be found in [31].
Usually, the classical dot product is applied and a arc equals 1. A change in the bilinear form and the weight leads to a different path drive.   The derivatives of the path parameter Eq. (2) are given by: Once the series terms are obtained, on evaluates the interval, namely the validity domain, in which the series converge. Hence the length of the solution branch is known and the end point of this branch is the start point of the next branch.
The Algorithm 1 summarizes the resolution of the non-linear analytical problem 1 with ANM and continuation. Padé approximants have been introduced in the context of ANM by [32]. Their efficiency has been shwon in terms of extending the validity domain and consequently on reducing the number of branch computation steps [33,29].
Let (U k ) k=0,...,n and (k k ) k=0,...,n be series obtained by solving an analytical non-linear problem using ANM. From these series, a validity domain can be a posteriori computed. Using Padé approximants allows to increase the length of the solution branch by building up an orthonormal basis via the Gram-Schmidt procedure. Starting from solutions given in the form: one can obtain a new representation by rational fraction of the solution path: where the expression of D i (t) is given in [33,29]. Accounting for Padé approximants results in adding an external plug after a series computation. The increase of computation time due to this plug is, in most of cases, much smaller than the benefit due to the increase of the length of the solution branch. The global Algorithm 1 is then replaced by Algorithm 2.
Compute high order term R nl

Global architecture
A previous work on Object Oriented implementation of a pathfollowing method devoted to structural mechanics has already been proposed [34]. Our new library, called MANITOO, applies to a wider kind of analytical non-linear problems. Moreover we attempt to define a more generic implementation of ANM which is also a path-following method.
The MANITOO library needs to consider the overall ANM process: Taylor terms, validity domain computation and continuation technique.
From the previous section, it seems obvious to split the ANM in several components: Path-following process Branch solution building Linear system solving Non-linear problem definition Derivatives computation (optional).
Only the first two parts are really specific to the ANM. The way of computing high order derivatives, defining the model and solving linear systems should not influence the ANM process for a given kind of non-linear problems. Then a generic ANM algorithm can be defined by Algorithms 3 and 4.
The path-following flow consists in computing successive solution branches up to a given number of steps or an optional loopbreaking condition (Algorithm 3) while the solution branch is built by Algorithm 4. So we distinguish the analysis (ANM specific features) from the model (non-linear problem expression), the way of differentiation (derivatives terms) and the linear solver. This leads to four packages related as in Fig. 1.

Finite Element Method
In the context of solid mechanics, the Finite Element Method (FEM) is one of the most popular numerical method to approximate the solution of partial differential equations. This method is here used to model an approximation of the non-linear problem.
It is quite impossible to establish an exhaustive list of all existing Object Oriented (OO) libraries dealing with Finite Element Methods (FEMs) or other discretization methods. A bibliography has been realized until 2003 in [24] and shows the amount of works of FEM and OO.
Thus, behind the term of ''Finite Element library'', one could find several kind of libraries. Some of them are related to the mesh building and entity organization (AOMD [35], libmesh [36]), some of them are just related to the implementation of the computational method (OOfem [37], getfem++ [38]), others provide an upper interface and pseudo-language so that the user only has to define his problem on a strong or variational form (freefem++ [39], Fenics [40]), etc.
The library catalog of Object Oriented Finite Element is really wide so we hoped to find a convenient one for our applications. However, in the context of continuum mechanics, we want to deal with linear and quadratic 2D and 3D elements. This first requirement reduces the number of potential candidates. A second requirement is to use a really specific shell element from Büchter et al. works [41] which gives good results for very thin shells. With this second requirement, at our knowledge, the potential candidate space is a null one. Thus we decided to develop an Object Oriented Finite Element Model with basic features like iso-parametric interpolation, numerical integration, element to global assembly, etc. Our model may not offer the same level of advanced functionality as in [37,[42][43][44] and, consequently, will not be detailed in this paper. The 2D and 3D models are quite similar (it just results in changing a template argument in most cases) while dealing with the cited shell elements needs the introduction of new internal variables, condensation technique and tensor algebra. This results to a specific class considering models based on shell elements.

Linear solver
Solving a non-linear problem results in solving successive linear systems with the same tangent matrix and changing the right hand side. Wrappers to the following solvers have been implemented: JAMA [45] SuperLU [46] PetsC [47][48][49] Mumps [50,51].
As the tangent matrix does not change during the successive system resolutions, we use features of these solvers to optimize the computation cost of the linear system resolution (please refer to user manuals of these libraries to get more details).
Tnt library is used for all our validation tests with small degrees of freedom and a dense linear matrix.
Three kinds of sparse matrices storage are defined to use Super-LU, PetsC, and Mumps.

Automatic Differentiation (AD)
As noted in the ANM overview, solving a problem with this method results in computing high order derivatives in a peculiar direction. A good way to automate the ANM solver is to use automatic differentiation as a pocket calculator. In this sense, an approach, named DIAMANT, coupling AD and ANM has been proposed in [19] and applied to academic problems [20]. Then righthand side terms and tangent matrix occurring in ANM are obtained by application of AD technique. It was already an improvement avoiding the manual calculation of terms which could be a quite difficult task.
Following the definition of [52], AD can be viewed as the way to generate a code which computes the derivative of a function itself given by a code. A natural approach is to use the forward mode by applying operator overloading. This is notably well-suited with C++ language [53]. This language enables to overload classical operators in a simple way [54]. Notice that AD by operator overloading has already been implemented for a long time in other classical languages [55].
Interested readers may refer to [56] for a review of AD techniques. But in most cases, works are concerned with low order derivatives (first or second order). We are here concerned with high order derivatives in direct mode as in [57] except that we only need to compute the derivatives in a peculiar direction (i.e. the Taylor series terms of a given variable or function).
From an informatics point of view, forward differentiation by operator overloading is realized using template facilities of C++ language. First, the function to be differentiated is implemented and tested in the continuous domain, after that the template parameter is changed into the high order derivative one. Consequently, one single implementation allows the computation of residuals and derivatives. A specific procedure is applied in the case of local implicit functions because of the different nature between the continuous function (implicit) and the derivative function (explicit).
The AD feature is also applied to the computation of the linear tangent matrix, noted L u T in Section 2.1.1 arising from the first order differentiation of the non-linear problem.

Object oriented implementation
A non-exhaustive overview of the implemented classes, restricted to this paper focus and organized by packages, is given in Fig. 2.
In this section, we describe some of the techniques used in the library. The objective is to split the overall process into small size objects to obtain the higher flexibility and genericity.
As explained in Section 2.1.1, we want to apply an iterative process to the computation of algebraic functions. One important point is the way of considering mathematical/ mechanical functions in our library. As mentioned in [58], function objects (also called functors) allow better optimization than function pointers in terms of computation time. However these functors have to be small enough to be inlined. They should also not be virtual. Functors are intensively used in this work to implement small functions. Moreover, to keep a reasonable computation time cost, we use meta-programming and limit the use of virtual functions.
Note that the library is implemented in french, so names of components, classes, operators have been translated for this paper.
Note also that all the UML diagrams have been built with Bouml freeware (http://bouml.free.fr).

Global class diagram
An ANM solver is a way of computing output data, usually named results, from input data given by a Model following a fixed algorithm. From the package view in Fig. 1, we focus on the ANM package. This package is made-off three main classes.
The path-following class The branch solution class The path equation class.
An overall class diagram is illustrated in Fig. 3. Referring to the class diagram in Fig. 3, user-defined class T_ProblemNL describes the non-linear problem and the computation of its high order derivatives (Taylor series coefficients). It is defined as a functor.
Using a templated class, here called T_LinearSystem, allows us to apply any solver of a set of linear equations.
Moreover, it seems obvious that we need to store any Taylor series terms and optionally to define devoted algebra rules to these data. Hence we defined a new data type.

Data Types
This subsection is concerned with data types of the Taylor Series variables (i.e. derivatives data type). As indicated in Section 2.5, automatic differentiation in forward mode could be done by overloading classical operators. So we would design new data type containing Taylor series coefficients and the corresponding arithmetic operators (and usual mathematics functions). An illustration of automatic differentiation by operator overloading is given in Section 3.6.
Moreover, solving non-linear problems in the field of continuum mechanics leads to using scalar variables and some containers related to linear algebra as vector or matrix. Here the term ''vector'' is related to the linear algebra term and not the STL term.
We make an intensive use of template ability offered by the C++ language so that a change in the scalar data type and in the container type could be done with only a few change in the library. To do so, we define a traits policy [59] containing: Scalar data type Vector data type Dense local matrix data type for computation with matrix of small size Linear tangent matrix data type which is usually a sparse one. that could be implemented as: Note that the type of the large matrices depends on the external linear solver. Some of them require Compressed Sparse Row matrix while others need Compressed Sparse Column, skyline or full storage. Moreover, we use our own implementation of vector and matrix container. Replacing by more efficient containers would be considered in future improvement of the library.
A prior feature of the differentiate data type is to contain all the terms of the series expansion. It has been chosen to define a data using static array. As a computation needs many calls to constructors/destructors, it is obvious that using dynamic allocator would lead to a loss of performance. The size of the array is defined using a macro variable. Taylor expansions are computed up to a given order less than the maximal array size.
Moreover, algebra could be added to the data class with Operator Overloading technique; this corresponds to an automatic differentiation in forward mode. Three classes of high order derivatives are available. Performances and features of these are described in [20]. As an illustration, a brief view of the class DType1 corresponding to the first version of DIAMANT approach is given in Fig. 4.
Note that due to the size of the manipulated data, we avoid at most to make a copy of scalars, vectors and matrices. Only one tangent matrix is created.

Path-following
As mentioned before, continuation consists in computing the validity domain of the series expansion, building the end point of the branch, then using this point to start the computation of a new branch.
Moreover, the matrix computation object, the linear solver object and optionally the residual computation object are created at the initial step of continuation process.
The sequence Fig. 5 outlines the process of continuation.

Branch computation
Branch computation consists in finding terms of the Taylor series expansion (i.e. high order derivatives terms). Knowing the initial solution (U 0 , k 0 ), Taylor series terms up to an order n are obtained.
As for continuation, this class is perhaps the most generic. However it needs some specifications depending on the kind of non-linear problem. Then a virtual base class has been designed (Fig. 6). In term of computation time, this should not be dramatic. Indeed, on non-academic simulations with a large number of degrees of freedom, access to the virtual table, also called vtable, is really less time-consuming than solving each right hand side term at each order and linear tangent assembly.
The sequence diagram related to the computation of the series terms is given in Fig. 7. Applied class would be expressed in the next section dealing with examples.

Padé class
As highlighted in Section 2.1.2, Padé approximants could be designed as a computational class. A simplified version of this class is described in Fig. 8. Members named Orthonormalization, SerieSum and Solution_Pade are private ones used to compute the validity domain and coefficients related to Padé method.
The sequence Fig. 9 illustrates the easiness of including a posteriori computation of Padé approximants and the extended validity domain.

Formulations
Formulations describe the non-linear problem. It computes function results with data and parameters given in the model.
A formulation has to be defined as a functor class even if some of them are not so small. Then formulations classes need at least parenthesis () operator. However, we want to compute the function and its derivatives up to a given order. There are at least two ways of implementing the function operator. The first corresponds to the classical hand-written developments. Each algebraic operation is developed by hand and, sometimes, smart and really specific optimizations could be investigated before the computer implementation.
The second way is to apply operator overloading based on the fact that a function is a compound of elementary and well-known operators. This leads to compute high order terms of a given function using AD technique in implements only once adequate formulae. This ''pocket calculator'' ability allows us to obtain more genericity in our library. The user should only implement computation function and then, with just referring to the appropriate data-type, one could compute residual value or differentiate terms.
The designed ANM library is free from the way of computing high order terms so that each user could choose his own preference. One has just to define a formulation as a functor with a good data type. However, independently of the computational way of high order derivatives, it is strongly recommended to enable different kinds of data treated by the formulation. Hence, with changing the data type, one could compute high order derivatives or evaluate the residual of a given solution. As an illustration, the function f(x, y)=x Ã y could be naively implemented with the following templated ScalarProduct class: It is quite easy to compute the function value in the continuous domain with automatic differentiation or with the hand-coded expansion: and its derivatives: Switching between scalar and derivative computation amounts to changing the variable type. With AD, recurrence formulae are implemented in the T_Diff class with operator overloading ability. It seems obvious that a hand-coded version leads to separate operations on derivatives from data so that any formulations need to rewrite associated recurrence formulae.
Coupling the automatic computation of formulation with our quite generic version of ANM allows us to claim that we obtain a quite generic and automatic version of a non-linear solver based on an asymptotic numerical method. We briefly discuss the performances and illustrate the implementation of a new formulation in the incoming application part.

Applications
In this section, we propose to solve two kinds of problems with our library. Each problem has different mathematical features and illustrates flexibility, reusability and genericity of MANITOO.
Computations have been performed under a 64 bit Windows 7 operating system on a Lenovo x200 Tablet. The processor is an Intel Core2 Duo L9400 @ 1.86 GHz and 2.9 Go of RAM are available.

Problem formulation and discretization
This example was solved in the ANM context by [29]. It was the first problem solved with the here-concerned MANITOO object-oriented library. The initial formulation of the problem is given by: The second formulation of the Bratu problem has been introduced to simplify the expansion made by hand. Hence in [29], a new kind of generic problems has been designed with a main equation and a local one expressed in a differential form. Using AD enables to work with the first formulation without any introduction of new variable or specific smart expansion. In the following part, we will deal with the second formulation, as the first is less complicated.
Concerning the discretization and numerical approximation of second order derivative in the spatial domain, the finite difference method is used to compute approximate the solution (u h , v h , k h ).
The domain C is discretized with a constant step size h.
ANM is applied on this discretized problem given by: Computing the tangent matrix respect to u h is obvious:

ANM
The unknowns are sought in the form: t i k i and the path-parameter is defined as: This problem could then be solved using Algorithm 5. At each order, the right hand side comes from the contribution of non-linearities, and then, in this problem, to the computation of the product kv.
Note that the variable v is not really an unknown but an internal variable (also called condensation variable). This variable is used for continuation at each branch computation. The algorithm could be improved for a hand-written differentiation (refer to Section 2.2 of [29]) but here is not the purpose.

Implementation
To implement the Bratu problem, we do not have to worry about condensation technique.
ANM is compound by a continuation class, a branch solution class, a linear solver class, the specific BratuProblem and BratuMatrix classes devoted to the computation of the function R(u, k), the high order derivatives and the tangent matrix with respect to the u variable. In this case AD is not used to compute the linear matrix.
The continuation, branch solution and linear solver classes are generic ones available in the library. As noticed from subSection 4.1.1, the Bratu problem is of the form R(U, k) with non-linearities in (U, k). Hence we use the BranchSolution_U_Lambda class. For this example we decide to use the SuperLU linear solver. The corresponding sparse matrix storage is defined to build the tangent matrix.
In the following, we give an implementation example of Bratu-Problem and BratuMatrix classes described in the diagram Figs. 10 and 11. BratuProblem class is an implementation of the discretized problem (9) while BratuMatrix is an implementation of the tangent matrix computation (10) (without condensation technique). Results are obtained with a call to the operators parenthesis returning a reference on a T_Vector in case of equilibrium computation (BratuProblem): As T_Scalar and T_Vector are template parameters of the class BratuProblem, on can choose to.
Compute the constant value by using for example the type double for T_Scalar and Vecteurhdoublei for T_Vector Compute the differentiate terms (or Taylor series terms) by using for example the type DType1 for T_Scalar and Vec-teurhDType1i for T_Vector.

In case of matrix computation (BratuMatrix), the result is a reference to a T_Matrix:
where the template parameters T_Vector, T_Matrix can be defined to choose the vector container type (for example Vectorhdoublei) and the type of linear matrix. Latter has to be coherent with the linear system solver called from ANM.
Note that an object of the class BratuProblemhdouble, Vectorh doubleii allows to compute the residual for a given solution through its parenthesis operator.

Results
To test the efficiency of MANITOO, we compare our implementation with a Matlab code. This code corresponds to a hand-written asymptotic numerical method. As the operations involved in the Bratu problem are few and quite simple, it is expected that the Matlab code is optimized in term of runtime performance. Moreover we use this example to evaluate the performance of automatic differentiation compared to hand-written ANM.
The classical solution curve at the middle of the domain is given in Fig. 12. Each point on the curve illustrates the end of a solution branch and the beginning of the next one.
Obviously, the finite difference method applied to Bratu problem leads to a really sparse matrix. In this case, we use the superLU solver with sparse matrix to solve linear systems in the Manitoo library while the sparse keyword is used in Matlab. Table 1 describes computation time (in s). The total and linear-solver cpu-time of Matlab, hand-written ANM and ANM with automatic differentiation are reported. Total time is the time spent to solve the problem while linear-solver time is the time spent to solve for each branch at each solver the linear systems KU k = F k .
We observe that the Matlab computation time seems to be greater than the computation time needed by AD or classical ANM. As the matrix is really sparse, we observe that AD time is more than twice of the hand-written code. Due to the use of condensation technique in the Matlab code, it could be viewed similar to the second version of DIAMANT [20] (i.e. there is no recomputation of the previous order) while we are here using the first version of automatic differentiation: for a given order we recompute all terms of the series up to this order. This kind of AD avoids us to store all intermediate variables. Moreover we do not need specific AD techniques such as check-pointing or sequence recording. This leads to a better flexibility and genericity of the algorithm. Note also that as observed on this case, memory management is a critical point for engineering problems.

Formulation
Let us consider an elastic deformable body. We look for the displacement field, denoted by u, of this body subjected to an external force.
Usually, the Green-Lagrange strain tensor c is defined by where h(u)=r x u is the displacement gradient tensor. In elasticity framework, the second Piola-Kirchhoff stress tensor, denoted S,i s related to the Green-Lagrange strain tensor by the fourth order elastic tensor L: The problem to solve is expressed as: In this example, non-linearity only lies in the expression of the strain field.
The problem (13) is discretized using the finite element method with 2D, 3D or shell elements. A Gauss quadrature approximation is used to integrate the non-linear form over the discretized do-main. Hence the stress tensor is a local variable evaluated at each integration point (or Gauss point).

ANM
The unknown are sought in the form: t i k i and the path-parameter is defined as: t ¼huðtÞÀu 0 ; u 1 iþðkðtÞÀk 0 Þk 1 : One has also to differentiate the intermediate variables: and the global formulation: using This problem could then be solved using Algorithm 1. At each order, the right hand side comes from the contribution of geometric non-linearities. According to the Bratu example, one could condense intermediate variables, such as the stress tensor. This leads to the convenient quadratic formulation [29]. In Manitoo, we recall that we work with any kind of formulation.

Design of the formulation class
As for the Bratu example, ANM is compound by a continuation class, a branch solution class, a linear solver class and classes defining non-linear problem and matrix computation.
In this example, the problem (13) could be written as R(U) À kF =0. Hence we do not need to compute the tangent term L T k related to the k variable. This term is here always equals to À F. Then we are here using the BranchSolution class accounting for this specificity. The continuation process does not change from the Bratu problem.
The non-linear formulation consists in computing an approximation of R(U) over a discretized domain. As the process of finite element computation is always the same, we design a generic function that iterates on all the elements of a mesh, apply on each one the computation of a given function, and assembles the element results to the global result. This function is called by a class named GlobalEquilibriumFE. It computes a function overall a finite element domain through the parenthesis operator.
Then we have to implement a class, called for example FML::InternalWorkElt, computing the numerical approximation of the internal work on an element. Concerning the computation of tangent matrix, even if it could be easily obtained and computed by hand, as explained in [29], we use the automated computation of this matrix. This allows us to change for example the constitutive law in the formulation without any change in the implementation of tangent matrix computation. As for the computation of non-linear formulation, computation of tangent matrix is quite generic. On each element a local tangent matrix is computed and assembled into the global matrix. The local matrix is computed with successive evaluations of the first derivative of element equilibrium equation on the canonical basis directions.
We give hereafter a simplified example of the element internal work computation class in a 2D case. This is a non-generic implementation of this kind of class and it is not the one used for the following results subsection. However the purpose is here to illustrate the relative programming effort to implement a new formulation at an element level.

Results
MANITOO is used to solve an academic buckling example dealing with the snap-through of a thin curved roof as in [12]. We consider that the behavior of the deformable body is linear and elastic. Shell elements are used to discretize the roof as in [12].
The mesh used, called mesh22 Â 22, is made up of 484 elements, 1541 nodes and 9246 degrees of freedom.
On Fig. 13 we compared the results obtained with MANITOO, Abaqus and Eve (a Fortran77 implementation of the ANM). On the Manitoo and Eve curves, points indicate the end of a branch and the beginning of the next one. Consequently, each point corresponds to the computation of a tangent matrix and its factorization.
On the Abaqus curves, each point corresponds to an increment computation and consequently to a linear system solution. In this case, each point is linked to the next one by a straight line. Accuracy of Abaqus results depends on the size of the increments and then to the number of matrix inversions. From Fig. 13, we see that when we need 13 matrix computation and factorization with the ANM, Abaqus needs about 75 matrix inversions to obtain a quite accurate curve.
However, the comparison of these two softwares is not really appropriate because MANITOO has a specific goal: ANM was designed to build analytical solutions that would be post-processed.
As MANITOO has been designed and developed to replace the Eve code, we compare the computation time. Simulation with Eve requires 2254 s while Manitoo only needs 1957 s.
We observe a saving of about 13%. Comparing results on this academic mechanical case does not lead to a global conclusion. However, this comparison ensures us that at least on this basic case, object oriented C++ implementation coupled with automatic differentiation overperforms the available Fortran77 implementation of ANM. We have obtained a faster solver with more genericity and flexibility than the Eve version.

Concluding remarks
This paper deals with the object oriented design of MANITOO which is an automatic non-linear solver based on asymptotic numerical method. We try to build a quite generic and flexible library by using some concepts such as functors, traits and classical Object Oriented paradigms.
Due to the nature of the mechanical problems, design of MAN-ITOO integrates a wrapper to numerous linear solvers. Moreover we have quickly designed our own finite element library due to the lack in the opensource community of a library accouting for our specific shell elements.
As illustrations we have applied our library to solve two different kinds of academic problems. On the one hand, we note the quality of the obtained solution from the solver which is a necessary condition to go on with this library. On the other hand, by comparing our library with Eve, a Fortran77 implementation of ANM, on a simple mechanical problem, we note that MANITOO allows to decrease the computation time. For our application, the object oriented design and C++ implementation is not slower than the Fortran77 one.
As a final result, we achieve to design, at our knowledge, the first generic non-linear solver based on Asymptotic Numerical Method with performance at least equivalent to the Fortran implementation existing in the ANM community.
Of course more works have to be realized to obtain a more complete library. Firstly, coupling MANITOO with an automated Finite Element tool such as the Fenics project [60] could result in a really generic finite element solver and could be of great interest for the ANM community. Secondly, a deeper analysis of the used design patterns as in [27] could result in improving the Oriented Object structure and adding more flexibility related to the Model design. Thirdly, improving internal algorithms and data structures, without changing the interface for the final user, would lead to decreasing the computation time and improving the memory management. And finally, allowing parallel or multi-threads feature would enable us to solve cases close to industrial problems.