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