Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
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
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
00049
00050
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
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
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
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
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
00193 cellRegion_.reset(new regionSplit(mesh, blockedFace));
00194
00195 Info<< "Detected " << cellRegion_().nRegions() << " layers." << nl << endl;
00196
00197
00198 regionCount_ = regionSum(scalarField(mesh.nCells(), 1.0));
00199
00200
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
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
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
00282 calcLayeredRegions(mesh, startFaces);
00283 }
00284
00285
00286
00287
00288
00289
00290
00291
00292