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

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  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             if (queryMesh.findCell(outsidePoint, -1, false) == -1)
 178             {
 179                 FatalErrorIn("surfaceToCell::combine(topoSet&, const bool)")
 180                     << "outsidePoint " << outsidePoint
 181                     << " is not inside any cell"
 182                     << exit(FatalError);
 183             }
 184         }
 185
 186         // Cut faces with surface and classify cells
 187
 188         cellClassification cellType
 189         (
 190             mesh_,
 191             queryMesh,
 192             querySurf(),
 193             outsidePoints_
 194         );
 195
 196
 197         Info<< "    Marked inside/outside in = "
 198             << timer.cpuTimeIncrement() << " s" << endl << endl;
 199
 200
 201         forAll(cellType, cellI)
 202         {
 203             label cType = cellType[cellI];
 204
 205             if
 206             (
 207                 (
 208                     includeCut_
 209                  && (cType == cellClassification::CUT)
 210                 )
 211              || (
 212                     includeInside_
 213                  && (cType == cellClassification::INSIDE)
 214                 )
 215              || (
 216                     includeOutside_
 217                  && (cType == cellClassification::OUTSIDE)
 218                 )
 219             )
 220             {
 221                 addOrDelete(set, cellI, add);
 222             }
 223         }
 224     }
 225
 226
 227     if (nearDist_ > 0)
 228     {
 229         //
 230         // Determine distance to surface
 231         //
 232
 233         const pointField& ctrs = mesh_.cellCentres();
 234
 235         // Box dimensions to search in octree.
 236         const vector span(nearDist_, nearDist_, nearDist_);
 237
 238
 239         if (curvature_ < -1)
 240         {
 241             Info<< "    Selecting cells with cellCentre closer than "
 242                 << nearDist_ << " to surface" << endl;
 243
 244             // No need to test curvature. Insert near cells into set.
 245
 246             forAll(ctrs, cellI)
 247             {
 248                 const point& c = ctrs[cellI];
 249
 250                 pointIndexHit inter = querySurf().nearest(c, span);
 251
 252                 if (inter.hit() && (mag(inter.hitPoint() - c) < nearDist_))
 253                 {
 254                     addOrDelete(set, cellI, add);
 255                 }
 256             }
 257
 258             Info<< "    Determined nearest surface point in = "
 259                 << timer.cpuTimeIncrement() << " s" << endl << endl;
 260
 261         }
 262         else
 263         {
 264             // Test near cells for curvature
 265
 266             Info<< "    Selecting cells with cellCentre closer than "
 267                 << nearDist_ << " to surface and curvature factor"
 268                 << " less than " << curvature_ << endl;
 269
 270             // Cache for nearest surface triangle for a point
 271             Map<label> pointToNearest(mesh_.nCells()/10);
 272
 273             forAll(ctrs, cellI)
 274             {
 275                 const point& c = ctrs[cellI];
 276
 277                 pointIndexHit inter = querySurf().nearest(c, span);
 278
 279                 if (inter.hit() && (mag(inter.hitPoint() - c) < nearDist_))
 280                 {
 281                     if
 282                     (
 283                         differingPointNormals
 284                         (
 285                             querySurf(),
 286                             span,
 287                             cellI,
 288                             inter.index(),      // nearest surface triangle
 289                             pointToNearest
 290                         )
 291                     )
 292                     {
 293                         addOrDelete(set, cellI, add);
 294                     }
 295                 }
 296             }
 297
 298             Info<< "    Determined nearest surface point in = "
 299                 << timer.cpuTimeIncrement() << " s" << endl << endl;
 300         }
 301     }
 302 }
 303
 304
 305 void Foam::surfaceToCell::checkSettings() const
 306 {
 307     if
 308     (
 309         (nearDist_ < 0)
 310      && (curvature_ < -1)
 311      && (
 312             (includeCut_ && includeInside_ && includeOutside_)
 313          || (!includeCut_ && !includeInside_ && !includeOutside_)
 314         )
 315     )
 316     {
 317         FatalErrorIn
 318         (
 319             "surfaceToCell:checkSettings()"
 320         )   << "Illegal include cell specification."
 321             << " Result would be either all or no cells." << endl
 322             << "Please set one of includeCut, includeInside, includeOutside"
 323             << " to true, set nearDistance to a value > 0"
 324             << " or set curvature to a value -1 .. 1."
 325             << exit(FatalError);
 326     }
 327 }
 328
 329
 330 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 331
 332 // Construct from components
 333 Foam::surfaceToCell::surfaceToCell
 334 (
 335     const polyMesh& mesh,
 336     const fileName& surfName,
 337     const pointField& outsidePoints,
 338     const bool includeCut,
 339     const bool includeInside,
 340     const bool includeOutside,
 341     const scalar nearDist,
 342     const scalar curvature
 343 )
 344 :
 345     topoSetSource(mesh),
 346     surfName_(surfName),
 347     outsidePoints_(outsidePoints),
 348     includeCut_(includeCut),
 349     includeInside_(includeInside),
 350     includeOutside_(includeOutside),
 351     nearDist_(nearDist),
 352     curvature_(curvature),
 353     surfPtr_(new triSurface(surfName_)),
 354     querySurfPtr_(new triSurfaceSearch(*surfPtr_)),
 355     IOwnPtrs_(true)
 356 {
 357     checkSettings();
 358 }
 359
 360
 361 // Construct from components. Externally supplied surface.
 362 Foam::surfaceToCell::surfaceToCell
 363 (
 364     const polyMesh& mesh,
 365     const fileName& surfName,
 366     const triSurface& surf,
 367     const triSurfaceSearch& querySurf,
 368     const pointField& outsidePoints,
 369     const bool includeCut,
 370     const bool includeInside,
 371     const bool includeOutside,
 372     const scalar nearDist,
 373     const scalar curvature
 374 )
 375 :
 376     topoSetSource(mesh),
 377     surfName_(surfName),
 378     outsidePoints_(outsidePoints),
 379     includeCut_(includeCut),
 380     includeInside_(includeInside),
 381     includeOutside_(includeOutside),
 382     nearDist_(nearDist),
 383     curvature_(curvature),
 384     surfPtr_(&surf),
 385     querySurfPtr_(&querySurf),
 386     IOwnPtrs_(false)
 387 {
 388     checkSettings();
 389 }
 390
 391
 392 // Construct from dictionary
 393 Foam::surfaceToCell::surfaceToCell
 394 (
 395     const polyMesh& mesh,
 396     const dictionary& dict
 397 )
 398 :
 399     topoSetSource(mesh),
 400     surfName_(dict.lookup("file")),
 401     outsidePoints_(dict.lookup("outsidePoints")),
 402     includeCut_(readBool(dict.lookup("includeCut"))),
 403     includeInside_(readBool(dict.lookup("includeInside"))),
 404     includeOutside_(readBool(dict.lookup("includeOutside"))),
 405     nearDist_(readScalar(dict.lookup("nearDistance"))),
 406     curvature_(readScalar(dict.lookup("curvature"))),
 407     surfPtr_(new triSurface(surfName_)),
 408     querySurfPtr_(new triSurfaceSearch(*surfPtr_)),
 409     IOwnPtrs_(true)
 410 {
 411     checkSettings();
 412 }
 413
 414
 415 // Construct from Istream
 416 Foam::surfaceToCell::surfaceToCell
 417 (
 418     const polyMesh& mesh,
 419     Istream& is
 420 )
 421 :
 422     topoSetSource(mesh),
 423     surfName_(checkIs(is)),
 424     outsidePoints_(checkIs(is)),
 425     includeCut_(readBool(checkIs(is))),
 426     includeInside_(readBool(checkIs(is))),
 427     includeOutside_(readBool(checkIs(is))),
 428     nearDist_(readScalar(checkIs(is))),
 429     curvature_(readScalar(checkIs(is))),
 430     surfPtr_(new triSurface(surfName_)),
 431     querySurfPtr_(new triSurfaceSearch(*surfPtr_)),
 432     IOwnPtrs_(true)
 433 {
 434     checkSettings();
 435 }
 436
 437
 438 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 439
 440 Foam::surfaceToCell::~surfaceToCell()
 441 {
 442     if (IOwnPtrs_)
 443     {
 444         if (surfPtr_)
 445         {
 446             delete surfPtr_;
 447         }
 448         if (querySurfPtr_)
 449         {
 450             delete querySurfPtr_;
 451         }
 452     }
 453 }
 454
 455
 456 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 457
 458 void Foam::surfaceToCell::applyToSet
 459 (
 460     const topoSetSource::setAction action,
 461     topoSet& set
 462 ) const
 463 {
 464     if ( (action == topoSetSource::NEW) || (action == topoSetSource::ADD))
 465     {
 466         Info<< "    Adding cells in relation to surface " << surfName_
 467             << " ..." << endl;
 468
 469         combine(set, true);
 470     }
 471     else if (action == topoSetSource::DELETE)
 472     {
 473         Info<< "    Removing cells in relation to surface " << surfName_
 474             << " ..." << endl;
 475
 476         combine(set, false);
 477     }
 478 }
 479
 480
 481 // ************************************************************************* //