diff --git a/Surface_mesh_topology/examples/Surface_mesh_topology/path_homotopy_with_symbols_2.cpp b/Surface_mesh_topology/examples/Surface_mesh_topology/path_homotopy_with_symbols_2.cpp
index 5c8ef1f09fc..2d6fef267f5 100644
--- a/Surface_mesh_topology/examples/Surface_mesh_topology/path_homotopy_with_symbols_2.cpp
+++ b/Surface_mesh_topology/examples/Surface_mesh_topology/path_homotopy_with_symbols_2.cpp
@@ -13,12 +13,12 @@ int main()
   ps.add_facet("a b -a c"); // First facet, giving directly its sequence of edges
   ps.add_facet("d -c e -b"); // Second facet
 
-  ps.begin_facet(); // Third facet
+  ps.init_facet(); // Third facet
   ps.add_edges_to_facet("f"); // Here, each edge is added one at a time
   ps.add_edges_to_facet("-d");
   ps.add_edges_to_facet("-f");
   ps.add_edges_to_facet("-e");
-  ps.end_facet();
+  ps.finish_facet();
 
   ps.perforate_facet("f");
   
diff --git a/Surface_mesh_topology/include/CGAL/Path_on_surface.h b/Surface_mesh_topology/include/CGAL/Path_on_surface.h
index 24fa8a946c7..cd21fd2d13c 100644
--- a/Surface_mesh_topology/include/CGAL/Path_on_surface.h
+++ b/Surface_mesh_topology/include/CGAL/Path_on_surface.h
@@ -37,13 +37,16 @@
 #include <initializer_list>
 
 // A Path_on_surface contains two vectors of equal length n
-// The first one is a vector of darts called m_path and the second one a vector of booleans called m_flip
-// If n = 0, the path represented by those vectors is the empty path
+// The first one is a vector of darts called m_path and the second one a vector
+// of booleans called m_flip.
+// If n = 0, the path represented by those vectors is the empty path.
 // Else, it is the path represented by the n-1 first elements of both vectors,
