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