src/autoMesh/autoHexMesh/shellSurfaces/shellSurfaces.C

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