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