diff --git a/Electrostatics/microstrip.geo b/Electrostatics/microstrip.geo
index ee8c8dfc688e969c151b430d542de3b1fdf8d0b1..04ecbf515d566ca77ed14566b1be84e24bb52706 100644
--- a/Electrostatics/microstrip.geo
+++ b/Electrostatics/microstrip.geo
@@ -1,8 +1,9 @@
 /* -------------------------------------------------------------------
    File "microstrip.geo"
 
-   This file is the geometrical description used by Gmsh to produce
-   the mesh file "microstrip.msh" that is read by GetDP.
+   This file is the geometrical description used by Gmsh to produce the mesh
+   file "microstrip.msh", which is read by GetDP when analysing the file
+   "microstrip.pro".
    ------------------------------------------------------------------- */
 
 /* Definition of some parameters for geometrical dimensions, i.e.  h (height of
@@ -40,19 +41,18 @@ Line(11) = {5,7};
 
 /* Definition of geometrical surfaces */
 
-Line Loop(12) = {1, 2, -8, -9, 7};   Plane Surface(13) = {12};
-Line Loop(14) = {10,-11,8,3,4,5}; Plane Surface(15) = {14};
+Curve Loop(12) = {1, 2, -8, -9, 7}; Plane Surface(13) = {12};
+Curve Loop(14) = {10,-11,8,3,4,5};  Plane Surface(15) = {14};
 
-/* Definition of Physical entities (surfaces, lines).  The definition of
-   Physical entities (Surfaces and Lines) tells Gmsh the elements and associated
-   region numbers that have to be saved in the mesh file 'microstrip.msh'.  For
-   example, Region 111 is made of the triangle elements of the geometric surface
-   13, whereas Region 121 is made of line elements of the geometric lines 9, 10
-   and 11. */
+/* Definition of Physical entities.  The definition of Physical entities
+   (Surfaces and Curves) tells Gmsh the elements and associated region numbers
+   that have to be saved in the mesh file 'microstrip.msh'.  For example, Region
+   111 is made of the triangle elements of the geometric surface 13, whereas
+   Region 121 is made of line elements of the geometric curves 9, 10 and 11. */
 
-Physical Surface ("Air", 101) = {15};
-Physical Surface ("Dielectric", 111) = {13};
+Physical Surface("Air", 101) = {15};
+Physical Surface("Dielectric", 111) = {13};
 
-Physical Line ("Ground", 120) = {1} ;
-Physical Line ("Electrode", 121) = {9,10,11} ;
-Physical Line ("Surface infinity", 130) = {2,3,4} ;
+Physical Curve("Ground", 120) = {1} ;
+Physical Curve("Electrode", 121) = {9,10,11} ;
+Physical Curve("Surface infinity", 130) = {2,3,4} ;
diff --git a/Electrostatics/microstrip.pro b/Electrostatics/microstrip.pro
index e7b1a988a4604807c37026b49c7d6494d8dca42d..bc81df1a78115f26504a1ca6801bef632bcd5abe 100644
--- a/Electrostatics/microstrip.pro
+++ b/Electrostatics/microstrip.pro
@@ -16,48 +16,67 @@
        Run (button at the bottom of the left panel)
    ------------------------------------------------------------------- */
 
