diff --git a/doc/texinfo/gmsh.texi b/doc/texinfo/gmsh.texi
index 39ade95af5ea14c02f18638c5491088c31aec6bf..ffd65316186ae583771abc45a59d359e88d695d6 100644
--- a/doc/texinfo/gmsh.texi
+++ b/doc/texinfo/gmsh.texi
@@ -5126,8 +5126,8 @@ system}. Screencasts that show how to use the GUI are available on
 * t1.geo:: Geometry basics, elementary entities, physical groups
 * t2.geo:: Transformations, extruded geometries, volumes
 * t3.geo:: Extruded meshes, input parameters, options
-* t4.geo:: Built-in functions, surface holes, annotations, mesh colors
-* t5.geo:: Characteristic lengths, arrays of variables, macros, loops
+* t4.geo:: Built-in functions, holes in surfaces, annotations, entity colors
+* t5.geo:: Characteristic lengths, macros, loops, holes in volumes
 * t6.geo:: Transfinite meshes
 * t7.geo:: Background meshes
 * t8.geo:: Post-processing, scripting, animations, options
diff --git a/tutorial/README.txt b/tutorial/README.txt
index 9a4af38ef15face086a8b462893d6c260215866e..0444264564cf23b142f2d6184f2c55772ea2512a 100644
--- a/tutorial/README.txt
+++ b/tutorial/README.txt
@@ -12,11 +12,11 @@ the command line, run "gmsh t1.geo" (which will launch the GUI) or "gmsh t1.geo
    instructions on how to compile the app from source.
 
 The `c++', `c', `python' and `julia' subdirectories contain C++, C, Python and
-Julia versions of (some of these) tutorials, written using the Gmsh Application
+Julia versions of (some of) these tutorials, written using the Gmsh Application
 Programming Interface (API): you will need the Gmsh dynamic library and the
 associated header files (for C++ and C) or modules (for Python and Julia) to run
-them. Each subdirectory contains additional information on how to run the
-tutorials for each supported language, as well as a extended tutorials (starting
+them. Each subdirectory also contains additional information on how to run the
+tutorials for each supported language, as well as extended tutorials (starting
 with `x') introducing features available through the API but not available in
 `.geo' files.
 
diff --git a/tutorial/c++/README.txt b/tutorial/c++/README.txt
index d0a158c6900be0d3cbb3798f3c234125bdcd7b5f..7c72104f639b9e3a3c7a48a2a702d81ce19ef487 100644
--- a/tutorial/c++/README.txt
+++ b/tutorial/c++/README.txt
@@ -1,5 +1,4 @@
-This directory contains C++ versions of the tutorials, written using the Gmsh
-API.
+This directory contains the Gmsh C++ tutorials, written using the Gmsh C++ API.
 
 To compile and run the C++ tutorials, you need the Gmsh dynamic library and the
 associated header file (`gmsh.h'). These can be either obtained
diff --git a/tutorial/c++/t1.cpp b/tutorial/c++/t1.cpp
index 987cf2a10544e2a25a3a2741c9a44f686b4c439f..823c3b3232710b9c06ae0780b003fe4af9aad340 100644
--- a/tutorial/c++/t1.cpp
+++ b/tutorial/c++/t1.cpp
@@ -1,74 +1,121 @@
-// This file reimplements gmsh/tutorial/t1.geo in C++.
-
-// For all the elementary explanations about the general philosphy of entities
-// in Gmsh, see the comments in the .geo file. Comments here focus on the
-// specifics of the C++ API.
-
-// The Gmsh API is entirely defined in the <gmsh.h> header. Read this file: it
-// contains the documentation for all the functions in the API.
+/*******************************************************************************
+ *
+ *  Gmsh C++ tutorial 1
+ *
+ *  Elementary entities (points, curves, surfaces), physical groups (points,
+ *  curves, surfaces)
+ *
+ *******************************************************************************/
+
+// The Gmsh C++ API is entirely defined in the `gmsh.h' header (which contains
+// the full documentation of all the functions in the API):
 #include <gmsh.h>
 
 int main(int argc, char **argv)
 {
-  // Before using any functions in the C++ API, Gmsh must be initialized.
+  // Before using any functions in the C++ API, Gmsh must be initialized:
   gmsh::initialize();
 
   // By default Gmsh will not print out any messages: in order to output
-  // messages on the terminal, just set the standard Gmsh option
-  // "General.Terminal" (same format and meaning as in .geo files) using
-  // gmsh::option::setNumber():
+  // messages on the terminal, just set the "General.Terminal" option to 1:
   gmsh::option::setNumber("General.Terminal", 1);
 
-  // This adds a new model, named "t1". If gmsh::model::add() is not called, a
+  // We now add a new model, named "t1". If gmsh::model::add() is not called, a
   // new default (unnamed) model will be created on the fly, if necessary.
   gmsh::model::add("t1");
 
-  // The C++ API provides direct access to the internal CAD kernels. The
-  // built-in CAD kernel was used in t1.geo: the corresponding API functions
-  // live in the "gmsh::model::geo" namespace. To create geometrical points with
-  // the built-in CAD kernel, one thus uses gmsh::model::geo::addPoint():
-  //
+  // The C++ API provides direct access to each supported geometry kernel. The
+  // built-in kernel is used in this first tutorial: the corresponding API
+  // functions live in the "gmsh::model::geo" namespace.
+
+  // The first type of `elementary entity' in Gmsh is a `Point'. To create a
+  // point with the built-in geometry kernel, the C++ API function is
+  // gmsh::model::geo::addPoint():
   // - the first 3 arguments are the point coordinates (x, y, z)
-  //
-  // - the next (optional) argument is the target mesh size close to the point
-  //
-  // - the last (optional) argument is the point tag
+  // - the next (optional) argument is the target mesh size (the "characteristic
+  //   length") close to the point
+  // - the last (optional) argument is the point tag (a stricly positive integer
+  //   that uniquely identifies the point)
   double lc = 1e-2;
   gmsh::model::geo::addPoint(0, 0, 0, lc, 1);
+
+  // The distribution of the mesh element sizes will be obtained by
+  // interpolation of these characteristic lengths throughout the
+  // geometry. Another method to specify characteristic lengths is to use
+  // general mesh size Fields (see `t10.cpp'). A particular case is the use of a
+  // background mesh (see `t7.cpp').
+  //
+  // If no target mesh size of provided, a default uniform coarse size will be
+  // used for the model, based on the overall model size.
+  //
+  // We can then define some additional points. All points should have different
+  // tags:
   gmsh::model::geo::addPoint(.1, 0,  0, lc, 2);
   gmsh::model::geo::addPoint(.1, .3, 0, lc, 3);
-  gmsh::model::geo::addPoint(0,  .3, 0, lc, 4);
 
-  // The API to create lines with the built-in kernel follows the same
-  // conventions: the first 2 arguments are point tags, the last (optional one)
-  // is the line tag.
+  // If the tag is not provided explicitly, a new tag is automatically created,
+  // and returned by the function:
+  int p4 = gmsh::model::geo::addPoint(0, .3, 0, lc);
+
+  // Curves are Gmsh's second type of elementery entities, and, amongst curves,
+  // straight lines are the simplest. The API to create straight line segments
+  // with the built-in kernel follows the same conventions: the first 2
+  // arguments are point tags (the start and end points of the line), and the
+  // last (optional one) is the line tag.
+  //
+  // In the commands below, for example, the line 1 starts at point 1 and ends
+  // at point 2.
+  //
+  // Note that curve tags are separate from point tags - hence we can reuse tag
+  // `1' for our first curve. And as a general rule, elementary entity tags in
+  // Gmsh have to be unique per geometrical dimension.
   gmsh::model::geo::addLine(1, 2, 1);
   gmsh::model::geo::addLine(3, 2, 2);
-  gmsh::model::geo::addLine(3, 4, 3);
-  gmsh::model::geo::addLine(4, 1, 4);
-
-  // The philosophy to construct curve loops and surfaces is similar: the first
-  // argument is now a vector of integers.
+  gmsh::model::geo::addLine(3, p4, 3);
+  gmsh::model::geo::addLine(4, 1, p4);
+
+  // The third elementary entity is the surface. In order to define a simple
+  // rectangular surface from the four curves defined above, a curve loop has
+  // first to be defined. A curve loop is defined by an ordered list of
+  // connected curves, a sign being associated with each curve (depending on the
+  // orientation of the curve to form a loop). The API function to create curve
+  // loops takes a vector of integers as first argument, and the curve loop tag
+  // (which must ne unique amongst curve loops) as the second (optional)
+  // argument:
   gmsh::model::geo::addCurveLoop({4, 1, -2, 3}, 1);
-  gmsh::model::geo::addPlaneSurface({1}, 1);
-
-  // Physical groups are defined by providing the dimension of the group (0 for
-  // physical points, 1 for physical curves, 2 for physical surfaces and 3 for
-  // phsyical volumes) followed by a vector of entity tags. The last (optional)
-  // argument is the tag of the new group to create.
-  gmsh::model::addPhysicalGroup(0, {1, 2}, 1);
-  gmsh::model::addPhysicalGroup(1, {1, 2}, 2);
-  gmsh::model::addPhysicalGroup(2, {1}, 6);
 
-  // Physical names are also defined by providing the dimension and tag of the
-  // entity.
-  gmsh::model::setPhysicalName(2, 6, "My surface");
+  // We can then define the surface as a list of curve loops (only one here,
+  // representing the external contour, since there are no holes--see `t4.cpp'
+  // for an example of a surface with a hole):
+  gmsh::model::geo::addPlaneSurface({1}, 1);
 
-  // Before it can be meshed, the internal CAD representation must be
+  // At this level, Gmsh knows everything to display the rectangular surface 6
+  // and to mesh it. An optional step is needed if we want to group elementary
+  // geometrical entities into more meaningful groups, e.g. to define some
+  // mathematical ("domain", "boundary"), functional ("left wing", "fuselage")
+  // or material ("steel", "carbon") properties.
+  //
+  // Such groups are called "Physical Groups" in Gmsh. By default, if physical
+  // groups are defined, Gmsh will export in output files only mesh elements
+  // that belong to at least one physical group. (To force Gmsh to save all
+  // elements, whether they belong to physical groups or not, set the
+  // `Mesh.SaveAll' option to 1.) Physical groups are also identified by tags,
+  // i.e. stricly positive integers, that should be unique per dimension (0D,
+  // 1D, 2D or 3D). Physical groups can also be given names.
+  //
+  // Here we define a physical curve that groups the left, bottom and right
+  // curves in a single group (with prescribed tag 5); and a physical surface
+  // with name "My surface" (with an automatic tag) containing the geometrical
+  // surface 1:
+  gmsh::model::addPhysicalGroup(1, {1, 2, 4}, 5);
+  int ps = gmsh::model::addPhysicalGroup(2, {1});
+  gmsh::model::setPhysicalName(2, ps, "My surface");
+
+  // Before it can be meshed, the internal geometrical representation must be
   // synchronized with the Gmsh model, which will create the relevant Gmsh data
   // structures. This is achieved by the gmsh::model::geo::synchronize() API
-  // call for the built-in CAD kernel. Synchronizations can be called at any
-  // time, but they involve a non trivial amount of processing; so while you
+  // call for the built-in geometry kernel. Synchronizations can be called at
+  // any time, but they involve a non trivial amount of processing; so while you
   // could synchronize the internal CAD data after every CAD command, it is
   // usually better to minimize the number of synchronization points.
   gmsh::model::geo::synchronize();
@@ -85,10 +132,32 @@ int main(int argc, char **argv)
   //
   // gmsh::option::setNumber("Mesh.SaveAll", 1);
 
-  // We could run the graphical user interface with
+  // We could run the graphical user interface with:
   // gmsh::fltk::run();
 
-  // This should be called at the end:
+  // Note that starting with Gmsh 3.0, models can be built using other geometry
+  // kernels than the default "built-in" kernel. To use the OpenCASCADE geometry
+  // kernel instead of the built-in kernel, you should use the functions in the
+  // "gmsh::model::occ" namespace.
+  //
+  // Different geometry kernels have different features. With OpenCASCADE,
+  // instead of defining the surface by successively defining 4 points, 4 curves
+  // and 1 curve loop, one can define the rectangular surface directly with
+  //
+  // gmsh::model::occ::addRectangle(.2, 0, 0, .1, .3);
+  //
+  // Afther synchronization with the Gmsh model with
+  //
+  // gmsh::model::occ::synchronize();
+  //
+  // the underlying curves and points could be accessed with
+  // gmsh::model::getBoundary().
+  //
+  // See e.g `t16.cpp', `t18.cpp' or `t19.cpp' for complete examples based on
+  // OpenCASCADE, and `demos/api' for more.
+
+  // This should be called when you are done using the Gmsh C++ API:
   gmsh::finalize();
+
   return 0;
 }
diff --git a/tutorial/c++/t10.cpp b/tutorial/c++/t10.cpp
index 3bb292ba51faa370b92522d6d06b516ca506421d..bebefee5c5326694aec344530266e5a8b52ea525 100644
--- a/tutorial/c++/t10.cpp
+++ b/tutorial/c++/t10.cpp
@@ -1,4 +1,14 @@
-// This reimplements gmsh/tutorial/t10.geo in C++.
+/*******************************************************************************
+ *
+ *  Gmsh C++ tutorial 10
+ *
+ *  General mesh size fields
+ *
+ *******************************************************************************/
+
+// In addition to specifying target mesh sizes at the points of the geometry
+// (see `t1.cpp') or using a background mesh (see `t7.cppq'), you can use
+// general mesh size "Fields".
 
 #include <gmsh.h>
 #include <sstream>
@@ -11,8 +21,9 @@ int main(int argc, char **argv)
   gmsh::initialize(argc, argv);
   gmsh::option::setNumber("General.Terminal", 1);
 
-  model::add("t1");
+  model::add("t10");
 
+  // Let's create a simple rectangular geometry
   double lc = .15;
   factory::addPoint(0.0,0.0,0,lc, 1);
   factory::addPoint(1,0.0,0,lc, 2);
@@ -28,11 +39,29 @@ int main(int argc, char **argv)
   factory::addCurveLoop({1,2,3,4}, 5);
   factory::addPlaneSurface({5}, 6);
 
+  factory::synchronize();
+
+  // Say we would like to obtain mesh elements with size lc/30 near curve 2 and
+  // point 5, and size lc elsewhere. To achieve this, we can use two fields:
+  // "Distance", and "Threshold". We first define a Distance field (`Field[1]')
+  // on points 5 and on curve 2. This field returns the distance to point 5 and
+  // to (100 equidistant points on) curve 2.
   model::mesh::field::add("Distance", 1);
   model::mesh::field::setNumbers(1, "NodesList", {5});
   model::mesh::field::setNumber(1, "NNodesByEdge", 100);
   model::mesh::field::setNumbers(1, "EdgesList", {2});
 
+  // We then define a `Threshold' field, which uses the return value of the
+  // `Distance' field 1 in order to define a simple change in element size
+  // depending on the computed distances
+  //
+  // LcMax -                         /------------------
+  //                               /
+  //                             /
+  //                           /
+  // LcMin -o----------------/
+  //        |                |       |
+  //      Point           DistMin DistMax
   model::mesh::field::add("Threshold", 2);
   model::mesh::field::setNumber(2, "IField", 1);
   model::mesh::field::setNumber(2, "LcMin", lc / 30);
@@ -40,17 +69,27 @@ int main(int argc, char **argv)
   model::mesh::field::setNumber(2, "DistMin", 0.15);
   model::mesh::field::setNumber(2, "DistMax", 0.5);
 
