| 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 | #ifdef GEOGRAM_WITH_PDEL | ||
| 41 | |||
| 42 | #include <geogram/delaunay/parallel_delaunay_3d.h> | ||
| 43 | #include <geogram/delaunay/cavity.h> | ||
| 44 | #include <geogram/mesh/mesh_reorder.h> | ||
| 45 | #include <geogram/numerics/predicates.h> | ||
| 46 | #include <geogram/basic/geometry.h> | ||
| 47 | #include <geogram/basic/stopwatch.h> | ||
| 48 | #include <geogram/basic/command_line.h> | ||
| 49 | #include <geogram/basic/permutation.h> | ||
| 50 | #include <geogram/bibliography/bibliography.h> | ||
| 51 | |||
| 52 | #include <mutex> | ||
| 53 | #include <condition_variable> | ||
| 54 | |||
| 55 | // ParallelDelaunayThread class, declared locally, has | ||
| 56 | // no out-of-line virtual functions. It is not a | ||
| 57 | // problem since they are only visible from this translation | ||
| 58 | // unit, but clang will complain. | ||
| 59 | #ifdef __clang__ | ||
| 60 | #pragma GCC diagnostic ignored "-Wweak-vtables" | ||
| 61 | #endif | ||
| 62 | |||
| 63 | namespace { | ||
| 64 | using namespace GEO; | ||
| 65 | |||
| 66 | /** | ||
| 67 | * \brief Generates a random integer. | ||
| 68 | * \return a random integer between 0 and \p choices - 1 | ||
| 69 | * \param [in] choices_in number of possible choices for the | ||
| 70 | * random variable (maximum value + 1) | ||
| 71 | * \details The function is thread-safe, and uses one seed | ||
| 72 | * per thread. | ||
| 73 | */ | ||
| 74 | 5 | index_t thread_safe_random(index_t choices_in) { | |
| 75 | #ifdef GARGANTUA | ||
| 76 | typedef Numeric::int64 Int; | ||
| 77 | #else | ||
| 78 | typedef long int Int; | ||
| 79 | #endif | ||
| 80 | 5 | signed_index_t choices = signed_index_t(choices_in); | |
| 81 | static thread_local Int randomseed = 1l ; | ||
| 82 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
|
5 | if (choices >= 714025l) { |
| 83 | ✗ | Int newrandom = (randomseed * 1366l + 150889l) % 714025l; | |
| 84 | ✗ | randomseed = (newrandom * 1366l + 150889l) % 714025l; | |
| 85 | ✗ | newrandom = newrandom * (choices / 714025l) + randomseed; | |
| 86 | ✗ | if (newrandom >= choices) { | |
| 87 | ✗ | return index_t(newrandom - choices); | |
| 88 | } else { | ||
| 89 | ✗ | return index_t(newrandom); | |
| 90 | } | ||
| 91 | } else { | ||
| 92 | 5 | randomseed = (randomseed * 1366l + 150889l) % 714025l; | |
| 93 | 5 | return index_t(randomseed % choices); | |
| 94 | } | ||
| 95 | } | ||
| 96 | |||
| 97 | /** | ||
| 98 | * \brief Generates a random integer between 0 and 3. | ||
| 99 | * \return a random integer between 0 and 3 | ||
| 100 | * \details The function is thread-safe, and uses one seed | ||
| 101 | * per thread. | ||
| 102 | */ | ||
| 103 | 4119 | index_t thread_safe_random_4() { | |
| 104 | static thread_local long int randomseed = 1l ; | ||
| 105 | 4119 | randomseed = (randomseed * 1366l + 150889l) % 714025l; | |
| 106 | 4119 | return index_t(randomseed % 4); | |
| 107 | } | ||
| 108 | } | ||
| 109 | |||
| 110 | namespace GEO { | ||
| 111 | |||
| 112 | /** | ||
| 113 | * \brief One of the threads of the multi-threaded | ||
| 114 | * 3d Delaunay implementation. | ||
| 115 | */ | ||
| 116 | class Delaunay3dThread : public GEO::Thread { | ||
| 117 | public: | ||
| 118 | |||
| 119 | /** | ||
| 120 | * \brief Symbolic value for cell_status_[t] that | ||
| 121 | * indicates that no thread owns t. | ||
| 122 | */ | ||
| 123 | static constexpr index_t NO_THREAD = CellStatusArray::FREE_CELL; | ||
| 124 | |||
| 125 | /** | ||
| 126 | * \brief Creates a new Delaunay3dThread. | ||
| 127 | * \details Each Delaunay3dThread has an affected working | ||
| 128 | * zone, i.e. a range of tetrahedra indices in which the | ||
| 129 | * thread is allowed to create tetrahedra. | ||
| 130 | * \param[in] master a pointer to the ParallelDelaunay3d | ||
| 131 | * this thread belongs to | ||
| 132 | * \param[in] pool_begin first tetrahedron index of | ||
| 133 | * the working zone of this Delaunay3dThread | ||
| 134 | * \param[in] pool_end one position past the last tetrahedron | ||
| 135 | * index of the working zone of this Delaunay3dThread | ||
| 136 | */ | ||
| 137 | 20 | Delaunay3dThread( | |
| 138 | ParallelDelaunay3d* master, | ||
| 139 | index_t pool_begin, | ||
| 140 | index_t pool_end | ||
| 141 | 20 | ) : | |
| 142 | 20 | master_(master), | |
| 143 | 20 | cell_to_v_store_(master_->cell_to_v_store_), | |
| 144 | 20 | cell_to_cell_store_(master_->cell_to_cell_store_), | |
| 145 | 20 | cell_next_(master_->cell_next_), | |
| 146 |
1/2✓ Branch 8 taken 20 times.
✗ Branch 9 not taken.
|
20 | cell_status_(master_->cell_status_) |
| 147 | { | ||
| 148 | |||
| 149 | // max_used_t_ is initialized to 1 so that | ||
| 150 | // computing modulos does not trigger FPEs | ||
| 151 | // at the beginning. | ||
| 152 | 20 | max_used_t_ = 1; | |
| 153 | 20 | max_t_ = master_->cell_next_.size(); | |
| 154 | |||
| 155 | 20 | nb_vertices_ = master_->nb_vertices(); | |
| 156 |
1/2✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
|
20 | vertices_ = master_->vertex_ptr(0); |
| 157 | 20 | weighted_ = master_->weighted_; | |
| 158 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 20 times.
|
20 | heights_ = weighted_ ? master_->heights_.data() : nullptr; |
| 159 | 20 | dimension_ = master_->dimension(); | |
| 160 | 20 | vertex_stride_ = dimension_; | |
| 161 | 20 | reorder_ = master_->reorder_.data(); | |
| 162 | |||
| 163 | // Initialize free list in memory pool | ||
| 164 | 20 | first_free_ = pool_begin; | |
| 165 |
2/2✓ Branch 0 taken 28386 times.
✓ Branch 1 taken 20 times.
|
28406 | for(index_t t=pool_begin; t<pool_end-1; ++t) { |
| 166 |
1/2✓ Branch 1 taken 28386 times.
✗ Branch 2 not taken.
|
28386 | cell_next_[t] = t+1; |
| 167 | } | ||
| 168 |
1/2✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
|
20 | cell_next_[pool_end-1] = END_OF_LIST; |
| 169 | 20 | nb_free_ = pool_end - pool_begin; | |
| 170 | 20 | memory_overflow_ = false; | |
| 171 | |||
| 172 | 20 | work_begin_ = NO_INDEX; | |
| 173 | 20 | work_end_ = NO_INDEX; | |
| 174 | 20 | finished_ = false; | |
| 175 | 20 | b_hint_ = NO_TETRAHEDRON; | |
| 176 | 20 | e_hint_ = NO_TETRAHEDRON; | |
| 177 | 20 | direction_ = true; | |
| 178 | |||
| 179 | #ifdef GEO_DEBUG | ||
| 180 | 20 | nb_acquired_tets_ = 0; | |
| 181 | #endif | ||
| 182 | 20 | interfering_thread_ = NO_THREAD; | |
| 183 | |||
| 184 | 20 | nb_rollbacks_ = 0; | |
| 185 | 20 | nb_failed_locate_ = 0; | |
| 186 | |||
| 187 | 20 | nb_tets_to_create_ = 0; | |
| 188 | 20 | t_boundary_ = NO_TETRAHEDRON; | |
| 189 | 20 | f_boundary_ = NO_INDEX; | |
| 190 | |||
| 191 | 20 | v1_ = NO_INDEX; | |
| 192 | 20 | v2_ = NO_INDEX; | |
| 193 | 20 | v3_ = NO_INDEX; | |
| 194 | 20 | v4_ = NO_INDEX; | |
| 195 | 20 | } | |
| 196 | |||
| 197 | /** | ||
| 198 | * \brief Copies some variables from another thread. | ||
| 199 | * \param[in] rhs the thread from which variables should | ||
| 200 | * be copied | ||
| 201 | * \details copies v1_, v2_, v3_, v4_ (indices of the vertices | ||
| 202 | * of the first created tetrahedron), max_used_t_ (maximum | ||
| 203 | * used tetrahedron index) and max_t_ (maximum valid tetrahedron | ||
| 204 | * index). | ||
| 205 | */ | ||
| 206 | 15 | void initialize_from(const Delaunay3dThread* rhs) { | |
| 207 | 15 | max_used_t_ = rhs->max_used_t_; | |
| 208 | 15 | max_t_ = rhs->max_t_; | |
| 209 | 15 | v1_ = rhs->v1_; | |
| 210 | 15 | v2_ = rhs->v2_; | |
| 211 | 15 | v3_ = rhs->v3_; | |
| 212 | 15 | v4_ = rhs->v4_; | |
| 213 | 15 | } | |
| 214 | |||
| 215 | /** | ||
| 216 | * \brief Gets the number of rollbacks. | ||
| 217 | * \return the number of rollbacks | ||
| 218 | * \details rollbacks occur whenever a point | ||
| 219 | * could not be inserted, due to interferences | ||
| 220 | * from other threads | ||
| 221 | */ | ||
| 222 | 32 | index_t nb_rollbacks() const { | |
| 223 | 32 | return nb_rollbacks_; | |
| 224 | } | ||
| 225 | |||
| 226 | /** | ||
| 227 | * \brief Gets the number of failed locate() calls. | ||
| 228 | * \return the number of failed locate() calls | ||
| 229 | * \details locate() can fail when it cannot acquire | ||
| 230 | * the tetrahedra that are traversed, due to | ||
| 231 | * interferences from another thread. | ||
| 232 | */ | ||
| 233 | 32 | index_t nb_failed_locate() const { | |
| 234 | 32 | return nb_failed_locate_; | |
| 235 | } | ||
| 236 | |||
| 237 | /** | ||
| 238 | * \brief Sets the point index sequence that | ||
| 239 | * should be processed by this thread. | ||
| 240 | * \param[in] b index of the first point to insert | ||
| 241 | * \param[in] e one position past the index of the | ||
| 242 | * last point to insert | ||
| 243 | */ | ||
| 244 | 5 | void set_work(index_t b, index_t e) { | |
| 245 | 5 | work_begin_ = b; | |
| 246 | // e is one position past the last point index | ||
| 247 | // to insert. | ||
| 248 | 5 | work_end_ = e-1; | |
| 249 | 5 | } | |
| 250 | |||
| 251 | /** | ||
| 252 | * \brief Gets the number of remaining points to | ||
| 253 | * be inserted. | ||
| 254 | * \return the number of points to be inserted by | ||
| 255 | * this thread | ||
| 256 | */ | ||
| 257 | 20 | index_t work_size() const { | |
| 258 |
3/4✓ Branch 0 taken 15 times.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 15 times.
✗ Branch 3 not taken.
|
20 | if(work_begin_ == NO_INDEX && work_end_ == NO_INDEX) { |
| 259 | 15 | return 0; | |
| 260 | } | ||
| 261 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
5 | geo_debug_assert(work_begin_ != NO_INDEX); |
| 262 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
5 | geo_debug_assert(work_end_ != NO_INDEX); |
| 263 | 5 | return std::max(work_end_ - work_begin_ + 1, index_t(0)); | |
| 264 | } | ||
| 265 | |||
| 266 | /** | ||
| 267 | * \brief Gets the number of threads. | ||
| 268 | * \return the number of threads created by | ||
| 269 | * the master ParallelDelaunay3d of this thread. | ||
| 270 | */ | ||
| 271 | index_t nb_threads() const { | ||
| 272 | return index_t(master_->threads_.size()); | ||
| 273 | } | ||
| 274 | |||
| 275 | /** | ||
| 276 | * \brief Gets a thread by index | ||
| 277 | * \pre t < nb_threads() | ||
| 278 | * \param[in] t index of the thread | ||
| 279 | * \return a poiner to the \p t th thread | ||
| 280 | */ | ||
| 281 | ✗ | Delaunay3dThread* thread(index_t t) { | |
| 282 | ✗ | return static_cast<Delaunay3dThread*>(master_->threads_[t].get()); | |
| 283 | } | ||
| 284 | |||
| 285 | /** | ||
| 286 | * \brief Inserts the point sequence allocated to | ||
| 287 | * this thread. | ||
| 288 | * \details The point sequence was previously defined | ||
| 289 | * by set_work(). | ||
| 290 | */ | ||
| 291 | 25 | void run() override { | |
| 292 | |||
| 293 | 25 | finished_ = false; | |
| 294 | |||
| 295 |
3/4✓ Branch 0 taken 10 times.
✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 10 times.
|
25 | if(work_begin_ == NO_INDEX || work_end_ == NO_INDEX) { |
| 296 | 15 | return ; | |
| 297 | } | ||
| 298 | |||
| 299 | 10 | memory_overflow_ = false; | |
| 300 | |||
| 301 | // Current hint associated with b | ||
| 302 | 10 | b_hint_ = NO_TETRAHEDRON; | |
| 303 | |||
| 304 | // Current hint associated with e | ||
| 305 | 10 | e_hint_ = NO_TETRAHEDRON; | |
| 306 | |||
| 307 | // If true, insert in b->e order, | ||
| 308 | // else insert in e->b order | ||
| 309 | 10 | direction_ = true; | |
| 310 | |||
| 311 |
3/4✓ Branch 0 taken 4058 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 4058 times.
✗ Branch 3 not taken.
|
4068 | while(work_end_ >= work_begin_ && !memory_overflow_) { |
| 312 |
1/2✓ Branch 0 taken 4058 times.
✗ Branch 1 not taken.
|
4058 | index_t v = direction_ ? work_begin_ : work_end_ ; |
| 313 |
1/2✓ Branch 0 taken 4058 times.
✗ Branch 1 not taken.
|
4058 | index_t& hint = direction_ ? b_hint_ : e_hint_; |
| 314 | |||
| 315 | // Try to insert v and update hint | ||
| 316 | 4058 | bool success = insert(reorder_[v],hint); | |
| 317 | |||
| 318 | // Notify all threads that are waiting for | ||
| 319 | // this thread to release some tetrahedra. | ||
| 320 | 4058 | send_event(); | |
| 321 | |||
| 322 |
1/2✓ Branch 0 taken 4058 times.
✗ Branch 1 not taken.
|
4058 | if(success) { |
| 323 |
1/2✓ Branch 0 taken 4058 times.
✗ Branch 1 not taken.
|
4058 | if(direction_) { |
| 324 | 4058 | ++work_begin_; | |
| 325 | } else { | ||
| 326 | ✗ | --work_end_; | |
| 327 | } | ||
| 328 | } else { | ||
| 329 | ✗ | ++nb_rollbacks_; | |
| 330 | ✗ | if(interfering_thread_ != NO_THREAD) { | |
| 331 | ✗ | if(id() < interfering_thread_) { | |
| 332 | // If this thread has a higher priority than | ||
| 333 | // the one that interfered, wait for the | ||
| 334 | // interfering thread to release the tets that | ||
| 335 | // it holds (then the loop will retry to insert | ||
| 336 | // the same vertex). | ||
| 337 | ✗ | wait_for_event(interfering_thread_); | |
| 338 | } else { | ||
| 339 | // If this thread has a lower priority than | ||
| 340 | // the interfering thread, try inserting | ||
| 341 | // from the other end of the points sequence. | ||
| 342 | ✗ | direction_ = !direction_; | |
| 343 | } | ||
| 344 | } | ||
| 345 | } | ||
| 346 | } | ||
| 347 | 10 | finished_ = true; | |
| 348 | |||
| 349 | // Fix by Hiep Vu: wake up threads that potentially missed | ||
| 350 | // the previous wake ups. | ||
| 351 | 10 | mutex_.lock(); | |
| 352 | 10 | send_event(); | |
| 353 | 10 | mutex_.unlock(); | |
| 354 | } | ||
| 355 | |||
| 356 | /** | ||
| 357 | * \brief Symbolic constant for uninitialized hint. | ||
| 358 | * \details Locate functions can be accelerated by | ||
| 359 | * specifying a hint. This constant indicates that | ||
| 360 | * no hint is given. | ||
| 361 | */ | ||
| 362 | static constexpr index_t NO_TETRAHEDRON = NO_INDEX; | ||
| 363 | |||
| 364 | /** | ||
| 365 | * \brief Symbolic value for a vertex of a | ||
| 366 | * tetrahedron that indicates a virtual tetrahedron. | ||
| 367 | * \details The three other vertices then correspond to a | ||
| 368 | * facet on the convex hull of the points. | ||
| 369 | */ | ||
| 370 | static constexpr index_t VERTEX_AT_INFINITY = NO_INDEX; | ||
| 371 | |||
| 372 | |||
| 373 | /** | ||
| 374 | * \brief Maximum valid index for a tetrahedron. | ||
| 375 | * \details This includes not only real tetrahedra, | ||
| 376 | * but also the virtual ones on the border, the conflict | ||
| 377 | * list and the free list. | ||
| 378 | * \return the maximum valid index for a tetrahedron | ||
| 379 | */ | ||
| 380 | 12047424 | index_t max_t() const { | |
| 381 | 12047424 | return max_t_; | |
| 382 | } | ||
| 383 | |||
| 384 | /** | ||
| 385 | * \brief Tests whether a given tetrahedron | ||
| 386 | * is a finite one. | ||
| 387 | * \details Infinite tetrahedra are the ones | ||
| 388 | * that are incident to the infinite vertex | ||
| 389 | * (index -1) | ||
| 390 | * \param[in] t the index of the tetrahedron | ||
| 391 | * \retval true if \p t is finite | ||
| 392 | * \retval false otherwise | ||
| 393 | */ | ||
| 394 | 32526 | bool tet_is_finite(index_t t) const { | |
| 395 | return | ||
| 396 | 32526 | cell_to_v_store_[4 * t] != NO_INDEX && | |
| 397 |
2/2✓ Branch 1 taken 31603 times.
✓ Branch 2 taken 921 times.
|
32524 | cell_to_v_store_[4 * t + 1] != NO_INDEX && |
| 398 |
4/4✓ Branch 0 taken 32524 times.
✓ Branch 1 taken 2 times.
✓ Branch 3 taken 31109 times.
✓ Branch 4 taken 494 times.
|
96159 | cell_to_v_store_[4 * t + 2] != NO_INDEX && |
| 399 |
2/2✓ Branch 1 taken 30307 times.
✓ Branch 2 taken 802 times.
|
63635 | cell_to_v_store_[4 * t + 3] != NO_INDEX ; |
| 400 | } | ||
| 401 | |||
| 402 | /** | ||
| 403 | * \brief Tests whether a tetrahedron is | ||
| 404 | * a real one. | ||
| 405 | * \details Real tetrahedra are incident to | ||
| 406 | * four user-specified vertices (there are also | ||
| 407 | * virtual tetrahedra that are incident to the | ||
| 408 | * vertex at infinity, with index -1) | ||
| 409 | * \param[in] t index of the tetrahedron | ||
| 410 | * \retval true if tetrahedron \p t is a real one | ||
| 411 | * \retval false otherwise | ||
| 412 | */ | ||
| 413 | 29784 | bool tet_is_real(index_t t) const { | |
| 414 |
4/4✓ Branch 1 taken 8371 times.
✓ Branch 2 taken 21413 times.
✓ Branch 4 taken 8005 times.
✓ Branch 5 taken 366 times.
|
29784 | return !tet_is_free(t) && tet_is_finite(t); |
| 415 | } | ||
| 416 | |||
| 417 | /** | ||
| 418 | * \brief Tests whether a tetrahedron is | ||
| 419 | * in the free list. | ||
| 420 | * \details Deleted tetrahedra are recycled | ||
| 421 | * in a free list. | ||
| 422 | * \param[in] t index of the tetrahedron | ||
| 423 | * \retval true if tetrahedron \p t is in | ||
| 424 | * the free list | ||
| 425 | * \retval false otherwise | ||
| 426 | */ | ||
| 427 | 180020 | bool tet_is_free(index_t t) const { | |
| 428 | 180020 | return tet_is_in_list(t); | |
| 429 | } | ||
| 430 | |||
| 431 | /** | ||
| 432 | * \brief Tests whether a tetrahedron is contained | ||
| 433 | * by a given linked list. | ||
| 434 | * \details Used for debugging purposes. | ||
| 435 | * \param[in] t the tetrahedron | ||
| 436 | * \param[in] first the first element of the list | ||
| 437 | * \retval true if \p t is contained in the list starting | ||
| 438 | * at \p first | ||
| 439 | * \retval false otherwise | ||
| 440 | */ | ||
| 441 | bool tet_is_in_list(index_t t, index_t first) const { | ||
| 442 | for( | ||
| 443 | index_t cur = first; cur != END_OF_LIST; | ||
| 444 | cur = tet_next(cur) | ||
| 445 | ) { | ||
| 446 | if(cur == t) { | ||
| 447 | return true; | ||
| 448 | } | ||
| 449 | } | ||
| 450 | return false; | ||
| 451 | } | ||
| 452 | |||
| 453 | |||
| 454 | /** | ||
| 455 | * \brief Finds in the pointset a set of four non-coplanar | ||
| 456 | * points and creates a tetrahedron that connects them. | ||
| 457 | * \details This function is used to initiate the incremental | ||
| 458 | * Delaunay construction, it should be called only once. | ||
| 459 | * \retval the index of the created tetrahedron | ||
| 460 | * \retval NO_TETRAHEDRON if all points were coplanar | ||
| 461 | */ | ||
| 462 | 5 | index_t create_first_tetrahedron() { | |
| 463 | index_t iv0,iv1,iv2,iv3; | ||
| 464 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 5 times.
|
5 | if(nb_vertices() < 4) { |
| 465 | ✗ | return NO_TETRAHEDRON; | |
| 466 | } | ||
| 467 | |||
| 468 | 5 | iv0 = 0; | |
| 469 | |||
| 470 | 5 | iv1 = 1; | |
| 471 | 5 | while( | |
| 472 |
2/4✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 5 times.
|
10 | iv1 < nb_vertices() && |
| 473 |
4/8✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 5 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 5 times.
|
5 | PCK::points_are_identical_3d( |
| 474 | vertex_ptr(iv0), vertex_ptr(iv1) | ||
| 475 | ) | ||
| 476 | ) { | ||
| 477 | ✗ | ++iv1; | |
| 478 | } | ||
| 479 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 5 times.
|
5 | if(iv1 == nb_vertices()) { |
| 480 | ✗ | return NO_TETRAHEDRON; | |
| 481 | } | ||
| 482 | |||
| 483 | 5 | iv2 = iv1 + 1; | |
| 484 | 5 | while( | |
| 485 |
2/4✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 5 times.
|
10 | iv2 < nb_vertices() && |
| 486 |
5/10✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 5 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 5 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 5 times.
|
5 | PCK::points_are_colinear_3d( |
| 487 | vertex_ptr(iv0), vertex_ptr(iv1), vertex_ptr(iv2) | ||
| 488 | ) | ||
| 489 | ) { | ||
| 490 | ✗ | ++iv2; | |
| 491 | } | ||
| 492 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 5 times.
|
5 | if(iv2 == nb_vertices()) { |
| 493 | ✗ | return NO_TETRAHEDRON; | |
| 494 | } | ||
| 495 | |||
| 496 | 5 | iv3 = iv2 + 1; | |
| 497 | 5 | Sign s = ZERO; | |
| 498 | 5 | while( | |
| 499 |
3/4✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 31 times.
✓ Branch 4 taken 5 times.
|
72 | iv3 < nb_vertices() && |
| 500 |
7/12✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 36 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 36 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 36 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 31 times.
✓ Branch 16 taken 5 times.
|
36 | (s = PCK::orient_3d( |
| 501 | vertex_ptr(iv0), vertex_ptr(iv1), | ||
| 502 | vertex_ptr(iv2), vertex_ptr(iv3) | ||
| 503 | )) == ZERO | ||
| 504 | ) { | ||
| 505 | 31 | ++iv3; | |
| 506 | } | ||
| 507 | |||
| 508 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 5 times.
|
5 | if(iv3 == nb_vertices()) { |
| 509 | ✗ | return NO_TETRAHEDRON; | |
| 510 | } | ||
| 511 | |||
| 512 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
5 | geo_debug_assert(s != ZERO); |
| 513 | |||
| 514 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 3 times.
|
5 | if(s == NEGATIVE) { |
| 515 | 2 | std::swap(iv2, iv3); | |
| 516 | } | ||
| 517 | |||
| 518 | // Create the first tetrahedron | ||
| 519 |
1/2✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
|
5 | index_t t0 = new_tetrahedron(iv0, iv1, iv2, iv3); |
| 520 | |||
| 521 | // Create the first four virtual tetrahedra surrounding it | ||
| 522 | index_t t[4]; | ||
| 523 |
2/2✓ Branch 0 taken 20 times.
✓ Branch 1 taken 5 times.
|
25 | for(index_t f = 0; f < 4; ++f) { |
| 524 | // In reverse order since it is an adjacent tetrahedron | ||
| 525 |
2/4✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 20 times.
✗ Branch 5 not taken.
|
20 | index_t v1 = tet_vertex(t0, tet_facet_vertex(f,2)); |
| 526 |
2/4✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 20 times.
✗ Branch 5 not taken.
|
20 | index_t v2 = tet_vertex(t0, tet_facet_vertex(f,1)); |
| 527 |
2/4✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 20 times.
✗ Branch 5 not taken.
|
20 | index_t v3 = tet_vertex(t0, tet_facet_vertex(f,0)); |
| 528 |
1/2✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
|
20 | t[f] = new_tetrahedron(VERTEX_AT_INFINITY, v1, v2, v3); |
| 529 | } | ||
| 530 | |||
| 531 | // Connect the virtual tetrahedra to the real one | ||
| 532 |
2/2✓ Branch 0 taken 20 times.
✓ Branch 1 taken 5 times.
|
25 | for(index_t f=0; f<4; ++f) { |
| 533 |
1/2✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
|
20 | set_tet_adjacent(t[f], 0, t0); |
| 534 |
1/2✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
|
20 | set_tet_adjacent(t0, f, t[f]); |
| 535 | } | ||
| 536 | |||
| 537 | // Interconnect the four virtual tetrahedra along their common | ||
| 538 | // faces | ||
| 539 |
2/2✓ Branch 0 taken 20 times.
✓ Branch 1 taken 5 times.
|
25 | for(index_t f = 0; f < 4; ++f) { |
| 540 | // In reverse order since it is an adjacent tetrahedron | ||
| 541 |
1/2✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
|
20 | index_t lv1 = tet_facet_vertex(f,2); |
| 542 |
1/2✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
|
20 | index_t lv2 = tet_facet_vertex(f,1); |
| 543 |
1/2✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
|
20 | index_t lv3 = tet_facet_vertex(f,0); |
| 544 |
1/2✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
|
20 | set_tet_adjacent(t[f], 1, t[lv1]); |
| 545 |
1/2✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
|
20 | set_tet_adjacent(t[f], 2, t[lv2]); |
| 546 |
1/2✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
|
20 | set_tet_adjacent(t[f], 3, t[lv3]); |
| 547 | } | ||
| 548 | |||
| 549 | 5 | v1_ = iv0; | |
| 550 | 5 | v2_ = iv1; | |
| 551 | 5 | v3_ = iv2; | |
| 552 | 5 | v4_ = iv3; | |
| 553 | |||
| 554 |
1/2✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
|
5 | release_tets(); |
| 555 | |||
| 556 | 5 | return t0; | |
| 557 | } | ||
| 558 | |||
| 559 | |||
| 560 | /** | ||
| 561 | * \brief Creates a star of tetrahedra filling the conflict | ||
| 562 | * zone. | ||
| 563 | * \param[in] v the index of the point to be inserted | ||
| 564 | * \details This function is used when the Cavity computed | ||
| 565 | * when traversing the conflict zone is OK, that is to say | ||
| 566 | * when its array sizes were not exceeded. | ||
| 567 | * \return the index of one the newly created tetrahedron | ||
| 568 | */ | ||
| 569 | 4032 | index_t stellate_cavity(index_t v) { | |
| 570 | 4032 | index_t new_tet = NO_INDEX; | |
| 571 | |||
| 572 |
2/2✓ Branch 1 taken 101866 times.
✓ Branch 2 taken 4032 times.
|
105898 | for(index_t f=0; f<cavity_.nb_facets(); ++f) { |
| 573 | 101866 | index_t old_tet = cavity_.facet_tet(f); | |
| 574 | 101866 | index_t lf = cavity_.facet_facet(f); | |
| 575 | 101866 | index_t t_neigh = tet_adjacent(old_tet, lf); | |
| 576 | 101866 | index_t v1 = cavity_.facet_vertex(f,0); | |
| 577 | 101866 | index_t v2 = cavity_.facet_vertex(f,1); | |
| 578 | 101866 | index_t v3 = cavity_.facet_vertex(f,2); | |
| 579 | 101866 | new_tet = new_tetrahedron(v, v1, v2, v3); | |
| 580 | 101866 | set_tet_adjacent(new_tet, 0, t_neigh); | |
| 581 | 101866 | set_tet_adjacent( | |
| 582 | t_neigh, find_tet_adjacent(t_neigh,old_tet), new_tet | ||
| 583 | ); | ||
| 584 | 101866 | cavity_.set_facet_tet(f, new_tet); | |
| 585 | } | ||
| 586 | |||
| 587 |
2/2✓ Branch 1 taken 101866 times.
✓ Branch 2 taken 4032 times.
|
105898 | for(index_t f=0; f<cavity_.nb_facets(); ++f) { |
| 588 |
1/2✓ Branch 1 taken 101866 times.
✗ Branch 2 not taken.
|
101866 | new_tet = cavity_.facet_tet(f); |
| 589 | index_t neigh1, neigh2, neigh3; | ||
| 590 |
1/2✓ Branch 1 taken 101866 times.
✗ Branch 2 not taken.
|
101866 | cavity_.get_facet_neighbor_tets(f, neigh1, neigh2, neigh3); |
| 591 |
1/2✓ Branch 1 taken 101866 times.
✗ Branch 2 not taken.
|
101866 | set_tet_adjacent(new_tet, 1, neigh1); |
| 592 |
1/2✓ Branch 1 taken 101866 times.
✗ Branch 2 not taken.
|
101866 | set_tet_adjacent(new_tet, 2, neigh2); |
| 593 |
1/2✓ Branch 1 taken 101866 times.
✗ Branch 2 not taken.
|
101866 | set_tet_adjacent(new_tet, 3, neigh3); |
| 594 | } | ||
| 595 | |||
| 596 | 4032 | return new_tet; | |
| 597 | } | ||
| 598 | |||
| 599 | /** | ||
| 600 | * \brief Inserts a point in the triangulation. | ||
| 601 | * \param[in] v the index of the point to be inserted | ||
| 602 | * \param[in,out] hint the index of a tetrahedron as near as | ||
| 603 | * possible to \p v, or NO_TETRAHEDRON if unspecified. On | ||
| 604 | * exit, the index of one of the tetrahedra incident to | ||
| 605 | * point \p v | ||
| 606 | * \retval true if insertion was successful | ||
| 607 | * \retval false otherwise | ||
| 608 | */ | ||
| 609 | 4058 | bool insert(index_t v, index_t& hint) { | |
| 610 | |||
| 611 | // If v is one of the vertices of the | ||
| 612 | // first tetrahedron, nothing to do. | ||
| 613 | 4058 | if( | |
| 614 |
2/2✓ Branch 0 taken 4053 times.
✓ Branch 1 taken 5 times.
|
4058 | v == v1_ || |
| 615 |
2/2✓ Branch 0 taken 4048 times.
✓ Branch 1 taken 5 times.
|
4053 | v == v2_ || |
| 616 |
2/2✓ Branch 0 taken 4043 times.
✓ Branch 1 taken 5 times.
|
4048 | v == v3_ || |
| 617 |
2/2✓ Branch 0 taken 5 times.
✓ Branch 1 taken 4038 times.
|
4043 | v == v4_ |
| 618 | ) { | ||
| 619 | 20 | return true; | |
| 620 | } | ||
| 621 | |||
| 622 | Sign orient[4]; | ||
| 623 |
2/4✓ Branch 1 taken 4038 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4038 times.
✗ Branch 5 not taken.
|
4038 | index_t t = locate(vertex_ptr(v),hint,orient); |
| 624 | |||
| 625 | // locate() may fail due to tets already owned by | ||
| 626 | // other threads. | ||
| 627 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4038 times.
|
4038 | if(t == NO_TETRAHEDRON) { |
| 628 | ✗ | ++nb_failed_locate_; | |
| 629 | ✗ | geo_debug_assert(nb_acquired_tets_ == 0); | |
| 630 | ✗ | return false; | |
| 631 | } | ||
| 632 | |||
| 633 | // At this point, t is a valid tetrahedron, | ||
| 634 | // and this thread acquired a lock on it. | ||
| 635 | |||
| 636 | // Test whether the point already exists in | ||
| 637 | // the triangulation. The point already exists | ||
| 638 | // if it's located on three faces of the | ||
| 639 | // tetrahedron returned by locate(). | ||
| 640 | 4038 | int nb_zero = | |
| 641 | 4038 | (orient[0] == ZERO) + | |
| 642 | 4038 | (orient[1] == ZERO) + | |
| 643 | 4038 | (orient[2] == ZERO) + | |
| 644 | 4038 | (orient[3] == ZERO) ; | |
| 645 | |||
| 646 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4038 times.
|
4038 | if(nb_zero >= 3) { |
| 647 | ✗ | release_tet(t); | |
| 648 | ✗ | return true; | |
| 649 | } | ||
| 650 | |||
| 651 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 4038 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
4038 | geo_debug_assert(nb_acquired_tets_ == 1); |
| 652 |
5/14✓ Branch 0 taken 4038 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 4038 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 4038 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 4038 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 4038 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
|
4038 | geo_debug_assert(weighted_ || tet_is_in_conflict(t,vertex_ptr(v))); |
| 653 | |||
| 654 | 4038 | index_t t_bndry = NO_TETRAHEDRON; | |
| 655 | 4038 | index_t f_bndry = NO_INDEX; | |
| 656 | |||
| 657 | 4038 | cavity_.clear(); | |
| 658 | |||
| 659 |
1/2✓ Branch 1 taken 4038 times.
✗ Branch 2 not taken.
|
4038 | bool ok = find_conflict_zone(v,t,t_bndry,f_bndry); |
| 660 | |||
| 661 | // When in multithreading mode, we cannot allocate memory | ||
| 662 | // dynamically and we use a fixed pool. If the fixed pool | ||
| 663 | // is full, then we exit the thread (and the missing points | ||
| 664 | // are inserted after, in sequential mode). | ||
| 665 | 4038 | if( | |
| 666 |
3/4✓ Branch 0 taken 2429 times.
✓ Branch 1 taken 1609 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4038 times.
|
6467 | nb_tets_to_create_ > nb_free_ && |
| 667 |
2/4✓ Branch 1 taken 2429 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 2429 times.
|
2429 | Process::is_running_threads() |
| 668 | ) { | ||
| 669 | ✗ | memory_overflow_ = true; | |
| 670 | ✗ | ok = false; | |
| 671 | } | ||
| 672 | |||
| 673 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4038 times.
|
4038 | if(!ok) { |
| 674 | // At this point, this thread did not successfully | ||
| 675 | // acquire all the tets in the conflict zone, so | ||
| 676 | // we need to rollback. | ||
| 677 | ✗ | release_tets(); | |
| 678 | ✗ | geo_debug_assert(nb_acquired_tets_ == 0); | |
| 679 | ✗ | return false; | |
| 680 | } | ||
| 681 | |||
| 682 | // The conflict list can be empty if | ||
| 683 | // the triangulation is weighted and v is not visible | ||
| 684 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 4038 times.
|
4038 | if(tets_to_delete_.size() == 0) { |
| 685 | ✗ | release_tets(); | |
| 686 | ✗ | geo_debug_assert(nb_acquired_tets_ == 0); | |
| 687 | ✗ | return true; | |
| 688 | } | ||
| 689 | |||
| 690 | |||
| 691 |
1/6✗ Branch 2 not taken.
✓ Branch 3 taken 4038 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
4038 | geo_debug_assert( |
| 692 | nb_acquired_tets_ == | ||
| 693 | tets_to_delete_.size() + tets_to_release_.size() | ||
| 694 | ); | ||
| 695 | |||
| 696 | #ifdef GEO_DEBUG | ||
| 697 | // Sanity check: make sure this threads owns all the tets | ||
| 698 | // in conflict and their neighbors. | ||
| 699 |
2/2✓ Branch 1 taken 77597 times.
✓ Branch 2 taken 4038 times.
|
81635 | for(index_t i=0; i<tets_to_delete_.size(); ++i) { |
| 700 |
1/2✓ Branch 1 taken 77597 times.
✗ Branch 2 not taken.
|
77597 | index_t tdel = tets_to_delete_[i]; |
| 701 |
2/8✓ Branch 1 taken 77597 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 77597 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
77597 | geo_debug_assert(owns_tet(tdel)); |
| 702 |
2/2✓ Branch 0 taken 310388 times.
✓ Branch 1 taken 77597 times.
|
387985 | for(index_t lf=0; lf<4; ++lf) { |
| 703 |
2/8✓ Branch 1 taken 310388 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 310388 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
310388 | geo_debug_assert(tet_adjacent(tdel,lf) != NO_INDEX); |
| 704 |
3/10✓ Branch 1 taken 310388 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 310388 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 310388 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
|
310388 | geo_debug_assert(owns_tet(tet_adjacent(tdel,lf))); |
| 705 | } | ||
| 706 | } | ||
| 707 | #endif | ||
| 708 |
2/8✓ Branch 1 taken 4038 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 4038 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
4038 | geo_debug_assert(owns_tet(t_bndry)); |
| 709 |
3/10✓ Branch 1 taken 4038 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4038 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 4038 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
|
4038 | geo_debug_assert(owns_tet(tet_adjacent(t_bndry,f_bndry))); |
| 710 |
3/10✓ Branch 1 taken 4038 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4038 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 4038 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
|
4038 | geo_debug_assert( |
| 711 | !tet_is_marked_as_conflict(tet_adjacent(t_bndry,f_bndry)) | ||
| 712 | ); | ||
| 713 | |||
| 714 | // At this point, this threads owns all the tets in conflict and | ||
| 715 | // their neighbors, therefore no other thread can interfere, and | ||
| 716 | // we can update the triangulation. | ||
| 717 | |||
| 718 | 4038 | index_t new_tet = NO_INDEX; | |
| 719 |
2/2✓ Branch 1 taken 4032 times.
✓ Branch 2 taken 6 times.
|
4038 | if(cavity_.OK()) { |
| 720 |
1/2✓ Branch 1 taken 4032 times.
✗ Branch 2 not taken.
|
4032 | new_tet = stellate_cavity(v); |
| 721 | } else { | ||
| 722 |
1/2✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
|
6 | new_tet = stellate_conflict_zone_iterative(v,t_bndry,f_bndry); |
| 723 | } | ||
| 724 | |||
| 725 | |||
| 726 | // Recycle the tetrahedra of the conflict zone. | ||
| 727 |
2/2✓ Branch 1 taken 73559 times.
✓ Branch 2 taken 4038 times.
|
77597 | for(index_t i=0; i<tets_to_delete_.size()-1; ++i) { |
| 728 |
3/6✓ Branch 1 taken 73559 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 73559 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 73559 times.
✗ Branch 8 not taken.
|
73559 | cell_next_[tets_to_delete_[i]] = tets_to_delete_[i+1]; |
| 729 | } | ||
| 730 |
2/4✓ Branch 1 taken 4038 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4038 times.
✗ Branch 5 not taken.
|
4038 | cell_next_[tets_to_delete_[tets_to_delete_.size()-1]] = |
| 731 | 4038 | first_free_; | |
| 732 |
1/2✓ Branch 1 taken 4038 times.
✗ Branch 2 not taken.
|
4038 | first_free_ = tets_to_delete_[0]; |
| 733 |
1/2✓ Branch 1 taken 4038 times.
✗ Branch 2 not taken.
|
4038 | nb_free_ += nb_tets_in_conflict(); |
| 734 | |||
| 735 | // For debugging purposes. | ||
| 736 | #ifdef GEO_DEBUG | ||
| 737 |
2/2✓ Branch 1 taken 77597 times.
✓ Branch 2 taken 4038 times.
|
81635 | for(index_t i=0; i<tets_to_delete_.size(); ++i) { |
| 738 |
1/2✓ Branch 1 taken 77597 times.
✗ Branch 2 not taken.
|
77597 | index_t tdel = tets_to_delete_[i]; |
| 739 |
1/2✓ Branch 1 taken 77597 times.
✗ Branch 2 not taken.
|
77597 | set_tet_vertex(tdel,0,VERTEX_OF_DELETED_TET); |
| 740 |
1/2✓ Branch 1 taken 77597 times.
✗ Branch 2 not taken.
|
77597 | set_tet_vertex(tdel,1,VERTEX_OF_DELETED_TET); |
| 741 |
1/2✓ Branch 1 taken 77597 times.
✗ Branch 2 not taken.
|
77597 | set_tet_vertex(tdel,2,VERTEX_OF_DELETED_TET); |
| 742 |
1/2✓ Branch 1 taken 77597 times.
✗ Branch 2 not taken.
|
77597 | set_tet_vertex(tdel,3,VERTEX_OF_DELETED_TET); |
| 743 | } | ||
| 744 | #endif | ||
| 745 | |||
| 746 | // Return one of the newly created tets | ||
| 747 | 4038 | hint=new_tet; | |
| 748 | |||
| 749 |
1/2✓ Branch 1 taken 4038 times.
✗ Branch 2 not taken.
|
4038 | release_tets(); |
| 750 | |||
| 751 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 4038 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
4038 | geo_debug_assert(nb_acquired_tets_ == 0); |
| 752 | 4038 | return true; | |
| 753 | } | ||
| 754 | |||
| 755 | /** | ||
| 756 | * \brief Determines the list of tetrahedra in conflict | ||
| 757 | * with a given point. | ||
| 758 | * \param[in] v the index of the point to be inserted | ||
| 759 | * \param[in] t the index of a tetrahedron that contains | ||
| 760 | * \p p, as returned by locate() | ||
| 761 | * \param[out] t_bndry a tetrahedron adjacent to the | ||
| 762 | * boundary of the conflict zone | ||
| 763 | * \param[out] f_bndry the facet along which t_bndry is | ||
| 764 | * adjacent to the boundary of the conflict zone | ||
| 765 | * The other tetrahedra are linked, and can be traversed | ||
| 766 | * from \p first by using tet_next() until \p last or END_OF_LIST | ||
| 767 | * is reached. | ||
| 768 | * The conflict zone can be empty under two circumstances: | ||
| 769 | * - the vertex \p v already exists in the triangulation | ||
| 770 | * - the triangulation is weighted and \p v is not visible | ||
| 771 | * in either cases, both \p first and \p last contain END_OF_LIST | ||
| 772 | * \retval true if all the tetrahedra of the conflict zone and their | ||
| 773 | * neighbors could be acquired by this thread | ||
| 774 | * \retval false otherwise | ||
| 775 | */ | ||
| 776 | 4038 | bool find_conflict_zone( | |
| 777 | index_t v, index_t t, | ||
| 778 | index_t& t_bndry, index_t& f_bndry | ||
| 779 | ) { | ||
| 780 | 4038 | nb_tets_to_create_ = 0; | |
| 781 | |||
| 782 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 4038 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
4038 | geo_debug_assert(t != NO_TETRAHEDRON); |
| 783 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 4038 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
4038 | geo_debug_assert(owns_tet(t)); |
| 784 | |||
| 785 | // Pointer to the coordinates of the point to be inserted | ||
| 786 | 4038 | const double* p = vertex_ptr(v); | |
| 787 | |||
| 788 | // Weighted triangulations can have dangling | ||
| 789 | // vertices. Such vertices p are characterized by | ||
| 790 | // the fact that p is not in conflict with the | ||
| 791 | // tetrahedron returned by locate(). | ||
| 792 |
2/6✗ Branch 0 not taken.
✓ Branch 1 taken 4038 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 4038 times.
|
4038 | if(weighted_ && !tet_is_in_conflict(t,p)) { |
| 793 | ✗ | release_tet(t); | |
| 794 | ✗ | return true; | |
| 795 | } | ||
| 796 | |||
| 797 | 4038 | mark_tet_as_conflict(t); | |
| 798 | |||
| 799 | // Sanity check: the vertex to be inserted should | ||
| 800 | // not correspond to one of the vertices of t. | ||
| 801 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 4038 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
4038 | geo_debug_assert(v != tet_vertex(t,0)); |
| 802 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 4038 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
4038 | geo_debug_assert(v != tet_vertex(t,1)); |
| 803 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 4038 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
4038 | geo_debug_assert(v != tet_vertex(t,2)); |
| 804 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 4038 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
4038 | geo_debug_assert(v != tet_vertex(t,3)); |
| 805 | |||
| 806 | // Note: points on edges and on facets are | ||
| 807 | // handled by the way tet_is_in_conflict() | ||
| 808 | // is implemented, that naturally inserts | ||
| 809 | // the correct tetrahedra in the conflict list. | ||
| 810 | |||
| 811 | // Determine the conflict list by greedy propagation from t. | ||
| 812 | 4038 | bool result = find_conflict_zone_iterative(p,t); | |
| 813 | 4038 | t_bndry = t_boundary_; | |
| 814 | 4038 | f_bndry = f_boundary_; | |
| 815 | 4038 | return result; | |
| 816 | } | ||
| 817 | |||
| 818 | |||
| 819 | /** | ||
| 820 | * \brief This function is used to implement find_conflict_zone. | ||
| 821 | * \details This function detects the neighbors of \p t that are | ||
| 822 | * in the conflict zone and calls itself recursively on them. | ||
| 823 | * \param[in] p the point to be inserted | ||
| 824 | * \param[in] t_in index of a tetrahedron in the conflict zone | ||
| 825 | * \pre The tetrahedron \p t was alredy marked as | ||
| 826 | * conflict (tet_is_in_list(t)) | ||
| 827 | */ | ||
| 828 | 4038 | bool find_conflict_zone_iterative( | |
| 829 | const double* p, index_t t_in | ||
| 830 | ) { | ||
| 831 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 4038 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
4038 | geo_debug_assert(owns_tet(t_in)); |
| 832 | 4038 | S_.push_back(t_in); | |
| 833 | |||
| 834 |
2/2✓ Branch 1 taken 77597 times.
✓ Branch 2 taken 4038 times.
|
81635 | while(S_.size() != 0) { |
| 835 | 77597 | index_t t = *(S_.rbegin()); | |
| 836 | 77597 | S_.pop_back(); | |
| 837 | |||
| 838 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 77597 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
77597 | geo_debug_assert(owns_tet(t)); |
| 839 | |||
| 840 |
2/2✓ Branch 0 taken 310388 times.
✓ Branch 1 taken 77597 times.
|
387985 | for(index_t lf = 0; lf < 4; ++lf) { |
| 841 |
1/2✓ Branch 1 taken 310388 times.
✗ Branch 2 not taken.
|
310388 | index_t t2 = tet_adjacent(t, lf); |
| 842 | |||
| 843 | // If t2 is already owned by current thread, then | ||
| 844 | // its status was previously determined. | ||
| 845 |
3/4✓ Branch 1 taken 310388 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 148233 times.
✓ Branch 4 taken 162155 times.
|
310388 | if(owns_tet(t2)) { |
| 846 |
3/10✓ Branch 1 taken 148233 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 148233 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 148233 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
|
148233 | geo_debug_assert( |
| 847 | tet_is_marked_as_conflict(t2) == | ||
| 848 | tet_is_in_conflict(t2,p) | ||
| 849 | ); | ||
| 850 | |||
| 851 | // If t2 is not in conflict list, then t has a facet | ||
| 852 | // on the border of the conflict zone, and there is | ||
| 853 | // a tet to create. | ||
| 854 |
3/4✓ Branch 1 taken 148233 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 14270 times.
✓ Branch 4 taken 133963 times.
|
148233 | if(!tet_is_marked_as_conflict(t2)) { |
| 855 | 14270 | ++nb_tets_to_create_; | |
| 856 |
7/14✓ Branch 1 taken 14270 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 14270 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 14270 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 14270 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 14270 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 14270 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 14270 times.
✗ Branch 20 not taken.
|
14270 | cavity_.new_facet( |
| 857 | t, lf, | ||
| 858 | tet_vertex(t, tet_facet_vertex(lf,0)), | ||
| 859 | tet_vertex(t, tet_facet_vertex(lf,1)), | ||
| 860 | tet_vertex(t, tet_facet_vertex(lf,2)) | ||
| 861 | ); | ||
| 862 | } | ||
| 863 | 221792 | continue; | |
| 864 | } | ||
| 865 | |||
| 866 |
2/4✓ Branch 1 taken 162155 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 162155 times.
|
162155 | if(!acquire_tet(t2)) { |
| 867 | ✗ | S_.resize(0); | |
| 868 | ✗ | return false; | |
| 869 | } | ||
| 870 | |||
| 871 |
2/8✓ Branch 1 taken 162155 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 162155 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
162155 | geo_debug_assert(owns_tet(t2)); |
| 872 | |||
| 873 |
3/4✓ Branch 1 taken 162155 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 88596 times.
✓ Branch 4 taken 73559 times.
|
162155 | if(!tet_is_in_conflict(t2,p)) { |
| 874 |
1/2✓ Branch 1 taken 88596 times.
✗ Branch 2 not taken.
|
88596 | mark_tet_as_neighbor(t2); |
| 875 | // If t2 is not in conflict list, then t has a facet | ||
| 876 | // on the border of the conflict zone, and there is | ||
| 877 | // a tet to create. | ||
| 878 | 88596 | ++nb_tets_to_create_; | |
| 879 | } else { | ||
| 880 |
1/2✓ Branch 1 taken 73559 times.
✗ Branch 2 not taken.
|
73559 | mark_tet_as_conflict(t2); |
| 881 |
2/8✓ Branch 1 taken 73559 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 73559 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
73559 | geo_debug_assert(owns_tet(t2)); |
| 882 |
1/2✓ Branch 1 taken 73559 times.
✗ Branch 2 not taken.
|
73559 | S_.push_back(t2); |
| 883 | 73559 | continue; | |
| 884 | } | ||
| 885 | |||
| 886 | // At this point, t is in conflict | ||
| 887 | // and t2 is not in conflict. | ||
| 888 | // We keep a reference to a tet on the boundary | ||
| 889 | 88596 | t_boundary_ = t; | |
| 890 | 88596 | f_boundary_ = lf; | |
| 891 | 88596 | ++nb_tets_to_create_; | |
| 892 |
7/14✓ Branch 1 taken 88596 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 88596 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 88596 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 88596 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 88596 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 88596 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 88596 times.
✗ Branch 20 not taken.
|
88596 | cavity_.new_facet( |
| 893 | t, lf, | ||
| 894 | tet_vertex(t, tet_facet_vertex(lf,0)), | ||
| 895 | tet_vertex(t, tet_facet_vertex(lf,1)), | ||
| 896 | tet_vertex(t, tet_facet_vertex(lf,2)) | ||
| 897 | ); | ||
| 898 |
2/8✓ Branch 1 taken 88596 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 88596 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
88596 | geo_debug_assert(tet_adjacent(t,lf) == t2); |
| 899 |
2/8✓ Branch 1 taken 88596 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 88596 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
88596 | geo_debug_assert(owns_tet(t)); |
| 900 |
2/8✓ Branch 1 taken 88596 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 88596 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
88596 | geo_debug_assert(owns_tet(t2)); |
| 901 | } | ||
| 902 | } | ||
| 903 | 4038 | return true; | |
| 904 | } | ||
| 905 | |||
| 906 | /** | ||
| 907 | * \brief Gets the lifted coordinate of a point by its 3d coordinates. | ||
| 908 | * \param[in] p a pointer to the coordinates of one of the vertices | ||
| 909 | * of the triangulation. | ||
| 910 | * \return the lifted coordinate of \p p | ||
| 911 | */ | ||
| 912 | ✗ | double lifted_coordinate(const double* p) const { | |
| 913 | // Compute the index of the point from its address | ||
| 914 | index_t pindex = index_t( | ||
| 915 | ✗ | (p - vertex_ptr(0)) / int(vertex_stride_) | |
| 916 | ✗ | ); | |
| 917 | ✗ | return heights_[pindex]; | |
| 918 | } | ||
| 919 | |||
| 920 | |||
| 921 | /** | ||
| 922 | * \brief Tests whether a given tetrahedron is in conflict with | ||
| 923 | * a given 3d point. | ||
| 924 | * \details A real tetrahedron is in conflict with a point whenever | ||
| 925 | * the point is contained by its circumscribed sphere, and a | ||
| 926 | * virtual tetrahedron is in conflict with a point whenever the | ||
| 927 | * tetrahedron formed by its real face and with the point has | ||
| 928 | * positive orientation. | ||
| 929 | * \param[in] t the index of the tetrahedron | ||
| 930 | * \param[in] p a pointer to the coordinates of the point | ||
| 931 | * \retval true if point \p p is in conflict with tetrahedron \p t | ||
| 932 | * \retval false otherwise | ||
| 933 | */ | ||
| 934 | 804420 | bool tet_is_in_conflict(index_t t, const double* p) const { | |
| 935 | |||
| 936 | // Lookup tetrahedron vertices | ||
| 937 | const double* pv[4]; | ||
| 938 |
2/2✓ Branch 0 taken 3217680 times.
✓ Branch 1 taken 804420 times.
|
4022100 | for(index_t i=0; i<4; ++i) { |
| 939 |
1/2✓ Branch 1 taken 3217680 times.
✗ Branch 2 not taken.
|
3217680 | index_t v = tet_vertex(t,i); |
| 940 |
3/4✓ Branch 0 taken 3038879 times.
✓ Branch 1 taken 178801 times.
✓ Branch 3 taken 3038879 times.
✗ Branch 4 not taken.
|
3217680 | pv[i] = (v == NO_INDEX) ? nullptr : vertex_ptr(v); |
| 941 | } | ||
| 942 | |||
| 943 | // Check for virtual tetrahedra (then in_sphere() | ||
| 944 | // is replaced with orient3d()) | ||
| 945 |
2/2✓ Branch 0 taken 3030651 times.
✓ Branch 1 taken 625619 times.
|
3656270 | for(index_t lf = 0; lf < 4; ++lf) { |
| 946 | |||
| 947 |
2/2✓ Branch 0 taken 178801 times.
✓ Branch 1 taken 2851850 times.
|
3030651 | if(pv[lf] == nullptr) { |
| 948 | |||
| 949 | // Facet of a virtual tetrahedron opposite to | ||
| 950 | // infinite vertex corresponds to | ||
| 951 | // the triangle on the convex hull of the points. | ||
| 952 | // Orientation is obtained by replacing vertex lf | ||
| 953 | // with p. | ||
| 954 | 178801 | pv[lf] = p; | |
| 955 |
1/2✓ Branch 1 taken 178801 times.
✗ Branch 2 not taken.
|
178801 | Sign sign = PCK::orient_3d(pv[0],pv[1],pv[2],pv[3]); |
| 956 | |||
| 957 |
2/2✓ Branch 0 taken 9667 times.
✓ Branch 1 taken 169134 times.
|
178801 | if(sign > 0) { |
| 958 | 9667 | return true; | |
| 959 | } | ||
| 960 | |||
| 961 |
2/2✓ Branch 0 taken 124522 times.
✓ Branch 1 taken 44612 times.
|
169134 | if(sign < 0) { |
| 962 | 124522 | return false; | |
| 963 | } | ||
| 964 | |||
| 965 | // If sign is zero, we check the real tetrahedron | ||
| 966 | // adjacent to the facet on the convex hull. | ||
| 967 |
2/8✓ Branch 1 taken 44612 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 44612 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
44612 | geo_debug_assert(tet_adjacent(t, lf) != NO_INDEX); |
| 968 |
1/2✓ Branch 1 taken 44612 times.
✗ Branch 2 not taken.
|
44612 | index_t t2 = tet_adjacent(t, lf); |
| 969 |
2/8✓ Branch 1 taken 44612 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 44612 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
44612 | geo_debug_assert(!tet_is_virtual(t2)); |
| 970 | |||
| 971 | // If t2 was already visited by this thread, then | ||
| 972 | // it is in conflict if it is already marked. | ||
| 973 |
3/4✓ Branch 1 taken 44612 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4374 times.
✓ Branch 4 taken 40238 times.
|
44612 | if(owns_tet(t2)) { |
| 974 |
1/2✓ Branch 1 taken 4374 times.
✗ Branch 2 not taken.
|
4374 | return tet_is_marked_as_conflict(t2); |
| 975 | } | ||
| 976 | |||
| 977 | // If t2 was not already visited, then we need to | ||
| 978 | // switch to the in_circum_circle_3d() predicate. | ||
| 979 | |||
| 980 | 40238 | const double* q0 = pv[(lf+1)%4]; | |
| 981 | 40238 | const double* q1 = pv[(lf+2)%4]; | |
| 982 | 40238 | const double* q2 = pv[(lf+3)%4]; | |
| 983 | |||
| 984 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 40238 times.
|
40238 | if(weighted_) { |
| 985 | return ( | ||
| 986 | ✗ | PCK::in_circle_3dlifted_SOS( | |
| 987 | q0, q1, q2, p, | ||
| 988 | lifted_coordinate(q0), | ||
| 989 | lifted_coordinate(q1), | ||
| 990 | lifted_coordinate(q2), | ||
| 991 | lifted_coordinate(p) | ||
| 992 | ) > 0 | ||
| 993 | ✗ | ); | |
| 994 | } else { | ||
| 995 | return ( | ||
| 996 |
1/2✓ Branch 1 taken 40238 times.
✗ Branch 2 not taken.
|
40238 | PCK::in_circle_3d_SOS( |
| 997 | q0,q1,q2,p | ||
| 998 | ) > 0 | ||
| 999 | 40238 | ); | |
| 1000 | } | ||
| 1001 | } | ||
| 1002 | } | ||
| 1003 | |||
| 1004 | // If the tetrahedron is a finite one, it is in conflict | ||
| 1005 | // if its circumscribed sphere contains the point (this is | ||
| 1006 | // the standard case). | ||
| 1007 | |||
| 1008 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 625619 times.
|
625619 | if(weighted_) { |
| 1009 | ✗ | double h0 = heights_[finite_tet_vertex(t, 0)]; | |
| 1010 | ✗ | double h1 = heights_[finite_tet_vertex(t, 1)]; | |
| 1011 | ✗ | double h2 = heights_[finite_tet_vertex(t, 2)]; | |
| 1012 | ✗ | double h3 = heights_[finite_tet_vertex(t, 3)]; | |
| 1013 | ✗ | double h = lifted_coordinate(p); | |
| 1014 | ✗ | return (PCK::orient_3dlifted_SOS( | |
| 1015 | pv[0],pv[1],pv[2],pv[3],p,h0,h1,h2,h3,h | ||
| 1016 | ✗ | ) > 0) ; | |
| 1017 | } | ||
| 1018 | |||
| 1019 |
1/2✓ Branch 1 taken 625619 times.
✗ Branch 2 not taken.
|
625619 | return (PCK::in_sphere_3d_SOS(pv[0], pv[1], pv[2], pv[3], p) > 0); |
| 1020 | } | ||
| 1021 | |||
| 1022 | |||
| 1023 | /** | ||
| 1024 | * \brief Finds the tetrahedron that contains a point. | ||
| 1025 | * \details The tetrahedron is acquired by this thread. If the | ||
| 1026 | * tetrahedron could not be acquired, then NO_TETRAHEDRON is returned. | ||
| 1027 | * If the point is on a face, edge or vertex, | ||
| 1028 | * the function returns one of the tetrahedra incident | ||
| 1029 | * to that face, edge or vertex. | ||
| 1030 | * \param[in] p a pointer to the coordinates of the point | ||
| 1031 | * \param[out] orient a pointer to an array of four Sign%s | ||
| 1032 | * or nullptr. If non-nullptr, returns the orientation with respect | ||
| 1033 | * to the four facets of the tetrahedron that contains \p p. | ||
| 1034 | * \retval the index of a tetrahedron that contains \p p. | ||
| 1035 | * If the point is outside the convex hull of | ||
| 1036 | * the inserted so-far points, then the returned tetrahedron | ||
| 1037 | * is a virtual one (first vertex is the "vertex at infinity" | ||
| 1038 | * of index -1) | ||
| 1039 | * \retval NO_TETRAHEDRON if the tetrahedron could not be | ||
| 1040 | * acquired by this thread, or if the virtual tetrahedra | ||
| 1041 | * were previously removed | ||
| 1042 | */ | ||
| 1043 | 4038 | index_t locate( | |
| 1044 | const double* p, index_t hint = NO_TETRAHEDRON, | ||
| 1045 | Sign* orient = nullptr | ||
| 1046 | ) { | ||
| 1047 | // Try improving the hint by using the | ||
| 1048 | // inexact locate function. This gains | ||
| 1049 | // (a little bit) performance (a few | ||
| 1050 | // percent in total Delaunay computation | ||
| 1051 | // time), but it is better than nothing... | ||
| 1052 | // Note: there is a maximum number of tets | ||
| 1053 | // traversed by locate_inexact() (2500) | ||
| 1054 | // since there exists configurations in which | ||
| 1055 | // locate_inexact() loops forever ! | ||
| 1056 | |||
| 1057 | { | ||
| 1058 |
1/2✓ Branch 1 taken 4038 times.
✗ Branch 2 not taken.
|
4038 | index_t new_hint = locate_inexact(p, hint, 2500); |
| 1059 | |||
| 1060 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4038 times.
|
4038 | if(new_hint == NO_TETRAHEDRON) { |
| 1061 | ✗ | return NO_TETRAHEDRON; | |
| 1062 | } | ||
| 1063 | |||
| 1064 | 4038 | hint = new_hint; | |
| 1065 | } | ||
| 1066 | |||
| 1067 | // If no hint specified, find a tetrahedron randomly | ||
| 1068 | |||
| 1069 |
1/2✓ Branch 0 taken 4038 times.
✗ Branch 1 not taken.
|
4038 | if(hint != NO_TETRAHEDRON) { |
| 1070 |
2/4✓ Branch 1 taken 4038 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 4038 times.
|
4038 | if(tet_is_free(hint)) { |
| 1071 | ✗ | hint = NO_TETRAHEDRON; | |
| 1072 | } else { | ||
| 1073 |
5/10✓ Branch 1 taken 4038 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4038 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 4038 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 4038 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 4038 times.
|
4038 | if( !owns_tet(hint) && !acquire_tet(hint) ) { |
| 1074 | ✗ | hint = NO_TETRAHEDRON; | |
| 1075 | } | ||
| 1076 |
4/8✓ Branch 0 taken 4038 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 4038 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 4038 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 4038 times.
|
4038 | if((hint != NO_TETRAHEDRON) && tet_is_free(hint)) { |
| 1077 | ✗ | release_tet(hint); | |
| 1078 | ✗ | hint = NO_TETRAHEDRON; | |
| 1079 | } | ||
| 1080 | } | ||
| 1081 | } | ||
| 1082 | |||
| 1083 | do { | ||
| 1084 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4038 times.
|
4038 | if(hint == NO_TETRAHEDRON) { |
| 1085 | ✗ | hint = thread_safe_random(max_used_t_); | |
| 1086 | } | ||
| 1087 | 4038 | if( | |
| 1088 |
3/6✓ Branch 1 taken 4038 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4038 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 4038 times.
|
8076 | tet_is_free(hint) || |
| 1089 |
2/8✓ Branch 1 taken 4038 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 4038 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
4038 | (!owns_tet(hint) && !acquire_tet(hint)) |
| 1090 | ) { | ||
| 1091 | ✗ | if(owns_tet(hint)) { | |
| 1092 | ✗ | release_tet(hint); | |
| 1093 | } | ||
| 1094 | ✗ | hint = NO_TETRAHEDRON; | |
| 1095 | } else { | ||
| 1096 |
2/2✓ Branch 0 taken 15114 times.
✓ Branch 1 taken 3107 times.
|
18221 | for(index_t f=0; f<4; ++f) { |
| 1097 |
3/4✓ Branch 1 taken 15114 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 931 times.
✓ Branch 4 taken 14183 times.
|
15114 | if(tet_vertex(hint,f) == VERTEX_AT_INFINITY) { |
| 1098 |
1/2✓ Branch 1 taken 931 times.
✗ Branch 2 not taken.
|
931 | index_t new_hint = tet_adjacent(hint,f); |
| 1099 | 931 | if( | |
| 1100 |
3/6✓ Branch 1 taken 931 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 931 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 931 times.
|
1862 | tet_is_free(new_hint) || |
| 1101 |
2/4✓ Branch 1 taken 931 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 931 times.
|
931 | !acquire_tet(new_hint) |
| 1102 | ) { | ||
| 1103 | ✗ | new_hint = NO_TETRAHEDRON; | |
| 1104 | } | ||
| 1105 |
1/2✓ Branch 1 taken 931 times.
✗ Branch 2 not taken.
|
931 | release_tet(hint); |
| 1106 | 931 | hint = new_hint; | |
| 1107 | 931 | break; | |
| 1108 | } | ||
| 1109 | } | ||
| 1110 | } | ||
| 1111 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4038 times.
|
4038 | } while(hint == NO_TETRAHEDRON) ; |
| 1112 | |||
| 1113 | 4038 | index_t t = hint; | |
| 1114 | 4038 | index_t t_pred = NO_TETRAHEDRON; | |
| 1115 | Sign orient_local[4]; | ||
| 1116 |
1/2✓ Branch 0 taken 4038 times.
✗ Branch 1 not taken.
|
4038 | if(orient == nullptr) { |
| 1117 | ✗ | orient = orient_local; | |
| 1118 | } | ||
| 1119 | |||
| 1120 | |||
| 1121 | 4038 | still_walking: | |
| 1122 | { | ||
| 1123 |
2/2✓ Branch 0 taken 81 times.
✓ Branch 1 taken 4038 times.
|
4119 | if(t_pred != NO_TETRAHEDRON) { |
| 1124 |
1/2✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
|
81 | release_tet(t_pred); |
| 1125 | } | ||
| 1126 | |||
| 1127 |
2/4✓ Branch 1 taken 4119 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 4119 times.
|
4119 | if(tet_is_free(t)) { |
| 1128 | 931 | return NO_TETRAHEDRON; | |
| 1129 | } | ||
| 1130 | |||
| 1131 |
6/10✓ Branch 1 taken 4119 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 81 times.
✓ Branch 4 taken 4038 times.
✓ Branch 6 taken 81 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 81 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 4119 times.
|
4119 | if(!owns_tet(t) && !acquire_tet(t)) { |
| 1132 | ✗ | return NO_TETRAHEDRON; | |
| 1133 | } | ||
| 1134 | |||
| 1135 | |||
| 1136 |
2/4✓ Branch 1 taken 4119 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 4119 times.
|
4119 | if(!tet_is_real(t)) { |
| 1137 | ✗ | release_tet(t); | |
| 1138 | ✗ | return NO_TETRAHEDRON; | |
| 1139 | } | ||
| 1140 | |||
| 1141 | const double* pv[4]; | ||
| 1142 |
2/4✓ Branch 1 taken 4119 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4119 times.
✗ Branch 5 not taken.
|
4119 | pv[0] = vertex_ptr(finite_tet_vertex(t,0)); |
| 1143 |
2/4✓ Branch 1 taken 4119 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4119 times.
✗ Branch 5 not taken.
|
4119 | pv[1] = vertex_ptr(finite_tet_vertex(t,1)); |
| 1144 |
2/4✓ Branch 1 taken 4119 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4119 times.
✗ Branch 5 not taken.
|
4119 | pv[2] = vertex_ptr(finite_tet_vertex(t,2)); |
| 1145 |
2/4✓ Branch 1 taken 4119 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4119 times.
✗ Branch 5 not taken.
|
4119 | pv[3] = vertex_ptr(finite_tet_vertex(t,3)); |
| 1146 | |||
| 1147 | // Start from a random facet | ||
| 1148 | 4119 | index_t f0 = thread_safe_random_4(); | |
| 1149 |
2/2✓ Branch 0 taken 14601 times.
✓ Branch 1 taken 3107 times.
|
17708 | for(index_t df = 0; df < 4; ++df) { |
| 1150 | 14601 | index_t f = (f0 + df) % 4; | |
| 1151 | |||
| 1152 |
1/2✓ Branch 1 taken 14601 times.
✗ Branch 2 not taken.
|
14601 | index_t t_next = tet_adjacent(t,f); |
| 1153 | |||
| 1154 | // If the opposite tet is -1, then it means that | ||
| 1155 | // we are trying to locate() (e.g. called from | ||
| 1156 | // nearest_vertex) within a tetrahedralization | ||
| 1157 | // from which the infinite tets were removed. | ||
| 1158 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 14601 times.
|
14601 | if(t_next == NO_INDEX) { |
| 1159 | ✗ | release_tet(t); | |
| 1160 | ✗ | return NO_TETRAHEDRON; | |
| 1161 | } | ||
| 1162 | |||
| 1163 | // If the candidate next tetrahedron is the | ||
| 1164 | // one we came from, then we know already that | ||
| 1165 | // the orientation is positive, thus we examine | ||
| 1166 | // the next candidate (or exit the loop if they | ||
| 1167 | // are exhausted). | ||
| 1168 |
2/2✓ Branch 0 taken 37 times.
✓ Branch 1 taken 14564 times.
|
14601 | if(t_next == t_pred) { |
| 1169 | 37 | orient[f] = POSITIVE ; | |
| 1170 | 37 | continue ; | |
| 1171 | } | ||
| 1172 | |||
| 1173 | // To test the orientation of p w.r.t. the facet f of | ||
| 1174 | // t, we replace vertex number f with p in t (same | ||
| 1175 | // convention as in CGAL). | ||
| 1176 | // This is equivalent to tet_facet_point_orient3d(t,f,p) | ||
| 1177 | // (but less costly, saves a couple of lookups) | ||
| 1178 | 14564 | const double* pv_bkp = pv[f]; | |
| 1179 | 14564 | pv[f] = p; | |
| 1180 |
1/2✓ Branch 1 taken 14564 times.
✗ Branch 2 not taken.
|
14564 | orient[f] = PCK::orient_3d(pv[0], pv[1], pv[2], pv[3]); |
| 1181 | |||
| 1182 | // If the orientation is not negative, then we cannot | ||
| 1183 | // walk towards t_next, and examine the next candidate | ||
| 1184 | // (or exit the loop if they are exhausted). | ||
| 1185 |
2/2✓ Branch 0 taken 13552 times.
✓ Branch 1 taken 1012 times.
|
14564 | if(orient[f] != NEGATIVE) { |
| 1186 | 13552 | pv[f] = pv_bkp; | |
| 1187 | 13552 | continue; | |
| 1188 | } | ||
| 1189 | |||
| 1190 | // If the opposite tet is a virtual tet, then | ||
| 1191 | // the point has a positive orientation relative | ||
| 1192 | // to the facet on the border of the convex hull, | ||
| 1193 | // thus t_next is a tet in conflict and we are | ||
| 1194 | // done. | ||
| 1195 |
3/4✓ Branch 1 taken 1012 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 931 times.
✓ Branch 4 taken 81 times.
|
1012 | if(tet_is_virtual(t_next)) { |
| 1196 |
1/2✓ Branch 1 taken 931 times.
✗ Branch 2 not taken.
|
931 | release_tet(t); |
| 1197 |
2/4✓ Branch 1 taken 931 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 931 times.
|
931 | if(!acquire_tet(t_next)) { |
| 1198 | ✗ | return NO_TETRAHEDRON; | |
| 1199 | } | ||
| 1200 |
2/2✓ Branch 0 taken 3724 times.
✓ Branch 1 taken 931 times.
|
4655 | for(index_t lf = 0; lf < 4; ++lf) { |
| 1201 | 3724 | orient[lf] = POSITIVE; | |
| 1202 | } | ||
| 1203 | 931 | return t_next; | |
| 1204 | } | ||
| 1205 | |||
| 1206 | // If we reach this point, then t_next is a valid | ||
| 1207 | // successor, thus we are still walking. | ||
| 1208 | 81 | t_pred = t; | |
| 1209 | 81 | t = t_next; | |
| 1210 | 81 | goto still_walking; | |
| 1211 | } | ||
| 1212 | } | ||
| 1213 | |||
| 1214 | // If we reach this point, we did not find a valid successor | ||
| 1215 | // for walking (a face for which p has negative orientation), | ||
| 1216 | // thus we reached the tet for which p has all positive | ||
| 1217 | // face orientations (i.e. the tet that contains p). | ||
| 1218 | |||
| 1219 | #ifdef GEO_DEBUG | ||
| 1220 |
2/8✓ Branch 1 taken 3107 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 3107 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
3107 | geo_debug_assert(tet_is_real(t)); |
| 1221 | |||
| 1222 | const double* pv[4]; | ||
| 1223 | Sign signs[4]; | ||
| 1224 |
2/4✓ Branch 1 taken 3107 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3107 times.
✗ Branch 5 not taken.
|
3107 | pv[0] = vertex_ptr(finite_tet_vertex(t,0)); |
| 1225 |
2/4✓ Branch 1 taken 3107 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3107 times.
✗ Branch 5 not taken.
|
3107 | pv[1] = vertex_ptr(finite_tet_vertex(t,1)); |
| 1226 |
2/4✓ Branch 1 taken 3107 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3107 times.
✗ Branch 5 not taken.
|
3107 | pv[2] = vertex_ptr(finite_tet_vertex(t,2)); |
| 1227 |
2/4✓ Branch 1 taken 3107 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3107 times.
✗ Branch 5 not taken.
|
3107 | pv[3] = vertex_ptr(finite_tet_vertex(t,3)); |
| 1228 |
2/2✓ Branch 0 taken 12428 times.
✓ Branch 1 taken 3107 times.
|
15535 | for(index_t f=0; f<4; ++f) { |
| 1229 | 12428 | const double* pv_bkp = pv[f]; | |
| 1230 | 12428 | pv[f] = p; | |
| 1231 |
1/2✓ Branch 1 taken 12428 times.
✗ Branch 2 not taken.
|
12428 | signs[f] = PCK::orient_3d(pv[0], pv[1], pv[2], pv[3]); |
| 1232 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 12428 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
12428 | geo_debug_assert(signs[f] >= 0); |
| 1233 | 12428 | pv[f] = pv_bkp; | |
| 1234 | } | ||
| 1235 | #endif | ||
| 1236 | |||
| 1237 | 3107 | return t; | |
| 1238 | } | ||
| 1239 | |||
| 1240 | |||
| 1241 | protected: | ||
| 1242 | |||
| 1243 | /** | ||
| 1244 | * \brief Tests whether a tetrahedron was marked as conflict. | ||
| 1245 | * \pre owns_tet(t) | ||
| 1246 | * \param[in] t the index of the tetrahedron to be tested | ||
| 1247 | * \retval true if \p t was marked as conflict | ||
| 1248 | * \retval false otherwise | ||
| 1249 | */ | ||
| 1250 | 388261 | bool tet_is_marked_as_conflict(index_t t) const { | |
| 1251 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 388261 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
388261 | geo_debug_assert(owns_tet(t)); |
| 1252 | 388261 | return cell_status_.cell_is_marked_as_conflict(t); | |
| 1253 | } | ||
| 1254 | |||
| 1255 | |||
| 1256 | /** | ||
| 1257 | * \brief Gets the number of tetrahedra in conflict. | ||
| 1258 | * \return the number of tetrahedra in conflict, | ||
| 1259 | * specified by mark_tet_as_conflict() | ||
| 1260 | */ | ||
| 1261 | 4038 | index_t nb_tets_in_conflict() const { | |
| 1262 | 4038 | return tets_to_delete_.size(); | |
| 1263 | } | ||
| 1264 | |||
| 1265 | /** | ||
| 1266 | * \brief Marks a tetrahedron as conflict. | ||
| 1267 | * \details The index of the tetrahedron is also | ||
| 1268 | * stored it in the list of conflict tetrahedra. | ||
| 1269 | * \param[in] t index of the tetrahedron to mark | ||
| 1270 | * \pre owns_tet(t) | ||
| 1271 | */ | ||
| 1272 | 77597 | void mark_tet_as_conflict(index_t t) { | |
| 1273 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 77597 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
77597 | geo_debug_assert(owns_tet(t)); |
| 1274 | 77597 | tets_to_delete_.push_back(t); | |
| 1275 | 77597 | cell_status_.mark_cell_as_conflict(t); | |
| 1276 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 77597 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
77597 | geo_debug_assert(owns_tet(t)); |
| 1277 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 77597 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
77597 | geo_debug_assert(tet_is_marked_as_conflict(t)); |
| 1278 | 77597 | } | |
| 1279 | |||
| 1280 | /** | ||
| 1281 | * \brief Marks a tetrahedron as neighbor of the conflict zone. | ||
| 1282 | * \details The index of the tetrahedron is also | ||
| 1283 | * stored it in the list of tetrahedra to release. | ||
| 1284 | * \param[in] t index of the tetrahedron to mark | ||
| 1285 | * \pre owns_tet(t) | ||
| 1286 | */ | ||
| 1287 | 88596 | void mark_tet_as_neighbor(index_t t) { | |
| 1288 | // Note: nothing to change in cell_status_[t] | ||
| 1289 | // since LSB=0 means neigbhor tet. | ||
| 1290 | 88596 | tets_to_release_.push_back(t); | |
| 1291 | 88596 | } | |
| 1292 | |||
| 1293 | /** | ||
| 1294 | * \brief Acquires a lock on a tetrahedron and keep | ||
| 1295 | * it in the list of acquired tetrahedra. | ||
| 1296 | * \param[in] t index of the tetrahedron to acquire | ||
| 1297 | */ | ||
| 1298 | 102891 | void acquire_and_mark_tet_as_created(index_t t) { | |
| 1299 | // The tet was created in this thread's tet pool, | ||
| 1300 | // therefore there is no need to use sync | ||
| 1301 | // primitives to acquire a lock on it. | ||
| 1302 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 102891 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
102891 | geo_debug_assert(cell_status_.cell_thread(t) == NO_THREAD); |
| 1303 | 102891 | cell_status_.set_cell_status( | |
| 1304 | 102891 | t, CellStatusArray::thread_index_t(id()) | |
| 1305 | ); | ||
| 1306 | #ifdef GEO_DEBUG | ||
| 1307 | 102891 | ++nb_acquired_tets_; | |
| 1308 | #endif | ||
| 1309 | 102891 | tets_to_release_.push_back(t); | |
| 1310 | 102891 | } | |
| 1311 | |||
| 1312 | |||
| 1313 | /** | ||
| 1314 | * \brief Releases all the tetrahedron locks that were | ||
| 1315 | * acquired using mark_tet_as_neighbor(), | ||
| 1316 | * acquire_and_mark_tet_as_created() and mark_as_conflict(). | ||
| 1317 | */ | ||
| 1318 | 4043 | void release_tets() { | |
| 1319 |
2/2✓ Branch 1 taken 191487 times.
✓ Branch 2 taken 4043 times.
|
195530 | for(index_t i=0; i<tets_to_release_.size(); ++i) { |
| 1320 | 191487 | release_tet(tets_to_release_[i]); | |
| 1321 | } | ||
| 1322 | 4043 | tets_to_release_.resize(0); | |
| 1323 |
2/2✓ Branch 1 taken 77597 times.
✓ Branch 2 taken 4043 times.
|
81640 | for(index_t i=0; i<tets_to_delete_.size(); ++i) { |
| 1324 | 77597 | release_tet(tets_to_delete_[i]); | |
| 1325 | } | ||
| 1326 | 4043 | tets_to_delete_.resize(0); | |
| 1327 | 4043 | } | |
| 1328 | |||
| 1329 | /** | ||
| 1330 | * \brief Atomically acquires a lock on a tetrahedron. | ||
| 1331 | * \details When the lock could not be acquired, interfering_thread_ | ||
| 1332 | * contains the id of the thread that owns the lock. | ||
| 1333 | * \param[in] t the index of the tetrahedron to acquire | ||
| 1334 | * \retval true if the lock was successfully acquired | ||
| 1335 | * \retval false otherwise | ||
| 1336 | */ | ||
| 1337 | 168136 | bool acquire_tet(index_t t) { | |
| 1338 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 168136 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
168136 | geo_debug_assert(t < max_t()); |
| 1339 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 168136 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
168136 | geo_debug_assert(!owns_tet(t)); |
| 1340 | |||
| 1341 | 168136 | interfering_thread_ = cell_status_.acquire_cell( | |
| 1342 | 168136 | t, CellStatusArray::thread_index_t(id()) | |
| 1343 | ); | ||
| 1344 | |||
| 1345 |
1/2✓ Branch 0 taken 168136 times.
✗ Branch 1 not taken.
|
168136 | if(interfering_thread_ == NO_THREAD) { |
| 1346 |
3/10✓ Branch 0 taken 168136 times.
✗ Branch 1 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 168136 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 168136 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
|
168136 | geo_debug_assert(t == first_free_ || !tet_is_in_list(t)); |
| 1347 | #ifdef GEO_DEBUG | ||
| 1348 | 168136 | ++nb_acquired_tets_; | |
| 1349 | #endif | ||
| 1350 | 168136 | return true; | |
| 1351 | } | ||
| 1352 | ✗ | return false; | |
| 1353 | } | ||
| 1354 | |||
| 1355 | /** | ||
| 1356 | * \brief Releases a lock on a tetrahedron, making it | ||
| 1357 | * available to the other threads. | ||
| 1358 | */ | ||
| 1359 | 271027 | void release_tet(index_t t) { | |
| 1360 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 271027 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
271027 | geo_debug_assert(t < max_t()); |
| 1361 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 271027 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
271027 | geo_debug_assert(owns_tet(t)); |
| 1362 | #ifdef GEO_DEBUG | ||
| 1363 | 271027 | --nb_acquired_tets_; | |
| 1364 | #endif | ||
| 1365 | 271027 | cell_status_.release_cell(t); | |
| 1366 | 271027 | } | |
| 1367 | |||
| 1368 | |||
| 1369 | /** | ||
| 1370 | * \brief Tests whether this thread owns a tetrahedron. | ||
| 1371 | * \param[in] t index of the tetrahedron | ||
| 1372 | * \retval true if this thread owns t | ||
| 1373 | * \retval false otherwise | ||
| 1374 | */ | ||
| 1375 | 3689592 | bool owns_tet(index_t t) const { | |
| 1376 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 3689592 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
3689592 | geo_debug_assert(t < max_t()); |
| 1377 | return ( | ||
| 1378 | 3689592 | cell_status_.cell_thread(t) == | |
| 1379 | 3689592 | CellStatusArray::thread_index_t(id()) | |
| 1380 | 3689592 | ); | |
| 1381 | } | ||
| 1382 | |||
| 1383 | /** | ||
| 1384 | * \brief Finds the tetrahedron that (approximately) | ||
| 1385 | * contains a point using inexact predicates. | ||
| 1386 | * \details The result of this function can be used as a hint | ||
| 1387 | * for locate(). It accelerates locate as compared to calling | ||
| 1388 | * it directly. This technique is referred to as "structural | ||
| 1389 | * filtering". | ||
| 1390 | * \param[in] p a pointer to the coordinates of the point | ||
| 1391 | * \param[in] max_iter maximum number of traversed tets | ||
| 1392 | * \return the index of a tetrahedron that (approximately) | ||
| 1393 | * contains \p p. | ||
| 1394 | * If the point is outside the convex hull of | ||
| 1395 | * the inserted so-far points, then the returned tetrahedron | ||
| 1396 | * is a virtual one (first vertex is the "vertex at infinity" | ||
| 1397 | * of index -1) or NO_TETRAHEDRON if the virtual tetrahedra | ||
| 1398 | * were previously removed. | ||
| 1399 | */ | ||
| 1400 | 4038 | index_t locate_inexact( | |
| 1401 | const double* p, index_t hint, index_t max_iter | ||
| 1402 | ) const { | ||
| 1403 | // If no hint was specified, or the specified hint refers to a | ||
| 1404 | // tetrahedron that another thread freed/recycled in the meanwhile, | ||
| 1405 | // find a tetrahedron randomly. | ||
| 1406 |
4/6✓ Branch 0 taken 4033 times.
✓ Branch 1 taken 5 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 4033 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 4038 times.
|
4038 | if(hint != NO_TETRAHEDRON && tet_is_free(hint)) { |
| 1407 | ✗ | hint = NO_TETRAHEDRON; | |
| 1408 | } | ||
| 1409 |
2/2✓ Branch 0 taken 5 times.
✓ Branch 1 taken 4038 times.
|
4043 | while(hint == NO_TETRAHEDRON) { |
| 1410 | 5 | hint = thread_safe_random(max_used_t_); | |
| 1411 |
3/6✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 5 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 5 times.
|
5 | if(tet_is_free(hint) || tet_thread(hint) != NO_THREAD) { |
| 1412 | ✗ | hint = NO_TETRAHEDRON; | |
| 1413 | } | ||
| 1414 | } | ||
| 1415 | |||
| 1416 | // Always start from a real tet. If the tet is virtual, | ||
| 1417 | // find its real neighbor (always opposite to the | ||
| 1418 | // infinite vertex) | ||
| 1419 |
2/2✓ Branch 1 taken 663 times.
✓ Branch 2 taken 3375 times.
|
4038 | if(tet_is_virtual(hint)) { |
| 1420 |
1/2✓ Branch 0 taken 1941 times.
✗ Branch 1 not taken.
|
1941 | for(index_t lf = 0; lf < 4; ++lf) { |
| 1421 |
2/2✓ Branch 1 taken 663 times.
✓ Branch 2 taken 1278 times.
|
1941 | if(tet_vertex(hint, lf) == VERTEX_AT_INFINITY) { |
| 1422 | 663 | hint = tet_adjacent(hint, lf); | |
| 1423 | |||
| 1424 | // Yes, this can happen if the tetrahedron was | ||
| 1425 | // modified by another thread in the meanwhile. | ||
| 1426 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 663 times.
|
663 | if(hint == NO_TETRAHEDRON) { |
| 1427 | ✗ | return NO_TETRAHEDRON; | |
| 1428 | } | ||
| 1429 | |||
| 1430 | 663 | break; | |
| 1431 | } | ||
| 1432 | } | ||
| 1433 | } | ||
| 1434 | |||
| 1435 | 4038 | index_t t = hint; | |
| 1436 | 4038 | index_t t_pred = NO_TETRAHEDRON; | |
| 1437 | |||
| 1438 | 28605 | still_walking: | |
| 1439 | { | ||
| 1440 | |||
| 1441 | // Lookup the vertices of the current tetrahedron. | ||
| 1442 | const double* pv[4]; | ||
| 1443 |
2/2✓ Branch 0 taken 114420 times.
✓ Branch 1 taken 28605 times.
|
143025 | for(index_t lv=0; lv<4; ++lv) { |
| 1444 |
1/2✓ Branch 1 taken 114420 times.
✗ Branch 2 not taken.
|
114420 | index_t iv = tet_vertex(t,lv); |
| 1445 | |||
| 1446 | |||
| 1447 | // Since we did not acquire any lock, | ||
| 1448 | // it is possible that another threads made | ||
| 1449 | // this tetrahedron virtual (iv == NO_INDEX) or | ||
| 1450 | // deleted this tetrahedron (iv == VERTEX_OF_DELETED_TET) | ||
| 1451 | // (in both cases we exit immediately). | ||
| 1452 |
2/4✓ Branch 0 taken 114420 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 114420 times.
|
114420 | if(iv == NO_INDEX || iv == VERTEX_OF_DELETED_TET) { |
| 1453 | 931 | return NO_TETRAHEDRON; | |
| 1454 | } | ||
| 1455 |
1/2✓ Branch 1 taken 114420 times.
✗ Branch 2 not taken.
|
114420 | pv[lv] = vertex_ptr(iv); |
| 1456 | } | ||
| 1457 | |||
| 1458 |
2/2✓ Branch 0 taken 72628 times.
✓ Branch 1 taken 3107 times.
|
75735 | for(index_t f = 0; f < 4; ++f) { |
| 1459 | |||
| 1460 |
1/2✓ Branch 1 taken 72628 times.
✗ Branch 2 not taken.
|
72628 | index_t t_next = tet_adjacent(t,f); |
| 1461 | |||
| 1462 | // If the opposite tet is -1, then it means that | ||
| 1463 | // we are trying to locate() (e.g. called from | ||
| 1464 | // nearest_vertex) within a tetrahedralization | ||
| 1465 | // from which the infinite tets were removed. | ||
| 1466 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 72628 times.
|
72628 | if(t_next == NO_INDEX) { |
| 1467 | ✗ | return NO_TETRAHEDRON; | |
| 1468 | } | ||
| 1469 | |||
| 1470 | // If the candidate next tetrahedron is the | ||
| 1471 | // one we came from, then we know already that | ||
| 1472 | // the orientation is positive, thus we examine | ||
| 1473 | // the next candidate (or exit the loop if they | ||
| 1474 | // are exhausted). | ||
| 1475 |
2/2✓ Branch 0 taken 12246 times.
✓ Branch 1 taken 60382 times.
|
72628 | if(t_next == t_pred) { |
| 1476 | 12246 | continue ; | |
| 1477 | } | ||
| 1478 | |||
| 1479 | // To test the orientation of p w.r.t. the facet f of | ||
| 1480 | // t, we replace vertex number f with p in t (same | ||
| 1481 | // convention as in CGAL). | ||
| 1482 | 60382 | const double* pv_bkp = pv[f]; | |
| 1483 | 60382 | pv[f] = p; | |
| 1484 |
1/2✓ Branch 1 taken 60382 times.
✗ Branch 2 not taken.
|
60382 | Sign ori = PCK::orient_3d_inexact( |
| 1485 | pv[0], pv[1], pv[2], pv[3] | ||
| 1486 | ); | ||
| 1487 | |||
| 1488 | // If the orientation is not negative, then we cannot | ||
| 1489 | // walk towards t_next, and examine the next candidate | ||
| 1490 | // (or exit the loop if they are exhausted). | ||
| 1491 |
2/2✓ Branch 0 taken 34884 times.
✓ Branch 1 taken 25498 times.
|
60382 | if(ori != NEGATIVE) { |
| 1492 | 34884 | pv[f] = pv_bkp; | |
| 1493 | 34884 | continue; | |
| 1494 | } | ||
| 1495 | |||
| 1496 | // If the opposite tet is a virtual tet, then | ||
| 1497 | // the point has a positive orientation relative | ||
| 1498 | // to the facet on the border of the convex hull, | ||
| 1499 | // thus t_next is a tet in conflict and we are | ||
| 1500 | // done. | ||
| 1501 |
3/4✓ Branch 1 taken 25498 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 931 times.
✓ Branch 4 taken 24567 times.
|
25498 | if(tet_is_virtual(t_next)) { |
| 1502 | 931 | return t_next; | |
| 1503 | } | ||
| 1504 | |||
| 1505 | // If we reach this point, then t_next is a valid | ||
| 1506 | // successor, thus we are still walking. | ||
| 1507 | 24567 | t_pred = t; | |
| 1508 | 24567 | t = t_next; | |
| 1509 |
1/2✓ Branch 0 taken 24567 times.
✗ Branch 1 not taken.
|
24567 | if(--max_iter != 0) { |
| 1510 | 24567 | goto still_walking; | |
| 1511 | } | ||
| 1512 | } | ||
| 1513 | } | ||
| 1514 | |||
| 1515 | // If we reach this point, we did not find a valid successor | ||
| 1516 | // for walking (a face for which p has negative orientation), | ||
| 1517 | // thus we reached the tet for which p has all positive | ||
| 1518 | // face orientations (i.e. the tet that contains p). | ||
| 1519 | |||
| 1520 | 3107 | return t; | |
| 1521 | } | ||
| 1522 | |||
| 1523 | |||
| 1524 | /** | ||
| 1525 | * \brief Tests whether a tetrahedron is | ||
| 1526 | * a virtual one. | ||
| 1527 | * \details Virtual tetrahedra are tetrahedra | ||
| 1528 | * incident to the vertex at infinity. | ||
| 1529 | * \param[in] t index of the tetrahedron | ||
| 1530 | * \retval true if tetrahedron \p t is virtual | ||
| 1531 | * \retval false otherwise | ||
| 1532 | */ | ||
| 1533 | 75160 | bool tet_is_virtual(index_t t) const { | |
| 1534 | return | ||
| 1535 |
1/2✓ Branch 1 taken 75160 times.
✗ Branch 2 not taken.
|
150320 | !tet_is_free(t) && ( |
| 1536 |
2/2✓ Branch 1 taken 75143 times.
✓ Branch 2 taken 17 times.
|
75160 | cell_to_v_store_[4 * t] == VERTEX_AT_INFINITY || |
| 1537 |
2/2✓ Branch 1 taken 74078 times.
✓ Branch 2 taken 1065 times.
|
75143 | cell_to_v_store_[4 * t + 1] == VERTEX_AT_INFINITY || |
| 1538 |
2/2✓ Branch 1 taken 73521 times.
✓ Branch 2 taken 557 times.
|
74078 | cell_to_v_store_[4 * t + 2] == VERTEX_AT_INFINITY || |
| 1539 |
2/2✓ Branch 1 taken 886 times.
✓ Branch 2 taken 72635 times.
|
148681 | cell_to_v_store_[4 * t + 3] == VERTEX_AT_INFINITY) ; |
| 1540 | } | ||
| 1541 | |||
| 1542 | |||
| 1543 | /** | ||
| 1544 | * \brief Returns the local index of a vertex by | ||
| 1545 | * facet and by local vertex index in the facet. | ||
| 1546 | * \details | ||
| 1547 | * tet facet vertex is such that the tetrahedron | ||
| 1548 | * formed with: | ||
| 1549 | * - vertex lv | ||
| 1550 | * - tet_facet_vertex(lv,0) | ||
| 1551 | * - tet_facet_vertex(lv,1) | ||
| 1552 | * - tet_facet_vertex(lv,2) | ||
| 1553 | * has the same orientation as the original tetrahedron for | ||
| 1554 | * any vertex lv. | ||
| 1555 | * \param[in] f local facet index, in (0,1,2,3) | ||
| 1556 | * \param[in] v local vertex index, in (0,1,2) | ||
| 1557 | * \return the local tetrahedron vertex index of | ||
| 1558 | * vertex \p v in facet \p f | ||
| 1559 | */ | ||
| 1560 | 366054 | static index_t tet_facet_vertex(index_t f, index_t v) { | |
| 1561 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 366054 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
366054 | geo_debug_assert(f < 4); |
| 1562 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 366054 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
366054 | geo_debug_assert(v < 3); |
| 1563 | 366054 | return index_t(tet_facet_vertex_[f][v]); | |
| 1564 | } | ||
| 1565 | |||
| 1566 | /** | ||
| 1567 | * \brief Gets the index of a vertex of a tetrahedron | ||
| 1568 | * \param[in] t index of the tetrahedron | ||
| 1569 | * \param[in] lv local vertex (0,1,2 or 3) index in \p t | ||
| 1570 | * \return the global index of the \p lv%th vertex of tetrahedron \p t | ||
| 1571 | * or -1 if the vertex is at infinity | ||
| 1572 | */ | ||
| 1573 | 3796345 | index_t tet_vertex(index_t t, index_t lv) const { | |
| 1574 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 3796345 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
3796345 | geo_debug_assert(t < max_t()); |
| 1575 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 3796345 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
3796345 | geo_debug_assert(lv < 4); |
| 1576 | 3796345 | return cell_to_v_store_[4 * t + lv]; | |
| 1577 | } | ||
| 1578 | |||
| 1579 | /** | ||
| 1580 | * \brief Finds the index of the vertex in a tetrahedron. | ||
| 1581 | * \param[in] t the tetrahedron | ||
| 1582 | * \param[in] v the vertex | ||
| 1583 | * \return iv such that tet_vertex(t,v)==iv | ||
| 1584 | * \pre \p t is incident to \p v | ||
| 1585 | */ | ||
| 1586 | 1500 | index_t find_tet_vertex(index_t t, index_t v) const { | |
| 1587 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 1500 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
1500 | geo_debug_assert(t < max_t()); |
| 1588 | // Find local index of v in tetrahedron t vertices. | ||
| 1589 | 1500 | const index_t* T = &(cell_to_v_store_[4 * t]); | |
| 1590 | 1500 | return find_4(T,v); | |
| 1591 | } | ||
| 1592 | |||
| 1593 | |||
| 1594 | /** | ||
| 1595 | * \brief Gets the index of a vertex of a tetrahedron | ||
| 1596 | * \param[in] t index of the tetrahedron | ||
| 1597 | * \param[in] lv local vertex (0,1,2 or 3) index in \p t | ||
| 1598 | * \return the global index of the \p lv%th vertex of tetrahedron \p t | ||
| 1599 | * \pre Vertex \p lv of tetrahedron \p t is not at infinity | ||
| 1600 | */ | ||
| 1601 | 28904 | index_t finite_tet_vertex(index_t t, index_t lv) const { | |
| 1602 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 28904 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
28904 | geo_debug_assert(t < max_t()); |
| 1603 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 28904 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
28904 | geo_debug_assert(lv < 4); |
| 1604 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 28904 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
28904 | geo_debug_assert(cell_to_v_store_[4 * t + lv] != NO_INDEX); |
| 1605 | 28904 | return cell_to_v_store_[4 * t + lv]; | |
| 1606 | } | ||
| 1607 | |||
| 1608 | /** | ||
| 1609 | * \brief Sets a tetrahedron-to-vertex adjacency. | ||
| 1610 | * \param[in] t index of the tetrahedron | ||
| 1611 | * \param[in] lv local vertex index (0,1,2 or 3) in \p t | ||
| 1612 | * \param[in] v global index of the vertex | ||
| 1613 | */ | ||
| 1614 | 311388 | void set_tet_vertex(index_t t, index_t lv, index_t v) { | |
| 1615 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 311388 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
311388 | geo_debug_assert(t < max_t()); |
| 1616 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 311388 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
311388 | geo_debug_assert(lv < 4); |
| 1617 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 311388 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
311388 | geo_debug_assert(owns_tet(t)); |
| 1618 | 311388 | cell_to_v_store_[4 * t + lv] = v; | |
| 1619 | 311388 | } | |
| 1620 | |||
| 1621 | /** | ||
| 1622 | * \brief Gets the index of a tetrahedron adjacent to another one. | ||
| 1623 | * \param[in] t index of the tetrahedron | ||
| 1624 | * \param[in] lf local facet (0,1,2 or 3) index in \p t | ||
| 1625 | * \return the tetrahedron adjacent to \p t accorss facet \p lf | ||
| 1626 | */ | ||
| 1627 | 1711871 | index_t tet_adjacent(index_t t, index_t lf) const { | |
| 1628 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 1711871 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
1711871 | geo_debug_assert(t < max_t()); |
| 1629 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 1711871 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
1711871 | geo_debug_assert(lf < 4); |
| 1630 | 1711871 | index_t result = cell_to_cell_store_[4 * t + lf]; | |
| 1631 | 1711871 | return result; | |
| 1632 | } | ||
| 1633 | |||
| 1634 | /** | ||
| 1635 | * \brief Sets a tetrahedron-to-tetrahedron adjacency. | ||
| 1636 | * \param[in] t1 index of the first tetrahedron | ||
| 1637 | * \param[in] lf1 local facet index (0,1,2 or 3) in t1 | ||
| 1638 | * \param[in] t2 index of the tetrahedron | ||
| 1639 | * adjacent to \p t1 accros \p lf1 | ||
| 1640 | */ | ||
| 1641 | 514430 | void set_tet_adjacent(index_t t1, index_t lf1, index_t t2) { | |
| 1642 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 514430 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
514430 | geo_debug_assert(t1 < max_t()); |
| 1643 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 514430 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
514430 | geo_debug_assert(t2 < max_t()); |
| 1644 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 514430 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
514430 | geo_debug_assert(lf1 < 4); |
| 1645 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 514430 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
514430 | geo_debug_assert(owns_tet(t1)); |
| 1646 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 514430 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
514430 | geo_debug_assert(owns_tet(t2)); |
| 1647 | 514430 | cell_to_cell_store_[4 * t1 + lf1] = t2; | |
| 1648 | 514430 | } | |
| 1649 | |||
| 1650 | /** | ||
| 1651 | * \brief Finds the index of the facet accros which t1 is | ||
| 1652 | * adjacent to t2. | ||
| 1653 | * \param[in] t1 first tetrahedron | ||
| 1654 | * \param[in] t2 second tetrahedron | ||
| 1655 | * \return f such that tet_adjacent(t1,f)==t2 | ||
| 1656 | * \pre \p t1 and \p t2 are adjacent | ||
| 1657 | */ | ||
| 1658 | 102866 | index_t find_tet_adjacent(index_t t1, index_t t2) const { | |
| 1659 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 102866 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
102866 | geo_debug_assert(t1 < max_t()); |
| 1660 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 102866 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
102866 | geo_debug_assert(t2 < max_t()); |
| 1661 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 102866 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
102866 | geo_debug_assert(t1 != t2); |
| 1662 | |||
| 1663 | // Find local index of t2 in tetrahedron t1 adajcent tets. | ||
| 1664 | 102866 | const index_t* T = &(cell_to_cell_store_[4 * t1]); | |
| 1665 | 102866 | index_t result = find_4(T,t2); | |
| 1666 | |||
| 1667 | // Sanity check: make sure that t1 is adjacent to t2 | ||
| 1668 | // only once! | ||
| 1669 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 102866 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
102866 | geo_debug_assert(tet_adjacent(t1,(result+1)%4) != t2); |
| 1670 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 102866 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
102866 | geo_debug_assert(tet_adjacent(t1,(result+2)%4) != t2); |
| 1671 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 102866 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
102866 | geo_debug_assert(tet_adjacent(t1,(result+3)%4) != t2); |
| 1672 | 102866 | return result; | |
| 1673 | } | ||
| 1674 | |||
| 1675 | |||
| 1676 | /** | ||
| 1677 | * Gets the local facet index incident to an | ||
| 1678 | * oriented halfedge. | ||
| 1679 | * \param[in] t index of the tetrahedron | ||
| 1680 | * \param[in] v1 global index of the first extremity | ||
| 1681 | * \param[in] v2 global index of the second extremity | ||
| 1682 | * \return the local index of the facet incident to | ||
| 1683 | * the oriented edge \p v1, \p v2. | ||
| 1684 | */ | ||
| 1685 | 2286 | index_t get_facet_by_halfedge( | |
| 1686 | index_t t, index_t v1, index_t v2 | ||
| 1687 | ) const { | ||
| 1688 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 2286 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
2286 | geo_debug_assert(t < max_t()); |
| 1689 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 2286 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
2286 | geo_debug_assert(v1 != v2); |
| 1690 | // Find local index of v1 and v2 in tetrahedron t | ||
| 1691 | 2286 | const index_t* T = &(cell_to_v_store_[4 * t]); | |
| 1692 | 2286 | index_t lv1 = find_4(T,v1); | |
| 1693 | 2286 | index_t lv2 = find_4(T,v2); | |
| 1694 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 2286 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
2286 | geo_debug_assert(lv1 != lv2); |
| 1695 | 2286 | return index_t(halfedge_facet_[lv1][lv2]); | |
| 1696 | } | ||
| 1697 | |||
| 1698 | |||
| 1699 | /** | ||
| 1700 | * Gets the local facet indices incident to an | ||
| 1701 | * oriented halfedge. | ||
| 1702 | * \param[in] t index of the tetrahedron | ||
| 1703 | * \param[in] v1 global index of the first extremity | ||
| 1704 | * \param[in] v2 global index of the second extremity | ||
| 1705 | * \param[out] f12 the local index of the facet | ||
| 1706 | * indicent to the halfedge [v1,v2] | ||
| 1707 | * \param[out] f21 the local index of the facet | ||
| 1708 | * indicent to the halfedge [v2,v1] | ||
| 1709 | */ | ||
| 1710 | 1500 | void get_facets_by_halfedge( | |
| 1711 | index_t t, index_t v1, index_t v2, | ||
| 1712 | index_t& f12, index_t& f21 | ||
| 1713 | ) const { | ||
| 1714 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 1500 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
1500 | geo_debug_assert(t < max_t()); |
| 1715 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 1500 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
1500 | geo_debug_assert(v1 != v2); |
| 1716 | |||
| 1717 | // Find local index of v1 and v2 in tetrahedron t | ||
| 1718 | // The following expression is 10% faster than using | ||
| 1719 | // if() statements (multiply by boolean result of test). | ||
| 1720 | // Thank to Laurent Alonso for this idea. | ||
| 1721 | 1500 | const index_t* T = &(cell_to_v_store_[4 * t]); | |
| 1722 | |||
| 1723 | 1500 | index_t lv1 = | |
| 1724 |
4/4✓ Branch 0 taken 366 times.
✓ Branch 1 taken 1134 times.
✓ Branch 2 taken 408 times.
✓ Branch 3 taken 1092 times.
|
1500 | index_t((T[1] == v1) | ((T[2] == v1) * 2) | ((T[3] == v1) * 3)); |
| 1725 | |||
| 1726 | 1500 | index_t lv2 = | |
| 1727 |
4/4✓ Branch 0 taken 394 times.
✓ Branch 1 taken 1106 times.
✓ Branch 2 taken 276 times.
✓ Branch 3 taken 1224 times.
|
1500 | index_t((T[1] == v2) | ((T[2] == v2) * 2) | ((T[3] == v2) * 3)); |
| 1728 | |||
| 1729 |
3/8✓ Branch 0 taken 306 times.
✓ Branch 1 taken 1194 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 306 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
1500 | geo_debug_assert(lv1 != 0 || T[0] == v1); |
| 1730 |
3/8✓ Branch 0 taken 566 times.
✓ Branch 1 taken 934 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 566 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
1500 | geo_debug_assert(lv2 != 0 || T[0] == v2); |
| 1731 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 1500 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
1500 | geo_debug_assert(lv1 != NO_INDEX); |
| 1732 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 1500 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
1500 | geo_debug_assert(lv2 != NO_INDEX); |
| 1733 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 1500 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
1500 | geo_debug_assert(lv1 != lv2); |
| 1734 | |||
| 1735 | 1500 | f12 = index_t(halfedge_facet_[lv1][lv2]); | |
| 1736 | 1500 | f21 = index_t(halfedge_facet_[lv2][lv1]); | |
| 1737 | 1500 | } | |
| 1738 | |||
| 1739 | /** | ||
| 1740 | * \brief Symbolic value of the cell_next_ field | ||
| 1741 | * that indicates the end of list in a linked | ||
| 1742 | * list of tetrahedra. | ||
| 1743 | */ | ||
| 1744 | static constexpr index_t END_OF_LIST = NO_INDEX; | ||
| 1745 | |||
| 1746 | |||
| 1747 | /** | ||
| 1748 | * \brief Symbolic value of the cell_next_ field | ||
| 1749 | * for a tetrahedron that is not in a list. | ||
| 1750 | */ | ||
| 1751 | static constexpr index_t NOT_IN_LIST = index_t(-2); | ||
| 1752 | |||
| 1753 | /** | ||
| 1754 | * \brief Symbolic value for t2v_[] indicating a deleted tetrahedron. | ||
| 1755 | */ | ||
| 1756 | static constexpr index_t VERTEX_OF_DELETED_TET = index_t(-2); | ||
| 1757 | |||
| 1758 | /** | ||
| 1759 | * \brief Gets the number of vertices. | ||
| 1760 | * \return the number of vertices in this Delaunay | ||
| 1761 | */ | ||
| 1762 | 4185724 | index_t nb_vertices() const { | |
| 1763 | 4185724 | return nb_vertices_; | |
| 1764 | } | ||
| 1765 | |||
| 1766 | /** | ||
| 1767 | * \brief Gets a pointer to a vertex by its global index. | ||
| 1768 | * \param[in] i global index of the vertex | ||
| 1769 | * \return a pointer to vertex \p i | ||
| 1770 | */ | ||
| 1771 | 3684480 | const double* vertex_ptr(index_t i) const { | |
| 1772 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 3684480 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
3684480 | geo_debug_assert(i < nb_vertices()); |
| 1773 | 3684480 | return vertices_ + vertex_stride_ * i; | |
| 1774 | } | ||
| 1775 | |||
| 1776 | /** | ||
| 1777 | * \brief Tests whether a tetrahedron belongs to a linked | ||
| 1778 | * list. | ||
| 1779 | * \details Tetrahedra can be linked, it is used to manage | ||
| 1780 | * both the free list that recycles deleted tetrahedra, | ||
| 1781 | * the conflict region and the list of newly created | ||
| 1782 | * tetrahedra. | ||
| 1783 | * \param[in] t the index of the tetrahedron | ||
| 1784 | * \retval true if tetrahedron \p t belongs to a linked list | ||
| 1785 | * \retval false otherwise | ||
| 1786 | */ | ||
| 1787 | 568204 | bool tet_is_in_list(index_t t) const { | |
| 1788 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 568204 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
568204 | geo_debug_assert(t < max_t()); |
| 1789 | 568204 | return (cell_next_[t] != NOT_IN_LIST); | |
| 1790 | } | ||
| 1791 | |||
| 1792 | /** | ||
| 1793 | * \brief Gets the index of a successor of a tetrahedron. | ||
| 1794 | * \details Tetrahedra can be linked, it is used to manage | ||
| 1795 | * both the free list that recycles deleted tetrahedra. | ||
| 1796 | * \param[in] t the index of the tetrahedron | ||
| 1797 | * \retval END_OF_LIST if the end of the list is reached | ||
| 1798 | * \retval the index of the successor of | ||
| 1799 | * tetrahedron \t otherwise | ||
| 1800 | * \pre tet_is_in_list(t) | ||
| 1801 | */ | ||
| 1802 | 102891 | index_t tet_next(index_t t) const { | |
| 1803 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 102891 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
102891 | geo_debug_assert(t < max_t()); |
| 1804 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 102891 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
102891 | geo_debug_assert(tet_is_in_list(t)); |
| 1805 | 102891 | return cell_next_[t]; | |
| 1806 | } | ||
| 1807 | |||
| 1808 | |||
| 1809 | 5 | index_t tet_thread(index_t t) const { | |
| 1810 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 5 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
5 | geo_debug_assert(t < max_t()); |
| 1811 | 5 | return cell_status_.cell_thread(t); | |
| 1812 | } | ||
| 1813 | |||
| 1814 | /** | ||
| 1815 | * \brief Adds a tetrahedron to a linked list. | ||
| 1816 | * \details Tetrahedra can be linked, it is used to manage | ||
| 1817 | * the free list that recycles deleted tetrahedra. | ||
| 1818 | * \param[in] t the index of the tetrahedron | ||
| 1819 | * \param[in,out] first first item of the list or END_OF_LIST if | ||
| 1820 | * the list is empty | ||
| 1821 | * \param[in,out] last last item of the list or END_OF_LIST if | ||
| 1822 | * the list is empty | ||
| 1823 | */ | ||
| 1824 | void add_tet_to_list(index_t t, index_t& first, index_t& last) { | ||
| 1825 | geo_debug_assert(t < max_t()); | ||
| 1826 | geo_debug_assert(!tet_is_in_list(t)); | ||
| 1827 | geo_debug_assert(owns_tet(t)); | ||
| 1828 | if(last == END_OF_LIST) { | ||
| 1829 | geo_debug_assert(first == END_OF_LIST); | ||
| 1830 | first = last = t; | ||
| 1831 | cell_next_[t] = END_OF_LIST; | ||
| 1832 | } else { | ||
| 1833 | cell_next_[t] = first; | ||
| 1834 | first = t; | ||
| 1835 | } | ||
| 1836 | } | ||
| 1837 | |||
| 1838 | /** | ||
| 1839 | * \brief Removes a tetrahedron from the linked list it | ||
| 1840 | * belongs to. | ||
| 1841 | * \details Tetrahedra can be linked, it is used to manage | ||
| 1842 | * the free list that recycles deleted tetrahedra. | ||
| 1843 | * \param[in] t the index of the tetrahedron | ||
| 1844 | */ | ||
| 1845 | 102891 | void remove_tet_from_list(index_t t) { | |
| 1846 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 102891 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
102891 | geo_debug_assert(t < max_t()); |
| 1847 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 102891 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
102891 | geo_debug_assert(tet_is_in_list(t)); |
| 1848 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 102891 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
102891 | geo_debug_assert(owns_tet(t)); |
| 1849 | 102891 | cell_next_[t] = NOT_IN_LIST; | |
| 1850 | 102891 | } | |
| 1851 | |||
| 1852 | |||
| 1853 | /** | ||
| 1854 | * \brief Creates a new tetrahedron. | ||
| 1855 | * \details Uses either a tetrahedron recycled | ||
| 1856 | * from the free list, or creates a new one by | ||
| 1857 | * expanding the two indices arrays. | ||
| 1858 | * \return the index of the newly created tetrahedron | ||
| 1859 | */ | ||
| 1860 | 102891 | index_t new_tetrahedron() { | |
| 1861 | |||
| 1862 | // If the memory pool is full, then we expand it. | ||
| 1863 | // This cannot be done when running multiple threads. | ||
| 1864 |
2/2✓ Branch 0 taken 18301 times.
✓ Branch 1 taken 84590 times.
|
102891 | if(first_free_ == END_OF_LIST) { |
| 1865 |
2/8✓ Branch 1 taken 18301 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 18301 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
18301 | geo_debug_assert(!Process::is_running_threads()); |
| 1866 | 36602 | master_->cell_to_v_store_.resize( | |
| 1867 |
1/2✓ Branch 2 taken 18301 times.
✗ Branch 3 not taken.
|
18301 | master_->cell_to_v_store_.size() + 4, NO_INDEX |
| 1868 | ); | ||
| 1869 | 36602 | master_->cell_to_cell_store_.resize( | |
| 1870 |
1/2✓ Branch 2 taken 18301 times.
✗ Branch 3 not taken.
|
18301 | master_->cell_to_cell_store_.size() + 4, NO_INDEX |
| 1871 | ); | ||
| 1872 |
1/2✓ Branch 1 taken 18301 times.
✗ Branch 2 not taken.
|
18301 | master_->cell_next_.push_back(END_OF_LIST); |
| 1873 |
1/2✓ Branch 1 taken 18301 times.
✗ Branch 2 not taken.
|
18301 | master_->cell_status_.grow(); |
| 1874 | 18301 | ++nb_free_; | |
| 1875 | 18301 | ++max_t_; | |
| 1876 | 18301 | first_free_ = master_->cell_status_.size() - 1; | |
| 1877 | } | ||
| 1878 | |||
| 1879 |
1/2✓ Branch 1 taken 102891 times.
✗ Branch 2 not taken.
|
102891 | acquire_and_mark_tet_as_created(first_free_); |
| 1880 | 102891 | index_t result = first_free_; | |
| 1881 | |||
| 1882 |
1/2✓ Branch 1 taken 102891 times.
✗ Branch 2 not taken.
|
102891 | first_free_ = tet_next(first_free_); |
| 1883 |
1/2✓ Branch 1 taken 102891 times.
✗ Branch 2 not taken.
|
102891 | remove_tet_from_list(result); |
| 1884 | |||
| 1885 |
1/2✓ Branch 1 taken 102891 times.
✗ Branch 2 not taken.
|
102891 | cell_to_cell_store_[4 * result] = NO_INDEX; |
| 1886 |
1/2✓ Branch 1 taken 102891 times.
✗ Branch 2 not taken.
|
102891 | cell_to_cell_store_[4 * result + 1] = NO_INDEX; |
| 1887 |
1/2✓ Branch 1 taken 102891 times.
✗ Branch 2 not taken.
|
102891 | cell_to_cell_store_[4 * result + 2] = NO_INDEX; |
| 1888 |
1/2✓ Branch 1 taken 102891 times.
✗ Branch 2 not taken.
|
102891 | cell_to_cell_store_[4 * result + 3] = NO_INDEX; |
| 1889 | |||
| 1890 | 102891 | max_used_t_ = std::max(max_used_t_, result); | |
| 1891 | |||
| 1892 | 102891 | --nb_free_; | |
| 1893 | 102891 | return result; | |
| 1894 | } | ||
| 1895 | |||
| 1896 | /** | ||
| 1897 | * \brief Creates a new tetrahedron. | ||
| 1898 | * \details Sets the vertices. Adjacent tetrahedra index are | ||
| 1899 | * left uninitialized. Uses either a tetrahedron recycled | ||
| 1900 | * from the free list, or creates a new one by | ||
| 1901 | * expanding the two indices arrays. | ||
| 1902 | * \param[in] v1 index of the first vertex | ||
| 1903 | * \param[in] v2 index of the second vertex | ||
| 1904 | * \param[in] v3 index of the third vertex | ||
| 1905 | * \param[in] v4 index of the fourth vertex | ||
| 1906 | * \return the index of the newly created tetrahedron | ||
| 1907 | */ | ||
| 1908 | 102891 | index_t new_tetrahedron(index_t v1, index_t v2, index_t v3, index_t v4) { | |
| 1909 | 102891 | index_t result = new_tetrahedron(); | |
| 1910 | 102891 | cell_to_v_store_[4 * result] = v1; | |
| 1911 | 102891 | cell_to_v_store_[4 * result + 1] = v2; | |
| 1912 | 102891 | cell_to_v_store_[4 * result + 2] = v3; | |
| 1913 | 102891 | cell_to_v_store_[4 * result + 3] = v4; | |
| 1914 | 102891 | return result; | |
| 1915 | } | ||
| 1916 | |||
| 1917 | /** | ||
| 1918 | * \brief Finds the index of an integer in an array of four integers. | ||
| 1919 | * \param[in] T a const pointer to an array of four integers | ||
| 1920 | * \param[in] v the integer to retrieve in \p T | ||
| 1921 | * \return the index (0,1,2 or 3) of \p v in \p T | ||
| 1922 | * \pre The four entries of \p T are different and one of them is | ||
| 1923 | * equal to \p v. | ||
| 1924 | */ | ||
| 1925 | 108938 | static index_t find_4(const index_t* T, index_t v) { | |
| 1926 | // The following expression is 10% faster than using | ||
| 1927 | // if() statements. This uses the C++ norm, that | ||
| 1928 | // ensures that the 'true' boolean value converted to | ||
| 1929 | // an int is always 1. With most compilers, this avoids | ||
| 1930 | // generating branching instructions. | ||
| 1931 | // Thank to Laurent Alonso for this idea. | ||
| 1932 | // Note: Laurent also has this version: | ||
| 1933 | // (T[0] != v)+(T[2]==v)+2*(T[3]==v) | ||
| 1934 | // that avoids a *3 multiply, but it is not faster in | ||
| 1935 | // practice. | ||
| 1936 | 108938 | index_t result = index_t( | |
| 1937 |
4/4✓ Branch 0 taken 30631 times.
✓ Branch 1 taken 78307 times.
✓ Branch 2 taken 29481 times.
✓ Branch 3 taken 79457 times.
|
108938 | (T[1] == v) | ((T[2] == v) * 2) | ((T[3] == v) * 3) |
| 1938 | ); | ||
| 1939 | // Sanity check, important if it was T[0], not explicitly | ||
| 1940 | // tested (detects input that does not meet the precondition). | ||
| 1941 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 108938 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
108938 | geo_debug_assert(T[result] == v); |
| 1942 | 108938 | return result; | |
| 1943 | } | ||
| 1944 | |||
| 1945 | /** | ||
| 1946 | * \brief Wakes up all the threads that are waiting for | ||
| 1947 | * this thread. | ||
| 1948 | */ | ||
| 1949 | 4068 | void send_event() { | |
| 1950 | 4068 | cond_.notify_all(); | |
| 1951 | 4068 | } | |
| 1952 | |||
| 1953 | /** | ||
| 1954 | * \brief Waits for a thread. | ||
| 1955 | * \details Sleeps until thread \p t calls send_event(). | ||
| 1956 | * \param[in] t index of the thread | ||
| 1957 | * \pre t < nb_threads() | ||
| 1958 | */ | ||
| 1959 | ✗ | void wait_for_event(index_t t) { | |
| 1960 | // Fixed by Hiep Vu: enlarged critical section (contains | ||
| 1961 | // now the test (!thrd->finished) | ||
| 1962 | ✗ | Delaunay3dThread* thrd = thread(t); | |
| 1963 | // RAII: ctor locks, dtor unlocks | ||
| 1964 | ✗ | std::unique_lock<std::mutex> L(thrd->mutex_); | |
| 1965 | ✗ | if(!thrd->finished_) { | |
| 1966 | ✗ | thrd->cond_.wait(L); | |
| 1967 | } | ||
| 1968 | ✗ | } | |
| 1969 | |||
| 1970 | /****** iterative stellate_conflict_zone *****************/ | ||
| 1971 | |||
| 1972 | /** | ||
| 1973 | * \brief Used to represent the stack in the | ||
| 1974 | * (de-recursified) stellate_conflict_zone_iterative() | ||
| 1975 | * function. | ||
| 1976 | */ | ||
| 1977 | class StellateConflictStack { | ||
| 1978 | public: | ||
| 1979 | |||
| 1980 | /** | ||
| 1981 | * \brief Pushes a new frame onto the stack. | ||
| 1982 | * \details This also creates the local variables (they are | ||
| 1983 | * left uninitialized). | ||
| 1984 | * \param[in] t1 index of a tetrahedron on the border of | ||
| 1985 | * the conflict zone | ||
| 1986 | * \param[in] t1fbord index of the facet of \p t1 that is | ||
| 1987 | * on the border of the conflict zone | ||
| 1988 | * \param[in] t1fprev index of the facet of \p t1 that we | ||
| 1989 | * come from, or NO_INDEX if \p t1 is the first tetrahedron | ||
| 1990 | */ | ||
| 1991 | 1000 | void push(index_t t1, index_t t1fbord, index_t t1fprev) { | |
| 1992 | 1000 | store_.resize(store_.size()+1); | |
| 1993 | 1000 | top().t1 = t1; | |
| 1994 | 1000 | top().t1fbord = Numeric::uint8(t1fbord); | |
| 1995 | 1000 | top().t1fprev = Numeric::uint8(t1fprev); | |
| 1996 | 1000 | } | |
| 1997 | |||
| 1998 | /** | ||
| 1999 | * \brief Saves local variables into the current stack frame. | ||
| 2000 | * \param[in] new_t the index of the newly created tetrahedron | ||
| 2001 | * \param[in] t1ft2 the facet of t1 that is adjacent to t2 | ||
| 2002 | * \param[in] t2ft1 the facet of t2 that is adjacent to t1 | ||
| 2003 | */ | ||
| 2004 | 994 | void save_locals(index_t new_t, index_t t1ft2, index_t t2ft1) { | |
| 2005 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 994 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
994 | geo_debug_assert(!empty()); |
| 2006 | 994 | top().new_t = new_t; | |
| 2007 | 994 | top().t1ft2 = Numeric::uint8(t1ft2); | |
| 2008 | 994 | top().t2ft1 = Numeric::uint8(t2ft1); | |
| 2009 | 994 | } | |
| 2010 | |||
| 2011 | /** | ||
| 2012 | * \brief Gets the parameters from the current stack frame. | ||
| 2013 | * \param[out] t1 index of a tetrahedron on the border of | ||
| 2014 | * the conflict zone | ||
| 2015 | * \param[out] t1fbord index of the facet of \p t1 that is | ||
| 2016 | * on the border of the conflict zone | ||
| 2017 | * \param[out] t1fprev index of the facet of \p t1 that we | ||
| 2018 | * come from, or NO_INDEX if \p t1 is the first tetrahedron | ||
| 2019 | */ | ||
| 2020 | 1994 | void get_parameters( | |
| 2021 | index_t& t1, index_t& t1fbord, index_t& t1fprev | ||
| 2022 | ) const { | ||
| 2023 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 1994 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
1994 | geo_debug_assert(!empty()); |
| 2024 | 1994 | t1 = top().t1; | |
| 2025 | 1994 | t1fbord = index_t(top().t1fbord); | |
| 2026 | 1994 | t1fprev = index_t(top().t1fprev); | |
| 2027 | 1994 | } | |
| 2028 | |||
| 2029 | |||
| 2030 | /** | ||
| 2031 | * \brief Gets the local variables from the current stack frame. | ||
| 2032 | * \param[out] new_t the index of the newly created tetrahedron | ||
| 2033 | * \param[out] t1ft2 the facet of t1 that is adjacent to t2 | ||
| 2034 | * \param[out] t2ft1 the facet of t2 that is adjacent to t1 | ||
| 2035 | */ | ||
| 2036 | 994 | void get_locals( | |
| 2037 | index_t& new_t, index_t& t1ft2, index_t& t2ft1 | ||
| 2038 | ) const { | ||
| 2039 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 994 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
994 | geo_debug_assert(!empty()); |
| 2040 | 994 | new_t = top().new_t; | |
| 2041 | 994 | t1ft2 = index_t(top().t1ft2); | |
| 2042 | 994 | t2ft1 = index_t(top().t2ft1); | |
| 2043 | 994 | } | |
| 2044 | |||
| 2045 | /** | ||
| 2046 | * \brief Pops a stack frame. | ||
| 2047 | */ | ||
| 2048 | 1000 | void pop() { | |
| 2049 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 1000 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
1000 | geo_debug_assert(!empty()); |
| 2050 | 1000 | store_.pop_back(); | |
| 2051 | 1000 | } | |
| 2052 | |||
| 2053 | /** | ||
| 2054 | * \brief Tests whether the stack is empty. | ||
| 2055 | * \retval true if the stack is empty | ||
| 2056 | * \retval false otherwise | ||
| 2057 | */ | ||
| 2058 | 20928 | bool empty() const { | |
| 2059 | 20928 | return store_.empty(); | |
| 2060 | } | ||
| 2061 | |||
| 2062 | private: | ||
| 2063 | |||
| 2064 | /** | ||
| 2065 | * \brief The parameters and local | ||
| 2066 | * variables stored in a stack frame. | ||
| 2067 | */ | ||
| 2068 | struct Frame { | ||
| 2069 | // Parameters | ||
| 2070 | index_t t1; | ||
| 2071 | index_t new_t; | ||
| 2072 | Numeric::uint8 t1fbord ; | ||
| 2073 | |||
| 2074 | // Local variables | ||
| 2075 | Numeric::uint8 t1fprev ; | ||
| 2076 | Numeric::uint8 t1ft2 ; | ||
| 2077 | Numeric::uint8 t2ft1 ; | ||
| 2078 | }; | ||
| 2079 | |||
| 2080 | /** | ||
| 2081 | * \brief Gets the top of the stack. | ||
| 2082 | * \return a modifiable reference to the Frame on | ||
| 2083 | * the top of the stack | ||
| 2084 | * \pre !empty() | ||
| 2085 | */ | ||
| 2086 | 5982 | Frame& top() { | |
| 2087 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 5982 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
5982 | geo_debug_assert(!empty()); |
| 2088 | 5982 | return *store_.rbegin(); | |
| 2089 | } | ||
| 2090 | |||
| 2091 | /** | ||
| 2092 | * \brief Gets the top of the stack. | ||
| 2093 | * \return a const reference to the Frame on | ||
| 2094 | * the top of the stack | ||
| 2095 | * \pre !empty() | ||
| 2096 | */ | ||
| 2097 | 8964 | const Frame& top() const { | |
| 2098 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 8964 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
8964 | geo_debug_assert(!empty()); |
| 2099 | 8964 | return *store_.rbegin(); | |
| 2100 | } | ||
| 2101 | |||
| 2102 | std::vector<Frame> store_; | ||
| 2103 | }; | ||
| 2104 | |||
| 2105 | /** | ||
| 2106 | * \brief Creates a star of tetrahedra filling the conflict | ||
| 2107 | * zone. | ||
| 2108 | * \details For each tetrahedron facet on the border of the | ||
| 2109 | * conflict zone, a new tetrahedron is created, resting on | ||
| 2110 | * the facet and incident to vertex \p v. The function is | ||
| 2111 | * called recursively until the entire conflict zone is filled. | ||
| 2112 | * \param[in] v the index of the point to be inserted | ||
| 2113 | * \param[in] t1 index of a tetrahedron on the border | ||
| 2114 | * of the conflict zone. | ||
| 2115 | * \param[in] t1fbord index of the facet along which \p t_bndry | ||
| 2116 | * is incident to the border of the conflict zone | ||
| 2117 | * \param[in] t1fprev the facet of \p t_bndry connected to the | ||
| 2118 | * tetrahedron that \p t_bndry was reached from, or NO_INDEX | ||
| 2119 | * if it is the first tetrahedron. | ||
| 2120 | * \return the index of one the newly created tetrahedron | ||
| 2121 | */ | ||
| 2122 | 6 | index_t stellate_conflict_zone_iterative( | |
| 2123 | index_t v, index_t t1, index_t t1fbord, | ||
| 2124 | index_t t1fprev = NO_INDEX | ||
| 2125 | ) { | ||
| 2126 | // This function is de-recursified because some degenerate | ||
| 2127 | // inputs can cause stack overflow (system stack is limited to | ||
| 2128 | // a few megs). For instance, it can happen when a large number | ||
| 2129 | // of points are on the same sphere exactly. | ||
| 2130 | |||
| 2131 | // To de-recursify, it uses class StellateConflictStack | ||
| 2132 | // that emulates system's stack for storing functions's | ||
| 2133 | // parameters and local variables in all the nested stack | ||
| 2134 | // frames. | ||
| 2135 | |||
| 2136 |
1/2✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
|
6 | S2_.push(t1, t1fbord, t1fprev); |
| 2137 | |||
| 2138 | index_t new_t; // the newly created tetrahedron. | ||
| 2139 | |||
| 2140 | index_t t1ft2; // traverses the 4 facets of t1. | ||
| 2141 | |||
| 2142 | index_t t2; // the tetrahedron on the border of | ||
| 2143 | // the conflict zone that shares an | ||
| 2144 | // edge with t1 along t1ft2. | ||
| 2145 | |||
| 2146 | index_t t2fbord; // the facet of t2 on the border of | ||
| 2147 | // the conflict zone. | ||
| 2148 | |||
| 2149 | index_t t2ft1; // the facet of t2 that is incident to t1. | ||
| 2150 | |||
| 2151 | 1000 | entry_point: | |
| 2152 |
1/2✓ Branch 1 taken 1000 times.
✗ Branch 2 not taken.
|
1000 | S2_.get_parameters(t1, t1fbord, t1fprev); |
| 2153 | |||
| 2154 | |||
| 2155 |
2/8✓ Branch 1 taken 1000 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 1000 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
1000 | geo_debug_assert(owns_tet(t1)); |
| 2156 |
2/8✓ Branch 1 taken 1000 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 1000 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
1000 | geo_debug_assert(tet_adjacent(t1,t1fbord) != NO_INDEX); |
| 2157 |
3/10✓ Branch 1 taken 1000 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1000 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1000 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
|
1000 | geo_debug_assert(owns_tet(tet_adjacent(t1,t1fbord))); |
| 2158 |
2/8✓ Branch 1 taken 1000 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 1000 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
1000 | geo_debug_assert(tet_is_marked_as_conflict(t1)); |
| 2159 |
3/10✓ Branch 1 taken 1000 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1000 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1000 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
|
1000 | geo_debug_assert( |
| 2160 | !tet_is_marked_as_conflict(tet_adjacent(t1,t1fbord)) | ||
| 2161 | ); | ||
| 2162 | |||
| 2163 | // Create new tetrahedron with same vertices as t_bndry | ||
| 2164 |
5/10✓ Branch 1 taken 1000 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1000 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1000 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1000 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 1000 times.
✗ Branch 14 not taken.
|
1000 | new_t = new_tetrahedron( |
| 2165 | tet_vertex(t1,0), | ||
| 2166 | tet_vertex(t1,1), | ||
| 2167 | tet_vertex(t1,2), | ||
| 2168 | tet_vertex(t1,3) | ||
| 2169 | ); | ||
| 2170 | |||
| 2171 | // Replace in new_t the vertex opposite to t1fbord with v | ||
| 2172 |
1/2✓ Branch 1 taken 1000 times.
✗ Branch 2 not taken.
|
1000 | set_tet_vertex(new_t, t1fbord, v); |
| 2173 | |||
| 2174 | // Connect new_t with t1's neighbor accros t1fbord | ||
| 2175 | { | ||
| 2176 |
1/2✓ Branch 1 taken 1000 times.
✗ Branch 2 not taken.
|
1000 | index_t tbord = tet_adjacent(t1,t1fbord); |
| 2177 |
1/2✓ Branch 1 taken 1000 times.
✗ Branch 2 not taken.
|
1000 | set_tet_adjacent(new_t, t1fbord, tbord); |
| 2178 |
2/4✓ Branch 1 taken 1000 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1000 times.
✗ Branch 5 not taken.
|
1000 | set_tet_adjacent(tbord, find_tet_adjacent(tbord,t1), new_t); |
| 2179 | } | ||
| 2180 | |||
| 2181 | // Lookup new_t's neighbors accros its three other | ||
| 2182 | // facets and connect them | ||
| 2183 |
2/2✓ Branch 0 taken 4000 times.
✓ Branch 1 taken 1000 times.
|
5000 | for(t1ft2=0; t1ft2<4; ++t1ft2) { |
| 2184 | |||
| 2185 |
7/8✓ Branch 0 taken 3006 times.
✓ Branch 1 taken 994 times.
✓ Branch 3 taken 3006 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1506 times.
✓ Branch 6 taken 1500 times.
✓ Branch 7 taken 2500 times.
✓ Branch 8 taken 1500 times.
|
4000 | if(t1ft2 == t1fprev || tet_adjacent(new_t,t1ft2) != NO_INDEX) { |
| 2186 | 2500 | continue; | |
| 2187 | } | ||
| 2188 | |||
| 2189 | // Get t1's neighbor along the border of the conflict zone | ||
| 2190 |
3/4✓ Branch 1 taken 1500 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 994 times.
✓ Branch 4 taken 506 times.
|
1500 | if(!get_neighbor_along_conflict_zone_border( |
| 2191 | t1,t1fbord,t1ft2, t2,t2fbord,t2ft1 | ||
| 2192 | )) { | ||
| 2193 | // If t1's neighbor is not a new tetrahedron, | ||
| 2194 | // create a new tetrahedron through a recursive call. | ||
| 2195 |
1/2✓ Branch 1 taken 994 times.
✗ Branch 2 not taken.
|
994 | S2_.save_locals(new_t, t1ft2, t2ft1); |
| 2196 |
1/2✓ Branch 1 taken 994 times.
✗ Branch 2 not taken.
|
994 | S2_.push(t2, t2fbord, t2ft1); |
| 2197 | 994 | goto entry_point; | |
| 2198 | |||
| 2199 | 1000 | return_point: | |
| 2200 | // This is the return value of the called function. | ||
| 2201 | 1000 | index_t result = new_t; | |
| 2202 |
1/2✓ Branch 1 taken 1000 times.
✗ Branch 2 not taken.
|
1000 | S2_.pop(); |
| 2203 | |||
| 2204 | // Special case: we were in the outermost frame, | ||
| 2205 | // then we (truly) return from the function. | ||
| 2206 |
2/2✓ Branch 1 taken 6 times.
✓ Branch 2 taken 994 times.
|
1000 | if(S2_.empty()) { |
| 2207 | 6 | return result; | |
| 2208 | } | ||
| 2209 | |||
| 2210 |
1/2✓ Branch 1 taken 994 times.
✗ Branch 2 not taken.
|
994 | S2_.get_parameters(t1, t1fbord, t1fprev); |
| 2211 |
1/2✓ Branch 1 taken 994 times.
✗ Branch 2 not taken.
|
994 | S2_.get_locals(new_t, t1ft2, t2ft1); |
| 2212 | 994 | t2 = result; | |
| 2213 | } | ||
| 2214 | |||
| 2215 |
1/2✓ Branch 1 taken 1500 times.
✗ Branch 2 not taken.
|
1500 | set_tet_adjacent(t2, t2ft1, new_t); |
| 2216 |
1/2✓ Branch 1 taken 1500 times.
✗ Branch 2 not taken.
|
1500 | set_tet_adjacent(new_t, t1ft2, t2); |
| 2217 | } | ||
| 2218 | |||
| 2219 | // Except for the initial call (see "Special case" above), | ||
| 2220 | // the nested calls all come from the same location, | ||
| 2221 | // thus there is only one possible return point | ||
| 2222 | // (no need to push any return address). | ||
| 2223 | 1000 | goto return_point; | |
| 2224 | } | ||
| 2225 | |||
| 2226 | /** | ||
| 2227 | * \brief Finds the neighbor of a tetrahedron on the border of the | ||
| 2228 | * conflict zone. | ||
| 2229 | * \details This function is used by stellate_conflict_zone_iterative() | ||
| 2230 | * \param[in] t1 a tetrahedron on the border of the conflict zone | ||
| 2231 | * \param[in] t1fborder the local facet index of \p t1 along which it | ||
| 2232 | * is on the border of the conflict zone | ||
| 2233 | * \param[in] t1ft2 the local facet index of \p t1 that will be | ||
| 2234 | * traversed | ||
| 2235 | * \param[out] t2 a tetrahedron on the border of the conflict zone, | ||
| 2236 | * with an edge common to facets \p t1fborder and \p t1ft2 of | ||
| 2237 | * tetrahedron \p t1 | ||
| 2238 | * \param[out] t2fborder the local facet index of \p t2 along which it | ||
| 2239 | * is on the border of the conflict zone | ||
| 2240 | * \param[out] t2ft1 the local index of the facet of \p t2 that has a | ||
| 2241 | * common edge with facets \p t1fborder and \p t1ft2 of tetrahedron | ||
| 2242 | * \p t1 | ||
| 2243 | * \retval true if \p t2 is a newly created tetrahedron | ||
| 2244 | * \retval false if \p t2 is an old tetrahedron in conflict | ||
| 2245 | */ | ||
| 2246 | 1500 | bool get_neighbor_along_conflict_zone_border( | |
| 2247 | index_t t1, | ||
| 2248 | index_t t1fborder, | ||
| 2249 | index_t t1ft2, | ||
| 2250 | index_t& t2, | ||
| 2251 | index_t& t2fborder, | ||
| 2252 | index_t& t2ft1 | ||
| 2253 | ) const { | ||
| 2254 | |||
| 2255 | // Note: this function is a bit long for an inline function, | ||
| 2256 | // but I observed a (modest) performance gain doing so. | ||
| 2257 | |||
| 2258 | // Find two vertices that are both on facets new_f and f1 | ||
| 2259 | // (the edge around which we are turning) | ||
| 2260 | // This uses duality as follows: | ||
| 2261 | // Primal form (not used here): | ||
| 2262 | // halfedge_facet_[v1][v2] returns a facet that is incident | ||
| 2263 | // to both v1 and v2. | ||
| 2264 | // Dual form (used here): | ||
| 2265 | // halfedge_facet_[f1][f2] returns a vertex that both | ||
| 2266 | // f1 and f2 are incident to. | ||
| 2267 | index_t ev1 = | ||
| 2268 |
1/2✓ Branch 1 taken 1500 times.
✗ Branch 2 not taken.
|
1500 | tet_vertex(t1, index_t(halfedge_facet_[t1ft2][t1fborder])); |
| 2269 | index_t ev2 = | ||
| 2270 |
1/2✓ Branch 1 taken 1500 times.
✗ Branch 2 not taken.
|
1500 | tet_vertex(t1, index_t(halfedge_facet_[t1fborder][t1ft2])); |
| 2271 | |||
| 2272 | // Turn around edge [ev1,ev2] inside the conflict zone | ||
| 2273 | // until we reach again the boundary of the conflict zone. | ||
| 2274 | // Traversing inside the conflict zone is faster (as compared | ||
| 2275 | // to outside) since it traverses a smaller number of tets. | ||
| 2276 | 1500 | index_t cur_t = t1; | |
| 2277 | 1500 | index_t cur_f = t1ft2; | |
| 2278 |
1/2✓ Branch 1 taken 1500 times.
✗ Branch 2 not taken.
|
1500 | index_t next_t = tet_adjacent(cur_t,cur_f); |
| 2279 |
3/4✓ Branch 1 taken 3786 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2286 times.
✓ Branch 4 taken 1500 times.
|
3786 | while(tet_is_marked_as_conflict(next_t)) { |
| 2280 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 2286 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
2286 | geo_debug_assert(next_t != t1); |
| 2281 | 2286 | cur_t = next_t; | |
| 2282 |
1/2✓ Branch 1 taken 2286 times.
✗ Branch 2 not taken.
|
2286 | cur_f = get_facet_by_halfedge(cur_t,ev1,ev2); |
| 2283 |
1/2✓ Branch 1 taken 2286 times.
✗ Branch 2 not taken.
|
2286 | next_t = tet_adjacent(cur_t, cur_f); |
| 2284 | } | ||
| 2285 | |||
| 2286 | // At this point, cur_t is in conflict zone and | ||
| 2287 | // next_t is outside the conflict zone. | ||
| 2288 | index_t f12,f21; | ||
| 2289 |
1/2✓ Branch 1 taken 1500 times.
✗ Branch 2 not taken.
|
1500 | get_facets_by_halfedge(next_t, ev1, ev2, f12, f21); |
| 2290 |
1/2✓ Branch 1 taken 1500 times.
✗ Branch 2 not taken.
|
1500 | t2 = tet_adjacent(next_t,f21); |
| 2291 |
1/2✓ Branch 1 taken 1500 times.
✗ Branch 2 not taken.
|
1500 | index_t v_neigh_opposite = tet_vertex(next_t,f12); |
| 2292 |
1/2✓ Branch 1 taken 1500 times.
✗ Branch 2 not taken.
|
1500 | t2ft1 = find_tet_vertex(t2, v_neigh_opposite); |
| 2293 | 1500 | t2fborder = cur_f; | |
| 2294 | |||
| 2295 | // Test whether the found neighboring tet was created | ||
| 2296 | // (then return true) or is an old tet in conflict | ||
| 2297 | // (then return false). | ||
| 2298 | 1500 | return(t2 != cur_t); | |
| 2299 | } | ||
| 2300 | |||
| 2301 | /** | ||
| 2302 | * \brief Used by the (de-recursified) | ||
| 2303 | * stellate_conflict_zone_iterative() function. | ||
| 2304 | */ | ||
| 2305 | StellateConflictStack S2_; | ||
| 2306 | |||
| 2307 | /*************************** debugging ************************/ | ||
| 2308 | |||
| 2309 | /** | ||
| 2310 | * \brief For debugging purposes, displays a tetrahedron adjacency. | ||
| 2311 | * \param[in] t index of the tetrahedron to display. | ||
| 2312 | * \param[in] lf local index (0,1,2 or 3) of the tetrahedron | ||
| 2313 | * facet adjacenty to display. | ||
| 2314 | */ | ||
| 2315 | 19112 | void show_tet_adjacent(index_t t, index_t lf) const { | |
| 2316 | 19112 | index_t adj = tet_adjacent(t, lf); | |
| 2317 |
2/2✓ Branch 0 taken 9488 times.
✓ Branch 1 taken 9624 times.
|
19112 | if(adj != NO_INDEX) { |
| 2318 |
2/2✓ Branch 1 taken 228 times.
✓ Branch 2 taken 9260 times.
|
9488 | std::cerr << (tet_is_in_list(adj) ? '*' : ' '); |
| 2319 | } | ||
| 2320 | 19112 | std::cerr << adj; | |
| 2321 | 19112 | std::cerr << ' '; | |
| 2322 | 19112 | } | |
| 2323 | |||
| 2324 | |||
| 2325 | /** | ||
| 2326 | * \brief For debugging purposes, displays a tetrahedron. | ||
| 2327 | * \param[in] t index of the tetrahedron to display. | ||
| 2328 | */ | ||
| 2329 | 4778 | void show_tet(index_t t) const { | |
| 2330 | std::cerr << "tet" | ||
| 2331 |
2/2✓ Branch 2 taken 2488 times.
✓ Branch 3 taken 2290 times.
|
4778 | << (tet_is_in_list(t) ? '*' : ' ') |
| 2332 | 4778 | << t | |
| 2333 | 4778 | << ", v=[" | |
| 2334 | 4778 | << tet_vertex(t, 0) | |
| 2335 | 4778 | << ' ' | |
| 2336 | 4778 | << tet_vertex(t, 1) | |
| 2337 | 4778 | << ' ' | |
| 2338 | 4778 | << tet_vertex(t, 2) | |
| 2339 | 4778 | << ' ' | |
| 2340 | 4778 | << tet_vertex(t, 3) | |
| 2341 | 4778 | << "] adj=["; | |
| 2342 | 4778 | show_tet_adjacent(t, 0); | |
| 2343 | 4778 | show_tet_adjacent(t, 1); | |
| 2344 | 4778 | show_tet_adjacent(t, 2); | |
| 2345 | 4778 | show_tet_adjacent(t, 3); | |
| 2346 | 4778 | std::cerr << "] "; | |
| 2347 | |||
| 2348 |
2/2✓ Branch 0 taken 19112 times.
✓ Branch 1 taken 4778 times.
|
23890 | for(index_t f = 0; f < 4; ++f) { |
| 2349 | 19112 | std::cerr << 'f' << f << ':'; | |
| 2350 |
2/2✓ Branch 0 taken 57336 times.
✓ Branch 1 taken 19112 times.
|
76448 | for(index_t v = 0; v < 3; ++v) { |
| 2351 | 57336 | std::cerr << tet_vertex(t, tet_facet_vertex(f,v)) | |
| 2352 | 57336 | << ','; | |
| 2353 | } | ||
| 2354 | 19112 | std::cerr << ' '; | |
| 2355 | } | ||
| 2356 | 4778 | std::cerr << std::endl; | |
| 2357 | 4778 | } | |
| 2358 | |||
| 2359 | public: | ||
| 2360 | |||
| 2361 | /** | ||
| 2362 | * \brief For debugging purposes, tests some combinatorial properties. | ||
| 2363 | */ | ||
| 2364 | 4 | void check_combinatorics(bool verbose) const { | |
| 2365 |
1/2✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
|
4 | if(verbose) { |
| 2366 |
1/2✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
|
4 | std::cerr << std::endl; |
| 2367 | } | ||
| 2368 | 4 | bool ok = true; | |
| 2369 |
1/2✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
|
4 | std::vector<bool> v_has_tet(nb_vertices(), false); |
| 2370 |
2/2✓ Branch 1 taken 4778 times.
✓ Branch 2 taken 4 times.
|
4782 | for(index_t t = 0; t < max_t(); ++t) { |
| 2371 |
3/4✓ Branch 1 taken 4778 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2488 times.
✓ Branch 4 taken 2290 times.
|
4778 | if(tet_is_free(t)) { |
| 2372 |
1/2✓ Branch 0 taken 2488 times.
✗ Branch 1 not taken.
|
2488 | if(verbose) { |
| 2373 |
1/2✓ Branch 1 taken 2488 times.
✗ Branch 2 not taken.
|
2488 | std::cerr << "-Deleted tet: "; |
| 2374 |
1/2✓ Branch 1 taken 2488 times.
✗ Branch 2 not taken.
|
2488 | show_tet(t); |
| 2375 | } | ||
| 2376 | } else { | ||
| 2377 |
1/2✓ Branch 0 taken 2290 times.
✗ Branch 1 not taken.
|
2290 | if(verbose) { |
| 2378 |
1/2✓ Branch 1 taken 2290 times.
✗ Branch 2 not taken.
|
2290 | std::cerr << "Checking tet: "; |
| 2379 |
1/2✓ Branch 1 taken 2290 times.
✗ Branch 2 not taken.
|
2290 | show_tet(t); |
| 2380 | } | ||
| 2381 |
2/2✓ Branch 0 taken 9160 times.
✓ Branch 1 taken 2290 times.
|
11450 | for(index_t lf = 0; lf < 4; ++lf) { |
| 2382 |
2/4✓ Branch 1 taken 9160 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 9160 times.
|
9160 | if(tet_adjacent(t, lf) == NO_INDEX) { |
| 2383 | ✗ | std::cerr << lf << ":Missing adjacent tet" | |
| 2384 | ✗ | << std::endl; | |
| 2385 | ✗ | ok = false; | |
| 2386 |
2/4✓ Branch 1 taken 9160 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 9160 times.
|
9160 | } else if(tet_adjacent(t, lf) == t) { |
| 2387 | ✗ | std::cerr << lf << ":Tet is adjacent to itself" | |
| 2388 | ✗ | << std::endl; | |
| 2389 | ✗ | ok = false; | |
| 2390 | } else { | ||
| 2391 |
1/2✓ Branch 1 taken 9160 times.
✗ Branch 2 not taken.
|
9160 | index_t t2 = tet_adjacent(t, lf); |
| 2392 | 9160 | bool found = false; | |
| 2393 |
2/2✓ Branch 0 taken 36640 times.
✓ Branch 1 taken 9160 times.
|
45800 | for(index_t lf2 = 0; lf2 < 4; ++lf2) { |
| 2394 |
3/4✓ Branch 1 taken 36640 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 9160 times.
✓ Branch 4 taken 27480 times.
|
36640 | if(tet_adjacent(t2, lf2) == t) { |
| 2395 | 9160 | found = true; | |
| 2396 | } | ||
| 2397 | } | ||
| 2398 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 9160 times.
|
9160 | if(!found) { |
| 2399 | std::cerr | ||
| 2400 | ✗ | << lf | |
| 2401 | ✗ | << ":Adjacent link is not bidirectional" | |
| 2402 | ✗ | << std::endl; | |
| 2403 | ✗ | ok = false; | |
| 2404 | } | ||
| 2405 | } | ||
| 2406 | } | ||
| 2407 | 2290 | index_t nb_infinite = 0; | |
| 2408 |
2/2✓ Branch 0 taken 9160 times.
✓ Branch 1 taken 2290 times.
|
11450 | for(index_t lv = 0; lv < 4; ++lv) { |
| 2409 |
3/4✓ Branch 1 taken 9160 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 732 times.
✓ Branch 4 taken 8428 times.
|
9160 | if(tet_vertex(t, lv) == NO_INDEX) { |
| 2410 | 732 | ++nb_infinite; | |
| 2411 | } | ||
| 2412 | } | ||
| 2413 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2290 times.
|
2290 | if(nb_infinite > 1) { |
| 2414 | ✗ | ok = false; | |
| 2415 | ✗ | std::cerr << "More than one infinite vertex" | |
| 2416 | ✗ | << std::endl; | |
| 2417 | } | ||
| 2418 | } | ||
| 2419 |
2/2✓ Branch 0 taken 19112 times.
✓ Branch 1 taken 4778 times.
|
23890 | for(index_t lv = 0; lv < 4; ++lv) { |
| 2420 |
1/2✓ Branch 1 taken 19112 times.
✗ Branch 2 not taken.
|
19112 | index_t v = tet_vertex(t, lv); |
| 2421 |
4/4✓ Branch 0 taken 8756 times.
✓ Branch 1 taken 10356 times.
✓ Branch 2 taken 8428 times.
✓ Branch 3 taken 328 times.
|
19112 | if(v != NO_INDEX && v != VERTEX_OF_DELETED_TET) { |
| 2422 | 8428 | v_has_tet[v] = true; | |
| 2423 | } | ||
| 2424 | } | ||
| 2425 | } | ||
| 2426 |
2/2✓ Branch 1 taken 458 times.
✓ Branch 2 taken 4 times.
|
462 | for(index_t v = 0; v < nb_vertices(); ++v) { |
| 2427 |
1/2✗ Branch 2 not taken.
✓ Branch 3 taken 458 times.
|
458 | if(!v_has_tet[v]) { |
| 2428 | ✗ | if(verbose) { | |
| 2429 | ✗ | std::cerr << "Vertex " << v | |
| 2430 | ✗ | << " is isolated (duplicated ?)" << std::endl; | |
| 2431 | } | ||
| 2432 | } | ||
| 2433 | } | ||
| 2434 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
4 | geo_assert(ok); |
| 2435 |
1/2✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
|
4 | if(verbose) { |
| 2436 |
1/2✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
|
4 | std::cerr << std::endl; |
| 2437 | } | ||
| 2438 |
3/6✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
|
4 | std::cerr << std::endl << "Delaunay Combi OK" << std::endl; |
| 2439 | 4 | } | |
| 2440 | |||
| 2441 | |||
| 2442 | /** | ||
| 2443 | * \brief For debugging purposes, test some geometrical properties. | ||
| 2444 | */ | ||
| 2445 | 4 | void check_geometry(bool verbose) const { | |
| 2446 | 4 | bool ok = true; | |
| 2447 |
2/2✓ Branch 1 taken 4778 times.
✓ Branch 2 taken 4 times.
|
4782 | for(index_t t = 0; t < max_t(); ++t) { |
| 2448 |
2/2✓ Branch 1 taken 2290 times.
✓ Branch 2 taken 2488 times.
|
4778 | if(!tet_is_free(t)) { |
| 2449 | 2290 | index_t v0 = tet_vertex(t, 0); | |
| 2450 | 2290 | index_t v1 = tet_vertex(t, 1); | |
| 2451 | 2290 | index_t v2 = tet_vertex(t, 2); | |
| 2452 | 2290 | index_t v3 = tet_vertex(t, 3); | |
| 2453 |
2/2✓ Branch 1 taken 498422 times.
✓ Branch 2 taken 2290 times.
|
500712 | for(index_t v = 0; v < nb_vertices(); ++v) { |
| 2454 |
8/8✓ Branch 0 taken 496134 times.
✓ Branch 1 taken 2288 times.
✓ Branch 2 taken 494140 times.
✓ Branch 3 taken 1994 times.
✓ Branch 4 taken 492014 times.
✓ Branch 5 taken 2126 times.
✓ Branch 6 taken 2020 times.
✓ Branch 7 taken 489994 times.
|
498422 | if(v == v0 || v == v1 || v == v2 || v == v3) { |
| 2455 | 8428 | continue; | |
| 2456 | } | ||
| 2457 |
1/2✗ Branch 2 not taken.
✓ Branch 3 taken 489994 times.
|
489994 | if(tet_is_in_conflict(t, vertex_ptr(v))) { |
| 2458 | ✗ | ok = false; | |
| 2459 | ✗ | if(verbose) { | |
| 2460 | ✗ | std::cerr << "Tet " << t << | |
| 2461 | ✗ | " is in conflict with vertex " << v | |
| 2462 | ✗ | << std::endl; | |
| 2463 | |||
| 2464 | ✗ | std::cerr << " offending tet: "; | |
| 2465 | ✗ | show_tet(t); | |
| 2466 | } | ||
| 2467 | } | ||
| 2468 | } | ||
| 2469 | } | ||
| 2470 | } | ||
| 2471 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
4 | geo_assert(ok); |
| 2472 | 4 | std::cerr << std::endl << "Delaunay Geo OK" << std::endl; | |
| 2473 | 4 | } | |
| 2474 | |||
| 2475 | private: | ||
| 2476 | ParallelDelaunay3d* master_; | ||
| 2477 | index_t nb_vertices_; | ||
| 2478 | const double* vertices_; | ||
| 2479 | const double* heights_; | ||
| 2480 | index_t* reorder_; | ||
| 2481 | index_t dimension_; | ||
| 2482 | index_t vertex_stride_; | ||
| 2483 | bool weighted_; | ||
| 2484 | index_t max_t_; | ||
| 2485 | index_t max_used_t_; | ||
| 2486 | |||
| 2487 | vector<index_t>& cell_to_v_store_; | ||
| 2488 | vector<index_t>& cell_to_cell_store_; | ||
| 2489 | vector<index_t>& cell_next_; | ||
| 2490 | CellStatusArray& cell_status_; | ||
| 2491 | |||
| 2492 | index_t first_free_; | ||
| 2493 | index_t nb_free_; | ||
| 2494 | bool memory_overflow_; | ||
| 2495 | |||
| 2496 | index_t v1_,v2_,v3_,v4_; // The first four vertices | ||
| 2497 | |||
| 2498 | vector<index_t> S_; | ||
| 2499 | index_t nb_tets_to_create_; | ||
| 2500 | index_t t_boundary_; // index of a tet,facet on the bndry | ||
| 2501 | index_t f_boundary_; // of the conflict zone. | ||
| 2502 | |||
| 2503 | bool direction_; | ||
| 2504 | index_t work_begin_; | ||
| 2505 | index_t work_end_; | ||
| 2506 | index_t b_hint_; | ||
| 2507 | index_t e_hint_; | ||
| 2508 | bool finished_; | ||
| 2509 | |||
| 2510 | // Whenever acquire_tet() is unsuccessful, contains | ||
| 2511 | // the index of the thread that was interfering | ||
| 2512 | // (shifted to the left by 1 !!) | ||
| 2513 | CellStatusArray::thread_index_t interfering_thread_; | ||
| 2514 | |||
| 2515 | #ifdef GEO_DEBUG | ||
| 2516 | index_t nb_acquired_tets_; | ||
| 2517 | #endif | ||
| 2518 | |||
| 2519 | vector<index_t> tets_to_delete_; | ||
| 2520 | vector<index_t> tets_to_release_; | ||
| 2521 | |||
| 2522 | index_t nb_rollbacks_; | ||
| 2523 | index_t nb_failed_locate_; | ||
| 2524 | |||
| 2525 | std::condition_variable cond_; | ||
| 2526 | std::mutex mutex_; | ||
| 2527 | |||
| 2528 | /** | ||
| 2529 | * \brief Gives the indexing of tetrahedron facet | ||
| 2530 | * vertices. | ||
| 2531 | * \details tet_facet_vertex[lf][lv] gives the | ||
| 2532 | * local vertex index (in 0,1,2,3) from a | ||
| 2533 | * local facet index lf (in 0,1,2,3) and a | ||
| 2534 | * local vertex index within the facet (in 0,1,2). | ||
| 2535 | */ | ||
| 2536 | static char tet_facet_vertex_[4][3]; | ||
| 2537 | |||
| 2538 | /** | ||
| 2539 | * \brief Gives a local facet index by | ||
| 2540 | * halfedge extremities local indices. | ||
| 2541 | */ | ||
| 2542 | static char halfedge_facet_[4][4]; | ||
| 2543 | |||
| 2544 | /** | ||
| 2545 | * \brief Stores the triangles on the boundary | ||
| 2546 | * of the cavity, for faster generation of the | ||
| 2547 | * new tetrahedra. | ||
| 2548 | */ | ||
| 2549 | Cavity cavity_; | ||
| 2550 | }; | ||
| 2551 | |||
| 2552 | |||
| 2553 | char Delaunay3dThread::halfedge_facet_[4][4] = { | ||
| 2554 | {4, 2, 3, 1}, | ||
| 2555 | {3, 4, 0, 2}, | ||
| 2556 | {1, 3, 4, 0}, | ||
| 2557 | {2, 0, 1, 4} | ||
| 2558 | }; | ||
| 2559 | |||
| 2560 | // tet facet vertex is such that the tetrahedron | ||
| 2561 | // formed with: | ||
| 2562 | // vertex lv | ||
| 2563 | // tet_facet_vertex[lv][0] | ||
| 2564 | // tet_facet_vertex[lv][1] | ||
| 2565 | // tet_facet_vertex[lv][2] | ||
| 2566 | // has the same orientation as the original tetrahedron for | ||
| 2567 | // any vertex lv. | ||
| 2568 | |||
| 2569 | char Delaunay3dThread::tet_facet_vertex_[4][3] = { | ||
| 2570 | {1, 2, 3}, | ||
| 2571 | {0, 3, 2}, | ||
| 2572 | {3, 0, 1}, | ||
| 2573 | {1, 0, 2} | ||
| 2574 | }; | ||
| 2575 | |||
| 2576 | |||
| 2577 | /*************************************************************************/ | ||
| 2578 | |||
| 2579 | 5 | ParallelDelaunay3d::ParallelDelaunay3d( | |
| 2580 | coord_index_t dimension | ||
| 2581 | 5 | ) : Delaunay(dimension) { | |
| 2582 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
5 | if(dimension != 3 && dimension != 4) { |
| 2583 | ✗ | throw InvalidDimension(dimension, "Delaunay3d", "3 or 4"); | |
| 2584 | } | ||
| 2585 | |||
| 2586 |
1/2✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
|
5 | geo_cite_with_info( |
| 2587 | "DBLP:journals/cj/Bowyer81", | ||
| 2588 | "One of the two initial references to the algorithm, " | ||
| 2589 | "discovered independently and simultaneously by Bowyer and Watson." | ||
| 2590 | ); | ||
| 2591 |
1/2✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
|
5 | geo_cite_with_info( |
| 2592 | "journals/cj/Watson81", | ||
| 2593 | "One of the two initial references to the algorithm, " | ||
| 2594 | "discovered independently and simultaneously by Bowyer and Watson." | ||
| 2595 | ); | ||
| 2596 |
1/2✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
|
5 | geo_cite_with_info( |
| 2597 | "DBLP:conf/compgeom/AmentaCR03", | ||
| 2598 | "Using spatial sorting has a dramatic impact on the performances." | ||
| 2599 | ); | ||
| 2600 |
1/2✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
|
5 | geo_cite_with_info( |
| 2601 | "DBLP:journals/comgeo/FunkeMN05", | ||
| 2602 | "Initializing \\verb|locate()| with a non-exact version " | ||
| 2603 | " (structural filtering) gains (a bit of) performance." | ||
| 2604 | ); | ||
| 2605 |
1/2✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
|
5 | geo_cite_with_info( |
| 2606 | "DBLP:journals/comgeo/BoissonnatDPTY02", | ||
| 2607 | "The idea of traversing the cavity from inside " | ||
| 2608 | " used in GEOGRAM is inspired by the implementation of " | ||
| 2609 | " \\verb|Delaunay_triangulation_3| in CGAL." | ||
| 2610 | ); | ||
| 2611 |
1/2✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
|
5 | geo_cite_with_info( |
| 2612 | "DBLP:conf/imr/Si06", | ||
| 2613 | "The triangulation data structure used in GEOGRAM is inspired " | ||
| 2614 | "by Tetgen." | ||
| 2615 | ); | ||
| 2616 |
1/2✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
|
5 | geo_cite_with_info( |
| 2617 | "DBLP:journals/ijfcs/DevillersPT02", | ||
| 2618 | "Analysis of the different versions of the line walk algorithm " | ||
| 2619 | " used by \\verb|locate()|." | ||
| 2620 | ); | ||
| 2621 | |||
| 2622 | 5 | weighted_ = (dimension == 4); | |
| 2623 | // In weighted mode, vertices are 4d but combinatorics is 3d. | ||
| 2624 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
|
5 | if(weighted_) { |
| 2625 | ✗ | cell_size_ = 4; | |
| 2626 | ✗ | cell_v_stride_ = 4; | |
| 2627 | ✗ | cell_neigh_stride_ = 4; | |
| 2628 | } | ||
| 2629 |
2/4✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5 times.
✗ Branch 5 not taken.
|
10 | debug_mode_ = CmdLine::get_arg_bool("dbg:delaunay"); |
| 2630 |
2/4✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5 times.
✗ Branch 5 not taken.
|
5 | verbose_debug_mode_ = CmdLine::get_arg_bool("dbg:delaunay_verbose"); |
| 2631 |
3/4✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 4 times.
✓ Branch 3 taken 1 times.
|
5 | debug_mode_ = (debug_mode_ || verbose_debug_mode_); |
| 2632 |
2/4✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5 times.
✗ Branch 5 not taken.
|
5 | benchmark_mode_ = CmdLine::get_arg_bool("dbg:delaunay_benchmark"); |
| 2633 | 5 | } | |
| 2634 | |||
| 2635 | 5 | void ParallelDelaunay3d::set_vertices( | |
| 2636 | index_t nb_vertices, const double* vertices | ||
| 2637 | ) { | ||
| 2638 | 5 | Stopwatch* W = nullptr ; | |
| 2639 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1 times.
|
5 | if(benchmark_mode_) { |
| 2640 |
3/8✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 4 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
|
12 | W = new Stopwatch("DelInternal"); |
| 2641 | } | ||
| 2642 | |||
| 2643 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
|
5 | if(weighted_) { |
| 2644 | ✗ | heights_.resize(nb_vertices); | |
| 2645 | ✗ | for(index_t i = 0; i < nb_vertices; ++i) { | |
| 2646 | // Client code uses 4d embedding with ti = sqrt(W - wi) | ||
| 2647 | // where W = max(wi) | ||
| 2648 | // We recompute the standard "shifted" lifting on | ||
| 2649 | // the paraboloid from it. | ||
| 2650 | // (we use wi - W, everything is shifted by W, but | ||
| 2651 | // we do not care since the power diagram is invariant | ||
| 2652 | // by a translation of all weights). | ||
| 2653 | ✗ | double w = -geo_sqr(vertices[4 * i + 3]); | |
| 2654 | ✗ | heights_[i] = -w + | |
| 2655 | ✗ | geo_sqr(vertices[4 * i]) + | |
| 2656 | ✗ | geo_sqr(vertices[4 * i + 1]) + | |
| 2657 | ✗ | geo_sqr(vertices[4 * i + 2]); | |
| 2658 | } | ||
| 2659 | } | ||
| 2660 | 5 | Delaunay::set_vertices(nb_vertices, vertices); | |
| 2661 | |||
| 2662 | 5 | index_t expected_tetra = nb_vertices * 7; | |
| 2663 | |||
| 2664 | // Allocate the tetrahedra | ||
| 2665 | 5 | cell_to_v_store_.assign(expected_tetra * 4,NO_INDEX); | |
| 2666 | 5 | cell_to_cell_store_.assign(expected_tetra * 4,NO_INDEX); | |
| 2667 | 5 | cell_next_.assign(expected_tetra,NO_INDEX); | |
| 2668 | 5 | cell_status_.resize(expected_tetra); | |
| 2669 | |||
| 2670 | // Reorder the points | ||
| 2671 |
1/2✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
|
5 | if(do_reorder_) { |
| 2672 | 5 | compute_BRIO_order( | |
| 2673 | 5 | nb_vertices, vertex_ptr(0), reorder_, | |
| 2674 | 5 | 3, dimension(), | |
| 2675 | 64, 0.125, | ||
| 2676 | &levels_ | ||
| 2677 | ); | ||
| 2678 | } else { | ||
| 2679 | ✗ | reorder_.resize(nb_vertices); | |
| 2680 | ✗ | for(index_t i = 0; i < nb_vertices; ++i) { | |
| 2681 | ✗ | reorder_[i] = i; | |
| 2682 | } | ||
| 2683 | ✗ | geo_debug_assert(levels_[0] == 0); | |
| 2684 | ✗ | geo_debug_assert(levels_[levels_.size()-1] == nb_vertices); | |
| 2685 | } | ||
| 2686 | |||
| 2687 | 5 | double sorting_time = 0; | |
| 2688 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1 times.
|
5 | if(benchmark_mode_) { |
| 2689 | 4 | sorting_time = W->elapsed_time(); | |
| 2690 |
3/6✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
|
12 | Logger::out("DelInternal1") << "BRIO sorting:" |
| 2691 |
1/2✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
|
4 | << sorting_time |
| 2692 |
1/2✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
|
4 | << std::endl; |
| 2693 | } | ||
| 2694 | |||
| 2695 | // Create the threads | ||
| 2696 | // The maximum number of threads is limited by the number of bits used by | ||
| 2697 | // cell_status_ (see delaunay_sync.h) | ||
| 2698 | 5 | index_t nb_threads = std::min( | |
| 2699 |
1/2✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
|
5 | Process::maximum_concurrent_threads(), |
| 2700 | CellStatusArray::MAX_THREADS | ||
| 2701 | 5 | ); | |
| 2702 | 5 | index_t pool_size = expected_tetra / nb_threads; | |
| 2703 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
|
5 | if (pool_size == 0) { |
| 2704 | // There are more threads than expected_tetra | ||
| 2705 | ✗ | pool_size = 1; | |
| 2706 | ✗ | nb_threads = expected_tetra; | |
| 2707 | } | ||
| 2708 | 5 | index_t pool_begin = 0; | |
| 2709 | 5 | threads_.clear(); | |
| 2710 |
2/2✓ Branch 0 taken 20 times.
✓ Branch 1 taken 5 times.
|
25 | for(index_t t=0; t<nb_threads; ++t) { |
| 2711 | 20 | index_t pool_end = | |
| 2712 |
2/2✓ Branch 0 taken 15 times.
✓ Branch 1 taken 5 times.
|
20 | (t == nb_threads - 1) ? expected_tetra : pool_begin + pool_size; |
| 2713 |
1/2✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
|
20 | threads_.push_back( |
| 2714 |
3/8✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 20 times.
✗ Branch 5 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 20 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
|
40 | new Delaunay3dThread(this, pool_begin, pool_end) |
| 2715 | ); | ||
| 2716 | 20 | pool_begin = pool_end; | |
| 2717 | } | ||
| 2718 | |||
| 2719 | |||
| 2720 | // Create first tetrahedron and triangulate first set of points | ||
| 2721 | // in sequential mode. | ||
| 2722 | |||
| 2723 | |||
| 2724 | 5 | index_t lvl = 1; | |
| 2725 |
5/6✓ Branch 1 taken 4 times.
✓ Branch 2 taken 5 times.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 4 times.
✓ Branch 7 taken 5 times.
|
9 | while(lvl < (levels_.size() - 1) && levels_[lvl] < 1000) { |
| 2726 | 4 | ++lvl; | |
| 2727 | } | ||
| 2728 | |||
| 2729 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1 times.
|
5 | if(benchmark_mode_) { |
| 2730 |
2/4✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
|
12 | Logger::out("PDEL") |
| 2731 |
4/8✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 4 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 4 times.
✗ Branch 12 not taken.
|
4 | << "Using " << levels_.size()-1 << " levels" << std::endl; |
| 2732 |
2/4✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
|
12 | Logger::out("PDEL") |
| 2733 |
2/4✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
|
4 | << "Levels 0 - " << lvl-1 |
| 2734 |
1/2✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
|
4 | << ": bootstraping with first levels in sequential mode" |
| 2735 |
1/2✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
|
4 | << std::endl; |
| 2736 | } | ||
| 2737 | Delaunay3dThread* thread0 = | ||
| 2738 | 5 | static_cast<Delaunay3dThread*>(threads_[0].get()); | |
| 2739 | 5 | thread0->create_first_tetrahedron(); | |
| 2740 | 5 | thread0->set_work(levels_[0], levels_[lvl]); | |
| 2741 | 5 | thread0->run(); | |
| 2742 | |||
| 2743 | 5 | index_t first_lvl = lvl; | |
| 2744 | |||
| 2745 | // Insert points in all BRIO levels | ||
| 2746 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 5 times.
|
5 | for(; lvl<levels_.size()-1; ++lvl) { |
| 2747 | |||
| 2748 | ✗ | if(benchmark_mode_) { | |
| 2749 | ✗ | Logger::out("PDEL") << "Level " | |
| 2750 | ✗ | << lvl << " : start" << std::endl; | |
| 2751 | } | ||
| 2752 | |||
| 2753 | ✗ | index_t lvl_b = levels_[lvl]; | |
| 2754 | ✗ | index_t lvl_e = levels_[lvl+1]; | |
| 2755 | ✗ | index_t work_size = (lvl_e - lvl_b)/index_t(threads_.size()); | |
| 2756 | |||
| 2757 | // Initialize threads | ||
| 2758 | ✗ | index_t b = lvl_b; | |
| 2759 | ✗ | for(index_t t=0; t<threads_.size(); ++t) { | |
| 2760 | ✗ | index_t e = t == threads_.size()-1 ? lvl_e : b+work_size; | |
| 2761 | Delaunay3dThread* thread = | ||
| 2762 | ✗ | static_cast<Delaunay3dThread*>(threads_[t].get()); | |
| 2763 | |||
| 2764 | // Copy the indices of the first created tetrahedron | ||
| 2765 | // and the maximum valid tetrahedron index max_t_ | ||
| 2766 | ✗ | if(lvl == first_lvl && t!=0) { | |
| 2767 | ✗ | thread->initialize_from(thread0); | |
| 2768 | } | ||
| 2769 | ✗ | thread->set_work(b,e); | |
| 2770 | ✗ | b = e; | |
| 2771 | } | ||
| 2772 | ✗ | Process::run_threads(threads_); | |
| 2773 | } | ||
| 2774 | |||
| 2775 | |||
| 2776 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1 times.
|
5 | if(benchmark_mode_) { |
| 2777 | 4 | index_t tot_rollbacks = 0 ; | |
| 2778 | 4 | index_t tot_failed_locate = 0 ; | |
| 2779 |
2/2✓ Branch 1 taken 16 times.
✓ Branch 2 taken 4 times.
|
20 | for(index_t t=0; t<threads_.size(); ++t) { |
| 2780 | Delaunay3dThread* thread = | ||
| 2781 | 16 | static_cast<Delaunay3dThread*>(threads_[t].get()); | |
| 2782 |
2/4✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 16 times.
✗ Branch 5 not taken.
|
32 | Logger::out("PDEL") |
| 2783 |
2/4✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
✓ Branch 6 taken 16 times.
✗ Branch 7 not taken.
|
16 | << "thread " << std::setw(3) << t << " : " |
| 2784 |
1/2✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
|
16 | << std::setw(3) |
| 2785 |
1/2✓ Branch 2 taken 16 times.
✗ Branch 3 not taken.
|
16 | << thread->nb_rollbacks() << " rollbacks " |
| 2786 |
1/2✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
|
16 | << std::setw(3) |
| 2787 |
2/4✓ Branch 2 taken 16 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 16 times.
✗ Branch 6 not taken.
|
16 | << thread->nb_failed_locate() << " restarted locate" |
| 2788 |
1/2✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
|
16 | << std::endl; |
| 2789 | 16 | tot_rollbacks += thread->nb_rollbacks(); | |
| 2790 | 16 | tot_failed_locate += thread->nb_failed_locate(); | |
| 2791 | } | ||
| 2792 |
4/8✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 4 times.
✗ Branch 11 not taken.
|
8 | Logger::out("PDEL") << "------------------" << std::endl; |
| 2793 |
3/6✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
|
12 | Logger::out("PDEL") << "total: " |
| 2794 |
2/4✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
|
4 | << tot_rollbacks << " rollbacks " |
| 2795 |
2/4✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
|
4 | << tot_failed_locate << " restarted locate" |
| 2796 |
1/2✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
|
4 | << std::endl; |
| 2797 | } | ||
| 2798 | |||
| 2799 | // Run threads sequentialy, to insert missing points if | ||
| 2800 | // memory overflow was encountered (in sequential mode, | ||
| 2801 | // dynamic memory growing works) | ||
| 2802 | |||
| 2803 | 5 | index_t nb_sequential_points = 0; | |
| 2804 |
2/2✓ Branch 1 taken 20 times.
✓ Branch 2 taken 5 times.
|
25 | for(index_t t=0; t<threads_.size(); ++t) { |
| 2805 | Delaunay3dThread* t1 = | ||
| 2806 | 20 | static_cast<Delaunay3dThread*>(threads_[t].get()); | |
| 2807 | |||
| 2808 | 20 | nb_sequential_points += t1->work_size(); | |
| 2809 | |||
| 2810 |
2/2✓ Branch 0 taken 15 times.
✓ Branch 1 taken 5 times.
|
20 | if(t != 0) { |
| 2811 | // We need to copy max_t_ from previous thread, | ||
| 2812 | // since the memory pool may have grown. | ||
| 2813 | Delaunay3dThread* t2 = | ||
| 2814 | 15 | static_cast<Delaunay3dThread*>(threads_[t-1].get()); | |
| 2815 | 15 | t1->initialize_from(t2); | |
| 2816 | } | ||
| 2817 | 20 | t1->run(); | |
| 2818 | } | ||
| 2819 | |||
| 2820 | // If some tetrahedra were created in sequential mode, then | ||
| 2821 | // the maximum valid tetrahedron index was increased by all | ||
| 2822 | // the threads in increasing number, so we copy it from the | ||
| 2823 | // last thread into thread0 since we use thread0 afterwards | ||
| 2824 | // to do the "compaction" afterwards. | ||
| 2825 | |||
| 2826 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
|
5 | if(nb_sequential_points != 0) { |
| 2827 | Delaunay3dThread* t0 = | ||
| 2828 | ✗ | static_cast<Delaunay3dThread*>(threads_[0].get()); | |
| 2829 | Delaunay3dThread* tn = | ||
| 2830 | static_cast<Delaunay3dThread*>( | ||
| 2831 | ✗ | threads_[threads_.size()-1].get() | |
| 2832 | ); | ||
| 2833 | ✗ | t0->initialize_from(tn); | |
| 2834 | } | ||
| 2835 | |||
| 2836 | |||
| 2837 | |||
| 2838 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1 times.
|
5 | if(benchmark_mode_) { |
| 2839 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
|
4 | if(nb_sequential_points != 0) { |
| 2840 | ✗ | Logger::out("PDEL") << "Local thread memory overflow occurred:" | |
| 2841 | ✗ | << std::endl; | |
| 2842 | ✗ | Logger::out("PDEL") << nb_sequential_points | |
| 2843 | ✗ | << " points inserted in sequential mode" | |
| 2844 | ✗ | << std::endl; | |
| 2845 | } else { | ||
| 2846 |
2/4✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
|
12 | Logger::out("PDEL") |
| 2847 |
1/2✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
|
4 | << "All the points were inserted in parallel mode" |
| 2848 |
1/2✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
|
4 | << std::endl; |
| 2849 | } | ||
| 2850 | } | ||
| 2851 | |||
| 2852 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1 times.
|
5 | if(benchmark_mode_) { |
| 2853 |
3/6✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
|
12 | Logger::out("DelInternal2") << "Core insertion algo:" |
| 2854 |
2/4✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
|
4 | << W->elapsed_time() - sorting_time |
| 2855 |
1/2✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
|
4 | << std::endl; |
| 2856 | } | ||
| 2857 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1 times.
|
5 | delete W; |
| 2858 | |||
| 2859 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1 times.
|
5 | if(debug_mode_) { |
| 2860 |
2/2✓ Branch 1 taken 16 times.
✓ Branch 2 taken 4 times.
|
20 | for(index_t i=0; i<threads_.size(); ++i) { |
| 2861 | 16 | std::cerr << i << " : " << | |
| 2862 | 16 | static_cast<Delaunay3dThread*>(threads_[i].get()) | |
| 2863 | 16 | ->max_t() << std::endl; | |
| 2864 | } | ||
| 2865 | |||
| 2866 | 4 | thread0->check_combinatorics(verbose_debug_mode_); | |
| 2867 | 4 | thread0->check_geometry(verbose_debug_mode_); | |
| 2868 | } | ||
| 2869 | |||
| 2870 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1 times.
|
5 | if(benchmark_mode_) { |
| 2871 |
3/8✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 4 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
|
12 | W = new Stopwatch("DelCompress"); |
| 2872 | } | ||
| 2873 | |||
| 2874 | // Compress cell_to_v_store_ and cell_to_cell_store_ | ||
| 2875 | // (remove free and virtual tetrahedra). | ||
| 2876 | // Since cell_next_ is not used at this point, | ||
| 2877 | // we reuse it for storing the conversion array that | ||
| 2878 | // maps old tet indices to new tet indices | ||
| 2879 | // Note: tet_is_real() uses the previous value of | ||
| 2880 | // cell_next(), but we are processing indices | ||
| 2881 | // in increasing order and since old2new[t] is always | ||
| 2882 | // smaller or equal to t, we never overwrite a value | ||
| 2883 | // before needing it. | ||
| 2884 | |||
| 2885 | 5 | vector<index_t>& old2new = cell_next_; | |
| 2886 | 5 | index_t nb_tets = 0; | |
| 2887 | 5 | index_t nb_tets_to_delete = 0; | |
| 2888 | |||
| 2889 | { | ||
| 2890 |
2/2✓ Branch 1 taken 46707 times.
✓ Branch 2 taken 5 times.
|
46712 | for(index_t t = 0; t < thread0->max_t(); ++t) { |
| 2891 | 46707 | if( | |
| 2892 |
6/6✓ Branch 0 taken 44318 times.
✓ Branch 1 taken 2389 times.
✓ Branch 3 taken 20169 times.
✓ Branch 4 taken 24149 times.
✓ Branch 5 taken 24928 times.
✓ Branch 6 taken 21779 times.
|
69265 | (keep_infinite_ && !thread0->tet_is_free(t)) || |
| 2893 |
2/2✓ Branch 1 taken 779 times.
✓ Branch 2 taken 21779 times.
|
22558 | thread0->tet_is_real(t) |
| 2894 | ) { | ||
| 2895 |
2/2✓ Branch 0 taken 18587 times.
✓ Branch 1 taken 6341 times.
|
24928 | if(t != nb_tets) { |
| 2896 | 18587 | Memory::copy( | |
| 2897 | 18587 | &cell_to_v_store_[nb_tets * 4], | |
| 2898 | 18587 | &cell_to_v_store_[t * 4], | |
| 2899 | 4 * sizeof(index_t) | ||
| 2900 | ); | ||
| 2901 | 18587 | Memory::copy( | |
| 2902 | 18587 | &cell_to_cell_store_[nb_tets * 4], | |
| 2903 | 18587 | &cell_to_cell_store_[t * 4], | |
| 2904 | 4 * sizeof(index_t) | ||
| 2905 | ); | ||
| 2906 | } | ||
| 2907 | 24928 | old2new[t] = nb_tets; | |
| 2908 | 24928 | ++nb_tets; | |
| 2909 | } else { | ||
| 2910 | 21779 | old2new[t] = NO_INDEX; | |
| 2911 | 21779 | ++nb_tets_to_delete; | |
| 2912 | } | ||
| 2913 | } | ||
| 2914 | |||
| 2915 | 5 | cell_to_v_store_.resize(4 * nb_tets); | |
| 2916 | 5 | cell_to_cell_store_.resize(4 * nb_tets); | |
| 2917 |
2/2✓ Branch 0 taken 99712 times.
✓ Branch 1 taken 5 times.
|
99717 | for(index_t i = 0; i < 4 * nb_tets; ++i) { |
| 2918 | 99712 | index_t t = cell_to_cell_store_[i]; | |
| 2919 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 99712 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
99712 | geo_debug_assert(t != NO_INDEX); |
| 2920 | 99712 | t = old2new[t]; | |
| 2921 | // Note: t can be equal to -1 when a real tet is | ||
| 2922 | // adjacent to a virtual one (and this is how the | ||
| 2923 | // rest of Vorpaline expects to see tets on the | ||
| 2924 | // border). | ||
| 2925 |
3/8✓ Branch 0 taken 96596 times.
✓ Branch 1 taken 3116 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 96596 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
99712 | geo_debug_assert(!(keep_infinite_ && t == NO_INDEX)); |
| 2926 | 99712 | cell_to_cell_store_[i] = t; | |
| 2927 | } | ||
| 2928 | } | ||
| 2929 | |||
| 2930 | // In "keep_infinite" mode, we reorder the cells in such | ||
| 2931 | // a way that finite cells have indices [0..nb_finite_cells_-1] | ||
| 2932 | // and infinite cells have indices [nb_finite_cells_ .. nb_cells_-1] | ||
| 2933 | |||
| 2934 |
2/2✓ Branch 0 taken 3 times.
✓ Branch 1 taken 2 times.
|
5 | if(keep_infinite_) { |
| 2935 | 3 | nb_finite_cells_ = 0; | |
| 2936 | 3 | index_t finite_ptr = 0; | |
| 2937 | 3 | index_t infinite_ptr = nb_tets - 1; | |
| 2938 | for(;;) { | ||
| 2939 |
2/2✓ Branch 1 taken 20703 times.
✓ Branch 2 taken 1599 times.
|
22302 | while(thread0->tet_is_finite(finite_ptr)) { |
| 2940 | 20703 | old2new[finite_ptr] = finite_ptr; | |
| 2941 | 20703 | ++finite_ptr; | |
| 2942 | 20703 | ++nb_finite_cells_; | |
| 2943 | } | ||
| 2944 |
2/2✓ Branch 1 taken 254 times.
✓ Branch 2 taken 1599 times.
|
1853 | while(!thread0->tet_is_finite(infinite_ptr)) { |
| 2945 | 254 | old2new[infinite_ptr] = infinite_ptr; | |
| 2946 | 254 | --infinite_ptr; | |
| 2947 | } | ||
| 2948 |
2/2✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1596 times.
|
1599 | if(finite_ptr > infinite_ptr) { |
| 2949 | 3 | break; | |
| 2950 | } | ||
| 2951 | 1596 | old2new[finite_ptr] = infinite_ptr; | |
| 2952 | 1596 | old2new[infinite_ptr] = finite_ptr; | |
| 2953 | 1596 | ++nb_finite_cells_; | |
| 2954 |
2/2✓ Branch 0 taken 6384 times.
✓ Branch 1 taken 1596 times.
|
7980 | for(index_t lf=0; lf<4; ++lf) { |
| 2955 | 6384 | std::swap( | |
| 2956 | 6384 | cell_to_cell_store_[4*finite_ptr + lf], | |
| 2957 | 6384 | cell_to_cell_store_[4*infinite_ptr + lf] | |
| 2958 | ); | ||
| 2959 | } | ||
| 2960 |
2/2✓ Branch 0 taken 6384 times.
✓ Branch 1 taken 1596 times.
|
7980 | for(index_t lv=0; lv<4; ++lv) { |
| 2961 | 6384 | std::swap( | |
| 2962 | 6384 | cell_to_v_store_[4*finite_ptr + lv], | |
| 2963 | 6384 | cell_to_v_store_[4*infinite_ptr + lv] | |
| 2964 | ); | ||
| 2965 | } | ||
| 2966 | 1596 | ++finite_ptr; | |
| 2967 | 1596 | --infinite_ptr; | |
| 2968 | 1596 | } | |
| 2969 |
2/2✓ Branch 0 taken 96596 times.
✓ Branch 1 taken 3 times.
|
96599 | for(index_t i = 0; i < 4 * nb_tets; ++i) { |
| 2970 | 96596 | index_t t = cell_to_cell_store_[i]; | |
| 2971 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 96596 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
96596 | geo_debug_assert(t != NO_INDEX); |
| 2972 | 96596 | t = old2new[t]; | |
| 2973 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 96596 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
96596 | geo_debug_assert(t != NO_INDEX); |
| 2974 | 96596 | cell_to_cell_store_[i] = t; | |
| 2975 | } | ||
| 2976 | } | ||
| 2977 | |||
| 2978 | |||
| 2979 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1 times.
|
5 | if(benchmark_mode_) { |
| 2980 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
|
4 | if(keep_infinite_) { |
| 2981 |
2/4✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
|
6 | Logger::out("DelCompress") |
| 2982 |
2/4✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
|
2 | << "Removed " << nb_tets_to_delete |
| 2983 |
2/4✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
|
2 | << " tets (free list)" << std::endl; |
| 2984 | } else { | ||
| 2985 |
2/4✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
|
6 | Logger::out("DelCompress") |
| 2986 |
2/4✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
|
2 | << "Removed " << nb_tets_to_delete |
| 2987 |
2/4✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
|
2 | << " tets (free list and infinite)" << std::endl; |
| 2988 | } | ||
| 2989 | } | ||
| 2990 | |||
| 2991 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1 times.
|
5 | delete W; |
| 2992 | |||
| 2993 | 5 | set_arrays( | |
| 2994 | nb_tets, | ||
| 2995 | 5 | cell_to_v_store_.data(), | |
| 2996 | 5 | cell_to_cell_store_.data() | |
| 2997 | ); | ||
| 2998 | 5 | } | |
| 2999 | |||
| 3000 | ✗ | index_t ParallelDelaunay3d::nearest_vertex(const double* p) const { | |
| 3001 | // TODO | ||
| 3002 | ✗ | return Delaunay::nearest_vertex(p); | |
| 3003 | } | ||
| 3004 | |||
| 3005 | ✗ | void ParallelDelaunay3d::set_BRIO_levels(const vector<index_t>& levels) { | |
| 3006 | ✗ | levels_ = levels; | |
| 3007 | ✗ | } | |
| 3008 | |||
| 3009 | } | ||
| 3010 | |||
| 3011 | #endif | ||
| 3012 |