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

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