applications/utilities/postProcessing/miscellaneous/postChannel/channelIndex.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 "channelIndex.H"
  28 #include <OpenFOAM/boolList.H>
  29 #include <OpenFOAM/syncTools.H>
  30 #include <OpenFOAM/OFstream.H>
  31 #include <meshTools/meshTools.H>
  32 #include <OpenFOAM/Time.H>
  33 #include <OpenFOAM/SortableList.H>
  34
  35 // * * * * * * * * * * * * * Static Member Data  * * * * * * * * * * * * * * //
  36
  37 template<>
  38 const char* Foam::NamedEnum<Foam::vector::components, 3>::names[] =
  39 {
  40     "x",
  41     "y",
  42     "z"
  43 };
  44
  45 const Foam::NamedEnum<Foam::vector::components, 3>
  46     Foam::channelIndex::vectorComponentsNames_;
  47
  48
  49 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  50
  51 // Determines face blocking
  52 void Foam::channelIndex::walkOppositeFaces
  53 (
  54     const polyMesh& mesh,
  55     const labelList& startFaces,
  56     boolList& blockedFace
  57 )
  58 {
  59     const cellList& cells = mesh.cells();
  60     const faceList& faces = mesh.faces();
  61     label nBnd = mesh.nFaces() - mesh.nInternalFaces();
  62
  63     DynamicList<label> frontFaces(startFaces);
  64     forAll(frontFaces, i)
  65     {
  66         label faceI = frontFaces[i];
  67         blockedFace[faceI] = true;
  68     }
  69
  70     while (returnReduce(frontFaces.size(), sumOp<label>()) > 0)
  71     {
  72         // Transfer across.
  73         boolList isFrontBndFace(nBnd, false);
  74         forAll(frontFaces, i)
  75         {
  76             label faceI = frontFaces[i];
  77
  78             if (!mesh.isInternalFace(faceI))
  79             {
  80                 isFrontBndFace[faceI-mesh.nInternalFaces()] = true;
  81             }
  82         }
  83         syncTools::swapBoundaryFaceList(mesh, isFrontBndFace, false);
  84
  85         // Add
  86         forAll(isFrontBndFace, i)
  87         {
  88             label faceI = mesh.nInternalFaces()+i;
  89             if (isFrontBndFace[i] && !blockedFace[faceI])
  90             {
  91                 blockedFace[faceI] = true;
  92                 frontFaces.append(faceI);
  93             }
  94         }
  95
  96         // Transfer across cells
  97         DynamicList<label> newFrontFaces(frontFaces.size());
  98
  99         forAll(frontFaces, i)
 100         {
 101             label faceI = frontFaces[i];
 102
 103             {
 104                 const cell& ownCell = cells[mesh.faceOwner()[faceI]];
 105
 106                 label oppositeFaceI = ownCell.opposingFaceLabel(faceI, faces);
 107
 108                 if (oppositeFaceI == -1)
 109                 {
 110                     FatalErrorIn("channelIndex::walkOppositeFaces(..)")
 111                         << "Face:" << faceI << " owner cell:" << ownCell
 112                         << " is not a hex?" << abort(FatalError);
 113                 }
 114                 else
 115                 {
 116                     if (!blockedFace[oppositeFaceI])
 117                     {
 118                         blockedFace[oppositeFaceI] = true;
 119                         newFrontFaces.append(oppositeFaceI);
 120                     }
 121                 }
 122             }
 123
 124             if (mesh.isInternalFace(faceI))
 125             {
 126                 const cell& neiCell = mesh.cells()[mesh.faceNeighbour()[faceI]];
 127
 128                 label oppositeFaceI = neiCell.opposingFaceLabel(faceI, faces);
 129
 130                 if (oppositeFaceI == -1)
 131                 {
 132                     FatalErrorIn("channelIndex::walkOppositeFaces(..)")
 133                         << "Face:" << faceI << " neighbour cell:" << neiCell
 134                         << " is not a hex?" << abort(FatalError);
 135                 }
 136                 else
 137                 {
 138                     if (!blockedFace[oppositeFaceI])
 139                     {
 140                         blockedFace[oppositeFaceI] = true;
 141                         newFrontFaces.append(oppositeFaceI);
 142                     }
 143                 }
 144             }
 145         }
 146
 147         frontFaces.transfer(newFrontFaces);
 148     }
 149 }
 150
 151
 152 // Calculate regions.
 153 void Foam::channelIndex::calcLayeredRegions
 154 (
 155     const polyMesh& mesh,
 156     const labelList& startFaces
 157 )
 158 {
 159     boolList blockedFace(mesh.nFaces(), false);
 160     walkOppositeFaces
 161     (
 162         mesh,
 163         startFaces,
 164         blockedFace
 165     );
 166
 167
 168     if (false)
 169     {
 170         OFstream str(mesh.time().path()/"blockedFaces.obj");
 171         label vertI = 0;
 172         forAll(blockedFace, faceI)
 173         {
 174             if (blockedFace[faceI])
 175             {
 176                 const face& f = mesh.faces()[faceI];
 177                 forAll(f, fp)
 178                 {
 179                     meshTools::writeOBJ(str, mesh.points()[f[fp]]);
 180                 }
 181                 str<< 'f';
 182                 forAll(f, fp)
 183                 {
 184                     str << ' ' << vertI+fp+1;
 185                 }
 186                 str << nl;
 187                 vertI += f.size();
 188             }
 189         }
 190     }
 191
 192
 193     // Do analysis for connected regions
 194     cellRegion_.reset(new regionSplit(mesh, blockedFace));
 195
 196     Info<< "Detected " << cellRegion_().nRegions() << " layers." << nl << endl;
 197
 198     // Sum number of entries per region
 199     regionCount_ = regionSum(scalarField(mesh.nCells(), 1.0));
 200
 201     // Average cell centres to determine ordering.
 202     pointField regionCc
 203     (
 204         regionSum(mesh.cellCentres())
 205       / regionCount_
 206     );
 207
 208     SortableList<scalar> sortComponent(regionCc.component(dir_));
 209
 210     sortMap_ = sortComponent.indices();
 211
 212     y_ = sortComponent;
 213
 214     if (symmetric_)
 215     {
 216         y_.setSize(cellRegion_().nRegions()/2);
 217     }
 218 }
 219
 220
 221 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 222
 223 Foam::channelIndex::channelIndex
 224 (
 225     const polyMesh& mesh,
 226     const dictionary& dict
 227 )
 228 :
 229     symmetric_(readBool(dict.lookup("symmetric"))),
 230     dir_(vectorComponentsNames_.read(dict.lookup("component")))
 231 {
 232     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 233
 234     const wordList patchNames(dict.lookup("patches"));
 235
 236     label nFaces = 0;
 237
 238     forAll(patchNames, i)
 239     {
 240         label patchI = patches.findPatchID(patchNames[i]);
 241
 242         if (patchI == -1)
 243         {
 244             FatalErrorIn("channelIndex::channelIndex(const polyMesh&)")
 245                 << "Illegal patch " << patchNames[i]
 246                 << ". Valid patches are " << patches.name()
 247                 << exit(FatalError);
 248         }
 249
 250         nFaces += patches[patchI].size();
 251     }
 252
 253     labelList startFaces(nFaces);
 254     nFaces = 0;
 255
 256     forAll(patchNames, i)
 257     {
 258         const polyPatch& pp = patches[patches.findPatchID(patchNames[i])];
 259
 260         forAll(pp, j)
 261         {
 262             startFaces[nFaces++] = pp.start()+j;
 263         }
 264     }
 265
 266     // Calculate regions.
 267     calcLayeredRegions(mesh, startFaces);
 268 }
 269
 270
 271 Foam::channelIndex::channelIndex
 272 (
 273     const polyMesh& mesh,
 274     const labelList& startFaces,
 275     const bool symmetric,
 276     const direction dir
 277 )
 278 :
 279     symmetric_(symmetric),
 280     dir_(dir)
 281 {
 282     // Calculate regions.
 283     calcLayeredRegions(mesh, startFaces);
 284 }
 285
 286
 287 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 288
 289
 290 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 291
 292
 293 // ************************ vim: set sw=4 sts=4 et: ************************ //