-// at the one we add the m_path[n-1] dart if m_flip[n-1] is false and the opposite of this dart if m_flip[n-1] is true
-// i.e. if m_flip[i] is true means that the i-th dart m_path[i] has to be flipped
-// We use flips because sometimes opposite darts doesn't exist on surfaces with boundaries
-// But if m_flip[i] is true doesn't necesary mean that m_path[i] is 2-free
+// at the one we add the m_path[n-1] dart if m_flip[n-1] is false and the
+// opposite of this dart if m_flip[n-1] is true i.e. if m_flip[i] is true means
+// that the i-th dart m_path[i] has to be flipped.
+// We use flips because sometimes opposite darts doesn't exist on surfaces with
+// boundaries. But if m_flip[i] is true doesn't necesary mean that
+// m_path[i] is 2-free
 
 namespace CGAL {
 namespace Surface_mesh_topology {
@@ -722,13 +725,19 @@ public:
   /// If this face does not exist (if it is a boundary) then replace the edge
   /// by the face on the other side. Problem of complexity when used many times
   /// (like in update_path_randomly).
-  void push_around_face(std::size_t i, bool update_isclosed=true)
+  bool push_around_face(std::size_t i, bool update_isclosed=true)
   {
     CGAL_assertion(i<length());
+
+    // It is not possible to push around a perforated face since it changes
+    // the homotopy of the path.
+    if (get_map().is_perforated(get_ith_dart(i))) { return false; }
+
     Self p2(get_mesh());
 
+    // 1) We add in p2 the part of the path which is pushed.
     if (get_ith_flip(i))
-    {// in this case the face of the ith dart doesn't exist
+    {
       Dart_const_handle dh=get_map().template beta<1>(get_ith_dart(i));
       do
       {
@@ -748,16 +757,20 @@ public:
       while(dh!=get_ith_dart(i));
     }
 
+    // 2) We copy the end of the path.
     p2.m_path.reserve(p2.length()+length()-i);
     for (std::size_t j=i+1; j<length(); ++j)
     { p2.push_back(get_ith_dart(j), get_ith_flip(j), false); }
 
+    // 3) We cut this path to keep the first i darts.
     cut(i, false);
     m_path.reserve(length()+p2.length());
     for (std::size_t j=0; j<p2.length(); ++j)
     { push_back(p2[j], p2.get_ith_flip(j), false); }
 
     if (update_isclosed) { update_is_closed(); }
+    return true;
+
     //CGAL_assertion(is_valid());
   }
 
@@ -768,8 +781,8 @@ public:
   {
     if (is_empty()) return;
 
-    for (unsigned int i=0; i<nb; ++i)
-    { push_around_face(random.get_int(0, length()), false); }
+    for (unsigned int i=0; i<nb; )
+    { if (push_around_face(random.get_int(0, length()), false)) { ++i; } }
     if (update_isclosed) { update_is_closed(); }
   }
 
diff --git a/Surface_mesh_topology/include/CGAL/Polygonal_schema.h b/Surface_mesh_topology/include/CGAL/Polygonal_schema.h
index 597cfff64ca..4fb9b1eb0f4 100644
--- a/Surface_mesh_topology/include/CGAL/Polygonal_schema.h
+++ b/Surface_mesh_topology/include/CGAL/Polygonal_schema.h
@@ -28,11 +28,14 @@
 #include <cstddef>
 #include <string>
 #include <initializer_list>
+#include <algorithm>
+#include <random>
 #include <CGAL/assertions.h>
 #include <CGAL/memory.h>
 #include <CGAL/Polygonal_schema_min_items.h>
 #include <CGAL/Combinatorial_map.h>
 #include <CGAL/Generalized_map.h>
+#include <CGAL/Random.h>
 
 namespace CGAL {
 namespace Surface_mesh_topology {
@@ -223,7 +226,7 @@ namespace Surface_mesh_topology {
     }
     
     /// Start a new facet.
-    void begin_facet()
+    void init_facet()
     {
       if (facet_started)
       {
@@ -291,9 +294,9 @@ namespace Surface_mesh_topology {
                  <<" but the previous facet is not yet ended."<<std::endl;
         return;
       }
-      begin_facet();
+      init_facet();
       add_edges_to_facet(s);
-      end_facet();
+      finish_facet();
     }
 
     /// add edges to the current facet,
@@ -321,13 +324,13 @@ namespace Surface_mesh_topology {
                  <<" but the previous facet is not yet ended."<<std::endl;
         return;
       }
-      begin_facet();
+      init_facet();
       add_edges_to_facet(l);
-      end_facet();
+      finish_facet();
     }
     
     /// End of the facet. Return the first dart of this facet.
-    Dart_handle end_facet()
+    Dart_handle finish_facet()
     {
       if (!facet_started)
       {
@@ -429,39 +432,6 @@ namespace Surface_mesh_topology {
     std::size_t get_perforated_mark() const
     { return mark_perforated; }
 
-    using Base::is_free;
-
-    /// same thing but using a label instead of a dart
-    template<unsigned int i>
-    bool is_free(const std::string& s) const
-    {
-      auto ite=edge_label_to_dart.find(s);
-      if (ite==edge_label_to_dart.end())
-      {// maybe there is no need to put an error message
-        std::cerr<<"Polygonal_schema ERROR: "
-                 <<"you ask if label "<<s<<" represents a dart border"
-                 <<" but this label does not exist yet."<<std::endl;
-        return false;
-      }
-
-      return Self::template is_free<i>(ite->second);
-    }
-
-    /// Non templated versions
-    bool is_free(const std::string& s, unsigned int i) const
-    {
-      auto ite=edge_label_to_dart.find(s);
-      if (ite==edge_label_to_dart.end())
-      {// maybe there is no need to put an error message
-        std::cerr<<"Polygonal_schema ERROR: "
-                 <<"you ask if label "<<s<<" represents a dart border"
-                 <<" but this label does not exist yet."<<std::endl;
-        return false;
-      }
-
-      return is_free(ite->second, i);
-    }
-
     void display_perforated_darts() const
     {
       std::cout<<"labels is_free<2> is_perforated"<<std::endl;
@@ -474,8 +444,8 @@ namespace Surface_mesh_topology {
     }
     
   protected:
-    // For each edge label, its corresponding dart. Stores both association a -a, to allow
-    // users to start to add either a or -a.
+    // For each edge label, its corresponding dart. Stores both association
+    // a -a, to allow users to start to add either a or -a.
     std::unordered_map<std::string, Dart_handle> edge_label_to_dart;
     std::size_t mark_perforated; // mark for perforated facets.
 
@@ -485,6 +455,7 @@ namespace Surface_mesh_topology {
     bool        facet_started;
   };
 
+  /// Polygonal schema with combinatorial map.
   template <class Items_=Polygonal_schema_min_items,
             class Alloc_=CGAL_ALLOCATOR(int),
             class Storage_= Combinatorial_map_storage_1<2, Items_, Alloc_> >
@@ -545,6 +516,7 @@ namespace Surface_mesh_topology {
     {}
   };
   
+  /// Polygonal schema with generalized map.
   template <class Items_=Polygonal_schema_min_items,
             class Alloc_=CGAL_ALLOCATOR(int),
             class Storage_= Generalized_map_storage_1<2, Items_, Alloc_> >
@@ -605,6 +577,53 @@ namespace Surface_mesh_topology {
     {}
   };
   
+  /// Generate a random polygonal schema ps.
+  /// @param nb_labels the number of labels used to generate ps.
+  /// @param seed the seed used for random
+  /// @param max_dart_per_face maximal number of darts per face
+  /// @param closed if true generates a closed polygonal schema.
+  /// @param percentage_of_perforated percentage of perforated faces. If 0
+  ///         no perforated faces.
+  template<typename PS>
+  void generate_random_polygonal_schema(PS& ps, std::size_t nb_labels,
+                                        unsigned int seed
+                                        =static_cast<unsigned int>(std::time(nullptr)),
+                                        std::size_t max_dart_per_face=30,
+                                        bool closed=true,
+                                        std::size_t percentage_of_perforated=20)
+  {
+    CGAL::Random random(seed);
+    std::vector<std::string> all_labels(nb_labels*2);
+    for (std::size_t i=0; i<nb_labels; ++i)
+    {
+      all_labels[2*i]=std::to_string(static_cast<int>(i)+1);
+      all_labels[(2*i)+1]=std::to_string(-(static_cast<int>(i)+1));
+    }
+
+    std::shuffle(all_labels.begin(), all_labels.end(),
+                 std::default_random_engine(seed));
+
+    std::size_t endlabel=all_labels.size();
+    if (!closed)
+    { endlabel-=(all_labels.size()/10); } // We remove 10% of labels.
+
+    for (std::size_t i=0; i<endlabel; )
+    {
+      ps.init_facet();
+      for (std::size_t j=0,
+           nb=static_cast<std::size_t>
+           (random.get_int(1, static_cast<int>(max_dart_per_face)));
+           i<endlabel && j<nb; ++i, ++j)
+      { ps.add_edges_to_facet(all_labels[i]); }
+      typename PS::Dart_handle dh=ps.finish_facet();
+
+      if (rand()%101<percentage_of_perforated)
+      { ps.perforate_facet(dh); }
+    }
+
+    ps.keep_biggest_connected_component(); // We keep only the biggest cc.
+  }
+
 } //namespace Surface_mesh_topology
 } //namespace CGAL
 
diff --git a/Surface_mesh_topology/include/CGAL/Surface_mesh_topology/internal/Minimal_quadrangulation.h b/Surface_mesh_topology/include/CGAL/Surface_mesh_topology/internal/Minimal_quadrangulation.h
index daf57cf76b1..149c69f52bb 100644
--- a/Surface_mesh_topology/include/CGAL/Surface_mesh_topology/internal/Minimal_quadrangulation.h
+++ b/Surface_mesh_topology/include/CGAL/Surface_mesh_topology/internal/Minimal_quadrangulation.h
@@ -697,7 +697,8 @@ protected:
           if (get_original_map().is_perforated(it))
           { get_reduced_map().mark(d1, m_mark_perforated); }
         }
-        else if (it<get_original_map().template beta<2>(it))
+        else if (Original_dart_const_handle(it)<
+                 get_original_map().template beta<2>(it))
         {
           d1=get_reduced_map().create_dart();
           d2=get_reduced_map().create_dart();
diff --git a/Surface_mesh_topology/test/Surface_mesh_topology/test_homotopy_with_polygonal_schema.cpp b/Surface_mesh_topology/test/Surface_mesh_topology/test_homotopy_with_polygonal_schema.cpp
index be36118bb48..612a7273db0 100644
--- a/Surface_mesh_topology/test/Surface_mesh_topology/test_homotopy_with_polygonal_schema.cpp
+++ b/Surface_mesh_topology/test/Surface_mesh_topology/test_homotopy_with_polygonal_schema.cpp
@@ -2,8 +2,6 @@
 #include <CGAL/Path_on_surface.h>
 #include <CGAL/Curves_on_surface_topology.h>
 #include <iostream>
-#include <vector>
-#include <algorithm>
 #include <ctime>
 #include <cstdlib>
 
@@ -14,55 +12,21 @@ using namespace CGAL::Surface_mesh_topology;
 
 static unsigned int starting_seed;
 
-///////////////////////////////////////////////////////////////////////////////
-int myrandom (int i) { return std::rand()%i;}
-//-----------------------------------------------------------------------------
-void generate_random_map(PS& ps, std::size_t nb_labels,
-                         std::size_t max_dart_per_face=30,
-                         bool closed=true,
-                         std::size_t percentage_of_perforated=20) // 0 => no perforated faces
-{
-  std::vector<std::string> all_labels(nb_labels*2);
-  for (std::size_t i=0; i<nb_labels; ++i)
-  {
-    all_labels[2*i]=std::to_string(static_cast<int>(i)+1);
-    all_labels[(2*i)+1]=std::to_string(-(static_cast<int>(i)+1));
-  }
-
-  std::random_shuffle(all_labels.begin(), all_labels.end(), myrandom);
-
-  std::size_t endlabel=all_labels.size();
-  if (!closed)
-  { endlabel-=(all_labels.size()/10); } // We remove 10% of labels.
-
-  for (std::size_t i=0; i<endlabel; )
-  {
-    ps.begin_facet();
-    for (std::size_t j=0,
-         nb=1+(static_cast<std::size_t>(rand())%max_dart_per_face);
-         i<endlabel && j<nb; ++i, ++j)
-    {
-      ps.add_edges_to_facet(all_labels[i]);
-    }
-    Dart_handle dh=ps.end_facet();
-
-    if (rand()%101<percentage_of_perforated)
-    { ps.perforate_facet(dh); }
-  }
-
-  ps.keep_biggest_connected_component();
-}
 ///////////////////////////////////////////////////////////////////////////////
 bool test_two_random_paths(const PS& ps,
                            const Curves_on_surface_topology<PS>& cst,
                            unsigned int nbtests,
-                           unsigned int lmin=5,
-                           unsigned int lmax=20)
+                           unsigned int& seed,
+                           int lmin=5,
+                           int lmax=20)
 {
   CGAL_assertion(lmin>0 && lmin<lmax);
-  CGAL::Random random(starting_seed+nbtests);
+  CGAL::Random random(seed++);
   Path_on_surface<PS> p1(ps), p2(ps);
-  internal::generate_random_closed_path(p1, random.get_int(lmin, lmax), random);
+  internal::generate_random_closed_path(p1,
+                                        static_cast<std::size_t>
+                                        (random.get_int(lmin, lmax)),
+                                        random);
 
   p2=p1;
   p2.update_path_randomly(100, random);
@@ -72,8 +36,8 @@ bool test_two_random_paths(const PS& ps,
   if (!cst.are_freely_homotopic(p1, p2))
   {
     res=false;
-    std::cout<<"ERROR for path number "<<nbtests<<" (seed="<<random.get_seed()
-             <<")"<<std::endl; res=false;
+    std::cout<<std::endl<<"ERROR for path number "<<nbtests<<" (seed="
+             <<random.get_seed()<<")"<<std::endl;
   }
   
   return res;
@@ -82,21 +46,21 @@ bool test_two_random_paths(const PS& ps,
 bool run_n_random_paths_tests(const PS& ps,
                               unsigned int n,
                               const Curves_on_surface_topology<PS>& cst,
-                              unsigned int lmin=5,
-                              unsigned int lmax=20)
+                              unsigned int& seed,
+                              int lmin=5,
+                              int lmax=20)
 {
   static std::size_t nbtest=0;
-  std::cout<<"Test "<<nbtest++<<std::flush;
+  std::cout<<"***** Test "<<nbtest++<<" "<<std::flush;
+  ps.display_characteristics(std::cout)<<"   ";
   bool res=true;
   for (unsigned int i=0; i<n; ++i)
   {
     std::cout<<"."<<std::flush;
-    if (!test_two_random_paths(ps, cst, i, lmin, lmax))
+    if (!test_two_random_paths(ps, cst, i, seed, lmin, lmax))
     { res=false; }
   }
-  if (res) { std::cout<<"OK."; }
-   else { std::cout<<"FAILED."; }
-  std::cout<<std::endl;
+  std::cout<<(res?" OK.":" FAILED.")<<std::endl;
   return res;
 }
 ///////////////////////////////////////////////////////////////////////////////
@@ -105,7 +69,8 @@ bool run_n_random_paths_tests(const PS& ps,
   std::cout<<"usage: "<<argv[0]<<" [-seed S]"<<std::endl
            <<"   Run different homotopy on 2-maps built using polygonal schema."
            <<std::endl
-           <<"   -seed S: uses S as seed of random generator. Otherwise use a different seed at each run (based on time)."<<std::endl
+           <<"   -seed S: uses S as seed of random generator. Otherwise use a"
+           <<" different seed at each run (based on time)."<<std::endl
            <<std::endl;
   exit(EXIT_FAILURE);
 }
@@ -141,27 +106,28 @@ void process_command_line(int argc, char** argv,
 int main(int argc, char** argv)
 {  
   process_command_line(argc, argv, starting_seed);
+  std::cout<<"Start tests="<<starting_seed<<"."<<std::endl;
 
-  std::cout<<"Start tests (seed="<<starting_seed<<")"<<std::endl;
-  srand(starting_seed);
-  CGAL::Random r(starting_seed);
-
-  PS ps1, ps2, ps3, ps4;
-  generate_random_map(ps1, 5000, 30, true, 0);   // Closed, no perforated faces
-  generate_random_map(ps2, 5000, 30, true, 15);  // Closed, 15% of perforated faces
-  generate_random_map(ps3, 5000, 30, false, 0);  // Open, no perforated faces
-  generate_random_map(ps4, 5000, 30, false, 15); // Open, 15% of perforated faces
-
-  Curves_on_surface_topology<PS> cst1(ps1);
-  Curves_on_surface_topology<PS> cst2(ps2);
-  Curves_on_surface_topology<PS> cst3(ps3);
-  Curves_on_surface_topology<PS> cst4(ps4);
-
+  bool closed=true;
+  std::size_t percent=0;
   bool res=true;
-  //res&=run_n_random_paths_tests(ps1, 20, cst1);
-  //res&=run_n_random_paths_tests(ps2, 20, cst2);
-  //res&=run_n_random_paths_tests(ps3, 20, cst3);
-  res&=run_n_random_paths_tests(ps4, 20, cst4);
+  for (int i=0; i<3; ++i)
+  {
+    for (int j=0; j<5; ++j)
+    { // 5 tests for the same "type" of polygonal schema
+      //    (open/closed; with  or without perforated)
+      PS ps;
+      Curves_on_surface_topology<PS> cst(ps);
+      generate_random_polygonal_schema(ps, 5000, starting_seed++,  30, closed,  percent);   // Closed, no perforated faces
+      res&=run_n_random_paths_tests(ps, 20, cst, starting_seed);
+    }
+    switch(i)
+    {
+      case 0: closed=true; percent=15; break;
+      case 1: closed=false; percent=0; break;
+      case 2: closed=false; percent=15; break;
+    }
+  }
 
   if (res)
   { std::cout<<"test_homotopy_with_polygonal_schema ALL TESTS OK."<<std::endl; }