src/dynamicMesh/polyMeshModifiers/layerAdditionRemoval/addCellLayer.C

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