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