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