applications/utilities/postProcessing/velocityField/streamFunction/streamFunction.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation, either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  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 \*---------------------------------------------------------------------------*/
  31
  32 #include "fvCFD.H"
  33 #include "pointFields.H"
  34 #include "emptyPolyPatch.H"
  35 #include "symmetryPolyPatch.H"
  36 #include "wedgePolyPatch.H"
  37 #include "OSspecific.H"
  38
  39 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  40 //  Main program:
  41
  42 int main(int argc, char *argv[])
  43 {
  44     timeSelector::addOptions();
  45
  46 #   include "setRootCase.H"
  47 #   include "createTime.H"
  48
  49     instantList timeDirs = timeSelector::select0(runTime, args);
  50
  51 #   include "createMeshNoClear.H"
  52
  53     pointMesh pMesh(mesh);
  54
  55     forAll(timeDirs, timeI)
  56     {
  57         runTime.setTime(timeDirs[timeI], timeI);
  58
  59         Info<< nl << "Time: " << runTime.timeName() << endl;
  60
  61         IOobject phiHeader
  62         (
  63             "phi",
  64             runTime.timeName(),
  65             mesh,
  66             IOobject::NO_READ
  67         );
  68
  69         if (phiHeader.headerOk())
  70         {
  71             mesh.readUpdate();
  72
  73             Info<< nl << "Reading field phi" << endl;
  74
  75             surfaceScalarField phi
  76             (
  77                 IOobject
  78                 (
  79                     "phi",
  80                     runTime.timeName(),
  81                     mesh,
  82                     IOobject::MUST_READ,
  83                     IOobject::NO_WRITE
  84                 ),
  85                 mesh
  86             );
  87
  88             pointScalarField streamFunction
  89             (
  90                 IOobject
  91                 (
  92                     "streamFunction",
  93                     runTime.timeName(),
  94                     mesh,
  95                     IOobject::NO_READ,
  96                     IOobject::NO_WRITE
  97                 ),
  98                 pMesh,
  99                 dimensionedScalar("zero", phi.dimensions(), 0.0)
 100             );
 101
 102             labelList visitedPoint(mesh.nPoints());
 103             forAll (visitedPoint, pointI)
 104             {
 105                 visitedPoint[pointI] = 0;
 106             }
 107             label nVisited = 0;
 108             label nVisitedOld = 0;
 109
 110             const unallocFaceList& faces = mesh.faces();
 111             const pointField& points = mesh.points();
 112
 113             label nInternalFaces = mesh.nInternalFaces();
 114
 115             vectorField unitAreas = mesh.faceAreas();
 116             unitAreas /= mag(unitAreas);
 117
 118             const polyPatchList& patches = mesh.boundaryMesh();
 119
 120             bool finished = true;
 121
 122             // Find the boundary face with zero flux. set the stream function
 123             // to zero on that face
 124             bool found = false;
 125
 126             do
 127             {
 128                 found = false;
 129
 130                 forAll (patches, patchI)
 131                 {
 132                     const primitivePatch& bouFaces = patches[patchI];
 133
 134                     if (!isType<emptyPolyPatch>(patches[patchI]))
 135                     {
 136                         forAll (bouFaces, faceI)
 137                         {
 138                             if
 139                             (
 140                                 magSqr(phi.boundaryField()[patchI][faceI])
 141                               < SMALL
 142                             )
 143                             {
 144                                 const labelList& zeroPoints = bouFaces[faceI];
 145
 146                                 // Zero flux face found
 147                                 found = true;
 148
 149                                 forAll (zeroPoints, pointI)
 150                                 {
 151                                     if (visitedPoint[zeroPoints[pointI]] == 1)
 152                                     {
 153                                         found = false;
 154                                         break;
 155                                     }
 156                                 }
 157
 158                                 if (found)
 159                                 {
 160                                     Info<< "Zero face: patch: " << patchI
 161                                         << "    face: " << faceI << endl;
 162
 163                                     forAll (zeroPoints, pointI)
 164                                     {
 165                                         streamFunction[zeroPoints[pointI]] = 0;
 166                                         visitedPoint[zeroPoints[pointI]] = 1;
 167                                         nVisited++;
 168                                     }
 169
 170                                     break;
 171                                 }
 172                             }
 173                         }
 174                     }
 175
 176                     if (found) break;
 177                 }
 178
 179                 if (!found)
 180                 {
 181                     Info<< "zero flux boundary face not found. "
 182                         << "Using cell as a reference."
 183                         << endl;
 184
 185                     const cellList& c = mesh.cells();
 186
 187                     forAll (c, cI)
 188                     {
 189                         labelList zeroPoints = c[cI].labels(mesh.faces());
 190
 191                         bool found = true;
 192
 193                         forAll (zeroPoints, pointI)
 194                         {
 195                             if (visitedPoint[zeroPoints[pointI]] == 1)
 196                             {
 197                                 found = false;
 198                                 break;
 199                             }
 200                         }
 201
 202                         if (found)
 203                         {
 204                             forAll (zeroPoints, pointI)
 205                             {
 206                                 streamFunction[zeroPoints[pointI]] = 0.0;
 207                                 visitedPoint[zeroPoints[pointI]] = 1;
 208                                 nVisited++;
 209                             }
 210
 211                             break;
 212                         }
 213                         else
 214                         {
 215                             FatalErrorIn(args.executable())
 216                                 << "Cannot find initialisation face or a cell."
 217                                 << abort(FatalError);
 218                         }
 219                     }
 220                 }
 221
 222                 // Loop through all faces. If one of the points on
 223                 // the face has the streamfunction value different
 224                 // from -1, all points with -1 ont that face have the
 225                 // streamfunction value equal to the face flux in
 226                 // that point plus the value in the visited point
 227                 do
 228                 {
 229                      finished = true;
 230
 231                      for
 232                      (
 233                          label faceI = nInternalFaces;
 234                          faceI<faces.size();
 235                          faceI++
 236                      )
 237                      {
 238                          const labelList& curBPoints = faces[faceI];
 239                          bool bPointFound = false;
 240
 241                          scalar currentBStream = 0.0;
 242                          vector currentBStreamPoint(0, 0, 0);
 243
 244                          forAll (curBPoints, pointI)
 245                          {
 246                              // Check if the point has been visited
 247                              if (visitedPoint[curBPoints[pointI]] == 1)
 248                              {
 249                                  // The point has been visited
 250                                  currentBStream =
 251                                      streamFunction[curBPoints[pointI]];
 252                                  currentBStreamPoint =
 253                                      points[curBPoints[pointI]];
 254
 255                                  bPointFound = true;
 256
 257                                  break;
 258                              }
 259                          }
 260
 261                          if (bPointFound)
 262                          {
 263                              // Sort out other points on the face
 264                              forAll (curBPoints, pointI)
 265                              {
 266                                  // Check if the point has been visited
 267                                  if (visitedPoint[curBPoints[pointI]] == 0)
 268                                  {
 269                                      label patchNo =
 270                                          mesh.boundaryMesh().whichPatch(faceI);
 271
 272                                      if
 273                                      (
 274                                         !isType<emptyPolyPatch>
 275                                          (patches[patchNo])
 276                                      && !isType<symmetryPolyPatch>
 277                                          (patches[patchNo])
 278                                      && !isType<wedgePolyPatch>
 279                                          (patches[patchNo])
 280                                      )
 281                                      {
 282                                          label faceNo =
 283                                              mesh.boundaryMesh()[patchNo]
 284                                              .whichFace(faceI);
 285
 286                                          vector edgeHat =
 287                                              points[curBPoints[pointI]]
 288                                              - currentBStreamPoint;
 289                                          edgeHat.replace(vector::Z, 0);
 290                                          edgeHat /= mag(edgeHat);
 291
 292                                          vector nHat = unitAreas[faceI];
 293
 294                                          if (edgeHat.y() > VSMALL)
 295                                          {
 296                                              visitedPoint[curBPoints[pointI]] =
 297                                                  1;
 298                                              nVisited++;
 299
 300                                              streamFunction[curBPoints[pointI]]
 301                                                  =
 302                                                  currentBStream
 303                                                + phi.boundaryField()
 304                                                  [patchNo][faceNo]
 305                                                  *sign(nHat.x());
 306                                          }
 307                                          else if (edgeHat.y() < -VSMALL)
 308                                          {
 309                                              visitedPoint[curBPoints[pointI]] =
 310                                                  1;
 311                                              nVisited++;
 312
 313                                              streamFunction[curBPoints[pointI]]
 314                                                  =
 315                                                  currentBStream
 316                                                - phi.boundaryField()
 317                                                  [patchNo][faceNo]
 318                                                  *sign(nHat.x());
 319                                          }
 320                                          else
 321                                          {
 322                                              if (edgeHat.x() > VSMALL)
 323                                              {
 324                                                  visitedPoint
 325                                                      [curBPoints[pointI]] = 1;
 326                                                  nVisited++;
 327
 328                                                  streamFunction
 329                                                      [curBPoints[pointI]] =
 330                                                      currentBStream
 331                                                    + phi.boundaryField()
 332                                                      [patchNo][faceNo]
 333                                                      *sign(nHat.y());
 334                                              }
 335                                              else if (edgeHat.x() < -VSMALL)
 336                                              {
 337                                                  visitedPoint
 338                                                      [curBPoints[pointI]] = 1;
 339                                                  nVisited++;
 340
 341                                                  streamFunction
 342                                                      [curBPoints[pointI]] =
 343                                                      currentBStream
 344                                                    - phi.boundaryField()
 345                                                      [patchNo][faceNo]
 346                                                      *sign(nHat.y());
 347                                              }
 348                                          }
 349                                      }
 350                                  }
 351                              }
 352                          }
 353                          else
 354                          {
 355                              finished = false;
 356                          }
 357                      }
 358
 359                      for (label faceI=0; faceI<nInternalFaces; faceI++)
 360                      {
 361                          // Get the list of point labels for the face
 362                          const labelList& curPoints = faces[faceI];
 363
 364                          bool pointFound = false;
 365
 366                          scalar currentStream = 0.0;
 367                          point currentStreamPoint(0, 0, 0);
 368
 369                          forAll (curPoints, pointI)
 370                          {
 371                              // Check if the point has been visited
 372                              if (visitedPoint[curPoints[pointI]] == 1)
 373                              {
 374                                  // The point has been visited
 375                                  currentStream =
 376                                      streamFunction[curPoints[pointI]];
 377                                  currentStreamPoint =
 378                                      points[curPoints[pointI]];
 379                                  pointFound = true;
 380
 381                                  break;
 382                              }
 383                          }
 384
 385                          if (pointFound)
 386                          {
 387                              // Sort out other points on the face
 388                              forAll (curPoints, pointI)
 389                              {
 390                                  // Check if the point has been visited
 391                                  if (visitedPoint[curPoints[pointI]] == 0)
 392                                  {
 393                                      vector edgeHat =
 394                                          points[curPoints[pointI]]
 395                                        - currentStreamPoint;
 396
 397                                      edgeHat.replace(vector::Z, 0);
 398                                      edgeHat /= mag(edgeHat);
 399
 400                                      vector nHat = unitAreas[faceI];
 401
 402                                      if (edgeHat.y() > VSMALL)
 403                                      {
 404                                          visitedPoint[curPoints[pointI]] = 1;
 405                                          nVisited++;
 406
 407                                          streamFunction[curPoints[pointI]] =
 408                                              currentStream
 409                                            + phi[faceI]*sign(nHat.x());
 410                                      }
 411                                      else if (edgeHat.y() < -VSMALL)
 412                                      {
 413                                          visitedPoint[curPoints[pointI]] = 1;
 414                                          nVisited++;
 415
 416                                          streamFunction[curPoints[pointI]] =
 417                                              currentStream
 418                                            - phi[faceI]*sign(nHat.x());
 419                                      }
 420                                  }
 421                              }
 422                          }
 423                          else
 424                          {
 425                              finished = false;
 426                          }
 427                      }
 428
 429                      Info<< ".";
 430 //                     Info<< "One pass, n visited = " << nVisited << endl;
 431
 432                      if (nVisited == nVisitedOld)
 433                      {
 434                          // Find new seed.  This must be a
 435                          // multiply connected domain
 436                          Info<< nl << "Exhausted a seed. Looking for new seed "
 437                              << "(this is correct for multiply connected "
 438                              << "domains).";
 439
 440                          break;
 441                      }
 442                      else
 443                      {
 444                          nVisitedOld = nVisited;
 445                      }
 446                 } while (!finished);
 447
 448                 Info << endl;
 449             } while (!finished);
 450
 451             streamFunction.boundaryField() = 0.0;
 452             streamFunction.write();
 453         }
 454         else
 455         {
 456             WarningIn(args.executable())
 457                 << "Flux field does not exist."
 458                 << " Stream function not calculated" << endl;
 459         }
 460     }
 461
 462     Info<< "\nEnd\n" << endl;
 463
 464     return 0;
 465 }
 466
 467
 468 // ************************************************************************* //