src/meshTools/searchableSurface/triSurfaceMesh.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software; you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM; if not, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*---------------------------------------------------------------------------*/
  26
  27 #include "triSurfaceMesh.H"
  28 #include "Random.H"
  29 #include "addToRunTimeSelectionTable.H"
  30 #include "EdgeMap.H"
  31 #include "triSurfaceFields.H"
  32 #include "Time.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 labelledTri& 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 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 222
 223 Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io, const triSurface& s)
 224 :
 225     searchableSurface(io),
 226     objectRegistry
 227     (
 228         IOobject
 229         (
 230             io.name(),
 231             io.instance(),
 232             io.local(),
 233             io.db(),
 234             io.readOpt(),
 235             io.writeOpt(),
 236             false       // searchableSurface already registered under name
 237         )
 238     ),
 239     triSurface(s),
 240     tolerance_(indexedOctree<treeDataTriSurface>::perturbTol()),
 241     surfaceClosed_(-1)
 242 {}
 243
 244
 245 Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io)
 246 :
 247     // Find instance for triSurfaceMesh
 248     searchableSurface(io),
 249     //(
 250     //    IOobject
 251     //    (
 252     //        io.name(),
 253     //        io.time().findInstance(io.local(), word::null),
 254     //        io.local(),
 255     //        io.db(),
 256     //        io.readOpt(),
 257     //        io.writeOpt(),
 258     //        io.registerObject()
 259     //    )
 260     //),
 261     // Reused found instance in objectRegistry
 262     objectRegistry
 263     (
 264         IOobject
 265         (
 266             io.name(),
 267             static_cast<const searchableSurface&>(*this).instance(),
 268             io.local(),
 269             io.db(),
 270             io.readOpt(),
 271             io.writeOpt(),
 272             false       // searchableSurface already registered under name
 273         )
 274     ),
 275     triSurface
 276     (
 277         checkFile
 278         (
 279             searchableSurface::filePath(),
 280             searchableSurface::objectPath()
 281         )
 282     ),
 283     tolerance_(indexedOctree<treeDataTriSurface>::perturbTol()),
 284     surfaceClosed_(-1)
 285 {}
 286
 287
 288 Foam::triSurfaceMesh::triSurfaceMesh
 289 (
 290     const IOobject& io,
 291     const dictionary& dict
 292 )
 293 :
 294     searchableSurface(io),
 295     //(
 296     //    IOobject
 297     //    (
 298     //        io.name(),
 299     //        io.time().findInstance(io.local(), word::null),
 300     //        io.local(),
 301     //        io.db(),
 302     //        io.readOpt(),
 303     //        io.writeOpt(),
 304     //        io.registerObject()
 305     //    )
 306     //),
 307     // Reused found instance in objectRegistry
 308     objectRegistry
 309     (
 310         IOobject
 311         (
 312             io.name(),
 313             static_cast<const searchableSurface&>(*this).instance(),
 314             io.local(),
 315             io.db(),
 316             io.readOpt(),
 317             io.writeOpt(),
 318             false       // searchableSurface already registered under name
 319         )
 320     ),
 321     triSurface
 322     (
 323         checkFile
 324         (
 325             searchableSurface::filePath(),
 326             searchableSurface::objectPath()
 327         )
 328     ),
 329     tolerance_(indexedOctree<treeDataTriSurface>::perturbTol()),
 330     surfaceClosed_(-1)
 331 {
 332     scalar scaleFactor = 0;
 333
 334     // allow rescaling of the surface points
 335     // eg, CAD geometries are often done in millimeters
 336     if (dict.readIfPresent("scale", scaleFactor) && scaleFactor > 0)
 337     {
 338         Info<< searchableSurface::name() << " : using scale " << scaleFactor
 339             << endl;
 340         triSurface::scalePoints(scaleFactor);
 341     }
 342
 343
 344     // Have optional non-standard search tolerance for gappy surfaces.
 345     if (dict.readIfPresent("tolerance", tolerance_) && tolerance_ > 0)
 346     {
 347         Info<< searchableSurface::name() << " : using intersection tolerance "
 348             << tolerance_ << endl;
 349     }
 350 }
 351
 352
 353 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 354
 355 Foam::triSurfaceMesh::~triSurfaceMesh()
 356 {
 357     clearOut();
 358 }
 359
 360
 361 void Foam::triSurfaceMesh::clearOut()
 362 {
 363     tree_.clear();
 364     edgeTree_.clear();
 365     triSurface::clearOut();
 366 }
 367
 368
 369 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 370
 371 void Foam::triSurfaceMesh::movePoints(const pointField& newPoints)
 372 {
 373     tree_.clear();
 374     edgeTree_.clear();
 375     triSurface::movePoints(newPoints);
 376 }
 377
 378
 379 const Foam::indexedOctree<Foam::treeDataTriSurface>&
 380     Foam::triSurfaceMesh::tree() const
 381 {
 382     if (tree_.empty())
 383     {
 384         // Random number generator. Bit dodgy since not exactly random ;-)
 385         Random rndGen(65431);
 386
 387         // Slightly extended bb. Slightly off-centred just so on symmetric
 388         // geometry there are less face/edge aligned items.
 389         treeBoundBox bb
 390         (
 391             treeBoundBox(points(), meshPoints()).extend(rndGen, 1E-4)
 392         );
 393         bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
 394         bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
 395
 396         scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
 397         indexedOctree<treeDataTriSurface>::perturbTol() = tolerance_;
 398
 399         tree_.reset
 400         (
 401             new indexedOctree<treeDataTriSurface>
 402             (
 403                 treeDataTriSurface(*this),
 404                 bb,
 405                 10,     // maxLevel
 406                 10,     // leafsize
 407                 3.0     // duplicity
 408             )
 409         );
 410
 411         indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
 412     }
 413
 414     return tree_();
 415 }
 416
 417
 418 const Foam::indexedOctree<Foam::treeDataEdge>&
 419  Foam::triSurfaceMesh::edgeTree() const
 420 {
 421     if (edgeTree_.empty())
 422     {
 423         // Boundary edges
 424         labelList bEdges
 425         (
 426             identity
 427             (
 428                 nEdges()
 429                -nInternalEdges()
 430             )
 431           + nInternalEdges()
 432         );
 433
 434         // Random number generator. Bit dodgy since not exactly random ;-)
 435         Random rndGen(65431);
 436
 437         // Slightly extended bb. Slightly off-centred just so on symmetric
 438         // geometry there are less face/edge aligned items.
 439         treeBoundBox bb
 440         (
 441             treeBoundBox(points(), meshPoints()).extend(rndGen, 1E-4)
 442         );
 443         bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
 444         bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
 445
 446         edgeTree_.reset
 447         (
 448             new indexedOctree<treeDataEdge>
 449             (
 450                 treeDataEdge
 451                 (
 452                     false,          // cachebb
 453                     edges(),        // edges
 454                     localPoints(),  // points
 455                     bEdges          // selected edges
 456                 ),
 457                 bb,     // bb
 458                 8,      // maxLevel
 459                 10,     // leafsize
 460                 3.0     // duplicity
 461             )
 462         );
 463     }
 464     return edgeTree_();
 465 }
 466
 467
 468 const Foam::wordList& Foam::triSurfaceMesh::regions() const
 469 {
 470     if (regions_.empty())
 471     {
 472         regions_.setSize(patches().size());
 473         forAll(regions_, regionI)
 474         {
 475             regions_[regionI] = patches()[regionI].name();
 476         }
 477     }
 478     return regions_;
 479 }
 480
 481
 482 // Find out if surface is closed.
 483 bool Foam::triSurfaceMesh::hasVolumeType() const
 484 {
 485     if (surfaceClosed_ == -1)
 486     {
 487         if (isSurfaceClosed())
 488         {
 489             surfaceClosed_ = 1;
 490         }
 491         else
 492         {
 493             surfaceClosed_ = 0;
 494         }
 495     }
 496
 497     return surfaceClosed_ == 1;
 498 }
 499
 500
 501 void Foam::triSurfaceMesh::findNearest
 502 (
 503     const pointField& samples,
 504     const scalarField& nearestDistSqr,
 505     List<pointIndexHit>& info
 506 ) const
 507 {
 508     const indexedOctree<treeDataTriSurface>& octree = tree();
 509
 510     info.setSize(samples.size());
 511
 512     scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
 513     indexedOctree<treeDataTriSurface>::perturbTol() = tolerance_;
 514
 515     forAll(samples, i)
 516     {
 517         static_cast<pointIndexHit&>(info[i]) = octree.findNearest
 518         (
 519             samples[i],
 520             nearestDistSqr[i]
 521         );
 522     }
 523
 524     indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
 525 }
 526
 527
 528 void Foam::triSurfaceMesh::findLine
 529 (
 530     const pointField& start,
 531     const pointField& end,
 532     List<pointIndexHit>& info
 533 ) const
 534 {
 535     const indexedOctree<treeDataTriSurface>& octree = tree();
 536
 537     info.setSize(start.size());
 538
 539     scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
 540     indexedOctree<treeDataTriSurface>::perturbTol() = tolerance_;
 541
 542     forAll(start, i)
 543     {
 544         static_cast<pointIndexHit&>(info[i]) = octree.findLine
 545         (
 546             start[i],
 547             end[i]
 548         );
 549     }
 550
 551     indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
 552 }
 553
 554
 555 void Foam::triSurfaceMesh::findLineAny
 556 (
 557     const pointField& start,
 558     const pointField& end,
 559     List<pointIndexHit>& info
 560 ) const
 561 {
 562     const indexedOctree<treeDataTriSurface>& octree = tree();
 563
 564     info.setSize(start.size());
 565
 566     scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
 567     indexedOctree<treeDataTriSurface>::perturbTol() = tolerance_;
 568
 569     forAll(start, i)
 570     {
 571         static_cast<pointIndexHit&>(info[i]) = octree.findLineAny
 572         (
 573             start[i],
 574             end[i]
 575         );
 576     }
 577
 578     indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
 579 }
 580
 581
 582 void Foam::triSurfaceMesh::findLineAll
 583 (
 584     const pointField& start,
 585     const pointField& end,
 586     List<List<pointIndexHit> >& info
 587 ) const
 588 {
 589     const indexedOctree<treeDataTriSurface>& octree = tree();
 590
 591     info.setSize(start.size());
 592
 593     scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
 594     indexedOctree<treeDataTriSurface>::perturbTol() = tolerance_;
 595
 596     // Work array
 597     DynamicList<pointIndexHit, 1, 1> hits;
 598
 599     // Tolerances:
 600     // To find all intersections we add a small vector to the last intersection
 601     // This is chosen such that
 602     // - it is significant (SMALL is smallest representative relative tolerance;
 603     //   we need something bigger since we're doing calculations)
 604     // - if the start-end vector is zero we still progress
 605     const vectorField dirVec(end-start);
 606     const scalarField magSqrDirVec(magSqr(dirVec));
 607     const vectorField smallVec
 608     (
 609         indexedOctree<treeDataTriSurface>::perturbTol()*dirVec
 610       + vector(ROOTVSMALL,ROOTVSMALL,ROOTVSMALL)
 611     );
 612
 613     forAll(start, pointI)
 614     {
 615         // See if any intersection between pt and end
 616         pointIndexHit inter = octree.findLine(start[pointI], end[pointI]);
 617
 618         if (inter.hit())
 619         {
 620             hits.clear();
 621             hits.append(inter);
 622
 623             point pt = inter.hitPoint() + smallVec[pointI];
 624
 625             while (((pt-start[pointI])&dirVec[pointI]) <= magSqrDirVec[pointI])
 626             {
 627                 // See if any intersection between pt and end
 628                 pointIndexHit inter = octree.findLine(pt, end[pointI]);
 629
 630                 // Check for not hit or hit same triangle as before (can happen
 631                 // if vector along surface of triangle)
 632                 if
 633                 (
 634                     !inter.hit()
 635                  || (inter.index() == hits[hits.size()-1].index())
 636                 )
 637                 {
 638                     break;
 639                 }
 640                 hits.append(inter);
 641
 642                 pt = inter.hitPoint() + smallVec[pointI];
 643             }
 644
 645             info[pointI].transfer(hits);
 646         }
 647         else
 648         {
 649             info[pointI].clear();
 650         }
 651     }
 652
 653     indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
 654 }
 655
 656
 657 void Foam::triSurfaceMesh::getRegion
 658 (
 659     const List<pointIndexHit>& info,
 660     labelList& region
 661 ) const
 662 {
 663     region.setSize(info.size());
 664     forAll(info, i)
 665     {
 666         if (info[i].hit())
 667         {
 668             region[i] = triSurface::operator[](info[i].index()).region();
 669         }
 670         else
 671         {
 672             region[i] = -1;
 673         }
 674     }
 675 }
 676
 677
 678 void Foam::triSurfaceMesh::getNormal
 679 (
 680     const List<pointIndexHit>& info,
 681     vectorField& normal
 682 ) const
 683 {
 684     normal.setSize(info.size());
 685
 686     forAll(info, i)
 687     {
 688         if (info[i].hit())
 689         {
 690             label triI = info[i].index();
 691             //- Cached:
 692             //normal[i] = faceNormals()[triI];
 693
 694             //- Uncached
 695             normal[i] = triSurface::operator[](triI).normal(points());
 696             normal[i] /= mag(normal[i]) + VSMALL;
 697         }
 698         else
 699         {
 700             // Set to what?
 701             normal[i] = vector::zero;
 702         }
 703     }
 704 }
 705
 706
 707 void Foam::triSurfaceMesh::getField
 708 (
 709     const word& fieldName,
 710     const List<pointIndexHit>& info,
 711     labelList& values
 712 ) const
 713 {
 714     const triSurfaceLabelField& fld = lookupObject<triSurfaceLabelField>
 715     (
 716         fieldName
 717     );
 718
 719     values.setSize(info.size());
 720     forAll(info, i)
 721     {
 722         if (info[i].hit())
 723         {
 724             values[i] = fld[info[i].index()];
 725         }
 726     }
 727 }
 728
 729
 730 void Foam::triSurfaceMesh::getVolumeType
 731 (
 732     const pointField& points,
 733     List<volumeType>& volType
 734 ) const
 735 {
 736     volType.setSize(points.size());
 737
 738     scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
 739     indexedOctree<treeDataTriSurface>::perturbTol() = tolerance_;
 740
 741     forAll(points, pointI)
 742     {
 743         const point& pt = points[pointI];
 744
 745         // - use cached volume type per each tree node
 746         // - cheat conversion since same values
 747         volType[pointI] = static_cast<volumeType>(tree().getVolumeType(pt));
 748     }
 749
 750     indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
 751 }
 752
 753
 754 //- Write using given format, version and compression
 755 bool Foam::triSurfaceMesh::writeObject
 756 (
 757     IOstream::streamFormat fmt,
 758     IOstream::versionNumber ver,
 759     IOstream::compressionType cmp
 760 ) const
 761 {
 762     fileName fullPath(searchableSurface::objectPath());
 763
 764     if (!mkDir(fullPath.path()))
 765     {
 766         return false;
 767     }
 768
 769     triSurface::write(fullPath);
 770
 771     if (!isFile(fullPath))
 772     {
 773         return false;
 774     }
 775
 776     //return objectRegistry::writeObject(fmt, ver, cmp);
 777     return true;
 778 }
 779
 780
 781 // ************************************************************************* //