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