PeriDyno 1.0.0
Loading...
Searching...
No Matches
SimplexSet.cu
Go to the documentation of this file.
1#include "SimplexSet.h"
2
3#include <thrust/sort.h>
4
5namespace dyno
6{
7 template<typename TDataType>
8 SimplexSet<TDataType>::SimplexSet()
9 : PointSet<TDataType>()
10 {
11 }
12
13 template<typename TDataType>
14 SimplexSet<TDataType>::~SimplexSet()
15 {
16 mEdgeIndex.clear();
17 mTriangleIndex.clear();
18 mTetrahedronIndex.clear();
19 }
20
21 template<typename TDataType>
22 void SimplexSet<TDataType>::copyFrom(SimplexSet<TDataType>& simplex)
23 {
24 PointSet<TDataType>::copyFrom(simplex);
25
26 mEdgeIndex.assign(simplex.mEdgeIndex);
27 mTriangleIndex.assign(simplex.mTriangleIndex);
28 mTetrahedronIndex.assign(simplex.mTetrahedronIndex);
29 }
30
31 template<typename TDataType>
32 bool SimplexSet<TDataType>::isEmpty()
33 {
34 bool empty = true;
35 empty |= mEdgeIndex.size() == 0;
36 empty |= mTriangleIndex.size() == 0;
37 empty |= mTetrahedronIndex.size() == 0;
38
39 return empty;
40 }
41
42 template<typename TDataType>
43 void SimplexSet<TDataType>::updateTopology()
44 {
45
46 }
47
48 template<typename Edge>
49 __global__ void Simplex_SetEdgeIndexIds(
50 DArray<uint> ids,
51 DArray<Edge> edges)
52 {
53 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
54 if (tId >= edges.size()) return;
55
56 Edge edge = edges[tId];
57
58 ids[2 * tId] = edge[0];
59 ids[2 * tId + 1] = edge[1];
60 }
61
62 __global__ void Simplex_CountVertexNumber(
63 DArray<int> counter,
64 DArray<uint> vIds)
65 {
66 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
67 if (tId >= vIds.size()) return;
68
69 if (tId == 0 || vIds[tId] != vIds[tId - 1])
70 counter[tId] = 1;
71 else
72 counter[tId] = 0;
73 }
74
75 template<typename Coord>
76 __global__ void Simplex_SetupEdgeVertices(
77 DArray<Coord> vertices,
78 DArray<Coord> originVertices,
79 DArray<uint> vertexIdMapper,
80 DArray<uint> vIds,
81 DArray<int> radix)
82 {
83 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
84 if (tId >= radix.size()) return;
85
86 int shift = radix[tId];
87 if (tId == 0 || radix[tId] != radix[tId - 1]) {
88 uint vId = vIds[tId];
89 vertices[shift] = originVertices[vId];
90 vertexIdMapper[vId] = shift;
91 }
92 }
93
94 template<typename Edge>
95 __global__ void Simplex_UpdateEdgeIndex(
96 DArray<Edge> newEdgeIndices,
97 DArray<Edge> oldEdgeIndices,
98 DArray<uint> vertexIdMapper)
99 {
100 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
101 if (tId >= oldEdgeIndices.size()) return;
102
103 Edge edge = oldEdgeIndices[tId];
104
105 uint v0 = vertexIdMapper[edge[0]];
106 uint v1 = vertexIdMapper[edge[1]];
107
108 newEdgeIndices[tId] = Edge(v0, v1);
109 }
110
111 template<typename TDataType>
112 void SimplexSet<TDataType>::extractSimplex1D(EdgeSet<TDataType>& es)
113 {
114 es.clear();
115
116 uint eNum = mEdgeIndex.size();
117
118 DArray<uint> vIds(eNum * 2);
119
120 cuExecute(eNum,
121 Simplex_SetEdgeIndexIds,
122 vIds,
123 mEdgeIndex);
124
125 thrust::sort(thrust::device, vIds.begin(), vIds.begin() + vIds.size());
126
127 DArray<int> radix(eNum * 2);
128
129 cuExecute(vIds.size(),
130 Simplex_CountVertexNumber,
131 radix,
132 vIds);
133
134 int vNum = thrust::reduce(thrust::device, radix.begin(), radix.begin() + radix.size());
135 thrust::exclusive_scan(thrust::device, radix.begin(), radix.begin() + radix.size(), radix.begin());
136
137 DArray<Coord> vertices(vNum);
138 DArray<Edge> edgeIndices(mEdgeIndex.size());
139 DArray<uint> vertexIdMapper(PointSet<TDataType>::mCoords.size());
140
141 cuExecute(radix.size(),
142 Simplex_SetupEdgeVertices,
143 vertices,
144 PointSet<TDataType>::mCoords,
145 vertexIdMapper,
146 vIds,
147 radix);
148
149 cuExecute(eNum,
150 Simplex_UpdateEdgeIndex,
151 edgeIndices,
152 mEdgeIndex,
153 vertexIdMapper);
154
155 es.setPoints(vertices);
156 es.setEdges(edgeIndices);
157 es.update();
158
159 vIds.clear();
160 radix.clear();
161 vertices.clear();
162 vertexIdMapper.clear();
163 edgeIndices.clear();
164 }
165
166 template<typename Triangle>
167 __global__ void Simplex_SetTriangleIndexIds(
168 DArray<uint> ids,
169 DArray<Triangle> triangles)
170 {
171 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
172 if (tId >= triangles.size()) return;
173
174 Triangle tri = triangles[tId];
175
176 ids[3 * tId] = tri[0];
177 ids[3 * tId + 1] = tri[1];
178 ids[3 * tId + 2] = tri[2];
179 }
180
181 template<typename Triangle>
182 __global__ void Simplex_UpdateTriangleIndex(
183 DArray<Triangle> newTriangleIndices,
184 DArray<Triangle> oldTriangleIndices,
185 DArray<uint> vertexIdMapper)
186 {
187 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
188 if (tId >= oldTriangleIndices.size()) return;
189
190 Triangle tri = oldTriangleIndices[tId];
191
192 uint v0 = vertexIdMapper[tri[0]];
193 uint v1 = vertexIdMapper[tri[1]];
194 uint v2 = vertexIdMapper[tri[2]];
195
196 newTriangleIndices[tId] = Triangle(v0, v1, v2);
197 }
198
199 template<typename TDataType>
200 void SimplexSet<TDataType>::extractSimplex2D(TriangleSet<TDataType>& ts)
201 {
202 ts.clear();
203
204 uint tNum = mTriangleIndex.size();
205
206 DArray<uint> vIds(tNum * 3);
207
208 cuExecute(tNum,
209 Simplex_SetTriangleIndexIds,
210 vIds,
211 mTriangleIndex);
212
213 thrust::sort(thrust::device, vIds.begin(), vIds.begin() + vIds.size());
214
215 DArray<int> radix(tNum * 3);
216
217 cuExecute(vIds.size(),
218 Simplex_CountVertexNumber,
219 radix,
220 vIds);
221
222 int vNum = thrust::reduce(thrust::device, radix.begin(), radix.begin() + radix.size());
223 thrust::exclusive_scan(thrust::device, radix.begin(), radix.begin() + radix.size(), radix.begin());
224
225 DArray<Coord> vertices(vNum);
226 DArray<Triangle> triangleIndex(mTriangleIndex.size());
227 DArray<uint> vertexIdMapper(PointSet<TDataType>::mCoords.size());
228
229 cuExecute(radix.size(),
230 Simplex_SetupEdgeVertices,
231 vertices,
232 PointSet<TDataType>::mCoords,
233 vertexIdMapper,
234 vIds,
235 radix);
236
237 cuExecute(tNum,
238 Simplex_UpdateTriangleIndex,
239 triangleIndex,
240 mTriangleIndex,
241 vertexIdMapper);
242
243 ts.setPoints(vertices);
244 ts.setTriangles(triangleIndex);
245 ts.update();
246
247 vIds.clear();
248 radix.clear();
249 vertices.clear();
250 vertexIdMapper.clear();
251 triangleIndex.clear();
252 }
253
254 template<typename Tetrahedron>
255 __global__ void Simplex_SetTetrahedronIndexIds(
256 DArray<uint> ids,
257 DArray<Tetrahedron> tetrahedrons)
258 {
259 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
260 if (tId >= tetrahedrons.size()) return;
261
262 Tetrahedron tet = tetrahedrons[tId];
263
264 ids[4 * tId] = tet[0];
265 ids[4 * tId + 1] = tet[1];
266 ids[4 * tId + 2] = tet[2];
267 ids[4 * tId + 3] = tet[3];
268 }
269
270 template<typename Tetrahedron>
271 __global__ void Simplex_UpdateTetrahedronIndex(
272 DArray<Tetrahedron> newTetIndices,
273 DArray<Tetrahedron> oldTetIndices,
274 DArray<uint> vertexIdMapper)
275 {
276 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
277 if (tId >= oldTetIndices.size()) return;
278
279 Tetrahedron tet = oldTetIndices[tId];
280
281 uint v0 = vertexIdMapper[tet[0]];
282 uint v1 = vertexIdMapper[tet[1]];
283 uint v2 = vertexIdMapper[tet[2]];
284 uint v3 = vertexIdMapper[tet[3]];
285
286 newTetIndices[tId] = Tetrahedron(v0, v1, v2, v3);
287 }
288
289 template<typename TDataType>
290 void SimplexSet<TDataType>::extractSimplex3D(TetrahedronSet<TDataType>& ts)
291 {
292 ts.clear();
293
294 uint tNum = mTetrahedronIndex.size();
295
296 DArray<uint> vIds(tNum * 4);
297
298 cuExecute(tNum,
299 Simplex_SetTetrahedronIndexIds,
300 vIds,
301 mTetrahedronIndex);
302
303 thrust::sort(thrust::device, vIds.begin(), vIds.begin() + vIds.size());
304
305 DArray<int> radix(tNum * 4);
306
307 cuExecute(vIds.size(),
308 Simplex_CountVertexNumber,
309 radix,
310 vIds);
311
312 int vNum = thrust::reduce(thrust::device, radix.begin(), radix.begin() + radix.size());
313 thrust::exclusive_scan(thrust::device, radix.begin(), radix.begin() + radix.size(), radix.begin());
314
315 DArray<Coord> vertices(vNum);
316 DArray<Tetrahedron> tetrahedronIndex(mTetrahedronIndex.size());
317 DArray<uint> vertexIdMapper(PointSet<TDataType>::mCoords.size());
318
319 cuExecute(radix.size(),
320 Simplex_SetupEdgeVertices,
321 vertices,
322 PointSet<TDataType>::mCoords,
323 vertexIdMapper,
324 vIds,
325 radix);
326
327 cuExecute(tNum,
328 Simplex_UpdateTetrahedronIndex,
329 tetrahedronIndex,
330 mTetrahedronIndex,
331 vertexIdMapper);
332
333 ts.setPoints(vertices);
334 ts.setTetrahedrons(tetrahedronIndex);
335 ts.update();
336
337 vIds.clear();
338 radix.clear();
339 vertices.clear();
340 vertexIdMapper.clear();
341 tetrahedronIndex.clear();
342 }
343
344 template<typename TDataType>
345 void SimplexSet<TDataType>::extractPointSet(PointSet<TDataType>& ps)
346 {
347 ps.clear();
348
349 ps.setPoints(PointSet<TDataType>::mCoords);
350 ps.update();
351 }
352
353 template<typename Edge, typename Triangle, typename Tetrahedron>
354 __global__ void Simplex_SetupEdgeKeyAndValue(
355 DArray<EKey> keys,
356 DArray<Edge> values,
357 DArray<Edge> edgeIndices,
358 DArray<Triangle> triIndices,
359 DArray<Tetrahedron> tetIndices)
360 {
361 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
362
363 uint total_num = edgeIndices.size() + triIndices.size() + tetIndices.size();
364
365 if (tId >= total_num) return;
366
367 if (tId < edgeIndices.size())
368 {
369 Edge eIdx = edgeIndices[tId];
370
371 uint v0 = eIdx[0];
372 uint v1 = eIdx[1];
373
374 EKey key(v0, v1);
375
376 keys[tId] = key;
377 values[tId] = eIdx;
378 }
379 else if (tId < edgeIndices.size() + triIndices.size())
380 {
381 uint i = tId - edgeIndices.size();
382
383 Triangle tIdx = triIndices[i];
384
385 uint v0 = tIdx[0];
386 uint v1 = tIdx[1];
387 uint v2 = tIdx[2];
388
389 uint offset = edgeIndices.size() + 3 * i;
390
391 keys[offset] = EKey(v0, v1);
392 values[offset] = Edge(v0, v1);
393
394 keys[offset + 1] = EKey(v1, v2);
395 values[offset + 1] = Edge(v1, v2);
396
397 keys[offset + 2] = EKey(v2, v0);
398 values[offset + 2] = Edge(v2, v0);
399 }
400 else
401 {
402 uint i = tId - edgeIndices.size() - triIndices.size();
403
404 Tetrahedron tIdx = tetIndices[i];
405
406 uint v0 = tIdx[0];
407 uint v1 = tIdx[1];
408 uint v2 = tIdx[2];
409 uint v3 = tIdx[3];
410
411 uint offset = edgeIndices.size() + 3 * triIndices.size() + 6 * i;
412
413 keys[offset] = EKey(v0, v1);
414 values[offset] = Edge(v0, v1);
415
416 keys[offset + 1] = EKey(v1, v2);
417 values[offset + 1] = Edge(v1, v2);
418
419 keys[offset + 2] = EKey(v2, v0);
420 values[offset + 2] = Edge(v2, v0);
421
422 keys[offset + 3] = EKey(v0, v3);
423 values[offset + 3] = Edge(v0, v3);
424
425 keys[offset + 4] = EKey(v1, v3);
426 values[offset + 4] = Edge(v1, v3);
427
428 keys[offset + 5] = EKey(v2, v3);
429 values[offset + 5] = Edge(v2, v3);
430 }
431 }
432
433 __global__ void Simplex_CountEdgeNumber(
434 DArray<int> counter,
435 DArray<EKey> vIds)
436 {
437 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
438 if (tId >= vIds.size()) return;
439
440 if (tId == 0 || vIds[tId] != vIds[tId - 1])
441 counter[tId] = 1;
442 else
443 counter[tId] = 0;
444 }
445
446 template<typename Edge>
447 __global__ void Simplex_SetupEdges(
448 DArray<Edge> newEdges,
449 DArray<Edge> oldEdges,
450 DArray<EKey> keys,
451 DArray<int> radix)
452 {
453 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
454 if (tId >= radix.size()) return;
455
456 int shift = radix[tId];
457 if (tId == 0 || radix[tId] != radix[tId - 1]) {
458 newEdges[shift] = oldEdges[tId];
459 }
460 }
461
462 template<typename TDataType>
463 void SimplexSet<TDataType>::extractEdgeSet(EdgeSet<TDataType>& es)
464 {
465 es.clear();
466
467 uint priNum = mEdgeIndex.size() + mTriangleIndex.size() + mTetrahedronIndex.size();
468 uint eNum = mEdgeIndex.size() + 3 * mTriangleIndex.size() + 6 * mTetrahedronIndex.size();
469
470 DArray<EKey> keys(eNum);
471 DArray<Edge> edges(eNum);
472
473 cuExecute(priNum,
474 Simplex_SetupEdgeKeyAndValue,
475 keys,
476 edges,
477 mEdgeIndex,
478 mTriangleIndex,
479 mTetrahedronIndex);
480
481 thrust::sort_by_key(thrust::device, keys.begin(), keys.begin() + keys.size(), edges.begin());
482
483 DArray<int> radix(eNum);
484
485 cuExecute(eNum,
486 Simplex_CountEdgeNumber,
487 radix,
488 keys);
489
490 int N = thrust::reduce(thrust::device, radix.begin(), radix.begin() + radix.size());
491 thrust::exclusive_scan(thrust::device, radix.begin(), radix.begin() + radix.size(), radix.begin());
492
493 DArray<Edge> distinctEdges(N);
494
495 cuExecute(eNum,
496 Simplex_SetupEdges,
497 distinctEdges,
498 edges,
499 keys,
500 radix);
501
502 es.setPoints(PointSet<TDataType>::mCoords);
503 es.setEdges(distinctEdges);
504 es.update();
505
506 keys.clear();
507 edges.clear();
508 radix.clear();
509 distinctEdges.clear();
510 }
511
512 template<typename TDataType>
513 void SimplexSet<TDataType>::extractTriangleSet(TriangleSet<TDataType>& ts)
514 {
515
516 }
517
518 DEFINE_CLASS(SimplexSet);
519}
Generated by doxygen 1.13.2