surfaceIntersection.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) 2011-2016 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 "surfaceIntersection.H"
27 #include "triSurfaceSearch.H"
28 #include "labelPairLookup.H"
29 #include "OFstream.H"
30 #include "HashSet.H"
31 #include "triSurface.H"
32 #include "pointIndexHit.H"
33 #include "mergePoints.H"
34 #include "plane.H"
35 #include "edgeIntersections.H"
36 
37 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
38 
39 namespace Foam
40 {
41 defineTypeNameAndDebug(surfaceIntersection, 0);
42 }
43 
44 
45 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
46 
47 // Checks if there exists a special topological situation that causes
48 // edge and the face it hit not to be recognized.
49 //
50 // For now if the face shares a point with the edge
51 bool Foam::surfaceIntersection::excludeEdgeHit
52 (
53  const triSurface& surf,
54  const label edgeI,
55  const label facei,
56  const scalar
57 )
58 {
59  const triSurface::FaceType& f = surf.localFaces()[facei];
60  const edge& e = surf.edges()[edgeI];
61 
62  forAll(f, fp)
63  {
64  if (f[0] == e.start() || f[0] == e.end())
65  {
66  return true;
67  }
68  }
69 
70 // {
71 // // Get edge vector
72 // vector eVec = e.vec(surf.localPoints());
73 // eVec /= mag(eVec) + VSMALL;
74 //
75 // const labelList& eLabels = surf.faceEdges()[facei];
76 //
77 // // Get edge vector of 0th edge of face
78 // vector e0Vec = surf.edges()[eLabels[0]].vec(surf.localPoints());
79 // e0Vec /= mag(e0Vec) + VSMALL;
80 //
81 // vector n = e0Vec ^ eVec;
82 //
83 // if (mag(n) < SMALL)
84 // {
85 // // e0 is aligned with e. Choose next edge of face.
86 // vector e1Vec = surf.edges()[eLabels[1]].vec(surf.localPoints());
87 // e1Vec /= mag(e1Vec) + VSMALL;
88 //
89 // n = e1Vec ^ eVec;
90 //
91 // if (mag(n) < SMALL)
92 // {
93 // // Problematic triangle. Two edges aligned with edgeI. Give
94 // // up.
95 // return true;
96 // }
97 // }
98 //
99 // // Check if same as faceNormal
100 // if (mag(n & surf.faceNormals()[facei]) > 1-tol)
101 // {
102 //
103 // Pout<< "edge:" << e << " face:" << facei
104 // << " e0Vec:" << e0Vec << " n:" << n
105 // << " normalComponent:" << (n & surf.faceNormals()[facei])
106 // << " tol:" << tol << endl;
107 //
108 // return true;
109 // }
110 // else
111 // {
112 // return false;
113 // }
114 // }
115 
116  return false;
117 }
118 
119 
123 //Foam::pointIndexHit Foam::surfaceIntersection::faceEdgeIntersection
124 //(
125 // const triSurface& surf,
126 // const label hitFacei,
127 //
128 // const vector& n,
129 // const point& eStart,
130 // const point& eEnd
131 //)
132 //{
133 // pointIndexHit pInter;
134 //
135 // const pointField& points = surf.points();
136 //
137 // const triSurface::FaceType& f = surf.localFaces()[hitFacei];
138 //
139 // // Plane for intersect test.
140 // plane pl(eStart, n);
141 //
142 // forAll(f, fp)
143 // {
144 // label fp1 = f.fcIndex(fp);
145 //
146 // const point& start = points[f[fp]];
147 // const point& end = points[f[fp1]];
148 //
149 // vector eVec(end - start);
150 //
151 // scalar s = pl.normalIntersect(start, eVec);
152 //
153 // if (s < 0 || s > 1)
154 // {
155 // pInter.setPoint(start + s*eVec);
156 //
157 // // Check if is correct one: orientation walking
158 // // eStart - eEnd - hitPoint should be opposite n
159 // vector n2(triPointRef(start, end, pInter.hitPoint()).normal());
160 //
161 // Pout<< "plane normal:" << n
162 // << " start:" << start << " end:" << end
163 // << " hit at:" << pInter.hitPoint()
164 // << " resulting normal:" << n2 << endl;
165 //
166 // if ((n2 & n) < 0)
167 // {
168 // pInter.setHit();
169 //
170 // // Find corresponding edge between f[fp] f[fp1]
171 // label edgeI =
172 // meshTools::findEdge
173 // (
174 // surf.edges(),
175 // surf.faceEdges()[hitFacei],
176 // f[fp],
177 // f[fp1]
178 // );
179 //
180 // pInter.setIndex(edgeI);
181 //
182 // return pInter;
183 // }
184 // }
185 // }
186 //
187 // FatalErrorInFunction
188 // << "Did not find intersection of plane " << pl
189 // << " with edges of face " << hitFacei << " verts:" << f
190 // << abort(FatalError);
191 //
192 // return pInter;
193 //}
194 
195 
196 void Foam::surfaceIntersection::storeIntersection
197 (
198  const bool isFirstSurf,
199  const labelList& facesA,
200  const label faceB,
201  DynamicList<edge>& allCutEdges,
202  DynamicList<point>& allCutPoints
203 )
204 {
205  forAll(facesA, facesAI)
206  {
207  label faceA = facesA[facesAI];
208 
209  // Combine two faces. Always make sure the face from the first surface
210  // is element 0.
211  FixedList<label, 2> twoFaces;
212  if (isFirstSurf)
213  {
214  twoFaces[0] = faceA;
215  twoFaces[1] = faceB;
216  }
217  else
218  {
219  twoFaces[0] = faceB;
220  twoFaces[1] = faceA;
221  }
222 
223  labelPairLookup::const_iterator iter = facePairToVertex_.find(twoFaces);
224 
225  if (iter == facePairToVertex_.end())
226  {
227  // New intersection. Store face-face intersection.
228  facePairToVertex_.insert(twoFaces, allCutPoints.size()-1);
229  }
230  else
231  {
232  // Second occurrence of surf1-surf2 intersection.
233  // Or rather the face on surf1 intersects a face on
234  // surface2 twice -> we found edge.
235 
236  // Check whether perhaps degenerate
237  const point& prevHit = allCutPoints[*iter];
238  const point& thisHit = allCutPoints.last();
239 
240  if (mag(prevHit - thisHit) < SMALL)
241  {
243  << "Encountered degenerate edge between face "
244  << twoFaces[0] << " on first surface"
245  << " and face " << twoFaces[1] << " on second surface"
246  << endl
247  << "Point on first surface:" << prevHit << endl
248  << "Point on second surface:" << thisHit << endl
249  << endl;
250  }
251  else
252  {
253  allCutEdges.append(edge(*iter, allCutPoints.size()-1));
254 
255  // Remember face on surf
256  facePairToEdge_.insert(twoFaces, allCutEdges.size()-1);
257  }
258  }
259  }
260 }
261 
262 
263 // Classify cut of edge of surface1 with surface2:
264 // 1- point of edge hits point on surface2
265 // 2- edge pierces point on surface2
266 // 3- point of edge hits edge on surface2
267 // 4- edge pierces edge on surface2
268 // 5- point of edge hits face on surface2
269 // 6- edge pierces face on surface2
270 //
271 // Note that handling of 2 and 4 should be the same but with surface1 and
272 // surface2 reversed.
273 void Foam::surfaceIntersection::classifyHit
274 (
275  const triSurface& surf1,
276  const scalarField& surf1PointTol,
277  const triSurface& surf2,
278  const bool isFirstSurf,
279  const label edgeI,
280  const scalar tolDim,
281  const pointIndexHit& pHit,
282 
283  DynamicList<edge>& allCutEdges,
284  DynamicList<point>& allCutPoints,
285  List<DynamicList<label>>& surfEdgeCuts
286 )
287 {
288  const edge& e = surf1.edges()[edgeI];
289 
290  const labelList& facesA = surf1.edgeFaces()[edgeI];
291 
292  // Label of face on surface2 edgeI intersected
293  label surf2Facei = pHit.index();
294 
295  // Classify point on surface2
296 
297  const triSurface::FaceType& f2 = surf2.localFaces()[surf2Facei];
298  const pointField& surf2Pts = surf2.localPoints();
299 
300  label nearType, nearLabel;
301 
302  f2.nearestPointClassify(pHit.hitPoint(), surf2Pts, nearType, nearLabel);
303 
304  // Classify points on edge of surface1
305  label edgeEnd =
306  classify
307  (
308  surf1PointTol[e.start()],
309  surf1PointTol[e.end()],
310  pHit.hitPoint(),
311  e,
312  surf1.localPoints()
313  );
314 
315  if (nearType == triPointRef::POINT)
316  {
317  if (edgeEnd >= 0)
318  {
319  // 1. Point hits point. Do nothing.
320  if (debug & 2)
321  {
322  Pout<< pHit.hitPoint() << " is surf1:"
323  << " end point of edge " << e
324  << " surf2: vertex " << f2[nearLabel]
325  << " coord:" << surf2Pts[f2[nearLabel]] << endl;
326  }
327  }
328  else
329  {
330  // 2. Edge hits point. Cut edge with new point.
331  if (debug & 2)
332  {
333  Pout<< pHit.hitPoint() << " is surf1:"
334  << " somewhere on edge " << e
335  << " surf2: vertex " << f2[nearLabel]
336  << " coord:" << surf2Pts[f2[nearLabel]] << endl;
337  }
338 
339  allCutPoints.append(pHit.hitPoint());
340  surfEdgeCuts[edgeI].append(allCutPoints.size()-1);
341 
342  const labelList& facesB = surf2.pointFaces()[f2[nearLabel]];
343 
344  forAll(facesB, faceBI)
345  {
346  storeIntersection
347  (
348  isFirstSurf,
349  facesA,
350  facesB[faceBI],
351  allCutEdges,
352  allCutPoints
353  );
354  }
355  }
356  }
357  else if (nearType == triPointRef::EDGE)
358  {
359  if (edgeEnd >= 0)
360  {
361  // 3. Point hits edge. Do nothing on this side. Reverse
362  // is handled by 2 (edge hits point)
363  label edge2I = getEdge(surf2, surf2Facei, nearLabel);
364  const edge& e2 = surf2.edges()[edge2I];
365 
366  if (debug&2)
367  {
368  Pout<< pHit.hitPoint() << " is surf1:"
369  << " end point of edge " << e
370  << " surf2: edge " << e2
371  << " coords:" << surf2Pts[e2.start()]
372  << surf2Pts[e2.end()] << endl;
373  }
374  }
375  else
376  {
377  // 4. Edge hits edge.
378 
379  // Cut edge with new point (creates duplicates when
380  // doing the surf2 with surf1 intersection but these
381  // are merged later on)
382 
383  label edge2I = getEdge(surf2, surf2Facei, nearLabel);
384  const edge& e2 = surf2.edges()[edge2I];
385 
386  if (debug&2)
387  {
388  Pout<< pHit.hitPoint() << " is surf1:"
389  << " somewhere on edge " << e
390  << " surf2: edge " << e2
391  << " coords:" << surf2Pts[e2.start()]
392  << surf2Pts[e2.end()] << endl;
393  }
394 
395  allCutPoints.append(pHit.hitPoint());
396  surfEdgeCuts[edgeI].append(allCutPoints.size()-1);
397 
398  // edge hits all faces on surf2 connected to the edge
399 
400  if (isFirstSurf)
401  {
402  // edge-edge intersection is symmetric, store only
403  // once.
404  // edge hits all faces on surf2 connected to the
405  // edge
406 
407  const labelList& facesB = surf2.edgeFaces()[edge2I];
408 
409  forAll(facesB, faceBI)
410  {
411  storeIntersection
412  (
413  isFirstSurf,
414  facesA,
415  facesB[faceBI],
416  allCutEdges,
417  allCutPoints
418  );
419  }
420  }
421  }
422  }
423  else
424  {
425  if (edgeEnd >= 0)
426  {
427  // 5. Point hits face. Do what? Introduce
428  // point & triangulation in face?
429  if (debug&2)
430  {
431  Pout<< pHit.hitPoint() << " is surf1:"
432  << " end point of edge " << e
433  << " surf2: face " << surf2Facei
434  << endl;
435  }
436 
437  //
438  // Look exactly at what side (of surf2) edge is. Leave out ones on
439  // inside of surf2 (i.e. on opposite side of normal)
440  //
441 
442  // Vertex on/near surf2
443  label nearVert = -1;
444 
445  if (edgeEnd == 0)
446  {
447  nearVert = e.start();
448  }
449  else
450  {
451  nearVert = e.end();
452  }
453 
454  const point& nearPt = surf1.localPoints()[nearVert];
455 
456  // Vertex away from surf2
457  label otherVert = e.otherVertex(nearVert);
458 
459  const point& otherPt = surf1.localPoints()[otherVert];
460 
461 
462  if (debug)
463  {
464  Pout
465  << pHit.hitPoint() << " is surf1:"
466  << " end point of edge " << e << " coord:"
467  << surf1.localPoints()[nearVert]
468  << " surf2: face " << surf2Facei << endl;
469  }
470 
471  vector eVec = otherPt - nearPt;
472 
473  if ((surf2.faceNormals()[surf2Facei] & eVec) > 0)
474  {
475  // otherVert on outside of surf2
476 
477  // Shift hitPoint a bit along edge.
478  //point hitPt = nearPt + 0.1*eVec;
479  point hitPt = nearPt;
480 
481  if (debug&2)
482  {
483  Pout<< "Shifted " << pHit.hitPoint()
484  << " to " << hitPt
485  << " along edge:" << e
486  << " coords:" << surf1.localPoints()[e.start()]
487  << surf1.localPoints()[e.end()] << endl;
488  }
489 
490  // Reclassify as normal edge-face pierce (see below)
491 
492  allCutPoints.append(hitPt);
493  surfEdgeCuts[edgeI].append(allCutPoints.size()-1);
494 
495  // edge hits single face only
496  storeIntersection
497  (
498  isFirstSurf,
499  facesA,
500  surf2Facei,
501  allCutEdges,
502  allCutPoints
503  );
504  }
505  else
506  {
507  if (debug&2)
508  {
509  Pout<< "Discarding " << pHit.hitPoint()
510  << " since edge " << e << " on inside of surf2."
511  << " surf2 normal:" << surf2.faceNormals()[surf2Facei]
512  << endl;
513  }
514  }
515  }
516  else
517  {
518  // 6. Edge pierces face. 'Normal' situation.
519  if (debug&2)
520  {
521  Pout<< pHit.hitPoint() << " is surf1:"
522  << " somewhere on edge " << e
523  << " surf2: face " << surf2Facei
524  << endl;
525  }
526 
527  // edgeI intersects surf2. Store point.
528  allCutPoints.append(pHit.hitPoint());
529  surfEdgeCuts[edgeI].append(allCutPoints.size()-1);
530 
531  // edge hits single face only
532  storeIntersection
533  (
534  isFirstSurf,
535  facesA,
536  surf2Facei,
537  allCutEdges,
538  allCutPoints
539  );
540  }
541  }
542  if (debug&2)
543  {
544  Pout<< endl;
545  }
546 }
547 
548 
549 // Cut all edges of surf1 with surf2. Sets
550 // - cutPoints : coordinates of cutPoints
551 // - cutEdges : newly created edges between cutPoints
552 // - facePairToVertex : hash from face1I and face2I to cutPoint
553 // - facePairToEdge : hash from face1I and face2I to cutEdge
554 // - surfEdgeCuts : gives for each edge the cutPoints
555 // (in order from start to end)
556 //
557 void Foam::surfaceIntersection::doCutEdges
558 (
559  const triSurface& surf1,
560  const triSurfaceSearch& querySurf2,
561  const bool isFirstSurf,
562  const bool isSelfIntersection,
563 
564  DynamicList<edge>& allCutEdges,
565  DynamicList<point>& allCutPoints,
566  List<DynamicList<label>>& surfEdgeCuts
567 )
568 {
569  scalar oldTol = intersection::setPlanarTol(1e-3);
570 
571  const pointField& surf1Pts = surf1.localPoints();
572 
573  // Calculate local (to point) tolerance based on min edge length.
574  scalarField surf1PointTol(surf1Pts.size());
575 
576  forAll(surf1PointTol, pointi)
577  {
578  surf1PointTol[pointi] =
580  * minEdgeLen(surf1, pointi);
581  }
582 
583  const triSurface& surf2 = querySurf2.surface();
584 
585  forAll(surf1.edges(), edgeI)
586  {
587  const edge& e = surf1.edges()[edgeI];
588 
589  point pStart = surf1Pts[e.start()];
590  const point& pEnd = surf1Pts[e.end()];
591 
592  const point tolVec = intersection::planarTol()*(pEnd-pStart);
593  const scalar tolDim = mag(tolVec);
594 
595  bool doTrack = false;
596  do
597  {
598  pointIndexHit pHit = querySurf2.tree().findLine(pStart, pEnd);
599 
600  if (pHit.hit())
601  {
602  if (isSelfIntersection)
603  {
604  // Skip all intersections which are hit at endpoints of
605  // edge.
606  // Problem is that if faces are almost coincident the
607  // intersection point will be calculated quite incorrectly
608  // The error might easily be larger than 1% of the edge
609  // length.
610  // So what we do here is to exclude hit faces if our edge
611  // is in their plane and they share a point with the edge.
612 
613  // Label of face on surface2 edgeI intersected
614  label hitFacei = pHit.index();
615 
616  if
617  (
618  !excludeEdgeHit
619  (
620  surf1,
621  edgeI,
622  hitFacei,
623  0.1 // 1-cos of angle between normals
624  )
625  )
626  {
627  // Classify point on surface1
628  label edgeEnd = classify
629  (
630  surf1PointTol[e.start()],
631  surf1PointTol[e.end()],
632  pHit.hitPoint(),
633  e,
634  surf1Pts
635  );
636 
637  if (edgeEnd < 0)
638  {
639  if (debug)
640  {
641  Pout<< "edge:" << edgeI << " vertices:" << e
642  << " start:" << surf1Pts[e.start()]
643  << " end:" << surf1Pts[e.end()]
644  << " hit:" << pHit.hitPoint()
645  << " tolDim:" << tolDim
646  << " planarTol:"
648  << endl;
649  }
650  allCutPoints.append(pHit.hitPoint());
651  surfEdgeCuts[edgeI].append(allCutPoints.size()-1);
652  }
653  }
654  }
655  else
656  {
657  classifyHit
658  (
659  surf1,
660  surf1PointTol,
661  surf2,
662  isFirstSurf,
663  edgeI,
664  tolDim,
665  pHit,
666 
667  allCutEdges,
668  allCutPoints,
669  surfEdgeCuts
670  );
671  }
672 
673  if (mag(pHit.hitPoint() - pEnd) < tolDim)
674  {
675  doTrack = false;
676  }
677  else
678  {
679  pStart = pHit.hitPoint() + tolVec;
680 
681  doTrack = true;
682  }
683  }
684  else
685  {
686  doTrack = false;
687  }
688  }
689  while (doTrack);
690  }
691 
693 }
694 
695 
696 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
697 
698 // Null constructor
700 :
701  cutPoints_(0),
702  cutEdges_(0),
703  facePairToVertex_(0),
704  facePairToEdge_(0),
705  surf1EdgeCuts_(0),
706  surf2EdgeCuts_(0)
707 {}
708 
709 
710 // Construct from two surfaces
712 (
713  const triSurfaceSearch& query1,
714  const triSurfaceSearch& query2
715 )
716 :
717  cutPoints_(0),
718  cutEdges_(0),
719  facePairToVertex_(2*max(query1.surface().size(), query2.surface().size())),
720  facePairToEdge_(2*max(query1.surface().size(), query2.surface().size())),
721  surf1EdgeCuts_(0),
722  surf2EdgeCuts_(0)
723 {
724  const triSurface& surf1 = query1.surface();
725  const triSurface& surf2 = query2.surface();
726 
727  //
728  // Cut all edges of surf1 with surf2.
729  //
730  if (debug)
731  {
732  Pout<< "Cutting surf1 edges" << endl;
733  }
734 
735 
736  DynamicList<edge> allCutEdges(surf1.nEdges()/20);
737  DynamicList<point> allCutPoints(surf1.nPoints()/20);
738 
739 
740  // From edge to cut index on surface1
741  List<DynamicList<label>> edgeCuts1(query1.surface().nEdges());
742 
743  doCutEdges
744  (
745  surf1,
746  query2,
747  true, // is first surface; construct labelPair in correct
748  // order
749  false, // not self intersection
750 
751  allCutEdges,
752  allCutPoints,
753  edgeCuts1
754  );
755  // Transfer to straight labelListList
756  transfer(edgeCuts1, surf1EdgeCuts_);
757 
758 
759  //
760  // Cut all edges of surf2 with surf1.
761  //
762  if (debug)
763  {
764  Pout<< "Cutting surf2 edges" << endl;
765  }
766 
767  // From edge to cut index
768  List<DynamicList<label>> edgeCuts2(query2.surface().nEdges());
769 
770  doCutEdges
771  (
772  surf2,
773  query1,
774  false, // is second surface
775  false, // not self intersection
776 
777  allCutEdges,
778  allCutPoints,
779  edgeCuts2
780  );
781 
782  // Transfer to straight label(List)List
783  transfer(edgeCuts2, surf2EdgeCuts_);
784  cutEdges_.transfer(allCutEdges);
785  cutPoints_.transfer(allCutPoints);
786 
787 
788  if (debug)
789  {
790  Pout<< "surfaceIntersection : Intersection generated:"
791  << endl
792  << " points:" << cutPoints_.size() << endl
793  << " edges :" << cutEdges_.size() << endl;
794 
795  Pout<< "surfaceIntersection : Writing intersection to intEdges.obj"
796  << endl;
797 
798  OFstream intStream("intEdges.obj");
799  writeOBJ(cutPoints_, cutEdges_, intStream);
800 
801  // Dump all cut edges to files
802  Pout<< "Dumping cut edges of surface1 to surf1EdgeCuts.obj" << endl;
803  OFstream edge1Stream("surf1EdgeCuts.obj");
804  writeIntersectedEdges(surf1, surf1EdgeCuts_, edge1Stream);
805 
806  Pout<< "Dumping cut edges of surface2 to surf2EdgeCuts.obj" << endl;
807  OFstream edge2Stream("surf2EdgeCuts.obj");
808  writeIntersectedEdges(surf2, surf2EdgeCuts_, edge2Stream);
809  }
810 }
811 
812 
813 // Construct from full intersection Poutrmation
815 (
816  const triSurface& surf1,
817  const edgeIntersections& intersections1,
818  const triSurface& surf2,
819  const edgeIntersections& intersections2
820 )
821 :
822  cutPoints_(0),
823  cutEdges_(0),
824  facePairToVertex_(2*max(surf1.size(), surf2.size())),
825  facePairToEdge_(2*max(surf1.size(), surf2.size())),
826  surf1EdgeCuts_(0),
827  surf2EdgeCuts_(0)
828 {
829 
830  // All intersection Pout (so for both surfaces)
831  DynamicList<edge> allCutEdges((surf1.nEdges() + surf2.nEdges())/20);
832  DynamicList<point> allCutPoints((surf1.nPoints() + surf2.nPoints())/20);
833 
834 
835  // Cut all edges of surf1 with surf2
836  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
837 
838  if (debug)
839  {
840  Pout<< "Storing surf1 intersections" << endl;
841  }
842 
843  {
844  // From edge to cut index on surface1
845  List<DynamicList<label>> edgeCuts1(surf1.nEdges());
846 
847  forAll(intersections1, edgeI)
848  {
849  const List<pointIndexHit>& intersections = intersections1[edgeI];
850 
851  forAll(intersections, i)
852  {
853  const pointIndexHit& pHit = intersections[i];
854 
855  // edgeI intersects surf2. Store point.
856  allCutPoints.append(pHit.hitPoint());
857  edgeCuts1[edgeI].append(allCutPoints.size()-1);
858 
859  storeIntersection
860  (
861  true, // is first surface
862  surf1.edgeFaces()[edgeI],
863  pHit.index(), // surf2Facei
864  allCutEdges,
865  allCutPoints
866  );
867  }
868  }
869 
870  // Transfer to straight labelListList
871  transfer(edgeCuts1, surf1EdgeCuts_);
872  }
873 
874 
875 
876  // Cut all edges of surf2 with surf1
877  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
878 
879  if (debug)
880  {
881  Pout<< "Storing surf2 intersections" << endl;
882  }
883 
884  {
885  // From edge to cut index on surface2
886  List<DynamicList<label>> edgeCuts2(surf2.nEdges());
887 
888  forAll(intersections2, edgeI)
889  {
890  const List<pointIndexHit>& intersections = intersections2[edgeI];
891 
892  forAll(intersections, i)
893  {
894  const pointIndexHit& pHit = intersections[i];
895 
896  // edgeI intersects surf1. Store point.
897  allCutPoints.append(pHit.hitPoint());
898  edgeCuts2[edgeI].append(allCutPoints.size()-1);
899 
900  storeIntersection
901  (
902  false, // is second surface
903  surf2.edgeFaces()[edgeI],
904  pHit.index(), // surf2Facei
905  allCutEdges,
906  allCutPoints
907  );
908  }
909  }
910 
911  // Transfer to surf2EdgeCuts_ (straight labelListList)
912  transfer(edgeCuts2, surf2EdgeCuts_);
913  }
914 
915 
916  // Transfer to straight label(List)List
917  cutEdges_.transfer(allCutEdges);
918  cutPoints_.transfer(allCutPoints);
919 
920 
921  if (debug)
922  {
923  Pout<< "surfaceIntersection : Intersection generated:"
924  << endl
925  << " points:" << cutPoints_.size() << endl
926  << " edges :" << cutEdges_.size() << endl;
927 
928  Pout<< "surfaceIntersection : Writing intersection to intEdges.obj"
929  << endl;
930 
931  OFstream intStream("intEdges.obj");
932  writeOBJ(cutPoints_, cutEdges_, intStream);
933 
934  // Dump all cut edges to files
935  Pout<< "Dumping cut edges of surface1 to surf1EdgeCuts.obj" << endl;
936  OFstream edge1Stream("surf1EdgeCuts.obj");
937  writeIntersectedEdges(surf1, surf1EdgeCuts_, edge1Stream);
938 
939  Pout<< "Dumping cut edges of surface2 to surf2EdgeCuts.obj" << endl;
940  OFstream edge2Stream("surf2EdgeCuts.obj");
941  writeIntersectedEdges(surf2, surf2EdgeCuts_, edge2Stream);
942  }
943 
944 
945  // Debugging stuff
946  {
947  // Check all facePairToVertex is used.
948  labelHashSet usedPoints;
949 
950  forAllConstIter(labelPairLookup, facePairToEdge_, iter)
951  {
952  label edgeI = iter();
953 
954  const edge& e = cutEdges_[edgeI];
955 
956  usedPoints.insert(e[0]);
957  usedPoints.insert(e[1]);
958  }
959 
960  forAllConstIter(labelPairLookup, facePairToVertex_, iter)
961  {
962  label pointi = iter();
963 
964  if (!usedPoints.found(pointi))
965  {
967  << "Problem: cut point:" << pointi
968  << " coord:" << cutPoints_[pointi]
969  << " not used by any edge" << endl;
970  }
971  }
972  }
973 }
974 
975 
976 // Construct from single surface. Used to test for self-intersection.
978 (
979  const triSurfaceSearch& query1
980 )
981 :
982  cutPoints_(0),
983  cutEdges_(0),
984  facePairToVertex_(2*query1.surface().size()),
985  facePairToEdge_(2*query1.surface().size()),
986  surf1EdgeCuts_(0),
987  surf2EdgeCuts_(0)
988 {
989  const triSurface& surf1 = query1.surface();
990 
991  //
992  // Cut all edges of surf1 with surf1 itself.
993  //
994  if (debug)
995  {
996  Pout<< "Cutting surf1 edges" << endl;
997  }
998 
999  DynamicList<edge> allCutEdges;
1000  DynamicList<point> allCutPoints;
1001 
1002  // From edge to cut index on surface1
1003  List<DynamicList<label>> edgeCuts1(query1.surface().nEdges());
1004 
1005  doCutEdges
1006  (
1007  surf1,
1008  query1,
1009  true, // is first surface; construct labelPair in correct
1010  // order
1011  true, // self intersection
1012 
1013 
1014  allCutEdges,
1015  allCutPoints,
1016  edgeCuts1
1017  );
1018 
1019  // Transfer to straight label(List)List
1020  transfer(edgeCuts1, surf1EdgeCuts_);
1021  cutEdges_.transfer(allCutEdges);
1022  cutPoints_.transfer(allCutPoints);
1023 
1024  // Shortcut.
1025  if (cutPoints_.empty() && cutEdges_.empty())
1026  {
1027  if (debug)
1028  {
1029  Pout<< "Empty intersection" << endl;
1030  }
1031  return;
1032  }
1033 
1034  //
1035  // Remove duplicate points (from edge-point or edge-edge cutting)
1036  //
1037 
1038  // Get typical dimension.
1039  scalar minEdgeLen = GREAT;
1040  forAll(surf1.edges(), edgeI)
1041  {
1042  minEdgeLen = min
1043  (
1044  minEdgeLen,
1045  surf1.edges()[edgeI].mag(surf1.localPoints())
1046  );
1047  }
1048 
1049  // Merge points
1050  labelList pointMap;
1051  pointField newPoints;
1052 
1053  bool hasMerged = mergePoints
1054  (
1055  cutPoints_,
1056  minEdgeLen*intersection::planarTol(),
1057  false,
1058  pointMap,
1059  newPoints
1060  );
1061 
1062  if (hasMerged)
1063  {
1064  if (debug)
1065  {
1066  Pout<< "Merged:" << hasMerged
1067  << " mergeDist:" << minEdgeLen*intersection::planarTol()
1068  << " cutPoints:" << cutPoints_.size()
1069  << " newPoints:" << newPoints.size()
1070  << endl;
1071  }
1072 
1073  // Copy points
1074  cutPoints_.transfer(newPoints);
1075 
1076  // Renumber vertices referenced by edges
1077  forAll(cutEdges_, edgeI)
1078  {
1079  edge& e = cutEdges_[edgeI];
1080 
1081  e.start() = pointMap[e.start()];
1082  e.end() = pointMap[e.end()];
1083 
1084  if (e.mag(cutPoints_) < minEdgeLen*intersection::planarTol())
1085  {
1086  if (debug)
1087  {
1088  Pout<< "Degenerate cut:" << edgeI << " vertices:" << e
1089  << " coords:" << cutPoints_[e.start()] << ' '
1090  << cutPoints_[e.end()] << endl;
1091  }
1092  }
1093  }
1094 
1095  // Renumber vertices referenced by edgeCut lists. Remove duplicates.
1096  forAll(surf1EdgeCuts_, edgeI)
1097  {
1098  // Get indices of cutPoints this edge is cut by
1099  labelList& cutVerts = surf1EdgeCuts_[edgeI];
1100 
1101  removeDuplicates(pointMap, cutVerts);
1102  }
1103  }
1104 
1105  if (debug)
1106  {
1107  Pout<< "surfaceIntersection : Intersection generated and compressed:"
1108  << endl
1109  << " points:" << cutPoints_.size() << endl
1110  << " edges :" << cutEdges_.size() << endl;
1111 
1112 
1113  Pout<< "surfaceIntersection : Writing intersection to intEdges.obj"
1114  << endl;
1115 
1116  OFstream intStream("intEdges.obj");
1117  writeOBJ(cutPoints_, cutEdges_, intStream);
1118  }
1119 
1120  if (debug)
1121  {
1122  // Dump all cut edges to files
1123  Pout<< "Dumping cut edges of surface1 to surf1EdgeCuts.obj" << endl;
1124  OFstream edge1Stream("surf1EdgeCuts.obj");
1125  writeIntersectedEdges(surf1, surf1EdgeCuts_, edge1Stream);
1126  }
1127 }
1128 
1129 
1130 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
1131 
1133 {
1134  return cutPoints_;
1135 }
1136 
1137 
1139 {
1140  return cutEdges_;
1141 }
1142 
1143 
1145 {
1146  return facePairToVertex_;
1147 }
1148 
1149 
1151 {
1152  return facePairToEdge_;
1153 }
1154 
1155 
1158  const bool isFirstSurf
1159 ) const
1160 {
1161  if (isFirstSurf)
1162  {
1163  return surf1EdgeCuts_;
1164  }
1165  else
1166  {
1167  return surf2EdgeCuts_;
1168  }
1169 }
1170 
1171 
1173 {
1174  return surf1EdgeCuts_;
1175 }
1176 
1177 
1179 {
1180  return surf2EdgeCuts_;
1181 }
1182 
1183 
1184 // ************************************************************************* //
label end() const
Return end vertex label.
Definition: edgeI.H:92
surfaceIntersection()
Construct null.
label nPoints() const
Return number of points supporting patch faces.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
const double e
Elementary charge.
Definition: doubleFloat.H:78
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
const triSurface & surface() const
Return reference to the surface.
bool empty() const
Return true if the UList is empty (ie, size() is zero)
Definition: UListI.H:313
const labelListList & surf2EdgeCuts() const
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
Output to file stream.
Definition: OFstream.H:81
const labelListList & edgeCuts(const bool) const
Access either surf1EdgeCuts (isFirstSurface = true) or.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
const labelPairLookup & facePairToEdge() const
const pointField & cutPoints() const
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
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
bool insert(const Key &key)
Insert a new entry.
Definition: HashSet.H:116
const labelPairLookup & facePairToVertex() const
static scalar setPlanarTol(const scalar t)
Set the planar tolerance, returning the previous value.
Definition: intersection.H:91
const labelListList & surf1EdgeCuts() const
Helper class to search on triSurface.
static scalar planarTol()
Return planar tolerance.
Definition: intersection.H:85
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
Definition: DynamicList.H:56
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
friend class const_iterator
Declare friendship with the const_iterator.
Definition: HashTable.H:192
edgeList edges() const
Return edges in face point ordering,.
Definition: triFaceI.H:318
const edgeList & edges() const
Return list of edges, address into LOCAL point list.
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:292
List< label > labelList
A List of labels.
Definition: labelList.H:56
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:29
bool found(const Key &) const
Return true if hashedEntry is found in table.
Definition: HashTable.C:109
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53
label start() const
Return start vertex label.
Definition: edgeI.H:81
scalar mag(const pointField &) const
Return scalar magnitude.
Definition: edgeI.H:163
defineTypeNameAndDebug(combustionModel, 0)
Merge points. See below.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
label nEdges() const
Return number of edges in patch.
vector point
Point is a vector.
Definition: point.H:41
#define WarningInFunction
Report a warning using Foam::Warning.
label index() const
Return index.
label mergePoints(const UList< Type > &points, const scalar mergeTol, const bool verbose, labelList &pointMap, const Type &origin=Type::zero)
Sorts and merges points. All points closer than/equal mergeTol get merged.
const Point & hitPoint() const
Return hit point.
dimensioned< scalar > mag(const dimensioned< Type > &)
const labelListList & edgeFaces() const
Return edge-face addressing.
const edgeList & cutEdges() const
Triangulated surface description with patch information.
Definition: triSurface.H:65
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
Definition: List.C:365
Holder of intersections of edges of a surface with another surface. Optionally shuffles around points...
const Field< PointType > & localPoints() const
Return pointField of points in patch.
Namespace for OpenFOAM.