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