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