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