src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.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 "surfaceFeatures.H"
  27 #include "triSurface.H"
  28 #include "indexedOctree.H"
  29 #include "treeDataEdge.H"
  30 #include "treeDataPoint.H"
  31 #include "meshTools.H"
  32 #include "linePointRef.H"
  33 #include "OFstream.H"
  34 #include "IFstream.H"
  35 #include "unitConversion.H"
  36
  37 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  38
  39 defineTypeNameAndDebug(Foam::surfaceFeatures, 0);
  40
  41
  42 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  43
  44 //- Get nearest point on edge and classify position on edge.
  45 Foam::pointIndexHit Foam::surfaceFeatures::edgeNearest
  46 (
  47     const point& start,
  48     const point& end,
  49     const point& sample
  50 )
  51 {
  52     pointHit eHit = linePointRef(start, end).nearestDist(sample);
  53
  54     // Classification of position on edge.
  55     label endPoint;
  56
  57     if (eHit.hit())
  58     {
  59         // Nearest point is on edge itself.
  60         // Note: might be at or very close to endpoint. We should use tolerance
  61         // here.
  62         endPoint = -1;
  63     }
  64     else
  65     {
  66         // Nearest point has to be one of the end points. Find out
  67         // which one.
  68         if
  69         (
  70             mag(eHit.rawPoint() - start)
  71           < mag(eHit.rawPoint() - end)
  72         )
  73         {
  74             endPoint = 0;
  75         }
  76         else
  77         {
  78             endPoint = 1;
  79         }
  80     }
  81
  82     return pointIndexHit(eHit.hit(), eHit.rawPoint(), endPoint);
  83 }
  84
  85
  86 // Go from selected edges only to a value for every edge
  87 Foam::List<Foam::surfaceFeatures::edgeStatus> Foam::surfaceFeatures::toStatus()
  88  const
  89 {
  90     List<edgeStatus> edgeStat(surf_.nEdges(), NONE);
  91
  92     // Region edges first
  93     for (label i = 0; i < externalStart_; i++)
  94     {
  95         edgeStat[featureEdges_[i]] = REGION;
  96     }
  97
  98     // External edges
  99     for (label i = externalStart_; i < internalStart_; i++)
 100     {
 101         edgeStat[featureEdges_[i]] = EXTERNAL;
 102     }
 103
 104     // Internal edges
 105     for (label i = internalStart_; i < featureEdges_.size(); i++)
 106     {
 107         edgeStat[featureEdges_[i]] = INTERNAL;
 108     }
 109
 110     return edgeStat;
 111 }
 112
 113
 114 // Set from value for every edge
 115 void Foam::surfaceFeatures::setFromStatus(const List<edgeStatus>& edgeStat)
 116 {
 117     // Count
 118
 119     label nRegion = 0;
 120     label nExternal = 0;
 121     label nInternal = 0;
 122
 123     forAll(edgeStat, edgeI)
 124     {
 125         if (edgeStat[edgeI] == REGION)
 126         {
 127             nRegion++;
 128         }
 129         else if (edgeStat[edgeI] == EXTERNAL)
 130         {
 131             nExternal++;
 132         }
 133         else if (edgeStat[edgeI] == INTERNAL)
 134         {
 135             nInternal++;
 136         }
 137     }
 138
 139     externalStart_ = nRegion;
 140     internalStart_ = externalStart_ + nExternal;
 141
 142
 143     // Copy
 144
 145     featureEdges_.setSize(internalStart_ + nInternal);
 146
 147     label regionI = 0;
 148     label externalI = externalStart_;
 149     label internalI = internalStart_;
 150
 151     forAll(edgeStat, edgeI)
 152     {
 153         if (edgeStat[edgeI] == REGION)
 154         {
 155             featureEdges_[regionI++] = edgeI;
 156         }
 157         else if (edgeStat[edgeI] == EXTERNAL)
 158         {
 159             featureEdges_[externalI++] = edgeI;
 160         }
 161         else if (edgeStat[edgeI] == INTERNAL)
 162         {
 163             featureEdges_[internalI++] = edgeI;
 164         }
 165     }
 166
 167     // Calculate consistent feature points
 168     calcFeatPoints(edgeStat);
 169 }
 170
 171
 172 //construct feature points where more than 2 feature edges meet
 173 void Foam::surfaceFeatures::calcFeatPoints(const List<edgeStatus>& edgeStat)
 174 {
 175     DynamicList<label> featurePoints(surf_.nPoints()/1000);
 176
 177     const labelListList& pointEdges = surf_.pointEdges();
 178
 179     forAll(pointEdges, pointI)
 180     {
 181         const labelList& pEdges = pointEdges[pointI];
 182
 183         label nFeatEdges = 0;
 184
 185         forAll(pEdges, i)
 186         {
 187             if (edgeStat[pEdges[i]] != NONE)
 188             {
 189                 nFeatEdges++;
 190             }
 191         }
 192
 193         if (nFeatEdges > 2)
 194         {
 195             featurePoints.append(pointI);
 196         }
 197     }
 198
 199     featurePoints_.transfer(featurePoints);
 200 }
 201
 202
 203 // Returns next feature edge connected to pointI with correct value.
 204 Foam::label Foam::surfaceFeatures::nextFeatEdge
 205 (
 206     const List<edgeStatus>& edgeStat,
 207     const labelList& featVisited,
 208     const label unsetVal,
 209     const label prevEdgeI,
 210     const label vertI
 211 ) const
 212 {
 213     const labelList& pEdges = surf_.pointEdges()[vertI];
 214
 215     label nextEdgeI = -1;
 216
 217     forAll(pEdges, i)
 218     {
 219         label edgeI = pEdges[i];
 220
 221         if
 222         (
 223             edgeI != prevEdgeI
 224          && edgeStat[edgeI] != NONE
 225          && featVisited[edgeI] == unsetVal
 226         )
 227         {
 228             if (nextEdgeI == -1)
 229             {
 230                 nextEdgeI = edgeI;
 231             }
 232             else
 233             {
 234                 // More than one feature edge to choose from. End of segment.
 235                 return -1;
 236             }
 237         }
 238     }
 239
 240     return nextEdgeI;
 241 }
 242
 243
 244 // Finds connected feature edges by walking from prevEdgeI in direction of
 245 // prevPointI. Marks feature edges visited in featVisited by assigning them
 246 // the current feature line number. Returns cumulative length of edges walked.
 247 // Works in one of two modes:
 248 // - mark : step to edges with featVisited = -1.
 249 //          Mark edges visited with currentFeatI.
 250 // - clear : step to edges with featVisited = currentFeatI
 251 //           Mark edges visited with -2 and erase from feature edges.
 252 Foam::surfaceFeatures::labelScalar Foam::surfaceFeatures::walkSegment
 253 (
 254     const bool mark,
 255     const List<edgeStatus>& edgeStat,
 256     const label startEdgeI,
 257     const label startPointI,
 258     const label currentFeatI,
 259     labelList& featVisited
 260 )
 261 {
 262     label edgeI = startEdgeI;
 263
 264     label vertI = startPointI;
 265
 266     scalar visitedLength = 0.0;
 267
 268     label nVisited = 0;
 269
 270     if (findIndex(featurePoints_, startPointI) >= 0)
 271     {
 272         // Do not walk across feature points
 273
 274         return labelScalar(nVisited, visitedLength);
 275     }
 276
 277     //
 278     // Now we have:
 279     //    edgeI : first edge on this segment
 280     //    vertI : one of the endpoints of this segment
 281     //
 282     // Start walking, marking off edges (in featVisited)
 283     // as we go along.
 284     //
 285
 286     label unsetVal;
 287
 288     if (mark)
 289     {
 290         unsetVal = -1;
 291     }
 292     else
 293     {
 294         unsetVal = currentFeatI;
 295     }
 296
 297     do
 298     {
 299         // Step to next feature edge with value unsetVal
 300         edgeI = nextFeatEdge(edgeStat, featVisited, unsetVal, edgeI, vertI);
 301
 302         if (edgeI == -1 || edgeI == startEdgeI)
 303         {
 304             break;
 305         }
 306
 307         // Mark with current value. If in clean mode also remove feature edge.
 308
 309         if (mark)
 310         {
 311             featVisited[edgeI] = currentFeatI;
 312         }
 313         else
 314         {
 315             featVisited[edgeI] = -2;
 316         }
 317
 318
 319         // Step to next vertex on edge
 320         const edge& e = surf_.edges()[edgeI];
 321
 322         vertI = e.otherVertex(vertI);
 323
 324
 325         // Update cumulative length
 326         visitedLength += e.mag(surf_.localPoints());
 327
 328         nVisited++;
 329
 330         if (nVisited > surf_.nEdges())
 331         {
 332             Warning<< "walkSegment : reached iteration limit in walking "
 333                 << "feature edges on surface from edge:" << startEdgeI
 334                 << " vertex:" << startPointI << nl
 335                 << "Returning with large length" << endl;
 336
 337             return labelScalar(nVisited, GREAT);
 338         }
 339     }
 340     while (true);
 341
 342     return labelScalar(nVisited, visitedLength);
 343 }
 344
 345
 346 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 347
 348 Foam::surfaceFeatures::surfaceFeatures(const triSurface& surf)
 349 :
 350     surf_(surf),
 351     featurePoints_(0),
 352     featureEdges_(0),
 353     externalStart_(0),
 354     internalStart_(0)
 355 {}
 356
 357
 358 // Construct from components.
 359 Foam::surfaceFeatures::surfaceFeatures
 360 (
 361     const triSurface& surf,
 362     const labelList& featurePoints,
 363     const labelList& featureEdges,
 364     const label externalStart,
 365     const label internalStart
 366 )
 367 :
 368     surf_(surf),
 369     featurePoints_(featurePoints),
 370     featureEdges_(featureEdges),
 371     externalStart_(externalStart),
 372     internalStart_(externalStart)
 373 {}
 374
 375
 376 // Construct from surface, angle and min length
 377 Foam::surfaceFeatures::surfaceFeatures
 378 (
 379     const triSurface& surf,
 380     const scalar includedAngle,
 381     const scalar minLen,
 382     const label minElems
 383 )
 384 :
 385     surf_(surf),
 386     featurePoints_(0),
 387     featureEdges_(0),
 388     externalStart_(0),
 389     internalStart_(0)
 390 {
 391     findFeatures(includedAngle);
 392
 393     if (minLen > 0 || minElems > 0)
 394     {
 395         trimFeatures(minLen, minElems);
 396     }
 397 }
 398
 399
 400 //- Construct from dictionary
 401 Foam::surfaceFeatures::surfaceFeatures
 402 (
 403     const triSurface& surf,
 404     const dictionary& featInfoDict
 405 )
 406 :
 407     surf_(surf),
 408     featurePoints_(featInfoDict.lookup("featurePoints")),
 409     featureEdges_(featInfoDict.lookup("featureEdges")),
 410     externalStart_(readLabel(featInfoDict.lookup("externalStart"))),
 411     internalStart_(readLabel(featInfoDict.lookup("internalStart")))
 412 {}
 413
 414
 415 //- Construct from file
 416 Foam::surfaceFeatures::surfaceFeatures
 417 (
 418     const triSurface& surf,
 419     const fileName& fName
 420 )
 421 :
 422     surf_(surf),
 423     featurePoints_(0),
 424     featureEdges_(0),
 425     externalStart_(0),
 426     internalStart_(0)
 427 {
 428     IFstream str(fName);
 429
 430     dictionary featInfoDict(str);
 431
 432     featureEdges_ = labelList(featInfoDict.lookup("featureEdges"));
 433     featurePoints_ = labelList(featInfoDict.lookup("featurePoints"));
 434     externalStart_ = readLabel(featInfoDict.lookup("externalStart"));
 435     internalStart_ = readLabel(featInfoDict.lookup("internalStart"));
 436 }
 437
 438
 439 //- Construct as copy
 440 Foam::surfaceFeatures::surfaceFeatures(const surfaceFeatures& sf)
 441 :
 442     surf_(sf.surface()),
 443     featurePoints_(sf.featurePoints()),
 444     featureEdges_(sf.featureEdges()),
 445     externalStart_(sf.externalStart()),
 446     internalStart_(sf.internalStart())
 447 {}
 448
 449
 450 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 451
 452 Foam::labelList Foam::surfaceFeatures::selectFeatureEdges
 453 (
 454     const bool regionEdges,
 455     const bool externalEdges,
 456     const bool internalEdges
 457 ) const
 458 {
 459     DynamicList<label> selectedEdges;
 460
 461     if (regionEdges)
 462     {
 463         selectedEdges.setCapacity(selectedEdges.size() + nRegionEdges());
 464
 465         for (label i = 0; i < externalStart_; i++)
 466         {
 467             selectedEdges.append(featureEdges_[i]);
 468         }
 469     }
 470
 471     if (externalEdges)
 472     {
 473         selectedEdges.setCapacity(selectedEdges.size() + nExternalEdges());
 474
 475         for (label i = externalStart_; i < internalStart_; i++)
 476         {
 477             selectedEdges.append(featureEdges_[i]);
 478         }
 479     }
 480
 481     if (internalEdges)
 482     {
 483         selectedEdges.setCapacity(selectedEdges.size() + nInternalEdges());
 484
 485         for (label i = internalStart_; i < featureEdges_.size(); i++)
 486         {
 487             selectedEdges.append(featureEdges_[i]);
 488         }
 489     }
 490
 491     return selectedEdges.shrink();
 492 }
 493
 494
 495 void Foam::surfaceFeatures::findFeatures(const scalar includedAngle)
 496 {
 497     scalar minCos = Foam::cos(degToRad(180.0 - includedAngle));
 498
 499     const labelListList& edgeFaces = surf_.edgeFaces();
 500     const vectorField& faceNormals = surf_.faceNormals();
 501     const pointField& points = surf_.points();
 502
 503     // Per edge whether is feature edge.
 504     List<edgeStatus> edgeStat(surf_.nEdges(), NONE);
 505
 506     forAll(edgeFaces, edgeI)
 507     {
 508         const labelList& eFaces = edgeFaces[edgeI];
 509
 510         if (eFaces.size() != 2)
 511         {
 512             // Non-manifold. What to do here? Is region edge? external edge?
 513             edgeStat[edgeI] = REGION;
 514         }
 515         else
 516         {
 517             label face0 = eFaces[0];
 518             label face1 = eFaces[1];
 519
 520             if (surf_[face0].region() != surf_[face1].region())
 521             {
 522                 edgeStat[edgeI] = REGION;
 523             }
 524             else
 525             {
 526                 if ((faceNormals[face0] & faceNormals[face1]) < minCos)
 527                 {
 528
 529                     // Check if convex or concave by looking at angle
 530                     // between face centres and normal
 531                     vector f0Tof1 =
 532                         surf_[face1].centre(points)
 533                         - surf_[face0].centre(points);
 534
 535                     if ((f0Tof1 & faceNormals[face0]) > 0.0)
 536                     {
 537                         edgeStat[edgeI] = INTERNAL;
 538                     }
 539                     else
 540                     {
 541                         edgeStat[edgeI] = EXTERNAL;
 542                     }
 543                 }
 544             }
 545         }
 546     }
 547
 548     setFromStatus(edgeStat);
 549 }
 550
 551
 552 // Remove small strings of edges. First assigns string number to
 553 // every edge and then works out the length of them.
 554 void Foam::surfaceFeatures::trimFeatures
 555 (
 556     const scalar minLen,
 557     const label minElems
 558 )
 559 {
 560     // Convert feature edge list to status per edge.
 561     List<edgeStatus> edgeStat(toStatus());
 562
 563
 564     // Mark feature edges according to connected lines.
 565     // -1: unassigned
 566     // -2: part of too small a feature line
 567     // >0: feature line number
 568     labelList featLines(surf_.nEdges(), -1);
 569
 570     // Current featureline number.
 571     label featI = 0;
 572
 573     // Current starting edge
 574     label startEdgeI = 0;
 575
 576     do
 577     {
 578         // Find unset featureline
 579         for (; startEdgeI < edgeStat.size(); startEdgeI++)
 580         {
 581             if
 582             (
 583                 edgeStat[startEdgeI] != NONE
 584              && featLines[startEdgeI] == -1
 585             )
 586             {
 587                 // Found unassigned feature edge.
 588                 break;
 589             }
 590         }
 591
 592         if (startEdgeI == edgeStat.size())
 593         {
 594             // No unset feature edge found. All feature edges have line number
 595             // assigned.
 596             break;
 597         }
 598
 599         featLines[startEdgeI] = featI;
 600
 601         const edge& startEdge = surf_.edges()[startEdgeI];
 602
 603         // Walk 'left' and 'right' from startEdge.
 604         labelScalar leftPath =
 605             walkSegment
 606             (
 607                 true,           // 'mark' mode
 608                 edgeStat,
 609                 startEdgeI,     // edge, not yet assigned to a featureLine
 610                 startEdge[0],   // start of edge
 611                 featI,          // Mark value
 612                 featLines       // Mark for all edges
 613             );
 614
 615         labelScalar rightPath =
 616             walkSegment
 617             (
 618                 true,
 619                 edgeStat,
 620                 startEdgeI,
 621                 startEdge[1],   // end of edge
 622                 featI,
 623                 featLines
 624             );
 625
 626         if
 627         (
 628             (
 629                 leftPath.len_
 630               + rightPath.len_
 631               + startEdge.mag(surf_.localPoints())
 632               < minLen
 633             )
 634          || (leftPath.n_ + rightPath.n_ + 1 < minElems)
 635         )
 636         {
 637             // Rewalk same route (recognizable by featLines == featI)
 638             // to reset featLines.
 639
 640             featLines[startEdgeI] = -2;
 641
 642             walkSegment
 643             (
 644                 false,          // 'clean' mode
 645                 edgeStat,
 646                 startEdgeI,     // edge, not yet assigned to a featureLine
 647                 startEdge[0],   // start of edge
 648                 featI,          // Unset value
 649                 featLines       // Mark for all edges
 650             );
 651
 652             walkSegment
 653             (
 654                 false,
 655                 edgeStat,
 656                 startEdgeI,
 657                 startEdge[1],   // end of edge
 658                 featI,
 659                 featLines
 660             );
 661         }
 662         else
 663         {
 664             featI++;
 665         }
 666     }
 667     while (true);
 668
 669     // Unmark all feature lines that have featLines=-2
 670     forAll(featureEdges_, i)
 671     {
 672         label edgeI = featureEdges_[i];
 673
 674         if (featLines[edgeI] == -2)
 675         {
 676             edgeStat[edgeI] = NONE;
 677         }
 678     }
 679
 680     // Convert back to edge labels
 681     setFromStatus(edgeStat);
 682 }
 683
 684
 685 void Foam::surfaceFeatures::writeDict(Ostream& writeFile) const
 686 {
 687
 688     dictionary featInfoDict;
 689     featInfoDict.add("externalStart", externalStart_);
 690     featInfoDict.add("internalStart", internalStart_);
 691     featInfoDict.add("featureEdges", featureEdges_);
 692     featInfoDict.add("featurePoints", featurePoints_);
 693
 694     featInfoDict.write(writeFile);
 695 }
 696
 697
 698 void Foam::surfaceFeatures::write(const fileName& fName) const
 699 {
 700     OFstream str(fName);
 701
 702     writeDict(str);
 703 }
 704
 705
 706 void Foam::surfaceFeatures::writeObj(const fileName& prefix) const
 707 {
 708     OFstream regionStr(prefix + "_regionEdges.obj");
 709     Pout<< "Writing region edges to " << regionStr.name() << endl;
 710
 711     label verti = 0;
 712     for (label i = 0; i < externalStart_; i++)
 713     {
 714         const edge& e = surf_.edges()[featureEdges_[i]];
 715
 716         meshTools::writeOBJ(regionStr, surf_.localPoints()[e[0]]); verti++;
 717         meshTools::writeOBJ(regionStr, surf_.localPoints()[e[1]]); verti++;
 718         regionStr << "l " << verti-1 << ' ' << verti << endl;
 719     }
 720
 721
 722     OFstream externalStr(prefix + "_externalEdges.obj");
 723     Pout<< "Writing external edges to " << externalStr.name() << endl;
 724
 725     verti = 0;
 726     for (label i = externalStart_; i < internalStart_; i++)
 727     {
 728         const edge& e = surf_.edges()[featureEdges_[i]];
 729
 730         meshTools::writeOBJ(externalStr, surf_.localPoints()[e[0]]); verti++;
 731         meshTools::writeOBJ(externalStr, surf_.localPoints()[e[1]]); verti++;
 732         externalStr << "l " << verti-1 << ' ' << verti << endl;
 733     }
 734
 735     OFstream internalStr(prefix + "_internalEdges.obj");
 736     Pout<< "Writing internal edges to " << internalStr.name() << endl;
 737
 738     verti = 0;
 739     for (label i = internalStart_; i < featureEdges_.size(); i++)
 740     {
 741         const edge& e = surf_.edges()[featureEdges_[i]];
 742
 743         meshTools::writeOBJ(internalStr, surf_.localPoints()[e[0]]); verti++;
 744         meshTools::writeOBJ(internalStr, surf_.localPoints()[e[1]]); verti++;
 745         internalStr << "l " << verti-1 << ' ' << verti << endl;
 746     }
 747
 748     OFstream pointStr(prefix + "_points.obj");
 749     Pout<< "Writing feature points to " << pointStr.name() << endl;
 750
 751     forAll(featurePoints_, i)
 752     {
 753         label pointI = featurePoints_[i];
 754
 755         meshTools::writeOBJ(pointStr, surf_.localPoints()[pointI]);
 756     }
 757 }
 758
 759
 760 // Get nearest vertex on patch for every point of surf in pointSet.
 761 Foam::Map<Foam::label> Foam::surfaceFeatures::nearestSamples
 762 (
 763     const labelList& pointLabels,
 764     const pointField& samples,
 765     const scalarField& maxDistSqr
 766 ) const
 767 {
 768     // Build tree out of all samples.
 769
 770     // Note: cannot be done one the fly - gcc4.4 compiler bug.
 771     treeBoundBox bb(samples);
 772
 773     indexedOctree<treeDataPoint> ppTree
 774     (
 775         treeDataPoint(samples),   // all information needed to do checks
 776         bb,                       // overall search domain
 777         8,      // maxLevel
 778         10,     // leafsize
 779         3.0     // duplicity
 780     );
 781
 782     // From patch point to surface point
 783     Map<label> nearest(2*pointLabels.size());
 784
 785     const pointField& surfPoints = surf_.localPoints();
 786
 787     forAll(pointLabels, i)
 788     {
 789         label surfPointI = pointLabels[i];
 790
 791         const point& surfPt = surfPoints[surfPointI];
 792
 793         pointIndexHit info = ppTree.findNearest
 794         (
 795             surfPt,
 796             maxDistSqr[i]
 797         );
 798
 799         if (!info.hit())
 800         {
 801             FatalErrorIn("surfaceFeatures::nearestSamples")
 802                 << "Problem for point "
 803                 << surfPointI << " in tree " << ppTree.bb()
 804                 << abort(FatalError);
 805         }
 806
 807         label sampleI = info.index();
 808
 809         if (magSqr(samples[sampleI] - surfPt) < maxDistSqr[sampleI])
 810         {
 811             nearest.insert(sampleI, surfPointI);
 812         }
 813     }
 814
 815
 816     if (debug)
 817     {
 818         //
 819         // Dump to obj file
 820         //
 821
 822         Pout<< endl
 823             << "Dumping nearest surface feature points to nearestSamples.obj"
 824             << endl
 825             << "View this Lightwave-OBJ file with e.g. javaview" << endl
 826             << endl;
 827
 828         OFstream objStream("nearestSamples.obj");
 829
 830         label vertI = 0;
 831         forAllConstIter(Map<label>, nearest, iter)
 832         {
 833             meshTools::writeOBJ(objStream, samples[iter.key()]); vertI++;
 834             meshTools::writeOBJ(objStream, surfPoints[iter()]); vertI++;
 835             objStream<< "l " << vertI-1 << ' ' << vertI << endl;
 836         }
 837     }
 838
 839     return nearest;
 840 }
 841
 842
 843 // Get nearest sample point for regularly sampled points along
 844 // selected edges. Return map from sample to edge label
 845 Foam::Map<Foam::label> Foam::surfaceFeatures::nearestSamples
 846 (
 847     const labelList& selectedEdges,
 848     const pointField& samples,
 849     const scalarField& sampleDist,
 850     const scalarField& maxDistSqr,
 851     const scalar minSampleDist
 852 ) const
 853 {
 854     const pointField& surfPoints = surf_.localPoints();
 855     const edgeList& surfEdges = surf_.edges();
 856
 857     scalar maxSearchSqr = max(maxDistSqr);
 858
 859     //Note: cannot be done one the fly - gcc4.4 compiler bug.
 860     treeBoundBox bb(samples);
 861
 862     indexedOctree<treeDataPoint> ppTree
 863     (
 864         treeDataPoint(samples),   // all information needed to do checks
 865         bb,                         // overall search domain
 866         8,      // maxLevel
 867         10,     // leafsize
 868         3.0     // duplicity
 869     );
 870
 871     // From patch point to surface edge.
 872     Map<label> nearest(2*selectedEdges.size());
 873
 874     forAll(selectedEdges, i)
 875     {
 876         label surfEdgeI = selectedEdges[i];
 877
 878         const edge& e = surfEdges[surfEdgeI];
 879
 880         if (debug && (i % 1000) == 0)
 881         {
 882             Pout<< "looking at surface feature edge " << surfEdgeI
 883                 << " verts:" << e << " points:" << surfPoints[e[0]]
 884                 << ' ' << surfPoints[e[1]] << endl;
 885         }
 886
 887         // Normalized edge vector
 888         vector eVec = e.vec(surfPoints);
 889         scalar eMag = mag(eVec);
 890         eVec /= eMag;
 891
 892
 893         //
 894         // Sample along edge
 895         //
 896
 897         bool exit = false;
 898
 899         // Coordinate along edge (0 .. eMag)
 900         scalar s = 0.0;
 901
 902         while (true)
 903         {
 904             point edgePoint(surfPoints[e.start()] + s*eVec);
 905
 906             pointIndexHit info = ppTree.findNearest
 907             (
 908                 edgePoint,
 909                 maxSearchSqr
 910             );
 911
 912             if (!info.hit())
 913             {
 914                 // No point close enough to surface edge.
 915                 break;
 916             }
 917
 918             label sampleI = info.index();
 919
 920             if (magSqr(info.hitPoint() - edgePoint) < maxDistSqr[sampleI])
 921             {
 922                 nearest.insert(sampleI, surfEdgeI);
 923             }
 924
 925             if (exit)
 926             {
 927                 break;
 928             }
 929
 930             // Step to next sample point using local distance.
 931             // Truncate to max 1/minSampleDist samples per feature edge.
 932             s += max(minSampleDist*eMag, sampleDist[sampleI]);
 933
 934             if (s >= (1-minSampleDist)*eMag)
 935             {
 936                 // Do one more sample, at edge endpoint
 937                 s = eMag;
 938                 exit = true;
 939             }
 940         }
 941     }
 942
 943
 944
 945     if (debug)
 946     {
 947         // Dump to obj file
 948
 949         Pout<< "Dumping nearest surface edges to nearestEdges.obj\n"
 950             << "View this Lightwave-OBJ file with e.g. javaview\n" << endl;
 951
 952         OFstream objStream("nearestEdges.obj");
 953
 954         label vertI = 0;
 955         forAllConstIter(Map<label>, nearest, iter)
 956         {
 957             const label sampleI = iter.key();
 958
 959             meshTools::writeOBJ(objStream, samples[sampleI]); vertI++;
 960
 961             const edge& e = surfEdges[iter()];
 962
 963             point nearPt =
 964                 e.line(surfPoints).nearestDist(samples[sampleI]).rawPoint();
 965
 966             meshTools::writeOBJ(objStream, nearPt); vertI++;
 967
 968             objStream<< "l " << vertI-1 << ' ' << vertI << endl;
 969         }
 970     }
 971
 972     return nearest;
 973 }
 974
 975
 976 // Get nearest edge on patch for regularly sampled points along the feature
 977 // edges. Return map from patch edge to feature edge.
 978 //
 979 // Q: using point-based sampleDist and maxDist (distance to look around
 980 // each point). Should they be edge-based e.g. by averaging or max()?
 981 Foam::Map<Foam::pointIndexHit> Foam::surfaceFeatures::nearestEdges
 982 (
 983     const labelList& selectedEdges,
 984     const edgeList& sampleEdges,
 985     const labelList& selectedSampleEdges,
 986     const pointField& samplePoints,
 987     const scalarField& sampleDist,
 988     const scalarField& maxDistSqr,
 989     const scalar minSampleDist
 990 ) const
 991 {
 992     // Build tree out of selected sample edges.
 993     indexedOctree<treeDataEdge> ppTree
 994     (
 995         treeDataEdge
 996         (
 997             false,
 998             sampleEdges,
 999             samplePoints,
1000             selectedSampleEdges
1001         ),                          // geometric info container for edges
1002         treeBoundBox(samplePoints), // overall search domain
1003         8,      // maxLevel
1004         10,     // leafsize
1005         3.0     // duplicity
1006     );
1007
1008     const pointField& surfPoints = surf_.localPoints();
1009     const edgeList& surfEdges = surf_.edges();
1010
1011     scalar maxSearchSqr = max(maxDistSqr);
1012
1013     Map<pointIndexHit> nearest(2*sampleEdges.size());
1014
1015     //
1016     // Loop over all selected edges. Sample at regular intervals. Find nearest
1017     // sampleEdges (using octree)
1018     //
1019
1020     forAll(selectedEdges, i)
1021     {
1022         label surfEdgeI = selectedEdges[i];
1023
1024         const edge& e = surfEdges[surfEdgeI];
1025
1026         if (debug && (i % 1000) == 0)
1027         {
1028             Pout<< "looking at surface feature edge " << surfEdgeI
1029                 << " verts:" << e << " points:" << surfPoints[e[0]]
1030                 << ' ' << surfPoints[e[1]] << endl;
1031         }
1032
1033         // Normalized edge vector
1034         vector eVec = e.vec(surfPoints);
1035         scalar eMag = mag(eVec);
1036         eVec /= eMag;
1037
1038
1039         //
1040         // Sample along edge
1041         //
1042
1043         bool exit = false;
1044
1045         // Coordinate along edge (0 .. eMag)
1046         scalar s = 0.0;
1047
1048         while (true)
1049         {
1050             point edgePoint(surfPoints[e.start()] + s*eVec);
1051
1052             pointIndexHit info = ppTree.findNearest
1053             (
1054                 edgePoint,
1055                 maxSearchSqr
1056             );
1057
1058             if (!info.hit())
1059             {
1060                 // No edge close enough to surface edge.
1061                 break;
1062             }
1063
1064             label index = info.index();
1065
1066             label sampleEdgeI = ppTree.shapes().edgeLabels()[index];
1067
1068             const edge& e = sampleEdges[sampleEdgeI];
1069
1070             if (magSqr(info.hitPoint() - edgePoint) < maxDistSqr[e.start()])
1071             {
1072                 nearest.insert
1073                 (
1074                     sampleEdgeI,
1075                     pointIndexHit(true, edgePoint, surfEdgeI)
1076                 );
1077             }
1078
1079             if (exit)
1080             {
1081                 break;
1082             }
1083
1084             // Step to next sample point using local distance.
1085             // Truncate to max 1/minSampleDist samples per feature edge.
1086             // s += max(minSampleDist*eMag, sampleDist[e.start()]);
1087             s += 0.01*eMag;
1088
1089             if (s >= (1-minSampleDist)*eMag)
1090             {
1091                 // Do one more sample, at edge endpoint
1092                 s = eMag;
1093                 exit = true;
1094             }
1095         }
1096     }
1097
1098
1099     if (debug)
1100     {
1101         // Dump to obj file
1102
1103         Pout<< "Dumping nearest surface feature edges to nearestEdges.obj\n"
1104             << "View this Lightwave-OBJ file with e.g. javaview\n" << endl;
1105
1106         OFstream objStream("nearestEdges.obj");
1107
1108         label vertI = 0;
1109         forAllConstIter(Map<pointIndexHit>, nearest, iter)
1110         {
1111             const label sampleEdgeI = iter.key();
1112
1113             const edge& sampleEdge = sampleEdges[sampleEdgeI];
1114
1115             // Write line from edgeMid to point on feature edge
1116             meshTools::writeOBJ(objStream, sampleEdge.centre(samplePoints));
1117             vertI++;
1118
1119             meshTools::writeOBJ(objStream, iter().rawPoint());
1120             vertI++;
1121
1122             objStream<< "l " << vertI-1 << ' ' << vertI << endl;
1123         }
1124     }
1125
1126     return nearest;
1127 }
1128
1129
1130 // Get nearest surface edge for every sample. Return in form of
1131 // labelLists giving surfaceEdge label&intersectionpoint.
1132 void Foam::surfaceFeatures::nearestSurfEdge
1133 (
1134     const labelList& selectedEdges,
1135     const pointField& samples,
1136     scalar searchSpanSqr,   // Search span
1137     labelList& edgeLabel,
1138     labelList& edgeEndPoint,
1139     pointField& edgePoint
1140 ) const
1141 {
1142     edgeLabel.setSize(samples.size());
1143     edgeEndPoint.setSize(samples.size());
1144     edgePoint.setSize(samples.size());
1145
1146     const pointField& localPoints = surf_.localPoints();
1147
1148     indexedOctree<treeDataEdge> ppTree
1149     (
1150         treeDataEdge
1151         (
1152             false,
1153             surf_.edges(),
1154             localPoints,
1155             selectedEdges
1156         ),          // all information needed to do geometric checks
1157         treeBoundBox(localPoints),  // overall search domain
1158         8,      // maxLevel
1159         10,     // leafsize
1160         3.0     // duplicity
1161     );
1162
1163     forAll(samples, i)
1164     {
1165         const point& sample = samples[i];
1166
1167         pointIndexHit info = ppTree.findNearest
1168         (
1169             sample,
1170             searchSpanSqr
1171         );
1172
1173         if (!info.hit())
1174         {
1175             edgeLabel[i] = -1;
1176         }
1177         else
1178         {
1179             edgeLabel[i] = selectedEdges[info.index()];
1180
1181             // Need to recalculate to classify edgeEndPoint
1182             const edge& e = surf_.edges()[edgeLabel[i]];
1183
1184             pointIndexHit pHit = edgeNearest
1185             (
1186                 localPoints[e.start()],
1187                 localPoints[e.end()],
1188                 sample
1189             );
1190
1191             edgePoint[i] = pHit.rawPoint();
1192             edgeEndPoint[i] = pHit.index();
1193         }
1194     }
1195 }
1196
1197
1198 // Get nearest point on nearest feature edge for every sample (is edge)
1199 void Foam::surfaceFeatures::nearestSurfEdge
1200 (
1201     const labelList& selectedEdges,
1202
1203     const edgeList& sampleEdges,
1204     const labelList& selectedSampleEdges,
1205     const pointField& samplePoints,
1206
1207     const vector& searchSpan,   // Search span
1208     labelList& edgeLabel,       // label of surface edge or -1
1209     pointField& pointOnEdge,    // point on above edge
1210     pointField& pointOnFeature  // point on sample edge
1211 ) const
1212 {
1213     edgeLabel.setSize(selectedSampleEdges.size());
1214     pointOnEdge.setSize(selectedSampleEdges.size());
1215     pointOnFeature.setSize(selectedSampleEdges.size());
1216
1217     indexedOctree<treeDataEdge> ppTree
1218     (
1219         treeDataEdge
1220         (
1221             false,
1222             surf_.edges(),
1223             surf_.localPoints(),
1224             selectedEdges
1225         ),          // all information needed to do geometric checks
1226         treeBoundBox(surf_.localPoints()),  // overall search domain
1227         8,      // maxLevel
1228         10,     // leafsize
1229         3.0     // duplicity
1230     );
1231
1232     forAll(selectedSampleEdges, i)
1233     {
1234         const edge& e = sampleEdges[selectedSampleEdges[i]];
1235
1236         linePointRef edgeLine = e.line(samplePoints);
1237
1238         point eMid(edgeLine.centre());
1239
1240         treeBoundBox tightest(eMid - searchSpan, eMid + searchSpan);
1241
1242         pointIndexHit info = ppTree.findNearest
1243         (
1244             edgeLine,
1245             tightest,
1246             pointOnEdge[i]
1247         );
1248
1249         if (!info.hit())
1250         {
1251             edgeLabel[i] = -1;
1252         }
1253         else
1254         {
1255             edgeLabel[i] = selectedEdges[info.index()];
1256
1257             pointOnFeature[i] = info.hitPoint();
1258         }
1259     }
1260 }
1261
1262
1263 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
1264
1265 void Foam::surfaceFeatures::operator=(const surfaceFeatures& rhs)
1266 {
1267     // Check for assignment to self
1268     if (this == &rhs)
1269     {
1270         FatalErrorIn
1271         (
1272             "Foam::surfaceFeatures::operator=(const Foam::surfaceFeatures&)"
1273         )   << "Attempted assignment to self"
1274             << abort(FatalError);
1275     }
1276
1277     if (&surf_ != &rhs.surface())
1278     {
1279         FatalErrorIn
1280         (
1281             "Foam::surfaceFeatures::operator=(const Foam::surfaceFeatures&)"
1282         )   << "Operating on different surfaces"
1283             << abort(FatalError);
1284     }
1285
1286     featurePoints_ = rhs.featurePoints();
1287     featureEdges_ = rhs.featureEdges();
1288     externalStart_ = rhs.externalStart();
1289     internalStart_ = rhs.internalStart();
1290 }
1291
1292
1293 // ************************************************************************* //