src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C

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