src/dynamicMesh/layerAdditionRemoval/addCellLayer.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2010 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
  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 "layerAdditionRemoval.H"
  27 #include "polyMesh.H"
  28 #include "primitiveMesh.H"
  29 #include "polyTopoChange.H"
  30 #include "polyTopoChanger.H"
  31 #include "polyAddPoint.H"
  32 #include "polyAddCell.H"
  33 #include "polyAddFace.H"
  34 #include "polyModifyFace.H"
  35
  36 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  37
  38 //- Calculate the offset to the next layer
  39 Foam::tmp<Foam::vectorField> Foam::layerAdditionRemoval::extrusionDir() const
  40 {
  41     const polyMesh& mesh = topoChanger().mesh();
  42     const primitiveFacePatch& masterFaceLayer =
  43         mesh.faceZones()[faceZoneID_.index()]();
  44
  45     const pointField& points = mesh.points();
  46     const labelList& mp = masterFaceLayer.meshPoints();
  47
  48     tmp<vectorField> textrusionDir(new vectorField(mp.size()));
  49     vectorField& extrusionDir = textrusionDir();
  50
  51     if (setLayerPairing())
  52     {
  53         if (debug)
  54         {
  55             Pout<< "void layerAdditionRemoval::extrusionDir() const "
  56                 << " for object " << name() << " : "
  57                 << "Using edges for point insertion" << endl;
  58         }
  59
  60         // Detected a valid layer.  Grab the point and face collapse mapping
  61         const labelList& ptc = pointsPairing();
  62
  63         forAll (extrusionDir, mpI)
  64         {
  65             extrusionDir[mpI] = points[ptc[mpI]] - points[mp[mpI]];
  66         }
  67     }
  68     else
  69     {
  70         if (debug)
  71         {
  72             Pout<< "void layerAdditionRemoval::extrusionDir() const "
  73                 << " for object " << name() << " : "
  74                 << "A valid layer could not be found in front of "
  75                 << "the addition face layer.  Using face-based "
  76                 << "point normals for point addition"
  77                 << endl;
  78         }
  79
  80         extrusionDir = minLayerThickness_*masterFaceLayer.pointNormals();
  81     }
  82     return textrusionDir;
  83 }
  84
  85
  86 void Foam::layerAdditionRemoval::addCellLayer
  87 (
  88     polyTopoChange& ref
  89 ) const
  90 {
  91     // Insert the layer addition instructions into the topological change
  92
  93     // Algorithm:
  94     // 1.  For every point in the master zone add a new point
  95     // 2.  For every face in the master zone add a cell
  96     // 3.  Go through all the edges of the master zone.  For all internal edges,
  97     //     add a face with the correct orientation and make it internal.
  98     //     For all boundary edges, find the patch it belongs to and add the face
  99     //     to this patch
 100     // 4.  Create a set of new faces on the patch image and assign them to be
 101     //     between the old master cells and the newly created cells
 102     // 5.  Modify all the faces in the patch such that they are located
 103     //     between the old slave cells and newly created cells
 104
 105     if (debug)
 106     {
 107         Pout<< "void layerAdditionRemoval::addCellLayer("
 108             << "polyTopoChange& ref) const for object " << name() << " : "
 109             << "Adding cell layer" << endl;
 110     }
 111
 112     // Create the points
 113
 114     const polyMesh& mesh = topoChanger().mesh();
 115     const primitiveFacePatch& masterFaceLayer =
 116         mesh.faceZones()[faceZoneID_.index()]();
 117
 118     const pointField& points = mesh.points();
 119     const labelList& mp = masterFaceLayer.meshPoints();
 120
 121     // Get the extrusion direction for the added points
 122
 123     tmp<vectorField> tpointOffsets = extrusionDir();
 124
 125     // Add the new points
 126     labelList addedPoints(mp.size());
 127
 128     forAll (mp, pointI)
 129     {
 130         // Add the point nominal thickness away from the master zone point
 131         // and grab the label
 132         addedPoints[pointI] =
 133             ref.setAction
 134             (
 135                 polyAddPoint
 136                 (
 137                     points[mp[pointI]]                  // point
 138                   + addDelta_*tpointOffsets()[pointI],
 139                     mp[pointI],                         // master point
 140                     -1,                                 // zone for point
 141                     true                                // supports a cell
 142                 )
 143             );
 144     }
 145
 146 // Pout << "mp: " << mp << " addedPoints: " << addedPoints << endl;
 147     // Create the cells
 148
 149     const labelList& mc =
 150         mesh.faceZones()[faceZoneID_.index()].masterCells();
 151     const labelList& sc =
 152         mesh.faceZones()[faceZoneID_.index()].slaveCells();
 153 // Pout << "mc: " << mc << " sc: " << sc << endl;
 154     const labelList& mf = mesh.faceZones()[faceZoneID_.index()];
 155     const boolList& mfFlip = mesh.faceZones()[faceZoneID_.index()].flipMap();
 156
 157     const labelList& own = mesh.faceOwner();
 158     const labelList& nei = mesh.faceNeighbour();
 159
 160     labelList addedCells(mf.size());
 161
 162     forAll (mf, faceI)
 163     {
 164         addedCells[faceI] =
 165             ref.setAction
 166             (
 167                 polyAddCell
 168                 (
 169                     -1,           // master point
 170                     -1,           // master edge
 171                     mf[faceI],    // master face
 172                     -1,           // master cell
 173                     -1            // zone for cell
 174                 )
 175             );
 176     }
 177
 178     // Create the new faces
 179
 180     // Grab the local faces from the master zone
 181     const faceList& zoneFaces = masterFaceLayer.localFaces();
 182
 183     // Create the new faces by copying the master zone.  All the added
 184     // faces need to point into the newly created cells, which means
 185     // all the faces from the master layer are flipped.  The flux flip
 186     // is determined from the original master layer face and the face
 187     // owner: if the master cell is equal to the face owner the flux
 188     // remains the same; otherwise it is flipped
 189
 190     forAll (zoneFaces, faceI)
 191     {
 192         const face oldFace = zoneFaces[faceI].reverseFace();
 193
 194         face newFace(oldFace.size());
 195
 196         forAll (oldFace, pointI)
 197         {
 198             newFace[pointI] = addedPoints[oldFace[pointI]];
 199         }
 200
 201         bool flipFaceFlux = false;
 202
 203         // Flip the face as necessary
 204         if
 205         (
 206             mc[faceI] == nei[mf[faceI]]
 207          || !mesh.isInternalFace(mf[faceI])
 208         )
 209         {
 210             flipFaceFlux = true;
 211         }
 212
 213         // Add the face
 214         ref.setAction
 215         (
 216             polyAddFace
 217             (
 218                 newFace,               // face
 219                 mc[faceI],             // owner
 220                 addedCells[faceI],     // neighbour
 221                 -1,                    // master point
 222                 -1,                    // master edge
 223                 mf[faceI],             // master face for addition
 224                 flipFaceFlux,          // flux flip
 225                 -1,                    // patch for face
 226                 -1,                    // zone for face
 227                 false                  // face zone flip
 228             )
 229         );
 230
 231 // Pout << "adding face: " << newFace << " own: " << mc[faceI] << " nei: " << addedCells[faceI] << endl;
 232     }
 233
 234     // Modify the faces from the master zone for the new neighbour
 235
 236     const faceList& faces = mesh.faces();
 237 // Pout << "mfFlip: " << mfFlip << endl;
 238     forAll (mf, faceI)
 239     {
 240         const label curfaceID = mf[faceI];
 241
 242         // If the face is internal, modify its owner to be the newly
 243         // created cell.  No flip is necessary
 244         if (!mesh.isInternalFace(curfaceID))
 245         {
 246             ref.setAction
 247             (
 248                 polyModifyFace
 249                 (
 250                     faces[curfaceID],            // modified face
 251                     curfaceID,                   // label of face being modified
 252                     addedCells[faceI],           // owner
 253                     -1,                          // neighbour
 254                     false,                       // face flip
 255                     mesh.boundaryMesh().whichPatch(curfaceID),// patch for face
 256                     false,                       // remove from zone
 257                     faceZoneID_.index(),         // zone for face
 258                     mfFlip[faceI]                // face flip in zone
 259                 )
 260             );
 261 // Pout << "Modifying a boundary face. Face: " << curfaceID << " flip: " << mfFlip[faceI] << endl;
 262         }
 263         // If slave cell is owner, the face remains the same (but with
 264         // a new neighbour - the newly created cell).  Otherwise, the
 265         // face is flipped.
 266         else if (sc[faceI] == own[curfaceID])
 267         {
 268             // Orientation is good, just change neighbour
 269             ref.setAction
 270             (
 271                 polyModifyFace
 272                 (
 273                     faces[curfaceID],            // modified face
 274                     curfaceID,                   // label of face being modified
 275                     own[curfaceID],              // owner
 276                     addedCells[faceI],           // neighbour
 277                     false,                       // face flip
 278                     mesh.boundaryMesh().whichPatch(curfaceID),// patch for face
 279                     false,                       // remove from zone
 280                     faceZoneID_.index(),         // zone for face
 281                     mfFlip[faceI]                // face flip in zone
 282                 )
 283             );
 284
 285 // Pout << "modify face, no flip " << curfaceID << " own: " << own[curfaceID] << " nei: " << addedCells[faceI] << endl;
 286         }
 287         else
 288         {
 289             // Change in orientation; flip face
 290             ref.setAction
 291             (
 292                 polyModifyFace
 293                 (
 294                     faces[curfaceID].reverseFace(), // modified face
 295                     curfaceID,                   // label of face being modified
 296                     nei[curfaceID],                 // owner
 297                     addedCells[faceI],              // neighbour
 298                     true,                           // face flip
 299                     mesh.boundaryMesh().whichPatch(curfaceID), // patch for face
 300                     false,                          // remove from zone
 301                     faceZoneID_.index(),            // zone for face
 302                     !mfFlip[faceI]                  // face flip in zone
 303                 )
 304             );
 305 // Pout << "modify face, with flip " << curfaceID << " own: " << own[curfaceID] << " nei: " << addedCells[faceI] << endl;
 306         }
 307     }
 308
 309     // Create new faces from edges
 310
 311     const edgeList& zoneLocalEdges = masterFaceLayer.edges();
 312
 313     const labelListList& edgeFaces = masterFaceLayer.edgeFaces();
 314
 315     label nInternalEdges = masterFaceLayer.nInternalEdges();
 316
 317     const labelList& meshEdges =
 318         mesh.faceZones()[faceZoneID_.index()].meshEdges();
 319
 320     // Do all internal edges
 321     for (label curEdgeID = 0; curEdgeID < nInternalEdges; curEdgeID++)
 322     {
 323         face newFace(4);
 324
 325         newFace[0] = mp[zoneLocalEdges[curEdgeID].start()];
 326         newFace[1] = mp[zoneLocalEdges[curEdgeID].end()];
 327         newFace[2] = addedPoints[zoneLocalEdges[curEdgeID].end()];
 328         newFace[3] = addedPoints[zoneLocalEdges[curEdgeID].start()];
 329
 330         ref.setAction
 331         (
 332             polyAddFace
 333             (
 334                 newFace,                               // face
 335                 addedCells[edgeFaces[curEdgeID][0]],   // owner
 336                 addedCells[edgeFaces[curEdgeID][1]],   // neighbour
 337                 -1,                                    // master point
 338                 meshEdges[curEdgeID],                  // master edge
 339                 -1,                                    // master face
 340                 false,                                 // flip flux
 341                 -1,                                    // patch for face
 342                 -1,                                    // zone for face
 343                 false                                  // face zone flip
 344             )
 345         );
 346
 347 // Pout << "Add internal face off edge: " << newFace << " own: " << addedCells[edgeFaces[curEdgeID][0]] << " nei: " << addedCells[edgeFaces[curEdgeID][1]] << endl;
 348     }
 349
 350     // Prepare creation of faces from boundary edges.
 351     // Note:
 352     // A tricky part of the algorithm is finding the patch into which the
 353     // newly created face will be added.  For this purpose, take the edge
 354     // and grab all the faces coming from it.  From the set of faces
 355     // eliminate the internal faces and find the boundary face from the rest.
 356     //  If there is more than one boundary face (excluding the ones in
 357     // the master zone), the patch with the lower label is selected.
 358
 359     const labelListList& meshEdgeFaces = mesh.edgeFaces();
 360
 361     const faceZoneMesh& zoneMesh = mesh.faceZones();
 362
 363     // Do all boundary edges
 364
 365     for
 366     (
 367         label curEdgeID = nInternalEdges;
 368         curEdgeID < zoneLocalEdges.size();
 369         curEdgeID++
 370     )
 371     {
 372         face newFace(4);
 373         newFace[0] = mp[zoneLocalEdges[curEdgeID].start()];
 374         newFace[1] = mp[zoneLocalEdges[curEdgeID].end()];
 375         newFace[2] = addedPoints[zoneLocalEdges[curEdgeID].end()];
 376         newFace[3] = addedPoints[zoneLocalEdges[curEdgeID].start()];
 377
 378         // Determine the patch for insertion
 379         const labelList& curFaces = meshEdgeFaces[meshEdges[curEdgeID]];
 380
 381         label patchID = -1;
 382         label zoneID = -1;
 383
 384         forAll (curFaces, faceI)
 385         {
 386             const label cf = curFaces[faceI];
 387
 388             if (!mesh.isInternalFace(cf))
 389             {
 390                 // Face not internal. Check if it is in the zone
 391                 if (zoneMesh.whichZone(cf) != faceZoneID_.index())
 392                 {
 393                     // Found the face in a boundary patch which is not in zone
 394                     patchID = mesh.boundaryMesh().whichPatch(cf);
 395                     zoneID = mesh.faceZones().whichZone(cf);
 396
 397                     break;
 398                 }
 399             }
 400         }
 401
 402         if (patchID < 0)
 403         {
 404             FatalErrorIn
 405             (
 406                 "void Foam::layerAdditionRemoval::setRefinement("
 407                 "polyTopoChange& ref)"
 408             )   << "Cannot find patch for edge " << meshEdges[curEdgeID]
 409                 << ". Edge: " << mesh.edges()[meshEdges[curEdgeID]]
 410                 << abort(FatalError);
 411         }
 412
 413         ref.setAction
 414         (
 415             polyAddFace
 416             (
 417                 newFace,                               // face
 418                 addedCells[edgeFaces[curEdgeID][0]],   // owner
 419                 -1,                                    // neighbour
 420                 -1,                                    // master point
 421                 meshEdges[curEdgeID],                  // master edge
 422                 -1,                                    // master face
 423                 false,                                 // flip flux
 424                 patchID,                               // patch for face
 425                 zoneID,                                // zone
 426                 false                                  // zone face flip
 427             )
 428         );
 429 // Pout << "add boundary face: " << newFace << " into patch " << patchID << " own: " << addedCells[edgeFaces[curEdgeID][0]] << endl;
 430     }
 431
 432     // Modify the remaining faces of the master cells to reconnect to the new
 433     // layer of faces.
 434     // Algorithm: Go through all the cells of the master zone and make
 435     // a map of faces to avoid duplicates.  Do not insert the faces in
 436     // the master patch (as they have already been dealt with).  Make
 437     // a master layer point renumbering map, which for every point in
 438     // the master layer gives its new label. Loop through all faces in
 439     // the map and attempt to renumber them using the master layer
 440     // point renumbering map.  Once the face is renumbered, compare it
 441     // with the original face; if they are the same, the face has not
 442     // changed; if not, modify the face but keep all of its old
 443     // attributes (apart from the vertex numbers).
 444
 445     // Create the map of faces in the master cell layer
 446     labelHashSet masterCellFaceMap(primitiveMesh::facesPerCell_*mc.size());
 447
 448     const cellList& cells = mesh.cells();
 449
 450     forAll (mc, cellI)
 451     {
 452         const labelList& curFaces = cells[mc[cellI]];
 453
 454         forAll (curFaces, faceI)
 455         {
 456             // Check if the face belongs to the master zone; if not add it
 457             if (zoneMesh.whichZone(curFaces[faceI]) != faceZoneID_.index())
 458             {
 459                 masterCellFaceMap.insert(curFaces[faceI]);
 460             }
 461         }
 462     }
 463
 464     // Create the master layer point map
 465     Map<label> masterLayerPointMap(2*mp.size());
 466
 467     forAll (mp, pointI)
 468     {
 469         masterLayerPointMap.insert
 470         (
 471             mp[pointI],
 472             addedPoints[pointI]
 473         );
 474     }
 475
 476     // Grab the list of faces of the master layer
 477     const labelList masterCellFaces = masterCellFaceMap.toc();
 478
 479     forAll (masterCellFaces, faceI)
 480     {
 481         // Attempt to renumber the face using the masterLayerPointMap.
 482         // Missing point remain the same
 483
 484         const label curFaceID = masterCellFaces[faceI];
 485
 486         const face& oldFace = faces[curFaceID];
 487
 488         face newFace(oldFace.size());
 489
 490         bool changed = false;
 491
 492         forAll (oldFace, pointI)
 493         {
 494             if (masterLayerPointMap.found(oldFace[pointI]))
 495             {
 496                 changed = true;
 497
 498                 newFace[pointI] = masterLayerPointMap.find(oldFace[pointI])();
 499             }
 500             else
 501             {
 502                 newFace[pointI] = oldFace[pointI];
 503             }
 504         }
 505
 506         // If the face has changed, create a modification entry
 507         if (changed)
 508         {
 509             // Get face zone and its flip
 510             label modifiedFaceZone = mesh.faceZones().whichZone(curFaceID);
 511             bool modifiedFaceZoneFlip = false;
 512
 513             if (modifiedFaceZone >= 0)
 514             {
 515                 modifiedFaceZoneFlip =
 516                     mesh.faceZones()[modifiedFaceZone].flipMap()
 517                     [
 518                         mesh.faceZones()[modifiedFaceZone].whichFace(curFaceID)
 519                     ];
 520             }
 521
 522             if (mesh.isInternalFace(curFaceID))
 523             {
 524                 ref.setAction
 525                 (
 526                     polyModifyFace
 527                     (
 528                         newFace,                // modified face
 529                         curFaceID,              // label of face being modified
 530                         own[curFaceID],         // owner
 531                         nei[curFaceID],         // neighbour
 532                         false,                  // face flip
 533                         -1,                     // patch for face
 534                         false,                  // remove from zone
 535                         modifiedFaceZone,       // zone for face
 536                         modifiedFaceZoneFlip    // face flip in zone
 537                     )
 538                 );
 539 // Pout << "modifying stick-out face. Internal Old face: " << oldFace << " new face: " << newFace << " own: " << own[curFaceID] << " nei: " << nei[curFaceID] << endl;
 540             }
 541             else
 542             {
 543                 ref.setAction
 544                 (
 545                     polyModifyFace
 546                     (
 547                         newFace,                // modified face
 548                         curFaceID,              // label of face being modified
 549                         own[curFaceID],         // owner
 550                         -1,                     // neighbour
 551                         false,                  // face flip
 552                         mesh.boundaryMesh().whichPatch(curFaceID), // patch for face
 553                         false,                  // remove from zone
 554                         modifiedFaceZone,       // zone for face
 555                         modifiedFaceZoneFlip    // face flip in zone
 556                     )
 557                 );
 558 // Pout << "modifying stick-out face. Boundary Old face: " << oldFace << " new face: " << newFace << " own: " << own[curFaceID] << " patch: " << mesh.boundaryMesh().whichPatch(curFaceID) << endl;
 559             }
 560         }
 561     }
 562
 563     if (debug)
 564     {
 565         Pout<< "void layerAdditionRemoval::addCellLayer("
 566             << "polyTopoChange& ref) const "
 567             << " for object " << name() << " : "
 568             << "Finished adding cell layer" << endl;
 569     }
 570 }
 571
 572
 573 // ************************************************************************* //