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