diff --git a/cmake/PugsDoc.cmake b/cmake/PugsDoc.cmake
index 099eda180d7b37d3824fbfac370151c018718952..7be934abae82d06b1e5dcf7265c7710521e02c8e 100644
--- a/cmake/PugsDoc.cmake
+++ b/cmake/PugsDoc.cmake
@@ -78,6 +78,8 @@ if (EMACS)
       DEPENDS
       "${PUGS_SOURCE_DIR}/doc/userdoc.org"
       "${PUGS_SOURCE_DIR}/doc/lisp/userdoc-pdf.el"
+      "${PUGS_SOURCE_DIR}/tools/pgs-pygments.sh"
+      "${PUGS_SOURCE_DIR}/tools/pgs-pygments.py"
       pugs
       pugsdoc-dir
       pugsdoc-download-elpa
diff --git a/doc/lisp/build-doc-config.el b/doc/lisp/build-doc-config.el
index 4cf83f69251b96aaa0f9473a4c6b0ba8fa8a1e51..1da3a5c9fecf3540dfd701380a340722379642cf 100644
--- a/doc/lisp/build-doc-config.el
+++ b/doc/lisp/build-doc-config.el
@@ -43,8 +43,7 @@
 
 (custom-set-variables
  '(org-export-html-postamble nil)
- ;; By now pugs is not known by Pygments raw text is better
- '(org-latex-minted-langs '((pugs "Text") (pugs-error "Text")))
+ '(org-latex-minted-langs '((pugs "Pugs") (pugs-error "Pugs")))
  '(org-confirm-babel-evaluate nil)
  '(org-html-validation-link nil)
  '(org-src-fontify-natively t)
diff --git a/doc/userdoc.org b/doc/userdoc.org
index b550e3315e0a29551624bd5af91345f64cb706c9..6538215304ff3da89bd385e8b7918569490cc85e 100644
--- a/doc/userdoc.org
+++ b/doc/userdoc.org
@@ -43,6 +43,9 @@
 #+LATEX_HEADER_EXTRA: \BeforeBeginEnvironment{warning}{\begin{mdframed}[linecolor=red,backgroundcolor=red!10]}
 #+LATEX_HEADER_EXTRA: \AfterEndEnvironment{warning}{\end{mdframed}}
 
