rewrite/refactor a lot of the code

STL_Extension/include/CGAL/Has_timestamp.h

@@ -36,6 +36,9 @@ namespace internal {
     // when T does not have a Has_timestamp tag
   {};
 
+  template <typename T>
+  constexpr bool has_timestamp_v = Has_timestamp<T>::value;
+
 } // end namespace internal
 } // end namespace CGAL
 

Stream_support/include/CGAL/IO/io.h

@@ -216,6 +216,19 @@ public:
   }
 };
 
+template <class Func>
+class Output_rep<Func, IO_manip_tag>
+{
+  Func f;
+
+public:
+  Output_rep(Func f) : f(f) {}
+  std::ostream& operator()(std::ostream& os) const
+  {
+    return f(os);
+  }
+};
+
 /*!
   \relates Output_rep
   \brief stream output of the \c Output_rep calls its \c operator().

Stream_support/include/CGAL/IO/io_tags.h

@@ -51,7 +51,7 @@ template<> struct Io_traits<float> { typedef io_Read_write Io_tag; };
 template<> struct Io_traits<double> { typedef io_Read_write Io_tag; };
 template<> struct Io_traits<long double> { typedef io_Read_write Io_tag; };
 
-
+struct IO_manip_tag{};
 
 } //namespace CGAL
 

Triangulation_2/include/CGAL/Constrained_triangulation_plus_2.h

@@ -378,57 +378,69 @@ public:
   }
 
 
+
+  // Removes pos from the constraint cid.
+  // Returns the iterator to vertex that was just after pos (or end())
+  // Writes the modified faces to out
   template <class OutputIterator>
   Vertices_in_constraint_iterator
   remove_vertex_from_constraint(Constraint_id cid, Vertices_in_constraint_iterator pos,
                                 OutputIterator out)
   {
     if(pos == vertices_in_constraint_begin(cid)){
-      ++pos;
+      // cid is [P, A, ..., B] -> split to aux=[P, A] and cid=[A...B]
-      Constraint_id aux = hierarchy.split2(cid,pos);
+      Constraint_id aux = hierarchy.split2(cid, std::next(pos));
       remove_constraint(aux, out);
-      return pos;
+      return vertices_in_constraint_begin(cid);
     }
 
-    Vertices_in_constraint_iterator it = vertices_in_constraint_end(cid);
+    if(pos == std::prev(vertices_in_constraint_end(cid))){
-    it--;
+      // cid is [A, ..., B, P] -> split to cid=[A...B] and aux=[B,P]
-    if(pos == it){
+      Constraint_id aux = hierarchy.split(cid, std::prev(pos));
-      --pos;
-      Constraint_id aux = hierarchy.split(cid, pos);
       remove_constraint(aux, out);
       return vertices_in_constraint_end(cid);
     }
 
-    Vertices_in_constraint_iterator pp = pos;
+    const auto second_vertex_of_cid = std::next(vertices_in_constraint_begin(cid));
-    --pp;
+    const auto next_to_last_vertex_of_cid = std::prev(vertices_in_constraint_end(cid), 2);
-    Vertex_handle a = *pp;
-    pp = pos;
-    ++pp;
-    Vertex_handle b = *pp;
-    --it;
-    Vertices_in_constraint_iterator beg = vertices_in_constraint_begin(cid);
-    ++beg;
-    Face_container<Constrained_triangulation_plus_2> fc(*this);
 
     Constraint_id head = nullptr, tail = nullptr;
-    if(pos != beg){
+    if(pos != second_vertex_of_cid){
+      // cid is [A, ..., B, P, C, ..., D]
+      // split to:
+      //    head = [A...B] and,
+      //    cid =  [B, P, C...D]
       // split off head
-      --pos;
+      head = hierarchy.split2(cid, std::prev(pos));
-      head = hierarchy.split2(cid, pos);
-      ++pos;
     }
-    if(pos != it){
+    if(pos != next_to_last_vertex_of_cid){
+      // cid is now [B, P, C, ..., D]
+      // split to:
+      //    cid = [B, P, C] and,
+      //    tail = [C...D]
       // split off tail
-      ++pos;
+      tail = hierarchy.split(cid,std::next(pos));
-      tail = hierarchy.split(cid,pos);
     }
 
-    Constraint_id aux = insert_constraint(a, b, std::back_inserter(fc));
+    // now:
-    pos = vertices_in_constraint_end(aux);
+    //   cid is [B, P, C]
-    --pos;
+    //   head is null or [A...B]
-    --pos; // this is not necessarily == vertices_in_constraint_begin(aux);
+    //   tail is null or [C...D]
+    // Let create insert [C,D] and conditionnaly concatenate head and tail,
+    // and return the iterator to C
+
+    Vertex_handle b = *std::prev(pos);
+    Vertex_handle c = *std::next(pos);
+
+    Face_container<Constrained_triangulation_plus_2> fc(*this);
+    Constraint_id aux = insert_constraint(b, c, std::back_inserter(fc));
+
+    auto pos_before_c = std::prev(vertices_in_constraint_end(aux), 2);
+    // `pos_before_c` is not necessarily == vertices_in_constraint_begin(aux)
+    // there might have been intersecting constraints
+
     hierarchy.swap(cid, aux);
-    remove_constraint(aux, std::back_inserter(fc));
+    remove_constraint(aux, std::back_inserter(fc)); // removes [B, P, C]
 
     if(head != nullptr){
       hierarchy.concatenate2(head, cid);
@@ -438,8 +450,9 @@ public:
       hierarchy.concatenate(cid, tail);
     }
     fc.write_faces(out);
-    ++pos; // we went one too far back because the last vertex gets removed by concatenate
+
-    return pos;
+    // we went one too far back because the last vertex `c` gets removed by concatenate
+    return std::next(pos_before_c);
   }
 
   // Inserts vh before pos
@@ -461,61 +474,56 @@ public:
     // Insertion after the last vertex
     if(pos == vertices_in_constraint_end(cid)){
       //std::cout << "insertion after last vertex" << std::endl;
-      pos--;
+      Constraint_id tail = insert_constraint(*std::prev(pos), vh, out);
-      Constraint_id tail = insert_constraint(*pos, vh, out);
+      auto returned_pos = std::prev(vertices_in_constraint_end(tail));
-      pos = vertices_in_constraint_end(tail);
-      --pos;
       hierarchy.concatenate(cid, tail);
-      return pos;
+      return returned_pos;
     }
+    Vertices_in_constraint_iterator pred = std::prev(pos);
+    Vertices_in_constraint_iterator last = std::prev(vertices_in_constraint_end(cid));
+    Vertex_handle a = *pred;
     Vertex_handle b = *pos;
-    --pos;
+
-    Vertex_handle a = *pos;
+    // cid is [..., A, B, ...] and M=*vh will be inserted between A and B
-    ++pos;
+
     Face_container<Constrained_triangulation_plus_2> fc(*this);
-    Vertices_in_constraint_iterator beg = vertices_in_constraint_begin(cid), vcit;
+    Constraint_id aux1 = insert_constraint(a, vh, std::back_inserter(fc));
-    ++beg;
+    Constraint_id aux2 = insert_constraint(vh, b, std::back_inserter(fc));
-    vcit = beg;
+    auto returned_pos = vertices_in_constraint_begin(aux2);
-    ++beg;
+    concatenate(aux1, aux2);
+    // here:
+    //   aux1 is [A, M, B]
+    //   aux2 is empty
+    // and returned_pos is the iterator to M
+
+    const auto second_vertex_of_cid = std::next(vertices_in_constraint_begin(cid));
     // If the constraint consists only of a segment, and we want to insert
-    // in the middle
+    // in the middle: cid is just the segment [A, B]
-    if((pos == vcit) && (beg == vertices_in_constraint_end(cid))){
+    if((pos == second_vertex_of_cid) && (second_vertex_of_cid == last)){
       //std::cout << "insertion in constraint which is a segment" << std::endl;
-      Constraint_id aux1 = insert_constraint(a, vh, std::back_inserter(fc));
-      Constraint_id aux2 = insert_constraint(vh, b, std::back_inserter(fc));
-      pos = vertices_in_constraint_begin(aux2);
-      concatenate(aux1, aux2);
       hierarchy.swap(cid, aux1);
       remove_constraint(aux1, std::back_inserter(fc));
       fc.write_faces(out);
-      return pos;
+      return returned_pos;
-
     }
+
     Constraint_id head = nullptr, tail = nullptr;
-    Vertices_in_constraint_iterator bit = vertices_in_constraint_begin(cid);
+    if(pos != second_vertex_of_cid){
-    Vertices_in_constraint_iterator pred = pos;
-    --pred;
-    ++bit;
-    if(pos != bit){
       //std::cout << "split head" << std::endl;
       head = split(cid, pred);
       std::swap(head,cid); // split2 does the job
       pred = vertices_in_constraint_begin(cid);
-      pos = pred;
+      pos = std::next(pred);
-      ++pos;
     }
-    Vertices_in_constraint_iterator eit = vertices_in_constraint_end(cid);
+    // head is now [..., A] or null
-    --eit;
+    //  cid is now [A, B, ...]
-    if(pos != eit){
+    if(pos != last){
       //std::cout << "split tail" << std::endl;
       tail = split(cid, pos);
     }
-
+    // head is now [..., A] or null
-    // make the new constraint
+    //  cid is now [A, B]
-    Constraint_id aux1 = insert_constraint(a, vh, std::back_inserter(fc));
+    // tail is now [B, ...] or null
-    Constraint_id aux2 = insert_constraint(vh, b, std::back_inserter(fc));
-    pos = vertices_in_constraint_begin(aux2);
-    concatenate(aux1, aux2);
 
     if(head != nullptr){
       //std::cout << "concatenate head" << std::endl;
@@ -988,7 +996,7 @@ insert_subconstraint(Vertex_handle vaa,
 
   //to debug
 public:
-  void print_hierarchy() { hierarchy.print(); }
+  void print_hierarchy(std::ostream& os = std::cout) { hierarchy.print(os); }
 
   //template member functions
 public:

Triangulation_2/include/CGAL/Triangulation_2/internal/Polyline_constraint_hierarchy_2.h

@@ -28,6 +28,8 @@
 #include <CGAL/Skiplist.h>
 #include <CGAL/Iterator_range.h>
 #include <CGAL/assertions.h>
+#include <CGAL/Has_timestamp.h>
+#include <CGAL/IO/io.h>
 
 #ifdef CGAL_CDT_2_DEBUG_INTERSECTIONS
 #  include <CGAL/IO/io.h>
@@ -62,7 +64,6 @@ private:
   private:
     Vertex_handle vertex_;
     Point point_;
-  public:
     bool input_;
   };
 
@@ -213,6 +214,7 @@ public:
   typedef typename Constraints_set::iterator Constraint_iterator;
   typedef const Constraints_set& Constraints;
   typedef typename Sc_to_c_map::const_iterator    Sc_iterator;
+  typedef typename Sc_to_c_map::iterator Sc_it;
   typedef Sc_iterator Subconstraint_and_contexts_iterator;
   typedef const Sc_to_c_map& Subconstraints_and_contexts;
 
@@ -447,7 +449,53 @@ public:
   void swap(Polyline_constraint_hierarchy_2& ch);
 
 private:
+  template <typename F>
+  void for_context_lists_of_all_subconstraints(Constraint_id cid, const F& f)
+  {
+    auto vl = cid.vl_ptr();
+    for(Vertex_it it = vl->skip_begin(), succ = it, end = vl->skip_end(); ++succ != end; ++it) {
+      auto scit = sc_to_c_map.find(sorted_pair(*it, *succ));
+      CGAL_assertion(scit != sc_to_c_map.end());
+      Context_list* hcl = scit->second;
+      f(hcl, it, scit);
+    }
+  }
+
+  static void replace_first_in_context_list(Context_list* hcl, Constraint_id old_id, Constraint_id new_id)
+  {
+    // std::find_if is a sort of std::for_each with a break
+    [[maybe_unused]] auto it = std::find_if(hcl->begin(), hcl->end(), [&](Context& ctxt) {
+      if(ctxt.enclosing == old_id) {
+        ctxt.enclosing = new_id;
+        return true;
+      }
+      return false;
+    });
+  }
+
+  static void update_first_context_position(Context_list* hcl, Constraint_id id, Vertex_it new_pos)
+  {
+    [[maybe_unused]] auto it = std::find_if(hcl->begin(), hcl->end(), [&](Context& ctxt) {
+      if(ctxt.enclosing == id) {
+        ctxt.pos = new_pos;
+        return true;
+      }
+      return false;
+    });
+  }
+
+  static void remove_first_in_context_list(Context_list* hcl, Constraint_id id)
+  {
+    auto it = std::find_if(hcl->begin(), hcl->end(), [&](Context& ctxt) {
+      return ctxt.enclosing == id;
+    });
+    if(it != hcl->end()) {
+      hcl->erase(it);
+    }
+  }
+
   Constraint_id new_constraint_id() const {
+    // TODO: handle ids
     auto id = number_of_constraints() == 0 ? 0 : constraints_set.rbegin()->id + 1;
     return Constraint_id(new Vertex_list, id);
   }
@@ -464,7 +512,7 @@ private:
 
   //to_debug
 public:
-  void   print() const;
+  void   print(std::ostream& os = std::cout) const;
 };
 
 template <class T, class Compare, class Point>
@@ -645,29 +693,18 @@ template <class T, class Compare, class Point>
 void
 Polyline_constraint_hierarchy_2<T,Compare,Point>::
 swap(Constraint_id constr_a, Constraint_id constr_b) {
-  auto substitute_enclosing_in_vertex_list = [this](Vertex_list_ptr vl, Constraint_id old_id, Constraint_id new_id) {
+  auto substitute_enclosing_in_vertex_list = [this](Constraint_id cid, Constraint_id old_id, Constraint_id new_id) {
-    // We have to look at all subconstraints
+    for_context_lists_of_all_subconstraints(cid, [&](Context_list* hcl, Vertex_it, Sc_it) {
-    for(Vertex_it it = vl->skip_begin(), succ = it, end = vl->skip_end(); ++succ != end; ++it) {
+      replace_first_in_context_list(hcl, old_id, new_id);
-      typename Sc_to_c_map::iterator scit = this->sc_to_c_map.find(sorted_pair(*it, *succ));
+    });
-      CGAL_assertion(scit != this->sc_to_c_map.end());
-      Context_list* hcl = scit->second;
-
-      // and replace the context of the constraint
-      for(Context_iterator ctit = hcl->begin(); ctit != hcl->end(); ctit++) {
-        if(ctit->enclosing == old_id) {
-          ctit->enclosing = new_id;
-          break;
-        }
-      }
-    }
   };
 
   Vertex_list_ptr constr_a_vl = constr_a.vl_ptr();
   Vertex_list_ptr constr_b_vl = constr_b.vl_ptr();
 
-  substitute_enclosing_in_vertex_list(constr_a_vl, constr_a, nullptr);
+  substitute_enclosing_in_vertex_list(constr_a, constr_a, nullptr);
-  substitute_enclosing_in_vertex_list(constr_b_vl, constr_b, constr_a);
+  substitute_enclosing_in_vertex_list(constr_b, constr_b, constr_a);
-  substitute_enclosing_in_vertex_list(constr_a_vl, nullptr, constr_b);
+  substitute_enclosing_in_vertex_list(constr_a, nullptr, constr_b);
 
   constr_a_vl->swap(*constr_b_vl);
 }
@@ -675,33 +712,22 @@ swap(Constraint_id constr_a, Constraint_id constr_b) {
 template <class T, class Compare, class Point>
 void
 Polyline_constraint_hierarchy_2<T,Compare,Point>::
-remove_constraint(Constraint_id cid){
+remove_constraint(Constraint_id cid)
+{
   constraints_set.erase(cid);
 
-  // We have to look at all subconstraints
+  for_context_lists_of_all_subconstraints(cid, [&](Context_list* hcl, Vertex_it, Sc_it scit) {
-  for(Vertex_it it = cid.begin(), succ = it, end = cid.end();
+    remove_first_in_context_list(hcl, cid);
-      ++succ != end;
-      ++it){
-    typename Sc_to_c_map::iterator scit = sc_to_c_map.find(sorted_pair(*it,*succ));
-    CGAL_assertion(scit != sc_to_c_map.end());
-    Context_list* hcl = scit->second;
 
-    // and remove the context of the constraint
-    for(Context_iterator ctit=hcl->begin(); ctit != hcl->end(); ctit++) {
-      if(ctit->enclosing == cid){
-            hcl->erase(ctit);
-                break;
-      }
-    }
     // If the constraint passes several times through the same subconstraint,
-    // the above loop maybe removes them in the wrong order
+    // the above call to `remove_first_in_context_list` may remove them in the wrong order.
 
     // If this was the only context in the list, delete the context list
-    if(hcl->empty()){
+    if(hcl->empty()) {
       sc_to_c_map.erase(scit);
       delete hcl;
     }
-  }
+  });
   delete cid.vl_ptr();
 }
 
@@ -802,53 +828,31 @@ template <class T, class Compare, class Point>
 typename Polyline_constraint_hierarchy_2<T,Compare,Point>::Constraint_id
 Polyline_constraint_hierarchy_2<T,Compare,Point>::concatenate(Constraint_id constr_a, Constraint_id constr_b)
 {
+  // constr_a is [A, ..., M]
+  // constr_b is [M, ..., B]
+  // we want:
+  //   constr_a = [A, ..., M, ..., B]
+  //   constr_b = []
+
   // std::cerr << std::format("concatenate({}, {}) ", constr_a.id, constr_b.id) << std::endl;
   Vertex_list_ptr constr_a_vl = constr_a.vl_ptr();
   Vertex_list_ptr constr_b_vl = constr_b.vl_ptr();
-  constraints_set.erase(constr_a);
   constraints_set.erase(constr_b);
-  // We have to look at all subconstraints
-  for(Vertex_it it = constr_b_vl->skip_begin(), succ = it, end = constr_b_vl->skip_end();
-      ++succ != end;
-      ++it){
-    typename Sc_to_c_map::iterator scit = sc_to_c_map.find(sorted_pair(*it,*succ));
-    CGAL_assertion(scit != sc_to_c_map.end());
-    Context_list* hcl = scit->second;
 
-    // and replace the context of the constraint
+  for_context_lists_of_all_subconstraints(constr_b, [&](Context_list* hcl, Vertex_it, Sc_it) {
-    for(Context_iterator ctit=hcl->begin(); ctit != hcl->end(); ctit++) {
+    replace_first_in_context_list(hcl, constr_b, constr_a);
-      if(ctit->enclosing == constr_b){
+  });
-        ctit->enclosing = constr_a;
+
-        break;
-      }
-    }
-  }
   // now we really concatenate the vertex lists
   // Note that all iterators pointing into constr_a remain valid.
   // This concerns user code, as well as  the data member "pos" of the Context class
+  CGAL_assertion(constr_a_vl->back().vertex() == constr_b_vl->front().vertex());
   constr_a_vl->pop_back(); // because it is the same as constr_b_vl.front()
-  Vertex_it back_it = constr_a_vl->skip_end();
+  constr_a_vl->splice(constr_a_vl->skip_end(), *constr_b_vl, constr_b_vl->skip_begin(), constr_b_vl->skip_end());
-  --back_it;
-  constr_a_vl->splice(constr_a_vl->skip_end(), *(constr_b_vl), constr_b_vl->skip_begin(), constr_b_vl->skip_end());
 
-  // Note that for VC8 with iterator debugging the iterators pointing into constr_b
+  for_context_lists_of_all_subconstraints(constr_a, [&](Context_list* hcl, Vertex_it it, Sc_it) {
-  // are NOT valid      So we have to update them
+    update_first_context_position(hcl, constr_a, it);
-  for(Vertex_it it = back_it, succ = it, end = constr_a_vl->skip_end();
+  });
-      ++succ != end;
-      ++it){
-    typename Sc_to_c_map::iterator scit = sc_to_c_map.find(sorted_pair(*it,*succ));
-    CGAL_assertion(scit != sc_to_c_map.end());
-    Context_list* hcl = scit->second;
-
-    // and update pos in the context of the constraint
-    for(Context_iterator ctit=hcl->begin(); ctit != hcl->end(); ctit++) {
-      if(ctit->enclosing == constr_a){
-        ctit->pos = it;
-        break;
-      }
-    }
-  }
-  constraints_set.insert(constr_a);
 
   delete constr_b_vl;
   return constr_a;
@@ -858,54 +862,34 @@ template <class T, class Compare, class Point>
 typename Polyline_constraint_hierarchy_2<T,Compare,Point>::Constraint_id
 Polyline_constraint_hierarchy_2<T,Compare,Point>::concatenate2(Constraint_id constr_a, Constraint_id constr_b)
 {
+  // constr_a is [A, ..., M]
+  // constr_b is [M, ..., B]
+  // we want:
+  //   constr_a =
+  //   constr_b = [A, ..., M, ..., B]
+
   Vertex_list_ptr constr_a_vl = constr_a.vl_ptr();
   Vertex_list_ptr constr_b_vl = constr_b.vl_ptr();
-  constraints_set.erase(constr_a);
+  constraints_set.erase(constr_a); // DIFF
-  constraints_set.erase(constr_b);
+
-  // We have to look at all subconstraints
+  for_context_lists_of_all_subconstraints(constr_a, [&](Context_list* hcl, Vertex_it, Sc_it) { // DIFF
-  for(Vertex_it it = constr_a_vl->skip_begin(), succ = it, end = constr_a_vl->skip_end();
+    replace_first_in_context_list(hcl, constr_a, constr_b); // DIFF
-      ++succ != end;
+  });
-      ++it){
-    typename Sc_to_c_map::iterator scit = sc_to_c_map.find(sorted_pair(*it,*succ));
-    CGAL_assertion(scit != sc_to_c_map.end());
-    Context_list* hcl = scit->second;
 
-    // and replace the context of the constraint
-    for(Context_iterator ctit=hcl->begin(); ctit != hcl->end(); ctit++) {
-      if(ctit->enclosing == constr_a){
-        ctit->enclosing = constr_b;
-        break;
-      }
-    }
-  }
   // now we really concatenate the vertex lists
   // Note that all iterators pointing into constr_b remain valid.
+  CGAL_assertion(constr_a_vl->back().vertex() == constr_b_vl->front().vertex());
   constr_a_vl->pop_back(); // because it is the same as constr_b_vl.front()
-  Vertex_it back_it = constr_b_vl->skip_begin();
+  constr_b_vl->splice(constr_b_vl->skip_begin(), *constr_a_vl, constr_a_vl->skip_begin(), constr_a_vl->skip_end()); // DIFF
-
-  constr_b_vl->splice(constr_b_vl->skip_begin(), *(constr_a_vl), constr_a_vl->skip_begin(), constr_a_vl->skip_end());
 
   // Note that for VC8 with iterator debugging the iterators pointing into constr_a
   // are NOT valid      So we have to update them
-  for(Vertex_it it = constr_b_vl->skip_begin(), succ = it, end = back_it;
+  for_context_lists_of_all_subconstraints(constr_b /*DIFF*/, [&](Context_list* hcl, Vertex_it it, Sc_it) {
-      ++succ != end;
+    update_first_context_position(hcl, constr_b, it); // DIFF
-      ++it){
+  });
-    typename Sc_to_c_map::iterator scit = sc_to_c_map.find(sorted_pair(*it,*succ));
-    CGAL_assertion(scit != sc_to_c_map.end());
-    Context_list* hcl = scit->second;
 
