FreeFOAM The Cross-Platform CFD Toolkit
Hosted by SourceForge:
Get FreeFOAM at SourceForge.net.
            Fast, secure and Free Open Source software downloads

channelIndex.C

Go to the documentation of this file.
00001 /*---------------------------------------------------------------------------*\
00002   =========                 |
00003   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
00004    \\    /   O peration     |
00005     \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
00006      \\/     M anipulation  |
00007 -------------------------------------------------------------------------------
00008 License
00009     This file is part of OpenFOAM.
00010 
00011     OpenFOAM is free software: you can redistribute it and/or modify it
00012     under the terms of the GNU General Public License as published by
00013     the Free Software Foundation, either version 3 of the License, or
00014     (at your option) any later version.
00015 
00016     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
00017     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00019     for more details.
00020 
00021     You should have received a copy of the GNU General Public License
00022     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
00023 
00024 \*---------------------------------------------------------------------------*/
00025 
00026 #include "channelIndex.H"
00027 #include <OpenFOAM/boolList.H>
00028 #include <OpenFOAM/syncTools.H>
00029 #include <OpenFOAM/OFstream.H>
00030 #include <meshTools/meshTools.H>
00031 #include <OpenFOAM/Time.H>
00032 #include <OpenFOAM/SortableList.H>
00033 
00034 // * * * * * * * * * * * * * Static Member Data  * * * * * * * * * * * * * * //
00035 
00036 template<>
00037 const char* Foam::NamedEnum<Foam::vector::components, 3>::names[] =
00038 {
00039     "x",
00040     "y",
00041     "z"
00042 };
00043 
00044 const Foam::NamedEnum<Foam::vector::components, 3>
00045     Foam::channelIndex::vectorComponentsNames_;
00046 
00047 
00048 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
00049 
00050 // Determines face blocking
00051 void Foam::channelIndex::walkOppositeFaces
00052 (
00053     const polyMesh& mesh,
00054     const labelList& startFaces,
00055     boolList& blockedFace
00056 )
00057 {
00058     const cellList& cells = mesh.cells();
00059     const faceList& faces = mesh.faces();
00060     label nBnd = mesh.nFaces() - mesh.nInternalFaces();
00061 
00062     DynamicList<label> frontFaces(startFaces);
00063     forAll(frontFaces, i)
00064     {
00065         label faceI = frontFaces[i];
00066         blockedFace[faceI] = true;
00067     }
00068 
00069     while (returnReduce(frontFaces.size(), sumOp<label>()) > 0)
00070     {
00071         // Transfer across.
00072         boolList isFrontBndFace(nBnd, false);
00073         forAll(frontFaces, i)
00074         {
00075             label faceI = frontFaces[i];
00076 
00077             if (!mesh.isInternalFace(faceI))
00078             {
00079                 isFrontBndFace[faceI-mesh.nInternalFaces()] = true;
00080             }
00081         }
00082         syncTools::swapBoundaryFaceList(mesh, isFrontBndFace, false);
00083 
00084         // Add 
00085         forAll(isFrontBndFace, i)
00086         {
00087             label faceI = mesh.nInternalFaces()+i;
00088             if (isFrontBndFace[i] && !blockedFace[faceI])
00089             {
00090                 blockedFace[faceI] = true;
00091                 frontFaces.append(faceI);
00092             }
00093         }
00094 
00095         // Transfer across cells
00096         DynamicList<label> newFrontFaces(frontFaces.size());
00097 
00098         forAll(frontFaces, i)
00099         {
00100             label faceI = frontFaces[i];
00101 
00102             {
00103                 const cell& ownCell = cells[mesh.faceOwner()[faceI]];
00104 
00105                 label oppositeFaceI = ownCell.opposingFaceLabel(faceI, faces);
00106 
00107                 if (oppositeFaceI == -1)
00108                 {
00109                     FatalErrorIn("channelIndex::walkOppositeFaces(..)")
00110                         << "Face:" << faceI << " owner cell:" << ownCell
00111                         << " is not a hex?" << abort(FatalError);
00112                 }
00113                 else
00114                 {
00115                     if (!blockedFace[oppositeFaceI])
00116                     {
00117                         blockedFace[oppositeFaceI] = true;
00118                         newFrontFaces.append(oppositeFaceI);
00119                     }
00120                 }
00121             }
00122 
00123             if (mesh.isInternalFace(faceI))
00124             {
00125                 const cell& neiCell = mesh.cells()[mesh.faceNeighbour()[faceI]];
00126 
00127                 label oppositeFaceI = neiCell.opposingFaceLabel(faceI, faces);
00128 
00129                 if (oppositeFaceI == -1)
00130                 {
00131                     FatalErrorIn("channelIndex::walkOppositeFaces(..)")
00132                         << "Face:" << faceI << " neighbour cell:" << neiCell
00133                         << " is not a hex?" << abort(FatalError);
00134                 }
00135                 else
00136                 {
00137                     if (!blockedFace[oppositeFaceI])
00138                     {
00139                         blockedFace[oppositeFaceI] = true;
00140                         newFrontFaces.append(oppositeFaceI);
00141                     }
00142                 }
00143             }
00144         }
00145 
00146         frontFaces.transfer(newFrontFaces);
00147     }
00148 }
00149 
00150 
00151 // Calculate regions.
00152 void Foam::channelIndex::calcLayeredRegions
00153 (
00154     const polyMesh& mesh,
00155     const labelList& startFaces
00156 )
00157 {
00158     boolList blockedFace(mesh.nFaces(), false);
00159     walkOppositeFaces
00160     (
00161         mesh,
00162         startFaces,
00163         blockedFace
00164     );
00165 
00166 
00167     if (false)
00168     {
00169         OFstream str(mesh.time().path()/"blockedFaces.obj");
00170         label vertI = 0;
00171         forAll(blockedFace, faceI)
00172         {
00173             if (blockedFace[faceI])
00174             {
00175                 const face& f = mesh.faces()[faceI];
00176                 forAll(f, fp)
00177                 {
00178                     meshTools::writeOBJ(str, mesh.points()[f[fp]]);
00179                 }
00180                 str<< 'f';
00181                 forAll(f, fp)
00182                 {
00183                     str << ' ' << vertI+fp+1;
00184                 }
00185                 str << nl;
00186                 vertI += f.size();
00187             }
00188         }
00189     }
00190 
00191 
00192     // Do analysis for connected regions
00193     cellRegion_.reset(new regionSplit(mesh, blockedFace));
00194 
00195     Info<< "Detected " << cellRegion_().nRegions() << " layers." << nl << endl;
00196 
00197     // Sum number of entries per region
00198     regionCount_ = regionSum(scalarField(mesh.nCells(), 1.0));
00199 
00200     // Average cell centres to determine ordering.
00201     pointField regionCc
00202     (
00203         regionSum(mesh.cellCentres())
00204       / regionCount_
00205     );
00206 
00207     SortableList<scalar> sortComponent(regionCc.component(dir_));
00208 
00209     sortMap_ = sortComponent.indices();
00210 
00211     y_ = sortComponent;
00212 
00213     if (symmetric_)
00214     {
00215         y_.setSize(cellRegion_().nRegions()/2);
00216     }
00217 }
00218 
00219 
00220 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
00221 
00222 Foam::channelIndex::channelIndex
00223 (
00224     const polyMesh& mesh,
00225     const dictionary& dict
00226 )
00227 :
00228     symmetric_(readBool(dict.lookup("symmetric"))),
00229     dir_(vectorComponentsNames_.read(dict.lookup("component")))
00230 {
00231     const polyBoundaryMesh& patches = mesh.boundaryMesh();
00232 
00233     const wordList patchNames(dict.lookup("patches"));
00234 
00235     label nFaces = 0;
00236 
00237     forAll(patchNames, i)
00238     {
00239         label patchI = patches.findPatchID(patchNames[i]);
00240 
00241         if (patchI == -1)
00242         {
00243             FatalErrorIn("channelIndex::channelIndex(const polyMesh&)")
00244                 << "Illegal patch " << patchNames[i]
00245                 << ". Valid patches are " << patches.name()
00246                 << exit(FatalError);
00247         }
00248 
00249         nFaces += patches[patchI].size();
00250     }
00251 
00252     labelList startFaces(nFaces);
00253     nFaces = 0;
00254 
00255     forAll(patchNames, i)
00256     {
00257         const polyPatch& pp = patches[patches.findPatchID(patchNames[i])];
00258 
00259         forAll(pp, j)
00260         {
00261             startFaces[nFaces++] = pp.start()+j;
00262         }
00263     }
00264 
00265     // Calculate regions.
00266     calcLayeredRegions(mesh, startFaces);
00267 }
00268 
00269 
00270 Foam::channelIndex::channelIndex
00271 (
00272     const polyMesh& mesh,
00273     const labelList& startFaces,
00274     const bool symmetric,
00275     const direction dir
00276 )
00277 :
00278     symmetric_(symmetric),
00279     dir_(dir)
00280 {
00281     // Calculate regions.
00282     calcLayeredRegions(mesh, startFaces);
00283 }
00284 
00285 
00286 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
00287 
00288 
00289 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
00290 
00291 
00292 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines