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/mesh_inspector.h>
41			#include <exploragram/hexdom/basic.h>
42			#include <exploragram/hexdom/extra_connectivity.h>
43			#include <geogram/mesh/mesh_tetrahedralize.h>
44			#include <geogram/delaunay/delaunay.h>
45
46
47			namespace GEO {
48
49		✗	bool volume_boundary_is_manifold(Mesh* m, std::string& msg) {
50		✗	m->cells.compute_borders();
51		✗	if (!surface_is_manifold(m, msg)) {
52			return false;
53			}
54		✗	m->facets.clear();
55		✗	return true;
56			}
57
58
59		✗	bool have_negative_tet_volume(Mesh*m) {
60		✗	FOR(c, m->cells.nb()) {
61		✗	vec3 A = X(m)[m->cells.vertex(c, 0)];
62		✗	vec3 B = X(m)[m->cells.vertex(c, 1)];
63		✗	vec3 C = X(m)[m->cells.vertex(c, 2)];
64		✗	vec3 D = X(m)[m->cells.vertex(c, 3)];
65			double vol = dot(cross(B - A, C - A), D - A);
66		✗	if (vol < 0) {
67		✗	Attribute<double> signed_volume(m->cells.attributes(), "signed_volume");
68		✗	signed_volume[c] = vol;
69			return true;
70			}
71			}
72			return false;
73			}
74
75		✗	bool surface_is_tetgenifiable(Mesh* m) {
76		✗	Mesh copy;
77		✗	copy.copy(*m);
78		✗	create_non_manifold_facet_adjacence(&copy);
79		✗	copy.facets.triangulate();
80			try {
81			mesh_tetrahedralize(copy, false, false, 1.);
82			}
83		✗	catch (const GEO::Delaunay::InvalidInput& error_report) {
84		✗	FOR(i, error_report.invalid_facets.size()) {
85		✗	plop(error_report.invalid_facets[i]);
86			}
87			return false;
88		✗	}
89		✗	return true;
90		✗	}
91
92		✗	bool volume_is_tetgenifiable(Mesh* m) {
93		✗	Mesh copy;
94		✗	copy.copy(*m);
95		✗	copy.edges.clear();
96		✗	copy.cells.compute_borders();
97		✗	copy.cells.clear();
98		✗	return surface_is_tetgenifiable(&copy);
99		✗	}
100		✗	bool surface_is_manifold(Mesh* m, std::string& msg) {
101		✗	if (m->facets.nb() == 0) return true;
102			{
103			// check for duplicated corners around a face
104		✗	FOR(f, m->facets.nb()) FOR(fc, m->facets.nb_corners(f))
105		✗	if (m->facets.vertex(f, fc) == m->facets.vertex(f, next_mod(fc, m->facets.nb_corners(f)))) {
106		✗	msg = "duplicated corner detected on (face = " + String::to_string(f) + " , local corner = " + String::to_string(fc) + " , vertex = " +
107		✗	String::to_string(m->facets.vertex(f, fc));
108		✗	return false;
109			}
110			// output the type of surface
111			index_t nb_edges_par_facets = m->facets.nb_corners(0);
112		✗	FOR(f, m->facets.nb()) if (m->facets.nb_corners(f) != nb_edges_par_facets) nb_edges_par_facets = index_t(-1);
113		✗	if (nb_edges_par_facets != index_t(-1)) plop(nb_edges_par_facets);
114
115			// check if the mesh is manifold
116
117		✗	Attribute<int> nb_opp(m->facet_corners.attributes(), "nb_opp");
118		✗	Attribute<int> nb_occ(m->facet_corners.attributes(), "nb_occ");
119		✗	FOR(h, m->facet_corners.nb()) { nb_opp[h] = 0; nb_occ[h] = 0; }
120
121			// edge connectivity
122		✗	FacetsExtraConnectivity fec(m);
123		✗	int nb_0_opp = 0;
124			//int nb_1_opp = 0;
125		✗	int nb_multiple_opp = 0;
126		✗	int nb_duplicated_edge = 0;
127		✗	FOR(h, m->facet_corners.nb()) {
128			index_t cir = h;
129			index_t result = NOT_AN_ID; // not found
130			do {
131		✗	index_t candidate = fec.prev(cir);
132		✗	if ((fec.org(candidate) == fec.dest(h)) && (fec.dest(candidate) == fec.org(h))) {
133		✗	nb_opp[h]++;
134		✗	if (result == NOT_AN_ID) result = candidate;
135		✗	else nb_multiple_opp++;
136			}
137		✗	if (cir != h && fec.dest(h) == fec.dest(cir)) {
138		✗	nb_duplicated_edge++;
139		✗	nb_occ[h]++;
140			}
141		✗	cir = fec.c2c[cir];
142		✗	} while (cir != h);
143		✗	if (result == NOT_AN_ID)nb_0_opp++;
144			//else nb_1_opp++;
145			}
146
147
148		✗	if (nb_0_opp > 0) {
149		✗	msg = "surface have halfedges without opposite, nb= " + String::to_string(nb_0_opp);
150		✗	return false;
151			}
152		✗	if (nb_multiple_opp > 0) {
153		✗	msg = "surface have halfedge with more than 2 opposites, nb= " + String::to_string(nb_multiple_opp);
154		✗	return false;
155			}
156		✗	if (nb_duplicated_edge > 0) {
157		✗	msg = "halfedge appears in more than one facet, nb= " + String::to_string(nb_duplicated_edge);
158		✗	return false;
159			}
160
161			// check for non manifold vertices
162		✗	Attribute<bool> nonmanifold(m->vertices.attributes(), "nonmanifold");
163		✗	FOR(v, m->vertices.nb()) nonmanifold[v] = false;
164
165		✗	FOR(h, m->facet_corners.nb()) {
166		✗	if (nb_opp[h] != 1 || nb_occ[h] != 0)
167		✗	nonmanifold[fec.org(h)] = true;
168			}
169		✗	vector<int> val(m->vertices.nb(), 0);
170		✗	FOR(f, m->facets.nb()) FOR(lc, m->facets.nb_vertices(f)) val[m->facets.vertex(f, lc)]++;
171		✗	FOR(h, m->facet_corners.nb()) {
172			int nb = 0;
173			index_t cir = h;
174			do {
175		✗	nb++;
176		✗	cir = fec.next_around_vertex(cir);// fec.c2c[cir];
177		✗	} while (cir != h);
178		✗	if (nb != val[fec.org(h)]) {
179		✗	msg = "Vertex " + String::to_string(fec.org(h)) + " is non manifold ";
180			return false;
181			}
182			}
183		✗	}
184		✗	m->vertices.attributes().delete_attribute_store("nonmanifold");
185		✗	m->facet_corners.attributes().delete_attribute_store("nb_opp");
186		✗	m->facet_corners.attributes().delete_attribute_store("nb_occ");
187		✗	return true;
188			}
189		✗	void get_facet_stats(Mesh* m, const char * msg, bool export_attribs) {
190			geo_argused(export_attribs);
191		✗	GEO::Logger::out("HexDom") << "-----------------------------------------" << std::endl;
192		✗	GEO::Logger::out("HexDom") << "get_facet_stats " << msg << std::endl;
193		✗	GEO::Logger::out("HexDom") << "-----------------------------------------" << std::endl;
194
195			{
196		✗	Attribute<int> nb_opp(m->facet_corners.attributes(), "nb_opp");
197		✗	Attribute<int> nb_occ(m->facet_corners.attributes(), "nb_occ");
198		✗	FOR(h, m->facet_corners.nb()) nb_opp[h] = 0;
199
200			// edge connectivity
201		✗	FacetsExtraConnectivity fec(m);
202			int nb_0_opp = 0;
203			int nb_1_opp = 0;
204			int nb_multiple_opp = 0;
205			int nb_duplicated_edge = 0;
206		✗	FOR(h, m->facet_corners.nb()) {
207			index_t cir = h;
208			index_t result = NOT_AN_ID; // not found
209			do {
210		✗	index_t candidate = fec.prev(cir);
211		✗	if ((fec.org(candidate) == fec.dest(h)) && (fec.dest(candidate) == fec.org(h))) {
212		✗	nb_opp[h]++;
213		✗	if (result == NOT_AN_ID) result = candidate;
214		✗	else nb_multiple_opp++;
215			}
216		✗	if (cir != h && fec.dest(h) == fec.dest(cir)) {
217		✗	nb_duplicated_edge++;
218		✗	nb_occ[h]++;
219			}
220		✗	cir = fec.c2c[cir];
221		✗	} while (cir != h);
222		✗	if (result == NOT_AN_ID)nb_1_opp++;
223		✗	else nb_0_opp++;
224			}
225
226
227
228		✗	FOR(f, m->facets.nb()) FOR(fc, m->facets.nb_corners(f))
229		✗	if (m->facets.vertex(f, fc) == m->facets.vertex(f, next_mod(fc, m->facets.nb_corners(f))))
230		✗	GEO::Logger::out("HexDom") << "Duplicated vertex found at facet #" << f << ", local corner= " << fc << " and vertex is " << m->facets.vertex(f, fc) << std::endl;
231
232		✗	plop(nb_0_opp);
233		✗	plop(nb_1_opp);
234		✗	plop(nb_multiple_opp);
235		✗	plop(nb_duplicated_edge);
236
237			// check for non manifold vertices
238		✗	Attribute<bool> nonmanifold(m->vertices.attributes(), "nonmanifold");
239		✗	FOR(v, m->vertices.nb()) nonmanifold[v] = false;
240
241		✗	FOR(h, m->facet_corners.nb()) {
242		✗	if (nb_opp[h] != 1 || nb_occ[h] != 0)
243		✗	nonmanifold[fec.org(h)] = true;
244
245			}
246
247		✗	vector<int> val(m->vertices.nb(), 0);
248		✗	FOR(f, m->facets.nb()) FOR(lc, m->facets.nb_vertices(f)) val[m->facets.vertex(f, lc)]++;
249		✗	FOR(h, m->facet_corners.nb()) {
250			int nb = 0;
251			index_t cir = h;
252			do {
253		✗	nb++;
254		✗	cir = fec.c2c[cir];
255		✗	} while (cir != h);
256		✗	if (nb != val[fec.org(h)]) {
257		✗	GEO::Logger::out("HexDom") << "Vertex " << fec.org(h) << " is non-manifold !!" << std::endl;
258		✗	nonmanifold[fec.org(h)] = true;
259			}
260			}
261		✗	}
262		✗	m->vertices.attributes().delete_attribute_store("nonmanifold");
263		✗	m->facet_corners.attributes().delete_attribute_store("nb_opp");
264		✗	m->facet_corners.attributes().delete_attribute_store("nb_occ");
265		✗	}
266		✗	double tet_vol(vec3 A, vec3 B, vec3 C, vec3 D) {
267			B = B - A;
268			C = C - A;
269			D = D - A;
270		✗	return (1. / 6.)*dot(D, cross(B, C));
271			}
272		✗	void get_hex_proportion(Mesh*m, double &nb_hex_prop, double &vol_hex_prop) {
273			int nb_tets = 0;
274			int nb_hexs = 0;
275			double vol_tets = 0;
276			double vol_hexs = 0;
277		✗	FOR(c, m->cells.nb()) if (m->cells.nb_facets(c) == 4) {
278		✗	vol_tets += tet_vol(X(m)[m->cells.vertex(c, 0)], X(m)[m->cells.vertex(c, 1)], X(m)[m->cells.vertex(c, 2)], X(m)[m->cells.vertex(c, 3)]);
279		✗	nb_tets++;
280			}
281		✗	FOR(c, m->cells.nb()) if (m->cells.nb_facets(c) == 6) {
282			vector<vec3> P(8);
283		✗	FOR(cv, 8) P[cv] = X(m)[m->cells.vertex(c, cv)];
284
285		✗	vol_hexs += tet_vol(P[0], P[3], P[2], P[6]);
286		✗	vol_hexs += tet_vol(P[0], P[7], P[3], P[6]);
287		✗	vol_hexs += tet_vol(P[0], P[7], P[6], P[4]);
288
289		✗	vol_hexs += tet_vol(P[0], P[1], P[3], P[7]);
290		✗	vol_hexs += tet_vol(P[0], P[1], P[7], P[5]);
291		✗	vol_hexs += tet_vol(P[0], P[5], P[7], P[4]);
292		✗	nb_hexs++;
293			}
294		✗	if (nb_hexs + nb_tets>0) nb_hex_prop = double(nb_hexs) / double(nb_hexs + nb_tets);
295		✗	if (nb_hexs + nb_tets>0) vol_hex_prop = double(vol_hexs) / double(vol_hexs + vol_tets);
296
297		✗	}
298			}
299