@@ -7,62 +7,59 @@ Extract visible tetrahedra - those intersecting with the planes
77Return corresponding points and facets for each region for drawing as mesh (Makie,MeshCat)
88or trisurf (pyplot)
99"""
10- function extract_visible_cells3D (coord,cellnodes,cellregions,nregions,xyzcut;
11- primepoints= zeros (0 ,0 ),Tp = SVector{3 ,Float32},Tf = SVector{3 ,Int32})
12- all_lt= ones (Bool,3 )
13- all_gt= ones (Bool,3 )
10+ function extract_visible_cells3D (coord, cellnodes, cellregions, nregions, xyzcut;
11+ primepoints = zeros (0 , 0 ), Tp = SVector{3 , Float32}, Tf = SVector{3 , Int32})
12+ all_lt = ones (Bool, 3 )
13+ all_gt = ones (Bool, 3 )
1414
15- function take (coord,simplex,xyzcut,all_lt,all_gt)
16- for idim= 1 : 3
17- all_lt[idim]= true
18- all_gt[idim]= true
19- for inode= 1 : 4
20- c= coord[idim,simplex[inode]]- xyzcut[idim]
21- all_lt[idim]= all_lt[idim] && (c< 0.0 )
22- all_gt[idim]= all_gt[idim] && (c> 0.0 )
15+ function take (coord, simplex, xyzcut, all_lt, all_gt)
16+ for idim = 1 : 3
17+ all_lt[idim] = true
18+ all_gt[idim] = true
19+ for inode = 1 : 4
20+ c = coord[idim, simplex[inode]] - xyzcut[idim]
21+ all_lt[idim] = all_lt[idim] && (c < 0.0 )
22+ all_gt[idim] = all_gt[idim] && (c > 0.0 )
2323 end
2424 end
25- tke= false
26- tke= tke || (! all_lt[1 ]) && (! all_gt[1 ]) && (! all_gt[2 ]) && (! all_gt[3 ])
27- tke= tke || (! all_lt[2 ]) && (! all_gt[2 ]) && (! all_gt[1 ]) && (! all_gt[3 ])
28- tke= tke || (! all_lt[3 ]) && (! all_gt[3 ]) && (! all_gt[1 ]) && (! all_gt[2 ])
25+ tke = false
26+ tke = tke || (! all_lt[1 ]) && (! all_gt[1 ]) && (! all_gt[2 ]) && (! all_gt[3 ])
27+ tke = tke || (! all_lt[2 ]) && (! all_gt[2 ]) && (! all_gt[1 ]) && (! all_gt[3 ])
28+ tke = tke || (! all_lt[3 ]) && (! all_gt[3 ]) && (! all_gt[1 ]) && (! all_gt[2 ])
2929 end
3030
31-
32- faces= [Vector {Tf} (undef,0 ) for iregion= 1 : nregions]
33- points= [Vector {Tp} (undef,0 ) for iregion= 1 : nregions]
34-
35- for iregion= 1 : nregions
36- for iprime= 1 : size (primepoints,2 )
37- @views push! (points[iregion],Tp (primepoints[:,iprime]))
31+ faces = [Vector {Tf} (undef, 0 ) for iregion = 1 : nregions]
32+ points = [Vector {Tp} (undef, 0 ) for iregion = 1 : nregions]
33+
34+ for iregion = 1 : nregions
35+ for iprime = 1 : size (primepoints, 2 )
36+ @views push! (points[iregion], Tp (primepoints[:, iprime]))
3837 end
3938 end
40- tet= zeros (Int32,4 )
41-
42- for itet= 1 : size (cellnodes,2 )
43- iregion= cellregions[itet]
44- for i= 1 : 4
45- tet[i]= cellnodes[i,itet]
39+ tet = zeros (Int32, 4 )
40+
41+ for itet = 1 : size (cellnodes, 2 )
42+ iregion = cellregions[itet]
43+ for i = 1 : 4
44+ tet[i] = cellnodes[i, itet]
4645 end
47- if take (coord,tet,xyzcut,all_lt,all_gt)
48- npts= size (points[iregion],1 )
46+ if take (coord, tet, xyzcut, all_lt, all_gt)
47+ npts = size (points[iregion], 1 )
4948 @views begin
50- push! (points[iregion],coord[:,cellnodes[1 ,itet]])
51- push! (points[iregion],coord[:,cellnodes[2 ,itet]])
52- push! (points[iregion],coord[:,cellnodes[3 ,itet]])
53- push! (points[iregion],coord[:,cellnodes[4 ,itet]])
54- push! (faces[iregion],(npts+ 1 , npts+ 2 , npts+ 3 ))
55- push! (faces[iregion],(npts+ 1 , npts+ 2 , npts+ 4 ))
56- push! (faces[iregion],(npts+ 2 , npts+ 3 , npts+ 4 ))
57- push! (faces[iregion],(npts+ 3 , npts+ 1 , npts+ 4 ))
49+ push! (points[iregion], coord[:, cellnodes[1 , itet]])
50+ push! (points[iregion], coord[:, cellnodes[2 , itet]])
51+ push! (points[iregion], coord[:, cellnodes[3 , itet]])
52+ push! (points[iregion], coord[:, cellnodes[4 , itet]])
53+ push! (faces[iregion], (npts + 1 , npts + 2 , npts + 3 ))
54+ push! (faces[iregion], (npts + 1 , npts + 2 , npts + 4 ))
55+ push! (faces[iregion], (npts + 2 , npts + 3 , npts + 4 ))
56+ push! (faces[iregion], (npts + 3 , npts + 1 , npts + 4 ))
5857 end
5958 end
6059 end
61- points,faces
60+ points, faces
6261end
6362
64-
65-
6663"""
6764$(SIGNATURES)
6865
@@ -72,57 +69,53 @@ Extract visible boundary faces - those not cut off by the planes
7269Return corresponding points and facets for each region for drawing as mesh (Makie,MeshCat)
7370or trisurf (pyplot)
7471"""
75- function extract_visible_bfaces3D (coord,bfacenodes,bfaceregions, nbregions, xyzcut;
76- primepoints= zeros (0 ,0 ), Tp= SVector{3 ,Float32},Tf= SVector{3 ,Int32})
72+ function extract_visible_bfaces3D (coord, bfacenodes, bfaceregions, nbregions, xyzcut;
73+ primepoints = zeros (0 , 0 ), Tp = SVector{3 , Float32}, Tf = SVector{3 , Int32})
74+ nbfaces = size (bfacenodes, 2 )
75+ cutcoord = zeros (3 )
7776
78-
79- nbfaces= size (bfacenodes,2 )
80- cutcoord= zeros (3 )
81-
82- function take (coord,simplex,xyzcut)
83- for idim= 1 : 3
84- all_gt= true
85- for inode= 1 : 3
86- c= coord[idim,simplex[inode]]- xyzcut[idim]
87- all_gt= all_gt && c> 0
77+ function take (coord, simplex, xyzcut)
78+ for idim = 1 : 3
79+ all_gt = true
80+ for inode = 1 : 3
81+ c = coord[idim, simplex[inode]] - xyzcut[idim]
82+ all_gt = all_gt && c > 0
8883 end
8984 if all_gt
9085 return false
9186 end
9287 end
9388 return true
9489 end
95-
9690
97- Tc= SVector{3 ,eltype (coord)}
98- xcoord= reinterpret (Tc,reshape (coord,(length (coord),)))
99-
100-
101- faces= [Vector {Tf} (undef,0 ) for iregion= 1 : nbregions]
102- points= [Vector {Tp} (undef,0 ) for iregion= 1 : nbregions]
103- for iregion= 1 : nbregions
104- for iprime= 1 : size (primepoints,2 )
105- @views push! (points[iregion],Tp (primepoints[:,iprime]))
91+ Tc = SVector{3 , eltype (coord)}
92+ xcoord = reinterpret (Tc, reshape (coord, (length (coord),)))
93+
94+ faces = [Vector {Tf} (undef, 0 ) for iregion = 1 : nbregions]
95+ points = [Vector {Tp} (undef, 0 ) for iregion = 1 : nbregions]
96+ for iregion = 1 : nbregions
97+ for iprime = 1 : size (primepoints, 2 )
98+ @views push! (points[iregion], Tp (primepoints[:, iprime]))
10699 end
107100 end
108101
109102 # remove some type instability here
110- function collct (points,faces)
111- trinodes= [1 ,2 , 3 ]
112- for i= 1 : nbfaces
113- iregion= bfaceregions[i]
114- trinodes[1 ]= bfacenodes[1 ,i]
115- trinodes[2 ]= bfacenodes[2 ,i]
116- trinodes[3 ]= bfacenodes[3 ,i]
117- if take (coord,trinodes,xyzcut)
118- npts= size (points[iregion],1 )
119- @views push! (points[iregion],xcoord[trinodes[1 ]])
120- @views push! (points[iregion],xcoord[trinodes[2 ]])
121- @views push! (points[iregion],xcoord[trinodes[3 ]])
122- @views push! (faces[iregion],(npts+ 1 , npts+ 2 , npts+ 3 ))
103+ function collct (points, faces)
104+ trinodes = [1 , 2 , 3 ]
105+ for i = 1 : nbfaces
106+ iregion = bfaceregions[i]
107+ trinodes[1 ] = bfacenodes[1 , i]
108+ trinodes[2 ] = bfacenodes[2 , i]
109+ trinodes[3 ] = bfacenodes[3 , i]
110+ if take (coord, trinodes, xyzcut)
111+ npts = size (points[iregion], 1 )
112+ @views push! (points[iregion], xcoord[trinodes[1 ]])
113+ @views push! (points[iregion], xcoord[trinodes[2 ]])
114+ @views push! (points[iregion], xcoord[trinodes[3 ]])
115+ @views push! (faces[iregion], (npts + 1 , npts + 2 , npts + 3 ))
123116 end
124117 end
125118 end
126- collct (points,faces)
127- points,faces
119+ collct (points, faces)
120+ points, faces
128121end
0 commit comments