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) 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 "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(average(fvm, *pointFieldPtr_).ptr());
 167             volFieldPtr_ = storedVolFieldPtr_.operator->();
 168         }
 169
 170
 171         if (debug)
 172         {
 173             Info<< "sampledIsoSurface::getIsoField() : volField "
 174                 << volFieldPtr_->name() << " min:" << min(*volFieldPtr_).value()
 175                 << " max:" << max(*volFieldPtr_).value() << endl;
 176             Info<< "sampledIsoSurface::getIsoField() : pointField "
 177                 << pointFieldPtr_->name()
 178                 << " min:" << gMin(pointFieldPtr_->internalField())
 179                 << " max:" << gMax(pointFieldPtr_->internalField()) << endl;
 180         }
 181     }
 182     else
 183     {
 184         // Get subMesh variants
 185         const fvMesh& subFvm = subMeshPtr_().subMesh();
 186
 187         // Either lookup on the submesh or subset the whole-mesh volField
 188
 189         if (subFvm.foundObject<volScalarField>(isoField_))
 190         {
 191             if (debug)
 192             {
 193                 Info<< "sampledIsoSurface::getIsoField() :"
 194                     << " submesh lookup volField "
 195                     << isoField_ << endl;
 196             }
 197             storedVolSubFieldPtr_.clear();
 198             volSubFieldPtr_ = &subFvm.lookupObject<volScalarField>(isoField_);
 199         }
 200         else
 201         {
 202             if (debug)
 203             {
 204                 Info<< "sampledIsoSurface::getIsoField() : subsetting volField "
 205                     << isoField_ << endl;
 206             }
 207             storedVolSubFieldPtr_.reset
 208             (
 209                 subMeshPtr_().interpolate
 210                 (
 211                     *volFieldPtr_
 212                 ).ptr()
 213             );
 214             storedVolSubFieldPtr_->checkOut();
 215             volSubFieldPtr_ = storedVolSubFieldPtr_.operator->();
 216         }
 217
 218
 219         // Pointfield on submesh
 220
 221         word pointFldName =
 222             "volPointInterpolate("
 223           + volSubFieldPtr_->name()
 224           + ')';
 225
 226         if (subFvm.foundObject<pointScalarField>(pointFldName))
 227         {
 228             if (debug)
 229             {
 230                 Info<< "sampledIsoSurface::getIsoField() :"
 231                     << " submesh lookup pointField " << pointFldName << endl;
 232             }
 233             storedPointSubFieldPtr_.clear();
 234             pointSubFieldPtr_ = &subFvm.lookupObject<pointScalarField>
 235             (
 236                 pointFldName
 237             );
 238         }
 239         else
 240         {
 241             if (debug)
 242             {
 243                 Info<< "sampledIsoSurface::getIsoField() :"
 244                     << " interpolating submesh volField "
 245                     << volSubFieldPtr_->name()
 246                     << " to get submesh pointField " << pointFldName << endl;
 247             }
 248             storedPointSubFieldPtr_.reset
 249             (
 250                 volPointInterpolation::New
 251                 (
 252                     subFvm
 253                 ).interpolate(*volSubFieldPtr_).ptr()
 254             );
 255             storedPointSubFieldPtr_->checkOut();
 256             pointSubFieldPtr_ = storedPointSubFieldPtr_.operator->();
 257         }
 258
 259
 260
 261         // If averaging redo the volField. Can only be done now since needs the
 262         // point field.
 263         if (average_)
 264         {
 265             storedVolSubFieldPtr_.reset
 266             (
 267                 average(subFvm, *pointSubFieldPtr_).ptr()
 268             );
 269             volSubFieldPtr_ = storedVolSubFieldPtr_.operator->();
 270         }
 271
 272
 273         if (debug)
 274         {
 275             Info<< "sampledIsoSurface::getIsoField() : volSubField "
 276                 << volSubFieldPtr_->name()
 277                 << " min:" << min(*volSubFieldPtr_).value()
 278                 << " max:" << max(*volSubFieldPtr_).value() << endl;
 279             Info<< "sampledIsoSurface::getIsoField() : pointSubField "
 280                 << pointSubFieldPtr_->name()
 281                 << " min:" << gMin(pointSubFieldPtr_->internalField())
 282                 << " max:" << gMax(pointSubFieldPtr_->internalField()) << endl;
 283         }
 284     }
 285 }
 286
 287
 288 Foam::tmp<Foam::volScalarField> Foam::sampledIsoSurface::average
 289 (
 290     const fvMesh& mesh,
 291     const pointScalarField& pfld
 292 ) const
 293 {
 294     tmp<volScalarField> tcellAvg
 295     (
 296         new volScalarField
 297         (
 298             IOobject
 299             (
 300                 "cellAvg",
 301                 mesh.time().timeName(),
 302                 mesh,
 303                 IOobject::NO_READ,
 304                 IOobject::NO_WRITE,
 305                 false
 306             ),
 307             mesh,
 308             dimensionedScalar("zero", dimless, scalar(0.0))
 309         )
 310     );
 311     volScalarField& cellAvg = tcellAvg();
 312
 313     labelField nPointCells(mesh.nCells(), 0);
 314     {
 315         for (label pointI = 0; pointI < mesh.nPoints(); pointI++)
 316         {
 317             const labelList& pCells = mesh.pointCells(pointI);
 318
 319             forAll(pCells, i)
 320             {
 321                 label cellI = pCells[i];
 322
 323                 cellAvg[cellI] += pfld[pointI];
 324                 nPointCells[cellI]++;
 325             }
 326         }
 327     }
 328     forAll(cellAvg, cellI)
 329     {
 330         cellAvg[cellI] /= nPointCells[cellI];
 331     }
 332     // Give value to calculatedFvPatchFields
 333     cellAvg.correctBoundaryConditions();
 334
 335     return tcellAvg;
 336 }
 337
 338
 339 Foam::tmp<Foam::pointScalarField> Foam::sampledIsoSurface::average
 340 (
 341     const pointMesh& pMesh,
 342     const volScalarField& fld
 343 ) const
 344 {
 345     tmp<pointScalarField> tpointAvg
 346     (
 347         new pointScalarField
 348         (
 349             IOobject
 350             (
 351                 "pointAvg",
 352                 fld.time().timeName(),
 353                 fld.db(),
 354                 IOobject::NO_READ,
 355                 IOobject::NO_WRITE,
 356                 false
 357             ),
 358             pMesh,
 359             dimensionedScalar("zero", dimless, scalar(0.0))
 360         )
 361     );
 362     pointScalarField& pointAvg = tpointAvg();
 363
 364     for (label pointI = 0; pointI < fld.mesh().nPoints(); pointI++)
 365     {
 366         const labelList& pCells = fld.mesh().pointCells(pointI);
 367
 368         forAll(pCells, i)
 369         {
 370             pointAvg[pointI] += fld[pCells[i]];
 371         }
 372         pointAvg[pointI] /= pCells.size();
 373     }
 374     // Give value to calculatedFvPatchFields
 375     pointAvg.correctBoundaryConditions();
 376
 377     return tpointAvg;
 378 }
 379
 380
 381 bool Foam::sampledIsoSurface::updateGeometry() const
 382 {
 383     const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
 384
 385     // no update needed
 386     if (fvm.time().timeIndex() == prevTimeIndex_)
 387     {
 388         return false;
 389     }
 390
 391     // Get any subMesh
 392     if (zoneID_.index() != -1 && !subMeshPtr_.valid())
 393     {
 394         const polyBoundaryMesh& patches = mesh().boundaryMesh();
 395
 396         // Patch to put exposed internal faces into
 397         label exposedPatchI = patches.findPatchID(exposedPatchName_);
 398
 399         if (debug)
 400         {
 401             Info<< "Allocating subset of size "
 402                 << mesh().cellZones()[zoneID_.index()].size()
 403                 << " with exposed faces into patch "
 404                 << patches[exposedPatchI].name() << endl;
 405         }
 406
 407         subMeshPtr_.reset
 408         (
 409             new fvMeshSubset(fvm)
 410         );
 411         subMeshPtr_().setLargeCellSubset
 412         (
 413             labelHashSet(mesh().cellZones()[zoneID_.index()]),
 414             exposedPatchI
 415         );
 416     }
 417
 418
 419     prevTimeIndex_ = fvm.time().timeIndex();
 420     getIsoFields();
 421
 422     // Clear any stored topo
 423     surfPtr_.clear();
 424     facesPtr_.clear();
 425
 426     // Clear derived data
 427     clearGeom();
 428
 429     if (subMeshPtr_.valid())
 430     {
 431         surfPtr_.reset
 432         (
 433             new isoSurface
 434             (
 435                 *volSubFieldPtr_,
 436                 *pointSubFieldPtr_,
 437                 isoVal_,
 438                 regularise_,
 439                 mergeTol_
 440             )
 441         );
 442     }
 443     else
 444     {
 445         surfPtr_.reset
 446         (
 447             new isoSurface
 448             (
 449                 *volFieldPtr_,
 450                 *pointFieldPtr_,
 451                 isoVal_,
 452                 regularise_,
 453                 mergeTol_
 454             )
 455         );
 456     }
 457
 458
 459     if (debug)
 460     {
 461         Pout<< "sampledIsoSurface::updateGeometry() : constructed iso:"
 462             << nl
 463             << "    regularise     : " << regularise_ << nl
 464             << "    average        : " << average_ << nl
 465             << "    isoField       : " << isoField_ << nl
 466             << "    isoValue       : " << isoVal_ << nl;
 467         if (subMeshPtr_.valid())
 468         {
 469             Pout<< "    zone size      : " << subMeshPtr_().subMesh().nCells()
 470                 << nl;
 471         }
 472         Pout<< "    points         : " << points().size() << nl
 473             << "    tris           : " << surface().size() << nl
 474             << "    cut cells      : " << surface().meshCells().size()
 475             << endl;
 476     }
 477
 478     return true;
 479 }
 480
 481
 482 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 483
 484 Foam::sampledIsoSurface::sampledIsoSurface
 485 (
 486     const word& name,
 487     const polyMesh& mesh,
 488     const dictionary& dict
 489 )
 490 :
 491     sampledSurface(name, mesh, dict),
 492     isoField_(dict.lookup("isoField")),
 493     isoVal_(readScalar(dict.lookup("isoValue"))),
 494     mergeTol_(dict.lookupOrDefault("mergeTol", 1E-6)),
 495     regularise_(dict.lookupOrDefault("regularise", true)),
 496     average_(dict.lookupOrDefault("average", false)),
 497     zoneID_(dict.lookupOrDefault("zone", word::null), mesh.cellZones()),
 498     exposedPatchName_(word::null),
 499     surfPtr_(NULL),
 500     facesPtr_(NULL),
 501     prevTimeIndex_(-1),
 502     storedVolFieldPtr_(NULL),
 503     volFieldPtr_(NULL),
 504     storedPointFieldPtr_(NULL),
 505     pointFieldPtr_(NULL)
 506 {
 507     if (!sampledSurface::interpolate())
 508     {
 509         FatalIOErrorIn
 510         (
 511             "sampledIsoSurface::sampledIsoSurface"
 512             "(const word&, const polyMesh&, const dictionary&)",
 513             dict
 514         )   << "Non-interpolated iso surface not supported since triangles"
 515             << " span across cells." << exit(FatalIOError);
 516     }
 517
 518     if (zoneID_.index() != -1)
 519     {
 520         dict.lookup("exposedPatchName") >> exposedPatchName_;
 521
 522         if (mesh.boundaryMesh().findPatchID(exposedPatchName_) == -1)
 523         {
 524             FatalIOErrorIn
 525             (
 526                 "sampledIsoSurface::sampledIsoSurface"
 527                 "(const word&, const polyMesh&, const dictionary&)",
 528                 dict
 529             )   << "Cannot find patch " << exposedPatchName_
 530                 << " in which to put exposed faces." << endl
 531                 << "Valid patches are " << mesh.boundaryMesh().names()
 532                 << exit(FatalIOError);
 533         }
 534
 535         if (debug && zoneID_.index() != -1)
 536         {
 537             Info<< "Restricting to cellZone " << zoneID_.name()
 538                 << " with exposed internal faces into patch "
 539                 << exposedPatchName_ << endl;
 540         }
 541     }
 542 }
 543
 544
 545 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 546
 547 Foam::sampledIsoSurface::~sampledIsoSurface()
 548 {}
 549
 550
 551 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 552
 553 bool Foam::sampledIsoSurface::needsUpdate() const
 554 {
 555     const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
 556
 557     return fvm.time().timeIndex() != prevTimeIndex_;
 558 }
 559
 560
 561 bool Foam::sampledIsoSurface::expire()
 562 {
 563     surfPtr_.clear();
 564     facesPtr_.clear();
 565     subMeshPtr_.clear();
 566
 567     // Clear derived data
 568     clearGeom();
 569
 570     // already marked as expired
 571     if (prevTimeIndex_ == -1)
 572     {
 573         return false;
 574     }
 575
 576     // force update
 577     prevTimeIndex_ = -1;
 578     return true;
 579 }
 580
 581
 582 bool Foam::sampledIsoSurface::update()
 583 {
 584     return updateGeometry();
 585 }
 586
 587
 588 Foam::tmp<Foam::scalarField> Foam::sampledIsoSurface::sample
 589 (
 590     const volScalarField& vField
 591 ) const
 592 {
 593     return sampleField(vField);
 594 }
 595
 596
 597 Foam::tmp<Foam::vectorField> Foam::sampledIsoSurface::sample
 598 (
 599     const volVectorField& vField
 600 ) const
 601 {
 602     return sampleField(vField);
 603 }
 604
 605
 606 Foam::tmp<Foam::sphericalTensorField> Foam::sampledIsoSurface::sample
 607 (
 608     const volSphericalTensorField& vField
 609 ) const
 610 {
 611     return sampleField(vField);
 612 }
 613
 614
 615 Foam::tmp<Foam::symmTensorField> Foam::sampledIsoSurface::sample
 616 (
 617     const volSymmTensorField& vField
 618 ) const
 619 {
 620     return sampleField(vField);
 621 }
 622
 623
 624 Foam::tmp<Foam::tensorField> Foam::sampledIsoSurface::sample
 625 (
 626     const volTensorField& vField
 627 ) const
 628 {
 629     return sampleField(vField);
 630 }
 631
 632
 633 Foam::tmp<Foam::scalarField> Foam::sampledIsoSurface::interpolate
 634 (
 635     const interpolation<scalar>& interpolator
 636 ) const
 637 {
 638     return interpolateField(interpolator);
 639 }
 640
 641
 642 Foam::tmp<Foam::vectorField> Foam::sampledIsoSurface::interpolate
 643 (
 644     const interpolation<vector>& interpolator
 645 ) const
 646 {
 647     return interpolateField(interpolator);
 648 }
 649
 650 Foam::tmp<Foam::sphericalTensorField> Foam::sampledIsoSurface::interpolate
 651 (
 652     const interpolation<sphericalTensor>& interpolator
 653 ) const
 654 {
 655     return interpolateField(interpolator);
 656 }
 657
 658
 659 Foam::tmp<Foam::symmTensorField> Foam::sampledIsoSurface::interpolate
 660 (
 661     const interpolation<symmTensor>& interpolator
 662 ) const
 663 {
 664     return interpolateField(interpolator);
 665 }
 666
 667
 668 Foam::tmp<Foam::tensorField> Foam::sampledIsoSurface::interpolate
 669 (
 670     const interpolation<tensor>& interpolator
 671 ) const
 672 {
 673     return interpolateField(interpolator);
 674 }
 675
 676
 677 void Foam::sampledIsoSurface::print(Ostream& os) const
 678 {
 679     os  << "sampledIsoSurface: " << name() << " :"
 680         << "  field   :" << isoField_
 681         << "  value   :" << isoVal_;
 682         //<< "  faces:" << faces().size()       // note: possibly no geom yet
 683         //<< "  points:" << points().size();
 684 }
 685
 686
 687 // ************************************************************************* //