1 // Copyright (c) 2009 INRIA Sophia-Antipolis (France).
2 // All rights reserved.
3 //
4 // This file is part of CGAL (www.cgal.org).
5 //
6 // $URL: https://github.com/CGAL/cgal/blob/v5.3/Mesh_3/include/CGAL/Mesh_3/C3T3_helpers.h $
7 // $Id: C3T3_helpers.h bdd4efe 2021-01-15T10:06:56+01:00 Sébastien Loriot
8 // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
9 //
10 //
11 // Author(s) : Stephane Tayeb
12 //
13 //******************************************************************************
14 //
15 //******************************************************************************
16
17 #ifndef CGAL_MESH_3_C3T3_HELPERS_H
18 #define CGAL_MESH_3_C3T3_HELPERS_H
19
20 #include <CGAL/license/Mesh_3.h>
21
22 #include <CGAL/disable_warnings.h>
23
24 #include <CGAL/Mesh_3/config.h>
25 #include <CGAL/use.h>
26
27 #include <CGAL/linear_least_squares_fitting_3.h>
28 #include <CGAL/Mesh_3/Triangulation_helpers.h>
29 #include <CGAL/Time_stamper.h>
30 #include <CGAL/tuple.h>
31 #include <CGAL/iterator.h>
32 #include <CGAL/array.h>
33 #include <CGAL/Handle_hash_function.h>
34
35 #ifdef CGAL_MESH_3_PROFILING
36 #include <CGAL/Mesh_3/Profiling_tools.h>
37 #endif
38
39 #include <boost/range/begin.hpp>
40 #include <boost/range/end.hpp>
41 #include <boost/optional.hpp>
42 #include <CGAL/boost/iterator/transform_iterator.hpp>
43 #include <boost/iterator/function_output_iterator.hpp>
44 #include <boost/type_traits/is_same.hpp>
45 #include <boost/type_traits/is_convertible.hpp>
46 #include <boost/unordered_set.hpp>
47
48 #ifdef CGAL_LINKED_WITH_TBB
49 # include <tbb/parallel_do.h>
50 # include <mutex>
51 #endif
52
53 #include <functional>
54 #include <vector>
55 #include <set>
56
57 namespace CGAL {
58 namespace Mesh_3 {
59
60 #ifdef CGAL_INTRUSIVE_LIST
61 template <typename Type>
62 class Intrusive_list {
63 public:
64
65 typedef Type Type_handle;
66 typedef Type_handle& reference;
67 typedef const Type_handle& const_reference;
68 typedef Type_handle value_type;
69
Intrusive_list()70 Intrusive_list()
71 : f(), b(), n(0)
72 {}
73
~Intrusive_list()74 ~Intrusive_list()
75 {
76 clear();
77 }
78
79
Intrusive_list(const Intrusive_list &)80 Intrusive_list(const Intrusive_list& )
81 {
82 CGAL_assertion(false);
83 }
84
85 #ifdef CGAL_CONSTRUCT_INTRUSIVE_LIST_RANGE_CONSTRUCTOR
86 template <typename IT>
Intrusive_list(IT first,IT last)87 Intrusive_list(IT first, IT last)
88 : f(), b(), n(0)
89 {
90 if(first == last){
91 return;
92 }
93
94 f = *first;
95 Type_handle ch = f;
96 ++n;
97 ++first;
98 while(first != last){
99 if((ch != Type(*first)) && ((*first)->next_intrusive()==Type_handle())){
100 // not yet inserted
101 ch->set_next_intrusive(*first);
102 (*first)->set_previous_intrusive(ch);
103 ch = *first;
104 ++n;
105 }
106 ++first;
107 }
108 b = ch;
109 b->set_next_intrusive(f);
110 f->set_previous_intrusive(b);
111 }
112 #endif
113
114 bool
is_valid()115 is_valid() const
116 {
117 if(n < 0){
118 std::cerr << "n < 0" << std::endl;
119 return false;
120 }
121 if(n == 0){
122 if (f != Type_handle()){
123 std::cerr << "n==0, but f!= Type_handle()" << std::endl;
124 return false;
125 }
126 if (b != Type_handle()){
127 std::cerr << "n==0, but b!= Type_handle()" << std::endl;
128 return false;
129 }
130 }else{
131 if(f->previous_intrusive() != b){
132 std::cerr << "f->previous_intrusive() != b" << std::endl;
133 return false;
134 }
135 if(b->next_intrusive() != f){
136 std::cerr << "b->next_intrusive() != f" << std::endl;
137 return false;
138 }
139
140
141 Type_handle ch = f;
142 for(std::size_t i = 1; i < n; i++){
143 if(ch->next_intrusive()->previous_intrusive() != ch){
144 std::cerr << "ch->next_intrusive()->previous_intrusive() != ch" << std::endl;
145 return false;
146 }
147 ch = ch->next_intrusive();
148 }
149 if(ch != b){
150 std::cerr << "ch!= b)" << std::endl;
151 return false;
152 }
153 }
154 return true;
155 }
156
157
clear()158 void clear()
159 {
160 if(!empty()){
161 while( f!= b ){
162 Type_handle h = f;
163 f=f->next_intrusive();
164 h->set_previous_intrusive(Type_handle());
165 h->set_next_intrusive(Type_handle());
166 }
167 b->set_previous_intrusive(Type_handle());
168 b->set_next_intrusive(Type_handle());
169 f = b = Type_handle();
170 }
171 n = 0;
172 }
173
size()174 std::size_t size() const
175 {
176 return n;
177 }
178
179
180 struct iterator {
181 Type_handle pos, b;
182
183 typedef Type_handle value_type;
184 typedef const Type_handle* pointer;
185 typedef const Type_handle& reference;
186 typedef std::size_t size_type;
187 typedef std::ptrdiff_t difference_type;
188 typedef std::forward_iterator_tag iterator_category;
189
iteratoriterator190 iterator(Type_handle f, Type_handle b)
191 : pos(f), b(b)
192 {}
193
iteratoriterator194 iterator()
195 : pos()
196 {}
197
198 iterator& operator++()
199 {
200 if(pos != Type_handle()){
201 if(pos == b){
202 pos = Type_handle(); // past the end
203 }else {
204 pos = pos->next_intrusive();
205 }
206 }
207 return *this;
208 }
209
210 iterator operator++(int)
211 {
212 iterator tmp(*this);
213 ++(*this);
214 return tmp;
215 }
216
217 bool operator==(const iterator& i) const
218 {
219 return pos == i.pos;
220 }
221
222 bool operator!=(const iterator& i) const
223 {
224 return !(*this == i);
225 }
226
227 reference operator*() const
228 {
229 return pos;
230 }
231
232 pointer operator->() const
233 {
234 return pos;
235 }
236 }; // struct iterator
237
238
begin()239 iterator begin()
240 {
241 return iterator(f,b);
242 }
243
end()244 iterator end()
245 {
246 return iterator();
247 }
248
249
front()250 Type_handle front() const
251 {
252 return f;
253 }
254
front()255 Type_handle& front()
256 {
257 return f;
258 }
259
260
back()261 Type_handle back() const
262 {
263 return b;
264 }
265
back()266 Type_handle& back()
267 {
268 return b;
269 }
270
insert(iterator,const Type_handle & ch)271 iterator insert(iterator /* position */,
272 const Type_handle& ch)
273 {
274 CGAL_assertion( (ch->next_intrusive() == Type_handle() && ch->previous_intrusive() == Type_handle()) ||
275 (ch->next_intrusive() != Type_handle() && ch->previous_intrusive() != Type_handle()) );
276 CGAL_expensive_assertion(is_valid());
277
278 if(ch->next_intrusive() != Type_handle()){
279 return iterator(ch->next_intrusive()/*first*/, ch/*last*/);
280 }
281 else{
282 insert(ch);
283 return iterator(ch->next_intrusive()/*first*/, ch/*last*/);
284 }
285 }
286
insert(Type_handle ch)287 void insert(Type_handle ch)
288 {
289 CGAL_assertion( (ch->next_intrusive() == Type_handle() && ch->previous_intrusive() == Type_handle()) ||
290 (ch->next_intrusive() != Type_handle() && ch->previous_intrusive() != Type_handle()) );
291 CGAL_expensive_assertion(is_valid());
292
293 if(ch->next_intrusive() != Type_handle()){
294 return;
295 }
296 if(empty()){
297 f = b = ch;
298 ch->set_next_intrusive(ch);
299 ch->set_previous_intrusive(ch);
300 } else {
301 ch->set_next_intrusive(f);
302 ch->set_previous_intrusive(b);
303 f->set_previous_intrusive(ch);
304 b->set_next_intrusive(ch);
305 b = ch;
306 }
307 n++;
308 }
309
erase(Type_handle ch)310 void erase(Type_handle ch)
311 {
312 CGAL_assertion( (ch->next_intrusive() == Type_handle() && ch->previous_intrusive() == Type_handle()) ||
313 (ch->next_intrusive() != Type_handle() && ch->previous_intrusive() != Type_handle()) );
314 CGAL_expensive_assertion(is_valid());
315 if(ch->next_intrusive() == Type_handle()){
316 return;
317 }
318 if(f == b){ // only 1 element in the list
319 CGAL_assertion(f == ch);
320 CGAL_assertion(n == 1);
321
322 f = b = Type_handle();
323 } else {
324 if(f == ch){
325 f = f->next_intrusive();
326 }
327 if(b == ch){
328 b = b->previous_intrusive();
329 }
330 Type_handle p = ch->previous_intrusive(), n = ch->next_intrusive();
331 p->set_next_intrusive(n);
332 n->set_previous_intrusive(p);
333 }
334 ch->set_next_intrusive(Type_handle());
335 ch->set_previous_intrusive(Type_handle());
336 CGAL_assertion(ch->next_intrusive() == Type_handle());
337 CGAL_assertion(ch->previous_intrusive() == Type_handle());
338 n--;
339 }
340
empty()341 bool empty() const
342 {
343 if(f == Type_handle()){
344 CGAL_assertion(b == Type_handle());
345 CGAL_assertion(n == 0);
346 }
347 return f == Type_handle();
348 }
349
contains(Type_handle th)350 bool contains(Type_handle th) const
351 {
352 if(th->next_intrusive() == Type_handle())
353 {
354 CGAL_assertion(th->previous_intrusive() == Type_handle());
355 return true;
356 }
357 else return false;
358 }
359
push_back(Type_handle ch)360 void push_back(Type_handle ch)
361 {
362 insert(ch);
363 }
364
365 private:
366 Type_handle f,b;
367 std::size_t n;
368 };
369 #endif // #ifdef CGAL_INTRUSIVE_LIST
370
371
372 /************************************************
373 // Class C3T3_helpers_base
374 // Two versions: sequential / parallel
375 ************************************************/
376
377 // Sequential
378 template <typename Tr, typename Concurrency_tag>
379 class C3T3_helpers_base
380 {
381 protected:
382 typedef typename Tr::Geom_traits Gt;
383 typedef typename Tr::Bare_point Bare_point;
384 typedef typename Tr::Weighted_point Weighted_point;
385 typedef typename Gt::FT FT;
386 typedef typename Tr::Vertex_handle Vertex_handle;
387 typedef typename Tr::Cell_handle Cell_handle;
388 typedef typename Tr::Facet Facet;
389 typedef typename Tr::Lock_data_structure Lock_data_structure;
390
C3T3_helpers_base(Lock_data_structure *)391 C3T3_helpers_base(Lock_data_structure *) {}
392
get_lock_data_structure()393 Lock_data_structure *get_lock_data_structure() const
394 {
395 return 0;
396 }
397
398 public:
399 // Dummy locks/unlocks
400 bool try_lock_point(const Weighted_point &, int = 0) const
401 {
402 return true;
403 }
404
405 bool try_lock_vertex(Vertex_handle, int = 0) const
406 {
407 return true;
408 }
409
410 bool try_lock_point_no_spin(const Weighted_point &, int = 0) const
411 {
412 return true;
413 }
414
415 bool try_lock_vertex_no_spin(Vertex_handle, int = 0) const
416 {
417 return true;
418 }
419
420 bool try_lock_element(Cell_handle, int = 0) const
421 {
422 return true;
423 }
424
425 bool try_lock_element(const Facet &, int = 0) const
426 {
427 return true;
428 }
429
430
is_point_locked_by_this_thread(const Weighted_point &)431 bool is_point_locked_by_this_thread(const Weighted_point &) const
432 { return false; }
433
is_cell_locked_by_this_thread(const Cell_handle &)434 bool is_cell_locked_by_this_thread(const Cell_handle &) const
435 { return false; }
436
unlock_all_elements()437 void unlock_all_elements() const {}
438
439 // Dummy locks/unlocks
lock_outdated_cells()440 void lock_outdated_cells() const {}
unlock_outdated_cells()441 void unlock_outdated_cells() const {}
lock_moving_vertices()442 void lock_moving_vertices() const {}
unlock_moving_vertices()443 void unlock_moving_vertices() const {}
lock_vertex_to_proj()444 void lock_vertex_to_proj() const {}
unlock_vertex_to_proj()445 void unlock_vertex_to_proj() const {}
446 };
447
448 #ifdef CGAL_LINKED_WITH_TBB
449 // Parallel
450 template <typename Tr>
451 class C3T3_helpers_base<Tr, Parallel_tag>
452 {
453 protected:
454 typedef typename Tr::Geom_traits Gt;
455 typedef typename Tr::Bare_point Bare_point;
456 typedef typename Tr::Weighted_point Weighted_point;
457 typedef typename Tr::Vertex_handle Vertex_handle;
458 typedef typename Tr::Cell_handle Cell_handle;
459 typedef typename Tr::Facet Facet;
460 typedef typename Tr::Lock_data_structure Lock_data_structure;
461
C3T3_helpers_base(Lock_data_structure * lock_ds)462 C3T3_helpers_base(Lock_data_structure *lock_ds)
463 : m_lock_ds(lock_ds) {}
464
465
466 public:
467 // LOCKS (CONCURRENCY)
468
469 /*Lock_data_structure *get_lock_data_structure() const
470 {
471 return m_lock_ds;
472 }*/
473
474 bool try_lock_point(const Weighted_point &p, int lock_radius = 0) const
475 {
476 if (m_lock_ds)
477 {
478 return m_lock_ds->try_lock(p, lock_radius);
479 }
480 return true;
481 }
482
483 bool try_lock_vertex(Vertex_handle vh, int lock_radius = 0) const
484 {
485 if (m_lock_ds)
486 {
487 return m_lock_ds->try_lock(vh->point(), lock_radius);
488 }
489 return true;
490 }
491
492 bool try_lock_point_no_spin(const Weighted_point &p, int lock_radius = 0) const
493 {
494 if (m_lock_ds)
495 {
496 return m_lock_ds->template try_lock<true>(p, lock_radius);
497 }
498 return true;
499 }
500
501 bool try_lock_vertex_no_spin(Vertex_handle vh, int lock_radius = 0) const
502 {
503 return try_lock_point_no_spin(vh->point(), lock_radius);
504 }
505
506 bool try_lock_element(Cell_handle cell_handle, int lock_radius = 0) const
507 {
508 bool success = true;
509
510 // Lock the element area on the grid
511 for (int iVertex = 0 ; success && iVertex < 4 ; ++iVertex)
512 {
513 Vertex_handle vh = cell_handle->vertex(iVertex);
514 success = try_lock_vertex(vh, lock_radius);
515 }
516
517 return success;
518 }
519
520 bool try_lock_element(const Facet &facet, int lock_radius = 0) const
521 {
522 bool success = true;
523
524 // Lock the element area on the grid
525 Cell_handle cell = facet.first;
526 for (int iVertex = (facet.second+1)&3 ;
527 success && iVertex != facet.second ; iVertex = (iVertex+1)&3)
528 {
529 Vertex_handle vh = cell->vertex(iVertex);
530 success = try_lock_vertex(vh, lock_radius);
531 }
532
533 return success;
534 }
535
is_point_locked_by_this_thread(const Weighted_point & p)536 bool is_point_locked_by_this_thread(const Weighted_point &p) const
537 {
538 bool locked = true;
539 if (m_lock_ds)
540 {
541 locked = m_lock_ds->is_locked_by_this_thread(p);
542 }
543 return locked;
544 }
545
is_cell_locked_by_this_thread(const Cell_handle & cell_handle)546 bool is_cell_locked_by_this_thread(const Cell_handle &cell_handle) const
547 {
548 bool locked = true;
549 if (m_lock_ds)
550 {
551 for (int iVertex = 0 ; locked && iVertex < 4 ; ++iVertex)
552 {
553 locked = m_lock_ds->is_locked_by_this_thread(
554 cell_handle->vertex(iVertex)->point());
555 }
556 }
557 return locked;
558 }
559
unlock_all_elements()560 void unlock_all_elements() const
561 {
562 if (m_lock_ds)
563 {
564 m_lock_ds->unlock_all_points_locked_by_this_thread();
565 }
566 }
567
lock_outdated_cells()568 void lock_outdated_cells() const
569 {
570 m_mut_outdated_cells.lock();
571 }
unlock_outdated_cells()572 void unlock_outdated_cells() const
573 {
574 m_mut_outdated_cells.unlock();
575 }
576
lock_moving_vertices()577 void lock_moving_vertices() const
578 {
579 m_mut_moving_vertices.lock();
580 }
unlock_moving_vertices()581 void unlock_moving_vertices() const
582 {
583 m_mut_moving_vertices.unlock();
584 }
585
lock_vertex_to_proj()586 void lock_vertex_to_proj() const
587 {
588 m_mut_vertex_to_proj.lock();
589 }
unlock_vertex_to_proj()590 void unlock_vertex_to_proj() const
591 {
592 m_mut_vertex_to_proj.unlock();
593 }
594
595 protected:
596 Lock_data_structure *m_lock_ds;
597
598 typedef std::mutex Mutex_type;
599 mutable Mutex_type m_mut_outdated_cells;
600 mutable Mutex_type m_mut_moving_vertices;
601 mutable Mutex_type m_mut_vertex_to_proj;
602 };
603 #endif // CGAL_LINKED_WITH_TBB
604
605
606 /************************************************
607 *
608 * C3T3_helpers class
609 *
610 ************************************************/
611
612 template <typename C3T3,
613 typename MeshDomain>
614 class C3T3_helpers
615 : public C3T3_helpers_base<typename C3T3::Triangulation,
616 typename C3T3::Concurrency_tag>
617 {
618 // -----------------------------------
619 // Private types
620 // -----------------------------------
621 typedef C3T3_helpers<C3T3, MeshDomain> Self;
622 typedef C3T3_helpers_base<typename C3T3::Triangulation,
623 typename C3T3::Concurrency_tag> Base;
624 typedef typename C3T3::Concurrency_tag Concurrency_tag;
625
626 typedef typename Base::Lock_data_structure Lock_data_structure;
627 typedef typename C3T3::Triangulation Tr;
628 typedef Tr Triangulation;
629 typedef typename Tr::Geom_traits Gt;
630
631 typedef typename Gt::FT FT;
632 typedef typename Tr::Bare_point Bare_point;
633 typedef typename Tr::Weighted_point Weighted_point;
634 typedef typename Gt::Vector_3 Vector_3;
635 typedef typename Gt::Plane_3 Plane_3;
636 typedef typename Gt::Tetrahedron_3 Tetrahedron;
637
638 typedef typename Tr::Vertex_handle Vertex_handle;
639 typedef typename Tr::Facet Facet;
640 typedef typename Tr::Cell Cell;
641 typedef typename Tr::Cell_handle Cell_handle;
642
643 typedef typename C3T3::Surface_patch_index Surface_patch_index;
644 typedef typename C3T3::Subdomain_index Subdomain_index;
645 typedef typename C3T3::Index Index;
646
647 typedef boost::optional<Surface_patch_index> Surface_patch;
648 typedef boost::optional<Subdomain_index> Subdomain;
649
650 typedef std::vector<Cell_handle> Cell_vector;
651 typedef std::set<Cell_handle> Cell_set;
652 typedef std::vector<Tetrahedron> Tet_vector;
653
654 typedef std::vector<Facet> Facet_vector;
655 typedef std::vector<Vertex_handle> Vertex_vector;
656
657 typedef CGAL::Hash_handles_with_or_without_timestamps Hash_fct;
658 typedef boost::unordered_set<Vertex_handle, Hash_fct> Vertex_set;
659
660 #ifdef CGAL_INTRUSIVE_LIST
661 typedef Intrusive_list<Cell_handle> Outdated_cell_set;
662 #else
663 typedef Cell_set Outdated_cell_set;
664 #endif //CGAL_INTRUSIVE_LIST
665
666 #ifdef CGAL_INTRUSIVE_LIST
667 typedef Intrusive_list<Vertex_handle> Moving_vertices_set;
668 #else
669 typedef Vertex_set Moving_vertices_set;
670 #endif //CGAL_INTRUSIVE_LIST
671
672 private:
673 typedef std::pair<Vertex_handle, Vertex_handle> Ordered_edge;
674 // Vertex_data is the dimension and Index of the third vertex of the facet.
675 typedef std::pair<int, Index> Vertex_data;
676 typedef std::pair<Surface_patch_index, Vertex_data> Facet_topology_description;
677 // Facet_boundary stores the edges, of the boundary of surface facets, with meta-data.
678 typedef std::map<Ordered_edge,
679 Facet_topology_description> Facet_boundary;
680
681 typedef Triangulation_helpers<Tr> Th;
682
683 public:
684 // -----------------------------------
685 // Public interface
686 // -----------------------------------
687 typedef boost::optional<Vertex_handle> Update_mesh;
688
689 using Base::try_lock_point;
690 using Base::try_lock_vertex;
691 using Base::try_lock_point_no_spin;
692 using Base::try_lock_vertex_no_spin;
693 using Base::try_lock_element;
694 using Base::is_point_locked_by_this_thread;
695 using Base::is_cell_locked_by_this_thread;
696 using Base::unlock_all_elements;
697 using Base::lock_outdated_cells;
698 using Base::unlock_outdated_cells;
699 using Base::lock_moving_vertices;
700 using Base::unlock_moving_vertices;
701 using Base::lock_vertex_to_proj;
702 using Base::unlock_vertex_to_proj;
703
704 /**
705 * Constructor
706 */
707 C3T3_helpers(C3T3& c3t3, const MeshDomain& domain,
708 Lock_data_structure *lock_ds = nullptr)
Base(lock_ds)709 : Base(lock_ds)
710 , c3t3_(c3t3)
711 , tr_(c3t3.triangulation())
712 , domain_(domain)
713 { }
714
715 /**
716 * @brief tries to move \c old_vertex to \c new_position in the mesh
717 * @param old_vertex the old vertex
718 * @param move the translation from the old position to the new
719 * @param new_position the new position of \c old_vertex
720 * @param criterion the criterion which will be used to verify the new
721 * position is ok. c3t3 minimal value of new criterion shall not decrease.
722 * @param modified_vertices contains the vertices incident to cells which
723 * may have been impacted by relocation
724 * @return a pair which contains:
725 * - a bool which is \c true if the move has been done.
726 * - a Vertex_handle which is always filled and may be the new vertex (if
727 * the move is a success), or the vertex which lies at \c v's position in
728 * the updated c3t3.
729 */
730 template <typename SliverCriterion, typename OutputIterator>
731 std::pair<bool,Vertex_handle>
732 update_mesh(const Vertex_handle& old_vertex,
733 const Vector_3& move,
734 const Weighted_point& new_position,
735 const SliverCriterion& criterion,
736 OutputIterator modified_vertices,
737 bool *could_lock_zone = nullptr);
738
739 /** @brief tries to move \c old_vertex to \c new_position in the mesh
740 *
741 * Same as update_mesh, but with the precondition that
742 * Th().no_topological_change(tr_, old_vertex, new_position,
743 * incident_cells_) return false.
744 */
745 template <typename SliverCriterion, typename OutputIterator>
746 std::pair<bool,Vertex_handle>
747 update_mesh_topo_change(const Vertex_handle& old_vertex,
748 const Weighted_point& new_position,
749 const SliverCriterion& criterion,
750 OutputIterator modified_vertices,
751 bool *could_lock_zone = nullptr);
752
753 /**
754 * Updates mesh moving vertex \c old_vertex to \c new_position. Returns the
755 * new vertex of the triangulation.
756 *
757 * Insert into modified vertices the vertices which are impacted by to move.
758 */
759 template <typename OutputIterator>
760 Vertex_handle update_mesh(const Vertex_handle& old_vertex,
761 const Vector_3& new_position,
762 OutputIterator modified_vertices,
763 bool fill_modified_vertices = true);
764
765 /**
766 * Updates mesh moving vertex \c old_vertex to \c new_position. Returns the
767 * new vertex of the triangulation.
768 */
update_mesh(const Vertex_handle & old_vertex,const Vector_3 & move)769 Vertex_handle update_mesh(const Vertex_handle& old_vertex,
770 const Vector_3& move)
771 {
772 return update_mesh(old_vertex, move, Emptyset_iterator(), false);
773 }
774
775 /**
776 * Rebuilds restricted Delaunay
777 */
778 template <typename ForwardIterator>
rebuild_restricted_delaunay(ForwardIterator first_cell,ForwardIterator last_cell)779 void rebuild_restricted_delaunay(ForwardIterator first_cell,
780 ForwardIterator last_cell)
781 {
782 Moving_vertices_set mvs;
783 return rebuild_restricted_delaunay(first_cell, last_cell, mvs);
784 }
785
786 /**
787 * Rebuilds restricted Delaunay
788 */
789 template <typename ForwardIterator>
790 void rebuild_restricted_delaunay(ForwardIterator first_cell,
791 ForwardIterator last_cell,
792 Moving_vertices_set& moving_vertices);
793
794 void update_restricted_facets();
795
796 #ifdef CGAL_INTRUSIVE_LIST
797 template <typename OutdatedCells>
798 void rebuild_restricted_delaunay(OutdatedCells& outdated_cells,
799 Moving_vertices_set& moving_vertices);
800 #endif
801
802 /**
803 * @brief Project \c p on surface, using incident facets of \c v
804 * @param v The vertex from which p was moved
805 * @param p The point to project
806 * @param index The index of the surface patch where v lies, if known.
807 * @return the projected point
808 *
809 * \c p is projected as follows using normal of least square fitting plane
810 * on \c v incident surface points. If \c index is specified, only
811 * surface points that are on the same surface patch are used to compute
812 * the fitting plane.
813 */
814 Bare_point
815 project_on_surface(const Vertex_handle& v, const Bare_point& p,
816 Surface_patch_index index = Surface_patch_index()) const;
817
818 /**
819 * Returns the minimum value for criterion for incident cells of \c vh
820 */
821 template <typename SliverCriterion>
822 FT min_incident_value(const Vertex_handle& vh,
823 const SliverCriterion& criterion) const;
824
825 /**
826 * Moves \c old_vertex to \c new_position
827 * Stores the cells which have to be updated in \c outdated_cells
828 * Updates the Vertex_handle old_vertex to its new value in \c moving_vertices
829 * The second one (with the could_lock_zone param) is for the parallel version
830 */
831 Vertex_handle move_point(const Vertex_handle& old_vertex,
832 const Vector_3& move,
833 Outdated_cell_set& outdated_cells_set,
834 Moving_vertices_set& moving_vertices) const;
835
836 Vertex_handle move_point(const Vertex_handle& old_vertex,
837 const Vector_3& move,
838 Outdated_cell_set& outdated_cells_set,
839 Moving_vertices_set& moving_vertices,
840 bool *could_lock_zone) const;
841
842 /**
843 * Try to lock ALL the incident cells and return in \c cells the ones
844 * whose \c filter says "true".
845 * Return value:
846 * - false: everything is unlocked and \c cells is empty
847 * - true: ALL incident cells are locked and \c cells is filled
848 */
849 template <typename Filter>
850 bool
851 try_lock_and_get_incident_cells(const Vertex_handle& v,
852 Cell_vector &cells,
853 const Filter &filter) const;
854
855 /**
856 * Try to lock ALL the incident cells and return in \c cells the slivers
857 * Return value:
858 * - false: everything is unlocked and \c cells is empty
859 * - true: incident cells are locked and \c cells contains all slivers
860 */
861 template <typename SliverCriterion>
862 bool
863 try_lock_and_get_incident_slivers(const Vertex_handle& v,
864 const SliverCriterion& criterion,
865 const FT& sliver_bound,
866 Cell_vector &cells) const;
867
868 template <typename SliverCriterion>
869 void
870 get_incident_slivers_without_using_tds_data(const Vertex_handle& v,
871 const SliverCriterion& criterion,
872 const FT& sliver_bound,
873 Cell_vector &slivers) const;
874
875 /**
876 * Outputs to out the sliver (wrt \c criterion and \c sliver_bound) incident
877 * to \c v
878 */
879 template <typename SliverCriterion, typename OutputIterator>
880 OutputIterator
881 incident_slivers(const Vertex_handle& v,
882 const SliverCriterion& criterion,
883 const FT& sliver_bound,
884 OutputIterator out) const;
885
886
887 template <typename SliverCriterion, typename OutputIterator>
888 OutputIterator
889 new_incident_slivers(const Vertex_handle& v,
890 const SliverCriterion& criterion,
891 const FT& sliver_bound,
892 OutputIterator out) const;
893
894 /**
895 * Returns the sliver (wrt \c criterion and \c sliver_bound) incident to \c v
896 */
897 template <typename SliverCriterion>
898 Cell_vector
incident_slivers(const Vertex_handle & v,const SliverCriterion & criterion,const FT & sliver_bound)899 incident_slivers(const Vertex_handle& v,
900 const SliverCriterion& criterion,
901 const FT& sliver_bound) const
902 {
903 Cell_vector slivers;
904 incident_slivers(v, criterion, sliver_bound, std::back_inserter(slivers));
905 return slivers;
906 }
907
908 template <typename SliverCriterion>
909 Cell_vector
new_incident_slivers(const Vertex_handle & v,const SliverCriterion & criterion,const FT & sliver_bound)910 new_incident_slivers(const Vertex_handle& v,
911 const SliverCriterion& criterion,
912 const FT& sliver_bound) const
913 {
914 Cell_vector slivers;
915 new_incident_slivers(v, criterion, sliver_bound, std::back_inserter(slivers));
916 return slivers;
917 }
918
919
920 /**
921 * Returns the number of slivers incident to \c v
922 */
923 template <typename SliverCriterion>
924 std::size_t
925 number_of_incident_slivers(const Vertex_handle& v,
926 const SliverCriterion& criterion,
927 const FT& sliver_bound) const;
928
929 template <typename SliverCriterion>
930 bool
931 is_sliver(const Cell_handle& ch,
932 const SliverCriterion& criterion,
933 const FT& sliver_bound) const;
934
935 /**
936 * Returns the minimum criterion value of cells contained in \c cells
937 * Precondition: cells of \c cells must not be infinite.
938 * Warning: Here we don't check if cells are in c3t3
939 */
940 template <typename SliverCriterion>
941 FT min_sliver_value(const Cell_vector& cells,
942 const SliverCriterion& criterion,
943 const bool use_cache = true) const;
944
945 /**
946 * Reset cache validity of all cells of c3t3_
947 */
reset_cache()948 void reset_cache() const
949 {
950 for(typename C3T3::Cells_in_complex_iterator it = c3t3_.cells_in_complex_begin();
951 it != c3t3_.cells_in_complex_end(); ++it)
952 it->reset_cache_validity();
953 }
954
955 private:
956 // -----------------------------------
957 // Useful Functors
958 // -----------------------------------
959 /**
960 * @class Get_all_facets
961 *
962 * A functor which adds to an output iterator canonical facets of a cell
963 */
964 template <typename OutputIterator>
965 class Get_all_facets
966 {
967 public:
Get_all_facets(const Triangulation & tr,OutputIterator out)968 Get_all_facets(const Triangulation& tr, OutputIterator out)
969 : tr_(tr)
970 , out_(out) {}
971
operator()972 void operator()(const Cell_handle& cell)
973 {
974 #ifndef CGAL_MESH_3_NEW_GET_FACETS
975 for ( int i=0 ; i<4 ; ++i )
976 if ( !tr_.is_infinite(cell,i) )
977 *out_++ = canonical_facet(cell,i);
978 #else
979 // Instead of iterating over the facets we iterate over the vertices
980 // If a vertex is infinite we report only the facet opposite to it and return
981 // If all vertices are finite we report all facets
982 // This approach makes less tests if vertices are infinite
983 int i=0;
984 for ( ; i<4 ; ++i ){
985 if ( tr_.is_infinite(cell->vertex(i)) ){
986 *out_++ = canonical_facet(cell,i);
987 return;
988 }
989 }
990 for ( i=0 ; i<4 ; ++i ){
991 *out_++ = canonical_facet(cell,i);
992 }
993 #endif
994 }
995
996 private:
canonical_facet(const Cell_handle & c,const int i)997 Facet canonical_facet(const Cell_handle& c, const int i) const
998 {
999 #ifndef CGAL_MESH_3_NEW_GET_FACETS
1000 Facet facet(c,i);
1001 Facet mirror = tr_.mirror_facet(facet);
1002 return ( (mirror<facet)?mirror:facet );
1003 #else
1004 Cell_handle n = c->neighbor(i);
1005 if(c < n){
1006 return Facet(c,i);
1007 }else{
1008 return Facet(n,n->index(c));
1009 }
1010 #endif
1011 }
1012
1013 private:
1014 const Triangulation& tr_;
1015 OutputIterator out_;
1016 };
1017
1018
1019 /**
1020 * @class Is_in_c3t3
1021 *
1022 * A functor which returns true if a given handle is in c3t3
1023 */
1024 template <typename Handle>
1025 class Is_in_c3t3 : public CGAL::cpp98::unary_function<Handle, bool>
1026 {
1027 public:
Is_in_c3t3(const C3T3 & c3t3)1028 Is_in_c3t3(const C3T3& c3t3) : c3t3_(c3t3) { }
operator()1029 bool operator()(const Handle& h) const { return c3t3_.is_in_complex(h); }
1030
1031 private:
1032 const C3T3& c3t3_;
1033 };
1034
1035
1036 /**
1037 * @class Is_sliver
1038 *
1039 * A functor which answers true if a Cell_handle is a sliver
1040 */
1041 template <typename SliverCriterion>
1042 struct Is_sliver : public CGAL::cpp98::unary_function<Cell_handle,bool>
1043 {
1044 Is_sliver(const C3T3& c3t3,
1045 const SliverCriterion& criterion,
1046 const FT& bound = 0)
c3t3_Is_sliver1047 : c3t3_(c3t3)
1048 , criterion_(criterion)
1049 , bound_(bound) { }
1050
operatorIs_sliver1051 bool operator()(const Cell_handle& c) const
1052 {
1053 if ( c3t3_.is_in_complex(c) )
1054 {
1055 CGAL_assertion(!c3t3_.triangulation().is_infinite(c));
1056
1057 if ( ! c->is_cache_valid() )
1058 {
1059 Sliver_criterion_value<SliverCriterion> sc_value(c3t3_.triangulation(),
1060 criterion_);
1061 (void) sc_value(c); // 'sc_value::operator()' updates the cache of 'c'
1062 }
1063 else
1064 {
1065 CGAL_expensive_assertion(c->sliver_value() == criterion_(c));
1066 }
1067 if(bound_ > 0)
1068 return ( c->sliver_value() <= bound_ );
1069 else
1070 return ( c->sliver_value() <= criterion_.sliver_bound() );
1071 }
1072 else
1073 return false;
1074 }
1075
1076 private:
1077 const C3T3& c3t3_;
1078 const SliverCriterion& criterion_;
1079 const FT bound_;
1080 };
1081
1082
1083 /**
1084 * @class Update_c3t3
1085 *
1086 * A functor which updates c3t3 w.r.t. the domain.
1087 */
1088 class Update_c3t3
1089 {
1090 public:
Update_c3t3(const MeshDomain & domain,C3T3 & c3t3)1091 Update_c3t3(const MeshDomain& domain, C3T3& c3t3)
1092 : domain_(domain)
1093 , c3t3_(c3t3) {}
1094
1095 /**
1096 * @brief Updates facet \c facet in c3t3
1097 * @param facet the facet to update
1098 * @param update if set to \c false, checking only is done
1099 * @return true if \c facet is in c3t3
1100 */
operator()1101 Surface_patch operator()(const Facet& facet, const bool update = true) const
1102 {
1103 return this->operator()(facet, update, update);
1104 }
1105
1106 /**
1107 * @brief Updates facet \c facet in c3t3
1108 * @param facet the facet to update
1109 * @param update_c3t3 if set to \c false, checking only is done
1110 * @param update_surface_center if set to \c true, the facet surface
1111 * center is updated.
1112 * @return true if \c facet is in c3t3
1113 *
1114 * By default, \c update_c3t3 is \c true, and \c update_surface_center
1115 * is equal to \c update_c3t3.
1116 */
operator()1117 Surface_patch operator()(const Facet& facet,
1118 const bool update_c3t3,
1119 const bool update_surface_center) const
1120 {
1121 typedef typename C3T3::Triangulation::Geom_traits Gt;
1122 typedef typename Gt::Segment_3 Segment_3;
1123 typedef typename Gt::Ray_3 Ray_3;
1124 typedef typename Gt::Line_3 Line_3;
1125
1126 // Nothing to do for infinite facets
1127 if ( c3t3_.triangulation().is_infinite(facet) )
1128 return Surface_patch();
1129
1130 // Functors
1131 typename Gt::Is_degenerate_3 is_degenerate =
1132 c3t3_.triangulation().geom_traits().is_degenerate_3_object();
1133
1134 // Get dual of facet
1135 Object dual = c3t3_.triangulation().dual(facet);
1136
1137 // The dual is a segment, a ray or a line
1138 if ( const Segment_3* p_segment = object_cast<Segment_3>(&dual) )
1139 {
1140 if (is_degenerate(*p_segment))
1141 return Surface_patch();
1142
1143 return dual_intersect(*p_segment,facet,
1144 update_c3t3,
1145 update_surface_center);
1146 }
1147 else if ( const Ray_3* p_ray = object_cast<Ray_3>(&dual) )
1148 {
1149 if (is_degenerate(*p_ray))
1150 return Surface_patch();
1151
1152 return dual_intersect(*p_ray,facet,update_c3t3,
1153 update_surface_center);
1154 }
1155 else if ( const Line_3* p_line = object_cast<Line_3>(&dual) )
1156 {
1157 return dual_intersect(*p_line,facet,update_c3t3,
1158 update_surface_center);
1159 }
1160
1161 CGAL_error_msg("This should not happen");
1162 return Surface_patch();
1163 }
1164
1165 /**
1166 * @brief Updates cell \c ch in c3t3
1167 * @param ch the cell to update
1168 * @param update if set to \c false, checking only is done
1169 * @return true if \c ch is in c3t3
1170 */
operator()1171 Subdomain operator()(const Cell_handle& ch, const bool update = true) const
1172 {
1173 if ( c3t3_.triangulation().is_infinite(ch) )
1174 return Subdomain();
1175
1176 // treat cell
1177 const Subdomain subdomain =
1178 domain_.is_in_domain_object()(c3t3_.triangulation().dual(ch));
1179 // function dual(cell) updates the circumcenter cache if there is one
1180
1181 if ( subdomain && update )
1182 {
1183 c3t3_.add_to_complex(ch,*subdomain);
1184 }
1185 else if(update)
1186 {
1187 c3t3_.remove_from_complex(ch);
1188 }
1189
1190 return subdomain;
1191 }
1192
1193 private:
1194
1195 // Returns true if query intersects the surface.
1196 template <typename Query>
dual_intersect(const Query & dual,const Facet & facet,const bool update_c3t3,const bool update_surface_center)1197 Surface_patch dual_intersect(const Query& dual,
1198 const Facet& facet,
1199 const bool update_c3t3,
1200 const bool update_surface_center) const
1201 {
1202 typedef typename MeshDomain::Intersection Intersection;
1203
1204 typename MeshDomain::Construct_intersection construct_intersection =
1205 domain_.construct_intersection_object();
1206
1207 #ifndef CGAL_MESH_3_NO_LONGER_CALLS_DO_INTERSECT_3
1208
1209 typename MeshDomain::Do_intersect_surface do_intersect_surface =
1210 domain_.do_intersect_surface_object();
1211 Surface_patch surface = do_intersect_surface(dual);
1212
1213 #else // CGAL_MESH_3_NO_LONGER_CALLS_DO_INTERSECT_3
1214
1215 Intersection intersection = construct_intersection(dual);
1216 Surface_patch surface =
1217 (std::get<2>(intersection) == 0) ? Surface_patch() :
1218 Surface_patch(
1219 domain_.surface_patch_index(std::get<1>(intersection)));
1220
1221 #endif // CGAL_MESH_3_NO_LONGER_CALLS_DO_INTERSECT_3
1222
1223 // Update if needed
1224 if(update_c3t3)
1225 {
1226 // Update status in c3t3
1227 if((bool)surface)
1228 c3t3_.add_to_complex(facet, *surface);
1229 else
1230 c3t3_.remove_from_complex(facet);
1231 }
1232
1233 if(update_surface_center)
1234 {
1235 Facet facet_m = c3t3_.triangulation().mirror_facet(facet);
1236 if(surface)
1237 {
1238 #ifndef CGAL_MESH_3_NO_LONGER_CALLS_DO_INTERSECT_3
1239 Intersection intersection = construct_intersection(dual);
1240 #endif // NOT CGAL_MESH_3_NO_LONGER_CALLS_DO_INTERSECT_3
1241
1242 // Update facet surface center
1243 Bare_point surface_center = std::get<0>(intersection);
1244 facet.first->set_facet_surface_center(facet.second, surface_center);
1245 facet_m.first->set_facet_surface_center(facet_m.second, surface_center);
1246 }
1247 }
1248
1249 return surface ? surface : Surface_patch();
1250 }
1251
1252
1253 private:
1254 const MeshDomain& domain_;
1255 C3T3& c3t3_;
1256 }; //end class Update_c3t3
1257
1258 class Facet_updater
1259 {
1260 const Self& m_c3t3_helpers;
1261 Vertex_set& vertex_to_proj;
1262 C3T3& c3t3_;
1263 Update_c3t3& c3t3_updater_;
1264
1265 public:
1266 typedef Facet& reference;
1267 typedef const Facet& const_reference;
1268
Facet_updater(const Self & c3t3_helpers,C3T3 & c3t3,Vertex_set & vertex_to_proj,Update_c3t3 & c3t3_updater_)1269 Facet_updater(const Self& c3t3_helpers,
1270 C3T3& c3t3, Vertex_set& vertex_to_proj, Update_c3t3& c3t3_updater_)
1271 : m_c3t3_helpers(c3t3_helpers),
1272 vertex_to_proj(vertex_to_proj), c3t3_(c3t3), c3t3_updater_(c3t3_updater_)
1273 {}
1274
1275 void
operator()1276 operator()(const Facet& f)
1277 {
1278 // Update facet
1279 c3t3_.remove_from_complex(f);
1280 c3t3_updater_(f);
1281
1282 // Update vertex_to_proj
1283 if ( c3t3_.is_in_complex(f) )
1284 {
1285 // Iterate on vertices
1286 int k = f.second;
1287 for ( int i=1 ; i<4 ; ++i )
1288 {
1289 const Vertex_handle& v = f.first->vertex((k+i)&3);
1290 if ( c3t3_.in_dimension(v) > 2 )
1291 {
1292 m_c3t3_helpers.lock_vertex_to_proj();
1293 vertex_to_proj.insert(v);
1294 m_c3t3_helpers.unlock_vertex_to_proj();
1295 }
1296 }
1297 }
1298 }
1299
1300 }; // end class Facet_updater
1301
1302
1303 /**
1304 * @class Sliver_criterion_value
1305 *
1306 * A functor which returns sliver criterion value for a Cell_handle
1307 * and updates its cache value
1308 */
1309 template <typename SliverCriterion>
1310 class Sliver_criterion_value
1311 : public CGAL::cpp98::unary_function<Cell_handle, double>
1312 {
1313 public:
Sliver_criterion_value(const Tr & tr,const SliverCriterion & criterion)1314 Sliver_criterion_value(const Tr& tr,
1315 const SliverCriterion& criterion)
1316 : p_tr_(&tr)
1317 , criterion_(criterion) {}
1318
operator()1319 FT operator()(const Cell_handle& ch) const
1320 {
1321 CGAL_precondition(!p_tr_->is_infinite(ch));
1322
1323 if ( ! ch->is_cache_valid() )
1324 {
1325 double sliver_value = criterion_(ch);
1326 ch->set_sliver_value(sliver_value);
1327 }
1328 else
1329 {
1330 CGAL_expensive_assertion(ch->sliver_value() == criterion_(ch));
1331 }
1332 return ch->sliver_value();
1333 }
1334
1335 private:
1336 // '=' is used, so p_tr_ must be a pointer ...
1337 const Tr* p_tr_;
1338 SliverCriterion criterion_;
1339 };
1340
1341 /**
1342 * to be used by the perturber
1343 */
1344 class Cell_from_ids
1345 {
1346 public:
Cell_from_ids(const C3T3 & c3t3,const Cell_handle & c)1347 Cell_from_ids(const C3T3& c3t3, const Cell_handle& c)
1348 : infinite_(c3t3.triangulation().is_infinite(c))
1349 , vertices_()
1350 , sorted_vertices_()
1351 {
1352 for(int i = 0; i < 4; ++i)
1353 {
1354 if (c3t3.triangulation().is_infinite(c->vertex(i)))
1355 continue;
1356 //the Id is set with an int by Sliver_perturber,
1357 // in initialize_vertices_id
1358 int id = static_cast<int>(c->vertex(i)->meshing_info());
1359 vertices_.push_back(id);
1360 }
1361 sorted_vertices_ = vertices_;//makes a copy of each element
1362 std::sort(sorted_vertices_.begin(), sorted_vertices_.end());
1363 CGAL_assertion((infinite_ && vertices_.size() == 3)
1364 || vertices_.size() == 4);
1365 }
1366
vertex_id(const std::size_t & i)1367 std::size_t vertex_id(const std::size_t& i) const
1368 {
1369 CGAL_precondition((infinite_ && i < 3) || i < 4);
1370 return vertices_[i];
1371 }
1372
1373 bool operator<(const Cell_from_ids& c) const
1374 {
1375 //std::array operator< compares lhs and rhs lexicographically
1376 return sorted_vertices_ < c.sorted_vertices_;
1377 }
1378
1379 private:
1380 bool infinite_;
1381 // vertices IDs, not sorted, to keep the ordering of the Cell_handle id's
1382 std::vector<int> vertices_;
1383 // vertices IDs, sorted, to be found in a std::set<Cell_from_ids>
1384 std::vector<int> sorted_vertices_;
1385 };
1386
1387 class Cell_data_backup
1388 {
1389 public:
1390 Cell_data_backup(const C3T3& c3t3,
1391 const Cell_handle& c,
1392 const bool do_backup = true)
cell_ids_(c3t3,c)1393 : cell_ids_(c3t3, c)
1394 {
1395 //backup is not done when constructor is called to
1396 //convert a newly created cell (has nothing to backup)
1397 //to a Cell_data_backup
1398 if(do_backup)
1399 {
1400 if (!c3t3.triangulation().is_infinite(c))
1401 backup_finite_cell(c);
1402 else
1403 backup_infinite_cell(c, c3t3);
1404 }
1405 }
1406
1407 private:
backup_finite_cell(const Cell_handle & c)1408 void backup_finite_cell(const Cell_handle& c)
1409 {
1410 if(c->is_cache_valid())
1411 sliver_value_ = c->sliver_value();
1412 else
1413 sliver_value_ = 0.;
1414
1415 subdomain_index_ = c->subdomain_index();
1416 for(std::size_t i = 0; i < 4; ++i)
1417 {
1418 const int ii = static_cast<int>(i);//avoid warnings
1419 surface_index_table_[i] = c->surface_patch_index(ii);
1420 facet_surface_center_[i] = c->get_facet_surface_center(ii);
1421 surface_center_index_table_[i] = c->get_facet_surface_center_index(ii);
1422 }
1423 //note c->next_intrusive() and c->previous_intrusive()
1424 //are lost by 'backup' and 'restore',
1425 //because all cells are changing during the move
1426 //they are not used in update_mesh functions involving a Sliver_criterion
1427 }
1428
backup_infinite_cell(const Cell_handle & c,const C3T3 & c3t3)1429 void backup_infinite_cell(const Cell_handle& c,
1430 const C3T3& c3t3)
1431 {
1432 for (int ii = 0; ii < 4; ++ii)
1433 {
1434 if (c3t3.triangulation().is_infinite(c->vertex(ii)))
1435 {
1436 surface_index_table_[0] = c->surface_patch_index(ii);
1437 facet_surface_center_[0] = c->get_facet_surface_center(ii);
1438 surface_center_index_table_[0] = c->get_facet_surface_center_index(ii);
1439 break;
1440 }
1441 }
1442 }
1443
1444 public:
1445 bool operator<(const Cell_data_backup& cb) const
1446 {
1447 return cell_ids_ < cb.cell_ids_;
1448 }
1449
1450 /**
1451 * new_cell has the same vertices as cell_ids_
1452 * (checked before function is called)
1453 * resets new_cell's meta-data to its back-uped values
1454 */
restore(Cell_handle new_cell,C3T3 & c3t3)1455 void restore(Cell_handle new_cell, C3T3& c3t3)
1456 {
1457 if (c3t3.triangulation().is_infinite(new_cell))
1458 return restore_infinite_cell(new_cell, c3t3);
1459
1460 IndexMap new_to_old_indices;
1461 CGAL_assertion_code(unsigned int nbv_found = 0);
1462 for(int i = 0; i < 4; ++i)
1463 {
1464 std::size_t new_vi_index =
1465 static_cast<std::size_t>(new_cell->vertex(i)->meshing_info());
1466 for(std::size_t j = 0; j < 4; ++j)
1467 {
1468 if(new_vi_index == cell_ids_.vertex_id(j))
1469 {
1470 new_to_old_indices[static_cast<std::size_t>(i)] = j;
1471 CGAL_assertion_code(++nbv_found);
1472 break;//loop on j
1473 }
1474 }//end loop j
1475 }//end loop i
1476 CGAL_assertion(nbv_found == 4);
1477
1478 restore(new_cell, new_to_old_indices, c3t3);
1479 }
1480
1481 private:
1482 typedef std::array<std::size_t, 4> IndexMap;
1483
restore(Cell_handle c,const IndexMap & index_map,C3T3 & c3t3)1484 void restore(Cell_handle c,
1485 const IndexMap& index_map,//new_to_old_indices
1486 C3T3& c3t3)
1487 {
1488 if(sliver_value_ > 0.)
1489 c->set_sliver_value(sliver_value_);
1490
1491 for(int i = 0; i < 4; ++i)
1492 c->reset_visited(i);
1493 //we don't need to store 'visited' information because it is
1494 //reset and used locally where it is needed
1495
1496 //add_to_complex sets the index, and updates the cell counter
1497 //if c should be in the c3t3, add_to_complex has to be used
1498 //to increment the nb of cells and facets in c3t3
1499 if(!( Subdomain_index() == subdomain_index_ ))
1500 c3t3.add_to_complex(c, subdomain_index_);
1501 else
1502 c3t3.remove_from_complex(c);
1503
1504 for(int i = 0; i < 4; ++i)
1505 {
1506 std::size_t old_i = index_map.at(static_cast<std::size_t>(i));
1507 Surface_patch_index index = surface_index_table_[old_i];
1508 //add_to_complex sets the index, and updates the facet counter
1509 if(!( Surface_patch_index() == index ))
1510 c3t3.add_to_complex(Facet(c, i), index);
1511 else
1512 c3t3.remove_from_complex(Facet(c,i));
1513
1514 c->set_facet_surface_center(i, facet_surface_center_[old_i]);
1515 const Facet mirror = c3t3.triangulation().mirror_facet(Facet(c, i));
1516 mirror.first->set_facet_surface_center(mirror.second, facet_surface_center_[old_i]);
1517 }
1518 }
1519
restore_infinite_cell(Cell_handle c,C3T3 & c3t3)1520 void restore_infinite_cell(Cell_handle c,
1521 C3T3& c3t3)
1522 {
1523 c3t3.remove_from_complex(c);//infinite
1524 for (unsigned int i = 0; i < 4; ++i)
1525 {
1526 if (!c3t3.triangulation().is_infinite(Facet(c,i)))
1527 {
1528 Surface_patch_index index = surface_index_table_[0];
1529 if (!( Surface_patch_index() == index ))
1530 c3t3.add_to_complex(Facet(c, i), index);
1531 else
1532 c3t3.remove_from_complex(Facet(c, i));
1533
1534 c->set_facet_surface_center(i, facet_surface_center_[0]);
1535 const Facet mirror = c3t3.triangulation().mirror_facet(Facet(c, i));
1536 mirror.first->set_facet_surface_center(mirror.second, facet_surface_center_[0]);
1537 return;
1538 }
1539 }
1540 }
1541
1542 private:
1543 typedef typename Tr::Cell::Subdomain_index Subdomain_index;
1544 typedef typename Tr::Cell::Surface_patch_index Surface_patch_index;
1545 typedef typename Tr::Cell::Index Index;
1546
1547 Cell_from_ids cell_ids_;
1548 FT sliver_value_;
1549 Subdomain_index subdomain_index_;
1550 std::array<Surface_patch_index, 4> surface_index_table_;
1551 std::array<Bare_point, 4> facet_surface_center_;
1552 std::array<Index, 4> surface_center_index_table_;
1553 };
1554
1555 private:
1556 // -----------------------------------
1557 // Private methods
1558 // -----------------------------------
1559 /**
1560 * Returns the minimum criterion value of c3t3 cells contained in \c cells.
1561 */
1562 template <typename SliverCriterion>
1563 FT min_sliver_in_c3t3_value(const Cell_vector& cells,
1564 const SliverCriterion& criterion,
1565 const bool use_cache = true) const
1566 {
1567 // Get complex cells only
1568 Cell_vector c3t3_cells_ = c3t3_cells(cells);
1569 return min_sliver_value(c3t3_cells_, criterion, use_cache);
1570 }
1571
c3t3_cells(const Cell_vector & cells)1572 Cell_vector c3t3_cells(const Cell_vector& cells) const
1573 {
1574 Cell_vector c3t3_cells;
1575 #ifdef CGAL_CXX17
1576 std::remove_copy_if(cells.begin(),
1577 cells.end(),
1578 std::back_inserter(c3t3_cells),
1579 std::not_fn(Is_in_c3t3<Cell_handle>(c3t3_)));
1580 #else
1581 std::remove_copy_if(cells.begin(),
1582 cells.end(),
1583 std::back_inserter(c3t3_cells),
1584 std::not1(Is_in_c3t3<Cell_handle>(c3t3_)) );
1585 #endif
1586 return c3t3_cells;
1587 }
1588
1589 /**
1590 * Removes objects of [begin,end[ range from \c c3t3_
1591 */
1592 template<typename ForwardIterator>
remove_from_c3t3(ForwardIterator begin,ForwardIterator end)1593 void remove_from_c3t3(ForwardIterator begin, ForwardIterator end) const
1594 {
1595 while ( begin != end )
1596 c3t3_.remove_from_complex(*begin++);
1597 }
1598
1599 /**
1600 * Remove cells and facets of \c cells from c3t3
1601 */
1602 template < typename ForwardIterator >
remove_cells_and_facets_from_c3t3(ForwardIterator cells_begin,ForwardIterator cells_end)1603 void remove_cells_and_facets_from_c3t3(ForwardIterator cells_begin,
1604 ForwardIterator cells_end) const
1605 {
1606 Facet_vector facets = get_facets_not_inplace(cells_begin,cells_end);
1607 remove_from_c3t3(facets.begin(), facets.end());
1608 remove_from_c3t3(cells_begin, cells_end);
1609 }
1610
1611 /**
1612 * Insert into \c out the vertices of range [cells_begin,cells_end[
1613 */
1614 template <typename InputIterator, typename OutputIterator>
1615 void fill_modified_vertices(InputIterator cells_begin,
1616 InputIterator cells_end,
1617 const Vertex_handle& vertex,
1618 OutputIterator out) const;
1619
1620 /**
1621 * Backup cells meta-data to a vector of Cell_data_backup
1622 */
1623 template <typename CellsVector, typename CellDataSet>
1624 void fill_cells_backup(const CellsVector& cells,
1625 CellDataSet& cells_backup) const;
1626
1627 template <typename CellsVector, typename CellDataSet>
1628 void restore_from_cells_backup(const CellsVector& cells,
1629 CellDataSet& cells_backup) const;
1630
1631
1632 /**
1633 * Update mesh iff sliver criterion value does not decrease.
1634 */
1635 template <typename SliverCriterion, typename OutputIterator>
1636 std::pair<bool,Vertex_handle>
1637 update_mesh_no_topo_change(const Vertex_handle& old_vertex,
1638 const Vector_3& move,
1639 const Weighted_point& new_position,
1640 const SliverCriterion& criterion,
1641 OutputIterator modified_vertices,
1642 const Cell_vector& conflict_cells);
1643
1644 /**
1645 * Move point and returns the set of cells that are not valid anymore, and
1646 * the set of cells which have been deleted by the move process.
1647 */
1648 template < typename OutdatedCellsOutputIterator,
1649 typename DeletedCellsOutputIterator >
1650 Vertex_handle move_point(const Vertex_handle& old_vertex,
1651 const Vector_3& move,
1652 OutdatedCellsOutputIterator outdated_cells,
1653 DeletedCellsOutputIterator deleted_cells) const;
1654
1655 Vertex_handle
1656 move_point_topo_change(const Vertex_handle& old_vertex,
1657 const Weighted_point& new_position,
1658 Outdated_cell_set& outdated_cells_set,
1659 bool *could_lock_zone = nullptr) const;
1660
1661 template < typename OutdatedCellsOutputIterator,
1662 typename DeletedCellsOutputIterator >
1663 Vertex_handle
1664 move_point_topo_change(const Vertex_handle& old_vertex,
1665 const Weighted_point& new_position,
1666 OutdatedCellsOutputIterator outdated_cells,
1667 DeletedCellsOutputIterator deleted_cells) const;
1668
1669 Vertex_handle move_point_topo_change(const Vertex_handle& old_vertex,
1670 const Weighted_point& new_position) const;
1671
1672 template < typename OutdatedCellsOutputIterator >
1673 Vertex_handle
1674 move_point_no_topo_change(const Vertex_handle& old_vertex,
1675 const Vector_3& move,
1676 const Weighted_point& new_position,
1677 OutdatedCellsOutputIterator outdated_cells) const;
1678
1679 Vertex_handle
1680 move_point_no_topo_change(const Vertex_handle& old_vertex,
1681 const Vector_3& move,
1682 const Weighted_point& new_position) const;
1683
1684 /**
1685 * Returns the least square plane from v, using adjacent surface points
1686 */
1687 boost::optional<Plane_3>
1688 get_least_square_surface_plane(const Vertex_handle& v,
1689 Bare_point& ref_point,
1690 Surface_patch_index index = Surface_patch_index()) const;
1691
1692 /**
1693 * @brief Project \c p on surface, using incident facets of \c v
1694 * @param v The vertex from which p was moved
1695 * @param p The point to project
1696 * @param index The index of the surface patch where v lies, if known.
1697 * @return a `boost::optional` with the projected point if the projection
1698 * was possible, or `boost::none`.
1699 *
1700 * \c p is projected using the normal of least square fitting plane
1701 * on \c v incident surface points. If \c index is specified, only
1702 * surface points that are on the same surface patch are used to compute
1703 * the fitting plane.
1704 */
1705 boost::optional<Bare_point>
1706 project_on_surface_if_possible(const Vertex_handle& v,
1707 const Bare_point& p,
1708 Surface_patch_index index = Surface_patch_index()) const;
1709
1710 /**
1711 * @brief Returns the projection of \c p, using direction of
1712 * \c projection_vector
1713 */
1714 Bare_point
1715 project_on_surface_aux(const Bare_point& p,
1716 const Bare_point& ref_point,
1717 const Vector_3& projection_vector) const;
1718
1719 /**
1720 * Reverts the move from \c old_point to \c new_vertex. Returns the inserted
1721 * vertex located at \c old_point
1722 * and an output iterator on outdated cells
1723 */
1724 template<typename OutputIterator>
revert_move(const Vertex_handle & new_vertex,const Weighted_point & old_point,OutputIterator outdated_cells)1725 Vertex_handle revert_move(const Vertex_handle& new_vertex,
1726 const Weighted_point& old_point,
1727 OutputIterator outdated_cells)
1728 {
1729 // Move vertex
1730 Vertex_handle revert_vertex =
1731 move_point_topo_change(new_vertex,
1732 old_point,
1733 outdated_cells,
1734 CGAL::Emptyset_iterator()); //deleted cells
1735 CGAL_assertion(Vertex_handle() != revert_vertex);
1736
1737 return revert_vertex;
1738 }
1739
1740 /**
1741 * Returns the boundary of restricted facets of \c facets,
1742 and the list of vertices of all restricted facets,
1743 which should not contain the vertex that is moving
1744 */
1745 Facet_boundary
1746 get_surface_boundary(const Vertex_handle& moving_vertex,
1747 const Facet_vector& facets,
1748 Vertex_set& incident_surface_vertices) const;
1749
1750 /**
1751 * Returns the boundary of restricted facets of \c cells
1752 and the list of vertices of all restricted facets.
1753 */
1754 Facet_boundary
get_surface_boundary(const Vertex_handle & moving_vertex,const Cell_vector & cells,Vertex_set & incident_surface_vertices)1755 get_surface_boundary(const Vertex_handle& moving_vertex,
1756 const Cell_vector& cells,
1757 Vertex_set& incident_surface_vertices) const
1758 {
1759 return get_surface_boundary(moving_vertex,
1760 get_facets(cells),
1761 incident_surface_vertices);
1762 }
1763
1764 /**
1765 * Returns false if there is a vertex belonging to one facet of \c facets
1766 * which has not his dimension < 3
1767 */
1768 bool check_no_inside_vertices(const Facet_vector& facets) const;
1769
1770 /**
1771 * Returns the impacted cells when moving \c vertex to \c conflict_point
1772 */
1773 template <typename OutputIterator>
1774 OutputIterator
1775 get_conflict_zone_no_topo_change(const Vertex_handle& vertex,
1776 OutputIterator conflict_cells) const;
1777
1778 template <typename OutputIterator>
1779 OutputIterator
1780 get_conflict_zone_topo_change(const Vertex_handle& vertex,
1781 const Weighted_point& conflict_point,
1782 OutputIterator conflict_cells) const;
1783
1784 template <typename CellsOutputIterator,
1785 typename FacetsOutputIterator>
1786 void
1787 get_conflict_zone_topo_change(const Vertex_handle& v,
1788 const Weighted_point& conflict_point,
1789 CellsOutputIterator insertion_conflict_cells,
1790 FacetsOutputIterator insertion_conflict_boundary,
1791 CellsOutputIterator removal_conflict_cells,
1792 bool *could_lock_zone = nullptr) const;
1793
1794
1795 template < typename ConflictCellsInputIterator,
1796 typename OutdatedCellsOutputIterator,
1797 typename DeletedCellsOutputIterator >
1798 Vertex_handle
1799 move_point_topo_change_conflict_zone_known(const Vertex_handle& old_vertex,
1800 const Weighted_point& new_position,
1801 const Facet& insertion_boundary_facet,
1802 ConflictCellsInputIterator insertion_conflict_cells_begin,
1803 ConflictCellsInputIterator insertion_conflict_cells_end,
1804 ConflictCellsInputIterator removal_conflict_cells_begin,
1805 ConflictCellsInputIterator removal_conflict_cells_end,
1806 OutdatedCellsOutputIterator outdated_cells,
1807 DeletedCellsOutputIterator deleted_cells) const;
1808
1809 /**
1810 * Updates \c boundary wrt \c edge: if edge is already in boundary we remove
1811 * it, else we add it.
1812 */
update_boundary(Facet_boundary & boundary,const Ordered_edge & edge,const Vertex_handle third_vertex,const Surface_patch_index & surface_index)1813 void update_boundary(Facet_boundary& boundary,
1814 const Ordered_edge& edge,
1815 const Vertex_handle third_vertex,
1816 const Surface_patch_index& surface_index) const
1817 {
1818 const typename Facet_boundary::value_type x =
1819 std::make_pair(edge,
1820 std::make_pair(surface_index,
1821 std::make_pair(c3t3_.in_dimension(third_vertex),
1822 c3t3_.index(third_vertex)
1823 )
1824 )
1825 );
1826 typename Facet_boundary::iterator boundary_it =
1827 boundary.find(edge);
1828
1829 if ( boundary_it != boundary.end() )
1830 boundary.erase(boundary_it);
1831 else
1832 boundary.insert(x);
1833 }
1834
1835 /**
1836 * Returns the facets of \c cells (returns each facet only once i.e. use
1837 * canonical facet)
1838 */
get_facets(const Cell_vector & cells)1839 Facet_vector get_facets(const Cell_vector& cells) const
1840 {
1841 return get_facets(cells.begin(),cells.end());
1842 }
1843
1844 // TODO: write get_facets so that it uses update_facets with a FacetUpdater that calls push_back
1845
1846 #if defined(CGAL_MESH_3_GET_FACETS_USING_INTRUSIVE_LIST) && defined(CGAL_INTRUSIVE_LIST)
1847 template <typename ForwardIterator>
get_facets(ForwardIterator first_cell,ForwardIterator last_cell)1848 Facet_vector get_facets(ForwardIterator first_cell,
1849 ForwardIterator last_cell) const
1850 {
1851 Facet_vector result; // AF: todo: resize?
1852 #ifdef CGAL_CONSTRUCT_INTRUSIVE_LIST_RANGE_CONSTRUCTOR
1853 Intrusive_list<Cell_handle> outdated_cells(first_cell, last_cell);
1854 #else
1855 Intrusive_list<Cell_handle> outdated_cells;
1856 for( ;first_cell!= last_cell; ++first_cell){
1857 outdated_cells.insert(*first_cell);
1858 }
1859 #endif
1860
1861 for(typename Intrusive_list<Cell_handle>::iterator it = outdated_cells.begin();
1862 it != outdated_cells.end();
1863 ++it){
1864 Cell_handle cell = *it;
1865 int i=0;
1866 bool inf = false;
1867 for ( ; i<4 && (!inf) ; ++i ){
1868 if ( tr_.is_infinite(cell->vertex(i)) ){
1869 inf = true;
1870 Cell_handle n = cell->neighbor(i);
1871 if(n->next_intrusive() != Cell_handle()){// the neighbor is also outdated
1872 if(cell < n){ // otherwise n will report it later
1873 result.push_back(Facet(cell,i));
1874 }
1875 } else { // report it now or never
1876 if(cell < n){
1877 result.push_back(Facet(cell,i));
1878 }else {
1879 result.push_back(Facet(n,n->index(cell)));
1880 }
1881 }
1882 }
1883 }
1884 if(! inf){
1885 for ( i=0 ; i<4 ; ++i ){
1886 Cell_handle n = cell->neighbor(i);
1887 if(n->next_intrusive() != Cell_handle()){// the neighbor is also outdated
1888 if(cell < n){ // otherwise n will report it later
1889 result.push_back(Facet(cell,i));
1890 }
1891 } else { // report it now or never
1892 if(cell < n){
1893 result.push_back(Facet(cell,i));
1894 }else {
1895 result.push_back(Facet(n,n->index(cell)));
1896 }
1897 }
1898 }
1899 }
1900 }
1901 return result;
1902 }
1903 #else
1904 /**
1905 * Returns the facets of \c cells (returns each facet only once i.e. use
1906 * canonical facet)
1907 */
1908 template <typename ForwardIterator>
get_facets(ForwardIterator first_cell,ForwardIterator last_cell)1909 Facet_vector get_facets(ForwardIterator first_cell,
1910 ForwardIterator last_cell) const
1911 {
1912 // Get all facets
1913 typedef Get_all_facets<std::back_insert_iterator<Facet_vector> > Get_facets;
1914
1915 Facet_vector all_facets;
1916 all_facets.reserve(64);
1917 std::for_each(first_cell,
1918 last_cell,
1919 Get_facets(tr_,std::back_inserter(all_facets)));
1920
1921 std::sort(all_facets.begin(), all_facets.end());
1922
1923 // Keep one copy of each facet (maybe copy could be avoided)
1924 // typename Facet_vector::iterator all_facets_end =
1925 // std::unique(all_facets.begin(), all_facets.end());
1926 Facet_vector facets;
1927 facets.reserve(64);
1928 std::unique_copy(all_facets.begin(),
1929 all_facets.end(),
1930 std::back_inserter(facets));
1931
1932 return facets;
1933 }
1934 #endif
1935
1936 #ifdef CGAL_INTRUSIVE_LIST
1937 template <typename FacetUpdater>
update_facets(Intrusive_list<Cell_handle> & outdated_cells,FacetUpdater updater)1938 void update_facets(Intrusive_list<Cell_handle>& outdated_cells, FacetUpdater updater)
1939 {
1940 # ifdef CGAL_LINKED_WITH_TBB
1941 // Parallel
1942 if (boost::is_convertible<Concurrency_tag, Parallel_tag>::value)
1943 {
1944 tbb::parallel_do(
1945 outdated_cells.begin(), outdated_cells.end(),
1946 Update_cell_facets<Self, FacetUpdater>(tr_, updater));
1947 }
1948 // Sequential
1949 else
1950 # endif // CGAL_LINKED_WITH_TBB
1951 {
1952 typename Intrusive_list<Cell_handle>::iterator it;
1953 for(it = outdated_cells.begin();
1954 it != outdated_cells.end();
1955 ++it)
1956 {
1957 Update_cell_facets<Self, FacetUpdater> ucf(tr_, updater);
1958 ucf(*it);
1959 }
1960 }
1961 }
1962 #endif //CGAL_INTRUSIVE_LIST
1963
1964 // Used by the parallel version
1965 template <typename FacetUpdater>
update_facets(std::vector<Cell_handle> & outdated_cells_vector,FacetUpdater updater)1966 void update_facets(std::vector<Cell_handle>& outdated_cells_vector, FacetUpdater updater)
1967 {
1968 # ifdef CGAL_LINKED_WITH_TBB
1969 // Parallel
1970 if (boost::is_convertible<Concurrency_tag, Parallel_tag>::value)
1971 {
1972 tbb::parallel_for
1973 (
1974 tbb::blocked_range<size_t>(0, outdated_cells_vector.size()),
1975 Update_cell_facets_for_parallel_for<Self, FacetUpdater>(
1976 tr_, updater, outdated_cells_vector)
1977 );
1978 }
1979 // Sequential
1980 else
1981 # endif // CGAL_LINKED_WITH_TBB
1982 {
1983 typename std::vector<Cell_handle>::iterator it;
1984 for(it = outdated_cells_vector.begin();
1985 it != outdated_cells_vector.end();
1986 ++it)
1987 {
1988 Update_cell_facets<Self, FacetUpdater> ucf(tr_, updater);
1989 ucf(*it);
1990 }
1991 }
1992 }
1993
1994
1995 /**
1996 * Returns the facets of \c cells (returns each facet only once i.e. use
1997 * canonical facet)
1998 */
1999 template <typename ForwardIterator>
get_facets_not_inplace(ForwardIterator first_cell,ForwardIterator last_cell)2000 Facet_vector get_facets_not_inplace(ForwardIterator first_cell,
2001 ForwardIterator last_cell) const
2002 {
2003 typedef Get_all_facets<std::back_insert_iterator<Facet_vector> > Get_facets;
2004
2005 Facet_vector all_facets;
2006 all_facets.reserve(64);
2007 std::for_each(first_cell,
2008 last_cell,
2009 Get_facets(tr_,std::back_inserter(all_facets)));
2010
2011 std::sort(all_facets.begin(), all_facets.end());
2012
2013 // Keep one copy of each facet (maybe copy could be avoided)
2014 // typename Facet_vector::iterator all_facets_end =
2015 // std::unique(all_facets.begin(), all_facets.end());
2016 Facet_vector facets;
2017 facets.reserve(64);
2018 std::unique_copy(all_facets.begin(),
2019 all_facets.end(),
2020 std::back_inserter(facets));
2021 CGAL_HISTOGRAM_PROFILER("|facets|",
2022 static_cast<unsigned int>(facets.size()));
2023 return facets;
2024 }
2025
2026
2027 /**
2028 * Returns false iff a surface facet of `cells` has entered or left the
2029 * restricted Delaunay, or if its surface patch index has changed
2030 *
2031 * That function does not modify the c3t3, but it does update the facet
2032 * surface centers. The function is only called by
2033 * `update_mesh_no_topo_change()`.
2034 */
verify_surface(const Cell_vector & cells)2035 bool verify_surface(const Cell_vector& cells) const
2036 {
2037 // Naive implementation.
2038 // Todo: improve this (maybe we don't have to check if no facet is on surface)
2039 Facet_vector facets = get_facets(cells);
2040 Facet_vector surface_facets;
2041
2042 // Check that nothing changed
2043 Update_c3t3 checker(domain_,c3t3_);
2044 for ( typename Facet_vector::iterator fit = facets.begin() ;
2045 fit != facets.end() ;
2046 ++fit )
2047 {
2048 if ( c3t3_.is_in_complex(*fit) )
2049 {
2050 surface_facets.push_back(*fit);
2051 }
2052 const Surface_patch sp = checker(*fit,
2053 false, /* do not update c3t3 */
2054 true); /* update surface centers */
2055 // false means "do not update the c3t3"
2056 if ( c3t3_.is_in_complex(*fit) != (bool)sp ||
2057 ((bool)sp && !(c3t3_.surface_patch_index(*fit) == *sp) ) )
2058 return false;
2059 }
2060
2061 return true;
2062 }
2063
2064 /**
2065 * Restore mesh for cells and facets of \c cells, using domain_
2066 */
2067 template <typename ForwardIterator>
restore_mesh(ForwardIterator first_cell,ForwardIterator last_cell)2068 void restore_mesh(ForwardIterator first_cell, ForwardIterator last_cell)
2069 {
2070 Facet_vector facets = get_facets(first_cell, last_cell);
2071 restore_mesh(first_cell, last_cell, facets.begin(), facets.end());
2072 }
2073
2074 /**
2075 * Restore mesh for cells of \c cells and facets of \c facets, using domain_
2076 */
2077 template <typename CellForwardIterator, typename FacetForwardIterator>
restore_mesh(CellForwardIterator first_cell,CellForwardIterator last_cell,FacetForwardIterator first_facet,FacetForwardIterator last_facet)2078 void restore_mesh(CellForwardIterator first_cell,
2079 CellForwardIterator last_cell,
2080 FacetForwardIterator first_facet,
2081 FacetForwardIterator last_facet)
2082 {
2083 // Update mesh
2084 Update_c3t3 updater(domain_,c3t3_);
2085 std::for_each(first_facet, last_facet, updater);
2086 std::for_each(first_cell, last_cell, updater);
2087 }
2088
2089 /**
2090 * Returns true if facets of \c facets have the same boundary as
2091 * \c old_boundary, and if the list of vertices has not changed.
2092 */
check_surface_mesh(const Vertex_handle & moving_vertex,const Facet_vector & facets,const Facet_boundary & old_boundary,const Vertex_set & old_incident_surface_vertices)2093 bool check_surface_mesh(const Vertex_handle& moving_vertex,
2094 const Facet_vector& facets,
2095 const Facet_boundary& old_boundary,
2096 const Vertex_set& old_incident_surface_vertices) const
2097 {
2098 Vertex_set incident_surface_vertices;
2099 Facet_boundary new_boundary = get_surface_boundary(moving_vertex,
2100 facets,
2101 incident_surface_vertices);
2102 return ( old_boundary.size() == new_boundary.size() &&
2103 old_incident_surface_vertices == incident_surface_vertices &&
2104 std::equal(new_boundary.begin(),
2105 new_boundary.end(),
2106 old_boundary.begin()) );
2107 }
2108
2109
set_facet_visited(const Facet & facet)2110 void set_facet_visited(const Facet& facet)
2111 {
2112 facet.first->set_facet_visited(facet.second);
2113 const Facet mirror_facet = tr_.mirror_facet(facet);
2114 mirror_facet.first->set_facet_visited(mirror_facet.second);
2115 }
2116
2117 /**
2118 * Orders handles \c h1, \c h2 & \c h3
2119 */
2120 template <typename Handle>
order_handles(Handle & h1,Handle & h2,Handle & h3)2121 void order_handles(Handle& h1, Handle& h2, Handle& h3) const
2122 {
2123 // Function used in get_surface_boundary
2124
2125 if ( h2 < h1 )
2126 std::swap(h1,h2);
2127
2128 if ( h3 < h2 )
2129 {
2130 std::swap(h2,h3);
2131
2132 if ( h2 < h1 ) // don't need to compare h2 & h1 if h2 didn't change
2133 std::swap(h1,h2);
2134 }
2135 }
2136
2137 template <typename CellRange>
reset_sliver_cache(CellRange & cell_range)2138 void reset_sliver_cache(CellRange& cell_range) const
2139 {
2140 reset_sliver_cache(boost::begin(cell_range),
2141 boost::end(cell_range));
2142 }
2143
2144 template <typename CellForwardIterator>
reset_sliver_cache(CellForwardIterator cells_begin,CellForwardIterator cells_end)2145 void reset_sliver_cache(CellForwardIterator cells_begin,
2146 CellForwardIterator cells_end) const
2147 {
2148 while(cells_begin != cells_end) {
2149 (*cells_begin)->reset_cache_validity();
2150 ++cells_begin;
2151 }
2152 }
2153
2154 template <typename CellRange>
reset_circumcenter_cache(CellRange & cell_range)2155 void reset_circumcenter_cache(CellRange& cell_range) const
2156 {
2157 reset_circumcenter_cache(boost::begin(cell_range),
2158 boost::end(cell_range));
2159 }
2160
2161 template <typename CellForwardIterator>
reset_circumcenter_cache(CellForwardIterator cells_begin,CellForwardIterator cells_end)2162 void reset_circumcenter_cache(CellForwardIterator cells_begin,
2163 CellForwardIterator cells_end) const
2164 {
2165 while(cells_begin != cells_end) {
2166 (*cells_begin)->invalidate_weighted_circumcenter_cache();
2167 ++cells_begin;
2168 }
2169 }
2170
2171
2172 private:
2173
2174 // Functor for update_facets function (base)
2175 template <typename C3T3_helpers_, typename FacetUpdater_>
2176 class Update_cell_facets
2177 {
2178 Tr & m_tr;
2179 FacetUpdater_ & m_facet_updater;
2180
2181 protected:
2182
update(const Cell_handle & cell)2183 void update(const Cell_handle& cell) const
2184 {
2185 Cell_handle null_cell;
2186 bool inf = false;
2187 for (int i=0 ; i<4 && (!inf) ; ++i ){
2188 if ( m_tr.is_infinite(cell->vertex(i)) ){
2189 inf = true;
2190 Cell_handle n = cell->neighbor(i);
2191 if(n->next_intrusive() != null_cell){// the neighbor is also outdated
2192 if(cell < n){ // otherwise n will report it later
2193 Facet f(cell,i);
2194 m_facet_updater(f);
2195 }
2196 } else { // report it now or never
2197 if(cell < n){
2198 Facet f(cell,i);
2199 m_facet_updater(f);
2200 }else {
2201 Facet f(n,n->index(cell));
2202 m_facet_updater(f);
2203 }
2204 }
2205 }
2206 }
2207 if(! inf){
2208 for ( int i=0 ; i<4 ; ++i ){
2209 Cell_handle n = cell->neighbor(i);
2210 if(n->next_intrusive() != null_cell){// the neighbor is also outdated
2211 if(cell < n){ // otherwise n will report it later
2212 Facet f(cell,i);
2213 m_facet_updater(f);
2214 }
2215 } else { // report it now or never
2216 if(cell < n){
2217 Facet f(cell,i);
2218 m_facet_updater(f);
2219 }else {
2220 Facet f(n,n->index(cell));
2221 m_facet_updater(f);
2222 }
2223 }
2224 }
2225 }
2226 }
2227
2228 public:
2229 // Constructor
Update_cell_facets(Tr & tr,FacetUpdater_ & fu)2230 Update_cell_facets(Tr &tr,
2231 FacetUpdater_& fu)
2232 : m_tr(tr), m_facet_updater(fu)
2233 {}
2234
2235 // Constructor
Update_cell_facets(const Update_cell_facets & uf)2236 Update_cell_facets(const Update_cell_facets &uf)
2237 : m_tr(uf.m_tr), m_facet_updater(uf.m_facet_updater)
2238 {}
2239
2240 // operator()
operator()2241 void operator()(const Cell_handle& cell) const
2242 {
2243 update(cell);
2244 }
2245 };
2246
2247 #ifdef CGAL_LINKED_WITH_TBB
2248 // Same functor: special version for tbb:parallel_for
2249 template <typename C3T3_helpers_, typename FacetUpdater_>
2250 class Update_cell_facets_for_parallel_for
2251 : Update_cell_facets<C3T3_helpers, FacetUpdater_>
2252 {
2253 typedef Update_cell_facets<C3T3_helpers, FacetUpdater_> Base;
2254 using Base::update;
2255
2256 const std::vector<Cell_handle> & m_outdated_cells;
2257
2258 public:
2259 // Constructor
Update_cell_facets_for_parallel_for(Tr & tr,FacetUpdater_ & fu,const std::vector<Cell_handle> & oc)2260 Update_cell_facets_for_parallel_for(
2261 Tr &tr,
2262 FacetUpdater_& fu,
2263 const std::vector<Cell_handle> &oc)
2264 : Base(tr, fu), m_outdated_cells(oc)
2265 {}
2266
2267 // Constructor
Update_cell_facets_for_parallel_for(const Update_cell_facets_for_parallel_for & uf)2268 Update_cell_facets_for_parallel_for(
2269 const Update_cell_facets_for_parallel_for &uf)
2270 : Base(uf), m_outdated_cells(uf.m_outdated_cells)
2271 {}
2272
2273 // operator()
operator()2274 void operator()(const tbb::blocked_range<size_t>& r) const
2275 {
2276 for( size_t i = r.begin() ; i != r.end() ; ++i)
2277 update(m_outdated_cells[i]);
2278 }
2279 };
2280 #endif //CGAL_LINKED_WITH_TBB
2281
2282 // -----------------------------------
2283 // -----------------------------------
2284 // -----------------------------------
2285
2286 // Functor for rebuild_restricted_delaunay function
2287 template <typename C3T3_, typename Update_c3t3_>
2288 class Update_cell
2289 {
2290 C3T3 & m_c3t3;
2291 Update_c3t3_ & m_updater;
2292
2293 protected:
2294 typedef typename C3T3_::Cell_handle Cell_handle;
2295
update(const Cell_handle & cell)2296 void update(const Cell_handle& cell) const
2297 {
2298 m_c3t3.remove_from_complex(cell);
2299 m_updater(cell);
2300 }
2301
2302 public:
2303 // Constructor
Update_cell(C3T3_ & c3t3,Update_c3t3_ & updater)2304 Update_cell(C3T3_ &c3t3, Update_c3t3_& updater)
2305 : m_c3t3(c3t3), m_updater(updater)
2306 {}
2307
2308 // Constructor
Update_cell(const Update_cell & uc)2309 Update_cell(const Update_cell &uc)
2310 : m_c3t3(uc.m_c3t3), m_updater(uc.m_updater)
2311 {}
2312
2313 // operator()
operator()2314 void operator()(const Cell_handle& cell) const
2315 {
2316 update(cell);
2317 }
2318 };
2319
2320 #ifdef CGAL_LINKED_WITH_TBB
2321 // Same functor: special version for tbb:parallel_for
2322 template <typename C3T3_, typename Update_c3t3_>
2323 class Update_cell_for_parallel_for
2324 : Update_cell<C3T3_, Update_c3t3_>
2325 {
2326 typedef Update_cell<C3T3_, Update_c3t3_> Base;
2327 using Base::update;
2328
2329 const std::vector<Cell_handle> & m_outdated_cells;
2330
2331 public:
2332 // Constructor
Update_cell_for_parallel_for(C3T3_ & c3t3,Update_c3t3_ & updater,const std::vector<Cell_handle> & oc)2333 Update_cell_for_parallel_for(
2334 C3T3_ &c3t3,
2335 Update_c3t3_& updater,
2336 const std::vector<Cell_handle> &oc)
2337 : Base(c3t3, updater), m_outdated_cells(oc)
2338 {}
2339
2340 // Constructor
Update_cell_for_parallel_for(const Update_cell_for_parallel_for & uc)2341 Update_cell_for_parallel_for(
2342 const Update_cell_for_parallel_for &uc)
2343 : Base(uc), m_outdated_cells(uc.m_outdated_cells)
2344 {}
2345
2346 // operator()
operator()2347 void operator()(const tbb::blocked_range<size_t>& r) const
2348 {
2349 for( size_t i = r.begin() ; i != r.end() ; ++i)
2350 update(m_outdated_cells[i]);
2351 }
2352 };
2353 #endif //CGAL_LINKED_WITH_TBB
2354
2355 // -----------------------------------
2356 // -----------------------------------
2357 // -----------------------------------
2358
2359 // Functor for rebuild_restricted_delaunay function
2360 #ifdef CGAL_LINKED_WITH_TBB
2361 template <typename C3T3_helpers_, typename C3T3_, typename Update_c3t3_,
2362 typename Vertex_to_proj_set_>
2363 class Update_facet
2364 {
2365 const C3T3_helpers_ & m_c3t3_helpers;
2366 C3T3_ & m_c3t3;
2367 Update_c3t3_ & m_updater;
2368 Vertex_to_proj_set_ & m_vertex_to_proj;
2369
2370 typedef typename C3T3_::Vertex_handle Vertex_handle;
2371 typedef typename C3T3_::Cell_handle Cell_handle;
2372 typedef typename C3T3_::Facet Facet;
2373 typedef typename C3T3::Surface_patch_index Surface_patch_index;
2374
2375 public:
2376 // Constructor
Update_facet(const C3T3_helpers_ & c3t3_helpers,C3T3_ & c3t3,Update_c3t3_ & updater,Vertex_to_proj_set_ & vertex_to_proj)2377 Update_facet(const C3T3_helpers_ & c3t3_helpers,
2378 C3T3_ &c3t3, Update_c3t3_& updater,
2379 Vertex_to_proj_set_ &vertex_to_proj)
2380 : m_c3t3_helpers(c3t3_helpers), m_c3t3(c3t3), m_updater(updater),
2381 m_vertex_to_proj(vertex_to_proj)
2382 {}
2383
2384 // Constructor
Update_facet(const Update_facet & uc)2385 Update_facet(const Update_facet &uc)
2386 : m_c3t3_helpers(uc.m_c3t3_helpers), m_c3t3(uc.m_c3t3),
2387 m_updater(uc.m_updater), m_vertex_to_proj(uc.m_vertex_to_proj)
2388 {}
2389
2390 // operator()
operator()2391 void operator()( const Facet& facet ) const
2392 {
2393 // Update facet
2394 m_c3t3.remove_from_complex(facet);
2395 m_updater(facet);
2396
2397 // Update m_vertex_to_proj
2398 if ( m_c3t3.is_in_complex(facet) )
2399 {
2400 // Iterate on vertices
2401 int k = facet.second;
2402 for ( int i=1 ; i<4 ; ++i )
2403 {
2404 const Vertex_handle& v = facet.first->vertex((k+i)&3);
2405 if ( m_c3t3.in_dimension(v) > 2 )
2406 {
2407 std::pair<Vertex_handle, Surface_patch_index> p
2408 = std::make_pair(v, m_c3t3.surface_patch_index(facet));
2409 m_c3t3_helpers.lock_vertex_to_proj();
2410 m_vertex_to_proj.insert(p);
2411 m_c3t3_helpers.unlock_vertex_to_proj();
2412 }
2413 }
2414 }
2415 }
2416 };
2417 #endif
2418
2419 // -----------------------------------
2420 // Private data
2421 // -----------------------------------
2422 C3T3& c3t3_;
2423 Tr& tr_;
2424 const MeshDomain& domain_;
2425 }; // class C3T3_helpers
2426
2427
2428 template <typename C3T3, typename MD>
2429 template <typename SliverCriterion, typename OutputIterator>
2430 std::pair<bool,typename C3T3_helpers<C3T3,MD>::Vertex_handle>
2431 C3T3_helpers<C3T3,MD>::
update_mesh(const Vertex_handle & old_vertex,const Vector_3 & move,const Weighted_point & new_position,const SliverCriterion & criterion,OutputIterator modified_vertices,bool * could_lock_zone)2432 update_mesh(const Vertex_handle& old_vertex,
2433 const Vector_3& move,
2434 const Weighted_point& new_position,
2435 const SliverCriterion& criterion,
2436 OutputIterator modified_vertices,
2437 bool *could_lock_zone)
2438 {
2439 // std::cerr << "\nupdate_mesh[v1](" << (void*)(&*old_vertex)
2440 // << "=" << tr_.point(old_vertex) << ",\n"
2441 // << " " << move << ",\n"
2442 // << " " << new_position << ")\n";
2443
2444 if (could_lock_zone)
2445 *could_lock_zone = true;
2446
2447 Cell_vector incident_cells_;
2448 incident_cells_.reserve(64);
2449 tr_.incident_cells(old_vertex, std::back_inserter(incident_cells_));
2450 if ( Th().no_topological_change(tr_, old_vertex, move, new_position, incident_cells_) )
2451 {
2452 return update_mesh_no_topo_change(old_vertex, move, new_position, criterion,
2453 modified_vertices, incident_cells_);
2454 }
2455 else
2456 {
2457 return update_mesh_topo_change(old_vertex, new_position, criterion,
2458 modified_vertices, could_lock_zone);
2459 }
2460 }
2461
2462 template <typename C3T3, typename MD>
2463 template <typename SliverCriterion, typename OutputIterator>
2464 std::pair<bool,typename C3T3_helpers<C3T3,MD>::Vertex_handle>
2465 C3T3_helpers<C3T3,MD>::
update_mesh_no_topo_change(const Vertex_handle & old_vertex,const Vector_3 & move,const Weighted_point & new_position,const SliverCriterion & criterion,OutputIterator modified_vertices,const Cell_vector & conflict_cells)2466 update_mesh_no_topo_change(const Vertex_handle& old_vertex,
2467 const Vector_3& move,
2468 const Weighted_point& new_position,
2469 const SliverCriterion& criterion,
2470 OutputIterator modified_vertices,
2471 const Cell_vector& conflict_cells )
2472 {
2473 #ifdef CGAL_MESH_3_C3T3_HELPERS_VERBOSE
2474 std::cerr << "update_mesh_no_topo_change("
2475 << (void*)(&*old_vertex) << " = " << tr_.point(old_vertex) << ",\n"
2476 << " " << move << ",\n"
2477 << " " << new_position << ")" << std::endl;
2478 #endif
2479
2480 typename Gt::Construct_opposite_vector_3 cov = tr_.geom_traits().construct_opposite_vector_3_object();
2481
2482 //backup metadata
2483 std::set<Cell_data_backup> cells_backup;
2484 fill_cells_backup(conflict_cells, cells_backup);
2485
2486 // Get old values
2487 criterion.before_move(c3t3_cells(conflict_cells));
2488 Weighted_point old_position = tr_.point(old_vertex); // intentional copy
2489
2490 // Move point
2491 reset_circumcenter_cache(conflict_cells);
2492 reset_sliver_cache(conflict_cells);
2493 move_point_no_topo_change(old_vertex, move, new_position);
2494
2495 // Check that surface mesh is still valid
2496 // and Get new criterion value (conflict_zone did not change)
2497 // warnings : valid_move updates caches
2498 // verify_surface does not change c3t3 when returns false,
2499 // but it does change circumcenters
2500 if( verify_surface(conflict_cells)
2501 && criterion.valid_move(c3t3_cells(conflict_cells)))
2502 {
2503 fill_modified_vertices(conflict_cells.begin(), conflict_cells.end(),
2504 old_vertex, modified_vertices);
2505 return std::make_pair(true,old_vertex);
2506 }
2507 else // revert move
2508 {
2509 #ifdef CGAL_MESH_3_C3T3_HELPERS_VERBOSE
2510 std::cerr << "update_mesh_no_topo_change: revert move to "
2511 << old_position << "\n";
2512 #endif
2513 reset_circumcenter_cache(conflict_cells);
2514
2515 // Sliver caches have been updated by valid_move
2516 reset_sliver_cache(conflict_cells);
2517 move_point_no_topo_change(old_vertex, cov(move), old_position);
2518
2519 // Restore meta-data (cells should have same connectivity as before move)
2520 // cells_backup does not contain infinite cells so they can be fewer
2521 CGAL_assertion(conflict_cells.size() >= cells_backup.size());
2522 restore_from_cells_backup(conflict_cells, cells_backup);
2523
2524 return std::make_pair(false, old_vertex);
2525 }
2526 }
2527
2528
2529 template <typename C3T3, typename MD>
2530 template <typename SliverCriterion, typename OutputIterator>
2531 std::pair<bool,typename C3T3_helpers<C3T3,MD>::Vertex_handle>
2532 C3T3_helpers<C3T3,MD>::
update_mesh_topo_change(const Vertex_handle & old_vertex,const Weighted_point & new_position,const SliverCriterion & criterion,OutputIterator modified_vertices,bool * could_lock_zone)2533 update_mesh_topo_change(const Vertex_handle& old_vertex,
2534 const Weighted_point& new_position,
2535 const SliverCriterion& criterion,
2536 OutputIterator modified_vertices,
2537 bool *could_lock_zone)
2538 {
2539 #ifdef CGAL_MESH_3_C3T3_HELPERS_VERBOSE
2540 std::cerr << "update_mesh_topo_change(" << (void*)(&*old_vertex)
2541 << "=" << tr_.point(old_vertex)
2542 << " " << new_position << ",\n"
2543 << ")" << std::endl;
2544 #endif
2545
2546 Cell_set insertion_conflict_cells;
2547 Cell_set removal_conflict_cells;
2548 Facet_vector insertion_conflict_boundary;
2549 insertion_conflict_boundary.reserve(64);
2550 get_conflict_zone_topo_change(old_vertex, new_position,
2551 std::inserter(insertion_conflict_cells,insertion_conflict_cells.end()),
2552 std::back_inserter(insertion_conflict_boundary),
2553 std::inserter(removal_conflict_cells, removal_conflict_cells.end()),
2554 could_lock_zone);
2555
2556 if (could_lock_zone && *could_lock_zone == false)
2557 return std::make_pair(false, Vertex_handle());
2558
2559 if(insertion_conflict_boundary.empty())
2560 return std::make_pair(false, old_vertex); // new_location is a vertex already
2561
2562 Cell_vector conflict_cells;
2563 conflict_cells.reserve(insertion_conflict_cells.size()+removal_conflict_cells.size());
2564 std::set_union(insertion_conflict_cells.begin(), insertion_conflict_cells.end(),
2565 removal_conflict_cells.begin(), removal_conflict_cells.end(),
2566 std::back_inserter(conflict_cells));
2567
2568 // Backup metadata
2569 std::set<Cell_data_backup> cells_backup;
2570 fill_cells_backup(conflict_cells, cells_backup);
2571 CGAL_assertion(conflict_cells.size() == cells_backup.size());
2572
2573 criterion.before_move(c3t3_cells(conflict_cells));
2574 Weighted_point old_position = tr_.point(old_vertex); // intentional copy
2575
2576 // Keep old boundary
2577 Vertex_set old_incident_surface_vertices;
2578 Facet_boundary old_surface_boundary =
2579 get_surface_boundary(old_vertex, conflict_cells, old_incident_surface_vertices);
2580
2581 reset_circumcenter_cache(conflict_cells);
2582 reset_sliver_cache(conflict_cells);
2583
2584 Cell_vector outdated_cells;
2585 outdated_cells.reserve(64);
2586 Vertex_handle new_vertex =
2587 move_point_topo_change_conflict_zone_known(old_vertex, new_position,
2588 insertion_conflict_boundary[0],
2589 insertion_conflict_cells.begin(),
2590 insertion_conflict_cells.end(),
2591 removal_conflict_cells.begin(),
2592 removal_conflict_cells.end(),
2593 std::back_inserter(outdated_cells),
2594 CGAL::Emptyset_iterator());
2595
2596 // If nothing changed, return
2597 if ( old_position == tr_.point(new_vertex) )
2598 {
2599 #ifdef CGAL_MESH_3_C3T3_HELPERS_VERBOSE
2600 std::cerr << "update_mesh_topo_change: no move!\n";
2601 #endif
2602 return std::make_pair(false,old_vertex);
2603 }
2604
2605 restore_mesh(outdated_cells.begin(),outdated_cells.end());
2606
2607 // Check that surface boundary does not change.
2608 // This check ensures that vertices which are inside c3t3 stay inside.
2609 if (criterion.valid_move(c3t3_cells(outdated_cells))
2610 && check_surface_mesh(new_vertex,
2611 get_facets(outdated_cells),
2612 old_surface_boundary,
2613 old_incident_surface_vertices) )
2614 {
2615 fill_modified_vertices(outdated_cells.begin(), outdated_cells.end(),
2616 new_vertex, modified_vertices);
2617 return std::make_pair(true,new_vertex);
2618 }
2619 else
2620 {
2621 // Removing from c3t3 cells which will be destroyed by revert_move
2622 // is done by move_point_topo_change_conflict_zone_known, called by revert_move
2623
2624 #ifdef CGAL_MESH_3_C3T3_HELPERS_VERBOSE
2625 std::cerr << "update_mesh_topo_change: revert move to "
2626 << old_position << "\n";
2627 #endif
2628
2629 //reset caches in case cells are re-used by the compact container
2630 reset_circumcenter_cache(outdated_cells);
2631 reset_sliver_cache(outdated_cells);
2632 outdated_cells.clear();
2633
2634 // Revert move
2635 Vertex_handle revert_vertex = revert_move(new_vertex, old_position,
2636 std::inserter(outdated_cells, outdated_cells.end()));
2637
2638 //restore meta-data (cells should have same connectivity as before move)
2639 //cells should be the same (connectivity-wise) as before initial move
2640 CGAL_assertion(outdated_cells.size() == cells_backup.size());
2641 restore_from_cells_backup(outdated_cells, cells_backup);
2642
2643 return std::make_pair(false,revert_vertex);
2644 }
2645 }
2646
2647 template <typename C3T3, typename MD>
2648 template <typename OutputIterator>
2649 typename C3T3_helpers<C3T3,MD>::Vertex_handle
2650 C3T3_helpers<C3T3,MD>::
update_mesh(const Vertex_handle & old_vertex,const Vector_3 & move,OutputIterator modified_vertices,bool fill_vertices)2651 update_mesh(const Vertex_handle& old_vertex,
2652 const Vector_3& move,
2653 OutputIterator modified_vertices,
2654 bool fill_vertices)
2655 {
2656 #ifdef CGAL_MESH_3_C3T3_HELPERS_VERBOSE
2657 std::cerr << "\nupdate_mesh[v2](" << (void*)(&*old_vertex)
2658 << "=" << tr_.point(old_vertex) << ",\n"
2659 << " " << move << ")\n";
2660 #endif
2661
2662 Cell_vector outdated_cells;
2663 Vertex_handle new_vertex = move_point(old_vertex, move,
2664 std::back_inserter(outdated_cells),
2665 CGAL::Emptyset_iterator());
2666 // move_point has invalidated caches
2667
2668 restore_mesh(outdated_cells.begin(), outdated_cells.end());
2669
2670 // Fill modified vertices
2671 if ( fill_vertices
2672 && !(boost::is_same<OutputIterator,CGAL::Emptyset_iterator>::value))
2673 {
2674 fill_modified_vertices(outdated_cells.begin(), outdated_cells.end(),
2675 new_vertex, modified_vertices);
2676 }
2677
2678 return new_vertex;
2679 }
2680
2681 #ifdef CGAL_INTRUSIVE_LIST
2682 template <typename C3T3, typename MD>
2683 template <typename OutdatedCells>
2684 void
2685 C3T3_helpers<C3T3,MD>::
rebuild_restricted_delaunay(OutdatedCells & outdated_cells,Moving_vertices_set & moving_vertices)2686 rebuild_restricted_delaunay(OutdatedCells& outdated_cells,
2687 Moving_vertices_set& moving_vertices)
2688 {
2689 typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object();
2690
2691 typename OutdatedCells::iterator first_cell = outdated_cells.begin();
2692 typename OutdatedCells::iterator last_cell = outdated_cells.end();
2693 Update_c3t3 updater(domain_,c3t3_);
2694
2695 # ifdef CGAL_MESH_3_PROFILING
2696 std::cerr << std::endl << " Updating cells...";
2697 WallClockTimer t;
2698 size_t num_cells = c3t3_.number_of_cells_in_complex();
2699 # endif
2700
2701 // Updates cells
2702 // Note: ~58% of rebuild_restricted_delaunay time
2703
2704 Vertex_set vertex_to_proj;
2705
2706 # ifdef CGAL_LINKED_WITH_TBB
2707 // Parallel
2708 if (boost::is_convertible<Concurrency_tag, Parallel_tag>::value)
2709 {
2710 std::vector<Cell_handle> outdated_cells_vector;
2711 outdated_cells_vector.reserve(outdated_cells.size());
2712 for ( ; first_cell != last_cell ; ++first_cell)
2713 {
2714 outdated_cells_vector.push_back(*first_cell);
2715 }
2716
2717 tbb::parallel_for(
2718 tbb::blocked_range<size_t>(0, outdated_cells_vector.size()),
2719 Update_cell_for_parallel_for<C3T3, Update_c3t3>(
2720 c3t3_, updater, outdated_cells_vector));
2721
2722 # ifdef CGAL_MESH_3_PROFILING
2723 std::cerr << " done in " << t.elapsed() << " seconds (#cells from "
2724 << num_cells << " to " << c3t3_.number_of_cells_in_complex() << ")."
2725 << std::endl;
2726 std::cerr << " Updating facets...";
2727 t.reset();
2728 # endif
2729
2730 // Get facets (returns each canonical facet only once)
2731 // Note: ~42% of rebuild_restricted_delaunay time
2732 // Facet_vector facets;
2733 lock_vertex_to_proj();
2734 Facet_updater facet_updater(*this, c3t3_,vertex_to_proj, updater);
2735 unlock_vertex_to_proj();
2736 update_facets(outdated_cells_vector, facet_updater);
2737
2738 // now we can clear
2739 outdated_cells.clear();
2740 }
2741 // Sequential
2742 else
2743 # endif // CGAL_LINKED_WITH_TBB
2744 {
2745 while ( first_cell != last_cell )
2746 {
2747 Cell_handle cell = *first_cell++;
2748 c3t3_.remove_from_complex(cell);
2749 updater(cell);
2750 }
2751
2752 # ifdef CGAL_MESH_3_PROFILING
2753 std::cerr << " done in " << t.elapsed() << " seconds (#cells from "
2754 << num_cells << " to " << c3t3_.number_of_cells_in_complex() << ")."
2755 << std::endl;
2756 std::cerr << " Updating facets...";
2757 t.reset();
2758 # endif
2759
2760 // Get facets (returns each canonical facet only once)
2761 // Note: ~42% of rebuild_restricted_delaunay time
2762 // Facet_vector facets;
2763 Facet_updater facet_updater(*this, c3t3_,vertex_to_proj, updater);
2764 update_facets(outdated_cells, facet_updater);
2765
2766 // now we can clear
2767 outdated_cells.clear();
2768 }
2769
2770 # ifdef CGAL_MESH_3_PROFILING
2771 std::cerr << " done in " << t.elapsed() << " seconds ("
2772 << vertex_to_proj.size() << " vertices to project)." << std::endl;
2773 std::cerr << " Projecting interior vertices...";
2774 t.reset();
2775 # endif
2776
2777 CGAL_HISTOGRAM_PROFILER("|vertex_to_proj|=",
2778 static_cast<unsigned int>(vertex_to_proj.size()));
2779
2780 // Project interior vertices
2781 // Note: ~0% of rebuild_restricted_delaunay time
2782 // TODO : iterate to be sure no interior vertice become on the surface
2783 // because of move ?
2784 for ( typename Vertex_set::iterator it = vertex_to_proj.begin() ;
2785 it != vertex_to_proj.end() ;
2786 ++it )
2787 {
2788 const Weighted_point& initial_position = tr_.point(*it);
2789 boost::optional<Bare_point> opt_new_pos = project_on_surface(*it, cp(initial_position));
2790
2791 if ( opt_new_pos )
2792 {
2793 //freezing needs 'erase' to be done before the vertex is actually destroyed
2794 // Update moving vertices (it becomes new_vertex)
2795 moving_vertices.erase(*it);
2796
2797 Vector_3 move(cp(initial_position), *opt_new_pos);
2798 Vertex_handle new_vertex = update_mesh(*it, move);
2799 c3t3_.set_dimension(new_vertex, 2);
2800
2801 moving_vertices.insert(new_vertex);
2802 }
2803 }
2804
2805 # ifdef CGAL_MESH_3_PROFILING
2806 std::cerr << " done in " << t.elapsed() << " seconds." << std::endl;
2807 # endif
2808 }
2809 #endif //CGAL_INTRUSIVE_LIST
2810
2811 template <typename C3T3, typename MD>
2812 template <typename ForwardIterator>
2813 void
2814 C3T3_helpers<C3T3,MD>::
rebuild_restricted_delaunay(ForwardIterator first_cell,ForwardIterator last_cell,Moving_vertices_set & moving_vertices)2815 rebuild_restricted_delaunay(ForwardIterator first_cell,
2816 ForwardIterator last_cell,
2817 Moving_vertices_set& moving_vertices)
2818 {
2819 typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object();
2820 typename Gt::Construct_vector_3 vector = tr_.geom_traits().construct_vector_3_object();
2821 typename Gt::Equal_3 equal = tr_.geom_traits().equal_3_object();
2822
2823 Update_c3t3 updater(domain_,c3t3_);
2824
2825 // Get facets (returns each canonical facet only once)
2826 Facet_vector facets = get_facets(first_cell, last_cell);
2827
2828 // Updates cells
2829 // Note: ~58% of rebuild_restricted_delaunay time
2830 #ifdef CGAL_LINKED_WITH_TBB
2831 // Parallel
2832 if (boost::is_convertible<Concurrency_tag, Parallel_tag>::value)
2833 {
2834 tbb::parallel_do(first_cell, last_cell,
2835 Update_cell<C3T3, Update_c3t3>(c3t3_, updater));
2836 }
2837 // Sequential
2838 else
2839 #endif // CGAL_LINKED_WITH_TBB
2840 {
2841 while ( first_cell != last_cell )
2842 {
2843 Update_cell<C3T3, Update_c3t3> uc(c3t3_, updater);
2844 uc(*first_cell++);
2845 }
2846 }
2847
2848 // Updates facets
2849 typedef std::set<std::pair<Vertex_handle, Surface_patch_index> > Vertex_to_proj_set;
2850 Vertex_to_proj_set vertex_to_proj;
2851
2852 #ifdef CGAL_LINKED_WITH_TBB
2853 // Parallel
2854 if (boost::is_convertible<Concurrency_tag, Parallel_tag>::value)
2855 {
2856 tbb::parallel_do(
2857 facets.begin(), facets.end(),
2858 Update_facet<Self, C3T3, Update_c3t3, Vertex_to_proj_set>(
2859 *this, c3t3_, updater, vertex_to_proj)
2860 );
2861 }
2862 // Sequential
2863 else
2864 #endif // CGAL_LINKED_WITH_TBB
2865 {
2866 for ( typename Facet_vector::iterator fit = facets.begin() ;
2867 fit != facets.end() ;
2868 ++fit )
2869 {
2870 // Update facet
2871 c3t3_.remove_from_complex(*fit);
2872 updater(*fit);
2873
2874 // Update vertex_to_proj
2875 if ( c3t3_.is_in_complex(*fit) )
2876 {
2877 // Iterate on vertices
2878 int k = fit->second;
2879 for ( int i=1 ; i<4 ; ++i )
2880 {
2881 const Vertex_handle& v = fit->first->vertex((k+i)&3);
2882 if ( c3t3_.in_dimension(v) > 2 )
2883 {
2884 vertex_to_proj.insert(
2885 std::make_pair(v, c3t3_.surface_patch_index(*fit)));
2886 }
2887 }
2888 }
2889 }
2890 }
2891
2892 // Project interior vertices
2893 // TODO : iterate to be sure no interior vertice become on the surface
2894 // because of move ?
2895 typename Vertex_to_proj_set::iterator it = vertex_to_proj.begin(),
2896 end = vertex_to_proj.end();
2897 for ( ; it != end; ++it )
2898 {
2899 Vertex_handle vh = it->first;
2900 const Weighted_point& initial_position = tr_.point(vh);
2901 boost::optional<Bare_point> opt_new_pos = project_on_surface(vh, cp(initial_position), it->second);
2902
2903 if ( opt_new_pos )
2904 {
2905 //freezing needs 'erase' to be done before the vertex is actually destroyed
2906 // Update moving vertices (it becomes new_vertex)
2907 moving_vertices.erase(vh);
2908
2909 Vertex_handle new_vertex = update_mesh(vh, vector(cp(initial_position), *opt_new_pos));
2910 c3t3_.set_dimension(new_vertex, 2);
2911
2912 moving_vertices.insert(new_vertex);
2913 }
2914 }
2915 }
2916
2917 template <typename C3T3, typename MD>
2918 void
2919 C3T3_helpers<C3T3, MD>::
update_restricted_facets()2920 update_restricted_facets()
2921 {
2922 Update_c3t3 updater(domain_, c3t3_);
2923 for (typename C3T3::Triangulation::Finite_facets_iterator
2924 fit = tr_.finite_facets_begin();
2925 fit != tr_.finite_facets_end(); ++fit)
2926 updater(*fit);
2927 }
2928
2929 template <typename C3T3, typename MD>
2930 template <typename OutdatedCellsOutputIterator,
2931 typename DeletedCellsOutputIterator>
2932 typename C3T3_helpers<C3T3,MD>::Vertex_handle
2933 C3T3_helpers<C3T3,MD>::
move_point(const Vertex_handle & old_vertex,const Vector_3 & move,OutdatedCellsOutputIterator outdated_cells,DeletedCellsOutputIterator deleted_cells)2934 move_point(const Vertex_handle& old_vertex,
2935 const Vector_3& move,
2936 OutdatedCellsOutputIterator outdated_cells,
2937 DeletedCellsOutputIterator deleted_cells) const
2938 {
2939 #ifdef CGAL_MESH_3_C3T3_HELPERS_VERBOSE
2940 std::cerr << "C3T3_helpers::move_point[v2](" << (void*)(&*old_vertex)
2941 << " = " << tr_.point(old_vertex) << ",\n"
2942 << " " << move << ")\n";
2943 #endif
2944
2945 typename Gt::Construct_translated_point_3 translate = tr_.geom_traits().construct_translated_point_3_object();
2946 typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object();
2947 typename Gt::Construct_weighted_point_3 cwp = tr_.geom_traits().construct_weighted_point_3_object();
2948
2949 Cell_vector incident_cells_;
2950 incident_cells_.reserve(64);
2951
2952 # ifdef CGAL_LINKED_WITH_TBB
2953 // Parallel
2954 if (boost::is_convertible<Concurrency_tag, Parallel_tag>::value)
2955 {
2956 tr_.incident_cells_threadsafe(old_vertex, std::back_inserter(incident_cells_));
2957 }
2958 // Sequential
2959 else
2960 # endif // CGAL_LINKED_WITH_TBB
2961 {
2962 tr_.incident_cells(old_vertex, std::back_inserter(incident_cells_));
2963 }
2964
2965 const Weighted_point& position = tr_.point(old_vertex);
2966 const Weighted_point& new_position = cwp(translate(cp(position), move));
2967
2968 if ( Th().no_topological_change(tr_, old_vertex, move, new_position, incident_cells_) )
2969 {
2970 reset_circumcenter_cache(incident_cells_);
2971 reset_sliver_cache(incident_cells_);
2972 std::copy(incident_cells_.begin(),incident_cells_.end(), outdated_cells);
2973 return move_point_no_topo_change(old_vertex, move, new_position);
2974 }
2975 else
2976 {
2977 return move_point_topo_change(old_vertex, new_position, outdated_cells, deleted_cells);
2978 }
2979 }
2980
2981 // Sequential
2982 template <typename C3T3, typename MD>
2983 typename C3T3_helpers<C3T3,MD>::Vertex_handle
2984 C3T3_helpers<C3T3,MD>::
move_point(const Vertex_handle & old_vertex,const Vector_3 & move,Outdated_cell_set & outdated_cells_set,Moving_vertices_set & moving_vertices)2985 move_point(const Vertex_handle& old_vertex,
2986 const Vector_3& move,
2987 Outdated_cell_set& outdated_cells_set,
2988 Moving_vertices_set& moving_vertices) const
2989 {
2990 #ifdef CGAL_MESH_3_C3T3_HELPERS_VERBOSE
2991 std::cerr << "C3T3_helpers::move_point(" << (void*)(&*old_vertex)
2992 << " = " << tr_.point(old_vertex) << ",\n"
2993 << " " << move << ")\n";
2994 #endif
2995
2996 typename Gt::Construct_translated_point_3 translate = tr_.geom_traits().construct_translated_point_3_object();
2997 typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object();
2998 typename Gt::Construct_weighted_point_3 cwp = tr_.geom_traits().construct_weighted_point_3_object();
2999
3000 Cell_vector incident_cells_;
3001 incident_cells_.reserve(64);
3002 tr_.incident_cells(old_vertex, std::back_inserter(incident_cells_));
3003
3004 const Weighted_point& position = tr_.point(old_vertex);
3005 const Weighted_point& new_position = cwp(translate(cp(position), move));
3006
3007 if ( Th().no_topological_change(tr_, old_vertex, move, new_position, incident_cells_) )
3008 {
3009 reset_circumcenter_cache(incident_cells_);
3010 reset_sliver_cache(incident_cells_);
3011 std::copy(incident_cells_.begin(),incident_cells_.end(),
3012 std::inserter(outdated_cells_set, outdated_cells_set.end()));
3013 return move_point_no_topo_change(old_vertex, move, new_position);
3014 }
3015 else
3016 {
3017 moving_vertices.erase(old_vertex);
3018
3019 Vertex_handle new_vertex = move_point_topo_change(old_vertex, new_position, outdated_cells_set);
3020
3021 moving_vertices.insert(new_vertex);
3022 return new_vertex;
3023 }
3024 }
3025
3026 // Parallel
3027 // In case of success (could_lock_zone = true), the zone is locked after the call
3028 // ==> the caller needs to call "unlock_all_elements" by itself
3029 // In case of failure (could_lock_zone = false), the zone is unlocked by this function
3030 template <typename C3T3, typename MD>
3031 typename C3T3_helpers<C3T3,MD>::Vertex_handle
3032 C3T3_helpers<C3T3,MD>::
move_point(const Vertex_handle & old_vertex,const Vector_3 & move,Outdated_cell_set & outdated_cells_set,Moving_vertices_set & moving_vertices,bool * could_lock_zone)3033 move_point(const Vertex_handle& old_vertex,
3034 const Vector_3& move,
3035 Outdated_cell_set& outdated_cells_set,
3036 Moving_vertices_set& moving_vertices,
3037 bool *could_lock_zone) const
3038 {
3039 CGAL_assertion(could_lock_zone != nullptr);
3040 *could_lock_zone = true;
3041
3042 if (!try_lock_vertex(old_vertex)) // LOCK
3043 {
3044 *could_lock_zone = false;
3045 unlock_all_elements();
3046 return Vertex_handle();
3047 }
3048
3049 //======= Get incident cells ==========
3050 Cell_vector incident_cells_;
3051 incident_cells_.reserve(64);
3052 if (tr_.try_lock_and_get_incident_cells(old_vertex, incident_cells_) == false)
3053 {
3054 *could_lock_zone = false;
3055 unlock_all_elements();
3056 return Vertex_handle();
3057 }
3058 //======= /Get incident cells ==========
3059
3060 typename Gt::Construct_translated_point_3 translate = tr_.geom_traits().construct_translated_point_3_object();
3061 typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object();
3062 typename Gt::Construct_weighted_point_3 cwp = tr_.geom_traits().construct_weighted_point_3_object();
3063
3064 const Weighted_point& position = tr_.point(old_vertex);
3065 const Weighted_point& new_position = cwp(translate(cp(position), move));
3066
3067 if (!try_lock_point(new_position)) // LOCK
3068 {
3069 *could_lock_zone = false;
3070 unlock_all_elements();
3071 return Vertex_handle();
3072 }
3073
3074 if ( Th().no_topological_change(tr_, old_vertex, move, new_position, incident_cells_) )
3075 {
3076 reset_circumcenter_cache(incident_cells_);
3077 reset_sliver_cache(incident_cells_);
3078
3079 lock_outdated_cells();
3080 std::copy(incident_cells_.begin(),incident_cells_.end(),
3081 std::inserter(outdated_cells_set, outdated_cells_set.end()));
3082
3083 Vertex_handle new_vertex =
3084 move_point_no_topo_change(old_vertex, move, new_position);
3085
3086 unlock_outdated_cells();
3087
3088 // Don't "unlock_all_elements" here, the caller may need it to do it himself
3089 return new_vertex;
3090 }
3091 else
3092 {
3093 //moving_vertices.erase(old_vertex); MOVED BELOW
3094
3095 Vertex_handle new_vertex =
3096 move_point_topo_change(old_vertex, new_position, outdated_cells_set,
3097 could_lock_zone);
3098
3099 if (*could_lock_zone == false)
3100 {
3101 unlock_all_elements();
3102 return Vertex_handle();
3103 }
3104
3105 lock_moving_vertices();
3106 moving_vertices.erase(old_vertex);
3107 moving_vertices.insert(new_vertex);
3108 unlock_moving_vertices();
3109
3110 // Don't "unlock_all_elements" here, the caller may need it to do it himself
3111 return new_vertex;
3112 }
3113 }
3114
3115 template <typename C3T3, typename MD>
3116 typename C3T3_helpers<C3T3,MD>::Vertex_handle
3117 C3T3_helpers<C3T3,MD>::
move_point_topo_change(const Vertex_handle & old_vertex,const Weighted_point & new_position,Outdated_cell_set & outdated_cells_set,bool * could_lock_zone)3118 move_point_topo_change(const Vertex_handle& old_vertex,
3119 const Weighted_point& new_position,
3120 Outdated_cell_set& outdated_cells_set,
3121 bool *could_lock_zone) const
3122 {
3123 #ifdef CGAL_MESH_3_C3T3_HELPERS_VERBOSE
3124 std::cerr << "move_point_topo_change(" << (void*)(&*old_vertex)
3125 << "=" << tr_.point(old_vertex)
3126 << " " << new_position << ",\n"
3127 << ")" << std::endl;
3128 #endif
3129
3130 Cell_set insertion_conflict_cells;
3131 Cell_set removal_conflict_cells;
3132 Facet_vector insertion_conflict_boundary;
3133 insertion_conflict_boundary.reserve(64);
3134
3135 get_conflict_zone_topo_change(old_vertex, new_position,
3136 std::inserter(insertion_conflict_cells,insertion_conflict_cells.end()),
3137 std::back_inserter(insertion_conflict_boundary),
3138 std::inserter(removal_conflict_cells, removal_conflict_cells.end()),
3139 could_lock_zone);
3140 if (insertion_conflict_cells.empty())
3141 return old_vertex;//new_position coincides with an existing vertex (not old_vertex)
3142 //and old_vertex should not be removed of the nb_vertices will change
3143 reset_circumcenter_cache(removal_conflict_cells);
3144 reset_sliver_cache(removal_conflict_cells);
3145 reset_circumcenter_cache(insertion_conflict_cells);
3146 reset_sliver_cache(insertion_conflict_cells);
3147
3148 if (could_lock_zone && *could_lock_zone == false)
3149 return Vertex_handle();
3150
3151 lock_outdated_cells();
3152 for(typename Cell_set::iterator it = insertion_conflict_cells.begin();
3153 it != insertion_conflict_cells.end(); ++it)
3154 outdated_cells_set.erase(*it);
3155 for(typename Cell_set::iterator it = removal_conflict_cells.begin();
3156 it != removal_conflict_cells.end(); ++it)
3157 outdated_cells_set.erase(*it);
3158 unlock_outdated_cells();
3159
3160 Cell_vector outdated_cells;
3161 Vertex_handle nv = move_point_topo_change_conflict_zone_known(old_vertex, new_position,
3162 insertion_conflict_boundary[0],
3163 insertion_conflict_cells.begin(),
3164 insertion_conflict_cells.end(),
3165 removal_conflict_cells.begin(),
3166 removal_conflict_cells.end(),
3167 std::back_inserter(outdated_cells),
3168 CGAL::Emptyset_iterator()); // deleted_cells
3169
3170 lock_outdated_cells();
3171 for(typename Cell_vector::iterator it = outdated_cells.begin();
3172 it != outdated_cells.end(); ++it)
3173 outdated_cells_set.insert(*it);
3174 unlock_outdated_cells();
3175
3176 return nv;
3177 }
3178
3179 template <typename C3T3, typename MD>
3180 template <typename OutdatedCellsOutputIterator,
3181 typename DeletedCellsOutputIterator>
3182 typename C3T3_helpers<C3T3,MD>::Vertex_handle
3183 C3T3_helpers<C3T3,MD>::
move_point_topo_change(const Vertex_handle & old_vertex,const Weighted_point & new_position,OutdatedCellsOutputIterator outdated_cells,DeletedCellsOutputIterator deleted_cells)3184 move_point_topo_change(const Vertex_handle& old_vertex,
3185 const Weighted_point& new_position,
3186 OutdatedCellsOutputIterator outdated_cells,
3187 DeletedCellsOutputIterator deleted_cells) const
3188 {
3189 #ifdef CGAL_MESH_3_C3T3_HELPERS_VERBOSE
3190 std::cerr << "move_point_topo_change with delcells(" << (void*)(&*old_vertex)
3191 << "=" << tr_.point(old_vertex)
3192 << " " << new_position << ",\n"
3193 << ")" << std::endl;
3194 #endif
3195
3196 Cell_set insertion_conflict_cells;
3197 Cell_set removal_conflict_cells;
3198 Facet_vector insertion_conflict_boundary;
3199 insertion_conflict_boundary.reserve(64);
3200
3201 get_conflict_zone_topo_change(old_vertex, new_position,
3202 std::inserter(insertion_conflict_cells,insertion_conflict_cells.end()),
3203 std::back_inserter(insertion_conflict_boundary),
3204 std::inserter(removal_conflict_cells, removal_conflict_cells.end()));
3205 reset_circumcenter_cache(removal_conflict_cells);
3206 reset_sliver_cache(removal_conflict_cells);
3207 reset_circumcenter_cache(insertion_conflict_cells);
3208 reset_sliver_cache(insertion_conflict_cells);
3209
3210 Vertex_handle nv = move_point_topo_change_conflict_zone_known(old_vertex, new_position,
3211 insertion_conflict_boundary[0],
3212 insertion_conflict_cells.begin(),
3213 insertion_conflict_cells.end(),
3214 removal_conflict_cells.begin(),
3215 removal_conflict_cells.end(),
3216 outdated_cells,
3217 deleted_cells);
3218
3219 return nv;
3220 }
3221
3222
3223 template <typename C3T3, typename MD>
3224 template < typename ConflictCellsInputIterator,
3225 typename OutdatedCellsOutputIterator,
3226 typename DeletedCellsOutputIterator >
3227 typename C3T3_helpers<C3T3,MD>::Vertex_handle
3228 C3T3_helpers<C3T3,MD>::
move_point_topo_change_conflict_zone_known(const Vertex_handle & old_vertex,const Weighted_point & new_position,const Facet & insertion_boundary_facet,ConflictCellsInputIterator insertion_conflict_cells_begin,ConflictCellsInputIterator insertion_conflict_cells_end,ConflictCellsInputIterator removal_conflict_cells_begin,ConflictCellsInputIterator removal_conflict_cells_end,OutdatedCellsOutputIterator outdated_cells,DeletedCellsOutputIterator deleted_cells)3229 move_point_topo_change_conflict_zone_known(
3230 const Vertex_handle& old_vertex,
3231 const Weighted_point& new_position,
3232 const Facet& insertion_boundary_facet,
3233 ConflictCellsInputIterator insertion_conflict_cells_begin,//ordered
3234 ConflictCellsInputIterator insertion_conflict_cells_end,
3235 ConflictCellsInputIterator removal_conflict_cells_begin,//ordered
3236 ConflictCellsInputIterator removal_conflict_cells_end,
3237 OutdatedCellsOutputIterator outdated_cells,
3238 DeletedCellsOutputIterator deleted_cells)//warning : this should not be an iterator to Intrusive_list
3239 //o.w. deleted_cells will point to null pointer or so and crash
3240 const
3241 {
3242 Weighted_point old_position = tr_.point(old_vertex); // intentional copy
3243 // make one set with conflict zone
3244 Cell_set conflict_zone;
3245 std::set_union(insertion_conflict_cells_begin, insertion_conflict_cells_end,
3246 removal_conflict_cells_begin, removal_conflict_cells_end,
3247 std::inserter(conflict_zone, conflict_zone.end()));
3248
3249 // Remove conflict zone cells from c3t3 (they will be deleted by insert/remove)
3250 remove_cells_and_facets_from_c3t3(conflict_zone.begin(), conflict_zone.end());
3251
3252 // Start Move point // Insert new_vertex, remove old_vertex
3253 int dimension = c3t3_.in_dimension(old_vertex);
3254 Index vertex_index = c3t3_.index(old_vertex);
3255 FT meshing_info = old_vertex->meshing_info();
3256
3257 // insert new point
3258 Vertex_handle new_vertex = tr_.insert_in_hole(new_position,
3259 insertion_conflict_cells_begin,
3260 insertion_conflict_cells_end,
3261 insertion_boundary_facet.first,
3262 insertion_boundary_facet.second);
3263
3264 // If new_position is hidden, update what should be and return default constructed handle
3265 if ( Vertex_handle() == new_vertex )
3266 {
3267 std::copy(conflict_zone.begin(), conflict_zone.end(), outdated_cells);
3268 return old_vertex;
3269 }
3270
3271 // remove old point
3272 tr_.remove(old_vertex);
3273
3274 c3t3_.set_dimension(new_vertex,dimension);
3275 c3t3_.set_index(new_vertex,vertex_index);
3276 new_vertex->set_meshing_info(meshing_info);
3277 // End Move point
3278
3279 //// Fill outdated_cells
3280 // Get conflict zone in new triangulation and set cells outdated
3281 Cell_vector new_conflict_cells;
3282 new_conflict_cells.reserve(64);
3283 get_conflict_zone_topo_change(new_vertex, old_position,
3284 std::back_inserter(new_conflict_cells));
3285 std::copy(new_conflict_cells.begin(),new_conflict_cells.end(),outdated_cells);
3286
3287 // Fill deleted_cells
3288 if(! boost::is_same<DeletedCellsOutputIterator,CGAL::Emptyset_iterator>::value)
3289 std::copy(conflict_zone.begin(), conflict_zone.end(), deleted_cells);
3290
3291 return new_vertex;
3292 }
3293
3294
3295 template <typename C3T3, typename MD>
3296 typename C3T3_helpers<C3T3,MD>::Vertex_handle
3297 C3T3_helpers<C3T3,MD>::
move_point_topo_change(const Vertex_handle & old_vertex,const Weighted_point & new_position)3298 move_point_topo_change(const Vertex_handle& old_vertex,
3299 const Weighted_point& new_position) const
3300 {
3301 // Insert new_vertex, remove old_vertex
3302 int dimension = c3t3_.in_dimension(old_vertex);
3303 Index vertex_index = c3t3_.index(old_vertex);
3304 FT meshing_info = old_vertex->meshing_info();
3305
3306 // insert new point
3307 Vertex_handle new_vertex = tr_.insert(new_position, old_vertex->cell());
3308 // If new_position is hidden, return default constructed handle
3309 if ( Vertex_handle() == new_vertex ) { return Vertex_handle(); }
3310 // remove old point
3311 tr_.remove(old_vertex);
3312
3313 c3t3_.set_dimension(new_vertex,dimension);
3314 c3t3_.set_index(new_vertex,vertex_index);
3315 new_vertex->set_meshing_info(meshing_info);
3316
3317 return new_vertex;
3318 }
3319
3320
3321 template <typename C3T3, typename MD>
3322 template <typename OutdatedCellsOutputIterator>
3323 typename C3T3_helpers<C3T3,MD>::Vertex_handle
3324 C3T3_helpers<C3T3,MD>::
move_point_no_topo_change(const Vertex_handle & old_vertex,const Vector_3 & move,const Weighted_point & new_position,OutdatedCellsOutputIterator outdated_cells)3325 move_point_no_topo_change(const Vertex_handle& old_vertex,
3326 const Vector_3& move,
3327 const Weighted_point& new_position,
3328 OutdatedCellsOutputIterator outdated_cells) const
3329 {
3330
3331 lock_outdated_cells();
3332 get_conflict_zone_no_topo_change(old_vertex, outdated_cells);
3333 unlock_outdated_cells();
3334
3335 return move_point_no_topo_change(old_vertex, move, new_position);
3336 }
3337
3338
3339 template <typename C3T3, typename MD>
3340 typename C3T3_helpers<C3T3,MD>::Vertex_handle
3341 C3T3_helpers<C3T3,MD>::
move_point_no_topo_change(const Vertex_handle & old_vertex,const Vector_3 & move,const Weighted_point & new_position)3342 move_point_no_topo_change(const Vertex_handle& old_vertex,
3343 const Vector_3& move,
3344 const Weighted_point& new_position) const
3345 {
3346 // Change vertex position
3347 tr_.set_point(old_vertex, move, new_position);
3348 CGAL_expensive_postcondition(tr_.is_valid());
3349 return old_vertex;
3350 }
3351
3352
3353 /**
3354 * @brief Returns the projection of \c p, using direction of
3355 * \c projection_vector
3356 */
3357 template <typename C3T3, typename MD>
3358 typename C3T3_helpers<C3T3,MD>::Bare_point
3359 C3T3_helpers<C3T3,MD>::
project_on_surface_aux(const Bare_point & p,const Bare_point & ref_point,const Vector_3 & projection_vector)3360 project_on_surface_aux(const Bare_point& p,
3361 const Bare_point& ref_point,
3362 const Vector_3& projection_vector) const
3363 {
3364 typedef typename Gt::Segment_3 Segment_3;
3365
3366 // Build a segment directed as projection_direction,
3367 typename Gt::Compute_squared_distance_3 sq_distance = tr_.geom_traits().compute_squared_distance_3_object();
3368 typename Gt::Compute_squared_length_3 sq_length = tr_.geom_traits().compute_squared_length_3_object();
3369 typename Gt::Construct_scaled_vector_3 scale = tr_.geom_traits().construct_scaled_vector_3_object();
3370 typename Gt::Is_degenerate_3 is_degenerate = tr_.geom_traits().is_degenerate_3_object();
3371 typename Gt::Construct_translated_point_3 translate = tr_.geom_traits().construct_translated_point_3_object();
3372
3373 typename MD::Construct_intersection construct_intersection =
3374 domain_.construct_intersection_object();
3375
3376 const FT sq_dist = sq_distance(p, ref_point);
3377 const FT sq_proj_length = sq_length(projection_vector);
3378
3379 if ( CGAL_NTS is_zero(sq_proj_length) )
3380 return ref_point;
3381
3382 const Vector_3 projection_scaled_vector =
3383 scale(projection_vector, CGAL::sqrt(sq_dist / sq_proj_length));
3384
3385 const Bare_point source = translate(p, projection_scaled_vector);
3386 const Bare_point target = translate(p, - projection_scaled_vector);
3387
3388 const Segment_3 proj_segment(source, target);
3389
3390 if ( is_degenerate(proj_segment) )
3391 return ref_point;
3392
3393 #ifndef CGAL_MESH_3_NO_LONGER_CALLS_DO_INTERSECT_3
3394
3395 typename MD::Do_intersect_surface do_intersect =
3396 domain_.do_intersect_surface_object();
3397
3398 if ( do_intersect(proj_segment) )
3399 return std::get<0>(construct_intersection(proj_segment));
3400 else
3401 return ref_point;
3402
3403 #else // CGAL_MESH_3_NO_LONGER_CALLS_DO_INTERSECT_3
3404
3405 typedef typename MD::Intersection Intersection;
3406 Intersection intersection = construct_intersection(proj_segment);
3407 if(std::get<2>(intersection) == 2)
3408 return std::get<0>(intersection);
3409 else
3410 return ref_point;
3411
3412 #endif // CGAL_MESH_3_NO_LONGER_CALLS_DO_INTERSECT_3
3413 }
3414
3415
3416 template <typename C3T3, typename MD>
3417 boost::optional<typename C3T3_helpers<C3T3,MD>::Plane_3>
3418 C3T3_helpers<C3T3,MD>::
get_least_square_surface_plane(const Vertex_handle & v,Bare_point & reference_point,Surface_patch_index patch_index)3419 get_least_square_surface_plane(const Vertex_handle& v,
3420 Bare_point& reference_point,
3421 Surface_patch_index patch_index) const
3422 {
3423 typedef typename C3T3::Triangulation::Triangle Triangle;
3424 typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object();
3425
3426 // Get incident facets
3427 Facet_vector facets;
3428 # ifdef CGAL_LINKED_WITH_TBB
3429 // Parallel
3430 if (boost::is_convertible<Concurrency_tag, Parallel_tag>::value)
3431 {
3432 tr_.finite_incident_facets_threadsafe(v, std::back_inserter(facets));
3433 }
3434 // Sequential
3435 else
3436 # endif // CGAL_LINKED_WITH_TBB
3437 {
3438 tr_.finite_incident_facets(v, std::back_inserter(facets));
3439 }
3440
3441 const Weighted_point& position = tr_.point(v);
3442
3443 // Get adjacent surface points
3444 std::vector<Triangle> triangles;
3445 typename C3T3::Facet ref_facet;
3446
3447 for (typename C3T3::Facet f : facets)
3448 {
3449 if ( c3t3_.is_in_complex(f) &&
3450 (patch_index == Surface_patch_index() ||
3451 c3t3_.surface_patch_index(f) == patch_index) )
3452 {
3453 ref_facet = f;
3454
3455 // In the case of a periodic triangulation, the incident facets of a point
3456 // do not necessarily have the same offsets. Worse, the surface centers
3457 // might not have the same offset as their facet. Thus, no solution except
3458 // calling a function 'get_closest_triangle(p, t)' that simply returns t
3459 // for a non-periodic triangulation, and checks all possible offsets for
3460 // periodic triangulations
3461
3462 Triangle t = c3t3_.triangulation().triangle(f);
3463 Triangle ct = tr_.get_closest_triangle(cp(position), t);
3464 triangles.push_back(ct);
3465 }
3466 }
3467
3468 // In some cases point is not a real surface point
3469 if ( triangles.empty() )
3470 return boost::none;
3471
3472 // Compute least square fitting plane
3473 Plane_3 plane;
3474 Bare_point point;
3475
3476 CGAL::linear_least_squares_fitting_3(triangles.begin(),
3477 triangles.end(),
3478 plane,
3479 point,
3480 Dimension_tag<2>(),
3481 tr_.geom_traits(),
3482 Default_diagonalize_traits<FT, 3>());
3483
3484 reference_point = ref_facet.first->get_facet_surface_center(ref_facet.second);
3485
3486 return plane;
3487 }
3488
3489
3490
3491 template <typename C3T3, typename MD>
3492 typename C3T3_helpers<C3T3,MD>::Bare_point
3493 C3T3_helpers<C3T3,MD>::
project_on_surface(const Vertex_handle & v,const Bare_point & p,Surface_patch_index index)3494 project_on_surface(const Vertex_handle& v,
3495 const Bare_point& p,
3496 Surface_patch_index index) const
3497 {
3498 boost::optional<Bare_point> opt_point =
3499 project_on_surface_if_possible(v, p, index);
3500 if(opt_point) return *opt_point;
3501 else return p;
3502 }
3503
3504
3505 template <typename C3T3, typename MD>
3506 boost::optional<typename C3T3_helpers<C3T3,MD>::Bare_point>
3507 C3T3_helpers<C3T3,MD>::
project_on_surface_if_possible(const Vertex_handle & v,const Bare_point & p,Surface_patch_index index)3508 project_on_surface_if_possible(const Vertex_handle& v,
3509 const Bare_point& p,
3510 Surface_patch_index index) const
3511 {
3512 // return domain_.project_on_surface(p);
3513
3514 typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object();
3515 typename Gt::Equal_3 equal = tr_.geom_traits().equal_3_object();
3516
3517 // Get plane
3518 Bare_point reference_point;
3519 boost::optional<Plane_3> opt_plane = get_least_square_surface_plane(v, reference_point, index);
3520
3521 if(!opt_plane) return boost::none;
3522
3523 // Project
3524 const Weighted_point& position = tr_.point(v);
3525 if ( ! equal(p, cp(position)) )
3526 return project_on_surface_aux(p, cp(position), opt_plane->orthogonal_vector());
3527 else
3528 return project_on_surface_aux(p, reference_point, opt_plane->orthogonal_vector());
3529 }
3530
3531
3532
3533 template <typename C3T3, typename MD>
3534 template <typename SliverCriterion>
3535 typename C3T3_helpers<C3T3,MD>::FT
3536 C3T3_helpers<C3T3,MD>::
min_incident_value(const Vertex_handle & vh,const SliverCriterion & criterion)3537 min_incident_value(const Vertex_handle& vh,
3538 const SliverCriterion& criterion) const
3539 {
3540 Cell_vector incident_cells_;
3541 tr_.finite_incident_cells(vh,std::back_inserter(incident_cells_));
3542
3543 return min_sliver_in_c3t3_value(incident_cells_, criterion);
3544 }
3545
3546 template <typename OutputIterator, typename CH, typename Fct>
3547 struct Filter {
3548
3549 mutable OutputIterator out;
3550 const Fct& fct;
3551
FilterFilter3552 Filter(OutputIterator out, const Fct& fct)
3553 : out(out), fct(fct)
3554 {}
3555
operatorFilter3556 void operator()(CH cell_handle) const
3557 {
3558 if(fct(cell_handle)){
3559 *out++ = cell_handle;
3560 }
3561 }
3562
3563 };
3564
3565 template <typename CH, typename Fct>
3566 struct Counter {
3567
3568 const Fct& fct;
3569 std::size_t& count;
3570
CounterCounter3571 Counter(const Fct& fct, std::size_t& count)
3572 : fct(fct), count(count)
3573 {}
3574
operatorCounter3575 void operator()(CH cell_handle)
3576 {
3577 if(fct(cell_handle)){
3578 ++count;
3579 }
3580 }
3581
3582 };
3583
3584 template <typename C3T3, typename MD>
3585 template <typename SliverCriterion>
3586 void
3587 C3T3_helpers<C3T3,MD>::
get_incident_slivers_without_using_tds_data(const Vertex_handle & v,const SliverCriterion & criterion,const FT & sliver_bound,Cell_vector & slivers)3588 get_incident_slivers_without_using_tds_data(const Vertex_handle& v,
3589 const SliverCriterion& criterion,
3590 const FT& sliver_bound,
3591 Cell_vector &slivers) const
3592 {
3593 typedef SliverCriterion Sc;
3594 typedef std::back_insert_iterator<Cell_vector> OutputIt;
3595 typedef Filter<OutputIt, Cell_handle, Is_sliver<Sc> > F;
3596 OutputIt slivers_it = std::back_inserter(slivers);
3597 Is_sliver<Sc> i_s(c3t3_, criterion, sliver_bound);
3598 F f(slivers_it, i_s);
3599 tr_.incident_cells_threadsafe(v, boost::make_function_output_iterator(f));
3600 }
3601
3602 template <typename C3T3, typename MD>
3603 template <typename Filter>
3604 bool
3605 C3T3_helpers<C3T3,MD>::
try_lock_and_get_incident_cells(const Vertex_handle & v,Cell_vector & cells,const Filter & filter)3606 try_lock_and_get_incident_cells(const Vertex_handle& v,
3607 Cell_vector &cells,
3608 const Filter &filter) const
3609 {
3610 std::vector<Cell_handle> tmp_cells;
3611 tmp_cells.reserve(64);
3612 bool ret = tr_.try_lock_and_get_incident_cells(v, tmp_cells);
3613 if (ret)
3614 {
3615 for(Cell_handle ch : tmp_cells)
3616 {
3617 if (filter(ch))
3618 cells.push_back(ch);
3619 }
3620 }
3621 else
3622 tr_.unlock_all_elements();
3623 return ret;
3624 }
3625
3626 template <typename C3T3, typename MD>
3627 template <typename SliverCriterion>
3628 bool
3629 C3T3_helpers<C3T3,MD>::
try_lock_and_get_incident_slivers(const Vertex_handle & v,const SliverCriterion & criterion,const FT & sliver_bound,Cell_vector & slivers)3630 try_lock_and_get_incident_slivers(const Vertex_handle& v,
3631 const SliverCriterion& criterion,
3632 const FT& sliver_bound,
3633 Cell_vector &slivers) const
3634 {
3635 Is_sliver<SliverCriterion> i_s(c3t3_, criterion, sliver_bound);
3636 return try_lock_and_get_incident_cells(v, slivers, i_s);
3637 }
3638
3639
3640 template <typename C3T3, typename MD>
3641 template <typename SliverCriterion, typename OutputIterator>
3642 OutputIterator
3643 C3T3_helpers<C3T3,MD>::
incident_slivers(const Vertex_handle & v,const SliverCriterion & criterion,const FT & sliver_bound,OutputIterator out)3644 incident_slivers(const Vertex_handle& v,
3645 const SliverCriterion& criterion,
3646 const FT& sliver_bound,
3647 OutputIterator out) const
3648 {
3649 typedef SliverCriterion Sc;
3650
3651 std::vector<Cell_handle> incident_cells_;
3652 tr_.incident_cells(v, std::back_inserter(incident_cells_));
3653
3654 #ifdef CGAL_CXX17
3655 std::remove_copy_if(incident_cells_.begin(),
3656 incident_cells_.end(),
3657 out,
3658 std::not_fn(Is_sliver<Sc>(c3t3_, criterion, sliver_bound)));
3659 #else
3660 std::remove_copy_if(incident_cells_.begin(),
3661 incident_cells_.end(),
3662 out,
3663 std::not1(Is_sliver<Sc>(c3t3_,criterion,sliver_bound)));
3664 #endif
3665
3666 return out;
3667 }
3668
3669 template <typename C3T3, typename MD>
3670 template <typename SliverCriterion, typename OutputIterator>
3671 OutputIterator
3672 C3T3_helpers<C3T3,MD>::
new_incident_slivers(const Vertex_handle & v,const SliverCriterion & criterion,const FT & sliver_bound,OutputIterator out)3673 new_incident_slivers(const Vertex_handle& v,
3674 const SliverCriterion& criterion,
3675 const FT& sliver_bound,
3676 OutputIterator out) const
3677 {
3678 typedef SliverCriterion Sc;
3679 typedef Filter<OutputIterator,Cell_handle,Is_sliver<Sc> > F;
3680
3681 Is_sliver<Sc> i_s(c3t3_, criterion, sliver_bound);
3682 F f(out, i_s);
3683 tr_.incident_cells(v,boost::make_function_output_iterator(f));
3684
3685 return f.out;
3686 }
3687
3688 template <typename C3T3, typename MD>
3689 template <typename SliverCriterion>
3690 bool
3691 C3T3_helpers<C3T3,MD>::
is_sliver(const Cell_handle & ch,const SliverCriterion & criterion,const FT & sliver_bound)3692 is_sliver(const Cell_handle& ch,
3693 const SliverCriterion& criterion,
3694 const FT& sliver_bound) const
3695 {
3696 Is_sliver<SliverCriterion> iss(c3t3_,criterion,sliver_bound);
3697 return iss(ch);
3698 }
3699
3700 template <typename C3T3, typename MD>
3701 template <typename SliverCriterion>
3702 std::size_t
3703 C3T3_helpers<C3T3,MD>::
number_of_incident_slivers(const Vertex_handle & v,const SliverCriterion & criterion,const FT & sliver_bound)3704 number_of_incident_slivers(const Vertex_handle& v,
3705 const SliverCriterion& criterion,
3706 const FT& sliver_bound) const
3707 {
3708 typedef SliverCriterion Sc;
3709 typedef Counter<Cell_handle,Is_sliver<Sc> > C;
3710
3711 std::size_t count = 0;
3712 Is_sliver<Sc> is_sliver(c3t3_,criterion,sliver_bound);
3713 C c(is_sliver, count);
3714 tr_.incident_cells(v, boost::make_function_output_iterator(c));
3715
3716 return count;
3717 }
3718
3719
3720 template <typename C3T3, typename MD>
3721 template <typename SliverCriterion>
3722 typename C3T3_helpers<C3T3,MD>::FT
3723 C3T3_helpers<C3T3,MD>::
min_sliver_value(const Cell_vector & cells,const SliverCriterion & criterion,const bool use_cache)3724 min_sliver_value(const Cell_vector& cells,
3725 const SliverCriterion& criterion,
3726 const bool use_cache) const
3727 {
3728 using boost::make_transform_iterator;
3729
3730 if ( cells.empty() )
3731 return criterion.get_max_value();
3732
3733 if ( ! use_cache )
3734 {
3735 reset_sliver_cache(cells.begin(),cells.end());
3736 }
3737
3738 // Return min dihedral angle
3739 //Sliver_criterion_value<SliverCriterion> sc_value(tr_,criterion);
3740 //
3741 //return *(std::min_element(make_transform_iterator(cells.begin(),sc_value),
3742 // make_transform_iterator(cells.end(),sc_value)));
3743 FT min_value = criterion.get_max_value();
3744 for(typename Cell_vector::const_iterator it = cells.begin();
3745 it != cells.end();
3746 ++it)
3747 {
3748 min_value = (std::min)(criterion(*it), min_value);
3749 }
3750 return min_value;
3751 }
3752
3753
3754 template <typename C3T3, typename MD>
3755 template <typename InputIterator, typename OutputIterator>
3756 void
3757 C3T3_helpers<C3T3,MD>::
fill_modified_vertices(InputIterator cells_begin,InputIterator cells_end,const Vertex_handle & vertex,OutputIterator out)3758 fill_modified_vertices(InputIterator cells_begin,
3759 InputIterator cells_end,
3760 const Vertex_handle& vertex,
3761 OutputIterator out) const
3762 {
3763 Vertex_set already_inserted_vertices;
3764 // Dont insert vertex in out
3765 already_inserted_vertices.insert(vertex);
3766
3767 for ( InputIterator it = cells_begin ; it != cells_end ; ++it )
3768 {
3769 for ( int i=0 ; i<4 ; ++i )
3770 {
3771 // Insert vertices if not already inserted
3772 const Vertex_handle& current_vertex = (*it)->vertex(i);
3773 if ( !tr_.is_infinite(current_vertex)
3774 && already_inserted_vertices.insert(current_vertex).second )
3775 {
3776 *out++ = current_vertex;
3777 }
3778 }
3779 }
3780 }
3781
3782
3783 template <typename C3T3, typename MD>
3784 template <typename CellsVector, typename CellDataSet>
3785 void
3786 C3T3_helpers<C3T3,MD>::
fill_cells_backup(const CellsVector & cells,CellDataSet & cells_backup)3787 fill_cells_backup(const CellsVector& cells,
3788 CellDataSet& cells_backup) const
3789 {
3790 typedef typename CellDataSet::value_type Cell_data;
3791 typename CellsVector::const_iterator cit;
3792 for(cit = cells.begin(); cit != cells.end(); ++cit)
3793 {
3794 cells_backup.insert(Cell_data(c3t3_,*cit));
3795 }
3796 }
3797
3798 template <typename C3T3, typename MD>
3799 template <typename CellsVector, typename CellDataSet>
3800 void
3801 C3T3_helpers<C3T3,MD>::
restore_from_cells_backup(const CellsVector & cells,CellDataSet & cells_backup)3802 restore_from_cells_backup(const CellsVector& cells,
3803 CellDataSet& cells_backup) const
3804 {
3805 for(typename CellsVector::const_iterator cit = cells.begin();
3806 cit != cells.end();
3807 ++cit)
3808 {
3809 typename CellDataSet::const_iterator cd_it
3810 = cells_backup.find(Cell_data_backup(c3t3_, *cit, false/*don't backup*/));
3811 if(cd_it != cells_backup.end())
3812 {
3813 typename CellDataSet::value_type cell_data = *cd_it;
3814 cell_data.restore(*cit, c3t3_);
3815 cells_backup.erase(cd_it);
3816 }
3817 else CGAL_error();
3818 }
3819 CGAL_assertion(cells_backup.empty());
3820 }
3821
3822 template <typename C3T3, typename MD>
3823 template <typename OutputIterator>
3824 OutputIterator
3825 C3T3_helpers<C3T3,MD>::
get_conflict_zone_no_topo_change(const Vertex_handle & vertex,OutputIterator conflict_cells)3826 get_conflict_zone_no_topo_change(const Vertex_handle& vertex,
3827 OutputIterator conflict_cells) const
3828 {
3829 return tr_.incident_cells(vertex, conflict_cells);
3830 }
3831
3832 template <typename C3T3, typename MD>
3833 template <typename CellsOutputIterator,
3834 typename FacetsOutputIterator>
3835 void
3836 C3T3_helpers<C3T3,MD>::
get_conflict_zone_topo_change(const Vertex_handle & v,const Weighted_point & conflict_point,CellsOutputIterator insertion_conflict_cells,FacetsOutputIterator insertion_conflict_boundary,CellsOutputIterator removal_conflict_cells,bool * could_lock_zone)3837 get_conflict_zone_topo_change(const Vertex_handle& v,
3838 const Weighted_point& conflict_point,
3839 CellsOutputIterator insertion_conflict_cells,
3840 FacetsOutputIterator insertion_conflict_boundary,
3841 CellsOutputIterator removal_conflict_cells,
3842 bool *could_lock_zone) const
3843 {
3844 // Get triangulation_vertex incident cells : removal conflict zone
3845 // TODO: hasn't it already been computed in "perturb_vertex" (when getting the slivers)?
3846 // We don't try to lock the incident cells since they've already been locked
3847
3848 # ifdef CGAL_LINKED_WITH_TBB
3849 // Parallel
3850 if (boost::is_convertible<Concurrency_tag, Parallel_tag>::value)
3851 {
3852 tr_.incident_cells_threadsafe(v, removal_conflict_cells);
3853 }
3854 // Sequential
3855 else
3856 # endif // CGAL_LINKED_WITH_TBB
3857 {
3858 tr_.incident_cells(v, removal_conflict_cells);
3859 }
3860
3861 // Get conflict_point conflict zone
3862 int li=0;
3863 int lj=0;
3864 typename Tr::Locate_type lt;
3865 Cell_handle cell = tr_.locate(
3866 conflict_point, lt, li, lj, v->cell(), could_lock_zone);
3867
3868 if (could_lock_zone && *could_lock_zone == false)
3869 return;
3870
3871 if ( lt == Tr::VERTEX ) // Vertex removal is forbidden
3872 return;
3873
3874 // Find conflict zone
3875 tr_.find_conflicts(conflict_point,
3876 cell,
3877 insertion_conflict_boundary,
3878 insertion_conflict_cells,
3879 could_lock_zone);
3880 }
3881
3882 template <typename C3T3, typename MD>
3883 template <typename OutputIterator>
3884 OutputIterator
3885 C3T3_helpers<C3T3,MD>::
get_conflict_zone_topo_change(const Vertex_handle & vertex,const Weighted_point & conflict_point,OutputIterator conflict_cells)3886 get_conflict_zone_topo_change(const Vertex_handle& vertex,
3887 const Weighted_point& conflict_point,
3888 OutputIterator conflict_cells) const
3889 {
3890 // Get triangulation_vertex incident cells
3891 Cell_vector incident_cells_;
3892 incident_cells_.reserve(64);
3893 tr_.incident_cells(vertex, std::back_inserter(incident_cells_));
3894
3895 // Get conflict_point conflict zone
3896 Cell_vector deleted_cells;
3897 deleted_cells.reserve(64);
3898
3899 // Vertex removal is forbidden
3900 int li=0;
3901 int lj=0;
3902 typename Tr::Locate_type locate_type;
3903 Cell_handle cell = tr_.locate(conflict_point, locate_type, li, lj, vertex->cell());
3904
3905 if ( Tr::VERTEX == locate_type )
3906 return conflict_cells;
3907
3908 // Find conflict zone
3909 tr_.find_conflicts(conflict_point,
3910 cell,
3911 CGAL::Emptyset_iterator(),
3912 std::back_inserter(deleted_cells),
3913 CGAL::Emptyset_iterator());
3914
3915 // Compute union of conflict_point conflict zone and triangulation_vertex
3916 // incident cells
3917 std::sort(deleted_cells.begin(),deleted_cells.end());
3918 std::sort(incident_cells_.begin(),incident_cells_.end());
3919
3920 std::set_union(deleted_cells.begin(), deleted_cells.end(),
3921 incident_cells_.begin(), incident_cells_.end(),
3922 conflict_cells);
3923
3924 return conflict_cells;
3925 }
3926
3927
3928 template <typename C3T3, typename MD>
3929 typename C3T3_helpers<C3T3,MD>::Facet_boundary
3930 C3T3_helpers<C3T3,MD>::
get_surface_boundary(const Vertex_handle & moving_vertex,const Facet_vector & facets,Vertex_set & incident_surface_vertices)3931 get_surface_boundary(const Vertex_handle& moving_vertex,
3932 const Facet_vector& facets,
3933 Vertex_set& incident_surface_vertices) const
3934 {
3935 Facet_boundary boundary;
3936 typename Facet_vector::const_iterator fit = facets.begin();
3937 for ( ; fit != facets.end() ; ++fit )
3938 {
3939 if ( c3t3_.is_in_complex(*fit) )
3940 {
3941 const Surface_patch_index surface_index = c3t3_.surface_patch_index(*fit);
3942 const int k = fit->second;
3943 Vertex_handle v1 = fit->first->vertex((k+1)&3);
3944 Vertex_handle v2 = fit->first->vertex((k+2)&3);
3945 Vertex_handle v3 = fit->first->vertex((k+3)&3);
3946 incident_surface_vertices.insert(v1);
3947 incident_surface_vertices.insert(v2);
3948 incident_surface_vertices.insert(v3);
3949
3950 // Check that each vertex is a surface one
3951 // This is a trick to ensure that in_domain vertices stay inside
3952 if ( c3t3_.in_dimension(v1) > 2
3953 || c3t3_.in_dimension(v2) > 2
3954 || c3t3_.in_dimension(v3) > 2 )
3955 {
3956 boundary.clear();
3957 return boundary; // return an empty boundary, that cannot be equal
3958 // to a real boundary
3959 }
3960
3961 order_handles(v1,v2,v3);
3962
3963 CGAL_assertion(v1<v2);
3964 CGAL_assertion(v2<v3);
3965
3966 update_boundary(boundary, Ordered_edge(v1,v2), v3, surface_index);
3967 update_boundary(boundary, Ordered_edge(v1,v3), v2, surface_index);
3968 update_boundary(boundary, Ordered_edge(v2,v3), v1, surface_index);
3969 }
3970
3971 incident_surface_vertices.erase(moving_vertex);
3972 }
3973
3974 return boundary;
3975 }
3976
3977 template <typename C3T3, typename MD>
3978 bool
3979 C3T3_helpers<C3T3,MD>::
check_no_inside_vertices(const Facet_vector & facets)3980 check_no_inside_vertices(const Facet_vector& facets) const
3981 {
3982 typename Facet_vector::const_iterator fit = facets.begin();
3983 for ( ; fit != facets.end() ; ++fit )
3984 {
3985 if ( c3t3_.is_in_complex(*fit) )
3986 {
3987 const int k = fit->second;
3988 const Vertex_handle& v1 = fit->first->vertex((k+1)&3);
3989 const Vertex_handle& v2 = fit->first->vertex((k+2)&3);
3990 const Vertex_handle& v3 = fit->first->vertex((k+3)&3);
3991
3992 // Check that each vertex is a surface one
3993 if ( c3t3_.in_dimension(v1) > 2
3994 || c3t3_.in_dimension(v2) > 2
3995 || c3t3_.in_dimension(v3) > 2 )
3996 {
3997 return false;
3998 }
3999 }
4000 }
4001
4002 return true;
4003 }
4004
4005 } // end namespace Mesh_3
4006 } // end namespace CGAL
4007
4008 #include <CGAL/enable_warnings.h>
4009
4010 #endif // CGAL_MESH_3_C3T3_HELPERS_H
4011