src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.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 "sampledCuttingPlane.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(sampledCuttingPlane, 0);
  38     addNamedToRunTimeSelectionTable
  39     (
  40         sampledSurface,
  41         sampledCuttingPlane,
  42         word,
  43         cuttingPlane
  44     );
  45 }
  46
  47 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  48
  49 void Foam::sampledCuttingPlane::createGeometry()
  50 {
  51     if (debug)
  52     {
  53         Pout<< "sampledCuttingPlane::createGeometry :updating geometry."
  54             << endl;
  55     }
  56
  57     // Clear any stored topologies
  58     facesPtr_.clear();
  59     isoSurfPtr_.ptr();
  60     pointDistance_.clear();
  61     cellDistancePtr_.clear();
  62
  63     // Clear derived data
  64     clearGeom();
  65
  66     // Get any subMesh
  67     if (zoneID_.index() != -1 && !subMeshPtr_.valid())
  68     {
  69         const polyBoundaryMesh& patches = mesh().boundaryMesh();
  70
  71         // Patch to put exposed internal faces into
  72         const label exposedPatchI = patches.findPatchID(exposedPatchName_);
  73
  74         if (debug)
  75         {
  76             Info<< "Allocating subset of size "
  77                 << mesh().cellZones()[zoneID_.index()].size()
  78                 << " with exposed faces into patch "
  79                 << patches[exposedPatchI].name() << endl;
  80         }
  81
  82         subMeshPtr_.reset
  83         (
  84             new fvMeshSubset(static_cast<const fvMesh&>(mesh()))
  85         );
  86         subMeshPtr_().setLargeCellSubset
  87         (
  88             labelHashSet(mesh().cellZones()[zoneID_.index()]),
  89             exposedPatchI
  90         );
  91     }
  92
  93
  94     // Select either the submesh or the underlying mesh
  95     const fvMesh& fvm =
  96     (
  97         subMeshPtr_.valid()
  98       ? subMeshPtr_().subMesh()
  99       : static_cast<const fvMesh&>(mesh())
 100     );
 101
 102
 103     // Distance to cell centres
 104     // ~~~~~~~~~~~~~~~~~~~~~~~~
 105
 106     cellDistancePtr_.reset
 107     (
 108         new volScalarField
 109         (
 110             IOobject
 111             (
 112                 "cellDistance",
 113                 fvm.time().timeName(),
 114                 fvm.time(),
 115                 IOobject::NO_READ,
 116                 IOobject::NO_WRITE,
 117                 false
 118             ),
 119             fvm,
 120             dimensionedScalar("zero", dimLength, 0)
 121         )
 122     );
 123     volScalarField& cellDistance = cellDistancePtr_();
 124
 125     // Internal field
 126     {
 127         const pointField& cc = fvm.cellCentres();
 128         scalarField& fld = cellDistance.internalField();
 129
 130         forAll(cc, i)
 131         {
 132             // Signed distance
 133             fld[i] = (cc[i] - plane_.refPoint()) & plane_.normal();
 134         }
 135     }
 136
 137     // Patch fields
 138     {
 139         forAll(cellDistance.boundaryField(), patchI)
 140         {
 141             if
 142             (
 143                 isA<emptyFvPatchScalarField>
 144                 (
 145                     cellDistance.boundaryField()[patchI]
 146                 )
 147             )
 148             {
 149                 cellDistance.boundaryField().set
 150                 (
 151                     patchI,
 152                     new calculatedFvPatchScalarField
 153                     (
 154                         fvm.boundary()[patchI],
 155                         cellDistance
 156                     )
 157                 );
 158
 159                 const polyPatch& pp = fvm.boundary()[patchI].patch();
 160                 pointField::subField cc = pp.patchSlice(fvm.faceCentres());
 161
 162                 fvPatchScalarField& fld = cellDistance.boundaryField()[patchI];
 163                 fld.setSize(pp.size());
 164                 forAll(fld, i)
 165                 {
 166                     fld[i] = (cc[i] - plane_.refPoint()) & plane_.normal();
 167                 }
 168             }
 169             else
 170             {
 171                 const pointField& cc = fvm.C().boundaryField()[patchI];
 172                 fvPatchScalarField& fld = cellDistance.boundaryField()[patchI];
 173
 174                 forAll(fld, i)
 175                 {
 176                     fld[i] = (cc[i] - plane_.refPoint()) & plane_.normal();
 177                 }
 178             }
 179         }
 180     }
 181
 182
 183     // On processor patches the mesh.C() will already be the cell centre
 184     // on the opposite side so no need to swap cellDistance.
 185
 186
 187     // Distance to points
 188     pointDistance_.setSize(fvm.nPoints());
 189     {
 190         const pointField& pts = fvm.points();
 191
 192         forAll(pointDistance_, i)
 193         {
 194             pointDistance_[i] = (pts[i] - plane_.refPoint()) & plane_.normal();
 195         }
 196     }
 197
 198
 199     if (debug)
 200     {
 201         Pout<< "Writing cell distance:" << cellDistance.objectPath() << endl;
 202         cellDistance.write();
 203         pointScalarField pDist
 204         (
 205             IOobject
 206             (
 207                 "pointDistance",
 208                 fvm.time().timeName(),
 209                 fvm.time(),
 210                 IOobject::NO_READ,
 211                 IOobject::NO_WRITE,
 212                 false
 213             ),
 214             pointMesh::New(fvm),
 215             dimensionedScalar("zero", dimLength, 0)
 216         );
 217         pDist.internalField() = pointDistance_;
 218
 219         Pout<< "Writing point distance:" << pDist.objectPath() << endl;
 220         pDist.write();
 221     }
 222
 223
 224     //- Direct from cell field and point field.
 225     isoSurfPtr_.reset
 226     (
 227         new isoSurface
 228         (
 229             cellDistance,
 230             pointDistance_,
 231             0.0,
 232             regularise_,
 233             mergeTol_
 234         )
 235         //new isoSurfaceCell
 236         //(
 237         //    fvm,
 238         //    cellDistance,
 239         //    pointDistance_,
 240         //    0.0,
 241         //    regularise_,
 242         //    mergeTol_
 243         //)
 244     );
 245
 246     if (debug)
 247     {
 248         print(Pout);
 249         Pout<< endl;
 250     }
 251 }
 252
 253
 254 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 255
 256 Foam::sampledCuttingPlane::sampledCuttingPlane
 257 (
 258     const word& name,
 259     const polyMesh& mesh,
 260     const dictionary& dict
 261 )
 262 :
 263     sampledSurface(name, mesh, dict),
 264     plane_(dict),
 265     mergeTol_(dict.lookupOrDefault("mergeTol", 1E-6)),
 266     regularise_(dict.lookupOrDefault("regularise", true)),
 267     average_(dict.lookupOrDefault("average", false)),
 268     zoneID_(dict.lookupOrDefault("zone", word::null), mesh.cellZones()),
 269     exposedPatchName_(word::null),
 270     needsUpdate_(true),
 271     subMeshPtr_(NULL),
 272     cellDistancePtr_(NULL),
 273     isoSurfPtr_(NULL),
 274     facesPtr_(NULL)
 275 {
 276     if (zoneID_.index() != -1)
 277     {
 278         dict.lookup("exposedPatchName") >> exposedPatchName_;
 279
 280         if (mesh.boundaryMesh().findPatchID(exposedPatchName_) == -1)
 281         {
 282             FatalErrorIn
 283             (
 284                 "sampledCuttingPlane::sampledCuttingPlane"
 285                 "(const word&, const polyMesh&, const dictionary&)"
 286             )   << "Cannot find patch " << exposedPatchName_
 287                 << " in which to put exposed faces." << endl
 288                 << "Valid patches are " << mesh.boundaryMesh().names()
 289                 << exit(FatalError);
 290         }
 291
 292         if (debug && zoneID_.index() != -1)
 293         {
 294             Info<< "Restricting to cellZone " << zoneID_.name()
 295                 << " with exposed internal faces into patch "
 296                 << exposedPatchName_ << endl;
 297         }
 298     }
 299 }
 300
 301
 302 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 303
 304 Foam::sampledCuttingPlane::~sampledCuttingPlane()
 305 {}
 306
 307
 308 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 309
 310 bool Foam::sampledCuttingPlane::needsUpdate() const
 311 {
 312     return needsUpdate_;
 313 }
 314
 315
 316 bool Foam::sampledCuttingPlane::expire()
 317 {
 318     if (debug)
 319     {
 320         Pout<< "sampledCuttingPlane::expire :"
 321             << " have-facesPtr_:" << facesPtr_.valid()
 322             << " needsUpdate_:" << needsUpdate_ << endl;
 323     }
 324
 325     // Clear any stored topologies
 326     facesPtr_.clear();
 327
 328     // Clear derived data
 329     clearGeom();
 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::sampledCuttingPlane::update()
 343 {
 344     if (debug)
 345     {
 346         Pout<< "sampledCuttingPlane::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::sampledCuttingPlane::sample
 365 (
 366     const volScalarField& vField
 367 ) const
 368 {
 369     return sampleField(vField);
 370 }
 371
 372
 373 Foam::tmp<Foam::vectorField>
 374 Foam::sampledCuttingPlane::sample
 375 (
 376     const volVectorField& vField
 377 ) const
 378 {
 379     return sampleField(vField);
 380 }
 381
 382
 383 Foam::tmp<Foam::sphericalTensorField>
 384 Foam::sampledCuttingPlane::sample
 385 (
 386     const volSphericalTensorField& vField
 387 ) const
 388 {
 389     return sampleField(vField);
 390 }
 391
 392
 393 Foam::tmp<Foam::symmTensorField>
 394 Foam::sampledCuttingPlane::sample
 395 (
 396     const volSymmTensorField& vField
 397 ) const
 398 {
 399     return sampleField(vField);
 400 }
 401
 402
 403 Foam::tmp<Foam::tensorField>
 404 Foam::sampledCuttingPlane::sample
 405 (
 406     const volTensorField& vField
 407 ) const
 408 {
 409     return sampleField(vField);
 410 }
 411
 412
 413 Foam::tmp<Foam::scalarField>
 414 Foam::sampledCuttingPlane::interpolate
 415 (
 416     const interpolation<scalar>& interpolator
 417 ) const
 418 {
 419     return interpolateField(interpolator);
 420 }
 421
 422
 423 Foam::tmp<Foam::vectorField>
 424 Foam::sampledCuttingPlane::interpolate
 425 (
 426     const interpolation<vector>& interpolator
 427 ) const
 428 {
 429     return interpolateField(interpolator);
 430 }
 431
 432 Foam::tmp<Foam::sphericalTensorField>
 433 Foam::sampledCuttingPlane::interpolate
 434 (
 435     const interpolation<sphericalTensor>& interpolator
 436 ) const
 437 {
 438     return interpolateField(interpolator);
 439 }
 440
 441
 442 Foam::tmp<Foam::symmTensorField>
 443 Foam::sampledCuttingPlane::interpolate
 444 (
 445     const interpolation<symmTensor>& interpolator
 446 ) const
 447 {
 448     return interpolateField(interpolator);
 449 }
 450
 451
 452 Foam::tmp<Foam::tensorField>
 453 Foam::sampledCuttingPlane::interpolate
 454 (
 455     const interpolation<tensor>& interpolator
 456 ) const
 457 {
 458     return interpolateField(interpolator);
 459 }
 460
 461
 462 void Foam::sampledCuttingPlane::print(Ostream& os) const
 463 {
 464     os  << "sampledCuttingPlane: " << name() << " :"
 465         << "  plane:" << plane_
 466         << "  faces:" << faces().size()
 467         << "  points:" << points().size();
 468 }
 469
 470
 471 // ************************************************************************* //