tcontactsearch.cuh - sphere - GPU-based 3D discrete element method algorithm with optional fluid coupling
HTML git clone git://src.adamsgaard.dk/sphere
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
tcontactsearch.cuh (31217B)
---
1 #ifndef CONTACTSEARCH_CUH_
2 #define CONTACTSEARCH_CUH_
3
4 // contactsearch.cuh
5 // Functions for identifying contacts and processing boundaries
6
7
8 // Calculate the distance modifier between a contact pair. The modifier
9 // accounts for periodic boundaries. If the target cell lies outside
10 // the grid, it returns -1.
11 // Function is called from overlapsInCell() and findContactsInCell().
12 __device__ int findDistMod(int3* targetCell, Float3* distmod)
13 {
14 // Check whether x- and y boundaries are to be treated as periodic
15 // 1: x- and y boundaries periodic
16 // 2: x boundaries periodic
17 if (devC_grid.periodic == 1) {
18
19 // Periodic x-boundary
20 if (targetCell->x < 0) {
21 //targetCell->x = devC_grid.num[0] - 1;
22 targetCell->x += devC_grid.num[0];
23 *distmod += MAKE_FLOAT3(devC_grid.L[0], 0.0, 0.0);
24 }
25 if (targetCell->x >= devC_grid.num[0]) {
26 //targetCell->x = 0;
27 targetCell->x -= devC_grid.num[0];
28 *distmod -= MAKE_FLOAT3(devC_grid.L[0], 0.0, 0.0);
29 }
30
31 // Periodic y-boundary
32 if (targetCell->y < 0) {
33 //targetCell->y = devC_grid.num[1] - 1;
34 targetCell->y += devC_grid.num[1];
35 *distmod += MAKE_FLOAT3(0.0, devC_grid.L[1], 0.0);
36 }
37 if (targetCell->y >= devC_grid.num[1]) {
38 //targetCell->y = 0;
39 targetCell->y -= devC_grid.num[1];
40 *distmod -= MAKE_FLOAT3(0.0, devC_grid.L[1], 0.0);
41 }
42
43
44 }
45 // Only x-boundaries are periodic
46 if (devC_grid.periodic == 2) {
47
48 // Periodic x-boundary
49 if (targetCell->x < 0) {
50 //targetCell->x = devC_grid.num[0] - 1;
51 targetCell->x += devC_grid.num[0];
52 *distmod += MAKE_FLOAT3(devC_grid.L[0], 0.0, 0.0);
53 }
54 if (targetCell->x >= devC_grid.num[0]) {
55 //targetCell->x = 0;
56 targetCell->x -= devC_grid.num[0];
57 *distmod -= MAKE_FLOAT3(devC_grid.L[0], 0.0, 0.0);
58 }
59
60 // Hande out-of grid cases on y-axis
61 if (targetCell->y < 0 || targetCell->y >= devC_grid.num[1])
62 return -1;
63
64
65 }
66 // No periodic boundaries
67 if (devC_grid.periodic > 2 || devC_grid.periodic < 1) {
68
69 // Don't modify targetCell or distmod.
70
71 // Hande out-of grid cases on x- and y-axes
72 if (targetCell->x < 0 || targetCell->x >= devC_grid.num[0])
73 return -1;
74 if (targetCell->y < 0 || targetCell->y >= devC_grid.num[1])
75 return -1;
76 }
77
78 // Handle out-of-grid cases on z-axis
79 if (targetCell->z < 0 || targetCell->z >= devC_grid.num[2])
80 return -1;
81 else
82 return 0;
83 }
84
85
86 // Find overlaps between particle no. 'idx' and particles in cell 'gridpos'.
87 // Contacts are processed as soon as they are identified.
88 // Used for contactmodel=1, where contact history is not needed.
89 // Kernel executed on device, and callable from device only.
90 // Function is called from interact().
91 __device__ void findAndProcessContactsInCell(
92 int3 targetCell,
93 const unsigned int idx_a,
94 const Float4 x_a,
95 const Float radius_a,
96 Float3* F,
97 Float3* T,
98 Float* es_dot,
99 Float* ev_dot,
100 Float* p,
101 const Float4* __restrict__ dev_x_sorted,
102 const Float4* __restrict__ dev_vel_sorted,
103 const Float4* __restrict__ dev_angvel_sorted,
104 const unsigned int* __restrict__ dev_cellStart,
105 const unsigned int* __restrict__ dev_cellEnd,
106 const Float4* __restrict__ dev_walls_nx,
107 Float4* __restrict__ dev_walls_mvfd)
108 //uint4 bonds)
109 {
110
111 // Get distance modifier for interparticle
112 // vector, if it crosses a periodic boundary
113 Float3 distmod = MAKE_FLOAT3(0.0f, 0.0f, 0.0f);
114 if (findDistMod(&targetCell, &distmod) != -1) {
115
116
117 //// Check and process particle-particle collisions
118
119 // Calculate linear cell ID
120 unsigned int cellID = targetCell.x + targetCell.y * devC_grid.num[0]
121 + (devC_grid.num[0] * devC_grid.num[1]) * targetCell.z;
122
123 // Lowest particle index in cell
124 unsigned int startIdx = dev_cellStart[cellID];
125
126 // Make sure cell is not empty
127 if (startIdx != 0xffffffff) {
128
129 // Highest particle index in cell + 1
130 unsigned int endIdx = dev_cellEnd[cellID];
131
132 // Iterate over cell particles
133 for (unsigned int idx_b = startIdx; idx_b<endIdx; ++idx_b) {
134 if (idx_b != idx_a) { // Do not collide particle with itself
135
136 // Fetch position and velocity of particle B.
137 Float4 x_b = dev_x_sorted[idx_b];
138 Float radius_b = x_b.w;
139 Float kappa = devC_params.kappa;
140
141 // Distance between particle centers (Float4 -> Float3)
142 Float3 x_ab = MAKE_FLOAT3(x_a.x - x_b.x,
143 x_a.y - x_b.y,
144 x_a.z - x_b.z);
145
146 // Adjust interparticle vector if periodic boundary/boundaries
147 // are crossed
148 x_ab += distmod;
149 Float x_ab_length = length(x_ab);
150
151 // Distance between particle perimeters
152 Float delta_ab = x_ab_length - (radius_a + radius_b);
153
154 // Check for particle overlap
155 if (delta_ab < 0.0f) {
156 contactLinearViscous(F, T, es_dot, ev_dot, p,
157 idx_a, idx_b,
158 dev_vel_sorted,
159 dev_angvel_sorted,
160 radius_a, radius_b,
161 x_ab, x_ab_length,
162 delta_ab, kappa);
163 } else if (delta_ab < devC_params.db) {
164 // Check wether particle distance satisfies the
165 // capillary bond distance
166 capillaryCohesion_exp(F, radius_a, radius_b, delta_ab,
167 x_ab, x_ab_length, kappa);
168 }
169
170 // Check wether particles are bonded together
171 /*if (bonds.x == idx_b || bonds.y == idx_b ||
172 bonds.z == idx_b || bonds.w == idx_b) {
173 bondLinear(F, T, es_dot, p, % ev_dot missing
174 idx_a, idx_b,
175 dev_x_sorted, dev_vel_sorted,
176 dev_angvel_sorted,
177 radius_a, radius_b,
178 x_ab, x_ab_length,
179 delta_ab);
180 }*/
181
182 } // Do not collide particle with itself end
183 } // Iterate over cell particles end
184 } // Check wether cell is empty end
185 } // Periodic boundary and distance adjustment end
186 } // End of overlapsInCell(...)
187
188
189 // Find overlaps between particle no. 'idx' and particles in cell 'gridpos'
190 // Write the indexes of the overlaps in array contacts[].
191 // Used for contactmodel=2, where bookkeeping of contact history is necessary.
192 // Kernel executed on device, and callable from device only.
193 // Function is called from topology().
194 __device__ void findContactsInCell(
195 int3 targetCell,
196 const unsigned int idx_a,
197 const Float4 x_a,
198 const Float radius_a,
199 const Float4* __restrict__ dev_x_sorted,
200 const unsigned int* __restrict__ dev_cellStart,
201 const unsigned int* __restrict__ dev_cellEnd,
202 const unsigned int* __restrict__ dev_gridParticleIndex,
203 int* nc,
204 unsigned int* __restrict__ dev_contacts,
205 Float4* __restrict__ dev_distmod)
206 {
207 // Get distance modifier for interparticle
208 // vector, if it crosses a periodic boundary
209 Float3 distmod = MAKE_FLOAT3(0.0f, 0.0f, 0.0f);
210 if (findDistMod(&targetCell, &distmod) != -1) {
211
212 //// Check and process particle-particle collisions
213
214 // Calculate linear cell ID
215 unsigned int cellID = targetCell.x + targetCell.y * devC_grid.num[0]
216 + (devC_grid.num[0] * devC_grid.num[1]) * targetCell.z;
217
218 // Lowest particle index in cell
219 unsigned int startIdx = dev_cellStart[cellID];
220
221 // Make sure cell is not empty
222 if (startIdx != 0xffffffff) {
223
224 // Highest particle index in cell + 1
225 unsigned int endIdx = dev_cellEnd[cellID];
226
227 // Read the original index of particle A
228 unsigned int idx_a_orig = dev_gridParticleIndex[idx_a];
229
230 // Iterate over cell particles
231 for (unsigned int idx_b = startIdx; idx_b<endIdx; ++idx_b) {
232 if (idx_b != idx_a) { // Do not collide particle with itself
233
234
235 // Fetch position and radius of particle B.
236 Float4 x_b = dev_x_sorted[idx_b];
237 Float radius_b = x_b.w;
238
239 // Read the original index of particle B
240 unsigned int idx_b_orig = dev_gridParticleIndex[idx_b];
241
242 //__syncthreads();
243
244 // Distance between particle centers (Float4 -> Float3)
245 Float3 x_ab = MAKE_FLOAT3(x_a.x - x_b.x,
246 x_a.y - x_b.y,
247 x_a.z - x_b.z);
248
249 // Adjust interparticle vector if periodic boundary/boundaries
250 // are crossed
251 x_ab += distmod;
252
253 Float x_ab_length = length(x_ab);
254
255 // Distance between particle perimeters
256 Float delta_ab = x_ab_length - (radius_a + radius_b);
257
258 // Check for particle overlap
259 if (delta_ab < 0.0f) {
260
261 // If the particle is not yet registered in the contact list,
262 // use the next position in the array
263 int cpos = *nc;
264 unsigned int cidx;
265
266 // Find out, if particle is already registered in contact list
267 for (int i=0; i < devC_nc; ++i) {
268 __syncthreads();
269 cidx = dev_contacts[(unsigned int)(idx_a_orig*devC_nc+i)];
270 if (cidx == devC_np) // Write to position of now-deleted contact
271 cpos = i;
272 else if (cidx == idx_b_orig) { // Write to position of same contact
273 cpos = i;
274 break;
275 }
276 }
277
278 __syncthreads();
279
280 // Write the particle index to the relevant position,
281 // no matter if it already is there or not (concurrency of write)
282 dev_contacts[(unsigned int)(idx_a_orig*devC_nc+cpos)] = idx_b_orig;
283
284
285 // Write the interparticle vector and radius of particle B
286 //dev_x_ab_r_b[(unsigned int)(idx_a_orig*devC_nc+cpos)] = make_Float4(x_ab, radius_b);
287 dev_distmod[(unsigned int)(idx_a_orig*devC_nc+cpos)] = MAKE_FLOAT4(distmod.x, distmod.y, distmod.z, radius_b);
288
289 // Increment contact counter
290 ++*nc;
291
292 // Check if the number of contacts of particle A
293 // exceeds the max. number of contacts per particle
294 if (*nc > devC_nc)
295 return; // I would like to throw an error, but it doesn't seem possible...
296
297 }
298
299 // Write the inter-particle position vector correction and radius of particle B
300 //dev_distmod[(unsigned int)(idx_a_orig*devC_nc+cpos)] = make_Float4(distmod, radius_b);
301
302 // Check wether particles are bonded together
303 /*if (bonds.x == idx_b || bonds.y == idx_b ||
304 bonds.z == idx_b || bonds.w == idx_b) {
305 bondLinear(F, T, es_dot, p, % ev_dot missing
306 idx_a, idx_b,
307 dev_x_sorted, dev_vel_sorted,
308 dev_angvel_sorted,
309 radius_a, radius_b,
310 x_ab, x_ab_length,
311 delta_ab);
312 }*/
313
314 } // Do not collide particle with itself end
315 } // Iterate over cell particles end
316 } // Check wether cell is empty end
317 } // Periodic boundary and distmod end
318 } // End of findContactsInCell(...)
319
320
321 // For a single particle:
322 // Search for neighbors to particle 'idx' inside the 27 closest cells,
323 // and save the contact pairs in global memory.
324 // Function is called from mainGPU loop.
325 __global__ void topology(
326 const unsigned int* __restrict__ dev_cellStart,
327 const unsigned int* __restrict__ dev_cellEnd,
328 const unsigned int* __restrict__ dev_gridParticleIndex,
329 const Float4* __restrict__ dev_x_sorted,
330 unsigned int* __restrict__ dev_contacts,
331 Float4* __restrict__ dev_distmod)
332 {
333 // Thread index equals index of particle A
334 unsigned int idx_a = blockIdx.x * blockDim.x + threadIdx.x;
335 if (idx_a < devC_np) {
336 // Fetch particle data in global read
337 Float4 x_a = dev_x_sorted[idx_a];
338 Float radius_a = x_a.w;
339
340 // Count the number of contacts in this time step
341 int nc = 0;
342
343 // Grid position of host and neighbor cells in uniform, cubic grid
344 int3 gridPos;
345 int3 targetPos;
346
347 // Calculate cell address in grid from position of particle
348 gridPos.x = floor((x_a.x - devC_grid.origo[0]) / (devC_grid.L[0]/devC_grid.num[0]));
349 gridPos.y = floor((x_a.y - devC_grid.origo[1]) / (devC_grid.L[1]/devC_grid.num[1]));
350 gridPos.z = floor((x_a.z - devC_grid.origo[2]) / (devC_grid.L[2]/devC_grid.num[2]));
351
352 // Find overlaps between particle no. idx and all particles
353 // from its own cell + 26 neighbor cells
354 for (int z_dim=-1; z_dim<2; ++z_dim) { // z-axis
355 for (int y_dim=-1; y_dim<2; ++y_dim) { // y-axis
356 for (int x_dim=-1; x_dim<2; ++x_dim) { // x-axis
357 targetPos = gridPos + make_int3(x_dim, y_dim, z_dim);
358 findContactsInCell(targetPos, idx_a, x_a, radius_a,
359 dev_x_sorted,
360 dev_cellStart, dev_cellEnd,
361 dev_gridParticleIndex,
362 &nc, dev_contacts, dev_distmod);
363 }
364 }
365 }
366 }
367 } // End of topology(...)
368
369
370 // For a single particle:
371 // If contactmodel=1:
372 // Search for neighbors to particle 'idx' inside the 27 closest cells,
373 // and compute the resulting force and torque on it.
374 // If contactmodel=2:
375 // Process contacts saved in dev_contacts by topology(), and compute
376 // the resulting force and torque on it.
377 // For all contactmodels:
378 // Collide with top- and bottom walls, save resulting force on upper wall.
379 // Kernel is executed on device, and is callable from host only.
380 // Function is called from mainGPU loop.
381 __global__ void interact(
382 const unsigned int* __restrict__ dev_gridParticleIndex, // in
383 const unsigned int* __restrict__ dev_cellStart, // in
384 const unsigned int* __restrict__ dev_cellEnd, // in
385 const Float4* __restrict__ dev_x, // in
386 const Float4* __restrict__ dev_x_sorted, // in
387 const Float4* __restrict__ dev_vel_sorted, // in
388 const Float4* __restrict__ dev_angvel_sorted, // in
389 const Float4* __restrict__ dev_vel, // in
390 const Float4* __restrict__ dev_angvel, // in
391 Float4* __restrict__ dev_force, // out
392 Float4* __restrict__ dev_torque, // out
393 Float* __restrict__ dev_es_dot, // out
394 Float* __restrict__ dev_ev_dot, // out
395 Float* __restrict__ dev_es, // out
396 Float* __restrict__ dev_ev, // out
397 Float* __restrict__ dev_p, // out
398 const Float4* __restrict__ dev_walls_nx, // in
399 Float4* __restrict__ dev_walls_mvfd, // in
400 Float* __restrict__ dev_walls_force_pp, // out
401 unsigned int* __restrict__ dev_contacts, // out
402 const Float4* __restrict__ dev_distmod, // in
403 Float4* __restrict__ dev_delta_t) // out
404 {
405 // Thread index equals index of particle A
406 unsigned int idx_a = blockIdx.x * blockDim.x + threadIdx.x;
407
408 if (idx_a < devC_np) {
409
410 // Fetch particle data in global read
411 unsigned int idx_a_orig = dev_gridParticleIndex[idx_a];
412 Float4 x_a = dev_x_sorted[idx_a];
413 Float radius_a = x_a.w;
414
415 // Fetch world dimensions in constant memory read
416 Float3 origo = MAKE_FLOAT3(devC_grid.origo[0],
417 devC_grid.origo[1],
418 devC_grid.origo[2]);
419 Float3 L = MAKE_FLOAT3(devC_grid.L[0],
420 devC_grid.L[1],
421 devC_grid.L[2]);
422
423 // Fetch wall data in global read
424 Float4 w_0_nx, w_1_nx, w_2_nx, w_3_nx, w_4_nx;
425 Float4 w_0_mvfd, w_1_mvfd, w_2_mvfd, w_3_mvfd, w_4_mvfd;
426
427 // default wall normals and positions
428 w_0_nx = MAKE_FLOAT4( 0.0f, 0.0f,-1.0f, L.z);
429 w_1_nx = MAKE_FLOAT4(-1.0f, 0.0f, 0.0f, L.x);
430 w_2_nx = MAKE_FLOAT4( 1.0f, 0.0f, 0.0f, 0.0f);
431 w_3_nx = MAKE_FLOAT4( 0.0f,-1.0f, 0.0f, L.y);
432 w_4_nx = MAKE_FLOAT4( 0.0f, 1.0f, 0.0f, 0.0f);
433
434 // default wall mass, vel, force, and sigma0
435 w_0_mvfd = MAKE_FLOAT4(0.0f, 0.0f, 0.0f, 0.0f);
436 w_1_mvfd = MAKE_FLOAT4(0.0f, 0.0f, 0.0f, 0.0f);
437 w_2_mvfd = MAKE_FLOAT4(0.0f, 0.0f, 0.0f, 0.0f);
438 w_3_mvfd = MAKE_FLOAT4(0.0f, 0.0f, 0.0f, 0.0f);
439 w_4_mvfd = MAKE_FLOAT4(0.0f, 0.0f, 0.0f, 0.0f);
440
441 // fetch data for dynamic walls
442 if (devC_nw >= 1) {
443 w_0_nx = dev_walls_nx[0];
444 w_0_mvfd = dev_walls_mvfd[0];
445 if (devC_nw >= 2) {
446 w_1_nx = dev_walls_nx[1];
447 w_1_mvfd = dev_walls_mvfd[1];
448 }
449 if (devC_nw >= 3) {
450 w_2_nx = dev_walls_nx[2];
451 w_2_mvfd = dev_walls_mvfd[2];
452 }
453 if (devC_nw >= 4) {
454 w_3_nx = dev_walls_nx[3];
455 w_3_mvfd = dev_walls_mvfd[3];
456 }
457 if (devC_nw >= 5) {
458 w_4_nx = dev_walls_nx[4];
459 w_4_mvfd = dev_walls_mvfd[4];
460 }
461 }
462
463 // Index of particle which is bonded to particle A.
464 // The index is equal to the particle no (p.np)
465 // if particle A is bond-less.
466 //uint4 bonds = dev_bonds_sorted[idx_a];
467
468 // Initiate shear friction loss rate at 0.0
469 Float es_dot = 0.0;
470 Float ev_dot = 0.0;
471
472 // Initiate pressure on particle at 0.0
473 Float p = 0.0;
474
475 // Allocate memory for temporal force/torque vector values
476 Float3 F = MAKE_FLOAT3(0.0, 0.0, 0.0);
477 Float3 T = MAKE_FLOAT3(0.0, 0.0, 0.0);
478
479 // Apply linear elastic, frictional contact model to registered contacts
480 if (devC_params.contactmodel == 2 || devC_params.contactmodel == 3) {
481 unsigned int idx_b_orig, mempos;
482 Float delta_n, x_ab_length, radius_b;
483 Float3 x_ab;
484 Float4 x_b, distmod;
485 Float4 vel_a = dev_vel_sorted[idx_a];
486 Float4 angvel4_a = dev_angvel_sorted[idx_a];
487 Float3 angvel_a = MAKE_FLOAT3(
488 angvel4_a.x, angvel4_a.y, angvel4_a.z);
489
490 // Loop over all possible contacts, and remove contacts
491 // that no longer are valid (delta_n > 0.0)
492 for (int i = 0; i<devC_nc; ++i) {
493 mempos = (unsigned int)(idx_a_orig * devC_nc + i);
494 __syncthreads();
495 idx_b_orig = dev_contacts[mempos];
496
497 if (idx_b_orig != (unsigned int)devC_np) {
498
499 // Read inter-particle distance correction vector
500 distmod = dev_distmod[mempos];
501
502 // Read particle b position and radius
503 x_b = dev_x[idx_b_orig];
504 radius_b = x_b.w;
505
506 // Inter-particle vector, corrected for periodic boundaries
507 x_ab = MAKE_FLOAT3(x_a.x - x_b.x + distmod.x,
508 x_a.y - x_b.y + distmod.y,
509 x_a.z - x_b.z + distmod.z);
510
511 x_ab_length = length(x_ab);
512 delta_n = x_ab_length - (radius_a + radius_b);
513
514 // Process collision if the particles are overlapping
515 if (delta_n < 0.0) {
516 if (devC_params.contactmodel == 2) {
517 contactLinear(&F, &T, &es_dot, &ev_dot, &p,
518 idx_a_orig,
519 idx_b_orig,
520 vel_a,
521 dev_vel,
522 angvel_a,
523 dev_angvel,
524 radius_a, radius_b,
525 x_ab, x_ab_length,
526 delta_n, dev_delta_t,
527 mempos);
528 } else if (devC_params.contactmodel == 3) {
529 contactHertz(&F, &T, &es_dot, &ev_dot, &p,
530 idx_a_orig,
531 idx_b_orig,
532 vel_a,
533 dev_vel,
534 angvel_a,
535 dev_angvel,
536 radius_a, radius_b,
537 x_ab, x_ab_length,
538 delta_n, dev_delta_t,
539 mempos);
540 }
541 } else {
542 __syncthreads();
543 // Remove this contact (there is no particle with
544 // index=np)
545 dev_contacts[mempos] = devC_np;
546 // Zero sum of shear displacement in this position
547 dev_delta_t[mempos] = MAKE_FLOAT4(0.0, 0.0, 0.0, 0.0);
548 }
549
550 } else { // if dev_contacts[mempos] == devC_np
551 __syncthreads();
552 // Zero sum of shear displacement in this position
553 dev_delta_t[mempos] = MAKE_FLOAT4(0.0, 0.0, 0.0, 0.0);
554 }
555 } // Contact loop end
556
557
558 // Find contacts and process collisions immidiately for
559 // contactmodel 1 (visco-frictional).
560 } else if (devC_params.contactmodel == 1) {
561
562 int3 gridPos;
563 int3 targetPos;
564
565 // Calculate address in grid from position
566 gridPos.x = floor((x_a.x - devC_grid.origo[0])
567 / (devC_grid.L[0]/devC_grid.num[0]));
568 gridPos.y = floor((x_a.y - devC_grid.origo[1])
569 / (devC_grid.L[1]/devC_grid.num[1]));
570 gridPos.z = floor((x_a.z - devC_grid.origo[2])
571 / (devC_grid.L[2]/devC_grid.num[2]));
572
573 // Find overlaps between particle no. idx and all particles
574 // from its own cell + 26 neighbor cells.
575 // Calculate resulting normal- and shear-force components and
576 // torque for the particle on the base of contactLinearViscous()
577 for (int z_dim=-1; z_dim<2; ++z_dim) { // z-axis
578 for (int y_dim=-1; y_dim<2; ++y_dim) { // y-axis
579 for (int x_dim=-1; x_dim<2; ++x_dim) { // x-axis
580 targetPos = gridPos + make_int3(x_dim, y_dim, z_dim);
581 findAndProcessContactsInCell(targetPos, idx_a,
582 x_a, radius_a,
583 &F, &T, &es_dot,
584 &ev_dot, &p,
585 dev_x_sorted,
586 dev_vel_sorted,
587 dev_angvel_sorted,
588 dev_cellStart,
589 dev_cellEnd,
590 dev_walls_nx,
591 dev_walls_mvfd);
592 }
593 }
594 }
595
596 }
597
598 //// Interact with walls
599 Float delta_w; // Overlap distance
600 Float3 w_n; // Wall surface normal
601 Float w_0_force = 0.0; // Force on wall 0 from particle A
602 Float w_1_force = 0.0; // Force on wall 1 from particle A
603 Float w_2_force = 0.0; // Force on wall 2 from particle A
604 Float w_3_force = 0.0; // Force on wall 3 from particle A
605 Float w_4_force = 0.0; // Force on wall 4 from particle A
606
607 // Upper wall (idx 0)
608 delta_w = w_0_nx.w - (x_a.z + radius_a);
609 w_n = MAKE_FLOAT3(w_0_nx.x, w_0_nx.y, w_0_nx.z);
610 if (delta_w < 0.0f) {
611 w_0_force = contactLinear_wall(&F, &T, &es_dot, &ev_dot, &p, idx_a,
612 radius_a, dev_vel_sorted,
613 dev_angvel_sorted, w_n,
614 delta_w, w_0_mvfd.y);
615 }
616
617 // Lower wall (force on wall not stored)
618 delta_w = x_a.z - radius_a - origo.z;
619 w_n = MAKE_FLOAT3(0.0f, 0.0f, 1.0f);
620 if (delta_w < 0.0f) {
621 (void)contactLinear_wall(&F, &T, &es_dot, &ev_dot, &p, idx_a,
622 radius_a, dev_vel_sorted,
623 dev_angvel_sorted, w_n, delta_w,
624 0.0);
625 }
626
627
628 if (devC_grid.periodic == 0) { // no periodic walls
629
630 // Right wall (idx 1)
631 delta_w = w_1_nx.w - (x_a.x + radius_a);
632 w_n = MAKE_FLOAT3(w_1_nx.x, w_1_nx.y, w_1_nx.z);
633 if (delta_w < 0.0f) {
634 w_1_force = contactLinear_wall(&F, &T, &es_dot, &ev_dot, &p,
635 idx_a, radius_a,
636 dev_vel_sorted,
637 dev_angvel_sorted, w_n,
638 delta_w, w_1_mvfd.y);
639 }
640
641 // Left wall (idx 2)
642 delta_w = x_a.x - radius_a - w_2_nx.w;
643 w_n = MAKE_FLOAT3(w_2_nx.x, w_2_nx.y, w_2_nx.z);
644 if (delta_w < 0.0f) {
645 w_2_force = contactLinear_wall(&F, &T, &es_dot, &ev_dot, &p,
646 idx_a, radius_a,
647 dev_vel_sorted,
648 dev_angvel_sorted, w_n,
649 delta_w, w_2_mvfd.y);
650 }
651
652 // Back wall (idx 3)
653 delta_w = w_3_nx.w - (x_a.y + radius_a);
654 w_n = MAKE_FLOAT3(w_3_nx.x, w_3_nx.y, w_3_nx.z);
655 if (delta_w < 0.0f) {
656 w_3_force = contactLinear_wall(&F, &T, &es_dot, &ev_dot, &p,
657 idx_a, radius_a,
658 dev_vel_sorted,
659 dev_angvel_sorted, w_n,
660 delta_w, w_3_mvfd.y);
661 }
662
663 // Front wall (idx 4)
664 delta_w = x_a.y - radius_a - w_4_nx.w;
665 w_n = MAKE_FLOAT3(w_4_nx.x, w_4_nx.y, w_4_nx.z);
666 if (delta_w < 0.0f) {
667 w_4_force = contactLinear_wall(&F, &T, &es_dot, &ev_dot, &p,
668 idx_a, radius_a,
669 dev_vel_sorted,
670 dev_angvel_sorted, w_n,
671 delta_w, w_4_mvfd.y);
672 }
673
674 } else if (devC_grid.periodic == 2) { // right and left walls periodic
675
676 // Back wall (idx 3)
677 delta_w = w_3_nx.w - (x_a.y + radius_a);
678 w_n = MAKE_FLOAT3(w_3_nx.x, w_3_nx.y, w_3_nx.z);
679 if (delta_w < 0.0f) {
680 w_3_force = contactLinear_wall(&F, &T, &es_dot, &ev_dot, &p,
681 idx_a, radius_a,
682 dev_vel_sorted,
683 dev_angvel_sorted, w_n,
684 delta_w, w_3_mvfd.y);
685 }
686
687 // Front wall (idx 4)
688 delta_w = x_a.y - radius_a - w_4_nx.w;
689 w_n = MAKE_FLOAT3(w_4_nx.x, w_4_nx.y, w_4_nx.z);
690 if (delta_w < 0.0f) {
691 w_4_force = contactLinear_wall(&F, &T, &es_dot, &ev_dot, &p,
692 idx_a, radius_a,
693 dev_vel_sorted,
694 dev_angvel_sorted, w_n,
695 delta_w, w_4_mvfd.y);
696 }
697 }
698
699 // Hold threads for coalesced write
700 __syncthreads();
701
702 // Write force to unsorted position
703 unsigned int orig_idx = dev_gridParticleIndex[idx_a];
704 dev_force[orig_idx] = MAKE_FLOAT4(F.x, F.y, F.z, 0.0f);
705 dev_torque[orig_idx] = MAKE_FLOAT4(T.x, T.y, T.z, 0.0f);
706 dev_es_dot[orig_idx] = es_dot;
707 dev_ev_dot[orig_idx] = ev_dot;
708 dev_es[orig_idx] += es_dot * devC_dt;
709 dev_ev[orig_idx] += ev_dot * devC_dt;
710 dev_p[orig_idx] = p;
711 if (devC_nw > 0)
712 dev_walls_force_pp[orig_idx] = w_0_force;
713 if (devC_nw > 1)
714 dev_walls_force_pp[orig_idx+devC_np] = w_1_force;
715 if (devC_nw > 2)
716 dev_walls_force_pp[orig_idx+devC_np*2] = w_2_force;
717 if (devC_nw > 3)
718 dev_walls_force_pp[orig_idx+devC_np*3] = w_3_force;
719 if (devC_nw > 4)
720 dev_walls_force_pp[orig_idx+devC_np*4] = w_4_force;
721 }
722 } // End of interact(...)
723
724
725 #endif
726 // vim: tabstop=8 expandtab shiftwidth=4 softtabstop=4