src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   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 "sampledIsoSurface.H"
  27 #include "dictionary.H"
  28 #include "volFields.H"
  29 #include "volPointInterpolation.H"
  30 #include "addToRunTimeSelectionTable.H"
  31
  32 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  33
  34 namespace Foam
  35 {
  36     defineTypeNameAndDebug(sampledIsoSurface, 0);
  37     addNamedToRunTimeSelectionTable
  38     (
  39         sampledSurface,
  40         sampledIsoSurface,
  41         word,
  42         isoSurface
  43     );
  44 }
  45
  46 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  47
  48 void Foam::sampledIsoSurface::getIsoFields() const
  49 {
  50     const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
  51
  52     // Get volField
  53     // ~~~~~~~~~~~~
  54
  55     if (fvm.foundObject<volScalarField>(isoField_))
  56     {
  57         if (debug)
  58         {
  59             Info<< "sampledIsoSurface::getIsoField() : lookup volField "
  60                 << isoField_ << endl;
  61         }
  62         storedVolFieldPtr_.clear();
  63         volFieldPtr_ = &fvm.lookupObject<volScalarField>(isoField_);
  64     }
  65     else
  66     {
  67         // Bit of a hack. Read field and store.
  68
  69         if (debug)
  70         {
  71             Info<< "sampledIsoSurface::getIsoField() : checking "
  72                 << isoField_ << " for same time " << fvm.time().timeName()
  73                 << endl;
  74         }
  75
  76         if
  77         (
  78             storedVolFieldPtr_.empty()
  79          || (fvm.time().timeName() != storedVolFieldPtr_().instance())
  80         )
  81         {
  82             if (debug)
  83             {
  84                 Info<< "sampledIsoSurface::getIsoField() : reading volField "
  85                     << isoField_ << " from time " << fvm.time().timeName()
  86                     << endl;
  87             }
  88
  89             storedVolFieldPtr_.reset
  90             (
  91                 new volScalarField
  92                 (
  93                     IOobject
  94                     (
  95                         isoField_,
  96                         fvm.time().timeName(),
  97                         fvm,
  98                         IOobject::MUST_READ,
  99                         IOobject::NO_WRITE,
 100                         false
 101                     ),
 102                     fvm
 103                 )
 104             );
 105             volFieldPtr_ = storedVolFieldPtr_.operator->();
 106         }
 107     }
 108
 109
 110
 111     // Get pointField
 112     // ~~~~~~~~~~~~~~
 113
 114     if (!subMeshPtr_.valid())
 115     {
 116         word pointFldName = "volPointInterpolate(" + isoField_ + ')';
 117
 118         if (fvm.foundObject<pointScalarField>(pointFldName))
 119         {
 120             if (debug)
 121             {
 122                 Info<< "sampledIsoSurface::getIsoField() : lookup pointField "
 123                     << pointFldName << endl;
 124             }
 125             pointFieldPtr_ = &fvm.lookupObject<pointScalarField>(pointFldName);
 126         }
 127         else
 128         {
 129             // Not in registry. Interpolate.
 130
 131             if (debug)
 132             {
 133                 Info<< "sampledIsoSurface::getIsoField() : checking pointField "
 134                     << pointFldName << " for same time "
 135                     << fvm.time().timeName() << endl;
 136             }
 137
 138             if
 139             (
 140                 storedPointFieldPtr_.empty()
 141              || (fvm.time().timeName() != storedPointFieldPtr_().instance())
 142             )
 143             {
 144                 if (debug)
 145                 {
 146                     Info<< "sampledIsoSurface::getIsoField() :"
 147                         << " interpolating volField " << volFieldPtr_->name()
 148                         << " to get pointField " << pointFldName << endl;
 149                 }
 150
 151                 storedPointFieldPtr_.reset
 152                 (
 153                     volPointInterpolation::New(fvm)
 154                     .interpolate(*volFieldPtr_).ptr()
 155                 );
 156                 storedPointFieldPtr_->checkOut();
 157                 pointFieldPtr_ = storedPointFieldPtr_.operator->();
 158             }
 159         }
 160
 161
 162         // If averaging redo the volField. Can only be done now since needs the
 163         // point field.
 164         if (average_)
 165         {
 166             storedVolFieldPtr_.reset
 167             (
 168                 pointAverage(*pointFieldPtr_).ptr()
 169             );
 170             volFieldPtr_ = storedVolFieldPtr_.operator->();
 171         }
 172
 173
 174         if (debug)
 175         {
 176             Info<< "sampledIsoSurface::getIsoField() : volField "
 177                 << volFieldPtr_->name() << " min:" << min(*volFieldPtr_).value()
 178                 << " max:" << max(*volFieldPtr_).value() << endl;
 179             Info<< "sampledIsoSurface::getIsoField() : pointField "
 180                 << pointFieldPtr_->name()
 181                 << " min:" << gMin(pointFieldPtr_->internalField())
 182                 << " max:" << gMax(pointFieldPtr_->internalField()) << endl;
 183         }
 184     }
 185     else
 186     {
 187         // Get subMesh variants
 188         const fvMesh& subFvm = subMeshPtr_().subMesh();
 189
 190         // Either lookup on the submesh or subset the whole-mesh volField
 191
 192         if (subFvm.foundObject<volScalarField>(isoField_))
 193         {
 194             if (debug)
 195             {
 196                 Info<< "sampledIsoSurface::getIsoField() :"
 197                     << " submesh lookup volField "
 198                     << isoField_ << endl;
 199             }
 200             storedVolSubFieldPtr_.clear();
 201             volSubFieldPtr_ = &subFvm.lookupObject<volScalarField>(isoField_);
 202         }
 203         else
 204         {
 205             if (debug)
 206             {
 207                 Info<< "sampledIsoSurface::getIsoField() : subsetting volField "
 208                     << isoField_ << endl;
 209             }
 210             storedVolSubFieldPtr_.reset
 211             (
 212                 subMeshPtr_().interpolate
 213                 (
 214                     *volFieldPtr_
 215                 ).ptr()
 216             );
 217             storedVolSubFieldPtr_->checkOut();
 218             volSubFieldPtr_ = storedVolSubFieldPtr_.operator->();
 219         }
 220
 221
 222         // Pointfield on submesh
 223
 224         word pointFldName =
 225             "volPointInterpolate("
 226           + volSubFieldPtr_->name()
 227           + ')';
 228
 229         if (subFvm.foundObject<pointScalarField>(pointFldName))
 230         {
 231             if (debug)
 232             {
 233                 Info<< "sampledIsoSurface::getIsoField() :"
 234                     << " submesh lookup pointField " << pointFldName << endl;
 235             }
 236             storedPointSubFieldPtr_.clear();
 237             pointSubFieldPtr_ = &subFvm.lookupObject<pointScalarField>
 238             (
 239                 pointFldName
 240             );
 241         }
 242         else
 243         {
 244             if (debug)
 245             {
 246                 Info<< "sampledIsoSurface::getIsoField() :"
 247                     << " interpolating submesh volField "
 248                     << volSubFieldPtr_->name()
 249                     << " to get submesh pointField " << pointFldName << endl;
 250             }
 251             storedPointSubFieldPtr_.reset
 252             (
 253                 volPointInterpolation::New
 254                 (
 255                     subFvm
 256                 ).interpolate(*volSubFieldPtr_).ptr()
 257             );
 258             storedPointSubFieldPtr_->checkOut();
 259             pointSubFieldPtr_ = storedPointSubFieldPtr_.operator->();
 260         }
 261
 262
 263
 264         // If averaging redo the volField. Can only be done now since needs the
 265         // point field.
 266         if (average_)
 267         {
 268             storedVolSubFieldPtr_.reset
 269             (
 270                 pointAverage(*pointSubFieldPtr_).ptr()
 271             );
 272             volSubFieldPtr_ = storedVolSubFieldPtr_.operator->();
 273         }
 274
 275
 276         if (debug)
 277         {
 278             Info<< "sampledIsoSurface::getIsoField() : volSubField "
 279                 << volSubFieldPtr_->name()
 280                 << " min:" << min(*volSubFieldPtr_).value()
 281                 << " max:" << max(*volSubFieldPtr_).value() << endl;
 282             Info<< "sampledIsoSurface::getIsoField() : pointSubField "
 283                 << pointSubFieldPtr_->name()
 284                 << " min:" << gMin(pointSubFieldPtr_->internalField())
 285                 << " max:" << gMax(pointSubFieldPtr_->internalField()) << endl;
 286         }
 287     }
 288 }
 289
 290
 291 bool Foam::sampledIsoSurface::updateGeometry() const
 292 {
 293     const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
 294
 295     // no update needed
 296     if (fvm.time().timeIndex() == prevTimeIndex_)
 297     {
 298         return false;
 299     }
 300
 301     // Get any subMesh
 302     if (zoneID_.index() != -1 && !subMeshPtr_.valid())
 303     {
 304         const polyBoundaryMesh& patches = mesh().boundaryMesh();
 305
 306         // Patch to put exposed internal faces into
 307         const label exposedPatchI = patches.findPatchID(exposedPatchName_);
 308
 309         if (debug)
 310         {
 311             Info<< "Allocating subset of size "
 312                 << mesh().cellZones()[zoneID_.index()].size()
 313                 << " with exposed faces into patch "
 314                 << patches[exposedPatchI].name() << endl;
 315         }
 316
 317         subMeshPtr_.reset
 318         (
 319             new fvMeshSubset(fvm)
 320         );
 321         subMeshPtr_().setLargeCellSubset
 322         (
 323             labelHashSet(mesh().cellZones()[zoneID_.index()]),
 324             exposedPatchI
 325         );
 326     }
 327
 328
 329     prevTimeIndex_ = fvm.time().timeIndex();
 330     getIsoFields();
 331
 332     // Clear any stored topo
 333     surfPtr_.clear();
 334     facesPtr_.clear();
 335
 336     // Clear derived data
 337     clearGeom();
 338
 339     if (subMeshPtr_.valid())
 340     {
 341         surfPtr_.reset
 342         (
 343             new isoSurface
 344             (
 345                 *volSubFieldPtr_,
 346                 *pointSubFieldPtr_,
 347                 isoVal_,
 348                 regularise_,
 349                 mergeTol_
 350             )
 351         );
 352     }
 353     else
 354     {
 355         surfPtr_.reset
 356         (
 357             new isoSurface
 358             (
 359                 *volFieldPtr_,
 360                 *pointFieldPtr_,
 361                 isoVal_,
 362                 regularise_,
 363                 mergeTol_
 364             )
 365         );
 366     }
 367
 368
 369     if (debug)
 370     {
 371         Pout<< "sampledIsoSurface::updateGeometry() : constructed iso:"
 372             << nl
 373             << "    regularise     : " << regularise_ << nl
 374             << "    average        : " << average_ << nl
 375             << "    isoField       : " << isoField_ << nl
 376             << "    isoValue       : " << isoVal_ << nl;
 377         if (subMeshPtr_.valid())
 378         {
 379             Pout<< "    zone size      : " << subMeshPtr_().subMesh().nCells()
 380                 << nl;
 381         }
 382         Pout<< "    points         : " << points().size() << nl
 383             << "    tris           : " << surface().size() << nl
 384             << "    cut cells      : " << surface().meshCells().size()
 385             << endl;
 386     }
 387
 388     return true;
 389 }
 390
 391
 392 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 393
 394 Foam::sampledIsoSurface::sampledIsoSurface
 395 (
 396     const word& name,
 397     const polyMesh& mesh,
 398     const dictionary& dict
 399 )
 400 :
 401     sampledSurface(name, mesh, dict),
 402     isoField_(dict.lookup("isoField")),
 403     isoVal_(readScalar(dict.lookup("isoValue"))),
 404     mergeTol_(dict.lookupOrDefault("mergeTol", 1E-6)),
 405     regularise_(dict.lookupOrDefault("regularise", true)),
 406     average_(dict.lookupOrDefault("average", false)),
 407     zoneID_(dict.lookupOrDefault("zone", word::null), mesh.cellZones()),
 408     exposedPatchName_(word::null),
 409     surfPtr_(NULL),
 410     facesPtr_(NULL),
 411     prevTimeIndex_(-1),
 412     storedVolFieldPtr_(NULL),
 413     volFieldPtr_(NULL),
 414     storedPointFieldPtr_(NULL),
 415     pointFieldPtr_(NULL)
 416 {
 417     if (!sampledSurface::interpolate())
 418     {
 419         FatalIOErrorIn
 420         (
 421             "sampledIsoSurface::sampledIsoSurface"
 422             "(const word&, const polyMesh&, const dictionary&)",
 423             dict
 424         )   << "Non-interpolated iso surface not supported since triangles"
 425             << " span across cells." << exit(FatalIOError);
 426     }
 427
 428     if (zoneID_.index() != -1)
 429     {
 430         dict.lookup("exposedPatchName") >> exposedPatchName_;
 431
 432         if (mesh.boundaryMesh().findPatchID(exposedPatchName_) == -1)
 433         {
 434             FatalIOErrorIn
 435             (
 436                 "sampledIsoSurface::sampledIsoSurface"
 437                 "(const word&, const polyMesh&, const dictionary&)",
 438                 dict
 439             )   << "Cannot find patch " << exposedPatchName_
 440                 << " in which to put exposed faces." << endl
 441                 << "Valid patches are " << mesh.boundaryMesh().names()
 442                 << exit(FatalIOError);
 443         }
 444
 445         if (debug && zoneID_.index() != -1)
 446         {
 447             Info<< "Restricting to cellZone " << zoneID_.name()
 448                 << " with exposed internal faces into patch "
 449                 << exposedPatchName_ << endl;
 450         }
 451     }
 452 }
 453
 454
 455 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 456
 457 Foam::sampledIsoSurface::~sampledIsoSurface()
 458 {}
 459
 460
 461 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 462
 463 bool Foam::sampledIsoSurface::needsUpdate() const
 464 {
 465     const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
 466
 467     return fvm.time().timeIndex() != prevTimeIndex_;
 468 }
 469
 470
 471 bool Foam::sampledIsoSurface::expire()
 472 {
 473     surfPtr_.clear();
 474     facesPtr_.clear();
 475     subMeshPtr_.clear();
 476
 477     // Clear derived data
 478     clearGeom();
 479
 480     // already marked as expired
 481     if (prevTimeIndex_ == -1)
 482     {
 483         return false;
 484     }
 485
 486     // force update
 487     prevTimeIndex_ = -1;
 488     return true;
 489 }
 490
 491
 492 bool Foam::sampledIsoSurface::update()
 493 {
 494     return updateGeometry();
 495 }
 496
 497
 498 Foam::tmp<Foam::scalarField> Foam::sampledIsoSurface::sample
 499 (
 500     const volScalarField& vField
 501 ) const
 502 {
 503     return sampleField(vField);
 504 }
 505
 506
 507 Foam::tmp<Foam::vectorField> Foam::sampledIsoSurface::sample
 508 (
 509     const volVectorField& vField
 510 ) const
 511 {
 512     return sampleField(vField);
 513 }
 514
 515
 516 Foam::tmp<Foam::sphericalTensorField> Foam::sampledIsoSurface::sample
 517 (
 518     const volSphericalTensorField& vField
 519 ) const
 520 {
 521     return sampleField(vField);
 522 }
 523
 524
 525 Foam::tmp<Foam::symmTensorField> Foam::sampledIsoSurface::sample
 526 (
 527     const volSymmTensorField& vField
 528 ) const
 529 {
 530     return sampleField(vField);
 531 }
 532
 533
 534 Foam::tmp<Foam::tensorField> Foam::sampledIsoSurface::sample
 535 (
 536     const volTensorField& vField
 537 ) const
 538 {
 539     return sampleField(vField);
 540 }
 541
 542
 543 Foam::tmp<Foam::scalarField> Foam::sampledIsoSurface::interpolate
 544 (
 545     const interpolation<scalar>& interpolator
 546 ) const
 547 {
 548     return interpolateField(interpolator);
 549 }
 550
 551
 552 Foam::tmp<Foam::vectorField> Foam::sampledIsoSurface::interpolate
 553 (
 554     const interpolation<vector>& interpolator
 555 ) const
 556 {
 557     return interpolateField(interpolator);
 558 }
 559
 560 Foam::tmp<Foam::sphericalTensorField> Foam::sampledIsoSurface::interpolate
 561 (
 562     const interpolation<sphericalTensor>& interpolator
 563 ) const
 564 {
 565     return interpolateField(interpolator);
 566 }
 567
 568
 569 Foam::tmp<Foam::symmTensorField> Foam::sampledIsoSurface::interpolate
 570 (
 571     const interpolation<symmTensor>& interpolator
 572 ) const
 573 {
 574     return interpolateField(interpolator);
 575 }
 576
 577
 578 Foam::tmp<Foam::tensorField> Foam::sampledIsoSurface::interpolate
 579 (
 580     const interpolation<tensor>& interpolator
 581 ) const
 582 {
 583     return interpolateField(interpolator);
 584 }
 585
 586
 587 void Foam::sampledIsoSurface::print(Ostream& os) const
 588 {
 589     os  << "sampledIsoSurface: " << name() << " :"
 590         << "  field   :" << isoField_
 591         << "  value   :" << isoVal_;
 592         //<< "  faces:" << faces().size()       // note: possibly no geom yet
 593         //<< "  points:" << points().size();
 594 }
 595
 596
 597 // ************************************************************************* //