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