Mesh_2: make it deterministic

Mesh_2/include/CGAL/Delaunay_mesh_face_base_2.h

@@ -17,6 +17,7 @@
 
 
 #include <CGAL/Constrained_Delaunay_triangulation_face_base_2.h>
+#include <CGAL/Has_timestamp.h>
 
 namespace CGAL {
 
@@ -67,6 +68,14 @@ public:
   /** compatibility with CGAL-3.2 */
   inline
   void set_marked(const bool b) { in_domain=b; }
+
+  typedef Tag_true Has_timestamp;
+
+  std::size_t time_stamp() const { return time_stamp_; }
+
+  void set_time_stamp(const std::size_t& ts) { time_stamp_ = ts; }
+
+  std::size_t time_stamp_;
 };
 
 } // namespace CGAL

Mesh_2/test/Mesh_2/reproductibility.cpp

@@ -0,0 +1,110 @@
+//#define CGAL_MESH_2_DEBUG_BAD_FACES
+//#define CGAL_MESH_2_DEBUG_CLUSTERS
+
+#include <CGAL/Constrained_Delaunay_triangulation_2.h>
+#include <CGAL/Delaunay_mesh_size_criteria_2.h>
+#include <CGAL/Delaunay_mesh_face_base_2.h>
+#include <CGAL/Delaunay_mesher_2.h>
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+#include <fstream>
+#include <sstream>
+#include <iostream>
+#include <string>
+
+using K = CGAL::Exact_predicates_inexact_constructions_kernel;
+using Vb = CGAL::Triangulation_vertex_base_2<K>;
+using Fb = CGAL::Delaunay_mesh_face_base_2<K>;
+using Tds = CGAL::Triangulation_data_structure_2<Vb, Fb>;
+using CDT = CGAL::Constrained_Delaunay_triangulation_2<K, Tds, CGAL::Exact_predicates_tag>;
+using Criteria = CGAL::Delaunay_mesh_size_criteria_2<CDT>;
+using Mesher =  CGAL::Delaunay_mesher_2<CDT, Criteria>;
+
+using Vertex_handle = CDT::Vertex_handle;
+using Point = CDT::Point;
+
+int main(int argc, char* argv[])
+{
+  std::string path = argv[0];
+  path = path.substr(0, path.rfind('/') + 1);
+
+  std::cout << "Current dir:" << path << std::endl;
+
+  auto triangulate = [&path](int index)
+  {
+    CDT cdt;
+
+    auto write_tr = [&](const std::string& filename)
+    {
+//       std::ofstream file(path + filename + "_" + std::to_string(index) + ".off");
+//
+//       cdt.file_output(file);
+//       file.close();
+    };
+
+    Vertex_handle va = cdt.insert(Point(-0.74397572, -0.54545455));
+    Vertex_handle vb = cdt.insert(Point(-0.13526831, -1));
+    Vertex_handle vc = cdt.insert(Point(0.067634156, -1));
+    Vertex_handle vd = cdt.insert(Point(0.33817078, -0.54545455));
+    Vertex_handle ve = cdt.insert(Point(0.74397572, 0.27272727));
+    Vertex_handle vf = cdt.insert(Point(0.74397572, 0.54545455));
+    Vertex_handle vg = cdt.insert(Point(0.067634156, 1));
+    Vertex_handle vh = cdt.insert(Point(-0.13526831, 1));
+    Vertex_handle vi = cdt.insert(Point(-0.74397572, -0.18181818));
+
+    cdt.insert_constraint(va, vb);
+    cdt.insert_constraint(vb, vc);
+    cdt.insert_constraint(vc, vd);
+    cdt.insert_constraint(vd, ve);
+    cdt.insert_constraint(ve, vf);
+    cdt.insert_constraint(vf, vg);
+    cdt.insert_constraint(vg, vh);
+    cdt.insert_constraint(vh, vi);
+    cdt.insert_constraint(vi, va);
+
+    const std::vector<Point> points{
+      Point(0.65605132, 0.43821259),
+      Point(0.23073753, -0.4476739),
+      Point(-0.037496007, -0.93636364),
+      Point(-0.00095596601, 0.88181818),
+      Point(-0.62452925, -0.30720903),
+      Point(-0.69663181, -0.45045525),
+    };
+
+    cdt.insert(points.cbegin(), points.cend());
+
+    std::cout << "Meshing: " << index << std::endl;
+
+    std::cout << "Number of vertices before: " << cdt.number_of_vertices() << std::endl;
+
+    write_tr("before_refine");
+
+    Mesher mesher(cdt);
+    mesher.set_criteria(Criteria(0.125, 0.05*std::sqrt(2)));
+
+//     mesher.clear_seeds();
+//     mesher.init();
+
+    mesher.refine_mesh();
+
+    write_tr("after_refine");
+
+    std::cout << "Number of vertices after: " << cdt.number_of_vertices() << std::endl;
+
+    std::stringstream ss;
+    ss << cdt;
+
+    return ss.str();
+  };
+
+  const std::string ref_cdts = triangulate(0);
+
+  for (int i = 1; i < 20; ++i)
+  {
+    const std::string cdts = triangulate(i);
+    if (ref_cdts != cdts)
+      return 1;
+  }
+
+  return 0;
+}

