parallel lloyd_optimize_mesh_2

Mesh_2/examples/Mesh_2/CMakeLists.txt

@@ -3,6 +3,7 @@ project(Mesh_2_Examples)
 
 find_package(CGAL REQUIRED)
 
+find_package(Threads)
 # create a target per cppfile
 file(
   GLOB cppfiles
@@ -11,3 +12,6 @@ file(
 foreach(cppfile ${cppfiles})
   create_single_source_cgal_program("${cppfile}")
 endforeach()
+
+target_link_libraries(test_parallel_Mesh_global_optimizer PRIVATE Threads::Threads)
+target_link_libraries(mesh_optimization PRIVATE Threads::Threads)

Mesh_2/examples/Mesh_2/test_parallel_Mesh_global_optimizer.cpp

@@ -0,0 +1,124 @@
+#define CGAL_MESH_2_OPTIMIZER_VERBOSE
+#define CGAL_MESH_2_OPTIMIZERS_DEBUG
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Constrained_Delaunay_triangulation_2.h>
+#include <CGAL/Delaunay_mesher_2.h>
+#include <CGAL/Delaunay_mesh_face_base_2.h>
+#include <CGAL/Delaunay_mesh_vertex_base_2.h>
+
+#include <CGAL/lloyd_optimize_mesh_2.h>
+#include <CGAL/point_generators_2.h>
+#include <CGAL/random_selection.h>
+
+#include <CGAL/draw_triangulation_2.h>
+#include <limits.h>
+
+
+#include <CGAL/Delaunay_triangulation_2.h>
+#include <CGAL/Delaunay_triangulation_adaptation_traits_2.h>
+#include <CGAL/Delaunay_triangulation_adaptation_policies_2.h>
+
+
+#include <iostream>
+
+
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+typedef CGAL::Delaunay_mesh_vertex_base_2<K>                Vb;
+typedef CGAL::Delaunay_mesh_face_base_2<K>                  Fb;
+typedef CGAL::Triangulation_data_structure_2<Vb, Fb>        Tds;
+typedef CGAL::Constrained_Delaunay_triangulation_2<K, Tds>  CDT;
+typedef CGAL::Delaunay_triangulation_2<K>                   DT;
+
+typedef CDT::Point Point;
+typedef CGAL::Creator_uniform_2<double,Point>            Creator;
+
+typedef CDT::Vertex_handle Vertex_handle;
+typedef CDT::Vertex_iterator Vertex_iterator;
+
+#define RADIUS 5.5
+#define RADIUS_BIRTH 5
+#define RADIUS_DEATH 5
+
+
+void cell_death(CDT &cdt)
+{
+  CDT::Finite_vertices_iterator it;
+  // CDT::Vertex_handle v; // equivalent à l'iterateur
+
+  //  std::vector<CDT::Vertex_handle> vector_vertices;
+  std::vector<CDT::Vertex_handle> vector_vertices;
+  std::vector<CDT::Vertex_handle> dead_cells;
+  std::vector<CDT::Vertex_handle>::iterator it_cells;
+  int count=0;
+  for (it = cdt.finite_vertices_begin();
+  	   it!= cdt.finite_vertices_end();
+  	   ++it)
+  	{
+	  // on ne garde que les points à l'intérieur des contraintes
+	  // delete only points without constraints
+	  if (!cdt.are_there_incident_constraints(it))
+	    {
+	      vector_vertices.push_back(it);
+	      count++;
+	    }
+  	}
+
+  int n =10;
+
+  random_selection(vector_vertices.begin(), vector_vertices.end(), n,
+		   std::back_inserter(dead_cells));
+
+
+  for (it_cells = dead_cells.begin(); it_cells!= dead_cells.end(); ++it_cells)
+    cdt.remove(*it_cells);
+
+}
+
+void cell_birth(CDT &cdt)
+{
+  // insert random points
+  CGAL::Random_points_in_disc_2<Point,Creator> g(RADIUS_BIRTH);
+  int new_cells = 10;
+  for (int i=0; i<new_cells; ++i)
+    {
+      cdt.insert(*g++);
+    }
+
+}
+
+int startSimu(int nb_points, int nb_iter)
+{
+
+  CDT cdt;
+  // insert random points at first
+  CGAL::Random_points_in_disc_2<Point,Creator> g(RADIUS);
+  for (int i=0; i<nb_points; ++i)
+    {
+      cdt.insert(*g++);
+    }
+
+
+  // loop on birth death and lloyd optimization
+  for (int n=0; n<nb_iter; ++n)
+    {
+      lloyd_optimize_mesh_2(cdt,
+      			    CGAL::parameters::max_iteration_number = nb_iter,
+      			    CGAL::parameters::mark = false);
+
+    }
+
+
+  return 0;
+}
+
+/// main function
+int main(int argc, char**argv)
+{
+  // arguments
+  int nb_points = 20; // number of points
+  int nb_iter=5; // nb of iterations
+
+
+  return startSimu(nb_points, nb_iter);
+}

Mesh_2/include/CGAL/Mesh_2/Lloyd_move_2.h

@@ -20,6 +20,8 @@
 #include <CGAL/Mesh_2/Uniform_sizing_field_2.h>
 #include <CGAL/Constrained_voronoi_diagram_2.h>
 
+#include <thread>
+
 namespace CGAL
 {
 namespace Mesh_2

Mesh_2/include/CGAL/Mesh_2/Mesh_global_optimizer_2.h

@@ -38,6 +38,9 @@
 #include <boost/format.hpp>
 #include <boost/math/constants/constants.hpp>
 
+#include <thread>
+#include <future>
+
 namespace CGAL {
 
 namespace Mesh_2 {
@@ -147,7 +150,9 @@ class Mesh_global_optimizer_2
       this->before_move();
 
       // Compute move for each vertex
-      Moves_vector moves = compute_moves(moving_vertices);
+      // Previous single thread version
+      //Moves_vector moves = compute_moves(moving_vertices);
+      Moves_vector moves = compute_moves_threaded(moving_vertices);
 
       //Pb with Freeze : sometimes a few vertices continue moving indefinitely
       //if the nb of moving vertices is < 2% of total nb AND does not decrease
@@ -221,6 +226,126 @@ class Mesh_global_optimizer_2
   }
 
 private:
+  // just for simplify notations
+  typedef typename Vertex_set::const_iterator VERTEXSET_ITERATOR;
+
+  std::mutex mtx;           // mutex for critical section
+
+  /// computes move for each thread
+  ///
+  /// tid: thread id
+  /// begin, end: iterator for locating vertices on the set. Deletions are
+  /// critical while computing the moves
+  /// moving_vertices: set of vertices that may move
+  Moves_vector threaded_move(int tid,
+			     VERTEXSET_ITERATOR begin, VERTEXSET_ITERATOR end,
+			     Vertex_set &moving_vertices)
+  {
+    typename Gt::Construct_translated_point_2 translate =
+      Gt().construct_translated_point_2_object();
+
+    Moves_vector moves;
+    // Get move for each moving vertex
+    VERTEXSET_ITERATOR vit;
+    for ( vit = begin ;
+	  vit != end ;
+	  )
+      {
+	Vertex_handle oldv = *vit;
+
+	Vector_2 move = compute_move(oldv);
+	vit++;
+
+	if ( CGAL::NULL_VECTOR != move )
+	  {
+
+	    Point_2 new_position = translate(oldv->point(), move);
+
+	    moves.push_back(std::make_pair(oldv, new_position));
+	  }
+	else if(sq_freeze_ratio_ > 0.) //freezing ON
+	  moving_vertices.erase(oldv);
+
+	// Stop if time_limit_ is reached
+	if ( is_time_limit_reached() )
+	  break;
+      }
+    return moves;
+  }
+
+  /**
+   * Split the vertices into subsets in order to dispatch the computation to multiple threads
+   *
+   */
+  Moves_vector compute_moves_threaded(Vertex_set moving_vertices)
+  {
+
+    // hard coded for testing purposes
+    int n_threads = 2;
+
+    typedef std::future<Moves_vector > TYPE_THREAD;
+    std::vector<TYPE_THREAD> vectorOfThreads;
+
+    typename Gt::Construct_translated_point_2 translate =
+      Gt().construct_translated_point_2_object();
+
+    // Store new location of points which have to move
+    Moves_vector moves;
+    moves.reserve(moving_vertices.size());
+
+    // reset worst_move list
+    std::fill(big_moves_.begin(), big_moves_.end(), FT(0));
+
+    // nb of vertices by thread
+    int n_vertices = moving_vertices.size();
+
+    int n_vertices_by_thread=0;
+    if (n_threads < n_vertices)
+      n_vertices_by_thread = n_vertices / n_threads + 1;
+    else
+      {
+	n_vertices_by_thread = n_vertices;
+	n_threads = 1;
+      }
+
+    // split vertices into different threads
+    for (int t=0; t<n_threads; ++t)
+      {
+	int begin = t*n_vertices_by_thread;
+	int end   = begin + n_vertices_by_thread;
+	if (end >= n_vertices)
+	  end = n_vertices;
+
+	VERTEXSET_ITERATOR itBegin = std::next(moving_vertices.begin(), begin);
+	VERTEXSET_ITERATOR itEnd   = std::next(moving_vertices.begin(), end);
+	if (end == n_vertices)
+	  {
+	    itEnd = moving_vertices.end();
+	  }
+
+	// create and start thread
+	TYPE_THREAD th = std::async(&Mesh_global_optimizer_2::threaded_move, this,
+				    t, std::ref(itBegin), std::ref(itEnd),
+				    std::ref(moving_vertices));
+
+	// store the thread for later use
+	vectorOfThreads.push_back(std::move(th));
+      }
+
+    // Get results in random order, wait for each thread to complete
+    for (auto &th : vectorOfThreads )
+      {
+    	Moves_vector moves_th = th.get(); // wait
+
+    	// agregation
+    	moves.insert(moves.end(), moves_th.begin(), moves_th.end());
+
+      }
+
+    return moves;
+  }
+
+
   /**
    * Returns moves for vertices of set \c moving_vertices
    */