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-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 "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                     zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
 328                 }
 329
 330                 //Pout<< "Putting exposed internal face " << faceI
 331                 //    << " fc:" << mesh_.faceCentres()[faceI]
 332                 //    << " into patch " << newPatchID[faceI] << endl;
 333
 334                 meshMod.setAction
 335                 (
 336                     polyModifyFace
 337                     (
 338                         f.reverseFace(),        // modified face
 339                         faceI,                  // label of face being modified
 340                         nei,                    // owner
 341                         -1,                     // neighbour
 342                         false,                  // face flip
 343                         newPatchID[faceI],      // patch for face
 344                         false,                  // remove from zone
 345                         zoneID,                 // zone for face
 346                         zoneFlip                // face flip in zone
 347                     )
 348                 );
 349             }
 350         }
 351         else if (removedCell[nei])
 352         {
 353             if (newPatchID[faceI] == -1)
 354             {
 355                 FatalErrorIn
 356                 (
 357                     "removeCells::setRefinement(const labelList&"
 358                     ", const labelList&, const labelList&"
 359                     ", polyTopoChange&)"
 360                 )   << "No patchID provided for exposed face " << faceI
 361                     << " on cell " << own << nl
 362                     << "Did you provide patch IDs for all exposed faces?"
 363                     << abort(FatalError);
 364             }
 365
 366             //Pout<< "Putting exposed internal face " << faceI
 367             //    << " fc:" << mesh_.faceCentres()[faceI]
 368             //    << " into patch " << newPatchID[faceI] << endl;
 369
 370             // own is remaining cell. FaceI becomes external cell.
 371             label zoneID = faceZones.whichZone(faceI);
 372             bool zoneFlip = false;
 373
 374             if (zoneID >= 0)
 375             {
 376                 const faceZone& fZone = faceZones[zoneID];
 377                 zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
 378             }
 379
 380             meshMod.setAction
 381             (
 382                 polyModifyFace
 383                 (
 384                     f,                      // modified face
 385                     faceI,                  // label of face being modified
 386                     own,                    // owner
 387                     -1,                     // neighbour
 388                     false,                  // face flip
 389                     newPatchID[faceI],      // patch for face
 390                     false,                  // remove from zone
 391                     zoneID,                 // zone for face
 392                     zoneFlip                // face flip in zone
 393                 )
 394             );
 395         }
 396     }
 397
 398     forAll(patches, patchI)
 399     {
 400         const polyPatch& pp = patches[patchI];
 401
 402         if (pp.coupled())
 403         {
 404             label faceI = pp.start();
 405
 406             forAll(pp, i)
 407             {
 408                 if (newPatchID[faceI] != -1)
 409                 {
 410                     //Pout<< "Putting uncoupled coupled face " << faceI
 411                     //    << " fc:" << mesh_.faceCentres()[faceI]
 412                     //    << " into patch " << newPatchID[faceI] << endl;
 413
 414                     label zoneID = faceZones.whichZone(faceI);
 415                     bool zoneFlip = false;
 416
 417                     if (zoneID >= 0)
 418                     {
 419                         const faceZone& fZone = faceZones[zoneID];
 420                         zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
 421                     }
 422
 423                     meshMod.setAction
 424                     (
 425                         polyModifyFace
 426                         (
 427                             faces[faceI],           // modified face
 428                             faceI,                  // label of face
 429                             faceOwner[faceI],       // owner
 430                             -1,                     // neighbour
 431                             false,                  // face flip
 432                             newPatchID[faceI],      // patch for face
 433                             false,                  // remove from zone
 434                             zoneID,                 // zone for face
 435                             zoneFlip                // face flip in zone
 436                         )
 437                     );
 438                 }
 439                 else if (removedCell[faceOwner[faceI]])
 440                 {
 441                     // Face no longer used
 442                     //Pout<< "Removing boundary face " << faceI
 443                     //    << " fc:" << mesh_.faceCentres()[faceI]
 444                     //    << endl;
 445
 446                     meshMod.setAction(polyRemoveFace(faceI));
 447                     uncount(faces[faceI], nFacesUsingPoint);
 448                 }
 449
 450                 faceI++;
 451             }
 452         }
 453         else
 454         {
 455             label faceI = pp.start();
 456
 457             forAll(pp, i)
 458             {
 459                 if (newPatchID[faceI] != -1)
 460                 {
 461                     FatalErrorIn
 462                     (
 463                         "removeCells::setRefinement(const labelList&"
 464                         ", const labelList&, const labelList&"
 465                         ", polyTopoChange&)"
 466                     )   << "new patchID provided for boundary face " << faceI
 467                         << " even though it is not on a coupled face."
 468                         << abort(FatalError);
 469                 }
 470
 471                 if (removedCell[faceOwner[faceI]])
 472                 {
 473                     // Face no longer used
 474                     //Pout<< "Removing boundary face " << faceI
 475                     //    << " fc:" << mesh_.faceCentres()[faceI]
 476                     //    << endl;
 477
 478                     meshMod.setAction(polyRemoveFace(faceI));
 479                     uncount(faces[faceI], nFacesUsingPoint);
 480                 }
 481
 482                 faceI++;
 483             }
 484         }
 485     }
 486
 487
 488     // Remove points that are no longer used.
 489     // Loop rewritten to not use pointFaces.
 490
 491     forAll(nFacesUsingPoint, pointI)
 492     {
 493         if (nFacesUsingPoint[pointI] == 0)
 494         {
 495             //Pout<< "Removing unused point " << pointI
 496             //    << " at:" << mesh_.points()[pointI] << endl;
 497
 498             meshMod.setAction(polyRemovePoint(pointI));
 499         }
 500         else if (nFacesUsingPoint[pointI] == 1)
 501         {
 502             WarningIn
 503             (
 504                 "removeCells::setRefinement(const labelList&"
 505                 ", const labelList&, const labelList&"
 506                 ", polyTopoChange&)"
 507             )   << "point " << pointI << " at coordinate "
 508                 << mesh_.points()[pointI]
 509                 << " is only used by 1 face after removing cells."
 510                 << " This probably results in an illegal mesh."
 511                 << endl;
 512         }
 513     }
 514 }
 515
 516
 517 // ************************************************************************* //