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