PeriDyno 1.0.0
Loading...
Searching...
No Matches
SemiAnalyticalIncompressibilityModule.cu
Go to the documentation of this file.
1/**
2 * @author : Yue Chang (yuechang@pku.edu.cn)
3 * @date : 2021-08-06
4 * @description: Implemendation of SemiAnalyticalIncompressibilityModule class, implemendation of semi-analytical perojection-based fluids
5 * introduced in the paper <Semi-analytical Solid Boundary Conditions for Free Surface Flows>
6 * @version : 1.1
7 */
8#include <cuda_runtime.h>
9#include "SemiAnalyticalIncompressibilityModule.h"
10
11#include "ParticleSystem/Module/SummationDensity.h"
12#include "ParticleSystem/Module/Kernel.h"
13
14#include "Collision/Attribute.h"
15#include "IntersectionArea.h"
16#include "Algorithm/Function2Pt.h"
17
18#include "Node.h"
19
20namespace dyno
21{
22__device__ inline float kernWeight(const float r, const float h)
23{
24 const float q = r / h;
25 if (q > 1.0f)
26 return 0.0f;
27 else
28 {
29 const float d = 1.0f - q;
30 const float hh = h * h;
31 // return 45.0f / ((float)M_PI * hh*h) *d*d;
32 return (1.0 - pow(q, 4.0f));
33 // return (1.0 - q)*(1.0 - q)*h*h;
34 }
35}
36
37__device__ inline float kernWeightMesh(const float r, const float h)
38{
39 if (r > h)
40 return 0.0f;
41 return 4.0 / 21.0 * pow(h, 3.0f) - pow(r, 3.0f) / 3.0 + pow(r, 7.0f) / 7.0 / pow(h, 4.0f);
42 //G(r) in eq11 for weighted kernel function
43}
44
45__device__ inline float kernWR(const float r, const float h)
46{
47 float w = kernWeight(r, h);
48 const float q = r / h;
49 if (q < 0.4f)
50 {
51 return w / (0.4f * h);
52 }
53 return w / r;
54}
55
56__device__ inline float kernWRMesh(const float r, const float h)
57{
58 const float q = r / h;
59 const float h04 = 0.4f * h;
60 if (q < 0.4f)
61 return //(pow(h04, 3.0f) - pow(r, 3.0f)) / 3.0 - (pow(h, 7.0f) - pow(r, 7.0f)) / 7.0 / pow(h, 4.0f) + 0.28 * h * h;
62 h * h / 3.0f - (h04 * h04 / 2.0f - pow(h04, 6.0f) / (6.0f * pow(h, 4.0f)))
63 + ((pow(h04, 3.0f) - pow(r, 3.0f)) / 3.0f - (pow(h04, 7.0f) - pow(r, 7.0f)) / 7.0f / pow(h, 4.0f)) / h04;
64 else
65 return h * h / 3.0f - (r * r / 2.0f - pow(r, 6.0f) / (6.0f * pow(h, 4.0f))); //(h * h - r * r) / 3.0;
66 //G(r) in eq11 for weighted kernel function's gradient
67}
68
69__device__ inline float kernWRR(const float r, const float h)
70{
71 float w = kernWeight(r, h);
72 const float q = r / h;
73 if (q < 0.4f)
74 {
75 return w / (0.16f * h * h);
76 }
77 return w / r / r;
78}
79
80__device__ inline float kernWRRMesh(const float r, const float h)
81{
82 const float q = r / h;
83 const float h04 = 0.4f * h;
84 if (q < 0.4f)
85 return 0.6f * h - 1.0f / 5.0f / pow(h, 4.0f) * (pow(h, 5.0f) - pow(h04, 5.0f)) + 1.0 / (0.16 * h * h) * (1.0 / 3.0 * (pow(h04, 3.0f) - pow(r, 3.0f)) - 1.0 / (7.0 * pow(h, 4.0f)) * (pow(h04, 7.0f) - pow(r, 7.0f)));
86 else
87 return h - r - (1.0 / 5.0 / pow(h, 4.0f) * (pow(h, 5.0f) - pow(r, 5.0f)));
88 //G(r) in eq11 for weighted kernel function's second-order gradient
89}
90//
91// //Triangle Clustering
92// template <typename Coord>
93// __global__ void VC_Sort_Neighbors(
94// DArray<Coord> position,
95// DArray<TopologyModule::Triangle> m_triangle_index,
96// DArray<Coord> positionTri,
97// DArrayList<int> neighborsTri)
98// {
99// int pId = threadIdx.x + (blockIdx.x * blockDim.x);
100// if (pId >= position.size())
101// return;
102// int nbSizeTri = neighborsTri.getNeighborSize(pId);
103//
104// Coord pos_i = position[pId];
105//
106// for (int ne = nbSizeTri / 2 - 1; ne >= 0; ne--)
107// {
108// int start = ne;
109// int end = nbSizeTri - 1;
110// int c = start;
111// int l = 2 * c + 1;
112// int tmp = neighborsTri.getElement(pId, c);
113// for (; l <= end; c = l, l = 2 * l + 1)
114// {
115// if (l < end)
116// {
117// bool judge = false;
118// {
119// int idx1, idx2;
120// idx1 = neighborsTri.getElement(pId, l);
121// idx2 = neighborsTri.getElement(pId, l + 1);
122// if (neighborsTri.getElement(pId, l) < 0 && neighborsTri.getElement(pId, l + 1) < 0)
123// {
124// idx1 *= -1;
125// idx1--;
126// idx2 *= -1;
127// idx2--;
128// Triangle3D t3d1(positionTri[m_triangle_index[idx1][0]], positionTri[m_triangle_index[idx1][1]], positionTri[m_triangle_index[idx1][2]]);
129// Triangle3D t3d2(positionTri[m_triangle_index[idx2][0]], positionTri[m_triangle_index[idx2][1]], positionTri[m_triangle_index[idx2][2]]);
130//
131// Coord normal1 = t3d1.normal().normalize();
132// Coord normal2 = t3d2.normal().normalize();
133//
134// Point3D p3d(pos_i);
135//
136// Plane3D PL1(positionTri[m_triangle_index[idx1][0]], t3d1.normal());
137// Plane3D PL2(positionTri[m_triangle_index[idx2][0]], t3d2.normal());
138//
139// Real dis1 = p3d.distance(PL1);
140// Real dis2 = p3d.distance(PL2);
141//
142// if (abs(dis1 - dis2) < EPSILON && abs(normal1[0] - normal2[0]) < EPSILON && abs(normal1[1] - normal2[1]) < EPSILON)
143// judge = normal1[2] < normal2[2] ? true : false;
144// else if (abs(dis1 - dis2) < EPSILON && abs(normal1[0] - normal2[0]) < EPSILON)
145// judge = normal1[1] < normal2[1] ? true : false;
146// else if (abs(dis1 - dis2) < EPSILON)
147// judge = normal1[0] < normal2[0] ? true : false;
148// else
149// judge = dis1 <= dis2 ? true : false;
150// }
151// else
152// judge = neighborsTri.getElement(pId, l) < neighborsTri.getElement(pId, l + 1) ? true : false;
153// }
154// if (judge)
155// l++;
156// }
157// bool judge = false;
158// {
159// int idx1, idx2;
160// idx1 = neighborsTri.getElement(pId, l);
161// idx2 = tmp;
162// if (neighborsTri.getElement(pId, l) < 0 && tmp < 0)
163// {
164// idx1 *= -1;
165// idx1--;
166// idx2 *= -1;
167// idx2--;
168// Triangle3D t3d1(positionTri[m_triangle_index[idx1][0]], positionTri[m_triangle_index[idx1][1]], positionTri[m_triangle_index[idx1][2]]);
169// Triangle3D t3d2(positionTri[m_triangle_index[idx2][0]], positionTri[m_triangle_index[idx2][1]], positionTri[m_triangle_index[idx2][2]]);
170//
171// Coord normal1 = t3d1.normal().normalize();
172// Coord normal2 = t3d2.normal().normalize();
173//
174// Point3D p3d(pos_i);
175//
176// Plane3D PL1(positionTri[m_triangle_index[idx1][0]], t3d1.normal());
177// Plane3D PL2(positionTri[m_triangle_index[idx2][0]], t3d2.normal());
178//
179// Real dis1 = p3d.distance(PL1);
180// Real dis2 = p3d.distance(PL2);
181//
182// if (abs(dis1 - dis2) < EPSILON && abs(normal1[0] - normal2[0]) < EPSILON && abs(normal1[1] - normal2[1]) < EPSILON)
183// judge = normal1[2] <= normal2[2] ? true : false;
184// else if (abs(dis1 - dis2) < EPSILON && abs(normal1[0] - normal2[0]) < EPSILON)
185// judge = normal1[1] <= normal2[1] ? true : false;
186// else if (abs(dis1 - dis2) < EPSILON)
187// judge = normal1[0] <= normal2[0] ? true : false;
188// else
189// judge = dis1 <= dis2 ? true : false;
190// }
191// else
192// judge = neighborsTri.getElement(pId, l) <= tmp ? true : false;
193// }
194// if (judge)
195// break;
196// else
197// {
198// neighborsTri.setElement(pId, c, neighborsTri.getElement(pId, l));
199// neighborsTri.setElement(pId, l, tmp);
200// }
201// }
202// }
203// for (int ne = nbSizeTri - 1; ne > 0; ne--)
204// {
205// int swap_tmp = neighborsTri.getElement(pId, 0);
206// neighborsTri.setElement(pId, 0, neighborsTri.getElement(pId, ne));
207// neighborsTri.setElement(pId, ne, swap_tmp);
208// int start = 0;
209// int end = ne - 1;
210// int c = start;
211// int l = 2 * c + 1;
212// int tmp = neighborsTri.getElement(pId, c);
213// for (; l <= end; c = l, l = 2 * l + 1)
214// {
215// if (l < end)
216// {
217// bool judge = false;
218// {
219// int idx1, idx2;
220// idx1 = neighborsTri.getElement(pId, l);
221// idx2 = neighborsTri.getElement(pId, l + 1);
222// if (neighborsTri.getElement(pId, l) < 0 && neighborsTri.getElement(pId, l + 1) < 0)
223// {
224// idx1 *= -1;
225// idx1--;
226// idx2 *= -1;
227// idx2--;
228// Triangle3D t3d1(positionTri[m_triangle_index[idx1][0]], positionTri[m_triangle_index[idx1][1]], positionTri[m_triangle_index[idx1][2]]);
229// Triangle3D t3d2(positionTri[m_triangle_index[idx2][0]], positionTri[m_triangle_index[idx2][1]], positionTri[m_triangle_index[idx2][2]]);
230//
231// Coord normal1 = t3d1.normal().normalize();
232// Coord normal2 = t3d2.normal().normalize();
233//
234// Point3D p3d(pos_i);
235//
236// Plane3D PL1(positionTri[m_triangle_index[idx1][0]], t3d1.normal());
237// Plane3D PL2(positionTri[m_triangle_index[idx2][0]], t3d2.normal());
238//
239// Real dis1 = p3d.distance(PL1);
240// Real dis2 = p3d.distance(PL2);
241//
242// if (abs(dis1 - dis2) < EPSILON && abs(normal1[0] - normal2[0]) < EPSILON && abs(normal1[1] - normal2[1]) < EPSILON)
243// judge = normal1[2] < normal2[2] ? true : false;
244// else if (abs(dis1 - dis2) < EPSILON && abs(normal1[0] - normal2[0]) < EPSILON)
245// judge = normal1[1] < normal2[1] ? true : false;
246// else if (abs(dis1 - dis2) < EPSILON)
247// judge = normal1[0] < normal2[0] ? true : false;
248// else
249// judge = dis1 < dis2 ? true : false;
250// }
251// else
252// judge = neighborsTri.getElement(pId, l) < neighborsTri.getElement(pId, l + 1) ? true : false;
253// }
254// if (judge)
255// l++;
256// }
257// bool judge = false;
258// {
259// int idx1, idx2;
260// idx1 = neighborsTri.getElement(pId, l);
261// idx2 = tmp;
262// if (neighborsTri.getElement(pId, l) < 0 && tmp < 0)
263// {
264// idx1 *= -1;
265// idx1--;
266// idx2 *= -1;
267// idx2--;
268// Triangle3D t3d1(positionTri[m_triangle_index[idx1][0]], positionTri[m_triangle_index[idx1][1]], positionTri[m_triangle_index[idx1][2]]);
269// Triangle3D t3d2(positionTri[m_triangle_index[idx2][0]], positionTri[m_triangle_index[idx2][1]], positionTri[m_triangle_index[idx2][2]]);
270//
271// Plane3D PL1(positionTri[m_triangle_index[idx1][0]], t3d1.normal());
272// Plane3D PL2(positionTri[m_triangle_index[idx2][0]], t3d2.normal());
273//
274// Coord normal1 = t3d1.normal().normalize();
275// Coord normal2 = t3d2.normal().normalize();
276//
277// Point3D p3d(pos_i);
278//
279// Real dis1 = p3d.distance(PL1);
280// Real dis2 = p3d.distance(PL2);
281//
282// if (abs(dis1 - dis2) < EPSILON && abs(normal1[0] - normal2[0]) < EPSILON && abs(normal1[1] - normal2[1]) < EPSILON)
283// judge = normal1[2] <= normal2[2] ? true : false;
284// else if (abs(dis1 - dis2) < EPSILON && abs(normal1[0] - normal2[0]) < EPSILON)
285// judge = normal1[1] <= normal2[1] ? true : false;
286// else if (abs(dis1 - dis2) < EPSILON)
287// judge = normal1[0] <= normal2[0] ? true : false;
288// else
289// judge = dis1 <= dis2 ? true : false;
290// }
291// else
292// judge = neighborsTri.getElement(pId, l) <= tmp ? true : false;
293// }
294// if (judge)
295// break;
296// else
297// {
298// neighborsTri.setElement(pId, c, neighborsTri.getElement(pId, l));
299// neighborsTri.setElement(pId, l, tmp);
300// }
301// }
302// }
303// }
304
305template <typename Real, typename Coord>
306__global__ void VC_ComputeAlphaTmp(
307 DArray<Real> alpha,
308 DArray<Real> rho_alpha,
309 DArray<Real> mass,
310 DArray<Coord> position,
311 DArray<TopologyModule::Triangle> m_triangle_index,
312 DArray<Coord> positionTri,
313 DArray<Attribute> attribute,
314 DArrayList<int> neighbors,
315 DArrayList<int> neighborsTri,
316 Real smoothingLength,
317 Real m_sampling_distance,
318 DArray<int> flip)
319{
320 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
321
322 if (pId >= position.size())
323 return;
324 if (!attribute[pId].isDynamic())
325 return;
326
327 Coord pos_i = position[pId];
328 Real alpha_i = 0.0f;
329 Real ra = 0.0f;
330 List<int>& list_i = neighbors[pId];
331 int nbSize = list_i.size();
332
333 List<int>& triList_i = neighborsTri[pId];
334 int nbSizeTri = triList_i.size();
335
336 Real debug_r = smoothingLength;
337 bool tmp = false;
338 Real AreaB;
339 for (int ne = 0; ne < nbSizeTri; ne++)
340 {
341 int j = triList_i[ne];
342 if (j >= 0)
343 continue;
344 j *= -1;
345 j--;
346 // Real m_sampling_distance = 0.015;
347
348 Triangle3D t3d(positionTri[m_triangle_index[j][0]], positionTri[m_triangle_index[j][1]], positionTri[m_triangle_index[j][2]]);
349 Plane3D PL(positionTri[m_triangle_index[j][0]], t3d.normal());
350 Point3D p3d(pos_i);
351 //Point3D nearest_pt = p3d.project(t3d);
352 Point3D nearest_pt = p3d.project(PL);
353 Real r = (nearest_pt.origin - pos_i).norm();
354 tmp = true;
355 float d = p3d.distance(PL);
356
357 Real AreaSum = calculateIntersectionArea(p3d, t3d, smoothingLength);
358 Real MinDistance = abs(p3d.distance(t3d));
359 Coord Min_Pt = (p3d.project(t3d)).origin - pos_i;
360 if (ne < nbSizeTri - 1 && triList_i[ne + 1] < 0)
361 {
362 //triangle clustering
363 int jn;
364 do
365 {
366 jn = triList_i[ne + 1];
367 if (jn > 0)
368 break;
369 jn *= -1;
370 jn--;
371
372 Triangle3D t3d_n(positionTri[m_triangle_index[jn][0]], positionTri[m_triangle_index[jn][1]], positionTri[m_triangle_index[jn][2]]);
373 if ((t3d.normal().cross(t3d_n.normal())).norm() > EPSILON)
374 break;
375
376 AreaSum += calculateIntersectionArea(p3d, t3d_n, smoothingLength);
377
378 if (abs(p3d.distance(t3d_n)) < MinDistance)
379 {
380 MinDistance = abs(p3d.distance(t3d_n));
381 Min_Pt = (p3d.project(t3d_n)).origin - pos_i;
382 }
383 //printf("%d %d\n", j, jn);
384 ne++;
385 } while (ne < nbSizeTri - 1);
386 }
387 Min_Pt /= Min_Pt.norm();
388
389 d = abs(d);
390 if (smoothingLength - d > EPSILON && smoothingLength * smoothingLength - d * d > EPSILON && d > EPSILON)
391 {
392 Coord n_PL = -t3d.normal();
393 if (flip[pId] < 0)
394 n_PL *= -1;
395 Coord n_TR = (p3d.project(t3d)).origin - pos_i;
396 n_PL = n_PL / n_PL.norm();
397 n_TR = n_TR / n_TR.norm();
398
399 AreaB = M_PI * (smoothingLength * smoothingLength - d * d);
400 //equation 12
401 Real a_ij =
402 kernWeightMesh(r, smoothingLength)
403 * 2.0 * (M_PI) * (1 - d / smoothingLength)
404 //* p3d.areaTriangle(t3d, smoothingLength) * n_PL.dot(n_TR)
405 * AreaSum * n_PL.dot(Min_Pt)
406 / ((M_PI) * (smoothingLength * smoothingLength - d * d)) / (m_sampling_distance * m_sampling_distance * m_sampling_distance);
407 debug_r = min(r, debug_r);
408 alpha_i += a_ij;
409 }
410 }
411
412 Real alpha_solid = alpha_i;
413 Real GT_alpha = Real(0);
414 for (int ne = 0; ne < nbSize; ne++)
415 {
416 int j = list_i[ne];
417 Real r = (pos_i - position[j]).norm();
418 ;
419
420 //if (r > EPSILON)
421 if (r > EPSILON && attribute[j].isDynamic())
422 {
423 Real a_ij = kernWeight(r, smoothingLength);
424 alpha_i += a_ij;
425 }
426 else if (r > EPSILON)
427 {
428 //printf("Alpha r: %.3lf posj: %.3lf %.3lf %.3lf \n", r,position[j][0], position[j][1],position[j][2]);
429 Real a_ij = kernWeight(r, smoothingLength);
430 GT_alpha += a_ij;
431 }
432 }
433
434 alpha[pId] = alpha_i;
435 rho_alpha[pId] = ra;
436}
437
438template <typename Real>
439__global__ void VC_CorrectAlphaTmp(
440 DArray<Real> alpha,
441 DArray<Real> rho_alpha,
442 DArray<Real> mass,
443 Real maxAlpha)
444{
445 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
446 if (pId >= alpha.size())
447 return;
448
449 Real alpha_i = alpha[pId];
450 Real ra = rho_alpha[pId];
451 if (alpha_i < maxAlpha)
452 {
453 ra += (maxAlpha - alpha_i) * mass[pId];
454 alpha_i = maxAlpha;
455 }
456 alpha[pId] = alpha_i;
457 rho_alpha[pId] = ra;
458}
459
460template <typename Real, typename Coord>
461__global__ void VC_ComputeDiagonalElementTmp(
462 DArray<Real> AiiFluid,
463 DArray<Real> AiiTotal,
464 DArray<Real> alpha,
465 DArray<Coord> position,
466 DArray<TopologyModule::Triangle> m_triangle_index,
467 DArray<Coord> positionTri,
468 DArray<Attribute> attribute,
469 DArrayList<int> neighbors,
470 DArrayList<int> neighborsTri,
471 Real smoothingLength,
472 Real m_sampling_distance,
473 DArray<int> flip)
474{
475 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
476 if (pId >= position.size())
477 return;
478
479 //printf("Yes\n");
480
481 if (!attribute[pId].isDynamic())
482 return;
483
484 Real invAlpha = 1.0f / alpha[pId];
485
486 Real diaA_total = 0.0f;
487 Real diaA_fluid = 0.0f;
488 Real diaA_test = 0.0f;
489 Coord pos_i = position[pId];
490
491 bool bNearWall = false;
492 List<int>& list_i = neighbors[pId];
493 int nbSize = list_i.size();
494
495 int nbSizeTri;
496 List<int>& triList_i = neighborsTri[pId];
497 nbSizeTri = triList_i.size();
498 float dd;
499 for (int ne = 0; ne < nbSizeTri; ne++)
500 {
501 int j = triList_i[ne];
502 if (j >= 0)
503 continue;
504 j *= -1;
505 j--;
506 //Real m_sampling_distance = 0.015;
507
508 Triangle3D t3d(positionTri[m_triangle_index[j][0]], positionTri[m_triangle_index[j][1]], positionTri[m_triangle_index[j][2]]);
509 Plane3D PL(positionTri[m_triangle_index[j][0]], t3d.normal());
510 Point3D p3d(pos_i);
511 //Point3D nearest_pt = p3d.project(t3d);
512 Point3D nearest_pt = p3d.project(PL);
513 Real r = (nearest_pt.origin - pos_i).norm();
514 float d = p3d.distance(PL);
515 //if (d < 0) continue;
516 d = abs(d);
517 dd = d;
518
519 Real AreaSum = calculateIntersectionArea(p3d, t3d, smoothingLength);
520 Real MinDistance = abs(p3d.distance(t3d));
521 Coord Min_Pt = (p3d.project(t3d)).origin - pos_i;
522 if (ne < nbSizeTri - 1 && triList_i[ne + 1] < 0)
523 {
524 //triangle clustering
525 int jn;
526 do
527 {
528 jn = triList_i[ne + 1];
529 if (jn > 0)
530 break;
531 jn *= -1;
532 jn--;
533
534 Triangle3D t3d_n(positionTri[m_triangle_index[jn][0]], positionTri[m_triangle_index[jn][1]], positionTri[m_triangle_index[jn][2]]);
535 if ((t3d.normal().cross(t3d_n.normal())).norm() > EPSILON)
536 break;
537
538 AreaSum += calculateIntersectionArea(p3d, t3d_n, smoothingLength);
539
540 if (abs(p3d.distance(t3d_n)) < MinDistance)
541 {
542 MinDistance = abs(p3d.distance(t3d_n));
543 Min_Pt = (p3d.project(t3d_n)).origin - pos_i;
544 }
545 //printf("%d %d\n", j, jn);
546 ne++;
547 } while (ne < nbSizeTri - 1);
548 }
549 Min_Pt /= Min_Pt.norm();
550 if (smoothingLength - d > EPSILON && smoothingLength * smoothingLength - d * d > EPSILON && d > EPSILON)
551 {
552 //equation 12
553
554 Coord n_PL = -t3d.normal();
555 if (flip[pId] < 0)
556 n_PL *= -1;
557 Coord n_TR = (p3d.project(t3d)).origin - pos_i;
558 n_PL = n_PL / n_PL.norm();
559 n_TR = n_TR / n_TR.norm();
560 Real wrr_ij = invAlpha * kernWRRMesh(r, smoothingLength)
561 * 2.0 * (M_PI) * (1 - d / smoothingLength)
562 //* p3d.areaTriangle(t3d, smoothingLength) * n_PL.dot(n_TR)
563 * AreaSum * n_PL.dot(Min_Pt)
564 / ((M_PI) * (smoothingLength * smoothingLength - d * d))
565 / (m_sampling_distance * m_sampling_distance * m_sampling_distance);
566 // printf("++++++++++++++++++++++++++++++++++++%.3lf\n", wrr_ij);
567 if (abs(wrr_ij) > 3000)
568 {
569 // printf("WRR: %.3lf ANGLE: %.7lf distance:%.3lf weight:%.3lf\n", kernWRRMesh(r, smoothingLength), (smoothingLength * smoothingLength - d * d),
570 // d,wrr_ij);
571 }
572 diaA_total += 2.0f * wrr_ij;
573 diaA_test += 2.0f * wrr_ij;
574 }
575 }
576 Real diaA_GT = Real(0);
577 for (int ne = 0; ne < nbSize; ne++)
578 {
579 int j = list_i[ne];
580 Real r = (pos_i - position[j]).norm();
581
582 Attribute att_j = attribute[j];
583 if (r > EPSILON)
584 {
585 Real wrr_ij = invAlpha * kernWRR(r, smoothingLength);
586 if (att_j.isDynamic())
587 {
588 diaA_total += wrr_ij;
589 diaA_fluid += wrr_ij;
590 atomicAdd(&AiiFluid[j], wrr_ij);
591 atomicAdd(&AiiTotal[j], wrr_ij);
592 }
593 else
594 {
595 //diaA_total += 2.0f * wrr_ij;
596 diaA_GT += 2.0f * wrr_ij;
597 }
598 }
599 }
600
601 atomicAdd(&AiiFluid[pId], diaA_fluid);
602 atomicAdd(&AiiTotal[pId], diaA_total);
603 //if (abs(diaA_test) > EPSILON)
604 // printf("========diaA_FOR_TEST: %.3lf GT:%.3lf DISTANCE: %.3lf\n",diaA_test, diaA_GT, dd);
605}
606
607template <typename Real, typename Coord>
608__global__ void VC_ComputeDiagonalElementTmp(
609 DArray<Real> diaA,
610 DArray<Real> alpha,
611 DArray<Coord> position,
612 DArray<Attribute> attribute,
613 DArrayList<int> neighbors,
614 Real smoothingLength)
615{
616 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
617 if (pId >= position.size())
618 return;
619 if (!attribute[pId].isDynamic())
620 return;
621
622 Coord pos_i = position[pId];
623 Real invAlpha_i = 1.0f / alpha[pId];
624 Real A_i = 0.0f;
625
626 List<int>& list_i = neighbors[pId];
627 int nbSize = list_i.size();
628
629 for (int ne = 0; ne < nbSize; ne++)
630 {
631 int j = list_i[ne];
632 Real r = (pos_i - position[j]).norm();
633
634 if (r > EPSILON && attribute[j].isDynamic())
635 {
636 Real wrr_ij = invAlpha_i * kernWRR(r, smoothingLength);
637 A_i += wrr_ij;
638 atomicAdd(&diaA[j], wrr_ij);
639 }
640 }
641
642 atomicAdd(&diaA[pId], A_i);
643}
644
645template <typename Real, typename Coord>
646__global__ void VC_DetectSurfaceTmp(
647 DArray<Real> Aii,
648 DArray<bool> bSurface,
649 DArray<Real> AiiFluid,
650 DArray<Real> AiiTotal,
651 DArray<Coord> position,
652 DArray<Attribute> attribute,
653 DArrayList<int> neighbors,
654 Real smoothingLength,
655 Real maxA)
656{
657 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
658 if (pId >= position.size())
659 return;
660 if (!attribute[pId].isDynamic())
661 return;
662
663 Real total_weight = 0.0f;
664 Coord div_i = Coord(0);
665
666 SmoothKernel<Real> kernSmooth;
667
668 Coord pos_i = position[pId];
669 List<int>& list_i = neighbors[pId];
670 int nbSize = list_i.size();
671 bool bNearWall = false;
672 for (int ne = 0; ne < nbSize; ne++)
673 {
674 int j = list_i[ne];
675 Real r = (pos_i - position[j]).norm();
676
677 if (r > EPSILON && attribute[j].isDynamic())
678 {
679 float weight = -kernSmooth.Gradient(r, smoothingLength);
680 total_weight += weight;
681 div_i += (position[j] - pos_i) * (weight / r);
682 }
683
684 if (!attribute[j].isDynamic())
685 {
686 //bNearWall = true;
687 }
688 }
689
690 total_weight = total_weight < EPSILON ? 1.0f : total_weight;
691 Real absDiv = div_i.norm() / total_weight;
692
693 bool bSurface_i = false;
694 Real diagF_i = AiiFluid[pId];
695 Real diagT_i = AiiTotal[pId];
696
697 if (abs(diagT_i - diagF_i) > EPSILON)
698 bNearWall = true;
699
700 Real aii = diagF_i;
701 Real eps = 0.001f;
702 Real diagS_i = diagT_i - diagF_i;
703 Real threshold = 0.0f;
704 if (bNearWall && diagT_i < maxA * (1.0f - threshold))
705 {
706 bSurface_i = true;
707 aii = maxA - (diagT_i - diagF_i);
708 }
709
710 if (!bNearWall && diagF_i < maxA * (1.0f - threshold))
711 {
712 bSurface_i = true;
713 aii = maxA;
714 }
715 bSurface[pId] = bSurface_i;
716 Aii[pId] = aii;
717}
718
719template <typename Real, typename Coord>
720__global__ void VC_ComputeDivergenceTmp(
721 DArray<Real> divergence,
722 DArray<Real> alpha,
723 DArray<Real> density,
724 DArray<Coord> position,
725 DArray<Coord> velocity,
726 DArray<Coord> velocityTri,
727 DArray<TopologyModule::Triangle> m_triangle_index,
728 DArray<Coord> positionTri,
729 DArray<bool> bSurface,
730 DArray<Attribute> attribute,
731 DArray<Real> mass,
732 DArray<Real> m_triangle_vertex_mass,
733 DArrayList<int> neighbors,
734 DArrayList<int> neighborsTri,
735 Real separation,
736 Real tangential,
737 Real restDensity,
738 Real smoothingLength,
739 Real m_sampling_distance,
740 Real dt,
741 DArray<int> flip)
742{
743 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
744 if (pId >= position.size())
745 return;
746 if (!attribute[pId].isDynamic())
747 return;
748
749 Coord pos_i = position[pId];
750 Coord vel_i = velocity[pId];
751
752 Real div_vi = 0.0f;
753
754 Real invAlpha_i = 1.0f / alpha[pId];
755 Real mass_i = mass[pId];
756
757 //mass_i = 10.0;
758 Real mass_i0 = mass_i;
759
760 int nbSizeTri;
761
762 List<int>& triList_i = neighborsTri[pId];
763 nbSizeTri = triList_i.size();
764
765 float div_debug, ddb;
766 div_debug = float(0.0);
767 int pop = 0;
768
769 Real DIV_GT, DIV_FL;
770 DIV_GT = Real(0);
771 DIV_FL = Real(0);
772
773 Real sum_weight_norm = Real(0);
774 Coord average_normal_j = Coord(0);
775 for (int ne = 0; ne < nbSizeTri; ne++)
776 {
777 int j = triList_i[ne];
778
779 if (j >= 0)
780 continue;
781 j *= -1;
782 j--;
783
784 Triangle3D t3d(positionTri[m_triangle_index[j][0]], positionTri[m_triangle_index[j][1]], positionTri[m_triangle_index[j][2]]);
785 Plane3D PL(positionTri[m_triangle_index[j][0]], t3d.normal());
786 Point3D p3d(pos_i);
787 //Point3D nearest_pt = p3d.project(t3d);
788 Point3D nearest_pt = p3d.project(PL);
789 Point3D plnpt = p3d.project(PL);
790 Real r = (nearest_pt.origin - pos_i).norm();
791 float d = p3d.distance(PL);
792 d = abs(d);
793 ddb = d;
794 if (r < EPSILON)
795 continue;
796 if (smoothingLength - d > EPSILON && smoothingLength * smoothingLength - d * d > EPSILON && d > EPSILON)
797 {
798
799 pop = 1;
800 Coord n_PL = nearest_pt.origin - pos_i;
801 Coord n_TR = (p3d.project(t3d)).origin - pos_i;
802 n_PL = n_PL / n_PL.norm();
803 n_TR = n_TR / n_TR.norm();
804
805 float wr_ij =
806 kernWeightMesh(r, smoothingLength)
807 * 2.0 * (M_PI) * (1 - d / smoothingLength)
808 * calculateIntersectionArea(p3d, t3d, smoothingLength) //* n_PL.dot(n_TR)
809 / ((M_PI) * (smoothingLength * smoothingLength - d * d))
810 / (m_sampling_distance * m_sampling_distance * m_sampling_distance);
811 Coord normal_j = t3d.normal().normalize();
812 //sum_weight_norm += wr_ij;
813 average_normal_j += wr_ij * normal_j;
814 }
815 }
816 average_normal_j = average_normal_j.normalize();
817
818 for (int ne = 0; ne < nbSizeTri; ne++)
819 {
820 int j = triList_i[ne];
821 //printf("%d\n",j);
822 if (j >= 0)
823 continue;
824 j *= -1;
825 j--;
826 //Real m_sampling_distance = 0.015;
827 //printf("YESSSSSSSSSSSSSS\n");
828 Triangle3D t3d(positionTri[m_triangle_index[j][0]], positionTri[m_triangle_index[j][1]], positionTri[m_triangle_index[j][2]]);
829 Plane3D PL(positionTri[m_triangle_index[j][0]], t3d.normal());
830 Point3D p3d(pos_i);
831 //Point3D nearest_pt = p3d.project(t3d);
832 Point3D nearest_pt = p3d.project(PL);
833 Point3D plnpt = p3d.project(PL);
834 Real r = (nearest_pt.origin - pos_i).norm();
835 float d = p3d.distance(PL);
836 d = abs(d);
837 ddb = d;
838
839 Real AreaSum = calculateIntersectionArea(p3d, t3d, smoothingLength);
840 Real MinDistance = abs(p3d.distance(t3d));
841 Coord Min_Pt = (p3d.project(t3d)).origin - pos_i;
842 if (ne < nbSizeTri - 1 && triList_i[ne + 1] < 0)
843 {
844 //triangle clustering
845 int jn;
846 do
847 {
848 jn = triList_i[ne + 1];
849 if (jn > 0)
850 break;
851 jn *= -1;
852 jn--;
853
854 Triangle3D t3d_n(positionTri[m_triangle_index[jn][0]], positionTri[m_triangle_index[jn][1]], positionTri[m_triangle_index[jn][2]]);
855 if ((t3d.normal().cross(t3d_n.normal())).norm() > EPSILON)
856 break;
857
858 AreaSum += calculateIntersectionArea(p3d, t3d_n, smoothingLength);
859
860 if (abs(p3d.distance(t3d_n)) < MinDistance)
861 {
862 MinDistance = abs(p3d.distance(t3d_n));
863 Min_Pt = (p3d.project(t3d_n)).origin - pos_i;
864 }
865 //printf("%d %d\n", j, jn);
866 ne++;
867 } while (ne < nbSizeTri - 1);
868 }
869 Min_Pt /= Min_Pt.norm();
870
871 if (r < EPSILON)
872 continue;
873 if (smoothingLength - d > EPSILON && smoothingLength * smoothingLength - d * d > EPSILON && d > EPSILON)
874 {
875
876 //equation 12
877 pop = 1;
878
879 //Coord n_PL = nearest_pt.origin - pos_i;
880 Coord n_PL = -t3d.normal();
881 if (flip[pId] < 0)
882 n_PL *= -1;
883 Coord n_TR = (p3d.project(t3d)).origin - pos_i;
884 n_PL = n_PL / n_PL.norm();
885 n_TR = n_TR / n_TR.norm();
886
887 float wr_ij =
888 kernWRMesh(r, smoothingLength)
889 * 2.0 * (M_PI) * (1 - d / smoothingLength)
890 // * p3d.areaTriangle(t3d, smoothingLength) * n_PL.dot(n_TR)
891 * AreaSum * n_PL.dot(Min_Pt)
892 / ((M_PI) * (smoothingLength * smoothingLength - d * d))
893 / (m_sampling_distance * m_sampling_distance * m_sampling_distance);
894 //printf("WRIJ IN DIV: %.3lf DISTANCE: %.3lf Kern:%.3lf areaTri: %.3lf\n", wr_ij,d, kernWRMesh(r, smoothingLength), p3d.areaTriangle(t3d, smoothingLength));
895
896 Coord g = -invAlpha_i * (pos_i - nearest_pt.origin) * wr_ij * (1.0f / r);
897 if (r < EPSILON)
898 g = -invAlpha_i * wr_ij * t3d.normal().normalize();
899
900 if (!((g.norm()) < 100000000000.0))
901 {
902 printf("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ %.10lf kermMesh:%.3lf areaTri: %.3lf AreaSphere%.3lf~~~~~~~~~~~~~~\n",
903 wr_ij,
904 kernWRMesh(r, smoothingLength),
905 calculateIntersectionArea(p3d, t3d, smoothingLength),
906 (smoothingLength * smoothingLength - d * d));
907 }
908
909 Coord Velj = (velocityTri[m_triangle_index[j][0]] + velocityTri[m_triangle_index[j][1]] + velocityTri[m_triangle_index[j][2]]) / 3.0;
910 //printf("J: %d Velj: %.3lf %.3lf %.3lf Pos:%.3lf %.3lf %.3lf\n",j ,positionTri[m_triangle_index[j][0]][0], positionTri[m_triangle_index[j][0]][1], positionTri[m_triangle_index[j][0]][2],
911 // position[pId][0], position[pId][1], position[pId][2]);
912
913 Real mass_j = m_triangle_vertex_mass[m_triangle_index[j][0]];
914
915 Coord normal_j = t3d.normal();
916 normal_j = normal_j.normalize();
917 //normal_j = average_normal_j;
918 //if (normal_j.dot(pos_i - nearest_pt.origin) < 0)
919 if (flip[pId] < 0)
920 n_PL *= -1;
921 normal_j *= -1;
922 //
923 // printf("NORMAL_J: %.3lf %.3lf %.3lf nij: %.3lf %.3lf %.3lf\n",normal_j[0], normal_j[1], normal_j[2], nij[0], nij[1], nij[2]);
924
925 Coord dVel = vel_i - Velj;
926 Real magNVel = dVel.dot(normal_j);
927 Coord nVel = magNVel * normal_j;
928 Coord tVel = dVel - nVel;
929
930 if (magNVel < -EPSILON)
931 {
932 Real div_ij = g.dot(2.0f * (nVel + tangential * tVel)) * restDensity * mass_i / dt;
933 atomicAdd(&divergence[pId], div_ij);
934 div_debug += div_ij;
935 // if (pos_i[1] < 0.035 && pId % 300 == 0)
936 // printf("down, %.3lf %.3lf %.3lf \n", nVel[0], nVel[1], nVel[2]);
937 }
938 else
939 {
940 Real div_ij = g.dot(2.0f * (( separation )*nVel + tangential * tVel)) * mass_i * restDensity / dt;
941 atomicAdd(&divergence[pId], div_ij);
942 div_debug += div_ij;
943 // if (pos_i[1] < 0.035 && pId % 300 == 0)
944 // printf("up, %.3lf %.3lf %.3lf \n",nVel[0],nVel[1],nVel[2]);
945 }
946 }
947 }
948
949 List<int>& list_i = neighbors[pId];
950 int nbSize = list_i.size();
951 for (int ne = 0; ne < nbSize; ne++)
952 {
953 int j = list_i[ne];
954 Real r = (pos_i - position[j]).norm();
955 Real mass_j = mass[j];
956
957 if (r > EPSILON)
958 {
959 Real wr_ij = kernWR(r, smoothingLength);
960 Coord g = -invAlpha_i * (pos_i - position[j]) * wr_ij * (1.0f / r);
961 //mass_i = min(mass_i0, mass_j);
962 if ((attribute[j].isDynamic() && attribute[j].isRigid() && attribute[pId].isRigid()))
963 {
964 }
965 //else if(attribute[j].IsDynamic())
966 else if (attribute[j].isFluid() && attribute[pId].isFluid())
967 {
968 Real div_ij = 0.5f * (vel_i - velocity[j]).dot(g) * ( mass_i )*restDensity / dt;
969
970 atomicAdd(&divergence[pId], div_ij);
971 atomicAdd(&divergence[j], div_ij);
972 }
973 }
974 }
975
976 if ((!(abs(div_debug) < EPSILON)))
977 {
978 // printf("FROM DIV: %.3lf GT: %.3lf\n", div_debug, divergence[pId] - div_debug);
979 }
980}
981
982template <typename Real, typename Coord>
983__global__ void VC_CompensateSourceTmp(
984 DArray<Real> divergence,
985 DArray<Real> density,
986 DArray<Attribute> attribute,
987 DArray<Coord> position,
988 Real restDensity,
989 Real dt)
990{
991 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
992 if (pId >= density.size())
993 return;
994 if (!attribute[pId].isDynamic())
995 return;
996
997 Coord pos_i = position[pId];
998 if (density[pId] > restDensity)
999 {
1000 Real ratio = (density[pId] - restDensity) / restDensity;
1001 atomicAdd(&divergence[pId], 1000000.0f * ratio / dt);
1002 }
1003}
1004
1005// compute Ax;
1006template <typename Real, typename Coord>
1007__global__ void VC_ComputeAxTmp(
1008 DArray<Real> residual,
1009 DArray<Real> pressure,
1010 DArray<Real> aiiSymArr,
1011 DArray<Real> alpha,
1012 DArray<Coord> position,
1013 DArray<Attribute> attribute,
1014 DArrayList<int> neighbor,
1015 Real smoothingLength)
1016{
1017 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
1018 if (pId >= position.size())
1019 return;
1020 if (!attribute[pId].isDynamic())
1021 return;
1022
1023 Coord pos_i = position[pId];
1024 Real invAlpha_i = 1.0f / alpha[pId];
1025
1026 atomicAdd(&residual[pId], aiiSymArr[pId] * pressure[pId]);
1027 Real con1 = 1.0f; // PARAMS.mass / PARAMS.restDensity / PARAMS.restDensity;
1028
1029 List<int>& list_i = neighbor[pId];
1030 int nbSize = list_i.size();
1031 for (int ne = 0; ne < nbSize; ne++)
1032 {
1033 int j = list_i[ne];
1034 Real r = (pos_i - position[j]).norm();
1035
1036 if (r > EPSILON && attribute[j].isDynamic())
1037 {
1038 Real wrr_ij = kernWRR(r, smoothingLength);
1039 Real a_ij = -invAlpha_i * wrr_ij;
1040 // residual += con1*a_ij*preArr[j];
1041 atomicAdd(&residual[pId], con1 * a_ij * pressure[j]);
1042 atomicAdd(&residual[j], con1 * a_ij * pressure[pId]);
1043 }
1044 }
1045}
1046
1047template <typename Real>
1048__global__ void VC_InitAttrTmp(
1049 DArray<Attribute> attribute,
1050 DArray<Real> mass)
1051{
1052 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
1053 if (pId >= attribute.size())
1054 return;
1055 attribute[pId].setDynamic();
1056 attribute[pId].setFluid();
1057 mass[pId] = 10.0;
1058}
1059
1060template <typename Coord, typename Real>
1061__global__ void VC_TriVelTmp(
1062 DArray<Coord> oldPos,
1063 DArray<Coord> newPos,
1064 DArray<Coord> vels,
1065 Real dt)
1066{
1067 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
1068 //printf("YES\n");
1069 if (pId >= oldPos.size())
1070 return;
1071 vels[pId] = (newPos[pId] - oldPos[pId]) / dt;
1072}
1073
1074template <typename Real, typename Coord>
1075__global__ void VC_UpdateVelocityBoundaryCorrectedTmp(
1076 DArray<Real> pressure,
1077 DArray<Real> alpha,
1078 DArray<bool> bSurface,
1079 DArray<Coord> position,
1080 DArray<Coord> velocity,
1081 DArray<Coord> velocityTri,
1082 DArray<TopologyModule::Triangle> m_triangle_index,
1083 DArray<Coord> positionTri,
1084 DArray<Attribute> attribute,
1085 DArray<Real> mass,
1086 DArray<Real> m_triangle_vertex_mass,
1087 DArrayList<int> neighbor,
1088 DArrayList<int> neighborTri,
1089 Real restDensity,
1090 Real airPressure,
1091 Real sliding,
1092 Real separation,
1093 Real smoothingLength,
1094 Real m_sampling_distance,
1095 Real dt,
1096 DArray<int> flip)
1097{
1098 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
1099 if (pId >= position.size())
1100 return;
1101 //printf("%d\n", position.size());
1102
1103 if (attribute[pId].isDynamic())
1104 {
1105 Attribute att_i = attribute[pId];
1106
1107 Coord pos_i = position[pId];
1108 Real p_i = pressure[pId];
1109
1110 List<int>& list_i = neighbor[pId];
1111 int nbSize = list_i.size();
1112
1113 List<int>& triList_i = neighborTri[pId];
1114 int nbSizeTri = triList_i.size();
1115 Real total_weight = 0.0f;
1116
1117 Real ceo = 1.6f;
1118
1119 Real invAlpha = 1.0f / alpha[pId];
1120 Real invAlpha_i = invAlpha;
1121
1122 Real mass_i = mass[pId];
1123
1124 Real mass_i0 = mass_i;
1125 Coord vel_i = velocity[pId];
1126 Coord dv_i(0.0f);
1127 Real scale = 1.0f * dt / restDensity;
1128 Real acuP = 0.0f;
1129 total_weight = 0.0f;
1130
1131 Real sum_weight_norm = Real(0);
1132
1133 bool tmp = false;
1134 for (int ne = 0; ne < nbSizeTri; ne++)
1135 {
1136 int j = triList_i[ne];
1137 if (j >= 0)
1138 continue;
1139 j *= -1;
1140 j--;
1141
1142 Triangle3D t3d(positionTri[m_triangle_index[j][0]], positionTri[m_triangle_index[j][1]], positionTri[m_triangle_index[j][2]]);
1143 Plane3D PL(positionTri[m_triangle_index[j][0]], t3d.normal());
1144 Point3D p3d(pos_i);
1145 Point3D nearest_ptt = p3d.project(t3d);
1146 Point3D nearest_pt = p3d.project(PL);
1147 Point3D plnpt = p3d.project(PL);
1148 Real r = (nearest_pt.origin - pos_i).norm();
1149
1150 if (r < EPSILON)
1151 continue;
1152 float d = p3d.distance(PL);
1153 //if (d < 0) continue;
1154 d = abs(d);
1155
1156 Real AreaSum = calculateIntersectionArea(p3d, t3d, smoothingLength);
1157 Real MinDistance = abs(p3d.distance(t3d));
1158 Coord Min_Pt = (p3d.project(t3d)).origin - pos_i;
1159 if (ne < nbSizeTri - 1 && triList_i[ne + 1] < 0)
1160 {
1161 //triangle clustering
1162 int jn;
1163 do
1164 {
1165 jn = triList_i[ne + 1];
1166 if (jn > 0)
1167 break;
1168 jn *= -1;
1169 jn--;
1170
1171 Triangle3D t3d_n(positionTri[m_triangle_index[jn][0]], positionTri[m_triangle_index[jn][1]], positionTri[m_triangle_index[jn][2]]);
1172 if ((t3d.normal().cross(t3d_n.normal())).norm() > EPSILON)
1173 break;
1174
1175 AreaSum += calculateIntersectionArea(p3d, t3d_n, smoothingLength);
1176
1177 if (abs(p3d.distance(t3d_n)) < MinDistance)
1178 {
1179 MinDistance = abs(p3d.distance(t3d_n));
1180 Min_Pt = (p3d.project(t3d_n)).origin - pos_i;
1181 }
1182 //printf("%d %d\n", j, jn);
1183 ne++;
1184 } while (ne < nbSizeTri - 1);
1185 }
1186 Min_Pt /= Min_Pt.norm();
1187
1188 //Coord n_PL = nearest_pt.origin - pos_i;
1189 Coord n_PL = -t3d.normal();
1190 //if (flip[pId] < 0) n_PL *= -1;
1191 Coord n_TR = (p3d.project(t3d)).origin - pos_i;
1192 n_PL = n_PL / n_PL.norm();
1193 n_TR = n_TR / n_TR.norm();
1194
1195 // if (pos_i[2] > 0.98 && pos_i[0] < 0.9)
1196 // printf("%.3lf %.3lf %.3lf %d\n", t3d.normal()[0], t3d.normal()[1], t3d.normal()[2], flip[pId]);
1197
1198 Real weight = -invAlpha * kernWRMesh(r, smoothingLength)
1199 * 2.0 * (M_PI) * (1 - d / smoothingLength)
1200 //* p3d.areaTriangle(t3d, smoothingLength) * n_PL.dot(n_TR)
1201 * AreaSum * n_PL.dot(Min_Pt)
1202 / ((M_PI) * (smoothingLength * smoothingLength - d * d))
1203 / (m_sampling_distance * m_sampling_distance * m_sampling_distance);
1204
1205 Coord dnij = (pos_i - nearest_pt.origin) * (1.0f / r);
1206 if (r < EPSILON)
1207 dnij = t3d.normal().normalize();
1208
1209 Coord corrected = dnij;
1210 if (corrected.norm() > EPSILON)
1211 {
1212 corrected = corrected.normalize();
1213 }
1214
1215 Real mass_j = m_triangle_vertex_mass[m_triangle_index[j][0]];
1216 corrected = -scale * weight * corrected / (mass_i);
1217 d = p3d.distance(PL);
1218 //if (d < 0) continue;
1219 d = abs(d);
1220 if (smoothingLength - d > EPSILON && smoothingLength * smoothingLength - d * d > EPSILON && d > EPSILON)
1221 {
1222 //equation 12
1223 Coord dvii = 2.0f * (pressure[pId]) * corrected;
1224 if (bSurface[pId])
1225 {
1226 dv_i += dvii;
1227 }
1228 float weight = 2.0f * invAlpha
1229 * kernWeightMesh(r, smoothingLength)
1230 * 2.0 * (M_PI) * (1 - d / smoothingLength)
1231 * AreaSum * n_PL.dot(Min_Pt)
1232 / ((M_PI) * (smoothingLength * smoothingLength - d * d))
1233 / (m_sampling_distance * m_sampling_distance * m_sampling_distance);
1234
1235 Coord nij = (pos_i - nearest_pt.origin);
1236 if (nij.norm() > EPSILON)
1237 {
1238 nij = nij.normalize();
1239 }
1240 else
1241 nij = t3d.normal().normalize();
1242
1243 Coord normal_j = t3d.normal();
1244 normal_j = normal_j.normalize();
1245 //if (flip[pId] < 0) normal_j *= -1;
1246
1247 Coord Velj = Coord(0); //(velocityTri[m_triangle_index[j][0]] + velocityTri[m_triangle_index[j][1]] + velocityTri[m_triangle_index[j][2]]) / 3.0;
1248
1249 Coord dVel = Velj - vel_i; //!!!!!!!!!!
1250 Real magNVel = dVel.dot(normal_j);
1251 Coord nVel = magNVel * normal_j;
1252 Coord tVel = dVel - nVel;
1253
1254 //if (pos_i[2] > 0.98 && pos_i[0] < 0.9 && pos_i[2] < 1.0)
1255 // printf("%.3lf %.3lf %.3lf %.3lf\n", nij[0], nij[1], nij[2], weight);
1256
1257 if (magNVel > EPSILON)
1258 {
1259 dv_i += weight * nij.dot(nVel + sliding * tVel) * nij;
1260 //dv_i += weight * (nVel + sliding * tVel);
1261 }
1262 else
1263 {
1264 dv_i += weight * nij.dot(separation * nVel + sliding * tVel) * nij;
1265 //dv_i += weight * (separation * nVel + sliding * tVel);
1266 }
1267 }
1268 }
1269 for (int ne = 0; ne < nbSize; ne++)
1270 {
1271 int j = list_i[ne];
1272 Real r = (pos_i - position[j]).norm();
1273
1274 Attribute att_j = attribute[j];
1275
1276 Real mass_j = mass[j];
1277
1278 if (r > EPSILON)
1279 {
1280 Real weight = -invAlpha * kernWR(r, smoothingLength);
1281 Coord dnij = (pos_i - position[j]) * (1.0f / r);
1282 Coord corrected = dnij;
1283 if (corrected.norm() > EPSILON)
1284 {
1285 corrected = corrected.normalize();
1286 }
1287
1288 Real mass_j = mass[j];
1289 corrected = -scale * weight * corrected / (mass_i);
1290
1291 Coord dvij = (pressure[j] - pressure[pId]) * corrected;
1292 Coord dvjj = (pressure[j] + airPressure) * corrected;
1293
1294 Coord dvij_sym = 0.5f * (pressure[pId] + pressure[j]) * corrected;
1295
1296 Coord dVel = (velocity[j] - vel_i);
1297
1298 Real ratio_mass = mass_j / (mass_i + mass_j);
1299 float weight_m = invAlpha * kernWeight(r, smoothingLength);
1300
1301 if (att_j.isDynamic())
1302 {
1303 if (bSurface[pId])
1304 dv_i += dvjj;
1305 else
1306 dv_i += dvij;
1307
1308 if (bSurface[j])
1309 {
1310 Coord dvii = -(pressure[pId] + airPressure) * corrected;
1311 atomicAdd(&velocity[j][0], ceo * dvii[0]);
1312 atomicAdd(&velocity[j][1], ceo * dvii[1]);
1313 atomicAdd(&velocity[j][2], ceo * dvii[2]);
1314 }
1315 else
1316 {
1317 atomicAdd(&velocity[j][0], ceo * dvij[0]);
1318 atomicAdd(&velocity[j][1], ceo * dvij[1]);
1319 atomicAdd(&velocity[j][2], ceo * dvij[2]);
1320 }
1321 }
1322 }
1323 }
1324
1325 dv_i *= ceo;
1326 if (attribute[pId].isFluid())
1327 {
1328 atomicAdd(&velocity[pId][0], dv_i[0]);
1329 atomicAdd(&velocity[pId][1], dv_i[1]);
1330 atomicAdd(&velocity[pId][2], dv_i[2]);
1331 }
1332 }
1333}
1334
1335template <typename TDataType>
1336SemiAnalyticalIncompressibilityModule<TDataType>::SemiAnalyticalIncompressibilityModule()
1337 : ConstraintModule()
1338 , m_airPressure(Real(00.0))
1339 , m_reduce(NULL)
1340 , m_arithmetic(NULL)
1341{
1342 m_smoothing_length.setValue(Real(0.011));
1343}
1344
1345template <typename TDataType>
1346SemiAnalyticalIncompressibilityModule<TDataType>::~SemiAnalyticalIncompressibilityModule()
1347{
1348 m_alpha.clear();
1349 m_Aii.clear();
1350 m_AiiFluid.clear();
1351 m_AiiTotal.clear();
1352 m_pressure.clear();
1353 //m_pressure2.release();
1354 m_divergence.clear();
1355 m_bSurface.clear();
1356
1357 m_y.clear();
1358 m_r.clear();
1359 m_p.clear();
1360
1361 //m_pressure.release();
1362
1363 if (m_reduce)
1364 {
1365 delete m_reduce;
1366 }
1367 if (m_arithmetic)
1368 {
1369 delete m_arithmetic;
1370 }
1371}
1372
1373template <typename TDataType>
1374void SemiAnalyticalIncompressibilityModule<TDataType>::constrain()
1375{
1376 //return true;
1377
1378 Real dt = getParentNode()->getDt();
1379
1380 // int start_f = Start.getValue();
1381 // cudaMemcpy(m_velocityAll.getValue().getDataPtr() + start_f, m_particle_velocity.getValue().getDataPtr(), num_f * sizeof(Coord), cudaMemcpyDeviceToDevice);
1382
1383 std::cout << "Element Count: " << m_particle_position.size() << std::endl;
1384
1385 uint pDims = cudaGridSize(m_particle_position.size(), BLOCK_SIZE);
1386
1387 //compute alpha_i = sigma w_j and A_i = sigma w_ij / r_ij / r_ij
1388
1389 printf("inside VC constraint NEW %d %d %d %d\n", m_particle_velocity.getData().size(), m_particle_mass.size(), m_particle_attribute.size(), m_particle_position.size());
1390
1391 int numTri = m_triangle_vertex.size();
1392 uint pDimsT = cudaGridSize(numTri, BLOCK_SIZE);
1393
1394 if (!m_particle_position.isEmpty())
1395 {
1396 printf("warning from second step!");
1397 int num = m_particle_position.size();
1398 uint pDims = cudaGridSize(num, BLOCK_SIZE);
1399
1400 //m_particle_attribute.resize(num);
1401 //m_particle_mass.resize(num);
1402
1403 // m_particle_normal.resize(num);
1404 //if (m_particle_attribute.isEmpty())printf("???\n");
1405 //m_particle_attribute.getReference()->reset();
1406 //VC_InitAttrTmp << <pDims, BLOCK_SIZE >> > (
1407 // m_particle_attribute.getValue(),
1408 // m_particle_mass.getValue()
1409 // );
1410 }
1411
1412 int num = m_particle_position.size();
1413 if (m_AiiFluid.isEmpty() || m_AiiFluid.size() != num)
1414 {
1415 //printf("RESIZW\n");
1416 m_alpha.resize(num);
1417 Rho_alpha.resize(num);
1418 m_Aii.resize(num);
1419 m_AiiFluid.resize(num);
1420 m_AiiTotal.resize(num);
1421 m_pressure.resize(num);
1422 m_divergence.resize(num);
1423 m_bSurface.resize(num);
1424 //m_density.resize(num);
1425
1426 m_y.resize(num);
1427 m_r.resize(num);
1428 m_p.resize(num);
1429
1430 m_flip.resize(num);
1431
1432 m_reduce = Reduction<float>::Create(num);
1433 m_arithmetic = Arithmetic<float>::Create(num);
1434 }
1435
1436 m_flip.getData().reset();
1437 m_alpha.reset();
1438 //printf("sampling_distance = %.10lf; smoothing_length = %.10lf\n", m_sampling_distance.getValue(), m_smoothing_length.getValue());
1439
1440 VC_TriVelTmp<<<pDimsT, BLOCK_SIZE>>>(
1441 m_triangle_vertex_old.getData(),
1442 m_triangle_vertex.getData(),
1443 m_meshVel,
1444 dt);
1445
1446// VC_Sort_Neighbors<<<pDims, BLOCK_SIZE>>>(
1447// m_particle_position.getData(),
1448// m_triangle_index.getData(),
1449// m_triangle_vertex.getData(),
1450// m_neighborhood_triangles.getData());
1451
1452 VC_ComputeAlphaTmp<<<pDims, BLOCK_SIZE>>>(
1453 m_alpha,
1454 Rho_alpha,
1455 m_particle_mass.getData(),
1456 m_particle_position.getData(),
1457 m_triangle_index.getData(),
1458 m_triangle_vertex.getData(),
1459 m_particle_attribute.getData(),
1460 this->inNeighborParticleIds()->getData(),
1461 this->inNeighborTriangleIds()->getData(),
1462 m_smoothing_length.getData(),
1463 m_sampling_distance.getData(),
1464 m_flip.getData());
1465 cuSynchronize();
1466
1467 //Real m_maxAlpha2 = m_reduce->maximum(m_alpha.getDataPtr(), m_alpha.size());
1468 //m_maxAlpha = max(m_maxAlpha2, m_maxAlpha);
1469
1470 VC_CorrectAlphaTmp<<<pDims, BLOCK_SIZE>>>(
1471 m_alpha,
1472 Rho_alpha,
1473 m_particle_mass.getData(),
1474 m_maxAlpha);
1475
1476 //compute the diagonal elements of the coefficient matrix
1477 m_AiiFluid.reset();
1478 m_AiiTotal.reset();
1479 VC_ComputeDiagonalElementTmp<<<pDims, BLOCK_SIZE>>>(
1480 m_AiiFluid,
1481 m_AiiTotal,
1482 m_alpha,
1483 m_particle_position.getData(),
1484 m_triangle_index.getData(),
1485 m_triangle_vertex.getData(),
1486 m_particle_attribute.getData(),
1487 this->inNeighborParticleIds()->getData(),
1488 this->inNeighborTriangleIds()->getData(),
1489 m_smoothing_length.getData(),
1490 m_sampling_distance.getData(),
1491 m_flip.getData());
1492
1493 m_bSurface.reset();
1494 m_Aii.reset();
1495 VC_DetectSurfaceTmp<<<pDims, BLOCK_SIZE>>>(
1496 m_Aii,
1497 m_bSurface,
1498 m_AiiFluid,
1499 m_AiiTotal,
1500 m_particle_position.getData(),
1501 m_particle_attribute.getData(),
1502 this->inNeighborParticleIds()->getData(),
1503 m_smoothing_length.getData(),
1504 m_maxA);
1505
1506 int itor = 0;
1507
1508 //compute the source term
1509 //m_densitySum->compute(m_density);
1510 m_densitySum->compute();
1511 //m_density = m_densitySum->outDensity()->getValue();
1512
1513 m_divergence.reset();
1514 VC_ComputeDivergenceTmp<<<pDims, BLOCK_SIZE>>>(
1515 m_divergence,
1516 m_alpha,
1517 m_densitySum->outDensity()->getData(),
1518 m_particle_position.getData(),
1519 m_particle_velocity.getData(),
1520 m_meshVel,
1521 m_triangle_index.getData(),
1522 m_triangle_vertex.getData(),
1523 m_bSurface,
1524 m_particle_attribute.getData(),
1525 m_particle_mass.getData(),
1526 m_triangle_vertex_mass.getData(),
1527 this->inNeighborParticleIds()->getData(),
1528 this->inNeighborTriangleIds()->getData(),
1529 m_separation,
1530 m_tangential,
1531 m_restDensity,
1532 m_smoothing_length.getData(),
1533 m_sampling_distance.getData(),
1534 dt,
1535 m_flip.getData());
1536
1537 VC_CompensateSourceTmp<<<pDims, BLOCK_SIZE>>>( // no need
1538 m_divergence,
1539 m_densitySum->outDensity()->getData(),
1540 m_particle_attribute.getData(),
1541 m_particle_position.getData(),
1542 m_restDensity,
1543 dt);
1544
1545 //solve the linear system of equations with a conjugate gradient method.
1546 m_y.reset();
1547 m_pressure.reset();
1548 VC_ComputeAxTmp<<<pDims, BLOCK_SIZE>>>(
1549 m_y,
1550 m_pressure,
1551 m_Aii,
1552 m_alpha,
1553 m_particle_position.getData(),
1554 m_particle_attribute.getData(),
1555 this->inNeighborParticleIds()->getData(),
1556 m_smoothing_length.getData());
1557
1558 m_r.reset();
1559 Function2Pt::subtract(m_r, m_divergence, m_y);
1560 m_p.assign(m_r);
1561 Real rr = m_arithmetic->Dot(m_r, m_r);
1562 Real err = sqrt(rr / m_r.size());
1563
1564 printf("ERROR: %.10lf\n", err);
1565
1566 while (itor < 1000 && err > 1.0f)
1567 {
1568 m_y.reset();
1569 //VC_ComputeAx << <pDims, BLOCK_SIZE >> > (*yArr, *pArr, *aiiArr, *alphaArr, *posArr, *attArr, *neighborArr);
1570 VC_ComputeAxTmp<<<pDims, BLOCK_SIZE>>>(
1571 m_y,
1572 m_p,
1573 m_Aii,
1574 m_alpha,
1575 m_particle_position.getData(),
1576 m_particle_attribute.getData(),
1577 this->inNeighborParticleIds()->getData(),
1578 m_smoothing_length.getData());
1579
1580 float alpha = rr / m_arithmetic->Dot(m_p, m_y);
1581 Function2Pt::saxpy(m_pressure, m_p, m_pressure, alpha);
1582 Function2Pt::saxpy(m_r, m_y, m_r, -alpha);
1583
1584 Real rr_old = rr;
1585
1586 rr = m_arithmetic->Dot(m_r, m_r);
1587
1588 Real beta = rr / rr_old;
1589 Function2Pt::saxpy(m_p, m_p, m_r, beta);
1590
1591 err = sqrt(rr / m_r.size());
1592 printf("err: %.3lf\n:", err);
1593 itor++;
1594 }
1595 //return true;
1596 //update the each particle's velocity
1597 VC_UpdateVelocityBoundaryCorrectedTmp<<<pDims, BLOCK_SIZE>>>(
1598 m_pressure,
1599 m_alpha,
1600 m_bSurface,
1601 m_particle_position.getData(),
1602 m_particle_velocity.getData(),
1603 m_meshVel,
1604 m_triangle_index.getData(),
1605 m_triangle_vertex.getData(),
1606 m_particle_attribute.getData(),
1607 m_particle_mass.getData(),
1608 m_triangle_vertex_mass.getData(),
1609 this->inNeighborParticleIds()->getData(),
1610 this->inNeighborTriangleIds()->getData(),
1611 m_restDensity,
1612 m_airPressure,
1613 m_tangential,
1614 m_separation,
1615 m_smoothing_length.getData(),
1616 m_sampling_distance.getData(),
1617 dt,
1618 m_flip.getData());
1619
1620 //printf("pressure size: %d\n", m_pressure.size());
1621 //cudaMemcpy(m_pressure2.getValue().getDataPtr(), m_pressure.getDataPtr() + start_f, num_f * sizeof(Real), cudaMemcpyDeviceToDevice);
1622 // cudaMemcpy(m_particle_velocity.getValue().getDataPtr(), m_velocityAll.getValue().getDataPtr() + start_f, num_f * sizeof(Coord), cudaMemcpyDeviceToDevice);
1623 // cudaMemcpy(PressureFluid.getValue().getDataPtr(), m_pressure.getDataPtr() + start_f, num_f * sizeof(Real), cudaMemcpyDeviceToDevice);
1624 //cudaMemcpy(PressureFluid.getValue().getDataPtr(), m_divergence_Tri.getDataPtr() + start_f, num_f * sizeof(Real), cudaMemcpyDeviceToDevice);
1625}
1626
1627template <typename TDataType>
1628bool SemiAnalyticalIncompressibilityModule<TDataType>::initializeImpl()
1629{
1630 return false;
1631
1632 first_step = true;
1633
1634 if (m_particle_position.isEmpty())
1635 printf("OKKKKKKKKKK\n");
1636 m_sampling_distance.setValue(0.005);
1637
1638 // printf("%d %d %d %d\n", m_position.size(), m_position1.size(), m_velocity.size(), m_velocity2.size());
1639
1640 if (!m_particle_position.isEmpty())
1641 {
1642 printf("warning from second step!");
1643 int num = m_particle_position.size();
1644 uint pDims = cudaGridSize(num, BLOCK_SIZE);
1645
1646 m_particle_attribute.resize(num);
1647 m_particle_mass.resize(num);
1648
1649 // m_particle_normal.resize(num);
1650 //if (m_particle_attribute.isEmpty())printf("???\n");
1651 m_particle_attribute.getDataPtr()->reset();
1652 VC_InitAttrTmp<<<pDims, BLOCK_SIZE>>>(
1653 m_particle_attribute.getData(),
1654 m_particle_mass.getData());
1655 }
1656 else
1657 {
1658 printf("YES~ m_triangle_index Size: %d\n", m_triangle_index.size());
1659 }
1660
1661 //TODO: input the time step size
1662 Real dt = 0.001f;
1663 int numt = m_triangle_vertex.size();
1664
1665 m_meshVel.resize(numt);
1666 uint pDims = cudaGridSize(numt, BLOCK_SIZE);
1667
1668 VC_TriVelTmp<<<pDims, BLOCK_SIZE>>>(
1669 m_triangle_vertex_old.getData(),
1670 m_triangle_vertex.getData(),
1671 m_meshVel,
1672 dt);
1673
1674 //m_neighborhood->initialize();
1675 /*
1676 m_densitySum = std::make_shared<DensitySummation<TDataType>>();
1677 m_smoothing_length.connect(&m_densitySum->m_smoothingLength);
1678 m_particle_position.connect(&m_densitySum->m_position);
1679 m_neighborhood_particles.connect(&m_densitySum->m_neighborhood);
1680 m_densitySum->initialize();
1681 */
1682
1683 m_densitySum = std::make_shared<SummationDensity<TDataType>>();
1684 m_smoothing_length.connect(m_densitySum->inSmoothingLength());
1685 m_particle_position.connect(m_densitySum->inPosition());
1686 this->inNeighborParticleIds()->connect(m_densitySum->inNeighborIds());
1687 m_densitySum->initialize();
1688
1689 int num = m_particle_position.size();
1690
1691 if (num == 0)
1692 {
1693 m_maxAlpha = 36.8;
1694 m_maxA = 35205.6;
1695 return true;
1696 }
1697
1698 num_f = m_particle_position.size();
1699 //int start_f = 11286;//5130;
1700
1701 m_alpha.resize(num);
1702 Rho_alpha.resize(num);
1703 m_Aii.resize(num);
1704 m_AiiFluid.resize(num);
1705 m_AiiTotal.resize(num);
1706 m_pressure.resize(num);
1707 m_divergence.resize(num);
1708 //m_divergence_Tri.resize(num);
1709 m_bSurface.resize(num);
1710 m_density.resize(num);
1711
1712 m_y.resize(num);
1713 m_r.resize(num);
1714 m_p.resize(num);
1715
1716 // m_pressure.resize(num);
1717
1718 m_reduce = Reduction<float>::Create(num);
1719 m_arithmetic = Arithmetic<float>::Create(num);
1720
1721 pDims = cudaGridSize(num, BLOCK_SIZE);
1722 /*
1723 VC_Sort_Neighbors << <pDims, BLOCK_SIZE >> > (
1724 m_particle_position.getValue(),
1725 m_triangle_index.getValue(),
1726 m_triangle_vertex.getValue(),
1727 m_neighborhood_triangles.getValue()
1728 );
1729
1730 VC_Calc_SolidAngle << <pDims, BLOCK_SIZE >> > (
1731 m_particle_position.getValue(),
1732 m_triangle_index.getValue(),
1733 m_triangle_vertex.getValue(),
1734 m_neighborhood_triangles.getValue(),
1735 m_smoothing_length.getValue()
1736 );
1737*/
1738
1739 //printf("TRI:%d\n", m_triangle_index.getValue().ge);
1740 // printf("NEI1:%d\n", m_neighborhood.isEmpty());
1741 m_alpha.reset();
1742 printf("FLIP: %d\n", m_flip.getData().size());
1743 VC_ComputeAlphaTmp<<<pDims, BLOCK_SIZE>>>(
1744 m_alpha,
1745 Rho_alpha,
1746 m_particle_mass.getData(),
1747 m_particle_position.getData(),
1748 m_triangle_index.getData(),
1749 m_triangle_vertex.getData(),
1750 m_particle_attribute.getData(),
1751 this->inNeighborParticleIds()->getData(),
1752 this->inNeighborTriangleIds()->getData(),
1753 m_smoothing_length.getData(),
1754 m_sampling_distance.getData(),
1755 m_flip.getData());
1756
1757 m_maxAlpha = m_reduce->maximum(m_alpha.begin(), m_alpha.size());
1758
1759 VC_CorrectAlphaTmp<<<pDims, BLOCK_SIZE>>>(
1760 m_alpha,
1761 Rho_alpha,
1762 m_particle_mass.getData(),
1763 m_maxAlpha);
1764
1765 m_AiiFluid.reset();
1766 VC_ComputeDiagonalElementTmp<<<pDims, BLOCK_SIZE>>>(
1767 m_AiiFluid,
1768 m_alpha,
1769 m_particle_position.getData(),
1770 m_particle_attribute.getData(),
1771 this->inNeighborParticleIds()->getData(),
1772 m_smoothing_length.getData());
1773
1774 m_maxA = m_reduce->maximum(m_AiiFluid.begin(), m_AiiFluid.size());
1775
1776 std::cout << "Max alpha: " << m_maxAlpha << std::endl;
1777 printf("%.10lf\n", m_maxAlpha);
1778 std::cout << "Max A: " << m_maxA << std::endl;
1779
1780 return true;
1781}
1782
1783DEFINE_CLASS(SemiAnalyticalIncompressibilityModule)
1784
1785}
Generated by doxygen 1.13.2