-    // and update pos in the context of the constraint
+  delete constr_a_vl; // DIFF
-    for(Context_iterator ctit=hcl->begin(); ctit != hcl->end(); ctit++) {
+  return constr_b; // DIFF
-      if(ctit->enclosing == constr_b){
-        ctit->pos = it;
-        break;
-      }
-    }
-  }
-  constraints_set.insert(constr_b);
-
-  delete constr_a_vl;
-  return constr_b;
 }
 
 
@@ -916,35 +900,35 @@ template <class T, class Compare, class Point>
 typename Polyline_constraint_hierarchy_2<T,Compare,Point>::Constraint_id
 Polyline_constraint_hierarchy_2<T,Compare,Point>::split(Constraint_id constr, Vertex_it vcit)
 {
-  Constraint_id new_constr = new_constraint_id();
+  // constr is [A, ..., B], vcit points to M in [A, ..., B]
-  constraints_set.erase(constr);
+
-  Vertex_list_ptr new_vl = new_constr.vl_ptr();
+  Constraint_id new_constr = new_constraint_id();
-  Vertex_list_ptr constr_vl = constr.vl_ptr();
+  constraints_set.insert(new_constr);
-  new_vl->splice(new_vl->skip_end(), *(constr_vl), vcit.base(), constr_vl->skip_end());
+
-  constr_vl->push_back(new_vl->front()); // Duplicate the common vertex
+  Vertex_list_ptr constr_vl = constr.vl_ptr();
-  Vertex_it vit = new_vl->skip_begin();
+  Vertex_list_ptr new_vl = new_constr.vl_ptr();
-  vit.input() = true;
+
-  vit = constr_vl->skip_end();
+  vcit.input() = true;
-  --vit;
+  auto middle_node = Node(*vcit, true);
-  vit.input() = true;
+
-  constraints_set.insert(constr);
+  // Let's split, that way:
-  constraints_set.insert(new_constr);
+  //       constr = [A, ..., M]
- // We have to look at all subconstraints
+  //   new_constr = [M, ..., B]
-  for(Vertex_it it = new_vl->skip_begin(), succ = it, end = new_vl->skip_end();
+
-      ++succ != end;
+  // The splice does:
-      ++it){
+  //       constr = [A, ..., M[ (without M)
-    typename Sc_to_c_map::iterator scit = sc_to_c_map.find(sorted_pair(*it,*succ));
+  //   new_constr = [M, ..., B]
-    CGAL_assertion(scit != sc_to_c_map.end());
+  new_vl->splice(new_vl->skip_end(), *constr_vl, vcit.base(), constr_vl->skip_end());
-    Context_list* hcl = scit->second;
+  constr_vl->push_back(middle_node); // add M to the end of constr
+
+  CGAL_assertion(vcit.base() == new_vl->skip_begin());
+  CGAL_assertion(new_vl->front().input() == true);
+  CGAL_assertion(constr_vl->back().input() == true);
+
+  for_context_lists_of_all_subconstraints(new_constr, [&](Context_list* hcl, Vertex_it, Sc_it) {
+    replace_first_in_context_list(hcl, constr, new_constr);
+  });
 
-    // and replace the context of the constraint
-    for(Context_iterator ctit=hcl->begin(); ctit != hcl->end(); ctit++) {
-      if(ctit->enclosing == constr){
-        ctit->enclosing = new_constr;
-        break;
-      }
-    }
-  }
   return new_constr;
 }
 
