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