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