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