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