+  // Say we want to modulate the mesh element sizes using a mathematical
+  // function of the spatial coordinates. We can do this with the MathEval
+  // field:
   model::mesh::field::add("MathEval", 3);
   model::mesh::field::setString(3, "F", "Cos(4*3.14*x) * Sin(4*3.14*y) / 10 + 0.101");
 
+  // We could also combine MathEval with values coming from other fields. For
+  // example, let's define a `Distance' field around point 1
   model::mesh::field::add("Distance", 4);
   model::mesh::field::setNumbers(4, "NodesList", {1});
 
+  // We can then create a `MathEval' field with a function that depends on the
+  // return value of the `Distance' field 4, i.e., depending on the distance to
+  // point 1 (here using a cubic law, with minimum element size = lc / 100)
   model::mesh::field::add("MathEval", 5);
   std::stringstream stream;
   stream << "F4^3 + " << lc / 100;
   model::mesh::field::setString(5, "F", stream.str());
 
+  // We could also use a `Box' field to impose a step change in element sizes
+  // inside a box
   model::mesh::field::add("Box", 6);
   model::mesh::field::setNumber(6, "VIn", lc / 15);
   model::mesh::field::setNumber(6, "VOut", lc);
@@ -59,12 +98,45 @@ int main(int argc, char **argv)
   model::mesh::field::setNumber(6, "YMin", 0.3);
   model::mesh::field::setNumber(6, "YMax", 0.6);
 
+  // Many other types of fields are available: see the reference manual for a
+  // complete list. You can also create fields directly in the graphical user
+  // interface by selecting `Define->Fields' in the `Mesh' module.
+
+  // Finally, let's use the minimum of all the fields as the background mesh
+  // field
   model::mesh::field::add("Min", 7);
   model::mesh::field::setNumbers(7, "FieldsList", {2, 3, 5, 6});
 
   model::mesh::field::setAsBackgroundMesh(7);
 
-  factory::synchronize();
+  // To determine the size of mesh elements, Gmsh locally computes the minimum
+  // of
+  //
+  // 1) the size of the model bounding box;
+  // 2) if `Mesh.CharacteristicLengthFromPoints' is set, the mesh size specified
+  //    at geometrical points;
+  // 3) if `Mesh.CharacteristicLengthFromCurvature' is set, the mesh size based
+  //    on the curvature and `Mesh.MinimumCirclePoints';
+  // 4) the background mesh field;
+  // 5) any per-entity mesh size constraint.
+  //
+  // This value is then constrained in the interval
+  // [`Mesh.CharacteristicLengthMin', `MeshCharacteristicLengthMax'] and
+  // multiplied by `Mesh.CharacteristicLengthFactor'.  In addition, boundary
+  // mesh sizes (on curves or surfaces) are interpolated inside the enclosed
+  // entity (surface or volume, respectively) if the option
+  // `Mesh.CharacteristicLengthExtendFromBoundary' is set (which is the case by
+  // default).
+  //
+  // When the element size is fully specified by a background mesh (as it is in
+  // this example), it is thus often desirable to set
+
+  gmsh::option::setNumber("Mesh.CharacteristicLengthExtendFromBoundary", 0);
+  gmsh::option::setNumber("Mesh.CharacteristicLengthFromPoints", 0);
+  gmsh::option::setNumber("Mesh.CharacteristicLengthFromCurvature", 0);
+
+  // This will prevent over-refinement due to small mesh sizes on the boundary.
+
   model::mesh::generate(2);
   gmsh::write("t10.msh");
   gmsh::finalize();
diff --git a/tutorial/c++/t11.cpp b/tutorial/c++/t11.cpp
index 40ed8eac6bf9d0569243d1731d4d96dfdc8559ab..18889509e5e0f63dd3013c6695203eac7c61c640 100644
--- a/tutorial/c++/t11.cpp
+++ b/tutorial/c++/t11.cpp
@@ -1,4 +1,10 @@
-// This file reimplements gmsh/tutorial/t11.geo in C++.
+/*******************************************************************************
+ *
+ *  Gmsh C++ tutorial 11
+ *
+ *  Unstructured quadrangular meshes
+ *
+ *******************************************************************************/
 
 #include <gmsh.h>
 
@@ -9,6 +15,11 @@ int main(int argc, char **argv)
 
   gmsh::model::add("t11");
 
+  // We have seen in tutorials `t3.cpp' and `t6.cpp' that extruded and
+  // transfinite meshes can be "recombined" into quads/prisms/hexahedra by using
+  // the "Recombine" keyword. Unstructured meshes can be recombined in the same
+  // way. Let's define a simple geometry with an analytical mesh size field:
+
   int p1 = gmsh::model::geo::addPoint(-1.25, -.5, 0);
   int p2 = gmsh::model::geo::addPoint(1.25, -.5, 0);
   int p3 = gmsh::model::geo::addPoint(1.25, 1.25, 0);
@@ -24,34 +35,61 @@ int main(int argc, char **argv)
 
   gmsh::model::geo::synchronize();
 
-  // add an analytical size field with tag 1
   gmsh::model::mesh::field::add("MathEval", 1);
   gmsh::model::mesh::field::setString
     (1, "F", "0.01*(1.0+30.*(y-x*x)*(y-x*x) + (1-x)*(1-x))");
   gmsh::model::mesh::field::setAsBackgroundMesh(1);
 
-  // to generate quadrangles instead of triangles for the plane surface, add the
-  // recombine constraint before meshing
+  // To generate quadrangles instead of triangles, we can simply add
   gmsh::model::mesh::setRecombine(2, pl);
-  gmsh::model::mesh::generate(2);
 
-  // you could also set the option "Mesh.RecombineAll"
+  // If we'd had several surfaces, we could have used the global option
+  // "Mesh.RecombineAll":
+  //
   // gmsh::option::setNumber("Mesh.RecombineAll", 1);
-  // gmsh::model::mesh::generate(2);
 
-  // You could also apply the recombination algo after meshing
-  // gmsh::model::mesh::generate(2);
-  // gmsh::model::mesh::recombine();
+  // The default recombination algorithm is called "Blossom": it uses a minimum
+  // cost perfect matching algorithm to generate fully quadrilateral meshes from
+  // triangulations. More details about the algorithm can be found in the
+  // following paper: J.-F. Remacle, J. Lambrechts, B. Seny, E. Marchandise,
+  // A. Johnen and C. Geuzaine, "Blossom-Quad: a non-uniform quadrilateral mesh
+  // generator using a minimum cost perfect matching algorithm", International
+  // Journal for Numerical Methods in Engineering 89, pp. 1102-1119, 2012.
 
-  // Better 2D planar quadrilateral meshes with the Frontal-Delaunay for quads
-  // algorithm:
+  // For even better 2D (planar) quadrilateral meshes, you can try the
+  // experimental "Frontal-Delaunay for quads" meshing algorithm, which is a
+  // triangulation algorithm that enables to create right triangles almost
+  // everywhere: J.-F. Remacle, F. Henrotte, T. Carrier-Baudouin, E. Bechet,
+  // E. Marchandise, C. Geuzaine and T. Mouton. A frontal Delaunay quad mesh
+  // generator using the L^inf norm. International Journal for Numerical Methods
+  // in Engineering, 94, pp. 494-512, 2013. Uncomment the following line to try
+  // the Frontal-Delaunay algorithms for quads:
+  //
   // gmsh::option::setNumber("Mesh.Algorithm", 8);
-  // gmsh::model::mesh::generate(2);
 
-  // Force a full-quad mesh with either option
-  // gmsh::option::setNumber("Mesh.SubdivisionAlgorithm", 1);
+  // The default recombination algorithm might leave some triangles in the mesh,
+  // if recombining all the triangles leads to badly shaped quads. In such cases,
+  // to generate full-quad meshes, you can either subdivide the resulting hybrid
+  // mesh (with Mesh.SubdivisionAlgorithm = 1), or use the full-quad recombination
+  // algorithm, which will automatically perform a coarser mesh followed by
+  // recombination, smoothing and subdivision. Uncomment the following line to try
+  // the full-quad algorithm:
+  //
   // gmsh::option::setNumber("Mesh.RecombinationAlgorithm", 2); // or 3
+
+  // You can also set the subdivision step alone, with
+  //
+  // gmsh::option::setNumber("Mesh.SubdivisionAlgorithm", 1);
+
+  gmsh::model::mesh::generate(2);
+
+  // Note that you could also apply the recombination algorithm and/or the
+  // subdivision step explicitly after meshing, as follows:
+  //
   // gmsh::model::mesh::generate(2);
+  // gmsh::model::mesh::recombine();
+  // gmsh::option::setNumber("Mesh.SubdivisionAlgorithm", 1);
+  // gmsh::model::mesh::refine();
 
   gmsh::finalize();
   return 0;
diff --git a/tutorial/c++/t12.cpp b/tutorial/c++/t12.cpp
index a9a73f5cd835f5d71c57a8ed3bc84186a3a2e211..9ce5596deb19f2e9f83de1ea224e8eb481f776fc 100644
--- a/tutorial/c++/t12.cpp
+++ b/tutorial/c++/t12.cpp
@@ -1,4 +1,39 @@
-// This file reimplements gmsh/tutorial/t12.geo in C++.
+/*******************************************************************************
+ *
+ *  Gmsh C++ tutorial 12
+ *
+ *  Cross-patch meshing with compounds
+ *
+ *******************************************************************************/
+
+// "Compound" meshing constraints allow to generate meshes across surface
+// boundaries, which can be useful e.g. for imported CAD models (e.g. STEP) with
+// undesired small features.
+
+// When a `setCompound()' meshing constraint is given, at mesh generation time
+// Gmsh
+//  1. meshes the underlying elementary geometrical entities, individually
+//  2. creates a discrete entity that combines all the individual meshes
+//  3. computes a discrete parametrization (i.e. a piece-wise linear mapping)
+//     on this discrete entity
+//  4. meshes the discrete entity using this discrete parametrization instead
+//     of the underlying geometrical description of the underlying elementary
+//     entities making up the compound
+//  5. optionally, reclassifies the mesh elements and nodes on the original
+//     entities
+
+// Step 3. above can only be performed if the mesh resulting from the
+// combination of the individual meshes can be reparametrized, i.e. if the shape
+// is "simple enough". If the shape is not amenable to reparametrization, you
+// should create a full mesh of the geometry and first re-classify it to
+// generate patches amenable to reparametrization (see `t13.cpp').
+
+// The mesh of the individual entities performed in Step 1. should usually be
+// finer than the desired final mesh; this can be controlled with the
+// `Mesh.CompoundCharacteristicLengthFactor' option.
+
+// The optional reclassification on the underlying elementary entities in Step
+// 5. is governed by the `Mesh.CompoundClassify' option.
 
 #include <gmsh.h>
 
@@ -45,12 +80,15 @@ int main(int argc, char **argv)
 
   factory::synchronize();
 
-  // set compound curve (dim 1); of curves 2, 3, 4
+  // Treat curves 2, 3 and 4 as a single curve when meshing (i.e. mesh across
+  // points 6 and 7)
   gmsh::model::mesh::setCompound(1, {2, 3, 4});
-  // set compound curve of curves 6, 7, 8
+
+  // Idem with curves 6, 7 and 8
   gmsh::model::mesh::setCompound(1, {6, 7, 8});
 
-  // set compound surface from surfaces 1, 5, 10
+  // Treat surfaces 1, 5 and 10 as a single surface when meshing (i.e. mesh
+  // across curves 9 and 10)
   gmsh::model::mesh::setCompound(2, {1, 5, 10});
 
   gmsh::finalize();
diff --git a/tutorial/c++/t13.cpp b/tutorial/c++/t13.cpp
index 8d25e5937e0b5015b72d56363d2d9d738b23cd5a..2a9d78623dbfd3da81fced5094268581e8ed8d33 100644
--- a/tutorial/c++/t13.cpp
+++ b/tutorial/c++/t13.cpp
@@ -1,4 +1,10 @@
-// This file reimplements gmsh/tutorial/t13.geo in C++.
+/*******************************************************************************
+ *
+ *  Gmsh C++ tutorial 13
+ *
+ *  Remeshing without an underlying CAD model
+ *
+ *******************************************************************************/
 
 #include <gmsh.h>
 #include <math.h>
@@ -26,14 +32,23 @@ int main(int argc, char **argv)
   // We first classify ("color") the surfaces by splitting the original surface
   // along sharp geometrical features. This will create new discrete surfaces,
   // curves and points.
-
   double angle = 40; // Angle for surface detection
-  bool forceParametrizablePatches = false; // Create surfaces guaranteed to be parametrizable?
+
+  // For complex geometries, patches can be too complex, too elongated or too
+  // large to be parametrized; setting the following option will force the
+  // creation of patches that are amenable to reparametrization:
+  bool forceParametrizablePatches = false;
+
+  // For open surfaces include the boundary edges in the classification process.
   bool includeBoundary = true;
 
+  // Force curves to be split on given angle:
+  double curveAngle = 180;
+
   gmsh::model::mesh::classifySurfaces(angle * M_PI / 180.,
                                       includeBoundary,
-                                      forceParametrizablePatches);
+                                      forceParametrizablePatches,
+                                      curveAngle);
 
   // Create a geometry for all the discrete curves and surfaces in the mesh, by
   // computing a parametrization for each one
@@ -42,16 +57,14 @@ int main(int argc, char **argv)
   // Create a volume from all the surfaces
   std::vector<std::pair<int, int> > s;
   gmsh::model::getEntities(s, 2);
-
   std::vector<int> sl;
   for(std::size_t i = 0; i < s.size(); i++) sl.push_back(s[i].second);
   int l = gmsh::model::geo::addSurfaceLoop(sl);
-
   gmsh::model::geo::addVolume({l});
 
   gmsh::model::geo::synchronize();
 
-  // element size imposed by a size field, just because we can :-)
+  // We specify element sizes imposed by a size field, just because we can :-)
   bool funny = true;//false;
   int f = gmsh::model::mesh::field::add("MathEval");
   if(funny)
diff --git a/tutorial/c++/t14.cpp b/tutorial/c++/t14.cpp
index b4d676b30bcc5faa2d8b482a5db1bd0fba2b4204..c051fac582bd3ddcd44880a1b52496b40cd42749 100644
--- a/tutorial/c++/t14.cpp
+++ b/tutorial/c++/t14.cpp
@@ -1,4 +1,10 @@
-// This file reimplements gmsh/tutorial/t14.geo in Python.
+/*******************************************************************************
+ *
+ *  Gmsh C++ tutorial 14
+ *
+ *  Homology and cohomology computation
+ *
+ *******************************************************************************/
 
 // Homology computation in Gmsh finds representative chains of (relative)
 // (co)homology space bases using a mesh of a model.  The representative basis
@@ -17,6 +23,8 @@ int main(int argc, char **argv)
 
   gmsh::model::add("t14");
 
+  // Create an example geometry
+
   double m = 0.5; // mesh characteristic length
   double h = 2; // geometry height in the z-direction
 
@@ -128,6 +136,10 @@ int main(int argc, char **argv)
   // Generate the mesh and perform the requested homology computations
   gmsh::model::mesh::generate(3);
 
+  // For more information, see M. Pellikka, S. Suuriniemi, L. Kettunen and
+  // C. Geuzaine. Homology and cohomology computation in finite element
+  // modeling. SIAM Journal on Scientific Computing 35(5), pp. 1195-1214, 2013.
+
   gmsh::finalize();
   return 0;
 }
