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.h>
41			#include <exploragram/hexdom/PGP.h>
42			#include <exploragram/hexdom/basic.h>
43			#include <exploragram/hexdom/extra_connectivity.h>
44			#include <geogram/numerics/matrix_util.h>
45			#include <geogram/basic/permutation.h>
46			#include <algorithm>
47			#include <geogram/mesh/triangle_intersection.h>
48			#include <geogram/mesh/mesh_tetrahedralize.h>
49			#include <geogram/delaunay/delaunay.h>
50
51			#include <geogram/points/nn_search.h>
52			#include <geogram/points/colocate.h>
53			#include <queue>
54
55			#include <exploragram/hexdom/mesh_inspector.h>
56			#include <exploragram/hexdom/intersect_tools.h>
57			#include <exploragram/hexdom/polygon.h>
58			#include <exploragram/hexdom/time_log.h>
59			namespace GEO {
60
61
62
63		✗	static void snap_vertices_to_ref_vertices(Mesh* m, Mesh* ref, double eps = 1e-10) {
64		✗	if (ref->vertices.nb() == 0) return;
65		✗	NearestNeighborSearch_var NN = NearestNeighborSearch::create(3, "default");
66		✗	NN->set_points(ref->vertices.nb(), (double*)ref->vertices.point_ptr(0));
67
68		✗	FOR(v, m->vertices.nb()) {
69		✗	index_t v_ref = NN->get_nearest_neighbor((double*)&(X(m)[v]));
70		✗	if ((X(m)[v] - X(ref)[v_ref]).length() < eps)
71		✗	X(m)[v] = X(ref)[v_ref];
72			}
73			}
74
75
76
77
78
79
80		✗	void kill_intersecting_hexes(Mesh* hex) {
81		✗	hex->cells.connect();
82			index_t nb_intersecting_hex = 0;
83		✗	vector<index_t> to_kill(hex->cells.nb(), false);
84
85			// alternativeway of finding intersections
86
87		✗	Mesh quadset;
88		✗	Attribute<index_t> cell_id(quadset.facets.attributes(), "cellid");
89			quadset.vertices.create_vertices(hex->vertices.nb());
90		✗	FOR(v, hex->vertices.nb()) X(&quadset)[v] = X(hex)[v];
91		✗	quadset.facets.create_quads(6 * hex->cells.nb());
92		✗	FOR(f, 6 * hex->cells.nb()) cell_id[f] = f / 6;
93		✗	FOR(c, hex->cells.nb()) FOR(lf, 6) FOR(lv, 4) quadset.facets.set_vertex(6 * c+lf, lv, hex->cells.facet_vertex(c,lf,lv));
94
95			// kill adjacency faces
96		✗	vector<index_t> to_kill_f(quadset.facets.nb(), false);
97		✗	FOR(c, hex->cells.nb()) FOR(lf, 6) if (hex->cells.adjacent(c, lf) != NOT_AN_ID) to_kill_f[6 * c + lf] = true;
98		✗	quadset.facets.delete_elements(to_kill_f);
99
100			//plop("gre");
101
102
103			//mesh_save(quadset, "C:\\DATA\\quadeset.meshb");
104		✗	vector<index_t> intersections = get_intersecting_faces(&quadset);
105		✗	plop(intersections.size());
106			nb_intersecting_hex = intersections.size();
107		✗	FOR(i, intersections.size()) to_kill[cell_id[intersections[i] ]] = true;
108		✗	FOR(i, intersections.size()) plop(intersections[i]);
109
110
111
112			//vector<BBox> inboxes(6 * hex->cells.nb());
113			//FOR(f, inboxes.size())FOR(cfv, 4)
114			// inboxes[f].add(X(hex)[hex->cells.facet_vertex(f / 6, f % 6, cfv)]);
115			//HBoxes hb(inboxes);
116
117			//FOR(c, hex->cells.nb()) FOR(cf, 6) {
118			// if (hex->cells.adjacent(c, cf) != NOT_AN_ID) continue; // just check on boundary
119			// BBox b;
120			// vector<vec3> Q(4);
121			// FOR(cfv, 4) {
122			// Q[cfv] = X(hex)[hex->cells.facet_vertex(c, cf, cfv)];
123			// b.add(Q[cfv]);
124			// }
125			// vector<index_t> primitives;
126			// hb.intersect(b, primitives);
127
128			// FOR(i, primitives.size()) {
129			// index_t other_c = primitives[i] / 6;
130			// if (other_c == c) continue;
131			// if (to_kill[other_c]) continue;
132			// index_t other_cf = primitives[i] % 6;
133			// if (hex->cells.adjacent(other_c, other_cf) != NOT_AN_ID) continue; // just check on boundary
134			// vector<vec3> P(4);
135			// FOR(other_cfv, 4) P[other_cfv] = X(hex)[hex->cells.facet_vertex(other_c, other_cf, other_cfv)];
136
137			// // check if opposite
138			// bool is_opp = false;
139			// FOR(v, 4) {
140			// index_t it = 0;
141			// while (it < 4 && (Q[it] - P[(v + it) % 4]).length2() == 0) it++;
142			// is_opp = is_opp || (it == 4);
143			// }
144
145			// if (is_opp) continue;
146
147			// // check intersection
148			// vector<TriangleIsect> result;
149			// FOR(c0, 2)FOR(c1, 2) {
150			// if (triangles_intersections(
151			// P[quad_rand_split[c0][0]], P[quad_rand_split[c0][1]], P[quad_rand_split[c0][2]],
152			// Q[quad_rand_split[c1][0]], Q[quad_rand_split[c1][1]], Q[quad_rand_split[c1][2]],
153			// result
154			// )) {
155			// plop("auto intersection found");
156			// to_kill[c] = true;
157			// nb_intersecting_hex++;
158			// }
159			// }
160			// }
161			//}
162		✗	logt.add_value("nb_intersecting_hex", nb_intersecting_hex);
163		✗	plop(nb_intersecting_hex);
164		✗	hex->cells.delete_elements(to_kill);
165		✗	if (nb_intersecting_hex > 0) kill_intersecting_hexes(hex);
166		✗	}
167
168
169		✗	void hex_set_2_hex_mesh(Mesh* hex, Mesh* quadtri) {
170		✗	logt.add_value("gna", 3);
171		✗	if (hex->cells.nb() == 0) return;
172
173			// merge vertices
174		✗	double eps = (1e-3)*get_cell_average_edge_size(hex);
175			{
176		✗	vector<index_t> to_kill(hex->vertices.nb(), 0);
177			vector<index_t> old2new(hex->vertices.nb());
178		✗	Geom::colocate(hex->vertices.point_ptr(0), 3, hex->vertices.nb(), old2new, eps);
179		✗	FOR(c, hex->cells.nb()) FOR(cv, 8) hex->cells.set_vertex(c, cv, old2new[hex->cells.vertex(c, cv)]);
180		✗	FOR(v, hex->vertices.nb()) if (old2new[v] != v) to_kill[v] = NOT_AN_ID;
181		✗	hex->vertices.delete_elements(to_kill);
182			}
183		✗	snap_vertices_to_ref_vertices(hex, quadtri, eps);
184
185			// check that there is no duplicated hex
186			{
187			int nb_duplicated_hex = 0;
188		✗	vector<vec3> sumVpos(hex->cells.nb(), vec3(0, 0, 0));
189		✗	vector<index_t> to_kill(hex->cells.nb(), 0);
190		✗	FOR(c, hex->cells.nb()) FOR(lv, 8) sumVpos[c] = sumVpos[c] + hex->vertices.point(hex->cells.vertex(c, lv));
191			vector<index_t> bary_old2new(hex->cells.nb());
192		✗	Geom::colocate((double*)(sumVpos.data()), 3, hex->cells.nb(), bary_old2new, eps);
193		✗	FOR(c, hex->cells.nb()) {
194			//geo_assert(bary_old2new[c] == c);
195		✗	if (bary_old2new[c] != c) {
196		✗	to_kill[c] = NOT_AN_ID;
197		✗	nb_duplicated_hex++;
198			}
199			}
200		✗	logt.add_value("nb_duplicated_hex", nb_duplicated_hex);
201		✗	plop(nb_duplicated_hex);
202		✗	hex->cells.delete_elements(to_kill);
203		✗	plop(hex->cells.nb());
204
205			}
206
207
208			// remove bad shaped hex
209			{
210			int nb_bad_shaped_hex = 0;
211		✗	vector<index_t> to_kill(hex->cells.nb(), false);
212		✗	FOR(c, hex->cells.nb()) {
213		✗	FOR(cf, 6)FOR(cfv, 4) {
214		✗	vec3 C = X(hex)[hex->cells.facet_vertex(c, cf, prev_mod(cfv, 4))];
215		✗	vec3 A = X(hex)[hex->cells.facet_vertex(c, cf, cfv)];
216		✗	vec3 B = X(hex)[hex->cells.facet_vertex(c, cf, next_mod(cfv, 4))];
217		✗	if (std::abs(dot(normalize(B - A), normalize(C - A))) > .6) {
218		✗	to_kill[c] = true;
219		✗	nb_bad_shaped_hex++;
220			}
221			}
222			}
223
224		✗	logt.add_value("nb_bad_shaped_hex", nb_bad_shaped_hex);
225		✗	plop(nb_bad_shaped_hex);
226		✗	hex->cells.delete_elements(to_kill);
227		✗	plop(hex->cells.nb());
228
229			}
230
231			// remove hexes that are linked by 3 vertices
232			{
233			index_t nb_hex_linked_by_3_vertices = 0;
234		✗	vector<index_t> to_kill(hex->cells.nb(), false);
235			vector<vector<index_t> > v2hexface(hex->vertices.nb());
236		✗	FOR(c, hex->cells.nb()) FOR(cf, 6) FOR(cfv, 4) v2hexface[hex->cells.facet_vertex(c, cf, cfv)].push_back(6 * c + cf);
237		✗	FOR(v, hex->vertices.nb()) FOR(h0, v2hexface[v].size())FOR(h1, h0) {
238		✗	index_t hexf0 = v2hexface[v][h0];
239		✗	index_t hexf1 = v2hexface[v][h1];
240		✗	if (to_kill[hexf0 / 6] || to_kill[hexf1 / 6]) continue;
241			index_t f[2][4];
242		✗	FOR(lc, 4) f[0][lc] = hex->cells.facet_vertex(hexf0 / 6, hexf0 % 6, lc);
243		✗	FOR(lc, 4) f[1][lc] = hex->cells.facet_vertex(hexf1 / 6, hexf1 % 6, lc);
244			index_t i0 = index_t(-1), i1 = index_t(-1);
245		✗	FOR(lc, 4) if (f[0][lc] == v) i0 = lc;
246		✗	FOR(lc, 4) if (f[1][lc] == v) i1 = lc;
247		✗	geo_assert(i0 != index_t(-1) && i1 != index_t(-1));
248			int nb_shared = 0;
249		✗	FOR(lc, 4) if (f[0][(i0 + lc) % 4] == f[1][(i1 + 4 - lc) % 4])nb_shared++;
250		✗	if (nb_shared == 3) {
251		✗	to_kill[hexf0 / 6] = true;
252		✗	nb_hex_linked_by_3_vertices++;
253			}
254			}
255		✗	logt.add_value("nb_hex_linked_by_3_vertices", nb_hex_linked_by_3_vertices);
256		✗	plop(nb_hex_linked_by_3_vertices);
257		✗	hex->cells.delete_elements(to_kill);
258		✗	plop(hex->cells.nb());
259
260			}
261		✗	kill_intersecting_hexes(hex);
262
263			// remove hex that are incompatible with quadtri
264			int sum_nb_hex_incompatible_with_quadtri = 0;
265		✗	if (quadtri != nullptr) {
266
267			int nb_hex_incompatible_with_quadtri = 0;
268
269			// remove all hex that touch border
270		✗	vector<index_t> to_kill(hex->cells.nb(), true);
271
272			vector<BBox> inboxes(quadtri->facets.nb());
273		✗	FOR(f, quadtri->facets.nb())
274		✗	FOR(fv, quadtri->facets.nb_vertices(f))
275		✗	inboxes[f].add(X(quadtri)[quadtri->facets.vertex(f, fv)]);
276		✗	HBoxes hb(inboxes);
277
278
279			bool may_have_intersections = true;
280		✗	while (may_have_intersections) {
281			nb_hex_incompatible_with_quadtri = 0;
282			may_have_intersections = false;
283			// check for intersection
284		✗	hex->cells.connect();
285		✗	FOR(c, hex->cells.nb()) {
286		✗	to_kill[c] = false;
287		✗	FOR(cf, 6) {
288		✗	if (hex->cells.adjacent(c, cf) != NOT_AN_ID) continue;
289			// Q contains the face of the hex PLUS two extra vertices to make a diamon shape
290			vector<vec3> Q;
291		✗	Q.reserve(6);
292		✗	FOR(cfv, 4) Q.push_back(X(hex)[hex->cells.facet_vertex(c, cf, cfv)]);
293		✗	vec3 G = Poly3d(Q).barycenter();
294		✗	vec3 n = Poly3d(Q).normal();
295			double decal = 0;
296		✗	FOR(cfv, 4) decal += .2*.25*(Q[cfv] - Q[next_mod(cfv, 4)]).length(); // .2*ave edge length
297		✗	Q.push_back(G + decal*n);
298		✗	Q.push_back(G - decal*n);
299
300			BBox b;
301		✗	FOR(cfv, 6) b.add(Q[cfv]);
302
303			vector<index_t> primitives;
304		✗	hb.intersect(b, primitives);
305
306			bool have_intersection = false;
307			bool have_tri_quad_intersect = false;
308		✗	FOR(i, primitives.size()) {
309			bool is_quatri_face = false;
310		✗	index_t f = primitives[i];
311			vector<vec3> P;
312		✗	FOR(fv, quadtri->facets.nb_vertices(f))
313		✗	P.push_back(X(quadtri)[quadtri->facets.vertex(f, fv)]);
314		✗	geo_assert(P.size() < 5);
315
316			// check "is_quatri_face"
317		✗	if (P.size() == 4) {
318			int nb_match = 0;
319		✗	FOR(v, 4) {
320		✗	FOR(it, 4) if ((Q[it] - P[(v + it) % 4]).length2() == 0)nb_match++;
321			//index_t it = 0;
322			//while (it < 4 && (Q[it] - P[(v + it) % 4]).length2() == 0) it++;
323			//is_quatri_face = is_quatri_face || (it == 4);
324			}
325		✗	if (nb_match == 4) is_quatri_face = true;
326		✗	if (nb_match == 3) have_tri_quad_intersect = true;
327			}
328			else {
329		✗	FOR(v, 3) {
330			index_t it = 0;
331		✗	while (it < 3 && (Q[it] - P[(v + it) % 3]).length2() == 0) it++;
332		✗	have_tri_quad_intersect = have_tri_quad_intersect || (it == 3);
333			}
334			}
335
336
337			// check for non degenerated intersections
338		✗	if (!is_quatri_face) {
339			TriangleIsects result;
340		✗	FOR(diam, 12) {
341		✗	if (P.size() == 3)
342		✗	have_intersection = have_intersection || triangles_intersections(
343			P[0], P[1], P[2],
344		✗	Q[diamon_split[diam][0]], Q[diamon_split[diam][1]], Q[diamon_split[diam][2]], result
345			);
346			else {
347		✗	FOR(qu, 4)
348		✗	have_intersection = have_intersection || triangles_intersections(
349		✗	P[quad_split[qu][0]], P[quad_split[qu][1]], P[quad_split[qu][2]],
350		✗	Q[diamon_split[diam][0]], Q[diamon_split[diam][1]], Q[diamon_split[diam][2]], result
351			);
352			}
353
354			}
355			}
356
357			}
358
359		✗	if (have_intersection) to_kill[c] = true;
360		✗	if (have_tri_quad_intersect) to_kill[c] = true;
361		✗	if (have_intersection || have_tri_quad_intersect)
362		✗	nb_hex_incompatible_with_quadtri++;
363			}
364			}
365
366		✗	sum_nb_hex_incompatible_with_quadtri += nb_hex_incompatible_with_quadtri;
367
368		✗	plop(nb_hex_incompatible_with_quadtri);
369		✗	may_have_intersections = (nb_hex_incompatible_with_quadtri > 0);
370		✗	hex->cells.delete_elements(to_kill);
371		✗	plop(hex->cells.nb());
372			}
373		✗	logt.add_value("nb_hex_incompatible_with_quadtri", sum_nb_hex_incompatible_with_quadtri);
374		✗	}
375
376		✗	hex->cells.connect();
377		✗	hex->cells.compute_borders();
378
379			}
380
381
382
383
384
385			}
386