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