GCC Code Coverage Report

Directory: ./
File: lib/geogram/mesh/mesh.cpp
Date: 2026-09-07 02:25:23
Exec Total Coverage
Lines: 483 1009 47.9%
Functions: 68 95 71.6%
Branches: 281 1159 24.2%
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 <geogram/mesh/mesh.h>
41 #include <geogram/basic/permutation.h>
42 #include <geogram/basic/logger.h>
43 #include <geogram/basic/algorithm.h>
44 #include <geogram/basic/string.h>
45
46 namespace GEO {
47
48 5733 MeshSubElementsStore::MeshSubElementsStore(Mesh& mesh) :
49 5733 mesh_(mesh),
50 5733 nb_(0) {
51 5733 }
52
53 11466 MeshSubElementsStore::~MeshSubElementsStore() {
54 11466 }
55
56 3147 void MeshSubElementsStore::clear_store(
57 bool keep_attributes, bool keep_memory
58 ) {
59 3147 attributes_.clear(keep_attributes, keep_memory);
60 3147 nb_ = 0;
61 3147 }
62
63 1278 void MeshSubElementsStore::resize_store(index_t new_size) {
64 1278 attributes_.resize(new_size);
65 1278 nb_ = new_size;
66 1278 }
67
68 /*************************************************************************/
69
70 3276 MeshElements::MeshElements() {
71 3276 }
72
73 6552 MeshElements::~MeshElements() {
74 6552 }
75
76 /*************************************************************************/
77
78 819 MeshVertices::MeshVertices(Mesh& mesh) :
79 MeshSubElementsStore(mesh),
80 819 edges_(mesh.edges),
81 819 facet_corners_(mesh.facet_corners),
82
1/2
✓ Branch 2 taken 819 times.
✗ Branch 3 not taken.
 819 cell_corners_(mesh.cell_corners) {
83 819 }
84
85
1/2
✓ Branch 0 taken 819 times.
✗ Branch 1 not taken.
 1638 MeshVertices::~MeshVertices() {
86 if(point_.is_bound()) {
87 1638 point_.unbind();
88 }
89 if(point_fp32_.is_bound()) {
90 ✗ point_fp32_.unbind();
91 }
92 1638 }
93
94
95
1/2
✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
 3 void MeshVertices::set_double_precision() {
96 if(double_precision()) {
97 return;
98 }
99
100 index_t dim = dimension();
101
102 ✗ point_.create_vector_attribute(
103 ✗ attributes(), "point", dim
104 );
105
106 ✗ for(index_t i=0; i<point_.nb_elements(); ++i) {
107 ✗ point_[i] = double(point_fp32_[i]);
108 }
109
110 ✗ point_fp32_.destroy();
111 }
112
113 ✗ void MeshVertices::set_single_precision() {
114 if(single_precision()) {
115 return;
116 }
117
118 index_t dim = dimension();
119
120 ✗ point_fp32_.create_vector_attribute(
121 ✗ attributes(), "point_fp32", dim
122 );
123
124 ✗ for(index_t i=0; i<point_.nb_elements(); ++i) {
125 ✗ point_fp32_[i] = float(point_[i]);
126 }
127
128 ✗ point_.destroy();
129 }
130
131
132
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 372 times.
 372 void MeshVertices::clear(bool keep_attributes, bool keep_memory) {
133 bool singlep = single_precision();
134 index_t dim = dimension();
135
136 // We need to unbind point attributes
137 // because it is not correct to clear the
138 // AttributesManager when an attribute is
139 // still bound.
140
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 372 times.
 372 if(!keep_attributes) {
141 if(point_.is_bound()) {
142 ✗ point_.unbind();
143 }
144 if(point_fp32_.is_bound()) {
145 ✗ point_fp32_.unbind();
146 }
147 }
148
149 372 clear_store(keep_attributes, keep_memory);
150
151 // Now we can re-create the point attributes.
152
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 372 times.
 372 if(!keep_attributes) {
153 ✗ bind_point_attribute(dim,singlep);
154 }
155 372 }
156
157 372 void MeshVertices::clear_store(
158 bool keep_attributes, bool keep_memory
159 ) {
160 372 MeshSubElementsStore::clear_store(keep_attributes, keep_memory);
161 372 }
162
163 173 void MeshVertices::resize_store(index_t new_size) {
164 173 MeshSubElementsStore::resize_store(new_size);
165 173 }
166
167 472 void MeshVertices::delete_elements(
168 vector<index_t>& to_delete,
169 bool remove_isolated_vertices
170 ) {
171 // to_delete is reused for re-indexing
172 vector<index_t>& old2new = to_delete;
173 geo_argused(remove_isolated_vertices);
174
2/2
✓ Branch 0 taken 170 times.
✓ Branch 1 taken 302 times.
 472 if(has_non_zero(to_delete)) {
175 index_t cur=0;
176
2/2
✓ Branch 0 taken 229749 times.
✓ Branch 1 taken 170 times.
 229919 for(index_t i=0; i<old2new.size(); ++i) {
177
2/2
✓ Branch 0 taken 142969 times.
✓ Branch 1 taken 86780 times.
 229749 if(old2new[i] == 0) {
178 142969 old2new[i] = cur;
179 142969 ++cur;
180 } else {
181 86780 old2new[i] = NO_INDEX;
182 }
183 }
184 170 attributes_.compress(old2new);
185 // cur now contains the new size.
186 170 resize_store(cur);
187
188
2/2
✓ Branch 0 taken 528 times.
✓ Branch 1 taken 170 times.
 698 for(index_t e=0; e<edges_.nb(); ++e) {
189
2/2
✓ Branch 0 taken 1056 times.
✓ Branch 1 taken 528 times.
 1584 for(index_t lv=0; lv<2; ++lv) {
190 index_t v = edges_.vertex(e,lv);
191 1056 v = old2new[v];
192 edges_.set_vertex(e,lv,v);
193 }
194 }
195
196
2/2
✓ Branch 0 taken 838086 times.
✓ Branch 1 taken 170 times.
 838256 for(index_t c=0; c<facet_corners_.nb(); ++c) {
197 index_t v = facet_corners_.vertex(c);
198 838086 v = old2new[v];
199 facet_corners_.set_vertex(c,v);
200 }
201
202
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 170 times.
 170 for(index_t c=0; c<cell_corners_.nb(); ++c) {
203 index_t v = cell_corners_.vertex(c);
204 // Cells can have padding
205 ✗ if(v == NO_VERTEX) {
206 ✗ continue;
207 }
208 ✗ v = old2new[v];
209 cell_corners_.set_vertex(c,v);
210 }
211 }
212 472 }
213
214 10 void MeshVertices::permute_elements(vector<index_t>& permutation) {
215 10 attributes_.apply_permutation(permutation);
216 10 Permutation::invert(permutation);
217
218
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
 10 for(index_t e=0; e<edges_.nb(); ++e) {
219 ✗ for(index_t lv=0; lv<2; ++lv) {
220 index_t v = edges_.vertex(e,lv);
221 ✗ v = permutation[v];
222 edges_.set_vertex(e,lv,v);
223 }
224 }
225
226
2/2
✓ Branch 0 taken 25698 times.
✓ Branch 1 taken 10 times.
 25708 for(index_t c=0; c<facet_corners_.nb(); ++c) {
227 index_t v = facet_corners_.vertex(c);
228 25698 v = permutation[v];
229 facet_corners_.set_vertex(c,v);
230 }
231
232
2/2
✓ Branch 0 taken 6144 times.
✓ Branch 1 taken 10 times.
 6154 for(index_t c=0; c<cell_corners_.nb(); ++c) {
233 index_t v = cell_corners_.vertex(c);
234 // Cells can have padding
235
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 6144 times.
 6144 if(v == NO_VERTEX) {
236 ✗ continue;
237 }
238 6144 v = permutation[v];
239 cell_corners_.set_vertex(c,v);
240 }
241 10 }
242
243 430 void MeshVertices::remove_isolated() {
244 430 vector<index_t> to_delete(nb(),1);
245
246
2/2
✓ Branch 0 taken 4241 times.
✓ Branch 1 taken 430 times.
 4671 for(index_t e=0; e<mesh_.edges.nb(); ++e) {
247
2/2
✓ Branch 0 taken 8482 times.
✓ Branch 1 taken 4241 times.
 12723 for(index_t lv=0; lv<2; ++lv) {
248 index_t v = mesh_.edges.vertex(e,lv);
249 8482 to_delete[v] = 0;
250 }
251 }
252
253
2/2
✓ Branch 0 taken 749855 times.
✓ Branch 1 taken 430 times.
 750285 for(index_t f=0; f<mesh_.facets.nb(); ++f) {
254 2270025 for(index_t co=mesh_.facets.corners_begin(f);
255
2/2
✓ Branch 0 taken 2270025 times.
✓ Branch 1 taken 749855 times.
 3019880 co<mesh_.facets.corners_end(f); ++co
256 ) {
257 index_t v = mesh_.facet_corners.vertex(co);
258 2270025 to_delete[v] = 0;
259 }
260 }
261
262
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 430 times.
 430 for(index_t c=0; c<mesh_.cells.nb(); ++c) {
263 ✗ for(
264 index_t co=mesh_.cells.corners_begin(c);
265 ✗ co < mesh_.cells.corners_end(c); ++co
266 ) {
267 ✗ index_t v = mesh_.cell_corners.vertex(co);
268 ✗ to_delete[v] = 0;
269 }
270 }
271
272
1/2
✓ Branch 1 taken 430 times.
✗ Branch 2 not taken.
 430 delete_elements(to_delete);
273 430 }
274
275 819 void MeshVertices::bind_point_attribute(
276 index_t dim, bool single_precision
277 ) {
278
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 819 times.
 819 if(single_precision) {
279 ✗ point_fp32_.create_vector_attribute(
280 ✗ attributes(), "point_fp32", dim
281 );
282 } else {
283
1/2
✓ Branch 1 taken 819 times.
✗ Branch 2 not taken.
 819 point_.create_vector_attribute(
284 1638 attributes(), "point", dim
285 );
286 }
287 819 }
288
289
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 15 times.
 15 void MeshVertices::assign_points(
290 vector<double>& points, index_t dim, bool steal_arg
291 ) {
292 // TODO: implement steal_arg
293 geo_argused(steal_arg);
294 15 index_t nb_pts = points.size()/dim;
295
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 15 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 15 geo_assert(dim*nb_pts == points.size());
296 15 assign_points(points.data(), dim, nb_pts);
297 15 }
298
299
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 15 times.
 15 void MeshVertices::assign_points(
300 const double* points, index_t dim, index_t nb_pts
301 ) {
302 ✗ geo_assert(!single_precision());
303
2/4
✓ Branch 0 taken 15 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 15 times.
✗ Branch 3 not taken.
 15 if(dim != dimension() || nb_pts != nb()) {
304 15 clear(true,false);
305 15 set_dimension(dim);
306 create_vertices(nb_pts);
307 }
308 15 Memory::copy(
309 15 point_ptr(0), points, nb_pts*dim*sizeof(double)
310 );
311 15 }
312
313 ✗ void MeshVertices::pop() {
314 geo_debug_assert(nb() != 0);
315 ✗ --nb_;
316 ✗ }
317
318 /**************************************************************************/
319
320
1/2
✓ Branch 2 taken 819 times.
✗ Branch 3 not taken.
 819 MeshEdges::MeshEdges(Mesh& mesh) : MeshSubElementsStore(mesh) {
321 819 }
322
323
2/2
✓ Branch 0 taken 254 times.
✓ Branch 1 taken 565 times.
 1638 MeshEdges::~MeshEdges() {
324 1638 }
325
326 3 void MeshEdges::delete_elements(
327 vector<index_t>& to_delete, bool remove_isolated_vertices
328 ) {
329 geo_debug_assert(to_delete.size() == nb());
330
331 // "Fast track" if no element should be deleted
332
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
 3 if(!has_non_zero(to_delete)) {
333 ✗ if(remove_isolated_vertices) {
334 ✗ mesh_.vertices.remove_isolated();
335 }
336 ✗ return;
337 }
338
339 // to_delete is used for both indicating
340 // which edges should be deleted and
341 // for storing the re-numbering map
342 vector<index_t>& edges_old2new = to_delete;
343 index_t new_nb_edges = 0;
344
345
2/2
✓ Branch 0 taken 64 times.
✓ Branch 1 taken 3 times.
 67 for(index_t e = 0; e < nb(); ++e) {
346
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 58 times.
 64 if(edges_old2new[e] != 0) {
347 6 edges_old2new[e] = NO_EDGE;
348 } else {
349 58 edges_old2new[e] = new_nb_edges;
350
2/2
✓ Branch 0 taken 20 times.
✓ Branch 1 taken 38 times.
 58 if(new_nb_edges != e) {
351 20 edge_vertex_[2*new_nb_edges] = edge_vertex_[2*e];
352 20 edge_vertex_[2*new_nb_edges+1] = edge_vertex_[2*e+1];
353 }
354 58 ++new_nb_edges;
355 }
356 }
357
358 // Manage facets store and attributes
359 3 attributes().compress(edges_old2new);
360 3 resize_store(new_nb_edges);
361
362
1/2
✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
 3 if(remove_isolated_vertices) {
363 3 mesh_.vertices.remove_isolated();
364 }
365 }
366
367 ✗ void MeshEdges::permute_elements(vector<index_t>& permutation) {
368 ✗ attributes_.apply_permutation(permutation);
369 ✗ Permutation::apply(
370 edge_vertex_.data(),
371 permutation,
372 index_t(sizeof(index_t) * 2)
373 );
374 ✗ }
375
376 624 void MeshEdges::clear(bool keep_attributes, bool keep_memory) {
377 624 clear_store(keep_attributes, keep_memory);
378 624 }
379
380 ✗ void MeshEdges::pop() {
381 geo_debug_assert(nb() != 0);
382 ✗ resize_store(nb()-1);
383 ✗ }
384
385 ✗ void MeshEdges::flip(index_t e) {
386 geo_debug_assert(e < nb());
387 index_t v1 = vertex(e,0);
388 index_t v2 = vertex(e,1);
389 set_vertex(e,0,v2);
390 set_vertex(e,1,v1);
391 ✗ }
392
393 624 void MeshEdges::clear_store(
394 bool keep_attributes, bool keep_memory
395 ) {
396
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 624 times.
 624 if(keep_memory) {
397 ✗ edge_vertex_.resize(0);
398 } else {
399 edge_vertex_.clear();
400 }
401 624 MeshSubElementsStore::clear_store(keep_attributes, keep_memory);
402 624 }
403
404 3 void MeshEdges::resize_store(index_t new_size) {
405 3 edge_vertex_.resize(new_size*2,NO_VERTEX);
406 3 MeshSubElementsStore::resize_store(new_size);
407 3 }
408
409
410 /**************************************************************************/
411
412 819 MeshFacetsStore::MeshFacetsStore(Mesh& mesh) :
413 MeshSubElementsStore(mesh),
414
1/2
✓ Branch 2 taken 819 times.
✗ Branch 3 not taken.
 819 is_simplicial_(true) {
415
1/4
✓ Branch 1 taken 819 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 819 facet_ptr_.push_back(0);
416 819 }
417
418 540 void MeshFacetsStore::clear_store(
419 bool keep_attributes, bool keep_memory
420 ) {
421
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 540 times.
 540 if(keep_memory) {
422 ✗ facet_ptr_.resize(0);
423 } else {
424 facet_ptr_.clear();
425 }
426 540 MeshSubElementsStore::clear_store(keep_attributes, keep_memory);
427 540 facet_ptr_.push_back(0);
428 540 }
429
430 543 void MeshFacetsStore::resize_store(index_t new_size) {
431
2/2
✓ Branch 0 taken 11 times.
✓ Branch 1 taken 532 times.
 543 if(!is_simplicial_) {
432 11 facet_ptr_.resize(new_size+1);
433 }
434 543 MeshSubElementsStore::resize_store(new_size);
435 543 }
436
437 /**************************************************************************/
438
439 819 MeshFacetCornersStore::MeshFacetCornersStore(Mesh& mesh) :
440 MeshSubElementsStore(mesh),
441 819 vertices_(mesh.vertices),
442 819 facets_(mesh.facets) {
443 819 }
444
445 540 void MeshFacetCornersStore::clear_store(
446 bool keep_attributes, bool keep_memory
447 ) {
448
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 540 times.
 540 if(keep_memory) {
449 ✗ corner_vertex_.resize(0);
450 ✗ corner_adjacent_facet_.resize(0);
451 } else {
452 corner_vertex_.clear();
453 corner_adjacent_facet_.clear();
454 }
455 540 MeshSubElementsStore::clear_store(keep_attributes, keep_memory);
456 540 }
457
458 543 void MeshFacetCornersStore::resize_store(index_t new_size) {
459 543 corner_vertex_.resize(new_size);
460 543 corner_adjacent_facet_.resize(new_size);
461 543 MeshSubElementsStore::resize_store(new_size);
462 543 }
463
464 /**************************************************************************/
465
466 819 MeshFacets::MeshFacets(Mesh& mesh) :
467 MeshFacetsStore(mesh),
468 819 vertices_(mesh.vertices),
469
1/2
✓ Branch 2 taken 819 times.
✗ Branch 3 not taken.
 819 facet_corners_(mesh.facet_corners) {
470 819 }
471
472 540 void MeshFacets::clear(bool keep_attributes, bool keep_memory) {
473 540 facet_corners_.clear_store(keep_attributes, keep_memory);
474 540 clear_store(keep_attributes, keep_memory);
475 is_simplicial();
476 540 }
477
478 287 void MeshFacets::delete_elements(
479 vector<index_t>& to_delete,
480 bool remove_isolated_vertices
481 ) {
482 geo_debug_assert(to_delete.size() == nb());
483
484 // "Fast track" if no element should be deleted
485
2/2
✓ Branch 0 taken 77 times.
✓ Branch 1 taken 210 times.
 287 if(!has_non_zero(to_delete)) {
486
1/2
✓ Branch 0 taken 77 times.
✗ Branch 1 not taken.
 77 if(remove_isolated_vertices) {
487 77 mesh_.vertices.remove_isolated();
488 }
489 77 return;
490 }
491
492 // to_delete is used for both indicating
493 // which facets should be deleted and
494 // for storing the re-numbering map
495 vector<index_t>& facets_old2new = to_delete;
496
497
2/2
✓ Branch 0 taken 115 times.
✓ Branch 1 taken 95 times.
 210 vector<index_t>& corner_vertex = facet_corners_.corner_vertex_;
498 vector<index_t>& corner_adjacent_facet =
499 facet_corners_.corner_adjacent_facet_;
500
501 index_t new_nb_facets = 0;
502 index_t new_nb_corners = 0;
503
504 // If there are some corner attributes, we need
505 // to compute the index mapping for them.
506 vector<index_t> corners_old2new;
507
2/2
✓ Branch 0 taken 115 times.
✓ Branch 1 taken 95 times.
 210 if(facet_corners_.attributes().nb() != 0) {
508
1/2
✓ Branch 1 taken 115 times.
✗ Branch 2 not taken.
 115 corners_old2new.resize(facet_corners_.nb(), NO_INDEX);
509 }
510
511
2/2
✓ Branch 0 taken 923711 times.
✓ Branch 1 taken 210 times.
 923921 for(index_t f = 0; f < nb(); ++f) {
512
2/2
✓ Branch 0 taken 378866 times.
✓ Branch 1 taken 544845 times.
 923711 if(facets_old2new[f] != 0) {
513 378866 facets_old2new[f] = NO_FACET;
514 } else {
515 544845 facets_old2new[f] = new_nb_facets;
516
2/2
✓ Branch 0 taken 20742 times.
✓ Branch 1 taken 524103 times.
 544845 if(!is_simplicial_) {
517 20742 facet_ptr_[new_nb_facets] = new_nb_corners;
518 }
519
2/2
✓ Branch 0 taken 1654995 times.
✓ Branch 1 taken 544845 times.
 2199840 for(index_t co = corners_begin(f); co != corners_end(f); ++co) {
520
2/2
✓ Branch 0 taken 652164 times.
✓ Branch 1 taken 1002831 times.
 1654995 if(corners_old2new.size() != 0) {
521 652164 corners_old2new[co] = new_nb_corners;
522 }
523
2/2
✓ Branch 0 taken 1425846 times.
✓ Branch 1 taken 229149 times.
 1654995 if(co != new_nb_corners) {
524 1425846 corner_vertex[new_nb_corners] =
525 corner_vertex[co];
526 1425846 corner_adjacent_facet[new_nb_corners] =
527 corner_adjacent_facet[co];
528 }
529 1654995 new_nb_corners++;
530 }
531 544845 new_nb_facets++;
532 }
533 }
534
535
2/2
✓ Branch 0 taken 11 times.
✓ Branch 1 taken 199 times.
 210 if(!is_simplicial_) {
536 11 facet_ptr_[new_nb_facets] = new_nb_corners;
537 }
538
539 // Map adjacent facets indices
540
2/2
✓ Branch 0 taken 2796832 times.
✓ Branch 1 taken 210 times.
 2797042 for(index_t c = 0; c < facet_corners_.nb(); ++c) {
541 2796832 index_t f = corner_adjacent_facet[c];
542
2/2
✓ Branch 0 taken 1467021 times.
✓ Branch 1 taken 1329811 times.
 2796832 if(f != NO_FACET) {
543 1467021 corner_adjacent_facet[c] = facets_old2new[f];
544 }
545 }
546
547 // Manage facets store and attributes
548
1/2
✓ Branch 1 taken 210 times.
✗ Branch 2 not taken.
 210 attributes().compress(facets_old2new);
549
1/2
✓ Branch 1 taken 210 times.
✗ Branch 2 not taken.
 210 resize_store(new_nb_facets);
550
551 // Manage corners store and attributes
552
2/2
✓ Branch 0 taken 115 times.
✓ Branch 1 taken 95 times.
 210 if(corners_old2new.size() != 0) {
553 // corners index mapping is computed only if there
554 // were some corner attributes.
555
1/2
✓ Branch 1 taken 115 times.
✗ Branch 2 not taken.
 115 facet_corners_.attributes().compress(corners_old2new);
556 }
557
1/2
✓ Branch 1 taken 210 times.
✗ Branch 2 not taken.
 210 facet_corners_.resize_store(new_nb_corners);
558
559
1/2
✓ Branch 0 taken 210 times.
✗ Branch 1 not taken.
 210 if(remove_isolated_vertices) {
560
1/2
✓ Branch 1 taken 210 times.
✗ Branch 2 not taken.
 210 mesh_.vertices.remove_isolated();
561 }
562 }
563
564 10 void MeshFacets::permute_elements(vector<index_t>& permutation) {
565 10 attributes_.apply_permutation(permutation);
566
567
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 6 times.
 10 vector<index_t>& corner_vertex = facet_corners_.corner_vertex_;
568 vector<index_t>& corner_adjacent_facet =
569 facet_corners_.corner_adjacent_facet_;
570
571
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 6 times.
 10 if(facet_corners_.attributes().nb() != 0) {
572 vector<index_t> facet_corners_permutation;
573
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 facet_corners_permutation.reserve(facet_corners_.nb());
574
575
2/2
✓ Branch 0 taken 4352 times.
✓ Branch 1 taken 4 times.
 4356 for(index_t new_f=0; new_f<nb(); ++new_f) {
576
1/2
✓ Branch 0 taken 4352 times.
✗ Branch 1 not taken.
 4352 index_t old_f = permutation[new_f];
577 13056 for(
578 4352 index_t old_c=corners_begin(old_f);
579
2/2
✓ Branch 0 taken 13056 times.
✓ Branch 1 taken 4352 times.
 17408 old_c<corners_end(old_f); ++old_c) {
580 facet_corners_permutation.push_back(old_c);
581 }
582 }
583
584
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 facet_corners_.attributes().apply_permutation(
585 facet_corners_permutation
586 );
587 }
588
589
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
 10 if(is_simplicial_) {
590 // If the surface is triangulated,
591 // everything can be done in-place (great !!)
592
593 10 Permutation::apply(
594 corner_vertex.data(),
595 permutation,
596 index_t(sizeof(index_t) * 3)
597 );
598
599 10 Permutation::apply(
600 corner_adjacent_facet.data(),
601 permutation,
602 index_t(sizeof(index_t) * 3)
603 );
604
605 10 Permutation::invert(permutation);
606
607
2/2
✓ Branch 0 taken 25698 times.
✓ Branch 1 taken 10 times.
 51416 for(index_t c = 0; c < corner_adjacent_facet.size(); ++c) {
608
2/2
✓ Branch 0 taken 25520 times.
✓ Branch 1 taken 178 times.
 25698 if(corner_adjacent_facet[c] != NO_FACET) {
609 25520 corner_adjacent_facet[c] =
610 permutation[corner_adjacent_facet[c]];
611 }
612 }
613
614 } else {
615
616 {
617 vector<index_t> new_corner_vertex;
618 ✗ new_corner_vertex.reserve(corner_vertex.size());
619 vector<index_t> new_corner_adjacent_facet;
620 ✗ new_corner_adjacent_facet.reserve(corner_adjacent_facet.size());
621 vector<index_t> new_facet_ptr;
622 ✗ new_facet_ptr.reserve(nb()+1);
623
624 ✗ new_facet_ptr.push_back(0);
625 ✗ for(index_t new_f=0; new_f<nb(); ++new_f) {
626 ✗ index_t old_f = permutation[new_f];
627 ✗ for(
628 index_t old_c = corners_begin(old_f);
629 ✗ old_c < corners_end(old_f); ++old_c
630 ) {
631 new_corner_vertex.push_back(
632 ✗ mesh_.facet_corners.vertex(old_c)
633 );
634 new_corner_adjacent_facet.push_back(
635 ✗ mesh_.facet_corners.adjacent_facet(old_c)
636 );
637 }
638 new_facet_ptr.push_back(
639 ✗ new_facet_ptr[new_facet_ptr.size()-1] +
640 nb_vertices(old_f)
641 );
642 }
643
644 corner_vertex.swap(new_corner_vertex);
645 corner_adjacent_facet.swap(new_corner_adjacent_facet);
646 facet_ptr_.swap(new_facet_ptr);
647 }
648
649
650 ✗ Permutation::invert(permutation);
651
652 ✗ for(index_t c = 0; c < corner_adjacent_facet.size(); ++c) {
653 ✗ if(corner_adjacent_facet[c] != NO_FACET) {
654 ✗ corner_adjacent_facet[c] =
655 permutation[corner_adjacent_facet[c]];
656 }
657 }
658 }
659 10 }
660
661 741 void MeshFacets::connect() {
662 741 connect(0, nb());
663 741 }
664
665 758 void MeshFacets::connect(index_t f_begin, index_t f_end) {
666
667
2/2
✓ Branch 0 taken 40 times.
✓ Branch 1 taken 718 times.
 758 if(f_begin == f_end) {
668 40 return;
669 }
670
671 // Sanity check: no facet is incident to same vertex
672 // several times
673 #ifdef GEO_DEBUG
674 {
675 for(index_t f = f_begin; f != f_end; ++f) {
676 for(index_t lv1=0; lv1<nb_vertices(f); ++lv1) {
677 for(index_t lv2=lv1+1; lv2<nb_vertices(f); ++lv2) {
678 geo_debug_assert(vertex(f,lv1) != vertex(f,lv2));
679 }
680 }
681 }
682 }
683 #endif
684
685 // Get facet corners slice
686 index_t c_begin = corners_begin(f_begin);
687 index_t c_end = corners_end(f_end-1);
688
689 // Get vertices slices indexed by facets in slice
690 index_t v_begin = NO_INDEX;
691 index_t v_end = NO_INDEX;
692
693
2/4
✓ Branch 0 taken 718 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 718 times.
✗ Branch 3 not taken.
 718 if(c_begin == 0 && c_end == facet_corners_.nb()) {
694 718 v_begin = 0;
695 718 v_end = vertices_.nb();
696 } else {
697 ✗ v_begin = facet_corners_.vertex(c_begin);
698 ✗ v_end = v_begin;
699 ✗ for(index_t c=c_begin; c!=c_end; ++c) {
700 index_t v = facet_corners_.vertex(c);
701 ✗ v_begin = std::min(v_begin,v);
702 ✗ v_end = std::max(v_end,v);
703 }
704 ✗ ++v_end;
705 }
706
707 // Gives for each corner the facet incident to it
708 // (or use c/3 if the surface is triangulated).
709 vector<index_t> c2f;
710
2/2
✓ Branch 0 taken 149 times.
✓ Branch 1 taken 569 times.
 718 if(!is_simplicial_) {
711
1/2
✓ Branch 1 taken 149 times.
✗ Branch 2 not taken.
 149 c2f.assign(c_end - c_begin, NO_FACET);
712
2/2
✓ Branch 0 taken 355670 times.
✓ Branch 1 taken 149 times.
 355819 for(index_t f = f_begin; f < f_end; ++f) {
713
2/2
✓ Branch 0 taken 1430227 times.
✓ Branch 1 taken 355670 times.
 1785897 for(index_t c = corners_begin(f); c < corners_end(f); ++c) {
714 geo_debug_assert(c >= c_begin);
715 1430227 c2f[c-c_begin] = f;
716 }
717 }
718 }
719
720
2/2
✓ Branch 0 taken 5358046 times.
✓ Branch 1 taken 718 times.
 5358764 for(index_t c = c_begin; c < c_end; ++c) {
721 5358046 facet_corners_.set_adjacent_facet(c, NO_FACET);
722 }
723
724 // Gives for each vertex a corner incident to it.
725
1/2
✓ Branch 1 taken 718 times.
✗ Branch 2 not taken.
 718 vector<index_t> v2c(v_end - v_begin, NO_CORNER);
726
727 // Chains the corners around each vertex.
728
1/2
✓ Branch 1 taken 718 times.
✗ Branch 2 not taken.
 718 vector<index_t> next_corner_around_vertex(c_end - c_begin, NO_CORNER);
729
730 // Step 1: chain corners around vertices and compute v2c
731
2/2
✓ Branch 0 taken 1664943 times.
✓ Branch 1 taken 718 times.
 1665661 for(index_t f = f_begin; f < f_end; ++f) {
732
2/2
✓ Branch 0 taken 5358046 times.
✓ Branch 1 taken 1664943 times.
 7022989 for(index_t c = corners_begin(f); c < corners_end(f); ++c) {
733 5358046 index_t v = facet_corners_.vertex(c);
734 5358046 next_corner_around_vertex[c - c_begin] = v2c[v - v_begin];
735 5358046 v2c[v - v_begin] = c;
736 }
737 }
738
739 // Step 2: connect
740
2/2
✓ Branch 0 taken 1664943 times.
✓ Branch 1 taken 718 times.
 1665661 for(index_t f1 = f_begin; f1 < f_end; ++f1) {
741
2/2
✓ Branch 0 taken 5358046 times.
✓ Branch 1 taken 1664943 times.
 12381035 for(index_t c1 = corners_begin(f1); c1 < corners_end(f1); ++c1) {
742
2/2
✓ Branch 0 taken 2693344 times.
✓ Branch 1 taken 2664702 times.
 5358046 if(facet_corners_.adjacent_facet(c1) == NO_FACET) {
743
744 index_t nb_candidates = 0;
745 index_t c_candidate = NO_CORNER;
746
747 index_t v1 = facet_corners_.vertex(c1);
748 index_t v2 = facet_corners_.vertex(
749 next_corner_around_facet(f1, c1)
750 );
751
752 // Traverse all the corners c2 incident to v1, and
753 // find among them the one(s) that is opposite to c1
754 2693344 for(
755 2693344 index_t c2 = v2c[v1 - v_begin];
756
2/2
✓ Branch 0 taken 33718738 times.
✓ Branch 1 taken 2693344 times.
 36412082 c2 != NO_CORNER;
757 33718738 c2 = next_corner_around_vertex[c2 - c_begin]
758 ) {
759
2/2
✓ Branch 0 taken 31025394 times.
✓ Branch 1 taken 2693344 times.
 33718738 if(c2 != c1) {
760 index_t f2 =
761
2/2
✓ Branch 0 taken 28790127 times.
✓ Branch 1 taken 2235267 times.
 31025394 is_simplicial_ ? c2/3 : c2f[c2 - c_begin];
762
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 31025394 times.
 31025394 index_t c2_prev = prev_corner_around_facet(f2, c2);
763
764 index_t v3 = facet_corners_.vertex(c2);
765 index_t v4 = facet_corners_.vertex(c2_prev);
766
767
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 31025394 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 31025394 geo_assert(v1 == v3);
768
769 if(
770
4/4
✓ Branch 0 taken 2665886 times.
✓ Branch 1 taken 28359508 times.
✓ Branch 2 taken 2665822 times.
✓ Branch 3 taken 64 times.
 31025394 v4 == v2 && (
771 facet_corners_.adjacent_facet(c2_prev) ==
772 NO_FACET
773 )
774 ) {
775 c_candidate = c2_prev;
776 2665822 ++nb_candidates;
777 }
778 }
779 }
780 // If there were more than 1 candidate, do not connect.
781
2/2
✓ Branch 0 taken 2664702 times.
✓ Branch 1 taken 28642 times.
 2693344 if(nb_candidates == 1) {
782 index_t c2 = c_candidate;
783
2/2
✓ Branch 0 taken 1955272 times.
✓ Branch 1 taken 709430 times.
 2664702 index_t f2 = is_simplicial_ ? (c2/3) : c2f[c2 - c_begin];
784 facet_corners_.set_adjacent_facet(c1,f2);
785 facet_corners_.set_adjacent_facet(c2,f1);
786 }
787 }
788 }
789 }
790 }
791
792 110 void MeshFacets::triangulate() {
793
2/2
✓ Branch 0 taken 91 times.
✓ Branch 1 taken 19 times.
 110 if(is_simplicial_) {
794 91 return;
795 }
796 index_t nb_triangles = 0;
797
2/2
✓ Branch 0 taken 13151 times.
✓ Branch 1 taken 19 times.
 13170 for(index_t f = 0; f < nb(); f++) {
798 13151 nb_triangles += (nb_vertices(f) - 2);
799 }
800 vector<index_t> new_corner_vertex_index;
801
1/2
✓ Branch 1 taken 19 times.
✗ Branch 2 not taken.
 19 new_corner_vertex_index.reserve(nb_triangles * 3);
802
2/2
✓ Branch 0 taken 13151 times.
✓ Branch 1 taken 19 times.
 13170 for(index_t f = 0; f < nb(); f++) {
803
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 13151 times.
 13151 index_t v0 = facet_corners_.vertex(corners_begin(f));
804 13151 for(index_t c = corners_begin(f) + 1;
805
3/4
✗ Branch 0 not taken.
✓ Branch 1 taken 38637 times.
✓ Branch 2 taken 25486 times.
✓ Branch 3 taken 13151 times.
 77274 c + 1 < corners_end(f); ++c
806 ) {
807 new_corner_vertex_index.push_back(v0);
808 new_corner_vertex_index.push_back(
809
1/2
✓ Branch 1 taken 25486 times.
✗ Branch 2 not taken.
 25486 facet_corners_.vertex(c)
810 );
811 new_corner_vertex_index.push_back(
812
1/4
✓ Branch 1 taken 25486 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 25486 facet_corners_.vertex(c + 1)
813 );
814 }
815 }
816
1/2
✓ Branch 1 taken 19 times.
✗ Branch 2 not taken.
 19 assign_triangle_mesh(new_corner_vertex_index, true);
817 }
818
819
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 718548 times.
 718548 void MeshFacets::flip(index_t f) {
820 index_t d = nb_vertices(f);
821
822 // Allocated on the stack (more multithread-friendly
823 // and no need to free)
824 index_t* corner_vertex_index =
825 718548 (index_t*) alloca(sizeof(index_t) * d);
826
827 index_t* corner_adjacent_facet =
828
1/2
✓ Branch 0 taken 718548 times.
✗ Branch 1 not taken.
 718548 (index_t*) alloca(sizeof(index_t) * d);
829
830 index_t c0 = corners_begin(f);
831
2/2
✓ Branch 0 taken 2155644 times.
✓ Branch 1 taken 718548 times.
 2874192 for(index_t i = 0; i < d; i++) {
832 2155644 corner_vertex_index[i] = facet_corners_.vertex(c0 + i);
833 2155644 corner_adjacent_facet[i] = facet_corners_.adjacent_facet(c0 + i);
834 }
835
2/2
✓ Branch 0 taken 2155644 times.
✓ Branch 1 taken 718548 times.
 2874192 for(index_t i = 0; i < d; i++) {
836 2155644 index_t i_v = d - 1 - i;
837
2/2
✓ Branch 0 taken 718548 times.
✓ Branch 1 taken 1437096 times.
 2155644 index_t i_f = (i_v == 0) ? d - 1 : i_v - 1;
838 2155644 facet_corners_.set_vertex(c0 + i, corner_vertex_index[i_v]);
839 2155644 facet_corners_.set_adjacent_facet(
840 2155644 c0 + i, corner_adjacent_facet[i_f]
841 );
842 }
843
2/2
✓ Branch 0 taken 718548 times.
✓ Branch 1 taken 718548 times.
 1437096 for(index_t i=0; i<d/2; i++) {
844 718548 mesh_.facet_corners.attributes().swap_items(c0+i,c0+d-1-i);
845 }
846 718548 }
847
848 144 void MeshFacets::compute_borders() {
849 144 mesh_.edges.clear();
850
2/2
✓ Branch 0 taken 4167 times.
✓ Branch 1 taken 144 times.
 4311 for(index_t f=0; f<nb(); ++f) {
851
2/2
✓ Branch 0 taken 12501 times.
✓ Branch 1 taken 4167 times.
 29169 for(index_t c1=corners_begin(f); c1!=corners_end(f); ++c1) {
852
2/2
✓ Branch 0 taken 4159 times.
✓ Branch 1 taken 8342 times.
 12501 if(mesh_.facet_corners.adjacent_facet(c1) == NO_FACET) {
853 index_t c2 = next_corner_around_facet(f,c1);
854 4159 mesh_.edges.create_edge(
855 mesh_.facet_corners.vertex(c1),
856 mesh_.facet_corners.vertex(c2)
857 );
858 }
859 }
860 }
861 144 }
862
863 11 void MeshFacets::assign_triangle_mesh(
864 coord_index_t dim,
865 vector<double>& vertices,
866 vector<index_t>& triangles,
867 bool steal_args
868 ) {
869 11 vertices_.assign_points(vertices, dim, steal_args);
870 11 assign_triangle_mesh(triangles, steal_args);
871 11 }
872
873
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 32 times.
 32 void MeshFacets::assign_triangle_mesh(
874 vector<index_t>& triangles,
875 bool steal_args
876 ) {
877
2/2
✓ Branch 0 taken 19 times.
✓ Branch 1 taken 13 times.
 32 index_t nb_triangles = triangles.size()/3;
878 is_simplicial();
879 facet_ptr_.clear();
880 32 resize_store(nb_triangles);
881
1/2
✓ Branch 0 taken 32 times.
✗ Branch 1 not taken.
 32 if(steal_args) {
882 32 facet_corners_.corner_vertex_.swap(triangles);
883 } else {
884 ✗ facet_corners_.corner_vertex_ = triangles;
885 }
886 32 facet_corners_.resize_store(nb_triangles*3);
887 32 facet_corners_.corner_adjacent_facet_.assign(
888 nb_triangles*3, NO_FACET
889 );
890 32 attributes().zero();
891 32 facet_corners_.attributes().zero();
892 32 }
893
894 185 void MeshFacets::pop() {
895 geo_debug_assert(nb() != 0);
896 index_t new_nb_corners =
897
1/2
✓ Branch 0 taken 185 times.
✗ Branch 1 not taken.
 185 is_simplicial_ ? 3*(nb()-1) : facet_ptr_[nb()-1];
898 185 resize_store(nb()-1);
899 185 facet_corners_.resize_store(new_nb_corners);
900 185 }
901
902 /**************************************************************************/
903
904 namespace MeshCellDescriptors {
905
906 GEOGRAM_API CellDescriptor tet_descriptor = {
907 4, // nb_vertices
908 4, // nb_facets
909 {3,3,3,3}, // nb_vertices in facet
910 { // facets
911 {1,3,2},
912 {0,2,3},
913 {3,1,0},
914 {0,1,2}
915 },
916 6, // nb_edges
917 { // edges
918 {1,2}, {2,3}, {3,1}, {0,1}, {0,2}, {0,3}
919 },
920 { // edges adjacent facets
921 {0,3}, {0,1}, {0,2}, {2,3}, {3,1}, {1,2}
922 }
923 };
924
925
926 GEOGRAM_API CellDescriptor hex_descriptor = {
927 8, // nb_vertices
928 6, // nb_facets
929 {4,4,4,4,4,4}, // nb_vertices in facet
930 { // facets
931 {0,2,6,4},
932 {3,1,5,7},
933 {1,0,4,5},
934 {2,3,7,6},
935 {1,3,2,0},
936 {4,6,7,5}
937 },
938 12, // nb_edges
939 { // edges
940 {0,1},{1,3},{3,2},{2,0},{4,5},{5,7},
941 {7,6},{6,4},{0,4},{1,5},{3,7},{2,6}
942 },
943 { // edges adjacent facets
944 {4,2},{4,1},{4,3},{4,0},{2,5},{1,5},
945 {3,5},{0,5},{2,0},{1,2},{3,1},{0,3}
946 }
947 };
948
949 GEOGRAM_API CellDescriptor prism_descriptor = {
950 6, // nb_vertices
951 5, // nb_facets
952 {3,3,4,4,4}, // nb_vertices in facet
953 { // facets
954 {0,1,2},
955 {3,5,4},
956 {0,3,4,1},
957 {0,2,5,3},
958 {1,4,5,2}
959 },
960 9, // nb_edges
961 { // edges
962 {0,1},{1,2},{2,0},{3,4},{4,5},{5,3},{0,3},{1,4},{2,5}
963 },
964 { // edges adjacent facets
965 {2,0},{4,0},{3,0},{1,2},{1,4},{1,3},{3,2},{2,4},{4,3}
966 }
967 };
968
969
970 GEOGRAM_API CellDescriptor pyramid_descriptor = {
971 5, // nb_vertices
972 5, // nb_facets
973 {4,3,3,3,3}, // nb_vertices in facet
974 { // facets
975 {0,1,2,3},
976 {0,4,1},
977 {0,3,4},
978 {2,4,3},
979 {2,1,4}
980 },
981 8, // nb_edges
982 { // edges
983 {0,1},{1,2},{2,3},{3,0},{0,4},{1,4},{2,4},{3,4}
984 },
985 { // edges adjacent facets
986 {1,0},{4,0},{3,0},{2,0},{2,1},{1,4},{4,3},{3,2}
987 }
988 };
989
990 GEOGRAM_API CellDescriptor connector_descriptor = {
991 4, // nb_vertices
992 3, // nb_facets
993 {4,3,3}, // nb_vertices in facet
994 { // facets
995 {0,1,2,3},
996 {2,1,0},
997 {3,2,0}
998 },
999 5, // nb_edges
1000 { // edges
1001 {0,1},{1,2},{2,3},{3,0},{0,2}
1002 },
1003 { // edges adjacent facets
1004 {1,0},{1,0},{2,0},{2,0},{2,1}
1005 }
1006 };
1007
1008 GEOGRAM_API CellDescriptor* cell_type_to_cell_descriptor[5] = {
1009 &tet_descriptor,
1010 &hex_descriptor,
1011 &prism_descriptor,
1012 &pyramid_descriptor,
1013 &connector_descriptor
1014 };
1015
1016 }
1017
1018 /********************************************************************/
1019
1020 819 MeshCellsStore::MeshCellsStore(Mesh& mesh) :
1021 MeshSubElementsStore(mesh),
1022
1/2
✓ Branch 2 taken 819 times.
✗ Branch 3 not taken.
 819 is_simplicial_(true) {
1023
1/4
✓ Branch 1 taken 819 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 819 cell_ptr_.push_back(0);
1024 819 }
1025
1026 357 void MeshCellsStore::clear_store(
1027 bool keep_attributes, bool keep_memory
1028 ) {
1029
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 357 times.
 357 if(keep_memory) {
1030 ✗ cell_ptr_.resize(0);
1031 ✗ cell_type_.resize(0);
1032 } else {
1033 cell_ptr_.clear();
1034 cell_type_.clear();
1035 }
1036 357 cell_ptr_.push_back(0);
1037 357 MeshSubElementsStore::clear_store(keep_attributes, keep_memory);
1038 357 }
1039
1040 4 void MeshCellsStore::resize_store(index_t new_size) {
1041
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 4 if(!is_simplicial_) {
1042 ✗ cell_ptr_.resize(new_size+1);
1043 ✗ cell_type_.resize(new_size);
1044 }
1045 4 MeshSubElementsStore::resize_store(new_size);
1046 4 }
1047
1048
1049 10368 const CellDescriptor& MeshCellsStore::descriptor(index_t c) const {
1050 geo_debug_assert(c < nb());
1051 10368 return is_simplicial_ ? MeshCellDescriptors::tet_descriptor :
1052 *(
1053 MeshCellDescriptors::cell_type_to_cell_descriptor[
1054 cell_type_[c]
1055 ✗ ]
1056
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10368 times.
 10368 );
1057 }
1058
1059 4 const CellDescriptor& MeshCellsStore::cell_type_to_cell_descriptor(
1060 MeshCellType t
1061 ) {
1062 geo_debug_assert(t < GEO::MESH_NB_CELL_TYPES);
1063 4 return *(MeshCellDescriptors::cell_type_to_cell_descriptor[t]);
1064 }
1065
1066 /**************************************************************************/
1067
1068 819 MeshCellCornersStore::MeshCellCornersStore(Mesh& mesh) :
1069 MeshSubElementsStore(mesh),
1070 819 vertices_(mesh.vertices) {
1071 819 }
1072
1073 357 void MeshCellCornersStore::clear_store(
1074 bool keep_attributes, bool keep_memory
1075 ) {
1076
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 357 times.
 357 if(keep_memory) {
1077 ✗ corner_vertex_.resize(0);
1078 } else {
1079 corner_vertex_.clear();
1080 }
1081 357 MeshSubElementsStore::clear_store(keep_attributes, keep_memory);
1082 357 }
1083
1084 4 void MeshCellCornersStore::resize_store(index_t new_size) {
1085 4 corner_vertex_.resize(new_size);
1086 4 MeshSubElementsStore::resize_store(new_size);
1087 4 }
1088
1089 /**************************************************************************/
1090
1091 819 MeshCellFacetsStore::MeshCellFacetsStore(Mesh& mesh) :
1092 MeshSubElementsStore(mesh),
1093 819 vertices_(mesh.vertices),
1094 819 cells_(mesh.cells) {
1095 819 }
1096
1097 357 void MeshCellFacetsStore::clear_store(
1098 bool keep_attributes, bool keep_memory
1099 ) {
1100
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 357 times.
 357 if(keep_memory) {
1101 ✗ adjacent_cell_.resize(0);
1102 } else {
1103 adjacent_cell_.clear();
1104 }
1105 357 MeshSubElementsStore::clear_store(keep_attributes, keep_memory);
1106 357 }
1107
1108 8 void MeshCellFacetsStore::resize_store(index_t new_size) {
1109 8 adjacent_cell_.resize(new_size);
1110 8 MeshSubElementsStore::resize_store(new_size);
1111 8 }
1112
1113
1114 /**************************************************************************/
1115
1116 819 MeshCells::MeshCells(Mesh& mesh) :
1117 MeshCellsStore(mesh),
1118 819 vertices_(mesh.vertices),
1119 819 cell_corners_(mesh.cell_corners),
1120
1/2
✓ Branch 2 taken 819 times.
✗ Branch 3 not taken.
 819 cell_facets_(mesh.cell_facets) {
1121 819 }
1122
1123 357 void MeshCells::clear(bool keep_attributes, bool keep_memory) {
1124 357 cell_corners_.clear_store(keep_attributes, keep_memory);
1125 357 cell_facets_.clear_store(keep_attributes, keep_memory);
1126 357 clear_store(keep_attributes, keep_memory);
1127 357 is_simplicial_ = true;
1128 357 }
1129
1130
1131 ✗ void MeshCells::delete_elements(
1132 vector<index_t>& to_delete,
1133 bool remove_isolated_vertices
1134 ) {
1135 // "Fast track" if no element should be deleted
1136 ✗ if(!has_non_zero(to_delete)) {
1137 ✗ if(remove_isolated_vertices) {
1138 ✗ mesh_.vertices.remove_isolated();
1139 }
1140 ✗ return;
1141 }
1142
1143 // to_delete is used for both indicating
1144 // which facets should be deleted and
1145 // for storing the re-numbering map
1146 vector<index_t>& cells_old2new = to_delete;
1147
1148 ✗ vector<index_t>& corner_vertex = cell_corners_.corner_vertex_;
1149 ✗ vector<index_t>& adjacent_cell = cell_facets_.adjacent_cell_;
1150
1151 index_t new_nb_cells = 0;
1152 index_t new_nb_corner_facets = 0;
1153
1154 // If there are some corners or facets
1155 // attributes, we need to compute the index
1156 // mapping for them.
1157 vector<index_t> corner_facets_old2new;
1158 if(
1159 ✗ cell_corners_.attributes().nb() != 0 ||
1160 cell_facets_.attributes().nb() != 0) {
1161 ✗ corner_facets_old2new.resize(cell_corners_.nb(), NO_INDEX);
1162 }
1163
1164 ✗ for(index_t c=0; c<nb(); ++c) {
1165 ✗ if(cells_old2new[c] != 0) {
1166 ✗ cells_old2new[c] = NO_CELL;
1167 } else {
1168 ✗ cells_old2new[c] = new_nb_cells;
1169
1170 ✗ if(!is_simplicial_) {
1171 ✗ cell_ptr_[new_nb_cells] = new_nb_corner_facets;
1172 ✗ cell_type_[new_nb_cells] = cell_type_[c];
1173 }
1174
1175 index_t b,e;
1176 ✗ if(is_simplicial_) {
1177 ✗ b = 4*c;
1178 ✗ e = b+4;
1179 } else {
1180 ✗ b = cell_ptr_[c];
1181 ✗ e = cell_ptr_[c+1];
1182 }
1183
1184 ✗ for(index_t cof=b; cof<e; ++cof) {
1185 ✗ if(corner_facets_old2new.size() != 0) {
1186 ✗ corner_facets_old2new[cof] = new_nb_corner_facets;
1187 }
1188 ✗ if(cof != new_nb_corner_facets) {
1189 ✗ corner_vertex[new_nb_corner_facets] =
1190 corner_vertex[cof];
1191 ✗ adjacent_cell[new_nb_corner_facets] =
1192 adjacent_cell[cof];
1193 }
1194 ✗ ++new_nb_corner_facets;
1195 }
1196 ✗ ++new_nb_cells;
1197 }
1198 }
1199
1200 ✗ if(!is_simplicial_) {
1201 ✗ cell_ptr_[new_nb_cells] = new_nb_corner_facets;
1202 }
1203
1204 // Map adjacent cell indices
1205 ✗ for(index_t f=0; f<cell_facets_.nb(); ++f) {
1206 ✗ index_t c = adjacent_cell[f];
1207 ✗ if(c != NO_CELL) {
1208 ✗ adjacent_cell[f] = cells_old2new[c];
1209 }
1210 }
1211
1212 // Manage cell store and attributes
1213 ✗ attributes().compress(cells_old2new);
1214 ✗ resize_store(new_nb_cells);
1215
1216 // Manage corners/facets store and attributes
1217 ✗ if(corner_facets_old2new.size() != 0) {
1218 // Corner/facet index mapping is only computed
1219 // if there ware some corner or facet attributes
1220 ✗ cell_corners_.attributes().compress(corner_facets_old2new);
1221 ✗ cell_facets_.attributes().compress(corner_facets_old2new);
1222 }
1223 ✗ cell_corners_.resize_store(new_nb_corner_facets);
1224 ✗ cell_facets_.resize_store(new_nb_corner_facets);
1225
1226 ✗ if(remove_isolated_vertices) {
1227 ✗ mesh_.vertices.remove_isolated();
1228 }
1229 }
1230
1231 4 void MeshCells::permute_elements(vector<index_t>& permutation) {
1232 4 attributes_.apply_permutation(permutation);
1233
1234 if(
1235
2/4
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 4 times.
 4 cell_corners_.attributes().nb() != 0 ||
1236
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 4 cell_facets_.attributes().nb() != 0
1237 ) {
1238 vector<index_t> cell_corner_facets_permutation;
1239 ✗ cell_corner_facets_permutation.reserve(cell_corners_.nb());
1240
1241 ✗ for(index_t new_cell = 0; new_cell<nb(); ++new_cell) {
1242 ✗ index_t old_cell = permutation[new_cell];
1243 index_t cell_size =
1244 ✗ std::max(nb_vertices(old_cell), nb_facets(old_cell));
1245 ✗ for(index_t i=0; i<cell_size; ++i) {
1246 cell_corner_facets_permutation.push_back(
1247 ✗ corners_begin(old_cell)+i
1248 );
1249 }
1250 }
1251
1252 ✗ if(cell_corners_.attributes().nb() != 0) {
1253 ✗ cell_corners_.attributes().apply_permutation(
1254 cell_corner_facets_permutation
1255 );
1256 }
1257 ✗ if(cell_facets_.attributes().nb() != 0) {
1258 ✗ cell_facets_.attributes().apply_permutation(
1259 cell_corner_facets_permutation
1260 );
1261 }
1262 }
1263
1264 4 vector<index_t>& corner_vertex = cell_corners_.corner_vertex_;
1265 4 vector<index_t>& facet_adjacent_cell = cell_facets_.adjacent_cell_;
1266
1267
1/2
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
 4 if(is_simplicial_) {
1268 // in-place permutation !
1269
1270 4 Permutation::apply(
1271 corner_vertex.data(),
1272 permutation,
1273 index_t(sizeof(index_t) * 4)
1274 );
1275
1276 4 Permutation::apply(
1277 facet_adjacent_cell.data(),
1278 permutation,
1279 index_t(sizeof(index_t) * 4)
1280 );
1281
1282 4 Permutation::invert(permutation);
1283
1284
2/2
✓ Branch 0 taken 6144 times.
✓ Branch 1 taken 4 times.
 12296 for(index_t f = 0; f < facet_adjacent_cell.size(); ++f) {
1285
2/2
✓ Branch 0 taken 5376 times.
✓ Branch 1 taken 768 times.
 6144 if(facet_adjacent_cell[f] != NO_CELL) {
1286 5376 facet_adjacent_cell[f] =
1287 permutation[facet_adjacent_cell[f]];
1288 }
1289 }
1290 } else {
1291 // we need to do some copies
1292
1293 ✗ vector<index_t> new_cell_ptr(nb()+1);
1294 ✗ vector<index_t> new_corner_vertex(cell_corners_.nb());
1295 ✗ vector<index_t> new_facet_adjacent_cell(cell_facets_.nb());
1296
1297 index_t new_ptr = 0;
1298 ✗ for(index_t new_c=0; new_c<nb(); ++new_c) {
1299 ✗ index_t old_c = permutation[new_c];
1300 ✗ index_t old_ptr = cell_ptr_[old_c];
1301 index_t cell_size = std::max(
1302 ✗ nb_vertices(old_c), nb_facets(old_c)
1303 );
1304 ✗ new_cell_ptr[new_c] = new_ptr;
1305 ✗ for(index_t i=0; i<cell_size; ++i) {
1306 ✗ new_corner_vertex[new_ptr+i] = corner_vertex[old_ptr+i];
1307 ✗ new_facet_adjacent_cell[new_ptr+i] =
1308 facet_adjacent_cell[old_ptr+i];
1309 }
1310 ✗ new_ptr += cell_size;
1311 }
1312 ✗ new_cell_ptr[nb()] = new_ptr;
1313
1314 ✗ Permutation::apply(
1315 cell_type_.data(), permutation, index_t(sizeof(Numeric::uint8))
1316 );
1317
1318 ✗ Permutation::invert(permutation);
1319
1320 ✗ for(index_t f = 0; f < new_facet_adjacent_cell.size(); ++f) {
1321 ✗ if(new_facet_adjacent_cell[f] != NO_CELL) {
1322 ✗ new_facet_adjacent_cell[f] =
1323 permutation[new_facet_adjacent_cell[f]];
1324 }
1325 }
1326
1327 corner_vertex.swap(new_corner_vertex);
1328 facet_adjacent_cell.swap(new_facet_adjacent_cell);
1329 cell_ptr_.swap(new_cell_ptr);
1330 }
1331 4 }
1332
1333
1334
1335 246 void MeshCells::connect_tets() {
1336
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 246 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 246 geo_assert(is_simplicial_);
1337
2/2
✓ Branch 0 taken 242 times.
✓ Branch 1 taken 4 times.
 246 if(nb() == 0) {
1338 242 return;
1339 }
1340 4 cell_facets_.resize_store(nb() * 4);
1341
2/2
✓ Branch 0 taken 6144 times.
✓ Branch 1 taken 4 times.
 6148 for(index_t f=0; f<cell_facets_.nb(); ++f) {
1342 cell_facets_.set_adjacent_cell(f,NO_CELL);
1343 }
1344
1345 GEO::vector<index_t> next_tet_corner_around_vertex(
1346 nb() * 4, NO_CORNER
1347 4 );
1348
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 GEO::vector<index_t> v2c(vertices_.nb(), NO_CORNER);
1349
1350 // Step 1: chain tet corners around vertices and compute v2c
1351
2/2
✓ Branch 0 taken 1536 times.
✓ Branch 1 taken 4 times.
 1540 for(index_t t = 0; t < nb(); ++t) {
1352
2/2
✓ Branch 0 taken 6144 times.
✓ Branch 1 taken 1536 times.
 7680 for(index_t lv = 0; lv < 4; ++lv) {
1353 index_t v = vertex(t, lv);
1354 6144 next_tet_corner_around_vertex[4 * t + lv] = v2c[v];
1355 6144 v2c[v] = 4 * t + lv;
1356 }
1357 }
1358
1359 // Step 2: connect tets
1360
2/2
✓ Branch 0 taken 1536 times.
✓ Branch 1 taken 4 times.
 1540 for(index_t t1 = 0; t1 < nb(); ++t1) {
1361
2/2
✓ Branch 0 taken 6144 times.
✓ Branch 1 taken 1536 times.
 7680 for(index_t lf1 = 0; lf1 < 4; ++lf1) {
1362
2/2
✓ Branch 0 taken 3456 times.
✓ Branch 1 taken 2688 times.
 6144 if(adjacent(t1, lf1) == NO_CELL) {
1363
1/2
✓ Branch 1 taken 3456 times.
✗ Branch 2 not taken.
 3456 index_t v1 = facet_vertex(t1, lf1, 0);
1364
1/2
✓ Branch 1 taken 3456 times.
✗ Branch 2 not taken.
 3456 index_t v2 = facet_vertex(t1, lf1, 1);
1365
1/2
✓ Branch 1 taken 3456 times.
✗ Branch 2 not taken.
 3456 index_t v3 = facet_vertex(t1, lf1, 2);
1366 3456 for(
1367
2/2
✓ Branch 0 taken 24100 times.
✓ Branch 1 taken 768 times.
 24868 index_t c2 = v2c[v1]; c2 != NO_CORNER;
1368 21412 c2 = next_tet_corner_around_vertex[c2]
1369 ) {
1370 24100 index_t t2 = c2/4;
1371 24100 index_t lf2 = find_tet_facet(t2, v3, v2, v1);
1372
2/2
✓ Branch 0 taken 2688 times.
✓ Branch 1 taken 21412 times.
 24100 if(lf2 != NO_FACET) {
1373 set_adjacent(t1, lf1, t2);
1374 set_adjacent(t2, lf2, t1);
1375 break;
1376 }
1377 }
1378 }
1379 }
1380 }
1381 }
1382
1383 ✗ bool MeshCells::facets_match(
1384 index_t c1, index_t f1, index_t c2, index_t f2
1385 ) const {
1386 index_t nbv = facet_nb_vertices(c1,f1);
1387 ✗ if(facet_nb_vertices(c2,f2) != nbv) {
1388 return false;
1389 }
1390 ✗ for(index_t offset=0; offset<nbv; ++offset) {
1391 bool match=true;
1392 ✗ for(index_t v1=0; v1<nbv; ++v1) {
1393 ✗ index_t v2 = (nbv-v1+offset)%nbv;
1394 ✗ if(
1395 ✗ facet_vertex(c1,f1,v1) !=
1396 ✗ facet_vertex(c2,f2,v2)
1397 ) {
1398 match=false;
1399 break;
1400 }
1401 }
1402 ✗ if(match) {
1403 return true;
1404 }
1405 }
1406 return false;
1407 }
1408
1409 /**
1410 * \brief Tests whether two indices triplets match
1411 * up to a circular permutation.
1412 * \param[in] v1 index of the first vertex of the first triangle
1413 * \param[in] v2 index of the second vertex of the first triangle
1414 * \param[in] v3 index of the third vertex of the first triangle
1415 * \param[in] w1 index of the first vertex of the second triangle
1416 * \param[in] w2 index of the second vertex of the second triangle
1417 * \param[in] w3 index of the third vertex of the second triangle
1418 * \retval true if (\p v1, \p v2, \p v3) = (\p w1, \p w2, \p w3)
1419 * up to a circular permutation
1420 * \retval false otherwise
1421 */
1422 inline bool triangles_equal(
1423 index_t v1, index_t v2, index_t v3,
1424 index_t w1, index_t w2, index_t w3
1425 ) {
1426 return (
1427 ✗ (v1 == w1 && v2 == w2 && v3 == w3) ||
1428 ✗ (v1 == w2 && v2 == w3 && v3 == w1) ||
1429 ✗ (v1 == w3 && v2 == w1 && v3 == w2)
1430 );
1431 }
1432
1433 ✗ bool MeshCells::triangular_facet_matches_quad_facet(
1434 index_t c1, index_t lf1,
1435 index_t c2, index_t lf2
1436 ) const {
1437 geo_debug_assert(facet_nb_vertices(c1,lf1) == 3);
1438 geo_debug_assert(facet_nb_vertices(c2,lf2) == 4);
1439
1440 ✗ index_t v1 = facet_vertex(c1,lf1,0);
1441 ✗ index_t v2 = facet_vertex(c1,lf1,1);
1442 ✗ index_t v3 = facet_vertex(c1,lf1,2);
1443 ✗ index_t w1 = facet_vertex(c2,lf2,0);
1444 ✗ index_t w2 = facet_vertex(c2,lf2,1);
1445 ✗ index_t w3 = facet_vertex(c2,lf2,2);
1446 ✗ index_t w4 = facet_vertex(c2,lf2,3);
1447
1448 // Note: subtriangles in (w1,w2,w3,w4) are
1449 // in reverse order since two facets can be
1450 // connected only if they have opposite
1451 // orientations.
1452 return (
1453 triangles_equal(v1,v2,v3,w4,w3,w2) ||
1454 triangles_equal(v1,v2,v3,w3,w2,w1) ||
1455 triangles_equal(v1,v2,v3,w2,w1,w4) ||
1456 triangles_equal(v1,v2,v3,w1,w4,w3)
1457 ✗ ) ;
1458 }
1459
1460 ✗ bool MeshCells::triangular_facets_have_common_edge(
1461 index_t c1, index_t f1,
1462 index_t c2, index_t f2,
1463 index_t& e1, index_t& e2
1464 ) const {
1465 geo_debug_assert(facet_nb_vertices(c1,f1) == 3);
1466 geo_debug_assert(facet_nb_vertices(c2,f2) == 3);
1467 ✗ for(e1=0; e1<3; ++e1) {
1468 ✗ for(e2=0; e2<3; ++e2) {
1469 if(
1470 ✗ facet_vertex(c1, f1, (e1+1)%3) ==
1471 ✗ facet_vertex(c2, f2, (e2+2)%3) &&
1472 ✗ facet_vertex(c1, f1, (e1+2)%3) ==
1473 ✗ facet_vertex(c2, f2, (e2+1)%3)
1474 ) {
1475 return true;
1476 }
1477 }
1478 }
1479 ✗ e1 = NO_EDGE;
1480 ✗ e2 = NO_EDGE;
1481 ✗ return false;
1482 }
1483
1484 ✗ bool MeshCells::create_connector(
1485 index_t c1, index_t lf1,
1486 const std::vector< std::pair<index_t, index_t> >& matches
1487 ) {
1488 ✗ if(matches.size() == 0) {
1489 return false;
1490 }
1491
1492 ✗ if(matches.size() == 1) {
1493 ✗ GEO::Logger::warn("Mesh")
1494 << "Found only one triangular facet adjacent to a quad facet"
1495 << std::endl;
1496 ✗ Attribute<bool> weird(attributes(),"weird");
1497 weird[c1] = true;
1498 ✗ for(index_t i=0; i<matches.size(); ++i) {
1499 ✗ weird[matches[i].first] = true;
1500 }
1501 return false;
1502 }
1503
1504 // Find among the matches two facets that have an edge in common
1505 // Yes, there can be more than two candidate facets with three
1506 // vertices in common with the quad facet ! But among them,
1507 // only two of them have an edge in common.
1508 index_t adj_c1 = NO_CELL;
1509 index_t adj_lf1 = NO_FACET;
1510 index_t adj_c2 = NO_CELL;
1511 index_t adj_lf2 = NO_FACET;
1512 index_t e1 = NO_EDGE;
1513 index_t e2 = NO_EDGE;
1514
1515 index_t nb_found=0;
1516 ✗ for(index_t i=0; i<index_t(matches.size()); ++i) {
1517 ✗ for(index_t j=i+1; j<index_t(matches.size()); ++j) {
1518 ✗ index_t cur_e1 = NO_EDGE;
1519 ✗ index_t cur_e2 = NO_EDGE;
1520 ✗ if(triangular_facets_have_common_edge(
1521 ✗ matches[i].first, matches[i].second,
1522 ✗ matches[j].first, matches[j].second,
1523 cur_e1, cur_e2
1524 )) {
1525 ✗ adj_c1 = matches[i].first;
1526 ✗ adj_lf1 = matches[i].second;
1527 ✗ adj_c2 = matches[j].first;
1528 ✗ adj_lf2 = matches[j].second;
1529 ✗ e1 = cur_e1;
1530 ✗ e2 = cur_e2;
1531 ✗ ++nb_found;
1532 }
1533 }
1534 }
1535
1536 // Sanity check: make sure that we only found a single pair
1537 // of triangular facets with a common edge that matches the quad.
1538 ✗ if(nb_found > 2) {
1539 ✗ GEO::Logger::warn("Mesh")
1540 << "Found more than two triangular facets adjacent to a quad"
1541 << " ( got " << nb_found << ")"
1542 << std::endl;
1543 ✗ Attribute<bool> weird(attributes(),"weird");
1544 weird[c1] = true;
1545 ✗ for(index_t i=0; i<matches.size(); ++i) {
1546 ✗ weird[matches[i].first] = true;
1547 }
1548
1549 return false;
1550 }
1551
1552 ✗ if(nb_found == 0) {
1553 ✗ GEO::Logger::warn("Mesh")
1554 << "Triangular facets adjacent to a quad have no common edge"
1555 << std::endl;
1556 ✗ return false;
1557 }
1558
1559 // Sanity check: make sure the triangular facets
1560 // are on the border.
1561 if(
1562 ✗ adjacent(adj_c1, adj_lf1) != NO_CELL ||
1563 adjacent(adj_c2, adj_lf2) != NO_CELL
1564 ) {
1565 /*
1566 GEO::Logger::warn("Mesh")
1567 << "Matching tet facets are not on border (\"thick sliver\")"
1568 << std::endl;
1569 */
1570 return false;
1571 }
1572
1573 // v1 and v2 are on the common edge
1574 ✗ index_t v1 = facet_vertex(
1575 ✗ adj_c1, adj_lf1, (e1+1)%3
1576 );
1577
1578 ✗ index_t v2 = facet_vertex(
1579 ✗ adj_c1, adj_lf1, (e1+2)%3
1580 );
1581
1582 // w1 and w2 are the opposite vertices
1583 ✗ index_t w1 = facet_vertex(adj_c1, adj_lf1, e1);
1584 ✗ index_t w2 = facet_vertex(adj_c2, adj_lf2, e2);
1585
1586 // Create the connector
1587 ✗ index_t conn = create_connector(
1588 v1, w2, v2, w1,
1589 c1, adj_c1, adj_c2
1590 );
1591
1592 // Connect the cells with the connector
1593 set_adjacent(c1, lf1, conn);
1594 set_adjacent(adj_c1, adj_lf1, conn);
1595 set_adjacent(adj_c2, adj_lf2, conn);
1596
1597 ✗ return true;
1598 }
1599
1600 246 void MeshCells::connect(bool remove_trivial_slivers, bool verbose_if_OK) {
1601 // "Fast track" for simplicial mesh
1602
1/2
✓ Branch 0 taken 246 times.
✗ Branch 1 not taken.
 246 if(is_simplicial_) {
1603 246 connect_tets();
1604 246 return;
1605 }
1606
1607 ✗ for(index_t f=0; f<cell_facets_.nb(); ++f) {
1608 cell_facets_.set_adjacent_cell(f,NO_CELL);
1609 }
1610
1611 vector<index_t> next_cell_around_vertex(
1612 ✗ cell_corners_.nb(), NO_CELL
1613 );
1614 ✗ vector<index_t> v2cell(vertices_.nb(), NO_CELL);
1615
1616 // Step 1: chain cells around vertices and compute v2cell
1617 ✗ for(index_t c = 0; c < nb(); ++c) {
1618 ✗ for(index_t lv = 0; lv < nb_vertices(c); ++lv) {
1619 index_t v = vertex(c, lv);
1620 ✗ next_cell_around_vertex[corners_begin(c) + lv] =
1621 v2cell[v];
1622 ✗ v2cell[v] = c;
1623 }
1624 }
1625
1626 // Step 2: connect cells
1627 // (c1,lf1) traverse all the cell facets
1628 ✗ for(index_t c1 = 0; c1 < nb(); ++c1) {
1629 ✗ for(index_t lf1 = 0; lf1 < nb_facets(c1); ++lf1) {
1630
1631 // If (c1,lf1) is on the border, try to connect it
1632 ✗ if(adjacent(c1, lf1) == NO_CELL) {
1633
1634 // v1 is one of the vertices of (c1,lf1)
1635 ✗ index_t v1 = facet_vertex(c1,lf1,0);
1636
1637 // c2 traverses all the cells incident to v1
1638 ✗ for(
1639 ✗ index_t c2 = v2cell[v1]; c2 != NO_CELL;
1640 ✗ c2 = next_cell_around_vertex[
1641 corners_begin(c2) +
1642 ✗ find_cell_vertex(c2,v1)
1643 ✗ ]
1644 ) {
1645
1646 // If we find a cell facet lf2 compatible with (c1,lf1)
1647 // in c2, then connect (c1,lf1) to c2 and
1648 // (c2,lf2) to c1.
1649 ✗ index_t lf2 = find_cell_facet(c2, c1, lf1);
1650 ✗ if(lf2 != NO_FACET) {
1651 set_adjacent(c1, lf1, c2);
1652 set_adjacent(c2, lf2, c1);
1653 break;
1654 }
1655 }
1656 }
1657 }
1658 }
1659
1660 // Step 3: Create connectors, i.e. artificial cells that represent
1661 // non-conformal connections between two triangular facets and
1662 // a quadrangular facet.
1663
1664 // Backup nb_cells since we are creating new cells (connectors)
1665 // during this loop.
1666 index_t nb_cells0 = nb();
1667
1668 // Keep track of the number of invalid configurations (does
1669 // not seem to happen anymore, but I keep the code just in case).
1670 index_t weird=0;
1671
1672 // For each quadrangular face, we compute the list of candidate
1673 // triangular faces to be connected with it (a vector of
1674 // (cell index, facet index) pairs).
1675 std::vector< std::pair<index_t, index_t> > matches;
1676
1677 // If remove_trivial_slivers is set, we also detect the trivial
1678 // slivers, i.e. the slivers that are glued on a quadrilateral facet.
1679
1680 std::vector<index_t> trivial_slivers;
1681
1682 // (c1,f1) traverse all quadrangular cell facets on the border
1683 ✗ for(index_t c1=0; c1 < nb_cells0; ++c1) {
1684 ✗ if(type(c1) == MESH_TET) {
1685 ✗ continue;
1686 }
1687 ✗ for(index_t lf1=0; lf1<nb_facets(c1); ++lf1) {
1688 ✗ if(
1689 ✗ facet_nb_vertices(c1,lf1) != 4 ||
1690 adjacent(c1,lf1) != NO_CELL
1691 ) {
1692 ✗ continue;
1693 }
1694
1695 // Now c2 traverses all the cells incident to one of
1696 // the vertex of (c1,lf1)
1697 ✗ matches.resize(0);
1698 ✗ for(index_t lv1=0; lv1<facet_nb_vertices(c1,lf1); ++lv1) {
1699 ✗ index_t v1 = facet_vertex(c1,lf1,lv1);
1700 ✗ for(
1701 ✗ index_t c2 = v2cell[v1]; c2 != NO_CELL;
1702 c2 = next_cell_around_vertex[
1703 corners_begin(c2) +
1704 ✗ find_cell_vertex(c2,v1)
1705 ✗ ]
1706 ) {
1707 geo_debug_assert(find_cell_vertex(c2,v1) != NO_VERTEX);
1708 ✗ if(c2 == c1 || type(c2) == MESH_HEX) {
1709 ✗ continue;
1710 }
1711
1712 // Among all the triangular facets of c2, find the ones
1713 // that can be connected to (c1,lf1)
1714 ✗ for(index_t lf2=0; lf2<nb_facets(c2); ++lf2) {
1715 ✗ if(facet_nb_vertices(c2,lf2) != 3) {
1716 ✗ continue;
1717 }
1718 ✗ if(triangular_facet_matches_quad_facet(
1719 c2,lf2,c1,lf1
1720 )) {
1721 ✗ matches.push_back(std::make_pair(c2,lf2));
1722 }
1723 }
1724 }
1725 }
1726
1727 // Make sure we get each match once only
1728 ✗ GEO::sort_unique(matches);
1729
1730 // This should not happen, but we keep this
1731 // sanity check and notify the user if some
1732 // connectors could not be created.
1733 if(
1734 ✗ matches.size() != 0 &&
1735 ✗ !create_connector(c1,lf1,matches)
1736 ) {
1737 ✗ ++weird;
1738 }
1739
1740 ✗ if(remove_trivial_slivers) {
1741 ✗ for(index_t i=0; i<matches.size(); ++i) {
1742 ✗ if(type(matches[i].first) != MESH_TET) {
1743 ✗ continue;
1744 }
1745 ✗ for(index_t j=i+1; j<matches.size(); ++j) {
1746 ✗ if(matches[j].first == matches[i].first) {
1747 ✗ trivial_slivers.push_back(matches[i].first);
1748 }
1749 }
1750 }
1751 }
1752 }
1753 }
1754 ✗ if(weird != 0) {
1755 ✗ GEO::Logger::warn("Mesh") << "Encountered "
1756 << weird
1757 << " invalid connector configurations"
1758 << std::endl;
1759 } else {
1760 ✗ if(verbose_if_OK) {
1761 ✗ GEO::Logger::out("Mesh") << "All connectors are OK"
1762 << std::endl;
1763 }
1764 }
1765 ✗ if(remove_trivial_slivers && trivial_slivers.size() != 0) {
1766 ✗ GEO::Logger::warn("Mesh") << "Removing "
1767 << trivial_slivers.size()
1768 << " trivial sliver(s)" << std::endl;
1769
1770 next_cell_around_vertex.clear();
1771 v2cell.clear();
1772
1773 ✗ vector<index_t> delete_c(nb(),0);
1774 ✗ for(index_t i=0; i<trivial_slivers.size(); ++i) {
1775 ✗ delete_c[trivial_slivers[i]] = 1;
1776 }
1777 // We need to remove the previously generated connectors,
1778 // some of them may be wrong if adjacent to a sliver that
1779 // was removed.
1780 ✗ for(index_t c=0; c<nb(); ++c) {
1781 ✗ if(type(c) == MESH_CONNECTOR) {
1782 ✗ delete_c[c] = 1;
1783 }
1784 }
1785 ✗ delete_elements(delete_c);
1786
1787 ✗ GEO::Logger::warn("Mesh")
1788 << "Re-trying to connect cells" << std::endl;
1789 ✗ connect(false,true);
1790 }
1791 }
1792
1793
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 void MeshCells::compute_borders() {
1794 Attribute<index_t> facet_cell;
1795
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 compute_borders(facet_cell);
1796 4 }
1797
1798 4 void MeshCells::compute_borders(Attribute<index_t>& facet_cell) {
1799 4 mesh_.facets.clear(true,false);
1800
1/2
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
 4 if(is_simplicial_) {
1801
2/2
✓ Branch 0 taken 1536 times.
✓ Branch 1 taken 4 times.
 1540 for(index_t t=0; t<nb(); ++t) {
1802
2/2
✓ Branch 0 taken 6144 times.
✓ Branch 1 taken 1536 times.
 7680 for(index_t f=0; f<4; ++f) {
1803
2/2
✓ Branch 0 taken 768 times.
✓ Branch 1 taken 5376 times.
 6144 if(adjacent(t,f) == NO_CELL) {
1804 768 index_t new_f = mesh_.facets.create_triangle(
1805 tet_facet_vertex(t,f,0),
1806 tet_facet_vertex(t,f,1),
1807 tet_facet_vertex(t,f,2)
1808 );
1809 if(facet_cell.is_bound()) {
1810 ✗ facet_cell[new_f] = t;
1811 }
1812 }
1813 }
1814 }
1815 } else {
1816 ✗ for(index_t c=0; c<nb(); ++c) {
1817 ✗ for(index_t f=0; f<nb_facets(c); ++f) {
1818 ✗ if(adjacent(c,f) == NO_CELL) {
1819 index_t new_f = NO_INDEX;
1820 ✗ switch(facet_nb_vertices(c,f)) {
1821 ✗ case 3:
1822 ✗ new_f = mesh_.facets.create_triangle(
1823 facet_vertex(c,f,0),
1824 facet_vertex(c,f,1),
1825 facet_vertex(c,f,2)
1826 );
1827 ✗ break;
1828 ✗ case 4:
1829 ✗ new_f = mesh_.facets.create_quad(
1830 facet_vertex(c,f,0),
1831 facet_vertex(c,f,1),
1832 facet_vertex(c,f,2),
1833 facet_vertex(c,f,3)
1834 );
1835 ✗ break;
1836 default:
1837 ✗ geo_assert_not_reached;
1838 }
1839 if(facet_cell.is_bound()) {
1840 ✗ facet_cell[new_f] = c;
1841 }
1842 }
1843 }
1844 }
1845 }
1846 4 mesh_.facets.connect();
1847 4 }
1848
1849 4 void MeshCells::assign_tet_mesh(
1850 coord_index_t dim,
1851 vector<double>& vertices,
1852 vector<index_t>& tets,
1853 bool steal_args
1854 ) {
1855 4 vertices_.assign_points(vertices, dim, steal_args);
1856 4 assign_tet_mesh(tets, steal_args);
1857 4 }
1858
1859
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 4 void MeshCells::assign_tet_mesh(
1860 vector<index_t>& tets,
1861 bool steal_args
1862 ) {
1863 4 index_t nb_tets = tets.size()/4;
1864
1/2
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
 4 is_simplicial_ = true;
1865 cell_ptr_.clear();
1866 cell_type_.clear();
1867
1/2
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
 4 if(steal_args) {
1868 4 cell_corners_.corner_vertex_.swap(tets);
1869 } else {
1870 ✗ cell_corners_.corner_vertex_ = tets;
1871 }
1872 4 resize_store(nb_tets);
1873 4 cell_corners_.resize_store(nb_tets*4);
1874 4 cell_facets_.resize_store(nb_tets*4);
1875 4 cell_facets_.adjacent_cell_.assign(
1876 nb_tets*4, NO_CELL
1877 );
1878 4 attributes().zero();
1879 4 cell_corners_.attributes().zero();
1880 4 cell_facets_.attributes().zero();
1881 4 }
1882
1883 ✗ void MeshCells::pop() {
1884 geo_debug_assert(nb() != 0);
1885 ✗ index_t corners_facets_new_size = cell_ptr_[nb()-1];
1886 ✗ cell_corners_.resize_store(corners_facets_new_size);
1887 ✗ cell_facets_.resize_store(corners_facets_new_size);
1888 ✗ resize_store(nb()-1);
1889 ✗ }
1890
1891 /**************************************************************************/
1892
1893 819 Mesh::Mesh(index_t dimension, bool single_precision)
1894 819 : vertices(*this),
1895
1/2
✓ Branch 1 taken 819 times.
✗ Branch 2 not taken.
 819 edges(*this),
1896
1/2
✓ Branch 1 taken 819 times.
✗ Branch 2 not taken.
 819 facets(*this),
1897
1/2
✓ Branch 1 taken 819 times.
✗ Branch 2 not taken.
 819 facet_corners(*this),
1898
1/2
✓ Branch 1 taken 819 times.
✗ Branch 2 not taken.
 819 cells(*this),
1899
1/2
✓ Branch 1 taken 819 times.
✗ Branch 2 not taken.
 819 cell_corners(*this),
1900
1/2
✓ Branch 2 taken 819 times.
✗ Branch 3 not taken.
 1638 cell_facets(*this)
1901 {
1902
1/2
✓ Branch 1 taken 819 times.
✗ Branch 2 not taken.
 819 vertices.bind_point_attribute(dimension, single_precision);
1903 819 }
1904
1905 1638 Mesh::~Mesh() {
1906 4914 }
1907
1908 357 void Mesh::clear(bool keep_attributes, bool keep_memory) {
1909 357 vertices.clear(keep_attributes, keep_memory);
1910 357 edges.clear(keep_attributes, keep_memory);
1911 357 facets.clear(keep_attributes, keep_memory);
1912 357 cells.clear(keep_attributes, keep_memory);
1913 357 }
1914
1915 60 void Mesh::copy(
1916 const Mesh& rhs,
1917 bool copy_attributes,
1918 MeshElementsFlags what
1919 ) {
1920
1/2
✓ Branch 0 taken 60 times.
✗ Branch 1 not taken.
 60 if(&rhs == this) {
1921 return;
1922 }
1923
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 60 times.
 60 if((what & MESH_VERTICES) == 0) {
1924 ✗ clear(false,false);
1925 ✗ return;
1926 }
1927 60 vertices.copy(rhs.vertices, copy_attributes);
1928
1/2
✓ Branch 0 taken 60 times.
✗ Branch 1 not taken.
 60 if((what & MESH_EDGES) != 0) {
1929 edges.copy(rhs.edges, copy_attributes);
1930 } else {
1931 ✗ edges.clear(false,false);
1932 }
1933
1/2
✓ Branch 0 taken 60 times.
✗ Branch 1 not taken.
 60 if((what & MESH_FACETS) != 0) {
1934 facets.copy(rhs.facets, copy_attributes);
1935 60 facet_corners.copy(rhs.facet_corners, copy_attributes);
1936 } else {
1937 ✗ facets.clear(false,false);
1938 }
1939
1/2
✓ Branch 0 taken 60 times.
✗ Branch 1 not taken.
 60 if((what & MESH_CELLS) != 0) {
1940 60 cells.copy(rhs.cells, copy_attributes);
1941 cell_corners.copy(rhs.cell_corners, copy_attributes);
1942 cell_facets.copy(rhs.cell_facets, copy_attributes);
1943 } else {
1944 ✗ cells.clear(false,false);
1945 }
1946 }
1947
1948 378 void Mesh::show_stats(const std::string& tag) const {
1949 index_t nb_borders = 0;
1950
2/2
✓ Branch 0 taken 3366114 times.
✓ Branch 1 taken 378 times.
 3366492 for(index_t co = 0; co < facet_corners.nb(); ++co) {
1951
2/2
✓ Branch 0 taken 33140 times.
✓ Branch 1 taken 3332974 times.
 3366114 if(facet_corners.adjacent_facet(co) == NO_FACET) {
1952 33140 nb_borders++;
1953 }
1954 }
1955
1956 378 Logger::out(tag)
1957 << (vertices.single_precision() ? "(FP32)" : "(FP64)")
1958 378 << " nb_v:" << vertices.nb()
1959 << " nb_e:" << edges.nb()
1960 << " nb_f:" << facets.nb()
1961 << " nb_b:" << nb_borders
1962 << " tri:" << facets.are_simplices()
1963 << " dim:" << vertices.dimension()
1964 << std::endl;
1965
1966
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 370 times.
 378 if(cells.nb() != 0) {
1967
1/2
✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
 8 if(cells.are_simplices()) {
1968 8 Logger::out(tag) << " nb_tets:"
1969 << cells.nb() << std::endl;
1970 } else {
1971
1972 index_t nb_cells_by_type[GEO::MESH_NB_CELL_TYPES];
1973 ✗ for(index_t i=0; i<GEO::MESH_NB_CELL_TYPES; ++i) {
1974 ✗ nb_cells_by_type[i] = 0;
1975 }
1976
1977 ✗ for(index_t c=0; c<cells.nb(); ++c) {
1978 geo_debug_assert(cells.type(c) < GEO::MESH_NB_CELL_TYPES);
1979 ✗ ++nb_cells_by_type[cells.type(c)];
1980 }
1981
1982 ✗ Logger::out(tag) << " Hybrid - nb_cells:"
1983 << cells.nb() << " "
1984 ✗ << " Tet:" << nb_cells_by_type[0]
1985 ✗ << " Hex:" << nb_cells_by_type[1]
1986 ✗ << " Psm:" << nb_cells_by_type[2]
1987 ✗ << " Pmd:" << nb_cells_by_type[3]
1988 ✗ << " Cnx:" << nb_cells_by_type[4]
1989 << std::endl;
1990 }
1991 }
1992
1993
1/2
✓ Branch 2 taken 378 times.
✗ Branch 3 not taken.
 378 display_attributes(tag, "vertices", vertices);
1994
1/2
✓ Branch 2 taken 378 times.
✗ Branch 3 not taken.
 378 display_attributes(tag, "edges", edges);
1995
1/2
✓ Branch 2 taken 378 times.
✗ Branch 3 not taken.
 378 display_attributes(tag, "facets", facets);
1996
1/2
✓ Branch 2 taken 378 times.
✗ Branch 3 not taken.
 378 display_attributes(tag, "facet_corners", facet_corners);
1997
1/2
✓ Branch 2 taken 378 times.
✗ Branch 3 not taken.
 378 display_attributes(tag, "cells", cells);
1998
1/2
✓ Branch 2 taken 378 times.
✗ Branch 3 not taken.
 378 display_attributes(tag, "cell_corners", cell_corners);
1999
1/2
✓ Branch 2 taken 378 times.
✗ Branch 3 not taken.
 378 display_attributes(tag, "cell_facets", cell_facets);
2000 378 }
2001
2002 ✗ void Mesh::assert_is_valid() {
2003 ✗ for(index_t f=0; f<facets.nb(); ++f) {
2004 ✗ for(
2005 index_t c=facets.corners_begin(f);
2006 ✗ c<facets.corners_end(f); ++c
2007 ) {
2008 ✗ geo_assert(facet_corners.vertex(c) < vertices.nb());
2009 index_t f2 = facet_corners.adjacent_facet(c);
2010 ✗ geo_assert(f2 == NO_FACET || f2 < facets.nb());
2011 }
2012 }
2013
2014 ✗ for(index_t c=0; c<cells.nb(); ++c) {
2015 ✗ for(index_t lv=0; lv<cells.nb_vertices(c); ++lv) {
2016 ✗ geo_assert(cells.vertex(c,lv) < vertices.nb());
2017 }
2018 ✗ for(index_t lf=0; lf<cells.nb_facets(c); ++lf) {
2019 index_t c2 = cells.adjacent(c,lf);
2020 ✗ geo_assert(c2 == NO_CELL || c2 < cells.nb());
2021 }
2022 }
2023 ✗ }
2024
2025
2/2
✓ Branch 0 taken 611 times.
✓ Branch 1 taken 2035 times.
 2646 void Mesh::display_attributes(
2026 const std::string& tag, const std::string& subelement_name,
2027 const MeshSubElementsStore& subelements
2028 ) const {
2029
2/2
✓ Branch 0 taken 611 times.
✓ Branch 1 taken 2035 times.
 2646 if(subelements.attributes().nb() != 0) {
2030 vector<std::string> names;
2031
1/2
✓ Branch 1 taken 611 times.
✗ Branch 2 not taken.
 611 subelements.attributes().list_attribute_names(names);
2032 std::string names_str;
2033
2/2
✓ Branch 0 taken 631 times.
✓ Branch 1 taken 611 times.
 2484 for(index_t i=0; i<names.size(); ++i) {
2034
2/2
✓ Branch 0 taken 20 times.
✓ Branch 1 taken 611 times.
 631 if(i != 0) {
2035
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 40 names_str = names_str + ",";
2036 }
2037
1/2
✓ Branch 2 taken 631 times.
✗ Branch 3 not taken.
 1262 names_str = names_str + names[i];
2038 AttributeStore* store =
2039
1/2
✓ Branch 1 taken 631 times.
✗ Branch 2 not taken.
 631 subelements.attributes().find_attribute_store(names[i]);
2040 631 index_t dim = store->dimension();
2041
2/2
✓ Branch 0 taken 532 times.
✓ Branch 1 taken 99 times.
 631 if(dim != 1) {
2042
1/4
✓ Branch 1 taken 532 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 1064 names_str += ("[" + String::to_string(dim) + "]");
2043 }
2044 }
2045
1/2
✓ Branch 1 taken 611 times.
✗ Branch 2 not taken.
 611 Logger::out(tag) << "Attributes on " << subelement_name
2046 << ": " << names_str << std::endl;
2047 }
2048 2646 }
2049
2050 ✗ index_t Mesh::nb_subelements_types() const {
2051 ✗ return 7;
2052 }
2053
2054 ✗ MeshSubElementsStore& Mesh::get_subelements_by_index(
2055 index_t i
2056 ) {
2057 ✗ switch(i) {
2058 ✗ case 0:
2059 ✗ return vertices;
2060 ✗ case 1:
2061 ✗ return edges;
2062 ✗ case 2:
2063 ✗ return facets;
2064 ✗ case 3:
2065 ✗ return facet_corners;
2066 ✗ case 4:
2067 ✗ return cells;
2068 ✗ case 5:
2069 ✗ return cell_corners;
2070 ✗ case 6:
2071 ✗ return cell_facets;
2072 default:
2073 ✗ geo_assert_not_reached;
2074 }
2075 }
2076
2077 ✗ const MeshSubElementsStore& Mesh::get_subelements_by_index(
2078 index_t i
2079 ) const {
2080 ✗ switch(i) {
2081 ✗ case 0:
2082 ✗ return vertices;
2083 ✗ case 1:
2084 ✗ return edges;
2085 ✗ case 2:
2086 ✗ return facets;
2087 ✗ case 3:
2088 ✗ return facet_corners;
2089 ✗ case 4:
2090 ✗ return cells;
2091 ✗ case 5:
2092 ✗ return cell_corners;
2093 ✗ case 6:
2094 ✗ return cell_facets;
2095 default:
2096 ✗ geo_assert_not_reached;
2097 }
2098 }
2099
2100 ✗ MeshSubElementsStore& Mesh::get_subelements_by_type(
2101 MeshElementsFlags what
2102 ) {
2103 ✗ switch(what) {
2104 ✗ case MESH_VERTICES:
2105 ✗ return vertices;
2106 ✗ case MESH_EDGES:
2107 ✗ return edges;
2108 ✗ case MESH_FACETS:
2109 ✗ return facets;
2110 ✗ case MESH_FACET_CORNERS:
2111 ✗ return facet_corners;
2112 ✗ case MESH_CELLS:
2113 ✗ return cells;
2114 ✗ case MESH_CELL_CORNERS:
2115 ✗ return cell_corners;
2116 ✗ case MESH_CELL_FACETS:
2117 ✗ return cell_facets;
2118 case MESH_NONE:
2119 case MESH_ALL_ELEMENTS:
2120 case MESH_ALL_SUBELEMENTS:
2121 ✗ geo_assert_not_reached;
2122 }
2123 return *(MeshSubElementsStore*)nullptr;
2124 }
2125
2126 ✗ const MeshSubElementsStore& Mesh::get_subelements_by_type(
2127 MeshElementsFlags what
2128 ) const {
2129 ✗ switch(what) {
2130 ✗ case MESH_VERTICES:
2131 ✗ return vertices;
2132 ✗ case MESH_EDGES:
2133 ✗ return edges;
2134 ✗ case MESH_FACETS:
2135 ✗ return facets;
2136 ✗ case MESH_FACET_CORNERS:
2137 ✗ return facet_corners;
2138 ✗ case MESH_CELLS:
2139 ✗ return cells;
2140 ✗ case MESH_CELL_CORNERS:
2141 ✗ return cell_corners;
2142 ✗ case MESH_CELL_FACETS:
2143 ✗ return cell_facets;
2144 case MESH_NONE:
2145 case MESH_ALL_ELEMENTS:
2146 case MESH_ALL_SUBELEMENTS:
2147 ✗ geo_assert_not_reached;
2148 }
2149 return *(MeshSubElementsStore*)nullptr;
2150 }
2151
2152 ✗ std::string Mesh::subelements_type_to_name(MeshElementsFlags what) {
2153 std::string result;
2154 ✗ switch(what) {
2155 case MESH_VERTICES:
2156 result = "vertices";
2157 break;
2158 case MESH_EDGES:
2159 result = "edges";
2160 break;
2161 case MESH_FACETS:
2162 result = "facets";
2163 break;
2164 case MESH_FACET_CORNERS:
2165 result = "facet_corners";
2166 break;
2167 case MESH_CELLS:
2168 result = "cells";
2169 break;
2170 case MESH_CELL_CORNERS:
2171 result = "cell_corners";
2172 break;
2173 case MESH_CELL_FACETS:
2174 result = "cell_facets";
2175 break;
2176 case MESH_NONE:
2177 case MESH_ALL_ELEMENTS:
2178 case MESH_ALL_SUBELEMENTS:
2179 ✗ geo_assert_not_reached;
2180 }
2181 ✗ return result;
2182 }
2183
2184 ✗ MeshElementsFlags Mesh::name_to_subelements_type(const std::string& name) {
2185 ✗ if(name == "vertices") {
2186 return MESH_VERTICES;
2187 ✗ } else if(name == "edges") {
2188 return MESH_EDGES;
2189 ✗ } else if(name == "facets") {
2190 return MESH_FACETS;
2191 ✗ } else if(name == "facet_corners") {
2192 return MESH_FACET_CORNERS;
2193 ✗ } else if(name == "cells") {
2194 return MESH_CELLS;
2195 ✗ } else if(name == "cell_corners") {
2196 return MESH_CELL_CORNERS;
2197 ✗ } else if(name == "cell_facets") {
2198 ✗ return MESH_CELL_FACETS;
2199 }
2200 return MESH_NONE;
2201 }
2202
2203 /**************************************************************************/
2204
2205 ✗ bool Mesh::parse_attribute_name(
2206 const std::string& full_attribute_name,
2207 MeshElementsFlags& where,
2208 std::string& attribute_name,
2209 index_t& component
2210 ) {
2211
2212 ✗ size_t pos1 = full_attribute_name.find('.');
2213 ✗ if(pos1 == std::string::npos) {
2214 return false;
2215 }
2216
2217 {
2218 ✗ std::string where_name = full_attribute_name.substr(0,pos1);
2219 ✗ where = Mesh::name_to_subelements_type(where_name);
2220 ✗ if(where == MESH_NONE) {
2221 return false;
2222 }
2223 }
2224
2225 ✗ attribute_name = full_attribute_name.substr(
2226 ✗ pos1+1, full_attribute_name.length()-pos1-1
2227 ✗ );
2228
2229 ✗ size_t pos2 = attribute_name.find('[');
2230 ✗ if(pos2 == std::string::npos) {
2231 ✗ component = 0;
2232 } else {
2233 ✗ if(attribute_name[attribute_name.length()-1] != ']') {
2234 ✗ return false;
2235 }
2236 std::string component_str = attribute_name.substr(
2237 ✗ pos2+1, attribute_name.length()-pos2-2
2238 ✗ );
2239 ✗ attribute_name = attribute_name.substr(0, pos2);
2240 try {
2241 ✗ component = String::to_uint(component_str);
2242 ✗ } catch(...) {
2243 return false;
2244 ✗ }
2245 }
2246
2247 return true;
2248 }
2249
2250 /**************************************************************************/
2251 }
2252
2253 namespace {
2254
2255 using namespace GEO;
2256
2257 /**
2258 * \brief Gets the names of all scalar attributes from an AttributeManager
2259 * \param[in] attributes a const reference to the attribute manager
2260 * \param[in] prefix a const rerefenre to a string to be prepended to
2261 * all attribute names
2262 * \return a ';'-separated list of all the scalar attributes
2263 */
2264 ✗ std::string get_scalar_attributes_impl(
2265 const AttributesManager& attributes,
2266 const std::string& prefix
2267 ) {
2268 std::string result;
2269 vector<std::string> attribute_names;
2270 ✗ attributes.list_attribute_names(attribute_names);
2271
2272 ✗ for(index_t i=0; i<attribute_names.size(); ++i) {
2273 const AttributeStore* store = attributes.
2274 ✗ find_attribute_store(attribute_names[i]);
2275 ✗ if(ReadOnlyScalarAttributeAdapter::can_be_bound_to(store)) {
2276 index_t dim =
2277 ✗ ReadOnlyScalarAttributeAdapter::nb_scalar_elements_per_item(
2278 store
2279 );
2280 ✗ if(dim == 1) {
2281 ✗ if(result != "") {
2282 result += ";";
2283 }
2284 ✗ result += prefix + "." + attribute_names[i];
2285 } else {
2286 ✗ for(index_t j=0; j<dim; ++j) {
2287 ✗ if(result != "") {
2288 result += ";";
2289 }
2290 result +=
2291 ✗ prefix + "." + attribute_names[i] +
2292 ✗ "[" + String::to_string(j) + "]";
2293 }
2294 }
2295 }
2296 }
2297 ✗ return result;
2298 }
2299
2300 /**
2301 * \brief Gets the names of all attributes from an AttributeManager
2302 * \param[in] attributes a const reference to the attribute manager
2303 * \param[in] prefix a const rerefenre to a string to be prepended to
2304 * all attribute names
2305 * \return a ';'-separated list of all the attributes
2306 */
2307 ✗ std::string get_attributes_impl(
2308 const AttributesManager& attributes,
2309 const std::string& prefix
2310 ) {
2311 std::string result;
2312 vector<std::string> attribute_names;
2313 ✗ attributes.list_attribute_names(attribute_names);
2314
2315 ✗ for(index_t i=0; i<attribute_names.size(); ++i) {
2316 ✗ if(result != "") {
2317 result += ";";
2318 }
2319 ✗ result += prefix + "." + attribute_names[i];
2320 }
2321 ✗ return result;
2322 }
2323
2324 /**
2325 * \brief Gets the names of all vector attributes from an AttributeManager
2326 * \param[in] attributes a const reference to the attribute manager
2327 * \param[in] prefix a const rerefenre to a string to be prepended to
2328 * all attribute names
2329 * \param[in] max_dim if non-zero, only return vector attributes with
2330 * dimension lower than max_dim
2331 * \return a ';'-separated list of all the vector attributes
2332 */
2333 ✗ std::string get_vector_attributes_impl(
2334 const AttributesManager& attributes,
2335 const std::string& prefix,
2336 index_t max_dim = 0
2337 ) {
2338 std::string result;
2339 vector<std::string> attribute_names;
2340 ✗ attributes.list_attribute_names(attribute_names);
2341
2342 ✗ for(index_t i=0; i<attribute_names.size(); ++i) {
2343 const AttributeStore* store = attributes.
2344 ✗ find_attribute_store(attribute_names[i]);
2345 if(
2346 ✗ store->dimension() >= 2 &&
2347 ✗ (max_dim == 0 || store->dimension() <= max_dim))
2348 {
2349 ✗ if(result != "") {
2350 result += ";";
2351 }
2352 ✗ result += prefix + "." + attribute_names[i];
2353 }
2354 ✗ if(
2355 ✗ store->elements_type_matches(typeid(vec2).name()) &&
2356 (max_dim == 0 || 2 <= max_dim)
2357 ) {
2358 ✗ if(result != "") {
2359 result += ";";
2360 }
2361 ✗ result += prefix + "." + attribute_names[i];
2362 }
2363 ✗ if(
2364 ✗ store->elements_type_matches(typeid(vec3).name()) &&
2365 ✗ (max_dim == 0 || 3 <= max_dim)
2366 ) {
2367 ✗ if(result != "") {
2368 result += ";";
2369 }
2370 ✗ result += prefix + "." + attribute_names[i];
2371 }
2372 }
2373 ✗ return result;
2374 }
2375
2376
2377 /**
2378 * \brief Appends a string to another one, with ';' delimiters.
2379 * \details If a is non-empty, a ';' delimiter is inserted.
2380 * \param[in,out] a a string
2381 * \param[in] b a string to be appended to a
2382 */
2383 ✗ static void strappend(std::string& a, const std::string& b) {
2384 ✗ if(b != "") {
2385 ✗ if(a != "") {
2386 a += ";";
2387 }
2388 a += b;
2389 }
2390 ✗ }
2391 }
2392
2393 namespace GEO {
2394
2395 ✗ std::string Mesh::get_attributes() const {
2396 std::string result;
2397 ✗ strappend(
2398 ✗ result,get_attributes_impl(vertices.attributes(),"vertices")
2399 );
2400 ✗ strappend(
2401 ✗ result,get_attributes_impl(edges.attributes(),"edges")
2402 );
2403 ✗ strappend(
2404 ✗ result,get_attributes_impl(facets.attributes(),"facets")
2405 );
2406 ✗ strappend(
2407 ✗ result,get_attributes_impl(
2408 ✗ facet_corners.attributes(),"facet_corners"
2409 )
2410 );
2411 ✗ strappend(
2412 ✗ result,get_attributes_impl(cells.attributes(),"cells")
2413 );
2414 ✗ strappend(
2415 ✗ result,get_attributes_impl(
2416 ✗ cell_corners.attributes(),"cell_corners"
2417 )
2418 );
2419 ✗ strappend(result,get_attributes_impl(
2420 ✗ cell_facets.attributes(),"cell_facets")
2421 );
2422 ✗ return result;
2423 }
2424
2425 ✗ std::string Mesh::get_scalar_attributes() const {
2426 std::string result;
2427 ✗ strappend(
2428 ✗ result,get_scalar_attributes_impl(vertices.attributes(),"vertices")
2429 );
2430 ✗ strappend(
2431 ✗ result,get_scalar_attributes_impl(edges.attributes(),"edges")
2432 );
2433 ✗ strappend(
2434 ✗ result,get_scalar_attributes_impl(facets.attributes(),"facets")
2435 );
2436 ✗ strappend(result,get_scalar_attributes_impl(
2437 ✗ facet_corners.attributes(),"facet_corners"
2438 )
2439 );
2440 ✗ strappend(
2441 ✗ result,get_scalar_attributes_impl(cells.attributes(),"cells")
2442 );
2443 ✗ strappend(
2444 ✗ result,get_scalar_attributes_impl(
2445 ✗ cell_corners.attributes(),"cell_corners"
2446 )
2447 );
2448 ✗ strappend(result,get_scalar_attributes_impl(
2449 ✗ cell_facets.attributes(),"cell_facets")
2450 );
2451 ✗ return result;
2452 }
2453
2454
2455 ✗ std::string Mesh::get_vector_attributes(index_t max_dim) const {
2456 std::string result;
2457 ✗ strappend(
2458 ✗ result,get_vector_attributes_impl(vertices.attributes(),"vertices",max_dim)
2459 );
2460 ✗ strappend(
2461 ✗ result,get_vector_attributes_impl(edges.attributes(),"edges",max_dim)
2462 );
2463 ✗ strappend(
2464 ✗ result,get_vector_attributes_impl(facets.attributes(),"facets",max_dim)
2465 );
2466 ✗ strappend(result,get_vector_attributes_impl(
2467 ✗ facet_corners.attributes(),"facet_corners",max_dim
2468 )
2469 );
2470 ✗ strappend(
2471 ✗ result,get_vector_attributes_impl(cells.attributes(),"cells",max_dim)
2472 );
2473 ✗ strappend(
2474 ✗ result,get_vector_attributes_impl(
2475 ✗ cell_corners.attributes(),"cell_corners",max_dim
2476 )
2477 );
2478 ✗ strappend(result,get_vector_attributes_impl(
2479 ✗ cell_facets.attributes(),"cell_facets",max_dim)
2480 );
2481 ✗ return result;
2482 }
2483
2484
2485 }
2486