diff --git a/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/clip.h b/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/clip.h
index 288220056162..2e7176b2636a 100644
--- a/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/clip.h
+++ b/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/clip.h
@@ -988,7 +988,7 @@ clip(TriangleMesh& tm,
   *   \cgalParamNEnd
   *
   *    \cgalParamNBegin{use_convex_specialization}
-  *      \cgalParamDescription{If set to `true`, a faster implementation specialized for convex meshes is used. The input mesh must be convex to guarantee a correct execution and results.}
+  *      \cgalParamDescription{If set to `true`, a faster implementation specialized for convex meshes is used. The input mesh must be convex and have no flat vertices to guarantee a correct execution and results.}
   *      \cgalParamType{Boolean}
   *      \cgalParamDefault{`false`}
   *    \cgalParamNEnd
diff --git a/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/internal/clip_convex.h b/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/internal/clip_convex.h
index c7d667ed6e76..7e25228a889e 100644
--- a/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/internal/clip_convex.h
+++ b/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/internal/clip_convex.h
@@ -33,6 +33,29 @@ namespace CGAL {
 namespace Polygon_mesh_processing {
 namespace internal {
 
+template<class Visitor, class PolygonMesh>
+struct Wrap_visitor_for_triangulation: public Triangulate_faces::Default_visitor<PolygonMesh>
+{
+  Wrap_visitor_for_triangulation(const PolygonMesh &pm, Visitor &v): pm(pm), visitor(v){}
+
+  typedef typename boost::graph_traits<PolygonMesh>::face_descriptor face_descriptor;
+  typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
+
+  void before_subface_creations(face_descriptor f_old) {
+    visitor.before_subface_creations(f_old, pm);
+  }
+  void after_subface_creations() {
+    visitor.after_subface_creations(pm);
+  }
+  void after_subface_created(face_descriptor f_new) {
+    visitor.after_subface_created(f_new, pm);
+  }
+
+private:
+  const PolygonMesh &pm;
+  Visitor &visitor;
+};
+
 /**
   * Given a convex mesh and a plane, return an halfedge crossing the plane from a vertex on positive side to a vertex on the plane or the negative side.
   * Return a null halfedge if the mesh and the plane do not intersect.
@@ -145,6 +168,157 @@ find_crossing_edge(PolygonMesh& pm,
   return halfedge(src, trg, pm).first;
 }
 
+/**
+  * Given a convex mesh, a plane and an halfedge crossing the plane from positive side, refine the mesh with the plane
+  * Go in the opposite direction, use only for input with boundaries
+  */
+template <class PolygonMesh, class Plane_3, class NamedParameters =  parameters::Default_named_parameters>
+std::vector<typename boost::graph_traits<PolygonMesh>::halfedge_descriptor>
+refine_convex_with_plane_using_opposite_direction(PolygonMesh& pm,
+                                                  const Plane_3& plane,
+                                                  typename boost::graph_traits<PolygonMesh>::halfedge_descriptor h,
+                                                  const NamedParameters& np = parameters::default_values())
+{
+  using parameters::choose_parameter;
+  using parameters::get_parameter;
+  using parameters::get_parameter_reference;
+
+  // graph typedefs
+  using BGT = boost::graph_traits<PolygonMesh>;
+  using face_descriptor = typename BGT::face_descriptor;
+  using halfedge_descriptor = typename BGT::halfedge_descriptor;
+  using vertex_descriptor = typename BGT::vertex_descriptor;
+
+  // np typedefs
+  using Default_visitor = Corefinement::Default_visitor<PolygonMesh>;
+  using Visitor_ref = typename internal_np::Lookup_named_param_def<internal_np::visitor_t, NamedParameters, Default_visitor>::reference;
+  using GT = typename GetGeomTraits<PolygonMesh, NamedParameters>::type;
+  GT traits = choose_parameter<GT>(get_parameter(np, internal_np::geom_traits));
+
+  using Default_f2f_map = boost::associative_property_map<std::map<face_descriptor, face_descriptor>>;
+  constexpr bool is_f2f_map = !parameters::is_default_parameter<NamedParameters, internal_np::face_to_face_map_t>::value;
+  auto f2f = choose_parameter<Default_f2f_map>(get_parameter(np, internal_np::face_to_face_map));
+
+  Default_visitor default_visitor;
+  Visitor_ref visitor = choose_parameter(get_parameter_reference(np, internal_np::visitor), default_visitor);
+  // constexpr bool has_visitor = !std::is_same_v<Default_visitor, std::remove_cv_t<std::remove_reference_t<Visitor_ref>>>;
+
+  auto vpm = choose_parameter(get_parameter(np, internal_np::vertex_point),
+                              get_property_map(vertex_point, pm));
+
+  // config flags
+  bool triangulate = !choose_parameter(get_parameter(np, internal_np::do_not_triangulate_faces), true);
+
+  auto oriented_side = traits.oriented_side_3_object();
+  auto intersection_point = traits.construct_plane_line_intersection_point_3_object();
+
+  // ____ Cut the convex along the plane by marching along crossing edges starting from the previous edge _____
+  std::vector<halfedge_descriptor> boundaries;
+
+  if(oriented_side(plane, get(vpm, target(h,pm)))!=ON_ORIENTED_BOUNDARY)
+  {
+    // split the first edge
+    auto pts = make_sorted_pair(get(vpm, source(h, pm)), get(vpm, target(h, pm)));
+    typename GT::Point_3 ip = intersection_point(plane, pts.first, pts.second);
+    visitor.before_edge_split(h, pm);
+    h = CGAL::Euler::split_edge(h, pm);
+    put(vpm, target(h, pm), ip);
+    visitor.new_vertex_added(std::numeric_limits<std::size_t>::max(), target(h,pm), pm);
+    visitor.edge_split(h, pm);
+    visitor.after_edge_split();
+  }
+
+  vertex_descriptor v_start = target(h, pm);
+  halfedge_descriptor h_start=h;
+  do{
+    halfedge_descriptor h_previous = h;
+    CGAL_assertion(oriented_side(plane, get(vpm, source(h,pm)))==ON_NEGATIVE_SIDE);
+    CGAL_assertion(oriented_side(plane, get(vpm, target(h,pm)))==ON_ORIENTED_BOUNDARY);
+
+    h = next(h, pm);
+    Oriented_side side_trg = oriented_side(plane, get(vpm, target(h,pm)));
+
+    if(side_trg != ON_POSITIVE_SIDE){
+      // The face does not cross the plane, we look for a face incident to the same vertex than h that crosses the plane.
+      while(side_trg == ON_ORIENTED_BOUNDARY){
+        // The edge is along the plane, add it to boundaries
+        boundaries.emplace_back(h);
+        set_halfedge(target(h, pm), h, pm);
+        h = next(h, pm);
+        side_trg=oriented_side(plane, get(vpm, target(h,pm)));
+      }
+      // continue on next face
+      h = opposite(h, pm);
+
+      // if h is a border, end the function
+      if(is_border(h, pm))
+        return boundaries;
+      continue;
+    }
+
+    // Search a crossing edge
+    h = next(h, pm);
+    side_trg=oriented_side(plane, get(vpm, target(h,pm)));
+    while(side_trg == ON_POSITIVE_SIDE){
+      h = next(h,pm);
+      side_trg=oriented_side(plane, get(vpm, target(h,pm)));
+    }
+
+    if(side_trg != ON_ORIENTED_BOUNDARY){
+      // Split the edge
+      auto pts = make_sorted_pair(get(vpm, source(h,pm)), get(vpm, target(h,pm)));
+      typename GT::Point_3 ip = intersection_point(plane, pts.first, pts.second);
+      visitor.before_edge_split(h, pm);
+      h = CGAL::Euler::split_edge(h, pm);
+      put(vpm, target(h, pm), ip);
+      visitor.new_vertex_added(std::numeric_limits<std::size_t>::max(), target(h,pm), pm);
+      visitor.edge_split(h, pm);
+      visitor.after_edge_split();
+    }
+
+    // Split the face
+    visitor.before_subface_created(pm);
+    visitor.before_subface_creations(face(h, pm), pm);
+    halfedge_descriptor sh = CGAL::Euler::split_face(h_previous, h, pm);
+    if constexpr(is_f2f_map)
+      put(f2f, face(opposite(sh, pm), pm), get(f2f, face(sh, pm)));
+
+    visitor.after_subface_created(face(h, pm), pm);
+    visitor.after_subface_creations(pm);
+    boundaries.emplace_back(sh);
+    set_halfedge(target(sh, pm), sh, pm);
+    visitor.add_retriangulation_edge(sh, pm);
+
+    // The input is assumed to be triangulated
+    if (triangulate){
+      halfedge_descriptor sh_opp = opposite(sh, pm);
+      if(!is_triangle(sh_opp, pm)){
+        visitor.before_subface_created(pm);
+        visitor.before_subface_creations(face(sh_opp, pm), pm);
+        halfedge_descriptor newh = CGAL::Euler::split_face(sh_opp, next(next(sh_opp, pm), pm), pm);
+        if constexpr(is_f2f_map)
+          put(f2f, face(opposite(newh, pm), pm), get(f2f, face(newh, pm)));
+        visitor.after_subface_created(face(opposite(newh, pm), pm), pm);
+        visitor.after_subface_creations(pm);
+        visitor.add_retriangulation_edge(newh, pm);
+      }
+    }
+
+    CGAL_assertion(target(sh, pm) == target(h, pm));
+    h = opposite(next(sh,pm), pm);
+
+    // if h is a border, end the function
+    if(is_border(h, pm))
+      return boundaries;
+
+    // During the loop, if a mesh vertex lies on the plane, we look for a face incident to that vertex that crosses the plane.
+    // The second part of the while-condition ensures we don't exit prematurely
+  } while(target(h, pm)!=v_start || (boundaries.empty() && h!=h_start));
+
+  CGAL_assertion(is_valid_polygon_mesh(pm));
+  return boundaries;
+}
+
 /**
   * Given a convex mesh, a plane and an halfedge crossing the plane from positive side, refine the mesh with the plane
   */
@@ -166,7 +340,6 @@ refine_convex_with_plane(PolygonMesh& pm,
   using vertex_descriptor = typename BGT::vertex_descriptor;
 
   // np typedefs
-  // using Default_ecm = Static_boolean_property_map<edge_descriptor, false>;
   using Default_visitor = Corefinement::Default_visitor<PolygonMesh>;
   using Visitor_ref = typename internal_np::Lookup_named_param_def<internal_np::visitor_t, NamedParameters, Default_visitor>::reference;
   using GT = typename GetGeomTraits<PolygonMesh, NamedParameters>::type;
@@ -194,13 +367,13 @@ refine_convex_with_plane(PolygonMesh& pm,
 
   if(oriented_side(plane, get(vpm, target(h,pm)))!=ON_ORIENTED_BOUNDARY)
   {
-    //split the first edge
+    // split the first edge
     auto pts = make_sorted_pair(get(vpm, source(h, pm)), get(vpm, target(h, pm)));
     typename GT::Point_3 ip = intersection_point(plane, pts.first, pts.second);
     visitor.before_edge_split(h, pm);
     h = CGAL::Euler::split_edge(h, pm);
     put(vpm, target(h, pm), ip);
-    // visitor.new_vertex_added(vpm.size()-1, target(h,pm), pm);
+    visitor.new_vertex_added(std::numeric_limits<std::size_t>::max(), target(h,pm), pm);
     visitor.edge_split(h, pm);
     visitor.after_edge_split();
   }
@@ -226,6 +399,14 @@ refine_convex_with_plane(PolygonMesh& pm,
       }
       // continue on next face
       h = opposite(h, pm);
+
+      // Specific case with boundaries, we explore in the opposite direction
+      if(is_border(h, pm)){
+        auto other_direction_boundaries = refine_convex_with_plane_using_opposite_direction(pm, plane, opposite(h_start, pm), np);
+        std::reverse(other_direction_boundaries.begin(), other_direction_boundaries.end());
+        other_direction_boundaries.insert(other_direction_boundaries.begin(), boundaries.begin(), boundaries.end());
+        return boundaries;
+      }
       continue;
     }
 
@@ -244,28 +425,35 @@ refine_convex_with_plane(PolygonMesh& pm,
       visitor.before_edge_split(h, pm);
       h = CGAL::Euler::split_edge(h, pm);
       put(vpm, target(h, pm), ip);
-      // visitor.new_vertex_added(vpm.size()-1, target(h,pm), pm);
+      visitor.new_vertex_added(std::numeric_limits<std::size_t>::max(), target(h,pm), pm);
       visitor.edge_split(h, pm);
       visitor.after_edge_split();
     }
 
     // Split the face
     visitor.before_subface_created(pm);
+    visitor.before_subface_creations(face(h, pm), pm);
     halfedge_descriptor sh = CGAL::Euler::split_face(h_previous, h, pm);
     if constexpr(is_f2f_map)
       put(f2f, face(opposite(sh, pm), pm), get(f2f, face(sh, pm)));
 
     visitor.after_subface_created(face(h, pm), pm);
+    visitor.after_subface_creations(pm);
     boundaries.emplace_back(sh);
     set_halfedge(target(sh, pm), sh, pm);
     visitor.add_retriangulation_edge(sh, pm);
 
+    // The input is assumed to be triangulated
     if (triangulate){
       halfedge_descriptor sh_opp = opposite(sh, pm);
       if(!is_triangle(sh_opp, pm)){
         visitor.before_subface_created(pm);
+        visitor.before_subface_creations(face(sh_opp, pm), pm);
         halfedge_descriptor newh = CGAL::Euler::split_face(sh_opp, next(next(sh_opp, pm), pm), pm);
+        if constexpr(is_f2f_map)
+          put(f2f, face(opposite(newh, pm), pm), get(f2f, face(newh, pm)));
         visitor.after_subface_created(face(opposite(newh, pm), pm), pm);
+        visitor.after_subface_creations(pm);
         visitor.add_retriangulation_edge(newh, pm);
       }
     }
@@ -273,6 +461,14 @@ refine_convex_with_plane(PolygonMesh& pm,
     CGAL_assertion(target(sh, pm) == target(h, pm));
     h = opposite(next(sh,pm), pm);
 
+    // Specific case with boundaries, we explore in the opposite direction
+    if(is_border(h, pm)){
+      auto other_direction_boundaries = refine_convex_with_plane_using_opposite_direction(pm, plane, opposite(h_start, pm), np);
+      std::reverse(other_direction_boundaries.begin(), other_direction_boundaries.end());
+      other_direction_boundaries.insert(other_direction_boundaries.begin(), boundaries.begin(), boundaries.end());
+      return boundaries;
+    }
+
     // During the loop, if a mesh vertex lies on the plane, we look for a face incident to that vertex that crosses the plane.
     // The second part of the while-condition ensures we don't exit prematurely
   } while(target(h, pm)!=v_start || (boundaries.empty() && h!=h_start));
@@ -295,6 +491,7 @@ remove_bounded_region_and_fill(PolygonMesh& pm,
   using parameters::choose_parameter;
   using parameters::is_default_parameter;
   using parameters::get_parameter;
+  using parameters::get_parameter_reference;
 
   // graph typedefs
   using BGT = boost::graph_traits<PolygonMesh>;
@@ -304,9 +501,8 @@ remove_bounded_region_and_fill(PolygonMesh& pm,
   using vertex_descriptor = typename BGT::vertex_descriptor;
 
   // np typedefs
-  // using Default_ecm = Static_boolean_property_map<edge_descriptor, false>;
-  // using Default_visitor = Corefinement::Default_visitor<PolygonMesh>;
-  // using Visitor_ref = typename internal_np::Lookup_named_param_def<internal_np::visitor_t, NamedParameters, Default_visitor>::reference;
+  using Default_visitor = Corefinement::Default_visitor<PolygonMesh>;
+  using Visitor_ref = typename internal_np::Lookup_named_param_def<internal_np::visitor_t, NamedParameters, Default_visitor>::reference;
   using GT = typename GetGeomTraits<PolygonMesh, NamedParameters>::type;
   // GT traits = choose_parameter<GT>(get_parameter(np, internal_np::geom_traits));
 
@@ -318,9 +514,9 @@ remove_bounded_region_and_fill(PolygonMesh& pm,
   if constexpr (update_bbox)
     bbox_pointer = get_parameter(np, internal_np::bounding_box);
 
-  // Default_visitor default_visitor;
-  // Visitor_ref visitor = choose_parameter(get_parameter_reference(np, internal_np::visitor), default_visitor);
-  // constexpr bool has_visitor = !std::is_same_v<Default_visitor, std::remove_cv_t<std::remove_reference_t<Visitor_ref>>>;
+  Default_visitor default_visitor;
+  Visitor_ref visitor = choose_parameter(get_parameter_reference(np, internal_np::visitor), default_visitor);
+  constexpr bool has_visitor = !std::is_same_v<Default_visitor, std::remove_cv_t<std::remove_reference_t<Visitor_ref>>>;
 
   // Used only if do_triangulate_faces or bounding_box check
   auto vpm = choose_parameter(get_parameter(np, internal_np::vertex_point),
@@ -421,6 +617,7 @@ remove_bounded_region_and_fill(PolygonMesh& pm,
 
   // Fill the hole
   if (clip_volume && faces(pm).size()>1){ // If there is only one face, the output is flat does not need to be closed
+    visitor.before_face_copy(BGT::null_face(), pm, pm);
     face_descriptor f=add_face(pm);
     for(auto h: boundaries){
       set_face(h, f, pm);
@@ -431,8 +628,16 @@ remove_bounded_region_and_fill(PolygonMesh& pm,
       auto f2f = choose_parameter<Default_f2f_map>(get_parameter(np, internal_np::face_to_face_map));
       put(f2f, f, plane_fd);
     }
-    if(triangulate)
-      triangulate_face(f, pm, parameters::vertex_point_map(vpm));
+    visitor.after_face_copy(BGT::null_face(), pm, f, pm);
+
+    if(triangulate){
+      if constexpr(has_visitor){
+        Wrap_visitor_for_triangulation wrap_visitor(pm, visitor);
+        triangulate_face(f, pm, parameters::vertex_point_map(vpm).visitor(wrap_visitor));
+      } else {
+        triangulate_face(f, pm, parameters::vertex_point_map(vpm));
+      }
+    }
   } else {
     for(auto h: boundaries){
       set_face(h, BGT::null_face(), pm);
@@ -483,13 +688,13 @@ clip_convex(PolygonMesh& pm,
 
   Default_visitor default_visitor;
   Visitor_ref visitor = choose_parameter(get_parameter_reference(np, internal_np::visitor), default_visitor);
-  // constexpr bool has_visitor = !std::is_same_v<Default_visitor, std::remove_cv_t<std::remove_reference_t<Visitor_ref>>>;
+  constexpr bool has_visitor = !std::is_same_v<Default_visitor, std::remove_cv_t<std::remove_reference_t<Visitor_ref>>>;
 
   auto vpm = choose_parameter(get_parameter(np, internal_np::vertex_point),
                               get_property_map(vertex_point, pm));
 
   // config flags
-  // bool triangulate = !choose_parameter(get_parameter(np, internal_np::do_not_triangulate_faces), true);
+  bool triangulate = !choose_parameter(get_parameter(np, internal_np::do_not_triangulate_faces), true);
 
   auto oriented_side = traits.oriented_side_3_object();
   auto intersection_point = traits.construct_plane_line_intersection_point_3_object();
@@ -532,7 +737,7 @@ clip_convex(PolygonMesh& pm,
     return *vertices(pm).begin();
 
   } else if(faces(pm).size()==1){
-    // Dimension == 2
+    // Dimension == 2 with only one face
     halfedge_descriptor h_start = halfedge(*faces(pm).begin(), pm);
     // halfedge_descriptor h = next(h_start, pm);
     halfedge_descriptor h = h_start;
@@ -590,6 +795,14 @@ clip_convex(PolygonMesh& pm,
       set_next(src_split, split_edge, pm);
       set_next(opposite(split_edge, pm), opposite(src_split, pm), pm);
       set_halfedge(*faces(pm).begin(), split_edge, pm);
+      if(triangulate){
+        if constexpr(has_visitor){
+          Wrap_visitor_for_triangulation wrap_visitor(pm, visitor);
+          triangulate_face(face(split_edge, pm), pm, parameters::vertex_point_map(vpm).visitor(wrap_visitor));
+        } else {
+          triangulate_face(face(split_edge, pm), pm, parameters::vertex_point_map(vpm));
+        }
+      }
     }
 
     if(vertices(pm).size() < 3){
diff --git a/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/refine_with_plane.h b/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/refine_with_plane.h
index 62d2ca97bf43..1d0ca1c95d9a 100644
--- a/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/refine_with_plane.h
+++ b/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/refine_with_plane.h
@@ -211,7 +211,7 @@ struct Orthogonal_cut_plane_traits
  *    \cgalParamNEnd
  *
  *    \cgalParamNBegin{use_convex_specialization}
- *      \cgalParamDescription{If set to `true`, a faster implementation specialized for convex meshes is used. The input mesh must be convex to guarantee a correct execution and results.}
+ *      \cgalParamDescription{If set to `true`, a faster implementation specialized for convex meshes is used. The input mesh must be convex and have no flat vertices to guarantee a correct execution and results.}
  *      \cgalParamType{Boolean}
  *      \cgalParamDefault{`false`}
  *      \cgalParamExtra{convex specialization is only used if `edge_is_constrained_map`, `edge_is_marked_map` and `vertex_oriented_side_map` are unused.}
diff --git a/PMP_Boolean_operations/test/PMP_Boolean_operations/test_mesh_kernel.cpp b/PMP_Boolean_operations/test/PMP_Boolean_operations/test_mesh_kernel.cpp
index c0346c0af1b8..e95b94743d08 100644
--- a/PMP_Boolean_operations/test/PMP_Boolean_operations/test_mesh_kernel.cpp
+++ b/PMP_Boolean_operations/test/PMP_Boolean_operations/test_mesh_kernel.cpp
@@ -26,6 +26,197 @@ rotation(double a, double b, double c)
   return aff;
 }
 
+using Surface_mesh = SM;
+struct Test_visitor: public CGAL::Polygon_mesh_processing::Corefinement::Default_visitor<Surface_mesh>
+{
+  using halfedge_descriptor = typename boost::graph_traits<Surface_mesh>::halfedge_descriptor;
+  using face_descriptor = typename boost::graph_traits<Surface_mesh>::face_descriptor;
+  using vertex_descriptor = typename boost::graph_traits<Surface_mesh>::vertex_descriptor;
+
+  Surface_mesh& sm;
+  Surface_mesh::Property_map<Surface_mesh::Face_index, int> fid_map;
+  Surface_mesh::Property_map<Surface_mesh::Vertex_index, int> vid_map;
+  Surface_mesh::Property_map<Surface_mesh::Halfedge_index, int> hid_map;
+
+  int fid=-1;
+  int hid=-1;
+  std::size_t nbf=0;
+  std::size_t nbe=0;
+  std::size_t nb_input_v=sm.number_of_vertices();
+  std::size_t nb_input_v_on=0;
+  std::size_t nb_new_v_on=0;
+
+  Test_visitor(Surface_mesh& sm)
+    : sm(sm)
+  {
+    bool is_new=false;
+    std::tie(fid_map, is_new)=sm.add_property_map<Surface_mesh::Face_index,int>("f:id", -1);
+    assert(is_new);
+    int i=0; for (auto f : faces(sm)) put(fid_map, f, i++);
+    std::tie(vid_map, is_new)=sm.add_property_map<Surface_mesh::Vertex_index,int>("v:id", -1);
+    assert(is_new);
+    i=0; for (auto v : vertices(sm)) put(vid_map, v, i++);
+    std::tie(hid_map, is_new)=sm.add_property_map<Surface_mesh::Halfedge_index,int>("h:id", -1);
+    assert(is_new);
+    i=0;
+    for (auto e : edges(sm))
+    {
+      auto h = halfedge(e, sm);
+      put(hid_map, opposite(h, sm), i);
+      put(hid_map, h, i++);
+    }
+  }
+
+  void before_subface_creations(face_descriptor f_split, const Surface_mesh& pm)
+  {
+    fid=get(fid_map, f_split);
+    assert(fid!=-1);
+    assert(&pm==&sm);
+  }
+  void after_subface_creations(const Surface_mesh& pm)
+  {
+    assert(&pm==&sm);
+  }
+  void before_subface_created(const Surface_mesh& pm)
+  {
+    nbf=sm.number_of_faces();
+    assert(&pm==&sm);
+  }
+  void after_subface_created(face_descriptor f_new, const Surface_mesh& pm)
+  {
+    assert(&pm==&sm);
+    assert(get(fid_map, f_new)==-1);
+    put(fid_map, f_new, fid);
+    assert(nbf+1==sm.number_of_faces());
+  }
+
+  void before_edge_split(halfedge_descriptor h, Surface_mesh& pm)
+  {
+    hid=get(hid_map, h);
+    assert(hid!=-1);
+    assert(&pm==&sm);
+    nbe=sm.number_of_edges();
+  }
+  void edge_split(halfedge_descriptor hnew, Surface_mesh& pm)
+  {
+    assert(&pm==&sm);
+    assert(get(hid_map, hnew)==-1);
+    assert(hid!=-1);
+    put(hid_map, hnew, hid);
+    put(hid_map, opposite(hnew, sm), hid);
+    assert(nbe+1==sm.number_of_edges());
+  }
+  void after_edge_split()
+  {
+    assert(nbe+1==sm.number_of_edges());
+  }
+  void add_retriangulation_edge(halfedge_descriptor hnew, const Surface_mesh& pm)
+  {
+    assert(get(hid_map, hnew)==-1);
+    put(hid_map, hnew, -2);
+    put(hid_map, opposite(hnew, sm), -2);
+    assert(&pm==&sm);
+  }
+
+  void intersection_point_detected(std::size_t /* i_id */,
+                                   int /* sdim */,
+                                   halfedge_descriptor h_e,
+                                   halfedge_descriptor h_f,
+                                   const Surface_mesh& tm_e,
+                                   const Surface_mesh& tm_f,
+                                   bool is_target_coplanar,
+                                   bool is_source_coplanar)
+  {
+    assert(is_source_coplanar==false);
+    assert(h_f==boost::graph_traits<Surface_mesh>::null_halfedge());
+    assert(h_e!=boost::graph_traits<Surface_mesh>::null_halfedge());
+    if (!is_target_coplanar)
+      ++nb_new_v_on;
+    else
+      ++nb_input_v_on;
+    assert(&tm_e==&sm);
+    assert(&tm_f==&sm);
+  }
+
+  void new_vertex_added(std::size_t id, vertex_descriptor v, const Surface_mesh& pm)
+  {
+    assert(&pm==&sm);
+    assert(get(vid_map, v)==-1);
+    put(vid_map, v, nb_input_v+id-nb_input_v_on);
+  }
+
+  void before_face_copy(face_descriptor f, const Surface_mesh& src, const Surface_mesh& tgt)
+  {
+    assert(f==boost::graph_traits<Surface_mesh>::null_face());
+    assert(&src==&sm);
+    assert(&tgt==&sm);
+  }
+
+  void after_face_copy(face_descriptor fsrc, const Surface_mesh& src, face_descriptor ftgt, const Surface_mesh& tgt)
+  {
+    assert(fsrc==boost::graph_traits<Surface_mesh>::null_face());
+    assert(ftgt!=boost::graph_traits<Surface_mesh>::null_face());
+    //assert(get(fid_map, ftgt)==-1);
+    put(fid_map, ftgt, -2);
+    assert(&src==&sm);
+    assert(&tgt==&sm);
+    for (auto h : halfedges_around_face(halfedge(ftgt, sm), sm))
+    {
+      if (get(hid_map, h)==-1)
+        put(hid_map, h, -2);
+    }
+  }
+
+  void before_edge_duplicated(halfedge_descriptor h, Surface_mesh& tm)
+  {
+    hid = get(hid_map, h);
+    assert(hid!=-1);
+    assert(&tm==&sm);
+  }
+
+  void after_edge_duplicated(halfedge_descriptor h, halfedge_descriptor new_hedge, Surface_mesh& tm)
+  {
+    assert(hid==get(hid_map, h));
+    assert(&tm==&sm);
+    put(hid_map, new_hedge, hid);
+    put(hid_map, opposite(new_hedge, sm), hid);
+  }
+
+  void before_vertex_copy(vertex_descriptor v, Surface_mesh& src, Surface_mesh& tgt)
+  {
+    assert(&src==&sm);
+    assert(&tgt==&sm);
+    assert(get(vid_map, v)!=-1);
+  }
+  void after_vertex_copy(vertex_descriptor v, Surface_mesh& src, vertex_descriptor nv, Surface_mesh& tgt)
+  {
+    assert(&src==&sm);
+    assert(&tgt==&sm);
+    assert(get(vid_map, v)!=-1);
+    put(vid_map, nv, get(vid_map, v));
+  }
+
+
+
+  void check()
+  {
+    for (auto f :faces(sm))
+    {
+      if (get(fid_map, f)==-1) std::cout << sm.point(source(halfedge(f, sm), sm)) << " " << sm.point(target(halfedge(f, sm), sm)) << " " << sm.point(target(next(halfedge(f, sm), sm), sm)) << "\n";
+      assert(get(fid_map, f)!=-1);
+    }
+    for (auto h :halfedges(sm))
+    {
+      if (get(hid_map, h)==-1)
+        std::cout << sm.point(source(h, sm)) << " " << sm.point(target(h,sm)) << "\n";
+      assert(get(hid_map, h)!=-1);
+    }
+    std::size_t nbv_max=nb_input_v+nb_new_v_on;
+    for (auto v :vertices(sm))
+      assert(get(vid_map, v)!=-1 && get(vid_map, v)<(int)nbv_max);
+  }
+
+};
 
 std::size_t i=0;
 template<class Mesh, class Plane_3>
@@ -41,6 +232,16 @@ void test_clip_convex_on_mesh(const Mesh &m, const Plane_3 &pl, std::size_t expe
   assert(edges(m_copy).size()    == expected_nb_edges);
   assert(faces(m_copy).size()    == expected_nb_faces);
 
+  if constexpr(std::is_same_v<Mesh, SM>){
+    m_copy = m;
+    Test_visitor visitor(m_copy);
+    PMP::internal::clip_convex(m_copy, Pl(pl[0], pl[1], pl[2]), CGAL::parameters::visitor(visitor));
+
+    assert(vertices(m_copy).size() == expected_nb_vertices);
+    assert(edges(m_copy).size()    == expected_nb_edges);
+    assert(faces(m_copy).size()    == expected_nb_faces);
+  }
+
   m_copy = m;
   PMP::internal::clip_convex(m_copy, Pl(pl[0], pl[1], pl[2]));
 
@@ -49,6 +250,29 @@ void test_clip_convex_on_mesh(const Mesh &m, const Plane_3 &pl, std::size_t expe
   assert(faces(m_copy).size()    == expected_nb_faces);
 }
 
+template<class Mesh, class Plane_3>
+void test_clip_convex_with_triangulation_on_mesh(const Mesh &m, const Plane_3 &pl){
+  using K =typename CGAL::Kernel_traits<typename Plane_3::value_type>::Kernel;
+  using Pl = typename K::Plane_3;
+
+  auto m_copy = m;
+  PMP::triangulate_faces(m_copy);
+  PMP::internal::clip_convex(m_copy, Pl(pl[0], pl[1], pl[2]));
+  PMP::triangulate_faces(m_copy);
+
+  std::size_t expected_nb_vertices = vertices(m_copy).size();
+  std::size_t expected_nb_edges    = edges(m_copy).size();
+  std::size_t expected_nb_faces    = faces(m_copy).size();
+
+  m_copy = m;
+  PMP::triangulate_faces(m_copy);
+  PMP::internal::clip_convex(m_copy, Pl(pl[0], pl[1], pl[2]), CGAL::parameters::do_not_triangulate_faces(false));
+
+  assert(vertices(m_copy).size() == expected_nb_vertices);
+  assert(edges(m_copy).size()    == expected_nb_edges);
+  assert(faces(m_copy).size()    == expected_nb_faces);
+}
+
 template<class Mesh>
 void test_kernel_on_mesh(const Mesh &input, std::size_t expected_nb_vertices, std::size_t expected_nb_edges, std::size_t expected_nb_faces, double expected_volume = 0){
   assert(PMP::has_empty_kernel(input, CGAL::parameters::allow_open_input(true)) == (expected_nb_vertices == 0));