// Copyright (c) 2009 INRIA Sophia-Antipolis (France). // All rights reserved. // // This file is part of CGAL (www.cgal.org). // // $URL: https://github.com/CGAL/cgal/blob/v5.2/Mesh_3/include/CGAL/Mesh_3/Mesh_global_optimizer.h $ // $Id: Mesh_global_optimizer.h 4fc2f59 2020-07-31T16:17:56+02:00 Laurent Rineau // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial // // // Author(s) : Stephane Tayeb // //****************************************************************************** // File Description : //****************************************************************************** #ifndef CGAL_MESH_3_MESH_GLOBAL_OPTIMIZER_H #define CGAL_MESH_3_MESH_GLOBAL_OPTIMIZER_H #include #include #include #include #include #include #include #include #include #include #include #include #include #ifdef CGAL_MESH_3_PROFILING #include #endif #include #include #include #include #ifdef CGAL_LINKED_WITH_TBB # include # include # include # include #endif namespace CGAL { namespace Mesh_3 { /************************************************ // Class Mesh_global_optimizer_base // Two versions: sequential / parallel ************************************************/ // Sequential template class Mesh_global_optimizer_base { protected: typedef typename Tr::Geom_traits Gt; typedef typename Gt::FT FT; typedef typename Gt::Vector_3 Vector_3; typedef typename Tr::Lock_data_structure Lock_data_structure; // The sizing field info is stored inside the move vector because it is computed // when the move is computed. This is because the parallel version uses the threadsafe // version of incident_cells (which thus requires points to not be moving yet) typedef std::vector > Moves_vector; typedef unsigned int Nb_frozen_points_type; Mesh_global_optimizer_base(const Bbox_3 &, int) : big_moves_size_(0) {} void update_big_moves(const FT& new_sq_move) { if (big_moves_.size() < big_moves_size_ ) big_moves_.insert(new_sq_move); else { FT smallest = *(big_moves_.begin()); if( new_sq_move > smallest ) { big_moves_.erase(big_moves_.begin()); big_moves_.insert(new_sq_move); } } } void clear_big_moves() { big_moves_.clear(); } Lock_data_structure *get_lock_data_structure() { return 0; } void unlock_all_elements() {} protected: std::size_t big_moves_size_; std::multiset big_moves_; }; #ifdef CGAL_LINKED_WITH_TBB // Parallel template class Mesh_global_optimizer_base { protected: typedef typename Tr::Geom_traits Gt; typedef typename Gt::FT FT; typedef typename Gt::Vector_3 Vector_3; typedef typename Tr::Lock_data_structure Lock_data_structure; typedef tbb::concurrent_vector > Moves_vector; typedef std::atomic Nb_frozen_points_type ; Mesh_global_optimizer_base(const Bbox_3 &bbox, int num_grid_cells_per_axis) : big_moves_size_(0) , m_lock_ds(bbox, num_grid_cells_per_axis) { big_moves_current_size_ = 0; big_moves_smallest_ = (std::numeric_limits::max)(); } void update_big_moves(const FT& new_sq_move) { if (++big_moves_current_size_ <= big_moves_size_ ) { std::lock_guard lock(m_big_moves_mutex); typename std::multiset::const_iterator it = big_moves_.insert(new_sq_move); // New smallest move of all big moves? if (it == big_moves_.begin()) big_moves_smallest_ = new_sq_move; } else { --big_moves_current_size_; if( new_sq_move > big_moves_smallest_ ) { std::lock_guard lock(m_big_moves_mutex); // Test it again since it may have been modified by another // thread in the meantime if( new_sq_move > big_moves_smallest_ ) { big_moves_.erase(big_moves_.begin()); typename std::multiset::const_iterator it = big_moves_.insert(new_sq_move); // New smallest move of all big moves? if (it == big_moves_.begin()) big_moves_smallest_ = new_sq_move; } } } } void clear_big_moves() { big_moves_current_size_ = 0; big_moves_smallest_ = (std::numeric_limits::max)(); big_moves_.clear(); } Lock_data_structure *get_lock_data_structure() { return &m_lock_ds; } void unlock_all_elements() { m_lock_ds.unlock_all_points_locked_by_this_thread(); } public: protected: std::atomic big_moves_current_size_; std::atomic big_moves_smallest_; std::size_t big_moves_size_; std::multiset big_moves_; std::mutex m_big_moves_mutex; /// Lock data structure Lock_data_structure m_lock_ds; }; #endif // CGAL_LINKED_WITH_TBB /************************************************ // Class Mesh_global_optimizer ************************************************/ template > class Mesh_global_optimizer : public Mesh_global_optimizer_base { // Types typedef typename C3T3::Concurrency_tag Concurrency_tag; typedef Mesh_global_optimizer Self; typedef Mesh_global_optimizer_base< typename C3T3::Triangulation, typename C3T3::Concurrency_tag> Base; using Base::get_lock_data_structure; using Base::big_moves_; using Base::big_moves_size_; typedef typename C3T3::Triangulation Tr; typedef typename Tr::Geom_traits Gt; typedef typename Tr::Bare_point Bare_point; typedef typename Tr::Weighted_point Weighted_point; typedef typename Tr::Cell_handle Cell_handle; typedef typename Tr::Vertex_handle Vertex_handle; typedef typename Tr::Edge Edge; typedef typename Tr::Vertex Vertex; typedef typename Gt::FT FT; typedef typename Gt::Vector_3 Vector_3; typedef typename std::vector Cell_vector; typedef typename std::vector Vertex_vector; typedef Hash_handles_with_or_without_timestamps Hash_fct; typedef typename boost::unordered_set Vertex_set; typedef typename Base::Moves_vector Moves_vector; typedef typename Base::Nb_frozen_points_type Nb_frozen_points_type; #ifdef CGAL_INTRUSIVE_LIST typedef Intrusive_list Outdated_cell_set; #else typedef std::set Outdated_cell_set; #endif //CGAL_INTRUSIVE_LIST #ifdef CGAL_INTRUSIVE_LIST typedef Intrusive_list Moving_vertices_set; #else typedef Vertex_set Moving_vertices_set; #endif typedef typename MoveFunction::Sizing_field Sizing_field; typedef class C3T3_helpers C3T3_helpers; // Visitor class // Should define: // - after_compute_moves() // - after_move_points() // - after_rebuild_restricted_delaunay() // - end_of_iteration(int iteration_number) typedef Visitor_ Visitor; public: /** * Constructor */ Mesh_global_optimizer(C3T3& c3t3, const MeshDomain& domain, const FT& freeze_ratio, const bool do_freeze, const FT& convergence_ratio, const MoveFunction move_function = MoveFunction()); /** * Launch optimization process * * @param nb_interations maximum number of iterations */ Mesh_optimization_return_code operator()(int nb_iterations, Visitor v = Visitor()); /** * collects all vertices of the triangulation in moving_vertices * (even the frozen ones) */ void collect_all_vertices(Moving_vertices_set& moving_vertices); /// Time accessors void set_time_limit(double time) { time_limit_ = time; } double time_limit() const { return time_limit_; } private: /** * Returns moves for vertices of set \c moving_vertices */ Moves_vector compute_moves(Moving_vertices_set& moving_vertices); /** * Returns the move for vertex \c v * \warning This function should be called only on moving vertices * even for frozen vertices, it could return a non-zero vector */ Vector_3 compute_move(const Vertex_handle& v); /** * Updates mesh using moves of \c moves vector. Updates moving_vertices with * the new set of moving vertices after the move. */ void update_mesh(const Moves_vector& moves, Moving_vertices_set& moving_vertices, Visitor& visitor); /** * Fill sizing field using sizes (avg circumradius) contained in tr_ */ void fill_sizing_field(); /** * Returns `true` if convergence is reached */ bool check_convergence() const; /** * Returns the average circumradius length of cells incident to \c v */ FT average_circumradius_length(const Vertex_handle& v) const; /** * Returns the minimum cicumradius length of cells incident to \c v */ FT min_circumradius_sq_length(const Vertex_handle& v, const Cell_vector& incident_cells) const; /** * Returns the squared circumradius length of cell \c cell */ FT sq_circumradius_length(const Cell_handle& cell, const Vertex_handle& v) const; /** * Returns `true` if time_limit is reached */ bool is_time_limit_reached() const { return ( (time_limit() > 0) && (running_time_.time() > time_limit()) ); } private: #ifdef CGAL_LINKED_WITH_TBB // Functor for compute_moves function template class Compute_move { typedef tbb::concurrent_vector Vertex_conc_vector; MGO & m_mgo; const Sizing_field_ & m_sizing_field; Moves_vector_ & m_moves; bool m_do_freeze; Vertex_conc_vector & m_vertices_not_moving_any_more; const Gt & m_gt; public: // Constructor Compute_move(MGO &mgo, const Sizing_field_ &sizing_field, Moves_vector_ &moves, bool do_freeze, Vertex_conc_vector &vertices_not_moving_any_more, const Gt >) : m_mgo(mgo), m_sizing_field(sizing_field), m_moves(moves), m_do_freeze(do_freeze), m_vertices_not_moving_any_more(vertices_not_moving_any_more), m_gt(gt) {} // Constructor Compute_move(const Compute_move &cm) : m_mgo(cm.m_mgo), m_sizing_field(cm.m_sizing_field), m_moves(cm.m_moves), m_do_freeze(cm.m_do_freeze), m_vertices_not_moving_any_more(cm.m_vertices_not_moving_any_more), m_gt(cm.m_gt) {} // operator() void operator()(const Vertex_handle& oldv) const { typename Gt::Construct_point_3 cp = m_gt.construct_point_3_object(); typename Gt::Construct_translated_point_3 translate = m_gt.construct_translated_point_3_object(); Vector_3 move = m_mgo.compute_move(oldv); if ( CGAL::NULL_VECTOR != move ) { FT size = 0.; if(MGO::Sizing_field::is_vertex_update_needed) { const Weighted_point& position = m_mgo.tr_.point(oldv); Bare_point new_position = translate(cp(position), move); size = m_sizing_field(new_position, oldv); } // typedef Triangulation_helpers Th; //if( !Th().inside_protecting_balls(tr_, oldv, new_position)) //note : this is not happening for Lloyd and ODT so it's commented // maybe for a new global optimizer it should be de-commented m_moves.push_back(std::make_tuple(oldv, move, size)); } else // CGAL::NULL_VECTOR == move { if(m_do_freeze) { m_vertices_not_moving_any_more.push_back(oldv); } } if ( m_mgo.is_time_limit_reached() ) m_mgo.cancel_group_execution(); } }; // Functor for fill_sizing_field function template class Compute_sizing_field_value { MGO & m_mgo; const Gt & m_gt; Local_list_ & m_local_lists; public: // Constructor Compute_sizing_field_value(MGO &mgo, const Gt >, Local_list_ &local_lists) : m_mgo(mgo), m_gt(gt), m_local_lists(local_lists) {} // Constructor Compute_sizing_field_value(const Compute_sizing_field_value &csfv) : m_mgo(csfv.m_mgo), m_gt(csfv.m_gt), m_local_lists(csfv.m_local_lists) {} // operator() void operator()(Vertex& v) const { typename Gt::Construct_point_3 cp = m_gt.construct_point_3_object(); Vertex_handle vh = Tr::Triangulation_data_structure::Vertex_range::s_iterator_to(v); const Weighted_point& position = m_mgo.tr_.point(vh); m_local_lists.local().push_back( std::make_pair(cp(position), m_mgo.average_circumradius_length(vh))); } }; // Functor for update_mesh function template class Move_vertex { MGO & m_mgo; const Helper & m_helper; const Moves_vector_ & m_moves; Moving_vertices_set_ & m_moving_vertices; Outdated_cell_set_ & m_outdated_cells; typedef typename Tr_::Bare_point Bare_point; typedef typename Tr_::Vertex_handle Vertex_handle; public: // Constructor Move_vertex(MGO &mgo, const Helper &helper, const Moves_vector_ &moves, Moving_vertices_set_ &moving_vertices, Outdated_cell_set_ &outdated_cells) : m_mgo(mgo), m_helper(helper), m_moves(moves), m_moving_vertices(moving_vertices), m_outdated_cells(outdated_cells) {} // Constructor Move_vertex(const Move_vertex &mv) : m_mgo(mv.m_mgo), m_helper(mv.m_helper), m_moves(mv.m_moves), m_moving_vertices(mv.m_moving_vertices), m_outdated_cells(mv.m_outdated_cells) {} // operator() void operator()( const tbb::blocked_range& r ) const { for( size_t i = r.begin() ; i != r.end() ; ++i) { const Vertex_handle& v = std::get<0>(m_moves[i]); const Vector_3& move = std::get<1>(m_moves[i]); // Get size at new position if ( MGO::Sizing_field::is_vertex_update_needed ) { FT size = std::get<2>(m_moves[i]); // Move point bool could_lock_zone; Vertex_handle new_v = m_helper.move_point( v, move, m_outdated_cells, m_moving_vertices, &could_lock_zone); while (could_lock_zone == false) { new_v = m_helper.move_point( v, move, m_outdated_cells, m_moving_vertices, &could_lock_zone); } // Restore size in meshing_info data new_v->set_meshing_info(size); } else // Move point { bool could_lock_zone; do { m_helper.move_point( v, move, m_outdated_cells, m_moving_vertices, &could_lock_zone); } while (!could_lock_zone); } m_mgo.unlock_all_elements(); // Stop if time_limit_ is reached, here we can't return without rebuilding // restricted delaunay if ( m_mgo.is_time_limit_reached() ) { m_mgo.cancel_group_execution(); break; } } } }; void cancel_group_execution() { tbb_task_group_context.cancel_group_execution(); } #endif // CGAL_LINKED_WITH_TBB // ----------------------------------- // Private data // ----------------------------------- C3T3& c3t3_; Tr& tr_; const MeshDomain& domain_; FT sq_freeze_ratio_; FT convergence_ratio_; C3T3_helpers helper_; MoveFunction move_function_; Sizing_field sizing_field_; double time_limit_; CGAL::Real_timer running_time_; bool do_freeze_; mutable Nb_frozen_points_type nb_frozen_points_; #ifdef CGAL_MESH_3_OPTIMIZER_VERBOSE mutable FT sum_moves_; #endif #ifdef CGAL_LINKED_WITH_TBB tbb::task_group_context tbb_task_group_context; #endif }; template Mesh_global_optimizer:: Mesh_global_optimizer(C3T3& c3t3, const Md& domain, const FT& freeze_ratio, const bool do_freeze, const FT& convergence_ratio, const Mf move_function) : Base(c3t3.bbox(), Concurrent_mesher_config::get().locking_grid_num_cells_per_axis) , c3t3_(c3t3) , tr_(c3t3_.triangulation()) , domain_(domain) , sq_freeze_ratio_(freeze_ratio*freeze_ratio) , convergence_ratio_(convergence_ratio) , helper_(c3t3_,domain_, get_lock_data_structure()) , move_function_(move_function) , sizing_field_(c3t3.triangulation()) , time_limit_(-1) , running_time_() , do_freeze_(do_freeze) #ifdef CGAL_MESH_3_OPTIMIZER_VERBOSE , sum_moves_(0) #endif // CGAL_MESH_3_OPTIMIZER_VERBOSE { nb_frozen_points_ = 0; // We put it here in case it's an "atomic" // If we're multi-thread tr_.set_lock_data_structure(get_lock_data_structure()); #ifdef CGAL_MESH_3_OPTIMIZER_VERBOSE std::cerr << "Fill sizing field..."; CGAL::Real_timer timer; timer.start(); #endif fill_sizing_field(); #ifdef CGAL_MESH_3_OPTIMIZER_VERBOSE std::cerr << "done (" << timer.time() << "s)\n"; #endif } template Mesh_optimization_return_code Mesh_global_optimizer:: operator()(int nb_iterations, Visitor visitor) { running_time_.reset(); running_time_.start(); #ifdef CGAL_MESH_3_PROFILING WallClockTimer t; #endif // Fill set containing moving vertices // first, take them all #ifdef CGAL_CONSTRUCT_INTRUSIVE_LIST_RANGE_CONSTRUCTOR typedef CGAL::Prevent_deref It; Moving_vertices_set moving_vertices(It(tr_.finite_vertices_begin()), It(tr_.finite_vertices_end())); #else Moving_vertices_set moving_vertices; collect_all_vertices(moving_vertices); #endif std::size_t initial_vertices_nb = moving_vertices.size(); #ifdef CGAL_MESH_3_OPTIMIZER_VERBOSE double step_begin = running_time_.time(); std::cerr << "Running " << Mf::name() << "-smoothing (" << initial_vertices_nb << " vertices)" << std::endl; #endif //CGAL_MESH_3_OPTIMIZER_VERBOSE // Initialize big moves (stores the largest moves) this->clear_big_moves(); big_moves_size_ = (std::max)(std::size_t(1), std::size_t(moving_vertices.size()/500)); std::size_t nb_vertices_moved = (std::numeric_limits::max)(); bool convergence_stop = false; // Iterate int i = -1; while ( ++i < nb_iterations && ! is_time_limit_reached() ) { if(!do_freeze_) nb_frozen_points_ = 0; else nb_vertices_moved = moving_vertices.size(); // Compute move for each vertex Moves_vector moves = compute_moves(moving_vertices); visitor.after_compute_moves(); // Pb with Freeze : sometimes a few vertices continue moving indefinitely // if the nb of moving vertices is < 1% of total nb AND does not decrease if(do_freeze_ && double(nb_vertices_moved) < 0.005 * double(initial_vertices_nb) && nb_vertices_moved == moving_vertices.size()) { // we should stop because we are probably entering an infinite instable loop convergence_stop = true; break; } // Stop if time_limit is reached if ( is_time_limit_reached() ) break; // Update mesh with those moves update_mesh(moves, moving_vertices, visitor); visitor.end_of_iteration(i); #ifdef CGAL_MESH_3_OPTIMIZER_VERBOSE std::size_t moving_vertices_size = moving_vertices.size(); std::cerr << boost::format("\r \r" "end iter.%1%, %2% / %3%, last step:%4$.2fs, step avg:%5$.2fs, avg big moves:%6$.3f ") % (i+1) % moving_vertices_size % initial_vertices_nb % (running_time_.time() - step_begin) % (running_time_.time() / (i+1)) % sum_moves_; step_begin = running_time_.time(); #endif if (do_freeze_ && nb_frozen_points_ == initial_vertices_nb ) break; if(check_convergence()) break; } #ifdef CGAL_MESH_3_OPTIMIZER_VERBOSE std::cerr << std::endl; #endif running_time_.stop(); #ifdef CGAL_MESH_3_PROFILING double optim_time = t.elapsed(); # ifdef CGAL_MESH_3_EXPORT_PERFORMANCE_DATA CGAL_MESH_3_SET_PERFORMANCE_DATA(std::string(Mf::name()) + "_optim_time", optim_time); # endif #endif #ifdef CGAL_MESH_3_OPTIMIZER_VERBOSE if ( do_freeze_ && nb_frozen_points_ == initial_vertices_nb ) std::cerr << "All vertices frozen" << std::endl; else if ( do_freeze_ && convergence_stop ) std::cerr << "Can't improve anymore" << std::endl; else if ( is_time_limit_reached() ) std::cerr << "Time limit reached" << std::endl; else if ( check_convergence() ) std::cerr << "Convergence reached" << std::endl; else if( i >= nb_iterations ) std::cerr << "Max iteration number reached" << std::endl; std::cerr << "Total optimization time: " << running_time_.time() << "s" << std::endl << std::endl; #endif #ifdef CGAL_MESH_3_PROFILING std::cerr << std::endl << "Total optimization 'wall-clock' time: " << optim_time << "s" << std::endl; #endif if ( do_freeze_ && nb_frozen_points_ == initial_vertices_nb ) return ALL_VERTICES_FROZEN; else if ( do_freeze_ && convergence_stop ) return CANT_IMPROVE_ANYMORE; else if ( is_time_limit_reached() ) return TIME_LIMIT_REACHED; else if ( check_convergence() ) return CONVERGENCE_REACHED; return MAX_ITERATION_NUMBER_REACHED; } template void Mesh_global_optimizer:: collect_all_vertices(Moving_vertices_set& moving_vertices) { typename Tr::Finite_vertices_iterator vit = tr_.finite_vertices_begin(); for(; vit != tr_.finite_vertices_end(); ++vit ) moving_vertices.insert(vit); } template typename Mesh_global_optimizer::Moves_vector Mesh_global_optimizer:: compute_moves(Moving_vertices_set& moving_vertices) { // Store new position of points which have to move Moves_vector moves; moves.reserve(moving_vertices.size()); // reset worst_move list this->clear_big_moves(); #ifdef CGAL_MESH_3_PROFILING std::cerr << "Computing moves..."; WallClockTimer t; #endif #ifdef CGAL_LINKED_WITH_TBB // Parallel if (boost::is_convertible::value) { tbb::concurrent_vector vertices_not_moving_any_more; // Get move for each moving vertex tbb::parallel_do( moving_vertices.begin(), moving_vertices.end(), Compute_move( *this, sizing_field_, moves, do_freeze_, vertices_not_moving_any_more, tr_.geom_traits()) ); typename tbb::concurrent_vector::const_iterator it = vertices_not_moving_any_more.begin(); typename tbb::concurrent_vector::const_iterator it_end = vertices_not_moving_any_more.end(); for ( ; it != it_end ; ++it) { moving_vertices.erase(*it); } } // Sequential else #endif // CGAL_LINKED_WITH_TBB { typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object(); typename Gt::Construct_translated_point_3 translate = tr_.geom_traits().construct_translated_point_3_object(); // Get move for each moving vertex typename Moving_vertices_set::iterator vit = moving_vertices.begin(); for ( ; vit != moving_vertices.end() ; ) { Vertex_handle oldv = *vit; ++vit; Vector_3 move = compute_move(oldv); if ( CGAL::NULL_VECTOR != move ) { FT size = 0.; if(Sizing_field::is_vertex_update_needed) { const Weighted_point& position = tr_.point(oldv); Bare_point new_position = translate(cp(position), move); size = sizing_field_(new_position, oldv); } moves.push_back(std::make_tuple(oldv, move, size)); } else // CGAL::NULL_VECTOR == move { if(do_freeze_) { // TODO: if non-intrusive, we can optimize since we have the iterator. // Don't forget to do "vit = mv.erase(vit)" instead ++vit. moving_vertices.erase(oldv); } } // Stop if time_limit_ is reached if ( is_time_limit_reached() ) break; } } #ifdef CGAL_MESH_3_PROFILING std::cerr << "done in " << t.elapsed() << " seconds." << std::endl; #endif return moves; } template typename Mesh_global_optimizer::Vector_3 Mesh_global_optimizer:: compute_move(const Vertex_handle& v) { typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object(); typename Gt::Compute_squared_length_3 sq_length = tr_.geom_traits().compute_squared_length_3_object(); typename Gt::Construct_translated_point_3 translate = tr_.geom_traits().construct_translated_point_3_object(); typename Gt::Construct_vector_3 vector = tr_.geom_traits().construct_vector_3_object(); Cell_vector incident_cells; incident_cells.reserve(64); #ifdef CGAL_LINKED_WITH_TBB // Parallel if (boost::is_convertible::value) { tr_.incident_cells_threadsafe(v, std::back_inserter(incident_cells)); } else #endif //CGAL_LINKED_WITH_TBB { tr_.incident_cells(v, std::back_inserter(incident_cells)); } // Get move from move function Vector_3 move = move_function_(v, incident_cells, c3t3_, sizing_field_); FT local_sq_size = min_circumradius_sq_length(v, incident_cells); if ( FT(0) == local_sq_size ) return CGAL::NULL_VECTOR; // Project surface vertex if ( c3t3_.in_dimension(v) == 2 ) { const Weighted_point& position = tr_.point(v); Bare_point new_position = translate(cp(position), move); Bare_point projected_new_position = helper_.project_on_surface(v, new_position); move = vector(cp(position), projected_new_position); } FT local_move_sq_ratio = sq_length(move) / local_sq_size; // Move point only if the displacement is big enough w.r.t. the local size if ( local_move_sq_ratio < sq_freeze_ratio_ ) { ++nb_frozen_points_; return CGAL::NULL_VECTOR; } // Update big moves this->update_big_moves(local_move_sq_ratio); return move; } template void Mesh_global_optimizer:: update_mesh(const Moves_vector& moves, Moving_vertices_set& moving_vertices, Visitor& visitor) { // Cells which have to be updated Outdated_cell_set outdated_cells; #ifdef CGAL_MESH_3_PROFILING std::cerr << "Moving vertices..."; WallClockTimer t; #endif #ifdef CGAL_LINKED_WITH_TBB // Parallel if (boost::is_convertible::value) { // Apply moves in triangulation tbb::parallel_for(tbb::blocked_range(0, moves.size()), Move_vertex< Self, C3T3_helpers, Tr, Moves_vector, Moving_vertices_set, Outdated_cell_set>( *this, helper_, moves, moving_vertices, outdated_cells) , tbb_task_group_context ); } // Sequential else #endif // CGAL_LINKED_WITH_TBB { // Apply moves in triangulation for ( typename Moves_vector::const_iterator it = moves.begin() ; it != moves.end() ; ++it ) { const Vertex_handle& v = std::get<0>(*it); const Vector_3& move = std::get<1>(*it); // Get size at new position if ( Sizing_field::is_vertex_update_needed ) { FT size = std::get<2>(*it); #ifdef CGAL_MESH_3_OPTIMIZER_VERBOSE std::cout << "Moving #" << it - moves.begin() << " addr: " << &*v << " pt: " << tr_.point(v) << " move: " << move << std::endl; #endif // Move point Vertex_handle new_v = helper_.move_point(v, move, outdated_cells, moving_vertices); // Restore size in meshing_info data new_v->set_meshing_info(size); } else // Move point { helper_.move_point(v, move, outdated_cells, moving_vertices); } // Stop if time_limit_ is reached, here we can't return without rebuilding // restricted delaunay if ( is_time_limit_reached() ) break; } } visitor.after_move_points(); #ifdef CGAL_MESH_3_PROFILING std::cerr << "done in " << t.elapsed() << " seconds." << std::endl; #endif #ifdef CGAL_MESH_3_PROFILING std::cerr << "Updating C3T3 (rebuilding restricted Delaunay)..."; t.reset(); #endif // Update c3t3 #ifdef CGAL_INTRUSIVE_LIST // AF: rebuild_restricted_delaunay does more cell insertion/removal // which clashes with our inplace list // That's why we hand it in, and call clear() when we no longer need it helper_.rebuild_restricted_delaunay(outdated_cells, moving_vertices); #else helper_.rebuild_restricted_delaunay(outdated_cells.begin(), outdated_cells.end(), moving_vertices); #endif visitor.after_rebuild_restricted_delaunay(); #ifdef CGAL_MESH_3_PROFILING std::cerr << "Updating C3T3 done in " << t.elapsed() << " seconds." << std::endl; #endif } template void Mesh_global_optimizer:: fill_sizing_field() { std::map value_map; #ifdef CGAL_LINKED_WITH_TBB // Parallel if (boost::is_convertible::value) { typedef tbb::enumerable_thread_specific< std::vector< std::pair > > Local_list; Local_list local_lists; tbb::parallel_do( tr_.finite_vertices_begin(), tr_.finite_vertices_end(), Compute_sizing_field_value(*this, tr_.geom_traits(), local_lists) ); for(typename Local_list::iterator it_list = local_lists.begin() ; it_list != local_lists.end() ; ++it_list ) { value_map.insert(it_list->begin(), it_list->end()); } } else #endif //CGAL_LINKED_WITH_TBB { typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object(); // Fill map with local size for(typename Tr::Finite_vertices_iterator vit = tr_.finite_vertices_begin(); vit != tr_.finite_vertices_end(); ++vit) { const Weighted_point& position = tr_.point(vit); value_map.insert(std::make_pair(cp(position), average_circumradius_length(vit))); } } // Fill the sizing field sizing_field_.fill(value_map); } template bool Mesh_global_optimizer:: check_convergence() const { FT sum(0); for( typename std::multiset::const_iterator it = big_moves_.begin(), end = big_moves_.end() ; it != end ; ++it ) { sum += CGAL::sqrt(*it); } // Even if set is not full, divide by the max size so that if only 1 point moves, // then it goes to 0. FT average_move = sum / FT(big_moves_size_); #ifdef CGAL_MESH_3_OPTIMIZER_VERBOSE sum_moves_ = average_move; #endif return ( average_move < convergence_ratio_ ); } template typename Mesh_global_optimizer::FT Mesh_global_optimizer:: average_circumradius_length(const Vertex_handle& v) const { Cell_vector incident_cells; incident_cells.reserve(64); #ifdef CGAL_LINKED_WITH_TBB // Parallel if (boost::is_convertible::value) { tr_.incident_cells_threadsafe(v, std::back_inserter(incident_cells)); } else #endif //CGAL_LINKED_WITH_TBB { tr_.incident_cells(v, std::back_inserter(incident_cells)); } FT sum_len (0); unsigned int nb = 0; for ( typename Cell_vector::iterator cit = incident_cells.begin() ; cit != incident_cells.end() ; ++cit) { if ( c3t3_.is_in_complex(*cit) ) { sum_len += CGAL::sqrt(sq_circumradius_length(*cit,v)); ++nb; } } // nb == 0 could happen if there is an isolated point. if ( 0 != nb ) { return sum_len/nb; } else { // Use outside cells to compute size of point for ( typename Cell_vector::iterator cit = incident_cells.begin() ; cit != incident_cells.end() ; ++cit) { if ( !tr_.is_infinite(*cit) ) { sum_len += CGAL::sqrt(sq_circumradius_length(*cit,v)); ++nb; } } CGAL_assertion(nb!=0); CGAL_assertion(sum_len!=0); return sum_len/nb; } } template typename Mesh_global_optimizer::FT Mesh_global_optimizer:: min_circumradius_sq_length(const Vertex_handle& v, const Cell_vector& incident_cells) const { // Get first cell sq_circumradius_length typename Cell_vector::const_iterator cit = incident_cells.begin(); while ( incident_cells.end() != cit && !c3t3_.is_in_complex(*cit) ) { ++cit; } // if vertex is isolated ... if ( incident_cells.end() == cit ) return FT(0); // Initialize min FT min_sq_len = sq_circumradius_length(*cit++,v); // Find the minimum value for ( ; cit != incident_cells.end() ; ++cit ) { if ( !c3t3_.is_in_complex(*cit) ) continue; min_sq_len = (std::min)(min_sq_len,sq_circumradius_length(*cit,v)); } return min_sq_len; } template typename Mesh_global_optimizer::FT Mesh_global_optimizer:: sq_circumradius_length(const Cell_handle& cell, const Vertex_handle& v) const { typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object(); typename Gt::Compute_squared_distance_3 sq_distance = tr_.geom_traits().compute_squared_distance_3_object(); const Bare_point circumcenter = tr_.dual(cell); const Weighted_point& position = tr_.point(cell, cell->index(v)); return ( sq_distance(cp(position), circumcenter) ); } } // end namespace Mesh_3 } //namespace CGAL #include #endif // CGAL_MESH_3_MESH_GLOBAL_OPTIMIZER_H