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-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 "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 ? i-1 : size-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         point thisPt = points[f[i]];
  64         point nextPt = points[f[f.fcIndex(i)]];
  65
  66         vector vec(nextPt - thisPt);
  67         vec /= mag(vec) + VSMALL;
  68
  69         edges[i] = vec;
  70     }
  71
  72     return tedges;
  73 }
  74
  75
  76 // Calculates half angle components of angle from e0 to e1
  77 void Foam::faceTriangulation::calcHalfAngle
  78 (
  79     const vector& normal,
  80     const vector& e0,
  81     const vector& e1,
  82     scalar& cosHalfAngle,
  83     scalar& sinHalfAngle
  84 )
  85 {
  86     // truncate cos to +-1 to prevent negative numbers
  87     scalar cos = max(-1, min(1, e0 & e1));
  88
  89     scalar sin = (e0 ^ e1) & normal;
  90
  91     if (sin < -ROOTVSMALL)
  92     {
  93         // 3rd or 4th quadrant
  94         cosHalfAngle = - Foam::sqrt(0.5*(1 + cos));
  95         sinHalfAngle = Foam::sqrt(0.5*(1 - cos));
  96     }
  97     else
  98     {
  99         // 1st or 2nd quadrant
 100         cosHalfAngle = Foam::sqrt(0.5*(1 + cos));
 101         sinHalfAngle = Foam::sqrt(0.5*(1 - cos));
 102     }
 103 }
 104
 105
 106 // Calculate intersection point between edge p1-p2 and ray (in 2D).
 107 // Return true and intersection point if intersection between p1 and p2.
 108 Foam::pointHit Foam::faceTriangulation::rayEdgeIntersect
 109 (
 110     const vector& normal,
 111     const point& rayOrigin,
 112     const vector& rayDir,
 113     const point& p1,
 114     const point& p2,
 115     scalar& posOnEdge
 116 )
 117 {
 118     // Start off from miss
 119     pointHit result(p1);
 120
 121     // Construct plane normal to rayDir and intersect
 122     const vector y = normal ^ rayDir;
 123
 124     posOnEdge = plane(rayOrigin, y).normalIntersect(p1, (p2-p1));
 125
 126     // Check intersection to left of p1 or right of p2
 127     if ((posOnEdge < 0) || (posOnEdge > 1))
 128     {
 129         // Miss
 130     }
 131     else
 132     {
 133         // Check intersection behind rayOrigin
 134         point intersectPt = p1 + posOnEdge * (p2 - p1);
 135
 136         if (((intersectPt - rayOrigin) & rayDir) < 0)
 137         {
 138             // Miss
 139         }
 140         else
 141         {
 142             // Hit
 143             result.setHit();
 144             result.setPoint(intersectPt);
 145             result.setDistance(mag(intersectPt - rayOrigin));
 146         }
 147     }
 148     return result;
 149 }
 150
 151
 152 // Return true if triangle given its three points (anticlockwise ordered)
 153 // contains point
 154 bool Foam::faceTriangulation::triangleContainsPoint
 155 (
 156     const vector& n,
 157     const point& p0,
 158     const point& p1,
 159     const point& p2,
 160     const point& pt
 161 )
 162 {
 163     scalar area01Pt = triPointRef(p0, p1, pt).normal() & n;
 164     scalar area12Pt = triPointRef(p1, p2, pt).normal() & n;
 165     scalar area20Pt = triPointRef(p2, p0, pt).normal() & n;
 166
 167     if ((area01Pt > 0) && (area12Pt > 0) && (area20Pt > 0))
 168     {
 169         return true;
 170     }
 171     else if ((area01Pt < 0) && (area12Pt < 0) && (area20Pt < 0))
 172     {
 173         FatalErrorIn("triangleContainsPoint") << abort(FatalError);
 174         return false;
 175     }
 176     else
 177     {
 178         return false;
 179     }
 180 }
 181
 182
 183 // Starting from startIndex find diagonal. Return in index1, index2.
 184 // Index1 always startIndex except when convex polygon
 185 void Foam::faceTriangulation::findDiagonal
 186 (
 187     const pointField& points,
 188     const face& f,
 189     const vectorField& edges,
 190     const vector& normal,
 191     const label startIndex,
 192     label& index1,
 193     label& index2
 194 )
 195 {
 196     const point& startPt = points[f[startIndex]];
 197
 198     // Calculate angle at startIndex
 199     const vector& rightE = edges[right(f.size(), startIndex)];
 200     const vector leftE = -edges[left(f.size(), startIndex)];
 201
 202     // Construct ray which bisects angle
 203     scalar cosHalfAngle = GREAT;
 204     scalar sinHalfAngle = GREAT;
 205     calcHalfAngle(normal, rightE, leftE, cosHalfAngle, sinHalfAngle);
 206     vector rayDir
 207     (
 208         cosHalfAngle*rightE
 209       + sinHalfAngle*(normal ^ rightE)
 210     );
 211     // rayDir should be normalized already but is not due to rounding errors
 212     // so normalize.
 213     rayDir /= mag(rayDir) + VSMALL;
 214
 215
 216     //
 217     // Check all edges (apart from rightE and leftE) for nearest intersection
 218     //
 219
 220     label faceVertI = f.fcIndex(startIndex);
 221
 222     pointHit minInter(false, vector::zero, GREAT, true);
 223     label minIndex = -1;
 224     scalar minPosOnEdge = GREAT;
 225
 226     for (label i = 0; i < f.size() - 2; i++)
 227     {
 228         scalar posOnEdge;
 229         pointHit inter =
 230             rayEdgeIntersect
 231             (
 232                 normal,
 233                 startPt,
 234                 rayDir,
 235                 points[f[faceVertI]],
 236                 points[f[f.fcIndex(faceVertI)]],
 237                 posOnEdge
 238             );
 239
 240         if (inter.hit() && inter.distance() < minInter.distance())
 241         {
 242             minInter = inter;
 243             minIndex = faceVertI;
 244             minPosOnEdge = posOnEdge;
 245         }
 246
 247         faceVertI = f.fcIndex(faceVertI);
 248     }
 249
 250
 251     if (minIndex == -1)
 252     {
 253         //WarningIn("faceTriangulation::findDiagonal")
 254         //    << "Could not find intersection starting from " << f[startIndex]
 255         //    << " for face " << f << endl;
 256
 257         index1 = -1;
 258         index2 = -1;
 259         return;
 260     }
 261
 262     const label leftIndex = minIndex;
 263     const label rightIndex = f.fcIndex(minIndex);
 264
 265     // Now ray intersects edge from leftIndex to rightIndex.
 266     // Check for intersection being one of the edge points. Make sure never
 267     // to return two consecutive points.
 268
 269     if (mag(minPosOnEdge) < edgeRelTol && f.fcIndex(startIndex) != leftIndex)
 270     {
 271         index1 = startIndex;
 272         index2 = leftIndex;
 273
 274         return;
 275     }
 276     if
 277     (
 278         mag(minPosOnEdge - 1) < edgeRelTol
 279      && f.fcIndex(rightIndex) != startIndex
 280     )
 281     {
 282         index1 = startIndex;
 283         index2 = rightIndex;
 284
 285         return;
 286     }
 287
 288     // Select visible vertex that minimizes
 289     // angle to bisection. Visibility checking by checking if inside triangle
 290     // formed by startIndex, leftIndex, rightIndex
 291
 292     const point& leftPt = points[f[leftIndex]];
 293     const point& rightPt = points[f[rightIndex]];
 294
 295     minIndex = -1;
 296     scalar maxCos = -GREAT;
 297
 298     // all vertices except for startIndex and ones to left and right of it.
 299     faceVertI = f.fcIndex(f.fcIndex(startIndex));
 300     for (label i = 0; i < f.size() - 3; i++)
 301     {
 302         const point& pt = points[f[faceVertI]];
 303
 304         if
 305         (
 306             (faceVertI == leftIndex)
 307          || (faceVertI == rightIndex)
 308          || (triangleContainsPoint(normal, startPt, leftPt, rightPt, pt))
 309         )
 310         {
 311             // pt inside triangle (so perhaps visible)
 312             // Select based on minimal angle (so guaranteed visible).
 313             vector edgePt0 = pt - startPt;
 314             edgePt0 /= mag(edgePt0);
 315
 316             scalar cos = rayDir & edgePt0;
 317             if (cos > maxCos)
 318             {
 319                 maxCos = cos;
 320                 minIndex = faceVertI;
 321             }
 322         }
 323         faceVertI = f.fcIndex(faceVertI);
 324     }
 325
 326     if (minIndex == -1)
 327     {
 328         // no vertex found. Return startIndex and one of the intersected edge
 329         // endpoints.
 330         index1 = startIndex;
 331
 332         if (f.fcIndex(startIndex) != leftIndex)
 333         {
 334             index2 = leftIndex;
 335         }
 336         else
 337         {
 338             index2 = rightIndex;
 339         }
 340
 341         return;
 342     }
 343
 344     index1 = startIndex;
 345     index2 = minIndex;
 346 }
 347
 348
 349 // Find label of vertex to start splitting from. Is:
 350 //     1] flattest concave angle
 351 //     2] flattest convex angle if no concave angles.
 352 Foam::label Foam::faceTriangulation::findStart
 353 (
 354     const face& f,
 355     const vectorField& edges,
 356     const vector& normal
 357 )
 358 {
 359     const label size = f.size();
 360
 361     scalar minCos = GREAT;
 362     label minIndex = -1;
 363
 364     forAll(f, fp)
 365     {
 366         const vector& rightEdge = edges[right(size, fp)];
 367         const vector leftEdge = -edges[left(size, fp)];
 368
 369         if (((rightEdge ^ leftEdge) & normal) < ROOTVSMALL)
 370         {
 371             scalar cos = rightEdge & leftEdge;
 372             if (cos < minCos)
 373             {
 374                 minCos = cos;
 375                 minIndex = fp;
 376             }
 377         }
 378     }
 379
 380     if (minIndex == -1)
 381     {
 382         // No concave angle found. Get flattest convex angle
 383         minCos = GREAT;
 384
 385         forAll(f, fp)
 386         {
 387             const vector& rightEdge = edges[right(size, fp)];
 388             const vector leftEdge = -edges[left(size, fp)];
 389
 390             scalar cos = rightEdge & leftEdge;
 391             if (cos < minCos)
 392             {
 393                 minCos = cos;
 394                 minIndex = fp;
 395             }
 396         }
 397     }
 398
 399     return minIndex;
 400 }
 401
 402
 403 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 404
 405 // Split face f into triangles. Handles all simple (convex & concave)
 406 // polygons.
 407 bool Foam::faceTriangulation::split
 408 (
 409     const bool fallBack,
 410     const pointField& points,
 411     const face& f,
 412     const vector& normal,
 413     label& triI
 414 )
 415 {
 416     const label size = f.size();
 417
 418     if (size <= 2)
 419     {
 420         WarningIn
 421         (
 422             "split(const bool, const pointField&, const face&"
 423             ", const vector&, label&)"
 424         )   << "Illegal face:" << f
 425             << " with points " << UIndirectList<point>(points, f)()
 426             << endl;
 427
 428         return false;
 429     }
 430     else if (size == 3)
 431     {
 432         // Triangle. Just copy.
 433         triFace& tri = operator[](triI++);
 434         tri[0] = f[0];
 435         tri[1] = f[1];
 436         tri[2] = f[2];
 437
 438         return true;
 439     }
 440     else
 441     {
 442         // General case. Start splitting for -flattest concave angle
 443         // -or flattest convex angle if no concave angles.
 444
 445         tmp<vectorField> tedges(calcEdges(f, points));
 446         const vectorField& edges = tedges();
 447
 448         label startIndex = findStart(f, edges, normal);
 449
 450         // Find diagonal to split face across
 451         label index1 = -1;
 452         label index2 = -1;
 453
 454         for (label iter = 0; iter < f.size(); iter++)
 455         {
 456             findDiagonal
 457             (
 458                 points,
 459                 f,
 460                 edges,
 461                 normal,
 462                 startIndex,
 463                 index1,
 464                 index2
 465             );
 466
 467             if (index1 != -1 && index2 != -1)
 468             {
 469                 // Found correct diagonal
 470                 break;
 471             }
 472
 473             // Try splitting from next startingIndex.
 474             startIndex = f.fcIndex(startIndex);
 475         }
 476
 477         if (index1 == -1 || index2 == -1)
 478         {
 479             if (fallBack)
 480             {
 481                 // Do naive triangulation. Find smallest angle to start
 482                 // triangulating from.
 483                 label maxIndex = -1;
 484                 scalar maxCos = -GREAT;
 485
 486                 forAll(f, fp)
 487                 {
 488                     const vector& rightEdge = edges[right(size, fp)];
 489                     const vector leftEdge = -edges[left(size, fp)];
 490
 491                     scalar cos = rightEdge & leftEdge;
 492                     if (cos > maxCos)
 493                     {
 494                         maxCos = cos;
 495                         maxIndex = fp;
 496                     }
 497                 }
 498
 499                 WarningIn
 500                 (
 501                     "split(const bool, const pointField&, const face&"
 502                     ", const vector&, label&)"
 503                 )   << "Cannot find valid diagonal on face " << f
 504                     << " with points " << UIndirectList<point>(points, f)()
 505                     << nl
 506                     << "Returning naive triangulation starting from "
 507                     << f[maxIndex] << " which might not be correct for a"
 508                     << " concave or warped face" << endl;
 509
 510
 511                 label fp = f.fcIndex(maxIndex);
 512
 513                 for (label i = 0; i < size-2; i++)
 514                 {
 515                     label nextFp = f.fcIndex(fp);
 516
 517                     triFace& tri = operator[](triI++);
 518                     tri[0] = f[maxIndex];
 519                     tri[1] = f[fp];
 520                     tri[2] = f[nextFp];
 521
 522                     fp = nextFp;
 523                 }
 524
 525                 return true;
 526             }
 527             else
 528             {
 529                 WarningIn
 530                 (
 531                     "split(const bool, const pointField&, const face&"
 532                     ", const vector&, label&)"
 533                 )   << "Cannot find valid diagonal on face " << f
 534                     << " with points " << UIndirectList<point>(points, f)()
 535                     << 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 " << UIndirectList<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 // ************************************************************************* //