00001 #include "anisodamrot.h"
00002 #include "matrix.h"
00003 #include "vector.h"
00004 #include "elastisomat.h"
00005 #include "global.h"
00006 #include "intpoints.h"
00007 #include "tensor.h"
00008 #include "vecttens.h"
00009 #include <math.h>
00010 #include <stdlib.h>
00011
00012 #define nijac 20
00013 #define limit 1.0e-16
00014 #ifndef M_PI
00015 #define M_PI 3.14159265358979323846
00016 #endif
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027 anisodamrot::anisodamrot (void)
00028 {
00029 cde = corr_off;
00030 ac = 0.0; bc = 0.0; y0c = 0.0; at = 0.0; bt = 0.0; y0t = 0.0;
00031 a = 0.0; b = 0.0; y0 = 0.0; gfc = 0.0; gft = 0.0; gf = 0.0;
00032 ft = uft = 0.0;
00033 fc = ufc = 0.0;
00034 }
00035
00036
00037
00038
00039
00040
00041
00042
00043 anisodamrot::~anisodamrot (void)
00044 {
00045
00046 }
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061 void anisodamrot::read (XFILE *in)
00062 {
00063 xfscanf (in, "%k%m", "dam_evol_f", &dam_evolfunc_kwdset, (int *)&damevf);
00064 switch(damevf)
00065 {
00066 case brittle:
00067 xfscanf (in,"%k%lf %k%lf", "f_t", &ft, "u_ft", &uft);
00068 xfscanf (in,"%k%lf %k%lf", "f_c", &fc, "u_fc", &ufc);
00069 xfscanf (in,"%k%m", "corr_dissip", &corr_disip_en_kwdset, (int *)(&cde));
00070 if (cde != corr_off)
00071 sra.read (in);
00072 break;
00073 case quasi_brittle:
00074 xfscanf (in, "%k%m", "corr_dissip", &corr_disip_en_kwdset, (int *)(&cde));
00075 switch (cde)
00076 {
00077 case corr_off:
00078 xfscanf (in, "%k%lf %k%lf %k%lf", "A_c", &ac, "B_c", &bc, "Y_0c", &y0c);
00079 xfscanf (in, "%k%lf %k%lf %k%lf", "A_t", &at, "B_t", &bt, "Y_0t", &y0t);
00080 xfscanf (in, "%k%lf %k%lf %k%lf", "a", &a, "b", &b, "y0", &y0);
00081 break;
00082 case corr_on:
00083 xfscanf (in, "%k%lf %k%lf %k%lf", "G_fc", &gfc, "B_c", &bc, "Y_0c", &y0c);
00084 xfscanf (in, "%k%lf %k%lf %k%lf", "G_ft", &gft, "B_t", &bt, "Y_0t", &y0t);
00085 xfscanf (in, "%k%lf %k%lf %k%lf", "G_f", &gf, "b", &b, "y0", &y0);
00086 break;
00087 default:
00088 print_err(" unknown type of correction of disipated energy is required", __FILE__, __LINE__,__func__);
00089 }
00090 break;
00091 default:
00092 print_err(" unknown type of damage evolution function is required", __FILE__, __LINE__,__func__);
00093 }
00094 }
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108 double anisodamrot::damdrvforce_vol(double nu, vector &peps)
00109 {
00110 long i;
00111 double ev = 0.0;
00112 double k0_prime, g0_prime, ret;
00113
00114 for (i = 0; i < 3; i++)
00115 ev += peps(i);
00116 ev /= 3.0;
00117
00118
00119 if (ev <= 0.0)
00120 return 0.0;
00121
00122
00123 k0_prime = 1.0/(1.0-2.0*nu);
00124 g0_prime = 1.0/(2.0*(1.0+nu));
00125 ret = 3.0/2.0*(k0_prime-2.0*g0_prime)*ev*ev;
00126 return ret;
00127 }
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144 void anisodamrot::damdrvforce_dev(double nu, vector &peps, vector &pyt, vector &pyc)
00145 {
00146 long i;
00147 double g0_prime;
00148
00149
00150 g0_prime = 1.0/(2.0*(1.0+nu));
00151 fillv(0.0, pyt);
00152 fillv(0.0, pyc);
00153 for (i = 0; i < 3; i++)
00154 {
00155 if (peps(i) > 0.0)
00156 pyt[i] = g0_prime*peps[i]*peps[i];
00157 if (peps(i) < 0.0)
00158 pyc[i] = g0_prime*peps[i]*peps[i];
00159 }
00160 }
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178 double anisodamrot::loadfuncvol(long ipp, double nu, vector &peps, double d, double aa)
00179 {
00180 double y, f, kappa;
00181 double tmp1, tmp2;
00182 double e, aft, auft;
00183
00184 y = damdrvforce_vol(nu, peps);
00185 switch(damevf)
00186 {
00187 case brittle:
00188 kappa = sqrt(y);
00189 if (kappa <= y0)
00190 f = -1.0;
00191 else
00192 {
00193 e = Mm->give_actual_ym(ipp);
00194
00195 aft = Mm->give_actual_ft(ipp);
00196
00197 auft = uft/(aft/ft);
00198 tmp1 = -(kappa-aft/e)/(auft-aft/e);
00199 f = 1.0 - aft/(e*kappa)*exp(tmp1) - d;
00200 }
00201 break;
00202 case quasi_brittle:
00203 if (y <= y0)
00204 f = -1.0;
00205 else
00206 {
00207
00208 tmp1 = y - y0;
00209 tmp2 = log(tmp1);
00210 tmp1 = exp(b*tmp2);
00211
00212 f = (1.0 - d)*(1.0+aa*tmp1) - 1.0;
00213 }
00214 break;
00215 default:
00216 f = 0.0;
00217 print_err(" unknown type of damage evolution function is required", __FILE__, __LINE__,__func__);
00218 }
00219 return f;
00220 }
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244 void anisodamrot::loadfuncdev(long ipp, double nu, vector &peps, vector &damt, vector &damc,
00245 double aat, double aac, vector &lft, vector &lfc)
00246 {
00247 long i;
00248 double e, af, auf, kappa, tmp1, tmp2;
00249 vector pyt(3), pyc(3);
00250
00251 damdrvforce_dev(nu, peps, pyt, pyc);
00252 switch(damevf)
00253 {
00254 case brittle:
00255 for (i=0; i<3; i++)
00256 {
00257
00258 kappa = sqrt(pyt(i));
00259 if (kappa <= y0t)
00260 lft(i) = -1.0;
00261 else
00262 {
00263 e = Mm->give_actual_ym(ipp);
00264
00265 af = Mm->give_actual_ft(ipp);
00266
00267 auf = uft/(af/ft);
00268 tmp1 = -(kappa-af/e)/(auf-af/e);
00269 lft(i) = 1.0 - af/(e*kappa)*exp(tmp1) - damt(i);
00270 }
00271
00272 kappa = sqrt(pyc(i));
00273 if (kappa <= y0c)
00274 lfc(i) = -1.0;
00275 else
00276 {
00277 e = Mm->give_actual_ym(ipp);
00278
00279 af = Mm->give_actual_fc(ipp);
00280
00281 auf = ufc/(af/fc);
00282 tmp1 = -(kappa-af/e)/(auf-af/e);
00283 lfc(i) = 1.0 - af/(e*kappa)*exp(tmp1) - damc(i);
00284 }
00285 }
00286 break;
00287 case quasi_brittle:
00288 for (i=0; i<3; i++)
00289 {
00290
00291 if (pyt(i) <= y0t)
00292 lft(i) = -1.0;
00293 else
00294 {
00295
00296 tmp1 = pyt(i) - y0t;
00297 tmp2 = log(tmp1);
00298 tmp1 = exp(bt*tmp2);
00299
00300 lft(i) = (1.0 - damt(i))*(1.0+aat*tmp1) - 1.0;
00301 }
00302
00303 if (pyc(i) <= y0c)
00304 lfc(i) = -1.0;
00305 else
00306 {
00307
00308 tmp1 = pyc(i) - y0c;
00309 tmp2 = log(tmp1);
00310 tmp1 = exp(bc*tmp2);
00311
00312 lfc(i) = (1.0 - damc(i))*(1.0+aac*tmp1) - 1.0;
00313 }
00314 }
00315 break;
00316 default:
00317 print_err(" unknown type of damage evolution function is required", __FILE__, __LINE__,__func__);
00318 }
00319 return;
00320 }
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338 double anisodamrot::dam_vol(long ipp, double y, double aa, double dvo)
00339 {
00340 double e, tmp, dam = 0.0;
00341 double aft, auft;
00342
00343 switch(damevf)
00344 {
00345 case brittle:
00346 e = Mm->give_actual_ym(ipp);
00347
00348 aft = Mm->give_actual_ft(ipp);
00349
00350 auft = uft/(aft/ft);
00351 dam = brittle_damage(ipp, y, e, aft, auft, dvo);
00352 break;
00353 case quasi_brittle:
00354 if (y < y0)
00355 break;
00356 tmp = log(y-y0);
00357 dam = aa*exp(tmp*b);
00358 tmp = 1.0/(1.0+dam);
00359 dam = 1.0-tmp;
00360 if (dam < 0.0)
00361 {
00362 fprintf(stdout, "\nDamage < 0.0 (ipp=%ld)\n", ipp);
00363 abort();
00364 }
00365 break;
00366 default:
00367 print_err(" unknown type of damage evolution function is required", __FILE__, __LINE__,__func__);
00368 }
00369 return dam;
00370 }
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391 void anisodamrot::pdam_dev(long ipp, vector &pyc, vector &pyt, double aat, double aac, vector &pdamt, vector &pdamc)
00392 {
00393 long i;
00394 double tmp, ltmp;
00395 double e, af, auf;
00396
00397 fillv(0.0, pdamt);
00398 fillv(0.0, pdamc);
00399
00400 switch(damevf)
00401 {
00402 case brittle:
00403 e = Mm->give_actual_ym(ipp);
00404 for(i=0; i<3; i++)
00405 {
00406
00407
00408 af = Mm->give_actual_ft(ipp);
00409
00410 auf = uft/(af/ft);
00411 pdamt(i) = brittle_damage(ipp, pyt(i), e, af, auf, pdamt(i));
00412
00413
00414 af = Mm->give_actual_fc(ipp);
00415
00416 auf = ufc/(af/fc);
00417 pdamc(i) = brittle_damage(ipp, pyc(i), e, af, auf, pdamc(i));
00418 }
00419 break;
00420 case quasi_brittle:
00421 for(i=0; i<3; i++)
00422 {
00423
00424 if (pyt(i) > y0t)
00425 {
00426 ltmp = log(pyt(i) - y0t);
00427 pdamt(i) = aat*exp(ltmp*bt);
00428 tmp = 1.0/(1.0+pdamt(i));
00429 pdamt(i) = 1.0 - tmp;
00430 }
00431 if (pdamt(i) < 0.0)
00432 {
00433 fprintf(stdout, "\nDamage < 0.0 (ipp=%ld)\n", ipp);
00434 abort();
00435 }
00436
00437 if (pyc(i) > y0c)
00438 {
00439 ltmp = log(pyc(i) - y0c);
00440 pdamc(i) = aac*exp(ltmp*bc);
00441 tmp = 1.0/(1.0+pdamc(i));
00442 pdamc(i) = 1.0-tmp;
00443 }
00444 if (pdamc(i) < 0.0)
00445 {
00446 fprintf(stdout, "\nDamage < 0.0 (ipp=%ld)\n", ipp);
00447 abort();
00448 }
00449 }
00450 break;
00451 default:
00452 print_err(" unknown type of damage evolution function is required", __FILE__, __LINE__,__func__);
00453 }
00454 }
00455
00456
00457
00458
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474 double anisodamrot::brittle_damage(long ipp, double y, double e, double f, double uf, double omegao)
00475 {
00476 double err,omega = 0.0, dtmp, tmp, h, indam, kappa;
00477 long i,ni, eid;
00478
00479 ni=sra.give_ni ();
00480 err=sra.give_err ();
00481
00482 indam = f/e;
00483
00484 kappa = sqrt(y);
00485 if (kappa < indam)
00486
00487 return 0.0;
00488 if ((Mm->ip[ipp].hmt & 2) > 0)
00489 {
00490 omega = 1.0-(f/(e*kappa))*exp(-(kappa-f/e)/(uf-f/e));
00491 return omega;
00492 }
00493
00494 eid = Mm->elip[ipp];
00495 switch (Mt->give_dimension(eid))
00496 {
00497 case 1:
00498 h = Mt->give_length(eid);
00499 break;
00500 case 2:
00501 h = sqrt(Mt->give_area(eid));
00502 break;
00503 case 3:
00504 h = pow(Mt->give_volume(eid), 1.0/3.0);
00505 break;
00506 default:
00507 print_err("unknown dimension of element", __FILE__, __LINE__, __func__);
00508 }
00509 if (cde == corr_off)
00510
00511 {
00512 omega = 1.0-(f/(e*kappa))*exp(-(kappa-f/e)/(uf-f/e));
00513 return omega;
00514 }
00515
00516 omega = 1.0;
00517 for (i = 0; i < ni; i++)
00518 {
00519 dtmp = -e*kappa+f/uf*h*kappa*exp(-h*omega*kappa/uf);
00520 tmp = (1.0-omega)*e*kappa-f*exp(-h*omega*kappa/uf);
00521 if (fabs(tmp) < ft*err)
00522 {
00523 omega += - tmp/dtmp;
00524 break;
00525 }
00526 omega += - tmp/dtmp;
00527 }
00528 if ((omega > 1.0) || (omega < 0.0) || (fabs(tmp) > ft*err))
00529
00530
00531 {
00532 omega = omegao;
00533 for (i = 0; i < ni; i++)
00534 {
00535 dtmp = -e*kappa+ft/uf*h*kappa*exp(-h*omega*kappa/uf);
00536 tmp = (1.0-omega)*e*kappa-ft*exp(-h*omega*kappa/uf);
00537 if (fabs(tmp) < ft*err)
00538 {
00539 omega += - tmp/dtmp;
00540 break;
00541 }
00542 omega += - tmp/dtmp;
00543 }
00544 }
00545
00546 if (fabs(tmp) > ft*err)
00547 {
00548 print_err("iteration of omega doesn't converge with required error\n"
00549 " elem %ld, ip %ld:\n"
00550 " kappa=%le, ft=%le, e=%le, uf=%le, h=%le",
00551 __FILE__, __LINE__, __func__, Mm->elip[ipp]+1, ipp+1, kappa, ft, e, uf, h);
00552 abort();
00553 }
00554
00555 if ((omega > 1.0) || (omega < 0.0))
00556 {
00557 print_err("computed omega=%le is out of prescribed range <0,1>\n"
00558 " elem %ld, ip %ld:\n"
00559 " kappa=%le, ft=%le, e=%le, uf=%le, h=%le",
00560 __FILE__, __LINE__, __func__, Mm->elip[ipp]+1, ipp+1, kappa, ft, e, uf, h);
00561 abort();
00562 }
00563
00564 return omega;
00565 }
00566
00567
00568
00569
00570
00571
00572
00573
00574
00575
00576
00577
00578
00579
00580
00581
00582
00583
00584
00585
00586 void anisodamrot::give_actual_param_a(long ipp, long ido, double &aa, double &aat, double &aac)
00587 {
00588 long eid;
00589 double e, h, tmp;
00590
00591 aa = Mm->ip[ipp].eqother[ido];
00592 aat = Mm->ip[ipp].eqother[ido+1];
00593 aac = Mm->ip[ipp].eqother[ido+2];
00594 if ((aa == 0.0) || (aac == 0.0) || (aat == 0.0))
00595 {
00596 if (cde == corr_off)
00597
00598 {
00599 aa = a;
00600 aat = at;
00601 aac = ac;
00602 Mm->ip[ipp].eqother[ido] = aa;
00603 Mm->ip[ipp].eqother[ido+1] = aat;
00604 Mm->ip[ipp].eqother[ido+2] = aac;
00605 return;
00606 }
00607 eid = Mm->elip[ipp];
00608 switch (Mt->give_dimension(eid))
00609 {
00610 case 1:
00611 h = Mt->give_length(eid);
00612 break;
00613 case 2:
00614 h = sqrt(Mt->give_area(eid));
00615 break;
00616 case 3:
00617 h = pow(Mt->give_volume(eid), 1.0/3.0);
00618 break;
00619 default:
00620 print_err("cannot determine dimension of element", __FILE__, __LINE__, __func__);
00621 }
00622 e = Mm->give_actual_ym(ipp);
00623
00624 tmp = (gf/h - y0/e)*b*e*sin(M_PI/b)/M_PI;
00625 if (tmp < 0.0)
00626 print_err("cannot express variable softening for volumetric damage", __FILE__, __LINE__, __func__);
00627 aa = pow(tmp, -b);
00628
00629 tmp = (gft/h - y0t/e)*bt*e*sin(M_PI/bt)/M_PI;
00630 if (tmp < 0.0)
00631 print_err("cannot express variable softening for tensile damage", __FILE__, __LINE__, __func__);
00632 aat = pow(tmp, -bt);
00633
00634 tmp = (gfc/h - y0c/e)*bc*e*sin(M_PI/bc)/M_PI;
00635 if (tmp < 0.0)
00636 print_err("cannot express variable softening for compressive damage", __FILE__, __LINE__, __func__);
00637 aac = pow(tmp, -bc);
00638 Mm->ip[ipp].eqother[ido] = aa;
00639 Mm->ip[ipp].eqother[ido+1] = aat;
00640 Mm->ip[ipp].eqother[ido+2] = aac;
00641 }
00642 }
00643
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653
00654
00655
00656 double anisodamrot::give_actual_ft(long ipp)
00657 {
00658 return ft;
00659 }
00660
00661
00662
00663
00664
00665
00666
00667
00668
00669
00670
00671
00672
00673 double anisodamrot::give_actual_fc(long ipp)
00674 {
00675 return fc;
00676 }
00677
00678
00679
00680
00681
00682
00683
00684
00685
00686
00687
00688
00689
00690
00691 void anisodamrot::initvalues(long ipp, long ido)
00692 {
00693 double aa, aat, aac;
00694 give_actual_param_a(ipp, ido, aa, aac, aat);
00695 }
00696
00697
00698
00699
00700
00701
00702
00703
00704
00705
00706
00707
00708
00709
00710 void anisodamrot::matstiff (matrix &d,long ipp,long ido)
00711 {
00712
00713
00714 switch (Mp->nlman->stmat)
00715 {
00716 case initial_stiff:
00717 elmatstiff(d, ipp);
00718 break;
00719 case tangent_stiff:
00720 elmatstiff(d, ipp);
00721
00722 break;
00723 default:
00724 print_err("unknown type of stiffness matrix is required", __FILE__, __LINE__, __func__);
00725 }
00726 }
00727
00728
00729
00730
00731
00732
00733
00734
00735
00736
00737
00738
00739
00740 void anisodamrot::elmatstiff (matrix &d,long ipp)
00741 {
00742 double e, e0;
00743 long idem = Mm->ip[ipp].gemid();
00744
00745 Mm->elmatstiff (d,ipp);
00746 e = Mm->give_actual_ym(ipp);
00747
00748 if (Mm->ip[ipp].tm[idem] != elisomat)
00749 print_err("invalid type of elastic material is required", __FILE__, __LINE__, __func__);
00750
00751 e0 = Mm->eliso[Mm->ip[ipp].idm[idem]].e;
00752 cmulm(e/e0, d);
00753 }
00754
00755
00756
00757
00758
00759
00760
00761
00762
00763
00764
00765
00766
00767
00768
00769
00770
00771 void anisodamrot::nlstresses (long ipp, long im, long ido)
00772 {
00773 long i;
00774 long ncompstr=Mm->ip[ipp].ncompstr;
00775
00776
00777 double aa, aat, aac;
00778
00779 double dvo;
00780
00781 double dv;
00782
00783 matrix dt(3,3), dc(3,3);
00784
00785 matrix eps_t(3,3), eps_c(3,3);
00786
00787 double yn;
00788 vector pyt(3), pyc(3);
00789
00790 vector pdt(3), pdc(3);
00791
00792 vector epsn(ncompstr);
00793 vector peps(3);
00794
00795 matrix epst(3,3), t(3,3);
00796
00797 matrix omegat(3,3), omegac(3,3);
00798
00799 matrix tmp1(3,3), tmp2(3,3);
00800
00801 matrix tomega(3,3);
00802
00803 matrix sigt(3,3);
00804 vector sigma(ncompstr);
00805
00806 double k0, g0, epsv, e, nu;
00807
00808 e = Mm->give_actual_ym(ipp);
00809 nu = Mm->give_actual_nu(ipp);
00810 if (Mm->ip[ipp].ssst == planestress)
00811 Mm->ip[ipp].strain[3] = -nu / (1.0 - nu) * (Mm->ip[ipp].strain[0]+Mm->ip[ipp].strain[1]);
00812
00813 give_actual_param_a(ipp, ido, aa, aat, aac);
00814
00815
00816 for (i=0;i<ncompstr;i++)
00817 epsn[i] = Mm->ip[ipp].strain[i];
00818
00819
00820 dvo = Mm->ip[ipp].eqother[ido+3];
00821
00822 dt[0][0] = Mm->ip[ipp].eqother[ido+3+1+0];
00823 dt[1][1] = Mm->ip[ipp].eqother[ido+3+1+1];
00824 dt[2][2] = Mm->ip[ipp].eqother[ido+3+1+2];
00825 dt[1][2] = dt[2][1] = Mm->ip[ipp].eqother[ido+3+1+3];
00826 dt[0][2] = dt[2][0] = Mm->ip[ipp].eqother[ido+3+1+4];
00827 dt[0][1] = dt[1][0] = Mm->ip[ipp].eqother[ido+3+1+5];
00828
00829 dc[0][0] = Mm->ip[ipp].eqother[ido+3+1+6+0];
00830 dc[1][1] = Mm->ip[ipp].eqother[ido+3+1+6+1];
00831 dc[2][2] = Mm->ip[ipp].eqother[ido+3+1+6+2];
00832 dc[1][2] = dc[2][1] = Mm->ip[ipp].eqother[ido+3+1+6+3];
00833 dc[0][2] = dc[2][0] = Mm->ip[ipp].eqother[ido+3+1+6+4];
00834 dc[0][1] = dc[1][0] = Mm->ip[ipp].eqother[ido+3+1+6+5];
00835
00836
00837 vector_tensor(epsn, epst, Mm->ip[ipp].ssst, strain);
00838 princ_val (epst,peps,t,nijac,limit,Mp->zero,3,1);
00839 yn = damdrvforce_vol(nu, peps);
00840 damdrvforce_dev(nu, peps, pyt, pyc);
00841
00842
00843 glmatrixtransf(dt, t);
00844 glmatrixtransf(dc, t);
00845
00846
00847 dv = dam_vol(ipp, yn, aa, dvo);
00848
00849 for (i=0; i<3; i++)
00850 {
00851 pdt[i] = dt[i][i];
00852 pdc[i] = dc[i][i];
00853 }
00854 pdam_dev(ipp, pyc, pyt, aat, aac, pdt, pdc);
00855
00856
00857 for (i = 0; i < 3; i++)
00858 {
00859 if (dt[i][i] < pdt[i])
00860 dt[i][i] = pdt[i];
00861 if (dc[i][i] < pdc[i])
00862 dc[i][i] = pdc[i];
00863 }
00864
00865 lgmatrixtransf(dt, t);
00866 lgmatrixtransf(dc, t);
00867
00868
00869 if (dv < dvo)
00870 dv = dvo;
00871
00872
00873 princ_val(dt, pdt, tomega, nijac, limit, Mp->zero, 3, 1);
00874 fillm(0.0, omegat);
00875 for(i=0; i<3; i++)
00876 {
00877 if ((pdt[i] < 0.0) || (pdt[i] > 1.0))
00878 {
00879 print_err("principal value of damage tensor for tension is out of range <0,1>.", __FILE__, __LINE__, __func__);
00880 if (pdt[i] < 0.0)
00881 pdt[i] = 0.0;
00882 else
00883 pdt[i] = 1.0;
00884 }
00885 omegat[i][i] = sqrt(1.0-pdt[i]);
00886 }
00887 lgmatrixtransf(omegat, tomega);
00888
00889
00890 princ_val(dc, pdc, tomega, nijac, limit, Mp->zero, 3, 1);
00891 fillm(0.0, omegac);
00892 for(i=0; i<3; i++)
00893 {
00894 if ((pdc[i] < 0.0) || (pdc[i] > 1.0))
00895 {
00896 print_err("principal value of damage tensor for compression is out of range <0,1>.", __FILE__, __LINE__, __func__);
00897 if (pdc[i] < 0.0)
00898 pdc[i] = 0.0;
00899 else
00900 pdc[i] = 1.0;
00901 }
00902 omegac[i][i] = sqrt(1.0-pdc[i]);
00903 }
00904 lgmatrixtransf(omegac, tomega);
00905
00906
00907 k0 = e;
00908 k0 /= (1-2.0*nu);
00909 g0 = e;
00910 g0 /= 2.0*(1.0+nu);
00911
00912 epsv = peps[0]+peps[1]+peps[2];
00913 epsv /= 3.0;
00914
00915 fillm(0.0, eps_t);
00916 for (i = 0; i < 3; i++)
00917 {
00918 if (peps[i] > 0.0)
00919 eps_t[i][i] = peps[i];
00920 }
00921 lgmatrixtransf(eps_t, t);
00922
00923 fillm(0.0, eps_c);
00924 for (i = 0; i < 3; i++)
00925 {
00926 if (peps[i] < 0.0)
00927 eps_c[i][i] = peps[i];
00928 }
00929 lgmatrixtransf(eps_c, t);
00930
00931 sigt[0][0] = k0 - 2.0*g0;
00932 if (epsv > 0.0)
00933 epsv -= dv*epsv;
00934 sigt[0][0] *= epsv;
00935 sigt[1][1] = sigt[2][2] = sigt[0][0];
00936
00937 mxm(omegat, eps_t, tmp1);
00938 mxm(tmp1, omegat, tmp2);
00939 cmulm(2.0*g0, tmp2);
00940 addm(sigt, tmp2, sigt);
00941
00942 mxm(omegac, eps_c, tmp1);
00943 mxm(tmp1, omegac, tmp2);
00944 cmulm(2.0*g0, tmp2);
00945 addm(sigt, tmp2, sigt);
00946 tensor_vector(sigma, sigt, Mm->ip[ipp].ssst, stress);
00947 for (i=0;i<ncompstr;i++)
00948 Mm->ip[ipp].stress[i]=sigma[i];
00949
00950
00951 Mm->ip[ipp].other[ido+3] = dv;
00952
00953 Mm->ip[ipp].other[ido+3+1+0] = dt[0][0];
00954 Mm->ip[ipp].other[ido+3+1+1] = dt[1][1];
00955 Mm->ip[ipp].other[ido+3+1+2] = dt[2][2];
00956 Mm->ip[ipp].other[ido+3+1+3] = dt[1][2];
00957 Mm->ip[ipp].other[ido+3+1+4] = dt[0][2];
00958 Mm->ip[ipp].other[ido+3+1+5] = dt[0][1];
00959
00960 Mm->ip[ipp].other[ido+3+1+6+0] = dc[0][0];
00961 Mm->ip[ipp].other[ido+3+1+6+1] = dc[1][1];
00962 Mm->ip[ipp].other[ido+3+1+6+2] = dc[2][2];
00963 Mm->ip[ipp].other[ido+3+1+6+3] = dc[1][2];
00964 Mm->ip[ipp].other[ido+3+1+6+4] = dc[0][2];
00965 Mm->ip[ipp].other[ido+3+1+6+5] = dc[0][1];
00966
00967 Mm->ip[ipp].other[ido+3+1+6+6+0] = pdt[0];
00968 Mm->ip[ipp].other[ido+3+1+6+6+1] = pdt[1];
00969 Mm->ip[ipp].other[ido+3+1+6+6+2] = pdt[2];
00970
00971 Mm->ip[ipp].other[ido+3+1+6+6+3+0] = pdc[0];
00972 Mm->ip[ipp].other[ido+3+1+6+6+3+1] = pdc[1];
00973 Mm->ip[ipp].other[ido+3+1+6+6+3+2] = pdc[2];
00974 }
00975
00976
00977
00978
00979
00980
00981
00982
00983
00984
00985
00986
00987
00988
00989
00990 void anisodamrot::updateval (long ipp,long im,long ido)
00991 {
00992 long i;
00993 long n = Mm->givencompeqother(ipp,im);
00994
00995 for (i=0;i<n;i++)
00996 Mm->ip[ipp].eqother[ido+i]=Mm->ip[ipp].other[ido+i];
00997 }