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