src/sampling/probes/probes.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 \*---------------------------------------------------------------------------*/
  25
  26 #include "probes.H"
  27 #include "volFields.H"
  28 #include "dictionary.H"
  29 #include "foamTime.H"
  30 #include "IOmanip.H"
  31
  32 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  33
  34 namespace Foam
  35 {
  36     defineTypeNameAndDebug(probes, 0);
  37 }
  38
  39 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  40
  41 void Foam::probes::findCells(const fvMesh& mesh)
  42 {
  43     if (cellList_.empty())
  44     {
  45         Info<< "Searching for probe point locations" << endl;
  46
  47         cellList_.setSize(probeLocations_.size());
  48
  49         forAll(probeLocations_, probeI)
  50         {
  51             cellList_[probeI] = mesh.findCell(probeLocations_[probeI]);
  52
  53             if (debug && cellList_[probeI] != -1)
  54             {
  55                 Pout<< "probes : found point " << probeLocations_[probeI]
  56                     << " in cell " << cellList_[probeI] << endl;
  57             }
  58         }
  59
  60
  61         // Check if all probes have been found.
  62         forAll(cellList_, probeI)
  63         {
  64             label cellI = cellList_[probeI];
  65
  66             // Check at least one processor with cell.
  67             reduce(cellI, maxOp<label>());
  68
  69             if (cellI == -1)
  70             {
  71                 if (Pstream::master())
  72                 {
  73                     WarningIn("probes::read()")
  74                         << "Did not find location " << probeLocations_[probeI]
  75                         << " in any cell. Skipping location." << endl;
  76                 }
  77             }
  78             else
  79             {
  80                 // Make sure location not on two domains.
  81                 if (cellList_[probeI] != -1 && cellList_[probeI] != cellI)
  82                 {
  83                     WarningIn("probes::read()")
  84                         << "Location " << probeLocations_[probeI]
  85                         << " seems to be on multiple domains:"
  86                         << " cell " << cellList_[probeI]
  87                         << " on my domain " << Pstream::myProcNo()
  88                         << " and cell " << cellI << " on some other domain."
  89                         << endl
  90                         << "This might happen if the probe location is on"
  91                         << " a processor patch. Change the location slightly"
  92                         << " to prevent this." << endl;
  93                 }
  94             }
  95         }
  96     }
  97 }
  98
  99
 100 bool Foam::probes::checkFieldTypes()
 101 {
 102     wordList fieldTypes(fieldNames_.size());
 103
 104     // check files for a particular time
 105     if (loadFromFiles_)
 106     {
 107         forAll(fieldNames_, fieldI)
 108         {
 109             IOobject io
 110             (
 111                 fieldNames_[fieldI],
 112                 obr_.time().timeName(),
 113                 refCast<const polyMesh>(obr_),
 114                 IOobject::MUST_READ,
 115                 IOobject::NO_WRITE,
 116                 false
 117             );
 118
 119             if (io.headerOk())
 120             {
 121                 fieldTypes[fieldI] = io.headerClassName();
 122             }
 123             else
 124             {
 125                 fieldTypes[fieldI] = "(notFound)";
 126             }
 127         }
 128     }
 129     else
 130     {
 131         // check objectRegistry
 132         forAll(fieldNames_, fieldI)
 133         {
 134             objectRegistry::const_iterator iter =
 135                 obr_.find(fieldNames_[fieldI]);
 136
 137             if (iter != obr_.end())
 138             {
 139                 fieldTypes[fieldI] = iter()->type();
 140             }
 141             else
 142             {
 143                 fieldTypes[fieldI] = "(notFound)";
 144             }
 145         }
 146     }
 147
 148
 149     label nFields = 0;
 150
 151     // classify fieldTypes
 152     nFields += countFields(scalarFields_, fieldTypes);
 153     nFields += countFields(vectorFields_, fieldTypes);
 154     nFields += countFields(sphericalTensorFields_, fieldTypes);
 155     nFields += countFields(symmTensorFields_, fieldTypes);
 156     nFields += countFields(tensorFields_, fieldTypes);
 157
 158     // concatenate all the lists into foundFields
 159     wordList foundFields(nFields);
 160
 161     label fieldI = 0;
 162     forAll(scalarFields_, i)
 163     {
 164         foundFields[fieldI++] = scalarFields_[i];
 165     }
 166     forAll(vectorFields_, i)
 167     {
 168         foundFields[fieldI++] = vectorFields_[i];
 169     }
 170     forAll(sphericalTensorFields_, i)
 171     {
 172         foundFields[fieldI++] = sphericalTensorFields_[i];
 173     }
 174     forAll(symmTensorFields_, i)
 175     {
 176         foundFields[fieldI++] = symmTensorFields_[i];
 177     }
 178     forAll(tensorFields_, i)
 179     {
 180         foundFields[fieldI++] = tensorFields_[i];
 181     }
 182
 183     if (Pstream::master())
 184     {
 185         fileName probeDir;
 186
 187         fileName probeSubDir = name_;
 188
 189         if (obr_.name() != polyMesh::defaultRegion)
 190         {
 191             probeSubDir = probeSubDir/obr_.name();
 192         }
 193         probeSubDir = probeSubDir/obr_.time().timeName();
 194
 195         if (Pstream::parRun())
 196         {
 197             // Put in undecomposed case
 198             // (Note: gives problems for distributed data running)
 199             probeDir = obr_.time().path()/".."/probeSubDir;
 200         }
 201         else
 202         {
 203             probeDir = obr_.time().path()/probeSubDir;
 204         }
 205
 206         // Close the file if any fields have been removed.
 207         forAllIter(HashPtrTable<OFstream>, probeFilePtrs_, iter)
 208         {
 209             if (findIndex(foundFields, iter.key()) == -1)
 210             {
 211                 if (debug)
 212                 {
 213                     Pout<< "close stream: " << iter()->name() << endl;
 214                 }
 215
 216                 delete probeFilePtrs_.remove(iter);
 217             }
 218         }
 219
 220         // Open new files for new fields. Keep existing files.
 221
 222         probeFilePtrs_.resize(2*foundFields.size());
 223
 224         forAll(foundFields, fieldI)
 225         {
 226             const word& fldName = foundFields[fieldI];
 227
 228             // Check if added field. If so open a stream for it.
 229
 230             if (!probeFilePtrs_.found(fldName))
 231             {
 232                 // Create directory if does not exist.
 233                 mkDir(probeDir);
 234
 235                 OFstream* sPtr = new OFstream(probeDir/fldName);
 236
 237                 if (debug)
 238                 {
 239                     Pout<< "open  stream: " << sPtr->name() << endl;
 240                 }
 241
 242                 probeFilePtrs_.insert(fldName, sPtr);
 243
 244                 unsigned int w = IOstream::defaultPrecision() + 7;
 245
 246                 for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
 247                 {
 248                     *sPtr<< '#' << setw(IOstream::defaultPrecision() + 6)
 249                         << vector::componentNames[cmpt];
 250
 251                     forAll(probeLocations_, probeI)
 252                     {
 253                         *sPtr<< ' ' << setw(w) << probeLocations_[probeI][cmpt];
 254                     }
 255                     *sPtr << endl;
 256                 }
 257
 258                 *sPtr<< '#' << setw(IOstream::defaultPrecision() + 6)
 259                     << "Time" << endl;
 260             }
 261         }
 262
 263         if (debug)
 264         {
 265             Pout<< "Probing fields:" << foundFields << nl
 266                 << "Probing locations:" << probeLocations_ << nl
 267                 << endl;
 268         }
 269     }
 270
 271
 272     return nFields > 0;
 273 }
 274
 275
 276 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 277
 278 Foam::probes::probes
 279 (
 280     const word& name,
 281     const objectRegistry& obr,
 282     const dictionary& dict,
 283     const bool loadFromFiles
 284 )
 285 :
 286     name_(name),
 287     obr_(obr),
 288     loadFromFiles_(loadFromFiles),
 289     fieldNames_(),
 290     probeLocations_(),
 291     scalarFields_(),
 292     vectorFields_(),
 293     sphericalTensorFields_(),
 294     symmTensorFields_(),
 295     tensorFields_(),
 296     cellList_(),
 297     probeFilePtrs_()
 298 {
 299     read(dict);
 300 }
 301
 302
 303 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 304
 305 Foam::probes::~probes()
 306 {}
 307
 308
 309 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 310
 311 void Foam::probes::execute()
 312 {
 313 }
 314
 315
 316 void Foam::probes::end()
 317 {
 318     // Do nothing - only valid on write
 319 }
 320
 321
 322 void Foam::probes::write()
 323 {
 324     // Check if the mesh is changing and if so, resample
 325     const fvMesh& mesh = refCast<const fvMesh>(obr_);
 326
 327     if (mesh.moving() || mesh.changing())
 328     {
 329         cellList_.clear();
 330         findCells(mesh);
 331     }
 332
 333     if (probeLocations_.size() && checkFieldTypes())
 334     {
 335         sampleAndWrite(scalarFields_);
 336         sampleAndWrite(vectorFields_);
 337         sampleAndWrite(sphericalTensorFields_);
 338         sampleAndWrite(symmTensorFields_);
 339         sampleAndWrite(tensorFields_);
 340     }
 341 }
 342
 343
 344 void Foam::probes::read(const dictionary& dict)
 345 {
 346     dict.lookup("fields") >> fieldNames_;
 347     dict.lookup("probeLocations") >> probeLocations_;
 348
 349     // Force all cell locations to be redetermined
 350     cellList_.clear();
 351     findCells(refCast<const fvMesh>(obr_));
 352     checkFieldTypes();
 353 }
 354
 355
 356 // ************************************************************************* //