00001 #include "anisodam.h"
00002 #include "bmatrix.h"
00003 #include "matrix.h"
00004 #include "vector.h"
00005 #include "elastisomat.h"
00006 #include "global.h"
00007 #include "intpoints.h"
00008 #include "tensor.h"
00009 #include "vecttens.h"
00010 #include <math.h>
00011 #include <stdlib.h>
00012
00013 #define nijac 200
00014 #define limit 1.0e-15
00015 #ifndef M_PI
00016 #define M_PI 3.14159265358979323846
00017 #endif
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 anisodam::anisodam (void)
00030 {
00031 cde = corr_off;
00032 ac = 0.0; bc = 0.0; y0c = 0.0; at = 0.0; bt = 0.0; y0t = 0.0;
00033 a = 0.0; b = 0.0; y0 = 0.0; gfc = 0.0; gft = 0.0; gf = 0.0;
00034 ft = uft = 0.0;
00035 fc = ufc = 0.0;
00036 }
00037
00038
00039
00040
00041
00042
00043
00044
00045 anisodam::~anisodam (void)
00046 {
00047
00048 }
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063 void anisodam::read (XFILE *in)
00064 {
00065 xfscanf (in, "%k%m", "dam_evol_f", &dam_evolfunc_kwdset, (int *)&damevf);
00066 switch(damevf)
00067 {
00068 case brittle:
00069 xfscanf (in,"%k%lf %k%lf", "f_t", &ft, "u_ft", &uft);
00070 xfscanf (in,"%k%lf %k%lf", "f_c", &fc, "u_fc", &ufc);
00071 xfscanf (in,"%k%m", "corr_dissip", &corr_disip_en_kwdset, (int *)(&cde));
00072 if (cde != corr_off)
00073 sra.read (in);
00074 break;
00075 case quasi_brittle:
00076 xfscanf (in, "%k%m", "corr_dissip", &corr_disip_en_kwdset, (int *)(&cde));
00077 switch (cde)
00078 {
00079 case corr_off:
00080 xfscanf (in, "%k%lf %k%lf %k%lf", "A_c", &ac, "B_c", &bc, "Y_0c", &y0c);
00081 xfscanf (in, "%k%lf %k%lf %k%lf", "A_t", &at, "B_t", &bt, "Y_0t", &y0t);
00082 xfscanf (in, "%k%lf %k%lf %k%lf", "a", &a, "b", &b, "y0", &y0);
00083 break;
00084 case corr_on:
00085 xfscanf (in, "%k%lf %k%lf %k%lf", "G_fc", &gfc, "B_c", &bc, "Y_0c", &y0c);
00086 xfscanf (in, "%k%lf %k%lf %k%lf", "G_ft", &gft, "B_t", &bt, "Y_0t", &y0t);
00087 xfscanf (in, "%k%lf %k%lf %k%lf", "G_f", &gf, "b", &b, "y0", &y0);
00088 break;
00089 default:
00090 print_err(" unknown type of correction of disipated energy is required", __FILE__, __LINE__,__func__);
00091 }
00092 break;
00093 default:
00094 print_err(" unknown type of damage evolution function is required", __FILE__, __LINE__,__func__);
00095 }
00096 }
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110 double anisodam::damdrvforce_vol(double nu, vector &peps)
00111 {
00112 long i;
00113 double ev = 0.0, ret = 0.0;
00114 double k0_prime;
00115 double g0_prime;
00116
00117 for (i = 0; i < 3; i++)
00118 ev += peps(i);
00119 ev /= 3.0;
00120
00121
00122 if (ev <= 0.0)
00123 return 0.0;
00124
00125 k0_prime = 1.0/(1.0-2.0*nu);
00126 g0_prime = 1.0/(2.0*(1.0+nu));
00127 ret = 3.0/2.0*(k0_prime-2.0*g0_prime)*ev*ev;
00128 return ret;
00129 }
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202 void anisodam::damdrvforce_dev(double nu, vector &peps, vector &pyt, vector &pyc)
00203 {
00204 long i;
00205 double g0_prime;
00206
00207
00208 g0_prime = 1.0/(2.0*(1.0+nu));
00209 fillv(0.0, pyt);
00210 fillv(0.0, pyc);
00211 for (i = 0; i < 3; i++)
00212 {
00213 if (peps(i) > 0.0)
00214 pyt[i] = g0_prime*peps[i]*peps[i];
00215 if (peps(i) < 0.0)
00216 pyc[i] = g0_prime*peps[i]*peps[i];
00217 }
00218 }
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235 double anisodam::loadfuncvol(long ipp, double nu, vector &peps, double d, double aa)
00236 {
00237 double y, f, kappa;
00238 double tmp1, tmp2;
00239 double e, aft, auft;
00240
00241 y = damdrvforce_vol(nu, peps);
00242 switch(damevf)
00243 {
00244 case brittle:
00245 kappa = sqrt(y);
00246 if (kappa <= y0)
00247 f = -1.0;
00248 else
00249 {
00250 e = Mm->give_actual_ym(ipp);
00251
00252 aft = Mm->give_actual_ft(ipp);
00253
00254 auft = uft/(aft/ft);
00255 tmp1 = -(kappa-aft/e)/(auft-aft/e);
00256 f = 1.0 - aft/(e*kappa)*exp(tmp1) - d;
00257 }
00258 break;
00259 case quasi_brittle:
00260 if (y <= y0)
00261 f = -1.0;
00262 else
00263 {
00264
00265 tmp1 = y - y0;
00266 tmp2 = log(tmp1);
00267 tmp1 = exp(b*tmp2);
00268
00269 f = (1.0 - d)*(1.0+aa*tmp1) - 1.0;
00270 }
00271 break;
00272 default:
00273 f = 0.0;
00274 print_err(" unknown type of damage evolution function is required", __FILE__, __LINE__,__func__);
00275 }
00276 return f;
00277 }
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300 void anisodam::loadfuncdev(long ipp, double nu, vector &peps, vector &damt, vector &damc,
00301 double aat, double aac, vector &lft, vector &lfc)
00302 {
00303 long i;
00304 double e, af, auf, kappa, tmp1, tmp2;
00305 vector pyt(3), pyc(3);
00306
00307 damdrvforce_dev(nu, peps, pyt, pyc);
00308 switch(damevf)
00309 {
00310 case brittle:
00311 for (i=0; i<3; i++)
00312 {
00313
00314 kappa = sqrt(pyt(i));
00315 if (kappa <= y0t)
00316 lft(i) = -1.0;
00317 else
00318 {
00319 e = Mm->give_actual_ym(ipp);
00320
00321 af = Mm->give_actual_ft(ipp);
00322
00323 auf = uft/(af/ft);
00324 tmp1 = -(kappa-af/e)/(auf-af/e);
00325 lft(i) = 1.0 - af/(e*kappa)*exp(tmp1) - damt(i);
00326 }
00327
00328 kappa = sqrt(pyc(i));
00329 if (kappa <= y0c)
00330 lfc(i) = -1.0;
00331 else
00332 {
00333 e = Mm->give_actual_ym(ipp);
00334
00335 af = Mm->give_actual_fc(ipp);
00336
00337 auf = ufc/(af/fc);
00338 tmp1 = -(kappa-af/e)/(auf-af/e);
00339 lfc(i) = 1.0 - af/(e*kappa)*exp(tmp1) - damc(i);
00340 }
00341 }
00342 break;
00343 case quasi_brittle:
00344 for (i=0; i<3; i++)
00345 {
00346
00347 if (pyt(i) <= y0t)
00348 lft(i) = -1.0;
00349 else
00350 {
00351
00352 tmp1 = pyt(i) - y0t;
00353 tmp2 = log(tmp1);
00354 tmp1 = exp(bt*tmp2);
00355
00356 lft(i) = (1.0 - damt(i))*(1.0+aat*tmp1) - 1.0;
00357 }
00358
00359 if (pyc(i) <= y0c)
00360 lfc(i) = -1.0;
00361 else
00362 {
00363
00364 tmp1 = pyc(i) - y0c;
00365 tmp2 = log(tmp1);
00366 tmp1 = exp(bc*tmp2);
00367
00368 lfc(i) = (1.0 - damc(i))*(1.0+aac*tmp1) - 1.0;
00369 }
00370 }
00371 break;
00372 default:
00373 print_err(" unknown type of damage evolution function is required", __FILE__, __LINE__,__func__);
00374 }
00375 return;
00376 }
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395 double anisodam::daminc_vol(long ipp, double y, double dy, double aa, double lf)
00396 {
00397 double e = Mm->give_actual_ym(ipp);
00398 double aft, auft, tmp, kappa, ddam = 0.0;
00399
00400 switch(damevf)
00401 {
00402 case brittle:
00403 if ((lf > 0.0) && (dy > 0.0))
00404 {
00405 kappa = sqrt(y);
00406
00407 aft = Mm->give_actual_ft(ipp);
00408
00409 auft = uft/(aft/ft);
00410 tmp = -(kappa-aft/e)/(auft-aft/e);
00411 ddam = (aft/e)*(kappa-1.0)*exp(tmp)/(kappa*kappa);
00412 }
00413 break;
00414 case quasi_brittle:
00415 if ((lf > 0.0) && (dy > 0.0))
00416 {
00417 tmp = log(y-y0);
00418 ddam = aa*b*exp(tmp*(b-1.0));
00419 tmp = 1.0+a*exp(tmp*b);
00420 ddam /= e*tmp*tmp;
00421 ddam *= dy;
00422 }
00423 break;
00424 default:
00425 print_err(" unknown type of damage evolution function is required", __FILE__, __LINE__,__func__);
00426 }
00427 return ddam;
00428 }
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444
00445
00446
00447
00448
00449
00450
00451
00452
00453 void anisodam::pdaminc_dev(long ipp, vector &pyc, vector &pyt, vector &dyt, vector &dyc,
00454 double aat, double aac, vector &lft, vector &lfc, vector &dpdamt, vector &dpdamc)
00455 {
00456 long i;
00457 double e = Mm->give_actual_ym(ipp);
00458 double kappa, tmp, ltmp;
00459 double af, auf;
00460
00461 fillv(0.0, dpdamt);
00462 fillv(0.0, dpdamc);
00463
00464 switch(damevf)
00465 {
00466 case brittle:
00467 for(i=0; i<3; i++)
00468 {
00469
00470 kappa = sqrt(pyt(i));
00471 if ((dyt[i] > 0.0) && (lft[i] > 0.0))
00472 {
00473
00474 af = Mm->give_actual_ft(ipp);
00475
00476 auf = uft/(af/ft);
00477 tmp = -(kappa-af/e)/(auf-af/e);
00478 dpdamt(i) = (af/e)*(kappa-1.0)*exp(tmp)/(kappa*kappa);
00479 }
00480
00481 kappa = sqrt(pyc(i));
00482 if ((dyc[i] > 0.0) && (lfc[i] > 0.0))
00483 {
00484
00485 af = Mm->give_actual_fc(ipp);
00486
00487 auf = ufc/(af/fc);
00488 tmp = -(kappa-af/e)/(auf-af/e);
00489 dpdamc(i) = (af/e)*(kappa-1.0)*exp(tmp)/(kappa*kappa);
00490 }
00491 }
00492 break;
00493 case quasi_brittle:
00494 for(i=0; i<3; i++)
00495 {
00496
00497 if ((dyt[i] > 0.0) && (lft[i] > 0.0))
00498 {
00499 ltmp = log(pyt(i) - y0t);
00500 dpdamt(i) = aat*bt*exp(ltmp*(bt-1.0));
00501 tmp = 1.0+aat*exp(ltmp*bt);
00502 dpdamt(i) /= e*(1.0+tmp*tmp);
00503 dpdamt(i) *= dyt(i);
00504 }
00505
00506 if ((dyc[i] > 0.0) && (lfc[i] > 0.0))
00507 {
00508 ltmp = log(pyc(i) - y0c);
00509 dpdamc(i) = aac*bc*exp(ltmp*(bc-1.0));
00510 tmp = 1.0+aac*exp(ltmp*bc);
00511 dpdamc(i) /= e*(1.0+tmp*tmp);
00512 dpdamc(i) *= dyc(i);
00513 }
00514 }
00515 break;
00516 default:
00517 print_err(" unknown type of damage evolution function is required", __FILE__, __LINE__,__func__);
00518 }
00519 }
00520
00521
00522
00523
00524
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535
00536
00537
00538
00539 double anisodam::brittle_damage(long ipp, double y, double e, double f, double uf, double omegao)
00540 {
00541 double err,omega = 0.0, dtmp, tmp, h, indam, kappa;
00542 long i,ni, eid;
00543
00544 ni=sra.give_ni ();
00545 err=sra.give_err ();
00546
00547 indam = f/e;
00548
00549 kappa = sqrt(y);
00550 if (kappa < indam)
00551
00552 return 0.0;
00553 if ((Mm->ip[ipp].hmt & 2) > 0)
00554 {
00555 omega = 1.0-(f/(e*kappa))*exp(-(kappa-f/e)/(uf-f/e));
00556 return omega;
00557 }
00558
00559 eid = Mm->elip[ipp];
00560 switch (Mt->give_dimension(eid))
00561 {
00562 case 1:
00563 h = Mt->give_length(eid);
00564 break;
00565 case 2:
00566 h = sqrt(Mt->give_area(eid));
00567 break;
00568 case 3:
00569 h = pow(Mt->give_volume(eid), 1.0/3.0);
00570 break;
00571 default:
00572 print_err("unknown dimension of element", __FILE__, __LINE__, __func__);
00573 }
00574 if (cde == corr_off)
00575
00576 {
00577 omega = 1.0-(f/(e*kappa))*exp(-(kappa-f/e)/(uf-f/e));
00578 return omega;
00579 }
00580
00581 omega = 1.0;
00582 for (i = 0; i < ni; i++)
00583 {
00584 dtmp = -e*kappa+f/uf*h*kappa*exp(-h*omega*kappa/uf);
00585 tmp = (1.0-omega)*e*kappa-f*exp(-h*omega*kappa/uf);
00586 if (fabs(tmp) < ft*err)
00587 {
00588 omega += - tmp/dtmp;
00589 break;
00590 }
00591 omega += - tmp/dtmp;
00592 }
00593 if ((omega > 1.0) || (omega < 0.0) || (fabs(tmp) > ft*err))
00594
00595
00596 {
00597 omega = omegao;
00598 for (i = 0; i < ni; i++)
00599 {
00600 dtmp = -e*kappa+ft/uf*h*kappa*exp(-h*omega*kappa/uf);
00601 tmp = (1.0-omega)*e*kappa-ft*exp(-h*omega*kappa/uf);
00602 if (fabs(tmp) < ft*err)
00603 {
00604 omega += - tmp/dtmp;
00605 break;
00606 }
00607 omega += - tmp/dtmp;
00608 }
00609 }
00610
00611 if (fabs(tmp) > ft*err)
00612 {
00613 print_err("iteration of omega doesn't converge with required error\n"
00614 " elem %ld, ip %ld:\n"
00615 " kappa=%le, ft=%le, e=%le, uf=%le, h=%le",
00616 __FILE__, __LINE__, __func__, Mm->elip[ipp]+1, ipp+1, kappa, ft, e, uf, h);
00617 abort();
00618 }
00619
00620 if ((omega > 1.0) || (omega < 0.0))
00621 {
00622 print_err("computed omega=%le is out of prescribed range <0,1>\n"
00623 " elem %ld, ip %ld:\n"
00624 " kappa=%le, ft=%le, e=%le, uf=%le, h=%le",
00625 __FILE__, __LINE__, __func__, Mm->elip[ipp]+1, ipp+1, kappa, ft, e, uf, h);
00626 abort();
00627 }
00628
00629 return omega;
00630 }
00631
00632
00633
00634
00635
00636
00637
00638
00639
00640
00641
00642
00643
00644
00645
00646
00647
00648 double anisodam::dam_vol(long ipp, double y, double aa, double dvo)
00649 {
00650 double e, tmp, dam = 0.0;
00651 double aft, auft;
00652
00653 switch(damevf)
00654 {
00655 case brittle:
00656 e = Mm->give_actual_ym(ipp);
00657
00658 aft = Mm->give_actual_ft(ipp);
00659
00660 auft = uft/(aft/ft);
00661 dam = brittle_damage(ipp, y, e, aft, auft, dvo);
00662 break;
00663 case quasi_brittle:
00664 if (y < y0)
00665 break;
00666 tmp = log(y-y0);
00667 dam = aa*exp(tmp*b);
00668 tmp = 1.0/(1.0+dam);
00669 dam = 1.0-tmp;
00670 if (dam < 0.0)
00671 {
00672 fprintf(stdout, "\nDamage < 0.0 (ipp=%ld)\n", ipp);
00673 abort();
00674 }
00675 break;
00676 default:
00677 print_err(" unknown type of damage evolution function is required", __FILE__, __LINE__,__func__);
00678 }
00679 return dam;
00680 }
00681
00682
00683
00684
00685
00686
00687
00688
00689
00690
00691
00692
00693
00694
00695
00696
00697
00698
00699
00700
00701
00702 void anisodam::pdam_dev(long ipp, vector &pyc, vector &pyt, double aat, double aac, vector &pdamt, vector &pdamc)
00703 {
00704 long i;
00705 double tmp, ltmp;
00706 double e, af, auf;
00707
00708 fillv(0.0, pdamt);
00709 fillv(0.0, pdamc);
00710
00711 switch(damevf)
00712 {
00713 case brittle:
00714 e = Mm->give_actual_ym(ipp);
00715 for(i=0; i<3; i++)
00716 {
00717
00718
00719 af = Mm->give_actual_ft(ipp);
00720
00721 auf = uft/(af/ft);
00722 pdamt(i) = brittle_damage(ipp, pyt(i), e, af, auf, pdamt(i));
00723
00724
00725 af = Mm->give_actual_fc(ipp);
00726
00727 auf = ufc/(af/fc);
00728 pdamc(i) = brittle_damage(ipp, pyc(i), e, af, auf, pdamc(i));
00729 }
00730 break;
00731 case quasi_brittle:
00732 for(i=0; i<3; i++)
00733 {
00734
00735 if (pyt(i) > y0t)
00736 {
00737 ltmp = log(pyt(i) - y0t);
00738 pdamt(i) = aat*exp(ltmp*bt);
00739 tmp = 1.0/(1.0+pdamt(i));
00740 pdamt(i) = 1.0 - tmp;
00741 }
00742 if (pdamt(i) < 0.0)
00743 {
00744 fprintf(stdout, "\nDamage < 0.0 (ipp=%ld)\n", ipp);
00745 abort();
00746 }
00747
00748 if (pyc(i) > y0c)
00749 {
00750 ltmp = log(pyc(i) - y0c);
00751 pdamc(i) = aac*exp(ltmp*bc);
00752 tmp = 1.0/(1.0+pdamc(i));
00753 pdamc(i) = 1.0-tmp;
00754 }
00755 if (pdamc(i) < 0.0)
00756 {
00757 fprintf(stdout, "\nDamage < 0.0 (ipp=%ld)\n", ipp);
00758 abort();
00759 }
00760 }
00761 break;
00762 default:
00763 print_err(" unknown type of damage evolution function is required", __FILE__, __LINE__,__func__);
00764 }
00765 }
00766
00767
00768
00769
00770
00771
00772
00773
00774
00775
00776
00777
00778
00779
00780
00781
00782
00783
00784
00785
00786
00787 void anisodam::give_actual_param_a(long ipp, long ido, double &aa, double &aat, double &aac)
00788 {
00789 long eid;
00790 double e, h, tmp;
00791
00792 aa = Mm->ip[ipp].eqother[ido];
00793 aat = Mm->ip[ipp].eqother[ido+1];
00794 aac = Mm->ip[ipp].eqother[ido+2];
00795 if ((aa == 0.0) || (aac == 0.0) || (aat == 0.0))
00796 {
00797 if (cde == corr_off)
00798
00799 {
00800 aa = a;
00801 aat = at;
00802 aac = ac;
00803 Mm->ip[ipp].eqother[ido] = aa;
00804 Mm->ip[ipp].eqother[ido+1] = aat;
00805 Mm->ip[ipp].eqother[ido+2] = aac;
00806 return;
00807 }
00808 eid = Mm->elip[ipp];
00809 switch (Mt->give_dimension(eid))
00810 {
00811 case 1:
00812 h = Mt->give_length(eid);
00813 break;
00814 case 2:
00815 h = sqrt(Mt->give_area(eid));
00816 break;
00817 case 3:
00818 h = pow(Mt->give_volume(eid), 1.0/3.0);
00819 break;
00820 default:
00821 print_err("cannot determine dimension of element", __FILE__, __LINE__, __func__);
00822 }
00823 e = Mm->give_actual_ym(ipp);
00824
00825 tmp = (gf/h - y0/e)*b*e*sin(M_PI/b)/M_PI;
00826 if (tmp < 0.0)
00827 print_err("cannot express variable softening for volumetric damage", __FILE__, __LINE__, __func__);
00828 aa = pow(tmp, -b);
00829
00830 tmp = (gft/h - y0t/e)*bt*e*sin(M_PI/bt)/M_PI;
00831 if (tmp < 0.0)
00832 print_err("cannot express variable softening for tensile damage", __FILE__, __LINE__, __func__);
00833 aat = pow(tmp, -bt);
00834
00835 tmp = (gfc/h - y0c/e)*bc*e*sin(M_PI/bc)/M_PI;
00836 if (tmp < 0.0)
00837 print_err("cannot express variable softening for compressive damage", __FILE__, __LINE__, __func__);
00838 aac = pow(tmp, -bc);
00839 Mm->ip[ipp].eqother[ido] = aa;
00840 Mm->ip[ipp].eqother[ido+1] = aat;
00841 Mm->ip[ipp].eqother[ido+2] = aac;
00842 }
00843 }
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857 double anisodam::give_actual_ft(long ipp)
00858 {
00859 return ft;
00860 }
00861
00862
00863
00864
00865
00866
00867
00868
00869
00870
00871
00872
00873
00874 double anisodam::give_actual_fc(long ipp)
00875 {
00876 return fc;
00877 }
00878
00879
00880
00881
00882
00883
00884
00885
00886
00887
00888
00889
00890
00891
00892 void anisodam::initvalues(long ipp, long ido)
00893 {
00894 double aa, aat, aac;
00895 give_actual_param_a(ipp, ido, aa, aac, aat);
00896 }
00897
00898
00899
00900
00901
00902
00903
00904
00905
00906
00907
00908
00909
00910
00911 void anisodam::matstiff (matrix &d,long ipp,long ido)
00912 {
00913
00914
00915 switch (Mp->nlman->stmat)
00916 {
00917 case initial_stiff:
00918 elmatstiff(d, ipp);
00919 break;
00920 case tangent_stiff:
00921
00922 tmatstiff(d, ipp, ido);
00923 break;
00924 default:
00925 print_err("unknown type of stiffness matrix is required", __FILE__, __LINE__, __func__);
00926 }
00927 }
00928
00929
00930
00931
00932
00933
00934
00935
00936
00937
00938
00939
00940
00941 void anisodam::elmatstiff (matrix &d,long ipp)
00942 {
00943 double e, e0;
00944 long idem = Mm->ip[ipp].gemid();
00945
00946 Mm->elmatstiff (d,ipp);
00947 e = Mm->give_actual_ym(ipp);
00948
00949 if (Mm->ip[ipp].tm[idem] != elisomat)
00950 print_err("invalid type of elastic material is required", __FILE__, __LINE__, __func__);
00951
00952 e0 = Mm->eliso[Mm->ip[ipp].idm[idem]].e;
00953 cmulm(e/e0, d);
00954 }
00955
00956
00957
00958
00959
00960
00961
00962
00963
00964
00965
00966
00967
00968
00969 void anisodam::tmatstiff (matrix &d,long ipp, long ido)
00970 {
00971 double e, nu, g, k0, dv, tmpt, tmpc;
00972 long i, j, k, l, m, idem, ime, ncompstr = Mm->ip[ipp].ncompstr;
00973 vector epsn(ncompstr), peps(3), pdt(3), pdc(3);
00974 matrix epst(3,3), t(3,3);
00975 matrix dt(3,3), dc(3,3), omega_t(3,3), omega_c(3,3), dpdi_deps(3,3), dpdj_deps(3,3);
00976 bmatrix d_c(3,3), d_t(3,3), tmp_dc(3,3), tmp_dt(3,3);
00977 matrix d_tot(6,6), c_tot(6,6), c_red(3,3);
00978 double epsv;
00979
00980 e = Mm->give_actual_ym(ipp);
00981 idem = Mm->ip[ipp].gemid();
00982 ime = Mm->ip[ipp].idm[idem];
00983 nu = Mm->eliso[ime].nu;
00984 if (Mm->ip[ipp].ssst == planestress)
00985 Mm->ip[ipp].strain[3] = -nu / (1.0 - nu) * (Mm->ip[ipp].strain[0]+Mm->ip[ipp].strain[1]);
00986
00987 if (Mm->ip[ipp].tm[idem] != elisomat)
00988 print_err("invalid type of elastic material is required", __FILE__, __LINE__, __func__);
00989
00990 nu = Mm->eliso[Mm->ip[ipp].idm[idem]].nu;
00991 g = e/(2.0*(1.0+nu));
00992 k0 = e/(1.0-2.0*nu);
00993 dv = Mm->ip[ipp].other[ido+3];
00994 for (i=0; i<3; i++)
00995 {
00996 pdt[i] = Mm->ip[ipp].eqother[ido+3+1+i];
00997 pdc[i] = Mm->ip[ipp].eqother[ido+3+1+3+i];
00998 }
00999
01000 for (i=0;i<ncompstr;i++)
01001 epsn[i] = Mm->ip[ipp].strain[i];
01002 vector_tensor(epsn, epst, Mm->ip[ipp].ssst, strain);
01003 princ_val (epst,peps,t,nijac,limit,Mp->zero,3,1);
01004 epsv = peps[0] + peps[1] + peps[2];
01005 epsv /= 3.0;
01006 fillm(0.0, d);
01007
01008 for(i=0; i<3; i++)
01009 {
01010 d_tot[i][i] = (1/3.0);
01011 if (epsv > 0.0)
01012 d_tot[i][i] -= dv*1/3.0;
01013 d_tot[i][i] *= (k0-2.0*g);
01014 }
01015 for(i=0; i<3; i++)
01016 {
01017 dt[i][i] = 1.0-pdt[i];
01018 dc[i][i] = 1.0-pdc[i];
01019 }
01020 lgmatrixtransf(dt, t);
01021 lgmatrixtransf(dc, t);
01022 f_tensor(dt, &sqrt, t, omega_t);
01023 f_tensor(dc, &sqrt, t, omega_c);
01024
01025 for(i=0; i<3; i++)
01026 {
01027 for(j=0; j<3; j++)
01028 {
01029 allocm(3, 3, d_c[i][j]);
01030 allocm(3, 3, d_t[i][j]);
01031 allocm(3, 3, tmp_dc[i][j]);
01032 allocm(3, 3, tmp_dt[i][j]);
01033 }
01034 }
01035 d_c.gen_indices();
01036 d_t.gen_indices();
01037 tmp_dc.gen_indices();
01038 tmp_dt.gen_indices();
01039
01040
01041 for (m=0; m<3; m++)
01042 {
01043 for(i=0; i<3; i++)
01044 {
01045 for(j=0; j<3; j++)
01046 {
01047 for(k=0; k<3; k++)
01048 {
01049 for(l=0; l<3; l++)
01050 {
01051 if (dpdir_da(t, peps, m, i, dpdi_deps))
01052 {
01053 elmatstiff(d, ipp);
01054 return;
01055 }
01056 dpdir_da(t, peps, m, j, dpdj_deps);
01057 {
01058 elmatstiff(d, ipp);
01059 return;
01060 }
01061 if (peps[m] >= 0.0)
01062 {
01063 tmpt = peps[m]*(dpdj_deps[k][l]*t[i][m]+dpdi_deps[k][l]*t[j][m]);
01064 tmpt += t[i][m]*t[j][m]*t[k][m]*t[l][m];
01065 }
01066 if (peps[m] < 0.0)
01067 {
01068 tmpc = peps[m]*(dpdj_deps[k][l]*t[i][m]+dpdi_deps[k][l]*t[j][m]);
01069 tmpc += t[i][m]*t[j][m]*t[k][m]*t[l][m];
01070 }
01071 d_c[i][j][k][l] += 2.0*g*tmpc;
01072 d_t[i][j][k][l] += 2.0*g*tmpt;
01073 }
01074 }
01075 }
01076 }
01077 }
01078
01079 fillm(0.0, tmp_dc);
01080 fillm(0.0, tmp_dt);
01081 for(i=0; i<3; i++)
01082 {
01083 for(j=0; j<3; j++)
01084 {
01085 for(k=0; k<3; k++)
01086 {
01087 for(l=0; l<3; l++)
01088 {
01089 for(m=0; m<3; m++)
01090 {
01091 tmp_dc[i][j][k][l] += omega_c[i][m] * d_c[m][j][k][l];
01092 tmp_dt[i][j][k][l] += omega_t[i][m] * d_t[m][j][k][l];
01093 }
01094 }
01095 }
01096 }
01097 }
01098 fillm(0.0, d_c);
01099 fillm(0.0, d_t);
01100 for(i=0; i<3; i++)
01101 {
01102 for(j=0; j<3; j++)
01103 {
01104 for(k=0; k<3; k++)
01105 {
01106 for(l=0; l<3; l++)
01107 {
01108 for(m=0; m<3; m++)
01109 {
01110 d_c[i][j][k][l] += omega_c[m][j] * tmp_dc[i][m][k][l];
01111 d_t[i][j][k][l] += omega_t[m][j] * tmp_dt[i][m][k][l];
01112 }
01113 }
01114 }
01115 }
01116 }
01117
01118 for(i=0; i<3; i++)
01119 {
01120 for(j=0; j<3; j++)
01121 d_tot[i][j] += d_c[i][i][j][j] + d_t[i][i][j][j];
01122 }
01123 for(i=0; i<3; i++)
01124 {
01125 d_tot[i][3] += d_c[i][i][1][2] + d_t[i][i][1][2];
01126 d_tot[i][4] += d_c[i][i][0][2] + d_t[i][i][0][2];
01127 d_tot[i][5] += d_c[i][i][0][1] + d_t[i][i][0][1];
01128
01129 d_tot[3][i] += d_c[1][2][i][i] + d_t[1][2][i][i];
01130 d_tot[4][i] += d_c[0][2][i][i] + d_t[0][2][i][i];
01131 d_tot[5][i] += d_c[0][1][i][i] + d_t[0][1][i][i];
01132 }
01133 d_tot[3][3] += d_c[1][2][1][2] + d_t[1][2][1][2];
01134 d_tot[3][4] += d_c[1][2][0][2] + d_t[1][2][0][2];
01135 d_tot[3][5] += d_c[1][2][0][1] + d_t[1][2][0][1];
01136
01137 d_tot[4][3] += d_c[0][2][1][2] + d_t[0][2][1][2];
01138 d_tot[4][4] += d_c[0][2][0][2] + d_t[0][2][0][2];
01139 d_tot[4][5] += d_c[0][2][0][1] + d_t[0][2][0][1];
01140
01141 d_tot[5][3] += d_c[0][1][1][2] + d_t[0][1][1][2];
01142 d_tot[5][4] += d_c[0][1][0][2] + d_t[0][1][0][2];
01143 d_tot[5][5] += d_c[0][1][0][1] + d_t[0][1][0][1];
01144
01145 switch (Mm->ip[ipp].ssst)
01146 {
01147 case planestress:
01148 {
01149 invm(d_tot, c_tot, 1e-8);
01150 for(i=0; i<2; i++)
01151 {
01152 for(j=0; j<2; j++)
01153 c_red[i][j] = c_tot[i][j];
01154 }
01155 c_red[0][2] = c_tot[0][5];
01156 c_red[1][2] = c_tot[1][5];
01157 c_red[2][2] = c_tot[5][5];
01158
01159 c_red[2][0] = c_tot[5][0];
01160 c_red[2][1] = c_tot[5][1];
01161 invm(c_red, d, 1e-8);
01162 break;
01163 }
01164 case planestrain:
01165 {
01166 for(i=0; i<2; i++)
01167 {
01168 for(j=0; j<2; j++)
01169 d[i][j] = d_tot[i][j];
01170 }
01171 d[0][2] = d_tot[0][5];
01172 d[1][2] = d_tot[1][5];
01173 d[2][2] = d_tot[5][5];
01174
01175 d[2][0] = d_tot[5][0];
01176 d[2][1] = d_tot[5][1];
01177 break;
01178 }
01179 case spacestress:
01180 copym(d_tot, d);
01181 break;
01182 default:
01183 print_err("Unknown type of stress/strain state is required", __FILE__, __LINE__, __func__);
01184 abort();
01185 }
01186 }
01187
01188
01189
01190
01191
01192
01193
01194
01195
01196
01197
01198
01199
01200
01201
01202
01203 void anisodam::nlstresses (long ipp, long im, long ido)
01204 {
01205 long ime, idem, i;
01206 long ncompstr=Mm->ip[ipp].ncompstr;
01207
01208
01209 double aa, aat, aac;
01210
01211
01212 vector pyto(3), pyco(3);
01213
01214 double dvo;
01215 vector pdto(3), pdco(3);
01216
01217
01218 vector pdyt(3), pdyc(3);
01219
01220 double yn;
01221 vector pytn(3), pycn(3);
01222
01223 double dv;
01224 vector pdt(3), pdc(3);
01225
01226 vector epsn(ncompstr);
01227
01228 matrix epst(3,3), t(3,3);
01229
01230 matrix sigt(3,3);
01231 vector sigma(ncompstr);
01232
01233 double k0, g0, epsv, e, nu;
01234
01235
01236 vector peps(3), psig(3), lft(3), lfc(3),pepsp(3);
01237 double tmp, h;
01238
01239 long eid;
01240
01241
01242
01243
01244
01245
01246 give_actual_param_a(ipp, ido, aa, aat, aac);
01247
01248 e = Mm->give_actual_ym(ipp);
01249 idem = Mm->ip[ipp].gemid();
01250 ime = Mm->ip[ipp].idm[idem];
01251 nu = Mm->eliso[ime].nu;
01252 if (Mm->ip[ipp].ssst == planestress)
01253 Mm->ip[ipp].strain[3] = -nu / (1.0 - nu) * (Mm->ip[ipp].strain[0]+Mm->ip[ipp].strain[1]);
01254
01255
01256 for (i=0;i<ncompstr;i++)
01257 epsn[i] = Mm->ip[ipp].strain[i];
01258
01259
01260 dvo = Mm->ip[ipp].eqother[ido+3];
01261 for (i=0; i<3; i++)
01262 {
01263
01264 pdto[i] = Mm->ip[ipp].eqother[ido+3+1+i];
01265 pdco[i] = Mm->ip[ipp].eqother[ido+3+1+3+i];
01266 }
01267
01268
01269 vector_tensor(epsn, epst, Mm->ip[ipp].ssst, strain);
01270 princ_val (epst,peps,t,nijac,limit,Mp->zero,3, 1);
01271
01272
01273
01274
01275
01276
01277
01278
01279
01280
01281
01282
01283
01284
01285
01286
01287
01288
01289
01290 yn = damdrvforce_vol(nu, peps);
01291 damdrvforce_dev(nu, peps, pytn, pycn);
01292
01293
01294 dv = dam_vol(ipp, yn, aa, dvo);
01295 if (dv < dvo)
01296 dv = dvo;
01297 copyv(pdto, pdt);
01298 copyv(pdco, pdc);
01299 pdam_dev(ipp, pycn, pytn, aat, aac, pdt, pdc);
01300 for (i=0; i<3; i++)
01301 {
01302 if (pdt[i] < pdto[i])
01303 pdt[i] = pdto[i];
01304 if (pdc[i] < pdco[i])
01305 pdc[i] = pdco[i];
01306 }
01307
01308
01309 k0 = e;
01310 k0 /= (1-2.0*nu);
01311 g0 = e;
01312 g0 /= 2.0*(1.0+nu);
01313 epsv = peps[0]+peps[1]+peps[2];
01314 epsv /= 3.0;
01315 psig[0] = k0 - 2.0*g0;
01316 if (epsv > 0.0)
01317 epsv -= dv*epsv;
01318
01319 psig[0] *= epsv;
01320 psig[1] = psig[2] = psig[0];
01321 for (i=0; i<3; i++)
01322 {
01323 tmp = 1.0;
01324 if (peps[i] > 0.0)
01325 tmp -= pdt[i];
01326 if (peps[i] < 0.0)
01327 tmp -= pdc[i];
01328 psig[i] += 2.0*g0*tmp*peps[i];
01329 }
01330
01331
01332 fillm(0.0, sigt);
01333 for (i = 0; i < 3; i++)
01334 sigt[i][i] = psig[i];
01335 lgmatrixtransf(sigt, t);
01336 tensor_vector(sigma, sigt, Mm->ip[ipp].ssst, stress);
01337 if (Mm->ip[ipp].ssst == planestress)
01338 sigma[3] = 0.0;
01339
01340
01341 for (i=0;i<ncompstr;i++)
01342 Mm->ip[ipp].stress[i]=sigma[i];
01343 Mm->ip[ipp].other[ido+3] = dv;
01344 for (i=0; i<3; i++)
01345 {
01346
01347 Mm->ip[ipp].other[ido+3+1+i] = pdt[i];
01348 Mm->ip[ipp].other[ido+3+1+3+i] = pdc[i];
01349 }
01350
01351
01352 fillm(0.0, sigt);
01353 for (i = 0; i < 3; i++)
01354 sigt[i][i] = pdt[i];
01355 lgmatrixtransf(sigt, t);
01356 for (i = 0; i < 3; i++)
01357 Mm->ip[ipp].other[ido+10+i] = sigt[i][i];
01358 Mm->ip[ipp].other[ido+13] = sigt[1][2];
01359 Mm->ip[ipp].other[ido+14] = sigt[0][2];
01360 Mm->ip[ipp].other[ido+15] = sigt[0][1];
01361
01362 for (i = 0; i < 3; i++)
01363 Mm->ip[ipp].other[ido+16+i] = sigt[i][i]+dv;
01364 Mm->ip[ipp].other[ido+19] = sigt[1][2];
01365 Mm->ip[ipp].other[ido+20] = sigt[0][2];
01366 Mm->ip[ipp].other[ido+21] = sigt[0][1];
01367
01368 eid = Mm->elip[ipp];
01369 switch (Mt->give_dimension(eid))
01370 {
01371 case 1:
01372 h = Mt->give_length(eid);
01373 break;
01374 case 2:
01375 h = sqrt(Mt->give_area(eid));
01376 break;
01377 case 3:
01378 h = pow(Mt->give_volume(eid), 1.0/3.0);
01379 break;
01380 default:
01381 print_err("unknown dimension of element", __FILE__, __LINE__, __func__);
01382 }
01383
01384 Mm->ip[ipp].other[ido+22+i] = (sqrt(yn)-(psig[0]+psig[1]+psig[2])/(3.0*e))*h;
01385
01386 for (i = 0; i < 3; i++)
01387 Mm->ip[ipp].other[ido+23+i] = (sqrt(pytn[i])-psig[i]/e)*h;
01388 }
01389
01390
01391
01392
01393
01394
01395
01396
01397
01398
01399
01400
01401
01402
01403
01404 void anisodam::updateval (long ipp,long im,long ido)
01405 {
01406 long i;
01407 long n = Mm->givencompeqother(ipp,im);
01408
01409 for (i=0;i<n;i++)
01410 Mm->ip[ipp].eqother[ido+i]=Mm->ip[ipp].other[ido+i];
01411 }