src/meshTools/triSurface/triangleFuncs/triangleFuncs.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "triangleFuncs.H"
  27 #include "pointField.H"
  28 #include "treeBoundBox.H"
  29 #include "SortableList.H"
  30 #include "boolList.H"
  31
  32 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  33
  34 void Foam::triangleFuncs::setIntersection
  35 (
  36     const point& oppositeSidePt,
  37     const scalar oppositeSign,
  38
  39     const point& thisSidePt,
  40     const scalar thisSign,
  41
  42     const scalar tol,
  43
  44     point& pt
  45 )
  46 {
  47     scalar denom = oppositeSign - thisSign;
  48
  49     if (mag(denom) < tol)
  50     {
  51         // If almost does not cut choose one which certainly cuts.
  52         pt = oppositeSidePt;
  53     }
  54     else
  55     {
  56         pt = oppositeSidePt + oppositeSign/denom*(thisSidePt - oppositeSidePt);
  57     }
  58 }
  59
  60
  61 void Foam::triangleFuncs::selectPt
  62 (
  63     const bool select0,
  64     const point& p0,
  65     const point& p1,
  66     point& min
  67 )
  68 {
  69     if (select0)
  70     {
  71         min = p0;
  72     }
  73     else
  74     {
  75         min = p1;
  76     }
  77 }
  78
  79
  80 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  81
  82 // Intersect triangle with parallel edges aligned with axis i0.
  83 // Returns true (and intersection in pInter) if any of them intersects triangle.
  84 bool Foam::triangleFuncs::intersectAxesBundle
  85 (
  86     const point& V0,
  87     const point& V10,
  88     const point& V20,
  89     const label i0,
  90     const pointField& origin,
  91     const scalar maxLength,
  92     point& pInter
  93 )
  94 {
  95     // Based on Graphics Gems - Fast Ray Triangle intersection.
  96     // Since direction is coordinate axis there is no need to do projection,
  97     // we can directly check u,v components for inclusion in triangle.
  98
  99     // Get other components
 100     label i1 = (i0 + 1) % 3;
 101     label i2 = (i1 + 1) % 3;
 102
 103     scalar u1 = V10[i1];
 104     scalar v1 = V10[i2];
 105
 106     scalar u2 = V20[i1];
 107     scalar v2 = V20[i2];
 108
 109     scalar localScale = mag(u1)+mag(v1)+mag(u2)+mag(v2);
 110
 111     scalar det = v2*u1 - u2*v1;
 112
 113     // Fix for  V0:(-31.71428 0 -15.10714)
 114     //          V10:(-1.285715 8.99165e-16 -1.142858)
 115     //          V20:(0 0 -1.678573)
 116     //          i0:0
 117     if (localScale < VSMALL || Foam::mag(det)/localScale < SMALL)
 118     {
 119         // Triangle parallel to dir
 120         return false;
 121     }
 122
 123     forAll(origin, originI)
 124     {
 125         const point& P = origin[originI];
 126
 127         scalar u0 = P[i1] - V0[i1];
 128         scalar v0 = P[i2] - V0[i2];
 129
 130         scalar alpha = 0;
 131         scalar beta = 0;
 132         bool inter = false;
 133
 134         if (Foam::mag(u1) < ROOTVSMALL)
 135         {
 136             beta = u0/u2;
 137             if ((beta >= 0) && (beta <= 1))
 138             {
 139                 alpha = (v0 - beta*v2)/v1;
 140                 inter = ((alpha >= 0) && ((alpha + beta) <= 1));
 141             }
 142         }
 143         else
 144         {
 145             beta = (v0*u1 - u0*v1)/det;
 146             if ((beta >= 0) && (beta <= 1))
 147             {
 148                 alpha = (u0 - beta*u2)/u1;
 149                 inter = ((alpha >= 0) && ((alpha + beta) <= 1));
 150             }
 151         }
 152
 153         if (inter)
 154         {
 155             pInter = V0 + alpha*V10 + beta*V20;
 156             scalar s = (pInter - origin[originI])[i0];
 157
 158             if ((s >= 0) && (s <= maxLength))
 159             {
 160                 return true;
 161             }
 162         }
 163     }
 164     return false;
 165 }
 166
 167
 168 // Intersect triangle with bounding box. Return true if
 169 // any of the faces of bb intersect triangle.
 170 // Note: so returns false if triangle inside bb.
 171 bool Foam::triangleFuncs::intersectBb
 172 (
 173     const point& p0,
 174     const point& p1,
 175     const point& p2,
 176     const treeBoundBox& cubeBb
 177 )
 178 {
 179     const vector p10 = p1 - p0;
 180     const vector p20 = p2 - p0;
 181
 182     // cubeBb points; counted as if cell with vertex0 at cubeBb.min().
 183     const point& min = cubeBb.min();
 184     const point& max = cubeBb.max();
 185
 186     const point& cube0 = min;
 187     const point  cube1(min.x(), min.y(), max.z());
 188     const point  cube2(max.x(), min.y(), max.z());
 189     const point  cube3(max.x(), min.y(), min.z());
 190
 191     const point  cube4(min.x(), max.y(), min.z());
 192     const point  cube5(min.x(), max.y(), max.z());
 193     const point  cube7(max.x(), max.y(), min.z());
 194
 195     //
 196     // Intersect all 12 edges of cube with triangle
 197     //
 198
 199     point pInter;
 200     pointField origin(4);
 201     // edges in x direction
 202     origin[0] = cube0;
 203     origin[1] = cube1;
 204     origin[2] = cube5;
 205     origin[3] = cube4;
 206
 207     scalar maxSx = max.x() - min.x();
 208
 209     if (intersectAxesBundle(p0, p10, p20, 0, origin, maxSx, pInter))
 210     {
 211         return true;
 212     }
 213
 214     // edges in y direction
 215     origin[0] = cube0;
 216     origin[1] = cube1;
 217     origin[2] = cube2;
 218     origin[3] = cube3;
 219
 220     scalar maxSy = max.y() - min.y();
 221
 222     if (intersectAxesBundle(p0, p10, p20, 1, origin, maxSy, pInter))
 223     {
 224         return true;
 225     }
 226
 227     // edges in z direction
 228     origin[0] = cube0;
 229     origin[1] = cube3;
 230     origin[2] = cube7;
 231     origin[3] = cube4;
 232
 233     scalar maxSz = max.z() - min.z();
 234
 235     if (intersectAxesBundle(p0, p10, p20, 2, origin, maxSz, pInter))
 236     {
 237         return true;
 238     }
 239
 240
 241     // Intersect triangle edges with bounding box
 242     if (cubeBb.intersects(p0, p1, pInter))
 243     {
 244         return true;
 245     }
 246     if (cubeBb.intersects(p1, p2, pInter))
 247     {
 248         return true;
 249     }
 250     if (cubeBb.intersects(p2, p0, pInter))
 251     {
 252         return true;
 253     }
 254
 255     return false;
 256 }
 257
 258
 259 //// Intersect triangle with bounding box. Return true if
 260 //// any of the faces of bb intersect triangle.
 261 //// Note: so returns false if triangle inside bb.
 262 //bool Foam::triangleFuncs::intersectBbExact
 263 //(
 264 //    const point& p0,
 265 //    const point& p1,
 266 //    const point& p2,
 267 //    const treeBoundBox& cubeBb
 268 //)
 269 //{
 270 //    const point& min = cubeBb.min();
 271 //    const point& max = cubeBb.max();
 272 //
 273 //    const point& cube0 = min;
 274 //    const point  cube1(min.x(), min.y(), max.z());
 275 //    const point  cube2(max.x(), min.y(), max.z());
 276 //    const point  cube3(max.x(), min.y(), min.z());
 277 //
 278 //    const point  cube4(min.x(), max.y(), min.z());
 279 //    const point  cube5(min.x(), max.y(), max.z());
 280 //    const point& cube6 = max;
 281 //    const point  cube7(max.x(), max.y(), min.z());
 282 //
 283 //    // Test intersection of triangle with twelve edges of box.
 284 //    {
 285 //        triPointRef tri(p0, p1, p2);
 286 //        if (tri.intersectionExact(cube0, cube1).hit())
 287 //        {
 288 //            return true;
 289 //        }
 290 //        if (tri.intersectionExact(cube1, cube2).hit())
 291 //        {
 292 //            return true;
 293 //        }
 294 //        if (tri.intersectionExact(cube2, cube3).hit())
 295 //        {
 296 //            return true;
 297 //        }
 298 //        if (tri.intersectionExact(cube3, cube0).hit())
 299 //        {
 300 //            return true;
 301 //        }
 302 //
 303 //        if (tri.intersectionExact(cube4, cube5).hit())
 304 //        {
 305 //            return true;
 306 //        }
 307 //        if (tri.intersectionExact(cube5, cube6).hit())
 308 //        {
 309 //            return true;
 310 //        }
 311 //        if (tri.intersectionExact(cube6, cube7).hit())
 312 //        {
 313 //            return true;
 314 //        }
 315 //        if (tri.intersectionExact(cube7, cube4).hit())
 316 //        {
 317 //            return true;
 318 //        }
 319 //
 320 //        if (tri.intersectionExact(cube0, cube4).hit())
 321 //        {
 322 //            return true;
 323 //        }
 324 //        if (tri.intersectionExact(cube1, cube5).hit())
 325 //        {
 326 //            return true;
 327 //        }
 328 //        if (tri.intersectionExact(cube2, cube6).hit())
 329 //        {
 330 //            return true;
 331 //        }
 332 //        if (tri.intersectionExact(cube3, cube7).hit())
 333 //        {
 334 //            return true;
 335 //        }
 336 //    }
 337 //    // Test intersection of triangle edges with bounding box
 338 //    {
 339 //        triPointRef tri(cube0, cube1, cube2);
 340 //        if (tri.intersectionExact(p0, p1).hit())
 341 //        {
 342 //            return true;
 343 //        }
 344 //        if (tri.intersectionExact(p1, p2).hit())
 345 //        {
 346 //            return true;
 347 //        }
 348 //        if (tri.intersectionExact(p2, p0).hit())
 349 //        {
 350 //            return true;
 351 //        }
 352 //    }
 353 //    {
 354 //        triPointRef tri(cube2, cube3, cube0);
 355 //        if (tri.intersectionExact(p0, p1).hit())
 356 //        {
 357 //            return true;
 358 //        }
 359 //        if (tri.intersectionExact(p1, p2).hit())
 360 //        {
 361 //            return true;
 362 //        }
 363 //        if (tri.intersectionExact(p2, p0).hit())
 364 //        {
 365 //            return true;
 366 //        }
 367 //    }
 368 //    {
 369 //        triPointRef tri(cube4, cube5, cube6);
 370 //        if (tri.intersectionExact(p0, p1).hit())
 371 //        {
 372 //            return true;
 373 //        }
 374 //        if (tri.intersectionExact(p1, p2).hit())
 375 //        {
 376 //            return true;
 377 //        }
 378 //        if (tri.intersectionExact(p2, p0).hit())
 379 //        {
 380 //            return true;
 381 //        }
 382 //    }
 383 //    {
 384 //        triPointRef tri(cube6, cube7, cube4);
 385 //        if (tri.intersectionExact(p0, p1).hit())
 386 //        {
 387 //            return true;
 388 //        }
 389 //        if (tri.intersectionExact(p1, p2).hit())
 390 //        {
 391 //            return true;
 392 //        }
 393 //        if (tri.intersectionExact(p2, p0).hit())
 394 //        {
 395 //            return true;
 396 //        }
 397 //    }
 398 //
 399 //
 400 //    {
 401 //        triPointRef tri(cube4, cube5, cube1);
 402 //        if (tri.intersectionExact(p0, p1).hit())
 403 //        {
 404 //            return true;
 405 //        }
 406 //        if (tri.intersectionExact(p1, p2).hit())
 407 //        {
 408 //            return true;
 409 //        }
 410 //        if (tri.intersectionExact(p2, p0).hit())
 411 //        {
 412 //            return true;
 413 //        }
 414 //    }
 415 //    {
 416 //        triPointRef tri(cube1, cube0, cube4);
 417 //        if (tri.intersectionExact(p0, p1).hit())
 418 //        {
 419 //            return true;
 420 //        }
 421 //        if (tri.intersectionExact(p1, p2).hit())
 422 //        {
 423 //            return true;
 424 //        }
 425 //        if (tri.intersectionExact(p2, p0).hit())
 426 //        {
 427 //            return true;
 428 //        }
 429 //    }
 430 //    {
 431 //        triPointRef tri(cube7, cube6, cube2);
 432 //        if (tri.intersectionExact(p0, p1).hit())
 433 //        {
 434 //            return true;
 435 //        }
 436 //        if (tri.intersectionExact(p1, p2).hit())
 437 //        {
 438 //            return true;
 439 //        }
 440 //        if (tri.intersectionExact(p2, p0).hit())
 441 //        {
 442 //            return true;
 443 //        }
 444 //    }
 445 //    {
 446 //        triPointRef tri(cube2, cube3, cube7);
 447 //        if (tri.intersectionExact(p0, p1).hit())
 448 //        {
 449 //            return true;
 450 //        }
 451 //        if (tri.intersectionExact(p1, p2).hit())
 452 //        {
 453 //            return true;
 454 //        }
 455 //        if (tri.intersectionExact(p2, p0).hit())
 456 //        {
 457 //            return true;
 458 //        }
 459 //    }
 460 //
 461 //    {
 462 //        triPointRef tri(cube0, cube4, cube7);
 463 //        if (tri.intersectionExact(p0, p1).hit())
 464 //        {
 465 //            return true;
 466 //        }
 467 //        if (tri.intersectionExact(p1, p2).hit())
 468 //        {
 469 //            return true;
 470 //        }
 471 //        if (tri.intersectionExact(p2, p0).hit())
 472 //        {
 473 //            return true;
 474 //        }
 475 //    }
 476 //    {
 477 //        triPointRef tri(cube7, cube3, cube0);
 478 //        if (tri.intersectionExact(p0, p1).hit())
 479 //        {
 480 //            return true;
 481 //        }
 482 //        if (tri.intersectionExact(p1, p2).hit())
 483 //        {
 484 //            return true;
 485 //        }
 486 //        if (tri.intersectionExact(p2, p0).hit())
 487 //        {
 488 //            return true;
 489 //        }
 490 //    }
 491 //    {
 492 //        triPointRef tri(cube1, cube5, cube6);
 493 //        if (tri.intersectionExact(p0, p1).hit())
 494 //        {
 495 //            return true;
 496 //        }
 497 //        if (tri.intersectionExact(p1, p2).hit())
 498 //        {
 499 //            return true;
 500 //        }
 501 //        if (tri.intersectionExact(p2, p0).hit())
 502 //        {
 503 //            return true;
 504 //        }
 505 //    }
 506 //    {
 507 //        triPointRef tri(cube6, cube2, cube1);
 508 //        if (tri.intersectionExact(p0, p1).hit())
 509 //        {
 510 //            return true;
 511 //        }
 512 //        if (tri.intersectionExact(p1, p2).hit())
 513 //        {
 514 //            return true;
 515 //        }
 516 //        if (tri.intersectionExact(p2, p0).hit())
 517 //        {
 518 //            return true;
 519 //        }
 520 //    }
 521 //    return false;
 522 //}
 523
 524
 525 bool Foam::triangleFuncs::intersect
 526 (
 527     const point& va0,
 528     const point& va10,
 529     const point& va20,
 530
 531     const point& base,
 532     const point& normal,
 533
 534     point& pInter0,
 535     point& pInter1
 536 )
 537 {
 538     // Get triangle normal
 539     vector na = va10 ^ va20;
 540     scalar magArea = mag(na);
 541     na/magArea;
 542
 543     if (mag(na & normal) > (1 - SMALL))
 544     {
 545         // Parallel
 546         return false;
 547     }
 548
 549     const point va1 = va0 + va10;
 550     const point va2 = va0 + va20;
 551
 552     // Find the triangle point on the other side.
 553     scalar sign0 = (va0 - base) & normal;
 554     scalar sign1 = (va1 - base) & normal;
 555     scalar sign2 = (va2 - base) & normal;
 556
 557     label oppositeVertex = -1;
 558
 559     if (sign0 < 0)
 560     {
 561         if (sign1 < 0)
 562         {
 563             if (sign2 < 0)
 564             {
 565                 // All on same side of plane
 566                 return false;
 567             }
 568             else    // sign2 >= 0
 569             {
 570                 // 2 on opposite side.
 571                 oppositeVertex = 2;
 572             }
 573         }
 574         else    // sign1 >= 0
 575         {
 576             if (sign2 < 0)
 577             {
 578                 // 1 on opposite side.
 579                 oppositeVertex = 1;
 580             }
 581             else
 582             {
 583                 // 0 on opposite side.
 584                 oppositeVertex = 0;
 585             }
 586         }
 587     }
 588     else    // sign0 >= 0
 589     {
 590         if (sign1 < 0)
 591         {
 592             if (sign2 < 0)
 593             {
 594                 // 0 on opposite side.
 595                 oppositeVertex = 0;
 596             }
 597             else    // sign2 >= 0
 598             {
 599                 // 1 on opposite side.
 600                 oppositeVertex = 1;
 601             }
 602         }
 603         else    // sign1 >= 0
 604         {
 605             if (sign2 < 0)
 606             {
 607                 // 2 on opposite side.
 608                 oppositeVertex = 2;
 609             }
 610             else    // sign2 >= 0
 611             {
 612                 // All on same side of plane
 613                 return false;
 614             }
 615         }
 616     }
 617
 618     scalar tol = SMALL*Foam::sqrt(magArea);
 619
 620     if (oppositeVertex == 0)
 621     {
 622         // 0 on opposite side. Cut edges 01 and 02
 623         setIntersection(va0, sign0, va1, sign1, tol, pInter0);
 624         setIntersection(va0, sign0, va2, sign2, tol, pInter1);
 625     }
 626     else if (oppositeVertex == 1)
 627     {
 628         // 1 on opposite side. Cut edges 10 and 12
 629         setIntersection(va1, sign1, va0, sign0, tol, pInter0);
 630         setIntersection(va1, sign1, va2, sign2, tol, pInter1);
 631     }
 632     else // oppositeVertex == 2
 633     {
 634         // 2 on opposite side. Cut edges 20 and 21
 635         setIntersection(va2, sign2, va0, sign0, tol, pInter0);
 636         setIntersection(va2, sign2, va1, sign1, tol, pInter1);
 637     }
 638
 639     return true;
 640 }
 641
 642
 643 bool Foam::triangleFuncs::intersect
 644 (
 645     const point& va0,
 646     const point& va10,
 647     const point& va20,
 648
 649     const point& vb0,
 650     const point& vb10,
 651     const point& vb20,
 652
 653     point& pInter0,
 654     point& pInter1
 655 )
 656 {
 657     // Get triangle normals
 658     vector na = va10 ^ va20;
 659     na/mag(na);
 660
 661     vector nb = vb10 ^ vb20;
 662     nb/mag(nb);
 663
 664     // Calculate intersection of triangle a with plane of b
 665     point planeB0;
 666     point planeB1;
 667     if (!intersect(va0, va10, va20, vb0, nb, planeB0, planeB1))
 668     {
 669         return false;
 670     }
 671
 672     //       ,,  triangle b with plane of a
 673     point planeA0;
 674     point planeA1;
 675     if (!intersect(vb0, vb10, vb20, va0, na, planeA0, planeA1))
 676     {
 677         return false;
 678     }
 679
 680     // Now check if intersections overlap (w.r.t. intersection of the two
 681     // planes)
 682
 683     vector intersection(na ^ nb);
 684
 685     scalar coordB0 = planeB0 & intersection;
 686     scalar coordB1 = planeB1 & intersection;
 687
 688     scalar coordA0 = planeA0 & intersection;
 689     scalar coordA1 = planeA1 & intersection;
 690
 691     // Put information in indexable form for sorting.
 692     List<point*> pts(4);
 693     boolList isFromB(4);
 694     SortableList<scalar> sortCoords(4);
 695
 696     pts[0] = &planeB0;
 697     isFromB[0] = true;
 698     sortCoords[0] = coordB0;
 699
 700     pts[1] = &planeB1;
 701     isFromB[1] = true;
 702     sortCoords[1] = coordB1;
 703
 704     pts[2] = &planeA0;
 705     isFromB[2] = false;
 706     sortCoords[2] = coordA0;
 707
 708     pts[3] = &planeA1;
 709     isFromB[3] = false;
 710     sortCoords[3] = coordA1;
 711
 712     sortCoords.sort();
 713
 714     const labelList& indices = sortCoords.indices();
 715
 716     if (isFromB[indices[0]] == isFromB[indices[1]])
 717     {
 718         // Entry 0 and 1 are of same region (both a or both b). Hence that
 719         // region does not overlap.
 720         return false;
 721     }
 722     else
 723     {
 724         // Different type. Start of overlap at indices[1], end at indices[2]
 725         pInter0 = *pts[indices[1]];
 726         pInter1 = *pts[indices[2]];
 727
 728         return true;
 729     }
 730 }
 731
 732
 733 // ************************************************************************* //