applications/solvers/incompressible/boundaryFoam/createFields.H

   1     Info<< "Reading field U\n" << endl;
   2     volVectorField U
   3     (
   4         IOobject
   5         (
   6             "U",
   7             runTime.timeName(),
   8             mesh,
   9             IOobject::MUST_READ,
  10             IOobject::AUTO_WRITE
  11         ),
  12         mesh
  13     );
  14
  15
  16     Info<< "Creating face flux\n" << endl;
  17     surfaceScalarField phi
  18     (
  19         IOobject
  20         (
  21             "phi",
  22             runTime.timeName(),
  23             mesh,
  24             IOobject::NO_READ,
  25             IOobject::NO_WRITE
  26         ),
  27         mesh,
  28         dimensionedScalar("zero", mesh.Sf().dimensions()*U.dimensions(), 0.0)
  29     );
  30
  31
  32     singlePhaseTransportModel laminarTransport(U, phi);
  33
  34     autoPtr<incompressible::RASModel> turbulence
  35     (
  36         incompressible::RASModel::New(U, phi, laminarTransport)
  37     );
  38
  39
  40     IOdictionary transportProperties
  41     (
  42         IOobject
  43         (
  44             "transportProperties",
  45             runTime.constant(),
  46             mesh,
  47             IOobject::MUST_READ,
  48             IOobject::NO_WRITE
  49         )
  50     );
  51
  52     dimensionedVector Ubar
  53     (
  54         transportProperties.lookup("Ubar")
  55     );
  56
  57     vector flowDirection = (Ubar/mag(Ubar)).value();
  58     tensor flowMask = sqr(flowDirection);
  59
  60
  61     // Search for wall patches faces and store normals
  62
  63     scalar nWallFaces(0);
  64     vector wallNormal(vector::zero);
  65
  66     const fvPatchList& patches = mesh.boundary();
  67
  68     forAll(patches, patchi)
  69     {
  70         const fvPatch& currPatch = patches[patchi];
  71
  72         if (isA<wallFvPatch>(currPatch))
  73         {
  74             forAll(currPatch, facei)
  75             {
  76                 nWallFaces++;
  77
  78                 if (nWallFaces == 1)
  79                 {
  80                     wallNormal =
  81                       - mesh.Sf().boundaryField()[patchi][facei]
  82                         /mesh.magSf().boundaryField()[patchi][facei];
  83                 }
  84                 else if (nWallFaces == 2)
  85                 {
  86                     vector wallNormal2 =
  87                         mesh.Sf().boundaryField()[patchi][facei]
  88                         /mesh.magSf().boundaryField()[patchi][facei];
  89
  90                     //- Check that wall faces are parallel
  91                     if
  92                     (
  93                         mag(wallNormal & wallNormal2) > 1.01
  94                       ||mag(wallNormal & wallNormal2) < 0.99
  95                     )
  96                     {
  97                         Info<< "boundaryFoam: wall faces are not parallel"
  98                             << endl
  99                             << abort(FatalError);
 100                     }
 101                 }
 102                 else
 103                 {
 104                     Info<< "boundaryFoam: number of wall faces > 2"
 105                         << endl
 106                         << abort(FatalError);
 107                 }
 108             }
 109         }
 110     }
 111
 112
 113     //- create position array for graph generation
 114     scalarField y = wallNormal & mesh.C().internalField();
 115
 116
 117     dimensionedVector gradP
 118     (
 119         "gradP",
 120         dimensionSet(0, 1, -2, 0, 0),
 121         vector(0, 0, 0)
 122     );