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) 2016 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 "streamFunction.H"
27 #include "surfaceFields.H"
28 #include "pointFields.H"
29 #include "emptyPolyPatch.H"
30 #include "symmetryPlanePolyPatch.H"
31 #include "symmetryPolyPatch.H"
32 #include "wedgePolyPatch.H"
33 #include "addToRunTimeSelectionTable.H"
34 
35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
36 
37 namespace Foam
38 {
39 namespace functionObjects
40 {
41  defineTypeNameAndDebug(streamFunction, 0);
42 
43  addToRunTimeSelectionTable
44  (
45  functionObject,
46  streamFunction,
47  dictionary
48  );
49 }
50 }
51 
52 
53 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
54 
55 Foam::tmp<Foam::pointScalarField> Foam::functionObjects::streamFunction::calc
56 (
57  const surfaceScalarField& phi
58 ) const
59 {
60  Log << " functionObjects::" << type() << " " << name()
61  << " calculating steam-function" << endl;
62 
63  Vector<label> slabNormal((Vector<label>::one - mesh_.geometricD())/2);
64  const direction slabDir
65  (
66  slabNormal
67  & Vector<label>(Vector<label>::X, Vector<label>::Y, Vector<label>::Z)
68  );
69 
70  scalar thickness = vector(slabNormal) & mesh_.bounds().span();
71 
72  const pointMesh& pMesh = pointMesh::New(mesh_);
73 
74  tmp<pointScalarField> tstreamFunction
75  (
76  new pointScalarField
77  (
78  IOobject
79  (
80  "streamFunction",
81  time_.timeName(),
82  mesh_
83  ),
84  pMesh,
85  dimensionedScalar("zero", phi.dimensions(), 0.0)
86  )
87  );
88  pointScalarField& streamFunction = tstreamFunction.ref();
89 
90  labelList visitedPoint(mesh_.nPoints());
91  forAll(visitedPoint, pointi)
92  {
93  visitedPoint[pointi] = 0;
94  }
95  label nVisited = 0;
96  label nVisitedOld = 0;
97 
98  const faceUList& faces = mesh_.faces();
99  const pointField& points = mesh_.points();
100 
101  label nInternalFaces = mesh_.nInternalFaces();
102 
103  vectorField unitAreas(mesh_.faceAreas());
104  unitAreas /= mag(unitAreas);
105 
106  const polyPatchList& patches = mesh_.boundaryMesh();
107 
108  bool finished = true;
109 
110  // Find the boundary face with zero flux. set the stream function
111  // to zero on that face
112  bool found = false;
113 
114  do
115  {
116  found = false;
117 
118  forAll(patches, patchi)
119  {
120  const primitivePatch& bouFaces = patches[patchi];
121 
122  if (!isType<emptyPolyPatch>(patches[patchi]))
123  {
124  forAll(bouFaces, facei)
125  {
126  if
127  (
128  magSqr(phi.boundaryField()[patchi][facei]) < SMALL
129  )
130  {
131  const labelList& zeroPoints = bouFaces[facei];
132 
133  // Zero flux face found
134  found = true;
135 
136  forAll(zeroPoints, pointi)
137  {
138  if (visitedPoint[zeroPoints[pointi]] == 1)
139  {
140  found = false;
141  break;
142  }
143  }
144 
145  if (found)
146  {
147  Log << " Zero face: patch: " << patchi
148  << " face: " << facei << endl;
149 
150  forAll(zeroPoints, pointi)
151  {
152  streamFunction[zeroPoints[pointi]] = 0;
153  visitedPoint[zeroPoints[pointi]] = 1;
154  nVisited++;
155  }
156 
157  break;
158  }
159  }
160  }
161  }
162 
163  if (found) break;
164  }
165 
166  if (!found)
167  {
168  Log << " Zero flux boundary face not found. "
169  << "Using cell as a reference."
170  << endl;
171 
172  const cellList& c = mesh_.cells();
173 
174  forAll(c, cI)
175  {
176  labelList zeroPoints = c[cI].labels(mesh_.faces());
177 
178  bool 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  forAll(zeroPoints, pointi)
192  {
193  streamFunction[zeroPoints[pointi]] = 0.0;
194  visitedPoint[zeroPoints[pointi]] = 1;
195  nVisited++;
196  }
197 
198  break;
199  }
200  else
201  {
202  FatalErrorInFunction
203  << "Cannot find initialisation face or a cell."
204  << exit(FatalError);
205  }
206  }
207  }
208 
209  // Loop through all faces. If one of the points on
210  // the face has the streamfunction value different
211  // from -1, all points with -1 ont that face have the
212  // streamfunction value equal to the face flux in
213  // that point plus the value in the visited point
214  do
215  {
216  finished = true;
217 
218  for (label facei = nInternalFaces; facei<faces.size(); facei++)
219  {
220  const labelList& curBPoints = faces[facei];
221  bool bPointFound = false;
222 
223  scalar currentBStream = 0.0;
224  vector currentBStreamPoint(0, 0, 0);
225 
226  forAll(curBPoints, pointi)
227  {
228  // Check if the point has been visited
229  if (visitedPoint[curBPoints[pointi]] == 1)
230  {
231  // The point has been visited
232  currentBStream = streamFunction[curBPoints[pointi]];
233  currentBStreamPoint = points[curBPoints[pointi]];
234 
235  bPointFound = true;
236 
237  break;
238  }
239  }
240 
241  if (bPointFound)
242  {
243  // Sort out other points on the face
244  forAll(curBPoints, pointi)
245  {
246  // Check if the point has been visited
247  if (visitedPoint[curBPoints[pointi]] == 0)
248  {
249  label patchNo =
250  mesh_.boundaryMesh().whichPatch(facei);
251 
252  if
253  (
254  !isType<emptyPolyPatch>(patches[patchNo])
255  && !isType<symmetryPlanePolyPatch>
256  (patches[patchNo])
257  && !isType<symmetryPolyPatch>(patches[patchNo])
258  && !isType<wedgePolyPatch>(patches[patchNo])
259  )
260  {
261  label faceNo =
262  mesh_.boundaryMesh()[patchNo]
263  .whichFace(facei);
264 
265  vector edgeHat =
266  points[curBPoints[pointi]]
267  - currentBStreamPoint;
268  edgeHat.replace(slabDir, 0);
269  edgeHat /= mag(edgeHat);
270 
271  vector nHat = unitAreas[facei];
272 
273  if (edgeHat.y() > VSMALL)
274  {
275  visitedPoint[curBPoints[pointi]] = 1;
276  nVisited++;
277 
278  streamFunction[curBPoints[pointi]] =
279  currentBStream
280  + phi.boundaryField()[patchNo][faceNo]
281  *sign(nHat.x());
282  }
283  else if (edgeHat.y() < -VSMALL)
284  {
285  visitedPoint[curBPoints[pointi]] = 1;
286  nVisited++;
287 
288  streamFunction[curBPoints[pointi]] =
289  currentBStream
290  - phi.boundaryField()[patchNo][faceNo]
291  *sign(nHat.x());
292  }
293  else
294  {
295  if (edgeHat.x() > VSMALL)
296  {
297  visitedPoint[curBPoints[pointi]] = 1;
298  nVisited++;
299 
300  streamFunction[curBPoints[pointi]] =
301  currentBStream
302  + phi.boundaryField()[patchNo][faceNo]
303  *sign(nHat.y());
304  }
305  else if (edgeHat.x() < -VSMALL)
306  {
307  visitedPoint[curBPoints[pointi]] = 1;
308  nVisited++;
309 
310  streamFunction[curBPoints[pointi]] =
311  currentBStream
312  - phi.boundaryField()[patchNo][faceNo]
313  *sign(nHat.y());
314  }
315  }
316  }
317  }
318  }
319  }
320  else
321  {
322  finished = false;
323  }
324  }
325 
326  for (label facei=0; facei<nInternalFaces; facei++)
327  {
328  // Get the list of point labels for the face
329  const labelList& curPoints = faces[facei];
330 
331  bool pointFound = false;
332 
333  scalar currentStream = 0.0;
334  point currentStreamPoint(0, 0, 0);
335 
336  forAll(curPoints, pointi)
337  {
338  // Check if the point has been visited
339  if (visitedPoint[curPoints[pointi]] == 1)
340  {
341  // The point has been visited
342  currentStream = streamFunction[curPoints[pointi]];
343  currentStreamPoint = points[curPoints[pointi]];
344  pointFound = true;
345 
346  break;
347  }
348  }
349 
350  if (pointFound)
351  {
352  // Sort out other points on the face
353  forAll(curPoints, pointi)
354  {
355  // Check if the point has been visited
356  if (visitedPoint[curPoints[pointi]] == 0)
357  {
358  vector edgeHat =
359  points[curPoints[pointi]] - currentStreamPoint;
360 
361  edgeHat.replace(slabDir, 0);
362  edgeHat /= mag(edgeHat);
363 
364  vector nHat = unitAreas[facei];
365 
366  if (edgeHat.y() > VSMALL)
367  {
368  visitedPoint[curPoints[pointi]] = 1;
369  nVisited++;
370 
371  streamFunction[curPoints[pointi]] =
372  currentStream
373  + phi[facei]*sign(nHat.x());
374  }
375  else if (edgeHat.y() < -VSMALL)
376  {
377  visitedPoint[curPoints[pointi]] = 1;
378  nVisited++;
379 
380  streamFunction[curPoints[pointi]] =
381  currentStream
382  - phi[facei]*sign(nHat.x());
383  }
384  }
385  }
386  }
387  else
388  {
389  finished = false;
390  }
391  }
392 
393  if (nVisited == nVisitedOld)
394  {
395  // Find new seed. This must be a
396  // multiply connected domain
397  Log << " Exhausted a seed, looking for new seed "
398  << "(this is correct for multiply connected domains).";
399 
400  break;
401  }
402  else
403  {
404  nVisitedOld = nVisited;
405  }
406  } while (!finished);
407  } while (!finished);
408 
409  // Normalise the stream-function by the 2D mesh thickness
410  streamFunction /= thickness;
411  streamFunction.boundaryFieldRef() = 0.0;
412 
413  return tstreamFunction;
414 }
415 
416 
417 bool Foam::functionObjects::streamFunction::calc()
418 {
419  if (foundObject<surfaceScalarField>(fieldName_))
420  {
421  const surfaceScalarField& phi =
422  mesh_.lookupObject<surfaceScalarField>(fieldName_);
423 
424  return store(resultName_, calc(phi));
425  }
426  else
427  {
428  return false;
429  }
430 }
431 
432 
433 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
434 
435 Foam::functionObjects::streamFunction::streamFunction
436 (
437  const word& name,
438  const Time& runTime,
439  const dictionary& dict
440 )
441 :
442  fieldExpression(name, runTime, dict, "phi")
443 {
444  setResultName("streamFunction", "phi");
445 
446  label nD = mesh_.nGeometricD();
447 
448  if (nD != 2)
449  {
450  FatalErrorInFunction
451  << "Case is not 2D, stream-function cannot be computed"
452  << exit(FatalError);
453  }
454 }
455 
456 
457 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
458 
459 Foam::functionObjects::streamFunction::~streamFunction()
460 {}
461 
462 
463 // ************************************************************************* //
Foam::pointScalarField
GeometricField< scalar, pointPatchField, pointMesh > pointScalarField
Definition: pointFields.H:49
surfaceFields.H
Foam::surfaceFields.
Foam::polyPatchList
PtrList< polyPatch > polyPatchList
container classes for polyPatch
Definition: polyPatchList.H:45
Foam::sign
dimensionedScalar sign(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:179
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
Foam::functionObjects::streamFunction::streamFunction
streamFunction(const word &name, const Time &runTime, const dictionary &dict)
Construct for given objectRegistry and dictionary.
Definition: streamFunction.C:436
Foam::label
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
Foam::exit
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
Foam::FatalError
error FatalError
Foam::dictionary
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
Foam::direction
uint8_t direction
Definition: direction.H:45
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
Foam::vector
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
Foam::VectorSpace::replace
void replace(const direction, const Cmpt &)
Definition: VectorSpaceI.H:159
Foam::Time
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
Definition: Time.H:68
addToRunTimeSelectionTable.H
Macros for easy insertion into run-time selection tables.
Foam::MeshObject< polyMesh, UpdateableMeshObject, pointMesh >::New
static const pointMesh & New(const polyMesh &mesh)
Definition: MeshObject.C:44
Foam::pointField
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
Foam::faceUList
UList< face > faceUList
Definition: faceListFwd.H:41
Foam::word
A class for handling words, derived from string.
Definition: word.H:59
Foam::labelList
List< label > labelList
A List of labels.
Definition: labelList.H:56
Foam::magSqr
dimensioned< scalar > magSqr(const dimensioned< Type > &)
Foam::primitivePatch
PrimitivePatch< face, SubList, const pointField & > primitivePatch
Addressing for a faceList slice.
Definition: primitivePatch.H:45
Foam::name
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
Foam::type
fileName::Type type(const fileName &, const bool followLink=true)
Return the file type: DIRECTORY or FILE.
Definition: POSIX.C:485
Foam::GeometricField::ref
Internal & ref()
Return a reference to the dimensioned internal field.
Definition: GeometricField.C:743
patchi
label patchi
Definition: getPatchFieldScalar.H:1
Foam::point
vector point
Point is a vector.
Definition: point.H:41
Foam::dimensionedScalar
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
Definition: dimensionedScalarFwd.H:41
Log
#define Log
Report write to Foam::Info if the local log switch is true.
Definition: messageStream.H:284
Foam::mag
dimensioned< scalar > mag(const dimensioned< Type > &)
Foam::vectorField
Field< vector > vectorField
Specialisation of Field<T> for vector.
Definition: primitiveFieldsFwd.H:49
Foam::functionObjects::defineTypeNameAndDebug
defineTypeNameAndDebug(fvMeshFunctionObject, 0)
Foam::tmp
A class for managing temporary objects.
Definition: PtrList.H:53
Foam::surfaceScalarField
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
Definition: surfaceFieldsFwd.H:52
Foam::cellList
List< cell > cellList
list of cells
Definition: cellList.H:42
Foam::functionObjects::addToRunTimeSelectionTable
addToRunTimeSelectionTable(functionObject, add, dictionary)
Foam
Namespace for OpenFOAM.
Definition: combustionModel.C:30
Generated by   doxygen 1.8.13