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) 2004-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
  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         )
 234         //new isoSurfaceCell
 235         //(
 236         //    fvm,
 237         //    cellDistance,
 238         //    pointDistance_,
 239         //    0.0,
 240         //    regularise_
 241         //)
 242     );
 243
 244     if (debug)
 245     {
 246         print(Pout);
 247         Pout<< endl;
 248     }
 249 }
 250
 251
 252 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 253
 254 Foam::sampledCuttingPlane::sampledCuttingPlane
 255 (
 256     const word& name,
 257     const polyMesh& mesh,
 258     const dictionary& dict
 259 )
 260 :
 261     sampledSurface(name, mesh, dict),
 262     plane_(dict),
 263     mergeTol_(dict.lookupOrDefault("mergeTol", 1E-6)),
 264     regularise_(dict.lookupOrDefault("regularise", true)),
 265     average_(dict.lookupOrDefault("average", false)),
 266     zoneID_(dict.lookupOrDefault("zone", word::null), mesh.cellZones()),
 267     exposedPatchName_(word::null),
 268     needsUpdate_(true),
 269     subMeshPtr_(NULL),
 270     cellDistancePtr_(NULL),
 271     isoSurfPtr_(NULL),
 272     facesPtr_(NULL)
 273 {
 274     if (zoneID_.index() != -1)
 275     {
 276         dict.lookup("exposedPatchName") >> exposedPatchName_;
 277
 278         if (mesh.boundaryMesh().findPatchID(exposedPatchName_) == -1)
 279         {
 280             FatalErrorIn
 281             (
 282                 "sampledCuttingPlane::sampledCuttingPlane"
 283                 "(const word&, const polyMesh&, const dictionary&)"
 284             )   << "Cannot find patch " << exposedPatchName_
 285                 << " in which to put exposed faces." << endl
 286                 << "Valid patches are " << mesh.boundaryMesh().names()
 287                 << exit(FatalError);
 288         }
 289
 290         if (debug && zoneID_.index() != -1)
 291         {
 292             Info<< "Restricting to cellZone " << zoneID_.name()
 293                 << " with exposed internal faces into patch "
 294                 << exposedPatchName_ << endl;
 295         }
 296     }
 297 }
 298
 299
 300 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 301
 302 Foam::sampledCuttingPlane::~sampledCuttingPlane()
 303 {}
 304
 305
 306 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 307
 308 bool Foam::sampledCuttingPlane::needsUpdate() const
 309 {
 310     return needsUpdate_;
 311 }
 312
 313
 314 bool Foam::sampledCuttingPlane::expire()
 315 {
 316     if (debug)
 317     {
 318         Pout<< "sampledCuttingPlane::expire :"
 319             << " have-facesPtr_:" << facesPtr_.valid()
 320             << " needsUpdate_:" << needsUpdate_ << endl;
 321     }
 322
 323     // Clear any stored topologies
 324     facesPtr_.clear();
 325
 326     // Clear derived data
 327     clearGeom();
 328
 329     // already marked as expired
 330     if (needsUpdate_)
 331     {
 332         return false;
 333     }
 334
 335     needsUpdate_ = true;
 336     return true;
 337 }
 338
 339
 340 bool Foam::sampledCuttingPlane::update()
 341 {
 342     if (debug)
 343     {
 344         Pout<< "sampledCuttingPlane::update :"
 345             << " have-facesPtr_:" << facesPtr_.valid()
 346             << " needsUpdate_:" << needsUpdate_ << endl;
 347     }
 348
 349     if (!needsUpdate_)
 350     {
 351         return false;
 352     }
 353
 354     createGeometry();
 355
 356     needsUpdate_ = false;
 357     return true;
 358 }
 359
 360
 361 Foam::tmp<Foam::scalarField>
 362 Foam::sampledCuttingPlane::sample
 363 (
 364     const volScalarField& vField
 365 ) const
 366 {
 367     return sampleField(vField);
 368 }
 369
 370
 371 Foam::tmp<Foam::vectorField>
 372 Foam::sampledCuttingPlane::sample
 373 (
 374     const volVectorField& vField
 375 ) const
 376 {
 377     return sampleField(vField);
 378 }
 379
 380
 381 Foam::tmp<Foam::sphericalTensorField>
 382 Foam::sampledCuttingPlane::sample
 383 (
 384     const volSphericalTensorField& vField
 385 ) const
 386 {
 387     return sampleField(vField);
 388 }
 389
 390
 391 Foam::tmp<Foam::symmTensorField>
 392 Foam::sampledCuttingPlane::sample
 393 (
 394     const volSymmTensorField& vField
 395 ) const
 396 {
 397     return sampleField(vField);
 398 }
 399
 400
 401 Foam::tmp<Foam::tensorField>
 402 Foam::sampledCuttingPlane::sample
 403 (
 404     const volTensorField& vField
 405 ) const
 406 {
 407     return sampleField(vField);
 408 }
 409
 410
 411 Foam::tmp<Foam::scalarField>
 412 Foam::sampledCuttingPlane::interpolate
 413 (
 414     const interpolation<scalar>& interpolator
 415 ) const
 416 {
 417     return interpolateField(interpolator);
 418 }
 419
 420
 421 Foam::tmp<Foam::vectorField>
 422 Foam::sampledCuttingPlane::interpolate
 423 (
 424     const interpolation<vector>& interpolator
 425 ) const
 426 {
 427     return interpolateField(interpolator);
 428 }
 429
 430 Foam::tmp<Foam::sphericalTensorField>
 431 Foam::sampledCuttingPlane::interpolate
 432 (
 433     const interpolation<sphericalTensor>& interpolator
 434 ) const
 435 {
 436     return interpolateField(interpolator);
 437 }
 438
 439
 440 Foam::tmp<Foam::symmTensorField>
 441 Foam::sampledCuttingPlane::interpolate
 442 (
 443     const interpolation<symmTensor>& interpolator
 444 ) const
 445 {
 446     return interpolateField(interpolator);
 447 }
 448
 449
 450 Foam::tmp<Foam::tensorField>
 451 Foam::sampledCuttingPlane::interpolate
 452 (
 453     const interpolation<tensor>& interpolator
 454 ) const
 455 {
 456     return interpolateField(interpolator);
 457 }
 458
 459
 460 void Foam::sampledCuttingPlane::print(Ostream& os) const
 461 {
 462     os  << "sampledCuttingPlane: " << name() << " :"
 463         << "  plane:" << plane_
 464         << "  faces:" << faces().size()
 465         << "  points:" << points().size();
 466 }
 467
 468
 469 // ************************************************************************* //