GCC Code Coverage Report

Line	Branch	Exec	Source
1			/*
2			* Copyright (c) 2000-2022 Inria
3			* All rights reserved.
4			*
5			* Redistribution and use in source and binary forms, with or without
6			* modification, are permitted provided that the following conditions are met:
7			*
8			* * Redistributions of source code must retain the above copyright notice,
9			* this list of conditions and the following disclaimer.
10			* * Redistributions in binary form must reproduce the above copyright notice,
11			* this list of conditions and the following disclaimer in the documentation
12			* and/or other materials provided with the distribution.
13			* * Neither the name of the ALICE Project-Team nor the names of its
14			* contributors may be used to endorse or promote products derived from this
15			* software without specific prior written permission.
16			*
17			* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18			* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19			* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20			* ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
21			* LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22			* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23			* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24			* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25			* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26			* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27			* POSSIBILITY OF SUCH DAMAGE.
28			*
29			* Contact: Bruno Levy
30			*
31			* https://www.inria.fr/fr/bruno-levy
32			*
33			* Inria,
34			* Domaine de Voluceau,
35			* 78150 Le Chesnay - Rocquencourt
36			* FRANCE
37			*
38			*/
39
40			#include <exploragram/hexdom/hex_cruncher.h>
41			#include <exploragram/hexdom/intersect_tools.h>
42			#include <exploragram/hexdom/polygon.h>
43			#include <exploragram/hexdom/mesh_inspector.h>
44			#include <geogram/basic/geometry.h>
45			#include <geogram/points/colocate.h>
46			#include <geogram/mesh/triangle_intersection.h>
47			#include <geogram/mesh/mesh_repair.h> //used in bourrin subdivide hex
48			#define FPG_UNCERTAIN_VALUE 0
49			#include <geogram/numerics/predicates/orient3d.h>
50
51			#ifdef GEO_COMPILER_MSVC
52			#include <intrin.h>
53			void nico_assert(bool b) {
54			if (!b) __debugbreak();
55			geo_assert(b);
56			}
57			#else
58			#define nico_assert(b) geo_assert(b)
59			#endif
60
61			namespace GEO {
62
63
64			// Extra connectivity dedicaed to surfaces with triangles and quads only.
65			// Halfedges are indiced by 4* facet + local_id => requires padding for triangles
66			struct QuadExtraConnectivity {
67		✗	void init(Mesh* p_m) {
68		✗	m = p_m;
69		✗	FOR(f, m->facets.nb()) {
70		✗	nico_assert(m->facets.nb_vertices(f) == 4 || m->facets.nb_vertices(f) == 3);// check that surface is quadragulated
71
72			}
73		✗	opp_h.clear();
74		✗	opp_h.resize(4 * m->facets.nb(), NOT_AN_ID); // NOT facet_corners !!!
75		✗	create_non_manifold_facet_adjacence(m);
76		✗	FOR(f, m->facets.nb()) FOR(lc, m->facets.nb_vertices(f)) {
77		✗	index_t opp_f = m->facets.adjacent(f, lc);
78		✗	if (opp_f == NOT_AN_ID) { plop("bad adjacency detected "); continue; }
79
80		✗	index_t opp_lc = NOT_AN_ID;
81		✗	FOR(opp_lc_it, m->facets.nb_vertices(opp_f))
82		✗	if (f == m->facets.adjacent(opp_f, opp_lc_it)
83		✗	&& m->facets.vertex(f, (lc + 1) % m->facets.nb_vertices(f)) == m->facets.vertex(opp_f, opp_lc_it)
84		✗	) opp_lc = opp_lc_it;
85		✗	if (opp_lc == NOT_AN_ID) plop("not a symetric opposite !");
86		✗	set_opp(4 * f + lc, 4 * opp_f + opp_lc);
87		✗	if (vertex(4 * f + lc) != vertex(next(opp(4 * f + lc)))) plop("wrong opposites");
88		✗	if (vertex(next(4 * f + lc)) != vertex(opp(4 * f + lc))) plop("wrong opposites");
89			}
90
91		✗	}
92
93
94			void check_integrity() {
95			FOR(h, 4 * m->facets.nb()) if (valid(h)) {
96			FOR(d, 5) nico_assert(valid(next(h, d)));
97			nico_assert(valid(opp(h)));
98			nico_assert(valid(next_around_vertex(h)));
99			nico_assert(valid(next_around_vertex(next_around_vertex(h))));
100			}
101			}
102
103		✗	void debug_export_adjacence() {
104		✗	FOR(f, m->facets.nb()) FOR(lc, m->facets.nb_vertices(f)) {
105		✗	if (opp(4 * f + lc) != NOT_AN_ID)
106		✗	m->facets.set_adjacent(f, lc, face(opp(4 * f + lc)));
107		✗	else m->facets.set_adjacent(f, lc, NOT_AN_ID);
108			}
109		✗	}
110
111		✗	bool valid(index_t h) { return h < 4 * m->facets.nb() && (h%4)<m->facets.nb_vertices(h/4); }
112
113		✗	index_t fsize(index_t e){ nico_assert(valid(e)); return m->facets.nb_vertices(face(e)); }
114		✗	void set_opp(index_t i, index_t j) { nico_assert(valid(i) && valid(j)); opp_h[i] = j; opp_h[j] = i; }
115		✗	index_t face(index_t e) { nico_assert(valid(e)); return e / 4; }
116		✗	index_t local_id(index_t e) { nico_assert(valid(e)); return e % 4; }
117		✗	index_t next(index_t e, index_t nb = 1) { nico_assert(valid(e)); return 4 * face(e) + ((e%4 + nb) % fsize(e)); }
118		✗	index_t opp(index_t e) { nico_assert(valid(e)); return opp_h[e]; }
119		✗	index_t vertex(index_t e) { nico_assert(valid(e)); return m->facets.vertex(face(e), local_id(e)); }
120		✗	index_t corner(index_t e) { nico_assert(valid(e)); return m->facets.corner(face(e), local_id(e)); }
121		✗	index_t next_around_vertex(index_t e) { nico_assert(valid(e)); return opp(next(e, fsize(e)-1)); }
122
123		✗	void set_vertex(index_t e, index_t v) { nico_assert(valid(e) && v < m->vertices.nb()); m->facets.set_vertex(face(e), local_id(e), v); }
124
125
126		✗	bool has_valid_one_ring(index_t e) {
127		✗	index_t cir = e;
128		✗	int count = 0;
129			do {
130		✗	if (cir == NOT_AN_ID) return false;
131		✗	if (count++ == 1000) return false;
132		✗	cir = next_around_vertex(cir);
133		✗	} while (cir != e);
134		✗	return true;
135			}
136		✗	int valence(index_t e) {
137		✗	nico_assert(has_valid_one_ring(e));
138		✗	index_t cir = e;
139		✗	int count = 0;
140			do {
141		✗	count++;
142		✗	cir = next_around_vertex(cir);
143		✗	} while (cir != e);
144		✗	return count;
145			}
146
147			bool is_closed() {
148			FOR(h, opp_h.size()) if (opp_h[h] == NOT_AN_ID) {
149			Attribute<double> deb(m->vertices.attributes(), "debug");
150			deb[vertex(h)] = 10;
151			deb[vertex(next(h))] = 10;
152			return false;
153			}
154			FOR(v, m->vertices.nb()) if (!has_valid_one_ring(v)) {
155			Attribute<double> deb(m->vertices.attributes(), "debug");
156			deb[v] = 10;
157			return false;
158			}
159			return true;
160			}
161
162
163			Mesh* m;
164			vector<index_t> opp_h;
165			};
166
167
168
169			struct CutSingularity {
170			Mesh* m; // ;)
171			vector<bool> visited; // prevents iterating more than once on the same cut
172			vec3 N; // normal to the current cut
173			QuadExtraConnectivity qem;
174			vector<index_t> border;
175			vector<vec3> pts;
176			vector<index_t> quads;
177
178		✗	CutSingularity(Mesh* p_m) {
179		✗	m = p_m;
180		✗	qem.init(m);
181		✗	visited.resize(4 * m->facets.nb(), false);
182		✗	}
183
184		✗	bool create_edge_loop(index_t h, vector<index_t>& test_border) {
185		✗	double sigma_angle = 0;
186		✗	index_t cir = h;
187			do {
188		✗	visited[cir] = true;
189		✗	if (qem.fsize(cir) == 3) return false;
190		✗	vec3 Nup = facet_normal(m, qem.face(cir));
191		✗	if (qem.opp(cir) == NOT_AN_ID) { test_border.clear(); break; }
192		✗	vec3 Ndown = facet_normal(m, qem.face(qem.opp(cir)));
193		✗	if (dot(N, Nup) > .5 || dot(N, Ndown) < -.5) return false;
194
195
196
197		✗	vec3 cir_dir = X(m)[qem.vertex(qem.next(cir))] - X(m)[qem.vertex(cir)];
198		✗	cir_dir = normalize(cir_dir);
199
200		✗	test_border.push_back(cir);
201		✗	index_t next = NOT_AN_ID;
202		✗	index_t in_cir = qem.next(cir);
203
204
205			//double best_dot = -1e20;
206		✗	double best_angle = -M_PI;
207			do {
208		✗	if (qem.fsize(in_cir)!=4 || qem.opp(in_cir) == NOT_AN_ID) {
209		✗	next = NOT_AN_ID;
210		✗	return false;
211			}
212
213		✗	vec3 in_cir_dir = X(m)[qem.vertex(qem.next(in_cir))] - X(m)[qem.vertex(in_cir)];
214		✗	in_cir_dir = normalize(in_cir_dir);
215
216		✗	if (fabs(dot(in_cir_dir, N)) < sin(M_PI / 8.) // stay in the cut plane (orthogonal to N)
217		✗	&& in_cir != qem.opp(cir) // do not go back
218			) {
219		✗	double newangle = atan2(dot(N, cross(cir_dir, in_cir_dir)), dot(in_cir_dir, cir_dir));
220		✗	if (best_angle < newangle) {
221		✗	best_angle = newangle;
222		✗	next = in_cir;
223			}
224			}
225		✗	in_cir = qem.next(qem.opp(in_cir));
226		✗	} while (in_cir != qem.next(cir));
227		✗	sigma_angle += best_angle;
228
229		✗	if (next == NOT_AN_ID || test_border.size() > 30) {
230		✗	return false;
231			}
232		✗	cir = next;
233		✗	} while (cir != h);
234		✗	if (sigma_angle < 0) return false;
235		✗	return true;
236			}
237
238		✗	bool cut_separates_2_shorts_closed_quads_strips(index_t h) {
239			// check that it separates two smalls and close quads strips
240		✗	index_t quad_strip_start[2] = { qem.next(h), qem.next(qem.opp(h)) };
241		✗	FOR(q, 2) {
242		✗	index_t it = quad_strip_start[q];
243		✗	int i = 0;
244			for (;;) {
245		✗	if (qem.fsize(it) == 3) return false;
246		✗	if (i == 30 || qem.opp(it) == NOT_AN_ID)
247		✗	return false;
248		✗	it = qem.next(qem.opp(it), 2);
249		✗	if (it == quad_strip_start[q]) break;
250		✗	i++;
251			}
252			}
253		✗	return true;
254			}
255
256
257		✗	bool can_easily_discard_current_edge_loop(vector<index_t>& test_border) {
258		✗	if (test_border.size() < 3) return true;
259		✗	if (!border.empty() && test_border.size() >= border.size()) return true;
260		✗	if (test_border.size() == 4
261		✗	&& (qem.face(qem.opp(test_border[0])) == qem.face(qem.opp(test_border[2]))
262		✗	|| qem.face(test_border[0]) == qem.face(test_border[2])
263			)
264		✗	) return true;
265		✗	if (test_border.size() % 2 != 0) return true;
266		✗	return false;
267			}
268
269
270
271		✗	bool apply() {
272		✗	vec3 best_N;
273		✗	double best_cut_area = 1e20;
274
275			// double ave_edge_length = get_facet_average_edge_size(m); // BL: unused.
276			//0;
277			//FOR(f, m->facets.nb()) FOR(v, 4) ave_edge_length += (X(m)[m->facets.vertex(f, v)] - X(m)[m->facets.vertex(f, (v + 1) % 4)]).length();
278			//ave_edge_length /= 4.*double(m->facets.nb());
279
280
281			// STEP 1: determine a valid cut along halfedges "border", its quadrangulation "quads", with vertices "pts[i]" (starting with vertices of "border")
282
283
284		✗	int nb_border_tried = 0;
285		✗	int nb_border_success = 0;
286		✗	Attribute<int> fail_c(m->facet_corners.attributes(), "fail_c");
287		✗	FOR(h, 4 * m->facets.nb()) {
288		✗	if (!qem.valid(h)) continue;
289		✗	fail_c[m->facets.corner(h / 4, h % 4)] = 0;
290		✗	fail_c[m->facets.corner(h/4,h%4)] = 10;
291
292		✗	if (visited[h]) continue;
293		✗	fail_c[m->facets.corner(h / 4, h % 4)] = 20;
294
295			// index_t f = qem.face(h); // [BL: unused]
296		✗	if (qem.fsize(h)!=4) continue;
297		✗	fail_c[m->facets.corner(h / 4, h % 4)] = 30;
298
299		✗	if (qem.opp(h) == NOT_AN_ID) continue;
300		✗	fail_c[m->facets.corner(h / 4, h % 4)] = 40;
301
302		✗	if (!cut_separates_2_shorts_closed_quads_strips(h)) continue;
303		✗	fail_c[m->facets.corner(h / 4, h % 4)] = 50;
304
305		✗	N = X(m)[qem.vertex(qem.next(h, 3))] - X(m)[qem.vertex(h)];
306		✗	N = normalize(N);
307		✗	vector<index_t> test_border;
308
309
310			// STEP 1.1 create the edge loop starting from h
311		✗	if (!create_edge_loop(h, test_border)) continue;
312
313		✗	fail_c[m->facets.corner(h / 4, h % 4)] = 60;
314
315
316
317
318			// STEP 1.2 check if the edge loop starting from h is valid and better than previous loop
319		✗	if (can_easily_discard_current_edge_loop(test_border)) continue;
320		✗	fail_c[m->facets.corner(h / 4, h % 4)] = 70;
321
322		✗	vector<vec3> test_pts;
323		✗	vector<index_t> test_quads;
324		✗	FOR(e, test_border.size()) test_pts.push_back(X(m)[qem.vertex(test_border[e])]);
325
326
327		✗	double cut_area = 0;
328		✗	FOR(e, test_pts.size()) cut_area -= dot(cross(test_pts[e], test_pts[(e + 1) % test_pts.size()]), N);
329		✗	if (best_cut_area < cut_area)continue;
330
331
332
333		✗	Poly3d p3(test_pts);
334		✗	nb_border_tried++;
335
336
337
338		✗	if (!p3.try_quadrangulate(test_quads)) {
339		✗	continue;
340			} else {
341		✗	nb_border_success++;
342		✗	bool cut_will_intersect = false;
343
344		✗	FacetIntersect finter(m);
345		✗	FOR(q, test_quads.size() / 4) {
346		✗	vector<vec3> quad;
347		✗	FOR(lc, 4) quad.push_back(test_pts[test_quads[4 * q + lc]]);
348		✗	cut_will_intersect = cut_will_intersect || finter.get_intersections(quad).size()>0;
349		✗	}
350
351		✗	if (cut_will_intersect) {
352		✗	continue;
353			}
354		✗	}
355
356			// STEP 1.3 validate the new loop
357			{
358		✗	border.swap(test_border);
359		✗	pts.swap(test_pts);
360		✗	test_quads.swap(quads);
361		✗	best_cut_area = cut_area;
362		✗	best_N = N;
363			}
364		✗	}
365
366		✗	plop(nb_border_tried);
367		✗	plop(nb_border_success);
368
369		✗	if (border.empty()) return false;
370
371		✗	vector<index_t> upper_v;
372		✗	vector<index_t> lower_v;
373			{
374		✗	index_t off_v = m->vertices.create_vertices(border.size());
375		✗	FOR(e, border.size()) {
376		✗	upper_v.push_back(qem.vertex(border[e]));
377		✗	pts.push_back(X(m)[upper_v[e]]);
378		✗	lower_v.push_back(off_v + e);
379		✗	X(m)[off_v + e] = pts[e];
380			};
381			}
382		✗	vector<vector<index_t> > opp_fan(border.size());
383		✗	FOR(e, border.size()) {
384		✗	index_t cir = qem.opp(border[e]);
385			do {
386		✗	opp_fan[e].push_back(cir);
387		✗	nico_assert(qem.fsize(cir) == 4);
388		✗	cir = qem.opp(qem.next(cir, 3));
389		✗	} while (cir != border[next_mod(e, border.size())]);
390			}
391		✗	FOR(e, border.size())
392		✗	FOR(v, opp_fan[e].size())
393		✗	qem.set_vertex(opp_fan[e][v], lower_v[next_mod(e, border.size())]);
394
395		✗	if (pts.size() > border.size()) {
396		✗	index_t off_v = m->vertices.create_vertices(2 * (pts.size() - border.size()));
397		✗	FOR(i, pts.size() - border.size()) {
398		✗	FOR(d,2) X(m)[off_v + 2 * i + d] = pts[border.size() + i];
399
400		✗	upper_v.push_back(off_v + 2 * i);
401		✗	lower_v.push_back(off_v + 2 * i + 1);
402			}
403			}
404		✗	FOR(q, quads.size() / 4) {
405		✗	m->facets.create_quad(
406		✗	upper_v[quads[4 * q + 0]],
407		✗	upper_v[quads[4 * q + 3]],
408		✗	upper_v[quads[4 * q + 2]],
409		✗	upper_v[quads[4 * q + 1]]
410			) ;
411		✗	m->facets.create_quad(
412		✗	lower_v[quads[4 * q + 0]],
413		✗	lower_v[quads[4 * q + 1]],
414		✗	lower_v[quads[4 * q + 2]],
415		✗	lower_v[quads[4 * q + 3]]
416			) ;
417			}
418			// debug output
419
420		✗	Attribute<int> date(m->edges.attributes(), "date");
421		✗	if (!border.empty()) {
422		✗	index_t off_e = m->edges.create_edges(border.size());
423		✗	FOR(e, border.size()) {
424		✗	date[off_e + e] = int(off_e);
425		✗	FOR(extr, 2)
426		✗	m->edges.set_vertex(off_e + e, extr, qem.vertex(border[(e + extr) % border.size()]));
427			}
428			}
429		✗	return true;
430
431		✗	}
432			};
433
434
435
436
437		✗	inline double cos_corner(vec3 B, vec3 A, vec3 C) {
438		✗	return dot(normalize(B - A), normalize(C - A));
439			}
440
441			struct VertexPuncher {
442			Mesh* m;
443			Mesh* newhex;
444			QuadExtraConnectivity qem;
445			DynamicHBoxes hb; // -> a static BBox tree
446			FacetIntersect finter;
447			double ave_edge_length;
448			Attribute<bool> dead_face;
449
450			int nb_punchs;
451			index_t punch_v;
452			index_t nv_punch_v;
453			index_t H[3][4];
454			index_t oppH[3][4];
455			vec3 old_vertex_position;
456			vec3 new_vertex_position;
457			vector<int> v2nb_facets; // -> facets that have moved
458			index_t via_facet[3];
459
460
461			Attribute<double> failt; ///DEBUG
462			Attribute<bool> isquad;
463
464
465		✗	VertexPuncher(Mesh* p_m, Mesh* p_newhex): finter(p_m) {
466		✗	m = p_m;
467		✗	newhex = p_newhex;
468		✗	dead_face.bind(m->facets.attributes(), "dead_face");
469		✗	}
470
471
472
473		✗	void unglue_duplicates() {
474		✗	FOR(lf, 3) {
475		✗	index_t f = via_facet[lf];
476		✗	if (f == NOT_AN_ID) continue;
477			index_t cir_h[4];
478			index_t cir_opp[4];
479		✗	index_t opp_f = NOT_AN_ID;
480		✗	FOR(lh, 4) {
481		✗	index_t h = 4 * f + lh;
482		✗	index_t h_opp = qem.opp(h);
483		✗	if (qem.vertex(qem.next(h_opp, 2)) == qem.vertex(qem.next(h, 3))
484		✗	&& qem.vertex(qem.next(h_opp, 3)) == qem.vertex(qem.next(h, 2))) {
485		✗	FOR(i, 4) {
486		✗	cir_h[i] = 4 * f + (lh + i) % 4;
487		✗	cir_opp[i] = 4 * qem.face(h_opp) + (qem.local_id(h_opp) + 4 - i) % 4;
488			}
489		✗	opp_f = qem.face(h_opp);
490		✗	continue;
491		✗	}
492			}
493		✗	if (opp_f == NOT_AN_ID) continue;
494
495		✗	dead_face[f] = true;
496		✗	dead_face[opp_f] = true;
497		✗	FOR(i, 4) {
498			//std::cerr << qem.vertex(cir_h[i]) << " " << qem.vertex(qem.next(cir_opp[i])) << " \n";
499		✗	if (qem.opp(cir_h[i]) == qem.opp(cir_opp[i])) continue;
500		✗	qem.set_opp(qem.opp(cir_opp[i]), qem.opp(cir_h[i]));
501		✗	qem.set_opp(cir_h[i], cir_opp[i]);
502			};
503			}
504		✗	}
505
506
507
508
509		✗	BBox facet_bbox(index_t f) {
510		✗	BBox res;
511		✗	FOR(fv, m->facets.nb_vertices(f)) res.add(X(m)[m->facets.vertex(f, fv)]);
512		✗	return res;
513			}
514
515		✗	void init() {
516		✗	isquad.bind(m->facets.attributes(), "isquad");
517
518		✗	qem.init(m);
519		✗	nb_punchs = 0;
520		✗	v2nb_facets.clear();
521			//moved_facets.clear();
522
523		✗	v2nb_facets.resize(m->vertices.nb(), 0);
524		✗	FOR(f, m->facets.nb()) FOR(v, m->facets.nb_vertices(f)) v2nb_facets[m->facets.vertex(f, v)]++;
525
526
527
528		✗	FOR(f, m->facets.nb())dead_face[f] = false;
529
530			// mesh resolution
531		✗	ave_edge_length = 0;
532		✗	int nb_samples = 0;
533		✗	FOR(h, 4 * m->facets.nb()) if (qem.valid(h)) {
534		✗	ave_edge_length += (X(m)[qem.vertex(h)] - X(m)[qem.vertex(qem.next(h))]).length();
535		✗	nb_samples++;
536			}
537		✗	ave_edge_length /= double(nb_samples);
538
539			// structures to find facets
540		✗	vector<BBox> inboxes(m->facets.nb());
541		✗	FOR(f, m->facets.nb()) inboxes[f] = facet_bbox(f);//FOR(fv, m->facets.nb_vertices(f)) inboxes[f].add(X(m)[m->facets.vertex(f, fv)]);
542		✗	hb.init(inboxes);
543
544
545		✗	}
546
547
548
549		✗	bool init_one_ring(index_t h) {
550			// check there is no triangles involved
551		✗	index_t cir = h;
552		✗	FOR(i, 3) {
553		✗	if (qem.fsize(cir) != 4) return false;
554		✗	cir = qem.next_around_vertex(cir);
555			}
556			// init H
557		✗	FOR(f, 3) {
558		✗	FOR(v, 4) {
559		✗	H[f][v] = qem.next(h, v);
560		✗	if (H[f][v] == NOT_AN_ID) return false;// TO REMOVE if m is closed
561			}
562			//if (!isquad[qem.face(H[f][0])]) return false;
563		✗	h = qem.next_around_vertex(h);
564			}
565
566			// check valence 3
567		✗	if (qem.next_around_vertex(H[2][0]) != H[0][0]) return false;
568
569			// init the opposites
570		✗	FOR(f, 3) FOR(e, 4) {
571		✗	oppH[f][e] = qem.opp(H[f][e]);
572		✗	if (oppH[f][e] == NOT_AN_ID) return false;// TO REMOVE if m is closed
573
574			}
575		✗	return true;
576			}
577			// create new hex
578		✗	void create_new_hex() {
579		✗	index_t off_v = newhex->vertices.create_vertices(8);
580		✗	Attribute<double> init(newhex->vertices.attributes(), "init");
581		✗	init[off_v] = -100; FOR(i, 7)init[off_v + i + 1] = nb_punchs;
582		✗	X(newhex)[off_v + 0] = old_vertex_position;
583		✗	X(newhex)[off_v + 1] = X(m)[qem.vertex(H[0][1])];
584		✗	X(newhex)[off_v + 2] = X(m)[qem.vertex(H[0][3])];
585		✗	X(newhex)[off_v + 3] = X(m)[qem.vertex(H[0][2])];
586		✗	X(newhex)[off_v + 4] = X(m)[qem.vertex(H[2][1])];
587		✗	X(newhex)[off_v + 5] = X(m)[qem.vertex(H[2][2])];
588		✗	X(newhex)[off_v + 6] = X(m)[qem.vertex(H[1][2])];
589		✗	X(newhex)[off_v + 7] = X(m)[nv_punch_v];
590		✗	newhex->cells.create_hex(off_v + 0, off_v + 1, off_v + 2, off_v + 3, off_v + 4, off_v + 5, off_v + 6, off_v + 7);
591		✗	}
592
593		✗	bool new_hex_geometry_is_crappy() {
594
595
596		✗	FOR(front, 2) {// front=0 for actual faces, front =1 for new faces
597		✗	FOR(fid, 3) {
598		✗	vector<vec3> v(4);
599		✗	if (front == 0) {
600		✗	v[0] = old_vertex_position;
601		✗	v[1] = X(m)[qem.vertex(H[fid][1])];
602		✗	v[2] = X(m)[qem.vertex(H[fid][2])];
603		✗	v[3] = X(m)[qem.vertex(H[fid][3])];
604			}
605			else {
606		✗	v[0] = X(m)[nv_punch_v];
607		✗	v[1] = X(m)[qem.vertex(H[next_mod(fid, 3)][2])];
608		✗	v[2] = X(m)[qem.vertex(H[next_mod(fid, 3)][1])];
609		✗	v[3] = X(m)[qem.vertex(H[fid][2])];
610			}
611		✗	vec3 n = Poly3d(v).normal();
612		✗	FOR(lv, 4) if (dot(n, cross(normalize(v[(lv + 2) % 4] - v[(lv + 1) % 4]), normalize(v[(lv) % 4] - v[(lv + 1) % 4]))) < .1)
613		✗	return true;
614
615		✗	}
616			}
617		✗	return false;
618
619			}
620
621			//topo_punch
622		✗	void topo_punch() {
623		✗	FOR(f, 3) qem.set_vertex(H[f][0], nv_punch_v);
624		✗	FOR(f, 3) qem.set_vertex(H[f][1], qem.vertex(oppH[(f + 1) % 3][1]));
625		✗	FOR(f, 3) qem.set_vertex(H[f][2], qem.vertex(oppH[(f + 1) % 3][2]));
626		✗	FOR(f, 3) qem.set_vertex(H[f][3], qem.vertex(oppH[(f + 2) % 3][1]));
627		✗	FOR(f, 3) qem.set_opp(H[f][1], oppH[(f + 1) % 3][2]);
628		✗	FOR(f, 3) qem.set_opp(H[f][2], oppH[(f + 2) % 3][1]);
629		✗	}
630		✗	bool topo_can_punch() {
631		✗	vector<index_t> neigh;
632			// check if a vertex is duplicated
633		✗	FOR(ring, 3) FOR(lv, 2) neigh.push_back(H[ring][lv + 1]);
634		✗	FOR(n0, 6) if (qem.opp(neigh[n0]) == NOT_AN_ID) return false;
635		✗	FOR(n0, 6)for (index_t n1 = 0; n1 < n0; n1++)
636		✗	if (qem.vertex(neigh[n0]) == qem.vertex(neigh[n1])) return false;;
637
638
639			// check is two boundary edges are opposite
640		✗	FOR(n0, 6)for (index_t n1 = 0; n1 < n0; n1++)
641		✗	if (qem.opp(neigh[n0]) == neigh[n1]) return false;
642		✗	return true;
643		✗	}
644
645
646			void produce_diamon(vector<vec3>& P, vec3 A, vec3 B, vec3 C, vec3 D, double h) {
647			P.reserve(6);
648			P.resize(4);
649			P[0] = A; P[1] = B; P[2] = C; P[3] = D;
650			vec3 G = Poly3d(P).barycenter();
651			vec3 n = Poly3d(P).normal();
652			P.push_back(G + h*n);
653			P.push_back(G - h*n);
654			}
655
656
657		✗	bool punch_will_produce_intersection() {
658			// check for geometric intersections
659			index_t Qv[3][4];
660		✗	FOR(f, 3) Qv[f][0] =nv_punch_v;
661		✗	FOR(f, 3) Qv[f][1] = qem.vertex(oppH[(f + 1) % 3][1]);
662		✗	FOR(f, 3) Qv[f][2] = qem.vertex(oppH[(f + 1) % 3][2]);
663		✗	FOR(f, 3) Qv[f][3] = qem.vertex(oppH[(f + 2) % 3][1]);
664
665
666		✗	FOR(f0, 3)FOR(f1, 3) {
667		✗	if (f0 >= f1) continue;
668		✗	vector<vec3> Q0(4), Q1(4);
669		✗	FOR(i, 4) Q0[i]= X(m)[Qv[f0][i]];
670		✗	FOR(i, 4) Q1[i] = X(m)[Qv[f1][i]];
671		✗	if (polyintersect(Q0, Q1)) return true;
672		✗	}
673		✗	FOR(fid, 3) {
674		✗	vector<vec3> Q;
675		✗	FOR(i, 4) Q.push_back(X(m)[Qv[fid][i]]);
676		✗	if (finter.get_intersections(Q).size() > 0) return true;
677		✗	continue;
678		✗	}
679		✗	return false;
680			}
681
682
683		✗	bool apply(int& nbmaxpunch) {
684		✗	if (m->facets.nb() == 0) return false;
685		✗	init();
686		✗	bool finished = false;
687		✗	failt.bind(m->vertices.attributes(), "failt");
688		✗	FOR(v, m->vertices.nb()) failt[v] = 0;
689
690
691
692		✗	while (!finished) {
693
694		✗	bool found_vertex_to_punch = false;
695		✗	FOR(seed, 4 * m->facets.nb()) {
696		✗	if (!qem.valid(seed)) continue;
697		✗	if (qem.next_around_vertex(seed) < seed || qem.next_around_vertex(qem.next_around_vertex(seed)) < seed) continue;
698		✗	punch_v = qem.vertex(seed);
699		✗	if (!init_one_ring(seed)) { failt[punch_v] = std::max(failt[punch_v], 10.); continue; }
700		✗	if (!topo_can_punch()) { failt[punch_v] = std::max(failt[punch_v], 20.); continue; }
701		✗	old_vertex_position = X(m)[punch_v];
702		✗	nv_punch_v = punch_v;
703
704		✗	if (tet_vol(X(m)[punch_v],
705		✗	X(m)[qem.vertex(H[0][1])],
706		✗	X(m)[qem.vertex(H[1][1])],
707		✗	X(m)[qem.vertex(H[2][1])]
708		✗	) > 0) continue;
709
710			// do we already have 4+ faces of the hex ?
711		✗	bool bad_config = false;
712		✗	index_t punch_cand[3] = { NOT_AN_ID, NOT_AN_ID, NOT_AN_ID };
713		✗	FOR(i, 3) via_facet[i] = NOT_AN_ID;
714		✗	index_t punch_cand_ref = NOT_AN_ID;
715		✗	FOR(i, 3) if (qem.valence(oppH[i][2]) == 3) {
716		✗	punch_cand[i] = qem.vertex(qem.next(oppH[i][2], 2));
717		✗	punch_cand_ref = punch_cand[i];
718		✗	via_facet[i] = qem.face(oppH[i][2]);
719			}
720		✗	FOR(i, 3) {
721		✗	if (punch_cand[i] != NOT_AN_ID && punch_cand[i] != punch_cand_ref)
722		✗	bad_config = true;
723		✗	if (qem.valence(H[i][2]) == 3)
724		✗	if ((via_facet[i] == NOT_AN_ID) != (via_facet[(i + 2) % 3] == NOT_AN_ID))
725		✗	bad_config = true;
726		✗	if (via_facet[i] != NOT_AN_ID && m->facets.nb_vertices(via_facet[i])==3)
727		✗	bad_config = true;
728			}
729		✗	if (bad_config) { failt[punch_v] = std::max(failt[punch_v], 20.); continue; }
730
731		✗	if (punch_cand_ref != NOT_AN_ID) {
732		✗	nv_punch_v = punch_cand_ref;
733		✗	FOR(i, 3) if (via_facet[i] != NOT_AN_ID) dead_face[via_facet[i]] = true; // the face will be cancelled by its opposite
734			}
735
736
737
738		✗	if (nv_punch_v == punch_v) {
739			//if (search_for_existing_vertex_to_punch) continue;
740			// check geometry of existing faces
741		✗	bool have_bad_angle = false;
742		✗	FOR(ring, 3) FOR(lv, 4)
743		✗	have_bad_angle = have_bad_angle || std::abs(cos_corner(
744		✗	X(m)[qem.vertex(H[ring][lv])],
745		✗	X(m)[qem.vertex(H[ring][(lv + 1) % 4])],
746		✗	X(m)[qem.vertex(H[ring][(lv + 2) % 4])]
747			)) > .8;
748		✗	if (have_bad_angle) { failt[punch_v] = std::max(failt[punch_v], 30.); continue; }
749
750			//-------------------------
751			{
752		✗	vec3 cubebary(0, 0, 0);
753		✗	FOR(i, 3) FOR(e, 2) cubebary = cubebary + X(m)[qem.vertex(H[i][1 + e])];
754		✗	cubebary = (1. / 6.) *cubebary;
755		✗	vec3 decal = cubebary - old_vertex_position;
756		✗	vec3 n[3];
757		✗	FOR(i, 3) n[i] = facet_normal(m, qem.face(H[i][0]));
758		✗	if (dot(decal, n[0] + n[1] + n[2]) > 0) {
759		✗	failt[punch_v] = std::max(failt[punch_v], 40.);
760		✗	continue;
761			}
762		✗	new_vertex_position = 2. * cubebary - old_vertex_position;
763			}
764			// new way to compute new position
765		✗	vec3 n[3];
766		✗	mat3 mat;
767			index_t tri[3][3];// 3 triplet of vertices that miss the last point
768		✗	FOR(f, 3) {
769		✗	index_t next_f = next_mod(f, 3);
770		✗	tri[f][0] = qem.vertex(H[f][2]);
771		✗	tri[f][1] = qem.vertex(H[next_f][1]);
772		✗	tri[f][2] = qem.vertex(H[next_f][2]);
773			}
774
775		✗	FOR(f, 3) {
776		✗	have_bad_angle = have_bad_angle
777		✗	|| std::abs(cos_corner(X(m)[tri[f][0]], X(m)[tri[f][1]], X(m)[tri[f][2]])) > cos(M_PI / 8.);
778			}
779		✗	if (have_bad_angle) { failt[punch_v] = std::max(failt[punch_v], 50.); continue; }
780
781		✗	FOR(f, 3) {
782		✗	n[f] = normalize(cross(X(m)[tri[f][1]] - X(m)[tri[f][0]], X(m)[tri[f][2]] - X(m)[tri[f][0]]));
783		✗	FOR(j, 3) mat(f, j) = n[f][j];
784			}
785		✗	mat3 inv = mat.inverse();
786		✗	vec3 b;
787		✗	FOR(i, 3) b[i] = dot(n[i], X(m)[tri[i][0]]);
788		✗	mult(inv, b.data(), new_vertex_position.data());
789		✗	X(m)[nv_punch_v] = new_vertex_position;
790			}
791
792			// check that punch vertex is convex
793		✗	if (new_hex_geometry_is_crappy()) {
794		✗	X(m)[punch_v] = old_vertex_position;
795		✗	failt[punch_v] = std::max(failt[punch_v], 60.);
796		✗	continue;
797			}
798
799
800
801		✗	if (punch_will_produce_intersection()) {
802		✗	X(m)[punch_v] = old_vertex_position;
803		✗	continue;
804			}
805
806			// split vertex if needed (non manifold)
807		✗	if (punch_v == nv_punch_v && v2nb_facets[punch_v] != 3) {
808		✗	index_t nvv = m->vertices.create_vertex();
809		✗	v2nb_facets[punch_v] -= 3;
810		✗	v2nb_facets.push_back(3);
811		✗	X(m)[punch_v] = old_vertex_position;
812
813		✗	X(m)[nvv] = new_vertex_position;
814		✗	punch_v = nvv;
815		✗	nv_punch_v = nvv;
816			}
817
818
819		✗	FOR(f, 3) { v2nb_facets[qem.vertex(H[f][1])]--; v2nb_facets[qem.vertex(H[f][2])]++; }
820		✗	plop("create_new_hex()");
821		✗	create_new_hex();
822		✗	topo_punch();
823		✗	FOR(lf, 3) hb.update_bbox(qem.face(H[lf][0]), facet_bbox(qem.face(H[lf][0])));
824		✗	FOR(lf, 3) finter.hb.update_bbox(qem.face(H[lf][0]), facet_bbox(qem.face(H[lf][0])));
825		✗	found_vertex_to_punch = true;
826
827		✗	unglue_duplicates();
828
829		✗	nb_punchs++;
830		✗	if (!(nb_punchs % 100)) plop(nb_punchs);
831		✗	if (-1 == --nbmaxpunch) {
832		✗	plop(nbmaxpunch);
833		✗	goto cleanup;// return true;
834			}// debug only... to be removed
835			}
836		✗	plop("done");
837		✗	finished = !found_vertex_to_punch;
838
839			}
840
841		✗	cleanup:
842
843		✗	qem.debug_export_adjacence();
844		✗	vector<index_t> to_kill(m->facets.nb());
845		✗	FOR(f, m->facets.nb()) to_kill[f]= (m->facets.nb_vertices(f) == 4) ;
846		✗	FOR(f, m->facets.nb()) {
847		✗	index_t opp = qem.face(qem.opp(4 * f));
848		✗	if (m->facets.nb_vertices(f) != 4) continue;
849		✗	if (m->facets.nb_vertices(opp) != 4) { to_kill[f] = false; continue; }
850		✗	for (index_t h = 4 * f; h < 4 * (f + 1); h++)
851		✗	if (qem.face(qem.opp(h)) != opp) {
852		✗	to_kill[f] = false;
853		✗	to_kill[opp] = false;
854			}
855			}
856		✗	m->facets.delete_elements(to_kill);
857			//check_no_intersecting_faces(m, true);
858		✗	plop(nb_punchs);
859		✗	return nb_punchs > 0;
860			}
861
862
863			};
864
865
866
867
868
869
870		✗	static void remove_scabs(Mesh* m, Mesh* newhex) {
871		✗	GEO::Logger::out("HexDom") << "try to remove_scabs" << std::endl;
872
873		✗	Attribute<int> ft(m->facets.attributes(), "ft"); // facet type
874		✗	vector<index_t> to_kill(m->facets.nb(), false); // ;( faces may be hurt and even killed in this fonction
875		✗	vector<index_t> local_id(m->vertices.nb(), NOT_AN_ID); // ids in the array of new vertices (projected onto the boundary)
876
877
878		✗	QuadExtraConnectivity qem;
879		✗	qem.init(m);
880		✗	FOR(f, m->facets.nb()) ft[f] = 0;
881		✗	FOR(h, 4 * m->facets.nb()) {
882
883		✗	if (!qem.valid(h)) continue;
884		✗	if (to_kill[qem.face(h)]) continue;
885		✗	if (qem.opp(h) == NOT_AN_ID) continue;
886		✗	if (ft[qem.face(qem.opp(h))] != 0) continue;
887
888		✗	vector<index_t> contour;
889		✗	bool fail = false;
890		✗	{index_t cir = h;
891			do {
892		✗	contour.push_back(cir);
893		✗	if (qem.opp(cir) == NOT_AN_ID) { fail = true; break; }
894		✗	if (qem.opp(qem.next(cir, 3)) == NOT_AN_ID) { fail = true; break; }
895		✗	if (qem.fsize(cir) != 4) { fail = true; break; }
896		✗	if (qem.fsize(qem.opp(qem.next(cir, 1))) != 4) { fail = true; break; }
897		✗	if (qem.fsize(qem.opp(qem.next(cir, 3))) == 4) { fail = true; break; }
898		✗	cir = qem.next(qem.opp(cir), 2);
899		✗	} while (cir != h);
900			}
901		✗	if (fail) continue;
902
903
904
905
906		✗	FOR(c, contour.size()) ft[qem.face(contour[c])] = 1;
907
908		✗	vector<index_t> in_facets;
909		✗	vector<index_t> out_facets;
910		✗	in_facets.push_back(qem.face(qem.opp(qem.next(h))));
911		✗	out_facets.push_back(qem.face(qem.opp(qem.next(h, 3))));
912
913		✗	ft[in_facets[0]] = 2;
914		✗	ft[out_facets[0]] = 3;
915
916		✗	vector<index_t> stack;
917		✗	stack.push_back(in_facets[0]);
918		✗	stack.push_back(out_facets[0]);
919		✗	while (!stack.empty() && !fail) {
920		✗	index_t f = stack.back();
921		✗	stack.pop_back();
922		✗	for (index_t cir = 4 * f; cir < 4 * f + m->facets.nb_vertices(f); cir++) {
923			// need to have an opposite
924		✗	if (qem.opp(cir) == NOT_AN_ID) {
925		✗	fail = true;
926		✗	break;
927			}
928			// do not link triangle/quads OR quads with the contour
929		✗	index_t oppf = qem.face(qem.opp(cir));
930		✗	if (ft[oppf] != 1 && m->facets.nb_vertices(oppf) != m->facets.nb_vertices(f)) {
931		✗	fail = true;
932		✗	break;
933			}
934
935
936
937
938		✗	if (ft[oppf] == 0) {
939			// linked by a reasonably flat edge---only for inside until outside have better mesh quality
940		✗	if (dot(facet_normal(m, f), facet_normal(m, oppf)) < .8 && ft[f] == 2) {
941		✗	fail = true;
942		✗	break;
943			}
944		✗	ft[oppf] = ft[f];
945		✗	if (ft[f] == 2) in_facets.push_back(oppf);
946		✗	if (ft[f] == 3) out_facets.push_back(oppf);
947		✗	stack.push_back(oppf);
948			}
949			}
950			}
951		✗	if (fail) {
952		✗	FOR(f, m->facets.nb()) ft[f] = 0;
953		✗	continue;
954		✗	}
955			// check that in_facets and outfacets are topo disks
956		✗	index_t in_facets_nb_neig = 0;
957		✗	FOR(i, in_facets.size()) {
958		✗	nico_assert(m->facets.nb_vertices(in_facets[i]) == 4);
959		✗	for (index_t cir = 4 * in_facets[i]; cir < 4 * in_facets[i] + 4; cir++) {
960		✗	if (ft[qem.face(qem.opp(cir))] != 2)
961		✗	in_facets_nb_neig++;
962			}
963			}
964		✗	if (in_facets_nb_neig != contour.size()) {
965		✗	FOR(f, m->facets.nb()) ft[f] = 0;
966		✗	continue;
967		✗	}
968		✗	index_t out_facets_nb_neig = 0;
969		✗	FOR(i, out_facets.size()) {
970		✗	nico_assert(m->facets.nb_vertices(out_facets[i]) != 4);
971		✗	for (index_t cir = 4 * out_facets[i]; cir < 4 * out_facets[i] + 3; cir++) {
972		✗	if (ft[qem.face(qem.opp(cir))] != 3)
973		✗	out_facets_nb_neig++;
974			}
975			}
976		✗	if (out_facets_nb_neig != contour.size()) {
977		✗	FOR(f, m->facets.nb()) ft[f] = 0;
978		✗	continue;
979		✗	}
980
981
982			//generate new vertices
983		✗	index_t nbv = 0;
984		✗	vector<vec3> packed_v; // for each vertex of infacets, stores the 3-uplet (pos, normal, projected pos)
985		✗	FOR(i, in_facets.size()) {
986		✗	for (index_t cir = 4 * in_facets[i]; cir < 4 * in_facets[i] + 4; cir++) {
987		✗	index_t v = qem.vertex(cir);
988		✗	if (local_id[v] == NOT_AN_ID) {
989		✗	local_id[v] = nbv;
990		✗	nbv++;
991		✗	packed_v.push_back(X(m)[v]);
992		✗	packed_v.push_back(vec3(0, 0, 0));
993		✗	packed_v.push_back(vec3(0, 0, 0));
994			}
995			// accumulate facet normals into vertex normal
996		✗	packed_v[3 * local_id[v] + 1] = packed_v[3 * local_id[v] + 1] + facet_normal(m, in_facets[i]);
997			}
998			}
999
1000
1001
1002			//return ;
1003
1004		✗	double ave_decal_length = 0;
1005		✗	int nb_intersections = 0;
1006		✗	FOR(v, nbv) {
1007		✗	vec3 O = packed_v[3 * v];
1008		✗	vec3 O2 = packed_v[3 * v] + normalize(packed_v[3 * v + 1]);
1009		✗	FOR(i, out_facets.size()) {
1010			//FOR(tr, 2)
1011			{
1012			//index_t cir = 4 * out_facets[i] + 2 * tr; // each quad is decomposed into 2 triangles
1013		✗	vec3 P[3];
1014		✗	FOR(p, 3) P[p] = X(m)[m->facets.vertex(out_facets[i], p)];// qem.vertex(qem.next(cir, p))];
1015			double sign[3];
1016		✗	FOR(p, 3) sign[p] = tetra_volume_sign(P[p], P[(p + 1) % 3], O, O2);
1017		✗	if (!same_sign(sign[0], sign[1]) || !same_sign(sign[0], sign[2])) continue;
1018		✗	double c0 = Geom::tetra_volume(P[0], P[1], P[2], O);
1019		✗	double c1 = Geom::tetra_volume(P[0], P[1], P[2], O2);
1020		✗	plop("intersection found");
1021		✗	packed_v[3 * v + 2] = O + (c0 / (c0 - c1))*(O2 - O);
1022		✗	ave_decal_length += (packed_v[3 * v + 2] - O).length();
1023		✗	nb_intersections++;
1024		✗	break;
1025			}
1026
1027			// if (packed_v[3 * v + 2].length2() != 0) break; //[Bruno: never executed, there is the "break" before].
1028			}
1029			}
1030		✗	ave_decal_length /= double(nb_intersections);
1031
1032		✗	FOR(v, nbv) {
1033		✗	if (packed_v[3 * v + 2].length2() == 0) {
1034		✗	plop("panic mode no intersection found");
1035		✗	packed_v[3 * v + 2] = packed_v[3 * v + 0] - ave_decal_length * normalize(packed_v[3 * v + 1]);
1036			}
1037			}
1038
1039
1040
1041			// vertices on border must match the contour
1042		✗	FOR(c, contour.size()) packed_v[3 * local_id[qem.vertex(qem.next(contour[c]))] + 2] = X(m)[qem.vertex(contour[c])];
1043
1044
1045			// generate hexes
1046		✗	index_t off_c = newhex->cells.create_hexes(in_facets.size());
1047		✗	index_t off_v = newhex->vertices.create_vertices(2 * nbv);
1048		✗	FOR(lv, nbv) {
1049		✗	X(newhex)[off_v + 2 * lv] = packed_v[3 * lv];
1050		✗	X(newhex)[off_v + 2 * lv + 1] = packed_v[3 * lv + 2];
1051			}
1052		✗	FOR(i, in_facets.size()) {
1053		✗	index_t f = in_facets[i];
1054		✗	newhex->cells.set_vertex(off_c + i, 0, off_v + 2 * local_id[qem.vertex(4 * f)]);
1055		✗	newhex->cells.set_vertex(off_c + i, 1, off_v + 2 * local_id[qem.vertex(4 * f + 1)]);
1056		✗	newhex->cells.set_vertex(off_c + i, 2, off_v + 2 * local_id[qem.vertex(4 * f + 3)]);
1057		✗	newhex->cells.set_vertex(off_c + i, 3, off_v + 2 * local_id[qem.vertex(4 * f + 2)]);
1058		✗	newhex->cells.set_vertex(off_c + i, 4, off_v + 1 + 2 * local_id[qem.vertex(4 * f)]);
1059		✗	newhex->cells.set_vertex(off_c + i, 5, off_v + 1 + 2 * local_id[qem.vertex(4 * f + 1)]);
1060		✗	newhex->cells.set_vertex(off_c + i, 6, off_v + 1 + 2 * local_id[qem.vertex(4 * f + 3)]);
1061		✗	newhex->cells.set_vertex(off_c + i, 7, off_v + 1 + 2 * local_id[qem.vertex(4 * f + 2)]);
1062			}
1063		✗	plop("gna");
1064		✗	FOR(c, contour.size()) to_kill[qem.face(contour[c])] = true;
1065		✗	FOR(i, in_facets.size()) to_kill[in_facets[i]] = true;
1066		✗	FOR(i, out_facets.size()) to_kill[out_facets[i]] = true;
1067		✗	FOR(v, m->vertices.nb()) local_id[v] = NOT_AN_ID; //not clearly needed, but brainless
1068
1069
1070		✗	}
1071		✗	m->facets.delete_elements(to_kill);
1072		✗	}
1073
1074
1075			/* [BL unused]
1076			static void evaluate_edges_valence(Mesh* m) {
1077			Attribute<int> edge_angle(m->facet_corners.attributes(), "edgeangle");
1078			FOR(f, m->facets.nb()) FOR(h, 4) edge_angle[m->facets.corner(f, h)] = rand() % 3;
1079
1080			}
1081			*/
1082
1083		✗	static bool next_crunch(Mesh* m, Mesh* newhex, int nb_max_punch) {
1084		✗	GEO::Logger::out("HexDom") << "next_crunch" << std::endl;
1085			//newhex->clear();
1086
1087		✗	if (m->facets.nb() == 0) return false;
1088			//static index_t nb_splits = 0;
1089
1090			static index_t iter = index_t(-1);
1091
1092		✗	if (iter == index_t(-1)) {
1093		✗	iter = 0;
1094			}
1095			else {
1096		✗	iter++;
1097			}
1098
1099			//if (iter == 0) evaluate_edges_valence(m);
1100
1101
1102		✗	plop(iter);
1103
1104		✗	bool did_something = false;
1105		✗	while (nb_max_punch > 0) {
1106
1107		✗	plop(nb_max_punch);
1108		✗	GEO::Logger::out("HexDom") << "Try to Punch" << std::endl;
1109		✗	VertexPuncher punch(m, newhex);
1110		✗	if (punch.apply(nb_max_punch)) {
1111			//check_no_intersecting_faces(m,true);
1112		✗	did_something = true;
1113		✗	if (nb_max_punch <= 0) return true;
1114			}
1115			else {
1116			//plop(did_something);
1117		✗	if (did_something) return true;
1118		✗	break;
1119			}
1120		✗	}
1121
1122
1123
1124		✗	GEO::Logger::out("HexDom") << "Try to cut" << std::endl;
1125			static int cutit = 0;
1126		✗	CutSingularity cut(m);
1127		✗	if (cut.apply()) {
1128		✗	GEO::Logger::out("HexDom") << "------------------------" << std::endl;
1129		✗	GEO::Logger::out("HexDom") << "CUT success... doing a small laplacian smoothing to separate collocated vertices" << std::endl;
1130		✗	plop(cutit);
1131		✗	cutit++;
1132			//return false; // DEBUG TO REMOVE
1133			//Attribute<vec3> real_geometry(m->vertices.attributes(), "real_geometry");
1134			//FOR(v, m->vertices.nb()) geo_assert((real_geometry[v] - X(m)[v]).length2()<1e-15);
1135			//return false;
1136			//check_no_intersecting_faces(m,true);
1137		✗	return true;
1138			}
1139
1140		✗	Attribute<vec3> real_geometry(m->vertices.attributes(), "real_geometry");
1141		✗	FOR(v, m->vertices.nb()) real_geometry[v] = X(m)[v];
1142
1143		✗	GEO::Logger::out("HexDom") << "------------------------" << std::endl;
1144		✗	GEO::Logger::out("HexDom") << "CUT FAILED" << std::endl;
1145		✗	remove_scabs(m, newhex);
1146		✗	return false;
1147		✗	}
1148
1149
1150
1151		✗	void punch_and_cut(Mesh* m, Mesh* newhex, int nb_iter, int nb_punch_per_iter, bool check_validity) {
1152		✗	FOR(i, nb_iter) {
1153		✗	m->edges.clear();
1154		✗	GEO::Logger::out("HexDom") << "iteration " << i << " / " << nb_iter << " with " << nb_punch_per_iter << " punch per iter" << std::endl;
1155		✗	if (!next_crunch(m, newhex, nb_punch_per_iter))break;
1156		✗	if (check_validity) {
1157		✗	m->edges.clear();
1158			}
1159			}
1160		✗	}
1161
1162
1163
1164		✗	void bourrin_quadrangulate_facets(Mesh* m) {
1165		✗	if (m->facets.nb() == 0) return;
1166		✗	Attribute<bool> isquad(m->facets.attributes(), "isquad");
1167		✗	index_t nb_facets = m->facets.nb();
1168		✗	FOR(f, nb_facets) {
1169		✗	index_t nbv = m->facets.nb_vertices(f);
1170		✗	vec3 bary(0, 0, 0);
1171		✗	FOR(fv, nbv) bary = bary + (1. / double(nbv))*X(m)[m->facets.vertex(f, fv)];
1172		✗	index_t nb_v = m->facets.nb_corners(f);
1173		✗	FOR(lc, nb_v) {
1174			vec3 v[3] = {
1175		✗	X(m)[m->facets.vertex(f, prev_mod(lc, nb_v))],
1176		✗	X(m)[m->facets.vertex(f, lc)],
1177		✗	X(m)[m->facets.vertex(f, next_mod(lc, nb_v))]
1178		✗	};
1179		✗	index_t off_v = m->vertices.create_vertices(4);
1180		✗	X(m)[off_v] = bary;
1181		✗	X(m)[off_v + 1] = 0.5*(v[0] + v[1]);
1182		✗	X(m)[off_v + 2] = v[1];
1183		✗	X(m)[off_v + 3] = 0.5*(v[1] + v[2]);
1184		✗	isquad[m->facets.create_quad(off_v + 0, off_v + 1, off_v + 2, off_v + 3)] = (nb_v == 4);
1185			}
1186			}
1187			{
1188		✗	vector<index_t> to_kill(m->facets.nb(), false);
1189		✗	FOR(f, nb_facets) to_kill[f] = true;
1190		✗	m->facets.delete_elements(to_kill);
1191		✗	}
1192
1193			// merge vertices
1194			{
1195		✗	vector<index_t> to_kill(m->vertices.nb(), 0);
1196		✗	vector<index_t> old2new(m->vertices.nb());
1197		✗	Geom::colocate(m->vertices.point_ptr(0), 3, m->vertices.nb(), old2new, 1e-15);
1198		✗	FOR(f, m->facets.nb()) FOR(fv, 4) m->facets.set_vertex(f, fv, old2new[m->facets.vertex(f, fv)]);
1199		✗	FOR(v, m->vertices.nb()) if (old2new[v] != v) to_kill[v] = NOT_AN_ID;
1200		✗	m->vertices.delete_elements(to_kill);
1201		✗	}
1202
1203
1204		✗	}
1205
1206		✗	void bourrin_subdivide_hexes(Mesh* m) {
1207		✗	if (m->cells.nb() == 0) return;
1208		✗	index_t nb_cells = m->cells.nb();
1209		✗	index_t off_v = m->vertices.create_vertices(64 * nb_cells);
1210		✗	index_t off_c = m->cells.create_hexes(8 * nb_cells);
1211		✗	FOR(c, nb_cells) FOR(nc, 8)FOR(nv, 8)
1212		✗	m->cells.set_vertex(off_c + 8 * c + nc, nv, off_v + 64 * c + 8 * nc + nv);
1213		✗	FOR(c, nb_cells) {
1214		✗	vec3 pts[8];
1215		✗	FOR(cv, 8) pts[cv] = X(m)[m->cells.vertex(c, cv)];
1216
1217		✗	FOR(nci, 2)FOR(ncj, 2)FOR(nck, 2) { // new cell position in cell c
1218
1219		✗	FOR(nvi, 2)FOR(nvj, 2)FOR(nvk, 2) { // new vertex position in new cell
1220		✗	index_t nc = nci + 2 * ncj + 4 * nck;
1221		✗	index_t nv = nvi + 2 * nvj + 4 * nvk;
1222		✗	vec3 coeff[2];
1223		✗	coeff[0] = vec3(double(nci + nvi) / 2., double(ncj + nvj) / 2., double(nck + nvk) / 2.);
1224		✗	coeff[1] = vec3(1, 1, 1) - coeff[0];
1225		✗	vec3 pos = vec3(0, 0, 0);
1226		✗	FOR(i, 2)FOR(j, 2)FOR(k, 2) {
1227		✗	pos += coeff[i][0] * coeff[j][1] * coeff[k][2] * pts[i + 2 * j + 4 * k];
1228			}
1229		✗	X(m)[off_v + 64 * c + 8 * nc + nv] = pos;
1230			}
1231
1232			}
1233			}
1234
1235
1236		✗	vector<index_t> to_kill(m->cells.nb(), false);
1237		✗	FOR(c, nb_cells) to_kill[c] = true;
1238		✗	m->cells.delete_elements(to_kill);
1239
1240		✗	mesh_repair(*m, MESH_REPAIR_COLOCATE, 1e-15);
1241
1242		✗	}
1243
1244			struct Polyline {
1245		✗	void compute_param() {
1246		✗	dist_to_org.resize(P.size());
1247		✗	FOR(i, P.size()) {
1248		✗	if (i == 0) {
1249		✗	dist_to_org[0] = 0;
1250			} else {
1251		✗	dist_to_org[i] = dist_to_org[i - 1] + (P[i - 1] - P[i]).length();
1252			}
1253			}
1254		✗	}
1255		✗	double length() {
1256		✗	geo_assert(!P.empty());
1257		✗	if (dist_to_org.size() != P.size()) {
1258		✗	compute_param();
1259			}
1260		✗	return dist_to_org.back();
1261			}
1262		✗	vec3 interpolate(double prop) {
1263		✗	double d = prop*length();
1264		✗	FOR(i, P.size()-1) {
1265		✗	if (dist_to_org[i + 1] - dist_to_org[i] < 1e-8) {
1266		✗	return P[i];
1267			}
1268		✗	double c = (d - dist_to_org[i]) / (dist_to_org[i + 1] - dist_to_org[i]);
1269		✗	if(c >= 0 && c <= 1.0) {
1270		✗	return c*P[i] + (1. - c)*P[i + 1];
1271			}
1272			}
1273		✗	return P.back();
1274			}
1275			vector<vec3> P;
1276			vector<double> dist_to_org;
1277			};
1278
1279
1280
1281
1282
1283		✗	void quadrangulate_easy_boundary(Mesh* m) {
1284		✗	plop("quadrangulate_easy_boundary");
1285		✗	QuadExtraConnectivity qem;
1286		✗	qem.init(m);
1287		✗	double ave_edge_length = get_facet_average_edge_size(m);
1288
1289
1290			//mark charts
1291		✗	index_t nb_charts = 0;
1292		✗	Attribute<index_t> chart(m->facets.attributes(), "chart");
1293		✗	FOR(seed, m->facets.nb()) chart[seed] = index_t(-1);
1294		✗	FOR(seed, m->facets.nb()) {
1295		✗	if (chart[seed] != index_t(-1)) continue;
1296		✗	if (m->facets.nb_vertices(seed)==4) continue;
1297		✗	chart[seed] = nb_charts;
1298		✗	vector<index_t> stack;
1299		✗	stack.push_back(seed);
1300		✗	while (!stack.empty()) {
1301		✗	index_t f = stack.back();
1302		✗	stack.pop_back();
1303		✗	FOR(e, 3) {
1304		✗	index_t h = 4 * f + e;
1305		✗	index_t fopp = qem.face (qem.opp(h));
1306		✗	if (qem.fsize(qem.opp(h)) == 4) continue;
1307		✗	if (chart[fopp] != index_t(-1)) continue;
1308		✗	stack.push_back(fopp);
1309		✗	chart[fopp] = nb_charts;
1310			}
1311			}
1312		✗	nb_charts++;
1313		✗	}
1314
1315		✗	vector<vector<int> > chart_to_subcharts(nb_charts);
1316		✗	vector<int> subchart_to_charts;
1317
1318			//mark sub charts
1319		✗	int nb_subcharts = 0;
1320		✗	Attribute<int> subchart(m->facets.attributes(), "subchart");
1321		✗	FOR(seed, m->facets.nb()) subchart[seed] = -1;
1322		✗	FOR(seed, m->facets.nb()) {
1323		✗	if (subchart[seed] != -1) continue;
1324		✗	if (m->facets.nb_vertices(seed) == 4) continue;
1325		✗	subchart[seed] = nb_subcharts;
1326		✗	chart_to_subcharts[chart[seed]] .push_back(subchart[seed]);
1327		✗	subchart_to_charts.push_back(int(chart[seed]));
1328
1329		✗	vector<index_t> stack;
1330		✗	stack.push_back(seed);
1331		✗	while (!stack.empty()) {
1332		✗	index_t f = stack.back();
1333		✗	stack.pop_back();
1334		✗	vec3 n_f = facet_normal(m, f);
1335		✗	FOR(e, 3) {
1336		✗	index_t h = 4 * f + e;
1337		✗	index_t fopp = qem.face(qem.opp(h));
1338		✗	vec3 n_fopp = facet_normal(m, fopp);
1339		✗	if (qem.fsize(qem.opp(h)) == 4) continue;
1340		✗	if (subchart[fopp] != -1) continue;
1341		✗	if (dot(n_f, n_fopp) < cos(M_PI / 4.)) continue;// hard edge detected
1342		✗	stack.push_back(fopp);
1343		✗	subchart[fopp] = nb_subcharts;
1344			}
1345			}
1346		✗	nb_subcharts++;
1347		✗	}
1348			// try to quadrangulate each subchart
1349		✗	vector<vector<vector<index_t> > > quads(nb_charts);
1350		✗	vector<vector<vector<vec3> > > pts(nb_charts);
1351		✗	FOR(c, nb_charts) quads[c].resize(chart_to_subcharts[c].size());
1352		✗	FOR(c, nb_charts) pts[c].resize(chart_to_subcharts[c].size());
1353
1354
1355		✗	Attribute<int> order(m->facet_corners.attributes(), "order");
1356		✗	FOR(sub, nb_subcharts) {
1357
1358		✗	vector<index_t> border;
1359			// gather halfedges in border
1360			{
1361		✗	FOR(h, 4 * m->facets.nb()) if (qem.valid(h))
1362		✗	if (subchart[qem.face(h)] == int(sub)
1363		✗	&& subchart[qem.face(qem.opp(h))] != int(sub))
1364		✗	border.push_back(h);
1365		✗	if (border.size() < 4) continue;
1366			}
1367			// order border halfedges
1368			{
1369			// do not start with an hardedge
1370		✗	FOR(cur, border.size()) if (qem.fsize(qem.opp(border[cur])) == 4) std::swap(border[0], border[cur]);
1371		✗	if (qem.fsize(qem.opp(border[0])) !=4) continue;
1372			//link them
1373		✗	for (index_t cur = 0; cur < border.size() - 1; cur++)
1374		✗	for (index_t next = cur + 1; next < border.size(); next++)
1375		✗	if (qem.vertex(qem.next(border[cur])) == qem.vertex(border[next]))
1376		✗	std::swap(border[cur + 1], border[next]);
1377			}
1378			// check that we have a loop
1379			{
1380		✗	bool topodisk = true;
1381		✗	for (index_t cur = 0; cur < border.size() - 1; cur++)
1382		✗	topodisk = topodisk && (qem.vertex(qem.next(border[cur])) == qem.vertex(border[cur + 1]));
1383		✗	if (!topodisk) { plop(topodisk); continue; }
1384			}
1385		✗	/*DEBUG*/FOR(cur,border.size()) order[qem.corner(border[cur])] = int(cur+10);
1386
1387			// constuct the problem.
1388		✗	int chartid = subchart_to_charts[sub];
1389		✗	vector<vec3> l_pts;
1390
1391		✗	vector<vec3> edge_org;
1392		✗	vector<bool> hardedge;
1393		✗	index_t offset = 0;
1394		✗	index_t N = border.size();
1395		✗	FOR(cur, N) {
1396		✗	edge_org.push_back(X(m)[qem.vertex(border[cur])]);
1397		✗	if (cur == 0) {
1398		✗	hardedge.push_back(false);
1399		✗	continue;
1400			}
1401		✗	hardedge.push_back(qem.fsize(qem.opp(border[cur])) != 4
1402		✗	&& subchart[qem.face(qem.opp(border[cur]))] == subchart[qem.face(qem.opp(border[cur - 1]))]);
1403		✗	if (!hardedge.back()) offset = cur;
1404			}
1405		✗	{vector<vec3> tmp(N); FOR(cur, N) tmp[cur] = edge_org[(cur + offset) % N]; edge_org = tmp; }
1406		✗	{vector<bool> tmp(N); FOR(cur, N) tmp[cur] = hardedge[(cur + offset) % N]; hardedge = tmp; }
1407
1408		✗	index_t i = 0;
1409		✗	while (i < N) {
1410			//plop(i);
1411		✗	l_pts.push_back(edge_org[i]);
1412		✗	if (!hardedge[i])
1413			{
1414		✗	i++;
1415			} else {
1416		✗	Polyline poly;
1417		✗	while (i < N+1 && hardedge[i%N]) {
1418
1419		✗	poly.P.push_back(edge_org[i%N]);
1420		✗	i++;
1421			}
1422		✗	plop(i);
1423		✗	plop(poly.length());
1424		✗	plop(ave_edge_length);
1425		✗	int nb_seg = int(2 * nint(0.5*poly.length() / ave_edge_length));
1426		✗	nb_seg = std::max(1, nb_seg);
1427		✗	FOR(s, nb_seg - 1) plop(poly.interpolate(double(s + 1) / double(nb_seg) - .001));
1428
1429		✗	}
1430			}
1431			//plop("border completed");
1432
1433		✗	vector<index_t> l_quads;
1434		✗	if (Poly3d(l_pts).try_quadrangulate(l_quads)) {
1435		✗	int loc_sub = -1;
1436		✗	FOR(l, chart_to_subcharts[chartid].size())
1437		✗	if (chart_to_subcharts[chartid][l] == int(sub))
1438		✗	loc_sub = int(l);
1439
1440		✗	FOR(p, l_pts.size()) plop(l_pts[p]);
1441		✗	FOR(p, l_quads.size()) plop(l_quads[p]);
1442
1443		✗	pts[chartid][loc_sub] = l_pts;
1444		✗	quads[chartid][loc_sub] = l_quads;
1445
1446			}
1447		✗	}
1448
1449			// replace charts that have been sucessfully remeshed
1450
1451		✗	vector<bool> chart_remesh_is_ok(nb_charts, true);
1452		✗	FOR(chartid, nb_charts) FOR(lsub, pts[chartid].size())
1453		✗	chart_remesh_is_ok[chartid] = chart_remesh_is_ok[chartid] && !pts[chartid][lsub].empty();
1454		✗	plop("gna");
1455
1456		✗	vector<index_t> to_kill(m->facets.nb(),false);
1457		✗	FOR(f, m->facets.nb()) if (chart[f]!=NOT_AN_ID) to_kill[f] = chart_remesh_is_ok[chart[f]];
1458		✗	plop("gna");
1459
1460		✗	Attribute<int> kill2(m->facets.attributes(), "kill2");
1461		✗	FOR(f, m->facets.nb()) kill2[f] = int(to_kill[f]);
1462
1463		✗	FOR(chartid, nb_charts) {
1464		✗	if (!chart_remesh_is_ok[chartid]) continue;
1465		✗	FOR(lsub, pts[chartid].size()) {
1466		✗	index_t off_v = m->vertices.create_vertices(pts[chartid][lsub].size());
1467		✗	FOR(p, pts[chartid][lsub].size())
1468		✗	X(m)[off_v+p] = pts[chartid][lsub][p];
1469		✗	index_t off_f = m->facets.create_quads(quads[chartid][lsub].size() / 4);
1470		✗	FOR(q, quads[chartid][lsub].size()/4)
1471		✗	FOR(lv,4) m->facets.set_vertex(off_f + q,lv,off_v+quads[chartid][lsub][4*q+lv]);
1472			}
1473			}
1474
1475
1476
1477		✗	to_kill.resize(m->facets.nb(),false);
1478		✗	m->facets.delete_elements(to_kill);
1479
1480			// merge vertices
1481		✗	vector<index_t> to_kill_v(m->vertices.nb(), 0);
1482		✗	vector<index_t> old2new(m->vertices.nb());
1483		✗	Geom::colocate(m->vertices.point_ptr(0), 3, m->vertices.nb(), old2new, 1e-15);
1484
1485
1486			//FOR(f, m->facets.nb()) FOR(lv, m->facets.nb_vertices(f))
1487			// if(old2new[m->facets.vertex(f, lv)]
1488			// == old2new[m->facets.vertex(f, (lv + 1) % m->facets.nb_vertices(f))])
1489			//{ Attribute<int> crush(m->facets.attributes(), "crush");
1490			//crush[f] = 1;
1491			//plop(f); plop(m->facets.vertex(f, lv)); return;
1492			//}
1493
1494
1495		✗	FOR(f, m->facets.nb()) FOR(fv, m->facets.nb_vertices(f)) m->facets.set_vertex(f, fv, old2new[m->facets.vertex(f, fv)]);
1496		✗	FOR(v, m->vertices.nb()) if (old2new[v] != v) to_kill_v[v] = NOT_AN_ID;
1497		✗	m->vertices.delete_elements(to_kill_v);
1498
1499		✗	}
1500
1501		✗	void prepare_crunch(Mesh* m, bool subdivide, bool revert) {
1502		✗	if (subdivide) bourrin_quadrangulate_facets(m);
1503		✗	if (revert) FOR(f, m->facets.nb()) {
1504		✗	index_t save_v = m->facets.vertex(f, 0);
1505		✗	m->facets.set_vertex(f, 0, m->facets.vertex(f, 2));
1506		✗	m->facets.set_vertex(f, 2, save_v);
1507			}
1508		✗	quadrangulate_easy_boundary(m);
1509		✗	}
1510		✗	void hex_crunch(Mesh* m, Mesh* hex) {
1511		✗	geo_argused(hex);
1512		✗	check_no_intersecting_faces(m, false);
1513		✗	prepare_crunch(m, false, true);
1514		✗	plop("prepare_crunch done");
1515		✗	plop(get_intersecting_faces(m).size());
1516		✗	return;
1517			}
1518
1519
1520			}
1521