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-2010 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
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*----------------------------------------------------------------------------*/
  25
  26 #include "polyMeshAdder.H"
  27 #include "mapAddedPolyMesh.H"
  28 #include "IOobject.H"
  29 #include "faceCoupleInfo.H"
  30 #include "processorPolyPatch.H"
  31 #include "SortableList.H"
  32 #include "Time.H"
  33 #include "globalMeshData.H"
  34 #include "mergePoints.H"
  35 #include "polyModifyFace.H"
  36 #include "polyRemovePoint.H"
  37 #include "polyTopoChange.H"
  38
  39 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  40
  41 //- Append all mapped elements of a list to a DynamicList
  42 void Foam::polyMeshAdder::append
  43 (
  44     const labelList& map,
  45     const labelList& lst,
  46     DynamicList<label>& dynLst
  47 )
  48 {
  49     dynLst.setCapacity(dynLst.size() + lst.size());
  50
  51     forAll(lst, i)
  52     {
  53         label newElem = map[lst[i]];
  54
  55         if (newElem != -1)
  56         {
  57             dynLst.append(newElem);
  58         }
  59     }
  60 }
  61
  62
  63 //- Append all mapped elements of a list to a DynamicList
  64 void Foam::polyMeshAdder::append
  65 (
  66     const labelList& map,
  67     const labelList& lst,
  68     const SortableList<label>& sortedLst,
  69     DynamicList<label>& dynLst
  70 )
  71 {
  72     dynLst.setSize(dynLst.size() + lst.size());
  73
  74     forAll(lst, i)
  75     {
  76         label newElem = map[lst[i]];
  77
  78         if (newElem != -1 && findSortedIndex(sortedLst, newElem) == -1)
  79         {
  80             dynLst.append(newElem);
  81         }
  82     }
  83 }
  84
  85
  86 // Get index of patch in new set of patchnames/types
  87 Foam::label Foam::polyMeshAdder::patchIndex
  88 (
  89     const polyPatch& p,
  90     DynamicList<word>& allPatchNames,
  91     DynamicList<word>& allPatchTypes
  92 )
  93 {
  94     // Find the patch name on the list.  If the patch is already there
  95     // and patch types match, return index
  96     const word& pType = p.type();
  97     const word& pName = p.name();
  98
  99     label patchI = findIndex(allPatchNames, pName);
 100
 101     if (patchI == -1)
 102     {
 103         // Patch not found. Append to the list
 104         allPatchNames.append(pName);
 105         allPatchTypes.append(pType);
 106
 107         return allPatchNames.size() - 1;
 108     }
 109     else if (allPatchTypes[patchI] == pType)
 110     {
 111         // Found name and types match
 112         return patchI;
 113     }
 114     else
 115     {
 116         // Found the name, but type is different
 117
 118         // Duplicate name is not allowed.  Create a composite name from the
 119         // patch name and case name
 120         const word& caseName = p.boundaryMesh().mesh().time().caseName();
 121
 122         allPatchNames.append(pName + "_" + caseName);
 123         allPatchTypes.append(pType);
 124
 125         Pout<< "label patchIndex(const polyPatch& p) : "
 126             << "Patch " << p.index() << " named "
 127             << pName << " in mesh " << caseName
 128             << " already exists, but patch types"
 129             << " do not match.\nCreating a composite name as "
 130             << allPatchNames[allPatchNames.size() - 1] << endl;
 131
 132         return allPatchNames.size() - 1;
 133     }
 134 }
 135
 136
 137 // Get index of zone in new set of zone names
 138 Foam::label Foam::polyMeshAdder::zoneIndex
 139 (
 140     const word& curName,
 141     DynamicList<word>& names
 142 )
 143 {
 144     label zoneI = findIndex(names, curName);
 145
 146     if (zoneI != -1)
 147     {
 148         return zoneI;
 149     }
 150     else
 151     {
 152         // Not found.  Add new name to the list
 153         names.append(curName);
 154
 155         return names.size() - 1;
 156     }
 157 }
 158
 159
 160 void Foam::polyMeshAdder::mergePatchNames
 161 (
 162     const polyBoundaryMesh& patches0,
 163     const polyBoundaryMesh& patches1,
 164
 165     DynamicList<word>& allPatchNames,
 166     DynamicList<word>& allPatchTypes,
 167
 168     labelList& from1ToAllPatches,
 169     labelList& fromAllTo1Patches
 170 )
 171 {
 172     // Insert the mesh0 patches and zones
 173     append(patches0.names(), allPatchNames);
 174     append(patches0.types(), allPatchTypes);
 175
 176
 177     // Patches
 178     // ~~~~~~~
 179     // Patches from 0 are taken over as is; those from 1 get either merged
 180     // (if they share name and type) or appended.
 181     // Empty patches are filtered out much much later on.
 182
 183     // Add mesh1 patches and build map both ways.
 184     from1ToAllPatches.setSize(patches1.size());
 185
 186     forAll(patches1, patchI)
 187     {
 188         from1ToAllPatches[patchI] = patchIndex
 189         (
 190             patches1[patchI],
 191             allPatchNames,
 192             allPatchTypes
 193         );
 194     }
 195     allPatchTypes.shrink();
 196     allPatchNames.shrink();
 197
 198     // Invert 1 to all patch map
 199     fromAllTo1Patches.setSize(allPatchNames.size());
 200     fromAllTo1Patches = -1;
 201
 202     forAll(from1ToAllPatches, i)
 203     {
 204         fromAllTo1Patches[from1ToAllPatches[i]] = i;
 205     }
 206 }
 207
 208
 209 Foam::labelList Foam::polyMeshAdder::getPatchStarts
 210 (
 211     const polyBoundaryMesh& patches
 212 )
 213 {
 214     labelList patchStarts(patches.size());
 215     forAll(patches, patchI)
 216     {
 217         patchStarts[patchI] = patches[patchI].start();
 218     }
 219     return patchStarts;
 220 }
 221
 222
 223 Foam::labelList Foam::polyMeshAdder::getPatchSizes
 224 (
 225     const polyBoundaryMesh& patches
 226 )
 227 {
 228     labelList patchSizes(patches.size());
 229     forAll(patches, patchI)
 230     {
 231         patchSizes[patchI] = patches[patchI].size();
 232     }
 233     return patchSizes;
 234 }
 235
 236
 237 Foam::List<Foam::polyPatch*> Foam::polyMeshAdder::combinePatches
 238 (
 239     const polyMesh& mesh0,
 240     const polyMesh& mesh1,
 241     const polyBoundaryMesh& allBoundaryMesh,
 242     const label nAllPatches,
 243     const labelList& fromAllTo1Patches,
 244
 245     const label nInternalFaces,
 246     const labelList& nFaces,
 247
 248     labelList& from0ToAllPatches,
 249     labelList& from1ToAllPatches
 250 )
 251 {
 252     const polyBoundaryMesh& patches0 = mesh0.boundaryMesh();
 253     const polyBoundaryMesh& patches1 = mesh1.boundaryMesh();
 254
 255
 256     // Compacted new patch list.
 257     DynamicList<polyPatch*> allPatches(nAllPatches);
 258
 259
 260     // Map from 0 to all patches (since gets compacted)
 261     from0ToAllPatches.setSize(patches0.size());
 262     from0ToAllPatches = -1;
 263
 264     label startFaceI = nInternalFaces;
 265
 266     // Copy patches0 with new sizes. First patches always come from
 267     // mesh0 and will always be present.
 268     for (label patchI = 0; patchI < patches0.size(); patchI++)
 269     {
 270         // Originates from mesh0. Clone with new size & filter out empty
 271         // patch.
 272         label filteredPatchI;
 273
 274         if (nFaces[patchI] == 0 && isA<processorPolyPatch>(patches0[patchI]))
 275         {
 276             //Pout<< "Removing zero sized mesh0 patch "
 277             //    << patches0[patchI].name() << endl;
 278             filteredPatchI = -1;
 279         }
 280         else
 281         {
 282             filteredPatchI = allPatches.size();
 283
 284             allPatches.append
 285             (
 286                 patches0[patchI].clone
 287                 (
 288                     allBoundaryMesh,
 289                     filteredPatchI,
 290                     nFaces[patchI],
 291                     startFaceI
 292                 ).ptr()
 293             );
 294             startFaceI += nFaces[patchI];
 295         }
 296
 297         // Record new index in allPatches
 298         from0ToAllPatches[patchI] = filteredPatchI;
 299
 300         // Check if patch was also in mesh1 and update its addressing if so.
 301         if (fromAllTo1Patches[patchI] != -1)
 302         {
 303             from1ToAllPatches[fromAllTo1Patches[patchI]] = filteredPatchI;
 304         }
 305     }
 306
 307     // Copy unique patches of mesh1.
 308     forAll(from1ToAllPatches, patchI)
 309     {
 310         label allPatchI = from1ToAllPatches[patchI];
 311
 312         if (allPatchI >= patches0.size())
 313         {
 314             // Patch has not been merged with any mesh0 patch.
 315
 316             label filteredPatchI;
 317
 318             if
 319             (
 320                 nFaces[allPatchI] == 0
 321              && isA<processorPolyPatch>(patches1[patchI])
 322             )
 323             {
 324                 //Pout<< "Removing zero sized mesh1 patch "
 325                 //    << patches1[patchI].name() << endl;
 326                 filteredPatchI = -1;
 327             }
 328             else
 329             {
 330                 filteredPatchI = allPatches.size();
 331
 332                 allPatches.append
 333                 (
 334                     patches1[patchI].clone
 335                     (
 336                         allBoundaryMesh,
 337                         filteredPatchI,
 338                         nFaces[allPatchI],
 339                         startFaceI
 340                     ).ptr()
 341                 );
 342                 startFaceI += nFaces[allPatchI];
 343             }
 344
 345             from1ToAllPatches[patchI] = filteredPatchI;
 346         }
 347     }
 348
 349     allPatches.shrink();
 350
 351     return allPatches;
 352 }
 353
 354
 355 Foam::labelList Foam::polyMeshAdder::getFaceOrder
 356 (
 357     const cellList& cells,
 358     const label nInternalFaces,
 359     const labelList& owner,
 360     const labelList& neighbour
 361 )
 362 {
 363     labelList oldToNew(owner.size(), -1);
 364
 365     // Leave boundary faces in order
 366     for (label faceI = nInternalFaces; faceI < owner.size(); faceI++)
 367     {
 368         oldToNew[faceI] = faceI;
 369     }
 370
 371     // First unassigned face
 372     label newFaceI = 0;
 373
 374     forAll(cells, cellI)
 375     {
 376         const labelList& cFaces = cells[cellI];
 377
 378         SortableList<label> nbr(cFaces.size());
 379
 380         forAll(cFaces, i)
 381         {
 382             label faceI = cFaces[i];
 383
 384             label nbrCellI = neighbour[faceI];
 385
 386             if (nbrCellI != -1)
 387             {
 388                 // Internal face. Get cell on other side.
 389                 if (nbrCellI == cellI)
 390                 {
 391                     nbrCellI = owner[faceI];
 392                 }
 393
 394                 if (cellI < nbrCellI)
 395                 {
 396                     // CellI is master
 397                     nbr[i] = nbrCellI;
 398                 }
 399                 else
 400                 {
 401                     // nbrCell is master. Let it handle this face.
 402                     nbr[i] = -1;
 403                 }
 404             }
 405             else
 406             {
 407                 // External face. Do later.
 408                 nbr[i] = -1;
 409             }
 410         }
 411
 412         nbr.sort();
 413
 414         forAll(nbr, i)
 415         {
 416             if (nbr[i] != -1)
 417             {
 418                 oldToNew[cFaces[nbr.indices()[i]]] = newFaceI++;
 419             }
 420         }
 421     }
 422
 423
 424     // Check done all faces.
 425     forAll(oldToNew, faceI)
 426     {
 427         if (oldToNew[faceI] == -1)
 428         {
 429             FatalErrorIn
 430             (
 431                 "polyMeshAdder::getFaceOrder"
 432                 "(const cellList&, const label, const labelList&"
 433                 ", const labelList&) const"
 434             )   << "Did not determine new position"
 435                 << " for face " << faceI
 436                 << abort(FatalError);
 437         }
 438     }
 439
 440     return oldToNew;
 441 }
 442
 443
 444 // Extends face f with split points. cutEdgeToPoints gives for every
 445 // edge the points introduced inbetween the endpoints.
 446 void Foam::polyMeshAdder::insertVertices
 447 (
 448     const edgeLookup& cutEdgeToPoints,
 449     const Map<label>& meshToMaster,
 450     const labelList& masterToCutPoints,
 451     const face& masterF,
 452
 453     DynamicList<label>& workFace,
 454     face& allF
 455 )
 456 {
 457     workFace.clear();
 458
 459     // Check any edge for being cut (check on the cut so takes account
 460     // for any point merging on the cut)
 461
 462     forAll(masterF, fp)
 463     {
 464         label v0 = masterF[fp];
 465         label v1 = masterF.nextLabel(fp);
 466
 467         // Copy existing face point
 468         workFace.append(allF[fp]);
 469
 470         // See if any edge between v0,v1
 471
 472         Map<label>::const_iterator v0Fnd = meshToMaster.find(v0);
 473         if (v0Fnd != meshToMaster.end())
 474         {
 475             Map<label>::const_iterator v1Fnd = meshToMaster.find(v1);
 476             if (v1Fnd != meshToMaster.end())
 477             {
 478                 // Get edge in cutPoint numbering
 479                 edge cutEdge
 480                 (
 481                     masterToCutPoints[v0Fnd()],
 482                     masterToCutPoints[v1Fnd()]
 483                 );
 484
 485                 edgeLookup::const_iterator iter = cutEdgeToPoints.find(cutEdge);
 486
 487                 if (iter != cutEdgeToPoints.end())
 488                 {
 489                     const edge& e = iter.key();
 490                     const labelList& addedPoints = iter();
 491
 492                     // cutPoints first in allPoints so no need for renumbering
 493                     if (e[0] == cutEdge[0])
 494                     {
 495                         forAll(addedPoints, i)
 496                         {
 497                             workFace.append(addedPoints[i]);
 498                         }
 499                     }
 500                     else
 501                     {
 502                         forAllReverse(addedPoints, i)
 503                         {
 504                             workFace.append(addedPoints[i]);
 505                         }
 506                     }
 507                 }
 508             }
 509         }
 510     }
 511
 512     if (workFace.size() != allF.size())
 513     {
 514         allF.transfer(workFace);
 515     }
 516 }
 517
 518
 519 // Adds primitives (cells, faces, points)
 520 // Cells:
 521 //  - all of mesh0
 522 //  - all of mesh1
 523 // Faces:
 524 //  - all uncoupled of mesh0
 525 //  - all coupled faces
 526 //  - all uncoupled of mesh1
 527 // Points:
 528 //  - all coupled
 529 //  - all uncoupled of mesh0
 530 //  - all uncoupled of mesh1
 531 void Foam::polyMeshAdder::mergePrimitives
 532 (
 533     const polyMesh& mesh0,
 534     const polyMesh& mesh1,
 535     const faceCoupleInfo& coupleInfo,
 536
 537     const label nAllPatches,                // number of patches in the new mesh
 538     const labelList& fromAllTo1Patches,
 539     const labelList& from1ToAllPatches,
 540
 541     pointField& allPoints,
 542     labelList& from0ToAllPoints,
 543     labelList& from1ToAllPoints,
 544
 545     faceList& allFaces,
 546     labelList& allOwner,
 547     labelList& allNeighbour,
 548     label& nInternalFaces,
 549     labelList& nFacesPerPatch,
 550     label& nCells,
 551
 552     labelList& from0ToAllFaces,
 553     labelList& from1ToAllFaces,
 554     labelList& from1ToAllCells
 555 )
 556 {
 557     const polyBoundaryMesh& patches0 = mesh0.boundaryMesh();
 558     const polyBoundaryMesh& patches1 = mesh1.boundaryMesh();
 559
 560     const primitiveFacePatch& cutFaces = coupleInfo.cutFaces();
 561     const indirectPrimitivePatch& masterPatch = coupleInfo.masterPatch();
 562     const indirectPrimitivePatch& slavePatch = coupleInfo.slavePatch();
 563
 564
 565     // Points
 566     // ~~~~~~
 567
 568     // Storage for new points
 569     allPoints.setSize(mesh0.nPoints() + mesh1.nPoints());
 570     label allPointI = 0;
 571
 572     from0ToAllPoints.setSize(mesh0.nPoints());
 573     from0ToAllPoints = -1;
 574     from1ToAllPoints.setSize(mesh1.nPoints());
 575     from1ToAllPoints = -1;
 576
 577     // Copy coupled points (on cut)
 578     {
 579         const pointField& cutPoints = coupleInfo.cutPoints();
 580
 581         //const labelListList& cutToMasterPoints =
 582         //   coupleInfo.cutToMasterPoints();
 583         labelListList cutToMasterPoints
 584         (
 585             invertOneToMany
 586             (
 587                 cutPoints.size(),
 588                 coupleInfo.masterToCutPoints()
 589             )
 590         );
 591
 592         //const labelListList& cutToSlavePoints =
 593         //    coupleInfo.cutToSlavePoints();
 594         labelListList cutToSlavePoints
 595         (
 596             invertOneToMany
 597             (
 598                 cutPoints.size(),
 599                 coupleInfo.slaveToCutPoints()
 600             )
 601         );
 602
 603         forAll(cutPoints, i)
 604         {
 605             allPoints[allPointI] = cutPoints[i];
 606
 607             // Mark all master and slave points referring to this point.
 608
 609             const labelList& masterPoints = cutToMasterPoints[i];
 610
 611             forAll(masterPoints, j)
 612             {
 613                 label mesh0PointI = masterPatch.meshPoints()[masterPoints[j]];
 614                 from0ToAllPoints[mesh0PointI] = allPointI;
 615             }
 616
 617             const labelList& slavePoints = cutToSlavePoints[i];
 618
 619             forAll(slavePoints, j)
 620             {
 621                 label mesh1PointI = slavePatch.meshPoints()[slavePoints[j]];
 622                 from1ToAllPoints[mesh1PointI] = allPointI;
 623             }
 624             allPointI++;
 625         }
 626     }
 627
 628     // Add uncoupled mesh0 points
 629     forAll(mesh0.points(), pointI)
 630     {
 631         if (from0ToAllPoints[pointI] == -1)
 632         {
 633             allPoints[allPointI] = mesh0.points()[pointI];
 634             from0ToAllPoints[pointI] = allPointI;
 635             allPointI++;
 636         }
 637     }
 638
 639     // Add uncoupled mesh1 points
 640     forAll(mesh1.points(), pointI)
 641     {
 642         if (from1ToAllPoints[pointI] == -1)
 643         {
 644             allPoints[allPointI] = mesh1.points()[pointI];
 645             from1ToAllPoints[pointI] = allPointI;
 646             allPointI++;
 647         }
 648     }
 649
 650     allPoints.setSize(allPointI);
 651
 652
 653     // Faces
 654     // ~~~~~
 655
 656     // Sizes per patch
 657     nFacesPerPatch.setSize(nAllPatches);
 658     nFacesPerPatch = 0;
 659
 660     // Storage for faces and owner/neighbour
 661     allFaces.setSize(mesh0.nFaces() + mesh1.nFaces());
 662     allOwner.setSize(allFaces.size());
 663     allOwner = -1;
 664     allNeighbour.setSize(allFaces.size());
 665     allNeighbour = -1;
 666     label allFaceI = 0;
 667
 668     from0ToAllFaces.setSize(mesh0.nFaces());
 669     from0ToAllFaces = -1;
 670     from1ToAllFaces.setSize(mesh1.nFaces());
 671     from1ToAllFaces = -1;
 672
 673     // Copy mesh0 internal faces (always uncoupled)
 674     for (label faceI = 0; faceI < mesh0.nInternalFaces(); faceI++)
 675     {
 676         allFaces[allFaceI] = renumber(from0ToAllPoints, mesh0.faces()[faceI]);
 677         allOwner[allFaceI] = mesh0.faceOwner()[faceI];
 678         allNeighbour[allFaceI] = mesh0.faceNeighbour()[faceI];
 679         from0ToAllFaces[faceI] = allFaceI++;
 680     }
 681
 682     // Copy coupled faces. Every coupled face has an equivalent master and
 683     // slave. Also uncount as boundary faces all the newly coupled faces.
 684     const labelList& cutToMasterFaces = coupleInfo.cutToMasterFaces();
 685     const labelList& cutToSlaveFaces = coupleInfo.cutToSlaveFaces();
 686
 687     forAll(cutFaces, i)
 688     {
 689         label masterFaceI = cutToMasterFaces[i];
 690
 691         label mesh0FaceI = masterPatch.addressing()[masterFaceI];
 692
 693         if (from0ToAllFaces[mesh0FaceI] == -1)
 694         {
 695             // First occurrence of face
 696             from0ToAllFaces[mesh0FaceI] = allFaceI;
 697
 698             // External face becomes internal so uncount
 699             label patch0 = patches0.whichPatch(mesh0FaceI);
 700             nFacesPerPatch[patch0]--;
 701         }
 702
 703         label slaveFaceI = cutToSlaveFaces[i];
 704
 705         label mesh1FaceI = slavePatch.addressing()[slaveFaceI];
 706
 707         if (from1ToAllFaces[mesh1FaceI] == -1)
 708         {
 709             from1ToAllFaces[mesh1FaceI] = allFaceI;
 710
 711             label patch1 = patches1.whichPatch(mesh1FaceI);
 712             nFacesPerPatch[from1ToAllPatches[patch1]]--;
 713         }
 714
 715         // Copy cut face (since cutPoints are copied first no renumbering
 716         // nessecary)
 717         allFaces[allFaceI] = cutFaces[i];
 718         allOwner[allFaceI] = mesh0.faceOwner()[mesh0FaceI];
 719         allNeighbour[allFaceI] = mesh1.faceOwner()[mesh1FaceI] + mesh0.nCells();
 720
 721         allFaceI++;
 722     }
 723
 724     // Copy mesh1 internal faces (always uncoupled)
 725     for (label faceI = 0; faceI < mesh1.nInternalFaces(); faceI++)
 726     {
 727         allFaces[allFaceI] = renumber(from1ToAllPoints, mesh1.faces()[faceI]);
 728         allOwner[allFaceI] = mesh1.faceOwner()[faceI] + mesh0.nCells();
 729         allNeighbour[allFaceI] = mesh1.faceNeighbour()[faceI] + mesh0.nCells();
 730         from1ToAllFaces[faceI] = allFaceI++;
 731     }
 732
 733     nInternalFaces = allFaceI;
 734
 735     // Copy (unmarked/uncoupled) external faces in new order.
 736     for (label allPatchI = 0; allPatchI < nAllPatches; allPatchI++)
 737     {
 738         if (allPatchI < patches0.size())
 739         {
 740             // Patch is present in mesh0
 741             const polyPatch& pp = patches0[allPatchI];
 742
 743             nFacesPerPatch[allPatchI] += pp.size();
 744
 745             label faceI = pp.start();
 746
 747             forAll(pp, i)
 748             {
 749                 if (from0ToAllFaces[faceI] == -1)
 750                 {
 751                     // Is uncoupled face since has not yet been dealt with
 752                     allFaces[allFaceI] = renumber
 753                     (
 754                         from0ToAllPoints,
 755                         mesh0.faces()[faceI]
 756                     );
 757                     allOwner[allFaceI] = mesh0.faceOwner()[faceI];
 758                     allNeighbour[allFaceI] = -1;
 759
 760                     from0ToAllFaces[faceI] = allFaceI++;
 761                 }
 762                 faceI++;
 763             }
 764         }
 765         if (fromAllTo1Patches[allPatchI] != -1)
 766         {
 767             // Patch is present in mesh1
 768             const polyPatch& pp = patches1[fromAllTo1Patches[allPatchI]];
 769
 770             nFacesPerPatch[allPatchI] += pp.size();
 771
 772             label faceI = pp.start();
 773
 774             forAll(pp, i)
 775             {
 776                 if (from1ToAllFaces[faceI] == -1)
 777                 {
 778                     // Is uncoupled face
 779                     allFaces[allFaceI] = renumber
 780                     (
 781                         from1ToAllPoints,
 782                         mesh1.faces()[faceI]
 783                     );
 784                     allOwner[allFaceI] =
 785                         mesh1.faceOwner()[faceI]
 786                       + mesh0.nCells();
 787                     allNeighbour[allFaceI] = -1;
 788
 789                     from1ToAllFaces[faceI] = allFaceI++;
 790                 }
 791                 faceI++;
 792             }
 793         }
 794     }
 795     allFaces.setSize(allFaceI);
 796     allOwner.setSize(allFaceI);
 797     allNeighbour.setSize(allFaceI);
 798
 799
 800     // So now we have all ok for one-to-one mapping.
 801     // For split slace faces:
 802     // - mesh consistent with slave side
 803     // - mesh not consistent with owner side. It is not zipped up, the
 804     //   original faces need edges split.
 805
 806     // Use brute force to prevent having to calculate addressing:
 807     // - get map from master edge to split edges.
 808     // - check all faces to find any edge that is split.
 809     {
 810         // From two cut-points to labels of cut-points inbetween.
 811         // (in order: from e[0] to e[1]
 812         const edgeLookup& cutEdgeToPoints = coupleInfo.cutEdgeToPoints();
 813
 814         // Get map of master face (in mesh labels) that are in cut. These faces
 815         // do not need to be renumbered.
 816         labelHashSet masterCutFaces(cutToMasterFaces.size());
 817         forAll(cutToMasterFaces, i)
 818         {
 819             label meshFaceI = masterPatch.addressing()[cutToMasterFaces[i]];
 820
 821             masterCutFaces.insert(meshFaceI);
 822         }
 823
 824         DynamicList<label> workFace(100);
 825
 826         forAll(from0ToAllFaces, face0)
 827         {
 828             if (!masterCutFaces.found(face0))
 829             {
 830                 label allFaceI = from0ToAllFaces[face0];
 831
 832                 insertVertices
 833                 (
 834                     cutEdgeToPoints,
 835                     masterPatch.meshPointMap(),
 836                     coupleInfo.masterToCutPoints(),
 837                     mesh0.faces()[face0],
 838
 839                     workFace,
 840                     allFaces[allFaceI]
 841                 );
 842             }
 843         }
 844
 845         // Same for slave face
 846
 847         labelHashSet slaveCutFaces(cutToSlaveFaces.size());
 848         forAll(cutToSlaveFaces, i)
 849         {
 850             label meshFaceI = slavePatch.addressing()[cutToSlaveFaces[i]];
 851
 852             slaveCutFaces.insert(meshFaceI);
 853         }
 854
 855         forAll(from1ToAllFaces, face1)
 856         {
 857             if (!slaveCutFaces.found(face1))
 858             {
 859                 label allFaceI = from1ToAllFaces[face1];
 860
 861                 insertVertices
 862                 (
 863                     cutEdgeToPoints,
 864                     slavePatch.meshPointMap(),
 865                     coupleInfo.slaveToCutPoints(),
 866                     mesh1.faces()[face1],
 867
 868                     workFace,
 869                     allFaces[allFaceI]
 870                 );
 871             }
 872         }
 873     }
 874
 875     // Now we have a full facelist and owner/neighbour addressing.
 876
 877
 878     // Cells
 879     // ~~~~~
 880
 881     from1ToAllCells.setSize(mesh1.nCells());
 882     from1ToAllCells = -1;
 883
 884     forAll(mesh1.cells(), i)
 885     {
 886         from1ToAllCells[i] = i + mesh0.nCells();
 887     }
 888
 889     // Make cells (= cell-face addressing)
 890     nCells = mesh0.nCells() + mesh1.nCells();
 891     cellList allCells(nCells);
 892     primitiveMesh::calcCells(allCells, allOwner, allNeighbour, nCells);
 893
 894     // Reorder faces for upper-triangular order.
 895     labelList oldToNew
 896     (
 897         getFaceOrder
 898         (
 899             allCells,
 900             nInternalFaces,
 901             allOwner,
 902             allNeighbour
 903         )
 904     );
 905
 906     inplaceReorder(oldToNew, allFaces);
 907     inplaceReorder(oldToNew, allOwner);
 908     inplaceReorder(oldToNew, allNeighbour);
 909     inplaceRenumber(oldToNew, from0ToAllFaces);
 910     inplaceRenumber(oldToNew, from1ToAllFaces);
 911 }
 912
 913
 914 void Foam::polyMeshAdder::mergePointZones
 915 (
 916     const pointZoneMesh& pz0,
 917     const pointZoneMesh& pz1,
 918     const labelList& from0ToAllPoints,
 919     const labelList& from1ToAllPoints,
 920
 921     DynamicList<word>& zoneNames,
 922     labelList& from1ToAll,
 923     List<DynamicList<label> >& pzPoints
 924 )
 925 {
 926     zoneNames.setCapacity(pz0.size() + pz1.size());
 927
 928     // Names
 929     append(pz0.names(), zoneNames);
 930
 931     from1ToAll.setSize(pz1.size());
 932
 933     forAll (pz1, zoneI)
 934     {
 935         from1ToAll[zoneI] = zoneIndex(pz1[zoneI].name(), zoneNames);
 936     }
 937     zoneNames.shrink();
 938
 939     // Point labels per merged zone
 940     pzPoints.setSize(zoneNames.size());
 941
 942     forAll(pz0, zoneI)
 943     {
 944         append(from0ToAllPoints, pz0[zoneI], pzPoints[zoneI]);
 945     }
 946
 947     // Get sorted zone contents for duplicate element recognition
 948     PtrList<SortableList<label> > pzPointsSorted(pzPoints.size());
 949     forAll(pzPoints, zoneI)
 950     {
 951         pzPointsSorted.set
 952         (
 953             zoneI,
 954             new SortableList<label>(pzPoints[zoneI])
 955         );
 956     }
 957
 958     // Now we have full addressing for points so do the pointZones of mesh1.
 959     forAll(pz1, zoneI)
 960     {
 961         // Relabel all points of zone and add to correct pzPoints.
 962         label allZoneI = from1ToAll[zoneI];
 963
 964         append
 965         (
 966             from1ToAllPoints,
 967             pz1[zoneI],
 968             pzPointsSorted[allZoneI],
 969             pzPoints[allZoneI]
 970         );
 971     }
 972
 973     forAll(pzPoints, i)
 974     {
 975         pzPoints[i].shrink();
 976     }
 977 }
 978
 979
 980 void Foam::polyMeshAdder::mergeFaceZones
 981 (
 982     const faceZoneMesh& fz0,
 983     const faceZoneMesh& fz1,
 984     const labelList& from0ToAllFaces,
 985     const labelList& from1ToAllFaces,
 986
 987     DynamicList<word>& zoneNames,
 988     labelList& from1ToAll,
 989     List<DynamicList<label> >& fzFaces,
 990     List<DynamicList<bool> >& fzFlips
 991 )
 992 {
 993     zoneNames.setCapacity(fz0.size() + fz1.size());
 994
 995     append(fz0.names(), zoneNames);
 996
 997     from1ToAll.setSize(fz1.size());
 998
 999     forAll (fz1, zoneI)
1000     {
1001         from1ToAll[zoneI] = zoneIndex(fz1[zoneI].name(), zoneNames);
1002     }
1003     zoneNames.shrink();
1004
1005
1006     // Create temporary lists for faceZones.
1007     fzFaces.setSize(zoneNames.size());
1008     fzFlips.setSize(zoneNames.size());
1009     forAll(fz0, zoneI)
1010     {
1011         DynamicList<label>& newZone = fzFaces[zoneI];
1012         DynamicList<bool>& newFlip = fzFlips[zoneI];
1013
1014         newZone.setCapacity(fz0[zoneI].size());
1015         newFlip.setCapacity(newZone.size());
1016
1017         const labelList& addressing = fz0[zoneI];
1018         const boolList& flipMap = fz0[zoneI].flipMap();
1019
1020         forAll(addressing, i)
1021         {
1022             label faceI = addressing[i];
1023
1024             if (from0ToAllFaces[faceI] != -1)
1025             {
1026                 newZone.append(from0ToAllFaces[faceI]);
1027                 newFlip.append(flipMap[i]);
1028             }
1029         }
1030     }
1031
1032     // Get sorted zone contents for duplicate element recognition
1033     PtrList<SortableList<label> > fzFacesSorted(fzFaces.size());
1034     forAll(fzFaces, zoneI)
1035     {
1036         fzFacesSorted.set
1037         (
1038             zoneI,
1039             new SortableList<label>(fzFaces[zoneI])
1040         );
1041     }
1042
1043     // Now we have full addressing for faces so do the faceZones of mesh1.
1044     forAll(fz1, zoneI)
1045     {
1046         label allZoneI = from1ToAll[zoneI];
1047
1048         DynamicList<label>& newZone = fzFaces[allZoneI];
1049         const SortableList<label>& newZoneSorted = fzFacesSorted[allZoneI];
1050         DynamicList<bool>& newFlip = fzFlips[allZoneI];
1051
1052         newZone.setCapacity(newZone.size() + fz1[zoneI].size());
1053         newFlip.setCapacity(newZone.size());
1054
1055         const labelList& addressing = fz1[zoneI];
1056         const boolList& flipMap = fz1[zoneI].flipMap();
1057
1058         forAll(addressing, i)
1059         {
1060             label faceI = addressing[i];
1061             label allFaceI = from1ToAllFaces[faceI];
1062
1063             if
1064             (
1065                 allFaceI != -1
1066              && findSortedIndex(newZoneSorted, allFaceI) == -1
1067             )
1068             {
1069                 newZone.append(allFaceI);
1070                 newFlip.append(flipMap[i]);
1071             }
1072         }
1073     }
1074
1075     forAll(fzFaces, i)
1076     {
1077         fzFaces[i].shrink();
1078         fzFlips[i].shrink();
1079     }
1080 }
1081
1082
1083 void Foam::polyMeshAdder::mergeCellZones
1084 (
1085     const cellZoneMesh& cz0,
1086     const cellZoneMesh& cz1,
1087     const labelList& from1ToAllCells,
1088
1089     DynamicList<word>& zoneNames,
1090     labelList& from1ToAll,
1091     List<DynamicList<label> >& czCells
1092 )
1093 {
1094     zoneNames.setCapacity(cz0.size() + cz1.size());
1095
1096     append(cz0.names(), zoneNames);
1097
1098     from1ToAll.setSize(cz1.size());
1099     forAll (cz1, zoneI)
1100     {
1101         from1ToAll[zoneI] = zoneIndex(cz1[zoneI].name(), zoneNames);
1102     }
1103     zoneNames.shrink();
1104
1105
1106     // Create temporary lists for cellZones.
1107     czCells.setSize(zoneNames.size());
1108     forAll(cz0, zoneI)
1109     {
1110         // Insert mesh0 cells
1111         append(cz0[zoneI], czCells[zoneI]);
1112     }
1113
1114
1115     // Cell mapping is trivial.
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             xferMove(allPoints),
1450             xferMove(allFaces),
1451             xferMove(allOwner),
1452             xferMove(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         xferMove(allPoints),
1754         xferMove(allFaces),
1755         xferMove(allOwner),
1756         xferMove(allNeighbour),
1757         patchSizes,     // size
1758         patchStarts,    // patchstarts
1759         validBoundary   // boundary valid?
1760     );
1761
1762     // Add zones to new mesh.
1763     mesh0.pointZones().clear();
1764     mesh0.faceZones().clear();
1765     mesh0.cellZones().clear();
1766     addZones
1767     (
1768         pointZoneNames,
1769         pzPoints,
1770
1771         faceZoneNames,
1772         fzFaces,
1773         fzFlips,
1774
1775         cellZoneNames,
1776         czCells,
1777         mesh0
1778     );
1779
1780     return mapPtr;
1781 }
1782
1783
1784 Foam::Map<Foam::label> Foam::polyMeshAdder::findSharedPoints
1785 (
1786     const polyMesh& mesh,
1787     const scalar mergeDist
1788 )
1789 {
1790     const labelList& sharedPointLabels = mesh.globalData().sharedPointLabels();
1791     const labelList& sharedPointAddr = mesh.globalData().sharedPointAddr();
1792
1793     // Because of adding the missing pieces e.g. when redistributing a mesh
1794     // it can be that there are multiple points on the same processor that
1795     // refer to the same shared point.
1796
1797     // Invert point-to-shared addressing
1798     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1799
1800     Map<labelList> sharedToMesh(sharedPointLabels.size());
1801
1802     label nMultiple = 0;
1803
1804     forAll(sharedPointLabels, i)
1805     {
1806         label pointI = sharedPointLabels[i];
1807
1808         label sharedI = sharedPointAddr[i];
1809
1810         Map<labelList>::iterator iter = sharedToMesh.find(sharedI);
1811
1812         if (iter != sharedToMesh.end())
1813         {
1814             // sharedI already used by other point. Add this one.
1815
1816             nMultiple++;
1817
1818             labelList& connectedPointLabels = iter();
1819
1820             label sz = connectedPointLabels.size();
1821
1822             // Check just to make sure.
1823             if (findIndex(connectedPointLabels, pointI) != -1)
1824             {
1825                 FatalErrorIn("polyMeshAdder::findSharedPoints(..)")
1826                     << "Duplicate point in sharedPoint addressing." << endl
1827                     << "When trying to add point " << pointI << " on shared "
1828                     << sharedI  << " with connected points "
1829                     << connectedPointLabels
1830                     << abort(FatalError);
1831             }
1832
1833             connectedPointLabels.setSize(sz+1);
1834             connectedPointLabels[sz] = pointI;
1835         }
1836         else
1837         {
1838             sharedToMesh.insert(sharedI, labelList(1, pointI));
1839         }
1840     }
1841
1842
1843     // Assign single master for every shared with multiple geometric points
1844     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1845
1846     Map<label> pointToMaster(nMultiple);
1847
1848     forAllConstIter(Map<labelList>, sharedToMesh, iter)
1849     {
1850         const labelList& connectedPointLabels = iter();
1851
1852         //Pout<< "For shared:" << iter.key()
1853         //    << " found points:" << connectedPointLabels
1854         //    << " at coords:"
1855         //    <<  pointField(mesh.points(), connectedPointLabels) << endl;
1856
1857         if (connectedPointLabels.size() > 1)
1858         {
1859             const pointField connectedPoints
1860             (
1861                 mesh.points(),
1862                 connectedPointLabels
1863             );
1864
1865             labelList toMergedPoints;
1866             pointField mergedPoints;
1867             bool hasMerged = Foam::mergePoints
1868             (
1869                 connectedPoints,
1870                 mergeDist,
1871                 false,
1872                 toMergedPoints,
1873                 mergedPoints
1874             );
1875
1876             if (hasMerged)
1877             {
1878                 // Invert toMergedPoints
1879                 const labelListList mergeSets
1880                 (
1881                     invertOneToMany
1882                     (
1883                         mergedPoints.size(),
1884                         toMergedPoints
1885                     )
1886                 );
1887
1888                 // Find master for valid merges
1889                 forAll(mergeSets, setI)
1890                 {
1891                     const labelList& mergeSet = mergeSets[setI];
1892
1893                     if (mergeSet.size() > 1)
1894                     {
1895                         // Pick lowest numbered point
1896                         label masterPointI = labelMax;
1897
1898                         forAll(mergeSet, i)
1899                         {
1900                             label pointI = connectedPointLabels[mergeSet[i]];
1901
1902                             masterPointI = min(masterPointI, pointI);
1903                         }
1904
1905                         forAll(mergeSet, i)
1906                         {
1907                             label pointI = connectedPointLabels[mergeSet[i]];
1908
1909                             //Pout<< "Merging point " << pointI
1910                             //    << " at " << mesh.points()[pointI]
1911                             //    << " into master point "
1912                             //    << masterPointI
1913                             //    << " at " << mesh.points()[masterPointI]
1914                             //    << endl;
1915
1916                             pointToMaster.insert(pointI, masterPointI);
1917                         }
1918                     }
1919                 }
1920             }
1921         }
1922     }
1923
1924     //- Old: geometric merging. Causes problems for two close shared points.
1925     //labelList sharedToMerged;
1926     //pointField mergedPoints;
1927     //bool hasMerged = Foam::mergePoints
1928     //(
1929     //    pointField
1930     //    (
1931     //        mesh.points(),
1932     //        sharedPointLabels
1933     //    ),
1934     //    mergeDist,
1935     //    false,
1936     //    sharedToMerged,
1937     //    mergedPoints
1938     //);
1939     //
1940     //// Find out which sets of points get merged and create a map from
1941     //// mesh point to unique point.
1942     //
1943     //Map<label> pointToMaster(10*sharedToMerged.size());
1944     //
1945     //if (hasMerged)
1946     //{
1947     //    labelListList mergeSets
1948     //    (
1949     //        invertOneToMany
1950     //        (
1951     //            sharedToMerged.size(),
1952     //            sharedToMerged
1953     //        )
1954     //    );
1955     //
1956     //    label nMergeSets = 0;
1957     //
1958     //    forAll(mergeSets, setI)
1959     //    {
1960     //        const labelList& mergeSet = mergeSets[setI];
1961     //
1962     //        if (mergeSet.size() > 1)
1963     //        {
1964     //            // Take as master the shared point with the lowest mesh
1965     //            // point label. (rather arbitrarily - could use max or
1966     //            // any other one of the points)
1967     //
1968     //            nMergeSets++;
1969     //
1970     //            label masterI = labelMax;
1971     //
1972     //            forAll(mergeSet, i)
1973     //            {
1974     //                label sharedI = mergeSet[i];
1975     //
1976     //                masterI = min(masterI, sharedPointLabels[sharedI]);
1977     //            }
1978     //
1979     //            forAll(mergeSet, i)
1980     //            {
1981     //                label sharedI = mergeSet[i];
1982     //
1983     //                pointToMaster.insert(sharedPointLabels[sharedI], masterI);
1984     //            }
1985     //        }
1986     //    }
1987     //
1988     //    //if (debug)
1989     //    //{
1990     //    //    Pout<< "polyMeshAdder : merging:"
1991     //    //        << pointToMaster.size() << " into " << nMergeSets
1992     //    //        << " sets." << endl;
1993     //    //}
1994     //}
1995
1996     return pointToMaster;
1997 }
1998
1999
2000 void Foam::polyMeshAdder::mergePoints
2001 (
2002     const polyMesh& mesh,
2003     const Map<label>& pointToMaster,
2004     polyTopoChange& meshMod
2005 )
2006 {
2007     // Remove all non-master points.
2008     forAll(mesh.points(), pointI)
2009     {
2010         Map<label>::const_iterator iter = pointToMaster.find(pointI);
2011
2012         if (iter != pointToMaster.end())
2013         {
2014             if (iter() != pointI)
2015             {
2016                 meshMod.removePoint(pointI, iter());
2017             }
2018         }
2019     }
2020
2021     // Modify faces for points. Note: could use pointFaces here but want to
2022     // avoid addressing calculation.
2023     const faceList& faces = mesh.faces();
2024
2025     forAll(faces, faceI)
2026     {
2027         const face& f = faces[faceI];
2028
2029         bool hasMerged = false;
2030
2031         forAll(f, fp)
2032         {
2033             label pointI = f[fp];
2034
2035             Map<label>::const_iterator iter = pointToMaster.find(pointI);
2036
2037             if (iter != pointToMaster.end())
2038             {
2039                 if (iter() != pointI)
2040                 {
2041                     hasMerged = true;
2042                     break;
2043                 }
2044             }
2045         }
2046
2047         if (hasMerged)
2048         {
2049             face newF(f);
2050
2051             forAll(f, fp)
2052             {
2053                 label pointI = f[fp];
2054
2055                 Map<label>::const_iterator iter = pointToMaster.find(pointI);
2056
2057                 if (iter != pointToMaster.end())
2058                 {
2059                     newF[fp] = iter();
2060                 }
2061             }
2062
2063             label patchID = mesh.boundaryMesh().whichPatch(faceI);
2064             label nei = (patchID == -1 ? mesh.faceNeighbour()[faceI] : -1);
2065             label zoneID = mesh.faceZones().whichZone(faceI);
2066             bool zoneFlip = false;
2067
2068             if (zoneID >= 0)
2069             {
2070                 const faceZone& fZone = mesh.faceZones()[zoneID];
2071                 zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
2072             }
2073
2074             meshMod.setAction
2075             (
2076                 polyModifyFace
2077                 (
2078                     newF,                       // modified face
2079                     faceI,                      // label of face
2080                     mesh.faceOwner()[faceI],    // owner
2081                     nei,                        // neighbour
2082                     false,                      // face flip
2083                     patchID,                    // patch for face
2084                     false,                      // remove from zone
2085                     zoneID,                     // zone for face
2086                     zoneFlip                    // face flip in zone
2087                 )
2088             );
2089         }
2090     }
2091 }
2092
2093
2094 // ************************************************************************* //