src/dynamicMesh/layerAdditionRemoval/addCellLayer.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2004-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
 232         //     << " own: " << mc[faceI]
 233         //     << " nei: " << addedCells[faceI]
 234         //     << endl;
 235     }
 236
 237     // Modify the faces from the master zone for the new neighbour
 238
 239     const faceList& faces = mesh.faces();
 240
 241     // Pout<< "mfFlip: " << mfFlip << endl;
 242
 243     forAll(mf, faceI)
 244     {
 245         const label curfaceID = mf[faceI];
 246
 247         // If the face is internal, modify its owner to be the newly
 248         // created cell.  No flip is necessary
 249         if (!mesh.isInternalFace(curfaceID))
 250         {
 251             ref.setAction
 252             (
 253                 polyModifyFace
 254                 (
 255                     faces[curfaceID],            // modified face
 256                     curfaceID,                   // label of face being modified
 257                     addedCells[faceI],           // owner
 258                     -1,                          // neighbour
 259                     false,                       // face flip
 260                     mesh.boundaryMesh().whichPatch(curfaceID),// patch for face
 261                     false,                       // remove from zone
 262                     faceZoneID_.index(),         // zone for face
 263                     mfFlip[faceI]                // face flip in zone
 264                 )
 265             );
 266
 267             // Pout<< "Modifying a boundary face. Face: " << curfaceID
 268             //     << " flip: " << mfFlip[faceI]
 269             //     << endl;
 270         }
 271         // If slave cell is owner, the face remains the same (but with
 272         // a new neighbour - the newly created cell).  Otherwise, the
 273         // face is flipped.
 274         else if (sc[faceI] == own[curfaceID])
 275         {
 276             // Orientation is good, just change neighbour
 277             ref.setAction
 278             (
 279                 polyModifyFace
 280                 (
 281                     faces[curfaceID],            // modified face
 282                     curfaceID,                   // label of face being modified
 283                     own[curfaceID],              // owner
 284                     addedCells[faceI],           // neighbour
 285                     false,                       // face flip
 286                     mesh.boundaryMesh().whichPatch(curfaceID),// patch for face
 287                     false,                       // remove from zone
 288                     faceZoneID_.index(),         // zone for face
 289                     mfFlip[faceI]                // face flip in zone
 290                 )
 291             );
 292
 293             // Pout<< "modify face, no flip " << curfaceID
 294             //     << " own: " << own[curfaceID]
 295             //     << " nei: " << addedCells[faceI]
 296             //     << endl;
 297         }
 298         else
 299         {
 300             // Change in orientation; flip face
 301             ref.setAction
 302             (
 303                 polyModifyFace
 304                 (
 305                     faces[curfaceID].reverseFace(), // modified face
 306                     curfaceID,                   // label of face being modified
 307                     nei[curfaceID],                 // owner
 308                     addedCells[faceI],              // neighbour
 309                     true,                           // face flip
 310                     mesh.boundaryMesh().whichPatch(curfaceID), // patch for face
 311                     false,                          // remove from zone
 312                     faceZoneID_.index(),            // zone for face
 313                     !mfFlip[faceI]                  // face flip in zone
 314                 )
 315             );
 316
 317             // Pout<< "modify face, with flip " << curfaceID
 318             //     << " own: " << own[curfaceID]
 319             //     << " nei: " << addedCells[faceI]
 320             //     << endl;
 321         }
 322     }
 323
 324     // Create new faces from edges
 325
 326     const edgeList& zoneLocalEdges = masterFaceLayer.edges();
 327
 328     const labelListList& edgeFaces = masterFaceLayer.edgeFaces();
 329
 330     label nInternalEdges = masterFaceLayer.nInternalEdges();
 331
 332     const labelList& meshEdges =
 333         mesh.faceZones()[faceZoneID_.index()].meshEdges();
 334
 335     // Do all internal edges
 336     for (label curEdgeID = 0; curEdgeID < nInternalEdges; curEdgeID++)
 337     {
 338         face newFace(4);
 339
 340         newFace[0] = mp[zoneLocalEdges[curEdgeID].start()];
 341         newFace[1] = mp[zoneLocalEdges[curEdgeID].end()];
 342         newFace[2] = addedPoints[zoneLocalEdges[curEdgeID].end()];
 343         newFace[3] = addedPoints[zoneLocalEdges[curEdgeID].start()];
 344
 345         ref.setAction
 346         (
 347             polyAddFace
 348             (
 349                 newFace,                               // face
 350                 addedCells[edgeFaces[curEdgeID][0]],   // owner
 351                 addedCells[edgeFaces[curEdgeID][1]],   // neighbour
 352                 -1,                                    // master point
 353                 meshEdges[curEdgeID],                  // master edge
 354                 -1,                                    // master face
 355                 false,                                 // flip flux
 356                 -1,                                    // patch for face
 357                 -1,                                    // zone for face
 358                 false                                  // face zone flip
 359             )
 360         );
 361
 362         // Pout<< "Add internal face off edge: " << newFace
 363         //     << " own: " << addedCells[edgeFaces[curEdgeID][0]]
 364         //     << " nei: " << addedCells[edgeFaces[curEdgeID][1]]
 365         //     << endl;
 366     }
 367
 368     // Prepare creation of faces from boundary edges.
 369     // Note:
 370     // A tricky part of the algorithm is finding the patch into which the
 371     // newly created face will be added.  For this purpose, take the edge
 372     // and grab all the faces coming from it.  From the set of faces
 373     // eliminate the internal faces and find the boundary face from the rest.
 374     //  If there is more than one boundary face (excluding the ones in
 375     // the master zone), the patch with the lower label is selected.
 376
 377     const labelListList& meshEdgeFaces = mesh.edgeFaces();
 378
 379     const faceZoneMesh& zoneMesh = mesh.faceZones();
 380
 381     // Do all boundary edges
 382
 383     for
 384     (
 385         label curEdgeID = nInternalEdges;
 386         curEdgeID < zoneLocalEdges.size();
 387         curEdgeID++
 388     )
 389     {
 390         face newFace(4);
 391         newFace[0] = mp[zoneLocalEdges[curEdgeID].start()];
 392         newFace[1] = mp[zoneLocalEdges[curEdgeID].end()];
 393         newFace[2] = addedPoints[zoneLocalEdges[curEdgeID].end()];
 394         newFace[3] = addedPoints[zoneLocalEdges[curEdgeID].start()];
 395
 396         // Determine the patch for insertion
 397         const labelList& curFaces = meshEdgeFaces[meshEdges[curEdgeID]];
 398
 399         label patchID = -1;
 400         label zoneID = -1;
 401
 402         forAll(curFaces, faceI)
 403         {
 404             const label cf = curFaces[faceI];
 405
 406             if (!mesh.isInternalFace(cf))
 407             {
 408                 // Face not internal. Check if it is in the zone
 409                 if (zoneMesh.whichZone(cf) != faceZoneID_.index())
 410                 {
 411                     // Found the face in a boundary patch which is not in zone
 412                     patchID = mesh.boundaryMesh().whichPatch(cf);
 413                     zoneID = mesh.faceZones().whichZone(cf);
 414
 415                     break;
 416                 }
 417             }
 418         }
 419
 420         if (patchID < 0)
 421         {
 422             FatalErrorIn
 423             (
 424                 "void Foam::layerAdditionRemoval::setRefinement("
 425                 "polyTopoChange& ref)"
 426             )   << "Cannot find patch for edge " << meshEdges[curEdgeID]
 427                 << ". Edge: " << mesh.edges()[meshEdges[curEdgeID]]
 428                 << abort(FatalError);
 429         }
 430
 431         ref.setAction
 432         (
 433             polyAddFace
 434             (
 435                 newFace,                               // face
 436                 addedCells[edgeFaces[curEdgeID][0]],   // owner
 437                 -1,                                    // neighbour
 438                 -1,                                    // master point
 439                 meshEdges[curEdgeID],                  // master edge
 440                 -1,                                    // master face
 441                 false,                                 // flip flux
 442                 patchID,                               // patch for face
 443                 zoneID,                                // zone
 444                 false                                  // zone face flip
 445             )
 446         );
 447
 448         // Pout<< "add boundary face: " << newFace
 449         //     << " into patch " << patchID
 450         //     << " own: " << addedCells[edgeFaces[curEdgeID][0]]
 451         //     << endl;
 452     }
 453
 454     // Modify the remaining faces of the master cells to reconnect to the new
 455     // layer of faces.
 456     // Algorithm: Go through all the cells of the master zone and make
 457     // a map of faces to avoid duplicates.  Do not insert the faces in
 458     // the master patch (as they have already been dealt with).  Make
 459     // a master layer point renumbering map, which for every point in
 460     // the master layer gives its new label. Loop through all faces in
 461     // the map and attempt to renumber them using the master layer
 462     // point renumbering map.  Once the face is renumbered, compare it
 463     // with the original face; if they are the same, the face has not
 464     // changed; if not, modify the face but keep all of its old
 465     // attributes (apart from the vertex numbers).
 466
 467     // Create the map of faces in the master cell layer
 468     labelHashSet masterCellFaceMap(primitiveMesh::facesPerCell_*mc.size());
 469
 470     const cellList& cells = mesh.cells();
 471
 472     forAll(mc, cellI)
 473     {
 474         const labelList& curFaces = cells[mc[cellI]];
 475
 476         forAll(curFaces, faceI)
 477         {
 478             // Check if the face belongs to the master zone; if not add it
 479             if (zoneMesh.whichZone(curFaces[faceI]) != faceZoneID_.index())
 480             {
 481                 masterCellFaceMap.insert(curFaces[faceI]);
 482             }
 483         }
 484     }
 485
 486     // Create the master layer point map
 487     Map<label> masterLayerPointMap(2*mp.size());
 488
 489     forAll(mp, pointI)
 490     {
 491         masterLayerPointMap.insert
 492         (
 493             mp[pointI],
 494             addedPoints[pointI]
 495         );
 496     }
 497
 498     // Grab the list of faces of the master layer
 499     const labelList masterCellFaces = masterCellFaceMap.toc();
 500
 501     forAll(masterCellFaces, faceI)
 502     {
 503         // Attempt to renumber the face using the masterLayerPointMap.
 504         // Missing point remain the same
 505
 506         const label curFaceID = masterCellFaces[faceI];
 507
 508         const face& oldFace = faces[curFaceID];
 509
 510         face newFace(oldFace.size());
 511
 512         bool changed = false;
 513
 514         forAll(oldFace, pointI)
 515         {
 516             if (masterLayerPointMap.found(oldFace[pointI]))
 517             {
 518                 changed = true;
 519
 520                 newFace[pointI] = masterLayerPointMap.find(oldFace[pointI])();
 521             }
 522             else
 523             {
 524                 newFace[pointI] = oldFace[pointI];
 525             }
 526         }
 527
 528         // If the face has changed, create a modification entry
 529         if (changed)
 530         {
 531             // Get face zone and its flip
 532             label modifiedFaceZone = mesh.faceZones().whichZone(curFaceID);
 533             bool modifiedFaceZoneFlip = false;
 534
 535             if (modifiedFaceZone >= 0)
 536             {
 537                 modifiedFaceZoneFlip =
 538                     mesh.faceZones()[modifiedFaceZone].flipMap()
 539                     [
 540                         mesh.faceZones()[modifiedFaceZone].whichFace(curFaceID)
 541                     ];
 542             }
 543
 544             if (mesh.isInternalFace(curFaceID))
 545             {
 546                 ref.setAction
 547                 (
 548                     polyModifyFace
 549                     (
 550                         newFace,                // modified face
 551                         curFaceID,              // label of face being modified
 552                         own[curFaceID],         // owner
 553                         nei[curFaceID],         // neighbour
 554                         false,                  // face flip
 555                         -1,                     // patch for face
 556                         false,                  // remove from zone
 557                         modifiedFaceZone,       // zone for face
 558                         modifiedFaceZoneFlip    // face flip in zone
 559                     )
 560                 );
 561
 562                 // Pout<< "modifying stick-out face. Internal Old face: "
 563                 //     << oldFace
 564                 //     << " new face: " << newFace
 565                 //     << " own: " << own[curFaceID]
 566                 //     << " nei: " << nei[curFaceID]
 567                 //     << endl;
 568             }
 569             else
 570             {
 571                 ref.setAction
 572                 (
 573                     polyModifyFace
 574                     (
 575                         newFace,                // modified face
 576                         curFaceID,              // label of face being modified
 577                         own[curFaceID],         // owner
 578                         -1,                     // neighbour
 579                         false,                  // face flip
 580                         mesh.boundaryMesh().whichPatch(curFaceID),
 581                                                 // patch for face
 582                         false,                  // remove from zone
 583                         modifiedFaceZone,       // zone for face
 584                         modifiedFaceZoneFlip    // face flip in zone
 585                     )
 586                 );
 587
 588                 // Pout<< "modifying stick-out face. Boundary Old face: "
 589                 //     << oldFace
 590                 //     << " new face: " << newFace
 591                 //     << " own: " << own[curFaceID]
 592                 //     << " patch: "
 593                 //     << mesh.boundaryMesh().whichPatch(curFaceID)
 594                 //     << endl;
 595             }
 596         }
 597     }
 598
 599     if (debug)
 600     {
 601         Pout<< "void layerAdditionRemoval::addCellLayer("
 602             << "polyTopoChange& ref) const "
 603             << " for object " << name() << " : "
 604             << "Finished adding cell layer" << endl;
 605     }
 606 }
 607
 608
 609 // ************************************************************************* //