src/dynamicMesh/meshCut/meshModifiers/meshCutter/meshCutter.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 "meshCutter.H"
  27 #include "polyMesh.H"
  28 #include "polyTopoChange.H"
  29 #include "cellCuts.H"
  30 #include "mapPolyMesh.H"
  31 #include "meshTools.H"
  32 #include "polyModifyFace.H"
  33 #include "polyAddPoint.H"
  34 #include "polyAddFace.H"
  35 #include "polyAddCell.H"
  36
  37 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  38
  39 defineTypeNameAndDebug(Foam::meshCutter, 0);
  40
  41
  42 // * * * * * * * * * * * * * Private Static Functions  * * * * * * * * * * * //
  43
  44 // Returns true if lst1 and lst2 share elements
  45 bool Foam::meshCutter::uses(const labelList& elems1, const labelList& elems2)
  46 {
  47     forAll(elems1, elemI)
  48     {
  49         if (findIndex(elems2, elems1[elemI]) != -1)
  50         {
  51             return true;
  52         }
  53     }
  54     return false;
  55 }
  56
  57
  58 // Check if twoCuts at two consecutive position in cuts.
  59 bool Foam::meshCutter::isIn
  60 (
  61     const edge& twoCuts,
  62     const labelList& cuts
  63 )
  64 {
  65     label index = findIndex(cuts, twoCuts[0]);
  66
  67     if (index == -1)
  68     {
  69         return false;
  70     }
  71
  72     return
  73     (
  74         cuts[cuts.fcIndex(index)] == twoCuts[1]
  75      || cuts[cuts.rcIndex(index)] == twoCuts[1]
  76     );
  77 }
  78
  79
  80 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  81
  82 // Returns the cell in cellLabels that is cut. Or -1.
  83 Foam::label Foam::meshCutter::findCutCell
  84 (
  85     const cellCuts& cuts,
  86     const labelList& cellLabels
  87 ) const
  88 {
  89     forAll(cellLabels, labelI)
  90     {
  91         label cellI = cellLabels[labelI];
  92
  93         if (cuts.cellLoops()[cellI].size())
  94         {
  95             return cellI;
  96         }
  97     }
  98     return -1;
  99 }
 100
 101
 102 //- Returns first pointI in pointLabels that uses an internal
 103 //  face. Used to find point to inflate cell/face from (has to be
 104 //  connected to internal face). Returns -1 (so inflate from nothing) if
 105 //  none found.
 106 Foam::label Foam::meshCutter::findInternalFacePoint
 107 (
 108     const labelList& pointLabels
 109 ) const
 110 {
 111     forAll(pointLabels, labelI)
 112     {
 113         label pointI = pointLabels[labelI];
 114
 115         const labelList& pFaces = mesh().pointFaces()[pointI];
 116
 117         forAll(pFaces, pFaceI)
 118         {
 119             label faceI = pFaces[pFaceI];
 120
 121             if (mesh().isInternalFace(faceI))
 122             {
 123                 return pointI;
 124             }
 125         }
 126     }
 127
 128     if (pointLabels.empty())
 129     {
 130         FatalErrorIn("meshCutter::findInternalFacePoint(const labelList&)")
 131             << "Empty pointLabels" << abort(FatalError);
 132     }
 133
 134     return -1;
 135 }
 136
 137
 138 // Get new owner and neighbour of face. Checks anchor points to see if
 139 // need to get original or added cell.
 140 void Foam::meshCutter::faceCells
 141 (
 142     const cellCuts& cuts,
 143     const label faceI,
 144     label& own,
 145     label& nei
 146 ) const
 147 {
 148     const labelListList& anchorPts = cuts.cellAnchorPoints();
 149     const labelListList& cellLoops = cuts.cellLoops();
 150
 151     const face& f = mesh().faces()[faceI];
 152
 153     own = mesh().faceOwner()[faceI];
 154
 155     if (cellLoops[own].size() && uses(f, anchorPts[own]))
 156     {
 157         own = addedCells_[own];
 158     }
 159
 160     nei = -1;
 161
 162     if (mesh().isInternalFace(faceI))
 163     {
 164         nei = mesh().faceNeighbour()[faceI];
 165
 166         if (cellLoops[nei].size() && uses(f, anchorPts[nei]))
 167         {
 168             nei = addedCells_[nei];
 169         }
 170     }
 171 }
 172
 173
 174 void Foam::meshCutter::getFaceInfo
 175 (
 176     const label faceI,
 177     label& patchID,
 178     label& zoneID,
 179     label& zoneFlip
 180 ) const
 181 {
 182     patchID = -1;
 183
 184     if (!mesh().isInternalFace(faceI))
 185     {
 186         patchID = mesh().boundaryMesh().whichPatch(faceI);
 187     }
 188
 189     zoneID = mesh().faceZones().whichZone(faceI);
 190
 191     zoneFlip = false;
 192
 193     if (zoneID >= 0)
 194     {
 195         const faceZone& fZone = mesh().faceZones()[zoneID];
 196
 197         zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
 198     }
 199 }
 200
 201
 202 // Adds a face on top of existing faceI.
 203 void Foam::meshCutter::addFace
 204 (
 205     polyTopoChange& meshMod,
 206     const label faceI,
 207     const face& newFace,
 208     const label own,
 209     const label nei
 210 )
 211 {
 212     label patchID, zoneID, zoneFlip;
 213
 214     getFaceInfo(faceI, patchID, zoneID, zoneFlip);
 215
 216     if ((nei == -1) || (own < nei))
 217     {
 218         // Ordering ok.
 219         if (debug & 2)
 220         {
 221             Pout<< "Adding face " << newFace
 222                 << " with new owner:" << own
 223                 << " with new neighbour:" << nei
 224                 << " patchID:" << patchID
 225                 << " zoneID:" << zoneID
 226                 << " zoneFlip:" << zoneFlip
 227                 << endl;
 228         }
 229
 230         meshMod.setAction
 231         (
 232             polyAddFace
 233             (
 234                 newFace,                    // face
 235                 own,                        // owner
 236                 nei,                        // neighbour
 237                 -1,                         // master point
 238                 -1,                         // master edge
 239                 faceI,                      // master face for addition
 240                 false,                      // flux flip
 241                 patchID,                    // patch for face
 242                 zoneID,                     // zone for face
 243                 zoneFlip                    // face zone flip
 244             )
 245         );
 246     }
 247     else
 248     {
 249         // Reverse owner/neighbour
 250         if (debug & 2)
 251         {
 252             Pout<< "Adding (reversed) face " << newFace.reverseFace()
 253                 << " with new owner:" << nei
 254                 << " with new neighbour:" << own
 255                 << " patchID:" << patchID
 256                 << " zoneID:" << zoneID
 257                 << " zoneFlip:" << zoneFlip
 258                 << endl;
 259         }
 260
 261         meshMod.setAction
 262         (
 263             polyAddFace
 264             (
 265                 newFace.reverseFace(),      // face
 266                 nei,                        // owner
 267                 own,                        // neighbour
 268                 -1,                         // master point
 269                 -1,                         // master edge
 270                 faceI,                      // master face for addition
 271                 false,                      // flux flip
 272                 patchID,                    // patch for face
 273                 zoneID,                     // zone for face
 274                 zoneFlip                    // face zone flip
 275             )
 276         );
 277     }
 278 }
 279
 280
 281 // Modifies existing faceI for either new owner/neighbour or new face points.
 282 void Foam::meshCutter::modFace
 283 (
 284     polyTopoChange& meshMod,
 285     const label faceI,
 286     const face& newFace,
 287     const label own,
 288     const label nei
 289 )
 290 {
 291     label patchID, zoneID, zoneFlip;
 292
 293     getFaceInfo(faceI, patchID, zoneID, zoneFlip);
 294
 295     if
 296     (
 297         (own != mesh().faceOwner()[faceI])
 298      || (
 299             mesh().isInternalFace(faceI)
 300          && (nei != mesh().faceNeighbour()[faceI])
 301         )
 302      || (newFace != mesh().faces()[faceI])
 303     )
 304     {
 305         if (debug & 2)
 306         {
 307             Pout<< "Modifying face " << faceI
 308                 << " old vertices:" << mesh().faces()[faceI]
 309                 << " new vertices:" << newFace
 310                 << " new owner:" << own
 311                 << " new neighbour:" << nei
 312                 << " new zoneID:" << zoneID
 313                 << " new zoneFlip:" << zoneFlip
 314                 << endl;
 315         }
 316
 317         if ((nei == -1) || (own < nei))
 318         {
 319             meshMod.setAction
 320             (
 321                 polyModifyFace
 322                 (
 323                     newFace,            // modified face
 324                     faceI,              // label of face being modified
 325                     own,                // owner
 326                     nei,                // neighbour
 327                     false,              // face flip
 328                     patchID,            // patch for face
 329                     false,              // remove from zone
 330                     zoneID,             // zone for face
 331                     zoneFlip            // face flip in zone
 332                 )
 333             );
 334         }
 335         else
 336         {
 337             meshMod.setAction
 338             (
 339                 polyModifyFace
 340                 (
 341                     newFace.reverseFace(),  // modified face
 342                     faceI,                  // label of face being modified
 343                     nei,                    // owner
 344                     own,                    // neighbour
 345                     false,                  // face flip
 346                     patchID,                // patch for face
 347                     false,                  // remove from zone
 348                     zoneID,                 // zone for face
 349                     zoneFlip                // face flip in zone
 350                 )
 351             );
 352         }
 353     }
 354 }
 355
 356
 357 // Copies face starting from startFp up to and including endFp.
 358 void Foam::meshCutter::copyFace
 359 (
 360     const face& f,
 361     const label startFp,
 362     const label endFp,
 363     face& newFace
 364 ) const
 365 {
 366     label fp = startFp;
 367
 368     label newFp = 0;
 369
 370     while (fp != endFp)
 371     {
 372         newFace[newFp++] = f[fp];
 373
 374         fp = (fp + 1) % f.size();
 375     }
 376     newFace[newFp] = f[fp];
 377 }
 378
 379
 380 // Actually split face in two along splitEdge v0, v1 (the two vertices in new
 381 // vertex numbering). Generates faces in same ordering
 382 // as original face. Replaces cutEdges by the points introduced on them
 383 // (addedPoints_).
 384 void Foam::meshCutter::splitFace
 385 (
 386     const face& f,
 387     const label v0,
 388     const label v1,
 389
 390     face& f0,
 391     face& f1
 392 ) const
 393 {
 394     // Check if we find any new vertex which is part of the splitEdge.
 395     label startFp = findIndex(f, v0);
 396
 397     if (startFp == -1)
 398     {
 399         FatalErrorIn
 400         (
 401             "meshCutter::splitFace"
 402             ", const face&, const label, const label, face&, face&)"
 403         )   << "Cannot find vertex (new numbering) " << v0
 404             << " on face " << f
 405             << abort(FatalError);
 406     }
 407
 408     label endFp = findIndex(f, v1);
 409
 410     if (endFp == -1)
 411     {
 412         FatalErrorIn
 413         (
 414             "meshCutter::splitFace("
 415             ", const face&, const label, const label, face&, face&)"
 416         )   << "Cannot find vertex (new numbering) " << v1
 417             << " on face " << f
 418             << abort(FatalError);
 419     }
 420
 421
 422     f0.setSize((endFp + 1 + f.size() - startFp) % f.size());
 423     f1.setSize(f.size() - f0.size() + 2);
 424
 425     copyFace(f, startFp, endFp, f0);
 426     copyFace(f, endFp, startFp, f1);
 427 }
 428
 429
 430 // Adds additional vertices (from edge cutting) to face. Used for faces which
 431 // are not split but still might use edge that has been cut.
 432 Foam::face Foam::meshCutter::addEdgeCutsToFace(const label faceI) const
 433 {
 434     const face& f = mesh().faces()[faceI];
 435
 436     face newFace(2 * f.size());
 437
 438     label newFp = 0;
 439
 440     forAll(f, fp)
 441     {
 442         // Duplicate face vertex .
 443         newFace[newFp++] = f[fp];
 444
 445         // Check if edge has been cut.
 446         label fp1 = f.fcIndex(fp);
 447
 448         HashTable<label, edge, Hash<edge> >::const_iterator fnd =
 449             addedPoints_.find(edge(f[fp], f[fp1]));
 450
 451         if (fnd != addedPoints_.end())
 452         {
 453             // edge has been cut. Introduce new vertex.
 454             newFace[newFp++] = fnd();
 455         }
 456     }
 457
 458     newFace.setSize(newFp);
 459
 460     return newFace;
 461 }
 462
 463
 464 // Walk loop (loop of cuts) across circumference of cellI. Returns face in
 465 // new vertices.
 466 // Note: tricky bit is that it can use existing edges which have been split.
 467 Foam::face Foam::meshCutter::loopToFace
 468 (
 469     const label cellI,
 470     const labelList& loop
 471 ) const
 472 {
 473     face newFace(2*loop.size());
 474
 475     label newFaceI = 0;
 476
 477     forAll(loop, fp)
 478     {
 479         label cut = loop[fp];
 480
 481         if (isEdge(cut))
 482         {
 483             label edgeI = getEdge(cut);
 484
 485             const edge& e = mesh().edges()[edgeI];
 486
 487             label vertI = addedPoints_[e];
 488
 489             newFace[newFaceI++] = vertI;
 490         }
 491         else
 492         {
 493             // cut is vertex.
 494             label vertI = getVertex(cut);
 495
 496             newFace[newFaceI++] = vertI;
 497
 498             label nextCut = loop[loop.fcIndex(fp)];
 499
 500             if (!isEdge(nextCut))
 501             {
 502                 // From vertex to vertex -> cross cut only if no existing edge.
 503                 label nextVertI = getVertex(nextCut);
 504
 505                 label edgeI = meshTools::findEdge(mesh(), vertI, nextVertI);
 506
 507                 if (edgeI != -1)
 508                 {
 509                     // Existing edge. Insert split-edge point if any.
 510                     HashTable<label, edge, Hash<edge> >::const_iterator fnd =
 511                         addedPoints_.find(mesh().edges()[edgeI]);
 512
 513                     if (fnd != addedPoints_.end())
 514                     {
 515                         newFace[newFaceI++] = fnd();
 516                     }
 517                 }
 518             }
 519         }
 520     }
 521     newFace.setSize(newFaceI);
 522
 523     return newFace;
 524 }
 525
 526
 527 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 528
 529 // Construct from components
 530 Foam::meshCutter::meshCutter(const polyMesh& mesh)
 531 :
 532     edgeVertex(mesh),
 533     addedCells_(),
 534     addedFaces_(),
 535     addedPoints_()
 536
 537 {}
 538
 539
 540 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 541
 542 Foam::meshCutter::~meshCutter()
 543 {}
 544
 545
 546 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 547
 548 void Foam::meshCutter::setRefinement
 549 (
 550     const cellCuts& cuts,
 551     polyTopoChange& meshMod
 552 )
 553 {
 554     // Clear and size maps here since mesh size will change.
 555     addedCells_.clear();
 556     addedCells_.resize(cuts.nLoops());
 557
 558     addedFaces_.clear();
 559     addedFaces_.resize(cuts.nLoops());
 560
 561     addedPoints_.clear();
 562     addedPoints_.resize(cuts.nLoops());
 563
 564     if (cuts.nLoops() == 0)
 565     {
 566         return;
 567     }
 568
 569     const labelListList& anchorPts = cuts.cellAnchorPoints();
 570     const labelListList& cellLoops = cuts.cellLoops();
 571
 572     //
 573     // Add new points along cut edges.
 574     //
 575
 576     forAll(cuts.edgeIsCut(), edgeI)
 577     {
 578         if (cuts.edgeIsCut()[edgeI])
 579         {
 580             const edge& e = mesh().edges()[edgeI];
 581
 582             // Check if there is any cell using this edge.
 583             if (debug && findCutCell(cuts, mesh().edgeCells()[edgeI]) == -1)
 584             {
 585                 FatalErrorIn
 586                 (
 587                     "meshCutter::setRefinement(const cellCuts&"
 588                     ", polyTopoChange&)"
 589                 )   << "Problem: cut edge but none of the cells using it is\n"
 590                     << "edge:" << edgeI << " verts:" << e
 591                     << abort(FatalError);
 592             }
 593
 594             // One of the edge end points should be master point of nbCellI.
 595             label masterPointI = e.start();
 596
 597             const point& v0 = mesh().points()[e.start()];
 598             const point& v1 = mesh().points()[e.end()];
 599
 600             scalar weight = cuts.edgeWeight()[edgeI];
 601
 602             point newPt = weight*v1 + (1.0-weight)*v0;
 603
 604             label addedPointI =
 605                 meshMod.setAction
 606                 (
 607                     polyAddPoint
 608                     (
 609                         newPt,              // point
 610                         masterPointI,       // master point
 611                         -1,                 // zone for point
 612                         true                // supports a cell
 613                     )
 614                 );
 615
 616             // Store on (hash of) edge.
 617             addedPoints_.insert(e, addedPointI);
 618
 619             if (debug & 2)
 620             {
 621                 Pout<< "Added point " << addedPointI
 622                     << " to vertex "
 623                     << masterPointI << " of edge " << edgeI
 624                     << " vertices " << e << endl;
 625             }
 626         }
 627     }
 628
 629     //
 630     // Add cells (on 'anchor' side of cell)
 631     //
 632
 633     forAll(cellLoops, cellI)
 634     {
 635         if (cellLoops[cellI].size())
 636         {
 637             // Add a cell to the existing cell
 638             label addedCellI =
 639                 meshMod.setAction
 640                 (
 641                     polyAddCell
 642                     (
 643                         -1,                 // master point
 644                         -1,                 // master edge
 645                         -1,                 // master face
 646                         cellI,              // master cell
 647                         -1                  // zone for cell
 648                     )
 649                 );
 650
 651             addedCells_.insert(cellI, addedCellI);
 652
 653             if (debug & 2)
 654             {
 655                 Pout<< "Added cell " << addedCells_[cellI] << " to cell "
 656                     << cellI << endl;
 657             }
 658         }
 659     }
 660
 661
 662     //
 663     // For all cut cells add an internal face
 664     //
 665
 666     forAll(cellLoops, cellI)
 667     {
 668         const labelList& loop = cellLoops[cellI];
 669
 670         if (loop.size())
 671         {
 672             //
 673             // Convert loop (=list of cuts) into proper face.
 674             // Orientation should already be ok. (done by cellCuts)
 675             //
 676             face newFace(loopToFace(cellI, loop));
 677
 678             // Pick any anchor point on cell
 679             label masterPointI = findInternalFacePoint(anchorPts[cellI]);
 680
 681             label addedFaceI =
 682                 meshMod.setAction
 683                 (
 684                     polyAddFace
 685                     (
 686                         newFace,                // face
 687                         cellI,                  // owner
 688                         addedCells_[cellI],     // neighbour
 689                         masterPointI,           // master point
 690                         -1,                     // master edge
 691                         -1,                     // master face for addition
 692                         false,                  // flux flip
 693                         -1,                     // patch for face
 694                         -1,                     // zone for face
 695                         false                   // face zone flip
 696                     )
 697                 );
 698
 699             addedFaces_.insert(cellI, addedFaceI);
 700
 701             if (debug & 2)
 702             {
 703                 // Gets edgeweights of loop
 704                 scalarField weights(loop.size());
 705                 forAll(loop, i)
 706                 {
 707                     label cut = loop[i];
 708
 709                     weights[i] =
 710                     (
 711                         isEdge(cut)
 712                       ? cuts.edgeWeight()[getEdge(cut)]
 713                       : -GREAT
 714                     );
 715                 }
 716
 717                 Pout<< "Added splitting face " << newFace << " index:"
 718                     << addedFaceI
 719                     << " to owner " << cellI
 720                     << " neighbour " << addedCells_[cellI]
 721                     << " from Loop:";
 722                 writeCuts(Pout, loop, weights);
 723                 Pout<< endl;
 724             }
 725         }
 726     }
 727
 728
 729     //
 730     // Modify faces on the outside and create new ones
 731     // (in effect split old faces into two)
 732     //
 733
 734     // Maintain whether face has been updated (for -split edges
 735     // -new owner/neighbour)
 736     boolList faceUptodate(mesh().nFaces(), false);
 737
 738     const Map<edge>& faceSplitCuts = cuts.faceSplitCut();
 739
 740     forAllConstIter(Map<edge>, faceSplitCuts, iter)
 741     {
 742         label faceI = iter.key();
 743
 744         // Renumber face to include split edges.
 745         face newFace(addEdgeCutsToFace(faceI));
 746
 747         // Edge splitting the face. Convert cuts to new vertex numbering.
 748         const edge& splitEdge = iter();
 749
 750         label cut0 = splitEdge[0];
 751
 752         label v0;
 753         if (isEdge(cut0))
 754         {
 755             label edgeI = getEdge(cut0);
 756             v0 = addedPoints_[mesh().edges()[edgeI]];
 757         }
 758         else
 759         {
 760             v0 = getVertex(cut0);
 761         }
 762
 763         label cut1 = splitEdge[1];
 764         label v1;
 765         if (isEdge(cut1))
 766         {
 767             label edgeI = getEdge(cut1);
 768             v1 = addedPoints_[mesh().edges()[edgeI]];
 769         }
 770         else
 771         {
 772             v1 = getVertex(cut1);
 773         }
 774
 775         // Split face along the elements of the splitEdge.
 776         face f0, f1;
 777         splitFace(newFace, v0, v1, f0, f1);
 778
 779         label own = mesh().faceOwner()[faceI];
 780
 781         label nei = -1;
 782
 783         if (mesh().isInternalFace(faceI))
 784         {
 785             nei = mesh().faceNeighbour()[faceI];
 786         }
 787
 788         if (debug & 2)
 789         {
 790             Pout<< "Split face " << mesh().faces()[faceI]
 791                 << " own:" << own << " nei:" << nei
 792                 << " into f0:" << f0
 793                 << " and f1:" << f1 << endl;
 794         }
 795
 796         // Check which face uses anchorPoints (connects to addedCell)
 797         // and which one doesn't (connects to original cell)
 798
 799         // Bit tricky. We have to know whether this faceSplit splits owner/
 800         // neighbour or both. Even if cell is cut we have to make sure this is
 801         // the one that cuts it (this face cut might not be the one splitting
 802         // the cell)
 803
 804         const face& f = mesh().faces()[faceI];
 805
 806         label f0Owner = -1;
 807         label f1Owner = -1;
 808
 809         if (cellLoops[own].empty())
 810         {
 811             f0Owner = own;
 812             f1Owner = own;
 813         }
 814         else if (isIn(splitEdge, cellLoops[own]))
 815         {
 816             // Owner is cut by this splitCut. See which of f0, f1 gets
 817             // owner, which gets addedCells_[owner]
 818             if (uses(f0, anchorPts[own]))
 819             {
 820                 f0Owner = addedCells_[own];
 821                 f1Owner = own;
 822             }
 823             else
 824             {
 825                 f0Owner = own;
 826                 f1Owner = addedCells_[own];
 827             }
 828         }
 829         else
 830         {
 831             // Owner not cut by this splitCut. Check on original face whether
 832             // use anchorPts.
 833             if (uses(f, anchorPts[own]))
 834             {
 835                 label newCellI = addedCells_[own];
 836                 f0Owner = newCellI;
 837                 f1Owner = newCellI;
 838             }
 839             else
 840             {
 841                 f0Owner = own;
 842                 f1Owner = own;
 843             }
 844         }
 845
 846
 847         label f0Neighbour = -1;
 848         label f1Neighbour = -1;
 849
 850         if (nei != -1)
 851         {
 852             if (cellLoops[nei].empty())
 853             {
 854                 f0Neighbour = nei;
 855                 f1Neighbour = nei;
 856             }
 857             else if (isIn(splitEdge, cellLoops[nei]))
 858             {
 859                 // Neighbour is cut by this splitCut. See which of f0, f1
 860                 // gets which neighbour/addedCells_[neighbour]
 861                 if (uses(f0, anchorPts[nei]))
 862                 {
 863                     f0Neighbour = addedCells_[nei];
 864                     f1Neighbour = nei;
 865                 }
 866                 else
 867                 {
 868                     f0Neighbour = nei;
 869                     f1Neighbour = addedCells_[nei];
 870                 }
 871             }
 872             else
 873             {
 874                 // neighbour not cut by this splitCut. Check on original face
 875                 // whether use anchorPts.
 876                 if (uses(f, anchorPts[nei]))
 877                 {
 878                     label newCellI = addedCells_[nei];
 879                     f0Neighbour = newCellI;
 880                     f1Neighbour = newCellI;
 881                 }
 882                 else
 883                 {
 884                     f0Neighbour = nei;
 885                     f1Neighbour = nei;
 886                 }
 887             }
 888         }
 889
 890         // f0 is the added face, f1 the modified one
 891         addFace(meshMod, faceI, f0, f0Owner, f0Neighbour);
 892
 893         modFace(meshMod, faceI, f1, f1Owner, f1Neighbour);
 894
 895         faceUptodate[faceI] = true;
 896     }
 897
 898
 899     //
 900     // Faces that have not been split but just appended to. Are guaranteed
 901     // to be reachable from an edgeCut.
 902     //
 903
 904     const boolList& edgeIsCut = cuts.edgeIsCut();
 905
 906     forAll(edgeIsCut, edgeI)
 907     {
 908         if (edgeIsCut[edgeI])
 909         {
 910             const labelList& eFaces = mesh().edgeFaces()[edgeI];
 911
 912             forAll(eFaces, i)
 913             {
 914                 label faceI = eFaces[i];
 915
 916                 if (!faceUptodate[faceI])
 917                 {
 918                     // Renumber face to include split edges.
 919                     face newFace(addEdgeCutsToFace(faceI));
 920
 921                     if (debug & 2)
 922                     {
 923                         Pout<< "Added edge cuts to face " << faceI
 924                             << " f:" << mesh().faces()[faceI]
 925                             << " newFace:" << newFace << endl;
 926                     }
 927
 928                     // Get (new or original) owner and neighbour of faceI
 929                     label own, nei;
 930                     faceCells(cuts, faceI, own, nei);
 931
 932                     modFace(meshMod, faceI, newFace, own, nei);
 933
 934                     faceUptodate[faceI] = true;
 935                 }
 936             }
 937         }
 938     }
 939
 940
 941     //
 942     // Correct any original faces on split cell for new neighbour/owner
 943     //
 944
 945     forAll(cellLoops, cellI)
 946     {
 947         if (cellLoops[cellI].size())
 948         {
 949             const labelList& cllFaces = mesh().cells()[cellI];
 950
 951             forAll(cllFaces, cllFaceI)
 952             {
 953                 label faceI = cllFaces[cllFaceI];
 954
 955                 if (!faceUptodate[faceI])
 956                 {
 957                     // Update face with new owner/neighbour (if any)
 958                     const face& f = mesh().faces()[faceI];
 959
 960                     if (debug && (f != addEdgeCutsToFace(faceI)))
 961                     {
 962                         FatalErrorIn
 963                         (
 964                             "meshCutter::setRefinement(const cellCuts&"
 965                             ", polyTopoChange&)"
 966                         )   << "Problem: edges added to face which does "
 967                             << " not use a marked cut" << endl
 968                             << "faceI:" << faceI << endl
 969                             << "face:" << f << endl
 970                             << "newFace:" << addEdgeCutsToFace(faceI)
 971                             << abort(FatalError);
 972                     }
 973
 974                     // Get (new or original) owner and neighbour of faceI
 975                     label own, nei;
 976                     faceCells(cuts, faceI, own, nei);
 977
 978                     modFace
 979                     (
 980                         meshMod,
 981                         faceI,
 982                         f,
 983                         own,
 984                         nei
 985                     );
 986
 987                     faceUptodate[faceI] = true;
 988                 }
 989             }
 990         }
 991     }
 992
 993     if (debug)
 994     {
 995         Pout<< "meshCutter:" << nl
 996             << "    cells split:" << addedCells_.size() << nl
 997             << "    faces added:" << addedFaces_.size() << nl
 998             << "    points added on edges:" << addedPoints_.size() << nl
 999             << endl;
1000     }
1001 }
1002
1003
1004 void Foam::meshCutter::updateMesh(const mapPolyMesh& morphMap)
1005 {
1006     // Update stored labels for mesh change.
1007
1008     {
1009         // Create copy since new label might (temporarily) clash with existing
1010         // key.
1011         Map<label> newAddedCells(addedCells_.size());
1012
1013         forAllConstIter(Map<label>, addedCells_, iter)
1014         {
1015             label cellI = iter.key();
1016             label newCellI = morphMap.reverseCellMap()[cellI];
1017
1018             label addedCellI = iter();
1019
1020             label newAddedCellI = morphMap.reverseCellMap()[addedCellI];
1021
1022             if (newCellI >= 0 && newAddedCellI >= 0)
1023             {
1024                 if
1025                 (
1026                     (debug & 2)
1027                  && (newCellI != cellI || newAddedCellI != addedCellI)
1028                 )
1029                 {
1030                     Pout<< "meshCutter::updateMesh :"
1031                         << " updating addedCell for cell " << cellI
1032                         << " from " << addedCellI
1033                         << " to " << newAddedCellI << endl;
1034                 }
1035                 newAddedCells.insert(newCellI, newAddedCellI);
1036             }
1037         }
1038
1039         // Copy
1040         addedCells_.transfer(newAddedCells);
1041     }
1042
1043     {
1044         Map<label> newAddedFaces(addedFaces_.size());
1045
1046         forAllConstIter(Map<label>, addedFaces_, iter)
1047         {
1048             label cellI = iter.key();
1049             label newCellI = morphMap.reverseCellMap()[cellI];
1050
1051             label addedFaceI = iter();
1052
1053             label newAddedFaceI = morphMap.reverseFaceMap()[addedFaceI];
1054
1055             if ((newCellI >= 0) && (newAddedFaceI >= 0))
1056             {
1057                 if
1058                 (
1059                     (debug & 2)
1060                  && (newCellI != cellI || newAddedFaceI != addedFaceI)
1061                 )
1062                 {
1063                     Pout<< "meshCutter::updateMesh :"
1064                         << " updating addedFace for cell " << cellI
1065                         << " from " << addedFaceI
1066                         << " to " << newAddedFaceI
1067                         << endl;
1068                 }
1069                 newAddedFaces.insert(newCellI, newAddedFaceI);
1070             }
1071         }
1072
1073         // Copy
1074         addedFaces_.transfer(newAddedFaces);
1075     }
1076
1077     {
1078         HashTable<label, edge, Hash<edge> > newAddedPoints(addedPoints_.size());
1079
1080         for
1081         (
1082             HashTable<label, edge, Hash<edge> >::const_iterator iter =
1083                 addedPoints_.begin();
1084             iter != addedPoints_.end();
1085             ++iter
1086         )
1087         {
1088             const edge& e = iter.key();
1089
1090             label newStart = morphMap.reversePointMap()[e.start()];
1091
1092             label newEnd = morphMap.reversePointMap()[e.end()];
1093
1094             label addedPointI = iter();
1095
1096             label newAddedPointI = morphMap.reversePointMap()[addedPointI];
1097
1098             if ((newStart >= 0) && (newEnd >= 0) && (newAddedPointI >= 0))
1099             {
1100                 edge newE = edge(newStart, newEnd);
1101
1102                 if
1103                 (
1104                     (debug & 2)
1105                  && (e != newE || newAddedPointI != addedPointI)
1106                 )
1107                 {
1108                     Pout<< "meshCutter::updateMesh :"
1109                         << " updating addedPoints for edge " << e
1110                         << " from " << addedPointI
1111                         << " to " << newAddedPointI
1112                         << endl;
1113                 }
1114
1115                 newAddedPoints.insert(newE, newAddedPointI);
1116             }
1117         }
1118
1119         // Copy
1120         addedPoints_.transfer(newAddedPoints);
1121     }
1122 }
1123
1124
1125 // ************************************************************************* //