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_utils.h>
41			#include <exploragram/hexdom/geometry.h>
42			#include <exploragram/hexdom/mesh_inspector.h>
43			#include <exploragram/hexdom/extra_connectivity.h>
44			#include <geogram/points/colocate.h>
45			#include <geogram/mesh/mesh_tetrahedralize.h>
46			#include <geogram/delaunay/delaunay.h>
47
48			namespace GEO {
49
50		✗	void compute_3D_edge_cot_w(Mesh* m, Attribute<index_t>& v2e, double anisoZ_cotW){
51
52		✗	Attribute<double> cot_w(m->edges.attributes(), "cot_w");
53		✗	FOR(e, m->edges.nb()) cot_w[e] = 0;
54		✗	FOR(c, m->cells.nb()){
55		✗	vec3 pts[4];
56		✗	FOR(cv, 4)pts[cv] = X(m)[m->cells.vertex(c, cv)];
57		✗	double aniso[6] = { 1, 1, anisoZ_cotW, 0, 0, 0 };
58		✗	CoTan3D cot(pts, aniso);
59		✗	FOR(cv, 4){
60		✗	index_t v = m->cells.vertex(c, cv);
61		✗	index_t start = v2e[v];
62		✗	index_t end = m->edges.nb();
63		✗	if (v + 1 < m->vertices.nb()) end = v2e[v + 1];
64		✗	for (index_t e = start; e < end; e++){
65		✗	geo_assert(v == m->edges.vertex(e, 0));
66		✗	FOR(cv2, 4){
67		✗	if (cv == cv2) continue;
68		✗	index_t v2 = m->cells.vertex(c, cv2);
69		✗	if (v2 == m->edges.vertex(e, 1)){
70		✗	FOR(cot_e, 6){
71		✗	if ((cot.org(cot_e) != cv || cot.dest(cot_e) != cv2)
72		✗	&& (cot.dest(cot_e) != cv || cot.org(cot_e) != cv2))
73		✗	continue;
74		✗	cot_w[e] += cot.w[cot_e];
75			}
76			}
77			}
78			}
79			}
80			}
81			// normalize a bit
82		✗	double sum = 0;
83		✗	FOR(e, m->edges.nb()) sum += cot_w[e];
84		✗	FOR(e, m->edges.nb()) cot_w[e] *= double(m->edges.nb()) / sum;
85		✗	}
86
87		✗	void kill_isolated_vertices(Mesh* m){
88		✗	vector<index_t> to_kill(m->vertices.nb(), NOT_AN_ID);
89		✗	FOR(e, m->edges.nb()) FOR(ev, 2) to_kill[m->edges.vertex(e, ev)] = 0;
90		✗	FOR(f, m->facets.nb()) FOR(fv, m->facets.nb_vertices(f)) to_kill[m->facets.vertex(f, fv)] = 0;
91		✗	FOR(c, m->cells.nb()) FOR(cv, m->cells.nb_vertices(c)) to_kill[m->cells.vertex(c, cv)] = 0;
92		✗	m->vertices.delete_elements(to_kill);
93		✗	}
94
95
96		✗	void merge_vertices(Mesh* m, double eps){
97		✗	vector<index_t> to_kill(m->vertices.nb(), 0);
98		✗	vector<index_t> old2new(m->vertices.nb());
99		✗	Geom::colocate(m->vertices.point_ptr(0), 3, m->vertices.nb(), old2new, eps);
100		✗	FOR(e, m->edges.nb()) FOR(ev, 2) m->edges.set_vertex(e, ev, old2new[m->edges.vertex(e, ev)]);
101		✗	FOR(f, m->facets.nb()) FOR(fv, m->facets.nb_vertices(f)) m->facets.set_vertex(f, fv, old2new[m->facets.vertex(f, fv)]);
102		✗	FOR(c, m->cells.nb()) FOR(cv, m->cells.nb_vertices(c)) m->cells.set_vertex(c, cv, old2new[m->cells.vertex(c, cv)]);
103		✗	FOR(v, m->vertices.nb()) if (old2new[v] != v) to_kill[v] = NOT_AN_ID;
104		✗	m->vertices.delete_elements(to_kill);
105		✗	}
106
107		✗	void facets_smooth_geom(Mesh* m, std::vector<bool>& lock_v, double fit_coeff) {
108		✗	vector<vec3> P(m->vertices.nb(), vec3(0, 0, 0));
109		✗	vector<index_t> val(m->vertices.nb(), 0);
110		✗	FOR(f, m->facets.nb()) {
111		✗	vec3 bary = facet_bary(m, f);
112		✗	FOR(fv, m->facets.nb_vertices(f)) {
113		✗	index_t v = m->facets.vertex(f, fv);
114		✗	P[v] = P[v] + bary;
115		✗	val[v]++;
116			}
117			}
118		✗	FOR(v, m->vertices.nb()) if (!lock_v[v]) X(m)[v] = fit_coeff*X(m)[v] + (1. - fit_coeff) * (1. / double(val[v]))*P[v];
119		✗	}
120
121		✗	void cells_smooth_geom(Mesh* m, std::vector<bool>& lock_v, double fit_coeff) {
122		✗	vector<vec3> P(m->vertices.nb(), vec3(0, 0, 0));
123		✗	vector<index_t> val(m->vertices.nb(), 0);
124		✗	FOR(c, m->cells.nb()) {
125		✗	vec3 bary = cell_bary(m, c);
126		✗	FOR(cv, m->cells.nb_vertices(c)) {
127		✗	index_t v = m->cells.vertex(c, cv);
128		✗	P[v] = P[v] + bary;
129		✗	val[v]++;
130			}
131			}
132		✗	FOR(v, m->vertices.nb()) if (!lock_v[v]) X(m)[v] = fit_coeff*X(m)[v] + (1. - fit_coeff) * (1. / double(val[v]))*P[v];
133		✗	}
134
135		✗	void facets_smooth_geom(Mesh* m, double fit_coeff) {
136		✗	std::vector<bool> lock_v(m->vertices.nb(), false);
137		✗	facets_smooth_geom(m, lock_v, fit_coeff);
138		✗	}
139
140		✗	void cells_smooth_geom(Mesh* m, double fit_coeff) {
141		✗	std::vector<bool> lock_v(m->vertices.nb(), false);
142		✗	cells_smooth_geom(m, lock_v, fit_coeff);
143		✗	}
144
145
146
147		✗	void create_non_manifold_facet_adjacence(Mesh* m) {
148		✗	FacetsExtraConnectivity fec(m);
149		✗	FOR(f, m->facets.nb())FOR(lh, m->facets.nb_vertices(f)) m->facets.set_adjacent(f, lh, NOT_AN_ID);
150		✗	FOR(h, m->facet_corners.nb()) {
151		✗	if (m->facets.adjacent(fec.facet(h), fec.local_id(h)) != NOT_AN_ID) continue;
152		✗	index_t cir = h;
153		✗	index_t best_candidate = NOT_AN_ID;
154		✗	double bestdot = 2;
155			do {
156		✗	index_t candidate = fec.prev(cir);
157		✗	if (candidate != NOT_AN_ID
158		✗	&& m->facets.adjacent(fec.facet(candidate), fec.local_id(candidate)) == NOT_AN_ID)
159		✗	if ((fec.org(candidate) == fec.dest(h)) && (fec.dest(candidate) == fec.org(h))) {
160		✗	double curdot = dot(facet_normal(m, fec.facet(h)), -facet_normal(m, fec.facet(candidate)));
161		✗	if (curdot < bestdot) {
162		✗	best_candidate = candidate;
163		✗	bestdot = curdot;
164			}
165			}
166		✗	cir = fec.c2c[cir];
167		✗	} while (cir != h);
168		✗	if (best_candidate != NOT_AN_ID) {
169		✗	m->facets.set_adjacent(fec.facet(h), fec.local_id(h), fec.facet(best_candidate));
170		✗	m->facets.set_adjacent(fec.facet(best_candidate), fec.local_id(best_candidate), fec.facet(h));
171
172		✗	geo_assert(fec.org(h) == fec.dest(best_candidate));
173		✗	geo_assert(fec.dest(h) == fec.org(best_candidate));
174			}
175			}
176		✗	}
177
178		✗	double get_cell_average_edge_size( Mesh* mesh) {
179		✗	double sum = 0;
180		✗	int nb = 0;
181		✗	FOR(c, mesh->cells.nb()) FOR(lf, mesh->cells.nb_facets(c)) FOR(lv, mesh->cells.facet_nb_vertices(c, lf))
182			{
183		✗	index_t v0 = mesh->cells.facet_vertex(c, lf, lv);
184		✗	index_t v1 = mesh->cells.facet_vertex(c, lf, (lv + 1) % mesh->cells.facet_nb_vertices(c, lf));
185		✗	sum += (mesh->vertices.point(v0) - mesh->vertices.point(v1)).length();
186		✗	nb++;
187			}
188		✗	geo_assert(nb > 0);
189		✗	return sum / double(nb);
190			}
191		✗	double get_facet_average_edge_size( Mesh* m) {
192		✗	geo_assert(m->facet_corners.nb() > 0);
193		✗	double ave_edge_length = 0;
194		✗	FOR(f, m->facets.nb()) FOR(v, m->facets.nb_vertices(f)) ave_edge_length += (X(m)[m->facets.vertex(f, v)] - X(m)[m->facets.vertex(f, (v + 1) % m->facets.nb_vertices(f))]).length();
195		✗	return ave_edge_length / double(m->facet_corners.nb());
196			}
197
198		✗	vec3 tet_facet_cross(Mesh* m, index_t c, index_t lf){
199		✗	vec3 pt[3];
200		✗	for (index_t v = 0; v < 3; v++) pt[v] = m->vertices.point(m->cells.facet_vertex(c, lf, v));
201		✗	return cross(normalize(pt[1] - pt[0]), normalize(pt[2] - pt[0]));
202			}
203
204			/**
205			* HalfedgeToTriangleInTet[cv1][cv2] is the local facet (cf) associated to the halfedge going from cv1 to cv2
206			*/
207			static index_t HalfedgeToTriangleInTet[4][4] = {
208			{NOT_AN_ID, 2, 3, 1},
209			{ 3, NOT_AN_ID, 0, 2 },
210			{ 1, 3, NOT_AN_ID, 0 },
211			{ 2, 0, 1, NOT_AN_ID }
212			};
213
214		✗	index_t next_cell_around_oriented_edge(Mesh* m, index_t cell_id, index_t v_org, index_t v_dest){
215		✗	index_t cv_org = NOT_AN_ID;
216		✗	index_t cv_dest = NOT_AN_ID;
217		✗	FOR(lv, 4) {
218		✗	index_t v = m->cells.vertex(cell_id, lv);
219		✗	if (v == v_org) cv_org = lv;
220		✗	if (v == v_dest) cv_dest = lv;
221			}
222		✗	geo_assert(cv_org != NOT_AN_ID);
223		✗	geo_assert(cv_dest != NOT_AN_ID);
224		✗	geo_assert(cv_org != cv_dest);
225		✗	return m->cells.adjacent(cell_id, HalfedgeToTriangleInTet[cv_org][cv_dest]);
226			}
227
228
229			/* __ __ _ _ _______ _
230			* | \/ | | | (_) |__ __| | |
231			* | \ / | __ _ _ __ ___| |__ _ _ __ __ _ | | ___| |_ ___
232			* | |\/| |/ _` | '__/ __| '_ \| | '_ \ / _` | | |/ _ \ __/ __|
233			* | | | | (_| | | | (__| | | | | | | | (_| | | | __/ |_\__ \
234			* |_| |_|\__,_|_| \___|_| |_|_|_| |_|\__, | |_|\___|\__|___/
235			* __/ |
236			* |___/
237			*/
238
239			const index_t tet_edge_vertices[6][2] = { { 0, 1 }, { 0, 2 }, { 0, 3 }, { 1, 2 }, { 1, 3 }, { 2, 3 } };
240
241			const index_t MTN = index_t(-1);
242			const index_t MT[16][4] = {
243			{MTN, MTN, MTN, MTN}, //0 0 0 0
244			{ 0, 2, 1, MTN },
245			{ 0, 3, 4, MTN },
246			{ 1, 3, 4, 2 },
247			{ 1, 5, 3, MTN }, //0 1 0 0
248			{ 0, 2, 5, 3 },
249			{ 1, 0, 4, 5 },
250			{ 2, 5, 4, MTN },
251			{ 2, 4, 5, MTN }, //1 0 0 0
252			{ 0, 4, 5, 1 },
253			{ 3, 5, 2, 0 },
254			{ 3, 5, 1, MTN },
255			{ 1, 2, 4, 3 }, //1 1 0 0
256			{ 0, 4, 3, MTN },
257			{ 2, 0, 1, MTN },
258			{ MTN,MTN,MTN,MTN }
259			};
260
261
262			}
263