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