@@ -952,35 +936,35 @@ template <class T, class Compare, class Point>
 typename Polyline_constraint_hierarchy_2<T,Compare,Point>::Constraint_id
 Polyline_constraint_hierarchy_2<T,Compare,Point>::split2(Constraint_id constr, Vertex_it vcit)
 {
-  Constraint_id new_constr = new_constraint_id();
+  // constr is [A, ..., B], vcit points to M in [A, ..., B]
-  constraints_set.erase(constr);
+
-  Vertex_list_ptr new_vl = new_constr.vl_ptr();
+  Constraint_id new_constr = new_constraint_id();
-  Vertex_list_ptr constr_vl = constr.vl_ptr();
+  constraints_set.insert(new_constr);
-  new_vl->splice(new_vl->skip_end(), *constr_vl, constr_vl->skip_begin(), vcit.base());
+
-  new_vl->push_back(constr_vl->front()); // Duplicate the common vertex
+  Vertex_list_ptr constr_vl = constr.vl_ptr();
-  Vertex_it vit = new_vl->skip_end();
+  Vertex_list_ptr new_vl = new_constr.vl_ptr();
-  --vit;
+
-  vit.input() = true;
+  vcit.input() = true;
-  vit = constr_vl->skip_begin();
+  auto middle_node = Node(*vcit, true);
-  vit.input() = true;
+
-  constraints_set.insert(constr);
+  // Let's split, that way:
-  constraints_set.insert(new_constr);
+  //       constr = [M, ..., B]
- // We have to look at all subconstraints
+  //   new_constr = [A, ..., M]
-  for(Vertex_it it = new_vl->skip_begin(), succ = it, end = new_vl->skip_end();
+
-      ++succ != end;
+  // The splice does:
-      ++it){
+  //       constr = [M, ..., B]
-    typename Sc_to_c_map::iterator scit = sc_to_c_map.find(sorted_pair(*it,*succ));
+  //   new_constr = [A, ..., M[  (without M)
-    CGAL_assertion(scit != sc_to_c_map.end());
+  new_vl->splice(new_vl->skip_end(), *constr_vl, constr_vl->skip_begin(), vcit.base());
-    Context_list* hcl = scit->second;
+  new_vl->push_back(middle_node); // add M to the end of new_constr
+
+  CGAL_assertion(vcit.base() == constr_vl->skip_begin());
+  CGAL_assertion(constr_vl->front().input() == true);
+  CGAL_assertion(new_vl->back().input() == true);
+
+  for_context_lists_of_all_subconstraints(new_constr, [&](Context_list* hcl, Vertex_it, Sc_it) {
+    replace_first_in_context_list(hcl, constr, new_constr);
+  });
 
-    // and replace the context of the constraint
-    for(Context_iterator ctit=hcl->begin(); ctit != hcl->end(); ctit++) {
-      if(ctit->enclosing == constr){
-        ctit->enclosing = new_constr;
-        break;
-      }
-    }
-  }
   return new_constr;
 }
 
