src/meshLibrary/utilities/boundaryLayers/boundaryLayersCheckTopologyOfBndFaces.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 "boundaryLayers.H"
  29 #include "meshSurfaceEngine.H"
  30 #include "decomposeCells.H"
  31 #include "helperFunctions.H"
  32 #include "HashSet.H"
  33
  34 #include <set>
  35
  36 //#define DEBUGLayer
  37
  38 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  39
  40 namespace Foam
  41 {
  42
  43 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  44
  45 void boundaryLayers::checkTopologyOfBoundaryFaces(const labelList& patchLabels)
  46 {
  47     if( !patchWiseLayers_ )
  48         return;
  49
  50     Info << "Checking topology of boundary faces" << endl;
  51
  52     labelHashSet usedPatches;
  53     forAll(patchLabels, i)
  54         usedPatches.insert(patchLabels[i]);
  55
  56     //- create a set of patch pairs. These are pairs at which the layers
  57     //- shall be terminated
  58     std::set<std::pair<label, label> > terminatedPairs;
  59     forAll(treatPatchesWithPatch_, patchI)
  60     {
  61         const DynList<label>& otherPatches = treatPatchesWithPatch_[patchI];
  62
  63         forAll(otherPatches, patchJ)
  64         {
  65             if( patchI == otherPatches[patchJ] )
  66                 continue;
  67
  68             terminatedPairs.insert
  69             (
  70                 std::make_pair
  71                 (
  72                     Foam::min(patchI, otherPatches[patchJ]),
  73                     Foam::max(patchI, otherPatches[patchJ])
  74                 )
  75             );
  76         }
  77     }
  78
  79     bool changed;
  80     label nDecomposed(0);
  81     boolList decomposeCell(mesh_.cells().size(), false);
  82
  83     do
  84     {
  85         changed = false;
  86
  87         const meshSurfaceEngine& mse = this->surfaceEngine();
  88         const faceList::subList& bFaces = mse.boundaryFaces();
  89         const labelList& faceOwner = mse.faceOwners();
  90         const labelList& facePatches = mse.boundaryFacePatches();
  91         const VRWGraph& faceEdges = mse.faceEdges();
  92         const VRWGraph& edgeFaces = mse.edgeFaces();
  93
  94         const Map<label>& otherProcPatch = mse.otherEdgeFacePatch();
  95
  96         VRWGraph newBoundaryFaces;
  97         labelLongList newBoundaryOwners;
  98         labelLongList newBoundaryPatches;
  99
 100         forAll(bFaces, bfI)
 101         {
 102             const face& bf = bFaces[bfI];
 103             const label fPatch = facePatches[bfI];
 104
 105             if( !usedPatches.found(fPatch) )
 106                 continue;
 107
 108             //- find patches of neighbour faces
 109             labelList neiPatches(bf.size());
 110             forAll(bf, eI)
 111             {
 112                 const label beI = faceEdges(bfI, eI);
 113
 114                 if( edgeFaces.sizeOfRow(beI) == 2 )
 115                 {
 116                     label neiFace = edgeFaces(beI, 0);
 117                     if( neiFace == bfI )
 118                         neiFace = edgeFaces(beI, 1);
 119
 120                     neiPatches[eI] = facePatches[neiFace];
 121                 }
 122                 else if( edgeFaces.sizeOfRow(beI) == 1 )
 123                 {
 124                     //- edge is at a parallel boundary
 125                     neiPatches[eI] = otherProcPatch[beI];
 126                 }
 127             }
 128
 129             //- find feature edges and check if the patches meeting there
 130             //- shall be treated together.
 131             bool storedFace(false);
 132             forAll(neiPatches, eI)
 133             {
 134                 if( neiPatches[eI] == fPatch )
 135                     continue;
 136
 137                 std::pair<label, label> pp
 138                 (
 139                     Foam::min(fPatch, neiPatches[eI]),
 140                     Foam::max(fPatch, neiPatches[eI])
 141                 );
 142
 143                 if( terminatedPairs.find(pp) == terminatedPairs.end() )
 144                     continue;
 145
 146                 //- create a new face from this edge and the neighbouring edges
 147                 bool usePrev(false), useNext(false);
 148                 if( neiPatches[neiPatches.rcIndex(eI)] == fPatch )
 149                 {
 150                     usePrev = true;
 151                 }
 152                 else
 153                 {
 154                     std::pair<label, label> ppPrev
 155                     (
 156                         Foam::min(fPatch, neiPatches[neiPatches.rcIndex(eI)]),
 157                         Foam::max(fPatch, neiPatches[neiPatches.rcIndex(eI)])
 158                     );
 159
 160                     if( terminatedPairs.find(ppPrev) == terminatedPairs.end() )
 161                         usePrev = true;
 162                 }
 163
 164                 if( neiPatches[neiPatches.fcIndex(eI)] == fPatch )
 165                 {
 166                     useNext = true;
 167                 }
 168                 else
 169                 {
 170                     std::pair<label, label> ppNext
 171                     (
 172                         Foam::min(fPatch, neiPatches[neiPatches.fcIndex(eI)]),
 173                         Foam::max(fPatch, neiPatches[neiPatches.fcIndex(eI)])
 174                     );
 175
 176                     if( terminatedPairs.find(ppNext) == terminatedPairs.end() )
 177                         useNext = true;
 178                 }
 179
 180                 DynList<edge> removeEdges;
 181                 if( useNext && usePrev )
 182                 {
 183                     removeEdges.setSize(3);
 184                     removeEdges[0] = bf.faceEdge(neiPatches.rcIndex(eI));
 185                     removeEdges[1] = bf.faceEdge(eI);
 186                     removeEdges[2] = bf.faceEdge(neiPatches.fcIndex(eI));
 187                 }
 188                 else if( useNext )
 189                 {
 190                     removeEdges.setSize(2);
 191                     removeEdges[0] = bf.faceEdge(neiPatches.fcIndex(eI));
 192                     removeEdges[1] = bf.faceEdge(eI);
 193                 }
 194                 else if( usePrev )
 195                 {
 196                     removeEdges.setSize(2);
 197                     removeEdges[0] = bf.faceEdge(neiPatches.rcIndex(eI));
 198                     removeEdges[1] = bf.faceEdge(eI);
 199                 }
 200
 201                 const face cutFace = help::removeEdgesFromFace(bf, removeEdges);
 202                 if( cutFace.size() > 2 )
 203                 {
 204                     newBoundaryFaces.appendList(cutFace);
 205                     newBoundaryOwners.append(faceOwner[bfI]);
 206                     newBoundaryPatches.append(fPatch);
 207                 }
 208                 const face rFace = help::createFaceFromRemovedPart(bf, cutFace);
 209                 if( rFace.size() > 2 )
 210                 {
 211                     newBoundaryFaces.appendList(rFace);
 212                     newBoundaryOwners.append(faceOwner[bfI]);
 213                     newBoundaryPatches.append(fPatch);
 214                 }
 215
 216                 if( (cutFace.size() > 2) && (rFace.size() > 2) )
 217                 {
 218                     decomposeCell[faceOwner[bfI]] = true;
 219                     changed = true;
 220                     ++nDecomposed;
 221                 }
 222
 223                 storedFace = true;
 224
 225                 break;
 226             }
 227
 228             if( !storedFace )
 229             {
 230                 newBoundaryFaces.appendList(bf);
 231                 newBoundaryOwners.append(faceOwner[bfI]);
 232                 newBoundaryPatches.append(fPatch);
 233             }
 234         }
 235
 236         //- Finally, replace the boundary faces
 237         reduce(changed, maxOp<bool>());
 238
 239         if( changed )
 240         {
 241             polyMeshGenModifier(mesh_).replaceBoundary
 242             (
 243                 patchNames_,
 244                 newBoundaryFaces,
 245                 newBoundaryOwners,
 246                 newBoundaryPatches
 247             );
 248
 249             clearOut();
 250         }
 251
 252     } while( changed );
 253
 254     //- decompose owner cells adjacent to the decomposed faces
 255     reduce(nDecomposed, sumOp<label>());
 256
 257     if( nDecomposed != 0 )
 258     {
 259         FatalError << "Critical. Not tested" << exit(FatalError);
 260         decomposeCells dc(mesh_);
 261         dc.decomposeMesh(decomposeCell);
 262
 263         clearOut();
 264     }
 265
 266     mesh_.write();
 267     Info << "Finished checking topology" << endl;
 268 }
 269
 270 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 271
 272 } // End namespace Foam
 273
 274 // ************************************************************************* //