diff --git a/tutorial/c++/t15.cpp b/tutorial/c++/t15.cpp
index 1c79a733a4529b7956566b662bcd0ae2c04eb3f2..300d103b2502cc5909c401528e92003300e565cf 100644
--- a/tutorial/c++/t15.cpp
+++ b/tutorial/c++/t15.cpp
@@ -1,4 +1,18 @@
-// This file reimplements gmsh/tutorial/t15.geo in C++.
+/*******************************************************************************
+ *
+ *  Gmsh C++ tutorial 15
+ *
+ *  Embedded points, lines and surfaces
+ *
+ *******************************************************************************/
+
+// By default, across geometrical dimensions meshes generated by Gmsh are only
+// conformal if lower dimensional entities are on the boundary of higher
+// dimensional ones (i.e. if points, curves or surfaces are part of the boundary
+// of volumes).
+
+// Embedding constraints allow to force a mesh to be conformal to other lower
+// dimensional entities.
 
 #include <gmsh.h>
 
@@ -12,7 +26,7 @@ int main(int argc, char **argv)
 
   gmsh::model::add("t15");
 
-  // Copied from t1.cpp...
+  // Copied from t1.cpp:
   double lc = 1e-2;
   factory::addPoint(0, 0, 0, lc, 1);
   factory::addPoint(.1, 0,  0, lc, 2);
@@ -28,23 +42,23 @@ int main(int argc, char **argv)
   model::addPhysicalGroup(1, {1, 2}, 2);
   model::addPhysicalGroup(2, {1}, 6);
   model::setPhysicalName(2, 6, "My surface");
-  // ...end of copy
 
   // We change the mesh size to generate a coarser mesh
   lc *=  4;
   factory::mesh::setSize({{0, 1}, {0, 2}, {0, 3}, {0, 4}}, lc);
 
-  // Define a new point and embed it in a surface
+  // We define a new point
   factory::addPoint(0.02, 0.02, 0., lc, 5);
 
-  // We have to synchronize before embedding entites.
-  // Otherwise, we get an error like "Point 5 does not exist"
+  // We have to synchronize before embedding entites:
   factory::synchronize();
 
-  // embed the point (dim 0) in the surface (dim 2)
+  // One can force this point to be included ("embedded") in the 2D mesh, using
+  // the `embed()' function:
   model::mesh::embed(0, {5}, 2, 1);
 
-  // We can also use embed to embed a curve in the 2D mesh
+  // In the same way, one use `embed()' to force a curve to be embedded in the
+  // 2D mesh:
   factory::addPoint(0.02, 0.12, 0., lc, 6);
   factory::addPoint(0.04, 0.18, 0., lc, 7);
   factory::addLine(6, 7, 5);
@@ -66,7 +80,7 @@ int main(int argc, char **argv)
   factory::synchronize();
   model::mesh::embed(1, {l}, 3, 1);
 
-  // finally, we can embed a surface in a volume
+  // Finally, we can also embed a surface in a volume:
   factory::addPoint(0.02, 0.12, 0.05, lc, p+2);
   factory::addPoint(0.04, 0.12, 0.05, lc, p+3);
   factory::addPoint(0.04, 0.18, 0.05, lc, p+4);
@@ -83,7 +97,6 @@ int main(int argc, char **argv)
   factory::synchronize();
   model::mesh::embed(2, {s}, 3, 1);
 
-  // create and show the mesh
   model::mesh::generate(3);
 
   gmsh::write("t15.msh");
diff --git a/tutorial/c++/t16.cpp b/tutorial/c++/t16.cpp
index b267b3ae5f807d9f7955d8f40745b1bea1fc3f5d..b249fc20957d0b7d1e6eeff6f9349b3651a0d7c1 100644
--- a/tutorial/c++/t16.cpp
+++ b/tutorial/c++/t16.cpp
@@ -1,4 +1,14 @@
-// This file reimplements gmsh/tutorial/t16.geo in C++.
+/*******************************************************************************
+ *
+ *  Gmsh C++ tutorial 16
+ *
+ *  Constructive Solid Geometry, OpenCASCADE geometry kernel
+ *
+ *******************************************************************************/
+
+// Instead of constructing a model in a bottom-up fashion with Gmsh's built-in
+// geometry kernel, starting with version 3 Gmsh allows you to directly use
+// alternative geometry kernels. Here we will use the OpenCASCADE kernel.
 
 #include <iostream>
 #include <gmsh.h>
@@ -13,21 +23,32 @@ int main(int argc, char **argv)
 
   model::add("t16");
 
-  //gmsh::logger::start();
+  // Let's build the same model as in `t5.geo', but using constructive solid
+  // geometry.
 
-  std::vector<std::pair<int, int> > ov;
-  std::vector<std::vector<std::pair<int, int> > > ovv;
+  // We can log all messages for further processing with:
+  gmsh::logger::start();
 
+  // We first create two cubes:
   try {
     factory::addBox(0,0,0, 1,1,1, 1);
+    factory::addBox(0,0,0, 0.5,0.5,0.5, 2);
   }
   catch(...) {
     gmsh::logger::write("Could not create OpenCASCADE shapes: bye!");
     return 0;
   }
 
-  factory::addBox(0,0,0, 0.5,0.5,0.5, 2);
+  // We apply a boolean difference to create the "cube minus one eigth" shape:
+  std::vector<std::pair<int, int> > ov;
+  std::vector<std::vector<std::pair<int, int> > > ovv;
   factory::cut({{3,1}}, {{3,2}}, ov, ovv, 3);
+
+  // Boolean operations with OpenCASCADE always create new entities. By default
+  // the extra arguments `removeObject' and `removeTool' in `occ::cut()' are set
+  // to `true', which will delete the original entities.
+
+  // We then create the five spheres:
   double x = 0, y = 0.75, z = 0, r = 0.09 ;
   std::vector<std::pair<int, int> > holes;
   for(int t = 1; t <= 5; t++){
@@ -36,20 +57,52 @@ int main(int argc, char **argv)
     factory::addSphere(x,y,z,r, 3 + t);
     holes.push_back({3, 3 + t});
   }
+
+  // If we had wanted five empty holes we would have used `cut()' again. Here we
+  // want five spherical inclusions, whose mesh should be conformal with the
+  // mesh of the cube: we thus use `fragment()', which intersects all volumes in
+  // a conformal manner (without creating duplicate interfaces):
   factory::fragment({{3,3}}, holes, ov, ovv);
 
   factory::synchronize();
 
+  // When the boolean operation leads to simple modifications of entities, and
+  // if one deletes the original entities, Gmsh tries to assign the same tag to
+  // the new entities. (This behavior is governed by the
+  // `Geometry.OCCBooleanPreserveNumbering' option.)
+
+  // Here the `Physical Volume' definitions can thus be made for the 5 spheres
+  // directly, as the five spheres (volumes 4, 5, 6, 7 and 8), which will be
+  // deleted by the fragment operations, will be recreated identically (albeit
+  // with new surfaces) with the same tags.:
+  for(int i = 1; i <= 5; i++)
+    gmsh::model::addPhysicalGroup(3, {3 + i}, i);
+
+  // The tag of the cube will change though, so we need to access it
+  // programmatically:
+  gmsh::model::addPhysicalGroup(3, {ov[ov.size() - 1].second}, 10);
+
+  // Creating entities using constructive solid geometry is very powerful, but can
+  // lead to small pratical issues for e.g. setting mesh sizes at points, or
+  // identifying boundaries.
+
+  // To identify point or other bounding entities you can take advantage of the
+  // `getEntities()', `getBoundary()' and `getEntitiesInBoundingBox()'
+  // functions:
+
   double lcar1 = .1;
   double lcar2 = .0005;
   double lcar3 = .055;
 
+  // Assign a mesh size to all the points:
   model::getEntities(ov, 0);
   model::mesh::setSize(ov, lcar1);
 
+  // Override this constraint on the points of the five spheres:
   model::getBoundary(holes, ov, false, false, true);
   model::mesh::setSize(ov, lcar3);
 
+  // Select the corner point by searching for it geometrically:
   double eps = 1e-3;
   model::getEntitiesInBoundingBox(0.5-eps, 0.5-eps, 0.5-eps,
                                   0.5+eps, 0.5+eps, 0.5+eps, ov, 0);
@@ -59,10 +112,13 @@ int main(int argc, char **argv)
 
   gmsh::write("t16.msh");
 
-  //std::vector<std::string> log;
-  //gmsh::logger::get(log);
-  //std::cout << "Logger has recorded " << log.size() << " lines" << std::endl;
-  //gmsh::logger::stop();
+  // Additional examples created with the OpenCASCADE geometry kernel are
+  // available in `t18.cpp', `t19.cpp' as well as in the `demos/api' directory.
+
+  std::vector<std::string> log;
+  gmsh::logger::get(log);
+  std::cout << "Logger has recorded " << log.size() << " lines" << std::endl;
+  gmsh::logger::stop();
 
   gmsh::finalize();
   return 0;
diff --git a/tutorial/c++/t17.cpp b/tutorial/c++/t17.cpp
index bf13760195c5b721d2dc5e13e7a6fa3f1aace60e..5bcd7c4ea7f6d97faf7df94da452d1db232755dc 100644
--- a/tutorial/c++/t17.cpp
+++ b/tutorial/c++/t17.cpp
@@ -1,4 +1,18 @@
-// This file reimplements gmsh/tutorial/t17.geo in C++.
+/*******************************************************************************
+ *
+ *  Gmsh C++ tutorial 17
+ *
+ *  Anisotropic background mesh
+ *
+ *******************************************************************************/
+
+// Characteristic lengths can be specified very accurately by providing a
+// background mesh, i.e., a post-processing view that contains the target mesh
+// sizes.
+
+// Here, the background mesh is represented as a metric tensor field defined on
+// a square. One should use bamg as 2d mesh generator to enable anisotropic
+// meshes in 2D.
 
 #include <gmsh.h>
 #include <math.h>
@@ -11,10 +25,13 @@ int main(int argc, char **argv)
   gmsh::initialize();
   gmsh::option::setNumber("General.Terminal", 1);
 
+  gmsh::model::add("t17");
+
+  // Create a square
   gmsh::model::occ::addRectangle(-1, -1, 0, 2, 2);
   gmsh::model::occ::synchronize();
 
-  // add a post-processing view to use as a size field
+  // Merge a post-processing view containing the target anisotropic mesh sizes
   try {
     gmsh::merge("../t17_bgmesh.pos");
   }
@@ -24,10 +41,11 @@ int main(int argc, char **argv)
     return 0;
   }
 
+  // Apply the view as the current background mesh
   int bg_field = gmsh::model::mesh::field::add("PostView");
   gmsh::model::mesh::field::setAsBackgroundMesh(bg_field);
 
-  // use bamg
+  // Use bamg
   gmsh::option::setNumber("Mesh.SmoothRatio", 3);
   gmsh::option::setNumber("Mesh.AnisoMax", 1000);
   gmsh::option::setNumber("Mesh.Algorithm", 7);
