src/sampling/probes/probes.C

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