applications/utilities/mesh/advanced/modifyMesh/cellSplitter.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Description
  25
  26 \*---------------------------------------------------------------------------*/
  27
  28 #include "cellSplitter.H"
  29 #include "polyMesh.H"
  30 #include "directTopoChange.H"
  31 #include "polyAddCell.H"
  32 #include "polyAddFace.H"
  33 #include "polyAddPoint.H"
  34 #include "polyModifyFace.H"
  35 #include "mapPolyMesh.H"
  36 #include "meshTools.H"
  37
  38
  39 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  40
  41 namespace Foam
  42 {
  43
  44 defineTypeNameAndDebug(cellSplitter, 0);
  45
  46 }
  47
  48 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  49
  50 void Foam::cellSplitter::getFaceInfo
  51 (
  52     const label faceI,
  53     label& patchID,
  54     label& zoneID,
  55     label& zoneFlip
  56 ) const
  57 {
  58     patchID = -1;
  59
  60     if (!mesh_.isInternalFace(faceI))
  61     {
  62         patchID = mesh_.boundaryMesh().whichPatch(faceI);
  63     }
  64
  65     zoneID = mesh_.faceZones().whichZone(faceI);
  66
  67     zoneFlip = false;
  68
  69     if (zoneID >= 0)
  70     {
  71         const faceZone& fZone = mesh_.faceZones()[zoneID];
  72
  73         zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
  74     }
  75 }
  76
  77
  78 // Find the new owner of faceI (since the original cell has been split into
  79 // newCells
  80 Foam::label Foam::cellSplitter::newOwner
  81 (
  82     const label faceI,
  83     const Map<labelList>& cellToCells
  84 ) const
  85 {
  86     label oldOwn = mesh_.faceOwner()[faceI];
  87
  88     Map<labelList>::const_iterator fnd = cellToCells.find(oldOwn);
  89
  90     if (fnd == cellToCells.end())
  91     {
  92         // Unsplit cell
  93         return oldOwn;
  94     }
  95     else
  96     {
  97         // Look up index of face in the cells' faces.
  98
  99         const labelList& newCells = fnd();
 100
 101         const cell& cFaces = mesh_.cells()[oldOwn];
 102
 103         return newCells[findIndex(cFaces, faceI)];
 104     }
 105 }
 106
 107
 108 Foam::label Foam::cellSplitter::newNeighbour
 109 (
 110     const label faceI,
 111     const Map<labelList>& cellToCells
 112 ) const
 113 {
 114     label oldNbr = mesh_.faceNeighbour()[faceI];
 115
 116     Map<labelList>::const_iterator fnd = cellToCells.find(oldNbr);
 117
 118     if (fnd == cellToCells.end())
 119     {
 120         // Unsplit cell
 121         return oldNbr;
 122     }
 123     else
 124     {
 125         // Look up index of face in the cells' faces.
 126
 127         const labelList& newCells = fnd();
 128
 129         const cell& cFaces = mesh_.cells()[oldNbr];
 130
 131         return newCells[findIndex(cFaces, faceI)];
 132     }
 133 }
 134
 135
 136 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 137
 138 // Construct from components
 139 Foam::cellSplitter::cellSplitter(const polyMesh& mesh)
 140 :
 141     mesh_(mesh),
 142     addedPoints_()
 143 {}
 144
 145
 146 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 147
 148 Foam::cellSplitter::~cellSplitter()
 149 {}
 150
 151
 152 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 153
 154 void Foam::cellSplitter::setRefinement
 155 (
 156     const Map<point>& cellToMidPoint,
 157     directTopoChange& meshMod
 158 )
 159 {
 160     addedPoints_.clear();
 161     addedPoints_.resize(cellToMidPoint.size());
 162
 163
 164     //
 165     // Introduce cellToMidPoints.
 166     //
 167
 168     forAllConstIter(Map<point>, cellToMidPoint, iter)
 169     {
 170         label cellI = iter.key();
 171
 172         label anchorPoint = mesh_.cellPoints()[cellI][0];
 173
 174         label addedPointI =
 175             meshMod.setAction
 176             (
 177                 polyAddPoint
 178                 (
 179                     iter(),         // point
 180                     anchorPoint,    // master point
 181                     -1,             // zone for point
 182                     true            // supports a cell
 183                 )
 184             );
 185         addedPoints_.insert(cellI, addedPointI);
 186
 187         //Pout<< "Added point " << addedPointI
 188         //    << iter() << " in cell " << cellI << " with centre "
 189         //    << mesh_.cellCentres()[cellI] << endl;
 190     }
 191
 192
 193     //
 194     // Add cells (first one is modified original cell)
 195     //
 196
 197     Map<labelList> cellToCells(cellToMidPoint.size());
 198
 199     forAllConstIter(Map<point>, cellToMidPoint, iter)
 200     {
 201         label cellI = iter.key();
 202
 203         const cell& cFaces = mesh_.cells()[cellI];
 204
 205         // Cells created for this cell.
 206         labelList newCells(cFaces.size());
 207
 208         // First pyramid is the original cell
 209         newCells[0] = cellI;
 210
 211         // Add other pyramids
 212         for (label i = 1; i < cFaces.size(); i++)
 213         {
 214             label addedCellI =
 215                 meshMod.setAction
 216                 (
 217                     polyAddCell
 218                     (
 219                         -1,     // master point
 220                         -1,     // master edge
 221                         -1,     // master face
 222                         cellI,  // master cell
 223                         -1      // zone
 224                     )
 225                 );
 226
 227             newCells[i] = addedCellI;
 228         }
 229
 230         cellToCells.insert(cellI, newCells);
 231
 232         //Pout<< "Split cell " << cellI
 233         //    << " with centre " << mesh_.cellCentres()[cellI] << nl
 234         //    << " faces:" << cFaces << nl
 235         //    << " into :" << newCells << endl;
 236     }
 237
 238
 239     //
 240     // Introduce internal faces. These go from edges of the cell to the mid
 241     // point.
 242     //
 243
 244     forAllConstIter(Map<point>, cellToMidPoint, iter)
 245     {
 246         label cellI = iter.key();
 247
 248         label midPointI = addedPoints_[cellI];
 249
 250         const cell& cFaces = mesh_.cells()[cellI];
 251
 252         const labelList& cEdges = mesh_.cellEdges()[cellI];
 253
 254         forAll(cEdges, i)
 255         {
 256             label edgeI = cEdges[i];
 257             const edge& e = mesh_.edges()[edgeI];
 258
 259             // Get the faces on the cell using the edge
 260             label face0, face1;
 261             meshTools::getEdgeFaces(mesh_, cellI, edgeI, face0, face1);
 262
 263             // Get the cells on both sides of the face by indexing into cFaces.
 264             // (since newly created cells are stored in cFaces order)
 265             const labelList& newCells = cellToCells[cellI];
 266
 267             label cell0 = newCells[findIndex(cFaces, face0)];
 268             label cell1 = newCells[findIndex(cFaces, face1)];
 269
 270             if (cell0 < cell1)
 271             {
 272                 // Construct face to midpoint that is pointing away from
 273                 // (pyramid split off from) cellI
 274
 275                 const face& f0 = mesh_.faces()[face0];
 276
 277                 label index = findIndex(f0, e[0]);
 278
 279                 bool edgeInFaceOrder = (f0[f0.fcIndex(index)] == e[1]);
 280
 281                 // Check if cellI is the face owner
 282
 283                 face newF(3);
 284                 if (edgeInFaceOrder == (mesh_.faceOwner()[face0] == cellI))
 285                 {
 286                     // edge used in face order.
 287                     newF[0] = e[1];
 288                     newF[1] = e[0];
 289                     newF[2] = midPointI;
 290                 }
 291                 else
 292                 {
 293                     newF[0] = e[0];
 294                     newF[1] = e[1];
 295                     newF[2] = midPointI;
 296                 }
 297
 298                 // Now newF points away from cell0
 299                 meshMod.setAction
 300                 (
 301                     polyAddFace
 302                     (
 303                         newF,                       // face
 304                         cell0,                      // owner
 305                         cell1,                      // neighbour
 306                         -1,                         // master point
 307                         -1,                         // master edge
 308                         face0,                      // master face for addition
 309                         false,                      // flux flip
 310                         -1,                         // patch for face
 311                         -1,                         // zone for face
 312                         false                       // face zone flip
 313                     )
 314                 );
 315             }
 316             else
 317             {
 318                 // Construct face to midpoint that is pointing away from
 319                 // (pyramid split off from) cellI
 320
 321                 const face& f1 = mesh_.faces()[face1];
 322
 323                 label index = findIndex(f1, e[0]);
 324
 325                 bool edgeInFaceOrder = (f1[f1.fcIndex(index)] == e[1]);
 326
 327                 // Check if cellI is the face owner
 328
 329                 face newF(3);
 330                 if (edgeInFaceOrder == (mesh_.faceOwner()[face1] == cellI))
 331                 {
 332                     // edge used in face order.
 333                     newF[0] = e[1];
 334                     newF[1] = e[0];
 335                     newF[2] = midPointI;
 336                 }
 337                 else
 338                 {
 339                     newF[0] = e[0];
 340                     newF[1] = e[1];
 341                     newF[2] = midPointI;
 342                 }
 343
 344                 // Now newF points away from cell1
 345                 meshMod.setAction
 346                 (
 347                     polyAddFace
 348                     (
 349                         newF,                       // face
 350                         cell1,                      // owner
 351                         cell0,                      // neighbour
 352                         -1,                         // master point
 353                         -1,                         // master edge
 354                         face0,                      // master face for addition
 355                         false,                      // flux flip
 356                         -1,                         // patch for face
 357                         -1,                         // zone for face
 358                         false                       // face zone flip
 359                     )
 360                 );
 361             }
 362         }
 363     }
 364
 365
 366     //
 367     // Update all existing faces for split owner or neighbour.
 368     //
 369
 370
 371     // Mark off affected face.
 372     boolList faceUpToDate(mesh_.nFaces(), true);
 373
 374     forAllConstIter(Map<point>, cellToMidPoint, iter)
 375     {
 376         label cellI = iter.key();
 377
 378         const cell& cFaces = mesh_.cells()[cellI];
 379
 380         forAll(cFaces, i)
 381         {
 382             label faceI = cFaces[i];
 383
 384             faceUpToDate[faceI] = false;
 385         }
 386     }
 387
 388     forAll(faceUpToDate, faceI)
 389     {
 390         if (!faceUpToDate[faceI])
 391         {
 392             const face& f = mesh_.faces()[faceI];
 393
 394             if (mesh_.isInternalFace(faceI))
 395             {
 396                 label newOwn = newOwner(faceI, cellToCells);
 397                 label newNbr = newNeighbour(faceI, cellToCells);
 398
 399                 if (newOwn < newNbr)
 400                 {
 401                     meshMod.setAction
 402                     (
 403                         polyModifyFace
 404                         (
 405                             f,
 406                             faceI,
 407                             newOwn,         // owner
 408                             newNbr,         // neighbour
 409                             false,          // flux flip
 410                             -1,             // patch for face
 411                             false,          // remove from zone
 412                             -1,             // zone for face
 413                             false           // face zone flip
 414                         )
 415                     );
 416                 }
 417                 else
 418                 {
 419                     meshMod.setAction
 420                     (
 421                         polyModifyFace
 422                         (
 423                             f.reverseFace(),
 424                             faceI,
 425                             newNbr,         // owner
 426                             newOwn,         // neighbour
 427                             false,          // flux flip
 428                             -1,             // patch for face
 429                             false,          // remove from zone
 430                             -1,             // zone for face
 431                             false           // face zone flip
 432                         )
 433                     );
 434                 }
 435
 436             }
 437             else
 438             {
 439                 label newOwn = newOwner(faceI, cellToCells);
 440
 441                 label patchID, zoneID, zoneFlip;
 442                 getFaceInfo(faceI, patchID, zoneID, zoneFlip);
 443
 444                 meshMod.setAction
 445                 (
 446                     polyModifyFace
 447                     (
 448                         mesh_.faces()[faceI],
 449                         faceI,
 450                         newOwn,         // owner
 451                         -1,             // neighbour
 452                         false,          // flux flip
 453                         patchID,        // patch for face
 454                         false,          // remove from zone
 455                         zoneID,         // zone for face
 456                         zoneFlip        // face zone flip
 457                     )
 458                 );
 459             }
 460
 461             faceUpToDate[faceI] = true;
 462         }
 463     }
 464 }
 465
 466
 467 void Foam::cellSplitter::updateMesh(const mapPolyMesh& morphMap)
 468 {
 469     // Create copy since we're deleting entries. Only if both cell and added
 470     // point get mapped do they get inserted.
 471     Map<label> newAddedPoints(addedPoints_.size());
 472
 473     forAllConstIter(Map<label>, addedPoints_, iter)
 474     {
 475         label oldCellI = iter.key();
 476
 477         label newCellI = morphMap.reverseCellMap()[oldCellI];
 478
 479         label oldPointI = iter();
 480
 481         label newPointI = morphMap.reversePointMap()[oldPointI];
 482
 483         if (newCellI >= 0 && newPointI >= 0)
 484         {
 485             newAddedPoints.insert(newCellI, newPointI);
 486         }
 487     }
 488
 489     // Copy
 490     addedPoints_.transfer(newAddedPoints);
 491 }
 492
 493
 494 // ************************************************************************* //