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