src/mesh/autoMesh/autoHexMesh/shellSurfaces/shellSurfaces.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 "shellSurfaces.H"
  27 #include "searchableSurface.H"
  28 #include "boundBox.H"
  29 #include "triSurfaceMesh.H"
  30 #include "refinementSurfaces.H"
  31 #include "searchableSurfaces.H"
  32 #include "orientedSurface.H"
  33 #include "pointIndexHit.H"
  34
  35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  36
  37
  38 namespace Foam
  39 {
  40
  41 template<>
  42 const char*
  43 NamedEnum<shellSurfaces::refineMode, 3>::
  44 names[] =
  45 {
  46     "inside",
  47     "outside",
  48     "distance"
  49 };
  50
  51 const NamedEnum<shellSurfaces::refineMode, 3> shellSurfaces::refineModeNames_;
  52
  53 } // End namespace Foam
  54
  55
  56
  57 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  58
  59 void Foam::shellSurfaces::setAndCheckLevels
  60 (
  61     const label shellI,
  62     const List<Tuple2<scalar, label> >& distLevels
  63 )
  64 {
  65     if (modes_[shellI] != DISTANCE && distLevels.size() != 1)
  66     {
  67         FatalErrorIn
  68         (
  69             "shellSurfaces::shellSurfaces"
  70             "(const searchableSurfaces&, const dictionary&)"
  71         )   << "For refinement mode "
  72             << refineModeNames_[modes_[shellI]]
  73             << " specify only one distance+level."
  74             << " (its distance gets discarded)"
  75             << exit(FatalError);
  76     }
  77     // Extract information into separate distance and level
  78     distances_[shellI].setSize(distLevels.size());
  79     levels_[shellI].setSize(distLevels.size());
  80
  81     forAll(distLevels, j)
  82     {
  83         distances_[shellI][j] = distLevels[j].first();
  84         levels_[shellI][j] = distLevels[j].second();
  85
  86         // Check in incremental order
  87         if (j > 0)
  88         {
  89             if
  90             (
  91                 (distances_[shellI][j] <= distances_[shellI][j-1])
  92              || (levels_[shellI][j] > levels_[shellI][j-1])
  93             )
  94             {
  95                 FatalErrorIn
  96                 (
  97                     "shellSurfaces::shellSurfaces"
  98                     "(const searchableSurfaces&, const dictionary&)"
  99                 )   << "For refinement mode "
 100                     << refineModeNames_[modes_[shellI]]
 101                     << " : Refinement should be specified in order"
 102                     << " of increasing distance"
 103                     << " (and decreasing refinement level)." << endl
 104                     << "Distance:" << distances_[shellI][j]
 105                     << " refinementLevel:" << levels_[shellI][j]
 106                     << exit(FatalError);
 107             }
 108         }
 109     }
 110
 111     const searchableSurface& shell = allGeometry_[shells_[shellI]];
 112
 113     if (modes_[shellI] == DISTANCE)
 114     {
 115         Info<< "Refinement level according to distance to "
 116             << shell.name() << endl;
 117         forAll(levels_[shellI], j)
 118         {
 119             Info<< "    level " << levels_[shellI][j]
 120                 << " for all cells within " << distances_[shellI][j]
 121                 << " meter." << endl;
 122         }
 123     }
 124     else
 125     {
 126         if (!allGeometry_[shells_[shellI]].hasVolumeType())
 127         {
 128             FatalErrorIn
 129             (
 130                 "shellSurfaces::shellSurfaces"
 131                 "(const searchableSurfaces&"
 132                 ", const PtrList<dictionary>&)"
 133             )   << "Shell " << shell.name()
 134                 << " does not support testing for "
 135                 << refineModeNames_[modes_[shellI]] << endl
 136                 << "Probably it is not closed."
 137                 << exit(FatalError);
 138         }
 139
 140         if (modes_[shellI] == INSIDE)
 141         {
 142             Info<< "Refinement level " << levels_[shellI][0]
 143                 << " for all cells inside " << shell.name() << endl;
 144         }
 145         else
 146         {
 147             Info<< "Refinement level " << levels_[shellI][0]
 148                 << " for all cells outside " << shell.name() << endl;
 149         }
 150     }
 151 }
 152
 153
 154 // Specifically orient triSurfaces using a calculated point outside.
 155 // Done since quite often triSurfaces not of consistent orientation which
 156 // is (currently) necessary for sideness calculation
 157 void Foam::shellSurfaces::orient()
 158 {
 159     // Determine outside point.
 160     boundBox overallBb = boundBox::invertedBox;
 161
 162     bool hasSurface = false;
 163
 164     forAll(shells_, shellI)
 165     {
 166         const searchableSurface& s = allGeometry_[shells_[shellI]];
 167
 168         if (modes_[shellI] != DISTANCE && isA<triSurfaceMesh>(s))
 169         {
 170             const triSurfaceMesh& shell = refCast<const triSurfaceMesh>(s);
 171
 172             if (shell.triSurface::size())
 173             {
 174                 const pointField& points = shell.points();
 175
 176                 hasSurface = true;
 177
 178                 boundBox shellBb(points[0], points[0]);
 179                 // Assume surface is compact!
 180                 forAll(points, i)
 181                 {
 182                     const point& pt = points[i];
 183                     shellBb.min() = min(shellBb.min(), pt);
 184                     shellBb.max() = max(shellBb.max(), pt);
 185                 }
 186
 187                 overallBb.min() = min(overallBb.min(), shellBb.min());
 188                 overallBb.max() = max(overallBb.max(), shellBb.max());
 189             }
 190         }
 191     }
 192
 193     if (hasSurface)
 194     {
 195         const point outsidePt = overallBb.max() + overallBb.span();
 196
 197         //Info<< "Using point " << outsidePt << " to orient shells" << endl;
 198
 199         forAll(shells_, shellI)
 200         {
 201             const searchableSurface& s = allGeometry_[shells_[shellI]];
 202
 203             if (modes_[shellI] != DISTANCE && isA<triSurfaceMesh>(s))
 204             {
 205                 triSurfaceMesh& shell = const_cast<triSurfaceMesh&>
 206                 (
 207                     refCast<const triSurfaceMesh>(s)
 208                 );
 209
 210                 // Flip surface so outsidePt is outside.
 211                 bool anyFlipped = orientedSurface::orient
 212                 (
 213                     shell,
 214                     outsidePt,
 215                     true
 216                 );
 217
 218                 if (anyFlipped)
 219                 {
 220                     // orientedSurface will have done a clearOut of the surface.
 221                     // we could do a clearout of the triSurfaceMeshes::trees()
 222                     // but these aren't affected by orientation
 223                     // (except for cached
 224                     // sideness which should not be set at this point.
 225                     // !!Should check!)
 226
 227                     Info<< "shellSurfaces : Flipped orientation of surface "
 228                         << s.name()
 229                         << " so point " << outsidePt << " is outside." << endl;
 230                 }
 231             }
 232         }
 233     }
 234 }
 235
 236
 237 // Find maximum level of a shell.
 238 void Foam::shellSurfaces::findHigherLevel
 239 (
 240     const pointField& pt,
 241     const label shellI,
 242     labelList& maxLevel
 243 ) const
 244 {
 245     const labelList& levels = levels_[shellI];
 246
 247     if (modes_[shellI] == DISTANCE)
 248     {
 249         // Distance mode.
 250
 251         const scalarField& distances = distances_[shellI];
 252
 253         // Collect all those points that have a current maxLevel less than
 254         // (any of) the shell. Also collect the furthest distance allowable
 255         // to any shell with a higher level.
 256
 257         pointField candidates(pt.size());
 258         labelList candidateMap(pt.size());
 259         scalarField candidateDistSqr(pt.size());
 260         label candidateI = 0;
 261
 262         forAll(maxLevel, pointI)
 263         {
 264             forAllReverse(levels, levelI)
 265             {
 266                 if (levels[levelI] > maxLevel[pointI])
 267                 {
 268                     candidates[candidateI] = pt[pointI];
 269                     candidateMap[candidateI] = pointI;
 270                     candidateDistSqr[candidateI] = sqr(distances[levelI]);
 271                     candidateI++;
 272                     break;
 273                 }
 274             }
 275         }
 276         candidates.setSize(candidateI);
 277         candidateMap.setSize(candidateI);
 278         candidateDistSqr.setSize(candidateI);
 279
 280         // Do the expensive nearest test only for the candidate points.
 281         List<pointIndexHit> nearInfo;
 282         allGeometry_[shells_[shellI]].findNearest
 283         (
 284             candidates,
 285             candidateDistSqr,
 286             nearInfo
 287         );
 288
 289         // Update maxLevel
 290         forAll(nearInfo, candidateI)
 291         {
 292             if (nearInfo[candidateI].hit())
 293             {
 294                 // Check which level it actually is in.
 295                 label minDistI = findLower
 296                 (
 297                     distances,
 298                     mag(nearInfo[candidateI].hitPoint()-candidates[candidateI])
 299                 );
 300
 301                 label pointI = candidateMap[candidateI];
 302
 303                 // pt is inbetween shell[minDistI] and shell[minDistI+1]
 304                 maxLevel[pointI] = levels[minDistI+1];
 305             }
 306         }
 307     }
 308     else
 309     {
 310         // Inside/outside mode
 311
 312         // Collect all those points that have a current maxLevel less than the
 313         // shell.
 314
 315         pointField candidates(pt.size());
 316         labelList candidateMap(pt.size());
 317         label candidateI = 0;
 318
 319         forAll(maxLevel, pointI)
 320         {
 321             if (levels[0] > maxLevel[pointI])
 322             {
 323                 candidates[candidateI] = pt[pointI];
 324                 candidateMap[candidateI] = pointI;
 325                 candidateI++;
 326             }
 327         }
 328         candidates.setSize(candidateI);
 329         candidateMap.setSize(candidateI);
 330
 331         // Do the expensive nearest test only for the candidate points.
 332         List<searchableSurface::volumeType> volType;
 333         allGeometry_[shells_[shellI]].getVolumeType(candidates, volType);
 334
 335         forAll(volType, i)
 336         {
 337             label pointI = candidateMap[i];
 338
 339             if
 340             (
 341                 (
 342                     modes_[shellI] == INSIDE
 343                  && volType[i] == searchableSurface::INSIDE
 344                 )
 345              || (
 346                     modes_[shellI] == OUTSIDE
 347                  && volType[i] == searchableSurface::OUTSIDE
 348                 )
 349             )
 350             {
 351                 maxLevel[pointI] = levels[0];
 352             }
 353         }
 354     }
 355 }
 356
 357
 358 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 359
 360 Foam::shellSurfaces::shellSurfaces
 361 (
 362     const searchableSurfaces& allGeometry,
 363     const PtrList<dictionary>& shellDicts
 364 )
 365 :
 366     allGeometry_(allGeometry)
 367 {
 368     shells_.setSize(shellDicts.size());
 369     modes_.setSize(shellDicts.size());
 370     distances_.setSize(shellDicts.size());
 371     levels_.setSize(shellDicts.size());
 372
 373     forAll(shellDicts, shellI)
 374     {
 375         const dictionary& dict = shellDicts[shellI];
 376         const word name = dict.lookup("name");
 377         const word type = dict.lookup("type");
 378
 379         shells_[shellI] = allGeometry_.findSurfaceID(name);
 380
 381         if (shells_[shellI] == -1)
 382         {
 383             FatalErrorIn
 384             (
 385                 "shellSurfaces::shellSurfaces"
 386                 "(const searchableSurfaces&, const PtrList<dictionary>&)"
 387             )   << "No surface called " << name << endl
 388                 << "Valid surfaces are " << allGeometry_.names()
 389                 << exit(FatalError);
 390         }
 391
 392         modes_[shellI] = refineModeNames_.read(dict.lookup("refineMode"));
 393
 394         // Read pairs of distance+level
 395         setAndCheckLevels(shellI, dict.lookup("levels"));
 396     }
 397
 398     // Orient shell surfaces before any searching is done. Note that this
 399     // only needs to be done for inside or outside. Orienting surfaces
 400     // constructs lots of addressing which we want to avoid.
 401     orient();
 402 }
 403
 404
 405 Foam::shellSurfaces::shellSurfaces
 406 (
 407     const searchableSurfaces& allGeometry,
 408     const dictionary& shellsDict
 409 )
 410 :
 411     allGeometry_(allGeometry)
 412 {
 413     shells_.setSize(shellsDict.size());
 414     modes_.setSize(shellsDict.size());
 415     distances_.setSize(shellsDict.size());
 416     levels_.setSize(shellsDict.size());
 417
 418     label shellI = 0;
 419     forAllConstIter(dictionary, shellsDict, iter)
 420     {
 421         shells_[shellI] = allGeometry_.findSurfaceID(iter().keyword());
 422
 423         if (shells_[shellI] == -1)
 424         {
 425             FatalErrorIn
 426             (
 427                 "shellSurfaces::shellSurfaces"
 428                 "(const searchableSurfaces&, const dictionary>&"
 429             )   << "No surface called " << iter().keyword() << endl
 430                 << "Valid surfaces are " << allGeometry_.names()
 431                 << exit(FatalError);
 432         }
 433         const dictionary& dict = shellsDict.subDict(iter().keyword());
 434
 435         modes_[shellI] = refineModeNames_.read(dict.lookup("mode"));
 436
 437         // Read pairs of distance+level
 438         setAndCheckLevels(shellI, dict.lookup("levels"));
 439
 440         shellI++;
 441     }
 442
 443     // Orient shell surfaces before any searching is done. Note that this
 444     // only needs to be done for inside or outside. Orienting surfaces
 445     // constructs lots of addressing which we want to avoid.
 446     orient();
 447 }
 448
 449
 450 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 451
 452 // Highest shell level
 453 Foam::label Foam::shellSurfaces::maxLevel() const
 454 {
 455     label overallMax = 0;
 456     forAll(levels_, shellI)
 457     {
 458         overallMax = max(overallMax, max(levels_[shellI]));
 459     }
 460     return overallMax;
 461 }
 462
 463
 464 void Foam::shellSurfaces::findHigherLevel
 465 (
 466     const pointField& pt,
 467     const labelList& ptLevel,
 468     labelList& maxLevel
 469 ) const
 470 {
 471     // Maximum level of any shell. Start off with level of point.
 472     maxLevel = ptLevel;
 473
 474     forAll(shells_, shellI)
 475     {
 476         findHigherLevel(pt, shellI, maxLevel);
 477     }
 478 }
 479
 480
 481 // ************************************************************************* //