getfem-commits
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Getfem-commits] (no subject)


From: Tetsuo Koyama
Subject: [Getfem-commits] (no subject)
Date: Thu, 14 Feb 2019 06:23:52 -0500 (EST)

branch: fixmisspell
commit a79e4ad543aa49fc8b02cfd1fc3a0a9a2e7e3fd2
Author: Tetsuo Koyama <address@hidden>
Date:   Thu Feb 14 20:22:57 2019 +0900

    fix sphinx compile warning
---
 doc/sphinx/source/about.rst                        |   4 +-
 doc/sphinx/source/biblio.rst                       |  16 +-
 doc/sphinx/source/bugs.rst                         |  12 +-
 doc/sphinx/source/copyright.rst                    |   2 +-
 doc/sphinx/source/glossary.rst                     |   4 +-
 doc/sphinx/source/gmm/blas.rst                     |   8 +-
 doc/sphinx/source/gmm/blas_interface.rst           | 170 +++++++++----------
 doc/sphinx/source/gmm/catch.rst                    |  16 +-
 doc/sphinx/source/gmm/denseqr.rst                  |   8 +-
 doc/sphinx/source/gmm/first-step.rst               |   8 +-
 doc/sphinx/source/gmm/index.rst                    |   2 +-
 doc/sphinx/source/gmm/inside.rst                   |  56 +++----
 doc/sphinx/source/gmm/install.rst                  |   4 +-
 doc/sphinx/source/gmm/iter.rst                     |  20 +--
 doc/sphinx/source/gmm/matrix.rst                   |   6 +-
 doc/sphinx/source/gmm/sub-matrix.rst               |   8 +-
 doc/sphinx/source/gmm/triangular.rst               |   2 +-
 doc/sphinx/source/install/install_linux.rst        |   2 +-
 doc/sphinx/source/install/install_mac.rst          |   4 +-
 doc/sphinx/source/install/install_windows.rst      |  12 +-
 doc/sphinx/source/links.rst                        |  10 +-
 doc/sphinx/source/lists.rst                        |   2 +-
 doc/sphinx/source/matlab/examples.rst              | 142 ++++++++--------
 doc/sphinx/source/matlab/install.rst               |   2 +-
 doc/sphinx/source/matlab/mlabgf.rst                |   4 +-
 doc/sphinx/source/matlab/oocmd.rst                 |   2 +-
 doc/sphinx/source/matlab/plotcmdref.rst            |  46 ++---
 doc/sphinx/source/project/contribute.rst           |   4 +-
 doc/sphinx/source/project/femdesc.rst              |  32 ++--
 doc/sphinx/source/project/intro.rst                |  20 +--
 doc/sphinx/source/project/libdesc.rst              |   2 +-
 doc/sphinx/source/project/libdesc_cont.rst         |   4 +-
 doc/sphinx/source/project/libdesc_fem.rst          |   2 +-
 .../source/project/libdesc_high_gen_assemb.rst     |  18 +-
 doc/sphinx/source/project/libdesc_im.rst           |   6 +-
 doc/sphinx/source/project/libdesc_interface.rst    |  70 ++++----
 doc/sphinx/source/project/libdesc_mesh.rst         |   2 +-
 doc/sphinx/source/project/libdesc_model.rst        |   2 +-
 doc/sphinx/source/python/examples.rst              | 186 ++++++++++-----------
 doc/sphinx/source/python/howtos.rst                |   2 +-
 doc/sphinx/source/python/intro.rst                 |   6 +-
 doc/sphinx/source/scilab/install.rst               |   2 +-
 doc/sphinx/source/scilab/intro.rst                 |   6 +-
 doc/sphinx/source/scilab/plotcmdref.rst            |  46 ++---
 doc/sphinx/source/scilab/scilabgf.rst              |   4 +-
 doc/sphinx/source/screenshots/shots.rst            |  90 +++++-----
 doc/sphinx/source/tutorial/basic_usage.rst         |   2 +-
 doc/sphinx/source/tutorial/install.rst             |   4 +-
 doc/sphinx/source/tutorial/thermo_coupling.rst     | 106 ++++++------
 doc/sphinx/source/tutorial/wheel.rst               |   8 +-
 doc/sphinx/source/userdoc/appendixA.rst            |  10 +-
 doc/sphinx/source/userdoc/appendixB.rst            |   2 +-
 doc/sphinx/source/userdoc/bfem.rst                 |   2 +-
 doc/sphinx/source/userdoc/bmesh.rst                |   8 +-
 doc/sphinx/source/userdoc/computeL2H1.rst          |   2 +-
 doc/sphinx/source/userdoc/convect.rst              |  20 +--
 doc/sphinx/source/userdoc/gasm_high.rst            |  80 +++++----
 doc/sphinx/source/userdoc/index.rst                |   2 +-
 doc/sphinx/source/userdoc/interMM.rst              |  16 +-
 doc/sphinx/source/userdoc/intro.rst                |   2 +-
 doc/sphinx/source/userdoc/model.rst                |   4 +-
 doc/sphinx/source/userdoc/model_ALE_rotating.rst   |  16 +-
 doc/sphinx/source/userdoc/model_Nitsche.rst        |  12 +-
 doc/sphinx/source/userdoc/model_bilaplacian.rst    |   6 +-
 .../source/userdoc/model_contact_friction.rst      |  18 +-
 .../model_contact_friction_large_sliding.rst       |  22 +--
 doc/sphinx/source/userdoc/model_continuation.rst   |  12 +-
 doc/sphinx/source/userdoc/model_dirichlet.rst      |  10 +-
 .../source/userdoc/model_generic_assembly.rst      |   6 +-
 .../source/userdoc/model_nonlinear_elasticity.rst  |   8 +-
 doc/sphinx/source/userdoc/model_object.rst         |  44 ++---
 .../userdoc/model_plasticity_small_strain.rst      |  18 +-
 .../source/userdoc/model_time_integration.rst      |  30 ++--
 doc/sphinx/source/userdoc/parallel.rst             |   6 +-
 doc/sphinx/source/userdoc/rmesh.rst                |  12 +-
 doc/sphinx/source/userdoc/xfem.rst                 |   2 +-
 doc/sphinx/source/whatsnew/1.7.rst                 |   2 +-
 doc/sphinx/source/whatsnew/2.0.rst                 |   2 +-
 doc/sphinx/source/whatsnew/3.0.rst                 |   2 +-
 doc/sphinx/source/whatsnew/4.1.1.rst               |   2 +-
 doc/sphinx/source/whatsnew/4.1.rst                 |   2 +-
 doc/sphinx/source/whatsnew/4.2.rst                 |   2 +-
 doc/sphinx/source/whatsnew/4.3.rst                 |   4 +-
 doc/sphinx/source/whatsnew/5.0.rst                 |   2 +-
 doc/sphinx/source/whatsnew/5.1.rst                 |   2 +-
 doc/sphinx/source/whatsnew/5.3.rst                 |   2 +-
 86 files changed, 791 insertions(+), 797 deletions(-)

diff --git a/doc/sphinx/source/about.rst b/doc/sphinx/source/about.rst
index 7b68393..ce69297 100644
--- a/doc/sphinx/source/about.rst
+++ b/doc/sphinx/source/about.rst
@@ -3,7 +3,7 @@ About these documents
 =====================
 
 These documents are generated from `reStructuredText
-<http://docutils.sourceforge.net/rst.html>`_ sources by *Sphinx*, a 
+<http://docutils.sourceforge.net/rst.html>`_ sources by *Sphinx*, a
 document processor specifically written for the Python documentation.
 
 In the online version of these documents, you can submit comments and
@@ -20,7 +20,7 @@ Many thanks go to:
 * The `Docutils <http://docutils.sourceforge.net/>`_ project for creating
   reStructuredText and the Docutils suite;
 * Fredrik Lundh for his `Alternative Python Reference
-  <http://effbot.org/zone/pyref.htm>`_ project from which Sphinx got many 
+  <http://effbot.org/zone/pyref.htm>`_ project from which Sphinx got many
   good ideas.
 
 See :ref:`reporting-bugs` for information how to report bugs in GetFEM++
diff --git a/doc/sphinx/source/biblio.rst b/doc/sphinx/source/biblio.rst
index 17d1ae8..d1089c2 100644
--- a/doc/sphinx/source/biblio.rst
+++ b/doc/sphinx/source/biblio.rst
@@ -30,12 +30,12 @@ References
 .. [BE-CO-DU2010] M. Bergot, G. Cohen, M. |durufle|.
    *Higher-order finite elements for hybrid meshes using new nodal pyramidal 
elements*
    J. Sci. Comput., 42, 345-381, 2010.
-   
+
 .. [br-ba-fo1989] F. Brezzi, K.J. Bathe, M. Fortin.
    *Mixed-interpolated element for Reissner-Mindlin plates*. Internat. J. 
Numer. Methods Engrg., 28, 1787-1801, 1989.
 
 .. [bu-ha2010] E. Burman, P. Hansbo.
-   *Fictitious domain finite element methods using cut elements: I. A 
stabilized Lagrange multiplier method*. Computer Methods in Applied Mechanics, 
199:41-44, 2680-2686, 2010. 
+   *Fictitious domain finite element methods using cut elements: I. A 
stabilized Lagrange multiplier method*. Computer Methods in Applied Mechanics, 
199:41-44, 2680-2686, 2010.
 
 .. [ca-re-so1994] D. Calvetti, L. Reichel and D.C. Sorensen.
    *An implicitely restarted Lanczos method for large symmetric eigenvalue 
problems*. Electronic Transaction on Numerical Analysis}. 2:1-21, 1994.
@@ -44,7 +44,7 @@ References
    *Crack-tip enrichment in the Xfem method using a cut-off function*. Int. J. 
Numer. Meth. Engng., 75(6):629-646, 2008.
 
 .. [CH-LA-RE2011] E. Chahine, P. Laborde, Y. Renard.
-   *A non-conformal eXtended Finite Element approach: Integral matching Xfem*. 
Applied Numerical Mathematics, 61:322-343, 2011. 
+   *A non-conformal eXtended Finite Element approach: Integral matching Xfem*. 
Applied Numerical Mathematics, 61:322-343, 2011.
 
 .. [ciarlet1978] P.G. Ciarlet.
    *The finite element method for elliptic problems*. Studies in Mathematics 
and its Applications vol. 4, North-Holland, 1978.
@@ -85,7 +85,7 @@ References
 .. [GR-GH1999] R.D. Graglia, I.-L. Gheorma.
    *Higher order interpolatory vector bases on pyramidal elements*
    IEEE transactions on antennas and propagation, 47:5, 775-782, 1999.
-   
+
 .. [GR-ST2015] D. Grandi, U. Stefanelli.
    *The Souza-Auricchio model for shape-memory alloys*
    Discrete and Continuous Dynamical Systems, Series S, 8(4):723-747, 2015.
@@ -118,10 +118,10 @@ References
    *An unconstrained integral approximation of large sliding frictional 
contact between deformable solids*. Computers and Structures, 153:75-90, 2015.
 
 .. [LA-RE2006] P. Laborde, Y. Renard.
-   *Fixed point strategies for elastostatic frictional contact problems*. 
Math. Meth. Appl. Sci., 31:415-441, 2008. 
+   *Fixed point strategies for elastostatic frictional contact problems*. 
Math. Meth. Appl. Sci., 31:415-441, 2008.
 
 .. [Li-Re2014] T. |ligursky| and Y. Renard.
-   *A Continuation Problem for Computing Solutions of Discretised Evolution 
Problems with Application to Plane Quasi-Static Contact Problems with 
Friction*. Comput. Methods Appl. Mech. Engrg. 280, 222-262, 2014. 
+   *A Continuation Problem for Computing Solutions of Discretised Evolution 
Problems with Application to Plane Quasi-Static Contact Problems with 
Friction*. Comput. Methods Appl. Mech. Engrg. 280, 222-262, 2014.
 
 .. [Li-Re2014hal] T. |ligursky| and Y. Renard.
    *Bifurcations in Piecewise-Smooth Steady-State Problems: Abstract Study and 
Application to Plane Contact Problems with Friction*. Computational Mechanics, 
56:1:39-62, 2015.
@@ -183,11 +183,11 @@ References
 
 .. |moes| unicode:: Mo U+00EB s
 .. |peric| unicode:: Peri U+0107
-.. |dolezel| unicode:: Dole U+017E el 
+.. |dolezel| unicode:: Dole U+017E el
    :rtrim:
 .. |ligursky| unicode:: Ligursk U+00FD
 .. |durufle| unicode:: Durufl U+00E9
-.. |solin| unicode:: U+0160 ol U+00ED n 
+.. |solin| unicode:: U+0160 ol U+00ED n
    :rtrim:
 
 
diff --git a/doc/sphinx/source/bugs.rst b/doc/sphinx/source/bugs.rst
index 3d34c7e..b3177c1 100644
--- a/doc/sphinx/source/bugs.rst
+++ b/doc/sphinx/source/bugs.rst
@@ -25,12 +25,12 @@ registration procedure. Otherwise, if you're not logged in, 
select
 report anonymously.
 
 Being now logged in, you can submit a bug.  Go to
-https://savannah.nongnu.org/projects/getfem and select the "Submit a new item" 
link 
+https://savannah.nongnu.org/projects/getfem and select the "Submit a new item" 
link
 (in the "Bug Tracker" table) to open the bug reporting form.
 
-The submission form has a number of fields.  For the "Summary" field, 
-enter a *very* short description of the problem; less than ten words is 
-good.  In the "Severity" field, select the severity of your problem; 
+The submission form has a number of fields.  For the "Summary" field,
+enter a *very* short description of the problem; less than ten words is
+good.  In the "Severity" field, select the severity of your problem;
 also select the "Privacy" and "Category" to which the bug relates.
 
 In the "Original Submission" field, describe the problem in detail,
@@ -45,7 +45,7 @@ each time action is taken on the bug.
 
 
 .. seealso::
-   
+
    * `How to Report Bugs Effectively
      <http://www.chiark.greenend.org.uk/~sgtatham/bugs.html>`_
      Article which goes into some detail about how to create a useful
@@ -54,6 +54,6 @@ each time action is taken on the bug.
 
    * `Bug Writing Guidelines
      <https://developer.mozilla.org/en/Bug_writing_guidelines>`_
-     Information about writing a good bug report.  Some of this is 
+     Information about writing a good bug report.  Some of this is
      specific to the Mozilla project, but describes general good
      practices.
diff --git a/doc/sphinx/source/copyright.rst b/doc/sphinx/source/copyright.rst
index 23a593a..2c90aa2 100644
--- a/doc/sphinx/source/copyright.rst
+++ b/doc/sphinx/source/copyright.rst
@@ -15,7 +15,7 @@ The text of this website and the documentations are available 
for modification a
 
 Publication director : Yves Renard, INSA Lyon, France.
 
-Editorial directors : Konstantinos Poulios, Technical University of Denmark  
and Franz Chouly, |Universite| de Franche |Comte|, |Besancon|, France. 
+Editorial directors : Konstantinos Poulios, Technical University of Denmark  
and Franz Chouly, |Universite| de Franche |Comte|, |Besancon|, France.
 
 The Website is hosted by GdS 2754 Mathrice, CNRS, France.
 
diff --git a/doc/sphinx/source/glossary.rst b/doc/sphinx/source/glossary.rst
index 0f843a7..d31009a 100644
--- a/doc/sphinx/source/glossary.rst
+++ b/doc/sphinx/source/glossary.rst
@@ -24,7 +24,7 @@ Glossary
       The degrees of freedom for a finite element method is the coefficients
       which multiply the shape functions in order to describe a
       (scalar or vector) field. Generally, they are the unknowns of the
-      problem in general. 
+      problem in general.
 
    Element
       An element is a small piece of a domain with a special shape (a segment,
@@ -66,4 +66,4 @@ Glossary
       For instance, the reference segment in |gf| is the segment [0,1].
       The reference triangle is the triangle (0,0), (0,1), (1,0). etc.
 
- 
+
diff --git a/doc/sphinx/source/gmm/blas.rst b/doc/sphinx/source/gmm/blas.rst
index 405aa31..d4c9aa7 100644
--- a/doc/sphinx/source/gmm/blas.rst
+++ b/doc/sphinx/source/gmm/blas.rst
@@ -34,13 +34,13 @@ transposition
 imaginary and real part
 -----------------------
 
-For a complex matrix ``M`` or a complex vector ``V``, 
+For a complex matrix ``M`` or a complex vector ``V``,
 ``gmm::real_part(M)``, ``gmm::real_part(V)``, ``gmm::imag_part(M)`` or 
``gmm::imag_part(V)`` give a possibily modifiable reference on the real or 
imaginary part of the matrix or vector (for instance 
``gmm::clear(gmm::imag_part(M))`` will set to zero the imaginary part of a 
matrix ``M``). These functions cannot be applied to real matrices or vectors.
 
 conjugate
 ---------
 
-For a matrix ``M`` or a vector ``V``, 
+For a matrix ``M`` or a vector ``V``,
 ``gmm::conjugated(M)`` and ``gmm::conjugated(V)`` give a constant reference on 
the conjugated vector or matrix. Of course, for a real vectors this has no 
effect (and no cost at all). Note : ``gmm::conjugated(M)`` transposes the 
matrix ``M`` so that this is the hermitian conjugate of ``M``. If you need only 
the conjugate of each component you have to use both transposition and 
conjugate with ``gmm::conjugated(gmm::transposed(M))`` or equivalently  
``gmm::transposed(gmm::conjugated(M))``.
 
 
@@ -90,9 +90,9 @@ Matrix-vector or matrix-matrix multiplication. Again, all the 
matrices and vecto
   gmm::col_matrix< gmm::wsvector<double> > M2(10, 10);
   gmm::col_matrix< gmm::vsvector<double> > M3(10, 10);
   ...
-  
+
   gmm::mult(M1, M2, M3); // M1 * M2 ---> M3
-  
+
   gmm::mult(gmm::sub_matrix(M1, sub_interval(0, 3)),
             gmm::sub_matrix(M2, sub_interval(4, 3)),
             gmm::sub_matrix(M3, sub_interval(2, 3)));
diff --git a/doc/sphinx/source/gmm/blas_interface.rst 
b/doc/sphinx/source/gmm/blas_interface.rst
index 0de1682..f391821 100644
--- a/doc/sphinx/source/gmm/blas_interface.rst
+++ b/doc/sphinx/source/gmm/blas_interface.rst
@@ -14,8 +14,8 @@ For better performance on dense matrices, it is possible to 
interface some opera
 
 to use this interface you have first to define ``GMM_USES_LAPACK`` before 
including |gmm| \ files::
 
-  \#define GMM_USES_LAPACK
-  \#include <gmm/gmm.h>
+  #define GMM_USES_LAPACK
+  #include <gmm/gmm.h>
 
   ... your code
 
@@ -35,97 +35,97 @@ Ask your system administrator if this configuration does 
not work.
 
 The following operations are interfaced::
 
-  vect_norm2(std::vector<T>)                                           
-                                                                        
-  vect_sp(std::vector<T>, std::vector<T>)                               
-  vect_sp(scaled(std::vector<T>), std::vector<T>)                       
-  vect_sp(std::vector<T>, scaled(std::vector<T>))                       
-  vect_sp(scaled(std::vector<T>), scaled(std::vector<T>))               
-                                                                        
-  vect_hp(std::vector<T>, std::vector<T>)                               
-  vect_hp(scaled(std::vector<T>), std::vector<T>)                       
-  vect_hp(std::vector<T>, scaled(std::vector<T>))                       
-  vect_hp(scaled(std::vector<T>), scaled(std::vector<T>))               
-                                                                        
-  add(std::vector<T>, std::vector<T>)                                   
-  add(scaled(std::vector<T>, a), std::vector<T>)                         
-
-  mult(dense_matrix<T>, dense_matrix<T>, dense_matrix<T>)               
-  mult(transposed(dense_matrix<T>), dense_matrix<T>, dense_matrix<T>)   
-  mult(dense_matrix<T>, transposed(dense_matrix<T>), dense_matrix<T>)   
-  mult(transposed(dense_matrix<T>), transposed(dense_matrix<T>),        
-       dense_matrix<T>)                                                 
-  mult(conjugated(dense_matrix<T>), dense_matrix<T>, dense_matrix<T>)   
-  mult(dense_matrix<T>, conjugated(dense_matrix<T>), dense_matrix<T>)   
-  mult(conjugated(dense_matrix<T>), conjugated(dense_matrix<T>),        
-       dense_matrix<T>)                                                 
-                                                                        
-  mult(dense_matrix<T>, std::vector<T>, std::vector<T>)                 
-  mult(transposed(dense_matrix<T>), std::vector<T>, std::vector<T>)     
-  mult(conjugated(dense_matrix<T>), std::vector<T>, std::vector<T>)     
-  mult(dense_matrix<T>, scaled(std::vector<T>), std::vector<T>)         
-  mult(transposed(dense_matrix<T>), scaled(std::vector<T>),             
-       std::vector<T>)                                                  
-  mult(conjugated(dense_matrix<T>), scaled(std::vector<T>),             
+  vect_norm2(std::vector<T>)
+
+  vect_sp(std::vector<T>, std::vector<T>)
+  vect_sp(scaled(std::vector<T>), std::vector<T>)
+  vect_sp(std::vector<T>, scaled(std::vector<T>))
+  vect_sp(scaled(std::vector<T>), scaled(std::vector<T>))
+
+  vect_hp(std::vector<T>, std::vector<T>)
+  vect_hp(scaled(std::vector<T>), std::vector<T>)
+  vect_hp(std::vector<T>, scaled(std::vector<T>))
+  vect_hp(scaled(std::vector<T>), scaled(std::vector<T>))
+
+  add(std::vector<T>, std::vector<T>)
+  add(scaled(std::vector<T>, a), std::vector<T>)
+
+  mult(dense_matrix<T>, dense_matrix<T>, dense_matrix<T>)
+  mult(transposed(dense_matrix<T>), dense_matrix<T>, dense_matrix<T>)
+  mult(dense_matrix<T>, transposed(dense_matrix<T>), dense_matrix<T>)
+  mult(transposed(dense_matrix<T>), transposed(dense_matrix<T>),
+       dense_matrix<T>)
+  mult(conjugated(dense_matrix<T>), dense_matrix<T>, dense_matrix<T>)
+  mult(dense_matrix<T>, conjugated(dense_matrix<T>), dense_matrix<T>)
+  mult(conjugated(dense_matrix<T>), conjugated(dense_matrix<T>),
+       dense_matrix<T>)
+
+  mult(dense_matrix<T>, std::vector<T>, std::vector<T>)
+  mult(transposed(dense_matrix<T>), std::vector<T>, std::vector<T>)
+  mult(conjugated(dense_matrix<T>), std::vector<T>, std::vector<T>)
+  mult(dense_matrix<T>, scaled(std::vector<T>), std::vector<T>)
+  mult(transposed(dense_matrix<T>), scaled(std::vector<T>),
+       std::vector<T>)
+  mult(conjugated(dense_matrix<T>), scaled(std::vector<T>),
        std::vector<T>)
 
-  mult_add(dense_matrix<T>, std::vector<T>, std::vector<T>)             
-  mult_add(transposed(dense_matrix<T>), std::vector<T>, std::vector<T>) 
-  mult_add(conjugated(dense_matrix<T>), std::vector<T>, std::vector<T>) 
-  mult_add(dense_matrix<T>, scaled(std::vector<T>), std::vector<T>)     
-  mult_add(transposed(dense_matrix<T>), scaled(std::vector<T>),         
-           std::vector<T>)                                              
-  mult_add(conjugated(dense_matrix<T>), scaled(std::vector<T>),         
-           std::vector<T>)                                              
-                                                                        
-  mult(dense_matrix<T>, std::vector<T>, std::vector<T>, std::vector<T>) 
-  mult(transposed(dense_matrix<T>), std::vector<T>, std::vector<T>,     
-       std::vector<T>)                                                  
-  mult(conjugated(dense_matrix<T>), std::vector<T>, std::vector<T>,     
-       std::vector<T>)                                                  
-  mult(dense_matrix<T>, scaled(std::vector<T>), std::vector<T>,         
-       std::vector<T>)                                                  
-  mult(transposed(dense_matrix<T>), scaled(std::vector<T>),             
-       std::vector<T>, std::vector<T>)                                  
-  mult(conjugated(dense_matrix<T>), scaled(std::vector<T>),             
-       std::vector<T>, std::vector<T>)                                  
-  mult(dense_matrix<T>, std::vector<T>, scaled(std::vector<T>),         
-       std::vector<T>)                                                  
-  mult(transposed(dense_matrix<T>), std::vector<T>,                     
-       scaled(std::vector<T>), std::vector<T>)                          
-  mult(conjugated(dense_matrix<T>), std::vector<T>,                     
-       scaled(std::vector<T>), std::vector<T>)                          
-  mult(dense_matrix<T>, scaled(std::vector<T>), scaled(std::vector<T>), 
-    std::vector<T>)                                                     
-  mult(transposed(dense_matrix<T>), scaled(std::vector<T>),             
-       scaled(std::vector<T>), std::vector<T>)                          
-  mult(conjugated(dense_matrix<T>), scaled(std::vector<T>),             
-       scaled(std::vector<T>), std::vector<T>)                          
-                                                                        
-  lower_tri_solve(dense_matrix<T>, std::vector<T>, k, b)                
-  upper_tri_solve(dense_matrix<T>, std::vector<T>, k, b)                
-  lower_tri_solve(transposed(dense_matrix<T>), std::vector<T>, k, b)    
-  upper_tri_solve(transposed(dense_matrix<T>), std::vector<T>, k, b)    
-  lower_tri_solve(conjugated(dense_matrix<T>), std::vector<T>, k, b)    
-  upper_tri_solve(conjugated(dense_matrix<T>), std::vector<T>, k, b)    
-                                                                        
-  lu_factor(dense_matrix<T>, std::vector<int>)                          
-  lu_solve(dense_matrix<T>, std::vector<T>, std::vector<T>)             
-  lu_solve(dense_matrix<T>, std::vector<int>, std::vector<T>,           
-           std::vector<T>)                                              
+  mult_add(dense_matrix<T>, std::vector<T>, std::vector<T>)
+  mult_add(transposed(dense_matrix<T>), std::vector<T>, std::vector<T>)
+  mult_add(conjugated(dense_matrix<T>), std::vector<T>, std::vector<T>)
+  mult_add(dense_matrix<T>, scaled(std::vector<T>), std::vector<T>)
+  mult_add(transposed(dense_matrix<T>), scaled(std::vector<T>),
+           std::vector<T>)
+  mult_add(conjugated(dense_matrix<T>), scaled(std::vector<T>),
+           std::vector<T>)
+
+  mult(dense_matrix<T>, std::vector<T>, std::vector<T>, std::vector<T>)
+  mult(transposed(dense_matrix<T>), std::vector<T>, std::vector<T>,
+       std::vector<T>)
+  mult(conjugated(dense_matrix<T>), std::vector<T>, std::vector<T>,
+       std::vector<T>)
+  mult(dense_matrix<T>, scaled(std::vector<T>), std::vector<T>,
+       std::vector<T>)
+  mult(transposed(dense_matrix<T>), scaled(std::vector<T>),
+       std::vector<T>, std::vector<T>)
+  mult(conjugated(dense_matrix<T>), scaled(std::vector<T>),
+       std::vector<T>, std::vector<T>)
+  mult(dense_matrix<T>, std::vector<T>, scaled(std::vector<T>),
+       std::vector<T>)
+  mult(transposed(dense_matrix<T>), std::vector<T>,
+       scaled(std::vector<T>), std::vector<T>)
+  mult(conjugated(dense_matrix<T>), std::vector<T>,
+       scaled(std::vector<T>), std::vector<T>)
+  mult(dense_matrix<T>, scaled(std::vector<T>), scaled(std::vector<T>),
+    std::vector<T>)
+  mult(transposed(dense_matrix<T>), scaled(std::vector<T>),
+       scaled(std::vector<T>), std::vector<T>)
+  mult(conjugated(dense_matrix<T>), scaled(std::vector<T>),
+       scaled(std::vector<T>), std::vector<T>)
+
+  lower_tri_solve(dense_matrix<T>, std::vector<T>, k, b)
+  upper_tri_solve(dense_matrix<T>, std::vector<T>, k, b)
+  lower_tri_solve(transposed(dense_matrix<T>), std::vector<T>, k, b)
+  upper_tri_solve(transposed(dense_matrix<T>), std::vector<T>, k, b)
+  lower_tri_solve(conjugated(dense_matrix<T>), std::vector<T>, k, b)
+  upper_tri_solve(conjugated(dense_matrix<T>), std::vector<T>, k, b)
+
+  lu_factor(dense_matrix<T>, std::vector<int>)
+  lu_solve(dense_matrix<T>, std::vector<T>, std::vector<T>)
+  lu_solve(dense_matrix<T>, std::vector<int>, std::vector<T>,
+           std::vector<T>)
   lu_solve_transposed(dense_matrix<T>, std::vector<int>, std::vector<T>,
-           std::vector<T>)                                              
-  lu_inverse(dense_matrix<T>)                                           
-  lu_inverse(dense_matrix<T>, std::vector<int>, dense_matrix<T>)        
-                                                                        
-  qr_factor(dense_matrix<T>, dense_matrix<T>, dense_matrix<T>) 
+           std::vector<T>)
+  lu_inverse(dense_matrix<T>)
+  lu_inverse(dense_matrix<T>, std::vector<int>, dense_matrix<T>)
+
+  qr_factor(dense_matrix<T>, dense_matrix<T>, dense_matrix<T>)
 
   implicit_qr_algorithm(dense_matrix<T>, std::vector<T>)
   implicit_qr_algorithm(dense_matrix<T>, std::vector<T>,
-                        dense_matrix<T>)                               
+                        dense_matrix<T>)
   implicit_qr_algorithm(dense_matrix<T>, std::vector<std::complex<T> >)
   implicit_qr_algorithm(dense_matrix<T>, std::vector<std::complex<T> >,
-                        dense_matrix<T>)                               
+                        dense_matrix<T>)
 
 
 Of course, it is not difficult to interface another operation if needed.
diff --git a/doc/sphinx/source/gmm/catch.rst b/doc/sphinx/source/gmm/catch.rst
index fe76a1f..ade2a86 100644
--- a/doc/sphinx/source/gmm/catch.rst
+++ b/doc/sphinx/source/gmm/catch.rst
@@ -11,16 +11,16 @@ Catch errors
 ============================
 
 
-Errors used in |gmm| are defined in the file ``gmm/gmm\_except.h``. In order 
to make easier  the error catching all errors derive from the type 
``std::logic\_error`` defined in the file `` stdexcept`` of the S.T.L.
+Errors used in |gmm| are defined in the file ``gmm/gmm_except.h``. In order to 
make easier  the error catching all errors derive from the type 
``std::logic_error`` defined in the file `` stdexcept`` of the S.T.L.
 
-A standard procedure, ``GMM\_STANDARD\_CATCH\_ERROR``, is defined in 
``gmm/gmm\_except.h``. This procedure catches all errors and print the error 
message when an error occurs. It can be used in the main procedure of the 
program as follows::
+A standard procedure, ``GMM_STANDARD_CATCH_ERROR``, is defined in 
``gmm/gmm_except.h``. This procedure catches all errors and print the error 
message when an error occurs. It can be used in the main procedure of the 
program as follows::
 
-  int main(void) \{ 
-    try \{ 
-      ... main program ... 
-        \} 
-     GMM\_STANDARD\_CATCH\_ERROR;
-  \}
+  int main(void) {
+    try {
+      ... main program ...
+        }
+     GMM_STANDARD_CATCH_ERROR;
+  }
 
 
 It is highly recommended to catch the errors at least in the main function, 
because if you do not so, you will not be able to see error messages.
diff --git a/doc/sphinx/source/gmm/denseqr.rst 
b/doc/sphinx/source/gmm/denseqr.rst
index 4d9891d..2e7f3d4 100644
--- a/doc/sphinx/source/gmm/denseqr.rst
+++ b/doc/sphinx/source/gmm/denseqr.rst
@@ -19,10 +19,10 @@ The following procedures are available in the file 
``gmm/gmm\_dense\_qr.h`` for
   implicit_qr_algorithm(M, eigval, double tol = 1E-16) // compute the
      // eigenvalues of M using the implicit QR factorisation (Householder and
      // Francis QR step version). eigval should be a vector of appropriate size
-     // in which the eigenvalues will be computed. If the matrix have 
+     // in which the eigenvalues will be computed. If the matrix have
      // complex eigenvalues, please use a complex vector.
 
-  implicit_qr_algorithm(M, eigval, shvect, double tol = 1E-16) // idem, 
+  implicit_qr_algorithm(M, eigval, shvect, double tol = 1E-16) // idem,
      // compute additionally the schur vectors in the matrix shvect.
 
   symmetric_qr_algorithm(M, eigval, double tol = 1E-16) // idem for symmetric
@@ -33,8 +33,8 @@ The following procedures are available in the file 
``gmm/gmm\_dense\_qr.h`` for
 
 
 