@@ -1244,9 +1228,9 @@ oriented_end(T va, T vb, T& vc) const
 template <class T, class Compare, class Point>
 void
 Polyline_constraint_hierarchy_2<T,Compare,Point>::
-print() const
+print(std::ostream& os) const
 {
-  std::map<T,int>  vertex_num_map;
+  std::map<Vertex_handle,int>  vertex_num_map;
   int num = 0;
   for(const auto& cid : constraints()) {
     for (const auto& node : cid.elements()){
@@ -1254,51 +1238,51 @@ print() const
     }
   }
 
-  struct Vertex_num {
+  auto disp_vertex = [&vertex_num_map](Vertex_handle v) {
-    std::map<T,int>& vertex_num_map;
+    return CGAL::IO::oformat(
-    int operator[](T v) const {
+        [v, &vertex_num_map](auto& os) -> decltype(os)& {
-#ifndef CGAL_CDT_2_DEBUG_INTERSECTIONS
+          constexpr bool star_v_has_timestamp =
-      auto it = vertex_num_map.find(v);
+              internal::has_timestamp_v<CGAL::cpp20::remove_cvref_t<decltype(*v)>>;
-      if(it == vertex_num_map.end()) return -1;
+          if constexpr(star_v_has_timestamp) {
-      return it->second;
+            return os << '#' << v->time_stamp();
-#else
+          } else {
-      return v->time_stamp();
+            auto it = vertex_num_map.find(v);
-#endif
+            auto n = (it == vertex_num_map.end()) ? -1 : it->second;
-    }
+            return os << n;
-  } vertex_num{vertex_num_map};
+          }
-//  typename std::map<T,int>::iterator vnit = vertex_num.begin();
+        },
-//  for(; vnit != vertex_num.end(); vnit++) {
+        IO_manip_tag{});
-//    vnit->second = ++num;
+  };
-//    std::cerr << "vertex num " << num  << " " << vnit->first->point()
-//              << std::endl;
-//  }
 
+  os << "# number of constraints: " << number_of_constraints() << std::endl;
   for(const auto& cid : constraints()) {
-    std::cout << std::endl;
+    os << "constraint(" << cid.id << ") ";
-    std::cout << "constraint(" << cid.id << ") ";
+    os << cid.vl_ptr();
-    std::cout << cid.vl_ptr();
+    os << "\n  vertex list ";
-    std::cout << "  subconstraints ";
     for(const auto& node : cid.elements()) {
-      std::cout << vertex_num[node.vertex()] << " ";
+      os << disp_vertex(node.vertex()) << " ";
     }
-    std::cout << ":  ";
+    os << "\n  corresponding points:  ";
     for(const auto& node : cid.elements()) {
-      std::cout << node.point() << "   ";
+      os << node.point() << "   ";
     }
+    os << std::endl;
   }
-  std::cout << std::endl;
+  os << std::endl;
+  os << "# number of subconstraints: " << sc_to_c_map.size() << std::endl;
   for(const auto& subconstraint : subconstraints()) {
-    std::cout << "subconstraint ";
+    os << "subconstraint (";
-    std::cout << vertex_num[subconstraint.first] << " " << vertex_num[subconstraint.second];
+    os << disp_vertex(subconstraint.first) << ", "
+              << disp_vertex(subconstraint.second) << ")";
     Context_iterator cb, ce;
     get_contexts(subconstraint.first, subconstraint.second, cb, ce);
 
-    std::cout << "  enclosing ";
+    os << "  enclosing:  ";
     for(; cb != ce; cb++) {
-      std::cout << "(" << cb->id().id << ") " << cb->id().vl_ptr();
+      os << "(cid " << cb->id().id << ") " << cb->id().vl_ptr();
-      std::cout << "   ";
+      os << ", pos: " << std::distance(cb->vertices_begin(), cb->pos) << "   ";
     }
-    std::cout << std::endl;
+    os << std::endl;
   }
   return;
 }

