diff --git a/doc/manual/manual-bibliography.bib b/doc/manual/manual-bibliography.bib
index db4e9d8af..4ebe2f8ff 100644
--- a/doc/manual/manual-bibliography.bib
+++ b/doc/manual/manual-bibliography.bib
@@ -1,576 +1,582 @@
 % This file was created with JabRef 2.10b2.
 % Encoding: UTF-8
 
 
 @Article{aifantis84a,
   Title                    = {On the microstructural origin of certain inelastic models},
   Author                   = {E. C. Aifantis},
   Journal                  = {Journal of Engineering Materials and Technology},
   Year                     = {1984},
   Pages                    = {326 - 330},
   Volume                   = {106}
 }
 
 @Article{Aragon:2013d,
   Title                    = {A hierarchical detection framework for computational contact mechanics},
   Author                   = {Alejandro M. Arag{\'o}n and Jean-Fran{\c c}ois Molinari},
   Journal                  = {Computer Methods in Applied Mechanics and Engineering},
   Year                     = {2014},
   Number                   = {0},
   Pages                    = {574 - 588},
   Volume                   = {268},
 
   Bdsk-url-1               = {http://www.sciencedirect.com/science/article/pii/S0045782513002533},
   Bdsk-url-2               = {http://dx.doi.org/10.1016/j.cma.2013.10.001},
   Date-added               = {2014-05-22 09:14:48 +0000},
   Date-modified            = {2014-05-22 09:14:48 +0000},
   Doi                      = {10.1016/j.cma.2013.10.001},
   ISSN                     = {0045-7825},
   Url                      = {http://dx.doi.org/10.1016/j.cma.2013.10.001}
 }
 
 @Article{bazant86a,
   Title                    = {Mechanics of distributed cracking},
   Author                   = {Zden\v{e}k P. Ba\v{z}ant},
   Journal                  = {Applied Mechanics Reviews},
   Year                     = {1986},
   Number                   = {5},
   Pages                    = {675 - 705},
   Volume                   = {39}
 }
 
 @Article{bazant76a,
   Title                    = {Instability, Ductility, and Size Effect in Strain-Softening Concrete},
   Author                   = {Zden\v{e}k P. Ba\v{z}ant},
   Journal                  = {Journal of the Engineering Mechanics Division},
   Year                     = {1976},
   Number                   = {5},
   Pages                    = {331 - 344},
   Volume                   = {102}
 }
 
 @Article{bazant85a,
   Title                    = {Wave Propagation in a Strain-Softening Bar: Exact Solution},
   Author                   = {Zden\v{e}k P. Ba\v{z}ant and Ted Belytschko},
   Journal                  = {Journal of Engineering Mechanics},
   Year                     = {1985},
   Number                   = {3},
   Pages                    = {381 - 389},
   Volume                   = {111}
 }
 
 @Article{bazant88a,
   Title                    = {{Nonlocal Continuum Damage, Localization Instability and Convergence}},
   Author                   = {Zden\v{e}k P. Ba\v{z}ant and Gilles Pijaudier-Cabot},
   Journal                  = {Journal of Applied Mechanics},
   Year                     = {1988},
   Number                   = {2},
   Pages                    = {287-293},
   Volume                   = {55},
 
   Optdoi                   = {10.1115/1.3173674},
   Publisher                = {ASME}
 }
 
 @Article{belytschko03a,
   Title                    = {Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment},
   Author                   = {T Belytschko and H Chen and J Xu and G Zi},
   Journal                  = {International Journal for Numerical Methods in Engineering},
   Year                     = {2003},
   Pages                    = {1875-1905},
   Volume                   = {58}
 }
 
 @Article{camacho_computational_1996,
   Title                    = {Computational modelling of impact damage in brittle materials},
   Author                   = {Camacho, {G.T.} and Ortiz, M.},
   Journal                  = {International Journal of Solids and Structures},
   Year                     = {1996},
 
   Month                    = aug,
   Number                   = {20{\textendash}22},
   Pages                    = {2899--2938},
   Volume                   = {33},
 
   Urldate                  = {2012-01-17}
 }
 
 @Article{caughey1960,
   Title                    = {Classical normal modes in damped linear dynamic systems},
   Author                   = {Caughey, TK},
   Journal                  = {Journal of Applied Mechanics},
   Year                     = {1960},
   Number                   = {2},
   Pages                    = {269--271},
   Volume                   = {27},
 
   Publisher                = {American Society of Mechanical Engineers}
 }
 
 @Book{courant56a,
   Title                    = {On the partial difference equations of mathematical physics},
   Author                   = {Courant, Richard and Friedrichs, Kurt Otto and Lewy, H.},
   Publisher                = {Courant Institute of Mathematical Sciences, New York University},
   Year                     = {1956},
 
   Address                  = {New York}
 }
 
 @Book{curnier92a,
   Title                    = {M{\'e}thodes num{\'e}riques en m{\'e}canique des solides},
   Author                   = {Curnier, A.},
   Publisher                = {EPFL},
   Year                     = {1992},
 
   Bdsk-url-1               = {http://books.google.ch/books?id=-Z2zygAACAAJ},
   Url                      = {http://books.google.ch/books?id=-Z2zygAACAAJ}
 }
 
 @Article{youn10a,
   Title                    = {Studies of dynamic crack propagation and crack branching with peridynamics},
   Author                   = {Youn Doh-Ha and Florin Bobaru},
   Journal                  = {International Journal of Fracture},
   Year                     = {2010},
   Pages                    = {229-1244},
   Volume                   = {162}
 }
 
 @Book{frey2009,
   Title                    = {M{\'e}thodes des {\'e}l{\'e}ments finis: Analyse des structures et milieux continus},
   Author                   = {Frey, F. and Jirousek, J.},
   Publisher                = {PPUR},
   Year                     = {2009},
   Series                   = {Trait{\'e} de g{\'e}nie civil de l'Ecole Polytechnique F{\'e}d{\'e}rale de Lausanne},
 
   Bdsk-url-1               = {http://books.google.fr/books?id=bCBtQgAACAAJ},
   ISBN                     = {9782880748524},
   Url                      = {http://books.google.fr/books?id=bCBtQgAACAAJ}
 }
 
 @Article{gmsh,
   Title                    = {Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities},
   Author                   = {Geuzaine, Christophe and Remacle, Jean-Fran{\c c}ois},
   Journal                  = {International Journal for Numerical Methods in Engineering},
   Year                     = {2009},
   Number                   = {11},
   Pages                    = {1309--1331},
   Volume                   = {79},
 
   Bdsk-url-1               = {http://dx.doi.org/10.1002/nme.2579},
   Doi                      = {10.1002/nme.2579},
   ISSN                     = {1097-0207},
   Keywords                 = {computer-aided design, mesh generation, post-processing, finite element method, open-source software},
   Publisher                = {John Wiley \& Sons, Ltd.},
   Url                      = {http://dx.doi.org/10.1002/nme.2579}
 }
 
 @Article{giry11a,
   Title                    = {Stress-based nonlocal damage model},
   Author                   = {C. Giry and F. Dufour and J. Mazars},
   Journal                  = {International Journal of Solids and Structures},
   Year                     = {2011},
   Pages                    = {3431 - 3443},
   Volume                   = {48}
 }
 
 @Article{hughes-83a,
   Title                    = {A precis of developments in computational methods for transient analysis},
   Author                   = {Hughes, T.J.R. and Belytschko, T.},
   Journal                  = {Journal. of Applied Mechanics (ASME)},
   Year                     = {1983},
   Pages                    = {1033-1041},
   Volume                   = {50}
 }
 
 @Article{hughes83a,
   Title                    = {{A Pr{\'e}cis of Developments in Computational Methods for Transient Analysis}},
   Author                   = {{Hughes}, T.~J.~R. and {Belytschko}, T.},
   Journal                  = {Journal of Applied Mechanics},
   Year                     = {1983},
 
   Month                    = {December},
   Pages                    = {1033},
   Volume                   = {50},
 
   Optdoi                   = {10.1115/1.3167186}
 }
 
 @Article{jirasek07a,
   Title                    = {Mathematical analysis of strain localization},
   Author                   = {Milan Jir\'asek},
   Journal                  = {Revue Europeenne de Genie Civil},
   Year                     = {2007},
   Pages                    = {977 - 991},
   Volume                   = {11}
 }
 
 @Article{jirasek07b,
   Title                    = {Nonlocal damage mechanics},
   Author                   = {Milan Jir\'asek},
   Journal                  = {Revue Europeenne de Genie Civil},
   Year                     = {2007},
   Pages                    = {993 - 1021},
   Volume                   = {11}
 }
 
 @Article{jirasek98a,
   Title                    = {Nonlocal models for damage and fracture: comparison of approaches},
   Author                   = {Milan Jir\'asek},
   Journal                  = {International Journal of Solids and Structures},
   Year                     = {1998},
   Pages                    = {4133 - 4145},
   Volume                   = {35}
 }
 
 @Article{jirasek03a,
   Title                    = {Comparison of integral-type nonlocal plasticity models for strain-softening materials},
   Author                   = {Milan Jir\'asek and Simon Rolshoven},
   Journal                  = {International Journal of Engineering Science},
   Year                     = {2003},
   Pages                    = {1553-1602},
   Volume                   = {41}
 }
 
 @Article{jirasek04a,
   Title                    = {Size effect on fracture energy induced by non-locality},
   Author                   = {Milan Jir\'asek and S Rolshoven and P Grassl},
   Journal                  = {International Journal for numerical and analytical methods in geomechanics},
   Year                     = {2004},
   Pages                    = {653 - 670},
   Volume                   = {28}
 }
 
 @Article{ladeveze92a,
   Title                    = {A damage computational method for composite structures},
   Author                   = {P. Ladeveze},
   Journal                  = {Computational \& Structures},
   Year                     = {1992},
   Pages                    = {79-87},
   Volume                   = {44}
 }
 
 @Book{Laursen:2002,
   Title                    = {Computational Contact and Impact Mechanics: Fundamentals of Modeling Interfacial Phenomena in Nonlinear Finite Element Analysis},
   Author                   = {Laursen, T.A.},
   Publisher                = {Springer},
   Year                     = {2002},
   Series                   = {Engineering online library},
 
   Bdsk-url-1               = {http://books.google.ch/books?id=umzsErNuyFgC},
   Date-added               = {2014-03-21 13:24:56 +0000},
   Date-modified            = {2014-03-21 13:24:56 +0000},
   ISBN                     = {9783540429067},
   Lccn                     = {2002511027}
 }
 
 @Book{lemaitre96a,
   Title                    = {A Course on Damage Mechanics},
   Author                   = {Lemaitre, Jean},
   Publisher                = {Springer Berlin Heidelberg},
   Year                     = {1996},
 
   ISBN                     = {978-3-540-60980-3}
 }
 
 @PhdThesis{levy10a,
   Title                    = {Exploring the {P}hysics behind {D}ynamic {F}ragmentation through {P}arallel {S}imulations},
   Author                   = {Levy, Sarah},
   School                   = {ENAC},
   Year                     = {2010},
 
   Address                  = {Lausanne},
 
   Affiliation              = {EPFL},
   Doctoral                 = {EDME},
   Institute                = {IIC},
   Original-unit            = {LSMS},
   Publisher                = {EPFL},
   Unit                     = {LSMS}
 }
 
 @Article{marigo81a,
   Title                    = {{Formulation d'une loi d'endommagement d'un mat\'eriau \'elastique}},
   Author                   = {Marigo, Jean-Jacques},
   Journal                  = {{C. R. Acad. Sci., Paris, S\'er. II}},
   Year                     = {1981},
   Number                   = {19},
   Pages                    = {1309-1312},
   Volume                   = {292},
 
   ISSN                     = {0249-6305}
 }
 
 @PhdThesis{mazars84a,
   Title                    = {{Application de la m\'ecanique de l'endommagement au comportement non lin\'eaire et \`a la rupture du b\'eton de structure}},
   Author                   = {Mazars, Jacky},
   School                   = {Universit\'e Paris 6},
   Year                     = {1984}
 }
 
 @Article{needleman88a,
   Title                    = {Material rate dependence and mesh sensitivity in localization problems},
   Author                   = {A. Needleman},
   Journal                  = {Computer Methods in Applied Mechanics and Engineering},
   Year                     = {1988},
   Number                   = {1},
   Pages                    = {69 - 85},
   Volume                   = {67},
 
   ISSN                     = {0045-7825},
   Optdoi                   = {10.1016/0045-7825(88)90069-2}
 }
 
 @Article{newmark59a,
   Title                    = {{A Method of Computation for Structural Dynamics}},
   Author                   = {Nathan M. Newmark},
   Journal                  = {Journal of the Engineering Mechanics Division},
   Year                     = {1959},
 
   Month                    = {July},
   Number                   = {3},
   Pages                    = {67-94},
   Volume                   = {85},
 
   Editor                   = {ASCE}
 }
 
 @Article{nguyen2001,
   Title                    = {A cohesive model of fatigue crack growth},
   Author                   = {Nguyen, O. and Repetto, E. A. and Ortiz, M. and Radovitzky, R. A.},
   Journal                  = {International Journal of Fracture},
   Year                     = {2001},
 
   Month                    = aug,
   Number                   = {4},
   Pages                    = {351--369},
   Volume                   = {110},
 
   Doi                      = {10.1023/A:1010839522926},
   ISSN                     = {0376-9429, 1573-2673},
   Language                 = {en},
   Urldate                  = {2015-02-17}
 }
 
 @TechReport{Omohundro:1989,
   Title                    = {Five Balltree Construction Algorithms},
   Author                   = {Stephen M. Omohundro},
   Institution              = {International Computer Science Institute, University of California at Berkeley},
   Year                     = {1989},
   Number                   = {TR-89-063},
   Type                     = {Technical Report},
 
   Bdsk-url-1               = {http://ftp.icsi.berkeley.edu/ftp/pub/techreports/1989/tr-89-063.pdf},
   Date-added               = {2014-05-22 09:47:58 +0000},
   Date-modified            = {2014-05-22 09:47:58 +0000},
   Url                      = {http://ftp.icsi.berkeley.edu/ftp/pub/techreports/1989/tr-89-063.pdf}
 }
 
 @Article{ortiz1999,
   Title                    = {Finite-deformation irreversible cohesive elements for three-dimensional crack-propagation analysis},
   Author                   = {Ortiz, M. and Pandolfi, A.},
   Journal                  = {International Journal for Numerical Methods in Engineering (IJNME)},
   Year                     = {1999},
   Pages                    = {1267-1282},
   Volume                   = {44}
 }
 
 @Article{ortiz99a,
   Title                    = {Finite-deformation irreversible cohesive elements for three-dimensional crack-propagation analysis},
   Author                   = {M Ortiz and A Pandolfi},
   Journal                  = {International Journal for Numerical Methods in Engineering},
   Year                     = {1999},
   Pages                    = {1267-1282},
   Volume                   = {44}
 }
 
 @Article{pandolfi12a,
   Title                    = {An eigenerosion approach to brittle fracture},
   Author                   = {A Pandolfi and M Ortiz},
   Journal                  = {International Journal for Numerical Methods in Engineering},
   Year                     = {2012},
   Pages                    = {694-714},
   Volume                   = {92}
 }
 
 @Article{patzak01a,
   Title                    = {Parallel explicit finite element dynamics with nonlocal constitutive models},
   Author                   = {B. Patz{\'a}k and D. Rypl and Z. Bittnar},
   Journal                  = {{Computers \& Structures}},
   Year                     = {2001},
   Number                   = {26--28},
   Pages                    = {2287 - 2297},
   Volume                   = {79},
 
   ISSN                     = {0045-7949}
 }
 
 @PhdThesis{Pietrzak:1997,
   Title                    = {Continuum mechanics modelling and augmented Lagrangian formulation of large deformation frictional contact problems},
   Author                   = {Pietrzak, Grzegorz},
   School                   = {{\'E}cole {P}olytechnique {F}{\'e}d{\'e}rale de {L}ausanne},
   Year                     = {1997},
 
   Date-added               = {2014-09-16 12:23:01 +0000},
   Date-modified            = {2014-09-16 12:41:09 +0000},
   Doi                      = {10.5075/epfl-thesis-1656}
 }
 
 @Article{pijaudier87a,
   Title                    = {Nonlocal damage theory},
   Author                   = {Gilles Pijaudier-Cabot and Zden\v{e}k P. Ba\v{z}ant},
   Journal                  = {Journal of Engineering Mechanics},
   Year                     = {1987},
   Number                   = {10},
   Pages                    = {1512 - 1533},
   Volume                   = {113}
 }
 
 @Article{rice76a,
   Title                    = {The Localisation of Plastic Deformation},
   Author                   = {James R. Rice},
   Journal                  = {Theoretical and Applied Mechanics (14th International Congress on Theoretical and Applied Mechanics, Delft, 1976, ed. W.T. Koiter)},
   Year                     = {1976},
   Pages                    = {207 - 220},
   Volume                   = {1}
 }
 
 @Article{sharon96a,
   Title                    = {Microbranching instability and the dynamic fracture of brittle materials},
   Author                   = {Eran Sharon and Jay Fineberg},
   Journal                  = {Physical Review B},
   Year                     = {1996},
   Number                   = {10},
   Pages                    = {7128 - 7139},
   Volume                   = {54}
 }
 
 @Article{silling00a,
   Title                    = {Reformulation of elasticity theory for discontinuities and long-range forces},
   Author                   = {S A Silling},
   Journal                  = {Journal of the Mechanics and Physics of Solids},
   Year                     = {2000},
   Pages                    = {175-209},
   Volume                   = {48}
 }
 
 @Book{simo92,
   Title                    = {Computational Inelasticity},
   Author                   = {Simo, J.C. and Hughes, T.J.R.},
   Publisher                = {Springer},
   Year                     = {1992}
 }
 
 @Article{snozzi_cohesive_2013,
   Title                    = {A cohesive element model for mixed mode loading with frictional contact capability},
   Author                   = {Snozzi, Leonardo and Molinari, Jean-Francois},
   Journal                  = {International Journal for Numerical Methods in Engineering},
   Year                     = {2013},
 
   Month                    = feb,
   Number                   = {5},
   Pages                    = {510--526},
   Volume                   = {93}
 }
 
 @Book{Belytschko:2000,
   Title                    = {Nonlinear Finite Elements for Continua and Structures},
   Author                   = {Ted Belytschko, Wing Kam Liu, Brian Moran},
   Publisher                = {Wiley},
   Year                     = {2000}
 }
 
 @Article{devree95,
   Title                    = {Comparison of nonlocal approaches in continuum damage mechanics},
   Author                   = {J.H.P. de Vree and W.A.M. Brekelmans and M.A.J. van Gils},
   Journal                  = {Computers \&; Structures},
   Year                     = {1995},
   Number                   = {4},
   Pages                    = {581 - 588},
   Volume                   = {55},
 
   ISSN                     = {0045-7949},
   Optdoi                   = {10.1016/0045-7949(94)00501-S}
 }
 
 @Unpublished{vocialta15,
   Title                    = {3D dynamic fragmentation with parallel dynamic insertion of cohesive elements},
   Author                   = {M. Vocialta and N. Richart and J.-F. Molinari},
   Note                     = {Submitted to IJNME},
   Year                     = {2015}
 }
 
 @Unpublished{wolff14a,
   Title                    = {A non-local continuum damage approach to model dynamic crack branching},
   Author                   = {C. Wolff and N. Richart and J.-F. Molinari},
   Note                     = {Submitted to IJNME},
   Year                     = {2014}
 }
 
 @Article{Zhou_Molinari_2004,
   Title                    = {Dynamic crack propagation with cohesive elements: a methodology to address mesh dependency},
   Author                   = {F. Zhou and J. F. Molinari},
   Journal                  = {International Journal for Numerical Methods in Engineering},
   Year                     = {2004},
 
   Timestamp                = {2015.07.30}
 }
 
 @Misc{abaqus,
   Title                    = {Simulia ABAQUS FEA},
 
   Bdsk-url-1               = {http://www.3ds.com/products-services/simulia/portfolio/abaqus/},
   Key                      = {Unified FEA},
   Url                      = {\url{http://www.3ds.com/products-services/simulia/portfolio/abaqus/}}
 }
 
 @Misc{cmake,
   Title                    = {CMake - Cross Platform Make},
 
   Bdsk-url-1               = {http://www.cmake.org/},
   Url                      = {\url{http://www.cmake.org/}}
 }
 
 @Misc{diana,
   Title                    = {TNO DIANA},
 
   Bdsk-url-1               = {http://tnodiana.com/content/DIANA},
   Key                      = {FEM},
   Url                      = {\url{http://tnodiana.com/content/DIANA}}
 }
 
 @Misc{mayavi,
   Title                    = {The MayaVi Data Visualizer},
 
   Bdsk-url-1               = {http://mayavi.sourceforge.net/},
   Url                      = {\url{http://mayavi.sourceforge.net/}}
 }
 
 @Misc{mumps,
   Title                    = {MUMPS : a parallel sparse direct solver},
 
   Bdsk-url-1               = {http://graal.ens-lyon.fr/MUMPS/},
   Key                      = {sparse matrix, direct solver, parallelisme},
   Url                      = {\url{http://graal.ens-lyon.fr/MUMPS/}}
 }
 
+@Misc{numpy,
+  Title                    = {NumPy - Fundamental package for scientific computing with Python},
+
+  Bdsk-url-1               = {http://www.numpy.org/},
+  Url                      = {\url{http://www.numpy.org/}}
+}  
 @Misc{paraview,
   Title                    = {ParaView - Open Source Scientific Visualization},
 
   Bdsk-url-1               = {http://www.paraview.org/},
   Url                      = {\url{http://www.paraview.org/}}
 }
 
 @Misc{scotch,
   Title                    = {SCOTCH: Static Mapping, Graph, Mesh and Hypergraph Partitioning},
 
   Bdsk-url-1               = {http://www.labri.fr/perso/pelegrin/scotch/},
   Url                      = {\url{http://www.labri.fr/perso/pelegrin/scotch/}}
 }
 
 @Misc{visit,
   Title                    = {VisIt Visualization Tool},
 
   Bdsk-url-1               = {http://wci.llnl.gov/codes/visit/},
   Url                      = {\url{http://wci.llnl.gov/codes/visit/}}
 }
 
diff --git a/doc/manual/manual-cohesive_elements_insertion.tex b/doc/manual/manual-cohesive_elements_insertion.tex
index a13692960..20e251bab 100644
--- a/doc/manual/manual-cohesive_elements_insertion.tex
+++ b/doc/manual/manual-cohesive_elements_insertion.tex
@@ -1,81 +1,81 @@
 For cohesive material, \akantu has a pre-defined material selector to assign 
 the first cohesive material by default to the cohesive elements which is called 
 \code{DefaultMaterialCohesiveSelector} and it inherits its properties from 
 \code{DefaultMaterialSelector}. Multiple cohesive materials can be assigned 
 using mesh data information (for more details, see \ref{intrinsic_insertion}).
 
 \subsection{Insertion of Cohesive Elements}
 <<<<<<< Updated upstream
 Cohesive elements are currently compatible only with static simulation
 and dynamic simulation with an explicit time integration scheme (see
 section~\ref{ssect:smm:expl-time-integr}). They do not have to be inserted when the mesh is generated (intrinsic) but can be added during the simulation (extrinsic). At any time during the simulation, it is possible to access the following energies with the relative function:
 \begin{cpp}
   Real Ed = model.getEnergy("dissipated");
   Real Er = model.getEnergy("reversible");
   Real Ec = model.getEnergy("contact");
 \end{cpp}
 
 A new model have to be call in a very similar way that the solid mechanics model:
 \begin{cpp}
   SolidMechanicsModelCohesive model(mesh);
   model.initFull(SolidMechanicsModelCohesiveOptions(_explicit_lumped_mass, true));
 \end{cpp} 
 
 
 \subsubsection{Extrinsic approach}
 The dynamic insertion of extrinsic cohesive elements should be initialized 
 in the following way:
 \begin{cpp}
   model.updateAutomaticInsertion();
 \end{cpp} 
 During the simulation, stress has to be checked along each facet in order to insert cohesive elements where thestress criterion is reached.
 This check is performed by calling the method \code{checkCohesiveStress}, as 
 example before each step resolution:
 \begin{cpp}
   model.checkCohesiveStress();
   model.solveStep();
 \end{cpp}
 The area where stresses are checked and cohesive elements inserted can be limited 
 using the method \code{limitInsertion} during initialization. As example, to 
 limit insertion in the range $[-1.5, 1.5]$ in the $x$ direction: 
 \begin{cpp}
   model.limitInsertion(_x, -1.5, 1.5);
   model.updateAutomaticInsertion();
 \end{cpp} 
 Additional restrictions with respect to $y$ and $z$ directions can be added as well.
 
 \subsubsection{Intrinsic approach \label{intrinsic_insertion}}
 Intrinsic cohesive elements are inserted in the mesh with the method 
 \code{insertIntrinsicElements}. Similarly, the range of insertion can me limited 
 with \code{limitInsertion}. As example with a static simulation,
 \begin{cpp}
   model.limitInsertion(_x, -1.5, 1.5);
   model.insertIntrinsicElements();
 \end{cpp} 
 Mesh data information becomes vital to the insertion of cohesive elements along 
 surface with more sophisticated geometry or when multiple cohesive materials are 
 wanted. To do so, cohesive elements can be inserted along a 
 specific group of surface elements identified in a GMSH geometry file. This can 
 be achieved with the material selector (see section~\ref{sect:smm:materialselector}), in the input file specify the name of these physical groups in the corresponding cohesive materials, and call these material in the \textit{mesh parameters} section. As example, with two physical surfaces named 
 \textit{weak\_interface} and \textit{strong\_interface} defined in the GMSH 
 geometry file:
 \begin{cpp}
 ...
-  material %\emph{cohesive constitutive\_law}% [
+  material %\emph{cohesive\_constitutive\_law}% [
      name = weak_interface
      sigma_c = $value$
      ...
   ]
-  material %\emph{cohesive constitutive\_law}% [
+  material %\emph{cohesive\_constitutive\_law}% [
      name = strong_interface
      sigma_c = $value$
      ...
   ]
   mesh parameters [
      	cohesive_surfaces = weak_interface,strong_interface
   ]
 \end{cpp}
 
 In this case, there is no need to call \code{insertIntrinsicElements} anymore 
 since the insertion of cohesive elements along physical surfaces is performed 
 automatically during \code{initFull} call.    
diff --git a/doc/manual/manual-constitutive-laws-non_local.tex b/doc/manual/manual-constitutive-laws-non_local.tex
index 469f5bb9c..9cd40c0eb 100644
--- a/doc/manual/manual-constitutive-laws-non_local.tex
+++ b/doc/manual/manual-constitutive-laws-non_local.tex
@@ -1,29 +1,29 @@
 \section{Non-Local Constitutive Laws \label{sect:smm:CLNL}}\index{Material}
 
 Continuum damage modeling of quasi-brittle materials undergo significant softening after the onset of damage. This fast growth of damage causes a loss of ellipticity of partial differential equations of equilibrium. Therefore, the numerical simulation results won't be objective anymore, because the dissipated energy will depend on mesh size used in the simulation. One way to avoid this effect is the use of non-local damage formulations. In this approach a local quantity such as the strain is replaced by its non-local average, where the size of the domain, over which the quantitiy is averaged, depends on the underlying material microstructure. 
-\akantu provides non-local versions of many constitutive laws for damage. Examples are for instance the material Mazar and the material Marigo, that can be used in a non-local context. In order to use the corresponding non-local formulation the user has to define the non-local material he wishes to use in the material file:
+\akantu provides non-local versions of many constitutive laws for damage. Examples are for instance the material Mazar and the material Marigo, that can be used in a non-local context. In order to use the corresponding non-local formulation the user has to define the non-local material he wishes to use in the text input file:
 \begin{cpp}
   material %\emph{constitutive\_law\_non\_local}% [
      name = %\emph{material\_name}
      rho = $value$
      ...
   ]
 \end{cpp}
 where \emph{constitutive\_law\_non\_local} is the name of the non-local consitutive law, \textit{e.g.} \emph{marigo\_non\_local}.
 In addition to the material the non-local neighborhood, that should be used for the averaging process needs to be defined in the material file as well: 
 \begin{cpp}
   non_local %\emph{neighborhood\_name}%  %\emph{weight\_function\_type}% [
      radius = $value$
      ...
       weight_function weight_parameter [
         damage_limit = $value$
         ...
      ]
   ]
 \end{cpp}
 for the non-local averaging, \textit{e.g.} \emph{base\_wf}, followed by the properties of the non-local neighborhood, such as the radius, and the weight function parameters. It is important to notice that the non-local neighborhood must have the same name as the material to which the neighborhood belongs!
 The following two sections list the non-local constitutive laws and different type of weight functions available in \akantu.
 \subsection{Non-local constitutive laws}
 \textbf{Description to be added!!!}
 \subsection{Non-local weight functions}
  \textbf{Description to be added!!!}
\ No newline at end of file
diff --git a/doc/manual/manual-constitutive-laws.tex b/doc/manual/manual-constitutive-laws.tex
index 220cff2f7..b2041c04f 100644
--- a/doc/manual/manual-constitutive-laws.tex
+++ b/doc/manual/manual-constitutive-laws.tex
@@ -1,543 +1,521 @@
 \section{Constitutive Laws \label{sect:smm:CL}}\index{Material}
 In order to compute an element's response to deformation, one needs to
 use an appropriate constitutive relationship. The constitutive law is
 used to compute the element's stresses from the element's strains.
 
 In the finite-element discretization, the constitutive formulation is
 applied to every quadrature point of each element. When the implicit
 formulation is used, the tangent matrix has to be computed.
 
 The chosen materials for the simulation have to be specified in the
 mesh file or, as an alternative, they can be assigned using the
 \code{element\_material} vector.  For every material assigned to the
 problem one has to specify the material characteristics (constitutive
-behavior and material properties) in a text file (\eg material.dat) as
-follows:
-\begin{cpp}
-  material %\emph{constitutive\_law}% %\emph{<optional flavor>}% [
-     name = $value$
-     rho = $value$
-     ...
-  ]
-\end{cpp}
-\index{Constitutive\_laws} where \emph{constitutive\_law} is the adopted
-constitutive law, followed by the material properties listed one by line in the
-bracket (\eg \code{name} and density \code{rho}). Some constitutive laws can
-also have an \emph{optional flavor}. The file needs to be loaded in \akantu
-using the \code{initialize} method of \akantu (as shown in
-Section~\ref{sec:writing_main})
-\begin{cpp}
-  initialize("material.dat", argc, argv);
-\end{cpp}
-% or, alternatively, the \code{initFull} method.
-% \begin{cpp}
-%   model.initFull("material.dat");
-% \end{cpp}
-
+behavior and material properties) using the text input file (see \ref{sect:io:material}).\\
 In order to conveniently store values at each quadrature in a material
 point \akantu provides a special data structure, the
 \code{InternalField}. The internal fields are inheriting from the
 \code{ElementTypeMapArray}.  Furthermore, it provides several functions for
 initialization, auto-resizing and auto removal of quadrature points.
 
 Sometimes it is also desired to generate random distributions of
 internal parameters. An example might be the critical stress at which the
-material fails. To generate such a field, in the material input file,
+material fails. To generate such a field, in the text input file,
 a random quantity needs be added to the base value:
 \begin{cpp}
   sigma_c = $base$
   sigma_c = $base$ uniform [$min$, $max$]
   sigma_c = $base$ weibull [$\lambda$, $m$]
 \end{cpp}
 
 All parameters are real numbers. For the uniform distribution, minimum
 and maximum values have to be specified.
 Random parameters are defined as a $base$ value to which we add a random number
 that follows the chosen distribution.
 
 The
 \href{http://en.wikipedia.org/wiki/Uniform\_distribution\_(continuous)}{\emph{Uniform}}
 distribution is gives a random values between in $[min, max)$. The
 \href{http://en.wikipedia.org/wiki/Weibull\_distribution}{\emph{Weibull}}
 distribution is characterized by the following cumulative distribution
 function:
 \begin{equation}
   F(x) = 1- e^{-\left({x/\lambda}\right)^m}
 \end{equation}
 which depends on  $m$ and $\lambda$, which are the shape parameter and the scale
 parameter. These random distributions are different each time the code
 is executed. In order to obtain always the same one, it possible to
 manually set the \emph{seed} that is the number from which these
 pseudo-random distributions are created. This can be done by adding
 the following line to the input file \emph{outside} the material
 parameters environments:
 \begin{cpp}
   seed = 1.0
 \end{cpp}
 where the value 1 can be substituted with any number. Currently
 \akantu is can reproduce always the same distribution when the seed is
 specified \emph{only} in serial.
 
 The following sections describe the constitutive models implemented in
 \akantu. In Appendix~\ref{app:material-parameters} a summary of the
 parameters for all materials of \akantu is provided.
 
 
 \subsection{Elasticity}\index{Material!Elastic}
 
 The elastic law is a commonly used constitutive relationship that can be used
 for a wide range of engineering materials (\eg metals, concrete, rock, wood,
 glass, rubber, etc.) provided that the strains remain small (\ie small
 deformation and stress lower than yield strength).
 
 The elastic laws are often expressed as $\mat{\sigma} =
 \mat{C}:\mat{\varepsilon}$ with where $\mat{\sigma}$ is the Cauchy stress tensor,
 $\mat{\varepsilon}$ represents the infinitesimal strain tensor and $\mat{C}$ is the
 elastic modulus tensor.
 
 \subsubsection{Linear isotropic\matlabel{ssect:smm:linear-elastic-isotropic}}
 
 The linear isotropic elastic behavior is described by Hooke's law, which states
 that the stress is linearly proportional to the applied strain (material behaves
 like an ideal spring), as illustrated in Figure~\ref{fig:smm:cl:elastic}.
 \begin{figure}[!htb]
   \begin{center}
 
     \subfloat[]{
       \begin{tikzpicture}
 	\draw[thick,latex-latex] (0,5) node[left] {$\sigma$} |- (5,0) node (x) [right, below] {$\varepsilon$};
 	\draw[thin] (1.5,1.5) -- (2.5,1.5) -- (2.5,2.5) node [midway, right] {E};
 	\draw[very thick,color=red] (0,0) -- (4,4);
 	\draw[very thick,latex-latex,color=red] (1,1) -- (3,3);
       \end{tikzpicture}
       \label{fig:smm:cl:elastic:stress_strain} }
     \hspace{0.05\textwidth} \subfloat[]{
       \raisebox{0.125\textwidth}{\includegraphics[width=0.25\textwidth,keepaspectratio=true]{figures/hooke_law.pdf}}
       \label{fig:smm:cl:elastic:hooke} }
     \caption{(a) Stress-strain curve for elastic material and (b)
       schematic representation of Hooke's law, denoted as a spring.}
     \label{fig:smm:cl:elastic}
   \end{center}
 \end{figure}
 The equation that relates the strains to the
 displacements is: % First the strain is computed (at every gauss
 point) from the displacements as follows:
 \begin{equation}
   \label{eqn:smm:strain_inf}
   \mat{\varepsilon} =
   \frac{1}{2} \left[ \nabla_0 \vec{u}+\nabla_0 \vec{u}^T \right]
 \end{equation}
 where $\mat{\varepsilon}$ represents the infinitesimal strain tensor,
 $\nabla_{0}\vec{u}$ the displacement gradient
 tensor according to the initial configuration. The constitutive equation
 for isotropic homogeneous media can be expressed as:
 \begin{equation}
   \label{eqn:smm:material:constitutive_elastic}
   \mat{\sigma } =\lambda\mathrm{tr}(\mat{\varepsilon})\mat{I}+2 \mu\mat{\varepsilon}
 \end{equation}
 where $\mat{\sigma}$ is the Cauchy stress tensor
 ($\lambda$ and $\mu$ are the the first and second Lame's
 coefficients).
 
 In Voigt notation this correspond to
 \begin{align}
   \left[\begin{array}{c}
       \sigma_{11}\\
       \sigma_{22}\\
       \sigma_{33}\\
       \sigma_{23}\\
       \sigma_{13}\\
       \sigma_{12}\\
     \end{array}\right]
   &= \frac{E}{(1+\nu)(1-2\nu)}\left[
     \begin{array}{cccccc}
       1-\nu & \nu   & \nu   & 0 & 0 & 0\\
       \nu   & 1-\nu & \nu   & 0 & 0 & 0\\
       \nu   & \nu   & 1-\nu & 0 & 0 & 0\\
       0     &  0    &  0    & \frac{1-2\nu}{2} & 0 & 0 \\
       0     &  0    &  0    & 0 & \frac{1-2\nu}{2} & 0 \\
       0     &  0    &  0    & 0 & 0 & \frac{1-2\nu}{2} \\
     \end{array}\right]
   \left[\begin{array}{c}
       \varepsilon_{11}\\
       \varepsilon_{22}\\
       \varepsilon_{33}\\
       2\varepsilon_{23}\\
       2\varepsilon_{13}\\
       2\varepsilon_{12}\\
     \end{array}\right]
 \end{align}
 
 \subsubsection{Linear anisotropic\matlabel{ssect:smm:linear-elastic-anisotropic}}
 This formulation is not sufficient to represent all elastic material
 behavior. Some materials have characteristic orientation that have to be taken
 into account. To represent this anisotropy a more general stress-strain law has
 to be used. For this we define the elastic modulus tensor as follow:
 
 \begin{align}
   \left[\begin{array}{c}
       \sigma_{11}\\
       \sigma_{22}\\
       \sigma_{33}\\
       \sigma_{23}\\
       \sigma_{13}\\
       \sigma_{12}\\
     \end{array}\right]
   &= \left[
     \begin{array}{cccccc}
       c_{11} & c_{12} & c_{13} & c_{14} & c_{15} & c_{16}\\
       c_{21} & c_{22} & c_{23} & c_{24} & c_{25} & c_{26}\\
       c_{31} & c_{32} & c_{33} & c_{34} & c_{35} & c_{36}\\
       c_{41} & c_{42} & c_{43} & c_{44} & c_{45} & c_{46}\\
       c_{51} & c_{52} & c_{53} & c_{54} & c_{55} & c_{56}\\
       c_{61} & c_{62} & c_{63} & c_{64} & c_{65} & c_{66}\\
     \end{array}\right]
   \left[\begin{array}{c}
       \varepsilon_{11}\\
       \varepsilon_{22}\\
       \varepsilon_{33}\\
       2\varepsilon_{23}\\
       2\varepsilon_{13}\\
       2\varepsilon_{12}\\
     \end{array}\right]
 \end{align}
 
 \begin{figure}[h]
   \centering
   \begin{tikzpicture}
     \draw[thick,latex-latex] (90:3) node[left] {$\vec{e_2}$} |- (0:3) node [right, below] {$\vec{e_1}$};
     \draw[ultra thick,latex-latex] (150:3) node[left] {$\vec{n_2}$} -- (0,0) -- (20:3) node [right] {$\vec{n_1}$};
   \end{tikzpicture}
   \caption{Material basis}
 \end{figure}
 
 To simplify the writing of input files the \mat{C} tensor is expressed in the
 material basis. And this basis as to be given too. This basis $\Omega_{\st{mat}}
 = \{\vec{n_1}, \vec{n_2}, \vec{n_3}\}$ is used to define the rotation $R_{ij} =
 \vec{n_j} . \vec{e_i}$. And $\mat{C}$ can be rotated in the global basis $\Omega
 = \{\vec{e_1}, \vec{e_2}, \vec{e_3}\}$ as follow:
 
 
 \begin{align}
 \mat{C}_{\Omega} &= \mat{R}_1 \mat{C}_{\Omega_{\st{mat}}} \mat{R}_2\\
 \mat{R}_1  &= \left[
   \begin{array}{cccccc}
     R_{11} R_{11} & R_{12} R_{12} & R_{13} R_{13} & R_{12} R_{13} & R_{11} R_{13} & R_{11} R_{12}\\
     R_{21} R_{21} & R_{22} R_{22} & R_{23} R_{23} & R_{22} R_{23} & R_{21} R_{23} & R_{21} R_{22}\\
     R_{31} R_{31} & R_{32} R_{32} & R_{33} R_{33} & R_{32} R_{33} & R_{31} R_{33} & R_{31} R_{32}\\
     R_{21} R_{31} & R_{22} R_{32} & R_{23} R_{33} & R_{22} R_{33} & R_{21} R_{33} & R_{21} R_{32}\\
     R_{11} R_{31} & R_{12} R_{32} & R_{13} R_{33} & R_{12} R_{33} & R_{11} R_{33} & R_{11} R_{32}\\
     R_{11} R_{21} & R_{12} R_{22} & R_{13} R_{23} & R_{12} R_{23} & R_{11} R_{23} & R_{11} R_{22}\\
   \end{array}\right]\\
 \mat{R}_2  &= \left[
   \begin{array}{cccccc}
     R_{11} R_{11} & R_{21} R_{21} & R_{31} R_{31} & R_{21} R_{31} & R_{11} R_{31} & R_{11} R_{21}\\
     R_{12} R_{12} & R_{22} R_{22} & R_{32} R_{32} & R_{22} R_{32} & R_{12} R_{32} & R_{12} R_{22}\\
     R_{13} R_{13} & R_{23} R_{23} & R_{33} R_{33} & R_{23} R_{33} & R_{13} R_{33} & R_{13} R_{23}\\
     R_{12} R_{13} & R_{22} R_{23} & R_{32} R_{33} & R_{22} R_{33} & R_{12} R_{33} & R_{12} R_{23}\\
     R_{11} R_{13} & R_{21} R_{23} & R_{31} R_{33} & R_{21} R_{33} & R_{11} R_{33} & R_{11} R_{23}\\
     R_{11} R_{12} & R_{21} R_{22} & R_{31} R_{32} & R_{21} R_{32} & R_{11} R_{32} & R_{11} R_{22}\\
   \end{array}\right]\\
 \end{align}
 
 \subsubsection{Linear orthotropic\matlabel{ssect:smm:linear-elastic-orthotropic}}
 
 A particular case of anisotropy is when the material basis is orthogonal in which case the elastic modulus tensor can be simplified and rewritten in terms of 9 independents material parameters.
 
 \begin{align}
   \left[\begin{array}{c}
       \sigma_{11}\\
       \sigma_{22}\\
       \sigma_{33}\\
       \sigma_{23}\\
       \sigma_{13}\\
       \sigma_{12}\\
     \end{array}\right]
   &= \left[
     \begin{array}{cccccc}
       c_{11} & c_{12} & c_{13} &   0   &   0   &   0  \\
             & c_{22} & c_{23} &   0   &   0   &   0  \\
             &       & c_{33} &   0   &   0   &   0  \\
             &       &       & c_{44} &   0   &   0  \\
             &  \multicolumn{2}{l}{\text{sym.}}       &       & c_{55} &   0  \\
             &       &       &       &       & c_{66}\\
     \end{array}\right]
   \left[\begin{array}{c}
       \varepsilon_{11}\\
       \varepsilon_{22}\\
       \varepsilon_{33}\\
       2\varepsilon_{23}\\
       2\varepsilon_{13}\\
       2\varepsilon_{12}\\
     \end{array}\right]
 \end{align}
 
 \begin{align}
   c_{11} &= E_1 (1 - \nu_{23}\nu_{32})\Gamma \qquad c_{22} = E_2 (1 - \nu_{13}\nu_{31})\Gamma \qquad c_{33} = E_3 (1 - \nu_{12}\nu_{21})\Gamma\\
   c_{12} &= E_1 (\nu_{21} - \nu_{31}\nu_{23})\Gamma = E_2 (\nu_{12} - \nu_{32}\nu_{13})\Gamma\\
   c_{13} &= E_1 (\nu_{31} - \nu_{21}\nu_{32})\Gamma = E_2 (\nu_{13} - \nu_{21}\nu_{23})\Gamma\\
   c_{23} &= E_2 (\nu_{32} - \nu_{12}\nu_{31})\Gamma = E_3 (\nu_{23} - \nu_{21}\nu_{13})\Gamma\\
   c_{44} &= \mu_{23} \qquad  c_{55} = \mu_{13} \qquad  c_{66} = \mu_{12} \\
   \Gamma &= \frac{1}{1 - \nu_{12} \nu_{21} - \nu_{13} \nu_{31} - \nu_{32} \nu_{23} - 2 \nu_{21} \nu_{32} \nu_{13}}
 \end{align}
 
 The Poisson ratios follow the rule $\nu_{ij} = \nu_{ji} E_i / E_j$.
 
 \subsection{Neo-Hookean\matlabel{ssect:smm:cl:neohookean}}\index{Material!Neohookean}
 The hyperelastic Neo-Hookean constitutive law results from an
 extension of the linear elastic relationship (Hooke's Law) for large
 deformation. Thus, the model predicts nonlinear stress-strain behavior
 for bodies undergoing large deformations.
 
 \begin{figure}[!htb]
   \begin{center}
     \includegraphics[width=0.4\textwidth,keepaspectratio=true]{figures/stress_strain_neo.pdf}
     \caption{Neo-hookean Stress-strain curve.}
     \label{fig:smm:cl:neo_hookean}
   \end{center}
 \end{figure}
 
 As illustrated in Figure~\ref{fig:smm:cl:neo_hookean}, the behavior is initially
 linear and the mechanical behavior is very close to the corresponding linear
 elastic material. This constitutive relationship, which accounts for compressibility,
 is a modified version of the one proposed by Ronald Rivlin \cite{Belytschko:2000}.
 
 The strain energy stored in the material is given by:
 \begin{equation}\label{eqn:smm:constitutive:neohookean_potential}
   \Psi(\mat{C}) = \frac{1}{2}\lambda_0\left(\ln J\right)^2-\mu_0\ln J+\frac{1}{2}
   \mu_0\left(\mathrm{tr}(\mat{C})-3\right)
 \end{equation}
 \noindent where $\lambda_0$ and $\mu_0$ are, respectively, Lam\'e's first parameter
 and the shear modulus at the initial configuration. $J$ is the jacobian of the deformation
 gradient ($\mat{F}=\nabla_{\!\!\vec{X}}\vec{x}$): $J=\text{det}(\mat{F})$. Finally $\mat{C}$ is the right Cauchy-Green
 deformation tensor.
 
 Since this kind of material is used for large deformation problems, a
 finite deformation framework should be used. Therefore, the Cauchy
 stress ($\mat{\sigma}$) should be computed through the second
 Piola-Kirchhoff stress tensor $\mat{S}$:
 
 \begin{equation}
   \mat{\sigma } = \frac{1}{J}\mat{F}\mat{S}\mat{F}^T
 \end{equation}
 
 Finally the second Piola-Kirchhoff stress tensor is given by:
 
 \begin{equation}
   \mat{S}  = 2\frac{\partial\Psi}{\partial\mat{C}} = \lambda_0\ln J
   \mat{C}^{-1}+\mu_0\left(\mat{I}-\mat{C}^{-1}\right)
 \end{equation}
 
 The parameters to indicate in the material file are the same
 as those for the elastic case: \code{E} (Young's modulus), \code{nu} (Poisson's
 ratio).
 
 
 \subsection{Visco-Elasticity\matlabel{ssect:smm:cl:sls}}
 % Standard Solid rheological model, see [] J.C. Simo, T.J.R. Hughes,
 % "Computational Inelasticity", Springer (1998), see Sections 10.2 and 10.3
 Visco-elasticity is characterized by strain rate dependent
 behavior. Moreover, when such a material undergoes a deformation it
 dissipates energy. This dissipation results in a hysteresis loop in
 the stress-strain curve at every loading cycle (see
 Figure~\ref{fig:smm:cl:visco-elastic:hyst}). In principle, it can be
 applied to many materials, since all materials exhibit a visco-elastic
 behavior if subjected to particular conditions (such as high
 temperatures).
 \begin{figure}[!htb]
   \begin{center}
 
     \subfloat[]{
       \includegraphics[width=0.4\textwidth,keepaspectratio=true]{figures/stress_strain_visco.pdf}
       \label{fig:smm:cl:visco-elastic:hyst}
     }
     \hspace{0.05\textwidth}
     \subfloat[]{
       \raisebox{0.025\textwidth}{\includegraphics[width=0.3\textwidth,keepaspectratio=true]{figures/visco_elastic_law.pdf}}
       \label{fig:smm:cl:visco-elastic:model}
     }
     \caption{(a) Characteristic stress-strain behavior of a visco-elastic material with hysteresis loop and (b) schematic representation of the standard rheological linear solid visco-elastic model.}
     \label{fig:smm:cl:visco-elastic}
   \end{center}
 \end{figure}
 The standard rheological linear solid model (see Sections 10.2 and 10.3
 of~\cite{simo92}) has been implemented in \akantu. This model results from the
 combination of a spring mounted in parallel with a spring and a dashpot
 connected in series, as illustrated in
 Figure~\ref{fig:smm:cl:visco-elastic:model}. The advantage of this model is that
 it allows to account for creep or stress relaxation. The equation that relates
 the stress to the strain is (in 1D):
 \begin{equation}
   \frac{d\varepsilon(t)}{dt} = \left ( E + E_V \right ) ^ {-1} \cdot \left [ \frac{d\sigma(t)}{dt} + \frac{E_V}{\eta}\sigma(t) - \frac{EE_V}{\eta}\varepsilon(t) \right ]
 \end{equation}
 where $\eta$ is the viscosity. The equilibrium condition is unique and
 is attained in the limit, as $t \to \infty $. At this stage, the
 response is elastic and depends on the Young's modulus $E$.  The
 mandatory parameters for the material file are the following:
 \code{rho} (density), \code{E} (Young's modulus), \code{nu} (Poisson's
 ratio), \code{Plane\_Stress} (if set to zero plane strain, otherwise
 plane stress), \code{eta} (dashpot viscosity) and \code{Ev} (stiffness
 of the viscous element).
 
 Note that the current standard linear solid model is applied only on the deviatoric part of the strain tensor. The spheric part of the strain tensor affects the stress tensor like an linear elastic material.
 
 \subsection{Small-Deformation Plasticity\matlabel{ssect:smm:cl:plastic}}\index{Material!Small-deformation Plasticity}
 The small-deformation plasticity is a simple plasticity material
 formulation which accounts for the additive decomposition of strain
 into elastic and plastic strain components. This formulation is
 applicable to infinitesimal deformation where the additive
 decomposition of the strain is a valid approximation. In this
 formulation, plastic strain is a shearing process where hydrostatic
 stress has no contribution to plasticity and consequently plasticity
 does not lead to volume change. Figure~\ref{fig:smm:cl:Lin-strain-hard}
 shows the linear strain hardening elasto-plastic behavior according to
 the additive decomposition of strain into the elastic and plastic
 parts in infinitesimal deformation as
 \begin{align}
   \mat{\varepsilon} &= \mat{\varepsilon}^e +\mat{\varepsilon}^p\\
   {\mat{\sigma}} &= 2G(\mat{\varepsilon}^e) + \lambda  \mathrm{tr}(\mat{\varepsilon}^e)\mat{I}
 \end{align}
 
 \begin{figure}[htp]
   \centering
   {\includegraphics[scale=0.4, clip]{figures/isotropic_hardening_plasticity.pdf}}
   \caption{
     Stress-strain curve for the small-deformation plasticity with linear isotropic hardening.
   }
   \label{fig:smm:cl:Lin-strain-hard}
 \end{figure}
 
 \noindent In this class, the von Mises yield criterion is used. In the von Mises yield criterion, the yield is independent of the hydrostatic stress. Other yielding criteria such as Tresca and Gurson can be easily implemented in this class as well.
 
 In the von Mises yield criterion, the hydrostatic stresses have no effect on the plasticity and consequently the yielding occurs when a critical elastic shear energy is achieved.
 \begin{equation} \label{eqn:smm:constitutive:von Mises}
   f = \sigma_{\st{eff}} - \sigma_y = \left(\frac{3}{2} {\mat{\sigma}}^{\st{tr}} : {\mat{\sigma}}^{\st{tr}}\right)^\frac{1}{2}-\sigma_y (\mat{\varepsilon}^p)
 \end{equation}
 \begin{equation} \label{eqn:smm:constitutive:yielding}
   f < 0 \quad \textrm{Elastic deformation,} \qquad f = 0 \quad  \textrm{Plastic deformation}
 \end{equation}
 where $\sigma_y$ is the yield strength of the material which can be function of plastic strain in case of hardening type of materials and ${\mat{\sigma}}^{\st{tr}}$ is the deviatoric part of stress given by
 \begin{equation} \label{eqn:smm:constitutive:deviatoric stress}
   {\mat{\sigma}}^{\st{tr}}=\mat{\sigma} - \frac{1}{3} \mathrm{tr}(\mat{\sigma}) \mat {I}
 \end{equation}
 
 After yielding $(f = 0)$, the normality hypothesis of plasticity determines the direction of plastic flow which is normal to the tangent to the yielding surface at the load point. Then, the tensorial form of the plastic constitutive equation using the von Mises yielding criterion (see equation 4.34) may be written as
 \begin{equation} \label{eqn:smm:constitutive:plastic contitutive equation}
   \Delta {\mat{\varepsilon}}^p = \Delta p \frac {\partial{f}}{\partial{\mat \sigma}}=\frac{3}{2} \Delta p \frac{{\mat{\sigma}}^{\st{tr}}}{\sigma_{\st{eff}}}
 \end{equation}
 
 In these expressions, the direction of the plastic strain increment (or equivalently, plastic strain rate) is given by $\frac{{\mat{\sigma}}^{\st{tr}}}{\sigma_{\st{eff}}}$ while the magnitude is defined by the plastic multiplier $\Delta p$. This can be obtained using the \emph{consistency condition} which impose the requirement for the load point to remain on the yielding surface in the plastic regime.
 
 Here, we summarize the implementation procedures for the
 small-deformation plasticity with linear isotropic hardening:
 \begin{enumerate}
 \item Compute the trial stress:
   \begin{equation}
     {\mat{\sigma}}^{\st{tr}} = {\mat{\sigma}}_t + 2G\Delta \mat{\varepsilon} + \lambda \mathrm{tr}(\Delta \mat{\varepsilon})\mat{I}
   \end{equation}
 \item Check the Yielding criteria:
   \begin{equation}
     f = (\frac{3}{2} {\mat{\sigma}}^{\st{tr}} : {\mat{\sigma}}^{\st{tr}})^{1/2}-\sigma_y (\mat{\varepsilon}^p)
   \end{equation}
 \item Compute the Plastic multiplier:
   \begin{align}
     d \Delta p &= \frac{\sigma^{tr}_{eff} - 3G \Delta P^{(k)}- \sigma_y^{(k)}}{3G + h}\\
     \Delta p^{(k+1)} &= \Delta p^{(k)}+ d\Delta p\\
     \sigma_y^{(k+1)} &= (\sigma_y)_t+ h\Delta p
   \end{align}
 \item Compute the plastic strain increment:
   \begin{equation}
     \Delta {\mat{\varepsilon}}^p = \frac{3}{2} \Delta p \frac{{\mat{\sigma}}^{\st{tr}}}{\sigma_{\st{eff}}}
   \end{equation}
 \item Compute the stress increment:
   \begin{equation}
     {\Delta \mat{\sigma}} = 2G(\Delta \mat{\varepsilon}-\Delta \mat{\varepsilon}^p) + \lambda  \mathrm{tr}(\Delta \mat{\varepsilon}-\Delta \mat{\varepsilon}^p)\mat{I}
   \end{equation}
 \item Update the variables:
   \begin{align}
     {\mat{\varepsilon^p}} &= {\mat{\varepsilon}}^p_t+{\Delta {\mat{\varepsilon}}^p}\\
     {\mat{\sigma}} &= {\mat{\sigma}}_t+{\Delta \mat{\sigma}}
   \end{align}
 \end{enumerate}
 
 We use an implicit integration technique called \emph{the radial
   return method} to obtain the plastic multiplier. This method has the
 advantage of being unconditionally stable, however, the accuracy
 remains dependent on the step size. The plastic parameters to indicate
 in the material file are: \code{$\sigma_y$} (Yield stress) and
 \code{h} (Hardening modulus). In addition, the elastic parameters need
 to be defined as previously mentioned: \code{E} (Young's modulus),
 \code{nu} (Poisson's ratio).
 
 \subsection{Damage}
 
 In the  simplified case of a  linear elastic and brittle  material, isotropic
 damage can be represented by a scalar variable $d$, which varies from $0$ to $1$
 for  no  damage  to  fully  broken  material  respectively.  The  stress-strain
 relationship then becomes:
 \begin{equation*}
   \mat{\sigma} = (1-d)\, \mat{C}:\mat{\varepsilon}
 \end{equation*}
 
 where  $\mat{\sigma}$,  $\mat{\varepsilon}$ are  the  Cauchy  stress and  strain
 tensors, and $\mat{C}$ is the elastic stiffness tensor. This formulation relies
 on the definition of an evolution law for the damage variable. In \akantu, many
 possibilities exist and they are listed below.
 
 \subsubsection{Marigo\matlabel{ssect:smm:cl:damage-marigo}}
 
 This damage evolution law is energy based as defined by Marigo \cite{marigo81a,
   lemaitre96a}. It is an isotropic damage law.
 \begin{align}
   Y &= \frac{1}{2}\mat{\varepsilon}:\mat{C}:\mat{\varepsilon}\\
   F &= Y - Y_d - S d\\
   d &= \left\{
     \begin{array}{l l}
       \mathrm{min}\left(\frac{Y-Y_d}{S},\;1\right) & \mathrm{if}\; F > 0\\
       \mathrm{unchanged} & \mathrm{otherwise}
     \end{array}
   \right.
 \end{align}
 In this formulation, $Y$ is the strain energy release rate, $Y_d$ the
 rupture criterion and $S$ the damage energy.  The non-local version of
 this damage evolution law is constructed by averaging the energy $Y$.
 
 \subsubsection{Mazars\matlabel{ssect:smm:cl:damage-mazars}}
 
 This law introduced by Mazars \cite{mazars84a} is a behavioral model to
 represent damage evolution in concrete. This model does not rely on the computation of the tangent stiffness, the damage is directly evaluated from the strain.
 
 The governing variable in this damage
 law is the equivalent strain $\varepsilon_{\st{eq}} =
 \sqrt{<\mat{\varepsilon}>_+:<\mat{\varepsilon}>_+}$, with $<.>_+$ the positive
 part of the tensor. This part is defined in the principal coordinates (I, II, III) as $\varepsilon_{\st{eq}} =
 \sqrt{<\mat{\varepsilon_I}>_+^2 + <\mat{\varepsilon_{II}}>_+^2 + <\mat{\varepsilon_{III}}>_+^2}$.
 The damage is defined as:
 \begin{align}
   D &= \alpha_t^\beta D_t + (1-\alpha_t)^\beta D_c\\
   D_t &= 1 - \frac{\kappa_0 (1- A_t)}{\varepsilon_{\st{eq}}} - A_t \exp^{-B_t(\varepsilon_{\st{eq}}-\kappa_0)}\\
   D_c &= 1 - \frac{\kappa_0 (1- A_c)}{\varepsilon_{\st{eq}}} - A_c
   \exp^{-B_c(\varepsilon_{\st{eq}}-\kappa_0)}\\
   \alpha_t &= \frac{\sum_{i=1}^3<\varepsilon_i>_+\varepsilon_{\st{nd}\;i}}{\varepsilon_{\st{eq}}^2}
 \end{align}
 With $\kappa_0$ the damage threshold, $A_t$ and $B_t$ the damage parameter in
 traction, $A_c$ and $B_c$ the damage parameter in compression, $\beta$ is the
 shear parameter. $\alpha_t$ is the coupling parameter between traction and
 compression, the $\varepsilon_i$ are the eigenstrain and the
 $\varepsilon_{\st{nd}\;i}$ are the eigenvalues of the strain if the material
 were undamaged.
 
 The coefficients $A$ and $B$ are the post-peak asymptotic
 value and the decay shape parameters.
 
 \IfFileExists{manual-constitutive-laws-non_local.tex}{\input{manual-constitutive-laws-non_local.tex}}{}
 
 \IfFileExists{manual-extra_materials.tex}{\input{manual-extra_materials}}{}
 
 \IfFileExists{manual-cohesive_laws.tex}{\input{manual-cohesive_laws}}{}
 
 
 %%% Local Variables:
 %%% mode: latex
 %%% TeX-master: "manual"
 %%% End:
diff --git a/doc/manual/manual-gettingstarted.tex b/doc/manual/manual-gettingstarted.tex
index a64ff3856..c0225df9a 100644
--- a/doc/manual/manual-gettingstarted.tex
+++ b/doc/manual/manual-gettingstarted.tex
@@ -1,434 +1,434 @@
 \chapter{Getting Started}
 \section{Downloading the Code}
 
 The \akantu source code can be requested using the form accessible at the URL
 \url{http://lsms.epfl.ch/akantu}.  There, you will be asked to accept the LGPL
 license terms.
 
 \section{Compiling \akantu}
 
 \akantu is a \code{cmake} project, so to configure it, you can either
 follow the usual way:
 \begin{command}
   > cd akantu
   > mkdir build
   > cd build
   > ccmake ..
   [ Set the options that you need ]
   > make
   > make install
 \end{command}
 
 \noindent Or, use the \code{Makefile} we added for your convenience to
 handle the \code{cmake} configuration
 
 \begin{command}
   > cd akantu
   > make config
   > make
   > make install
 \end{command}
 
 \noindent All the \akantu options are documented in Appendix
 \ref{app:package-dependencies}.
 
 
 \section{Writing a \texttt{main} Function\label{sect:common:main}}
 \label{sec:writing_main}
 
 First of all, \akantu needs to be initialized.  The memory management
 included in the core library handles the correct allocation and
 de-allocation of vectors, structures and/or objects. Moreover, in
 parallel computations, the initialization procedure performs the
 communication setup. This is achieved by a pair of functions
 (\code{initialize} and \code{finalize}) that are used as follows:
 \begin{cpp}
 #include "aka_common.hh"
 #include "..."
 
 using namespace akantu;
 
 int main(int argc, char *argv[]) {
-  initialize("material.dat", argc, argv);
+  initialize("input_file.dat", argc, argv);
 
   // your code
   ...
 
   finalize();
 }
 \end{cpp}
-The \code{initialize} function takes the material file and the program
-parameters which can be interpreted by \akantu in due form. Obviously
-it is necessary to include all files needed in main. In this manual
+The \code{initialize} function takes the text inpute file and the program
+parameters which can be parsed by \akantu in due form (see \ref{sect:parser}). 
+Obviously it is necessary to include all files needed in main. In this manual
 all provided code implies the usage of \code{akantu} as namespace.
 
 \section{Creating and Loading a Mesh\label{sect:common:mesh}}
 
 In its current state, \akantu supports three types of meshes: Gmsh~\cite{gmsh},
 Abaqus~\cite{abaqus} and Diana~\cite{diana}. Once a \code{Mesh} object is
 created with a given spatial dimension, it can be filled by reading a mesh input file.
 The method \code{read} of the class \code{Mesh} infers the mesh type from the file extension. If a non-standard file extension is used, the mesh type has to be specified.
 
 \begin{cpp}
 UInt spatial_dimension = 2;
 Mesh mesh(spatial_dimension);
 
 // Reading Gmsh files
 mesh.read("my_gmsh_mesh.msh");
 mesh.read("my_gmsh_mesh", _miot_gmsh);
 
 
 // Reading Abaqus files
 mesh.read("my_abaqus_mesh.inp");
 mesh.read("my_abaqus_mesh", _miot_abaqus);
 
 // Reading Diana files
 mesh.read("my_diana_mesh.dat");
 mesh.read("my_diana_mesh", _miot_diana);
 \end{cpp}
 
 The Gmsh reader adds the geometrical and physical tags as mesh data. The physical
 values are stored as a \code{UInt} data called \code{tag\_0}, if a string
 name is provided it is stored as a \code{std::string} data named
 \code{physical\_names}. The geometrical tag is stored as a \code{UInt} data named
 \code{tag\_1}.
 
 
 The Abaqus reader stores the \code{ELSET} in ElementGroups and the \code{NSET}
 in NodeGroups. The material assignment can be retrieved from the
 \code{std::string} mesh data named \code{abaqus\_material}.
 
 % \akantu supports meshes generated with Gmsh~\cite{gmsh}, a free
 % software available at \url{http://geuz.org/gmsh/} where a detailed
 % documentation can be found. Consequently, this manual will not provide
 % Gmsh usage directions. Gmsh outputs meshes in \code{.msh} format that can be read
 % by \akantu. In order to import a mesh, it is necessary to create
 % a \code{Mesh} object through the following function calls:
 % \begin{cpp}
 %   UInt spatial_dimension = 2;
 %   Mesh mesh(spatial_dimension);
 % \end{cpp}
 % The only parameter that has to be specified by the user is the spatial
 % dimension of the problem. Now it is possible to read a \code{.msh} file with
 % a \code{MeshIOMSH} object that takes care of loading a mesh to memory.
 % This step is carried out by:
 % \begin{cpp}
 %   mesh.read("square.msh");
 % \end{cpp}
 % where the \code{MeshIOMSH} object is first created before being
 % used to read the \code{.msh} file. The mesh file name as well as the \code{Mesh}
 % object must be specified by the user.
 
 % The \code{MeshIOMSH} object can also write mesh files. This feature
 % is useful to save a mesh that has been modified during a
 % simulation. The \code{write} method takes care of it:
 % \begin{cpp}
 %   mesh_io.write("square_modified.msh", mesh);
 % \end{cpp}
 % which works exactly like the \code{read} method.
 
 % \akantu supports also meshes generated by
 % DIANA (\url{http://tnodiana.com}), but only in reading mode. A similar
 % procedure applies where the only
 % difference is that the \code{MeshIODiana} object should be used
 % instead of the \code{MeshIOMSH} one. Additional mesh readers can be
 % introduced into \akantu by coding new \code{MeshIO} classes.
 
 \section{Using \texttt{Arrays}}
 
 Data in \akantu can be stored in data containers implemented by
 the \code{Array} class. In its most basic usage, the \code{Array} class
 implemented in \akantu is similar to the \code{vector} class of
 the Standard Template Library (STL) for C++. A simple \code{Array}
 containing a sequence of \code{nb\_element} values (of a given type) can be generated with:
 \begin{cpp}
   Array<type> example_array(nb_element);
 \end{cpp}
 where \code{type} usually is \code{Real}, \code{Int}, \code{UInt} or
 \code{bool}. Each value is associated to an index, so that data can be
 accessed by typing:
 
 \begin{cpp}
   type & val = example_array(index)
 \end{cpp}
 
 \code{Arrays} can also contain tuples of values for each index. In
 that case, the number of components per tuple must be specified at the
 \code{Array} creation.  For example, if we want to create an
 \code{Array} to store the coordinates (sequences of three values) of
 ten nodes, the appropriate code is the following:
 \begin{cpp}
   UInt nb_nodes = 10;
   UInt spatial_dimension = 3;
 
   Array<Real> position(nb_nodes, spatial_dimension);
 \end{cpp}
 In this case the $x$ position of the eighth node number will be given by
 \code{position(7, 0)} (in C++, numbering starts at 0 and not
 1). If the number of components for the sequences is not specified, the
 default value of 1 is used.
 
 It is very common in \akantu to loop over arrays to perform a specific
 treatment. This ranges from geometric calculation on nodal quantities
 to tensor algebra (in constitutive laws for example).
 The \code{Array} object has the possibility to request iterators
 in order to make the writing of loops easier and enhance readability.
 For instance, a loop over the nodal coordinates can be performed like:
 \begin{cpp}
   //accessing the nodal coordinates Array (spatial_dimension components)
   Array<Real> nodes = mesh.getNodes();
 
   //creating the iterators
   Array<Real>::vector_iterator it  = nodes.begin(spatial_dimension);
   Array<Real>::vector_iterator end = nodes.end(spatial_dimension);
 
   for (; it != end; ++it){
     Vector<Real> & coords = (*it);
 
     //do what you need
     ....
 
   }
 \end{cpp}
 In that example, each \code{Vector<Real>} is a geometrical array of
 size \code{spatial\_dimension} and the iteration is conveniently
 performed by the \code{Array} iterator.
 
 The \code{Array} object is intensively used to store second order
 tensor values.  In that case, it should be specified that the returned
 object type is a matrix when constructing the iterator. This is done
 when calling the \code{begin} function. For instance, assuming that we
 have a \code{Array} storing stresses, we can loop over the stored
 tensors by:
 
 \begin{cpp}
   //creating the iterators
   Array<Real>::matrix_iterator it  = stresses.begin(spatial_dimension,spatial_dimension);
   Array<Real>::matrix_iterator end = stresses.end(spatial_dimension,spatial_dimension);
 
   for (; it != end; ++it){
     Matrix<Real> & stress = (*it);
 
     //do what you need
     ....
 
   }
 \end{cpp}
 In that last example, the \code{Matrix} objects are
 \code{spatial\_dimension} $\times$ \code{spatial\_dimension} matrices.
 The light objects \code{Matrix} and \code{Vector} can be used and
 combined to do most common linear algebra. If the number of component
 is 1, it is possible to use a scalar\_iterator rather than the
 vector/matrix one.
 
 
 In general, a mesh consists of several kinds of elements.
 Consequently, the amount of data to be stored can differ for each
 element type.  The straightforward example is the connectivity array,
 namely the sequences of nodes belonging to each element (linear
 triangular elements have fewer nodes than, say, rectangular quadratic
 elements etc.).  A particular data structure called
 \code{ElementTypeMapArray} is provided to easily manage this kind of
 data.  It consists of a group of \code{Arrays}, each associated to an
 element type.  The following code can retrieve the
 \code{ElementTypeMapArray} which stores the connectivity arrays for a
 mesh:
 \begin{cpp}
   ElementTypeMapArray<UInt> & connectivities = mesh.getConnectivities();
 \end{cpp}
 Then, the specific array associated to a given element type can be obtained by
 \begin{cpp}
   Array<UInt> & connectivity_triangle = connectivities(_triangle_3);
 \end{cpp}
 where the first order 3-node triangular element was used in the presented piece
 of code.
 
 \subsection{Vector \& Matrix}
 
 The \code{Array} iterators as presented in the previous section can be shaped as
 \code{Vector} or \code{Matrix}. This objects represent $1^{st}$ and $2^{nd}$
 order tensors. As such they come with some functionalities that we will present
 a bit more into detail in this here.
 
 \subsubsection{\texttt{Vector<T>}}
 
 \begin{enumerate}
 \item Accessors:
   \begin{itemize}
   \item \code{v(i)} gives the $i^{th}$ component of the vector \code{v}
   \item \code{v[i]} gives the $i^{th}$ component of the vector \code{v}
   \item \code{v.size()} gives the number of component
   \end{itemize}
 \item Level 1: (results are scalars)
   \begin{itemize}
   \item \code{v.norm()} returns the geometrical norm ($L_2$)
   \item \code{v.norm<N>()} returns the $L_N$ norm defined as $\left(\sum_i
     |\code{v(i)}|^N\right)^{1/N}$. N can take any positive
     integer value. There are also some particular values for the most commonly
     used norms, \code{L\_1} for the Manhattan norm, \code{L\_2} for the
     geometrical norm and \code{L\_inf} for the norm infinity.
   \item \code{v.dot(x)} return the dot product of \code{v} and \code{x}
   \item \code{v.distance(x)} return the geometrical norm of $\code{v} -
     \code{x}$
   \end{itemize}
 \item Level 2: (results are vectors)
   \begin{itemize}
   \item \code{v += s}, \code{v -= s}, \code{v *= s}, \code{v /= s} those are
     element-wise operators that sum, substract, multiply or divide all the
     component of \code{v} by the scalar \code{s}
   \item \code{v += x}, \code{v -= x} sums or substracts the vector \code{x}
     to/from \code{v}
   \item \code{v.mul(A, x, alpha)} stores the result of $\alpha \mat{A} \vec{x}$ in \code{v},
     $\alpha$ is equal to 1 by default
   \item \code{v.solve(A, b)} stores the result of the resolution of the system $\mat{A} \vec{x} =
     \vec{b}$ in \code{v}
   \item \code{v.crossProduct(v1, v2)} computes the cross product of \code{v1} and \code{v2} and
     stores the result in \code{v}
   \end{itemize}
 \end{enumerate}
 
 \subsubsection{\texttt{Matrix<T>}}
 \begin{enumerate}
 \item Accessors:
   \begin{itemize}
   \item \code{A(i, j)} gives the component $A_{ij}$ of the matrix \code{A}
   \item \code{A(i)} gives the $i^{th}$ column of the matrix as a \code{Vector}
   \item \code{A[k]} gives the $k^{th}$ component of the matrix, matrices are
     stored in a column major way, which means that to access $A_{ij}$, $k = i +
     j M$
   \item \code{A.rows()} gives the number of rows of \code{A} ($M$)
   \item \code{A.cols()} gives the number of columns of \code{A} ($N$)
   \item \code{A.size()} gives the number of component in the matrix ($M \times N$)
   \end{itemize}
 \item Level 1: (results are scalars)
 
   \begin{itemize}
   \item \code{A.norm<N>()} returns the $L_N$ norm defined as $\left(\sum_i
     |\code{A(i)}|^N\right)^{1/N}$. N can take any positive
     integer value.
   \item \code{A.norm()} is equivalent to \code{A.norm<L\_2>()}
   \item \code{A.trace()} return the trace of \code{A}
   \item \code{A.det()} return the determinant of \code{A}
   \item \code{A.doubleDot(B)} return the double dot product of \code{A} and
     \code{B}, $\mat{A}:\mat{B}$
   \end{itemize}
 \item Level 3: (results are matrices)
   \begin{itemize}
   \item \code{A.eye(s)}, \code{Matrix<T>::eye(s)} fills/creates a matrix with
     the $s\mat{I}$ with $\mat{I}$ the identity matrix
   \item \code{A.inverse(B)} stores $\mat{B}^{-1}$ in \code{A}
   \item \code{A.transpose()} returns  $\mat{A}^{t}$
   \item \code{A.outerProduct(v1, v2)} stores $\vec{v_1} \vec{v_2}^{t}$ in
     \code{A}
   \item \code{C.mul<t\_A, t\_B>(A, B, alpha)}: stores the result of the product of
     \code{A} and code{B} time the scalar \code{alpha} in \code{C}. \code{t\_A}
     and \code{t\_B} are boolean defining if \code{A} and \code{B} should be
     transposed or not.
 
     \begin{tabular}{ccl}
       \toprule
       \code{t\_A} & \code{t\_B} & result \\
       \midrule
       false & false & $\mat{C} = \alpha \mat{A} \mat{B}$\\
       false & true  & $\mat{C} = \alpha \mat{A} \mat{B}^t$\\
       true  & false & $\mat{C} = \alpha \mat{A}^t \mat{B}$\\
       true  & true  & $\mat{C} = \alpha \mat{A}^t \mat{B}^t$\\
       \bottomrule
     \end{tabular}
   \item \code{A.eigs(d, V)} this method computes the eigenvalues and
     eigenvectors of \code{A} and store the results in \code{d} and \code{V} such
     that $\code{d(i)} = \lambda_i$ and $\code{V(i)} = \vec{v_i}$ with
     $\mat{A}\vec{v_i} = \lambda_i\vec{v_i}$ and $\lambda_1 > ... > \lambda_i >
     ... > \lambda_N$
   \end{itemize}
 \end{enumerate}
 
 \subsubsection{\texttt{Tensor3<T>}}
 Accessors:
 \begin{itemize}
   \item \code{t(i, j, k)} gives the component $T_{ijk}$ of the tensor \code{t}
   \item \code{t(k)} gives the $k^{th}$ two-dimensional tensor as a \code{Matrix}
   \item \code{t[k]} gives the $k^{th}$ two-dimensional tensor as a \code{Matrix}
 \end{itemize}
 
 
 \section{Manipulating group of nodes and/or elements\label{sect:common:groups}}
 \akantu provides the possibility to manipulate
 subgroups of elements and nodes.
 Any \code{ElementGroup} and/or \code{NodeGroup} must be managed
 by a \code{GroupManager}. Such a manager has the role to
 associate group objects to names. This is a useful feature,
 in particular for the application of the boundary conditions,
 as will be demonstrated in section \ref{sect:smm:boundary}.
 To most general group manager is the \code{Mesh} class
 which inheritates from the \code{GroupManager} class.
 
 For instance, the following code shows how to request
 an element group to a mesh:
 
 \begin{cpp}
   // request creation of a group of nodes
   NodeGroup & my_node_group = mesh.createNodeGroup("my_node_group");
   // request creation of a group of elements
   ElementGroup & my_element_group = mesh.createElementGroup("my_element_group");
 
   /* fill and use the groups */
 \end{cpp}
 
 
 \subsection{The \texttt{NodeGroup} object}
 A group of nodes is stored in \code{NodeGroup} objects.
 They are quite simple objects which store the indexes
 of the selected nodes in a \code{Array<UInt>}.
 Nodes are selected by adding them when calling
 \code{NodeGroup::add}. For instance you can select nodes
 having a positive $X$ coordinate with the following code:
 \begin{cpp}
   Array<Real> & nodes = mesh.getNodes();
   NodeGroup & group = mesh.createNodeGroup("XpositiveNode");
 
   Array<Real>::const_vector_iterator it  = nodes.begin(spatial_dimension);
   Array<Real>::const_vector_iterator end = nodes.end(spatial_dimension);
 
   UInt index = 0;
 
   for (; it != end ; ++it , ++index){
     const Vector<Real> & position = *it;
     if (position(0) > 0) group.add(index);
   }
 \end{cpp}
 
 \subsection{The \texttt{ElementGroup} object}
 A group of elements is stored in \code{ElementGroup} objects.
 Since a group can contain elements of various types
 the \code{ElementGroup} object stores indexes in
 a \code{ElementTypeMapArray<UInt>} object.
 Then elements can be added to the group by calling \code{addElement}.
 
 For instance, selecting the elements for which the barycenter of the nodes
 has a positive $X$ coordinate can be made with:
 
 \begin{cpp}
   ElementGroup & group = mesh.createElementGroup("XpositiveElement");
 
   Mesh::type_iterator it  = mesh.firstType();
   Mesh::type_iterator end = mesh.lastType();
 
   Vector<Real> barycenter(spatial_dimension);
 
   for(; it != end; ++it){
     UInt nb_element  = mesh.getNbElement(*it);
     for(UInt e = 0; e < nb_element; ++e) {
       ElementType type = *it;
       mesh.getBarycenter(e, type, barycenter.storage());
       if (barycenter(0) > 0) group.add(type,e);
     }
   }
 \end{cpp}
 
 
 
 %%% Local Variables:
 %%% mode: latex
 %%% TeX-master: "manual"
 %%% End:
diff --git a/doc/manual/manual-io.tex b/doc/manual/manual-io.tex
index 07100480e..a053b08f8 100644
--- a/doc/manual/manual-io.tex
+++ b/doc/manual/manual-io.tex
@@ -1,200 +1,270 @@
 \chapter{Input/Output}\index{I\/O}
 
-\section{Generic data}
+\section{Input file \label{sect:parser}}
+
+The text input file of a simulation should be precised using the method \code{initialize} which will instantiate the static \code{Parser} object of \akantu. This section explains how to manipulate \code{Parser} objects to input data in \akantu.      
+\begin{cpp}
+int main(int argc, char *argv[]) {
+  initialize("input_files.dat", argc, argv);
+  ...
+\end{cpp}
+
+\subsection{Akantu Parser}
+
+\akantu file parser has a tree organization.
+\begin{itemize}
+\item \code{Parser}, the root of the tree, can be accessed using
+\begin{cpp}
+Parser & parser = getStaticParser();
+\end{cpp}
+\item \code{ParserSection}, branch of the tree, contains map a of sub-sections (\code{SectionType}, \code{ParserSection}) and a \code{ParserSection *} pointing to the parent section. The user section of the input file can directly be accessed by
+\begin{cpp}
+const ParserSection & usersect = getUserParser();
+\end{cpp} 
+\item \code{ParserParameter}, the leaf of the tree, carries data of the input file which can be casted to the correct type with
+\begin{cpp}
+Real mass = usersect.getParameter("mass");
+\end{cpp}
+or used directly within an expression
+\begin{cpp}
+Real dead_weight = 9.81 * usersect.getParameterValue<Real>("mass");
+\end{cpp}
+\end{itemize}
+
+\subsection{Grammar}
+
+The structure of text input files consists of different sections containing a list of parameters. As example, the file parsed in the previous section will look like
+\begin{cpp}
+  user parameters [
+     mass = 10.5
+  ]
+\end{cpp}  
+Basically every standard arithmetic operations can be used inside of input files as well as the constant \code{pi} and \code{e} and the exponent operator \code{\^{}}. Operations between \code{ParserParameter} are also possible with the convention that only parameters of the current and the parent sections are available. \code{Vector} and \code{Matrix} can also be read according to the \code{NumPy}\cite{numpy} writing convention (a.e. cauchy$\_$stress$\_$tensor = [[$\sigma_{xx}$, $\sigma_{xy}$],[$\sigma_{yx}$,$\sigma_{yy}$]]). 
+An example illustrating how to parse the following input file can be found in \code{example$\backslash$io$\backslash$parser$\backslash$example$\_$parser.cc}.
+\begin{cpp}
+user parameters [
+   spatial$\_$dimension = 2
+   mesh$\_$file = swiss$\_$cheese.msh
+   inner$\_$holes = holes
+   outter$\_$crust = crust
+   lactostatic$\_$p = 30e3
+   stress = [[lactostatic$\_$p,0],[0,lactostatic$\_$p]]
+   max$\_$nb$\_$iterations = 100
+   precision = 1e-9
+]
+\end{cpp}
+\subsection{Material section \label{sect:io:material}}
+The input file should also be used to specify material characteristics (constitutive behavior and material properties). The dedicated material section is then read by \code{initFull} method of \code{SolidMechanicsModel} which initializes the different materials specified with the following convention:
+\begin{cpp}
+  material %\emph{constitutive\_law}% %\emph{<optional flavor>}% [
+     name = $value$
+     rho = $value$
+     ...
+  ]
+\end{cpp}
+\index{Constitutive\_laws} where \emph{constitutive\_law} is the adopted
+constitutive law, followed by the material properties listed one by line in the
+bracket (\eg \code{name} and density \code{rho}). Some constitutive laws can
+also have an \emph{optional flavor}. More information can be found in sections relative to material
+constitutive laws \ref{sect:smm:CL} or in Appendix \ref{app:material-parameters}.
+
+\section{Output data}
+
+\subsection{Generic data}
 In this chapter, we address ways to get the internal data in human-readable formats.
 The models in \akantu handle data associated to the
 mesh, but this data can be split into several \code{Arrays}. For example, the
 data stored per element type in a \code{ElementTypeMapArray} is composed of as
 many \code{Array}s as types in the mesh.
 
 In order to get this data in a visualization software, the models contain a
 object to dump \code{VTK} files. These files can be visualized in software such
 as \code{ParaView}\cite{paraview}, \code{ViSit}\cite{visit} or \code{Mayavi}\cite{mayavi}.
 
 The internal dumper of the model can be configured to specify which data fields
 are to be written. This is done with the
 \code{addDumpField}\index{I\/O!addDumpField} method. By default all the files
 are generated in a folder called \code{paraview/}
 
 \begin{cpp}
   model.setBaseName("output"); // prefix for all generated files
 
   model.addDumpField("displacement");
   model.addDumpField("stress");
   ...
 
   model.dump()
 \end{cpp}
 
 The fields are dumped with the number of components of the memory. For example, in 2D, the memory has 
 \code{Vector}s of 2 components, or the $2^{nd}$ order tensors with $2\times2$ components.  
 This memory can be dealt with \code{addDumpFieldVector}\index{I\/O!addDumpFieldVector} which always dumps
 \code{Vector}s with 3 components or \code{addDumpFieldTensor}\index{I\/O!addDumpFieldTensor} which dumps $2^{nd}$
 order tensors with $3\times3$ components respectively. The routines \code{addDumpFieldVector}\index{I\/O!addDumpFieldVector} and
 \code{addDumpFieldTensor}\index{I\/O!addDumpFieldTensor} were introduced because of \code{ParaView} which mostly manipulate 3D data.
 
 Those fields which are stored by quadrature point are modified to be seen in the
 \code{VTK} file as elemental data. To do this, the default is to average the
 values of all the quadrature points.
 
 The list of fields depends on the models (for
 \code{SolidMechanicsModel} see table~\ref{tab:io:smm_field_list}).
 
 \begin{table}
   \centering
   \begin{tabular}{llll}
     \toprule
     key          &    type      & support \\
     \midrule
     displacement & Vector<Real> & nodes  \\
     mass         & Vector<Real> & nodes  \\
     velocity     & Vector<Real> & nodes  \\
     acceleration & Vector<Real> & nodes  \\
     force	       & Vector<Real> & nodes  \\
     residual     & Vector<Real> & nodes  \\
     increment     & Vector<Real> & nodes  \\
     {blocked\_dofs}  & Vector<bool> & nodes  \\    
     partitions   & Real         & elements \\
     material\_index & variable  & elements\\    
     strain & Matrix<Real> & quadrature points  \\
     Green strain & Matrix<Real> & quadrature points  \\
     principal strain & Vector<Real> & quadrature points  \\
     principal Green strain & Vector<Real> & quadrature points  \\
     grad\_u & Matrix<Real> & quadrature points  \\    
     stress & Matrix<Real> & quadrature points  \\
     Von Mises stress & Real & quadrature points  \\        
     material\_index & variable  & quadrature points \\
     \bottomrule
   \end{tabular}
   \caption{List of dumpable fields for \code{SolidMechanicsModel}.}
   \label{tab:io:smm_field_list}
 \end{table}
 
-\section{Cohesive elements' data}
+\subsection{Cohesive elements' data}
 Cohesive elements and their relative data can be easily dumped thanks
 to a specific dumper contained in
 \code{SolidMechanicsModelCohesive}. In order to use it, one has just
 to add the string \code{"cohesive elements"} when calling each method
 already illustrated. Here is an example on how to dump displacement
 and damage:
 \begin{cpp}
   model.setBaseNameToDumper("cohesive elements", "cohesive_elements_output");
   model.addDumpFieldVectorToDumper("cohesive elements", "displacement");
   model.addDumpFieldToDumper("cohesive elements", "damage");
   ...
 
   model.dump("cohesive elements");
 \end{cpp}
 
-\subsection{Fragmentation data}
+\subsubsection{Fragmentation data}
 
 Whenever the \code{SolidMechanicsModelCohesive} is used, it is
 possible to dump additional data about the fragments that get formed
 in the simulation both in serial and parallel. This task is carried
 out by the \code{FragmentManager} class, that takes care of computing
 the following quantities for each fragment:
 \begin{itemize}
 \item index;
 \item mass;
 \item moments of inertia;
 \item velocity;
 \item number of elements.
 \end{itemize}
 These computations can be realized at once by calling the function
 \code{computeAllData}, or individually by calling the other public
 functions of the class. The data can be dumped to be visualized in
 Paraview, or can be accessed within the simulation. An example of
 usage is:
 \begin{cpp}
   FragmentManager fragment_manager(model);
   fragment_manager.buildAllData();
   ...
 
   model.addDumpField("fragments");       // this field contains the indices
   model.addDumpField("fragments mass");
   model.addDumpField("moments of inertia");
   model.addDumpField("elements per fragment");
   ...
 
   for (UInt step = 1; step <= total_steps; ++step) {
     ...
 
     fragment_manager.buildAllData();
     model.dump();
   }
   ...
 
   const Array<Real> & fragment_velocities = fragment_manager.getVelocity();
   ...
 \end{cpp}
 At the end of this example the velocities of the fragments are
 accessed with a reference to a \code{const Array<Real>}. The size of
 this array is the number of fragments, and its number of components is
 the spatial dimension in this case.
 
 
-\section{Advanced dumping}
+\subsection{Advanced dumping}
 
-\subsection{Arbitrary fields}
+\subsubsection{Arbitrary fields}
 In addition to the predetermined fields from the models and materials, the user
 can add any data to a dumper as long as the support is the same. That is to say
 data that have the size of the full mesh on if the dumper is dumping the mesh,
 or of the size of an element group if it is a filtered dumper.
 
 For this the easiest is to use the ``external'' fields register functions\index{I\/O!addDumpFieldExternal}
 
 The simple case force nodal and elemental data are to pass directly the data
 container itself if it as the good size.
 \begin{itemize}
 \item For nodal fields :
 \begin{cpp}
   Array<T> nodal_data(nb_nodes, nb_component);
   ...
   model.addDumpFieldExternal("my_field", nodal_data);
 \end{cpp}
 
 \item For elemental fields :
 \begin{cpp}
   ElementTypeMapArray<T> elem_data;
   ...
   model.addDumpFieldExternal("my_field", elem_data);
 \end{cpp}
 \end{itemize}
 
 If some changes have to be applied on the data as for example a padding for
 \code{ParaView} vectors, this can be done by using the
 field interface.
 
 \begin{cpp}
   model.addDumpFieldExternal(const std::string & field_name,
                              dumper::Field * field);
 \end{cpp}
 
 An example of code presenting this interface is present in the
 \shellcode{\examplesdir/io/dumper}. This interface is part of the
 \code{Dumpable} class from which the \code{Mesh} inherits.
 
-\subsection{Creating a new dumper}
+\subsubsection{Creating a new dumper}
 
 You can also create you own dumpers, \akantu uses a third-party library in order
 to write the output files, \code{IOHelper}. \akantu supports the \code{ParaView}
 format and a Text format defined by \code{IOHelper}.
 
 This two files format are handled by the classes
 \code{DumperParaview}\index{I\/O!DumperParaview} and
 \code{DumperText}\index{I\/O!DumperText}.
 
 In order to use them you can instantiate on of this object in your code. This
 dumper have a simple interface. You can register a mesh
 \code{registerMesh}\index{I\/O!registerMesh},
 \code{registerFilteredMesh}\index{I\/O!registerFilteredMesh} or a field,
 \code{registerField}\index{I\/O!registerField}.
 
 An example of code presenting this low level interface is present in the
 \shellcode{\examplesdir/io/dumper}. The different types of \code{Field} that can
 be created are present in the source folder \shellcode{src/io/dumper}.
 
 %%% Local Variables:
 %%% mode: latex
 %%% TeX-master: "manual"
 %%% End: