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