Triangulation_2/test/Triangulation_2/issue_4025.cpp

@@ -1,89 +1,110 @@
-#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+// This code was first submitted in an issue:
-#include <CGAL/Triangulation_2.h>
+//    https://github.com/CGAL/cgal/issues/4025
+// and then rewrote a lot, keeping the observed behavior.
+
+#include <CGAL/Base_with_time_stamp.h>
 #include <CGAL/Constrained_Delaunay_triangulation_2.h>
 #include <CGAL/Constrained_triangulation_plus_2.h>
 
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 
-typedef CGAL::Exact_predicates_inexact_constructions_kernel                 K;
+#include <string_view>
-typedef CGAL::Polygon_2<K>                                                Polygon_2;
+#include <iostream>
-typedef CGAL::Exact_intersections_tag                                     Itag_;
+#include <iterator>
-typedef CGAL::Constrained_Delaunay_triangulation_2<K,CGAL::Default, Itag_> CDT;
-typedef CGAL::Constrained_triangulation_plus_2<CDT>                       CDTP;
 
-typedef CDTP::Point                                                       Point;
+using K = CGAL::Exact_predicates_inexact_constructions_kernel;
-typedef CDTP::Constraint_id                                               Cid;
+using Polygon_2 = CGAL::Polygon_2<K>;
-typedef CDTP::Vertex_handle                                               Vertex_handle;
+using Itag_ = CGAL::Exact_intersections_tag;
-typedef CDTP::Constraint_id                                               Constraint_id;
+using Vb = CGAL::Base_with_time_stamp<CGAL::Triangulation_vertex_base_2<K>>;
-typedef CDTP::Vertices_in_constraint_iterator                             Vertices_in_constraint_iterator;
+using Cb = CGAL::Base_with_time_stamp<CGAL::Constrained_triangulation_face_base_2<K>>;
+using Tds = CGAL::Triangulation_data_structure_2<Vb, Cb>;
+using CDT = CGAL::Constrained_Delaunay_triangulation_2<K, Tds, Itag_>;
+using CDTP = CGAL::Constrained_triangulation_plus_2<CDT>;
 
