src/dynamicMesh/polyMeshAdder/polyMeshAdder.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 "polyMeshAdder.H"
  27 #include "mapAddedPolyMesh.H"
  28 #include "IOobject.H"
  29 #include "faceCoupleInfo.H"
  30 #include "processorPolyPatch.H"
  31 #include "SortableList.H"
  32 #include "Time.H"
  33 #include "globalMeshData.H"
  34 #include "mergePoints.H"
  35 #include "polyModifyFace.H"
  36 #include "polyRemovePoint.H"
  37 #include "polyTopoChange.H"
  38
  39 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  40
  41 // Append all mapped elements of a list to a DynamicList
  42 void Foam::polyMeshAdder::append
  43 (
  44     const labelList& map,
  45     const labelList& lst,
  46     DynamicList<label>& dynLst
  47 )
  48 {
  49     dynLst.setCapacity(dynLst.size() + lst.size());
  50
  51     forAll(lst, i)
  52     {
  53         const label newElem = map[lst[i]];
  54
  55         if (newElem != -1)
  56         {
  57             dynLst.append(newElem);
  58         }
  59     }
  60 }
  61
  62
  63 // Append all mapped elements of a list to a DynamicList
  64 void Foam::polyMeshAdder::append
  65 (
  66     const labelList& map,
  67     const labelList& lst,
  68     const SortableList<label>& sortedLst,
  69     DynamicList<label>& dynLst
  70 )
  71 {
  72     dynLst.setSize(dynLst.size() + lst.size());
  73
  74     forAll(lst, i)
  75     {
  76         const label newElem = map[lst[i]];
  77
  78         if (newElem != -1 && findSortedIndex(sortedLst, newElem) == -1)
  79         {
  80             dynLst.append(newElem);
  81         }
  82     }
  83 }
  84
  85
  86 // Get index of patch in new set of patchnames/types
  87 Foam::label Foam::polyMeshAdder::patchIndex
  88 (
  89     const polyPatch& p,
  90     DynamicList<word>& allPatchNames,
  91     DynamicList<word>& allPatchTypes
  92 )
  93 {
  94     // Find the patch name on the list.  If the patch is already there
  95     // and patch types match, return index
  96     const word& pType = p.type();
  97     const word& pName = p.name();
  98
  99     label patchI = findIndex(allPatchNames, pName);
 100
 101     if (patchI == -1)
 102     {
 103         // Patch not found. Append to the list
 104         allPatchNames.append(pName);
 105         allPatchTypes.append(pType);
 106
 107         return allPatchNames.size() - 1;
 108     }
 109     else if (allPatchTypes[patchI] == pType)
 110     {
 111         // Found name and types match
 112         return patchI;
 113     }
 114     else
 115     {
 116         // Found the name, but type is different
 117
 118         // Duplicate name is not allowed.  Create a composite name from the
 119         // patch name and case name
 120         const word& caseName = p.boundaryMesh().mesh().time().caseName();
 121
 122         allPatchNames.append(pName + "_" + caseName);
 123         allPatchTypes.append(pType);
 124
 125         Pout<< "label patchIndex(const polyPatch& p) : "
 126             << "Patch " << p.index() << " named "
 127             << pName << " in mesh " << caseName
 128             << " already exists, but patch types"
 129             << " do not match.\nCreating a composite name as "
 130             << allPatchNames.last() << endl;
 131
 132         return allPatchNames.size() - 1;
 133     }
 134 }
 135
 136
 137 // Get index of zone in new set of zone names
 138 Foam::label Foam::polyMeshAdder::zoneIndex
 139 (
 140     const word& curName,
 141     DynamicList<word>& names
 142 )
 143 {
 144     label zoneI = findIndex(names, curName);
 145
 146     if (zoneI != -1)
 147     {
 148         return zoneI;
 149     }
 150     else
 151     {
 152         // Not found.  Add new name to the list
 153         names.append(curName);
 154
 155         return names.size() - 1;
 156     }
 157 }
 158
 159
 160 void Foam::polyMeshAdder::mergePatchNames
 161 (
 162     const polyBoundaryMesh& patches0,
 163     const polyBoundaryMesh& patches1,
 164
 165     DynamicList<word>& allPatchNames,
 166     DynamicList<word>& allPatchTypes,
 167
 168     labelList& from1ToAllPatches,
 169     labelList& fromAllTo1Patches
 170 )
 171 {
 172     // Insert the mesh0 patches and zones
 173     allPatchNames.append(patches0.names());
 174     allPatchTypes.append(patches0.types());
 175
 176
 177     // Patches
 178     // ~~~~~~~
 179     // Patches from 0 are taken over as is; those from 1 get either merged
 180     // (if they share name and type) or appended.
 181     // Empty patches are filtered out much much later on.
 182
 183     // Add mesh1 patches and build map both ways.
 184     from1ToAllPatches.setSize(patches1.size());
 185
 186     forAll(patches1, patchI)
 187     {
 188         from1ToAllPatches[patchI] = patchIndex
 189         (
 190             patches1[patchI],
 191             allPatchNames,
 192             allPatchTypes
 193         );
 194     }
 195     allPatchTypes.shrink();
 196     allPatchNames.shrink();
 197
 198     // Invert 1 to all patch map
 199     fromAllTo1Patches.setSize(allPatchNames.size());
 200     fromAllTo1Patches = -1;
 201
 202     forAll(from1ToAllPatches, i)
 203     {
 204         fromAllTo1Patches[from1ToAllPatches[i]] = i;
 205     }
 206 }
 207
 208
 209 Foam::labelList Foam::polyMeshAdder::getPatchStarts
 210 (
 211     const polyBoundaryMesh& patches
 212 )
 213 {
 214     labelList patchStarts(patches.size());
 215     forAll(patches, patchI)
 216     {
 217         patchStarts[patchI] = patches[patchI].start();
 218     }
 219     return patchStarts;
 220 }
 221
 222
 223 Foam::labelList Foam::polyMeshAdder::getPatchSizes
 224 (
 225     const polyBoundaryMesh& patches
 226 )
 227 {
 228     labelList patchSizes(patches.size());
 229     forAll(patches, patchI)
 230     {
 231         patchSizes[patchI] = patches[patchI].size();
 232     }
 233     return patchSizes;
 234 }
 235
 236
 237 Foam::List<Foam::polyPatch*> Foam::polyMeshAdder::combinePatches
 238 (
 239     const polyMesh& mesh0,
 240     const polyMesh& mesh1,
 241     const polyBoundaryMesh& allBoundaryMesh,
 242     const label nAllPatches,
 243     const labelList& fromAllTo1Patches,
 244
 245     const label nInternalFaces,
 246     const labelList& nFaces,
 247
 248     labelList& from0ToAllPatches,
 249     labelList& from1ToAllPatches
 250 )
 251 {
 252     const polyBoundaryMesh& patches0 = mesh0.boundaryMesh();
 253     const polyBoundaryMesh& patches1 = mesh1.boundaryMesh();
 254
 255
 256     // Compacted new patch list.
 257     DynamicList<polyPatch*> allPatches(nAllPatches);
 258
 259
 260     // Map from 0 to all patches (since gets compacted)
 261     from0ToAllPatches.setSize(patches0.size());
 262     from0ToAllPatches = -1;
 263
 264     label startFaceI = nInternalFaces;
 265
 266     // Copy patches0 with new sizes. First patches always come from
 267     // mesh0 and will always be present.
 268     forAll(patches0, patchI)
 269     {
 270         // Originates from mesh0. Clone with new size & filter out empty
 271         // patch.
 272         label filteredPatchI;
 273
 274         if (nFaces[patchI] == 0 && isA<processorPolyPatch>(patches0[patchI]))
 275         {
 276             //Pout<< "Removing zero sized mesh0 patch "
 277             //    << patches0[patchI].name() << endl;
 278             filteredPatchI = -1;
 279         }
 280         else
 281         {
 282             filteredPatchI = allPatches.size();
 283
 284             allPatches.append
 285             (
 286                 patches0[patchI].clone
 287                 (
 288                     allBoundaryMesh,
 289                     filteredPatchI,
 290                     nFaces[patchI],
 291                     startFaceI
 292                 ).ptr()
 293             );
 294             startFaceI += nFaces[patchI];
 295         }
 296
 297         // Record new index in allPatches
 298         from0ToAllPatches[patchI] = filteredPatchI;
 299
 300         // Check if patch was also in mesh1 and update its addressing if so.
 301         if (fromAllTo1Patches[patchI] != -1)
 302         {
 303             from1ToAllPatches[fromAllTo1Patches[patchI]] = filteredPatchI;
 304         }
 305     }
 306
 307     // Copy unique patches of mesh1.
 308     forAll(from1ToAllPatches, patchI)
 309     {
 310         label allPatchI = from1ToAllPatches[patchI];
 311
 312         if (allPatchI >= patches0.size())
 313         {
 314             // Patch has not been merged with any mesh0 patch.
 315
 316             label filteredPatchI;
 317
 318             if
 319             (
 320                 nFaces[allPatchI] == 0
 321              && isA<processorPolyPatch>(patches1[patchI])
 322             )
 323             {
 324                 //Pout<< "Removing zero sized mesh1 patch "
 325                 //    << patches1[patchI].name() << endl;
 326                 filteredPatchI = -1;
 327             }
 328             else
 329             {
 330                 filteredPatchI = allPatches.size();
 331
 332                 allPatches.append
 333                 (
 334                     patches1[patchI].clone
 335                     (
 336                         allBoundaryMesh,
 337                         filteredPatchI,
 338                         nFaces[allPatchI],
 339                         startFaceI
 340                     ).ptr()
 341                 );
 342                 startFaceI += nFaces[allPatchI];
 343             }
 344
 345             from1ToAllPatches[patchI] = filteredPatchI;
 346         }
 347     }
 348
 349     allPatches.shrink();
 350
 351     return allPatches;
 352 }
 353
 354
 355 Foam::labelList Foam::polyMeshAdder::getFaceOrder
 356 (
 357     const cellList& cells,
 358     const label nInternalFaces,
 359     const labelList& owner,
 360     const labelList& neighbour
 361 )
 362 {
 363     labelList oldToNew(owner.size(), -1);
 364
 365     // Leave boundary faces in order
 366     for (label faceI = nInternalFaces; faceI < owner.size(); ++faceI)
 367     {
 368         oldToNew[faceI] = faceI;
 369     }
 370
 371     // First unassigned face
 372     label newFaceI = 0;
 373
 374     forAll(cells, cellI)
 375     {
 376         const labelList& cFaces = cells[cellI];
 377
 378         SortableList<label> nbr(cFaces.size());
 379
 380         forAll(cFaces, i)
 381         {
 382             label faceI = cFaces[i];
 383
 384             label nbrCellI = neighbour[faceI];
 385
 386             if (nbrCellI != -1)
 387             {
 388                 // Internal face. Get cell on other side.
 389                 if (nbrCellI == cellI)
 390                 {
 391                     nbrCellI = owner[faceI];
 392                 }
 393
 394                 if (cellI < nbrCellI)
 395                 {
 396                     // CellI is master
 397                     nbr[i] = nbrCellI;
 398                 }
 399                 else
 400                 {
 401                     // nbrCell is master. Let it handle this face.
 402                     nbr[i] = -1;
 403                 }
 404             }
 405             else
 406             {
 407                 // External face. Do later.
 408                 nbr[i] = -1;
 409             }
 410         }
 411
 412         nbr.sort();
 413
 414         forAll(nbr, i)
 415         {
 416             if (nbr[i] != -1)
 417             {
 418                 oldToNew[cFaces[nbr.indices()[i]]] = newFaceI++;
 419             }
 420         }
 421     }
 422
 423
 424     // Check done all faces.
 425     forAll(oldToNew, faceI)
 426     {
 427         if (oldToNew[faceI] == -1)
 428         {
 429             FatalErrorIn
 430             (
 431                 "polyMeshAdder::getFaceOrder"
 432                 "(const cellList&, const label, const labelList&"
 433                 ", const labelList&) const"
 434             )   << "Did not determine new position"
 435                 << " for face " << faceI
 436                 << abort(FatalError);
 437         }
 438     }
 439
 440     return oldToNew;
 441 }
 442
 443
 444 // Extends face f with split points. cutEdgeToPoints gives for every
 445 // edge the points introduced inbetween the endpoints.
 446 void Foam::polyMeshAdder::insertVertices
 447 (
 448     const edgeLookup& cutEdgeToPoints,
 449     const Map<label>& meshToMaster,
 450     const labelList& masterToCutPoints,
 451     const face& masterF,
 452
 453     DynamicList<label>& workFace,
 454     face& allF
 455 )
 456 {
 457     workFace.clear();
 458
 459     // Check any edge for being cut (check on the cut so takes account
 460     // for any point merging on the cut)
 461
 462     forAll(masterF, fp)
 463     {
 464         label v0 = masterF[fp];
 465         label v1 = masterF.nextLabel(fp);
 466
 467         // Copy existing face point
 468         workFace.append(allF[fp]);
 469
 470         // See if any edge between v0,v1
 471
 472         Map<label>::const_iterator v0Fnd = meshToMaster.find(v0);
 473         if (v0Fnd != meshToMaster.end())
 474         {
 475             Map<label>::const_iterator v1Fnd = meshToMaster.find(v1);
 476             if (v1Fnd != meshToMaster.end())
 477             {
 478                 // Get edge in cutPoint numbering
 479                 edge cutEdge
 480                 (
 481                     masterToCutPoints[v0Fnd()],
 482                     masterToCutPoints[v1Fnd()]
 483                 );
 484
 485                 edgeLookup::const_iterator iter = cutEdgeToPoints.find(cutEdge);
 486
 487                 if (iter != cutEdgeToPoints.end())
 488                 {
 489                     const edge& e = iter.key();
 490                     const labelList& addedPoints = iter();
 491
 492                     // cutPoints first in allPoints so no need for renumbering
 493                     if (e[0] == cutEdge[0])
 494                     {
 495                         forAll(addedPoints, i)
 496                         {
 497                             workFace.append(addedPoints[i]);
 498                         }
 499                     }
 500                     else
 501                     {
 502                         forAllReverse(addedPoints, i)
 503                         {
 504                             workFace.append(addedPoints[i]);
 505                         }
 506                     }
 507                 }
 508             }
 509         }
 510     }
 511
 512     if (workFace.size() != allF.size())
 513     {
 514         allF.transfer(workFace);
 515     }
 516 }
 517
 518
 519 // Adds primitives (cells, faces, points)
 520 // Cells:
 521 //  - all of mesh0
 522 //  - all of mesh1
 523 // Faces:
 524 //  - all uncoupled of mesh0
 525 //  - all coupled faces
 526 //  - all uncoupled of mesh1
 527 // Points:
 528 //  - all coupled
 529 //  - all uncoupled of mesh0
 530 //  - all uncoupled of mesh1
 531 void Foam::polyMeshAdder::mergePrimitives
 532 (
 533     const polyMesh& mesh0,
 534     const polyMesh& mesh1,
 535     const faceCoupleInfo& coupleInfo,
 536
 537     const label nAllPatches,                // number of patches in the new mesh
 538     const labelList& fromAllTo1Patches,
 539     const labelList& from1ToAllPatches,
 540
 541     pointField& allPoints,
 542     labelList& from0ToAllPoints,
 543     labelList& from1ToAllPoints,
 544
 545     faceList& allFaces,
 546     labelList& allOwner,
 547     labelList& allNeighbour,
 548     label& nInternalFaces,
 549     labelList& nFacesPerPatch,
 550     label& nCells,
 551
 552     labelList& from0ToAllFaces,
 553     labelList& from1ToAllFaces,
 554     labelList& from1ToAllCells
 555 )
 556 {
 557     const polyBoundaryMesh& patches0 = mesh0.boundaryMesh();
 558     const polyBoundaryMesh& patches1 = mesh1.boundaryMesh();
 559
 560     const primitiveFacePatch& cutFaces = coupleInfo.cutFaces();
 561     const indirectPrimitivePatch& masterPatch = coupleInfo.masterPatch();
 562     const indirectPrimitivePatch& slavePatch = coupleInfo.slavePatch();
 563
 564
 565     // Points
 566     // ~~~~~~
 567
 568     // Storage for new points
 569     allPoints.setSize(mesh0.nPoints() + mesh1.nPoints());
 570     label allPointI = 0;
 571
 572     from0ToAllPoints.setSize(mesh0.nPoints());
 573     from0ToAllPoints = -1;
 574     from1ToAllPoints.setSize(mesh1.nPoints());
 575     from1ToAllPoints = -1;
 576
 577     // Copy coupled points (on cut)
 578     {
 579         const pointField& cutPoints = coupleInfo.cutPoints();
 580
 581         //const labelListList& cutToMasterPoints =
 582         //   coupleInfo.cutToMasterPoints();
 583         labelListList cutToMasterPoints
 584         (
 585             invertOneToMany
 586             (
 587                 cutPoints.size(),
 588                 coupleInfo.masterToCutPoints()
 589             )
 590         );
 591
 592         //const labelListList& cutToSlavePoints =
 593         //    coupleInfo.cutToSlavePoints();
 594         labelListList cutToSlavePoints
 595         (
 596             invertOneToMany
 597             (
 598                 cutPoints.size(),
 599                 coupleInfo.slaveToCutPoints()
 600             )
 601         );
 602
 603         forAll(cutPoints, i)
 604         {
 605             allPoints[allPointI] = cutPoints[i];
 606
 607             // Mark all master and slave points referring to this point.
 608
 609             const labelList& masterPoints = cutToMasterPoints[i];
 610
 611             forAll(masterPoints, j)
 612             {
 613                 label mesh0PointI = masterPatch.meshPoints()[masterPoints[j]];
 614                 from0ToAllPoints[mesh0PointI] = allPointI;
 615             }
 616
 617             const labelList& slavePoints = cutToSlavePoints[i];
 618
 619             forAll(slavePoints, j)
 620             {
 621                 label mesh1PointI = slavePatch.meshPoints()[slavePoints[j]];
 622                 from1ToAllPoints[mesh1PointI] = allPointI;
 623             }
 624             allPointI++;
 625         }
 626     }
 627
 628     // Add uncoupled mesh0 points
 629     forAll(mesh0.points(), pointI)
 630     {
 631         if (from0ToAllPoints[pointI] == -1)
 632         {
 633             allPoints[allPointI] = mesh0.points()[pointI];
 634             from0ToAllPoints[pointI] = allPointI;
 635             allPointI++;
 636         }
 637     }
 638
 639     // Add uncoupled mesh1 points
 640     forAll(mesh1.points(), pointI)
 641     {
 642         if (from1ToAllPoints[pointI] == -1)
 643         {
 644             allPoints[allPointI] = mesh1.points()[pointI];
 645             from1ToAllPoints[pointI] = allPointI;
 646             allPointI++;
 647         }
 648     }
 649
 650     allPoints.setSize(allPointI);
 651
 652
 653     // Faces
 654     // ~~~~~
 655
 656     // Sizes per patch
 657     nFacesPerPatch.setSize(nAllPatches);
 658     nFacesPerPatch = 0;
 659
 660     // Storage for faces and owner/neighbour
 661     allFaces.setSize(mesh0.nFaces() + mesh1.nFaces());
 662     allOwner.setSize(allFaces.size());
 663     allOwner = -1;
 664     allNeighbour.setSize(allFaces.size());
 665     allNeighbour = -1;
 666     label allFaceI = 0;
 667
 668     from0ToAllFaces.setSize(mesh0.nFaces());
 669     from0ToAllFaces = -1;
 670     from1ToAllFaces.setSize(mesh1.nFaces());
 671     from1ToAllFaces = -1;
 672
 673     // Copy mesh0 internal faces (always uncoupled)
 674     for (label faceI = 0; faceI < mesh0.nInternalFaces(); faceI++)
 675     {
 676         allFaces[allFaceI] = renumber(from0ToAllPoints, mesh0.faces()[faceI]);
 677         allOwner[allFaceI] = mesh0.faceOwner()[faceI];
 678         allNeighbour[allFaceI] = mesh0.faceNeighbour()[faceI];
 679         from0ToAllFaces[faceI] = allFaceI++;
 680     }
 681
 682     // Copy coupled faces. Every coupled face has an equivalent master and
 683     // slave. Also uncount as boundary faces all the newly coupled faces.
 684     const labelList& cutToMasterFaces = coupleInfo.cutToMasterFaces();
 685     const labelList& cutToSlaveFaces = coupleInfo.cutToSlaveFaces();
 686
 687     forAll(cutFaces, i)
 688     {
 689         label masterFaceI = cutToMasterFaces[i];
 690
 691         label mesh0FaceI = masterPatch.addressing()[masterFaceI];
 692
 693         if (from0ToAllFaces[mesh0FaceI] == -1)
 694         {
 695             // First occurrence of face
 696             from0ToAllFaces[mesh0FaceI] = allFaceI;
 697
 698             // External face becomes internal so uncount
 699             label patch0 = patches0.whichPatch(mesh0FaceI);
 700             nFacesPerPatch[patch0]--;
 701         }
 702
 703         label slaveFaceI = cutToSlaveFaces[i];
 704
 705         label mesh1FaceI = slavePatch.addressing()[slaveFaceI];
 706
 707         if (from1ToAllFaces[mesh1FaceI] == -1)
 708         {
 709             from1ToAllFaces[mesh1FaceI] = allFaceI;
 710
 711             label patch1 = patches1.whichPatch(mesh1FaceI);
 712             nFacesPerPatch[from1ToAllPatches[patch1]]--;
 713         }
 714
 715         // Copy cut face (since cutPoints are copied first no renumbering
 716         // nessecary)
 717         allFaces[allFaceI] = cutFaces[i];
 718         allOwner[allFaceI] = mesh0.faceOwner()[mesh0FaceI];
 719         allNeighbour[allFaceI] = mesh1.faceOwner()[mesh1FaceI] + mesh0.nCells();
 720
 721         allFaceI++;
 722     }
 723
 724     // Copy mesh1 internal faces (always uncoupled)
 725     for (label faceI = 0; faceI < mesh1.nInternalFaces(); faceI++)
 726     {
 727         allFaces[allFaceI] = renumber(from1ToAllPoints, mesh1.faces()[faceI]);
 728         allOwner[allFaceI] = mesh1.faceOwner()[faceI] + mesh0.nCells();
 729         allNeighbour[allFaceI] = mesh1.faceNeighbour()[faceI] + mesh0.nCells();
 730         from1ToAllFaces[faceI] = allFaceI++;
 731     }
 732
 733     nInternalFaces = allFaceI;
 734
 735     // Copy (unmarked/uncoupled) external faces in new order.
 736     for (label allPatchI = 0; allPatchI < nAllPatches; allPatchI++)
 737     {
 738         if (allPatchI < patches0.size())
 739         {
 740             // Patch is present in mesh0
 741             const polyPatch& pp = patches0[allPatchI];
 742
 743             nFacesPerPatch[allPatchI] += pp.size();
 744
 745             label faceI = pp.start();
 746
 747             forAll(pp, i)
 748             {
 749                 if (from0ToAllFaces[faceI] == -1)
 750                 {
 751                     // Is uncoupled face since has not yet been dealt with
 752                     allFaces[allFaceI] = renumber
 753                     (
 754                         from0ToAllPoints,
 755                         mesh0.faces()[faceI]
 756                     );
 757                     allOwner[allFaceI] = mesh0.faceOwner()[faceI];
 758                     allNeighbour[allFaceI] = -1;
 759
 760                     from0ToAllFaces[faceI] = allFaceI++;
 761                 }
 762                 faceI++;
 763             }
 764         }
 765         if (fromAllTo1Patches[allPatchI] != -1)
 766         {
 767             // Patch is present in mesh1
 768             const polyPatch& pp = patches1[fromAllTo1Patches[allPatchI]];
 769
 770             nFacesPerPatch[allPatchI] += pp.size();
 771
 772             label faceI = pp.start();
 773
 774             forAll(pp, i)
 775             {
 776                 if (from1ToAllFaces[faceI] == -1)
 777                 {
 778                     // Is uncoupled face
 779                     allFaces[allFaceI] = renumber
 780                     (
 781                         from1ToAllPoints,
 782                         mesh1.faces()[faceI]
 783                     );
 784                     allOwner[allFaceI] =
 785                         mesh1.faceOwner()[faceI]
 786                       + mesh0.nCells();
 787                     allNeighbour[allFaceI] = -1;
 788
 789                     from1ToAllFaces[faceI] = allFaceI++;
 790                 }
 791                 faceI++;
 792             }
 793         }
 794     }
 795     allFaces.setSize(allFaceI);
 796     allOwner.setSize(allFaceI);
 797     allNeighbour.setSize(allFaceI);
 798
 799
 800     // So now we have all ok for one-to-one mapping.
 801     // For split slace faces:
 802     // - mesh consistent with slave side
 803     // - mesh not consistent with owner side. It is not zipped up, the
 804     //   original faces need edges split.
 805
 806     // Use brute force to prevent having to calculate addressing:
 807     // - get map from master edge to split edges.
 808     // - check all faces to find any edge that is split.
 809     {
 810         // From two cut-points to labels of cut-points inbetween.
 811         // (in order: from e[0] to e[1]
 812         const edgeLookup& cutEdgeToPoints = coupleInfo.cutEdgeToPoints();
 813
 814         // Get map of master face (in mesh labels) that are in cut. These faces
 815         // do not need to be renumbered.
 816         labelHashSet masterCutFaces(cutToMasterFaces.size());
 817         forAll(cutToMasterFaces, i)
 818         {
 819             label meshFaceI = masterPatch.addressing()[cutToMasterFaces[i]];
 820
 821             masterCutFaces.insert(meshFaceI);
 822         }
 823
 824         DynamicList<label> workFace(100);
 825
 826         forAll(from0ToAllFaces, face0)
 827         {
 828             if (!masterCutFaces.found(face0))
 829             {
 830                 label allFaceI = from0ToAllFaces[face0];
 831
 832                 insertVertices
 833                 (
 834                     cutEdgeToPoints,
 835                     masterPatch.meshPointMap(),
 836                     coupleInfo.masterToCutPoints(),
 837                     mesh0.faces()[face0],
 838
 839                     workFace,
 840                     allFaces[allFaceI]
 841                 );
 842             }
 843         }
 844
 845         // Same for slave face
 846
 847         labelHashSet slaveCutFaces(cutToSlaveFaces.size());
 848         forAll(cutToSlaveFaces, i)
 849         {
 850             label meshFaceI = slavePatch.addressing()[cutToSlaveFaces[i]];
 851
 852             slaveCutFaces.insert(meshFaceI);
 853         }
 854
 855         forAll(from1ToAllFaces, face1)
 856         {
 857             if (!slaveCutFaces.found(face1))
 858             {
 859                 label allFaceI = from1ToAllFaces[face1];
 860
 861                 insertVertices
 862                 (
 863                     cutEdgeToPoints,
 864                     slavePatch.meshPointMap(),
 865                     coupleInfo.slaveToCutPoints(),
 866                     mesh1.faces()[face1],
 867
 868                     workFace,
 869                     allFaces[allFaceI]
 870                 );
 871             }
 872         }
 873     }
 874
 875     // Now we have a full facelist and owner/neighbour addressing.
 876
 877
 878     // Cells
 879     // ~~~~~
 880
 881     from1ToAllCells.setSize(mesh1.nCells());
 882     from1ToAllCells = -1;
 883
 884     forAll(mesh1.cells(), i)
 885     {
 886         from1ToAllCells[i] = i + mesh0.nCells();
 887     }
 888
 889     // Make cells (= cell-face addressing)
 890     nCells = mesh0.nCells() + mesh1.nCells();
 891     cellList allCells(nCells);
 892     primitiveMesh::calcCells(allCells, allOwner, allNeighbour, nCells);
 893
 894     // Reorder faces for upper-triangular order.
 895     labelList oldToNew
 896     (
 897         getFaceOrder
 898         (
 899             allCells,
 900             nInternalFaces,
 901             allOwner,
 902             allNeighbour
 903         )
 904     );
 905
 906     inplaceReorder(oldToNew, allFaces);
 907     inplaceReorder(oldToNew, allOwner);
 908     inplaceReorder(oldToNew, allNeighbour);
 909     inplaceRenumber(oldToNew, from0ToAllFaces);
 910     inplaceRenumber(oldToNew, from1ToAllFaces);
 911 }
 912
 913
 914 void Foam::polyMeshAdder::mergePointZones
 915 (
 916     const pointZoneMesh& pz0,
 917     const pointZoneMesh& pz1,
 918     const labelList& from0ToAllPoints,
 919     const labelList& from1ToAllPoints,
 920
 921     DynamicList<word>& zoneNames,
 922     labelList& from1ToAll,
 923     List<DynamicList<label> >& pzPoints
 924 )
 925 {
 926     zoneNames.setCapacity(pz0.size() + pz1.size());
 927     zoneNames.append(pz0.names());
 928
 929     from1ToAll.setSize(pz1.size());
 930
 931     forAll(pz1, zoneI)
 932     {
 933         from1ToAll[zoneI] = zoneIndex(pz1[zoneI].name(), zoneNames);
 934     }
 935     zoneNames.shrink();
 936
 937     // Point labels per merged zone
 938     pzPoints.setSize(zoneNames.size());
 939
 940     forAll(pz0, zoneI)
 941     {
 942         append(from0ToAllPoints, pz0[zoneI], pzPoints[zoneI]);
 943     }
 944
 945     // Get sorted zone contents for duplicate element recognition
 946     PtrList<SortableList<label> > pzPointsSorted(pzPoints.size());
 947     forAll(pzPoints, zoneI)
 948     {
 949         pzPointsSorted.set
 950         (
 951             zoneI,
 952             new SortableList<label>(pzPoints[zoneI])
 953         );
 954     }
 955
 956     // Now we have full addressing for points so do the pointZones of mesh1.
 957     forAll(pz1, zoneI)
 958     {
 959         // Relabel all points of zone and add to correct pzPoints.
 960         const label allZoneI = from1ToAll[zoneI];
 961
 962         append
 963         (
 964             from1ToAllPoints,
 965             pz1[zoneI],
 966             pzPointsSorted[allZoneI],
 967             pzPoints[allZoneI]
 968         );
 969     }
 970
 971     forAll(pzPoints, i)
 972     {
 973         pzPoints[i].shrink();
 974     }
 975 }
 976
 977
 978 void Foam::polyMeshAdder::mergeFaceZones
 979 (
 980     const faceZoneMesh& fz0,
 981     const faceZoneMesh& fz1,
 982     const labelList& from0ToAllFaces,
 983     const labelList& from1ToAllFaces,
 984
 985     DynamicList<word>& zoneNames,
 986     labelList& from1ToAll,
 987     List<DynamicList<label> >& fzFaces,
 988     List<DynamicList<bool> >& fzFlips
 989 )
 990 {
 991     zoneNames.setCapacity(fz0.size() + fz1.size());
 992     zoneNames.append(fz0.names());
 993
 994     from1ToAll.setSize(fz1.size());
 995
 996     forAll(fz1, zoneI)
 997     {
 998         from1ToAll[zoneI] = zoneIndex(fz1[zoneI].name(), zoneNames);
 999     }
1000     zoneNames.shrink();
1001
1002
1003     // Create temporary lists for faceZones.
1004     fzFaces.setSize(zoneNames.size());
1005     fzFlips.setSize(zoneNames.size());
1006     forAll(fz0, zoneI)
1007     {
1008         DynamicList<label>& newZone = fzFaces[zoneI];
1009         DynamicList<bool>& newFlip = fzFlips[zoneI];
1010
1011         newZone.setCapacity(fz0[zoneI].size());
1012         newFlip.setCapacity(newZone.size());
1013
1014         const labelList& addressing = fz0[zoneI];
1015         const boolList& flipMap = fz0[zoneI].flipMap();
1016
1017         forAll(addressing, i)
1018         {
1019             label faceI = addressing[i];
1020
1021             if (from0ToAllFaces[faceI] != -1)
1022             {
1023                 newZone.append(from0ToAllFaces[faceI]);
1024                 newFlip.append(flipMap[i]);
1025             }
1026         }
1027     }
1028
1029     // Get sorted zone contents for duplicate element recognition
1030     PtrList<SortableList<label> > fzFacesSorted(fzFaces.size());
1031     forAll(fzFaces, zoneI)
1032     {
1033         fzFacesSorted.set
1034         (
1035             zoneI,
1036             new SortableList<label>(fzFaces[zoneI])
1037         );
1038     }
1039
1040     // Now we have full addressing for faces so do the faceZones of mesh1.
1041     forAll(fz1, zoneI)
1042     {
1043         label allZoneI = from1ToAll[zoneI];
1044
1045         DynamicList<label>& newZone = fzFaces[allZoneI];
1046         const SortableList<label>& newZoneSorted = fzFacesSorted[allZoneI];
1047         DynamicList<bool>& newFlip = fzFlips[allZoneI];
1048
1049         newZone.setCapacity(newZone.size() + fz1[zoneI].size());
1050         newFlip.setCapacity(newZone.size());
1051
1052         const labelList& addressing = fz1[zoneI];
1053         const boolList& flipMap = fz1[zoneI].flipMap();
1054
1055         forAll(addressing, i)
1056         {
1057             label faceI = addressing[i];
1058             label allFaceI = from1ToAllFaces[faceI];
1059
1060             if
1061             (
1062                 allFaceI != -1
1063              && findSortedIndex(newZoneSorted, allFaceI) == -1
1064             )
1065             {
1066                 newZone.append(allFaceI);
1067                 newFlip.append(flipMap[i]);
1068             }
1069         }
1070     }
1071
1072     forAll(fzFaces, i)
1073     {
1074         fzFaces[i].shrink();
1075         fzFlips[i].shrink();
1076     }
1077 }
1078
1079
1080 void Foam::polyMeshAdder::mergeCellZones
1081 (
1082     const cellZoneMesh& cz0,
1083     const cellZoneMesh& cz1,
1084     const labelList& from1ToAllCells,
1085
1086     DynamicList<word>& zoneNames,
1087     labelList& from1ToAll,
1088     List<DynamicList<label> >& czCells
1089 )
1090 {
1091     zoneNames.setCapacity(cz0.size() + cz1.size());
1092     zoneNames.append(cz0.names());
1093
1094     from1ToAll.setSize(cz1.size());
1095     forAll(cz1, zoneI)
1096     {
1097         from1ToAll[zoneI] = zoneIndex(cz1[zoneI].name(), zoneNames);
1098     }
1099     zoneNames.shrink();
1100
1101
1102     // Create temporary lists for cellZones.
1103     czCells.setSize(zoneNames.size());
1104     forAll(cz0, zoneI)
1105     {
1106         // Insert mesh0 cells
1107         czCells[zoneI].append(cz0[zoneI]);
1108     }
1109
1110
1111     // Cell mapping is trivial.
1112     forAll(cz1, zoneI)
1113     {
1114         const label allZoneI = from1ToAll[zoneI];
1115
1116         append(from1ToAllCells, cz1[zoneI], czCells[allZoneI]);
1117     }
1118
1119     forAll(czCells, i)
1120     {
1121         czCells[i].shrink();
1122     }
1123 }
1124
1125
1126 void Foam::polyMeshAdder::mergeZones
1127 (
1128     const polyMesh& mesh0,
1129     const polyMesh& mesh1,
1130     const labelList& from0ToAllPoints,
1131     const labelList& from0ToAllFaces,
1132
1133     const labelList& from1ToAllPoints,
1134     const labelList& from1ToAllFaces,
1135     const labelList& from1ToAllCells,
1136
1137     DynamicList<word>& pointZoneNames,
1138     List<DynamicList<label> >& pzPoints,
1139
1140     DynamicList<word>& faceZoneNames,
1141     List<DynamicList<label> >& fzFaces,
1142     List<DynamicList<bool> >& fzFlips,
1143
1144     DynamicList<word>& cellZoneNames,
1145     List<DynamicList<label> >& czCells
1146 )
1147 {
1148     labelList from1ToAllPZones;
1149     mergePointZones
1150     (
1151         mesh0.pointZones(),
1152         mesh1.pointZones(),
1153         from0ToAllPoints,
1154         from1ToAllPoints,
1155
1156         pointZoneNames,
1157         from1ToAllPZones,
1158         pzPoints
1159     );
1160
1161     labelList from1ToAllFZones;
1162     mergeFaceZones
1163     (
1164         mesh0.faceZones(),
1165         mesh1.faceZones(),
1166         from0ToAllFaces,
1167         from1ToAllFaces,
1168
1169         faceZoneNames,
1170         from1ToAllFZones,
1171         fzFaces,
1172         fzFlips
1173     );
1174
1175     labelList from1ToAllCZones;
1176     mergeCellZones
1177     (
1178         mesh0.cellZones(),
1179         mesh1.cellZones(),
1180         from1ToAllCells,
1181
1182         cellZoneNames,
1183         from1ToAllCZones,
1184         czCells
1185     );
1186 }
1187
1188
1189 void Foam::polyMeshAdder::addZones
1190 (
1191     const DynamicList<word>& pointZoneNames,
1192     const List<DynamicList<label> >& pzPoints,
1193
1194     const DynamicList<word>& faceZoneNames,
1195     const List<DynamicList<label> >& fzFaces,
1196     const List<DynamicList<bool> >& fzFlips,
1197
1198     const DynamicList<word>& cellZoneNames,
1199     const List<DynamicList<label> >& czCells,
1200
1201     polyMesh& mesh
1202 )
1203 {
1204     List<pointZone*> pZones(pzPoints.size());
1205     forAll(pZones, i)
1206     {
1207         pZones[i] = new pointZone
1208         (
1209             pointZoneNames[i],
1210             pzPoints[i],
1211             i,
1212             mesh.pointZones()
1213         );
1214     }
1215
1216     List<faceZone*> fZones(fzFaces.size());
1217     forAll(fZones, i)
1218     {
1219         fZones[i] = new faceZone
1220         (
1221             faceZoneNames[i],
1222             fzFaces[i],
1223             fzFlips[i],
1224             i,
1225             mesh.faceZones()
1226         );
1227     }
1228
1229     List<cellZone*> cZones(czCells.size());
1230     forAll(cZones, i)
1231     {
1232         cZones[i] = new cellZone
1233         (
1234             cellZoneNames[i],
1235             czCells[i],
1236             i,
1237             mesh.cellZones()
1238         );
1239     }
1240
1241     mesh.addZones
1242     (
1243         pZones,
1244         fZones,
1245         cZones
1246     );
1247 }
1248
1249
1250
1251 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
1252
1253 // Returns new mesh and sets
1254 // - map from new cell/face/point/patch to either mesh0 or mesh1
1255 //
1256 // mesh0Faces/mesh1Faces: corresponding faces on both meshes.
1257 Foam::autoPtr<Foam::polyMesh> Foam::polyMeshAdder::add
1258 (
1259     const IOobject& io,
1260     const polyMesh& mesh0,
1261     const polyMesh& mesh1,
1262     const faceCoupleInfo& coupleInfo,
1263     autoPtr<mapAddedPolyMesh>& mapPtr
1264 )
1265 {
1266     const polyBoundaryMesh& patches0 = mesh0.boundaryMesh();
1267     const polyBoundaryMesh& patches1 = mesh1.boundaryMesh();
1268
1269
1270     DynamicList<word> allPatchNames(patches0.size() + patches1.size());
1271     DynamicList<word> allPatchTypes(allPatchNames.size());
1272
1273     // Patch maps
1274     labelList from1ToAllPatches(patches1.size());
1275     labelList fromAllTo1Patches(allPatchNames.size(), -1);
1276
1277     mergePatchNames
1278     (
1279         patches0,
1280         patches1,
1281         allPatchNames,
1282         allPatchTypes,
1283         from1ToAllPatches,
1284         fromAllTo1Patches
1285     );
1286
1287
1288     // New points
1289     pointField allPoints;
1290
1291     // Map from mesh0/1 points to allPoints.
1292     labelList from0ToAllPoints(mesh0.nPoints(), -1);
1293     labelList from1ToAllPoints(mesh1.nPoints(), -1);
1294
1295     // New faces
1296     faceList allFaces;
1297     label nInternalFaces;
1298
1299     // New cells
1300     labelList allOwner;
1301     labelList allNeighbour;
1302     label nCells;
1303
1304     // Sizes per patch
1305     labelList nFaces(allPatchNames.size(), 0);
1306
1307     // Maps
1308     labelList from0ToAllFaces(mesh0.nFaces(), -1);
1309     labelList from1ToAllFaces(mesh1.nFaces(), -1);
1310
1311     // Map
1312     labelList from1ToAllCells(mesh1.nCells(), -1);
1313
1314     mergePrimitives
1315     (
1316         mesh0,
1317         mesh1,
1318         coupleInfo,
1319
1320         allPatchNames.size(),
1321         fromAllTo1Patches,
1322         from1ToAllPatches,
1323
1324         allPoints,
1325         from0ToAllPoints,
1326         from1ToAllPoints,
1327
1328         allFaces,
1329         allOwner,
1330         allNeighbour,
1331         nInternalFaces,
1332         nFaces,
1333         nCells,
1334
1335         from0ToAllFaces,
1336         from1ToAllFaces,
1337         from1ToAllCells
1338     );
1339
1340
1341     // Zones
1342     // ~~~~~
1343
1344     DynamicList<word> pointZoneNames;
1345     List<DynamicList<label> > pzPoints;
1346
1347     DynamicList<word> faceZoneNames;
1348     List<DynamicList<label> > fzFaces;
1349     List<DynamicList<bool> > fzFlips;
1350
1351     DynamicList<word> cellZoneNames;
1352     List<DynamicList<label> > czCells;
1353
1354     mergeZones
1355     (
1356         mesh0,
1357         mesh1,
1358
1359         from0ToAllPoints,
1360         from0ToAllFaces,
1361
1362         from1ToAllPoints,
1363         from1ToAllFaces,
1364         from1ToAllCells,
1365
1366         pointZoneNames,
1367         pzPoints,
1368
1369         faceZoneNames,
1370         fzFaces,
1371         fzFlips,
1372
1373         cellZoneNames,
1374         czCells
1375     );
1376
1377
1378     // Patches
1379     // ~~~~~~~
1380
1381     // Map from 0 to all patches (since gets compacted)
1382     labelList from0ToAllPatches(patches0.size(), -1);
1383
1384     List<polyPatch*> allPatches
1385     (
1386         combinePatches
1387         (
1388             mesh0,
1389             mesh1,
1390             patches0,           // Should be boundaryMesh() on new mesh.
1391             allPatchNames.size(),
1392             fromAllTo1Patches,
1393             mesh0.nInternalFaces()
1394           + mesh1.nInternalFaces()
1395           + coupleInfo.cutFaces().size(),
1396             nFaces,
1397
1398             from0ToAllPatches,
1399             from1ToAllPatches
1400         )
1401     );
1402
1403
1404     // Map information
1405     // ~~~~~~~~~~~~~~~
1406
1407     mapPtr.reset
1408     (
1409         new mapAddedPolyMesh
1410         (
1411             mesh0.nPoints(),
1412             mesh0.nFaces(),
1413             mesh0.nCells(),
1414
1415             mesh1.nPoints(),
1416             mesh1.nFaces(),
1417             mesh1.nCells(),
1418
1419             from0ToAllPoints,
1420             from0ToAllFaces,
1421             identity(mesh0.nCells()),
1422
1423             from1ToAllPoints,
1424             from1ToAllFaces,
1425             from1ToAllCells,
1426
1427             from0ToAllPatches,
1428             from1ToAllPatches,
1429             getPatchSizes(patches0),
1430             getPatchStarts(patches0)
1431         )
1432     );
1433
1434
1435
1436     // Now we have extracted all information from all meshes.
1437     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1438
1439     // Construct mesh
1440     autoPtr<polyMesh> tmesh
1441     (
1442         new polyMesh
1443         (
1444             io,
1445             xferMove(allPoints),
1446             xferMove(allFaces),
1447             xferMove(allOwner),
1448             xferMove(allNeighbour)
1449         )
1450     );
1451     polyMesh& mesh = tmesh();
1452
1453     // Add zones to new mesh.
1454     addZones
1455     (
1456         pointZoneNames,
1457         pzPoints,
1458
1459         faceZoneNames,
1460         fzFaces,
1461         fzFlips,
1462
1463         cellZoneNames,
1464         czCells,
1465         mesh
1466     );
1467
1468     // Add patches to new mesh
1469     mesh.addPatches(allPatches);
1470
1471     return tmesh;
1472 }
1473
1474
1475 // Inplace add mesh1 to mesh0
1476 Foam::autoPtr<Foam::mapAddedPolyMesh> Foam::polyMeshAdder::add
1477 (
1478     polyMesh& mesh0,
1479     const polyMesh& mesh1,
1480     const faceCoupleInfo& coupleInfo,
1481     const bool validBoundary
1482 )
1483 {
1484     const polyBoundaryMesh& patches0 = mesh0.boundaryMesh();
1485     const polyBoundaryMesh& patches1 = mesh1.boundaryMesh();
1486
1487     DynamicList<word> allPatchNames(patches0.size() + patches1.size());
1488     DynamicList<word> allPatchTypes(allPatchNames.size());
1489
1490     // Patch maps
1491     labelList from1ToAllPatches(patches1.size());
1492     labelList fromAllTo1Patches(allPatchNames.size(), -1);
1493
1494     mergePatchNames
1495     (
1496         patches0,
1497         patches1,
1498         allPatchNames,
1499         allPatchTypes,
1500         from1ToAllPatches,
1501         fromAllTo1Patches
1502     );
1503
1504
1505     // New points
1506     pointField allPoints;
1507
1508     // Map from mesh0/1 points to allPoints.
1509     labelList from0ToAllPoints(mesh0.nPoints(), -1);
1510     labelList from1ToAllPoints(mesh1.nPoints(), -1);
1511
1512     // New faces
1513     faceList allFaces;
1514     labelList allOwner;
1515     labelList allNeighbour;
1516     label nInternalFaces;
1517     // Sizes per patch
1518     labelList nFaces(allPatchNames.size(), 0);
1519     label nCells;
1520
1521     // Maps
1522     labelList from0ToAllFaces(mesh0.nFaces(), -1);
1523     labelList from1ToAllFaces(mesh1.nFaces(), -1);
1524     // Map
1525     labelList from1ToAllCells(mesh1.nCells(), -1);
1526
1527     mergePrimitives
1528     (
1529         mesh0,
1530         mesh1,
1531         coupleInfo,
1532
1533         allPatchNames.size(),
1534         fromAllTo1Patches,
1535         from1ToAllPatches,
1536
1537         allPoints,
1538         from0ToAllPoints,
1539         from1ToAllPoints,
1540
1541         allFaces,
1542         allOwner,
1543         allNeighbour,
1544         nInternalFaces,
1545         nFaces,
1546         nCells,
1547
1548         from0ToAllFaces,
1549         from1ToAllFaces,
1550         from1ToAllCells
1551     );
1552
1553
1554     // Zones
1555     // ~~~~~
1556
1557     DynamicList<word> pointZoneNames;
1558     List<DynamicList<label> > pzPoints;
1559
1560     DynamicList<word> faceZoneNames;
1561     List<DynamicList<label> > fzFaces;
1562     List<DynamicList<bool> > fzFlips;
1563
1564     DynamicList<word> cellZoneNames;
1565     List<DynamicList<label> > czCells;
1566
1567     mergeZones
1568     (
1569         mesh0,
1570         mesh1,
1571
1572         from0ToAllPoints,
1573         from0ToAllFaces,
1574
1575         from1ToAllPoints,
1576         from1ToAllFaces,
1577         from1ToAllCells,
1578
1579         pointZoneNames,
1580         pzPoints,
1581
1582         faceZoneNames,
1583         fzFaces,
1584         fzFlips,
1585
1586         cellZoneNames,
1587         czCells
1588     );
1589
1590
1591     // Patches
1592     // ~~~~~~~
1593
1594
1595     // Store mesh0 patch info before modifying patches0.
1596     labelList mesh0PatchSizes(getPatchSizes(patches0));
1597     labelList mesh0PatchStarts(getPatchStarts(patches0));
1598
1599     // Map from 0 to all patches (since gets compacted)
1600     labelList from0ToAllPatches(patches0.size(), -1);
1601
1602     // Inplace extend mesh0 patches (note that patches0.size() now also
1603     // has changed)
1604     polyBoundaryMesh& allPatches =
1605         const_cast<polyBoundaryMesh&>(mesh0.boundaryMesh());
1606     allPatches.setSize(allPatchNames.size());
1607
1608     label startFaceI = nInternalFaces;
1609
1610     // Copy patches0 with new sizes. First patches always come from
1611     // mesh0 and will always be present.
1612     label allPatchI = 0;
1613
1614     forAll(from0ToAllPatches, patch0)
1615     {
1616         // Originates from mesh0. Clone with new size & filter out empty
1617         // patch.
1618
1619         if (nFaces[patch0] == 0 && isA<processorPolyPatch>(allPatches[patch0]))
1620         {
1621             //Pout<< "Removing zero sized mesh0 patch " << allPatchNames[patch0]
1622             //    << endl;
1623             from0ToAllPatches[patch0] = -1;
1624             // Check if patch was also in mesh1 and update its addressing if so.
1625             if (fromAllTo1Patches[patch0] != -1)
1626             {
1627                 from1ToAllPatches[fromAllTo1Patches[patch0]] = -1;
1628             }
1629         }
1630         else
1631         {
1632             // Clone.
1633             allPatches.set
1634             (
1635                 allPatchI,
1636                 allPatches[patch0].clone
1637                 (
1638                     allPatches,
1639                     allPatchI,
1640                     nFaces[patch0],
1641                     startFaceI
1642                 )
1643             );
1644
1645             // Record new index in allPatches
1646             from0ToAllPatches[patch0] = allPatchI;
1647
1648             // Check if patch was also in mesh1 and update its addressing if so.
1649             if (fromAllTo1Patches[patch0] != -1)
1650             {
1651                 from1ToAllPatches[fromAllTo1Patches[patch0]] = allPatchI;
1652             }
1653
1654             startFaceI += nFaces[patch0];
1655
1656             allPatchI++;
1657         }
1658     }
1659
1660     // Copy unique patches of mesh1.
1661     forAll(from1ToAllPatches, patch1)
1662     {
1663         label uncompactAllPatchI = from1ToAllPatches[patch1];
1664
1665         if (uncompactAllPatchI >= from0ToAllPatches.size())
1666         {
1667             // Patch has not been merged with any mesh0 patch.
1668
1669             if
1670             (
1671                 nFaces[uncompactAllPatchI] == 0
1672              && isA<processorPolyPatch>(patches1[patch1])
1673             )
1674             {
1675                 //Pout<< "Removing zero sized mesh1 patch "
1676                 //    << allPatchNames[uncompactAllPatchI] << endl;
1677                 from1ToAllPatches[patch1] = -1;
1678             }
1679             else
1680             {
1681                 // Clone.
1682                 allPatches.set
1683                 (
1684                     allPatchI,
1685                     patches1[patch1].clone
1686                     (
1687                         allPatches,
1688                         allPatchI,
1689                         nFaces[uncompactAllPatchI],
1690                         startFaceI
1691                     )
1692                 );
1693
1694                 // Record new index in allPatches
1695                 from1ToAllPatches[patch1] = allPatchI;
1696
1697                 startFaceI += nFaces[uncompactAllPatchI];
1698
1699                 allPatchI++;
1700             }
1701         }
1702     }
1703
1704     allPatches.setSize(allPatchI);
1705
1706
1707     // Construct map information before changing mesh0 primitives
1708     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1709
1710     autoPtr<mapAddedPolyMesh> mapPtr
1711     (
1712         new mapAddedPolyMesh
1713         (
1714             mesh0.nPoints(),
1715             mesh0.nFaces(),
1716             mesh0.nCells(),
1717
1718             mesh1.nPoints(),
1719             mesh1.nFaces(),
1720             mesh1.nCells(),
1721
1722             from0ToAllPoints,
1723             from0ToAllFaces,
1724             identity(mesh0.nCells()),
1725
1726             from1ToAllPoints,
1727             from1ToAllFaces,
1728             from1ToAllCells,
1729
1730             from0ToAllPatches,
1731             from1ToAllPatches,
1732
1733             mesh0PatchSizes,
1734             mesh0PatchStarts
1735         )
1736     );
1737
1738
1739
1740     // Now we have extracted all information from all meshes.
1741     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1742
1743     labelList patchSizes(getPatchSizes(allPatches));
1744     labelList patchStarts(getPatchStarts(allPatches));
1745
1746     mesh0.resetMotion();    // delete any oldPoints.
1747     mesh0.resetPrimitives
1748     (
1749         xferMove(allPoints),
1750         xferMove(allFaces),
1751         xferMove(allOwner),
1752         xferMove(allNeighbour),
1753         patchSizes,     // size
1754         patchStarts,    // patchstarts
1755         validBoundary   // boundary valid?
1756     );
1757
1758     // Add zones to new mesh.
1759     mesh0.pointZones().clear();
1760     mesh0.faceZones().clear();
1761     mesh0.cellZones().clear();
1762     addZones
1763     (
1764         pointZoneNames,
1765         pzPoints,
1766
1767         faceZoneNames,
1768         fzFaces,
1769         fzFlips,
1770
1771         cellZoneNames,
1772         czCells,
1773         mesh0
1774     );
1775
1776     return mapPtr;
1777 }
1778
1779
1780 Foam::Map<Foam::label> Foam::polyMeshAdder::findSharedPoints
1781 (
1782     const polyMesh& mesh,
1783     const scalar mergeDist
1784 )
1785 {
1786     const labelList& sharedPointLabels = mesh.globalData().sharedPointLabels();
1787     const labelList& sharedPointAddr = mesh.globalData().sharedPointAddr();
1788
1789     // Because of adding the missing pieces e.g. when redistributing a mesh
1790     // it can be that there are multiple points on the same processor that
1791     // refer to the same shared point.
1792
1793     // Invert point-to-shared addressing
1794     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1795
1796     Map<labelList> sharedToMesh(sharedPointLabels.size());
1797
1798     label nMultiple = 0;
1799
1800     forAll(sharedPointLabels, i)
1801     {
1802         label pointI = sharedPointLabels[i];
1803
1804         label sharedI = sharedPointAddr[i];
1805
1806         Map<labelList>::iterator iter = sharedToMesh.find(sharedI);
1807
1808         if (iter != sharedToMesh.end())
1809         {
1810             // sharedI already used by other point. Add this one.
1811
1812             nMultiple++;
1813
1814             labelList& connectedPointLabels = iter();
1815
1816             label sz = connectedPointLabels.size();
1817
1818             // Check just to make sure.
1819             if (findIndex(connectedPointLabels, pointI) != -1)
1820             {
1821                 FatalErrorIn("polyMeshAdder::findSharedPoints(..)")
1822                     << "Duplicate point in sharedPoint addressing." << endl
1823                     << "When trying to add point " << pointI << " on shared "
1824                     << sharedI  << " with connected points "
1825                     << connectedPointLabels
1826                     << abort(FatalError);
1827             }
1828
1829             connectedPointLabels.setSize(sz+1);
1830             connectedPointLabels[sz] = pointI;
1831         }
1832         else
1833         {
1834             sharedToMesh.insert(sharedI, labelList(1, pointI));
1835         }
1836     }
1837
1838
1839     // Assign single master for every shared with multiple geometric points
1840     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1841
1842     Map<label> pointToMaster(nMultiple);
1843
1844     forAllConstIter(Map<labelList>, sharedToMesh, iter)
1845     {
1846         const labelList& connectedPointLabels = iter();
1847
1848         //Pout<< "For shared:" << iter.key()
1849         //    << " found points:" << connectedPointLabels
1850         //    << " at coords:"
1851         //    <<  pointField(mesh.points(), connectedPointLabels) << endl;
1852
1853         if (connectedPointLabels.size() > 1)
1854         {
1855             const pointField connectedPoints
1856             (
1857                 mesh.points(),
1858                 connectedPointLabels
1859             );
1860
1861             labelList toMergedPoints;
1862             pointField mergedPoints;
1863             bool hasMerged = Foam::mergePoints
1864             (
1865                 connectedPoints,
1866                 mergeDist,
1867                 false,
1868                 toMergedPoints,
1869                 mergedPoints
1870             );
1871
1872             if (hasMerged)
1873             {
1874                 // Invert toMergedPoints
1875                 const labelListList mergeSets
1876                 (
1877                     invertOneToMany
1878                     (
1879                         mergedPoints.size(),
1880                         toMergedPoints
1881                     )
1882                 );
1883
1884                 // Find master for valid merges
1885                 forAll(mergeSets, setI)
1886                 {
1887                     const labelList& mergeSet = mergeSets[setI];
1888
1889                     if (mergeSet.size() > 1)
1890                     {
1891                         // Pick lowest numbered point
1892                         label masterPointI = labelMax;
1893
1894                         forAll(mergeSet, i)
1895                         {
1896                             label pointI = connectedPointLabels[mergeSet[i]];
1897
1898                             masterPointI = min(masterPointI, pointI);
1899                         }
1900
1901                         forAll(mergeSet, i)
1902                         {
1903                             label pointI = connectedPointLabels[mergeSet[i]];
1904
1905                             //Pout<< "Merging point " << pointI
1906                             //    << " at " << mesh.points()[pointI]
1907                             //    << " into master point "
1908                             //    << masterPointI
1909                             //    << " at " << mesh.points()[masterPointI]
1910                             //    << endl;
1911
1912                             pointToMaster.insert(pointI, masterPointI);
1913                         }
1914                     }
1915                 }
1916             }
1917         }
1918     }
1919
1920     //- Old: geometric merging. Causes problems for two close shared points.
1921     //labelList sharedToMerged;
1922     //pointField mergedPoints;
1923     //bool hasMerged = Foam::mergePoints
1924     //(
1925     //    pointField
1926     //    (
1927     //        mesh.points(),
1928     //        sharedPointLabels
1929     //    ),
1930     //    mergeDist,
1931     //    false,
1932     //    sharedToMerged,
1933     //    mergedPoints
1934     //);
1935     //
1936     //// Find out which sets of points get merged and create a map from
1937     //// mesh point to unique point.
1938     //
1939     //Map<label> pointToMaster(10*sharedToMerged.size());
1940     //
1941     //if (hasMerged)
1942     //{
1943     //    labelListList mergeSets
1944     //    (
1945     //        invertOneToMany
1946     //        (
1947     //            sharedToMerged.size(),
1948     //            sharedToMerged
1949     //        )
1950     //    );
1951     //
1952     //    label nMergeSets = 0;
1953     //
1954     //    forAll(mergeSets, setI)
1955     //    {
1956     //        const labelList& mergeSet = mergeSets[setI];
1957     //
1958     //        if (mergeSet.size() > 1)
1959     //        {
1960     //            // Take as master the shared point with the lowest mesh
1961     //            // point label. (rather arbitrarily - could use max or
1962     //            // any other one of the points)
1963     //
1964     //            nMergeSets++;
1965     //
1966     //            label masterI = labelMax;
1967     //
1968     //            forAll(mergeSet, i)
1969     //            {
1970     //                label sharedI = mergeSet[i];
1971     //
1972     //                masterI = min(masterI, sharedPointLabels[sharedI]);
1973     //            }
1974     //
1975     //            forAll(mergeSet, i)
1976     //            {
1977     //                label sharedI = mergeSet[i];
1978     //
1979     //                pointToMaster.insert(sharedPointLabels[sharedI], masterI);
1980     //            }
1981     //        }
1982     //    }
1983     //
1984     //    //if (debug)
1985     //    //{
1986     //    //    Pout<< "polyMeshAdder : merging:"
1987     //    //        << pointToMaster.size() << " into " << nMergeSets
1988     //    //        << " sets." << endl;
1989     //    //}
1990     //}
1991
1992     return pointToMaster;
1993 }
1994
1995
1996 void Foam::polyMeshAdder::mergePoints
1997 (
1998     const polyMesh& mesh,
1999     const Map<label>& pointToMaster,
2000     polyTopoChange& meshMod
2001 )
2002 {
2003     // Remove all non-master points.
2004     forAll(mesh.points(), pointI)
2005     {
2006         Map<label>::const_iterator iter = pointToMaster.find(pointI);
2007
2008         if (iter != pointToMaster.end())
2009         {
2010             if (iter() != pointI)
2011             {
2012                 meshMod.removePoint(pointI, iter());
2013             }
2014         }
2015     }
2016
2017     // Modify faces for points. Note: could use pointFaces here but want to
2018     // avoid addressing calculation.
2019     const faceList& faces = mesh.faces();
2020
2021     forAll(faces, faceI)
2022     {
2023         const face& f = faces[faceI];
2024
2025         bool hasMerged = false;
2026
2027         forAll(f, fp)
2028         {
2029             label pointI = f[fp];
2030
2031             Map<label>::const_iterator iter = pointToMaster.find(pointI);
2032
2033             if (iter != pointToMaster.end())
2034             {
2035                 if (iter() != pointI)
2036                 {
2037                     hasMerged = true;
2038                     break;
2039                 }
2040             }
2041         }
2042
2043         if (hasMerged)
2044         {
2045             face newF(f);
2046
2047             forAll(f, fp)
2048             {
2049                 label pointI = f[fp];
2050
2051                 Map<label>::const_iterator iter = pointToMaster.find(pointI);
2052
2053                 if (iter != pointToMaster.end())
2054                 {
2055                     newF[fp] = iter();
2056                 }
2057             }
2058
2059             label patchID = mesh.boundaryMesh().whichPatch(faceI);
2060             label nei = (patchID == -1 ? mesh.faceNeighbour()[faceI] : -1);
2061             label zoneID = mesh.faceZones().whichZone(faceI);
2062             bool zoneFlip = false;
2063
2064             if (zoneID >= 0)
2065             {
2066                 const faceZone& fZone = mesh.faceZones()[zoneID];
2067                 zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
2068             }
2069
2070             meshMod.setAction
2071             (
2072                 polyModifyFace
2073                 (
2074                     newF,                       // modified face
2075                     faceI,                      // label of face
2076                     mesh.faceOwner()[faceI],    // owner
2077                     nei,                        // neighbour
2078                     false,                      // face flip
2079                     patchID,                    // patch for face
2080                     false,                      // remove from zone
2081                     zoneID,                     // zone for face
2082                     zoneFlip                    // face flip in zone
2083                 )
2084             );
2085         }
2086     }
2087 }
2088
2089
2090 // ************************************************************************* //