00001 #include <math.h>
00002 #include "chen.h"
00003 #include "global.h"
00004 #include "genfile.h"
00005 #include "intpoints.h"
00006 #include "vecttens.h"
00007 #include "alias.h"
00008 #include "stochdriver.h"
00009 #include "matrix.h"
00010
00011 chen::chen (void)
00012 {
00013 fyc=0.0;
00014 fyt=0.0;
00015 fybc=0.0;
00016 fc=0.0;
00017 ft=0.0;
00018 fbc=0.0;
00019 epsu=0.0;
00020
00021 ay=0.0; au=0.0;
00022 ky=0.0; ku=0.0;
00023 alpha=0.0; beta=0.0;
00024
00025
00026
00027
00028 hard=1;
00029 }
00030
00031 chen::~chen (void)
00032 {
00033
00034 }
00035
00036
00037
00038
00039
00040
00041 void chen::read (FILE *in)
00042 {
00043 fscanf (in,"%lf %lf %lf %lf %lf %lf %lf",&fyc,&fyt,&fybc,&fc,&ft,&fbc,&epsu);
00044 sra.read (in);
00045 }
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055 void chen::matstiff (matrix &d,long ipp,long ido)
00056 {
00057 if (Mp->nlman->stmat==initial_stiff){
00058
00059 Mm->elmatstiff (d,ipp);
00060 }
00061 if (Mp->nlman->stmat==tangent_stiff){
00062
00063
00064
00065 matrix ad(d.m,d.n);
00066 Mm->elmatstiff (ad,ipp);
00067 tangentstiff (ad,d,ipp,ido);
00068 }
00069 }
00070
00071
00072 void chen::tangentstiff (matrix &d,matrix &td,long ipp,long ido)
00073 {
00074 long ncomp=Mm->ip[ipp].ncompstr;
00075 double denom,gamma;
00076 vector str,av(d.m),q(1);
00077 matrix sig(3,3),am(d.m,d.n);
00078
00079 gamma=Mm->ip[ipp].eqother[ido+ncomp];
00080 if (gamma<1.0e-10){
00081 copym (d,td);
00082 }
00083 else{
00084
00085 allocv (ncomp,str);
00086
00087 Mm->givestress (0,ipp,str);
00088 vector_tensor (str,sig,Mm->ip[ipp].ssst,stress);
00089 deryieldfdsigma (sig,q,sig);
00090 tensor_vector (str,sig,Mm->ip[ipp].ssst,stress);
00091
00092 if (Mm->ip[ipp].ssst==planestress){
00093 vector auxstr(3);
00094 auxstr[0]=str[0];auxstr[1]=str[1];auxstr[2]=str[2];
00095 destrv (str);
00096 allocv (d.m,str);
00097 str[0]=auxstr[0];str[1]=auxstr[1];str[2]=auxstr[2];
00098 }
00099
00100 mxv (d,str,av);
00101 scprd (av,str,denom);
00102
00103
00104 q[0] = Mm->ip[ipp].eqother[ido+ncomp+1];
00105 denom+= plasmodscalar(q);
00106
00107 if (fabs(denom)<1.0e-10){
00108 copym (d,td);
00109 }
00110 else{
00111 vxv (str,str,am);
00112 mxm (d,am,td);
00113 mxm (td,d,am);
00114
00115 cmulm (1.0/denom,am);
00116
00117 subm (d,am,td);
00118 }
00119 }
00120
00121 }
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131 double chen::yieldfunction (matrix &sig,vector &q)
00132 {
00133 double f,invar,j2,zone,kappa;
00134 matrix dev(3,3);
00135
00136 deviator (sig,dev);
00137 invar=first_invar (sig);
00138 j2=second_invar(dev);
00139
00140 zone=sqrt(j2)+invar/sqrt(3.0);
00141
00142
00143 if(invar<0.0 && zone<0.0)
00144 {
00145 ay = (fybc*fybc-fyc*fyc)/(2.0*fybc-fyc);
00146 au = (fbc*fbc-fc*fc)/(2.0*fbc-fc);
00147 ky = fyc*fybc*(2*fyc-fybc)/3.0/(2*fybc-fyc);
00148 ku = fc*fbc*(2*fc-fbc)/3.0/(2*fbc-fc);
00149 alpha = (au-ay)/(ku-ky);
00150 beta = (ay*ku-au*ky)/(ku-ky);
00151
00152 kappa = sqrt(ky + q[0]*(ku-ky));
00153 if (kappa>sqrt(ku))
00154 kappa=sqrt(ku);
00155
00156 if (hard==0){
00157
00158 f = j2 + ay/3.0*invar-ky;
00159 }
00160 if (hard==1){
00161
00162 f = j2 + beta/3.0*invar-kappa*kappa*(1.0-alpha/3.0*invar);
00163 }
00164 }
00165
00166 else
00167 {
00168 ay = (fyc-fyt)/2.0;
00169 au = (fc-ft)/2.0;
00170 ky = fyc*fyt/6.0;
00171 ku = fc*ft/6.0;
00172 alpha = (au-ay)/(ku-ky);
00173 beta = (ay*ku-au*ky)/(ku-ky);
00174
00175 kappa = sqrt(ky + q[0]*(ku-ky));
00176 if (kappa>sqrt(ku))
00177 kappa=sqrt(ku);
00178
00179 if (hard==0){
00180
00181 f = j2 - (invar*invar)/6.0+ay/3.0*invar-ky;
00182 }
00183 if (hard==1){
00184
00185 f = j2 - (invar*invar)/6.0 + beta/3.0*invar-kappa*kappa*(1.0-alpha/3.0*invar);
00186 }
00187
00188 }
00189
00190
00191
00192
00193 return f;
00194 }
00195
00196
00197
00198
00199
00200
00201
00202 void chen::deryieldfdsigma (matrix &sig,vector &q,matrix &dfds)
00203 {
00204 double invar,j2,zone,kappa;
00205 matrix dev(3,3);
00206
00207 deviator (sig,dev);
00208 invar=first_invar (sig);
00209 j2=second_invar(dev);
00210
00211 zone=sqrt(j2)+invar/sqrt(3.0);
00212
00213
00214 if(invar<0.0 && zone<0.0)
00215 {
00216 ay = (fybc*fybc-fyc*fyc)/(2.0*fybc-fyc);
00217 au = (fbc*fbc-fc*fc)/(2.0*fbc-fc);
00218 ky = fyc*fybc*(2*fyc-fybc)/3.0/(2*fybc-fyc);
00219 ku = fc*fbc*(2*fc-fbc)/3.0/(2*fbc-fc);
00220 alpha = (au-ay)/(ku-ky);
00221 beta = (ay*ku-au*ky)/(ku-ky);
00222
00223 kappa = sqrt(ky + q[0]*(ku-ky));
00224 if (kappa>sqrt(ku))
00225 kappa=sqrt(ku);
00226
00227 if (hard==0){
00228 dfds[0][0]=1.0/3.0*(2.0*sig[0][0]-sig[1][1]-sig[2][2]) + beta/3.0+alpha/3.0*kappa*kappa;
00229 dfds[1][1]=1.0/3.0*(2.0*sig[1][1]-sig[0][0]-sig[2][2]) + beta/3.0+alpha/3.0*kappa*kappa;
00230 dfds[2][2]=1.0/3.0*(2.0*sig[2][2]-sig[0][0]-sig[1][1]) + beta/3.0+alpha/3.0*kappa*kappa;
00231 }
00232 if (hard==1){
00233 dfds[0][0]=1.0/3.0*(2.0*sig[0][0]-sig[1][1]-sig[2][2]) + ay/3.0;
00234 dfds[1][1]=1.0/3.0*(2.0*sig[1][1]-sig[0][0]-sig[2][2]) + ay/3.0;
00235 dfds[2][2]=1.0/3.0*(2.0*sig[2][2]-sig[0][0]-sig[1][1]) + ay/3.0;
00236 }
00237
00238 dfds[0][1]=2*sig[0][1];
00239 dfds[1][0]=2*sig[0][1];
00240 dfds[0][2]=2*sig[0][2];
00241 dfds[2][0]=2*sig[0][2];
00242 dfds[1][2]=2*sig[1][2];
00243 dfds[2][1]=2*sig[1][2];
00244 }
00245
00246
00247
00248
00249 else
00250 {
00251 ay = (fyc-fyt)/2.0;
00252 au = (fc-ft)/2.0;
00253 ky = fyc*fyt/6.0;
00254 ku = fc*ft/6.0;
00255 alpha = (au-ay)/(ku-ky);
00256 beta = (ay*ku-au*ky)/(ku-ky);
00257
00258 kappa = sqrt(ky + q[0]*(ku-ky));
00259 if (kappa>sqrt(ku))
00260 kappa=sqrt(ku);
00261
00262 if (hard==0){
00263 dfds[0][0]=1.0/3.0*(2.0*sig[0][0]-sig[1][1]-sig[2][2]) + beta/3.0+alpha/3.0*kappa*kappa;
00264 dfds[1][1]=1.0/3.0*(2.0*sig[1][1]-sig[0][0]-sig[2][2]) + beta/3.0+alpha/3.0*kappa*kappa;
00265 dfds[2][2]=1.0/3.0*(2.0*sig[2][2]-sig[0][0]-sig[1][1]) + beta/3.0+alpha/3.0*kappa*kappa;
00266 }
00267 if (hard==1){
00268 dfds[0][0]=1.0/3.0*(2.0*sig[0][0]-sig[1][1]-sig[2][2]) + ay/3.0;
00269 dfds[1][1]=1.0/3.0*(2.0*sig[1][1]-sig[0][0]-sig[2][2]) + ay/3.0;
00270 dfds[2][2]=1.0/3.0*(2.0*sig[2][2]-sig[0][0]-sig[1][1]) + ay/3.0;
00271 }
00272
00273 dfds[0][1]=2*sig[0][1];
00274 dfds[1][0]=2*sig[0][1];
00275 dfds[0][2]=2*sig[0][2];
00276 dfds[2][0]=2*sig[0][2];
00277 dfds[1][2]=2*sig[1][2];
00278 dfds[2][1]=2*sig[1][2];
00279
00280 dfds[0][0]-=1.0/3.0*invar;
00281 dfds[1][1]-=1.0/3.0*invar;
00282 dfds[2][2]-=1.0/3.0*invar;
00283
00284
00285 }
00286
00287 }
00288
00289
00290
00291
00292
00293
00294 void chen::deryieldfdsigma_old (matrix &sig,matrix &dfds)
00295 {
00296 double invar,j2,zone;
00297 matrix dev(3,3);
00298
00299 deviator (sig,dev);
00300 invar=first_invar (sig);
00301 j2=second_invar(dev);
00302
00303 zone=sqrt(j2)+invar/sqrt(3.0);
00304
00305
00306 if(invar<0.0 && zone<0.0)
00307 {
00308 ay = (fybc*fybc-fyc*fyc)/(2.0*fybc-fyc);
00309
00310 dfds[0][0]=1.0/3.0*(2.0*sig[0][0]-sig[1][1]-sig[2][2])+ay/3.0;
00311 dfds[1][1]=1.0/3.0*(2.0*sig[1][1]-sig[0][0]-sig[2][2])+ay/3.0;
00312 dfds[2][2]=1.0/3.0*(2.0*sig[2][2]-sig[0][0]-sig[1][1])+ay/3.0;
00313
00314 dfds[0][1]=2*sig[0][1];
00315 dfds[1][0]=2*sig[0][1];
00316 dfds[0][2]=2*sig[0][2];
00317 dfds[2][0]=2*sig[0][2];
00318 dfds[1][2]=2*sig[1][2];
00319 dfds[2][1]=2*sig[1][2];
00320 }
00321
00322
00323
00324
00325 else
00326 {
00327 ay = (fyc-fyt)/2.0;
00328
00329 dfds[0][0]=1.0/3.0*(2.0*sig[0][0]-sig[1][1]-sig[2][2])+ay/3.0;
00330 dfds[1][1]=1.0/3.0*(2.0*sig[1][1]-sig[0][0]-sig[2][2])+ay/3.0;
00331 dfds[2][2]=1.0/3.0*(2.0*sig[2][2]-sig[0][0]-sig[1][1])+ay/3.0;
00332
00333 dfds[0][1]=2*sig[0][1];
00334 dfds[1][0]=2*sig[0][1];
00335 dfds[0][2]=2*sig[0][2];
00336 dfds[2][0]=2*sig[0][2];
00337 dfds[1][2]=2*sig[1][2];
00338 dfds[2][1]=2*sig[1][2];
00339
00340 dfds[0][0]-=1.0/3.0*invar;
00341 dfds[1][1]-=1.0/3.0*invar;
00342 dfds[2][2]-=1.0/3.0*invar;
00343
00344
00345 }
00346
00347 }
00348
00349
00350
00351
00352
00353
00354 void chen::deryieldfdsigmadsigma (matrix &sig,matrix &dfdsds)
00355 {
00356 double invar,j2,zone;
00357 matrix dev(3,3);
00358
00359 deviator (sig,dev);
00360 invar=first_invar (sig);
00361 j2=second_invar(dev);
00362
00363 fillm (0.0,dfdsds);
00364
00365 zone = sqrt(j2)+invar/sqrt(3.0);
00366
00367
00368 if(invar<0.0 && zone<0.0){
00369 dfdsds[0][0] = 2.0/3.0;
00370 dfdsds[0][1] = -1.0/3.0;
00371 dfdsds[0][2] = -1.0/3.0;
00372
00373 dfdsds[1][0] = -1.0/3.0;
00374 dfdsds[1][1] = 2.0/3.0;
00375 dfdsds[1][2] = -1.0/3.0;
00376
00377 dfdsds[2][0] = -1.0/3.0;
00378 dfdsds[2][1] = -1.0/3.0;
00379 dfdsds[2][2] = 2.0/3.0;
00380
00381 dfdsds[3][3] = 2.0;
00382 dfdsds[4][4] = 2.0;
00383 dfdsds[5][5] = 2.0;
00384 }
00385
00386
00387 else{
00388 dfdsds[0][0] = 1.0/3.0;
00389 dfdsds[0][1] = -2.0/3.0;
00390 dfdsds[0][2] = -2.0/3.0;
00391
00392 dfdsds[1][0] = -2.0/3.0;
00393 dfdsds[1][1] = 1.0/3.0;
00394 dfdsds[1][2] = -2.0/3.0;
00395
00396 dfdsds[2][0] = -2.0/3.0;
00397 dfdsds[2][1] = -2.0/3.0;
00398 dfdsds[2][2] = 1.0/3.0;
00399
00400 dfdsds[3][3] = 2.0;
00401 dfdsds[4][4] = 2.0;
00402 dfdsds[5][5] = 2.0;
00403 }
00404
00405 }
00406
00407
00408
00409
00410
00411
00412 void chen::deryieldfdsigma_old_old (matrix &sig,matrix &dfds)
00413 {
00414 double invar,j2,a0;
00415 matrix dev(3,3);
00416
00417 deviator (sig,dev);
00418 invar=first_invar (sig);
00419 normedtensor (dev,dfds);
00420 j2=second_invar(dev);
00421
00422
00423 if(invar<0.0 && (sqrt(j2)+invar/sqrt(3.0))<0.0)
00424 {
00425 a0 = (fybc*fybc-fyc*fyc)/(2.0*fybc-fyc);
00426 dfds[0][0]+=a0/3.0;
00427 dfds[1][1]+=a0/3.0;
00428 dfds[2][2]+=a0/3.0;
00429 }
00430
00431
00432
00433
00434 else
00435 {
00436 a0 = (fyc-fyt)/2.0;
00437 dfds[0][0]+=a0/3.0-1.0/3.0*invar;
00438 dfds[1][1]+=a0/3.0-1.0/3.0*invar;
00439 dfds[2][2]+=a0/3.0-1.0/3.0*invar;
00440 }
00441
00442 }
00443
00444
00445
00446
00447
00448
00449 void chen::deryieldfdq (matrix &sig,vector &q,vector &dfdq)
00450 {
00451 double invar,j2,zone,kappa;
00452 matrix dev(3,3);
00453
00454 if (hard==0){
00455 nullv (dfdq.a,dfdq.n);
00456 }
00457 if (hard==1){
00458 deviator (sig,dev);
00459 invar=first_invar (sig);
00460 j2=second_invar(dev);
00461
00462 zone=sqrt(j2)+invar/sqrt(3.0);
00463
00464
00465 if(invar<0.0 && zone<0.0)
00466 {
00467 ay = (fybc*fybc-fyc*fyc)/(2.0*fybc-fyc);
00468 au = (fbc*fbc-fc*fc)/(2.0*fbc-fc);
00469 ky = fyc*fybc*(2*fyc-fybc)/3.0/(2*fybc-fyc);
00470 ku = fc*fbc*(2*fc-fbc)/3.0/(2*fbc-fc);
00471 alpha = (au-ay)/(ku-ky);
00472 beta = (ay*ku-au*ky)/(ku-ky);
00473 }
00474
00475 else{
00476 ay = (fyc-fyt)/2.0;
00477 au = (fc-ft)/2.0;
00478 ky = fyc*fyt/6.0;
00479 ku = fc*ft/6.0;
00480 alpha = (au-ay)/(ku-ky);
00481 beta = (ay*ku-au*ky)/(ku-ky);
00482 }
00483
00484 kappa = sqrt(ky + q[0]*(ku-ky));
00485 if (kappa>sqrt(ku))
00486 kappa=sqrt(ku);
00487
00488 dfdq[0]=-2.0*kappa*(1.0-alpha/3.0*invar);
00489 }
00490 }
00491
00492
00493
00494
00495
00496
00497 void chen::deryieldfdqdq (matrix &sig,vector &q,matrix &dfdqdq)
00498 {
00499 double invar,j2,zone,kappa;
00500 matrix dev(3,3);
00501
00502 if (hard==0){
00503 dfdqdq[0][0]=0.0;
00504 }
00505 if (hard==1){
00506 deviator (sig,dev);
00507 invar=first_invar (sig);
00508 j2=second_invar(dev);
00509
00510 zone=sqrt(j2)+invar/sqrt(3.0);
00511
00512
00513 if(invar<0.0 && zone<0.0)
00514 {
00515 ay = (fybc*fybc-fyc*fyc)/(2.0*fybc-fyc);
00516 au = (fbc*fbc-fc*fc)/(2.0*fbc-fc);
00517 ky = fyc*fybc*(2*fyc-fybc)/3.0/(2*fybc-fyc);
00518 ku = fc*fbc*(2*fc-fbc)/3.0/(2*fbc-fc);
00519 alpha = (au-ay)/(ku-ky);
00520 beta = (ay*ku-au*ky)/(ku-ky);
00521 }
00522 else{
00523 ay = (fyc-fyt)/2.0;
00524 au = (fc-ft)/2.0;
00525 ky = fyc*fyt/6.0;
00526 ku = fc*ft/6.0;
00527 alpha = (au-ay)/(ku-ky);
00528 beta = (ay*ku-au*ky)/(ku-ky);
00529 }
00530
00531 dfdqdq[0][0]=-2.0*(1.0-alpha/3.0*invar);
00532 }
00533 }
00534
00535
00536 void chen::deryieldfdsigmadq (matrix &sig,vector &q,matrix &dfdsdq)
00537 {
00538 double invar,j2,zone,kappa;
00539 matrix dev(3,3);
00540
00541 if (hard==0){
00542 fillm (0.0,dfdsdq);
00543 }
00544 if (hard==1){
00545 deviator (sig,dev);
00546 invar=first_invar (sig);
00547 j2=second_invar(dev);
00548
00549 zone=sqrt(j2)+invar/sqrt(3.0);
00550
00551
00552 if(invar<0.0 && zone<0.0)
00553 {
00554 ay = (fybc*fybc-fyc*fyc)/(2.0*fybc-fyc);
00555 au = (fbc*fbc-fc*fc)/(2.0*fbc-fc);
00556 ky = fyc*fybc*(2*fyc-fybc)/3.0/(2*fybc-fyc);
00557 ku = fc*fbc*(2*fc-fbc)/3.0/(2*fbc-fc);
00558 alpha = (au-ay)/(ku-ky);
00559 beta = (ay*ku-au*ky)/(ku-ky);
00560
00561 kappa = sqrt(ky + q[0]*(ku-ky));
00562 if (kappa>sqrt(ku))
00563 kappa=sqrt(ku);
00564
00565 dfdsdq[0][0]=2.0*alpha/3.0*kappa;
00566 dfdsdq[1][0]=2.0*alpha/3.0*kappa;
00567 dfdsdq[2][0]=2.0*alpha/3.0*kappa;
00568
00569 dfdsdq[3][0]=0.0;
00570 dfdsdq[4][0]=0.0;
00571 dfdsdq[5][0]=0.0;
00572
00573 }
00574
00575
00576
00577
00578 else
00579 {
00580 ay = (fyc-fyt)/2.0;
00581 au = (fc-ft)/2.0;
00582 ky = fyc*fyt/6.0;
00583 ku = fc*ft/6.0;
00584 alpha = (au-ay)/(ku-ky);
00585 beta = (ay*ku-au*ky)/(ku-ky);
00586
00587 kappa = sqrt(ky + q[0]*(ku-ky));
00588 if (kappa>sqrt(ku))
00589 kappa=sqrt(ku);
00590
00591 dfdsdq[0][0]=2.0*alpha/3.0*kappa;
00592 dfdsdq[1][0]=2.0*alpha/3.0*kappa;
00593 dfdsdq[2][0]=2.0*alpha/3.0*kappa;
00594
00595 dfdsdq[3][0]=0.0;
00596 dfdsdq[4][0]=0.0;
00597 dfdsdq[5][0]=0.0;
00598 }
00599 }
00600
00601 }
00602
00603
00604 void chen::plasmod (long ipp,vector &epsp,matrix &sig,matrix &h)
00605 {
00606 long ncompstr=epsp.n;
00607 double s,eq;
00608 matrix epspt(3,3);
00609
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620
00621
00622
00623
00624
00625
00626
00627
00628
00629
00630
00631
00632
00633
00634
00635
00636
00637
00638
00639
00640
00641
00642
00643 double f,invar,j2,zone;
00644 matrix dev(3,3);
00645
00646 deviator (sig,dev);
00647 invar=first_invar (sig);
00648 j2=second_invar(dev);
00649
00650 zone=sqrt(j2)+invar/sqrt(3.0);
00651
00652
00653 if(invar<0.0 && zone<0.0){
00654 ky = fyc*fybc*(2*fyc-fybc)/3.0/(2*fybc-fyc);
00655 ku = fc*fbc*(2*fc-fbc)/3.0/(2*fbc-fc);
00656 }
00657
00658 else{
00659 ky = fyc*fyt/6.0;
00660 ku = fc*ft/6.0;
00661 }
00662
00663 h[0][0]=(ku-ky)/epsu*epsu;
00664
00665
00666
00667
00668
00669
00670
00671
00672
00673
00674
00675
00676 }
00677
00678 double chen::plasmodscalar(vector &qtr)
00679 {
00680
00681
00682
00683
00684
00685
00686
00687
00688
00689
00690
00691 }
00692
00693 void chen::updateq(long ipp,double dgamma,vector &epsp,matrix &sig,vector &q)
00694 {
00695
00696
00697
00698
00699
00700
00701
00702
00703
00704
00705
00706
00707
00708
00709
00710
00711 }
00712
00713
00714
00715
00716
00717
00718
00719
00720
00721
00722
00723 void chen::nlstresses (long ipp, long im, long ido)
00724 {
00725 long i,ni,n,nhard;
00726 double gamma,err;
00727
00728
00729 n = Mm->ip[ipp].ncompstr;
00730
00731 vector epsn(n),epsp(n);
00732
00733
00734
00735
00736
00737
00738
00739
00740 vector q(1);
00741
00742
00743 for (i=0;i<n;i++){
00744
00745 epsn[i]=Mm->ip[ipp].strain[i];
00746
00747 epsp[i]=Mm->ip[ipp].eqother[ido+i];
00748 }
00749
00750 gamma = Mm->ip[ipp].eqother[ido+n];
00751
00752
00753 if (hard==0)
00754 q[0]=0.0;
00755 if (hard==1){
00756
00757 double f,invar,j2,zone,kappa;
00758 vector sig (epsp.n);
00759 matrix sigt(3,3),dev(3,3);
00760
00761
00762 vector_tensor (sig,sigt,Mm->ip[ipp].ssst,stress);
00763 deviator (sigt,dev);
00764 invar=first_invar (sigt);
00765 j2=second_invar(dev);
00766
00767 zone=sqrt(j2)+invar/sqrt(3.0);
00768
00769
00770 if(invar<0.0 && zone<0.0){
00771 ky = fyc*fybc*(2*fyc-fybc)/3.0/(2*fybc-fyc);
00772 ku = fc*fbc*(2*fc-fbc)/3.0/(2*fbc-fc);
00773 }
00774
00775 else{
00776 ky = fyc*fyt/6.0;
00777 ku = fc*ft/6.0;
00778 }
00779
00780 q[0] = Mm->ip[ipp].eqother[ido+n+1]+ky;
00781
00782 }
00783
00784
00785 if (sra.give_tsra () == cp){
00786 ni=sra.give_ni ();
00787 err=sra.give_err ();
00788
00789
00790
00791 Mm->newton_stress_return_2 (ipp,im,ido,gamma,epsn,epsp,q,ni,err);
00792
00793 }
00794 else{
00795 fprintf (stderr,"\n\n wrong type of stress return algorithm is required in nlstresses (file %s, line %d).\n",__FILE__,__LINE__);
00796 abort ();
00797 }
00798
00799
00800 for (i=0;i<n;i++){
00801 Mm->ip[ipp].other[ido+i]=epsp[i];
00802 }
00803 Mm->ip[ipp].other[ido+n]=gamma;
00804 if (hard==1){
00805 Mm->ip[ipp].other[ido+n+1]=q[0];
00806 }
00807
00808 }
00809
00810 void chen::nonloc_nlstresses (long ipp, long im, long ido)
00811 {
00812 long i,ni, n = Mm->ip[ipp].ncompstr;
00813 double gamma,err;
00814 vector epsn(n),epsp(n),q(1);
00815
00816
00817 for (i=0;i<n;i++){
00818 epsn[i]=Mm->ip[ipp].strain[i];
00819 epsp[i]=Mm->ip[ipp].nonloc[i];
00820 }
00821 gamma = Mm->ip[ipp].eqother[ido+n];
00822 q[0] = Mm->ip[ipp].eqother[ido+n+1];
00823
00824
00825 if (sra.give_tsra () == cp){
00826 ni=sra.give_ni ();
00827 err=sra.give_err ();
00828
00829 Mm->cutting_plane (ipp,im,ido,gamma,epsn,epsp,q,ni,err);
00830 }
00831 else{
00832 fprintf (stderr,"\n\n wrong type of stress return algorithm is required in nlstresses (file %s, line %d).\n",__FILE__,__LINE__);
00833 abort ();
00834 }
00835
00836
00837 for (i=0;i<n;i++){
00838 Mm->ip[ipp].other[ido+i]=epsp[i];
00839 }
00840 Mm->ip[ipp].other[ido+n]=gamma;
00841 Mm->ip[ipp].other[ido+n+1]=q[0];
00842 }
00843
00844
00845
00846
00847
00848
00849 void chen::updateval (long ipp, long im, long ido)
00850 {
00851 long i,n = Mm->givencompeqother(ipp, im);
00852
00853 for (i=0;i<n;i++){
00854 Mm->ip[ipp].eqother[ido+i]=Mm->ip[ipp].other[ido+i];
00855 }
00856 }
00857
00858
00859
00860
00861
00862
00863 void chen::changeparam (atsel &atm,vector &val)
00864 {
00865 long i;
00866
00867 for (i=0;i<atm.num;i++){
00868 switch (atm.atrib[i]){
00869 case 0:{
00870 fyc=val[i];
00871 break;
00872 }
00873 case 1:{
00874 fyt=val[i];
00875 break;
00876 }
00877 case 2:{
00878 fybc=val[i];
00879 break;
00880 }
00881 case 3:{
00882 fc=val[i];
00883 break;
00884 }
00885 case 4:{
00886 ft=val[i];
00887 break;
00888 }
00889 case 5:{
00890 fbc=val[i];
00891 break;
00892 }
00893 default:{
00894 fprintf (stderr,"\n\n wrong number of atribute in function changeparam (file %s, line %d).\n",__FILE__,__LINE__);
00895 }
00896 }
00897 }
00898 }
00899
00900 double chen::give_consparam (long ipp,long ido)
00901 {
00902 long ncompstr;
00903 double gamma;
00904
00905 ncompstr=Mm->ip[ipp].ncompstr;
00906 gamma = Mm->ip[ipp].eqother[ido+ncompstr];
00907
00908 return gamma;
00909 }
00910
00911 long chen::give_num_interparam ()
00912 {
00913 return 1;
00914 }
00915
00916 void chen::give_interparam (long ipp,long ido,vector &q)
00917 {
00918 long ncompstr=Mm->ip[ipp].ncompstr;
00919
00920 q[0]=Mm->ip[ipp].eqother[ido+ncompstr+1];
00921 }