FreeFOAM The Cross-Platform CFD Toolkit
booleanSurface.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 Description
25 
26 \*---------------------------------------------------------------------------*/
27 
28 #include "booleanSurface.H"
32 #include <OpenFOAM/OFstream.H>
33 #include <meshTools/treeBoundBox.H>
34 #include <meshTools/meshTools.H>
35 
36 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
37 
39 
40 
41 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
42 
43 // Check whether at least one of faces connected to the intersection has been
44 // marked.
45 void Foam::booleanSurface::checkIncluded
46 (
47  const intersectedSurface& surf,
48  const labelList& faceZone,
49  const label includedFace
50 )
51 {
52  forAll(surf.intersectionEdges(), intEdgeI)
53  {
54  label edgeI = surf.intersectionEdges()[intEdgeI];
55 
56  const labelList& myFaces = surf.edgeFaces()[edgeI];
57 
58  bool usesIncluded = false;
59 
60  forAll(myFaces, myFaceI)
61  {
62  if (faceZone[myFaces[myFaceI]] == faceZone[includedFace])
63  {
64  usesIncluded = true;
65 
66  break;
67  }
68  }
69 
70  if (!usesIncluded)
71  {
73  (
74  "booleanSurface::checkIncluded(const intersectedSurface&"
75  ", const labelList&, const label)"
76  ) << "None of the faces reachable from face " << includedFace
77  << " connects to the intersection."
78  << exit(FatalError);
79  }
80  }
81 }
82 
83 
84 // Linear lookup
85 Foam::label Foam::booleanSurface::index
86 (
87  const labelList& elems,
88  const label elem
89 )
90 {
91  forAll(elems, elemI)
92  {
93  if (elems[elemI] == elem)
94  {
95  return elemI;
96  }
97  }
98  return -1;
99 }
100 
101 
103 (
104  const edgeList& edges,
105  const labelList& edgeLabels,
106  const edge& e
107 )
108 {
109  forAll(edgeLabels, edgeLabelI)
110  {
111  if (edges[edgeLabels[edgeLabelI]] == e)
112  {
113  return edgeLabels[edgeLabelI];
114  }
115  }
117  (
118  "booleanSurface::findEdge(const edgeList&, const labelList&"
119  ", const edge&)"
120  ) << "Cannot find edge " << e << " in edges " << edgeLabels
121  << abort(FatalError);
122 
123  return -1;
124 }
125 
126 
127 // Generate combined patchList (returned). Sets patchMap to map from surf2
128 // region numbers into combined patchList
129 Foam::geometricSurfacePatchList Foam::booleanSurface::mergePatches
130 (
131  const triSurface& surf1,
132  const triSurface& surf2,
133  labelList& patchMap2
134 )
135 {
136  // Size too big.
137  geometricSurfacePatchList combinedPatches
138  (
139  surf1.patches().size()
140  + surf2.patches().size()
141  );
142 
143  // Copy all patches of surf1
144  label combinedPatchI = 0;
145  forAll(surf1.patches(), patchI)
146  {
147  combinedPatches[combinedPatchI++] = surf1.patches()[patchI];
148  }
149 
150  // (inefficiently) add unique patches from surf2
151  patchMap2.setSize(surf2.patches().size());
152 
153  forAll(surf2.patches(), patch2I)
154  {
155  label index = -1;
156 
157  forAll(surf1.patches(), patch1I)
158  {
159  if (surf1.patches()[patch1I] == surf2.patches()[patch2I])
160  {
161  index = patch1I;
162 
163  break;
164  }
165  }
166 
167  if (index == -1)
168  {
169  combinedPatches[combinedPatchI] = surf2.patches()[patch2I];
170  patchMap2[patch2I] = combinedPatchI;
171  combinedPatchI++;
172  }
173  else
174  {
175  patchMap2[patch2I] = index;
176  }
177  }
178 
179  combinedPatches.setSize(combinedPatchI);
180 
181  return combinedPatches;
182 }
183 
184 
185 void Foam::booleanSurface::propagateEdgeSide
186 (
187  const triSurface& surf,
188  const label prevVert0,
189  const label prevFaceI,
190  const label prevState,
191  const label edgeI,
192  labelList& side
193 )
194 {
195  const labelList& eFaces = surf.sortedEdgeFaces()[edgeI];
196 
197  // Simple case. Propagate side.
198  if (eFaces.size() == 2)
199  {
200  forAll(eFaces, eFaceI)
201  {
202  propagateSide
203  (
204  surf,
205  prevState,
206  eFaces[eFaceI],
207  side
208  );
209  }
210  }
211 
212 
213  if (((eFaces.size() % 2) == 1) && (eFaces.size() != 1))
214  {
216  (
217  "booleanSurface::propagateEdgeSide(const triSurface&,"
218  "const label, const label, const label, const label,"
219  " labelList&)"
220  ) << "Don't know how to handle edges with odd number of faces"
221  << endl
222  << "edge:" << edgeI << " vertices:" << surf.edges()[edgeI]
223  << " coming from face:" << prevFaceI
224  << " edgeFaces:" << eFaces << abort(FatalError);
225  }
226 
227 
228  // Get position of face in edgeFaces
229  label ind = index(eFaces, prevFaceI);
230 
231  // Determine orientation of faces around edge prevVert0
232  // (might be opposite of edge)
233  const edge& e = surf.edges()[edgeI];
234 
235  // Get next face to include
236  label nextInd;
237  label prevInd;
238 
239  if (e.start() == prevVert0)
240  {
241  // Edge (and hence eFaces) in same order as prevVert0.
242  // Take next face from sorted list
243  nextInd = eFaces.fcIndex(ind);
244  prevInd = eFaces.rcIndex(ind);
245  }
246  else
247  {
248  // Take previous face from sorted neighbours
249  nextInd = eFaces.rcIndex(ind);
250  prevInd = eFaces.fcIndex(ind);
251  }
252 
253 
254  if (prevState == OUTSIDE)
255  {
256  // Coming from outside. nextInd is outside, rest is inside.
257 
258  forAll(eFaces, eFaceI)
259  {
260  if (eFaceI != ind)
261  {
262  label nextState;
263 
264  if (eFaceI == nextInd)
265  {
266  nextState = OUTSIDE;
267  }
268  else
269  {
270  nextState = INSIDE;
271  }
272 
273  propagateSide
274  (
275  surf,
276  nextState,
277  eFaces[eFaceI],
278  side
279  );
280  }
281  }
282  }
283  else
284  {
285  // Coming from inside. prevInd is inside as well, rest is outside.
286 
287  forAll(eFaces, eFaceI)
288  {
289  if (eFaceI != ind)
290  {
291  label nextState;
292 
293  if (eFaceI == prevInd)
294  {
295  nextState = INSIDE;
296  }
297  else
298  {
299  nextState = OUTSIDE;
300  }
301 
302  propagateSide
303  (
304  surf,
305  nextState,
306  eFaces[eFaceI],
307  side
308  );
309  }
310  }
311  }
312 }
313 
314 
315 // Face-edge walk. Determines inside/outside for all faces connected to an edge.
316 void Foam::booleanSurface::propagateSide
317 (
318  const triSurface& surf,
319  const label prevState,
320  const label faceI,
321  labelList& side
322 )
323 {
324  if (side[faceI] == UNVISITED)
325  {
326  side[faceI] = prevState;
327 
328  const labelledTri& tri = surf[faceI];
329 
330  // Get copy of face labels
331  label a = tri[0];
332  label b = tri[1];
333  label c = tri[2];
334 
335  // Go and visit my edges' face-neighbours.
336 
337  const labelList& myEdges = surf.faceEdges()[faceI];
338 
339  label edgeAB = findEdge(surf.edges(), myEdges, edge(a, b));
340 
341  propagateEdgeSide
342  (
343  surf,
344  a,
345  faceI,
346  prevState,
347  edgeAB,
348  side
349  );
350 
351  label edgeBC = findEdge(surf.edges(), myEdges, edge(b, c));
352 
353  propagateEdgeSide
354  (
355  surf,
356  b,
357  faceI,
358  prevState,
359  edgeBC,
360  side
361  );
362 
363  label edgeCA = findEdge(surf.edges(), myEdges, edge(c, a));
364 
365  propagateEdgeSide
366  (
367  surf,
368  c,
369  faceI,
370  prevState,
371  edgeCA,
372  side
373  );
374  }
375 }
376 
377 
378 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
379 
380 // Null constructor
382 :
383  triSurface()
384 {}
385 
386 
387 // Construct from surfaces and face to include for every surface
389 (
390  const triSurface& surf1,
391  const triSurface& surf2,
392  const surfaceIntersection& inter,
393  const label includeFace1,
394  const label includeFace2
395 )
396 :
397  triSurface(),
398  faceMap_()
399 {
400  if (debug)
401  {
402  Pout<< "booleanSurface : Generating intersected surface for surf1"
403  << endl;
404  }
405 
406  // Add intersection to surface1 (retriangulates cut faces)
407  intersectedSurface cutSurf1(surf1, true, inter);
408 
409 
410  if (debug)
411  {
412  Pout<< "booleanSurface : Generated cutSurf1: " << endl;
413  cutSurf1.writeStats(Pout);
414 
415  Pout<< "Writing to file cutSurf1.ftr" << endl;
416  OFstream cutSurf1Stream("cutSurf1.ftr");
417  cutSurf1.write(cutSurf1Stream);
418  }
419 
420  if (debug)
421  {
422  Pout<< "booleanSurface : Generating intersected surface for surf2"
423  << endl;
424  }
425 
426  // Add intersection to surface2
427  intersectedSurface cutSurf2(surf2, false, inter);
428 
429  if (debug)
430  {
431  Pout<< "booleanSurface : Generated cutSurf2: " << endl;
432  cutSurf2.writeStats(Pout);
433 
434  Pout<< "Writing to file cutSurf2.ftr" << endl;
435  OFstream cutSurf2Stream("cutSurf2.ftr");
436  cutSurf2.write(cutSurf2Stream);
437  }
438 
439 
440  // Find (first) face of cutSurf1 that originates from includeFace1
441  label cutSurf1FaceI = index(cutSurf1.faceMap(), includeFace1);
442 
443  if (debug)
444  {
445  Pout<< "cutSurf1 : starting to fill from face:" << cutSurf1FaceI
446  << endl;
447  }
448 
449  if (cutSurf1FaceI == -1)
450  {
452  (
453  "booleanSurface(const triSurfaceSearch&"
454  ", const label, const triSurfaceSearch&, const label)"
455  ) << "Did not find face with label " << includeFace1
456  << " in intersectedSurface."
457  << exit(FatalError);
458  }
459 
460  // Find (first) face of cutSurf2 that originates from includeFace1
461  label cutSurf2FaceI = index(cutSurf2.faceMap(), includeFace2);
462 
463  if (debug)
464  {
465  Pout<< "cutSurf2 : starting to fill from face:" << cutSurf2FaceI
466  << endl;
467  }
468  if (cutSurf2FaceI == -1)
469  {
471  (
472  "booleanSurface(const triSurfaceSearch&"
473  ", const label, const triSurfaceSearch&, const label)"
474  ) << "Did not find face with label " << includeFace2
475  << " in intersectedSurface."
476  << exit(FatalError);
477  }
478 
479 
480  //
481  // Mark faces of cutSurf1 that need to be kept by walking from includeFace1
482  // without crossing any edges of the intersection.
483  //
484 
485  // Mark edges on intersection
486  const labelList& int1Edges = cutSurf1.intersectionEdges();
487 
488  boolList isIntersectionEdge1(cutSurf1.nEdges(), false);
489  forAll(int1Edges, intEdgeI)
490  {
491  label edgeI = int1Edges[intEdgeI];
492  isIntersectionEdge1[edgeI] = true;
493  }
494 
495  labelList faceZone1;
496  (void)cutSurf1.markZones(isIntersectionEdge1, faceZone1);
497 
498 
499  // Check whether at least one of sides of intersection has been marked.
500  checkIncluded(cutSurf1, faceZone1, cutSurf1FaceI);
501 
502  // Subset zone which includes cutSurf2FaceI
503  boolList includedFaces1(cutSurf1.size(), false);
504 
505  forAll(faceZone1, faceI)
506  {
507  if (faceZone1[faceI] == faceZone1[cutSurf1FaceI])
508  {
509  includedFaces1[faceI] = true;
510  }
511  }
512 
513  // Subset to include only interesting part
514  labelList pointMap1;
515  labelList faceMap1;
516 
517  triSurface subSurf1
518  (
519  cutSurf1.subsetMesh
520  (
521  includedFaces1,
522  pointMap1,
523  faceMap1
524  )
525  );
526 
527 
528  //
529  // Mark faces of cutSurf2 that need to be kept by walking from includeFace2
530  // without crossing any edges of the intersection.
531  //
532 
533  // Mark edges and points on intersection
534  const labelList& int2Edges = cutSurf2.intersectionEdges();
535 
536  boolList isIntersectionEdge2(cutSurf2.nEdges(), false);
537  forAll(int2Edges, intEdgeI)
538  {
539  label edgeI = int2Edges[intEdgeI];
540  isIntersectionEdge2[edgeI] = true;
541  }
542 
543  labelList faceZone2;
544  (void)cutSurf2.markZones(isIntersectionEdge2, faceZone2);
545 
546 
547  // Check whether at least one of sides of intersection has been marked.
548  checkIncluded(cutSurf2, faceZone2, cutSurf2FaceI);
549 
550  // Subset zone which includes cutSurf2FaceI
551  boolList includedFaces2(cutSurf2.size(), false);
552 
553  forAll(faceZone2, faceI)
554  {
555  if (faceZone2[faceI] == faceZone2[cutSurf2FaceI])
556  {
557  includedFaces2[faceI] = true;
558  }
559  }
560 
561  labelList pointMap2;
562  labelList faceMap2;
563 
564  triSurface subSurf2
565  (
566  cutSurf2.subsetMesh
567  (
568  includedFaces2,
569  pointMap2,
570  faceMap2
571  )
572  );
573 
574 
575  //
576  // Now match up the corresponding points on the intersection. The
577  // intersectedSurfaces will have the points resulting from the
578  // intersection last in their points and in the same
579  // order so we can use the pointMaps from the subsets to find them.
580  //
581  // We keep the vertices on the first surface and renumber those on the
582  // second one.
583 
584 
585  //
586  // points
587  //
588  pointField combinedPoints
589  (
590  subSurf1.nPoints()
591  + subSurf2.nPoints()
592  - (cutSurf2.nPoints() - cutSurf2.nSurfacePoints())
593  );
594 
595  // Copy points from subSurf1 and remember the labels of the ones in
596  // the intersection
597  labelList intersectionLabels(cutSurf1.nPoints() - cutSurf1.nSurfacePoints());
598 
599  label combinedPointI = 0;
600 
601  forAll(subSurf1.points(), pointI)
602  {
603  // Label in cutSurf
604  label cutSurfPointI = pointMap1[pointI];
605 
606  if (!cutSurf1.isSurfacePoint(cutSurfPointI))
607  {
608  // Label in original intersection is equal to the cutSurfPointI
609 
610  // Remember label in combinedPoints for intersection point.
611  intersectionLabels[cutSurfPointI] = combinedPointI;
612  }
613 
614  // Copy point
615  combinedPoints[combinedPointI++] = subSurf1.points()[pointI];
616  }
617 
618  // Append points from subSurf2 (if they are not intersection points)
619  // and construct mapping
620  labelList pointMap(subSurf2.nPoints());
621 
622  forAll(subSurf2.points(), pointI)
623  {
624  // Label in cutSurf
625  label cutSurfPointI = pointMap2[pointI];
626 
627  if (!cutSurf2.isSurfacePoint(cutSurfPointI))
628  {
629  // Lookup its label in combined point list.
630  pointMap[pointI] = intersectionLabels[cutSurfPointI];
631  }
632  else
633  {
634  pointMap[pointI] = combinedPointI;
635 
636  combinedPoints[combinedPointI++] = subSurf2.points()[pointI];
637  }
638  }
639 
640 
641  //
642  // patches
643  //
644 
645  labelList patchMap2;
646 
647  geometricSurfacePatchList combinedPatches
648  (
649  mergePatches
650  (
651  surf1,
652  surf2,
653  patchMap2
654  )
655  );
656 
657 
658  //
659  // faces
660  //
661 
662  List<labelledTri> combinedFaces(subSurf1.size() + subSurf2.size());
663 
664  faceMap_.setSize(combinedFaces.size());
665 
666  // Copy faces from subSurf1. No need for renumbering.
667  label combinedFaceI = 0;
668  forAll(subSurf1, faceI)
669  {
670  faceMap_[combinedFaceI] = faceMap1[faceI];
671  combinedFaces[combinedFaceI++] = subSurf1[faceI];
672  }
673 
674  // Copy and renumber faces from subSurf2.
675  forAll(subSurf2, faceI)
676  {
677  const labelledTri& f = subSurf2[faceI];
678 
679  faceMap_[combinedFaceI] = -faceMap2[faceI]-1;
680 
681  combinedFaces[combinedFaceI++] =
683  (
684  pointMap[f[0]],
685  pointMap[f[1]],
686  pointMap[f[2]],
687  patchMap2[f.region()]
688  );
689  }
690 
692  (
693  triSurface
694  (
695  combinedFaces,
696  combinedPatches,
697  combinedPoints
698  )
699  );
700 }
701 
702 
703 // Construct from surfaces and boolean operation
705 (
706  const triSurface& surf1,
707  const triSurface& surf2,
708  const surfaceIntersection& inter,
709  const label booleanOp
710 )
711 :
712  triSurface(),
713  faceMap_()
714 {
715  if (debug)
716  {
717  Pout<< "booleanSurface : Testing surf1 and surf2" << endl;
718 
719  {
720  const labelListList& edgeFaces = surf1.edgeFaces();
721 
722  forAll(edgeFaces, edgeI)
723  {
724  const labelList& eFaces = edgeFaces[edgeI];
725 
726  if (eFaces.size() == 1)
727  {
728  WarningIn("booleanSurface::booleanSurface")
729  << "surf1 is open surface at edge " << edgeI
730  << " verts:" << surf1.edges()[edgeI]
731  << " connected to faces " << eFaces << endl;
732  }
733  }
734  }
735  {
736  const labelListList& edgeFaces = surf2.edgeFaces();
737 
738  forAll(edgeFaces, edgeI)
739  {
740  const labelList& eFaces = edgeFaces[edgeI];
741 
742  if (eFaces.size() == 1)
743  {
744  WarningIn("booleanSurface::booleanSurface")
745  << "surf2 is open surface at edge " << edgeI
746  << " verts:" << surf2.edges()[edgeI]
747  << " connected to faces " << eFaces << endl;
748  }
749  }
750  }
751  }
752 
753 
754  //
755  // Surface 1
756  //
757 
758  if (debug)
759  {
760  Pout<< "booleanSurface : Generating intersected surface for surf1"
761  << endl;
762  }
763 
764  // Add intersection to surface1 (retriangulates cut faces)
765  intersectedSurface cutSurf1(surf1, true, inter);
766 
767  if (debug)
768  {
769  Pout<< "booleanSurface : Generated cutSurf1: " << endl;
770  cutSurf1.writeStats(Pout);
771 
772  Pout<< "Writing to file cutSurf1.ftr" << endl;
773  OFstream cutSurf1Stream("cutSurf1.ftr");
774  cutSurf1.write(cutSurf1Stream);
775  }
776 
777 
778  //
779  // Surface 2
780  //
781 
782  if (debug)
783  {
784  Pout<< "booleanSurface : Generating intersected surface for surf2"
785  << endl;
786  }
787 
788  // Add intersection to surface2
789  intersectedSurface cutSurf2(surf2, false, inter);
790 
791  if (debug)
792  {
793  Pout<< "booleanSurface : Generated cutSurf2: " << endl;
794  cutSurf2.writeStats(Pout);
795 
796  Pout<< "Writing to file cutSurf2.ftr" << endl;
797  OFstream cutSurf2Stream("cutSurf2.ftr");
798  cutSurf2.write(cutSurf2Stream);
799  }
800 
801 
802  //
803  // patches
804  //
805 
806  labelList patchMap2;
807 
808  geometricSurfacePatchList combinedPatches
809  (
810  mergePatches
811  (
812  surf1,
813  surf2,
814  patchMap2
815  )
816  );
817 
818 
819  //
820  // Now match up the corresponding points on the intersection. The
821  // intersectedSurfaces will have the points resulting from the
822  // intersection first in their points and in the same
823  // order
824  //
825  // We keep the vertices on the first surface and renumber those on the
826  // second one.
827 
828  pointField combinedPoints(cutSurf1.nPoints() + cutSurf2.nSurfacePoints());
829 
830  // Copy all points from 1 and non-intersection ones from 2.
831 
832  label combinedPointI = 0;
833 
834  forAll(cutSurf1.points(), pointI)
835  {
836  combinedPoints[combinedPointI++] = cutSurf1.points()[pointI];
837  }
838 
839  for
840  (
841  label pointI = 0;
842  pointI < cutSurf2.nSurfacePoints();
843  pointI++
844  )
845  {
846  combinedPoints[combinedPointI++] = cutSurf2.points()[pointI];
847  }
848 
849  // Point order is now
850  // - 0.. cutSurf1.nSurfacePoints : original surf1 points
851  // - .. cutSurf1.nPoints : intersection points
852  // - .. cutSurf2.nSurfacePoints : original surf2 points
853 
854  if (debug)
855  {
856  Pout<< "booleanSurface : generated points:" << nl
857  << " 0 .. " << cutSurf1.nSurfacePoints()-1
858  << " : original surface1"
859  << nl
860  << " " << cutSurf1.nSurfacePoints()
861  << " .. " << cutSurf1.nPoints()-1
862  << " : intersection points"
863  << nl
864  << " " << cutSurf1.nPoints() << " .. "
865  << cutSurf2.nSurfacePoints()-1
866  << " : surface2 points"
867  << nl
868  << endl;
869  }
870 
871  // Copy faces. Faces from surface 1 keep vertex numbering and region info.
872  // Faces from 2 get vertices and region renumbered.
873  List<labelledTri> combinedFaces(cutSurf1.size() + cutSurf2.size());
874 
875  label combinedFaceI = 0;
876 
877  forAll(cutSurf1, faceI)
878  {
879  combinedFaces[combinedFaceI++] = cutSurf1[faceI];
880  }
881 
882  forAll(cutSurf2, faceI)
883  {
884  labelledTri& combinedTri = combinedFaces[combinedFaceI++];
885 
886  const labelledTri& tri = cutSurf2[faceI];
887 
888  forAll(tri, fp)
889  {
890  if (cutSurf2.isSurfacePoint(tri[fp]))
891  {
892  // Renumber. Surface2 points are after ones from surf 1.
893  combinedTri[fp] = tri[fp] + cutSurf1.nPoints();
894  }
895  else
896  {
897  // Is intersection.
898  combinedTri[fp] =
899  tri[fp]
900  - cutSurf2.nSurfacePoints()
901  + cutSurf1.nSurfacePoints();
902  }
903  }
904  combinedTri.region() = patchMap2[tri.region()];
905  }
906 
907 
908  // Now we have surface in combinedFaces and combinedPoints. Use
909  // booleanOp to determine which part of what to keep.
910 
911  // Construct addressing for whole part.
912  triSurface combinedSurf
913  (
914  combinedFaces,
915  combinedPatches,
916  combinedPoints
917  );
918 
919  if (debug)
920  {
921  Pout<< "booleanSurface : Generated combinedSurf: " << endl;
922  combinedSurf.writeStats(Pout);
923 
924  Pout<< "Writing to file combinedSurf.ftr" << endl;
925  OFstream combinedSurfStream("combinedSurf.ftr");
926  combinedSurf.write(combinedSurfStream);
927  }
928 
929 
930  if (booleanOp == booleanSurface::ALL)
931  {
932  // Special case: leave surface multiply connected
933 
934  faceMap_.setSize(combinedSurf.size());
935 
936  label combinedFaceI = 0;
937 
938  forAll(cutSurf1, faceI)
939  {
940  faceMap_[combinedFaceI++] = cutSurf1.faceMap()[faceI];
941  }
942  forAll(cutSurf2, faceI)
943  {
944  faceMap_[combinedFaceI++] = -cutSurf2.faceMap()[faceI] - 1;
945  }
946 
947  triSurface::operator=(combinedSurf);
948 
949  return;
950  }
951 
952 
953  // Get outside point.
954  point outsidePoint = 2 * treeBoundBox(combinedSurf.localPoints()).span();
955 
956  //
957  // Linear search for nearest point on surface.
958  //
959 
960  const pointField& pts = combinedSurf.points();
961 
962  label minFaceI = -1;
963  pointHit minHit(false, vector::zero, GREAT, true);
964 
965  forAll(combinedSurf, faceI)
966  {
967  const labelledTri& f = combinedSurf[faceI];
968 
969  pointHit curHit =
971  (
972  pts[f[0]],
973  pts[f[1]],
974  pts[f[2]]
975  ).nearestPoint(outsidePoint);
976 
977  if (curHit.distance() < minHit.distance())
978  {
979  minHit = curHit;
980 
981  minFaceI = faceI;
982  }
983  }
984 
985  if (debug)
986  {
987  Pout<< "booleanSurface : found for point:" << outsidePoint
988  << " nearest face:" << minFaceI
989  << " nearest point:" << minHit.rawPoint()
990  << endl;
991  }
992 
993  // Visibility/side of face:
994  // UNVISITED: unvisited
995  // OUTSIDE: visible from outside
996  // INSIDE: invisible from outside
997  labelList side(combinedSurf.size(), UNVISITED);
998 
999  // Walk face-edge-face and propagate inside/outside status.
1000  propagateSide(combinedSurf, OUTSIDE, minFaceI, side);
1001 
1002 
1003  // Depending on operation include certain faces.
1004  // AND: faces on inside of 1 and of 2
1005  // OR: faces on outside of 1 and of 2
1006  // XOR: faces on outside of 1 and inside of 2
1007 
1008  boolList include(combinedSurf.size(), false);
1009 
1010  forAll(side, faceI)
1011  {
1012  if (side[faceI] == UNVISITED)
1013  {
1014  FatalErrorIn
1015  (
1016  "booleanSurface::booleanSurface"
1017  "(const triSurfaceSearch&, const triSurfaceSearch&"
1018  ", const label booleanOp)"
1019  ) << "Face " << faceI << " has not been reached by walking from"
1020  << " nearest point " << minHit.rawPoint()
1021  << " nearest face " << minFaceI << exit(FatalError);
1022  }
1023  else if (side[faceI] == OUTSIDE)
1024  {
1025  if (booleanOp == booleanSurface::OR)
1026  {
1027  include[faceI] = true;
1028  }
1029  else if (booleanOp == booleanSurface::AND)
1030  {
1031  include[faceI] = false;
1032  }
1033  else // xor
1034  {
1035  include[faceI] = (faceI < cutSurf1.size()); // face from surf1
1036  }
1037  }
1038  else // inside
1039  {
1040  if (booleanOp == booleanSurface::OR)
1041  {
1042  include[faceI] = false;
1043  }
1044  else if (booleanOp == booleanSurface::AND)
1045  {
1046  include[faceI] = true;
1047  }
1048  else // xor
1049  {
1050  include[faceI] = (faceI >= cutSurf1.size()); // face from surf2
1051  }
1052  }
1053  }
1054 
1055  // Create subsetted surface
1056  labelList subToCombinedPoint;
1057  labelList subToCombinedFace;
1058  triSurface subSurf
1059  (
1060  combinedSurf.subsetMesh
1061  (
1062  include,
1063  subToCombinedPoint,
1064  subToCombinedFace
1065  )
1066  );
1067 
1068  // Create face map
1069  faceMap_.setSize(subSurf.size());
1070 
1071  forAll(subToCombinedFace, faceI)
1072  {
1073  // Get label in combinedSurf
1074  label combinedFaceI = subToCombinedFace[faceI];
1075 
1076  // First faces in combinedSurf come from cutSurf1.
1077 
1078  if (combinedFaceI < cutSurf1.size())
1079  {
1080  label cutSurf1Face = combinedFaceI;
1081 
1082  faceMap_[faceI] = cutSurf1.faceMap()[cutSurf1Face];
1083  }
1084  else
1085  {
1086  label cutSurf2Face = combinedFaceI - cutSurf1.size();
1087 
1088  faceMap_[faceI] = - cutSurf2.faceMap()[cutSurf2Face] - 1;
1089  }
1090  }
1091 
1092  // Orient outwards
1093  orientedSurface outSurf(subSurf);
1094 
1095  // Assign.
1096  triSurface::operator=(outSurf);
1097 }
1098 
1099 
1100 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
1101 
1102 
1103 // ************************ vim: set sw=4 sts=4 et: ************************ //