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 \*---------------------------------------------------------------------------*/
  26
  27 #include "surfaceToCell.H"
  28 #include "polyMesh.H"
  29 #include "meshSearch.H"
  30 #include "triSurface.H"
  31 #include "triSurfaceSearch.H"
  32 #include "cellClassification.H"
  33
  34 #include "addToRunTimeSelectionTable.H"
  35
  36 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  37
  38 namespace Foam
  39 {
  40
  41 defineTypeNameAndDebug(surfaceToCell, 0);
  42
  43 addToRunTimeSelectionTable(topoSetSource, surfaceToCell, word);
  44
  45 addToRunTimeSelectionTable(topoSetSource, surfaceToCell, istream);
  46
  47 }
  48
  49
  50 Foam::topoSetSource::addToUsageTable Foam::surfaceToCell::usage_
  51 (
  52     surfaceToCell::typeName,
  53     "\n    Usage: surfaceToCell"
  54     "<surface> <outsidePoints> <cut> <inside> <outside> <near> <curvature>\n\n"
  55     "    <surface> name of triSurface\n"
  56     "    <outsidePoints> list of points that define outside\n"
  57     "    <cut> boolean whether to include cells cut by surface\n"
  58     "    <inside>   ,,                 ,,       inside surface\n"
  59     "    <outside>  ,,                 ,,       outside surface\n"
  60     "    <near> scalar; include cells with centre <= near to surface\n"
  61     "    <curvature> scalar; include cells close to strong curvature"
  62     " on surface\n"
  63     "    (curvature defined as difference in surface normal at nearest"
  64     " point on surface for each vertex of cell)\n\n"
  65 );
  66
  67
  68 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  69
  70 Foam::label Foam::surfaceToCell::getNearest
  71 (
  72     const triSurfaceSearch& querySurf,
  73     const label pointI,
  74     const point& pt,
  75     const vector& span,
  76     Map<label>& cache
  77 )
  78 {
  79     Map<label>::const_iterator iter = cache.find(pointI);
  80
  81     if (iter != cache.end())
  82     {
  83         // Found cached answer
  84         return iter();
  85     }
  86     else
  87     {
  88         pointIndexHit inter = querySurf.nearest(pt, span);
  89
  90         // Triangle label (can be -1)
  91         label triI = inter.index();
  92
  93         // Store triangle on point
  94         cache.insert(pointI, triI);
  95
  96         return triI;
  97     }
  98 }
  99
 100
 101 // Return true if nearest surface to points on cell makes largish angle
 102 // with nearest surface to cell centre. Returns false otherwise. Points visited
 103 // are cached in pointToNearest
 104 bool Foam::surfaceToCell::differingPointNormals
 105 (
 106     const triSurfaceSearch& querySurf,
 107
 108     const vector& span,         // current search span
 109     const label cellI,
 110     const label cellTriI,       // nearest (to cell centre) surface triangle
 111
 112     Map<label>& pointToNearest  // cache for nearest triangle to point
 113 ) const
 114 {
 115     const triSurface& surf = querySurf.surface();
 116     const vectorField& normals = surf.faceNormals();
 117
 118     const faceList& faces = mesh().faces();
 119     const pointField& points = mesh().points();
 120
 121     const labelList& cFaces = mesh().cells()[cellI];
 122
 123     forAll(cFaces, cFaceI)
 124     {
 125         const face& f = faces[cFaces[cFaceI]];
 126
 127         forAll(f, fp)
 128         {
 129             label pointI = f[fp];
 130
 131             label pointTriI =
 132                 getNearest
 133                 (
 134                     querySurf,
 135                     pointI,
 136                     points[pointI],
 137                     span,
 138                     pointToNearest
 139                 );
 140
 141             if (pointTriI != -1 && pointTriI != cellTriI)
 142             {
 143                 scalar cosAngle = normals[pointTriI] & normals[cellTriI];
 144
 145                 if (cosAngle < 0.9)
 146                 {
 147                     return true;
 148                 }
 149             }
 150         }
 151     }
 152     return false;
 153 }
 154
 155
 156 void Foam::surfaceToCell::combine(topoSet& set, const bool add) const
 157 {
 158     cpuTime timer;
 159
 160     if (includeCut_ || includeInside_ || includeOutside_)
 161     {
 162         //
 163         // Cut cells with surface and classify cells
 164         //
 165
 166
 167         // Construct search engine on mesh
 168
 169         meshSearch queryMesh(mesh_, true);
 170
 171
 172         // Check all 'outside' points
 173         forAll(outsidePoints_, outsideI)
 174         {
 175             const point& outsidePoint = outsidePoints_[outsideI];
 176
 177             // Find cell point is in. Linear search.
 178             if (queryMesh.findCell(outsidePoint, -1, false) == -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 // ************************************************************************* //