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