src/meshTools/searchableSurface/triSurfaceMesh.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   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 "triSurfaceMesh.H"
  27 #include "Random.H"
  28 #include "addToRunTimeSelectionTable.H"
  29 #include "EdgeMap.H"
  30 #include "triSurfaceFields.H"
  31 #include "Time.H"
  32 #include "PackedBoolList.H"
  33
  34 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  35
  36 namespace Foam
  37 {
  38
  39 defineTypeNameAndDebug(triSurfaceMesh, 0);
  40 addToRunTimeSelectionTable(searchableSurface, triSurfaceMesh, dict);
  41
  42 }
  43
  44 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  45
  46 //// Special version of Time::findInstance that does not check headerOk
  47 //// to search for instances of raw files
  48 //Foam::word Foam::triSurfaceMesh::findRawInstance
  49 //(
  50 //    const Time& runTime,
  51 //    const fileName& dir,
  52 //    const word& name
  53 //)
  54 //{
  55 //    // Check current time first
  56 //    if (isFile(runTime.path()/runTime.timeName()/dir/name))
  57 //    {
  58 //        return runTime.timeName();
  59 //    }
  60 //    instantList ts = runTime.times();
  61 //    label instanceI;
  62 //
  63 //    for (instanceI = ts.size()-1; instanceI >= 0; --instanceI)
  64 //    {
  65 //        if (ts[instanceI].value() <= runTime.timeOutputValue())
  66 //        {
  67 //            break;
  68 //        }
  69 //    }
  70 //
  71 //    // continue searching from here
  72 //    for (; instanceI >= 0; --instanceI)
  73 //    {
  74 //        if (isFile(runTime.path()/ts[instanceI].name()/dir/name))
  75 //        {
  76 //            return ts[instanceI].name();
  77 //        }
  78 //    }
  79 //
  80 //
  81 //    // not in any of the time directories, try constant
  82 //
  83 //    // Note. This needs to be a hard-coded constant, rather than the
  84 //    // constant function of the time, because the latter points to
  85 //    // the case constant directory in parallel cases
  86 //
  87 //    if (isFile(runTime.path()/runTime.constant()/dir/name))
  88 //    {
  89 //        return runTime.constant();
  90 //    }
  91 //
  92 //    FatalErrorIn
  93 //    (
  94 //        "searchableSurfaces::findRawInstance"
  95 //        "(const Time&, const fileName&, const word&)"
  96 //    )   << "Cannot find file \"" << name << "\" in directory "
  97 //        << runTime.constant()/dir
  98 //        << exit(FatalError);
  99 //
 100 //    return runTime.constant();
 101 //}
 102
 103
 104 //- Check file existence
 105 const Foam::fileName& Foam::triSurfaceMesh::checkFile
 106 (
 107     const fileName& fName,
 108     const fileName& objectName
 109 )
 110 {
 111     if (fName.empty())
 112     {
 113         FatalErrorIn
 114         (
 115             "triSurfaceMesh::checkFile(const fileName&, const fileName&)"
 116         )   << "Cannot find triSurfaceMesh starting from "
 117             << objectName << exit(FatalError);
 118     }
 119     return fName;
 120 }
 121
 122
 123 bool Foam::triSurfaceMesh::addFaceToEdge
 124 (
 125     const edge& e,
 126     EdgeMap<label>& facesPerEdge
 127 )
 128 {
 129     EdgeMap<label>::iterator eFnd = facesPerEdge.find(e);
 130     if (eFnd != facesPerEdge.end())
 131     {
 132         if (eFnd() == 2)
 133         {
 134             return false;
 135         }
 136         eFnd()++;
 137     }
 138     else
 139     {
 140         facesPerEdge.insert(e, 1);
 141     }
 142     return true;
 143 }
 144
 145
 146 bool Foam::triSurfaceMesh::isSurfaceClosed() const
 147 {
 148     // Construct pointFaces. Let's hope surface has compact point
 149     // numbering ...
 150     labelListList pointFaces;
 151     invertManyToMany(points().size(), *this, pointFaces);
 152
 153     // Loop over all faces surrounding point. Count edges emanating from point.
 154     // Every edge should be used by two faces exactly.
 155     // To prevent doing work twice per edge only look at edges to higher
 156     // point
 157     EdgeMap<label> facesPerEdge(100);
 158     forAll(pointFaces, pointI)
 159     {
 160         const labelList& pFaces = pointFaces[pointI];
 161
 162         facesPerEdge.clear();
 163         forAll(pFaces, i)
 164         {
 165             const triSurface::FaceType& f = triSurface::operator[](pFaces[i]);
 166             label fp = findIndex(f, pointI);
 167
 168             // Something weird: if I expand the code of addFaceToEdge in both
 169             // below instances it gives a segmentation violation on some
 170             // surfaces. Compiler (4.3.2) problem?
 171
 172
 173             // Forward edge
 174             label nextPointI = f[f.fcIndex(fp)];
 175
 176             if (nextPointI > pointI)
 177             {
 178                 bool okFace = addFaceToEdge
 179                 (
 180                     edge(pointI, nextPointI),
 181                     facesPerEdge
 182                 );
 183
 184                 if (!okFace)
 185                 {
 186                     return false;
 187                 }
 188             }
 189             // Reverse edge
 190             label prevPointI = f[f.rcIndex(fp)];
 191
 192             if (prevPointI > pointI)
 193             {
 194                 bool okFace = addFaceToEdge
 195                 (
 196                     edge(pointI, prevPointI),
 197                     facesPerEdge
 198                 );
 199
 200                 if (!okFace)
 201                 {
 202                     return false;
 203                 }
 204             }
 205         }
 206
 207         // Check for any edges used only once.
 208         forAllConstIter(EdgeMap<label>, facesPerEdge, iter)
 209         {
 210             if (iter() != 2)
 211             {
 212                 return false;
 213             }
 214         }
 215     }
 216
 217     return true;
 218 }
 219
 220
 221 // Gets all intersections after initial one. Adds smallVec and starts tracking
 222 // from there.
 223 void Foam::triSurfaceMesh::getNextIntersections
 224 (
 225     const indexedOctree<treeDataTriSurface>& octree,
 226     const point& start,
 227     const point& end,
 228     const vector& smallVec,
 229     DynamicList<pointIndexHit, 1, 1>& hits
 230 )
 231 {
 232     const vector dirVec(end-start);
 233     const scalar magSqrDirVec(magSqr(dirVec));
 234
 235     // Initial perturbation amount
 236     vector perturbVec(smallVec);
 237
 238     while (true)
 239     {
 240         // Start tracking from last hit.
 241         point pt = hits.last().hitPoint() + perturbVec;
 242
 243         if (((pt-start)&dirVec) > magSqrDirVec)
 244         {
 245             return;
 246         }
 247
 248         // See if any intersection between pt and end
 249         pointIndexHit inter = octree.findLine(pt, end);
 250
 251         if (!inter.hit())
 252         {
 253             return;
 254         }
 255
 256         // Check if already found this intersection
 257         bool duplicateHit = false;
 258         forAllReverse(hits, i)
 259         {
 260             if (hits[i].index() == inter.index())
 261             {
 262                 duplicateHit = true;
 263                 break;
 264             }
 265         }
 266
 267
 268         if (duplicateHit)
 269         {
 270             // Hit same triangle again. Increase perturbVec and try again.
 271             perturbVec *= 2;
 272         }
 273         else
 274         {
 275             // Proper hit
 276             hits.append(inter);
 277             // Restore perturbVec
 278             perturbVec = smallVec;
 279         }
 280     }
 281 }
 282
 283
 284 void Foam::triSurfaceMesh::calcBounds(boundBox& bb, label& nPoints) const
 285 {
 286     // Unfortunately nPoints constructs meshPoints() so do compact version
 287     // ourselves
 288
 289     const triSurface& s = static_cast<const triSurface&>(*this);
 290
 291     PackedBoolList pointIsUsed(points().size());
 292
 293     nPoints = 0;
 294     bb = boundBox::invertedBox;
 295
 296     forAll(s, faceI)
 297     {
 298         const triSurface::FaceType& f = s[faceI];
 299
 300         forAll(f, fp)
 301         {
 302             label pointI = f[fp];
 303             if (pointIsUsed.set(pointI, 1u))
 304             {
 305                 bb.min() = ::Foam::min(bb.min(), points()[pointI]);
 306                 bb.max() = ::Foam::max(bb.max(), points()[pointI]);
 307                 nPoints++;
 308             }
 309         }
 310     }
 311 }
 312
 313
 314 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 315
 316 Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io, const triSurface& s)
 317 :
 318     searchableSurface(io),
 319     objectRegistry
 320     (
 321         IOobject
 322         (
 323             io.name(),
 324             io.instance(),
 325             io.local(),
 326             io.db(),
 327             io.readOpt(),
 328             io.writeOpt(),
 329             false       // searchableSurface already registered under name
 330         )
 331     ),
 332     triSurface(s),
 333     tolerance_(indexedOctree<treeDataTriSurface>::perturbTol()),
 334     minQuality_(-1),
 335     maxTreeDepth_(10),
 336     surfaceClosed_(-1)
 337 {}
 338
 339
 340 Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io)
 341 :
 342     // Find instance for triSurfaceMesh
 343     searchableSurface(io),
 344     //(
 345     //    IOobject
 346     //    (
 347     //        io.name(),
 348     //        io.time().findInstance(io.local(), word::null),
 349     //        io.local(),
 350     //        io.db(),
 351     //        io.readOpt(),
 352     //        io.writeOpt(),
 353     //        io.registerObject()
 354     //    )
 355     //),
 356     // Reused found instance in objectRegistry
 357     objectRegistry
 358     (
 359         IOobject
 360         (
 361             io.name(),
 362             static_cast<const searchableSurface&>(*this).instance(),
 363             io.local(),
 364             io.db(),
 365             io.readOpt(),
 366             io.writeOpt(),
 367             false       // searchableSurface already registered under name
 368         )
 369     ),
 370     triSurface
 371     (
 372         checkFile
 373         (
 374             searchableSurface::filePath(),
 375             searchableSurface::objectPath()
 376         )
 377     ),
 378     tolerance_(indexedOctree<treeDataTriSurface>::perturbTol()),
 379     minQuality_(-1),
 380     maxTreeDepth_(10),
 381     surfaceClosed_(-1)
 382 {}
 383
 384
 385 Foam::triSurfaceMesh::triSurfaceMesh
 386 (
 387     const IOobject& io,
 388     const dictionary& dict
 389 )
 390 :
 391     searchableSurface(io),
 392     //(
 393     //    IOobject
 394     //    (
 395     //        io.name(),
 396     //        io.time().findInstance(io.local(), word::null),
 397     //        io.local(),
 398     //        io.db(),
 399     //        io.readOpt(),
 400     //        io.writeOpt(),
 401     //        io.registerObject()
 402     //    )
 403     //),
 404     // Reused found instance in objectRegistry
 405     objectRegistry
 406     (
 407         IOobject
 408         (
 409             io.name(),
 410             static_cast<const searchableSurface&>(*this).instance(),
 411             io.local(),
 412             io.db(),
 413             io.readOpt(),
 414             io.writeOpt(),
 415             false       // searchableSurface already registered under name
 416         )
 417     ),
 418     triSurface
 419     (
 420         checkFile
 421         (
 422             searchableSurface::filePath(),
 423             searchableSurface::objectPath()
 424         )
 425     ),
 426     tolerance_(indexedOctree<treeDataTriSurface>::perturbTol()),
 427     minQuality_(-1),
 428     maxTreeDepth_(10),
 429     surfaceClosed_(-1)
 430 {
 431     scalar scaleFactor = 0;
 432
 433     // allow rescaling of the surface points
 434     // eg, CAD geometries are often done in millimeters
 435     if (dict.readIfPresent("scale", scaleFactor) && scaleFactor > 0)
 436     {
 437         Info<< searchableSurface::name() << " : using scale " << scaleFactor
 438             << endl;
 439         triSurface::scalePoints(scaleFactor);
 440     }
 441
 442
 443     // Have optional non-standard search tolerance for gappy surfaces.
 444     if (dict.readIfPresent("tolerance", tolerance_) && tolerance_ > 0)
 445     {
 446         Info<< searchableSurface::name() << " : using intersection tolerance "
 447             << tolerance_ << endl;
 448     }
 449
 450     // Have optional minimum quality for normal calculation
 451     if (dict.readIfPresent("minQuality", minQuality_) && minQuality_ > 0)
 452     {
 453         Info<< searchableSurface::name()
 454             << " : ignoring triangles with quality < "
 455             << minQuality_ << " for normals calculation." << endl;
 456     }
 457
 458     // Have optional non-standard tree-depth to limit storage.
 459     if (dict.readIfPresent("maxTreeDepth", maxTreeDepth_) && maxTreeDepth_ > 0)
 460     {
 461         Info<< searchableSurface::name() << " : using maximum tree depth "
 462             << maxTreeDepth_ << endl;
 463     }
 464 }
 465
 466
 467 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 468
 469 Foam::triSurfaceMesh::~triSurfaceMesh()
 470 {
 471     clearOut();
 472 }
 473
 474
 475 void Foam::triSurfaceMesh::clearOut()
 476 {
 477     tree_.clear();
 478     edgeTree_.clear();
 479     triSurface::clearOut();
 480 }
 481
 482
 483 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 484
 485 Foam::pointField Foam::triSurfaceMesh::coordinates() const
 486 {
 487     // Use copy to calculate face centres so they don't get stored
 488     return PrimitivePatch<triSurface::FaceType, SubList, const pointField&>
 489     (
 490         SubList<triSurface::FaceType>(*this, triSurface::size()),
 491         triSurface::points()
 492     ).faceCentres();
 493 }
 494
 495
 496 void Foam::triSurfaceMesh::movePoints(const pointField& newPoints)
 497 {
 498     tree_.clear();
 499     edgeTree_.clear();
 500     triSurface::movePoints(newPoints);
 501 }
 502
 503
 504 const Foam::indexedOctree<Foam::treeDataTriSurface>&
 505 Foam::triSurfaceMesh::tree() const
 506 {
 507     if (tree_.empty())
 508     {
 509         // Calculate bb without constructing local point numbering.
 510         treeBoundBox bb;
 511         label nPoints;
 512         calcBounds(bb, nPoints);
 513
 514         if (nPoints != points().size())
 515         {
 516             WarningIn("triSurfaceMesh::tree() const")
 517                 << "Surface " << searchableSurface::name()
 518                 << " does not have compact point numbering."
 519                 << " Of " << points().size() << " only " << nPoints
 520                 << " are used. This might give problems in some routines."
 521                 << endl;
 522         }
 523
 524
 525         // Random number generator. Bit dodgy since not exactly random ;-)
 526         Random rndGen(65431);
 527
 528         // Slightly extended bb. Slightly off-centred just so on symmetric
 529         // geometry there are less face/edge aligned items.
 530         bb = bb.extend(rndGen, 1E-4);
 531         bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
 532         bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
 533
 534         scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
 535         indexedOctree<treeDataTriSurface>::perturbTol() = tolerance_;
 536
 537         tree_.reset
 538         (
 539             new indexedOctree<treeDataTriSurface>
 540             (
 541                 treeDataTriSurface(*this, tolerance_),
 542                 bb,
 543                 maxTreeDepth_,  // maxLevel
 544                 10,             // leafsize
 545                 3.0             // duplicity
 546             )
 547         );
 548
 549         indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
 550     }
 551
 552     return tree_();
 553 }
 554
 555
 556 const Foam::indexedOctree<Foam::treeDataEdge>&
 557 Foam::triSurfaceMesh::edgeTree() const
 558 {
 559     if (edgeTree_.empty())
 560     {
 561         // Boundary edges
 562         labelList bEdges
 563         (
 564             identity
 565             (
 566                 nEdges()
 567                -nInternalEdges()
 568             )
 569           + nInternalEdges()
 570         );
 571
 572         treeBoundBox bb;
 573         label nPoints;
 574         calcBounds(bb, nPoints);
 575
 576         // Random number generator. Bit dodgy since not exactly random ;-)
 577         Random rndGen(65431);
 578
 579         // Slightly extended bb. Slightly off-centred just so on symmetric
 580         // geometry there are less face/edge aligned items.
 581
 582         bb = bb.extend(rndGen, 1E-4);
 583         bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
 584         bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
 585
 586         scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
 587         indexedOctree<treeDataEdge>::perturbTol() = tolerance_;
 588
 589         edgeTree_.reset
 590         (
 591             new indexedOctree<treeDataEdge>
 592             (
 593                 treeDataEdge
 594                 (
 595                     false,          // cachebb
 596                     edges(),        // edges
 597                     localPoints(),  // points
 598                     bEdges          // selected edges
 599                 ),
 600                 bb,                 // bb
 601                 maxTreeDepth_,      // maxLevel
 602                 10,                 // leafsize
 603                 3.0                 // duplicity
 604             )
 605         );
 606
 607         indexedOctree<treeDataEdge>::perturbTol() = oldTol;
 608     }
 609     return edgeTree_();
 610 }
 611
 612
 613 const Foam::wordList& Foam::triSurfaceMesh::regions() const
 614 {
 615     if (regions_.empty())
 616     {
 617         regions_.setSize(patches().size());
 618         forAll(regions_, regionI)
 619         {
 620             regions_[regionI] = patches()[regionI].name();
 621         }
 622     }
 623     return regions_;
 624 }
 625
 626
 627 // Find out if surface is closed.
 628 bool Foam::triSurfaceMesh::hasVolumeType() const
 629 {
 630     if (surfaceClosed_ == -1)
 631     {
 632         if (isSurfaceClosed())
 633         {
 634             surfaceClosed_ = 1;
 635         }
 636         else
 637         {
 638             surfaceClosed_ = 0;
 639         }
 640     }
 641
 642     return surfaceClosed_ == 1;
 643 }
 644
 645
 646 void Foam::triSurfaceMesh::findNearest
 647 (
 648     const pointField& samples,
 649     const scalarField& nearestDistSqr,
 650     List<pointIndexHit>& info
 651 ) const
 652 {
 653     const indexedOctree<treeDataTriSurface>& octree = tree();
 654
 655     info.setSize(samples.size());
 656
 657     scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
 658     indexedOctree<treeDataTriSurface>::perturbTol() = tolerance_;
 659
 660     forAll(samples, i)
 661     {
 662         static_cast<pointIndexHit&>(info[i]) = octree.findNearest
 663         (
 664             samples[i],
 665             nearestDistSqr[i]
 666         );
 667     }
 668
 669     indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
 670 }
 671
 672
 673 void Foam::triSurfaceMesh::findLine
 674 (
 675     const pointField& start,
 676     const pointField& end,
 677     List<pointIndexHit>& info
 678 ) const
 679 {
 680     const indexedOctree<treeDataTriSurface>& octree = tree();
 681
 682     info.setSize(start.size());
 683
 684     scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
 685     indexedOctree<treeDataTriSurface>::perturbTol() = tolerance_;
 686
 687     forAll(start, i)
 688     {
 689         static_cast<pointIndexHit&>(info[i]) = octree.findLine
 690         (
 691             start[i],
 692             end[i]
 693         );
 694     }
 695
 696     indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
 697 }
 698
 699
 700 void Foam::triSurfaceMesh::findLineAny
 701 (
 702     const pointField& start,
 703     const pointField& end,
 704     List<pointIndexHit>& info
 705 ) const
 706 {
 707     const indexedOctree<treeDataTriSurface>& octree = tree();
 708
 709     info.setSize(start.size());
 710
 711     scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
 712     indexedOctree<treeDataTriSurface>::perturbTol() = tolerance_;
 713
 714     forAll(start, i)
 715     {
 716         static_cast<pointIndexHit&>(info[i]) = octree.findLineAny
 717         (
 718             start[i],
 719             end[i]
 720         );
 721     }
 722
 723     indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
 724 }
 725
 726
 727 void Foam::triSurfaceMesh::findLineAll
 728 (
 729     const pointField& start,
 730     const pointField& end,
 731     List<List<pointIndexHit> >& info
 732 ) const
 733 {
 734     const indexedOctree<treeDataTriSurface>& octree = tree();
 735
 736     info.setSize(start.size());
 737
 738     scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
 739     indexedOctree<treeDataTriSurface>::perturbTol() = tolerance_;
 740
 741     // Work array
 742     DynamicList<pointIndexHit, 1, 1> hits;
 743
 744     // Tolerances:
 745     // To find all intersections we add a small vector to the last intersection
 746     // This is chosen such that
 747     // - it is significant (SMALL is smallest representative relative tolerance;
 748     //   we need something bigger since we're doing calculations)
 749     // - if the start-end vector is zero we still progress
 750     const vectorField dirVec(end-start);
 751     const vectorField smallVec
 752     (
 753         indexedOctree<treeDataTriSurface>::perturbTol()*dirVec
 754       + vector(ROOTVSMALL,ROOTVSMALL,ROOTVSMALL)
 755     );
 756
 757     forAll(start, pointI)
 758     {
 759         // See if any intersection between pt and end
 760         pointIndexHit inter = octree.findLine(start[pointI], end[pointI]);
 761
 762         if (inter.hit())
 763         {
 764             hits.clear();
 765             hits.append(inter);
 766
 767             getNextIntersections
 768             (
 769                 octree,
 770                 start[pointI],
 771                 end[pointI],
 772                 smallVec[pointI],
 773                 hits
 774             );
 775
 776             info[pointI].transfer(hits);
 777         }
 778         else
 779         {
 780             info[pointI].clear();
 781         }
 782     }
 783
 784     indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
 785 }
 786
 787
 788 void Foam::triSurfaceMesh::getRegion
 789 (
 790     const List<pointIndexHit>& info,
 791     labelList& region
 792 ) const
 793 {
 794     region.setSize(info.size());
 795     forAll(info, i)
 796     {
 797         if (info[i].hit())
 798         {
 799             region[i] = triSurface::operator[](info[i].index()).region();
 800         }
 801         else
 802         {
 803             region[i] = -1;
 804         }
 805     }
 806 }
 807
 808
 809 void Foam::triSurfaceMesh::getNormal
 810 (
 811     const List<pointIndexHit>& info,
 812     vectorField& normal
 813 ) const
 814 {
 815     normal.setSize(info.size());
 816
 817     if (minQuality_ >= 0)
 818     {
 819         // Make sure we don't use triangles with low quality since
 820         // normal is not reliable.
 821
 822         const triSurface& s = static_cast<const triSurface&>(*this);
 823         const labelListList& faceFaces = s.faceFaces();
 824
 825         forAll(info, i)
 826         {
 827             if (info[i].hit())
 828             {
 829                 label faceI = info[i].index();
 830
 831                 scalar qual = s[faceI].tri(points()).quality();
 832
 833                 if (qual < minQuality_)
 834                 {
 835                     // Search neighbouring triangles
 836                     const labelList& fFaces = faceFaces[faceI];
 837
 838                     forAll(fFaces, j)
 839                     {
 840                         label nbrI = fFaces[j];
 841                         scalar nbrQual = s[nbrI].tri(points()).quality();
 842                         if (nbrQual > qual)
 843                         {
 844                             qual = nbrQual;
 845                             normal[i] = s[nbrI].normal(points());
 846                         }
 847                     }
 848                 }
 849                 else
 850                 {
 851                     normal[i] = s[faceI].normal(points());
 852                 }
 853                 normal[i] /= mag(normal[i]);
 854             }
 855             else
 856             {
 857                 // Set to what?
 858                 normal[i] = vector::zero;
 859             }
 860         }
 861     }
 862     else
 863     {
 864         forAll(info, i)
 865         {
 866             if (info[i].hit())
 867             {
 868                 label faceI = info[i].index();
 869                 //- Cached:
 870                 //normal[i] = faceNormals()[faceI];
 871
 872                 //- Uncached
 873                 normal[i] = triSurface::operator[](faceI).normal(points());
 874                 normal[i] /= mag(normal[i]) + VSMALL;
 875             }
 876             else
 877             {
 878                 // Set to what?
 879                 normal[i] = vector::zero;
 880             }
 881         }
 882     }
 883 }
 884
 885
 886 void Foam::triSurfaceMesh::setField(const labelList& values)
 887 {
 888     autoPtr<triSurfaceLabelField> fldPtr
 889     (
 890         new triSurfaceLabelField
 891         (
 892             IOobject
 893             (
 894                 "values",
 895                 objectRegistry::time().timeName(),  // instance
 896                 "triSurface",                       // local
 897                 *this,
 898                 IOobject::NO_READ,
 899                 IOobject::AUTO_WRITE
 900             ),
 901             *this,
 902             dimless,
 903             labelField(values)
 904         )
 905     );
 906
 907     // Store field on triMesh
 908     fldPtr.ptr()->store();
 909 }
 910
 911
 912 void Foam::triSurfaceMesh::getField
 913 (
 914     const List<pointIndexHit>& info,
 915     labelList& values
 916 ) const
 917 {
 918     if (foundObject<triSurfaceLabelField>("values"))
 919     {
 920         values.setSize(info.size());
 921
 922         const triSurfaceLabelField& fld = lookupObject<triSurfaceLabelField>
 923         (
 924             "values"
 925         );
 926
 927         forAll(info, i)
 928         {
 929             if (info[i].hit())
 930             {
 931                 values[i] = fld[info[i].index()];
 932             }
 933         }
 934     }
 935 }
 936
 937
 938 void Foam::triSurfaceMesh::getVolumeType
 939 (
 940     const pointField& points,
 941     List<volumeType>& volType
 942 ) const
 943 {
 944     volType.setSize(points.size());
 945
 946     scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
 947     indexedOctree<treeDataTriSurface>::perturbTol() = tolerance_;
 948
 949     forAll(points, pointI)
 950     {
 951         const point& pt = points[pointI];
 952
 953         // - use cached volume type per each tree node
 954         // - cheat conversion since same values
 955         volType[pointI] = static_cast<volumeType>(tree().getVolumeType(pt));
 956     }
 957
 958     indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
 959 }
 960
 961
 962 //- Write using given format, version and compression
 963 bool Foam::triSurfaceMesh::writeObject
 964 (
 965     IOstream::streamFormat fmt,
 966     IOstream::versionNumber ver,
 967     IOstream::compressionType cmp
 968 ) const
 969 {
 970     fileName fullPath(searchableSurface::objectPath());
 971
 972     if (!mkDir(fullPath.path()))
 973     {
 974         return false;
 975     }
 976
 977     triSurface::write(fullPath);
 978
 979     if (!isFile(fullPath))
 980     {
 981         return false;
 982     }
 983
 984     //return objectRegistry::writeObject(fmt, ver, cmp);
 985     return true;
 986 }
 987
 988
 989 // ************************************************************************* //