1 // Copyright (c) 2009-2014 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/Mesher_3.h $
7 // $Id: Mesher_3.h 59a0da4 2021-05-19T17:23:53+02:00 Laurent Rineau
8 // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
9 //
10 //
11 // Author(s) : Laurent Rineau, Stephane Tayeb, Clement Jamin
12 //
13 //******************************************************************************
14 // File Description :
15 // Implements a mesher_3 with two mesher levels : one for facets and one for
16 // tets.
17 //******************************************************************************
18
19 #ifndef CGAL_MESH_3_MESHER_3_H
20 #define CGAL_MESH_3_MESHER_3_H
21
22 #include <CGAL/license/Mesh_3.h>
23
24 #include <CGAL/disable_warnings.h>
25
26 #include <CGAL/Mesh_3/config.h>
27
28 #if CGAL_MESH_3_USE_INTEL_ITT
29 # include <ittnotify.h>
30 # define CGAL_MESH_3_TASK_BEGIN(task_handle) __itt_task_begin(mesh_3_domain, __itt_null, __itt_null, task_handle);
31 # define CGAL_MESH_3_TASK_END(task_handle) __itt_task_end(mesh_3_domain);
32 #else
33 # define CGAL_MESH_3_TASK_BEGIN(task_handle)
34 # define CGAL_MESH_3_TASK_END(task_handle)
35 #endif
36
37 #include <CGAL/Mesh_error_code.h>
38
39 #include <CGAL/Mesh_3/Dump_c3t3.h>
40
41 #include <CGAL/Mesh_3/Refine_facets_3.h>
42 #include <CGAL/Mesh_3/Refine_facets_manifold_base.h>
43 #include <CGAL/Mesh_3/Refine_cells_3.h>
44 #include <CGAL/Mesh_3/Refine_tets_visitor.h>
45 #include <CGAL/Mesher_level_visitors.h>
46 #include <CGAL/Kernel_traits.h>
47 #include <CGAL/point_generators_3.h>
48
49 #ifdef CGAL_MESH_3_USE_OLD_SURFACE_RESTRICTED_DELAUNAY_UPDATE
50 #include <CGAL/Surface_mesher/Surface_mesher_visitor.h>
51 #endif
52
53 #include <CGAL/Mesh_3/Concurrent_mesher_config.h>
54 #include <CGAL/Real_timer.h>
55
56 #ifdef CGAL_MESH_3_PROFILING
57 #include <CGAL/Mesh_3/Profiling_tools.h>
58 #endif
59
60 #ifdef CGAL_LINKED_WITH_TBB
61 # include <thread>
62 #endif
63
64 #include <boost/format.hpp>
65 #include <boost/type_traits/is_convertible.hpp>
66 #include <string>
67 #include <atomic>
68
69 namespace CGAL {
70 namespace Mesh_3 {
71
72 /************************************************
73 // Class Mesher_3_base
74 // Two versions: sequential / parallel
75 ************************************************/
76
77 // Sequential
78 template <typename Tr, typename Concurrency_tag>
79 class Mesher_3_base
80 {
81 protected:
82 typedef typename Tr::Lock_data_structure Lock_data_structure;
83
Mesher_3_base(const Bbox_3 &,int)84 Mesher_3_base(const Bbox_3 &, int) {}
85
get_lock_data_structure()86 Lock_data_structure *get_lock_data_structure() const { return 0; }
get_worksharing_data_structure()87 WorksharingDataStructureType *get_worksharing_data_structure() const { return 0; }
set_bbox(const Bbox_3 &)88 void set_bbox(const Bbox_3 &) {}
89 };
90
91 #ifdef CGAL_LINKED_WITH_TBB
92 // Parallel
93 template <typename Tr>
94 class Mesher_3_base<Tr, Parallel_tag>
95 {
96 protected:
97 typedef typename Tr::Lock_data_structure Lock_data_structure;
98
Mesher_3_base(const Bbox_3 & bbox,int num_grid_cells_per_axis)99 Mesher_3_base(const Bbox_3 &bbox, int num_grid_cells_per_axis)
100 : m_lock_ds(bbox, num_grid_cells_per_axis),
101 m_worksharing_ds(bbox)
102 {}
103
get_lock_data_structure()104 Lock_data_structure *get_lock_data_structure()
105 {
106 return &m_lock_ds;
107 }
get_worksharing_data_structure()108 WorksharingDataStructureType *get_worksharing_data_structure()
109 {
110 return &m_worksharing_ds;
111 }
get_worksharing_data_structure()112 const WorksharingDataStructureType *get_worksharing_data_structure() const
113 {
114 return &m_worksharing_ds;
115 }
116
set_bbox(const Bbox_3 & bbox)117 void set_bbox(const Bbox_3 &bbox)
118 {
119 m_lock_ds.set_bbox(bbox);
120 m_worksharing_ds.set_bbox(bbox);
121 }
122
123 /// Lock data structure
124 Lock_data_structure m_lock_ds;
125 /// Worksharing data structure
126 WorksharingDataStructureType m_worksharing_ds;
127 };
128 #endif // CGAL_LINKED_WITH_TBB
129
130
131 /************************************************
132 *
133 * Mesher_3 class
134 *
135 ************************************************/
136
137 template<class C3T3, class MeshCriteria, class MeshDomain>
138 class Mesher_3
139 : public Mesher_3_base<
140 typename C3T3::Triangulation,
141 #ifdef CGAL_DEBUG_FORCE_SEQUENTIAL_MESH_REFINEMENT
142 Sequential_tag
143 #else
144 typename C3T3::Concurrency_tag
145 #endif
146 >
147 {
148 public:
149 #ifdef CGAL_DEBUG_FORCE_SEQUENTIAL_MESH_REFINEMENT
150 typedef Sequential_tag Concurrency_tag;
151 #else
152 typedef typename C3T3::Concurrency_tag Concurrency_tag;
153 #endif
154 typedef typename C3T3::Triangulation Triangulation;
155 typedef typename Triangulation::Point Point;
156 typedef typename Triangulation::Bare_point Bare_point;
157 typedef typename Kernel_traits<Point>::Kernel Kernel;
158 typedef typename Kernel::Vector_3 Vector;
159 typedef typename MeshDomain::Index Index;
160
161 // Self
162 typedef Mesher_3<C3T3, MeshCriteria, MeshDomain> Self;
163 #ifdef CGAL_DEBUG_FORCE_SEQUENTIAL_MESH_REFINEMENT
164 typedef Mesher_3_base<Triangulation, Sequential_tag> Base;
165 #else
166 typedef Mesher_3_base<Triangulation, typename C3T3::Concurrency_tag> Base;
167 #endif
168
169 using Base::get_lock_data_structure;
170
171 //-------------------------------------------------------
172 // Mesher_levels
173 //-------------------------------------------------------
174 // typedef Refine_facets_manifold_base<Rf_base> Rf_manifold_base;
175
176 /// Facets mesher level
177 typedef Refine_facets_3<
178 Triangulation,
179 typename MeshCriteria::Facet_criteria,
180 MeshDomain,
181 C3T3,
182 Null_mesher_level,
183 Concurrency_tag,
184 Refine_facets_manifold_base
185 > Facets_level;
186
187 /// Cells mesher level
188 typedef Refine_cells_3<
189 Triangulation,
190 typename MeshCriteria::Cell_criteria,
191 MeshDomain,
192 C3T3,
193 Facets_level,
194 Concurrency_tag> Cells_level;
195
196 //-------------------------------------------------------
197 // Visitors
198 //-------------------------------------------------------
199 /// Facets visitor : to update cells when refining surface
200 typedef tets::Refine_facets_visitor<
201 Triangulation,
202 Cells_level,
203 Null_mesh_visitor> Facets_visitor;
204
205 #ifndef CGAL_MESH_3_USE_OLD_SURFACE_RESTRICTED_DELAUNAY_UPDATE
206 /// Cells visitor : it just need to know previous level
207 typedef Null_mesh_visitor_level<Facets_visitor> Cells_visitor;
208 #else
209 /// Cells visitor : to update surface (restore restricted Delaunay)
210 /// when refining cells
211 typedef Surface_mesher::Visitor<
212 Triangulation,
213 Facets_level,
214 Facets_visitor> Cells_visitor;
215 #endif
216
217 /// Constructor
218 Mesher_3(C3T3& c3t3,
219 const MeshDomain& domain,
220 const MeshCriteria& criteria,
221 int mesh_topology = 0,
222 std::size_t maximal_number_of_vertices = 0,
223 Mesh_error_code* error_code = 0
224 #ifndef CGAL_NO_ATOMIC
225 , std::atomic<bool>* stop_ptr = 0
226 #endif
227 );
228
229 /// Destructor
~Mesher_3()230 ~Mesher_3()
231 {
232 // The lock data structure is going to be destroyed
233 r_c3t3_.triangulation().set_lock_data_structure(nullptr);
234 }
235
236 /// Launch mesh refinement
237 double refine_mesh(std::string dump_after_refine_surface_prefix = "");
238
239 /// Debug
240 std::string debug_info() const;
241 std::string debug_info_header() const;
242
243 // Step-by-step methods
244 void initialize();
245 void fix_c3t3();
246 void display_number_of_bad_elements();
247 void one_step();
248 bool is_algorithm_done();
249
250 #ifdef CGAL_MESH_3_MESHER_STATUS_ACTIVATED
251 struct Mesher_status
252 {
253 std::size_t vertices, facet_queue, cells_queue;
254
Mesher_statusMesher_status255 Mesher_status(std::size_t v, std::size_t f, std::size_t c)
256 : vertices(v), facet_queue(f), cells_queue(c) {}
257 };
258
259 Mesher_status status() const;
260 #endif
261
262 private:
263 /// The oracle
264 const MeshDomain& r_oracle_;
265
266 /// Meshers
267 Null_mesher_level null_mesher_;
268 Facets_level facets_mesher_;
269 Cells_level cells_mesher_;
270
271 /// Visitors
272 Null_mesh_visitor null_visitor_;
273 Facets_visitor facets_visitor_;
274 Cells_visitor cells_visitor_;
275
276 /// The container of the resulting mesh
277 C3T3& r_c3t3_;
278
279 enum Refinement_stage {
280 NOT_INITIALIZED,
281 REFINE_FACETS,
282 REFINE_FACETS_AND_EDGES,
283 REFINE_FACETS_AND_EDGES_AND_VERTICES,
284 REFINE_ALL
285 };
286
287 Refinement_stage refinement_stage;
288
289 /// Maximal number of vertices
290 std::size_t maximal_number_of_vertices_;
291
292 /// Pointer to the error code
293 Mesh_error_code* const error_code_;
294
295 #ifndef CGAL_NO_ATOMIC
296 /// Pointer to the atomic Boolean that can stop the process
297 std::atomic<bool>* const stop_ptr;
298 #endif
299
300 #ifdef CGAL_LINKED_WITH_TBB
approximate_number_of_vertices(CGAL::Parallel_tag)301 std::size_t approximate_number_of_vertices(CGAL::Parallel_tag) const {
302 # if CGAL_CONCURRENT_COMPACT_CONTAINER_APPROXIMATE_SIZE
303 return r_c3t3_.triangulation().tds().vertices().approximate_size();
304 # else // not CGAL_CONCURRENT_COMPACT_CONTAINER_APPROXIMATE_SIZE
305 CGAL_error_msg(
306 "If you want to use the Mesh_3 feature \"maximal_number_of_vertices\"\n"
307 "with CGAL::Parallel_tag then you need to recompile the code with the\n"
308 "preprocessor macro CGAL_CONCURRENT_COMPACT_CONTAINER_APPROXIMATE_SIZE\n"
309 "set to 1. That will induce a performance loss of 3%.\n");
310 # endif // not CGAL_CONCURRENT_COMPACT_CONTAINER_APPROXIMATE_SIZE
311 }
312 #endif // CGAL_LINKED_WITH_TBB
313
approximate_number_of_vertices(CGAL::Sequential_tag)314 std::size_t approximate_number_of_vertices(CGAL::Sequential_tag) const {
315 return r_c3t3_.triangulation().number_of_vertices();
316 }
317
forced_stop()318 bool forced_stop() const {
319 #ifndef CGAL_NO_ATOMIC
320 if(stop_ptr != 0 &&
321 stop_ptr->load(std::memory_order_acquire) == true)
322 {
323 if(error_code_ != 0) *error_code_ = CGAL_MESH_3_STOPPED;
324 return true;
325 }
326 #endif // not defined CGAL_NO_ATOMIC
327 if(maximal_number_of_vertices_ != 0 &&
328 approximate_number_of_vertices(Concurrency_tag()) >= maximal_number_of_vertices_)
329 {
330 if(error_code_ != 0) {
331 *error_code_ = CGAL_MESH_3_MAXIMAL_NUMBER_OF_VERTICES_REACHED;
332 }
333 return true;
334 }
335 return false;
336 }
337
338 private:
339 // Disabled copy constructor
340 Mesher_3(const Self& src);
341 // Disabled assignment operator
342 Self& operator=(const Self& src);
343
344 }; // end class Mesher_3
345
346
347
348 template<class C3T3, class MC, class MD>
Mesher_3(C3T3 & c3t3,const MD & domain,const MC & criteria,int mesh_topology,std::size_t maximal_number_of_vertices,Mesh_error_code * error_code,std::atomic<bool> * stop_ptr)349 Mesher_3<C3T3,MC,MD>::Mesher_3(C3T3& c3t3,
350 const MD& domain,
351 const MC& criteria,
352 int mesh_topology,
353 std::size_t maximal_number_of_vertices,
354 Mesh_error_code* error_code
355 #ifndef CGAL_NO_ATOMIC
356 , std::atomic<bool>* stop_ptr
357 #endif
358 )
359 : Base(c3t3.bbox(),
360 Concurrent_mesher_config::get().locking_grid_num_cells_per_axis)
361 , r_oracle_(domain)
362 , null_mesher_()
363 , facets_mesher_(c3t3.triangulation(),
364 criteria.facet_criteria_object(),
365 domain,
366 null_mesher_,
367 c3t3,
368 mesh_topology,
369 maximal_number_of_vertices
370 #ifndef CGAL_NO_ATOMIC
371 , stop_ptr
372 #endif
373 )
374 , cells_mesher_(c3t3.triangulation(),
375 criteria.cell_criteria_object(),
376 domain,
377 facets_mesher_,
378 c3t3,
379 maximal_number_of_vertices
380 #ifndef CGAL_NO_ATOMIC
381 , stop_ptr
382 #endif
383 )
384 , null_visitor_()
385 , facets_visitor_(&cells_mesher_, &null_visitor_)
386 #ifndef CGAL_MESH_3_USE_OLD_SURFACE_RESTRICTED_DELAUNAY_UPDATE
387 , cells_visitor_(facets_visitor_)
388 #else
389 , cells_visitor_(&facets_mesher_, &facets_visitor_)
390 #endif
391 , r_c3t3_(c3t3)
392 , maximal_number_of_vertices_(maximal_number_of_vertices)
393 , error_code_(error_code)
394 #ifndef CGAL_NO_ATOMIC
395 , stop_ptr(stop_ptr)
396 #endif
397 {
398 facets_mesher_.set_lock_ds(this->get_lock_data_structure());
399 facets_mesher_.set_worksharing_ds(this->get_worksharing_data_structure());
400 cells_mesher_.set_lock_ds(this->get_lock_data_structure());
401 cells_mesher_.set_worksharing_ds(this->get_worksharing_data_structure());
402 #ifndef CGAL_NO_ATOMIC
403 cells_mesher_.set_stop_pointer(stop_ptr);
404 facets_mesher_.set_stop_pointer(stop_ptr);
405 #endif
406 }
407
408
409
410 template<class C3T3, class MC, class MD>
411 double
412 Mesher_3<C3T3,MC,MD>::
refine_mesh(std::string dump_after_refine_surface_prefix)413 refine_mesh(std::string dump_after_refine_surface_prefix)
414 {
415 CGAL::Real_timer timer;
416 timer.start();
417 double elapsed_time = 0.;
418 #if CGAL_MESH_3_USE_INTEL_ITT
419 auto mesh_3_domain = __itt_domain_create("org.cgal.Mesh_3.refine_mesh");
420 auto initialize_task_handle = __itt_string_handle_create("Mesher_3::initialize");
421 auto refine_surface_mesh_task_handle = __itt_string_handle_create("Mesher_3 refine surface mesh");
422 auto scan_cells_task_handle = __itt_string_handle_create("Mesher_3 scan triangulation for bad cells");
423 auto refine_volume_mesh_task_handle = __itt_string_handle_create("Mesher_3 refine volume mesh");
424 #endif // CGAL_MESH_3_USE_INTEL_ITT
425
426 // First surface mesh could modify c3t3 without notifying cells_mesher
427 // So we have to ensure that no old cell will be left in c3t3
428 // Second, the c3t3 object could have been corrupted since the last call
429 // to `refine_mesh`, for example by inserting new vertices in the
430 // triangulation.
431 r_c3t3_.clear_cells_and_facets_from_c3t3();
432
433 const Triangulation& r_tr = r_c3t3_.triangulation();
434 CGAL_USE(r_tr);
435
436 #ifndef CGAL_MESH_3_VERBOSE
437 // Scan surface and refine it
438 CGAL_MESH_3_TASK_BEGIN(initialize_task_handle);
439 initialize();
440 CGAL_MESH_3_TASK_END(initialize_task_handle);
441
442 #ifdef CGAL_MESH_3_PROFILING
443 std::cerr << "Refining facets..." << std::endl;
444 WallClockTimer t;
445 #endif
446 CGAL_MESH_3_TASK_BEGIN(refine_surface_mesh_task_handle);
447 facets_mesher_.refine(facets_visitor_);
448 facets_mesher_.scan_edges();
449 refinement_stage = REFINE_FACETS_AND_EDGES;
450 facets_mesher_.refine(facets_visitor_);
451 facets_mesher_.scan_vertices();
452 refinement_stage = REFINE_FACETS_AND_EDGES_AND_VERTICES;
453 facets_mesher_.refine(facets_visitor_);
454 CGAL_MESH_3_TASK_END(refine_surface_mesh_task_handle);
455 #ifdef CGAL_MESH_3_PROFILING
456 double facet_ref_time = t.elapsed();
457 std::cerr << "==== Facet refinement: " << facet_ref_time << " seconds ===="
458 << std::endl << std::endl;
459 # ifdef CGAL_MESH_3_EXPORT_PERFORMANCE_DATA
460 // If it's parallel but the refinement is forced to sequential, we don't
461 // output the value
462 # ifndef CGAL_DEBUG_FORCE_SEQUENTIAL_MESH_REFINEMENT
463 CGAL_MESH_3_SET_PERFORMANCE_DATA("Facets_time", facet_ref_time);
464 # endif
465 # endif
466 #endif
467
468 #if defined(CHECK_AND_DISPLAY_THE_NUMBER_OF_BAD_ELEMENTS_IN_THE_END)
469 std::cerr << std::endl
470 << "===============================================================" << std::endl
471 << "=== CHECK_AND_DISPLAY_THE_NUMBER_OF_BAD_ELEMENTS_IN_THE_END ===" << std::endl;
472 display_number_of_bad_elements();
473 std::cerr
474 << "==============================================================="
475 << std::endl << std::endl;
476 #endif
477
478 // Then activate facet to surface visitor (surface could be
479 // refined again if it is encroached)
480 facets_visitor_.activate();
481
482 dump_c3t3(r_c3t3_, dump_after_refine_surface_prefix);
483
484 if(!forced_stop())
485 {
486 // Then scan volume and refine it
487 CGAL_MESH_3_TASK_BEGIN(scan_cells_task_handle);
488 cells_mesher_.scan_triangulation();
489 CGAL_MESH_3_TASK_END(scan_cells_task_handle);
490 refinement_stage = REFINE_ALL;
491 #ifdef CGAL_MESH_3_PROFILING
492 std::cerr << "Refining cells..." << std::endl;
493 t.reset();
494 #endif
495 CGAL_MESH_3_TASK_BEGIN(refine_volume_mesh_task_handle);
496 cells_mesher_.refine(cells_visitor_);
497 CGAL_MESH_3_TASK_END(refine_volume_mesh_task_handle);
498 #ifdef CGAL_MESH_3_PROFILING
499 double cell_ref_time = t.elapsed();
500 std::cerr << "==== Cell refinement: " << cell_ref_time << " seconds ===="
501 << std::endl << std::endl;
502 # ifdef CGAL_MESH_3_EXPORT_PERFORMANCE_DATA
503 // If it's parallel but the refinement is forced to sequential, we don't
504 // output the value
505 # ifndef CGAL_DEBUG_FORCE_SEQUENTIAL_MESH_REFINEMENT
506 CGAL_MESH_3_SET_PERFORMANCE_DATA("Cells_refin_time", cell_ref_time);
507 # endif
508 # endif
509 #endif
510
511 #if defined(CGAL_MESH_3_VERBOSE) || defined(CGAL_MESH_3_PROFILING)
512 std::cerr
513 << "Vertices: " << r_c3t3_.triangulation().number_of_vertices() << std::endl
514 << "Facets : " << r_c3t3_.number_of_facets_in_complex() << std::endl
515 << "Tets : " << r_c3t3_.number_of_cells_in_complex() << std::endl;
516 #endif
517 } // end test of `maximal_number_of_vertices`
518 #else // ifdef CGAL_MESH_3_VERBOSE
519 std::cerr << "Start surface scan...";
520 CGAL_MESH_3_TASK_BEGIN(initialize_task_handle);
521 initialize();
522 CGAL_MESH_3_TASK_END(initialize_task_handle);
523 std::cerr << "end scan. [Bad facets:" << facets_mesher_.size() << "]";
524 std::cerr << std::endl << std::endl;
525 elapsed_time += timer.time();
526 timer.stop(); timer.reset(); timer.start();
527
528 int nbsteps = 0;
529
530 std::cerr << "Refining Surface...\n";
531 std::cerr << "Legend of the following line: "
532 << "(#vertices,#steps," << cells_mesher_.debug_info_header()
533 << ")\n";
534
535 std::cerr << "(" << r_tr.number_of_vertices() << ","
536 << nbsteps << "," << cells_mesher_.debug_info() << ")";
537
538 CGAL_MESH_3_TASK_BEGIN(refine_surface_mesh_task_handle);
539 while ( ! facets_mesher_.is_algorithm_done() &&
540 ! forced_stop() )
541 {
542 facets_mesher_.one_step(facets_visitor_);
543 std::cerr
544 << boost::format("\r \r"
545 "(%1%,%2%,%3%) (%|4$.1f| vertices/s)")
546 % r_tr.number_of_vertices()
547 % nbsteps % cells_mesher_.debug_info()
548 % (nbsteps / timer.time());
549 if(refinement_stage == REFINE_FACETS &&
550 ! forced_stop() &&
551 facets_mesher_.is_algorithm_done())
552 {
553 facets_mesher_.scan_edges();
554 refinement_stage = REFINE_FACETS_AND_EDGES;
555 }
556 if(refinement_stage == REFINE_FACETS_AND_EDGES &&
557 ! forced_stop() &&
558 facets_mesher_.is_algorithm_done())
559 {
560 facets_mesher_.scan_vertices();
561 refinement_stage = REFINE_FACETS_AND_EDGES_AND_VERTICES;
562 }
563 ++nbsteps;
564 }
565 CGAL_MESH_3_TASK_END(refine_surface_mesh_task_handle);
566 std::cerr << std::endl;
567 std::cerr << "Total refining surface time: " << timer.time() << "s" << std::endl;
568 std::cerr << std::endl;
569
570 CGAL_triangulation_postcondition(r_tr.is_valid());
571
572 elapsed_time += timer.time();
573 timer.stop(); timer.reset(); timer.start();
574 nbsteps = 0;
575
576 facets_visitor_.activate();
577 dump_c3t3(r_c3t3_, dump_after_refine_surface_prefix);
578 std::cerr << "Start volume scan...";
579 CGAL_MESH_3_TASK_BEGIN(scan_cells_task_handle);
580 cells_mesher_.scan_triangulation();
581 CGAL_MESH_3_TASK_END(scan_cells_task_handle);
582 refinement_stage = REFINE_ALL;
583 std::cerr << "end scan. [Bad tets:" << cells_mesher_.size() << "]";
584 std::cerr << std::endl << std::endl;
585 elapsed_time += timer.time();
586 timer.stop(); timer.reset(); timer.start();
587
588 std::cerr << "Refining...\n";
589 std::cerr << "Legend of the following line: "
590 << "(#vertices,#steps," << cells_mesher_.debug_info_header()
591 << ")\n";
592 std::cerr << "(" << r_tr.number_of_vertices() << ","
593 << nbsteps << "," << cells_mesher_.debug_info() << ")";
594
595 CGAL_MESH_3_TASK_BEGIN(refine_volume_mesh_task_handle);
596 while ( ! cells_mesher_.is_algorithm_done() &&
597 ! forced_stop() )
598 {
599 cells_mesher_.one_step(cells_visitor_);
600 std::cerr
601 << boost::format("\r \r"
602 "(%1%,%2%,%3%) (%|4$.1f| vertices/s)")
603 % r_tr.number_of_vertices()
604 % nbsteps % cells_mesher_.debug_info()
605 % (nbsteps / timer.time());
606 ++nbsteps;
607 }
608 CGAL_MESH_3_TASK_END(refine_volume_mesh_task_handle);
609 std::cerr << std::endl;
610
611 std::cerr << "Total refining volume time: " << timer.time() << "s" << std::endl;
612 std::cerr << "Total refining time: " << timer.time()+elapsed_time << "s" << std::endl;
613 std::cerr << std::endl;
614
615 CGAL_triangulation_postcondition(r_tr.is_valid());
616 #endif
617
618 (void)(forced_stop()); // sets *error_code
619
620 timer.stop();
621 elapsed_time += timer.time();
622
623 #if defined(CHECK_AND_DISPLAY_THE_NUMBER_OF_BAD_ELEMENTS_IN_THE_END) \
624 || defined(SHOW_REMAINING_BAD_ELEMENT_IN_RED)
625 std::cerr << std::endl
626 << "===============================================================" << std::endl
627 << "=== CHECK_AND_DISPLAY_THE_NUMBER_OF_BAD_ELEMENTS_IN_THE_END ===" << std::endl;
628 display_number_of_bad_elements();
629 std::cerr
630 << "==============================================================="
631 << std::endl << std::endl;
632 #endif
633
634 return elapsed_time;
635 }
636
637
638 template<class C3T3, class MC, class MD>
639 void
640 Mesher_3<C3T3,MC,MD>::
initialize()641 initialize()
642 {
643 #ifdef CGAL_MESH_3_PROFILING
644 std::cerr << "Initializing... ";
645 WallClockTimer t;
646 #endif
647 //=====================================
648 // Bounding box estimation
649 //=====================================
650 #if defined(CGAL_LINKED_WITH_TBB) || \
651 defined(CGAL_SEQUENTIAL_MESH_3_ADD_OUTSIDE_POINTS_ON_A_FAR_SPHERE)
652
653 #ifndef CGAL_SEQUENTIAL_MESH_3_ADD_OUTSIDE_POINTS_ON_A_FAR_SPHERE
654 if(boost::is_convertible<Concurrency_tag, Parallel_tag>::value)
655 #endif // If that macro is defined, then estimated_bbox must be initialized
656 {
657 Base::set_bbox(r_oracle_.bbox());
658 }
659 #endif // CGAL_LINKED_WITH_TBB||"sequential use far sphere"
660
661 //========================================
662 // Initialization: parallel or sequential
663 //========================================
664
665 #ifdef CGAL_LINKED_WITH_TBB
666 // Parallel
667 if (boost::is_convertible<Concurrency_tag, Parallel_tag>::value)
668 {
669 // we're not multi-thread, yet
670 r_c3t3_.triangulation().set_lock_data_structure(0);
671
672 # ifndef CGAL_PARALLEL_MESH_3_DO_NOT_ADD_OUTSIDE_POINTS_ON_A_FAR_SPHERE
673
674 if (r_c3t3_.number_of_far_points() == 0 && r_c3t3_.number_of_facets() == 0)
675 {
676 const Bbox_3 &bbox = r_oracle_.bbox();
677
678 // Compute radius for far sphere
679 const double xdelta = bbox.xmax()-bbox.xmin();
680 const double ydelta = bbox.ymax()-bbox.ymin();
681 const double zdelta = bbox.zmax()-bbox.zmin();
682 const double radius = 5. * std::sqrt(xdelta*xdelta +
683 ydelta*ydelta +
684 zdelta*zdelta);
685 const Vector center(
686 bbox.xmin() + 0.5*xdelta,
687 bbox.ymin() + 0.5*ydelta,
688 bbox.zmin() + 0.5*zdelta);
689 # ifdef CGAL_CONCURRENT_MESH_3_VERBOSE
690 std::cerr << "Adding points on a far sphere (radius = " << radius <<")...";
691 # endif
692 Random_points_on_sphere_3<Bare_point> random_point(radius);
693 const int NUM_PSEUDO_INFINITE_VERTICES = static_cast<int>(
694 float(std::thread::hardware_concurrency())
695 * Concurrent_mesher_config::get().num_pseudo_infinite_vertices_per_core);
696 for (int i = 0 ; i < NUM_PSEUDO_INFINITE_VERTICES ; ++i, ++random_point)
697 r_c3t3_.add_far_point(r_c3t3_.triangulation().geom_traits().construct_weighted_point_3_object()
698 (r_c3t3_.triangulation().geom_traits().construct_translated_point_3_object()(*random_point, center)));
699
700 # ifdef CGAL_CONCURRENT_MESH_3_VERBOSE
701 std::cerr << "done." << std::endl;
702 # endif
703 }
704 # endif // CGAL_PARALLEL_MESH_3_DO_NOT_ADD_OUTSIDE_POINTS_ON_A_FAR_SPHERE
705
706 #ifdef CGAL_MESH_3_PROFILING
707 double init_time = t.elapsed();
708 std::cerr << "done in " << init_time << " seconds." << std::endl;
709 #endif
710
711 // Scan triangulation
712 facets_mesher_.scan_triangulation();
713 refinement_stage = REFINE_FACETS;
714
715 // From now on, we're multi-thread
716 r_c3t3_.triangulation().set_lock_data_structure(get_lock_data_structure());
717 }
718 // Sequential
719 else
720 #endif // CGAL_LINKED_WITH_TBB
721 {
722 if (r_c3t3_.number_of_far_points() == 0 &&
723 r_c3t3_.number_of_facets() == 0 &&
724 (r_c3t3_.triangulation().dimension() < 3
725 #ifdef CGAL_SEQUENTIAL_MESH_3_ADD_OUTSIDE_POINTS_ON_A_FAR_SPHERE
726 || true // in sequential mode, one wants the far points only
727 // if the macro
728 // `CGAL_SEQUENTIAL_MESH_3_ADD_OUTSIDE_POINTS_ON_A_FAR_SPHERE`
729 // is defined, or if the dimension of the
730 // triangulation is 2.
731 #endif
732 )) {
733 const Bbox_3 &bbox = r_oracle_.bbox();
734
735 // Compute radius for far sphere
736 const double xdelta = bbox.xmax()-bbox.xmin();
737 const double ydelta = bbox.ymax()-bbox.ymin();
738 const double zdelta = bbox.zmax()-bbox.zmin();
739 const double radius = 5. * std::sqrt(xdelta*xdelta +
740 ydelta*ydelta +
741 zdelta*zdelta);
742 const Vector center(
743 bbox.xmin() + 0.5*xdelta,
744 bbox.ymin() + 0.5*ydelta,
745 bbox.zmin() + 0.5*zdelta);
746 # ifdef CGAL_MESH_3_VERBOSE
747 std::cerr << "Adding points on a far sphere (radius = " << radius << ")...";
748 # endif
749 Random_points_on_sphere_3<Bare_point> random_point(radius);
750 const int NUM_PSEUDO_INFINITE_VERTICES = 12*2;
751 for (int i = 0 ; i < NUM_PSEUDO_INFINITE_VERTICES ; ++i, ++random_point)
752 r_c3t3_.add_far_point(r_c3t3_.triangulation().geom_traits().construct_weighted_point_3_object()
753 (r_c3t3_.triangulation().geom_traits().construct_translated_point_3_object()(*random_point, center)));
754 # ifdef CGAL_MESH_3_VERBOSE
755 std::cerr << "done." << std::endl;
756 # endif
757 }
758
759 #ifdef CGAL_MESH_3_PROFILING
760 double init_time = t.elapsed();
761 std::cerr << "done in " << init_time << " seconds." << std::endl;
762 #endif
763
764 // Scan triangulation
765 facets_mesher_.scan_triangulation();
766 refinement_stage = REFINE_FACETS;
767 }
768 }
769
770
771 template<class C3T3, class MC, class MD>
772 void
773 Mesher_3<C3T3,MC,MD>::
fix_c3t3()774 fix_c3t3()
775 {
776 if ( ! facets_visitor_.is_active() )
777 {
778 cells_mesher_.scan_triangulation();
779 }
780 }
781
782
783 template<class C3T3, class MC, class MD>
784 void
785 Mesher_3<C3T3,MC,MD>::
display_number_of_bad_elements()786 display_number_of_bad_elements()
787 {
788 int nf = facets_mesher_.number_of_bad_elements();
789 int nc = cells_mesher_.number_of_bad_elements();
790 std::cerr << "Bad facets: " << nf << " - Bad cells: " << nc << std::endl;
791 }
792
793 template<class C3T3, class MC, class MD>
794 void
795 Mesher_3<C3T3,MC,MD>::
one_step()796 one_step()
797 {
798 if ( ! facets_visitor_.is_active() )
799 {
800 facets_mesher_.one_step(facets_visitor_);
801
802 if ( facets_mesher_.is_algorithm_done() )
803 {
804 if(forced_stop()) {
805 return;
806 }
807
808 switch(refinement_stage) {
809 case REFINE_FACETS:
810 facets_mesher_.scan_edges();
811 refinement_stage = REFINE_FACETS_AND_EDGES;
812 if (!facets_mesher_.is_algorithm_done()) break;
813 CGAL_FALLTHROUGH;
814 case REFINE_FACETS_AND_EDGES:
815 facets_mesher_.scan_vertices();
816 refinement_stage = REFINE_FACETS_AND_EDGES_AND_VERTICES;
817 if (!facets_mesher_.is_algorithm_done()) break;
818 CGAL_FALLTHROUGH;
819 default:
820 facets_visitor_.activate();
821 cells_mesher_.scan_triangulation();
822 refinement_stage = REFINE_ALL;
823 }
824 }
825 }
826 else
827 {
828 CGAL_assertion(refinement_stage == REFINE_ALL);
829 cells_mesher_.one_step(cells_visitor_);
830 }
831 }
832
833 template<class C3T3, class MC, class MD>
834 bool
835 Mesher_3<C3T3,MC,MD>::
is_algorithm_done()836 is_algorithm_done()
837 {
838 return cells_mesher_.is_algorithm_done();
839 }
840
841
842 #ifdef CGAL_MESH_3_MESHER_STATUS_ACTIVATED
843 template<class C3T3, class MC, class MD>
844 typename Mesher_3<C3T3,MC,MD>::Mesher_status
845 Mesher_3<C3T3,MC,MD>::
status()846 status() const
847 {
848 #ifdef CGAL_LINKED_WITH_TBB
849 if(boost::is_convertible<Concurrency_tag, Parallel_tag>::value) {
850 return Mesher_status(
851 # if CGAL_CONCURRENT_COMPACT_CONTAINER_APPROXIMATE_SIZE
852 approximate_number_of_vertices(Concurrency_tag()),
853 #else
854 // not thread-safe, but that is not important
855 approximate_number_of_vertices(CGAL::Sequential_tag()),
856 #endif
857 0,
858 0);
859 }
860 else
861 #endif // with TBB
862 {
863 return Mesher_status(r_c3t3_.triangulation().number_of_vertices(),
864 facets_mesher_.queue_size(),
865 cells_mesher_.queue_size());
866 }
867 }
868 #endif
869
870 template<class C3T3, class MC, class MD>
871 inline
872 std::string
debug_info()873 Mesher_3<C3T3,MC,MD>::debug_info() const
874 {
875 return cells_mesher_.debug_info();
876 }
877
878 template<class C3T3, class MC, class MD>
879 inline
880 std::string
debug_info_header()881 Mesher_3<C3T3,MC,MD>::debug_info_header() const
882 {
883 return cells_mesher_.debug_info_header();
884 }
885
886 } // end namespace Mesh_3
887
888 } // end namespace CGAL
889
890 #include <CGAL/enable_warnings.h>
891
892 #endif // CGAL_MESH_3_MESHER_3_H
893