diff --git a/tutorial/t13.geo b/tutorial/t13.geo
index 5f9bb7097ab3ad8c077cd77258ec8d1dfc70844c..16f2acb6eca27b884b0e9750b5cc4bd7d453e90e 100644
--- a/tutorial/t13.geo
+++ b/tutorial/t13.geo
@@ -1,19 +1,51 @@
-cl1 = 0.01;
-Point(1) = {0, 0, 0, 0.01};
-Point(2) = {0.1, 0, 0, 0.01};
-Point(3) = {0.1, 0.3, 0, 0.01};
-Point(4) = {0, 0.3, 0, 0.01};
-Line(1) = {1, 2};
-Line(2) = {3, 2};
-Line(3) = {3, 4};
-Line(4) = {4, 1};
-Line Loop(6) = {4, 1, -2, 3};
-Plane Surface(6) = {6};
-Physical Point(1) = {1, 2};
-MyLine = 99;
-Physical Line(MyLine) = {1, 2, 4};
-Physical Surface("My fancy surface label") = {6};
-Field[1] = PostView;
-Field[1].CropNegativeValues = 1;
-Field[1].IView = 0;
+/*********************************************************************
+ *
+ *  Gmsh tutorial 13
+ *
+ *  Remeshing STL with compounds
+ *
+ *********************************************************************/
+
+// Since compound geometrical compute a new parametrization, one can
+// also use them to remesh STL files, even if in this case there's
+// usually only a single elementary geometrical entity per compound.
+
+// Let's merge the mesh that we would like to remesh. This mesh was
+// reclassified ("colored") from an initial STL triangulation using
+// the "Reclassify 2D" tool in Gmsh, so that we could split it along
+// sharp geometrical features.
+Merge "t13_data.msh";
+
+// Since the original mesh is a bit coarse, we refine it once
+RefineMesh;
+
+// Create the topology of the discrete model
+CreateTopology;
+
+// We can now define a compound line (resp. surface) for each discrete
+// line (resp. surface) in the model
+ll[] = Line "*";
+For j In {0 : #ll[]-1}
+  Compound Line(newl) = ll[j];
+EndFor
+ss[] = Surface "*";
+s = news;
+For i In {0 : #ss[]-1}
+  Compound Surface(s+i) = ss[i];
+EndFor
+
+// And we can create the volume based on the new compound entities
+Surface Loop(1) = {s : s + #ss[]-1};
+Volume(1) = {1};
+
+Physical Surface(1) = {s : s + #ss[]-1};
+Physical Volume(1) = 1;
+
+// Apply a funny mesh size field, just because we can :-)
+Field[1] = MathEval;
+Field[1].F = "2*Sin((x+y)/5) + 3";
 Background Field = 1;
+
+Mesh.RemeshAlgorithm = 1; // (0) no split (1) automatic (2) automatic only with metis
+Mesh.RemeshParametrization = 7; // (0) harmonic (1) conformal spectral (7) conformal finite element
+Geometry.HideCompounds = 0; // don't hide the compound entities