intersectedSurface.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-2012 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 "intersectedSurface.H"
27 #include "surfaceIntersection.H"
28 #include "faceList.H"
29 #include "faceTriangulation.H"
30 #include "treeBoundBox.H"
31 #include "OFstream.H"
32 #include "error.H"
33 #include "meshTools.H"
34 #include "edgeSurface.H"
35 #include "DynamicList.H"
36 #include "transform.H"
37 
38 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
39 
40 namespace Foam
41 {
42 defineTypeNameAndDebug(intersectedSurface, 0);
43 }
44 
45 
49 const Foam::label Foam::intersectedSurface::BOTH = STARTTOEND | ENDTOSTART;
50 
51 
52 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
53 
54 // Write whole pointField and edges to stream
55 void Foam::intersectedSurface::writeOBJ
56 (
57  const pointField& points,
58  const edgeList& edges,
59  Ostream& os
60 )
61 {
62  forAll(points, pointI)
63  {
64  const point& pt = points[pointI];
65 
66  os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << nl;
67  }
68  forAll(edges, edgeI)
69  {
70  const edge& e = edges[edgeI];
71 
72  os << "l " << e.start()+1 << ' ' << e.end()+1 << nl;
73  }
74 }
75 
76 
77 // Write whole pointField and selected edges to stream
78 void Foam::intersectedSurface::writeOBJ
79 (
80  const pointField& points,
81  const edgeList& edges,
82  const labelList& faceEdges,
83  Ostream& os
84 )
85 {
86  forAll(points, pointI)
87  {
88  const point& pt = points[pointI];
89 
90  os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << nl;
91  }
92  forAll(faceEdges, i)
93  {
94  const edge& e = edges[faceEdges[i]];
95 
96  os << "l " << e.start()+1 << ' ' << e.end()+1 << nl;
97  }
98 }
99 
100 
101 // write local points and edges to stream
102 void Foam::intersectedSurface::writeLocalOBJ
103 (
104  const pointField& points,
105  const edgeList& edges,
106  const labelList& faceEdges,
107  const fileName& fName
108 )
109 {
110  OFstream os(fName);
111 
112  labelList pointMap(points.size(), -1);
113 
114  label maxVertI = 0;
115 
116  forAll(faceEdges, i)
117  {
118  const edge& e = edges[faceEdges[i]];
119 
120  forAll(e, i)
121  {
122  label pointI = e[i];
123 
124  if (pointMap[pointI] == -1)
125  {
126  const point& pt = points[pointI];
127 
128  os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << nl;
129 
130  pointMap[pointI] = maxVertI++;
131  }
132  }
133  }
134 
135  forAll(faceEdges, i)
136  {
137  const edge& e = edges[faceEdges[i]];
138 
139  os << "l " << pointMap[e.start()]+1 << ' ' << pointMap[e.end()]+1
140  << nl;
141  }
142 }
143 
144 
145 // Write whole pointField and face to stream
146 void Foam::intersectedSurface::writeOBJ
147 (
148  const pointField& points,
149  const face& f,
150  Ostream& os
151 )
152 {
153  forAll(points, pointI)
154  {
155  const point& pt = points[pointI];
156 
157  os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << nl;
158  }
159 
160  os << 'f';
161 
162  forAll(f, fp)
163  {
164  os << ' ' << f[fp]+1;
165  }
166  os << nl;
167 }
168 
169 
170 // Print current visited state.
171 void Foam::intersectedSurface::printVisit
172 (
173  const edgeList& edges,
174  const labelList& edgeLabels,
175  const Map<label>& visited
176 )
177 {
178  Pout<< "Visited:" << nl;
179  forAll(edgeLabels, i)
180  {
181  label edgeI = edgeLabels[i];
182 
183  const edge& e = edges[edgeI];
184 
185  label stat = visited[edgeI];
186 
187  if (stat == UNVISITED)
188  {
189  Pout<< " edge:" << edgeI << " verts:" << e
190  << " unvisited" << nl;
191  }
192  else if (stat == STARTTOEND)
193  {
194  Pout<< " edge:" << edgeI << " from " << e[0]
195  << " to " << e[1] << nl;
196  }
197  else if (stat == ENDTOSTART)
198  {
199  Pout<< " edge:" << edgeI << " from " << e[1]
200  << " to " << e[0] << nl;
201  }
202  else
203  {
204  Pout<< " edge:" << edgeI << " both " << e
205  << nl;
206  }
207  }
208  Pout<< endl;
209 }
210 
211 
212 // Check if the two vertices that f0 and f1 share are in the same order on
213 // both faces.
214 bool Foam::intersectedSurface::sameEdgeOrder
215 (
216  const labelledTri& fA,
217  const labelledTri& fB
218 )
219 {
220  forAll(fA, fpA)
221  {
222  label fpB = findIndex(fB, fA[fpA]);
223 
224  if (fpB != -1)
225  {
226  // Get prev/next vertex on fA
227  label vA1 = fA[fA.fcIndex(fpA)];
228  label vAMin1 = fA[fA.rcIndex(fpA)];
229 
230  // Get prev/next vertex on fB
231  label vB1 = fB[fB.fcIndex(fpB)];
232  label vBMin1 = fB[fB.rcIndex(fpB)];
233 
234  if (vA1 == vB1 || vAMin1 == vBMin1)
235  {
236  return true;
237  }
238  else if (vA1 == vBMin1 || vAMin1 == vB1)
239  {
240  // shared vertices in opposite order.
241  return false;
242  }
243  else
244  {
245  FatalErrorIn("intersectedSurface::sameEdgeOrder")
246  << "Triangle:" << fA << " and triangle:" << fB
247  << " share a point but not an edge"
248  << abort(FatalError);
249  }
250  }
251  }
252 
253  FatalErrorIn("intersectedSurface::sameEdgeOrder")
254  << "Triangle:" << fA << " and triangle:" << fB
255  << " do not share an edge"
256  << abort(FatalError);
257 
258  return false;
259 }
260 
261 
262 void Foam::intersectedSurface::incCount
263 (
264  Map<label>& visited,
265  const label key,
266  const label offset
267 )
268 {
269  Map<label>::iterator iter = visited.find(key);
270 
271  if (iter == visited.end())
272  {
273  visited.insert(key, offset);
274  }
275  else
276  {
277  iter() += offset;
278  }
279 }
280 
281 
282 // Calculate point to edge addressing for the face given by the edge
283 // subset faceEdges. Constructs facePointEdges which for every point
284 // gives a list of edge labels connected to it.
286 Foam::intersectedSurface::calcPointEdgeAddressing
287 (
288  const edgeSurface& eSurf,
289  const label faceI
290 )
291 {
292  const pointField& points = eSurf.points();
293  const edgeList& edges = eSurf.edges();
294 
295  const labelList& fEdges = eSurf.faceEdges()[faceI];
296 
297  Map<DynamicList<label> > facePointEdges(4*fEdges.size());
298 
299  forAll(fEdges, i)
300  {
301  label edgeI = fEdges[i];
302 
303  const edge& e = edges[edgeI];
304 
305  // Add e.start to point-edges
306  Map<DynamicList<label> >::iterator iter =
307  facePointEdges.find(e.start());
308 
309  if (iter == facePointEdges.end())
310  {
311  DynamicList<label> oneEdge;
312  oneEdge.append(edgeI);
313  facePointEdges.insert(e.start(), oneEdge);
314  }
315  else
316  {
317  iter().append(edgeI);
318  }
319 
320  // Add e.end to point-edges
321  Map<DynamicList<label> >::iterator iter2 =
322  facePointEdges.find(e.end());
323 
324  if (iter2 == facePointEdges.end())
325  {
326  DynamicList<label> oneEdge;
327  oneEdge.append(edgeI);
328  facePointEdges.insert(e.end(), oneEdge);
329  }
330  else
331  {
332  iter2().append(edgeI);
333  }
334  }
335 
336  // Shrink it
337  forAllIter(Map< DynamicList<label> >, facePointEdges, iter)
338  {
339  iter().shrink();
340 
341  // Check on dangling points.
342  if (iter().empty())
343  {
345  (
346  "intersectedSurface::calcPointEdgeAddressing"
347  "(const edgeSurface&, const label)"
348  ) << "Point:" << iter.key() << " used by too few edges:"
349  << iter() << abort(FatalError);
350  }
351  }
352 
353  if (debug & 2)
354  {
355  // Print facePointEdges
356  Pout<< "calcPointEdgeAddressing: face consisting of edges:" << endl;
357  forAll(fEdges, i)
358  {
359  label edgeI = fEdges[i];
360  const edge& e = edges[edgeI];
361  Pout<< " " << edgeI << ' ' << e
362  << points[e.start()]
363  << points[e.end()] << endl;
364  }
365 
366  Pout<< " Constructed point-edge adressing:" << endl;
367  forAllConstIter(Map< DynamicList<label> >, facePointEdges, iter)
368  {
369  Pout<< " vertex " << iter.key() << " is connected to edges "
370  << iter() << endl;
371  }
372  Pout<<endl;
373  }
374 
375  return facePointEdges;
376 }
377 
378 
379 // Find next (triangle or cut) edge label coming from point prevVertI on
380 // prevEdgeI doing a right handed walk (i.e. following right wall).
381 // Note: normal is provided externally. Could be deducted from angle between
382 // two triangle edges but these could be in line.
383 Foam::label Foam::intersectedSurface::nextEdge
384 (
385  const edgeSurface& eSurf,
386  const Map<label>& visited,
387  const label faceI,
388  const vector& n,
389  const Map<DynamicList<label> >& facePointEdges,
390  const label prevEdgeI,
391  const label prevVertI
392 )
393 {
394  const pointField& points = eSurf.points();
395  const edgeList& edges = eSurf.edges();
396  const labelList& fEdges = eSurf.faceEdges()[faceI];
397 
398 
399  // Edges connected to prevVertI
400  const DynamicList<label>& connectedEdges = facePointEdges[prevVertI];
401 
402  if (connectedEdges.size() <= 1)
403  {
404  // Problem. Point not connected.
405  {
406  Pout<< "Writing face:" << faceI << " to face.obj" << endl;
407  OFstream str("face.obj");
408  writeOBJ(points, edges, fEdges, str);
409 
410  Pout<< "Writing connectedEdges edges to faceEdges.obj" << endl;
411  writeLocalOBJ(points, edges, connectedEdges, "faceEdges.obj");
412  }
413 
415  (
416  "intersectedSurface::nextEdge(const pointField&, const edgeList&"
417  ", const vector&, Map<DynamicList<label> >, const label"
418  ", const label)"
419  ) << "Problem: prevVertI:" << prevVertI << " on edge " << prevEdgeI
420  << " has less than 2 connected edges."
421  << " connectedEdges:" << connectedEdges << abort(FatalError);
422 
423  return -1;
424  }
425 
426  if (connectedEdges.size() == 2)
427  {
428  // Simple case. Take other edge
429  if (connectedEdges[0] == prevEdgeI)
430  {
431  if (debug & 2)
432  {
433  Pout<< "Stepped from edge:" << edges[prevEdgeI]
434  << " to edge:" << edges[connectedEdges[1]]
435  << " chosen from candidates:" << connectedEdges << endl;
436  }
437  return connectedEdges[1];
438  }
439  else
440  {
441  if (debug & 2)
442  {
443  Pout<< "Stepped from edge:" << edges[prevEdgeI]
444  << " to edge:" << edges[connectedEdges[0]]
445  << " chosen from candidates:" << connectedEdges << endl;
446  }
447  return connectedEdges[0];
448  }
449  }
450 
451 
452  // Multiple choices. Look at angle between edges.
453 
454  const edge& prevE = edges[prevEdgeI];
455 
456  // x-axis of coordinate system
457  vector e0 = n ^ (points[prevE.otherVertex(prevVertI)] - points[prevVertI]);
458  e0 /= mag(e0) + VSMALL;
459 
460  // Get y-axis of coordinate system
461  vector e1 = n ^ e0;
462 
463  if (mag(mag(e1) - 1) > SMALL)
464  {
465  {
466  Pout<< "Writing face:" << faceI << " to face.obj" << endl;
467  OFstream str("face.obj");
468  writeOBJ(points, edges, fEdges, str);
469 
470  Pout<< "Writing connectedEdges edges to faceEdges.obj" << endl;
471  writeLocalOBJ(points, edges, connectedEdges, "faceEdges.obj");
472  }
473 
474  FatalErrorIn("intersectedSurface::nextEdge")
475  << "Unnormalized normal e1:" << e1
476  << " formed from cross product of e0:" << e0 << " n:" << n
477  << abort(FatalError);
478  }
479 
480 
481  //
482  // Determine maximum angle over all connected (unvisited) edges.
483  //
484 
485  scalar maxAngle = -GREAT;
486  label maxEdgeI = -1;
487 
488  forAll(connectedEdges, connI)
489  {
490  label edgeI = connectedEdges[connI];
491 
492  if (edgeI != prevEdgeI)
493  {
494  label stat = visited[edgeI];
495 
496  const edge& e = edges[edgeI];
497 
498  // Find out whether walk of edge from prevVert would be acceptible.
499  if
500  (
501  stat == UNVISITED
502  || (
503  stat == STARTTOEND
504  && prevVertI == e[1]
505  )
506  || (
507  stat == ENDTOSTART
508  && prevVertI == e[0]
509  )
510  )
511  {
512  // Calculate angle of edge with respect to base e0, e1
513  vector vec =
514  n ^ (points[e.otherVertex(prevVertI)] - points[prevVertI]);
515 
516  vec /= mag(vec) + VSMALL;
517 
518  scalar angle = pseudoAngle(e0, e1, vec);
519 
520  if (angle > maxAngle)
521  {
522  maxAngle = angle;
523  maxEdgeI = edgeI;
524  }
525  }
526  }
527  }
528 
529 
530  if (maxEdgeI == -1)
531  {
532  // No unvisited edge found
533  {
534  Pout<< "Writing face:" << faceI << " to face.obj" << endl;
535  OFstream str("face.obj");
536  writeOBJ(points, edges, fEdges, str);
537 
538  Pout<< "Writing connectedEdges edges to faceEdges.obj" << endl;
539  writeLocalOBJ(points, edges, connectedEdges, "faceEdges.obj");
540  }
541 
543  (
544  "intersectedSurface::nextEdge(const pointField&, const edgeList&"
545  ", const Map<label>&, const vector&"
546  ", const Map<DynamicList<label> >&"
547  ", const label, const label"
548  ) << "Trying to step from edge " << edges[prevEdgeI]
549  << ", vertex " << prevVertI
550  << " but cannot find 'unvisited' edges among candidates:"
551  << connectedEdges
552  << abort(FatalError);
553  }
554 
555  if (debug & 2)
556  {
557  Pout<< "Stepped from edge:" << edges[prevEdgeI]
558  << " to edge:" << maxEdgeI << " angle:" << edges[maxEdgeI]
559  << " with angle:" << maxAngle
560  << endl;
561  }
562 
563  return maxEdgeI;
564 }
565 
566 
567 // Create (polygonal) face by walking full circle starting from startVertI on
568 // startEdgeI.
569 // Uses nextEdge(..) to do the walking. Returns face. Updates visited with
570 // the labels of visited edges.
571 Foam::face Foam::intersectedSurface::walkFace
572 (
573  const edgeSurface& eSurf,
574  const label faceI,
575  const vector& n,
576  const Map<DynamicList<label> >& facePointEdges,
577 
578  const label startEdgeI,
579  const label startVertI,
580 
581  Map<label>& visited
582 )
583 {
584  const pointField& points = eSurf.points();
585  const edgeList& edges = eSurf.edges();
586 
587  // Overestimate size of face
588  face f(eSurf.faceEdges()[faceI].size());
589 
590  label fp = 0;
591 
592  label vertI = startVertI;
593  label edgeI = startEdgeI;
594 
595  while (true)
596  {
597  const edge& e = edges[edgeI];
598 
599  if (debug & 2)
600  {
601  Pout<< "Now at:" << endl
602  << " edge:" << edgeI << " vertices:" << e
603  << " positions:" << points[e.start()] << ' ' << points[e.end()]
604  << " vertex:" << vertI << endl;
605  }
606 
607  // Mark edge as visited
608  if (e[0] == vertI)
609  {
610  visited[edgeI] |= STARTTOEND;
611  }
612  else
613  {
614  visited[edgeI] |= ENDTOSTART;
615  }
616 
617 
618  // Store face vertex
619  f[fp++] = vertI;
620 
621  // step to other vertex
622  vertI = e.otherVertex(vertI);
623 
624  if (vertI == startVertI)
625  {
626  break;
627  }
628 
629  // step from vertex to next edge
630  edgeI = nextEdge
631  (
632  eSurf,
633  visited,
634  faceI,
635  n,
636  facePointEdges,
637  edgeI,
638  vertI
639  );
640  }
641 
642  f.setSize(fp);
643 
644  return f;
645 }
646 
647 
648 void Foam::intersectedSurface::findNearestVisited
649 (
650  const edgeSurface& eSurf,
651  const label faceI,
652  const Map<DynamicList<label> >& facePointEdges,
653  const Map<label>& pointVisited,
654  const point& pt,
655  const label excludePointI,
656 
657  label& minVertI,
658  scalar& minDist
659 )
660 {
661  minVertI = -1;
662  minDist = GREAT;
663 
664  forAllConstIter(Map<label>, pointVisited, iter)
665  {
666  label pointI = iter.key();
667 
668  if (pointI != excludePointI)
669  {
670  label nVisits = iter();
671 
672  if (nVisits == 2*facePointEdges[pointI].size())
673  {
674  // Fully visited (i.e. both sides of all edges)
675  scalar dist = mag(eSurf.points()[pointI] - pt);
676 
677  if (dist < minDist)
678  {
679  minDist = dist;
680  minVertI = pointI;
681  }
682  }
683  }
684  }
685 
686  if (minVertI == -1)
687  {
688  const labelList& fEdges = eSurf.faceEdges()[faceI];
689 
690  SeriousErrorIn("intersectedSurface::findNearestVisited")
691  << "Dumping face edges to faceEdges.obj" << endl;
692 
693  writeLocalOBJ(eSurf.points(), eSurf.edges(), fEdges, "faceEdges.obj");
694 
695  FatalErrorIn("intersectedSurface::findNearestVisited")
696  << "No fully visited edge found for pt " << pt
697  << abort(FatalError);
698  }
699 }
700 
701 
702 // Sometimes there are bunches of edges that are not connected to the
703 // other edges in the face. This routine tries to connect the loose
704 // edges up to the 'proper' edges by adding two extra edges between a
705 // properly connected edge and an unconnected one. Since at this level the
706 // adressing is purely in form of points and a cloud of edges this can
707 // be easily done.
708 // (edges are to existing points. Don't want to introduce new vertices here
709 // since then also neighbouring face would have to be split)
710 Foam::faceList Foam::intersectedSurface::resplitFace
711 (
712  const triSurface& surf,
713  const label faceI,
714  const Map<DynamicList<label> >& facePointEdges,
715  const Map<label>& visited,
716  edgeSurface& eSurf
717 )
718 {
719  {
720  // Dump face for debugging.
721  Pout<< "Writing face:" << faceI << " to face.obj" << endl;
722  OFstream str("face.obj");
723  writeOBJ(eSurf.points(), eSurf.edges(), eSurf.faceEdges()[faceI], str);
724  }
725 
726 
727  // Count the number of times point has been visited so we
728  // can compare number to facePointEdges.
729  Map<label> pointVisited(2*facePointEdges.size());
730 
731 
732  {
733  OFstream str0("visitedNone.obj");
734  OFstream str1("visitedOnce.obj");
735  OFstream str2("visitedTwice.obj");
736  forAll(eSurf.points(), pointI)
737  {
738  const point& pt = eSurf.points()[pointI];
739 
740  str0 << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << nl;
741  str1 << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << nl;
742  str2 << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << nl;
743  }
744 
745 
746  forAllConstIter(Map<label>, visited, iter)
747  {
748  label edgeI = iter.key();
749 
750  const edge& e = eSurf.edges()[edgeI];
751 
752  label stat = iter();
753 
754  if (stat == STARTTOEND || stat == ENDTOSTART)
755  {
756  incCount(pointVisited, e[0], 1);
757  incCount(pointVisited, e[1], 1);
758 
759  str1 << "l " << e[0]+1 << ' ' << e[1]+1 << nl;
760  }
761  else if (stat == BOTH)
762  {
763  incCount(pointVisited, e[0], 2);
764  incCount(pointVisited, e[1], 2);
765 
766  str2 << "l " << e[0]+1 << ' ' << e[1]+1 << nl;
767  }
768  else if (stat == UNVISITED)
769  {
770  incCount(pointVisited, e[0], 0);
771  incCount(pointVisited, e[1], 0);
772 
773  str0 << "l " << e[0]+1 << ' ' << e[1]+1 << nl;
774  }
775  }
776  }
777 
778 
779  {
780  forAllConstIter(Map<label>, pointVisited, iter)
781  {
782  label pointI = iter.key();
783 
784  label nVisits = iter();
785 
786  Pout<< "point:" << pointI << " nVisited:" << nVisits
787  << " pointEdges:" << facePointEdges[pointI].size() << endl;
788  }
789  }
790 
791 
792  // Find nearest point pair where one is not fully visited and
793  // the other is.
794 
795  label visitedVert0 = -1;
796  label unvisitedVert0 = -1;
797 
798  {
799  scalar minDist = GREAT;
800 
801  forAllConstIter(Map<DynamicList<label> >, facePointEdges, iter)
802  {
803  label pointI = iter.key();
804 
805  label nVisits = pointVisited[pointI];
806 
807  const DynamicList<label>& pEdges = iter();
808 
809  if (nVisits < 2*pEdges.size())
810  {
811  // Not fully visited. Find nearest fully visited.
812 
813  scalar nearDist;
814  label nearVertI;
815 
816  findNearestVisited
817  (
818  eSurf,
819  faceI,
820  facePointEdges,
821  pointVisited,
822  eSurf.points()[pointI],
823  -1, // Do not exclude vertex
824  nearVertI,
825  nearDist
826  );
827 
828 
829  if (nearDist < minDist)
830  {
831  minDist = nearDist;
832  visitedVert0 = nearVertI;
833  unvisitedVert0 = pointI;
834  }
835  }
836  }
837  }
838 
839 
840  // Find second intersection.
841  label visitedVert1 = -1;
842  label unvisitedVert1 = -1;
843  {
844  scalar minDist = GREAT;
845 
846  forAllConstIter(Map<DynamicList<label> >, facePointEdges, iter)
847  {
848  label pointI = iter.key();
849 
850  if (pointI != unvisitedVert0)
851  {
852  label nVisits = pointVisited[pointI];
853 
854  const DynamicList<label>& pEdges = iter();
855 
856  if (nVisits < 2*pEdges.size())
857  {
858  // Not fully visited. Find nearest fully visited.
859 
860  scalar nearDist;
861  label nearVertI;
862 
863  findNearestVisited
864  (
865  eSurf,
866  faceI,
867  facePointEdges,
868  pointVisited,
869  eSurf.points()[pointI],
870  visitedVert0, // vertex to exclude
871  nearVertI,
872  nearDist
873  );
874 
875 
876  if (nearDist < minDist)
877  {
878  minDist = nearDist;
879  visitedVert1 = nearVertI;
880  unvisitedVert1 = pointI;
881  }
882  }
883  }
884  }
885  }
886 
887 
888  Pout<< "resplitFace : adding intersection from " << visitedVert0
889  << " to " << unvisitedVert0 << endl
890  << " and from " << visitedVert1
891  << " to " << unvisitedVert1 << endl;
892 
893 
894 // // Add the new intersection edges to the edgeSurface
895 // edgeList additionalEdges(2);
896 // additionalEdges[0] = edge(visitedVert0, unvisitedVert0);
897 // additionalEdges[1] = edge(visitedVert1, unvisitedVert1);
898 
899  // Add the new intersection edges to the edgeSurface
900  edgeList additionalEdges(1);
901  additionalEdges[0] = edge(visitedVert0, unvisitedVert0);
902 
903  eSurf.addIntersectionEdges(faceI, additionalEdges);
904 
905  fileName newFName("face_" + Foam::name(faceI) + "_newEdges.obj");
906  Pout<< "Dumping face:" << faceI << " to " << newFName << endl;
907  writeLocalOBJ
908  (
909  eSurf.points(),
910  eSurf.edges(),
911  eSurf.faceEdges()[faceI],
912  newFName
913  );
914 
915  // Retry splitFace. Use recursion since is rare situation.
916  return splitFace(surf, faceI, eSurf);
917 }
918 
919 
920 Foam::faceList Foam::intersectedSurface::splitFace
921 (
922  const triSurface& surf,
923  const label faceI,
924  edgeSurface& eSurf
925 )
926 {
927  // Alias
928  const pointField& points = eSurf.points();
929  const edgeList& edges = eSurf.edges();
930  const labelList& fEdges = eSurf.faceEdges()[faceI];
931 
932  // Create local (for the face only) point-edge connectivity.
933  Map<DynamicList<label> > facePointEdges
934  (
935  calcPointEdgeAddressing
936  (
937  eSurf,
938  faceI
939  )
940  );
941 
942  // Order in which edges have been walked. Initialize outside edges.
943  Map<label> visited(fEdges.size()*2);
944 
945  forAll(fEdges, i)
946  {
947  label edgeI = fEdges[i];
948 
949  if (eSurf.isSurfaceEdge(edgeI))
950  {
951  // Edge is edge from original surface so an outside edge for
952  // the current face.
953  label surfEdgeI = eSurf.parentEdge(edgeI);
954 
955  label owner = surf.edgeOwner()[surfEdgeI];
956 
957  if
958  (
959  owner == faceI
960  || sameEdgeOrder
961  (
962  surf.localFaces()[owner],
963  surf.localFaces()[faceI]
964  )
965  )
966  {
967  // Edge is in same order as current face.
968  // Mark off the opposite order.
969  visited.insert(edgeI, ENDTOSTART);
970  }
971  else
972  {
973  // Edge is in same order as current face.
974  // Mark off the opposite order.
975  visited.insert(edgeI, STARTTOEND);
976  }
977  }
978  else
979  {
980  visited.insert(edgeI, UNVISITED);
981  }
982  }
983 
984 
985 
986  // Storage for faces
987  DynamicList<face> faces(fEdges.size());
988 
989  while (true)
990  {
991  // Find starting edge:
992  // - unvisted triangle edge
993  // - once visited intersection edge
994  // Give priority to triangle edges.
995  label startEdgeI = -1;
996  label startVertI = -1;
997 
998  forAll(fEdges, i)
999  {
1000  label edgeI = fEdges[i];
1001 
1002  const edge& e = edges[edgeI];
1003 
1004  label stat = visited[edgeI];
1005 
1006  if (stat == STARTTOEND)
1007  {
1008  startEdgeI = edgeI;
1009  startVertI = e[1];
1010 
1011  if (eSurf.isSurfaceEdge(edgeI))
1012  {
1013  break;
1014  }
1015  }
1016  else if (stat == ENDTOSTART)
1017  {
1018  startEdgeI = edgeI;
1019  startVertI = e[0];
1020 
1021  if (eSurf.isSurfaceEdge(edgeI))
1022  {
1023  break;
1024  }
1025  }
1026  }
1027 
1028  if (startEdgeI == -1)
1029  {
1030  break;
1031  }
1032 
1033  //Pout<< "splitFace: starting walk from edge:" << startEdgeI
1034  // << ' ' << edges[startEdgeI] << " vertex:" << startVertI << endl;
1035 
1037  //printVisit(eSurf.edges(), fEdges, visited);
1038 
1039  //{
1040  // Pout<< "Writing face:" << faceI << " to face.obj" << endl;
1041  // OFstream str("face.obj");
1042  // writeOBJ(eSurf.points(), eSurf.edges(), fEdges, str);
1043  //}
1044 
1045  faces.append
1046  (
1047  walkFace
1048  (
1049  eSurf,
1050  faceI,
1051  surf.faceNormals()[faceI],
1052  facePointEdges,
1053 
1054  startEdgeI,
1055  startVertI,
1056 
1057  visited
1058  )
1059  );
1060  }
1061 
1062 
1063  // Check if any unvisited edges left.
1064  forAll(fEdges, i)
1065  {
1066  label edgeI = fEdges[i];
1067 
1068  label stat = visited[edgeI];
1069 
1070  if (eSurf.isSurfaceEdge(edgeI) && stat != BOTH)
1071  {
1072  SeriousErrorIn("Foam::intersectedSurface::splitFace")
1073  << "Dumping face edges to faceEdges.obj" << endl;
1074 
1075  writeLocalOBJ(points, edges, fEdges, "faceEdges.obj");
1076 
1077  FatalErrorIn("intersectedSurface::splitFace")
1078  << "Problem: edge " << edgeI << " vertices "
1079  << edges[edgeI] << " on face " << faceI
1080  << " has visited status " << stat << " from a "
1081  << "righthanded walk along all"
1082  << " of the triangle edges. Are the original surfaces"
1083  << " closed and non-intersecting?"
1084  << abort(FatalError);
1085  }
1086  else if (stat != BOTH)
1087  {
1088  //{
1089  // Pout<< "Dumping faces so far to faces.obj" << nl
1090  // << faces << endl;
1091  //
1092  // OFstream str("faces.obj");
1093  //
1094  // forAll(faces, i)
1095  // {
1096  // writeOBJ(points, faces[i], str);
1097  // }
1098  //}
1099 
1100  Pout<< "** Resplitting **" << endl;
1101 
1102  // Redo face after introducing extra edge. Edge introduced
1103  // should be one nearest to any fully visited edge.
1104  return resplitFace
1105  (
1106  surf,
1107  faceI,
1108  facePointEdges,
1109  visited,
1110  eSurf
1111  );
1112  }
1113  }
1114 
1115 
1116  // See if normal needs flipping.
1117  faces.shrink();
1118 
1119  vector n = faces[0].normal(eSurf.points());
1120 
1121  if ((n & surf.faceNormals()[faceI]) < 0)
1122  {
1123  forAll(faces, i)
1124  {
1125  reverse(faces[i]);
1126  }
1127  }
1128 
1129  return faces;
1130 }
1131 
1132 
1133 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
1134 
1135 // Null constructor
1137 :
1138  triSurface(),
1139  intersectionEdges_(0),
1140  faceMap_(0),
1141  nSurfacePoints_(0)
1142 {}
1143 
1144 
1145 // Construct from components
1147 :
1148  triSurface(surf),
1149  intersectionEdges_(0),
1150  faceMap_(0),
1151  nSurfacePoints_(0)
1152 {}
1153 
1154 
1155 // Construct from surface and intersection
1158  const triSurface& surf,
1159  const bool isFirstSurface,
1160  const surfaceIntersection& inter
1161 )
1162 :
1163  triSurface(),
1164  intersectionEdges_(0),
1165  faceMap_(0),
1166  nSurfacePoints_(surf.nPoints())
1167 {
1168  if (inter.cutPoints().empty() && inter.cutEdges().empty())
1169  {
1170  // No intersection. Make straight copy.
1171  triSurface::operator=(surf);
1172 
1173  // Identity for face map
1174  faceMap_.setSize(size());
1175 
1176  forAll(faceMap_, faceI)
1177  {
1178  faceMap_[faceI] = faceI;
1179  }
1180  return;
1181  }
1182 
1183 
1184  // Calculate face-edge addressing for surface + intersection.
1185  edgeSurface eSurf(surf, isFirstSurface, inter);
1186 
1187  // Update point information for any removed points. (when are they removed?
1188  // - but better make sure)
1189  nSurfacePoints_ = eSurf.nSurfacePoints();
1190 
1191  // Now we have full points, edges and edge addressing for surf. Start
1192  // extracting faces and triangulate them.
1193 
1194  // Storage for new triangles of surface.
1195  DynamicList<labelledTri> newTris(eSurf.edges().size()/2);
1196 
1197  // Start in newTris for decomposed face.
1198  labelList startTriI(surf.size(), 0);
1199 
1200  forAll(surf, faceI)
1201  {
1202  startTriI[faceI] = newTris.size();
1203 
1204  if (eSurf.faceEdges()[faceI].size() != surf.faceEdges()[faceI].size())
1205  {
1206  // Face has been cut by intersection.
1207  // Cut face into multiple subfaces. Use faceEdge information
1208  // from edgeSurface only. (original triSurface 'surf' is used for
1209  // faceNormals and regionnumber only)
1210  faceList newFaces
1211  (
1212  splitFace
1213  (
1214  surf,
1215  faceI, // current triangle
1216  eSurf // face-edge description of surface
1217  // + intersection
1218  )
1219  );
1220  forAll(newFaces, newFaceI)
1221  {
1222  const face& newF = newFaces[newFaceI];
1223 
1224 // {
1225 // fileName fName
1226 // (
1227 // "face_"
1228 // + Foam::name(faceI)
1229 // + "_subFace_"
1230 // + Foam::name(newFaceI)
1231 // + ".obj"
1232 // );
1233 // Pout<< "Writing original face:" << faceI << " subFace:"
1234 // << newFaceI << " to " << fName << endl;
1235 //
1236 // OFstream str(fName);
1237 //
1238 // forAll(newF, fp)
1239 // {
1240 // meshTools::writeOBJ(str, eSurf.points()[newF[fp]]);
1241 // }
1242 // str << 'l';
1243 // forAll(newF, fp)
1244 // {
1245 // str << ' ' << fp+1;
1246 // }
1247 // str<< " 1" << nl;
1248 // }
1249 
1250 
1251  const vector& n = surf.faceNormals()[faceI];
1252  const label region = surf[faceI].region();
1253 
1254  faceTriangulation tris(eSurf.points(), newF, n);
1255 
1256  forAll(tris, triI)
1257  {
1258  const triFace& t = tris[triI];
1259 
1260  forAll(t, i)
1261  {
1262  if (t[i] < 0 || t[i] >= eSurf.points().size())
1263  {
1264  FatalErrorIn
1265  (
1266  "intersectedSurface::intersectedSurface"
1267  ) << "Face triangulation of face " << faceI
1268  << " uses points outside range 0.."
1269  << eSurf.points().size()-1 << endl
1270  << "Triangulation:"
1271  << tris << abort(FatalError);
1272  }
1273  }
1274 
1275  newTris.append(labelledTri(t[0], t[1], t[2], region));
1276  }
1277  }
1278  }
1279  else
1280  {
1281  // Face has not been cut at all. No need to renumber vertices since
1282  // eSurf keeps surface vertices first.
1283  newTris.append(surf.localFaces()[faceI]);
1284  }
1285  }
1286 
1287  newTris.shrink();
1288 
1289 
1290  // Construct triSurface. Note that addressing will be compact since
1291  // edgeSurface is compact.
1293  (
1294  triSurface
1295  (
1296  newTris,
1297  surf.patches(),
1298  eSurf.points()
1299  )
1300  );
1301 
1302  // Construct mapping back into original surface
1303  faceMap_.setSize(size());
1304 
1305  for (label faceI = 0; faceI < surf.size()-1; faceI++)
1306  {
1307  for (label triI = startTriI[faceI]; triI < startTriI[faceI+1]; triI++)
1308  {
1309  faceMap_[triI] = faceI;
1310  }
1311  }
1312  for (label triI = startTriI[surf.size()-1]; triI < size(); triI++)
1313  {
1314  faceMap_[triI] = surf.size()-1;
1315  }
1316 
1317 
1318  // Find edges on *this which originate from 'cuts'. (i.e. newEdgeI >=
1319  // nSurfaceEdges) Renumber edges into local triSurface numbering.
1320 
1321  intersectionEdges_.setSize(eSurf.edges().size() - eSurf.nSurfaceEdges());
1322 
1323  label intersectionEdgeI = 0;
1324 
1325  for
1326  (
1327  label edgeI = eSurf.nSurfaceEdges();
1328  edgeI < eSurf.edges().size();
1329  edgeI++
1330  )
1331  {
1332  // Get edge vertices in triSurface local numbering
1333  const edge& e = eSurf.edges()[edgeI];
1334  label surfStartI = meshPointMap()[e.start()];
1335  label surfEndI = meshPointMap()[e.end()];
1336 
1337  // Find edge connected to surfStartI which also uses surfEndI.
1338  label surfEdgeI = -1;
1339 
1340  const labelList& pEdges = pointEdges()[surfStartI];
1341 
1342  forAll(pEdges, i)
1343  {
1344  const edge& surfE = edges()[pEdges[i]];
1345 
1346  // Edge already connected to surfStart for sure. See if also
1347  // connects to surfEnd
1348  if (surfE.start() == surfEndI || surfE.end() == surfEndI)
1349  {
1350  surfEdgeI = pEdges[i];
1351 
1352  break;
1353  }
1354  }
1355 
1356  if (surfEdgeI != -1)
1357  {
1358  intersectionEdges_[intersectionEdgeI++] = surfEdgeI;
1359  }
1360  else
1361  {
1362  FatalErrorIn("intersectedSurface::intersectedSurface")
1363  << "Cannot find edge among candidates " << pEdges
1364  << " which uses points " << surfStartI
1365  << " and " << surfEndI
1366  << abort(FatalError);
1367  }
1368  }
1369 }
1370 
1371 
1372 // ************************************************************************* //
const edgeList & edges() const
Definition: edgeSurface.H:147
#define SeriousErrorIn(functionName)
Report an error message using Foam::SeriousError.
const pointField & points() const
Definition: edgeSurface.H:137
vector point
Point is a vector.
Definition: point.H:41
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurence of given element and return index,.
Description of surface in form of &#39;cloud of edges&#39;.
Definition: edgeSurface.H:73
dimensioned< scalar > mag(const dimensioned< Type > &)
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
labelList f(nPoints)
Basic surface-surface intersection description. Constructed from two surfaces it creates a descriptio...
bool empty() const
Return true if the UList is empty (ie, size() is zero).
Definition: UListI.H:313
const Map< label > & meshPointMap() const
Mesh point map. Given the global point index find its.
A HashTable to objects of type <T> with a label key.
Definition: Map.H:49
#define forAllIter(Container, container, iter)
Definition: UList.H:440
Triangulated surface description with patch information.
Definition: triSurface.H:57
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
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
triSurface()
Construct null.
Definition: triSurface.C:611
const labelListList & faceEdges() const
From face to our edges_.
Definition: edgeSurface.H:183
label nSurfacePoints() const
Definition: edgeSurface.H:142
friend Ostream & operator(Ostream &, const UList< T > &)
intersectedSurface()
Construct null.
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
label nPoints() const
Return number of points supporting patch faces.
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
Namespace for OpenFOAM.
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:39
label n
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
Triangle with additional region number.
Definition: labelledTri.H:49
static const char nl
Definition: Ostream.H:260
void setSize(const label)
Reset size of List.
Definition: List.C:318
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
const geometricSurfacePatchList & patches() const
Definition: triSurface.H:306
static const label ENDTOSTART
#define forAll(list, i)
Definition: UList.H:421
const edgeList & edges() const
Return list of edges, address into LOCAL point list.
const pointField & cutPoints() const
const Cmpt & x() const
Definition: VectorI.H:65
Triangulation of faces. Handles concave polygons as well (inefficiently)
void operator=(const triSurface &)
Definition: triSurface.C:1125
errorManip< error > abort(error &err)
Definition: errorManip.H:131
const Field< PointType > & faceNormals() const
Return face normals for patch.
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
List< edge > edgeList
Definition: edgeList.H:38
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
Definition: DynamicList.H:56
error FatalError
const labelListList & pointEdges() const
Return point-edge addressing.
label end() const
Return end vertex label.
Definition: edgeI.H:92
List< label > labelList
A List of labels.
Definition: labelList.H:56
scalar pseudoAngle(const vector &e0, const vector &e1, const vector &vec)
Estimate angle of vec in coordinate system (e0, e1, e0^e1).
Definition: transform.H:223
HashTable< label, label, Hash< label > >::iterator iterator
Definition: Map.H:56
const labelListList & faceEdges() const
Return face-edge addressing.
3D tensor transformation operations.
static const label UNVISITED
label size() const
Return the number of elements in the UList.
static const label STARTTOEND
A triangular face using a FixedList of labels corresponding to mesh vertices.
Definition: triFace.H:68
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
void reverse(UList< T > &, const label n)
Definition: UListI.H:322
const List< Face > & localFaces() const
Return patch faces addressing into local point list.
label walkFace(const primitiveMesh &, const label faceI, const label startEdgeI, const label startVertI, const label nEdges)
Returns label of edge nEdges away from startEdge (in the direction.
Definition: meshTools.C:608
label start() const
Return start vertex label.
Definition: edgeI.H:81
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of point.
Definition: meshTools.C:209
defineTypeNameAndDebug(combustionModel, 0)
const edgeList & cutEdges() const
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53
label nSurfaceEdges() const
Definition: edgeSurface.H:152