PeriDyno 1.0.0
Loading...
Searching...
No Matches
SemiAnalyticalPBD.cu
Go to the documentation of this file.
1#include "SemiAnalyticalPBD.h"
2
3#include "SemiAnalyticalSummationDensity.h"
4#include "IntersectionArea.h"
5
6namespace dyno
7{
8 IMPLEMENT_TCLASS(SemiAnalyticalPBD, TDataType)
9
10 template <typename TDataType>
11 SemiAnalyticalPBD<TDataType>::SemiAnalyticalPBD()
12 : ConstraintModule()
13 {
14 mCalculateDensity = std::make_shared<SemiAnalyticalSummationDensity<TDataType>>();
15 this->inSmoothingLength()->connect(mCalculateDensity->inSmoothingLength());
16 this->inSamplingDistance()->connect(mCalculateDensity->inSamplingDistance());
17 this->inPosition()->connect(mCalculateDensity->inPosition());
18 this->inNeighborParticleIds()->connect(mCalculateDensity->inNeighborIds());
19 this->inNeighborTriangleIds()->connect(mCalculateDensity->inNeighborTriIds());
20 this->inTriangleSet()->connect(mCalculateDensity->inTriangleSet());
21 }
22
23 template <typename TDataType>
24 SemiAnalyticalPBD<TDataType>::~SemiAnalyticalPBD()
25 {
26 mLamda.clear();
27 mDeltaPos.clear();
28 mPosBuffer.clear();
29 }
30
31 template <typename TDataType>
32 void SemiAnalyticalPBD<TDataType>::constrain()
33 {
34 int num = this->inPosition()->size();
35
36 if (mLamda.size() != num) {
37 mLamda.resize(num);
38 mDeltaPos.resize(num);
39 mPosBuffer.resize(num);
40 }
41
42 mPosBuffer.assign(this->inPosition()->getData());
43
44 int iter = 0;
45 auto maxIter = this->varInterationNumber()->getData();
46 while (iter++ < maxIter) {
47 takeOneIteration();
48 }
49
50 updateVelocity();
51 }
52
53 __device__ inline float kernGradientMeshPBD(const float r, const float h)
54 {
55 const Real q = r / h;
56 if (q > 1.0f)
57 return 0.0;
58 else
59 {
60 //G(r) in equation 6
61 const Real d = 1.0 - q;
62 const Real hh = h * h;
63 return -45.0f / ((Real)M_PI * hh * h) * (1.0f / 3.0f * (hh * h - r * r * r) - 1.0f / 2.0f / h * (hh * hh - r * r * r * r) + 1.0f / 5.0f / hh * (hh * hh * h - r * r * r * r * r));
64 }
65 }
66
67 template <typename Real, typename Coord>
68 __global__ void K_ComputeLambdasMesh(
69 DArray<Real> lambdaArr,
70 DArray<Real> rhoArr,
71 DArray<Coord> posArr,
72 DArray<TopologyModule::Triangle> Tri,
73 DArray<Coord> positionTri,
74 DArrayList<int> neighbors,
75 DArrayList<int> neighborsTri,
76 SpikyKernel<Real> kern,
77 Real smoothingLength,
78 Real samplingDistance)
79 {
80 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
81 if (pId >= posArr.size())
82 return;
83
84 List<int>& triList_i = neighborsTri[pId];
85 int nbSizeTri = triList_i.size();
86 Real dis_n = REAL_MAX;
87
88 //semi-analytical boundary integration
89 Real lamda_i = Real(0);
90 Coord pos_i = posArr[pId];
91 Coord grad_ci(0);
92 for (int ne = 0; ne < nbSizeTri; ne++)
93 {
94 int j = triList_i[ne];
95 Triangle3D t3d(positionTri[Tri[j][0]], positionTri[Tri[j][1]], positionTri[Tri[j][2]]);
96 Plane3D PL(positionTri[Tri[j][0]], t3d.normal());
97 Point3D p3d(pos_i);
98 Point3D nearest_pt = p3d.project(PL);
99 Real r = (nearest_pt.origin - pos_i).norm();
100
101 Real AreaSum = calculateIntersectionArea(p3d, t3d, smoothingLength); //A_s in equation 10
102 Real MinDistance = (p3d.distance(t3d)); //d_n (scalar) in equation 10
103 Coord Min_Pt = (p3d.project(t3d)).origin - pos_i; //d_n (vector) in equation 10
104 Coord Min_Pos = p3d.project(t3d).origin;
105 if (ne < nbSizeTri - 1 && triList_i[ne + 1] < 0)
106 {
107 //triangle clustering
108 int jn;
109 do
110 {
111 jn = triList_i[ne + 1];
112
113 Triangle3D t3d_n(positionTri[Tri[jn][0]], positionTri[Tri[jn][1]], positionTri[Tri[jn][2]]);
114 if ((t3d.normal().cross(t3d_n.normal())).norm() > EPSILON)
115 break;
116
117 AreaSum += calculateIntersectionArea(p3d, t3d_n, smoothingLength);
118
119 if (abs(p3d.distance(t3d_n)) < abs(MinDistance))
120 {
121 MinDistance = (p3d.distance(t3d_n));
122 Min_Pt = (p3d.project(t3d_n)).origin - pos_i;
123 Min_Pos = (p3d.project(t3d_n)).origin;
124 }
125 //printf("%d %d\n", j, jn);
126 ne++;
127 } while (ne < nbSizeTri - 1);
128 }
129 if (abs(MinDistance) < abs(dis_n))
130 dis_n = MinDistance;
131 Min_Pt /= (-Min_Pt.norm());
132
133 float d = p3d.distance(PL);
134 d = abs(d);
135
136 // equation 6
137 if (smoothingLength - d > EPSILON && smoothingLength * smoothingLength - d * d > EPSILON && d > EPSILON)
138 {
139
140 Real a_ij =
141 kernGradientMeshPBD(r, smoothingLength)
142 / (samplingDistance * samplingDistance * samplingDistance)
143 * 2.0 * (M_PI) * (1 - d / smoothingLength) //eq 11
144 * AreaSum //p3d.areaTriangle(t3d, smoothingLength)
145 / ((M_PI) * (smoothingLength * smoothingLength - d * d)) //eq 11
146 * t3d.normal().dot(Min_Pt) / t3d.normal().norm();
147
148 {
149 Coord g = a_ij * (pos_i - nearest_pt.origin) / r;
150 grad_ci += g;
151 lamda_i += g.dot(g);
152 }
153 }
154 }
155 grad_ci *= (dis_n / abs(dis_n));
156 lamda_i *= (dis_n / abs(dis_n));
157
158 //traditional integration position based fluids
159 List<int>& parList_i = neighbors[pId];
160 int nbSize = parList_i.size();
161 for (int ne = 0; ne < nbSize; ne++)
162 {
163 int j = parList_i[ne];
164
165 Real r = (pos_i - posArr[j]).norm();
166 if (r > EPSILON)
167 {
168 Coord g = kern.Gradient(r, smoothingLength) * (pos_i - posArr[j]) * (1.0f / r);
169 grad_ci += g;
170 lamda_i += g.dot(g);
171 }
172 }
173
174 lamda_i += grad_ci.dot(grad_ci);
175
176 Real rho_i = rhoArr[pId];
177
178 lamda_i = -(rho_i - 1000.0f) / (lamda_i + 0.1f);
179
180 lambdaArr[pId] = lamda_i > 0.0f ? 0.0f : lamda_i;
181 }
182
183 template <typename Real, typename Coord>
184 __global__ void K_ComputeDisplacementMesh(
185 DArray<Coord> dPos,
186 DArray<Real> lambdas,
187 DArray<Coord> posArr,
188 DArray<TopologyModule::Triangle> Tri,
189 DArray<Coord> positionTri,
190 DArrayList<int> neighbors,
191 DArrayList<int> neighborsTri,
192 SpikyKernel<Real> kern,
193 Real smoothingLength,
194 Real dt,
195 Real sampling_distance)
196 {
197 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
198 if (pId >= posArr.size())
199 return;
200
201 Coord pos_i = posArr[pId];
202 Real lamda_i = lambdas[pId];
203
204 Coord dP_i(0);
205 List<int>& list_i = neighbors[pId];
206 int nbSize = list_i.size();
207
208 List<int>& triList_i = neighborsTri[pId];
209 int nbSizeTri = triList_i.size();
210 Real dis_n = REAL_MAX;
211 for (int ne = 0; ne < nbSizeTri; ne++)
212 {
213 int j = triList_i[ne];
214 Triangle3D t3d(positionTri[Tri[j][0]], positionTri[Tri[j][1]], positionTri[Tri[j][2]]);
215 Plane3D PL(positionTri[Tri[j][0]], t3d.normal());
216 Point3D p3d(pos_i);
217 if (abs(p3d.distance(t3d)) < abs(dis_n))
218 {
219 dis_n = p3d.distance(t3d);
220 }
221 }
222 //semi-analytical boundary integration
223 for (int ne = 0; ne < nbSizeTri; ne++)
224 {
225 int j = triList_i[ne];
226 Triangle3D t3d(positionTri[Tri[j][0]], positionTri[Tri[j][1]], positionTri[Tri[j][2]]);
227 Plane3D PL(positionTri[Tri[j][0]], t3d.normal());
228 Point3D p3d(pos_i);
229 //Point3D nearest_pt = p3d.project(t3d);
230 Point3D nearest_pt = p3d.project(PL);
231 Real r = (nearest_pt.origin - pos_i).norm();
232
233 Real tmp = 1.0;
234 float d = p3d.distance(PL);
235
236 Coord ttm = PL.normal;
237 Real AreaSum = calculateIntersectionArea(p3d, t3d, smoothingLength); //A_s in equation 10
238 Real MinDistance = abs(p3d.distance(t3d)); //d_n (scalar) in equation 10
239 Coord Min_Pt = (p3d.project(t3d)).origin - pos_i; //d_n (vector) in equation 10
240 Coord Min_Pos = p3d.project(t3d).origin;
241 if (ne < nbSizeTri - 1 && triList_i[ne + 1] < 0)
242 {
243 //triangle clustering
244 int jn;
245 do
246 {
247 jn = triList_i[ne + 1];
248 Triangle3D t3d_n(positionTri[Tri[jn][0]], positionTri[Tri[jn][1]], positionTri[Tri[jn][2]]);
249 if ((t3d.normal().cross(t3d_n.normal())).norm() > EPSILON)
250 break;
251
252 AreaSum += calculateIntersectionArea(p3d, t3d_n, smoothingLength);
253
254 if (abs(p3d.distance(t3d_n)) < abs(MinDistance))
255 {
256 MinDistance = (p3d.distance(t3d_n));
257 Min_Pt = (p3d.project(t3d_n)).origin - pos_i;
258 Min_Pos = (p3d.project(t3d_n)).origin;
259 }
260 //printf("%d %d\n", j, jn);
261 ne++;
262 } while (ne < nbSizeTri - 1);
263 }
264
265 Min_Pt /= (-Min_Pt.norm());
266
267 d = abs(d);
268 //r = max((r - sampling_distance / 2.0), 0.0);
269 if (smoothingLength - d > EPSILON && smoothingLength * smoothingLength - d * d > EPSILON && d > EPSILON)
270 {
271 //equaltion 6
272 Real a_ij =
273 kernGradientMeshPBD(r, smoothingLength)
274 * 2.0 * (M_PI) * (1 - d / smoothingLength) //eq 11
275 * AreaSum //p3d.areaTriangle(t3d, smoothingLength)
276 / ((M_PI) * (smoothingLength * smoothingLength - d * d)) // eq11
277 * t3d.normal().dot(Min_Pt) / t3d.normal().norm() // / (p3d.project(t3d).origin - p3d.origin).norm()
278 / (sampling_distance * sampling_distance * sampling_distance);
279 //a_ij *= (dis_n / abs(dis_n));
280
281 Coord dp_ij = 40.0f * (pos_i - nearest_pt.origin) * (lamda_i)*a_ij * (1.0 / (pos_i - nearest_pt.origin).norm());
282
283 //if (a_ij < 0)
284 {
285 dp_ij *= tmp;
286 dP_i += dp_ij;
287 atomicAdd(&dPos[pId][0], dp_ij[0]);
288
289 if (Coord::dims() >= 2)
290 atomicAdd(&dPos[pId][1], dp_ij[1]);
291
292 if (Coord::dims() >= 3)
293 atomicAdd(&dPos[pId][2], dp_ij[2]);
294 }
295 }
296 }
297
298 //traditional integration position based fluids
299 for (int ne = 0; ne < nbSize; ne++)
300 {
301 int j = list_i[ne];
302 Real r = (pos_i - posArr[j]).norm();
303 if (r > EPSILON)
304 {
305 Coord dp_ij = 10.0f * (pos_i - posArr[j]) * (lamda_i + lambdas[j]) * kern.Gradient(r, smoothingLength) * (1.0 / r);
306 dP_i += dp_ij;
307
308 atomicAdd(&dPos[pId][0], dp_ij[0]);
309 atomicAdd(&dPos[j][0], -dp_ij[0]);
310
311 if (Coord::dims() >= 2)
312 {
313 atomicAdd(&dPos[pId][1], dp_ij[1]);
314 atomicAdd(&dPos[j][1], -dp_ij[1]);
315 }
316
317 if (Coord::dims() >= 3)
318 {
319 atomicAdd(&dPos[pId][2], dp_ij[2]);
320 atomicAdd(&dPos[j][2], -dp_ij[2]);
321 }
322 }
323 }
324 }
325
326 template <typename Coord>
327 __global__ void K_UpdatePositionMesh(
328 DArray<Coord> posArr,
329 DArray<Coord> dPos)
330 {
331 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
332 if (pId >= posArr.size())
333 return;
334
335 posArr[pId] += dPos[pId];
336 }
337
338 template <typename Real, typename Coord>
339 __global__ void DP_UpdateVelocityMesh(
340 DArray<Coord> velArr,
341 DArray<Coord> prePos,
342 DArray<Coord> curPos,
343 Real dt)
344 {
345 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
346 if (pId >= velArr.size())
347 return;
348
349 velArr[pId] += (curPos[pId] - prePos[pId]) / dt;
350 }
351
352 template <typename TDataType>
353 void SemiAnalyticalPBD<TDataType>::takeOneIteration()
354 {
355 auto& inPos = this->inPosition()->getData();
356
357 auto ts = this->inTriangleSet()->constDataPtr();
358 auto& triVertex = ts->getPoints();
359 auto& triIndex = ts->getTriangles();
360
361 mDeltaPos.reset();
362 mCalculateDensity->update();
363
364 int num = inPos.size();
365
366 cuExecute(num,
367 K_ComputeLambdasMesh,
368 mLamda,
369 mCalculateDensity->outDensity()->getData(),
370 inPos,
371 triIndex,
372 triVertex,
373 this->inNeighborParticleIds()->getData(),
374 this->inNeighborTriangleIds()->getData(),
375 m_kernel,
376 this->inSmoothingLength()->getData(),
377 this->inSamplingDistance()->getData());
378
379 cuExecute(num,
380 K_ComputeDisplacementMesh,
381 mDeltaPos,
382 mLamda,
383 inPos,
384 triIndex,
385 triVertex,
386 this->inNeighborParticleIds()->getData(),
387 this->inNeighborTriangleIds()->getData(),
388 m_kernel,
389 this->inSmoothingLength()->getData(),
390 this->inTimeStep()->getData(),
391 this->inSamplingDistance()->getData());
392
393 cuExecute(num,
394 K_UpdatePositionMesh,
395 inPos,
396 mDeltaPos);
397 }
398
399 template <typename TDataType>
400 void SemiAnalyticalPBD<TDataType>::updateVelocity()
401 {
402 int num = this->inPosition()->size();
403 cuExecute(num,
404 DP_UpdateVelocityMesh,
405 this->inVelocity()->getData(),
406 mPosBuffer,
407 this->inPosition()->getData(),
408 this->inTimeStep()->getData());
409 }
410
411 DEFINE_CLASS(SemiAnalyticalPBD);
412}
Generated by doxygen 1.13.2