diff --git a/Triangulation_3/include/CGAL/Delaunay_triangulation_3.h b/Triangulation_3/include/CGAL/Delaunay_triangulation_3.h index e9afaef555d..5b9c2449d6d 100644 --- a/Triangulation_3/include/CGAL/Delaunay_triangulation_3.h +++ b/Triangulation_3/include/CGAL/Delaunay_triangulation_3.h @@ -413,13 +413,94 @@ private: spatial_sort(indices.begin(),indices.end(),Search_traits(&(points[0]),geom_traits())); - Vertex_handle hint; - for (typename std::vector::const_iterator - it = indices.begin(), end = indices.end(); - it != end; ++it){ - hint = insert(points[*it], hint); - if (hint!=Vertex_handle()) hint->info()=infos[*it]; - } +#ifdef CGAL_LINKED_WITH_TBB + if(this->is_parallel()){ + + size_t num_points = points.size(); + Vertex_handle hint; + std::vector far_sphere_vertices; + +#ifdef CGAL_CONCURRENT_TRIANGULATION_3_ADD_TEMPORARY_POINTS_ON_FAR_SPHERE + const size_t MIN_NUM_POINTS_FOR_FAR_SPHERE_POINTS = 1000000; + if (num_points >= MIN_NUM_POINTS_FOR_FAR_SPHERE_POINTS) + { + // Add temporary vertices on a "far sphere" to reduce contention on + // the infinite vertex + + // Get bbox + const Bbox_3 &bbox = *this->get_bbox(); + // Compute radius for far sphere + const double& xdelta = bbox.xmax() - bbox.xmin(); + const double& ydelta = bbox.ymax() - bbox.ymin(); + const double& zdelta = bbox.zmax() - bbox.zmin(); + const double radius = 1.3 * 0.5 * std::sqrt(xdelta*xdelta + + ydelta*ydelta + + zdelta*zdelta); + // WARNING - TODO: this code has to be fixed because Vector_3 is not + // required by the traits concept + const typename Gt::Vector_3 center( + bbox.xmin() + 0.5*xdelta, + bbox.ymin() + 0.5*ydelta, + bbox.zmin() + 0.5*zdelta); + Random_points_on_sphere_3 random_point(radius); + const int NUM_PSEUDO_INFINITE_VERTICES = static_cast( + tbb::task_scheduler_init::default_num_threads() * 3.5); + std::vector points_on_far_sphere; + for (int i = 0 ; i < NUM_PSEUDO_INFINITE_VERTICES ; ++i, ++random_point) + points_on_far_sphere.push_back(*random_point + center); + + spatial_sort(points_on_far_sphere.begin(), + points_on_far_sphere.end(), + geom_traits()); + + std::vector::const_iterator it_p = points_on_far_sphere.begin(); + std::vector::const_iterator it_p_end = points_on_far_sphere.end(); + for ( ; it_p != it_p_end ; ++it_p) + { + hint = insert(*it_p, hint); + far_sphere_vertices.push_back(hint); + } + } +#endif // CGAL_CONCURRENT_TRIANGULATION_3_ADD_TEMPORARY_POINTS_ON_FAR_SPHERE + + size_t i = 0; + // Insert "num_points_seq" points sequentially + // (or more if dim < 3 after that) + size_t num_points_seq = (std::min)(num_points, (size_t)100); + while (dimension() < 3 || i < num_points_seq) + { + hint = insert(points[indices[i]], hint); + if (hint != Vertex_handle()) hint->info() = infos[indices[i]]; + ++i; + } + + tbb::enumerable_thread_specific tls_hint(hint); + tbb::parallel_for( + tbb::blocked_range( i, num_points ), + Insert_point_with_info(*this, points, infos, indices, tls_hint) + ); + +#ifdef CGAL_CONCURRENT_TRIANGULATION_3_ADD_TEMPORARY_POINTS_ON_FAR_SPHERE + if (num_points >= MIN_NUM_POINTS_FOR_FAR_SPHERE_POINTS) + { + // Remove the temporary vertices on far sphere + remove(far_sphere_vertices.begin(), far_sphere_vertices.end()); + } +#endif // CGAL_CONCURRENT_TRIANGULATION_3_ADD_TEMPORARY_POINTS_ON_FAR_SPHERE + + } + // Sequential + else +#endif + { + Vertex_handle hint; + for (typename std::vector::const_iterator + it = indices.begin(), end = indices.end(); + it != end; ++it) { + hint = insert(points[*it], hint); + if (hint != Vertex_handle()) hint->info() = infos[*it]; + } + } return number_of_vertices() - n; } @@ -909,6 +990,88 @@ protected: { m_dt.unlock_all_elements(); +#ifdef CGAL_CONCURRENT_TRIANGULATION_3_PROFILING + bcounter.increment_branch_2(); // THIS is an early withdrawal! +#endif + } + } + + } + } + }; + + // Functor for parallel insert(begin, end) function with info + template + class Insert_point_with_info + { + typedef typename DT::Point Point; + typedef typename DT::Triangulation_data_structure::Vertex::Info Info; + typedef typename DT::Vertex_handle Vertex_handle; + + DT & m_dt; + const std::vector & m_points; + const std::vector & m_infos; + const std::vector & m_indices; + tbb::enumerable_thread_specific & m_tls_hint; + + public: + // Constructor + Insert_point_with_info(DT & dt, + const std::vector & points, + const std::vector & infos, + const std::vector & indices, + tbb::enumerable_thread_specific & tls_hint) + : m_dt(dt), m_points(points), m_infos(infos), m_indices(indices), m_tls_hint(tls_hint) + {} + + // Constructor + Insert_point_with_info(const Insert_point_with_info &ip) + : m_dt(ip.m_dt), m_points(ip.m_points), m_infos(ip.m_infos), m_indices(ip.m_indices), m_tls_hint(ip.m_tls_hint) + {} + + // operator() + void operator()( const tbb::blocked_range& r ) const + { +#ifdef CGAL_CONCURRENT_TRIANGULATION_3_PROFILING + static Profile_branch_counter_3 bcounter( + "early withdrawals / late withdrawals / successes [Delaunay_tri_3::insert]"); +#endif + + Vertex_handle &hint = m_tls_hint.local(); + for( std::size_t i_idx = r.begin() ; i_idx != r.end() ; ++i_idx) + { + bool success = false; + std::ptrdiff_t i_point = m_indices[i_idx]; + while(!success) + { + if (m_dt.try_lock_vertex(hint) && m_dt.try_lock_point(m_points[i_point])) + { + bool could_lock_zone; + Vertex_handle new_hint = m_dt.insert( + m_points[i_point], hint, &could_lock_zone); + + m_dt.unlock_all_elements(); + + if (could_lock_zone) + { + hint = new_hint; + if (hint!=Vertex_handle()) hint->info()=m_infos[i_point]; + success = true; +#ifdef CGAL_CONCURRENT_TRIANGULATION_3_PROFILING + ++bcounter; +#endif + } +#ifdef CGAL_CONCURRENT_TRIANGULATION_3_PROFILING + else + { + bcounter.increment_branch_1(); // THIS is a late withdrawal! + } +#endif + } + else + { + m_dt.unlock_all_elements(); + #ifdef CGAL_CONCURRENT_TRIANGULATION_3_PROFILING bcounter.increment_branch_2(); // THIS is an early withdrawal! #endif