PeriDyno 1.0.0
Loading...
Searching...
No Matches
CapillaryWave.cu
Go to the documentation of this file.
1#include "CapillaryWave.h"
2
3#include "SceneGraph.h"
4
5#include "Mapping/HeightFieldToTriangleSet.h"
6#include "GLSurfaceVisualModule.h"
7
8namespace dyno
9{
10 template<typename TDataType>
11 CapillaryWave<TDataType>::CapillaryWave()
12 : Node()
13 {
14 auto heights = std::make_shared<HeightField<TDataType>>();
15 this->stateHeightField()->setDataPtr(heights);
16
17 auto mapper = std::make_shared<HeightFieldToTriangleSet<DataType3f>>();
18 this->stateHeightField()->connect(mapper->inHeightField());
19 this->graphicsPipeline()->pushModule(mapper);
20
21 // mapper->varScale()->setValue(0.1);
22 // mapper->varTranslation()->setValue(Vec3f(-2, 0.2, -2));
23
24 auto sRender = std::make_shared<GLSurfaceVisualModule>();
25 sRender->setColor(Color(0, 0.2, 1.0));
26 mapper->outTriangleSet()->connect(sRender->inTriangleSet());
27 this->graphicsPipeline()->pushModule(sRender);
28 }
29
30 template<typename TDataType>
31 CapillaryWave<TDataType>::~CapillaryWave()
32 {
33 mDeviceGrid.clear();
34 mDeviceGridNext.clear();
35 }
36
37 template <typename Coord3D, typename Coord4D>
38 __global__ void CW_UpdateHeightDisp(
39 DArray2D<Coord3D> displacement,
40 DArray2D<Coord4D> dis)
41 {
42 unsigned int i = blockIdx.x * blockDim.x + threadIdx.x;
43 unsigned int j = blockIdx.y * blockDim.y + threadIdx.y;
44 if (i < displacement.nx() && j < displacement.ny())
45 {
46 displacement(i, j).y = dis(i, j).x;
47 }
48 }
49
50 template <typename Coord>
51 __device__ float C_GetU(Coord gp)
52 {
53 Real h = max(gp.x, 0.0f);
54 Real uh = gp.y;
55
56 Real h4 = h * h * h * h;
57 return sqrtf(2.0f) * h * uh / (sqrtf(h4 + max(h4, EPSILON)));
58 }
59
60 template <typename Coord>
61 __device__ Real C_GetV(Coord gp)
62 {
63 Real h = max(gp.x, 0.0f);
64 Real vh = gp.z;
65
66 Real h4 = h * h * h * h;
67 return sqrtf(2.0f) * h * vh / (sqrtf(h4 + max(h4, EPSILON)));
68 }
69
70 template <typename Coord>
71 __global__ void CW_MoveSimulatedRegion(
72 DArray2D<Coord> grid_next,
73 DArray2D<Coord> grid,
74 int width,
75 int height,
76 int dx,
77 int dy,
78 Real level)
79 {
80 int i = threadIdx.x + blockIdx.x * blockDim.x;
81 int j = threadIdx.y + blockIdx.y * blockDim.y;
82 if (i < width && j < height)
83 {
84 int gx = i + 1;
85 int gy = j + 1;
86
87 Coord gp = grid(gx, gy);
88 Coord gp_init = Coord(level, 0.0f, 0.0f, gp.w);
89
90 int new_i = i - dx;
91 int new_j = j - dy;
92
93 if (new_i < 0 || new_i >= width) gp = gp_init;
94
95 new_i = new_i % width;
96 if (new_i < 0) new_i = width + new_i;
97
98 if (new_j < 0 || new_j >= height) gp = gp_init;
99
100 new_j = new_j % height;
101 if (new_j < 0) new_j = height + new_j;
102
103 grid(new_i + 1, new_j + 1) = gp;
104 }
105 }
106
107 template<typename TDataType>
108 void CapillaryWave<TDataType>::moveDynamicRegion(int nx, int ny)
109 {
110 auto res = this->varResolution()->getValue();
111
112 auto level = this->varWaterLevel()->getValue();
113
114 int extNx = res + 2;
115 int extNy = res + 2;
116
117 cuExecute2D(make_uint2(extNx, extNy),
118 CW_MoveSimulatedRegion,
119 mDeviceGridNext,
120 mDeviceGrid,
121 res,
122 res,
123 nx,
124 ny,
125 level);
126
127 mOriginX += nx;
128 mOriginY += ny;
129 }
130
131 template<typename TDataType>
132 void CapillaryWave<TDataType>::resetStates()
133 {
134 int res = this->varResolution()->getValue();
135 Real length = this->varLength()->getValue();
136
137 Real level = this->varWaterLevel()->getValue();
138
139 mRealGridSize = length / res;
140
141 int extNx = res + 2;
142 int extNy = res + 2;
143
144 mDeviceGrid.resize(extNx, extNy);
145 mDeviceGridNext.resize(extNx, extNy);
146 this->stateHeight()->resize(res, res);
147
148 //init grid with initial values
149 cuExecute2D(make_uint2(extNx, extNy),
150 InitDynamicRegion,
151 mDeviceGrid,
152 extNx,
153 extNy,
154 level);
155
156 //init grid_next with initial values
157 cuExecute2D(make_uint2(extNx, extNy),
158 InitDynamicRegion,
159 mDeviceGridNext,
160 extNx,
161 extNy,
162 level);
163
164 auto topo = this->stateHeightField()->getDataPtr();
165 topo->setExtents(res, res);
166 topo->setGridSpacing(mRealGridSize);
167 topo->setOrigin(Coord3D(-0.5 * mRealGridSize * topo->width(), 0, -0.5 * mRealGridSize * topo->height()));
168
169 auto& disp = topo->getDisplacement();
170
171 uint2 extent;
172 extent.x = disp.nx();
173 extent.y = disp.ny();
174
175 cuExecute2D(extent,
176 CW_InitHeightDisp,
177 this->stateHeight()->getData(),
178 disp,
179 mDeviceGrid,
180 level);
181 }
182
183 template<typename TDataType>
184 void CapillaryWave<TDataType>::updateStates()
185 {
186 Real dt = this->stateTimeStep()->getValue();
187
188 Real level = this->varWaterLevel()->getValue();
189
190 uint res = this->varResolution()->getValue();
191
192 int extNx = res + 2;
193 int extNy = res + 2;
194
195 int nStep = 1;
196 float timestep = dt / nStep;
197
198 auto scn = this->getSceneGraph();
199 auto GRAVITY = scn->getGravity().norm();
200
201 for (int iter = 0; iter < nStep; iter++)
202 {
203 cuExecute2D(make_uint2(extNx, extNy),
204 CW_ImposeBC,
205 mDeviceGridNext,
206 mDeviceGrid,
207 extNx,
208 extNy);
209
210 cuExecute2D(make_uint2(res, res),
211 CW_OneWaveStep,
212 mDeviceGrid,
213 mDeviceGridNext,
214 res,
215 res,
216 GRAVITY,
217 timestep);
218 }
219
220 cuExecute2D(make_uint2(res, res),
221 CW_InitHeights,
222 this->stateHeight()->getData(),
223 mDeviceGrid,
224 res,
225 mRealGridSize);
226
227 cuExecute2D(make_uint2(res, res),
228 CW_InitHeightGrad,
229 this->stateHeight()->getData(),
230 res);
231
232 //Update topology
233 auto topo = this->stateHeightField()->getDataPtr();
234
235 auto& disp = topo->getDisplacement();
236
237 uint2 extent;
238 extent.x = disp.nx();
239 extent.y = disp.ny();
240
241 cuExecute2D(extent,
242 CW_UpdateHeightDisp,
243 disp,
244 this->stateHeight()->getData());
245 }
246
247 template <typename Coord4D>
248 __global__ void InitDynamicRegion(DArray2D<Coord4D> grid, int gridwidth, int gridheight, float level)
249 {
250 int x = threadIdx.x + blockIdx.x * blockDim.x;
251 int y = threadIdx.y + blockIdx.y * blockDim.y;
252 if (x < gridwidth && y < gridheight)
253 {
254 Coord4D gp;
255 gp.x = level;
256 gp.y = 0.0f;
257 gp.z = 0.0f;
258 gp.w = 0.0f;
259
260 grid(x, y) = gp;
261// if ((x - 256) * (x - 256) + (y - 256) * (y - 256) <= 2500) grid(x, y).x = level;
262 }
263 }
264
265 template <typename Coord4D>
266 __global__ void CW_ImposeBC(
267 DArray2D<Coord4D> grid_next,
268 DArray2D<Coord4D> grid,
269 int width,
270 int height)
271 {
272 int x = threadIdx.x + blockIdx.x * blockDim.x;
273 int y = threadIdx.y + blockIdx.y * blockDim.y;
274 if (x < width && y < height)
275 {
276 if (x == 0)
277 {
278 Coord4D a = grid(1, y);
279 grid_next(x, y) = a;
280 }
281 else if (x == width - 1)
282 {
283 Coord4D a = grid(width - 2, y);
284 grid_next(x, y) = a;
285 }
286 else if (y == 0)
287 {
288 Coord4D a = grid(x, 1);
289 grid_next(x, y) = a;
290 }
291 else if (y == height - 1)
292 {
293 Coord4D a = grid(x, height - 2);
294 grid_next(x, y) = a;
295 }
296 else
297 {
298 Coord4D a = grid(x, y);
299 grid_next(x, y) = a;
300 }
301 }
302 }
303
304 template <typename Coord>
305 __host__ __device__ void CW_FixShore(Coord& l, Coord& c, Coord& r)
306 {
307
308 if (r.x < 0.0f || l.x < 0.0f || c.x < 0.0f)
309 {
310 c.x = c.x + l.x + r.x;
311 c.x = max(0.0f, c.x);
312 l.x = 0.0f;
313 r.x = 0.0f;
314 }
315 float h = c.x;
316 float h4 = h * h * h * h;
317 float v = sqrtf(2.0f) * h * c.y / (sqrtf(h4 + max(h4, EPSILON)));
318 float u = sqrtf(2.0f) * h * c.z / (sqrtf(h4 + max(h4, EPSILON)));
319
320 c.y = u * h;
321 c.z = v * h;
322 }
323
324 template <typename Coord>
325 __host__ __device__ Coord CW_VerticalPotential(Coord gp, float GRAVITY)
326 {
327 float h = max(gp.x, 0.0f);
328 float uh = gp.y;
329 float vh = gp.z;
330
331 float h4 = h * h * h * h;
332 float v = sqrtf(2.0f) * h * vh / (sqrtf(h4 + max(h4, EPSILON)));
333
334 Coord G;
335 G.x = v * h;
336 G.y = uh * v;
337 G.z = vh * v + GRAVITY * h * h;
338 G.w = 0.0f;
339 return G;
340 }
341
342 template <typename Coord>
343 __device__ Coord CW_HorizontalPotential(Coord gp, float GRAVITY)
344 {
345 float h = max(gp.x, 0.0f);
346 float uh = gp.y;
347 float vh = gp.z;
348
349 float h4 = h * h * h * h;
350 float u = sqrtf(2.0f) * h * uh / (sqrtf(h4 + max(h4, EPSILON)));
351
352 Coord F;
353 F.x = u * h;
354 F.y = uh * u + GRAVITY * h * h;
355 F.z = vh * u;
356 F.w = 0.0f;
357 return F;
358 }
359
360 template <typename Coord>
361 __device__ Coord CW_SlopeForce(Coord c, Coord n, Coord e, Coord s, Coord w, float GRAVITY)
362 {
363 float h = max(c.x, 0.0f);
364
365 Coord H;
366 H.x = 0.0f;
367 H.y = -GRAVITY * h * (e.w - w.w);
368 H.z = -GRAVITY * h * (s.w - n.w);
369 H.w = 0.0f;
370 return H;
371 }
372
373 template <typename Coord4D>
374 __global__ void CW_OneWaveStep(
375 DArray2D<Coord4D> grid_next,
376 DArray2D<Coord4D> grid,
377 int width,
378 int height,
379 float GRAVITY,
380 float timestep)
381 {
382 int x = threadIdx.x + blockIdx.x * blockDim.x;
383 int y = threadIdx.y + blockIdx.y * blockDim.y;
384
385 if (x < width && y < height)
386 {
387 int gridx = x + 1;
388 int gridy = y + 1;
389
390 Coord4D center = grid(gridx, gridy);
391
392 Coord4D north = grid(gridx, gridy - 1);
393
394 Coord4D west = grid(gridx - 1, gridy);
395
396 Coord4D south = grid(gridx, gridy + 1);
397
398 Coord4D east = grid(gridx + 1, gridy);
399
400 CW_FixShore(west, center, east);
401 CW_FixShore(north, center, south);
402
403 Coord4D u_south = 0.5f * (south + center) - timestep * (CW_VerticalPotential(south, GRAVITY) - CW_VerticalPotential(center, GRAVITY));
404 Coord4D u_north = 0.5f * (north + center) - timestep * (CW_VerticalPotential(center, GRAVITY) - CW_VerticalPotential(north, GRAVITY));
405 Coord4D u_west = 0.5f * (west + center) - timestep * (CW_HorizontalPotential(center, GRAVITY) - CW_HorizontalPotential(west, GRAVITY));
406 Coord4D u_east = 0.5f * (east + center) - timestep * (CW_HorizontalPotential(east, GRAVITY) - CW_HorizontalPotential(center, GRAVITY));
407
408 Coord4D u_center = center + timestep * CW_SlopeForce(center, north, east, south, west, GRAVITY) - timestep * (CW_HorizontalPotential(u_east, GRAVITY) - CW_HorizontalPotential(u_west, GRAVITY)) - timestep * (CW_VerticalPotential(u_south, GRAVITY) - CW_VerticalPotential(u_north, GRAVITY));
409 u_center.x = max(0.0f, u_center.x);
410
411 grid_next(gridx, gridy) = u_center;
412 }
413 }
414
415 template <typename Coord>
416 __global__ void CW_InitHeights(
417 DArray2D<Coord> height,
418 DArray2D<Coord> grid,
419 int patchSize,
420 float realSize)
421 {
422 int i = threadIdx.x + blockIdx.x * blockDim.x;
423 int j = threadIdx.y + blockIdx.y * blockDim.y;
424 if (i < patchSize && j < patchSize)
425 {
426 int gridx = i + 1;
427 int gridy = j + 1;
428
429 Coord gp = grid(gridx, gridy);
430 height(i, j) = gp;
431 }
432 }
433
434 template <typename Coord4D>
435 __global__ void CW_InitHeightGrad(
436 DArray2D<Coord4D> height,
437 int patchSize)
438 {
439 int i = threadIdx.x + blockIdx.x * blockDim.x;
440 int j = threadIdx.y + blockIdx.y * blockDim.y;
441 if (i < patchSize && j < patchSize)
442 {
443 int i_minus_one = (i - 1 + patchSize) % patchSize;
444 int i_plus_one = (i + 1) % patchSize;
445 int j_minus_one = (j - 1 + patchSize) % patchSize;
446 int j_plus_one = (j + 1) % patchSize;
447
448 Coord4D Dx = (height(i_plus_one, j) - height(i_minus_one, j)) / 2;
449 Coord4D Dz = (height(i, j_plus_one) - height(i, j_minus_one)) / 2;
450
451 height(i, j).z = Dx.y;
452 height(i, j).w = Dz.y;
453 }
454 }
455
456 template <typename Real, typename Coord3D, typename Coord4D>
457 __global__ void CW_InitHeightDisp(
458 DArray2D<Coord4D> heights,
459 DArray2D<Coord3D> displacement,
460 DArray2D<Coord4D> grid,
461 Real horizon)
462 {
463 unsigned int i = blockIdx.x * blockDim.x + threadIdx.x;
464 unsigned int j = blockIdx.y * blockDim.y + threadIdx.y;
465 if (i < displacement.nx() && j < displacement.ny())
466 {
467 int gridx = i + 1;
468 int gridy = j + 1;
469
470 Coord4D gij = grid(gridx, gridy);
471
472 displacement(i, j).x = 0;
473 displacement(i, j).y = gij.x + gij.w;
474 displacement(i, j).z = 0;
475
476 heights(i, j) = gij;
477 }
478 }
479
480 DEFINE_CLASS(CapillaryWave);
481}
Generated by doxygen 1.13.2