src/meshTools/searchableSurface/triSurfaceMesh.C

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