src/OpenFOAM/meshes/meshShapes/face/face.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 "face.H"
  28 #include "triPointRef.H"
  29 #include "mathematicalConstants.H"
  30
  31 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  32
  33 const char* const Foam::face::typeName = "face";
  34
  35 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  36
  37 Foam::tmp<Foam::vectorField>
  38 Foam::face::calcEdges(const pointField& points) const
  39 {
  40     tmp<vectorField> tedges(new vectorField(size()));
  41     vectorField& edges = tedges();
  42
  43     forAll(*this, i)
  44     {
  45         label ni = fcIndex(i);
  46
  47         point thisPt = points[operator[](i)];
  48         point nextPt = points[operator[](ni)];
  49
  50         vector vec(nextPt - thisPt);
  51         vec /= Foam::mag(vec) + VSMALL;
  52
  53         edges[i] = vec;
  54     }
  55
  56     return tedges;
  57 }
  58
  59
  60 Foam::scalar Foam::face::edgeCos
  61 (
  62     const vectorField& edges,
  63     const label index
  64 ) const
  65 {
  66     label leftEdgeI = left(index);
  67     label rightEdgeI = right(index);
  68
  69     // note negate on left edge to get correct left-pointing edge.
  70     return -(edges[leftEdgeI] & edges[rightEdgeI]);
  71 }
  72
  73
  74 Foam::label Foam::face::mostConcaveAngle
  75 (
  76     const pointField& points,
  77     const vectorField& edges,
  78     scalar& maxAngle
  79 ) const
  80 {
  81     vector n(normal(points));
  82
  83     label index = 0;
  84     maxAngle = -GREAT;
  85
  86     forAll(edges, i)
  87     {
  88         label leftEdgeI = left(i);
  89         label rightEdgeI = right(i);
  90
  91         vector edgeNormal = edges[rightEdgeI] ^ edges[leftEdgeI];
  92
  93         scalar edgeCos = edges[leftEdgeI] & edges[rightEdgeI];
  94         scalar edgeAngle = acos(max(-1.0, min(1.0, edgeCos)));
  95
  96         scalar angle;
  97
  98         if ((edgeNormal & n) > 0)
  99         {
 100             // Concave angle.
 101             angle = mathematicalConstant::pi + edgeAngle;
 102         }
 103         else
 104         {
 105             // Convex angle. Note '-' to take into account that rightEdge
 106             // and leftEdge are head-to-tail connected.
 107             angle = mathematicalConstant::pi - edgeAngle;
 108         }
 109
 110         if (angle > maxAngle)
 111         {
 112             maxAngle = angle;
 113             index = i;
 114         }
 115     }
 116
 117     return index;
 118 }
 119
 120
 121 void Foam::face::split
 122 (
 123     const face::splitMode mode,
 124     const pointField& points,
 125     label& triI,
 126     label& quadI,
 127     faceList& triFaces,
 128     faceList& quadFaces
 129 ) const
 130 {
 131     if (size() <= 2)
 132     {
 133         FatalErrorIn
 134         (
 135             "face::split"
 136             "(const face::splitMode, const pointField&, label&, label&"
 137             ", faceList&, faceList&)"
 138         )
 139             << "Serious problem: asked to split a face with < 3 vertices"
 140             << abort(FatalError);
 141     }
 142     if (size() == 3)
 143     {
 144         // Triangle. Just copy.
 145         if ((mode == COUNTQUAD) || (mode == COUNTTRIANGLE))
 146         {
 147             triI++;
 148         }
 149         else
 150         {
 151             triFaces[triI++] = *this;
 152         }
 153     }
 154     else if (size() == 4)
 155     {
 156         if (mode == COUNTTRIANGLE)
 157         {
 158             triI += 2;
 159         }
 160         if (mode == COUNTQUAD)
 161         {
 162             quadI++;
 163         }
 164         else if (mode == SPLITTRIANGLE)
 165         {
 166             //  Start at point with largest internal angle.
 167             const vectorField edges(calcEdges(points));
 168
 169             scalar minAngle;
 170             label startIndex = mostConcaveAngle(points, edges, minAngle);
 171
 172             label nextIndex = fcIndex(startIndex);
 173             label splitIndex = fcIndex(nextIndex);
 174
 175             // Create triangles
 176             face triFace(3);
 177             triFace[0] = operator[](startIndex);
 178             triFace[1] = operator[](nextIndex);
 179             triFace[2] = operator[](splitIndex);
 180
 181             triFaces[triI++] = triFace;
 182
 183             triFace[0] = operator[](splitIndex);
 184             triFace[1] = operator[](fcIndex(splitIndex));
 185             triFace[2] = operator[](startIndex);
 186
 187             triFaces[triI++] = triFace;
 188         }
 189         else
 190         {
 191             quadFaces[quadI++] = *this;
 192         }
 193     }
 194     else
 195     {
 196         // General case. Like quad: search for largest internal angle.
 197
 198         const vectorField edges(calcEdges(points));
 199
 200         scalar minAngle = 1;
 201         label startIndex = mostConcaveAngle(points, edges, minAngle);
 202
 203         scalar bisectAngle = minAngle/2;
 204         vector rightEdge = edges[right(startIndex)];
 205
 206         //
 207         // Look for opposite point which as close as possible bisects angle
 208         //
 209
 210         // split candidate starts two points away.
 211         label index = fcIndex(fcIndex(startIndex));
 212
 213         label minIndex = index;
 214         scalar minDiff = Foam::mathematicalConstant::pi;
 215
 216         for(label i = 0; i < size() - 3; i++)
 217         {
 218             vector splitEdge
 219             (
 220                 points[operator[](index)]
 221               - points[operator[](startIndex)]
 222             );
 223             splitEdge /= Foam::mag(splitEdge) + VSMALL;
 224
 225             const scalar splitCos = splitEdge & rightEdge;
 226             const scalar splitAngle = acos(max(-1.0, min(1.0, splitCos)));
 227             const scalar angleDiff = fabs(splitAngle - bisectAngle);
 228
 229             if (angleDiff < minDiff)
 230             {
 231                 minDiff = angleDiff;
 232                 minIndex = index;
 233             }
 234
 235             // go to next candidate
 236             index = fcIndex(index);
 237         }
 238
 239
 240         // Split into two subshapes.
 241         //     face1: startIndex to minIndex
 242         //     face2: minIndex to startIndex
 243
 244         // Get sizes of the two subshapes
 245         label diff = 0;
 246         if (minIndex > startIndex)
 247         {
 248             diff = minIndex - startIndex;
 249         }
 250         else
 251         {
 252             // folded round
 253             diff = minIndex + size() - startIndex;
 254         }
 255
 256         label nPoints1 = diff + 1;
 257         label nPoints2 = size() - diff + 1;
 258
 259         // Collect face1 points
 260         face face1(nPoints1);
 261
 262         index = startIndex;
 263         for (label i = 0; i < nPoints1; i++)
 264         {
 265             face1[i] = operator[](index);
 266             index = fcIndex(index);
 267         }
 268
 269         // Collect face2 points
 270         face face2(nPoints2);
 271
 272         index = minIndex;
 273         for (label i = 0; i < nPoints2; i++)
 274         {
 275             face2[i] = operator[](index);
 276             index = fcIndex(index);
 277         }
 278
 279         // Split faces
 280         face1.split(mode, points, triI, quadI, triFaces, quadFaces);
 281         face2.split(mode, points, triI, quadI, triFaces, quadFaces);
 282     }
 283 }
 284
 285
 286 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
 287
 288
 289 // return
 290 //   0: no match
 291 //  +1: identical
 292 //  -1: same face, but different orientation
 293 int Foam::face::compare(const face& a, const face& b)
 294 {
 295     // Basic rule: we assume that the sequence of labels in each list
 296     // will be circular in the same order (but not necessarily in the
 297     // same direction or from the same starting point).
 298
 299     // Trivial reject: faces are different size
 300     label sizeA = a.size();
 301     label sizeB = b.size();
 302
 303     if (sizeA != sizeB)
 304     {
 305         return 0;
 306     }
 307
 308
 309     // Full list comparison
 310     const label firstA = a[0];
 311     label Bptr = -1;
 312
 313     forAll (b, i)
 314     {
 315         if (b[i] == firstA)
 316         {
 317             Bptr = i;        // 'found match' at element 'i'
 318             break;
 319         }
 320     }
 321
 322     // If no match was found, return 0
 323     if (Bptr < 0)
 324     {
 325         return 0;
 326     }
 327
 328     // Now we must look for the direction, if any
 329     label secondA = a[1];
 330
 331     if (sizeA > 1 && (secondA == firstA || firstA == a[sizeA - 1]))
 332     {
 333         face ca = a;
 334         ca.collapse();
 335
 336         face cb = b;
 337         cb.collapse();
 338
 339         return face::compare(ca, cb);
 340     }
 341
 342     int dir = 0;
 343
 344     // Check whether at top of list
 345     Bptr++;
 346     if (Bptr == b.size())
 347     {
 348         Bptr = 0;
 349     }
 350
 351     // Test whether upward label matches second A label
 352     if (b[Bptr] == secondA)
 353     {
 354         // Yes - direction is 'up'
 355         dir = 1;
 356     }
 357     else
 358     {
 359         // No - so look downwards, checking whether at bottom of list
 360         Bptr -= 2;
 361
 362         if (Bptr < 0)
 363         {
 364             // wraparound
 365             Bptr += b.size();
 366         }
 367
 368         // Test whether downward label matches second A label
 369         if (b[Bptr] == secondA)
 370         {
 371             // Yes - direction is 'down'
 372             dir = -1;
 373         }
 374     }
 375
 376     // Check whether a match was made at all, and exit 0 if not
 377     if (dir == 0)
 378     {
 379         return 0;
 380     }
 381
 382     // Decrement size by 2 to account for first searches
 383     sizeA -= 2;
 384
 385     // We now have both direction of search and next element
 386     // to search, so we can continue search until no more points.
 387     label Aptr = 1;
 388     if (dir > 0)
 389     {
 390         while (sizeA--)
 391         {
 392             Aptr++;
 393             if (Aptr >= a.size())
 394             {
 395                 Aptr = 0;
 396             }
 397
 398             Bptr++;
 399             if (Bptr >= b.size())
 400             {
 401                 Bptr = 0;
 402             }
 403
 404             if (a[Aptr] != b[Bptr])
 405             {
 406                 return 0;
 407             }
 408         }
 409     }
 410     else
 411     {
 412         while (sizeA--)
 413         {
 414             Aptr++;
 415             if (Aptr >= a.size())
 416             {
 417                 Aptr = 0;
 418             }
 419
 420             Bptr--;
 421             if (Bptr < 0)
 422             {
 423                 Bptr = b.size() - 1;
 424             }
 425
 426             if (a[Aptr] != b[Bptr])
 427             {
 428                 return 0;
 429             }
 430         }
 431     }
 432
 433     // They must be equal - return direction
 434     return dir;
 435 }
 436
 437
 438 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 439
 440
 441 Foam::label Foam::face::collapse()
 442 {
 443     if (size() > 1)
 444     {
 445         label ci = 0;
 446         for (label i=1; i<size(); i++)
 447         {
 448             if (operator[](i) != operator[](ci))
 449             {
 450                 operator[](++ci) = operator[](i);
 451             }
 452         }
 453
 454         if (operator[](ci) != operator[](0))
 455         {
 456             ci++;
 457         }
 458
 459         setSize(ci);
 460     }
 461
 462     return size();
 463 }
 464
 465
 466 Foam::point Foam::face::centre(const pointField& meshPoints) const
 467 {
 468     // Calculate the centre by breaking the face into triangles and
 469     // area-weighted averaging their centres
 470
 471     // If the face is a triangle, do a direct calculation
 472     if (size() == 3)
 473     {
 474         return
 475             (1.0/3.0)
 476            *(
 477                meshPoints[operator[](0)]
 478              + meshPoints[operator[](1)]
 479              + meshPoints[operator[](2)]
 480             );
 481     }
 482
 483     label nPoints = size();
 484
 485     point centrePoint = point::zero;
 486     for (register label pI=0; pI<nPoints; pI++)
 487     {
 488         centrePoint += meshPoints[operator[](pI)];
 489     }
 490     centrePoint /= nPoints;
 491
 492     scalar sumA = 0;
 493     vector sumAc = vector::zero;
 494
 495     for (register label pI=0; pI<nPoints; pI++)
 496     {
 497         const point& nextPoint = meshPoints[operator[]((pI + 1) % nPoints)];
 498
 499         // Calculate 3*triangle centre
 500         vector ttc
 501         (
 502             meshPoints[operator[](pI)]
 503           + nextPoint
 504           + centrePoint
 505         );
 506
 507         // Calculate 2*triangle area
 508         scalar ta = Foam::mag
 509         (
 510             (meshPoints[operator[](pI)] - centrePoint)
 511           ^ (nextPoint - centrePoint)
 512         );
 513
 514         sumA += ta;
 515         sumAc += ta*ttc;
 516     }
 517
 518     if (sumA > VSMALL)
 519     {
 520         return sumAc/(3*sumA);
 521     }
 522     else
 523     {
 524         return centrePoint;
 525     }
 526 }
 527
 528
 529 Foam::vector Foam::face::normal(const pointField& meshPoints) const
 530 {
 531     // Calculate the normal by summing the face triangle normals.
 532     // Changed to deal with small concavity by using a central decomposition
 533     //
 534
 535     // If the face is a triangle, do a direct calculation to avoid round-off
 536     // error-related problems
 537     //
 538     if (size() == 3)
 539     {
 540         return triPointRef
 541         (
 542             meshPoints[operator[](0)],
 543             meshPoints[operator[](1)],
 544             meshPoints[operator[](2)]
 545         ).normal();
 546     }
 547
 548     vector n = vector::zero;
 549
 550     point centrePoint = Foam::average(points(meshPoints));
 551
 552     label nPoints = size();
 553
 554     point nextPoint = centrePoint;
 555
 556     register label pI;
 557
 558     for (pI = 0; pI < nPoints; pI++)
 559     {
 560         if (pI < nPoints - 1)
 561         {
 562             nextPoint = meshPoints[operator[](pI + 1)];
 563         }
 564         else
 565         {
 566             nextPoint = meshPoints[operator[](0)];
 567         }
 568
 569         // Note: for best accuracy, centre point always comes last
 570         //
 571         n += triPointRef
 572         (
 573             meshPoints[operator[](pI)],
 574             nextPoint,
 575             centrePoint
 576         ).normal();
 577     }
 578
 579     return n;
 580 }
 581
 582
 583 Foam::face Foam::face::reverseFace() const
 584 {
 585     // reverse the label list and return
 586     // Changed to make sure that the starting point of the original
 587     // and the reverse face is identical.
 588     //
 589
 590     const labelList& myList = *this;
 591
 592     labelList newList(size());
 593
 594     newList[0] = myList[0];
 595
 596     for (label pointI = 1; pointI < newList.size(); pointI++)
 597     {
 598         newList[pointI] = myList[size() - pointI];
 599     }
 600
 601     return face(newList);
 602 }
 603
 604
 605 Foam::label Foam::face::which(const label globalIndex) const
 606 {
 607     label pointInFace = -1;
 608     const labelList& f = *this;
 609
 610     forAll (f, i)
 611     {
 612         if (f[i] == globalIndex)
 613         {
 614             pointInFace = i;
 615             break;
 616         }
 617     }
 618
 619     return pointInFace;
 620 }
 621
 622
 623 Foam::scalar Foam::face::sweptVol
 624 (
 625     const pointField& oldPoints,
 626     const pointField& newPoints
 627 ) const
 628 {
 629     scalar sv = 0;
 630
 631     // Calculate the swept volume by breaking the face into triangles and
 632     // summing their swept volumes.
 633     // Changed to deal with small concavity by using a central decomposition
 634
 635     point centreOldPoint = centre(oldPoints);
 636     point centreNewPoint = centre(newPoints);
 637
 638     label nPoints = size();
 639
 640     point nextOldPoint = centreOldPoint;
 641
 642     point nextNewPoint = centreNewPoint;
 643
 644     register label pI;
 645
 646     for (pI = 0; pI < nPoints; pI++)
 647     {
 648         if (pI < nPoints - 1)
 649         {
 650             nextOldPoint = oldPoints[operator[](pI + 1)];
 651             nextNewPoint = newPoints[operator[](pI + 1)];
 652         }
 653         else
 654         {
 655             nextOldPoint = oldPoints[operator[](0)];
 656             nextNewPoint = newPoints[operator[](0)];
 657         }
 658
 659         // Note: for best accuracy, centre point always comes last
 660         sv += triPointRef
 661               (
 662                   centreOldPoint,
 663                   oldPoints[operator[](pI)],
 664                   nextOldPoint
 665               ).sweptVol
 666               (
 667                   triPointRef
 668                   (
 669                       centreNewPoint,
 670                       newPoints[operator[](pI)],
 671                       nextNewPoint
 672                   )
 673               );
 674     }
 675
 676     return sv;
 677 }
 678
 679
 680 Foam::edgeList Foam::face::edges() const
 681 {
 682     const labelList& points = *this;
 683
 684     edgeList e(points.size());
 685
 686     label pointI;
 687
 688     for (pointI = 0; pointI < points.size() - 1; pointI++)
 689     {
 690         e[pointI] = edge(points[pointI], points[pointI + 1]);
 691     }
 692
 693     // add last edge
 694     e[points.size() - 1] = edge(points[points.size() - 1], points[0]);
 695
 696     return e;
 697 }
 698
 699
 700 int Foam::face::edgeDirection(const edge& e) const
 701 {
 702     if (size() > 2)
 703     {
 704         edge found(-1,-1);
 705
 706         // find start/end points - this breaks down for degenerate faces
 707         forAll (*this, i)
 708         {
 709             if (operator[](i) == e.start())
 710             {
 711                 found.start() = i;
 712             }
 713             else if (operator[](i) == e.end())
 714             {
 715                 found.end() = i;
 716             }
 717         }
 718
 719         label diff = found.end() - found.start();
 720         if (!diff || found.start() < 0 || found.end() < 0)
 721         {
 722             return 0;
 723         }
 724
 725         // forward direction
 726         if (diff == 1 || diff == 1 - size())
 727         {
 728             return 1;
 729         }
 730         // reverse direction
 731         if (diff ==  -1 || diff == -1 + size())
 732         {
 733             return -1;
 734         }
 735     }
 736
 737     return 0;
 738 }
 739
 740
 741 // Number of triangles directly known from number of vertices
 742 Foam::label Foam::face::nTriangles
 743 (
 744     const pointField&
 745 ) const
 746 {
 747     return size() - 2;
 748 }
 749
 750
 751 void Foam::face::triangles
 752 (
 753     const pointField& points,
 754     label& triI,
 755     faceList& triFaces
 756 ) const
 757 {
 758     faceList quadFaces;
 759     label quadI = 0;
 760
 761     split(SPLITTRIANGLE, points, triI, quadI, triFaces, quadFaces);
 762 }
 763
 764
 765 void Foam::face::nTrianglesQuads
 766 (
 767     const pointField& points,
 768     label& triI,
 769     label& quadI
 770 ) const
 771 {
 772     faceList triFaces;
 773     faceList quadFaces;
 774
 775     split(COUNTQUAD, points, triI, quadI, triFaces, quadFaces);
 776 }
 777
 778
 779 void Foam::face::trianglesQuads
 780 (
 781     const pointField& points,
 782     label& triI,
 783     label& quadI,
 784     faceList& triFaces,
 785     faceList& quadFaces
 786 ) const
 787 {
 788     split(SPLITQUAD, points, triI, quadI, triFaces, quadFaces);
 789 }
 790
 791
 792 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 793
 794
 795 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 796
 797 // ************************************************************************* //
 798