PeriDyno 1.0.0
Loading...
Searching...
No Matches
EnergyDensityFunction.h
Go to the documentation of this file.
1
16#pragma once
17#include "Platform.h"
18#include "Matrix.h"
19
20namespace dyno
21{
35
36 template<typename Real, typename Matrix>
38 {
39 public:
40 //DYN_FUNC HyperelasticityModel() {};
41
42 DYN_FUNC virtual Real getEnergy(Real lambda1, Real lambda2, Real lambda3) = 0;
43 DYN_FUNC virtual Matrix getStressTensorPositive(Real lambda1, Real lambda2, Real lambda3) = 0;
44 DYN_FUNC virtual Matrix getStressTensorNegative(Real lambda1, Real lambda2, Real lambda3) = 0;
45
47 };
48
49 //dynamic initialization is not supported for a __constant__ variable
50 template<typename Real>
52 {
53 public:
54 //dynamic initialization is not supported for a __constant__ variable
55 DYN_FUNC LinearModel()
56 {
57// density = Real(1000);
58// s1 = Real(48000);
59// s0 = Real(12000);
60 }
61
62 DYN_FUNC LinearModel(Real _s0, Real _s1)
63 {
64 s0 = _s0;
65 s1 = _s1;
66 density = Real(1000);
67 }
68
69 DYN_FUNC Real getEnergy(Real lambda1, Real lambda2, Real lambda3)
70 {
71 return s0 * (lambda1 - 1)*(lambda1 - 1) + s0 * (lambda2 - 1)*(lambda2 - 1) + s0 * (lambda3 - 1)*(lambda3 - 1);
72 }
73
74 DYN_FUNC SquareMatrix<Real, 3> getStressTensorPositive(Real lambda1, Real lambda2, Real lambda3)
75 {
77 }
78
79 DYN_FUNC SquareMatrix<Real, 3> getStressTensorNegative(Real lambda1, Real lambda2, Real lambda3)
80 {
82 }
83
86
88 };
89
90
91 template<typename Real>
93 {
94 public:
95 //dynamic initialization is not supported for a __constant__ variable
96 DYN_FUNC StVKModel()
97 {
98// density = Real(1000);
99// s1 = Real(48000);
100// s0 = Real(12000);
101 }
102
103 DYN_FUNC StVKModel(Real _s0, Real _s1)
104 {
105 s0 = _s0;
106 s1 = _s1;
107 density = Real(1000);
108 }
109
110 DYN_FUNC Real getEnergy(Real lambda1, Real lambda2, Real lambda3)
111 {
112 Real I = lambda1 * lambda1 + lambda2 * lambda2 + lambda3 * lambda3;
113 Real sq1 = lambda1 * lambda1;
114 Real sq2 = lambda2 * lambda2;
115 Real sq3 = lambda3 * lambda3;
116 Real II = sq1 * sq1 + sq2 * sq2 + sq3 * sq3;
117 return 0.5*s0*(I - 3)*(I - 3) + 0.25*s1*(II - 2 * I + 3);
118 }
119
120 DYN_FUNC SquareMatrix<Real, 3> getStressTensorPositive(Real lambda1, Real lambda2, Real lambda3)
121 {
122 Real I = lambda1 * lambda1 + lambda2 * lambda2 + lambda3 * lambda3;
123
124 Real D1 = 2 * s0*I + s1 * lambda1*lambda1;
125 Real D2 = 2 * s0*I + s1 * lambda2*lambda2;
126 Real D3 = 2 * s0*I + s1 * lambda3*lambda3;
127
129 D(0, 0) = D1;
130 D(1, 1) = D2;
131 D(2, 2) = D3;
132 return D;
133 }
134
135 DYN_FUNC SquareMatrix<Real, 3> getStressTensorNegative(Real lambda1, Real lambda2, Real lambda3)
136 {
138 D(0, 0) *= lambda1;
139 D(1, 1) *= lambda2;
140 D(2, 2) *= lambda3;
141 return D;
142 }
143
146
148 };
149
150 template<typename Real>
152 {
153 public:
154 //dynamic initialization is not supported for a __constant__ variable
156 {
157// s0 = Real(12000);
158// s1 = Real(48000);
159// density = Real(1000);
160 }
161
162 DYN_FUNC NeoHookeanModel(Real _s0, Real _s1)
163 {
164 s0 = _s0;
165 s1 = _s1;
166 density = Real(1000);
167 }
168
169 DYN_FUNC Real getEnergy(Real lambda1, Real lambda2, Real lambda3)
170 {
171 Real I = lambda1 * lambda1 + lambda2 * lambda2 + lambda3 * lambda3;
172 Real sqrtIII = lambda1 * lambda2*lambda3;
173 return s0 * (I - 3 - 2 * log(sqrtIII)) + s1 * (sqrtIII - 1)*(sqrtIII - 1);
174 }
175
176 DYN_FUNC SquareMatrix<Real, 3> getStressTensorPositive(Real lambda1, Real lambda2, Real lambda3)
177 {
178 Real sq1 = lambda1 * lambda1;
179 Real sq2 = lambda2 * lambda2;
180 Real sq3 = lambda3 * lambda3;
181
182 Real D1 = 2 * s0 + 2 * s1*sq2*sq3;
183 Real D2 = 2 * s0 + 2 * s1*sq3*sq1;
184 Real D3 = 2 * s0 + 2 * s1*sq1*sq2;
185
187 D(0, 0) = D1;
188 D(1, 1) = D2;
189 D(2, 2) = D3;
190 return D;
191 }
192
193 DYN_FUNC SquareMatrix<Real, 3> getStressTensorNegative(Real lambda1, Real lambda2, Real lambda3)
194 {
195 Real sqrtIII = lambda1* lambda2* lambda3;
196
197 Real sq1 = lambda1 * lambda1;
198 Real sq2 = lambda2 * lambda2;
199 Real sq3 = lambda3 * lambda3;
200
202 D(0, 0) = 2 * s0 / sq1 + 2 * s1 * sqrtIII / sq1;
203 D(1, 1) = 2 * s0 / sq2 + 2 * s1 * sqrtIII / sq2;
204 D(2, 2) = 2 * s0 / sq3 + 2 * s1 * sqrtIII / sq3;
205
206 return D;
207 }
208
211
213 };
214
215 template<typename Real, int n>
217 {
218 public:
220 {
221 density = Real(1000);
222
223 for (int i = 0; i <= n; i++)
224 {
225 for (int j = 0; j <= n; j++)
226 {
227 C[i][j] = Real(12000);
228 }
229 }
230 }
231
232 DYN_FUNC Real getEnergy(Real lambda1, Real lambda2, Real lambda3)
233 {
234 Real sq1 = lambda1 * lambda1;
235 Real sq2 = lambda2 * lambda2;
236 Real sq3 = lambda3 * lambda3;
237 Real I1 = sq1 + sq2 + sq3;
238 Real I2 = sq1 * sq2 + sq2 * sq3 + sq3 * sq1;
239
240 Real totalE = Real(0);
241 for (int i = 0; i <= n; i++)
242 {
243 for (int j = 0; j <= n; j++)
244 {
245 totalE += C[i][j] * pow(I1 - 3, Real(i))*pow(I2 - 3, Real(j));
246 }
247 }
248
249 return totalE;
250 }
251
252 DYN_FUNC SquareMatrix<Real, 3> getStressTensorPositive(Real lambda1, Real lambda2, Real lambda3)
253 {
254 Real sq1 = lambda1 * lambda1;
255 Real sq2 = lambda2 * lambda2;
256 Real sq3 = lambda3 * lambda3;
257 Real I1 = sq1 + sq2 + sq3;
258 Real I2 = sq1 * sq2 + sq2 * sq3 + sq3 * sq1;
259
260 Real D1 = Real(0);
261 Real D2 = Real(0);
262 Real D3 = Real(0);
263 for (int i = 0; i <= n; i++)
264 {
265 for (int j = 0; j <= n; j++)
266 {
267 int i_minus_one = i - 1 < 0 ? 0 : i - 1;
268 int j_minus_one = j - 1 < 0 ? 0 : j - 1;
269
270 int i_minus_one_even = 2 * floor(i_minus_one / 2);
271 int j_minus_one_even = 2 * floor(j_minus_one / 2);
272 int i_even = 2 * floor(i / 2);
273 int j_even = 2 * floor(j / 2);
274
275 Real C1_positive;
276 if (i_minus_one_even == i_minus_one && j_even == j)
277 {
278 C1_positive = Real(1);
279 }
280 else if (i_minus_one_even != i_minus_one && j_even == j)
281 {
282 C1_positive = I1;
283 }
284 else if (i_minus_one_even == i_minus_one && j_even != j)
285 {
286 C1_positive = I2;
287 }
288 else
289 {
290 C1_positive = I1 * I2 + 9;
291 }
292
293 Real C2_positive;
294 if (i_even == i && j_minus_one_even == j_minus_one)
295 {
296 C2_positive = Real(1);
297 }
298 else if (i_even != i && j_minus_one_even == j_minus_one)
299 {
300 C2_positive = I1;
301 }
302 else if (i_even == i && j_minus_one_even != j_minus_one)
303 {
304 C2_positive = I2;
305 }
306 else
307 {
308 C2_positive = I1 * I2 + 9;
309 }
310
311 int exp_i = i < 0 ? 1 : i;
312 int exp_j = j < 0 ? 1 : j;
313 D1 += 2 * exp_i * C[i][j] * pow(I1 - 3, Real(i_minus_one_even))*pow(I2 - 3, Real(j_even)) * C1_positive + 4 * exp_j * C[i][j] * pow(I1 - 3, Real(i_even))*pow(I2 - 3, Real(j_minus_one_even))*C2_positive*sq1;
314 D2 += 2 * exp_i * C[i][j] * pow(I1 - 3, Real(i_minus_one_even))*pow(I2 - 3, Real(j_even)) * C1_positive + 4 * exp_j * C[i][j] * pow(I1 - 3, Real(i_even))*pow(I2 - 3, Real(j_minus_one_even))*C2_positive*sq2;
315 D3 += 2 * exp_i * C[i][j] * pow(I1 - 3, Real(i_minus_one_even))*pow(I2 - 3, Real(j_even)) * C1_positive + 4 * exp_j * C[i][j] * pow(I1 - 3, Real(i_even))*pow(I2 - 3, Real(j_minus_one_even))*C2_positive*sq3;
316
317 }
318 }
319
321 D(0, 0) = D1 * lambda1;
322 D(1, 1) = D2 * lambda2;
323 D(2, 2) = D3 * lambda3;
324
325 return D;
326 }
327
328 DYN_FUNC SquareMatrix<Real, 3> getStressTensorNegative(Real lambda1, Real lambda2, Real lambda3)
329 {
330 Real sq1 = lambda1 * lambda1;
331 Real sq2 = lambda2 * lambda2;
332 Real sq3 = lambda3 * lambda3;
333 Real I1 = sq1 + sq2 + sq3;
334 Real I2 = sq1 * sq2 + sq2 * sq3 + sq3 * sq1;
335
336 Real D1 = Real(0);
337 Real D2 = Real(0);
338 Real D3 = Real(0);
339 for (int i = 0; i <= n; i++)
340 {
341 for (int j = 0; j <= n; j++)
342 {
343 int i_minus_one = i - 1 < 0 ? 0 : i - 1;
344 int j_minus_one = j - 1 < 0 ? 0 : j - 1;
345
346 int i_minus_one_even = 2 * floor(i_minus_one / 2);
347 int j_minus_one_even = 2 * floor(j_minus_one / 2);
348 int i_even = 2 * floor(i / 2);
349 int j_even = 2 * floor(j / 2);
350
351 Real C1_positive;
352 if (i_minus_one_even == i_minus_one && j_even == j)
353 {
354 C1_positive = 0;
355 }
356 else if (i_minus_one_even != i_minus_one && j_even == j)
357 {
358 C1_positive = 3;
359 }
360 else if (i_minus_one_even == i_minus_one && j_even != j)
361 {
362 C1_positive = 3;
363 }
364 else
365 {
366 C1_positive = 3 * I1 + 3 * I2;
367 }
368
369 Real C2_positive;
370 if (i_even == i && j_minus_one_even == j_minus_one)
371 {
372 C2_positive = 0;
373 }
374 else if (i_even != i && j_minus_one_even == j_minus_one)
375 {
376 C2_positive = 3;
377 }
378 else if (i_even == i && j_minus_one_even != j_minus_one)
379 {
380 C2_positive = 3;
381 }
382 else
383 {
384 C2_positive = 3 * I1 + 3 * I2;
385 }
386
387 int exp_i = i < 0 ? 1 : i;
388 int exp_j = j < 0 ? 1 : j;
389 D1 += 2 * exp_i * C[i][j] * pow(I1 - 3, Real(i_minus_one_even))*pow(I2 - 3, Real(j_even)) * C1_positive + 4 * exp_j * C[i][j] * pow(I1 - 3, Real(i_even))*pow(I2 - 3, Real(j_minus_one_even))*C2_positive*sq1;
390 D2 += 2 * exp_i * C[i][j] * pow(I1 - 3, Real(i_minus_one_even))*pow(I2 - 3, Real(j_even)) * C1_positive + 4 * exp_j * C[i][j] * pow(I1 - 3, Real(i_even))*pow(I2 - 3, Real(j_minus_one_even))*C2_positive*sq2;
391 D3 += 2 * exp_i * C[i][j] * pow(I1 - 3, Real(i_minus_one_even))*pow(I2 - 3, Real(j_even)) * C1_positive + 4 * exp_j * C[i][j] * pow(I1 - 3, Real(i_even))*pow(I2 - 3, Real(j_minus_one_even))*C2_positive*sq3;
392
393 }
394 }
395
397 D(0, 0) = D1;
398 D(1, 1) = D2;
399 D(2, 2) = D3;
400
401 return D;
402 }
403
404 Real C[n + 1][n + 1];
405
407 };
408
409
410 template<typename Real>
412 {
413 public:
414 //dynamic initialization is not supported for a __constant__ variable
415 DYN_FUNC XuModel()
416 {
417// density = Real(1000);
418// s0 = Real(48000);
419 }
420
421 DYN_FUNC XuModel(Real _s0)
422 {
423 s0 = _s0;
424 density = Real(1000);
425 }
426
427 DYN_FUNC Real getEnergy(Real lambda1, Real lambda2, Real lambda3)
428 {
429 Real sq1 = lambda1 * lambda1;
430 Real sq2 = lambda2 * lambda2;
431 Real sq3 = lambda3 * lambda3;
432 Real E1 = ((sq1*sq1 - 1) / 4 + (1 / sq1 - 1) / 2) / 3;
433 Real E2 = ((sq2*sq2 - 1) / 4 + (1 / sq2 - 1) / 2) / 3;
434 Real E3 = ((sq3*sq3 - 1) / 4 + (1 / sq3 - 1) / 2) / 3;
435
436 return s0 * (E1 + E2 + E3);
437 }
438
439 DYN_FUNC SquareMatrix<Real, 3> getStressTensorPositive(Real lambda1, Real lambda2, Real lambda3)
440 {
442 D(0, 0) = s0 * (lambda1*lambda1) / 3;
443 D(1, 1) = s0 * (lambda2*lambda2) / 3;
444 D(2, 2) = s0 * (lambda3*lambda3) / 3;
445 return D;
446 }
447
448 DYN_FUNC SquareMatrix<Real, 3> getStressTensorNegative(Real lambda1, Real lambda2, Real lambda3)
449 {
451 D(0, 0) = lambda1 * s0 / (lambda1*lambda1*lambda1*lambda1) / 3;
452 D(1, 1) = lambda2 * s0 / (lambda2*lambda2*lambda2*lambda2) / 3;
453 D(2, 2) = lambda3 * s0 / (lambda3*lambda3*lambda3*lambda3) / 3;
454 return D;
455 }
456
458
460 };
461
462 template<typename Real>
464 {
465 public:
467 {
468 }
469
470 DYN_FUNC MooneyRivlinModel(Real _s0, Real _s1, Real _s2)
471 {
472 s0 = _s0;
473 s1 = _s1;
474 s2 = _s2;
475 density = Real(1000);
476 }
477
478 DYN_FUNC Real getEnergy(Real lambda1, Real lambda2, Real lambda3)
479 {
480 Real I = lambda1 * lambda1 + lambda2 * lambda2 + lambda3 * lambda3;
481 Real sq1 = lambda1 * lambda1;
482 Real sq2 = lambda2 * lambda2;
483 Real sq3 = lambda3 * lambda3;
484 Real II = sq1 * sq1 + sq2 * sq2 + sq3 * sq3;
485 Real J = lambda1 * lambda2 * lambda3;
486 Real sqrtIII = 1 / (lambda1 * lambda2 * lambda3);
487 Real cbrtIII = 1 / cbrt(sq1 * sq2 * sq3);
488 return s0 * (cbrtIII * I - 3) + s1 * (J - 1) * (J - 1) + s2 * (0.5 * cbrtIII * cbrtIII * (I * I - II) - 3);
489 }
490
491 DYN_FUNC SquareMatrix<Real, 3> getStressTensorPositive(Real lambda1, Real lambda2, Real lambda3)
492 {
493 Real sq1 = lambda1 * lambda1;
494 Real sq2 = lambda2 * lambda2;
495 Real sq3 = lambda3 * lambda3;
496
497 Real I = lambda1 * lambda1 + lambda2 * lambda2 + lambda3 * lambda3;
498 Real II = sq1 * sq1 + sq2 * sq2 + sq3 * sq3;
499 Real J = lambda1 * lambda2 * lambda3;
500 Real sqrtIII = 1 / (lambda1 * lambda2 * lambda3);
501 Real cbrtIII = 1 / cbrt(sq1 * sq2 * sq3);
502
503 Real D1 = 2 * s0 * cbrtIII + 2.0 * s1 * J * J + 2.0 / 3.0 * s2 / lambda1 * II * cbrtIII * cbrtIII + 2 * s2 * I * cbrtIII * cbrtIII;
504 Real D2 = 2 * s0 * cbrtIII + 2.0 * s1 * J * J + 2.0 / 3.0 * s2 / lambda2 * II * cbrtIII * cbrtIII + 2 * s2 * I * cbrtIII * cbrtIII;
505 Real D3 = 2 * s0 * cbrtIII + 2.0 * s1 * J * J + 2.0 / 3.0 * s2 / lambda3 * II * cbrtIII * cbrtIII + 2 * s2 * I * cbrtIII * cbrtIII;
506
508 D(0, 0) = D1;
509 D(1, 1) = D2;
510 D(2, 2) = D3;
511 return D;
512 }
513
514 DYN_FUNC SquareMatrix<Real, 3> getStressTensorNegative(Real lambda1, Real lambda2, Real lambda3)
515 {
516 Real sq1 = lambda1 * lambda1;
517 Real sq2 = lambda2 * lambda2;
518 Real sq3 = lambda3 * lambda3;
519
520 Real I = lambda1 * lambda1 + lambda2 * lambda2 + lambda3 * lambda3;
521 Real II = sq1 * sq1 + sq2 * sq2 + sq3 * sq3;
522 Real J = lambda1 * lambda2 * lambda3;
523 Real sqrtIII = 1 / (lambda1 * lambda2 * lambda3);
524 Real cbrtIII = 1 / cbrt(sq1 * sq2 * sq3);
525
527 D(0, 0) = 2.0 / 3.0 * s0 / lambda1 * cbrtIII * I + 2 * s1 / lambda1 * J + 2.0 / 3.0 * s2 * cbrtIII * cbrtIII / lambda1 * I * I + 2 * s2 * lambda1 * lambda1 * lambda1 * cbrtIII * cbrtIII;
528 D(1, 1) = 2.0 / 3.0 * s0 / lambda2 * cbrtIII * I + 2 * s1 / lambda2 * J + 2.0 / 3.0 * s2 * cbrtIII * cbrtIII / lambda2 * I * I + 2 * s2 * lambda2 * lambda2 * lambda2 * cbrtIII * cbrtIII;
529 D(2, 2) = 2.0 / 3.0 * s0 / lambda3 * cbrtIII * I + 2 * s1 / lambda3 * J + 2.0 / 3.0 * s2 * cbrtIII * cbrtIII / lambda3 * I * I + 2 * s2 * lambda3 * lambda3 * lambda3 * cbrtIII * cbrtIII;
530
531 return D;
532 }
533
538 };
539
540 template<typename Real>
542 {
543 public:
544 DYN_FUNC FungModel()
545 {
546 }
547
548 DYN_FUNC FungModel(Real _s0, Real _s1, Real _s2, Real _s3)
549 {
550 s0 = _s0;
551 s1 = _s1;
552 s2 = _s2;
553 s3 = _s3;
554 density = Real(1000);
555 }
556
557 DYN_FUNC Real getEnergy(Real lambda1, Real lambda2, Real lambda3)
558 {
559 Real I = lambda1 * lambda1 + lambda2 * lambda2 + lambda3 * lambda3;
560 Real sq1 = lambda1 * lambda1;
561 Real sq2 = lambda2 * lambda2;
562 Real sq3 = lambda3 * lambda3;
563 Real II = sq1 * sq1 + sq2 * sq2 + sq3 * sq3;
564 Real J = lambda1 * lambda2 * lambda3;
565 Real sqrtIII = 1 / (lambda1 * lambda2 * lambda3);
566 Real cbrtIII = 1 / cbrt(sq1 * sq2 * sq3);
567 return s0 * (cbrtIII * I - 3) + s1 * (J - 1) * (J - 1) + s2 * exp(s3 * (cbrtIII * I - 3) - 1);
568 }
569
570 DYN_FUNC SquareMatrix<Real, 3> getStressTensorPositive(Real lambda1, Real lambda2, Real lambda3)
571 {
572 Real sq1 = lambda1 * lambda1;
573 Real sq2 = lambda2 * lambda2;
574 Real sq3 = lambda3 * lambda3;
575
576 Real I = lambda1 * lambda1 + lambda2 * lambda2 + lambda3 * lambda3;
577 Real II = sq1 * sq1 + sq2 * sq2 + sq3 * sq3;
578 Real J = lambda1 * lambda2 * lambda3;
579 Real sqrtIII = 1 / (lambda1 * lambda2 * lambda3);
580 Real cbrtIII = 1 / cbrt(sq1 * sq2 * sq3);
581
582 Real D1 = 2 * s0 * cbrtIII + 2.0 * s1 * J * J + 2 * s2 * s3 * cbrtIII * exp(s3 * (cbrtIII * I - 3));
583 Real D2 = 2 * s0 * cbrtIII + 2.0 * s1 * J * J + 2 * s2 * s3 * cbrtIII * exp(s3 * (cbrtIII * I - 3));
584 Real D3 = 2 * s0 * cbrtIII + 2.0 * s1 * J * J + 2 * s2 * s3 * cbrtIII * exp(s3 * (cbrtIII * I - 3));
585
587 D(0, 0) = D1;
588 D(1, 1) = D2;
589 D(2, 2) = D3;
590 return D;
591 }
592
593 DYN_FUNC SquareMatrix<Real, 3> getStressTensorNegative(Real lambda1, Real lambda2, Real lambda3)
594 {
595 Real sq1 = lambda1 * lambda1;
596 Real sq2 = lambda2 * lambda2;
597 Real sq3 = lambda3 * lambda3;
598
599 Real I = lambda1 * lambda1 + lambda2 * lambda2 + lambda3 * lambda3;
600 Real II = sq1 * sq1 + sq2 * sq2 + sq3 * sq3;
601 Real J = lambda1 * lambda2 * lambda3;
602 Real sqrtIII = 1 / (lambda1 * lambda2 * lambda3);
603 Real cbrtIII = 1 / cbrt(sq1 * sq2 * sq3);
604
606 D(0, 0) = 2.0 / 3.0 * s0 / lambda1 * cbrtIII * I + 2 * s1 / lambda1 * J + 2.0 / 3.0 * s2 * s3 * cbrtIII / lambda1 * I * exp(s3 * (cbrtIII * I - 3));
607 D(1, 1) = 2.0 / 3.0 * s0 / lambda2 * cbrtIII * I + 2 * s1 / lambda2 * J + 2.0 / 3.0 * s2 * s3 * cbrtIII / lambda2 * I * exp(s3 * (cbrtIII * I - 3));
608 D(2, 2) = 2.0 / 3.0 * s0 / lambda3 * cbrtIII * I + 2 * s1 / lambda2 * J + 2.0 / 3.0 * s2 * s3 * cbrtIII / lambda3 * I * exp(s3 * (cbrtIII * I - 3));
609
610 return D;
611 }
612
618 };
619
620 template<typename Real>
622 {
623 public:
625 {
626 }
627
628 DYN_FUNC ArrudaBoyceModel(Real _s0, Real _s1, Real _s2)
629 {
630 s0 = _s0;
631 s1 = _s1;
632 s2 = _s2;
633 density = Real(1000);
634 }
635
636 DYN_FUNC Real getEnergy(Real lambda1, Real lambda2, Real lambda3)
637 {
638 Real I = lambda1 * lambda1 + lambda2 * lambda2 + lambda3 * lambda3;
639 Real sq1 = lambda1 * lambda1;
640 Real sq2 = lambda2 * lambda2;
641 Real sq3 = lambda3 * lambda3;
642 Real II = sq1 * sq1 + sq2 * sq2 + sq3 * sq3;
643 //Real J = lambda1 * lambda2 * lambda3;
644 //Real sqrtIII = 1 / (lambda1 * lambda2 * lambda3);
645 //Real cbrtIII = 1 / cbrt(sq1 * sq2 * sq3);
646 return 0.5 * s0 * (II - 3) + 0.25 * s0 * s1 * (II * II - 9) + s0 * s2 / 6.0 * (II * II * II - 27);
647 }
648
649 DYN_FUNC SquareMatrix<Real, 3> getStressTensorPositive(Real lambda1, Real lambda2, Real lambda3)
650 {
651 Real sq1 = lambda1 * lambda1;
652 Real sq2 = lambda2 * lambda2;
653 Real sq3 = lambda3 * lambda3;
654
655 Real II = sq1 * sq1 + sq2 * sq2 + sq3 * sq3;
656
657 Real D1 = 2 * s0 * lambda1 * lambda1 + 2.0 * lambda1 * lambda1 * s0 * s1 * II + 2.0 * s0 * s2 * lambda1 * lambda1 * II * II;
658 Real D2 = 2 * s0 * lambda2 * lambda2 + 2.0 * lambda2 * lambda2 * s0 * s1 * II + 2.0 * s0 * s2 * lambda2 * lambda2 * II * II;
659 Real D3 = 2 * s0 * lambda3 * lambda3 + 2.0 * lambda3 * lambda3 * s0 * s1 * II + 2.0 * s0 * s2 * lambda3 * lambda3 * II * II;
660
662 D(0, 0) = D1;
663 D(1, 1) = D2;
664 D(2, 2) = D3;
665 return D;
666 }
667
668 DYN_FUNC SquareMatrix<Real, 3> getStressTensorNegative(Real lambda1, Real lambda2, Real lambda3)
669 {
670
672 D(0, 0) = 0.;
673 D(1, 1) = 0.;
674 D(2, 2) = 0.;
675
676 return D;
677 }
678
683 };
684
685 template<typename Real>
687 {
688 public:
689 DYN_FUNC OgdenModel()
690 {
691 }
692
693 DYN_FUNC OgdenModel(Real _s0, Real _s1)
694 {
695 s0 = _s0;
696 s1 = _s1;
697 density = Real(1000);
698 }
699
700 DYN_FUNC Real getEnergy(Real lambda1, Real lambda2, Real lambda3)
701 {
702 Real I = lambda1 * lambda1 + lambda2 * lambda2 + lambda3 * lambda3;
703 Real sq1 = lambda1 * lambda1;
704 Real sq2 = lambda2 * lambda2;
705 Real sq3 = lambda3 * lambda3;
706 Real II = sq1 * sq1 + sq2 * sq2 + sq3 * sq3;
707 Real J = lambda1 * lambda2 * lambda3;
708 Real sqrtIII = 1 / (lambda1 * lambda2 * lambda3);
709 Real cbrtIII = 1 / cbrt(sq1 * sq2 * sq3);
710 return s0 * (II - 3) + 2.0 * s0 * log(J) + s1 * (J - 1) * (J - 1);
711 }
712
713 DYN_FUNC SquareMatrix<Real, 3> getStressTensorPositive(Real lambda1, Real lambda2, Real lambda3)
714 {
715 Real sq1 = lambda1 * lambda1;
716 Real sq2 = lambda2 * lambda2;
717 Real sq3 = lambda3 * lambda3;
718
719 Real II = sq1 * sq1 + sq2 * sq2 + sq3 * sq3;
720 Real J = lambda1 * lambda2 * lambda3;
721 Real sqrtIII = 1 / (lambda1 * lambda2 * lambda3);
722 Real cbrtIII = 1 / cbrt(sq1 * sq2 * sq3);
723
724 Real D1 = 4.0 * s0 * lambda1 * lambda1 * lambda1 + 2.0 * s1 * J * J / lambda1;
725 Real D2 = 4.0 * s0 * lambda2 * lambda2 * lambda2 + 2.0 * s1 * J * J / lambda2;
726 Real D3 = 4.0 * s0 * lambda3 * lambda3 * lambda3 + 2.0 * s1 * J * J / lambda3;
727
729 D(0, 0) = D1 / lambda1;
730 D(1, 1) = D2 / lambda2;
731 D(2, 2) = D3 / lambda3;
732 return D;
733 }
734
735 DYN_FUNC SquareMatrix<Real, 3> getStressTensorNegative(Real lambda1, Real lambda2, Real lambda3)
736 {
737 Real sq1 = lambda1 * lambda1;
738 Real sq2 = lambda2 * lambda2;
739 Real sq3 = lambda3 * lambda3;
740
741 Real II = sq1 * sq1 + sq2 * sq2 + sq3 * sq3;
742 Real J = lambda1 * lambda2 * lambda3;
743 Real sqrtIII = 1 / (lambda1 * lambda2 * lambda3);
744 Real cbrtIII = 1 / cbrt(sq1 * sq2 * sq3);
746 D(0, 0) = 4.0 * s0 / lambda1 + 2.0 * s1 * J / lambda1;
747 D(1, 1) = 4.0 * s0 / lambda2 + 2.0 * s1 * J / lambda2;
748 D(2, 2) = 4.0 * s0 / lambda3 + 2.0 * s1 * J / lambda3;
749
750 return D;
751 }
752
756 };
757
758 template<typename Real>
760 {
761 public:
762 DYN_FUNC YeohModel()
763 {
764 }
765
766 DYN_FUNC YeohModel(Real _s0, Real _s1, Real _s2)
767 {
768 s0 = _s0;
769 s1 = _s1;
770 s2 = _s2;
771 density = Real(1000);
772 }
773
774 DYN_FUNC Real getEnergy(Real lambda1, Real lambda2, Real lambda3)
775 {
776 Real I = lambda1 * lambda1 + lambda2 * lambda2 + lambda3 * lambda3;
777 return s0 * (I - 3) + s1 * (I - 3) * (I - 3)/* + s2 * (I - 3) * (I - 3) * (I - 3)*/;
778 }
779
780 DYN_FUNC SquareMatrix<Real, 3> getStressTensorPositive(Real lambda1, Real lambda2, Real lambda3)
781 {
782
783
784 Real I = lambda1 * lambda1 + lambda2 * lambda2 + lambda3 * lambda3;
785
786
787 Real D1 = 2.0 * s0 + 4 * s1 * I + 6 * lambda1 * s2 * I * I + 54 * s2;
788 Real D2 = 2.0 * s0 + 4 * s1 * I + 6 * lambda2 * s2 * I * I + 54 * s2;
789 Real D3 = 2.0 * s0 + 4 * s1 * I + 6 * lambda3 * s2 * I * I + 54 * s2;
790
792 D(0, 0) = D1;
793 D(1, 1) = D2;
794 D(2, 2) = D3;
795 return D;
796 }
797
798 DYN_FUNC SquareMatrix<Real, 3> getStressTensorNegative(Real lambda1, Real lambda2, Real lambda3)
799 {
800
801
802 Real I = lambda1 * lambda1 + lambda2 * lambda2 + lambda3 * lambda3;
803
805 D(0, 0) = 12.0 * lambda1 * s1 + 36.0 * lambda1 * s2 * I;
806 D(1, 1) = 12.0 * lambda2 * s1 + 36.0 * lambda2 * s2 * I;
807 D(2, 2) = 12.0 * lambda3 * s1 + 36.0 * lambda3 * s2 * I;
808
809 return D;
810 }
811
816 };
817
818 //dynamic initialization is not supported for a __constant__ variable
819 template<typename Real>
821 {
822 public:
823 //dynamic initialization is not supported for a __constant__ variable
824 DYN_FUNC FiberModel()
825 {
826 // density = Real(1000);
827 /* s[0] = 1;
828 s[1] = 1;
829 s[2] = 1;
830 */
831 }
832
833 DYN_FUNC FiberModel(Real _s0, Real _s1, Real _s2)
834 {
835 s[0] = _s0;
836 s[1] = _s1;
837 s[2] = _s2;
838 density = Real(1);
839 }
840
841 DYN_FUNC inline SquareMatrix<Real, 3> getF_TF(Real lambda1, Real lambda2, Real lambda3, SquareMatrix<Real, 3> V) {
843 res(0, 0) = lambda1 * lambda1;
844 res(1, 1) = lambda2 * lambda2;
845 res(2, 2) = lambda3 * lambda3;
846 res = V * res * V.transpose();
847 return res;
848
849 }
850
851 DYN_FUNC inline SquareMatrix<Real, 3> get_diag(Real l1, Real l2, Real l3) {
853 res(0, 0) = l1;
854 res(1, 1) = l2;
855 res(2, 2) = l3;
856 return res;
857 }
858 DYN_FUNC inline SquareMatrix<Real, 3> get_diag_bool(bool l1, bool l2, bool l3) {
860 res(0, 0) = l1 ? (Real) 1 : (Real) 0;
861 res(1, 1) = l2 ? (Real) 1 : (Real) 0;
862 res(2, 2) = l3 ? (Real) 1 : (Real) 0;
863 return res;
864
865 }
866 DYN_FUNC inline Real getJ(Real l1, Real l2, Real l3) {
867 return l1 * l2 * l3;
868 }
869 DYN_FUNC inline SquareMatrix<Real, 3> get_entry(int i, int j) {
871 if (i >= 3 || j >= 3 || i < 0 || j < 0)
872 return res;
873 res(i, j) = (Real) 1.0;
874 res(j, i) = (Real) 1.0;
875 return res;
876 }
877
878 DYN_FUNC Real getEnergy(Real lambda1, Real lambda2, Real lambda3, SquareMatrix<Real,3> V)
879 {
880 auto FTF = getF_TF(lambda1, lambda2, lambda3, V);
881 Real energy = 0.0;
882 for (int i = 0; i < 3; i++) {
883 energy += pow(s[i] * (FTF(i, i) - 1), 2);
884 energy += s[i] * s[(i + 1) % 3] * pow(FTF(i, (i + 1) % 3), 2) / (FTF(i, i) * FTF((i + 1) % 3, (i + 1) % 3));
885 }
886 energy += s[0] * s[1] * s[2] * pow(getJ(lambda1, lambda2, lambda3) -1,2);
887 return density * energy;
888 }
889
891 {
893 auto FTF = getF_TF(lambda1, lambda2, lambda3, V);
894 auto J = getJ(lambda1, lambda2, lambda3);
895 for (int i = 0; i < 3; i++) {
896 if (FTF(i, i) >= (1.0+0.05)) {
897 res += 400 * (FTF(i, i) - 1) * pow(s[i], 2) * get_diag_bool(i==0, i==1, i==2);
898 }
899 /*
900 if (FTF(i, (i + 1) % 3) >= 0.05)
901 res += 2 * s[i] * s[(i + 1) % 3] * FTF(i, (i + 1) % 3) / (FTF(i, i) * FTF((i + 1) % 3, (i + 1) % 3)) * get_entry(i, (i + 1) % 3); */
902 }
903 res = get_diag(lambda1, lambda2, lambda3) * V.transpose() * res * V * get_diag(1 / lambda1, 1 / lambda2, 1 / lambda3);
904 if (J >= (1.0 + 0.05)) {
905 res += 20 * J * (J-1) * s[0] * s[1] * s[2] * get_diag(1/pow(lambda1,2), 1/pow(lambda2,2), 1/pow(lambda3,2));
906 }
907
908 return density * res;
909 }
910
912 {
914 auto FTF = getF_TF(lambda1, lambda2, lambda3, V);
915 auto J = getJ(lambda1, lambda2, lambda3);
916
917 for (int i = 0; i < 3; i++) {
918 if (FTF(i, i) < (1.0 -0.05) )
919 {
920 res += 400 * (1 - FTF(i, i)) * pow(s[i], 2) * get_diag_bool(i == 0, i == 1, i == 2);
921 }
922
923 /*if (FTF(i, (i + 1) % 3) <= -0.05)*/
924 res += 200 * s[i] * s[(i + 1) % 3] * FTF(i, (i + 1) % 3) / (FTF(i, i) * FTF((i + 1) % 3, (i + 1) % 3)) * get_entry(i, (i + 1) % 3);
925 res += 200 * s[i] * s[(i + 1) % 3] * pow(FTF(i, (i + 1) % 3),2) / (FTF(i, i) * FTF((i + 1) % 3, (i + 1) % 3)) * (
926 get_diag_bool(i==0, i==1, i==2) / FTF(i, i) +
927 get_diag_bool((i+1)%3 == 0, (i+1)%3 == 1, (i+1)%3 == 2) / FTF((i + 1) % 3, (i + 1) % 3) );
928 }
929
930 res = get_diag(lambda1, lambda2, lambda3) * V.transpose() * res * V;
931
932 if (J < (1.0 - 0.05)) {
933 res += 20 * J * (1 - J) * s[0] * s[1] * s[2] * get_diag(1/lambda1, 1/lambda2, 1/lambda3);
934 }
935
936 return density * res;
937 }
938 DYN_FUNC void getInfo(Real lambda1, Real lambda2, Real lambda3, SquareMatrix<Real, 3> V)
939 {
940 /*
941 auto FTF = getF_TF(lambda1, lambda2, lambda3, V);
942 printf("=========================FTF=========================\nrow1:%f,\t%f,\t%f\nrow2:%f,\t%f,\t%f\nrow3:%f,\t%f,\t%f\n", FTF(0, 0), FTF(0, 1), FTF(0, 2),FTF(1, 0), FTF(1, 1), FTF(1, 2), FTF(2, 0), FTF(2, 1), FTF(2, 2));
943 */
944 }
945 Real s[3];
946
948 };
949
950 template<typename Real>
964}
double Real
Definition Typedef.inl:23
DYN_FUNC SquareMatrix< Real, 3 > getStressTensorNegative(Real lambda1, Real lambda2, Real lambda3)
DYN_FUNC Real getEnergy(Real lambda1, Real lambda2, Real lambda3)
DYN_FUNC SquareMatrix< Real, 3 > getStressTensorPositive(Real lambda1, Real lambda2, Real lambda3)
DYN_FUNC ArrudaBoyceModel(Real _s0, Real _s1, Real _s2)
DYN_FUNC SquareMatrix< Real, 3 > get_diag(Real l1, Real l2, Real l3)
DYN_FUNC Real getEnergy(Real lambda1, Real lambda2, Real lambda3, SquareMatrix< Real, 3 > V)
DYN_FUNC Real getJ(Real l1, Real l2, Real l3)
DYN_FUNC SquareMatrix< Real, 3 > get_diag_bool(bool l1, bool l2, bool l3)
DYN_FUNC SquareMatrix< Real, 3 > get_entry(int i, int j)
DYN_FUNC SquareMatrix< Real, 3 > getF_TF(Real lambda1, Real lambda2, Real lambda3, SquareMatrix< Real, 3 > V)
DYN_FUNC SquareMatrix< Real, 3 > getStressTensorNegative(Real lambda1, Real lambda2, Real lambda3, SquareMatrix< Real, 3 > V)
DYN_FUNC SquareMatrix< Real, 3 > getStressTensorPositive(Real lambda1, Real lambda2, Real lambda3, SquareMatrix< Real, 3 > V)
DYN_FUNC void getInfo(Real lambda1, Real lambda2, Real lambda3, SquareMatrix< Real, 3 > V)
DYN_FUNC FiberModel(Real _s0, Real _s1, Real _s2)
DYN_FUNC SquareMatrix< Real, 3 > getStressTensorNegative(Real lambda1, Real lambda2, Real lambda3)
DYN_FUNC SquareMatrix< Real, 3 > getStressTensorPositive(Real lambda1, Real lambda2, Real lambda3)
DYN_FUNC FungModel(Real _s0, Real _s1, Real _s2, Real _s3)
DYN_FUNC Real getEnergy(Real lambda1, Real lambda2, Real lambda3)
virtual DYN_FUNC Matrix getStressTensorPositive(Real lambda1, Real lambda2, Real lambda3)=0
virtual DYN_FUNC Real getEnergy(Real lambda1, Real lambda2, Real lambda3)=0
virtual DYN_FUNC Matrix getStressTensorNegative(Real lambda1, Real lambda2, Real lambda3)=0
DYN_FUNC SquareMatrix< Real, 3 > getStressTensorPositive(Real lambda1, Real lambda2, Real lambda3)
DYN_FUNC LinearModel(Real _s0, Real _s1)
DYN_FUNC SquareMatrix< Real, 3 > getStressTensorNegative(Real lambda1, Real lambda2, Real lambda3)
DYN_FUNC Real getEnergy(Real lambda1, Real lambda2, Real lambda3)
DYN_FUNC Real getEnergy(Real lambda1, Real lambda2, Real lambda3)
DYN_FUNC SquareMatrix< Real, 3 > getStressTensorPositive(Real lambda1, Real lambda2, Real lambda3)
DYN_FUNC SquareMatrix< Real, 3 > getStressTensorNegative(Real lambda1, Real lambda2, Real lambda3)
DYN_FUNC MooneyRivlinModel(Real _s0, Real _s1, Real _s2)
DYN_FUNC NeoHookeanModel(Real _s0, Real _s1)
DYN_FUNC Real getEnergy(Real lambda1, Real lambda2, Real lambda3)
DYN_FUNC SquareMatrix< Real, 3 > getStressTensorNegative(Real lambda1, Real lambda2, Real lambda3)
DYN_FUNC SquareMatrix< Real, 3 > getStressTensorPositive(Real lambda1, Real lambda2, Real lambda3)
DYN_FUNC SquareMatrix< Real, 3 > getStressTensorNegative(Real lambda1, Real lambda2, Real lambda3)
DYN_FUNC OgdenModel(Real _s0, Real _s1)
DYN_FUNC Real getEnergy(Real lambda1, Real lambda2, Real lambda3)
DYN_FUNC SquareMatrix< Real, 3 > getStressTensorPositive(Real lambda1, Real lambda2, Real lambda3)
DYN_FUNC Real getEnergy(Real lambda1, Real lambda2, Real lambda3)
DYN_FUNC SquareMatrix< Real, 3 > getStressTensorPositive(Real lambda1, Real lambda2, Real lambda3)
DYN_FUNC SquareMatrix< Real, 3 > getStressTensorNegative(Real lambda1, Real lambda2, Real lambda3)
DYN_FUNC StVKModel(Real _s0, Real _s1)
DYN_FUNC Real getEnergy(Real lambda1, Real lambda2, Real lambda3)
DYN_FUNC SquareMatrix< Real, 3 > getStressTensorNegative(Real lambda1, Real lambda2, Real lambda3)
DYN_FUNC SquareMatrix< Real, 3 > getStressTensorPositive(Real lambda1, Real lambda2, Real lambda3)
DYN_FUNC XuModel(Real _s0)
DYN_FUNC SquareMatrix< Real, 3 > getStressTensorNegative(Real lambda1, Real lambda2, Real lambda3)
DYN_FUNC Real getEnergy(Real lambda1, Real lambda2, Real lambda3)
DYN_FUNC SquareMatrix< Real, 3 > getStressTensorPositive(Real lambda1, Real lambda2, Real lambda3)
DYN_FUNC SquareMatrix< Real, 3 > getStressTensorPositive(Real lambda1, Real lambda2, Real lambda3)
DYN_FUNC Real getEnergy(Real lambda1, Real lambda2, Real lambda3)
DYN_FUNC YeohModel(Real _s0, Real _s1, Real _s2)
DYN_FUNC SquareMatrix< Real, 3 > getStressTensorNegative(Real lambda1, Real lambda2, Real lambda3)
#define V(a, b, c)
This is an implementation of AdditiveCCD based on peridyno.
Definition Array.h:25
DYN_FUNC Complex< Real > log(const Complex< Real > &)
Definition Complex.inl:305
DYN_FUNC Complex< Real > pow(const Complex< Real > &, const Real &)
Definition Complex.inl:370
DYN_FUNC Complex< Real > exp(const Complex< Real > &)
Definition Complex.inl:338
FungModel< Real > fungModel
FiberModel< Real > fiberModel
OgdenModel< Real > ogdenModel
LinearModel< Real > linearModel
YeohModel< Real > yeohModel
StVKModel< Real > stvkModel
MooneyRivlinModel< Real > mrModel
ArrudaBoyceModel< Real > abModel
NeoHookeanModel< Real > neohookeanModel