PeriDyno 1.0.0
Loading...
Searching...
No Matches
DualParticleIsphModule.cu
Go to the documentation of this file.
1#include <cuda_runtime.h>
2#include "DualParticleIsphModule.h"
3#include "Node.h"
4#include "Field.h"
5#include "ParticleSystem/Module/SummationDensity.h"
6#include "Collision/Attribute.h"
7#include "ParticleSystem/Module/Kernel.h"
8#include "Algorithm/Function2Pt.h"
9#include "Algorithm/CudaRand.h"
10
11
12namespace dyno
13{
14 IMPLEMENT_TCLASS(DualParticleIsphModule, TDataType)
15
16 template<typename TDataType>
17 DualParticleIsphModule<TDataType>::DualParticleIsphModule()
18 : ConstraintModule()
19 , m_airPressure(Real(0))
20 , m_reduce(NULL)
21 , m_arithmetic(NULL)
22 {
23 this->inParticleAttribute()->tagOptional(true);
24 this->inBoundaryNorm()->tagOptional(true);
25
26 this->varSamplingDistance()->setValue(Real(0.005));
27 this->varRestDensity()->setValue(Real(1000));
28 this->varSmoothingLength()->setValue(Real( 0.0125));
29 this->varPpeSmoothingLength()->setValue(Real( 0.0125));
30
31 m_summation = std::make_shared<SummationDensity<TDataType>>();
32 this->varRestDensity()->connect(m_summation->varRestDensity());
33 this->varSmoothingLength()->connect(m_summation->inSmoothingLength());
34 this->varSamplingDistance()->connect(m_summation->inSamplingDistance());
35 this->inRPosition()->connect(m_summation->inPosition());
36 this->inNeighborIds()->connect(m_summation->inNeighborIds());
37
38 m_vr_summation = std::make_shared<SummationDensity<TDataType>>();
39 this->varRestDensity()->connect(m_vr_summation->varRestDensity());
40 this->varSmoothingLength()->connect(m_vr_summation->inSmoothingLength());
41 this->varSamplingDistance()->connect(m_vr_summation->inSamplingDistance());
42 this->inVPosition()->connect(m_vr_summation->inPosition());
43 this->inRPosition()->connect(m_vr_summation->inOther());
44 this->inVRNeighborIds()->connect(m_vr_summation->inNeighborIds());
45
46 m_vv_summation = std::make_shared<SummationDensity<TDataType>>();
47 this->varRestDensity()->connect(m_vv_summation->varRestDensity());
48 this->varSmoothingLength()->connect(m_vv_summation->inSmoothingLength());
49 this->varSamplingDistance()->connect(m_vv_summation->inSamplingDistance());
50 this->inVPosition()->connect(m_vv_summation->inPosition());
51 this->inVVNeighborIds()->connect(m_vv_summation->inNeighborIds());
52
53 }
54
55 template<typename TDataType>
56 DualParticleIsphModule<TDataType>::~DualParticleIsphModule()
57 {
58
59 m_r.clear();
60 m_Aii.clear();
61 m_virtualAirFlag.clear();
62 m_virtualAirWeight.clear();
63 m_pressure.clear();
64 m_solidVirtualPaticleFlag.clear();
65 m_virtualVelocity.clear();
66 m_source.clear();
67 m_Ax.clear();
68 m_r.clear();
69 m_p.clear();
70 m_residual.clear();
71 m_Gp.clear();
72 m_GpNearSolid.clear();
73
74
75 if (m_reduce)
76 {
77 delete m_reduce;
78 }
79 if (m_arithmetic)
80 {
81 delete m_arithmetic;
82 }
83
84 }
85
86 /*
87 * MatrixDotVector
88 */
89 template <typename Vec4, typename Matrix4x4>
90 __device__ inline Vec4 MatrixDotVector(const Matrix4x4 M, const Vec4 V)
91 {
92 return Vec4(
93 V[0] * M(0, 0) + V[1] * M(0, 1) + V[2] * M(0, 2) + V[3] * M(0, 3),
94 V[0] * M(1, 0) + V[1] * M(1, 1) + V[2] * M(1, 2) + V[3] * M(1, 3),
95 V[0] * M(2, 0) + V[1] * M(2, 1) + V[2] * M(2, 2) + V[3] * M(2, 3),
96 V[0] * M(3, 0) + V[1] * M(3, 1) + V[2] * M(3, 2) + V[3] * M(3, 3)
97 );
98 }
99
100 /*
101 * Name : UpdateVelocity;
102 * Function : Update velocities with the pressure gradient field;
103 * Formular : vi = vi + dt * Gp / rho;
104 */
105 template <typename Real, typename Coord>
106 __global__ void DualParticle_UpdateVelocity(
107 DArray<Coord> gradientPress,
108 DArray<Coord> velocity,
109 DArray<Real> density,
110 Real mass,
111 Real dt
112 )
113 {
114 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
115 if (pId >= velocity.size()) return;
116
117 Coord temp = dt * gradientPress[pId] / density[pId];
118 velocity[pId] = velocity[pId] - temp;
119 }
120
121
122 template<typename TDataType>
123 bool DualParticleIsphModule<TDataType>::virtualArraysResize()
124 {
125 int num = this->inRPosition()->size();
126 int num_v = this->inVPosition()->size();
127
128 if (m_pressure.size() != num_v)
129 m_pressure.resize(num_v);
130
131 if (m_Ax.size() != num_v)
132 m_Ax.resize(num_v);
133 if (m_Aii.size() != num_v)
134 m_Aii.resize(num_v);
135 if (m_r.size() != num_v)
136 m_r.resize(num_v);
137 if (m_p.size() != num_v)
138 m_p.resize(num_v);
139 if (m_virtualAirFlag.size() != num_v)
140 m_virtualAirFlag.resize(num_v);
141 if (m_source.size() != num_v)
142 m_source.resize(num_v);
143 if (m_virtualVelocity.size() != num_v)
144 m_virtualVelocity.resize(num_v);
145 if (m_solidVirtualPaticleFlag.size() != num_v)
146 m_solidVirtualPaticleFlag.resize(num_v);
147 if (this->outVirtualBool()->isEmpty())
148 this->outVirtualBool()->allocate();
149 if (this->outVirtualWeight()->isEmpty())
150 this->outVirtualWeight()->allocate();
151 if (this->outVirtualBool()->size() != num_v)
152 this->outVirtualBool()->resize(num_v);
153 if (this->outVirtualWeight()->size() != num_v)
154 this->outVirtualWeight()->resize(num_v);
155
156 if (m_reduce)
157 {
158 delete m_reduce;
159 m_reduce = Reduction<float>::Create(num_v);
160 }
161 else
162 {
163 m_reduce = Reduction<float>::Create(num_v);
164 }
165
166 if (m_arithmetic)
167 {
168 delete m_arithmetic;
169 m_arithmetic = Arithmetic<float>::Create(num_v);
170 }
171 else
172 {
173 m_arithmetic = Arithmetic<float>::Create(num_v);
174 }
175 return true;
176 }
177
178 template<typename TDataType>
179 bool DualParticleIsphModule<TDataType>::realArraysResize()
180 {
181 int num = this->inRPosition()->size();
182
183 if (m_Gp.size() != num)
184 m_Gp.resize(num);
185
186 if (m_GpNearSolid.size() != num)
187 m_GpNearSolid.resize(num);
188
189 //if (m_arithmetic_r)
190 //{
191 // delete m_arithmetic_r;
192 // m_arithmetic_r = Arithmetic<float>::Create(num);
193 //}
194 //else
195 //{
196 // m_arithmetic_r = Arithmetic<float>::Create(num);
197 //}
198 return true;
199 }
200
201
202
203 /*
204 * Name : DualParticle_SolidVirtualParticleDetect
205 * Function : A virtual particle is Neumann boundary particle or not?
206 */
207 template <typename Real, typename Coord>
208 __global__ void DualParticle_VirtualSolidParticleDetection(
209 DArray<bool> solidVirtual,
210 DArray<Real> fluidFraction,
211 DArray<Coord> virtualPosition,
212 DArray<Coord> realPosition,
213 DArray<Attribute> attribute,
214 DArray<Real> density,
215 DArrayList<int> vr_neighbors,
216 DArray<bool> virtualAir,
217 CubicKernel<Real> kernel,
218 Real mass,
219 Real threshold,
220 Real h
221 )
222 {
223 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
224 if (pId >= solidVirtual.size()) return;
225
226 solidVirtual[pId] = false;
227
228 List<int> & list_i_vr = vr_neighbors[pId];
229 int nbSize = list_i_vr.size();
230 Real c = 0.0f;
231 Real w = 0.0f;
232 for (int ne = 0; ne < nbSize; ne++)
233 {
234 int j = list_i_vr[ne];
235 Real r_ij = (virtualPosition[pId] - realPosition[j]).norm();
236 if (!attribute[j].isFluid())
237 {
238 c += kernel.Weight(r_ij, h) * mass / density[j];
239 }
240 w += kernel.Weight(r_ij, h) * mass / density[j];
241 //w = 1.0f;
242 }
243 if (w < EPSILON) w = EPSILON;
244 fluidFraction[pId] = c / w;
245
246
247
248 if (fluidFraction[pId] > threshold)
249 {
250 solidVirtual[pId] = true;
251 }
252 else
253 {
254 solidVirtual[pId] = false;
255
256 }
257
258
259
260 }
261
262
263 template <typename Real, typename Coord>
264 __global__ void DualParticle_SolidBoundaryDivergenceCompsate(
265 DArray<Real> source,
266 DArray<bool> solidVirtual,
267 DArray<bool> airVirtual,
268 DArray<Coord> virtualVelocity,
269 DArray<Coord> realVelocity,
270 DArray<Coord> virtualPosition,
271 DArray<Coord> realPosition,
272 DArray<Attribute> attribute,
273 DArray<Coord> norm,
274 DArray<Real> density,
275 DArrayList<int> vr_neighbors,
276 CubicKernel<Real> kernel,
277 Real restDensity,
278 Real mass,
279 Real dt,
280 Real h
281 )
282 {
283 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
284 if (pId >= virtualPosition.size()) return;
285 if (solidVirtual[pId] == true) return;
286 if (airVirtual[pId] == true) return;
287
288 List<int> & list_i_vr = vr_neighbors[pId];
289 int nbSize_vr = list_i_vr.size();
290
291 Real value(0);
292 for (int ne = 0; ne < nbSize_vr; ne++)
293 {
294 int j = list_i_vr[ne];
295 if (!attribute[j].isFluid())
296 {
297
298 Real r_ij = (virtualPosition[pId] - realPosition[j]).norm();
299 Coord Gwij(0);
300
301 if (r_ij > EPSILON && !attribute[j].isFluid())
302 {
303
304 Coord dVel = virtualVelocity[pId] - realVelocity[j];
305 Real magNVel = dVel.dot(norm[j]);
306 Coord nVel = magNVel * norm[j];
307 Coord tVel = dVel - nVel;
308
309
310 Gwij = kernel.Gradient(r_ij, h) * (virtualPosition[pId] - realPosition[j]) / r_ij;
311
312
313 if (magNVel < -EPSILON)
314 {
315 value += 2 * (nVel + 0.01*tVel).dot(Gwij) * mass / density[pId];
316 }
317 else
318 {
319 value += 2 * (0.1 * nVel + 0.01*tVel).dot(Gwij) * mass / density[pId];
320 }
321
322 }
323
324
325 }
326 }
327
328 source[pId] -= -value * restDensity / dt;
329
330 }
331
332 /*
333 * @briedf Use CSPH method to calculate the virtual particle velocity;
334 * Virtual-air particle velocity should be zero;
335 *
336 */
337
338 template <typename Real, typename Coord>
339 __global__ void DualParticle_SmoothVirtualVelocity(
340 DArray<Coord> virtualVelocity,
341 DArray<Coord> velocity,
342 DArray<Coord> virtualPosition,
343 DArray<Coord> realPosition,
344 DArrayList<int> vr_neighbors,
345 DArray<bool> virtualAir,
346 DArray<bool> virtualSolid,
347 DArray<Real> density,
348 CubicKernel<Real> kernel,
349 Real mass,
350 Real h
351 )
352 {
353 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
354 if (pId >= virtualPosition.size()) return;
355 if (virtualSolid[pId] == true) return;
356 if (virtualAir[pId] == true)
357 {
358 virtualVelocity[pId] = Coord(0.0f);
359 return;
360 }
361
362 List<int> & list_i_vr = vr_neighbors[pId];
363 int nbSize_vr = list_i_vr.size();
364 Real total_w(0);
365 Coord total_vw(0);
366 for (int ne = 0; ne < nbSize_vr; ne++)
367 {
368 int j = list_i_vr[ne];
369 Real r_ij = (virtualPosition[pId] - realPosition[j]).norm();
370 Real wij = kernel.Weight(r_ij ,h);
371 total_w += wij;
372 total_vw += velocity[j] * wij;
373 }
374
375 if (total_w > EPSILON)
376 {
377 if (nbSize_vr == 0) printf("1-ERROR: pId %d \r\n", pId);
378 virtualVelocity[pId] = total_vw / total_w;
379 }
380 else
381 {
382 if (nbSize_vr != 0) printf("0-ERROR: pId %d \r\n", pId);
383 virtualVelocity[pId] = Coord(0.0f);
384 }
385
386 }
387
388 template <typename Real, typename Coord>
389 __global__ void DualParticle_SourceTerm(
390 DArray<Real> source,
391 DArray<Coord> virtualVelocity,
392 DArray<Coord> velocity,
393 DArray<Coord> virtualPosition,
394 DArray<Coord> realPosition,
395 DArray<Real> density,
396 DArray<Real> vdensity,
397 DArrayList<int> vr_neighbors,
398 DArray<Attribute> attribute,
399 DArray<bool> virtualAir,
400 DArray<bool> virtualSolid,
401 DArray<Coord> boundaryNorm,
402 CubicKernel<Real> kernel,
403 Real restDensity,
404 Real mass,
405 Real dt,
406 Real h
407 )
408 {
409 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
410 if (pId >= source.size()) return;
411 if (virtualSolid[pId] == true)
412 {
413 source[pId] = 0.0f;
414 return;
415 }
416
417 if (virtualAir[pId] == true)
418 {
419 source[pId] = 0.0f;
420 }
421 else
422 {
423 List<int> & list_i_vr = vr_neighbors[pId];
424 int nbSize_vr = list_i_vr.size();
425
426 Real value(0);
427 for (int ne = 0; ne < nbSize_vr; ne++)
428 {
429 int j = list_i_vr[ne];
430 Real r_ij = (virtualPosition[pId] - realPosition[j]).norm();
431 Coord Gwij(0);
432 if(r_ij > EPSILON & attribute[j].isFluid())
433 {
434 Gwij = kernel.Gradient(r_ij, h) * (virtualPosition[pId] - realPosition[j]) / r_ij;
435 }
436 value += (velocity[j] - virtualVelocity[pId]).dot(Gwij) * mass / vdensity[pId];
437 }
438 source[pId] = -value * restDensity / dt;
439 }
440 }
441
442 /*
443 * Name: VirtualAirParticleDetection
444 * Return: Virtual-air particle flag
445 * Function: rho_v(virtual particles'real density) < threashold ? true: false
446 */
447 template <typename Real>
448 __global__ void DualParticle_VirtualAirParticleDetection
449 (
450 DArray<bool> m_virtualAirFlag,
451 DArrayList<int> vr_neighbors,
452 DArray<Real> rho_v,
453 Real threshold
454 )
455 {
456 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
457 if (pId >= m_virtualAirFlag.size()) return;
458
459 List<int> & list_i_vr = vr_neighbors[pId];
460 int nbSize_vr = list_i_vr.size();
461
462
463 if ((rho_v[pId] > threshold) && (nbSize_vr > 0))
464 {
465 m_virtualAirFlag[pId] = false;
466
467 }
468 else
469 {
470 m_virtualAirFlag[pId] = true;
471 }
472 }
473
474
475
476
477 template <typename Real, typename Coord>
478 __global__ void DualParticle_LaplacianPressure
479 (
480 DArray<Real> Ax,
481 DArray<Real> Aii,
482 DArray<Real> pressure,
483 DArray<Coord> v_pos,
484 DArray<bool> virtualAir,
485 DArray<bool> virtualSolid,
486 DArrayList<int> vv_neighbors,
487 DArray<Real> rho_vv,
488 CubicKernel<Real> kernel,
489 Real mass,
490 Real h
491 )
492 {
493 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
494 if (pId >= Ax.size()) return;
495 if (virtualSolid[pId] == true)
496 {
497 return;
498 }
499
500
501 if (virtualAir[pId] == true)
502 {
503 Ax[pId] = 0.0f;
504 pressure[pId] = 0.0f;
505 return;
506 }
507 List<int> & list_i = vv_neighbors[pId];
508 int nbSize_i = list_i.size();
509
510 Real temp = 0;
511 for (int ne = 0; ne < nbSize_i; ne++)
512 {
513 int j = list_i[ne];
514
515 Real rij = (v_pos[pId] - v_pos[j]).norm();
516
517
518
519 if (rij > EPSILON && virtualAir[j] == false && virtualSolid[j] == false)
520 //if (rij > EPSILON && virtualAir[j] == false)
521 {
522 Real dwij = kernel.Gradient(rij, h) / rij;
523 temp += 8 * mass * dwij / pow((rho_vv[pId] + rho_vv[j]), 2) * pressure[j];
524 }
525
526 }
527 Ax[pId] = Aii[pId] * pressure[pId] + temp;
528 }
529
530
531
532
533 template <typename Real, typename Coord>
534 __global__ void DualParticle_AiiInLaplacian
535 (
536 DArray<Real> Aii,
537 DArray<Coord> v_pos,
538 DArray<bool> virtualAir,
539 DArray<bool> virtualSolid,
540 DArrayList<int > vv_neighbors,
541 DArray<Real> rho_vv,
542 CubicKernel<Real> kernel,
543 Real mass,
544 Real h
545 )
546 {
547 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
548 if (pId >= Aii.size()) return;
549 if (virtualSolid[pId] == true) return;
550 if (virtualAir[pId] == true) {
551 Aii[pId] = EPSILON;
552 return;
553 }
554 List<int> & list_i = vv_neighbors[pId];
555 int nbSize_i = list_i.size();
556
557 Real temp = 0;
558
559 for (int ne = 0; ne < nbSize_i; ne++)
560 {
561 int j = list_i[ne];
562
563 Real rij = (v_pos[pId] - v_pos[j]).norm();
564
565 //if (rij > EPSILON && virtualAir[j] == false && virtualSolid[j] == false)
566 if (rij > EPSILON && virtualAir[j] == false)
567 {
568 Real dwij = kernel.Gradient(rij, h) / rij;
569 temp += 8 * mass * dwij / pow((rho_vv[pId] + rho_vv[j]), 2);
570 }
571 }
572 Aii[pId] = -temp;
573
574 }
575
576 /*
577 * @brief: Using SPH aproach to calculate Gradient pressures of fluid particles.
578 */
579 template <typename Real, typename Coord>
580 __global__ void DualParticle_GradientPressure
581 (
582 DArray<Coord> gradient,
583 DArray<Real> pressure,
584 DArray<Coord> velocity,
585 DArray<Coord> v_pos,
586 DArray<Coord> r_pos,
587 DArray<Attribute> attribute,
588 DArray<bool> virtualAir,
589 DArray<bool> virtualSolid,
590 DArrayList<int> vr_neighbors,
591 DArrayList<int> rv_neighbors,
592 DArray<Real> rho_r,
593 DArray<Real> rho_v,
594 CubicKernel<Real> kernel,
595 Real rho_0,
596 Real dt,
597 Real mass,
598 Real h
599 )
600 {
601 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
602
603 if (pId >= gradient.size()) return;
604 if (!attribute[pId].isFluid()) return;
605
606 Coord value(0);
607 List<int> & list_i_rv = rv_neighbors[pId];
608 int nbSize_i = list_i_rv.size();
609
610 for (int ne = 0; ne < nbSize_i; ne++)
611 {
612 int j = list_i_rv[ne];
613 Real r_ij = (r_pos[pId] - v_pos[j]).norm();
614 if ((r_ij > EPSILON) && (virtualAir[j] != true) && (virtualSolid[j] != true))
615 {
616 value += (mass / rho_0) * pressure[j] * (r_pos[pId] - v_pos[j]) * kernel.Gradient(r_ij, h) / r_ij;
617 }
618 }
619 value = value * dt / rho_0;
620 gradient[pId] = value;
621 //Coord v = velocity[pId];
622 velocity[pId] -= gradient[pId] ;
623
624 }
625
626
627 /*
628 * Gradient correction Martrix
629 * @Paper: Fang et al 2009.
630 * @Paper: Improved SPH methods for simulating free surface flows of viscous fluids
631 * @Paper: Applied Numerical Mathematics 59 (2009) 251¨C271
632 */
633
634 template <typename Real, typename Coord>
635 __global__ void DualParticle_ImprovedGradient(
636 DArray<Coord> improvedGradient,
637 DArrayList<int> rv_neighbors,
638 DArray<Coord> vpos,
639 DArray<Coord> rpos,
640 DArray<Real> vrho,
641 DArray<bool> virtualAir,
642 DArray<bool> virtualSolid,
643 CubicKernel<Real> kernel,
644 Real mass,
645 Real rho_0,
646 Real h
647 )
648 {
649 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
650 if (pId >= improvedGradient.size()) return;
651
652 Real V = mass / rho_0;
653 glm::mat4x4 C(0);
654 List<int> & list_i_rr = rv_neighbors[pId];
655 int nbSize_i = list_i_rr.size();
656 for (int ne = 0; ne < nbSize_i; ne++)
657 {
658 int j = list_i_rr[ne];
659 Coord xij = rpos[pId] - vpos[j];
660 Real rij = xij.norm();
661 Coord nij = xij / rij;
662 Real weight = kernel.Weight(rij, h);
663 Real dweight = kernel.Gradient(rij, h);
664 C[0][0] += weight;
665 C[0][1] += xij[0] * weight;
666 C[0][2] += xij[1] * weight;
667 C[0][3] += xij[2] * weight;
668
669 C[1][0] += nij[0] * dweight;
670 C[1][1] += xij[0] * nij[0] * dweight;
671 C[1][2] += xij[1] * nij[0] * dweight;
672 C[1][3] += xij[2] * nij[0] * dweight;
673
674 C[2][0] += nij[1] * dweight;
675 C[2][1] += xij[0] * nij[1] * dweight;
676 C[2][2] += xij[1] * nij[1] * dweight;
677 C[2][3] += xij[2] * nij[1] * dweight;
678
679 C[3][0] += nij[2] * dweight;
680 C[3][1] += xij[0] * nij[2] * dweight;
681 C[3][2] += xij[1] * nij[2] * dweight;
682 C[3][3] += xij[2] * nij[2] * dweight;
683
684 }
685 C = V * C;
686 //glm::mat4x4 i;
687 C = glm::inverse(C);
688 if (pId == 132)
689 {
690 printf("GLM:INVERSE:");
691 for (int ci = 0; ci < 4; ci++)
692 {
693 for (int cj = 0; cj < 4; cj++)
694 {
695 printf("%f, ", C[ci][cj]);
696 }
697 printf(" | ");
698 }
699 printf("\r\n");
700 }
701
702 }
703
704 /*
705 * @brief: Modify gradient pressures if the fluid particle near solids.
706 * @Paper: Bridson. Fluid simulation for computer graphics, Second Edition. (Section 5.1, The Discrete Pressure Gradient)
707 */
708
709 template <typename Real, typename Coord>
710 __global__ void DualParticle_GradientNearSolid
711 (
712 DArray<Coord> gradidentComp,
713 DArray<Coord> velocity,
714 DArray<Coord> v_pos,
715 DArray<Coord> r_pos,
716 DArray<Attribute> attribute,
717 DArray<Coord> norm,
718 DArray<bool> virtualAir,
719 DArray<bool> virtualSolid,
720 DArrayList<int> rr_neighbors,
721 DArray<Real> rho_r,
722 DArray<Real> rho_v,
723 CubicKernel<Real> kernel,
724 Real rho_0,
725 Real dt,
726 Real mass,
727 Real h
728 )
729 {
730 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
731 if (pId >= gradidentComp.size()) return;
732
733 gradidentComp[pId] = Coord(0.0f);
734
735 List<int> & list_i_rr = rr_neighbors[pId];
736 int nbSize_i = list_i_rr.size();
737 Coord & vel_i = velocity[pId];
738 Coord dvij(0.0f);
739 for (int ne = 0; ne < nbSize_i; ne++)
740 {
741 int j = list_i_rr[ne];
742 Real r_ij = (r_pos[pId] - r_pos[j]).norm();
743 Real weight = kernel.Weight(r_ij, h);
744 if (!attribute[j].isFluid())
745 {
746
747 Coord nij = (r_pos[pId] - r_pos[j]);
748 if(nij.norm() > EPSILON)
749 {
750 nij = nij.normalize();
751 }
752 else
753 {
754 nij = Coord(1.0f, 0.0f, 0.0f);
755 }
756
757 Coord normal_j = norm[j];
758 Coord dVel = vel_i - velocity[j];
759 Real magNVel = dVel.dot(normal_j);
760 Coord nVel = magNVel * normal_j;
761 Coord tVel = dVel - nVel;
762 if (magNVel < -EPSILON)
763 {
764 dvij += nij.dot(nVel + 0.01 * dVel) * weight * nij;
765 }
766 else
767 {
768 dvij += nij.dot(0.1 * nVel + 0.01 * dVel) * weight * nij;
769 }
770
771 }
772 }
773
774 gradidentComp[pId] = 2 * dt * dvij / rho_0;
775 velocity[pId] -= gradidentComp[pId];
776 }
777
778 /*
779 * @brief: Semi-analytical Dirichlet boundary.
780 * @Paper: Nair and Tomar, Computers & Fluids, 2014, 10(102): 304-314, An improved free surface modeling for incompressible SPH.
781 */
782 template <typename Real>
783 __global__ void DualParticle_CorrectAii
784 (
785 DArray<Real> Aii,
786 Real max_Aii,
787 Real c
788 )
789 {
790 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
791 if (pId >= Aii.size()) return;
792 max_Aii = max_Aii * c;
793 Real temp = Aii[pId];
794 if (temp < max_Aii)
795 {
796 temp = max_Aii;
797 }
798 Aii[pId] = temp;
799 }
800
801 /*
802 * @brief: Modified density drifting error (Liu, et al. Tog 2024, Equation 13.)
803 * @Paper: Abbas Khayyer and Hitoshi Gotoh. J. Comput. Phys. 230, 8 (2011). Enhancement of stability and accuracy of the moving particle semi-implicit method
804 */
805 template <typename Real>
806 __global__ void DualParticle_DensityCompensate
807 (
808 DArray<Real> source,
809 DArray<Real> rho,
810 DArray<bool> virtualAir,
811 DArray<bool> virtualSolid,
812 Real rho_0,
813 Real dt,
814 Real h
815 )
816 {
817 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
818 if (pId >= source.size()) return;
819 if (virtualSolid[pId] == true)
820 {
821 return;
822 }
823
824 if (virtualAir[pId] == true)
825 {
826 return;
827 }
828
829 if (rho[pId] > rho_0)
830 {
831 source[pId] += 1000000.0f * (rho[pId] - rho_0) / rho_0;
832 }
833 }
834
835
836 /*
837 * @brief: Neumann Boundary in Matrix of pressrue Poisson eqtution.
838 * @Paper: Bridson. Fluid simulation for computer graphics, Second Edition. (Section 5.3, The Pressure Equations.)
839 */
840 template <typename Real, typename Coord>
841 __global__ void DualParticle_AiiNeumannCorrect
842 (
843 DArray<Real> Aii,
844 DArray<bool> virtualAir,
845 DArray<bool> virtualSolid,
846 DArrayList<int > vv_neighbors,
847 DArray<Real> rho_vv,
848 DArray<Coord> v_pos,
849 CubicKernel<Real> kernel,
850 Real max_Aii,
851 Real mass,
852 Real h
853 )
854 {
855 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
856 if (pId >= Aii.size()) return;
857 if (virtualSolid[pId] == true) return;
858 if (virtualAir[pId] == true) return;
859
860 List<int> & list_i = vv_neighbors[pId];
861 int nbSize_i = list_i.size();
862
863 Real temp = 0;
864 int solidCounter = 0;
865
866 for (int ne = 0; ne < nbSize_i; ne++)
867 {
868 int j = list_i[ne];
869
870 Real rij = (v_pos[pId] - v_pos[j]).norm();
871
872 //if (rij > EPSILON && virtualAir[j] == false && virtualSolid[j] == true)
873 if (rij > EPSILON && virtualSolid[j] == true)
874 {
875 Real dwij = kernel.Gradient(rij, h) / rij;
876 temp += 8 * mass * dwij / pow((rho_vv[pId] + rho_vv[j]), 2);
877 solidCounter++;
878 }
879 }
880
881 Real value = Aii[pId] + temp;
882 Aii[pId] = value;
883 }
884
885 template <typename Coord>
886 __global__ void DualParticle_SolidVelocityReset
887 (
888 DArray<Coord> velocity,
889 DArray<Attribute> attribute
890 )
891 {
892 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
893 if (pId >= velocity.size()) return;
894
895 if (!attribute[pId].isFluid())
896 {
897 velocity[pId] = Coord(0.0f);
898 }
899
900 }
901
902
903
904 template <typename Real>
905 __global__ void DualParticle_VolumeTest
906 (
907 DArray<Real> volume,
908 DArray<Real> density,
909 Real mass
910 )
911 {
912 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
913 if (pId >= volume.size()) return;
914
915 volume[pId] = mass / density[pId];
916 }
917
918
919 template <typename Real>
920 __global__ void DualParticle_VirtualVolumeTest
921 (
922 DArray<Real> volume,
923 DArray<Real> density,
924 DArray<bool> air,
925 Real mass
926 )
927 {
928 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
929 if (pId >= volume.size()) return;
930 if (air[pId] == false)
931 volume[pId] = mass / density[pId];
932 else
933 volume[pId] = 0.0f;
934 }
935
936
937 /*
938 * @brief: Using euleric finite divegence aproach to calculate Gradient pressures
939 */
940 template <typename Real, typename Coord>
941 __global__ void DualParticle_VirtualGradientPressureFD(
942 DArray<Coord> vGp,
943 DArray<Coord> vpos,
944 DArrayList<int> vv_neighbors,
945 DArray<Real> pressure,
946 DArray<bool> virtualAir,
947 DArray<bool> virtualSolid,
948 Real dx
949 )
950 {
951 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
952 if (pId >= vpos.size()) return;
953
954 List<int> & list_i = vv_neighbors[pId];
955 int nbSize_i = list_i.size();
956 int count = 0;
957 Coord gp(0.0f);
958 for (int ne = 0; ne < nbSize_i; ne++)
959 {
960 int j = list_i[ne];
961 Real rij = (vpos[pId] - vpos[j]).norm();
962 if ((rij < dx)&&(rij > EPSILON))
963 {
964 count++;
965 gp += (pressure[pId] - pressure[j]) * (vpos[pId] - vpos[j]) / rij;
966 }
967 }
968
969 vGp[pId] = gp;
970
971 }
972
973 template <typename Real, typename Coord>
974 __global__ void DualParticle_VirtualGradientPressureMeshless(
975 DArray<Coord> vGp,
976 DArray<Coord> vpos,
977 DArrayList<int> vv_neighbors,
978 DArray<Real> pressure,
979 DArray<bool> virtualAir,
980 DArray<bool> virtualSolid,
981 CubicKernel<Real> kernel,
982 Real mass,
983 Real rho_0,
984 Real h
985 ) {
986
987 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
988 if (pId >= vpos.size()) return;
989
990 if (virtualAir[pId] == true)
991 {
992 vGp[pId] = Coord(0.0f);
993 return;
994 }
995
996 List<int> & list_i = vv_neighbors[pId];
997 int nbSize_i = list_i.size();
998 Real dwij(0.0f);
999 Coord value(0.0f);
1000 for (int ne = 0; ne < nbSize_i; ne++)
1001 {
1002 int j = list_i[ne];
1003 Real rij = (vpos[pId] - vpos[j]).norm();
1004 if(rij > EPSILON)
1005 {
1006 dwij = kernel.Gradient(rij, h);
1007 value += mass * (pressure[j] - pressure[pId]) * dwij * (vpos[pId] - vpos[j]) / (rij * rho_0);
1008 //value += mass * (pressure[j]) * dwij * (vpos[pId] - vpos[j]) / (rij * rho_0);
1009 }
1010 }
1011
1012 vGp[pId] = value;
1013 }
1014
1015
1016
1017 /*
1018 * @brief: Using moving least square aproach to calculate Gradient pressures
1019 */
1020 template <typename Real, typename Coord, typename Vector4, typename Matrix4x4>
1021 __global__ void DualParticle_MlsGradientPressure(
1022 DArray<Coord> rGp,
1023 DArray<Coord> velocity,
1024 DArray<Coord> vGp,
1025 DArray<Coord> vpos,
1026 DArray<Coord> rpos,
1027 DArrayList<int> rv_neighbors,
1028 DArray<Real> pressure,
1029 CubicKernel<Real> kernel,
1030 Real rho_0,
1031 Real dt,
1032 Real h
1033 )
1034 {
1035 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
1036 if (pId >= velocity.size()) return;
1037
1038
1039 List<int> & list_i = rv_neighbors[pId];
1040 int nbSize_i = list_i.size();
1041
1042 Vector4 P(0.0f); /* Basis Function: P((xi - xj)/h) */
1043 Vector4 P0(0.0f); /* Basis Function: P(0) */
1044 Matrix4x4 M(0.0f); /* Momentum Matrix of MLS*/
1045 Real dx_Ji(0.0f); /* (xi - xj)/h*/
1046 Real dy_Ji(0.0f); /* (yi - yj)/h*/
1047 Real dz_Ji(0.0f); /* (zi - zj)/h*/
1048 Real wij(0.0f); /*Weight*/
1049 Real rij(0.0f);
1050
1051
1052 /*Momentum Matrix*/
1053 for (int ne = 0; ne < nbSize_i; ne++)
1054 {
1055 int j = list_i[ne];
1056 rij = (rpos[pId] - vpos[j]).norm();
1057 wij = kernel.Weight(rij, h);
1058
1059 dx_Ji = (vpos[j] - rpos[pId])[0] / h;
1060 dy_Ji = (vpos[j] - rpos[pId])[1] / h;
1061 dz_Ji = (vpos[j] - rpos[pId])[2] / h;
1062
1063 M(0, 0) += wij; M(0, 1) += wij * dx_Ji; M(0, 2) += wij * dy_Ji; M(0, 3) += wij * dz_Ji;
1064 M(1, 0) += wij * dx_Ji; M(1, 1) += wij * dx_Ji * dx_Ji; M(1, 2) += wij * dx_Ji * dy_Ji; M(1, 3) += wij * dx_Ji * dz_Ji;
1065 M(2, 0) += wij * dy_Ji; M(2, 1) += wij * dx_Ji * dy_Ji; M(2, 2) += wij * dy_Ji * dy_Ji; M(2, 3) += wij * dy_Ji * dz_Ji;
1066 M(3, 0) += wij * dz_Ji; M(3, 1) += wij * dx_Ji * dz_Ji; M(3, 2) += wij * dy_Ji * dz_Ji; M(3, 3) += wij * dz_Ji * dz_Ji;
1067 }
1068
1069 Matrix4x4 M_inv = M.inverse();
1070
1071
1072 Coord tgp(0.0f);
1073 for (int ne = 0; ne < nbSize_i; ne++)
1074 {
1075 int j = list_i[ne];
1076 rij = (rpos[pId] - vpos[j]).norm();
1077 wij = kernel.Weight(rij, h);
1078
1079 dx_Ji = (vpos[j] - rpos[pId])[0] / h;
1080 dy_Ji = (vpos[j] - rpos[pId])[1] / h;
1081 dz_Ji = (vpos[j] - rpos[pId])[2] / h;
1082
1083 P[0] = 1.0f;
1084 P[1] = dx_Ji;
1085 P[2] = dy_Ji;
1086 P[3] = dz_Ji;
1087
1088 P0[0] = 1.0f;
1089 P0[1] = 0.0f;
1090 P0[2] = 0.0f;
1091 P0[3] = 0.0f;
1092
1093 tgp += P0.dot(MatrixDotVector(M_inv, P)) * wij * vGp[j];
1094
1095 }
1096
1097 rGp[pId] = dt * tgp / rho_0;
1098
1099 velocity[pId] -= rGp[pId];
1100
1101 }
1102
1103
1104 __global__ void DualParticle_AttributeInit
1105 (
1106 DArray<Attribute> atts
1107 )
1108 {
1109 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
1110 if (pId >= atts.size()) return;
1111
1112 atts[pId].setFluid();
1113 atts[pId].setDynamic();
1114 }
1115
1116 template<typename TDataType>
1117 void DualParticleIsphModule<TDataType>::constrain()
1118 {
1119 int num = this->inRPosition()->size();
1120 int num_v = this->inVPosition()->size();
1121
1122 cudaDeviceSynchronize();
1123
1124 if (m_Gp.size() != num)
1125 {
1126 realArraysResize();
1127 }
1128 if (m_Ax.size() != num_v)
1129 {
1130 virtualArraysResize();
1131 }
1132
1133 auto & m_virtualSolidFlag = this->outVirtualBool()->getData();
1134 Real dt = this->inTimeStep()->getData();
1135 Real h1 = this->varSmoothingLength()->getData();
1136 Real h2 = this->varPpeSmoothingLength()->getData();
1137
1138
1139 if (this->inParticleAttribute()->isEmpty()
1140 || this->inParticleAttribute()->size() != num
1141 || this->inBoundaryNorm()->size() != num
1142 )
1143 {
1144 this->inParticleAttribute()->allocate();
1145 this->inParticleAttribute()->resize(num);
1146 cuExecute(num, DualParticle_AttributeInit, this->inParticleAttribute()->getData());
1147 this->inBoundaryNorm()->resize(num);
1148 this->inBoundaryNorm()->reset();
1149 }
1150
1151
1152 m_summation->update();
1153 m_vv_summation->update();
1154 m_vr_summation->update();
1155
1156 m_particleMass = m_summation->getParticleMass();
1157 m_v_particleMass = m_vv_summation->getParticleMass();
1158
1159 Real restDensity= this->varRestDensity()->getValue();
1160
1161 Real MaxVirtualDensity =
1162 m_reduce->maximum(
1163 m_vr_summation->outDensity()->getData().begin(),
1164 m_vr_summation->outDensity()->getData().size()
1165 );
1166
1167 cuExecute(num_v, DualParticle_VirtualAirParticleDetection,
1168 m_virtualAirFlag,
1169 this->inVRNeighborIds()->getData(),
1170 m_vr_summation->outDensity()->getData(),
1171 MaxVirtualDensity * 0.1f
1172 );
1173
1174 cuExecute(num_v, DualParticle_VirtualSolidParticleDetection,
1175 m_virtualSolidFlag,
1176 this->outVirtualWeight()->getData(),
1177 this->inVPosition()->getData(),
1178 this->inRPosition()->getData(),
1179 this->inParticleAttribute()->getData(),
1180 m_summation->outDensity()->getData(),
1181 this->inVRNeighborIds()->getData(),
1182 m_virtualAirFlag,
1183 kernel,
1184 m_particleMass,
1185 0.999f,
1186 h1
1187 );
1188
1189 cuExecute(num_v, DualParticle_SmoothVirtualVelocity,
1190 m_virtualVelocity,
1191 this->inVelocity()->getData(),
1192 this->inVPosition()->getData(),
1193 this->inRPosition()->getData(),
1194 this->inVRNeighborIds()->getData(),
1195 m_virtualAirFlag,
1196 m_virtualSolidFlag,
1197 m_summation->outDensity()->getData(),
1198 kernel,
1199 m_particleMass,
1200 h1
1201 );
1202
1203 m_source.reset();
1204
1205 cuExecute(num, DualParticle_SolidVelocityReset,
1206 this->inVelocity()->getData(),
1207 this->inParticleAttribute()->getData()
1208 );
1209
1210 cuExecute(num_v, DualParticle_SourceTerm,
1211 m_source,
1212 m_virtualVelocity,
1213 this->inVelocity()->getData(),
1214 this->inVPosition()->getData(),
1215 this->inRPosition()->getData(),
1216 m_summation->outDensity()->getData(),
1217 m_vr_summation->outDensity()->getData(),
1218 this->inVRNeighborIds()->getData(),
1219 this->inParticleAttribute()->getData(),
1220 m_virtualAirFlag,
1221 m_virtualSolidFlag,
1222 this->inBoundaryNorm()->getData(),
1223 kernel,
1224 restDensity,
1225 m_particleMass,
1226 dt,
1227 h1
1228 );
1229
1230 cuExecute(num_v, DualParticle_SolidBoundaryDivergenceCompsate,
1231 m_source,
1232 m_virtualSolidFlag,
1233 m_virtualAirFlag,
1234 m_virtualVelocity,
1235 this->inVelocity()->getData(),
1236 this->inVPosition()->getData(),
1237 this->inRPosition()->getData(),
1238 this->inParticleAttribute()->getData(),
1239 this->inBoundaryNorm()->getData(),
1240 m_vr_summation->outDensity()->getData(),
1241 this->inVRNeighborIds()->getData(),
1242 kernel,
1243 restDensity,
1244 m_v_particleMass,
1245 dt,
1246 h1
1247 );
1248
1249 cuExecute(num_v, DualParticle_DensityCompensate,
1250 m_source,
1251 m_vr_summation->outDensity()->getData(),
1252 m_virtualAirFlag,
1253 m_virtualSolidFlag,
1254 restDensity,
1255 dt,
1256 h1
1257 );
1258
1259 m_r.reset();
1260 m_Ax.reset();
1261
1262 if (m_pressure.size() != virtualNumber_old)
1263 {
1264 m_pressure.reset();
1265 }
1266
1267 virtualNumber_old = m_pressure.size();
1268
1269 cuExecute(num_v, DualParticle_AiiInLaplacian,
1270 m_Aii,
1271 this->inVPosition()->getData(),
1272 m_virtualAirFlag,
1273 m_virtualSolidFlag,
1274 this->inVVNeighborIds()->getData(),
1275 m_vv_summation->outDensity()->getData(),
1276 kernel,
1277 m_v_particleMass,
1278 h2
1279 );
1280
1281 if ((frag_number <= 3) || abs(max_Aii) < EPSILON)
1282 {
1283 //if(m_reduce->maximum(m_Aii.begin(), m_Aii.size()) > 0)
1284 max_Aii = m_reduce->maximum(m_Aii.begin(), m_Aii.size());
1285 }
1286
1287 cuExecute(num_v, DualParticle_CorrectAii,
1288 m_Aii,
1289 max_Aii,
1290 1.0f
1291 );
1292
1293 cuExecute(num_v, DualParticle_AiiNeumannCorrect,
1294 m_Aii,
1295 m_virtualAirFlag,
1296 m_virtualSolidFlag,
1297 this->inVVNeighborIds()->getData(),
1298 m_vv_summation->outDensity()->getData(),
1299 this->inVPosition()->getData(),
1300 kernel,
1301 max_Aii,
1302 m_v_particleMass,
1303 h2
1304 );
1305
1306 cuExecute(num_v, DualParticle_LaplacianPressure,
1307 m_Ax,
1308 m_Aii,
1309 m_pressure,
1310 this->inVPosition()->getData(),
1311 m_virtualAirFlag,
1312 m_virtualSolidFlag,
1313 this->inVVNeighborIds()->getData(),
1314 m_vv_summation->outDensity()->getData(),
1315 kernel,
1316 m_v_particleMass,
1317 h2
1318 );
1319
1320 Function2Pt::subtract(m_r, m_source, m_Ax);
1321 m_p.assign(m_r);
1322
1323 Real rr = m_arithmetic->Dot(m_r, m_r);
1324 Real err = m_r.size() > 0 ? sqrt(rr / m_r.size()) : 0.0f;
1325 Real max_err = err;
1326 if (abs(max_err) < EPSILON) max_err = EPSILON;
1327 unsigned int iter = 0;
1328 Real threshold = this->varResidualThreshold()->getValue();
1329 while ((iter < 500) && (err / max_err > threshold) && (err > threshold))
1330 {
1331 iter++;
1332
1333 m_Ax.reset();
1334
1335 cuExecute(num_v, DualParticle_LaplacianPressure,
1336 m_Ax,
1337 m_Aii,
1338 m_p,
1339 this->inVPosition()->getData(),
1340 m_virtualAirFlag,
1341 m_virtualSolidFlag,
1342 this->inVVNeighborIds()->getData(),
1343 m_vv_summation->outDensity()->getData(),
1344 kernel,
1345 m_v_particleMass,
1346 h2
1347 );
1348
1349 float alpha = rr / m_arithmetic->Dot(m_p, m_Ax);
1350 Function2Pt::saxpy(m_pressure, m_p, m_pressure, alpha);
1351 Function2Pt::saxpy(m_r, m_Ax, m_r, -alpha);
1352
1353 Real rr_old = rr;
1354
1355 rr = m_arithmetic->Dot(m_r, m_r);
1356
1357 Real beta = rr / rr_old;
1358 Function2Pt::saxpy(m_p, m_p, m_r, beta);
1359 err = sqrt(rr / m_r.size());
1360 //std::cout<<"*DUAL-ISPH:: iter:"<< iter <<": Err-" << err << std::endl;
1361 }
1362 std::cout << "*DUAL-ISPH::Solver::Iteration:" << iter << "||RelativeError:" << err/ max_err * 100 <<"%" << std::endl;
1363
1364 m_GpNearSolid.reset();
1365
1366 cuExecute(num, DualParticle_GradientPressure,
1367 m_Gp,
1368 m_pressure,
1369 this->inVelocity()->getData(),
1370 this->inVPosition()->getData(),
1371 this->inRPosition()->getData(),
1372 this->inParticleAttribute()->getData(),
1373 m_virtualAirFlag,
1374 m_virtualSolidFlag,
1375 this->inVRNeighborIds()->getData(),
1376 this->inRVNeighborIds()->getData(),
1377 m_summation->outDensity()->getData(),
1378 m_vr_summation->outDensity()->getData(),
1379 kernel,
1380 restDensity,
1381 dt,
1382 m_v_particleMass,
1383 h1
1384 );
1385
1386 cuExecute(num, DualParticle_GradientNearSolid,
1387 m_GpNearSolid,
1388 this->inVelocity()->getData(),
1389 this->inVPosition()->getData(),
1390 this->inRPosition()->getData(),
1391 this->inParticleAttribute()->getData(),
1392 this->inBoundaryNorm()->getData(),
1393 m_virtualAirFlag,
1394 m_virtualSolidFlag,
1395 this->inNeighborIds()->getData(),
1396 m_summation->outDensity()->getData(),
1397 m_vr_summation->outDensity()->getData(),
1398 kernel,
1399 restDensity,
1400 dt,
1401 m_v_particleMass,
1402 h1
1403 );
1404
1405 /*
1406 * Gradient Pressure On Virtual Points (Finite Difference)
1407 */
1408 //DualParticle_VirtualGradientPressureFD << <pDims_r, BLOCK_SIZE >> > (
1409 // m_vGp,
1410 // this->inVPosition()->getData(), //DArray<Coord> vpos,
1411 // this->inVVNeighborIds()->getData(), //DArrayList<int> vv_neighbors,
1412 // m_pressure, //DArray<Real> pressure,
1413 // m_virtualAirFlag, //DArray<bool> virtualAir,
1414 // m_virtualSolidFlag, //DArray<bool> virtualSolid,
1415 // 0.0051f
1416 //);
1417
1418
1419 /*
1420 * Improved method for Gradient Pressure On Virtual Points
1421 */
1422 //DualParticle_ImprovedGradient << <pDims_r, BLOCK_SIZE >> > (
1423 // m_improvedGradient,
1424 // this->inRVNeighborIds()->getData(),
1425 // this->inVPosition()->getData(),
1426 // this->inRPosition()->getData(),
1427 // m_vv_summation->outDensity()->getData(),
1428 // m_virtualAirFlag,
1429 // m_virtualSolidFlag,
1430 // kernel,
1431 // m_v_particleMass,
1432 // 1000.0f,
1433 // h1
1434 // );
1435
1436
1437 /*
1438 * Gradient Pressure On Virtual Points (SPH)
1439 */
1440 //DualParticle_VirtualGradientPressureMeshless << <pDims_v, BLOCK_SIZE >> > (
1441 // m_vGp,
1442 // this->inVPosition()->getData(), //DArray<Coord> vpos,
1443 // this->inVVNeighborIds()->getData(), //DArrayList<int> vv_neighbors,
1444 // //pseudoPressure,
1445 // m_pressure, ////DArray<Real> pressure,
1446 // m_virtualAirFlag, //DArray<bool> virtualAir,
1447 // m_virtualSolidFlag, //DArray<bool> virtualSolid,
1448 // kernel,
1449 // m_v_particleMass,
1450 // 1000.0f,
1451 // 0.011f
1452 // );
1453
1454
1455 /*
1456 * Moving least square approgh for Gradient Pressure On Virtual Points (SPH)
1457 */
1458 //DualParticle_MlsGradientPressure <Real, Coord, Vector<Real, 4>, SquareMatrix<Real , 4>> << <pDims_r, BLOCK_SIZE >> > (
1459 // m_rGp, //DArray<Coord> rGp,
1460 // this->inVelocity()->getData(),
1461 // m_vGp, //DArray<Coord> vGp,
1462 // this->inVPosition()->getData(), //DArray<Coord> vpos,
1463 // this->inRPosition()->getData(), //DArray<Coord> rpos,
1464 // this->inRVNeighborIds()->getData(), //DArrayList<int> rv_neighbors,
1465 // m_pressure, //DArray<Real> pressure,
1466 // kernel,
1467 // 1000.0f,
1468 // dt,
1469 // 0.011 //Real h
1470 // );
1471
1472
1473
1474 }
1475
1476 template<typename TDataType>
1477 bool DualParticleIsphModule<TDataType>::initializeImpl()
1478 {
1479 cuSynchronize();
1480
1481 int num = this->inRPosition()->size();
1482 int num_v = this->inVPosition()->size();
1483
1484 m_reduce = Reduction<float>::Create(num_v);
1485 m_arithmetic = Arithmetic<float>::Create(num_v);
1486 // m_arithmetic_r = Arithmetic<float>::Create(num);
1487
1488 return true;
1489 }
1490
1491 DEFINE_CLASS(DualParticleIsphModule);
1492}
Generated by doxygen 1.13.2