PeriDyno 1.0.0
Loading...
Searching...
No Matches
CoSemiImplicitHyperelasticitySolver.cu
Go to the documentation of this file.
1#include "CoSemiImplicitHyperelasticitySolver.h"
2#include "Matrix/MatrixFunc.h"
3#include "ParticleSystem/Module/Kernel.h"
4#include "curand_kernel.h"
5#include "Algorithm/CudaRand.h"
6
7namespace dyno
8{
9 IMPLEMENT_TCLASS(CoSemiImplicitHyperelasticitySolver, TDataType);
10
11 __constant__ EnergyModels<Real> ENERGY_FUNC;
12
13 template<typename Real>
14 __device__ Real constantWeight(Real r, Real h)
15 {
16 Real d = h / r;
17 return d * d;
18 }
19
20 template<typename TDataType>
21 CoSemiImplicitHyperelasticitySolver<TDataType>::CoSemiImplicitHyperelasticitySolver()
22 : LinearElasticitySolver<TDataType>()
23 {
24 this->varIterationNumber()->setValue(30);
25
26 mContactRule = std::make_shared<ContactRule<TDataType>>();
27
28 this->inTriangularMesh()->connect(mContactRule->inTriangularMesh());
29 this->inOldPosition()->connect(mContactRule->inOldPosition());
30 this->inVelocity()->connect(mContactRule->inVelocity());
31 this->inTimeStep()->connect(mContactRule->inTimeStep());
32 this->inUnit()->connect(mContactRule->inUnit());
33 }
34
35 template<typename TDataType>
36 void CoSemiImplicitHyperelasticitySolver<TDataType>::initializeVolume()
37 {
38 int numOfParticles = this->inY()->getData().size();
39 uint pDims = cudaGridSize(numOfParticles, BLOCK_SIZE);
40 std::cout << "dev: " << numOfParticles << " particles\n";
41 HM_InitVolume << <pDims, BLOCK_SIZE >> > (m_volume, m_objectVolume, m_objectVolumeSet, m_particleVolume, m_particleVolumeSet);
42 }
43
44 template<typename TDataType>
45 CoSemiImplicitHyperelasticitySolver<TDataType>::~CoSemiImplicitHyperelasticitySolver()
46 {
47 this->mWeights.clear();
48 this->mDisplacement.clear();
49 this->mInvK.clear();
50 this->mF.clear();
51 this->mPosBuf.clear();
52 this->mPosBuf_March.clear();
53 }
54
55 template <typename Real, typename Coord, typename Matrix>
56 __global__ void HM_ComputeEnergy(
57 DArray<Real> energy,
58 DArray<Coord> eigens,
59 EnergyType type)
60 {
61 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
62 if (pId >= energy.size()) return;
63
64 Coord eigen_i = eigens[pId];
65
66 if (type == StVK) {
67 energy[pId] = ENERGY_FUNC.stvkModel.getEnergy(eigen_i[0], eigen_i[1], eigen_i[2]);
68 }
69 else if (type == NeoHooekean) {
70 energy[pId] = ENERGY_FUNC.neohookeanModel.getEnergy(eigen_i[0], eigen_i[1], eigen_i[2]);
71 }
72 else if (type == Linear) {
73 energy[pId] = ENERGY_FUNC.linearModel.getEnergy(eigen_i[0], eigen_i[1], eigen_i[2]);
74 }
75 else if (type == Xuetal) {
76 energy[pId] = ENERGY_FUNC.xuModel.getEnergy(eigen_i[0], eigen_i[1], eigen_i[2]);
77 }
78 }
79
80
81 template<typename TDataType>
82 void CoSemiImplicitHyperelasticitySolver<TDataType>::solveElasticity()
83 {
84 cudaMemcpyToSymbol(ENERGY_FUNC, this->inEnergyModels()->constDataPtr().get(), sizeof(EnergyModels<Real>));
85
86 enforceHyperelasticity();
87 }
88
89 template <typename Coord>
90 __global__ void HM_ComputeGradient(
91 DArray<Coord> grad,
92 DArray<Coord> y_pre,
93 DArray<Coord> y_next)
94 {
95 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
96 if (pId >= y_next.size()) return;
97
98 grad[pId] = y_next[pId] - y_pre[pId];
99 }
100
101 template <typename Real, typename Coord>
102 __global__ void HM_ComputeCurrentPosition(
103 DArray<Coord> grad,
104 DArray<Coord> y_current,
105 DArray<Coord> y_next,
106 DArray<Real> alpha)
107 {
108 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
109 if (pId >= y_next.size()) return;
110
111 y_next[pId] = y_current[pId] + alpha[pId] * grad[pId];
112 }
113
114 template <typename Real, typename Coord>
115 __global__ void HM_ComputeCurrentPosition(
116 DArray<Coord> grad,
117 DArray<Coord> y_current,
118 DArray<Coord> y_next,
119 Real alpha)
120 {
121 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
122 if (pId >= y_next.size()) return;
123
124 y_next[pId] = y_current[pId] + alpha * grad[pId];
125 }
126
127 template <typename Real, typename Coord, typename Matrix, typename Bond>
128 __global__ void HM_Compute1DEnergy(
129 DArray<Real> energy,
130 DArray<Coord> energyGradient,
131 DArray<Coord> X,
132 DArray<Coord> pos_current,
133 DArray<Matrix> F,
134 DArray<Real> volume,
135 DArray<bool> validOfK,
136 DArray<Coord> eigenValues,
137 DArrayList<Bond> bonds,
138 EnergyType type)
139 {
140 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
141 if (pId >= energy.size()) return;
142
143 Coord pos_current_i = pos_current[pId];
144
145 Real totalEnergy = 0.0f;
146 Coord totalEnergyGradient = Coord(0);
147 Real V_i = volume[pId];
148
149 int size_i = bonds[pId].size();
150
151 Coord x_i = X[pId];
152 Coord eigen_value_i = eigenValues[pId];
153 bool valid_i = validOfK[pId];
154
155 Matrix F_i = F[pId];
156
157 for (int ne = 0; ne < size_i; ne++)
158 {
159 Bond bond_ij = bonds[pId][ne];
160 int j = bond_ij.idx;
161 Coord pos_current_j = pos_current[j];
162 Coord x_j = X[j];
163 Real r = (x_j - x_i).norm();
164
165 Real V_j = volume[j];
166
167 if (r > EPSILON)
168 {
169 Real norm_ij = (pos_current_j - pos_current_i).norm();
170 Real lambda_ij = norm_ij / r;
171
172 Real deltaEnergy;
173 Coord deltaEnergyGradient;
174 Coord dir_ij = norm_ij < EPSILON ? Coord(0) : (pos_current_i - pos_current_j) / (r);
175
176 if (type == StVK) {
177 deltaEnergy = V_j * ENERGY_FUNC.stvkModel.getEnergy(lambda_ij, lambda_ij, lambda_ij);
178 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;
179 }
180 else if (type == NeoHooekean) {
181 deltaEnergy = V_j * ENERGY_FUNC.neohookeanModel.getEnergy(lambda_ij, lambda_ij, lambda_ij);
182 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;
183 }
184 else if (type == Linear) {
185 deltaEnergy = V_j * ENERGY_FUNC.linearModel.getEnergy(lambda_ij, lambda_ij, lambda_ij);
186 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;
187 }
188 else if (type == Xuetal) {
189 deltaEnergy = V_j * ENERGY_FUNC.xuModel.getEnergy(lambda_ij, lambda_ij, lambda_ij);
190 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;
191 }
192
193 totalEnergy += deltaEnergy;
194 totalEnergyGradient += deltaEnergyGradient;
195 }
196 }
197
198 energy[pId] = totalEnergy * V_i;
199 energyGradient[pId] = totalEnergyGradient * V_i;
200 }
201
202 template <typename Coord>
203 __global__ void HM_Chebyshev_Acceleration(DArray<Coord> next_X, DArray<Coord> X, DArray<Coord> prev_X, float omega)
204 {
205 int pId = blockDim.x * blockIdx.x + threadIdx.x;
206 if (pId >= prev_X.size()) return;
207
208 next_X[pId] = (next_X[pId] - X[pId]) * 0.666 + X[pId];
209
210 next_X[pId] = omega * (next_X[pId] - prev_X[pId]) + prev_X[pId];
211 }
212
213 template <typename Real, typename Coord, typename Matrix, typename Bond>
214 __global__ void HM_ComputeStepLength(
215 DArray<Real> stepLength,
216 DArray<Coord> gradient,
217 DArray<Coord> energyGradient,
218 DArray<Real> volume,
219 DArray<Matrix> A,
220 DArray<Real> energy,
221 DArrayList<Bond> bonds)
222 {
223 int pId = blockDim.x * blockIdx.x + threadIdx.x;
224 if (pId >= stepLength.size()) return;
225
226 Real mass_i = volume[pId] * 1000.0;
227 Real energy_i = energy[pId];
228
229 Real deltaE_i = abs(energyGradient[pId].dot(gradient[pId]));
230
231 Real alpha = deltaE_i < EPSILON || deltaE_i < energy_i ? Real(1) : energy_i / deltaE_i;
232
233 alpha /= Real(1 + bonds[pId].size());
234
235 stepLength[pId] = alpha;
236 }
237
238 template <typename Coord>
239 __global__ void HM_FixCoM(
240 DArray<Coord> positions,
241 DArray<Attribute> atts,
242 Coord CoM)
243 {
244 int pId = blockDim.x * blockIdx.x + threadIdx.x;
245 if (pId >= positions.size()) return;
246 if ((positions[pId] - CoM).norm() < 0.05)
247 atts[pId].setFixed();
248 }
249
250 template <typename Coord, typename Bond>
251 __global__ void K_UpdateRestShape(
252 DArrayList<Bond> shape,
253 DArrayList<int> nbr,
254 DArray<Coord> pos)
255 {
256 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
257 if (pId >= pos.size()) return;
258
259 Bond np;
260
261 List<Bond>& rest_shape_i = shape[pId];
262 List<int>& list_id_i = nbr[pId];
263 int nbSize = list_id_i.size();
264 for (int ne = 0; ne < nbSize; ne++)
265 {
266 int j = list_id_i[ne];
267 np.index = j;
268 np.pos = pos[j];
269 np.weight = 1;
270
271 rest_shape_i.insert(np);
272 if (pId == j)
273 {
274 Bond np_0 = rest_shape_i[0];
275 rest_shape_i[0] = np;
276 rest_shape_i[ne] = np_0;
277 }
278 }
279 }
280
281 template <typename Coord>
282 __global__ void test_HM_UpdatePosition(
283 DArray<Coord> position, //position target
284 DArray<Coord> y_next, //position reference
285 DArray<Attribute> att)
286 {
287 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
288 if (pId >= position.size()) return;
289
290 if (!att[pId].isFixed()) {
291 position[pId] = y_next[pId];
292 }
293 }
294
295 template <typename Coord>
296 __global__ void test_HM_UpdatePosition(
297 DArray<Coord> position,
298 DArray<Coord> velocity,
299 DArray<Coord> y_next,
300 DArray<Coord> position_old,
301 DArray<Attribute> att,
302 Real dt)
303 {
304 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
305 if (pId >= position.size()) return;
306
307 if (!att[pId].isFixed()) {
308 position[pId] = y_next[pId];
309
310 velocity[pId] += (position[pId] - position_old[pId]) / dt;
311 }
312 }
313
314 template <typename Real, typename Coord, typename Matrix, typename Bond>
315 __global__ void HM_ComputeF(
316 DArray<Matrix> F,
317 DArray<Coord> eigens,
318 DArray<Matrix> invK,
319 DArray<bool> validOfK,
320 DArray<Matrix> matU,
321 DArray<Matrix> matV,
322 DArray<Matrix> Rots,
323 DArray<Coord> X,
324 DArray<Coord> Y,
325 DArrayList<Bond> bonds,
326 Real horizon,
327 Real const strainLimit,
328 DArray<Coord> restNorm,
329 DArray<Coord> Norm)
330 {
331 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
332 if (pId >= Y.size()) return;
333
334 Coord x_i = X[pId];
335 int size_i = bonds[pId].size();
336 Real total_weight = Real(0);
337 Matrix matL_i(0);
338 Matrix matK_i(0);
339
340 Real t = 1;
341
342#ifdef DEBUG_INFO
343 if (pId == 497)
344 {
345 printf("Position in HM_ComputeF %d: %f %f %f \n", pId, Y[pId][0], Y[pId][1], Y[pId][2]);
346 }
347#endif // DEBUG_INFO
348 //printf("%d %d \n", pId, size_i);
349
350 Real maxDist = Real(0);
351 for (int ne = 0; ne < size_i; ne++)
352 {
353 Bond bond_ij = bonds[pId][ne];
354 int j = bond_ij.idx;
355 Coord y_j = X[j];
356 Real r = (x_i - y_j).norm();
357
358 maxDist = max(maxDist, r);
359 }
360 maxDist = maxDist < EPSILON ? Real(1) : maxDist;
361
362#ifdef DEBUG_INFO
363 printf("Max distance %d: %f \n", pId, maxDist);
364#endif // DEBUG_INFO
365
366 for (int ne = 0; ne < size_i; ne++)
367 {
368 Bond bond_ij = bonds[pId][ne];
369 int j = bond_ij.idx;
370 Coord x_j = X[j];
371 Real r = (x_i - x_j).norm();
372
373 if (r > EPSILON)
374 {
375 Real weight = Real(1);
376
377 Coord p = (Y[j] - Y[pId]) / maxDist;
378 Coord q = (x_j - x_i) / maxDist;
379
380
381 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;
382 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;
383 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;
384
385 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;
386 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;
387 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;
388
389 total_weight += weight;
390
391
392#ifdef DEBUG_INFO
393 if (pId == 497)
394 {
395 printf("%d Neighbor %d: %f %f %f \n", pId, j, Y[j][0], Y[j][1], Y[j][2]);
396 }
397#endif // DEBUG_INFO
398 }
399 }
400
401
402 Coord n = Norm[pId];
403 Coord nr = restNorm[pId];
404 matK_i(0, 0) += t*nr[0] * nr[0]; matK_i(0, 1) += t*nr[0] * nr[1]; matK_i(0, 2) += t*nr[0] * nr[2];
405 matK_i(1, 0) += t*nr[1] * nr[0]; matK_i(1, 1) += t*nr[1] * nr[1]; matK_i(1, 2) += t*nr[1] * nr[2];
406 matK_i(2, 0) += t*nr[2] * nr[0]; matK_i(2, 1) += t*nr[2] * nr[1]; matK_i(2, 2) += t*nr[2] * nr[2];
407
408
409
410 matL_i(0, 0) += t*n[0] * nr[0]; matL_i(0, 1) += t*n[0] * nr[1]; matL_i(0, 2) += t*n[0] * nr[2];
411 matL_i(1, 0) += t*n[1] * nr[0]; matL_i(1, 1) += t*n[1] * nr[1]; matL_i(1, 2) += t*n[1] * nr[2];
412 matL_i(2, 0) += t*n[2] * nr[0]; matL_i(2, 1) += t*n[2] * nr[1]; matL_i(2, 2) += t*n[2] * nr[2];
413
414 total_weight += 2 * t;
415 if (total_weight > EPSILON)
416 {
417 matL_i *= (1.0f / total_weight);
418 matK_i *= (1.0f / total_weight);
419 }
420
421
422 Matrix R, U, D, V;
423 polarDecomposition(matK_i, R, U, D, V);
424
425#ifdef DEBUG_INFO
426 if (pId == 497)
427 {
428 Matrix mat_out = matK_i;
429 printf("matK_i: \n %f %f %f \n %f %f %f \n %f %f %f \n\n",
430 mat_out(0, 0), mat_out(0, 1), mat_out(0, 2),
431 mat_out(1, 0), mat_out(1, 1), mat_out(1, 2),
432 mat_out(2, 0), mat_out(2, 1), mat_out(2, 2));
433
434 mat_out = matL_i;
435 printf("matL_i: \n %f %f %f \n %f %f %f \n %f %f %f \n\n",
436 mat_out(0, 0), mat_out(0, 1), mat_out(0, 2),
437 mat_out(1, 0), mat_out(1, 1), mat_out(1, 2),
438 mat_out(2, 0), mat_out(2, 1), mat_out(2, 2));
439
440 mat_out = U * D * V.transpose();
441 printf("matK polar: \n %f %f %f \n %f %f %f \n %f %f %f \n\n",
442 mat_out(0, 0), mat_out(0, 1), mat_out(0, 2),
443 mat_out(1, 0), mat_out(1, 1), mat_out(1, 2),
444 mat_out(2, 0), mat_out(2, 1), mat_out(2, 2));
445
446 printf("Horizon: %f; Det: %f \n", horizon, matK_i.determinant());
447 }
448#endif // DEBUG_INFO
449
450 Real maxK = maximum(abs(D(0, 0)), maximum(abs(D(1, 1)),abs(D(2,2))));
451 Real minK = minimum(abs(D(0, 0)), minimum(abs(D(1, 1)),abs(D(2,2))));
452
453
454 bool valid_K = (minK < EPSILON || maxK / minK > Real(1/(strainLimit * strainLimit))) ? false : true;
455 validOfK[pId] = valid_K;
456
457 Matrix F_i;
458 if (valid_K)
459 {
460 invK[pId] = matK_i.inverse();
461 F_i = matL_i * matK_i.inverse();
462
463 }
464 else
465 {
466 Real threshold = 0.0001f * horizon;
467 D(0, 0) = D(0, 0) > threshold ? 1.0 / D(0, 0) : 1.0;
468 D(1, 1) = D(1, 1) > threshold ? 1.0 / D(1, 1) : 1.0;
469 D(2, 2) = D(2, 2) > threshold ? 1.0 / D(2, 2) : 1.0;
470 invK[pId] = V * D * U.transpose();
471 F_i = matL_i * invK[pId];
472 }
473
474 polarDecomposition(F_i, R, U, D, V);
475
476#ifdef DEBUG_INFO
477 if (pId == 497)
478 {
479 Matrix mat_out = F_i;
480 printf("matF_i: \n %f %f %f \n %f %f %f \n %f %f %f \n\n",
481 mat_out(0, 0), mat_out(0, 1), mat_out(0, 2),
482 mat_out(1, 0), mat_out(1, 1), mat_out(1, 2),
483 mat_out(2, 0), mat_out(2, 1), mat_out(2, 2));
484
485 mat_out = D;
486 printf("matD: \n %f %f %f \n %f %f %f \n %f %f %f \n\n",
487 mat_out(0, 0), mat_out(0, 1), mat_out(0, 2),
488 mat_out(1, 0), mat_out(1, 1), mat_out(1, 2),
489 mat_out(2, 0), mat_out(2, 1), mat_out(2, 2));
490
491 printf("Horizon: %f; Det: %f \n", horizon, F_i.determinant());
492 }
493#endif // DEBUG_INFO
494 if (F_i.determinant() < maxDist * EPSILON)
495 {
496 R = Rots[pId];
497 }
498 else
499 Rots[pId] = R;
500
501 //strain limit
502 Coord n0 = Coord(1, 0, 0);
503 Coord n1 = Coord(0, 1, 0);
504 Coord n2 = Coord(0, 0, 1);
505
506#ifdef DEBUG_INFO
507 Coord n0_rot = R * n0;
508 Coord n1_rot = R * n1;
509 Coord n2_rot = R * n2;
510
511 Real un0 = n0_rot.dot(n);
512 Real un1 = n1_rot.dot(n);
513 Real un2 = n2_rot.dot(n);
514 printf("U * N: %f,%f,%f\n",un0,un1,un2);
515
516
517 /*
518 Real l0 = n0_rot.dot(F_i * n0);
519 Real l1 = n1_rot.dot(F_i * n1);
520 Real l2 = n2_rot.dot(F_i * n2);
521 */
522#endif
523
524 matU[pId] = U;
525 matV[pId] = V;
526 Real l0 = D(0, 0);
527 Real l1 = D(1, 1);
528 Real l2 = D(2, 2);
529
530 const Real slimit = min(t, strainLimit);
531 l0 = clamp(l0, slimit, 1/slimit);
532 l1 = clamp(l1, slimit, 1/slimit);
533 l2 = clamp(l2, slimit, 1/slimit);
534 /*
535 if (l0 == slimit)
536 printf("l0 got low clamp, %f\n",l0);
537 else if(l0 == 1/slimit)
538 printf("l0 got high clamp, %f\n",l0);
539 if (l1 == slimit)
540 printf("l1 got low clamp, %f\n",l1);
541 else if (l1 == 1 / slimit)
542 printf("l1 got high clamp, %f\n",l1);
543 if (l2 == slimit)
544 printf("l2 got low clamp, %f\n",l2);
545 else if (l2 == 1 / slimit)
546 printf("l2 got high clamp, %f\n",l2);
547 */
548
549 D(0, 0) = l0;
550 D(1, 1) = l1;
551 D(2, 2) = l2;
552
553 eigens[pId] = Coord(D(0,0), D(1,1), D(2,2));
554 F[pId] = U * D * V.transpose();
555
556 }
557
558
559 template <typename Real, typename Coord, typename Matrix, typename Bond>
560 __global__ void HM_JacobiStepNonsymmetric(
561 DArray<Coord> source,
562 DArray<Matrix> A,
563 DArray<Coord> X,
564 DArray<Coord> y_pre,
565 DArray<Matrix> matU,
566 DArray<Matrix> matV,
567 DArray<Matrix> matR,
568 DArray<Coord> eigen,
569 DArray<bool> validOfK,
570 DArray<Matrix> F,
571 Real k_bend,
572 DArrayList<Bond> bonds,
573 Real horizon,
574 DArray<Real> volume,
575 Real dt,
576 EnergyType type,
577 bool NeighborSearchingAdjacent)
578 {
579 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
580 if (pId >= y_pre.size()) return;
581
582 Coord x_i = X[pId];
583 int size_i = bonds[pId].size();
584
585 Real maxDist = Real(0);
586 for (int ne = 0; ne < size_i; ne++)
587 {
588 Bond bond_ij = bonds[pId][ne];
589 int j = bond_ij.idx;
590 Coord x_j = X[j];
591 Real r = (x_i - x_j).norm();
592
593 maxDist = max(maxDist, r);
594 }
595 maxDist = maxDist < EPSILON ? Real(1) : maxDist;
596
597 Real kappa = 16 / (15 * M_PI * maxDist * maxDist * maxDist * maxDist * maxDist);
598
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 Matrix S1_i;
606 Matrix S2_i;
607 Real vol = volume[pId];
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 == Fiber)
625 {
626 S1_i = ENERGY_FUNC.fiberModel.getStressTensorPositive(lambda_i1, lambda_i2, lambda_i3, V_i);
627 S2_i = ENERGY_FUNC.fiberModel.getStressTensorNegative(lambda_i1, lambda_i2, lambda_i3, V_i);
628 ENERGY_FUNC.fiberModel.getInfo(lambda_i1, lambda_i2, lambda_i3, V_i);
629 }
630
631 Matrix PK1_i = U_i * S1_i * U_i.transpose();
632 Matrix PK2_i = U_i * S2_i * V_i.transpose();
633
634 Real Vol_i = volume[pId];
635
636 Matrix F_i = F[pId];
637
638 Coord y_pre_i = y_pre[pId];
639
640 bool K_valid = validOfK[pId];
641
642 Matrix mat_i(0);
643 Coord source_i(0);
644
645 for (int ne = 0; ne < size_i; ne++)
646 {
647 Bond bond_ij = bonds[pId][ne];
648 int j = bond_ij.idx;
649 Coord y_pre_j = y_pre[j];
650 Coord x_j = X[j];
651 Real r = (x_j - x_i).norm();
652
653 if (r > EPSILON)
654 {
655 Real weight_ij = 1.;
656
657 Real Vol_j = volume[j];
658
659 Coord y_ij = y_pre_i - y_pre_j;
660
661
662 Real lambda = y_ij.norm() / r;
663
664 const Real scale = Vol_i * Vol_j * kappa;
665 const Real sw_ij = dt * dt * scale * weight_ij;
666
667 bool K_valid_ij = K_valid & validOfK[j];
668
669 Matrix S1_ij;
670 Matrix S2_ij;
671 if (type == StVK) {
672 S1_ij = ENERGY_FUNC.stvkModel.getStressTensorPositive(lambda, lambda, lambda);
673 S2_ij = ENERGY_FUNC.stvkModel.getStressTensorNegative(lambda, lambda, lambda);
674 }
675 else if (type == NeoHooekean) {
676 S1_ij = ENERGY_FUNC.neohookeanModel.getStressTensorPositive(lambda, lambda, lambda);
677 S2_ij = ENERGY_FUNC.neohookeanModel.getStressTensorNegative(lambda, lambda, lambda);
678 }
679 else if (type == Linear) {
680 S1_ij = ENERGY_FUNC.linearModel.getStressTensorPositive(lambda, lambda, lambda);
681 S2_ij = ENERGY_FUNC.linearModel.getStressTensorNegative(lambda, lambda, lambda);
682 }
683 else if (type == Xuetal) {
684 S1_ij = ENERGY_FUNC.xuModel.getStressTensorPositive(lambda, lambda, lambda);
685 S2_ij = ENERGY_FUNC.xuModel.getStressTensorNegative(lambda, lambda, lambda);
686 }
687 else if (type == Fiber) {
688 S1_ij = ENERGY_FUNC.fiberModel.getStressTensorPositive(lambda, lambda,lambda, V_i);
689 S2_ij = ENERGY_FUNC.fiberModel.getStressTensorNegative(lambda, lambda,lambda, V_i);
690 }
691
692 Matrix PK1_ij = K_valid_ij ? PK1_i : S1_ij;
693 Matrix PK2_ij = K_valid_ij ? PK2_i : S2_ij;
694
695 if (NeighborSearchingAdjacent)
696 {
697 Real standardVol = 0.0025f * 0.005f;
698
699 PK1_ij *= (x_i - x_j).normSquared() / standardVol;
700 PK2_ij *= (x_i - x_j).normSquared() / standardVol;
701 }
702
703 Coord dir_ij = lambda > EPSILON ? y_ij.normalize() : Coord(1, 0, 0);
704
705 Coord x_ij = K_valid_ij ? x_i - x_j : dir_ij * (x_i - x_j).norm();
706
707 Matrix F_i_1 = F_i.inverse();
708 Matrix F_i_T = F_i_1.transpose();
709
710 Matrix F_TF_1 = k_bend * F_i_T * F_i_1;
711
712 Coord ds_ij = sw_ij * (PK1_ij + F_TF_1) * y_pre_j + sw_ij * (PK2_ij + k_bend * F_i_1) * x_ij;
713 Coord ds_ji = sw_ij * (PK1_ij + F_TF_1) * y_pre_i - sw_ij * (PK2_ij + k_bend * F_i_1) * x_ij;
714
715 Matrix mat_ij = sw_ij * (PK1_ij + F_TF_1);
716
717
718 source_i += ds_ij;
719
720 mat_i += (mat_ij);
721
722
723 atomicAdd(&source[j][0], ds_ji[0]);
724 atomicAdd(&source[j][1], ds_ji[1]);
725 atomicAdd(&source[j][2], ds_ji[2]);
726
727 atomicAdd(&A[j](0, 0), mat_ij(0, 0));
728 atomicAdd(&A[j](0, 1), mat_ij(0, 1));
729 atomicAdd(&A[j](0, 2), mat_ij(0, 2));
730 atomicAdd(&A[j](1, 0), mat_ij(1, 0));
731 atomicAdd(&A[j](1, 1), mat_ij(1, 1));
732 atomicAdd(&A[j](1, 2), mat_ij(1, 2));
733 atomicAdd(&A[j](2, 0), mat_ij(2, 0));
734 atomicAdd(&A[j](2, 1), mat_ij(2, 1));
735 atomicAdd(&A[j](2, 2), mat_ij(2, 2));
736 }
737 }
738
739 atomicAdd(&source[pId][0], source_i[0]);
740 atomicAdd(&source[pId][1], source_i[1]);
741 atomicAdd(&source[pId][2], source_i[2]);
742
743 atomicAdd(&A[pId](0, 0), mat_i(0, 0));
744 atomicAdd(&A[pId](0, 1), mat_i(0, 1));
745 atomicAdd(&A[pId](0, 2), mat_i(0, 2));
746 atomicAdd(&A[pId](1, 0), mat_i(1, 0));
747 atomicAdd(&A[pId](1, 1), mat_i(1, 1));
748 atomicAdd(&A[pId](1, 2), mat_i(1, 2));
749 atomicAdd(&A[pId](2, 0), mat_i(2, 0));
750 atomicAdd(&A[pId](2, 1), mat_i(2, 1));
751 atomicAdd(&A[pId](2, 2), mat_i(2, 2));
752 }
753
754
755
756 template <typename Matrix>
757 __global__ void HM_InitRotation(
758 DArray<Matrix> U,
759 DArray<Matrix> V,
760 DArray<Matrix> R)
761 {
762 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
763 if (pId >= U.size()) return;
764
765 U[pId] = Matrix::identityMatrix();
766 V[pId] = Matrix::identityMatrix();
767 R[pId] = Matrix::identityMatrix();
768 }
769
770
771
772 template<typename Matrix, typename Coord, typename Real>
773 __global__ void HM_ComputeNextPosition(
774 DArray<Coord> y_next,
775 DArray<Coord> y_inter,
776 DArray<Real> volume,
777 DArray<Coord> source,
778 DArray<Matrix> A
779 )
780 {
781 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
782 if (pId >= y_next.size()) return;
783
784 Real mass_i = volume[pId] * 1000.0;
785 Matrix mat_i = A[pId] + mass_i * Matrix::identityMatrix();
786 Coord src_i = source[pId] + mass_i * y_inter[pId];
787 y_next[pId] = mat_i.inverse() * src_i;
788 }
789
790
791 template <typename Real>
792 __global__ void HM_InitVolume(
793 DArray<Real> volume,
794 Real objectVolume,
795 bool objectVolumeSet,
796 Real particleVolume,
797 bool particleVolumeSet
798 )
799 {
800 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
801 if (pId >= volume.size()) return;
802
803 if (objectVolumeSet && volume.size() != 0) volume[pId] = Real(objectVolume / volume.size());
804 else if (particleVolumeSet) volume[pId] = particleVolume;
805 else volume[pId] = 0.001;
806 }
807
808
809 template <typename Coord, typename Real>
810 __global__ void HM_ComputeGradient(
811 DArray<Coord> grad,
812 DArray<Real> grad_m,
813 DArray<Coord> y_pre,
814 DArray<Coord> y_next)
815 {
816 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
817 if (pId >= y_next.size()) return;
818
819 grad[pId] = y_next[pId] - y_pre[pId];
820 grad_m[pId] = sqrt(grad[pId].dot(grad[pId]));
821 }
822
823
824 template<typename TDataType>
825 void CoSemiImplicitHyperelasticitySolver<TDataType>::resizeAllFields()
826 {
827
828 uint num = this->inY()->getData().size();
829
830 if (m_F.size() == num)
831 return;
832
833 m_F.resize(num);
834
835 m_invF.resize(num);
836 m_invK.resize(num);
837 m_validOfK.resize(num);
838 m_validOfF.resize(num);
839
840 m_eigenValues.resize(num);
841 m_matU.resize(num);
842 m_matV.resize(num);
843 m_matR.resize(num);
844 m_volume.resize(num);
845
846 m_energy.resize(num);
847 m_alpha.resize(num);
848 m_gradient.resize(num);
849 mEnergyGradient.resize(num);
850
851 y_pre.resize(num);
852 y_next.resize(num);
853 y_current.resize(num);
854 y_gradC.resize(num);
855
856 m_source.resize(num);
857 m_A.resize(num);
858 m_gradientMagnitude.resize(num);
859 this->mPosBuf.resize(num);
860
861 m_fraction.resize(num);
862
863 m_bFixed.resize(num);
864 m_fixedPos.resize(num);
865 mPosBuf_March.resize(num);
866
867 initializeVolume();
868
869 cuExecute(m_matU.size(),
870 HM_InitRotation,
871 m_matU,
872 m_matV,
873 m_matR);
874
875 m_reduce = Reduction<Real>::Create(num);
876
877 auto triSet = TypeInfo::cast<TriangleSet<TDataType>>(this->inTriangularMesh()->getDataPtr());
878 triSet->getPoints().assign(this->inY()->getData());
879 triSet->updateAngleWeightedVertexNormal(this->inNorm()->getData());
880 }
881
882 template<typename TDataType>
883 void CoSemiImplicitHyperelasticitySolver<TDataType>::enforceHyperelasticity()
884 {
885 resizeAllFields();
886
887 int numOfParticles = this->inY()->getData().size();
888 uint pDims = cudaGridSize(numOfParticles, BLOCK_SIZE);
889
890 std::cout << "enforceElasticity Particles: " << numOfParticles << std::endl;
891
892 this->mWeights.reset();
893
894 Log::sendMessage(Log::User, "\n \n \n \n=================== solver start!!!===================");
895
896
897
898 /*====================================== Jacobi method ======================================*/
899 // initialize y_now, y_next_iter
900 y_current.assign(this->inY()->getData());
901 this->mPosBuf.assign(this->inY()->getData());
902
903
904 // do Jacobi method Loop
905 bool convergeFlag = false; // converge or not
906 int iterCount = 0;
907 Real omega;
908 Real alpha = 1.0f;
909 Real max_grad_mag = Real(1000.0);
910
911 if (true) {
912 while (iterCount < this->varIterationNumber()->getData() && max_grad_mag >= this->grad_res_eps) {
913
914 m_source.reset();
915 m_A.reset();
916
917 HM_ComputeF << <pDims, BLOCK_SIZE >> > (
918 m_F,
919 m_eigenValues,
920 m_invK,
921 m_validOfK,
922 m_matU,
923 m_matV,
924 m_matR,
925 this->inX()->getData(),
926 y_current,
927 this->inBonds()->getData(),
928 this->inHorizon()->getData(),
929 (Real const)0.3,
930 this->inRestNorm()->getData(),
931 this->inNorm()->getData());
932 cuSynchronize();
933
934 HM_JacobiStepNonsymmetric << <pDims, BLOCK_SIZE >> > (
935 m_source,
936 m_A,
937 this->inX()->getData(),
938 y_current,
939 m_matU,
940 m_matV,
941 m_matR,
942 m_eigenValues,
943 m_validOfK,
944 m_F,
945 this->k_bend,
946 this->inBonds()->getData(),
947 this->inHorizon()->getData(),
948 m_volume,
949 this->inTimeStep()->getData(),
950 this->inEnergyType()->getData(),
951 this->varNeighborSearchingAdjacent()->getData());
952 cuSynchronize();
953
954 cuExecute(y_current.size(),
955 HM_ComputeNextPosition,
956 y_next,
957 this->mPosBuf,
958 m_volume,
959 m_source,
960 m_A);
961
962 //sub steping
963 cuExecute(m_gradient.size(),
964 HM_ComputeGradient,
965 m_gradient,
966 m_gradientMagnitude,
967 y_current,
968 y_next);
969
970
971 Reduction<Real> reduce;
972 max_grad_mag = reduce.maximum(m_gradientMagnitude.begin(), m_gradientMagnitude.size());
973
974 if (this->m_alphaCompute) {
975 cuExecute(m_energy.size(),
976 HM_Compute1DEnergy,
977 m_energy,
978 mEnergyGradient,
979 this->inX()->getData(),
980 y_current,
981 m_F,
982 m_volume,
983 m_validOfK,
984 m_eigenValues,
985 this->inBonds()->getData(),
986 this->inEnergyType()->getData());
987
988 cuExecute(m_alpha.size(),
989 HM_ComputeStepLength,
990 m_alpha,
991 m_gradient,
992 mEnergyGradient,
993 m_volume,
994 m_A,
995 m_energy,
996 this->inBonds()->getData());
997
998 cuExecute(m_gradient.size(),
999 HM_ComputeCurrentPosition,
1000 m_gradient,
1001 y_current,
1002 y_next,
1003 m_alpha);
1004 }
1005 else {
1006 alpha = Real(1.0);
1007 cuExecute(m_gradient.size(),
1008 HM_ComputeCurrentPosition,
1009 m_gradient,
1010 y_current,
1011 y_next,
1012 alpha);
1013 }
1014
1015 iterCount++;
1016
1017 y_current.assign(y_next);
1018 }
1019
1020 // do Jacobi method Loop
1021 mContactRule->inNewPosition()->assign(y_next);
1022 mContactRule->initCCDBroadPhase();
1023
1024 mContactRule->inNewPosition()->assign(this->mPosBuf);
1025 convergeFlag = false; // converge or not
1026 iterCount = 0;
1027 alpha = 1.0f;
1028 max_grad_mag = 1e3;
1029 }
1030 mPosBuf_March.assign(this->mPosBuf);
1031
1032
1033 auto& cntPos = mContactRule->inNewPosition()->getData();
1034 while (iterCount < this->varIterationNumber()->getData() && max_grad_mag >= this->grad_res_eps) {
1035 m_source.reset();
1036 m_A.reset();
1037 y_current.assign(cntPos);
1038
1039 HM_ComputeF << <pDims, BLOCK_SIZE >> > (
1040 m_F,
1041 m_eigenValues,
1042 m_invK,
1043 m_validOfK,
1044 m_matU,
1045 m_matV,
1046 m_matR,
1047 this->inX()->getData(),
1048 y_current,
1049 this->inBonds()->getData(),
1050 this->inHorizon()->getData(),
1051 (Real const)0.3,
1052 this->inRestNorm()->getData(),
1053 this->inNorm()->getData());
1054 cuSynchronize();
1055
1056 HM_JacobiStepNonsymmetric << <pDims, BLOCK_SIZE >> > (
1057 m_source,
1058 m_A,
1059 this->inX()->getData(),
1060 y_current,
1061 m_matU,
1062 m_matV,
1063 m_matR,
1064 m_eigenValues,
1065 m_validOfK,
1066 m_F,
1067 this->k_bend,
1068 this->inBonds()->getData(),
1069 this->inHorizon()->getData(),
1070 m_volume,
1071 this->inTimeStep()->getData(),
1072 this->inEnergyType()->getData(),
1073 this->varNeighborSearchingAdjacent()->getData());
1074 cuSynchronize();
1075
1076 cuExecute(cntPos.size(),
1077 HM_ComputeNextPosition,
1078 cntPos,
1079 mPosBuf_March,
1080 m_volume,
1081 m_source,
1082 m_A);
1083
1084 //sub steping
1085 cuExecute(m_gradient.size(),
1086 HM_ComputeGradient,
1087 m_gradient,
1088 m_gradientMagnitude,
1089 y_current,
1090 cntPos);
1091
1092 Reduction<Real> reduce;
1093 max_grad_mag = reduce.maximum(m_gradientMagnitude.begin(), m_gradientMagnitude.size());
1094
1095 if (this->m_alphaCompute) {
1096 cuExecute(m_energy.size(),
1097 HM_Compute1DEnergy,
1098 m_energy,
1099 mEnergyGradient,
1100 this->inX()->getData(),
1101 y_current,
1102 m_F,
1103 m_volume,
1104 m_validOfK,
1105 m_eigenValues,
1106 this->inBonds()->getData(),
1107 this->inEnergyType()->getData());
1108
1109 cuExecute(m_alpha.size(),
1110 HM_ComputeStepLength,
1111 m_alpha,
1112 m_gradient,
1113 mEnergyGradient,
1114 m_volume,
1115 m_A,
1116 m_energy,
1117 this->inBonds()->getData());
1118
1119 cuExecute(m_gradient.size(),
1120 HM_ComputeCurrentPosition,
1121 m_gradient,
1122 y_current,
1123 cntPos,
1124 m_alpha);
1125 }
1126 else {
1127 alpha = Real(1.0);
1128
1129 cuExecute(m_gradient.size(),
1130 HM_ComputeCurrentPosition,
1131 m_gradient,
1132 y_current,
1133 cntPos,
1134 alpha);
1135 }
1136 if (this->selfContact) {
1137 if (this->acc) {
1138 mContactRule->update();
1139 }
1140 else {
1141 mContactRule->constrain();
1142 }
1143 }
1144
1145 iterCount++;
1146 }
1147
1148 /*========================= end of alg, marching time step==============================*/
1149 printf("========= Enforcement elastic run %d iteration =======\n", iterCount);
1150
1151
1152 cuExecute(this->inY()->getDataPtr()->size(),
1153 test_HM_UpdatePosition,
1154 this->inY()->getData(),
1155 this->inVelocity()->getData(),
1156 cntPos,
1157 this->mPosBuf,
1158 this->inAttribute()->getData(),
1159 this->inTimeStep()->getData());
1160
1161
1162 Log::sendMessage(Log::User, "\n==================== solver end!!!====================\n");
1163 }
1164
1165 DEFINE_CLASS(CoSemiImplicitHyperelasticitySolver);
1166}
Generated by doxygen 1.13.2