1 // Copyright (c) 2009, 2014 INRIA Sophia-Antipolis (France).
2 // Copyright (c) 2017 GeometryFactory (France).
3 // All rights reserved.
4 //
5 // This file is part of CGAL (www.cgal.org).
6 //
7 // $URL: https://github.com/CGAL/cgal/blob/v5.3/Periodic_3_mesh_3/include/CGAL/refine_periodic_3_mesh_3.h $
8 // $Id: refine_periodic_3_mesh_3.h 53d4c9b 2019-10-28T11:29:08+01:00 Mael Rouxel-Labbé
9 // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
10 //
11 // Author(s)     : Stephane Tayeb,
12 //                 Mikhail Bogdanov,
13 //                 Mael Rouxel-Labbé
14 //
15 #ifndef CGAL_REFINE_PERIODIC_3_MESH_3_H
16 #define CGAL_REFINE_PERIODIC_3_MESH_3_H
17 
18 #include <CGAL/license/Periodic_3_mesh_3.h>
19 
20 #include <CGAL/Mesh_3/config.h>
21 #include <CGAL/Periodic_3_mesh_3/config.h>
22 
23 #include <CGAL/Mesh_3/C3T3_helpers.h>
24 #include <CGAL/Mesh_3/Dump_c3t3.h>
25 #include <CGAL/Mesh_3/Triangulation_helpers.h>
26 #include <CGAL/refine_mesh_3.h>
27 #include <CGAL/Time_stamper.h>
28 
29 #include <boost/parameter/preprocessor.hpp>
30 #include <boost/unordered_set.hpp>
31 
32 #include <algorithm>
33 #include <iostream>
34 #include <iterator>
35 
36 namespace CGAL {
37 
38 namespace internal {
39 
40 template<class C3T3, class MeshDomain>
project_dummy_points_of_surface(C3T3 & c3t3,const MeshDomain & domain)41 void project_dummy_points_of_surface(C3T3& c3t3, const MeshDomain& domain)
42 {
43   typedef typename C3T3::Vertex_handle                     Vertex_handle;
44   typedef CGAL::Hash_handles_with_or_without_timestamps    Hash_fct;
45   typedef boost::unordered_set<Vertex_handle, Hash_fct>    Vertex_container;
46 
47   Vertex_container vertex_container;
48   find_points_to_project(c3t3, std::insert_iterator<Vertex_container>(vertex_container, vertex_container.begin()));
49 
50   project_points(c3t3, domain, vertex_container.begin(), vertex_container.end());
51 }
52 
53 template<class C3T3, class OutputIterator>
find_points_to_project(C3T3 & c3t3,OutputIterator vertices)54 void find_points_to_project(C3T3& c3t3, OutputIterator vertices)
55 {
56   typedef typename C3T3::Vertex_handle Vertex_handle;
57   typedef typename C3T3::Cell_handle Cell_handle;
58 
59   // Don't project all dummy points, but simply the ones that are involved in
60   // surface facets of the c3t3.
61   for(typename C3T3::Facets_in_complex_iterator face_it = c3t3.facets_in_complex_begin();
62                                                 face_it != c3t3.facets_in_complex_end();
63                                                 ++face_it)
64   {
65     int ind = face_it->second;
66     Cell_handle c = face_it->first;
67 
68     for(int i = 1; i < 4; i++) {
69       Vertex_handle v = c->vertex((ind+i)&3);
70 
71       typename C3T3::Index index = c3t3.index(v);
72       if(const int* i = boost::get<int>(&index))
73       {
74         if(*i == 0) // '0' is the index of dummies
75           *vertices++ = v;
76       }
77     }
78   }
79 }
80 
81 template<class C3T3, class MeshDomain, class InputIterator>
project_points(C3T3 & c3t3,const MeshDomain & domain,InputIterator vertex_begin,InputIterator vertex_end)82 void project_points(C3T3& c3t3,
83                     const MeshDomain& domain,
84                     InputIterator vertex_begin,
85                     InputIterator vertex_end)
86 {
87   typedef typename C3T3::Vertex_handle         Vertex_handle;
88 
89   typedef typename C3T3::Triangulation         Tr;
90 
91   typedef typename Tr::Bare_point              Bare_point;
92   typedef typename Tr::Weighted_point          Weighted_point;
93   typedef typename Tr::Geom_traits::Vector_3   Vector_3;
94   typedef typename Tr::Geom_traits::FT         FT;
95 
96   typename C3T3::Triangulation::Geom_traits::Construct_point_3 cp =
97     c3t3.triangulation().geom_traits().construct_point_3_object();
98 
99   CGAL::Mesh_3::C3T3_helpers<C3T3, MeshDomain> helper(c3t3, domain);
100   CGAL::Mesh_3::Triangulation_helpers<Tr> tr_helpers;
101 
102   for(InputIterator it = vertex_begin; it != vertex_end; ++it)
103   {
104     Vertex_handle vh = *it;
105 
106     const Weighted_point& vh_wp = c3t3.triangulation().point(vh);
107     const Bare_point& vh_p = cp(vh_wp);
108     const Bare_point new_point = helper.project_on_surface(vh, vh_p);
109 
110     const FT sq_d = CGAL::squared_distance(new_point, vh_p);
111 
112 #ifdef CGAL_PERIODIC_3_MESH_3_DEBUG_DUMMY_PROJECTION
113     std::cerr << "vh: " << &*vh << std::endl;
114     std::cerr << "vhp: " << vh_p << std::endl;
115     std::cerr << "projected: " << new_point << std::endl;
116     std::cerr << "squared distance from dummy to surface: " << sq_d << std::endl;
117 #endif
118 
119     // Skip tiny moves for efficiency
120     if(sq_d < 1e-10) // arbitrary value, maybe compare it to the surface distance criterium ?
121       continue;
122 
123     // Do not project if the projected point is in a protection ball
124     if(tr_helpers.inside_protecting_balls(c3t3.triangulation(), vh, new_point))
125       continue;
126 
127     const Vector_3 move(vh_p, new_point);
128     Vertex_handle new_vertex = helper.update_mesh(vh, move);
129     if(new_vertex != vh) // if the move has successfully been performed
130       c3t3.set_dimension(new_vertex, 2);
131   }
132 }
133 
134 } // namespace internal
135 
136 #if defined(BOOST_MSVC)
137 #  pragma warning(push)
138 #  pragma warning(disable:4003) // not enough actual parameters for macro
139 #endif
140 
141 // see <CGAL/config.h>
142 CGAL_PRAGMA_DIAG_PUSH
143 // see <CGAL/boost/parameter.h>
144 CGAL_IGNORE_BOOST_PARAMETER_NAME_WARNINGS
145 
146 BOOST_PARAMETER_FUNCTION(
147   (void),
148   refine_periodic_3_mesh_3,
149   parameters::tag,
150   (required (in_out(c3t3),*) (domain,*) (criteria,*) ) // nondeduced
151   (deduced
152     (optional
153       (exude_param, (parameters::internal::Exude_options), parameters::no_exude()) // another default parameter distinct from Mesh_3
154       (perturb_param, (parameters::internal::Perturb_options), parameters::no_perturb()) // another default parameter distinct from Mesh_3
155       (odt_param, (parameters::internal::Odt_options), parameters::no_odt())
156       (lloyd_param, (parameters::internal::Lloyd_options), parameters::no_lloyd())
157       (reset_param, (parameters::Reset), parameters::reset_c3t3())
158       (mesh_options_param, (parameters::internal::Mesh_3_options),
159        parameters::internal::Mesh_3_options())
160       (manifold_options_param, (parameters::internal::Manifold_options),
161        parameters::internal::Manifold_options())
162     )
163   )
164 )
165 {
166   return refine_periodic_3_mesh_3_impl(c3t3, domain, criteria,
167                                        exude_param,
168                                        perturb_param,
169                                        odt_param,
170                                        lloyd_param,
171                                        reset_param(),
172                                        mesh_options_param,
173                                        manifold_options_param);
174 }
175 
176 CGAL_PRAGMA_DIAG_POP
177 
178 #if defined(BOOST_MSVC)
179 #  pragma warning(pop)
180 #endif
181 
182 /**
183  * @brief This function refines the mesh c3t3 wrt domain & criteria
184  *
185  * @param c3t3 the mesh to be refined.
186  * @param domain the domain to be discretized
187  * @param criteria the criteria
188  * @param exude if \c true, an exudation step will be done at
189  *   the end of the Delaunay refinement process
190  * @param perturb if \c true, an explicit vertex perturbation step will be
191  *   done at the end of refinement process
192  * @param reset_c3t3 if \c true, a new C3T3 will be construct from param c3t3.
193  *   The new c3t3 keeps only the vertices (as NON-weighted points with their
194  *   dimension and Index) of the triangulation. That allows to refine a mesh
195  *   which has been exuded.
196  * @param mesh_3_options is a struct object used to pass non-documented options,
197  *   for debugging purpose.
198  */
199 template<class C3T3, class MeshDomain, class MeshCriteria>
200 void refine_periodic_3_mesh_3_impl(C3T3& c3t3,
201                                    const MeshDomain& domain,
202                                    const MeshCriteria& criteria,
203                                    const parameters::internal::Exude_options& exude,
204                                    const parameters::internal::Perturb_options& perturb,
205                                    const parameters::internal::Odt_options& odt,
206                                    const parameters::internal::Lloyd_options& lloyd,
207                                    bool reset_c3t3,
208                                    const parameters::internal::Mesh_3_options&
209                                      mesh_options = parameters::internal::Mesh_3_options(),
210                                    const parameters::internal::Manifold_options&
211                                      manifold_options = parameters::internal::Manifold_options())
212 {
213   // 'refine_periodic_3_mesh_3' does not insert dummy points, so the triangulation
214   // must already be 1-sheeted prior to calls to this function.
215   CGAL_precondition(c3t3.triangulation().is_1_cover());
216 
217   typedef Mesh_3::Mesher_3<C3T3, MeshCriteria, MeshDomain> Mesher;
218 
219   // Reset c3t3 (i.e. remove weights) if needed
220   if ( reset_c3t3 )
221   {
222     C3T3 tmp_c3t3;
223     std::for_each(c3t3.triangulation().finite_vertices_begin(),
224                   c3t3.triangulation().finite_vertices_end(),
225                   details::Insert_vertex_in_c3t3<C3T3>(tmp_c3t3));
226     // TODO: corners and edges are not restored
227     c3t3.swap(tmp_c3t3);
228   }
229 
230   dump_c3t3(c3t3, mesh_options.dump_after_init_prefix);
231 
232   // Build mesher and launch refinement process
233   Mesher mesher(c3t3, domain, criteria,
234                 manifold_options.mesh_topology,
235                 mesh_options.maximal_number_of_vertices,
236                 mesh_options.pointer_to_error_code);
237 
238   double refine_time = mesher.refine_mesh(mesh_options.dump_after_refine_surface_prefix);
239   c3t3.clear_manifold_info();
240 
241   CGAL_expensive_postcondition(c3t3.triangulation().tds().is_valid());
242   CGAL_expensive_postcondition(c3t3.triangulation().is_valid());
243   CGAL_expensive_postcondition(c3t3.is_valid());
244   dump_c3t3(c3t3, mesh_options.dump_after_init_prefix);
245 
246   // Project dummy points on the surface to lessen their influence on the output
247   internal::project_dummy_points_of_surface(c3t3, domain);
248 
249   CGAL_expensive_postcondition(c3t3.triangulation().tds().is_valid());
250   CGAL_expensive_postcondition(c3t3.triangulation().is_valid());
251   CGAL_expensive_postcondition(c3t3.is_valid());
252 
253   // Odt
254   if(odt)
255   {
256     odt_optimize_mesh_3(c3t3, domain,
257                         parameters::time_limit = odt.time_limit(),
258                         parameters::max_iteration_number = odt.max_iteration_number(),
259                         parameters::convergence = odt.convergence(),
260                         parameters::freeze_bound = odt.bound());
261   }
262 
263   // Lloyd
264   if(lloyd)
265   {
266     lloyd_optimize_mesh_3(c3t3, domain,
267                           parameters::time_limit = lloyd.time_limit(),
268                           parameters::max_iteration_number = lloyd.max_iteration_number(),
269                           parameters::convergence = lloyd.convergence(),
270                           parameters::freeze_bound = lloyd.bound());
271   }
272 
273   if(odt || lloyd)
274   {
275     dump_c3t3(c3t3, mesh_options.dump_after_glob_opt_prefix);
276   }
277 
278   // Perturbation
279   if(perturb)
280   {
281     double perturb_time_limit = refine_time;
282 
283     if(perturb.is_time_limit_set())
284       perturb_time_limit = perturb.time_limit();
285 
286     perturb_mesh_3(c3t3, domain,
287                    parameters::time_limit = perturb_time_limit,
288                    parameters::sliver_bound = perturb.bound());
289 
290     dump_c3t3(c3t3, mesh_options.dump_after_perturb_prefix);
291   }
292 
293   // Exudation
294   if(exude)
295   {
296     double exude_time_limit = refine_time;
297 
298     if(exude.is_time_limit_set())
299       exude_time_limit = exude.time_limit();
300 
301     exude_mesh_3(c3t3,
302                  parameters::time_limit = exude_time_limit,
303                  parameters::sliver_bound = exude.bound());
304 
305     dump_c3t3(c3t3, mesh_options.dump_after_perturb_prefix);
306   }
307 
308   CGAL_expensive_postcondition(c3t3.triangulation().tds().is_valid());
309   CGAL_expensive_postcondition(c3t3.triangulation().is_valid());
310   CGAL_expensive_postcondition(c3t3.is_valid());
311 }
312 
313 } // end namespace CGAL
314 
315 #endif // CGAL_REFINE_PERIODIC_3_MESH_3_H
316