PeriDyno 1.0.0
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
ElastoplasticityModule.cu
Go to the documentation of this file.
1#include "ElastoplasticityModule.h"
2#include "Matrix/MatrixFunc.h"
3
4namespace dyno
5{
6 IMPLEMENT_TCLASS(ElastoplasticityModule, TDataType)
7
8 template<typename TDataType>
9 ElastoplasticityModule<TDataType>::ElastoplasticityModule()
10 : LinearElasticitySolver<TDataType>()
11 {
12 this->varCohesion()->setRange(0.0f, 1.0f);
13 this->varFrictionAngle()->setRange(0.0f, 1.0f);
14
15 mDensityPBD = std::make_shared<IterativeDensitySolver<TDataType>>();
16 mDensityPBD->varIterationNumber()->setValue(1);
17 this->inTimeStep()->connect(mDensityPBD->inTimeStep());
18 this->inHorizon()->connect(mDensityPBD->inSmoothingLength());
19 this->inY()->connect(mDensityPBD->inPosition());
20 this->inVelocity()->connect(mDensityPBD->inVelocity());
21 this->inNeighborIds()->connect(mDensityPBD->inNeighborIds());
22 }
23
24 template<typename TDataType>
25 ElastoplasticityModule<TDataType>::~ElastoplasticityModule()
26 {
27 m_bYield.clear();
28 m_invF.clear();
29 m_yiled_I1.clear();
30 m_yield_J2.clear();
31 m_I1.clear();
32 }
33
34 __device__ Real Hardening(Real rho)
35 {
36 Real hardening = 1.0f;
37
38 return 1.0f;
39
40 if (rho >= 1000)
41 {
42 Real ratio = rho / 1000;
43 return pow((float)M_E, hardening*(ratio - 1.0f));
44 }
45 else
46 return 1.0f;
47 }
48
49 template <typename Real, typename Coord, typename Matrix, typename Bond>
50 __global__ void PM_ComputeInvariants(
51 DArray<bool> bYield,
52 DArray<Real> yield_I1,
53 DArray<Real> yield_J2,
54 DArray<Real> arrI1,
55 DArray<Coord> X,
56 DArray<Coord> Y,
57 DArray<Real> bulk_stiffiness,
58 DArrayList<Bond> bonds,
59 Real horizon,
60 Real A,
61 Real B,
62 Real mu,
63 Real lambda)
64 {
65 int i = threadIdx.x + (blockIdx.x * blockDim.x);
66 if (i >= Y.size()) return;
67
68 CorrectedKernel<Real> kernSmooth;
69
70 Real weaking = 1.0f;// Softening(rhoArr[i]);
71
72 Real s_A = weaking*A;
73
74 List<Bond>& bonds_i = bonds[i];
75 Coord x_i = X[i];
76 Coord y_i = Y[i];
77
78 Real I1_i = 0.0f;
79 Real J2_i = 0.0f;
80 //compute the first and second invariants of the deformation state, i.e., I1 and J2
81 int size_i = bonds_i.size();
82 Real total_weight = Real(0);
83 // Compute e^{iso} in Equation (16)
84 for (int ne = 0; ne < size_i; ne++)
85 {
86 Bond bond_ij = bonds_i[ne];
87 int j = bond_ij.idx;
88 Coord x_j = X[j];
89
90 Real r = (x_i - x_j).norm();
91
92 if (r > 0.01*horizon)
93 {
94 Real weight = kernSmooth.Weight(r, horizon);
95 Coord p = (Y[j] - y_i);
96 Real ratio_ij = p.norm() / r;
97
98 I1_i += weight*ratio_ij;
99
100 total_weight += weight;
101 }
102 }
103
104 if (total_weight > EPSILON)
105 {
106 I1_i /= total_weight;
107 }
108 else
109 {
110 I1_i = 1.0f;
111 }
112
113 // Compute e^{dev} in Equation (16)
114 for (int ne = 0; ne < size_i; ne++)
115 {
116 Bond bond_ij = bonds_i[ne];
117 int j = bond_ij.idx;
118 Coord x_j = X[j];
119 Real r = (x_i - x_j).norm();
120
121 if (r > 0.01*horizon)
122 {
123 Real weight = kernSmooth.Weight(r, horizon);
124 Vec3f p = (Y[j] - y_i);
125 Real ratio_ij = p.norm() / r;
126 J2_i = (ratio_ij - I1_i)*(ratio_ij - I1_i)*weight;
127 }
128 }
129 if (total_weight > EPSILON)
130 {
131 J2_i /= total_weight;
132 J2_i = sqrt(J2_i);
133 }
134 else
135 {
136 J2_i = 0.0f;
137 }
138
139 Real D1 = 1 - I1_i; //positive for compression and negative for stretching
140
141 Real yield_I1_i = 0.0f;
142 Real yield_J2_i = 0.0f;
143
144 //Equation (17)
145 Real s_J2 = J2_i*mu*bulk_stiffiness[i];
146 Real s_D1 = D1*lambda*bulk_stiffiness[i];
147
148 //Drucker-Prager yield criterion not violated
149 if (s_J2 <= s_A + B*s_D1)
150 {
151 yield_I1[i] = Real(0);
152 yield_J2[i] = Real(0);
153
154 bYield[i] = false;
155 }
156 else //Equation (18)
157 {
158 //bulk_stiffiness[i] = 0.0f;
159 if (s_A + B*s_D1 > 0.0f)
160 {
161 yield_I1_i = 0.0f;
162
163 yield_J2_i = (s_J2 - (s_A + B*s_D1)) / s_J2;
164 }
165 else
166 {
167 yield_I1_i = 1.0f;
168 if (s_A + B*s_D1 < -EPSILON)
169 {
170 yield_I1_i = (s_A + B*s_D1) / (B*s_D1);
171 }
172
173 yield_J2_i = 1.0f;
174 }
175
176 yield_I1[i] = yield_I1_i;
177 yield_J2[i] = yield_J2_i;
178 }
179
180 arrI1[i] = I1_i;
181 }
182
183 template <typename Real, typename Coord, typename Matrix, typename Bond>
184 __global__ void PM_ApplyYielding(
185 DArray<Real> yield_I1,
186 DArray<Real> yield_J2,
187 DArray<Real> arrI1,
188 DArray<Coord> X,
189 DArray<Coord> Y,
190 DArrayList<Bond> bonds)
191 {
192 int i = threadIdx.x + (blockIdx.x * blockDim.x);
193 if (i >= Y.size()) return;
194
195 List<Bond>& bonds_i = bonds[i];
196 Coord x_i = X[i];
197 Coord y_i = Y[i];
198
199 Real yield_I1_i = yield_I1[i];
200 Real yield_J2_i = yield_J2[i];
201 Real I1_i = arrI1[i];
202
203 //add permanent deformation
204 int size_i = bonds_i.size();
205 for (int ne = 0; ne < size_i; ne++)
206 {
207 Bond bond_ij = bonds_i[ne];
208 int j = bond_ij.idx;
209 Coord x_j = X[j];
210
211
212 Real yield_I1_j = yield_I1[j];
213 Real yield_J2_j = yield_J2[j];
214 Real I1_j = arrI1[j];
215
216 Real r = (x_i - x_j).norm();
217
218 Coord p = (Y[j] - y_i);
219 Coord q = (x_j - x_i);
220
221 //Coord new_q = q*I1_i;
222 Coord new_q = q*(I1_i + I1_j) / 2;
223 Coord D_iso = new_q - q;
224
225 Coord dir_q = q;
226 dir_q = dir_q.norm() > EPSILON ? dir_q.normalize() : Coord(0);
227
228 Coord D_dev = p.norm()*dir_q - new_q;
229 //Coord D_dev = p - new_q;
230
231 Bond newBond_ij;
232
233 //Coord new_rest_pos_j = rest_pos_j + yield_I1_i * D_iso + yield_J2_i * D_dev;
234 Coord new_rest_pos_j = x_j + (yield_I1_i + yield_I1_j) / 2 * D_iso + (yield_J2_i + yield_J2_j) / 2 * D_dev;
235
236 newBond_ij.xi = new_rest_pos_j - x_i;
237 newBond_ij.idx = j;
238 bonds_i[ne] = newBond_ij;
239 }
240
241 }
242
243 // int iter = 0;
244 template<typename TDataType>
245 void ElastoplasticityModule<TDataType>::constrain()
246 {
247 if (m_invF.size() != this->inY()->size())
248 {
249 m_invF.resize(this->inY()->size());
250 m_yiled_I1.resize(this->inY()->size());
251 m_yield_J2.resize(this->inY()->size());
252 m_I1.resize(this->inY()->size());
253 m_bYield.resize(this->inY()->size());
254
255 m_bYield.reset();
256 }
257
258 this->solveElasticity();
259 this->applyPlasticity();
260 }
261
262
263 template<typename TDataType>
264 void ElastoplasticityModule<TDataType>::solveElasticity()
265 {
266 this->mPosBuf.assign(this->inY()->getData());
267
268 this->computeInverseK();
269
270 int iter = 0;
271 int total = this->varIterationNumber()->getData();
272 while (iter < total)
273 {
274 this->enforceElasticity();
275 if (this->varIncompressible()->getValue() == true) {
276 mDensityPBD->update();
277 }
278
279 iter++;
280 }
281
282 this->updateVelocity();
283 }
284
285 template<typename TDataType>
286 void ElastoplasticityModule<TDataType>::applyPlasticity()
287 {
288 this->rotateRestShape();
289
290 this->computeMaterialStiffness();
291 this->applyYielding();
292
293 this->reconstructRestShape();
294 }
295
296 template<typename TDataType>
297 void ElastoplasticityModule<TDataType>::applyYielding()
298 {
299 int num = this->inY()->size();
300 uint pDims = cudaGridSize(num, BLOCK_SIZE);
301
302 Real A = computeA();
303 Real B = computeB();
304
305 PM_ComputeInvariants<Real, Coord, Matrix, Bond> << <pDims, BLOCK_SIZE >> > (
306 m_bYield,
307 m_yiled_I1,
308 m_yield_J2,
309 m_I1,
310 this->inX()->getData(),
311 this->inY()->getData(),
312 this->mBulkStiffness,
313 this->inBonds()->getData(),
314 this->inHorizon()->getData(),
315 A,
316 B,
317 this->varMu()->getData(),
318 this->varLambda()->getData());
319 cuSynchronize();
320 //
321 PM_ApplyYielding<Real, Coord, Matrix, Bond> << <pDims, BLOCK_SIZE >> > (
322 m_yiled_I1,
323 m_yield_J2,
324 m_I1,
325 this->inX()->getData(),
326 this->inY()->getData(),
327 this->inBonds()->getData());
328 cuSynchronize();
329 }
330
331 template <typename Real, typename Coord, typename Matrix, typename Bond>
332 __global__ void PM_ReconstructRestShape(
333 DArrayList<Bond> newBonds,
334 DArray<bool> bYield,
335 DArray<Coord> X,
336 DArray<Coord> Y,
337 DArray<Real> I1,
338 DArray<Real> I1_yield,
339 DArray<Real> J2_yield,
340 DArray<Matrix> invF,
341 DArrayList<int> neighborhood,
342 DArrayList<Bond> bonds,
343 Real horizon)
344 {
345 int i = threadIdx.x + (blockIdx.x * blockDim.x);
346 if (i >= newBonds.size()) return;
347
348 List<int>& list_i = neighborhood[i];
349 List<Bond>& bonds_i = bonds[i];
350 List<Bond>& newBonds_i = newBonds[i];
351
352 // update neighbors
353 if (!bYield[i])
354 {
355 int new_size = bonds_i.size();
356 for (int ne = 0; ne < new_size; ne++)
357 {
358 Bond pair = bonds_i[ne];
359 newBonds_i.insert(pair);
360 }
361 }
362 else
363 {
364 int nbSize = list_i.size();
365 Coord y_i = Y[i];
366
367 Matrix invF_i = invF[i];
368
369 Bond np;
370 for (int ne = 0; ne < nbSize; ne++)
371 {
372 int j = list_i[ne];
373 Matrix invF_j = invF[j];
374
375 np.idx = j;
376 np.xi = 0.5*(invF_i + invF_j)*(Y[j] - y_i);
377
378 newBonds_i.insert(np);
379
380// if (i == j)
381// {
382// Bond np_0 = new_list_np_i[0];
383// new_list_np_i[0] = np;
384// new_list_np_i[ne] = np_0;
385// }
386 }
387 }
388
389 bYield[i] = false;
390 }
391
392 template <typename Bond>
393 __global__ void PM_ReconfigureRestShape(
394 DArray<uint> nbSize,
395 DArray<bool> bYield,
396 DArrayList<int> neighborhood,
397 DArrayList<Bond> restShape)
398 {
399 int i = threadIdx.x + (blockIdx.x * blockDim.x);
400 if (i >= nbSize.size()) return;
401
402 if (bYield[i]) {
403 nbSize[i] = neighborhood[i].size();
404 }
405 else {
406 nbSize[i] = restShape[i].size();
407 }
408 }
409
410 template <typename Real, typename Coord, typename Matrix, typename Bond>
411 __global__ void PM_ComputeInverseDeformation(
412 DArray<Matrix> invF,
413 DArray<Coord> X,
414 DArray<Coord> Y,
415 DArrayList<Bond> bonds,
416 Real horizon)
417 {
418 int i = threadIdx.x + (blockIdx.x * blockDim.x);
419 if (i >= invF.size()) return;
420
421 CorrectedKernel<Real> kernSmooth;
422
423 //reconstruct the rest shape as the yielding condition is violated.
424 Real total_weight = 0.0f;
425 Matrix curM(0);
426 Matrix refM(0);
427
428 List<Bond>& bonds_i = bonds[i];
429 Coord x_i = X[i];
430 int size_i = bonds_i.size();
431 for (int ne = 0; ne < size_i; ne++)
432 {
433 Bond bond_ij = bonds_i[ne];
434 int j = bond_ij.idx;
435 Coord x_j = X[j];
436 Real r = (x_j - x_i).norm();
437
438 if (r > EPSILON)
439 {
440 Real weight = kernSmooth.Weight(r, horizon);
441
442 Coord p = (Y[j] - Y[i]) / horizon;
443 Coord q = (x_j - x_i) / horizon;
444
445 curM(0, 0) += p[0] * p[0] * weight; curM(0, 1) += p[0] * p[1] * weight; curM(0, 2) += p[0] * p[2] * weight;
446 curM(1, 0) += p[1] * p[0] * weight; curM(1, 1) += p[1] * p[1] * weight; curM(1, 2) += p[1] * p[2] * weight;
447 curM(2, 0) += p[2] * p[0] * weight; curM(2, 1) += p[2] * p[1] * weight; curM(2, 2) += p[2] * p[2] * weight;
448
449 refM(0, 0) += q[0] * p[0] * weight; refM(0, 1) += q[0] * p[1] * weight; refM(0, 2) += q[0] * p[2] * weight;
450 refM(1, 0) += q[1] * p[0] * weight; refM(1, 1) += q[1] * p[1] * weight; refM(1, 2) += q[1] * p[2] * weight;
451 refM(2, 0) += q[2] * p[0] * weight; refM(2, 1) += q[2] * p[1] * weight; refM(2, 2) += q[2] * p[2] * weight;
452
453 total_weight += weight;
454 }
455 }
456
457
458 if (total_weight < EPSILON)
459 {
460 total_weight = Real(1);
461 }
462 refM *= (1.0f / total_weight);
463 curM *= (1.0f / total_weight);
464
465 Real threshold = Real(0.00001);
466 Matrix curR, curU, curD, curV;
467
468 polarDecomposition(curM, curR, curU, curD, curV);
469
470 curD(0, 0) = curD(0, 0) > threshold ? 1.0 / curD(0, 0) : 1.0 / threshold;
471 curD(1, 1) = curD(1, 1) > threshold ? 1.0 / curD(1, 1) : 1.0 / threshold;
472 curD(2, 2) = curD(2, 2) > threshold ? 1.0 / curD(2, 2) : 1.0 / threshold;
473 refM *= curV*curD*curU.transpose();
474
475 // if (abs(refM.determinant() - 1) > 0.05f)
476 // {
477 // refM = Matrix::identityMatrix();
478 // }
479
480 if (refM.determinant() < EPSILON)
481 {
482 refM = Matrix::identityMatrix();
483 }
484
485 invF[i] = refM;
486 }
487
488 __global__ void PM_EnableAllReconstruction(
489 DArray<bool> bYield)
490 {
491 int i = threadIdx.x + (blockIdx.x * blockDim.x);
492 if (i >= bYield.size()) return;
493
494 bYield[i] = true;
495 }
496
497 template<typename TDataType>
498 void ElastoplasticityModule<TDataType>::reconstructRestShape()
499 {
500 //constructRestShape(m_neighborhood.getData(), m_position.getData());
501
502 auto& pts = this->inY()->getData();
503
504 if (this->varRenewNeighborhood()->getValue())
505 {
506 cuExecute(m_bYield.size(),
507 PM_EnableAllReconstruction,
508 m_bYield);
509 }
510
511 DArray<uint> index(this->inY()->getDataPtr()->size());
512
513 cuExecute(index.size(),
514 PM_ReconfigureRestShape,
515 index,
516 m_bYield,
517 this->inNeighborIds()->getData(),
518 this->inBonds()->getData());
519
520 DArrayList<Bond> newBonds;
521 newBonds.resize(index);
522
523 cuExecute(m_invF.size(),
524 PM_ComputeInverseDeformation,
525 m_invF,
526 this->inX()->getData(),
527 this->inY()->getData(),
528 this->inBonds()->getData(),
529 this->inHorizon()->getData());
530
531 cuExecute(newBonds.size(),
532 PM_ReconstructRestShape,
533 newBonds,
534 m_bYield,
535 this->inX()->getData(),
536 this->inY()->getData(),
537 m_I1,
538 m_yiled_I1,
539 m_yield_J2,
540 m_invF,
541 this->inNeighborIds()->getData(),
542 this->inBonds()->getData(),
543 this->inHorizon()->getData());
544
545 this->inBonds()->assign(newBonds);
546
547 newBonds.clear();
548 index.clear();
549 cuSynchronize();
550 }
551
552 template <typename Real, typename Coord, typename Matrix, typename Bond>
553 __global__ void EM_RotateRestShape(
554 DArray<Coord> X,
555 DArray<Coord> Y,
556 DArray<bool> bYield,
557 DArrayList<Bond> bonds,
558 Real smoothingLength)
559 {
560 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
561 if (pId >= Y.size()) return;
562
563 SmoothKernel<Real> kernSmooth;
564
565 List<Bond>& bonds_i = bonds[pId];
566 Coord x_i = X[pId];
567 int size_i = bonds_i.size();
568
569 // cout << i << " " << rids[shape_i.ids[shape_i.idx]] << endl;
570 Real total_weight = 0.0f;
571 Matrix mat_i(0);
572 Matrix invK_i(0);
573 for (int ne = 0; ne < size_i; ne++)
574 {
575 Bond bond_ij = bonds_i[ne];
576 int j = bond_ij.idx;
577 Coord x_j = X[j];
578 Real r = (x_i - x_j).norm();
579
580 if (r > EPSILON)
581 {
582 Real weight = kernSmooth.Weight(r, smoothingLength);
583
584 Coord p = (Y[j] - Y[pId]) / smoothingLength;
585 //Vec3f q = (shape_i.pos[ne] - rest_i)*(1.0f/r)*weight;
586 Coord q = (x_j - x_i) / smoothingLength;
587
588 mat_i(0, 0) += p[0] * q[0] * weight; mat_i(0, 1) += p[0] * q[1] * weight; mat_i(0, 2) += p[0] * q[2] * weight;
589 mat_i(1, 0) += p[1] * q[0] * weight; mat_i(1, 1) += p[1] * q[1] * weight; mat_i(1, 2) += p[1] * q[2] * weight;
590 mat_i(2, 0) += p[2] * q[0] * weight; mat_i(2, 1) += p[2] * q[1] * weight; mat_i(2, 2) += p[2] * q[2] * weight;
591
592 invK_i(0, 0) += q[0] * q[0] * weight; invK_i(0, 1) += q[0] * q[1] * weight; invK_i(0, 2) += q[0] * q[2] * weight;
593 invK_i(1, 0) += q[1] * q[0] * weight; invK_i(1, 1) += q[1] * q[1] * weight; invK_i(1, 2) += q[1] * q[2] * weight;
594 invK_i(2, 0) += q[2] * q[0] * weight; invK_i(2, 1) += q[2] * q[1] * weight; invK_i(2, 2) += q[2] * q[2] * weight;
595
596 total_weight += weight;
597 }
598 }
599
600 if (total_weight > EPSILON)
601 {
602 mat_i *= (1.0f / total_weight);
603 invK_i *= (1.0f / total_weight);
604 }
605
606 Matrix R, U, D, V;
607 polarDecomposition(invK_i, R, U, D, V);
608
609 Real threshold = 0.0001f*smoothingLength;
610 D(0, 0) = D(0, 0) > threshold ? 1.0 / D(0, 0) : 1.0;
611 D(1, 1) = D(1, 1) > threshold ? 1.0 / D(1, 1) : 1.0;
612 D(2, 2) = D(2, 2) > threshold ? 1.0 / D(2, 2) : 1.0;
613
614 invK_i = V*D*U.transpose();
615
616 mat_i *= invK_i;
617
618 polarDecomposition(mat_i, R, U, D, V);
619
620 for (int ne = 0; ne < size_i; ne++)
621 {
622 Bond bond_ij = bonds_i[ne];
623 Coord rest_pos_j = X[bond_ij.idx];
624
625 Coord new_rest_pos_j = x_i + R*(rest_pos_j - x_i);
626 bond_ij.xi = new_rest_pos_j - x_i;
627 bonds_i[ne] = bond_ij;
628 }
629 }
630
631
632 template<typename TDataType>
633 void ElastoplasticityModule<TDataType>::rotateRestShape()
634 {
635 int num = this->inY()->size();
636 uint pDims = cudaGridSize(num, BLOCK_SIZE);
637
638 EM_RotateRestShape <Real, Coord, Matrix, Bond> << <pDims, BLOCK_SIZE >> > (
639 this->inX()->getData(),
640 this->inY()->getData(),
641 m_bYield,
642 this->inBonds()->getData(),
643 this->inHorizon()->getData());
644 cuSynchronize();
645 }
646
647 DEFINE_CLASS(ElastoplasticityModule);
648}
Generated by doxygen 1.13.2