FreeFOAM The Cross-Platform CFD Toolkit
edgeIntersections.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 "edgeIntersections.H"
31 #include <OpenFOAM/OFstream.H>
32 #include <OpenFOAM/HashSet.H>
33 #include <triSurface/triSurface.H>
37 #include <meshTools/meshTools.H>
38 #include <OpenFOAM/plane.H>
39 #include <OpenFOAM/Random.H>
41 #include <meshTools/treeBoundBox.H>
42 
43 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
44 
46 
49 
50 
51 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
52 
53 void Foam::edgeIntersections::checkEdges(const triSurface& surf)
54 {
55  const pointField& localPoints = surf.localPoints();
56  const edgeList& edges = surf.edges();
57  const labelListList& edgeFaces = surf.edgeFaces();
58 
59  treeBoundBox bb(localPoints);
60 
61  scalar minSize = SMALL * bb.minDim();
62 
63  forAll(edges, edgeI)
64  {
65  const edge& e = edges[edgeI];
66 
67  scalar eMag = e.mag(localPoints);
68 
69  if (eMag < minSize)
70  {
71  WarningIn
72  (
73  "Foam::edgeIntersections::checkEdges(const triSurface& surf)"
74  ) << "Edge " << edgeI << " vertices " << e
75  << " coords:" << localPoints[e[0]] << ' '
76  << localPoints[e[1]] << " is very small compared to bounding"
77  << " box dimensions " << bb << endl
78  << "This might lead to problems in intersection"
79  << endl;
80  }
81 
82  if (edgeFaces[edgeI].size() == 1)
83  {
84  WarningIn
85  (
86  "Foam::edgeIntersections::checkEdges(const triSurface& surf)"
87  ) << "Edge " << edgeI << " vertices " << e
88  << " coords:" << localPoints[e[0]] << ' '
89  << localPoints[e[1]] << " has only one face connected to it:"
90  << edgeFaces[edgeI] << endl
91  << "This might lead to problems in intersection"
92  << endl;
93  }
94  }
95 }
96 
97 
98 // Update intersections for selected edges.
99 void Foam::edgeIntersections::intersectEdges
100 (
101  const triSurface& surf1,
102  const pointField& points1, // surf1 meshPoints (not localPoints!)
103  const triSurfaceSearch& querySurf2,
104  const scalarField& surf1PointTol, // surf1 tolerance per point
105  const labelList& edgeLabels
106 )
107 {
108  const triSurface& surf2 = querySurf2.surface();
109  const vectorField& normals2 = surf2.faceNormals();
110 
111  const labelList& meshPoints = surf1.meshPoints();
112 
113  if (debug)
114  {
115  Pout<< "Calculating intersection of " << edgeLabels.size() << " edges"
116  << " out of " << surf1.nEdges() << " with " << surf2.size()
117  << " triangles ..." << endl;
118  }
119 
120  // Construct octree.
121  const indexedOctree<treeDataTriSurface>& tree = querySurf2.tree();
122 
123 
124  label nHits = 0;
125 
126 
127  // Go through all edges, calculate intersections
128  forAll(edgeLabels, i)
129  {
130  label edgeI = edgeLabels[i];
131 
132  if (debug && (i % 1000 == 0))
133  {
134  Pout<< "Intersecting edge " << edgeI << " with surface" << endl;
135  }
136 
137  const edge& e = surf1.edges()[edgeI];
138 
139  const point& pStart = points1[meshPoints[e.start()]];
140  const point& pEnd = points1[meshPoints[e.end()]];
141 
142  const vector eVec(pEnd - pStart);
143  const scalar eMag = mag(eVec);
144  const vector n(eVec/(eMag + VSMALL));
145 
146  // Smallish length for intersection calculation errors.
147  const point tolVec = 1e-6*eVec;
148 
149  // Start tracking somewhat before pStart and upto somewhat after p1.
150  // Note that tolerances here are smaller than those used to classify
151  // hit below.
152  // This will cause this hit to be marked as degenerate and resolved
153  // later on.
154  point p0 = pStart - 0.5*surf1PointTol[e[0]]*n;
155  const point p1 = pEnd + 0.5*surf1PointTol[e[1]]*n;
156  const scalar maxS = mag(p1 - pStart);
157 
158  // Get all intersections of the edge with the surface
159 
160  DynamicList<pointIndexHit> currentIntersections(100);
161  DynamicList<label> currentIntersectionTypes(100);
162 
163  while (true)
164  {
165  pointIndexHit pHit = tree.findLine(p0, p1);
166 
167  if (pHit.hit())
168  {
169  nHits++;
170 
171  currentIntersections.append(pHit);
172 
173  // Classify point on surface1 edge.
174  label edgeEnd = -1;
175 
176  if (mag(pHit.hitPoint() - pStart) < surf1PointTol[e[0]])
177  {
178  edgeEnd = 0;
179  }
180  else if (mag(pHit.hitPoint() - pEnd) < surf1PointTol[e[1]])
181  {
182  edgeEnd = 1;
183  }
184  else if (mag(n & normals2[pHit.index()]) < alignedCos_)
185  {
186  Pout<< "Flat angle edge:" << edgeI
187  << " face:" << pHit.index()
188  << " cos:" << mag(n & normals2[pHit.index()])
189  << endl;
190  edgeEnd = 2;
191  }
192 
193  currentIntersectionTypes.append(edgeEnd);
194 
195  if (edgeEnd == 1)
196  {
197  // Close to end
198  break;
199  }
200  else
201  {
202  // Continue tracking. Shift by a small amount.
203  p0 = pHit.hitPoint() + tolVec;
204 
205  if (((p0-pStart) & n) >= maxS)
206  {
207  break;
208  }
209  }
210  }
211  else
212  {
213  // No hit.
214  break;
215  }
216  }
217 
218 
219  // Done current edge. Transfer all data into *this
220  operator[](edgeI).transfer(currentIntersections);
221  classification_[edgeI].transfer(currentIntersectionTypes);
222  }
223 
224  if (debug)
225  {
226  Pout<< "Found " << nHits << " intersections of edges with surface ..."
227  << endl;
228  }
229 
230 }
231 
232 
233 // If edgeI intersections are close to endpoint of edge shift endpoints
234 // slightly along edge
235 // Updates
236 // - points1 with new endpoint position
237 // - affectedEdges with all edges affected by moving the point
238 // Returns true if changed anything.
239 bool Foam::edgeIntersections::inlinePerturb
240 (
241  const triSurface& surf1,
242  const scalarField& surf1PointTol, // surf1 tolerance per point
243  const label edgeI,
244  Random& rndGen,
245  pointField& points1,
246  boolList& affectedEdges
247 ) const
248 {
249  bool hasPerturbed = false;
250 
251  // Check if edge close to endpoint. Note that we only have to check
252  // the intersection closest to the edge endpoints (i.e. first and last in
253  // edgeEdgs)
254 
255  const labelList& edgeEnds = classification_[edgeI];
256 
257  if (edgeEnds.size())
258  {
259  bool perturbStart = false;
260  bool perturbEnd = false;
261 
262  // Check first intersection.
263  if (edgeEnds[0] == 0)
264  {
265  perturbStart = true;
266  }
267 
268  if (edgeEnds[edgeEnds.size()-1] == 1)
269  {
270  perturbEnd = true;
271  }
272 
273 
274  if (perturbStart || perturbEnd)
275  {
276  const edge& e = surf1.edges()[edgeI];
277 
278  label v0 = surf1.meshPoints()[e[0]];
279  label v1 = surf1.meshPoints()[e[1]];
280 
281  vector eVec(points1[v1] - points1[v0]);
282  vector n = eVec/mag(eVec);
283 
284  if (perturbStart)
285  {
286  // Perturb with something (hopefully) larger than tolerance.
287  scalar t = 4.0*(rndGen.scalar01() - 0.5);
288  points1[v0] += t*surf1PointTol[e[0]]*n;
289 
290  const labelList& pEdges = surf1.pointEdges()[e[0]];
291 
292  forAll(pEdges, i)
293  {
294  affectedEdges[pEdges[i]] = true;
295  }
296  }
297  if (perturbEnd)
298  {
299  // Perturb with something larger than tolerance.
300  scalar t = 4.0*(rndGen.scalar01() - 0.5);
301  points1[v1] += t*surf1PointTol[e[1]]*n;
302 
303  const labelList& pEdges = surf1.pointEdges()[e[1]];
304 
305  forAll(pEdges, i)
306  {
307  affectedEdges[pEdges[i]] = true;
308  }
309  }
310  hasPerturbed = true;
311  }
312  }
313 
314  return hasPerturbed;
315 }
316 
317 
318 // Perturb single edge endpoint when perpendicular to face
319 bool Foam::edgeIntersections::rotatePerturb
320 (
321  const triSurface& surf1,
322  const scalarField& surf1PointTol, // surf1 tolerance per point
323  const label edgeI,
324 
325  Random& rndGen,
326  pointField& points1,
327  boolList& affectedEdges
328 ) const
329 {
330  const labelList& meshPoints = surf1.meshPoints();
331 
332  const labelList& edgeEnds = classification_[edgeI];
333 
334  bool hasPerturbed = false;
335 
336  forAll(edgeEnds, i)
337  {
338  if (edgeEnds[i] == 2)
339  {
340  const edge& e = surf1.edges()[edgeI];
341 
342  // Endpoint to modify. Choose either start or end.
343  label pointI = e[rndGen.bit()];
344  //label pointI = e[0];
345 
346  // Generate random vector slightly larger than tolerance.
347  vector rndVec = rndGen.vector01() - vector(0.5, 0.5, 0.5);
348 
349  // Make sure rndVec only perp to edge
350  vector n(points1[meshPoints[e[1]]] - points1[meshPoints[e[0]]]);
351  scalar magN = mag(n) + VSMALL;
352  n /= magN;
353 
354  rndVec -= n*(n & rndVec);
355 
356  // Normalize
357  rndVec /= mag(rndVec) + VSMALL;
358 
359  // Scale to be moved by tolerance.
360  rndVec *= 0.01*magN;
361 
362  Pout<< "rotating: shifting endpoint " << meshPoints[pointI]
363  << " of edge:" << edgeI << " verts:"
364  << points1[meshPoints[e[0]]] << ' '
365  << points1[meshPoints[e[1]]]
366  << " by " << rndVec
367  << " tol:" << surf1PointTol[pointI] << endl;
368 
369  points1[meshPoints[pointI]] += rndVec;
370 
371  // Mark edges affected by change to point
372  const labelList& pEdges = surf1.pointEdges()[pointI];
373 
374  forAll(pEdges, i)
375  {
376  affectedEdges[pEdges[i]] = true;
377  }
378 
379  hasPerturbed = true;
380 
381  // Enough done for current edge; no need to test other intersections
382  // of this edge.
383  break;
384  }
385  }
386 
387  return hasPerturbed;
388 }
389 
390 
391 // Perturb edge when close to face
392 bool Foam::edgeIntersections::offsetPerturb
393 (
394  const triSurface& surf1,
395  const triSurface& surf2,
396  const label edgeI,
397 
398  Random& rndGen,
399  pointField& points1,
400  boolList& affectedEdges
401 ) const
402 {
403  const labelList& meshPoints = surf1.meshPoints();
404 
405  const edge& e = surf1.edges()[edgeI];
406 
407  const List<pointIndexHit>& hits = operator[](edgeI);
408 
409 
410  bool hasPerturbed = false;
411 
412  // For all hits on edge
413  forAll(hits, i)
414  {
415  const pointIndexHit& pHit = hits[i];
416 
417  // Classify point on face of surface2
418  label surf2FaceI = pHit.index();
419 
420  const labelledTri& f2 = surf2.localFaces()[surf2FaceI];
421 
422  const pointField& surf2Pts = surf2.localPoints();
423 
424  label nearType;
425  label nearLabel;
426 
427  triPointRef tri
428  (
429  surf2Pts[f2[0]],
430  surf2Pts[f2[1]],
431  surf2Pts[f2[2]]
432  );
433 
434  point ctr = tri.centre();
435 
436  // Get measure for tolerance.
437  scalar tolDim = 0.001*mag(tri.a() - ctr);
438 
439  tri.classify(pHit.hitPoint(), tolDim, nearType, nearLabel);
440 
441  if (nearType == triPointRef::POINT || nearType == triPointRef::EDGE)
442  {
443  // Shift edge towards tri centre
444  vector offset = 0.01*rndGen.scalar01()*(ctr - pHit.hitPoint());
445 
446  // shift e[0]
447  points1[meshPoints[e[0]]] += offset;
448 
449  // Mark edges affected by change to e0
450  const labelList& pEdges0 = surf1.pointEdges()[e[0]];
451 
452  forAll(pEdges0, i)
453  {
454  affectedEdges[pEdges0[i]] = true;
455  }
456 
457  // shift e[1]
458  points1[meshPoints[e[1]]] += offset;
459 
460  // Mark edges affected by change to e1
461  const labelList& pEdges1 = surf1.pointEdges()[e[1]];
462 
463  forAll(pEdges1, i)
464  {
465  affectedEdges[pEdges1[i]] = true;
466  }
467 
468  hasPerturbed = true;
469 
470  // No need to test any other hits on this edge
471  break;
472  }
473  }
474 
475  return hasPerturbed;
476 }
477 
478 
479 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
480 
481 // Construct null
483 :
484  List<List<pointIndexHit> >(),
485  classification_()
486 {}
487 
488 
489 // Construct from surface and tolerance
491 (
492  const triSurface& surf1,
493  const triSurfaceSearch& query2,
494  const scalarField& surf1PointTol
495 )
496 :
497  List<List<pointIndexHit> >(surf1.nEdges()),
498  classification_(surf1.nEdges())
499 {
500  checkEdges(surf1);
501  checkEdges(query2.surface());
502 
503  // Current set of edges to test
504  labelList edgesToTest(surf1.nEdges());
505 
506  // Start off with all edges
507  forAll(edgesToTest, i)
508  {
509  edgesToTest[i] = i;
510  }
511 
512  // Determine intersections for edgesToTest
513  intersectEdges
514  (
515  surf1,
516  surf1.points(), // surf1 meshPoints (not localPoints!)
517  query2,
518  surf1PointTol,
519  edgesToTest
520  );
521 }
522 
523 
524 // Construct from components
526 (
527  const List<List<pointIndexHit> >& intersections,
528  const labelListList& classification
529 )
530 :
531  List<List<pointIndexHit> >(intersections),
532  classification_(classification)
533 {}
534 
535 
536 // * * * * * * * * * * * * * * * Static Functions * * * * * * * * * * * * * //
537 
539 {
540  const pointField& localPoints = surf.localPoints();
541  const labelListList& pointEdges = surf.pointEdges();
542  const edgeList& edges = surf.edges();
543 
544  scalarField minLen(localPoints.size());
545 
546  forAll(minLen, pointI)
547  {
548  const labelList& pEdges = pointEdges[pointI];
549 
550  scalar minDist = GREAT;
551 
552  forAll(pEdges, i)
553  {
554  minDist = min(minDist, edges[pEdges[i]].mag(localPoints));
555  }
556 
557  minLen[pointI] = minDist;
558  }
559  return minLen;
560 }
561 
562 
563 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
564 
566 (
567  const label nIters,
568  const triSurface& surf1,
569  const triSurfaceSearch& query2,
570  const scalarField& surf1PointTol,
571  pointField& points1
572 )
573 {
574  const triSurface& surf2 = query2.surface();
575 
576  Random rndGen(356574);
577 
578  // Current set of edges to (re)test
579  labelList edgesToTest(surf1.nEdges());
580 
581  // Start off with all edges
582  forAll(edgesToTest, i)
583  {
584  edgesToTest[i] = i;
585  }
586 
587 
588  label iter = 0;
589 
590  for (; iter < nIters; iter++)
591  {
592  // Go through all edges to (re)test and perturb points if they are
593  // degenerate hits. Mark off edges that need to be recalculated.
594 
595  boolList affectedEdges(surf1.nEdges(), false);
596  label nShifted = 0;
597  label nRotated = 0;
598  label nOffset = 0;
599 
600  forAll(edgesToTest, i)
601  {
602  label edgeI = edgesToTest[i];
603 
604  // If edge not already marked for retesting
605  if (!affectedEdges[edgeI])
606  {
607  // 1. Check edges close to endpoint and perturb if nessecary.
608 
609  bool shiftedEdgeEndPoints =
610  inlinePerturb
611  (
612  surf1,
613  surf1PointTol,
614  edgeI,
615  rndGen,
616  points1,
617  affectedEdges
618  );
619 
620  nShifted += (shiftedEdgeEndPoints ? 1 : 0);
621 
622  if (!shiftedEdgeEndPoints)
623  {
624  bool rotatedEdge =
625  rotatePerturb
626  (
627  surf1,
628  surf1PointTol,
629  edgeI,
630  rndGen,
631  points1,
632  affectedEdges
633  );
634 
635  nRotated += (rotatedEdge ? 1 : 0);
636 
637  if (!rotatedEdge)
638  {
639  // 2. we're sure now that the edge actually pierces the
640  // face. Now check the face for intersections close its
641  // points/edges
642 
643  bool offsetEdgePoints =
644  offsetPerturb
645  (
646  surf1,
647  surf2,
648  edgeI,
649  rndGen,
650  points1,
651  affectedEdges
652  );
653 
654  nOffset += (offsetEdgePoints ? 1 : 0);
655  }
656  }
657  }
658  }
659 
660  if (debug)
661  {
662  Pout<< "Edges to test : " << nl
663  << " total:" << edgesToTest.size() << nl
664  << " resolved by:" << nl
665  << " shifting : " << nShifted << nl
666  << " rotating : " << nRotated << nl
667  << " offsetting : " << nOffset << nl
668  << endl;
669  }
670 
671 
672  if (nShifted == 0 && nRotated == 0 && nOffset == 0)
673  {
674  // Nothing changed in current iteration. Current hit pattern is
675  // without any degenerates.
676  break;
677  }
678 
679  // Repack affected edges
680  labelList newEdgesToTest(surf1.nEdges());
681  label newEdgeI = 0;
682 
683  forAll(affectedEdges, edgeI)
684  {
685  if (affectedEdges[edgeI])
686  {
687  newEdgesToTest[newEdgeI++] = edgeI;
688  }
689  }
690  newEdgesToTest.setSize(newEdgeI);
691 
692  if (debug)
693  {
694  Pout<< "Edges to test:" << nl
695  << " was : " << edgesToTest.size() << nl
696  << " is : " << newEdgesToTest.size() << nl
697  << endl;
698  }
699 
700  // Transfer and test.
701  edgesToTest.transfer(newEdgesToTest);
702 
703  if (edgesToTest.empty())
704  {
705  FatalErrorIn("perturb") << "oops" << abort(FatalError);
706  }
707 
708  // Re intersect moved edges.
709  intersectEdges
710  (
711  surf1,
712  points1, // surf1 meshPoints (not localPoints!)
713  query2,
714  surf1PointTol,
715  edgesToTest
716  );
717  }
718 
719  return iter;
720 }
721 
722 
723 // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
724 
725 
726 // * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * //
727 
728 
729 // * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
730 
731 
732 // ************************ vim: set sw=4 sts=4 et: ************************ //