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