FreeFOAM The Cross-Platform CFD Toolkit
streamFunction.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 Application
25  streamFunction
26 
27 Description
28  Calculates and writes the stream function of velocity field U at each time
29 
30 Usage
31 
32  - streamFunction [OPTIONS]
33 
34  @param -noZero \n
35  Ignore timestep 0.
36 
37  @param -constant \n
38  Include the constant directory.
39 
40  @param -time <time>\n
41  Apply only to specific time.
42 
43  @param -latestTime \n
44  Only apply to latest time step.
45 
46  @param -case <dir>\n
47  Case directory.
48 
49  @param -parallel \n
50  Run in parallel.
51 
52  @param -help \n
53  Display help message.
54 
55  @param -doc \n
56  Display Doxygen API documentation page for this application.
57 
58  @param -srcDoc \n
59  Display Doxygen source documentation page for this application.
60 
61 \*---------------------------------------------------------------------------*/
62 
63 #include <finiteVolume/fvCFD.H>
64 #include <OpenFOAM/pointFields.H>
68 #include <OpenFOAM/OSspecific.H>
69 
70 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
71 // Main program:
72 
73 int main(int argc, char *argv[])
74 {
75  timeSelector::addOptions();
76 
77 # include <OpenFOAM/setRootCase.H>
78 # include <OpenFOAM/createTime.H>
79 
80  instantList timeDirs = timeSelector::select0(runTime, args);
81 
83 
84  pointMesh pMesh(mesh);
85 
86  forAll(timeDirs, timeI)
87  {
88  runTime.setTime(timeDirs[timeI], timeI);
89 
90  Info<< nl << "Time: " << runTime.timeName() << endl;
91 
92  IOobject phiHeader
93  (
94  "phi",
95  runTime.timeName(),
96  mesh,
97  IOobject::NO_READ
98  );
99 
100  if (phiHeader.headerOk())
101  {
102  mesh.readUpdate();
103 
104  Info<< nl << "Reading field phi" << endl;
105 
107  (
108  IOobject
109  (
110  "phi",
111  runTime.timeName(),
112  mesh,
113  IOobject::MUST_READ,
114  IOobject::NO_WRITE
115  ),
116  mesh
117  );
118 
119  pointScalarField streamFunction
120  (
121  IOobject
122  (
123  "streamFunction",
124  runTime.timeName(),
125  mesh,
126  IOobject::NO_READ,
127  IOobject::NO_WRITE
128  ),
129  pMesh,
130  dimensionedScalar("zero", phi.dimensions(), 0.0)
131  );
132 
133  labelList visitedPoint(mesh.nPoints());
134  forAll (visitedPoint, pointI)
135  {
136  visitedPoint[pointI] = 0;
137  }
138  label nVisited = 0;
139  label nVisitedOld = 0;
140 
141  const unallocFaceList& faces = mesh.faces();
142  const pointField& points = mesh.points();
143 
144  label nInternalFaces = mesh.nInternalFaces();
145 
146  vectorField unitAreas = mesh.faceAreas();
147  unitAreas /= mag(unitAreas);
148 
150 
151  bool finished = true;
152 
153  // Find the boundary face with zero flux. set the stream function
154  // to zero on that face
155  bool found = false;
156 
157  do
158  {
159  found = false;
160 
161  forAll (patches, patchI)
162  {
163  const primitivePatch& bouFaces = patches[patchI];
164 
165  if (!isType<emptyPolyPatch>(patches[patchI]))
166  {
167  forAll (bouFaces, faceI)
168  {
169  if
170  (
171  magSqr(phi.boundaryField()[patchI][faceI])
172  < SMALL
173  )
174  {
175  const labelList& zeroPoints = bouFaces[faceI];
176 
177  // Zero flux face found
178  found = true;
179 
180  forAll (zeroPoints, pointI)
181  {
182  if (visitedPoint[zeroPoints[pointI]] == 1)
183  {
184  found = false;
185  break;
186  }
187  }
188 
189  if (found)
190  {
191  Info<< "Zero face: patch: " << patchI
192  << " face: " << faceI << endl;
193 
194  forAll (zeroPoints, pointI)
195  {
196  streamFunction[zeroPoints[pointI]] = 0;
197  visitedPoint[zeroPoints[pointI]] = 1;
198  nVisited++;
199  }
200 
201  break;
202  }
203  }
204  }
205  }
206 
207  if (found) break;
208  }
209 
210  if (!found)
211  {
212  Info<< "zero flux boundary face not found. "
213  << "Using cell as a reference."
214  << endl;
215 
216  const cellList& c = mesh.cells();
217 
218  forAll (c, cI)
219  {
220  labelList zeroPoints = c[cI].labels(mesh.faces());
221 
222  bool found = true;
223 
224  forAll (zeroPoints, pointI)
225  {
226  if (visitedPoint[zeroPoints[pointI]] == 1)
227  {
228  found = false;
229  break;
230  }
231  }
232 
233  if (found)
234  {
235  forAll (zeroPoints, pointI)
236  {
237  streamFunction[zeroPoints[pointI]] = 0.0;
238  visitedPoint[zeroPoints[pointI]] = 1;
239  nVisited++;
240  }
241 
242  break;
243  }
244  else
245  {
247  << "Cannot find initialisation face or a cell."
248  << abort(FatalError);
249  }
250  }
251  }
252 
253  // Loop through all faces. If one of the points on
254  // the face has the streamfunction value different
255  // from -1, all points with -1 ont that face have the
256  // streamfunction value equal to the face flux in
257  // that point plus the value in the visited point
258  do
259  {
260  finished = true;
261 
262  for
263  (
264  label faceI = nInternalFaces;
265  faceI<faces.size();
266  faceI++
267  )
268  {
269  const labelList& curBPoints = faces[faceI];
270  bool bPointFound = false;
271 
272  scalar currentBStream = 0.0;
273  vector currentBStreamPoint(0, 0, 0);
274 
275  forAll (curBPoints, pointI)
276  {
277  // Check if the point has been visited
278  if (visitedPoint[curBPoints[pointI]] == 1)
279  {
280  // The point has been visited
281  currentBStream =
282  streamFunction[curBPoints[pointI]];
283  currentBStreamPoint =
284  points[curBPoints[pointI]];
285 
286  bPointFound = true;
287 
288  break;
289  }
290  }
291 
292  if (bPointFound)
293  {
294  // Sort out other points on the face
295  forAll (curBPoints, pointI)
296  {
297  // Check if the point has been visited
298  if (visitedPoint[curBPoints[pointI]] == 0)
299  {
300  label patchNo =
301  mesh.boundaryMesh().whichPatch(faceI);
302 
303  if
304  (
305  !isType<emptyPolyPatch>
306  (patches[patchNo])
307  && !isType<symmetryPolyPatch>
308  (patches[patchNo])
309  && !isType<wedgePolyPatch>
310  (patches[patchNo])
311  )
312  {
313  label faceNo =
314  mesh.boundaryMesh()[patchNo]
315  .whichFace(faceI);
316 
317  vector edgeHat =
318  points[curBPoints[pointI]]
319  - currentBStreamPoint;
320  edgeHat.replace(vector::Z, 0);
321  edgeHat /= mag(edgeHat);
322 
323  vector nHat = unitAreas[faceI];
324 
325  if (edgeHat.y() > VSMALL)
326  {
327  visitedPoint[curBPoints[pointI]] =
328  1;
329  nVisited++;
330 
331  streamFunction[curBPoints[pointI]]
332  =
333  currentBStream
334  + phi.boundaryField()
335  [patchNo][faceNo]
336  *sign(nHat.x());
337  }
338  else if (edgeHat.y() < -VSMALL)
339  {
340  visitedPoint[curBPoints[pointI]] =
341  1;
342  nVisited++;
343 
344  streamFunction[curBPoints[pointI]]
345  =
346  currentBStream
347  - phi.boundaryField()
348  [patchNo][faceNo]
349  *sign(nHat.x());
350  }
351  else
352  {
353  if (edgeHat.x() > VSMALL)
354  {
355  visitedPoint
356  [curBPoints[pointI]] = 1;
357  nVisited++;
358 
359  streamFunction
360  [curBPoints[pointI]] =
361  currentBStream
362  + phi.boundaryField()
363  [patchNo][faceNo]
364  *sign(nHat.y());
365  }
366  else if (edgeHat.x() < -VSMALL)
367  {
368  visitedPoint
369  [curBPoints[pointI]] = 1;
370  nVisited++;
371 
372  streamFunction
373  [curBPoints[pointI]] =
374  currentBStream
375  - phi.boundaryField()
376  [patchNo][faceNo]
377  *sign(nHat.y());
378  }
379  }
380  }
381  }
382  }
383  }
384  else
385  {
386  finished = false;
387  }
388  }
389 
390  for (label faceI=0; faceI<nInternalFaces; faceI++)
391  {
392  // Get the list of point labels for the face
393  const labelList& curPoints = faces[faceI];
394 
395  bool pointFound = false;
396 
397  scalar currentStream = 0.0;
398  point currentStreamPoint(0, 0, 0);
399 
400  forAll (curPoints, pointI)
401  {
402  // Check if the point has been visited
403  if (visitedPoint[curPoints[pointI]] == 1)
404  {
405  // The point has been visited
406  currentStream =
407  streamFunction[curPoints[pointI]];
408  currentStreamPoint =
409  points[curPoints[pointI]];
410  pointFound = true;
411 
412  break;
413  }
414  }
415 
416  if (pointFound)
417  {
418  // Sort out other points on the face
419  forAll (curPoints, pointI)
420  {
421  // Check if the point has been visited
422  if (visitedPoint[curPoints[pointI]] == 0)
423  {
424  vector edgeHat =
425  points[curPoints[pointI]]
426  - currentStreamPoint;
427 
428  edgeHat.replace(vector::Z, 0);
429  edgeHat /= mag(edgeHat);
430 
431  vector nHat = unitAreas[faceI];
432 
433  if (edgeHat.y() > VSMALL)
434  {
435  visitedPoint[curPoints[pointI]] = 1;
436  nVisited++;
437 
438  streamFunction[curPoints[pointI]] =
439  currentStream
440  + phi[faceI]*sign(nHat.x());
441  }
442  else if (edgeHat.y() < -VSMALL)
443  {
444  visitedPoint[curPoints[pointI]] = 1;
445  nVisited++;
446 
447  streamFunction[curPoints[pointI]] =
448  currentStream
449  - phi[faceI]*sign(nHat.x());
450  }
451  }
452  }
453  }
454  else
455  {
456  finished = false;
457  }
458  }
459 
460  Info<< ".";
461 // Info<< "One pass, n visited = " << nVisited << endl;
462 
463  if (nVisited == nVisitedOld)
464  {
465  // Find new seed. This must be a
466  // multiply connected domain
467  Info<< nl << "Exhausted a seed. Looking for new seed "
468  << "(this is correct for multiply connected "
469  << "domains).";
470 
471  break;
472  }
473  else
474  {
475  nVisitedOld = nVisited;
476  }
477  } while (!finished);
478 
479  Info << endl;
480  } while (!finished);
481 
482  streamFunction.boundaryField() = 0.0;
483  streamFunction.write();
484  }
485  else
486  {
488  << "Flux field does not exist."
489  << " Stream function not calculated" << endl;
490  }
491  }
492 
493  Info<< "\nEnd\n" << endl;
494 
495  return 0;
496 }
497 
498 
499 // ************************ vim: set sw=4 sts=4 et: ************************ //