src/sampling/sampledSurface/distanceSurface/distanceSurface.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 "distanceSurface.H"
  27 #include "dictionary.H"
  28 #include "volFields.H"
  29 #include "volPointInterpolation.H"
  30 #include "addToRunTimeSelectionTable.H"
  31 #include "fvMesh.H"
  32
  33 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  34
  35 namespace Foam
  36 {
  37     defineTypeNameAndDebug(distanceSurface, 0);
  38     addToRunTimeSelectionTable(sampledSurface, distanceSurface, word);
  39 }
  40
  41 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  42
  43 void Foam::distanceSurface::createGeometry()
  44 {
  45     if (debug)
  46     {
  47         Pout<< "distanceSurface::createGeometry :updating geometry." << endl;
  48     }
  49
  50     // Clear any stored topologies
  51     facesPtr_.clear();
  52
  53     // Clear derived data
  54     clearGeom();
  55
  56     const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
  57
  58     // Distance to cell centres
  59     // ~~~~~~~~~~~~~~~~~~~~~~~~
  60
  61     cellDistancePtr_.reset
  62     (
  63         new volScalarField
  64         (
  65             IOobject
  66             (
  67                 "cellDistance",
  68                 fvm.time().timeName(),
  69                 fvm.time(),
  70                 IOobject::NO_READ,
  71                 IOobject::NO_WRITE,
  72                 false
  73             ),
  74             fvm,
  75             dimensionedScalar("zero", dimLength, 0)
  76         )
  77     );
  78     volScalarField& cellDistance = cellDistancePtr_();
  79
  80     // Internal field
  81     {
  82         const pointField& cc = fvm.C();
  83         scalarField& fld = cellDistance.internalField();
  84
  85         List<pointIndexHit> nearest;
  86         surfPtr_().findNearest
  87         (
  88             cc,
  89             scalarField(cc.size(), GREAT),
  90             nearest
  91         );
  92
  93         if (signed_)
  94         {
  95             List<searchableSurface::volumeType> volType;
  96
  97             surfPtr_().getVolumeType(cc, volType);
  98
  99             forAll(volType, i)
 100             {
 101                 searchableSurface::volumeType vT = volType[i];
 102
 103                 if (vT == searchableSurface::OUTSIDE)
 104                 {
 105                     fld[i] = Foam::mag(cc[i] - nearest[i].hitPoint());
 106                 }
 107                 else if (vT == searchableSurface::INSIDE)
 108                 {
 109                     fld[i] = -Foam::mag(cc[i] - nearest[i].hitPoint());
 110                 }
 111                 else
 112                 {
 113                     FatalErrorIn
 114                     (
 115                         "void Foam::distanceSurface::createGeometry()"
 116                     )   << "getVolumeType failure, neither INSIDE or OUTSIDE"
 117                         << exit(FatalError);
 118                 }
 119             }
 120         }
 121         else
 122         {
 123             forAll(nearest, i)
 124             {
 125                 fld[i] = Foam::mag(cc[i] - nearest[i].hitPoint());
 126             }
 127         }
 128     }
 129
 130     // Patch fields
 131     {
 132         forAll(fvm.C().boundaryField(), patchI)
 133         {
 134             const pointField& cc = fvm.C().boundaryField()[patchI];
 135             fvPatchScalarField& fld = cellDistance.boundaryField()[patchI];
 136
 137             List<pointIndexHit> nearest;
 138             surfPtr_().findNearest
 139             (
 140                 cc,
 141                 scalarField(cc.size(), GREAT),
 142                 nearest
 143             );
 144
 145             if (signed_)
 146             {
 147                 List<searchableSurface::volumeType> volType;
 148
 149                 surfPtr_().getVolumeType(cc, volType);
 150
 151                 forAll(volType, i)
 152                 {
 153                     searchableSurface::volumeType vT = volType[i];
 154
 155                     if (vT == searchableSurface::OUTSIDE)
 156                     {
 157                         fld[i] = Foam::mag(cc[i] - nearest[i].hitPoint());
 158                     }
 159                     else if (vT == searchableSurface::INSIDE)
 160                     {
 161                         fld[i] = -Foam::mag(cc[i] - nearest[i].hitPoint());
 162                     }
 163                     else
 164                     {
 165                         FatalErrorIn
 166                         (
 167                             "void Foam::distanceSurface::createGeometry()"
 168                         )   << "getVolumeType failure, "
 169                             << "neither INSIDE or OUTSIDE"
 170                             << exit(FatalError);
 171                     }
 172                 }
 173             }
 174             else
 175             {
 176                 forAll(nearest, i)
 177                 {
 178                     fld[i] = Foam::mag(cc[i] - nearest[i].hitPoint());
 179                 }
 180             }
 181         }
 182     }
 183
 184
 185     // On processor patches the mesh.C() will already be the cell centre
 186     // on the opposite side so no need to swap cellDistance.
 187
 188
 189     // Distance to points
 190     pointDistance_.setSize(fvm.nPoints());
 191     {
 192         const pointField& pts = fvm.points();
 193
 194         List<pointIndexHit> nearest;
 195         surfPtr_().findNearest
 196         (
 197             pts,
 198             scalarField(pts.size(), GREAT),
 199             nearest
 200         );
 201
 202         if (signed_)
 203         {
 204             List<searchableSurface::volumeType> volType;
 205
 206             surfPtr_().getVolumeType(pts, volType);
 207
 208             forAll(volType, i)
 209             {
 210                 searchableSurface::volumeType vT = volType[i];
 211
 212                 if (vT == searchableSurface::OUTSIDE)
 213                 {
 214                     pointDistance_[i] =
 215                         Foam::mag(pts[i] - nearest[i].hitPoint());
 216                 }
 217                 else if (vT == searchableSurface::INSIDE)
 218                 {
 219                     pointDistance_[i] =
 220                         -Foam::mag(pts[i] - nearest[i].hitPoint());
 221                 }
 222                 else
 223                 {
 224                     FatalErrorIn
 225                     (
 226                         "void Foam::distanceSurface::createGeometry()"
 227                     )   << "getVolumeType failure, neither INSIDE or OUTSIDE"
 228                         << exit(FatalError);
 229                 }
 230             }
 231         }
 232         else
 233         {
 234             forAll(nearest, i)
 235             {
 236                 pointDistance_[i] = Foam::mag(pts[i]-nearest[i].hitPoint());
 237             }
 238         }
 239     }
 240
 241
 242     if (debug)
 243     {
 244         Pout<< "Writing cell distance:" << cellDistance.objectPath() << endl;
 245         cellDistance.write();
 246         pointScalarField pDist
 247         (
 248             IOobject
 249             (
 250                 "pointDistance",
 251                 fvm.time().timeName(),
 252                 fvm.time(),
 253                 IOobject::NO_READ,
 254                 IOobject::NO_WRITE,
 255                 false
 256             ),
 257             pointMesh::New(fvm),
 258             dimensionedScalar("zero", dimLength, 0)
 259         );
 260         pDist.internalField() = pointDistance_;
 261
 262         Pout<< "Writing point distance:" << pDist.objectPath() << endl;
 263         pDist.write();
 264     }
 265
 266
 267     //- Direct from cell field and point field.
 268     isoSurfPtr_.reset
 269     (
 270         new isoSurface
 271         (
 272             cellDistance,
 273             pointDistance_,
 274             distance_,
 275             regularise_
 276         )
 277         //new isoSurfaceCell
 278         //(
 279         //    fvm,
 280         //    cellDistance,
 281         //    pointDistance_,
 282         //    distance_,
 283         //    regularise_
 284         //)
 285     );
 286
 287     if (debug)
 288     {
 289         print(Pout);
 290         Pout<< endl;
 291     }
 292 }
 293
 294
 295 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 296
 297 Foam::distanceSurface::distanceSurface
 298 (
 299     const word& name,
 300     const polyMesh& mesh,
 301     const dictionary& dict
 302 )
 303 :
 304     sampledSurface(name, mesh, dict),
 305     surfPtr_
 306     (
 307         searchableSurface::New
 308         (
 309             dict.lookup("surfaceType"),
 310             IOobject
 311             (
 312                 dict.lookupOrDefault("surfaceName", name),  // name
 313                 mesh.time().constant(),                     // directory
 314                 "triSurface",                               // instance
 315                 mesh.time(),                                // registry
 316                 IOobject::MUST_READ,
 317                 IOobject::NO_WRITE
 318             ),
 319             dict
 320         )
 321     ),
 322     distance_(readScalar(dict.lookup("distance"))),
 323     signed_(readBool(dict.lookup("signed"))),
 324     regularise_(dict.lookupOrDefault("regularise", true)),
 325     average_(dict.lookupOrDefault("average", false)),
 326     zoneKey_(keyType::null),
 327     needsUpdate_(true),
 328     isoSurfPtr_(NULL),
 329     facesPtr_(NULL)
 330 {
 331 //    dict.readIfPresent("zone", zoneKey_);
 332 //
 333 //    if (debug && zoneKey_.size() && mesh.cellZones().findZoneID(zoneKey_) < 0)
 334 //    {
 335 //        Info<< "cellZone " << zoneKey_
 336 //            << " not found - using entire mesh" << endl;
 337 //    }
 338 }
 339
 340
 341 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 342
 343 Foam::distanceSurface::~distanceSurface()
 344 {}
 345
 346
 347 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 348
 349 bool Foam::distanceSurface::needsUpdate() const
 350 {
 351     return needsUpdate_;
 352 }
 353
 354
 355 bool Foam::distanceSurface::expire()
 356 {
 357     if (debug)
 358     {
 359         Pout<< "distanceSurface::expire :"
 360             << " have-facesPtr_:" << facesPtr_.valid()
 361             << " needsUpdate_:" << needsUpdate_ << endl;
 362     }
 363
 364     // Clear any stored topologies
 365     facesPtr_.clear();
 366
 367     // Clear derived data
 368     clearGeom();
 369
 370     // already marked as expired
 371     if (needsUpdate_)
 372     {
 373         return false;
 374     }
 375
 376     needsUpdate_ = true;
 377     return true;
 378 }
 379
 380
 381 bool Foam::distanceSurface::update()
 382 {
 383     if (debug)
 384     {
 385         Pout<< "distanceSurface::update :"
 386             << " have-facesPtr_:" << facesPtr_.valid()
 387             << " needsUpdate_:" << needsUpdate_ << endl;
 388     }
 389
 390     if (!needsUpdate_)
 391     {
 392         return false;
 393     }
 394
 395     createGeometry();
 396
 397     needsUpdate_ = false;
 398     return true;
 399 }
 400
 401
 402 Foam::tmp<Foam::scalarField> Foam::distanceSurface::sample
 403 (
 404     const volScalarField& vField
 405 ) const
 406 {
 407     return sampleField(vField);
 408 }
 409
 410
 411 Foam::tmp<Foam::vectorField> Foam::distanceSurface::sample
 412 (
 413     const volVectorField& vField
 414 ) const
 415 {
 416     return sampleField(vField);
 417 }
 418
 419
 420 Foam::tmp<Foam::sphericalTensorField> Foam::distanceSurface::sample
 421 (
 422     const volSphericalTensorField& vField
 423 ) const
 424 {
 425     return sampleField(vField);
 426 }
 427
 428
 429 Foam::tmp<Foam::symmTensorField> Foam::distanceSurface::sample
 430 (
 431     const volSymmTensorField& vField
 432 ) const
 433 {
 434     return sampleField(vField);
 435 }
 436
 437
 438 Foam::tmp<Foam::tensorField> Foam::distanceSurface::sample
 439 (
 440     const volTensorField& vField
 441 ) const
 442 {
 443     return sampleField(vField);
 444 }
 445
 446
 447 Foam::tmp<Foam::scalarField> Foam::distanceSurface::interpolate
 448 (
 449     const interpolation<scalar>& interpolator
 450 ) const
 451 {
 452     return interpolateField(interpolator);
 453 }
 454
 455
 456 Foam::tmp<Foam::vectorField> Foam::distanceSurface::interpolate
 457 (
 458     const interpolation<vector>& interpolator
 459 ) const
 460 {
 461     return interpolateField(interpolator);
 462 }
 463
 464 Foam::tmp<Foam::sphericalTensorField> Foam::distanceSurface::interpolate
 465 (
 466     const interpolation<sphericalTensor>& interpolator
 467 ) const
 468 {
 469     return interpolateField(interpolator);
 470 }
 471
 472
 473 Foam::tmp<Foam::symmTensorField> Foam::distanceSurface::interpolate
 474 (
 475     const interpolation<symmTensor>& interpolator
 476 ) const
 477 {
 478     return interpolateField(interpolator);
 479 }
 480
 481
 482 Foam::tmp<Foam::tensorField> Foam::distanceSurface::interpolate
 483 (
 484     const interpolation<tensor>& interpolator
 485 ) const
 486 {
 487     return interpolateField(interpolator);
 488 }
 489
 490
 491 void Foam::distanceSurface::print(Ostream& os) const
 492 {
 493     os  << "distanceSurface: " << name() << " :"
 494         << "  surface:" << surfPtr_().name()
 495         << "  distance:" << distance_
 496         << "  faces:" << faces().size()
 497         << "  points:" << points().size();
 498 }
 499
 500
 501 // ************************************************************************* //