streamFunction.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2016-2020 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  pointScalarField::New
77  (
78  "streamFunction",
79  pMesh,
80  dimensionedScalar(phi.dimensions(), 0)
81  )
82  );
83  pointScalarField& streamFunction = tstreamFunction.ref();
84 
85  labelList visitedPoint(mesh_.nPoints());
86  forAll(visitedPoint, pointi)
87  {
88  visitedPoint[pointi] = 0;
89  }
90  label nVisited = 0;
91  label nVisitedOld = 0;
92 
93  const faceUList& faces = mesh_.faces();
94  const pointField& points = mesh_.points();
95 
96  label nInternalFaces = mesh_.nInternalFaces();
97 
98  vectorField unitAreas(mesh_.faceAreas());
99  unitAreas /= mag(unitAreas);
100 
101  const polyPatchList& patches = mesh_.boundaryMesh();
102 
103  bool finished = true;
104 
105  // Find the boundary face with zero flux. set the stream function
106  // to zero on that face
107  bool found = false;
108 
109  do
110  {
111  found = false;
112 
113  forAll(patches, patchi)
114  {
115  const primitivePatch& bouFaces = patches[patchi];
116 
117  if (!isType<emptyPolyPatch>(patches[patchi]))
118  {
119  forAll(bouFaces, facei)
120  {
121  if
122  (
123  magSqr(phi.boundaryField()[patchi][facei]) < small
124  )
125  {
126  const labelList& zeroPoints = bouFaces[facei];
127 
128  // Zero flux face found
129  found = true;
130 
131  forAll(zeroPoints, pointi)
132  {
133  if (visitedPoint[zeroPoints[pointi]] == 1)
134  {
135  found = false;
136  break;
137  }
138  }
139 
140  if (found)
141  {
142  Log << " Zero face: patch: " << patchi
143  << " face: " << facei << endl;
144 
145  forAll(zeroPoints, pointi)
146  {
147  streamFunction[zeroPoints[pointi]] = 0;
148  visitedPoint[zeroPoints[pointi]] = 1;
149  nVisited++;
150  }
151 
152  break;
153  }
154  }
155  }
156  }
157 
158  if (found) break;
159  }
160 
161  if (!found)
162  {
163  Log << " Zero flux boundary face not found. "
164  << "Using cell as a reference."
165  << endl;
166 
167  const cellList& c = mesh_.cells();
168 
169  forAll(c, cI)
170  {
171  labelList zeroPoints = c[cI].labels(mesh_.faces());
172 
173  bool found = true;
174 
175  forAll(zeroPoints, pointi)
176  {
177  if (visitedPoint[zeroPoints[pointi]] == 1)
178  {
179  found = false;
180  break;
181  }
182  }
183 
184  if (found)
185  {
186  forAll(zeroPoints, pointi)
187  {
188  streamFunction[zeroPoints[pointi]] = 0.0;
189  visitedPoint[zeroPoints[pointi]] = 1;
190  nVisited++;
191  }
192 
193  break;
194  }
195  else
196  {
197  FatalErrorInFunction
198  << "Cannot find initialisation face or a cell."
199  << exit(FatalError);
200  }
201  }
202  }
203 
204  // Loop through all faces. If one of the points on
205  // the face has the streamfunction value different
206  // from -1, all points with -1 ont that face have the
207  // streamfunction value equal to the face flux in
208  // that point plus the value in the visited point
209  do
210  {
211  finished = true;
212 
213  for (label facei = nInternalFaces; facei<faces.size(); facei++)
214  {
215  const labelList& curBPoints = faces[facei];
216  bool bPointFound = false;
217 
218  scalar currentBStream = 0.0;
219  vector currentBStreamPoint(0, 0, 0);
220 
221  forAll(curBPoints, pointi)
222  {
223  // Check if the point has been visited
224  if (visitedPoint[curBPoints[pointi]] == 1)
225  {
226  // The point has been visited
227  currentBStream = streamFunction[curBPoints[pointi]];
228  currentBStreamPoint = points[curBPoints[pointi]];
229 
230  bPointFound = true;
231 
232  break;
233  }
234  }
235 
236  if (bPointFound)
237  {
238  // Sort out other points on the face
239  forAll(curBPoints, pointi)
240  {
241  // Check if the point has been visited
242  if (visitedPoint[curBPoints[pointi]] == 0)
243  {
244  label patchNo =
245  mesh_.boundaryMesh().whichPatch(facei);
246 
247  if
248  (
249  !isType<emptyPolyPatch>(patches[patchNo])
250  && !isType<symmetryPlanePolyPatch>
251  (patches[patchNo])
252  && !isType<symmetryPolyPatch>(patches[patchNo])
253  && !isType<wedgePolyPatch>(patches[patchNo])
254  )
255  {
256  label faceNo =
257  mesh_.boundaryMesh()[patchNo]
258  .whichFace(facei);
259 
260  vector edgeHat =
261  points[curBPoints[pointi]]
262  - currentBStreamPoint;
263  edgeHat.replace(slabDir, 0);
264  edgeHat /= mag(edgeHat);
265 
266  vector nHat = unitAreas[facei];
267 
268  if (edgeHat.y() > vSmall)
269  {
270  visitedPoint[curBPoints[pointi]] = 1;
271  nVisited++;
272 
273  streamFunction[curBPoints[pointi]] =
274  currentBStream
275  + phi.boundaryField()[patchNo][faceNo]
276  *sign(nHat.x());
277  }
278  else if (edgeHat.y() < -vSmall)
279  {
280  visitedPoint[curBPoints[pointi]] = 1;
281  nVisited++;
282 
283  streamFunction[curBPoints[pointi]] =
284  currentBStream
285  - phi.boundaryField()[patchNo][faceNo]
286  *sign(nHat.x());
287  }
288  else
289  {
290  if (edgeHat.x() > vSmall)
291  {
292  visitedPoint[curBPoints[pointi]] = 1;
293  nVisited++;
294 
295  streamFunction[curBPoints[pointi]] =
296  currentBStream
297  + phi.boundaryField()[patchNo][faceNo]
298  *sign(nHat.y());
299  }
300  else if (edgeHat.x() < -vSmall)
301  {
302  visitedPoint[curBPoints[pointi]] = 1;
303  nVisited++;
304 
305  streamFunction[curBPoints[pointi]] =
306  currentBStream
307  - phi.boundaryField()[patchNo][faceNo]
308  *sign(nHat.y());
309  }
310  }
311  }
312  }
313  }
314  }
315  else
316  {
317  finished = false;
318  }
319  }
320 
321  for (label facei=0; facei<nInternalFaces; facei++)
322  {
323  // Get the list of point labels for the face
324  const labelList& curPoints = faces[facei];
325 
326  bool pointFound = false;
327 
328  scalar currentStream = 0.0;
329  point currentStreamPoint(0, 0, 0);
330 
331  forAll(curPoints, pointi)
332  {
333  // Check if the point has been visited
334  if (visitedPoint[curPoints[pointi]] == 1)
335  {
336  // The point has been visited
337  currentStream = streamFunction[curPoints[pointi]];
338  currentStreamPoint = points[curPoints[pointi]];
339  pointFound = true;
340 
341  break;
342  }
343  }
344 
345  if (pointFound)
346  {
347  // Sort out other points on the face
348  forAll(curPoints, pointi)
349  {
350  // Check if the point has been visited
351  if (visitedPoint[curPoints[pointi]] == 0)
352  {
353  vector edgeHat =
354  points[curPoints[pointi]] - currentStreamPoint;
355 
356  edgeHat.replace(slabDir, 0);
357  edgeHat /= mag(edgeHat);
358 
359  vector nHat = unitAreas[facei];
360 
361  if (edgeHat.y() > vSmall)
362  {
363  visitedPoint[curPoints[pointi]] = 1;
364  nVisited++;
365 
366  streamFunction[curPoints[pointi]] =
367  currentStream
368  + phi[facei]*sign(nHat.x());
369  }
370  else if (edgeHat.y() < -vSmall)
371  {
372  visitedPoint[curPoints[pointi]] = 1;
373  nVisited++;
374 
375  streamFunction[curPoints[pointi]] =
376  currentStream
377  - phi[facei]*sign(nHat.x());
378  }
379  }
380  }
381  }
382  else
383  {
384  finished = false;
385  }
386  }
387 
388  if (nVisited == nVisitedOld)
389  {
390  // Find new seed. This must be a
391  // multiply connected domain
392  Log << " Exhausted a seed, looking for new seed "
393  << "(this is correct for multiply connected domains).";
394 
395  break;
396  }
397  else
398  {
399  nVisitedOld = nVisited;
400  }
401  } while (!finished);
402  } while (!finished);
403 
404  // Normalise the stream-function by the 2D mesh thickness
405  streamFunction /= thickness;
406  streamFunction.boundaryFieldRef() = 0.0;
407 
408  return tstreamFunction;
409 }
410 
411 
412 bool Foam::functionObjects::streamFunction::calc()
413 {
414  if (foundObject<surfaceScalarField>(fieldName_))
415  {
416  const surfaceScalarField& phi =
417  mesh_.lookupObject<surfaceScalarField>(fieldName_);
418 
419  return store(resultName_, calc(phi));
420  }
421  else
422  {
423  cannotFindObject<surfaceScalarField>(fieldName_);
424 
425  return false;
426  }
427 }
428 
429 
430 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
431 
432 Foam::functionObjects::streamFunction::streamFunction
433 (
434  const word& name,
435  const Time& runTime,
436  const dictionary& dict
437 )
438 :
439  fieldExpression(name, runTime, dict, typeName, "phi")
440 {
441  label nD = mesh_.nGeometricD();
442 
443  if (nD != 2)
444  {
445  FatalErrorInFunction
446  << "Case is not 2D, stream-function cannot be computed"
447  << exit(FatalError);
448  }
449 }
450 
451 
452 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
453 
454 Foam::functionObjects::streamFunction::~streamFunction()
455 {}
456 
457 
458 // ************************************************************************* //
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:434
Foam::functionObjects::streamFunction::streamFunction
streamFunction(const word &name, const Time &runTime, const dictionary &dict)
Construct for given objectRegistry and dictionary.
Definition: streamFunction.C:433
Foam::label
FvWallInfoData< WallInfo, label > label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: FvWallInfoData.H:189
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:156
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
Foam::direction
uint8_t direction
Definition: direction.H:45
Foam::functionObjects::addToRunTimeSelectionTable
addToRunTimeSelectionTable(functionObject, Qdot, dictionary)
Foam::GeometricField< scalar, pointPatchField, pointMesh >::New
static tmp< GeometricField< scalar, pointPatchField, pointMesh > > New(const word &name, const Internal &, const PtrList< pointPatchField< scalar >> &)
Return a temporary field constructed from name,.
Definition: GeometricField.C:766
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
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:149
Foam::Time
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
Definition: Time.H:69
addToRunTimeSelectionTable.H
Macros for easy insertion into run-time selection tables.
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::pointScalarField
GeometricField< scalar, pointPatchField, pointMesh > pointScalarField
Definition: pointFieldsFwd.H:58
Foam::magSqr
dimensioned< scalar > magSqr(const dimensioned< Type > &)
Foam::primitivePatch
PrimitivePatch< SubList< face >, const pointField & > primitivePatch
Addressing for a faceList slice.
Definition: primitivePatch.H:44
Foam::name
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
Foam::functionObjects::defineTypeNameAndDebug
defineTypeNameAndDebug(Qdot, 0)
Foam::GeometricField::ref
Internal & ref()
Return a reference to the dimensioned internal field.
Definition: GeometricField.C:1015
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
Foam::type
fileType type(const fileName &, const bool checkVariants=true, const bool followLink=true)
Return the file type: directory or file.
Definition: POSIX.C:488
Log
#define Log
Report write to Foam::Info if the local log switch is true.
Definition: messageStream.H:280
Foam::mag
dimensioned< scalar > mag(const dimensioned< Type > &)
Foam::vectorField
Field< vector > vectorField
Specialisation of Field<T> for vector.
Definition: primitiveFieldsFwd.H:49
Foam::tmp
A class for managing temporary objects.
Definition: PtrList.H:53
Foam::surfaceScalarField
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
Definition: surfaceFieldsFwd.H:58
Foam::cellList
List< cell > cellList
list of cells
Definition: cellList.H:42
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.H:213