+/* In this first tutorial we consider the solution of electric fields given a
+   static distribution of electric potentials. This corresponds to a so-called
+   "electrostatic" physical model, obtained by combining the time-invariant
+   Maxwell-Ampere equation (Curl e = 0, with e the electric field) with Gauss'
+   law (Div d = rho, with d the displacement field and rho the charge density)
+   and the dielectric constitutive law (d = epsilon e, with epsilon the
+   dielectric permittivity).
+
+   Since Curl e = 0, e can be derived from a scalar electric potential v, such
+   that e = -Grad v. Plugging this potential in Gauss' law and using the
+   constitutive law leads to a scalar Poisson equation in terms of the electric
+   potential: -Div(epsilon Grad v) = rho.
+
+   We consider here the special case where rho = 0, to model a conducting
+   microstrip line on top of a dielectric substrate. A Dirichlet boundary
+   condition sets the potential to 1 mV on the boundary of the line (called
+   "Electrode" below) and 0 V on the ground. A homogeneous Neumann boundary
+   condition (zero flux of the displacement field, i.e. n.d = 0) is imposed on a
+   surface truncating the modelling domain. */
+
 Group {
-  /* One starts by giving explicit meaningful names to
-     the Physical regions defined in the "microstrip.msh" mesh file.
-     There are 2 volume regions and 3 surface regions in this model. */
+  /* One starts by giving explicit meaningful names to the Physical regions
+     defined in the "microstrip.msh" mesh file. There are 2 volume regions and 3
+     surface regions in this model. */
 
   Air = Region[101];
   Diel1 = Region[111];
 
   Ground = Region[120];
   Electrode = Region[121];
-  SurfInf = Region[130];
 
-  /* We now define abstract regions to be used below
-     in the definition of the scalar electric potential formulation:
+  /* We now define abstract regions to be used below in the definition of the
+     scalar electric potential formulation:
 
-     Vol_Dielectric_Ele : dielectric volume regions where
-                          "Div ( epsr[] Grad v)" is solved
-     Sur_Dir_Ele        : Dirichlet boundary condition (v imposed)
-     Sur_Neu_Ele        : Neumann bondary condition ( epsr[] n.Grad v = 0 )
+     Vol_Ele        : volume where -Div(epsilon Grad v) = 0 is solved
+     SurNeumann_Ele : surface where non homogeneous Neumann boundary conditions
+                      (on n.d == epsilon n.Grad v) are imposed
 
-     Vol_xxx groups contain only volume elements of the mesh (triangles here).
-     Sur_xxx groups contain only surface elements of the mesh (lines here).
-  */
+     Vol* groups contain only volume elements of the mesh (triangles here).
+     Sur* groups contain only surface elements of the mesh (lines here). */
 
-  Vol_Dielectric_Ele = Region[ {Air, Diel1} ];
-  Sur_Dir_Ele = Region[ {Ground, Electrode} ];
-  Sur_Neu_Ele = Region[ {SurfInf} ];
+  Vol_Ele = Region[ {Air, Diel1} ];
+  SurNeumann_Ele = Region[ {} ]; // empty for this model
 }
 
 Function {
-  /* Material laws (here the relative permittivity)
-     are defined piecewise in terms of the above defined physical regions */
-
-  epsr[Air] = 1.;
-  epsr[Diel1] = 9.8;
+  /* Constants are evaluated (only once) when the .pro file is parsed by GetDP,
+     before any finite element calculation takes place. */
+  eps0 = 8.854187818e-12; // permittivity of empty space
+
+  /* Material laws (here the dielectric permittivity) are defined as piecewise
+     functions (note the square brackets), in terms of the above defined
+     physical regions. */
+  epsilon[Air] = 1. * eps0;
+  epsilon[Diel1] = 9.8 * eps0;
 }
 
 Constraint {
-  /* As for material laws, the Dirichlet boundary condition
-     is defined piecewise.
-     The constraint "Dirichlet_Ele" is invoked in the FunctionSpace below */
-
+  /* As for material laws, the Dirichlet boundary condition is defined
+     piecewise.  The constraint "Dirichlet_Ele" is invoked in the FunctionSpace
+     below. */
   { Name Dirichlet_Ele; Type Assign;
     Case {
       { Region Ground; Value 0.; }
@@ -67,53 +86,60 @@ Constraint {
 }
 
 Group{
-  /* The domain of definition of a FunctionSpace lists all regions
-     on which a field is defined.
+  /* The domain of definition of a FunctionSpace lists all regions on which a
+     field is defined.
 
-     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.*/
+     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.*/
 
-  Dom_Hgrad_v_Ele =  Region[ {Vol_Dielectric_Ele, Sur_Dir_Ele, Sur_Neu_Ele} ];
+  Dom_Hgrad_v_Ele =  Region[ {Vol_Ele, SurNeumann_Ele} ];
 }
 
 FunctionSpace {
-  /* The function space in which we shall pick
-     the electric scalar potential "v" solution
-     is definied by
-     - a domain of definition ("Dom_Hgrad_v_Ele")
+  /* The function space in which we shall pick the electric scalar potential "v"
+     solution is definied by
+     - a domain of definition (the "Support": "Dom_Hgrad_v_Ele")
      - a type ("Form0" means scalar field)
      - a set of scalar basis functions ("BF_Node" means nodal basis functions)
+     - a set of entities to which the basis functions are associated ("Entity":
+       here all the nodes of the domain of definition "NodesOf[All]")
      - a constraint (here the Dirichlet boundary conditions)
 
      The FE expansion of the unknown field "v" reads
 
-     v = Sum_k vn_k sn_k
+     v(x,y) = Sum_k vn_k sn_k(x,y)
 
-     where the "vn_k" are the nodal values (connectors)
-     and "sn_k" the basis functions.
-     Not all connectors are unknowns of the FE problem,
-     due to the "Constraint", which assigns particular values
-     to the nodes of the region "Sur_Dir_Ele".
-     GetDP deals with that automatically
-     on basis of the definition of the FunctionSpace.
+     where the "vn_k" are the nodal values (connectors) and "sn_k(x,y)" the
+     nodal basis functions.  Not all connectors are unknowns of the FE problem,
+     due to the "Constraint", which assigns particular values to the nodes of
+     the Ground and Electrode regions. GetDP deals with that automatically on
+     basis of the definition of the FunctionSpace.
   */
-
   { Name Hgrad_v_Ele; Type Form0;
     BasisFunction {
       { 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
+      // all nodes
     }
     Constraint {
-      { NameOfCoef vn; EntityType NodesOf;
-        NameOfConstraint Dirichlet_Ele; }
+      { NameOfCoef vn; EntityType NodesOf; NameOfConstraint Dirichlet_Ele; }
     }
   }
 }
 
 Jacobian {
+  /* Jacobians are used to specify the mapping between elements in the mesh and
+     the reference elements (in a standard unit-type cell) over which
+     integration is performed. "Vol" represents the classical 1-to-1 mapping
+     between elements of identical geometrical dimension, i.e. in this case
+     triangles/quadrangles in the plane onto the reference triangle/quadrangle
+     (2D to 2D). "Sur" would be used to map triangles/quadrangles in a 3D volume
+     onto the reference triangle/quadrangle (3D to 2D). "Lin" would be used to
+     map segments in a volume to the unit segment (3D to 1D). */
   { Name Vol ;
     Case {
       { Region All ; Jacobian Vol ; }
@@ -122,95 +148,104 @@ Jacobian {
 }
 
 Integration {
+  /* A Gauss quadrature rule with 4 points is used for all integrations below. */
   { Name Int ;
-    Case { {Type Gauss ;
-            Case { { GeoElement Triangle    ; NumberOfPoints  4 ; }
-                   { GeoElement Quadrangle  ; NumberOfPoints  4 ; } }
+    Case { { Type Gauss ;
+             Case { { GeoElement Triangle    ; NumberOfPoints  4 ; }
+                    { GeoElement Quadrangle  ; NumberOfPoints  4 ; } }
       }
     }
   }
 }
 
 Formulation {
-  /* The syntax of the Formulation{} section is a hard nut to crack.
-     So let's deal with it carefully.
+  /* The syntax of the Formulation section is a harder nut to crack. So let's
+     deal with it carefully.
+
+     A GetDP Formulation encodes a so-called weak formulation of the original
+     partial differential equation, i.e. of -Div(epsilon Grad v) = 0. This weak
+     formulation involves finding v such that
+
+     (-Div(epsilon Grad v) , v')_Vol_Ele = 0
+
+     holds for all so-called "test-functions" v', where (.,.) denotes an inner
+     product. If the test-functions v' are differentiable, integration by parts
+     using Green's identity leads to finding v such that
+
+     (epsilon Grad v, Grad v')_Vol_Ele + (epsilon n.Grad v, v')_SurNeumann_Ele = 0
 
-     The weak form of the electrostatic problem is particularly simple
-     in this model, as it has only one term:
+     holds for all v'.
 
-     (epsr[] Grad v, Grad vp)_Vol_Dielectric_Ele = 0
-     for all vp in S0(Vol_Dielectric_Ele).
+     In this particular example the Neumann boundary term vanishes and we are
+     looking for functions v in the function space Hgrad_v_Ele such that
 
-     The corresponding Euler-Lagrange equations are:
-     * Div ( epsr[] Grad v) = 0  on Vol_Dielectric_Ele
-     * epsr[] n.Grad v = 0       on Sur_Neu_Ele
+     (epsilon Grad v, Grad v')_Vol_Ele = 0
+
+     holds for all v', where the test functions v' are the same basis functions
+     ("sn_k") as the ones used to interpolate the unknown potential v.
+
+     The Galerkin statement in the Formulation is a symbolic representation of
+     this weak formulation.
 
-     The Galerkin{} statement is a symbolic representation
-     of this weak formulation term.
      It has got 4 semicolon separated arguments:
-     * the density [,] to be integrated,
+     * the density [.,.] to be integrated, with the test-functions (always) on
+       the right side of the comma
      * the domain of integration,
-     * the jacobian of the transformation reference element -> real element,
+     * the Jacobian of the transformation reference element <-> real element,
      * the integration method to be used.
-     The symbol "d" represents the exterior derivative,
-     and it is a synonym of "Grad" when applied to a scalar function.
-     The expression "d v" stands thus for the gradient of the v field,
-     i.e., for the electric field up to a minus sign.
 
-     What is a bit confusing is that the two comma-separated terms
-     of the bracket [,], in the first argument,
-     are not interpreted the same way. Let us unravel this in detail.
+     In the density, braces around a quantity (such as "{v}") indicate that this
+     quantity belongs to a FunctionSpace. Differential operators can be applied
+     within braces (such as "{Grad v}"); in particular the symbol "d" represents
+     the exterior derivative, and it is a synonym of "Grad" when applied to a
+     scalar function (such as "{d v}").
 
-     As the Galerkin method uses as trial (test) functions
-     the basis functions "sn_k" of the unknown field "v",
-     the density should be something like this
+     What can be a bit confusing is that the two comma-separated terms of the
+     density are not interpreted exactly in the same way. Let us unravel this in
+     detail.  As the Galerkin method uses as test functions the basis functions
+     "sn_k" of the unknown field "v", the density should be something like this
 
-     [ epsr[] * {d v} , basis_functions_of {d v} ].
+     [ epsilon[] * {d v} , basis_functions_of {d v} ].
 
-     Since the second argument is devoted to the trial 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}".
+     Since the second argument 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}".
 
-     The first argument, on the other hand, can contain
-     a much wider variety of expressions than the second one.
-     In this problem, it contains
+     The first argument, on the other hand, can contain a much wider variety of
+     expressions than the second one. In this problem, it contains
 
-     epsr[] * {d v} = Sum_k   vn_k  epsr[]*{d sn_k}
+     epsilon[] * {d v} = Sum_k vn_k epsilon[]*{d sn_k}
 
-     which is the eletric displacement vector, up to a sign.
-     Here, we have two valid syntaxes, with very different meanings.
+     which is the eletric displacement vector, up to a sign.  Here, we have two
+     valid syntaxes, with very different meanings.
 
-     The first argument can be expressed in terms of
-     the FE expansion of "v" at the present system solution.
-     This is indicated by invoking the Dof{} operator:
+     The first argument can be expressed in terms of the FE expansion of "v" at
+     the present system solution, i.e. when the coefficients vn_k are unknown.
+     This is indicated by prefixing the braces with "Dof":
 
-     [ epsr[] * Dof{d v} , {d v} ].
+     [ epsilon[] * Dof{d v} , {d v} ].
 
-     Another option, which would not work here,
-     is to evaluate the first argument
-     with the last available already computed solution.
-     For this the Dof{} operator is omitted:
+     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 operator would be omitted:
 
-     [ epsr[] * {d v} , {d v} ].
+     [ epsilon[] * {d v} , {d v} ].
 
-     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 come
+     back on this often in subsequent tutorials.
 
-     To express the stiffness matrix of the electrostatic problem at hand,
-     we have to take the first option.
-     Hence the final expression of the density below.
-
-   */
+     To express the stiffness matrix of the electrostatic problem at hand, we
+     have to take the first option.  Hence the final expression of the density
+     below. */
   { Name Electrostatics_v; Type FemEquation;
     Quantity {
       { Name v; Type Local; NameOfSpace Hgrad_v_Ele; }
     }
     Equation {
-      Galerkin { [ epsr[] * Dof{d v} , {d v} ];
-	In Vol_Dielectric_Ele;
-	Jacobian Vol; Integration Int; }
+      Galerkin { [ epsilon[] * Dof{d v} , {d v} ];
+	In Vol_Ele; Jacobian Vol; Integration Int; }
     }
   }
 }
@@ -226,27 +261,23 @@ Resolution {
   }
 }
 
-
-eps0 = 8.854187818e-12;  // permittivity of empty space
-
 /* Post-processing is done in two parts.
-   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 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 vector 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
+   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,
    - 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 generated post-processing files are automatically displayed by Gmsh if
+   the "Merge result automatically" option is enabled (which is the default). */
 
 PostProcessing {
   { Name EleSta_v; NameOfFormulation Electrostatics_v;
@@ -254,7 +285,7 @@ PostProcessing {
       { Name v;
         Value {
           Local { [ {v} ];
-	    In Dom_Hgrad_v_Ele; Jacobian Vol; }
+            In Dom_Hgrad_v_Ele; Jacobian Vol; }
         }
       }
       { Name e;
@@ -265,7 +296,7 @@ PostProcessing {
       }
       { Name d;
         Value {
-          Local { [ -eps0*epsr[] * {d v} ];
+          Local { [ -epsilon[] * {d v} ];
 	    In Dom_Hgrad_v_Ele; Jacobian Vol; }
         }
       }
@@ -281,10 +312,6 @@ PostOperation {
      Operation {
        Print [ v, OnElementsOf Dom_Hgrad_v_Ele, File "mStrip_v.pos" ];
        Print [ e, OnElementsOf Dom_Hgrad_v_Ele, File "mStrip_e.pos" ];
-       Echo[ StrCat["l=PostProcessing.NbViews-1;",
-		    "View[l].IntervalsType = 3;",
-		    "View[l].NbIso = 40;"],
-	     File "tmp.geo", LastTimeStepOnly] ;
        Print [ e, OnLine {{e,h,0}{14.e-3,h,0}}{60}, File "Cut_e.pos" ];
      }
   }