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