surfaceFeatures.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "surfaceFeatures.H"
27 #include "triSurface.H"
28 #include "indexedOctree.H"
29 #include "treeDataEdge.H"
30 #include "treeDataPoint.H"
31 #include "meshTools.H"
32 #include "linePointRef.H"
33 #include "OFstream.H"
34 #include "IFstream.H"
35 #include "unitConversion.H"
36 #include "EdgeMap.H"
37 
38 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
39 
40 namespace Foam
41 {
42  defineTypeNameAndDebug(surfaceFeatures, 0);
43 
44  const scalar surfaceFeatures::parallelTolerance = sin(degToRad(1.0));
45 }
46 
47 
48 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
49 
50 Foam::pointIndexHit Foam::surfaceFeatures::edgeNearest
51 (
52  const point& start,
53  const point& end,
54  const point& sample
55 )
56 {
57  pointHit eHit = linePointRef(start, end).nearestDist(sample);
58 
59  // Classification of position on edge.
60  label endPoint;
61 
62  if (eHit.hit())
63  {
64  // Nearest point is on edge itself.
65  // Note: might be at or very close to endpoint. We should use tolerance
66  // here.
67  endPoint = -1;
68  }
69  else
70  {
71  // Nearest point has to be one of the end points. Find out
72  // which one.
73  if
74  (
75  mag(eHit.rawPoint() - start)
76  < mag(eHit.rawPoint() - end)
77  )
78  {
79  endPoint = 0;
80  }
81  else
82  {
83  endPoint = 1;
84  }
85  }
86 
87  return pointIndexHit(eHit.hit(), eHit.rawPoint(), endPoint);
88 }
89 
90 
91 // Go from selected edges only to a value for every edge
93  const
94 {
95  List<edgeStatus> edgeStat(surf_.nEdges(), NONE);
96 
97  // Region edges first
98  for (label i = 0; i < externalStart_; i++)
99  {
100  edgeStat[featureEdges_[i]] = REGION;
101  }
102 
103  // External edges
104  for (label i = externalStart_; i < internalStart_; i++)
105  {
106  edgeStat[featureEdges_[i]] = EXTERNAL;
107  }
108 
109  // Internal edges
110  for (label i = internalStart_; i < featureEdges_.size(); i++)
111  {
112  edgeStat[featureEdges_[i]] = INTERNAL;
113  }
114 
115  return edgeStat;
116 }
117 
118 
119 // Set from value for every edge
121 (
122  const List<edgeStatus>& edgeStat,
123  const scalar includedAngle
124 )
125 {
126  // Count
127 
128  label nRegion = 0;
129  label nExternal = 0;
130  label nInternal = 0;
131 
132  forAll(edgeStat, edgeI)
133  {
134  if (edgeStat[edgeI] == REGION)
135  {
136  nRegion++;
137  }
138  else if (edgeStat[edgeI] == EXTERNAL)
139  {
140  nExternal++;
141  }
142  else if (edgeStat[edgeI] == INTERNAL)
143  {
144  nInternal++;
145  }
146  }
147 
148  externalStart_ = nRegion;
149  internalStart_ = externalStart_ + nExternal;
150 
151 
152  // Copy
153 
154  featureEdges_.setSize(internalStart_ + nInternal);
155 
156  label regionI = 0;
157  label externalI = externalStart_;
158  label internalI = internalStart_;
159 
160  forAll(edgeStat, edgeI)
161  {
162  if (edgeStat[edgeI] == REGION)
163  {
164  featureEdges_[regionI++] = edgeI;
165  }
166  else if (edgeStat[edgeI] == EXTERNAL)
167  {
168  featureEdges_[externalI++] = edgeI;
169  }
170  else if (edgeStat[edgeI] == INTERNAL)
171  {
172  featureEdges_[internalI++] = edgeI;
173  }
174  }
175 
176  const scalar minCos = Foam::cos(degToRad(180.0 - includedAngle));
177 
178  calcFeatPoints(edgeStat, minCos);
179 }
180 
181 
182 //construct feature points where more than 2 feature edges meet
183 void Foam::surfaceFeatures::calcFeatPoints
184 (
185  const List<edgeStatus>& edgeStat,
186  const scalar minCos
187 )
188 {
189  DynamicList<label> featurePoints(surf_.nPoints()/1000);
190 
191  const labelListList& pointEdges = surf_.pointEdges();
192  const edgeList& edges = surf_.edges();
193  const pointField& localPoints = surf_.localPoints();
194 
195  forAll(pointEdges, pointi)
196  {
197  const labelList& pEdges = pointEdges[pointi];
198 
199  label nFeatEdges = 0;
200 
201  forAll(pEdges, i)
202  {
203  if (edgeStat[pEdges[i]] != NONE)
204  {
205  nFeatEdges++;
206  }
207  }
208 
209  if (nFeatEdges > 2)
210  {
211  featurePoints.append(pointi);
212  }
213  else if (nFeatEdges == 2)
214  {
215  // Check the angle between the two edges
216  DynamicList<vector> edgeVecs(2);
217 
218  forAll(pEdges, i)
219  {
220  const label edgeI = pEdges[i];
221 
222  if (edgeStat[edgeI] != NONE)
223  {
224  edgeVecs.append(edges[edgeI].vec(localPoints));
225  edgeVecs.last() /= mag(edgeVecs.last());
226  }
227  }
228 
229  if (mag(edgeVecs[0] & edgeVecs[1]) < minCos)
230  {
231  featurePoints.append(pointi);
232  }
233  }
234  }
235 
236  featurePoints_.transfer(featurePoints);
237 }
238 
239 
240 void Foam::surfaceFeatures::classifyFeatureAngles
241 (
242  const labelListList& edgeFaces,
243  List<edgeStatus>& edgeStat,
244  const scalar minCos,
245  const bool geometricTestOnly
246 ) const
247 {
248  const vectorField& faceNormals = surf_.faceNormals();
249  const pointField& points = surf_.points();
250 
251  // Special case: minCos=1
252  bool selectAll = (mag(minCos-1.0) < small);
253 
254  forAll(edgeFaces, edgeI)
255  {
256  const labelList& eFaces = edgeFaces[edgeI];
257 
258  if (eFaces.size() != 2)
259  {
260  // Non-manifold. What to do here? Is region edge? external edge?
261  edgeStat[edgeI] = REGION;
262  }
263  else
264  {
265  label face0 = eFaces[0];
266  label face1 = eFaces[1];
267 
268  if
269  (
270  !geometricTestOnly
271  && surf_[face0].region() != surf_[face1].region()
272  )
273  {
274  edgeStat[edgeI] = REGION;
275  }
276  else if
277  (
278  selectAll
279  || ((faceNormals[face0] & faceNormals[face1]) < minCos)
280  )
281  {
282  // Check if convex or concave by looking at angle
283  // between face centres and normal
284  vector f0Tof1 =
285  surf_[face1].centre(points)
286  - surf_[face0].centre(points);
287 
288  if ((f0Tof1 & faceNormals[face0]) >= 0.0)
289  {
290  edgeStat[edgeI] = INTERNAL;
291  }
292  else
293  {
294  edgeStat[edgeI] = EXTERNAL;
295  }
296  }
297  }
298  }
299 }
300 
301 
302 // Returns next feature edge connected to pointi with correct value.
303 Foam::label Foam::surfaceFeatures::nextFeatEdge
304 (
305  const List<edgeStatus>& edgeStat,
306  const labelList& featVisited,
307  const label unsetVal,
308  const label prevEdgeI,
309  const label vertI
310 ) const
311 {
312  const labelList& pEdges = surf_.pointEdges()[vertI];
313 
314  label nextEdgeI = -1;
315 
316  forAll(pEdges, i)
317  {
318  label edgeI = pEdges[i];
319 
320  if
321  (
322  edgeI != prevEdgeI
323  && edgeStat[edgeI] != NONE
324  && featVisited[edgeI] == unsetVal
325  )
326  {
327  if (nextEdgeI == -1)
328  {
329  nextEdgeI = edgeI;
330  }
331  else
332  {
333  // More than one feature edge to choose from. End of segment.
334  return -1;
335  }
336  }
337  }
338 
339  return nextEdgeI;
340 }
341 
342 
343 // Finds connected feature edges by walking from prevEdgeI in direction of
344 // prevPointi. Marks feature edges visited in featVisited by assigning them
345 // the current feature line number. Returns cumulative length of edges walked.
346 // Works in one of two modes:
347 // - mark : step to edges with featVisited = -1.
348 // Mark edges visited with currentFeatI.
349 // - clear : step to edges with featVisited = currentFeatI
350 // Mark edges visited with -2 and erase from feature edges.
351 Foam::surfaceFeatures::labelScalar Foam::surfaceFeatures::walkSegment
352 (
353  const bool mark,
354  const List<edgeStatus>& edgeStat,
355  const label startEdgeI,
356  const label startPointi,
357  const label currentFeatI,
358  labelList& featVisited
359 )
360 {
361  label edgeI = startEdgeI;
362 
363  label vertI = startPointi;
364 
365  scalar visitedLength = 0.0;
366 
367  label nVisited = 0;
368 
369  if (findIndex(featurePoints_, startPointi) >= 0)
370  {
371  // Do not walk across feature points
372 
373  return labelScalar(nVisited, visitedLength);
374  }
375 
376  //
377  // Now we have:
378  // edgeI : first edge on this segment
379  // vertI : one of the endpoints of this segment
380  //
381  // Start walking, marking off edges (in featVisited)
382  // as we go along.
383  //
384 
385  label unsetVal;
386 
387  if (mark)
388  {
389  unsetVal = -1;
390  }
391  else
392  {
393  unsetVal = currentFeatI;
394  }
395 
396  do
397  {
398  // Step to next feature edge with value unsetVal
399  edgeI = nextFeatEdge(edgeStat, featVisited, unsetVal, edgeI, vertI);
400 
401  if (edgeI == -1 || edgeI == startEdgeI)
402  {
403  break;
404  }
405 
406  // Mark with current value. If in clean mode also remove feature edge.
407 
408  if (mark)
409  {
410  featVisited[edgeI] = currentFeatI;
411  }
412  else
413  {
414  featVisited[edgeI] = -2;
415  }
416 
417 
418  // Step to next vertex on edge
419  const edge& e = surf_.edges()[edgeI];
420 
421  vertI = e.otherVertex(vertI);
422 
423 
424  // Update cumulative length
425  visitedLength += e.mag(surf_.localPoints());
426 
427  nVisited++;
428 
429  if (nVisited > surf_.nEdges())
430  {
431  Warning<< "walkSegment : reached iteration limit in walking "
432  << "feature edges on surface from edge:" << startEdgeI
433  << " vertex:" << startPointi << nl
434  << "Returning with large length" << endl;
435 
436  return labelScalar(nVisited, great);
437  }
438  }
439  while (true);
440 
441  return labelScalar(nVisited, visitedLength);
442 }
443 
444 
445 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
446 
448 :
449  surf_(surf),
450  featurePoints_(0),
451  featureEdges_(0),
452  externalStart_(0),
453  internalStart_(0)
454 {}
455 
456 
457 // Construct from components.
459 (
460  const triSurface& surf,
461  const labelList& featurePoints,
462  const labelList& featureEdges,
463  const label externalStart,
464  const label internalStart
465 )
466 :
467  surf_(surf),
468  featurePoints_(featurePoints),
469  featureEdges_(featureEdges),
470  externalStart_(externalStart),
471  internalStart_(externalStart)
472 {}
473 
474 
475 // Construct from surface, angle and min length
477 (
478  const triSurface& surf,
479  const scalar includedAngle,
480  const scalar minLen,
481  const label minElems,
482  const bool geometricTestOnly
483 )
484 :
485  surf_(surf),
486  featurePoints_(0),
487  featureEdges_(0),
488  externalStart_(0),
489  internalStart_(0)
490 {
491  findFeatures(includedAngle, geometricTestOnly);
492 
493  if (minLen > 0 || minElems > 0)
494  {
495  trimFeatures(minLen, minElems, includedAngle);
496  }
497 }
498 
499 
501 (
502  const triSurface& surf,
503  const dictionary& featInfoDict
504 )
505 :
506  surf_(surf),
507  featurePoints_(featInfoDict.lookup("featurePoints")),
508  featureEdges_(featInfoDict.lookup("featureEdges")),
509  externalStart_(featInfoDict.lookup<label>("externalStart")),
510  internalStart_(featInfoDict.lookup<label>("internalStart"))
511 {}
512 
513 
515 (
516  const triSurface& surf,
517  const fileName& fName
518 )
519 :
520  surf_(surf),
521  featurePoints_(0),
522  featureEdges_(0),
523  externalStart_(0),
524  internalStart_(0)
525 {
526  IFstream str(fName);
527 
528  dictionary featInfoDict(str);
529 
530  featureEdges_ = labelList(featInfoDict.lookup("featureEdges"));
531  featurePoints_ = labelList(featInfoDict.lookup("featurePoints"));
532  externalStart_ = featInfoDict.lookup<label>("externalStart");
533  internalStart_ = featInfoDict.lookup<label>("internalStart");
534 }
535 
536 
538 (
539  const triSurface& surf,
540  const pointField& points,
541  const edgeList& edges,
542  const scalar mergeTol,
543  const bool geometricTestOnly
544 )
545 :
546  surf_(surf),
547  featurePoints_(0),
548  featureEdges_(0),
549  externalStart_(0),
550  internalStart_(0)
551 {
552  // Match edge mesh edges with the triSurface edges
553 
554  const labelListList& surfEdgeFaces = surf_.edgeFaces();
555  const edgeList& surfEdges = surf_.edges();
556 
557  scalar mergeTolSqr = sqr(mergeTol);
558 
559  EdgeMap<label> dynFeatEdges(2*edges.size());
560  DynamicList<labelList> dynFeatureEdgeFaces(edges.size());
561 
562  labelList edgeLabel;
563 
565  (
566  edges,
567  points,
568  mergeTolSqr,
569  edgeLabel // label of surface edge or -1
570  );
571 
572  label count = 0;
573  forAll(edgeLabel, sEdgeI)
574  {
575  const label sEdge = edgeLabel[sEdgeI];
576 
577  if (sEdge == -1)
578  {
579  continue;
580  }
581 
582  dynFeatEdges.insert(surfEdges[sEdge], count++);
583  dynFeatureEdgeFaces.append(surfEdgeFaces[sEdge]);
584  }
585 
586  // Find whether an edge is external or internal
587  List<edgeStatus> edgeStat(dynFeatEdges.size(), NONE);
588 
589  classifyFeatureAngles
590  (
591  dynFeatureEdgeFaces,
592  edgeStat,
593  great,
594  geometricTestOnly
595  );
596 
597  // Transfer the edge status to a list encompassing all edges in the surface
598  // so that calcFeatPoints can be used.
599  List<edgeStatus> allEdgeStat(surf_.nEdges(), NONE);
600 
601  forAll(allEdgeStat, eI)
602  {
603  EdgeMap<label>::const_iterator iter = dynFeatEdges.find(surfEdges[eI]);
604 
605  if (iter != dynFeatEdges.end())
606  {
607  allEdgeStat[eI] = edgeStat[iter()];
608  }
609  }
610 
611  edgeStat.clear();
612  dynFeatEdges.clear();
613 
614  setFromStatus(allEdgeStat, great);
615 }
616 
617 
619 :
620  surf_(sf.surface()),
621  featurePoints_(sf.featurePoints()),
622  featureEdges_(sf.featureEdges()),
623  externalStart_(sf.externalStart()),
624  internalStart_(sf.internalStart())
625 {}
626 
627 
628 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
629 
631 (
632  const bool regionEdges,
633  const bool externalEdges,
634  const bool internalEdges
635 ) const
636 {
637  DynamicList<label> selectedEdges;
638 
639  if (regionEdges)
640  {
641  selectedEdges.setCapacity(selectedEdges.size() + nRegionEdges());
642 
643  for (label i = 0; i < externalStart_; i++)
644  {
645  selectedEdges.append(featureEdges_[i]);
646  }
647  }
648 
649  if (externalEdges)
650  {
651  selectedEdges.setCapacity(selectedEdges.size() + nExternalEdges());
652 
653  for (label i = externalStart_; i < internalStart_; i++)
654  {
655  selectedEdges.append(featureEdges_[i]);
656  }
657  }
658 
659  if (internalEdges)
660  {
661  selectedEdges.setCapacity(selectedEdges.size() + nInternalEdges());
662 
663  for (label i = internalStart_; i < featureEdges_.size(); i++)
664  {
665  selectedEdges.append(featureEdges_[i]);
666  }
667  }
668 
669  return selectedEdges.shrink();
670 }
671 
672 
674 (
675  const scalar includedAngle,
676  const bool geometricTestOnly
677 )
678 {
679  scalar minCos = Foam::cos(degToRad(180.0 - includedAngle));
680 
681  // Per edge whether is feature edge.
682  List<edgeStatus> edgeStat(surf_.nEdges(), NONE);
683 
684  classifyFeatureAngles
685  (
686  surf_.edgeFaces(),
687  edgeStat,
688  minCos,
689  geometricTestOnly
690  );
691 
692  setFromStatus(edgeStat, includedAngle);
693 }
694 
695 
696 // Remove small strings of edges. First assigns string number to
697 // every edge and then works out the length of them.
699 (
700  const scalar minLen,
701  const label minElems,
702  const scalar includedAngle
703 )
704 {
705  // Convert feature edge list to status per edge.
706  List<edgeStatus> edgeStat(toStatus());
707 
708 
709  // Mark feature edges according to connected lines.
710  // -1: unassigned
711  // -2: part of too small a feature line
712  // >0: feature line number
713  labelList featLines(surf_.nEdges(), -1);
714 
715  // Current featureline number.
716  label featI = 0;
717 
718  // Current starting edge
719  label startEdgeI = 0;
720 
721  do
722  {
723  // Find unset featureline
724  for (; startEdgeI < edgeStat.size(); startEdgeI++)
725  {
726  if
727  (
728  edgeStat[startEdgeI] != NONE
729  && featLines[startEdgeI] == -1
730  )
731  {
732  // Found unassigned feature edge.
733  break;
734  }
735  }
736 
737  if (startEdgeI == edgeStat.size())
738  {
739  // No unset feature edge found. All feature edges have line number
740  // assigned.
741  break;
742  }
743 
744  featLines[startEdgeI] = featI;
745 
746  const edge& startEdge = surf_.edges()[startEdgeI];
747 
748  // Walk 'left' and 'right' from startEdge.
749  labelScalar leftPath =
750  walkSegment
751  (
752  true, // 'mark' mode
753  edgeStat,
754  startEdgeI, // edge, not yet assigned to a featureLine
755  startEdge[0], // start of edge
756  featI, // Mark value
757  featLines // Mark for all edges
758  );
759 
760  labelScalar rightPath =
761  walkSegment
762  (
763  true,
764  edgeStat,
765  startEdgeI,
766  startEdge[1], // end of edge
767  featI,
768  featLines
769  );
770 
771  if
772  (
773  (
774  leftPath.len_
775  + rightPath.len_
776  + startEdge.mag(surf_.localPoints())
777  < minLen
778  )
779  || (leftPath.n_ + rightPath.n_ + 1 < minElems)
780  )
781  {
782  // Rewalk same route (recognisable by featLines == featI)
783  // to reset featLines.
784 
785  featLines[startEdgeI] = -2;
786 
787  walkSegment
788  (
789  false, // 'clean' mode
790  edgeStat,
791  startEdgeI, // edge, not yet assigned to a featureLine
792  startEdge[0], // start of edge
793  featI, // Unset value
794  featLines // Mark for all edges
795  );
796 
797  walkSegment
798  (
799  false,
800  edgeStat,
801  startEdgeI,
802  startEdge[1], // end of edge
803  featI,
804  featLines
805  );
806  }
807  else
808  {
809  featI++;
810  }
811  }
812  while (true);
813 
814  // Unmark all feature lines that have featLines=-2
815  forAll(featureEdges_, i)
816  {
817  label edgeI = featureEdges_[i];
818 
819  if (featLines[edgeI] == -2)
820  {
821  edgeStat[edgeI] = NONE;
822  }
823  }
824 
825  // Convert back to edge labels
826  setFromStatus(edgeStat, includedAngle);
827 
828  return featLines;
829 }
830 
831 
833 {
834 
835  dictionary featInfoDict;
836  featInfoDict.add("externalStart", externalStart_);
837  featInfoDict.add("internalStart", internalStart_);
838  featInfoDict.add("featureEdges", featureEdges_);
839  featInfoDict.add("featurePoints", featurePoints_);
840 
841  featInfoDict.write(writeFile);
842 }
843 
844 
845 void Foam::surfaceFeatures::write(const fileName& fName) const
846 {
847  OFstream str(fName);
848 
849  writeDict(str);
850 }
851 
852 
853 void Foam::surfaceFeatures::writeObj(const fileName& prefix) const
854 {
855  OFstream regionStr(prefix + "_regionEdges.obj");
856  Pout<< "Writing region edges to " << regionStr.name() << endl;
857 
858  label verti = 0;
859  for (label i = 0; i < externalStart_; i++)
860  {
861  const edge& e = surf_.edges()[featureEdges_[i]];
862 
863  meshTools::writeOBJ(regionStr, surf_.localPoints()[e[0]]); verti++;
864  meshTools::writeOBJ(regionStr, surf_.localPoints()[e[1]]); verti++;
865  regionStr << "l " << verti-1 << ' ' << verti << endl;
866  }
867 
868 
869  OFstream externalStr(prefix + "_externalEdges.obj");
870  Pout<< "Writing external edges to " << externalStr.name() << endl;
871 
872  verti = 0;
873  for (label i = externalStart_; i < internalStart_; i++)
874  {
875  const edge& e = surf_.edges()[featureEdges_[i]];
876 
877  meshTools::writeOBJ(externalStr, surf_.localPoints()[e[0]]); verti++;
878  meshTools::writeOBJ(externalStr, surf_.localPoints()[e[1]]); verti++;
879  externalStr << "l " << verti-1 << ' ' << verti << endl;
880  }
881 
882  OFstream internalStr(prefix + "_internalEdges.obj");
883  Pout<< "Writing internal edges to " << internalStr.name() << endl;
884 
885  verti = 0;
886  for (label i = internalStart_; i < featureEdges_.size(); i++)
887  {
888  const edge& e = surf_.edges()[featureEdges_[i]];
889 
890  meshTools::writeOBJ(internalStr, surf_.localPoints()[e[0]]); verti++;
891  meshTools::writeOBJ(internalStr, surf_.localPoints()[e[1]]); verti++;
892  internalStr << "l " << verti-1 << ' ' << verti << endl;
893  }
894 
895  OFstream pointStr(prefix + "_points.obj");
896  Pout<< "Writing feature points to " << pointStr.name() << endl;
897 
898  forAll(featurePoints_, i)
899  {
900  label pointi = featurePoints_[i];
901 
902  meshTools::writeOBJ(pointStr, surf_.localPoints()[pointi]);
903  }
904 }
905 
906 
907 // Get nearest vertex on patch for every point of surf in pointSet.
909 (
910  const labelList& pointLabels,
911  const pointField& samples,
912  const scalarField& maxDistSqr
913 ) const
914 {
915  // Build tree out of all samples.
916 
917  // Note: cannot be done one the fly - gcc4.4 compiler bug.
918  treeBoundBox bb(samples);
919 
921  (
922  treeDataPoint(samples), // all information needed to do checks
923  bb, // overall search domain
924  8, // maxLevel
925  10, // leafsize
926  3.0 // duplicity
927  );
928 
929  // From patch point to surface point
930  Map<label> nearest(2*pointLabels.size());
931 
932  const pointField& surfPoints = surf_.localPoints();
933 
934  forAll(pointLabels, i)
935  {
936  label surfPointi = pointLabels[i];
937 
938  const point& surfPt = surfPoints[surfPointi];
939 
940  pointIndexHit info = ppTree.findNearest
941  (
942  surfPt,
943  maxDistSqr[i]
944  );
945 
946  if (!info.hit())
947  {
949  << "Problem for point "
950  << surfPointi << " in tree " << ppTree.bb()
951  << abort(FatalError);
952  }
953 
954  label sampleI = info.index();
955 
956  if (magSqr(samples[sampleI] - surfPt) < maxDistSqr[sampleI])
957  {
958  nearest.insert(sampleI, surfPointi);
959  }
960  }
961 
962 
963  if (debug)
964  {
965  //
966  // Dump to obj file
967  //
968 
969  Pout<< endl
970  << "Dumping nearest surface feature points to nearestSamples.obj"
971  << endl
972  << "View this Lightwave-OBJ file with e.g. javaview" << endl
973  << endl;
974 
975  OFstream objStream("nearestSamples.obj");
976 
977  label vertI = 0;
978  forAllConstIter(Map<label>, nearest, iter)
979  {
980  meshTools::writeOBJ(objStream, samples[iter.key()]); vertI++;
981  meshTools::writeOBJ(objStream, surfPoints[iter()]); vertI++;
982  objStream<< "l " << vertI-1 << ' ' << vertI << endl;
983  }
984  }
985 
986  return nearest;
987 }
988 
989 
990 // Get nearest sample point for regularly sampled points along
991 // selected edges. Return map from sample to edge label
993 (
994  const labelList& selectedEdges,
995  const pointField& samples,
996  const scalarField& sampleDist,
997  const scalarField& maxDistSqr,
998  const scalar minSampleDist
999 ) const
1000 {
1001  const pointField& surfPoints = surf_.localPoints();
1002  const edgeList& surfEdges = surf_.edges();
1003 
1004  scalar maxSearchSqr = max(maxDistSqr);
1005 
1006  // Note: cannot be done one the fly - gcc4.4 compiler bug.
1007  treeBoundBox bb(samples);
1008 
1010  (
1011  treeDataPoint(samples), // all information needed to do checks
1012  bb, // overall search domain
1013  8, // maxLevel
1014  10, // leafsize
1015  3.0 // duplicity
1016  );
1017 
1018  // From patch point to surface edge.
1019  Map<label> nearest(2*selectedEdges.size());
1020 
1021  forAll(selectedEdges, i)
1022  {
1023  label surfEdgeI = selectedEdges[i];
1024 
1025  const edge& e = surfEdges[surfEdgeI];
1026 
1027  if (debug && (i % 1000) == 0)
1028  {
1029  Pout<< "looking at surface feature edge " << surfEdgeI
1030  << " verts:" << e << " points:" << surfPoints[e[0]]
1031  << ' ' << surfPoints[e[1]] << endl;
1032  }
1033 
1034  // Normalised edge vector
1035  vector eVec = e.vec(surfPoints);
1036  scalar eMag = mag(eVec);
1037  eVec /= eMag;
1038 
1039 
1040  //
1041  // Sample along edge
1042  //
1043 
1044  bool exit = false;
1045 
1046  // Coordinate along edge (0 .. eMag)
1047  scalar s = 0.0;
1048 
1049  while (true)
1050  {
1051  point edgePoint(surfPoints[e.start()] + s*eVec);
1052 
1053  pointIndexHit info = ppTree.findNearest
1054  (
1055  edgePoint,
1056  maxSearchSqr
1057  );
1058 
1059  if (!info.hit())
1060  {
1061  // No point close enough to surface edge.
1062  break;
1063  }
1064 
1065  label sampleI = info.index();
1066 
1067  if (magSqr(info.hitPoint() - edgePoint) < maxDistSqr[sampleI])
1068  {
1069  nearest.insert(sampleI, surfEdgeI);
1070  }
1071 
1072  if (exit)
1073  {
1074  break;
1075  }
1076 
1077  // Step to next sample point using local distance.
1078  // Truncate to max 1/minSampleDist samples per feature edge.
1079  s += max(minSampleDist*eMag, sampleDist[sampleI]);
1080 
1081  if (s >= (1-minSampleDist)*eMag)
1082  {
1083  // Do one more sample, at edge endpoint
1084  s = eMag;
1085  exit = true;
1086  }
1087  }
1088  }
1089 
1090 
1091 
1092  if (debug)
1093  {
1094  // Dump to obj file
1095 
1096  Pout<< "Dumping nearest surface edges to nearestEdges.obj\n"
1097  << "View this Lightwave-OBJ file with e.g. javaview\n" << endl;
1098 
1099  OFstream objStream("nearestEdges.obj");
1100 
1101  label vertI = 0;
1102  forAllConstIter(Map<label>, nearest, iter)
1103  {
1104  const label sampleI = iter.key();
1105 
1106  meshTools::writeOBJ(objStream, samples[sampleI]); vertI++;
1107 
1108  const edge& e = surfEdges[iter()];
1109 
1110  point nearPt =
1111  e.line(surfPoints).nearestDist(samples[sampleI]).rawPoint();
1112 
1113  meshTools::writeOBJ(objStream, nearPt); vertI++;
1114 
1115  objStream<< "l " << vertI-1 << ' ' << vertI << endl;
1116  }
1117  }
1118 
1119  return nearest;
1120 }
1121 
1122 
1123 // Get nearest edge on patch for regularly sampled points along the feature
1124 // edges. Return map from patch edge to feature edge.
1125 //
1126 // Q: using point-based sampleDist and maxDist (distance to look around
1127 // each point). Should they be edge-based e.g. by averaging or max()?
1130  const labelList& selectedEdges,
1131  const edgeList& sampleEdges,
1132  const labelList& selectedSampleEdges,
1133  const pointField& samplePoints,
1134  const scalarField& sampleDist,
1135  const scalarField& maxDistSqr,
1136  const scalar minSampleDist
1137 ) const
1138 {
1139  // Build tree out of selected sample edges.
1141  (
1142  treeDataEdge
1143  (
1144  false,
1145  sampleEdges,
1146  samplePoints,
1147  selectedSampleEdges
1148  ), // geometric info container for edges
1149  treeBoundBox(samplePoints), // overall search domain
1150  8, // maxLevel
1151  10, // leafsize
1152  3.0 // duplicity
1153  );
1154 
1155  const pointField& surfPoints = surf_.localPoints();
1156  const edgeList& surfEdges = surf_.edges();
1157 
1158  scalar maxSearchSqr = max(maxDistSqr);
1159 
1160  Map<pointIndexHit> nearest(2*sampleEdges.size());
1161 
1162  //
1163  // Loop over all selected edges. Sample at regular intervals. Find nearest
1164  // sampleEdges (using octree)
1165  //
1166 
1167  forAll(selectedEdges, i)
1168  {
1169  label surfEdgeI = selectedEdges[i];
1170 
1171  const edge& e = surfEdges[surfEdgeI];
1172 
1173  if (debug && (i % 1000) == 0)
1174  {
1175  Pout<< "looking at surface feature edge " << surfEdgeI
1176  << " verts:" << e << " points:" << surfPoints[e[0]]
1177  << ' ' << surfPoints[e[1]] << endl;
1178  }
1179 
1180  // Normalised edge vector
1181  vector eVec = e.vec(surfPoints);
1182  scalar eMag = mag(eVec);
1183  eVec /= eMag;
1184 
1185 
1186  //
1187  // Sample along edge
1188  //
1189 
1190  bool exit = false;
1191 
1192  // Coordinate along edge (0 .. eMag)
1193  scalar s = 0.0;
1194 
1195  while (true)
1196  {
1197  point edgePoint(surfPoints[e.start()] + s*eVec);
1198 
1199  pointIndexHit info = ppTree.findNearest
1200  (
1201  edgePoint,
1202  maxSearchSqr
1203  );
1204 
1205  if (!info.hit())
1206  {
1207  // No edge close enough to surface edge.
1208  break;
1209  }
1210 
1211  label index = info.index();
1212 
1213  label sampleEdgeI = ppTree.shapes().edgeLabels()[index];
1214 
1215  const edge& e = sampleEdges[sampleEdgeI];
1216 
1217  if (magSqr(info.hitPoint() - edgePoint) < maxDistSqr[e.start()])
1218  {
1219  nearest.insert
1220  (
1221  sampleEdgeI,
1222  pointIndexHit(true, edgePoint, surfEdgeI)
1223  );
1224  }
1225 
1226  if (exit)
1227  {
1228  break;
1229  }
1230 
1231  // Step to next sample point using local distance.
1232  // Truncate to max 1/minSampleDist samples per feature edge.
1233  // s += max(minSampleDist*eMag, sampleDist[e.start()]);
1234  s += 0.01*eMag;
1235 
1236  if (s >= (1-minSampleDist)*eMag)
1237  {
1238  // Do one more sample, at edge endpoint
1239  s = eMag;
1240  exit = true;
1241  }
1242  }
1243  }
1244 
1245 
1246  if (debug)
1247  {
1248  // Dump to obj file
1249 
1250  Pout<< "Dumping nearest surface feature edges to nearestEdges.obj\n"
1251  << "View this Lightwave-OBJ file with e.g. javaview\n" << endl;
1252 
1253  OFstream objStream("nearestEdges.obj");
1254 
1255  label vertI = 0;
1256  forAllConstIter(Map<pointIndexHit>, nearest, iter)
1257  {
1258  const label sampleEdgeI = iter.key();
1259 
1260  const edge& sampleEdge = sampleEdges[sampleEdgeI];
1261 
1262  // Write line from edgeMid to point on feature edge
1263  meshTools::writeOBJ(objStream, sampleEdge.centre(samplePoints));
1264  vertI++;
1265 
1266  meshTools::writeOBJ(objStream, iter().rawPoint());
1267  vertI++;
1268 
1269  objStream<< "l " << vertI-1 << ' ' << vertI << endl;
1270  }
1271  }
1272 
1273  return nearest;
1274 }
1275 
1276 
1277 // Get nearest surface edge for every sample. Return in form of
1278 // labelLists giving surfaceEdge label&intersectionpoint.
1281  const labelList& selectedEdges,
1282  const pointField& samples,
1283  scalar searchSpanSqr, // Search span
1284  labelList& edgeLabel,
1285  labelList& edgeEndPoint,
1286  pointField& edgePoint
1287 ) const
1288 {
1289  edgeLabel.setSize(samples.size());
1290  edgeEndPoint.setSize(samples.size());
1291  edgePoint.setSize(samples.size());
1292 
1293  const pointField& localPoints = surf_.localPoints();
1294 
1295  treeBoundBox searchDomain(localPoints);
1296  searchDomain.inflate(0.1);
1297 
1299  (
1300  treeDataEdge
1301  (
1302  false,
1303  surf_.edges(),
1304  localPoints,
1305  selectedEdges
1306  ), // all information needed to do geometric checks
1307  searchDomain, // overall search domain
1308  8, // maxLevel
1309  10, // leafsize
1310  3.0 // duplicity
1311  );
1312 
1313  forAll(samples, i)
1314  {
1315  const point& sample = samples[i];
1316 
1317  pointIndexHit info = ppTree.findNearest
1318  (
1319  sample,
1320  searchSpanSqr
1321  );
1322 
1323  if (!info.hit())
1324  {
1325  edgeLabel[i] = -1;
1326  }
1327  else
1328  {
1329  edgeLabel[i] = selectedEdges[info.index()];
1330 
1331  // Need to recalculate to classify edgeEndPoint
1332  const edge& e = surf_.edges()[edgeLabel[i]];
1333 
1334  pointIndexHit pHit = edgeNearest
1335  (
1336  localPoints[e.start()],
1337  localPoints[e.end()],
1338  sample
1339  );
1340 
1341  edgePoint[i] = pHit.rawPoint();
1342  edgeEndPoint[i] = pHit.index();
1343  }
1344  }
1345 }
1346 
1347 
1348 // Get nearest point on nearest feature edge for every sample (is edge)
1351  const labelList& selectedEdges,
1352 
1353  const edgeList& sampleEdges,
1354  const labelList& selectedSampleEdges,
1355  const pointField& samplePoints,
1356 
1357  const vector& searchSpan, // Search span
1358  labelList& edgeLabel, // label of surface edge or -1
1359  pointField& pointOnEdge, // point on above edge
1360  pointField& pointOnFeature // point on sample edge
1361 ) const
1362 {
1363  edgeLabel.setSize(selectedSampleEdges.size());
1364  pointOnEdge.setSize(selectedSampleEdges.size());
1365  pointOnFeature.setSize(selectedSampleEdges.size());
1366 
1367  treeBoundBox searchDomain(surf_.localPoints());
1368 
1370  (
1371  treeDataEdge
1372  (
1373  false,
1374  surf_.edges(),
1375  surf_.localPoints(),
1376  selectedEdges
1377  ), // all information needed to do geometric checks
1378  searchDomain, // overall search domain
1379  8, // maxLevel
1380  10, // leafsize
1381  3.0 // duplicity
1382  );
1383 
1384  forAll(selectedSampleEdges, i)
1385  {
1386  const edge& e = sampleEdges[selectedSampleEdges[i]];
1387 
1388  linePointRef edgeLine = e.line(samplePoints);
1389 
1390  point eMid(edgeLine.centre());
1391 
1392  treeBoundBox tightest(eMid - searchSpan, eMid + searchSpan);
1393 
1394  pointIndexHit info = ppTree.findNearest
1395  (
1396  edgeLine,
1397  tightest,
1398  pointOnEdge[i]
1399  );
1400 
1401  if (!info.hit())
1402  {
1403  edgeLabel[i] = -1;
1404  }
1405  else
1406  {
1407  edgeLabel[i] = selectedEdges[info.index()];
1408 
1409  pointOnFeature[i] = info.hitPoint();
1410  }
1411  }
1412 }
1413 
1414 
1417  const edgeList& edges,
1418  const pointField& points,
1419  scalar searchSpanSqr, // Search span
1420  labelList& edgeLabel
1421 ) const
1422 {
1423  edgeLabel = labelList(surf_.nEdges(), -1);
1424 
1425  treeBoundBox searchDomain(points);
1426  searchDomain.inflate(0.1);
1427 
1429  (
1430  treeDataEdge
1431  (
1432  false,
1433  edges,
1434  points,
1435  identity(edges.size())
1436  ), // all information needed to do geometric checks
1437  searchDomain, // overall search domain
1438  8, // maxLevel
1439  10, // leafsize
1440  3.0 // duplicity
1441  );
1442 
1443  const edgeList& surfEdges = surf_.edges();
1444  const pointField& surfLocalPoints = surf_.localPoints();
1445 
1446  forAll(surfEdges, edgeI)
1447  {
1448  const edge& sample = surfEdges[edgeI];
1449 
1450  const point& startPoint = surfLocalPoints[sample.start()];
1451  const point& midPoint = sample.centre(surfLocalPoints);
1452 
1453  pointIndexHit infoMid = ppTree.findNearest
1454  (
1455  midPoint,
1456  searchSpanSqr
1457  );
1458 
1459  if (infoMid.hit())
1460  {
1461  const vector surfEdgeDir = midPoint - startPoint;
1462 
1463  const edge& featEdge = edges[infoMid.index()];
1464  const vector featEdgeDir = featEdge.vec(points);
1465 
1466  // Check that the edges are nearly parallel
1467  if (mag(surfEdgeDir ^ featEdgeDir) < parallelTolerance)
1468  {
1469  edgeLabel[edgeI] = edgeI;
1470  }
1471  }
1472  }
1473 }
1474 
1475 
1476 // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
1477 
1479 {
1480  // Check for assignment to self
1481  if (this == &rhs)
1482  {
1484  << "Attempted assignment to self"
1485  << abort(FatalError);
1486  }
1487 
1488  if (&surf_ != &rhs.surface())
1489  {
1491  << "Operating on different surfaces"
1492  << abort(FatalError);
1493  }
1494 
1495  featurePoints_ = rhs.featurePoints();
1496  featureEdges_ = rhs.featureEdges();
1497  externalStart_ = rhs.externalStart();
1498  internalStart_ = rhs.internalStart();
1499 }
1500 
1501 
1502 // * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * //
1503 
1504 void Foam::selectBox
1506  const triSurface& surf,
1507  const boundBox& bb,
1508  const bool inside,
1510 )
1511 {
1512  forAll(edgeStat, edgei)
1513  {
1514  const point eMid = surf.edges()[edgei].centre(surf.localPoints());
1515 
1516  if (!inside ? bb.contains(eMid) : !bb.contains(eMid))
1517  {
1518  edgeStat[edgei] = surfaceFeatures::NONE;
1519  }
1520  }
1521 }
1522 
1523 
1526  const triSurface& surf,
1527  const plane& cutPlane,
1529 )
1530 {
1531  const pointField& points = surf.points();
1532  const labelList& meshPoints = surf.meshPoints();
1533 
1534  forAll(edgeStat, edgei)
1535  {
1536  const edge& e = surf.edges()[edgei];
1537  const point& p0 = points[meshPoints[e.start()]];
1538  const point& p1 = points[meshPoints[e.end()]];
1539  const linePointRef line(p0, p1);
1540 
1541  // If edge does not intersect the plane, delete.
1542  const scalar intersect = cutPlane.lineIntersect(line);
1543  if (intersect < 0 || intersect > 1)
1544  {
1545  edgeStat[edgei] = surfaceFeatures::NONE;
1546  }
1547  }
1548 }
1549 
1550 
1553  const triSurface& surf,
1554  const scalar tol,
1555  const scalar includedAngle,
1556  const label edgei
1557 )
1558 {
1559  const edge& e = surf.edges()[edgei];
1560  const labelList& eFaces = surf.edgeFaces()[edgei];
1561 
1562  // Bin according to normal
1563 
1564  DynamicList<Foam::vector> normals(2);
1565  DynamicList<labelList> bins(2);
1566 
1567  forAll(eFaces, eFacei)
1568  {
1569  const Foam::vector& n = surf.faceNormals()[eFaces[eFacei]];
1570 
1571  // Find the normal in normals
1572  label index = -1;
1573  forAll(normals, normalI)
1574  {
1575  if (mag(n&normals[normalI]) > (1-tol))
1576  {
1577  index = normalI;
1578  break;
1579  }
1580  }
1581 
1582  if (index != -1)
1583  {
1584  bins[index].append(eFacei);
1585  }
1586  else if (normals.size() >= 2)
1587  {
1588  // Would be third normal. Mark as feature.
1589  // Pout<< "** at edge:" << surf.localPoints()[e[0]]
1590  // << surf.localPoints()[e[1]]
1591  // << " have normals:" << normals
1592  // << " and " << n << endl;
1593  return surfaceFeatures::REGION;
1594  }
1595  else
1596  {
1597  normals.append(n);
1598  bins.append(labelList(1, eFacei));
1599  }
1600  }
1601 
1602  // Check resulting number of bins
1603  if (bins.size() == 1)
1604  {
1605  // Note: should check here whether they are two sets of faces
1606  // that are planar or indeed 4 faces al coming together at an edge.
1607  // Pout<< "** at edge:"
1608  // << surf.localPoints()[e[0]]
1609  // << surf.localPoints()[e[1]]
1610  // << " have single normal:" << normals[0]
1611  // << endl;
1612  return surfaceFeatures::NONE;
1613  }
1614  else
1615  {
1616  // Two bins. Check if normals make an angle
1617 
1618  // Pout<< "** at edge:"
1619  // << surf.localPoints()[e[0]]
1620  // << surf.localPoints()[e[1]] << nl
1621  // << " normals:" << normals << nl
1622  // << " bins :" << bins << nl
1623  // << endl;
1624 
1625  if (includedAngle >= 0)
1626  {
1627  scalar minCos = Foam::cos(degToRad(180.0 - includedAngle));
1628 
1629  forAll(eFaces, i)
1630  {
1631  const Foam::vector& ni = surf.faceNormals()[eFaces[i]];
1632  for (label j=i+1; j<eFaces.size(); j++)
1633  {
1634  const Foam::vector& nj = surf.faceNormals()[eFaces[j]];
1635  if (mag(ni & nj) < minCos)
1636  {
1637  // Pout<< "have sharp feature between normal:" << ni
1638  // << " and " << nj << endl;
1639 
1640  // Is feature. Keep as region or convert to
1641  // feature angle? For now keep as region.
1642  return surfaceFeatures::REGION;
1643  }
1644  }
1645  }
1646  }
1647 
1648 
1649  // So now we have two normals bins but need to make sure both
1650  // bins have the same regions in it.
1651 
1652  // 1. store + or - region number depending
1653  // on orientation of triangle in bins[0]
1654  const labelList& bin0 = bins[0];
1655  labelList regionAndNormal(bin0.size());
1656  forAll(bin0, i)
1657  {
1658  const labelledTri& t = surf.localFaces()[eFaces[bin0[i]]];
1659  int dir = t.edgeDirection(e);
1660 
1661  if (dir > 0)
1662  {
1663  regionAndNormal[i] = t.region()+1;
1664  }
1665  else if (dir == 0)
1666  {
1668  << exit(FatalError);
1669  }
1670  else
1671  {
1672  regionAndNormal[i] = -(t.region()+1);
1673  }
1674  }
1675 
1676  // 2. check against bin1
1677  const labelList& bin1 = bins[1];
1678  labelList regionAndNormal1(bin1.size());
1679  forAll(bin1, i)
1680  {
1681  const labelledTri& t = surf.localFaces()[eFaces[bin1[i]]];
1682  int dir = t.edgeDirection(e);
1683 
1684  label myRegionAndNormal;
1685  if (dir > 0)
1686  {
1687  myRegionAndNormal = t.region()+1;
1688  }
1689  else
1690  {
1691  myRegionAndNormal = -(t.region()+1);
1692  }
1693 
1694  regionAndNormal1[i] = myRegionAndNormal;
1695 
1696  label index = findIndex(regionAndNormal, -myRegionAndNormal);
1697  if (index == -1)
1698  {
1699  // Not found.
1700  // Pout<< "cannot find region " << myRegionAndNormal
1701  // << " in regions " << regionAndNormal << endl;
1702 
1703  return surfaceFeatures::REGION;
1704  }
1705  }
1706 
1707  return surfaceFeatures::NONE;
1708  }
1709 }
1710 
1711 
1714  const triSurface& surf,
1715  const scalar tol,
1716  const scalar includedAngle,
1718 )
1719 {
1720  forAll(edgeStat, edgei)
1721  {
1722  const labelList& eFaces = surf.edgeFaces()[edgei];
1723 
1724  if
1725  (
1726  eFaces.size() > 2
1727  && edgeStat[edgei] == surfaceFeatures::REGION
1728  && (eFaces.size() % 2) == 0
1729  )
1730  {
1731  edgeStat[edgei] = checkNonManifoldEdge
1732  (
1733  surf,
1734  tol,
1735  includedAngle,
1736  edgei
1737  );
1738  }
1739  }
1740 }
1741 
1742 
1743 // ************************************************************************* //
void write(Ostream &, const bool subDict=true) const
Write dictionary, normally with sub-dictionary formatting.
Definition: dictionaryIO.C:207
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
FvWallInfoData< WallInfo, label > label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
A line primitive.
Definition: line.H:56
A class for handling file names.
Definition: fileName.H:79
Holds (reference to) pointField. Encapsulation of data needed for octree searches. Used for searching for nearest point. No bounding boxes around points. Only overlaps and calcNearest are implemented, rest makes little sense.
Definition: treeDataPoint.H:59
label nRegionEdges() const
Return number of region edges.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
static iteratorEnd end()
iteratorEnd set to beyond the end of any HashTable
Definition: HashTable.H:112
void write(const fileName &fName) const
Write as dictionary to file.
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
#define forAllConstIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:477
void selectCutEdges(const triSurface &surf, const plane &cutPlane, List< surfaceFeatures::edgeStatus > &edgeStat)
Select edges that are intersected by the given plane.
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: HashTable.H:59
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Unit conversion functions.
Output to file stream.
Definition: OFstream.H:82
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
bool contains(const point &) const
Contains point? (inside or on edge)
Definition: boundBoxI.H:170
surfaceFeatures::edgeStatus checkNonManifoldEdge(const triSurface &surf, const scalar tol, const scalar includedAngle, const label edgei)
Divide into multiple normal bins.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
label count(const ListType &l, typename ListType::const_reference x)
Count the number of occurrences of a value in a list.
Holds data for octree to work on an edges subset.
Definition: treeDataEdge.H:53
label externalStart() const
Start of external edges.
Map< label > nearestSamples(const labelList &selectedPoints, const pointField &samples, const scalarField &maxDistSqr) const
Find nearest sample for selected surface points.
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
void selectBox(const triSurface &surf, const boundBox &bb, const bool removeInside, List< surfaceFeatures::edgeStatus > &edgeStat)
Select edges inside or outside bounding box.
A bounding box defined in terms of the points at its extremities.
Definition: boundBox.H:58
Map< pointIndexHit > nearestEdges(const labelList &selectedEdges, const edgeList &sampleEdges, const labelList &selectedSampleEdges, const pointField &samplePoints, const scalarField &sampleDist, const scalarField &maxDistSqr, const scalar minSampleDist=0.1) const
Like nearestSamples but now gets nearest point on.
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
Definition: PointIndexHit.H:53
PointIndexHit< point > pointIndexHit
Definition: pointIndexHit.H:42
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of point.
Definition: meshTools.C:203
const labelList & featureEdges() const
Return feature edge list.
label otherVertex(const label a) const
Given one vertex, return the other.
Definition: edgeI.H:140
Geometric class that creates a 2D plane and can return the intersection point between a line and the ...
Definition: plane.H:60
const Field< PointType > & localPoints() const
Return pointField of points in patch.
const fileName & name() const
Return the name of the stream.
Definition: OFstream.H:120
void writeObj(const fileName &prefix) const
Write to separate OBJ files (region, external, internal edges,.
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
const labelList & meshPoints() const
Return labelList of mesh points in patch. They are constructed.
point centre(const pointField &) const
Return centre (centroid)
Definition: edgeI.H:169
bool add(entry *, bool mergeEntry=false)
Add a new entry.
Definition: dictionary.C:1153
labelList selectFeatureEdges(const bool regionEdges, const bool externalEdges, const bool internalEdges) const
Helper function: select a subset of featureEdges_.
int edgeDirection(const edge &) const
Return the edge direction on the face.
Definition: triFaceI.H:352
List< edgeStatus > toStatus() const
From member feature edges to status per edge.
bool insert(const edge &, const T &newElmt)
Insert a new hashedEntry.
Definition: HashTableI.H:80
iterator find(const Key &)
Find and return an iterator set at the hashedEntry.
Definition: HashTable.C:142
const Vector< Cmpt > & centre(const Foam::List< Vector< Cmpt >> &) const
Return *this (used for point which is a typedef to Vector<scalar>.
Definition: VectorI.H:116
void selectManifoldEdges(const triSurface &surf, const scalar tol, const scalar includedAngle, List< surfaceFeatures::edgeStatus > &edgeStat)
Select manifold edges.
const Point & hitPoint() const
Return hit point.
labelList trimFeatures(const scalar minLen, const label minElems, const scalar includedAngle)
Delete small sets of edges. Edges are stringed up and any.
label region() const
Return region label.
Definition: labelledTriI.H:68
void operator=(const surfaceFeatures &)
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
linePointRef line(const pointField &) const
Return edge line.
Definition: edgeI.H:187
dimensionedScalar cos(const dimensionedScalar &ds)
line< point, const point & > linePointRef
Line using referred points.
Definition: linePointRef.H:45
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
Definition: DynamicList.H:56
Mid-point interpolation (weighting factors = 0.5) scheme class.
Definition: midPoint.H:50
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
void clear()
Clear the list, i.e. set size to zero.
Definition: ListI.H:125
void setFromStatus(const List< edgeStatus > &edgeStat, const scalar includedAngle)
Set from status per edge.
void append(const T &)
Append an element at the end of the list.
Definition: ListI.H:178
void setCapacity(const label)
Alter the size of the underlying storage.
Definition: DynamicListI.H:130
const triSurface & surface() const
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:296
bool hit() const
Is there a hit.
const Field< PointType > & points() const
Return reference to global points.
const Type & shapes() const
Reference to shape.
surfaceFeatures(const triSurface &)
Construct from surface.
const labelListList & edgeFaces() const
Return edge-face addressing.
const edgeList & edges() const
Return list of edges, address into LOCAL point list.
List< label > labelList
A List of labels.
Definition: labelList.H:56
scalar mag(const pointField &) const
Return scalar magnitude.
Definition: edgeI.H:181
Triangle with additional region number.
Definition: labelledTri.H:57
const Field< PointType > & faceNormals() const
Return face normals for patch.
const treeBoundBox & bb() const
Top bounding box.
errorManip< error > abort(error &err)
Definition: errorManip.H:131
const Point & rawPoint() const
Return point with no checking.
void inflate(const scalar s)
Inflate box by factor*mag(span) in all dimensions.
Definition: boundBox.C:210
dimensioned< scalar > magSqr(const dimensioned< Type > &)
An Ostream is an abstract base class for all output systems (streams, files, token lists...
Definition: Ostream.H:54
dimensionedScalar sin(const dimensionedScalar &ds)
vector vec(const pointField &) const
Return the vector (end - start)
Definition: edgeI.H:175
label nEdges() const
Return number of edges in patch.
DynamicList< T, SizeInc, SizeMult, SizeDiv > & shrink()
Shrink the allocated space to the number of elements used.
Definition: DynamicListI.H:252
static const char nl
Definition: Ostream.H:260
const labelList & featurePoints() const
Return feature point list.
Map from edge (expressed as its endpoints) to value.
Definition: EdgeMap.H:47
defineTypeNameAndDebug(combustionModel, 0)
Input from file stream.
Definition: IFstream.H:81
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
volScalarField sf(fieldObject, mesh)
messageStream Warning
scalar lineIntersect(const line< Point, PointRef > &l) const
Return the cutting point between the plane and.
Definition: plane.H:186
Point centre() const
Return centre (centroid)
Definition: lineI.H:73
void setSize(const label)
Reset size of List.
Definition: List.C:281
vector point
Point is a vector.
Definition: point.H:41
Non-pointer based hierarchical recursive searching.
Definition: treeDataEdge.H:47
void findFeatures(const scalar includedAngle, const bool geometricTestOnly)
Find feature edges using provided included angle.
label end() const
Return end vertex label.
Definition: edgeI.H:92
void nearestSurfEdge(const labelList &selectedEdges, const pointField &samples, scalar searchSpanSqr, labelList &edgeLabel, labelList &edgeEndPoint, pointField &edgePoint) const
Find nearest surface edge (out of selectedEdges) for.
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
const List< FaceType > & localFaces() const
Return patch faces addressing into local point list.
void writeDict(Ostream &) const
Write as dictionary.
void nearestFeatEdge(const edgeList &edges, const pointField &points, scalar searchSpanSqr, labelList &edgeLabel) const
Find nearest feature edge to each surface edge. Uses the.
Standard boundBox + extra functionality for use in octree.
Definition: treeBoundBox.H:87
dimensioned< scalar > mag(const dimensioned< Type > &)
label n
label nInternalEdges() const
Return number of internal edges.
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:105
PointHit< point > pointHit
Definition: pointHit.H:41
label internalStart() const
Start of internal edges.
T & last()
Return the last element of the list.
Definition: UListI.H:128
Triangulated surface description with patch information.
Definition: triSurface.H:66
label index() const
Return index.
PointHit< Point > nearestDist(const Point &p) const
Return nearest distance to line from a given point.
Definition: lineI.H:95
label start() const
Return start vertex label.
Definition: edgeI.H:81
Holds feature edges/points of surface.
label nExternalEdges() const
Return number of external edges.
Namespace for OpenFOAM.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:864
A HashTable to objects of type <T> with a label key.
Definition: Map.H:49