applications/utilities/mesh/manipulation/mergeMeshes/mergePolyMesh.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   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 "mergePolyMesh.H"
  28 #include "Time.H"
  29 #include "directTopoChange.H"
  30 #include "mapPolyMesh.H"
  31 #include "polyAddPoint.H"
  32 #include "polyAddCell.H"
  33 #include "polyAddFace.H"
  34
  35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  36
  37 defineTypeNameAndDebug(Foam::mergePolyMesh, 1);
  38
  39
  40 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  41
  42 Foam::label Foam::mergePolyMesh::patchIndex(const polyPatch& p)
  43 {
  44     // Find the patch name on the list.  If the patch is already there
  45     // and patch types match, return index
  46     const word& pType = p.type();
  47     const word& pName = p.name();
  48
  49     bool nameFound = false;
  50
  51     forAll (patchNames_, patchI)
  52     {
  53         if (patchNames_[patchI] == pName)
  54         {
  55             if (patchTypes_[patchI] == pType)
  56             {
  57                 // Found name and types match
  58                 return patchI;
  59             }
  60             else
  61             {
  62                 // Found the name, but type is different
  63                 nameFound = true;
  64             }
  65         }
  66     }
  67
  68     // Patch not found.  Append to the list
  69     patchTypes_.append(pType);
  70
  71     if (nameFound)
  72     {
  73         // Duplicate name is not allowed.  Create a composite name from the
  74         // patch name and case name
  75         const word& caseName = p.boundaryMesh().mesh().time().caseName();
  76
  77         patchNames_.append(pName + "_" + caseName);
  78
  79         Info<< "label patchIndex(const polyPatch& p) : "
  80             << "Patch " << p.index() << " named "
  81             << pName << " in mesh " << caseName
  82             << " already exists, but patch types "
  83             << " do not match.\nCreating a composite name as "
  84             << patchNames_[patchNames_.size() - 1] << endl;
  85     }
  86     else
  87     {
  88         patchNames_.append(pName);
  89     }
  90
  91     return patchNames_.size() - 1;
  92 }
  93
  94
  95 Foam::label Foam::mergePolyMesh::zoneIndex
  96 (
  97     DynamicList<word>& names,
  98     const word& curName
  99 )
 100 {
 101     forAll (names, zoneI)
 102     {
 103         if (names[zoneI] == curName)
 104         {
 105             return zoneI;
 106         }
 107     }
 108
 109     // Not found.  Add new name to the list
 110     names.append(curName);
 111
 112     return names.size() - 1;
 113 }
 114
 115
 116 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 117
 118 Foam::mergePolyMesh::mergePolyMesh(const IOobject& io)
 119 :
 120     polyMesh(io),
 121     meshMod_(*this),
 122     patchTypes_(2*boundaryMesh().size()),
 123     patchNames_(2*boundaryMesh().size()),
 124     pointZoneNames_(),
 125     faceZoneNames_(),
 126     cellZoneNames_(),
 127     pointZones_(),
 128     faceZones_(),
 129     faceZoneFlips_(),
 130     cellZones_()
 131 {
 132     // Insert the original patches into the list
 133     wordList curPatchTypes = boundaryMesh().types();
 134     wordList curPatchNames = boundaryMesh().names();
 135
 136     forAll (curPatchTypes, patchI)
 137     {
 138         patchTypes_.append(curPatchTypes[patchI]);
 139         patchNames_.append(curPatchNames[patchI]);
 140     }
 141
 142     // Insert point, face and cell zones into the list
 143
 144     // Point zones
 145     wordList curPointZoneNames = pointZones().names();
 146     if (curPointZoneNames.size())
 147     {
 148         pointZoneNames_.setCapacity(2*curPointZoneNames.size());
 149     }
 150
 151     forAll (curPointZoneNames, zoneI)
 152     {
 153         pointZoneNames_.append(curPointZoneNames[zoneI]);
 154     }
 155     const pointField& p = points();
 156     const pointZoneMesh& pz = pointZones();
 157     forAll(p, pointI)
 158     {
 159         pointZones_.append(pz.whichZone(pointI));
 160     }
 161
 162     // Face zones
 163     wordList curFaceZoneNames = faceZones().names();
 164
 165     if (curFaceZoneNames.size())
 166     {
 167         faceZoneNames_.setCapacity(2*curFaceZoneNames.size());
 168     }
 169     forAll (curFaceZoneNames, zoneI)
 170     {
 171         faceZoneNames_.append(curFaceZoneNames[zoneI]);
 172     }
 173     const faceList& f = faces();
 174     const faceZoneMesh& fz = faceZones();
 175     forAll(f, faceI)
 176     {
 177         label zoneID = fz.whichZone(faceI);
 178         faceZones_.append(zoneID);
 179
 180         if(zoneID > -1)
 181             faceZoneFlips_.append(fz[zoneID].flipMap()[faceI]);
 182         else
 183             faceZoneFlips_.append(false);
 184     }
 185     faceZones_.shrink();
 186     faceZoneFlips_.shrink();
 187
 188     // Cell zones
 189     wordList curCellZoneNames = cellZones().names();
 190
 191     if (curCellZoneNames.size())
 192     {
 193         cellZoneNames_.setCapacity(2*curCellZoneNames.size());
 194     }
 195     forAll (curCellZoneNames, zoneI)
 196     {
 197         cellZoneNames_.append(curCellZoneNames[zoneI]);
 198     }
 199     const cellList& c = cells();
 200     const cellZoneMesh& cz = cellZones();
 201     forAll(c, cellI)
 202     {
 203         cellZones_.append(cz.whichZone(cellI));
 204     }
 205     cellZones_.shrink();
 206 }
 207
 208
 209 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 210
 211
 212 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 213
 214 void Foam::mergePolyMesh::addMesh(const polyMesh& m)
 215 {
 216     // Add all the points, faces and cells of the new mesh
 217
 218     // Add points
 219
 220     label zoneID = -1;
 221
 222     const pointField& p = m.points();
 223     labelList renumberPoints(p.size());
 224
 225     const pointZoneMesh& pz = m.pointZones();
 226     labelList pointZoneIndices(pz.size());
 227
 228     forAll (pz, zoneI)
 229     {
 230         pointZoneIndices[zoneI] = zoneIndex(pointZoneNames_, pz[zoneI].name());
 231     }
 232
 233     forAll (p, pointI)
 234     {
 235         // Grab zone ID.  If a point is not in a zone, it will return -1
 236         zoneID = pz.whichZone(pointI);
 237
 238         if (zoneID > -1)
 239         {
 240             // Translate zone ID into the new index
 241             zoneID = pointZoneIndices[zoneID];
 242         }
 243
 244         renumberPoints[pointI] =
 245             meshMod_.setAction
 246             (
 247                 polyAddPoint
 248                 (
 249                     p[pointI],            // Point to add
 250                     -1,                   // Master point (straight addition)
 251                     zoneID,               // Zone for point
 252                     pointI < m.nPoints()  // Is in cell?
 253                 )
 254             );
 255
 256         pointZones_.append(zoneID);
 257
 258     }
 259
 260     // Add cells
 261
 262     const cellList& c = m.cells();
 263     labelList renumberCells(c.size());
 264
 265     const cellZoneMesh& cz = m.cellZones();
 266     labelList cellZoneIndices(cz.size());
 267
 268     forAll (cz, zoneI)
 269     {
 270         cellZoneIndices[zoneI] = zoneIndex(cellZoneNames_, cz[zoneI].name());
 271     }
 272
 273     forAll (c, cellI)
 274     {
 275         // Grab zone ID.  If a cell is not in a zone, it will return -1
 276         zoneID = cz.whichZone(cellI);
 277
 278         if (zoneID > -1)
 279         {
 280             // Translate zone ID into the new index
 281             zoneID = cellZoneIndices[zoneID];
 282         }
 283
 284         renumberCells[cellI] =
 285             meshMod_.setAction
 286             (
 287                 polyAddCell
 288                 (
 289                     -1,                   // Master point
 290                     -1,                   // Master edge
 291                     -1,                   // Master face
 292                     -1,                   // Master cell
 293                     zoneID                // Zone for cell
 294                 )
 295             );
 296
 297         cellZones_.append(zoneID);
 298     }
 299
 300     // Add faces
 301     const polyBoundaryMesh& bm = m.boundaryMesh();
 302
 303     // Gather the patch indices
 304     labelList patchIndices(bm.size());
 305
 306     forAll (patchIndices, patchI)
 307     {
 308         patchIndices[patchI] = patchIndex(bm[patchI]);
 309     }
 310
 311     // Temporary: update number of allowable patches. This should be
 312     // determined at the top - before adding anything.
 313     meshMod_.setNumPatches(patchNames_.size());
 314
 315
 316
 317     const faceZoneMesh& fz = m.faceZones();
 318     labelList faceZoneIndices(fz.size());
 319
 320     forAll (fz, zoneI)
 321     {
 322         faceZoneIndices[zoneI] = zoneIndex(faceZoneNames_, fz[zoneI].name());
 323     }
 324
 325     const faceList& f = m.faces();
 326     labelList renumberFaces(f.size());
 327
 328     const labelList& own = m.faceOwner();
 329     const labelList& nei = m.faceNeighbour();
 330
 331     label newOwn, newNei, newPatch, newZone;
 332     bool newZoneFlip;
 333
 334     forAll (f, faceI)
 335     {
 336         const face& curFace = f[faceI];
 337
 338         face newFace(curFace.size());
 339
 340         forAll (curFace, pointI)
 341         {
 342             newFace[pointI] = renumberPoints[curFace[pointI]];
 343         }
 344
 345         if (debug)
 346         {
 347             // Check that the face is valid
 348             if (min(newFace) < 0)
 349             {
 350                 FatalErrorIn("void mergePolyMesh::addMesh(const polyMesh&)")
 351                     << "Error in point mapping for face " << faceI
 352                     << ".  Old face: " << curFace << " New face: " << newFace
 353                     << abort(FatalError);
 354             }
 355         }
 356
 357         if (faceI < m.nInternalFaces() || faceI >= m.nFaces())
 358         {
 359             newPatch = -1;
 360         }
 361         else
 362         {
 363             newPatch = patchIndices[bm.whichPatch(faceI)];
 364         }
 365
 366         newOwn = own[faceI];
 367         if (newOwn > -1) newOwn = renumberCells[newOwn];
 368
 369         if (newPatch > -1)
 370         {
 371             newNei = -1;
 372         }
 373         else
 374         {
 375             newNei = nei[faceI];
 376             newNei = renumberCells[newNei];
 377         }
 378
 379
 380         newZone = fz.whichZone(faceI);
 381         newZoneFlip = false;
 382
 383         if (newZone > -1)
 384         {
 385             newZoneFlip = fz[newZone].flipMap()[fz[newZone].whichFace(faceI)];
 386
 387             // Grab the new zone
 388             newZone = faceZoneIndices[newZone];
 389         }
 390
 391         renumberFaces[faceI] =
 392             meshMod_.setAction
 393             (
 394                 polyAddFace
 395                 (
 396                     newFace,
 397                     newOwn,
 398                     newNei,
 399                     -1,
 400                     -1,
 401                     -1,
 402                     false,
 403                     newPatch,
 404                     newZone,
 405                     newZoneFlip
 406                 )
 407             );
 408         faceZones_.append(newZone);
 409         faceZoneFlips_.append(newZoneFlip);
 410     }
 411
 412     faceZones_.shrink();
 413     faceZoneFlips_.shrink();
 414 }
 415
 416
 417 void Foam::mergePolyMesh::merge()
 418 {
 419     Info<< "patch names: " << patchNames_ << nl
 420         << "patch types: " << patchTypes_ << nl
 421         << "point zone names: " << pointZoneNames_ << nl
 422         << "face zone names: " << faceZoneNames_ << nl
 423         << "cell zone names: " << cellZoneNames_ << endl;
 424
 425     // Add the patches if necessary
 426     if (patchNames_.size() != boundaryMesh().size())
 427     {
 428         Info << "Copying old patches" << endl;
 429
 430         List<polyPatch*> newPatches(patchNames_.size());
 431
 432         const polyBoundaryMesh& oldPatches = boundaryMesh();
 433
 434         // Note.  Re-using counter in two for loops
 435         label patchI = 0;
 436
 437         for (patchI = 0; patchI < oldPatches.size(); patchI++)
 438         {
 439             newPatches[patchI] = oldPatches[patchI].clone(oldPatches).ptr();
 440         }
 441
 442         Info << "Adding new patches. " << endl;
 443
 444         label endOfLastPatch =
 445             oldPatches[patchI - 1].start() + oldPatches[patchI - 1].size();
 446
 447         for (; patchI < patchNames_.size(); patchI++)
 448         {
 449             // Add a patch
 450             newPatches[patchI] =
 451             (
 452                 polyPatch::New
 453                 (
 454                     patchTypes_[patchI],
 455                     patchNames_[patchI],
 456                     0,
 457                     endOfLastPatch,
 458                     patchI,
 459                     oldPatches
 460                 ).ptr()
 461             );
 462         }
 463
 464         removeBoundary();
 465         addPatches(newPatches);
 466     }
 467
 468     // Add the zones if necessary
 469     if
 470     (
 471         pointZoneNames_.size() != pointZones().size()
 472      || faceZoneNames_.size() != faceZones().size()
 473      || cellZoneNames_.size() != cellZones().size()
 474     )
 475     {
 476
 477         List<pointZone*> pZones(pointZoneNames_.size());
 478         List<faceZone*> fZones(faceZoneNames_.size());
 479         List<cellZone*> cZones(cellZoneNames_.size());
 480
 481         if(pointZoneNames_.size() > 0)
 482         {
 483             Info << "Updating point zones" << endl;
 484
 485             List<DynamicList<label> > pzPoints(pointZoneNames_.size());
 486             forAll(pointZones_, pointI)
 487             {
 488                 label zoneID = pointZones_[pointI];
 489
 490                 if(zoneID > -1)
 491                     pzPoints[zoneID].append(pointI);
 492             }
 493
 494             forAll(pZones, i)
 495             {
 496                 pZones[i] = new pointZone
 497                     (
 498                         pointZoneNames_[i],
 499                         pzPoints[i],
 500                         i,
 501                         pointZones()
 502                     );
 503             }
 504         }
 505
 506         if(faceZoneNames_.size() > 0)
 507         {
 508             Info << "Updating face zones" << endl;
 509
 510             List<DynamicList<label> > fzFaces(faceZoneNames_.size());
 511             List<DynamicList<bool> >  fzFlips(faceZoneNames_.size());
 512
 513             forAll(faceZones_, faceI)
 514             {
 515                 label zoneID = faceZones_[faceI];
 516
 517                 if(zoneID > -1)
 518                 {
 519                     fzFaces[zoneID].append(faceI);
 520                     fzFlips[zoneID].append(faceZoneFlips_[faceI]);
 521                 }
 522             }
 523
 524             forAll(fZones, i)
 525             {
 526                 fzFaces[i].shrink();
 527                 fzFlips[i].shrink();
 528
 529                 fZones[i] = new faceZone
 530                     (
 531                         faceZoneNames_[i],
 532                         fzFaces[i],
 533                         fzFlips[i],
 534                         i,
 535                         faceZones()
 536                     );
 537             }
 538         }
 539
 540         if(cellZoneNames_.size() > 0)
 541         {
 542             Info << "Updating cell zones" << endl;
 543
 544             List<DynamicList<label> > czCells(cellZoneNames_.size());
 545             forAll(cellZones_, cellI)
 546             {
 547                 label zoneID = cellZones_[cellI];
 548
 549                 if(zoneID > -1)
 550                     czCells[zoneID].append(cellI);
 551             }
 552
 553             forAll(cZones, i)
 554             {
 555                 czCells[i].shrink();
 556
 557                 cZones[i] = new cellZone
 558                     (
 559                         cellZoneNames_[i],
 560                         czCells[i],
 561                         i,
 562                         cellZones()
 563                     );
 564             }
 565         }
 566
 567         removeZones();
 568         addZones(pZones, fZones, cZones);
 569         write();
 570     }
 571
 572     // Change mesh. No inflation
 573     meshMod_.changeMesh(*this, false);
 574
 575     // Clear topo change for the next operation
 576     meshMod_.clear();
 577 }
 578
 579
 580 // ************************************************************************* //