src/sampling/sampledSurface/distanceSurface/distanceSurface.C

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