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 | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2016 OpenFOAM Foundation
9 Copyright (C) 2019-2020 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27\*---------------------------------------------------------------------------*/
28
29#include "streamFunction.H"
30#include "surfaceFields.H"
31#include "pointFields.H"
32#include "emptyPolyPatch.H"
34#include "symmetryPolyPatch.H"
35#include "wedgePolyPatch.H"
37
38// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
39
40namespace Foam
41{
42namespace functionObjects
43{
46}
47}
48
49
50// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
51
52Foam::tmp<Foam::pointScalarField> Foam::functionObjects::streamFunction::calc
53(
54 const surfaceScalarField& phi
55) const
56{
57 Log << " functionObjects::" << type() << " " << name()
58 << " calculating stream-function" << endl;
59
60 Vector<label> slabNormal((Vector<label>::one - mesh_.geometricD())/2);
61 const direction slabDir
62 (
63 slabNormal
64 & Vector<label>(Vector<label>::X, Vector<label>::Y, Vector<label>::Z)
65 );
66
67 scalar thickness = vector(slabNormal) & mesh_.bounds().span();
68
69 const pointMesh& pMesh = pointMesh::New(mesh_);
70
71 auto tstreamFunction = tmp<pointScalarField>::New
72 (
73 IOobject
74 (
75 "streamFunction",
77 mesh_
78 ),
79 pMesh,
81 );
82 pointScalarField& streamFunction = tstreamFunction.ref();
83
84 labelList visitedPoint(mesh_.nPoints(), Zero);
85
86 label nVisited = 0;
87 label nVisitedOld = 0;
88
89 const faceUList& faces = mesh_.faces();
90 const pointField& points = mesh_.points();
91
92 label nInternalFaces = mesh_.nInternalFaces();
93
94 vectorField unitAreas(mesh_.faceAreas());
95 unitAreas /= mag(unitAreas);
96
98
99 bool finished = true;
100
101 // Find the boundary face with zero flux. Set the stream function
102 // to zero on that face
103 bool found = false;
104
105 do
106 {
107 found = false;
108
109 forAll(patches, patchi)
110 {
111 const primitivePatch& bouFaces = patches[patchi];
112
113 if (!isType<emptyPolyPatch>(patches[patchi]))
114 {
115 forAll(bouFaces, facei)
116 {
117 if (magSqr(phi.boundaryField()[patchi][facei]) < SMALL)
118 {
119 const labelList& zeroPoints = bouFaces[facei];
120
121 // Zero flux face found
122 found = true;
123
124 forAll(zeroPoints, pointi)
125 {
126 if (visitedPoint[zeroPoints[pointi]] == 1)
127 {
128 found = false;
129 break;
130 }
131 }
132
133 if (found)
134 {
135 Log << " Zero face: patch: " << patchi
136 << " face: " << facei << endl;
137
138 forAll(zeroPoints, pointi)
139 {
140 streamFunction[zeroPoints[pointi]] = 0;
141 visitedPoint[zeroPoints[pointi]] = 1;
142 nVisited++;
143 }
144
145 break;
146 }
147 }
148 }
149 }
150
151 if (found) break;
152 }
153
154 if (!found)
155 {
156 Log << " Zero flux boundary face not found. "
157 << "Using cell as a reference."
158 << endl;
159
160 const cellList& c = mesh_.cells();
161
162 forAll(c, ci)
163 {
164 labelList zeroPoints = c[ci].labels(mesh_.faces());
165
166 bool found = true;
167
168 forAll(zeroPoints, pointi)
169 {
170 if (visitedPoint[zeroPoints[pointi]] == 1)
171 {
172 found = false;
173 break;
174 }
175 }
176
177 if (found)
178 {
179 forAll(zeroPoints, pointi)
180 {
181 streamFunction[zeroPoints[pointi]] = 0.0;
182 visitedPoint[zeroPoints[pointi]] = 1;
183 nVisited++;
184 }
185
186 break;
187 }
188 else
189 {
191 << "Cannot find initialisation face or a cell."
192 << exit(FatalError);
193 }
194 }
195 }
196
197 // Loop through all faces. If one of the points on
198 // the face has the streamfunction value different
199 // from -1, all points with -1 ont that face have the
200 // streamfunction value equal to the face flux in
201 // that point plus the value in the visited point
202 do
203 {
204 finished = true;
205
206 for (label facei = nInternalFaces; facei<faces.size(); facei++)
207 {
208 const labelList& curBPoints = faces[facei];
209 bool bPointFound = false;
210
211 scalar currentBStream = 0.0;
212 vector currentBStreamPoint(0, 0, 0);
213
214 forAll(curBPoints, pointi)
215 {
216 // Check if the point has been visited
217 if (visitedPoint[curBPoints[pointi]] == 1)
218 {
219 // The point has been visited
220 currentBStream = streamFunction[curBPoints[pointi]];
221 currentBStreamPoint = points[curBPoints[pointi]];
222
223 bPointFound = true;
224
225 break;
226 }
227 }
228
229 if (bPointFound)
230 {
231 // Sort out other points on the face
232 forAll(curBPoints, pointi)
233 {
234 // Check if the point has been visited
235 if (visitedPoint[curBPoints[pointi]] == 0)
236 {
237 label patchNo =
239
240 if
241 (
242 !isType<emptyPolyPatch>(patches[patchNo])
243 && !isType<symmetryPlanePolyPatch>
244 (patches[patchNo])
245 && !isType<symmetryPolyPatch>(patches[patchNo])
246 && !isType<wedgePolyPatch>(patches[patchNo])
247 )
248 {
249 label faceNo =
250 mesh_.boundaryMesh()[patchNo]
251 .whichFace(facei);
252
253 vector edgeHat =
254 points[curBPoints[pointi]]
255 - currentBStreamPoint;
256 edgeHat.replace(slabDir, 0);
257 edgeHat.normalise();
258
259 vector nHat = unitAreas[facei];
260
261 if (edgeHat.y() > VSMALL)
262 {
263 visitedPoint[curBPoints[pointi]] = 1;
264 nVisited++;
265
266 streamFunction[curBPoints[pointi]] =
267 currentBStream
268 + phi.boundaryField()[patchNo][faceNo]
269 *sign(nHat.x());
270 }
271 else if (edgeHat.y() < -VSMALL)
272 {
273 visitedPoint[curBPoints[pointi]] = 1;
274 nVisited++;
275
276 streamFunction[curBPoints[pointi]] =
277 currentBStream
278 - phi.boundaryField()[patchNo][faceNo]
279 *sign(nHat.x());
280 }
281 else
282 {
283 if (edgeHat.x() > VSMALL)
284 {
285 visitedPoint[curBPoints[pointi]] = 1;
286 nVisited++;
287
288 streamFunction[curBPoints[pointi]] =
289 currentBStream
290 + phi.boundaryField()[patchNo][faceNo]
291 *sign(nHat.y());
292 }
293 else if (edgeHat.x() < -VSMALL)
294 {
295 visitedPoint[curBPoints[pointi]] = 1;
296 nVisited++;
297
298 streamFunction[curBPoints[pointi]] =
299 currentBStream
300 - phi.boundaryField()[patchNo][faceNo]
301 *sign(nHat.y());
302 }
303 }
304 }
305 }
306 }
307 }
308 else
309 {
310 finished = false;
311 }
312 }
313
314 for (label facei=0; facei<nInternalFaces; facei++)
315 {
316 // Get the list of point labels for the face
317 const labelList& curPoints = faces[facei];
318
319 bool pointFound = false;
320
321 scalar currentStream = 0.0;
322 point currentStreamPoint(0, 0, 0);
323
324 forAll(curPoints, pointi)
325 {
326 // Check if the point has been visited
327 if (visitedPoint[curPoints[pointi]] == 1)
328 {
329 // The point has been visited
330 currentStream = streamFunction[curPoints[pointi]];
331 currentStreamPoint = points[curPoints[pointi]];
332 pointFound = true;
333
334 break;
335 }
336 }
337
338 if (pointFound)
339 {
340 // Sort out other points on the face
341 forAll(curPoints, pointi)
342 {
343 // Check if the point has been visited
344 if (visitedPoint[curPoints[pointi]] == 0)
345 {
346 vector edgeHat =
347 points[curPoints[pointi]] - currentStreamPoint;
348
349 edgeHat.replace(slabDir, 0);
350 edgeHat.normalise();
351
352 vector nHat = unitAreas[facei];
353
354 if (edgeHat.y() > VSMALL)
355 {
356 visitedPoint[curPoints[pointi]] = 1;
357 nVisited++;
358
359 streamFunction[curPoints[pointi]] =
360 currentStream
361 + phi[facei]*sign(nHat.x());
362 }
363 else if (edgeHat.y() < -VSMALL)
364 {
365 visitedPoint[curPoints[pointi]] = 1;
366 nVisited++;
367
368 streamFunction[curPoints[pointi]] =
369 currentStream
370 - phi[facei]*sign(nHat.x());
371 }
372 }
373 }
374 }
375 else
376 {
377 finished = false;
378 }
379 }
380
381 if (nVisited == nVisitedOld)
382 {
383 // Find new seed. This must be a
384 // multiply connected domain
385 Log << " Exhausted a seed, looking for new seed "
386 << "(this is correct for multiply connected domains).";
387
388 break;
389 }
390 else
391 {
392 nVisitedOld = nVisited;
393 }
394 } while (!finished);
395 } while (!finished);
396
397 // Normalise the stream-function by the 2D mesh thickness
398 streamFunction /= thickness;
399 streamFunction.boundaryFieldRef() = 0.0;
400
401 return tstreamFunction;
402}
403
404
405bool Foam::functionObjects::streamFunction::calc()
406{
407 if (foundObject<surfaceScalarField>(fieldName_))
408 {
409 const surfaceScalarField& phi =
410 mesh_.lookupObject<surfaceScalarField>(fieldName_);
411
412 return store(resultName_, calc(phi));
413 }
414
415 return false;
416}
417
418
419// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
420
422(
423 const word& name,
424 const Time& runTime,
425 const dictionary& dict
426)
427:
429{
430 setResultName(typeName, "phi");
431
432 const label nD = mesh_.nGeometricD();
433
434 if (nD != 2)
435 {
437 << "Case is not 2D, stream-function cannot be computed"
438 << exit(FatalError);
439 }
440}
441
442
443// ************************************************************************* //
#define Log
Definition: PDRblock.C:35
bool found
Macros for easy insertion into run-time selection tables.
#define addToRunTimeSelectionTable(baseType, thisType, argNames)
Add to construction table with typeName as the key.
surfaceScalarField & phi
const dimensionSet & dimensions() const
Return dimensions.
Internal & ref(const bool updateAccessTime=true)
Return a reference to the dimensioned internal field.
const Boundary & boundaryField() const
Return const-reference to the boundary field.
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
Definition: Time.H:80
static autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
Definition: Time.C:717
static word timeName(const scalar t, const int precision=precision_)
Definition: Time.C:780
vector span() const
The bounding box span (from minimum to maximum)
Definition: boundBoxI.H:127
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
Abstract base-class for Time/database function objects.
const word & name() const noexcept
Return the name of this functionObject.
virtual const word & type() const =0
Runtime type information.
Intermediate class for handling field expression function objects (e.g. blendingFactor etc....
void setResultName(const word &typeName, const word &defaultArg)
Set the name of result field.
const fvMesh & mesh_
Reference to the fvMesh.
Computes the stream function (i.e. https://w.wiki/Ncm).
const Time & time_
Reference to the time database.
label whichPatch(const label faceIndex) const
Return patch index for a given face label.
label nGeometricD() const
Return the number of valid geometric dimensions in the mesh.
Definition: polyMesh.C:883
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1108
const boundBox & bounds() const
Return mesh bounding box.
Definition: polyMesh.H:462
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:456
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1083
const Vector< label > & geometricD() const
Return the vector of geometric directions in mesh.
Definition: polyMesh.C:872
label nInternalFaces() const noexcept
Number of internal faces.
label nPoints() const noexcept
Number of mesh points.
const vectorField & faceAreas() const
const cellList & cells() const
A class for managing temporary objects.
Definition: tmp.H:65
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
A class for handling words, derived from Foam::string.
Definition: word.H:68
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
const polyBoundaryMesh & patches
engineTime & runTime
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
const pointField & points
const dimensionedScalar c
Speed of light in a vacuum.
Namespace for OpenFOAM.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
PtrList< polyPatch > polyPatchList
Store lists of polyPatch as a PtrList.
Definition: polyPatch.H:63
GeometricField< scalar, pointPatchField, pointMesh > pointScalarField
dimensionedScalar sign(const dimensionedScalar &ds)
List< label > labelList
A List of labels.
Definition: List.H:66
List< cell > cellList
A List of cells.
Definition: cellListFwd.H:47
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
vector point
Point is a vector.
Definition: point.H:43
PrimitivePatch< SubList< face >, const pointField & > primitivePatch
A PrimitivePatch with a SubList addressing for the faces, const reference for the point field.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
UList< face > faceUList
A UList of faces.
Definition: faceListFwd.H:46
Field< vector > vectorField
Specialisation of Field<T> for vector.
uint8_t direction
Definition: direction.H:56
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
error FatalError
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
Vector< scalar > vector
Definition: vector.H:61
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
dictionary dict
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
Foam::surfaceFields.