-`Remark`: The computation of eigenvectors for non hermitian matrices is not 
yet implemented. You can use for the moment the functions 
``geev_interface_left`` and ``geev_interface_right`` from the LAPACK interface 
(see ``gmm/gmm_lapack_interface.h``. These LAPACK functions compute right and 
left eigen vectors. 
-                   
+`Remark`: The computation of eigenvectors for non hermitian matrices is not 
yet implemented. You can use for the moment the functions 
``geev_interface_left`` and ``geev_interface_right`` from the LAPACK interface 
(see ``gmm/gmm_lapack_interface.h``. These LAPACK functions compute right and 
left eigen vectors.
+
 
 The following function defined in the file ``gmm/gmm\_condition\_number.h``::
 
diff --git a/doc/sphinx/source/gmm/first-step.rst 
b/doc/sphinx/source/gmm/first-step.rst
index d3950e4..bdab584 100644
--- a/doc/sphinx/source/gmm/first-step.rst
+++ b/doc/sphinx/source/gmm/first-step.rst
@@ -21,8 +21,8 @@ It is not possible in |gmm| to invert all kind of matrices. 
For the moment, the
   gmm::scale(M, 2.0);                    // M = 2 * Id.
   M(1,2) = 1.0;
 
-  gmm::copy(M, M2);  
- 
+  gmm::copy(M, M2);
+
   gmm::lu_inverse(M);
 
   gmm::mult(M, M2, M3);
@@ -43,7 +43,7 @@ You have more than one possibility to solve a linear system. 
If you have a dense
 
   std::vector<double> X(3), B(3), Bagain(3);
   B[0] = 1.0; B[1] = 2.0; B[2] = 3.0;  // B = [1 2 3]
- 
+
   gmm::lu_solve(M, X, B);
 
   gmm::mult(M, X, Bagain);
@@ -52,7 +52,7 @@ You have more than one possibility to solve a linear system. 
If you have a dense
 
 
 If, now, you have a sparse system coming for example from a pde 
discretization, you have various iterative solvers, with or without 
preconditioners. This is an example with a precontionned GMRES::
- 
+
   int nbdof = 1000; // number of degrees of freedom.
   gmm::row_matrix< gmm::rsvector<double> > M(nbdof, nbdof); // a sparse matrix
   std::vector<double> X(nbdof), B(nbdof); // Unknown and left hand side.
diff --git a/doc/sphinx/source/gmm/index.rst b/doc/sphinx/source/gmm/index.rst
index a3abc38..e7fe148 100644
--- a/doc/sphinx/source/gmm/index.rst
+++ b/doc/sphinx/source/gmm/index.rst
@@ -31,4 +31,4 @@ Gmm++ Library
    qd
    first-step
    inside
-   noverif
\ No newline at end of file
+   noverif
diff --git a/doc/sphinx/source/gmm/inside.rst b/doc/sphinx/source/gmm/inside.rst
index fd1b0d6..92f60c8 100644
--- a/doc/sphinx/source/gmm/inside.rst
+++ b/doc/sphinx/source/gmm/inside.rst
@@ -31,7 +31,7 @@ For a vector type, the following informations are available::
   typename gmm::linalg_traits<V>::linalg_type    --> should be abstract_vector
   typename gmm::linalg_traits<V>::index_sorted    --> linalg_true or 
linalg_false
   typename gmm::linalg_traits<V>::const_iterator --> const iterator to iterate 
on the
-                                                     components of the vector 
in 
+                                                     components of the vector 
in
                                                      order to read them.
   typename gmm::linalg_traits<V>::iterator       --> iterator to iterate on the
                                                      components of the vector 
in
@@ -43,7 +43,7 @@ For a vector type, the following informations are available::
   typename gmm::linalg_traits<V>::origin_type    --> the type of vector itself
                                                      or the type of referenced
                                                      vector for a reference.
- 
+
   gmm::linalg_traits<V>::size(v)     --> a method which gives the size of the 
vector.
   gmm::linalg_traits<V>::begin(v)    --> a method which gives an iterator on 
the
                                          beginning of the vector
@@ -54,11 +54,11 @@ For a vector type, the following informations are 
available::
 
   gmm::linalg_traits<V>::access(o, it, ite, i) --> return the ith component or 
a
                                           reference on the ith component. o is 
a
-                                          pointer o type ``origin_type *'' or 
+                                          pointer o type ``origin_type *'' or
                                           ``const origin_type *''.
 
   gmm::linalg_traits<V>::clear(o, it, ite) --> clear the vector. o is a
-                                          pointer o type ``origin_type *'' or 
+                                          pointer o type ``origin_type *'' or
                                           ``const origin_type *''.
 
 and for a matrix type::
@@ -77,13 +77,13 @@ and for a matrix type::
                                                       row_and_col or 
col_and_row.
   typename gmm::linalg_traits<M>::sub_col_type      --> type of reference on a 
column
                                                       (if the matrix is not 
row_major)
-  typename gmm::linalg_traits<M>::const_sub_col_type --> type of const 
reference on a 
+  typename gmm::linalg_traits<M>::const_sub_col_type --> type of const 
reference on a
                                                        column
   typename gmm::linalg_traits<M>::col_iterator      --> iterator on the columns
   typename gmm::linalg_traits<M>::const_col_iterator --> const iterator on the 
columns
   typename gmm::linalg_traits<M>::sub_row_type      --> type of reference on a 
row
                                                       (if the matrix is not 
col_major)
-  typename gmm::linalg_traits<M>::const_sub_row_type --> type of const 
reference on a 
+  typename gmm::linalg_traits<M>::const_sub_row_type --> type of const 
reference on a
                                                        row
   typename gmm::linalg_traits<M>::const_row_iterator --> const iterator on the 
rows
   typename gmm::linalg_traits<M>::row_iterator       --> iterator on the rows
@@ -106,7 +106,7 @@ and for a matrix type::
                                           iterator  (if not row_major)
   gmm::linalg_traits<M>::origin(m)    --> gives a void pointer allowing to 
identify
                                           the matrix
-  gmm::linalg_traits<M>::access(it,i) --> return the ith component or a 
reference 
+  gmm::linalg_traits<M>::access(it,i) --> return the ith component or a 
reference
                                           on the ith component of the row or
                                           column pointed by it.
   gmm::linalg_traits<M>::do_clear(m)  --> make a clear on the matrix
@@ -120,17 +120,13 @@ How to iterate on the components of a vector
 
 Here is an example which accumulate the components of a vector. It is assumed 
that ``V`` is a vector type and ``v`` an instantiated vector::
 
-  
   typename gmm::linalg_traits<V>::value_type r(0); // scalar in which we 
accumulate
-  typename gmm::linalg_traits<V>::const_iterator it = vect_const_begin(v); // 
beginning 
-                                                                           // 
of v
+  typename gmm::linalg_traits<V>::const_iterator it = vect_const_begin(v); // 
beginning of v
   typename gmm::linalg_traits<V>::const_iterator ite = vect_const_end(v); // 
end of v
 
   for (; it != ite; ++it)  // loop on the components
     r += *it;              // accumulate the components
 
-
-
 This piece of code will work with every kind of interfaced vector.
 
 For sparse or skyline vectors, it is possible to obtain the index of the 
components pointed by the iterator with ``it.index()``. Here is the example of 
the scalar product of two sparse or skyline vectors, assuming ``V1`` and ``V2`` 
are two vector types and ``v1``, ``v2`` two corresponding instantiated vectors::
@@ -142,17 +138,17 @@ For sparse or skyline vectors, it is possible to obtain 
the index of the compone
    typename gmm::linalg_traits<V1>::value_type r(0); // it is assumed that V2 
have a
                                                 // compatible value_type
 
-   while (it1 != ite1 && it2 != ite2) \{  // loops on the components
-     if (it1.index() == it2.index()) \{
+   while (it1 != ite1 && it2 != ite2) {  // loops on the components
+     if (it1.index() == it2.index()) {
        res += (*it1) * (*it2));          // if the indices are equals 
accumulate
        ++it1;
        ++it2;
-     \}
+     }
      else if (it1.index() < it2.index())
        ++it1;
      else
        ++it2;
-   \}
+   }
 
 This algorithm use the fact that indices are increasing in a sparse vector. 
This code will not work for dense vectors because dense vector iterators do not 
have the method ``it.index()``.
 
@@ -163,27 +159,27 @@ You can iterate on the rows of a matrix if it is not a 
column major matrix and o
 
 If you need not to be optimal, you can use a basic loop like that::
 
-  for (size_t i = 0; i < gmm::mat_nrows(m); ++i) \{
+  for (size_t i = 0; i < gmm::mat_nrows(m); ++i) {
     typename gmm::linalg_traits<M>::const_sub_row_type row = mat_const_row(M, 
i);
 
     ...
 
     std::cout << "norm of row " << i << " : " << vect_norm2(row) << std::endl;
-  \}
+  }
 
 But you can also use iterators, like that::
 
   typename gmm::linalg_traits<M>::const_row_iterator it = 
mat_row_const_begin(m);
   typename gmm::linalg_traits<M>::const_row_iterator ite = 
mat_row_const_end(m);
 
-  for (; it != ite; ++it) \{
+  for (; it != ite; ++it) {
     typename gmm::linalg_traits<M>::const_sub_row_type
       row = gmm::linalg_traits<M>::row(it);
 
     ...
 
     std::cout << "norm of row " << i << " : " << vect_norm2(row) << std::endl;
-  \}
+  }
 
 
 How to make your algorithm working on all type of matrices
@@ -191,36 +187,36 @@ How to make your algorithm working on all type of matrices
 
 For this, you will generally have to specialize it. For instance, let us take 
a look at the code for ``gmm::nnz`` which count the number of stored components 
(in fact, the real ``gmm::nnz`` algorithm is specialized in most of the cases 
so that it does not count the components one by one)::
 
-  template <class L> inline size_type nnz(const L& l) \{
+  template <class L> inline size_type nnz(const L& l) {
     return nnz(l, typename linalg_traits<L>::linalg_type());
-  \}
+  }
 
-  template <class L> inline size_type nnz(const L& l, abstract_vector) \{ 
+  template <class L> inline size_type nnz(const L& l, abstract_vector) {
     typename linalg_traits<L>::const_iterator it = vect_const_begin(l);
     typename linalg_traits<L>::const_iterator ite = vect_const_end(l);
     size_type res(0);
     for (; it != ite; ++it) ++res;
     return res;
-  \}
+  }
 
-  template <class L> inline size_type nnz(const L& l, abstract_matrix) \{
+  template <class L> inline size_type nnz(const L& l, abstract_matrix) {
     return nnz(l,  typename principal_orientation_type<typename
                    linalg_traits<L>::sub_orientation>::potype());
-  \}
+  }
 
-  template <class L> inline size_type nnz(const L& l, row_major) \{
+  template <class L> inline size_type nnz(const L& l, row_major) {
     size_type res(0);
     for (size_type i = 0; i < mat_nrows(l); ++i)
       res += nnz(mat_const_row(l, i));
     return res;
-  \} 
+  }
 
-  template <class L> inline size_type nnz(const L& l, col_major) \{
+  template <class L> inline size_type nnz(const L& l, col_major) {
     size_type res(0);
     for (size_type i = 0; i < mat_ncols(l); ++i)
       res += nnz(mat_const_col(l, i));
     return res;
-  \}
+  }
 
 
 The first function dispatch on the second or the third function respectively 
if the parameter is a vector or a matrix. The third function dispatch again on 
the fourth and the fifth function respectively if the matrix is row_major or 
column major. Of course, as the function are declared ``inline``, at least the 
two dispatcher functions will not be implemented. Which means that this 
construction is not costly.
diff --git a/doc/sphinx/source/gmm/install.rst 
b/doc/sphinx/source/gmm/install.rst
index cbc396b..476a1d6 100644
--- a/doc/sphinx/source/gmm/install.rst
+++ b/doc/sphinx/source/gmm/install.rst
@@ -6,7 +6,7 @@
 
 .. _gmm-install:
 
-Installation 
+Installation
 ============
 
 Since we use standard GNU tools, the installation of the |gmm| library is 
somewhat standard.
@@ -47,4 +47,4 @@ the distribution.
 
 Now, to use |gmm| in you programs, the simpler manner is to include the file 
``gmm/gmm.h`` which includes all the template library. If the compilation time 
is too important, the minimum to be included is contained is the file 
``gmm/gmm\_kernel.h`` (vectors and matrix types, blas, sub vector and sub 
matrices).
 
-DO NOT FORGET to catch errors messages. See the corresponding section.
\ No newline at end of file
+DO NOT FORGET to catch errors messages. See the corresponding section.
diff --git a/doc/sphinx/source/gmm/iter.rst b/doc/sphinx/source/gmm/iter.rst
index 8673db9..67eabee 100644
--- a/doc/sphinx/source/gmm/iter.rst
+++ b/doc/sphinx/source/gmm/iter.rst
@@ -16,7 +16,7 @@ Most of the solvers provided in |gmm| come form ITL with 
slight modifications (g
 iterations
 ----------
 
-  The iteration object of |gmm| is a modification of the one in ITL. This is 
not a template type as in ITL. 
+  The iteration object of |gmm| is a modification of the one in ITL. This is 
not a template type as in ITL.
 
 The simplest initialization is::
 
@@ -27,14 +27,14 @@ Some possibilities::
 
   iter.set_noisy(n) // n = 0 : no output
                     // n = 1 : output of iterations on the standard output
-                    // n = 2 : output of iterations and sub-iterations 
+                    // n = 2 : output of iterations and sub-iterations
                     //         on the standard output
                     // ...
   iter.get_iteration() // after a computation, gives the number of
                        // iterations made.
   iter.converged()     // true if the method converged.
   iter.set_maxiter(n)  // Set the maximum of iterations.
-                       // A solver stops if the maximum of iteration is 
+                       // A solver stops if the maximum of iteration is
                        // reached, iter.converged() is then false.
 
 
@@ -51,7 +51,7 @@ Here is the list of available linear solvers::
   ...
   gmm::iteration iter(10E-9);// Iteration object with the max residu
   size_t restart = 50;       // restart parameter for GMRES
-  
+
   gmm::cg(A, X, B, PS, PR, iter); // Conjugate gradient
 
   gmm::bicgstab(A, X, B, PR, iter); // BICGSTAB BiConjugate Gradient Stabilized
@@ -72,15 +72,15 @@ Preconditioners
 
 The following preconditioners, to be used with linear solvers, are available::
 
-  gmm::identity_matrix P;   // No preconditioner 
+  gmm::identity_matrix P;   // No preconditioner
 
   gmm::diagonal_precond<matrix_type> P(SM); // diagonal preconditioner
- 
+
   gmm::mr_approx_inverse_precond<matrix_type> P(SM, 10, 10E-17);
                                                // preconditioner based on MR
                                                // iterations
 
-  gmm::ildlt_precond<matrix_type> P(SM); // incomplete (level 0) ldlt 
+  gmm::ildlt_precond<matrix_type> P(SM); // incomplete (level 0) ldlt
                                         // preconditioner. Fast to be
                                         // computed but less efficient than
                                         // gmm::ildltt_precond.
@@ -89,7 +89,7 @@ The following preconditioners, to be used with linear 
solvers, are available::
   // Efficient but could be costly.
   gmm::ildltt_precond<matrix_type> P(SM, k, threshold);
 
-  gmm::ilu_precond<matrix_type> P(SM);  // incomplete (level 0) ilu 
+  gmm::ilu_precond<matrix_type> P(SM);  // incomplete (level 0) ilu
                                         // preconditioner. Very fast to be
                                         // computed but less efficient than
                                         // gmm::ilut_precond.
@@ -100,7 +100,7 @@ The following preconditioners, to be used with linear 
solvers, are available::
   gmm::ilut_precond<matrix_type> P(SM, k, threshold);
 
   // incomplete LU with k fill-in, threshold and column pivoting 
preconditioner.
-  // Try it when ilut encounter too small pivots. 
+  // Try it when ilut encounter too small pivots.
   gmm::ilutp_precond<matrix_type> P(SM, k, threshold);
 
 
@@ -114,7 +114,7 @@ The additive Schwarz method is a decomposition domain 
method allowing the resolu
 For the moment, the method is not parallelized (this should be done ...). The 
call is the following::
 
  gmm::sequential_additive_schwarz(A, u, f, P, vB, iter, local_solver, 
global_solver)
-                           
+
 ``A`` is the matrix of the linear system. ``u`` is the unknown vector. ``f`` 
is the right hand side. ``P`` is an eventual preconditioner for the local 
solver. ``vB`` is a vector of rectangular sparse matrices (``of type const 
std::vector<vBMatrix>``, where ``vBMatrix`` is a sparse matrix type), each of 
these matrices is of size :math:`N \times N_i` where :math:`N` is the size of 
``A`` and :math:`N_i` the number of variables in the :math:`i^{th}` sub-domain 
; each column of the matrix is  [...]
 
 The test program ``schwarz_additive.C`` is the directory ``tests`` of GetFEM++ 
is an example of the resolution with the additive Schwarz method of an 
elastostatic problem with the use of coarse mesh to make a better 
preconditioning (i.e. one of the sub-domains represents in fact a coarser mesh).
diff --git a/doc/sphinx/source/gmm/matrix.rst b/doc/sphinx/source/gmm/matrix.rst
index a27fb22..2ba67aa 100644
--- a/doc/sphinx/source/gmm/matrix.rst
+++ b/doc/sphinx/source/gmm/matrix.rst
@@ -78,7 +78,7 @@ Those two type of matrices store an array of ``VECT`` so the 
memory is not conti
   gmm::row_matrix< std::vector<double> > M1(10, 10);  // dense row matrix
   gmm::col_matrix< gmm::wsvector<double> > M2(5, 20); // sparse column matrix
 
-Of course ``gmm::row_matrix<VECT>`` is a row matrix and it is impossible to 
access to a particular column of this matrix. 
+Of course ``gmm::row_matrix<VECT>`` is a row matrix and it is impossible to 
access to a particular column of this matrix.
 
 
 ``gmm::mat_nrows(M)`` gives the number of rows of a matrix and 
``gmm::mat_ncols(M)`` the number of columns.
@@ -98,7 +98,7 @@ sparse matrices
 ---------------
 
 Similarly, ``gmm::row_matrix< gmm::wsvector<double> >`` or ``gmm::col_matrix< 
gmm::rsvector<double> >`` represents some sparse matrices, but |gmm| provides 
also two types of classical sparse matrix types::
- 
+
   gmm::csr_matrix<T>
   gmm::csc_matrix<T>
 
@@ -119,4 +119,4 @@ Matrices ``gmm::csr_matrix<T>`` and ``gmm::csc_matrix<T>`` 
have the advantage to
   gmm::csc_matrix<double, 1> M1;
   gmm::csr_matrix<double, 1> M2;
 
-The ``1`` means that a shift will be done on all the indices.
\ No newline at end of file
+The ``1`` means that a shift will be done on all the indices.
diff --git a/doc/sphinx/source/gmm/sub-matrix.rst 
b/doc/sphinx/source/gmm/sub-matrix.rst
index e8da5a8..8bd0bd4 100644
--- a/doc/sphinx/source/gmm/sub-matrix.rst
+++ b/doc/sphinx/source/gmm/sub-matrix.rst
@@ -48,7 +48,7 @@ Now ``gmm::sub_vector(V, subi)`` gives a reference to a 
sub-vector::
 
   gmm::vsvector<double> V(10);
   V[5] = 3.0;
-  std::cout << gmm::sub_vector(V, gmm::sub_interval(2, 3)) << std::endl;  
+  std::cout << gmm::sub_vector(V, gmm::sub_interval(2, 3)) << std::endl;
 
 prints to the standard output ``V[2], V[3]`` and ``V[4]``.
 
@@ -57,17 +57,17 @@ prints to the standard output ``V[2], V[3]`` and ``V[4]``.
   gmm::col_matrix< gmm::wsvector<double> > M(5, 20);
   M(3, 2) = 5.0;
   std::cout << gmm::sub_matrix(M, gmm::sub_interval(2, 3), 
gmm::sub_interval(2, 3))
-            << std::endl;  
+            << std::endl;
 
 prints to the output a sub-matrix. If the two sub-indices are equal, it is 
possible to omit the second. For instance::
 
   gmm::col_matrix< gmm::wsvector<double> > M(5, 20);
   M(3, 2) = 5.0;
-  std::cout << gmm::sub_matrix(V, gmm::sub_interval(2, 3)) << std::endl;  
+  std::cout << gmm::sub_matrix(V, gmm::sub_interval(2, 3)) << std::endl;
 
 The reference on sub_matrix is writable if the corresponding matrix is 
writable (so you can copy on a sub_matrix, add sub-matrices ...).
 
 row and column of a matrix
 --------------------------
 
-``gmm::mat_row(M, i)`` gives a (possibly writable) reference to the row ``i`` 
of matrix ``M``, and ``gmm::mat_col(M, i)``  gives a (possibly writable) 
reference to the column ``i``. It is not possible to access to the rows if 
``M`` is a column matrix and to the columns if it is a row matrix. It is 
possible to use ``gmm::mat_const_row(M, i)`` and ``gmm::mat_const_col(M, i)`` 
to have constant references.
\ No newline at end of file
+``gmm::mat_row(M, i)`` gives a (possibly writable) reference to the row ``i`` 
of matrix ``M``, and ``gmm::mat_col(M, i)``  gives a (possibly writable) 
reference to the column ``i``. It is not possible to access to the rows if 
``M`` is a column matrix and to the columns if it is a row matrix. It is 
possible to use ``gmm::mat_const_row(M, i)`` and ``gmm::mat_const_col(M, i)`` 
to have constant references.
diff --git a/doc/sphinx/source/gmm/triangular.rst 
b/doc/sphinx/source/gmm/triangular.rst
index 058362b..0b63445 100644
--- a/doc/sphinx/source/gmm/triangular.rst
+++ b/doc/sphinx/source/gmm/triangular.rst
@@ -19,4 +19,4 @@ If ``M`` is a triangular matrix (upper or lower) and ``X`` a 
vector containing t
    gmm::lower_tri_solve(M, X, true)  // Solving a lower triangular system
                                      // assuming there is 1 on the diagonal
 
-components which are lower the diagonal are ignored by 
``gmm::upper_tri_solve`` and components which are upper the diagonal are 
ignored by ``gmm::lower_tri_solve``.
\ No newline at end of file
+components which are lower the diagonal are ignored by 
``gmm::upper_tri_solve`` and components which are upper the diagonal are 
ignored by ``gmm::lower_tri_solve``.
diff --git a/doc/sphinx/source/install/install_linux.rst 
b/doc/sphinx/source/install/install_linux.rst
index 71da831..5432483 100644
--- a/doc/sphinx/source/install/install_linux.rst
+++ b/doc/sphinx/source/install/install_linux.rst
@@ -157,7 +157,7 @@ If you want to use a different compiler than the one chosen 
automatically by the
 Once getfem is compiled:
 
   - Go to the scilab getfem++ interface install directory 
(interface/src/scilab if the installation is not done)
- 
+
   - launch scilab
 
   - load the getfem++ toolbox with:
diff --git a/doc/sphinx/source/install/install_mac.rst 
b/doc/sphinx/source/install/install_mac.rst
index e49a991..6a48759 100644
--- a/doc/sphinx/source/install/install_mac.rst
+++ b/doc/sphinx/source/install/install_mac.rst
@@ -38,7 +38,7 @@ For the sequential mumps,
 
    $ brew install mumps --without-mpi
 
-For the parallel one, just forget --without-mpi and install also mpi and 
metis. 
+For the parallel one, just forget --without-mpi and install also mpi and metis.
 
 For Qhull
 
@@ -245,7 +245,7 @@ If you want to use a different compiler than the one chosen 
automatically by the
 Once getfem is compiled:
 
   - Go to the scilab getfem++ interface install directory 
(interface/src/scilab if the installation is not done)
- 
+
   - launch scilab
 
   - load the getfem++ toolbox with:
diff --git a/doc/sphinx/source/install/install_windows.rst 
b/doc/sphinx/source/install/install_windows.rst
index 18e6b80..ed57106 100644
--- a/doc/sphinx/source/install/install_windows.rst
+++ b/doc/sphinx/source/install/install_windows.rst
@@ -70,7 +70,7 @@ possible with `Cygwin <https://www.cygwin.com/>`_.
       $ cd ..
 
     Do the same for LAPACK with::
-   
+
       $ cd lapack-?.?.?
       $ cp make.inc.example make.inc
       $ cp ../BLAS-?.?.?/blas_LINUX.a librefblas.a
@@ -81,14 +81,14 @@ possible with `Cygwin <https://www.cygwin.com/>`_.
   - You now need an installation of QHULL library in order to have access to
     the meshing and Xfem tools of |gf|
     (see `Qhull <http://www.qhull.org>`_ and
-    `Qhull install instructions <http://www.qhull.org/README.txt>`_). 
+    `Qhull install instructions <http://www.qhull.org/README.txt>`_).
     Download the sources of Qhull in you Msys home (similarly as what you have
     done for BLAS and LAPACK), untar them and enter into the Qhull source
     directory with Msys. You can compile and install Qhull simply with::
 
       $ make SO=dll
       $ make install
-  
+
   - Similarly, we will compile and install now a sequential version of
     the sparse linear solver `MUMPS <http://mumps.enseeiht.fr/>`_.
     Again, it is not strictly necessary since a version of
@@ -133,7 +133,7 @@ possible with `Cygwin <https://www.cygwin.com/>`_.
 Build with the Python interface
 *******************************
 
-Additionaly to build the Python interface, you will have first to install a 
64bits version of Python 2 or 3 on your system together with Numpy and Scipy 
packages. The simpler way is to install Anaconda2 or 3 (it already contains 
Numpy and Scipy packages which are necessary). For instance with Anaconda2  
+Additionaly to build the Python interface, you will have first to install a 
64bits version of Python 2 or 3 on your system together with Numpy and Scipy 
packages. The simpler way is to install Anaconda2 or 3 (it already contains 
Numpy and Scipy packages which are necessary). For instance with Anaconda2
 
   - Install Anaconda2 and add to windows path ``C:\install_dir\Anaconda2``
     and ``C:\install_dir\Anaconda2\Scripts`` where ``install_dir`` is the
@@ -198,13 +198,13 @@ Here follows the additional step to build the Matlab 
interface. You have first,
     path of the compiled interface (on the matlab command line) with::
 
       add_path('c:\msys\home\login\getfem-5.?\interface/tests/matlab')
-    
+
     and try the demo matlab programs of the interface in
     ``interface\tests\matlab``. In order not to have to call the ``addpath``
     command each time you open Matlab, you can add a Windows system variable
     ``MATLABPATH`` set to
     ``c:\msys\home\login\getfem-5.?\interface/tests/matlab``.
-    You can also move the ``interface\tests\matlab`` directory into your 
+    You can also move the ``interface\tests\matlab`` directory into your
     Matlab installation directory.
 
 
diff --git a/doc/sphinx/source/links.rst b/doc/sphinx/source/links.rst
index 1280aa4..e85943e 100644
--- a/doc/sphinx/source/links.rst
+++ b/doc/sphinx/source/links.rst
@@ -89,14 +89,14 @@ Examples of publications based on GetFEM++
 
 .. |link8_2| replace:: hal.archives-ouvertes.fr
 .. _link8_2: https://hal.archives-ouvertes.fr/hal-00960996/
-      
+
 * Fabre M., Pousin J. and Renard Y. A fictitious domain method for 
frictionless contact problems in elasticity using Nitsche's method. SMAI J. of 
Comput. Math., 2:19--50, 2016. |link8_2|_.
 
 .. |link8_3| replace:: hal.archives-ouvertes.fr
 .. _link8_3: https://hal.archives-ouvertes.fr/hal-00937569/
-  
+
 * Poulios K., Renard Y. An unconstrained integral approximation of large 
sliding frictional contact between deformable solids. Computers and Structures, 
153:75--90, 2015. |link8_3|_.
-      
+
 .. |link8_5| replace:: orbit.dtu.dk
 .. _link8_5: 
http://orbit.dtu.dk/ws/files/125153379/Postprint_Homogenization_of_long_fiber_reinforced_composites_including_fiber_bending_effects.pdf
 
@@ -104,7 +104,7 @@ Examples of publications based on GetFEM++
 
 .. |link10| replace:: sciencedirect.com
 .. _link10: https://www.sciencedirect.com/science/article/pii/S0950061818311292
-      
+
 * Lozinguez E., Barthélémy J.-F., Bouteiller V. and Desbois T., Contribution 
of Sacrificial Anode in reinforced concrete patch repair: Results of numerical 
simulations. Construction and Building Materials, 178 (2018) 405–417. |link10|_.
 
 .. |link11| replace:: link.springer.com
@@ -114,7 +114,7 @@ Examples of publications based on GetFEM++
 
 .. |link12| replace:: onlinelibrary.wiley.com
 .. _link12: https://onlinelibrary.wiley.com/doi/abs/10.1002/hbm.21479
-      
+
 * Windhoff M., Opitz A., and Thielscher A., Electric Field Calculations in 
Brain Stimulation Based on Finite Elements: An Optimized Processing Pipeline 
for the Generation and Usage of Accurate Individual Head Models. Human Brain 
Mapping, 34(4), 923-35, 2013. DOI: 10.1002/hbm.21479. |link12|_.
 
 
diff --git a/doc/sphinx/source/lists.rst b/doc/sphinx/source/lists.rst
index b8f7434..da4edc3 100644
--- a/doc/sphinx/source/lists.rst
+++ b/doc/sphinx/source/lists.rst
@@ -10,7 +10,7 @@ GetFEM++ is maintened on the `Savannah 
<http://Savannah.gnu.org>`_ collaborative
 
 
 
-The mailing lists of GetFEM++ are listed on the page `getfem mailing lists 
<https://savannah.nongnu.org/mail/?group=getfem>`_ 
+The mailing lists of GetFEM++ are listed on the page `getfem mailing lists 
<https://savannah.nongnu.org/mail/?group=getfem>`_
 
 The main mainling list is the `user one 
<https://lists.nongnu.org/mailman/listinfo/getfem-users>`_. All kind of 
problems or questions about using, install or improve GetFEM++ can be posted 
there. Don't forget to register to the list before to post a message.
 
diff --git a/doc/sphinx/source/matlab/examples.rst 
b/doc/sphinx/source/matlab/examples.rst
index 02b0441..276119a 100644
--- a/doc/sphinx/source/matlab/examples.rst
+++ b/doc/sphinx/source/matlab/examples.rst
@@ -14,15 +14,15 @@ Examples
 A step-by-step basic example
 ----------------------------
 
-This example shows the basic usage of getfem, on the über-canonical problem 
above 
-all others: solving the :envvar:`Laplacian`, :math:`-\Delta u = f` on a 
square, 
-with the Dirichlet condition :math:`u = g(x)` on the domain boundary. You can 
find 
-the **m-file** of this example under the name **demo_step_by_step.m** in the 
+This example shows the basic usage of getfem, on the über-canonical problem 
above
+all others: solving the :envvar:`Laplacian`, :math:`-\Delta u = f` on a square,
+with the Dirichlet condition :math:`u = g(x)` on the domain boundary. You can 
find
+the **m-file** of this example under the name **demo_step_by_step.m** in the
 directory ``interface/tests/matlab/`` of the |gf| distribution.
 
-The first step is to **create a mesh**. It is possible to create simple 
structured meshes or unstructured meshes for simple geometries (see 
``gf_mesh('generate', mesher_object mo, scalar h))``) or to rely on an external 
mesher (see ``gf_mesh('import', string 
-FORMAT, string FILENAME))``).  For this example, we 
-just consider a regular **cartesian mesh** whose nodes are 
+The first step is to **create a mesh**. It is possible to create simple 
structured meshes or unstructured meshes for simple geometries (see 
``gf_mesh('generate', mesher_object mo, scalar h))``) or to rely on an external 
mesher (see ``gf_mesh('import', string
+FORMAT, string FILENAME))``).  For this example, we
+just consider a regular **cartesian mesh** whose nodes are
 :math:`\{x_{i=0\ldots10,j=0..10}=(i/10,j/10)\}`::
 
   >> % creation of a simple cartesian mesh
@@ -31,37 +31,37 @@ just consider a regular **cartesian mesh** whose nodes are
        id: 0
       cid: 0
 
-If you try to look at the value of ``m``, you'll notice that it appears to be 
a 
-structure containing two integers. The first one is its identifier, the second 
one 
-is its class-id, i.e. an identifier of its type. This small structure is just 
an 
-"handle" or "descriptor" to the real object, which is stored in the |gf| 
memory 
-and cannot be represented via |Mlab| data structures. Anyway, you can still 
+If you try to look at the value of ``m``, you'll notice that it appears to be a
+structure containing two integers. The first one is its identifier, the second 
one
+is its class-id, i.e. an identifier of its type. This small structure is just 
an
+"handle" or "descriptor" to the real object, which is stored in the |gf| memory
+and cannot be represented via |Mlab| data structures. Anyway, you can still
 inspect the |gf| objects via the command ``gf_workspace('stats')``.
 
-Now we can try to have a **look at the mesh**, with its vertices numbering and 
the 
+Now we can try to have a **look at the mesh**, with its vertices numbering and 
the
 convexes numbering::
 
   >> % we enable vertices and convexes labels
   >> gf_plot_mesh(m, 'vertices', 'on', 'convexes', 'on');
 
-As you can see, the mesh is regular, and the numbering of its nodes and 
convexes 
-is also regular (this is guaranteed for cartesian meshes, but do not hope a 
+As you can see, the mesh is regular, and the numbering of its nodes and 
convexes
+is also regular (this is guaranteed for cartesian meshes, but do not hope a
 similar numbering for the degrees of freedom).
 
-The next step is to **create a mesh_fem object**. This one links a mesh with a 
set 
+The next step is to **create a mesh_fem object**. This one links a mesh with a 
set
 of FEM::
 
   >> % create a mesh_fem of for a field of dimension 1 (i.e. a scalar field)
   >> mf = gf_mesh_fem(m,1);
   >> gf_mesh_fem_set(mf,'fem',gf_fem('FEM_QK(2,2)'));
 
-The first instruction builds a new |mlab_mf| object, the second argument 
specifies 
-that this object will be used to interpolate scalar fields (since the unknown 
-:math:`u` is a scalar field). The second instruction assigns the :math:`Q^2` 
FEM 
+The first instruction builds a new |mlab_mf| object, the second argument 
specifies
+that this object will be used to interpolate scalar fields (since the unknown
+:math:`u` is a scalar field). The second instruction assigns the :math:`Q^2` 
FEM
 to every convex (each basis function is a polynomial of degree 4, remember that
-:math:`P^k\Rightarrow` polynomials of degree :math:`k`, while 
-:math:`Q^k\Rightarrow` polynomials of degree :math:`2k`). As :math:`Q^2` is a 
-polynomial FEM, you can view the expression of its basis functions on the 
+:math:`P^k\Rightarrow` polynomials of degree :math:`k`, while
+:math:`Q^k\Rightarrow` polynomials of degree :math:`2k`). As :math:`Q^2` is a
+polynomial FEM, you can view the expression of its basis functions on the
 reference convex::
 
   >> gf_fem_get(gf_fem('FEM_QK(2,2)'), 'poly_str');
@@ -76,11 +76,11 @@ reference convex::
       '-4*x*y + 4*x^2*y + 8*x*y^2 - 8*x^2*y^2'
       'x*y - 2*x^2*y - 2*x*y^2 + 4*x^2*y^2'
 
-It is also possible to make use of the "object oriented" features of |mlab|. 
As 
-you may have noticed, when a class "foo" is provided by the |gfi|, it is build 
-with the function ``gf_foo``, and manipulated with the functions 
``gf_foo_get`` 
-and ``gf_foo_set``. But (with matlab 6.x and better) you may also create the 
-object with the ``gfFoo`` constructor , and manipulated with the ``get(..)`` 
and 
+It is also possible to make use of the "object oriented" features of |mlab|. As
+you may have noticed, when a class "foo" is provided by the |gfi|, it is build
+with the function ``gf_foo``, and manipulated with the functions ``gf_foo_get``
+and ``gf_foo_set``. But (with matlab 6.x and better) you may also create the
+object with the ``gfFoo`` constructor , and manipulated with the ``get(..)`` 
and
 ``set(..)`` methods. For example, the previous steps could have been::
 
   >> gfFem('FEM_QK(2,2)');