-int countVertex(CDTP &cdtp, CDTP::Constraint_id id)
+using Point = CDTP::Point;
+using Vertex_handle = CDTP::Vertex_handle;
+using Constraint_id = CDTP::Constraint_id;
+using Vertices_in_constraint_iterator = CDTP::Vertices_in_constraint_iterator;
+
+auto nb_of_vertices(CDTP &cdtp, Constraint_id id)
 {
-    Vertices_in_constraint_iterator v=cdtp.vertices_in_constraint_begin(id);
+    return static_cast<std::size_t>(std::distance(cdtp.vertices_in_constraint_begin(id),
-
+                                                  cdtp.vertices_in_constraint_end(id)));
-    int count=0;
-    while(v!=cdtp.vertices_in_constraint_end(id))
-    {
-        count++;
-        v++;
-    }
-
-    return count;
 }
 
+auto value_check_expected =
+    [](auto&& value, [[maybe_unused]] const auto& expected) -> decltype(auto)
+    {
+        assert(value == expected);
+        return std::forward<decltype(value)>(value);
+    };
+
+auto oformat =
+    [](Vertex_handle vh)
+    {
+        return CGAL::IO::oformat(vh, CGAL::With_point_tag{});
+    };
 
 int main()
 {
     CDTP cdtp;
 
-    std::list<Point> pointsListCollinear;
+    auto print_cdtp = [&cdtp](std::string_view msg)
+    {
+        std::cout << msg << std::endl;
+        cdtp.print_hierarchy(std::cout);
+    };
 
-    pointsListCollinear.push_back(Point(0,0));
+    const std::array<Point, 6> collinear_points = {
-    pointsListCollinear.push_back(Point(0,1));
+        Point(0,0), Point(0,1), Point(0,2), Point(0,3), Point(0,4), Point(0,5)
-    pointsListCollinear.push_back(Point(0,2));
+    };
-    pointsListCollinear.push_back(Point(0,3));
-    pointsListCollinear.push_back(Point(0,4));
-    pointsListCollinear.push_back(Point(0,5));
 
-    std::list<Point> pointsListNoCollinear;
+    const std::array<Point, 6> non_collinear_points = {
+        Point(1,0), Point(2,1), Point(4,2), Point(2,3), Point(4,4), Point(1,5)
+    };
 
-    pointsListNoCollinear.push_back(Point(1,0));
+    Constraint_id collinear_cid = cdtp.insert_constraint(collinear_points.begin(), collinear_points.end());
-    pointsListNoCollinear.push_back(Point(2,1));
+    Constraint_id non_collinear_cid = cdtp.insert_constraint(non_collinear_points.begin(), non_collinear_points.end());
-    pointsListNoCollinear.push_back(Point(4,2));
-    pointsListNoCollinear.push_back(Point(2,3));
-    pointsListNoCollinear.push_back(Point(4,4));
-    pointsListNoCollinear.push_back(Point(1,5));
 
+    print_cdtp("Initial state");
 
-    Constraint_id ctIdCollinear=cdtp.insert_constraint(pointsListCollinear.begin(),pointsListCollinear.end());
+    Vertices_in_constraint_iterator vertex_it = std::next(cdtp.vertices_in_constraint_begin(collinear_cid), 2);
-    Constraint_id ctIdNoCollinear=cdtp.insert_constraint(pointsListNoCollinear.begin(),pointsListNoCollinear.end());
+    [[maybe_unused]] auto next_it = std::next(vertex_it);
+    std::cout << "\n-> attempt to remove vertex " << oformat(*vertex_it) << std::endl;
+    vertex_it = cdtp.remove_vertex_from_constraint(collinear_cid, vertex_it);
+    std::cout << "   cdtp.remove_vertex_from_constraint(collinear_cid, vertex_it) returned the vertex "
+              << oformat(*vertex_it) << std::endl;
+    assert(vertex_it == next_it);
 
+    print_cdtp("\nAfter removing third vertex from the collinear constraint");
 
-    //******************************* attempt with the collinear constraint
+    // The first constraint (ID `collinear_cid`) is collinear. `cdtp.remove_vertex_from_constraint`
-    Vertices_in_constraint_iterator vertexToRemoveCollinear=cdtp.vertices_in_constraint_begin(ctIdCollinear);
+    // cannot remove the third vertex from it, because it is collinear with the triangulation vertex
-    vertexToRemoveCollinear++;
+    // with the point (0, 2).
-    vertexToRemoveCollinear++;
 
+    std::cout << "\nnumber of subconstraints: "
+              << value_check_expected(cdtp.number_of_subconstraints(), 10U) << std::endl;
+    std::cout << "number of constraints: "
+              << value_check_expected(cdtp.number_of_constraints(), 2U) << std::endl;
+    std::cout << "number of vertex in collinear constraint: "
+              << value_check_expected(nb_of_vertices(cdtp, collinear_cid), 6U) << std::endl;
 
-        std::cout<<"attempt to remove vertex "<<(*vertexToRemoveCollinear)->point().x()<<" , "<<(*vertexToRemoveCollinear)->point().y() <<std::endl;
+    vertex_it = std::next(cdtp.vertices_in_constraint_begin(non_collinear_cid), 2);
-    cdtp.remove_vertex_from_constraint(ctIdCollinear,vertexToRemoveCollinear);
+    next_it = std::next(vertex_it);
+    std::cout << "\n-> attempt to remove vertex " << oformat(*vertex_it) << std::endl;
+    vertex_it = cdtp.remove_vertex_from_constraint(non_collinear_cid, vertex_it);
+    std::cout << "   cdtp.remove_vertex_from_constraint(non_collinear_cid, vertex_it) returned the vertex "
+              << oformat(*vertex_it) << std::endl;
+    assert(vertex_it == next_it);
 
-    std::cout<<"number of subconstraints "<<cdtp.number_of_subconstraints()<<std::endl; //--> 5, expected 4
+    print_cdtp("\nAfter removing third vertex from the non-collinear constraint");
-    std::cout<<"number of constraints "<<cdtp.number_of_constraints()<<std::endl; //--> 1
-    std::cout<<"number of vertex in constraint "<<countVertex(cdtp,ctIdCollinear)<<std::endl; //--> 6, expected 5
 
-
+    std::cout << "number of subconstraints: "
-    //******************************* attempt with the collinear constraint
+              << value_check_expected(cdtp.number_of_subconstraints(), 9U) << std::endl;
-    Vertices_in_constraint_iterator vertexToRemoveNoCollinear=cdtp.vertices_in_constraint_begin(ctIdNoCollinear);
+    std::cout << "number of constraints: "
-    vertexToRemoveNoCollinear++;
+              << value_check_expected(cdtp.number_of_constraints(), 2U) << std::endl;
-    vertexToRemoveNoCollinear++;
+    std::cout << "number of vertex in collinear constraint: "
-
+              << value_check_expected(nb_of_vertices(cdtp, non_collinear_cid), 5U) << std::endl;
-        std::cout<<"attempt to remove vertex "<<(*vertexToRemoveNoCollinear)->point().x()<<" , "<<(*vertexToRemoveNoCollinear)->point().y() << std::endl;
-    cdtp.remove_vertex_from_constraint(ctIdNoCollinear,vertexToRemoveNoCollinear);
-
-    std::cout<<"number of subconstraints "<<cdtp.number_of_subconstraints()<<std::endl; //--> 4, ok
-    std::cout<<"number of constraints "<<cdtp.number_of_constraints()<<std::endl;  //--> 1
-    std::cout<<"number of vertex in constraint "<<countVertex(cdtp,ctIdNoCollinear)<<std::endl;  //--> 5, ok
 
     return 0;
-
 }