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