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/preprocessing.h>
41			#include <geogram/mesh/mesh_tetrahedralize.h>
42
43			#include <exploragram/hexdom/quadmesher.h> // for debug output
44
45			namespace GEO {
46
47		✗	static vec3 triangle_normal(const Mesh& M, index_t f){
48		✗	vec3 pt[3];
49		✗	for (index_t v = 0; v < 3; v++) pt[v] = M.vertices.point(M.facets.vertex(f,v));
50		✗	return cross(normalize(pt[1] - pt[0]), normalize(pt[2] - pt[0]));
51			}
52
53		✗	static bool get_a_triangle_patch(
54			const Mesh& M,
55			index_t facet,
56			std::vector<bool>& tri_patch_flag,
57			index_t& patch_size,
58			double cos_angle,
59			index_t nb_tri_min_in_patch) {
60
61			/* Add facet in the patch */
62		✗	tri_patch_flag[facet] = true;
63		✗	patch_size += 1;
64
65			/* Check stop condition */
66		✗	if (patch_size >= nb_tri_min_in_patch) {
67			return true;
68			}
69
70			/* Recursive call */
71		✗	vec3 n = normalize(triangle_normal(M, facet));
72		✗	for (index_t le = 0; le < 3; ++le) {
73			index_t a = M.facets.adjacent(facet, le);
74		✗	if (a == GEO::NO_FACET) continue;
75		✗	vec3 n_a = normalize(triangle_normal(M, a));
76		✗	if (dot(n, n_a) > cos_angle) {
77		✗	if (!tri_patch_flag[a])
78		✗	get_a_triangle_patch(M, a, tri_patch_flag, patch_size, cos_angle, nb_tri_min_in_patch);
79			}
80			}
81		✗	return patch_size >= nb_tri_min_in_patch;
82			}
83
84		✗	static void generate_facet_is_in_patch_attribute(const Mesh& M) {
85		✗	geo_assert(M.facets.nb() > 0);
86
87			/* The input constrained are only computed from "valid" triangles
88			* - triangles are considered valid if they are in a smooth patch of at
89			* least nb_tri_min_in_patch triangles
90			* - two adjacent triangles are in the same patch if the dot product of their
91			* normals is superior to cos_angle */
92
93			/* Parameters */
94			index_t nb_tri_min_in_patch = 8;
95			const double cos_angle = 0.975;
96
97		✗	if (M.facets.nb() < 100) nb_tri_min_in_patch = 4; /* for small models */
98
99		✗	Attribute<int> is_in_patch(M.facets.attributes(), "is_valid");
100			is_in_patch.fill(0);
101
102		✗	std::vector<bool> is_flagged(M.facets.nb());
103		✗	for(index_t f = 0; f < M.facets.nb(); ++f) {
104		✗	if (is_in_patch[f]) continue;
105		✗	index_t patch_size = 0;
106		✗	std::fill(is_flagged.begin(), is_flagged.end(), false); /* reinitialize to zero */
107		✗	bool ok = get_a_triangle_patch(M, f, is_flagged, patch_size, cos_angle, nb_tri_min_in_patch);
108		✗	if (ok) { /* A patch starting at f has been found */
109		✗	for (index_t i = 0; i < is_flagged.size(); ++i) {
110			/* triangles in the patch are keep for imposing constraints */
111		✗	if (is_flagged[i]) is_in_patch[i] = 1;
112			}
113			}
114			}
115		✗	}
116
117		✗	static void compute_input_constraints(Mesh* m, bool relaxed = false) {
118
119		✗	Attribute<mat3> B(m->vertices.attributes(), "B");
120		✗	Attribute<vec3> lockB(m->vertices.attributes(), "lockB");// how many vectors are locked
121		✗	Attribute<vec3> U(m->vertices.attributes(), "U");
122		✗	Attribute<vec3> lockU(m->vertices.attributes(), "lockU");// how many dimensions are locked
123
124			// init all normal for each vertex
125			vector<vector<vec3> > normals(m->vertices.nb());
126			vector<vector<double> > weight(m->vertices.nb());
127
128			/* Compute the facet normals and store them at vertices */
129		✗	if (!relaxed) {
130		✗	FOR(c, m->cells.nb()) FOR(cf, 4) {
131		✗	if ((m->cells.adjacent(c, cf) != NO_CELL)) continue;
132		✗	vec3 n = tet_facet_cross(m, c, cf);
133		✗	if (n.length2() > 1e-10) {
134		✗	FOR(cfv, 3) {
135		✗	normals[m->cells.facet_vertex(c, cf, cfv)].push_back(normalize(n));
136		✗	weight[m->cells.facet_vertex(c, cf, cfv)].push_back(n.length());
137			}
138			}
139			}
140			} else {
141			/* The relaxation of constraints is achieved by storing only the normals of "valid"
142			* triangles.
143			* They are flagged via the facet attribute is_valid
144			* One possibility for flag them is to use generate_facet_is_in_patch_attribute()
145			* See the method for tweaking the parameters. */
146		✗	m->cells.compute_borders();
147		✗	generate_facet_is_in_patch_attribute(*m);
148		✗	Attribute<int> is_valid(m->facets.attributes(), "is_valid");
149
150		✗	for(index_t f = 0; f < m->facets.nb(); ++f) {
151		✗	if (!is_valid[f]) continue;
152		✗	vec3 n = facet_normal(m, f);
153		✗	if (n.length2() > 1e-10) {
154		✗	FOR(lv, 3) {
155		✗	normals[m->facets.vertex(f, lv)].push_back(normalize(n));
156		✗	weight[m->facets.vertex(f, lv)].push_back(n.length());
157			}
158			}
159			}
160			}
161			/* Build the constraints */
162		✗	FOR(v, m->vertices.nb()) {
163			vector<vec3>& n = normals[v];
164			vector<double>& w = weight[v];
165
166			B[v].load_identity();
167		✗	lockB[v] = vec3(0, 0, 0);
168		✗	lockU[v] = vec3(0, 0, 0);
169		✗	if (n.size() > 0) {
170		✗	B[v] = Frame::representative_frame(n, w);// rot_to_B(representative_frame(n, w));
171			AxisPermutation ap;
172		✗	ap.make_col2_equal_to_z(B[v], n[0]);
173		✗	B[v] = Frame(B[v]).apply_permutation(ap);
174		✗	FOR(i, n.size()) FOR(a, 3)
175		✗	if (std::abs(dot(col(B[v], a), n[i])) > .7) {
176		✗	lockU[v][a] = 1;
177		✗	lockB[v][a] = 1;
178			}
179		✗	if (lockB[v].length2() == 1) lockB[v] = vec3(0, 0, 1); // it may not always be true (happened once)
180		✗	if (lockB[v].length2() > 1) lockB[v] = vec3(1, 1, 1);
181		✗	U[v] = vec3(0, 0, 0);
182			}
183
184			//FOR(i,3) FOR(j,3) B[v](i,j) *= 5.;
185			}
186
187			if (relaxed) {
188			// m->facets.clear(false); // Keep the surface part for debugging
189			}
190		✗	}
191
192		✗	static void reorder_vertices_according_to_constraints(Mesh* m,bool hibert_sort) {
193		✗	Attribute<vec3> lockB(m->vertices.attributes(), "lockB");// how many vectors are locked
194
195			GEO::vector<index_t > ind_map(m->vertices.nb());
196			index_t n = 0;
197			index_t num_l_v = 0;
198			index_t num_ln_v = 0;
199		✗	FOR(v, m->vertices.nb()) if (lockB[v][0] == 1) { geo_assert(lockB[v][1] == 1); geo_assert(lockB[v][2] == 1); ind_map[n] = v; n++; }
200			num_l_v = n;
201		✗	FOR(v, m->vertices.nb()) if (lockB[v][0] == 0 && lockB[v][2] == 1) { geo_assert(lockB[v][1] == 0); ind_map[n] = v; n++; }
202			num_ln_v = n;
203		✗	FOR(v, m->vertices.nb()) if (lockB[v][2] == 0) { ind_map[n] = v; n++; }
204
205		✗	for (index_t i = 0; i < num_l_v; i++) geo_assert(lockB[ind_map[i]][0] == 1);
206		✗	for (index_t i = num_l_v; i < num_ln_v; i++) geo_assert(lockB[ind_map[i]][2] == 1);
207		✗	for (index_t i = num_ln_v; i < m->vertices.nb(); i++) geo_assert(lockB[ind_map[i]][2] == 0);
208
209		✗	geo_assert(num_l_v <= num_ln_v && num_ln_v <= m->vertices.nb());
210
211		✗	plop(hibert_sort);
212		✗	if (hibert_sort) {
213		✗	compute_Hilbert_order(m->vertices.nb(), m->vertices.point_ptr(0), ind_map, 0, num_l_v,3);
214		✗	compute_Hilbert_order(m->vertices.nb(), m->vertices.point_ptr(0), ind_map, num_l_v, num_ln_v, 3);
215		✗	compute_Hilbert_order(m->vertices.nb(), m->vertices.point_ptr(0), ind_map, num_ln_v, m->vertices.nb(), 3);
216			}
217		✗	m->vertices.permute_elements(ind_map); // note: it also updates the cell_corners.vertex... and invert ind_map :(
218		✗	}
219
220
221		✗	void produce_hexdom_input(Mesh* m,std::string& error_msg,bool hilbert_sort, bool relaxed) {
222
223		✗	m->edges.clear();
224		✗	m->vertices.remove_isolated();
225
226		✗	if (m->cells.nb() == 0) {
227		✗	if (m->facets.nb() == 0) throw ("mesh have no cells and no facets");
228			mesh_tetrahedralize(*m, true, true, .8);
229			}
230
231		✗	if (have_negative_tet_volume(m)) {
232		✗	throw ("contains tets with negative volume");
233			}
234
235		✗	if (!m->cells.are_simplices()) {
236		✗	throw ("cells contains non tet elements");
237			}
238
239		✗	if (!volume_boundary_is_manifold(m, error_msg)) {
240		✗	throw (error_msg.c_str());
241			}
242
243		✗	if (!volume_is_tetgenifiable(m)) {
244		✗	throw (" tetgen is not able to remesh the volume from its boundary");
245			}
246
247
248
249			// add some attributes
250		✗	compute_input_constraints(m, relaxed);
251
252
253			// compute scale
254		✗	double wanted_edge_length = get_cell_average_edge_size(m);
255
256		✗	Attribute<mat3> B(m->vertices.attributes(), "B");
257		✗	FOR(v, m->vertices.nb()) FOR(ij, 9) B[v].data()[ij] *= wanted_edge_length;
258		✗	reorder_vertices_according_to_constraints(m,hilbert_sort );
259
260		✗	}
261
262
263			}
264