src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2008 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 "refinementSurfaces.H"
  28 #include "orientedSurface.H"
  29 #include "Time.H"
  30 #include "searchableSurfaces.H"
  31 #include "shellSurfaces.H"
  32 #include "triSurfaceMesh.H"
  33 #include "labelPair.H"
  34 #include "searchableSurfacesQueries.H"
  35
  36
  37 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  38
  39
  40 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  41
  42
  43 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  44
  45 // Construct from components
  46 Foam::refinementSurfaces::refinementSurfaces
  47 (
  48     const searchableSurfaces& allGeometry,
  49     const PtrList<dictionary>& surfaceDicts
  50 )
  51 :
  52     allGeometry_(allGeometry),
  53     surfaces_(surfaceDicts.size()),
  54     names_(surfaceDicts.size()),
  55     faceZoneNames_(surfaceDicts.size()),
  56     cellZoneNames_(surfaceDicts.size()),
  57     zoneInside_(surfaceDicts.size()),
  58     regionOffset_(surfaceDicts.size())
  59 {
  60     labelList globalMinLevel(surfaceDicts.size(), 0);
  61     labelList globalMaxLevel(surfaceDicts.size(), 0);
  62     scalarField globalAngle(surfaceDicts.size(), -GREAT);
  63     List<Map<label> > regionMinLevel(surfaceDicts.size());
  64     List<Map<label> > regionMaxLevel(surfaceDicts.size());
  65     List<Map<scalar> > regionAngle(surfaceDicts.size());
  66
  67     //wordList globalPatchType(surfaceDicts.size());
  68     //List<HashTable<word> > regionPatchType(surfaceDicts.size());
  69     //List<HashTable<word> > regionPatchName(surfaceDicts.size());
  70
  71     forAll(surfaceDicts, surfI)
  72     {
  73         const dictionary& dict = surfaceDicts[surfI];
  74
  75         dict.lookup("name") >> names_[surfI];
  76
  77         surfaces_[surfI] = allGeometry_.findSurfaceID(names_[surfI]);
  78
  79         // Global refinement level
  80         globalMinLevel[surfI] = readLabel(dict.lookup("minRefinementLevel"));
  81         globalMaxLevel[surfI] = readLabel(dict.lookup("maxRefinementLevel"));
  82
  83         // Global zone names per surface
  84         if (dict.found("faceZone"))
  85         {
  86             dict.lookup("faceZone") >> faceZoneNames_[surfI];
  87             dict.lookup("cellZone") >> cellZoneNames_[surfI];
  88             dict.lookup("zoneInside") >> zoneInside_[surfI];
  89         }
  90
  91         // Global perpendicular angle
  92         if (dict.found("perpendicularAngle"))
  93         {
  94             globalAngle[surfI] = readScalar(dict.lookup("perpendicularAngle"));
  95         }
  96
  97         //// Global patch name per surface
  98         //if (dict.found("patchType"))
  99         //{
 100         //    dict.lookup("patchType") >> globalPatchType[surfI];
 101         //}
 102
 103
 104         if (dict.found("regions"))
 105         {
 106             PtrList<dictionary> regionDicts(dict.lookup("regions"));
 107
 108             const wordList& regionNames =
 109                 allGeometry_[surfaces_[surfI]].regions();
 110
 111             forAll(regionDicts, dictI)
 112             {
 113                 const dictionary& regionDict = regionDicts[dictI];
 114
 115                 const word regionName(regionDict.lookup("name"));
 116
 117                 label regionI = findIndex(regionNames, regionName);
 118
 119                 if (regionI == -1)
 120                 {
 121                     FatalErrorIn
 122                     (
 123                         "refinementSurfaces::refinementSurfaces"
 124                         "(const IOobject&, const PtrList<dictionary>&)"
 125                     )   << "No region called " << regionName << " on surface "
 126                         << allGeometry_[surfaces_[surfI]].name() << endl
 127                         << "Valid regions are " << regionNames
 128                         << exit(FatalError);
 129                 }
 130
 131
 132                 label min = readLabel(regionDict.lookup("minRefinementLevel"));
 133                 label max = readLabel(regionDict.lookup("maxRefinementLevel"));
 134
 135                 bool hasInserted = regionMinLevel[surfI].insert(regionI, min);
 136                 if (!hasInserted)
 137                 {
 138                     FatalErrorIn
 139                     (
 140                         "refinementSurfaces::refinementSurfaces"
 141                         "(const IOobject&, const PtrList<dictionary>&)"
 142                     )   << "Duplicate region name " << regionName
 143                         << " on surface " << names_[surfI]
 144                         << exit(FatalError);
 145                 }
 146                 regionMaxLevel[surfI].insert(regionI, max);
 147
 148                 if (regionDict.found("perpendicularAngle"))
 149                 {
 150                     regionAngle[surfI].insert
 151                     (
 152                         regionI,
 153                         readScalar(regionDict.lookup("perpendicularAngle"))
 154                     );
 155                 }
 156             }
 157         }
 158     }
 159
 160
 161     // Check for duplicate surface names
 162     {
 163         HashTable<label> surfaceNames(names_.size());
 164
 165         forAll(names_, surfI)
 166         {
 167             if (!surfaceNames.insert(names_[surfI], surfI))
 168             {
 169                 FatalErrorIn
 170                 (
 171                     "refinementSurfaces::refinementSurfaces"
 172                     "(const IOobject&, const PtrList<dictionary>&)"
 173                 )   << "Duplicate surface name " << names_[surfI] << endl
 174                     << "Previous occurrence of name at surface "
 175                     << surfaceNames[names_[surfI]]
 176                     << exit(FatalError);
 177             }
 178         }
 179     }
 180
 181     // Calculate local to global region offset
 182     label nRegions = 0;
 183
 184     forAll(surfaceDicts, surfI)
 185     {
 186         regionOffset_[surfI] = nRegions;
 187
 188         nRegions += allGeometry_[surfaces_[surfI]].regions().size();
 189     }
 190
 191     // Rework surface specific information into information per global region
 192     minLevel_.setSize(nRegions);
 193     minLevel_ = 0;
 194     maxLevel_.setSize(nRegions);
 195     maxLevel_ = 0;
 196     perpendicularAngle_.setSize(nRegions);
 197     perpendicularAngle_ = -GREAT;
 198     //patchName_.setSize(nRegions);
 199     //patchType_.setSize(nRegions);
 200
 201     forAll(surfaceDicts, surfI)
 202     {
 203         label nRegions = allGeometry_[surfaces_[surfI]].regions().size();
 204
 205         // Initialise to global (i.e. per surface)
 206         for (label i = 0; i < nRegions; i++)
 207         {
 208             minLevel_[regionOffset_[surfI] + i] = globalMinLevel[surfI];
 209             maxLevel_[regionOffset_[surfI] + i] = globalMaxLevel[surfI];
 210             perpendicularAngle_[regionOffset_[surfI] + i] = globalAngle[surfI];
 211         }
 212
 213         // Overwrite with region specific information
 214         forAllConstIter(Map<label>, regionMinLevel[surfI], iter)
 215         {
 216             label globalRegionI = regionOffset_[surfI] + iter.key();
 217
 218             minLevel_[globalRegionI] = iter();
 219             maxLevel_[globalRegionI] = regionMaxLevel[surfI][iter.key()];
 220
 221             // Check validity
 222             if
 223             (
 224                 minLevel_[globalRegionI] < 0
 225              || maxLevel_[globalRegionI] < minLevel_[globalRegionI]
 226             )
 227             {
 228                 FatalErrorIn
 229                 (
 230                     "refinementSurfaces::refinementSurfaces"
 231                     "(const IOobject&, const PtrList<dictionary>&)"
 232                 )   << "Illegal level or layer specification for surface "
 233                     << names_[surfI]
 234                     << " : minLevel:" << minLevel_[globalRegionI]
 235                     << " maxLevel:" << maxLevel_[globalRegionI]
 236                     << exit(FatalError);
 237             }
 238         }
 239         forAllConstIter(Map<scalar>, regionAngle[surfI], iter)
 240         {
 241             label globalRegionI = regionOffset_[surfI] + iter.key();
 242
 243             perpendicularAngle_[globalRegionI] = regionAngle[surfI][iter.key()];
 244         }
 245
 246         //// Optional patch names and patch types
 247         //forAllConstIter(HashTable<word>, regionPatchName[surfI], iter)
 248         //{
 249         //    label regionI = findIndex(regionNames, iter.key());
 250         //    label globalRegionI = regionOffset_[surfI] + regionI;
 251         //
 252         //    patchName_[globalRegionI] = iter();
 253         //    patchType_[globalRegionI] = regionPatchType[surfI][iter.key()];
 254         //}
 255     }
 256 }
 257
 258
 259 Foam::refinementSurfaces::refinementSurfaces
 260 (
 261     const searchableSurfaces& allGeometry,
 262     const dictionary& surfacesDict
 263 )
 264 :
 265     allGeometry_(allGeometry),
 266     surfaces_(surfacesDict.size()),
 267     names_(surfacesDict.size()),
 268     faceZoneNames_(surfacesDict.size()),
 269     cellZoneNames_(surfacesDict.size()),
 270     zoneInside_(surfacesDict.size()),
 271     regionOffset_(surfacesDict.size())
 272 {
 273     labelList globalMinLevel(surfacesDict.size(), 0);
 274     labelList globalMaxLevel(surfacesDict.size(), 0);
 275     scalarField globalAngle(surfacesDict.size(), -GREAT);
 276     List<Map<label> > regionMinLevel(surfacesDict.size());
 277     List<Map<label> > regionMaxLevel(surfacesDict.size());
 278     List<Map<scalar> > regionAngle(surfacesDict.size());
 279
 280     label surfI = 0;
 281     forAllConstIter(dictionary, surfacesDict, iter)
 282     {
 283         names_[surfI] = iter().keyword();
 284
 285         surfaces_[surfI] = allGeometry_.findSurfaceID(names_[surfI]);
 286
 287         if (surfaces_[surfI] == -1)
 288         {
 289             FatalErrorIn
 290             (
 291                 "refinementSurfaces::refinementSurfaces"
 292                 "(const searchableSurfaces&, const dictionary>&"
 293             )   << "No surface called " << iter().keyword() << endl
 294                 << "Valid surfaces are " << allGeometry_.names()
 295                 << exit(FatalError);
 296         }
 297         const dictionary& dict = surfacesDict.subDict(iter().keyword());
 298
 299         const labelPair refLevel(dict.lookup("level"));
 300         globalMinLevel[surfI] = refLevel[0];
 301         globalMaxLevel[surfI] = refLevel[1];
 302
 303         // Global zone names per surface
 304         if (dict.found("faceZone"))
 305         {
 306             dict.lookup("faceZone") >> faceZoneNames_[surfI];
 307             dict.lookup("cellZone") >> cellZoneNames_[surfI];
 308             dict.lookup("zoneInside") >> zoneInside_[surfI];
 309         }
 310
 311         // Global perpendicular angle
 312         if (dict.found("perpendicularAngle"))
 313         {
 314             globalAngle[surfI] = readScalar(dict.lookup("perpendicularAngle"));
 315         }
 316
 317         if (dict.found("regions"))
 318         {
 319             const dictionary& regionsDict = dict.subDict("regions");
 320             const wordList& regionNames =
 321                 allGeometry_[surfaces_[surfI]].regions();
 322
 323             forAllConstIter(dictionary, regionsDict, iter)
 324             {
 325                 const word& key = iter().keyword();
 326
 327                 if (regionsDict.isDict(key))
 328                 {
 329                     // Get the dictionary for region iter.keyword()
 330                     const dictionary& regionDict = regionsDict.subDict(key);
 331
 332                     label regionI = findIndex(regionNames, key);
 333                     if (regionI == -1)
 334                     {
 335                         FatalErrorIn
 336                         (
 337                             "refinementSurfaces::refinementSurfaces"
 338                             "(const searchableSurfaces&, const dictionary>&"
 339                         )   << "No region called " << key << " on surface "
 340                             << allGeometry_[surfaces_[surfI]].name() << endl
 341                             << "Valid regions are " << regionNames
 342                             << exit(FatalError);
 343                     }
 344
 345                     const labelPair refLevel(regionDict.lookup("level"));
 346
 347                     regionMinLevel[surfI].insert(regionI, refLevel[0]);
 348                     regionMaxLevel[surfI].insert(regionI, refLevel[1]);
 349
 350                     if (regionDict.found("perpendicularAngle"))
 351                     {
 352                         regionAngle[surfI].insert
 353                         (
 354                             regionI,
 355                             readScalar(regionDict.lookup("perpendicularAngle"))
 356                         );
 357                     }
 358                 }
 359             }
 360         }
 361         surfI++;
 362     }
 363
 364     // Calculate local to global region offset
 365     label nRegions = 0;
 366
 367     forAll(surfacesDict, surfI)
 368     {
 369         regionOffset_[surfI] = nRegions;
 370         nRegions += allGeometry_[surfaces_[surfI]].regions().size();
 371     }
 372
 373     // Rework surface specific information into information per global region
 374     minLevel_.setSize(nRegions);
 375     minLevel_ = 0;
 376     maxLevel_.setSize(nRegions);
 377     maxLevel_ = 0;
 378     perpendicularAngle_.setSize(nRegions);
 379     perpendicularAngle_ = -GREAT;
 380
 381
 382     forAll(globalMinLevel, surfI)
 383     {
 384         label nRegions = allGeometry_[surfaces_[surfI]].regions().size();
 385
 386         // Initialise to global (i.e. per surface)
 387         for (label i = 0; i < nRegions; i++)
 388         {
 389             minLevel_[regionOffset_[surfI] + i] = globalMinLevel[surfI];
 390             maxLevel_[regionOffset_[surfI] + i] = globalMaxLevel[surfI];
 391             perpendicularAngle_[regionOffset_[surfI] + i] = globalAngle[surfI];
 392         }
 393
 394         // Overwrite with region specific information
 395         forAllConstIter(Map<label>, regionMinLevel[surfI], iter)
 396         {
 397             label globalRegionI = regionOffset_[surfI] + iter.key();
 398
 399             minLevel_[globalRegionI] = iter();
 400             maxLevel_[globalRegionI] = regionMaxLevel[surfI][iter.key()];
 401
 402             // Check validity
 403             if
 404             (
 405                 minLevel_[globalRegionI] < 0
 406              || maxLevel_[globalRegionI] < minLevel_[globalRegionI]
 407             )
 408             {
 409                 FatalErrorIn
 410                 (
 411                     "refinementSurfaces::refinementSurfaces"
 412                     "(const searchableSurfaces&, const dictionary>&"
 413                 )   << "Illegal level or layer specification for surface "
 414                     << names_[surfI]
 415                     << " : minLevel:" << minLevel_[globalRegionI]
 416                     << " maxLevel:" << maxLevel_[globalRegionI]
 417                     << exit(FatalError);
 418             }
 419         }
 420
 421         forAllConstIter(Map<scalar>, regionAngle[surfI], iter)
 422         {
 423             label globalRegionI = regionOffset_[surfI] + iter.key();
 424
 425             perpendicularAngle_[globalRegionI] = regionAngle[surfI][iter.key()];
 426         }
 427     }
 428 }
 429
 430
 431 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 432
 433 // Get indices of named surfaces (surfaces with cellZoneName)
 434 Foam::labelList Foam::refinementSurfaces::getNamedSurfaces() const
 435 {
 436    labelList namedSurfaces(cellZoneNames_.size());
 437
 438     label namedI = 0;
 439     forAll(cellZoneNames_, surfI)
 440     {
 441         if (cellZoneNames_[surfI].size() > 0)
 442         {
 443             namedSurfaces[namedI++] = surfI;
 444         }
 445     }
 446     namedSurfaces.setSize(namedI);
 447
 448     return namedSurfaces;
 449 }
 450
 451
 452 // Get indices of closed named surfaces
 453 Foam::labelList Foam::refinementSurfaces::getClosedNamedSurfaces() const
 454 {
 455     labelList named(getNamedSurfaces());
 456
 457     labelList closed(named.size());
 458     label closedI = 0;
 459
 460     forAll(named, i)
 461     {
 462         label surfI = named[i];
 463
 464         if (allGeometry_[surfaces_[surfI]].hasVolumeType())
 465         {
 466             closed[closedI++] = surfI;
 467         }
 468     }
 469     closed.setSize(closedI);
 470
 471     return closed;
 472 }
 473
 474
 475 // Orient surface (if they're closed) before any searching is done.
 476 void Foam::refinementSurfaces::orientSurface
 477 (
 478     const point& keepPoint,
 479     triSurfaceMesh& s
 480 )
 481 {
 482     // Flip surface so keepPoint is outside.
 483     bool anyFlipped = orientedSurface::orient(s, keepPoint, true);
 484
 485     if (anyFlipped)
 486     {
 487         // orientedSurface will have done a clearOut of the surface.
 488         // we could do a clearout of the triSurfaceMeshes::trees()
 489         // but these aren't affected by orientation (except for cached
 490         // sideness which should not be set at this point. !!Should
 491         // check!)
 492
 493         Info<< "orientSurfaces : Flipped orientation of surface "
 494             << s.searchableSurface::name()
 495             << " so point " << keepPoint << " is outside." << endl;
 496     }
 497 }
 498
 499
 500 // Count number of triangles per surface region
 501 Foam::labelList Foam::refinementSurfaces::countRegions(const triSurface& s)
 502 {
 503     const geometricSurfacePatchList& regions = s.patches();
 504
 505     labelList nTris(regions.size(), 0);
 506
 507     forAll(s, triI)
 508     {
 509         nTris[s[triI].region()]++;
 510     }
 511     return nTris;
 512 }
 513
 514
 515 // Precalculate the refinement level for every element of the searchable
 516 // surface. This is currently hardcoded for triSurfaceMesh only.
 517 void Foam::refinementSurfaces::setMinLevelFields
 518 (
 519     const shellSurfaces& shells
 520 )
 521 {
 522     minLevelFields_.setSize(surfaces_.size());
 523
 524     forAll(surfaces_, surfI)
 525     {
 526         const searchableSurface& geom = allGeometry_[surfaces_[surfI]];
 527
 528         if (isA<triSurfaceMesh>(geom))
 529         {
 530             const triSurfaceMesh& triMesh = refCast<const triSurfaceMesh>(geom);
 531
 532             minLevelFields_.set
 533             (
 534                 surfI,
 535                 new triSurfaceLabelField
 536                 (
 537                     IOobject
 538                     (
 539                         triMesh.searchableSurface::name(),
 540                         triMesh.objectRegistry::time().constant(),// directory
 541                         "triSurface",               // instance
 542                         triMesh,
 543                         IOobject::NO_READ,
 544                         IOobject::AUTO_WRITE,
 545                         false
 546                     ),
 547                     triMesh,
 548                     dimless
 549                 )
 550             );
 551             triSurfaceLabelField& minLevelField = minLevelFields_[surfI];
 552
 553             const triSurface& s = static_cast<const triSurface&>(triMesh);
 554
 555             // Initialise fields to region wise minLevel
 556             forAll(s, triI)
 557             {
 558                 minLevelField[triI] = minLevel(surfI, s[triI].region());
 559             }
 560
 561             // Find out if triangle inside shell with higher level
 562             pointField fc(s.size());
 563             forAll(s, triI)
 564             {
 565                 fc[triI] = s[triI].centre(s.points());
 566             }
 567             // What level does shell want to refine fc to?
 568             labelList shellLevel;
 569             shells.findHigherLevel(fc, minLevelField, shellLevel);
 570
 571             forAll(minLevelField, triI)
 572             {
 573                 minLevelField[triI] = max
 574                 (
 575                     minLevelField[triI],
 576                     shellLevel[triI]
 577                 );
 578             }
 579         }
 580     }
 581 }
 582
 583
 584 // Find intersections of edge. Return -1 or first surface with higher minLevel
 585 // number.
 586 void Foam::refinementSurfaces::findHigherIntersection
 587 (
 588     const pointField& start,
 589     const pointField& end,
 590     const labelList& currentLevel,   // current cell refinement level
 591
 592     labelList& surfaces,
 593     labelList& surfaceLevel
 594 ) const
 595 {
 596     surfaces.setSize(start.size());
 597     surfaces = -1;
 598     surfaceLevel.setSize(start.size());
 599     surfaceLevel = -1;
 600
 601     if (surfaces_.size() == 0)
 602     {
 603         return;
 604     }
 605
 606     // Work arrays
 607     labelList hitMap(identity(start.size()));
 608     pointField p0(start);
 609     pointField p1(end);
 610     List<pointIndexHit> hitInfo(start.size());
 611
 612
 613     forAll(surfaces_, surfI)
 614     {
 615         allGeometry_[surfaces_[surfI]].findLineAny(p0, p1, hitInfo);
 616
 617         // Copy all hits into arguments, continue with misses
 618         label newI = 0;
 619         forAll(hitInfo, hitI)
 620         {
 621             // Get the minLevel for the point
 622             label minLocalLevel = -1;
 623
 624             if (hitInfo[hitI].hit())
 625             {
 626                 // Check if minLevelField for this surface.
 627                 if
 628                 (
 629                     minLevelFields_.set(surfI)
 630                  && minLevelFields_[surfI].size() > 0
 631                 )
 632                 {
 633                     minLocalLevel =
 634                         minLevelFields_[surfI][hitInfo[hitI].index()];
 635                 }
 636                 else
 637                 {
 638                     // Use the min level for the surface instead. Assume
 639                     // single region 0.
 640                     minLocalLevel = minLevel(surfI, 0);
 641                 }
 642             }
 643
 644             label pointI = hitMap[hitI];
 645
 646             if (minLocalLevel > currentLevel[pointI])
 647             {
 648                 surfaces[pointI] = surfI;
 649                 surfaceLevel[pointI] = minLocalLevel;
 650             }
 651             else
 652             {
 653                 if (hitI != newI)
 654                 {
 655                     hitMap[newI] = hitMap[hitI];
 656                     p0[newI] = p0[hitI];
 657                     p1[newI] = p1[hitI];
 658                 }
 659                 newI++;
 660             }
 661         }
 662
 663         // All done? Note that this decision should be synchronised
 664         if (newI == 0)
 665         {
 666             break;
 667         }
 668
 669         // Trim and continue
 670         hitMap.setSize(newI);
 671         p0.setSize(newI);
 672         p1.setSize(newI);
 673         hitInfo.setSize(newI);
 674     }
 675 }
 676
 677
 678 void Foam::refinementSurfaces::findAllHigherIntersections
 679 (
 680     const pointField& start,
 681     const pointField& end,
 682     const labelList& currentLevel,   // current cell refinement level
 683
 684     List<vectorList>& surfaceNormal,
 685     labelListList& surfaceLevel
 686 ) const
 687 {
 688     surfaceLevel.setSize(start.size());
 689     surfaceNormal.setSize(start.size());
 690
 691     if (surfaces_.size() == 0)
 692     {
 693         return;
 694     }
 695
 696     // Work arrays
 697     List<List<pointIndexHit> > hitInfo;
 698     labelList pRegions;
 699     vectorField pNormals;
 700
 701     forAll(surfaces_, surfI)
 702     {
 703         allGeometry_[surfaces_[surfI]].findLineAll(start, end, hitInfo);
 704
 705         forAll(hitInfo, pointI)
 706         {
 707             const List<pointIndexHit>& pHits = hitInfo[pointI];
 708             allGeometry_[surfaces_[surfI]].getRegion(pHits, pRegions);
 709             allGeometry_[surfaces_[surfI]].getNormal(pHits, pNormals);
 710
 711             // Extract those hits that are on higher levelled surfaces.
 712             // Note: should move extraction of region, normal outside of loop
 713             // below for if getRegion/getNormal have high overhead.
 714
 715             forAll(pHits, pHitI)
 716             {
 717                 label region = globalRegion(surfI, pRegions[pHitI]);
 718
 719                 if (maxLevel_[region] > currentLevel[pointI])
 720                 {
 721                     // Append to pointI info
 722                     label sz = surfaceNormal[pointI].size();
 723                     surfaceNormal[pointI].setSize(sz+1);
 724                     surfaceNormal[pointI][sz] = pNormals[pHitI];
 725
 726                     surfaceLevel[pointI].setSize(sz+1);
 727                     surfaceLevel[pointI][sz] = maxLevel_[region];
 728                 }
 729             }
 730         }
 731     }
 732 }
 733
 734
 735 void Foam::refinementSurfaces::findNearestIntersection
 736 (
 737     const labelList& surfacesToTest,
 738     const pointField& start,
 739     const pointField& end,
 740
 741     labelList& surface1,
 742     List<pointIndexHit>& hit1,
 743     labelList& region1,
 744     labelList& surface2,
 745     List<pointIndexHit>& hit2,
 746     labelList& region2
 747 ) const
 748 {
 749     // 1. intersection from start to end
 750     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 751
 752     // Initialize arguments
 753     surface1.setSize(start.size());
 754     surface1 = -1;
 755     hit1.setSize(start.size());
 756     region1.setSize(start.size());
 757
 758     // Current end of segment to test.
 759     pointField nearest(end);
 760     // Work array
 761     List<pointIndexHit> nearestInfo(start.size());
 762     labelList region;
 763
 764     forAll(surfacesToTest, testI)
 765     {
 766         label surfI = surfacesToTest[testI];
 767
 768         // See if any intersection between start and current nearest
 769         allGeometry_[surfaces_[surfI]].findLine
 770         (
 771             start,
 772             nearest,
 773             nearestInfo
 774         );
 775         allGeometry_[surfaces_[surfI]].getRegion
 776         (
 777             nearestInfo,
 778             region
 779         );
 780
 781         forAll(nearestInfo, pointI)
 782         {
 783             if (nearestInfo[pointI].hit())
 784             {
 785                 hit1[pointI] = nearestInfo[pointI];
 786                 surface1[pointI] = surfI;
 787                 region1[pointI] = region[pointI];
 788                 nearest[pointI] = hit1[pointI].hitPoint();
 789             }
 790         }
 791     }
 792
 793
 794     // 2. intersection from end to last intersection
 795     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 796
 797     // Find the nearest intersection from end to start. Note that we initialize
 798     // to the first intersection (if any).
 799     surface2 = surface1;
 800     hit2 = hit1;
 801     region2 = region1;
 802
 803     // Set current end of segment to test.
 804     forAll(nearest, pointI)
 805     {
 806         if (hit1[pointI].hit())
 807         {
 808             nearest[pointI] = hit1[pointI].hitPoint();
 809         }
 810         else
 811         {
 812             // Disable testing by setting to end.
 813             nearest[pointI] = end[pointI];
 814         }
 815     }
 816
 817     forAll(surfacesToTest, testI)
 818     {
 819         label surfI = surfacesToTest[testI];
 820
 821         // See if any intersection between end and current nearest
 822         allGeometry_[surfaces_[surfI]].findLine
 823         (
 824             end,
 825             nearest,
 826             nearestInfo
 827         );
 828         allGeometry_[surfaces_[surfI]].getRegion
 829         (
 830             nearestInfo,
 831             region
 832         );
 833
 834         forAll(nearestInfo, pointI)
 835         {
 836             if (nearestInfo[pointI].hit())
 837             {
 838                 hit2[pointI] = nearestInfo[pointI];
 839                 surface2[pointI] = surfI;
 840                 region2[pointI] = region[pointI];
 841                 nearest[pointI] = hit2[pointI].hitPoint();
 842             }
 843         }
 844     }
 845
 846
 847     // Make sure that if hit1 has hit something, hit2 will have at least the
 848     // same point (due to tolerances it might miss its end point)
 849     forAll(hit1, pointI)
 850     {
 851         if (hit1[pointI].hit() && !hit2[pointI].hit())
 852         {
 853             hit2[pointI] = hit1[pointI];
 854             surface2[pointI] = surface1[pointI];
 855             region2[pointI] = region1[pointI];
 856         }
 857     }
 858 }
 859
 860
 861 void Foam::refinementSurfaces::findAnyIntersection
 862 (
 863     const pointField& start,
 864     const pointField& end,
 865
 866     labelList& hitSurface,
 867     List<pointIndexHit>& hitInfo
 868 ) const
 869 {
 870     searchableSurfacesQueries::findAnyIntersection
 871     (
 872         allGeometry_,
 873         surfaces_,
 874         start,
 875         end,
 876         hitSurface,
 877         hitInfo
 878     );
 879 }
 880
 881
 882 void Foam::refinementSurfaces::findNearest
 883 (
 884     const labelList& surfacesToTest,
 885     const pointField& samples,
 886     const  scalarField& nearestDistSqr,
 887     labelList& hitSurface,
 888     List<pointIndexHit>& hitInfo
 889 ) const
 890 {
 891     labelList geometries(IndirectList<label>(surfaces_, surfacesToTest));
 892
 893     // Do the tests. Note that findNearest returns index in geometries.
 894     searchableSurfacesQueries::findNearest
 895     (
 896         allGeometry_,
 897         geometries,
 898         samples,
 899         nearestDistSqr,
 900         hitSurface,
 901         hitInfo
 902     );
 903
 904     // Rework the hitSurface to be surface (i.e. index into surfaces_)
 905     forAll(hitSurface, pointI)
 906     {
 907         if (hitSurface[pointI] != -1)
 908         {
 909             hitSurface[pointI] = surfacesToTest[hitSurface[pointI]];
 910         }
 911     }
 912 }
 913
 914
 915 //// Find intersection with max of edge. Return -1 or the surface
 916 //// with the highest maxLevel above currentLevel
 917 //Foam::label Foam::refinementSurfaces::findHighestIntersection
 918 //(
 919 //    const point& start,
 920 //    const point& end,
 921 //    const label currentLevel,   // current cell refinement level
 922 //
 923 //    pointIndexHit& maxHit
 924 //) const
 925 //{
 926 //    // surface with highest maxlevel
 927 //    label maxSurface = -1;
 928 //    // maxLevel of maxSurface
 929 //    label maxLevel = currentLevel;
 930 //
 931 //    forAll(*this, surfI)
 932 //    {
 933 //        pointIndexHit hit = operator[](surfI).findLineAny(start, end);
 934 //
 935 //        if (hit.hit())
 936 //        {
 937 //            const triSurface& s = operator[](surfI);
 938 //
 939 //            label region = globalRegion(surfI, s[hit.index()].region());
 940 //
 941 //            if (maxLevel_[region] > maxLevel)
 942 //            {
 943 //                maxSurface = surfI;
 944 //                maxLevel = maxLevel_[region];
 945 //                maxHit = hit;
 946 //            }
 947 //        }
 948 //    }
 949 //
 950 //    if (maxSurface == -1)
 951 //    {
 952 //        // maxLevel unchanged. No interesting surface hit.
 953 //        maxHit.setMiss();
 954 //    }
 955 //
 956 //    return maxSurface;
 957 //}
 958
 959
 960 void Foam::refinementSurfaces::findInside
 961 (
 962     const labelList& testSurfaces,
 963     const pointField& pt,
 964     labelList& insideSurfaces
 965 ) const
 966 {
 967     insideSurfaces.setSize(pt.size());
 968     insideSurfaces = -1;
 969
 970     forAll(testSurfaces, i)
 971     {
 972         label surfI = testSurfaces[i];
 973
 974         if (allGeometry_[surfaces_[surfI]].hasVolumeType())
 975         {
 976             List<searchableSurface::volumeType> volType;
 977             allGeometry_[surfaces_[surfI]].getVolumeType(pt, volType);
 978
 979             forAll(volType, pointI)
 980             {
 981                 if (insideSurfaces[pointI] == -1)
 982                 {
 983                     if
 984                     (
 985                         (
 986                             volType[pointI] == triSurfaceMesh::INSIDE
 987                          && zoneInside_[surfI]
 988                         )
 989                      || (
 990                             volType[pointI] == triSurfaceMesh::OUTSIDE
 991                          && !zoneInside_[surfI]
 992                         )
 993                     )
 994                     {
 995                         insideSurfaces[pointI] = surfI;
 996                     }
 997                 }
 998             }
 999         }
1000     }
1001 }
1002
1003
1004 // ************************************************************************* //