PeriDyno 1.0.0
Loading...
Searching...
No Matches
RigidSandCoupling.cu
Go to the documentation of this file.
1#include "RigidSandCoupling.h"
2
3#include "Math/Lerp.h"
4
5#include "Primitive/Primitive3D.h"
6
7namespace dyno
8{
9 template<typename TDataType>
10 RigidSandCoupling<TDataType>::RigidSandCoupling()
11 : Node()
12 {
13 }
14
15 template<typename TDataType>
16 RigidSandCoupling<TDataType>::~RigidSandCoupling()
17 {
18 }
19
20 template<typename TDataType>
21 void RigidSandCoupling<TDataType>::resetStates()
22 {
23 }
24
25 template<typename Real, typename Coord3D, typename Coord4D, typename Box3D>
26 __global__ void CQWS_CoupleSandWithBoxes(
27 DArray2D<Coord4D> grid,
28 DArray2D<Coord3D> displacements,
29 DArray<Box3D> boxes,
30 DArray<Coord3D> boxVel,
31 DArray<Coord3D> boxAngularVel,
32 Coord3D origin,
33 Real spacing,
34 uint offset)
35 {
36 int i = threadIdx.x + blockIdx.x * blockDim.x;
37 int j = threadIdx.y + blockIdx.y * blockDim.y;
38
39 uint nx = displacements.nx();
40 uint ny = displacements.ny();
41
42 if (i < nx && j < ny)
43 {
44 Coord3D v = origin + spacing * Coord3D(i, 0, j) + displacements(i, j);
45 Point3D p = Point3D(v);
46
47 uint nb = boxes.size();
48 int boxId = -1;
49 Real Dmin = 100000.0f;
50 for (uint k = 0; k < nb; k++)
51 {
52 Real d = p.distance(boxes[k]);
53 if (d < Dmin)
54 {
55 Dmin = d;
56 boxId = k;
57 }
58 }
59
60 if (Dmin < 0 && boxId != -1)
61 {
62 Coord4D gij = grid(i + 1, j + 1);
63
64 Coord3D c = boxes[boxId].center;
65
66 Coord3D v_box = boxVel[offset + boxId];
67 Coord3D omega_box = boxAngularVel[offset + boxId];
68
69 Coord3D vv = v_box + omega_box.cross(v - c);
70
71 Real duh = gij.x * vv.x;
72 Real dvh = gij.x * vv.z;
73 Real dw = abs(Dmin) * 5;
74 atomicAdd(&grid(i + 1, j + 1).y, duh);
75 atomicAdd(&grid(i + 1, j + 1).z, dvh);
76
77 atomicAdd(&grid(i + 2, j + 1).y, 0.25 * dw);
78 atomicAdd(&grid(i + 0, j + 1).y, -0.25 * dw);
79
80 atomicAdd(&grid(i + 1, j + 2).z, 0.25 * dw);
81 atomicAdd(&grid(i + 1, j + 0).z, -0.25 * dw);
82 }
83 }
84 }
85
86
87 template<typename Real, typename Coord3D, typename Coord4D, typename Sphere3D>
88 __global__ void CQWS_CoupleSandAndSpheres(
89 DArray2D<Coord4D> grid,
90 DArray2D<Coord3D> displacements,
91 DArray<Sphere3D> spheres,
92 DArray<Coord3D> ele_vel,
93 DArray<Coord3D> ele_vel_angular,
94 Coord3D origin,
95 Real spacing,
96 uint offset)
97 {
98 int i = threadIdx.x + blockIdx.x * blockDim.x;
99 int j = threadIdx.y + blockIdx.y * blockDim.y;
100
101 uint nx = displacements.nx();
102 uint ny = displacements.ny();
103
104 if (i < nx && j < ny)
105 {
106 Coord3D v = origin + spacing * Coord3D(i, 0, j) + displacements(i, j);
107 Point3D p = Point3D(v);
108
109 uint nb = spheres.size();
110 int sphereId = -1;
111 Real Dmin = 100000.0f;
112 for (uint k = 0; k < nb; k++)
113 {
114 Real d = p.distance(spheres[k]);
115 if (d < Dmin)
116 {
117 Dmin = d;
118 sphereId = k;
119 }
120 }
121
122 if (Dmin < 0 && sphereId != -1)
123 {
124 Coord4D gij = grid(i + 1, j + 1);
125
126 Coord3D c = spheres[sphereId].center;
127
128 Coord3D v_i = ele_vel[offset + sphereId];
129 Coord3D omega_i = ele_vel_angular[offset + sphereId];
130
131 Coord3D vv = v_i + omega_i.cross(v - c);
132
133 Real duh = gij.x * vv.x;
134 Real dvh = gij.x * vv.z;
135 Real dw = abs(Dmin) * 5;
136 atomicAdd(&grid(i + 1, j + 1).y, duh);
137 atomicAdd(&grid(i + 1, j + 1).z, dvh);
138
139 atomicAdd(&grid(i + 2, j + 1).y, 0.25 * dw);
140 atomicAdd(&grid(i + 0, j + 1).y, -0.25 * dw);
141
142 atomicAdd(&grid(i + 1, j + 2).z, 0.25 * dw);
143 atomicAdd(&grid(i + 1, j + 0).z, -0.25 * dw);
144 }
145 }
146 }
147
148 template<typename Real, typename Coord3D, typename Coord4D, typename Capsule3D>
149 __global__ void CQWS_CoupleSandAndCapsules(
150 DArray2D<Coord4D> grid,
151 DArray2D<Coord3D> displacements,
152 DArray<Capsule3D> capsules,
153 DArray<Coord3D> ele_vel,
154 DArray<Coord3D> ele_vel_angular,
155 Coord3D origin,
156 Real spacing,
157 uint offset)
158 {
159 int i = threadIdx.x + blockIdx.x * blockDim.x;
160 int j = threadIdx.y + blockIdx.y * blockDim.y;
161
162 uint nx = displacements.nx();
163 uint ny = displacements.ny();
164
165 if (i < nx && j < ny)
166 {
167 Coord3D v = origin + spacing * Coord3D(i, 0, j) + displacements(i, j);
168 Point3D p = Point3D(v);
169
170 uint nb = capsules.size();
171 int sphereId = -1;
172 Real Dmin = 100000.0f;
173 for (uint k = 0; k < nb; k++)
174 {
175 Real d = p.distance(capsules[k]);
176 if (d < Dmin)
177 {
178 Dmin = d;
179 sphereId = k;
180 }
181 }
182
183 if (Dmin < 0 && sphereId != -1)
184 {
185 Coord4D gij = grid(i + 1, j + 1);
186
187 Coord3D c = capsules[sphereId].center;
188
189 Coord3D v_i = ele_vel[offset + sphereId];
190 Coord3D omega_i = ele_vel_angular[offset + sphereId];
191
192 Coord3D vv = v_i + omega_i.cross(v - c);
193
194 Real duh = gij.x * vv.x;
195 Real dvh = gij.x * vv.z;
196 Real dw = abs(Dmin) * 5;
197 atomicAdd(&grid(i + 1, j + 1).y, duh);
198 atomicAdd(&grid(i + 1, j + 1).z, dvh);
199
200 atomicAdd(&grid(i + 2, j + 1).y, 0.25 * dw);
201 atomicAdd(&grid(i + 0, j + 1).y, -0.25 * dw);
202
203 atomicAdd(&grid(i + 1, j + 2).z, 0.25 * dw);
204 atomicAdd(&grid(i + 1, j + 0).z, -0.25 * dw);
205 }
206 }
207 }
208
209 template<typename TDataType>
210 void RigidSandCoupling<TDataType>::updateStates()
211 {
212 Real dt = this->stateTimeStep()->getData();
213
214 auto sand = this->getGranularMedia();
215 auto rigidbody = this->getRigidBodySystem();
216
217 auto heights = sand->stateHeightField()->constDataPtr();
218 auto elements = rigidbody->stateTopology()->constDataPtr();
219
220 auto& grid2d = sand->stateGrid()->getData();
221
222 auto& ele_vel = rigidbody->stateVelocity()->getData();
223 auto& ele_angular_vel = rigidbody->stateAngularVelocity()->getData();
224
225 auto& boxes = elements->boxesInGlobal();
226 auto& spheres = elements->spheresInGlobal();
227 auto& capsules = elements->capsulesInGlobal();
228
229 auto offset = elements->calculateElementOffset();
230
231 auto& disp = heights->getDisplacement();
232 auto origin = heights->getOrigin();
233 auto spacing = heights->getGridSpacing();
234
235 uint boxOffset = offset.boxIndex();
236 uint sphereOffset = offset.sphereIndex();
237 uint capsuleOffset = offset.capsuleIndex();
238
239 cuExecute2D(make_uint2(disp.nx(), disp.ny()),
240 CQWS_CoupleSandWithBoxes,
241 grid2d,
242 disp,
243 boxes,
244 ele_vel,
245 ele_angular_vel,
246 origin,
247 spacing,
248 boxOffset);
249
250 cuExecute2D(make_uint2(disp.nx(), disp.ny()),
251 CQWS_CoupleSandAndSpheres,
252 grid2d,
253 disp,
254 spheres,
255 ele_vel,
256 ele_angular_vel,
257 origin,
258 spacing,
259 sphereOffset);
260
261 cuExecute2D(make_uint2(disp.nx(), disp.ny()),
262 CQWS_CoupleSandAndCapsules,
263 grid2d,
264 disp,
265 capsules,
266 ele_vel,
267 ele_angular_vel,
268 origin,
269 spacing,
270 capsuleOffset);
271 }
272
273 DEFINE_CLASS(RigidSandCoupling);
274}
Generated by doxygen 1.13.2