SPH
MarchingCubes.cpp
Go to the documentation of this file.
1 #include "post/MarchingCubes.h"
2 #include "objects/finders/NeighbourFinder.h"
3 #include "post/Analysis.h"
4 #include "quantities/Storage.h"
5 #include "sph/kernel/Kernel.h"
6 #include "system/Factory.h"
7 #include "system/Profiler.h"
8 #include "system/Settings.h"
9 #include "thread/ThreadLocal.h"
10 
11 NAMESPACE_SPH_BEGIN
12 
15 
16 // clang-format off
17 const int MC_EDGES[256]= {
18  0x000, 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c, 0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00,
19  0x190, 0x099, 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c, 0x99c, 0x895, 0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90,
20  0x230, 0x339, 0x033, 0x13a, 0x636, 0x73f, 0x435, 0x53c, 0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30,
21  0x3a0, 0x2a9, 0x1a3, 0x0aa, 0x7a6, 0x6af, 0x5a5, 0x4ac, 0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0,
22  0x460, 0x569, 0x663, 0x76a, 0x066, 0x16f, 0x265, 0x36c, 0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, 0xa69, 0xb60,
23  0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0x0ff, 0x3f5, 0x2fc, 0xdfc, 0xcf5, 0xfff, 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0,
24  0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x055, 0x15c, 0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950,
25  0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, 0x2cf, 0x1c5, 0x0cc, 0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0,
26  0x8c0, 0x9c9, 0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc, 0x0cc, 0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9, 0x7c0,
27  0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c, 0x15c, 0x055, 0x35f, 0x256, 0x55a, 0x453, 0x759, 0x650,
28  0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc, 0x2fc, 0x3f5, 0x0ff, 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0,
29  0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f, 0xd65, 0xc6c, 0x36c, 0x265, 0x16f, 0x066, 0x76a, 0x663, 0x569, 0x460,
30  0xca0, 0xda9, 0xea3, 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac, 0x4ac, 0x5a5, 0x6af, 0x7a6, 0x0aa, 0x1a3, 0x2a9, 0x3a0,
31  0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c, 0x53c, 0x435, 0x73f, 0x636, 0x13a, 0x033, 0x339, 0x230,
32  0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c, 0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x099, 0x190,
33  0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c, 0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x000
34 };
35 
36 const int MC_TRIANGLES[256][16] =
37 {{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
38 {0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
39 {0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
40 {1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
41 {1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
42 {0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
43 {9, 2, 10, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
44 {2, 8, 3, 2, 10, 8, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1},
45 {3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
46 {0, 11, 2, 8, 11, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
47 {1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
48 {1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1},
49 {3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
50 {0, 10, 1, 0, 8, 10, 8, 11, 10, -1, -1, -1, -1, -1, -1, -1},
51 {3, 9, 0, 3, 11, 9, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1},
52 {9, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
53 {4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
54 {4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
55 {0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
56 {4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1},
57 {1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
58 {3, 4, 7, 3, 0, 4, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1},
59 {9, 2, 10, 9, 0, 2, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1},
60 {2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1, -1},
61 {8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
62 {11, 4, 7, 11, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1},
63 {9, 0, 1, 8, 4, 7, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1},
64 {4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1, -1, -1, -1},
65 {3, 10, 1, 3, 11, 10, 7, 8, 4, -1, -1, -1, -1, -1, -1, -1},
66 {1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1, -1, -1, -1},
67 {4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1, -1, -1, -1},
68 {4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1},
69 {9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
70 {9, 5, 4, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
71 {0, 5, 4, 1, 5, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
72 {8, 5, 4, 8, 3, 5, 3, 1, 5, -1, -1, -1, -1, -1, -1, -1},
73 {1, 2, 10, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
74 {3, 0, 8, 1, 2, 10, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1},
75 {5, 2, 10, 5, 4, 2, 4, 0, 2, -1, -1, -1, -1, -1, -1, -1},
76 {2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, -1, -1, -1, -1},
77 {9, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
78 {0, 11, 2, 0, 8, 11, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1},
79 {0, 5, 4, 0, 1, 5, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1},
80 {2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, -1, -1, -1, -1},
81 {10, 3, 11, 10, 1, 3, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1},
82 {4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, -1, -1, -1, -1},
83 {5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, -1, -1, -1, -1},
84 {5, 4, 8, 5, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1},
85 {9, 7, 8, 5, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
86 {9, 3, 0, 9, 5, 3, 5, 7, 3, -1, -1, -1, -1, -1, -1, -1},
87 {0, 7, 8, 0, 1, 7, 1, 5, 7, -1, -1, -1, -1, -1, -1, -1},
88 {1, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
89 {9, 7, 8, 9, 5, 7, 10, 1, 2, -1, -1, -1, -1, -1, -1, -1},
90 {10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, -1, -1, -1, -1},
91 {8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, -1, -1, -1, -1},
92 {2, 10, 5, 2, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1},
93 {7, 9, 5, 7, 8, 9, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1},
94 {9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, -1, -1, -1, -1},
95 {2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, -1, -1, -1, -1},
96 {11, 2, 1, 11, 1, 7, 7, 1, 5, -1, -1, -1, -1, -1, -1, -1},
97 {9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, -1, -1, -1, -1},
98 {5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, -1},
99 {11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, -1},
100 {11, 10, 5, 7, 11, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
101 {10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
102 {0, 8, 3, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
103 {9, 0, 1, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
104 {1, 8, 3, 1, 9, 8, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1},
105 {1, 6, 5, 2, 6, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
106 {1, 6, 5, 1, 2, 6, 3, 0, 8, -1, -1, -1, -1, -1, -1, -1},
107 {9, 6, 5, 9, 0, 6, 0, 2, 6, -1, -1, -1, -1, -1, -1, -1},
108 {5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, -1, -1, -1, -1},
109 {2, 3, 11, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
110 {11, 0, 8, 11, 2, 0, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1},
111 {0, 1, 9, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1},
112 {5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, -1, -1, -1, -1},
113 {6, 3, 11, 6, 5, 3, 5, 1, 3, -1, -1, -1, -1, -1, -1, -1},
114 {0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, -1, -1, -1, -1},
115 {3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, -1, -1, -1, -1},
116 {6, 5, 9, 6, 9, 11, 11, 9, 8, -1, -1, -1, -1, -1, -1, -1},
117 {5, 10, 6, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
118 {4, 3, 0, 4, 7, 3, 6, 5, 10, -1, -1, -1, -1, -1, -1, -1},
119 {1, 9, 0, 5, 10, 6, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1},
120 {10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, -1, -1, -1, -1},
121 {6, 1, 2, 6, 5, 1, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1},
122 {1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, -1, -1, -1, -1},
123 {8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, -1, -1, -1, -1},
124 {7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, -1},
125 {3, 11, 2, 7, 8, 4, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1},
126 {5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, -1, -1, -1, -1},
127 {0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1},
128 {9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, -1},
129 {8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, -1, -1, -1, -1},
130 {5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11, -1},
131 {0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, -1},
132 {6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, -1, -1, -1, -1},
133 {10, 4, 9, 6, 4, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
134 {4, 10, 6, 4, 9, 10, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1},
135 {10, 0, 1, 10, 6, 0, 6, 4, 0, -1, -1, -1, -1, -1, -1, -1},
136 {8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, -1, -1, -1, -1},
137 {1, 4, 9, 1, 2, 4, 2, 6, 4, -1, -1, -1, -1, -1, -1, -1},
138 {3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, -1, -1, -1, -1},
139 {0, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
140 {8, 3, 2, 8, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1},
141 {10, 4, 9, 10, 6, 4, 11, 2, 3, -1, -1, -1, -1, -1, -1, -1},
142 {0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, -1, -1, -1, -1},
143 {3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, -1, -1, -1, -1},
144 {6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1, -1},
145 {9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, -1, -1, -1, -1},
146 {8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, -1},
147 {3, 11, 6, 3, 6, 0, 0, 6, 4, -1, -1, -1, -1, -1, -1, -1},
148 {6, 4, 8, 11, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
149 {7, 10, 6, 7, 8, 10, 8, 9, 10, -1, -1, -1, -1, -1, -1, -1},
150 {0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, -1, -1, -1, -1},
151 {10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, -1, -1, -1, -1},
152 {10, 6, 7, 10, 7, 1, 1, 7, 3, -1, -1, -1, -1, -1, -1, -1},
153 {1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, -1, -1, -1, -1},
154 {2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, -1},
155 {7, 8, 0, 7, 0, 6, 6, 0, 2, -1, -1, -1, -1, -1, -1, -1},
156 {7, 3, 2, 6, 7, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
157 {2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, -1, -1, -1, -1},
158 {2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, -1},
159 {1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, -1},
160 {11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, -1, -1, -1, -1},
161 {8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, -1},
162 {0, 9, 1, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
163 {7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, -1, -1, -1, -1},
164 {7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
165 {7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
166 {3, 0, 8, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
167 {0, 1, 9, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
168 {8, 1, 9, 8, 3, 1, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1},
169 {10, 1, 2, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
170 {1, 2, 10, 3, 0, 8, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1},
171 {2, 9, 0, 2, 10, 9, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1},
172 {6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, -1, -1, -1, -1},
173 {7, 2, 3, 6, 2, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
174 {7, 0, 8, 7, 6, 0, 6, 2, 0, -1, -1, -1, -1, -1, -1, -1},
175 {2, 7, 6, 2, 3, 7, 0, 1, 9, -1, -1, -1, -1, -1, -1, -1},
176 {1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, -1, -1, -1, -1},
177 {10, 7, 6, 10, 1, 7, 1, 3, 7, -1, -1, -1, -1, -1, -1, -1},
178 {10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, -1, -1, -1, -1},
179 {0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, -1, -1, -1, -1},
180 {7, 6, 10, 7, 10, 8, 8, 10, 9, -1, -1, -1, -1, -1, -1, -1},
181 {6, 8, 4, 11, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
182 {3, 6, 11, 3, 0, 6, 0, 4, 6, -1, -1, -1, -1, -1, -1, -1},
183 {8, 6, 11, 8, 4, 6, 9, 0, 1, -1, -1, -1, -1, -1, -1, -1},
184 {9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, -1, -1, -1, -1},
185 {6, 8, 4, 6, 11, 8, 2, 10, 1, -1, -1, -1, -1, -1, -1, -1},
186 {1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, -1, -1, -1, -1},
187 {4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, -1, -1, -1, -1},
188 {10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3, -1},
189 {8, 2, 3, 8, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1},
190 {0, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
191 {1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, -1, -1, -1, -1},
192 {1, 9, 4, 1, 4, 2, 2, 4, 6, -1, -1, -1, -1, -1, -1, -1},
193 {8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, -1, -1, -1, -1},
194 {10, 1, 0, 10, 0, 6, 6, 0, 4, -1, -1, -1, -1, -1, -1, -1},
195 {4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, -1},
196 {10, 9, 4, 6, 10, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
197 {4, 9, 5, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
198 {0, 8, 3, 4, 9, 5, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1},
199 {5, 0, 1, 5, 4, 0, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1},
200 {11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, -1, -1, -1, -1},
201 {9, 5, 4, 10, 1, 2, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1},
202 {6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, -1, -1, -1, -1},
203 {7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, -1, -1, -1, -1},
204 {3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6, -1},
205 {7, 2, 3, 7, 6, 2, 5, 4, 9, -1, -1, -1, -1, -1, -1, -1},
206 {9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, -1, -1, -1, -1},
207 {3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, -1, -1, -1, -1},
208 {6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8, -1},
209 {9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7, -1, -1, -1, -1},
210 {1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4, -1},
211 {4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10, -1},
212 {7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10, -1, -1, -1, -1},
213 {6, 9, 5, 6, 11, 9, 11, 8, 9, -1, -1, -1, -1, -1, -1, -1},
214 {3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, -1, -1, -1, -1},
215 {0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, -1, -1, -1, -1},
216 {6, 11, 3, 6, 3, 5, 5, 3, 1, -1, -1, -1, -1, -1, -1, -1},
217 {1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, -1, -1, -1, -1},
218 {0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10, -1},
219 {11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5, -1},
220 {6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3, -1, -1, -1, -1},
221 {5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, -1, -1, -1, -1},
222 {9, 5, 6, 9, 6, 0, 0, 6, 2, -1, -1, -1, -1, -1, -1, -1},
223 {1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8, -1},
224 {1, 5, 6, 2, 1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
225 {1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6, -1},
226 {10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0, -1, -1, -1, -1},
227 {0, 3, 8, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
228 {10, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
229 {11, 5, 10, 7, 5, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
230 {11, 5, 10, 11, 7, 5, 8, 3, 0, -1, -1, -1, -1, -1, -1, -1},
231 {5, 11, 7, 5, 10, 11, 1, 9, 0, -1, -1, -1, -1, -1, -1, -1},
232 {10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, -1, -1, -1, -1},
233 {11, 1, 2, 11, 7, 1, 7, 5, 1, -1, -1, -1, -1, -1, -1, -1},
234 {0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, -1, -1, -1, -1},
235 {9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, -1, -1, -1, -1},
236 {7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2, -1},
237 {2, 5, 10, 2, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1},
238 {8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, -1, -1, -1, -1},
239 {9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, -1, -1, -1, -1},
240 {9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, -1},
241 {1, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
242 {0, 8, 7, 0, 7, 1, 1, 7, 5, -1, -1, -1, -1, -1, -1, -1},
243 {9, 0, 3, 9, 3, 5, 5, 3, 7, -1, -1, -1, -1, -1, -1, -1},
244 {9, 8, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
245 {5, 8, 4, 5, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1},
246 {5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, -1, -1, -1, -1},
247 {0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, -1, -1, -1, -1},
248 {10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4, -1},
249 {2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, -1, -1, -1, -1},
250 {0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11, -1},
251 {0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5, -1},
252 {9, 4, 5, 2, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
253 {2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, -1, -1, -1, -1},
254 {5, 10, 2, 5, 2, 4, 4, 2, 0, -1, -1, -1, -1, -1, -1, -1},
255 {3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9, -1},
256 {5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2, -1, -1, -1, -1},
257 {8, 4, 5, 8, 5, 3, 3, 5, 1, -1, -1, -1, -1, -1, -1, -1},
258 {0, 4, 5, 1, 0, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
259 {8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, -1, -1, -1, -1},
260 {9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
261 {4, 11, 7, 4, 9, 11, 9, 10, 11, -1, -1, -1, -1, -1, -1, -1},
262 {0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, -1, -1, -1, -1},
263 {1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, -1, -1, -1, -1},
264 {3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4, -1},
265 {4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, -1, -1, -1, -1},
266 {9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3, -1},
267 {11, 7, 4, 11, 4, 2, 2, 4, 0, -1, -1, -1, -1, -1, -1, -1},
268 {11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, -1, -1, -1, -1},
269 {2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, -1, -1, -1, -1},
270 {9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7, -1},
271 {3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10, -1},
272 {1, 10, 2, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
273 {4, 9, 1, 4, 1, 7, 7, 1, 3, -1, -1, -1, -1, -1, -1, -1},
274 {4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, -1, -1, -1, -1},
275 {4, 0, 3, 7, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
276 {4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
277 {9, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
278 {3, 0, 9, 3, 9, 11, 11, 9, 10, -1, -1, -1, -1, -1, -1, -1},
279 {0, 1, 10, 0, 10, 8, 8, 10, 11, -1, -1, -1, -1, -1, -1, -1},
280 {3, 1, 10, 11, 3, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
281 {1, 2, 11, 1, 11, 9, 9, 11, 8, -1, -1, -1, -1, -1, -1, -1},
282 {3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, -1, -1, -1, -1},
283 {0, 2, 11, 8, 0, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
284 {3, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
285 {2, 3, 8, 2, 8, 10, 10, 8, 9, -1, -1, -1, -1, -1, -1, -1},
286 {9, 10, 2, 0, 9, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
287 {2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, -1, -1, -1, -1},
288 {1, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
289 {1, 3, 8, 9, 1, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
290 {0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
291 {0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
292 {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}};
293 
294 // clang-format on
295 
296 
297 const Size IDXS1[12] = { 0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3 };
298 const Size IDXS2[12] = { 1, 2, 3, 0, 5, 6, 7, 4, 4, 5, 6, 7 };
299 
300 template <typename TFunctor>
301 bool MarchingCubes::iterateWithIndices(const Box& box, const Vector& step, TFunctor&& functor) {
302  MEASURE_SCOPE("MC - evaluating field");
303  SPH_ASSERT(box != Box::EMPTY());
304 
305  Size reportCnt = max(Size(box.size()[Z] / step[Z]), 1u);
306  Size reportStep = max(reportCnt / 100, 1u);
307  std::atomic_int counter{ 0 };
308  std::atomic_bool shouldContinue{ true };
309  auto task = [this, &step, &box, &functor, reportStep, reportCnt, &counter, &shouldContinue](
310  const Size k) {
311  const Float z = box.lower()[Z] + k * step[Z];
312  Size i = 0;
313  Size j = 0;
314  for (Float y = box.lower()[Y]; y <= box.upper()[Y]; y += step[Y], j++) {
315  i = 0;
316  for (Float x = box.lower()[X]; x <= box.upper()[X]; x += step[X], i++) {
317  functor(Indices(i, j, k), Vector(x, y, z));
318  }
319  }
320 
321  if (progressCallback && (++counter % reportStep == 0)) {
322  shouldContinue = shouldContinue && progressCallback(Float(counter) / Float(reportCnt));
323  }
324  };
325  parallelFor(scheduler, 0, Size(box.size()[Z] / step[Z]) + 1, task);
326  return shouldContinue;
327 }
328 
329 MarchingCubes::MarchingCubes(IScheduler& scheduler,
330  const Float surfaceLevel,
331  const SharedPtr<IScalarField>& field,
332  Function<bool(Float progress)> progressCallback)
333  : scheduler(scheduler)
334  , surfaceLevel(surfaceLevel)
335  , field(field)
336  , progressCallback(progressCallback) {}
337 
338 void MarchingCubes::addComponent(const Box& box, const Float gridResolution) {
339  MEASURE_SCOPE("MC addComponent");
340 
341  const Vector dr = min(Vector(gridResolution), box.size() * (1._f - EPS));
342  cached.phi.clear();
343  // multiply by (1 + EPS) to handle case where box size is divisible by dr
344  Indices cnts((1._f + EPS) * box.size() / dr);
345  SPH_ASSERT(cnts[X] >= 1 && cnts[Y] >= 1 && cnts[Z] >= 1);
346 
347  // find values of grid nodes
348  auto mapping = [&cnts](const Indices& idxs) {
349  SPH_ASSERT(idxs[X] >= 0 && idxs[X] <= cnts[X], idxs[X], cnts[X]);
350  SPH_ASSERT(idxs[Y] >= 0 && idxs[Y] <= cnts[Y], idxs[Y], cnts[Y]);
351  SPH_ASSERT(idxs[Z] >= 0 && idxs[Z] <= cnts[Z], idxs[Z], idxs[Z]);
352  return idxs[X] + (cnts[X] + 1) * idxs[Y] + (cnts[X] + 1) * (cnts[Y] + 1) * idxs[Z];
353  };
354  cached.phi.resize((cnts[X] + 1) * (cnts[Y] + 1) * (cnts[Z] + 1));
355  bool shouldContinue =
356  this->iterateWithIndices(box, dr, [this, &mapping](const Indices& idxs, const Vector& v) { //
357  cached.phi[mapping(idxs)] = (*field)(v);
358  });
359  if (!shouldContinue) {
360  return;
361  }
362 
363  // for each non-empty grid, find all intersecting triangles
364  Box boxWithoutLast(box.lower(), box.upper() - dr);
365  ThreadLocal<Array<Triangle>> tri(scheduler);
366 
367  auto intersect = [this, &dr, &mapping, &tri](const Indices& idxs, const Vector& v) {
368  Cell cell;
369  bool allBelow = true, allAbove = true;
370  Size i = 0;
371  for (Size z = 0; z <= 1; ++z) {
372  for (Size y = 0; y <= 1; ++y) {
373  for (Size x = 0; x <= 1; ++x) {
374  Size actX = (x + y) % 2; // the mapping convention is (0, 0) - (1, 0) - (1, 1) - (0, 1)
375  cell.node(i) = v + dr * Vector(actX, y, z);
376  cell.value(i) = cached.phi[mapping(idxs + Indices(actX, y, z))];
377  if (cell.value(i) > surfaceLevel) {
378  allBelow = false;
379  } else if (cell.value(i) < surfaceLevel) {
380  allAbove = false;
381  }
382  i++;
383  }
384  }
385  }
386 
387  if (!allBelow && !allAbove) {
388  this->intersectCell(cell, tri.local());
389  }
390  };
391  shouldContinue = this->iterateWithIndices(boxWithoutLast, dr, intersect);
392  if (!shouldContinue) {
393  return;
394  }
395 
396  for (Array<Triangle>& triTl : tri) {
397  triangles.pushAll(triTl);
398  }
399 }
400 
401 void MarchingCubes::intersectCell(Cell& cell, Array<Triangle>& tri) {
402  Size cubeIdx = 0;
403  for (Size i = 0; i < 8; ++i) {
404  if (cell.value(i) <= surfaceLevel) {
405  cubeIdx |= 1 << i;
406  }
407  }
408 
409  if (MC_EDGES[cubeIdx] == 0) {
410  // cube is entirely in/out of the surface
411  return;
412  }
413 
414  // find the vertices where the surface intersects the cube
415  StaticArray<Vector, 12> vertices;
416  for (Size i = 0; i < 12; ++i) {
417  if (MC_EDGES[cubeIdx] & (1 << i)) {
418  const Size k = IDXS1[i];
419  const Size l = IDXS2[i];
420  vertices[i] = this->interpolate(cell.node(k), cell.value(k), cell.node(l), cell.value(l));
421  } else {
422  vertices[i] = Vector(NAN);
423  }
424  }
425 
426  for (Size i = 0; MC_TRIANGLES[cubeIdx][i] != -1; i += 3) {
427  Triangle t;
428  t[0] = vertices[MC_TRIANGLES[cubeIdx][i + 0]];
429  t[1] = vertices[MC_TRIANGLES[cubeIdx][i + 1]];
430  t[2] = vertices[MC_TRIANGLES[cubeIdx][i + 2]];
431  if (!t.isValid()) {
432  // skip degenerated triangles
433  continue;
434  }
435  tri.push(t);
436  }
437 }
438 
439 INLINE Vector MarchingCubes::interpolate(const Vector& v1,
440  const Float p1,
441  const Vector& v2,
442  const Float p2) const {
443  if (almostEqual(p1, surfaceLevel)) {
444  return v1;
445  }
446  if (almostEqual(p2, surfaceLevel)) {
447  return v2;
448  }
449  if (almostEqual(p1, p2)) {
450  // small difference between values, just return the center to avoid instabilities
451  return 0.5_f * (v1 + v2);
452  }
453 
454  Float mu = (p1 - surfaceLevel) / (p1 - p2);
455  SPH_ASSERT(mu >= 0._f && mu <= 1._f);
456  return v1 + mu * (v2 - v1);
457 }
458 
459 namespace {
460 
461 class ColorField : public IScalarField {
462 private:
463  LutKernel<3> kernel;
464  AutoPtr<IBasicFinder> finder;
465 
466  ArrayView<const Vector> r;
467  ArrayView<const Float> m, rho;
468  ArrayView<const Size> flag;
469  ArrayView<const SymmetricTensor> G;
470  Float maxH = 0._f;
471 
472  ThreadLocal<Array<NeighbourRecord>> neighs;
473 
474 public:
486  ColorField(const Storage& storage,
487  IScheduler& scheduler,
488  const ArrayView<const Vector> r,
489  const ArrayView<const SymmetricTensor> aniso,
490  const Float maxH,
491  LutKernel<3>&& kernel,
492  AutoPtr<IBasicFinder>&& finder)
493  : kernel(std::move(kernel))
494  , finder(std::move(finder))
495  , r(r)
496  , G(aniso)
497  , maxH(maxH)
498  , neighs(scheduler) {
499  tie(m, rho) = storage.getValues<Float>(QuantityId::MASS, QuantityId::DENSITY);
500  flag = storage.getValue<Size>(QuantityId::FLAG);
501 
502  // we have to re-build the tree since we are using different positions (in general)
503  this->finder->build(scheduler, r);
504  }
505 
506  virtual Float operator()(const Vector& pos) override {
507  SPH_ASSERT(maxH > 0._f);
508  Array<NeighbourRecord>& neighsTl = neighs.local();
511  finder->findAll(pos, maxH * kernel.radius(), neighsTl);
512  Float phi = 0._f;
513 
514  // find average h of neighbours and the flag of the closest particle
515  Size closestFlag = 0;
516  Float flagDistSqr = INFTY;
517  for (NeighbourRecord& n : neighsTl) {
518  const Size j = n.index;
519  if (n.distanceSqr < flagDistSqr) {
520  closestFlag = flag[j];
521  flagDistSqr = n.distanceSqr;
522  }
523  }
524 
525  // interpolate values of neighbours
526  for (NeighbourRecord& n : neighsTl) {
527  const Size j = n.index;
528  if (flag[j] != closestFlag) {
529  continue;
530  }
531  phi += m[j] / rho[j] * G[j].determinant() * kernel.valueImpl(getSqrLength(G[j] * (pos - r[j])));
532  }
533  return phi;
534  }
535 };
536 
537 
538 class FallbackField : public IScalarField {
539 private:
540  LutKernel<3> kernel;
541  AutoPtr<IBasicFinder> finder;
542 
543  ArrayView<const Vector> r;
544  ArrayView<const SymmetricTensor> G;
545  Float maxH = 0._f;
546 
547  ThreadLocal<Array<NeighbourRecord>> neighs;
548 
549 public:
550  FallbackField(IScheduler& scheduler,
551  const ArrayView<const Vector> r,
552  const ArrayView<const SymmetricTensor> aniso,
553  const Float maxH,
554  LutKernel<3>&& kernel,
555  AutoPtr<IBasicFinder>&& finder)
556  : kernel(std::move(kernel))
557  , finder(std::move(finder))
558  , r(r)
559  , G(aniso)
560  , maxH(maxH)
561  , neighs(scheduler) {
562 
563  // we have to re-build the tree since we are using different positions (in general)
564  this->finder->build(scheduler, r);
565  }
566 
567  virtual Float operator()(const Vector& pos) override {
568  SPH_ASSERT(maxH > 0._f);
569  Array<NeighbourRecord>& neighsTl = neighs.local();
572  finder->findAll(pos, maxH * kernel.radius(), neighsTl);
573  Float phi = 0._f;
574 
575  // interpolate values of neighbours
576  for (NeighbourRecord& n : neighsTl) {
577  const Size j = n.index;
578  phi += sphereVolume(0.5_f * r[j][H]) * G[j].determinant() *
579  kernel.valueImpl(getSqrLength(G[j] * (pos - r[j])));
580  }
581  return phi;
582  }
583 };
584 
585 } // namespace
586 
587 INLINE Float weight(const Vector& r1, const Vector& r2) {
588  const Float lengthSqr = getSqrLength(r1 - r2);
589  // Eq. (11)
590  if (lengthSqr < sqr(2._f * r1[H])) {
591  return 1._f - pow<3>(sqrt(lengthSqr) / (2._f * r1[H]));
592  } else {
593  return 0._f;
594  }
595 }
596 
597 Array<Triangle> getSurfaceMesh(IScheduler& scheduler, const Storage& storage, const McConfig& config) {
598  MEASURE_SCOPE("getSurfaceMesh");
599 
600  // (according to http://www.cc.gatech.edu/~turk/my_papers/sph_surfaces.pdf)
601 
602  ArrayView<const Vector> r = storage.getValue<Vector>(QuantityId::POSITION);
603  RunSettings settings;
604  LutKernel<3> kernel = Factory::getKernel<3>(settings);
605  AutoPtr<IBasicFinder> finder = Factory::getFinder(settings);
606 
607  finder->build(scheduler, r);
608 
609  Array<Vector> r_bar(r.size());
610  Array<SymmetricTensor> G(r.size()); // anisotropy matrix
611 
612  ThreadLocal<Array<NeighbourRecord>> neighsData(scheduler);
613  parallelFor(scheduler, neighsData, 0, r.size(), [&](const Size i, Array<NeighbourRecord>& neighs) {
615  r_bar[i] = r[i];
616  r_bar[i][H] = r[i][H] * config.smoothingMult;
617 
618  if (config.useAnisotropicKernels) {
619  Vector r_center = Vector(0._f);
620  finder->findAll(r_bar[i], 2 * r_bar[i][H], neighs);
621  for (const NeighbourRecord& n : neighs) {
622  r_center += r_bar[n.index];
623  }
624  r_center /= neighs.size();
625  SymmetricTensor C = SymmetricTensor::null();
626  for (const NeighbourRecord& n : neighs) {
627  C += symmetricOuter(r[n.index] - r_center, r[n.index] - r_center);
628  }
629 
630  Svd svd = singularValueDecomposition(C);
631  const Float maxSigma = maxElement(svd.S);
632  for (Size i = 0; i < 3; ++i) {
633  svd.S[i] = 1._f / std::max(svd.S[i], 0.125_f * maxSigma);
634  }
635 
636  AffineMatrix sigma = convert<AffineMatrix>(SymmetricTensor(svd.S, Vector(0._f)));
637  G[i] = convert<SymmetricTensor>(svd.V * sigma * svd.U.transpose());
638  } else {
639  G[i] = SymmetricTensor(Vector(1._f / r[i][H]), Vector(0._f));
640  }
641  });
642  // 5. find bounding box and maximum h (we need to search neighbours of arbitrary point in space)
643 
644  Float maxH = 0._f;
645  for (Size i = 0; i < r_bar.size(); ++i) {
646  maxH = max(maxH, r_bar[i][H]);
647  }
648  SharedPtr<IScalarField> field;
649 
650  if (storage.has(QuantityId::MASS) && storage.has(QuantityId::DENSITY) && storage.has(QuantityId::FLAG)) {
651  field =
652  makeShared<ColorField>(storage, scheduler, r_bar, G, maxH, std::move(kernel), std::move(finder));
653  } else {
654  field = makeShared<FallbackField>(scheduler, r_bar, G, maxH, std::move(kernel), std::move(finder));
655  }
656 
657  MarchingCubes mc(scheduler, config.surfaceLevel, field, config.progressCallback);
658 
659  Array<Size> components;
660  const Size numComponents = Post::findComponents(storage, 2._f, Post::ComponentFlag::OVERLAP, components);
661 
662  // 6. find the surface using marching cubes for each component
663  Array<Box> boxes(numComponents);
664  Array<Size> counts(numComponents);
665  counts.fill(0);
666  for (Size j = 0; j < components.size(); ++j) {
667  const Vector padding(max(2._f * r_bar[j][H], 2._f * config.gridResolution));
668  boxes[components[j]].extend(r_bar[j] + padding);
669  boxes[components[j]].extend(r_bar[j] - padding);
670  counts[components[j]]++;
671  }
672  for (Size i = 0; i < numComponents; ++i) {
673  if (counts[i] > 10) {
674  mc.addComponent(boxes[i], config.gridResolution);
675  }
676  }
677 
678  return std::move(mc.getTriangles());
679 }
680 
681 
682 NAMESPACE_SPH_END
almostEqual
INLINE bool almostEqual(const AffineMatrix &m1, const AffineMatrix &m2, const Float eps=EPS)
Definition: AffineMatrix.h:278
Analysis.h
Various function for interpretation of the results of a simulation.
SPH_ASSERT
#define SPH_ASSERT(x,...)
Definition: Assert.h:94
NAMESPACE_SPH_BEGIN
NAMESPACE_SPH_BEGIN
Definition: BarnesHut.cpp:13
Size
uint32_t Size
Integral type used to index arrays (by default).
Definition: Globals.h:16
Float
double Float
Precision used withing the code. Use Float instead of float or double where precision is important.
Definition: Globals.h:13
Kernel.h
SPH kernels.
MC_TRIANGLES
const int MC_TRIANGLES[256][16]
Definition: MarchingCubes.cpp:36
getSurfaceMesh
Array< Triangle > getSurfaceMesh(IScheduler &scheduler, const Storage &storage, const McConfig &config)
Returns the triangle mesh of the body surface (or surfaces of bodies).
Definition: MarchingCubes.cpp:597
weight
INLINE Float weight(const Vector &r1, const Vector &r2)
Definition: MarchingCubes.cpp:587
MC_EDGES
NAMESPACE_SPH_BEGIN const int MC_EDGES[256]
Definition: MarchingCubes.cpp:17
IDXS2
const Size IDXS2[12]
Definition: MarchingCubes.cpp:298
IDXS1
const Size IDXS1[12]
Definition: MarchingCubes.cpp:297
max
constexpr INLINE T max(const T &f1, const T &f2)
Definition: MathBasic.h:20
min
NAMESPACE_SPH_BEGIN constexpr INLINE T min(const T &f1, const T &f2)
Minimum & Maximum value.
Definition: MathBasic.h:15
EPS
constexpr Float EPS
Definition: MathUtils.h:30
sqr
constexpr INLINE T sqr(const T &f) noexcept
Return a squared value.
Definition: MathUtils.h:67
sqrt
INLINE T sqrt(const T f)
Return a squared root of a value.
Definition: MathUtils.h:78
pow< 3 >
constexpr INLINE Float pow< 3 >(const Float v)
Definition: MathUtils.h:132
sphereVolume
constexpr INLINE Float sphereVolume(const Float radius)
Computes a volume of a sphere given its radius.
Definition: MathUtils.h:375
INFTY
constexpr Float INFTY
Definition: MathUtils.h:38
INLINE
#define INLINE
Macros for conditional compilation based on selected compiler.
Definition: Object.h:31
NAMESPACE_SPH_END
#define NAMESPACE_SPH_END
Definition: Object.h:12
Profiler.h
Tool to measure time spent in functions and profile the code.
MEASURE_SCOPE
#define MEASURE_SCOPE(name)
Definition: Profiler.h:70
QuantityId::FLAG
@ FLAG
ID of original body, used to implement discontinuities between bodies in SPH.
QuantityId::POSITION
@ POSITION
Positions (velocities, accelerations) of particles, always a vector quantity,.
QuantityId::DENSITY
@ DENSITY
Density, always a scalar quantity.
QuantityId::MASS
@ MASS
Paricles masses, always a scalar quantity.
parallelFor
INLINE void parallelFor(IScheduler &scheduler, const Size from, const Size to, TFunctor &&functor)
Executes a functor concurrently from all available threads.
Definition: Scheduler.h:98
tie
StaticArray< T0 &, sizeof...(TArgs)+1 > tie(T0 &t0, TArgs &... rest)
Creates a static array from a list of l-value references.
Definition: StaticArray.h:281
Storage.h
Container for storing particle quantities and materials.
ThreadLocal.h
Template for thread-local storage.
getSqrLength
INLINE Float getSqrLength(const Vector &v)
Definition: Vector.h:574
Vector
BasicVector< Float > Vector
Definition: Vector.h:539
H
@ H
Definition: Vector.h:25
Y
@ Y
Definition: Vector.h:23
X
@ X
Definition: Vector.h:22
Z
@ Z
Definition: Vector.h:24
ArrayView
Object providing safe access to continuous memory of data.
Definition: ArrayView.h:17
ArrayView::size
INLINE TCounter size() const
Definition: ArrayView.h:101
Array::push
INLINE void push(U &&u)
Adds new element to the end of the array, resizing the array if necessary.
Definition: Array.h:306
Array::size
INLINE TCounter size() const noexcept
Definition: Array.h:193
Array::pushAll
void pushAll(const TIter first, const TIter last)
Definition: Array.h:312
Box
Helper object defining three-dimensional interval (box).
Definition: Box.h:17
Box::lower
INLINE const Vector & lower() const
Returns lower bounds of the box.
Definition: Box.h:82
Box::upper
INLINE const Vector & upper() const
Returns upper bounds of the box.
Definition: Box.h:94
Box::size
INLINE Vector size() const
Returns box dimensions.
Definition: Box.h:106
Box::EMPTY
static Box EMPTY()
Syntactic sugar, returns a default-constructed (empty) box.
Definition: Box.h:41
Cell
Single cell used in mesh generation.
Definition: MarchingCubes.h:16
Cell::value
INLINE Float & value(const Size idx)
Definition: MarchingCubes.h:22
Cell::node
INLINE Vector & node(const Size idx)
Definition: MarchingCubes.h:26
IBasicFinder::build
void build(IScheduler &scheduler, ArrayView< const Vector > points)
Constructs the struct with an array of vectors.
Definition: NeighbourFinder.cpp:6
IBasicFinder::findAll
virtual Size findAll(const Size index, const Float radius, Array< NeighbourRecord > &neighbours) const =0
Finds all neighbours within given radius from the point given by index.
IScalarField
Inferface for a generic scalar field, returning a float for given position.:w.
Definition: MarchingCubes.h:32
IScalarField::operator()
virtual Float operator()(const Vector &pos)=0
Returns the value of the scalar field at given position.
IScheduler
Interface that allows unified implementation of sequential and parallelized versions of algorithms.
Definition: Scheduler.h:27
Indices
Helper object for storing three (possibly four) int or bool values.
Definition: Indices.h:16
LutKernel::valueImpl
INLINE Float valueImpl(const Float qSqr) const noexcept
Definition: Kernel.h:111
LutKernel::radius
INLINE Float radius() const noexcept
Definition: Kernel.h:107
MarchingCubes
Marching cubes algorithm for generation of mesh from iso-surface of given scalar field.
Definition: MarchingCubes.h:39
MarchingCubes::addComponent
void addComponent(const Box &box, const Float gridResolution)
Adds a triangle mesh representing the boundary of particles.
Definition: MarchingCubes.cpp:338
MarchingCubes::MarchingCubes
MarchingCubes(IScheduler &scheduler, const Float surfaceLevel, const SharedPtr< IScalarField > &field, Function< bool(Float progress)> progressCallback=nullptr)
Constructs the object using given scalar field.
Definition: MarchingCubes.cpp:329
StaticArray
Array with fixed number of allocated elements.
Definition: StaticArray.h:19
Storage
Container storing all quantities used within the simulations.
Definition: Storage.h:230
Storage::getValues
auto getValues(const QuantityId first, const QuantityId second, const TArgs... others)
Retrieves an array of quantities from the key.
Definition: Storage.h:359
Storage::has
bool has(const QuantityId key) const
Checks if the storage contains quantity with given key.
Definition: Storage.cpp:130
Storage::getValue
Array< TValue > & getValue(const QuantityId key)
Retrieves a quantity values from the storage, given its key and value type.
Definition: Storage.cpp:191
ThreadLocal
Template for storing a copy of a value for every thread in given scheduler.
Definition: ThreadLocal.h:36
ThreadLocal::local
INLINE Type & local()
Return a value for current thread.
Definition: ThreadLocal.h:88
Triangle
Triangle.h.
Definition: Triangle.h:14
Triangle::isValid
INLINE bool isValid() const
Definition: Triangle.h:70
Factory.h
Creating code components based on values from settings.
Settings.h
Generic storage and input/output routines of settings.
Factory::getFinder
AutoPtr< ISymmetricFinder > getFinder(const RunSettings &settings)
Definition: Factory.cpp:156
Post::ComponentFlag::OVERLAP
@ OVERLAP
Specifies that overlapping particles belong into the same component.
Post::findComponents
Size findComponents(const Storage &storage, const Float particleRadius, const Flags< ComponentFlag > flags, Array< Size > &indices)
Finds and marks connected components (a.k.a. separated bodies) in the array of vertices.
Definition: Analysis.cpp:86
std
Overload of std::swap for Sph::Array.
Definition: Array.h:578
NeighbourRecord
Holds information about a neighbour particles.
Definition: NeighbourFinder.h:18
Generated by doxygen 1.9.1