src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.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 \*---------------------------------------------------------------------------*/
  26
  27 #include "sampledSurfaces.H"
  28 #include "volFields.H"
  29 #include "dictionary.H"
  30 #include "Time.H"
  31 #include "IOmanip.H"
  32 #include "ListListOps.H"
  33 #include "mergePoints.H"
  34 #include "volPointInterpolation.H"
  35
  36 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  37
  38 namespace Foam
  39 {
  40     //- Used to offset faces in Pstream::combineOffset
  41     template <>
  42     class offsetOp<face>
  43     {
  44
  45     public:
  46
  47         face operator()
  48         (
  49             const face& x,
  50             const label offset
  51         ) const
  52         {
  53             face result(x.size());
  54
  55             forAll(x, xI)
  56             {
  57                 result[xI] = x[xI] + offset;
  58             }
  59             return result;
  60         }
  61     };
  62
  63
  64     defineTypeNameAndDebug(sampledSurfaces, 0);
  65 }
  66
  67
  68 bool Foam::sampledSurfaces::verbose_(false);
  69 Foam::scalar Foam::sampledSurfaces::mergeTol_(1e-10);
  70
  71 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  72
  73 Foam::label Foam::sampledSurfaces::classifyFieldTypes()
  74 {
  75     label nFields = 0;
  76
  77     scalarFields_.clear();
  78     vectorFields_.clear();
  79     sphericalTensorFields_.clear();
  80     symmTensorFields_.clear();
  81     tensorFields_.clear();
  82
  83     forAll(fieldNames_, fieldI)
  84     {
  85         const word& fieldName = fieldNames_[fieldI];
  86         word fieldType = "";
  87
  88         // check files for a particular time
  89         if (loadFromFiles_)
  90         {
  91             IOobject io
  92             (
  93                 fieldName,
  94                 mesh_.time().timeName(),
  95                 mesh_,
  96                 IOobject::MUST_READ,
  97                 IOobject::NO_WRITE,
  98                 false
  99             );
 100
 101             if (io.headerOk())
 102             {
 103                 fieldType = io.headerClassName();
 104             }
 105             else
 106             {
 107                 continue;
 108             }
 109         }
 110         else
 111         {
 112             // check objectRegistry
 113             objectRegistry::const_iterator iter = mesh_.find(fieldName);
 114
 115             if (iter != mesh_.objectRegistry::end())
 116             {
 117                 fieldType = iter()->type();
 118             }
 119             else
 120             {
 121                 continue;
 122             }
 123         }
 124
 125
 126         if (fieldType == volScalarField::typeName)
 127         {
 128             scalarFields_.append(fieldName);
 129             nFields++;
 130         }
 131         else if (fieldType == volVectorField::typeName)
 132         {
 133             vectorFields_.append(fieldName);
 134             nFields++;
 135         }
 136         else if (fieldType == volSphericalTensorField::typeName)
 137         {
 138             sphericalTensorFields_.append(fieldName);
 139             nFields++;
 140         }
 141         else if (fieldType == volSymmTensorField::typeName)
 142         {
 143             symmTensorFields_.append(fieldName);
 144             nFields++;
 145         }
 146         else if (fieldType == volTensorField::typeName)
 147         {
 148             tensorFields_.append(fieldName);
 149             nFields++;
 150         }
 151
 152     }
 153
 154     return nFields;
 155 }
 156
 157
 158 void Foam::sampledSurfaces::writeGeometry() const
 159 {
 160     // Write to time directory under outputPath_
 161     // skip surface without faces (eg, a failed cut-plane)
 162
 163     const fileName outputDir = outputPath_/mesh_.time().timeName();
 164
 165     forAll(*this, surfI)
 166     {
 167         const sampledSurface& s = operator[](surfI);
 168
 169         if (Pstream::parRun())
 170         {
 171             if (Pstream::master() && mergeList_[surfI].faces.size())
 172             {
 173                 genericFormatter_->write
 174                 (
 175                     outputDir,
 176                     s.name(),
 177                     mergeList_[surfI].points,
 178                     mergeList_[surfI].faces
 179                 );
 180             }
 181         }
 182         else if (s.faces().size())
 183         {
 184             genericFormatter_->write
 185             (
 186                 outputDir,
 187                 s.name(),
 188                 s.points(),
 189                 s.faces()
 190             );
 191         }
 192     }
 193 }
 194
 195
 196 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 197
 198 Foam::sampledSurfaces::sampledSurfaces
 199 (
 200     const word& name,
 201     const objectRegistry& obr,
 202     const dictionary& dict,
 203     const bool loadFromFiles
 204 )
 205 :
 206     PtrList<sampledSurface>(),
 207     name_(name),
 208     mesh_(refCast<const fvMesh>(obr)),
 209     loadFromFiles_(loadFromFiles),
 210     outputPath_(fileName::null),
 211     fieldNames_(),
 212     interpolationScheme_(word::null),
 213     writeFormat_(word::null),
 214     mergeList_(),
 215     genericFormatter_(NULL),
 216     scalarFields_(),
 217     vectorFields_(),
 218     sphericalTensorFields_(),
 219     symmTensorFields_(),
 220     tensorFields_()
 221 {
 222     if (Pstream::parRun())
 223     {
 224         outputPath_ = mesh_.time().path()/".."/name_;
 225     }
 226     else
 227     {
 228         outputPath_ = mesh_.time().path()/name_;
 229     }
 230
 231     read(dict);
 232 }
 233
 234
 235 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 236
 237 Foam::sampledSurfaces::~sampledSurfaces()
 238 {}
 239
 240
 241 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 242
 243 void Foam::sampledSurfaces::verbose(const bool verbosity)
 244 {
 245     verbose_ = verbosity;
 246 }
 247
 248
 249 void Foam::sampledSurfaces::execute()
 250 {
 251     // Do nothing - only valid on write
 252 }
 253
 254
 255 void Foam::sampledSurfaces::end()
 256 {
 257     // Do nothing - only valid on write
 258 }
 259
 260
 261 void Foam::sampledSurfaces::write()
 262 {
 263     if (size())
 264     {
 265         // finalize surfaces, merge points etc.
 266         update();
 267
 268         const label nFields = classifyFieldTypes();
 269
 270         if (Pstream::master())
 271         {
 272             if (debug)
 273             {
 274                 Pout<< "timeName = " << mesh_.time().timeName() << nl
 275                     << "scalarFields    " << scalarFields_ << nl
 276                     << "vectorFields    " << vectorFields_ << nl
 277                     << "sphTensorFields " << sphericalTensorFields_ << nl
 278                     << "symTensorFields " << symmTensorFields_ <<nl
 279                     << "tensorFields    " << tensorFields_ <<nl;
 280
 281                 Pout<< "Creating directory "
 282                     << outputPath_/mesh_.time().timeName() << nl << endl;
 283
 284             }
 285
 286             mkDir(outputPath_/mesh_.time().timeName());
 287         }
 288
 289         // write geometry first if required, or when no fields would otherwise
 290         // be written
 291         if (nFields == 0 || genericFormatter_->separateFiles())
 292         {
 293             writeGeometry();
 294         }
 295
 296         sampleAndWrite(scalarFields_);
 297         sampleAndWrite(vectorFields_);
 298         sampleAndWrite(sphericalTensorFields_);
 299         sampleAndWrite(symmTensorFields_);
 300         sampleAndWrite(tensorFields_);
 301     }
 302 }
 303
 304
 305 void Foam::sampledSurfaces::read(const dictionary& dict)
 306 {
 307     fieldNames_ = wordList(dict.lookup("fields"));
 308
 309     const label nFields = fieldNames_.size();
 310
 311     scalarFields_.reset(nFields);
 312     vectorFields_.reset(nFields);
 313     sphericalTensorFields_.reset(nFields);
 314     symmTensorFields_.reset(nFields);
 315     tensorFields_.reset(nFields);
 316
 317     interpolationScheme_ = dict.lookupOrDefault<word>
 318     (
 319         "interpolationScheme",
 320         "cell"
 321     );
 322     writeFormat_ = dict.lookupOrDefault<word>
 323     (
 324         "surfaceFormat",
 325         "null"
 326     );
 327
 328
 329     // define the generic (geometry) writer
 330     genericFormatter_ = surfaceWriter<bool>::New(writeFormat_);
 331
 332
 333     PtrList<sampledSurface> newList
 334     (
 335         dict.lookup("surfaces"),
 336         sampledSurface::iNew(mesh_)
 337     );
 338
 339     transfer(newList);
 340
 341     if (Pstream::parRun())
 342     {
 343         mergeList_.setSize(size());
 344     }
 345
 346     // ensure all surfaces and merge information are expired
 347     expire();
 348
 349     if (Pstream::master() && debug)
 350     {
 351         Pout<< "sample fields:" << fieldNames_ << nl
 352             << "sample surfaces:" << nl << "(" << nl;
 353
 354         forAll(*this, surfI)
 355         {
 356             Pout<< "  " << operator[](surfI) << endl;
 357         }
 358         Pout<< ")" << endl;
 359     }
 360 }
 361
 362
 363 void Foam::sampledSurfaces::updateMesh(const mapPolyMesh&)
 364 {
 365     expire();
 366 }
 367
 368
 369 void Foam::sampledSurfaces::movePoints(const pointField&)
 370 {
 371     expire();
 372 }
 373
 374
 375 void Foam::sampledSurfaces::readUpdate(const polyMesh::readUpdateState state)
 376 {
 377     if (state != polyMesh::UNCHANGED)
 378     {
 379         expire();
 380     }
 381 }
 382
 383
 384 bool Foam::sampledSurfaces::needsUpdate() const
 385 {
 386     forAll(*this, surfI)
 387     {
 388         if (operator[](surfI).needsUpdate())
 389         {
 390             return true;
 391         }
 392     }
 393
 394     return false;
 395 }
 396
 397
 398 bool Foam::sampledSurfaces::expire()
 399 {
 400     bool justExpired = false;
 401
 402     forAll(*this, surfI)
 403     {
 404         if (operator[](surfI).expire())
 405         {
 406             justExpired = true;
 407         }
 408
 409         // clear merge information
 410         if (Pstream::parRun())
 411         {
 412             mergeList_[surfI].clear();
 413         }
 414     }
 415
 416     // reset interpolation
 417     pointMesh::Delete(mesh_);
 418     volPointInterpolation::Delete(mesh_);
 419
 420     // true if any surfaces just expired
 421     return justExpired;
 422 }
 423
 424
 425 bool Foam::sampledSurfaces::update()
 426 {
 427     bool updated = false;
 428
 429     if (!needsUpdate())
 430     {
 431         return updated;
 432     }
 433
 434     // serial: quick and easy, no merging required
 435     if (!Pstream::parRun())
 436     {
 437         forAll(*this, surfI)
 438         {
 439             if (operator[](surfI).update())
 440             {
 441                 updated = true;
 442             }
 443         }
 444
 445         return updated;
 446     }
 447
 448     // dimension as fraction of mesh bounding box
 449     scalar mergeDim = mergeTol_ * mesh_.globalData().bb().mag();
 450
 451     if (Pstream::master() && debug)
 452     {
 453         Pout<< nl << "Merging all points within "
 454             << mergeDim << " meter" << endl;
 455     }
 456
 457     forAll(*this, surfI)
 458     {
 459         sampledSurface& s = operator[](surfI);
 460
 461         if (s.update())
 462         {
 463             updated = true;
 464         }
 465         else
 466         {
 467             continue;
 468         }
 469
 470
 471         // Collect points from all processors
 472         List<pointField> gatheredPoints(Pstream::nProcs());
 473         gatheredPoints[Pstream::myProcNo()] = s.points();
 474         Pstream::gatherList(gatheredPoints);
 475
 476         if (Pstream::master())
 477         {
 478             mergeList_[surfI].points = ListListOps::combine<pointField>
 479             (
 480                 gatheredPoints,
 481                 accessOp<pointField>()
 482             );
 483         }
 484
 485         // Collect faces from all processors and renumber using sizes of
 486         // gathered points
 487         List<faceList> gatheredFaces(Pstream::nProcs());
 488         gatheredFaces[Pstream::myProcNo()] = s.faces();
 489         Pstream::gatherList(gatheredFaces);
 490
 491         if (Pstream::master())
 492         {
 493             mergeList_[surfI].faces = static_cast<const faceList&>
 494             (
 495                 ListListOps::combineOffset<faceList>
 496                 (
 497                     gatheredFaces,
 498                     ListListOps::subSizes
 499                     (
 500                         gatheredPoints,
 501                         accessOp<pointField>()
 502                     ),
 503                     accessOp<faceList>(),
 504                     offsetOp<face>()
 505                 )
 506             );
 507         }
 508
 509         pointField newPoints;
 510         labelList oldToNew;
 511
 512         bool hasMerged = mergePoints
 513         (
 514             mergeList_[surfI].points,
 515             mergeDim,
 516             false,                  // verbosity
 517             oldToNew,
 518             newPoints
 519         );
 520
 521         if (hasMerged)
 522         {
 523             // Store point mapping
 524             mergeList_[surfI].pointsMap.transfer(oldToNew);
 525
 526             // Copy points
 527             mergeList_[surfI].points.transfer(newPoints);
 528
 529             // Relabel faces
 530             faceList& faces = mergeList_[surfI].faces;
 531
 532             forAll(faces, faceI)
 533             {
 534                 inplaceRenumber(mergeList_[surfI].pointsMap, faces[faceI]);
 535             }
 536
 537             if (Pstream::master() && debug)
 538             {
 539                 Pout<< "For surface " << surfI << " merged from "
 540                     << mergeList_[surfI].pointsMap.size() << " points down to "
 541                     << mergeList_[surfI].points.size()    << " points" << endl;
 542             }
 543         }
 544     }
 545
 546     return updated;
 547 }
 548
 549
 550 // ************************************************************************* //