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_candidates.h>
41			#include <exploragram/hexdom/mesh_utils.h>
42			#include <geogram/mesh/mesh.h>
43			#include <geogram/points/colocate.h>
44			#include <exploragram/hexdom/geometry.h>
45			namespace GEO {
46
47
48		✗	static vec3 change_tet_basis(vec3 in, vec3 P[4], vec3 P_img[4]) { // TODO sortir ceci + trglgrad dans un fichier apart
49			Matrix<12, double> M;
50			M.load_zero();
51			double RHS[12];
52			double X[12];
53
54			index_t cur_row = 0; // in fact cur_row can be replaced by v*3+dimx
55		✗	FOR(v, 4) FOR(dimx, 3) {
56		✗	FOR(dimu, 3)
57		✗	M(4 * dimx + dimu, cur_row) = P[v][dimu];
58		✗	M(4 * dimx + 3, cur_row) = 1.0;
59		✗	RHS[cur_row] = P_img[v][dimx];
60		✗	cur_row++;
61
62			}
63		✗	geo_assert(cur_row == 12);
64
65		✗	Matrix<12, double> inv = M.inverse();
66
67		✗	FOR(i, 12) X[i] = 0;
68		✗	FOR(i, 12) FOR(j, 12) X[i] += inv(j, i)*RHS[j];
69
70			vec3 res(0, 0, 0);
71		✗	FOR(dimx, 3) FOR(dimu, 3)
72		✗	res[dimx] += in[dimu] * X[4 * dimx + dimu];
73
74		✗	FOR(dimx, 3) res[dimx] += X[4 * dimx + 3];
75		✗	return res;
76			}
77
78		✗	static bool in_tet(vec3 test, vec3 P[4], double eps) {
79		✗	int find[4][3] = { { 1, 3, 2 },{ 3, 0, 2 },{ 0, 3, 1 },{ 0, 1, 2 } };
80		✗	if (std::abs(dot(P[3] - P[0], cross(P[2] - P[0], P[1] - P[0]))) < 1e-15) return false;// check for flat tet
81		✗	FOR(f, 4) {
82		✗	vec3 vA = P[find[f][0]];
83		✗	vec3 vB = P[find[f][1]];
84		✗	vec3 vC = P[find[f][2]];
85			vec3 n = cross(vB - vA, vC - vA);
86		✗	n = normalize(n);
87		✗	if (dot(test - vA, n) > eps) return false;
88			}
89			return true;
90			}
91
92		✗	static void get_grid_vertices(Mesh* m, Attribute<vec3>& UC, index_t t, std::vector<vec3>& psetX, std::vector<vec3>& psetU, bool dual) {
93			// compute a bbox in the parametric domain
94		✗	int bbox[2][3] = { { 1000000, 1000000, 1000000 },{ -1000000, -1000000, -1000000 } };
95		✗	FOR(i, 4) FOR(dim, 3) {
96		✗	bbox[0][dim] = std::min(bbox[0][dim], int(std::floor(UC[m->cells.corner(t, i)][dim])) - 1);
97		✗	bbox[1][dim] = std::max(bbox[1][dim], int(std::ceil(UC[m->cells.corner(t, i)][dim])) + 1);
98			}
99
100		✗	vec3 lX[4], lU[4];
101		✗	FOR(i, 4) {
102		✗	lX[i] = m->vertices.point(m->cells.vertex(t, i));
103		✗	lU[i] = UC[m->cells.corner(t, i)];
104			}
105
106			// raster the bbox and project to geometric space points included in the parametric space tet
107		✗	for (int i0 = bbox[0][0]; i0 <= bbox[1][0]; i0++) {
108		✗	for (int i1 = bbox[0][1]; i1 <= bbox[1][1]; i1++) {
109		✗	for (int i2 = bbox[0][2]; i2 <= bbox[1][2]; i2++) {
110		✗	vec3 testpt(i0, i1, i2);
111		✗	if (dual) testpt += vec3(.5, .5, .5);
112		✗	if (in_tet(testpt, lU, .001)) {
113		✗	psetX.push_back(change_tet_basis(testpt, lU, lX));
114			psetU.push_back(testpt);
115			}
116			}
117			}
118			}
119		✗	}
120
121
122
123
124		✗	void export_points(Mesh* m, Mesh* pts) {
125		✗	Attribute<bool> has_param(m->cell_facets.attributes(), "has_param");
126		✗	Attribute<vec3> UC(m->cell_corners.attributes(), "U");
127
128		✗	FOR(c, m->cells.nb()) {
129		✗	if (c % 100 == 0) GEO::Logger::out("HexDom") << " EXPORT point set: tet" << c << " / " << m->cells.nb() << std::endl;
130		✗	if (!has_param[m->cells.facet(c, 0)] || !has_param[m->cells.facet(c, 1)] || !has_param[m->cells.facet(c, 2)] || !has_param[m->cells.facet(c, 3)]) continue;
131
132			std::vector<vec3> seedU;
133			std::vector<vec3> seedX;
134		✗	get_grid_vertices(m, UC, c, seedX, seedU, false);
135		✗	index_t off = pts->vertices.create_vertices(index_t(seedX.size()));
136		✗	FOR(lv, seedX.size()) X(pts)[off + lv] = seedX[lv];
137			}
138
139			// remove colocated vertices
140		✗	double eps = (1e-2)*get_cell_average_edge_size(m);
141		✗	vector<index_t> to_kill(pts->vertices.nb(), 0);
142			vector<index_t> old2new(pts->vertices.nb());
143		✗	Geom::colocate(pts->vertices.point_ptr(0), 3, pts->vertices.nb(), old2new, eps);
144		✗	FOR(v, pts->vertices.nb()) if (old2new[v] != v) to_kill[v] = NOT_AN_ID;
145		✗	plop(pts->vertices.nb());
146		✗	pts->vertices.delete_elements(to_kill);
147		✗	plop(pts->vertices.nb());
148
149		✗	}
150		✗	inline bool intersect_unit_box(vec3 p_boxcenter, vec3 tri[3]) {
151		✗	float boxcenter[3] = { float(p_boxcenter[0]), float(p_boxcenter[1]), float(p_boxcenter[2]) };
152		✗	float boxhalfsize[3] = { .5, .5, .5 };
153			float triverts[3][3] = {
154		✗	{ float(tri[0][0]), float(tri[0][1]), float(tri[0][2]) },
155		✗	{ float(tri[1][0]), float(tri[1][1]), float(tri[1][2]) },
156		✗	{ float(tri[2][0]), float(tri[2][1]), float(tri[2][2]) }
157		✗	};
158		✗	return (triBoxOverlap(boxcenter, boxhalfsize, triverts) != 0);
159			}
160
161		✗	inline void init_centered_unit_cube_face_bary(vec3 *cube_face_bary) {
162		✗	vec3 cubeU[8];
163		✗	FOR(k, 2)FOR(j, 2)FOR(i, 2)
164		✗	cubeU[4 * i + 2 * j + k] = vec3(i, j, k) - vec3(.5, .5, .5);
165		✗	CellDescriptor hexdescr = MeshCellsStore::cell_type_to_cell_descriptor(MESH_HEX);
166		✗	FOR(lf, 6) {
167		✗	cube_face_bary[lf] = vec3(0, 0, 0);
168		✗	for (int lv = 0; lv < 4; lv++) cube_face_bary[lf] += .5 * cubeU[hexdescr.facet_vertex[lf][lv]];
169			}
170		✗	}
171
172			static bool cell_has_param(Mesh* m, Attribute<bool>& has_param, index_t c) {
173		✗	FOR(f, 4) if (!has_param[m->cells.facet(c, f)]) return false;
174			return true;
175			}
176
177		✗	static bool param_is_degenerated(Mesh* m, Attribute<vec3>& UC, index_t tet) {
178		✗	vec3 P[4];
179		✗	FOR(i, 4) P[i] = UC[m->cells.corner(tet, i)];
180		✗	return (dot(P[3] - P[0], cross(P[2] - P[0], P[1] - P[0])) == 0);
181			}
182
183
184		✗	static void make_neig_tet_compatible(Mesh* m, Attribute<vec3>& UC, index_t ref_tet, index_t cf) {
185			{// preconditions
186		✗	Attribute<bool> has_param(m->cell_facets.attributes(), "has_param");
187		✗	geo_assert(cell_has_param(m, has_param, ref_tet));
188		✗	geo_assert(m->cells.adjacent(ref_tet, cf) != NOT_AN_ID);
189		✗	geo_assert(cell_has_param(m, has_param, m->cells.adjacent(ref_tet, cf)));
190			}
191
192			index_t opp_tet = m->cells.adjacent(ref_tet, cf);
193			index_t ABC[3];
194		✗	FOR(cfv, 3) ABC[cfv] = m->cells.facet_vertex(ref_tet, cf, cfv);
195
196			index_t ref_c[3];
197			index_t opp_c[3];
198
199			index_t Dc = NOT_AN_ID;
200
201
202		✗	FOR(opp_corner, 4) {
203			index_t v = m->cells.vertex(opp_tet, opp_corner);
204			bool match_found = false;
205		✗	FOR(i, 3) if (ABC[i] == v) {
206		✗	opp_c[i] = m->cells.corner(opp_tet, opp_corner);
207			match_found = true;
208			}
209		✗	if (!match_found) {
210			Dc = m->cells.corner(opp_tet, opp_corner);
211			}
212			}
213		✗	geo_assert(Dc != NOT_AN_ID);
214		✗	FOR(i, 3) ref_c[i] = cell_facet_corner_id(m, ref_tet, cf, i);
215
216		✗	vec3 ref_xyz[3];
217		✗	ref_xyz[0] = normalize(UC[ref_c[1]] - UC[ref_c[0]]);
218		✗	ref_xyz[1] = normalize(UC[ref_c[2]] - UC[ref_c[0]]);
219		✗	ref_xyz[2] = normalize(cross(ref_xyz[0], ref_xyz[1]));
220		✗	ref_xyz[1] = normalize(cross(ref_xyz[2], ref_xyz[0]));
221
222		✗	vec3 opp_xyz[3];
223		✗	opp_xyz[0] = normalize(UC[opp_c[1]] - UC[opp_c[0]]);
224		✗	opp_xyz[1] = normalize(UC[opp_c[2]] - UC[opp_c[0]]);
225		✗	opp_xyz[2] = normalize(cross(opp_xyz[0], opp_xyz[1]));
226		✗	opp_xyz[1] = normalize(cross(opp_xyz[2], opp_xyz[0]));
227
228			vec3 local;
229		✗	FOR(d, 3)local[d] = dot(opp_xyz[d], UC[Dc] - UC[opp_c[0]]);
230
231		✗	vec3 save = UC[Dc];
232		✗	UC[Dc] = UC[ref_c[0]];
233		✗	FOR(d, 3) UC[Dc] = UC[Dc] + local[d] * ref_xyz[d];
234
235		✗	FOR(v, 3) UC[opp_c[v]] = UC[ref_c[v]];
236
237
238			//assumes that two coordinates are never (non integer and closer than 1e-8)
239		✗	FOR(d, 3) FOR(ds, 3) {
240		✗	if (std::abs((UC[Dc][d] - round(UC[Dc][d])) - (save[ds] - round(save[ds]))) < 1e-8) {
241		✗	UC[Dc][d] = round(UC[Dc][d]) + save[ds] - round(save[ds]);
242			}
243		✗	if (std::abs((UC[Dc][d] - round(UC[Dc][d])) + (save[ds] - round(save[ds]))) < 1e-8) {
244		✗	UC[Dc][d] = round(UC[Dc][d]) - save[ds] + round(save[ds]);
245			}
246			}
247
248			//FOR(d, 3) if (std::abs(UC[Dc][d] - round(UC[Dc][d])) <.05) UC[Dc][d] = round(UC[Dc][d]);
249
250		✗	}
251
252		✗	void export_hexes(Mesh* m, Mesh* hex) {
253		✗	Attribute<bool> has_param(m->cell_facets.attributes(), "has_param");
254		✗	Attribute<vec3> UC(m->cell_corners.attributes(), "U");
255
256		✗	vec3 cube_face_bary[6];
257		✗	init_centered_unit_cube_face_bary(cube_face_bary);
258
259
260
261		✗	FOR(c, m->cells.nb()) {
262
263		✗	if (c % 1000 == 0) GEO::Logger::out("HexDom") << " EXPORT HEXES tet = " << c << " / " << m->cells.nb() << std::endl;
264		✗	if (!cell_has_param(m, has_param, c)) continue;
265		✗	if (param_is_degenerated(m, UC, c)) continue;
266			std::vector<vec3> seedU; // contains centers of all cubes inside tet c
267			{// init seedU
268			std::vector<vec3> trash;
269		✗	get_grid_vertices(m, UC, c, trash, seedU, true);
270			}
271		✗	FOR(sid, seedU.size()) {
272		✗	vec3 ptsU[8]; // vertices of the cube centered in seedU[sid]
273		✗	vec3 ptsX[8];
274		✗	bool ptsdone[8] = {};
275			bool have_boundary_face[6] = {};
276		✗	FOR(k, 2) FOR(j, 2) FOR(i, 2)
277		✗	ptsU[4 * i + 2 * j + k] = seedU[sid] - vec3(.5, .5, .5) + vec3(i, j, k);
278
279
280			std::vector<index_t> tet_stack;
281			tet_stack.push_back(c);
282			std::vector<index_t> tet_done;
283			bool has_sing_tet = false;
284
285		✗	while (!tet_stack.empty()) {
286		✗	index_t curt = tet_stack.back(); tet_stack.pop_back();
287
288			// check if tet is singular
289		✗	if (!cell_has_param(m, has_param, curt)) { has_sing_tet = true; break; }
290
291
292			// create local geometry
293		✗	vec3 lX[4], lU[4];
294		✗	FOR(i, 4) {
295		✗	lX[i] = m->vertices.point(m->cells.vertex(curt, i));
296		✗	lU[i] = UC[m->cells.corner(curt, i)];
297			}
298
299			// find cubes corners located inside the current tet
300		✗	FOR(i, 8) if (in_tet(ptsU[i], lU, 1e-5)) {
301		✗	ptsX[i] = change_tet_basis(ptsU[i], lU, lX);
302		✗	ptsdone[i] = true;
303			}
304
305		✗	FOR(lf, 4) {
306			index_t corners[3] = {
307		✗	cell_facet_corner_id(m, curt, lf, 0),
308		✗	cell_facet_corner_id(m, curt, lf, 1),
309		✗	cell_facet_corner_id(m, curt, lf, 2)
310			};
311			//index_t corners[3];
312			//FOR(v, 3) FOR(corn, 4) if (verts[v] == m->cell_corners.vertex(m->cells.corner(curt, corn)))
313			// corners[v] = m->cells.corner(curt, corn);
314
315
316		✗	vec3 tri[3] = { UC[corners[0]], UC[corners[1]], UC[corners[2]] };
317
318		✗	index_t oppt = m->cells.adjacent(curt, lf);
319		✗	if (oppt == NO_CELL) {
320			// if the facet matches a unit cube facet in U coordinates, mark this unit cube facet
321		✗	vec3 face_normal = normalize(cross(tri[1] - tri[0], tri[2] - tri[0]));
322		✗	FOR(cubef, 6)
323			if ((face_normal - cube_face_bary[cubef]).length2() < 1e-10
324		✗	&& round(.5 + dot(cube_face_bary[cubef], tri[0] - seedU[sid])) == 1) {
325			have_boundary_face[cubef] = true;
326			// need to check that it is not a concave hardedge
327		✗	FOR(cubef_in, 6) {
328		✗	if (cubef == cubef_in) continue;
329
330			bool all_are_outside = true;
331		✗	FOR(d, 3) {
332			all_are_outside = all_are_outside &&
333			(.5 + dot(cube_face_bary[cubef_in], tri[d] - seedU[sid]) > 1 - 1e-5);
334			}
335			have_boundary_face[cubef] = have_boundary_face[cubef] && !all_are_outside;
336			}
337			}
338			continue;
339		✗	}
340
341			bool alreadydone = false;
342		✗	FOR(p, tet_done.size()) alreadydone = alreadydone || (tet_done[p] == oppt);
343		✗	if (alreadydone) continue;
344
345		✗	if (intersect_unit_box(seedU[sid], tri)) {
346		✗	if (cell_has_param(m, has_param, oppt) && !param_is_degenerated(m, UC, oppt)) {
347		✗	make_neig_tet_compatible(m, UC, curt, lf);
348			tet_stack.push_back(oppt);
349			}
350			else {
351			has_sing_tet = true;
352		✗	break;
353			}
354			}
355			}
356			tet_done.push_back(curt);
357			}
358			bool all_vertices_found = true;
359		✗	FOR(i, 8) all_vertices_found = all_vertices_found && ptsdone[i];
360		✗	if (all_vertices_found && !has_sing_tet) {
361
362			index_t off = hex->vertices.create_vertices(8);
363		✗	hex->cells.create_hex(off, off + 1, off + 2, off + 3, off + 4, off + 5, off + 6, off + 7);
364		✗	FOR(i, 8) hex->vertices.point(off + i) = ptsX[i];
365			}
366			}
367			}
368		✗	}
369
370			}
371