src/meshLibrary/utilities/decomposeCells/decomposeCellsDecomposition.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | cfMesh: A library for mesh generation
   4    \\    /   O peration     |
   5     \\  /    A nd           | Author: Franjo Juretic (franjo.juretic@c-fields.com)
   6      \\/     M anipulation  | Copyright (C) Creative Fields, Ltd.
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of cfMesh.
  10
  11     cfMesh 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     cfMesh 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 cfMesh.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Description
  25
  26 \*---------------------------------------------------------------------------*/
  27
  28 #include "decomposeCells.H"
  29 #include "demandDrivenData.H"
  30 #include "polyMeshGenAddressing.H"
  31 #include "meshSurfaceEngine.H"
  32 #include "decomposeFaces.H"
  33 #include "labelLongList.H"
  34
  35 //#define DEBUGDecompose
  36
  37 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  38
  39 namespace Foam
  40 {
  41
  42 void decomposeCells::decomposeMesh(const boolList& decomposeCell)
  43 {
  44     checkFaceConnections(decomposeCell);
  45
  46     createPointsAndCellFaces(decomposeCell);
  47
  48     storeBoundaryFaces(decomposeCell);
  49
  50     removeDecomposedCells(decomposeCell);
  51
  52     addNewCells();
  53
  54     # ifdef DEBUGDecompose
  55     mesh_.addressingData().checkMesh();
  56     # endif
  57 }
  58
  59 void decomposeCells::checkFaceConnections(const boolList& decomposeCell)
  60 {
  61     const faceListPMG& faces = mesh_.faces();
  62     const cellListPMG& cells = mesh_.cells();
  63
  64     boolList decomposeFace(faces.size(), false);
  65     forAll(cells, cellI)
  66     {
  67         if( decomposeCell[cellI] )
  68         {
  69             DynList<label, 32> vrt;
  70             DynList<edge, 64> edges;
  71             DynList<DynList<label, 8> > faceEdges;
  72             DynList<DynList<label, 2>, 64> edgeFaces;
  73
  74             findAddressingForCell(cellI, vrt, edges, faceEdges, edgeFaces);
  75
  76             forAll(faceEdges, fI)
  77             {
  78                 const DynList<label, 8>& fEdges = faceEdges[fI];
  79
  80                 labelHashSet neiFaces;
  81                 forAll(fEdges, feI)
  82                 {
  83                     label neiFace = edgeFaces[fEdges[feI]][0];
  84                     if( neiFace == fI )
  85                         neiFace = edgeFaces[fEdges[feI]][1];
  86
  87                     if( neiFaces.found(neiFace) )
  88                     {
  89                         decomposeFace[cells[cellI][fI]] = true;
  90                     }
  91                     else
  92                     {
  93                         neiFaces.insert(neiFace);
  94                     }
  95                 }
  96             }
  97         }
  98     }
  99
 100     if( Pstream::parRun() )
 101     {
 102         const PtrList<processorBoundaryPatch>& procBoundaries =
 103             mesh_.procBoundaries();
 104
 105         //- send information to the neighbour processor
 106         forAll(procBoundaries, patchI)
 107         {
 108             const label start = procBoundaries[patchI].patchStart();
 109             boolList decFace(procBoundaries[patchI].patchSize(), false);
 110             const label size = decFace.size();
 111
 112             for(label i=0;i<size;++i)
 113             {
 114                 if( decomposeFace[start+i] )
 115                     decFace[i] = true;
 116             }
 117
 118             OPstream toOtherProc
 119             (
 120                 Pstream::blocking,
 121                 procBoundaries[patchI].neiProcNo(),
 122                 decFace.byteSize()
 123             );
 124
 125             toOtherProc << decFace;
 126         }
 127
 128         //- receive information from the neighbour processor
 129         forAll(procBoundaries, patchI)
 130         {
 131             boolList decFace;
 132
 133             IPstream fromOtherProc
 134             (
 135                 Pstream::blocking,
 136                 procBoundaries[patchI].neiProcNo()
 137             );
 138
 139             fromOtherProc >> decFace;
 140
 141             const label start = procBoundaries[patchI].patchStart();
 142             forAll(decFace, i)
 143             {
 144                 if( decFace[i] )
 145                     decomposeFace[start+i] = true;
 146             }
 147         }
 148     }
 149
 150     //- decompose faces which would cause invalid connections
 151     decomposeFaces(mesh_).decomposeMeshFaces(decomposeFace);
 152 }
 153
 154 void decomposeCells::createPointsAndCellFaces(const boolList& decomposeCell)
 155 {
 156     facesOfNewCells_.clear();
 157
 158     forAll(decomposeCell, cI)
 159         if( decomposeCell[cI] )
 160         {
 161             decomposeCellIntoPyramids(cI);
 162         }
 163 }
 164
 165 void decomposeCells::storeBoundaryFaces(const boolList& decomposeCell)
 166 {
 167     meshSurfaceEngine mse(mesh_);
 168     const faceList::subList& bFaces = mse.boundaryFaces();
 169     const labelList& facePatch = mse.boundaryFacePatches();
 170
 171     forAll(bFaces, bfI)
 172     {
 173         newBoundaryFaces_.appendList(bFaces[bfI]);
 174         newBoundaryPatches_.append(facePatch[bfI]);
 175     }
 176 }
 177
 178 void decomposeCells::removeDecomposedCells(const boolList& decomposeCell)
 179 {
 180     # ifdef DEBUGDecompose
 181     Info << "Number of cells before removal " << mesh_.cells().size() << endl;
 182     # endif
 183
 184     polyMeshGenModifier meshModifier(mesh_);
 185     meshModifier.removeCells(decomposeCell, false);
 186
 187     # ifdef DEBUGDecompose
 188     Info << "Number of cells after removal " << mesh_.cells().size() << endl;
 189     # endif
 190 }
 191
 192 void decomposeCells::addNewCells()
 193 {
 194     Info << "Adding new cells " << endl;
 195     polyMeshGenModifier(mesh_).addCells(facesOfNewCells_);
 196     facesOfNewCells_.clear();
 197     Info << "Reordering bnd faces" << endl;
 198     polyMeshGenModifier(mesh_).reorderBoundaryFaces();
 199
 200     Info << "Finding bnd faces" << endl;
 201     const faceListPMG& faces = mesh_.faces();
 202     const labelList& owner = mesh_.owner();
 203     const VRWGraph& pointFaces = mesh_.addressingData().pointFaces();
 204
 205     labelLongList newBoundaryOwners;
 206
 207     forAll(newBoundaryFaces_, faceI)
 208     {
 209         face bf(newBoundaryFaces_.sizeOfRow(faceI));
 210         forAllRow(newBoundaryFaces_, faceI, pI)
 211             bf[pI] = newBoundaryFaces_(faceI, pI);
 212
 213         # ifdef DEBUGDecompose
 214         Info << "Finding cell for boundary face " << bf << endl;
 215         bool found(false);
 216         forAllRow(pointFaces, bf[0], pfI)
 217             if( bf == faces[pointFaces(bf[0], pfI)] )
 218                 found = true;
 219         if( !found )
 220             FatalErrorIn
 221             (
 222                 "void decomposeCells::addNewCells()"
 223             ) << "Face " << bf << " does not exist in the mesh"
 224                 << abort(FatalError);
 225         #endif
 226
 227         forAllRow(pointFaces, bf[0], pfI)
 228         {
 229             const label fLabel = pointFaces(bf[0], pfI);
 230             if( (mesh_.faceIsInPatch(fLabel) != -1) && (bf == faces[fLabel]) )
 231             {
 232                 # ifdef DEBUGDecompose
 233                 Info << "Boundary face " << bf << " is in cell "
 234                 << owner[fLabel]] << endl;
 235                 # endif
 236
 237                 newBoundaryOwners.append(owner[fLabel]);
 238             }
 239         }
 240     }
 241
 242     polyMeshGenModifier(mesh_).replaceBoundary
 243     (
 244         patchNames_,
 245         newBoundaryFaces_,
 246         newBoundaryOwners,
 247         newBoundaryPatches_
 248     );
 249
 250     polyMeshGenModifier(mesh_).removeUnusedVertices();
 251     polyMeshGenModifier(mesh_).clearAll();
 252 }
 253
 254 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *//
 255
 256 } // End namespace Foam
 257
 258 // ************************************************************************* //