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