src/meshTools/sets/cellSources/surfaceToCell/surfaceToCell.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2011 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 "surfaceToCell.H"
  27 #include "polyMesh.H"
  28 #include "meshSearch.H"
  29 #include "triSurface.H"
  30 #include "triSurfaceSearch.H"
  31 #include "cellClassification.H"
  32
  33 #include "addToRunTimeSelectionTable.H"
  34
  35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  36
  37 namespace Foam
  38 {
  39
  40 defineTypeNameAndDebug(surfaceToCell, 0);
  41
  42 addToRunTimeSelectionTable(topoSetSource, surfaceToCell, word);
  43
  44 addToRunTimeSelectionTable(topoSetSource, surfaceToCell, istream);
  45
  46 }
  47
  48
  49 Foam::topoSetSource::addToUsageTable Foam::surfaceToCell::usage_
  50 (
  51     surfaceToCell::typeName,
  52     "\n    Usage: surfaceToCell"
  53     "<surface> <outsidePoints> <cut> <inside> <outside> <near> <curvature>\n\n"
  54     "    <surface> name of triSurface\n"
  55     "    <outsidePoints> list of points that define outside\n"
  56     "    <cut> boolean whether to include cells cut by surface\n"
  57     "    <inside>   ,,                 ,,       inside surface\n"
  58     "    <outside>  ,,                 ,,       outside surface\n"
  59     "    <near> scalar; include cells with centre <= near to surface\n"
  60     "    <curvature> scalar; include cells close to strong curvature"
  61     " on surface\n"
  62     "    (curvature defined as difference in surface normal at nearest"
  63     " point on surface for each vertex of cell)\n\n"
  64 );
  65
  66
  67 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  68
  69 Foam::label Foam::surfaceToCell::getNearest
  70 (
  71     const triSurfaceSearch& querySurf,
  72     const label pointI,
  73     const point& pt,
  74     const vector& span,
  75     Map<label>& cache
  76 )
  77 {
  78     Map<label>::const_iterator iter = cache.find(pointI);
  79
  80     if (iter != cache.end())
  81     {
  82         // Found cached answer
  83         return iter();
  84     }
  85     else
  86     {
  87         pointIndexHit inter = querySurf.nearest(pt, span);
  88
  89         // Triangle label (can be -1)
  90         label triI = inter.index();
  91
  92         // Store triangle on point
  93         cache.insert(pointI, triI);
  94
  95         return triI;
  96     }
  97 }
  98
  99
 100 // Return true if nearest surface to points on cell makes largish angle
 101 // with nearest surface to cell centre. Returns false otherwise. Points visited
 102 // are cached in pointToNearest
 103 bool Foam::surfaceToCell::differingPointNormals
 104 (
 105     const triSurfaceSearch& querySurf,
 106
 107     const vector& span,         // current search span
 108     const label cellI,
 109     const label cellTriI,       // nearest (to cell centre) surface triangle
 110
 111     Map<label>& pointToNearest  // cache for nearest triangle to point
 112 ) const
 113 {
 114     const triSurface& surf = querySurf.surface();
 115     const vectorField& normals = surf.faceNormals();
 116
 117     const faceList& faces = mesh().faces();
 118     const pointField& points = mesh().points();
 119
 120     const labelList& cFaces = mesh().cells()[cellI];
 121
 122     forAll(cFaces, cFaceI)
 123     {
 124         const face& f = faces[cFaces[cFaceI]];
 125
 126         forAll(f, fp)
 127         {
 128             label pointI = f[fp];
 129
 130             label pointTriI =
 131                 getNearest
 132                 (
 133                     querySurf,
 134                     pointI,
 135                     points[pointI],
 136                     span,
 137                     pointToNearest
 138                 );
 139
 140             if (pointTriI != -1 && pointTriI != cellTriI)
 141             {
 142                 scalar cosAngle = normals[pointTriI] & normals[cellTriI];
 143
 144                 if (cosAngle < 0.9)
 145                 {
 146                     return true;
 147                 }
 148             }
 149         }
 150     }
 151     return false;
 152 }
 153
 154
 155 void Foam::surfaceToCell::combine(topoSet& set, const bool add) const
 156 {
 157     cpuTime timer;
 158
 159     if (includeCut_ || includeInside_ || includeOutside_)
 160     {
 161         //
 162         // Cut cells with surface and classify cells
 163         //
 164
 165
 166         // Construct search engine on mesh
 167
 168         meshSearch queryMesh(mesh_, true);
 169
 170
 171         // Check all 'outside' points
 172         forAll(outsidePoints_, outsideI)
 173         {
 174             const point& outsidePoint = outsidePoints_[outsideI];
 175
 176             // Find cell point is in. Linear search.
 177             label cellI = queryMesh.findCell(outsidePoint, -1, false);
 178             if (returnReduce(cellI, maxOp<label>()) == -1)
 179             {
 180                 FatalErrorIn("surfaceToCell::combine(topoSet&, const bool)")
 181                     << "outsidePoint " << outsidePoint
 182                     << " is not inside any cell"
 183                     << exit(FatalError);
 184             }
 185         }
 186
 187         // Cut faces with surface and classify cells
 188
 189         cellClassification cellType
 190         (
 191             mesh_,
 192             queryMesh,
 193             querySurf(),
 194             outsidePoints_
 195         );
 196
 197
 198         Info<< "    Marked inside/outside in = "
 199             << timer.cpuTimeIncrement() << " s" << endl << endl;
 200
 201
 202         forAll(cellType, cellI)
 203         {
 204             label cType = cellType[cellI];
 205
 206             if
 207             (
 208                 (
 209                     includeCut_
 210                  && (cType == cellClassification::CUT)
 211                 )
 212              || (
 213                     includeInside_
 214                  && (cType == cellClassification::INSIDE)
 215                 )
 216              || (
 217                     includeOutside_
 218                  && (cType == cellClassification::OUTSIDE)
 219                 )
 220             )
 221             {
 222                 addOrDelete(set, cellI, add);
 223             }
 224         }
 225     }
 226
 227
 228     if (nearDist_ > 0)
 229     {
 230         //
 231         // Determine distance to surface
 232         //
 233
 234         const pointField& ctrs = mesh_.cellCentres();
 235
 236         // Box dimensions to search in octree.
 237         const vector span(nearDist_, nearDist_, nearDist_);
 238
 239
 240         if (curvature_ < -1)
 241         {
 242             Info<< "    Selecting cells with cellCentre closer than "
 243                 << nearDist_ << " to surface" << endl;
 244
 245             // No need to test curvature. Insert near cells into set.
 246
 247             forAll(ctrs, cellI)
 248             {
 249                 const point& c = ctrs[cellI];
 250
 251                 pointIndexHit inter = querySurf().nearest(c, span);
 252
 253                 if (inter.hit() && (mag(inter.hitPoint() - c) < nearDist_))
 254                 {
 255                     addOrDelete(set, cellI, add);
 256                 }
 257             }
 258
 259             Info<< "    Determined nearest surface point in = "
 260                 << timer.cpuTimeIncrement() << " s" << endl << endl;
 261
 262         }
 263         else
 264         {
 265             // Test near cells for curvature
 266
 267             Info<< "    Selecting cells with cellCentre closer than "
 268                 << nearDist_ << " to surface and curvature factor"
 269                 << " less than " << curvature_ << endl;
 270
 271             // Cache for nearest surface triangle for a point
 272             Map<label> pointToNearest(mesh_.nCells()/10);
 273
 274             forAll(ctrs, cellI)
 275             {
 276                 const point& c = ctrs[cellI];
 277
 278                 pointIndexHit inter = querySurf().nearest(c, span);
 279
 280                 if (inter.hit() && (mag(inter.hitPoint() - c) < nearDist_))
 281                 {
 282                     if
 283                     (
 284                         differingPointNormals
 285                         (
 286                             querySurf(),
 287                             span,
 288                             cellI,
 289                             inter.index(),      // nearest surface triangle
 290                             pointToNearest
 291                         )
 292                     )
 293                     {
 294                         addOrDelete(set, cellI, add);
 295                     }
 296                 }
 297             }
 298
 299             Info<< "    Determined nearest surface point in = "
 300                 << timer.cpuTimeIncrement() << " s" << endl << endl;
 301         }
 302     }
 303 }
 304
 305
 306 void Foam::surfaceToCell::checkSettings() const
 307 {
 308     if
 309     (
 310         (nearDist_ < 0)
 311      && (curvature_ < -1)
 312      && (
 313             (includeCut_ && includeInside_ && includeOutside_)
 314          || (!includeCut_ && !includeInside_ && !includeOutside_)
 315         )
 316     )
 317     {
 318         FatalErrorIn
 319         (
 320             "surfaceToCell:checkSettings()"
 321         )   << "Illegal include cell specification."
 322             << " Result would be either all or no cells." << endl
 323             << "Please set one of includeCut, includeInside, includeOutside"
 324             << " to true, set nearDistance to a value > 0"
 325             << " or set curvature to a value -1 .. 1."
 326             << exit(FatalError);
 327     }
 328 }
 329
 330
 331 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 332
 333 // Construct from components
 334 Foam::surfaceToCell::surfaceToCell
 335 (
 336     const polyMesh& mesh,
 337     const fileName& surfName,
 338     const pointField& outsidePoints,
 339     const bool includeCut,
 340     const bool includeInside,
 341     const bool includeOutside,
 342     const scalar nearDist,
 343     const scalar curvature
 344 )
 345 :
 346     topoSetSource(mesh),
 347     surfName_(surfName),
 348     outsidePoints_(outsidePoints),
 349     includeCut_(includeCut),
 350     includeInside_(includeInside),
 351     includeOutside_(includeOutside),
 352     nearDist_(nearDist),
 353     curvature_(curvature),
 354     surfPtr_(new triSurface(surfName_)),
 355     querySurfPtr_(new triSurfaceSearch(*surfPtr_)),
 356     IOwnPtrs_(true)
 357 {
 358     checkSettings();
 359 }
 360
 361
 362 // Construct from components. Externally supplied surface.
 363 Foam::surfaceToCell::surfaceToCell
 364 (
 365     const polyMesh& mesh,
 366     const fileName& surfName,
 367     const triSurface& surf,
 368     const triSurfaceSearch& querySurf,
 369     const pointField& outsidePoints,
 370     const bool includeCut,
 371     const bool includeInside,
 372     const bool includeOutside,
 373     const scalar nearDist,
 374     const scalar curvature
 375 )
 376 :
 377     topoSetSource(mesh),
 378     surfName_(surfName),
 379     outsidePoints_(outsidePoints),
 380     includeCut_(includeCut),
 381     includeInside_(includeInside),
 382     includeOutside_(includeOutside),
 383     nearDist_(nearDist),
 384     curvature_(curvature),
 385     surfPtr_(&surf),
 386     querySurfPtr_(&querySurf),
 387     IOwnPtrs_(false)
 388 {
 389     checkSettings();
 390 }
 391
 392
 393 // Construct from dictionary
 394 Foam::surfaceToCell::surfaceToCell
 395 (
 396     const polyMesh& mesh,
 397     const dictionary& dict
 398 )
 399 :
 400     topoSetSource(mesh),
 401     surfName_(dict.lookup("file")),
 402     outsidePoints_(dict.lookup("outsidePoints")),
 403     includeCut_(readBool(dict.lookup("includeCut"))),
 404     includeInside_(readBool(dict.lookup("includeInside"))),
 405     includeOutside_(readBool(dict.lookup("includeOutside"))),
 406     nearDist_(readScalar(dict.lookup("nearDistance"))),
 407     curvature_(readScalar(dict.lookup("curvature"))),
 408     surfPtr_(new triSurface(surfName_)),
 409     querySurfPtr_(new triSurfaceSearch(*surfPtr_)),
 410     IOwnPtrs_(true)
 411 {
 412     checkSettings();
 413 }
 414
 415
 416 // Construct from Istream
 417 Foam::surfaceToCell::surfaceToCell
 418 (
 419     const polyMesh& mesh,
 420     Istream& is
 421 )
 422 :
 423     topoSetSource(mesh),
 424     surfName_(checkIs(is)),
 425     outsidePoints_(checkIs(is)),
 426     includeCut_(readBool(checkIs(is))),
 427     includeInside_(readBool(checkIs(is))),
 428     includeOutside_(readBool(checkIs(is))),
 429     nearDist_(readScalar(checkIs(is))),
 430     curvature_(readScalar(checkIs(is))),
 431     surfPtr_(new triSurface(surfName_)),
 432     querySurfPtr_(new triSurfaceSearch(*surfPtr_)),
 433     IOwnPtrs_(true)
 434 {
 435     checkSettings();
 436 }
 437
 438
 439 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 440
 441 Foam::surfaceToCell::~surfaceToCell()
 442 {
 443     if (IOwnPtrs_)
 444     {
 445         if (surfPtr_)
 446         {
 447             delete surfPtr_;
 448         }
 449         if (querySurfPtr_)
 450         {
 451             delete querySurfPtr_;
 452         }
 453     }
 454 }
 455
 456
 457 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 458
 459 void Foam::surfaceToCell::applyToSet
 460 (
 461     const topoSetSource::setAction action,
 462     topoSet& set
 463 ) const
 464 {
 465     if ( (action == topoSetSource::NEW) || (action == topoSetSource::ADD))
 466     {
 467         Info<< "    Adding cells in relation to surface " << surfName_
 468             << " ..." << endl;
 469
 470         combine(set, true);
 471     }
 472     else if (action == topoSetSource::DELETE)
 473     {
 474         Info<< "    Removing cells in relation to surface " << surfName_
 475             << " ..." << endl;
 476
 477         combine(set, false);
 478     }
 479 }
 480
 481
 482 // ************************************************************************* //