TDS_2/include/CGAL/Triangulation_data_structure_2.h

@@ -1032,7 +1032,9 @@ insert_dim_up(Vertex_handle w,  bool orient)
 
       for ( ; lfit != faces_list.end() ; ++lfit) {
         f = * lfit;
-        g = create_face(f); //calls copy constructor of face
+        g = create_face(f->vertex(0),f->vertex(1),f->vertex(2),
+                        f->neighbor(0),f->neighbor(1),f->neighbor(2));
+
         f->set_vertex(dim,v);
         g->set_vertex(dim,w);
         set_adjacency(f, dim, g, dim);

Triangulation_2/include/CGAL/Constrained_triangulation_2.h

@@ -148,6 +148,7 @@ public:
   using Triangulation::all_edges_begin;
   using Triangulation::all_edges_end;
   using Triangulation::mirror_index;
+  using Triangulation::mirror_edge;
   using Triangulation::orientation;
 #endif
 
@@ -656,6 +657,19 @@ insert(const Point& a, Locate_type lt, Face_handle loc, int li)
   Vertex_handle v1, v2;
   bool insert_in_constrained_edge = false;
 
+  std::list<std::pair<Vertex_handle,Vertex_handle> > constrained_edges;
+  bool one_dimensional = false;
+  if(dimension() == 1){
+    one_dimensional = true;
+    for(Finite_edges_iterator it = finite_edges_begin();
+        it != finite_edges_end();
+        ++it){
+      if(is_constrained(*it)){
+        constrained_edges.push_back(std::make_pair(it->first->vertex(cw(it->second)),
+                                                   it->first->vertex(ccw(it->second))));
+      }
+    }
+  }
   if ( lt == Triangulation::EDGE && loc->is_constrained(li) )
   {
     if(boost::is_same<Itag, No_constraint_intersection_tag>::value)
@@ -668,6 +682,18 @@ insert(const Point& a, Locate_type lt, Face_handle loc, int li)
 
   va = Triangulation::insert(a,lt,loc,li);
 
+  if(one_dimensional && (dimension() == 2)){
+    for(const std::pair<Vertex_handle,Vertex_handle> vp : constrained_edges){
+      Face_handle fh;
+      int i;
+      if(this->is_edge(vp.first, vp.second, fh,i)){
+        fh->set_constraint(i,true);
+        boost::tie(fh,i) = mirror_edge(Edge(fh,i));
+        fh->set_constraint(i,true);
+      }
+    }
+  }
+
   if (insert_in_constrained_edge)
     update_constraints_incident(va, v1,v2);
   else if(lt != Triangulation::VERTEX)