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