applications/test/nearWallDist-wave/testWallDist2.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2008 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 Description
  26     Wave propagation of nearwall distance through grid. Every iteration
  27     information goes through one layer of cells.
  28
  29 \*---------------------------------------------------------------------------*/
  30
  31 #include "fvCFD.H"
  32 #include "wallFvPatch.H"
  33 #include "MeshWave.H"
  34 #include "wallPoint.H"
  35
  36
  37 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  38 // Main program:
  39
  40 int main(int argc, char *argv[])
  41 {
  42 #   include "setRootCase.H"
  43 #   include "createTime.H"
  44 #   include "createMesh.H"
  45
  46     Info<< "Mesh read in = "
  47         << runTime.cpuTimeIncrement()
  48         << " s\n" << endl << endl;
  49
  50     Info<< "Creating field wDistNC\n" << endl;
  51     volScalarField wallDistUncorrected
  52     (
  53         IOobject
  54         (
  55             "wDistNC",
  56             runTime.timeName(),
  57             mesh,
  58             IOobject::NO_READ,
  59             IOobject::AUTO_WRITE
  60         ),
  61         mesh,
  62         dimensionedScalar
  63         (
  64             "wallDist",
  65             dimensionSet(0, 1, 0, 0, 0),
  66             0.0
  67         )
  68     );
  69
  70     //
  71     // Set initial changed faces: set wallPoint for wall faces to wall centre
  72     //
  73
  74     // Count walls
  75     label nWalls = 0;
  76     forAll (mesh.boundary(), patchI)
  77     {
  78         const fvPatch& patch = mesh.boundary()[patchI];
  79
  80         if (typeid(patch) == typeid(wallFvPatch))
  81         {
  82             nWalls += patch.size();
  83         }
  84     }
  85
  86     List<wallPoint> faceDist(nWalls);
  87     labelList changedFaces(nWalls);
  88
  89     label nChangedFaces = 0;
  90     forAll (mesh.boundary(), patchI)
  91     {
  92         const fvPatch& patch = mesh.boundary()[patchI];
  93
  94         if (typeid(patch) == typeid(wallFvPatch))
  95         {
  96             forAll (patch.Cf(), patchFaceI)
  97             {
  98                 const polyPatch& polyPatch = mesh.boundaryMesh()[patchI];
  99
 100                 label meshFaceI = polyPatch.start() + patchFaceI;
 101
 102                 changedFaces[nChangedFaces] = meshFaceI;
 103
 104                 faceDist[nChangedFaces] =
 105                     wallPoint(patch.Cf()[patchFaceI], 0.0);
 106
 107                 nChangedFaces++;
 108             }
 109         }
 110     }
 111
 112
 113     MeshWave<wallPoint> wallDistCalc
 114     (
 115         mesh,
 116         changedFaces,
 117         faceDist,
 118         0           // max iterations
 119     );
 120
 121     Info<< "\nStarting time loop\n" << endl;
 122
 123     for (runTime++; !runTime.end(); runTime++)
 124     {
 125         Info<< "Time = " << runTime.timeName() << endl;
 126
 127
 128         label nCells = wallDistCalc.faceToCell();
 129
 130         Info<< "    Total changed cells   : " << nCells << endl;
 131
 132         if (nCells == 0)
 133         {
 134             break;
 135         }
 136
 137
 138         label nFaces = wallDistCalc.cellToFace();
 139
 140         Info<< "    Total changed faces   : " << nFaces << endl;
 141
 142         if (nFaces == 0)
 143         {
 144             break;
 145         }
 146
 147
 148         //
 149         // Copy face and cell values into field
 150         //
 151
 152         const List<wallPoint>& cellInfo = wallDistCalc.allCellInfo();
 153         const List<wallPoint>& faceInfo = wallDistCalc.allFaceInfo();
 154
 155         label nIllegal = 0;
 156
 157         // Copy cell values
 158         forAll(cellInfo, cellI)
 159         {
 160             scalar dist = cellInfo[cellI].distSqr();
 161             if (cellInfo[cellI].valid())
 162             {
 163                 wallDistUncorrected[cellI] = Foam::sqrt(dist);
 164             }
 165             else
 166             {
 167                 wallDistUncorrected[cellI] = -1;
 168                 nIllegal++;
 169             }
 170         }
 171
 172         // Copy boundary values
 173         forAll (wallDistUncorrected.boundaryField(), patchI)
 174         {
 175             fvPatchScalarField& patchField =
 176                 wallDistUncorrected.boundaryField()[patchI];
 177
 178             forAll(patchField, patchFaceI)
 179             {
 180                 label meshFaceI =
 181                     patchField.patch().patch().start() + patchFaceI;
 182
 183                 scalar dist = faceInfo[meshFaceI].distSqr();
 184                 if (faceInfo[meshFaceI].valid())
 185                 {
 186                     patchField[patchFaceI] = Foam::sqrt(dist);
 187                 }
 188                 else
 189                 {
 190                     patchField[patchFaceI] = dist;
 191                     nIllegal++;
 192                 }
 193             }
 194         }
 195
 196         Info<< "nIllegal:" << nIllegal << endl;
 197
 198
 199         //
 200         // Write it
 201         //
 202
 203         wallDistUncorrected.write();
 204
 205         Info<< "ExecutionTime = "
 206             << runTime.elapsedCpuTime()
 207             << " s\n" << endl << endl;
 208     }
 209
 210     Info<< "End\n" << endl;
 211
 212     return 0;
 213 }
 214
 215
 216 // ************************************************************************* //