PeriDyno 1.0.0
Loading...
Searching...
No Matches
SemiImplicitHyperelasticitySolver.cu
Go to the documentation of this file.
1#include "SemiImplicitHyperelasticitySolver.h"
2
3#include "Matrix/MatrixFunc.h"
4#include "ParticleSystem/Module/Kernel.h"
5#include "curand_kernel.h"
6#include "Algorithm/CudaRand.h"
7
8namespace dyno
9{
10 IMPLEMENT_TCLASS(SemiImplicitHyperelasticitySolver, TDataType);
11
12 __constant__ EnergyModels<Real> ENERGY_FUNC;
13
14 template<typename Real>
15 __device__ Real constantWeight(Real r, Real h)
16 {
17 Real d = h / r;
18 return d * d;
19 }
20
21 template<typename Real>
22 __device__ Real D_Weight(Real r, Real h)
23 {
24 CorrectedKernel<Real> kernSmooth;
25 return kernSmooth.WeightRR(r, 4 * h);
26 }
27
28 template<typename TDataType>
29 SemiImplicitHyperelasticitySolver<TDataType>::SemiImplicitHyperelasticitySolver()
30 : LinearElasticitySolver<TDataType>()
31 {
32 this->varIterationNumber()->setValue(5);
33 }
34
35 template<typename TDataType>
36 SemiImplicitHyperelasticitySolver<TDataType>::~SemiImplicitHyperelasticitySolver()
37 {
38 this->mWeights.clear();
39 this->mDisplacement.clear();
40 this->mInvK.clear();
41 this->mF.clear();
42 this->mPosBuf.clear();
43 }
44
45 template <typename Real, typename Coord, typename Matrix>
46 __global__ void HM_ComputeEnergy(
47 DArray<Real> energy,
48 DArray<Coord> eigens,
49 EnergyType type)
50 {
51 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
52 if (pId >= energy.size()) return;
53
54 Coord eigen_i = eigens[pId];
55
56 if (type == StVK) {
57 energy[pId] = ENERGY_FUNC.stvkModel.getEnergy(eigen_i[0], eigen_i[1], eigen_i[2]);
58 }
59 else if (type == NeoHooekean) {
60 energy[pId] = ENERGY_FUNC.neohookeanModel.getEnergy(eigen_i[0], eigen_i[1], eigen_i[2]);
61 }
62 else if (type == Linear) {
63 energy[pId] = ENERGY_FUNC.linearModel.getEnergy(eigen_i[0], eigen_i[1], eigen_i[2]);
64 }
65 else if (type == Xuetal) {
66 energy[pId] = ENERGY_FUNC.xuModel.getEnergy(eigen_i[0], eigen_i[1], eigen_i[2]);
67 }
68 else if (type == MooneyRivlin) {
69 energy[pId] = ENERGY_FUNC.mrModel.getEnergy(eigen_i[0], eigen_i[1], eigen_i[2]);
70 }
71 else if (type == Fung) {
72 energy[pId] = ENERGY_FUNC.fungModel.getEnergy(eigen_i[0], eigen_i[1], eigen_i[2]);
73 }
74 else if (type == Ogden) {
75 energy[pId] = ENERGY_FUNC.ogdenModel.getEnergy(eigen_i[0], eigen_i[1], eigen_i[2]);
76 }
77 else if (type == Yeoh) {
78 energy[pId] = ENERGY_FUNC.yeohModel.getEnergy(eigen_i[0], eigen_i[1], eigen_i[2]);
79 }
80 else if (type == ArrudaBoyce) {
81 energy[pId] = ENERGY_FUNC.abModel.getEnergy(eigen_i[0], eigen_i[1], eigen_i[2]);
82 }
83 }
84
85
86 template<typename TDataType>
87 void SemiImplicitHyperelasticitySolver<TDataType>::solveElasticity()
88 {
89 cudaMemcpyToSymbol(ENERGY_FUNC, this->inEnergyModels()->constDataPtr().get(), sizeof(EnergyModels<Real>));
90
91 enforceHyperelasticity();
92 }
93
94 template <typename Coord>
95 __global__ void HM_ComputeGradient(
96 DArray<Coord> grad,
97 DArray<Coord> y_pre,
98 DArray<Coord> y_next)
99 {
100 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
101 if (pId >= y_next.size()) return;
102
103 grad[pId] = y_next[pId] - y_pre[pId];
104 }
105
106 template <typename Real, typename Coord>
107 __global__ void HM_ComputeCurrentPosition(
108 DArray<Coord> grad,
109 DArray<Coord> y_current,
110 DArray<Coord> y_next,
111 DArray<Attribute> attribute,
112 DArray<Real> alpha)
113 {
114 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
115 if (pId >= y_next.size()) return;
116
117 y_next[pId] = attribute[pId].isDynamic() ? y_current[pId] + alpha[pId] * grad[pId] : y_current[pId];
118 }
119
120 template <typename Real, typename Coord>
121 __global__ void HM_ComputeCurrentPosition(
122 DArray<Coord> grad,
123 DArray<Coord> y_current,
124 DArray<Coord> y_next,
125 Real alpha)
126 {
127 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
128 if (pId >= y_next.size()) return;
129
130 y_next[pId] = y_current[pId] + alpha * grad[pId];
131 }
132
133 template <typename Real, typename Coord, typename Matrix, typename Bond>
134 __global__ void HM_Compute1DEnergy(
135 DArray<Real> energy,
136 DArray<Coord> energyGradient,
137 DArray<Coord> X,
138 DArray<Coord> Y,
139 DArray<Matrix> F,
140 DArray<Real> volume,
141 DArray<bool> validOfK,
142 DArray<Coord> eigenValues,
143 DArrayList<Bond> bonds,
144 EnergyType type)
145 {
146 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
147 if (pId >= energy.size()) return;
148
149 Coord y_i = Y[pId];
150
151 Real totalEnergy = 0.0f;
152 Coord totalEnergyGradient = Coord(0);
153 Real V_i = volume[pId];
154
155 int size_i = bonds[pId].size();
156
157 Coord x_i = X[pId];
158 Coord eigen_value_i = eigenValues[pId];
159 bool valid_i = validOfK[pId];
160
161 Matrix F_i = F[pId];
162
163 for (int ne = 1; ne < size_i; ne++)
164 {
165 Bond bond_ij = bonds[pId][ne];
166 int j = bond_ij.idx;
167 Coord y_j = Y[j];
168 Real r = bond_ij.xi.norm();
169
170 Real V_j = volume[j];
171
172 if (r > EPSILON)
173 {
174 Real norm_ij = (y_j - y_i).norm();
175 Real lambda_ij = norm_ij / r;
176
177 Real deltaEnergy;
178 Coord deltaEnergyGradient;
179 Coord dir_ij = norm_ij < EPSILON ? Coord(0) : (y_i - y_j) / (r);
180
181 if (type == StVK) {
182 deltaEnergy = V_j * ENERGY_FUNC.stvkModel.getEnergy(lambda_ij, lambda_ij, lambda_ij);
183 deltaEnergyGradient = V_j * (ENERGY_FUNC.stvkModel.getStressTensorPositive(lambda_ij, lambda_ij, lambda_ij) - ENERGY_FUNC.stvkModel.getStressTensorNegative(lambda_ij, lambda_ij, lambda_ij)) * dir_ij;
184 }
185 else if (type == NeoHooekean) {
186 deltaEnergy = V_j * ENERGY_FUNC.neohookeanModel.getEnergy(lambda_ij, lambda_ij, lambda_ij);
187 deltaEnergyGradient = V_j * (ENERGY_FUNC.neohookeanModel.getStressTensorPositive(lambda_ij, lambda_ij, lambda_ij) - ENERGY_FUNC.neohookeanModel.getStressTensorNegative(lambda_ij, lambda_ij, lambda_ij)) * dir_ij;
188 }
189 else if (type == Linear) {
190 deltaEnergy = V_j * ENERGY_FUNC.linearModel.getEnergy(lambda_ij, lambda_ij, lambda_ij);
191 deltaEnergyGradient = V_j * (ENERGY_FUNC.linearModel.getStressTensorPositive(lambda_ij, lambda_ij, lambda_ij) - ENERGY_FUNC.linearModel.getStressTensorNegative(lambda_ij, lambda_ij, lambda_ij)) * dir_ij;
192 }
193 else if (type == Xuetal) {
194 deltaEnergy = V_j * ENERGY_FUNC.xuModel.getEnergy(lambda_ij, lambda_ij, lambda_ij);
195 deltaEnergyGradient = V_j * (ENERGY_FUNC.xuModel.getStressTensorPositive(lambda_ij, lambda_ij, lambda_ij) - ENERGY_FUNC.xuModel.getStressTensorNegative(lambda_ij, lambda_ij, lambda_ij)) * dir_ij;
196 }
197 else if (type == MooneyRivlin) {
198 deltaEnergy = V_j * ENERGY_FUNC.mrModel.getEnergy(lambda_ij, lambda_ij, lambda_ij);
199 deltaEnergyGradient = V_j * (ENERGY_FUNC.mrModel.getStressTensorPositive(lambda_ij, lambda_ij, lambda_ij) - ENERGY_FUNC.xuModel.getStressTensorNegative(lambda_ij, lambda_ij, lambda_ij)) * dir_ij;
200 }
201 else if (type == Fung) {
202 deltaEnergy = V_j * ENERGY_FUNC.fungModel.getEnergy(lambda_ij, lambda_ij, lambda_ij);
203 deltaEnergyGradient = V_j * (ENERGY_FUNC.fungModel.getStressTensorPositive(lambda_ij, lambda_ij, lambda_ij) - ENERGY_FUNC.xuModel.getStressTensorNegative(lambda_ij, lambda_ij, lambda_ij)) * dir_ij;
204 }
205 else if (type == Ogden) {
206 deltaEnergy = V_j * ENERGY_FUNC.ogdenModel.getEnergy(lambda_ij, lambda_ij, lambda_ij);
207 deltaEnergyGradient = V_j * (ENERGY_FUNC.ogdenModel.getStressTensorPositive(lambda_ij, lambda_ij, lambda_ij) - ENERGY_FUNC.xuModel.getStressTensorNegative(lambda_ij, lambda_ij, lambda_ij)) * dir_ij;
208 }
209 else if (type == Yeoh) {
210 deltaEnergy = V_j * ENERGY_FUNC.yeohModel.getEnergy(lambda_ij, lambda_ij, lambda_ij);
211 deltaEnergyGradient = V_j * (ENERGY_FUNC.yeohModel.getStressTensorPositive(lambda_ij, lambda_ij, lambda_ij) - ENERGY_FUNC.xuModel.getStressTensorNegative(lambda_ij, lambda_ij, lambda_ij)) * dir_ij;
212 }
213 else if (type == ArrudaBoyce) {
214 deltaEnergy = V_j * ENERGY_FUNC.abModel.getEnergy(lambda_ij, lambda_ij, lambda_ij);
215 deltaEnergyGradient = V_j * (ENERGY_FUNC.abModel.getStressTensorPositive(lambda_ij, lambda_ij, lambda_ij) - ENERGY_FUNC.xuModel.getStressTensorNegative(lambda_ij, lambda_ij, lambda_ij)) * dir_ij;
216 }
217
218 totalEnergy += deltaEnergy;
219 totalEnergyGradient += deltaEnergyGradient;
220 }
221 }
222
223 energy[pId] = totalEnergy * V_i;
224 energyGradient[pId] = totalEnergyGradient * V_i;
225 }
226
227 template <typename Coord>
228 __global__ void HM_Chebyshev_Acceleration(DArray<Coord> next_X, DArray<Coord> X, DArray<Coord> prev_X, float omega)
229 {
230 int pId = blockDim.x * blockIdx.x + threadIdx.x;
231 if (pId >= prev_X.size()) return;
232
233 next_X[pId] = (next_X[pId] - X[pId]) * 0.666 + X[pId];
234
235 next_X[pId] = omega * (next_X[pId] - prev_X[pId]) + prev_X[pId];
236 }
237
238 template <typename Real, typename Coord, typename Matrix, typename Bond>
239 __global__ void HM_ComputeStepLength(
240 DArray<Real> stepLength,
241 DArray<Coord> gradient,
242 DArray<Coord> energyGradient,
243 DArray<Real> volume,
244 DArray<Matrix> A,
245 DArray<Real> energy,
246 DArrayList<Bond> restShapes)
247 {
248 int pId = blockDim.x * blockIdx.x + threadIdx.x;
249 if (pId >= stepLength.size()) return;
250
251 //TODO: replace 1000 with an input
252 Real mass_i = volume[pId] * 1000.0;
253 Real energy_i = energy[pId];
254
255 Real deltaE_i = abs(energyGradient[pId].dot(gradient[pId]));
256
257 Real alpha = deltaE_i < EPSILON || deltaE_i < energy_i ? Real(1) : energy_i / deltaE_i;
258
259 alpha /= Real(1 + restShapes[pId].size());
260
261 stepLength[pId] = alpha;
262 }
263
264 template <typename Coord>
265 __global__ void HM_FixCoM(
266 DArray<Coord> positions,
267 DArray<Attribute> atts,
268 Coord CoM)
269 {
270 int pId = blockDim.x * blockIdx.x + threadIdx.x;
271 if (pId >= positions.size()) return;
272 if ((positions[pId] - CoM).norm() < 0.05)
273 atts[pId].setFixed();
274 }
275
276 template <typename Coord, typename Bond>
277 __global__ void K_UpdateRestShape(
278 DArrayList<Bond> shape,
279 DArrayList<int> nbr,
280 DArray<Coord> pos)
281 {
282 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
283 if (pId >= pos.size()) return;
284
285 Bond np;
286
287 List<Bond>& rest_shape_i = shape[pId];
288 List<int>& list_id_i = nbr[pId];
289 int nbSize = list_id_i.size();
290 for (int ne = 0; ne < nbSize; ne++)
291 {
292 int j = list_id_i[ne];
293 np.index = j;
294 np.pos = pos[j];
295 np.weight = 1;
296
297 rest_shape_i.insert(np);
298 if (pId == j)
299 {
300 Bond np_0 = rest_shape_i[0];
301 rest_shape_i[0] = np;
302 rest_shape_i[ne] = np_0;
303 }
304 }
305 }
306
307 template <typename Coord>
308 __global__ void test_HM_UpdatePosition(
309 DArray<Coord> position,
310 DArray<Coord> velocity,
311 DArray<Coord> y_next,
312 DArray<Coord> position_old,
313 DArray<Attribute> att,
314 Real dt)
315 {
316 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
317 if (pId >= position.size()) return;
318
319 if (att[pId].isDynamic()) {
320 position[pId] = y_next[pId];
321 velocity[pId] += (position[pId] - position_old[pId]) / dt;
322 }
323 }
324
325 template<typename Real>
326 __device__ Real LimitStrain(Real s, Real slimit)
327 {
328 Real ret;
329 // if (s < 0)
330 // {
331 // ret = 0;// s + 0.05 * (slimit - s);
332 // }
333 // // else if (s > 1.1 / slimit)
334 // // {
335 // // ret = s - 0.05 * (s - 1 / slimit);
336 // // }
337 // else
338 {
339 ret = clamp(s, Real(0.2), 1 / slimit);
340 }
341
342 return ret;
343 }
344
345 template <typename Real, typename Coord, typename Matrix, typename Bond>
346 __global__ void HM_ComputeF(
347 DArray<Matrix> F,
348 DArray<Coord> eigens,
349 DArray<Matrix> invK,
350 DArray<bool> validOfK,
351 DArray<Matrix> matU,
352 DArray<Matrix> matV,
353 DArray<Matrix> Rots,
354 DArray<Coord> X,
355 DArray<Coord> Y,
356 DArrayList<Bond> bonds,
357 Real strainLimiting,
358 Real horizon)
359 {
360 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
361 if (pId >= Y.size()) return;
362 /*if (pId % 100 == 0)
363 printf("%d %d \n", pId, pId);*/
364
365 Coord x_i = X[pId];
366
367 Real total_weight = Real(0);
368 Matrix matL_i(0);
369 Matrix matK_i(0);
370
371#ifdef DEBUG_INFO
372 if (pId == 497)
373 {
374 printf("Position in HM_ComputeF %d: %f %f %f \n", pId, Y[pId][0], Y[pId][1], Y[pId][2]);
375 }
376#endif // DEBUG_INFO
377
378 List<Bond>& bonds_i = bonds[pId];
379
380 Real maxDist = Real(0);
381 for (int ne = 0; ne < bonds_i.size(); ne++)
382 {
383 maxDist = max(maxDist, bonds_i[ne].xi.norm());
384 }
385 maxDist = maxDist < EPSILON ? Real(1) : maxDist;
386
387 //printf("maxdist %d: %f; \n\n", pId, maxDist);
388
389#ifdef DEBUG_INFO
390 printf("Max distance %d: %f \n", pId, maxDist);
391#endif // DEBUG_INFO
392
393 for (int ne = 0; ne < bonds_i.size(); ne++)
394 {
395 Bond bond_ij = bonds[pId][ne];
396 int j = bond_ij.idx;
397 Coord x_j = X[j];
398 Real r = (x_i - x_j).norm();
399
400 if (r > EPSILON)
401 {
402 Coord p = (Y[j] - Y[pId]) / maxDist;
403 Coord q = (x_j - x_i) / maxDist;
404 //Real weight = D_Weight(r, horizon);
405 Real weight = 1.;
406
407 matL_i(0, 0) += p[0] * q[0] * weight; matL_i(0, 1) += p[0] * q[1] * weight; matL_i(0, 2) += p[0] * q[2] * weight;
408 matL_i(1, 0) += p[1] * q[0] * weight; matL_i(1, 1) += p[1] * q[1] * weight; matL_i(1, 2) += p[1] * q[2] * weight;
409 matL_i(2, 0) += p[2] * q[0] * weight; matL_i(2, 1) += p[2] * q[1] * weight; matL_i(2, 2) += p[2] * q[2] * weight;
410
411 matK_i(0, 0) += q[0] * q[0] * weight; matK_i(0, 1) += q[0] * q[1] * weight; matK_i(0, 2) += q[0] * q[2] * weight;
412 matK_i(1, 0) += q[1] * q[0] * weight; matK_i(1, 1) += q[1] * q[1] * weight; matK_i(1, 2) += q[1] * q[2] * weight;
413 matK_i(2, 0) += q[2] * q[0] * weight; matK_i(2, 1) += q[2] * q[1] * weight; matK_i(2, 2) += q[2] * q[2] * weight;
414
415 total_weight += weight;
416
417#ifdef DEBUG_INFO
418 if (pId == 497)
419 {
420 printf("%d Neighbor %d: %f %f %f \n", pId, j, Y[j][0], Y[j][1], Y[j][2]);
421 }
422#endif // DEBUG_INFO
423 }
424 }
425
426 if (total_weight > EPSILON)
427 {
428 matL_i *= (1.0f / total_weight);
429 matK_i *= (1.0f / total_weight);
430 }
431
432 Matrix R, U, D, V;
433 polarDecomposition(matK_i, R, U, D, V);
434
435#ifdef DEBUG_INFO
436 if (pId == 497)
437 {
438 Matrix mat_out = matK_i;
439 printf("matK_i: \n %f %f %f \n %f %f %f \n %f %f %f \n\n",
440 mat_out(0, 0), mat_out(0, 1), mat_out(0, 2),
441 mat_out(1, 0), mat_out(1, 1), mat_out(1, 2),
442 mat_out(2, 0), mat_out(2, 1), mat_out(2, 2));
443
444 mat_out = matL_i;
445 printf("matL_i: \n %f %f %f \n %f %f %f \n %f %f %f \n\n",
446 mat_out(0, 0), mat_out(0, 1), mat_out(0, 2),
447 mat_out(1, 0), mat_out(1, 1), mat_out(1, 2),
448 mat_out(2, 0), mat_out(2, 1), mat_out(2, 2));
449
450 mat_out = U * D * V.transpose();
451 printf("matK polar: \n %f %f %f \n %f %f %f \n %f %f %f \n\n",
452 mat_out(0, 0), mat_out(0, 1), mat_out(0, 2),
453 mat_out(1, 0), mat_out(1, 1), mat_out(1, 2),
454 mat_out(2, 0), mat_out(2, 1), mat_out(2, 2));
455
456 printf("Horizon: %f; Det: %f \n", horizon, matK_i.determinant());
457 }
458#endif // DEBUG_INFO
459
460 Real maxK = maximum(abs(D(0, 0)), maximum(abs(D(1, 1)), abs(D(2, 2))));
461 Real minK = minimum(abs(D(0, 0)), minimum(abs(D(1, 1)), abs(D(2, 2))));
462
463 bool valid_K = (minK < EPSILON) ? false : true;
464
465 validOfK[pId] = valid_K;
466
467 Matrix F_i;
468 if (valid_K)
469 {
470 invK[pId] = matK_i.inverse();
471 F_i = matL_i * matK_i.inverse();
472 }
473 else
474 {
475 Real threshold = 0.0001f * horizon;
476 D(0, 0) = D(0, 0) > threshold ? 1.0 / D(0, 0) : 1.0;
477 D(1, 1) = D(1, 1) > threshold ? 1.0 / D(1, 1) : 1.0;
478 D(2, 2) = D(2, 2) > threshold ? 1.0 / D(2, 2) : 1.0;
479 invK[pId] = V * D * U.transpose();
480 F_i = matL_i * V * D * U.transpose();
481 }
482
483 polarDecomposition(F_i, R, U, D, V);
484
485 // Real minDist = Real(1000.0f);
486 // for (int ne = 0; ne < size_i; ne++)
487 // {
488 // Bond np_j = restShapes[pId][ne];
489 // int j = np_j.index;
490 // Coord rest_pos_j = np_j.pos;
491 // Real r = (rest_pos_i - rest_pos_j).norm();
492 //
493 // if(r > EPSILON)
494 // minDist = min(minDist, r);
495 // }
496 //
497 // printf("ID %d: min value %f \n", pId, minDist);
498
499 if (F_i.determinant() < maxDist * EPSILON)
500 {
501 }
502 else
503 {
504 matU[pId] = U;
505 matV[pId] = V;
506 }
507
508 // Coord n0 = Coord(1, 0, 0);
509 // Coord n1 = Coord(0, 1, 0);
510 // Coord n2 = Coord(0, 0, 1);
511 //
512 // Coord n0_rot = R * n0;
513 // Coord n1_rot = R * n1;
514 // Coord n2_rot = R * n2;
515
516 Real l0 = D(0, 0);
517 Real l1 = D(1, 1);
518 Real l2 = D(2, 2);
519
520 const Real slimit = strainLimiting;
521
522 l0 = clamp(l0, slimit, 1 / slimit);
523 l1 = clamp(l1, slimit, 1 / slimit);
524 l2 = clamp(l2, slimit, 1 / slimit);
525
526 D(0, 0) = l0;
527 D(1, 1) = l1;
528 D(2, 2) = l2;
529
530 eigens[pId] = Coord(l0, l1, l2);
531
532 F[pId] = U * D * V.transpose();
533
534 // if (F_i.determinant() < maxDist*EPSILON)
535 // {
536 // R = Rots[pId];
537 // }
538 // else
539 // {
540 // Rots[pId] = R;
541 // }
542 //
543 // Coord n0 = Coord(1, 0, 0);
544 // Coord n1 = Coord(0, 1, 0);
545 // Coord n2 = Coord(0, 0, 1);
546 //
547 // Coord n0_rot = R * n0;
548 // Coord n1_rot = R * n1;
549 // Coord n2_rot = R * n2;
550 //
551 // Real l0 = n0_rot.dot(F_i*n0);
552 // Real l1 = n1_rot.dot(F_i*n1);
553 // Real l2 = n2_rot.dot(F_i*n2);
554 //
555 // matU[pId] = R;
556 // matV[pId] = Matrix::identityMatrix();
557 //
558 // const Real slimit = strainLimiting;
559 //
560 // l0 = clamp(l0, slimit, 1 / slimit);
561 // l1 = clamp(l1, slimit, 1 / slimit);
562 // l2 = clamp(l2, slimit, 1 / slimit);
563 //
564 // D(0, 0) = l0;
565 // D(1, 1) = l1;
566 // D(2, 2) = l2;
567 //
568 // eigens[pId] = Coord(l0, l1, l2);
569 //
570 // F[pId] = R * D;
571 }
572
573 template <typename Real, typename Coord, typename Matrix, typename Bond>
574 __global__ void HM_JacobiStepNonsymmetric(
575 DArray<Coord> source,
576 DArray<Matrix> A,
577 DArray<Coord> y_pre,
578 DArray<Matrix> matU,
579 DArray<Matrix> matV,
580 DArray<Matrix> matR,
581 DArray<Coord> eigen,
582 DArray<bool> validOfK,
583 DArray<Matrix> F,
584 DArray<Coord> X,
585 DArrayList<Bond> bonds,
586 Real horizon,
587 DArray<Real> volume,
588 DArrayList<Real> volumePair,
589 Real dt,
590 EnergyType type)
591 {
592 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
593 if (pId >= y_pre.size()) return;
594
595 Coord x_i = X[pId];
596 int size_i = bonds[pId].size();
597
598 Real kappa = 4 / (3 * M_PI * horizon * horizon * horizon);
599 Real lambda_i1 = eigen[pId][0];
600 Real lambda_i2 = eigen[pId][1];
601 Real lambda_i3 = eigen[pId][2];
602
603 Matrix U_i = matU[pId];
604 Matrix V_i = matV[pId];
605
606 Matrix S1_i;
607 Matrix S2_i;
608 if (type == StVK) {
609 S1_i = ENERGY_FUNC.stvkModel.getStressTensorPositive(lambda_i1, lambda_i2, lambda_i3);
610 S2_i = ENERGY_FUNC.stvkModel.getStressTensorNegative(lambda_i1, lambda_i2, lambda_i3);
611 }
612 else if (type == NeoHooekean) {
613 S1_i = ENERGY_FUNC.neohookeanModel.getStressTensorPositive(lambda_i1, lambda_i2, lambda_i3);
614 S2_i = ENERGY_FUNC.neohookeanModel.getStressTensorNegative(lambda_i1, lambda_i2, lambda_i3);
615 }
616 else if (type == Linear) {
617 S1_i = ENERGY_FUNC.linearModel.getStressTensorPositive(lambda_i1, lambda_i2, lambda_i3);
618 S2_i = ENERGY_FUNC.linearModel.getStressTensorNegative(lambda_i1, lambda_i2, lambda_i3);
619 }
620 else if (type == Xuetal) {
621 S1_i = ENERGY_FUNC.xuModel.getStressTensorPositive(lambda_i1, lambda_i2, lambda_i3);
622 S2_i = ENERGY_FUNC.xuModel.getStressTensorNegative(lambda_i1, lambda_i2, lambda_i3);
623 }
624 else if (type == MooneyRivlin) {
625 S1_i = ENERGY_FUNC.mrModel.getStressTensorPositive(lambda_i1, lambda_i2, lambda_i3);
626 S2_i = ENERGY_FUNC.mrModel.getStressTensorNegative(lambda_i1, lambda_i2, lambda_i3);
627 }
628 else if (type == Fung) {
629 S1_i = ENERGY_FUNC.fungModel.getStressTensorPositive(lambda_i1, lambda_i2, lambda_i3);
630 S2_i = ENERGY_FUNC.fungModel.getStressTensorNegative(lambda_i1, lambda_i2, lambda_i3);
631 }
632 else if (type == Ogden) {
633 S1_i = ENERGY_FUNC.ogdenModel.getStressTensorPositive(lambda_i1, lambda_i2, lambda_i3);
634 S2_i = ENERGY_FUNC.ogdenModel.getStressTensorNegative(lambda_i1, lambda_i2, lambda_i3);
635 }
636 else if (type == Yeoh) {
637 S1_i = ENERGY_FUNC.yeohModel.getStressTensorPositive(lambda_i1, lambda_i2, lambda_i3);
638 S2_i = ENERGY_FUNC.yeohModel.getStressTensorNegative(lambda_i1, lambda_i2, lambda_i3);
639 }
640 else if (type == ArrudaBoyce) {
641 S1_i = ENERGY_FUNC.abModel.getStressTensorPositive(lambda_i1, lambda_i2, lambda_i3);
642 S2_i = ENERGY_FUNC.abModel.getStressTensorNegative(lambda_i1, lambda_i2, lambda_i3);
643 }
644
645 Matrix PK1_i = U_i * S1_i * U_i.transpose();
646 Matrix PK2_i = U_i * S2_i * V_i.transpose();
647
648 //Real Vol_i = volume[pId];
649
650 Matrix F_i = F[pId];
651
652 Coord y_pre_i = y_pre[pId];
653
654 bool K_valid = validOfK[pId];
655
656 Matrix mat_i(0);
657 Coord source_i(0);
658
659 for (int ne = 0; ne < size_i; ne++)
660 {
661 Bond bond_ij = bonds[pId][ne];
662 int j = bond_ij.idx;
663 Coord x_j = X[j];
664 Coord y_pre_j = y_pre[j];
665 Real r = bond_ij.xi.norm();
666
667 if (r > EPSILON)
668 {
669 Real weight_ij = 1.;
670 //Real Vol_j = volume[j];
671 Real V_ij = volumePair[pId][ne];
672
673 Coord y_ij = y_pre_i - y_pre_j;
674
675 Real lambda = y_ij.norm() / r;
676
677 //const Real scale = Vol_i * Vol_j*kappa;
678 const Real scale = V_ij * V_ij * kappa;
679 const Real sw_ij = dt * dt * scale * weight_ij / (r * r);
680
681 bool K_valid_ij = K_valid & validOfK[j];
682
683 Matrix S1_ij;
684 Matrix S2_ij;
685 if (type == StVK) {
686 S1_ij = ENERGY_FUNC.stvkModel.getStressTensorPositive(lambda, lambda, lambda);
687 S2_ij = ENERGY_FUNC.stvkModel.getStressTensorNegative(lambda, lambda, lambda);
688 }
689 else if (type == NeoHooekean) {
690 S1_ij = ENERGY_FUNC.neohookeanModel.getStressTensorPositive(lambda, lambda, lambda);
691 S2_ij = ENERGY_FUNC.neohookeanModel.getStressTensorNegative(lambda, lambda, lambda);
692 }
693 else if (type == Linear) {
694 S1_ij = ENERGY_FUNC.linearModel.getStressTensorPositive(lambda, lambda, lambda);
695 S2_ij = ENERGY_FUNC.linearModel.getStressTensorNegative(lambda, lambda, lambda);
696 }
697 else if (type == Xuetal) {
698 S1_ij = ENERGY_FUNC.xuModel.getStressTensorPositive(lambda, lambda, lambda);
699 S2_ij = ENERGY_FUNC.xuModel.getStressTensorNegative(lambda, lambda, lambda);
700 }
701 else if (type == MooneyRivlin) {
702 S1_ij = ENERGY_FUNC.mrModel.getStressTensorPositive(lambda, lambda, lambda);
703 S2_ij = ENERGY_FUNC.mrModel.getStressTensorNegative(lambda, lambda, lambda);
704 }
705 else if (type == Fung) {
706 S1_ij = ENERGY_FUNC.fungModel.getStressTensorPositive(lambda, lambda, lambda);
707 S2_ij = ENERGY_FUNC.fungModel.getStressTensorNegative(lambda, lambda, lambda);
708 }
709 else if (type == Ogden) {
710 S1_ij = ENERGY_FUNC.ogdenModel.getStressTensorPositive(lambda, lambda, lambda);
711 S2_ij = ENERGY_FUNC.ogdenModel.getStressTensorNegative(lambda, lambda, lambda);
712 }
713 else if (type == Yeoh) {
714 S1_ij = ENERGY_FUNC.yeohModel.getStressTensorPositive(lambda, lambda, lambda);
715 S2_ij = ENERGY_FUNC.yeohModel.getStressTensorNegative(lambda, lambda, lambda);
716 }
717 else if (type == ArrudaBoyce) {
718 S1_ij = ENERGY_FUNC.abModel.getStressTensorPositive(lambda, lambda, lambda);
719 S2_ij = ENERGY_FUNC.abModel.getStressTensorNegative(lambda, lambda, lambda);
720 }
721
722 Matrix PK1_ij = K_valid_ij ? PK1_i : S1_ij;
723 Matrix PK2_ij = K_valid_ij ? PK2_i : S2_ij;
724
725 Coord dir_ij = lambda > EPSILON ? y_ij.normalize() : Coord(1, 0, 0);
726
727 Coord x_ij = K_valid_ij ? x_i - x_j : dir_ij * (x_i - x_j).norm();
728
729 Coord ds_ij = sw_ij * PK1_ij * y_pre_j + sw_ij * PK2_ij * x_ij;
730 Coord ds_ji = sw_ij * PK1_ij * y_pre_i - sw_ij * PK2_ij * x_ij;
731
732 Matrix mat_ij = sw_ij * PK1_ij;
733
734 source_i += ds_ij;
735 mat_i += mat_ij;
736
737 atomicAdd(&source[j][0], ds_ji[0]);
738 atomicAdd(&source[j][1], ds_ji[1]);
739 atomicAdd(&source[j][2], ds_ji[2]);
740
741 atomicAdd(&A[j](0, 0), mat_ij(0, 0));
742 atomicAdd(&A[j](0, 1), mat_ij(0, 1));
743 atomicAdd(&A[j](0, 2), mat_ij(0, 2));
744 atomicAdd(&A[j](1, 0), mat_ij(1, 0));
745 atomicAdd(&A[j](1, 1), mat_ij(1, 1));
746 atomicAdd(&A[j](1, 2), mat_ij(1, 2));
747 atomicAdd(&A[j](2, 0), mat_ij(2, 0));
748 atomicAdd(&A[j](2, 1), mat_ij(2, 1));
749 atomicAdd(&A[j](2, 2), mat_ij(2, 2));
750 }
751 }
752
753 atomicAdd(&source[pId][0], source_i[0]);
754 atomicAdd(&source[pId][1], source_i[1]);
755 atomicAdd(&source[pId][2], source_i[2]);
756
757 atomicAdd(&A[pId](0, 0), mat_i(0, 0));
758 atomicAdd(&A[pId](0, 1), mat_i(0, 1));
759 atomicAdd(&A[pId](0, 2), mat_i(0, 2));
760 atomicAdd(&A[pId](1, 0), mat_i(1, 0));
761 atomicAdd(&A[pId](1, 1), mat_i(1, 1));
762 atomicAdd(&A[pId](1, 2), mat_i(1, 2));
763 atomicAdd(&A[pId](2, 0), mat_i(2, 0));
764 atomicAdd(&A[pId](2, 1), mat_i(2, 1));
765 atomicAdd(&A[pId](2, 2), mat_i(2, 2));
766 }
767
768 template <typename Real, typename Coord, typename Matrix>
769 __global__ void HM_ComputeNextPosition(
770 DArray<Coord> y_next,
771 DArray<Coord> y_pre,
772 DArray<Coord> y_old,
773 DArray<Coord> source,
774 DArray<Matrix> A,
775 DArray<Real> volume,
776 DArray<Attribute> attribute,
777 EnergyType type)
778 {
779 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
780 if (pId >= y_next.size()) return;
781
782 //TODO: replace 1000 with an input
783 Real mass_i = volume[pId] * 1000;
784
785 Matrix mat_i = A[pId] + mass_i * Matrix::identityMatrix();
786 Coord src_i = source[pId] + mass_i * y_old[pId];
787
788 y_next[pId] = attribute[pId].isDynamic() ? mat_i.inverse() * src_i : y_pre[pId];
789 }
790
791 template <typename Matrix>
792 __global__ void HM_InitRotation(
793 DArray<Matrix> U,
794 DArray<Matrix> V,
795 DArray<Matrix> R)
796 {
797 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
798 if (pId >= U.size()) return;
799
800 U[pId] = Matrix::identityMatrix();
801 V[pId] = Matrix::identityMatrix();
802 R[pId] = Matrix::identityMatrix();
803 }
804
805 template<typename TDataType>
806 void SemiImplicitHyperelasticitySolver<TDataType>::resizeAllFields()
807 {
808 uint num = this->inY()->size();
809
810 if (m_F.size() == num)
811 return;
812
813 m_F.resize(num);
814 m_invF.resize(num);
815 m_invK.resize(num);
816 m_validOfK.resize(num);
817 m_validOfF.resize(num);
818
819 m_eigenValues.resize(num);
820 m_matU.resize(num);
821 m_matV.resize(num);
822 m_matR.resize(num);
823
824 m_energy.resize(num);
825 m_alpha.resize(num);
826 m_gradient.resize(num);
827 mEnergyGradient.resize(num);
828
829 y_next.resize(num);
830 y_pre.resize(num);
831 y_current.resize(num);
832 y_residual.resize(num);
833 y_gradC.resize(num);
834
835 m_source.resize(num);
836 m_A.resize(num);
837
838 this->mPosBuf.resize(num);
839
840 m_fraction.resize(num);
841
842 cuExecute(m_matU.size(),
843 HM_InitRotation,
844 m_matU,
845 m_matV,
846 m_matR);
847
848 m_reduce = Reduction<Real>::Create(num);
849 }
850
851 template<typename TDataType>
852 void SemiImplicitHyperelasticitySolver<TDataType>::enforceHyperelasticity()
853 {
854 //std::cout << "enforceElasticity\n";
855
856 resizeAllFields();
857
858 int numOfParticles = this->inY()->size();
859 uint pDims = cudaGridSize(numOfParticles, BLOCK_SIZE);
860
861 std::cout << "enforceElasticity " << numOfParticles << std::endl;
862
863 this->mWeights.reset();
864
865 Log::sendMessage(Log::User, "\n \n \n \n *************solver start!!!***************");
866
867 /**************************** Jacobi method ************************************************/
868 // initialize y_now, y_next_iter
869 y_current.assign(this->inY()->getData());
870 y_next.assign(this->inY()->getData());
871 this->mPosBuf.assign(this->inY()->getData());
872
873 // do Jacobi method Loop
874 bool convergeFlag = false; // converge or not
875 int iterCount = 0;
876
877 Real omega;
878 Real alpha = 1.0f;
879 std::cout << "enforceElasticity 2 " << numOfParticles << std::endl;
880
881 while (iterCount < this->varIterationNumber()->getData()) {
882 //printf("%.3lf\n", inHorizon()->getData());
883
884
885 m_source.reset();
886 m_A.reset();
887 HM_ComputeF << <pDims, BLOCK_SIZE >> > (
888 m_F,
889 m_eigenValues,
890 m_invK,
891 m_validOfK,
892 m_matU,
893 m_matV,
894 m_matR,
895 this->inX()->getData(),
896 y_current,
897 this->inBonds()->getData(),
898 this->varStrainLimiting()->getData(),
899 this->inHorizon()->getData());
900 cuSynchronize();
901
902 HM_JacobiStepNonsymmetric << <pDims, BLOCK_SIZE >> > (
903 m_source,
904 m_A,
905 y_current,
906 m_matU,
907 m_matV,
908 m_matR,
909 m_eigenValues,
910 m_validOfK,
911 m_F,
912 this->inX()->getData(),
913 this->inBonds()->getData(),
914 this->inHorizon()->getData(),
915 this->inVolume()->getData(),
916 this->inVolumePair()->getData(),
917 this->inTimeStep()->getData(),
918 this->inEnergyType()->getData());
919 cuSynchronize();
920
921 cuExecute(y_next.size(),
922 HM_ComputeNextPosition,
923 y_next,
924 y_current,
925 this->mPosBuf,
926 m_source,
927 m_A,
928 this->inVolume()->getData(),
929 this->inAttribute()->getData(),
930 this->inEnergyType()->getData());
931
932 //sub steping
933 cuExecute(m_gradient.size(),
934 HM_ComputeGradient,
935 m_gradient,
936 y_current,
937 y_next);
938
939
940 cuExecute(m_energy.size(),
941 HM_Compute1DEnergy,
942 m_energy,
943 mEnergyGradient,
944 this->inX()->getData(),
945 y_current,
946 m_F,
947 this->inVolume()->getData(),
948 m_validOfK,
949 m_eigenValues,
950 this->inBonds()->getData(),
951 this->inEnergyType()->getData());
952
953 m_alphaCompute = this->varIsAlphaComputed()->getData();
954 if (this->m_alphaCompute) {
955 cuExecute(m_alpha.size(),
956 HM_ComputeStepLength,
957 m_alpha,
958 m_gradient,
959 mEnergyGradient,
960 this->inVolume()->getData(),
961 m_A,
962 m_energy,
963 this->inBonds()->getData());
964
965 cuExecute(m_gradient.size(),
966 HM_ComputeCurrentPosition,
967 m_gradient,
968 y_current,
969 y_next,
970 this->inAttribute()->getData(),
971 m_alpha);
972
973 }
974 else {
975 alpha = Real(0.15f);
976 cuExecute(m_gradient.size(),
977 HM_ComputeCurrentPosition,
978 m_gradient,
979 y_current,
980 y_next,
981 alpha);
982 }
983
984 y_current.assign(y_next);
985
986 iterCount++;
987 }
988
989 cuExecute(this->inY()->getDataPtr()->size(),
990 test_HM_UpdatePosition,
991 this->inY()->getData(),
992 this->inVelocity()->getData(),
993 y_next,
994 this->mPosBuf,
995 this->inAttribute()->getData(),
996 this->inTimeStep()->getData());
997 printf("outside\n");
998 }
999
1000 DEFINE_CLASS(SemiImplicitHyperelasticitySolver);
1001}
Generated by doxygen 1.13.2