diff --git a/Elasticity/wrench2D.pro b/Elasticity/wrench2D.pro
index 5751251be66825a17fbf0c4fcbe1c0ca66c2b668..a1f8dc128ec65847db1bbf821b4e924a4d4fb221 100644
--- a/Elasticity/wrench2D.pro
+++ b/Elasticity/wrench2D.pro
@@ -2,7 +2,7 @@
    Tutorial 3 : linear elastic model of a wrench
 
    Features:
-   - "Grad u" GetDP specific formulation for linear elasticity
+   - "grad u" GetDP specific formulation for linear elasticity
    - first and second order elements
    - triangular and quadrangular elements
 
@@ -16,15 +16,15 @@
    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
+   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.
+   relies on the operators grad, curl and div.
 
-   The "Grad u" formulation entails a small increase of assembly work but makes
+   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
@@ -142,7 +142,7 @@ Constraint {
 
    u . n = ux Cos [th] + uy Sin [th] = ... ;
 
-   are less naturally accounted for within the "Grad u" formulation; but they
+   are less naturally accounted for within the "grad u" formulation; but they
    could be easily implemented with a Lagrange multplier.
 
    Finite element shape (triangles or quadrangles) makes no difference in the
diff --git a/Electrostatics/microstrip.pro b/Electrostatics/microstrip.pro
index 3ac146210b73d2af440d6fafb4ae4953e9259fc0..43bf1dad6e653db345aa0bf64607ccfaf5d73fd3 100644
--- a/Electrostatics/microstrip.pro
+++ b/Electrostatics/microstrip.pro
@@ -19,15 +19,15 @@
 /* In this first tutorial we consider the calculation of the electric field
    given a static distribution of electric potential. This corresponds to an
    "electrostatic" physical model, obtained by combining the time-invariant
-   Faraday 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)
+   Faraday 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
+   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.
+   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
@@ -52,9 +52,9 @@ Group {
   /* We now define abstract regions to be used below in the definition of the
      scalar electric potential formulation:
 
-     Vol_Ele    : volume where -Div(epsilon Grad v) = 0 is solved
+     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
+                  (on n.d = -n . (epsilon 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).
@@ -172,24 +172,24 @@ Formulation {
      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
+     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
+     (-div(epsilon grad v) , v')_Vol_Ele = 0
 
      holds for all so-called "test-functions" v', where (.,.)_D denotes an inner
      product over a domain D. 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 + (n . (epsilon Grad v), v')_Bnd_Vol_Ele = 0
+     (epsilon grad v, grad v')_Vol_Ele + (n . (epsilon grad v), v')_Bnd_Vol_Ele = 0
 
      holds for all v', where Bnd_Vol_Ele is the boundary of Vol_Ele. In our
      microstrip example this surface term vanishes, as there is either no test
-     function v' (on the Dirichlet boundary), or "n. (epsilon Grad v)" is zero
+     function v' (on the Dirichlet boundary), or "n. (epsilon grad v)" is zero
      (on the homogeneous Neumann boundary). We are thus eventually looking for
      functions v in the function space Hgrad_v_Ele such that
 
-     (epsilon Grad v, Grad v')_Vol_Ele = 0
+     (epsilon grad v, grad v')_Vol_Ele = 0
 
      holds for all v'. Finally, our choice here is to use a Gakerkin method,
      where the test functions v' are the same basis functions ("sn_k") as the
diff --git a/ElectrostaticsFloating/floating.pro b/ElectrostaticsFloating/floating.pro
index 58f6b19aacd239fc116bbe20bfcdb57ec107ab4f..dbbfe60b1d561ac8fa653389a6e39e6510b94341 100644
--- a/ElectrostaticsFloating/floating.pro
+++ b/ElectrostaticsFloating/floating.pro
@@ -56,9 +56,9 @@ Group {
 
   /* Abstract regions:
 
-     Vol_Ele            : volume where -Div(epsilon Grad v) = 0 is solved
+     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
+                          (on n.d = -n . (epsilon grad v)) are imposed
      Sur_Electrodes_Ele : electrode regions */
 
   Vol_Ele = Region[ {Air, Diel1} ];
@@ -140,15 +140,15 @@ FunctionSpace {
      carried by that electrode. Indeed, let us consider the electrostatic weak
      formulation derived in Tutorial 1: find v in Hgradv_Ele such that
 
-     (epsilon Grad v, Grad v')_Vol_Ele + (n . (epsilon Grad v), v')_Bnd_Vol_Ele = 0
+     (epsilon grad v, grad v')_Vol_Ele + (n . (epsilon grad v), v')_Bnd_Vol_Ele = 0
 
      holds for all test functions v'. When the test-function v' is BF_electrode,
      the boundary term reduces to
 
-     (n . (epsilon Grad v), BF_electrode)_Sur_Electrodes_Ele.
+     (n . (epsilon grad v), BF_electrode)_Sur_Electrodes_Ele.
 
      Since BF_electrode == 1 on the electrode, the boundary term is actually
-     simply equal to the integral of (n . epsilon Grad v) on the electrode,
+     simply equal to the integral of (n . epsilon grad v) on the electrode,
      i.e. the flux of the displacement field, which is by definition the
      charge Q_electrode carried by the electrodes.
 
@@ -215,7 +215,7 @@ Formulation {
      satisfies (just consider the equation corresponding to the test function
      BF_electrode):
 
-     Q_electrode = (-epsilon Grad v, Grad BF_electrode)_Vol_Ele */
+     Q_electrode = (-epsilon grad v, grad BF_electrode)_Vol_Ele */
   { Name Electrostatics_v; Type FemEquation;
     Quantity {
       { Name v; Type Local; NameOfSpace Hgrad_v_Ele; }
diff --git a/Thermics/brick.pro b/Thermics/brick.pro
index f2720fb14d151c371eb675f194aa5d3dd2f93f36..78ca76cf29a4533c686afe980b3b222dc8e6ecea 100644
--- a/Thermics/brick.pro
+++ b/Thermics/brick.pro
@@ -15,149 +15,126 @@
        Run (button at the bottom of the left panel)
    ------------------------------------------------------------------- */
 
-/* This model is a rectangular brick with two windows,
-   where various kinds of thermal constraints can be set.
-   Dirichlet, Neumann and convection boundary conditions are imposed
-   on different parts of the surface of the brick.
-   The model is rather academic but it
-   demonstrates some useful high-level GetDP features.
+/* This model is a rectangular brick with two windows, where various kinds of
+   thermal constraints can be set.  Dirichlet, Neumann and convection boundary
+   conditions are imposed on different parts of the surface of the brick.  The
+   model is rather academic but it demonstrates some useful high-level GetDP
+   features.
 
    Governing eqations are
 
-   div ( -lambda[] grad T ) = Q    in Vol_Lambda_The
-   -lambda[] grad T . n = qn = 0   on Sur_Neumann_The
+   div ( -lambda grad T ) = Q    in Vol_The
+   -lambda grad T . n = qn = 0   on Sur_Neu_The
 
    Contact thermal resistance:
-   First, it is shown how to implement contact thermal resistances
-   with surface elements. The surface elements are associated with a
-   thickness and a thermal conductivity (typically much lower than that
-   of surrounding regions). The implementation takes advantage
-   of the powerful FunctionSpace definition in GetDP.
 
-   With the flag "Flag_Regularization", the contact surface
-   can be, for the sake of comparison, replaced
-   by a thin volume conducting region.
+   First, it is shown how to implement contact thermal resistances with surface
+   elements. The surface elements are associated with a thickness and a thermal
+   conductivity (typically much lower than that of surrounding regions). The
+   implementation takes advantage of the powerful FunctionSpace definition in
+   GetDP.
+
+   With the flag "Flag_Regularization", the contact surface can be, for the sake
+   of comparison, replaced by a thin volume conducting region.
 
    Thermal "electrode":
-   The floating potential idea (introduced in tutorial 4)
-   is reconsidered here in a thermal context to represent a region
-   with a very large thermal conductivity where, consequently,
-   the temperature field is uniform (exactly like the electric potential
-   is uniform on an electrode).
-   The dual quantity of this uniform temperature "T_electrode" [K]
-   (which is the "associated global quantity" in GetDP language)
-   is the heat flux "Q_electrode" [W] injected in the electrode
-   by the agent that maintains the temperature equal
+
+   The floating potential idea (introduced in tutorial 4) is reconsidered here
+   in a thermal context to represent a region with a very large thermal
+   conductivity where, consequently, the temperature field is uniform (exactly
+   like the electric potential is uniform on an electrode).  The dual quantity
+   of this uniform temperature "T_electrode" [K] (which is the "associated
+   global quantity" in GetDP language) is the heat flux "Q_electrode" [W]
+   injected in the electrode by the agent that maintains the temperature equal
    to the prescribed value.
 
-   The value of Q_electrode is a by-product of the system resolution
-   provided the term
+   The value of Q_electrode is a by-product of the system resolution provided
+   the term
 
    GlobalTerm { [-Dof{Q_electrode} , {T_electrode} ] ; In Tfloating_The ; }
 
-   is present in the "Resolution" section. This term triggers the writing
-   in the linear system of a supplementary equation associated
-   with the global basis function BF{T_electrode}.
-   All integrations are automatically done by getDP,
-   and the value of Q_electrode is obtained in postprocessing
-   with the PostOperation
+   is present in the "Resolution" section. This term triggers the writing in the
+   linear system of a supplementary equation associated with the global basis
+   function BF{T_electrode}. All integrations are automatically done by GetDP,
+   and the value of Q_electrode is obtained in postprocessing with the
+   PostOperation
 
    Print[ Q_electrode, OnRegion Tfloating_The, ... ]
 
    Heat flux through surfaces:
-   The purpose of a thermal simulation usually goes beyond
-   the mere calculation of a temperature distribution.
-   One is in general also interested in evaluating
-   the heat flux q(S) through some specific surface S:
-
-   q(S) = ( -lambda[] grad T . n )_S
-
-   This quantity cannot be computed from the temperature
-   distribution available on the surface S only.
-   As heat flux is related with the gradient of temperature
-   in the direction normal to the surface, its computation relies on
-   the temperature distribution in a neighborhood of the surface.
-   This means that volume elements in contact with the considered surface
-   need be involved in the computation.
-   To achieve this with getDP, a good method proceeds
-   by the definition a smooth auxiliary function g(S),
-   with g(S)=1 on S, and g(S)=0 outside a finite neighborhood of S.
-   Typically, g(S) is the sum of the shape functions of the nodes on S.
-   Let w(S) be the support of g(S),
-   and let dw(S) denote the boundary of w(S).
-   We then have, just adding and substracting dw(S)
-   to the surface of integration S
-
-   q(S) = ( -lambda[] grad T . n g(S) )_{ dw(S) - ( dw(S)-S ) }.
-
-   dw(S) being a boundary, Stokes theorem can be invoked and,
-   after an integration by part one ends up with
-
-   q(S) = ( -lambda[] grad T . grad g(S) )_w(S)
-        - ( -lambda[] grad T . n g(S) )_{dw(S)-S}.
+
+   The purpose of a thermal simulation usually goes beyond the mere calculation
+   of a temperature distribution.  One is in general also interested in
+   evaluating the heat flux q(S) through some specific surface S:
+
+   q(S) = ( -lambda grad T . n )_S
+
+   This quantity cannot be computed from the temperature distribution available
+   on the surface S only.  As heat flux is related with the gradient of
+   temperature in the direction normal to the surface, its computation relies on
+   the temperature distribution in a neighborhood of the surface.  This means
+   that volume elements in contact with the considered surface need be involved
+   in the computation.  To achieve this with getDP, a good method proceeds by
+   the definition a smooth auxiliary function g(S), with g(S)=1 on S, and g(S)=0
+   outside a finite neighborhood of S.  Typically, g(S) is the sum of the shape
+   functions of the nodes on S.  Let w(S) be the support of g(S), and let dw(S)
+   denote the boundary of w(S).  We then have, just adding and substracting
+   dw(S) to the surface of integration S
+
+   q(S) = ( -lambda grad T . n g(S) )_{ dw(S) - ( dw(S)-S ) }.
+
+   dw(S) being a boundary, Stokes theorem can be invoked and, after an
+   integration by part one ends up with
+
+   q(S) = ( -lambda grad T . grad g(S) )_w(S)
+        - ( -lambda grad T . n g(S) )_{dw(S)-S}.
 	+ ( Q g(S) )_w(S)
 
    Now, g(S) is zero on {dw(S)-S}, except maybe at some surface elements
-   adjacents to dS, but not in dS.
-   The second terme vanishes then if either S is closed,
-   or adjacent to a homogeneous Neumann boundary condition.
-
-   The third term also vanishes, except if a region
-   with a nonzero heatsource Q is in contact with the surface S.
+   adjacents to dS, but not in dS.  The second terme vanishes then if either S
+   is closed, or adjacent to a homogeneous Neumann boundary condition.
 
-   So we have nearly always the following practical formula
-   to evaluate the heat flux across a surface S,
-   in terms of a well-chosen auxiliary scalar function g(S).
+   The third term also vanishes, except if a region with a nonzero heatsource Q
+   is in contact with the surface S.
 
-   q(S) = ( -lambda[] grad T . grad g(S) )_{support of g(S)}
+   So we have nearly always the following practical formula to evaluate the heat
+   flux across a surface S, in terms of a well-chosen auxiliary scalar function
+   g(S).
 
+   q(S) = ( -lambda grad T . grad g(S) )_{support of g(S)}
 
    Particular cases:
 
-   For the heat flux through the boundary of a thermal electrode,
-   one uses g(S) = BF(T_electrode).
-   Note that this heat flux is equal to Q_electrode
-   in the stationary case.
-   This is the case for the heat flux through the boundary of Window2
+   For the heat flux through the boundary of a thermal electrode, one uses g(S)
+   = BF(T_electrode).  Note that this heat flux is equal to Q_electrode in the
+   stationary case.  This is the case for the heat flux through the boundary of
+   Window2
 
-   Auxiliary functions g(S) are also generated
-   for the surfaces named "Surface_i", i=1,2,3 in the model.
-   Note that the flux computed through Surface_3 is incorrect
-   because this surface is not adjacent to surfaces
-   with homogeneous Neumann boundary conditions.
-*/
+   Auxiliary functions g(S) are also generated for the surfaces named
+   "Surface_i", i=1,2,3 in the model.  Note that the flux computed through
+   Surface_3 is incorrect because this surface is not adjacent to surfaces with
+   homogeneous Neumann boundary conditions. */
 
 Include "brick_common.pro";
 
-QWindow1 =
-  DefineNumber[1e3, Name "Window1/Heat source [W]"];
-
-/* The user is given the choice of setting either the global temperature
-   or the global heat flux in Window2.
-   Check how the variable "Flag_ConstraintWin2" is used at different
-   places to alter the model according to that choice
-   and to manage the visibility of related input and output data.*/
-QWindow2 =
-  DefineNumber[1e3, Name "Window2/Heat source [W]",
-	       Visible Flag_ConstraintWin2];
-TWindow2 =
-  DefineNumber[50, Name "Window2/Temperature [degC]",
-	       Visible !Flag_ConstraintWin2];
-outQWindow2 =
-  DefineNumber[0, Name "Output/Q window 2 [degC]",
-	       Visible !Flag_ConstraintWin2, Highlight "Ivory"];
-outTWindow2 =
-  DefineNumber[0, Name "Output/T window 2 [degC]",
-	       Visible Flag_ConstraintWin2, Highlight "Ivory"];
-
-
-ConvectionCoef =
-  DefineNumber[1000, Name"Surface2/hconv",
-	       Label "Convection coefficient [W/(m^2K)]"];
-T_Ambiance =
-  DefineNumber[20, Name"Surface2/Ambiance temperature [degC]"];
-T_Dirichlet =
-  DefineNumber[20, Name"Surface1/Imposed temperature [degC]"];
+QWindow1 = DefineNumber[1e3, Name "Window1/Heat source [W]"];
+
+/* The user is given the choice of setting either the global temperature or the
+   global heat flux in Window2.  Check how the variable "Flag_ConstraintWin2" is
+   used at different places to alter the model according to that choice and to
+   manage the visibility of related input and output data.*/
+QWindow2 = DefineNumber[1e3, Name "Window2/Heat source [W]",
+                        Visible Flag_ConstraintWin2];
+TWindow2 = DefineNumber[50, Name "Window2/Temperature [degC]",
+	                Visible !Flag_ConstraintWin2];
+outQWindow2 = DefineNumber[0, Name "Output/Q window 2 [degC]",
+	                   Visible !Flag_ConstraintWin2, Highlight "Ivory"];
+outTWindow2 = DefineNumber[0, Name "Output/T window 2 [degC]",
+	                   Visible Flag_ConstraintWin2, Highlight "Ivory"];
+ConvectionCoef = DefineNumber[1000, Name"Surface2/hconv",
+	                      Label "Convection coefficient [W/(m^2K)]"];
+T_Ambiance = DefineNumber[20, Name"Surface2/Ambiance temperature [degC]"];
+T_Dirichlet = DefineNumber[20, Name"Surface1/Imposed temperature [degC]"];
 
 Group {
   /* Geometrical regions: */
@@ -175,17 +152,17 @@ Group {
 
   /* Abstract regions:
 
-     Vol_Lambda_Ele     : volume regions with a thermal conductivity lambda[]
-     Sur_Dirichlet_The  : Dirichlet bondary condition surface
-     Sur_Neu_Ele        : Neumann bondary condition surface
+     Vol_The            : volume regions with a thermal conductivity
+     Sur_Dir_The        : non-homogeneous Dirichlet boundary condition surface
+     Sur_Neu_The        : homogeneous Neumann boundary condition surface
      Sur_Convection_The : convective surface q.n = h ( T-Tinf )
      Vol_Qsource_The    : volume heat source regions
      Tfloating_The      : thermal electrodes
   */
 
-  Vol_Lambda_The = Region[ {Brick, Window1, Window2, LayerWindow1} ];
-  Sur_Dirichlet_The = Region[ Surface_1 ];
-  Sur_Neumann_The = Region[ Surface_3 ];
+  Vol_The = Region[ {Brick, Window1, Window2, LayerWindow1} ];
+  Sur_Dir_The = Region[ Surface_1 ];
+  Sur_Neu_The = Region[ Surface_3 ];
   Sur_Convection_The = Region[ Surface_2 ];
 
   Vol_Qsource_The = Region[ Window1 ];
@@ -200,7 +177,6 @@ Group {
   EndIf
 }
 
-
 Function {
   lambda_Brick = 50.; // steel
   lambda[ Brick ] = lambda_Brick;
@@ -226,7 +202,7 @@ Function {
 Constraint {
   { Name Dirichlet_The ;
     Case {
-      { Region Sur_Dirichlet_The ; Value T_Dirichlet ; }
+      { Region Sur_Dir_The ; Value T_Dirichlet ; }
     }
   }
   { Name T_Discontinuity ;
@@ -260,8 +236,7 @@ Constraint {
 Integration {
   { Name Int ;
     Case {
-      {
-        Type Gauss ;
+      { Type Gauss ;
         Case {
           { GeoElement Triangle   ; NumberOfPoints  4 ; }
           { GeoElement Quadrangle ; NumberOfPoints  4 ; }
@@ -286,12 +261,9 @@ Jacobian {
 }
 
 Group {
-  Dom_Hgrad_T = Region[ {Vol_Lambda_The,
-			 Sur_Convection_The,
-			 Sur_Tdisc_The} ];
+  Dom_Hgrad_T = Region[ {Vol_The, Sur_Convection_The, Sur_Tdisc_The} ];
   DomainWithSurf_TL_The =
-    ElementsOf[ {Vol_OneSide_The, Sur_Tdisc_The},
-		OnOneSideOf Sur_Tdisc_The ];
+    ElementsOf[ {Vol_OneSide_The, Sur_Tdisc_The}, OnOneSideOf Sur_Tdisc_The ];
 }
 
 FunctionSpace {
@@ -336,7 +308,6 @@ FunctionSpace {
     }
   }
   EndFor
-
 }
 
 Formulation {
@@ -350,36 +321,33 @@ Formulation {
       For i In {1:NbSurface}
         { Name un~{i} ; Type Local ; NameOfSpace FluxLayer~{i} ; }
       EndFor
-
     }
     Equation {
-
       Integral { [ lambda[]  * Dof{d T} , {d T} ] ;
-                 In Vol_Lambda_The; Integration Int ; Jacobian Vol ; }
+        In Vol_The; Integration Int ; Jacobian Vol ; }
 
       Integral { [ ( lambda[]/thickness[] ) * Dof{Tdisc} , {Tdisc} ] ;
-	         In Sur_Tdisc_The; Integration Int ; Jacobian Sur ; }
+        In Sur_Tdisc_The; Integration Int ; Jacobian Sur ; }
 
       Integral { [ -Qsource[] , {T} ] ;
-                 In Vol_Qsource_The ; Integration Int ; Jacobian Vol ; }
+        In Vol_Qsource_The ; Integration Int ; Jacobian Vol ; }
 
       Integral { [ h[] * Dof{T} , {T} ] ;
-                 In Sur_Convection_The ; Integration Int ; Jacobian Sur ; }
+        In Sur_Convection_The ; Integration Int ; Jacobian Sur ; }
 
       Integral { [ -h[] * Tinf[] , {T} ] ;
-                 In Sur_Convection_The ; Integration Int ; Jacobian Sur ; }
+        In Sur_Convection_The ; Integration Int ; Jacobian Sur ; }
 
       GlobalTerm { [-Dof{Qglob} , {Tglob} ] ; In Tfloating_The ; }
 
       For i In {1:NbSurface}
         Integral { [ 0 * Dof{un~{i}} , {un~{i}} ] ;
-          In Vol_Lambda_The ; Integration Int ; Jacobian Vol ; }
+          In Vol_The ; Integration Int ; Jacobian Vol ; }
       EndFor
     }
   }
 }
 
-
 Resolution {
   { Name Thermal_T ;
     System {
@@ -411,9 +379,9 @@ PostProcessing {
 
       For i In {1:NbSurface}
       { Name un~{i} ; Value { Local { [ {un~{i}} ] ;
-	    In Vol_Lambda_The ; Jacobian Vol ; } } }
+	    In Vol_The ; Jacobian Vol ; } } }
       { Name IFlux~{i} ; Value { Integral { [ -lambda[]*{d T} * {d un~{i}} ];
-	    In Vol_Lambda_The ; Jacobian Vol ; Integration Int ; } } }
+	    In Vol_The ; Jacobian Vol ; Integration Int ; } } }
       EndFor
 
     }
@@ -422,8 +390,8 @@ PostProcessing {
 
 PostOperation Map_T UsingPost Thermal_T {
   If( !Flag_Regularization )
-    Print[ Tcont, OnElementsOf Vol_Lambda_The, File "Tcont.pos"] ;
-    Print[ Tdisc, OnElementsOf Vol_Lambda_The, File "Tdisc.pos"] ;
+    Print[ Tcont, OnElementsOf Vol_The, File "Tcont.pos"] ;
+    Print[ Tdisc, OnElementsOf Vol_The, File "Tdisc.pos"] ;
   EndIf
 
   If(Flag_ConstraintWin2)
@@ -435,13 +403,13 @@ PostOperation Map_T UsingPost Thermal_T {
   EndIf
 
   For i In {1:NbSurface}
-    Print[un~{i}, OnElementsOf Vol_Lambda_The,
+    Print[un~{i}, OnElementsOf Vol_The,
     	  File Sprintf("FluxLayer_%g.pos",i)];
-    Print[ IFlux~{i}[Vol_Lambda_The], OnGlobal,
+    Print[ IFlux~{i}[Vol_The], OnGlobal,
 	   Format TimeTable, File > "Fluxes.dat", Color "Ivory",
 	   SendToServer Sprintf("Output/Heat flux surface %g [W]", i)];
   EndFor
-  Print[ T, OnElementsOf Vol_Lambda_The, File "T.pos" ] ;
+  Print[ T, OnElementsOf Vol_The, File "T.pos" ] ;
   Echo[ StrCat["l=PostProcessing.NbViews-1;",
 	       "View[l].IntervalsType = 3;",
 	       "View[l].NbIso = 30;",