73 int main(
int argc,
char *argv[])
75 timeSelector::addOptions();
88 runTime.setTime(timeDirs[timeI], timeI);
90 Info<<
nl <<
"Time: " << runTime.timeName() <<
endl;
100 if (phiHeader.headerOk())
134 forAll (visitedPoint, pointI)
136 visitedPoint[pointI] = 0;
139 label nVisitedOld = 0;
147 unitAreas /=
mag(unitAreas);
151 bool finished =
true;
165 if (!isType<emptyPolyPatch>(patches[patchI]))
171 magSqr(
phi.boundaryField()[patchI][faceI])
175 const labelList& zeroPoints = bouFaces[faceI];
180 forAll (zeroPoints, pointI)
182 if (visitedPoint[zeroPoints[pointI]] == 1)
191 Info<<
"Zero face: patch: " << patchI
192 <<
" face: " << faceI <<
endl;
194 forAll (zeroPoints, pointI)
196 streamFunction[zeroPoints[pointI]] = 0;
197 visitedPoint[zeroPoints[pointI]] = 1;
212 Info<<
"zero flux boundary face not found. "
213 <<
"Using cell as a reference."
224 forAll (zeroPoints, pointI)
226 if (visitedPoint[zeroPoints[pointI]] == 1)
235 forAll (zeroPoints, pointI)
237 streamFunction[zeroPoints[pointI]] = 0.0;
238 visitedPoint[zeroPoints[pointI]] = 1;
247 <<
"Cannot find initialisation face or a cell."
264 label faceI = nInternalFaces;
269 const labelList& curBPoints = faces[faceI];
270 bool bPointFound =
false;
272 scalar currentBStream = 0.0;
273 vector currentBStreamPoint(0, 0, 0);
275 forAll (curBPoints, pointI)
278 if (visitedPoint[curBPoints[pointI]] == 1)
282 streamFunction[curBPoints[pointI]];
283 currentBStreamPoint =
284 points[curBPoints[pointI]];
295 forAll (curBPoints, pointI)
298 if (visitedPoint[curBPoints[pointI]] == 0)
305 !isType<emptyPolyPatch>
307 && !isType<symmetryPolyPatch>
309 && !isType<wedgePolyPatch>
318 points[curBPoints[pointI]]
319 - currentBStreamPoint;
321 edgeHat /=
mag(edgeHat);
323 vector nHat = unitAreas[faceI];
325 if (edgeHat.
y() > VSMALL)
327 visitedPoint[curBPoints[pointI]] =
331 streamFunction[curBPoints[pointI]]
334 +
phi.boundaryField()
338 else if (edgeHat.
y() < -VSMALL)
340 visitedPoint[curBPoints[pointI]] =
344 streamFunction[curBPoints[pointI]]
347 -
phi.boundaryField()
353 if (edgeHat.
x() > VSMALL)
356 [curBPoints[pointI]] = 1;
360 [curBPoints[pointI]] =
362 +
phi.boundaryField()
366 else if (edgeHat.
x() < -VSMALL)
369 [curBPoints[pointI]] = 1;
373 [curBPoints[pointI]] =
375 -
phi.boundaryField()
390 for (label faceI=0; faceI<nInternalFaces; faceI++)
393 const labelList& curPoints = faces[faceI];
395 bool pointFound =
false;
397 scalar currentStream = 0.0;
398 point currentStreamPoint(0, 0, 0);
400 forAll (curPoints, pointI)
403 if (visitedPoint[curPoints[pointI]] == 1)
407 streamFunction[curPoints[pointI]];
409 points[curPoints[pointI]];
419 forAll (curPoints, pointI)
422 if (visitedPoint[curPoints[pointI]] == 0)
425 points[curPoints[pointI]]
426 - currentStreamPoint;
429 edgeHat /=
mag(edgeHat);
431 vector nHat = unitAreas[faceI];
433 if (edgeHat.
y() > VSMALL)
435 visitedPoint[curPoints[pointI]] = 1;
438 streamFunction[curPoints[pointI]] =
442 else if (edgeHat.
y() < -VSMALL)
444 visitedPoint[curPoints[pointI]] = 1;
447 streamFunction[curPoints[pointI]] =
463 if (nVisited == nVisitedOld)
467 Info<<
nl <<
"Exhausted a seed. Looking for new seed "
468 <<
"(this is correct for multiply connected "
475 nVisitedOld = nVisited;
482 streamFunction.boundaryField() = 0.0;
483 streamFunction.
write();
488 <<
"Flux field does not exist."
489 <<
" Stream function not calculated" << endl;