src/dynamicMesh/meshCut/meshModifiers/meshCutAndRemove/meshCutAndRemove.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 "meshCutAndRemove.H"
  27 #include "polyMesh.H"
  28 #include "polyTopoChange.H"
  29 #include "polyAddFace.H"
  30 #include "polyAddPoint.H"
  31 #include "polyRemovePoint.H"
  32 #include "polyRemoveFace.H"
  33 #include "polyModifyFace.H"
  34 #include "cellCuts.H"
  35 #include "mapPolyMesh.H"
  36 #include "meshTools.H"
  37
  38 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  39
  40 defineTypeNameAndDebug(Foam::meshCutAndRemove, 0);
  41
  42 // * * * * * * * * * * * * * Private Static Functions  * * * * * * * * * * * //
  43
  44 // Returns -1 or index in elems1 of first shared element.
  45 Foam::label Foam::meshCutAndRemove::firstCommon
  46 (
  47     const labelList& elems1,
  48     const labelList& elems2
  49 )
  50 {
  51     forAll(elems1, elemI)
  52     {
  53         label index1 = findIndex(elems2, elems1[elemI]);
  54
  55         if (index1 != -1)
  56         {
  57             return index1;
  58         }
  59     }
  60     return -1;
  61 }
  62
  63
  64 // Check if twoCuts at two consecutive position in cuts.
  65 bool Foam::meshCutAndRemove::isIn
  66 (
  67     const edge& twoCuts,
  68     const labelList& cuts
  69 )
  70 {
  71     label index = findIndex(cuts, twoCuts[0]);
  72
  73     if (index == -1)
  74     {
  75         return false;
  76     }
  77
  78     return
  79     (
  80         cuts[cuts.fcIndex(index)] == twoCuts[1]
  81      || cuts[cuts.rcIndex(index)] == twoCuts[1]
  82     );
  83 }
  84
  85
  86 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  87
  88 // Returns the cell in cellLabels that is cut. Or -1.
  89 Foam::label Foam::meshCutAndRemove::findCutCell
  90 (
  91     const cellCuts& cuts,
  92     const labelList& cellLabels
  93 ) const
  94 {
  95     forAll(cellLabels, labelI)
  96     {
  97         label cellI = cellLabels[labelI];
  98
  99         if (cuts.cellLoops()[cellI].size())
 100         {
 101             return cellI;
 102         }
 103     }
 104     return -1;
 105 }
 106
 107
 108 //- Returns first pointI in pointLabels that uses an internal
 109 //  face. Used to find point to inflate cell/face from (has to be
 110 //  connected to internal face). Returns -1 (so inflate from nothing) if
 111 //  none found.
 112 Foam::label Foam::meshCutAndRemove::findInternalFacePoint
 113 (
 114     const labelList& pointLabels
 115 ) const
 116 {
 117     forAll(pointLabels, labelI)
 118     {
 119         label pointI = pointLabels[labelI];
 120
 121         const labelList& pFaces = mesh().pointFaces()[pointI];
 122
 123         forAll(pFaces, pFaceI)
 124         {
 125             label faceI = pFaces[pFaceI];
 126
 127             if (mesh().isInternalFace(faceI))
 128             {
 129                 return pointI;
 130             }
 131         }
 132     }
 133
 134     if (pointLabels.empty())
 135     {
 136         FatalErrorIn
 137         (
 138             "meshCutAndRemove::findInternalFacePoint(const labelList&)"
 139         )
 140             << "Empty pointLabels" << abort(FatalError);
 141     }
 142
 143     return -1;
 144 }
 145
 146
 147 // Find point on face that is part of original mesh and that is point connected
 148 // to the patch
 149 Foam::label Foam::meshCutAndRemove::findPatchFacePoint
 150 (
 151     const face& f,
 152     const label exposedPatchI
 153 ) const
 154 {
 155     const labelListList& pointFaces = mesh().pointFaces();
 156     const polyBoundaryMesh& patches = mesh().boundaryMesh();
 157
 158     forAll(f, fp)
 159     {
 160         label pointI = f[fp];
 161
 162         if (pointI < mesh().nPoints())
 163         {
 164             const labelList& pFaces = pointFaces[pointI];
 165
 166             forAll(pFaces, i)
 167             {
 168                 if (patches.whichPatch(pFaces[i]) == exposedPatchI)
 169                 {
 170                     return pointI;
 171                 }
 172             }
 173         }
 174     }
 175     return -1;
 176 }
 177
 178
 179 // Get new owner and neighbour of face. Checks anchor points to see if
 180 // cells have been removed.
 181 void Foam::meshCutAndRemove::faceCells
 182 (
 183     const cellCuts& cuts,
 184     const label exposedPatchI,
 185     const label faceI,
 186     label& own,
 187     label& nei,
 188     label& patchID
 189 ) const
 190 {
 191     const labelListList& anchorPts = cuts.cellAnchorPoints();
 192     const labelListList& cellLoops = cuts.cellLoops();
 193
 194     const face& f = mesh().faces()[faceI];
 195
 196     own = mesh().faceOwner()[faceI];
 197
 198     if (cellLoops[own].size() && firstCommon(f, anchorPts[own]) == -1)
 199     {
 200         // owner has been split and this is the removed part.
 201         own = -1;
 202     }
 203
 204     nei = -1;
 205
 206     if (mesh().isInternalFace(faceI))
 207     {
 208         nei = mesh().faceNeighbour()[faceI];
 209
 210         if (cellLoops[nei].size() && firstCommon(f, anchorPts[nei]) == -1)
 211         {
 212             nei = -1;
 213         }
 214     }
 215
 216     patchID = mesh().boundaryMesh().whichPatch(faceI);
 217
 218     if (patchID == -1 && (own == -1 || nei == -1))
 219     {
 220         // Face was internal but becomes external
 221         patchID = exposedPatchI;
 222     }
 223 }
 224
 225
 226 void Foam::meshCutAndRemove::getZoneInfo
 227 (
 228     const label faceI,
 229     label& zoneID,
 230     bool& zoneFlip
 231 ) const
 232 {
 233     zoneID = mesh().faceZones().whichZone(faceI);
 234
 235     zoneFlip = false;
 236
 237     if (zoneID >= 0)
 238     {
 239         const faceZone& fZone = mesh().faceZones()[zoneID];
 240
 241         zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
 242     }
 243 }
 244
 245
 246 // Adds a face from point.
 247 void Foam::meshCutAndRemove::addFace
 248 (
 249     polyTopoChange& meshMod,
 250     const label faceI,
 251     const label masterPointI,
 252     const face& newFace,
 253     const label own,
 254     const label nei,
 255     const label patchID
 256 )
 257 {
 258     label zoneID;
 259     bool zoneFlip;
 260
 261     getZoneInfo(faceI, zoneID, zoneFlip);
 262
 263     if ((nei == -1) || (own != -1 && own < nei))
 264     {
 265         // Ordering ok.
 266         if (debug & 2)
 267         {
 268             Pout<< "Adding face " << newFace
 269                 << " with new owner:" << own
 270                 << " with new neighbour:" << nei
 271                 << " patchID:" << patchID
 272                 << " anchor:" << masterPointI
 273                 << " zoneID:" << zoneID
 274                 << " zoneFlip:" << zoneFlip
 275                 << endl;
 276         }
 277
 278         meshMod.setAction
 279         (
 280             polyAddFace
 281             (
 282                 newFace,                    // face
 283                 own,                        // owner
 284                 nei,                        // neighbour
 285                 masterPointI,               // master point
 286                 -1,                         // master edge
 287                 -1,                         // master face for addition
 288                 false,                      // flux flip
 289                 patchID,                    // patch for face
 290                 zoneID,                     // zone for face
 291                 zoneFlip                    // face zone flip
 292             )
 293         );
 294     }
 295     else
 296     {
 297         // Reverse owner/neighbour
 298         if (debug & 2)
 299         {
 300             Pout<< "Adding (reversed) face " << newFace.reverseFace()
 301                 << " with new owner:" << nei
 302                 << " with new neighbour:" << own
 303                 << " patchID:" << patchID
 304                 << " anchor:" << masterPointI
 305                 << " zoneID:" << zoneID
 306                 << " zoneFlip:" << zoneFlip
 307                 << endl;
 308         }
 309
 310         meshMod.setAction
 311         (
 312             polyAddFace
 313             (
 314                 newFace.reverseFace(),      // face
 315                 nei,                        // owner
 316                 own,                        // neighbour
 317                 masterPointI,               // master point
 318                 -1,                         // master edge
 319                 -1,                         // master face for addition
 320                 false,                      // flux flip
 321                 patchID,                    // patch for face
 322                 zoneID,                     // zone for face
 323                 zoneFlip                    // face zone flip
 324             )
 325         );
 326     }
 327 }
 328
 329
 330 // Modifies existing faceI for either new owner/neighbour or new face points.
 331 void Foam::meshCutAndRemove::modFace
 332 (
 333     polyTopoChange& meshMod,
 334     const label faceI,
 335     const face& newFace,
 336     const label own,
 337     const label nei,
 338     const label patchID
 339 )
 340 {
 341     label zoneID;
 342     bool zoneFlip;
 343
 344     getZoneInfo(faceI, zoneID, zoneFlip);
 345
 346     if
 347     (
 348         (own != mesh().faceOwner()[faceI])
 349      || (
 350             mesh().isInternalFace(faceI)
 351          && (nei != mesh().faceNeighbour()[faceI])
 352         )
 353      || (newFace != mesh().faces()[faceI])
 354     )
 355     {
 356         if (debug & 2)
 357         {
 358             Pout<< "Modifying face " << faceI
 359                 << " old vertices:" << mesh().faces()[faceI]
 360                 << " new vertices:" << newFace
 361                 << " new owner:" << own
 362                 << " new neighbour:" << nei
 363                 << " new patch:" << patchID
 364                 << " new zoneID:" << zoneID
 365                 << " new zoneFlip:" << zoneFlip
 366                 << endl;
 367         }
 368
 369         if ((nei == -1) || (own != -1 && own < nei))
 370         {
 371             meshMod.setAction
 372             (
 373                 polyModifyFace
 374                 (
 375                     newFace,            // modified face
 376                     faceI,              // label of face being modified
 377                     own,                // owner
 378                     nei,                // neighbour
 379                     false,              // face flip
 380                     patchID,            // patch for face
 381                     false,              // remove from zone
 382                     zoneID,             // zone for face
 383                     zoneFlip            // face flip in zone
 384                 )
 385             );
 386         }
 387         else
 388         {
 389             meshMod.setAction
 390             (
 391                 polyModifyFace
 392                 (
 393                     newFace.reverseFace(),  // modified face
 394                     faceI,                  // label of face being modified
 395                     nei,                    // owner
 396                     own,                    // neighbour
 397                     false,                  // face flip
 398                     patchID,                // patch for face
 399                     false,                  // remove from zone
 400                     zoneID,                 // zone for face
 401                     zoneFlip                // face flip in zone
 402                 )
 403             );
 404         }
 405     }
 406 }
 407
 408
 409 // Copies face starting from startFp up to and including endFp.
 410 void Foam::meshCutAndRemove::copyFace
 411 (
 412     const face& f,
 413     const label startFp,
 414     const label endFp,
 415     face& newFace
 416 ) const
 417 {
 418     label fp = startFp;
 419
 420     label newFp = 0;
 421
 422     while (fp != endFp)
 423     {
 424         newFace[newFp++] = f[fp];
 425
 426         fp = (fp + 1) % f.size();
 427     }
 428     newFace[newFp] = f[fp];
 429 }
 430
 431
 432 // Actually split face in two along splitEdge v0, v1 (the two vertices in new
 433 // vertex numbering). Generates faces in same ordering
 434 // as original face. Replaces cutEdges by the points introduced on them
 435 // (addedPoints_).
 436 void Foam::meshCutAndRemove::splitFace
 437 (
 438     const face& f,
 439     const label v0,
 440     const label v1,
 441
 442     face& f0,
 443     face& f1
 444 ) const
 445 {
 446     // Check if we find any new vertex which is part of the splitEdge.
 447     label startFp = findIndex(f, v0);
 448
 449     if (startFp == -1)
 450     {
 451         FatalErrorIn
 452         (
 453             "meshCutAndRemove::splitFace"
 454             ", const face&, const label, const label, face&, face&)"
 455         )   << "Cannot find vertex (new numbering) " << v0
 456             << " on face " << f
 457             << abort(FatalError);
 458     }
 459
 460     label endFp = findIndex(f, v1);
 461
 462     if (endFp == -1)
 463     {
 464         FatalErrorIn
 465         (
 466             "meshCutAndRemove::splitFace("
 467             ", const face&, const label, const label, face&, face&)"
 468         )   << "Cannot find vertex (new numbering) " << v1
 469             << " on face " << f
 470             << abort(FatalError);
 471     }
 472
 473
 474     f0.setSize((endFp + 1 + f.size() - startFp) % f.size());
 475     f1.setSize(f.size() - f0.size() + 2);
 476
 477     copyFace(f, startFp, endFp, f0);
 478     copyFace(f, endFp, startFp, f1);
 479 }
 480
 481
 482 // Adds additional vertices (from edge cutting) to face. Used for faces which
 483 // are not split but still might use edge that has been cut.
 484 Foam::face Foam::meshCutAndRemove::addEdgeCutsToFace(const label faceI) const
 485 {
 486     const face& f = mesh().faces()[faceI];
 487
 488     face newFace(2 * f.size());
 489
 490     label newFp = 0;
 491
 492     forAll(f, fp)
 493     {
 494         // Duplicate face vertex.
 495         newFace[newFp++] = f[fp];
 496
 497         // Check if edge has been cut.
 498         label fp1 = f.fcIndex(fp);
 499
 500         HashTable<label, edge, Hash<edge> >::const_iterator fnd =
 501             addedPoints_.find(edge(f[fp], f[fp1]));
 502
 503         if (fnd != addedPoints_.end())
 504         {
 505             // edge has been cut. Introduce new vertex.
 506             newFace[newFp++] = fnd();
 507         }
 508     }
 509
 510     newFace.setSize(newFp);
 511
 512     return newFace;
 513 }
 514
 515
 516 // Walk loop (loop of cuts) across circumference of cellI. Returns face in
 517 // new vertices.
 518 // Note: tricky bit is that it can use existing edges which have been split.
 519 Foam::face Foam::meshCutAndRemove::loopToFace
 520 (
 521     const label cellI,
 522     const labelList& loop
 523 ) const
 524 {
 525     face newFace(2*loop.size());
 526
 527     label newFaceI = 0;
 528
 529     forAll(loop, fp)
 530     {
 531         label cut = loop[fp];
 532
 533         if (isEdge(cut))
 534         {
 535             label edgeI = getEdge(cut);
 536
 537             const edge& e = mesh().edges()[edgeI];
 538
 539             label vertI = addedPoints_[e];
 540
 541             newFace[newFaceI++] = vertI;
 542         }
 543         else
 544         {
 545             // cut is vertex.
 546             label vertI = getVertex(cut);
 547
 548             newFace[newFaceI++] = vertI;
 549
 550             label nextCut = loop[loop.fcIndex(fp)];
 551
 552             if (!isEdge(nextCut))
 553             {
 554                 // From vertex to vertex -> cross cut only if no existing edge.
 555                 label nextVertI = getVertex(nextCut);
 556
 557                 label edgeI = meshTools::findEdge(mesh(), vertI, nextVertI);
 558
 559                 if (edgeI != -1)
 560                 {
 561                     // Existing edge. Insert split-edge point if any.
 562                     HashTable<label, edge, Hash<edge> >::const_iterator fnd =
 563                         addedPoints_.find(mesh().edges()[edgeI]);
 564
 565                     if (fnd != addedPoints_.end())
 566                     {
 567                         newFace[newFaceI++] = fnd();
 568                     }
 569                 }
 570             }
 571         }
 572     }
 573     newFace.setSize(newFaceI);
 574
 575     return newFace;
 576 }
 577
 578
 579 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 580
 581 // Construct from components
 582 Foam::meshCutAndRemove::meshCutAndRemove(const polyMesh& mesh)
 583 :
 584     edgeVertex(mesh),
 585     addedFaces_(),
 586     addedPoints_()
 587 {}
 588
 589
 590 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 591
 592 void Foam::meshCutAndRemove::setRefinement
 593 (
 594     const label exposedPatchI,
 595     const cellCuts& cuts,
 596     const labelList& cutPatch,
 597     polyTopoChange& meshMod
 598 )
 599 {
 600     // Clear and size maps here since mesh size will change.
 601     addedFaces_.clear();
 602     addedFaces_.resize(cuts.nLoops());
 603
 604     addedPoints_.clear();
 605     addedPoints_.resize(cuts.nLoops());
 606
 607     if (cuts.nLoops() == 0)
 608     {
 609         return;
 610     }
 611
 612     const labelListList& anchorPts = cuts.cellAnchorPoints();
 613     const labelListList& cellLoops = cuts.cellLoops();
 614     const polyBoundaryMesh& patches = mesh().boundaryMesh();
 615
 616     if (exposedPatchI < 0 || exposedPatchI >= patches.size())
 617     {
 618         FatalErrorIn
 619         (
 620             "meshCutAndRemove::setRefinement("
 621             ", const label, const cellCuts&, const labelList&"
 622             ", polyTopoChange&)"
 623         )   << "Illegal exposed patch " << exposedPatchI
 624             << abort(FatalError);
 625     }
 626
 627
 628     //
 629     // Add new points along cut edges.
 630     //
 631
 632     forAll(cuts.edgeIsCut(), edgeI)
 633     {
 634         if (cuts.edgeIsCut()[edgeI])
 635         {
 636             const edge& e = mesh().edges()[edgeI];
 637
 638             // Check if there is any cell using this edge.
 639             if (debug && findCutCell(cuts, mesh().edgeCells()[edgeI]) == -1)
 640             {
 641                 FatalErrorIn
 642                 (
 643                     "meshCutAndRemove::setRefinement("
 644                     ", const label, const cellCuts&, const labelList&"
 645                     ", polyTopoChange&)"
 646                 )   << "Problem: cut edge but none of the cells using it is\n"
 647                     << "edge:" << edgeI << " verts:" << e
 648                     << abort(FatalError);
 649             }
 650
 651             // One of the edge end points should be master point of nbCellI.
 652             label masterPointI = e.start();
 653
 654             const point& v0 = mesh().points()[e.start()];
 655             const point& v1 = mesh().points()[e.end()];
 656
 657             scalar weight = cuts.edgeWeight()[edgeI];
 658
 659             point newPt = weight*v1 + (1.0-weight)*v0;
 660
 661             label addedPointI =
 662                 meshMod.setAction
 663                 (
 664                     polyAddPoint
 665                     (
 666                         newPt,              // point
 667                         masterPointI,       // master point
 668                         -1,                 // zone for point
 669                         true                // supports a cell
 670                     )
 671                 );
 672
 673             // Store on (hash of) edge.
 674             addedPoints_.insert(e, addedPointI);
 675
 676             if (debug & 2)
 677             {
 678                 Pout<< "Added point " << addedPointI
 679                     << " to vertex "
 680                     << masterPointI << " of edge " << edgeI
 681                     << " vertices " << e << endl;
 682             }
 683         }
 684     }
 685
 686
 687     //
 688     // Remove all points that will not be used anymore
 689     //
 690     {
 691         boolList usedPoint(mesh().nPoints(), false);
 692
 693         forAll(cellLoops, cellI)
 694         {
 695             const labelList& loop = cellLoops[cellI];
 696
 697             if (loop.size())
 698             {
 699                 // Cell is cut. Uses only anchor points and loop itself.
 700                 forAll(loop, fp)
 701                 {
 702                     label cut = loop[fp];
 703
 704                     if (!isEdge(cut))
 705                     {
 706                         usedPoint[getVertex(cut)] = true;
 707                     }
 708                 }
 709
 710                 const labelList& anchors = anchorPts[cellI];
 711
 712                 forAll(anchors, i)
 713                 {
 714                     usedPoint[anchors[i]] = true;
 715                 }
 716             }
 717             else
 718             {
 719                 // Cell is not cut so use all its points
 720                 const labelList& cPoints = mesh().cellPoints()[cellI];
 721
 722                 forAll(cPoints, i)
 723                 {
 724                     usedPoint[cPoints[i]] = true;
 725                 }
 726             }
 727         }
 728
 729
 730         // Check
 731         const Map<edge>& faceSplitCut = cuts.faceSplitCut();
 732
 733         forAllConstIter(Map<edge>, faceSplitCut, iter)
 734         {
 735             const edge& fCut = iter();
 736
 737             forAll(fCut, i)
 738             {
 739                 label cut = fCut[i];
 740
 741                 if (!isEdge(cut))
 742                 {
 743                     label pointI = getVertex(cut);
 744
 745                     if (!usedPoint[pointI])
 746                     {
 747                         FatalErrorIn
 748                         (
 749                             "meshCutAndRemove::setRefinement("
 750                             ", const label, const cellCuts&, const labelList&"
 751                             ", polyTopoChange&)"
 752                         )   << "Problem: faceSplitCut not used by any loop"
 753                             << " or cell anchor point"
 754                             << "face:" << iter.key() << " point:" << pointI
 755                             << " coord:" << mesh().points()[pointI]
 756                             << abort(FatalError);
 757                     }
 758                 }
 759             }
 760         }
 761
 762         forAll(cuts.pointIsCut(), pointI)
 763         {
 764             if (cuts.pointIsCut()[pointI])
 765             {
 766                 if (!usedPoint[pointI])
 767                 {
 768                     FatalErrorIn
 769                     (
 770                         "meshCutAndRemove::setRefinement("
 771                         ", const label, const cellCuts&, const labelList&"
 772                         ", polyTopoChange&)"
 773                     )   << "Problem: point is marked as cut but"
 774                         << " not used by any loop"
 775                         << " or cell anchor point"
 776                         << "point:" << pointI
 777                         << " coord:" << mesh().points()[pointI]
 778                         << abort(FatalError);
 779                 }
 780             }
 781         }
 782
 783
 784         // Remove unused points.
 785         forAll(usedPoint, pointI)
 786         {
 787             if (!usedPoint[pointI])
 788             {
 789                 meshMod.setAction(polyRemovePoint(pointI));
 790
 791                 if (debug & 2)
 792                 {
 793                     Pout<< "Removing unused point " << pointI << endl;
 794                 }
 795             }
 796         }
 797     }
 798
 799
 800     //
 801     // For all cut cells add an internal or external face
 802     //
 803
 804     forAll(cellLoops, cellI)
 805     {
 806         const labelList& loop = cellLoops[cellI];
 807
 808         if (loop.size())
 809         {
 810             if (cutPatch[cellI] < 0 || cutPatch[cellI] >= patches.size())
 811             {
 812                 FatalErrorIn
 813                 (
 814                     "meshCutAndRemove::setRefinement("
 815                     ", const label, const cellCuts&, const labelList&"
 816                     ", polyTopoChange&)"
 817                 )   << "Illegal patch " << cutPatch[cellI]
 818                     << " provided for cut cell " << cellI
 819                     << abort(FatalError);
 820             }
 821
 822             //
 823             // Convert loop (=list of cuts) into proper face.
 824             // cellCuts sets orientation is towards anchor side so reverse.
 825             //
 826             face newFace(loopToFace(cellI, loop));
 827
 828             reverse(newFace);
 829
 830             // Pick any anchor point on cell
 831             label masterPointI = findPatchFacePoint(newFace, exposedPatchI);
 832
 833             label addedFaceI =
 834                 meshMod.setAction
 835                 (
 836                     polyAddFace
 837                     (
 838                         newFace,                // face
 839                         cellI,                  // owner
 840                         -1,                     // neighbour
 841                         masterPointI,           // master point
 842                         -1,                     // master edge
 843                         -1,                     // master face for addition
 844                         false,                  // flux flip
 845                         cutPatch[cellI],        // patch for face
 846                         -1,                     // zone for face
 847                         false                   // face zone flip
 848                     )
 849                 );
 850
 851             addedFaces_.insert(cellI, addedFaceI);
 852
 853             if (debug & 2)
 854             {
 855                 Pout<< "Added splitting face " << newFace << " index:"
 856                     << addedFaceI << " from masterPoint:" << masterPointI
 857                     << " to owner " << cellI << " with anchors:"
 858                     << anchorPts[cellI]
 859                     << " from Loop:";
 860
 861                 // Gets edgeweights of loop
 862                 scalarField weights(loop.size());
 863                 forAll(loop, i)
 864                 {
 865                     label cut = loop[i];
 866
 867                     weights[i] =
 868                     (
 869                         isEdge(cut)
 870                       ? cuts.edgeWeight()[getEdge(cut)]
 871                       : -GREAT
 872                     );
 873                 }
 874                 writeCuts(Pout, loop, weights);
 875                 Pout<< endl;
 876             }
 877         }
 878     }
 879
 880
 881     //
 882     // Modify faces to use only anchorpoints and loop points
 883     // (so throw away part without anchorpoints)
 884     //
 885
 886
 887     // Maintain whether face has been updated (for -split edges
 888     // -new owner/neighbour)
 889     boolList faceUptodate(mesh().nFaces(), false);
 890
 891     const Map<edge>& faceSplitCuts = cuts.faceSplitCut();
 892
 893     forAllConstIter(Map<edge>, faceSplitCuts, iter)
 894     {
 895         label faceI = iter.key();
 896
 897         // Renumber face to include split edges.
 898         face newFace(addEdgeCutsToFace(faceI));
 899
 900         // Edge splitting the face. Convert edge to new vertex numbering.
 901         const edge& splitEdge = iter();
 902
 903         label cut0 = splitEdge[0];
 904
 905         label v0;
 906         if (isEdge(cut0))
 907         {
 908             label edgeI = getEdge(cut0);
 909             v0 = addedPoints_[mesh().edges()[edgeI]];
 910         }
 911         else
 912         {
 913             v0 = getVertex(cut0);
 914         }
 915
 916         label cut1 = splitEdge[1];
 917         label v1;
 918         if (isEdge(cut1))
 919         {
 920             label edgeI = getEdge(cut1);
 921             v1 = addedPoints_[mesh().edges()[edgeI]];
 922         }
 923         else
 924         {
 925             v1 = getVertex(cut1);
 926         }
 927
 928         // Split face along the elements of the splitEdge.
 929         face f0, f1;
 930         splitFace(newFace, v0, v1, f0, f1);
 931
 932         label own = mesh().faceOwner()[faceI];
 933
 934         label nei = -1;
 935
 936         if (mesh().isInternalFace(faceI))
 937         {
 938             nei = mesh().faceNeighbour()[faceI];
 939         }
 940
 941         if (debug & 2)
 942         {
 943             Pout<< "Split face " << mesh().faces()[faceI]
 944                 << " own:" << own << " nei:" << nei
 945                 << " into f0:" << f0
 946                 << " and f1:" << f1 << endl;
 947         }
 948
 949
 950         // Check which cell using face uses anchorPoints (so is kept)
 951         // and which one doesn't (gets removed)
 952
 953         // Bit tricky. We have to know whether this faceSplit splits owner/
 954         // neighbour or both. Even if cell is cut we have to make sure this is
 955         // the one that cuts it (this face cut might not be the one splitting
 956         // the cell)
 957         // The face f gets split into two parts, f0 and f1.
 958         // Each of these can have a different owner and or neighbour.
 959
 960         const face& f = mesh().faces()[faceI];
 961
 962         label f0Own = -1;
 963         label f1Own = -1;
 964
 965         if (cellLoops[own].empty())
 966         {
 967             // Owner side is not split so keep both halves.
 968             f0Own = own;
 969             f1Own = own;
 970         }
 971         else if (isIn(splitEdge, cellLoops[own]))
 972         {
 973             // Owner is cut by this splitCut. See which of f0, f1 gets
 974             // preserved and becomes owner, and which gets removed.
 975             if (firstCommon(f0, anchorPts[own]) != -1)
 976             {
 977                 // f0 preserved so f1 gets deleted
 978                 f0Own = own;
 979                 f1Own = -1;
 980             }
 981             else
 982             {
 983                 f0Own = -1;
 984                 f1Own = own;
 985             }
 986         }
 987         else
 988         {
 989             // Owner not cut by this splitCut but by another.
 990             // Check on original face whether
 991             // use anchorPts.
 992             if (firstCommon(f, anchorPts[own]) != -1)
 993             {
 994                 // both f0 and f1 owner side preserved
 995                 f0Own = own;
 996                 f1Own = own;
 997             }
 998             else
 999             {
1000                 // both f0 and f1 owner side removed
1001                 f0Own = -1;
1002                 f1Own = -1;
1003             }
1004         }
1005
1006
1007         label f0Nei = -1;
1008         label f1Nei = -1;
1009
1010         if (nei != -1)
1011         {
1012             if (cellLoops[nei].empty())
1013             {
1014                 f0Nei = nei;
1015                 f1Nei = nei;
1016             }
1017             else if (isIn(splitEdge, cellLoops[nei]))
1018             {
1019                 // Neighbour is cut by this splitCut. So anchor part of it
1020                 // gets kept, non-anchor bit gets removed. See which of f0, f1
1021                 // connects to which part.
1022
1023                 if (firstCommon(f0, anchorPts[nei]) != -1)
1024                 {
1025                     f0Nei = nei;
1026                     f1Nei = -1;
1027                 }
1028                 else
1029                 {
1030                     f0Nei = -1;
1031                     f1Nei = nei;
1032                 }
1033             }
1034             else
1035             {
1036                 // neighbour not cut by this splitCut. Check on original face
1037                 // whether use anchorPts.
1038
1039                 if (firstCommon(f, anchorPts[nei]) != -1)
1040                 {
1041                     f0Nei = nei;
1042                     f1Nei = nei;
1043                 }
1044                 else
1045                 {
1046                     // both f0 and f1 on neighbour side removed
1047                     f0Nei = -1;
1048                     f1Nei = -1;
1049                 }
1050             }
1051         }
1052
1053
1054         if (debug & 2)
1055         {
1056             Pout<< "f0 own:" << f0Own << " nei:" << f0Nei
1057                 << "  f1 own:" << f1Own << " nei:" << f1Nei
1058                 << endl;
1059         }
1060
1061
1062         // If faces were internal but now become external set a patch.
1063         // If they were external already keep the patch.
1064         label patchID = patches.whichPatch(faceI);
1065
1066         if (patchID == -1)
1067         {
1068             patchID = exposedPatchI;
1069         }
1070
1071
1072         // Do as much as possible by modifying faceI. Delay any remove
1073         // face. Keep track of whether faceI has been used.
1074
1075         bool modifiedFaceI = false;
1076
1077         if (f0Own == -1)
1078         {
1079             if (f0Nei != -1)
1080             {
1081                 // f0 becomes external face (note:modFace will reverse face)
1082                 modFace(meshMod, faceI, f0, f0Own, f0Nei, patchID);
1083                 modifiedFaceI = true;
1084             }
1085         }
1086         else
1087         {
1088             if (f0Nei == -1)
1089             {
1090                 // f0 becomes external face
1091                 modFace(meshMod, faceI, f0, f0Own, f0Nei, patchID);
1092                 modifiedFaceI = true;
1093             }
1094             else
1095             {
1096                 // f0 stays internal face.
1097                 modFace(meshMod, faceI, f0, f0Own, f0Nei, -1);
1098                 modifiedFaceI = true;
1099             }
1100         }
1101
1102
1103         // f1 is added face (if at all)
1104
1105         if (f1Own == -1)
1106         {
1107             if (f1Nei == -1)
1108             {
1109                 // f1 not needed.
1110             }
1111             else
1112             {
1113                 // f1 becomes external face (note:modFace will reverse face)
1114                 if (!modifiedFaceI)
1115                 {
1116                     modFace(meshMod, faceI, f1, f1Own, f1Nei, patchID);
1117                     modifiedFaceI = true;
1118                 }
1119                 else
1120                 {
1121                     label masterPointI = findPatchFacePoint(f1, patchID);
1122
1123                     addFace
1124                     (
1125                         meshMod,
1126                         faceI,          // face for zone info
1127                         masterPointI,   // inflation point
1128                         f1,             // vertices of face
1129                         f1Own,
1130                         f1Nei,
1131                         patchID         // patch for new face
1132                     );
1133                 }
1134             }
1135         }
1136         else
1137         {
1138             if (f1Nei == -1)
1139             {
1140                 // f1 becomes external face
1141                 if (!modifiedFaceI)
1142                 {
1143                     modFace(meshMod, faceI, f1, f1Own, f1Nei, patchID);
1144                     modifiedFaceI = true;
1145                 }
1146                 else
1147                 {
1148                     label masterPointI = findPatchFacePoint(f1, patchID);
1149
1150                     addFace
1151                     (
1152                         meshMod,
1153                         faceI,
1154                         masterPointI,
1155                         f1,
1156                         f1Own,
1157                         f1Nei,
1158                         patchID
1159                     );
1160                 }
1161             }
1162             else
1163             {
1164                 // f1 is internal face.
1165                 if (!modifiedFaceI)
1166                 {
1167                     modFace(meshMod, faceI, f1, f1Own, f1Nei, -1);
1168                     modifiedFaceI = true;
1169                 }
1170                 else
1171                 {
1172                     label masterPointI = findPatchFacePoint(f1, -1);
1173
1174                     addFace(meshMod, faceI, masterPointI, f1, f1Own, f1Nei, -1);
1175                 }
1176             }
1177         }
1178
1179         if (f0Own == -1 && f0Nei == -1 && !modifiedFaceI)
1180         {
1181             meshMod.setAction(polyRemoveFace(faceI));
1182
1183             if (debug & 2)
1184             {
1185                 Pout<< "Removed face " << faceI << endl;
1186             }
1187         }
1188
1189         faceUptodate[faceI] = true;
1190     }
1191
1192
1193     //
1194     // Faces that have not been split but just appended to. Are guaranteed
1195     // to be reachable from an edgeCut.
1196     //
1197
1198     const boolList& edgeIsCut = cuts.edgeIsCut();
1199
1200     forAll(edgeIsCut, edgeI)
1201     {
1202         if (edgeIsCut[edgeI])
1203         {
1204             const labelList& eFaces = mesh().edgeFaces()[edgeI];
1205
1206             forAll(eFaces, i)
1207             {
1208                 label faceI = eFaces[i];
1209
1210                 if (!faceUptodate[faceI])
1211                 {
1212                     // So the face has not been split itself (i.e. its owner
1213                     // or neighbour have not been split) so it only
1214                     // borders by edge a cell which has been split.
1215
1216                     // Get (new or original) owner and neighbour of faceI
1217                     label own, nei, patchID;
1218                     faceCells(cuts, exposedPatchI, faceI, own, nei, patchID);
1219
1220
1221                     if (own == -1 && nei == -1)
1222                     {
1223                         meshMod.setAction(polyRemoveFace(faceI));
1224
1225                         if (debug & 2)
1226                         {
1227                             Pout<< "Removed face " << faceI << endl;
1228                         }
1229                     }
1230                     else
1231                     {
1232                         // Renumber face to include split edges.
1233                         face newFace(addEdgeCutsToFace(faceI));
1234
1235                         if (debug & 2)
1236                         {
1237                             Pout<< "Added edge cuts to face " << faceI
1238                                 << " f:" << mesh().faces()[faceI]
1239                                 << " newFace:" << newFace << endl;
1240                         }
1241
1242                         modFace
1243                         (
1244                             meshMod,
1245                             faceI,
1246                             newFace,
1247                             own,
1248                             nei,
1249                             patchID
1250                         );
1251                     }
1252
1253                     faceUptodate[faceI] = true;
1254                 }
1255             }
1256         }
1257     }
1258
1259
1260     //
1261     // Remove any faces on the non-anchor side of a split cell.
1262     // Note: could loop through all cut cells only and check their faces but
1263     //       looping over all faces is cleaner and probably faster for dense
1264     //       cut patterns.
1265
1266     const faceList& faces = mesh().faces();
1267
1268     forAll(faces, faceI)
1269     {
1270         if (!faceUptodate[faceI])
1271         {
1272             // Get (new or original) owner and neighbour of faceI
1273             label own, nei, patchID;
1274             faceCells(cuts, exposedPatchI, faceI, own, nei, patchID);
1275
1276             if (own == -1 && nei == -1)
1277             {
1278                 meshMod.setAction(polyRemoveFace(faceI));
1279
1280                 if (debug & 2)
1281                 {
1282                     Pout<< "Removed face " << faceI << endl;
1283                 }
1284             }
1285             else
1286             {
1287                 modFace(meshMod, faceI, faces[faceI], own, nei, patchID);
1288             }
1289
1290             faceUptodate[faceI] = true;
1291         }
1292     }
1293
1294     if (debug)
1295     {
1296         Pout<< "meshCutAndRemove:" << nl
1297             << "    cells split:" << cuts.nLoops() << nl
1298             << "    faces added:" << addedFaces_.size() << nl
1299             << "    points added on edges:" << addedPoints_.size() << nl
1300             << endl;
1301     }
1302 }
1303
1304
1305 void Foam::meshCutAndRemove::updateMesh(const mapPolyMesh& map)
1306 {
1307     // Update stored labels for mesh change.
1308     {
1309         Map<label> newAddedFaces(addedFaces_.size());
1310
1311         forAllConstIter(Map<label>, addedFaces_, iter)
1312         {
1313             label cellI = iter.key();
1314             label newCellI = map.reverseCellMap()[cellI];
1315
1316             label addedFaceI = iter();
1317
1318             label newAddedFaceI = map.reverseFaceMap()[addedFaceI];
1319
1320             if ((newCellI >= 0) && (newAddedFaceI >= 0))
1321             {
1322                 if
1323                 (
1324                     (debug & 2)
1325                  && (newCellI != cellI || newAddedFaceI != addedFaceI)
1326                 )
1327                 {
1328                     Pout<< "meshCutAndRemove::updateMesh :"
1329                         << " updating addedFace for cell " << cellI
1330                         << " from " << addedFaceI
1331                         << " to " << newAddedFaceI
1332                         << endl;
1333                 }
1334                 newAddedFaces.insert(newCellI, newAddedFaceI);
1335             }
1336         }
1337
1338         // Copy
1339         addedFaces_.transfer(newAddedFaces);
1340     }
1341
1342     {
1343         HashTable<label, edge, Hash<edge> > newAddedPoints(addedPoints_.size());
1344
1345         for
1346         (
1347             HashTable<label, edge, Hash<edge> >::const_iterator iter =
1348                 addedPoints_.begin();
1349             iter != addedPoints_.end();
1350             ++iter
1351         )
1352         {
1353             const edge& e = iter.key();
1354
1355             label newStart = map.reversePointMap()[e.start()];
1356
1357             label newEnd = map.reversePointMap()[e.end()];
1358
1359             label addedPointI = iter();
1360
1361             label newAddedPointI = map.reversePointMap()[addedPointI];
1362
1363             if ((newStart >= 0) && (newEnd >= 0) && (newAddedPointI >= 0))
1364             {
1365                 edge newE = edge(newStart, newEnd);
1366
1367                 if
1368                 (
1369                     (debug & 2)
1370                  && (e != newE || newAddedPointI != addedPointI)
1371                 )
1372                 {
1373                     Pout<< "meshCutAndRemove::updateMesh :"
1374                         << " updating addedPoints for edge " << e
1375                         << " from " << addedPointI
1376                         << " to " << newAddedPointI
1377                         << endl;
1378                 }
1379
1380                 newAddedPoints.insert(newE, newAddedPointI);
1381             }
1382         }
1383
1384         // Copy
1385         addedPoints_.transfer(newAddedPoints);
1386     }
1387 }
1388
1389
1390 // ************************************************************************* //