intersectedSurface.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "intersectedSurface.H"
27 #include "surfaceIntersection.H"
28 #include "faceList.H"
29 #include "polygonTriangulate.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 
46 const Foam::label Foam::intersectedSurface::UNVISITED = 0;
47 const Foam::label Foam::intersectedSurface::STARTTOEND = 1;
48 const Foam::label Foam::intersectedSurface::ENDTOSTART = 2;
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  {
246  << "Triangle:" << fA << " and triangle:" << fB
247  << " share a point but not an edge"
248  << abort(FatalError);
249  }
250  }
251  }
252 
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  << "Point:" << iter.key() << " used by too few edges:"
346  << iter() << abort(FatalError);
347  }
348  }
349 
350  if (debug & 2)
351  {
352  // Print facePointEdges
353  Pout<< "calcPointEdgeAddressing: face consisting of edges:" << endl;
354  forAll(fEdges, i)
355  {
356  label edgeI = fEdges[i];
357  const edge& e = edges[edgeI];
358  Pout<< " " << edgeI << ' ' << e
359  << points[e.start()]
360  << points[e.end()] << endl;
361  }
362 
363  Pout<< " Constructed point-edge addressing:" << endl;
364  forAllConstIter(Map<DynamicList<label>>, facePointEdges, iter)
365  {
366  Pout<< " vertex " << iter.key() << " is connected to edges "
367  << iter() << endl;
368  }
369  Pout<<endl;
370  }
371 
372  return facePointEdges;
373 }
374 
375 
376 // Find next (triangle or cut) edge label coming from point prevVertI on
377 // prevEdgeI doing a right handed walk (i.e. following right wall).
378 // Note: normal is provided externally. Could be deducted from angle between
379 // two triangle edges but these could be in line.
380 Foam::label Foam::intersectedSurface::nextEdge
381 (
382  const edgeSurface& eSurf,
383  const Map<label>& visited,
384  const label facei,
385  const vector& n,
386  const Map<DynamicList<label>>& facePointEdges,
387  const label prevEdgeI,
388  const label prevVertI
389 )
390 {
391  const pointField& points = eSurf.points();
392  const edgeList& edges = eSurf.edges();
393  const labelList& fEdges = eSurf.faceEdges()[facei];
394 
395 
396  // Edges connected to prevVertI
397  const DynamicList<label>& connectedEdges = facePointEdges[prevVertI];
398 
399  if (connectedEdges.size() <= 1)
400  {
401  // Problem. Point not connected.
402  {
403  Pout<< "Writing face:" << facei << " to face.obj" << endl;
404  OFstream str("face.obj");
405  writeOBJ(points, edges, fEdges, str);
406 
407  Pout<< "Writing connectedEdges edges to faceEdges.obj" << endl;
408  writeLocalOBJ(points, edges, connectedEdges, "faceEdges.obj");
409  }
410 
412  << "Problem: prevVertI:" << prevVertI << " on edge " << prevEdgeI
413  << " has less than 2 connected edges."
414  << " connectedEdges:" << connectedEdges << abort(FatalError);
415 
416  return -1;
417  }
418 
419  if (connectedEdges.size() == 2)
420  {
421  // Simple case. Take other edge
422  if (connectedEdges[0] == prevEdgeI)
423  {
424  if (debug & 2)
425  {
426  Pout<< "Stepped from edge:" << edges[prevEdgeI]
427  << " to edge:" << edges[connectedEdges[1]]
428  << " chosen from candidates:" << connectedEdges << endl;
429  }
430  return connectedEdges[1];
431  }
432  else
433  {
434  if (debug & 2)
435  {
436  Pout<< "Stepped from edge:" << edges[prevEdgeI]
437  << " to edge:" << edges[connectedEdges[0]]
438  << " chosen from candidates:" << connectedEdges << endl;
439  }
440  return connectedEdges[0];
441  }
442  }
443 
444 
445  // Multiple choices. Look at angle between edges.
446 
447  const edge& prevE = edges[prevEdgeI];
448 
449  // x-axis of coordinate system
450  vector e0 = n ^ (points[prevE.otherVertex(prevVertI)] - points[prevVertI]);
451  e0 /= mag(e0) + vSmall;
452 
453  // Get y-axis of coordinate system
454  vector e1 = n ^ e0;
455 
456  if (mag(mag(e1) - 1) > small)
457  {
458  {
459  Pout<< "Writing face:" << facei << " to face.obj" << endl;
460  OFstream str("face.obj");
461  writeOBJ(points, edges, fEdges, str);
462 
463  Pout<< "Writing connectedEdges edges to faceEdges.obj" << endl;
464  writeLocalOBJ(points, edges, connectedEdges, "faceEdges.obj");
465  }
466 
468  << "Unnormalised normal e1:" << e1
469  << " formed from cross product of e0:" << e0 << " n:" << n
470  << abort(FatalError);
471  }
472 
473 
474  //
475  // Determine maximum angle over all connected (unvisited) edges.
476  //
477 
478  scalar maxAngle = -great;
479  label maxEdgeI = -1;
480 
481  forAll(connectedEdges, connI)
482  {
483  label edgeI = connectedEdges[connI];
484 
485  if (edgeI != prevEdgeI)
486  {
487  label stat = visited[edgeI];
488 
489  const edge& e = edges[edgeI];
490 
491  // Find out whether walk of edge from prevVert would be acceptable.
492  if
493  (
494  stat == UNVISITED
495  || (
496  stat == STARTTOEND
497  && prevVertI == e[1]
498  )
499  || (
500  stat == ENDTOSTART
501  && prevVertI == e[0]
502  )
503  )
504  {
505  // Calculate angle of edge with respect to base e0, e1
506  vector vec =
507  n ^ (points[e.otherVertex(prevVertI)] - points[prevVertI]);
508 
509  vec /= mag(vec) + vSmall;
510 
511  scalar angle = pseudoAngle(e0, e1, vec);
512 
513  if (angle > maxAngle)
514  {
515  maxAngle = angle;
516  maxEdgeI = edgeI;
517  }
518  }
519  }
520  }
521 
522 
523  if (maxEdgeI == -1)
524  {
525  // No unvisited edge found
526  {
527  Pout<< "Writing face:" << facei << " to face.obj" << endl;
528  OFstream str("face.obj");
529  writeOBJ(points, edges, fEdges, str);
530 
531  Pout<< "Writing connectedEdges edges to faceEdges.obj" << endl;
532  writeLocalOBJ(points, edges, connectedEdges, "faceEdges.obj");
533  }
534 
536  << "Trying to step from edge " << edges[prevEdgeI]
537  << ", vertex " << prevVertI
538  << " but cannot find 'unvisited' edges among candidates:"
539  << connectedEdges
540  << abort(FatalError);
541  }
542 
543  if (debug & 2)
544  {
545  Pout<< "Stepped from edge:" << edges[prevEdgeI]
546  << " to edge:" << maxEdgeI << " angle:" << edges[maxEdgeI]
547  << " with angle:" << maxAngle
548  << endl;
549  }
550 
551  return maxEdgeI;
552 }
553 
554 
555 // Create (polygonal) face by walking full circle starting from startVertI on
556 // startEdgeI.
557 // Uses nextEdge(..) to do the walking. Returns face. Updates visited with
558 // the labels of visited edges.
559 Foam::face Foam::intersectedSurface::walkFace
560 (
561  const edgeSurface& eSurf,
562  const label facei,
563  const vector& n,
564  const Map<DynamicList<label>>& facePointEdges,
565 
566  const label startEdgeI,
567  const label startVertI,
568 
569  Map<label>& visited
570 )
571 {
572  const pointField& points = eSurf.points();
573  const edgeList& edges = eSurf.edges();
574 
575  // Overestimate size of face
576  face f(eSurf.faceEdges()[facei].size());
577 
578  label fp = 0;
579 
580  label vertI = startVertI;
581  label edgeI = startEdgeI;
582 
583  while (true)
584  {
585  const edge& e = edges[edgeI];
586 
587  if (debug & 2)
588  {
589  Pout<< "Now at:" << endl
590  << " edge:" << edgeI << " vertices:" << e
591  << " positions:" << points[e.start()] << ' ' << points[e.end()]
592  << " vertex:" << vertI << endl;
593  }
594 
595  // Mark edge as visited
596  if (e[0] == vertI)
597  {
598  visited[edgeI] |= STARTTOEND;
599  }
600  else
601  {
602  visited[edgeI] |= ENDTOSTART;
603  }
604 
605 
606  // Store face vertex
607  f[fp++] = vertI;
608 
609  // step to other vertex
610  vertI = e.otherVertex(vertI);
611 
612  if (vertI == startVertI)
613  {
614  break;
615  }
616 
617  // step from vertex to next edge
618  edgeI = nextEdge
619  (
620  eSurf,
621  visited,
622  facei,
623  n,
624  facePointEdges,
625  edgeI,
626  vertI
627  );
628  }
629 
630  f.setSize(fp);
631 
632  return f;
633 }
634 
635 
636 void Foam::intersectedSurface::findNearestVisited
637 (
638  const edgeSurface& eSurf,
639  const label facei,
640  const Map<DynamicList<label>>& facePointEdges,
641  const Map<label>& pointVisited,
642  const point& pt,
643  const label excludePointi,
644 
645  label& minVertI,
646  scalar& minDist
647 )
648 {
649  minVertI = -1;
650  minDist = great;
651 
652  forAllConstIter(Map<label>, pointVisited, iter)
653  {
654  label pointi = iter.key();
655 
656  if (pointi != excludePointi)
657  {
658  label nVisits = iter();
659 
660  if (nVisits == 2*facePointEdges[pointi].size())
661  {
662  // Fully visited (i.e. both sides of all edges)
663  scalar dist = mag(eSurf.points()[pointi] - pt);
664 
665  if (dist < minDist)
666  {
667  minDist = dist;
668  minVertI = pointi;
669  }
670  }
671  }
672  }
673 
674  if (minVertI == -1)
675  {
676  const labelList& fEdges = eSurf.faceEdges()[facei];
677 
679  << "Dumping face edges to faceEdges.obj" << endl;
680 
681  writeLocalOBJ(eSurf.points(), eSurf.edges(), fEdges, "faceEdges.obj");
682 
684  << "No fully visited edge found for pt " << pt
685  << abort(FatalError);
686  }
687 }
688 
689 
690 // Sometimes there are bunches of edges that are not connected to the
691 // other edges in the face. This routine tries to connect the loose
692 // edges up to the 'proper' edges by adding two extra edges between a
693 // properly connected edge and an unconnected one. Since at this level the
694 // addressing is purely in form of points and a cloud of edges this can
695 // be easily done.
696 // (edges are to existing points. Don't want to introduce new vertices here
697 // since then also neighbouring face would have to be split)
698 Foam::faceList Foam::intersectedSurface::resplitFace
699 (
700  const triSurface& surf,
701  const label facei,
702  const Map<DynamicList<label>>& facePointEdges,
703  const Map<label>& visited,
704  edgeSurface& eSurf
705 )
706 {
707  {
708  // Dump face for debugging.
709  Pout<< "Writing face:" << facei << " to face.obj" << endl;
710  OFstream str("face.obj");
711  writeOBJ(eSurf.points(), eSurf.edges(), eSurf.faceEdges()[facei], str);
712  }
713 
714 
715  // Count the number of times point has been visited so we
716  // can compare number to facePointEdges.
717  Map<label> pointVisited(2*facePointEdges.size());
718 
719 
720  {
721  OFstream str0("visitedNone.obj");
722  OFstream str1("visitedOnce.obj");
723  OFstream str2("visitedTwice.obj");
724  forAll(eSurf.points(), pointi)
725  {
726  const point& pt = eSurf.points()[pointi];
727 
728  str0 << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << nl;
729  str1 << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << nl;
730  str2 << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << nl;
731  }
732 
733 
734  forAllConstIter(Map<label>, visited, iter)
735  {
736  label edgeI = iter.key();
737 
738  const edge& e = eSurf.edges()[edgeI];
739 
740  label stat = iter();
741 
742  if (stat == STARTTOEND || stat == ENDTOSTART)
743  {
744  incCount(pointVisited, e[0], 1);
745  incCount(pointVisited, e[1], 1);
746 
747  str1 << "l " << e[0]+1 << ' ' << e[1]+1 << nl;
748  }
749  else if (stat == BOTH)
750  {
751  incCount(pointVisited, e[0], 2);
752  incCount(pointVisited, e[1], 2);
753 
754  str2 << "l " << e[0]+1 << ' ' << e[1]+1 << nl;
755  }
756  else if (stat == UNVISITED)
757  {
758  incCount(pointVisited, e[0], 0);
759  incCount(pointVisited, e[1], 0);
760 
761  str0 << "l " << e[0]+1 << ' ' << e[1]+1 << nl;
762  }
763  }
764  }
765 
766 
767  {
768  forAllConstIter(Map<label>, pointVisited, iter)
769  {
770  label pointi = iter.key();
771 
772  label nVisits = iter();
773 
774  Pout<< "point:" << pointi << " nVisited:" << nVisits
775  << " pointEdges:" << facePointEdges[pointi].size() << endl;
776  }
777  }
778 
779 
780  // Find nearest point pair where one is not fully visited and
781  // the other is.
782 
783  label visitedVert0 = -1;
784  label unvisitedVert0 = -1;
785 
786  {
787  scalar minDist = great;
788 
789  forAllConstIter(Map<DynamicList<label>>, facePointEdges, iter)
790  {
791  label pointi = iter.key();
792 
793  label nVisits = pointVisited[pointi];
794 
795  const DynamicList<label>& pEdges = iter();
796 
797  if (nVisits < 2*pEdges.size())
798  {
799  // Not fully visited. Find nearest fully visited.
800 
801  scalar nearDist;
802  label nearVertI;
803 
804  findNearestVisited
805  (
806  eSurf,
807  facei,
808  facePointEdges,
809  pointVisited,
810  eSurf.points()[pointi],
811  -1, // Do not exclude vertex
812  nearVertI,
813  nearDist
814  );
815 
816 
817  if (nearDist < minDist)
818  {
819  minDist = nearDist;
820  visitedVert0 = nearVertI;
821  unvisitedVert0 = pointi;
822  }
823  }
824  }
825  }
826 
827 
828  // Find second intersection.
829  label visitedVert1 = -1;
830  label unvisitedVert1 = -1;
831  {
832  scalar minDist = great;
833 
834  forAllConstIter(Map<DynamicList<label>>, facePointEdges, iter)
835  {
836  label pointi = iter.key();
837 
838  if (pointi != unvisitedVert0)
839  {
840  label nVisits = pointVisited[pointi];
841 
842  const DynamicList<label>& pEdges = iter();
843 
844  if (nVisits < 2*pEdges.size())
845  {
846  // Not fully visited. Find nearest fully visited.
847 
848  scalar nearDist;
849  label nearVertI;
850 
851  findNearestVisited
852  (
853  eSurf,
854  facei,
855  facePointEdges,
856  pointVisited,
857  eSurf.points()[pointi],
858  visitedVert0, // vertex to exclude
859  nearVertI,
860  nearDist
861  );
862 
863 
864  if (nearDist < minDist)
865  {
866  minDist = nearDist;
867  visitedVert1 = nearVertI;
868  unvisitedVert1 = pointi;
869  }
870  }
871  }
872  }
873  }
874 
875 
876  Pout<< "resplitFace : adding intersection from " << visitedVert0
877  << " to " << unvisitedVert0 << endl
878  << " and from " << visitedVert1
879  << " to " << unvisitedVert1 << endl;
880 
881 
882 // // Add the new intersection edges to the edgeSurface
883 // edgeList additionalEdges(2);
884 // additionalEdges[0] = edge(visitedVert0, unvisitedVert0);
885 // additionalEdges[1] = edge(visitedVert1, unvisitedVert1);
886 
887  // Add the new intersection edges to the edgeSurface
888  edgeList additionalEdges(1);
889  additionalEdges[0] = edge(visitedVert0, unvisitedVert0);
890 
891  eSurf.addIntersectionEdges(facei, additionalEdges);
892 
893  fileName newFName("face_" + Foam::name(facei) + "_newEdges.obj");
894  Pout<< "Dumping face:" << facei << " to " << newFName << endl;
895  writeLocalOBJ
896  (
897  eSurf.points(),
898  eSurf.edges(),
899  eSurf.faceEdges()[facei],
900  newFName
901  );
902 
903  // Retry splitFace. Use recursion since is rare situation.
904  return splitFace(surf, facei, eSurf);
905 }
906 
907 
908 Foam::faceList Foam::intersectedSurface::splitFace
909 (
910  const triSurface& surf,
911  const label facei,
912  edgeSurface& eSurf
913 )
914 {
915  // Alias
916  const pointField& points = eSurf.points();
917  const edgeList& edges = eSurf.edges();
918  const labelList& fEdges = eSurf.faceEdges()[facei];
919 
920  // Create local (for the face only) point-edge connectivity.
921  Map<DynamicList<label>> facePointEdges
922  (
923  calcPointEdgeAddressing
924  (
925  eSurf,
926  facei
927  )
928  );
929 
930  // Order in which edges have been walked. Initialise outside edges.
931  Map<label> visited(fEdges.size()*2);
932 
933  forAll(fEdges, i)
934  {
935  label edgeI = fEdges[i];
936 
937  if (eSurf.isSurfaceEdge(edgeI))
938  {
939  // Edge is edge from original surface so an outside edge for
940  // the current face.
941  label surfEdgeI = eSurf.parentEdge(edgeI);
942 
943  label owner = surf.edgeOwner()[surfEdgeI];
944 
945  if
946  (
947  owner == facei
948  || sameEdgeOrder
949  (
950  surf.localFaces()[owner],
951  surf.localFaces()[facei]
952  )
953  )
954  {
955  // Edge is in same order as current face.
956  // Mark off the opposite order.
957  visited.insert(edgeI, ENDTOSTART);
958  }
959  else
960  {
961  // Edge is in same order as current face.
962  // Mark off the opposite order.
963  visited.insert(edgeI, STARTTOEND);
964  }
965  }
966  else
967  {
968  visited.insert(edgeI, UNVISITED);
969  }
970  }
971 
972 
973 
974  // Storage for faces
975  DynamicList<face> faces(fEdges.size());
976 
977  while (true)
978  {
979  // Find starting edge:
980  // - unvisted triangle edge
981  // - once visited intersection edge
982  // Give priority to triangle edges.
983  label startEdgeI = -1;
984  label startVertI = -1;
985 
986  forAll(fEdges, i)
987  {
988  label edgeI = fEdges[i];
989 
990  const edge& e = edges[edgeI];
991 
992  label stat = visited[edgeI];
993 
994  if (stat == STARTTOEND)
995  {
996  startEdgeI = edgeI;
997  startVertI = e[1];
998 
999  if (eSurf.isSurfaceEdge(edgeI))
1000  {
1001  break;
1002  }
1003  }
1004  else if (stat == ENDTOSTART)
1005  {
1006  startEdgeI = edgeI;
1007  startVertI = e[0];
1008 
1009  if (eSurf.isSurfaceEdge(edgeI))
1010  {
1011  break;
1012  }
1013  }
1014  }
1015 
1016  if (startEdgeI == -1)
1017  {
1018  break;
1019  }
1020 
1021  // Pout<< "splitFace: starting walk from edge:" << startEdgeI
1022  // << ' ' << edges[startEdgeI] << " vertex:" << startVertI << endl;
1023 
1025  // printVisit(eSurf.edges(), fEdges, visited);
1026 
1027  //{
1028  // Pout<< "Writing face:" << facei << " to face.obj" << endl;
1029  // OFstream str("face.obj");
1030  // writeOBJ(eSurf.points(), eSurf.edges(), fEdges, str);
1031  //}
1032 
1033  faces.append
1034  (
1035  walkFace
1036  (
1037  eSurf,
1038  facei,
1039  surf.faceNormals()[facei],
1040  facePointEdges,
1041 
1042  startEdgeI,
1043  startVertI,
1044 
1045  visited
1046  )
1047  );
1048  }
1049 
1050 
1051  // Check if any unvisited edges left.
1052  forAll(fEdges, i)
1053  {
1054  label edgeI = fEdges[i];
1055 
1056  label stat = visited[edgeI];
1057 
1058  if (eSurf.isSurfaceEdge(edgeI) && stat != BOTH)
1059  {
1061  << "Dumping face edges to faceEdges.obj" << endl;
1062 
1063  writeLocalOBJ(points, edges, fEdges, "faceEdges.obj");
1064 
1066  << "Problem: edge " << edgeI << " vertices "
1067  << edges[edgeI] << " on face " << facei
1068  << " has visited status " << stat << " from a "
1069  << "righthanded walk along all"
1070  << " of the triangle edges. Are the original surfaces"
1071  << " closed and non-intersecting?"
1072  << abort(FatalError);
1073  }
1074  else if (stat != BOTH)
1075  {
1076  //{
1077  // Pout<< "Dumping faces so far to faces.obj" << nl
1078  // << faces << endl;
1079  //
1080  // OFstream str("faces.obj");
1081  //
1082  // forAll(faces, i)
1083  // {
1084  // writeOBJ(points, faces[i], str);
1085  // }
1086  //}
1087 
1088  Pout<< "** Resplitting **" << endl;
1089 
1090  // Redo face after introducing extra edge. Edge introduced
1091  // should be one nearest to any fully visited edge.
1092  return resplitFace
1093  (
1094  surf,
1095  facei,
1096  facePointEdges,
1097  visited,
1098  eSurf
1099  );
1100  }
1101  }
1102 
1103 
1104  // See if normal needs flipping.
1105  faces.shrink();
1106 
1107  vector a = faces[0].area(eSurf.points());
1108 
1109  if ((a & surf.faceNormals()[facei]) < 0)
1110  {
1111  forAll(faces, i)
1112  {
1113  reverse(faces[i]);
1114  }
1115  }
1116 
1117  return move(faces);
1118 }
1119 
1120 
1121 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
1122 
1123 // Null constructor
1125 :
1126  triSurface(),
1127  intersectionEdges_(0),
1128  faceMap_(0),
1129  nSurfacePoints_(0)
1130 {}
1131 
1132 
1133 // Construct from components
1135 :
1136  triSurface(surf),
1137  intersectionEdges_(0),
1138  faceMap_(0),
1139  nSurfacePoints_(0)
1140 {}
1141 
1142 
1143 // Construct from surface and intersection
1146  const triSurface& surf,
1147  const bool isFirstSurface,
1148  const surfaceIntersection& inter
1149 )
1150 :
1151  triSurface(),
1152  intersectionEdges_(0),
1153  faceMap_(0),
1154  nSurfacePoints_(surf.nPoints())
1155 {
1156  if (inter.cutPoints().empty() && inter.cutEdges().empty())
1157  {
1158  // No intersection. Make straight copy.
1159  triSurface::operator=(surf);
1160 
1161  // Identity for face map
1162  faceMap_.setSize(size());
1163 
1164  forAll(faceMap_, facei)
1165  {
1166  faceMap_[facei] = facei;
1167  }
1168  return;
1169  }
1170 
1171 
1172  // Calculate face-edge addressing for surface + intersection.
1173  edgeSurface eSurf(surf, isFirstSurface, inter);
1174 
1175  // Update point information for any removed points. (when are they removed?
1176  // - but better make sure)
1177  nSurfacePoints_ = eSurf.nSurfacePoints();
1178 
1179  // Now we have full points, edges and edge addressing for surf. Start
1180  // extracting faces and triangulate them.
1181 
1182  // Storage for new triangles of surface.
1183  DynamicList<labelledTri> newTris(eSurf.edges().size()/2);
1184 
1185  // Start in newTris for decomposed face.
1186  labelList startTriI(surf.size(), 0);
1187 
1188  // Create a triangulation engine
1189  polygonTriangulate triEngine;
1190 
1191  forAll(surf, facei)
1192  {
1193  startTriI[facei] = newTris.size();
1194 
1195  if (eSurf.faceEdges()[facei].size() != surf.faceEdges()[facei].size())
1196  {
1197  // Face has been cut by intersection.
1198  // Cut face into multiple subfaces. Use faceEdge information
1199  // from edgeSurface only. (original triSurface 'surf' is used for
1200  // faceNormals and regionnumber only)
1201  faceList newFaces
1202  (
1203  splitFace
1204  (
1205  surf,
1206  facei, // current triangle
1207  eSurf // face-edge description of surface
1208  // + intersection
1209  )
1210  );
1211  forAll(newFaces, newFacei)
1212  {
1213  const face& newF = newFaces[newFacei];
1214 
1215  const vector& n = surf.faceNormals()[facei];
1216  const label region = surf[facei].region();
1217 
1218  triEngine.triangulate
1219  (
1220  UIndirectList<point>(eSurf.points(), newF),
1221  n
1222  );
1223 
1224  forAll(triEngine.triPoints(), trii)
1225  {
1226  newTris.append
1227  (
1228  labelledTri(triEngine.triPoints(trii, newF), region)
1229  );
1230  }
1231  }
1232  }
1233  else
1234  {
1235  // Face has not been cut at all. No need to renumber vertices since
1236  // eSurf keeps surface vertices first.
1237  newTris.append(surf.localFaces()[facei]);
1238  }
1239  }
1240 
1241  newTris.shrink();
1242 
1243 
1244  // Construct triSurface. Note that addressing will be compact since
1245  // edgeSurface is compact.
1247  (
1248  triSurface
1249  (
1250  newTris,
1251  surf.patches(),
1252  eSurf.points()
1253  )
1254  );
1255 
1256  // Construct mapping back into original surface
1257  faceMap_.setSize(size());
1258 
1259  for (label facei = 0; facei < surf.size()-1; facei++)
1260  {
1261  for (label triI = startTriI[facei]; triI < startTriI[facei+1]; triI++)
1262  {
1263  faceMap_[triI] = facei;
1264  }
1265  }
1266  for (label triI = startTriI[surf.size()-1]; triI < size(); triI++)
1267  {
1268  faceMap_[triI] = surf.size()-1;
1269  }
1270 
1271 
1272  // Find edges on *this which originate from 'cuts'. (i.e. newEdgeI >=
1273  // nSurfaceEdges) Renumber edges into local triSurface numbering.
1274 
1275  intersectionEdges_.setSize(eSurf.edges().size() - eSurf.nSurfaceEdges());
1276 
1277  label intersectionEdgeI = 0;
1278 
1279  for
1280  (
1281  label edgeI = eSurf.nSurfaceEdges();
1282  edgeI < eSurf.edges().size();
1283  edgeI++
1284  )
1285  {
1286  // Get edge vertices in triSurface local numbering
1287  const edge& e = eSurf.edges()[edgeI];
1288  label surfStartI = meshPointMap()[e.start()];
1289  label surfEndI = meshPointMap()[e.end()];
1290 
1291  // Find edge connected to surfStartI which also uses surfEndI.
1292  label surfEdgeI = -1;
1293 
1294  const labelList& pEdges = pointEdges()[surfStartI];
1295 
1296  forAll(pEdges, i)
1297  {
1298  const edge& surfE = edges()[pEdges[i]];
1299 
1300  // Edge already connected to surfStart for sure. See if also
1301  // connects to surfEnd
1302  if (surfE.start() == surfEndI || surfE.end() == surfEndI)
1303  {
1304  surfEdgeI = pEdges[i];
1305 
1306  break;
1307  }
1308  }
1309 
1310  if (surfEdgeI != -1)
1311  {
1312  intersectionEdges_[intersectionEdgeI++] = surfEdgeI;
1313  }
1314  else
1315  {
1317  << "Cannot find edge among candidates " << pEdges
1318  << " which uses points " << surfStartI
1319  << " and " << surfEndI
1320  << abort(FatalError);
1321  }
1322  }
1323 }
1324 
1325 
1326 // ************************************************************************* //
const labelListList & pointEdges() const
Return point-edge addressing.
label nPoints() const
Return number of points supporting patch faces.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
bool empty() const
Return true if the UList is empty (ie, size() is zero)
Definition: UListI.H:325
label nSurfacePoints() const
Definition: edgeSurface.H:140
FvWallInfoData< WallInfo, label > label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
friend Ostream & operator(Ostream &, const UList< T > &)
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
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:306
error FatalError
void operator=(const triSurface &)
Definition: triSurface.C:1385
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
#define forAllConstIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:477
label nSurfaceEdges() const
Definition: edgeSurface.H:150
HashTable< label, label, Hash< label > >::iterator iterator
Definition: Map.H:56
#define forAllIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:459
Basic surface-surface intersection description. Constructed from two surfaces it creates a descriptio...
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
intersectedSurface()
Construct null.
const pointField & points() const
Definition: edgeSurface.H:135
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
Description of surface in form of &#39;cloud of edges&#39;.
Definition: edgeSurface.H:73
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of point.
Definition: meshTools.C:203
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
#define SeriousErrorInFunction
Report an error message using Foam::SeriousError.
const edgeList & cutEdges() const
const Map< label > & meshPointMap() const
Mesh point map. Given the global point index find its.
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
List< edge > edgeList
Definition: edgeList.H:38
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
3D tensor transformation operations.
const edgeList & edges() const
Definition: edgeSurface.H:145
const edgeList & edges() const
Return list of edges, address into LOCAL point list.
List< label > labelList
A List of labels.
Definition: labelList.H:56
Triangle with additional region number.
Definition: labelledTri.H:57
const Field< PointType > & faceNormals() const
Return face normals for patch.
errorManip< error > abort(error &err)
Definition: errorManip.H:131
const Cmpt & x() const
Definition: VectorI.H:75
void reverse(UList< T > &, const label n)
Definition: UListI.H:334
static const char nl
Definition: Ostream.H:260
defineTypeNameAndDebug(combustionModel, 0)
labelList f(nPoints)
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
const pointField & cutPoints() const
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
void setSize(const label)
Reset size of List.
Definition: List.C:281
const geometricSurfacePatchList & patches() const
Definition: triSurface.H:322
const labelListList & faceEdges() const
From face to our edges_.
Definition: edgeSurface.H:179
triSurface()
Construct null.
Definition: triSurface.C:534
vector point
Point is a vector.
Definition: point.H:41
Triangulation of three-dimensional polygons.
label end() const
Return end vertex label.
Definition: edgeI.H:92
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
const List< FaceType > & localFaces() const
Return patch faces addressing into local point list.
A List with indirect addressing.
Definition: fvMatrix.H:106
const labelListList & faceEdges() const
Return face-edge addressing.
dimensioned< scalar > mag(const dimensioned< Type > &)
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:587
Triangulated surface description with patch information.
Definition: triSurface.H:66
label size() const
Return the number of elements in the UList.
Definition: ListI.H:171
label start() const
Return start vertex label.
Definition: edgeI.H:81
Namespace for OpenFOAM.
A HashTable to objects of type <T> with a label key.
Definition: Map.H:49