src/dynamicMesh/meshCut/cellLooper/topoCellLooper.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 "topoCellLooper.H"
  27 #include "cellFeatures.H"
  28 #include "polyMesh.H"
  29 #include "unitConversion.H"
  30 #include "DynamicList.H"
  31 #include "ListOps.H"
  32 #include "meshTools.H"
  33 #include "hexMatcher.H"
  34
  35 #include "addToRunTimeSelectionTable.H"
  36
  37 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  38
  39 namespace Foam
  40 {
  41    defineTypeNameAndDebug(topoCellLooper, 0);
  42    addToRunTimeSelectionTable(cellLooper, topoCellLooper, word);
  43 }
  44
  45 // Angle for polys to be considered splitHexes.
  46 const Foam::scalar Foam::topoCellLooper::featureCos = Foam::cos(degToRad(10.0));
  47
  48
  49 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  50
  51 // In-memory truncate a list
  52 template <class T>
  53 void Foam::topoCellLooper::subsetList
  54 (
  55     const label startI,
  56     const label freeI,
  57     DynamicList<T>& lst
  58 )
  59 {
  60     if (startI == 0)
  61     {
  62         // Truncate (setSize decides itself not to do anything if nothing
  63         // changed)
  64         if (freeI < 0)
  65         {
  66             FatalErrorIn("topoCellLooper::subsetList")
  67                 << "startI:" << startI << "  freeI:" << freeI
  68                 << "  lst:" << lst << abort(FatalError);
  69         }
  70         lst.setCapacity(freeI);
  71     }
  72     else
  73     {
  74         // Shift elements down
  75         label newI = 0;
  76         for (label elemI = startI; elemI < freeI; elemI++)
  77         {
  78             lst[newI++] = lst[elemI];
  79         }
  80
  81         if ((freeI - startI) < 0)
  82         {
  83             FatalErrorIn("topoCellLooper::subsetList")
  84                 << "startI:" << startI << "  freeI:" << freeI
  85                 << "  lst:" << lst << abort(FatalError);
  86         }
  87
  88         lst.setCapacity(freeI - startI);
  89     }
  90 }
  91
  92
  93 // Returns label of edge nFeaturePts away from startEdge (in the direction of
  94 // startVertI) and not counting non-featurePoints.
  95 //
  96 // When stepping to this face it can happen in 3 different ways:
  97 //
  98 //  --|------
  99 //    |
 100 //  1 |   0
 101 //    |A
 102 //    |
 103 //    |
 104 //  --|------
 105 //
 106 //  A: jumping from face0 to face1 across edge A.
 107 //      startEdge != -1
 108 //      startVert == -1
 109 //
 110 //  --|------
 111 //    |
 112 //  1 |   0
 113 //    +B
 114 //    |
 115 //    |
 116 //  --|------
 117 //
 118 //  B: jumping from face0 to face1 across (non-feature) point B
 119 //      startEdge == -1
 120 //      startVert != -1
 121 //
 122 //  --|------
 123 //  0 |   1
 124 //    |C
 125 //  --+
 126 //    |
 127 //    |
 128 //  --|------
 129 //
 130 //  C: jumping from face0 to face1 across (feature) point C on edge.
 131 //      startEdge != -1
 132 //      startVert != -1
 133 //
 134 void Foam::topoCellLooper::walkFace
 135 (
 136     const cellFeatures& features,
 137     const label faceI,
 138     const label startEdgeI,
 139     const label startVertI,
 140     const label nFeaturePts,
 141
 142     label& edgeI,
 143     label& vertI
 144 ) const
 145 {
 146     const labelList& fEdges = mesh().faceEdges()[faceI];
 147
 148     edgeI = startEdgeI;
 149
 150     vertI = startVertI;
 151
 152     // Number of feature points crossed so far
 153     label nVisited = 0;
 154
 155     if (vertI == -1)
 156     {
 157         // Started on edge. Go to one of its endpoints.
 158         vertI = mesh().edges()[edgeI].start();
 159
 160         if (features.isFeatureVertex(faceI, vertI))
 161         {
 162             nVisited++;
 163         }
 164     }
 165
 166     if ((edgeI == -1) || !meshTools::edgeOnFace(mesh(), faceI, edgeI))
 167     {
 168         // Either edge is not set or not on current face.  Just take one of
 169         // the edges on this face as starting edge.
 170         edgeI = getFirstVertEdge(faceI, vertI);
 171     }
 172
 173     // Now we should have starting edge on face and a vertex on that edge.
 174
 175     do
 176     {
 177
 178         edgeI = meshTools::otherEdge(mesh(), fEdges, edgeI, vertI);
 179
 180         if (nVisited == nFeaturePts)
 181         {
 182             break;
 183         }
 184
 185         vertI = mesh().edges()[edgeI].otherVertex(vertI);
 186
 187         if (features.isFeatureVertex(faceI, vertI))
 188         {
 189             nVisited++;
 190         }
 191     }
 192     while (true);
 193 }
 194
 195
 196 // Returns list of vertices on 'superEdge' i.e. list of edges connected by
 197 // non-feature points. First and last are feature points, ones inbetween are
 198 // not.
 199 Foam::labelList Foam::topoCellLooper::getSuperEdge
 200 (
 201     const cellFeatures& features,
 202     const label faceI,
 203     const label startEdgeI,
 204     const label startVertI
 205 ) const
 206 {
 207     const labelList& fEdges = mesh().faceEdges()[faceI];
 208
 209     labelList superVerts(fEdges.size());
 210     label superVertI = 0;
 211
 212
 213     label edgeI = startEdgeI;
 214
 215     label vertI = startVertI;
 216
 217     superVerts[superVertI++] = vertI;
 218
 219     label prevEdgeI = -1;
 220
 221     do
 222     {
 223         vertI = mesh().edges()[edgeI].otherVertex(vertI);
 224
 225         superVerts[superVertI++] = vertI;
 226
 227         prevEdgeI = edgeI;
 228
 229         edgeI = meshTools::otherEdge(mesh(), fEdges, edgeI, vertI);
 230     }
 231     while (!features.isFeaturePoint(prevEdgeI, edgeI));
 232
 233     superVerts.setSize(superVertI);
 234
 235     return superVerts;
 236 }
 237
 238
 239 // Return non-feature edge from cells' edges that is most perpendicular
 240 // to refinement direction.
 241 Foam::label Foam::topoCellLooper::getAlignedNonFeatureEdge
 242 (
 243     const vector& refDir,
 244     const label cellI,
 245     const cellFeatures& features
 246 ) const
 247 {
 248     const labelList& cEdges = mesh().cellEdges()[cellI];
 249
 250     const point& ctr = mesh().cellCentres()[cellI];
 251
 252     label cutEdgeI = -1;
 253     scalar maxCos = -GREAT;
 254
 255     forAll(cEdges, cEdgeI)
 256     {
 257         label edgeI = cEdges[cEdgeI];
 258
 259         if (!features.isFeatureEdge(edgeI))
 260         {
 261             const edge& e = mesh().edges()[edgeI];
 262
 263             // Get plane spanned by e.start, e.end and cell centre.
 264             vector e0 = mesh().points()[e.start()] - ctr;
 265             vector e1 = mesh().points()[e.end()] - ctr;
 266
 267             vector n = e0 ^ e1;
 268             n /= mag(n);
 269
 270             scalar cosAngle = mag(refDir & n);
 271
 272             if (cosAngle > maxCos)
 273             {
 274                 maxCos = cosAngle;
 275
 276                 cutEdgeI = edgeI;
 277             }
 278         }
 279     }
 280
 281     return cutEdgeI;
 282 }
 283
 284
 285 // Starts from edge and vertex on edge on face (or neighbouring face)
 286 // and steps either to existing vertex (vertI != -1) or to edge (vertI == -1)
 287 // by walking point-edge and crossing nFeats featurePoints.
 288 void Foam::topoCellLooper::walkAcrossFace
 289 (
 290     const cellFeatures& features,
 291     const label faceI,
 292     const label startEdgeI,
 293     const label startVertI,
 294     const label nFeats,
 295
 296     label& edgeI,
 297     label& vertI
 298 ) const
 299 {
 300     label oppositeVertI = -1;
 301     label oppositeEdgeI = -1;
 302
 303     // Go to oppositeEdge and a vertex on that.
 304     walkFace
 305     (
 306         features,
 307         faceI,
 308         startEdgeI,
 309         startVertI,
 310         nFeats,
 311
 312         oppositeEdgeI,
 313         oppositeVertI
 314     );
 315
 316     // Loop over super edge to find internal points if there are any.
 317
 318     labelList superEdge =
 319         getSuperEdge
 320         (
 321             features,
 322             faceI,
 323             oppositeEdgeI,
 324             oppositeVertI
 325         );
 326
 327     label sz = superEdge.size();
 328
 329     if (sz == 2)
 330     {
 331         // No non-feature point inbetween feature points.
 332         // Mark edge.
 333
 334         vertI = -1;
 335         edgeI = oppositeEdgeI;
 336     }
 337     else if (sz == 3)
 338     {
 339         vertI = superEdge[1];
 340         edgeI = -1;
 341     }
 342     else
 343     {
 344         //Should choose acc. to geometry!
 345         label index = sz/2;
 346
 347         if (debug)
 348         {
 349             Pout<< "    Don't know what to do. Stepped to non-feature point "
 350                 << "at index " << index << " in superEdge:" << superEdge
 351                 << endl;
 352         }
 353
 354         vertI = superEdge[index];
 355         edgeI = -1;
 356     }
 357 }
 358
 359
 360 // Walks cell circumference. Updates face, edge, vertex.
 361 //
 362 // Position on face is given by:
 363 //
 364 //  vertI == -1, faceI != -1, edgeI != -1
 365 //      on edge of face. Cross edge to neighbouring face.
 366 //
 367 //  vertI != -1, edgeI != -1, faceI == -1
 368 //      coming from edge onto vertex vertI. Need to step to one
 369 //      of the faces not using edgeI.
 370 //
 371 //  vertI != -1, edgeI == -1, faceI != -1
 372 //      coming from vertex on side of face. Step to one of the faces
 373 //      using vertI but not faceI
 374 void Foam::topoCellLooper::walkSplitHex
 375 (
 376     const label cellI,
 377     const cellFeatures& features,
 378     const label fromFaceI,
 379     const label fromEdgeI,
 380     const label fromVertI,
 381
 382     DynamicList<label>& loop,
 383     DynamicList<scalar>& loopWeights
 384 ) const
 385 {
 386     // Work vars giving position on cell
 387     label faceI = fromFaceI;
 388     label edgeI = fromEdgeI;
 389     label vertI = fromVertI;
 390
 391     do
 392     {
 393         if (debug)
 394         {
 395             Pout<< "Entering walk with : cell:" << cellI << " face:" << faceI;
 396             if (faceI != -1)
 397             {
 398                 Pout<< " verts:" << mesh().faces()[faceI];
 399             }
 400             Pout<< " edge:" << edgeI;
 401             if (edgeI != -1)
 402             {
 403                 Pout<< " verts:" << mesh().edges()[edgeI];
 404             }
 405             Pout<< " vert:" << vertI << endl;
 406         }
 407
 408         label startLoop = -1;
 409
 410         if
 411         (
 412             (vertI != -1)
 413          && (
 414                 (startLoop =
 415                     findIndex
 416                     (
 417                         loop,
 418                         vertToEVert(vertI)
 419                     )
 420                 )
 421             != -1
 422             )
 423         )
 424         {
 425             // Breaking walk since vertI already cut
 426             label firstFree = loop.size();
 427
 428             subsetList(startLoop, firstFree, loop);
 429             subsetList(startLoop, firstFree, loopWeights);
 430
 431             break;
 432         }
 433         if
 434         (
 435             (edgeI != -1)
 436          && (
 437                 (startLoop =
 438                     findIndex
 439                     (
 440                         loop,
 441                         edgeToEVert(edgeI)
 442                     )
 443                 )
 444              != -1
 445             )
 446         )
 447         {
 448             // Breaking walk since edgeI already cut
 449             label firstFree = loop.size();
 450
 451             subsetList(startLoop, firstFree, loop);
 452             subsetList(startLoop, firstFree, loopWeights);
 453
 454             break;
 455         }
 456
 457
 458         if (vertI == -1)
 459         {
 460             // On edge
 461             if (edgeI == -1)
 462             {
 463                 FatalErrorIn("walkSplitHex") << "Illegal edge and vert"
 464                     << abort(FatalError);
 465             }
 466
 467             loop.append(edgeToEVert(edgeI));
 468             loopWeights.append(0.5);
 469
 470             // Cross edge to next face
 471             faceI = meshTools::otherFace(mesh(), cellI, faceI, edgeI);
 472
 473             if (debug)
 474             {
 475                 Pout<< "    stepped across edge " << mesh().edges()[edgeI]
 476                     << " to face " << faceI << " verts:"
 477                     << mesh().faces()[faceI] << endl;
 478             }
 479
 480             label nextEdgeI = -1;
 481             label nextVertI = -1;
 482
 483             // Walk face along its edges
 484             walkAcrossFace
 485             (
 486                 features,
 487                 faceI,
 488                 edgeI,
 489                 vertI,
 490                 2,
 491
 492                 nextEdgeI,
 493                 nextVertI
 494             );
 495
 496             edgeI = nextEdgeI;
 497             vertI = nextVertI;
 498         }
 499         else
 500         {
 501             // On vertex.
 502
 503             loop.append(vertToEVert(vertI));
 504             loopWeights.append(-GREAT);
 505
 506             if (edgeI == -1)
 507             {
 508                 // Normal vertex on edge of face. Get edges connected to it
 509                 // which are not on faceI.
 510                 labelList nextEdges = getVertEdgesNonFace
 511                 (
 512                     cellI,
 513                     faceI,
 514                     vertI
 515                 );
 516
 517                 if (nextEdges.empty())
 518                 {
 519                     // Cross to other face (there is only one since no edges)
 520                     const labelList& pFaces = mesh().pointFaces()[vertI];
 521
 522                     forAll(pFaces, pFaceI)
 523                     {
 524                         label thisFaceI = pFaces[pFaceI];
 525
 526                         if
 527                         (
 528                             (thisFaceI != faceI)
 529                          && meshTools::faceOnCell(mesh(), cellI, thisFaceI)
 530                         )
 531                         {
 532                             faceI = thisFaceI;
 533                             break;
 534                         }
 535                     }
 536
 537                     if (debug)
 538                     {
 539                         Pout<< "    stepped from non-edge vertex " << vertI
 540                             << " to face " << faceI << " verts:"
 541                             << mesh().faces()[faceI]
 542                             << " since candidate edges:" << nextEdges << endl;
 543                     }
 544
 545                     label nextEdgeI = -1;
 546                     label nextVertI = -1;
 547
 548                     walkAcrossFace
 549                     (
 550                         features,
 551                         faceI,
 552                         edgeI,
 553                         vertI,
 554                         2,          // 2 vertices to cross
 555
 556                         nextEdgeI,
 557                         nextVertI
 558                     );
 559
 560                     edgeI = nextEdgeI;
 561                     vertI = nextVertI;
 562                 }
 563                 else if (nextEdges.size() == 1)
 564                 {
 565                     // One edge. Go along this one.
 566                     edgeI = nextEdges[0];
 567
 568                     if (debug)
 569                     {
 570                         Pout<< "    stepped from non-edge vertex " << vertI
 571                             << " along edge " << edgeI << " verts:"
 572                             << mesh().edges()[edgeI]
 573                             << " out of candidate edges:"
 574                             << nextEdges << endl;
 575                     }
 576
 577                     vertI = mesh().edges()[edgeI].otherVertex(vertI);
 578
 579                     faceI = -1;
 580                 }
 581                 else
 582                 {
 583                     // Multiple edges to choose from. Pick any one.
 584                     // (ideally should be geometric)
 585
 586                     label index = nextEdges.size()/2;
 587
 588                     edgeI = nextEdges[index];
 589
 590                     if (debug)
 591                     {
 592                         Pout<< "    stepped from non-edge vertex " << vertI
 593                             << " along edge " << edgeI << " verts:"
 594                             << mesh().edges()[edgeI]
 595                             << " out of candidate edges:" << nextEdges << endl;
 596                     }
 597
 598                     vertI = mesh().edges()[edgeI].otherVertex(vertI);
 599
 600                     faceI = -1;
 601                 }
 602             }
 603             else
 604             {
 605                 // Get faces connected to startVertI but not startEdgeI
 606                 labelList nextFaces =
 607                     getVertFacesNonEdge
 608                     (
 609                         cellI,
 610                         edgeI,
 611                         vertI
 612                     );
 613
 614                 if (nextFaces.size() == 1)
 615                 {
 616                     // Only one face to cross.
 617                     faceI = nextFaces[0];
 618
 619                     label nextEdgeI = -1;
 620                     label nextVertI = -1;
 621
 622                     walkAcrossFace
 623                     (
 624                         features,
 625                         faceI,
 626                         edgeI,
 627                         vertI,
 628                         2,          // 2 vertices to cross
 629
 630                         nextEdgeI,
 631                         nextVertI
 632                     );
 633
 634                     edgeI = nextEdgeI;
 635                     vertI = nextVertI;
 636                 }
 637                 else if (nextFaces.size() == 2)
 638                 {
 639                     // Split face. Get edge inbetween.
 640                     faceI = -1;
 641
 642                     edgeI =
 643                         meshTools::getSharedEdge
 644                         (
 645                             mesh(),
 646                             nextFaces[0],
 647                             nextFaces[1]
 648                         );
 649
 650                     vertI = mesh().edges()[edgeI].otherVertex(vertI);
 651                 }
 652                 else
 653                 {
 654                     FatalErrorIn("walkFromVert") << "Not yet implemented"
 655                         << "Choosing from more than "
 656                         << "two candidates:" << nextFaces
 657                         << " when coming from vertex " << vertI << " on cell "
 658                         << cellI << abort(FatalError);
 659                 }
 660             }
 661         }
 662
 663         if (debug)
 664         {
 665             Pout<< "Walked to : face:" << faceI;
 666             if (faceI != -1)
 667             {
 668                 Pout<< " verts:" << mesh().faces()[faceI];
 669             }
 670             Pout<< " edge:" << edgeI;
 671             if (edgeI != -1)
 672             {
 673                 Pout<< " verts:" << mesh().edges()[edgeI];
 674             }
 675             Pout<< " vert:" << vertI << endl;
 676         }
 677     }
 678     while (true);
 679 }
 680
 681
 682 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 683
 684 // Construct from components
 685 Foam::topoCellLooper::topoCellLooper(const polyMesh& mesh)
 686 :
 687     hexCellLooper(mesh)
 688 {}
 689
 690
 691 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 692
 693 Foam::topoCellLooper::~topoCellLooper()
 694 {}
 695
 696
 697 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 698
 699 bool Foam::topoCellLooper::cut
 700 (
 701     const vector& refDir,
 702     const label cellI,
 703     const boolList& vertIsCut,
 704     const boolList& edgeIsCut,
 705     const scalarField& edgeWeight,
 706
 707     labelList& loop,
 708     scalarField& loopWeights
 709 ) const
 710 {
 711     if (mesh().cellShapes()[cellI].model() == hex_)
 712     {
 713         // Let parent handle hex case.
 714         return
 715             hexCellLooper::cut
 716             (
 717                 refDir,
 718                 cellI,
 719                 vertIsCut,
 720                 edgeIsCut,
 721                 edgeWeight,
 722                 loop,
 723                 loopWeights
 724             );
 725     }
 726     else
 727     {
 728         cellFeatures superCell(mesh(), featureCos, cellI);
 729
 730         if (hexMatcher().isA(superCell.faces()))
 731         {
 732             label edgeI =
 733                 getAlignedNonFeatureEdge
 734                 (
 735                     refDir,
 736                     cellI,
 737                     superCell
 738                 );
 739
 740             label vertI = -1;
 741
 742             label faceI = -1;
 743
 744             if (edgeI != -1)
 745             {
 746                 // Found non-feature edge. Start walking from vertex on edge.
 747                 vertI = mesh().edges()[edgeI].start();
 748             }
 749             else
 750             {
 751                 // No 'matching' non-feature edge found on cell. Get starting
 752                 // normal i.e. feature edge.
 753                 edgeI = getMisAlignedEdge(refDir, cellI);
 754
 755                 // Get any face using edge
 756                 label face0;
 757                 label face1;
 758                 meshTools::getEdgeFaces(mesh(), cellI, edgeI, face0, face1);
 759
 760                 faceI = face0;
 761             }
 762
 763             label nEstCuts = 2*mesh().cells()[cellI].size();
 764
 765             DynamicList<label> localLoop(nEstCuts);
 766             DynamicList<scalar> localLoopWeights(nEstCuts);
 767
 768             walkSplitHex
 769             (
 770                 cellI,
 771                 superCell,
 772                 faceI,
 773                 edgeI,
 774                 vertI,
 775
 776                 localLoop,
 777                 localLoopWeights
 778             );
 779
 780             if (localLoop.size() <=2)
 781             {
 782                 return false;
 783             }
 784             else
 785             {
 786                 loop.transfer(localLoop);
 787                 loopWeights.transfer(localLoopWeights);
 788
 789                 return true;
 790             }
 791         }
 792         else
 793         {
 794             // Let parent handle poly case.
 795             return hexCellLooper::cut
 796             (
 797                 refDir,
 798                 cellI,
 799                 vertIsCut,
 800                 edgeIsCut,
 801                 edgeWeight,
 802                 loop,
 803                 loopWeights
 804             );
 805         }
 806     }
 807 }
 808
 809
 810 bool Foam::topoCellLooper::cut
 811 (
 812     const plane& cutPlane,
 813     const label cellI,
 814     const boolList& vertIsCut,
 815     const boolList& edgeIsCut,
 816     const scalarField& edgeWeight,
 817
 818     labelList& loop,
 819     scalarField& loopWeights
 820 ) const
 821 {
 822     // Let parent handle cut with plane.
 823     return
 824         hexCellLooper::cut
 825         (
 826             cutPlane,
 827             cellI,
 828             vertIsCut,
 829             edgeIsCut,
 830             edgeWeight,
 831
 832             loop,
 833             loopWeights
 834         );
 835 }
 836
 837
 838 // ************************************************************************* //