1
+ function assign_coordinates_at_point_in_plane! (
2
+ ixcoord,
3
+ ixvalues,
4
+ coordinates,
5
+ function_values,
6
+ node_index,
7
+ intersection_count
8
+ )
9
+
10
+ @views ixcoord[:, intersection_count] .= coordinates[:, node_index]
11
+ ixvalues[intersection_count] = function_values[node_index]
12
+
13
+ return nothing
14
+
15
+ end
16
+
17
+ function assign_coordinates_at_intersection! (
18
+ ixcoord,
19
+ ixvalues,
20
+ coordinates,
21
+ function_values,
22
+ planeq_values_start,
23
+ planeq_values_end,
24
+ node_index_start,
25
+ node_index_end,
26
+ intersection_count
27
+ )
28
+
29
+ t = planeq_values_start / (planeq_values_start - planeq_values_end)
30
+ @views @. ixcoord[:, intersection_count] = coordinates[:, node_index_start] + t * (coordinates[:, node_index_end] - coordinates[:, node_index_start])
31
+ ixvalues[intersection_count] = function_values[node_index_start] + t * (function_values[node_index_end] - function_values[node_index_start])
32
+
33
+ return nothing
34
+ end
35
+
36
+
1
37
"""
2
38
$(SIGNATURES)
3
39
Calculate intersections between tetrahedron with given piecewise linear
10
46
and the plane.
11
47
12
48
Input:
13
- - pointlist : 3xN array of grid point coordinates
14
- - node_indices: 4 element array of node indices (pointing into pointlist and function_values)
49
+ - coordinates : 3xN array of grid point coordinates
50
+ - node_indices: 4 element array of node indices (pointing into coordinates and function_values)
15
51
- planeq_values: 4 element array of plane equation evaluated at the node coordinates
16
52
- function_values: N element array of function values
17
53
20
56
- ixvalues: 4 element array of function values at plane - tetdedge intersections
21
57
22
58
Returns:
23
- - nxs ,ixcoord,ixvalues
59
+ - amount_intersections ,ixcoord,ixvalues
24
60
25
61
This method can be used both for the evaluation of plane sections and for
26
62
the evaluation of function isosurfaces.
27
63
"""
28
- function tet_x_plane ! (
64
+ function calculate_plane_tetrahedron_intersection ! (
29
65
ixcoord,
30
66
ixvalues,
31
- pointlist ,
67
+ coordinates ,
32
68
node_indices,
33
69
planeq_values,
34
70
function_values;
35
71
tol = 0.0
36
72
)
37
73
38
- # If all nodes lie on one side of the plane, no intersection
39
- @fastmath if (
40
- mapreduce (a -> a < - tol, * , planeq_values) ||
41
- mapreduce (a -> a > tol, * , planeq_values)
42
- )
43
- return 0
44
- end
45
- # Interpolate coordinates and function_values according to
46
- # evaluation of the plane equation
47
- nxs = 0
48
- @inbounds @simd for n1 in 1 : 3
74
+
75
+ amount_intersections = 0
76
+
77
+ @inbounds for n1 in 1 : 4
49
78
N1 = node_indices[n1]
50
- @inbounds @fastmath @simd for n2 in (n1 + 1 ): 4
51
- N2 = node_indices[n2]
52
- if planeq_values[n1] != planeq_values[n2] &&
53
- planeq_values[n1] * planeq_values[n2] < tol
54
- nxs += 1
55
- t = planeq_values[n1] / (planeq_values[n1] - planeq_values[n2])
56
- ixcoord[1 , nxs] = pointlist[1 , N1] + t * (pointlist[1 , N2] - pointlist[1 , N1])
57
- ixcoord[2 , nxs] = pointlist[2 , N1] + t * (pointlist[2 , N2] - pointlist[2 , N1])
58
- ixcoord[3 , nxs] = pointlist[3 , N1] + t * (pointlist[3 , N2] - pointlist[3 , N1])
59
- ixvalues[nxs] = function_values[N1] + t * (function_values[N2] - function_values[N1])
79
+ if abs (planeq_values[n1]) < tol
80
+ amount_intersections += 1
81
+ assign_coordinates_at_point_in_plane! (ixcoord, ixvalues, coordinates, function_values, N1, amount_intersections)
82
+ else
83
+ for n2 in (n1 + 1 ): 4
84
+ N2 = node_indices[n2]
85
+ if (abs (planeq_values[n2]) < tol) # We do not allow the 2nd node to be in the plane
86
+ continue
87
+ end
88
+ if planeq_values[n1] * planeq_values[n2] < tol^ 2
89
+ amount_intersections += 1
90
+ assign_coordinates_at_intersection! (ixcoord, ixvalues, coordinates, function_values, planeq_values[n1], planeq_values[n2], N1, N2, amount_intersections)
91
+ end
60
92
end
61
93
end
62
94
end
63
- return nxs
95
+
96
+ if amount_intersections > 4
97
+ @warn " computed $(amount_intersections) intersection points of a tetrahedron and a plane. Expected at most 4."
98
+ end
99
+ return amount_intersections
64
100
end
65
101
66
102
"""
@@ -190,7 +226,7 @@ function marching_tetrahedra(
190
226
191
227
function pushtris (ns, ixcoord, ixvalues)
192
228
# number of intersection points can be 3 or 4
193
- return if ns >= 3
229
+ if ns >= 3
194
230
last_i = length (all_ixvalues)
195
231
for is in 1 : ns
196
232
@views push! (all_ixcoord, ixcoord[:, is])
@@ -201,6 +237,7 @@ function marching_tetrahedra(
201
237
push! (all_ixfaces, (last_i + 3 , last_i + 2 , last_i + 4 ))
202
238
end
203
239
end
240
+ return nothing
204
241
end
205
242
206
243
for igrid in 1 : length (allcoords)
@@ -221,7 +258,7 @@ function marching_tetrahedra(
221
258
planeq[2 ] = all_planeq[node_indices[2 ]]
222
259
planeq[3 ] = all_planeq[node_indices[3 ]]
223
260
planeq[4 ] = all_planeq[node_indices[4 ]]
224
- nxs = tet_x_plane ! (
261
+ nxs = calculate_plane_tetrahedron_intersection ! (
225
262
ixcoord,
226
263
ixvalues,
227
264
coord,
0 commit comments