00001 #include <math.h>
00002 #include <float.h>
00003 #include <stdlib.h>
00004 #include "shefplast.h"
00005 #include "global.h"
00006 #include "vecttens.h"
00007 #include "intpoints.h"
00008 #include "matrix.h"
00009 #include "vector.h"
00010 #include "tensor.h"
00011
00012 #define nijac 20
00013 #define limit 1.0e-8
00014
00015
00016
00017
00018 shefplast::shefplast (void)
00019 {
00020 rc = 0.0;
00021 rt = 0.0;
00022 gamma = 0.99;
00023 alpha = 0.5;
00024 p = 0.0;
00025 a = -1.0;
00026 k0 = 0.1;
00027 ah = bh = ch = 0.0;
00028 kh0 = 0.001;
00029 allocm(3,3,dev);
00030 khat = rhoc = r = m = xi = ft = 0.0;
00031 b0 = b1 = c0 = c1 = c2 = c3 = c4 = c5 = theta = d0 = d1 = d2 = 0.0;
00032 }
00033
00034
00035
00036
00037 shefplast::~shefplast (void)
00038 {
00039
00040 }
00041
00042
00043
00044
00045
00046
00047
00048
00049 void shefplast::read (XFILE *in)
00050 {
00051 xfscanf (in, "%lf %lf %lf %lf %lf %lf", &rc, &rt, &p, &ah, &bh, &ch);
00052 sra.read (in);
00053
00054 ft = rt/rc;
00055 m = 3.0*(1.0-pow(ft,2.0/gamma));
00056 m /= ft+2.0*pow(ft,1.0/gamma);
00057 }
00058
00059 void shefplast::compute_khat(matrix &sig, vector &q)
00060 {
00061 double k;
00062
00063 k = k0 + (1.0 - k0)*sqrt(q[0]*(2.0-q[0]));
00064 khat = pow(k,2.0*p)*(1.0 - xi*xi*(1.0-k)*(1.0-k)/(a*a));
00065 }
00066
00067 void shefplast::compute_rhoc(matrix &sig)
00068 {
00069 double tmp;
00070
00071 tmp = -m+sqrt(m*m-12.0*sqrt(3.0)*m*xi+36.0);
00072 rhoc = pow(1.0/6.0, gamma)*sqrt(2.0/3.0)*pow(tmp,gamma);
00073 }
00074
00075 void shefplast::compute_r(matrix &sig)
00076 {
00077 double tmp, rhoe;
00078 double pi=4.0*atan(1.0) ;
00079
00080 tmp = -m+sqrt(m*m-3.0*sqrt(3.0)*m*xi+9.0);
00081 rhoe = pow(1.0/3.0, gamma)*sqrt(2.0/3.0)*pow(tmp, gamma);
00082 b0 = rhoe/rhoc;
00083 if(b0<=0.5) {
00084 b0 = 0.5 + 1.0e-6 ;
00085 }
00086 else if(b0>1.0) {
00087 b0 = 1.0 ;
00088 }
00089 b1 = sqrt(3.0)*(1.0-alpha)*b0/(1.0+b0)+2.0*alpha*b0/sqrt(3.0);
00090 c0 = (2.0-sqrt(3.0)*b1)*(2.0*b0-sqrt(3.0)*b1) / ((b1*(1.0+b0)-sqrt(3.0)*b0)*(b1*(1.0+b0)-sqrt(3.0)*b0));
00091 c1 = 3.0-c0*(1.0+b0)*(1.0+b0);
00092 c2 = 1.0+3.0*c0*(1.0-b0)*(1.0-b0);
00093 c3 = 2.0*c0*sqrt(3.0)*(1.0-b0*b0);
00094 c4 = (1.0+b0)*(1.0-b0*c0);
00095 c5 = (1.0-b0)*(1.0-3.0*b0*c0);
00096 if((-sqrt(3.0)/2.0*3.0*j3s/pow(j2s,1.5)) > 1.0) {
00097 theta = pi/6.0 ;
00098 }
00099 else if((-sqrt(3.0)/2.0*3.0*j3s/pow(j2s,1.5)) < -1.0) {
00100 theta = -pi/6.0 ;
00101 }
00102 else
00103 theta = 1.0/3.0*asin(-sqrt(3.0)/2.0*3.0*j3s/pow(j2s,1.5));
00104 d0 = c1*cos(theta)*cos(theta) - c2*sin(theta)*sin(theta) + c3*sin(theta)*cos(theta);
00105 d1 = 2.0*(c4*sqrt(3.0)*cos(theta)-c5*sin(theta));
00106 d2 = b0*(4.0-3.0*b0*c0);
00107 r = 2.0*d0/(d1-sqrt(d1*d1-4.0*d0*d2));
00108 }
00109
00110 void shefplast::compute_classvar(matrix &sig, vector &q)
00111 {
00112 deviator(sig,dev);
00113 i1s = first_invar(sig);
00114 j2s = second_invar(dev);
00115 j3s = third_invar(dev);
00116 xi = i1s/(rc*sqrt(3.0));
00117 compute_khat(sig, q);
00118 compute_rhoc(sig);
00119 compute_r(sig);
00120 }
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132 double shefplast::yieldfunction (matrix &sig, vector &q)
00133 {
00134 double rho, f;
00135
00136 compute_classvar(sig, q);
00137 rho = pow(2.0*j2s,0.5)/rc;
00138 f = rho*rho - khat*rhoc*rhoc /(r*r);
00139
00140 return f;
00141 }
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154 void shefplast::derqdsigma (matrix &sig, matrix &drds)
00155 {
00156 long i, j, k;
00157 double dqdd0, dqdd1, dqdd2, dd0db0, dd1db0, dd2db0;
00158 double db1db0, dc0db0,dc1db0,dc2db0,dc3db0,dc4db0,dc5db0;
00159 double dqdb0, dqdth, db0dxi;
00160 double dd0dth, dd1dth;
00161 double dthdj2, dthdj3;
00162 double tmp = sqrt(d1*d1-4.0*d0*d2);
00163 double tmp1, tmp2, tmp3;
00164 double drhocdxi,drhoedxi,rhoe;
00165 double pi=4.0*atan(1.0) ;
00166 double epsi=1.0e-6 ;
00167 matrix dthds(3,3);
00168 matrix dj3ds(3,3);
00169 matrix dxids(3,3);
00170 matrix db0ds(3,3);
00171
00172 dqdd0 = 4*d2*d0/2.0/tmp - d1+tmp;
00173 dqdd0 /= (2.0*d0*d0);
00174 dqdd1 = (1.0-d1/tmp)/2.0/d0;
00175 dqdd2 = 1.0/tmp;
00176 db1db0 = sqrt(3.0)*(1.0-alpha)/(1.0+b0)/(1.0+b0)+2.0*alpha/sqrt(3.0);
00177
00178 tmp1 = -sqrt(3.0)*db1db0*(2.0*b0-sqrt(3.0)*b1)+(2.0-sqrt(3.0)*b1)*(2.0-sqrt(3.0)*db1db0);
00179 tmp2 = b1*(1.0+b0)-sqrt(3.0)*b0;
00180 tmp3 = 2.0*tmp2*(db1db0*(1.0+b0)+b1-sqrt(3.0))*(2.0-sqrt(3.0)*b1)*(2.0*b0-sqrt(3.0)*b1);
00181
00182 dc0db0 = (tmp1*tmp2*tmp2 - tmp3)/pow(tmp2, 4.0);
00183 dc1db0 = -dc0db0*(1.0+b0)*(1.0+b0) - c0*2.0*(1.0+b0);
00184 dc2db0 = 3.0*dc0db0*(1.0-b0)*(1.0-b0)-6.0*c0*(1.0-b0);
00185 dc3db0 = 2.0*sqrt(3.0)*dc0db0*(1.0-b0*b0)-4.0*sqrt(3.0)*c0*b0;
00186 dc4db0 = 1.0-b0*c0-(1.0+b0)*(c0+b0*dc0db0);
00187 dc5db0 = 3.0*b0*c0-1.0-3.0*(1.0-b0)*(c0+b0*dc0db0);
00188
00189 dd0db0 = dc1db0*cos(theta)*cos(theta)-dc2db0*sin(theta)*sin(theta)+dc3db0*sin(theta)*cos(theta);
00190 dd1db0 = 2.0*(dc4db0*sqrt(3.0)*cos(theta)-dc5db0*sin(theta));
00191 dd2db0 = 4.0 - 3.0*b0*c0 - 3.0*b0*(c0+b0*dc0db0);
00192
00193 dqdb0 = dqdd0*dd0db0 + dqdd1*dd1db0 + dqdd2*dd2db0;
00194
00195
00196 tmp1 = sqrt(m*m-12.0*sqrt(3.0)*m*xi+36.0);
00197 drhocdxi = -gamma*pow((tmp1-m)/6.0, gamma-1.0);
00198 drhocdxi *= m*sqrt(2.0)/tmp1;
00199
00200 tmp1 = sqrt(m*m-3.0*sqrt(3.0)*m*xi+9.0);
00201 drhoedxi = -gamma*pow((tmp1-m)/3.0, gamma-1.0);
00202 drhoedxi *= m/sqrt(2.0)/tmp1;
00203
00204 tmp = -m+sqrt(m*m-3.0*sqrt(3.0)*m*xi+9.0);
00205 rhoe = pow(1.0/3.0, gamma)*sqrt(2.0/3.0)*pow(tmp, gamma);
00206
00207 db0dxi = (drhoedxi*rhoc-drhocdxi*rhoe)/rhoc/rhoc;
00208 for(i=0; i < 3; i++)
00209 dxids[i][i] = 1.0/rc/sqrt(3.0);
00210 cmulm(db0dxi,dxids,db0ds);
00211
00212
00213 dd0dth = -2.0*(c1+c2)*cos(theta)*sin(theta)+c3*(cos(theta)*cos(theta)-sin(theta)*sin(theta));
00214 dd1dth = 2.0*(-c4*sqrt(3.0)*sin(theta)-c5*cos(theta));
00215 dqdth = dqdd0*dd0dth+dqdd1*dd1dth;
00216
00217
00218 for (i=0; i < 3; i++)
00219 {
00220 for (j=0; j < 3; j++)
00221 {
00222 for (k = 0; k < 3; k++)
00223 dj3ds[i][j] += dev[i][k]*dev[k][j];
00224 }
00225 dj3ds[i][i] -= 2.0/3.0*j2s;
00226 }
00227 if(fabs(fabs(theta/(pi/6.0))-1.0)<epsi) {
00228 for (i = 0; i < 3; i++) {
00229 for (j = 0; j < 3; j++)
00230 dthds[i][j] = 0.0 ;
00231 }
00232 }
00233 else {
00234 dthdj2 = 3.0 * sqrt(3.0)*j3s/4.0/cos(3.0*theta)/pow(j2s,2.5);
00235 dthdj3 = -sqrt(3.0)/2.0/cos(3.0*theta)/pow(j2s,1.5);
00236 for (i = 0; i < 3; i++) {
00237 for (j = 0; j < 3; j++)
00238 dthds[i][j] = dthdj2*dev[i][j] + dthdj3*dj3ds[i][j];
00239 }
00240 }
00241
00242
00243
00244 for (i = 0; i < 3; i++) {
00245 for (j = 0; j < 3; j++) {
00246 drds[i][j] = dqdb0*db0ds[i][j]+dqdth*dthds[i][j];
00247 }
00248 }
00249 cmulm((2.0/r)*khat*rhoc*rhoc,drds);
00250 }
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264 void shefplast::derkhatds (long ipp,matrix &sig, vector &q, matrix &dkhatds, long ido)
00265 {
00266 long i;
00267 double k, dkhatdxi;
00268
00269 matrix dxids(3,3);
00270
00271 k = k0 + (1.0-k0)*sqrt(q[0]*(2.0-q[0]));
00272
00273 dkhatdxi = -2.0*pow(k,2.0*p)*(1.0-k)*(1.0-k)*xi/a/a;
00274
00275 for(i=0; i < 3; i++)
00276 dxids[i][i] = 1.0/rc/sqrt(3.0);
00277
00278
00279 cmulm(dkhatdxi, dxids, dkhatds);
00280 cmulm(rhoc*rhoc/r/r, dkhatds);
00281 }
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295 void shefplast::deryieldfsigma (long ipp,matrix &sig, vector &q, matrix &dfds,long ido)
00296 {
00297 long i;
00298 double tmp, tmp1, rhoc;
00299 matrix drhods(3,3);
00300 matrix drhocds(3,3);
00301 matrix dqds(3,3);
00302 matrix dkhatds(3,3);
00303 matrix drhoeds(3,3);
00304
00305 compute_classvar (sig, q);
00306
00307
00308 cmulm(2.0/rc/rc,dev,drhods);
00309
00310 tmp = sqrt(m*m-12.0*sqrt(3.0)*m*xi+36.0);
00311 tmp1 = -gamma*pow((tmp-m)/6.0,gamma-1.0);
00312 tmp1 *= m*sqrt(2.0)/tmp;
00313 tmp1 /= (sqrt(3.0)*rc);
00314 for (i = 0; i < 3; i++)
00315 drhocds[i][i] = tmp1;
00316 rhoc = pow((1.0/6.0),gamma)*sqrt(2.0/3.0)*pow((tmp-m),gamma);
00317 cmulm(2.0*rhoc*khat/r/r,drhocds);
00318
00319
00320 derqdsigma(sig, dqds);
00321
00322 derkhatds(ipp,sig, q, dkhatds, ido);
00323
00324 addm(dkhatds, drhocds, dfds);
00325 addm(dfds, dqds, dfds);
00326 subm(drhods, dfds, dfds);
00327 }
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338 void shefplast::dderyieldfsigma (matrix &ddfds, strastrestate ssst)
00339 {
00340 }
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351 void shefplast::derpotsigma (long ipp,matrix &sig, vector &q, matrix &dgds,long ido)
00352 {
00353 deryieldfsigma (ipp,sig, q, dgds,ido);
00354 }
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367 void shefplast::deryieldfq(matrix &sig, vector &q, vector &dfq)
00368 {
00369 double k ;
00370 double dkdkh, dkhatdk ;
00371
00372 k = k0 + (1.0-k0)*sqrt(q[0]*(2.0-q[0]));
00373 dkhatdk = 2.0*p*pow(k,2.0*p-1.0)*(1.0-xi*xi*(1.0-k)*(1.0-k)/a/a)
00374 +2.0*(1.0-k)*(pow(k,p)*xi/a)*(pow(k,p)*xi/a) ;
00375
00376 dkdkh = (1.0-k0)*(1.0-q[0])/pow(q[0]*(2.0-q[0]),0.5) ;
00377
00378 dfq[0] = -dkhatdk*dkdkh*rhoc*rhoc/r/r;
00379 return;
00380 }
00381
00382
00383
00384
00385
00386
00387
00388
00389 void shefplast::matstiff (matrix &d, long ipp,long ido)
00390 {
00391 if (Mp->nlman->stmat==initial_stiff) {
00392
00393 Mm->elmatstiff (d,ipp);
00394 }
00395 if (Mp->nlman->stmat==tangent_stiff) {
00396
00397
00398 matrix ad(d.m,d.n);
00399 Mm->elmatstiff (ad,ipp);
00400 tangentstiff (ad,d,ipp,ido);
00401 }
00402 }
00403
00404
00405
00406
00407
00408
00409
00410 void shefplast::tangentstiff (matrix &d,matrix &td,long ipp,long ido)
00411 {
00412 long ncomp=Mm->ip[ipp].ncompstr;
00413 long i,j;
00414 double tmps,gamma,hf,dhdk;
00415 vector mxi,tmpv,str,dfdk(1),dhds,dfdsdk,dfds,av,q(1);
00416 matrix td7;
00417 matrix a,ainv,dinv,sig,dfdstens,dfdsds,am,tmp;
00418
00419 gamma=Mm->ip[ipp].eqother[ido+ncomp];
00420
00421 if (gamma<1.0e-10) {
00422 copym (d,td);
00423 }
00424 else {
00425 allocv (ncomp,str);
00426 allocv (7,mxi);
00427 allocv (7,tmpv);
00428 allocv (6,av);
00429 allocv (6,dhds);
00430 allocv (6,dfdsdk);
00431 allocv (6,dfds);
00432
00433 allocm (3,3,sig);
00434 allocm(6,6,dinv);
00435 allocm(7,7,a);
00436 allocm(7,7,tmp);
00437 allocm(7,7,ainv);
00438 allocm(3,3,dfdstens);
00439 allocm(6,6,dfdsds);
00440 allocm(6,6,am);
00441 allocm(d.m+1,d.n+1,td7);
00442
00443 invm(d,dinv,Mp->zero);
00444
00445 q[0] = Mm->ip[ipp].eqother[ido+ncomp+1];
00446
00447 Mm->givestress (0,ipp,str);
00448 vector_tensor (str,sig,Mm->ip[ipp].ssst,stress);
00449
00450 deryieldfsigma (ipp,sig,q,dfdstens,ido);
00451 tensor_vector (dfds,dfdstens,Mm->ip[ipp].ssst,stress);
00452 numdiff_dfdsdkc (ipp,sig,q,dfdsdk,ido);
00453 numdiff_dfdsdsc (ipp,sig,q,dfdsds,ido);
00454
00455 hf = hardening(ipp,sig,q,ido);
00456
00457 numdiff_dhdkc (ipp,sig,q,dhdk,ido);
00458 numdiff_dhdsc (ipp,sig,q,dhds,ido);
00459 deryieldfq (sig,q,dfdk);
00460
00461
00462 cmulm(gamma,dfdsds,am);
00463 addm(am,dinv,am);
00464
00465 for (i=0;i<6;i++) {
00466 for (j=0;j<6;j++) {
00467 a[i][j]=am[i][j];
00468 }
00469 }
00470
00471
00472 cmulv(gamma,dfdsdk,av);
00473
00474 for (i=0;i<6;i++) {
00475 a[i][6]=av[i];
00476 mxi[i] = dfds[i] ;
00477 }
00478 mxi[6] = dfdk[0] ;
00479
00480
00481 cmulv(-1.0*gamma,dhds,av);
00482
00483 for (i=0;i<6;i++) {
00484 a[6][i]=av[i];
00485 }
00486
00487
00488 a[6][6]=1.0-gamma*dhdk;
00489
00490 invm(a,ainv,Mp->zero);
00491 vxm(mxi,ainv,tmpv);
00492
00493
00494 mxi[6] = -1.0* hf ;
00495 scprd(tmpv,mxi,tmps) ;
00496 cmulv(1.0/tmps,tmpv,tmpv);
00497 vxv(mxi,tmpv,tmp);
00498 mxm(ainv,tmp,td7);
00499 subm(ainv,td7,td7);
00500
00501
00502 for (i=0;i<d.m;i++) {
00503 for (j=0;j<d.n;j++) {
00504 td[i][j] = td7[i][j] ;
00505 }
00506 }
00507 destrm(tmp);destrm(dinv);destrm(sig);destrm(dfdstens);destrm(dfdsds);destrm(am);
00508 destrm(a);destrm(ainv);destrm(td7);
00509
00510 destrv(dfds);destrv(str);destrv(dfdsdk);destrv(dhds);destrv(dfdk);
00511 destrv(av);destrv(mxi);destrv(tmpv);
00512 }
00513 }
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524 void shefplast::nlstresses (long ipp, long ido)
00525
00526 {
00527 long i,n=Mm->ip[ipp].ncompstr;
00528 double gamma;
00529 vector epsn(n),epsp(n),q(1);
00530
00531
00532 for (i=0; i<n; i++) {
00533 epsn[i]=Mm->ip[ipp].strain[i];
00534 epsp[i]=Mm->ip[ipp].eqother[ido+i];
00535 }
00536 gamma=Mm->ip[ipp].eqother[ido+n];
00537 q[0] = Mm->ip[ipp].eqother[ido+n+1];
00538 if (q[0] == 0.0)
00539 q[0] = kh0;
00540
00541
00542
00543 stress_return (ipp,gamma,q,epsn,epsp,ido);
00544
00545
00546 for (i=0; i<n; i++){
00547 Mm->ip[ipp].other[ido+i]=epsp[i];
00548 }
00549 Mm->ip[ipp].other[ido+n]=gamma;
00550 Mm->ip[ipp].other[ido+n+1]=q[0];
00551 }
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561 void shefplast::updateval (long ipp,long ido)
00562 {
00563 long i,n = Mm->ip[ipp].ncompstr;
00564
00565 for (i=0;i<n;i++){
00566 Mm->ip[ipp].eqother[ido+i]=Mm->ip[ipp].other[ido+i];
00567 }
00568
00569 Mm->ip[ipp].eqother[ido+n]=Mm->ip[ipp].other[ido+n];
00570 Mm->ip[ipp].eqother[ido+n+1]=Mm->ip[ipp].other[ido+n+1];
00571 }
00572
00573 void shefplast::stress_return (long ipp,double &lambda,vector &k,vector &eps,vector &epsp,long ido)
00574 {
00575 long i,ii,jj,stop,ni;
00576 double f,hf,nlambda,dlambda,dhdk,nory,err,tmp1,tmp2;
00577 vector dfds,dfdk,dfdsdk,dhds,nsig,dsig,mxi,tmp,x,y,av,sigtrial,nk,tmp6;
00578 matrix d,dinv,dfdsds,a,ainv,am,sigtens,nsigtens,dfdstens;
00579
00580
00581
00582
00583 ni=sra.give_ni ();
00584 err=sra.give_err ();
00585
00586 allocv (6,dfds);
00587 allocv (6,dfdsdk);
00588 allocv (6,dhds);
00589 allocv (6,nsig);
00590 allocv (6,dsig);
00591 allocv (6,sigtrial);
00592 allocv (7,mxi);
00593 allocv (7,tmp);
00594 allocv (6,tmp6);
00595 allocv (7,x);
00596 allocv (7,y);
00597 allocv (6,av);
00598 allocv (1,nk);
00599 allocv (1,dfdk);
00600
00601 allocm (6,6,d);
00602 allocm (6,6,dinv);
00603 allocm (6,6,dfdsds);
00604 allocm (7,7,a);
00605 allocm (7,7,ainv);
00606 allocm (6,6,am);
00607 allocm (3,3,sigtens);
00608 allocm (3,3,nsigtens);
00609 allocm (3,3,dfdstens);
00610
00611 nlambda = 0.0;
00612 nk[0] = k[0];
00613
00614
00615 Mm->elmatstiff (d,ipp);
00616
00617
00618
00619
00620
00621
00622 invm(d,dinv,Mp->zero);
00623
00624
00625 subv(eps,epsp,av);
00626
00627
00628 mxv(d,av,sigtrial);
00629
00630 vector_tensor(sigtrial,nsigtens,Mm->ip[ipp].ssst,stress);
00631 copyv(sigtrial,nsig);
00632
00633
00634
00635 stop=0;
00636 for (i=0;i<ni;i++) {
00637
00638 f=yieldfunction (nsigtens,nk);
00639
00640 if (f<err && i==0) {
00641 stop =1 ;
00642 break ;
00643 }
00644
00645
00646 deryieldfsigma (ipp,nsigtens,nk,dfdstens,ido);
00647 tensor_vector (dfds,dfdstens,Mm->ip[ipp].ssst,stress);
00648
00649
00650 deryieldfq (nsigtens,nk,dfdk);
00651
00652
00653 numdiff_dfdsdsc (ipp,nsigtens,nk,dfdsds,ido);
00654
00655
00656 numdiff_dfdsdkc (ipp,nsigtens,nk,dfdsdk,ido);
00657
00658
00659 hf = hardening(ipp,nsigtens,nk,ido);
00660
00661
00662 numdiff_dhdsc (ipp,nsigtens,nk,dhds,ido);
00663
00664
00665 numdiff_dhdkc (ipp,nsigtens,nk,dhdk,ido);
00666
00667
00668
00669
00670
00671
00672 cmulm(nlambda,dfdsds,am);
00673 addm(am,dinv,am);
00674
00675 for (ii=0;ii<6;ii++) {
00676 for (jj=0;jj<6;jj++) {
00677 a[ii][jj]=am[ii][jj];
00678 }
00679 }
00680
00681
00682 cmulv(nlambda,dfdsdk,av);
00683
00684 for (ii=0;ii<6;ii++) {
00685 a[ii][6]=av[ii];
00686 mxi[ii] = dfds[ii] ;
00687 }
00688 mxi[6] = dfdk[0] ;
00689
00690
00691 cmulv(-1.0*nlambda,dhds,av);
00692
00693 for (ii=0;ii<6;ii++) {
00694 a[6][ii]=av[ii];
00695 }
00696
00697
00698 a[6][6]=1.0-nlambda*dhdk;
00699
00700
00701
00702 cmulv(nlambda,dfds,av);
00703 mxv(dinv,nsig,tmp6);
00704 addv(tmp6,av,av);
00705
00706 mxv(dinv,sigtrial,tmp6);
00707 subv(av,tmp6,av);
00708
00709 for (ii=0;ii<6;ii++) {
00710 y[ii]=av[ii];
00711 }
00712
00713
00714 if(nk[0]==1.0)
00715 y[6] = 0.0 ;
00716 else
00717 y[6]=nk[0]-nlambda*hf-k[0];
00718
00719
00720 invm(a,ainv,Mp->zero);
00721 vxm(mxi,ainv,tmp);
00722 scprd(tmp,y,tmp1) ;
00723
00724 mxi[6] = -1.0* hf ;
00725 scprd(tmp,mxi,tmp2) ;
00726 dlambda = (f-tmp1)/tmp2 ;
00727
00728 cmulv(dlambda,mxi,mxi);
00729 addv(y,mxi,tmp);
00730 cmulv(-1.0,tmp,tmp);
00731 mxv(ainv,tmp,x) ;
00732
00733
00734
00735
00736 for (ii=0;ii<6;ii++) {
00737 dsig[ii]=x[ii];
00738 }
00739
00740
00741 addv(nsig,dsig,nsig);
00742 vector_tensor (nsig,nsigtens,Mm->ip[ipp].ssst,stress);
00743
00744
00745 nk[0]+=x[6];
00746 if(nk[0] < 0.0) {
00747 nk[0] = 0.0 ;
00748 }
00749 if(nk[0] > 1.0) {
00750 nk[0] = 1.0 ;
00751 }
00752
00753 nlambda+=dlambda;
00754
00755
00756 nory = fabs(y[0]) ;
00757 for (ii=1;ii<7;ii++) {
00758 nory = maxim(nory,fabs(y[ii]));
00759 }
00760 nory = maxim(nory,f);
00761
00762 if(ipp==0)
00763 fprintf(stderr,"local error %e/%e \n",nory,err);
00764 if (nory<err) {
00765 stop=1;
00766 break;
00767 }
00768 }
00769
00770 if (stop==1) {
00771 copyv (nk,k);
00772 Mm->storestress (0,ipp,nsig);
00773 lambda=nlambda;
00774 mxv(dinv,nsig,tmp6);
00775 subv(eps,tmp6,epsp);
00776 }
00777 if(i==ni) {
00778 fprintf (stderr,"\n\n stress return algo in %s is not successfull after %ld steps for ipp = %ld.\n",__FILE__,ni,ipp);
00779 }
00780
00781 destrm (dfdstens); destrm (nsigtens); destrm (sigtens);
00782 destrm (am); destrm (a); destrm (dfdsds);
00783 destrm (d); destrm (dinv);
00784
00785 destrv (dfdk); destrv (nk); destrv (av);
00786 destrv (y); destrv (x); destrv (sigtrial);
00787 destrv (nsig); destrv (dhds); destrv (dfdsdk);
00788 destrv (dfds);
00789 }
00790
00791
00792
00793
00794
00795
00796 void shefplast::compzeta()
00797 {
00798 if (xi>0.0)
00799
00800 zeta = 2.0e-5;
00801 else
00802 zeta = -ah + sqrt(ah*ah-bh*xi+ch);
00803 }
00804
00805 double shefplast::hardening(long ipp,matrix &sigtens,vector &q,long ido)
00806 {
00807 double h;
00808 matrix dfds;
00809
00810 if (q[0]<1.0) {
00811 allocm (3,3,dfds);
00812 deryieldfsigma (ipp,sigtens,q,dfds,ido);
00813 compzeta ();
00814
00815 h = sqrt(2.0/3.0)*tensornorm (dfds)/zeta;
00816 destrm (dfds);
00817 }
00818 else {
00819 h=0.0;
00820 }
00821
00822 return h;
00823 }
00824
00825
00826
00827
00828
00829
00830
00831
00832
00833
00834 void shefplast::numdiff_dfdsdsc(long ipp,matrix &sigtens,vector &q,matrix &dfdsds,long ido)
00835 {
00836 long i,j;
00837 double dh;
00838 long ncomp = Mm->ip[ipp].ncompstr;
00839 vector sig(ncomp),sigdh(ncomp),sigmdh(ncomp);
00840 matrix dfds(3,3),dfdsdh(3,3),sigtensdh(3,3),sigtensmdh(3,3);
00841
00842 tensor_vector (sig,sigtens,Mm->ip[ipp].ssst,stress);
00843
00844 for (i=0;i<6;i++) {
00845 copyv (sig,sigdh);
00846 copyv (sig,sigmdh);
00847 dh = normv(sig);
00848 dh *= sqrt(DBL_EPSILON);
00849
00850 sigdh[i]+=dh;
00851 sigmdh[i]-=dh;
00852 vector_tensor (sigdh,sigtensdh,Mm->ip[ipp].ssst,stress);
00853 vector_tensor (sigmdh,sigtensmdh,Mm->ip[ipp].ssst,stress);
00854
00855 deryieldfsigma (ipp,sigtensmdh,q,dfds,ido);
00856 deryieldfsigma (ipp,sigtensdh,q,dfdsdh,ido);
00857
00858 subm(dfdsdh,dfds,dfds);
00859 cmulm(0.5/dh,dfds,dfds);
00860
00861 tensor_vector (sigdh,dfds,Mm->ip[ipp].ssst,stress);
00862
00863 for (j=0;j<6;j++) {
00864 dfdsds[j][i]=sigdh[j];
00865 }
00866 }
00867 }
00868
00869
00870
00871
00872
00873
00874
00875
00876
00877
00878 void shefplast::numdiff_dfdsdkc(long ipp,matrix &sigtens,vector &q,vector &dfdsdk,long ido)
00879 {
00880 double dh;
00881 vector qdh,qmdh;
00882 matrix dfds,dfdsdh,dfdsmdh;
00883
00884 allocv (1,qdh);
00885 allocv (1,qmdh);
00886 allocm (3,3,dfds);
00887 allocm (3,3,dfdsdh);
00888 allocm (3,3,dfdsmdh);
00889
00890 dh = normv(q);
00891 dh *= sqrt(DBL_EPSILON);
00892
00893 qdh[0]=q[0]+dh;
00894 qmdh[0]=q[0]-dh;
00895
00896 deryieldfsigma (ipp,sigtens,qmdh,dfds,ido);
00897 deryieldfsigma (ipp,sigtens,qdh,dfdsdh,ido);
00898
00899 subm(dfdsdh,dfds,dfds);
00900 cmulm(0.5/dh,dfds,dfds);
00901
00902 tensor_vector (dfdsdk,dfds,Mm->ip[ipp].ssst,stress);
00903
00904 destrm (dfdsdh); destrm (dfdsmdh); destrm (dfds);
00905 destrv (qdh); destrv (qmdh);
00906 }
00907
00908
00909
00910
00911
00912
00913
00914 void shefplast::numdiff_dhdsc(long ipp,matrix &sigtens,vector &q,vector &dhds,long ido)
00915 {
00916 long i;
00917 long ncomp = Mm->ip[ipp].ncompstr;
00918 double h,hdh,dh;
00919 vector sig,sigdh,sigmdh;
00920 matrix sigtensdh;
00921
00922 allocv (ncomp,sig);
00923 allocv (ncomp,sigdh);
00924 allocv (ncomp,sigmdh);
00925 allocm (3,3,sigtensdh);
00926
00927 tensor_vector (sig,sigtens,Mm->ip[ipp].ssst,stress);
00928
00929 for (i=0;i<6;i++) {
00930 copyv (sig,sigdh);
00931 copyv (sig,sigmdh);
00932 dh = normv(sig);
00933 dh *= sqrt(DBL_EPSILON);
00934
00935 sigdh[i]+=dh;
00936 sigmdh[i]-=dh;
00937
00938 vector_tensor (sigdh,sigtensdh,Mm->ip[ipp].ssst,stress);
00939 vector_tensor (sigmdh,sigtens,Mm->ip[ipp].ssst,stress);
00940
00941 hdh = hardening (ipp,sigtensdh,q,ido);
00942 h = hardening (ipp,sigtens,q,ido);
00943
00944 dhds[i]=(hdh-h)*0.5/dh;
00945 }
00946
00947 destrm (sigtensdh);
00948 destrv (sigdh); destrv (sigmdh); destrv (sig);
00949 }
00950
00951
00952
00953
00954
00955
00956
00957 void shefplast::numdiff_dhdkc(long ipp,matrix &sigtens,vector &q,double &dhdk,long ido)
00958 {
00959 double h,hdh,dh;
00960 vector qdh,qmdh;
00961
00962 allocv (1,qdh);
00963 allocv (1,qmdh);
00964
00965 dh = normv(q);
00966 dh *= sqrt(DBL_EPSILON);
00967 qdh[0]=q[0]+dh;
00968 qmdh[0]=q[0]-dh;
00969
00970 hdh = hardening (ipp,sigtens,qdh,ido);
00971 h = hardening (ipp,sigtens,qmdh,ido);
00972
00973 dhdk = (hdh-h)*0.5/dh;
00974
00975 destrv (qdh);
00976 destrv (qmdh);
00977 }
00978
00979
00980 inline double
00981 shefplast::maxim (double a,double b)
00982 {
00983 return(a>b ? a:b);
00984 }