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-2008 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.setSize(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         workFace.shrink();
 516         allF.transfer(workFace);
 517     }
 518 }
 519
 520
 521 // Adds primitives (cells, faces, points)
 522 // Cells:
 523 //  - all of mesh0
 524 //  - all of mesh1
 525 // Faces:
 526 //  - all uncoupled of mesh0
 527 //  - all coupled faces
 528 //  - all uncoupled of mesh1
 529 // Points:
 530 //  - all coupled
 531 //  - all uncoupled of mesh0
 532 //  - all uncoupled of mesh1
 533 void Foam::polyMeshAdder::mergePrimitives
 534 (
 535     const polyMesh& mesh0,
 536     const polyMesh& mesh1,
 537     const faceCoupleInfo& coupleInfo,
 538
 539     const label nAllPatches,                // number of patches in the new mesh
 540     const labelList& fromAllTo1Patches,
 541     const labelList& from1ToAllPatches,
 542
 543     pointField& allPoints,
 544     labelList& from0ToAllPoints,
 545     labelList& from1ToAllPoints,
 546
 547     faceList& allFaces,
 548     labelList& allOwner,
 549     labelList& allNeighbour,
 550     label& nInternalFaces,
 551     labelList& nFacesPerPatch,
 552     label& nCells,
 553
 554     labelList& from0ToAllFaces,
 555     labelList& from1ToAllFaces,
 556     labelList& from1ToAllCells
 557 )
 558 {
 559     const polyBoundaryMesh& patches0 = mesh0.boundaryMesh();
 560     const polyBoundaryMesh& patches1 = mesh1.boundaryMesh();
 561
 562     const primitiveFacePatch& cutFaces = coupleInfo.cutFaces();
 563     const indirectPrimitivePatch& masterPatch = coupleInfo.masterPatch();
 564     const indirectPrimitivePatch& slavePatch = coupleInfo.slavePatch();
 565
 566
 567     // Points
 568     // ~~~~~~
 569
 570     // Storage for new points
 571     allPoints.setSize(mesh0.nPoints() + mesh1.nPoints());
 572     label allPointI = 0;
 573
 574     from0ToAllPoints.setSize(mesh0.nPoints());
 575     from0ToAllPoints = -1;
 576     from1ToAllPoints.setSize(mesh1.nPoints());
 577     from1ToAllPoints = -1;
 578
 579     // Copy coupled points (on cut)
 580     {
 581         const pointField& cutPoints = coupleInfo.cutPoints();
 582
 583         //const labelListList& cutToMasterPoints =
 584         //   coupleInfo.cutToMasterPoints();
 585         labelListList cutToMasterPoints
 586         (
 587             invertOneToMany
 588             (
 589                 cutPoints.size(),
 590                 coupleInfo.masterToCutPoints()
 591             )
 592         );
 593
 594         //const labelListList& cutToSlavePoints =
 595         //    coupleInfo.cutToSlavePoints();
 596         labelListList cutToSlavePoints
 597         (
 598             invertOneToMany
 599             (
 600                 cutPoints.size(),
 601                 coupleInfo.slaveToCutPoints()
 602             )
 603         );
 604
 605         forAll(cutPoints, i)
 606         {
 607             allPoints[allPointI] = cutPoints[i];
 608
 609             // Mark all master and slave points referring to this point.
 610
 611             const labelList& masterPoints = cutToMasterPoints[i];
 612
 613             forAll(masterPoints, j)
 614             {
 615                 label mesh0PointI = masterPatch.meshPoints()[masterPoints[j]];
 616                 from0ToAllPoints[mesh0PointI] = allPointI;
 617             }
 618
 619             const labelList& slavePoints = cutToSlavePoints[i];
 620
 621             forAll(slavePoints, j)
 622             {
 623                 label mesh1PointI = slavePatch.meshPoints()[slavePoints[j]];
 624                 from1ToAllPoints[mesh1PointI] = allPointI;
 625             }
 626             allPointI++;
 627         }
 628     }
 629
 630     // Add uncoupled mesh0 points
 631     forAll(mesh0.points(), pointI)
 632     {
 633         if (from0ToAllPoints[pointI] == -1)
 634         {
 635             allPoints[allPointI] = mesh0.points()[pointI];
 636             from0ToAllPoints[pointI] = allPointI;
 637             allPointI++;
 638         }
 639     }
 640
 641     // Add uncoupled mesh1 points
 642     forAll(mesh1.points(), pointI)
 643     {
 644         if (from1ToAllPoints[pointI] == -1)
 645         {
 646             allPoints[allPointI] = mesh1.points()[pointI];
 647             from1ToAllPoints[pointI] = allPointI;
 648             allPointI++;
 649         }
 650     }
 651
 652     allPoints.setSize(allPointI);
 653
 654
 655     // Faces
 656     // ~~~~~
 657
 658     // Sizes per patch
 659     nFacesPerPatch.setSize(nAllPatches);
 660     nFacesPerPatch = 0;
 661
 662     // Storage for faces and owner/neighbour
 663     allFaces.setSize(mesh0.nFaces() + mesh1.nFaces());
 664     allOwner.setSize(allFaces.size());
 665     allOwner = -1;
 666     allNeighbour.setSize(allFaces.size());
 667     allNeighbour = -1;
 668     label allFaceI = 0;
 669
 670     from0ToAllFaces.setSize(mesh0.nFaces());
 671     from0ToAllFaces = -1;
 672     from1ToAllFaces.setSize(mesh1.nFaces());
 673     from1ToAllFaces = -1;
 674
 675     // Copy mesh0 internal faces (always uncoupled)
 676     for (label faceI = 0; faceI < mesh0.nInternalFaces(); faceI++)
 677     {
 678         allFaces[allFaceI] = renumber(from0ToAllPoints, mesh0.faces()[faceI]);
 679         allOwner[allFaceI] = mesh0.faceOwner()[faceI];
 680         allNeighbour[allFaceI] = mesh0.faceNeighbour()[faceI];
 681         from0ToAllFaces[faceI] = allFaceI++;
 682     }
 683
 684     // Copy coupled faces. Every coupled face has an equivalent master and
 685     // slave. Also uncount as boundary faces all the newly coupled faces.
 686     const labelList& cutToMasterFaces = coupleInfo.cutToMasterFaces();
 687     const labelList& cutToSlaveFaces = coupleInfo.cutToSlaveFaces();
 688
 689     forAll(cutFaces, i)
 690     {
 691         label masterFaceI = cutToMasterFaces[i];
 692
 693         label mesh0FaceI = masterPatch.addressing()[masterFaceI];
 694
 695         if (from0ToAllFaces[mesh0FaceI] == -1)
 696         {
 697             // First occurrence of face
 698             from0ToAllFaces[mesh0FaceI] = allFaceI;
 699
 700             // External face becomes internal so uncount
 701             label patch0 = patches0.whichPatch(mesh0FaceI);
 702             nFacesPerPatch[patch0]--;
 703         }
 704
 705         label slaveFaceI = cutToSlaveFaces[i];
 706
 707         label mesh1FaceI = slavePatch.addressing()[slaveFaceI];
 708
 709         if (from1ToAllFaces[mesh1FaceI] == -1)
 710         {
 711             from1ToAllFaces[mesh1FaceI] = allFaceI;
 712
 713             label patch1 = patches1.whichPatch(mesh1FaceI);
 714             nFacesPerPatch[from1ToAllPatches[patch1]]--;
 715         }
 716
 717         // Copy cut face (since cutPoints are copied first no renumbering
 718         // nessecary)
 719         allFaces[allFaceI] = cutFaces[i];
 720         allOwner[allFaceI] = mesh0.faceOwner()[mesh0FaceI];
 721         allNeighbour[allFaceI] = mesh1.faceOwner()[mesh1FaceI] + mesh0.nCells();
 722
 723         allFaceI++;
 724     }
 725
 726     // Copy mesh1 internal faces (always uncoupled)
 727     for (label faceI = 0; faceI < mesh1.nInternalFaces(); faceI++)
 728     {
 729         allFaces[allFaceI] = renumber(from1ToAllPoints, mesh1.faces()[faceI]);
 730         allOwner[allFaceI] = mesh1.faceOwner()[faceI] + mesh0.nCells();
 731         allNeighbour[allFaceI] = mesh1.faceNeighbour()[faceI] + mesh0.nCells();
 732         from1ToAllFaces[faceI] = allFaceI++;
 733     }
 734
 735     nInternalFaces = allFaceI;
 736
 737     // Copy (unmarked/uncoupled) external faces in new order.
 738     for (label allPatchI = 0; allPatchI < nAllPatches; allPatchI++)
 739     {
 740         if (allPatchI < patches0.size())
 741         {
 742             // Patch is present in mesh0
 743             const polyPatch& pp = patches0[allPatchI];
 744
 745             nFacesPerPatch[allPatchI] += pp.size();
 746
 747             label faceI = pp.start();
 748
 749             forAll(pp, i)
 750             {
 751                 if (from0ToAllFaces[faceI] == -1)
 752                 {
 753                     // Is uncoupled face since has not yet been dealt with
 754                     allFaces[allFaceI] = renumber
 755                     (
 756                         from0ToAllPoints,
 757                         mesh0.faces()[faceI]
 758                     );
 759                     allOwner[allFaceI] = mesh0.faceOwner()[faceI];
 760                     allNeighbour[allFaceI] = -1;
 761
 762                     from0ToAllFaces[faceI] = allFaceI++;
 763                 }
 764                 faceI++;
 765             }
 766         }
 767         if (fromAllTo1Patches[allPatchI] != -1)
 768         {
 769             // Patch is present in mesh1
 770             const polyPatch& pp = patches1[fromAllTo1Patches[allPatchI]];
 771
 772             nFacesPerPatch[allPatchI] += pp.size();
 773
 774             label faceI = pp.start();
 775
 776             forAll(pp, i)
 777             {
 778                 if (from1ToAllFaces[faceI] == -1)
 779                 {
 780                     // Is uncoupled face
 781                     allFaces[allFaceI] = renumber
 782                     (
 783                         from1ToAllPoints,
 784                         mesh1.faces()[faceI]
 785                     );
 786                     allOwner[allFaceI] =
 787                         mesh1.faceOwner()[faceI]
 788                       + mesh0.nCells();
 789                     allNeighbour[allFaceI] = -1;
 790
 791                     from1ToAllFaces[faceI] = allFaceI++;
 792                 }
 793                 faceI++;
 794             }
 795         }
 796     }
 797     allFaces.setSize(allFaceI);
 798     allOwner.setSize(allFaceI);
 799     allNeighbour.setSize(allFaceI);
 800
 801
 802     // So now we have all ok for one-to-one mapping.
 803     // For split slace faces:
 804     // - mesh consistent with slave side
 805     // - mesh not consistent with owner side. It is not zipped up, the
 806     //   original faces need edges split.
 807
 808     // Use brute force to prevent having to calculate addressing:
 809     // - get map from master edge to split edges.
 810     // - check all faces to find any edge that is split.
 811     {
 812         // From two cut-points to labels of cut-points inbetween.
 813         // (in order: from e[0] to e[1]
 814         const edgeLookup& cutEdgeToPoints = coupleInfo.cutEdgeToPoints();
 815
 816         // Get map of master face (in mesh labels) that are in cut. These faces
 817         // do not need to be renumbered.
 818         labelHashSet masterCutFaces(cutToMasterFaces.size());
 819         forAll(cutToMasterFaces, i)
 820         {
 821             label meshFaceI = masterPatch.addressing()[cutToMasterFaces[i]];
 822
 823             masterCutFaces.insert(meshFaceI);
 824         }
 825
 826         DynamicList<label> workFace(100);
 827
 828         forAll(from0ToAllFaces, face0)
 829         {
 830             if (!masterCutFaces.found(face0))
 831             {
 832                 label allFaceI = from0ToAllFaces[face0];
 833
 834                 insertVertices
 835                 (
 836                     cutEdgeToPoints,
 837                     masterPatch.meshPointMap(),
 838                     coupleInfo.masterToCutPoints(),
 839                     mesh0.faces()[face0],
 840
 841                     workFace,
 842                     allFaces[allFaceI]
 843                 );
 844             }
 845         }
 846
 847         // Same for slave face
 848
 849         labelHashSet slaveCutFaces(cutToSlaveFaces.size());
 850         forAll(cutToSlaveFaces, i)
 851         {
 852             label meshFaceI = slavePatch.addressing()[cutToSlaveFaces[i]];
 853
 854             slaveCutFaces.insert(meshFaceI);
 855         }
 856
 857         forAll(from1ToAllFaces, face1)
 858         {
 859             if (!slaveCutFaces.found(face1))
 860             {
 861                 label allFaceI = from1ToAllFaces[face1];
 862
 863                 insertVertices
 864                 (
 865                     cutEdgeToPoints,
 866                     slavePatch.meshPointMap(),
 867                     coupleInfo.slaveToCutPoints(),
 868                     mesh1.faces()[face1],
 869
 870                     workFace,
 871                     allFaces[allFaceI]
 872                 );
 873             }
 874         }
 875     }
 876
 877     // Now we have a full facelist and owner/neighbour addressing.
 878
 879
 880     // Cells
 881     // ~~~~~
 882
 883     from1ToAllCells.setSize(mesh1.nCells());
 884     from1ToAllCells = -1;
 885
 886     forAll(mesh1.cells(), i)
 887     {
 888         from1ToAllCells[i] = i + mesh0.nCells();
 889     }
 890
 891     // Make cells (= cell-face addressing)
 892     nCells = mesh0.nCells() + mesh1.nCells();
 893     cellList allCells(nCells);
 894     primitiveMesh::calcCells(allCells, allOwner, allNeighbour, nCells);
 895
 896     // Reorder faces for upper-triangular order.
 897     labelList oldToNew
 898     (
 899         getFaceOrder
 900         (
 901             allCells,
 902             nInternalFaces,
 903             allOwner,
 904             allNeighbour
 905         )
 906     );
 907
 908     inplaceReorder(oldToNew, allFaces);
 909     inplaceReorder(oldToNew, allOwner);
 910     inplaceReorder(oldToNew, allNeighbour);
 911     inplaceRenumber(oldToNew, from0ToAllFaces);
 912     inplaceRenumber(oldToNew, from1ToAllFaces);
 913 }
 914
 915
 916 void Foam::polyMeshAdder::mergePointZones
 917 (
 918     const pointZoneMesh& pz0,
 919     const pointZoneMesh& pz1,
 920     const labelList& from0ToAllPoints,
 921     const labelList& from1ToAllPoints,
 922
 923     DynamicList<word>& zoneNames,
 924     labelList& from1ToAll,
 925     List<DynamicList<label> >& pzPoints
 926 )
 927 {
 928     zoneNames.setSize(pz0.size() + pz1.size());
 929
 930     // Names
 931     append(pz0.names(), zoneNames);
 932
 933     from1ToAll.setSize(pz1.size());
 934
 935     forAll(pz1, zoneI)
 936     {
 937         from1ToAll[zoneI] = zoneIndex(pz1[zoneI].name(), zoneNames);
 938     }
 939     zoneNames.shrink();
 940
 941     // Point labels per merged zone
 942     pzPoints.setSize(zoneNames.size());
 943
 944     forAll(pz0, zoneI)
 945     {
 946         append(from0ToAllPoints, pz0[zoneI], pzPoints[zoneI]);
 947     }
 948
 949     // Get sorted zone contents for duplicate element recognition
 950     PtrList<SortableList<label> > pzPointsSorted(pzPoints.size());
 951     forAll(pzPoints, zoneI)
 952     {
 953         pzPointsSorted.set
 954         (
 955             zoneI,
 956             new SortableList<label>(pzPoints[zoneI].shrink())
 957         );
 958     }
 959
 960     // Now we have full addressing for points so do the pointZones of mesh1.
 961     forAll(pz1, zoneI)
 962     {
 963         // Relabel all points of zone and add to correct pzPoints.
 964         label allZoneI = from1ToAll[zoneI];
 965
 966         append
 967         (
 968             from1ToAllPoints,
 969             pz1[zoneI],
 970             pzPointsSorted[allZoneI],
 971             pzPoints[allZoneI]
 972         );
 973     }
 974
 975     forAll(pzPoints, i)
 976     {
 977         pzPoints[i].shrink();
 978     }
 979 }
 980
 981
 982 void Foam::polyMeshAdder::mergeFaceZones
 983 (
 984     const faceZoneMesh& fz0,
 985     const faceZoneMesh& fz1,
 986     const labelList& from0ToAllFaces,
 987     const labelList& from1ToAllFaces,
 988
 989     DynamicList<word>& zoneNames,
 990     labelList& from1ToAll,
 991     List<DynamicList<label> >& fzFaces,
 992     List<DynamicList<bool> >& fzFlips
 993 )
 994 {
 995     zoneNames.setSize(fz0.size() + fz1.size());
 996
 997     append(fz0.names(), zoneNames);
 998
 999     from1ToAll.setSize(fz1.size());
1000
1001     forAll (fz1, zoneI)
1002     {
1003         from1ToAll[zoneI] = zoneIndex(fz1[zoneI].name(), zoneNames);
1004     }
1005     zoneNames.shrink();
1006
1007
1008     // Create temporary lists for faceZones.
1009     fzFaces.setSize(zoneNames.size());
1010     fzFlips.setSize(zoneNames.size());
1011     forAll(fz0, zoneI)
1012     {
1013         DynamicList<label>& newZone = fzFaces[zoneI];
1014         DynamicList<bool>& newFlip = fzFlips[zoneI];
1015
1016         newZone.setSize(fz0[zoneI].size());
1017         newFlip.setSize(newZone.size());
1018
1019         const labelList& addressing = fz0[zoneI];
1020         const boolList& flipMap = fz0[zoneI].flipMap();
1021
1022         forAll(addressing, i)
1023         {
1024             label faceI = addressing[i];
1025
1026             if (from0ToAllFaces[faceI] != -1)
1027             {
1028                 newZone.append(from0ToAllFaces[faceI]);
1029                 newFlip.append(flipMap[i]);
1030             }
1031         }
1032     }
1033
1034     // Get sorted zone contents for duplicate element recognition
1035     PtrList<SortableList<label> > fzFacesSorted(fzFaces.size());
1036     forAll(fzFaces, zoneI)
1037     {
1038         fzFacesSorted.set
1039         (
1040             zoneI,
1041             new SortableList<label>(fzFaces[zoneI].shrink())
1042         );
1043     }
1044
1045     // Now we have full addressing for faces so do the faceZones of mesh1.
1046     forAll(fz1, zoneI)
1047     {
1048         label allZoneI = from1ToAll[zoneI];
1049         DynamicList<label>& newZone = fzFaces[allZoneI];
1050         const SortableList<label>& newZoneSorted = fzFacesSorted[allZoneI];
1051         DynamicList<bool>& newFlip = fzFlips[allZoneI];
1052
1053         newZone.setSize(newZone.size() + fz1[zoneI].size());
1054         newFlip.setSize(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.setSize(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     // Cell mapping is trivial. Also no duplicate elements possible.
1116     forAll(cz1, zoneI)
1117     {
1118         label allZoneI = from1ToAll[zoneI];
1119
1120         append(from1ToAllCells, cz1[zoneI], czCells[allZoneI]);
1121     }
1122
1123     forAll(czCells, i)
1124     {
1125         czCells[i].shrink();
1126     }
1127 }
1128
1129
1130 void Foam::polyMeshAdder::mergeZones
1131 (
1132     const polyMesh& mesh0,
1133     const polyMesh& mesh1,
1134     const labelList& from0ToAllPoints,
1135     const labelList& from0ToAllFaces,
1136
1137     const labelList& from1ToAllPoints,
1138     const labelList& from1ToAllFaces,
1139     const labelList& from1ToAllCells,
1140
1141     DynamicList<word>& pointZoneNames,
1142     List<DynamicList<label> >& pzPoints,
1143
1144     DynamicList<word>& faceZoneNames,
1145     List<DynamicList<label> >& fzFaces,
1146     List<DynamicList<bool> >& fzFlips,
1147
1148     DynamicList<word>& cellZoneNames,
1149     List<DynamicList<label> >& czCells
1150 )
1151 {
1152     labelList from1ToAllPZones;
1153     mergePointZones
1154     (
1155         mesh0.pointZones(),
1156         mesh1.pointZones(),
1157         from0ToAllPoints,
1158         from1ToAllPoints,
1159
1160         pointZoneNames,
1161         from1ToAllPZones,
1162         pzPoints
1163     );
1164
1165     labelList from1ToAllFZones;
1166     mergeFaceZones
1167     (
1168         mesh0.faceZones(),
1169         mesh1.faceZones(),
1170         from0ToAllFaces,
1171         from1ToAllFaces,
1172
1173         faceZoneNames,
1174         from1ToAllFZones,
1175         fzFaces,
1176         fzFlips
1177     );
1178
1179     labelList from1ToAllCZones;
1180     mergeCellZones
1181     (
1182         mesh0.cellZones(),
1183         mesh1.cellZones(),
1184         from1ToAllCells,
1185
1186         cellZoneNames,
1187         from1ToAllCZones,
1188         czCells
1189     );
1190 }
1191
1192
1193 void Foam::polyMeshAdder::addZones
1194 (
1195     const DynamicList<word>& pointZoneNames,
1196     const List<DynamicList<label> >& pzPoints,
1197
1198     const DynamicList<word>& faceZoneNames,
1199     const List<DynamicList<label> >& fzFaces,
1200     const List<DynamicList<bool> >& fzFlips,
1201
1202     const DynamicList<word>& cellZoneNames,
1203     const List<DynamicList<label> >& czCells,
1204
1205     polyMesh& mesh
1206 )
1207 {
1208     List<pointZone*> pZones(pzPoints.size());
1209     forAll(pZones, i)
1210     {
1211         pZones[i] = new pointZone
1212         (
1213             pointZoneNames[i],
1214             pzPoints[i],
1215             i,
1216             mesh.pointZones()
1217         );
1218     }
1219
1220     List<faceZone*> fZones(fzFaces.size());
1221     forAll(fZones, i)
1222     {
1223         fZones[i] = new faceZone
1224         (
1225             faceZoneNames[i],
1226             fzFaces[i],
1227             fzFlips[i],
1228             i,
1229             mesh.faceZones()
1230         );
1231     }
1232
1233     List<cellZone*> cZones(czCells.size());
1234     forAll(cZones, i)
1235     {
1236         cZones[i] = new cellZone
1237         (
1238             cellZoneNames[i],
1239             czCells[i],
1240             i,
1241             mesh.cellZones()
1242         );
1243     }
1244
1245     mesh.addZones
1246     (
1247         pZones,
1248         fZones,
1249         cZones
1250     );
1251 }
1252
1253
1254
1255 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
1256
1257 // Returns new mesh and sets
1258 // - map from new cell/face/point/patch to either mesh0 or mesh1
1259 //
1260 // mesh0Faces/mesh1Faces: corresponding faces on both meshes.
1261 Foam::autoPtr<Foam::polyMesh> Foam::polyMeshAdder::add
1262 (
1263     const IOobject& io,
1264     const polyMesh& mesh0,
1265     const polyMesh& mesh1,
1266     const faceCoupleInfo& coupleInfo,
1267     autoPtr<mapAddedPolyMesh>& mapPtr
1268 )
1269 {
1270     const polyBoundaryMesh& patches0 = mesh0.boundaryMesh();
1271     const polyBoundaryMesh& patches1 = mesh1.boundaryMesh();
1272
1273
1274     DynamicList<word> allPatchNames(patches0.size() + patches1.size());
1275     DynamicList<word> allPatchTypes(allPatchNames.size());
1276
1277     // Patch maps
1278     labelList from1ToAllPatches(patches1.size());
1279     labelList fromAllTo1Patches(allPatchNames.size(), -1);
1280
1281     mergePatchNames
1282     (
1283         patches0,
1284         patches1,
1285         allPatchNames,
1286         allPatchTypes,
1287         from1ToAllPatches,
1288         fromAllTo1Patches
1289     );
1290
1291
1292     // New points
1293     pointField allPoints;
1294
1295     // Map from mesh0/1 points to allPoints.
1296     labelList from0ToAllPoints(mesh0.nPoints(), -1);
1297     labelList from1ToAllPoints(mesh1.nPoints(), -1);
1298
1299     // New faces
1300     faceList allFaces;
1301     label nInternalFaces;
1302
1303     // New cells
1304     labelList allOwner;
1305     labelList allNeighbour;
1306     label nCells;
1307
1308     // Sizes per patch
1309     labelList nFaces(allPatchNames.size(), 0);
1310
1311     // Maps
1312     labelList from0ToAllFaces(mesh0.nFaces(), -1);
1313     labelList from1ToAllFaces(mesh1.nFaces(), -1);
1314
1315     // Map
1316     labelList from1ToAllCells(mesh1.nCells(), -1);
1317
1318     mergePrimitives
1319     (
1320         mesh0,
1321         mesh1,
1322         coupleInfo,
1323
1324         allPatchNames.size(),
1325         fromAllTo1Patches,
1326         from1ToAllPatches,
1327
1328         allPoints,
1329         from0ToAllPoints,
1330         from1ToAllPoints,
1331
1332         allFaces,
1333         allOwner,
1334         allNeighbour,
1335         nInternalFaces,
1336         nFaces,
1337         nCells,
1338
1339         from0ToAllFaces,
1340         from1ToAllFaces,
1341         from1ToAllCells
1342     );
1343
1344
1345     // Zones
1346     // ~~~~~
1347
1348     DynamicList<word> pointZoneNames;
1349     List<DynamicList<label> > pzPoints;
1350
1351     DynamicList<word> faceZoneNames;
1352     List<DynamicList<label> > fzFaces;
1353     List<DynamicList<bool> > fzFlips;
1354
1355     DynamicList<word> cellZoneNames;
1356     List<DynamicList<label> > czCells;
1357
1358     mergeZones
1359     (
1360         mesh0,
1361         mesh1,
1362
1363         from0ToAllPoints,
1364         from0ToAllFaces,
1365
1366         from1ToAllPoints,
1367         from1ToAllFaces,
1368         from1ToAllCells,
1369
1370         pointZoneNames,
1371         pzPoints,
1372
1373         faceZoneNames,
1374         fzFaces,
1375         fzFlips,
1376
1377         cellZoneNames,
1378         czCells
1379     );
1380
1381
1382     // Patches
1383     // ~~~~~~~
1384
1385     // Map from 0 to all patches (since gets compacted)
1386     labelList from0ToAllPatches(patches0.size(), -1);
1387
1388     List<polyPatch*> allPatches
1389     (
1390         combinePatches
1391         (
1392             mesh0,
1393             mesh1,
1394             patches0,           // Should be boundaryMesh() on new mesh.
1395             allPatchNames.size(),
1396             fromAllTo1Patches,
1397             mesh0.nInternalFaces()
1398           + mesh1.nInternalFaces()
1399           + coupleInfo.cutFaces().size(),
1400             nFaces,
1401
1402             from0ToAllPatches,
1403             from1ToAllPatches
1404         )
1405     );
1406
1407
1408     // Map information
1409     // ~~~~~~~~~~~~~~~
1410
1411     mapPtr.reset
1412     (
1413         new mapAddedPolyMesh
1414         (
1415             mesh0.nPoints(),
1416             mesh0.nFaces(),
1417             mesh0.nCells(),
1418
1419             mesh1.nPoints(),
1420             mesh1.nFaces(),
1421             mesh1.nCells(),
1422
1423             from0ToAllPoints,
1424             from0ToAllFaces,
1425             identity(mesh0.nCells()),
1426
1427             from1ToAllPoints,
1428             from1ToAllFaces,
1429             from1ToAllCells,
1430
1431             from0ToAllPatches,
1432             from1ToAllPatches,
1433             getPatchSizes(patches0),
1434             getPatchStarts(patches0)
1435         )
1436     );
1437
1438
1439
1440     // Now we have extracted all information from all meshes.
1441     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1442
1443     // Construct mesh
1444     autoPtr<polyMesh> tmesh
1445     (
1446         new polyMesh
1447         (
1448             io,
1449             allPoints,
1450             allFaces,
1451             allOwner,
1452             allNeighbour
1453         )
1454     );
1455     polyMesh& mesh = tmesh();
1456
1457     // Add zones to new mesh.
1458     addZones
1459     (
1460         pointZoneNames,
1461         pzPoints,
1462
1463         faceZoneNames,
1464         fzFaces,
1465         fzFlips,
1466
1467         cellZoneNames,
1468         czCells,
1469         mesh
1470     );
1471
1472     // Add patches to new mesh
1473     mesh.addPatches(allPatches);
1474
1475     return tmesh;
1476 }
1477
1478
1479 // Inplace add mesh1 to mesh0
1480 Foam::autoPtr<Foam::mapAddedPolyMesh> Foam::polyMeshAdder::add
1481 (
1482     polyMesh& mesh0,
1483     const polyMesh& mesh1,
1484     const faceCoupleInfo& coupleInfo,
1485     const bool validBoundary
1486 )
1487 {
1488     const polyBoundaryMesh& patches0 = mesh0.boundaryMesh();
1489     const polyBoundaryMesh& patches1 = mesh1.boundaryMesh();
1490
1491     DynamicList<word> allPatchNames(patches0.size() + patches1.size());
1492     DynamicList<word> allPatchTypes(allPatchNames.size());
1493
1494     // Patch maps
1495     labelList from1ToAllPatches(patches1.size());
1496     labelList fromAllTo1Patches(allPatchNames.size(), -1);
1497
1498     mergePatchNames
1499     (
1500         patches0,
1501         patches1,
1502         allPatchNames,
1503         allPatchTypes,
1504         from1ToAllPatches,
1505         fromAllTo1Patches
1506     );
1507
1508
1509     // New points
1510     pointField allPoints;
1511
1512     // Map from mesh0/1 points to allPoints.
1513     labelList from0ToAllPoints(mesh0.nPoints(), -1);
1514     labelList from1ToAllPoints(mesh1.nPoints(), -1);
1515
1516     // New faces
1517     faceList allFaces;
1518     labelList allOwner;
1519     labelList allNeighbour;
1520     label nInternalFaces;
1521     // Sizes per patch
1522     labelList nFaces(allPatchNames.size(), 0);
1523     label nCells;
1524
1525     // Maps
1526     labelList from0ToAllFaces(mesh0.nFaces(), -1);
1527     labelList from1ToAllFaces(mesh1.nFaces(), -1);
1528     // Map
1529     labelList from1ToAllCells(mesh1.nCells(), -1);
1530
1531     mergePrimitives
1532     (
1533         mesh0,
1534         mesh1,
1535         coupleInfo,
1536
1537         allPatchNames.size(),
1538         fromAllTo1Patches,
1539         from1ToAllPatches,
1540
1541         allPoints,
1542         from0ToAllPoints,
1543         from1ToAllPoints,
1544
1545         allFaces,
1546         allOwner,
1547         allNeighbour,
1548         nInternalFaces,
1549         nFaces,
1550         nCells,
1551
1552         from0ToAllFaces,
1553         from1ToAllFaces,
1554         from1ToAllCells
1555     );
1556
1557
1558     // Zones
1559     // ~~~~~
1560
1561     DynamicList<word> pointZoneNames;
1562     List<DynamicList<label> > pzPoints;
1563
1564     DynamicList<word> faceZoneNames;
1565     List<DynamicList<label> > fzFaces;
1566     List<DynamicList<bool> > fzFlips;
1567
1568     DynamicList<word> cellZoneNames;
1569     List<DynamicList<label> > czCells;
1570
1571     mergeZones
1572     (
1573         mesh0,
1574         mesh1,
1575
1576         from0ToAllPoints,
1577         from0ToAllFaces,
1578
1579         from1ToAllPoints,
1580         from1ToAllFaces,
1581         from1ToAllCells,
1582
1583         pointZoneNames,
1584         pzPoints,
1585
1586         faceZoneNames,
1587         fzFaces,
1588         fzFlips,
1589
1590         cellZoneNames,
1591         czCells
1592     );
1593
1594
1595     // Patches
1596     // ~~~~~~~
1597
1598
1599     // Store mesh0 patch info before modifying patches0.
1600     labelList mesh0PatchSizes(getPatchSizes(patches0));
1601     labelList mesh0PatchStarts(getPatchStarts(patches0));
1602
1603     // Map from 0 to all patches (since gets compacted)
1604     labelList from0ToAllPatches(patches0.size(), -1);
1605
1606     // Inplace extend mesh0 patches (note that patches0.size() now also
1607     // has changed)
1608     polyBoundaryMesh& allPatches =
1609         const_cast<polyBoundaryMesh&>(mesh0.boundaryMesh());
1610     allPatches.setSize(allPatchNames.size());
1611
1612     label startFaceI = nInternalFaces;
1613
1614     // Copy patches0 with new sizes. First patches always come from
1615     // mesh0 and will always be present.
1616     label allPatchI = 0;
1617
1618     forAll(from0ToAllPatches, patch0)
1619     {
1620         // Originates from mesh0. Clone with new size & filter out empty
1621         // patch.
1622
1623         if (nFaces[patch0] == 0 && isA<processorPolyPatch>(allPatches[patch0]))
1624         {
1625             //Pout<< "Removing zero sized mesh0 patch " << allPatchNames[patch0]
1626             //    << endl;
1627             from0ToAllPatches[patch0] = -1;
1628             // Check if patch was also in mesh1 and update its addressing if so.
1629             if (fromAllTo1Patches[patch0] != -1)
1630             {
1631                 from1ToAllPatches[fromAllTo1Patches[patch0]] = -1;
1632             }
1633         }
1634         else
1635         {
1636             // Clone.
1637             allPatches.set
1638             (
1639                 allPatchI,
1640                 allPatches[patch0].clone
1641                 (
1642                     allPatches,
1643                     allPatchI,
1644                     nFaces[patch0],
1645                     startFaceI
1646                 )
1647             );
1648
1649             // Record new index in allPatches
1650             from0ToAllPatches[patch0] = allPatchI;
1651
1652             // Check if patch was also in mesh1 and update its addressing if so.
1653             if (fromAllTo1Patches[patch0] != -1)
1654             {
1655                 from1ToAllPatches[fromAllTo1Patches[patch0]] = allPatchI;
1656             }
1657
1658             startFaceI += nFaces[patch0];
1659
1660             allPatchI++;
1661         }
1662     }
1663
1664     // Copy unique patches of mesh1.
1665     forAll(from1ToAllPatches, patch1)
1666     {
1667         label uncompactAllPatchI = from1ToAllPatches[patch1];
1668
1669         if (uncompactAllPatchI >= from0ToAllPatches.size())
1670         {
1671             // Patch has not been merged with any mesh0 patch.
1672
1673             if
1674             (
1675                 nFaces[uncompactAllPatchI] == 0
1676              && isA<processorPolyPatch>(patches1[patch1])
1677             )
1678             {
1679                 //Pout<< "Removing zero sized mesh1 patch "
1680                 //    << allPatchNames[uncompactAllPatchI] << endl;
1681                 from1ToAllPatches[patch1] = -1;
1682             }
1683             else
1684             {
1685                 // Clone.
1686                 allPatches.set
1687                 (
1688                     allPatchI,
1689                     patches1[patch1].clone
1690                     (
1691                         allPatches,
1692                         allPatchI,
1693                         nFaces[uncompactAllPatchI],
1694                         startFaceI
1695                     )
1696                 );
1697
1698                 // Record new index in allPatches
1699                 from1ToAllPatches[patch1] = allPatchI;
1700
1701                 startFaceI += nFaces[uncompactAllPatchI];
1702
1703                 allPatchI++;
1704             }
1705         }
1706     }
1707
1708     allPatches.setSize(allPatchI);
1709
1710
1711     // Construct map information before changing mesh0 primitives
1712     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1713
1714     autoPtr<mapAddedPolyMesh> mapPtr
1715     (
1716         new mapAddedPolyMesh
1717         (
1718             mesh0.nPoints(),
1719             mesh0.nFaces(),
1720             mesh0.nCells(),
1721
1722             mesh1.nPoints(),
1723             mesh1.nFaces(),
1724             mesh1.nCells(),
1725
1726             from0ToAllPoints,
1727             from0ToAllFaces,
1728             identity(mesh0.nCells()),
1729
1730             from1ToAllPoints,
1731             from1ToAllFaces,
1732             from1ToAllCells,
1733
1734             from0ToAllPatches,
1735             from1ToAllPatches,
1736
1737             mesh0PatchSizes,
1738             mesh0PatchStarts
1739         )
1740     );
1741
1742
1743
1744     // Now we have extracted all information from all meshes.
1745     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1746
1747     labelList patchSizes(getPatchSizes(allPatches));
1748     labelList patchStarts(getPatchStarts(allPatches));
1749
1750     mesh0.resetMotion();    // delete any oldPoints.
1751     mesh0.resetPrimitives
1752     (
1753         allFaces.size(),
1754         allPoints,
1755         allFaces,
1756         allOwner,
1757         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     //        IndirectList<point>
1933     //        (
1934     //            mesh.points(),
1935     //            sharedPointLabels
1936     //        )()
1937     //    ),
1938     //    mergeDist,
1939     //    false,
1940     //    sharedToMerged,
1941     //    mergedPoints
1942     //);
1943     //
1944     //// Find out which sets of points get merged and create a map from
1945     //// mesh point to unique point.
1946     //
1947     //Map<label> pointToMaster(10*sharedToMerged.size());
1948     //
1949     //if (hasMerged)
1950     //{
1951     //    labelListList mergeSets
1952     //    (
1953     //        invertOneToMany
1954     //        (
1955     //            sharedToMerged.size(),
1956     //            sharedToMerged
1957     //        )
1958     //    );
1959     //
1960     //    label nMergeSets = 0;
1961     //
1962     //    forAll(mergeSets, setI)
1963     //    {
1964     //        const labelList& mergeSet = mergeSets[setI];
1965     //
1966     //        if (mergeSet.size() > 1)
1967     //        {
1968     //            // Take as master the shared point with the lowest mesh
1969     //            // point label. (rather arbitrarily - could use max or
1970     //            // any other one of the points)
1971     //
1972     //            nMergeSets++;
1973     //
1974     //            label masterI = labelMax;
1975     //
1976     //            forAll(mergeSet, i)
1977     //            {
1978     //                label sharedI = mergeSet[i];
1979     //
1980     //                masterI = min(masterI, sharedPointLabels[sharedI]);
1981     //            }
1982     //
1983     //            forAll(mergeSet, i)
1984     //            {
1985     //                label sharedI = mergeSet[i];
1986     //
1987     //                pointToMaster.insert(sharedPointLabels[sharedI], masterI);
1988     //            }
1989     //        }
1990     //    }
1991     //
1992     //    //if (debug)
1993     //    //{
1994     //    //    Pout<< "polyMeshAdder : merging:"
1995     //    //        << pointToMaster.size() << " into " << nMergeSets
1996     //    //        << " sets." << endl;
1997     //    //}
1998     //}
1999
2000     return pointToMaster;
2001 }
2002
2003
2004 void Foam::polyMeshAdder::mergePoints
2005 (
2006     const polyMesh& mesh,
2007     const Map<label>& pointToMaster,
2008     polyTopoChange& meshMod
2009 )
2010 {
2011     // Remove all non-master points.
2012     forAll(mesh.points(), pointI)
2013     {
2014         Map<label>::const_iterator iter = pointToMaster.find(pointI);
2015
2016         if (iter != pointToMaster.end())
2017         {
2018             if (iter() != pointI)
2019             {
2020                 //1.4.1: meshMod.removePoint(pointI, iter());
2021                 meshMod.setAction(polyRemovePoint(pointI));
2022             }
2023         }
2024     }
2025
2026     // Modify faces for points. Note: could use pointFaces here but want to
2027     // avoid addressing calculation.
2028     const faceList& faces = mesh.faces();
2029
2030     forAll(faces, faceI)
2031     {
2032         const face& f = faces[faceI];
2033
2034         bool hasMerged = false;
2035
2036         forAll(f, fp)
2037         {
2038             label pointI = f[fp];
2039
2040             Map<label>::const_iterator iter = pointToMaster.find(pointI);
2041
2042             if (iter != pointToMaster.end())
2043             {
2044                 if (iter() != pointI)
2045                 {
2046                     hasMerged = true;
2047                     break;
2048                 }
2049             }
2050         }
2051
2052         if (hasMerged)
2053         {
2054             face newF(f);
2055
2056             forAll(f, fp)
2057             {
2058                 label pointI = f[fp];
2059
2060                 Map<label>::const_iterator iter = pointToMaster.find(pointI);
2061
2062                 if (iter != pointToMaster.end())
2063                 {
2064                     newF[fp] = iter();
2065                 }
2066             }
2067
2068             label patchID = mesh.boundaryMesh().whichPatch(faceI);
2069             label nei = (patchID == -1 ? mesh.faceNeighbour()[faceI] : -1);
2070             label zoneID = mesh.faceZones().whichZone(faceI);
2071             bool zoneFlip = false;
2072
2073             if (zoneID >= 0)
2074             {
2075                 const faceZone& fZone = mesh.faceZones()[zoneID];
2076                 zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
2077             }
2078
2079             meshMod.setAction
2080             (
2081                 polyModifyFace
2082                 (
2083                     newF,                       // modified face
2084                     faceI,                      // label of face
2085                     mesh.faceOwner()[faceI],    // owner
2086                     nei,                        // neighbour
2087                     false,                      // face flip
2088                     patchID,                    // patch for face
2089                     false,                      // remove from zone
2090                     zoneID,                     // zone for face
2091                     zoneFlip                    // face flip in zone
2092                 )
2093             );
2094         }
2095     }
2096 }
2097
2098
2099 // ************************************************************************* //