src/dynamicMesh/polyTopoChange/polyTopoChange/removeCells.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 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 "removeCells.H"
  28 #include "polyMesh.H"
  29 #include "polyTopoChange.H"
  30 #include "polyRemoveCell.H"
  31 #include "polyRemoveFace.H"
  32 #include "polyModifyFace.H"
  33 #include "polyRemovePoint.H"
  34 #include "syncTools.H"
  35
  36 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  37
  38 namespace Foam
  39 {
  40
  41 defineTypeNameAndDebug(removeCells, 0);
  42
  43 }
  44
  45 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  46
  47 // Remove count of elements of f.
  48 void Foam::removeCells::uncount
  49 (
  50     const labelList& f,
  51     labelList& nUsage
  52 )
  53 {
  54     forAll(f, fp)
  55     {
  56         nUsage[f[fp]]--;
  57     }
  58 }
  59
  60
  61 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  62
  63 // Construct from mesh
  64 Foam::removeCells::removeCells
  65 (
  66     const polyMesh& mesh,
  67     const bool syncPar
  68 )
  69 :
  70     mesh_(mesh),
  71     syncPar_(syncPar)
  72 {}
  73
  74
  75 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  76
  77 //- Get labels of exposed faces. These are
  78 //  - internal faces that become boundary faces
  79 //  - coupled faces that become uncoupled (since on of the sides
  80 //    gets deleted)
  81 Foam::labelList Foam::removeCells::getExposedFaces
  82 (
  83     const labelList& cellLabels
  84 ) const
  85 {
  86     // Create list of cells to be removed
  87     boolList removedCell(mesh_.nCells(), false);
  88
  89     // Go from labelList of cells-to-remove to a boolList.
  90     forAll(cellLabels, i)
  91     {
  92         removedCell[cellLabels[i]] = true;
  93     }
  94
  95
  96     const labelList& faceOwner = mesh_.faceOwner();
  97     const labelList& faceNeighbour = mesh_.faceNeighbour();
  98
  99     // Count cells using face.
 100     labelList nCellsUsingFace(mesh_.nFaces(), 0);
 101
 102     for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
 103     {
 104         label own = faceOwner[faceI];
 105         label nei = faceNeighbour[faceI];
 106
 107         if (!removedCell[own])
 108         {
 109             nCellsUsingFace[faceI]++;
 110         }
 111         if (!removedCell[nei])
 112         {
 113             nCellsUsingFace[faceI]++;
 114         }
 115     }
 116
 117     for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
 118     {
 119         label own = faceOwner[faceI];
 120
 121         if (!removedCell[own])
 122         {
 123             nCellsUsingFace[faceI]++;
 124         }
 125     }
 126
 127     // Coupled faces: add number of cells using face across couple.
 128     if (syncPar_)
 129     {
 130         syncTools::syncFaceList
 131         (
 132             mesh_,
 133             nCellsUsingFace,
 134             plusEqOp<label>(),
 135             false
 136         );
 137     }
 138
 139     // Now nCellsUsingFace:
 140     // 0 : internal face whose both cells get deleted
 141     //     boundary face whose all cells get deleted
 142     // 1 : internal face that gets exposed
 143     //     unaffected (uncoupled) boundary face
 144     //     coupled boundary face that gets exposed ('uncoupled')
 145     // 2 : unaffected internal face
 146     //     unaffected coupled boundary face
 147
 148     DynamicList<label> exposedFaces(mesh_.nFaces()/10);
 149
 150     for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
 151     {
 152         if (nCellsUsingFace[faceI] == 1)
 153         {
 154             exposedFaces.append(faceI);
 155         }
 156     }
 157
 158     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 159
 160     forAll(patches, patchI)
 161     {
 162         const polyPatch& pp = patches[patchI];
 163
 164         if (pp.coupled())
 165         {
 166             label faceI = pp.start();
 167
 168             forAll(pp, i)
 169             {
 170                 label own = faceOwner[faceI];
 171
 172                 if (nCellsUsingFace[faceI] == 1 && !removedCell[own])
 173                 {
 174                     // My owner not removed but other side is so has to become
 175                     // normal, uncoupled, boundary face
 176                     exposedFaces.append(faceI);
 177                 }
 178
 179                 faceI++;
 180             }
 181         }
 182     }
 183
 184     return exposedFaces.shrink();
 185 }
 186
 187
 188 void Foam::removeCells::setRefinement
 189 (
 190     const labelList& cellLabels,
 191     const labelList& exposedFaceLabels,
 192     const labelList& exposedPatchIDs,
 193     polyTopoChange& meshMod
 194 ) const
 195 {
 196     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 197
 198     if (exposedFaceLabels.size() != exposedPatchIDs.size())
 199     {
 200         FatalErrorIn
 201         (
 202             "removeCells::setRefinement(const labelList&"
 203             ", const labelList&, const labelList&, polyTopoChange&)"
 204         )   << "Size of exposedFaceLabels " << exposedFaceLabels.size()
 205             << " differs from size of exposedPatchIDs "
 206             << exposedPatchIDs.size()
 207             << abort(FatalError);
 208     }
 209
 210     // List of new patchIDs
 211     labelList newPatchID(mesh_.nFaces(), -1);
 212
 213     forAll(exposedFaceLabels, i)
 214     {
 215         label patchI = exposedPatchIDs[i];
 216
 217         if (patchI < 0 || patchI >= patches.size())
 218         {
 219             FatalErrorIn
 220             (
 221                 "removeCells::setRefinement(const labelList&"
 222                 ", const labelList&, const labelList&, polyTopoChange&)"
 223             )   << "Invalid patch " << patchI
 224                 << " for exposed face " << exposedFaceLabels[i] << endl
 225                 << "Valid patches 0.." << patches.size()-1
 226                 << abort(FatalError);
 227         }
 228
 229         if (patches[patchI].coupled())
 230         {
 231             FatalErrorIn
 232             (
 233                 "removeCells::setRefinement(const labelList&"
 234                 ", const labelList&, const labelList&, polyTopoChange&)"
 235             )   << "Trying to put exposed face " << exposedFaceLabels[i]
 236                 << " into a coupled patch : " << patches[patchI].name()
 237                 << endl
 238                 << "This is illegal."
 239                 << abort(FatalError);
 240         }
 241
 242         newPatchID[exposedFaceLabels[i]] = patchI;
 243     }
 244
 245
 246     // Create list of cells to be removed
 247     boolList removedCell(mesh_.nCells(), false);
 248
 249     // Go from labelList of cells-to-remove to a boolList and remove all
 250     // cells mentioned.
 251     forAll(cellLabels, i)
 252     {
 253         label cellI = cellLabels[i];
 254
 255         removedCell[cellI] = true;
 256
 257         //Pout<< "Removing cell " << cellI
 258         //    << " cc:" << mesh_.cellCentres()[cellI] << endl;
 259
 260         meshMod.setAction(polyRemoveCell(cellI));
 261     }
 262
 263
 264     // Remove faces that are no longer used. Modify faces that
 265     // are used by one cell only.
 266
 267     const faceList& faces = mesh_.faces();
 268     const labelList& faceOwner = mesh_.faceOwner();
 269     const labelList& faceNeighbour = mesh_.faceNeighbour();
 270     const faceZoneMesh& faceZones = mesh_.faceZones();
 271
 272     // Count starting number of faces using each point. Keep up to date whenever
 273     // removing a face.
 274     labelList nFacesUsingPoint(mesh_.nPoints(), 0);
 275
 276     forAll(faces, faceI)
 277     {
 278         const face& f = faces[faceI];
 279
 280         forAll(f, fp)
 281         {
 282             nFacesUsingPoint[f[fp]]++;
 283         }
 284     }
 285
 286
 287     for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
 288     {
 289         const face& f = faces[faceI];
 290         label own = faceOwner[faceI];
 291         label nei = faceNeighbour[faceI];
 292
 293         if (removedCell[own])
 294         {
 295             if (removedCell[nei])
 296             {
 297                 // Face no longer used
 298                 //Pout<< "Removing internal face " << faceI
 299                 //    << " fc:" << mesh_.faceCentres()[faceI] << endl;
 300
 301                 meshMod.setAction(polyRemoveFace(faceI));
 302                 uncount(f, nFacesUsingPoint);
 303             }
 304             else
 305             {
 306                 if (newPatchID[faceI] == -1)
 307                 {
 308                     FatalErrorIn
 309                     (
 310                         "removeCells::setRefinement(const labelList&"
 311                         ", const labelList&, const labelList&"
 312                         ", polyTopoChange&)"
 313                     )   << "No patchID provided for exposed face " << faceI
 314                         << " on cell " << nei << nl
 315                         << "Did you provide patch IDs for all exposed faces?"
 316                         << abort(FatalError);
 317                 }
 318
 319                 // nei is remaining cell. FaceI becomes external cell
 320
 321                 label zoneID = faceZones.whichZone(faceI);
 322                 bool zoneFlip = false;
 323
 324                 if (zoneID >= 0)
 325                 {
 326                     const faceZone& fZone = faceZones[zoneID];
 327                     // Note: we reverse the owner/neighbour of the face
 328                     // so should also select the other side of the zone
 329                     zoneFlip = !fZone.flipMap()[fZone.whichFace(faceI)];
 330                 }
 331
 332                 //Pout<< "Putting exposed internal face " << faceI
 333                 //    << " fc:" << mesh_.faceCentres()[faceI]
 334                 //    << " into patch " << newPatchID[faceI] << endl;
 335
 336                 meshMod.setAction
 337                 (
 338                     polyModifyFace
 339                     (
 340                         f.reverseFace(),        // modified face
 341                         faceI,                  // label of face being modified
 342                         nei,                    // owner
 343                         -1,                     // neighbour
 344                         true,                   // face flip
 345                         newPatchID[faceI],      // patch for face
 346                         false,                  // remove from zone
 347                         zoneID,                 // zone for face
 348                         zoneFlip                // face flip in zone
 349                     )
 350                 );
 351             }
 352         }
 353         else if (removedCell[nei])
 354         {
 355             if (newPatchID[faceI] == -1)
 356             {
 357                 FatalErrorIn
 358                 (
 359                     "removeCells::setRefinement(const labelList&"
 360                     ", const labelList&, const labelList&"
 361                     ", polyTopoChange&)"
 362                 )   << "No patchID provided for exposed face " << faceI
 363                     << " on cell " << own << nl
 364                     << "Did you provide patch IDs for all exposed faces?"
 365                     << abort(FatalError);
 366             }
 367
 368             //Pout<< "Putting exposed internal face " << faceI
 369             //    << " fc:" << mesh_.faceCentres()[faceI]
 370             //    << " into patch " << newPatchID[faceI] << endl;
 371
 372             // own is remaining cell. FaceI becomes external cell.
 373             label zoneID = faceZones.whichZone(faceI);
 374             bool zoneFlip = false;
 375
 376             if (zoneID >= 0)
 377             {
 378                 const faceZone& fZone = faceZones[zoneID];
 379                 zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
 380             }
 381
 382             meshMod.setAction
 383             (
 384                 polyModifyFace
 385                 (
 386                     f,                      // modified face
 387                     faceI,                  // label of face being modified
 388                     own,                    // owner
 389                     -1,                     // neighbour
 390                     false,                  // face flip
 391                     newPatchID[faceI],      // patch for face
 392                     false,                  // remove from zone
 393                     zoneID,                 // zone for face
 394                     zoneFlip                // face flip in zone
 395                 )
 396             );
 397         }
 398     }
 399
 400     forAll(patches, patchI)
 401     {
 402         const polyPatch& pp = patches[patchI];
 403
 404         if (pp.coupled())
 405         {
 406             label faceI = pp.start();
 407
 408             forAll(pp, i)
 409             {
 410                 if (newPatchID[faceI] != -1)
 411                 {
 412                     //Pout<< "Putting uncoupled coupled face " << faceI
 413                     //    << " fc:" << mesh_.faceCentres()[faceI]
 414                     //    << " into patch " << newPatchID[faceI] << endl;
 415
 416                     label zoneID = faceZones.whichZone(faceI);
 417                     bool zoneFlip = false;
 418
 419                     if (zoneID >= 0)
 420                     {
 421                         const faceZone& fZone = faceZones[zoneID];
 422                         zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
 423                     }
 424
 425                     meshMod.setAction
 426                     (
 427                         polyModifyFace
 428                         (
 429                             faces[faceI],           // modified face
 430                             faceI,                  // label of face
 431                             faceOwner[faceI],       // owner
 432                             -1,                     // neighbour
 433                             false,                  // face flip
 434                             newPatchID[faceI],      // patch for face
 435                             false,                  // remove from zone
 436                             zoneID,                 // zone for face
 437                             zoneFlip                // face flip in zone
 438                         )
 439                     );
 440                 }
 441                 else if (removedCell[faceOwner[faceI]])
 442                 {
 443                     // Face no longer used
 444                     //Pout<< "Removing boundary face " << faceI
 445                     //    << " fc:" << mesh_.faceCentres()[faceI]
 446                     //    << endl;
 447
 448                     meshMod.setAction(polyRemoveFace(faceI));
 449                     uncount(faces[faceI], nFacesUsingPoint);
 450                 }
 451
 452                 faceI++;
 453             }
 454         }
 455         else
 456         {
 457             label faceI = pp.start();
 458
 459             forAll(pp, i)
 460             {
 461                 if (newPatchID[faceI] != -1)
 462                 {
 463                     FatalErrorIn
 464                     (
 465                         "removeCells::setRefinement(const labelList&"
 466                         ", const labelList&, const labelList&"
 467                         ", polyTopoChange&)"
 468                     )   << "new patchID provided for boundary face " << faceI
 469                         << " even though it is not on a coupled face."
 470                         << abort(FatalError);
 471                 }
 472
 473                 if (removedCell[faceOwner[faceI]])
 474                 {
 475                     // Face no longer used
 476                     //Pout<< "Removing boundary face " << faceI
 477                     //    << " fc:" << mesh_.faceCentres()[faceI]
 478                     //    << endl;
 479
 480                     meshMod.setAction(polyRemoveFace(faceI));
 481                     uncount(faces[faceI], nFacesUsingPoint);
 482                 }
 483
 484                 faceI++;
 485             }
 486         }
 487     }
 488
 489
 490     // Remove points that are no longer used.
 491     // Loop rewritten to not use pointFaces.
 492
 493     forAll(nFacesUsingPoint, pointI)
 494     {
 495         if (nFacesUsingPoint[pointI] == 0)
 496         {
 497             //Pout<< "Removing unused point " << pointI
 498             //    << " at:" << mesh_.points()[pointI] << endl;
 499
 500             meshMod.setAction(polyRemovePoint(pointI));
 501         }
 502         else if (nFacesUsingPoint[pointI] == 1)
 503         {
 504             WarningIn
 505             (
 506                 "removeCells::setRefinement(const labelList&"
 507                 ", const labelList&, const labelList&"
 508                 ", polyTopoChange&)"
 509             )   << "point " << pointI << " at coordinate "
 510                 << mesh_.points()[pointI]
 511                 << " is only used by 1 face after removing cells."
 512                 << " This probably results in an illegal mesh."
 513                 << endl;
 514         }
 515     }
 516 }
 517
 518
 519 // ************************************************************************* //