diff --git a/CircuitCoupling/R_circuit.geo b/CircuitCoupling/R_circuit.geo
index dad9c9a86d34162a7e00eb758a892b0d4e78c395..1d8e7e30e988e5bedc346d6e167130ff92a08a6c 100644
--- a/CircuitCoupling/R_circuit.geo
+++ b/CircuitCoupling/R_circuit.geo
@@ -2,7 +2,7 @@
 
 Include "R_circuit_common.pro";
 
-lc = 0.25;      // Cjaracteristic length
+lc = 0.25;      // Characteristic length
 
 Point(1) = {0, 0, 0, lc};
 Extrude {1, 0, 0} {
diff --git a/Elasticity/wrench2D.pro b/Elasticity/wrench2D.pro
index 76147ed1b51383cea809f24f305debff127c7ec0..cd70662385c43098b5a1786396c2d3849de74fa9 100644
--- a/Elasticity/wrench2D.pro
+++ b/Elasticity/wrench2D.pro
@@ -1,35 +1,38 @@
 /* -------------------------------------------------------------------
-   Tutorial 3 : linear elastic model of a wrench
+    Tutorial 3 : linear elastic model of a wrench
 
-   Features:
-   - "grad u" GetDP specific formulation for linear elasticity
-   - first and second order elements
-   - triangular and quadrangular elements
+    Features:
+    - "grad u" GetDP specific formulation for linear elasticity
+    - first and second order elements
+    - triangular and quadrangular elements
 
-   To compute the solution interactively from the Gmsh GUI:
+    To compute the solution in a terminal:
+       getdp wrench2D.pro -solve Elast_u -pos pos
+
+    To compute the solution interactively from the Gmsh GUI:
        File > Open > wrench.pro
        Run (button at the bottom of the left panel)
    ------------------------------------------------------------------- */
 
-/* Linear elasticity with GetDP:
-
-   GetDP has a peculiar way to deal with linear elasticity.  Instead of a vector
-   field "u = Vector[ ux, uy, uz ]", the displacement field is regarded as two
-   (2D case) or 3 (3D case) scalar fields. Unlike conventional formulations
-   then, GetDP's formulation is written in terms of the gradient "grad u" of the
-   displacement field, which is a non-symmetric tensor, and the needed
-   symmetrization (to define the strain tensor and relate it to the stress
-   tensor) is done through the constitutive relationship (Hooke law).  The
-   reason for this unusual formulation is to be able to use also for elastic
-   problems the powerful geometrical and homological kernel of GetDP, which
-   relies on the operators grad, curl and div.
-
-   The "grad u" formulation entails a small increase of assembly work but makes
-   in counterpart lots of geometrical features implemented in GetDP (change of
-   coordinates, function spaces, etc...)  applicable to elastic problems
-   out-of-the-box, since the scalar fields { ux, uy, uz } have exactly the same
-   geometrical properties as, e.g. a scalar electric potential or a temperature
-   field. */
+/* 
+    Particularities of linear elasticity in GetDP:
+
+    Instead of a vector field "u = Vector[ ux, uy, uz ]", the displacement field 
+    is regarded as two (2D case) or three (3D case) scalar fields. 
+    Unlike conventional formulations, GetDP formulation is then written in terms 
+    of the gradient "grad u" of the displacement field, which is a non-symmetric 
+    tensor, and the needed symmetrization (to define the strain tensor and relate 
+    it to the stress tensor) is done via the constitutive relationship (Hooke law). 
+    This unusual formulation allows to take advantage of the powerful geometrical 
+    and homological GetDP kernel, which relies on the operators grad, curl and div.
+    
+    The "grad u" formulation entails a small increase of assembly work but makes
+    in counterpart lots of geometrical features implemented in GetDP (change of
+    coordinates, function spaces, etc...) applicable to elastic problems
+    out-of-the-box, since the scalar fields { ux, uy, uz } have exactly the same
+    geometrical properties as, e.g. an electric scalar potential or a temperature
+    field. 
+*/
 
 Include "wrench2D_common.pro";
 
@@ -54,24 +57,30 @@ Group {
   Sur_Clamp_Mec = Region[ Grip ];
   Sur_Force_Mec = Region[ Force ];
 
-  /* Signification of the abstract regions:
+  /* Meaning of abstract regions:
      - Vol_Mec       : Elastic domain
-     - Vol_Force_Mec : Region with imposed volumic force
+     - Vol_Force_Mec : Region with imposed volume force
      - Sur_Force_Mec : Surface with imposed surface traction
      - Sur_Clamp_Mec : Surface with imposed zero displacements (all components) */
 }
 
 Function {
-  /* Material coefficients.
-     No need to define them regionwise here ( E[{Wrench}] = ... ; )
-     as there is only one region in this model. */
+  /* 
+    Material coefficients
+    No need of regionwise definition ( E[{Wrench}] = ... ; )
+    as this model comprises only one region. 
+    If there is more than one region and coefficients are not particularised, 
+    the same value holds for all of them. 
+ */
   E[] = Young;
   nu[] = Poisson;
-  /* Components of the volumic force applied to the region "Vol_Force_Mec"
-     Gravity could be defined here ( force_y[] = 7000*9.81; ) ; */
+  /* 
+    Volume force components applied to the region "Vol_Force_Mec"
+    Gravity could be defined here as well ( force_y[] = 7000*9.81; ) ; 
+  */
   force_x[] = 0;
   force_y[] = 0;
-  /* Components of the surface traction force applied to the region "Sur_Force_Mec" */
+  /* Surface traction force components applied to the region "Sur_Force_Mec" */
   pressure_x[] = 0;
   pressure_y[] = -AppliedForce/(SurfaceArea[]*Thickness); // downward vertical force
 }
@@ -82,7 +91,7 @@ Function {
 
    sigma_ij = C_ijkl epsilon_ij
 
-   is represented in 2D by 4 2x2 tensors C_ij[], i,j=1,2, depending on the Lamé
+   is represented in 2D by four 2x2 tensors C_ij[], i,j=1,2, depending on the Lamé
    coefficients of the isotropic linear material,
 
    lambda = E[]*nu[]/(1.+nu[])/(1.-2.*nu[])
@@ -94,7 +103,8 @@ Function {
    EPD:  a[] = lambda + 2 mu     b[] = mu     c[] = lambda
     3D:  a[] = lambda + 2 mu     b[] = mu     c[] = lambda
 
-   respectively for the 2D plane strain (EPD), 2D plane stress (EPS) and 3D cases. */
+   respectively for the 2D plane strain (EPD), 2D plane stress (EPS) and 3D cases. 
+*/
 
 Function {
   Flag_EPC = 1;
@@ -128,7 +138,7 @@ Constraint {
   }
 }
 
-/* As explained above, the displacement field is discretized as two scalar
+/* As mentioned above, the displacement field is discretized as two scalar
    fields "ux" and "uy", which are the spatial components of the vector field
    "u" in a fixed Cartesian coordinate system.
 
@@ -143,17 +153,17 @@ Constraint {
    u . n = ux Cos [th] + uy Sin [th] = ... ;
 
    are less naturally accounted for within the "grad u" formulation; but they
-   could be easily implemented with a Lagrange multplier.
+   could be easily implemented with e.g. a Lagrange multiplier.
 
-   Finite element shape (triangles or quadrangles) makes no difference in the
-   definition of the FunctionSpaces. The appropriate shape functions to be used
+   The finite element shape (e.g. triangles or quadrangles in 2D) has no influence 
+   in the definition of the FunctionSpaces. The appropriate shape functions
    are determined by GetDP at a much lower level on basis of the information
    contained in the *.msh file.
 
-   Second order elements, on the other hand, are implemented in the hierarchical
-   fashion by adding to the first order node-based shape functions a set of
-   second order edge-based functions to complete a basis for 2d order
-   polynomials on the reference element. */
+   Second order elements are hierarchically implemented by adding to the first 
+   order node-based shape functions a set of second order edge-based functions 
+   to complete a basis for 2D order polynomials on the reference element. 
+*/
 
 // Domain of definition of the "ux" and "uy" FunctionSpaces
 Group {
@@ -217,8 +227,8 @@ Jacobian {
   }
 }
 
-/* Adapt the number of Gauss points to the polynomial degree of the finite elements
-   is as simple as this: */
+/* Adapting the number of Gauss points to the polynomial degree of the finite elements
+   is simple: */
 Integration {
   { Name Gauss_v;
     Case {
diff --git a/Electrostatics/microstrip.pro b/Electrostatics/microstrip.pro
index a5b4f33089c0e6dd5fa7116243a88d4b4a06e573..91d96e44d411b2c77b8caf1c23d1462a2e4ed5b3 100644
--- a/Electrostatics/microstrip.pro
+++ b/Electrostatics/microstrip.pro
@@ -36,12 +36,12 @@
    boundary condition (zero flux of the displacement field, i.e. n.d = 0) is
    imposed on the left boundary of the domain to account for the symmetry of the
    problem, as well as on the top and right boundaries that truncate the
-   modelling domain. */
+   simulation domain. */
 
 Group {
   /* One starts by giving explicit meaningful names to the Physical regions
-     defined in the "microstrip.msh" mesh file. We only use 2 volume regions and
-     2 surface regions in this model. */
+     defined in the "microstrip.msh" mesh file. This model comprises only 
+     2 volume regions and 2 surface regions. */
 
   Air = Region[101];
   Diel1 = Region[111];
@@ -61,6 +61,10 @@ Group {
 
      Since there are no non-homogeneous Neumann conditions in this particular
      example, Sur_Neu_Ele is defined as empty.
+     
+     Note that volume elements are those that correspond to the higher dimension 
+     of the model at hand (2D elements here), surface elements correspond to the 
+     higher dimension of the model minus one (1D elements here).
      */
 
   Vol_Ele = Region[ {Air, Diel1} ];
@@ -78,8 +82,8 @@ Function {
 }
 
 Constraint {
-  /* As for material laws, the Dirichlet boundary condition is defined
-     piecewise.  The constraint "Dirichlet_Ele" is invoked in the FunctionSpace
+  /* The Dirichlet boundary condition is also defined piecewise.  
+     The constraint "Dirichlet_Ele" is invoked in the FunctionSpace
      below. */
 
   { Name Dirichlet_Ele; Type Assign;
@@ -97,13 +101,13 @@ Group{
      Contrary to Vol_xxx and Sur_xxx regions, which contain only volume or
      surface regions, resp., domains of definition Dom_xxx regions may contain
      both volume and surface regions. Hence the use of the prefixes Vol_, Sur_
-     and Dom_ to avoid confusions.*/
+     and Dom_ to avoid confusion.*/
 
   Dom_Hgrad_v_Ele =  Region[ {Vol_Ele, Sur_Neu_Ele} ];
 }
 
 FunctionSpace {
-  /* The function space in which we shall pick the electric scalar potential "v"
+  /* The function space in which we pick the electric scalar potential "v"
      solution is defined by
      - a domain of definition (the "Support": "Dom_Hgrad_v_Ele")
      - a type ("Form0" means scalar field)
@@ -127,7 +131,7 @@ FunctionSpace {
       { Name sn; NameOfCoef vn; Function BF_Node;
         Support Dom_Hgrad_v_Ele; Entity NodesOf[ All ]; }
       // using "NodesOf[All]" instead of "NodesOf[Dom_Hgrad_v_Ele]" is an
-      // optimization, which allows GetDP to not explicitly build the list of
+      // optimization, which avoids explicitly building the list of 
       // all the nodes
     }
     Constraint {
@@ -138,7 +142,7 @@ FunctionSpace {
 
 Jacobian {
   /* Jacobians are used to specify the mapping between elements in the mesh and
-     the reference elements (defined in standardized unit cells) over which
+     the reference elements (defined in standardized/reference unit cells) over which
      integration is performed. "Vol" represents the classical 1-to-1 mapping
      between identical spatial dimensions, i.e. in this case a reference
      triangle/quadrangle onto triangles/quadrangles in the z=0 plane (2D <->
@@ -173,7 +177,7 @@ Integration {
 }
 
 Formulation {
-  /* The syntax of the Formulation section is a harder nut to crack. So let's
+  /* The syntax of the Formulation section is a harder nut to crack. So let us
      deal with it carefully.
 
      A GetDP Formulation encodes a so-called weak formulation of the original
@@ -201,11 +205,11 @@ Formulation {
      ones used to interpolate the unknown potential v.
 
      The "Integral" statement in the Formulation is a symbolic representation of
-     this weak formulation. It has got 4 semicolon separated arguments:
+     this weak formulation. It has 4 semicolon separated arguments:
      * the density [.,.] to be integrated (note the square brackets instead of
-       the parentheses), with the test-functions (always) after the comma
-     * the domain of integration,
-     * the Jacobian of the transformation reference element <-> real element,
+       the parentheses), with the test-functions (always) after the comma;
+     * the domain of integration;
+     * the Jacobian of the transformation reference element <-> real element;
      * the integration method to be used.
 
      In the density, braces around a quantity (such as "{v}") indicate that this
@@ -223,33 +227,33 @@ Formulation {
      be something like this [ ... , basis_functions_of {d v} ]. However, since
      the second term is always devoted to test functions, the operator
      "basis_functions_of" would always be there.  It can therefore be made
-     implicit, and, according to the GetDP syntax, it is omitted. So, one writes
-     simply [ ... , {d v} ].
+     implicit and, according to the GetDP syntax, omitted. So, one simply writes
+     [ ... , {d v} ].
 
-     The first term, on the other hand, can contain a much wider variety of
+     On the other hand, the first term can contain a much wider variety of
      expressions than the second one. In our case it should be expressed in
      terms of the FE expansion of "v" at the present system solution, i.e. when
      the coefficients vn_k in the expansion of "v = Sum_k vn_k sn_k" are
-     unknown. This is indicated by prefixing the braces with "Dof", which leads
-     to the following density:
+     unknown. This is indicated by prefixing the braces with "Dof" (degrees of 
+     freedom), which leads to the following density:
 
      [ epsilon[] * Dof{d v} , {d v} ],
 
      a so-called bilinear term that will contribute to the stiffness matrix of
      the electrostatic problem at hand.
 
-     Another option, which would not work here, is to evaluate the first
-     argument with the last available already computed solution, i.e. simply
-     perform the interpolation with known coefficients vn_k. For this the Dof
-     prefix would be omitted and we would have:
+     Another option, which does not work here, is to evaluate the first
+     argument with the last available computed solution, i.e. simply
+     perform the interpolation with known coefficients vn_k. The Dof
+     prefix is then omitted and we have:
 
      [ epsilon[] * {d v} , {d v} ],
 
-     a so-called linear term that will contribute to the right-hand side of the
+     a so-called linear term that contributes to the right-hand side of the
      linear system.
 
-     Both choices are commonly used in different contexts, and we shall come
-     back on this often in subsequent tutorials. */
+     Both choices are commonly used in different contexts, and we shall often 
+     come back in subsequent tutorials. */
 
   { Name Electrostatics_v; Type FemEquation;
     Quantity {
@@ -259,20 +263,21 @@ Formulation {
       Integral { [ epsilon[] * Dof{d v} , {d v} ];
 	    In Vol_Ele; Jacobian Vol; Integration Int; }
 
-   /* Additional Integral terms could be added here. For example, the
-      following term would account for non-homogeneous Neumann boundary
+   /* Additional Integral terms can be added here. For example, the
+      following term may account for non-homogeneous Neumann boundary
       conditions, provided that the function nd[] is defined:
 
       Integral { [ nd[] , {v} ];
 	    In Sur_Neu_Ele; Jacobian Sur; Integration Int; }
 
-      All the terms are added, and an implicit "= 0" is considered at the end. */
+      All the terms in the Equation environment are added, 
+      and an implicit "= 0" is considered at the end. */
     }
   }
 }
 
-/* What to do with a weak formulation is specified in the Resolution: here we
-   simply generale a linear system, solve it and save the solution (.res file)
+/* In the Resolution environment we specify what to do with a weak formulation: 
+   here we simply generate a linear system, solve it and save the solution (.res file)
    to disk. */
 
 Resolution {
@@ -291,18 +296,18 @@ Resolution {
    The first part defines, in terms of the Formulation, which itself refers to
    the FunctionSpace, a number of quantities that can be evaluated at the
    postprocessing stage. The three quantities defined here are:
-   - the scalar electric potential,
-   - the electric field,
+   - the scalar electric potential;
+   - the electric field;
    - the electric displacement.
 
    The second part consists in defining groups of post-processing operations,
    which can be invoked separately. The first group is invoked by default when
    Gmsh is run interactively. Each Operation specifies
-   - a quantity to be eveluated,
-   - the region on which the evaluation is done,
+   - a quantity to evaluate;
+   - the region on which the evaluation is done;
    - the name of the output file.
    The generated post-processing files are automatically displayed by Gmsh if
-   the "Merge result automatically" option is enabled (which is the default). */
+   the "Merge result automatically" option is enabled (the default). */
 
 PostProcessing {
   { Name EleSta_v; NameOfFormulation Electrostatics_v;
diff --git a/ElectrostaticsFloating/floating.pro b/ElectrostaticsFloating/floating.pro
index 50daaa566e9883856e604b23c982fa03b46ffc35..cbe37d6f78f8d8018c90d2a0d05bbf4397f35213 100644
--- a/ElectrostaticsFloating/floating.pro
+++ b/ElectrostaticsFloating/floating.pro
@@ -12,7 +12,7 @@
 
    ------------------------------------------------------------------- */
 
-/* A thing GetDP is pretty good at is the management of global (non-local) basis
+/* GetDP is pretty good at the management of global (non-local) basis
    functions. Finite element expansions typically associate basis functions to
    individual nodes or edges in the mesh. But consider the situation where a
    scalar field is set to be uniform over a region of the problem (say a
@@ -20,12 +20,12 @@
    identical nodal value "v_electrode", a global (non-local) basis function
    "BF_electrode" is obtained as factor which is the sum of the shape functions
    of all the nodes in the electrode region. This basis function "BF_electrode"
-   - is a continuous function, scalar in this case,
-   - is equal to 1 at the nodes of the electrode, and to 0 at all other nodes,
+   - is a continuous function, scalar in this case;
+   - is equal to 1 at the nodes of the electrode, and to 0 at all other nodes;
    - decreases from 1 to 0 over the one-element-thick layer of elements sharing
      at least one node with the electrode region.
 
-   One such glabal basis function can be associated with each electrode in the
+   One such global basis function can be associated with each electrode in the
    system, so that the finite element expansion of the electric scalar potential
    reads:
 
@@ -34,13 +34,13 @@
    with the the sum_k running over all nodes except those of the electrode
    regions.
 
-   We show in this tutorial how GetDP takes advantage of global quantities and
+   This tutorial shows how GetDP takes advantage of global quantities and
    the associated global basis functions
-   - to reduce the number of unknowns
-   - to compute efficiently the electrode charges "Q_electrode", which are
-     precisely the energy duals of the global "v_electrode" quantities
+   - to reduce the number of unknowns;
+   - to efficiently compute the electrode charges "Q_electrode", which are
+     precisely the energy duals of the global "v_electrode" quantities;
    - to deal with floating potentials, which are the computed electrode
-     potential when the electrode charge is imposed
+     potentials when the electrode charge is imposed;
    - to provide output quantities (charges, armature voltages, capacitances,
      ...) that can be immediately used in a external circuit. */
 
@@ -57,8 +57,8 @@ Group {
   /* Abstract regions:
 
      Vol_Ele            : volume where -div(epsilon grad v) = 0 is solved
-     Sur_Neu_Ele        : surface where non homogeneous Neumann boundary conditions
-                          (on n.d = -n . (epsilon grad v)) are imposed
+     Sur_Neu_Ele        : surface with imposed non homogeneous Neumann boundary conditions
+                          (on n.d = -n . (epsilon grad v))
      Sur_Electrodes_Ele : electrode regions */
 
   Vol_Ele = Region[ {Air, Diel1} ];
@@ -66,7 +66,7 @@ Group {
   Sur_Electrodes_Ele = Region [ {Ground, Microstrip} ];
 }
 
-/* A number of ONELAB parameters are defined to define model parameters or model
+/* A number of ONELAB parameters are defined to provide model parameters or model
    options interactively. */
 
 MicrostripTypeBC = DefineNumber[0, Name "1Microstrip excitation/Type",
@@ -86,7 +86,7 @@ Function {
 
 Constraint {
   /* The Dirichlet boundary condition on the local electric potential is no
-     longer used. The microstrip and the ground are now treated as electrodes,
+     longer used. The microstrip and the ground are herein treated as electrodes,
      whose voltage is imposed with the "SetGlobalPotential" constraint below. */
   { Name Dirichlet_Ele; Type Assign;
     Case {
@@ -95,12 +95,12 @@ Constraint {
 
   { Name SetGlobalPotential; Type Assign;
     Case {
-      /* Define the imposed potential regionwise on the different parts of
+      /* Impose the potential regionwise on the different parts of
 	 "Sur_Electrodes_Ele". No voltage is imposed to the Microstrip electrode
 	 when the "Fixed charge" option is enabled (if MicrostripTypeBC != 0). */
       { Region Ground; Value 0; }
       If(!MicrostripTypeBC)
-	{ Region Microstrip; Value MicrostripValueBC; }
+	    { Region Microstrip; Value MicrostripValueBC; }
       EndIf
     }
   }
@@ -108,26 +108,26 @@ Constraint {
     Case {
       /* Impose the charge if MicrostripTypeBC != 0 */
       If(MicrostripTypeBC)
-	{ Region Microstrip; Value MicrostripValueBC; }
+	    { Region Microstrip; Value MicrostripValueBC; }
       EndIf
     }
   }
 }
 
 Group{
-  /* The domain of definition lists all regions on which the field "v" is
+  /* The domain of definition comprises all regions on which the field "v" is
      defined.*/
   Dom_Hgrad_v_Ele =  Region[ {Vol_Ele, Sur_Neu_Ele, Sur_Electrodes_Ele} ];
 }
 
 FunctionSpace {
-  /* The magic in the treatment of global quantitities by GetDP is in the fact
+  /* The magic in the treatment of global quantitities by GetDP lies in the fact
      that nearly all the work is done at the level of the FunctionSpace
      definition. The finite element expansion of "v" is
 
      v = Sum_k sn_k vn_k + Sum_electrode v_electrode BF_electrode
 
-     with the the sum_k running over all nodes except those of the electrode
+     with the sum_k running over all nodes except those of the electrode
      regions. This is exactly what one finds in the FunctionSpace definition
      below with "sf" standing for "BF_electrode" and "vf" for "v_electrode".
 
@@ -135,10 +135,10 @@ FunctionSpace {
      "GlobalQuantity" section; these names are used in the corresponding
      "GlobalTerm" in the Formulation. Such global terms are the equivalent of a
      "Integral" term, but where no integration needs to be performed. The
-     "AssociatedWith" statement manifests the fact that the global potential of
+     "AssociatedWith" statement refers to the fact that the global potential of
      an electrode is the (electrostatic) energy dual of the electric charge
      carried by that electrode. Indeed, let us consider the electrostatic weak
-     formulation derived in Tutorial 1: find v in Hgradv_Ele such that
+     formulation derived in Tutorial 1: find v in Hgrad_v_Ele such that
 
      (epsilon grad v, grad v')_Vol_Ele + (n . (epsilon grad v), v')_Bnd_Vol_Ele = 0
 
@@ -153,14 +153,14 @@ FunctionSpace {
      charge Q_electrode carried by the electrodes.
 
      By checking the "Display global basis functions" checkbox and running the
-     model, you can take a look on how the two "BF_electrode" basis functions in
+     model, you can take a look at how the two "BF_electrode" basis functions in
      this model look like.  Constraints can then be set on either component of
      the FunctionSpace.  Besides the usual Dirichlet boundary condition on the
-     local field, which is left here for the sake of completeness but is not
-     used in this model, there is the possibility to fix either the
-     GlobalPotential or the ArmatureCharge of each indidual electrode (not both,
-     of course).  When the ArmatureCharge is fixed, the computed GlobalPotential
-     computed for that electrode is the so-called floating potential. */
+     local field (left here for the sake of completeness but not
+     used in this model), one may fix either the
+     GlobalPotential or the ArmatureCharge of each indidual electrode (never both,
+     of course). When the ArmatureCharge is fixed, the computed GlobalPotential
+     for that electrode is the so-called floating potential. */
 
   { Name Hgrad_v_Ele; Type Form0;
     BasisFunction {
@@ -207,8 +207,8 @@ Integration {
 }
 
 Formulation {
-  /* The formulation only contains minor changes compared to formulation from
-     the first tutorial.  The global quantities are declared as "Global" in the
+  /* The formulation contains only minor changes compared to formulation from
+     the first tutorial. The global quantities are declared as "Global" in the
      "Quantity" section, and a "GlobalTerm" is added that triggers the assembly
      of the additional equation per electrode (the "pre-integrated" boundary
      term) in the system to compute the charge Q_electrode, which
@@ -290,8 +290,8 @@ PostProcessing {
   }
 }
 
-/* Various output results are generated, which are both displayed in the
-   graphical user interface, and stored in disk files.  In particular, global
+/* Several output results are generated, which are both displayed in the
+   graphical user interface, and stored in disk files. In particular, global
    quantities related results are stored in the "output.txt" file. A user option
    allows to chose to not overwrite the "output.txt" file when running a new
    simulation. */
diff --git a/MagneticForces/InfiniteBox.geo b/MagneticForces/InfiniteBox.geo
index 821ab74e6658a51568233e737ae7dc4466cda322..b3a13268c555f672cd1f75467cda7d56d46f39da 100644
--- a/MagneticForces/InfiniteBox.geo
+++ b/MagneticForces/InfiniteBox.geo
@@ -117,7 +117,7 @@ p6 = newp; Point (p6) = {xCnt+xInt, yCnt-yInt, zCnt+zInt, lc1inf};
 p7 = newp; Point (p7) = {xCnt+xInt, yCnt+yInt, zCnt+zInt, lc1inf};
 p8 = newp; Point (p8) = {xCnt-xInt, yCnt+yInt, zCnt+zInt, lc1inf};
 
-
+VolInf[]={}; // Define empty array in case Flag_InfiniteBox is not active
 If(Flag_InfiniteBox)
   xExt = xInt * ratioInf;
   yExt = yInt * ratioInf;
diff --git a/MagneticForces/magnets.pro b/MagneticForces/magnets.pro
index 7e73f2ddf8d7ef2f4a17929ed55d2b6c0c361282..b1703c2169ccf47356b33d7f270fac8ed44f9316 100644
--- a/MagneticForces/magnets.pro
+++ b/MagneticForces/magnets.pro
@@ -8,35 +8,39 @@
    - Maxwell stress tensor and rigid-body magnetic forces
 
    To compute the solution in a terminal:
-       getdp magnets.pro
+   First generate the (3D) mesh and then run getdp with the chosen resolution 
+       gmsh magnets.geo -3
+       getdp magnets.pro -solve MagSta_a
+       OR
+       getdp magnets.pro -solve MagSta_phi
 
    To compute the solution interactively from the Gmsh GUI:
        File > Open > magnets.pro
+       Resolution can be chosen from the menu on the left: 
+       MagSta_a (default) or MagSta_phi
        Run (button at the bottom of the left panel)
    ------------------------------------------------------------------- */
 
 /*
- This rather didactic tutorial solves the electromagnetic field 
+ This tutorial solves the electromagnetic field 
  and the rigid-body forces acting on a set of magnetic pieces
- of parallelepipedic or cylindrical shape.
+ of either parallelepipedic or cylindrical shape.
  Besides position and dimension, each piece is attributed 
- a (constant) magnetic permability and/or a remanence field.
- The pieces are all (a bit abusively) generically called "Magnet" 
- in the problem decription below, irresective of whether they are
- truly permanent magnets or ferromagnetic barrels. 
-
- The tutorial model proposes both dual 3D magnetostatic formulations:
- the magnetic vector potential formulation with spanning-tree gauging,
- and the scalar magnetic potential formulation.
- The later is rather simple in this case since, as there are no conductors,
- the known coercive field hc[] is the only source field "hs" one needs 
- to represent the magnetic field:
-  
+ a (constant) magnetic permeability and/or a remanence field.
+ Hereafter, the pieces are all, simply though imprecisely, referred to as "Magnet", 
+ irresective of whether they are truly permanent magnets or ferromagnetic barrels. 
+
+ The tutorial model proposes two dual 3D magnetostatic formulations:
+ - the magnetic vector potential formulation with spanning-tree gauging;
+ - the scalar magnetic potential formulation.
+ As there are no conductors, the later is rather simple. The source field "hs" is
+ directly the the known coercive field hc[]: 
+   
    h = hs - grad phi   ,  hs = -hc.
 
- As in tutorial 2 (magnetostatic field of an electromagnet), a shell
- of so-called infinite elements is used here to impose the exact
- zero-field boundary condition at infinity.
+ If the "Add infinite box" box is ticked, a transformation to infinity shell is 
+ used to impose the exact zero-field boundary condition at infinity. 
+ See also Tutorial 2: magnetostatic field of an electromagnet.
  The shell is generated automatically by including "InfiniteBox.geo"
  at the end of the geometrical description of the model. 
  It can be placed rather close of the magnets without loss of accuracy.
@@ -44,24 +48,23 @@
  The preferred way to compute electromagnetic forces in GetDP
  is as an explicit by-product of the Maxwell stress tensor "TM[{b}]",
  which is a material dependent function of the magnetic induction "b" field. 
- Exactly like we computed the heat flux "q(S)" through a surface "S"
- using a special auxiliary function "g(S)" associated with that surface 
- in "Tutorial 5 : thermal problem with contact resistances",
- the magnetic force acting on a rigid body in empty space can be evaluated
- as the flux of the Maxwell stress tensor through that surface.
- There is one auxiliary function "g(SkinMagnet~{i}) = un~{i}"
- for each magnet, and the resultant magnetic force acting on "Magnet~{i}"
- is given by the integral:
+ The magnetic force acting on a rigid body in empty space can be evaluated
+ as the flux of the Maxwell stress tensor through a surface "S" (surrounding the body).
+ A special auxiliary function "g(S)" linked "S" is defined for each magnet, i.e.
+ "g(SkinMagnet~{i}) = un~{i}".
+ The resultant magnetic force acting on "Magnet~{i}" is given by the integral:
 
  f~{i} = Integral [ TM[{b}] * {-grad un~{i}} ] ;
-
- It should be insisted on the fact that the Maxwell stress tensor
- is always discontinuous on material discontinuities, 
+ 
+ This approach is analogous to the computation of heat flux "q(S)" through a 
+ surface "S" described in "Tutorial 5: thermal problem with contact resistances".
+ 
+ Note that the Maxwell stress tensor is always discontinuous on material discontinuities, 
  and that magnetic forces acting on rigid bodies
- only depend on the Maxwell stress tensor of empty space, 
+ depend only on the Maxwell stress tensor in empty space, 
  and on the "b" and "h" field distribution, 
  on the external side of "SkinMagnet~{i}" 
- (the side of the surface in contact with air).
+ (side of the surface in contact with air).
 
  "{-grad un~{i}}" in the above formula can be regarded 
  as the normal vector to "SkinMagnet~{i}"
@@ -79,7 +82,7 @@
 Include "magnets_common.pro"
 
 DefineConstant[
-  // preset all getdp options and make them invisible
+  // preset all getdp options and make them (in)visible
   R_ = {"MagSta_a", Name "GetDP/1ResolutionChoices", Visible 1,
 	Choices {"MagSta_a", "MagSta_phi"}},
   C_ = {"-solve -v 5 -v2 -bin", Name "GetDP/9ComputeCommand", Visible 0}
@@ -198,7 +201,7 @@ Constraint {
 }
 
 FunctionSpace {
-  { Name Hgrad_phi ; Type Form0 ; // scalar magnetic potential
+  { Name Hgrad_phi ; Type Form0 ; // magnetic scalar potential
     BasisFunction {
       { Name sn ; NameOfCoef phin ; Function BF_Node ;
         Support Dom_Hgrad_phi ; Entity NodesOf[ All ] ; }
@@ -207,7 +210,7 @@ FunctionSpace {
       { NameOfCoef phin ; EntityType NodesOf ; NameOfConstraint phi ; }
     }
   }
-  { Name Hcurl_a; Type Form1; // vector magnetic potential
+  { Name Hcurl_a; Type Form1; // magnetic vector potential
     BasisFunction {
       { Name se;  NameOfCoef ae;  Function BF_Edge; 
 	Support Dom_Hcurl_a ;Entity EdgesOf[ All ]; }
diff --git a/Magnetodynamics/electromagnet.pro b/Magnetodynamics/electromagnet.pro
index b5e7138e09ad332ffd43db3f16ba3594bfc02bfd..a381e598b2fa389d36ee9badc349d12dc647c848 100644
--- a/Magnetodynamics/electromagnet.pro
+++ b/Magnetodynamics/electromagnet.pro
@@ -14,6 +14,8 @@
 
    To compute the solution interactively from the Gmsh GUI:
        File > Open > electromagnet.pro
+       You may choose the Resolution in the left panel: 
+       Magnetodynamics2D_av (default) or Magnetostatics2D_a
        Run (button at the bottom of the left panel)
    ------------------------------------------------------------------- */
 
@@ -32,8 +34,8 @@ Group {
   Vol_Mag = Region[{Air, Core, Ind, AirInf}]; // full magnetic domain
   Vol_C_Mag = Region[Core]; // massive conductors
   Vol_S_Mag = Region[Ind]; // stranded conductors (coils)
-  Vol_Inf_Mag = Region[AirInf]; // annulus for infinite shell transformation
-  Val_Rint = rInt; Val_Rext = rExt; // interior and exterior radii of annulus
+  Vol_Inf_Mag = Region[AirInf]; // ring-shaped shell for infinite transformation
+  Val_Rint = rInt; Val_Rext = rExt; // interior and exterior radii of ring
 }
 
 Function {
@@ -50,7 +52,7 @@ Function {
   sigma[ Ind ] = 5e7;
 
   Ns[ Ind ] = 1000 ; // number of turns in coil
-  Sc[ Ind ] = SurfaceArea[] ; // surface (cross section) of coil
+  Sc[ Ind ] = SurfaceArea[] ; // area of coil cross section
   // Current density in each coil portion for a unit current (will be multiplied
   // by the actual total current in the coil)
   js0[ Ind ] = Ns[]/Sc[] * Vector[0,0,-1];
@@ -66,7 +68,7 @@ Constraint {
   }
   { Name Current_2D;
     Case {
-      // represents the phasor amplitude for a dynamic analysis
+      // represents the phasor amplitude (peak to peak value) for a dynamic analysis
       { Region Ind; Value Current; }
     }
   }
diff --git a/Magnetostatics/electromagnet.pro b/Magnetostatics/electromagnet.pro
index 066ada9edb5cae64c536067f1433c678e26d1453..4261455468d6a29f9867713205cf8285d2d48b20 100644
--- a/Magnetostatics/electromagnet.pro
+++ b/Magnetostatics/electromagnet.pro
@@ -29,20 +29,19 @@
    the reluctivity.
 
    Electromagnetic fields expand to infinity. The corresponding boundary
-   condition can be imposed rigorously by means of a gometrical transformation
+   condition can be imposed rigorously by means of a geometrical transformation
    that maps a ring (or shell) of finite elements to the complementary of its
    interior.  As this is a mere geometric transformation, it is enough in the
    model description to attribute a special Jacobian to the ring region
    ("AirInf") - see the "Jacobian" section below.  With this information, GetDP
-   is able to deal with the correct transformation of all quantities involved in
-   the model.
+   is able to correctly transform all quantities involved in the model.
 
    The special Jacobian "VolSphShell" takes 2 parameters in this case,
    "Val_Rint" and "Val_Rext", which represent the inner and outer radii of the
-   transformed ring region and whose value must match those used in the
+   transformed ring region and the value of which must match those used in the
    geometrical description of the model (.geo file).  This is a typical case
    where Gmsh and GetDP must consistently share parameter values.  To ensure
-   consistency in all cases, common parameters are defined is a specific file
+   consistency, common parameters are defined is a specific file
    "electromagnet_common.pro", which is included in both the .geo and .pro file
    of the model.
 
@@ -72,8 +71,8 @@ Group {
 
      The abstract regions in this model have the following interpretation:
      - Vol_Mag     : full volume domain
-     - Vol_S_Mag   : region where the current source js is defined
-     - Vol_Inf_Mag : region where the infinite ring geometric transformation is applied
+     - Vol_S_Mag   : region with imposed current source js
+     - Vol_Inf_Mag : region where the infinite ring geometric transformation applies
      - Sur_Dir_Mag : part of the boundary with homogenous Dirichlet conditions
      - Sur_Neu_Mag : part of the boundary with non-homogeneous Neumann conditions
   */
@@ -85,7 +84,7 @@ Group {
 }
 
 /* The weak formulation for this problem is derived in a similar way to the
-   electrostatic weak formulation from Tutorial 1. The main difference is that
+   electrostatic weak formulation from Tutorial 1. The main differences are that
    the fields are now vector-valued, and that we have one linear (source) term
    in addition to the bilinear term. The weak formulation reads: find a such
    that
@@ -109,13 +108,13 @@ Group {
 Function {
   mu0 = 4.e-7 * Pi;
 
-  /* New ONELAB variables can then be defined using defineNumber, e.g.: */
+  /* New ONELAB variables can then be defined using DefineNumber, e.g.: */
   murCore = DefineNumber[100, Name "Model parameters/Mur core",
 			 Help "Magnetic relative permeability of Core"];
 
   /* When the script is run, if the parameter named "Model parameters/Mur core"
      has not been previously defined, it takes the value (100) provided in
-     defineNumber and is sent to the ONELAB server. The "/" character in the
+     DefineNumber and is sent to the ONELAB server. The "/" character in the
      variable name is interpreted as a path separator, and results in the
      creation of a sub-tree in the graphical user interface. If the script is
      re-run later, the value will be updated using the value from the server
@@ -150,7 +149,7 @@ Constraint {
 }
 
 /* In the 2D approximation, the magnetic vector potential a and the current
-   density js are vectors with a z-component only, i.e.:
+   density js are vectors with only the z-component, i.e.:
 
    a  = Vector [ 0, 0, az(x,y) ]
    js = Vector [ 0, 0, jsz(x,y) ]
@@ -162,7 +161,7 @@ Constraint {
    with "BF_PerpendicularEdge" basis functions associated to nodes of the mesh)
    and in the "Hregion_j_Mag_2D" FunctionSpace for js ("Vector" with
    "BF_RegionZ" basis functions associated with the region Vol_S_Mag). Without
-   giving all the details, a is thus node-based, whereas js is region-wise
+   providing explicit details, a is thus node-based, whereas js is region-wise
    constant. */
 
 Group {
@@ -232,7 +231,7 @@ Integration {
    Integral { [ - Dof{js} , {a} ];
      In Vol_S_Mag; Jacobian Vol; Integration Int; }
 
-   The chosen implementation below is however more effeicient as it avoids
+   The chosen implementation below is however more efficient as it avoids
    evaluating repeatedly the function js_fct[] during assembly.
 */