00001 #include "ortodamrot.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 // zero limit for principal values of tensors (Jacobi method)
00015 #define qlimit 1.0e-8 // zero limit for switch from quadratuc Bezier to pure quadrat
00016 #ifndef M_PI
00017 #define M_PI 3.14159265358979323846
00018 #endif
00019 #define sqr(x) (x)*(x)
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 ortodamrot::ortodamrot (void)
00030 {
00031 cde = corr_off;
00032 fat = fatigue_off;
00033 ac = 0.0; bc = 0.0; y0c = 0.0; at = 0.0; bt = 0.0; y0t = 0.0;
00034 gfc = 0.0; gft = 0.0;
00035 ft = uft = ult = 0.0;
00036 fc = ufc = ulc = 0.0;
00037 pqt = pqc = 0;
00038 betac = betat = 0.0;
00039 }
00040
00041
00042
00043
00044
00045
00046
00047
00048 ortodamrot::~ortodamrot (void)
00049 {
00050
00051 }
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066 void ortodamrot::read (XFILE *in)
00067 {
00068 xfscanf (in, "%k%m", "dam_evol_f", &dam_evolfunc_kwdset, (int *)&damevf);
00069 switch(damevf)
00070 {
00071 case brittle:
00072 xfscanf (in,"%k%lf %k%lf", "f_t", &ft, "u_ft", &uft);
00073 xfscanf (in,"%k%lf %k%lf", "f_c", &fc, "u_fc", &ufc);
00074 xfscanf (in,"%k%m", "corr_dissip", &corr_disip_en_kwdset, (int *)(&cde));
00075 if (cde != corr_off)
00076 sra.read (in);
00077 break;
00078 case quasi_brittle:
00079 xfscanf (in, "%k%m", "corr_dissip", &corr_disip_en_kwdset, (int *)(&cde));
00080 switch (cde)
00081 {
00082 case corr_off:
00083 xfscanf (in, "%k%lf %k%lf %k%lf", "A_c", &ac, "B_c", &bc, "Y_0c", &y0c);
00084 xfscanf (in, "%k%lf %k%lf %k%lf", "A_t", &at, "B_t", &bt, "Y_0t", &y0t);
00085 break;
00086 case corr_on:
00087 xfscanf (in, "%k%lf %k%lf %k%lf", "G_fc", &gfc, "B_c", &bc, "Y_0c", &y0c);
00088 xfscanf (in, "%k%lf %k%lf %k%lf", "G_ft", &gft, "B_t", &bt, "Y_0t", &y0t);
00089 break;
00090 default:
00091 print_err("unknown type of correction of disipated energy is required.", __FILE__, __LINE__,__func__);
00092 }
00093 break;
00094 case quadbezier:
00095 xfscanf (in,"%k%lf %k%lf %k%lf", "f_t", &ft, "u_ft", &uft, "u_limt", &ult);
00096 if (fabs(uft/ult-0.5) < qlimit)
00097 pqt = 1;
00098 xfscanf (in,"%k%lf %k%lf %k%lf", "f_c", &fc, "u_fc", &ufc, "u_limc", &ulc);
00099 if (fabs(ufc/ulc-0.5) < qlimit)
00100 pqc = 1;
00101 xfscanf (in,"%k%m", "corr_dissip", &corr_disip_en_kwdset, (int *)(&cde));
00102 if (cde != corr_off)
00103 sra.read (in);
00104 break;
00105 default:
00106 print_err("unknown type of damage evolution function is required.", __FILE__, __LINE__,__func__);
00107 }
00108 xfscanf (in,"%k%m", "fatigue", &fatigue_flag_kwdset, (int *)(&fat));
00109 if (fat == fatigue_on)
00110 xfscanf (in,"%k%lf %k%lf", "beta_t", &betat, "beta_c", &betac);
00111 }
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134 void ortodamrot::loadfunc(long ipp, double nu, vector &peps, vector &damt, vector &damc,
00135 double aat, double aac, vector &lft, vector &lfc)
00136 {
00137 long i;
00138 double e, af, auf, aul, kappa, tmp1, tmp2, a, b, h;
00139 vector pyt(3), pyc(3);
00140
00141 switch(damevf)
00142 {
00143 case brittle:
00144 for (i=0; i<3; i++)
00145 {
00146
00147 if (peps(i) > 0.0)
00148 {
00149 kappa = peps(i);
00150 if (kappa <= y0t)
00151 lft(i) = -1.0;
00152 else
00153 {
00154 e = Mm->give_actual_ym(ipp);
00155
00156 af = Mm->give_actual_ft(ipp);
00157
00158 auf = uft*(af/ft);
00159 tmp1 = -(kappa-af/e)/(auf-af/e);
00160 lft(i) = 1.0 - af/(e*kappa)*exp(tmp1) - damt(i);
00161 }
00162 }
00163
00164 if (peps(i) < 0.0)
00165 {
00166 kappa = fabs(peps(i));
00167 if (kappa <= y0c)
00168 lfc(i) = -1.0;
00169 else
00170 {
00171 e = Mm->give_actual_ym(ipp);
00172
00173 af = Mm->give_actual_fc(ipp);
00174
00175 auf = ufc*(af/fc);
00176 tmp1 = -(kappa-af/e)/(auf-af/e);
00177 lfc(i) = 1.0 - af/(e*kappa)*exp(tmp1) - damc(i);
00178 }
00179 }
00180 }
00181 break;
00182 case quasi_brittle:
00183 for (i=0; i<3; i++)
00184 {
00185
00186 if (peps(i) > 0.0)
00187 {
00188 kappa = peps(i);
00189 if (kappa <= y0t)
00190 lft(i) = -1.0;
00191 else
00192 {
00193
00194 tmp1 = kappa - y0t;
00195 tmp2 = log(tmp1);
00196 tmp1 = exp(bt*tmp2);
00197
00198 lft(i) = (1.0 - damt(i))*(1.0+aat*tmp1) - 1.0;
00199 }
00200 }
00201
00202 if (peps(i) < 0.0)
00203 {
00204 kappa = fabs(peps(i));
00205 if (kappa <= y0c)
00206 lfc(i) = -1.0;
00207 else
00208 {
00209
00210 tmp1 = kappa - y0c;
00211 tmp2 = log(tmp1);
00212 tmp1 = exp(bc*tmp2);
00213
00214 lfc(i) = (1.0 - damc(i))*(1.0+aac*tmp1) - 1.0;
00215 }
00216 }
00217 }
00218 break;
00219 case quadbezier:
00220 h = 1.0;
00221 for (i=0; i<3; i++)
00222 {
00223
00224 if (peps(i) > 0.0)
00225 {
00226 kappa = peps(i);
00227 if (kappa <= y0t)
00228 lft(i) = -1.0;
00229 else
00230 {
00231 e = Mm->give_actual_ym(ipp);
00232
00233 af = Mm->give_actual_ft(ipp);
00234
00235 auf = uft*(af/ft);
00236
00237 aul = ult*(af/ft);
00238 a = 4.0*auf*auf;
00239 b = (aul-2.0*auf);
00240 if (pqt)
00241 lft(i) = e*kappa-damt(i)*e*kappa-af*(1.0-damt(i)*kappa*h/auf+sqr(damt(i)*kappa*h)*(aul-auf)/(auf*aul*aul));
00242 else
00243
00244 lft(i) = e*kappa-damt(i)*e*kappa-af*(1.0-auf*(2.0*aul-3.0*auf)/(b*b)-2.0*(aul-auf)/(b*b)*sqrt(a+b*damt(i)*h*kappa)+damt(i)*h*kappa/b);
00245 }
00246 }
00247
00248 if (peps(i) < 0.0)
00249 {
00250 kappa = fabs(peps(i));
00251 if (kappa <= y0c)
00252 lfc(i) = -1.0;
00253 else
00254 {
00255 e = Mm->give_actual_ym(ipp);
00256
00257 af = Mm->give_actual_fc(ipp);
00258
00259 auf = ufc*(af/fc);
00260
00261 aul = ult*(af/ft);
00262 a = 4.0*auf*auf;
00263 b = (aul-2.0*auf);
00264 if (pqc)
00265 lfc(i) = e*kappa-damc(i)*e*kappa-af*(1.0-damc(i)*kappa*h/auf+sqr(damc(i)*kappa*h)*(aul-auf)/(auf*aul*aul));
00266 else
00267
00268 lfc(i) = e*kappa-damc(i)*e*kappa-af*(1.0-auf*(2.0*aul-3.0*auf)/(b*b)-2.0*(aul-auf)/(b*b)*sqrt(a+b*damc(i)*h*kappa)+damc(i)*h*kappa/b);
00269 }
00270 }
00271 }
00272 break;
00273 default:
00274 print_err("unknown type of damage evolution function is required.", __FILE__, __LINE__,__func__);
00275 }
00276 return;
00277 }
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297 double ortodamrot::brittle_damage(long ipp, double y, double e, double f, double uf, double omegao)
00298 {
00299 double err,omega = 0.0, dtmp, tmp, h, indam, kappa;
00300 long i,ni, eid;
00301
00302 ni=sra.give_ni ();
00303 err=sra.give_err ();
00304
00305 indam = f/e;
00306
00307 kappa = y;
00308 if (kappa < indam)
00309
00310 return 0.0;
00311 if ((Mm->ip[ipp].hmt & 2) > 0)
00312 {
00313 omega = 1.0-(f/(e*kappa))*exp(-(kappa-f/e)/(uf-f/e));
00314 return omega;
00315 }
00316
00317 eid = Mm->elip[ipp];
00318 switch (Mt->give_dimension(eid))
00319 {
00320 case 1:
00321 h = Mt->give_length(eid);
00322 break;
00323 case 2:
00324 h = sqrt(Mt->give_area(eid));
00325 break;
00326 case 3:
00327 h = pow(Mt->give_volume(eid), 1.0/3.0);
00328 break;
00329 default:
00330 print_err("unknown dimension of element is required.", __FILE__, __LINE__, __func__);
00331 }
00332 if (cde == corr_off)
00333
00334 {
00335 omega = 1.0-(f/(e*kappa))*exp(-(kappa-f/e)/(uf-f/e));
00336 return omega;
00337 }
00338
00339 omega = 1.0;
00340 for (i = 0; i < ni; i++)
00341 {
00342 dtmp = -e*kappa+f/uf*h*kappa*exp(-h*omega*kappa/uf);
00343 tmp = (1.0-omega)*e*kappa-f*exp(-h*omega*kappa/uf);
00344 if (fabs(tmp) < ft*err)
00345 {
00346 omega += - tmp/dtmp;
00347 break;
00348 }
00349 omega += - tmp/dtmp;
00350 }
00351 if ((omega > 1.0) || (omega < 0.0) || (fabs(tmp) > ft*err))
00352
00353
00354 {
00355 omega = omegao;
00356 for (i = 0; i < ni; i++)
00357 {
00358 dtmp = -e*kappa+ft/uf*h*kappa*exp(-h*omega*kappa/uf);
00359 tmp = (1.0-omega)*e*kappa-ft*exp(-h*omega*kappa/uf);
00360 if (fabs(tmp) < ft*err)
00361 {
00362 omega += - tmp/dtmp;
00363 break;
00364 }
00365 omega += - tmp/dtmp;
00366 }
00367 }
00368
00369 if (fabs(tmp) > ft*err)
00370 {
00371 print_err("iteration of omega doesn't converge with required error\n"
00372 " elem %ld, ip %ld:\n"
00373 " kappa=%le, ft=%le, e=%le, uf=%le, h=%le",
00374 __FILE__, __LINE__, __func__, Mm->elip[ipp]+1, ipp+1, kappa, ft, e, uf, h);
00375 abort();
00376 }
00377
00378 if ((omega > 1.0) || (omega < 0.0))
00379 {
00380 print_err("computed omega=%le is out of prescribed range <0,1>\n"
00381 " elem %ld, ip %ld:\n"
00382 " kappa=%le, ft=%le, e=%le, uf=%le, h=%le",
00383 __FILE__, __LINE__, __func__, Mm->elip[ipp]+1, ipp+1, kappa, ft, e, uf, h);
00384 abort();
00385 }
00386
00387 return omega;
00388 }
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410 double ortodamrot::qbezier_damage(long ipp, double y, double e, double f, double uf, double ul,double omegao,long pq)
00411 {
00412 double err,omega = 0.0, dtmp, tmp, h, indam, kappa, a, b;
00413 long i,ni, eid;
00414
00415 ni=sra.give_ni ();
00416 err=sra.give_err ();
00417
00418 indam = f/e;
00419
00420 kappa = y;
00421 if (kappa-indam <= Mp->zero)
00422
00423 return 0.0;
00424 a = uf*uf;
00425 b = (ul-2.0*uf);
00426 if ((Mm->ip[ipp].hmt & 2) > 0)
00427 {
00428 if (kappa >= ul)
00429 return 1.0;
00430 omega = 1.0-2.0*uf*(ul-uf)/(b*b)-2.0*(ul-uf)*sqrt(a+b*(kappa-f/e))/(b*b)+(kappa-f/e)/b;
00431 omega = 1.0-(f/(e*(kappa-f/e)))*omega;
00432 return omega;
00433 }
00434
00435 eid = Mm->elip[ipp];
00436 switch (Mt->give_dimension(eid))
00437 {
00438 case 1:
00439 h = Mt->give_length(eid);
00440 break;
00441 case 2:
00442 h = sqrt(Mt->give_area(eid));
00443 break;
00444 case 3:
00445 h = pow(Mt->give_volume(eid), 1.0/3.0);
00446 break;
00447 default:
00448 print_err("unknown dimension of element is required.", __FILE__, __LINE__, __func__);
00449 }
00450 if (cde == corr_off)
00451
00452 {
00453 if (kappa >= ul)
00454 return 1.0;
00455 if (pq)
00456 omega = (1.0-omega)*e*kappa-f*(1.0-(kappa-f/e)/uf+(ul-uf)/(uf*ul*ul)*(kappa-f/e)*(kappa-f/e));
00457 else
00458 {
00459 omega = 1.0-2.0*uf*(ul-uf)/(b*b)-2.0*(ul-uf)*sqrt(a+b*(kappa-f/e))/(b*b)+(kappa-f/e)/b;
00460 omega = 1.0-(f/(e*(kappa-f/e)))*omega;
00461 }
00462 return omega;
00463 }
00464 if (omega < 0.5)
00465 omega = 0.5;
00466 else
00467 omega = omegao;
00468 if (omega < 0.0)
00469 print_err("previous iteration of omega returned value < 0.", __FILE__, __LINE__, __func__);
00470
00471 if (kappa*omegao*h >= ul)
00472 return 1.0;
00473 for (i = 0; i < ni; i++)
00474 {
00475 if (pq)
00476 {
00477 dtmp = 2.0*f*(uf-ul)*sqr(kappa*h)*omega/(uf*ul*ul)+kappa*h*f/uf-e*kappa;
00478 tmp = (1.0-omega)*e*kappa-f*(1.0-omega*kappa*h/uf+(ul-uf)*sqr(omega*kappa*h)/(uf*ul*ul));
00479 }
00480 else
00481 {
00482 dtmp = f*(ul-uf)*kappa*h/(b*sqrt(a+omega*b*kappa*h))-(e+f*h/b)*kappa;
00483 tmp = (1.0-omega)*e*kappa-f*(1.0+(2.0*uf*(ul-uf))/(b*b)-(2.0*(ul-uf)*sqrt(a+b*omega*h*kappa))/(b*b)+omega*h*kappa/b);
00484 }
00485 omega += - tmp/dtmp;
00486 if (omega < 0.0)
00487 {
00488 print_err("during iteration, omega reaches value < 0.", __FILE__, __LINE__, __func__);
00489 abort();
00490 }
00491 if (fabs(tmp) <= ft*err)
00492 break;
00493 }
00494 if (fabs(tmp) > ft*err)
00495 print_err("iteration of omega doesn't converge with required error.", __FILE__, __LINE__, __func__);
00496
00497 if (omega > 1.0)
00498 {
00499 if ((omega-1.0) > err)
00500 print_err("iteration of omega returns value > 1.", __FILE__, __LINE__, __func__);
00501 else
00502 omega = 1.0;
00503 }
00504 if (omega < 0.0)
00505 {
00506 if (fabs(omega) > err)
00507 print_err("iteration of omega returns value < 0.", __FILE__, __LINE__, __func__);
00508 else
00509 omega = 0.0;
00510 }
00511
00512 return omega;
00513 }
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526
00527
00528
00529
00530
00531
00532 void ortodamrot::princ_dam(long ipp, vector &peps, double aat, double aac, vector &pdamt, vector &pdamc)
00533 {
00534 long i;
00535 double tmp, ltmp;
00536 double e, af, auf, aul;
00537
00538
00539 switch(damevf)
00540 {
00541 case brittle:
00542 e = Mm->give_actual_ym(ipp);
00543 for(i=0; i<3; i++)
00544 {
00545
00546 if (peps(i) > 0.0)
00547 {
00548
00549 af = Mm->give_actual_ft(ipp);
00550
00551 auf = uft*(af/ft);
00552 pdamt(i) = brittle_damage(ipp, peps(i), e, af, auf, pdamt(i));
00553 }
00554
00555 if (peps(i) < 0.0)
00556 {
00557
00558 af = Mm->give_actual_fc(ipp);
00559
00560 auf = ufc*(af/fc);
00561 pdamc(i) = brittle_damage(ipp, fabs(peps(i)), e, af, auf, pdamc(i));
00562 }
00563 }
00564 break;
00565 case quasi_brittle:
00566 for(i=0; i<3; i++)
00567 {
00568
00569 if (peps(i) > 0.0)
00570 {
00571 if (peps(i) > y0t)
00572 {
00573 ltmp = log(peps(i) - y0t);
00574 pdamt(i) = aat*exp(ltmp*bt);
00575 tmp = 1.0/(1.0+pdamt(i));
00576 pdamt(i) = 1.0 - tmp;
00577 }
00578 if (pdamt(i) < 0.0)
00579 {
00580 fprintf(stdout, "\nDamage < 0.0 (ipp=%ld)\n", ipp);
00581 abort();
00582 }
00583 }
00584
00585 if (peps(i) < 0.0)
00586 {
00587 if (fabs(peps(i)) > y0c)
00588 {
00589 ltmp = log(fabs(peps(i)) - y0c);
00590 pdamc(i) = aac*exp(ltmp*bc);
00591 tmp = 1.0/(1.0+pdamc(i));
00592 pdamc(i) = 1.0-tmp;
00593 }
00594 if (pdamc(i) < 0.0)
00595 {
00596 fprintf(stdout, "\nDamage < 0.0 (ipp=%ld)\n", ipp);
00597 abort();
00598 }
00599 }
00600 }
00601 break;
00602 case quadbezier:
00603 e = Mm->give_actual_ym(ipp);
00604 for(i=0; i<3; i++)
00605 {
00606
00607 if (peps(i) > 0.0)
00608 {
00609
00610 af = Mm->give_actual_ft(ipp);
00611
00612 auf = uft*(af/ft);
00613
00614 aul = ult*(af/ft);
00615 pdamt(i) = qbezier_damage(ipp, peps(i), e, af, auf, aul, pdamt(i),pqt);
00616 }
00617
00618 if (peps(i) < 0.0)
00619 {
00620
00621 af = Mm->give_actual_fc(ipp);
00622
00623 auf = ufc*(af/fc);
00624
00625 aul = ulc*(af/fc);
00626 pdamc(i) = qbezier_damage(ipp, fabs(peps(i)), e, af, auf, aul, pdamc(i),pqc);
00627 }
00628 }
00629 break;
00630 default:
00631 print_err("unknown type of damage evolution function is required.", __FILE__, __LINE__,__func__);
00632 }
00633 }
00634
00635
00636
00637
00638
00639
00640
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653
00654 void ortodamrot::give_actual_param_a(long ipp, long ido, double &aat, double &aac)
00655 {
00656 long eid;
00657 double e, h, tmp;
00658
00659 aat = Mm->ip[ipp].eqother[ido];
00660 aac = Mm->ip[ipp].eqother[ido+1];
00661 if ((aac == 0.0) || (aat == 0.0))
00662 {
00663 if (cde == corr_off)
00664
00665 {
00666 aat = at;
00667 aac = ac;
00668 Mm->ip[ipp].eqother[ido] = aat;
00669 Mm->ip[ipp].eqother[ido+1] = aac;
00670 return;
00671 }
00672 eid = Mm->elip[ipp];
00673 switch (Mt->give_dimension(eid))
00674 {
00675 case 1:
00676 h = Mt->give_length(eid);
00677 break;
00678 case 2:
00679 h = sqrt(Mt->give_area(eid));
00680 break;
00681 case 3:
00682 h = pow(Mt->give_volume(eid), 1.0/3.0);
00683 break;
00684 default:
00685 print_err("cannot determine dimension of element.", __FILE__, __LINE__, __func__);
00686 }
00687 e = Mm->give_actual_ym(ipp);
00688
00689 tmp = (gft/h - y0t/e)*bt*e*sin(M_PI/bt)/M_PI;
00690
00691 if (tmp < 0.0)
00692 print_err("cannot express variable softening for tensile damage.", __FILE__, __LINE__, __func__);
00693
00694 aat = pow(tmp, -bt);
00695
00696 tmp = (gfc/h - y0c/e)*bc*e*sin(M_PI/bc)/M_PI;
00697
00698 if (tmp < 0.0)
00699 print_err("cannot express variable softening for compressive damage.", __FILE__, __LINE__, __func__);
00700
00701 aac = pow(tmp, -bc);
00702 Mm->ip[ipp].eqother[ido] = aat;
00703 Mm->ip[ipp].eqother[ido+1] = aac;
00704 }
00705 }
00706
00707
00708
00709
00710
00711
00712
00713
00714
00715
00716
00717
00718
00719 double ortodamrot::give_actual_ft(long ipp)
00720 {
00721 return ft;
00722 }
00723
00724
00725
00726
00727
00728
00729
00730
00731
00732
00733
00734
00735
00736 double ortodamrot::give_actual_fc(long ipp)
00737 {
00738 return fc;
00739 }
00740
00741
00742
00743
00744
00745
00746
00747
00748
00749
00750
00751
00752
00753
00754 void ortodamrot::initvalues(long ipp, long ido)
00755 {
00756 double aat, aac;
00757 give_actual_param_a(ipp, ido, aac, aat);
00758 }
00759
00760
00761
00762
00763
00764
00765
00766
00767
00768
00769
00770
00771
00772
00773 void ortodamrot::matstiff (matrix &d,long ipp,long ido)
00774 {
00775
00776
00777 switch (Mp->nlman->stmat)
00778 {
00779 case initial_stiff:
00780 elmatstiff(d, ipp);
00781 break;
00782 case tangent_stiff:
00783 elmatstiff(d, ipp);
00784
00785 break;
00786 default:
00787 print_err("unknown type of stifness matrix is required.", __FILE__, __LINE__, __func__);
00788 }
00789 }
00790
00791
00792
00793
00794
00795
00796
00797
00798
00799
00800
00801
00802
00803 void ortodamrot::elmatstiff (matrix &d,long ipp)
00804 {
00805 double e, e0;
00806 long idem = Mm->ip[ipp].gemid();
00807 long idoem=Mm->givencompeqother(ipp,0);
00808 long tmp=Mm->givencompeqother(ipp,idem);
00809 idoem -= tmp;
00810
00811 Mm->elmatstiff (d,ipp);
00812 e = Mm->give_actual_ym(ipp);
00813 e0 = Mm->give_actual_ym(ipp,idem,idoem);
00814 cmulm(e/e0, d);
00815 }
00816
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826
00827
00828
00829
00830
00831
00832 void ortodamrot::pelmatstiff (long ipp, matrix &d, double e, double nu)
00833 {
00834 double b;
00835
00836 switch (Mm->ip[ipp].ssst)
00837 {
00838 case planestress:
00839 b = e/(1.0-nu*nu);
00840 d[0][0] = b; d[0][1] = b*nu; d[0][2] = b*nu;
00841 d[1][0] = b*nu; d[1][1] = b; d[1][2] = b*nu;
00842 d[2][0] = b*nu; d[2][1] = b*nu; d[2][2] = b;
00843 break;
00844 case planestrain:
00845 case spacestress:
00846 b = e /((1.0 + nu)*(1.0 - 2.0*nu));
00847 d[0][0] = d[1][1] = d[2][2] = 1.0 - nu;
00848 d[0][1] = d[0][2] = d[1][0] = d[1][2] = d[2][0] = d[2][1] = nu;
00849 cmulm(b, d);
00850 break;
00851 default:
00852 print_err("unknown type of stress/strain state is required.", __FILE__, __LINE__, __func__);
00853 break;
00854 }
00855 }
00856
00857
00858
00859
00860
00861
00862
00863
00864
00865
00866
00867
00868
00869
00870
00871 void ortodamrot::nlstresses (long ipp, long im, long ido)
00872 {
00873 long i;
00874 long ncompstr=Mm->ip[ipp].ncompstr;
00875
00876
00877 double aat, aac;
00878
00879 matrix dt(3,3), dc(3,3);
00880
00881 matrix d(ncompstr, ncompstr);
00882
00883 matrix eps_t(3,3), eps_c(3,3);
00884
00885 vector pyt(3), pyc(3);
00886
00887 vector pdt(3), pdc(3);
00888
00889 vector epsn(ncompstr);
00890 vector peps(3);
00891
00892 vector eps(ncompstr);
00893
00894 matrix epst(3,3), t(3,3);
00895
00896 vector epsp(ncompstr);
00897 matrix epspt(3,3);
00898
00899 matrix omegat(3,3), omegac(3,3);
00900
00901 matrix tmp(3,3);
00902
00903 matrix tomega(3,3);
00904
00905 matrix sigt(3,3);
00906
00907 vector sig(ncompstr);
00908
00909 vector sigma(ncompstr);
00910
00911 vector psig(3);
00912
00913 double e, nu, h;
00914
00915 long eid;
00916
00917 if ((Mp->matmodel == nonlocal) && (Mm->ip[ipp].hmt & 2))
00918
00919 {
00920 if (Mp->nonlocphase == 1)
00921 {
00922
00923
00924 return;
00925 }
00926 else
00927 {
00928
00929 for (i=0;i<ncompstr;i++)
00930 epsn[i] = Mm->ip[ipp].nonloc[i];
00931 }
00932 }
00933 else
00934
00935 {
00936
00937 for (i=0;i<ncompstr;i++)
00938 epsn[i] = Mm->ip[ipp].strain[i];
00939 }
00940
00941
00942 give_actual_param_a(ipp, ido, aat, aac);
00943 e = Mm->give_actual_ym(ipp);
00944 nu = Mm->give_actual_nu(ipp);
00945 if (Mm->ip[ipp].ssst == planestress)
00946 epsn[3] = -nu/(1.0-nu)*(epsn[0]+epsn[1]);
00947
00948
00949
00950
00951 dt[0][0] = Mm->ip[ipp].eqother[ido+2+0];
00952 dt[1][1] = Mm->ip[ipp].eqother[ido+2+1];
00953 dt[2][2] = Mm->ip[ipp].eqother[ido+2+2];
00954 dt[1][2] = dt[2][1] = Mm->ip[ipp].eqother[ido+2+3];
00955 dt[0][2] = dt[2][0] = Mm->ip[ipp].eqother[ido+2+4];
00956 dt[0][1] = dt[1][0] = Mm->ip[ipp].eqother[ido+2+5];
00957
00958 dc[0][0] = Mm->ip[ipp].eqother[ido+2+6+0];
00959 dc[1][1] = Mm->ip[ipp].eqother[ido+2+6+1];
00960 dc[2][2] = Mm->ip[ipp].eqother[ido+2+6+2];
00961 dc[1][2] = dc[2][1] = Mm->ip[ipp].eqother[ido+2+6+3];
00962 dc[0][2] = dc[2][0] = Mm->ip[ipp].eqother[ido+2+6+4];
00963 dc[0][1] = dc[1][0] = Mm->ip[ipp].eqother[ido+2+6+5];
00964
00965
00966 if (fat == fatigue_on)
00967 {
00968 epspt[0][0] = Mm->ip[ipp].eqother[ido+2+6+6+0];
00969 epspt[1][1] = Mm->ip[ipp].eqother[ido+2+6+6+1];
00970 epspt[2][2] = Mm->ip[ipp].eqother[ido+2+6+6+2];
00971 epspt[1][2] = epspt[2][1] = Mm->ip[ipp].eqother[ido+2+6+6+3];
00972 epspt[0][2] = epspt[2][0] = Mm->ip[ipp].eqother[ido+2+6+6+4];
00973 epspt[0][1] = epspt[1][0] = Mm->ip[ipp].eqother[ido+2+6+6+5];
00974 tensor_vector(epsp, epspt, Mm->ip[ipp].ssst, strain);
00975 subv(epsn,epsp,eps);
00976 }
00977 else
00978 copyv(epsn, eps);
00979
00980
00981 vector_tensor(eps, epst, Mm->ip[ipp].ssst, strain);
00982 princ_val(epst, peps, t, nijac, limit, Mp->zero, 3, 1);
00983
00984 glmatrixtransf(dt, t);
00985 glmatrixtransf(dc, t);
00986
00987 for (i=0; i<3; i++)
00988 {
00989 pdt[i] = dt[i][i];
00990 pdc[i] = dc[i][i];
00991 }
00992
00993
00994 princ_dam(ipp, peps, aat, aac, pdt, pdc);
00995 for (i=0; i<3; i++)
00996 {
00997 if (pdt[i] > dt[i][i])
00998 dt[i][i] = pdt[i];
00999 if (pdc[i] > dc[i][i])
01000 dc[i][i] = pdc[i];
01001 }
01002 lgmatrixtransf(dt, t);
01003 lgmatrixtransf(dc, t);
01004
01005 if ((Mp->matmodel == nonlocal) && (Mm->ip[ipp].hmt & 2))
01006
01007 {
01008
01009 for (i=0;i<ncompstr;i++)
01010 epsn[i] = Mm->ip[ipp].strain[i];
01011 }
01012
01013
01014 princ_val(dt, pdt, tomega, nijac, limit, Mp->zero, 3, 1);
01015 fillm(0.0, omegat);
01016 for(i=0; i<3; i++)
01017 {
01018 if ((pdt[i] < -1.0e-8) || (pdt[i] > 1.00000001))
01019 {
01020
01021 if (pdt[i] < 0.0)
01022 pdt[i] = 0.0;
01023 else
01024 pdt[i] = 1.0;
01025 }
01026 omegat[i][i] = sqrt(1.0-pdt[i]);
01027 }
01028 lgmatrixtransf(omegat, tomega);
01029
01030 if (fat == fatigue_on)
01031 {
01032 fillm(0.0, epspt);
01033 for (i=0; i<3; i++)
01034 epspt[i][i] = 0.5*betat*pdt[i]*pdt[i]/(1.0-pdt[i]);
01035 lgmatrixtransf(epspt, tomega);
01036 tensor_vector(epsp, epspt, Mm->ip[ipp].ssst, strain);
01037 subv(epsn,epsp,eps);
01038 }
01039 else
01040 copyv(epsn, eps);
01041
01042
01043 princ_val(dc, pdc, tomega, nijac, limit, Mp->zero, 3, 1);
01044 fillm(0.0, omegac);
01045 for(i=0; i<3; i++)
01046 {
01047 if ((pdc[i] < -1.0e-8) || (pdc[i] > 1.00000001))
01048 {
01049
01050 if (pdc[i] < 0.0)
01051 pdc[i] = 0.0;
01052 else
01053 pdc[i] = 1.0;
01054 }
01055 omegac[i][i] = sqrt(1.0-pdc[i]);
01056 }
01057 lgmatrixtransf(omegac, tomega);
01058
01059 if (fat == fatigue_on)
01060 {
01061 fillm(0.0, epspt);
01062 for (i=0; i<3; i++)
01063 epspt[i][i] -= 0.5*betac*pdc[i]*pdc[i]/(1.0-pdc[i]);
01064 lgmatrixtransf(epspt, tomega);
01065 tensor_vector(epsp, epspt, Mm->ip[ipp].ssst, strain);
01066 subv(epsn,epsp,eps);
01067 }
01068 else
01069 copyv(epsn, eps);
01070
01071
01072
01073 vector_tensor(eps, epst, Mm->ip[ipp].ssst, strain);
01074 princ_val(epst, peps, t, nijac, limit, Mp->zero, 3, 1);
01075
01076
01077 fillm(0.0, epst);
01078 for(i=0; i<3; i++)
01079 {
01080 if (peps[i] > 0.0)
01081 epst[i][i] = peps[i];
01082 }
01083 lgmatrixtransf(epst, t);
01084 tensor_vector(eps, epst, Mm->ip[ipp].ssst, strain);
01085
01086
01087 Mm->elmatstiff(d, ipp);
01088 mxv(d, eps, sig);
01089 vector_tensor(sig, sigt, Mm->ip[ipp].ssst, stress);
01090
01091 mxm(omegat, sigt, tmp);
01092 mxm(tmp, omegat, sigt);
01093 tensor_vector(sigma, sigt, Mm->ip[ipp].ssst, stress);
01094
01095
01096 fillm(0.0, epst);
01097 for(i=0; i<3; i++)
01098 {
01099 if (peps[i] < 0.0)
01100 epst[i][i] = peps[i];
01101 }
01102 lgmatrixtransf(epst, t);
01103 tensor_vector(eps, epst, Mm->ip[ipp].ssst, strain);
01104
01105
01106 Mm->elmatstiff(d, ipp);
01107 mxv(d, eps, sig);
01108 vector_tensor(sig, sigt, Mm->ip[ipp].ssst, stress);
01109
01110 mxm(omegat, sigt, tmp);
01111 mxm(tmp, omegat, sigt);
01112 tensor_vector(sig, sigt, Mm->ip[ipp].ssst, stress);
01113
01114
01115 addv(sigma, sig, sigma);
01116
01117
01118 for (i=0;i<ncompstr;i++)
01119 Mm->ip[ipp].stress[i]=sigma[i];
01120
01121 Mm->ip[ipp].other[ido+2+0] = dt[0][0];
01122 Mm->ip[ipp].other[ido+2+1] = dt[1][1];
01123 Mm->ip[ipp].other[ido+2+2] = dt[2][2];
01124 Mm->ip[ipp].other[ido+2+3] = dt[1][2];
01125 Mm->ip[ipp].other[ido+2+4] = dt[0][2];
01126 Mm->ip[ipp].other[ido+2+5] = dt[0][1];
01127
01128 Mm->ip[ipp].other[ido+2+6+0] = dc[0][0];
01129 Mm->ip[ipp].other[ido+2+6+1] = dc[1][1];
01130 Mm->ip[ipp].other[ido+2+6+2] = dc[2][2];
01131 Mm->ip[ipp].other[ido+2+6+3] = dc[1][2];
01132 Mm->ip[ipp].other[ido+2+6+4] = dc[0][2];
01133 Mm->ip[ipp].other[ido+2+6+5] = dc[0][1];
01134
01135
01136 if (fat == fatigue_on)
01137 {
01138 subv(epsn, eps, epsp);
01139 Mm->ip[ipp].other[ido+2+6+6+0] = epsp[0];
01140 Mm->ip[ipp].other[ido+2+6+6+1] = epsp[1];
01141 Mm->ip[ipp].other[ido+2+6+6+2] = epsp[2];
01142 Mm->ip[ipp].other[ido+2+6+6+3] = epsp[3];
01143 Mm->ip[ipp].other[ido+2+6+6+4] = epsp[4];
01144 Mm->ip[ipp].other[ido+2+6+6+5] = epsp[5];
01145 }
01146
01147
01148 eid = Mm->elip[ipp];
01149 switch (Mt->give_dimension(eid))
01150 {
01151 case 1:
01152 h = Mt->give_length(eid);
01153 break;
01154 case 2:
01155 h = sqrt(Mt->give_area(eid));
01156 break;
01157 case 3:
01158 h = pow(Mt->give_volume(eid), 1.0/3.0);
01159 break;
01160 default:
01161 print_err("unknown dimension of element", __FILE__, __LINE__, __func__);
01162 }
01163 vector_tensor(sigma, sigt, Mm->ip[ipp].ssst, stress);
01164 princ_val(sigt, psig, t, nijac, limit, Mp->zero, 3, 1);
01165 for (i = 0; i < 3; i++)
01166 {
01167 if (peps[i] > 0.0)
01168 Mm->ip[ipp].other[ido+2+6+6+6+i] = pdt[i]*(peps[i]-psig[i]/e)*h;
01169 else
01170 Mm->ip[ipp].other[ido+2+6+6+6+i] = -pdc[i]*(peps[i]-psig[i]/e)*h;
01171 }
01172 }
01173
01174
01175
01176
01177
01178
01179
01180
01181
01182
01183
01184
01185
01186
01187
01188 void ortodamrot::updateval (long ipp,long im,long ido)
01189 {
01190 long i;
01191 long n = Mm->givencompeqother(ipp,im);
01192
01193 for (i=0;i<n;i++)
01194 Mm->ip[ipp].eqother[ido+i]=Mm->ip[ipp].other[ido+i];
01195 }