surfaceBooleanFeatures.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) 2016-2022 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
27Application
28 surfaceBooleanFeatures
29
30Group
31 grpSurfaceUtilities
32
33Description
34 Generates the extendedFeatureEdgeMesh for the interface between a boolean
35 operation on two surfaces.
36
37 Assumes that the orientation of the surfaces is correct:
38 - if the operation is union or intersection, that both surface's normals
39 (n) have the same orientation with respect to a point, i.e. surfaces
40 A and B are orientated the same with respect to point x:
41
42 \verbatim
43 _______
44 | |--> n
45 | ___|___ x
46 |A | | |--> n
47 |___|___| B|
48 | |
49 |_______|
50
51 \endverbatim
52
53 - if the operation is a subtraction, the surfaces should be oppositely
54 oriented with respect to a point, i.e. for (A - B), then B's orientation
55 should be such that x is "inside", and A's orientation such that x is
56 "outside"
57
58 \verbatim
59 _______
60 | |--> n
61 | ___|___ x
62 |A | | |
63 |___|___| B|
64 | n <--|
65 |_______|
66
67 \endverbatim
68
69 When the operation is performed - for union, all of the edges generated
70 where one surfaces cuts another are all "internal" for union,
71 and "external" for intersection, (B - A) and (A - B).
72 This has been assumed, formal (dis)proof is invited.
73
74\*---------------------------------------------------------------------------*/
75
76#include "triSurface.H"
77#include "argList.H"
78#include "Time.H"
79#include "featureEdgeMesh.H"
81#include "triSurfaceSearch.H"
82#include "triSurfaceMesh.H"
83#include "OFstream.H"
84#include "OBJstream.H"
85#include "booleanSurface.H"
86#include "edgeIntersections.H"
87#include "meshTools.H"
88#include "DynamicField.H"
89#include "Enum.H"
90
91#ifndef NO_CGAL
92
93// Silence boost bind deprecation warnings (before CGAL-5.2.1)
94#include "CGAL/version.h"
95#if defined(CGAL_VERSION_NR) && (CGAL_VERSION_NR < 1050211000)
96#define BOOST_BIND_GLOBAL_PLACEHOLDERS
97#endif
98
99#include <CGAL/AABB_tree.h>
100#include <CGAL/AABB_traits.h>
101#include <CGAL/AABB_face_graph_triangle_primitive.h>
103#include "PolyhedronReader.H"
104typedef CGAL::AABB_face_graph_triangle_primitive
105<
106 Polyhedron, CGAL::Default, CGAL::Tag_false
107> Primitive;
108typedef CGAL::AABB_traits<K, Primitive> Traits;
109typedef CGAL::AABB_tree<Traits> Tree;
110
111typedef boost::optional
112<
113 Tree::Intersection_and_primitive_id<Segment>::Type
114> Segment_intersection;
115
116#endif // NO_CGAL
117
118
119using namespace Foam;
120
121// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
122
123bool intersectSurfaces
124(
125 triSurface& surf,
126 edgeIntersections& edgeCuts
127)
128{
129 bool hasMoved = false;
130
131 for (label iter = 0; iter < 10; iter++)
132 {
133 Info<< "Determining intersections of surface edges with itself" << endl;
134
135 // Determine surface edge intersections. Allow surface to be moved.
136
137 // Number of iterations needed to resolve degenerates
138 label nIters = 0;
139 {
140 triSurfaceSearch querySurf(surf);
141
142 scalarField surfPointTol
143 (
144 max(1e-3*edgeIntersections::minEdgeLength(surf), SMALL)
145 );
146
147 // Determine raw intersections
148 edgeCuts = edgeIntersections
149 (
150 surf,
151 querySurf,
152 surfPointTol
153 );
154
155 // Shuffle a bit to resolve degenerate edge-face hits
156 {
157 pointField points(surf.points());
158
159 nIters =
160 edgeCuts.removeDegenerates
161 (
162 5, // max iterations
163 surf,
164 querySurf,
165 surfPointTol,
166 points // work array
167 );
168
169 if (nIters != 0)
170 {
171 // Update geometric quantities
172 surf.movePoints(points);
173 hasMoved = true;
174 }
175 }
176 }
177 }
178
179 if (hasMoved)
180 {
181 fileName newFile("surf.obj");
182 Info<< "Surface has been moved. Writing to " << newFile << endl;
183 surf.write(newFile);
184 }
185
186 return hasMoved;
187}
188
189
190// Keep on shuffling surface points until no more degenerate intersections.
191// Moves both surfaces and updates set of edge cuts.
192bool intersectSurfaces
193(
194 triSurface& surf1,
195 edgeIntersections& edgeCuts1,
196 triSurface& surf2,
197 edgeIntersections& edgeCuts2
198)
199{
200 bool hasMoved1 = false;
201 bool hasMoved2 = false;
202
203 for (label iter = 0; iter < 10; iter++)
204 {
205 Info<< "Determining intersections of surf1 edges with surf2"
206 << " faces" << endl;
207
208 // Determine surface1 edge intersections. Allow surface to be moved.
209
210 // Number of iterations needed to resolve degenerates
211 label nIters1 = 0;
212 {
213 triSurfaceSearch querySurf2(surf2);
214
215 scalarField surf1PointTol
216 (
217 max(1e-3*edgeIntersections::minEdgeLength(surf1), SMALL)
218 );
219
220 // Determine raw intersections
221 edgeCuts1 = edgeIntersections
222 (
223 surf1,
224 querySurf2,
225 surf1PointTol
226 );
227
228 // Shuffle a bit to resolve degenerate edge-face hits
229 {
230 pointField points1(surf1.points());
231
232 nIters1 =
233 edgeCuts1.removeDegenerates
234 (
235 5, // max iterations
236 surf1,
237 querySurf2,
238 surf1PointTol,
239 points1 // work array
240 );
241
242 if (nIters1 != 0)
243 {
244 // Update geometric quantities
245 surf1.movePoints(points1);
246 hasMoved1 = true;
247 }
248 }
249 }
250
251 Info<< "Determining intersections of surf2 edges with surf1"
252 << " faces" << endl;
253
254 label nIters2 = 0;
255 {
256 triSurfaceSearch querySurf1(surf1);
257
258 scalarField surf2PointTol
259 (
260 max(1e-3*edgeIntersections::minEdgeLength(surf2), SMALL)
261 );
262
263 // Determine raw intersections
264 edgeCuts2 = edgeIntersections
265 (
266 surf2,
267 querySurf1,
268 surf2PointTol
269 );
270
271 // Shuffle a bit to resolve degenerate edge-face hits
272 {
273 pointField points2(surf2.points());
274
275 nIters2 =
276 edgeCuts2.removeDegenerates
277 (
278 5, // max iterations
279 surf2,
280 querySurf1,
281 surf2PointTol,
282 points2 // work array
283 );
284
285 if (nIters2 != 0)
286 {
287 // Update geometric quantities
288 surf2.movePoints(points2);
289 hasMoved2 = true;
290 }
291 }
292 }
293
294 if (nIters1 == 0 && nIters2 == 0)
295 {
296 Info<< "** Resolved all intersections to be proper edge-face pierce"
297 << endl;
298 break;
299 }
300 }
301
302 if (hasMoved1)
303 {
304 fileName newFile("surf1.obj");
305 Info<< "Surface 1 has been moved. Writing to " << newFile
306 << endl;
307 surf1.write(newFile);
308 }
309
310 if (hasMoved2)
311 {
312 fileName newFile("surf2.obj");
313 Info<< "Surface 2 has been moved. Writing to " << newFile
314 << endl;
315 surf2.write(newFile);
316 }
317
318 return hasMoved1 || hasMoved2;
319}
320
321
322label calcNormalDirection
323(
324 const vector& normal,
325 const vector& otherNormal,
326 const vector& edgeDir,
327 const vector& faceCentre,
328 const vector& pointOnEdge
329)
330{
331 const vector cross = normalised(normal ^ edgeDir);
332
333 const vector fC0tofE0 = normalised(faceCentre - pointOnEdge);
334
335 label nDir = ((cross & fC0tofE0) > 0.0 ? 1 : -1);
336
337 nDir *= ((otherNormal & fC0tofE0) > 0.0 ? -1 : 1);
338
339 return nDir;
340}
341
342
343void calcEdgeCuts
344(
345 triSurface& surf1,
346 triSurface& surf2,
347 const bool perturb,
348 edgeIntersections& edgeCuts1,
349 edgeIntersections& edgeCuts2
350)
351{
352 if (perturb)
353 {
354 intersectSurfaces
355 (
356 surf1,
357 edgeCuts1,
358 surf2,
359 edgeCuts2
360 );
361 }
362 else
363 {
364 triSurfaceSearch querySurf2(surf2);
365
366 Info<< "Determining intersections of surf1 edges with surf2 faces"
367 << endl;
368
369 edgeCuts1 = edgeIntersections
370 (
371 surf1,
372 querySurf2,
373 max(1e-3*edgeIntersections::minEdgeLength(surf1), SMALL)
374 );
375
376 triSurfaceSearch querySurf1(surf1);
377
378 Info<< "Determining intersections of surf2 edges with surf1 faces"
379 << endl;
380
381 edgeCuts2 = edgeIntersections
382 (
383 surf2,
384 querySurf1,
385 max(1e-3*edgeIntersections::minEdgeLength(surf2), SMALL)
386 );
387 }
388}
389
390
391// CGAL variants
392
393#ifndef NO_CGAL
394
395void visitPointRegion
396(
397 const triSurface& s,
398 const label zoneI,
399 const label pointI,
400 const label startEdgeI,
401 const label startFaceI,
402 labelList& pFacesZone
403)
404{
405 const labelList& eFaces = s.edgeFaces()[startEdgeI];
406
407 if (eFaces.size() == 2)
408 {
409 label nextFaceI;
410 if (eFaces[0] == startFaceI)
411 {
412 nextFaceI = eFaces[1];
413 }
414 else if (eFaces[1] == startFaceI)
415 {
416 nextFaceI = eFaces[0];
417 }
418 else
419 {
421 << "problem" << exit(FatalError);
422
423 nextFaceI = -1;
424 }
425
426
427
428 label index = s.pointFaces()[pointI].find(nextFaceI);
429
430 if (pFacesZone[index] == -1)
431 {
432 // Mark face as been visited.
433 pFacesZone[index] = zoneI;
434
435 // Step to next edge on face which is still using pointI
436 const labelList& fEdges = s.faceEdges()[nextFaceI];
437
438 label nextEdgeI = -1;
439
440 forAll(fEdges, i)
441 {
442 label edgeI = fEdges[i];
443 const edge& e = s.edges()[edgeI];
444
445 if (edgeI != startEdgeI && (e[0] == pointI || e[1] == pointI))
446 {
447 nextEdgeI = edgeI;
448
449 break;
450 }
451 }
452
453 if (nextEdgeI == -1)
454 {
456 << "Problem: cannot find edge out of " << fEdges
457 << "on face " << nextFaceI << " that uses point " << pointI
458 << " and is not edge " << startEdgeI << abort(FatalError);
459 }
460
461
462 visitPointRegion
463 (
464 s,
465 zoneI,
466 pointI,
467 nextEdgeI,
468 nextFaceI,
469 pFacesZone
470 );
471 }
472 }
473}
474
475
476label dupNonManifoldPoints(triSurface& s, labelList& pointMap)
477{
478 const labelListList& pf = s.pointFaces();
479 const labelListList& fe = s.faceEdges();
480 const edgeList& edges = s.edges();
481
482
483 DynamicField<point> newPoints(s.points());
484 // From dupSurf back to s.pointa
485 DynamicList<label> newPointMap(identity(newPoints.size()));
486 List<labelledTri> newFaces(s);
487 label nNonManifold = 0;
488
489 forAll(pf, pointI)
490 {
491 const labelList& pFaces = pf[pointI];
492
493 // Visited faces (as indices into pFaces)
494 labelList pFacesZone(pFaces.size(), -1);
495
496 label nZones = 0;
497 label index = 0;
498 for (; index < pFacesZone.size(); index++)
499 {
500 if (pFacesZone[index] == -1)
501 {
502 label zoneI = nZones++;
503 pFacesZone[index] = zoneI;
504
505 label faceI = pFaces[index];
506 const labelList& fEdges = fe[faceI];
507
508 // Find starting edge
509 forAll(fEdges, fEdgeI)
510 {
511 label edgeI = fEdges[fEdgeI];
512 const edge& e = edges[edgeI];
513
514 if (e[0] == pointI || e[1] == pointI)
515 {
516 visitPointRegion
517 (
518 s,
519 zoneI,
520 pointI,
521 edgeI,
522 faceI,
523 pFacesZone
524 );
525 }
526 }
527 }
528 }
529
530
531 // Subset
532 if (nZones > 1)
533 {
534 for (label zoneI = 1; zoneI < nZones; zoneI++)
535 {
536 label newPointI = newPoints.size();
537 newPointMap.append(s.meshPoints()[pointI]);
538 newPoints.append(s.points()[s.meshPoints()[pointI]]);
539
540 forAll(pFacesZone, index)
541 {
542 if (pFacesZone[index] == zoneI)
543 {
544 label faceI = pFaces[index];
545 const labelledTri& localF = s.localFaces()[faceI];
546 forAll(localF, fp)
547 {
548 if (localF[fp] == pointI)
549 {
550 newFaces[faceI][fp] = newPointI;
551 }
552 }
553 }
554 }
555 }
556 nNonManifold++;
557 }
558 }
559
560
561 Info<< "Duplicating " << nNonManifold << " points out of " << s.nPoints()
562 << endl;
563 if (nNonManifold > 0)
564 {
565 triSurface dupSurf(newFaces, s.patches(), newPoints, true);
566
567 // Create map from dupSurf localPoints to s.localPoints
568 const Map<label> mpm = s.meshPointMap();
569
570 const labelList& dupMp = dupSurf.meshPoints();
571
572 labelList dupPointMap(dupMp.size());
573 forAll(dupMp, pointI)
574 {
575 label dupMeshPointI = dupMp[pointI];
576 label meshPointI = newPointMap[dupMeshPointI];
577 dupPointMap[pointI] = mpm[meshPointI];
578 }
579
580
581 forAll(dupPointMap, pointI)
582 {
583 const point& dupPt = dupSurf.points()[dupMp[pointI]];
584 label sLocalPointI = dupPointMap[pointI];
585 label sMeshPointI = s.meshPoints()[sLocalPointI];
586 const point& sPt = s.points()[sMeshPointI];
587
588 if (mag(dupPt-sPt) > SMALL)
589 {
591 << "dupPt:" << dupPt
592 << " sPt:" << sPt
593 << exit(FatalError);
594 }
595 }
596
597
598 //s.transfer(dupSurf);
599 s = dupSurf;
600 pointMap = labelUIndList(pointMap, dupPointMap)();
601 }
602
603 return nNonManifold;
604}
605
606
607// Find intersections of surf1 by edges of surf2. Return number of degenerate
608// cuts (cuts through faces/edges/points)
609labelPair edgeIntersectionsCGAL
610(
611 const Tree& tree,
612 const triSurface& surf,
613 const pointField& points,
614 edgeIntersections& edgeCuts
615)
616{
617 const edgeList& edges = surf.edges();
618 const labelList& meshPoints = surf.meshPoints();
619
620 //Info<< "Intersecting CGAL surface ..." << endl;
621 List<List<pointIndexHit>> intersections(edges.size());
622 labelListList classifications(edges.size());
623
624 label nPoints = 0;
625 label nSegments = 0;
626
627 std::vector<Segment_intersection> segments;
628 forAll(edges, edgeI)
629 {
630 const edge& e = edges[edgeI];
631
632 const point& a = points[meshPoints[e[0]]];
633 const point& b = points[meshPoints[e[1]]];
634
635 K::Segment_3 segment_query
636 (
637 Point(a.x(), a.y(), a.z()),
638 Point(b.x(), b.y(), b.z())
639 );
640
641 segments.clear();
642 tree.all_intersections(segment_query, std::back_inserter(segments));
643
644
645 for (const Segment_intersection& intersect : segments)
646 {
647 // Get intersection object
648 if
649 (
650 const Point* ptPtr = boost::get<Point>(&(intersect->first))
651 )
652 {
653 point pt
654 (
655 CGAL::to_double(ptPtr->x()),
656 CGAL::to_double(ptPtr->y()),
657 CGAL::to_double(ptPtr->z())
658 );
659
660 #if defined (CGAL_VERSION_NR) && (CGAL_VERSION_NR < 1041400000)
661 Polyhedron::Face_handle f = (intersect->second);
662 #else
663 // 1.14 and later
664 Polyhedron::Face_handle f = (intersect->second).first;
665 #endif
666
667 intersections[edgeI].append
668 (
670 (
671 true,
672 pt,
673 f->index
674 )
675 );
676 // Intersection on edge interior
677 classifications[edgeI].append(-1);
678 ++nPoints;
679 }
680 else if
681 (
682 const Segment* sPtr = boost::get<Segment>(&(intersect->first))
683 )
684 {
685 #if defined (CGAL_VERSION_NR) && (CGAL_VERSION_NR < 1041400000)
686 Polyhedron::Face_handle f = (intersect->second);
687 #else
688 // 1.14 and later
689 Polyhedron::Face_handle f = (intersect->second).first;
690 #endif
691
692 //std::cout
693 // << "intersection object is a segment:" << sPtr->source()
694 // << " " << sPtr->target() << std::endl;
695
696 //std::cout<< "triangle:" << f->index
697 // << " region:" << f->region << std::endl;
698
699 const point source
700 (
701 CGAL::to_double(sPtr->source().x()),
702 CGAL::to_double(sPtr->source().y()),
703 CGAL::to_double(sPtr->source().z())
704 );
705
706 const point target
707 (
708 CGAL::to_double(sPtr->target().x()),
709 CGAL::to_double(sPtr->target().y()),
710 CGAL::to_double(sPtr->target().z())
711 );
712
713 // Make up some intersection point
714 intersections[edgeI].append
715 (
717 (
718 true,
719 0.5*(source+target),
720 f->index
721 )
722 );
723 // Intersection aligned with face. Tbd: enums
724 classifications[edgeI].append(2);
725 ++nSegments;
726 }
727 }
728 }
729
730 edgeCuts = edgeIntersections(intersections, classifications);
731
732 return labelPair(nPoints, nSegments);
733}
734
735
736// Intersect edges of surf1 until no more degenerate intersections. Return
737// number of degenerates
738labelPair edgeIntersectionsAndShuffleCGAL
739(
740 Random& rndGen,
741 const triSurface& surf2,
742 const scalarField& surf1PointTol,
743 triSurface& surf1,
744 edgeIntersections& edgeCuts1
745)
746{
747 //Info<< "Constructing CGAL surface ..." << endl;
749 PolyhedronReader(surf2, p);
750
751
752 //Info<< "Constructing CGAL tree ..." << endl;
753 const Tree tree(p.facets_begin(), p.facets_end(), p);
754
755
756 labelPair cuts(0, 0);
757
758 // Update surface1 points until no longer intersecting
759 pointField surf1Points(surf1.points());
760
761 for (label iter = 0; iter < 5; iter++)
762 {
763 // See which edges of 1 intersect 2
764 Info<< "Determining intersections of " << surf1.nEdges()
765 << " edges with surface of " << label(tree.size()) << " triangles"
766 << endl;
767 cuts = edgeIntersectionsCGAL
768 (
769 tree,
770 surf1,
771 surf1Points,
772 edgeCuts1
773 );
774 Info<< "Determined intersections:" << nl
775 << " edges : " << surf1.nEdges() << nl
776 << " non-degenerate cuts : " << cuts.first() << nl
777 << " degenerate cuts : " << cuts.second() << nl
778 << endl;
779
780 if (cuts.second() == 0)
781 {
782 break;
783 }
784
785 Info<< "Shuffling conflicting points" << endl;
786
787 const labelListList& edgeStat = edgeCuts1.classification();
788 const edgeList& edges = surf1.edges();
789 const labelList& mp = surf1.meshPoints();
790 const point p05(0.5, 0.5, 0.5);
791
792 forAll(edgeStat, edgeI)
793 {
794 const labelList& stat = edgeStat[edgeI];
795 forAll(stat, i)
796 {
797 if (stat[i] == 2)
798 {
799 const edge& e = edges[edgeI];
800 forAll(e, eI)
801 {
802 vector d = rndGen.sample01<vector>() - p05;
803 surf1Points[mp[e[eI]]] += surf1PointTol[e[eI]]*d;
804 }
805 }
806 }
807 }
808 }
809 surf1.movePoints(surf1Points);
810
811 return cuts;
812}
813
814
815// Return map from subSurf.edges() back to surf.edges()
816labelList matchEdges
817(
818 const triSurface& surf,
819 const triSurface& subSurf,
820 const labelList& pointMap
821)
822{
823 if (pointMap.size() != subSurf.nPoints())
824 {
826 << "problem" << exit(FatalError);
827 }
828
829 labelList edgeMap(subSurf.nEdges(), -1);
830
831 const edgeList& edges = surf.edges();
832 const labelListList& pointEdges = surf.pointEdges();
833
834 const edgeList& subEdges = subSurf.edges();
835
836
837 forAll(subEdges, subEdgeI)
838 {
839 const edge& subE = subEdges[subEdgeI];
840
841 // Match points on edge to those on surf
842 label start = pointMap[subE[0]];
843 label end = pointMap[subE[1]];
844 const labelList& pEdges = pointEdges[start];
845 forAll(pEdges, pEdgeI)
846 {
847 label edgeI = pEdges[pEdgeI];
848 const edge& e = edges[edgeI];
849
850 if (e.otherVertex(start) == end)
851 {
852 if (edgeMap[subEdgeI] == -1)
853 {
854 edgeMap[subEdgeI] = edgeI;
855 }
856 else if (edgeMap[subEdgeI] != edgeI)
857 {
859 << subE << " points:"
860 << subE.line(subSurf.localPoints())
861 << " matches to " << edgeI
862 << " points:" << e.line(surf.localPoints())
863 << " but also " << edgeMap[subEdgeI]
864 << " points:"
865 << edges[edgeMap[subEdgeI]].line(surf.localPoints())
866 << exit(FatalError);
867 }
868 break;
869 }
870 }
871
872 if (edgeMap[subEdgeI] == -1)
873 {
875 << subE << " at:" << subSurf.localPoints()[subE[0]]
876 << subSurf.localPoints()[subE[1]]
877 << exit(FatalError);
878 }
879 }
880
881 return edgeMap;
882}
883
884
885void calcEdgeCutsCGAL
886(
887 triSurface& surf1,
888 triSurface& surf2,
889 const bool perturb,
890 edgeIntersections& edgeCuts1,
891 edgeIntersections& edgeCuts2
892)
893{
894 if (!perturb)
895 {
896 // See which edges of 1 intersect 2
897 {
898 Info<< "Intersect surface 1 edges with surface 2:" << nl;
899 Info<< " constructing CGAL surface ..." << endl;
901 PolyhedronReader(surf2, p);
902
903 Info<< " constructing CGAL tree ..." << endl;
904 const Tree tree(p.facets_begin(), p.facets_end(), p);
905
906 edgeIntersectionsCGAL
907 (
908 tree,
909 surf1,
910 surf1.points(),
911 edgeCuts1
912 );
913 }
914 // See which edges of 2 intersect 1
915 {
916 Info<< "Intersect surface 2 edges with surface 1:" << nl;
917 Info<< " constructing CGAL surface ..." << endl;
919 PolyhedronReader(surf1, p);
920
921 Info<< " constructing CGAL tree ..." << endl;
922 const Tree tree(p.facets_begin(), p.facets_end(), p);
923
924 edgeIntersectionsCGAL
925 (
926 tree,
927 surf2,
928 surf2.points(),
929 edgeCuts2
930 );
931 }
932 Info<< endl;
933 }
934 else
935 {
936 const scalarField surf1PointTol
937 (
938 max(1e-8*edgeIntersections::minEdgeLength(surf1), SMALL)
939 );
940 const scalarField surf2PointTol
941 (
942 max(1e-8*edgeIntersections::minEdgeLength(surf2), SMALL)
943 );
944
945
946 Random rndGen(0);
947
948 labelPair cuts1;
949 labelPair cuts2;
950
951 for (label iter = 0; iter < 10; iter++)
952 {
953 // Find intersections of surf1 edges with surf2 triangles
954 cuts1 = edgeIntersectionsAndShuffleCGAL
955 (
956 rndGen,
957 surf2,
958 surf1PointTol,
959 surf1,
960 edgeCuts1
961 );
962
963 // Find intersections of surf2 edges with surf1 triangles
964 cuts2 = edgeIntersectionsAndShuffleCGAL
965 (
966 rndGen,
967 surf1,
968 surf2PointTol,
969 surf2,
970 edgeCuts2
971 );
972
973 if (cuts1.second() + cuts2.second() == 0)
974 {
975 break;
976 }
977 }
978 }
979}
980
981
982void calcEdgeCutsBitsCGAL
983(
984 triSurface& surf1,
985 triSurface& surf2,
986 const bool perturb,
987 edgeIntersections& edgeCuts1,
988 edgeIntersections& edgeCuts2
989)
990{
991 label nZones1 = 1;
992 labelList zone1;
993 {
994 labelHashSet orientationEdge(surf1.size()/1000);
995 PatchTools::checkOrientation(surf1, false, &orientationEdge);
996 nZones1 = PatchTools::markZones(surf1, orientationEdge, zone1);
997
998 Info<< "Surface triangles " << surf1.size()
999 << " number of zones (orientation or non-manifold):"
1000 << nZones1 << endl;
1001 }
1002
1003 label nZones2 = 1;
1004 labelList zone2;
1005 {
1006 labelHashSet orientationEdge(surf2.size()/1000);
1007 PatchTools::checkOrientation(surf2, false, &orientationEdge);
1008 nZones2 = PatchTools::markZones(surf2, orientationEdge, zone2);
1009
1010 Info<< "Surface triangles " << surf2.size()
1011 << " number of zones (orientation or non-manifold):"
1012 << nZones2 << endl;
1013 }
1014
1015
1016 if (nZones1 == 1 && nZones2 == 1)
1017 {
1018 calcEdgeCutsCGAL
1019 (
1020 surf1,
1021 surf2,
1022 perturb,
1023 edgeCuts1,
1024 edgeCuts2
1025 );
1026 }
1027 else
1028 {
1029 edgeCuts1 = edgeIntersections
1030 (
1032 labelListList(surf1.nEdges())
1033 );
1034 edgeCuts2 = edgeIntersections
1035 (
1037 labelListList(surf2.nEdges())
1038 );
1039
1040
1041 for (label zone1I = 0; zone1I < nZones1; zone1I++)
1042 {
1043 // Generate sub surface for zone1I
1044 boolList includeMap1(surf1.size(), false);
1045
1046 forAll(zone1, faceI)
1047 {
1048 if (zone1[faceI] == zone1I)
1049 {
1050 includeMap1[faceI] = true;
1051 }
1052 }
1053
1054 // Subset. Map from local points on subset to local points on
1055 // original
1056 labelList pointMap1;
1057 labelList faceMap1;
1058 triSurface subSurf1
1059 (
1060 surf1.subsetMesh
1061 (
1062 includeMap1,
1063 pointMap1,
1064 faceMap1
1065 )
1066 );
1067
1068 // Split non-manifold points; update pointMap
1069 dupNonManifoldPoints(subSurf1, pointMap1);
1070
1071 const boundBox subBb1
1072 (
1073 subSurf1.points(),
1074 subSurf1.meshPoints(),
1075 false
1076 );
1077
1078 const labelList edgeMap1
1079 (
1080 matchEdges
1081 (
1082 surf1,
1083 subSurf1,
1084 pointMap1
1085 )
1086 );
1087
1088
1089 for (label zone2I = 0; zone2I < nZones2; zone2I++)
1090 {
1091 // Generate sub surface for zone2I
1092 boolList includeMap2(surf2.size(), false);
1093
1094 forAll(zone2, faceI)
1095 {
1096 if (zone2[faceI] == zone2I)
1097 {
1098 includeMap2[faceI] = true;
1099 }
1100 }
1101
1102 labelList pointMap2;
1103 labelList faceMap2;
1104 triSurface subSurf2
1105 (
1106 surf2.subsetMesh
1107 (
1108 includeMap2,
1109 pointMap2,
1110 faceMap2
1111 )
1112 );
1113
1114
1115 const boundBox subBb2
1116 (
1117 subSurf2.points(),
1118 subSurf2.meshPoints(),
1119 false
1120 );
1121
1122 // Short-circuit expensive calculations
1123 if (!subBb2.overlaps(subBb1))
1124 {
1125 continue;
1126 }
1127
1128
1129 // Split non-manifold points; update pointMap
1130 dupNonManifoldPoints(subSurf2, pointMap2);
1131
1132 const labelList edgeMap2
1133 (
1134 matchEdges
1135 (
1136 surf2,
1137 subSurf2,
1138 pointMap2
1139 )
1140 );
1141
1142
1143 // Do cuts
1144 edgeIntersections subEdgeCuts1;
1145 edgeIntersections subEdgeCuts2;
1146 calcEdgeCutsCGAL
1147 (
1148 subSurf1,
1149 subSurf2,
1150 perturb,
1151 subEdgeCuts1,
1152 subEdgeCuts2
1153 );
1154
1155 // Move original surface
1156 {
1157 pointField points2(surf2.points());
1158 forAll(pointMap2, i)
1159 {
1160 label subMeshPointI = subSurf2.meshPoints()[i];
1161 label meshPointI = surf2.meshPoints()[pointMap2[i]];
1162 points2[meshPointI] = subSurf2.points()[subMeshPointI];
1163 }
1164 surf2.movePoints(points2);
1165 }
1166
1167 // Insert into main structure
1168 edgeCuts1.merge(subEdgeCuts1, edgeMap1, faceMap2);
1169 edgeCuts2.merge(subEdgeCuts2, edgeMap2, faceMap1);
1170 }
1171
1172
1173 // Move original surface
1174 {
1175 pointField points1(surf1.points());
1176 forAll(pointMap1, i)
1177 {
1178 label subMeshPointI = subSurf1.meshPoints()[i];
1179 label meshPointI = surf1.meshPoints()[pointMap1[i]];
1180 points1[meshPointI] = subSurf1.points()[subMeshPointI];
1181 }
1182 surf1.movePoints(points1);
1183 }
1184 }
1185 }
1186}
1187
1188#endif // NO_CGAL
1189
1190
1191//void calcFeaturePoints(const pointField& points, const edgeList& edges)
1192//{
1193// edgeMesh eMesh(points, edges);
1194//
1195// const labelListList& pointEdges = eMesh.pointEdges();
1196//
1197//
1198// // Get total number of feature points
1199// label nFeaturePoints = 0;
1200// forAll(pointEdges, pI)
1201// {
1202// const labelList& pEdges = pointEdges[pI];
1203//
1204// if (pEdges.size() == 1)
1205// {
1206// nFeaturePoints++;
1207// }
1208// }
1209//
1210//
1211// // Calculate addressing from feature point to cut point and cut edge
1212// labelList featurePointToCutPoint(nFeaturePoints);
1213// labelList featurePointToCutEdge(nFeaturePoints);
1214//
1215// label nFeatPts = 0;
1216// forAll(pointEdges, pI)
1217// {
1218// const labelList& pEdges = pointEdges[pI];
1219//
1220// if (pEdges.size() == 1)
1221// {
1222// featurePointToCutPoint[nFeatPts] = pI;
1223// featurePointToCutEdge[nFeatPts] = pEdges[0];
1224// nFeatPts++;
1225// }
1226// }
1227//
1228//
1229//
1230// label concaveStart = 0;
1231// label mixedStart = 0;
1232// label nonFeatureStart = nFeaturePoints;
1233//
1234//
1235// labelListList featurePointNormals(nFeaturePoints);
1236// labelListList featurePointEdges(nFeaturePoints);
1237// labelList regionEdges;
1238//}
1239
1240
1242(
1243 const IOobject& io,
1244 const booleanSurface::booleanOpType action,
1245 const bool surf1Baffle,
1246 const bool surf2Baffle,
1247 const bool invertedSpace,
1248 const triSurface& surf1,
1249 const edgeIntersections& edgeCuts1,
1250 const triSurface& surf2,
1251 const edgeIntersections& edgeCuts2
1252)
1253{
1254 // Determine intersection edges from the edge cuts
1256 (
1257 surf1,
1258 edgeCuts1,
1259 surf2,
1260 edgeCuts2
1261 );
1262
1263 label nFeatEds = inter.cutEdges().size();
1264
1265 DynamicList<vector> normals(2*nFeatEds);
1266 vectorField edgeDirections(nFeatEds, Zero);
1268 (
1269 2*nFeatEds
1270 );
1271 List<DynamicList<label>> edgeNormals(nFeatEds);
1272 List<DynamicList<label>> normalDirections(nFeatEds);
1273
1274
1275 const triSurface& s1 = surf1;
1276 const triSurface& s2 = surf2;
1277
1278 forAllConstIters(inter.facePairToEdgeId(), iter)
1279 {
1280 const labelPair& facePair = iter.key();
1281 const label cutEdgeI = iter.val();
1282
1283 const edge& fE = inter.cutEdges()[cutEdgeI];
1284
1285 const vector& norm1 = s1.faceNormals()[facePair.first()];
1286 const vector& norm2 = s2.faceNormals()[facePair.second()];
1287
1288 DynamicList<label>& eNormals = edgeNormals[cutEdgeI];
1289 DynamicList<label>& nDirections = normalDirections[cutEdgeI];
1290
1291 edgeDirections[cutEdgeI] = fE.vec(inter.cutPoints());
1292
1293 normals.append(norm1);
1294 eNormals.append(normals.size() - 1);
1295
1296 if (surf1Baffle)
1297 {
1298 normalVolumeTypes.append(extendedFeatureEdgeMesh::BOTH);
1299
1300 nDirections.append(1);
1301 }
1302 else
1303 {
1304 normalVolumeTypes.append(extendedFeatureEdgeMesh::INSIDE);
1305 nDirections.append
1306 (
1307 calcNormalDirection
1308 (
1309 norm1,
1310 norm2,
1311 edgeDirections[cutEdgeI],
1312 s1[facePair.first()].centre(s1.points()),
1313 inter.cutPoints()[fE.start()]
1314 )
1315 );
1316 }
1317
1318 normals.append(norm2);
1319 eNormals.append(normals.size() - 1);
1320
1321 if (surf2Baffle)
1322 {
1323 normalVolumeTypes.append(extendedFeatureEdgeMesh::BOTH);
1324
1325 nDirections.append(1);
1326 }
1327 else
1328 {
1329 normalVolumeTypes.append(extendedFeatureEdgeMesh::INSIDE);
1330
1331 nDirections.append
1332 (
1333 calcNormalDirection
1334 (
1335 norm2,
1336 norm1,
1337 edgeDirections[cutEdgeI],
1338 s2[facePair.second()].centre(s2.points()),
1339 inter.cutPoints()[fE.start()]
1340 )
1341 );
1342 }
1343
1344
1345 if (surf1Baffle)
1346 {
1347 normals.append(norm2);
1348
1349 if (surf2Baffle)
1350 {
1351 normalVolumeTypes.append(extendedFeatureEdgeMesh::BOTH);
1352
1353 nDirections.append(1);
1354 }
1355 else
1356 {
1357 normalVolumeTypes.append(extendedFeatureEdgeMesh::INSIDE);
1358
1359 nDirections.append
1360 (
1361 calcNormalDirection
1362 (
1363 norm2,
1364 norm1,
1365 edgeDirections[cutEdgeI],
1366 s2[facePair.second()].centre(s2.points()),
1367 inter.cutPoints()[fE.start()]
1368 )
1369 );
1370 }
1371
1372 eNormals.append(normals.size() - 1);
1373 }
1374
1375 if (surf2Baffle)
1376 {
1377 normals.append(norm1);
1378
1379 if (surf1Baffle)
1380 {
1381 normalVolumeTypes.append(extendedFeatureEdgeMesh::BOTH);
1382
1383 nDirections.append(1);
1384 }
1385 else
1386 {
1387 normalVolumeTypes.append(extendedFeatureEdgeMesh::INSIDE);
1388
1389 nDirections.append
1390 (
1391 calcNormalDirection
1392 (
1393 norm1,
1394 norm2,
1395 edgeDirections[cutEdgeI],
1396 s1[facePair.first()].centre(s1.points()),
1397 inter.cutPoints()[fE.start()]
1398 )
1399 );
1400 }
1401
1402 eNormals.append(normals.size() - 1);
1403 }
1404 }
1405
1406
1407 label internalStart = -1;
1408 label nIntOrExt = 0;
1409 label nFlat = 0;
1410 label nOpen = 0;
1411 label nMultiple = 0;
1412
1413 forAll(edgeNormals, eI)
1414 {
1415 label nEdNorms = edgeNormals[eI].size();
1416
1417 if (nEdNorms == 1)
1418 {
1419 nOpen++;
1420 }
1421 else if (nEdNorms == 2)
1422 {
1423 const vector& n0(normals[edgeNormals[eI][0]]);
1424 const vector& n1(normals[edgeNormals[eI][1]]);
1425
1426 if ((n0 & n1) > extendedFeatureEdgeMesh::cosNormalAngleTol_)
1427 {
1428 nFlat++;
1429 }
1430 else
1431 {
1432 nIntOrExt++;
1433 }
1434 }
1435 else if (nEdNorms > 2)
1436 {
1437 nMultiple++;
1438 }
1439 }
1440
1441 if (action == booleanSurface::UNION)
1442 {
1443 if (!invertedSpace)
1444 {
1445 // All edges are internal
1446 internalStart = 0;
1447 }
1448 else
1449 {
1450 // All edges are external
1451 internalStart = nIntOrExt;
1452 }
1453 }
1454 else if (action == booleanSurface::INTERSECTION)
1455 {
1456 if (!invertedSpace)
1457 {
1458 // All edges are external
1459 internalStart = nIntOrExt;
1460 }
1461 else
1462 {
1463 // All edges are internal
1464 internalStart = 0;
1465 }
1466 }
1467 else if (action == booleanSurface::DIFFERENCE)
1468 {
1469 // All edges are external
1470 internalStart = nIntOrExt;
1471 }
1472 else
1473 {
1475 << "Unsupported booleanSurface:booleanOpType and space "
1476 << action << " " << invertedSpace
1477 << abort(FatalError);
1478 }
1479
1480 // There are no feature points supported by surfaceIntersection
1481 // Flat, open or multiple edges are assumed to be impossible
1482 // Region edges are not explicitly supported by surfaceIntersection
1483
1484 vectorField normalsTmp(normals);
1486 (
1487 normalVolumeTypes
1488 );
1489 labelListList edgeNormalsTmp(edgeNormals.size());
1490 forAll(edgeNormalsTmp, i)
1491 {
1492 edgeNormalsTmp[i] = edgeNormals[i];
1493 }
1494 labelListList normalDirectionsTmp(normalDirections.size());
1495 forAll(normalDirectionsTmp, i)
1496 {
1497 normalDirectionsTmp[i] = normalDirections[i];
1498 }
1499
1500 //calcFeaturePoints(inter.cutPoints(), inter.cutEdges());
1501
1503 (
1504 io,
1505 inter.cutPoints(),
1506 inter.cutEdges(),
1507
1508 0, // concaveStart,
1509 0, // mixedStart,
1510 0, // nonFeatureStart,
1511
1512 internalStart, // internalStart,
1513 nIntOrExt, // flatStart,
1514 nIntOrExt + nFlat, // openStart,
1515 nIntOrExt + nFlat + nOpen, // multipleStart,
1516
1517 normalsTmp,
1518 normalVolumeTypesTmp,
1519 edgeDirections,
1520 normalDirectionsTmp,
1521 edgeNormalsTmp,
1522
1523 labelListList(), // featurePointNormals,
1524 labelListList(), // featurePointEdges,
1525 labelList() // regionEdges
1526 );
1527}
1528
1529int main(int argc, char *argv[])
1530{
1531 argList::addNote
1532 (
1533 "Generates a extendedFeatureEdgeMesh for the interface created by"
1534 " a boolean operation on two surfaces."
1535 #ifndef NO_CGAL
1536 " [Compiled with CGAL]"
1537 #else
1538 " [Compiled without CGAL]"
1539 #endif
1540 );
1541
1542 argList::noParallel();
1543 argList::addArgument
1544 (
1545 "action",
1546 "One of (intersection | union | difference)"
1547 );
1548 argList::addArgument("surface1", "The input surface file 1");
1549 argList::addArgument("surface2", "The input surface file 2");
1550
1551 argList::addOption
1552 (
1553 "scale",
1554 "factor",
1555 "Geometry scaling factor (both surfaces)"
1556 );
1557 argList::addBoolOption
1558 (
1559 "surf1Baffle",
1560 "Mark surface 1 as a baffle"
1561 );
1562
1563 argList::addBoolOption
1564 (
1565 "surf2Baffle",
1566 "Mark surface 2 as a baffle"
1567 );
1568
1569 argList::addBoolOption
1570 (
1571 "perturb",
1572 "Perturb surface points to escape degenerate intersections"
1573 );
1574
1575 argList::addBoolOption
1576 (
1577 "no-cgal",
1578 #ifndef NO_CGAL
1579 "Do not use CGAL algorithms"
1580 #else
1581 "Ignored, compiled without CGAL"
1582 #endif
1583 );
1584
1585 argList::addBoolOption
1586 (
1587 "invertedSpace",
1588 "Do the surfaces have inverted space orientation, "
1589 "i.e. a point at infinity is considered inside. "
1590 "This is only sensible for union and intersection."
1591 );
1592
1593 argList::addOption
1594 (
1595 "trim",
1596 "((surface1 volumeType) .. (surfaceN volumeType))",
1597 "Trim resulting intersection with additional surfaces;"
1598 " volumeType is 'inside' (keep (parts of) edges that are inside)"
1599 ", 'outside' (keep (parts of) edges that are outside) or"
1600 " 'mixed' (keep all)"
1601 );
1602
1603 #include "setRootCase.H"
1604 #include "createTime.H"
1605
1606 const word action(args[1]);
1607
1608 const Enum<booleanSurface::booleanOpType> validActions
1609 {
1610 { booleanSurface::INTERSECTION, "intersection" },
1611 { booleanSurface::UNION, "union" },
1612 { booleanSurface::DIFFERENCE, "difference" }
1613 };
1614
1615 if (!validActions.found(action))
1616 {
1618 << "Unsupported action " << action << endl
1619 << "Supported actions:" << validActions << nl
1620 << abort(FatalError);
1621 }
1622
1623
1624 List<Pair<word>> surfaceAndSide;
1625 if (args.readIfPresent("trim", surfaceAndSide))
1626 {
1627 Info<< "Trimming edges with " << surfaceAndSide << endl;
1628 }
1629
1630
1631 // Scale factor for both surfaces:
1632 const scalar scaleFactor = args.getOrDefault<scalar>("scale", -1);
1633
1634 const word surf1Name(args[2]);
1635 Info<< "Reading surface " << surf1Name << endl;
1636 triSurfaceMesh surf1
1637 (
1638 IOobject
1639 (
1640 surf1Name,
1641 runTime.constant(),
1642 triSurfaceMesh::meshSubDir,
1643 runTime
1644 )
1645 );
1646 if (scaleFactor > 0)
1647 {
1648 Info<< "Scaling : " << scaleFactor << nl;
1649 surf1.scalePoints(scaleFactor);
1650 }
1651
1652 Info<< surf1Name << " statistics:" << endl;
1653 surf1.writeStats(Info);
1654 Info<< endl;
1655
1656 const word surf2Name(args[3]);
1657 Info<< "Reading surface " << surf2Name << endl;
1658 triSurfaceMesh surf2
1659 (
1660 IOobject
1661 (
1662 surf2Name,
1663 runTime.constant(),
1664 triSurfaceMesh::meshSubDir,
1665 runTime
1666 )
1667 );
1668 if (scaleFactor > 0)
1669 {
1670 Info<< "Scaling : " << scaleFactor << nl;
1671 surf2.scalePoints(scaleFactor);
1672 }
1673
1674 Info<< surf2Name << " statistics:" << endl;
1675 surf2.writeStats(Info);
1676 Info<< endl;
1677
1678 const bool surf1Baffle = args.found("surf1Baffle");
1679 const bool surf2Baffle = args.found("surf2Baffle");
1680
1681 edgeIntersections edgeCuts1;
1682 edgeIntersections edgeCuts2;
1683
1684 const bool invertedSpace = args.found("invertedSpace");
1685
1686 if (invertedSpace && validActions[action] == booleanSurface::DIFFERENCE)
1687 {
1689 << "Inverted space only makes sense for union or intersection."
1690 << exit(FatalError);
1691 }
1692
1693
1694 // Calculate where edges are cut by the other surface
1695#ifndef NO_CGAL
1696 if (!args.found("no-cgal"))
1697 {
1698 calcEdgeCutsBitsCGAL
1699 (
1700 surf1,
1701 surf2,
1702 args.found("perturb"),
1703 edgeCuts1,
1704 edgeCuts2
1705 );
1706 }
1707 else
1708#endif // NO_CGAL
1709 {
1710 calcEdgeCuts
1711 (
1712 surf1,
1713 surf2,
1714 args.found("perturb"),
1715 edgeCuts1,
1716 edgeCuts2
1717 );
1718 }
1719
1720 const fileName sFeatFileName
1721 (
1722 fileName(surf1Name).nameLessExt()
1723 + "_"
1724 + fileName(surf2Name).nameLessExt()
1725 + "_"
1726 + action
1727 );
1728
1730 (
1731 createEdgeMesh
1732 (
1733 IOobject
1734 (
1735 sFeatFileName + ".extendedFeatureEdgeMesh",
1736 runTime.constant(),
1737 "extendedFeatureEdgeMesh",
1738 runTime,
1739 IOobject::NO_READ,
1740 IOobject::NO_WRITE
1741 ),
1742 booleanSurface::booleanOpTypeNames[action],
1743 surf1Baffle,
1744 surf2Baffle,
1745 invertedSpace,
1746 surf1,
1747 edgeCuts1,
1748 surf2,
1749 edgeCuts2
1750 )
1751 );
1752
1753
1754 // Trim intersections
1755 forAll(surfaceAndSide, trimI)
1756 {
1757 const word& trimName = surfaceAndSide[trimI].first();
1758 const volumeType trimType
1759 (
1760 volumeType::names[surfaceAndSide[trimI].second()]
1761 );
1762
1763 Info<< "Reading trim surface " << trimName << endl;
1764 const triSurfaceMesh trimSurf
1765 (
1766 IOobject
1767 (
1768 trimName,
1769 runTime.constant(),
1770 triSurfaceMesh::meshSubDir,
1771 runTime,
1772 IOobject::MUST_READ,
1773 IOobject::NO_WRITE
1774 )
1775 );
1776
1777 Info<< trimName << " statistics:" << endl;
1778 trimSurf.writeStats(Info);
1779 Info<< endl;
1780
1781 labelList pointMap;
1782 labelList edgeMap;
1783 feMeshPtr().trim
1784 (
1785 trimSurf,
1786 trimType,
1787 pointMap,
1788 edgeMap
1789 );
1790 }
1791
1792
1793 const extendedFeatureEdgeMesh& feMesh = feMeshPtr();
1794
1795 feMesh.writeStats(Info);
1796 Info<< nl << "Writing extendedFeatureEdgeMesh to "
1797 << feMesh.objectRelPath() << nl
1798 << endl;
1799 feMesh.write();
1800 feMesh.writeObj(feMesh.path()/sFeatFileName);
1801
1802 {
1803 // Write a featureEdgeMesh for backwards compatibility
1804 featureEdgeMesh bfeMesh
1805 (
1806 IOobject
1807 (
1808 sFeatFileName + ".eMesh", // name
1809 runTime.constant(), // instance
1810 triSurfaceMesh::meshSubDir,
1811 runTime, // registry
1812 IOobject::NO_READ,
1813 IOobject::NO_WRITE,
1814 false
1815 ),
1816 feMesh.points(),
1817 feMesh.edges()
1818 );
1819
1820 Info<< nl << "Writing featureEdgeMesh to "
1821 << bfeMesh.objectRelPath() << endl;
1822
1823 bfeMesh.regIOobject::write();
1824 }
1825
1826 Info << nl << "End\n" << endl;
1827
1828 return 0;
1829}
1830
1831
1832// ************************************************************************* //
CGAL data structures used for triSurface handling.
CGAL::Segment_3< K > Segment
CGAL::Polyhedron_3< K, My_items > Polyhedron
CGAL::Point_3< K > Point
Y[inertIndex] max(0.0)
Dynamically sized Field.
Definition: DynamicField.H:65
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:72
void append(const T &val)
Copy append an element to the end of this list.
Definition: DynamicListI.H:503
Enum is a wrapper around a list of names/values that represent particular enumeration (or int) values...
Definition: Enum.H:61
T & first() noexcept
The first element of the list, position [0].
Definition: FixedListI.H:207
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:170
fileName objectRelPath() const
The object path relative to the root.
Definition: IOobject.C:558
fileName path() const
The complete path.
Definition: IOobject.C:524
A HashTable to objects of type <T> with a label key.
Definition: Map.H:60
const T & second() const noexcept
Return second element, which is also the last element.
Definition: PairI.H:120
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
Definition: PointIndexHit.H:66
label nEdges() const
Number of edges in patch.
label nPoints() const
Number of points supporting patch faces.
const labelListList & pointEdges() const
Return point-edge addressing.
const edgeList & edges() const
Return list of edges, address into LOCAL point list.
const labelList & meshPoints() const
Return labelList of mesh points in patch.
const Field< point_type > & localPoints() const
Return pointField of points in patch.
const Field< point_type > & faceNormals() const
Return face unit normals for patch.
const Field< point_type > & points() const noexcept
Return reference to global points.
Random number generator.
Definition: Random.H:60
A List with indirect addressing. Like IndirectList but does not store addressing.
Definition: IndirectList.H:79
T & first()
Return the first element of the list.
Definition: UListI.H:202
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
const Cmpt & z() const
Access to the vector z component.
Definition: VectorI.H:85
const Cmpt & y() const
Access to the vector y component.
Definition: VectorI.H:79
const Cmpt & x() const
Access to the vector x component.
Definition: VectorI.H:73
bool found(const word &optName) const
Return true if the named option is found.
Definition: argListI.H:178
bool readIfPresent(const word &optName, T &val) const
Read a value from the named option if present.
Definition: argListI.H:323
T getOrDefault(const word &optName, const T &deflt) const
Get a value from the named option if present, or return default.
Definition: argListI.H:307
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: autoPtr.H:66
booleanOpType
Enumeration listing the possible volume operator types.
A bounding box defined in terms of min/max extrema points.
Definition: boundBox.H:64
Holder of intersections of edges of a surface with another surface. Optionally shuffles around points...
const labelListList & classification() const
For every intersection the classification status.
void merge(const edgeIntersections &, const labelList &edgeMap, const labelList &faceMap, const bool merge=true)
Merge (or override) edge intersection for a subset.
label removeDegenerates(const label nIters, const triSurface &surf1, const triSurfaceSearch &query2, const scalarField &surf1PointTol, pointField &points1)
Resolve ties. Shuffles points so all edge - face intersections.
const pointField & points() const noexcept
Return points.
Definition: edgeMeshI.H:99
const edgeList & edges() const noexcept
Return edges.
Definition: edgeMeshI.H:105
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:66
vector vec(const UList< point > &pts) const
Return the vector (end - start)
Definition: edgeI.H:417
linePointRef line(const UList< point > &pts) const
Return edge line.
Definition: edgeI.H:456
label start() const
Return start (first) vertex label.
Definition: edgeI.H:95
void writeObj(const fileName &prefix) const
Write all components of the extendedEdgeMesh as obj files.
virtual void writeStats(Ostream &os) const
Dump some information.
virtual bool write(const bool valid=true) const
Give precedence to the regIOobject write.
A class for handling file names.
Definition: fileName.H:76
A triFace with additional (region) index.
Definition: labelledTri.H:60
Basic surface-surface intersection description. Constructed from two surfaces it creates a descriptio...
IOoject and searching on triSurface.
Helper class to search on triSurface.
Triangulated surface description with patch information.
Definition: triSurface.H:79
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
virtual void movePoints(const pointField &pts)
Move points.
Definition: triSurface.C:612
virtual void scalePoints(const scalar scaleFactor)
Scale points. A non-positive factor is ignored.
Definition: triSurface.C:638
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
An enumeration wrapper for classification of a location as being inside/outside of a volume.
Definition: volumeType.H:61
A class for handling words, derived from Foam::string.
Definition: word.H:68
volScalarField & p
engineTime & runTime
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
const pointField & points
label nPoints
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
label nZones
IOobject io("surfaceFilmProperties", mesh.time().constant(), mesh, IOobject::READ_IF_PRESENT, IOobject::NO_WRITE, false)
const dimensionedScalar mp
Proton mass.
Namespace for OpenFOAM.
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i)
Definition: labelList.C:38
quaternion normalised(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
Definition: quaternionI.H:680
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
errorManip< error > abort(error &err)
Definition: errorManip.H:144
void cross(FieldField< Field1, typename crossProduct< Type1, Type2 >::type > &f, const FieldField< Field1, Type1 > &f1, const FieldField< Field2, Type2 > &f2)
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
constexpr auto end(C &c) -> decltype(c.end())
Return iterator to the end of the container c.
Definition: stdFoam.H:158
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints, -1);const cellModel &hex=cellModel::ref(cellModel::HEX);labelList hexLabels(8);label activeCells=0;labelList pointMap(nPoints);forAll(pointMap, i){ pointMap[i]=i;}for(label i=0;i< nPoints;i++){ if(f[i] > 0.0) { hexLabels[0]=i;hexLabels[1]=i1tab[i];hexLabels[2]=i3tab[i1tab[i]];hexLabels[3]=i3tab[i];hexLabels[4]=i8tab[i];hexLabels[5]=i1tab[i8tab[i]];hexLabels[6]=i3tab[i1tab[i8tab[i]]];hexLabels[7]=i3tab[i8tab[i]];cellShapes[activeCells].reset(hex, hexLabels);edgeList edges=cellShapes[activeCells].edges();forAll(edges, ei) { if(edges[ei].mag(points)< SMALL) { label start=pointMap[edges[ei].start()];while(start !=pointMap[start]) { start=pointMap[start];} label end=pointMap[edges[ei].end()];while(end !=pointMap[end]) { end=pointMap[end];} label minLabel=min(start, end);pointMap[start]=pointMap[end]=minLabel;} } cellZoning[activeCells]=idreg[i];activeCells++;}}cellShapes.setSize(activeCells);cellZoning.setSize(activeCells);forAll(cellShapes, celli){ cellShape &cs=cellShapes[celli];forAll(cs, i) { cs[i]=pointMap[cs[i]];} cs.collapse();}label bcIDs[11]={-1, 0, 2, 4, -1, 5, -1, 6, 7, 8, 9};const label nBCs=12;const word *kivaPatchTypes[nBCs]={ &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &symmetryPolyPatch::typeName, &wedgePolyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &symmetryPolyPatch::typeName, &oldCyclicPolyPatch::typeName};enum patchTypeNames{ PISTON, VALVE, LINER, CYLINDERHEAD, AXIS, WEDGE, INFLOW, OUTFLOW, PRESIN, PRESOUT, SYMMETRYPLANE, CYCLIC};const char *kivaPatchNames[nBCs]={ "piston", "valve", "liner", "cylinderHead", "axis", "wedge", "inflow", "outflow", "presin", "presout", "symmetryPlane", "cyclic"};List< SLList< face > > pFaces[nBCs]
Definition: readKivaGrid.H:235
labelList f(nPoints)
Foam::argList args(argc, argv)
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
#define forAllConstIters(container, iter)
Iterate across all elements of the container object with const access.
Definition: stdFoam.H:278
Random rndGen
Definition: createFields.H:23