@@ -96,28 +96,28 @@ object with the ``gfFoo`` constructor , and manipulated 
with the ``get(..)`` and
   gfMeshFem object: ID=1 [1316 bytes], qdim=1, nbdof=441,
     linked gfMesh object: dim=2, nbpts=121, nbcvs=100
 
-Now, in order to perform numerical integrations on ``mf``, we need to **build 
a 
+Now, in order to perform numerical integrations on ``mf``, we need to **build a
 mesh_im object**::
 
   >> % assign the same integration method on all convexes
   >> mim = gf_mesh_im(m, gf_integ('IM_EXACT_PARALLELEPIPED(2)'));
 
-The integration method will be used to compute the various integrals on each 
-element: here we choose to perform exact computations (no :envvar:`quadrature 
-formula`), which is possible since the geometric transformation of these 
convexes 
-from the reference convex is linear (this is true for all simplices, and this 
is 
-also true for the parallelepipeds of our regular mesh, but it is not true for 
-general quadrangles), and the chosen FEM is polynomial. Hence it is possible 
to 
-analytically integrate every basis function/product of basis 
-functions/gradients/etc. There are many alternative FEM methods and 
integration 
+The integration method will be used to compute the various integrals on each
+element: here we choose to perform exact computations (no :envvar:`quadrature
+formula`), which is possible since the geometric transformation of these 
convexes
+from the reference convex is linear (this is true for all simplices, and this 
is
+also true for the parallelepipeds of our regular mesh, but it is not true for
+general quadrangles), and the chosen FEM is polynomial. Hence it is possible to
+analytically integrate every basis function/product of basis
+functions/gradients/etc. There are many alternative FEM methods and integration
 methods (see :ref:`ud`).
 
-Note however that in the general case, approximate integration methods are a 
+Note however that in the general case, approximate integration methods are a
 better choice than exact integration methods.
 
-Now we have to **find the** ":envvar:`boundary`" **of the domain**, in order 
to 
-set a Dirichlet condition. A mesh object has the ability to store some sets of 
-convexes and convex faces. These sets (called "regions") are accessed via an 
+Now we have to **find the** ":envvar:`boundary`" **of the domain**, in order to
+set a Dirichlet condition. A mesh object has the ability to store some sets of
+convexes and convex faces. These sets (called "regions") are accessed via an
 integer #id::
 
   >> % detect the border of the mesh
@@ -126,7 +126,7 @@ integer #id::
   >> gf_mesh_set(m, 'region', 42, border);
   >> gf_plot_mesh(m, 'regions', [42]); % the boundary edges appears in red
 
-Here we find the faces of the convexes which are on the boundary of the mesh 
(i.e. 
+Here we find the faces of the convexes which are on the boundary of the mesh 
(i.e.
 the faces which are not shared by two convexes).
 
  Remark:
@@ -134,22 +134,22 @@ the faces which are not shared by two convexes).
    we could have used ``gf_mesh_get(m, 'OuTEr_faCes')``, as the interface is
    case-insensitive, and whitespaces can be replaced by underscores.
 
-The array ``border`` has two rows, on the first row is a convex number, on the 
-second row is a face number (which is local to the convex, there is no global 
+The array ``border`` has two rows, on the first row is a convex number, on the
+second row is a face number (which is local to the convex, there is no global
 numbering of faces). Then this set of faces is assigned to the region number 
42.
 
-At this point, we just have to desribe the model and run the solver to get the 
-solution! The ":envvar:`model`" is created with the ``gf_model`` (or 
``gfModel``) 
-constructor. A model is basically an object which build a global linear system 
-(tangent matrix for non-linear problems) and its associated right hand side.  
-Typical modifications are insertion of the stiffness matrix for the problem 
-considered (linear elasticity, laplacian, etc), handling of a set of 
contraints, 
-Dirichlet condition, addition of a source term to the right hand side etc. The 
-global tangent matrix and its right hand side are stored in the 
":envvar:`model`" 
+At this point, we just have to desribe the model and run the solver to get the
+solution! The ":envvar:`model`" is created with the ``gf_model`` (or 
``gfModel``)
+constructor. A model is basically an object which build a global linear system
+(tangent matrix for non-linear problems) and its associated right hand side.
+Typical modifications are insertion of the stiffness matrix for the problem
+considered (linear elasticity, laplacian, etc), handling of a set of 
contraints,
+Dirichlet condition, addition of a source term to the right hand side etc. The
+global tangent matrix and its right hand side are stored in the 
":envvar:`model`"
 structure.
 
-Let us build a problem with an easy solution: :math:`u=x(x-1)y(y-1)+x^5`, then 
we 
-have :math:`\Delta u=2(x^2+y^2)-2(x+y)+20x^3` (the FEM won't be able to catch 
the 
+Let us build a problem with an easy solution: :math:`u=x(x-1)y(y-1)+x^5`, then 
we
+have :math:`\Delta u=2(x^2+y^2)-2(x+y)+20x^3` (the FEM won't be able to catch 
the
 exact solution since we use a :math:`Q^2` method).
 
 We start with an empty real model::
@@ -157,16 +157,16 @@ We start with an empty real model::
   >> % empty real model
   >> md = gf_model('real');
 
-(a model is either ``'real'`` or ``'complex'``). And we declare that ``u`` is 
an 
+(a model is either ``'real'`` or ``'complex'``). And we declare that ``u`` is 
an
 unknown of the system on the finite element method `mf` by::
 
   >> % declare that "u" is an unknown of the system
   >> % on the finite element method `mf`
   >> gf_model_set(md, 'add fem variable', 'u', mf);
 
-Now, we add a "generic elliptic" brick, which handles 
:math:`-\nabla\cdot(A:\nabla 
-u) = \ldots` problems, where :math:`A` can be a scalar field, a matrix field, 
or 
-an order 4 tensor field. By default, :math:`A=1`. We add it on our main 
variable 
+Now, we add a "generic elliptic" brick, which handles 
:math:`-\nabla\cdot(A:\nabla
+u) = \ldots` problems, where :math:`A` can be a scalar field, a matrix field, 
or
+an order 4 tensor field. By default, :math:`A=1`. We add it on our main 
variable
 ``u`` with::
 
   >> % add generic elliptic brick on "u"
@@ -181,21 +181,21 @@ Next we add a Dirichlet condition on the domain boundary::
   >> gf_model_set(md, 'add Dirichlet condition with multipliers', mim, 'u', 
mf, 42, 'DirichletData');
 
 
-The two first lines defines a data of the model which represents the value of 
the 
-Dirichlet condition. The third one add a Dirichlet condition to the variable 
``u`` 
-on the boundary number ``42``. The dirichlet condition is imposed with 
lagrange 
-multipliers. Another possibility is to use a penalization. A |mlab_mf| 
argument is 
-also required, as the Dirichlet condition :math:`u=g` is imposed in a weak 
form 
-:math:`\int_\Gamma u(x)v(x) = \int_\Gamma g(x)v(x) ~ \forall v` where 
:math:`v` is 
+The two first lines defines a data of the model which represents the value of 
the
+Dirichlet condition. The third one add a Dirichlet condition to the variable 
``u``
+on the boundary number ``42``. The dirichlet condition is imposed with lagrange
+multipliers. Another possibility is to use a penalization. A |mlab_mf| 
argument is
+also required, as the Dirichlet condition :math:`u=g` is imposed in a weak form
+:math:`\int_\Gamma u(x)v(x) = \int_\Gamma g(x)v(x) ~ \forall v` where 
:math:`v` is
 taken in the space of multipliers given by here by ``mf``.
 
 
 .. topic:: Remark:
 
-   the polynomial expression was interpolated on ``mf``. It is possible only 
if 
-   ``mf`` is of Lagrange type. In this first example we use the same |mlab_mf| 
for 
-   the unknown and for the data such as ``g``, but in the general case, ``mf`` 
-   won't be Lagrangian and another (Lagrangian) |mf| will be used for the 
+   the polynomial expression was interpolated on ``mf``. It is possible only if
+   ``mf`` is of Lagrange type. In this first example we use the same |mlab_mf| 
for
+   the unknown and for the data such as ``g``, but in the general case, ``mf``
+   won't be Lagrangian and another (Lagrangian) |mf| will be used for the
    description of Dirichlet conditions, source terms etc.
 
 A source term can be added with the following lines::
@@ -205,13 +205,13 @@ A source term can be added with the following lines::
   >> gf_model_set(md, 'add initialized fem data', 'VolumicData', mf, f);
   >> gf_model_set(md, 'add source term brick', mim, 'u', 'VolumicData');
 
-It only remains now to launch the solver. The linear system is assembled and 
solve 
+It only remains now to launch the solver. The linear system is assembled and 
solve
 with the instruction::
 
   >> % solve the linear system
   >> gf_model_get(md, 'solve');
 
-The model now contains the solution (as well as other things, such as the 
linear 
+The model now contains the solution (as well as other things, such as the 
linear
 system which was solved). It is extracted, a display into a |mlab| figure::
 
   >> % extracted solution
@@ -329,9 +329,9 @@ Using Matlab Object-Oriented features
 
 The basic functions of the |gf| toolbox do not use any advanced |mlab| features
 (except that the handles to getfem objects are stored in a small |mlab|
-structure). But the toolbox comes with a set of |Mlab| objects, which 
encapsulate 
+structure). But the toolbox comes with a set of |Mlab| objects, which 
encapsulate
 the handles and make them look as real |mlab| objects. The aim is not to 
provide
-extra-functionalities, but to have a better integration of the toolbox with 
+extra-functionalities, but to have a better integration of the toolbox with
 |mlab|.
 
 Here is an example of its use::
diff --git a/doc/sphinx/source/matlab/install.rst 
b/doc/sphinx/source/matlab/install.rst
index 641a054..47de936 100644
--- a/doc/sphinx/source/matlab/install.rst
+++ b/doc/sphinx/source/matlab/install.rst
@@ -6,7 +6,7 @@
 
 .. _mlab-install:
 
-Installation 
+Installation
 ============
 
 The installation of the |gfi| toolbox can be somewhat tricky, since it 
combines a
diff --git a/doc/sphinx/source/matlab/mlabgf.rst 
b/doc/sphinx/source/matlab/mlabgf.rst
index 2061200..7fa6036 100644
--- a/doc/sphinx/source/matlab/mlabgf.rst
+++ b/doc/sphinx/source/matlab/mlabgf.rst
@@ -77,8 +77,8 @@ Functions
 Objects
 -------
 
-Various "objects" can be manipulated by the |gfm| toolbox, see fig. 
-:ref:`malb-fig-hierarchy`. The MESH and MESHFEM objects are the two most 
+Various "objects" can be manipulated by the |gfm| toolbox, see fig.
+:ref:`malb-fig-hierarchy`. The MESH and MESHFEM objects are the two most
 important objects.
 
 .. _malb-fig-hierarchy:
diff --git a/doc/sphinx/source/matlab/oocmd.rst 
b/doc/sphinx/source/matlab/oocmd.rst
index 9350748..46de8f1 100644
--- a/doc/sphinx/source/matlab/oocmd.rst
+++ b/doc/sphinx/source/matlab/oocmd.rst
@@ -91,7 +91,7 @@ induce a call to a mex-file, and additional |mlab| code::
 Hence you should always try to store data in |mlab| arrays instead of
 repetitively calling the getfem functions.
 
-Avalaible object types are :envvar:`gfCvStruct`, :envvar:`gfGeoTrans`, 
+Avalaible object types are :envvar:`gfCvStruct`, :envvar:`gfGeoTrans`,
 :envvar:`gfEltm`, :envvar:`gfInteg`, :envvar:`gfFem`, :envvar:`gfMesh`,
 :envvar:`gfMeshFem`, :envvar:`gfMeshIm`, :envvar:`gfMdBrick`,
 :envvar:`gfMdState`, :envvar:`gfModel`, :envvar:`gfSpmat`, :envvar:`gfPrecond`,
diff --git a/doc/sphinx/source/matlab/plotcmdref.rst 
b/doc/sphinx/source/matlab/plotcmdref.rst
index 0388312..e44e002 100644
--- a/doc/sphinx/source/matlab/plotcmdref.rst
+++ b/doc/sphinx/source/matlab/plotcmdref.rst
@@ -41,19 +41,19 @@ gf_plot
 
   The options are specified as pairs of "option name"/"option value"
 
-  'zplot',{'off'|'on'}       : values of ``U`` are mapped on the $z$-axis 
(only possible when qdim=1, mdim=2). 
+  'zplot',{'off'|'on'}       : values of ``U`` are mapped on the $z$-axis 
(only possible when qdim=1, mdim=2).
   'norm', {'off'|'on'}       : if qdim >= 2, color-plot the norm of the field
   'dir',[]                   : or the scalar product of the field with 'dir' 
(can be a vector, or 'x', 'y' etc..)
   'refine',8                 : nb of refinments for curved edges and surface 
plots
   'interpolated',{'off'|'on'}: if triangular patch are interpolated
   'pcolor',{'on'|'off'}      : if the field is scalar, a color plot of its 
values is plotted
-  'quiver',{'on'|'off'}      : if the field is vector, represent arrows        
       
+  'quiver',{'on'|'off'}      : if the field is vector, represent arrows        
   'quiver_density',50        : density of arrows in quiver plot
   'quiver_scale',1           : scaling of arrows (0=>no scaling)
   'mesh',{'off'|'on'}        : show the mesh ?
-  'meshopts',{cell(0)}          : cell array of options passed to 
gf_plot_slice for the mesh 
+  'meshopts',{cell(0)}          : cell array of options passed to 
gf_plot_slice for the mesh
   'deformed_mesh', {'off'|'on'} : shows the deformed mesh (only when qdim == 
mdim)
-  'deformed_meshopts', {cell(0)}: cell array of options passed to 
gf_plot_slice for the deformed mesh 
+  'deformed_meshopts', {cell(0)}: cell array of options passed to 
gf_plot_slice for the deformed mesh
   'deformation',[]           : plots on the deformed object (only when qdim == 
mdim)
   'deformation_mf',[]        : plots on the deformed object (only when qdim == 
mdim)
   'deformation_scale','10%'  : indicate the amplitude of the deformation. Can 
be a percentage of the mesh width if given as a string, or an absolute value if 
given as a number
@@ -81,14 +81,14 @@ gf_plot
   dof) labels.
 
   For example, plotting a scalar field on the border of a 3D mesh can be done 
with ::
-  
+
     % load the 'strange.mesh_fem' (found in the getfem_matlab/tests directory)
-    mf=gf_mesh_fem('load', 'strange.mesh_fem') 
+    mf=gf_mesh_fem('load', 'strange.mesh_fem')
     U=rand(1, gf_mesh_fem_get(mf, 'nbdof')); # random field that will be drawn
     gf_plot(mf, U, 'refine', 25, 'cvlst', gf_mesh_get(mf,'outer faces'), 
'mesh','on');
- 
 
- 
+
+
 
 gf_plot_1D
 -------------------------------------------
@@ -123,18 +123,18 @@ gf_plot_mesh
 
   gf_plot_mesh(m, ...)
 
-  'vertices', {'off'|'on'}    : displays also vertices numbers. 
-  'convexes', {'off'|'on'}    : displays also convexes numbers. 
+  'vertices', {'off'|'on'}    : displays also vertices numbers.
+  'convexes', {'off'|'on'}    : displays also convexes numbers.
   'dof',{'off'|'on'}          : displays also finite element nodes. In that 
case, ``m`` should be a ``mesh_fem`` identifier.
   'regions',BLST              : displays the boundaries listed in BLST.
   'cvlst',CVLST               : display only the listed convexes. If CVLST has 
two rows, display only the faces listed in the second row.
   'edges', {'on' | 'off'}     : display edges ?
   'faces', {'off'|'on'}       : fills each 2D-face of the mesh
   'curved', {'off'|'on'}      : displays curved edges
-  'refine',N                  : refine curved edges and filled faces N times  
+  'refine',N                  : refine curved edges and filled faces N times
   'deformation', Udef         : optionnal deformation applied to the mesh (M 
must be a mesh_fem object)
   'edges_color',[.6 .6 1]     : RGB values for the color of edges
-  'edges_width',1             : width of edges              
+  'edges_width',1             : width of edges
   'faces_color',[.75 .75 .75]): RGB values for the color of faces
   'quality',{ 'off' | 'on' }  : Display the quality of the mesh.
 
@@ -143,14 +143,14 @@ gf_plot_mesh
 
   This function is used to display a mesh.
 
-  Example :: 
+  Example ::
 
     % the mesh is in the tests directory of the distribution
     m=gf_mesh('import','gid','donut_with_quadratic_tetra_314_elements.msh');
-    gf_plot_mesh(m,'refine',15,'cvlst',gf_mesh_get(m,'outer 
faces'),'faces','on',\ldots, 'faces_color',[1. .9 
.2],'curved','on','edges_width',2); 
+    gf_plot_mesh(m,'refine',15,'cvlst',gf_mesh_get(m,'outer 
faces'),'faces','on',\ldots, 'faces_color',[1. .9 
.2],'curved','on','edges_width',2);
     camlight % turn on the light!
 
- 
+
 
 gf_plot_slice
 -------------------------------------------
@@ -178,7 +178,7 @@ gf_plot_slice
   pcolor, ['on']      : if the field is scalar, a color plot of its values is 
plotted
   quiver, ['on']      : if the field is vector, represent arrows
   quiver_density, 50  : density of arrows in quiver plot
-  quiver_scale, 1     : density of arrows in quiver plot 
+  quiver_scale, 1     : density of arrows in quiver plot
   tube, ['on']        : use tube plot for 'filar' (1D) parts of the slice
   tube_color, ['red'] : color of tubes (ignored if 'data' is not empty and 
'pcolor' is on)
   tube_radius, ['0.5%'] : tube radius; you can use a constant, or a percentage 
(of the mesh size) or a vector of nodal values
@@ -201,29 +201,29 @@ gf_plot_slice
     Usl=gf_compute(pde.mf_u,U,'interpolate on', sl);  % interpolate the 
solution on the slice
     % show the norm of the displacement on this slice
     
gf_plot_slice(sl,'mesh','on','data',sqrt(sum(Usl.^2,1)),'mesh_slice_edges','off');
-    
+
     % another slice: now we take the lower part of the mesh
     
sl=gf_slice(\{'boundary',\{'intersection',\{'planar',+1,[0;0;6],[0;0;-1]\},\ldots
                                             
\{'planar',+1,[0;0;0],[0;1;0]\}\}\},m,6);
     Usl=gf_compute(pde.mf_u,U,'interpolate on', sl);
     hold on;
     
