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-2023 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 {
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  {
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 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
1145 (
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 // ************************************************************************* //
label n
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
#define forAllIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:459
#define forAllConstIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:477
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:78
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:296
DynamicList< T, SizeInc, SizeMult, SizeDiv > & shrink()
Shrink the allocated space to the number of elements used.
Definition: DynamicListI.H:252
friend class iterator
Declare friendship with the iterator.
Definition: HashTable.H:194
label size() const
Return the number of elements in the UList.
Definition: ListI.H:171
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
void setSize(const label)
Reset size of List.
Definition: List.C:281
A HashTable to objects of type <T> with a label key.
Definition: Map.H:52
const labelListList & pointEdges() const
Return point-edge addressing.
const edgeList & edges() const
Return list of edges, address into LOCAL point list.
const Map< label > & meshPointMap() const
Mesh point map. Given the global point index find its.
const Field< PointType > & points() const
Return reference to global points.
const List< FaceType > & localFaces() const
Return patch faces addressing into local point list.
const labelListList & faceEdges() const
Return face-edge addressing.
const Field< PointType > & faceNormals() const
Return face normals for patch.
A List with indirect addressing.
Definition: UIndirectList.H:60
bool empty() const
Return true if the UList is empty (ie, size() is zero)
Definition: UListI.H:325
friend Ostream & operator(Ostream &, const UList< T > &)
const Cmpt & x() const
Definition: VectorI.H:75
Description of surface in form of 'cloud of edges'.
Definition: edgeSurface.H:74
label nSurfaceEdges() const
Definition: edgeSurface.H:150
const labelListList & faceEdges() const
From face to our edges_.
Definition: edgeSurface.H:179
label nSurfacePoints() const
Definition: edgeSurface.H:140
const edgeList & edges() const
Definition: edgeSurface.H:145
const pointField & points() const
Definition: edgeSurface.H:135
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:61
label end() const
Return end vertex label.
Definition: edgeI.H:92
label start() const
Return start vertex label.
Definition: edgeI.H:81
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:76
Given triSurface and intersection creates the intersected (properly triangulated) surface....
intersectedSurface()
Construct null.
static const label ENDTOSTART
static const label UNVISITED
static const label STARTTOEND
Triangle with additional region number.
Definition: labelledTri.H:60
Triangulation of three-dimensional polygons.
const UList< triFace > & triPoints() const
Get the triangles' points.
Basic surface-surface intersection description. Constructed from two surfaces it creates a descriptio...
const edgeList & cutEdges() const
const pointField & cutPoints() const
Triangulated surface description with patch information.
Definition: triSurface.H:69
triSurface()
Construct null.
Definition: triSurface.C:535
void operator=(const triSurface &)
Definition: triSurface.C:1386
const geometricSurfacePatchList & patches() const
Definition: triSurface.H:322
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:334
const pointField & points
label nPoints
#define SeriousErrorInFunction
Report an error message using Foam::SeriousError.
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
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of point.
Definition: meshTools.C:203
Namespace for OpenFOAM.
const doubleScalar e
Definition: doubleScalar.H:106
scalar minDist(const List< pointIndexHit > &hitList)
List< label > labelList
A List of labels.
Definition: labelList.H:56
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
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:257
word name(const bool)
Return a word representation of a bool.
Definition: boolIO.C:39
errorManip< error > abort(error &err)
Definition: errorManip.H:131
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
vector point
Point is a vector.
Definition: point.H:41
void reverse(UList< T > &, const label n)
Definition: UListI.H:334
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
dimensioned< scalar > mag(const dimensioned< Type > &)
defineTypeNameAndDebug(combustionModel, 0)
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
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
error FatalError
List< edge > edgeList
Definition: edgeList.H:38
void offset(label &lst, const label o)
static const char nl
Definition: Ostream.H:266
labelList f(nPoints)
3D tensor transformation operations.