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