gf_plot_slice(sl,'mesh','on','data',sqrt(sum(Usl.^2,1)),'mesh_slice_edges','off');
-    
+
     % this slice contains the transparent mesh faces displayed on the picture
     sl2=gf_slice(\{'boundary',\{'planar',+1,[0;0;0],[0;1;0]\}\},\ldots
                 m,6,setdiff(all_faces',TOPfaces','rows')');
-    gf_plot_slice(sl2,'mesh_faces','off','mesh','on','pcolor','off'); 
-    
+    gf_plot_slice(sl2,'mesh_faces','off','mesh','on','pcolor','off');
+
     % last step is to plot the streamlines
     hh=[1 5 9 12.5 16 19.5]; % vertical position of the different starting 
points of the streamlines
     H=[zeros(2,numel(hh));hh];
-    
+
     % compute the streamlines
     tsl=gf_slice('streamlines',pde.mf_u,U,H);
     Utsl=gf_compute(pde.mf_u,U,'interpolate on', tsl);
-    
+
     % render them with "tube plot"
-    
[a,h]=gf_plot_slice(tsl,'mesh','off','tube_radius',.2,'tube_color','white'); 
+    
[a,h]=gf_plot_slice(tsl,'mesh','off','tube_radius',.2,'tube_color','white');
     hold off;
     % use a nice colormap
     caxis([0 .7]);
diff --git a/doc/sphinx/source/project/contribute.rst 
b/doc/sphinx/source/project/contribute.rst
index 2fd533a..8e07ba9 100644
--- a/doc/sphinx/source/project/contribute.rst
+++ b/doc/sphinx/source/project/contribute.rst
@@ -96,7 +96,7 @@ You can now transfer your modifications to the Savannah 
repository with ::
 
   git push origin devel-me-rewrite-fem-kernel
 
-where of course *devel-me-rewrite-fem-kernel* is still the name of your 
branch. At this stage your modifications are registered in the branch 
*devel-me-rewrite-fem-kernel* of Savannah repository. 
+where of course *devel-me-rewrite-fem-kernel* is still the name of your 
branch. At this stage your modifications are registered in the branch 
*devel-me-rewrite-fem-kernel* of Savannah repository.
 Your role stops here, since you are not allowed to modify the master branch of 
|gf|.
 
 
@@ -136,7 +136,7 @@ The recommended way for new contributors to translate 
document is to join |tfweb
   cd doc/sphinx
   tx pull -l <lang>
 
-Set code for your native language to <lang> (see |cfvlang|_ ). 
+Set code for your native language to <lang> (see |cfvlang|_ ).
 
 .. warning::
 
diff --git a/doc/sphinx/source/project/femdesc.rst 
b/doc/sphinx/source/project/femdesc.rst
index 987fed4..5f9e6cb 100644
--- a/doc/sphinx/source/project/femdesc.rst
+++ b/doc/sphinx/source/project/femdesc.rst
@@ -103,9 +103,9 @@ The following functions build the descriptions:
 
    description of the parallelepiped of reference of dimension ``d``.
 
-The vertices correspond to the classical vertices for such reference element. 
For 
-instance the vertices for the triangle are :math:`(0, 0)`, :math:`(1, 0)` and 
-:math:`(0, 1)`. It corresponds to the configuration shown in Figure 
+The vertices correspond to the classical vertices for such reference element. 
For
+instance the vertices for the triangle are :math:`(0, 0)`, :math:`(1, 0)` and
+:math:`(0, 1)`. It corresponds to the configuration shown in Figure
 :ref:`dp-fig-elem`
 
 If ``p`` is of type |bg_pcr| then ``p->structure()`` is the corresponding 
convex
@@ -160,8 +160,8 @@ geometric nodes are denoted:
 
    g^i, i = 0, \ldots, n_g - 1.
 
-The geometric transformation is described thanks to a :math:`n_g` components 
-polynomial vector (In fact, as an extention, non polynomial geometric 
+The geometric transformation is described thanks to a :math:`n_g` components
+polynomial vector (In fact, as an extention, non polynomial geometric
 transformation can also be supported by |gf|, but this is very rarely used)
 
 .. math::
@@ -192,9 +192,9 @@ The derivative of :math:`\tau` is then
 
    \fbox{$K(\widehat{x}) := \nabla\tau(\widehat{x}) = G\cdot\nabla {\cal 
N}(\widehat{x})$,}
 
-where :math:`K(\widehat{x}) = \nabla\tau(\widehat{x})` is a :math:`N\times P` 
matrix and 
-:math:`\nabla {\cal N}(\widehat{x})` a :math:`n_g\times P` matrix. The 
(transposed) 
-pseudo-inverse of :math:`\nabla\tau(\widehat{x})` is a :math:`N\times P` 
matrix denoted 
+where :math:`K(\widehat{x}) = \nabla\tau(\widehat{x})` is a :math:`N\times P` 
matrix and
+:math:`\nabla {\cal N}(\widehat{x})` a :math:`n_g\times P` matrix. The 
(transposed)
+pseudo-inverse of :math:`\nabla\tau(\widehat{x})` is a :math:`N\times P` 
matrix denoted
 :math:`B(\widehat{x})`:
 
 .. math::
@@ -212,12 +212,12 @@ where ``"name of trans"`` can be chosen among the 
following list.
 
 * ``"GT_PK(n,k)"``
 
-  Description of the simplex transformation of dimension ``n`` and degree 
``k`` 
+  Description of the simplex transformation of dimension ``n`` and degree ``k``
   (Most of the time, the degree 1 is used).
 
 * ``"GT_QK(n,k)"``
 
-  Description of the parallelepiped transformation of dimension ``n`` and 
degree 
+  Description of the parallelepiped transformation of dimension ``n`` and 
degree
   ``k``.
 
 * ``"GT_PRISM(n,k)"``
@@ -230,9 +230,9 @@ where ``"name of trans"`` can be chosen among the following 
list.
 
 * ``"GT_LINEAR_PRODUCT(a,b)"``
 
-  Description of the direct product of the two transformations ``a`` and ``b`` 
-  keeping a linear transformation (this is a restriction of the previous 
-  function). This allows, for instance, to use exact integrations on regular 
+  Description of the direct product of the two transformations ``a`` and ``b``
+  keeping a linear transformation (this is a restriction of the previous
+  function). This allows, for instance, to use exact integrations on regular
   meshes with parallelograms.
 
 
@@ -282,8 +282,8 @@ one has
 
 where :math:`\alpha` is the vector whose ith component is :math:`\alpha_i`.
 
-A certain number of description of classical finite element method are defined 
in 
-the file :file:`getfem_fem.h`. See :ref:`ud-appendixa` for an exhaustive list 
of 
+A certain number of description of classical finite element method are defined 
in
+the file :file:`getfem_fem.h`. See :ref:`ud-appendixa` for an exhaustive list 
of
 available finite element methods.
 
 A pointer to the finite element descriptor of a method is obtained using the
@@ -291,5 +291,5 @@ function::
 
   getfem::pfem pfe = getfem::fem_descriptor("name of method");
 
-We refer to the file :file:`getfem_fem.cc` for how to define a new finite 
element 
+We refer to the file :file:`getfem_fem.cc` for how to define a new finite 
element
 method.
diff --git a/doc/sphinx/source/project/intro.rst 
b/doc/sphinx/source/project/intro.rst
index 76fb2ac..7021a5b 100644
--- a/doc/sphinx/source/project/intro.rst
+++ b/doc/sphinx/source/project/intro.rst
@@ -21,7 +21,7 @@ on Savannah |linktask|_.
 
 
 The |gf| project focuses on the development of an open source generic
-finite element library. 
+finite element library.
 The goal is to provide a finite element framework which allows to
 easily build numerical code for the modelisation of system described
 by partial differential equations (p.d.e.). A special attention is paid
@@ -34,12 +34,12 @@ codes, is the complete separation between the description 
of p.d.e.
 models and finite element methods. Moreover, a separation is made
 between integration methods (exact or approximated), geometric
 transformations (linear or not) and finite element methods of
-arbitrary degrees described on a reference element. |gf| can 
+arbitrary degrees described on a reference element. |gf| can
 be used to build very general finite elements codes, where the
 finite elements, integration methods, dimension of the meshes,
-are just some parameters that can 
-be changed very easily, thus allowing a large spectrum of experimentations. 
-Numerous examples are available in the :file:`tests` directory of the 
+are just some parameters that can
+be changed very easily, thus allowing a large spectrum of experimentations.
+Numerous examples are available in the :file:`tests` directory of the
 distribution.
 
 The goal is also to make the addition of new finite element method
@@ -47,14 +47,14 @@ as simple as possible. For standard method, a description 
of the
 finite element shape functions and the type of connection of degrees
 of freedom on the reference element is sufficient. Extensions are
 provided for Hermite elements, piecewise polynomial, non-polynomial,
-vectorial elements and XFem. Examples of predefined 
-available methods are :math:`P_k` on simplices in arbitrary degrees and 
+vectorial elements and XFem. Examples of predefined
+available methods are :math:`P_k` on simplices in arbitrary degrees and
 dimensions, :math:`Q_k` on parallelepipeds, :math:`P_1`, :math:`P_2`
 with bubble functions, Hermite elements, elements with hierarchic
 basis (for multigrid methods for instance),
 discontinuous :math:`P_k` or :math:`Q_k`, XFem, Argyris, HCT, Raviart-Thomas.
 
-The library also includes the usual tools for finite elements such as assembly 
+The library also includes the usual tools for finite elements such as assembly
 procedures for classical PDEs, interpolation methods, computation of norms,
 mesh operations, boundary conditions, post-processing tools such as
 extraction of slices from a mesh ...
@@ -64,8 +64,8 @@ extraction of slices from a mesh ...
 The aim of the |gf| project is not to provide a ready to use
 finite element code allowing for instance structural mechanics
 computations with a graphic interface. It is basically a library
-allowing the build of C++  finite element codes. 
-However, the Python, Scilab and matlab interfaces allows to easily build 
application 
+allowing the build of C++  finite element codes.
+However, the Python, Scilab and matlab interfaces allows to easily build 
application
 coupling the definition of the problem, the finite element methods
 selection and the graphical post-processing.
 
diff --git a/doc/sphinx/source/project/libdesc.rst 
b/doc/sphinx/source/project/libdesc.rst
index b8570cc..d2cda4a 100644
--- a/doc/sphinx/source/project/libdesc.rst
+++ b/doc/sphinx/source/project/libdesc.rst
@@ -9,7 +9,7 @@
 Description of the different parts of the library
 =================================================
 
-Figure :ref:`dp-fig-diagram` describes the diagram of the different modules of 
+Figure :ref:`dp-fig-diagram` describes the diagram of the different modules of
 the |gf| library. The current state and perspective for each module is
 described in section :ref:`dp-libdesc`.
 
diff --git a/doc/sphinx/source/project/libdesc_cont.rst 
b/doc/sphinx/source/project/libdesc_cont.rst
index 1e62f32..8ec382c 100644
--- a/doc/sphinx/source/project/libdesc_cont.rst
+++ b/doc/sphinx/source/project/libdesc_cont.rst
@@ -24,12 +24,12 @@ Files
    :widths: 8, 15
 
    :file:`getfem_continuation.h` and :file:`getfem_continuation.cc`, "The 
generic continuation and branching method"
-   
+
 
 State
 ^^^^^
 
-Have already generic and advanced functionalities. 
+Have already generic and advanced functionalities.
 
 Perspectives
 ^^^^^^^^^^^^
diff --git a/doc/sphinx/source/project/libdesc_fem.rst 
b/doc/sphinx/source/project/libdesc_fem.rst
index 51fe177..5b9f2f9 100644
--- a/doc/sphinx/source/project/libdesc_fem.rst
+++ b/doc/sphinx/source/project/libdesc_fem.rst
@@ -36,7 +36,7 @@ Files
    :file:`getfem_interpolated_fem.h` and :file:`getfem_interpolated_fem.cc`, 
"Dfines a fem which is the interpolation of a finite element space (represented 
by a mesh_fem) on a different mesh. Note that the high-generic assembly 
language offers also this functionality by means of the interpolated 
transformations."
 
 
-   
+
 
 State
 ^^^^^
diff --git a/doc/sphinx/source/project/libdesc_high_gen_assemb.rst 
b/doc/sphinx/source/project/libdesc_high_gen_assemb.rst
index 35e1697..bd01d9a 100644
--- a/doc/sphinx/source/project/libdesc_high_gen_assemb.rst
+++ b/doc/sphinx/source/project/libdesc_high_gen_assemb.rst
@@ -29,7 +29,7 @@ Files
    :file:`getfem_generic_assembly_function_and_operators.h` and 
:file:`getfem_generic_assembly_function_and_operators.cc`, "Definition of 
redefined function and nonlinear operator of the weak form language."
    :file:`getfem_generic_assembly_semantic.h` and 
:file:`getfem_generic_assembly_semantic.cc`, "Semantic analysis and enrichment 
of the syntax tree. Include some operations such as making the derivation of a 
tree with respect to a variable or computing the tree corresponding to the 
gradient of an expression."
    :file:`getfem_generic_assembly_workspace.cc`, "Methodes of the workspace 
object (defined in :file:`getfem_generic_assembly.h`)."
-   :file:`getfem_generic_assembly_compile_and_exec.h` and 
:file:`getfem_generic_assembly_compile_and_exec.cc`, "Definition of the 
optimized instructions, compilation into a sequel of optimize instructions and 
execution of the instructions on Gauss point/interpolation points."   
+   :file:`getfem_generic_assembly_compile_and_exec.h` and 
:file:`getfem_generic_assembly_compile_and_exec.cc`, "Definition of the 
optimized instructions, compilation into a sequel of optimize instructions and 
execution of the instructions on Gauss point/interpolation points."
    :file:`getfem_generic_assembly_interpolation.cc`, "Interpolation operations 
and interpolate transformations."
 
 
@@ -86,12 +86,12 @@ Assembly tensors
 
 Assembly tensors are represented on each node by a ``bgeot::tensor<double>`` 
object. However, there is a specific structure in 
:file:`src/getfem\_generic\_assembly.cc` for assembly tensors because there is 
several format of assembly tensors :
 
-- Normal tensor. The first and second indices may represent the test function 
local indices if the node represent a first or second order term. Remember that 
in |gf| all tensors are stored with a Fortran order. This means that for 
instance a for a :math:`N\times P\times Q` tensor one has ``t(i, j, k) = t[i + 
j*N + k*N*P]``. 
+- Normal tensor. The first and second indices may represent the test function 
local indices if the node represent a first or second order term. Remember that 
in |gf| all tensors are stored with a Fortran order. This means that for 
instance a for a :math:`N\times P\times Q` tensor one has ``t(i, j, k) = t[i + 
j*N + k*N*P]``.
 
 - Copied tensor. When a node is detected to have exactly the same expression 
compared to an already compiled one, the assembly tensor will contain a pointer 
to the assembly tensor of the already compiled node. The consequence is that no 
unnecessary copy is made.
 
 - Sparse tensor with a listed sparsity. When working with a vector field, the 
finite element method is applied on each component. This results on vector base 
functions having only one nonzero component and some components are duplicated. 
The tensor are fully represented because it would be difficul to gain in 
efficiency with that kind of small sparse tensor format. However, some 
operation can be optimized with the knoledge of a certain sparsity (and 
duplication). This can change the orde [...]
-  
+
   - 1: Vectorized base sparsity format: The tensor represent a vectorized
     value. Each value of the condensed tensor is repeated on :math:`Q`
     components of the vectorized tensor. The mesh dimensions is denoted
@@ -110,7 +110,7 @@ Assembly tensors are represented on each node by a 
``bgeot::tensor<double>`` obj
 
        ":math:`\bar{t}(i) = \varphi_i(x)`", ":math:`t(j,k) = \varphi_{j/Q}(x) 
\delta_{k, (j \mbox{ mod } Q)}`"
        ":math:`[\varphi_0(x), \varphi_1(x)]`", ":math:`[\varphi_0(x), 0, 
\varphi_1(x), 0, 0, \varphi_0(x), 0, \varphi_1(x)]`"
-  
+
   - 2: Grad condensed format
 
     .. csv-table::
@@ -119,22 +119,22 @@ Assembly tensors are represented on each node by a 
``bgeot::tensor<double>`` obj
 
        ":math:`\bar{t}(i,j) = \partial_j\varphi_i(x)`", ":math:`t(k,l,m) = 
\partial_m\varphi_{k/Q}(x) \delta_{l, (m \mbox{ mod } Q)}`"
        ":math:`[\partial_0\varphi_0(x), \partial_0\varphi_1(x),` 
:math:`\partial_1\varphi_0(x), \partial_1\varphi_1(x),` 
:math:`\partial_2\varphi_0(x), \partial_2\varphi_1(x)]`",""
-  
+
   - 3: Hessian condensed format
 
 
   - 10: Vectorized mass: the tensor represent a scalar product of two
     vectorised base functions. This means a tensor :math:`t(\cdot,\cdot)`
     where :math:`t(i*Q+k, j*Q+l) = 0` for :math:`k \ne l` and
-    :math:`t(i*Q+k, j*Q+k)` are equals for :math:`0 \le k < Q`. 
-   
+    :math:`t(i*Q+k, j*Q+k)` are equals for :math:`0 \le k < Q`.
+
+
 
-    
 
 Optimized instructions
 ^^^^^^^^^^^^^^^^^^^^^^
 
-Optimized instructions for variable evaluation, operations, vector and matrix 
assembly ... to be described. 
+Optimized instructions for variable evaluation, operations, vector and matrix 
assembly ... to be described.
 
 Predefined functions
 ^^^^^^^^^^^^^^^^^^^^
diff --git a/doc/sphinx/source/project/libdesc_im.rst 
b/doc/sphinx/source/project/libdesc_im.rst
index b1087ee..1012bdf 100644
--- a/doc/sphinx/source/project/libdesc_im.rst
+++ b/doc/sphinx/source/project/libdesc_im.rst
@@ -37,14 +37,14 @@ Files
    :file:`getfem_im_list.h`, "file generated by 
:file:`cubature/make_getfem_list` with the integration methods defined in 
cubature directory. This gives the possibility to define a new integration 
method just listing the Gauss points and weigth in a text file."
 
 
-State 
+State
 ^^^^^
 
 This module meets the present needs for the project and is considered as
-stabilized. The list of available cubature formulas is given in 
+stabilized. The list of available cubature formulas is given in
 :ref:`ud-appendixb`.
 
-Perspectives 
+Perspectives
 ^^^^^^^^^^^^
 
 No change needed for the moment. An effort could be done on the documentation 
to
diff --git a/doc/sphinx/source/project/libdesc_interface.rst 
b/doc/sphinx/source/project/libdesc_interface.rst
index 426366c..29a4c42 100644
--- a/doc/sphinx/source/project/libdesc_interface.rst
+++ b/doc/sphinx/source/project/libdesc_interface.rst
@@ -14,7 +14,7 @@ A simplified (but rather complete) interface of |gf| is 
provided, so that it is
 
 Description
 ^^^^^^^^^^^
- 
+
 All sources are located in the :file:`interface/src` directory. The interface 
is
 composed of one large library ``getfemint`` (which stands for getfem
 interaction), which acts as a layer above the |gf| library, and is used by
@@ -33,75 +33,75 @@ All the files in the directory :file:`interface\src`. A 
short description of mai
 
 * :file:`getfem_interface.cc`.
 
-  This is the bridge between the script language and the getfem interface. The 
-  function getfem_interface_main is exported as an ``extern "C"`` function, so 
-  this is a sort of c++ barrier between the script language and the getfem 
+  This is the bridge between the script language and the getfem interface. The
+  function getfem_interface_main is exported as an ``extern "C"`` function, so
+  this is a sort of c++ barrier between the script language and the getfem
   interface (exporting only a C interface avoids many compilation problems).
 
 * :file:`matlab/gfm_mex.c`.
 
-  The matlab interface. The only thing it knows about getfem is in 
+  The matlab interface. The only thing it knows about getfem is in
   :file:`getfem_interface.h`.
 
 * :file:`python/getfem_python.c`.
 
-  The python interface. The only thing it knows about getfem is in 
+  The python interface. The only thing it knows about getfem is in
   :file:`getfem_interface.h`.
 
 * :file:`gfi_array.h`, :file:`gfi_array.c`.
 
-  Both :file:`gfm_mex.c` and :file:`getfem_python.c` need a simple convention 
on 
-  how to send and receive arrays, and object handles, from 
+  Both :file:`gfm_mex.c` and :file:`getfem_python.c` need a simple convention 
on
+  how to send and receive arrays, and object handles, from
   ``getfem_interface_main()``. This file provide such functionnality.
 
 * :file:`getfemint_gsparse.h`, :file:`getfemint_precond.h`, etc.
 
-  Files specific to an interfaced object if needed. 
-  (getfemint_gsparse which export some kind of mutable sparse matrix that can 
+  Files specific to an interfaced object if needed.
+  (getfemint_gsparse which export some kind of mutable sparse matrix that can
   switch between different storage types, and real of complex elements).
 
 * :file:`gf_workspace.cc`, :file:`gf_delete.cc`.
 
   Memory management for getfem objects. There is a layer which handles the
   dependency between for example a ``mesh`` and a ``mesh_fem``.
-  It makes sure that no object 
+  It makes sure that no object
   will be destroyed while there is still another getfem_object using it.
-  The goal 
-  is to make sure that under no circumstances the user is able to crash getfem 
+  The goal
+  is to make sure that under no circumstances the user is able to crash getfem
   (and the host program, matlab, scilab or python) by passing incorrect
   argument to the getfem interface.
 
-  It also provides a kind of workspace stack, which was designed to simplify 
-  handling and cleaning of many getfem objects in matlab (since matlab does 
not 
+  It also provides a kind of workspace stack, which was designed to simplify
+  handling and cleaning of many getfem objects in matlab (since matlab does not
   have "object destructors").
 
 * :file:`getfemint.h`, :file:`getfemint.cc`.
 
-  Define the ``mexarg_in``, ``mexarg_out`` classes, which are used to parse 
the 
+  Define the ``mexarg_in``, ``mexarg_out`` classes, which are used to parse the
   list of input and output arguments to the getfem interface functions.
   The name  is not adequate anymore since any reference to "mex"
-  has been moved into 
+  has been moved into
   :file:`gfm_mex.c`.
 
 * :file:`gf_mesh.cc`, :file:`gf_mesh_get.cc`, :file:`gf_mesh_set.cc`,
   :file:`gf_fem.cc`, etc.
 
-  All the functions exported be the getfem interfaces, sorted by object type 
-  (``gf_mesh*``, ``gf_mesh_fem*``, ``gf_fem*``), and then organized as one for 
-  the object construction (``gf_mesh``), one for the object modification 
-  (``gf_mesh_set``), and one for the object inquiry (``gf_mesh_get``). Each of 
-  these files contain one main function, that receives a ``mexargs_in`` and 
-  ``mexargs_out`` stack of arguments. It parses then, and usually interprets 
the 
-  first argument as the name of a subfunction (``gf_mesh_get('nbpts')`` in 
+  All the functions exported be the getfem interfaces, sorted by object type
+  (``gf_mesh*``, ``gf_mesh_fem*``, ``gf_fem*``), and then organized as one for
+  the object construction (``gf_mesh``), one for the object modification
+  (``gf_mesh_set``), and one for the object inquiry (``gf_mesh_get``). Each of
+  these files contain one main function, that receives a ``mexargs_in`` and
+  ``mexargs_out`` stack of arguments. It parses then, and usually interprets 
the
+  first argument as the name of a subfunction (``gf_mesh_get('nbpts')`` in
   matlab, or ``Mesh.nbpts()`` in python).
 
 * :file:`matlab/gfm_rpx_mexint.c`.
 
-  An alternative to :file:`gfm_mex.c` which is used when the 
-  ``--enable-matlab-rpc`` is passed to the ``./configure`` script. The main 
use 
-  for that is debugging the interface, since in that case, the matlab 
interface 
-  communicates via sockets with a "getfem_server" program, so it is possible 
to 
-  debug that server program, and identify memory leaks or anything else 
without 
+  An alternative to :file:`gfm_mex.c` which is used when the
+  ``--enable-matlab-rpc`` is passed to the ``./configure`` script. The main use
+  for that is debugging the interface, since in that case, the matlab interface
+  communicates via sockets with a "getfem_server" program, so it is possible to
+  debug that server program, and identify memory leaks or anything else without
   having to mess with matlab (it is pain to debug).
 
 * :file:`python/getfem.py`.
@@ -112,7 +112,7 @@ All the files in the directory :file:`interface\src`. A 
short description of mai
 
 
 
-Objects, methods and functions of the interface 
+Objects, methods and functions of the interface
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
 The main concepts manipulated by the interface are a limited number of objects
@@ -217,7 +217,7 @@ To make all the things work automatically, a certain number 
of rules have to be
         - address@hidden  for  address@hidden
 
 
-  Three dots at the end of the parameter list (``...``) mean that 
+  Three dots at the end of the parameter list (``...``) mean that
   additional parameters are possible. Optional parameters can be described
   with brackets. For instance ``/address@hidden v = ('name'[, @int i])``. But
   be careful how it is interpreted by the :file:`extract_doc` script
@@ -236,12 +236,12 @@ To make all the things work automatically, a certain 
number of rules have to be
 
 .. _reStructuredText: http://docutils.sourceforge.net/rst.html
 
-* The parts of documentation included in the c++ commentaries should be in 
+* The parts of documentation included in the c++ commentaries should be in
   `reStructuredText`_ format. In particular, math formulas can be included
   with \:math\:\`f(x) = 3x^2+2x+4\` or with::
-  
+
     .. math::
- 
+
       f(x) = 3x^2+2x+4
 
   It is possible to refer to another method or function of the interface
@@ -298,7 +298,7 @@ For the management of the object, you have to declare the 
class at the begining
   id_type store_"name"_object(const std::shared_ptr<object_class> &shp);
   object_class *to_"name"_object(const mexarg_in &p);
 
-where "name" is the name of the object in the interface and ``object_class`` 
is the class name in getfem (for instance  ``getfem::mesh`` for the mesh 
object). Alternatively, for the object that are manipulated by a shared pointer 
in |gf|, the third function can return a shared pointer. 
+where "name" is the name of the object in the interface and ``object_class`` 
is the class name in getfem (for instance  ``getfem::mesh`` for the mesh 
object). Alternatively, for the object that are manipulated by a shared pointer 
in |gf|, the third function can return a shared pointer.
 
 IMPORTANT: In order to be interfaced, a |gf| object has to derive from 
``dal::static_stored_object``. However, if it is not the case, a wrapper class 
can be defined such as the one for ``bgeot::base_poly`` (see the end of 
:file:`getfemint.h`).
 
diff --git a/doc/sphinx/source/project/libdesc_mesh.rst 
b/doc/sphinx/source/project/libdesc_mesh.rst
index d0525eb..778b782 100644
--- a/doc/sphinx/source/project/libdesc_mesh.rst
+++ b/doc/sphinx/source/project/libdesc_mesh.rst
@@ -76,7 +76,7 @@ namespace (``getfem``).
 
 A  bibliographical review on how to efficiently store a mesh and implement the 
main operations (add a node, an element, deal with faces, find the neighbour 
elements, the isolated faces ...) would be interesting to make the mesh 
structure evolve.
 
-A senstive algorithm is the one (in bgeot_node_tab.cc) which identify the too 
much close nodes. More investigations (and documentation) are probably 
necessary. 
+A senstive algorithm is the one (in bgeot_node_tab.cc) which identify the too 
much close nodes. More investigations (and documentation) are probably 
necessary.
 
 
 
diff --git a/doc/sphinx/source/project/libdesc_model.rst 
b/doc/sphinx/source/project/libdesc_model.rst
index 7910665..162b47a 100644
--- a/doc/sphinx/source/project/libdesc_model.rst
+++ b/doc/sphinx/source/project/libdesc_model.rst
@@ -22,7 +22,7 @@ Files
    :header: "File(s)", "Description"
    :widths: 8, 15
 
-   :file:`getfem_models.h` and :file:`getfem_models.cc`, "Defines the object 
models, its internal and the standard bricks (linear elasticity, generic 
elliptic brick, Dirichlet boundary conditions ...)." 
+   :file:`getfem_models.h` and :file:`getfem_models.cc`, "Defines the object 
models, its internal and the standard bricks (linear elasticity, generic 
elliptic brick, Dirichlet boundary conditions ...)."
    :file:`getfem_model_solvers.h` and :file:`getfem_model_solvers.cc`, 
"Defines the the standard solvers for the model object."
    :file:`getfem_contact_and_friction_common.h` and 
:file:`getfem_contact_and_friction_common.cc`, "Common algorithms for 
contact/friction conditions on deformable bodies"
    :file:`getfem_contact_and_friction_integral.h` and 
:file:`getfem_contact_and_friction_integral.cc`, "Small sliding 
Contact/friction bricks of integral type."
diff --git a/doc/sphinx/source/python/examples.rst 
b/doc/sphinx/source/python/examples.rst
index a971d92..a600021 100644
--- a/doc/sphinx/source/python/examples.rst
+++ b/doc/sphinx/source/python/examples.rst
@@ -16,84 +16,84 @@ Examples
 A step-by-step basic example
 ----------------------------
 
-This example shows the basic usage of getfem, on the über-canonical problem 
above 
-all others: solving the :envvar:`Laplacian`, :math:`-\Delta u = f` on a 
square, 
-with the Dirichlet condition :math:`u = g(x)` on the domain boundary. You can 
find 
-the **py-file** of this example under the name **demo_step_by_step.py** in the 
+This example shows the basic usage of getfem, on the über-canonical problem 
above
+all others: solving the :envvar:`Laplacian`, :math:`-\Delta u = f` on a square,
+with the Dirichlet condition :math:`u = g(x)` on the domain boundary. You can 
find
+the **py-file** of this example under the name **demo_step_by_step.py** in the
 directory ``interface/tests/python/`` of the |gf| distribution.
 
-The first step is to **create a Mesh object**. It is possible to create simple 
structured meshes or unstructured meshes for simple geometries (see 
``getfem.Mesh('generate', mesher_object mo, scalar h))``) or to rely on an 
external mesher (see ``getfem.Mesh('import', 
-string FORMAT, string FILENAME)``), or use very simple meshes. For this 
example, 
-we just consider a regular mesh\index{cartesian mesh} whose nodes are 
+The first step is to **create a Mesh object**. It is possible to create simple 
structured meshes or unstructured meshes for simple geometries (see 
``getfem.Mesh('generate', mesher_object mo, scalar h))``) or to rely on an 
external mesher (see ``getfem.Mesh('import',
+string FORMAT, string FILENAME)``), or use very simple meshes. For this 
example,
+we just consider a regular mesh\index{cartesian mesh} whose nodes are
 :math:`\{x_{i=0\ldots10,j=0..10}=(i/10,j/10)\}`
 
 .. literalinclude:: code_samples/demo_step_by_step.py
    :linenos:
    :lines: 4-9
 
-The next step is to **create a MeshFem object**. This one links a mesh with a 
set 
+The next step is to **create a MeshFem object**. This one links a mesh with a 
set
 of FEM
 
 .. literalinclude:: code_samples/demo_step_by_step.py
    :linenos:
    :lines: 11-14
 
-The first instruction builds a new |py_mf| object, the second argument 
specifies 
-that this object will be used to interpolate scalar fields (since the unknown 
-:math:`u` is a scalar field). The second instruction assigns the :math:`Q^2` 
FEM 
-to every convex (each basis function is a polynomial of degree 4, remember 
that 
-:math:`P^k\Rightarrow` polynomials of degree :math:`k`, while 
-:math:`Q^k\Rightarrow` polynomials of degree :math:`2k`). As :math:`Q^2` is a 
-polynomial FEM, you can view the expression of its basis functions on the 
+The first instruction builds a new |py_mf| object, the second argument 
specifies
+that this object will be used to interpolate scalar fields (since the unknown
+:math:`u` is a scalar field). The second instruction assigns the :math:`Q^2` 
FEM
+to every convex (each basis function is a polynomial of degree 4, remember that
+:math:`P^k\Rightarrow` polynomials of degree :math:`k`, while
+:math:`Q^k\Rightarrow` polynomials of degree :math:`2k`). As :math:`Q^2` is a
+polynomial FEM, you can view the expression of its basis functions on the
 reference convex:
 
 .. literalinclude:: code_samples/demo_step_by_step.py
    :linenos:
    :lines: 16-17
 
-Now, in order to perform numerical integrations on ``mf``, we need to **build 
a 
+Now, in order to perform numerical integrations on ``mf``, we need to **build a
 MeshIm object**
 
 .. literalinclude:: code_samples/demo_step_by_step.py
    :linenos:
    :lines: 19-20
 
-The integration method will be used to compute the various integrals on each 
-element: here we choose to perform exact computations (no :envvar:`quadrature 
-formula`), which is possible since the geometric transformation of these 
convexes 
-from the reference convex is linear (this is true for all simplices, and this 
is 
-also true for the parallelepipeds of our regular mesh, but it is not true for 
-general quadrangles), and the chosen FEM is polynomial. Hence it is possible 
to 
-analytically integrate every basis function/product of basis 
-functions/gradients/etc. There are many alternative FEM methods and 
integration 
+The integration method will be used to compute the various integrals on each
+element: here we choose to perform exact computations (no :envvar:`quadrature
+formula`), which is possible since the geometric transformation of these 
convexes
+from the reference convex is linear (this is true for all simplices, and this 
is
+also true for the parallelepipeds of our regular mesh, but it is not true for
+general quadrangles), and the chosen FEM is polynomial. Hence it is possible to
+analytically integrate every basis function/product of basis
+functions/gradients/etc. There are many alternative FEM methods and integration
 methods (see :ref:`ud`).
 
-Note however that in the general case, approximate integration methods are a 
+Note however that in the general case, approximate integration methods are a
 better choice than exact integration methods.
 
-Now we have to **find the** <:envvar:`boundary`> **of the domain**, in order 
to 
-set a Dirichlet condition. A mesh object has the ability to store some sets of 
-convexes and convex faces. These sets (called <regions>) are accessed via an 
+Now we have to **find the** <:envvar:`boundary`> **of the domain**, in order to
+set a Dirichlet condition. A mesh object has the ability to store some sets of
+convexes and convex faces. These sets (called <regions>) are accessed via an
 integer *#id*
 
 .. literalinclude:: code_samples/demo_step_by_step.py
    :linenos:
    :lines: 22-25
 
-Here we find the faces of the convexes which are on the boundary of the mesh 
(i.e. 
+Here we find the faces of the convexes which are on the boundary of the mesh 
(i.e.
 the faces which are not shared by two convexes).
 
-The array ``border`` has two rows, on the first row is a convex number, on the 
-second row is a face number (which is local to the convex, there is no global 
+The array ``border`` has two rows, on the first row is a convex number, on the
+second row is a face number (which is local to the convex, there is no global
 numbering of faces). Then this set of faces is assigned to the region number 
42.
 
-At this point, we just have to describe the model and run the solver to get 
the 
-solution! The ":envvar:`model`" is created with the |py_md| constructor. A 
model 
-is basically an object which build a global linear system (tangent matrix for 
-non-linear problems) and its associated right hand side. Typical modifications 
are 
-insertion of the stiffness matrix for the problem considered (linear 
elasticity, 
-laplacian, etc), handling of a set of constraints, Dirichlet condition, 
addition of 
-a source term to the right hand side etc. The global tangent matrix and its 
right 
+At this point, we just have to describe the model and run the solver to get the
+solution! The ":envvar:`model`" is created with the |py_md| constructor. A 
model
+is basically an object which build a global linear system (tangent matrix for
+non-linear problems) and its associated right hand side. Typical modifications 
are
+insertion of the stiffness matrix for the problem considered (linear 
elasticity,
+laplacian, etc), handling of a set of constraints, Dirichlet condition, 
addition of
+a source term to the right hand side etc. The global tangent matrix and its 
right
 hand side are stored in the ":envvar:`model`" structure.
 
 Let us build a problem with an easy solution: :math:`u = x(x-1)-y(y-1)`, then
@@ -106,16 +106,16 @@ We start with an empty real model
    :linenos:
    :lines: 27-28
 
-(a model is either ``'real'`` or ``'complex'``). And we declare that ``u`` is 
an 
+(a model is either ``'real'`` or ``'complex'``). And we declare that ``u`` is 
an
 unknown of the system on the finite element method `mf` by
 
 .. literalinclude:: code_samples/demo_step_by_step.py
    :linenos:
    :lines: 30-32
 
-Now, we add a `generic elliptic` brick, which handles 
:math:`-\nabla\cdot(A:\nabla 
-u) = \ldots` problems, where :math:`A` can be a scalar field, a matrix field, 
or 
-an order 4 tensor field. By default, :math:`A=1`. We add it on our main 
variable 
+Now, we add a `generic elliptic` brick, which handles 
:math:`-\nabla\cdot(A:\nabla
+u) = \ldots` problems, where :math:`A` can be a scalar field, a matrix field, 
or
+an order 4 tensor field. By default, :math:`A=1`. We add it on our main 
variable
 ``u`` with
 
 .. literalinclude:: code_samples/demo_step_by_step.py
@@ -128,20 +128,20 @@ Next we add a Dirichlet condition on the domain boundary
    :linenos:
    :lines: 37-40
 
-The two first lines defines a data of the model which represents the value of 
the 
-Dirichlet condition. The third one add a Dirichlet condition to the variable 
``u`` 
-on the boundary number ``42``. The dirichlet condition is imposed with 
lagrange 
-multipliers. Another possibility is to use a penalization. A |py_mf| argument 
is 
-also required, as the Dirichlet condition :math:`u=g` is imposed in a weak 
form 
-:math:`\int_\Gamma u(x)v(x) = \int_\Gamma g(x)v(x)\ \forall v` where :math:`v` 
is 
+The two first lines defines a data of the model which represents the value of 
the
+Dirichlet condition. The third one add a Dirichlet condition to the variable 
``u``
+on the boundary number ``42``. The dirichlet condition is imposed with lagrange
+multipliers. Another possibility is to use a penalization. A |py_mf| argument 
is
+also required, as the Dirichlet condition :math:`u=g` is imposed in a weak form
+:math:`\int_\Gamma u(x)v(x) = \int_\Gamma g(x)v(x)\ \forall v` where :math:`v` 
is
 taken in the space of multipliers given by here by ``mf``.
 
 .. topic:: Remark:
 
-   the polynomial expression was interpolated on ``mf``. It is possible only 
if 
-   ``mf`` is of Lagrange type. In this first example we use the same |py_mf| 
for 
-   the unknown and for the data such as ``g``, but in the general case, ``mf`` 
-   won't be Lagrangian and another (Lagrangian) |py_mf| will be used for the 
+   the polynomial expression was interpolated on ``mf``. It is possible only if
+   ``mf`` is of Lagrange type. In this first example we use the same |py_mf| 
for
+   the unknown and for the data such as ``g``, but in the general case, ``mf``
+   won't be Lagrangian and another (Lagrangian) |py_mf| will be used for the
    description of Dirichlet conditions, source terms etc.
 
 A source term can be added with (uncommented) the following lines
@@ -150,14 +150,14 @@ A source term can be added with (uncommented) the 
following lines
    :linenos:
    :lines: 42-45
 
-It only remains now to launch the solver. The linear system is assembled and 
solve 
+It only remains now to launch the solver. The linear system is assembled and 
solve
 with the instruction
 
 .. literalinclude:: code_samples/demo_step_by_step.py
    :linenos:
    :lines: 47-48
 
-The model now contains the solution (as well as other things, such as the 
linear 
+The model now contains the solution (as well as other things, such as the 
linear
 system which was solved). It is extracted
 
 .. literalinclude:: code_samples/demo_step_by_step.py
@@ -183,12 +183,12 @@ and view with ``gmsh u.pos``, see figure 
:ref:`py-fig-sbs`.
 Another Laplacian with exact solution (source term)
 ---------------------------------------------------
 
-This example shows the basic usage of getfem, on the canonical problem: 
solving 
-the Laplacian, :math:`-\Delta u = f` on a square, with the Dirichlet condition 
-:math:`u = g(x)` on the domain boundary :math:`\Gamma_D` and the Neumann 
condition 
-:math:`\frac{\partial u}{\partial\eta} = h(x)` on the domain boundary 
-:math:`\Gamma_N`. You can find the **py-file** of this example under the name 
-**demo_laplacian.py** in the directory ``interface/tests/python/`` of the |gf| 
+This example shows the basic usage of getfem, on the canonical problem: solving
+the Laplacian, :math:`-\Delta u = f` on a square, with the Dirichlet condition
+:math:`u = g(x)` on the domain boundary :math:`\Gamma_D` and the Neumann 
condition
+:math:`\frac{\partial u}{\partial\eta} = h(x)` on the domain boundary
+:math:`\Gamma_N`. You can find the **py-file** of this example under the name
+**demo_laplacian.py** in the directory ``interface/tests/python/`` of the |gf|
 distribution.
 
 We create Mesh, MeshFem, MeshIm object and find the boundary of the domain in
@@ -210,7 +210,7 @@ and we bricked the problem as in the previous example
    :linenos:
    :lines: 77-102
 
-the only change is the add of `source term` bricks. Finally the solution of 
the 
+the only change is the add of `source term` bricks. Finally the solution of the
 problem is extracted and exported
 
 .. literalinclude:: code_samples/demo_laplacian.py
@@ -229,16 +229,16 @@ view with ``gmsh sol.pos``:
 Linear and non-linear elasticity
 --------------------------------
 
-This example uses a mesh that was generated with `GiD`_. The object is meshed 
-with quadratic tetrahedrons. You can find the **py-file** of this example 
under 
-the name :file:`demo_tripod.py` in the directory 
:file:`interface/tests/python/` 
+This example uses a mesh that was generated with `GiD`_. The object is meshed
+with quadratic tetrahedrons. You can find the **py-file** of this example under
+the name :file:`demo_tripod.py` in the directory 
:file:`interface/tests/python/`
 of the |gf| distribution.
 
 .. literalinclude:: code_samples/demo_tripod.py
    :linenos:
    :lines: 3-89
 
-Here is the final figure, displaying the :envvar:`Von Mises` stress and 
+Here is the final figure, displaying the :envvar:`Von Mises` stress and
 displacements norms:
 
 .. figure:: images/tripod.png
@@ -251,52 +251,52 @@ displacements norms:
 Avoiding the model framework
 ----------------------------
 
-The model bricks are very convenient, as they hide most of the details of the 
-assembly of the final linear systems. However it is also possible to stay at a 
-lower level, and handle the assembly of linear systems, and their resolution, 
-directly in |py|. For example, the demonstration :file:`demo_tripod_alt.py` is 
+The model bricks are very convenient, as they hide most of the details of the
+assembly of the final linear systems. However it is also possible to stay at a
+lower level, and handle the assembly of linear systems, and their resolution,
+directly in |py|. For example, the demonstration :file:`demo_tripod_alt.py` is
 very similar to the :file:`demo_tripod.py` except that the assembly is explicit
 
 .. literalinclude:: code_samples/demo_tripod_alt.py
    :lines: 21,23,25,27,49-51,53,58,62-64,70,73-81,83,85-99,113,118-
 
-In |gfi|, the assembly of vectors, and matrices is done via the ``gf.asm_*`` 
-functions. The Dirichlet condition :math:`h(x)u(x) = r(x)` is handled in the 
-weak form :math:`\int (h(x)u(x)).v(x) = \int r(x).v(x)\quad\forall v` (where 
-:math:`h(x)` is a :math:`3\times 3` matrix field -- here it is constant and 
-equal to the identity). The reduced system ``KK UU = FF`` is then built via 
the 
-elimination of Dirichlet constraints from the original system. Note that it 
-might be more efficient (and simpler) to deal with Dirichlet condition via a 
+In |gfi|, the assembly of vectors, and matrices is done via the ``gf.asm_*``
+functions. The Dirichlet condition :math:`h(x)u(x) = r(x)` is handled in the
+weak form :math:`\int (h(x)u(x)).v(x) = \int r(x).v(x)\quad\forall v` (where
+:math:`h(x)` is a :math:`3\times 3` matrix field -- here it is constant and
+equal to the identity). The reduced system ``KK UU = FF`` is then built via the
+elimination of Dirichlet constraints from the original system. Note that it
+might be more efficient (and simpler) to deal with Dirichlet condition via a
 penalization technique.
 
 Other examples
 --------------
 
-* the :file:`demo_refine.py` script shows a simple 2D or 3D bar whose 
extremity 
-  is clamped. An adaptative refinement is used to obtain a better 
approximation 
-  in the area where the stress is singular (the transition between the clamped 
+* the :file:`demo_refine.py` script shows a simple 2D or 3D bar whose extremity
+  is clamped. An adaptative refinement is used to obtain a better approximation
+  in the area where the stress is singular (the transition between the clamped
   area and the neumann boundary).
 
-* the :file:`demo_nonlinear_elasticity.py` script shows a 3D bar which is is 
-  bended and twisted. This is a quasi-static problem as the deformation is 
-  applied in many steps. At each step, a non-linear (large deformations) 
+* the :file:`demo_nonlinear_elasticity.py` script shows a 3D bar which is is
+  bended and twisted. This is a quasi-static problem as the deformation is
+  applied in many steps. At each step, a non-linear (large deformations)
   elasticity problem is solved.
 
-* the :file:`demo_stokes_3D_tank.py` script shows a Stokes (viscous fluid) 
-  problem in a tank. The :file:`demo_stokes_3D_tank_draw.py` shows how to draw 
-  a nice plot of the solution, with mesh slices and stream lines. Note that 
the 
-  :file:`demo_stokes_3D_tank_alt.py` is the old example, which uses the 
+* the :file:`demo_stokes_3D_tank.py` script shows a Stokes (viscous fluid)
+  problem in a tank. The :file:`demo_stokes_3D_tank_draw.py` shows how to draw
+  a nice plot of the solution, with mesh slices and stream lines. Note that the
+  :file:`demo_stokes_3D_tank_alt.py` is the old example, which uses the
   deprecated ``gf_solve`` function.
 
-* the :file:`demo_bilaplacian.py` script is just an adaption of the |gf| 
-  example :file:`tests/bilaplacian.cc`. Solve the bilaplacian (or a 
+* the :file:`demo_bilaplacian.py` script is just an adaption of the |gf|
+  example :file:`tests/bilaplacian.cc`. Solve the bilaplacian (or a
   Kirchhoff-Love plate model) on a square.
 
-* the :file:`demo_plasticity.py` script is an adaptation of the |gf| example 
-  :file:`tests/plasticity.cc`: a 2D or 3D bar is bended in many steps, and the 
-  plasticity of the material is taken into account (plastification occurs when 
+* the :file:`demo_plasticity.py` script is an adaptation of the |gf| example
+  :file:`tests/plasticity.cc`: a 2D or 3D bar is bended in many steps, and the
+  plasticity of the material is taken into account (plastification occurs when
   the material's Von Mises exceeds a given threshold).
 
-* the :file:`demo_wave2D.py` is a 2D scalar wave equation example (diffraction 
-  of a plane wave by a cylinder), with high order geometric transformations 
and 
+* the :file:`demo_wave2D.py` is a 2D scalar wave equation example (diffraction
+  of a plane wave by a cylinder), with high order geometric transformations and
   high order FEMs.
diff --git a/doc/sphinx/source/python/howtos.rst 
b/doc/sphinx/source/python/howtos.rst
index fe56984..a69bba3 100644
--- a/doc/sphinx/source/python/howtos.rst
+++ b/doc/sphinx/source/python/howtos.rst
@@ -29,7 +29,7 @@ regions. We can import that file (*quad.msh*) to getfem::
 with the second command we can see the *regions ids*. When we import the mesh,
 we might be warned with the following::
 
-  Level 3 Warning in getfem_import.cc, line 137: 
+  Level 3 Warning in getfem_import.cc, line 137:
     All regions must have different number!
 
 this means that the parametrization of the mesh in |gmsh| *.geo file* must
diff --git a/doc/sphinx/source/python/intro.rst 
b/doc/sphinx/source/python/intro.rst
index 78ff421..766307f 100644
--- a/doc/sphinx/source/python/intro.rst
+++ b/doc/sphinx/source/python/intro.rst
@@ -9,9 +9,9 @@
 Introduction
 ============
 
-This guide provides a reference about the |py| interface of |gf|. For a 
complete 
-reference of |gf|, please report to the `specific guides`_, but you should be 
able 
-to use the |gfi|'s without any particular knowledge of the |gf| internals, 
+This guide provides a reference about the |py| interface of |gf|. For a 
complete
+reference of |gf|, please report to the `specific guides`_, but you should be 
able
+to use the |gfi|'s without any particular knowledge of the |gf| internals,
 although a basic knowledge about Finite Elements is required. This 
documentation
 is however not self contained. You should in
 particular refer to the `user documentation`_ to have a more extensive
diff --git a/doc/sphinx/source/scilab/install.rst 
b/doc/sphinx/source/scilab/install.rst
index fb3c9bc..d742add 100644
--- a/doc/sphinx/source/scilab/install.rst
+++ b/doc/sphinx/source/scilab/install.rst
@@ -6,7 +6,7 @@
 
 .. _sci-install:
 
-Installation 
+Installation
 ============
 
 The installation of the |sci| |gf| toolbox can be somewhat tricky, since it 
combines a C++ compiler, libraries and |sci| interaction. In case of troubles 
with a
diff --git a/doc/sphinx/source/scilab/intro.rst 
b/doc/sphinx/source/scilab/intro.rst
index f9338fb..e3fa55f 100644
--- a/doc/sphinx/source/scilab/intro.rst
+++ b/doc/sphinx/source/scilab/intro.rst
@@ -4,9 +4,9 @@
 Introduction
 ============
 
-This guide provides a reference about the |sci| interface of |gf|. For a 
complete 
-reference of |gf|, please report to the `specific guides`_, but you should be 
-able to use the scilab interface without any particular knowledge of the |gf| 
+This guide provides a reference about the |sci| interface of |gf|. For a 
complete
+reference of |gf|, please report to the `specific guides`_, but you should be
+able to use the scilab interface without any particular knowledge of the |gf|
 internals, although a basic knowledge about Finite Elements is required.
 This documentation is however not self contained. You should in
 particular refer to the `user documentation`_ to have a more extensive
diff --git a/doc/sphinx/source/scilab/plotcmdref.rst 
b/doc/sphinx/source/scilab/plotcmdref.rst
index 563593a..173c0b2 100644
--- a/doc/sphinx/source/scilab/plotcmdref.rst
+++ b/doc/sphinx/source/scilab/plotcmdref.rst
@@ -36,19 +36,19 @@ gf_plot
 
   The options are specified as pairs of "option name"/"option value"
 
-  'zplot',{'off'|'on'}       : values of ``U`` are mapped on the $z$-axis 
(only possible when qdim=1, mdim=2). 
+  'zplot',{'off'|'on'}       : values of ``U`` are mapped on the $z$-axis 
(only possible when qdim=1, mdim=2).
   'norm', {'off'|'on'}       : if qdim >= 2, color-plot the norm of the field
   'dir',[]                   : or the scalar product of the field with 'dir' 
(can be a vector, or 'x', 'y' etc..)
   'refine',8                 : nb of refinments for curved edges and surface 
plots
   'interpolated',{'off'|'on'}: if triangular patch are interpolated
   'pcolor',{'on'|'off'}      : if the field is scalar, a color plot of its 
values is plotted
-  'quiver',{'on'|'off'}      : if the field is vector, represent arrows        
       
+  'quiver',{'on'|'off'}      : if the field is vector, represent arrows        
   'quiver_density',50        : density of arrows in quiver plot
   'quiver_scale',1           : scaling of arrows (0=>no scaling)
   'mesh',{'off'|'on'}        : show the mesh ?
-  'meshopts',{cell(0)}          : cell array of options passed to 
gf_plot_slice for the mesh 
+  'meshopts',{cell(0)}          : cell array of options passed to 
gf_plot_slice for the mesh
   'deformed_mesh', {'off'|'on'} : shows the deformed mesh (only when qdim == 
mdim)
-  'deformed_meshopts', {cell(0)}: cell array of options passed to 
gf_plot_slice for the deformed mesh 
+  'deformed_meshopts', {cell(0)}: cell array of options passed to 
gf_plot_slice for the deformed mesh
   'deformation',[]           : plots on the deformed object (only when qdim == 
mdim)
   'deformation_mf',[]        : plots on the deformed object (only when qdim == 
mdim)
   'deformation_scale','10%'  : indicate the amplitude of the deformation. Can 
be a percentage of the mesh width if given as a string, or an absolute value if 
given as a number
@@ -76,14 +76,14 @@ gf_plot
   dof) labels.
 
   For example, plotting a scalar field on the border of a 3D mesh can be done 
with ::
-  
+
     % load the 'strange.mesh_fem' (found in the getfem_matlab/tests directory)
-    mf=gf_mesh_fem('load', 'strange.mesh_fem') 
+    mf=gf_mesh_fem('load', 'strange.mesh_fem')
     U=rand(1, gf_mesh_fem_get(mf, 'nbdof')); # random field that will be drawn
     gf_plot(mf, U, 'refine', 25, 'cvlst', gf_mesh_get(mf,'outer faces'), 
'mesh','on');
- 
 
- 
+
+
 
 gf_plot_1D
 -------------------------------------------
@@ -117,18 +117,18 @@ gf_plot_mesh
 
   gf_plot_mesh(m, ...)
 
-  'vertices', {'off'|'on'}    : displays also vertices numbers. 
-  'convexes', {'off'|'on'}    : displays also convexes numbers. 
+  'vertices', {'off'|'on'}    : displays also vertices numbers.
+  'convexes', {'off'|'on'}    : displays also convexes numbers.
   'dof',{'off'|'on'}          : displays also finite element nodes. In that 
case, ``m`` should be a ``mesh_fem`` identifier.
   'regions',BLST              : displays the boundaries listed in BLST.
   'cvlst',CVLST               : display only the listed convexes. If CVLST has 
two rows, display only the faces listed in the second row.
   'edges', {'on' | 'off'}     : display edges ?
   'faces', {'off'|'on'}       : fills each 2D-face of the mesh
   'curved', {'off'|'on'}      : displays curved edges
-  'refine',N                  : refine curved edges and filled faces N times  
+  'refine',N                  : refine curved edges and filled faces N times
   'deformation', Udef         : optionnal deformation applied to the mesh (M 
must be a mesh_fem object)
   'edges_color',[.6 .6 1]     : RGB values for the color of edges
-  'edges_width',1             : width of edges              
+  'edges_width',1             : width of edges
   'faces_color',[.75 .75 .75]): RGB values for the color of faces
   'quality',{ 'off' | 'on' }  : Display the quality of the mesh.
 
@@ -137,14 +137,14 @@ gf_plot_mesh
 
   This function is used to display a mesh.
 
-  Example :: 
+  Example ::
 
     % the mesh is in the tests directory of the distribution
     m=gf_mesh('import','gid','donut_with_quadratic_tetra_314_elements.msh');
-    gf_plot_mesh(m,'refine',15,'cvlst',gf_mesh_get(m,'outer 
faces'),'faces','on',\ldots, 'faces_color',[1. .9 
.2],'curved','on','edges_width',2); 
+    gf_plot_mesh(m,'refine',15,'cvlst',gf_mesh_get(m,'outer 
faces'),'faces','on',\ldots, 'faces_color',[1. .9 
.2],'curved','on','edges_width',2);
     camlight % turn on the light!
 
- 
+
 
 gf_plot_slice
 -------------------------------------------
@@ -172,7 +172,7 @@ gf_plot_slice
   pcolor, ['on']      : if the field is scalar, a color plot of its values is 
plotted
   quiver, ['on']      : if the field is vector, represent arrows
   quiver_density, 50  : density of arrows in quiver plot
-  quiver_scale, 1     : density of arrows in quiver plot 
+  quiver_scale, 1     : density of arrows in quiver plot
   tube, ['on']        : use tube plot for 'filar' (1D) parts of the slice
   tube_color, ['red'] : color of tubes (ignored if 'data' is not empty and 
'pcolor' is on)
   tube_radius, ['0.5%'] : tube radius; you can use a constant, or a percentage 
(of the mesh size) or a vector of nodal values
@@ -195,29 +195,29 @@ gf_plot_slice
     Usl=gf_compute(pde.mf_u,U,'interpolate on', sl);  % interpolate the 
solution on the slice
     % show the norm of the displacement on this slice
     
gf_plot_slice(sl,'mesh','on','data',sqrt(sum(Usl.^2,1)),'mesh_slice_edges','off');
-    
+
     % another slice: now we take the lower part of the mesh
     
sl=gf_slice(\{'boundary',\{'intersection',\{'planar',+1,[0;0;6],[0;0;-1]\},\ldots
                                             
\{'planar',+1,[0;0;0],[0;1;0]\}\}\},m,6);
     Usl=gf_compute(pde.mf_u,U,'interpolate on', sl);
     hold on;
     
gf_plot_slice(sl,'mesh','on','data',sqrt(sum(Usl.^2,1)),'mesh_slice_edges','off');
-    
+
     % this slice contains the transparent mesh faces displayed on the picture
     sl2=gf_slice(\{'boundary',\{'planar',+1,[0;0;0],[0;1;0]\}\},\ldots
                 m,6,setdiff(all_faces',TOPfaces','rows')');
-    gf_plot_slice(sl2,'mesh_faces','off','mesh','on','pcolor','off'); 
-    
+    gf_plot_slice(sl2,'mesh_faces','off','mesh','on','pcolor','off');
+
     % last step is to plot the streamlines
     hh=[1 5 9 12.5 16 19.5]; % vertical position of the different starting 
points of the streamlines
     H=[zeros(2,numel(hh));hh];
-    
+
     % compute the streamlines
     tsl=gf_slice('streamlines',pde.mf_u,U,H);
     Utsl=gf_compute(pde.mf_u,U,'interpolate on', tsl);
-    
+
     % render them with "tube plot"
-    
[a,h]=gf_plot_slice(tsl,'mesh','off','tube_radius',.2,'tube_color','white'); 
+    
[a,h]=gf_plot_slice(tsl,'mesh','off','tube_radius',.2,'tube_color','white');
     hold off;
     % use a nice colormap
     caxis([0 .7]);
diff --git a/doc/sphinx/source/scilab/scilabgf.rst 
b/doc/sphinx/source/scilab/scilabgf.rst
index 41e4d81..1e21d62 100644
--- a/doc/sphinx/source/scilab/scilabgf.rst
+++ b/doc/sphinx/source/scilab/scilabgf.rst
@@ -80,8 +80,8 @@ Functions
 Objects
 -------
 
-Various "objects" can be manipulated by the |gfm| toolbox, see fig. 
-:ref:`scilab-fig-hierarchy`. The MESH and MESHFEM objects are the two most 
+Various "objects" can be manipulated by the |gfm| toolbox, see fig.
+:ref:`scilab-fig-hierarchy`. The MESH and MESHFEM objects are the two most
 important objects.
 
 .. _scilab-fig-hierarchy:
diff --git a/doc/sphinx/source/screenshots/shots.rst 
b/doc/sphinx/source/screenshots/shots.rst
index 33d365c..76794fc 100644
--- a/doc/sphinx/source/screenshots/shots.rst
+++ b/doc/sphinx/source/screenshots/shots.rst
@@ -9,8 +9,8 @@ GetFEM++ in action ...
 Generic mesh handling
 ---------------------
 
-The first images illustrate the general mesh handling of getfem. The `mesh 
description`_ 
-is hand-made, and involves many different element types and convex types, as 
you can see 
+The first images illustrate the general mesh handling of getfem. The `mesh 
description`_
+is hand-made, and involves many different element types and convex types, as 
you can see
 (the mesh, and a random function interpolated on the mesh):
 
 .. _mesh description: ../_static/strange.mesh_fem
@@ -24,17 +24,17 @@ is hand-made, and involves many different element types and 
convex types, as you
 
 .. centered:: |im1|_ |im2|_
 
-The mesh is 3D. There is a quadrangle, a curved quadrangle/triangle, a kind of 
curved 
-prism and hexahedron, and a very curved (geometrical transformation of degree 
3) 
+The mesh is 3D. There is a quadrangle, a curved quadrangle/triangle, a kind of 
curved
+prism and hexahedron, and a very curved (geometrical transformation of degree 
3)
 quadrangle.
 
 Linear elasticity
 -----------------
 
-A tripod is fixed on the ground and loaded with a vertical force on its top. 
The mesh was 
-generated with `GiD`_, using quadratic (i.e. curved) tetrahedrons. The 
solution is 
-computed on a P2 FEM (i.e. P2 isoparametric FEM). Below is the Von Mises 
stress, 
-represented on the deformed tripod. The source code of this example uses the 
matlab 
+A tripod is fixed on the ground and loaded with a vertical force on its top. 
The mesh was
+generated with `GiD`_, using quadratic (i.e. curved) tetrahedrons. The 
solution is
+computed on a P2 FEM (i.e. P2 isoparametric FEM). Below is the Von Mises 
stress,
+represented on the deformed tripod. The source code of this example uses the 
matlab
 interface, and can be found here: :ref:`tripod-source`.
 
 .. |im-tri| image:: images/tripodvonmiseswithmesh_small.*
@@ -43,7 +43,7 @@ interface, and can be found here: :ref:`tripod-source`.
 .. centered:: |im-tri|_
 
 
-If you want to see what is inside the tripod, download the following animation 
(mpeg-4 
+If you want to see what is inside the tripod, download the following animation 
(mpeg-4
 movie, 6MB, 45secs) `tripod_slice.avi`_
 
 .. _tripod_slice.avi: http://getfem.org/misc/tripod_slice.avi
@@ -52,9 +52,9 @@ movie, 6MB, 45secs) `tripod_slice.avi`_
 Stokes equation
 ---------------
 
-An incompressible viscous fluid flows in a 2D tube. The mesh is made of curved 
triangles, 
-and the solution is computed on a mixed P2+/P1 FEM (P2 with a cubic bubble for 
the 
-velocity field, and discontinuous P1 for the pressure field). The source code 
is here: 
+An incompressible viscous fluid flows in a 2D tube. The mesh is made of curved 
triangles,
+and the solution is computed on a mixed P2+/P1 FEM (P2 with a cubic bubble for 
the
+velocity field, and discontinuous P1 for the pressure field). The source code 
is here:
 :ref:`stokes-source`.
 
 .. |im-sto| image:: images/tube_small.*
@@ -62,7 +62,7 @@ velocity field, and discontinuous P1 for the pressure field). 
The source code is
 
 .. centered:: |im-sto|_
 
-The next example is still the Stokes problem, inside a 3D cylindrical tank. 
The picture show the norm of the fluid velocity, with some streamlines. 
+The next example is still the Stokes problem, inside a 3D cylindrical tank. 
The picture show the norm of the fluid velocity, with some streamlines.
 3D tank
 
 .. |im-cuve| image:: images/cuve_3D_streamlines_small.*
@@ -73,10 +73,10 @@ The next example is still the Stokes problem, inside a 3D 
cylindrical tank. The
 Helmholtz equation
 ------------------
 
-This is a basic 2D scattering example. An incoming plane wave is scaterred by 
a perfectly 
-reflective circular obstacle. The mesh is made of only 25 quadrangles whose 
geometric 
-transformations are polynomials of degree 6. Computations are done with a P10 
FEM, hence 
-it is possible to have 2 wavelength per element ! (with a P1 fem, the rule is 
at least 6 
+This is a basic 2D scattering example. An incoming plane wave is scaterred by 
a perfectly
+reflective circular obstacle. The mesh is made of only 25 quadrangles whose 
geometric
+transformations are polynomials of degree 6. Computations are done with a P10 
FEM, hence
+it is possible to have 2 wavelength per element ! (with a P1 fem, the rule is 
at least 6
 elements per wavelength). The source is here: :ref:`helmholtz-source`.
 
 .. |im-helm1| image:: images/helm_mesh_k7_P10_gt6.*
@@ -95,40 +95,40 @@ Eigenmodes of a structure (thanks to Paolo Bertolo)
 
 .. centered:: |im-paolo|
 
-eigenmode of a vibrating structure You can look at a small movie showing the 
24 first 
-modes of the structure: (mpeg1, 4MB) `oggetto_modes.mpeg`_ or (mpeg4, 8MB) 
+eigenmode of a vibrating structure You can look at a small movie showing the 
24 first
+modes of the structure: (mpeg1, 4MB) `oggetto_modes.mpeg`_ or (mpeg4, 8MB)
 `oggetto_modes.avi`_.
 
 .. _oggetto_modes.mpeg: http://getfem.org/misc/oggetto_modes.mpeg
 
 .. _oggetto_modes.avi: http://getfem.org/misc/oggetto_modes.avi
-  
+
 Contact with friction problem (Houari Khenous)
 ----------------------------------------------
 
-This example shows the deformation of a tire under its own weight. The tire is 
meshed 
-with one layer of regular hexahedric cells (384 cells), whose geometric 
transformation is 
-of order 2, and a Q2 FEM. This picture shows the Von Mises criterion on the 
deformed 
+This example shows the deformation of a tire under its own weight. The tire is 
meshed
+with one layer of regular hexahedric cells (384 cells), whose geometric 
transformation is
+of order 2, and a Q2 FEM. This picture shows the Von Mises criterion on the 
deformed
 tire.
- 
+
 .. |im-houari| image:: images/pneu_Q2_vonmises_small.*
 
 .. centered:: |im-houari|
 
-An animation of a (soft) elastic disk is also available (mpeg-4 movie, 4MB, 
12secs) 
-`disk_in_contact.avi`_ (mpeg1, 1MB) (A newmark scheme adapted for the 
unilateral contact 
+An animation of a (soft) elastic disk is also available (mpeg-4 movie, 4MB, 
12secs)
+`disk_in_contact.avi`_ (mpeg1, 1MB) (A newmark scheme adapted for the 
unilateral contact
 condition).
 
 .. _disk_in_contact.avi: http://getfem.org/misc/disk_in_contact.avi
 
- 
+
 Xfem cracks in a beam
 ---------------------
 
-Here we used XFem to handle cracks in a beam. XFem is an enrichment of the 
classical 
-finite element space (a P2 FEM was used for this example) with a discontinuous 
function. 
-Thanks to this function, the crack path does not have to follow the original 
mesh. Note 
-how the crack cross elements on the mesh below. Four singular functions, which 
form a 
+Here we used XFem to handle cracks in a beam. XFem is an enrichment of the 
classical
+finite element space (a P2 FEM was used for this example) with a discontinuous 
function.
+Thanks to this function, the crack path does not have to follow the original 
mesh. Note
+how the crack cross elements on the mesh below. Four singular functions, which 
form a
 basis for asymptotical solution to the linear elasticity problem near the 
crack tips.
 
 .. |im-crack| image:: images/xfembeammesh.*
@@ -144,7 +144,7 @@ basis for asymptotical solution to the linear elasticity 
problem near the crack
 A 3D crack, made via level-set
 ------------------------------
 
-In this example, the mesh was a simple cartesian mesh 20x20x1, and the crack 
geometry was 
+In this example, the mesh was a simple cartesian mesh 20x20x1, and the crack 
geometry was
 defined implicitely via a levelset.
 
 .. |im-crack3d| image:: images/fissure_3d_de_traviole.*
@@ -154,26 +154,26 @@ defined implicitely via a levelset.
 Large strain
 ------------
 
-In this example, a bar is twisted. Each step is solved with a Newton method. 
The material 
-law is a "Ciarlet Geymonat" one. A P2 FEM is used. The source code for this 
example can 
-be found in the `tests/nonlinear_elastostatic.cc` file of |gf| package. This 
picture was 
+In this example, a bar is twisted. Each step is solved with a Newton method. 
The material
+law is a "Ciarlet Geymonat" one. A P2 FEM is used. The source code for this 
example can
+be found in the `tests/nonlinear_elastostatic.cc` file of |gf| package. This 
picture was
 made with OpenDX.
 
 .. |im-largestrain| image:: images/torsion034.*
 
 .. centered:: |im-largestrain|
- 
-A short animation is also available: (mpeg-4 movie, 3MB) `torsion.avi`_. 
+
+A short animation is also available: (mpeg-4 movie, 3MB) `torsion.avi`_.
 
 .. _torsion.avi: http://getfem.org/misc/torsion.avi
 
 Shape and topological optimization
 ----------------------------------
 
-This images were obtained with the script 
-`interface/tests/matlab/demo_structural_optimization.m` (Alassane SY and Yves 
Renard). It 
-represents a shape optimization of a structure submitted to a vertical load at 
the right 
-and clambed at the left. A (Xfem like) fictitious domain approach is used 
together with 
+This images were obtained with the script
+`interface/tests/matlab/demo_structural_optimization.m` (Alassane SY and Yves 
Renard). It
+represents a shape optimization of a structure submitted to a vertical load at 
the right
+and clambed at the left. A (Xfem like) fictitious domain approach is used 
together with
 both a shape gradient and a topological gradient.
 
 .. |im-shape1| image:: images/shape1.*
@@ -182,9 +182,9 @@ both a shape gradient and a topological gradient.
 
 .. centered:: |im-shape1| |im-shape2|
 
-  
-The first image corresponds to an initial structure with pre-existing holes. 
For the 
-second one the holes are initiated by the topological optimization. The two 
following 
+
+The first image corresponds to an initial structure with pre-existing holes. 
For the
+second one the holes are initiated by the topological optimization. The two 
following
 images correspond to a 3D case.
 
 .. |im-shape3| image:: images/shape3.*
@@ -196,7 +196,7 @@ images correspond to a 3D case.
 3D planetary gears
 ------------------
 
-This image comes from the application developed by Konstantinos Poulios 
+This image comes from the application developed by Konstantinos Poulios
 which is freely available at |link23|_. It is based on |gf| and is intended to 
be a tool for easy, almost automatic, creation and calculation of gear 
transmissions.
 
 
diff --git a/doc/sphinx/source/tutorial/basic_usage.rst 
b/doc/sphinx/source/tutorial/basic_usage.rst
index 5a9111f..72f39cb 100644
--- a/doc/sphinx/source/tutorial/basic_usage.rst
+++ b/doc/sphinx/source/tutorial/basic_usage.rst
@@ -20,7 +20,7 @@ A program using getfem will often have the following 
structure ::
 
 
   ... define one or more Mesh
-  
+
   ... define one or more MeshFem
 
   ... define one or more MeshIm
diff --git a/doc/sphinx/source/tutorial/install.rst 
b/doc/sphinx/source/tutorial/install.rst
index 1b06631..1d98fb3 100644
--- a/doc/sphinx/source/tutorial/install.rst
+++ b/doc/sphinx/source/tutorial/install.rst
@@ -17,7 +17,7 @@ The main dependences of |Gf| on other libraries are
   from git version to get the latest changes.
 
 * Python development files (Python.h etc.) and also the |np| and |sp| packages 
if
-  you want to build the python interface. 
+  you want to build the python interface.
 
 * sequential MUMPS package (direct solver for sparse matrices) if you want to 
use it instead of the SuperLU version distributed along with |gf|.
 
@@ -29,7 +29,7 @@ The main dependences of |Gf| on other libraries are
 
 |gf| C++ library can be build on its own or together with the Python, Scilab 
and/or Matlab interface.
 
-You can also install the stable release of Getfem on linux distributions using 
the corresponding package management system. 
+You can also install the stable release of Getfem on linux distributions using 
the corresponding package management system.
 
 
 More specific information on how to build Getfem C++ library can be found on 
the `download and install <../download.html>`_ page.
diff --git a/doc/sphinx/source/tutorial/thermo_coupling.rst 
b/doc/sphinx/source/tutorial/thermo_coupling.rst
index 7659e22..d57573c 100644
--- a/doc/sphinx/source/tutorial/thermo_coupling.rst
+++ b/doc/sphinx/source/tutorial/thermo_coupling.rst
@@ -17,7 +17,7 @@ Let :math:`\Omega \subset \R^2` be the reference 
configuration of a 2D plate (se
 Thermal problem
 ***************
 
-The lateral faces of the plates are supposed to be in thermal insulation since 
the front and back faces of the plate are supposed to be in thermal exchange 
with the air (supposed at 20 |degreC|) with a heat transfer coefficient 
:math:`D`. 
+The lateral faces of the plates are supposed to be in thermal insulation since 
the front and back faces of the plate are supposed to be in thermal exchange 
with the air (supposed at 20 |degreC|) with a heat transfer coefficient 
:math:`D`.
 
 The equation on the temperature :math:`\theta` and boundary condition can be 
written as follows:
 
@@ -26,7 +26,7 @@ The equation on the temperature :math:`\theta` and boundary 
condition can be wri
   \left\{\begin{array}{l}
   -\mbox{div}(\varepsilon\kappa(\nabla \theta)) + 2D(\theta - T_0) - 
\varepsilon\sigma|\nabla V|^2 = 0 ~~ \mbox{ in } \Omega, \\
   \kappa\nabla \theta \cdot n = 0 ~~ \mbox{ on } \partial \Omega,
-  \end{array} \right. 
+  \end{array} \right.
 
 
 where the thermal conductivity is designed by :math:`\kappa`, :math:`T_0` is 
the temperature of the air, :math:`\partial \Omega` the boundary of the domain 
:math:`\Omega` and :math:`n` the outward unit normal vector to :math:`\Omega` 
on :math:`\partial \Omega`.
@@ -44,7 +44,7 @@ We consider a potential difference of :math:`0.1V` between 
the right and left la
   -\mbox{div}(\varepsilon\sigma(\nabla V)) = 0 ~~ \mbox{ in } \Omega, \\
   \sigma\nabla V \cdot n = 0 ~~ \mbox{ on the top an bottom lateral faces}, \\
   V = 0 ~~ \mbox{ on the right lateral face}, \\
-  V = 0.1 ~~ \mbox{ on the left lateral face}, \\  
+  V = 0.1 ~~ \mbox{ on the left lateral face}, \\
   \end{array} \right.
 
 where :math:`\sigma` is still the electrical conductivity. Moreover, we 
consider that :math:`\sigma` depends on the temperature as follows:
@@ -68,13 +68,13 @@ We consider the membrane small deformation of the plate 
under a force applied on
   -\mbox{div}(\bar{\sigma}(u)) = 0 ~~ \mbox{ in } \Omega, \\
   \bar{\sigma}\ n = 0 ~~ \mbox{ on the top an bottom lateral faces}, \\
   \bar{\sigma}\ n = F ~~ \mbox{ on the right lateral face}, \\
-  u = 0 ~~ \mbox{ on the left lateral face}, 
+  u = 0 ~~ \mbox{ on the left lateral face},
   \end{array} \right.
 
 where :math:`F` is the force density applied on the right lateral boundary and 
:math:`\bar{\sigma}(u)` is the Cauchy stress tensor defined by
 
 .. math::
- 
+
   \bar{\sigma}(u) = \lambda^* \mbox{div}(u) I + 2\mu \bar{\varepsilon}(u) + 
\beta(T_0-\theta) I,
 
 :math:`\bar{\varepsilon}(u) = (\nabla u + (\nabla u)^T)/2` being the 
linearized strain tensor, :math:`I` the identity second order tensor and 
:math:`\lambda^*, \mu` being the |Lame| coefficients defined  by
@@ -115,18 +115,18 @@ Let us now make a detailed presentation of the use of 
|gf| to approximate the pr
 Initialization
 **************
 
-First, in C++, ones has to include a certain number of headers for the model 
object, the generic assembly, the linear algebra interface (Gmm++), the 
experimental mesher and the export facilities. For Python, this is simpler, 
|gf| can be imported globally (numpy has also to be imported). For Scilab, the 
library has first to be loaded in the Scilab console (this is not described 
here) and for Matlab, nothing is necessary, except a `gf_workspace('clear 
all')` which allows to clear all |gf|  [...]
+First, in C++, ones has to include a certain number of headers for the model 
object, the generic assembly, the linear algebra interface (Gmm++), the 
experimental mesher and the export facilities. For Python, this is simpler, 
|gf| can be imported globally (numpy has also to be imported). For Scilab, the 
library has first to be loaded in the Scilab console (this is not described 
here) and for Matlab, nothing is necessary, except a `gf_workspace('clear 
all')` which allows to clear all |gf|  [...]
 
 
 ========== ================================================
-**C++**    .. code-block:: c++                             
+**C++**    .. code-block:: c++
 
-             #include "getfem/getfem_model_solvers.h" 
+             #include "getfem/getfem_model_solvers.h"
              #include "getfem/getfem_export.h"
              #include "gmm/gmm.h"
              #include "getfem/getfem_mesher.h"
              #include "getfem/getfem_generic_assembly.h"
-           
+
              using bgeot::size_type;
              using bgeot::base_node;
              using bgeot::base_small_vector;
@@ -134,16 +134,16 @@ First, in C++, ones has to include a certain number of 
headers for the model obj
 
              int main(void) {
 ---------- ------------------------------------------------
-**Python** .. code-block:: python                                     
+**Python** .. code-block:: python
 
              import getfem as gf
              import numpy as np
 ---------- ------------------------------------------------
-**Scilab** .. code-block:: matlab                             
+**Scilab** .. code-block:: matlab
 
              gf_workspace('clear all'); // In case of multiple runs
 ---------- ------------------------------------------------
-**Matlab** .. code-block:: matlab                                    
+**Matlab** .. code-block:: matlab
 
              gf_workspace('clear all'); % In case of multiple runs
 ========== ================================================
@@ -156,7 +156,7 @@ Let us now define the different physical and numerical 
parameters of the problem
 
 
 =========== ================================================
-**C++**     .. code-block:: c++                             
+**C++**     .. code-block:: c++
 
                 double epsilon = 1.; // Thickness of the plate (cm)
                 double E = 21E6;     // Young Modulus (N/cm^2)
@@ -176,7 +176,7 @@ Let us now define the different physical and numerical 
parameters of the problem
                 double h = 2.        // Approximate mesh size
                 bgeot::dim_type elements_degree = 2; // Degree of the finite 
element methods
 ----------- ------------------------------------------------
-**Scripts** .. code-block:: python                                     
+**Scripts** .. code-block:: python
 
                 epsilon = 1.; E = 21E6; nu = 0.3;
                 clambda = E*nu/((1+nu)*(1-2*nu));
@@ -187,7 +187,7 @@ Let us now define the different physical and numerical 
parameters of the problem
                 T0 = 20; rho_0 = 1.754E-8;
                 alpha = 0.0039;
 
-                h = 2; elements_degree = 2;     
+                h = 2; elements_degree = 2;
 
 =========== ================================================
 
@@ -203,7 +203,7 @@ The geometry of the domain is supposed to be a rectangle 
with three circular hol
 
 
 ========== 
===========================================================================
-**C++**    .. code-block:: c++                             
+**C++**    .. code-block:: c++
 
                 getfem::mesh mesh;
                getfem::pmesher_signed_distance
@@ -217,7 +217,7 @@ The geometry of the domain is supposed to be a rectangle 
with three circular hol
                 std::vector<getfem::base_node> fixed;
                 getfem::build_mesh(mesh, mo, h, fixed, 2, -2);
 ---------- 
---------------------------------------------------------------------------
-**Python** .. code-block:: python                                     
+**Python** .. code-block:: python
 
                 mo1 = gf.MesherObject('rectangle', [0., 0.], [100., 25.])
                 mo2 = gf.MesherObject('ball', [25., 12.5], 8.)
@@ -228,7 +228,7 @@ The geometry of the domain is supposed to be a rectangle 
with three circular hol
 
                 mesh = gf.Mesh('generate', mo, h, 2)
 ---------- 
---------------------------------------------------------------------------
-**Scilab** .. code-block:: matlab                             
+**Scilab** .. code-block:: matlab
 
                 mo1 = gf_mesher_object('rectangle', [0 0], [100 25]);
                 mo2 = gf_mesher_object('ball', [25 12.5], 8);
@@ -265,7 +265,7 @@ Since we have different boundary conditions on the 
different parts of the bounda
 
 
 ========== 
===========================================================================
-**C++**    .. code-block:: c++                             
+**C++**    .. code-block:: c++
 
                 getfem::mesh_region border_faces;
                 getfem::outer_faces_of_mesh(mesh, border_faces);
@@ -291,7 +291,7 @@ Since we have different boundary conditions on the 
different parts of the bounda
                 mesh.region(   TOP_BOUND) = getfem::mesh_region::subtract(fb4, 
fb1);
                 mesh.region(BOTTOM_BOUND) = getfem::mesh_region::subtract(fb5, 
fb1);
 ---------- 
---------------------------------------------------------------------------
-**Python** .. code-block:: python                                     
+**Python** .. code-block:: python
 
                 fb1 = mesh.outer_faces_in_box([1., 1.], [99., 24.])
                 fb2 = mesh.outer_faces_with_direction([ 1., 0.], 0.01)
@@ -311,7 +311,7 @@ Since we have different boundary conditions on the 
different parts of the bounda
                 mesh.region_subtract(   TOP_BOUND, HOLE_BOUND)
                 mesh.region_subtract(BOTTOM_BOUND, HOLE_BOUND)
 ---------- 
---------------------------------------------------------------------------
-**Scilab** .. code-block:: matlab                             
+**Scilab** .. code-block:: matlab
 
                 fb1 = gf_mesh_get(mesh, 'outer faces in box', [1 1], [99 24]);
                 fb2 = gf_mesh_get(mesh, 'outer faces with direction', [ 1 0], 
0.01);
@@ -357,7 +357,7 @@ Mesh draw
 In order to preview the mesh and to control its validity, the following 
instructions can be used:
 
 ========== 
===========================================================================
-**C++**    .. code-block:: c++                             
+**C++**    .. code-block:: c++
 
                 getfem::vtk_export exp("mesh.vtk", false);
                 exp.exporting(mesh);
@@ -365,7 +365,7 @@ In order to preview the mesh and to control its validity, 
the following instruct
                 // You can view the mesh for instance with
                 // mayavi2 -d mesh.vtk -f ExtractEdges -m Surface
 ---------- 
---------------------------------------------------------------------------
-**Python** .. code-block:: python                                     
+**Python** .. code-block:: python
 
                 mesh.export_to_vtk('mesh.vtk');
                 # You can view the mesh for instance with
@@ -396,7 +396,7 @@ Definition of finite element methods and integration method
 
 We will define three finite element methods. The first one, `mfu` is to 
approximate the displacement field. This is a vector field. This is defined in 
C++ by
 
-.. code-block:: c++   
+.. code-block:: c++
 
         getfem::mesh_fem mfu(mesh, 2);
         mfu.set_classical_finite_element(elements_degree);
@@ -413,7 +413,7 @@ The last thing to define is an integration method `mim`. 
There is no default int
 
 
 ========== 
===========================================================================
-**C++**    .. code-block:: c++                             
+**C++**    .. code-block:: c++
 
                 getfem::mesh_fem mfu(mesh, 2);
                 mfu.set_classical_finite_element(elements_degree);
@@ -425,7 +425,7 @@ The last thing to define is an integration method `mim`. 
There is no default int
                 getfem::mesh_im  mim(mesh);
                 
mim.set_integration_method(bgeot::dim_type(gmm::sqr(elements_degree)));
 ---------- 
---------------------------------------------------------------------------
-**Python** .. code-block:: python                                     
+**Python** .. code-block:: python
 
                 mfu = gf.MeshFem(mesh, 2)
                 mfu.set_classical_fem(elements_degree)
@@ -466,20 +466,20 @@ The model object in |gf| gather the variables of the 
models (the unknowns), the
 
 This is not strictly mandatory to use the model object since one may use 
directly the assembly procedures and build by it own the (tangent) linear 
system. The model object allows a rapid build of the model since most classical 
parts of a model are pre-programmed: standard boundary conditions, standard 
partial differential equations, use of multipliers to prescribe a constraint 
... Moreover, some bricks are designed to extend the possibilities of the 
standard bricks (generic assembly bric [...]
 
-There are two versions of the model: the real one and the complex one. Complex 
models have to be reserved for special applications (some electromagnetism 
problems for instance) where it is advantageous to solve a complex linear 
system. 
+There are two versions of the model: the real one and the complex one. Complex 
models have to be reserved for special applications (some electromagnetism 
problems for instance) where it is advantageous to solve a complex linear 
system.
 
 Let us declare a real model with the three variables corresponding to the 
three fields to be computed:
 
 
 ========== 
===========================================================================
-**C++**    .. code-block:: c++                             
+**C++**    .. code-block:: c++
 
                 getfem::model md;
                 md.add_fem_variable("u", mfu);
                 md.add_fem_variable("theta", mft);
                 md.add_fem_variable("V", mft);
 ---------- 
---------------------------------------------------------------------------
-**Python** .. code-block:: python                                     
+**Python** .. code-block:: python
 
                 md=gf.Model('real');
                 md.add_fem_variable('u', mfu)
@@ -491,7 +491,7 @@ Let us declare a real model with the three variables 
corresponding to the three
                 md=gf_model('real');
                 gf_model_set(md, 'add fem variable', 'u', mfu);
                 gf_model_set(md, 'add fem variable', 'theta', mft);
-                gf_model_set(md, 'add fem variable', 'V', mft);               
+                gf_model_set(md, 'add fem variable', 'V', mft);
 ---------- 
---------------------------------------------------------------------------
 **Matlab** .. code-block:: matlab
 
@@ -525,7 +525,7 @@ The following program allows to take into account the whole 
elastic deformation
 
 
 ========== 
================================================================================================================
-**C++**    .. code-block:: c++                             
+**C++**    .. code-block:: c++
 
                 md.add_initialized_scalar_data("cmu", cmu);
                 md.add_initialized_scalar_data("clambdastar", clambdastar);
@@ -540,7 +540,7 @@ The following program allows to take into account the whole 
elastic deformation
                 md.add_initialized_scalar_data("beta", alpha_th*E/(1-2*nu));
                 getfem::add_linear_term(md, mim, "beta*(T0-theta)*Div_Test_u");
 ---------- 
----------------------------------------------------------------------------------------------------------------
-**Python** .. code-block:: python                                     
+**Python** .. code-block:: python
 
                 md.add_initialized_data('cmu', [cmu])
                 md.add_initialized_data('clambdastar', [clambdastar])
@@ -589,7 +589,7 @@ Electric potential problem
 Similarly, the following program take into account the electric potential 
equation. Note the definition of the  electrical conductivity :math:`\sigma` 
and again the use of weak form language terms.
 
 ========== 
===========================================================================
-**C++**    .. code-block:: c++                             
+**C++**    .. code-block:: c++
 
                 std::string sigmaeps = "(eps/(rho_0*(1+alpha*(theta-T0))))";
                 md.add_initialized_scalar_data("eps", epsilon);
@@ -603,7 +603,7 @@ Similarly, the following program take into account the 
electric potential equati
                 getfem::add_Dirichlet_condition_with_multipliers
                   (md, mim, "V", bgeot::dim_type(elements_degree-1), 
LEFT_BOUND, "DdataV");
 ---------- 
---------------------------------------------------------------------------
-**Python** .. code-block:: python                                     
+**Python** .. code-block:: python
 
                 sigmaeps = '(eps/(rho_0*(1+alpha*(theta-T0))))'
                 md.add_initialized_data('eps', [epsilon])
@@ -643,7 +643,7 @@ Thermal problem
 Now, the program to take into account the thermal problem:
 
 ========== 
===========================================================================
-**C++**    .. code-block:: c++                             
+**C++**    .. code-block:: c++
 
                 md.add_initialized_scalar_data("kaeps", kappa*epsilon);
                 getfem::add_generic_elliptic_brick(md, mim, "theta", "kaeps");
@@ -655,7 +655,7 @@ Now, the program to take into account the thermal problem:
                 getfem::add_nonlinear_term
                   (md, mim, "-"+sigmaeps+"*Norm_sqr(Grad_V)*Test_theta");
 ---------- 
---------------------------------------------------------------------------
-**Python** .. code-block:: python                                     
+**Python** .. code-block:: python
 
                 md.add_initialized_data('kaeps', [kappa*epsilon])
                 md.add_generic_elliptic_brick(mim, 'theta', 'kaeps')
@@ -686,7 +686,7 @@ Now, the program to take into account the thermal problem:
                 gf_model_set(md, 'add mass brick', mim, 'theta', 'D2');
                 gf_model_set(md, 'add source term brick', mim, 'theta', 
'D2airt');
 
-                gf_model_set(md, 'add nonlinear term', mim, ['-' sigmaeps 
'*Norm_sqr(Grad_V)*Test_theta']);          
+                gf_model_set(md, 'add nonlinear term', mim, ['-' sigmaeps 
'*Norm_sqr(Grad_V)*Test_theta']);
 ========== 
===========================================================================
 
 
@@ -696,22 +696,22 @@ Model solve
 Once the model is correctly defined, we can simply solve it by:
 
 ========== 
===========================================================================
-**C++**    .. code-block:: c++                             
+**C++**    .. code-block:: c++
 
                 gmm::iteration iter(1E-9, 1, 100);
                 getfem::standard_solve(md, iter);
 ---------- 
---------------------------------------------------------------------------
-**Python** .. code-block:: python                                     
+**Python** .. code-block:: python
 
                 md.solve('max_res', 1E-9, 'max_iter', 100, 'noisy')
 ---------- 
---------------------------------------------------------------------------
 **Scilab** .. code-block:: matlab
 
-                gf_model_get(md, 'solve', 'max_res', 1E-9, 'max_iter', 100, 
'noisy'); 
+                gf_model_get(md, 'solve', 'max_res', 1E-9, 'max_iter', 100, 
'noisy');
 ---------- 
---------------------------------------------------------------------------
 **Matlab** .. code-block:: matlab
 
-                gf_model_get(md, 'solve', 'max_res', 1E-9, 'max_iter', 100, 
'noisy');  
+                gf_model_get(md, 'solve', 'max_res', 1E-9, 'max_iter', 100, 
'noisy');
 ========== 
===========================================================================
 
 Since the problem is globally nonlinear, a Newton method is used to 
iteratively solve the problem. It needs a few iterations (about 4 in that case).
@@ -723,7 +723,7 @@ Model solve with two steps
 Another option to solve the problem is to solve first the thermal and electric 
potential problems. Indeed, in our model, the thermal and  electric potential 
do not depend on the deformation. Once the  thermal and electric potential 
problem, we then solve the deformation problem. This can be done as follows:
 
 ========== 
===========================================================================
-**C++**    .. code-block:: c++                             
+**C++**    .. code-block:: c++
 
                 gmm::iteration iter(1E-9, 1, 100);
                 md.disable_variable("u");
@@ -732,9 +732,9 @@ Another option to solve the problem is to solve first the 
thermal and electric p
                 md.disable_variable("theta");
                 md.disable_variable("V");
                 iter.init();
-                getfem::standard_solve(md, iter);      
+                getfem::standard_solve(md, iter);
 ---------- 
---------------------------------------------------------------------------
-**Python** .. code-block:: python                                     
+**Python** .. code-block:: python
 
                 md.disable_variable('u')
                 md.solve('max_res', 1E-9, 'max_iter', 100, 'noisy')
@@ -770,7 +770,7 @@ Export/visualization of the solution
 The finite element problem is now solved. We can plot the solution as follows. 
Note that for the C++ and Python programs, it is necessary to use an external 
graphical post-processor. Note also that arbitrary quantities can be 
post-processed using the generic interpolation (see 
`ga_interpolation_Lagrange_fem` below). It is also possible to make complex 
exports and slices (see :ref:`ud-export`).
 
 ========== 
=====================================================================================================================================================
-**C++**    .. code-block:: c++                             
+**C++**    .. code-block:: c++
 
                 plain_vector U(mfu.nb_dof()); gmm::copy(md.real_variable("u"), 
U);
                 plain_vector V(mft.nb_dof()); gmm::copy(md.real_variable("V"), 
V);
@@ -780,14 +780,14 @@ The finite element problem is now solved. We can plot the 
solution as follows. N
                   (md, "u", "clambdastar", "cmu", mfvm, VM, false);
                 plain_vector CO(mfvm.nb_dof() * 2);
                 getfem::ga_interpolation_Lagrange_fem(md, 
"-"+sigmaeps+"*Grad_V",  mfvm, CO);
-  
+
                 getfem::vtk_export exp("displacement_with_von_mises.vtk", 
false);
                 exp.exporting(mfu);
                 exp.write_point_data(mfu, U, "elastostatic displacement");
                 exp.write_point_data(mfvm, VM, "Von Mises stress");
                 cout << "\nYou can view solutions with for 
instance:\n\nmayavi2 "
                   "-d displacement_with_von_mises.vtk -f WarpVector -m 
Surface\n" << endl;
-  
+
                 getfem::vtk_export exp2("temperature.vtk", false);
                 exp2.exporting(mft);
                 exp2.write_point_data(mft, THETA, "Temperature");
@@ -798,9 +798,9 @@ The finite element problem is now solved. We can plot the 
solution as follows. N
                 exp3.write_point_data(mft, V, "Electric potential");
                 cout << "mayavi2 -d electric_potential.vtk -f WarpScalar -m 
Surface\n"
                      << endl;
-                }      
+                }
 ---------- 
-----------------------------------------------------------------------------------------------------------------------------------------------------
-**Python** .. code-block:: python                                     
+**Python** .. code-block:: python
 
                 U = md.variable('u')
                 V = md.variable('V')
@@ -815,7 +815,7 @@ The finite element problem is now solved. We can plot the 
solution as follows. N
                 mft.export_to_vtk('temperature.vtk', mft, THETA, 'Temperature')
                 print ('mayavi2 -d temperature.vtk -f WarpScalar -m Surface')
                 mft.export_to_vtk('electric_potential.vtk', mft, V, 'Electric 
potential')
-                print ('mayavi2 -d electric_potential.vtk -f WarpScalar -m 
Surface')    
+                print ('mayavi2 -d electric_potential.vtk -f WarpScalar -m 
Surface')
 ---------- 
-----------------------------------------------------------------------------------------------------------------------------------------------------
 **Scilab** .. code-block:: matlab
 
@@ -824,7 +824,7 @@ The finite element problem is now solved. We can plot the 
solution as follows. N
                 THETA = gf_model_get(md, 'variable', 'theta');
                 VM = gf_model_get(md, 
'compute_isotropic_linearized_Von_Mises_or_Tresca', 'u', 'clambdastar', 'cmu', 
mfvm);
                 CO = matrix(gf_model_get(md, 'interpolation', 
'-'+sigmaeps+'*Grad_V', mfvm), [2 gf_mesh_fem_get(mfvm, 'nbdof')]);
-    
+
                 hh = scf(2);
                 hh.color_map = jetcolormap(255);
                 subplot(3,1,1);
@@ -852,7 +852,7 @@ The finite element problem is now solved. We can plot the 
solution as follows. N
                 THETA = gf_model_get(md, 'variable', 'theta');
                 VM = gf_model_get(md, 
'compute_isotropic_linearized_Von_Mises_or_Tresca', 'u', 'clambdastar', 'cmu', 
mfvm);
                 CO = reshape(gf_model_get(md, 'interpolation', ['-' sigmaeps 
'*Grad_V'], mfvm), [2 gf_mesh_fem_get(mfvm, 'nbdof')]);
-    
+
                 figure(2);
                 subplot(3,1,1);
                 gf_plot(mfvm, VM, 'mesh', 'off', 'deformed_mesh','off', 
'deformation', U, 'deformation_mf', mfu, 'deformation_scale', 100, 'refine', 8);
@@ -871,7 +871,7 @@ The finite element problem is now solved. We can plot the 
solution as follows. N
                 gf_plot(mft, THETA, 'mesh', 'off', 'deformed_mesh','off', 
'deformation', U, 'deformation_mf', mfu, 'deformation_scale', 100, 'refine', 8);
                 colorbar;
                 title('Temperature in ?C (on the deformed configuration, scale 
factor x100)');
-     
+
 ========== 
=====================================================================================================================================================
 
 
diff --git a/doc/sphinx/source/tutorial/wheel.rst 
b/doc/sphinx/source/tutorial/wheel.rst
index d1121db..d46b9fc 100644
--- a/doc/sphinx/source/tutorial/wheel.rst
+++ b/doc/sphinx/source/tutorial/wheel.rst
@@ -35,7 +35,7 @@ Let us begin by loading Getfem and fixing the parameters of 
the problem
 
   h = 1                      # Approximate mesh size
   elements_degree = 2        # Degree of the finite element methods
-  gamma0 = 1./E;             # Augmentation parameter for the augmented 
Lagrangian 
+  gamma0 = 1./E;             # Augmentation parameter for the augmented 
Lagrangian
 
 
 
@@ -51,7 +51,7 @@ We consider that the radius of the wheel is 15cm and the one 
of the rim 8cm and
   mo3 = gf.MesherObject('set minus', mo1, mo2)
   gf.util('trace level', 2)   # No trace for mesh generation
   mesh1 = gf.Mesh('generate', mo3, h, 2)
-  
+
   mesh2 = 
gf.Mesh('import','structured','GT="GT_PK(2,1)";SIZES=[30,10];NOISED=0;NSUBDIV=[%d,%d];'
 % (int(30/h)+1, int(10/h)+1));
   mesh2.translate([-15.,-10.])
 
@@ -143,7 +143,7 @@ Contact condition (use of interpolate transformations)
 
 Now, let us see how to prescribed the contact condition between the two 
structures. It is possible to use predefined bricks (see  
:ref:`ud-model-contact-friction` for small deformation/small sliding contact 
and :ref:`ud-model-contact-friction-large` for large deformation/large sliding 
contact). However, we will see here how to directly prescribe a contact 
condition using an augmented Lagrangian formulation and the interpolate 
transformations.
 
-For small deformation contact, the correspondence between points of one 
contact surface to the other have to be described on the reference 
configuration and is not evolving, which is of course simpler but is an 
approximation. 
+For small deformation contact, the correspondence between points of one 
contact surface to the other have to be described on the reference 
configuration and is not evolving, which is of course simpler but is an 
approximation.
 
 We consider that the contact boundary of the wheel is the slave one and we 
have to describe the transformation from the contact boundary of the wheel to 
the contact boundary of the foundation. This is quite simple here, since the 
contact boundary of the foundation corresponds to a vanishing vertical 
coordinate. So we define the transformation
 
@@ -158,7 +158,7 @@ where :math:`X` is the vector of coordinates of the point. 
We add this transform
   md.add_interpolate_transformation_from_expression('Proj1', mesh1, mesh2, 
'[X(1);0]')
 
 As a consequence, it will be possible to use this transformation, from the 
mesh of the wheel to the mesh of the foundation, into weak form language 
expressions. Notes that this is here a very simple constant expression. More 
complex expressions depending on the data or even the variables of the model 
can be used. If the expression of a transformation depends on the variable of 
the model, the tangent linear system will automatically takes into account this 
dependence (see :ref:`ud-gasm-hi [...]
- 
+
 Using the defined transformation, we can write an integral contact condition 
using an augmented Lagrangian formulation (see :ref:`ud-model-contact-friction` 
for more details). The corresponding term (to be added to the rest of the weak 
formulation) reads:
 
 .. math::
diff --git a/doc/sphinx/source/userdoc/appendixA.rst 
b/doc/sphinx/source/userdoc/appendixA.rst
index 2a4d1d4..b14ad67 100644
--- a/doc/sphinx/source/userdoc/appendixA.rst
+++ b/doc/sphinx/source/userdoc/appendixA.rst
@@ -1391,7 +1391,7 @@ the shape functions read:
      - No :math:`(Q = 1)`
      - Yes
      - No
-       
+
    * - :math:`1`
      - :math:`3`
      - :math:`5`
@@ -1409,7 +1409,7 @@ the shape functions read:
      - No
 
 
-   
+
 .. list-table:: Discontinuous Lagrange element of order 0, 1 or 2 
``"FEM_PYRAMID_DISCONTINUOUS_LAGRANGE(K)"``
    :widths: 10 10 10 10 10 10 10
    :header-rows: 1
@@ -1429,7 +1429,7 @@ the shape functions read:
      - No :math:`(Q = 1)`
      - Yes
      - No
-       
+
    * - :math:`1`
      - :math:`3`
      - :math:`5`
@@ -1446,8 +1446,8 @@ the shape functions read:
      - Yes
      - No
 
-       
-        
+
+       
 
 Elements with additional bubble functions
 +++++++++++++++++++++++++++++++++++++++++
diff --git a/doc/sphinx/source/userdoc/appendixB.rst 
b/doc/sphinx/source/userdoc/appendixB.rst
index 679841c..5443344 100644
--- a/doc/sphinx/source/userdoc/appendixB.rst
+++ b/doc/sphinx/source/userdoc/appendixB.rst
@@ -27,7 +27,7 @@ the degree of each integration method.
 Exact Integration methods
 -------------------------
 
-|gf| furnishes a set of exact integration methods. This means that polynomials 
are integrated exactly. However, their use is (very) limited and not 
recommended. The use of exact integration methods is limited to the low-level 
generic assembly for polynomial :math:`\tau`-equivalent elements with linear 
transformations and for linear terms. It is not possible to use them in the 
high-level generic assembly. 
+|gf| furnishes a set of exact integration methods. This means that polynomials 
are integrated exactly. However, their use is (very) limited and not 
recommended. The use of exact integration methods is limited to the low-level 
generic assembly for polynomial :math:`\tau`-equivalent elements with linear 
transformations and for linear terms. It is not possible to use them in the 
high-level generic assembly.
 
 The list of available exact integration methods is the following
 
diff --git a/doc/sphinx/source/userdoc/bfem.rst 
b/doc/sphinx/source/userdoc/bfem.rst
index 0f70539..f9329de 100644
--- a/doc/sphinx/source/userdoc/bfem.rst
+++ b/doc/sphinx/source/userdoc/bfem.rst
@@ -199,7 +199,7 @@ the one of the |mf| object have to match. To sum it up,
 Additionally, if the field to be represented is a tensor field instead of a 
vector field (for instance the stress or strain tensor field in elasticity), it 
is possible to specify the tensor dimensions with the methods::
 
   mf.set_qdim(dim_type M, dim_type N)
-  mf.set_qdim(dim_type M, dim_type N, dim_type O, dim_type P)  
+  mf.set_qdim(dim_type M, dim_type N, dim_type O, dim_type P)
   mf.set_qdim(const bgeot::multi_index &mii)
 
 respectively for a tensor field of order two, four and arbitrary (but limited 
to 6). For most of the operations, this is equivalent to declare a vector field 
of the size the product of the dimensions. However, the declared tensor 
dimensions are taken into account into the high level generic assembly. 
Remember that the components inside a tensor are stored in Fortran order.
diff --git a/doc/sphinx/source/userdoc/bmesh.rst 
b/doc/sphinx/source/userdoc/bmesh.rst
index 576c3a9..0f026de 100644
--- a/doc/sphinx/source/userdoc/bmesh.rst
+++ b/doc/sphinx/source/userdoc/bmesh.rst
@@ -52,9 +52,9 @@ The most basic function to add a new element to a mesh is::
 
   j = mymesh.add_convex(pgt, it);
 
-This is a template function, with ``pgt`` of type |bg_pgt| (basically a 
pointer 
-to an instance of type |bg_gt|) and ``it`` is an iterator on a list of indexes 
of 
-already existing points. For instance, if one needs to add a new triangle in a 
3D 
+This is a template function, with ``pgt`` of type |bg_pgt| (basically a pointer
+to an instance of type |bg_gt|) and ``it`` is an iterator on a list of indexes 
of
+already existing points. For instance, if one needs to add a new triangle in a 
3D
 mesh, one needs to define first an array with the indexes of the three points::
 
   std::vector<bgeot::size_type> ind(3);
@@ -108,7 +108,7 @@ Specialized functions exist also::
   mymesh.add_prism(N, it);
   mymesh.add_prism_by_points(N, itp);
 
-The order of the points in the array of points is not important for simplices 
+The order of the points in the array of points is not important for simplices
 (except if you care about the orientation of your simplices). For other 
elements, it is important to respect the vertex order shown in 
:ref:`ud-fig-elem` (first order elements).
 
 .. _ud-fig-elem:
diff --git a/doc/sphinx/source/userdoc/computeL2H1.rst 
b/doc/sphinx/source/userdoc/computeL2H1.rst
index ab06b36..910bbb2 100644
--- a/doc/sphinx/source/userdoc/computeL2H1.rst
+++ b/doc/sphinx/source/userdoc/computeL2H1.rst
@@ -20,7 +20,7 @@ the different norms::
 where ``mim`` is a |gf_mim| used for the integration, ``mf`` is a |gf_mf| and
 describes the finite element method on which the solution is defined, ``U`` is 
the
 vector of values of the solution on each degree of freedom of ``mf`` and 
``region`` is an optional parameter which specify the mesh region on which the 
norm is computed. The size of
-``U`` should be ``mf.nb_dof()``. 
+``U`` should be ``mf.nb_dof()``.
 
 In order to compare two solutions, it is often simpler and faster to use the
 following function than to interpolate one |mf| on another::
diff --git a/doc/sphinx/source/userdoc/convect.rst 
b/doc/sphinx/source/userdoc/convect.rst
index 2b0f94b..d8ca138 100644
--- a/doc/sphinx/source/userdoc/convect.rst
+++ b/doc/sphinx/source/userdoc/convect.rst
@@ -6,17 +6,17 @@
 
 .. _ud-convect:
 
-A pure convection method 
+A pure convection method
 ========================
 
-A method to compute a pure convection is defined in the file 
+A method to compute a pure convection is defined in the file
 :file:`getfem/getfem_convect.h`. The call of the function is::
 
   getfem::convect(mf, U, mf_v, V, dt, nt, option = CONVECT_EXTRAPOLATION);
 
-where ``mf`` is a variable of type |gf_mf|, ``U`` is a vector which represent 
the 
-field to be convected, ``mf_v`` is a |gf_mf| for the velocity field, ``V`` is 
the 
-dof vector for the velocity field, ``dt`` is the pseudo time of convection and 
+where ``mf`` is a variable of type |gf_mf|, ``U`` is a vector which represent 
the
+field to be convected, ``mf_v`` is a |gf_mf| for the velocity field, ``V`` is 
the
+dof vector for the velocity field, ``dt`` is the pseudo time of convection and
 ``nt`` the number of iterations for the computation of characteristics. 
``option`` is an option for the boundary condition where there is a re-entrant 
convection. The possibilities are getfem::CONVECT_EXTRAPOLATION (extrapolation 
of the field on the nearest element) or getfem::CONVECT_UNCHANGED (no change of 
the value on the boundary).
 
 The method integrate the partial differential equation
@@ -27,21 +27,21 @@ The method integrate the partial differential equation
 
 on the time intervall :math:`[0, dt]`.
 
-The method used is of Galerkin-Characteristic kind. It is a very simple 
version 
-which is inconditionnally stable but rather dissipative. See [ZT1989]_ and 
also the Freefem++ documentation on convect 
+The method used is of Galerkin-Characteristic kind. It is a very simple version
+which is inconditionnally stable but rather dissipative. See [ZT1989]_ and 
also the Freefem++ documentation on convect
 command.
 
-The defined method works only if ``mf`` is a pure Lagrange finite element 
method 
+The defined method works only if ``mf`` is a pure Lagrange finite element 
method
 for the moment. The principle is to convect backward the finite element nodes 
by solving the ordinary differential equation:
 
 .. math::
 
    \frac{d X}{d t} = -V(X),
 
-with an initial condition corresponding to each node. This convection is made 
with ``nt`` steps. Then the solution is interploated on 
+with an initial condition corresponding to each node. This convection is made 
with ``nt`` steps. Then the solution is interploated on
 the convected nodes.
 
-In order to make the extrapolation not too expensive, the product 
:math:`dt\times V` 
+In order to make the extrapolation not too expensive, the product 
:math:`dt\times V`
 should not be too large.
 
 Note that this method can be used to solve convection dominant problems 
coupling it with a splitting scheme.
diff --git a/doc/sphinx/source/userdoc/gasm_high.rst 
b/doc/sphinx/source/userdoc/gasm_high.rst
index 9119120..6c8d4e2 100644
--- a/doc/sphinx/source/userdoc/gasm_high.rst
+++ b/doc/sphinx/source/userdoc/gasm_high.rst
@@ -42,7 +42,7 @@ A specific weak form language has been developed to describe 
the weak formulatio
 
   - Some constants: ``pi``, ``meshdim`` (the dimension of the current mesh), 
``qdim(u)`` and ``qdims(u)`` the dimensions of the variable ``u`` (the size for 
fixed size variables and the dimension of the vector field for FEM variables), 
``Id(n)`` the identity :math:`n\times n` matrix.
 
-  - Parentheses can be used to change the operations order in a standard way. 
For instance ``(1+2)*4`` or ``(u+v)*Test_u`` are valid expressions. 
+  - Parentheses can be used to change the operations order in a standard way. 
For instance ``(1+2)*4`` or ``(u+v)*Test_u`` are valid expressions.
 
   - The access to a component of a vector/matrix/tensor can be done by 
following a term by a left parenthesis, the list of components and a right 
parenthesis. For instance ``[1,1,2](3)`` is correct and will return ``2``. Note 
that indices are assumed to begin by 1 (even in C++ and with the python 
interface). A colon can replace the value of an index in a Matlab like syntax.
 
@@ -60,7 +60,7 @@ A specific weak form language has been developed to describe 
the weak formulatio
 
   - A certain number of linear and nonlinear operators (``Trace``, ``Norm``, 
``Det``, ``Deviator``, ``Contract``, ...). The nonlinear operators cannot be 
applied to test functions.
 
-  - ``Diff(expression, variable)``: The possibility to explicit differentiate 
an expression with respect to a variable (symbolic differentiation). 
+  - ``Diff(expression, variable)``: The possibility to explicit differentiate 
an expression with respect to a variable (symbolic differentiation).
 
   - ``Diff(expression, variable, direction)``: computes the derivative of 
``expression`` with respect to ``variable`` in the direction ``direction``.
 
@@ -83,7 +83,7 @@ The weak formulation for the Poisson problem on a domain 
:math:`\Omega`
 
   -\mbox{div } \nabla u = f, \mbox{ in } \Omega,
 
-with Dirichlet boundary conditions :math:`u = 0` on :math:`\partial\Omega` is 
classically 
+with Dirichlet boundary conditions :math:`u = 0` on :math:`\partial\Omega` is 
classically
 
 
 .. math::
@@ -123,20 +123,19 @@ Another classical equation is linear elasticity:
 
 for :math:`u` a vector field and :math:`\sigma(u) = \lambda \mbox{div } u + 
\mu (\nabla u + (\nabla u)^T)` when isotropic linear elasticity is considered. 
The corresponding assembly string to describe the weak formulation can be 
written::
 
-  (lambda*Trace(Grad_u)*Id(qdim(u)) + mu*(Grad_u+Grad_u')):Grad_Test_u - 
my_f.Test_u
+  "(lambda*Trace(Grad_u)*Id(qdim(u)) + mu*(Grad_u+Grad_u')):Grad_Test_u - 
my_f.Test_u"
 
-or:: 
+or::
 
-  lambda*Div_u*Div_Test_u + mu*(Grad_u + Grad_u'):Grad_Test_u - my_f.Test_u
+  "lambda*Div_u*Div_Test_u + mu*(Grad_u + Grad_u'):Grad_Test_u - my_f.Test_u"
 
 Here again, the coefficients ``lambda`` and ``mu`` can be given constants, or 
scalar field or explicit expression or even expression coming from some other 
variables in order to couples some problems. For instance, if the coefficients 
depends on a temperature field one can write::
 
-  my_f1(theta)*Div_u*Div_Test_u
-  + my_f2(theta)*(Grad_u + Grad_u'):Grad_Test_u - my_f.Grad_Test_u
+  "my_f1(theta)*Div_u*Div_Test_u + my_f2(theta)*(Grad_u + Grad_u'):Grad_Test_u 
- my_f.Grad_Test_u"
 
 where ``theta`` is the temperature which can be the solution to a Poisson 
equation::
 
-  Grad_theta.Grad_Test_theta - my_f*Grad_Test_theta
+  "Grad_theta.Grad_Test_theta - my_f*Grad_Test_theta"
 
 and ``my_f1`` and ``my_f2`` are some given functions. Note that in that case, 
the problem is nonlinear due to the coupling, even if the two functions  
``my_f1`` and ``my_f2`` are linear.
 
@@ -171,15 +170,15 @@ with ``model`` a previously define ``getfem::model`` 
object. In that case the va
 In that case, the variable and constant have to be added to the workspace. 
This can be done thanks to the following methods::
 
   workspace.add_fem_variable(name, mf, I, V);
-  
+
   workspace.add_fixed_size_variable(name, I, V);
 
   workspace.add_fem_constant(name, mf, V);
-  
+
   workspace.add_fixed_size_constant(name, V);
 
   workspace.add_im_data(name, imd, V);
-  
+
 where ``name`` is the variable/constant name (see in the next sections the 
restriction on possible names), ``mf`` is the ``getfem::mesh_fem`` object 
describing the finite element method, ``I`` is an object of class 
``gmm::sub_interval`` indicating the interval of the variable on the assembled 
vector/matrix and ``V`` is a ``getfem::base_vector`` being the value of the 
variable/constant. The last method add a constant defined on an ``im_data`` 
object ``imd`` which allows to store scalar/ve [...]
 
 
@@ -230,7 +229,7 @@ As a first example, if one needs to perform the assembly of 
a Poisson problem
 
   -\mbox{div } \nabla u = f, \mbox{ in } \Omega,
 
-the stiffness matrix is given 
+the stiffness matrix is given
 
 .. math::
 
@@ -238,7 +237,6 @@ the stiffness matrix is given
 
 and will be assembled by the following code::
 
-
   getfem::ga_workspace workspace;
   getfem::size_type nbdof = mf.nb_dof();
   getfem::base_vector U(nbdof);
@@ -421,7 +419,7 @@ A certain number of predefined scalar functions can be 
used. The exhaustive list
   - ``pos_part(t)`` (:math:`tH(t)`)
   - ``reg_pos_part(t, eps)`` (:math:`(t-eps/2-t^2/(2eps))H(t-eps) + 
t^2H(t)/(2eps)`)
   - ``neg_part(t)`` (:math:`-tH(-t)`), ``max(t, u)``, ``min(t, u)``
-     
+
 A scalar function can be applied to a scalar expression, but also to a tensor 
one. If is is applied to a tensor expression, is is applied componentwise and 
the result is a tensor with the same dimensions. For functions having two 
arguments (pow(t,u), min(t,u) ...) if two non-scalar arguments are passed, the 
dimension have to be the same. For instance "max([1;2],[0;3])" will return 
"[0;3]".
 
 
@@ -462,7 +460,7 @@ Binary operations
 
 A certain number of binary operations between tensors are available:
 
-  
+
     - ``+`` and ``-`` are the standard addition and subtraction of scalar, 
vector, matrix or tensors.
 
     - ``*`` stands for the scalar, matrix-vector, matrix-matrix or (fourth 
order tensor)-matrix multiplication.
@@ -480,24 +478,24 @@ A certain number of binary operations between tensors are 
available:
     - address@hidden stands for the tensor product.
 
     - ``Contract(A, i, B, j)`` stands for the contraction of tensors A and B 
with respect to the ith index of A and jth index of B. The first index is 
numbered 1. For instance ``Contract(V,1,W,1)`` is equivalent to ``V.W`` for two 
vectors ``V`` and ``W``.
-      
+
     - ``Contract(A, i, j, B, k, l)`` stands for the double contraction of 
tensors A and B with respect to indices i,j of A and indices k,l of B. The 
first index is numbered 1. For instance ``Contract(A,1,2,B,1,2)`` is equivalent 
to ``A:B`` for two matrices ``A`` and ``B``.
-      
+
 
 Unary operators
 ---------------
- 
+
   - ``-`` the unary minus operator: change the sign of an expression.
-  
+
   - ``'`` stands for the transpose of a matrix or line view of a vector. It a 
tensor ``A`` is of order greater than two,``A'`` denotes the inversion of the 
two first indices.
-  
+
   - ``Contract(A, i, j)`` stands for the contraction of tensor A with respect 
to its ith and jth indices. The first index is numbered 1. For instance, 
``Contract(A, 1, 2)`` is equivalent to ``Trace(A)`` for a matrix ``A``.
 
   - ``Swap_indices(A, i, j)`` exchange indices number i and j. The first index 
is numbered 1. For instance ``Swap_indices(A, 1, 2)`` is equivalent to ``A'`` 
for a matrix ``A``.
 
   - ``Index_move_last(A, i)`` move the index number i in order to be the last 
one. For instance, if ``A`` is a fourth order tensor :math:`A_{i_1i_2i_3i_4}`, 
then the result of ``Index_move_last(A, 2)`` will be the tensor 
:math:`B_{i_1i_3i_4i_2} = A_{i_1i_2i_3i_4}`. For a matrix, ``Index_move_last(A, 
1)`` is equivalent to ``A'``.
 
-    
+
 Parentheses
 -----------
 
@@ -530,17 +528,17 @@ Constant expressions
 --------------------
 
   - Floating points with standards notations (for instance ``3``, ``1.456``, 
``1E-6``)
-  - ``pi``: the constant Pi. 
+  - ``pi``: the constant Pi.
   - ``meshdim``: the dimension of the current mesh (i.e. size of geometrical 
nodes)
-  - ``timestep``: the main time step of the model on which this assembly 
string is evaluated (defined by ``model.set_time_step(dt)``). Do not work on 
pure workspaces. 
+  - ``timestep``: the main time step of the model on which this assembly 
string is evaluated (defined by ``model.set_time_step(dt)``). Do not work on 
pure workspaces.
   - ``Id(n)``: the identity matrix of size :math:`n\times n`. `n` should be an 
integer expression. For instance ``Id(meshdim)`` is allowed.
   - ``qdim(u)``: the total dimension of the variable ``u`` (i.e. the  size for 
fixed size variables and the total dimension of the vector/tensor field for FEM 
variables)
   - ``qdims(u)``: the dimensions of the variable ``u`` (i.e. the size for 
fixed size variables and the vector of dimensions of the vector/tensor field 
for FEM variables)
 
-Special expressions linked to the current position 
+Special expressions linked to the current position
 --------------------------------------------------
 
-  - ``X`` is the current coordinate on the real element (i.e. the position on 
the mesh of the current Gauss point on which the expression is evaluated), 
``X(i)`` is its i-th component. For instance ``sin(X(1)+X(2))`` is a valid 
expression on a mesh of dimension greater or equal to two. 
+  - ``X`` is the current coordinate on the real element (i.e. the position on 
the mesh of the current Gauss point on which the expression is evaluated), 
``X(i)`` is its i-th component. For instance ``sin(X(1)+X(2))`` is a valid 
expression on a mesh of dimension greater or equal to two.
 
   - ``Normal`` the outward unit normal vector to a boundary when integration 
on a boundary is performed.
 
@@ -704,7 +702,7 @@ is equivalent to::
   Grad_u
 
 for a varible ``u``.
-  
+
 .. _ud-gasm-high-transf:
 
 Interpolate transformations
@@ -714,10 +712,10 @@ The ``Interpolate`` operation allows to compute integrals 
between quantities whi
 
 In order to use this functionality, the user have first to declare to the 
workspace or to the model object an interpolate transformation which described 
the map between the current integration point and the point lying on the same 
mesh or on another mesh.
 
-Different kind of transformations can be described. Several kinds of 
transformations has been implemented. The first one, described hereafter is a 
transformation described by an expression. A second one corresponds to the 
raytracing contact detection (see 
:ref:`ud-model-contact-friction_raytrace_inter_trans`). Some other 
transformations (neighbour element and element extrapolation) are describe in 
the next sections. 
+Different kind of transformations can be described. Several kinds of 
transformations has been implemented. The first one, described hereafter is a 
transformation described by an expression. A second one corresponds to the 
raytracing contact detection (see 
:ref:`ud-model-contact-friction_raytrace_inter_trans`). Some other 
transformations (neighbour element and element extrapolation) are describe in 
the next sections.
 
 The transformation defined by an expression can be added to the workspace or 
the model thanks to the command::
- 
+
   add_interpolate_transformation_from_expression
     (workspace, transname, source_mesh, target_mesh, expr);
 
@@ -792,12 +790,12 @@ A specific transformation (see previous section) is 
defined in order to allows t
   (workspace, transname, my_mesh, std::map<size_type, size_type> &elt_corr);
 
 The map elt_corr should contain the correspondances between the elements where 
the transformation is to be applied and the respective elements where the 
extrapolation has to be made. On the element not listed in the map, no 
transformation is applied and the evaluation is performed normally on the 
current element.
-  
+
 The following functions allow to change the element correspondance of a 
previously added element extrapolation transformation::
 
   set_element_extrapolation_correspondance
   (model, transname, std::map<size_type, size_type> &elt_corr);
-  
+
   set_element_extrapolation_correspondance
   (workspace, transname, std::map<size_type, size_type> &elt_corr);
 
@@ -853,24 +851,24 @@ In some very special cases, it can be interesting to 
compute an integral on the
 
 where :math:`k(x,y)` is a given kernel, :math:`u` a quantity defined on 
:math:`\Omega_1` and  :math:`v` a quantity defined on :math:`\Omega_2`, 
eventually with  :math:`\Omega_1` and :math:`\Omega_2` the same domain. This 
can be interesting either to compute such an integral or to define an 
interaction term between two variables defined on two different domains.
 
-CAUTION: Of course, this kind of term have to be used with great care, since 
it naturally leads to fully populated stiffness or tangent matrices. 
+CAUTION: Of course, this kind of term have to be used with great care, since 
it naturally leads to fully populated stiffness or tangent matrices.
 
 
 The weak form language of |gf| furnishes a mechanism to compute such a term. 
First, the secondary domain has to be declared in the workspace/model with its 
integration methods. The addition of a standard secondary domain can be done 
with one of the two following functions::
 
   add_standard_secondary_domain(model, domain_name, mim, region);
-  
+
   add_standard_secondary_domain(workspace, domain_name, mim, region);
 
 where ``model`` or ``workspace`` is the model or workspace where the secondary 
domain has to be declared, ``domain_name`` is a string for the identification 
of this domain together with the mesh region and integration method, ``mim`` 
the integration method and ``region`` a mesh region. Note that with these 
standard secondary domains, the integration is done on the whole region for 
each element of the primary domain. It can be interesting to implement specific 
secondary domains restrictin [...]
 
-Once a secondary domain has been declared, it can be specified that a weak 
form language expression has to be assembled on the direct product of a current 
domain and a secondary domain, adding the name of the secondary domain to the 
``add_expression`` method of the workspace object or using 
``add_linear_twodomain_term``, ``add_nonlinear_twodomain_term`` or 
``add_twodomain_source_term`` functions:: 
+Once a secondary domain has been declared, it can be specified that a weak 
form language expression has to be assembled on the direct product of a current 
domain and a secondary domain, adding the name of the secondary domain to the 
``add_expression`` method of the workspace object or using 
``add_linear_twodomain_term``, ``add_nonlinear_twodomain_term`` or 
``add_twodomain_source_term`` functions::
 
   workspace.add_expression(expr, mim, region, derivative_order, 
secondary_domain)
   add_twodomain_source_term(model, mim, expr, region, secondary_domain)
   add_linear_twodomain_term(model, mim, expr, region, secondary_domain)
   add_nonlinear_twodomain_term(model, mim, expr, region, secondary_domain)
-  
+
 For the utilisation with the Python/Scilab/Matlab interface, see the 
documentation on ``gf_asm`` command and the ``model`` object.
 
 
@@ -946,7 +944,7 @@ For |gf| 5.1. When using a fem cut by a level-set (using 
fem_level_set or mesh_f
   Xfem_plus(Test_Grad_u)
   Xfem_plus(Test_Div_u)
   Xfem_plus(Test_Hess_u)
-  
+
   Xfem_minus(u)
   Xfem_minus(Grad_u)
   Xfem_minus(Div_u)
@@ -955,7 +953,7 @@ For |gf| 5.1. When using a fem cut by a level-set (using 
fem_level_set or mesh_f
   Xfem_minus(Test_Grad_u)
   Xfem_minus(Test_Div_u)
   Xfem_minus(Test_Hess_u)
-  
+
 which are only available when the evaluation (integration) is made on the 
curve/surface separating two zones of continuity, i.e. on the zero level-set of 
a considered level-set function (using a ``mesh_im_level_set`` object). For 
instance, a jump in the variable ``u`` will be given by::
 
   Xfem_plus(u)-Xfem_minus(u)
@@ -964,10 +962,10 @@ and the average by::
 
   (Xfem_plus(u)+Xfem_minus(u))/2
 
-  The value ``Xfem_plus(u)`` is the value of ``u`` on the side where the 
corresponding level-set function is positive and ``Xfem_minus(u)`` the value of 
``u`` on the side where the level-set function is negative.
+The value ``Xfem_plus(u)`` is the value of ``u`` on the side where the 
corresponding level-set function is positive and ``Xfem_minus(u)`` the value of 
``u`` on the side where the level-set function is negative.
+
+Additionally, note that, when integrating on a level-set with a 
``mesh_im_level_set`` object, ``Normal`` stands for the normal unit vector to 
the level-set in the direction of the gradient of the level-set function.
 
-  Additionally, note that, when integrating on a level-set with a 
``mesh_im_level_set`` object, ``Normal`` stands for the normal unit vector to 
the level-set in the direction of the gradient of the level-set function.
-  
 Storage of sub-expressions in a getfem::im_data object during assembly
 ----------------------------------------------------------------------
 
@@ -994,7 +992,7 @@ the tangent system).
 If before = 0 (default), the assignement is done after the assembly terms.
 
 Additionally, In a model, the method::
-  
+
   model.clear_assembly_assignments()
 
 allows to cancel all the assembly assignments previously added.
diff --git a/doc/sphinx/source/userdoc/index.rst 
b/doc/sphinx/source/userdoc/index.rst
index 3fd75ba..801659e 100644
--- a/doc/sphinx/source/userdoc/index.rst
+++ b/doc/sphinx/source/userdoc/index.rst
@@ -37,7 +37,7 @@ User Documentation
    model_nonlinear_elasticity
    model_plasticity_small_strain
    model_ALE_rotating
-   
+
    appendixA
    appendixB
 
diff --git a/doc/sphinx/source/userdoc/interMM.rst 
b/doc/sphinx/source/userdoc/interMM.rst
index 556efdc..1600a29 100644
--- a/doc/sphinx/source/userdoc/interMM.rst
+++ b/doc/sphinx/source/userdoc/interMM.rst
@@ -14,21 +14,21 @@ Once a solution has been computed, it is quite easy to 
extract any quantity of i
 Basic interpolation
 *******************
 
-The file :file:`getfem/getfem_interpolation.h` defines the function 
-``getfem::interpolation(...)`` to interpolate a solution from a given 
mesh/finite 
+The file :file:`getfem/getfem_interpolation.h` defines the function
+``getfem::interpolation(...)`` to interpolate a solution from a given 
mesh/finite
 element method on another mesh and/or another Lagrange finite element method::
 
   getfem::interpolation(mf1, mf2, U, V, extrapolation = 0);
 
-where ``mf1`` is a variable of type |gf_mf| and describes the finite element 
-method on which the source field ``U`` is defined, ``mf2`` is the finite 
element 
-method on which ``U`` will be interpolated. ``extrapolation`` is an optional 
-parameter. The values are ``0`` not to allow the extrapolation, ``1`` for an 
-extrapolation of the exterior points near the boundary and ``2`` for the 
+where ``mf1`` is a variable of type |gf_mf| and describes the finite element
+method on which the source field ``U`` is defined, ``mf2`` is the finite 
element
+method on which ``U`` will be interpolated. ``extrapolation`` is an optional
+parameter. The values are ``0`` not to allow the extrapolation, ``1`` for an
+extrapolation of the exterior points near the boundary and ``2`` for the
 extrapolation of all exterior points (could be expensive).
 
 
-The dimension of ``U`` should be a multiple of ``mf1.nb_dof()``, and the 
+The dimension of ``U`` should be a multiple of ``mf1.nb_dof()``, and the
 interpolated data ``V`` should be correctly sized (multiple of 
``mf2.nb_dof()``).
 
 ... important::
diff --git a/doc/sphinx/source/userdoc/intro.rst 
b/doc/sphinx/source/userdoc/intro.rst
index 4b609a6..bfa7000 100644
--- a/doc/sphinx/source/userdoc/intro.rst
+++ b/doc/sphinx/source/userdoc/intro.rst
@@ -45,7 +45,7 @@ are just some parameters that can be changed very easily, thus
 allowing a large spectrum of experimentations. Numerous examples are
 available in the ``tests`` directory of the distribution.
 
-|gf| has only a (very) experimental meshing procedure (and produces regular 
meshes), hence it is generally 
+|gf| has only a (very) experimental meshing procedure (and produces regular 
meshes), hence it is generally
 necessary to import meshes. Imports formats currently known by |gf|
 are |gid|, |gmsh| and *emc2* mesh files. However, given a mesh, it
 is possible to refine it automatically.
diff --git a/doc/sphinx/source/userdoc/model.rst 
b/doc/sphinx/source/userdoc/model.rst
index 2602803..628a8f2 100644
--- a/doc/sphinx/source/userdoc/model.rst
+++ b/doc/sphinx/source/userdoc/model.rst
@@ -37,7 +37,7 @@ The kernel of the model description is contained in the file
 
 .. toctree::
    :maxdepth: 2
-  
+
    model_object
    model_generic_assembly
    model_generic_elliptic
@@ -57,4 +57,4 @@ The kernel of the model description is contained in the file
    model_time_integration
    model_contact_friction
    model_contact_friction_large_sliding
-  
+
diff --git a/doc/sphinx/source/userdoc/model_ALE_rotating.rst 
b/doc/sphinx/source/userdoc/model_ALE_rotating.rst
index 49cb811..833ce78 100644
--- a/doc/sphinx/source/userdoc/model_ALE_rotating.rst
+++ b/doc/sphinx/source/userdoc/model_ALE_rotating.rst
@@ -61,7 +61,7 @@ We have then
 .. math::
  \dot{r}(t,X) = \dot{\theta}AR(t)X
 
-If :math:`\varphi(t, X)` is the deformation of the body which maps the 
reference configuration :math:`\Omega^0` to the deformed configuration 
:math:`\Omega_t` at time :math:`t`, the ALE description consists in the 
decomposition of the deformation of the cylinder in 
+If :math:`\varphi(t, X)` is the deformation of the body which maps the 
reference configuration :math:`\Omega^0` to the deformed configuration 
:math:`\Omega_t` at time :math:`t`, the ALE description consists in the 
decomposition of the deformation of the cylinder in
 
 .. math::
    \varphi(t, X) = (\tau(t) \circ \bar{\varphi}(t) \circ r(t))(X) = 
\bar{\varphi}(t, r(t, X)) + Z(t)
@@ -72,7 +72,7 @@ With :math:`\bar{X} = R(t)X` the new considered deformation is
   \bar{\varphi}(t,\bar{X}) = \varphi(X) - Z(t)
 
 
-Thanks to the rotation symmetry of the reference configuration 
:math:`\Omega^0:`, we note that :math:`\bar{\Omega}^0 = r(t, \Omega^0)` is 
independant of :math:`t` and will serve as the new reference configuration. 
This is illustrated in the following figure: 
+Thanks to the rotation symmetry of the reference configuration 
:math:`\Omega^0:`, we note that :math:`\bar{\Omega}^0 = r(t, \Omega^0)` is 
independant of :math:`t` and will serve as the new reference configuration. 
This is illustrated in the following figure:
 
 .. _ud-fig-rotating_cylinder_conf:
 
@@ -111,8 +111,8 @@ This should not be forgotten that a correction has to be 
provided for each evolv
   \Frac{\partial u}{\partial t} = \Frac{\partial \bar{u}}{\partial t} + 
\dot{\theta} (I_d + \nabla \bar{u}) A \bar{X} + \dot{Z}(t),
 
   \Frac{\partial^2 u}{\partial t^2} = \Frac{\partial^2 \bar{u}}{\partial t^2} 
+ 2\dot{\theta} \nabla\Frac{\partial \bar{u}}{\partial t}A \bar{X} +  
\dot{\theta}^2\mbox{div}(((I_d + \nabla\bar{u})A \bar{X}) \otimes (A \bar{X}) ) 
 + \ddot{\theta} (I_d + \nabla\bar{u}) A \bar{X}  + \ddot{Z}(t).
-  
-  
+
+
 
 
 
@@ -125,7 +125,7 @@ Assuming :math:`\rho^0` the density in the reference 
configuration having a rota
    \int_{\Omega^0} \rho^0 \Frac{\partial^2 u}{\partial t^2}\cdot vdX =
 
    \int_{\bar{\Omega}^0} \rho^0 \left[\Frac{\partial^2 \bar{u}}{\partial t^2} 
+ 2\dot{\theta} \nabla\Frac{\partial \bar{u}}{\partial t}A \bar{X} +  
\dot{\theta}^2\mbox{div}(((I_d + \nabla\bar{u})A \bar{X}) \otimes (A \bar{X}) ) 
 + \ddot{\theta} (I_d + \nabla\bar{u}) A \bar{X}  + \ddot{Z}(t) \right] \cdot 
\bar{v} d\bar{X}.
-   
+
 The third term in the right hand side can be integrated by part as follows:
 
 .. math::
@@ -158,7 +158,7 @@ To be adapted  ::
 
   ind = getfem::brick_name(parmeters);
 
-where ``parameters`` are the parameters ... 
+where ``parameters`` are the parameters ...
 
 
 
@@ -209,7 +209,7 @@ the displacement on the intermediary configuration, then it 
is easy to check tha
 
 .. math::
    \Frac{\partial \varphi}{\partial t} = \Frac{\partial \bar{u}}{\partial t} - 
\nabla \bar{u} \dot{Z}
- 
+
    \Frac{\partial^2 \varphi}{\partial t^2} = \Frac{\partial^2 
\bar{u}}{\partial t^2} - \nabla\Frac{\partial \bar{u}}{\partial t}\dot{Z} + 
\Frac{\partial^2 \bar{u}}{\partial \dot{Z}^2} - \nabla\bar{u}\ddot{Z}.
 
 
@@ -226,7 +226,7 @@ Assuming :math:`\rho^0` the density in the reference being 
invariant with the co
    \int_{\bar{\Omega}^{0}} \rho^0 \left[\Frac{\partial^2 \bar{u}}{\partial 
t^2} - 2\nabla\Frac{\partial \bar{u}}{\partial t}\dot{Z} - 
\nabla\bar{u}\ddot{Z}\right]\cdot \bar{v}  - \rho^0 
(\nabla\bar{u}\dot{Z}).(\nabla\bar{v}\dot{Z}) d\bar{X} + \int_{\partial 
\bar{\Omega}^0} \rho^0 (\nabla\bar{u}\dot{Z}).\bar{v}(\dot{Z}.\bar{N}) 
d\bar{\Gamma},
 
 where :math:`\bar{N}` is the outward unit normal vector on :math:`\partial 
\bar{\Omega}^0`. Note that the last term vanishes on :math:`(\alpha, \beta) 
\times \partial \omega^0` but not necessarily on :math:`\{\alpha\} \times 
\omega^0` and :math:`\{\beta\} \times \omega^0`.
-   
+
 
 
 
diff --git a/doc/sphinx/source/userdoc/model_Nitsche.rst 
b/doc/sphinx/source/userdoc/model_Nitsche.rst
index 7575ca5..b0bee0e 100644
--- a/doc/sphinx/source/userdoc/model_Nitsche.rst
+++ b/doc/sphinx/source/userdoc/model_Nitsche.rst
@@ -33,7 +33,7 @@ If additionally a mixed incompressibility brick is added with 
a variable :math:`
 :math:`G = \sigma(u)n - pn`
 
 In order to allow a generic implementation in which the brick imposing 
Nitsche's method will work for every partial differential term applied to the 
concerned variables, each brick adding a partial differential term to a model 
is required to give its expression via an assembly string (weak form language).
- 
+
 These expressions are utilized in a special method of the model object::
 
   expr = md.Neumann_term(variable, region)
@@ -42,9 +42,9 @@ which allows to automatically derive an expression for the 
sum of all Neumann te
 Of course it is required that all volumic bricks were added to the model prior 
to the call of this method.
 The derivation of the Neumann term works only for second order partial 
differential equations.
 A generic implementation for higher order pde would be more complicated.
-   
 
-Generic Nitsche's method for a Dirichlet condition 
+
+Generic Nitsche's method for a Dirichlet condition
 ++++++++++++++++++++++++++++++++++++++++++++++++++
 
 Assume that the variable :math:`u` is considered and that one wants to 
prescribe the condition
@@ -71,7 +71,7 @@ The bricks adding a Dirichlet condition with Nitsche's method 
to a model are the
 
   getfem::add_Dirichlet_condition_with_Nitsche_method
      (model &md, const mesh_im &mim, const std::string &varname,
-      const std::string &Neumannterm, 
+      const std::string &Neumannterm,
       const std::string &gamma0name, size_type region,
       scalar_type theta = scalar_type(1),
       const std::string &dataname = std::string());
@@ -127,7 +127,7 @@ even for nonlinear problems. Returns the brick index in the 
model.
 
   getfem::add_generalized_Dirichlet_condition_with_Nitsche_method
      (model &md, const mesh_im &mim, const std::string &varname,
-      const std::string &Neumannterm, 
+      const std::string &Neumannterm,
       const std::string &gamma0name, size_type region, scalar_type theta,
       const std::string &dataname, const std::string &Hname);
 
@@ -160,7 +160,7 @@ or described on a scalar fem. Returns the brick index in 
the model.
 
 .. _nitsche_contact_small_def_section:
 
-Generic Nitsche's method for contact with friction condition 
+Generic Nitsche's method for contact with friction condition
 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
 We describe here the use of Nitsche's method to prescribe a contact with 
Coulomb friction condition in the small deformations framework. This 
corresponds to a weak integral contact condition which as some similarity with 
the ones which use Lagrange multipliers describe in the corresponding section, 
see :ref:`weak_integral_contact_section`
diff --git a/doc/sphinx/source/userdoc/model_bilaplacian.rst 
b/doc/sphinx/source/userdoc/model_bilaplacian.rst
index 40c20e9..e8128ce 100644
--- a/doc/sphinx/source/userdoc/model_bilaplacian.rst
+++ b/doc/sphinx/source/userdoc/model_bilaplacian.rst
@@ -18,7 +18,7 @@ The following function ::
                              region = size_type(-1));
 
 adds a bilaplacian brick on the variable `varname` and on the mesh region 
`region`. This represent a term :math:`\Delta(D \Delta u)`. where :math:`D(x)` 
is a coefficient determined by `dataname` which could be constant or described 
on a f.e.m. The corresponding weak form is :math:`\int D(x)\Delta u(x) \Delta 
v(x) dx`.
- 
+
 
 For the Kirchhoff-Love plate model, the weak form is a bit different (and more 
stable than the previous one). the function to add that term is ::
 
@@ -67,7 +67,7 @@ And a Dirichlet condition on the normal derivative can be 
prescribed thanks to t
        R_must_be_derivated = false);
 
 These bricks add a Dirichlet condition on the normal derivative of the variable
-`varname` and on the mesh region `region` (which should be a boundary. 
+`varname` and on the mesh region `region` (which should be a boundary.
 The general form is :math:`\int \partial_n u(x)v(x) = \int r(x)v(x) \forall v`
 where :math:`r(x)` is the right hand side for the Dirichlet condition (0 for
 homogeneous conditions) and :math:`v` is in a space of multipliers
@@ -91,6 +91,6 @@ The test program :file:`bilaplacian.cc` is a good example of 
the use of the prev
 
 
 
-   
+
 
 
diff --git a/doc/sphinx/source/userdoc/model_contact_friction.rst 
b/doc/sphinx/source/userdoc/model_contact_friction.rst
index 05c528c..eb83042 100644
--- a/doc/sphinx/source/userdoc/model_contact_friction.rst
+++ b/doc/sphinx/source/userdoc/model_contact_friction.rst
@@ -68,13 +68,13 @@ Suppose now that you approximate a linearized elasticity 
problem submitted to co
 .. math::
 
   u_N(a_i) = (B_N U)_i,
- 
+
   (\dot{u}_T(a_i))_k = (B_T \dot{U})_{(d-1)(i-1)+k},
 
 where :math:`d` is the dimension of the domain and :math:`k = 1..d-1`. The 
expression of the elasticity problem with contact with friction can be written 
as
 
 .. math::
- 
+
   K U = L + B_N^T \lambda_N + B_T^T \lambda_T,
 
   -\frac{1}{r\alpha_i}(\lambda_N^i - P_{]-\infty, 0]}(\lambda_N^i - \alpha_i r 
((B_N U)_i - \text{gap}_i))) = 0, ~~ i = 1..N_c,
@@ -102,7 +102,7 @@ where :math:`a_i`, :math:`~~i=1..N_c` are the finite 
element nodes corresponding
 
 .. math::
 
-  \int_{\Gamma_c} (\mu_N^h - \lambda_N^h) (u_N - \text{gap}) d\Gamma \ge 0 ~~ 
\forall \mu_N^h \in \Lambda_N^h. 
+  \int_{\Gamma_c} (\mu_N^h - \lambda_N^h) (u_N - \text{gap}) d\Gamma \ge 0 ~~ 
\forall \mu_N^h \in \Lambda_N^h.
 
 In that case, the component :math:`\lambda_N^i` is a contact stress 
(:math:`N/m^2`) and the matrix :math:`B_N` can be written
 
@@ -118,10 +118,10 @@ The matrix :math:`B_T` can also be written in a similar 
way. The friction condit
 
 where :math:`\Lambda_T^h({\mathscr F}\lambda_N^h)` is the discrete set of 
admissible friction stress.
 
-Finally, the expression of the direct nodal contact condition are recovered 
+Finally, the expression of the direct nodal contact condition are recovered
 
 .. math::
- 
+
   K U = L + B_N^T \lambda_N + B_T^T \lambda_T,
 
   -\frac{1}{r\alpha_i}(\lambda_N^i - P_{]-\infty, 0]}(\lambda_N^i - \alpha_i r 
((B_N U)_i - \text{gap}_i))) = 0, ~~ i = 1..N_c,
@@ -196,7 +196,7 @@ where :math:`a(\cdot, \cdot)` and :math:`\ell(v)` represent 
the remaining parts
   ~~~~~~ \displaystyle +\int_{\Gamma_c}({\mathscr F} 
D_{\rho}P_{B(\rho)}(\lambda^h_T - 
r\alpha(u^h_T-w^h_T))\delta_{u_N})\cdot\mu^h_T d\Gamma \\
   ~~~~~~ \displaystyle -\int_{\Gamma_c}(\frac{\mathscr F}{r} 
D_{\rho}P_{B(\rho)}(\lambda^h_T - 
r\alpha(u^h_T-w^h_T))\delta_{\lambda_N})\cdot\mu^h_T d\Gamma = \cdots ~~~ 
\forall \mu^h \in W^h,
   \end{array}\right.
-  
+
 where :math:`H(\cdot)` is the Heaviside function (0 for a negative argument 
and 1 for a non-negative argument), :math:`D_xP_{B(\rho)}(x)` and 
:math:`D_{\rho}P_{B(\rho)}(x)` are the derivatives of the projection on 
:math:`B(\rho)` (assumed to vanish for :math:`\rho \le 0`) and 
:math:`\delta_{\lambda}` and :math:`\delta_{u}` are the unknown corresponding 
to the tangent problem.
 
 
@@ -336,7 +336,7 @@ The contact condition is prescribed thank to a multiplier
 size should be the number of rows of ``B_T``.
 The parameter ``dataname_friction_coeff`` describes the friction
 coefficient. It could be a scalar or a vector describing the
-coefficient on each contact condition. 
+coefficient on each contact condition.
 The augmentation parameter ``r`` should be chosen in a range of acceptable 
values
 (see Getfem user documentation). ``dataname_gap`` is an
 optional parameter representing the initial gap. It can be a single value
@@ -369,7 +369,7 @@ are 'x', 'y' in 2D and 'x', 'y', 'z' in 3D. For instance, 
if the rigid
 obstacle correspond to :math:`z \le 0`, the corresponding signed distance will
 be simply 'z'. ``multname_n`` should be a fixed size variable whose size is
 the number of degrees of freedom on boundary ``region``. It represents the
-contact equivalent nodal forces. 
+contact equivalent nodal forces.
 The augmentation parameter ``r`` should be chosen in a
 range of acceptable values (close to the Young modulus of the elastic
 body, see Getfem user documentation). 1 for the non-symmetric Alart-Curnier 
augmented Lagrangian, 2 for the symmetric one, 3 for the unsymmetric method 
based on augmented multipliers.
@@ -392,7 +392,7 @@ are 'x', 'y' in 2D and 'x', 'y', 'z' in 3D. For instance, 
if the rigid
 obstacle correspond to :math:`z \le 0`, the corresponding signed distance will
 be simply 'z'. ``multname_n`` should be a fixed size variable whose size is
 the number of degrees of freedom on boundary ``region``. It represents the
-contact equivalent nodal forces. 
+contact equivalent nodal forces.
 ``multname_t`` should be a fixed size variable whose size is
 the number of degrees of freedom on boundary ``region`` multiplied by
 :math:`d-1` where :math:`d` is the domain dimension. It represents the
diff --git a/doc/sphinx/source/userdoc/model_contact_friction_large_sliding.rst 
b/doc/sphinx/source/userdoc/model_contact_friction_large_sliding.rst
index ec3a981..170a6ba 100644
--- a/doc/sphinx/source/userdoc/model_contact_friction_large_sliding.rst
+++ b/doc/sphinx/source/userdoc/model_contact_friction_large_sliding.rst
@@ -28,7 +28,7 @@ In order to incorporate the contact detection in the 
high-level generic assembly
    :align: center
    :scale: 45
 
-The slave surface is the "contactor" and the master one the "target". Rigid 
obstacle are also considered. They are always master surfaces.  The basic rule 
is that the contact is considered between a slave surface and a master one. 
However, the multi-contact frame object and the |gf| bricks allow multi-contact 
situations, including contact between two master surfaces, self-contact of a 
master surface and an arbitrary number of slave and master surfaces. 
+The slave surface is the "contactor" and the master one the "target". Rigid 
obstacle are also considered. They are always master surfaces.  The basic rule 
is that the contact is considered between a slave surface and a master one. 
However, the multi-contact frame object and the |gf| bricks allow multi-contact 
situations, including contact between two master surfaces, self-contact of a 
master surface and an arbitrary number of slave and master surfaces.
 
 Basically, in order to detect the contact pairs, Gauss points or f.e.m. nodes 
of slave surfaces are projected on master surfaces (see  
:ref:`figure<ud-fig-masterslave>`). If self-contact is considered, Gauss points 
or f.e.m. nodes of master surface are also projected on master surfaces.
 
@@ -53,13 +53,13 @@ The raytracing transformation is added without any slave or 
master contact bound
 where ``dispname`` is the variable name which represent the displacement on 
that contact
 boundary. The difference between master and slave contact boundary is that the 
contact detection is to be performed starting from a slave or master boundary 
toward a master boundary. The contact detection is not performed toward a slave 
boundary. Consequently, only the influence boxes of the elements of the master 
surfaces are computed and stored.
 
-It is also possible to add a rigid obstacle (considered as a master surface) 
thanks to the function:: 
+It is also possible to add a rigid obstacle (considered as a master surface) 
thanks to the function::
 
   add_rigid_obstacle_to_raytracing_transformation(model &md,
              const std::string &transname,
              const std::string &expr, size_type N)
 
-where ``expr`` is the expression of a signed distance to the obstacle using 
the syntax of the weak form language (``X`` being the current position, 
``X(0)``, ``X(1)`` ... the corresponding components). For instance an 
expression ``X(0) + 5`` will correspond to a flat obstacle lying on the right 
of the position ``-5`` of the first coordinate. Be aware that the expression 
have to be close to a signed distance, which in particular means that the 
gradient norm have to be close to 1. 
+where ``expr`` is the expression of a signed distance to the obstacle using 
the syntax of the weak form language (``X`` being the current position, 
``X(0)``, ``X(1)`` ... the corresponding components). For instance an 
expression ``X(0) + 5`` will correspond to a flat obstacle lying on the right 
of the position ``-5`` of the first coordinate. Be aware that the expression 
have to be close to a signed distance, which in particular means that the 
gradient norm have to be close to 1.
 
 In order to distinguish between non-contact situations and the occurence of a 
contact with another deformable body or with a rigid obstacle, the 
transformation returns an integer identifier which can be used by the 
`Interpolate_filter` command of the weak form language (see 
:ref:`ud-gasm-high-transf`). The different values:
 
@@ -116,12 +116,12 @@ Some details on the algorithm:
     of contact pairs. The influence boxes are stored in a region tree object
     in order to find the boxes containing a point with an algorithm having
     a mean complexity in :math:`O(log(N))`.
-  
+
   - **What is a potential contact pair.** A potential contact pair is a pair
     slave point - master element face which will be investigated.
     The projection of the slave point on the master surface will be done
     and criteria will be applied.
- 
+
   - **Projection algorithm.** The projection of the slave point onto a
     master element face is done by a parametrization of the surface on the
     reference element via the geometric transformation and the displacement
@@ -153,7 +153,7 @@ The list of criteria:
     When Newton's algorithms (and BFGS one for projection) used to compute the
     projection/raytrace of the slave point on the master element surface
     fails to converge, the pair is not considered. A warning is generated.
-    
+
   - **Criterion 3 : the projected point should be inside the element.**
     The slave point is projected on the surface of the master element
     without the constraint to remain inside the face
@@ -207,7 +207,7 @@ We denote by :math:`x_i = \varphi^h(X_i)` the corresponding 
node on the deformed
 
   g_i = n_y . (\varphi^h(X_i) - \varphi^h(Y_i)) = \|\varphi^h(X_i) - 
\varphi^h(Y_i)\| \text{Sign}(n_y . (\varphi^h(X_i) - \varphi^h(Y_i))),
 
-where :math:`n_y` is the outward unit normal vector of the master surface at 
:math:`y`. 
+where :math:`n_y` is the outward unit normal vector of the master surface at 
:math:`y`.
 
 Considering only stationnary rigid obstacles and applying the principle of 
Alart-Curnier augmented Lagrangian [AL-CU1991]_, the problem with nodal contact 
with friction condition can be expressed as follows in an unsymmetric version 
(see [renard2013]_ for the linear elasticity case)
 
@@ -246,7 +246,7 @@ The following nonlinear operators are defined in the weak 
form language (see :re
     .. math::
 
       \partial_{u} n_{trans}[\delta u] = -(I - n_{trans}\otimes n_{trans})(I+ 
\nabla u)^{-T}(\nabla \delta u)^T n_{trans}
-  
+
       \partial_{n} n_{trans}[\delta n] = \Frac{(I+ \nabla u)^{-T}\delta n - 
n_{trans}(n_{trans}\cdot \delta n)}{\|(I+\nabla u)^{-T} n\|}
 
   - ``Coulomb_friction_coupled_projection(lambda, n, Vs, g, f, r)``
@@ -271,7 +271,7 @@ The following nonlinear operators are defined in the weak 
form language (see :re
       \partial_q P_{B(n,\tau)}(q) =
       \left\{\begin{array}{cl}
       0 & \mbox{for } \tau \le 0 \\
-      \mathbf{T}_n & \mbox{for } \|q_{_T}\| \le \tau \\ 
+      \mathbf{T}_n & \mbox{for } \|q_{_T}\| \le \tau \\
       \Frac{\tau}{\|q_{_T}\|}
       \left(\mathbf{T}_n - \Frac{q_{_T}}{\|q_{_T}\|}\otimes 
\Frac{q_{_T}}{\|q_{_T}\|}
       \right) & \mbox{otherwise }
@@ -305,7 +305,7 @@ The following nonlinear operators are defined in the weak 
form language (see :re
       \left|\begin{array}{l} \partial_n P_{B(n,\tau)}
       +\partial_{\tau} P_{B(n,\tau)} \otimes \partial_n \tau \\
       \hspace*{3em}+H(-\lambda\cdot n - r\,g) ~
-      \left(n \otimes \lambda - 
+      \left(n \otimes \lambda -
       (2~\lambda\cdot n + r\,g)~n \otimes n +
       (\lambda\cdot n + r\,g)~\mathbf{I}\right),
       \end{array}\right.
@@ -334,7 +334,7 @@ Add of the brick::
 
   indbrick = add_integral_large_sliding_contact_brick_raytracing
     (model &md, const std::string &dataname_r,
-     scalar_type release_distance, 
+     scalar_type release_distance,
      const std::string &dataname_friction_coeff = "0",
      const std::string &dataname_alpha = "1");
 
diff --git a/doc/sphinx/source/userdoc/model_continuation.rst 
b/doc/sphinx/source/userdoc/model_continuation.rst
index 221654c..8ea2b2d 100644
--- a/doc/sphinx/source/userdoc/model_continuation.rst
+++ b/doc/sphinx/source/userdoc/model_continuation.rst
@@ -21,7 +21,7 @@ written in the form
    F(U) = 0.
 
 In what follows, we shall suppose that the model depends on an additional 
scalar
-parameter :math:`\lambda` so that :math:`F(U) = F(U, \lambda)`. 
+parameter :math:`\lambda` so that :math:`F(U) = F(U, \lambda)`.
 
 Numerical continuation
 ++++++++++++++++++++++
@@ -243,14 +243,14 @@ Then :math:`\tau_{\mathrm{BP}}(Y(s))` changes its sign at 
:math:`s = \bar{s}`.
 Obviously, if one takes :math:`B`, :math:`C` and :math:`d` randomly, it is
 highly possible that they satisfy the requirements above. Consequently, the
 numerical continuation method is able to detect bifurcation points by
-taking the vectors :math:`Y` and :math:`T` supplied by the correction at each 
+taking the vectors :math:`Y` and :math:`T` supplied by the correction at each
 continuation step and monitoring the signs of :math:`\tau_{\mathrm{BP}}`.
 
 Once a bifurcation point :math:`\bar{Y}` is detected by a sign change
 :math:`\tau_{\mathrm{BP}}(Y_{j}) \tau_{\mathrm{BP}}(Y_{j+1}) < 0`, it can be
 approximated more precisely by the predictor-corrector steps described above
 with a special step-length adaptation (see Section 8.1 in [Al-Ge1997]_). 
Namely,
-one can take the subsequent step lengths as 
+one can take the subsequent step lengths as
 
    .. math::
 
@@ -427,8 +427,8 @@ Optionally, parametrisation by a vector datum is then 
declared by::
 
   S.set_parametrised_data_names(initdata_name, finaldata_name, 
currentdata_name);
 
-Here, the data names ``initdata_name`` and ``finaldata_name`` should represent 
-:math:`P^{0}` and :math:`P^{1}`, respectively. Under ``currentdata_name``, the 
+Here, the data names ``initdata_name`` and ``finaldata_name`` should represent
+:math:`P^{0}` and :math:`P^{1}`, respectively. Under ``currentdata_name``, the
 values of :math:`P(\lambda)` have to be stored, that is, actual values of the
 datum the model depends on.
 
@@ -439,7 +439,7 @@ Next, the continuation is initialised by::
 where ``U`` should be a solution for the value of the parameter :math:`\lambda`
 equal to ``lambda`` so that :math:`Y_{0}=` (\ ``U``\ ,\ ``lambda``\ ). During
 this initialisation, an initial unit tangent :math:`T_{0}` corresponding to
-:math:`Y_{0}` is computed in accordance with the sign of the initial value 
+:math:`Y_{0}` is computed in accordance with the sign of the initial value
 ``T_lambda``, and it is returned in ``T_U``, ``T_lambda``. Moreover, ``h`` is
 set to the initial step size ``h_init``.
 
diff --git a/doc/sphinx/source/userdoc/model_dirichlet.rst 
b/doc/sphinx/source/userdoc/model_dirichlet.rst
index 0013ef0..082cc96 100644
--- a/doc/sphinx/source/userdoc/model_dirichlet.rst
+++ b/doc/sphinx/source/userdoc/model_dirichlet.rst
@@ -99,17 +99,17 @@ If `dataname` is ommited, an homogeneous Dirichlet 
condition is applied. If `dat
 Generalized Dirichlet condition brick
 -------------------------------------
 
-The generalized Dirichlet condition is a boundary condition of a vector field 
u of 
+The generalized Dirichlet condition is a boundary condition of a vector field 
u of
 the type
 
 .. math::
 
    H u  = r
 
-where :math:`H` is a matrix field. The functions adding the corresponding 
bricks 
-are similar to the ones of the standard Dirichlet condition except that they 
need 
-the supplementary parameter `Hname` which gives the name of the data 
corresponding 
-to :math:`H`. This data can be a matrix field described on a scalar fem or a 
+where :math:`H` is a matrix field. The functions adding the corresponding 
bricks
+are similar to the ones of the standard Dirichlet condition except that they 
need
+the supplementary parameter `Hname` which gives the name of the data 
corresponding
+to :math:`H`. This data can be a matrix field described on a scalar fem or a
 constant matrix. ::
 
 
diff --git a/doc/sphinx/source/userdoc/model_generic_assembly.rst 
b/doc/sphinx/source/userdoc/model_generic_assembly.rst
index 64af4fe..31f5b08 100644
--- a/doc/sphinx/source/userdoc/model_generic_assembly.rst
+++ b/doc/sphinx/source/userdoc/model_generic_assembly.rst
@@ -18,7 +18,7 @@ A mean to add a term either on one variable or on several 
ones is to directly us
    size_type getfem::add_nonlinear_term(md, mim, expr,
                          region = -1, is_sym = false, is_coercive = false);
 
-This adds a brick to the model ``md``, using the integration method ``mim``, 
the assembly string ``expr`` on the mesh region ``region``. If the result is 
symmetric, you can specify it on the 5th argument and if it is coercive on the 
6th argument. The latter indications of symmetry and coercivness are used to 
determine the right linear solver. If you are not so sure, it is preferable not 
to indicate anything. 
+This adds a brick to the model ``md``, using the integration method ``mim``, 
the assembly string ``expr`` on the mesh region ``region``. If the result is 
symmetric, you can specify it on the 5th argument and if it is coercive on the 
6th argument. The latter indications of symmetry and coercivness are used to 
determine the right linear solver. If you are not so sure, it is preferable not 
to indicate anything.
 
 However, this brick consider that the expression is nonlinear. This brick is 
especially indicated to obtain nonlinear coupled terms between several 
variables. This means in particular that the assembly of the term is performed 
at each call of the assembly of the model and that a Newton algorithm will be 
used to solve the problem. If the term is indeed linear, you should use 
instead::
 
@@ -40,9 +40,9 @@ where ``F`` is a pre-defined constant of the model 
representing the right hand s
 
   getfem::add_linear_term(md, mim, "Grad_u.Grad_Test_u", -1, true, true);
   getfem::add_source_term(md, mim, "F*Test_u");
- 
 
- 
+
+
 
 
 
diff --git a/doc/sphinx/source/userdoc/model_nonlinear_elasticity.rst 
b/doc/sphinx/source/userdoc/model_nonlinear_elasticity.rst
index 01c796b..383b5fa 100644
--- a/doc/sphinx/source/userdoc/model_nonlinear_elasticity.rst
+++ b/doc/sphinx/source/userdoc/model_nonlinear_elasticity.rst
@@ -225,7 +225,7 @@ Incompressible material.
   \intertext{with the additional constraint:}
   i_3( C) = 1
 
-The incompressibility constraint :math:`i_3( C) = 1` is handled with a 
Lagrange multiplier :math:`p` (the pressure) 
+The incompressibility constraint :math:`i_3( C) = 1` is handled with a 
Lagrange multiplier :math:`p` (the pressure)
 
 constraint: :math:`\sigma = -pI \Rightarrow {\hat{\hat{\sigma}}} = 
-p\nabla\Phi\nabla\Phi^{-T}\det\nabla\Phi`
 
@@ -238,7 +238,7 @@ constraint: :math:`\sigma = -pI \Rightarrow 
{\hat{\hat{\sigma}}} = -p\nabla\Phi\
 .. math::
 
   B &= -\int_{\Omega_0} p(\nabla\Phi)^{-T} \det \nabla\Phi : \nabla v  dX \\
-  K &= \int_{\Omega_0} \left( p(\nabla\Phi)^{-T}(\nabla 
h)^{T}(\nabla\Phi)^{-T}\det\nabla\Phi : \nabla v  dX - 
+  K &= \int_{\Omega_0} \left( p(\nabla\Phi)^{-T}(\nabla 
h)^{T}(\nabla\Phi)^{-T}\det\nabla\Phi : \nabla v  dX -
   p(\nabla\Phi)^{-T}(\det \nabla\Phi(\nabla\Phi)^{-T}:\nabla h) : \nabla v 
\right)  dX\\
   &= \int_{\Omega_0} p(\nabla h^T\nabla\Phi^{-T}):(\nabla\Phi^{-1}\nabla 
v)\det\nabla\Phi dX - \int_{\Omega_0} p(\nabla\Phi^{-T}:\nabla 
h)(\nabla\Phi^{-T}:\nabla v)\det\nabla\Phi dX
 
@@ -279,7 +279,7 @@ Since :math:`\frac{\partial}{\partial C} {W}(C) = 
\displaystyle\sum_{j}\frac{\pa
 ``Plane strain hyper-elasticity``
 <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 
-All previous models are valid in volumic domains. Corresponding plane strain 
2D models can be obtained by restricting the stress tensor and the fourth order 
tensor :math:`\mathcal{A}` to their plane components.  
+All previous models are valid in volumic domains. Corresponding plane strain 
2D models can be obtained by restricting the stress tensor and the fourth order 
tensor :math:`\mathcal{A}` to their plane components.
 
 
 
@@ -305,7 +305,7 @@ The Mooney-Rivlin law accepts two optional flags, the first 
one determines if th
 
 The plane strain hyperelastic law takes a pointer on a hyperelastic law as a 
parameter and performs a 2D plane strain approximation.
 
-``md`` is the model variable, ``mim`` the integration method, ``varname`` the 
string being the name of the variable on which the term is added, ``dataname`` 
the string being the name of the data in the model representing the 
coefficients of the law (can be constant or describe on a finite element 
method) and ``region`` is the region on which the term is considered (by 
default, all the mesh). 
+``md`` is the model variable, ``mim`` the integration method, ``varname`` the 
string being the name of the variable on which the term is added, ``dataname`` 
the string being the name of the data in the model representing the 
coefficients of the law (can be constant or describe on a finite element 
method) and ``region`` is the region on which the term is considered (by 
default, all the mesh).
 
 
 The program :file:`nonlinear_elastostatic.cc` in :file:`tests` directory and 
:file:`demo_nonlinear_elasticity.m` in :file:`interface/tests/matlab` directory 
are some examples of use of this brick with or without an incompressibility 
condition.
diff --git a/doc/sphinx/source/userdoc/model_object.rst 
b/doc/sphinx/source/userdoc/model_object.rst
index 97ee537..c79fcce 100644
--- a/doc/sphinx/source/userdoc/model_object.rst
+++ b/doc/sphinx/source/userdoc/model_object.rst
@@ -34,26 +34,26 @@ numbers.
 
    The (tangent) linear system
 
-There are different kinds of variables/data in the model. The variables are 
the 
-unknown of the model. They will be (generally) computed by solving the 
(tangent) 
-linear system built by the model. Generally, the model will have several 
-variables. Each variable has a certain size (number of degrees of freedom) and 
the 
-different variables are sorted in alphanumeric order to form the global 
unknown 
-(:math:`U` in Fig. :ref:`ud-fig-syslin`). Each variable will be associated to 
an 
-interval :math:`I = [n_1, n_2]` which will represent the degrees of freedom 
-indices corresponding to this variable in the global system. The model stores 
also 
-some data (in the same format as the variables). The difference between data 
-and variables is that data is not an unknown of the model. The value of the 
-data should be provided. In some cases (nonlinear models) some variables can 
be 
-considered as some data for certain terms. Variables and data are of two 
kinds. 
-They can have a fixed size, or they can depend on a finite element method (be 
the 
+There are different kinds of variables/data in the model. The variables are the
+unknown of the model. They will be (generally) computed by solving the 
(tangent)
+linear system built by the model. Generally, the model will have several
+variables. Each variable has a certain size (number of degrees of freedom) and 
the
+different variables are sorted in alphanumeric order to form the global unknown
+(:math:`U` in Fig. :ref:`ud-fig-syslin`). Each variable will be associated to 
an
+interval :math:`I = [n_1, n_2]` which will represent the degrees of freedom
+indices corresponding to this variable in the global system. The model stores 
also
+some data (in the same format as the variables). The difference between data
+and variables is that data is not an unknown of the model. The value of the
+data should be provided. In some cases (nonlinear models) some variables can be
+considered as some data for certain terms. Variables and data are of two kinds.
+They can have a fixed size, or they can depend on a finite element method (be 
the
 d.o.f. of a finite element method).
 
-For instance, in the situation described in Fig. :ref:`ud-fig-syslin`, there 
are four variables in the model, namely :math:`X, Y, V` and :math:`W`. The role 
of 
-the model object will be to assemble the linear system, i.e. to fill the sub 
-matrices corresponding to each variable (:math:`R_{X,X}, R_{Y,Y}, R_{V,V}`, 
and 
-:math:`R_{W,W}`) and the coupling terms between two variables (:math:`R_{X,Y}, 
-R_{X,V}, R_{W,V}, \cdots`). This different contributions will be given by the 
+For instance, in the situation described in Fig. :ref:`ud-fig-syslin`, there 
are four variables in the model, namely :math:`X, Y, V` and :math:`W`. The role 
of
+the model object will be to assemble the linear system, i.e. to fill the sub
+matrices corresponding to each variable (:math:`R_{X,X}, R_{Y,Y}, R_{V,V}`, and
+:math:`R_{W,W}`) and the coupling terms between two variables (:math:`R_{X,Y},
+R_{X,V}, R_{W,V}, \cdots`). This different contributions will be given by the
 different bricks added to the model.
 
 The main useful methods on a |mo| object are
@@ -328,10 +328,10 @@ to this condition prescribed with a Lagrange multiplier 
are
 
    \int_{\Gamma} u \mu\ d\Gamma = \int_{\Gamma} u_D \mu\ d\Gamma, \forall \mu 
\in M,
 
-where :math:`M` is an appropriate multiplier space. The contributions to the 
-global linear system can be viewed in Fig. :ref:`ud-fig-syslinDir`. The matrix 
-:math:`B` is the "mass matrix" between the finite element space of the 
variable 
-:math:`u` and the finite element space of the multiplier :math:`\mu`. 
+where :math:`M` is an appropriate multiplier space. The contributions to the
+global linear system can be viewed in Fig. :ref:`ud-fig-syslinDir`. The matrix
+:math:`B` is the "mass matrix" between the finite element space of the variable
+:math:`u` and the finite element space of the multiplier :math:`\mu`.
 :math:`L_{u}` is the right hand side corresponding to the data :math:`u_D`.
 
 .. _ud-fig-syslinDir:
diff --git a/doc/sphinx/source/userdoc/model_plasticity_small_strain.rst 
b/doc/sphinx/source/userdoc/model_plasticity_small_strain.rst
index fb21127..95c5dcc 100644
--- a/doc/sphinx/source/userdoc/model_plasticity_small_strain.rst
+++ b/doc/sphinx/source/userdoc/model_plasticity_small_strain.rst
@@ -131,7 +131,7 @@ A solution would be to solve the whole problem with all the 
unknows, that is :ma
 
    {\mathscr A} : (u_{n+1}, \zeta_{n}, \eta_n) \mapsto \alpha_{n+1}
 
-   
+
 
 with :math:`\eta_n, \zeta_{n}` the right hand side of equations 
:eq:`thetascheme1`, :eq:`thetascheme2`, i.e.
 
@@ -205,7 +205,7 @@ In the decomposition of plastic and elastic part of the 
strain tensor, we assume
 
 and
 
-.. math:: \varepsilon^e_{3,3} + \varepsilon^p_{3,3} = \varepsilon_{3,3} = 0. 
+.. math:: \varepsilon^e_{3,3} + \varepsilon^p_{3,3} = \varepsilon_{3,3} = 0.
 
 The adaptation to the plane strain approximation to plastic model is most of 
the time an  easy task. An isotropic linearized elastic response reads
 
@@ -232,7 +232,7 @@ The plane stress approximation describe generally the 2D 
membrane deformation of
 
 .. math:: \sigma = \left(\hspace{-0.5em}\begin{array}{ccc} \sigma_{1,1} & 
\sigma_{1,2} & 0 \\ \sigma_{1,2} & \sigma_{2,2} & 0 \\ 0 & 0 & 0 
\end{array}\hspace{-0.5em}\right).
 
-We will still denote 
+We will still denote
 
 .. math:: \bar{\sigma} =  \left(\hspace{-0.5em}\begin{array}{cc} \sigma_{1,1} 
& \sigma_{1,2} \\ \sigma_{1,2} & \sigma_{2,2} \end{array}\hspace{-0.5em}\right)
 
@@ -305,13 +305,13 @@ One has
 
 .. math:: \|\mbox{Dev}(\sigma_{n+1})\| = 
2\mu\|\mbox{Dev}(\varepsilon(u_{n+1})) -\varepsilon^p_{n+1}\| = 
\Frac{2\mu}{1+2\mu\theta\Delta t \xi_{n+1}}\|\mbox{Dev}(\varepsilon(u_{n+1})) - 
\zeta_n\|,
 
-Thus, denoting :math:`B = \mbox{Dev}(\varepsilon(u_{n+1})) - \zeta_n`, either 
+Thus, denoting :math:`B = \mbox{Dev}(\varepsilon(u_{n+1})) - \zeta_n`, either
 
 .. math:: 2\mu\|B\| \le \sqrt{\frac{2}{3}}\sigma_y,
 
 and :math:`\xi_{n+1} = 0`, i.e. we are in the elastic case, or  
:math:`\|\mbox{Dev}(\sigma_{n+1})\| =  \sqrt{\frac{2}{3}}` and one obtains
 
-.. math:: 1+2\mu\theta\Delta t \xi_{n+1} = 
\Frac{2\mu\|B\|}{\sqrt{\frac{2}{3}}\sigma_y}, 
+.. math:: 1+2\mu\theta\Delta t \xi_{n+1} = 
\Frac{2\mu\|B\|}{\sqrt{\frac{2}{3}}\sigma_y},
 
 and thus
 
@@ -334,7 +334,7 @@ The plane strain approximation has the same expression 
replacing the 3D strain t
 
 where :math:`\overline{\mbox{Dev}}(\bar{\varepsilon}) = \bar{\varepsilon} - 
\Frac{\mbox{tr}(\bar{\varepsilon})}{3} \bar{I}` is the 2D restriction of the 3D 
deviator.
 
-Moreover, for the yield condition, 
+Moreover, for the yield condition,
 
 .. math:: \|\mbox{Dev}(\sigma)\|^2 = 
4\mu^2\left(\|\overline{\mbox{Dev}}\bar{\varepsilon}(u) - 
\bar{\varepsilon}^p\|^2 + \left(\Frac{\mbox{tr}(\bar{\varepsilon}(u))}{3} 
-\mbox{tr}(\bar{\varepsilon}^p) \right)^2\right).
 
@@ -468,7 +468,7 @@ to be done ...
 
 
 
-Elasto-plasticity bricks 
+Elasto-plasticity bricks
 +++++++++++++++++++++++++
 
 See the test programs :file:`tests/plasticity.cc`, 
:file:`interface/tests/matlab/demo_plasticity.m`, 
:file:`interface/tests/matlab/demo_plasticity.py` and in 
:file:`contrib/test_plasticity`.
@@ -623,7 +623,7 @@ Additionaly, The function: ::
 computes the new stress constraint values supported by the material after a 
load or an unload (once a solve has been done earlier) and upload the variables 
``varname`` and ``datasigma`` as follows:
 
 .. math::
-   
+
    u^{n+1} \Rightarrow u^n \ \ \ \ \ \textrm{ and } \ \ \ \ \ \sigma^{n+1} 
\Rightarrow \sigma^n
 
 Then, :math:`u^n` and :math:`\sigma^n` contains the new values computed and 
one can restart the process.
@@ -642,7 +642,7 @@ The function: ::
       getfem::compute_plastic_part
           (md, mim, mf_pl, varname, previous_varname, ACP, datalambda, datamu, 
datathreshold, datasigma, Plast);
 
-computes on ``mf_pl`` the plastic part of the material, that could appear 
after a load and an unload, into the vector ``Plast``. 
+computes on ``mf_pl`` the plastic part of the material, that could appear 
after a load and an unload, into the vector ``Plast``.
 
 Note that ``datasigma`` should be the vector containing the new stress 
constraint values, i.e. after a load or an unload of the material.
 
diff --git a/doc/sphinx/source/userdoc/model_time_integration.rst 
b/doc/sphinx/source/userdoc/model_time_integration.rst
index 3e22428..80ce7ea 100644
--- a/doc/sphinx/source/userdoc/model_time_integration.rst
+++ b/doc/sphinx/source/userdoc/model_time_integration.rst
@@ -25,7 +25,7 @@ Although time integration scheme can be written directly 
using the model object
 
 * The data `t` represent the time parameter and can be used (either in the 
weak form language or as parameter of some bricks). Before the assembly of the 
system, the data `t` is automatically updated to the time step `n`.
 
-* The problem to be solved at each iteration correspond to the formulation of 
the transient problem in its natural (weak) formulation in which the velocity 
and the acceleration are expressed by the intermediate variables introduced. 
For instance, the translation into the weak form language of the problem 
+* The problem to be solved at each iteration correspond to the formulation of 
the transient problem in its natural (weak) formulation in which the velocity 
and the acceleration are expressed by the intermediate variables introduced. 
For instance, the translation into the weak form language of the problem
 
   .. math::
 
@@ -39,8 +39,8 @@ Although time integration scheme can be written directly 
using the model object
 
 * For all implemented one-step schemes, the time step can be changed from an 
iteration to another for both order one and order two in time problems (or even 
quasi-static problems).
 
-* A scheme for second order in time problem (resp. first order in time) can be 
applied to a second or first order in time or even to a quasi-static problem 
(resp. to a first order or quasi-static problem) without any problem except 
that the initial data corresponding to the velocity/displacement have to be 
initialized with respect ot the order of the scheme. Conversely, of course, a 
scheme for first order problem cannot be applied to a second order in time 
problem. 
- 
+* A scheme for second order in time problem (resp. first order in time) can be 
applied to a second or first order in time or even to a quasi-static problem 
(resp. to a first order or quasi-static problem) without any problem except 
that the initial data corresponding to the velocity/displacement have to be 
initialized with respect ot the order of the scheme. Conversely, of course, a 
scheme for first order problem cannot be applied to a second order in time 
problem.
+
 
 The implicit theta-method for first-order problems
 **************************************************
@@ -75,7 +75,7 @@ When applying this scheme to a variable "u" of the model, the 
following affine d
 
   "Dot_u"
 
-which represent the time derivative of the variable and can be used in some 
brick definition. 
+which represent the time derivative of the variable and can be used in some 
brick definition.
 
 The following data are also added::
 
@@ -153,7 +153,7 @@ When aplying this scheme to a variable "u" of the model, 
the following affine de
 
   "Dot_u", "Dot2_u"
 
-which represent the first and second order time derivative of the variable and 
can be used in some brick definition. 
+which represent the first and second order time derivative of the variable and 
can be used in some brick definition.
 
 The following data are also added::
 
@@ -207,7 +207,7 @@ When aplying this scheme to a variable "u" of the model, 
the following affine de
 
   "Dot_u", "Dot2_u"
 
-which represent the first and second order time derivative of the variable and 
can be used in some brick definition. 
+which represent the first and second order time derivative of the variable and 
can be used in some brick definition.
 
 The following data are also added::
 
@@ -243,7 +243,7 @@ Similarly, every existing model brick of |gf| can be 
applied to "Dot_u". This is
 
    getfem::add_mass_brick(model, mim, "Dot_u");
 
-which adds the same transient term. 
+which adds the same transient term.
 
 VERY IMPORTANT: When adding an existing model brick applied to an affine 
dependent variable such as "Dot_u", it is always assumed that the corresponding 
test function is the one of the corresponding original variable (i.e. "Test_u" 
here). In other words, "Test_Dot_u", the test variable corresponding to the 
velocity, is not used. This corresponds to the choice made to solve the problem 
in term of the original variable, so that the test function corresponds to the 
original variable.
 
@@ -251,7 +251,7 @@ Another example of model brick which can be used to account 
for a Kelvin-Voigt l
 
    getfem::add_isotropic_linearized_elasticity_brick(model, mim, "Dot_u", 
"lambda_viscosity", "mu_viscosity");
 
-when applied to an order two transient elasticity problem. 
+when applied to an order two transient elasticity problem.
 
 Computation on the sequence of time steps
 *****************************************
@@ -260,12 +260,12 @@ Typically, the solve on the different time steps will 
take the following form::
 
 
   for (scalar_type t = 0.; t < T; t += dt) { // time loop
-    
+
     // Eventually compute here some time dependent terms
 
     iter.init();
     getfem::standard_solve(model, iter);
-    
+
     // + Do something with the solution (plot or store it)
 
     model.shift_variables_for_time_integration();
@@ -306,12 +306,12 @@ Assuming that `mf_u` and `mim` are valid finite element 
and integration methods
   // transient part.
   getfem::add_theta_method_for_first_order(model, "u", theta);
   getfem::add_mass_brick(model, mim, "Dot_u");
-  
+
   gmm::iteration iter(residual, 0, 40000);
-  
+
   model.set_time(0.);        // Init time is 0 (not mandatory)
   model.set_time_step(dt);   // Init of the time step.
-  
+
   // Null initial value for the temperature.
   gmm::clear(model.set_real_variable("Previous_u"));
 
@@ -323,10 +323,10 @@ Assuming that `mf_u` and `mim` are valid finite element 
and integration methods
 
   // Iterations in time
   for (scalar_type t = 0.; t < T; t += dt) {
-    
+
     iter.init();
     getfem::standard_solve(model, iter);
-    
+
     // + Do something with the solution (plot or store it)
 
     // Copy the current variables "u" and "Dot_u" into "Previous_u"
diff --git a/doc/sphinx/source/userdoc/parallel.rst 
b/doc/sphinx/source/userdoc/parallel.rst
index a3e1e59..87a7045 100644
--- a/doc/sphinx/source/userdoc/parallel.rst
+++ b/doc/sphinx/source/userdoc/parallel.rst
@@ -16,7 +16,7 @@ using the mesh regions to parallelize assembly procedures.
 
 Nevertheless, the brick system offers a generic parallelization based on MPI
 (communication between processes),
-`METIS <http://glaros.dtc.umn.edu/gkhome/metis/metis/overview>`_ 
+`METIS <http://glaros.dtc.umn.edu/gkhome/metis/metis/overview>`_
 (partition of the mesh)
 and `MUMPS <http://graal.ens-lyon.fr/MUMPS>`_ (parallel sparse direct solver).
 It is available with the compiler option ``-D GETFEM_PARA_LEVEL=2``
@@ -90,7 +90,7 @@ Parallelization of getfem is still considered a "work in 
progress". A certain nu
   assembly procedures of standard bricks use a METIS partition of the
   meshes to distribute the assembly. The tangent/stiffness matrices
   remain distibuted and the standard solve call the parallel version
-  of MUMPS (which accept distributed matrices). 
+  of MUMPS (which accept distributed matrices).
 
   For the moment, the procedure ``actualize_sizes()`` of the model
   object remains sequential and is executed on each process.
@@ -105,7 +105,7 @@ Parallelization of getfem is still considered a "work in 
progress". A certain nu
   * The explicit rhs brick: the given vector is not considered to be
     distributed. Only the given vector on the master process is taken into
     account.
-  
+
   * Constraint brick: The given matrix and rhs are not considered to be
     distributed. Only the given matrix and vector on the master process are
     taken into account.
diff --git a/doc/sphinx/source/userdoc/rmesh.rst 
b/doc/sphinx/source/userdoc/rmesh.rst
index b60edca..4b866fa 100644
--- a/doc/sphinx/source/userdoc/rmesh.rst
+++ b/doc/sphinx/source/userdoc/rmesh.rst
@@ -9,16 +9,16 @@
 Mesh refinement
 ===============
 
-Mesh refinement with the Bank et all method (see [bank1983]_) is available in 
-dimension 1, 2 or 3 for simplex meshes (segments, triangles and tetrahedrons). 
+Mesh refinement with the Bank et all method (see [bank1983]_) is available in
+dimension 1, 2 or 3 for simplex meshes (segments, triangles and tetrahedrons).
 For a given object ``mymesh`` of type |gf_m|, the method::
 
   mymesh.Bank_refine(bv);
 
-refines the elements whose indices are stored in ``bv`` (a |dal_bv| object). 
The 
-conformity of the mesh is kept thanks to additional refinement (the so called 
-green triangles). Information about green triangles (in Figure 
-:ref:`ud-fig-refine`) is stored on the mesh object to gather them for further 
+refines the elements whose indices are stored in ``bv`` (a |dal_bv| object). 
The
+conformity of the mesh is kept thanks to additional refinement (the so called
+green triangles). Information about green triangles (in Figure
+:ref:`ud-fig-refine`) is stored on the mesh object to gather them for further
 refinements (see [bank1983]_).
 
 .. _ud-fig-refine:
diff --git a/doc/sphinx/source/userdoc/xfem.rst 
b/doc/sphinx/source/userdoc/xfem.rst
index 88b077c..1224a04 100644
--- a/doc/sphinx/source/userdoc/xfem.rst
+++ b/doc/sphinx/source/userdoc/xfem.rst
@@ -161,7 +161,7 @@ where "desc" is a string containing the description of the 
boolean operation whi
 
   "c-(a+b)" is the domain of the third level-set minus the union of
   the domains of the two others.
-      
+
   "!a" is the complementary of the domain of a (i.e. it is the
   domain where a(x)>0)
 
diff --git a/doc/sphinx/source/whatsnew/1.7.rst 
b/doc/sphinx/source/whatsnew/1.7.rst
index 111d253..9376c38 100644
--- a/doc/sphinx/source/whatsnew/1.7.rst
+++ b/doc/sphinx/source/whatsnew/1.7.rst
@@ -55,5 +55,5 @@ package). Note that, while it is `documented
 <http://getfem.org/python/index.html>`_ and
 working, the python interface is still considered a *work in
 progress*. You have to enable it explicitly with ``./configure
---enable-python``. An interface to 
+--enable-python``. An interface to
 the Gmm++ sparse matrices and solvers is also provided.
diff --git a/doc/sphinx/source/whatsnew/2.0.rst 
b/doc/sphinx/source/whatsnew/2.0.rst
index 6fb1e9f..786b2b9 100644
--- a/doc/sphinx/source/whatsnew/2.0.rst
+++ b/doc/sphinx/source/whatsnew/2.0.rst
@@ -26,7 +26,7 @@
 
    * Some news features have been introduced in this release:
 
-     * Introduction of level-set objects. Integration methods can be 
+     * Introduction of level-set objects. Integration methods can be
        cut with respect to these level-set and discontinuous
        elements across the level-set are provided.
 
diff --git a/doc/sphinx/source/whatsnew/3.0.rst 
b/doc/sphinx/source/whatsnew/3.0.rst
index e3b5bca..e8fd317 100644
--- a/doc/sphinx/source/whatsnew/3.0.rst
+++ b/doc/sphinx/source/whatsnew/3.0.rst
@@ -13,7 +13,7 @@ Not so many changes, but some of them are incompatible with 
|gf| 2.0:
      directives have to be updated:
 
      * ``#include "gmm_xxx.h"`` should be replaced with
-       ``#include "gmm/gmm_xxx.h"`` 
+       ``#include "gmm/gmm_xxx.h"``
 
      * ``#include "getfem_xxx.h"`` should be replaced with
        ``#include "getfem/getfem_xxx.h"``
diff --git a/doc/sphinx/source/whatsnew/4.1.1.rst 
b/doc/sphinx/source/whatsnew/4.1.1.rst
index ff6d83a..041753f 100644
--- a/doc/sphinx/source/whatsnew/4.1.1.rst
+++ b/doc/sphinx/source/whatsnew/4.1.1.rst
@@ -9,7 +9,7 @@ Minor release, 2010/10/09:
 
    * Scilab files in the archive.
 
-   * Some minor bugs fixed including correction of the Mooney Rivlin 
hyperelastic law. 
+   * Some minor bugs fixed including correction of the Mooney Rivlin 
hyperelastic law.
 
 
 
diff --git a/doc/sphinx/source/whatsnew/4.1.rst 
b/doc/sphinx/source/whatsnew/4.1.rst
index e6d7538..81663e1 100644
--- a/doc/sphinx/source/whatsnew/4.1.rst
+++ b/doc/sphinx/source/whatsnew/4.1.rst
@@ -32,7 +32,7 @@ The main changes are:
      to Luis Saavedra for is important work to re-write the main part of the
      documentations into the sphinx/rst format.
      Documentation for the interfaces are now fully automatic.
-     
+
    * A convection scheme based on Characteristic Galerkin method has been
      added. Useful to update the level-sets.
 
diff --git a/doc/sphinx/source/whatsnew/4.2.rst 
b/doc/sphinx/source/whatsnew/4.2.rst
index 0a4fc47..21a5f12 100644
--- a/doc/sphinx/source/whatsnew/4.2.rst
+++ b/doc/sphinx/source/whatsnew/4.2.rst
@@ -30,7 +30,7 @@ The main changes are:
 
    * The experimental mesher of Getfem is now (partially) interfaced with
      python/scilab/matlab.
-     
+
    * Some tools to verify the consistence (tangent term) of a brick or
      the whole model has been added.
 
diff --git a/doc/sphinx/source/whatsnew/4.3.rst 
b/doc/sphinx/source/whatsnew/4.3.rst
index e891428..9fa0d8b 100644
--- a/doc/sphinx/source/whatsnew/4.3.rst
+++ b/doc/sphinx/source/whatsnew/4.3.rst
@@ -4,7 +4,7 @@
   What's New in |gf| 4.3
 ****************************
 
-This release is a transitional one until next version 5.0. The new high-level 
generic assembly based on a weak form language is working and usable. However, 
the basic model bricks still use the old generic assembly and the new assembly 
is for the moment incompatible with a few things (Nitsche bricks and time 
dispatcher bricks). 
+This release is a transitional one until next version 5.0. The new high-level 
generic assembly based on a weak form language is working and usable. However, 
the basic model bricks still use the old generic assembly and the new assembly 
is for the moment incompatible with a few things (Nitsche bricks and time 
dispatcher bricks).
 
 Released version, 2014/07/14.
 
@@ -17,7 +17,7 @@ The main changes are:
    * The introduction of interpolate transformations in the weak form language
      to deal with the assembly of terms on different meshes or part of meshes.
      Example of applications : mortar methods, periodic boundary conditions,
-     large sliding contact conditions 
+     large sliding contact conditions
 
    * A large sliding contact with friction brick is now working
      (work of Konstantinos Poulios and Yves Renard) and will be extended soon.
diff --git a/doc/sphinx/source/whatsnew/5.0.rst 
b/doc/sphinx/source/whatsnew/5.0.rst
index 31059f7..3ac8eee 100644
--- a/doc/sphinx/source/whatsnew/5.0.rst
+++ b/doc/sphinx/source/whatsnew/5.0.rst
@@ -48,7 +48,7 @@ The main changes are:
    * New im_data object version to store and interpolate data on the
      gauss points of a boundary.
 
-   
+
 
 
 
diff --git a/doc/sphinx/source/whatsnew/5.1.rst 
b/doc/sphinx/source/whatsnew/5.1.rst
index 172768e..2d804e0 100644
--- a/doc/sphinx/source/whatsnew/5.1.rst
+++ b/doc/sphinx/source/whatsnew/5.1.rst
@@ -20,7 +20,7 @@ The main changes are:
    * Large strain plasticity bricks (Simo-Miehe model).
 
    * Redesign of the dof enumeration algorithm: a local sort instead of a 
global one (but it still remains to be parallelized !)
-     
+
    * Addition of a local projection generic function on discontinuous fems.
 
    * Addition of a specific transformation allowing inter-element computation 
in the weak form language (in order to compute inter-element jump of any 
quantity, average value, error estimator, ...)
diff --git a/doc/sphinx/source/whatsnew/5.3.rst 
b/doc/sphinx/source/whatsnew/5.3.rst
index 8e8c069..968f99d 100644
--- a/doc/sphinx/source/whatsnew/5.3.rst
+++ b/doc/sphinx/source/whatsnew/5.3.rst
@@ -16,7 +16,7 @@ The main changes are:
      the spatial gradient and Diff(expression, variable, direction)
      which performs a directional derivative. A few more tensor contraction
      operations has been added on this occasion.
-  
+
    * The support for pyramidal has been extended and stabilized.
 
    * Incomplete pyramidal (13-node) and incomplete prism (15-node)



reply via email to

[Prev in Thread] Current Thread [Next in Thread]