+#+LATEX_HEADER_EXTRA: \usemintedstyle{perldoc}
+#+LATEX_HEADER_EXTRA: \renewcommand{\MintedPygmentize}{${PUGS_SOURCE_DIR}/tools/pgs-pygments.sh}
+
 #+begin_src latex :exports results
   \let\OldTexttt\texttt
   \renewcommand{\texttt}[1]{\fcolorbox{gray!50}{gray!5}{\OldTexttt{#1}}}
@@ -149,7 +152,7 @@ already be discussed.
 - There is no predefined constant in ~pugs~. Here a value is provided
   for ~pi~.
 - There are two kinds of variable in ~pugs~: variables of basic types
-  and variables of high-level types. This two kinds of variable behave
+  and variables of high-level types. These two kinds of variable behave
   almost the same but one must know their differences to understand
   better the underlying mechanisms and choices that we made. See
   sections [[basic-types]] and [[high-level-types]] for details.
@@ -181,13 +184,13 @@ boundary conditions, equations of state, source terms for a specific
 model. Choosing a numerical method or even more, setting the model
 itself, is common in large codes.
 
-In ~pugs~, all these "parameters" are set through a
-DSL[fn:DSL-def]. Thus, when ~pugs~ is launched, it actually executes a
-provided script. A ~C++~ function is associated to each instruction of
-the script. The ~C++~ components of ~pugs~ are completely unaware one of
-the others. ~pugs~ interpreter is responsible for the data flow between
-the components: it manages the data transfer between those ~C++~
-components and ensures that the workflow is properly defined.
+In ~pugs~, all these "parameters" are set through a DSL. Thus, when ~pugs~
+is launched, it actually executes a provided script. A ~C++~ function is
+associated to each instruction of the script. The ~C++~ components of
+~pugs~ are completely unaware one of the others. ~pugs~ interpreter is
+responsible for the data flow between the components: it manages the
+data transfer between those ~C++~ components and ensures that the
+workflow is properly defined.
 
 **** Why?
 
@@ -200,9 +203,9 @@ There are lots of reasons not to use data files. By data file, we
 refer to a set of options that describe physical models, numerical
 methods or their settings.
 
-- Data files are not flexible. This implies in the one hand that
+- Data files are not flexible. This implies on the one hand that
   application scenarios must be known somehow precisely to reflect
-  possible option combinations and in the other hand even defining a
+  possible option combinations and on the other hand even defining a
   specific initial data may require the creation of a new option and
   its associated code (in ~C++~ for instance). \\
   Usually, the last point is addressed by adding a local interpreter
@@ -215,7 +218,7 @@ methods or their settings.
 - Generally data files become rapidly obsolete. An option was not the
   right one, or its type changed to allow other contexts... This puts
   pressure on the user.
-- Even worst, options meaning can depend on other
+- Even worst, option meanings can depend on other
   options. Unfortunately, this happens commonly. For instance, a
   global option can change implicitly the treatment associated to
   another one. This is dangerous since writing or reading the data
@@ -237,10 +240,10 @@ files or scripts), but it presents several drawbacks.
   things that should not be changed.
 - Again, one can easily have access to irrelevant options and it
   requires a great knowledge of the code to find important ones.
-- With that regard, defining a simulation properly can be a difficult
+- With this in mind, defining a simulation properly can be a difficult
   task. For instance, in the early developments of ~pugs~ (when it was
   just a raw ~C++~ code) it was tricky to change boundary conditions for
-  coupled physics.
+  multiphysics problems.
 - Another difficulty is related to the fact that code's internal API
   is likely to change permanently in a research code. Thus valid
   constructions or settings may become rapidly obsolete.  In other
@@ -259,7 +262,7 @@ solution to all problems. However, it offers some advantages.
 - It allows to structure the code in the sense that new developments
   have to be designed not only focusing on the result but also on the
   way it should be used (its interactions with the scripting language).
-- In the same idea, it provides a framework that drives the desired
+- In the same vein, it provides a framework that drives the desired
   principle of "do simple things and do them well".
 - There are no hidden dependencies between numerical options: the DSL
   code is easier to read (than data files) and is less likely to
@@ -292,10 +295,10 @@ developer, by reading it, the user will have a better understanding of
 the development choices and the underlying policy that the code
 follows.
 
-Actually the development framework imposed by the DSL tends to guide
-writing of new methods.
+Actually the development framework imposed by the DSL is a guideline
+for writing of new methods.
 
-- In the process of writing a *new numerical methods*, one must create
+- In the process of writing a *new numerical method*, one must create
   *new functions in the language*. Moreover, if a new method is close to
   an existing one, it is generally *better* to use completely new
   underlying ~C++~ code than to patch existing methods. Starting from a
@@ -307,7 +310,7 @@ writing of new methods.
     numerical method and the *design* of the code.
   - From the computer science point of view, early design for new
     numerical methods is generally wrong: usually one cannot
-    anticipate precisely enough eventual problems or method
+    anticipate precisely enough possible problems or method
     corrections.
   - It is much more difficult to introduce bugs in existing methods,
     since previously validated methods are unchanged!
@@ -324,7 +327,7 @@ writing of new methods.
     deteriorated. At this time, it is likely that the numerical method
     design is finished, thus (re)designing the source code makes more
     sense.
-- Another consequence is that utilities are not be developed again and
+- Another consequence is that utilities are not developed again and
   again.
   - This implies an /a priori/ safer code: utilities are well tested and
     validated.
@@ -332,12 +335,12 @@ writing of new methods.
   - The code of numerical methods is not polluted by environment
     instructions (data initialization, error calculation,
     post-processing,...)
-  - The counterpart is somehow classical. In the one hand, the
+  - The counterpart is somehow classical. On the one hand, the
     knowledge of existing utilities is required, this document tries
-    to address a part of it. In the other hand, if the developer
+    to address a part of it. On the other hand, if the developer
     requires a new utility, a good practice is to discuss with the
-    other ones to check if it could benefit to them. Then one can
-    determine if it should integrate rapidly or not the main
+    other developers to check if it could benefit to them. Then one
+    can determine if it should integrate rapidly or not the main
     development branch.
 
 ***** Why not python or any other scripting language?
@@ -348,7 +351,7 @@ too much freedom: it is not easy to protect data. For instance in the
 [[high-level-types]]). It is important since it prevents the user from
 modifying data in inconsistent ways. Also, one must keep in mind that
 constraining the expressiveness is actually a strength. As said
-before, one can warranty coherence of the data, perform calculations
+before, one can warranty consistency of the data, perform calculations
 without paying attention to the parallelism aspects,... Observe that
 it is not a limitation: if the DSL's field of application needs to be
 extended, it is always possible. But these extensions should never
@@ -363,7 +366,7 @@ Finally, python is ugly.
 *** A high-level language
 
 Following the previous discussion, the reader should now understand
-the motivations that drove the design choices that conduct to build
+the motivations that drove the design choices that conducted to build
 ~pugs~ as a ~C++~ toolbox driven by a user friendly language.
 
 #+begin_verse
@@ -409,7 +412,7 @@ define high-level optimizations.
   ~pugs~ script.
 
 Another benefit of not providing low-level instructions is that the
-scripts are more easy to write and read, and it is more difficult to
+scripts are easier to write and read, and it is more difficult to
 write errors.
 
 * Language
@@ -430,8 +433,8 @@ linefeed string (there is no character type in ~pugs~, just strings).
 
 Actually, ~cout~ is itself a variable, we will come to this later.
 
-~pugs~ is a strongly typed language. It means that a variable *cannot*
-change of type in its lifetime.
+~pugs~ is a strongly typed language. It means that the type of a
+variable *cannot* change in its lifetime.
 
 
 *** Declaration and affectation syntax
@@ -443,7 +446,7 @@ To declare a variable ~v~ of a given type ~V~, one simply writes
   let v:V;
 #+END_SRC
 
-This instruction is read as
+This instruction reads as
 #+begin_verse
 Let $v\in V$.
 #+end_verse
@@ -589,7 +592,7 @@ we give a few examples.
 #+END_SRC
 #+results: out-of-scope-variable-use
 
-**** Variable name *cannot* be reused in an enclosed scope
+**** A variable name *cannot* be reused in an enclosed scope
 #+NAME: nested-scope-variable-example
 #+BEGIN_SRC pugs-error :exports both :results output
   {
@@ -613,9 +616,9 @@ read.
 
 *** Basic types<<basic-types>>
 
-Basic types in ~pugs~ are boolean (~B~), natural integers (~N~), integers (~Z~),
-real (~R~), small vectors (~R^1~, ~R^2~ and ~R^3~), small matrices (~R^1x1~, ~R^2x2~
-and ~R^3x3~) and strings (~string~).
+Basic types in ~pugs~ are boolean (~B~), natural integers (~N~), integers
+(~Z~), real numbers (~R~), small vectors (~R^1~, ~R^2~ and ~R^3~), small
+matrices (~R^1x1~, ~R^2x2~ and ~R^3x3~) and strings (~string~).
 
 #+BEGIN_note
 Observe that while mathematically, obviously $\mathbb{R} = \mathbb{R}^1
@@ -624,14 +627,15 @@ in ~pugs~ and are *not implicitly* convertible from one to the other!
 
 This may sound strange but there are few reasons for that.
 - First, these are the reflect of internal ~pugs~ ~C++~-data types that
-  are used to write algorithms. In its core design pugs aim at writing
-  numerical methods generically with regard to the dimension. One of
-  the ingredients to achieve this purpose is to use dimension $1$
-  vectors and matrices when some algorithms reduce to dimension $1$
-  instead of ~double~ values. To avoid ambiguity that may arise in some
-  situations (this can lead to very tricky code), we decided to forbid
-  automatic conversions of these types with ~double~. When designing the
-  language, we adopted the same rule to avoid ambiguity.
+  are used to write algorithms. In its core design pugs aims at
+  writing numerical methods generically with regard to the
+  dimension. One of the ingredients to achieve this purpose is to use
+  dimension $1$ vectors and matrices when some algorithms reduce to
+  dimension $1$ instead of ~double~ values. To avoid ambiguity that may
+  arise in some situations (this can lead to very tricky code), we
+  decided to forbid automatic conversions of these types with
+  ~double~. When designing the language, we adopted the same rule to
+  avoid ambiguity.
 - A second reason is connected to the first one. Since ~pugs~ aims at
   providing numerical methods for problems in dimension $1$, $2$ or
   $3$, this allows to distinguish the nature of the underlying objects.
@@ -639,8 +643,8 @@ This may sound strange but there are few reasons for that.
     defining a mesh in dimension $d$ are elements of $\mathbb{R}^d$,
   - or that a velocity or a displacement are also defined as
     $\mathbb{R}^d$ values.
-  Thus using ~R^1~ in dimension $1$ for this kind data precise their
-  nature in some sense .
+  Thus using ~R^1~ in dimension $1$ for this kind of data makes precise
+  their nature in some sense .
 #+END_note
 
 **** Expression types
@@ -735,9 +739,9 @@ which is not a surprise. However, the use of the ~+=~ operator results
 in the modification of the stored value. There is no copy.
 
 Actually, this is not really important from the user point of
-view. One just have to keep in mind that, as it will be depicted
-after, high-level variables *are not mutable*: their values can be
-*replaced* by a new ones but *cannot be modified*.
+view. One just has to keep in mind that, as it will be depicted below,
+high-level variables *are not mutable*: their values can be *replaced* by
+new ones but *cannot be modified*.
 
 *** Implicit type conversions<<implicit-conversion>>
 
@@ -868,11 +872,11 @@ are sorted by type of left hand side variable.
   | ~R^1~                           |
 
 - ~R^2~: vector of dimension 2 ($\mathbb{R}^2$) left hand side variable.
-  | ~R^2 =~ allowed expression type               |
-  |---------------------------------------------|
-  | ~R^2~                                         |
-  | ~0~  (special value)                          |
-  | list of 2 scalar (~B~, ~N~, ~Z~ or ~R~) expressions |
+  | ~R^2 =~ allowed expression type                |
+  |----------------------------------------------|
+  | ~R^2~                                          |
+  | ~0~  (special value)                           |
+  | list of 2 scalars (~B~, ~N~, ~Z~ or ~R~) expressions |
   An example of initialization using an $\mathbb{R}^2$ value or the special value ~0~ is
   #+NAME: affectation-to-R2-from-list
   #+BEGIN_SRC pugs :exports both :results output
@@ -971,6 +975,27 @@ are sorted by type of left hand side variable.
   | ~R^2x2~                            |
   | ~R^3x3~                            |
   | ~string~                           |
+  The stored value is the same as the output value described
+  above. Here is an example.
+  #+NAME: affectation-to-string-example
+  #+BEGIN_SRC pugs :exports both :results output
+    let s_from_B:string,
+        s_from_B = 2>1;
+
+    let s_from_R2:string,
+        s_from_R2 = [ -3.5, 1.3];
+
+    let s_from_R3x3:string,
+        s_from_R3x3 = [[ -3, 2.5, 1E-2],
+                       [  2, 1.7,   -2],
+                       [1.2,   4,  2.3]];
+
+    cout << "s_from_B    = " << s_from_B << "\n";
+    cout << "s_from_R2   = " << s_from_R2 << "\n";
+    cout << "s_from_R3x3 = " << s_from_R3x3 << "\n";
+  #+END_SRC
+  the output is
+  #+RESULTS: affectation-to-string-example
 
 ***** List of defined operator ~+=~ for basic types.
 
@@ -1044,6 +1069,23 @@ are sorted by type of left hand side variable.
   | ~R^2x2~                             |
   | ~R^3x3~                             |
   | ~string~                            |
+  The concatenated value is the same as the output value described
+  above, for instance
+  #+NAME: concatenate-to-string-example
+  #+BEGIN_SRC pugs :exports both :results output
+    let s:string, s = "foo";
+
+    s += [ -3.5, 1.3, -2];
+    s += "_";
+    s += 1>2;
+    s += "_";
+    s += [[ -3, 2.5],
+          [  2, 1.7]];
+
+    cout << "s = " << s << "\n";
+  #+END_SRC
+  the output is
+  #+RESULTS: concatenate-to-string-example
 
 ***** List of defined operator ~-=~ for basic types.
 
@@ -1181,8 +1223,8 @@ are sorted by type of left hand side variable.
 Observe that for these small matrix types ($\mathbb{R}^{d\times d}$) the
 construction ~A *= B;~ where ~B~ is a matrix of the same type as ~A~ is not
 allowed. The main reason for that is that for $d>1$ this operation has
-no interests since it requires a temporary. One will see bellow that
-it is possible to write ~A = A*B;~ if needed.
+no interest since it requires a temporary. One will see bellow that it
+is possible to write ~A = A*B;~ if needed.
 #+END_note
 
 - ~string~: the ~*=~ operator is not defined for left hand side string variables.
@@ -1234,7 +1276,7 @@ The ~not~, ~+~ and ~-~ operators apply to the *expression* on their right. ~++~
 and ~--~ operators apply only to a *variable* that can be positioned
 before (pre increment/decrement) or after the token (post
 increment/decrement). These operators are also inspired from their ~C++~
-counterparts for commodity.
+counterparts for convenience.
 
 The ~+~ unary operator is a convenient operator that is *elided* when
 parsing the script.
@@ -1287,30 +1329,30 @@ This code produces
 **** Binary operators
 
 Syntax for binary operators follows again a classical structure: if
-~exp1~ and ~exp2~ denotes two expressions and if ~op~ denotes a binary
+~exp1~ and ~exp2~ denote two expressions and if ~op~ denotes a binary
 operator, one simply writes ~exp1 op exp2~.
 
 Here is the list of binary operators
-| keyword | operator              |
-|---------+-----------------------|
-| ~and~     | logic and             |
-| ~or~      | logic or              |
-| ~xor~     | logic exclusive or    |
-|---------+-----------------------|
-| ~==~      | equality              |
-| ~!=~      | non-equality          |
-| ~<~       | lower than            |
-| ~<=~      | lower or equal than   |
-| ~>~       | greater than          |
-| ~>=~      | greater or equal than |
-|---------+-----------------------|
-| ~<<~      | shift left            |
-| ~>>~      | shift right           |
-|---------+-----------------------|
-| ~+~       | sum                   |
-| ~-~       | difference            |
-| ~*~       | product               |
-| ~/~       | division              |
+| keyword | operator                 |
+|---------+--------------------------|
+| ~and~     | logic and                |
+| ~or~      | logic or                 |
+| ~xor~     | logic exclusive or       |
+|---------+--------------------------|
+| ~==~      | equality                 |
+| ~!=~      | non-equality             |
+| ~<~       | lower than               |
+| ~<=~      | lower than or equal to   |
+| ~>~       | greater than             |
+| ~>=~      | greater than or equal to |
+|---------+--------------------------|
+| ~<<~      | shift left               |
+| ~>>~      | shift right              |
+|---------+--------------------------|
+| ~+~       | sum                      |
+| ~-~       | difference               |
+| ~*~       | product                  |
+| ~/~       | division                 |
 
 Binary operators can be defined for high-level types. For basic types,
 they follow a few rules.
@@ -1329,7 +1371,7 @@ they follow a few rules.
     \end{equation*}
   #+end_src
 - Comparison operators ~==~, ~!=~, ~<~, ~<=~, ~>~ and ~>=~ are defined for all
-  basic scalar type and return a boolean value.
+  basic scalar types and return a boolean value.
   #+begin_src latex :results drawer :exports results
     \begin{equation*}
       \forall \mathbb{S}_1, \mathbb{S}_2 \in \{\mathbb{B},\mathbb{N},\mathbb{Z},\mathbb{R}\},
@@ -1386,7 +1428,7 @@ they follow a few rules.
   combinations of basic types. We classify them by their returned
   types.
 
-  - Their is no arithmetic operation that returns a boolean ~B~.
+  - There are no arithmetic operations that return a boolean ~B~.
 
   - Operators that return a natural integer ~N~.
     #+begin_src latex :results drawer :exports results
@@ -1528,7 +1570,7 @@ they follow a few rules.
 To avoid confusions, the operators precedence in ~pugs~ language follows
 the same rules as in ~C++~.
 
-This is summarized in the following table, where ~a~ and ~b~ denotes two
+This is summarized in the following table, where ~a~ and ~b~ denote two
 expressions.
 | Precedence | Operator |
 |------------+----------|
@@ -1591,7 +1633,7 @@ the output is
 *** High-level types<<high-level-types>>
 
 Aside from the basic types described in the previous section, ~pugs~
-also deals with "high-level" types. This term is more to understand as
+also deals with "high-level" types. This term is better understood as
 "non-basic types". The ~pugs~ language is not object oriented to keep it
 simple.
 
@@ -1610,7 +1652,7 @@ script).
 The second difference is that data of these types are *constant*. More
 precisely, the content of high-level variables can be replaced by a
 new one *but* it cannot be modified. For this reason, the following
-operators can never be applied to variables of this kind
+operators can never be applied to variables of these kinds
 | forbidden operators | description                 |
 |---------------------+-----------------------------|
 | ~++~                  | pre/post increment operator |
@@ -1620,7 +1662,7 @@ operators can never be applied to variables of this kind
 | ~*=~                  | assignment by product       |
 | ~/=~                  | assignment by quotient      |
 
-We conclude by stating that if access operator ~[]~ can eventually be
+We conclude by stating that if access operator ~[]~ can possibly be
 defined for high-level types, it should be done with care. It is not
 recommended.
 
@@ -1650,10 +1692,10 @@ In this example, we are dealing with 3 ~mesh~ variables.
   duplication of the mesh in memory.
 - Then, the variable ~m3~ is defined to refer to a *new* mesh (called
   $\mathcal{M}_b$. It is /similar/ to $\mathcal{M}_a$, but it is *not* the
-  same one! When ~m3~ is defined two meshes $\mathcal{M}_a$ and
+  same one! When ~m3~ is defined, two meshes $\mathcal{M}_a$ and
   $\mathcal{M}_b$ (and two distinct connectivities) are resident in
   memory.
-- Finally, the last instruction (~m3 = m1;~) sets ~m3~ to also refer
+- Finally, the last instruction (~m3 = m1;~) sets ~m3~ to also refer to
   $\mathcal{M}_a$. Since no other variable refers to $\mathcal{M}_b$,
   this mesh is destroyed (memory is freed). At the end of the program,
   all the variables ~m1~, ~m2~ and ~m3~ are referring to $\mathcal{M}_a$
@@ -1721,8 +1763,8 @@ It produces the following error
 While the variable ~x~ is defined *before* ~y~, this kind of construction is
 forbidden. From a technical point of view, this behavior would be easy
 to change (allow to use the fresh value of ~x~ in the definition of ~y~),
-but this make the code unclear and this is not the purpose of compound
-types.
+but this makes the code unclear and this is not the purpose of
+compound types.
 
 #+BEGIN_note
 Observe that there is no implicit conversion when dealing with
@@ -1755,9 +1797,9 @@ is the operator ~=~.
 **** Use of compound types
 
 Actually using compound types the way it is presented in this
-paragraph is not recommend. The purpose of compounds types and
-compound affectations is related to functions. As one will see in
-section [[functions]], functions can return compound values, thus compound
+paragraph is not recommend. The purpose of compound types and compound
+affectations is related to functions. As one will see in section
+[[functions]], functions can return compound values, thus compound
 affectations (or definitions) are needed to get returned values in
 that case.
 
@@ -1801,7 +1843,7 @@ Observe that the special value ~0~ is used there.
 
 **** Tuple purpose
 
-Tuples variables are just lists of data of the same type. In the ~pugs~
+Tuple variables are just lists of data of the same type. In the ~pugs~
 language, one cannot access to a specific value of the list nor alter
 one of them. This is not something that should ever change. Tuples are
 not arrays! The ~pugs~ language is not meant to allow low-level
@@ -1816,10 +1858,10 @@ boundary conditions to a method.
 The ~pugs~ language supports classical statements to control the data
 flow. For simplicity, these statements syntax follow their ~C++~
 counterpart. The only statement that is not implemented in ~pugs~ is the
-~switch...case~. This may change but in the one hand, up to now it was
-never necessary (up to now, we did not encountered the need to chain
-~if...else~ statements), and on the other hand, keeping the language as
-simple as possible remains the policy in ~pugs~ development.
+~switch...case~. This may change but on the one hand, up to now it has
+never been necessary (up to now, we did not encountered the need to
+chain ~if...else~ statements), and on the other hand, keeping the
+language as simple as possible remains the policy in ~pugs~ development.
 
 *** ~if...else~ statement
 
@@ -1902,11 +1944,11 @@ are optional. The ~condition~ argument, if it is present, *must* be a
 boolean value (type ~B~). If it is absent, the default value ~true~ is
 used.
 
-- The ~declarationinstruction~ is execute only *once* /before/ the beginning
+- The ~declarationinstruction~ is executed only *once* /before/ the beginning
   of the loop. The lifetime of the variable declared here is defined
   by the ~for~ instruction itself.
 
-- The ~condition~ is evaluate /before/ each loop.
+- The ~condition~ is evaluated /before/ each loop.
 
 - The ~postinstruction~ is executed /after/ each loop.
 
@@ -1986,7 +2028,7 @@ The ~statement~ is either a single instruction or a block of
 instructions. The ~condition~ is an expression of boolean value (type
 ~B~).
 
-This time is the ~condition~ is not satisfied (~false~ when reaching the
+This time if the ~condition~ is not satisfied (~false~ when reaching the
 ~while~ instruction), the ~statment~ is never executed.
 
 An example of the ~while~ loop is the following.
@@ -2087,10 +2129,10 @@ act as operators.
 
 #+BEGIN_note
 Actually these functions are not strictly /pure functions/ in the
-computer science sense. The reason is that they can eventually have
-side effects. As an example, it is possible to modify the random seed
-used by the code. In that case, the modified value is not a variable
-of the language itself but the internal random seed.
+computer science sense. The reason is that they may have side
+effects. As an example, it is possible to modify the random seed used
+by the code. In that case, the modified value is not a variable of the
+language itself but the internal random seed.
 #+END_note
 
 *** Implicit type conversion for parameters and returned values
@@ -2167,10 +2209,10 @@ Using compound types as input and output, one can write
 This meaningless example produces the following result.
 #+results: R22-R-string-to-R-string-function
 
-**** Lifetime of functions arguments
+**** Lifetime of function arguments
 
-The arguments used to define a function are *local* variables that
-exists only during the evaluation of the function.
+The arguments used to define a function are *local* variables that exist
+only during the evaluation of the function.
 
 Let us consider the following example
 #+NAME: lifetime-of-function-args
@@ -2187,7 +2229,7 @@ This gives the expected result: the value of the variable ~a~ is
 unchanged.
 #+results: lifetime-of-function-args
 
-**** Non-arguments variables in function expressions
+**** Non-argument variables in function expressions
 
 Here we discuss rapidly of using variables (which are not arguments)
 in function expressions.
@@ -2209,8 +2251,8 @@ in function expressions.
 Running the example, one gets
 #+results: non-arg-variables-in-functions
 While the function itself is a constant object, one sees that since
-the value of ~a~ is changed, the value function is implicitly
-modified. /This is a dangerous feature and should be avoid!/
+the value of ~a~ is changed, the function value is implicitly
+modified. /This is a dangerous feature and should be avoided!/
 
 Since functions themselves are variables one can use functions in
 function expressions.
@@ -2231,7 +2273,7 @@ output since ~cout~ does not handle compound types output. One gets
 **** Lifetime of user-defined functions
 
 Since functions are somehow variables, the lifetime of functions
-follows the similar rules.
+follows similar rules.
 
 Let us give an example
 #+NAME: functions-lifetime
@@ -2282,19 +2324,19 @@ produces the following compilation time error
 *** Builtin functions<<builtin-functions>>
 
 In ~pugs~ language, builtin functions are ~C++~ pieces of code that can be
-called in scripts. There usage is very similar to user-defined
+called in scripts. Their usage is very similar to user-defined
 functions. They differ from user-defined functions in three points.
-- Builtin functions can have no parameter or no returned value.
+- Builtin functions may have no parameter or no returned value.
 - Builtin functions are polymorphic. More precisely, this means that
   the signature of a builtin function is also defined by its expected
-  arguments types.
+  argument types.
 - Builtin functions can take user-defined functions as parameters.
   - user-defined functions cannot take functions as parameters
   - builtin functions cannot take builtin functions as parameters
     (actually, this is not a limitation since it is trivial to embed a
     builtin function into a user-defined one).
 
-Here is a simple example of builtin function embedding in a user
+Here is a simple example of builtin function embedded in a user
 function.
 #+NAME: builtin-function-embedding
 #+BEGIN_SRC pugs :exports both :results output
@@ -2323,7 +2365,7 @@ mathematical functions, one writes in the preamble of the script
 #+BEGIN_warning
 A work in progress
 - At the time of writing this documentation, one should note that
-  module inter-dependencies is still not implemented.
+  module inter-dependencies are still not implemented.
 - Also, (and especially with regard to the ~scheme~ module), module
   contents are likely to change and to be reorganized.
 - Finally it is almost sure that modules will be equipped with a
@@ -2332,22 +2374,21 @@ A work in progress
   be more natural.
 #+END_warning
 
-One can access to the list of available modules inside the language.
+One can access the list of available modules inside the language.
 #+NAME: get-available-modules
 #+BEGIN_SRC pugs :exports both :results output
   cout << getAvailableModules() << "\n";
 #+END_SRC
 The output lists all available modules
 #+RESULTS: get-available-modules
-Let us comment a bit this output. One notices that there are two kind
+Let us comment a bit this output. One notices that there are two kinds
 of modules. Modules that are automatically imported (tagged with a ~*~)
 and the other ones.
 
 In this section we will not describe exhaustively the whole module
-contents but will give the basic information that should allow the
-user to find his way. To do so, it is important to examine carefully
-the content of the ~core~ module, since it contains some helper
-functions.
+content but will give the basic information that should allow the user
+to find his way. To do so, it is important to examine carefully the
+content of the ~core~ module, since it contains some helper functions.
 
 *** The ~core~ module
 
@@ -2387,7 +2428,7 @@ operator ~<<~ is an ~ostream~, the result of the operation is also an
 One can overload the ~ostream <<~ construction for high-level types.
 
 Other variables of the type ~ostream~ can be created (in order to write
-to files for instance) as one should see below.
+to files for instance) as we will see below.
 
 **** ~core~ provided functions
 
@@ -2423,8 +2464,7 @@ name of the function and its input and output sets.
 
 #+BEGIN_warning
 Observe that this function does not provide the list of operators that
-are defined in the module (eventually associated to the defined
-types).
+are defined in the module (possibly associated to the defined types).
 #+END_warning
 
 ***** ~getPugsBuildInfo: void -> string~
@@ -2475,8 +2515,8 @@ Running this example gives
 
 ***** ~ofstream: string -> ostream~
 
-This function is used to create an ~ostream~ that actually write to the
-file which name is given by the ~string~ argument. One should notice
+This function is used to create an ~ostream~ that actually writes to the
+file whose name is given by the ~string~ argument. One should notice
 that the file is *created* at the function call. If a file with the same
 name existed, it is *erased*.
 #+NAME: ofstream-example
@@ -2484,10 +2524,8 @@ name existed, it is *erased*.
   let fout:ostream, fout = ofstream("filename.txt");
   fout << [1,2] << " is a vector of R^2\n";
 #+END_SRC
-Running this example produces no output
-#+RESULTS: ofstream-example
-But a file is created (in the execution directory), with the name
-~"filename.txt"~. Its content is
+Running this example produces no output but a file is created (in the
+execution directory), with the name ~"filename.txt"~. Its content is
 #+NAME: cat-filename-txt
 #+BEGIN_SRC shell :exports results :results output
   cat filename.txt
@@ -2558,15 +2596,15 @@ their ~C++~ man pages for details.
 #+END_note
 
 #+BEGIN_note
-Let us comment the use of the ~pow~ function. Actually one can wonder
+Let us comment on the use of the ~pow~ function. Actually one can wonder
 why we did not use a syntax like ~x^y~? The reason is that if
-mathematically ${x^y}^z = x^{(y^z)}$, many software treat it (by mistake)
+mathematically ${x^y}^z = x^{(y^z)}$, many softwares treat it (by mistake)
 as ${(x^y)}^z$. Thus, using the ~pow~ function avoids any confusion.
 #+END_note
 
 *** The ~mesh~ module
 
-This is an important module. It provides mesh utilities tools.
+This is an important module. It provides mesh utility tools.
 #+NAME: get-module-info-mesh
 #+BEGIN_SRC pugs :exports both :results output
   cout << getModuleInfo("mesh") << "\n";
@@ -2582,23 +2620,23 @@ a ~mesh~ that is either designated by an integer or by a ~string~.
 
 A ~boundary~ can designate a set of nodes, edges or faces. The ~boundary~
 (descriptor) itself is not related to any ~mesh~, thus the nature of the
-~boundary~ is precised when it is used with a particular ~mesh~.
+~boundary~ is made precise when it is used with a particular ~mesh~.
 
 #+BEGIN_warning
 A ~boundary~ *cannot* be used to refer to an interface (/ie/ an inner set of
-items)..
+items).
 #+END_warning
 
 ***** ~zone~
 
 Following the same idea, a ~zone~ is a descriptor of a set of cells. It
 can be either defined by an integer or by a ~string~. Its meaning is
-precised when it is associated with a ~mesh~.
+made precise when it is associated with a ~mesh~.
 
 ***** ~mesh~
 
 The type ~mesh~ is an *abstract* type that is used to store meshes. A
-variable of that type can refer for instance unstructured meshes of
+variable of that type can refer for instance to unstructured meshes of
 dimension 1, 2 or 3.
 
 The following binary operator is provided.
@@ -2634,7 +2672,7 @@ This type is used to designate kinds of items (cell, face, edge or node).
 
 The type ~item_value~ is an abstract type use to designate arrays of
 values defined on the entities of a ~mesh~. Actually, these values are
-not associated to the mesh itself but on the *connectivity*. The values
+not associated to the mesh itself but to the *connectivity*. The values
 on the entities can be of type ~B~, ~N~, ~Z~, ~R~, ~R^1~, ~R^2~, ~R^3~, ~R^1x1~, ~R^2x2~
 and ~R^3x3~. Entities themselves can be cells, faces, edges or nodes.
 
@@ -2645,17 +2683,17 @@ By now, no mathematical operation is defined on ~item_value~ variables.
 #+END_warning
 
 Moreover, ~item_value~ variables can be post processed. Observe that
-edges or faces values cannot be post processed since neither ~VTK~ nor
+edge or face values cannot be post processed since neither ~VTK~ nor
 ~Gnuplot~ can handle these data.
 
 ***** ~sub_item_value~
 
-This abstract type is handle values defined on the sub items of the
+This abstract type handles values defined on the sub items of the
 items of a ~mesh~. Similarly these values are attached to a *connectivity*
-and not to a mesh. Values can bed of type The values on the entities
-can be of type ~B~, ~N~, ~Z~, ~R~, ~R^1~, ~R^2~, ~R^3~, ~R^1x1~, ~R^2x2~ or ~R^3x3~. An
-example of ~sub_item_value~ is the $\mathbf{C}_{jr}$ vectors which are
-defined at each node of each cell.
+and not to a mesh. The values associated to the sub items can be of
+type ~B~, ~N~, ~Z~, ~R~, ~R^1~, ~R^2~, ~R^3~, ~R^1x1~, ~R^2x2~ or ~R^3x3~. An example of
+~sub_item_value~ is the $\mathbf{C}_{jr}$ vectors which are defined at each
+node of each cell.
 
 These variables are used to pass data from one function to
 another. They cannot be post processed.
@@ -2695,7 +2733,7 @@ Creates a zone descriptor from a ~string~ name
 
 ****** ~zoneTag: Z -> zone~
 
-Associates a zone descriptor from zone tag
+Associates a zone descriptor from a zone tag
 #+BEGIN_SRC pugs :exports both :results none
   import mesh;
 
@@ -2707,9 +2745,9 @@ Associates a zone descriptor from zone tag
 Creates a cartesian mesh of dimension $d\in\{1, 2, 3\}$. The produced
 cartesian grid is aligned with the axis and made of identical cells.
 
-The two first arguments are two opposite corners of the box (or
+The first two arguments are two opposite corners of the box (or of the
 segment in 1d) and the list of natural integers (type ~(N)~) sets the
-number of *cells* in each direction. Thus size of the list of ~N~ is $d$.
+number of *cells* in each direction. Thus the size of the list ~N~ is $d$.
 
 For instance one can write:
 #+BEGIN_SRC pugs :exports both :results none
@@ -2866,8 +2904,8 @@ The ~mesh~ is represented in Figure [[fig:gmsh-hybrid-2d]].
 #+RESULTS: diamond-dual-img
 
 #+BEGIN_note
-The mesh storage mechanisms in ~pugs~ is such that the diamond dual mesh
-is built only once. This means that is one writes for instance
+The mesh storage mechanism in ~pugs~ is such that the diamond dual mesh
+is built only once. This means that if one writes for instance
 #+BEGIN_SRC pugs :exports both :results none
   import mesh;
 
@@ -2919,7 +2957,7 @@ The ~mesh~ is represented in Figure [[fig:gmsh-hybrid-2d]].
 #+RESULTS: median-dual-img
 
 #+BEGIN_note
-In ~pugs~, the storage mechanisms of median dual meshes follows the same
+In ~pugs~, the storage mechanisms of median dual meshes follow the same
 rules as the diamond dual meshes. As long as the primary mesh lives
 and as long as the median dual mesh is referred, it is kept in memory,
 thus constructed only once.
@@ -2931,7 +2969,9 @@ available in parallel
 #+END_warning
 
 ***** Item types
-The following functions are used to refer to ~item_type~
+
+\\
+The following functions are used to designate a specific ~item_type~
 - ~cell: void -> item_type~
 - ~face: void -> item_type~
 - ~edge: void -> item_type~
@@ -2964,8 +3004,8 @@ In this example, we set three arrays defined at all nodes, all the
 This function allows to compute a new mesh as the transformation of a
 given mesh by displacing its nodes through a user defined function.
 
-For a mesh of dimension $d$, the mesh must be a function $\mathbb{R}^d
-\to\mathbb{R}^d$.
+For a mesh of dimension $d$, the transformation must be a function
+$\mathbb{R}^d \to\mathbb{R}^d$.
 
 #+BEGIN_SRC pugs :exports both :results none
 import mesh;
@@ -2983,7 +3023,7 @@ write_mesh(gnuplot_writer("transformed"), m1);
 
 #+BEGIN_note
 One should keep in mind that the mesh produced by the ~transform~
-function *shares* the same connectivity than the given mesh. This means
+function *shares* the same connectivity as the original mesh. This means
 that in ~pugs~ internals, there is only one connectivity object for
 these two meshes.
 #+END_note
@@ -3067,10 +3107,10 @@ different meshes produced in this example are displayed in Figure
   unset ytics
   set square
   set terminal png truecolor enhanced size 640,640
-  plot '<(sed "" $PUGS_SOURCE_DIR/doc/relax_example_m0.gnu)' lt rgb "green" w l, '<(sed "" $PUGS_SOURCE_DIR/doc/relax_example_m1.gnu)' lt rgb "blue"  w l, '<(sed "" $PUGS_SOURCE_DIR/doc/relax_example_m2.gnu)' lt rgb "black" w l
+  plot '<(sed "" $PUGS_SOURCE_DIR/doc/relax_example_m0.gnu)' lt rgb "black" w l, '<(sed "" $PUGS_SOURCE_DIR/doc/relax_example_m1.gnu)' lt rgb "blue"  w l, '<(sed "" $PUGS_SOURCE_DIR/doc/relax_example_m2.gnu)' lw 2 lt rgb "red" w l
 #+END_SRC
 
-#+CAPTION: Example of meshes relaxation. The relaxed mesh $\mathcal{M}_2$ (black) and the original meshes in green ($\mathcal{M}_0$) and blue ($\mathcal{M}_1$).
+#+CAPTION: Example of meshes relaxation. The relaxed mesh $\mathcal{M}_2$ (red) and the original meshes in black ($\mathcal{M}_0$) and blue ($\mathcal{M}_1$).
 #+NAME: fig:relax
 #+ATTR_LATEX: :width 0.38\textwidth
 #+ATTR_HTML: :width 300px;
@@ -3219,7 +3259,7 @@ operand.
 ****** ~R^1*Vh -> Vh~ and ~Vh*R^1 -> Vh~
 
 These functions are defined for $\mathbb{P}_0(\mathbb{R}^1)$ data and the
-return value is also a $\mathbb{P}_0(\mathbb{R})$ function.
+return value is a $\mathbb{P}_0(\mathbb{R})$ function.
 
 The following functions
 - ~dot: Rˆ1*Vh -> Vh~
@@ -3228,7 +3268,7 @@ The following functions
 ****** ~R^2*Vh -> Vh~ and ~Vh*R^2 -> Vh~
 
 These functions are defined for $\mathbb{P}_0(\mathbb{R}^2)$ data and the
-return value is also a $\mathbb{P}_0(\mathbb{R})$ function.
+return value is a $\mathbb{P}_0(\mathbb{R})$ function.
 
 The following functions
 - ~dot: Rˆ2*Vh -> Vh~
@@ -3237,7 +3277,7 @@ The following functions
 ****** ~R^3*Vh -> Vh~ and ~Vh*R^3 -> Vh~
 
 These functions are defined for $\mathbb{P}_0(\mathbb{R}^3)$ data and the
-return value is also a $\mathbb{P}_0(\mathbb{R})$ function.
+return value is a $\mathbb{P}_0(\mathbb{R})$ function.
 
 The following functions
 - ~dot: Rˆ3*Vh -> Vh~
@@ -3373,7 +3413,7 @@ data.
 
 ****** ~interpolate: mesh*Vh_type*(function) -> Vh~
 
-This functions takes a ~mesh~, a type of discrete function (~Vh_type~) and
+This function takes a ~mesh~, a type of discrete function (~Vh_type~) and
 a list of user functions as arguments. It returns a $\mathbb{P}_0$
 function defined at the mesh.
 
@@ -3448,11 +3488,11 @@ dimension 3.
 
 ****** ~interpolate: mesh*(zone)*Vh_type*(function) -> Vh~
 
-This function works exactly the same as the previous function. The
-additional parameter, the ~zone~ lists is used to define the cells where
-the user function (or the user function list) is interpolate. For
-cells that are not in the ~zone~ list, the discrete function is set to
-the value $0$.
+This function is similar to the previous function. The additional
+parameter, the ~zone~ list is used to define the cells where the user
+function (or the user function list) is interpolated. For cells that
+are not in the ~zone~ list, the discrete function is set to the value
+$0$.
 
 #+BEGIN_SRC pugs :exports both :results none
   import mesh;
@@ -3493,7 +3533,7 @@ dimension 3.
 ****** ~integrate: mesh*quadrature*function -> Vh~ <<integrate-classic>>
 
 This function integrates the given user function in each cell of a
-~mesh~ using a prescribed ~quadrature~. The result is of type as
+~mesh~ using a prescribed ~quadrature~. The result is of type
 $\mathbb{P}_0(\mathbb{R})$, $\mathbb{P}_0(\mathbb{R}^1)$,
 $\mathbb{P}_0(\mathbb{R}^2)$, $\mathbb{P}_0(\mathbb{R}^3)$,
 $\mathbb{P}_0(\mathbb{R}^{1\times1})$, $\mathbb{P}_0(\mathbb{R}^{2\times2})$
@@ -3511,8 +3551,8 @@ Let us consider the following example
   let U:Vh, U = integrate(m, Gauss(5), u);
 #+END_SRC
 Here, for each cell $j$, the value of the discrete function
-$\mathbf{F}_j$ is computed using a Gauss quadrature formula that is
-exact for polynomials of degree $5$, $\mathbf{F}_j \approx\int_j
+$\mathbf{U}_j$ is computed using a Gauss quadrature formula that is
+exact for polynomials of degree $5$, $\mathbf{U}_j \approx\int_j
 \mathbf{u}$. More details about quadrature formula will be given
 below.
 
@@ -3534,7 +3574,7 @@ cells.
 
 ****** ~integrate: mesh*quadrature*Vh_type*(function) -> Vh~ <<integrate-P1-vector>>
 
-This function behaves the same, the user function list size defines
+This function behaves similarly, the user function list size defines
 the dimension of the vector value of the produced
 $\vec{\mathbb{P}}_0(\mathbb{R})$ discrete function. Actually the
 ~Vh_type~ parameter is there to allow the construction of
@@ -3653,9 +3693,9 @@ sets of cells where to integrate the list of user functions.
 
 For numerical tests purpose, it is often useful to create meshes with
 random vertices positions. This is the aim of the functions that are
-described in this section. These function share some properties.
-- The generate mesh is always suitable for calculations in the sense
-  that cells volumes are warrantied to be positive.
+described in this section. These functions share some properties.
+- The generated mesh is always suitable for calculations in the sense
+  that cell volumes are warrantied to be positive.
 - Generated cells can be non-convex.
 - One has to specify boundary conditions to drive the mesh
   displacement on boundaries.
@@ -3666,7 +3706,7 @@ described in this section. These function share some properties.
 
 ****** ~randomizeMesh: mesh*(boundary_condition) -> mesh~
 
-This function creates a random mesh by displacing the nodes of a given
+This function creates a random mesh by moving the nodes of a given
 ~mesh~ and a list of ~boundary_condition~.
 
 The supported boundary conditions are the following:
@@ -3679,10 +3719,10 @@ One should refer to the section [[boundary-condition-descriptor]] for a
 documentation of the boundary condition descriptors.
 
 #+BEGIN_note
-Let us precise these boundary conditions behavior
+Let us make precise these boundary conditions behavior
 - In dimension 1, ~fixed~, ~axis~ and ~symmetry~ boundary conditions have
   the same effect.
-- In dimension 2, ~axis~ and ~symmetry~ behave the same. Thus, boundaries
+- In dimension 2, ~axis~ and ~symmetry~ behave similarly. Thus, boundaries
   supporting this kind of boundary conditions *must* form *straight*
   lines.
 - In dimension 3, boundaries describing ~axis~ conditions *must* be
@@ -3819,13 +3859,12 @@ functions may vary.
 #+END_note
 
 #+BEGIN_note
-There a three kind of boundaries are supported by ~pugs~, boundaries
-made
+There are three kinds of boundaries supported by ~pugs~, boundaries made
 - of sets of nodes,
 - of sets of edges, or
 - of sets of faces.
 
-In dimension 1, nodes, edges and faces denotes the same entities. In
+In dimension 1, nodes, edges and faces denote the same entities. In
 dimension 2, edges and faces refer the same entities.
 #+END_note
 
@@ -3837,8 +3876,8 @@ For instance, if an algorithm or a method requires a set of nodes to
 set some numerical treatment, it can be deduced from a set of faces.
 
 Obviously, these conversions can be meaningless, for instance, if one
-expects a *line* in 3d, cannot be defined by a set of faces. ~pugs~ will
-forbid this kind of conversion at runtime.
+expects a *line* in 3d, it cannot be defined by a set of faces. ~pugs~
+will forbid this kind of conversion at runtime.
 #+END_note
 
 #+BEGIN_note
@@ -3850,7 +3889,7 @@ that the given boundary is actually *straight* or *planar*.
 ~pugs~ checks that boundaries do not contain /inner/ items.
 #+END_note
 
-We regroup the boundary condition descriptors functions according to
+We regroup the boundary condition descriptor functions according to
 their arguments
 ****** ~boundary -> boundary_condition~
 - ~axis: boundary -> boundary_condition~\\
@@ -3887,7 +3926,8 @@ This function returns the quadrature descriptor associated to Gauss
 formulas for the given ~N~.
 
 In the following table, we summarize the *maximal degree* quadrature
-that are available in ~pugs~ for various elements.
+(exact for polynomials of a given degree) that are available in ~pugs~
+for various elements.
 | element type           | max. degree |
 |------------------------+-------------|
 | segment                |          23 |
@@ -3907,8 +3947,8 @@ degree given in argument.
 
 The maximum allowed degree is 23.
 
-For dimension 2 or 3 elements, Gauss-Legendre formulas are defined by
-tensorization. Conform transformations are used to map the cube
+For 2 or 3-dimensional elements, Gauss-Legendre formulas are defined
+by tensorization. Conform transformations are used to map the cube
 $]-1,1[^d$ to supported elements.
 
 ****** ~GaussLobatto: N -> quadrature~
@@ -3918,9 +3958,9 @@ degree given in argument.
 
 The maximum allowed degree is "only" 13.
 
-For dimension 2 or 3 elements, Gauss-Lobatto formulas are defined by
-tensorization. Conform transformations are used to map cube $]-1,1[^d$
-to supported elements.
+For 2 or 3-dimensional elements, Gauss-Lobatto formulas are defined by
+tensorization. Conform transformations are used to map the cube
+$]-1,1[^d$ to supported elements.
 
 ***** ~lagrangian: mesh*Vh -> Vh~
 
@@ -3994,7 +4034,7 @@ no effect.
 
 ***** Binary operators
 
-The supported binary operators for ~vh~ data types are arithmetic
+The supported binary operators for ~Vh~ data types are arithmetic
 operators.
 
 #+begin_src latex :results drawer :exports results
@@ -4091,7 +4131,7 @@ Let us consider the following example
   cout << "integral(uh)  = " << integral_of_R(uh) << "\n";
   cout << "integral(uh0) = " << integral_of_R(uh0) << "\n";
 #+END_SRC
-Here we substract the mean value of a discrete function.
+Here we subtract the mean value of a discrete function.
 #+results: substract-mean-value-to-Vh
 
 ****** Additional ~*~ operators
@@ -4105,7 +4145,7 @@ The following constructions are allowed for ~*~ operator.
   \end{equation*}
 #+end_src
 Obviously, if $\mathbb{S}=\mathbb{R}^{d\times d}$, for $d\in\{1,2,3\}$,
-the right operand must be have a compatible type, for instance
+the right operand must have a compatible type, for instance
 $\mathbb{P}_0(\mathbb{R}^{d\times d})$ or $\mathbb{P}_0(\mathbb{R}^d)$.
 
 Additionally these operators are defined
@@ -4151,8 +4191,8 @@ An important function is
 
 ****** ~getLSAvailable: void -> string~
 
-This shows the available options or libraries that are available. It
-depends on the compilation options of the code.
+This shows the available options or libraries. It depends on the
+compilation options of the code.
 
 #+NAME: get-ls-available-example
 #+BEGIN_SRC pugs :exports both :results output
@@ -4225,14 +4265,14 @@ different.
 ***** ~writer~
 
 Variables of this type manage outputs: which format is used and
-eventually the writing policy. This policy sets for instance the time
+possibly the writing policy. This policy sets for instance the time
 period for time-dependent post processing.
 
 **** ~writer~ provided functions
 
 ***** ~name_output: Vh*string -> output~
 
-This function give a name to a discrete function.
+This function gives a name to a discrete function.
 
 #+BEGIN_SRC pugs :exports both :results none
   import mesh;
@@ -4359,7 +4399,7 @@ of ~fh~.
 
 ***** ~name_output: item_value*string -> output~
 
-This function give a name to an ~item_value~.
+This function gives a name to an ~item_value~.
 
 #+BEGIN_SRC pugs :exports both :results none
   import mesh;
@@ -4395,14 +4435,14 @@ This function give a name to an ~item_value~.
 
 #+BEGIN_warning
 One observes that we supplied a mesh to the writer function. The
-reason for that is that ~item_value~ refers a connectivity. It requires
+reason for that is that ~item_value~ refers to a connectivity. It requires
 a ~mesh~ with the same connectivity to be written.
 
-However, is a discrete function (of type ~Vh~) build on the same
+However, if a discrete function (of type ~Vh~) build on the same
 connectivity is provided, there is no need to specify the ~mesh~.
 #+END_warning
 
-Let us illustrate it by an important second example.
+Let us give an example.
 
 #+BEGIN_SRC pugs :exports both :results none
   import mesh;
@@ -4456,15 +4496,15 @@ dimension 1 and 2 (~gnuplot_writer~).
 Both of these writers can be used for single output or time series
 outputs. In the case of single output, the filename is completed by
 adding the extension ~.gnu~, in the case of time series, the filename is
-extending by adding ~.abcd.gnu~, where ~abcd~ is the number of the output
-in the series.
+extending by adding ~.abcd.gnu~, where ~abcd~ is the number in the output
+series.
 
 #+BEGIN_note
 The ~gnuplot~ writers are implemented in parallel.
 
-The ~gnuplot~ post processing of produced files is the same whichever is
-the number of processors (as soon as the saved data is also the same,
-which is warrantied by ~pugs~ for explicit methods).
+The ~gnuplot~ post processing of produced files does not depend on the
+number of processors (as soon as the saved data is also the same,
+which is ensured by ~pugs~ for explicit methods).
 #+END_note
 
 For an obvious practical reason, each ~gnuplot~ file starts with a
@@ -4487,11 +4527,11 @@ Here is an example of preamble of a produced ~gnuplot~ file.
 
 ****** ~gnuplot_1d_writer~ functions
 
-This writer family make only sense in 1d.
+This writer family makes sense only in 1d.
 
 #+BEGIN_note
 In parallel, as soon as the saved data themselves are the same, the
-~gnuplot_1d_writer~ generates *exactly* the same files (whichever is the
+~gnuplot_1d_writer~ generates *exactly* the same files (whichever the
 number of processors) since the coordinates of the post processed data
 are sorted according to growing abscissa.
 #+END_note
@@ -4570,8 +4610,8 @@ A typical use of this writer is the following.
 ******* ~gnuplot_1d_writer: string*R -> writer~ <<gnuplot-1d-series>>
 
 This writer differs from the previous one by handling output
-series. The real value argument defines the period to respect between
-two outputs. It can be viewed as an helper to outputs.
+series. The real value argument defines the period between two
+outputs. It can be viewed as helper to outputs.
 
 Let us give an example to fix ideas.
 #+BEGIN_SRC pugs :exports both :results none
@@ -4613,8 +4653,8 @@ Running this example produces the following files
 #+END_SRC
 #+RESULTS: ls-produced-gp-1d-series
 
-Each of these file contains the numerical solution at following saving
-times:
+Each of these files contains the numerical solution at following
+saving times:
 #+NAME: times-in-gp-1d-series
 #+BEGIN_SRC shell :exports results :results output
   grep -n "# time = " gp_1d_exp_sin.*.gnu | cut -d '=' -f 2
@@ -4625,9 +4665,9 @@ times:
 
 These writers differ from the previous ones since it draws the cells
 and affects the cell value to its nodes. This produces larger files
-but allows 2d representations. Also, if the saved data is exactly the
-same in parallel, the order of the cells is generally different since
-they are written processor by processor.
+but allows for 2d representations. Also, if the saved data is exactly
+the same in parallel, the order of the cells is generally different
+since they are written processor by processor.
 
 Additionally, this writer allows to write 2d meshes, see paragraph
 [[write-mesh]].
@@ -4665,7 +4705,7 @@ Figure [[fig:writer-gp-sin]].
   plot '<(sed "" $PUGS_SOURCE_DIR/doc/gp_sin.gnu)' lw 2 w lp
 #+END_SRC
 
-#+CAPTION: Example of produced gnuplot results from the ~gnuplot_writer~. One can compare ths produced result to the one of the ~gnuplot_1d_writer~ given in Figure [[fig:writer-gp-1d-sin]]
+#+CAPTION: Example of produced gnuplot results from the ~gnuplot_writer~. One can compare the produced result to the one of the ~gnuplot_1d_writer~ given in Figure [[fig:writer-gp-1d-sin]]
 #+NAME: fig:writer-gp-sin
 #+ATTR_LATEX: :width 0.38\textwidth
 #+ATTR_HTML: :width 300px;
@@ -4713,7 +4753,7 @@ The gnuplot result is displayed on Figure [[fig:writer-gp-2d-cos-sin]].
 ******* ~gnuplot_writer: string*R -> writer~
 
 This is the time series function in the case of the ~gnuplot_writer~. It
-behaves the same as [[gnuplot-1d-series]].
+behaves similarly to [[gnuplot-1d-series]].
 
 ***** ~vtk~ writers
 
@@ -4721,8 +4761,8 @@ For more complex post processing (including 3d), ~pugs~ can generate ~vtk~
 outputs.
 
 The used format is one file in the ~vtu~ format for each parallel domain
-(and eventually each time). The output is done using binary data for
-performance reasons. For each time step a ~pvtu~ file is generate to
+(and each time). The output is produced using binary data for
+performance reasons. For each time step a ~pvtu~ file is generated to
 handle parallelism. And for a complete time series, a ~pvd~ file is
 produced. This is the file that should be loaded.
 
@@ -4730,7 +4770,7 @@ Observe that each of these files (~vtu~, ~pvtu~ and ~pvd~) contains a
 comment that stores the creation date and the version of ~pugs~ that was
 run.
 
-The use is exactly the same than for ~gnuplot~ writers so we do not
+The use is exactly the same as for ~gnuplot~ writers so we do not
 provide full examples.
 
 ~vtk~ writers are compatible with the ~write_mesh~ function, see paragraph
@@ -4744,10 +4784,9 @@ produced ~pvd~ file is built by adding ~.pvd~ to the provided ~string~.
 ****** ~vtk_writer: string*R -> writer~
 
 This function follows the same rule. One just specifies the output
-period. The generated ~pvd~ file is built the same way, one adds ~.pvd~ to
+period. The generated ~pvd~ file is built similarly, one adds ~.pvd~ to
 the provided ~string~.
 
-
 ***** ~write~, ~force_write~ and ~write_mesh~ functions
 
 Once a mesh writer has been defined, these functions are called to
@@ -4755,9 +4794,9 @@ effectively generate the post processing files.
 
 ****** ~write: writer*(output) -> void~
 
-As parameters, it takes a ~writer~ that is not a time series one and a
-list of ~output~ (which are named discrete functions that are defined on
-the *same* mesh).
+As parameters, it takes a ~writer~ (that is not a time series writer)
+and a list of ~output~ (which are named discrete functions that are
+defined on the *same* mesh).
 
 We have already shown a lot of examples of use of the ~write~
 function. So let us focus on improper calls.
@@ -4854,7 +4893,7 @@ share the same connectivity with the ~mesh~.
 One probably noticed that using the ~write~ function with a time series
 ~writer~, last time of the calculation may not be written (see section
 [[gnuplot-1d-series]]). The ~force_write~ function does not check that the
-saving time has been reached. It just checks that the current time has
+saving time has been reached. It only checks that the current time has
 not already been saved.
 
 Let us improve slightly the example given in section
@@ -4904,8 +4943,8 @@ Running this example produces the following files
 #+RESULTS: ls-produced-gp-1d-series-force
 One can see the additional file.
 
-Each of these file contains the numerical solution at following saving
-times:
+Each of these files contains the numerical solution at following
+saving times:
 #+NAME: times-in-gp-1d-series-force
 #+BEGIN_SRC shell :exports results :results output
   grep -n "# time = " gp_1d_exp_sin_force.*.gnu | cut -d '=' -f 2
diff --git a/src/analysis/CMakeLists.txt b/src/analysis/CMakeLists.txt
index 723083e448194938b0f0a54245bc53d9dd0176bb..ffab3ba49306f99fbb272da67322c189a1ff86b5 100644
--- a/src/analysis/CMakeLists.txt
+++ b/src/analysis/CMakeLists.txt
@@ -13,3 +13,10 @@ add_library(
   TensorialGaussLobattoQuadrature.cpp
   TetrahedronGaussQuadrature.cpp
   TriangleGaussQuadrature.cpp)
+
+if (CMAKE_CXX_COMPILER_ID STREQUAL "GNU")
+  if((CMAKE_CXX_COMPILER_VERSION VERSION_GREATER_EQUAL "12.0.0") AND (CMAKE_CXX_COMPILER_VERSION VERSION_LESS "13.0.0"))
+    # Deactivated since it produces false positive warning in this file only ...
+    set_source_files_properties(PyramidGaussQuadrature.cpp PROPERTIES COMPILE_FLAGS "-Wno-array-bounds")
+  endif()
+endif()
diff --git a/src/language/modules/BuiltinModule.hpp b/src/language/modules/BuiltinModule.hpp
index 33afb27d718be626cad0e14da1736761e5f8a624..e15dda2972686761a52211eebf2cd82fa87f8825 100644
--- a/src/language/modules/BuiltinModule.hpp
+++ b/src/language/modules/BuiltinModule.hpp
@@ -6,6 +6,7 @@
 
 #include <utils/Exceptions.hpp>
 
+#include <functional>
 #include <sstream>
 
 class IBuiltinFunctionEmbedder;
diff --git a/src/language/utils/ASTNodeDataType.hpp b/src/language/utils/ASTNodeDataType.hpp
index b7cc48fefaa6d849cf970641b37279a5a85ca15e..4ef3a80c188faa5898fd29e0bc35cce651c7f6c8 100644
--- a/src/language/utils/ASTNodeDataType.hpp
+++ b/src/language/utils/ASTNodeDataType.hpp
@@ -3,6 +3,7 @@
 
 #include <utils/PugsAssert.hpp>
 
+#include <array>
 #include <limits>
 #include <memory>
 #include <string>
@@ -119,7 +120,7 @@ class ASTNodeDataType
   }
 
   ASTNodeDataType& operator=(const ASTNodeDataType&) = default;
-  ASTNodeDataType& operator=(ASTNodeDataType&&) = default;
+  ASTNodeDataType& operator=(ASTNodeDataType&&)      = default;
 
   template <DataType data_type>
   [[nodiscard]] static ASTNodeDataType
diff --git a/tools/pgs-pygments.py b/tools/pgs-pygments.py
new file mode 100755
index 0000000000000000000000000000000000000000..a6622c0dc1b25aa8f998a2480af0f5d080c4a5d5
--- /dev/null
+++ b/tools/pgs-pygments.py
@@ -0,0 +1,42 @@
+from pygments.lexer import *
+from pygments.token import *
+
+class PugsLexer(RegexLexer):
+    name = 'Pugs'
+    aliases = ['pgs']
+    filenames = ['*.pgs']
+
+    tokens = {
+        'root': [
+            (r'([Ll]?)(")', bygroups(String.Affix, String.Double), 'string'),
+            (r'([Ll]?)(\')(\\[^\']+)(\')', bygroups(String.Affix, String.Char, String.Escape, String.Char)),
+            (r'\s+', Whitespace),
+            (r'/\*', Comment.Multiline, 'comment'),
+            (r'//.*?$', Comment.Singleline),
+            (r'/', Text),
+            (words(('B', 'N', 'Z', 'R', 'R^1', 'R^2','R^3', 'R^1x1', 'R^2x2', 'R^3x3', 'string'), suffix=r'\b'), Keyword.Type),
+            (words(('let', 'import', 'do', 'while', 'for', 'if', 'else' 'break'), suffix=r'\b'), Keyword.Reserved),
+            (words(('and', 'or', 'xor', 'not', 'true', 'false'), suffix=r'\b'), Keyword.Reserved),
+            (words(('cout', 'cerr', 'clog'), suffix=r'\b'), Name.Variable),
+            (r'[,~!%&*+=|?:<>/-]', Operator),
+            (r'[+-]?\d+(\.\d*)?[Ee][+-]?\d+', Number.Float),
+            (r'[+-]?(\d+\.\d*)|(\d*\.\d+)([Ee][+-]?\d+)?', Number.Float),
+            (r'\d+[LlUu]*', Number.Integer),
+            (r'[\(\)\[\];.]', Punctuation),
+            (r'[{}]', Punctuation),
+            (r'\w+', Name),
+        ],
+        'comment': [
+            (r'[^*/]+', Comment.Multiline),
+            (r'/\*', Comment.Multiline, '#push'),
+            (r'\*/', Comment.Multiline, '#pop'),
+            (r'[*/]', Comment.Multiline)
+        ],
+        'string': [
+            (r'"', String, '#pop'),
+            (r'\\([\\abfnrtv"\']|x[a-fA-F0-9]{2,4}|[0-7]{1,3})', String.Escape),
+            (r'\\\n', String.Escape),  # line continuation
+            (r'[^\\"]+', String),  # all other characters
+            (r'\\', String),  # stray backslash
+        ],
+    }
diff --git a/tools/pgs-pygments.sh b/tools/pgs-pygments.sh
new file mode 100755
index 0000000000000000000000000000000000000000..98f90c55e028546aa8fa49fc069f06b1fe3f057a
--- /dev/null
+++ b/tools/pgs-pygments.sh
@@ -0,0 +1,28 @@
+#!/bin/env bash
+
+arguments=()
+
+while [[ "$#" -ne 0 ]]; do
+  case "$1" in
+  -l)
+    arguments+=("$1")
+
+    shift
+    lexer=$1
+
+    if [[ "$lexer" == "Pugs" ]]; then
+      arguments+=("${PUGS_SOURCE_DIR}/tools/pgs-pygments.py:PugsLexer")
+      arguments+=("-x")
+    else
+      arguments+=("${lexer}")
+    fi
+
+    ;;
+  *)
+    arguments+=("$1")
+    ;;
+  esac
+  shift
+done
+
+pygmentize -Ostyle=colorful "${arguments[@]}"