diff --git a/tutorial/c++/t18.cpp b/tutorial/c++/t18.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..373af42a25d8b3e179e00876a1144b45f6da4443
--- /dev/null
+++ b/tutorial/c++/t18.cpp
@@ -0,0 +1,135 @@
+/*******************************************************************************
+ *
+ *  Gmsh GEO tutorial 18
+ *
+ *  Periodic meshes
+ *
+ *******************************************************************************/
+
+// Periodic meshing constraints can be imposed on surfaces and curves.
+
+#include <gmsh.h>
+
+int main(int argc, char **argv)
+{
+  gmsh::initialize(argc, argv);
+  gmsh::option::setNumber("General.Terminal", 1);
+
+  gmsh::model::add("t18");
+
+  // Let's use the OpenCASCADE geometry kernel to build two geometries.
+
+  // The first geometry is very simple: a unit cube with a non-uniform mesh size
+  // constraint (set on purpose to be able to verify visually that the
+  // periodicity constraint works!):
+
+  gmsh::model::occ::addBox(0, 0, 0, 1, 1, 1, 1);
+  gmsh::model::occ::synchronize();
+
+  std::vector<std::pair<int, int> > out;
+  gmsh::model::getEntities(out, 0);
+  gmsh::model::mesh::setSize(out, 0.1);
+  gmsh::model::mesh::setSize({{0,1}}, 0.02);
+
+  // To impose that the mesh on surface 2 (the right side of the cube) should
+  // match the mesh from surface 1 (the left side), the following periodicity
+  // constraint is set:
+  std::vector<double> translation({1,0,0,1, 0,1,0,0, 0,0,1,0, 0,0,0,0});
+  gmsh::model::mesh::setPeriodic(2, {2}, {1}, translation);
+
+  // The periodicity transform is provided as a 4x4 affine transformation
+  // matrix, given by row.
+
+  // During mesh generation, the mesh on surface 2 will be created by copying
+  // the mesh from surface 1.
+
+  // For more complicated cases, finding the corresponding surfaces by hand can
+  // be tedious, especially when geometries are created through solid
+  // modelling. Let's construct a slightly more complicated geometry.
+
+  // We start with a cube and some spheres:
+  gmsh::model::occ::addBox(2, 0, 0, 1, 1, 1, 10);
+  double x = 2-0.3, y = 0, z = 0;
+  gmsh::model::occ::addSphere(x, y, z, 0.35, 11);
+  gmsh::model::occ::addSphere(x+1, y, z, 0.35, 12);
+  gmsh::model::occ::addSphere(x, y+1, z, 0.35, 13);
+  gmsh::model::occ::addSphere(x, y, z+1, 0.35, 14);
+  gmsh::model::occ::addSphere(x+1, y+1, z, 0.35, 15);
+  gmsh::model::occ::addSphere(x, y+1, z+1, 0.35, 16);
+  gmsh::model::occ::addSphere(x+1, y, z+1, 0.35, 17);
+  gmsh::model::occ::addSphere(x+1, y+1, z+1, 0.35, 18);
+
+  // We first fragment all the volumes, which will leave parts of spheres
+  // protruding outside the cube:
+  std::vector<std::vector<std::pair<int, int> > > out_map;
+  std::vector<std::pair<int, int> > sph;
+  for(int i = 11; i <= 18; i++) sph.push_back(std::pair<int, int>(3, i));
+  gmsh::model::occ::fragment({{3,10}}, sph, out, out_map);
+  gmsh::model::occ::synchronize();
+
+  // Ask OpenCASCADE to compute more accurate bounding boxes of entities using
+  // the STL mesh:
+  gmsh::option::setNumber("Geometry.OCCBoundsUseStl", 1);
+
+  // We then retrieve all the volumes in the bounding box of the original cube,
+  // and delete all the parts outside it:
+  double eps = 1e-3;
+  std::vector<std::pair<int, int> > in;
+  gmsh::model::getEntitiesInBoundingBox(2-eps,-eps,-eps, 2+1+eps,1+eps,1+eps,
+                                        in, 3);
+  std::vector<int> boundary_tags, complement_tags;
+  for(auto it = in.begin(); it != in.end(); ++it) {
+    auto it2 = std::find(out.begin(), out.end(), *it);
+    if(it2 != out.end()) out.erase(it2);
+  }
+  gmsh::model::removeEntities(out, true); // Delete outside parts recursively
+
+  // We now set some a non-uniform mesh size constraint (again to check results
+  // visually):
+  std::vector<std::pair<int, int> > p;
+  gmsh::model::getBoundary(in, p, false, false, true); // Get all points
+  gmsh::model::mesh::setSize(p, 0.1);
+  gmsh::model::getEntitiesInBoundingBox(2-eps, -eps, -eps, 2+eps, eps, eps,
+                                        p, 0);
+  gmsh::model::mesh::setSize(p, 0.001);
+
+  // We now identify corresponding surfaces on the left and right sides of the
+  // geometry automatically.
+
+  // First we get all surfaces on the left:
+  std::vector<std::pair<int, int> > sxmin;
+  gmsh::model::getEntitiesInBoundingBox(2-eps, -eps, -eps, 2+eps, 1+eps, 1+eps,
+                                        sxmin, 2);
+  for(std::size_t i = 0; i < sxmin.size(); i++) {
+    // Then we get the bounding box of each left surface
+    double xmin, ymin, zmin, xmax, ymax, zmax;
+    gmsh::model::getBoundingBox(sxmin[i].first, sxmin[i].second,
+                                xmin, ymin, zmin, xmax, ymax, zmax);
+    // We translate the bounding box to the right and look for surfaces inside
+    // it:
+    std::vector<std::pair<int, int> > sxmax;
+    gmsh::model::getEntitiesInBoundingBox(xmin-eps+1, ymin-eps, zmin-eps,
+                                          xmax+eps+1, ymax+eps, zmax+eps,
+                                          sxmax, 2);
+    // For all the matches, we compare the corresponding bounding boxes...
+    for(std::size_t j = 0; j < sxmax.size(); j++) {
+      double xmin2, ymin2, zmin2, xmax2, ymax2, zmax2;
+      gmsh::model::getBoundingBox(sxmax[j].first, sxmax[j].second,
+                                  xmin2, ymin2, zmin2, xmax2, ymax2, zmax2);
+      xmin2 -= 1;
+      xmax2 -= 1;
+      // ...and if they match, we apply the periodicity constraint
+      if(std::abs(xmin2 - xmin) < eps && std::abs(xmax2 - xmax) < eps &&
+         std::abs(ymin2 - ymin) < eps && std::abs(ymax2 - ymax) < eps &&
+         std::abs(zmin2 - zmin) < eps && std::abs(zmax2 - zmax) < eps) {
+        gmsh::model::mesh::setPeriodic(2, {sxmax[j].second}, {sxmin[i].second},
+                                       translation);
+      }
+    }
+  }
+
+  // gmsh::fltk::run();
+
+  gmsh::finalize();
+  return 0;
+}
diff --git a/tutorial/c++/t19.cpp b/tutorial/c++/t19.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a5994e35828439a33927e7eed8ce9c03615da5b9
--- /dev/null
+++ b/tutorial/c++/t19.cpp
@@ -0,0 +1,102 @@
+/*******************************************************************************
+ *
+ *  Gmsh GEO tutorial 19
+ *
+ *  Thrusections, fillets, pipes, mesh size from curvature
+ *
+ *******************************************************************************/
+
+// The OpenCASCADE geometry kernel supports several useful features for solid
+// modelling.
+
+#include <cmath>
+#include <gmsh.h>
+
+int main(int argc, char **argv)
+{
+  gmsh::initialize(argc, argv);
+  gmsh::option::setNumber("General.Terminal", 1);
+
+  gmsh::model::add("t19");
+
+  // Volumes can be constructed from (closed) curve loops thanks to the
+  // `addThruSections()' function
+  gmsh::model::occ::addCircle(0,0,0, 0.5, 1);
+  gmsh::model::occ::addCurveLoop({1}, 1);
+  gmsh::model::occ::addCircle(0.1,0.05,1, 0.1, 2);
+  gmsh::model::occ::addCurveLoop({2}, 2);
+  gmsh::model::occ::addCircle(-0.1,-0.1,2, 0.3, 3);
+  gmsh::model::occ::addCurveLoop({3}, 3);
+  std::vector<std::pair<int, int> > out;
+  gmsh::model::occ::addThruSections({1,2,3}, out, 1);
+  gmsh::model::occ::synchronize();
+
+  // With `Ruled ThruSections' you can force the use of ruled surfaces:
+  gmsh::model::occ::addCircle(2+0,0,0, 0.5, 11);
+  gmsh::model::occ::addCurveLoop({11}, 11);
+  gmsh::model::occ::addCircle(2+0.1,0.05,1, 0.1, 12);
+  gmsh::model::occ::addCurveLoop({12}, 12);
+  gmsh::model::occ::addCircle(2-0.1,-0.1,2, 0.3, 13);
+  gmsh::model::occ::addCurveLoop({13}, 13);
+  gmsh::model::occ::addThruSections({11,12,13}, out, 11, true, true);
+  gmsh::model::occ::synchronize();
+
+  // We copy the first volume, and fillet all its edges:
+  gmsh::model::occ::copy({{3, 1}}, out);
+  gmsh::model::occ::translate(out, 4, 0, 0);
+  gmsh::model::occ::synchronize();
+  std::vector<std::pair<int, int> > f;
+  gmsh::model::getBoundary(out, f);
+  std::vector<std::pair<int, int> > e;
+  gmsh::model::getBoundary(f, e, false);
+  std::vector<int> c;
+  for(std::size_t i = 0; i < e.size(); i++) c.push_back(std::abs(e[i].second));
+  gmsh::model::occ::fillet({out[0].second}, c, {0.1}, out);
+  gmsh::model::occ::synchronize();
+
+  // OpenCASCADE also allows general extrusions along a smooth path. Let's first
+  // define a spline curve:
+  double nturns = 1;
+  int npts = 20;
+  double r = 1;
+  double h = 1 * nturns;
+  std::vector<int> p;
+  for(int i = 0; i < npts; i++) {
+    double theta = i * 2*M_PI*nturns/npts;
+    gmsh::model::occ::addPoint(r * cos(theta), r * sin(theta), i * h/npts, 1,
+                               1000 + i);
+    p.push_back(1000 + i);
+  }
+  gmsh::model::occ::addSpline(p, 1000);
+
+  // A wire is like a curve loop, but open:
+  gmsh::model::occ::addWire({1000}, 1000);
+
+  // We define the shape we would like to extrude along the spline (a disk):
+  gmsh::model::occ::addDisk(1,0,0, 0.2, 0.2, 1000);
+  gmsh::model::occ::rotate({{2,1000}}, 0, 0, 0, 1, 0, 0, M_PI/2);
+
+  // We extrude the disk along the spline to create a pipe:
+  gmsh::model::occ::addPipe({{2,1000}}, 1000, out);
+
+  // We delete the source surface, and increase the number of sub-edges for a
+  // nicer display of the geometry:
+  gmsh::model::occ::remove({{2, 1000}});
+  gmsh::option::setNumber("Geometry.NumSubEdges", 1000);
+
+  gmsh::model::occ::synchronize();
+
+  // We can activate the calculation of mesh element sizes based on curvature:
+  gmsh::option::setNumber("Mesh.CharacteristicLengthFromCurvature", 1);
+
+  // And we set the minimum number of elements per 2*Pi radians:
+  gmsh::option::setNumber("Mesh.MinimumElementsPerTwoPi", 20);
+
+  // We can constraints the min and max element sizes to stay within reasonnable
+  // values (see `t10.geo' for more details):
+  gmsh::option::setNumber("Mesh.CharacteristicLengthMin", 0.001);
+  gmsh::option::setNumber("Mesh.CharacteristicLengthMax", 0.3);
+
+  // gmsh::fltk::run();
+  return 0;
+}
diff --git a/tutorial/c++/t2.cpp b/tutorial/c++/t2.cpp
index 6b415e933fbef61b45deecf6ad2282539a8d28df..eccaa7ce2c9ad97bbd573e343138a82e231f8e8f 100644
--- a/tutorial/c++/t2.cpp
+++ b/tutorial/c++/t2.cpp
@@ -1,17 +1,24 @@
-// This file reimplements gmsh/tutorial/t2.geo in C++. Comments focus on the new
-// API functions used, compared to the ones introduced in t1.cpp.
+/*******************************************************************************
+ *
+ *  Gmsh C++ tutorial 2
+ *
+ *  Geometrical transformations, extruded geometries, elementary entities
+ *  (volumes), physical groups (volumes)
+ *
+ *******************************************************************************/
 
 #include <gmsh.h>
 
-// give nice shortcuts for some namespaces
+// We start by giving some nice shortcuts for some namespaces
 namespace model = gmsh::model;
 namespace factory = gmsh::model::geo;
 
 int main(int argc, char **argv)
 {
-  // If argc/argv are passed, Gmsh will parse the commandline in the same way as
-  // the standalone Gmsh app.
+  // If argc/argv are passed to gmsh::initialize(), Gmsh will parse the
+  // commandline in the same way as the standalone Gmsh app:
   gmsh::initialize(argc, argv);
+
   gmsh::option::setNumber("General.Terminal", 1);
 
   model::add("t2");
@@ -28,38 +35,53 @@ int main(int argc, char **argv)
   factory::addLine(4, 1, 4);
   factory::addCurveLoop({4, 1, -2, 3}, 1);
   factory::addPlaneSurface({1}, 1);
-  model::addPhysicalGroup(0, {1, 2}, 1);
-  model::addPhysicalGroup(1, {1, 2}, 2);
-  model::addPhysicalGroup(2, {1}, 6);
-  model::setPhysicalName(2, 6, "My surface");
-  // ...end of copy
+  model::addPhysicalGroup(1, {1, 2, 4}, 5);
+  int ps = model::addPhysicalGroup(2, {1});
+  model::setPhysicalName(2, ps, "My surface");
 
+  // We can then add new points and curves in the same way as we did in
+  // `t1.cpp':
   factory::addPoint(0, .4, 0, lc, 5);
   factory::addLine(4, 5, 5);
 
-  // Geometrical transformations take a vector of pairs of integers as first
-  // argument, which contains the list of entities, represented by (dimension,
-  // tag) pairs. Here we thus translate point 3 (dimension=0, tag=3), by
-  // dx=-0.05, dy=0, dz=0.
+  // But Gmsh also provides tools to transform (translate, rotate, etc.)
+  // elementary entities or copies of elementary entities.  Geometrical
+  // transformations take a vector of pairs of integers as first argument, which
+  // contains the list of entities, represented by (dimension, tag) pairs.  For
+  // example, the point 3 (dimension=0, tag=3) can be moved by 0.05 to the left
+  // (dx=-0.05, dy=0, dz=0) with
   factory::translate({{0, 3}}, -0.05, 0, 0);
 
-  // The "Duplicata" functionality in .geo files is handled by
+  // Note that there are no units in Gmsh: coordinates are just numbers - it's
+  // up to the user to associate a meaning to them.
+
+  // The resulting point can also be duplicated and translated by 0.1 along the
+  // y axis. The "Duplicata" functionality in .geo files is handled by
   // model::geo::copy(), which takes a vector of (dim, tag) pairs as input, and
-  // returns another vector of (dim, tag) pairs.
+  // returns another vector of (dim, tag) pairs:
   std::vector<std::pair<int, int> > ov;
   factory::copy({{0, 3}}, ov);
   factory::translate(ov, 0, 0.1, 0);
 
+  // The new point tag is available in ov[0].second, and can be used to create
+  // new lines:
   factory::addLine(3, ov[0].second, 7);
   factory::addLine(ov[0].second, 5, 8);
   factory::addCurveLoop({5,-8,-7,3}, 10);
   factory::addPlaneSurface({10}, 11);
 
+  // In the same way, we can translate copies of the two surfaces 6 and 11 to
+  // the right with the following command:
   factory::copy({{2, 1}, {2, 11}}, ov);
   factory::translate(ov, 0.12, 0, 0);
 
   std::printf("New surfaces '%d' and '%d'\n", ov[0].second, ov[1].second);
 
+  // Volumes are the fourth type of elementary entities in Gmsh. In the same way
+  // one defines curve loops to build surfaces, one has to define surface loops
+  // (i.e. `shells') to build volumes. The following volume does not have holes
+  // and thus consists of a single surface loop:
+
   factory::addPoint(0., 0.3, 0.13, lc, 100);
   factory::addPoint(0.08, 0.3, 0.1, lc, 101);
   factory::addPoint(0.08, 0.4, 0.1, lc, 102);
@@ -85,30 +107,58 @@ int main(int argc, char **argv)
   factory::addCurveLoop({115, 116, 117, 114}, 126);
   factory::addPlaneSurface({126}, 127);
 
-  // The API to create surface loops ("shells") and volumes is similar to the
-  // one used to create curve loops and surfaces.
   factory::addSurfaceLoop({127, 119, 121, 123, 125, 11}, 128);
   factory::addVolume({128}, 129);
 
-  // Extrusion works as expected, by providing a vector of (dim, tag) pairs as
-  // input, the translation vector, and a vector of (dim, tag) pairs as output.
+  // When a volume can be extruded from a surface, it is usually easier to use
+  // the `extrude()' function directly instead of creating all the points,
+  // curves and surfaces by hand. For example, the following command extrudes
+  // the surface 11 along the z axis and automatically creates a new volume (as
+  // well as all the needed points, curves and surfaces). As expected, the
+  // function takes a vector of (dim, tag) pairs as input as well as the
+  // translation vector, and returns a vector of (dim, tag) pairs as output:
   std::vector<std::pair<int, int> > ov2;
   factory::extrude({ov[1]}, 0, 0, 0.12, ov2);
 
   // Mesh sizes associated to geometrical points can be set by passing a vector
-  // of (dim, tag) pairs for the corresponding points.
+  // of (dim, tag) pairs for the corresponding points:
   factory::mesh::setSize({{0,103}, {0,105}, {0,109}, {0,102}, {0,28},
                           {0, 24}, {0,6}, {0,5}}, lc * 3);
 
-  model::addPhysicalGroup(3, {129,130}, 1);
+  // We finally group volumes 129 and 130 in a single physical group with tag
+  // `1' and name "The volume":
+  model::addPhysicalGroup(3, {129, 130}, 1);
   model::setPhysicalName(3, 1, "The volume");
 
+  // We finish by synchronizing the data from the built-in geometry kernel with
+  // the Gmsh model, and by generating and saving the mesh:
   factory::synchronize();
-
   model::mesh::generate(3);
-
   gmsh::write("t2.msh");
 
+  // Note that, if the transformation tools are handy to create complex
+  // geometries, it is also sometimes useful to generate the `flat' geometry, with
+  // an explicit representation of all the elementary entities.
+  //
+  // With the built-in geometry kernel, this can be achieved by saving the model
+  // in the `Gmsh Unrolled GEO' format:
+  //
+  // gmsh::write("t2.geo_unrolled");
+  //
+  // With the OpenCASCADE geometry kernel, unrolling the geometry can be
+  // achieved by exporting in the `OpenCASCADE BRep' format:
+  //
+  // gmsh::write("t2.brep");
+  //
+  // (OpenCASCADE geometries can also be exported as STEP files.)
+
+  // It is important to note that Gmsh never translates geometry data into a
+  // common representation: all the operations on a geometrical entity are
+  // performed natively with the associated geometry kernel. Consequently, one
+  // cannot export a geometry constructed with the built-in kernel as an
+  // OpenCASCADE BRep file; or export an OpenCASCADE model as an Unrolled GEO
+  // file.
+
   gmsh::finalize();
   return 0;
 }
diff --git a/tutorial/c++/t3.cpp b/tutorial/c++/t3.cpp
index 67031c8c204680c50a40b84ccd9999939c6fa04a..5fb4e9ed67a60dd1cc7fe3c0432dfeb1b5add844 100644
--- a/tutorial/c++/t3.cpp
+++ b/tutorial/c++/t3.cpp
@@ -1,4 +1,10 @@
-// This files reimplements gmsh/tutorial/t3.geo in C++.
+/*******************************************************************************
+ *
+ *  Gmsh C++ tutorial 3
+ *
+ *  Extruded meshes, options
+ *
+ *******************************************************************************/
 
 #include <cmath>
 #include <gmsh.h>
