GCC Code Coverage Report

Directory: ./
File: lib/geogram/mesh/mesh.h
Date: 2026-09-07 02:36:43
Exec Total Coverage
Lines: 510 716 71.2%
Functions: 127 162 78.4%
Branches: 224 874 25.6%
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 #ifndef GEOGRAM_MESH_MESH
41 #define GEOGRAM_MESH_MESH
42
43 #include <geogram/basic/common.h>
44 #include <geogram/basic/range.h>
45 #include <geogram/basic/attributes.h>
46 #include <geogram/basic/vector_attribute.h>
47 #include <geogram/basic/geometry.h>
48 #include <tuple>
49
50 // GOMGEN (Graphite's pre-processor) does not understand
51 // modern C++ (but it is not a big drama since it is just
52 // used to generate the GUI).
53 #ifdef GOMGEN
54 #define MESH_NO_SYNTAXIC_SUGAR
55 #endif
56
57 /**
58 * \file geogram/mesh/mesh.h
59 * \brief The class that represents a mesh.
60 */
61
62 namespace GEO {
63
64 class Mesh;
65
66 static constexpr index_t NO_VERTEX = NO_INDEX;
67 static constexpr index_t NO_EDGE = NO_INDEX;
68 static constexpr index_t NO_FACET = NO_INDEX;
69 static constexpr index_t NO_CELL = NO_INDEX;
70 static constexpr index_t NO_CORNER = NO_INDEX;
71
72 /**
73 * \brief Base class for mesh sub-element storage.
74 * \details Sub-elements are those that cannot exist
75 * independently (such as MeshFacetCorner, MeshCellCorner
76 * and MeshCellFacet).
77 * \relates Mesh
78 */
79 class GEOGRAM_API MeshSubElementsStore {
80 public:
81
82 /**
83 * \brief Constructs a new MeshSubElementStore.
84 * \param[in] mesh a reference to the mesh
85 * this MeshElementStore belongs to.
86 */
87 MeshSubElementsStore(Mesh& mesh);
88
89 /**
90 * \brief MeshElementStore destructor.
91 */
92 virtual ~MeshSubElementsStore();
93
94 /**
95 * \brief Gets the number of (sub-)elements.
96 * \return the number of (sub-)elements in this store
97 */
98 1702709470 index_t nb() const {
99 1702709470 return nb_;
100 }
101
102 /**
103 * \brief Gets the attributes manager.
104 * \details The returned reference is not a const one,
105 * so that attributes can be bound / unbound / accessed
106 * even if the mesh is const.
107 * \return a modifiable reference to the attributes manager.
108 */
109 2756599 AttributesManager& attributes() const {
110 2756599 return attributes_;
111 }
112
113 /**
114 * \brief Used by range-based for.
115 * \return The index of the first position.
116 */
117 95155 index_as_iterator begin() const {
118 95155 return index_as_iterator(0);
119 }
120
121 /**
122 * \brief Used by range-based for.
123 * \return The index of one position past the last position.
124 */
125 95155 index_as_iterator end() const {
126 95155 return index_as_iterator(nb());
127 }
128
129 protected:
130
131 /**
132 * \brief Removes all the elements and attributes.
133 * \param[in] keep_attributes if true, then all the
134 * existing attribute names / bindings are kept (but
135 * they are cleared). If false, they are destroyed.
136 * \param[in] keep_memory if true, then memory is
137 * kept and can be reused by subsequent mesh
138 * element creations.
139 */
140 virtual void clear_store(
141 bool keep_attributes, bool keep_memory = false
142 );
143
144 /**
145 * \brief Resizes this MeshSubElementsStore.
146 * \details On exit, nb() == new_size, elements are
147 * created or destroyed if needed.
148 * \param[in] new_size the desired size
149 */
150 virtual void resize_store(index_t new_size);
151
152
153 /**
154 * \brief Reserves space for new elements
155 * \param[in] nb_to_reserve the number of subelements to reserve
156 */
157 114 void reserve_store(index_t nb_to_reserve) {
158 114 index_t nb = this->nb();
159 114 resize_store(nb + nb_to_reserve);
160 114 resize_store(nb);
161 114 }
162
163 /**
164 * \brief Creates a contiguous chunk of attributes for sub-elements.
165 * \param[in] nb number of sub-elements to create
166 * \return the index of the first created sub-element
167 */
168 1182 index_t create_sub_elements(index_t nb) {
169 1182 index_t result = nb_;
170
2/2
✓ Branch 1 taken 1123 times.
✓ Branch 2 taken 59 times.
 1182 if(nb_ + nb > attributes_.size()) {
171 1123 index_t new_capacity=nb_ + nb;
172
2/2
✓ Branch 0 taken 639 times.
✓ Branch 1 taken 484 times.
 1123 if(nb < 128) {
173 639 new_capacity = std::max(index_t(16),attributes_.capacity());
174
2/2
✓ Branch 0 taken 528 times.
✓ Branch 1 taken 639 times.
 1167 while(new_capacity < nb_ + nb) {
175 528 new_capacity *= 2;
176 }
177 }
178 1123 attributes_.reserve(new_capacity);
179 }
180 1182 nb_ += nb;
181 1182 attributes_.resize(nb_);
182 1182 return result;
183 }
184
185 /**
186 * \brief Creates attributes for a sub-element
187 * \return the index of the created element
188 */
189 6136740 index_t create_sub_element() {
190 6136740 index_t result = nb_;
191 6136740 ++nb_;
192
2/2
✓ Branch 1 taken 6909 times.
✓ Branch 2 taken 6129831 times.
 6136740 if(attributes_.capacity() < nb_) {
193 index_t new_capacity =
194 6909 std::max(index_t(16),attributes_.capacity()*2);
195 6909 attributes_.reserve(new_capacity);
196 }
197 6136740 attributes_.resize(nb_);
198 6136740 return result;
199 }
200
201 /**
202 * \brief Makes the size of the store tightly match
203 * the number of the elements.
204 * \details When elements are created one by one, the
205 * system may allocate more memory than necessary, to
206 * amortize the cost of container reallocation. Once
207 * the object is constructed, this function may be called
208 * to release the memory that was not used.
209 */
210 void adjust_store() {
211 attributes_.resize(nb_);
212 }
213
214 /**
215 * \brief Copies a MeshSubElementsStore into
216 * this one.
217 * \param[in] rhs a const reference to the
218 * MeshSubElementsStore to be copied.
219 * \param[in] copy_attributes if true, copies
220 * also the attributes, else attributes are
221 * cleared.
222 */
223 420 void copy(
224 const MeshSubElementsStore& rhs,
225 bool copy_attributes = true
226 ) {
227 420 nb_ = rhs.nb();
228
1/2
✓ Branch 0 taken 420 times.
✗ Branch 1 not taken.
 420 if(copy_attributes) {
229 420 attributes_.copy(rhs.attributes_);
230 } else {
231 ✗ attributes_.clear(false,false);
232 ✗ attributes_.resize(rhs.attributes_.size());
233 }
234 420 }
235
236 protected:
237 Mesh& mesh_;
238 mutable AttributesManager attributes_;
239 index_t nb_;
240 };
241
242
243 /**************************************************************************/
244
245 /**
246 * \brief Base class for mesh elements.
247 * \details Mesh elements can be created / manipulated independantly,
248 * in contrast with sub-elements that cannot. Mesh elements are
249 * vertices, facets and cells.
250 * \relates Mesh
251 */
252 class GEOGRAM_API MeshElements {
253 public:
254 MeshElements();
255 virtual ~MeshElements();
256
257 /**
258 * \brief Deletes a set of elements.
259 * \param[in] to_delete a vector of size nb(). If to_delete[e]
260 * is different from 0, then element e will be destroyed, else
261 * it will be kept. On exit, to_delete is modified (it is used
262 * for internal bookkeeping).
263 * \param[in] remove_isolated_vertices if true, then the vertices
264 * that are no longer incident to any element are deleted.
265 */
266 virtual void delete_elements(
267 vector<index_t>& to_delete,
268 bool remove_isolated_vertices=true
269 ) = 0;
270
271 /**
272 * \brief Applies a permutation to the elements and their attributes.
273 * \details On exit, permutation is modified (used for internal
274 * bookkeeping). Applying a permutation \p permutation is equivalent
275 * to:
276 * \code
277 * for(i=0; i<permutation.size(); i++) {
278 * data2[i] = data[permutation[i]]
279 * }
280 * data = data2 ;
281 * \endcode
282 */
283 virtual void permute_elements(vector<index_t>& permutation) = 0;
284
285 /**
286 * \brief Removes all the elements and attributes.
287 * \param[in] keep_attributes if true, then all the
288 * existing attribute names / bindings are kept (but
289 * they are cleared). If false, they are destroyed.
290 * \param[in] keep_memory if true, then memory is
291 * kept and can be reused by subsequent mesh
292 * element creations.
293 */
294 virtual void clear(
295 bool keep_attributes=true, bool keep_memory=false
296 ) = 0;
297
298 /**
299 * \brief Removes the last element.
300 */
301 virtual void pop() = 0;
302
303 protected:
304 /**
305 * \brief Tests whether a vector contains a non-zero value.
306 * \details This function is used internally by delete_elements()
307 * \param[in] I a vector of signed integers
308 * \retval true if \p I contains at least a non-zero value
309 * \retval false otherwise
310 */
311 762 static bool has_non_zero(const GEO::vector<index_t>& I) {
312
2/2
✓ Branch 1 taken 567380 times.
✓ Branch 2 taken 379 times.
 567759 for(index_t i = 0; i < I.size(); i++) {
313
2/2
✓ Branch 1 taken 383 times.
✓ Branch 2 taken 566997 times.
 567380 if(I[i] != 0) {
314 383 return true;
315 }
316 }
317 379 return false;
318 }
319 };
320
321 /**************************************************************************/
322
323 class MeshEdges;
324 class MeshFacetCornersStore;
325 class MeshCellCornersStore;
326
327 /**
328 * \brief The vertices of a mesh.
329 * \relates Mesh
330 */
331 class GEOGRAM_API MeshVertices :
332 public MeshSubElementsStore, public MeshElements {
333 public:
334 MeshVertices(Mesh& mesh);
335 ~MeshVertices() override;
336
337 /**
338 * \brief Removes the vertices that have no mesh element
339 * incident to them.
340 */
341 void remove_isolated();
342
343 void delete_elements(
344 vector<index_t>& to_delete, bool remove_isolated_vertices=true
345 ) override;
346
347 void permute_elements(vector<index_t>& permutation) override;
348
349 /**
350 * \brief Creates a new vertex
351 * \return the index of the created vertex
352 */
353 674731 index_t create_vertex() {
354 674731 return MeshSubElementsStore::create_sub_element();
355 }
356
357 /**
358 * \brief Creates a new vertex
359 * \param[in] coords a pointer to dimension() coordinates
360 * \return the index of the created vertex
361 */
362 120062 index_t create_vertex(const double* coords) {
363 // Sanity check:
364 // It is not correct to call create_vertex(point_ptr(v)) since
365 // create_vertex may realloc the points coordinates vector, thus
366 // invalidate point_ptr(v).
367
6/12
✓ Branch 1 taken 119828 times.
✓ Branch 2 taken 234 times.
✓ Branch 4 taken 29655 times.
✓ Branch 5 taken 90173 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 29655 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 120062 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
 120062 geo_debug_assert(
368 nb() == 0 ||
369 coords < point_ptr(0) ||
370 coords >= point_ptr(0) + nb() * dimension()
371 );
372 120062 index_t result = create_vertex();
373
2/2
✓ Branch 1 taken 356933 times.
✓ Branch 2 taken 120062 times.
 476995 for(index_t c=0; c<dimension(); ++c) {
374 356933 point_ptr(result)[c] = coords[c];
375 }
376 120062 return result;
377 }
378
379 /**
380 * \brief Creates a vertex from a 3d point
381 * \param[in] p a const reference to the 3d point
382 * \return the index of the created vertex
383 * \pre dimension() == 3
384 */
385 7672 template <index_t DIM> index_t create_vertex(
386 const vecng<DIM,double>& p
387 ) {
388
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 4179 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 7672 geo_debug_assert(dimension() == DIM);
389 7672 return create_vertex(p.data());
390 }
391
392 /**
393 * \brief Creates a contiguous chunk of vertices.
394 * \param[in] nb number of sub-elements to create
395 * \return the index of the first created vertex
396 */
397 386 index_t create_vertices(index_t nb) {
398 386 return MeshSubElementsStore::create_sub_elements(nb);
399 }
400
401 void clear(
402 bool keep_attributes=true, bool keep_memory=false
403 ) override;
404
405 /**
406 * \brief Sets single precision mode.
407 * \details Single-precision mode is used for instance
408 * by Vorpaview, to make it more memory efficient.
409 * Existing point coordinates are copied and converted to
410 * single precision.
411 */
412 void set_single_precision();
413
414 /**
415 * \brief Sets double precision mode.
416 * \details Double precision mode is the default.
417 * Single-precision mode is used for instance
418 * by Vorpaview, to make it more memory efficient.
419 * Existing point coordinates are copied and converted to
420 * double precision.
421 */
422 void set_double_precision();
423
424 /**
425 * \brief Tests whether vertices are stored in
426 * single-precision mode.
427 * \details Single-precision mode is used for instance
428 * by Vorpaview, to make it more memory efficient.
429 */
430 595078145 bool single_precision() const {
431 595078145 return point_fp32_.is_bound();
432 }
433
434 /**
435 * \brief Tests whether vertices are stored in
436 * double-precision mode.
437 * \details Double precision mode is the default.
438 * Single-precision mode is used for instance
439 * by Vorpaview, to make it more memory efficient.
440 */
441 3 bool double_precision() const {
442 3 return point_.is_bound();
443 }
444
445 /**
446 * \brief Gets the dimension of the vertices.
447 * \return the number of coordinates in each vertex
448 */
449 246138488 index_t dimension() const {
450 return
451
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 246138488 times.
 246138488 single_precision() ?
452 ✗ point_fp32_.dimension() :
453 246138488 point_.dimension() ;
454 }
455
456 /**
457 * \brief Sets the dimension of the vertices
458 * \details Existing coordinates are kept, newly created
459 * coordinates are initialized to zero.
460 * \param[in] dim new dimension
461 */
462 700 void set_dimension(index_t dim) {
463
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 700 times.
 700 if(single_precision()) {
464 ✗ point_fp32_.redim(dim);
465 } else {
466 700 point_.redim(dim);
467 }
468 700 }
469
470 /**
471 * \brief Gets a point
472 * \param[in] v the index of the vertex
473 * \return a const pointer to the coordinates of the point
474 * that correspond to the vertex
475 * \pre !single_precision()
476 */
477 103216137 const double* point_ptr(index_t v) const {
478
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 103216137 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 103216137 geo_debug_assert(v < nb());
479
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 103216137 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 103216137 geo_debug_assert(!single_precision());
480 103216137 return &point_[v*point_.dimension()];
481 }
482
483 /**
484 * \brief Gets a point
485 * \param[in] v the vertex, in 0..nb()-1
486 * \return a pointer to the coordinates of the point
487 * that correspond to the vertex
488 * \pre !single_precision()
489 */
490 2448452 double* point_ptr(index_t v) {
491
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 2448452 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 2448452 geo_debug_assert(v < nb());
492
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 2448452 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 2448452 geo_debug_assert(!single_precision());
493 2448452 return &point_[v*point_.dimension()];
494 }
495
496
497 /**
498 * \brief Gets a point
499 * \param[in] v the vertex, in 0..nb()-1
500 * \return a modifiable reference to the point
501 * that corresponds to the vertex
502 * \pre !single_precision()
503 */
504 262728377 template<index_t DIM=3> vecng<DIM,double>& point(index_t v) {
505
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 239216625 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 262728377 geo_debug_assert(v < nb());
506
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 239216625 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 262728377 geo_debug_assert(!single_precision());
507
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 239216625 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 262728377 geo_debug_assert(dimension() >= DIM);
508 262728377 return Memory::pointer_as_reference<vecng<DIM,double>>(
509 262728377 &point_[v*point_.dimension()]
510 262728377 );
511 }
512
513 /**
514 * \brief Gets a point
515 * \param[in] v the vertex, in 0..nb()-1
516 * \return a const reference to the point
517 * that corresponds to the vertex
518 * \pre !single_precision()
519 */
520 3186434 template <index_t DIM=3> const vecng<DIM,double>& point(
521 index_t v
522 ) const {
523
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 3078794 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 3186434 geo_debug_assert(v < nb());
524
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 3078794 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 3186434 geo_debug_assert(!single_precision());
525
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 3078794 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 3186434 geo_debug_assert(dimension() >= DIM);
526 3186434 return Memory::pointer_as_reference<vecng<DIM,double>>(
527 3186434 &point_[v*point_.dimension()]
528 3186434 );
529 }
530
531 /**
532 * \brief Gets a (single-precision) point
533 * \param[in] v the index of the vertex
534 * \return a const pointer to the coordinates
535 * of the point that corresponds to the vertex
536 * \pre single_precision()
537 */
538 ✗ const float* single_precision_point_ptr(index_t v) const {
539 ✗ geo_debug_assert(v < nb());
540 ✗ geo_debug_assert(single_precision());
541 ✗ return &point_fp32_[v*point_fp32_.dimension()];
542 }
543
544 /**
545 * \brief Gets a (single-precision) point
546 * \param[in] v the index of the vertex
547 * \return a pointer to the coordinates of the point
548 * that corresponds to the vertex
549 * \pre single_precision()
550 */
551 ✗ float* single_precision_point_ptr(index_t v) {
552 ✗ geo_debug_assert(v < nb());
553 ✗ geo_debug_assert(single_precision());
554 ✗ return &point_fp32_[v*point_fp32_.dimension()];
555 }
556
557 /**
558 * \brief Gets the coordinates of all the points as a single vector.
559 * \details The vector contains nb() * dimension() coordinates
560 * [x0, y0, z0, x1, y1, z1, ...]. It is forbidden to change the
561 * size of the returned vector; the reference is invalidated by
562 * create_vertices() (same lifetime rule as point_ptr()).
563 * \return a const reference to the vector of all point coordinates.
564 * \pre !single_precision()
565 */
566 const vector<double>& point_coordinates() const {
567 geo_debug_assert(!single_precision());
568 return point_.get_vector();
569 }
570
571 /**
572 * \brief Gets the coordinates of all the points as a single vector.
573 * \return a modifiable reference to the vector of all point
574 * coordinates (see the const overload for the contract).
575 * \pre !single_precision()
576 */
577 vector<double>& point_coordinates() {
578 geo_debug_assert(!single_precision());
579 return point_.get_vector();
580 }
581
582
583 /**
584 * \brief Assigns all the points.
585 * \param[in] points a vector that contains all the coordinates
586 * of the points
587 * \param[in] dim the dimension of the points, i.e. number of
588 * coordinates per point
589 * \param[in] steal_arg if true, memory is stolen from \p points,
590 * using std::vector::swap (no memory copy).
591 */
592 void assign_points(
593 vector<double>& points, index_t dim, bool steal_arg
594 );
595
596 /**
597 * \brief Assigns all the points.
598 * \param[in] points a const pointer to the (\p dim * \p nb_pts)
599 * coordinates of al the points
600 * \param[in] dim the dimension of the points, i.e. number of
601 * coordinates per point
602 * \param[in] nb_pts number of points
603 */
604 void assign_points(
605 const double* points, index_t dim, index_t nb_pts
606 );
607
608 void pop() override;
609
610 #ifndef MESH_NO_SYNTAXIC_SUGAR
611
612 /**
613 * \brief Gets the 3D points of the mesh as an iterable sequence
614 * \details Each point is returned as a const reference
615 */
616 76 template <index_t DIM = 3> auto points() const {
617 typedef vecng<DIM,double> vecn;
618
1/2
✓ Branch 1 taken 76 times.
✗ Branch 2 not taken.
 76 return transform_range_ref(
619 152 index_range(0, nb()),
620 229229 [this](index_t v)->const vecn& {
621 // for MSVC that cannot chose among const/non-const versions
622 229229 return Memory::pointer_as_reference<vecn>(point_ptr(v));
623 // return point<DIM>(v); // MSVC does not understand this one
624 }
625 152 );
626 }
627
628 /**
629 * \brief Gets the 3D points of the mesh as an iterable sequence
630 * \details Each point is returned as a modifiable reference
631 */
632 9 template <index_t DIM = 3> auto points() {
633 typedef vecng<DIM,double> vecn;
634
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 9 return transform_range_ref(
635 18 index_range(0, nb()),
636 107537 [this](index_t v)->vecn& {
637 // for MSVC that cannot chose among const/non-const versions
638 107537 return Memory::pointer_as_reference<vecn>(point_ptr(v));
639 // return point<DIM>(v); //MSVC does not understand this one
640 }
641 18 );
642 }
643
644 #endif
645
646 protected:
647
648 void clear_store(
649 bool keep_attributes, bool keep_memory = false
650 ) override;
651
652 void resize_store(index_t new_size) override;
653
654 void bind_point_attribute(index_t dim, bool single_precision=false);
655
656 60 void copy(const MeshVertices& rhs, bool copy_attributes=true) {
657 60 index_t dim = rhs.dimension();
658
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 60 times.
 60 if(point_fp32_.is_bound()) {
659 ✗ point_fp32_.destroy();
660 }
661
1/2
✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
 60 if(point_.is_bound()) {
662 60 point_.destroy();
663 }
664 60 MeshSubElementsStore::copy(rhs, copy_attributes);
665
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 60 times.
 60 if(rhs.single_precision()) {
666 ✗ point_fp32_.bind_if_is_defined(attributes(),"point_fp32");
667 ✗ if(!point_fp32_.is_bound()) {
668 ✗ point_fp32_.create_vector_attribute(
669 ✗ attributes(), "point_fp32", dim
670 );
671 }
672 } else {
673
2/4
✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 60 times.
✗ Branch 6 not taken.
 120 point_.bind_if_is_defined(attributes(),"point");
674
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 60 times.
 60 if(!point_.is_bound()) {
675 ✗ point_.create_vector_attribute(
676 ✗ attributes(), "point", dim
677 );
678 }
679 }
680 // Even if we do not copy the attributes, we need at least
681 // to copy the coordinates of the points !!
682
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 60 times.
 60 if(!copy_attributes) {
683 ✗ if(rhs.single_precision()) {
684 ✗ Memory::copy(
685 ✗ single_precision_point_ptr(0),
686 ✗ rhs.single_precision_point_ptr(0),
687 ✗ rhs.dimension()*rhs.nb()*sizeof(float)
688 );
689 } else {
690 ✗ Memory::copy(
691 ✗ point_ptr(0),
692 ✗ rhs.point_ptr(0),
693 ✗ rhs.dimension()*rhs.nb()*sizeof(double)
694 );
695 }
696 }
697 60 }
698
699 MeshEdges& edges_;
700 MeshFacetCornersStore& facet_corners_;
701 MeshCellCornersStore& cell_corners_;
702 Attribute<double> point_;
703 Attribute<float> point_fp32_;
704
705 friend class Mesh;
706 friend class GeogramIOHandler;
707 };
708
709 /*************************************************************************/
710
711 /**
712 * \brief The edges of a mesh.
713 * \relates Mesh
714 */
715 class GEOGRAM_API MeshEdges :
716 public MeshSubElementsStore, public MeshElements {
717 public:
718 MeshEdges(Mesh& mesh);
719 ~MeshEdges() override;
720
721 /**
722 * \brief Gets the index of an edge vertex
723 * \param[in] e index of the edge
724 * \param[in] lv local index of the vertex, in {0,1}
725 * \return the global index of vertex \p lv in edge \p e
726 */
727 114220 index_t vertex(index_t e, index_t lv) const {
728
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 114220 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 114220 geo_debug_assert(e < nb());
729
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 114220 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 114220 geo_debug_assert(lv < 2);
730 114220 return edge_vertex_[2*e+lv];
731 }
732
733 /**
734 * \brief Sets a vertex of an edge
735 * \param[in] e index of the edge
736 * \param[in] lv local index of the vertex, in {0,1}
737 * \param[in] v global index of the vertex
738 */
739 20136 void set_vertex(index_t e, index_t lv, index_t v) {
740
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 20136 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 20136 geo_debug_assert(e < nb());
741
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 20136 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 20136 geo_debug_assert(lv < 2);
742 20136 edge_vertex_[2*e+lv] = v;
743 20136 }
744
745
746 /**
747 * \brief Gets a pointer to a vertex index by corner index
748 * \param[in] c corner index (2 * edge index + 0 or 1)
749 * \return a pointer to the index of the vertex.
750 * \note Normal uses do not call this function
751 */
752 index_t* vertex_index_ptr(index_t c) {
753 geo_debug_assert(c < 2*nb());
754 return &(edge_vertex_[c]);
755 }
756
757 /**
758 * \brief Gets a pointer to a vertex index by corner index
759 * \param[in] c corner index (2 * edge index + 0 or 1)
760 * \return a pointer to the index of the vertex.
761 * \note Normal uses do not call this function
762 */
763 const index_t* vertex_index_ptr(index_t c) const {
764 geo_debug_assert(c < 2*nb());
765 return &(edge_vertex_[c]);
766 }
767
768 /**
769 * \brief Creates a new edge
770 * \return the index of the created edge
771 */
772 8226 index_t create_edge() {
773 8226 return create_sub_element();
774 }
775
776 /**
777 * \brief Creates a batch of edges
778 * \param[in] nb number of edges to create
779 * \return the index of the first created edge
780 */
781 88 index_t create_edges(index_t nb) {
782 88 return create_sub_elements(nb);
783 }
784
785 /**
786 * \brief Creates a new edge
787 * \param[in] v1 , v2 global indices of the vertices of the edge
788 * \return the index of the created edge
789 */
790 8226 index_t create_edge(index_t v1, index_t v2) {
791 8226 index_t result = create_edge();
792 8226 set_vertex(result,0,v1);
793 8226 set_vertex(result,1,v2);
794 8226 return result;
795 }
796
797 void delete_elements(
798 vector<index_t>& to_delete, bool remove_isolated_vertices=true
799 ) override;
800
801 void permute_elements(vector<index_t>& permutation) override;
802
803 void clear(
804 bool keep_attributes=true, bool keep_memory=false
805 ) override;
806
807 void pop() override;
808
809 /**
810 * \brief Swaps both extremities of an edge
811 * \param[in] e the edge
812 */
813 void flip(index_t e);
814
815 protected:
816 void clear_store(
817 bool keep_attributes, bool keep_memory = false
818 ) override;
819
820 void resize_store(index_t new_size) override;
821
822 8226 index_t create_sub_element() {
823 8226 edge_vertex_.push_back(NO_VERTEX);
824 8226 edge_vertex_.push_back(NO_VERTEX);
825 8226 return MeshSubElementsStore::create_sub_element();
826 }
827
828 88 index_t create_sub_elements(index_t nb_in) {
829 88 edge_vertex_.resize(2*(nb()+nb_in),NO_VERTEX);
830 88 return MeshSubElementsStore::create_sub_elements(nb_in);
831 }
832
833 60 void copy(const MeshEdges& rhs, bool copy_attributes=true) {
834 60 MeshSubElementsStore::copy(rhs, copy_attributes);
835 60 edge_vertex_ = rhs.edge_vertex_;
836 60 }
837
838 vector<index_t> edge_vertex_;
839 friend class Mesh;
840 friend class GeogramIOHandler;
841 };
842
843 /**************************************************************************/
844
845 /**
846 * \brief Stores the facets of a mesh (low-level store)
847 * \relates MeshFacets
848 */
849 class GEOGRAM_API MeshFacetsStore : public MeshSubElementsStore {
850 public:
851 MeshFacetsStore(Mesh& mesh);
852
853 /**
854 * \brief Gets the first element for iterating over
855 * the corners of a facet
856 * \param[in] f the facet
857 * \return the first corner of the facet
858 */
859 301512359 index_t corners_begin(index_t f) const {
860
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 301512359 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 301512359 geo_debug_assert(f < nb());
861
2/2
✓ Branch 0 taken 276809426 times.
✓ Branch 1 taken 24702933 times.
 301512359 return (is_simplicial_ ? 3*f : facet_ptr_[f]);
862 }
863
864 /**
865 * \brief Gets the upper limit for iterating over the
866 * corners of a facet
867 * \param[in] f the facet
868 * \return one position past the last corner of the facet
869 */
870 188565844 index_t corners_end(index_t f) const {
871
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 188565844 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 188565844 geo_debug_assert(f < nb());
872
2/2
✓ Branch 0 taken 176412574 times.
✓ Branch 1 taken 12153270 times.
 188565844 return (is_simplicial_ ? 3*(f+1): facet_ptr_[f+1]);
873 }
874
875 /**
876 * \brief Gets the number of corners in a facet
877 * \param[in] f the facet
878 * \return the number of corners in facet \p f
879 */
880 142626374 index_t nb_corners(index_t f) const {
881
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 142626374 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 142626374 geo_debug_assert(f < nb());
882
2/2
✓ Branch 0 taken 129670788 times.
✓ Branch 1 taken 12955586 times.
 142626374 return (is_simplicial_ ? 3 : facet_ptr_[f+1] - facet_ptr_[f]);
883 }
884
885 /**
886 * \brief Gets a corner by facet and local vertex index
887 * \param[in] f the facet
888 * \param[in] lv the local index of the vertex in facet \p f
889 * \return the \p lv%th corner of facet \p f
890 * \pre lv < nb_corners(f)
891 */
892 101321696 index_t corner(index_t f, index_t lv) const {
893
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 101321696 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 101321696 geo_debug_assert(f < nb());
894
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 101321696 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 101321696 geo_debug_assert(lv < nb_corners(f));
895 101321696 return corners_begin(f)+lv;
896 }
897
898 /**
899 * \brief Tests whether all the facets are triangles
900 * \details When all the facets are triangles, storage
901 * and access is optimized.
902 * \retval true if all the facets are triangles
903 * \retval false otherwise
904 */
905 7006151 bool are_simplices() const {
906 7006151 return is_simplicial_;
907 }
908
909 /**
910 * \brief Gets a pointer to the first element for iterating over
911 * the corners of a facet
912 * \param[in] f the facet
913 * \return a pointer to the first corner of the facet
914 */
915 index_t* corners_begin_ptr(index_t f) {
916 geo_debug_assert(!is_simplicial_);
917 geo_debug_assert(f < nb());
918 return &facet_ptr_[f];
919 }
920
921 /**
922 * \brief Gets a pointer to the first element for iterating over
923 * the corners of a facet
924 * \param[in] f the facet
925 * \return a const pointer to the first corner of the facet
926 */
927 const index_t* corners_begin_ptr(index_t f) const {
928 geo_debug_assert(!is_simplicial_);
929 geo_debug_assert(f < nb());
930 return &facet_ptr_[f];
931 }
932
933
934 protected:
935 void clear_store(
936 bool keep_attributes, bool keep_memory = false
937 ) override;
938
939 void resize_store(index_t new_size) override;
940
941 1277766 index_t create_sub_element() {
942
2/2
✓ Branch 0 taken 335166 times.
✓ Branch 1 taken 942600 times.
 1277766 if(!is_simplicial_) {
943 335166 facet_ptr_.push_back(NO_CORNER);
944 }
945 1277766 return MeshSubElementsStore::create_sub_element();
946 }
947
948 348 index_t create_sub_elements(index_t nb) {
949
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 338 times.
 348 if(!is_simplicial_) {
950
2/2
✓ Branch 0 taken 20678 times.
✓ Branch 1 taken 10 times.
 20688 for(index_t i=0; i<nb; ++i) {
951 20678 facet_ptr_.push_back(NO_CORNER);
952 }
953 }
954 348 return MeshSubElementsStore::create_sub_elements(nb);
955 }
956
957 60 void copy(const MeshFacetsStore& rhs, bool copy_attributes=true) {
958 60 MeshSubElementsStore::copy(rhs,copy_attributes);
959 60 is_simplicial_ = rhs.is_simplicial_;
960 60 facet_ptr_ = rhs.facet_ptr_;
961 60 }
962
963 protected:
964 bool is_simplicial_;
965 vector<index_t> facet_ptr_;
966 friend class Mesh;
967 friend class GeogramIOHandler;
968 };
969
970 /*************************************************************************/
971
972 /**
973 * \brief Stores the facet corners of a mesh (low-level store)
974 * \relates MeshFacets
975 */
976 class GEOGRAM_API MeshFacetCornersStore : public MeshSubElementsStore {
977 public:
978 MeshFacetCornersStore(Mesh& mesh);
979
980 /**
981 * \brief Gets the vertex that a corner is incident to
982 * \param[in] c the corner
983 * \return the vertex that corner \p c is incident to
984 */
985 344376338 index_t vertex(index_t c) const {
986
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 344376338 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 344376338 geo_debug_assert(c < nb());
987 344376338 return corner_vertex_[c];
988 }
989
990 /**
991 * \brief Gets the facet that a corner is adjacent to
992 * \param[in] c the corner
993 * \return the facet that corner \p is adjacent to or
994 * NO_FACET if \p c is on the border
995 */
996 57086569 index_t adjacent_facet(index_t c) const {
997
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 57086569 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 57086569 geo_debug_assert(c < nb());
998 57086569 return corner_adjacent_facet_[c];
999 }
1000
1001 /**
1002 * \brief Gets a pointer to the the facet index
1003 * that a corner is adjacent to
1004 * \param[in] c the corner
1005 * \return a pointer to the the facet index
1006 * that corner \p is adjacent to.
1007 */
1008 const index_t* adjacent_facet_ptr(index_t c) const {
1009 geo_debug_assert(c < nb());
1010 return &corner_adjacent_facet_[c];
1011 }
1012
1013
1014 /**
1015 * \brief Gets a pointer to the the facet index
1016 * that a corner is adjacent to
1017 * \param[in] c the corner
1018 * \return a pointer to the the facet index
1019 * that corner \p is adjacent to.
1020 */
1021 index_t* adjacent_facet_ptr(index_t c) {
1022 geo_debug_assert(c < nb());
1023 return &corner_adjacent_facet_[c];
1024 }
1025
1026 /**
1027 * \brief Sets the vertex that a corner is incident to
1028 * \param[in] c the corner
1029 * \param[in] v the vertex that corner \p c is incident to
1030 * \pre v < mesh.vertices.nb()
1031 */
1032 10912063 void set_vertex(index_t c, index_t v) {
1033
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 10912063 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 10912063 geo_debug_assert(c < nb());
1034
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 10912063 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 10912063 geo_debug_assert(v < vertices_.nb());
1035 10912063 corner_vertex_[c] = v;
1036 10912063 }
1037
1038 /**
1039 * \brief Sets the vertex that a corner is incident to
1040 * \details Does not check whether \p v is a valid vertex
1041 * index. This function is useful for some algorithms that
1042 * need to create/update the facets before creating the
1043 * vertices.
1044 * \param[in] c the corner
1045 * \param[in] v the vertex that corner \p c is incident to
1046 * \note Normal uses do not call this function
1047 */
1048 2185176 void set_vertex_no_check(index_t c, index_t v) {
1049
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 2185176 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 2185176 geo_debug_assert(c < nb());
1050 2185176 corner_vertex_[c] = v;
1051 2185176 }
1052
1053 /**
1054 * \brief Sets the facet that a corner is adjacent to
1055 * \param[in] c the corner
1056 * \param[in] f the facet that corner \p is adjacent to or
1057 * NO_FACET if \p c is on the border
1058 */
1059 23665569 void set_adjacent_facet(index_t c, index_t f) {
1060
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 23665569 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 23665569 geo_debug_assert(c < nb());
1061
4/10
✓ Branch 0 taken 11479840 times.
✓ Branch 1 taken 12185729 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 11479840 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 23665569 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 23665569 geo_debug_assert(f == NO_FACET || f < facets_.nb());
1062 23665569 corner_adjacent_facet_[c] = f;
1063 23665569 }
1064
1065 /**
1066 * \brief Gets a pointer to the vertex that a corner is incident to
1067 * \param[in] c the corner
1068 * \return a pointer to the index of the vertex that this corner
1069 * is incident to
1070 * \note Normal uses do not call this function
1071 */
1072 index_t* vertex_index_ptr(index_t c) {
1073 geo_debug_assert(c < nb());
1074 return &(corner_vertex_[c]);
1075 }
1076
1077 /**
1078 * \brief Gets a pointer to the vertex that a corner is incident to
1079 * \param[in] c the corner
1080 * \return a pointer to the index of the vertex that this corner
1081 * is incident to
1082 * \note Normal uses do not call this function
1083 */
1084 const index_t* vertex_index_ptr(index_t c) const {
1085 geo_debug_assert(c < nb());
1086 return &(corner_vertex_[c]);
1087 }
1088
1089 /**
1090 * \brief Gets the vertex indices of all corners as a single vector.
1091 * \details On a triangulated mesh, this is the triangle list
1092 * (3 vertex indices per facet). It is forbidden to change the
1093 * size of the returned vector.
1094 * \return a const reference to the corner-to-vertex table.
1095 */
1096 const vector<index_t>& vertex_indices() const {
1097 return corner_vertex_;
1098 }
1099
1100 /**
1101 * \brief Gets the point associated with a corner
1102 * \param[in] c the corner
1103 * \return a reference to the DIM-d point associated with the corner
1104 * \pre vertices.dimension() >= DIM
1105 */
1106 ✗ template <index_t DIM=3> vecng<DIM,double>& point(index_t c) {
1107 ✗ geo_debug_assert(c < nb());
1108 ✗ return vertices_.point<DIM>(vertex(c));
1109 }
1110
1111 /**
1112 * \brief Gets the point associated with a corner
1113 * \param[in] c the corner
1114 * \return const a reference to the DIM-d point associated with
1115 * the corner
1116 * \pre vertices.dimension() >= DIM
1117 */
1118 102006 template <index_t DIM=3> const vecng<DIM,double>& point(
1119 index_t c
1120 ) const {
1121
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 102006 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 102006 geo_debug_assert(c < nb());
1122 102006 return vertices_.point(vertex(c));
1123 }
1124
1125 protected:
1126 void clear_store(
1127 bool keep_attributes, bool keep_memory = false
1128 ) override;
1129
1130 void resize_store(index_t new_size) override;
1131
1132 4176017 index_t create_sub_element(index_t v, index_t f = NO_FACET) {
1133 4176017 corner_vertex_.push_back(v);
1134 4176017 corner_adjacent_facet_.push_back(f);
1135 4176017 return MeshSubElementsStore::create_sub_element();
1136 }
1137
1138 348 index_t create_sub_elements(index_t nb) {
1139
2/2
✓ Branch 0 taken 2013704 times.
✓ Branch 1 taken 348 times.
 2014052 for(index_t i=0; i<nb; ++i) {
1140 2013704 corner_vertex_.push_back(NO_VERTEX);
1141 }
1142
2/2
✓ Branch 0 taken 2013704 times.
✓ Branch 1 taken 348 times.
 2014052 for(index_t i=0; i<nb; ++i) {
1143 2013704 corner_adjacent_facet_.push_back(NO_FACET);
1144 }
1145 348 return MeshSubElementsStore::create_sub_elements(nb);
1146 }
1147
1148 60 void copy(
1149 const MeshFacetCornersStore& rhs, bool copy_attributes=true
1150 ) {
1151 60 MeshSubElementsStore::copy(rhs, copy_attributes);
1152 60 corner_vertex_ = rhs.corner_vertex_;
1153 60 corner_adjacent_facet_ = rhs.corner_adjacent_facet_;
1154 60 }
1155
1156 protected:
1157 MeshVertices& vertices_;
1158 MeshFacetsStore& facets_;
1159 vector<index_t> corner_vertex_;
1160 vector<index_t> corner_adjacent_facet_;
1161
1162 friend class MeshFacets;
1163 friend class Mesh;
1164 friend class GeogramIOHandler;
1165 };
1166
1167 /*************************************************************************/
1168
1169 /**
1170 * \brief The facets of a mesh
1171 * \relates Mesh
1172 */
1173 class GEOGRAM_API MeshFacets : public MeshFacetsStore, public MeshElements {
1174 public:
1175
1176 /**
1177 * \brief MeshFacets constructor
1178 * \param[in] mesh a reference to the mesh
1179 * this MeshFacets is attached to
1180 */
1181 MeshFacets(Mesh& mesh);
1182
1183 /**
1184 * \brief Gets the number of vertices of a facet
1185 * \param[in] f the facet
1186 * \return the number of vertices of facet \p f
1187 */
1188 41304678 index_t nb_vertices(index_t f) const {
1189 41304678 return nb_corners(f);
1190 }
1191
1192 /**
1193 * \brief Gets a vertex by facet and local vertex index
1194 * \param[in] f the facet
1195 * \param[in] lv the local vertex index in \p f
1196 * \return the \p lv%th vertex of facet \p f
1197 * \pre lv < nb_vertices(f)
1198 */
1199 90659713 index_t vertex(index_t f, index_t lv) const {
1200 90659713 return facet_corners_.vertex(corner(f,lv));
1201 }
1202
1203
1204 /**
1205 * \brief Gets a point by facet and local vertex index
1206 * \param[in] f the facet
1207 * \param[in] lv the local vertex index in \p f
1208 * \return a const reference to the point corresponding to
1209 * the \p lv%th vertex of facet \p f
1210 * \pre lv < nb_vertices(f)
1211 */
1212 7681336 template <index_t DIM=3> const vecng<DIM,double>& point(
1213 index_t f, index_t lv
1214 ) const {
1215 7681336 return vertices_.point<DIM>(vertex(f,lv));
1216 }
1217
1218 /**
1219 * \brief Gets a point by facet and local vertex index
1220 * \param[in] f the facet
1221 * \param[in] lv the local vertex index in \p f
1222 * \return a reference to the point corresponding to
1223 * the \p lv%th vertex of facet \p f
1224 * \pre lv < nb_vertices(f)
1225 */
1226 1525938 template <index_t DIM=3> vecng<DIM,double>& point(
1227 index_t f, index_t lv
1228 ) {
1229 1525938 return vertices_.point<DIM>(vertex(f,lv));
1230 }
1231
1232 /**
1233 * \brief Sets a vertex by facet and local vertex index
1234 * \param[in] f the facet
1235 * \param[in] lv the local vertex index in \p f
1236 * \param[in] v specifies the \p lv%th vertex of facet \p f
1237 * \pre lv < nb_vertices(f)
1238 */
1239 4393696 void set_vertex(index_t f, index_t lv, index_t v) {
1240 4393696 facet_corners_.set_vertex(corner(f,lv),v);
1241 4393696 }
1242
1243 /**
1244 * \brief Gets the local index of a vertex in a facet.
1245 * \param[in] f a facet
1246 * \param[in] v a vertex
1247 * \return lv such that vertex(f,lv) == v or NO_VERTE if f is
1248 * not incident to v
1249 */
1250 ✗ index_t find_vertex(index_t f, index_t v) const {
1251 ✗ for(index_t lv=0; lv<nb_vertices(f); ++lv) {
1252 ✗ if(vertex(f,lv) == v) {
1253 ✗ return lv;
1254 }
1255 }
1256 ✗ return NO_VERTEX;
1257 }
1258
1259 /**
1260 * \brief finds a common vertex shared by two facets
1261 * \param[in] f1 , f2 the two facets
1262 * \return the local index in \p f1 of a vertex present in \p f2,
1263 * or NO_VERTEX if there is no such vertex.
1264 */
1265 ✗ index_t find_common_vertex(index_t f1, index_t f2) const {
1266 ✗ for(index_t lv=0; lv<nb_vertices(f1); ++lv) {
1267 ✗ index_t v = vertex(f1,lv);
1268 ✗ if(find_vertex(f2,v) != NO_VERTEX) {
1269 ✗ return lv;
1270 }
1271 }
1272 ✗ return NO_VERTEX;
1273 }
1274
1275 /**
1276 * \brief Gets an adjacent facet by facet and local edge index
1277 * \param[in] f the facet
1278 * \param[in] le the local index of an edge in facet \p f
1279 * \return the facet incident to \p f along edge \p le or
1280 * NO_FACET if \p le is on the border
1281 */
1282 4205277 index_t adjacent(index_t f, index_t le) const {
1283 4205277 return facet_corners_.adjacent_facet(corner(f,le));
1284 }
1285
1286 /**
1287 * \brief Gets the local index of a facet adjacent to another one.
1288 * \param[in] f a facet
1289 * \param[in] f2 another facet
1290 * \return le such that adjacent(f,le) == f2 or NO_INDEX if f and f2
1291 * are not adjacent.
1292 */
1293 267235 index_t find_adjacent(index_t f, index_t f2) const {
1294
1/2
✓ Branch 1 taken 536403 times.
✗ Branch 2 not taken.
 536403 for(index_t le=0; le<nb_vertices(f); ++le) {
1295
2/2
✓ Branch 1 taken 267235 times.
✓ Branch 2 taken 269168 times.
 536403 if(adjacent(f,le) == f2) {
1296 267235 return le;
1297 }
1298 }
1299 ✗ return NO_INDEX;
1300 }
1301
1302 /**
1303 * \brief Sets an adjacent facet by facet and local edge index
1304 * \param[in] f the facet
1305 * \param[in] le the local index of an edge in facet \p f
1306 * \param[in] f2 specifies the facet incident to \p f along edge
1307 * \p le or NO_FACET if \p le is on the border
1308 */
1309 356940 void set_adjacent(index_t f, index_t le, index_t f2) {
1310 356940 facet_corners_.set_adjacent_facet(corner(f,le),f2);
1311 356940 }
1312
1313 /**
1314 * \brief Gets the successor of a corner around a facet
1315 * \param[in] f the facet
1316 * \param[in] c the corner
1317 * \return the successor of corner \p c around facet \p f
1318 * \pre c >= corners_begin(f) && c < corners_end(f)
1319 */
1320 18387862 index_t next_corner_around_facet(index_t f, index_t c) const {
1321
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 18387862 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 18387862 geo_debug_assert(f < nb());
1322
3/10
✓ Branch 1 taken 18387862 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 18387862 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 18387862 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 18387862 geo_debug_assert(c >= corners_begin(f) && c < corners_end(f));
1323
2/2
✓ Branch 1 taken 5592798 times.
✓ Branch 2 taken 12795064 times.
 18387862 return c + 1 == corners_end(f) ? corners_begin(f) : c + 1;
1324 }
1325
1326 /**
1327 * \brief Gets the predecessor of a corner around a facet
1328 * \param[in] f the facet
1329 * \param[in] c the corner
1330 * \return the predecessor of corner \p c around facet \p f
1331 * \pre c >= corners_begin(f) && c < corners_end(f)
1332 */
1333 41056732 index_t prev_corner_around_facet(index_t f, index_t c) const {
1334
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 41056732 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 41056732 geo_debug_assert(f < nb());
1335
3/10
✓ Branch 1 taken 41056732 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 41056732 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 41056732 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 41056732 geo_debug_assert(c >= corners_begin(f) && c < corners_end(f));
1336
2/2
✓ Branch 1 taken 7535745 times.
✓ Branch 2 taken 33520987 times.
 41056732 return c == corners_begin(f) ? corners_end(f) - 1 : c - 1;
1337 }
1338
1339 /**
1340 * \brief Finds an edge by vertex indices
1341 * \param[in] f a facet
1342 * \param[in] v1 , v2 two vertex indices
1343 * \return the edge le such that vertex(f,le) = v1 and
1344 * vertex(f, (le+1)%nb_vertices(f)) == v2
1345 */
1346 370476 index_t find_edge(index_t f, index_t v1, index_t v2) const {
1347
1/2
✓ Branch 2 taken 740952 times.
✗ Branch 3 not taken.
 740952 for(index_t c1 = corners_begin(f); c1 != corners_end(f); ++c1) {
1348 740952 index_t c2 = next_corner_around_facet(f,c1);
1349 740952 if(
1350
4/4
✓ Branch 1 taken 370476 times.
✓ Branch 2 taken 370476 times.
✓ Branch 3 taken 370476 times.
✓ Branch 4 taken 370476 times.
 1111428 facet_corners_.vertex(c1) == v1 &&
1351
1/2
✓ Branch 1 taken 370476 times.
✗ Branch 2 not taken.
 370476 facet_corners_.vertex(c2) == v2
1352 ) {
1353 370476 return c1 - corners_begin(f);
1354 }
1355 }
1356 ✗ return NO_INDEX;
1357 }
1358
1359 void delete_elements(
1360 vector<index_t>& to_delete,
1361 bool remove_isolated_vertices=true
1362 ) override;
1363
1364 void permute_elements(vector<index_t>& permutation) override;
1365
1366 void clear(
1367 bool keep_attributes=true, bool keep_memory=false
1368 ) override;
1369
1370 /**
1371 * \brief Creates a contiguous chunk of facets
1372 * \param[in] nb_facets number of facets to create
1373 * \param[in] nb_vertices_per_polygon number of vertices
1374 * in each facet
1375 * \return the index of the first facet
1376 */
1377 348 index_t create_facets(
1378 index_t nb_facets, index_t nb_vertices_per_polygon
1379 ) {
1380
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 338 times.
 348 if(nb_vertices_per_polygon != 3) {
1381 10 is_not_simplicial();
1382 }
1383
1384 348 index_t first_facet = nb();
1385 348 index_t co = facet_corners_.nb();
1386 348 facet_corners_.create_sub_elements(
1387 nb_facets*nb_vertices_per_polygon
1388 );
1389 348 index_t result = create_sub_elements(nb_facets);
1390
1391
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 338 times.
 348 if(!is_simplicial_) {
1392
2/2
✓ Branch 0 taken 20688 times.
✓ Branch 1 taken 10 times.
 20698 for(index_t f=first_facet; f<=first_facet+nb_facets; ++f) {
1393 20688 facet_ptr_[f] = co;
1394 20688 co += nb_vertices_per_polygon;
1395 }
1396
1/6
✗ Branch 2 not taken.
✓ Branch 3 taken 10 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 10 geo_debug_assert(facet_ptr_.size() == nb()+1);
1397
1/6
✗ Branch 3 not taken.
✓ Branch 4 taken 10 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 10 geo_debug_assert(facet_ptr_[nb()] == facet_corners_.nb());
1398 }
1399 348 return result;
1400 }
1401
1402 /**
1403 * \brief Reserves space for new facets
1404 * \param[in] nb_to_reserve the number of facets to reserve
1405 * \details Does not change size
1406 */
1407 57 void reserve(index_t nb_to_reserve) {
1408 57 facet_corners_.reserve_store(nb_to_reserve*3);
1409 57 this->reserve_store(nb_to_reserve);
1410 57 }
1411
1412 /**
1413 * \brief Creates a contiguous chunk of triangles
1414 * \param[in] nb_triangles number of triangles to create
1415 * \return the index of the first triangle
1416 */
1417 338 index_t create_triangles(index_t nb_triangles) {
1418 338 return create_facets(nb_triangles, 3);
1419 }
1420
1421 /**
1422 * \brief Creates a contiguous chunk of quads
1423 * \param[in] nb_quads number of quads to create
1424 * \return the index of the first quad
1425 */
1426 10 index_t create_quads(index_t nb_quads) {
1427 10 return create_facets(nb_quads, 4);
1428 }
1429
1430 /**
1431 * \brief Creates a triangle
1432 * \param[in] v1 , v2 , v3 the vertices of the triangle
1433 * \return the index of the created triangle
1434 */
1435 659930 index_t create_triangle(index_t v1, index_t v2, index_t v3) {
1436
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 659930 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 659930 geo_debug_assert(v1 != v2);
1437
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 659930 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 659930 geo_debug_assert(v2 != v3);
1438
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 659930 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 659930 geo_debug_assert(v3 != v1);
1439 659930 facet_corners_.create_sub_element(v1);
1440 659930 facet_corners_.create_sub_element(v2);
1441 659930 facet_corners_.create_sub_element(v3);
1442 659930 index_t result = create_sub_element();
1443
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 659930 times.
 659930 if(!is_simplicial_) {
1444 ✗ facet_ptr_[result+1] = facet_corners_.nb();
1445 ✗ geo_debug_assert(facet_ptr_.size() == nb()+1);
1446 ✗ geo_debug_assert(facet_ptr_[nb()] == facet_corners_.nb());
1447 }
1448 659930 return result;
1449 }
1450
1451 /**
1452 * \brief Creates a quad
1453 * \param[in] v1 , v2 , v3 , v4 the vertices of the quad
1454 * \return the index of the created quad
1455 */
1456 108 index_t create_quad(index_t v1, index_t v2, index_t v3, index_t v4) {
1457 108 is_not_simplicial();
1458 108 facet_corners_.create_sub_element(v1);
1459 108 facet_corners_.create_sub_element(v2);
1460 108 facet_corners_.create_sub_element(v3);
1461 108 facet_corners_.create_sub_element(v4);
1462 108 index_t result = create_sub_element();
1463 108 facet_ptr_[result+1] = facet_corners_.nb();
1464
1/6
✗ Branch 2 not taken.
✓ Branch 3 taken 108 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 108 geo_debug_assert(facet_ptr_.size() == nb()+1);
1465
1/6
✗ Branch 3 not taken.
✓ Branch 4 taken 108 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 108 geo_debug_assert(facet_ptr_[nb()] == facet_corners_.nb());
1466 108 return result;
1467 }
1468
1469 /**
1470 * \brief Creates a polygonal facet
1471 * \param[in] nb_vertices number of vertices of the facet
1472 * \return the index of the created facet
1473 */
1474 617728 index_t create_polygon(index_t nb_vertices) {
1475
2/2
✓ Branch 0 taken 248598 times.
✓ Branch 1 taken 369130 times.
 617728 if(nb_vertices != 3) {
1476 248598 is_not_simplicial();
1477 }
1478
2/2
✓ Branch 0 taken 2195795 times.
✓ Branch 1 taken 617728 times.
 2813523 for(index_t i=0; i<nb_vertices; ++i) {
1479 2195795 facet_corners_.create_sub_element(NO_VERTEX);
1480 }
1481 617728 index_t result = create_sub_element();
1482
2/2
✓ Branch 0 taken 335058 times.
✓ Branch 1 taken 282670 times.
 617728 if(!is_simplicial_) {
1483 335058 facet_ptr_[result+1] = facet_corners_.nb();
1484
1/6
✗ Branch 2 not taken.
✓ Branch 3 taken 335058 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 335058 geo_debug_assert(facet_ptr_.size() == nb()+1);
1485
1/6
✗ Branch 3 not taken.
✓ Branch 4 taken 335058 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 335058 geo_debug_assert(facet_ptr_[nb()] == facet_corners_.nb());
1486 }
1487 617728 return result;
1488 }
1489
1490 /**
1491 * \brief Creates a polygonal facet
1492 * \param[in] nb_vertices number of vertices of the facet
1493 * \param[in] vertices a const pointer to the \p nb_vertices vertices
1494 * \return the index of the created facet
1495 */
1496 index_t create_polygon(index_t nb_vertices, const index_t* vertices) {
1497 if(nb_vertices != 3) {
1498 is_not_simplicial();
1499 }
1500 for(index_t i=0; i<nb_vertices; ++i) {
1501 facet_corners_.create_sub_element(vertices[i]);
1502 }
1503 index_t result = create_sub_element();
1504 if(!is_simplicial_) {
1505 facet_ptr_[result+1] = facet_corners_.nb();
1506 geo_debug_assert(facet_ptr_.size() == nb()+1);
1507 geo_debug_assert(facet_ptr_[nb()] == facet_corners_.nb());
1508 }
1509 return result;
1510 }
1511
1512 /**
1513 * \brief Creates a polygonal facet
1514 * \param[in] vertices a const reference to a vector that
1515 * contains the vertices
1516 * \return the index of the created facet
1517 */
1518 index_t create_polygon(const vector<index_t>& vertices) {
1519 return create_polygon(vertices.size(), vertices.data());
1520 }
1521
1522 /**
1523 * \brief Connects the facets
1524 * \details Finds the adjacent_facet() links based on vertex indices
1525 */
1526 void connect();
1527
1528 /**
1529 * \brief Connects a contiguous sequence of facets
1530 * \details Finds the adjacent_facet() links based on vertex indices
1531 * \param[in] facets_begin first facet to connect
1532 * \param[in] facets_end one position past the last facet to connect
1533 */
1534 void connect(index_t facets_begin, index_t facets_end);
1535
1536 /**
1537 * \brief Triangulates the facets
1538 * \note Attributes are zeroed
1539 */
1540 void triangulate();
1541
1542 /**
1543 * \brief Flips a facet
1544 * \details The order of the corners is reversed
1545 * \param[in] f the facet to be flipped
1546 */
1547 void flip(index_t f);
1548
1549 /**
1550 * \brief Replaces the edges of this mesh
1551 * with the borders of the surfacic part.
1552 */
1553 void compute_borders();
1554
1555 /**
1556 * \brief Copies a triangle mesh into this Mesh.
1557 * \details Facet adjacence are not computed.
1558 * Facet and corner attributes are zeroed.
1559 * \param[in] dim dimension of the vertices
1560 * \param[in] vertices coordinates of the vertices
1561 * \param[in] triangles facet to vertex links
1562 * \param[in] steal_args if set, vertices and triangles
1563 * are 'stolen' from the arguments
1564 * (using vector::swap).
1565 */
1566 void assign_triangle_mesh(
1567 coord_index_t dim,
1568 vector<double>& vertices,
1569 vector<index_t>& triangles,
1570 bool steal_args
1571 );
1572
1573 /*
1574 * \brief Copies a triangle mesh into this Mesh.
1575 * \details Facet adjacence are not computed.
1576 * Facet and corner attributes are zeroed.
1577 * \param[in] triangles facet to vertex links
1578 * \param[in] steal_args if set, vertices and triangles
1579 * are 'stolen' from the arguments
1580 * (using vector::swap).
1581 */
1582 void assign_triangle_mesh(
1583 vector<index_t>& triangles,
1584 bool steal_args
1585 );
1586
1587 void pop() override;
1588
1589 /**
1590 * \brief Gets the corners of a facet.
1591 * \param[in] f the index of the facet.
1592 * \return a range with all the corners of the facet.
1593 */
1594 30334065 index_range corners(index_t f) const {
1595
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 30334065 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 30334065 geo_debug_assert(f < nb());
1596 30334065 return index_range(
1597
1/2
✓ Branch 1 taken 30334065 times.
✗ Branch 2 not taken.
 30334065 index_as_iterator(corners_begin(f)),
1598
1/2
✓ Branch 1 taken 30334065 times.
✗ Branch 2 not taken.
 30334065 index_as_iterator(corners_end(f))
1599 30334065 );
1600 }
1601
1602 #ifndef MESH_NO_SYNTAXIC_SUGAR
1603
1604 /**
1605 * \brief Gets the vertices of a facet.
1606 * \param[in] f the index of the facet.
1607 * \return a range with all the vertices indices of the facet.
1608 */
1609 175 auto vertices(index_t f) const {
1610
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 175 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 175 geo_debug_assert(f < nb());
1611
1/2
✓ Branch 1 taken 175 times.
✗ Branch 2 not taken.
 175 return transform_range(
1612
1/2
✓ Branch 1 taken 175 times.
✗ Branch 2 not taken.
 350 corners(f), [this](index_t c)->index_t {
1613 525 return facet_corners_.vertex(c);
1614 }
1615 350 );
1616 }
1617
1618 /**
1619 * \brief Gets the facets adjacent to a given facet.
1620 * \param[in] f the index of the facet.
1621 * \return a range with all the indices of the facets adjacent to this
1622 * facet. It will output exactly one item per edge of the facet. Edges
1623 * on the border will output NO_INDEX (no adjacent facet).
1624 */
1625 2591187 auto adjacent(index_t f) const {
1626
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 2591187 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 2591187 geo_debug_assert(f < nb());
1627
1/2
✓ Branch 1 taken 2591187 times.
✗ Branch 2 not taken.
 2591187 return transform_range(
1628
1/2
✓ Branch 1 taken 2591187 times.
✗ Branch 2 not taken.
 5182374 corners(f), [this](index_t c)->index_t {
1629 8018604 return facet_corners_.adjacent_facet(c);
1630 }
1631 5182374 );
1632 }
1633
1634 /**
1635 * \brief Gets the points associated with the vertices of a facet.
1636 * \param[in] f the index of the facet.
1637 * \return a range with the points that correspond to the vertices of
1638 * the facet, returned as const references.
1639 */
1640 template <index_t DIM=3> auto points(index_t f) const {
1641 geo_debug_assert(f < nb());
1642 typedef vecng<DIM,double> vecn;
1643 return transform_range_ref(
1644 corners(f), [this](index_t c)->const vecn& {
1645 index_t v = facet_corners_.vertex(c);
1646 return vertices_.point<DIM>(v);
1647 }
1648 );
1649 }
1650
1651 /**
1652 * \brief Gets the points associated with the vertices of a facet.
1653 * \param[in] f the index of the facet.
1654 * \return a range with the points that correspond to the vertices of
1655 * the facet, returned as modifiable references.
1656 */
1657 385468 template <index_t DIM=3> auto points(index_t f) {
1658
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 192734 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 385468 geo_debug_assert(f < nb());
1659 typedef vecng<DIM,double> vecn;
1660
1/2
✓ Branch 1 taken 192734 times.
✗ Branch 2 not taken.
 385468 return transform_range_ref(
1661
1/2
✓ Branch 1 taken 192734 times.
✗ Branch 2 not taken.
 1349316 corners(f), [this](index_t c)->vecn& {
1662 578380 index_t v = facet_corners_.vertex(c);
1663 578380 return vertices_.point<DIM>(v);
1664 }
1665 770936 );
1666 }
1667
1668 /**
1669 * \brief Decomposes a facet into triangles
1670 * \param[in] f the index of the facet.
1671 * \return a range with a decomposition of the facet into triangles,
1672 * returned as std::tuple<index_t, index_t, index_t>
1673 */
1674 3684340 auto triangles(index_t f) const {
1675
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 3684340 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 3684340 geo_debug_assert(f < nb());
1676 3684340 index_t v0 = facet_corners_.vertex(corners_begin(f));
1677
1/2
✓ Branch 1 taken 3684340 times.
✗ Branch 2 not taken.
 3684340 return transform_range(
1678 3684340 index_range(
1679
1/2
✓ Branch 1 taken 3684340 times.
✗ Branch 2 not taken.
 3684340 index_as_iterator(corners_begin(f)+1),
1680
1/2
✓ Branch 1 taken 3684340 times.
✗ Branch 2 not taken.
 3684340 index_as_iterator(corners_end(f)-1)
1681 7368680 ),
1682 3700430 [this,v0](index_t c)->std::tuple<index_t, index_t, index_t> {
1683 return std::make_tuple(
1684 3700430 v0,
1685 ✗ facet_corners_.vertex(c),
1686
2/4
✓ Branch 1 taken 3700430 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3700430 times.
✗ Branch 5 not taken.
 7400860 facet_corners_.vertex(c+1)
1687
1/2
✓ Branch 1 taken 3700430 times.
✗ Branch 2 not taken.
 7400860 );
1688 }
1689 7368680 );
1690 }
1691
1692 /**
1693 * \brief Decomposes a facet into triangles
1694 * \param[in] f the index of the facet.
1695 * \return a range with a decomposition of the facet into triangles,
1696 * returned as std::tuple<const vecn&, const vecn&, const vecn&>
1697 * where vecn can be vec2, vec3, vec4 ... depending on the DIM
1698 * template argument
1699 */
1700 3144068 template <index_t DIM=3> auto triangle_points(index_t f) const {
1701
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 3144068 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 3144068 geo_debug_assert(f < nb());
1702 typedef vecng<DIM,double> vecn;
1703
1/2
✓ Branch 1 taken 3144068 times.
✗ Branch 2 not taken.
 3144068 return transform_range(
1704
1/2
✓ Branch 1 taken 3144068 times.
✗ Branch 2 not taken.
 6288136 triangles(f),
1705 3159428 [this](std::tuple<index_t, index_t, index_t> T)
1706 ->std::tuple<const vecn&, const vecn&, const vecn&> {
1707 return std::tuple<const vecn&, const vecn&, const vecn&>(
1708 3159428 vertices_.point<DIM>(std::get<0>(T)),
1709 3159428 vertices_.point<DIM>(std::get<1>(T)),
1710 3159428 vertices_.point<DIM>(std::get<2>(T))
1711 3159428 );
1712 }
1713 6288136 );
1714 }
1715
1716 /**
1717 * \brief Decomposes a facet into triangles
1718 * \param[in] f the index of the facet.
1719 * \return a range with a decomposition of the facet into triangles,
1720 * returned as std::tuple<vecn&, vecn&, vecn&>
1721 * where vecn can be vec2, vec3, vec4 ... depending on the DIM
1722 * template argument
1723 */
1724 48564 template <index_t DIM=3> auto triangle_points(index_t f) {
1725
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 48564 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 48564 geo_debug_assert(f < nb());
1726 typedef vecng<DIM,double> vecn;
1727
1/2
✓ Branch 1 taken 48564 times.
✗ Branch 2 not taken.
 48564 return transform_range(
1728
1/2
✓ Branch 1 taken 48564 times.
✗ Branch 2 not taken.
 97128 triangles(f),
1729 48564 [this](std::tuple<index_t, index_t, index_t> T)
1730 ->std::tuple<vecn&, vecn&, vecn&> {
1731 return std::tuple<vecn&, vecn&, vecn&>(
1732 48564 vertices_.point(std::get<0>(T)),
1733 48564 vertices_.point(std::get<1>(T)),
1734 48564 vertices_.point(std::get<2>(T))
1735 48564 );
1736 }
1737 97128 );
1738 }
1739
1740 #endif
1741
1742 protected:
1743
1744 /**
1745 * \brief Indicates that the stored elements are only triangles.
1746 */
1747 578 void is_simplicial() {
1748
2/2
✓ Branch 0 taken 43 times.
✓ Branch 1 taken 535 times.
 578 if(!is_simplicial_) {
1749 43 is_simplicial_ = true;
1750 43 facet_ptr_.resize(1);
1751 43 facet_ptr_[0] = 0;
1752 }
1753 578 }
1754
1755 /**
1756 * \brief Indicates that the stored elements are no
1757 * longer only triangles.
1758 * \details Creates the facet pointers for the pre-existing
1759 * triangles if any.
1760 */
1761 248716 void is_not_simplicial() {
1762
2/2
✓ Branch 0 taken 142 times.
✓ Branch 1 taken 248574 times.
 248716 if(is_simplicial_) {
1763 142 is_simplicial_ = false;
1764 142 facet_ptr_.resize(nb()+1);
1765
2/2
✓ Branch 1 taken 185 times.
✓ Branch 2 taken 142 times.
 327 for(index_t f=0; f<facet_ptr_.size(); ++f) {
1766 185 facet_ptr_[f] = 3*f;
1767 }
1768 }
1769 248716 }
1770
1771 protected:
1772 MeshVertices& vertices_;
1773 MeshFacetCornersStore& facet_corners_;
1774 friend class Mesh;
1775 friend class GeogramIOHandler;
1776 friend void GEOGRAM_API tessellate_facets(
1777 Mesh& M, index_t max_nb_vertices
1778 );
1779 };
1780
1781 /*************************************************************************/
1782
1783 enum MeshCellType {
1784 MESH_TET = 0,
1785 MESH_HEX = 1,
1786 MESH_PRISM = 2,
1787 MESH_PYRAMID = 3,
1788 MESH_CONNECTOR = 4,
1789 MESH_NB_CELL_TYPES = 5
1790 };
1791
1792 /**
1793 * \brief Lookup tables that describe the combinatorics
1794 * of each cell type.
1795 * \relates MeshCells
1796 */
1797 struct CellDescriptor {
1798 /** Number of vertices */
1799 index_t nb_vertices;
1800
1801 /** Number of facets */
1802 index_t nb_facets;
1803
1804 /** Number of vertices in each facet */
1805 index_t nb_vertices_in_facet[6];
1806
1807 /**
1808 * Cell vertex index by (facet index,facet vertex index).
1809 */
1810 index_t facet_vertex[6][4];
1811
1812 /**
1813 * Number of edges */
1814 index_t nb_edges;
1815
1816 /**
1817 * Cell vertex index by (edge index, edge vertex index).
1818 */
1819 index_t edge_vertex[12][2];
1820
1821 /**
1822 * Cell facet index by (edge index, adjacent facet index).
1823 */
1824 index_t edge_adjacent_facet[12][2];
1825 };
1826
1827
1828 /**
1829 * \brief Gathers declarations of global cell descriptors.
1830 * \details Cannot be declared as static variables in
1831 * MeshCellsStore, since visual C++ does not allows
1832 * exporting class static variables from a DLL.
1833 */
1834 namespace MeshCellDescriptors {
1835 /**
1836 * \brief Maps a cell type to the associated cell descriptor.
1837 */
1838 GEOGRAM_API extern CellDescriptor*
1839 cell_type_to_cell_descriptor[GEO::MESH_NB_CELL_TYPES];
1840
1841 GEOGRAM_API extern CellDescriptor tet_descriptor;
1842 GEOGRAM_API extern CellDescriptor hex_descriptor;
1843 GEOGRAM_API extern CellDescriptor prism_descriptor;
1844 GEOGRAM_API extern CellDescriptor pyramid_descriptor;
1845 GEOGRAM_API extern CellDescriptor connector_descriptor;
1846 }
1847
1848 /**
1849 * \brief Stores the cells of a mesh (low-level store)
1850 * \relates MeshCells
1851 */
1852 class GEOGRAM_API MeshCellsStore : public MeshSubElementsStore {
1853 public:
1854 MeshCellsStore(Mesh& mesh);
1855
1856 /**
1857 * \brief Tests whether all the cells are tetrahedra
1858 * \details Storage and access are optimized when all the
1859 * cells are tetrahedra
1860 * \retval true if all the cells are tetrahedra
1861 * \retval false otherwise
1862 */
1863 15 bool are_simplices() const {
1864 15 return is_simplicial_;
1865 }
1866
1867 /**
1868 * \brief Gets the type of a cell
1869 * \param[in] c the cell, in 0..nb()-1
1870 * \return one of
1871 * MESH_TET, MESH_HEX, MESH_PRISM, MESH_PYRAMID, MESH_CONNECTOR.
1872 */
1873 3116 MeshCellType type(index_t c) const {
1874
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 3116 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 3116 geo_debug_assert(c < nb());
1875
1/2
✓ Branch 0 taken 3116 times.
✗ Branch 1 not taken.
 3116 return is_simplicial_ ? MESH_TET : MeshCellType(cell_type_[c]);
1876 }
1877
1878 /**
1879 * \brief Gets the descriptor of a cell
1880 * \details The descriptor of a cell is a set of static
1881 * arrays that facilitate some accesses (most client
1882 * code do not need to use this function)
1883 * \param[in] c a cell, in 0..nb()-1
1884 * \return the descriptor of cell \p c
1885 */
1886 const CellDescriptor& descriptor(index_t c) const;
1887
1888 /**
1889 * \brief Gets a descriptor by cell type
1890 * \details The descriptor of a cell is a set of static
1891 * arrays that facilitate some accesses (most client
1892 * code do not need to use this function)
1893 * \param[in] t one of
1894 * MESH_TET, MESH_HEX, MESH_PRISM, MESH_PYRAMID, MESH_CONNECTOR
1895 * \return the descriptor of cell \p c
1896 */
1897 static const CellDescriptor& cell_type_to_cell_descriptor(
1898 MeshCellType t
1899 );
1900
1901 /**
1902 * \brief Gets the number of corners of a cell
1903 * \param[in] c a cell, in 0..nb()-1
1904 * \return the number of corners of cell \p c
1905 */
1906 317784 index_t nb_corners(index_t c) const {
1907
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 317784 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 317784 geo_debug_assert(c < nb());
1908 317784 return descriptor(c).nb_vertices;
1909 }
1910
1911 /**
1912 * \brief Gets the first element for iterating over
1913 * the corners of a cell
1914 * \param[in] c the cell, in 0..nb()-1
1915 * \return the first corner of the cell
1916 */
1917 317784 index_t corners_begin(index_t c) const {
1918
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 317784 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 317784 geo_debug_assert(c < nb());
1919
1/2
✓ Branch 0 taken 317784 times.
✗ Branch 1 not taken.
 317784 return is_simplicial_ ? 4*c : cell_ptr_[c];
1920 }
1921
1922 /**
1923 * \brief Gets the upper limit for iterating over the
1924 * corners of a cell
1925 * \param[in] c the cell, in 0..nb()-1
1926 * \return one position past the last corner of the cell
1927 */
1928 ✗ index_t corners_end(index_t c) const {
1929 ✗ geo_debug_assert(c < nb());
1930 ✗ return is_simplicial_ ? 4*(c+1) : cell_ptr_[c] + nb_corners(c);
1931 }
1932
1933 /**
1934 * \brief Gets a corner of a cell by local vertex index
1935 * \param[in] c the cell, in 0..nb()-1
1936 * \param[in] lv the local vertex index, in 0..nb_corners(c)-1
1937 * \return the corner incident to vertex \p lv in cell \p c
1938 */
1939 317784 index_t corner(index_t c, index_t lv) const {
1940
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 317784 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 317784 geo_debug_assert(c < nb());
1941 // There seems to be a linkage problem under MSVC for the
1942 // following assertion check...
1943 #ifndef GEO_OS_WINDOWS
1944
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 317784 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 317784 geo_debug_assert(lv < nb_corners(c));
1945 #endif
1946 317784 return corners_begin(c) + lv;
1947 }
1948
1949 /**
1950 * \brief Gets the number of facets of a cell
1951 * \param[in] c the cell, in 0..nb()-1
1952 * \return the number of facets of cell \p c
1953 */
1954 28032 index_t nb_facets(index_t c) const {
1955
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 28032 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 28032 geo_debug_assert(c < nb());
1956 28032 return descriptor(c).nb_facets;
1957 }
1958
1959 /**
1960 * \brief Gets the first element for iterating over
1961 * the facets of a cell
1962 * \param[in] c the cell, in 0..nb()-1
1963 * \return the first facet of the cell
1964 */
1965 17664 index_t facets_begin(index_t c) const {
1966
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 17664 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 17664 geo_debug_assert(c < nb());
1967
1/2
✓ Branch 0 taken 17664 times.
✗ Branch 1 not taken.
 17664 return is_simplicial_ ? 4*c : cell_ptr_[c];
1968 }
1969
1970 /**
1971 * \brief Gets the upper limit for iterating over the
1972 * facets of a cell
1973 * \param[in] c the cell, in 0..nb()-1
1974 * \return one position past the last facet of the facet
1975 */
1976 index_t facets_end(index_t c) const {
1977 geo_debug_assert(c < nb());
1978 return is_simplicial_ ? 4*(c+1) : cell_ptr_[c] + nb_facets(c);
1979 }
1980
1981 /**
1982 * \brief Gets a facet of a cell by local facet index
1983 * \param[in] c the cell, in 0..nb()-1
1984 * \param[in] lf the local facet index, in 0..nb_facets(c)-1
1985 * \return the facet \p lf in cell \p c
1986 */
1987 17664 index_t facet(index_t c, index_t lf) const {
1988
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 17664 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 17664 geo_debug_assert(c < nb());
1989
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 17664 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 17664 geo_debug_assert(lf < nb_facets(c));
1990 17664 return facets_begin(c) + lf;
1991 }
1992
1993 /**
1994 * \brief Gets the number of edges in a cell
1995 * \param[in] c the cell, in 0..nb()-1
1996 * \return the number of edges in cell \p c
1997 */
1998 index_t nb_edges(index_t c) const {
1999 return descriptor(c).nb_edges;
2000 }
2001
2002 /**
2003 * \brief Gets a pointer to a cell pointer index by cell index
2004 * \param[in] c cell index
2005 * \return a pointer to the cell pointer index
2006 * \note Normal uses do not call this function
2007 */
2008 index_t* cell_ptr_ptr(index_t c) {
2009 geo_debug_assert(!is_simplicial_);
2010 return &cell_ptr_[c];
2011 }
2012
2013 /**
2014 * \brief Gets a pointer to a cell pointer index by cell index
2015 * \param[in] c cell index
2016 * \return a pointer to the cell pointer index
2017 * \note Normal uses do not call this function
2018 */
2019 const index_t* cell_ptr_ptr(index_t c) const {
2020 geo_debug_assert(!is_simplicial_);
2021 return &cell_ptr_[c];
2022 }
2023
2024 /**
2025 * \brief Gets a pointer to a cell type by cell index
2026 * \param[in] c cell index
2027 * \return a pointer to the cell type
2028 * \note Normal uses do not call this function
2029 */
2030 Numeric::uint8* cell_type_ptr(index_t c) {
2031 geo_debug_assert(!is_simplicial_);
2032 return &cell_type_[c];
2033 }
2034
2035 /**
2036 * \brief Gets a pointer to a cell type by cell index
2037 * \param[in] c cell index
2038 * \return a const pointer to the cell type
2039 * \note Normal uses do not call this function
2040 */
2041 const Numeric::uint8* cell_type_ptr(index_t c) const {
2042 geo_debug_assert(!is_simplicial_);
2043 return &cell_type_[c];
2044 }
2045
2046
2047 protected:
2048 void clear_store(
2049 bool keep_attributes, bool keep_memory = false
2050 ) override;
2051
2052 void resize_store(index_t new_size) override;
2053
2054 ✗ index_t create_sub_element(MeshCellType type) {
2055 ✗ if(!is_simplicial_) {
2056 ✗ cell_ptr_.push_back(NO_CORNER);
2057 ✗ cell_type_.push_back(Numeric::uint8(type));
2058 }
2059 ✗ return MeshSubElementsStore::create_sub_element();
2060 }
2061
2062 4 index_t create_sub_elements(index_t nb, MeshCellType type) {
2063
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 4 if(!is_simplicial_) {
2064 ✗ for(index_t i=0; i<nb; ++i) {
2065 ✗ cell_ptr_.push_back(NO_CORNER);
2066 ✗ cell_type_.push_back(Numeric::uint8(type));
2067 }
2068 }
2069 4 return MeshSubElementsStore::create_sub_elements(nb);
2070 }
2071
2072 60 void copy(
2073 const MeshCellsStore& rhs, bool copy_attributes=true
2074 ) {
2075 60 MeshSubElementsStore::copy(rhs, copy_attributes);
2076 60 is_simplicial_ = rhs.is_simplicial_;
2077 60 cell_type_ = rhs.cell_type_;
2078 60 cell_ptr_ = rhs.cell_ptr_;
2079 60 }
2080
2081 protected:
2082 bool is_simplicial_;
2083 vector<Numeric::uint8> cell_type_;
2084 vector<index_t> cell_ptr_;
2085
2086 protected:
2087 friend class Mesh;
2088 friend class GeogramIOHandler;
2089 };
2090
2091 /*************************************************************************/
2092
2093 /**
2094 * \brief Stores the cell corners of a mesh (low-level store)
2095 * \relates MeshCells
2096 */
2097 class GEOGRAM_API MeshCellCornersStore : public MeshSubElementsStore {
2098 public:
2099 MeshCellCornersStore(Mesh& mesh);
2100
2101 /**
2102 * \brief Gets the vertex that a corner is incident to
2103 * \param[in] c the corner, in 0..nb()-1
2104 * \return the vertex that corner \p c is incident to
2105 */
2106 597408 index_t vertex(index_t c) const {
2107
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 597408 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 597408 geo_debug_assert(c < nb());
2108 597408 return corner_vertex_[c];
2109 }
2110
2111 /**
2112 * \brief Sets the vertex that a corner is incident to
2113 * \param[in] c the corner, in 0..nb()-1
2114 * \param[in] v specifies the vertex that corner \p c is incident to
2115 */
2116 12288 void set_vertex(index_t c, index_t v) {
2117
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 12288 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 12288 geo_debug_assert(c < nb());
2118
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 12288 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 12288 geo_debug_assert(v < vertices_.nb());
2119 12288 corner_vertex_[c] = v;
2120 12288 }
2121
2122 /**
2123 * \brief Gets a pointer to the vertex that a corner is incident to
2124 * \param[in] c the corner
2125 * \return a pointer to the index of the vertex that this corner
2126 * is incident to
2127 * \note Normal uses do not call this function
2128 */
2129 index_t* vertex_index_ptr(index_t c) {
2130 geo_debug_assert(c < nb());
2131 return &(corner_vertex_[c]);
2132 }
2133
2134 /**
2135 * \brief Gets a pointer to the vertex that a corner is incident to
2136 * \param[in] c the corner
2137 * \return a const pointer to the index of the vertex that this corner
2138 * is incident to
2139 * \note Normal uses do not call this function
2140 */
2141 const index_t* vertex_index_ptr(index_t c) const {
2142 geo_debug_assert(c < nb());
2143 return &(corner_vertex_[c]);
2144 }
2145
2146 /**
2147 * \brief Gets the point associated with a corner
2148 * \param[in] c the corner
2149 * \return a reference to the DIM-d point associated with the corner
2150 * \pre vertices.dimension() >= DIM
2151 */
2152 template <index_t DIM=3> vecng<DIM,double>& point(index_t c) {
2153 geo_debug_assert(c < nb());
2154 return vertices_.point<DIM>(vertex(c));
2155 }
2156
2157 /**
2158 * \brief Gets the point associated with a corner
2159 * \param[in] c the corner
2160 * \return const a reference to the DIM-d point associated with
2161 * the corner
2162 * \pre vertices.dimension() >= DIM
2163 */
2164 template <index_t DIM=3> const vecng<DIM,double>& point(
2165 index_t c
2166 ) const {
2167 geo_debug_assert(c < nb());
2168 return vertices_.point(vertex(c));
2169 }
2170
2171 protected:
2172 void clear_store(
2173 bool keep_attributes, bool keep_memory = false
2174 ) override;
2175
2176 void resize_store(index_t new_size) override;
2177
2178 ✗ index_t create_sub_element(index_t v) {
2179 ✗ corner_vertex_.push_back(v);
2180 ✗ return MeshSubElementsStore::create_sub_element();
2181 }
2182
2183 4 index_t create_sub_elements(index_t nb) {
2184
2/2
✓ Branch 0 taken 6144 times.
✓ Branch 1 taken 4 times.
 6148 for(index_t i=0; i<nb; ++i) {
2185 6144 corner_vertex_.push_back(NO_VERTEX);
2186 }
2187 4 return MeshSubElementsStore::create_sub_elements(nb);
2188 }
2189
2190 60 void copy(
2191 const MeshCellCornersStore& rhs, bool copy_attributes=true
2192 ) {
2193 60 MeshSubElementsStore::copy(rhs, copy_attributes);
2194 60 corner_vertex_ = rhs.corner_vertex_;
2195 60 }
2196
2197 protected:
2198 MeshVertices& vertices_;
2199 vector<index_t> corner_vertex_;
2200
2201 friend class MeshCells;
2202 friend class Mesh;
2203 friend class GeogramIOHandler;
2204 };
2205
2206 /*************************************************************************/
2207
2208 /**
2209 * \brief Stores the cell facets of a mesh (low-level store)
2210 * \relates MeshCells
2211 */
2212 class GEOGRAM_API MeshCellFacetsStore : public MeshSubElementsStore {
2213 public:
2214 /**
2215 * \brief MeshCellFacetsStore constructor
2216 * \param[in] mesh the mesh that this MeshCellFacetsStore is attached to
2217 */
2218 MeshCellFacetsStore(Mesh& mesh);
2219
2220 /**
2221 * \brief Gets a cell adjacent to a facet
2222 * \param[in] f the facet, in 0..nb()-1
2223 * \return the cell adjacent to facet \p f, or NO_FACET if \p f
2224 * is on the border
2225 */
2226 12288 index_t adjacent_cell(index_t f) const {
2227
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 12288 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 12288 geo_debug_assert(f < nb());
2228 12288 return adjacent_cell_[f];
2229 }
2230
2231 /**
2232 * \brief Sets a cell adjacent to a facet
2233 * \param[in] f the facet, in 0..nb()-1
2234 * \param[in] c specifies the cell adjacent
2235 * to facet \p f, or is set to NO_FACET
2236 * if \p f is on the border
2237 */
2238 11520 void set_adjacent_cell(index_t f, index_t c) {
2239
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 11520 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 11520 geo_debug_assert(f < nb());
2240
4/10
✓ Branch 0 taken 5376 times.
✓ Branch 1 taken 6144 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 5376 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 11520 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 11520 geo_debug_assert(c == NO_CELL || c < cells_.nb());
2241 11520 adjacent_cell_[f] = c;
2242 11520 }
2243
2244 /**
2245 * \brief Gets a const pointer to a cell adjacent to a facet
2246 * \param[in] f the facet, in 0..nb()-1
2247 * \return a const pointer to the cell adjacent to facet \p f,
2248 * or NO_FACET if \p f is on the border
2249 */
2250 const index_t* adjacent_cell_ptr(index_t f) const {
2251 geo_debug_assert(f < nb());
2252 return &adjacent_cell_[f];
2253 }
2254
2255 /**
2256 * \brief Gets a pointer to a cell adjacent to a facet
2257 * \param[in] f the facet, in 0..nb()-1
2258 * \return a pointer to the cell adjacent to facet \p f,
2259 * or NO_FACET if \p f is on the border
2260 */
2261 index_t* adjacent_cell_ptr(index_t f) {
2262 geo_debug_assert(f < nb());
2263 return &adjacent_cell_[f];
2264 }
2265
2266 protected:
2267 void clear_store(
2268 bool keep_attributes, bool keep_memory = false
2269 ) override;
2270
2271 void resize_store(index_t new_size) override;
2272
2273 ✗ index_t create_sub_element(index_t c = NO_CELL) {
2274 ✗ adjacent_cell_.push_back(c);
2275 ✗ return MeshSubElementsStore::create_sub_element();
2276 }
2277
2278 4 index_t create_sub_elements(index_t nb) {
2279
2/2
✓ Branch 0 taken 6144 times.
✓ Branch 1 taken 4 times.
 6148 for(index_t i=0; i<nb; ++i) {
2280 6144 adjacent_cell_.push_back(NO_CELL);
2281 }
2282 4 return MeshSubElementsStore::create_sub_elements(nb);
2283 }
2284
2285 60 void copy(
2286 const MeshCellFacetsStore& rhs, bool copy_attributes=true
2287 ) {
2288 60 MeshSubElementsStore::copy(rhs, copy_attributes);
2289 60 adjacent_cell_ = rhs.adjacent_cell_;
2290 60 }
2291
2292 protected:
2293 MeshVertices& vertices_;
2294 MeshCellsStore& cells_;
2295 vector<index_t> adjacent_cell_;
2296
2297 friend class MeshCells;
2298 friend class Mesh;
2299 friend class GeogramIOHandler;
2300 };
2301
2302 /*************************************************************************/
2303
2304 /**
2305 * \brief The cells of a mesh.
2306 * \relates Mesh
2307 */
2308 class GEOGRAM_API MeshCells : public MeshCellsStore, public MeshElements {
2309 public:
2310 /**
2311 * \brief MeshCells constructor
2312 * \param[in] mesh the mesh this MeshCells is attached to
2313 */
2314 MeshCells(Mesh& mesh);
2315
2316 /**
2317 * \brief Gets the number of vertices of a cell
2318 * \param[in] c a cell, in 0..nb()-1
2319 * \return the number of vertices of cell \p c
2320 */
2321 ✗ index_t nb_vertices(index_t c) const {
2322 ✗ return nb_corners(c);
2323 }
2324
2325 /**
2326 * \brief Gets a vertex of a cell by local vertex index
2327 * \param[in] c the cell, in 0..nb()-1
2328 * \param[in] lv local vertex index, in 0..nb_vertices(c)-1
2329 * \return the vertex \p lv of cell \p c
2330 */
2331 301272 index_t vertex(index_t c, index_t lv) const {
2332 301272 return cell_corners_.vertex(corner(c,lv));
2333 }
2334
2335 /**
2336 * \brief Sets a vertex of a cell by local vertex index
2337 * \param[in] c the cell, in 0..nb()-1
2338 * \param[in] lv local vertex index, in 0..nb_vertices(c)-1
2339 * \param[in] v specifies the vertex \p lv of cell \p c
2340 */
2341 6144 void set_vertex(index_t c, index_t lv, index_t v) {
2342 6144 cell_corners_.set_vertex(corner(c,lv),v);
2343 6144 }
2344
2345 /**
2346 * \brief Gets a point by cell and local vertex index
2347 * \param[in] c the cell
2348 * \param[in] lv the local vertex index in \p c
2349 * \return a const reference to the point corresponding to
2350 * the \p lv%th vertex of cell \p c
2351 * \pre lv < nb_vertices(c)
2352 */
2353 21504 template <index_t DIM=3> const vecng<DIM,double>& point(
2354 index_t c, index_t lv
2355 ) const {
2356 21504 return vertices_.point<DIM>(vertex(c,lv));
2357 }
2358
2359 /**
2360 * \brief Gets a point by cell and local vertex index
2361 * \param[in] c the cell
2362 * \param[in] lv the local vertex index in \p c
2363 * \return a reference to the point corresponding to
2364 * the \p lv%th vertex of cell \p c
2365 * \pre lv < nb_vertices(c)
2366 */
2367 template <index_t DIM=3> vecng<DIM,double>& point(
2368 index_t c, index_t lv
2369 ) {
2370 return vertices_.point<DIM>(vertex(c,lv));
2371 }
2372
2373 /**
2374 * \brief Gets a cell adjacent to another one by local facet index
2375 * \param[in] c the cell, in 0..nb()-1
2376 * \param[in] lf local facet index, in 0..nb_facets(c)-1
2377 * \return the cell adjacent to \p c along facet \p lf or NO_CELL
2378 * if no such cell exists
2379 */
2380 12288 index_t adjacent(index_t c, index_t lf) const {
2381 12288 return cell_facets_.adjacent_cell(facet(c,lf));
2382 }
2383
2384 /**
2385 * \brief Sets a cell adjacent to another one by local facet index
2386 * \param[in] c the cell, in 0..nb()-1
2387 * \param[in] lf local facet index, in 0..nb_facets(c)-1
2388 * \param[in] c2 specifies the cell adjacent to \p c along
2389 * facet \p lf or NO_CELL if no such cell exists
2390 */
2391 5376 void set_adjacent(index_t c, index_t lf, index_t c2) {
2392 5376 cell_facets_.set_adjacent_cell(facet(c,lf),c2);
2393 5376 }
2394
2395 /**
2396 * \brief Gets the number of vertices in a cell facet
2397 * \param[in] c the cell, in 0..nb()-1
2398 * \param[in] lf the local facet index, in 0..nb_facets(c)-1
2399 * \return the number of vertices of facet \p lf in cell \p c
2400 */
2401 10368 index_t facet_nb_vertices(index_t c, index_t lf) const {
2402
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 10368 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 10368 geo_debug_assert(lf < nb_facets(c));
2403 10368 return descriptor(c).nb_vertices_in_facet[lf];
2404 }
2405
2406 /**
2407 * \brief Gets a vertex of a cell by local facet index and
2408 * local vertex index in the facet
2409 * \param[in] c the cell, in 0..nb()-1
2410 * \param[in] lf the local facet index, in 0..nb_facets(c)-1
2411 * \param[in] lv the local vertex index, in 0..facet_nb_vertices(c,lf)-1
2412 * \return vertex \p lv of facet \p lf in cell \p c
2413 */
2414 10368 index_t facet_vertex(index_t c, index_t lf, index_t lv) const {
2415
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 10368 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 10368 geo_debug_assert(lv < facet_nb_vertices(c, lf));
2416 10368 return cell_corners_.vertex(
2417 10368 corner(c, descriptor(c).facet_vertex[lf][lv])
2418 10368 );
2419 }
2420 /**
2421 * \brief Gets a corner of a cell by local facet index and
2422 * local corner index in the facet
2423 * \param[in] c the cell, in 0..nb()-1
2424 * \param[in] lf the local facet index, in 0..nb_facets(c)-1
2425 * \param[in] lc the local corner index, in 0..facet_nb_vertices(c,lf)-1
2426 * \return corner \p lc of facet \p lf in cell \p c
2427 */
2428 index_t facet_corner(index_t c, index_t lf, index_t lc) const {
2429 geo_debug_assert(lc < facet_nb_vertices(c, lf));
2430 return corner(c, descriptor(c).facet_vertex[lf][lc]);
2431 }
2432
2433 /**
2434 * \brief Gets a cell vertex by local edge index and local
2435 * vertex index in the edge
2436 * \param[in] c the cell, in 0..nb()-1
2437 * \param[in] le the local edge index, in 0..nb_edges(c)-1
2438 * \param[in] lv the local index in the edge, one of 0,1
2439 * \return vertex \p lv of edge \p le in cell \p c
2440 */
2441 index_t edge_vertex(index_t c, index_t le, index_t lv) const {
2442 geo_debug_assert(le < nb_edges(c));
2443 geo_debug_assert(lv < 2);
2444 return cell_corners_.vertex(
2445 corner(c,descriptor(c).edge_vertex[le][lv])
2446 );
2447 }
2448
2449 /**
2450 * \brief Gets a cell local facet index by local edge index and local
2451 * facet index in the edge
2452 * \param[in] c the cell, in 0..nb()-1
2453 * \param[in] le the local edge index, in 0..nb_edges(c)-1
2454 * \param[in] lf the local index in the edge, one of 0,1
2455 * \return the local facet index of the facet adjacent to the edge.
2456 * - If \p lf=0, it gets the facet on the left of the oriented edge.
2457 * - If \p lf=1, it gets the facet on the right of the oriented edge.
2458 */
2459 index_t edge_adjacent_facet(index_t c, index_t le, index_t lf) const {
2460 geo_debug_assert(le < nb_edges(c));
2461 geo_debug_assert(lf < 2);
2462 return descriptor(c).edge_adjacent_facet[le][lf];
2463 }
2464
2465 /**
2466 * \brief Gets the corners of a cell.
2467 * \param[in] c the index of the cell.
2468 * \return a range with all the corners of the facet.
2469 */
2470 ✗ index_range corners(index_t c) const {
2471 ✗ geo_debug_assert(c < nb());
2472 ✗ return index_range(
2473 ✗ index_as_iterator(corners_begin(c)),
2474 ✗ index_as_iterator(corners_end(c))
2475 ✗ );
2476 }
2477
2478 /**
2479 * \brief Gets the facets of a cell.
2480 * \param[in] c the index of the cell.
2481 * \return a range with all the global indices of the cell facets.
2482 */
2483 index_range facets(index_t c) const {
2484 geo_debug_assert(c < nb());
2485 return index_range(
2486 index_as_iterator(facets_begin(c)),
2487 index_as_iterator(facets_end(c))
2488 );
2489 }
2490
2491 #ifndef MESH_NO_SYNTAXIC_SUGAR
2492
2493 /**
2494 * \brief Gets the points associated with the vertices of a cell.
2495 * \param[in] cell the index of the cell.
2496 * \return a range with the points that correspond to the vertices of
2497 * the cell, returned as const references.
2498 */
2499 ✗ template <index_t DIM=3> auto points(index_t cell) const {
2500 ✗ geo_debug_assert(cell < nb());
2501 typedef vecng<DIM,double> vecn;
2502 ✗ return transform_range_ref(
2503 ✗ corners(cell), [this](index_t c)->const vecn& {
2504 ✗ index_t v = cell_corners_.vertex(c);
2505 ✗ return vertices_.point<DIM>(v);
2506 }
2507 ✗ );
2508 }
2509
2510 /**
2511 * \brief Gets the points associated with the vertices of a cell.
2512 * \param[in] cell the index of the cell.
2513 * \return a range with the points that correspond to the vertices of
2514 * the cell, returned as modifiable references.
2515 */
2516 ✗ template <index_t DIM=3> auto points(index_t cell) {
2517 ✗ geo_debug_assert(cell < nb());
2518 typedef vecng<DIM,double> vecn;
2519 ✗ return transform_range_ref(
2520 ✗ corners(cell), [this](index_t c)->vecn& {
2521 ✗ index_t v = cell_corners_.vertex(c);
2522 ✗ return vertices_.point<DIM>(v);
2523 }
2524 ✗ );
2525 }
2526
2527 /**
2528 * \brief Gets the cells adjacent to a given cell.
2529 * \param[in] c the index of the cell.
2530 * \return a range with all the indices of the cels adjacent to this
2531 * cell. It will output exactly one item per facet of the cell. Facets
2532 * on the border will output NO_INDEX (no adjacent cell).
2533 */
2534 auto adjacent(index_t c) const {
2535 geo_debug_assert(c < nb());
2536 return transform_range(
2537 facets(c), [this](index_t f)->index_t {
2538 return cell_facets_.adjacent_cell(f);
2539 }
2540 );
2541 }
2542
2543 #endif
2544
2545 void clear(
2546 bool keep_attributes=true, bool keep_memory=false
2547 ) override;
2548
2549 void delete_elements(
2550 vector<index_t>& to_delete,
2551 bool remove_isolated_vertices=true
2552 ) override;
2553
2554 void permute_elements(vector<index_t>& permutation) override;
2555
2556 /**
2557 * \brief Creates a contiguous chunk of cells of the
2558 * same type
2559 * \param[in] nb_cells number of cells to create
2560 * \param[in] type type of the cells to create, one of
2561 * MESH_TET, MESH_HEX, MESH_PRISM, MESH_PYRAMID, MESH_CONNECTOR.
2562 * \return the first created cell
2563 */
2564 4 index_t create_cells(index_t nb_cells, MeshCellType type) {
2565
2566
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 4 if(nb_cells == 0) {
2567 ✗ return NO_CELL;
2568 }
2569
2570
2571
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 4 if(type != MESH_TET) {
2572 ✗ is_not_simplicial();
2573 }
2574
2575 4 const CellDescriptor& desc = cell_type_to_cell_descriptor(type);
2576
2577 // Note: there is padding, the same number of corners and
2578 // faces is created for each cell, so that a single cell
2579 // pointer is used for both.
2580
2581 4 index_t cell_size = std::max(desc.nb_vertices, desc.nb_facets);
2582 4 index_t first_cell = nb();
2583 4 index_t co = cell_corners_.nb();
2584
2585 4 cell_corners_.create_sub_elements(
2586 nb_cells*cell_size
2587 );
2588
2589 4 cell_facets_.create_sub_elements(
2590 nb_cells*cell_size
2591 );
2592
2593 4 index_t result = create_sub_elements(nb_cells, type);
2594
2595
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 4 if(!is_simplicial_) {
2596 ✗ for(index_t c=first_cell; c<=first_cell+nb_cells; ++c) {
2597 ✗ cell_ptr_[c] = co;
2598 ✗ co += cell_size;
2599 }
2600
2601 ✗ geo_debug_assert(cell_ptr_.size() == nb()+1);
2602 ✗ geo_debug_assert(cell_ptr_[nb()] == cell_corners_.nb());
2603 ✗ geo_debug_assert(cell_ptr_[nb()] == cell_facets_.nb());
2604 }
2605
2606 4 return result;
2607 }
2608
2609 /**
2610 * \brief Creates a contiguous chunk of tetrahedra
2611 * \param[in] nb_tets number of tetrahedra to create
2612 * \return the first created tetrahedron
2613 */
2614 4 index_t create_tets(index_t nb_tets) {
2615 4 return create_cells(nb_tets, MESH_TET);
2616 }
2617
2618 /**
2619 * \brief Creates a contiguous chunk of hexahedra
2620 * \param[in] nb_hexes number of hexahedra to create
2621 * \return the first created hexahedron
2622 */
2623 ✗ index_t create_hexes(index_t nb_hexes) {
2624 ✗ return create_cells(nb_hexes, MESH_HEX);
2625 }
2626
2627 /**
2628 * \brief Creates a contiguous chunk of prisms
2629 * \param[in] nb_prisms number of prisms to create
2630 * \return the first created prism
2631 */
2632 ✗ index_t create_prisms(index_t nb_prisms) {
2633 ✗ return create_cells(nb_prisms, MESH_PRISM);
2634 }
2635
2636 /**
2637 * \brief Creates a contiguous chunk of pyramids
2638 * \param[in] nb_pyramids number of pyramids to create
2639 * \return the first created pyramid
2640 */
2641 ✗ index_t create_pyramids(index_t nb_pyramids) {
2642 ✗ return create_cells(nb_pyramids, MESH_PYRAMID);
2643 }
2644
2645 /**
2646 * \brief Creates a tetrahedron
2647 * \param[in] v1 , v2 , v3 , v4 the vertices of the tetrahedron,
2648 * all in 0 .. mesh.vertices.nb()-1
2649 * \param[in] adj1 , adj2 , adj3 , adj4
2650 * adjacent cells, or NO_CELL if unspecified / on border
2651 * \return the created tetrahedron
2652 */
2653 ✗ index_t create_tet(
2654 index_t v1, index_t v2, index_t v3, index_t v4,
2655 index_t adj1 = NO_CELL,
2656 index_t adj2 = NO_CELL,
2657 index_t adj3 = NO_CELL,
2658 index_t adj4 = NO_CELL
2659 ) {
2660 ✗ cell_corners_.create_sub_element(v1);
2661 ✗ cell_corners_.create_sub_element(v2);
2662 ✗ cell_corners_.create_sub_element(v3);
2663 ✗ cell_corners_.create_sub_element(v4);
2664 ✗ cell_facets_.create_sub_element(adj1);
2665 ✗ cell_facets_.create_sub_element(adj2);
2666 ✗ cell_facets_.create_sub_element(adj3);
2667 ✗ cell_facets_.create_sub_element(adj4);
2668 ✗ index_t result = create_sub_element(MESH_TET);
2669 ✗ if(!is_simplicial_) {
2670 ✗ cell_ptr_[nb()] = cell_corners_.nb();
2671 }
2672 ✗ geo_debug_assert(cell_facets_.nb() == cell_corners_.nb());
2673 ✗ return result;
2674 }
2675
2676 /**
2677 * \brief Creates an hexahedron
2678 * \param[in] v1 , v2 , v3 , v4 , v5 , v6 , v7 , v8
2679 * the vertices of the hexahedron,
2680 * all in 0 .. mesh.vertices.nb()-1
2681 * \param[in] adj1 , adj2 , adj3 , adj4 , adj5 , adj6
2682 * adjacent cells, or NO_CELL if unspecified / on border
2683 * \return the created hexahedron
2684 */
2685 ✗ index_t create_hex(
2686 index_t v1, index_t v2, index_t v3, index_t v4,
2687 index_t v5, index_t v6, index_t v7, index_t v8,
2688 index_t adj1 = NO_CELL,
2689 index_t adj2 = NO_CELL,
2690 index_t adj3 = NO_CELL,
2691 index_t adj4 = NO_CELL,
2692 index_t adj5 = NO_CELL,
2693 index_t adj6 = NO_CELL
2694 ) {
2695 ✗ is_not_simplicial();
2696 ✗ cell_corners_.create_sub_element(v1);
2697 ✗ cell_corners_.create_sub_element(v2);
2698 ✗ cell_corners_.create_sub_element(v3);
2699 ✗ cell_corners_.create_sub_element(v4);
2700 ✗ cell_corners_.create_sub_element(v5);
2701 ✗ cell_corners_.create_sub_element(v6);
2702 ✗ cell_corners_.create_sub_element(v7);
2703 ✗ cell_corners_.create_sub_element(v8);
2704 ✗ cell_facets_.create_sub_element(adj1);
2705 ✗ cell_facets_.create_sub_element(adj2);
2706 ✗ cell_facets_.create_sub_element(adj3);
2707 ✗ cell_facets_.create_sub_element(adj4);
2708 ✗ cell_facets_.create_sub_element(adj5);
2709 ✗ cell_facets_.create_sub_element(adj6);
2710 ✗ cell_facets_.create_sub_element(NO_CELL); // padding
2711 ✗ cell_facets_.create_sub_element(NO_CELL); // padding
2712 ✗ index_t result = create_sub_element(MESH_HEX);
2713 ✗ cell_ptr_[nb()] = cell_corners_.nb();
2714 ✗ geo_debug_assert(cell_facets_.nb() == cell_corners_.nb());
2715 ✗ return result;
2716 }
2717
2718 /**
2719 * \brief Creates a prism
2720 * \param[in] v1 , v2 , v3 , v4 , v5 , v6
2721 * the vertices of the prism
2722 * all in 0 .. mesh.vertices.nb()-1
2723 * \param[in] adj1 , adj2 , adj3 , adj4 , adj5
2724 * adjacent cells, or NO_CELL if unspecified / on border
2725 * \return the created prism
2726 */
2727 ✗ index_t create_prism(
2728 index_t v1, index_t v2,
2729 index_t v3, index_t v4,
2730 index_t v5, index_t v6,
2731 index_t adj1 = NO_CELL,
2732 index_t adj2 = NO_CELL,
2733 index_t adj3 = NO_CELL,
2734 index_t adj4 = NO_CELL,
2735 index_t adj5 = NO_CELL
2736 ) {
2737 ✗ is_not_simplicial();
2738 ✗ cell_corners_.create_sub_element(v1);
2739 ✗ cell_corners_.create_sub_element(v2);
2740 ✗ cell_corners_.create_sub_element(v3);
2741 ✗ cell_corners_.create_sub_element(v4);
2742 ✗ cell_corners_.create_sub_element(v5);
2743 ✗ cell_corners_.create_sub_element(v6);
2744 ✗ cell_facets_.create_sub_element(adj1);
2745 ✗ cell_facets_.create_sub_element(adj2);
2746 ✗ cell_facets_.create_sub_element(adj3);
2747 ✗ cell_facets_.create_sub_element(adj4);
2748 ✗ cell_facets_.create_sub_element(adj5);
2749 ✗ cell_facets_.create_sub_element(NO_CELL); // padding
2750 ✗ index_t result = create_sub_element(MESH_PRISM);
2751 ✗ cell_ptr_[nb()] = cell_corners_.nb();
2752 ✗ geo_debug_assert(cell_facets_.nb() == cell_corners_.nb());
2753 ✗ return result;
2754 }
2755
2756 /**
2757 * \brief Creates a pyramid
2758 * \param[in] v1 , v2 , v3 , v4 , v5
2759 * the vertices of the pyramid
2760 * all in 0 .. mesh.vertices.nb()-1
2761 * \param[in] adj1 , adj2 , adj3 , adj4 , adj5
2762 * adjacent cells, or NO_CELL if unspecified / on border
2763 * \return the created pyramid
2764 */
2765 ✗ index_t create_pyramid(
2766 index_t v1, index_t v2, index_t v3, index_t v4, index_t v5,
2767 index_t adj1 = NO_CELL,
2768 index_t adj2 = NO_CELL,
2769 index_t adj3 = NO_CELL,
2770 index_t adj4 = NO_CELL,
2771 index_t adj5 = NO_CELL
2772 ) {
2773 ✗ is_not_simplicial();
2774 ✗ cell_corners_.create_sub_element(v1);
2775 ✗ cell_corners_.create_sub_element(v2);
2776 ✗ cell_corners_.create_sub_element(v3);
2777 ✗ cell_corners_.create_sub_element(v4);
2778 ✗ cell_corners_.create_sub_element(v5);
2779 ✗ cell_facets_.create_sub_element(adj1);
2780 ✗ cell_facets_.create_sub_element(adj2);
2781 ✗ cell_facets_.create_sub_element(adj3);
2782 ✗ cell_facets_.create_sub_element(adj4);
2783 ✗ cell_facets_.create_sub_element(adj5);
2784 ✗ index_t result = create_sub_element(MESH_PYRAMID);
2785 ✗ cell_ptr_[nb()] = cell_corners_.nb();
2786 ✗ geo_debug_assert(cell_facets_.nb() == cell_corners_.nb());
2787 ✗ return result;
2788 }
2789
2790 /**
2791 * \brief Creates a connector
2792 * \details Connector are automatically
2793 * created by connect() (most client codes
2794 * do not use this function)
2795 * \param[in] v1 , v2 , v3 , v4
2796 * the vertices of the connector
2797 * all in 0 .. mesh.vertices.nb()-1
2798 * \param[in] adj1 , adj2 , adj3
2799 * adjacent cells, or NO_CELL if unspecified / on border
2800 * \return the created connector
2801 */
2802 ✗ index_t create_connector(
2803 index_t v1, index_t v2, index_t v3, index_t v4,
2804 index_t adj1 = NO_CELL,
2805 index_t adj2 = NO_CELL,
2806 index_t adj3 = NO_CELL
2807 ) {
2808 ✗ is_not_simplicial();
2809 ✗ cell_corners_.create_sub_element(v1);
2810 ✗ cell_corners_.create_sub_element(v2);
2811 ✗ cell_corners_.create_sub_element(v3);
2812 ✗ cell_corners_.create_sub_element(v4);
2813 ✗ cell_facets_.create_sub_element(adj1);
2814 ✗ cell_facets_.create_sub_element(adj2);
2815 ✗ cell_facets_.create_sub_element(adj3);
2816 ✗ cell_facets_.create_sub_element(NO_CELL); // padding
2817 ✗ index_t result = create_sub_element(MESH_CONNECTOR);
2818 ✗ cell_ptr_[nb()] = cell_corners_.nb();
2819 ✗ geo_debug_assert(cell_facets_.nb() == cell_corners_.nb());
2820 ✗ return result;
2821 }
2822
2823 /**
2824 * \brief Connects the cells.
2825 * \details This creates as needed the connectors that represent
2826 * non-conformal connections between a quadrilateral facet and
2827 * two triangular facets.
2828 * \param[in] remove_trivial_slivers if set, this removes the
2829 * slivers that are adjacent to a quadrilateral facet.
2830 * \param[in] verbose_if_OK if set, says OK if no bad connector
2831 * configuration was detected.
2832 */
2833 void connect(
2834 bool remove_trivial_slivers = true, bool verbose_if_OK=false
2835 );
2836
2837 /**
2838 * \brief Replaces the surfacic part of this mesh
2839 * with the borders of the volumetric part.
2840 */
2841 void compute_borders();
2842
2843 /**
2844 * \brief Replaces the surfacic part of this mesh
2845 * with the borders of the volumetric part.
2846 * \param[out] facet_cell on exit, stores the
2847 * index of the cell adjacent to the facet
2848 * on the border.
2849 */
2850 void compute_borders(Attribute<index_t>& facet_cell);
2851
2852 /**
2853 * \brief Copies a tetrahedron mesh into this Mesh.
2854 * \details Tetrahedron adjacences are not computed.
2855 * \param[in] dim dimension of the vertices
2856 * \param[in] vertices coordinates of the vertices
2857 * \param[in] tets tetrahedron to vertex links
2858 * \param[in] steal_args if set, vertices and tets
2859 * are 'stolen' from the arguments
2860 * (using vector::swap).
2861 */
2862 void assign_tet_mesh(
2863 coord_index_t dim,
2864 vector<double>& vertices,
2865 vector<index_t>& tets,
2866 bool steal_args
2867 );
2868
2869 /**
2870 * \brief Copies a tetrahedron mesh into this Mesh.
2871 * \details Tetrahedron adjacences are not computed.
2872 * \param[in] tets tetrahedron to vertex links
2873 * \param[in] steal_args if set, vertices and tets
2874 * are 'stolen' from the arguments
2875 * (using vector::swap).
2876 */
2877 void assign_tet_mesh(
2878 vector<index_t>& tets,
2879 bool steal_args
2880 );
2881
2882 void pop() override;
2883
2884 3315888 index_t tet_adjacent(index_t t, index_t lf) const {
2885
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 3315888 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 3315888 geo_debug_assert(is_simplicial_);
2886
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 3315888 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 3315888 geo_debug_assert(t < nb());
2887
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 3315888 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 3315888 geo_debug_assert(lf < 4);
2888 3315888 return cell_facets_.adjacent_cell_[4*t+lf];
2889 }
2890
2891 35367 index_t find_tet_adjacent(index_t t, index_t t2) const {
2892
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 35367 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 35367 geo_debug_assert(is_simplicial_);
2893
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 35367 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 35367 geo_debug_assert(t < nb());
2894
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 35367 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 35367 geo_debug_assert(t2 < nb());
2895
1/2
✓ Branch 0 taken 88293 times.
✗ Branch 1 not taken.
 88293 for(index_t lf=0; lf<4; ++lf) {
2896
2/2
✓ Branch 1 taken 35367 times.
✓ Branch 2 taken 52926 times.
 88293 if(cell_facets_.adjacent_cell_[4*t+lf] == t2) {
2897 35367 return lf;
2898 }
2899 }
2900 ✗ return NO_FACET;
2901 }
2902
2903 5238022 index_t tet_vertex(index_t t, index_t lv) const {
2904
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 5238022 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 5238022 geo_debug_assert(is_simplicial_);
2905
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 5238022 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 5238022 geo_debug_assert(t < nb());
2906
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 5238022 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 5238022 geo_debug_assert(lv < 4);
2907 5238022 return cell_corners_.corner_vertex_[4*t+lv];
2908 }
2909
2910 index_t find_tet_vertex(index_t t, index_t v) const {
2911 geo_debug_assert(is_simplicial_);
2912 geo_debug_assert(t < nb());
2913 geo_debug_assert(v < vertices_.nb());
2914 for(index_t lv=0; lv<4; ++lv) {
2915 if(cell_corners_.corner_vertex_[4*t+lv] == v) {
2916 return lv;
2917 }
2918 }
2919 return NO_VERTEX;
2920 }
2921
2922 /**
2923 * \brief Gets a vertex of a tetrahedron by local facet
2924 * index and local vertex index in facet.
2925 * \param[in] t global index of the tetrahedron
2926 * \param[in] lf local facet index (0,1,2 or 3)
2927 * \param[in] lv local vertex index in facet (0,1 or 2)
2928 * \return the global index of vertex \p lv in facet \p lf of
2929 * tetrahedron \p t
2930 * \pre are_simplices()
2931 */
2932 279624 index_t tet_facet_vertex(
2933 index_t t, index_t lf, index_t lv
2934 ) const {
2935
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 279624 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 279624 geo_debug_assert(is_simplicial_);
2936
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 279624 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 279624 geo_debug_assert(t < nb());
2937
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 279624 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 279624 geo_debug_assert(lf < 4);
2938
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 279624 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 279624 geo_debug_assert(lv < 3);
2939 838872 return cell_corners_.vertex(
2940 279624 4 * t + local_tet_facet_vertex_index(lf,lv)
2941 279624 );
2942 }
2943
2944 /**
2945 * \brief Finds the local index of a facet in a tetrahedron
2946 * by the global indices of its vertices.
2947 * \param[in] t index of the tetrahedron
2948 * \param[in] v1 global index of the first vertex
2949 * \param[in] v2 global index of the second vertex
2950 * \param[in] v3 global index of the third vertex
2951 * \return the local index (0,1,2 or 3) of the facet of
2952 * \p t that has \p v1, \p v2, \p v3 as vertices modulo a
2953 * circular permutation, or NO_FACET if such a facet does not
2954 * exist in \p t.
2955 * \pre are_simplices()
2956 */
2957 24100 index_t find_tet_facet(
2958 index_t t, index_t v1, index_t v2, index_t v3
2959 ) const {
2960
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 24100 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 24100 geo_debug_assert(is_simplicial_);
2961
2/2
✓ Branch 0 taken 92440 times.
✓ Branch 1 taken 21412 times.
 113852 for(index_t lf = 0; lf < 4; ++lf) {
2962 92440 index_t w1 = tet_facet_vertex(t, lf, 0);
2963 92440 index_t w2 = tet_facet_vertex(t, lf, 1);
2964 92440 index_t w3 = tet_facet_vertex(t, lf, 2);
2965
2/2
✓ Branch 0 taken 6332 times.
✓ Branch 1 taken 86108 times.
 92440 if(
2966
6/6
✓ Branch 0 taken 1008 times.
✓ Branch 1 taken 5324 times.
✓ Branch 2 taken 196 times.
✓ Branch 3 taken 812 times.
✓ Branch 4 taken 5212 times.
✓ Branch 5 taken 86416 times.
 92440 (v1 == w1 && v2 == w2 && v3 == w3) ||
2967
6/6
✓ Branch 0 taken 832 times.
✓ Branch 1 taken 4380 times.
✓ Branch 2 taken 324 times.
✓ Branch 3 taken 508 times.
✓ Branch 4 taken 7336 times.
✓ Branch 5 taken 83784 times.
 91628 (v1 == w2 && v2 == w3 && v3 == w1) ||
2968
4/4
✓ Branch 0 taken 1616 times.
✓ Branch 1 taken 5720 times.
✓ Branch 2 taken 1368 times.
✓ Branch 3 taken 248 times.
 7336 (v1 == w3 && v2 == w1 && v3 == w2)
2969 ) {
2970 2688 return lf;
2971 }
2972 }
2973 21412 return NO_FACET;
2974 }
2975
2976 /**
2977 * \brief Gives the local index of a vertex in a
2978 * tetrahedron from its facet and vertex local indices.
2979 * \param[in] lf local facet index (0,1,2 or 3)
2980 * \param[in] lv local vertex index in \p lf (0,1 or 2)
2981 * \return the local vertex index (0,1,2 or 3) of the
2982 * \p lv%th vertex in facet \p lf
2983 */
2984 2017098 static index_t local_tet_facet_vertex_index(index_t lf, index_t lv) {
2985
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 2017098 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 2017098 geo_debug_assert(lf < 4);
2986
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 2017098 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 2017098 geo_debug_assert(lv < 3);
2987 2017098 return MeshCellDescriptors::tet_descriptor.facet_vertex[lf][lv];
2988 }
2989
2990 protected:
2991
2992 /**
2993 * \brief Indicates that the stored elements are no
2994 * longer only tetrahedra.
2995 * \details Creates the cell pointers and cell types
2996 * for the pre-existing cells if any.
2997 */
2998 ✗ void is_not_simplicial() {
2999 ✗ if(is_simplicial_) {
3000 ✗ is_simplicial_ = false;
3001 ✗ cell_ptr_.resize(nb()+1);
3002 ✗ cell_type_.assign(nb(), MESH_TET);
3003 ✗ for(index_t c=0; c<cell_ptr_.size(); ++c) {
3004 ✗ cell_ptr_[c] = 4*c;
3005 }
3006 }
3007 ✗ }
3008
3009 /**
3010 * \brief Tests whether two cell facets can be connected.
3011 * \details Two cell facets can be connected if they have the
3012 * same vertices in reverse order.
3013 * \param[in] c1 index of the first cell
3014 * \param[in] f1 index of the first facet in \p c1
3015 * \param[in] c2 index of the second cell
3016 * \param[in] f2 index of the second facet in \p c2
3017 * \retval true if \p c1 and \p c2 can be connected by \p f1 and \p f2
3018 * \retval false otherwise
3019 */
3020 bool facets_match(
3021 index_t c1, index_t f1, index_t c2, index_t f2
3022 ) const;
3023
3024 /**
3025 * \brief Finds the local index of a vertex in a cell.
3026 * \param[in] c index of the cell
3027 * \param[in] v global index of the vertex
3028 * \return the local index
3029 * (in 0..cell_nb_vertices(c)-1) of the vertex in
3030 * cell \p c or NO_VERTEX if \p c is not incident to \p v
3031 */
3032 ✗ index_t find_cell_vertex(index_t c, index_t v) const {
3033 ✗ geo_debug_assert(c < nb());
3034 ✗ geo_debug_assert(v < vertices_.nb());
3035 ✗ for(index_t lv=0; lv<nb_vertices(c); ++lv) {
3036 ✗ if(vertex(c,lv) == v) {
3037 ✗ return lv;
3038 }
3039 }
3040 ✗ return NO_VERTEX;
3041 }
3042
3043 /**
3044 * \brief Finds the local index of a facet in a cell
3045 * that can be connected to a facet of another cell
3046 * \param[in] c1 index of the cell
3047 * \param[in] c2 index of the other cell
3048 * \param[in] f2 facet of the other cell
3049 * \return the local index (in 0 .. cell_nb_facets(c1)) of the facet of
3050 * \p c1 that has the same vertices as \p f2 in \p c2 in reverse order,
3051 * modulo a circular permutation, or NO_FACET if such a facet does not
3052 * exist in \p c1.
3053 */
3054 ✗ index_t find_cell_facet(
3055 index_t c1, index_t c2, index_t f2
3056 ) const {
3057 ✗ for(index_t f1=0; f1<nb_facets(c1); ++f1) {
3058 ✗ if(facets_match(c1,f1,c2,f2)) {
3059 ✗ return f1;
3060 }
3061 }
3062 ✗ return NO_FACET;
3063 }
3064
3065 /**
3066 * \brief Tests whether a triangular facet matches a quad facet.
3067 * \details Used to detect non-conformal configurations that should
3068 * be resolved by a connector.
3069 * \param[in] c1 index of the first cell
3070 * \param[in] lf1 index of a triangular facet in \p c1
3071 * \param[in] c2 index of the second cell
3072 * \param[in] lf2 index of a quadrangular facet in \p c2
3073 * \retval true if the three vertices of \p f1 appear in \p f2
3074 * in reverse order
3075 * \retval false otherwise
3076 */
3077 bool triangular_facet_matches_quad_facet(
3078 index_t c1, index_t lf1,
3079 index_t c2, index_t lf2
3080 ) const;
3081
3082
3083 /**
3084 * \brief Tests whether two triangular cell facets have a common edge.
3085 * \param[in] c1 index of the first cell
3086 * \param[in] f1 index of a triangular facet of \p c1
3087 * \param[in] c2 index of the second cell
3088 * \param[in] f2 index of a triangular facet of \p c2
3089 * \param[out] e1 index of the common edge in \p f1
3090 * or NO_EDGE if no such edge exists
3091 * \param[out] e2 index of the common edge in \p f2
3092 * or NO_EDGE if no such edge exists
3093 * \retval true if \p f1 and \p f2 have a common edge
3094 * \retval false otherwise
3095 */
3096 bool triangular_facets_have_common_edge(
3097 index_t c1, index_t f1,
3098 index_t c2, index_t f2,
3099 index_t& e1, index_t& e2
3100 ) const;
3101
3102 /**
3103 * \brief Creates a connector between a quadrandular facet and two
3104 * triangular facets.
3105 * \details This function is used by connect_cells()
3106 * \param[in] c1 index of the cell that has the quadrangular facet
3107 * \param[in] lf1 index of the quadrangular facet in \p c1
3108 * \param[in] matches a const reference to a vector of
3109 * (cell index, facet index) pairs that are candidate triangles to be
3110 * connected to the quadrangular facet. Each of them
3111 * has three vertices in common with the quadrangular facet.
3112 * It may contain more than two (cell,facet) index pairs.
3113 * In this case, among them we select the pair of triangular facets
3114 * that have an edge in common.
3115 * \retval true if a connector was created. A connector is created if
3116 * among the candidate triangular facets there are exactly two facets
3117 * on the border with an edge in common.
3118 * \retval false otherwise
3119 */
3120 bool create_connector(
3121 index_t c1, index_t lf1,
3122 const std::vector< std::pair<index_t, index_t> >& matches
3123 );
3124
3125 /**
3126 * \brief Optimized implementation of connect() used
3127 * when the mesh is simplicial.
3128 */
3129 void connect_tets();
3130
3131 protected:
3132 MeshVertices& vertices_;
3133 MeshCellCornersStore& cell_corners_;
3134 MeshCellFacetsStore& cell_facets_;
3135 friend class Mesh;
3136 friend class GeogramIOHandler;
3137 };
3138
3139 /*************************************************************************/
3140
3141 /**
3142 * \brief Indicates the mesh elements (vertices, facets or cells)
3143 * present in a mesh.
3144 * \details The set of elements present in a mesh is represented
3145 * by a bitwise-or combination of the constants.
3146 * \relates Mesh
3147 */
3148 enum MeshElementsFlags {
3149 MESH_NONE = 0,
3150 MESH_VERTICES = 1,
3151 MESH_FACETS = 2,
3152 MESH_EDGES = 4,
3153 MESH_CELLS = 8,
3154 MESH_ALL_ELEMENTS = 15,
3155 MESH_FACET_CORNERS = 16,
3156 MESH_CELL_CORNERS = 32,
3157 MESH_CELL_FACETS = 64,
3158 MESH_ALL_SUBELEMENTS = 65
3159 };
3160
3161 /*************************************************************************/
3162
3163 /**
3164 * \brief Represents a mesh.
3165 * \details A mesh can have vertices, optionally facets and
3166 * optionally volumetric cells. Attributes can be attached
3167 * to all elements and sub-elements.
3168 */
3169 class GEOGRAM_API Mesh {
3170 public:
3171 MeshVertices vertices;
3172 MeshEdges edges;
3173 MeshFacets facets;
3174 MeshFacetCornersStore facet_corners;
3175 MeshCells cells;
3176 MeshCellCornersStore cell_corners;
3177 MeshCellFacetsStore cell_facets;
3178
3179 /**
3180 * \brief Mesh constructor
3181 * \param[in] dimension dimension of the vertices
3182 * \param[in] single_precision if true, vertices are
3183 * stored in single precision (float), else they are
3184 * stored as double precision (double).
3185 */
3186 Mesh(index_t dimension=3, bool single_precision=false);
3187
3188 /**
3189 * \brief Mesh destructor.
3190 */
3191 virtual ~Mesh();
3192
3193 /**
3194 * \brief Removes all the elements and attributes of
3195 * this mesh.
3196 * \param[in] keep_attributes if true, then all the
3197 * existing attribute names / bindings are kept (but
3198 * they are cleared). If false, they are destroyed.
3199 * \param[in] keep_memory if true, then memory is
3200 * kept and can be reused by subsequent mesh
3201 * element creations.
3202 */
3203 void clear(bool keep_attributes=true, bool keep_memory=false);
3204
3205 /**
3206 * \brief Displays number of vertices, facets and borders.
3207 */
3208
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 27 void show_stats(const std::string& tag = "Mesh") const;
3209
3210
3211 /**
3212 * \brief Does some validity checks.
3213 * \details Used for debugging. If the
3214 * validity checks are not satisfied,
3215 * then it crashes with an assertion
3216 * failure.
3217 */
3218 void assert_is_valid();
3219
3220
3221 /**
3222 * \brief Copies a mesh onto this one
3223 * \details All attributes of this mesh are deleted, and
3224 * Attribute instances connected to this mesh are unbound.
3225 * \param[in] rhs a const reference to the mesh to be copied
3226 * \param[in] copy_attributes if true, all the attributes are
3227 * copied.
3228 * \param[in] what a combination of MESH_VERTICES, MESH_EDGES,
3229 * MESH_FACETS, MESH_CELLS flags. Set to MESH_ALL_ELEMENTS
3230 * to copy everything (default). If MESH_VERTICES is not set,
3231 * then the mesh is cleared.
3232 */
3233 void copy(
3234 const Mesh& rhs,
3235 bool copy_attributes=true,
3236 MeshElementsFlags what=MESH_ALL_ELEMENTS
3237 );
3238
3239 /**
3240 * \brief Loads this mesh from a file.
3241 * \details The file format is deduced from the extension. Supported
3242 * formats are those registered in the mesh I/O handlers
3243 * (.obj, .off, .ply, .stl, .mesh/.meshb, .geogram, ...).
3244 * Equivalent to mesh_load(filename, *this); implemented in
3245 * mesh_io.cpp.
3246 * \param[in] filename the name of the file
3247 * \retval true on success
3248 */
3249 bool load(const std::string& filename);
3250
3251 /**
3252 * \brief Saves this mesh to a file.
3253 * \details The file format is deduced from the extension.
3254 * Equivalent to mesh_save(*this, filename).
3255 * \param[in] filename the name of the file
3256 * \retval true on success
3257 */
3258 bool save(const std::string& filename) const;
3259
3260 /**
3261 * \brief Gets the list of all attributes.
3262 * \return a ';'-separated list of all attributes.
3263 */
3264 std::string get_attributes() const;
3265
3266 /**
3267 * \brief Gets the list of all scalar attributes.
3268 * \return a ';'-separated list of all scalar attributes.
3269 * \details Whenever there is a vector attribute v of dim d,
3270 * it appends v[0];v[1];...v[d-1] to the list.
3271 */
3272 std::string get_scalar_attributes() const;
3273
3274 /**
3275 * \brief Gets the list of all vector attributes.
3276 * \param[in] max_dim if non-zero, only vector attributes of
3277 * dimension lower than \p max_dim are returned.
3278 * \return a ';'-separated list of all vector attributes.
3279 */
3280 std::string get_vector_attributes(index_t max_dim = 0) const;
3281
3282 /**
3283 * \brief Gets the number of subelements types.
3284 * \return the number of subelements types.
3285 */
3286 index_t nb_subelements_types() const;
3287
3288 /**
3289 * \brief Gets a MeshSubElementsStore by index.
3290 * \param[in] i index of the subelements
3291 * \return a reference to the corresponding MeshSubElementsStore
3292 * \pre i < nb_subelements_types()
3293 */
3294 MeshSubElementsStore& get_subelements_by_index(index_t i);
3295
3296 /**
3297 * \brief Gets a MeshSubElementsStore by index.
3298 * \param[in] i index of the subelements
3299 * \return a const reference to the corresponding MeshSubElementsStore
3300 * \pre i < nb_subelements_types()
3301 */
3302 const MeshSubElementsStore& get_subelements_by_index(index_t i) const;
3303
3304
3305 /**
3306 * \brief Gets a MeshSubElementsStore by subelements type.
3307 * \param[in] what one of MESH_VERTICES, MESH_EDGES, MESH_FACETS,
3308 * MESH_FACET_CORNERS, MESH_CELLS, MESH_CELL_CORNERS, MESH_CELL_FACETS
3309 * \return a reference to the corresponding MeshSubElementsStore
3310 */
3311 MeshSubElementsStore& get_subelements_by_type(MeshElementsFlags what);
3312
3313 /**
3314 * \brief Gets a MeshSubElementsStore by subelements type.
3315 * \param[in] what one of MESH_VERTICES, MESH_EDGES, MESH_FACETS,
3316 * MESH_FACET_CORNERS, MESH_CELLS, MESH_CELL_CORNERS, MESH_CELL_FACETS
3317 * \return a const reference to the corresponding MeshSubElementsStore
3318 */
3319 const MeshSubElementsStore& get_subelements_by_type(
3320 MeshElementsFlags what
3321 ) const;
3322
3323 /**
3324 * \brief Gets a subelement name by subelement type.
3325 * \param[in] what one of MESH_VERTICES, MESH_EDGES, MESH_FACETS,
3326 * MESH_FACET_CORNERS, MESH_CELLS, MESH_CELL_CORNERS, MESH_CELL_FACETS
3327 * \return a string with the name of the subelement.
3328 */
3329 static std::string subelements_type_to_name(MeshElementsFlags what);
3330
3331 /**
3332 * \brief Gets a subelement type by subelement name.
3333 * \param[in] name the name of the subelement as a string
3334 * \return one of MESH_VERTICES, MESH_EDGES, MESH_FACETS,
3335 * MESH_FACET_CORNERS, MESH_CELLS, MESH_CELL_CORNERS, MESH_CELL_FACETS
3336 * or MESH_NONE if the name is invalid
3337 */
3338 static MeshElementsFlags name_to_subelements_type(
3339 const std::string& name
3340 );
3341
3342 /**
3343 * \brief Extracts localisation, name and optional component from
3344 * an attribute name.
3345 * \param[in] full_attribute_name for instance, facets.density, or
3346 * vertices.normal[0]
3347 * \param[out] where one of MESH_VERTICES, MESH_EDGES, MESH_FACETS,
3348 * MESH_FACET_CORNERS, MESH_CELLS, MESH_CELL_FACETS, MESH_CELL_CORNERS
3349 * \param[out] attribute_name the name of the attribute, without the
3350 * localisation and without the component
3351 * \param[out] component the component (between square brackets in
3352 * \p full_attribute_name) or 0 if no component was specified
3353 * \retval true if the attribute name could be parsed
3354 * \retval false if the attribute name has invalid syntax
3355 */
3356 static bool parse_attribute_name(
3357 const std::string& full_attribute_name,
3358 MeshElementsFlags& where,
3359 std::string& attribute_name,
3360 index_t& component
3361 );
3362
3363 protected:
3364 /**
3365 * \brief Displays the list of attributes to the Logger.
3366 * \param[in] tag the tag to be sent to the Logger
3367 * \param[in] subelement_name the name of the subelement
3368 * (vertices, facets, facet_corners ...)
3369 * \param[in] subelements a const reference to the MeshSubElementsStore
3370 */
3371 void display_attributes(
3372 const std::string& tag, const std::string& subelement_name,
3373 const MeshSubElementsStore& subelements
3374 ) const;
3375
3376 /**
3377 * \brief Forbids copy.
3378 * \details This is to make sure that client code does
3379 * not unintentionlly copies a Mesh (for
3380 * instance by passing it by-value to a function).
3381 * Use copy() instead.
3382 */
3383 Mesh(const Mesh& rhs) = delete;
3384
3385 /**
3386 * \brief Forbids copy.
3387 * \details This is to make sure that client code does
3388 * not unintentionlly copies a Mesh (for
3389 * instance by passing it by-value to a function).
3390 * Use copy() instead.
3391 */
3392 const Mesh& operator=(const Mesh& rhs) = delete;
3393 };
3394
3395 /*************************************************************************/
3396 }
3397
3398 #endif
3399