src/triSurface/faceTriangulation/faceTriangulation.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 \*---------------------------------------------------------------------------*/
  26
  27 #include "faceTriangulation.H"
  28 #include "plane.H"
  29
  30
  31 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  32
  33 const Foam::scalar Foam::faceTriangulation::edgeRelTol = 1E-6;
  34
  35
  36 //- Edge to the right of face vertex i
  37 Foam::label Foam::faceTriangulation::right(const label, label i)
  38 {
  39     return i;
  40 }
  41
  42
  43 //- Edge to the left of face vertex i
  44 Foam::label Foam::faceTriangulation::left(const label size, label i)
  45 {
  46     return i == 0 ? size - 1 : (i - 1);
  47 }
  48
  49
  50 // Calculate (normalized) edge vectors.
  51 // edges[i] gives edge between point i+1 and i.
  52 Foam::tmp<Foam::vectorField> Foam::faceTriangulation::calcEdges
  53 (
  54     const face& f,
  55     const pointField& points
  56 )
  57 {
  58     tmp<vectorField> tedges(new vectorField(f.size()));
  59     vectorField& edges = tedges();
  60
  61     forAll(f, i)
  62     {
  63         label ni = f.fcIndex(i);
  64
  65         point thisPt = points[f[i]];
  66         point nextPt = points[f[ni]];
  67
  68         vector vec(nextPt - thisPt);
  69         vec /= mag(vec) + VSMALL;
  70
  71         edges[i] = vec;
  72     }
  73
  74     return tedges;
  75 }
  76
  77
  78 // Calculates half angle components of angle from e0 to e1
  79 void Foam::faceTriangulation::calcHalfAngle
  80 (
  81     const vector& normal,
  82     const vector& e0,
  83     const vector& e1,
  84     scalar& cosHalfAngle,
  85     scalar& sinHalfAngle
  86 )
  87 {
  88     // truncate cos to +-1 to prevent negative numbers
  89     scalar cos = max(-1, min(1, e0 & e1));
  90
  91     scalar sin = (e0 ^ e1) & normal;
  92
  93     if (sin < -SMALL)
  94     {
  95         // 3rd or 4th quadrant
  96         cosHalfAngle = - Foam::sqrt(0.5*(1 + cos));
  97         sinHalfAngle = Foam::sqrt(0.5*(1 - cos));
  98     }
  99     else
 100     {
 101         // 1st or 2nd quadrant
 102         cosHalfAngle = Foam::sqrt(0.5*(1 + cos));
 103         sinHalfAngle = Foam::sqrt(0.5*(1 - cos));
 104     }
 105 }
 106
 107
 108 // Calculate intersection point between edge p1-p2 and ray (in 2D).
 109 // Return true and intersection point if intersection between p1 and p2.
 110 Foam::pointHit Foam::faceTriangulation::rayEdgeIntersect
 111 (
 112     const vector& normal,
 113     const point& rayOrigin,
 114     const vector& rayDir,
 115     const point& p1,
 116     const point& p2,
 117     scalar& posOnEdge
 118 )
 119 {
 120     // Start off from miss
 121     pointHit result(p1);
 122
 123     // Construct plane normal to rayDir and intersect
 124     const vector y = normal ^ rayDir;
 125
 126     posOnEdge = plane(rayOrigin, y).normalIntersect(p1, (p2-p1));
 127
 128     // Check intersection to left of p1 or right of p2
 129     if ((posOnEdge < 0) || (posOnEdge > 1))
 130     {
 131         // Miss
 132     }
 133     else
 134     {
 135         // Check intersection behind rayOrigin
 136         point intersectPt = p1 + posOnEdge * (p2 - p1);
 137
 138         if (((intersectPt - rayOrigin) & rayDir) < 0)
 139         {
 140             // Miss
 141         }
 142         else
 143         {
 144             // Hit
 145             result.setHit();
 146             result.setPoint(intersectPt);
 147             result.setDistance(mag(intersectPt - rayOrigin));
 148         }
 149     }
 150     return result;
 151 }
 152
 153
 154 // Return true if triangle given its three points (anticlockwise ordered)
 155 // contains point
 156 bool Foam::faceTriangulation::triangleContainsPoint
 157 (
 158     const vector& n,
 159     const point& p0,
 160     const point& p1,
 161     const point& p2,
 162     const point& pt
 163 )
 164 {
 165     scalar area01Pt = triPointRef(p0, p1, pt).normal() & n;
 166     scalar area12Pt = triPointRef(p1, p2, pt).normal() & n;
 167     scalar area20Pt = triPointRef(p2, p0, pt).normal() & n;
 168
 169     if ((area01Pt > 0) && (area12Pt > 0) && (area20Pt > 0))
 170     {
 171         return true;
 172     }
 173     else if ((area01Pt < 0) && (area12Pt < 0) && (area20Pt < 0))
 174     {
 175         FatalErrorIn("triangleContainsPoint") << abort(FatalError);
 176         return false;
 177     }
 178     else
 179     {
 180         return false;
 181     }
 182 }
 183
 184
 185 // Starting from startIndex find diagonal. Return in index1, index2.
 186 // Index1 always startIndex except when convex polygon
 187 void Foam::faceTriangulation::findDiagonal
 188 (
 189     const pointField& points,
 190     const face& f,
 191     const vectorField& edges,
 192     const vector& normal,
 193     const label startIndex,
 194     label& index1,
 195     label& index2
 196 )
 197 {
 198     const point& startPt = points[f[startIndex]];
 199
 200     // Calculate angle at startIndex
 201     const vector& rightE = edges[right(f.size(), startIndex)];
 202     const vector leftE = -edges[left(f.size(), startIndex)];
 203
 204     // Construct ray which bisects angle
 205     scalar cosHalfAngle = GREAT;
 206     scalar sinHalfAngle = GREAT;
 207     calcHalfAngle(normal, rightE, leftE, cosHalfAngle, sinHalfAngle);
 208     vector rayDir
 209     (
 210         cosHalfAngle*rightE
 211       + sinHalfAngle*(normal ^ rightE)
 212     );
 213     // rayDir should be normalized already but is not due to rounding errors
 214     // so normalize.
 215     rayDir /= mag(rayDir) + VSMALL;
 216
 217
 218     //
 219     // Check all edges (apart from rightE and leftE) for nearest intersection
 220     //
 221
 222     label faceVertI = f.fcIndex(startIndex);
 223
 224     pointHit minInter(false, vector::zero, GREAT, true);
 225     label minIndex = -1;
 226     scalar minPosOnEdge = GREAT;
 227
 228     for(label i = 0; i < f.size() - 2; i++)
 229     {
 230         scalar posOnEdge;
 231         pointHit inter =
 232             rayEdgeIntersect
 233             (
 234                 normal,
 235                 startPt,
 236                 rayDir,
 237                 points[f[faceVertI]],
 238                 points[f[f.fcIndex(faceVertI)]],
 239                 posOnEdge
 240             );
 241
 242         if (inter.hit() && inter.distance() < minInter.distance())
 243         {
 244             minInter = inter;
 245             minIndex = faceVertI;
 246             minPosOnEdge = posOnEdge;
 247         }
 248
 249         faceVertI = f.fcIndex(faceVertI);
 250     }
 251
 252
 253     if (minIndex == -1)
 254     {
 255         //WarningIn("faceTriangulation::findDiagonal")
 256         //    << "Could not find intersection starting from " << f[startIndex]
 257         //    << " for face " << f << endl;
 258
 259         index1 = -1;
 260         index2 = -1;
 261         return;
 262     }
 263
 264     const label leftIndex = minIndex;
 265     const label rightIndex = f.fcIndex(minIndex);
 266
 267     // Now ray intersects edge from leftIndex to rightIndex.
 268     // Check for intersection being one of the edge points. Make sure never
 269     // to return two consecutive points.
 270
 271     if (mag(minPosOnEdge) < edgeRelTol && f.fcIndex(startIndex) != leftIndex)
 272     {
 273         index1 = startIndex;
 274         index2 = leftIndex;
 275
 276         return;
 277     }
 278     if
 279     (
 280         mag(minPosOnEdge - 1) < edgeRelTol
 281      && f.fcIndex(rightIndex) != startIndex
 282     )
 283     {
 284         index1 = startIndex;
 285         index2 = rightIndex;
 286
 287         return;
 288     }
 289
 290     // Select visible vertex that minimizes
 291     // angle to bisection. Visibility checking by checking if inside triangle
 292     // formed by startIndex, leftIndex, rightIndex
 293
 294     const point& leftPt = points[f[leftIndex]];
 295     const point& rightPt = points[f[rightIndex]];
 296
 297     minIndex = -1;
 298     scalar maxCos = -GREAT;
 299
 300     // all vertices except for startIndex and ones to left and right of it.
 301     faceVertI = f.fcIndex(f.fcIndex(startIndex));
 302     for(label i = 0; i < f.size() - 3; i++)
 303     {
 304         const point& pt = points[f[faceVertI]];
 305
 306         if
 307         (
 308             (faceVertI == leftIndex)
 309          || (faceVertI == rightIndex)
 310          || (triangleContainsPoint(normal, startPt, leftPt, rightPt, pt))
 311         )
 312         {
 313             // pt inside triangle (so perhaps visible)
 314             // Select based on minimal angle (so guaranteed visible).
 315             vector edgePt0 = pt - startPt;
 316             edgePt0 /= mag(edgePt0);
 317
 318             scalar cos = rayDir & edgePt0;
 319             if (cos > maxCos)
 320             {
 321                 maxCos = cos;
 322                 minIndex = faceVertI;
 323             }
 324         }
 325         faceVertI = f.fcIndex(faceVertI);
 326     }
 327
 328     if (minIndex == -1)
 329     {
 330         // no vertex found. Return startIndex and one of the intersected edge
 331         // endpoints.
 332         index1 = startIndex;
 333
 334         if (f.fcIndex(startIndex) != leftIndex)
 335         {
 336             index2 = leftIndex;
 337         }
 338         else
 339         {
 340             index2 = rightIndex;
 341         }
 342
 343         return;
 344     }
 345
 346     index1 = startIndex;
 347     index2 = minIndex;
 348 }
 349
 350
 351 // Find label of vertex to start splitting from. Is:
 352 //     1] flattest concave angle
 353 //     2] flattest convex angle if no concave angles.
 354 Foam::label Foam::faceTriangulation::findStart
 355 (
 356     const face& f,
 357     const vectorField& edges,
 358     const vector& normal
 359 )
 360 {
 361     const label size = f.size();
 362
 363     scalar minCos = GREAT;
 364     label minIndex = -1;
 365
 366     forAll(f, fp)
 367     {
 368         const vector& rightEdge = edges[right(size, fp)];
 369         const vector leftEdge = -edges[left(size, fp)];
 370
 371         if (((rightEdge ^ leftEdge) & normal) < SMALL)
 372         {
 373             scalar cos = rightEdge & leftEdge;
 374             if (cos < minCos)
 375             {
 376                 minCos = cos;
 377                 minIndex = fp;
 378             }
 379         }
 380     }
 381
 382     if (minIndex == -1)
 383     {
 384         // No concave angle found. Get flattest convex angle
 385         minCos = GREAT;
 386
 387         forAll(f, fp)
 388         {
 389             const vector& rightEdge = edges[right(size, fp)];
 390             const vector leftEdge = -edges[left(size, fp)];
 391
 392             scalar cos = rightEdge & leftEdge;
 393             if (cos < minCos)
 394             {
 395                 minCos = cos;
 396                 minIndex = fp;
 397             }
 398         }
 399     }
 400
 401     return minIndex;
 402 }
 403
 404
 405 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 406
 407 // Split face f into triangles. Handles all simple (convex & concave)
 408 // polygons.
 409 bool Foam::faceTriangulation::split
 410 (
 411     const bool fallBack,
 412     const pointField& points,
 413     const face& f,
 414     const vector& normal,
 415     label& triI
 416 )
 417 {
 418     const label size = f.size();
 419
 420     if (size <= 2)
 421     {
 422         WarningIn
 423         (
 424             "split(const bool, const pointField&, const face&"
 425             ", const vector&, label&)"
 426         )   << "Illegal face:" << f
 427             << " with points " << IndirectList<point>(points, f)()
 428             << endl;
 429
 430         return false;
 431     }
 432     else if (size == 3)
 433     {
 434         // Triangle. Just copy.
 435         triFace& tri = operator[](triI++);
 436         tri[0] = f[0];
 437         tri[1] = f[1];
 438         tri[2] = f[2];
 439
 440         return true;
 441     }
 442     else
 443     {
 444         // General case. Start splitting for -flattest concave angle
 445         // -or flattest convex angle if no concave angles.
 446
 447         tmp<vectorField> tedges(calcEdges(f, points));
 448         const vectorField& edges = tedges();
 449
 450         label startIndex = findStart(f, edges, normal);
 451
 452         // Find diagonal to split face across
 453         label index1 = -1;
 454         label index2 = -1;
 455
 456         for (label iter = 0; iter < f.size(); iter++)
 457         {
 458             findDiagonal
 459             (
 460                 points,
 461                 f,
 462                 edges,
 463                 normal,
 464                 startIndex,
 465                 index1,
 466                 index2
 467             );
 468
 469             if (index1 != -1 && index2 != -1)
 470             {
 471                 // Found correct diagonal
 472                 break;
 473             }
 474
 475             // Try splitting from next startingIndex.
 476             startIndex = f.fcIndex(startIndex);
 477         }
 478
 479         if (index1 == -1 || index2 == -1)
 480         {
 481             if (fallBack)
 482             {
 483                 // Do naive triangulation. Find smallest angle to start
 484                 // triangulating from.
 485                 label maxIndex = -1;
 486                 scalar maxCos = -GREAT;
 487
 488                 forAll(f, fp)
 489                 {
 490                     const vector& rightEdge = edges[right(size, fp)];
 491                     const vector leftEdge = -edges[left(size, fp)];
 492
 493                     scalar cos = rightEdge & leftEdge;
 494                     if (cos > maxCos)
 495                     {
 496                         maxCos = cos;
 497                         maxIndex = fp;
 498                     }
 499                 }
 500
 501                 WarningIn
 502                 (
 503                     "split(const bool, const pointField&, const face&"
 504                     ", const vector&, label&)"
 505                 )   << "Cannot find valid diagonal on face " << f
 506                     << " with points " << IndirectList<point>(points, f)() << nl
 507                     << "Returning naive triangulation starting from "
 508                     << f[maxIndex] << " which might not be correct for a"
 509                     << " concave or warped face" << endl;
 510
 511
 512                 label fp = f.fcIndex(maxIndex);
 513
 514                 for (label i = 0; i < size-2; i++)
 515                 {
 516                     label nextFp = f.fcIndex(fp);
 517
 518                     triFace& tri = operator[](triI++);
 519                     tri[0] = f[maxIndex];
 520                     tri[1] = f[fp];
 521                     tri[2] = f[nextFp];
 522
 523                     fp = nextFp;
 524                 }
 525
 526                 return true;
 527             }
 528             else
 529             {
 530                 WarningIn
 531                 (
 532                     "split(const bool, const pointField&, const face&"
 533                     ", const vector&, label&)"
 534                 )   << "Cannot find valid diagonal on face " << f
 535                     << " with points " << IndirectList<point>(points, f)() << nl
 536                     << "Returning empty triFaceList" << endl;
 537
 538                 return false;
 539             }
 540         }
 541
 542
 543         // Split into two subshapes.
 544         //     face1: index1 to index2
 545         //     face2: index2 to index1
 546
 547         // Get sizes of the two subshapes
 548         label diff = 0;
 549         if (index2 > index1)
 550         {
 551             diff = index2 - index1;
 552         }
 553         else
 554         {
 555             // folded round
 556             diff = index2 + size - index1;
 557         }
 558
 559         label nPoints1 = diff + 1;
 560         label nPoints2 = size - diff + 1;
 561
 562         if (nPoints1 == size || nPoints2 == size)
 563         {
 564             FatalErrorIn
 565             (
 566                 "split(const bool, const pointField&, const face&"
 567                 ", const vector&, label&)"
 568             )   << "Illegal split of face:" << f
 569                 << " with points " << IndirectList<point>(points, f)()
 570                 << " at indices " << index1 << " and " << index2
 571                 << abort(FatalError);
 572         }
 573
 574
 575         // Collect face1 points
 576         face face1(nPoints1);
 577
 578         label faceVertI = index1;
 579         for (int i = 0; i < nPoints1; i++)
 580         {
 581             face1[i] = f[faceVertI];
 582             faceVertI = f.fcIndex(faceVertI);
 583         }
 584
 585         // Collect face2 points
 586         face face2(nPoints2);
 587
 588         faceVertI = index2;
 589         for (int i = 0; i < nPoints2; i++)
 590         {
 591             face2[i] = f[faceVertI];
 592             faceVertI = f.fcIndex(faceVertI);
 593         }
 594
 595         // Decompose the split faces
 596         //Pout<< "Split face:" << f << " into " << face1 << " and " << face2
 597         //    << endl;
 598         //string oldPrefix(Pout.prefix());
 599         //Pout.prefix() = "  " + oldPrefix;
 600
 601         bool splitOk =
 602             split(fallBack, points, face1, normal, triI)
 603          && split(fallBack, points, face2, normal, triI);
 604
 605         //Pout.prefix() = oldPrefix;
 606
 607         return splitOk;
 608     }
 609 }
 610
 611
 612 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 613
 614 // Null constructor
 615 Foam::faceTriangulation::faceTriangulation()
 616 :
 617     triFaceList()
 618 {}
 619
 620
 621 // Construct from components
 622 Foam::faceTriangulation::faceTriangulation
 623 (
 624     const pointField& points,
 625     const face& f,
 626     const bool fallBack
 627 )
 628 :
 629     triFaceList(f.size()-2)
 630 {
 631     vector avgNormal = f.normal(points);
 632     avgNormal /= mag(avgNormal) + VSMALL;
 633
 634     label triI = 0;
 635
 636     bool valid = split(fallBack, points, f, avgNormal, triI);
 637
 638     if (!valid)
 639     {
 640         setSize(0);
 641     }
 642 }
 643
 644
 645 // Construct from components
 646 Foam::faceTriangulation::faceTriangulation
 647 (
 648     const pointField& points,
 649     const face& f,
 650     const vector& n,
 651     const bool fallBack
 652 )
 653 :
 654     triFaceList(f.size()-2)
 655 {
 656     label triI = 0;
 657
 658     bool valid = split(fallBack, points, f, n, triI);
 659
 660     if (!valid)
 661     {
 662         setSize(0);
 663     }
 664 }
 665
 666
 667 // Construct from Istream
 668 Foam::faceTriangulation::faceTriangulation(Istream& is)
 669 :
 670     triFaceList(is)
 671 {}
 672
 673
 674 // ************************************************************************* //