@@ -29,29 +35,76 @@ int main(int argc, char **argv)
   model::addPhysicalGroup(1, {1, 2}, 2);
   model::addPhysicalGroup(2, {1}, 6);
   model::setPhysicalName(2, 6, "My surface");
-  // ...end of copy
+
+  // As in `t2.cpp', we plan to perform an extrusion along the z axis.  But
+  // here, instead of only extruding the geometry, we also want to extrude the
+  // 2D mesh. This is done with the same `extrude()' function, but by specifying
+  // element 'Layers' (2 layers in this case, the first one with 8 subdivisions
+  // and the second one with 2 subdivisions, both with a height of h/2). The
+  // number of elements for each layer and the (end) height of each layer are
+  // specified in two vectors:
 
   double h = 0.1, angle = 90.;
   std::vector<std::pair<int, int> > ov;
-
-  // Extruding the mesh in addition to the geometry works as in .geo files: the
-  // number of elements for each layer and the (end) height of each layer are
-  // specified in two vectors.
   factory::extrude({{2,1}}, 0, 0, h, ov, {8,2}, {0.5,1});
 
-  // Rotational and twisted extrusions are available as well with the built-in
-  // CAD kernel. The last (optional) argument for the Extrude/Revolve/Twist
-  // commands specifies whether the extruded mesh should be recombined or not.
+  // The extrusion can also be performed with a rotation instead of a
+  // translation, and the resulting mesh can be recombined into prisms (we use
+  // only one layer here, with 7 subdivisions). All rotations are specified by
+  // an an axis point (-0.1, 0, 0.1), an axis direction (0, 1, 0), and a
+  // rotation angle (-Pi/2):
   factory::revolve({{2,28}}, -0.1,0,0.1, 0,1,0, -M_PI/2, ov, {7});
+
+  // Using the built-in geometry kernel, only rotations with angles < Pi are
+  // supported. To do a full turn, you will thus need to apply at least 3
+  // rotations. The OpenCASCADE geometry kernel does not have this limitation.
+
+  // A translation (-2*h, 0, 0) and a rotation ((0,0.15,0.25}, (1,0,0), Pi/2)
+  // can also be combined to form a "twist".  The last (optional) argument for
+  // the extrude() and twist() functions specifies whether the extruded mesh
+  // should be recombined or not.
   factory::twist({{2,50}}, 0,0.15,0.25, -2*h,0,0, 1,0,0, angle*M_PI/180.,
                  ov, {10}, {}, true);
 
   factory::synchronize();
 
+  // All the extrusion functions return a vector of extruded entities: the "top"
+  // of the extruded surface (in `ov[0]'), the newly created volume (in `ov[1]')
+  // and the tags of the lateral surfaces (in `ov[2]', `ov[3]', ...).
+
+  // We can then define a new physical volume (with tag 101) to group all the
+  // elementary volumes:
   model::addPhysicalGroup(3, {1, 2, ov[1].second}, 101);
 
   model::mesh::generate(3);
   gmsh::write("t3.msh");
+
+  // Let us now change some options... Since all interactive options are
+  // accessible through the API, we can for example make point tags visible or
+  // redefine some colors:
+
+  gmsh::option::setNumber("Geometry.PointNumbers", 1);
+  gmsh::option::setColor("Geometry.Points", 255, 165, 0);
+  gmsh::option::setColor("General.Text", 255, 255, 255);
+  gmsh::option::setColor("Mesh.Points", 255, 0, 0);
+
+  // Note that color options are special: setting a color option of "X.Y"
+  // actually sets the option "X.Color.Y".
+
+  int r, g, b, a;
+  gmsh::option::getColor("Geometry.Points", r, g, b, a);
+  gmsh::option::setColor("Geometry.Surfaces", r, g, b, a);
+
+  // Launch the GUI to see the effects of the color changes:
+
+  // gmsh::fltk::run();
+
+  // When the GUI is launched, you can use the `Help->Current options' menu to
+  // see the current values of all options. To save the options in a file, use
+  // `File->Export->Gmsh options', or through the api:
+
+  // gmsh::write("t3.opt");
+
   gmsh::finalize();
 
   return 0;
diff --git a/tutorial/c++/t4.cpp b/tutorial/c++/t4.cpp
index 82d28e3f0b7383b9642421a7f5b2d6037865ba28..26412f5a12f2586334929e315132f8b82bbecafe 100644
--- a/tutorial/c++/t4.cpp
+++ b/tutorial/c++/t4.cpp
@@ -1,4 +1,10 @@
-// This file reimplements gmsh/tutorial/t4.geo in C++.
+/*******************************************************************************
+ *
+ *  Gmsh C++ tutorial 4
+ *
+ *  Holes in surfaces, annotations, entity colors
+ *
+ *******************************************************************************/
 
 #include <math.h>
 #include <gmsh.h>
@@ -25,6 +31,7 @@ int main(int argc, char **argv)
   double ccos = (-h5*R1 + e2 * hypot(h5, hypot(e2, R1))) / (h5*h5 + e2*e2);
   double ssin = sqrt(1 - ccos*ccos);
 
+  // We start by defining some points and some lines:
   factory::addPoint(-e1-e2, 0    , 0, Lc1, 1);
   factory::addPoint(-e1-e2, h1   , 0, Lc1, 2);
   factory::addPoint(-e3-r , h1   , 0, Lc2, 3);
@@ -57,7 +64,16 @@ int main(int argc, char **argv)
   factory::addLine(1 , 17, 1);
   factory::addLine(17, 16, 2);
 
+  // Gmsh provides other curve primitives than straight lines: splines,
+  // B-splines, circle arcs, ellipse arcs, etc. Here we define a new circle arc,
+  // starting at point 14 and ending at point 16, with the circle's center being
+  // the point 15:
   factory::addCircleArc(14,15,16, 3);
+
+  // Note that, in Gmsh, circle arcs should always be smaller than Pi. The
+  // OpenCASCADE geometry kernel does not have this limitation.
+
+  // We can then define additional lines and circles, as well as a new surface:
   factory::addLine(14,13, 4);
   factory::addLine(13,12, 5);
   factory::addLine(12,11, 6);
@@ -78,29 +94,83 @@ int main(int argc, char **argv)
 
   factory::addCurveLoop({17,-15,18,19,-20,16}, 21);
   factory::addPlaneSurface({21}, 22);
-  factory::addCurveLoop({11,-12,13,14,1,2,-3,4,5,6,7,-8,9,10}, 23);
 
-  // A surface with one hole is specified using 2 curve loops:
+  // But we still need to define the exterior surface. Since this surface has a
+  // hole, its definition now requires two curves loops:
+  factory::addCurveLoop({11,-12,13,14,1,2,-3,4,5,6,7,-8,9,10}, 23);
   factory::addPlaneSurface({23,21}, 24);
 
-  // FIXME: this will be implemented through the gmshView API
-  /*
-  View "comments" {
-    T2(10, -10, 0){ StrCat("Created on ", Today, " with Gmsh") };
-    T3(0, 0.11, 0, TextAttributes("Align", "Center", "Font", "Helvetica")){ "Hole" };
-    T3(0, 0.09, 0, TextAttributes("Align", "Center")){ "file://image.png@0.01x0" };
-    T3(-0.01, 0.09, 0, 0){ "file://image.png@0.01x0,0,0,1,0,1,0" };
-    T3(0, 0.12, 0, TextAttributes("Align", "Center")){ "file://image.png@0.01x0#" };
-    T2(350, -7, 0){ "file://image.png@20x0" };
-  };
-  */
+  // As a general rule, if a surface has N holes, it is defined by N+1 curve loops:
+  // the first loop defines the exterior boundary; the other loops define the
+  // boundaries of the holes.
 
   factory::synchronize();
 
+  // Finally, we can add some comments by creating a post-processing view
+  // containing some strings:
+  int v = gmsh::view::add("comments");
+
+  // Add a text string in window coordinates, 10 pixels from the left and 10
+  // pixels from the bottom:
+  gmsh::view::addListDataString(v, {10, -10}, {"Created with Gmsh"});
+
+  // Add a text string in model coordinates centered at (X,Y,Z) = (0, 0.11, 0),
+  // with some style attributes:
+  gmsh::view::addListDataString(v, {0, 0.11, 0}, {"Hole"},
+                                {"Align", "Center", "Font", "Helvetica"});
+
+  // If a string starts with `file://', the rest is interpreted as an image
+  // file. For 3D annotations, the size in model coordinates can be specified
+  // after a `@' symbol in the form `widthxheight' (if one of `width' or
+  // `height' is zero, natural scaling is used; if both are zero, original image
+  // dimensions in pixels are used):
+  gmsh::view::addListDataString(v, {0, 0.09, 0},
+                                {"file://../t4_image.png@0.01x0"},
+                                {"Align", "Center"});
+
+  // The 3D orientation of the image can be specified by proving the direction
+  // of the bottom and left edge of the image in model space:
+  gmsh::view::addListDataString(v, {-0.01, 0.09, 0},
+                                {"file://../t4_image.png@0.01x0,0,0,1,0,1,0"});
+
+  // The image can also be drawn in "billboard" mode, i.e. always parallel to
+  // the camera, by using the `#' symbol:
+  gmsh::view::addListDataString(v, {0, 0.12, 0},
+                                {"file://../t4_image.png@0.01x0#"},
+                                {"Align", "Center"});
+
+  // The size of 2D annotations is given directly in pixels:
+  gmsh::view::addListDataString(v, {150, -7}, {"file://../t4_image.png@20x0"});
+
+  // These annotations are handled by a list-based post-processing view. For
+  // large post-processing datasets, that contain actual field values defined on
+  // a mesh, you should use model-based post-processing views instead, which
+  // allow to efficiently store continuous or discontinuous scalar, vector and
+  // tensor fields, or arbitrary polynomial order.
+
+  // Views and geometrical entities can be made to respond to double-click
+  // events, here to print some messages to the console:
+  gmsh::option::setString("View[0].DoubleClickedCommand",
+                          "Printf('View[0] has been double-clicked!');");
+  gmsh::option::setString("Geometry.DoubleClickedLineCommand",
+                          "Printf('Curve %g has been double-clicked!', "
+                          "Geometry.DoubleClickedEntityTag);");
+
+  // We can also change the color of some entities:
+  gmsh::model::setColor({{2, 22}}, 127, 127, 127); // Gray50
+  gmsh::model::setColor({{2, 24}}, 160, 32, 240); // Purple
+  for(int i = 1; i <= 14; i++)
+    gmsh::model::setColor({{1, i}}, 255, 0, 0); // Red
+  for(int i = 15; i <= 20; i++)
+    gmsh::model::setColor({{1, i}}, 255, 255, 0); // Yellow
+
   model::mesh::generate(2);
 
   gmsh::write("t4.msh");
 
+  // Launch the GUI to see the results:
+  // gmsh::fltk::run();
+
   gmsh::finalize();
   return 0;
 }
diff --git a/tutorial/c++/t5.cpp b/tutorial/c++/t5.cpp
index c515c403598ca3e18f4bec2817e612735dfc2214..2cede4f98b17630b230808d0c1a86de478319030 100644
--- a/tutorial/c++/t5.cpp
+++ b/tutorial/c++/t5.cpp
@@ -1,7 +1,10 @@
-// This file reimplements gmsh/tutorial/t5.geo in C++.
-
-// The same geometry created with the OpenCASCADE CAD kernel (instead of the
-// built-in kernel) is avaiable in t16.cpp.
+/*******************************************************************************
+ *
+ *  Gmsh C++ tutorial 5
+ *
+ *  Characteristic lengths, holes in volumes
+ *
+ *******************************************************************************/
 
 #include <gmsh.h>
 #include <cstdio>
