src/dynamicMesh/attachDetach/detachInterface.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 "attachDetach.H"
  27 #include "polyMesh.H"
  28 #include "primitiveMesh.H"
  29 #include "polyTopoChange.H"
  30 #include "polyTopoChanger.H"
  31 #include "polyAddPoint.H"
  32 #include "polyModifyFace.H"
  33 #include "polyAddFace.H"
  34
  35 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  36
  37 void Foam::attachDetach::detachInterface
  38 (
  39     polyTopoChange& ref
  40 ) const
  41 {
  42     // Algorithm:
  43     // 1. Create new points for points of the master face zone
  44     // 2. Modify all faces of the master zone, by putting them into the master
  45     //    patch (look for orientation) and their renumbered mirror images
  46     //    into the slave patch
  47     // 3. Create a point renumbering list, giving a new point index for original
  48     //    points in the face patch
  49     // 4. Grab all faces in cells on the master side and renumber them
  50     //    using the point renumbering list.  Exclude the ones that belong to
  51     //    the master face zone
  52     //
  53     // Note on point creation:
  54     // In order to take into account the issues related to partial
  55     // blocking in an attach/detach mesh modifier, special treatment
  56     // is required for the duplication of points on the edge of the
  57     // face zone.  Points which are shared only by internal edges need
  58     // not to be duplicated, as this would propagate the discontinuity
  59     // in the mesh beyond the face zone.  Therefore, before creating
  60     // the new points, check the external edge loop.  For each edge
  61     // check if the edge is internal (i.e. does not belong to a
  62     // patch); if so, exclude both of its points from duplication.
  63
  64     if (debug)
  65     {
  66         Pout<< "void attachDetach::detachInterface("
  67             << "polyTopoChange& ref) const "
  68             << " for object " << name() << " : "
  69             << "Detaching interface" << endl;
  70     }
  71
  72     const polyMesh& mesh = topoChanger().mesh();
  73     const faceZoneMesh& zoneMesh = mesh.faceZones();
  74
  75     // Check that zone is in increasing order (needed since adding faces
  76     // in same order - otherwise polyTopoChange face ordering will mess up
  77     // correspondence)
  78     if (debug)
  79     {
  80         const labelList& faceLabels = zoneMesh[faceZoneID_.index()];
  81         if (faceLabels.size() > 0)
  82         {
  83             for (label i = 1; i < faceLabels.size(); i++)
  84             {
  85                 if (faceLabels[i] <= faceLabels[i-1])
  86                 {
  87                     FatalErrorIn
  88                     (
  89                         "attachDetach::detachInterface"
  90                         "(polyTopoChange&) const"
  91                     )   << "faceZone " << zoneMesh[faceZoneID_.index()].name()
  92                         << " does not have mesh face labels in"
  93                         << " increasing order." << endl
  94                         << "Face label " << faceLabels[i]
  95                         << " at position " << i
  96                         << " is smaller than the previous value "
  97                         << faceLabels[i-1]
  98                         << exit(FatalError);
  99                 }
 100             }
 101         }
 102     }
 103
 104
 105
 106     const primitiveFacePatch& masterFaceLayer = zoneMesh[faceZoneID_.index()]();
 107     const pointField& points = mesh.points();
 108     const labelListList& meshEdgeFaces = mesh.edgeFaces();
 109
 110     const labelList& mp = masterFaceLayer.meshPoints();
 111     const edgeList& zoneLocalEdges = masterFaceLayer.edges();
 112
 113     const labelList& meshEdges = zoneMesh[faceZoneID_.index()].meshEdges();
 114
 115     // Create the points
 116
 117     labelList addedPoints(mp.size(), -1);
 118
 119     // Go through boundary edges of the master patch.  If all the faces from
 120     // this patch are internal, mark the points in the addedPoints lookup
 121     // with their original labels to stop duplication
 122     label nIntEdges = masterFaceLayer.nInternalEdges();
 123
 124     for (label curEdgeID = nIntEdges; curEdgeID < meshEdges.size(); curEdgeID++)
 125     {
 126         const labelList& curFaces = meshEdgeFaces[meshEdges[curEdgeID]];
 127
 128         bool edgeIsInternal = true;
 129
 130         forAll (curFaces, faceI)
 131         {
 132             if (!mesh.isInternalFace(curFaces[faceI]))
 133             {
 134                 // The edge belongs to a boundary face
 135                 edgeIsInternal = false;
 136                 break;
 137             }
 138         }
 139
 140         if (edgeIsInternal)
 141         {
 142             const edge& e = zoneLocalEdges[curEdgeID];
 143
 144             // Reset the point creation
 145             addedPoints[e.start()] = mp[e.start()];
 146             addedPoints[e.end()] = mp[e.end()];
 147         }
 148     }
 149 // Pout << "addedPoints before point creation: " << addedPoints << endl;
 150
 151     // Create new points for face zone
 152     forAll (addedPoints, pointI)
 153     {
 154         if (addedPoints[pointI] < 0)
 155         {
 156             addedPoints[pointI] =
 157                 ref.setAction
 158                 (
 159                     polyAddPoint
 160                     (
 161                         points[mp[pointI]],        // point
 162                         mp[pointI],                // master point
 163                         -1,                        // zone ID
 164                         true                       // supports a cell
 165                     )
 166                 );
 167             //Pout<< "Adding point " << addedPoints[pointI]
 168             //    << " coord1:" << points[mp[pointI]]
 169             //    << " coord2:" << masterFaceLayer.localPoints()[pointI]
 170             //    << " for original point " << mp[pointI] << endl;
 171         }
 172     }
 173
 174     // Modify faces in the master zone and duplicate for the slave zone
 175
 176     const labelList& mf = zoneMesh[faceZoneID_.index()];
 177     const boolList& mfFlip = zoneMesh[faceZoneID_.index()].flipMap();
 178     const faceList& zoneFaces = masterFaceLayer.localFaces();
 179
 180     const faceList& faces = mesh.faces();
 181     const labelList& own = mesh.faceOwner();
 182     const labelList& nei = mesh.faceNeighbour();
 183
 184     forAll (mf, faceI)
 185     {
 186         const label curFaceID = mf[faceI];
 187
 188         // Build the face for the slave patch by renumbering
 189         const face oldFace = zoneFaces[faceI].reverseFace();
 190
 191         face newFace(oldFace.size());
 192
 193         forAll (oldFace, pointI)
 194         {
 195             newFace[pointI] = addedPoints[oldFace[pointI]];
 196         }
 197
 198         if (mfFlip[faceI])
 199         {
 200             // Face needs to be flipped for the master patch
 201             ref.setAction
 202             (
 203                 polyModifyFace
 204                 (
 205                     faces[curFaceID].reverseFace(), // modified face
 206                     curFaceID,                   // label of face being modified
 207                     nei[curFaceID],                 // owner
 208                     -1,                             // neighbour
 209                     true,                           // face flip
 210                     masterPatchID_.index(),         // patch for face
 211                     false,                          // remove from zone
 212                     faceZoneID_.index(),            // zone for face
 213                     !mfFlip[faceI]                  // face flip in zone
 214                 )
 215             );
 216
 217             // Add renumbered face into the slave patch
 218             //label addedFaceI =
 219             ref.setAction
 220             (
 221                 polyAddFace
 222                 (
 223                     newFace,                        // face
 224                     own[curFaceID],                 // owner
 225                     -1,                             // neighbour
 226                     -1,                             // master point
 227                     -1,                             // master edge
 228                     curFaceID,                      // master face
 229                     false,                          // flip flux
 230                     slavePatchID_.index(),          // patch to add the face to
 231                     -1,                             // zone for face
 232                     false                           // zone flip
 233                 )
 234             );
 235             //{
 236             //    pointField newPts(ref.points());
 237             //Pout<< "Flip.  Modifying face: " << ref.faces()[curFaceID]
 238             //    << " fc:" <<  ref.faces()[curFaceID].centre(newPts)
 239             //    << " next to cell: " << nei[curFaceID]
 240             //    << " and adding face: " << newFace
 241             //    << " fc:" << ref.faces()[addedFaceI].centre(newPts)
 242             //    << " next to cell: " << own[curFaceID] << endl;
 243             //}
 244         }
 245         else
 246         {
 247             // No flip
 248             ref.setAction
 249             (
 250                 polyModifyFace
 251                 (
 252                     faces[curFaceID],         // modified face
 253                     curFaceID,                // label of face being modified
 254                     own[curFaceID],           // owner
 255                     -1,                       // neighbour
 256                     false,                    // face flip
 257                     masterPatchID_.index(),   // patch for face
 258                     false,                    // remove from zone
 259                     faceZoneID_.index(),      // zone for face
 260                     mfFlip[faceI]             // face flip in zone
 261                 )
 262             );
 263
 264             // Add renumbered face into the slave patch
 265             //label addedFaceI =
 266             ref.setAction
 267             (
 268                 polyAddFace
 269                 (
 270                     newFace,                        // face
 271                     nei[curFaceID],                 // owner
 272                     -1,                             // neighbour
 273                     -1,                             // master point
 274                     -1,                             // master edge
 275                     curFaceID,                      // master face
 276                     true,                           // flip flux
 277                     slavePatchID_.index(),          // patch to add the face to
 278                     -1,                             // zone for face
 279                     false                           // face flip in zone
 280                 )
 281             );
 282             //{
 283             //    pointField newPts(ref.points());
 284             //Pout<< "No flip.  Modifying face: " << ref.faces()[curFaceID]
 285             //    << " fc:" <<  ref.faces()[curFaceID].centre(newPts)
 286             //    << " next to cell: " << own[curFaceID]
 287             //    << " and adding face: " << newFace
 288             //    << " fc:" << ref.faces()[addedFaceI].centre(newPts)
 289             //    << " next to cell: " << nei[curFaceID] << endl;
 290             //}
 291         }
 292     }
 293
 294     // Modify the remaining faces of the master cells to reconnect to the new
 295     // layer of faces.
 296
 297     // Algorithm: Go through all the cells of the master zone and make
 298     // a map of faces to avoid duplicates.  Do not insert the faces in
 299     // the master patch (as they have already been dealt with).  Make
 300     // a master layer point renumbering map, which for every point in
 301     // the master layer gives its new label. Loop through all faces in
 302     // the map and attempt to renumber them using the master layer
 303     // point renumbering map.  Once the face is renumbered, compare it
 304     // with the original face; if they are the same, the face has not
 305     // changed; if not, modify the face but keep all of its old
 306     // attributes (apart from the vertex numbers).
 307
 308     // Create the map of faces in the master cell layer
 309     const labelList& mc =
 310         mesh.faceZones()[faceZoneID_.index()].masterCells();
 311
 312     labelHashSet masterCellFaceMap(6*mc.size());
 313
 314     const cellList& cells = mesh.cells();
 315
 316     forAll (mc, cellI)
 317     {
 318         const labelList& curFaces = cells[mc[cellI]];
 319
 320         forAll (curFaces, faceI)
 321         {
 322             // Check if the face belongs to the master patch; if not add it
 323             if (zoneMesh.whichZone(curFaces[faceI]) != faceZoneID_.index())
 324             {
 325                 masterCellFaceMap.insert(curFaces[faceI]);
 326             }
 327         }
 328     }
 329
 330     // Extend the map to include first neighbours of the master cells to
 331     // deal with multiple corners.
 332     { // Protection and memory management
 333         // Make a map of master cells for quick reject
 334         labelHashSet mcMap(2*mc.size());
 335
 336         forAll (mc, mcI)
 337         {
 338             mcMap.insert(mc[mcI]);
 339         }
 340
 341         // Go through all the faces in the masterCellFaceMap.  If the
 342         // cells around them are not already used, add all of their
 343         // faces to the map
 344         const labelList mcf = masterCellFaceMap.toc();
 345
 346         forAll (mcf, mcfI)
 347         {
 348             // Do the owner side
 349             const label ownCell = own[mcf[mcfI]];
 350
 351             if (!mcMap.found(ownCell))
 352             {
 353                 // Cell not found. Add its faces to the map
 354                 const cell& curFaces = cells[ownCell];
 355
 356                 forAll (curFaces, faceI)
 357                 {
 358                     masterCellFaceMap.insert(curFaces[faceI]);
 359                 }
 360             }
 361
 362             // Do the neighbour side if face is internal
 363             if (mesh.isInternalFace(mcf[mcfI]))
 364             {
 365                 const label neiCell = nei[mcf[mcfI]];
 366
 367                 if (!mcMap.found(neiCell))
 368                 {
 369                     // Cell not found. Add its faces to the map
 370                     const cell& curFaces = cells[neiCell];
 371
 372                     forAll (curFaces, faceI)
 373                     {
 374                         masterCellFaceMap.insert(curFaces[faceI]);
 375                     }
 376                 }
 377             }
 378         }
 379     }
 380
 381     // Create the master layer point map
 382     Map<label> masterLayerPointMap(2*mp.size());
 383
 384     forAll (mp, pointI)
 385     {
 386         masterLayerPointMap.insert
 387         (
 388             mp[pointI],
 389             addedPoints[pointI]
 390         );
 391     }
 392
 393     // Grab the list of faces of the master layer
 394     const labelList masterCellFaces = masterCellFaceMap.toc();
 395
 396     forAll (masterCellFaces, faceI)
 397     {
 398         // Attempt to renumber the face using the masterLayerPointMap.
 399         // Missing point remain the same
 400
 401         const label curFaceID = masterCellFaces[faceI];
 402
 403         const face& oldFace = faces[curFaceID];
 404
 405         face newFace(oldFace.size());
 406
 407         bool changed = false;
 408
 409         forAll (oldFace, pointI)
 410         {
 411             if (masterLayerPointMap.found(oldFace[pointI]))
 412             {
 413                 changed = true;
 414
 415                 newFace[pointI] = masterLayerPointMap.find(oldFace[pointI])();
 416             }
 417             else
 418             {
 419                 newFace[pointI] = oldFace[pointI];
 420             }
 421         }
 422
 423         // If the face has changed, create a modification entry
 424         if (changed)
 425         {
 426             if (mesh.isInternalFace(curFaceID))
 427             {
 428                 ref.setAction
 429                 (
 430                     polyModifyFace
 431                     (
 432                         newFace,                    // face
 433                         curFaceID,                  // master face
 434                         own[curFaceID],             // owner
 435                         nei[curFaceID],             // neighbour
 436                         false,                      // flip flux
 437                         -1,                         // patch for face
 438                         false,                      // remove from zone
 439                         -1,                         // zone for face
 440                         false                       // face zone flip
 441                     )
 442                 );
 443 // Pout << "modifying stick-out face. Internal Old face: " << oldFace << " new face: " << newFace << " own: " << own[curFaceID] << " nei: " << nei[curFaceID] << endl;
 444             }
 445             else
 446             {
 447                 ref.setAction
 448                 (
 449                     polyModifyFace
 450                     (
 451                         newFace,                     // face
 452                         curFaceID,                   // master face
 453                         own[curFaceID],              // owner
 454                         -1,                          // neighbour
 455                         false,                       // flip flux
 456                         mesh.boundaryMesh().whichPatch(curFaceID), // patch
 457                         false,                        // remove from zone
 458                         -1,                           // zone for face
 459                         false                         // face zone flip
 460                     )
 461                 );
 462 // Pout << "modifying stick-out face. Boundary Old face: " << oldFace << " new face: " << newFace << " own: " << own[curFaceID] << " patch: " << mesh.boundaryMesh().whichPatch(curFaceID) << endl;
 463             }
 464         }
 465     }
 466
 467     if (debug)
 468     {
 469         Pout<< "void attachDetach::detachInterface("
 470             << "polyTopoChange& ref) const "
 471             << " for object " << name() << " : "
 472             << "Finished detaching interface" << endl;
 473     }
 474 }
 475
 476
 477 // ************************************************************************* //