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/extra_connectivity.h>
41			#include <exploragram/hexdom/mesh_utils.h>
42
43			namespace GEO {
44
45		✗	FacetsExtraConnectivity::FacetsExtraConnectivity(Mesh * p_m){
46		✗	m = p_m;
47		✗	reset();
48		✗	}
49
50		✗	void FacetsExtraConnectivity::reset(){
51		✗	index_t nbc = m->facet_corners.nb();
52		✗	c2f.resize(nbc);
53		✗	c2c.resize(nbc);
54		✗	v2c.resize(m->vertices.nb());
55		✗	FOR(f, m->facets.nb()) FOR(fc, m->facets.nb_corners(f)){
56		✗	index_t c = m->facets.corner(f, fc);
57		✗	c2f[c] = f;
58		✗	c2c[c] = c;
59		✗	v2c[m->facets.vertex(f, fc)] = c;
60			}
61		✗	FOR(f, m->facets.nb()) FOR(fc, m->facets.nb_corners(f)){
62		✗	index_t c = m->facets.corner(f, fc);
63		✗	c2c[c] = v2c[m->facets.vertex(f, fc)];
64		✗	v2c[m->facets.vertex(f, fc)] = c;
65			}
66		✗	}
67
68		✗	index_t FacetsExtraConnectivity::org(index_t corner_id){ return m->facet_corners.vertex(corner_id); }
69		✗	index_t FacetsExtraConnectivity::dest(index_t corner_id){ return m->facet_corners.vertex(next(corner_id)); }
70
71		✗	index_t FacetsExtraConnectivity::opposite(index_t corner_id){
72		✗	index_t cir = corner_id;
73		✗	index_t result = NOT_AN_ID; // not found
74			do {
75		✗	index_t candidate = prev(cir);
76		✗	if ((org(candidate) == dest(corner_id)) && (dest(candidate) == org(corner_id))){
77		✗	if (result == NOT_AN_ID) result = candidate;
78		✗	else return NOT_AN_ID; // found more than one
79			}
80		✗	if (cir != corner_id && dest(corner_id) == dest(cir))
81		✗	return NOT_AN_ID; // the edge is non manifold
82		✗	cir = c2c[cir];
83		✗	} while (cir != corner_id);
84		✗	return result;
85			}
86		✗	index_t FacetsExtraConnectivity::next_around_vertex(index_t cir) { return opposite(prev(cir)); }
87
88		✗	index_t FacetsExtraConnectivity::facet(index_t corner_id) { return c2f[corner_id]; }
89		✗	index_t FacetsExtraConnectivity::local_id(index_t corner_id) { return corner_id - m->facets.corners_begin(c2f[corner_id]); }
90
91		✗	index_t FacetsExtraConnectivity::next(index_t corner_id) {
92		✗	index_t fc = local_id(corner_id);
93		✗	index_t offset = corner_id - fc;
94		✗	return offset + next_mod(fc, m->facets.nb_corners(c2f[corner_id]));
95			}
96		✗	index_t FacetsExtraConnectivity::prev(index_t corner_id) {
97		✗	index_t fc = local_id(corner_id);
98		✗	index_t offset = corner_id - fc;
99		✗	return offset + prev_mod(fc, m->facets.nb_corners(c2f[corner_id]));
100			}
101
102		✗	vec3 FacetsExtraConnectivity::geom(index_t corner_id){
103		✗	return X(m)[dest(corner_id)] - X(m)[org(corner_id)];
104			}
105
106
107
108
109
110
111
112
113
114		✗	FacetsExtraConnectivityWithInvalidFacets::FacetsExtraConnectivityWithInvalidFacets(Mesh * p_m) {
115		✗	m = p_m;
116		✗	facet_is_valid.bind(m->facets.attributes(), "is_valid");
117		✗	FOR(f, m->facets.nb()) facet_is_valid[f] = true;
118		✗	reset();
119		✗	}
120
121		✗	void FacetsExtraConnectivityWithInvalidFacets::reset() {
122			//plop(m->facets.nb());
123		✗	index_t nbc = m->facet_corners.nb();
124		✗	c2f.resize(nbc,NOT_AN_ID);
125		✗	c2c.resize(nbc, NOT_AN_ID);
126		✗	v2c.resize(m->vertices.nb(), NOT_AN_ID);
127		✗	FOR(f, m->facets.nb()) FOR(fc, m->facets.nb_corners(f)) {
128		✗	index_t c = m->facets.corner(f, fc);
129		✗	if (facet_is_valid[f]) { // everything is NOT_AN_ID for invalid facet and associated corner
130		✗	c2f[c] = f;
131		✗	c2c[c] = c;
132		✗	v2c[m->facets.vertex(f, fc)] = c;
133			}
134			}
135		✗	FOR(f, m->facets.nb()) {
136		✗	if (!facet_is_valid[f]) continue;
137		✗	FOR(fc, m->facets.nb_corners(f)) {
138		✗	index_t c = m->facets.corner(f, fc);
139		✗	c2c[c] = v2c[m->facets.vertex(f, fc)];
140		✗	v2c[m->facets.vertex(f, fc)] = c;
141			}
142			}
143		✗	}
144
145		✗	index_t FacetsExtraConnectivityWithInvalidFacets::org(index_t corner_id) { return m->facet_corners.vertex(corner_id); }
146		✗	index_t FacetsExtraConnectivityWithInvalidFacets::dest(index_t corner_id) { return m->facet_corners.vertex(next(corner_id)); }
147
148		✗	index_t FacetsExtraConnectivityWithInvalidFacets::opposite(index_t corner_id) {
149		✗	index_t cir = corner_id;
150		✗	index_t result = NOT_AN_ID; // not found
151			do {
152		✗	index_t candidate = prev(cir);
153		✗	if ((org(candidate) == dest(corner_id)) && (dest(candidate) == org(corner_id))) {
154		✗	if (result == NOT_AN_ID) result = candidate;
155		✗	else return NOT_AN_ID; // found more than one
156			}
157		✗	if (cir != corner_id && dest(corner_id) == dest(cir))
158		✗	return NOT_AN_ID; // the edge is non manifold
159		✗	cir = c2c[cir];
160		✗	} while (cir != corner_id);
161		✗	return result;
162			}
163			//index_t FacetsExtraConnectivityWithInvalidFacets::next_around_vertex(index_t cir) { return opposite(prev(cir)); }
164
165		✗	index_t FacetsExtraConnectivityWithInvalidFacets::facet(index_t corner_id) { return c2f[corner_id]; }
166		✗	index_t FacetsExtraConnectivityWithInvalidFacets::local_id(index_t corner_id) { return corner_id - m->facets.corners_begin(c2f[corner_id]); }
167
168		✗	index_t FacetsExtraConnectivityWithInvalidFacets::next(index_t corner_id) {
169		✗	index_t fc = local_id(corner_id);
170		✗	index_t offset = corner_id - fc;
171		✗	return offset + next_mod(fc, m->facets.nb_corners(c2f[corner_id]));
172			}
173		✗	index_t FacetsExtraConnectivityWithInvalidFacets::prev(index_t corner_id) {
174		✗	index_t fc = local_id(corner_id);
175		✗	index_t offset = corner_id - fc;
176		✗	return offset + prev_mod(fc, m->facets.nb_corners(c2f[corner_id]));
177			}
178
179		✗	vec3 FacetsExtraConnectivityWithInvalidFacets::geom(index_t corner_id) {
180		✗	return X(m)[dest(corner_id)] - X(m)[org(corner_id)];
181			}
182
183
184
185
186
187
188
189
190
191
192		✗	void halfedge_manip_example(Mesh* m){
193		✗	FacetsExtraConnectivity fec(m);
194		✗	FOR(c, m->facet_corners.nb()){// we can directly loop over each halfedges
195			// turning around a facet
196		✗	index_t cir = c;
197			do {
198		✗	cir = fec.next(cir);
199		✗	} while (cir != c);
200
201			//iterate on vertex incident
202		✗	cir = c;
203			do {
204		✗	cir = fec.c2c[cir];
205		✗	} while (cir != c);
206
207			//turning around a manifold vertex
208		✗	cir = c;
209			do {
210		✗	cir = fec.next_around_vertex(cir);
211		✗	} while (cir != c && cir != NOT_AN_ID);
212			}
213		✗	}
214
215
216
217
218		✗	void create_facet_adjacence(Mesh* m, bool has_border){
219		✗	FacetsExtraConnectivity qfec(m);
220		✗	FOR(h, m->facet_corners.nb()){
221		✗	index_t opp = qfec.opposite(h);
222		✗	if (!has_border && opp == NOT_AN_ID) GEO::Logger::out("HexDom") << "PANIC MODE, INPUT IS NON MANIFOLD !!! --- check if the surface has border" << std::endl;
223			else {
224		✗	m->facets.set_adjacent(qfec.facet(h), qfec.local_id(h), qfec.facet(opp));
225		✗	m->facets.set_adjacent(qfec.facet(opp), qfec.local_id(opp), qfec.facet(h));
226			}
227			}
228		✗	}
229
230
231
232
233
234		✗	void cell_edges_in_RCS(Mesh* m, vector<index_t>& offset_from_org, vector<index_t>& dest){
235		✗	geo_assert(m->cells.are_simplices());
236		✗	dest.clear();
237		✗	offset_from_org.clear();
238		✗	offset_from_org.reserve(m->vertices.nb());
239		✗	vector<index_t> v2cc(m->vertices.nb(), NOT_AN_ID); // v2cc maps vertices to cell corner
240		✗	vector<index_t> next(4 * m->cells.nb(), NOT_AN_ID); // next chains cell corners
241
242		✗	FOR(c, m->cells.nb()) FOR(cv, 4){
243		✗	index_t v = m->cells.vertex(c, cv);
244		✗	next[4 * c + cv] = v2cc[v];
245		✗	v2cc[v] = 4 * c + cv;
246			}
247		✗	FOR(v, m->vertices.nb()){
248		✗	offset_from_org.push_back(dest.size());
249		✗	for (index_t i = v2cc[v]; i != NOT_AN_ID; i = next[i])
250		✗	FOR(lc, 4) if (i % 4 != lc)
251		✗	dest.push_back(m->cells.vertex(i / 4, lc));
252		✗	std::sort(dest.begin() + int(offset_from_org.back()), dest.end());
253		✗	vector<index_t>::iterator last = std::unique(dest.begin() + int(offset_from_org.back()), dest.end());
254		✗	dest.resize(size_t(last - dest.begin()));
255			}
256		✗	offset_from_org.push_back(dest.size());
257		✗	}
258
259
260
261
262		✗	void compute_tet_edge_graph(Mesh* m,vector<index_t> & v2e,bool store_both_directions){
263		✗	geo_assert(m->cells.are_simplices());
264		✗	m->edges.clear();
265		✗	v2e.resize(m->vertices.nb());
266
267			//Attribute<index_t> v2e(m->vertices.attributes(), "v2e");
268
269		✗	vector<index_t> adj_;
270		✗	vector<index_t> adj_off_;
271		✗	adj_off_.reserve(m->vertices.nb() + 1);
272
273		✗	vector<index_t> v2c(m->vertices.nb(), NOT_AN_ID);
274		✗	vector<index_t> next(4 * m->cells.nb(), NOT_AN_ID);
275
276		✗	FOR(c, m->cells.nb()) FOR(cv, 4){
277		✗	index_t v = m->cells.vertex(c, cv);
278		✗	next[4 * c + cv] = v2c[v];
279		✗	v2c[v] = 4 * c + cv;
280			}
281		✗	FOR(v, m->vertices.nb()){
282		✗	adj_off_.push_back(adj_.size());
283		✗	for (index_t i = v2c[v]; i != NOT_AN_ID; i = next[i]) FOR(lc, 4) if (i % 4 != lc){
284		✗	index_t nv = m->cells.vertex(i / 4, lc);
285		✗	if (nv>v || store_both_directions) adj_.push_back(nv);
286			}
287		✗	std::sort(adj_.begin() + int(adj_off_.back()), adj_.end());
288		✗	vector<index_t>::iterator last = std::unique(adj_.begin() + int(adj_off_.back()), adj_.end());
289		✗	adj_.resize(size_t(last - adj_.begin()));
290			}
291		✗	adj_off_.push_back(adj_.size());
292
293		✗	index_t offe = m->edges.create_edges(adj_.size());
294		✗	FOR(v, m->vertices.nb()){
295		✗	v2e[v] = adj_off_[v];
296		✗	for (index_t e = adj_off_[v]; e < adj_off_[v + 1]; e++){
297		✗	m->edges.set_vertex(offe + e, 0, v);
298		✗	m->edges.set_vertex(offe + e, 1, adj_[e]);
299			}
300			}
301		✗	}
302
303		✗	void restore_v2e(Mesh* m, vector<index_t>& v2e) {
304		✗	v2e.resize(m->vertices.nb());
305		✗	FOR(v, m->vertices.nb())v2e[v] = NOT_AN_ID;
306		✗	index_t lastv = 0;
307		✗	FOR(e, m->edges.nb()) {
308		✗	index_t org = m->edges.vertex(e, 0);
309		✗	while (lastv <= org) { v2e[lastv] = e; lastv++; }
310			}
311		✗	while (lastv < m->vertices.nb()) { v2e[lastv] = m->edges.nb(); lastv++; }
312		✗	}
313
314
315
316
317			/**
318			* check that v2f[v] gives all facets that
319			* are adjacent to v
320			* are not to be killed
321			*/
322
323		✗	bool v2f_is_valid(Mesh* m, vector<vector<index_t> >& v2f, vector<index_t> &to_kill) {
324		✗	vector<vector<bool> > is_ref(m->vertices.nb());
325		✗	FOR(v, m->vertices.nb()) is_ref[v].resize(v2f[v].size(), false);
326			// check that all facets are referenced
327		✗	FOR(f, m->facets.nb()) if (!to_kill[f])
328		✗	FOR(lc, m->facets.nb_corners(f)) {
329		✗	index_t v = m->facets.vertex(f, lc);
330		✗	bool f_is_ref = false;
331		✗	FOR(i, v2f[v].size()) if (v2f[v][i] == f) {
332		✗	is_ref[v][i] = true;
333		✗	f_is_ref = true;
334			}
335		✗	if (!f_is_ref) {
336		✗	GEO::Logger::out("HexDom") << "facet " << f << " is not referenced in vertex " << v << " !!!" << std::endl;
337		✗	return false;
338			}
339			}
340			// check that no extra facet is referenced
341		✗	FOR(v, m->vertices.nb()) FOR(i, v2f[v].size()) if (!is_ref[v][i]) {
342		✗	GEO::Logger::out("HexDom") << "one facet of v2f does not exist!!!" << std::endl;
343		✗	return false;
344			}
345		✗	return true;
346		✗	}
347			}
348