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