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_dominant.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 <exploragram/hexdom/hex_cruncher.h>
52
53
54			#include <geogram/points/nn_search.h>
55			#include <geogram/points/colocate.h>
56			#include <queue>
57
58			#include <exploragram/hexdom/mesh_inspector.h>
59			#include <exploragram/hexdom/intersect_tools.h>
60			#include <exploragram/hexdom/polygon.h>
61			#include <exploragram/hexdom/frame.h>
62			#define FPG_UNCERTAIN_VALUE 0
63			#include <geogram/numerics/predicates/orient3d.h>
64			namespace GEO {
65
66		✗	static void fill_quad_tri_surface(Mesh* m, bool with_pyramids) {
67		✗	if (m->vertices.nb() == 0) return;
68		✗	m->edges.clear();
69
70		✗	if (with_pyramids) {
71
72			vector<index_t> pyrindex;
73			// Remark: spliting quads in 2 may be better on boundary... but it is not clear for constrained boundary...
74		✗	vector<index_t> to_kill(m->facets.nb(), 0);
75
76			index_t init_nb_facets = m->facets.nb();
77		✗	FOR(f, init_nb_facets) {
78		✗	if (m->facets.nb_vertices(f) != 4) continue;
79
80		✗	index_t nvv = m->vertices.create_vertex();
81		✗	vec3 n = facet_normal(m, f);
82			double d = 0;
83		✗	FOR(e, 4) d += .25*(X(m)[m->facets.vertex(f, e)] - X(m)[m->facets.vertex(f, (e + 1) % 4)]).length();
84		✗	X(m)[nvv] = facet_bary(m, f) +.2*d*n;
85		✗	to_kill[f] = 1;
86
87		✗	index_t off_f = m->facets.create_triangles(4);
88		✗	FOR(e, 4) {
89		✗	to_kill.push_back(0);
90		✗	m->facets.set_vertex(off_f + e, 0, m->facets.vertex(f, e));
91		✗	m->facets.set_vertex(off_f + e, 1, m->facets.vertex(f, (e + 1) % 4));
92			m->facets.set_vertex(off_f + e, 2, nvv);
93			}
94		✗	FOR(e, 4) pyrindex.push_back(m->facets.vertex(f, e));
95			pyrindex.push_back(nvv);
96			}
97		✗	m->facets.delete_elements(to_kill, false);
98
99		✗	m->facets.triangulate();
100		✗	create_non_manifold_facet_adjacence(m);
101			try {
102			mesh_tetrahedralize(*m, false, true, 1.);
103		✗	index_t off_c = m->cells.create_pyramids(pyrindex.size() / 5);
104		✗	FOR(p, pyrindex.size() / 5) FOR(lv, 5)
105		✗	m->cells.set_vertex(off_c + p, lv, pyrindex[5 * p + lv]);
106			}
107		✗	catch (const Delaunay::InvalidInput& error_report) {
108		✗	FOR(i, error_report.invalid_facets.size())
109		✗	plop(error_report.invalid_facets[i]);
110		✗	}
111
112			}
113			else {// !with_pyramids
114		✗	m->facets.triangulate();
115		✗	create_non_manifold_facet_adjacence(m);
116			try {
117			mesh_tetrahedralize(*m, false, true, 1.);
118			}
119		✗	catch (const Delaunay::InvalidInput& error_report) {
120		✗	FOR(i, error_report.invalid_facets.size())
121		✗	plop(error_report.invalid_facets[i]);
122		✗	}
123
124			}
125			}
126
127		✗	void fill_cavity_with_tetgen(Mesh* input, Mesh* tri, bool with_pyramid) {
128		✗	tri->copy(*input, false);
129		✗	fill_quad_tri_surface(tri, with_pyramid);
130		✗	}
131
132
133		✗	void add_hexes_to_tetmesh(Mesh* hex, Mesh* tet_mesh) {
134		✗	if (hex->cells.nb() == 0) return;
135
136
137			index_t off_v = tet_mesh->vertices.create_vertices(hex->vertices.nb());
138		✗	FOR(v, hex->vertices.nb()) X(tet_mesh)[off_v + v] = X(hex)[v];
139		✗	index_t off_c = tet_mesh->cells.create_hexes(hex->cells.nb());
140		✗	FOR(c, hex->cells.nb())FOR(cv, 8) {
141		✗	index_t vid = hex->cells.vertex(c, cv)+off_v;
142		✗	tet_mesh->cells.set_vertex(off_c + c, cv, vid);
143			}
144			// merge vertices
145		✗	double eps = (1e-3)*get_cell_average_edge_size(tet_mesh);
146			{
147		✗	vector<index_t> to_kill(tet_mesh->vertices.nb(), 0);
148			vector<index_t> old2new(tet_mesh->vertices.nb());
149		✗	Geom::colocate(tet_mesh->vertices.point_ptr(0), 3, tet_mesh->vertices.nb(), old2new, eps);
150		✗	FOR(c, tet_mesh->cells.nb()) FOR(cv, tet_mesh->cells.nb_vertices(c))
151		✗	tet_mesh->cells.set_vertex(c, cv, old2new[tet_mesh->cells.vertex(c, cv)]);
152		✗	FOR(v, tet_mesh->vertices.nb()) if (old2new[v] != v) to_kill[v] = NOT_AN_ID;
153		✗	tet_mesh->vertices.delete_elements(to_kill);
154			}
155
156
157			}
158
159			/*
160			* _____ _ ____ _ _
161			* / ____| (_) | _ \ | | (_)
162			* | | __ _ _ __ _ __ _ ___ _ __ | |_) | __ _ _ _ __| | ___ _ _ __
163			* | | / _` | '__| '__| |/ _ \ '__| | _ < / _` | | | |/ _` |/ _ \| | '_ \
164			* | |___| (_| | | | | | | __/ | | |_) | (_| | |_| | (_| | (_) | | | | |
165			* \_____\__,_|_| |_| |_|\___|_| |____/ \__,_|\__,_|\__,_|\___/|_|_| |_|
166			*/
167
168
169
170		✗	static bool in_volume(Mesh* surface, vec3 request) {
171			int accum = 0;
172			vec2 R(request[0], request[1]);
173		✗	FOR(f, surface->facets.nb()) {
174		✗	FOR(fan, surface->facets.nb_vertices(f) - 2) {
175			index_t v[3] = {
176			surface->facets.vertex(f,0),
177		✗	surface->facets.vertex(f,1 + fan),
178		✗	surface->facets.vertex(f,2 + fan)
179		✗	};
180		✗	vec2 P[3];
181		✗	FOR(vid, 3) P[vid] = vec2(X(surface)[v[vid]][0], X(surface)[v[vid]][1]) - R;
182			bool in_triangle = true;
183			double tr_orient = det(P[1] - P[0], P[2] - P[0]);
184		✗	if (tr_orient > 0)tr_orient = 1; else tr_orient = -1;
185		✗	FOR(vid, 3) in_triangle = in_triangle && ((det(P[vid], P[(vid + 1) % 3]) > 0) == (tr_orient >0)); ;
186		✗	if (!in_triangle) continue;
187
188		✗	int sign = orient_3d_filter(X(surface)[v[0]].data(), X(surface)[v[1]].data(), X(surface)[v[2]].data(), request.data());
189		✗	if (sign == 0) return false; // I don't want to deal with degenerate case
190		✗	accum += sign;
191			}
192			}
193		✗	return accum != 0;
194
195			}
196
197			struct GrowPt {
198		✗	GrowPt(vec3 p_pos, mat3 p_r) { pos = p_pos; r = p_r; }
199			vec3 pos;
200			mat3 r;
201			};
202
203		✗	void Baudoin_mesher(Mesh* m) {
204			geo_argused(m);
205			//return;
206			// init rot
207		✗	Attribute<mat3> B(m->vertices.attributes(), "B");
208			{
209			vector<vector<vec3> > dir(m->vertices.nb());
210		✗	FOR(f, m->facets.nb()) {
211		✗	if (m->facets.nb_vertices(f) != 4) continue;
212		✗	FOR(e, 4)
213		✗	dir[m->facets.vertex(f, e)].push_back(normalize(X(m)[m->facets.vertex(f, e)] - X(m)[m->facets.vertex(f, (e + 1) % 4)]));
214			}
215		✗	FOR(f, m->facets.nb()) {
216		✗	if (m->facets.nb_vertices(f)!=4) continue;
217		✗	FOR(e, 4) dir[m->facets.vertex(f, e)].clear();
218			}
219		✗	FOR(v, m->vertices.nb())
220		✗	if (dir[v].empty()) B[v].load_identity();
221		✗	else B[v] = Frame::representative_frame(dir[v]);
222			}
223
224			// compute ave QUAD edge length
225			double ave = 0;
226			double nb = 0;
227		✗	FOR(f, m->facets.nb()) {
228		✗	if (m->facets.nb_vertices(f) != 4) continue;
229		✗	FOR(e, 4) {
230		✗	ave += (X(m)[m->facets.vertex(f, e)] - X(m)[m->facets.vertex(f, (e + 1) % 4)]).length();
231		✗	nb += 1.;
232			}
233
234			}
235		✗	ave /= nb;
236
237
238
239			vector<vec3> nvvertices;
240
241			vector<GrowPt> g;
242			index_t cur = 0;
243			// add useless boundary vertices (to prevent intersections)
244		✗	FOR(v, m->vertices.nb()) if (B[v].is_identity()) g.push_back(GrowPt(X(m)[v], B[v]));
245			cur = g.size();
246		✗	FOR(v, m->vertices.nb()) if (!B[v].is_identity()) g.push_back(GrowPt(X(m)[v], B[v]));
247
248
249			// original paper
250			int nb_max_pts = 3000;
251		✗	while (cur < g.size() && nb_max_pts-->0) {
252		✗	FOR(d, 3) FOR(s, 2) {
253			vec3 cand = ave * col(g[cur].r,d);
254		✗	if (s > 0) cand *= -1;
255		✗	cand = cand + g[cur].pos;
256			bool fail = false;
257		✗	FOR(i, g.size())
258		✗	if ((g[i].pos - cand).length2() < pow(.5*ave, 2.))
259			fail = true;
260		✗	if (!fail && in_volume(m, cand)) {
261		✗	g.push_back(GrowPt(cand, g[cur].r));
262			nvvertices.push_back(cand);
263			}
264			}
265		✗	cur++;
266			}
267
268			index_t off_v = m->vertices.create_vertices(nvvertices.size());
269		✗	FOR(i, nvvertices.size()) X(m)[off_v + i] = nvvertices[i];
270
271		✗	}
272
273
274		✗	static void stuff_with_tets_and_pyramids(Mesh* m) {
275		✗	if (m->vertices.nb() == 0 || m->facets.nb()==0) return;
276		✗	m->edges.clear();
277		✗	check_no_intersecting_faces(m);
278		✗	double ave_edge_length = get_facet_average_edge_size(m);
279
280			vector<index_t> pyrindex;
281			vector<vec3> pyr_Z;
282			vector<index_t> pyr_top_index;
283
284			{// split quads into 4 triangles... and remember them to produce pyramids
285			// Remark: spliting quads in 2 may be better on boundary... but it is not clear for constrained boundary...
286		✗	vector<index_t> to_kill(m->facets.nb(), 0);
287
288
289			index_t init_nb_facets = m->facets.nb();
290		✗	FOR(f, init_nb_facets) {
291		✗	if (m->facets.nb_vertices(f) != 4) {
292		✗	geo_assert(m->facets.nb_vertices(f)==3);
293		✗	continue;
294			}
295		✗	index_t nvv = m->vertices.create_vertex();
296		✗	pyr_Z.push_back(facet_normal(m, f));
297
298			//double d = 0;
299			//FOR(e, 4) d += .25*(X(m)[m->facets.vertex(f, e)] - X(m)[m->facets.vertex(f, (e + 1) % 4)]).length();
300		✗	X(m)[nvv] = facet_bary(m, f);// +.2*d*n;
301		✗	to_kill[f] = 1;
302
303		✗	index_t off_f = m->facets.create_triangles(4);
304		✗	FOR(e, 4) {
305		✗	to_kill.push_back(0);
306		✗	m->facets.set_vertex(off_f + e, 0, m->facets.vertex(f, e));
307		✗	m->facets.set_vertex(off_f + e, 1, m->facets.vertex(f, (e + 1) % 4));
308			m->facets.set_vertex(off_f + e, 2, nvv);
309			}
310		✗	FOR(e, 4) pyrindex.push_back(m->facets.vertex(f, e));
311			pyrindex.push_back(nvv);
312			pyr_top_index.push_back(nvv);
313			}
314		✗	m->facets.delete_elements(to_kill, false);
315			}
316
317			{// mode the tip of pyramid in the normal direction to produce better shape (not flat)
318		✗	check_no_intersecting_faces(m);
319		✗	vector<double> max_pyr_top_coeff(pyr_top_index.size(), .5);
320			bool done = false;
321		✗	while (!done){
322		✗	Mesh copy;
323		✗	copy.copy(*m);
324		✗	create_non_manifold_facet_adjacence(&copy);
325
326		✗	FOR(v, pyr_top_index.size())
327		✗	X(&copy)[pyr_top_index[v]] = X(m)[pyr_top_index[v]]
328		✗	- max_pyr_top_coeff[v]*ave_edge_length * pyr_Z[v];
329
330		✗	vector<index_t> intersections = get_intersecting_faces(&copy);
331		✗	done = (intersections.size() == 0);
332		✗	vector<bool> vertexpb(copy.vertices.nb(), false);
333		✗	FOR(i, intersections.size()) FOR(lv, 3) vertexpb[copy.facets.vertex(intersections[i], lv)] = true;
334
335		✗	FOR(v, pyr_top_index.size()) if (vertexpb[pyr_top_index[v]]) {
336		✗	if (max_pyr_top_coeff[v] > .02)max_pyr_top_coeff[v] /= 2.;
337		✗	else max_pyr_top_coeff[v] = 0;
338			}
339		✗	}
340
341		✗	FOR(v, pyr_top_index.size()) X(m)[pyr_top_index[v]]
342		✗	-= std::max(0.0, .7*max_pyr_top_coeff[v])*ave_edge_length * pyr_Z[v];
343
344		✗	check_no_intersecting_faces(m);
345			}
346			// keep vertices geometry (indices are broken by the tetrahedrisation)
347			vector<vec3> pyrpos(pyrindex.size());
348		✗	FOR(nv, pyrindex.size()) pyrpos[nv] = X(m)[pyrindex[nv]];
349
350			// tetrahedrize inside
351			try {
352		✗	FOR(f, m->facets.nb()) geo_assert(m->facets.nb_vertices(f)==3);
353
354		✗	mesh_save(*m, "C:/DATA/debug/pretriangulate.geogram");
355		✗	m->facets.triangulate();
356		✗	mesh_save(*m, "C:/DATA/debug/posttrinagulate.geogram");
357		✗	create_non_manifold_facet_adjacence(m);
358			mesh_tetrahedralize(*m, false, true, 1.);
359			}
360		✗	catch (const Delaunay::InvalidInput& error_report) {
361		✗	FOR(i, error_report.invalid_facets.size())
362		✗	plop(error_report.invalid_facets[i]);
363		✗	Attribute<int> intersection(m->facets.attributes(), "intersection");
364		✗	FOR(f, m->facets.nb()) intersection[f] = 0;
365		✗	FOR(i, error_report.invalid_facets.size())
366		✗	intersection[error_report.invalid_facets[i]] = 1;
367		✗	mesh_save(*m, "C:/DATA/debug/intersectingsurface2.geogram");
368		✗	geo_assert_not_reached;
369		✗	}
370
371			// restore indices from geometry
372		✗	NearestNeighborSearch_var NN = NearestNeighborSearch::create(3);
373		✗	NN->set_points(m->vertices.nb(), m->vertices.point_ptr(0), 3);
374		✗	FOR(nv, pyrindex.size()) pyrindex[nv] = NN->get_nearest_neighbor(pyrpos[nv].data());
375
376
377			// produce pyramids
378		✗	index_t off_c = m->cells.create_pyramids(pyrindex.size() / 5);
379		✗	FOR(p, pyrindex.size() / 5) FOR(lv, 5)
380		✗	m->cells.set_vertex(off_c + p, lv, pyrindex[5 * p + lv]);
381
382			}
383
384
385
386		✗	void hex_dominant(Mesh* cavity, Mesh* hexahedrons, Mesh* result) {
387		✗	plop("HexDominant with vertex_puncher and pyramids");
388		✗	hex_crunch(cavity, hexahedrons);
389
390			//return;
391		✗	Mesh tets;
392		✗	tets.copy(*cavity);
393		✗	stuff_with_tets_and_pyramids(&tets);
394		✗	result->copy(tets);
395		✗	result->facets.clear();
396
397
398
399		✗	index_t off_c = result->cells.create_hexes(hexahedrons->cells.nb());
400			index_t off_v = result->vertices.create_vertices(hexahedrons->vertices.nb());
401		✗	FOR(v, hexahedrons->vertices.nb()) X(result)[off_v + v] = X(hexahedrons)[v];
402		✗	FOR(c, hexahedrons->cells.nb())FOR(cv, 8) {
403		✗	result->cells.set_vertex(off_c + c, cv, off_v + hexahedrons->cells.vertex(c, cv));
404			}
405
406
407			// merge vertices
408		✗	double eps = (1e-3)*get_cell_average_edge_size(result);
409			{
410		✗	vector<index_t> to_kill(result->vertices.nb(), 0);
411			vector<index_t> old2new(result->vertices.nb());
412		✗	Geom::colocate(result->vertices.point_ptr(0), 3, result->vertices.nb(), old2new, eps);
413		✗	FOR(c, result->cells.nb()) FOR(cv, result->cells.nb_vertices(c))
414		✗	result->cells.set_vertex(c, cv, old2new[result->cells.vertex(c, cv)]);
415		✗	FOR(v, result->vertices.nb()) if (old2new[v] != v) to_kill[v] = NOT_AN_ID;
416		✗	result->vertices.delete_elements(to_kill);
417			}
418
419		✗	result->cells.connect();
420		✗	result->facets.clear();
421			// result->cells.compute_borders();
422		✗	return;
423		✗	}
424
425			}
426