@@ -12,7 +15,9 @@ namespace factory = gmsh::model::geo;
 void cheeseHole(double x, double y, double z, double r, double lc,
                 std::vector<int> &shells, std::vector<int> &volumes)
 {
-  // When the tag is not specified, a new one is automatically provided
+  // This function will create a spherical hole in a volume. We don't specify
+  // tags manually, and let the functions return them automatically:
+
   int p1 = factory::addPoint(x,  y,  z,  lc);
   int p2 = factory::addPoint(x+r,y,  z,   lc);
   int p3 = factory::addPoint(x,  y+r,z,   lc);
@@ -43,6 +48,17 @@ void cheeseHole(double x, double y, double z, double r, double lc,
   int l7 = factory::addCurveLoop({-c2,-c7,-c12});
   int l8 = factory::addCurveLoop({-c6,-c9,c2});
 
+  // We need non-plane surfaces to define the spherical holes. Here we use the
+  // `gmsh::model::geo::addSurfaceFilling()' function, which can be used for
+  // surfaces with 3 or 4 curves on their boundary. With the he built-in kernel,
+  // if the curves are circle arcs, ruled surfaces are created; otherwise
+  // transfinite interpolation is used.
+  //
+  // With the OpenCASCADE kernel, `gmsh::model::occ::addSurfaceFilling()' uses a
+  // much more general generic surface filling algorithm, creating a BSpline
+  // surface passing through an arbitrary number of boundary curves. The
+  // `gmsh::model::geo::addThruSections()' allows to create ruled surfaces.
+
   int s1 = factory::addSurfaceFilling({l1});
   int s2 = factory::addSurfaceFilling({l2});
   int s3 = factory::addSurfaceFilling({l3});
@@ -60,13 +76,36 @@ void cheeseHole(double x, double y, double z, double r, double lc,
 
 int main(int argc, char **argv)
 {
-  gmsh::initialize();
+  gmsh::initialize(argc, argv);
   gmsh::option::setNumber("General.Terminal", 1);
 
   double lcar1 = .1;
   double lcar2 = .0005;
   double lcar3 = .055;
 
+  // If we wanted to change these mesh sizes globally (without changing the
+  // above definitions), we could give a global scaling factor for all
+  // characteristic lengths with e.g.
+  //
+  // gmsh::option::setNumber("Mesh.CharacteristicLengthFactor", 0.1);
+  //
+  // Since we pass `argc' and `argv' to `gmsh::initialize()', we can also give
+  // the option on the command line with the `-clscale' switch. For example,
+  // with:
+  //
+  // > ./t5.exe -clscale 1
+  //
+  // this tutorial produces a mesh of approximately 3000 nodes and 14,000
+  // tetrahedra. With
+  //
+  // > ./t5.exe -clscale 0.2
+  //
+  // the mesh counts approximately 231,000 nodes and 1,360,000 tetrahedra. You
+  // can check mesh statistics in the graphical user interface
+  // (gmsh::fltk::run()) with the `Tools->Statistics' menu.
+
+  // We proceed by defining some elementary entities describing a truncated cube:
+
   factory::addPoint(0.5,0.5,0.5, lcar2, 1);
   factory::addPoint(0.5,0.5,0, lcar1, 2);
   factory::addPoint(0,0.5,0.5, lcar1, 3);
@@ -128,6 +167,7 @@ int main(int argc, char **argv)
   int sl = factory::addSurfaceLoop({35,31,29,37,33,23,39,25,27});
   shells.push_back(sl);
 
+  // We create five holes in the cube:
   double x = 0, y = 0.75, z = 0, r = 0.09 ;
   for(int t = 1; t <= 5; t++){
     x += 0.166 ;
@@ -138,12 +178,52 @@ int main(int argc, char **argv)
                 t, x, y, z, r, volumes.back());
   }
 
-  factory::addVolume(shells, 186);
+  // The volume of the cube, without the 5 holes, is defined by 6 surface loops:
+  // the first surface loop defines the exterior surface; the surface loops
+  // other than the first one define holes:
+  int ve = factory::addVolume(shells);
+
+  // Note that using solid modelling with the OpenCASCADE geometry kernel, the
+  // same geometry could be built quite differently: see `t16.cpp'.
+
+  // We finally define a physical volume for the elements discretizing the cube,
+  // without the holes (for which physical groups were already defined in the
+  // `cheeseHole()' function):
+  model::addPhysicalGroup(3, {ve}, 10);
 
-  model::addPhysicalGroup(3, {186}, 10);
   factory::synchronize();
+
+  // We could make only part of the model visible to only mesh this subset:
+  // std::vector<std::pair<int, int> > ent;
+  // gmsh::model::getEntities(ent);
+  // gmsh::model::setVisibility(ent, false);
+  // gmsh::model::setVisibility({{3, 5}}, true, true);
+  // gmsh::option::setNumber("Mesh.MeshOnlyVisible", 1);
+
+  // Meshing algorithms can changed globally using options:
+  gmsh::option::setNumber("Mesh.Algorithm", 6); // 2D algorithm set to Frontal-Delaunay
+  gmsh::option::setNumber("Mesh.Algorithm3D", 1); // 3D algorithm set to Delaunay
+
+  // They can also be set for individual surfaces, e.g. for using `MeshAdapt' on
+  // surface 1:
+  gmsh::model::mesh::setAlgorithm(2, 33, 1);
+
+  // To generate a curvilinear mesh and optimize it to produce provably valid
+  // curved elements (see A. Johnen, J.-F. Remacle and C. Geuzaine. Geometric
+  // validity of curvilinear finite elements. Journal of Computational Physics
+  // 233, pp. 359-372, 2013; and T. Toulorge, C. Geuzaine, J.-F. Remacle,
+  // J. Lambrechts. Robust untangling of curvilinear meshes. Journal of
+  // Computational Physics 254, pp. 8-26, 2013), you can uncomment the following
+  // lines:
+  //
+  // gmsh::option::setNumber("Mesh.ElementOrder", 2);
+  // gmsh::option::setNumber("Mesh.HighOrderOptimize", 2);
+
   model::mesh::generate(3);
   gmsh::write("t5.msh");
+
+  // gmsh::fltk::run();
+
   gmsh::finalize();
 
   return 0;
diff --git a/tutorial/c++/t6.cpp b/tutorial/c++/t6.cpp
index 431d66dc93f33ff7065655e35363b145bd6efd75..9461debb7764cb11b75da2001099ddbb847c5039 100644
--- a/tutorial/c++/t6.cpp
+++ b/tutorial/c++/t6.cpp
@@ -1,4 +1,10 @@
-// This file reimplements gmsh/tutorial/t6.geo in C++.
+/*******************************************************************************
+ *
+ *  Gmsh C++ tutorial 6
+ *
+ *  Transfinite meshes
+ *
+ *******************************************************************************/
 
 #include <gmsh.h>
 
@@ -24,54 +30,51 @@ int main(int argc, char **argv)
   factory::addLine(4, 1, 4);
   factory::addCurveLoop({4, 1, -2, 3}, 1);
   factory::addPlaneSurface({1}, 1);
-  model::addPhysicalGroup(0, {1, 2}, 1);
-  model::addPhysicalGroup(1, {1, 2}, 2);
-  model::addPhysicalGroup(2, {1}, 6);
-  model::setPhysicalName(2, 6, "My surface");
-  // ...end of copy
 
-  // Delete surface 1 and left boundary (curve 4)
+  // Delete the surface and the left line, and replace the line with 3 new ones:
   factory::remove({{2,1}, {1,4}});
 
-  // Replace left boundary with 3 new lines
   int p1 = factory::addPoint(-0.05, 0.05, 0, lc);
   int p2 = factory::addPoint(-0.05, 0.1, 0, lc);
   int l1 = factory::addLine(1, p1);
   int l2 = factory::addLine(p1, p2);
   int l3 = factory::addLine(p2, 4);
 
-  // Recreate surface
+  // Create surface:
   factory::addCurveLoop({2, -1, l1, l2, l3, -3}, 2);
   factory::addPlaneSurface({-2}, 1);
 
-  // Put 20 points with a refinement toward the extremities on curve 2
+  // The `setTransfiniteCurve()' meshing constraints explicitly specifies the
+  // location of the nodes on the curve. For example, the following command
+  // forces 20 uniformly placed nodes on curve 2 (including the nodes on the two
+  // end points):
   factory::mesh::setTransfiniteCurve(2, 20, "Bump", 0.05);
 
-  // Put 20 points total on combination of curves l1, l2 and l3 (beware that the
-  // points p1 and p2 are shared by the curves, so we do not create 6 + 6 + 10 =
-  // 22 points, but 20!)
+  // Let's put 20 points total on combination of curves `l1', `l2' and `l3'
+  // (beware that the points `p1' and `p2' are shared by the curves, so we do
+  // not create 6 + 6 + 10 = 22 nodes, but 20!)
   factory::mesh::setTransfiniteCurve(l1, 6);
   factory::mesh::setTransfiniteCurve(l2, 6);
   factory::mesh::setTransfiniteCurve(l3, 10);
 
-  // Put 30 points following a geometric progression on curve 1 (reversed) and
-  // on curve 3
+  // Finally, we put 30 nodes following a geometric progression on curve 1
+  // (reversed) and on curve 3: Put 30 points following a geometric progression
   factory::mesh::setTransfiniteCurve(1, 30, "Progression", -1.2);
   factory::mesh::setTransfiniteCurve(3, 30, "Progression", 1.2);
 
-  // Define the Surface as transfinite, by specifying the four corners of the
-  // transfinite interpolation
+  // The `setTransfiniteSurface()' meshing constraint uses a transfinite
+  // interpolation algorithm in the parametric plane of the surface to connect
+  // the nodes on the boundary using a structured grid. If the surface has more
+  // than 4 corner points, the corners of the transfinite interpolation have to
+  // be specified by hand:
   factory::mesh::setTransfiniteSurface(1, "Left", {1,2,3,4});
 
-  // Recombine the triangles into quads
+  // To create quadrangles instead of triangles, one can use the `setRecombine'
+  // constraint:
   factory::mesh::setRecombine(2, 1);
 
-  // Apply an elliptic smoother to the grid
-  gmsh::option::setNumber("Mesh.Smoothing", 100);
-  model::addPhysicalGroup(2, {1}, 1);
-
-  // When the surface has only 3 or 4 control points, the transfinite constraint
-  // can be applied automatically (without specifying the corners explictly).
+  // When the surface has only 3 or 4 points on its boundary the list of corners
+  // can be omitted in the `setTransfiniteSurface()' call:
   factory::addPoint(0.2, 0.2, 0, 1.0, 7);
   factory::addPoint(0.2, 0.1, 0, 1.0, 8);
   factory::addPoint(0, 0.3, 0, 1.0, 9);
@@ -87,7 +90,14 @@ int main(int argc, char **argv)
     factory::mesh::setTransfiniteCurve(i, 10);
   factory::mesh::setTransfiniteSurface(15);
 
-  model::addPhysicalGroup(2, {15}, 2);
+  // The way triangles are generated can be controlled by specifying "Left",
+  // "Right" or "Alternate" in `setTransfiniteSurface()' command. Try e.g.
+  //
+  // factory::mesh::setTransfiniteSurface(15, "Alternate");
+
+  // Finally we apply an elliptic smoother to the grid to have a more regular
+  // mesh:
+  gmsh::option::setNumber("Mesh.Smoothing", 100);
 
   factory::synchronize();
   model::mesh::generate(2);
diff --git a/tutorial/c++/t7.cpp b/tutorial/c++/t7.cpp
index e337db595bc80d2dba79f300db3f12c4c55434e4..43914abd9fe392c1062109d2713a8aab07f98968 100644
--- a/tutorial/c++/t7.cpp
+++ b/tutorial/c++/t7.cpp
@@ -1,34 +1,39 @@
-// This file reimplements gmsh/tutorial/t7.geo in C++.
-//
-// Background mesh
+/*******************************************************************************
+ *
+ *  Gmsh C++ tutorial 7
+ *
+ *  Background mesh
+ *
+ *******************************************************************************/
 
 #include <gmsh.h>
 
-namespace model = gmsh::model;
-namespace factory = gmsh::model::geo;
+// Mesh sizes can be specified very accurately by providing a background mesh,
+// i.e., a post-processing view that contains the target characteristic lengths.
 
 int main(int argc, char **argv)
 {
   gmsh::initialize();
   gmsh::option::setNumber("General.Terminal", 1);
 
-  model::add("t7");
+  gmsh::model::add("t7");
 
-  // Copied from t1.cpp...
+  // Create a simple rectangular geometry
   double lc = 1e-2;
-  factory::addPoint(0, 0, 0, lc, 1);
-  factory::addPoint(.1, 0,  0, lc, 2);
-  factory::addPoint(.1, .3, 0, lc, 3);
-  factory::addPoint(0,  .3, 0, lc, 4);
-  factory::addLine(1, 2, 1);
-  factory::addLine(3, 2, 2);
-  factory::addLine(3, 4, 3);
-  factory::addLine(4, 1, 4);
-  factory::addCurveLoop({4, 1, -2, 3}, 1);
-  factory::addPlaneSurface({1}, 1);
-  factory::synchronize();
-
-  // add the background mesh file as a view
+  gmsh::model::geo::addPoint(0, 0, 0, lc, 1);
+  gmsh::model::geo::addPoint(.1, 0,  0, lc, 2);
+  gmsh::model::geo::addPoint(.1, .3, 0, lc, 3);
+  gmsh::model::geo::addPoint(0,  .3, 0, lc, 4);
+  gmsh::model::geo::addLine(1, 2, 1);
+  gmsh::model::geo::addLine(3, 2, 2);
+  gmsh::model::geo::addLine(3, 4, 3);
+  gmsh::model::geo::addLine(4, 1, 4);
+  gmsh::model::geo::addCurveLoop({4, 1, -2, 3}, 1);
+  gmsh::model::geo::addPlaneSurface({1}, 1);
+
+  gmsh::model::geo::synchronize();
+
+  // Merge a post-processing view containing the target mesh sizes
   try {
     gmsh::merge("../t7_bgmesh.pos");
   }
@@ -38,11 +43,13 @@ int main(int argc, char **argv)
     return 0;
   }
 
-  // add the post-processing view as a new size field
-  int bg_field = model::mesh::field::add("PostView");
-  model::mesh::field::setAsBackgroundMesh(bg_field);
+  // Add the post-processing view as a new size field
+  int bg_field = gmsh::model::mesh::field::add("PostView");
 
-  model::mesh::generate(2);
+  // Apply the view as the current background mesh
+  gmsh::model::mesh::field::setAsBackgroundMesh(bg_field);
+
+  gmsh::model::mesh::generate(2);
   gmsh::write("t7.msh");
 
   gmsh::finalize();
diff --git a/tutorial/c++/t8.cpp b/tutorial/c++/t8.cpp
index eabf507216ad903966630a43c82efaff58d7bab1..96944b5faeff170cdc3967e0c6d5e999e914b028 100644
--- a/tutorial/c++/t8.cpp
+++ b/tutorial/c++/t8.cpp
@@ -1,7 +1,16 @@
-// This file reimplements gmsh/tutorial/t8.geo in C++.
+/*******************************************************************************
+ *
+ *  Gmsh C++ tutorial 8
+ *
+ *  Post-processing, animations, options
+ *
+ *******************************************************************************/
 
 #include <gmsh.h>
 
+// In addition to creating geometries and meshes, the C++ API can also be used
+// to manipulate post-processing datasets (called "views" in Gmsh).
+
 namespace model = gmsh::model;
 namespace factory = gmsh::model::geo;
 namespace option = gmsh::option;
@@ -13,7 +22,7 @@ int main(int argc, char **argv)
 
   model::add("t8");
 
-  // Copied from t1.cpp...
+  // We first create a simple geometry
   double lc = 1e-2;
   factory::addPoint(0, 0, 0, lc, 1);
   factory::addPoint(.1, 0,  0, lc, 2);
@@ -28,7 +37,7 @@ int main(int argc, char **argv)
   factory::synchronize();
 
   try {
-    // add post-processing views to work on
+    // We merge some post-processing views to work on
     gmsh::merge("../view1.pos");
     gmsh::merge("../view1.pos");
     gmsh::merge("../view4.pos"); // contains 2 views inside
@@ -39,7 +48,14 @@ int main(int argc, char **argv)
     return 0;
   }
 
-  // set general options
+  // Gmsh can read post-processing views in various formats. Here the
+  // `view1.pos' and `view4.pos' files are in the Gmsh "parsed" format, which is
+  // interpreted by the GEO script parser. The parsed format should only be used
+  // for relatively small datasets of course: for larger datasets using e.g. MSH
+  // files is much more efficient. Post-processing views can also be created
+  // directly from the C++ API.
+
+  // We then set some general options:
   option::setNumber("General.Trackball", 0);
   option::setNumber("General.RotationX", 0);
   option::setNumber("General.RotationY", 0);
@@ -47,9 +63,6 @@ int main(int argc, char **argv)
 
   int white[3] = {255, 255, 255};
   int black[3] = {0, 0, 0};
-
-  // Color options are special: setting a color option of "X.Y" actually sets
-  // the option "X.Color.Y"
   option::setColor("General.Background", white[0], white[1], white[2]);
 
   // We could make our own shorter versions of repetitive methods
@@ -62,10 +75,10 @@ int main(int argc, char **argv)
   option::setNumber("General.Axes", 0);
   option::setNumber("General.SmallAxes", 0);
 
-  // show the GUI
+  // Show the GUI
   gmsh::fltk::initialize();
 
-  // set options for each view
+  // We also set some options for each post-processing view:
   option::setNumber("View[0].IntervalsType", 2);
   option::setNumber("View[0].OffsetZ", 0.05);
   option::setNumber("View[0].RaiseZ", 0);
@@ -74,7 +87,7 @@ int main(int argc, char **argv)
   option::setNumber("View[0].SmoothNormals", 1);
 
   option::setNumber("View[1].IntervalsType", 1);
-  // can't set ColorTable in API yet
+  // We can't yet set the ColorTable through the API
   // option::setColorTable("View[1].ColorTable", "{ Green, Blue }");
   option::setNumber("View[1].NbIso", 10);
   option::setNumber("View[1].ShowScale", 0);
@@ -92,14 +105,24 @@ int main(int argc, char **argv)
 
   option::setNumber("View[3].Visible", 0);
 
-  int t = 0;
+  // You can save an MPEG movie directly by selecting `File->Export' in the
+  // GUI. Several predefined animations are setup, for looping on all the time
+  // steps in views, or for looping between views.
+
+  // But the API can be used to build much more complex animations, by changing
+  // options at run-time and re-rendering the graphics. Each frame can then be
+  // saved to disk as an image, and multiple frames can be encoded to form a
+  // movie. Below is an example of such a custom animation.
+
+  int t = 0; // Initial step
 
-  // step through time for each view
   for(int num = 1; num <= 3; num++) {
 
     double nbt;
     option::getNumber("View[0].NbTimeStep", nbt);
     t = (t < nbt - 1) ? t + 1 : 0;
+
+    // Set time step
     option::setNumber("View[0].TimeStep", t);
     option::setNumber("View[1].TimeStep", t);
     option::setNumber("View[2].TimeStep", t);
@@ -118,6 +141,8 @@ int main(int argc, char **argv)
 
     int frames = 50;
     for(int num2 = 1; num2 <= frames; num2++) {
+
+      // Incrementally rotate the scene
       double rotx;
       option::getNumber("General.RotationX", rotx);
       option::setNumber("General.RotationX", rotx + 10);
@@ -126,11 +151,11 @@ int main(int argc, char **argv)
       option::getNumber("General.RotationZ", rotz);
       option::setNumber("General.RotationZ", rotz + 0.1);
 
+      // Draw the scene
       gmsh::graphics::draw();
 
-      // write out the graphics scene to an image file Gmsh will try to detect
-      // the file extension
       if(num == 3) {
+        // Uncomment the following lines to save each frame to an image file
         // gmsh::write("t2-" + std::to_string(num2) + ".gif");
         // gmsh::write("t2-" + std::to_string(num2) + ".ppm");
         // gmsh::write("t2-" + std::to_string(num2) + ".jpg");
diff --git a/tutorial/c++/t9.cpp b/tutorial/c++/t9.cpp
index 00ed8f88c8499b186a5cad6776bb4a192f771a3b..a2071b86708a23fa751f5b77f3ccc853fd956809 100644
--- a/tutorial/c++/t9.cpp
+++ b/tutorial/c++/t9.cpp
@@ -1,4 +1,19 @@
-// This file reimplements gmsh/tutorial/t9.geo in C++.
+/*******************************************************************************
+ *
+ *  Gmsh C++ tutorial 9
+ *
+ *  Post-processing plugins (levelsets, sections, annotations)
+ *
+ *******************************************************************************/
+
+// Plugins can be added to Gmsh in order to extend its capabilities. For
+// example, post-processing plugins can modify views, or create new views based
+// on previously loaded views. Several default plugins are statically linked
+// with Gmsh, e.g. Isosurface, CutPlane, CutSphere, Skin, Transform or Smooth.
+//
+// Plugins can be controlled through the API functions in the `gmsh::plugin'
+// namespace, or from the graphical interface (right click on the view button,
+// then `Plugins').
 
 #include <gmsh.h>
 
@@ -9,7 +24,7 @@ int main(int argc, char **argv)
 
   gmsh::model::add("t9");
 
-  // add a three-dimensional scalar view to work on:
+  // Let us for example include a three-dimensional scalar view:
   try {
     gmsh::merge("../view3.pos");
   }
@@ -19,12 +34,14 @@ int main(int argc, char **argv)
     return 0;
   }
 
-  // First plugin is Isosurface
-  gmsh::plugin::setNumber("Isosurface", "Value", 0.67);
-  gmsh::plugin::setNumber("Isosurface", "View", 0);
-  gmsh::plugin::run("Isosurface");
+  // We then set some options for the `Isosurface' plugin (which extracts an
+  // isosurface from a 3D scalar view), and run it:
+  gmsh::plugin::setNumber("Isosurface", "Value", 0.67); // Iso-value level
+  gmsh::plugin::setNumber("Isosurface", "View", 0); // Source view is View[0]
+  gmsh::plugin::run("Isosurface"); // Run the plugin!
 
-  // Second is CutPlane
+  // We also set some options for the `CutPlane' plugin (which computes a
+  // section of a 3D view using the plane A*x+B*y+C*z+D=0), and then run it:
   gmsh::plugin::setNumber("CutPlane", "A", 0);
   gmsh::plugin::setNumber("CutPlane", "B", 0.2);
   gmsh::plugin::setNumber("CutPlane", "C", 1);
@@ -32,9 +49,10 @@ int main(int argc, char **argv)
   gmsh::plugin::setNumber("CutPlane", "View", 0);
   gmsh::plugin::run("CutPlane");
 
-  // Third is Annotate
+  // Add a title (By convention, for window coordinates a value greater than
+  // 99999 represents the center. We could also use `General.GraphicsWidth / 2',
+  // but that would only center the string for the current window size.):
   gmsh::plugin::setString("Annotate", "Text", "A nice title");
-  // By convention, window coordinates larger than 99999 represent the center
   gmsh::plugin::setNumber("Annotate", "X", 1.e5);
   gmsh::plugin::setNumber("Annotate", "Y", 50);
   gmsh::plugin::setString("Annotate", "Font", "Times-BoldItalic");
@@ -49,7 +67,7 @@ int main(int argc, char **argv)
   gmsh::plugin::setNumber("Annotate", "FontSize", 12);
   gmsh::plugin::run("Annotate");
 
-  // set some general options
+  // We finish by setting some options:
   gmsh::option::setNumber("View[0].Light", 1);
   gmsh::option::setNumber("View[0].IntervalsType", 1);
   gmsh::option::setNumber("View[0].NbIso", 6);
diff --git a/tutorial/c++/x1.cpp b/tutorial/c++/x1.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c44ffe4486ddb48f8ae87170077fe357a43b892c
--- /dev/null
+++ b/tutorial/c++/x1.cpp
@@ -0,0 +1,113 @@
+/*******************************************************************************
+ *
+ *  Gmsh C++ extended tutorial 1
+ *
+ *  Accessing geometrical and mesh data
+ *
+ *******************************************************************************/
+
+// The C++ API allows to do much more than what can be done in .geo files. These
+// additional features are introduced gradually in the extended tutorials,
+// starting with `x1.cpp'.
+
+// In this first extended tutorial, we start by using the API to explore basic
+// geometrical and mesh data.
+
+#include <iostream>
+#include <gmsh.h>
+
+int main(int argc, char **argv)
+{
+  if(argc < 2){
+    std::cout << "Usage: " << argv[0] << " file" << std::endl;
+    return 0;
+  }
+
+  gmsh::initialize();
+  gmsh::option::setNumber("General.Terminal", 1);
+  gmsh::open(argv[1]);
+
+  // Geometrical data is made of elementary model `entities', called `points'
+  // (entities of dimension 0), `curves' (entities of dimension 1), `surfaces'
+  // (entities of dimension 2) and `volumes' (entities of dimension 3). As we
+  // have seen in the other C++ tutorials, elementary model entities are
+  // identified by their dimension and by a `tag': a strictly positive
+  // identification number. Model entities can be either CAD entities (from the
+  // built-in `geo' kernel or from the OpenCASCADE `occ' kernel) or `discrete'
+  // entities (defined by a mesh).
+
+  // Get all the entities in the model, as a vector of (dimension, tag) pairs:
+  std::vector<std::pair<int, int> > entities;
+  gmsh::model::getEntities(entities);
+
+  for(unsigned int i = 0; i < entities.size(); i++){
+    // Mesh data is made of `elements' (points, lines, triangles, ...), defined
+    // by an ordered list of their `nodes'. Elements and nodes are identified by
+    // `tags' as well (strictly positive identification numbers), and are stored
+    // ("classified") in the model entity they discretize.
+
+    // A model entity of dimension 0 (a geometrical point) will contain a mesh
+    // element of type point, as well as a mesh node. A model curve will contain
+    // line elements as well as its interior nodes, while its boundary nodes
+    // will be stored in the bounding model points. A model surface will contain
+    // triangular and/or quadrangular elements and all the nodes not classified
+    // on its boundary or on its embedded entities (curves and points). A model
+    // volume will contain tetrahedra, hexahedra, etc. and all the nodes not
+    // classified on its boundary or on its embedded entities (surfaces, curves
+    // and points).
+
+    // Get the mesh nodes for each entity:
+    std::vector<std::size_t> nodeTags;
+    std::vector<double> nodeCoords, nodeParams;
+    int dim = entities[i].first, tag = entities[i].second;
+    gmsh::model::mesh::getNodes(nodeTags, nodeCoords, nodeParams, dim, tag);
+
+    // Get the mesh elements for each entity:
+    std::vector<int> elemTypes;
+    std::vector<std::vector<std::size_t> > elemTags, elemNodeTags;
+    gmsh::model::mesh::getElements(elemTypes, elemTags, elemNodeTags, dim, tag);
+
+    // Print a summary of information available on the entity and its mesh:
+    int numElem = 0;
+    for(unsigned int i = 0; i < elemTags.size(); i++)
+      numElem += elemTags[i].size();
+    std::string type;
+    gmsh::model::getType(dim, tag, type);
+    std::cout << nodeTags.size() << " mesh nodes and "
+              << numElem << " mesh elements on entity ("
+              << dim << "," << tag << ") of type " << type << "\n";
+    std::vector<int> physicalTags;
+    gmsh::model::getPhysicalGroupsForEntity(dim, tag, physicalTags);
+    if(physicalTags.size()) {
+      std::cout << " - Physical group tag(s):";
+      for(unsigned int i = 0; i < physicalTags.size(); i++)
+        std::cout << " " << physicalTags[i];
+      std::cout << "\n";
+    }
+    std::vector<int> partitions;
+    gmsh::model::getPartitions(dim, tag, partitions);
+    if(partitions.size()){
+      std::cout << " - Partition tag(s):";
+      for(unsigned int i = 0; i < partitions.size(); i++)
+        std::cout << " " << partitions[i];
+      int parentDim, parentTag;
+      gmsh::model::getParent(dim, tag, parentDim, parentTag);
+      std::cout << " - parent entity (" << parentDim << "," << parentTag << ")\n";
+    }
+    for(unsigned int i = 0; i < elemTypes.size(); i++){
+      std::string name;
+      int d, order, numv, numpv;
+      std::vector<double> param;
+      gmsh::model::mesh::getElementProperties(elemTypes[i], name, d, order,
+                                              numv, param, numpv);
+      std::cout << " - Element type: " << name << ", order " << order << "\n";
+      std::cout << "   with " << numv << " nodes in param coord: (";
+      for(unsigned int j = 0; j < param.size(); j++)
+        std::cout << param[j] << " ";
+      std::cout << ")\n";
+    }
+  }
+
+  gmsh::finalize();
+  return 0;
+}
diff --git a/tutorial/c/t1.c b/tutorial/c/t1.c
index d26e6bc43bf33c94012e2da370afa58575daf15d..4be9ce9d13c069f0c28fa6ee2ee00df3ba6c4ce8 100644
--- a/tutorial/c/t1.c
+++ b/tutorial/c/t1.c
@@ -1,83 +1,128 @@
-/* This file reimplements gmsh/tutorial/t1.geo in C.
-
-   For all the elementary explanations about the general philosphy of entities
-   in Gmsh, see the comments in the .geo file. Comments here focus on the
-   specifics of the C API.
-
-   The Gmsh API is entirely defined in the <gmsh.h> header. Read this file: it
-   contains the documentation for all the functions in the API.
-*/
+/*******************************************************************************
+ *
+ *  Gmsh C++ tutorial 1
+ *
+ *  Elementary entities (points, curves, surfaces), physical groups (points,
+ *  curves, surfaces)
+ *
+ *******************************************************************************/
+
+/* The Gmsh C API is entirely defined in the `gmshc.h' header (which contains
+   the full documentation of all the functions in the API): */
 #include <gmshc.h>
 
 int main(int argc, char **argv)
 {
-  int ierr;
+  int ierr, p4;
   double lc = 1e-2;
-  int cl1[] = {4, 1, -2, 3}, ps1[] = {1};
-  int pg1[] = {1, 2}, pg2[] = {1,2}, pg6[] = {1};
+  int cl1[] = {4, 1, -2, 3}, s1[] = {1};
+  int g5[] = {1, 2, 4}, g6[] = {1}, ps;
 
   /* Before using any functions in the C API, Gmsh must be initialized. In the C
-     API the last argument of all functions return the error code, if any. */
+     API the last argument of all functions returns the error code, if any. */
   gmshInitialize(argc, argv, 1, &ierr);
 
   /* By default Gmsh will not print out any messages: in order to output
-     messages on the terminal, just set the standard Gmsh option
-     "General.Terminal" (same format and meaning as in .geo files) using
-     gmshOptionSetNumber() */
+     messages on the terminal, just set the "General.Terminal" option to 1: */
   gmshOptionSetNumber("General.Terminal", 1, &ierr);
 
-  /* This adds a new model, named "t1". If gmshModelAdd() is not called, a
+  /* We now add a new model, named "t1". If gmsh::model::add() is not called, a
      new default (unnamed) model will be created on the fly, if necessary. */
   gmshModelAdd("t1", &ierr);
 
-  /* The C API provides direct access to the internal CAD kernels. The
-     built-in CAD kernel was used in t1.geo: the corresponding API functions
-     live in the "gmshModelGeo" namespace. To create geometrical points with
-     the built-in CAD kernel, one thus uses gmshModelGeoAddPoint():
+  /* The C API provides direct access to each supported geometry kernel. The
+     built-in kernel is used in this first tutorial: the corresponding API
+     functions have the "gmshModelGeo" prefix.
 
+     The first type of `elementary entity' in Gmsh is a `Point'. To create a
+     point with the built-in geometry kernel, the C API function is
+     gmshModelGeoAddPoint():
      - the first 3 arguments are the point coordinates (x, y, z)
+     - the next argument is the target mesh size (the "characteristic length")
+       close to the point
+     - the last argument is the point tag (a stricly positive integer that
+       uniquely identifies the point); if the tag is set to -1, the function
+       will return a new tag */
+  gmshModelGeoAddPoint(0, 0, 0, lc, 1, &ierr);
 
-     - the next argument is the target mesh size close to the point
+  /* The distribution of the mesh element sizes will be obtained by
+     interpolation of these characteristic lengths throughout the
+     geometry. Another method to specify characteristic lengths is to use
+     general mesh size Fields. A particular case is the use of a background
+     mesh.
 
-     - the next argument is the point tag */
-  gmshModelGeoAddPoint(0, 0, 0, lc, 1, &ierr);
+     If no target mesh size of provided, a default uniform coarse size will be
+     used for the model, based on the overall model size.
+
+     We can then define some additional points. All points should have different
+     tags: */
   gmshModelGeoAddPoint(.1, 0,  0, lc, 2, &ierr);
   gmshModelGeoAddPoint(.1, .3, 0, lc, 3, &ierr);
-  gmshModelGeoAddPoint(0,  .3, 0, lc, 4, &ierr);
 
-  /* The API to create lines with the built-in kernel follows the same
-     conventions: the first 2 arguments are point tags, the last (optional one)
-     is the line tag. */
+  /* If the tag is not provided explicitly, a new tag is automatically created,
+     and returned by the function: */
+  p4 = gmshModelGeoAddPoint(0,  .3, 0, lc, -1, &ierr);
+
+  /* Curves are Gmsh's second type of elementery entities, and, amongst curves,
+     straight lines are the simplest. The API to create straight line segments
+     with the built-in kernel follows the same conventions: the first 2
+     arguments are point tags (the start and end points of the line), and the
+     last is the line tag.
+
+     In the commands below, for example, the line 1 starts at point 1 and ends
+     at point 2.
+
+     Note that curve tags are separate from point tags - hence we can reuse tag
+     `1' for our first curve. And as a general rule, elementary entity tags in
+     Gmsh have to be unique per geometrical dimension. */
   gmshModelGeoAddLine(1, 2, 1, &ierr);
   gmshModelGeoAddLine(3, 2, 2, &ierr);
-  gmshModelGeoAddLine(3, 4, 3, &ierr);
-  gmshModelGeoAddLine(4, 1, 4, &ierr);
-
-  /* The philosophy to construct curve loops and surfaces is similar: the first
-     argument is now a vector of integers. */
+  gmshModelGeoAddLine(3, p4, 3, &ierr);
+  gmshModelGeoAddLine(4, 1, p4, &ierr);
+
+  /* The third elementary entity is the surface. In order to define a simple
+     rectangular surface from the four curves defined above, a curve loop has
+     first to be defined. A curve loop is defined by an ordered list of
+     connected curves, a sign being associated with each curve (depending on the
+     orientation of the curve to form a loop). The API function to create curve
+     loops takes an array of integers as first argument, and the curve loop tag
+     (which must ne unique amongst curve loops) as the second argument: */
   gmshModelGeoAddCurveLoop(cl1, sizeof(cl1)/sizeof(cl1[0]), 1, &ierr);
-  gmshModelGeoAddPlaneSurface(ps1, sizeof(ps1)/sizeof(ps1[0]), 1, &ierr);
 
-  /* Physical groups are defined by providing the dimension of the group (0 for
-     physical points, 1 for physical curves, 2 for physical surfaces and 3 for
-     phsyical volumes) followed by a vector of entity tags. The last (optional)
-     argument is the tag of the new group to create. */
-  gmshModelAddPhysicalGroup(0, pg1, sizeof(pg1)/sizeof(pg1[0]), 1, &ierr);
-  gmshModelAddPhysicalGroup(1, pg2, sizeof(pg2)/sizeof(pg2[0]), 2, &ierr);
-  gmshModelAddPhysicalGroup(2, pg6, sizeof(pg6)/sizeof(pg6[0]), 6, &ierr);
-
-  /* Physical names are also defined by providing the dimension and tag of the
-     entity. */
-  gmshModelSetPhysicalName(2, 6, "My surface", &ierr);
+  /* We can then define the surface as a list of curve loops (only one here,
+     representing the external contour, since there are no holes): */
+  gmshModelGeoAddPlaneSurface(s1, sizeof(s1)/sizeof(s1[0]), 1, &ierr);
+
+  /* At this level, Gmsh knows everything to display the rectangular surface 6
+     and to mesh it. An optional step is needed if we want to group elementary
+     geometrical entities into more meaningful groups, e.g. to define some
+     mathematical ("domain", "boundary"), functional ("left wing", "fuselage")
+     or material ("steel", "carbon") properties.
+
+     Such groups are called "Physical Groups" in Gmsh. By default, if physical
+     groups are defined, Gmsh will export in output files only mesh elements
+     that belong to at least one physical group. (To force Gmsh to save all
+     elements, whether they belong to physical groups or not, set the
+     `Mesh.SaveAll' option to 1.) Physical groups are also identified by tags,
+     i.e. stricly positive integers, that should be unique per dimension (0D,
+     1D, 2D or 3D). Physical groups can also be given names.
+
+     Here we define a physical curve that groups the left, bottom and right
+     curves in a single group (with prescribed tag 5); and a physical surface
+     with name "My surface" (with an automatic tag) containing the geometrical
+     surface 1: */
+  gmshModelAddPhysicalGroup(1, g5, sizeof(g5)/sizeof(g5[0]), 5, &ierr);
+  ps = gmshModelAddPhysicalGroup(2, g6, sizeof(g6)/sizeof(g6[0]), -1, &ierr);
+  gmshModelSetPhysicalName(2, ps, "My surface", &ierr);
 
   /* Before it can be meshed, the internal CAD representation must be
      synchronized with the Gmsh model, which will create the relevant Gmsh data
-     structures. This is achieved by the gmsh::model::geo::synchronize() API
-     call for the built-in CAD kernel. Synchronizations can be called at any
+     structures. This is achieved by the `gmshModelGeoSynchronize()' API call
+     for the built-in geometry kernel. Synchronizations can be called at any
      time, but they involve a non trivial amount of processing; so while you
      could synchronize the internal CAD data after every CAD command, it is
      usually better to minimize the number of synchronization points. */
-  gmshModelGeoSynchronize( &ierr);
+  gmshModelGeoSynchronize(&ierr);
 
   /* We can then generate a 2D mesh... */
   gmshModelMeshGenerate(2, &ierr);
@@ -92,7 +137,29 @@ int main(int argc, char **argv)
      gmshOptionSetNumber("Mesh.SaveAll", 1);
   */
 
-  /* This should be called at the end: */
+  /* We could run the graphical user interface with:
+     gmsh::fltk::run();
+  */
+
+  /* Note that starting with Gmsh 3.0, models can be built using other geometry
+     kernels than the default "built-in" kernel. To use the OpenCASCADE geometry
+     kernel instead of the built-in kernel, you should use the functions with the
+     "gmshModelOcc" prefix.
+
+     Different geometry kernels have different features. With OpenCASCADE,
+     instead of defining the surface by successively defining 4 points, 4 curves
+     and 1 curve loop, one can define the rectangular surface directly with
+
+     gmshModelOccAddRectangle(.2, 0, 0, .1, .3, -1, 0, &ierr);
+
+     Afther synchronization with the Gmsh model with
+
+     gmshModelOccSynchronize();
+
+     the underlying curves and points could be accessed with
+     gmshModelGetBoundary(). */
+
+  /* This should be called when you are done using the Gmsh C API: */
   gmshFinalize(&ierr);
   return 0;
 }
diff --git a/tutorial/t1.geo b/tutorial/t1.geo
index 8c82d24634a633a3958a641c7b8fb9a99a6ad291..69061da7a0ba7bfd680fcf2eba07e138391ba986 100644
--- a/tutorial/t1.geo
+++ b/tutorial/t1.geo
@@ -115,5 +115,5 @@ Physical Surface("My surface") = {1};
 // The underlying curves and points could be accessed with the `Boundary' or
 // `CombinedBoundary' operators.
 //
-// See `t16.geo', `t17.geo' or `t18.geo' for complete examples based on
+// See e.g. `t16.geo', `t18.geo' or `t19.geo' for complete examples based on
 // OpenCASCADE, and `demos/boolean' for more.
diff --git a/tutorial/t11.geo b/tutorial/t11.geo
index 5b998a1042973721b2b6b6411caf613161bd6254..73f05f6db54525ea82cf00ea109f3ff7634b4bcb 100644
--- a/tutorial/t11.geo
+++ b/tutorial/t11.geo
@@ -24,7 +24,7 @@ Field[1].F = "0.01*(1.0+30.*(y-x*x)*(y-x*x) + (1-x)*(1-x))";
 Background Field = 1;
 
 // To generate quadrangles instead of triangles, we can simply add
-Recombine Surface{100};
+//Recombine Surface{100};
 
 // If we'd had several surfaces, we could have used `Recombine Surface {:};'.
 // Yet another way would be to specify the global option "Mesh.RecombineAll =
@@ -58,3 +58,11 @@ Recombine Surface{100};
 // the full-quad algorithm:
 //
 // Mesh.RecombinationAlgorithm = 2; // or 3
+
+// Note that you could also apply the recombination algorithm and/or the
+// subdivision step explicitly after meshing, as follows:
+//
+// Mesh 2;
+// RecombineMesh;
+// Mesh.SubdivisionAlgorithm = 1;
+// RefineMesh;
diff --git a/tutorial/t16.geo b/tutorial/t16.geo
index 1991322e639a6134b1eade5b491e6ff48e331bc9..596667cf92549626329c079a486b2cc2422d3326 100644
--- a/tutorial/t16.geo
+++ b/tutorial/t16.geo
@@ -8,11 +8,11 @@
 
 // Instead of constructing a model in a bottom-up fashion with Gmsh's built-in
 // geometry kernel, starting with version 3 Gmsh allows you to directly use
-// alternative geometry kernels. Let us use the OpenCASCADE kernel:
+// alternative geometry kernels. Here we use the OpenCASCADE kernel:
 
 SetFactory("OpenCASCADE");
 
-// And let's build the same model as in `t5.geo', but using constructive solid
+// Let's build the same model as in `t5.geo', but using constructive solid
 // geometry.
 
 // We first create two cubes:
@@ -26,7 +26,7 @@ BooleanDifference(3) = { Volume{1}; Delete; }{ Volume{2}; Delete; };
 // `Delete' in the arguments allows to automatically delete the original
 // entities.
 
-// We then create the five spheres
+// We then create the five spheres:
 x = 0 ; y = 0.75 ; z = 0 ; r = 0.09 ;
 For t In {1:5}
   x += 0.166 ;
@@ -42,7 +42,7 @@ EndFor
 v() = BooleanFragments{ Volume{3}; Delete; }{ Volume{3 + 1 : 3 + 5}; Delete; };
 
 // When the boolean operation leads to simple modifications of entities, and if
-// one deletes the original entities with `Delete', Gmsh tries to assign the 
+// one deletes the original entities with `Delete', Gmsh tries to assign the
 // same tag to the new entities. (This behavior is governed by the
 // `Geometry.OCCBooleanPreserveNumbering' option.)
 
@@ -79,4 +79,5 @@ p() = Point In BoundingBox{0.5-eps, 0.5-eps, 0.5-eps,
 Characteristic Length{ p() } = lcar2;
 
 // Additional examples created with the OpenCASCADE geometry kernel are
-// available in the `demos/boolean' directory.
+// available in `t18.geo', `t19.geo' as well as in the `demos/boolean'
+// directory.
diff --git a/tutorial/t17.geo b/tutorial/t17.geo
index 6f02728e3950ee359b3217ca15376e240a7bb0d2..c892dc9501a020fa449125c24527f23e15826191 100644
--- a/tutorial/t17.geo
+++ b/tutorial/t17.geo
@@ -16,7 +16,7 @@
 
 SetFactory("OpenCASCADE");
 
-// the square
+// Create a square
 Rectangle(1) = {-1, -1, 0, 2, 2};
 
 // Merge a post-processing view containing the target anisotropic mesh sizes
diff --git a/tutorial/t18.geo b/tutorial/t18.geo
index 50c94fe418700584f5aa2dd464086ebe22a7e468..02a454fc6e6e38d469de8f9b448eb45a0850dbb0 100644
--- a/tutorial/t18.geo
+++ b/tutorial/t18.geo
@@ -70,25 +70,23 @@ Characteristic Length {p()} = 0.001;
 // geometry automatically.
 
 // First we get all surfaces on the left:
-e = 1e-3;
-Sxmin() = Surface In BoundingBox{2-e, -e, -e, 2+e, 1+e, 1+e};
+Sxmin() = Surface In BoundingBox{2-eps, -eps, -eps, 2+eps, 1+eps, 1+eps};
 
 For i In {0:#Sxmin()-1}
   // Then we get the bounding box of each left surface
   bb() = BoundingBox Surface { Sxmin(i) };
-  // We translate the bounding box to the right...
-  bbe() = { bb(0)-e+1, bb(1)-e, bb(2)-e, bb(3)+e+1, bb(4)+e, bb(5)+e };
-  // ...and look for surfaces inside it
-  Sxmax() = Surface In BoundingBox { bbe() };
+  // We translate the bounding box to the right and look for surfaces inside i
+  Sxmax() = Surface In BoundingBox { bb(0)-eps+1, bb(1)-eps, bb(2)-eps,
+                                     bb(3)+eps+1, bb(4)+eps, bb(5)+eps };
   // For all the matches, we compare the corresponding bounding boxes...
   For j In {0:#Sxmax()-1}
     bb2() = BoundingBox Surface { Sxmax(j) };
     bb2(0) -= 1;
     bb2(3) -= 1;
     // ...and if they match, we apply the periodicity constraint
-    If(Fabs(bb2(0)-bb(0)) < e && Fabs(bb2(1)-bb(1)) < e &&
-       Fabs(bb2(2)-bb(2)) < e && Fabs(bb2(3)-bb(3)) < e &&
-       Fabs(bb2(4)-bb(4)) < e && Fabs(bb2(5)-bb(5)) < e)
+    If(Fabs(bb2(0)-bb(0)) < eps && Fabs(bb2(1)-bb(1)) < eps &&
+       Fabs(bb2(2)-bb(2)) < eps && Fabs(bb2(3)-bb(3)) < eps &&
+       Fabs(bb2(4)-bb(4)) < eps && Fabs(bb2(5)-bb(5)) < eps)
       Periodic Surface {Sxmax(j)} = {Sxmin(i)} Translate {1,0,0};
     EndIf
   EndFor
diff --git a/tutorial/t3.geo b/tutorial/t3.geo
index d443cc473e38cff4424c30c9acc6e4671d3719ee..baaac2b54fb90818c18f5e4c91e747df6c59038c 100644
--- a/tutorial/t3.geo
+++ b/tutorial/t3.geo
@@ -68,7 +68,7 @@ out[] = Extrude { {-2*h,0,0}, {1,0,0} , {0,0.15,0.25} , angle * Pi / 180 } {
 // programmatically by using the return value (a list) of the `Extrude'
 // command. This list contains the "top" of the extruded surface (in `out[0]'),
 // the newly created volume (in `out[1]') and the tags of the lateral surfaces
-// (in `out[2]', `out[3]', ...)
+// (in `out[2]', `out[3]', ...).
 
 // We can then define a new physical volume (with tag 101) to group all the
 // elementary volumes:
diff --git a/tutorial/t4.geo b/tutorial/t4.geo
index 5c57c15b73c1f8860b0b8a602290f5fc010c2096..595e88c4391c28a0b8e77b13520f6d1dde13711e 100644
--- a/tutorial/t4.geo
+++ b/tutorial/t4.geo
@@ -2,7 +2,7 @@
  *
  *  Gmsh GEO tutorial 4
  *
- *  Built-in functions, surface holes, annotations, mesh colors
+ *  Built-in functions, holes in surfaces, annotations, entity colors
  *
  *******************************************************************************/
 
@@ -121,7 +121,7 @@ View[0].DoubleClickedCommand = "Printf('View[0] has been double-clicked!');";
 Geometry.DoubleClickedLineCommand = "Printf('Curve %g has been double-clicked!',
   Geometry.DoubleClickedEntityTag);";
 
-// We can also change the color of some mesh entities:
+// We can also change the color of some entities:
 
 Color Grey50{ Surface{ 22 }; }
 Color Purple{ Surface{ 24 }; }
diff --git a/tutorial/t5.geo b/tutorial/t5.geo
index 4b1e3d8c8c5608b8fc8274eb91e2cf3634c61195..afb839969a0e833ffc4419325224120b7ef9b919 100644
--- a/tutorial/t5.geo
+++ b/tutorial/t5.geo
@@ -2,7 +2,7 @@
  *
  *  Gmsh GEO tutorial 5
  *
- *  Characteristic lengths, arrays of variables, macros, loops
+ *  Characteristic lengths, macros, loops, holes in volumes
  *
  *******************************************************************************/
 
@@ -169,8 +169,8 @@ Volume(186) = {theloops[]};
 // same geometry could be built quite differently: see `t16.geo'.
 
 // We finally define a physical volume for the elements discretizing the cube,
-// without the holes (whose elements were already tagged with numbers 1 to 5 in
-// the `For' loop):
+// without the holes (for which physical groups were already created in the
+// `For' loop):
 
 Physical Volume (10) = 186;
 
diff --git a/tutorial/t6.geo b/tutorial/t6.geo
index 76bbdad7c401c1e9cb2ff4d190d44de5d2e9a1ae..9accfee7ed2d211523bd2a57f31611052399c426 100644
--- a/tutorial/t6.geo
+++ b/tutorial/t6.geo
@@ -9,7 +9,7 @@
 // Let's use the geometry from the first tutorial as a basis for this one:
 Include "t1.geo";
 
-// Delete the left line and replace it with 3 new ones:
+// Delete the surface and the left line, and replace the line with 3 new ones:
 Delete{ Surface{1}; Curve{4}; }
 
 p1 = newp; Point(p1) = {-0.05, 0.05, 0, lc};