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