00001 #include <math.h>
00002 #include "glasgmech.h"
00003 #include "vector.h"
00004 #include "matrix.h"
00005 #include "vecttens.h"
00006 #include "tensor.h"
00007 #include "global.h"
00008 #include "intpoints.h"
00009 #include "elastisomat.h"
00010
00011
00012 #define nijac 20
00013 #define limit 1.0e-8
00014
00015 glasgmech::glasgmech (void)
00016 {
00017 tkappa0=0.0;
00018 st=0.0;
00019 k=0.0;
00020 gf0=0.0;
00021 lc=0.0;
00022 a=0.0;
00023 b=0.0;
00024 c=0.0;
00025 ft = paramf_type(0);
00026 }
00027
00028 glasgmech::~glasgmech (void)
00029 {
00030
00031 }
00032
00033
00034
00035
00036 void glasgmech::read (XFILE *in)
00037 {
00038 xfscanf (in,"%lf %lf %lf %lf %m", &tkappa0, &st, &gf0, &lc, ¶mf_type_kwdset, (int *)(&ft));
00039 xfscanf(in, "%lf", &k);
00040 xfscanf (in,"%lf %lf %lf",&a,&b,&c);
00041 }
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051 double glasgmech::betacoeff (long ipp,double t)
00052 {
00053 double tstar,tbar,beta;
00054
00055 tstar=(470.0-20.0)/100.0;
00056 tbar=(t-20.0)/100.0;
00057
00058 if (0.0 > tbar){
00059 fprintf (stderr,"\n\n Warning: normalized temperature t_bar is less than 0.0 at integration point %ld",ipp);
00060 return 0.0;
00061 }
00062 if (0.0<=tbar && tbar<=tstar){
00063 beta=(2.0*a*tbar+b)*0.01;
00064 }
00065 else{
00066 beta=(2.0*c*(tbar-tstar)+2.0*a*tstar+b)*0.01;
00067 }
00068
00069 return beta;
00070 }
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081 void glasgmech::damfuncpar(long ipp, vector &eps, vector &kappa)
00082 {
00083 vector poseps;
00084 vector princeps;
00085 matrix pvect;
00086 matrix epst;
00087 matrix dev;
00088 double i1e, j2e, nu, e, tmp, a, b;
00089 long ncomp, idem;
00090
00091 switch (ft)
00092 {
00093 case vonmises :
00094 ncomp=Mm->ip[ipp].ncompstr;
00095 allocm(3, 3, epst);
00096 allocm(3, 3, dev);
00097 vector_tensor(eps, epst, Mm->ip[ipp].ssst, strain);
00098 i1e = first_invar(epst);
00099 deviator (epst,dev);
00100 j2e = second_invar (dev);
00101 idem = Mm->ip[ipp].gemid();
00102 if (Mm->ip[ipp].tm[idem] != elisomat)
00103 {
00104 fprintf(stderr, "\n\n Invalid type of elastic material is required");
00105 fprintf(stderr, "\n in function glasgmech::damfuncpar, (file %s, line %d)\n", __FILE__, __LINE__);
00106 }
00107 e = Mm->eliso[idem].e;
00108 nu = Mm->eliso[idem].nu;
00109 a = (k-1.0)/(2.0*k)/(1.0-2.0*nu);
00110 b = 3.0/k/(1+nu);
00111 tmp = a*a*i1e*i1e + b*j2e;
00112 kappa[0] = a*i1e+sqrt(tmp);
00113 break;
00114 case normazar :
00115 allocm(3, 3, epst);
00116 allocm(3, 3, pvect);
00117 allocv(3, princeps);
00118 allocv(3,poseps);
00119 fillv(0.0, princeps);
00120 vector_tensor(eps, epst, Mm->ip[ipp].ssst, strain);
00121 princ_val (epst,princeps,pvect,nijac,limit,Mp->zero,3,1);
00122 extractposv (princeps, poseps);
00123 scprd(poseps, poseps, tmp);
00124 kappa[0] = sqrt(tmp);
00125 break;
00126 default :
00127 {
00128 fprintf(stderr, "\n\n Unknown damage parameter function type");
00129 fprintf(stderr, "\n in function glasgmech::damfuncpar, (file %s, line %d)\n", __FILE__, __LINE__);
00130 }
00131 }
00132 return;
00133 }
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145 double glasgmech::damfunction(long ipp,double tempr,double chi,vector &kappa)
00146 {
00147 double omega = 0.0, e0, kappa0, tn, gamma, gf,th;
00148 long idem = Mm->ip[ipp].gemid();
00149
00150 if (Mm->ip[ipp].tm[idem] != elisomat)
00151 {
00152 fprintf(stderr, "\n\nError - ip %ld has wrong type of elastic material combined with", ipp);
00153 fprintf(stderr, "\n scalar damge. The elastic isotropic material is requried");
00154 fprintf(stderr, "\n Function glasgmech::damfunction (file %s, line %d)\n", __FILE__, __LINE__);
00155 return 0.0;
00156 }
00157
00158
00159
00160 e0 = Mm->eliso[Mm->ip[ipp].idm[idem]].e;
00161
00162
00163 tn = (tempr-20.0)/100.0;
00164
00165 th = (kappa[1]-20.0)/100.0;
00166 if (tn>th){
00167 th=tn;
00168 }
00169
00170 if ((th < 0.0) || (th > 7.9)){
00171 fprintf(stderr, "\n\nWarning - largest value of normalised temperature\n");
00172 fprintf(stderr, " is out of suggested range 0.0-7.9 (ipp %ld)\n",ipp);
00173 }
00174
00175
00176 kappa0 = st/e0*(1.0 - 0.016*th*th)/(1.0 - 0.1*th)/(1.0 - 0.1*th);
00177
00178
00179 gf = gf0*(1.0+0.39*th-0.07*th*th);
00180
00181
00182 gamma = (1.0-chi)*st*lc/gf;
00183
00184 if (kappa[0] < kappa0)
00185
00186 return 0.0;
00187
00188
00189 omega = 1.0-(kappa0/kappa[0])*exp(-gamma*(kappa[0]-kappa0));
00190
00191 return omega;
00192 }
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203 double glasgmech::domegadkmd (long ipp,double tempr,double chi, vector &kappa)
00204 {
00205 double domega = 0.0, e0, kappa0, tn, gamma, gf,th;
00206 long idem = Mm->ip[ipp].gemid();
00207
00208 if (Mm->ip[ipp].tm[idem] != elisomat)
00209 {
00210 fprintf(stderr, "\n\nError - ip %ld has wrong type of elastic material combined with", ipp);
00211 fprintf(stderr, "\n scalar damge. The elastic isotropic material is requried");
00212 fprintf(stderr, "\n Function glasgmech::damfunction (file %s, line %d)\n", __FILE__, __LINE__);
00213 return 0.0;
00214 }
00215
00216
00217
00218 e0 = Mm->eliso[Mm->ip[ipp].idm[idem]].e;
00219
00220 tn = (tempr-20.0)/100.0;
00221
00222 th = (kappa[1]-20.0)/100.0;
00223 if (tn>th){
00224 th=tn;
00225 }
00226
00227 if ((th < 0.0) || (th > 7.9)){
00228 fprintf(stderr, "\n\nWarning - largest value of normalised temperature\n");
00229 fprintf(stderr, " is out of suggested range 0.0-7.9 (ipp %ld)\n",ipp);
00230 }
00231
00232
00233 kappa0 = st/e0*(1.0 - 0.016*th*th)/(1.0 - 0.1*th)/(1.0 - 0.1*th);
00234
00235
00236 gf = gf0*(1.0+0.39*th-0.07*th*th);
00237
00238
00239 gamma = (1.0-chi)*st*lc/gf;
00240
00241 if (kappa[0] < kappa0)
00242
00243 return 0.0;
00244
00245 domega = -kappa0*exp(-gamma*(kappa[0]-kappa0))/kappa[0];
00246 domega += kappa0*kappa[0]*gamma*exp(-gamma*(kappa[0]-kappa0));
00247 domega /= kappa[0]*kappa[0];
00248
00249 return domega;
00250 }
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262 double glasgmech::domegadt (long ipp,double tempr,double chi, vector &kappa)
00263 {
00264 double domega = 0.0, e0, kappa0, tn, gamma, gf,th, dkmd0dt, dgdt;
00265 long idem = Mm->ip[ipp].gemid();
00266
00267 if (Mm->ip[ipp].tm[idem] != elisomat)
00268 {
00269 fprintf(stderr, "\n\nError - ip %ld has wrong type of elastic material combined with", ipp);
00270 fprintf(stderr, "\n scalar damge. The elastic isotropic material is requried");
00271 fprintf(stderr, "\n Function glasgmech::damfunction (file %s, line %d)\n", __FILE__, __LINE__);
00272 return 0.0;
00273 }
00274
00275
00276
00277 e0 = Mm->eliso[Mm->ip[ipp].idm[idem]].e;
00278
00279 tn = (tempr-20.0)/100.0;
00280
00281 th = (kappa[1]-20.0)/100.0;
00282 if (tn>th){
00283 th=tn;
00284 }
00285
00286 if ((th < 0.0) || (th > 7.9)){
00287 fprintf(stderr, "\n\nWarning - largest value of normalised temperature\n");
00288 fprintf(stderr, " is out of suggested range 0.0-7.9 (ipp %ld)\n",ipp);
00289 }
00290
00291
00292 kappa0 = st/e0*(1.0 - 0.016*th*th)/(1.0 - 0.1*th)/(1.0 - 0.1*th);
00293
00294
00295 gf = gf0*(1.0+0.39*th-0.07*th*th);
00296
00297
00298 gamma = (1.0-chi)*st*lc/gf;
00299
00300 if (kappa[0] < kappa0)
00301
00302 return 0.0;
00303
00304
00305 dkmd0dt = st/e0/100.0*(-0.032*th*(1.0-0.1*th)*(1.0-0.1*th)+(1.0-0.016*th*th)*2.0*(1.0-0.1*th)*0.1);
00306 dkmd0dt /= (1.0-0.1*th)*(1.0-0.1*th)*(1.0-0.1*th)*(1.0-0.1*th);
00307
00308 dgdt = -(1-chi)*st*lc/(100.0*gf*gf)*(0.39-0.14*th)*gf0;
00309
00310 domega = dkmd0dt*exp(-gamma*(kappa[0]-kappa0))/kappa[0];
00311 domega += kappa0/kappa[0]*exp(-gamma*(kappa[0]-kappa0))*(-(dgdt*(kappa[0]-kappa0)-gamma*dkmd0dt));
00312 domega = -domega;
00313
00314 return domega;
00315 }
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325 void glasgmech::matstiff (matrix &d,long ipp,long ido)
00326 {
00327 double dp;
00328
00329 switch (Mp->nlman->stmat)
00330 {
00331 case initial_stiff :
00332 Mm->elmatstiff (d,ipp);
00333 break;
00334 case tangent_stiff :
00335 Mm->elmatstiff (d,ipp);
00336 dp=Mm->ip[ipp].eqother[ido+1];
00337 if (dp > 0.999999)
00338 dp = 0.999999;
00339 cmulm (1.0-dp,d);
00340 break;
00341 default :
00342 fprintf(stderr, "\n\nError - unknown type of stifness matrix");
00343 fprintf(stderr, "\n in function glasgmech::matstiff (%s, line %d)\n", __FILE__, __LINE__);
00344 }
00345 }
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355 void glasgmech::compute_thermdilat (double t,double dt,matrix &eps)
00356 {
00357 double auxt,alpha;
00358
00359 auxt = (t-20.0)/100.0;
00360 if (auxt>=0.0 && auxt<=6.0){
00361 alpha = 6.0e-5/(7.0-auxt);
00362 }
00363 else{ alpha=0.0; }
00364
00365 fillm (0.0,eps);
00366
00367 eps[0][0]=alpha*dt;
00368 eps[1][1]=alpha*dt;
00369 eps[2][2]=alpha*dt;
00370
00371 }
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384 double glasgmech::thermdamfunction (long ipp,double tempr,vector &kappa)
00385 {
00386 double chi = 0.0;
00387
00388 if (tempr > kappa[0])
00389 kappa[0] = tempr;
00390 if (kappa[0] > tkappa0)
00391 chi = 20.0*(kappa[0] - tkappa0)*(1.0-5.0*(kappa[0]-tkappa0));
00392 else
00393 chi = 0.0;
00394
00395 return chi;
00396 }
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409 double glasgmech::dtdfdt (long ipp,double tempr,vector &kappa)
00410 {
00411 double htd;
00412
00413 if (tempr > kappa[0])
00414 kappa[0] = tempr;
00415
00416 if (kappa[0] > tkappa0)
00417 htd = 20.0*(1.0-5.0*(kappa[0]-tkappa0)) - 5.0*20.0*(kappa[0] - tkappa0);
00418 else
00419 htd = 0.0;
00420
00421 return htd;
00422 }
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433 void glasgmech::compute_lits (long ipp,vector &epstm,matrix &sigt,double told,double tnew)
00434 {
00435 long idem;
00436 double dt,i1,beta,nu,fc;
00437 vector princsig(3);
00438 matrix pvect(3,3);
00439 strastrestate ssst;
00440
00441 idem=Mm->ip[ipp].gemid();
00442 if (Mm->ip[ipp].tm[idem] != elisomat)
00443 {
00444 fprintf(stderr, "\n\n Invalid type of elastic material is required");
00445 fprintf(stderr, "\n in function glasgmech::compute_lits, (file %s, line %d)\n", __FILE__, __LINE__);
00446 }
00447 nu = Mm->eliso[idem].nu;
00448
00449 fillv (0.0,epstm);
00450
00451
00452 dt=tnew-told;
00453 if (dt<=0.0){
00454 return;
00455 }
00456 else{
00457 princ_val (sigt,princsig,pvect,nijac,limit,Mp->zero,3,1);
00458 extractnegv (princsig,princsig);
00459 fillm (0.0,sigt);
00460 sigt[0][0]=princsig[0];
00461 sigt[1][1]=princsig[1];
00462 sigt[2][2]=princsig[2];
00463 lgmatrixtransf(sigt, pvect);
00464
00465 i1=first_invar (sigt);
00466
00467 beta = betacoeff (ipp,tnew);
00468
00469 fc = st*k;
00470
00471 cmulm ((1.0+nu),sigt);
00472 sigt[0][0]-=nu*i1;
00473 sigt[1][1]-=nu*i1;
00474 sigt[2][2]-=nu*i1;
00475 cmulm (beta/fc*dt,sigt);
00476
00477 ssst=Mm->ip[ipp].ssst;
00478 tensor_vector (epstm,sigt,ssst,strain);
00479 }
00480 }
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491 void glasgmech::depseqdepsel(long ipp, vector &eps, vector &deeqdeel)
00492 {
00493 double q1, q2, q3, e, nu, i1e, j2e, tmp;
00494 double r, af, bf;
00495 long i;
00496 long idem = Mm->ip[ipp].gemid();
00497 long ncomp = Mm->ip[ipp].ncompstr;
00498 vector pv(ncomp), tv(ncomp);
00499 matrix pm(ncomp,ncomp);
00500 matrix epst(3,3), dev(3,3);
00501
00502 fillv(0.0, pv);
00503 fillm(0.0, pm);
00504 fillv(0.0,deeqdeel);
00505
00506 vector_tensor(eps, epst, Mm->ip[ipp].ssst, strain);
00507 i1e = first_invar(epst);
00508 deviator (epst,dev);
00509 j2e = second_invar (dev);
00510 idem = Mm->ip[ipp].gemid();
00511 if (Mm->ip[ipp].tm[idem] != elisomat)
00512 {
00513 fprintf(stderr, "\n\n Invalid type of elastic material is required");
00514 fprintf(stderr, "\n in function glasgmech::depseqds, (file %s, line %d)\n", __FILE__, __LINE__);
00515 }
00516 e = Mm->eliso[idem].e;
00517 nu = Mm->eliso[idem].nu;
00518 af = (k-1.0)/(2.0*k)/(1.0-2.0*nu);
00519 bf = 3.0/k/(1+nu);
00520 r = af*af*i1e*i1e + bf*j2e;
00521 r = sqrt(r);
00522 if (r > 0.0)
00523 r = 0.5 / r;
00524 else
00525 r = 0.0;
00526 q3 = nu/(1.0-nu);
00527 q1 = 1.0 + q3 + q3*q3;
00528 q2 = 1.0 - 2.0*q3 - 2.0*q3*q3;
00529 switch (Mm->ip[ipp].ssst)
00530 {
00531 case planestrain:
00532 case axisymm:
00533 pv[0] = 1.0-q3;
00534 pv[1] = 1.0-q3;
00535 pv[2] = 0.0;
00536 pv[3] = 0.0;
00537 pm[0][0] = 2.0/3.0*q1;
00538 pm[1][1] = pm[0][0];
00539 pm[2][2] = 0.0;
00540 pm[3][3] = 2.0;
00541 pm[0][1] = -q2/3.0;
00542 pm[1][0] = pm[0][1];
00543 break;
00544 case spacestress:
00545 for (i = 0; i < 3; i++)
00546 {
00547 pv[i] = 1.0;
00548 pm[i][i] = 2.0/3.0;
00549 pm[i+3][i+3] = 2.0;
00550 }
00551 pm[0][1] = -q2/3.0;
00552 pm[0][2] = pm[0][1];
00553 pm[1][2] = pm[0][1];
00554 pm[2][0] = pm[0][1];
00555 pm[2][1] = pm[0][1];
00556 break;
00557 default:
00558 fprintf(stderr, "\n\n Invalid stress/strain state is required in function glasgmech::depseqds\n");
00559 fprintf(stderr, " in integeration point %ld, (file %s, line %d)\n", ipp, __FILE__, __LINE__);
00560 }
00561 for(i = ncomp / 2; i < ncomp; i++)
00562 eps[i] *= 0.5;
00563 tmp = af+2.0*af*af*i1e*r;
00564 cmulv(tmp, pv);
00565 mxv(pm, pv, tv);
00566 cmulv(bf*r, tv);
00567 addv(tv, pv, deeqdeel);
00568 for(i = ncomp / 2; i < ncomp; i++)
00569 {
00570 deeqdeel[i] *= 0.5;
00571 eps[i] *= 2.0;
00572 }
00573 }
00574
00575
00576
00577
00578
00579
00580
00581
00582
00583
00584
00585
00586 void glasgmech::nlstresses (long ipp,long im,long ido)
00587 {
00588 long i,ncompstr;
00589
00590
00591 ncompstr = Mm->ip[ipp].ncompstr;
00592
00593 double newt,oldt,dt,chi,htd,omega,domega, omegao, kappao, depseq;
00594
00595 vector epsn(ncompstr),epso(ncompstr),epse(ncompstr),epsft(ncompstr),epstm(ncompstr);
00596 vector epsc(ncompstr),epsirr(ncompstr),epstirr(ncompstr), deeqdeel(ncompstr), deps(ncompstr);
00597 vector sig(ncompstr),sigtrial(ncompstr);
00598 vector kappa(2),tkappa(1);
00599 vector deftdt(ncompstr), s(ncompstr);
00600
00601 matrix eps(3,3),d(ncompstr,ncompstr);
00602 matrix sigt(3,3);
00603
00604 i = 0;
00605
00606 if (i==1){
00607
00608
00609
00610
00611
00612 for (i=0;i<ncompstr;i++){
00613 epsn[i] = Mm->ip[ipp].strain[i];
00614 }
00615
00616 for (i=0;i<ncompstr;i++){
00617 epso[i] = Mm->ip[ipp].eqother[ncompstr+i];
00618 }
00619
00620
00621
00622
00623
00624
00625
00626 newt = Mm->givenonmechq(temperature, ipp);
00627
00628 oldt = Mm->ip[ipp].eqother[4*ncompstr];
00629
00630 dt = newt-oldt;
00631
00632
00633 compute_thermdilat (newt,dt,eps);
00634
00635 tensor_vector (epsft,eps,Mm->ip[ipp].ssst,strain);
00636
00637
00638
00639
00640 vector_tensor(sig, sigt, Mm->ip[ipp].ssst, stress);
00641 compute_lits (ipp, epstm, sigt, oldt, newt);
00642
00643
00644
00645
00646
00647
00648
00649 tkappa[0]=Mm->ip[ipp].eqother[4*ncompstr+3];
00650
00651
00652 chi = thermdamfunction (ipp,newt,tkappa);
00653
00654
00655
00656
00657
00658
00659 damfuncpar(ipp,epsn,kappa);
00660
00661
00662 kappao = Mm->ip[ipp].eqother[4*ncompstr+1];
00663
00664 omegao = Mm->ip[ipp].eqother[4*ncompstr+2];
00665
00666 kappa[1] = Mm->ip[ipp].eqother[4*ncompstr];
00667
00668 omega = damfunction(ipp,newt,chi,kappa);
00669 if (omega < kappao)
00670 {
00671 omega = omegao;
00672 domega = 0.0;
00673 }
00674 else
00675 {
00676
00677
00678
00679 depseqdepsel(ipp, epsn, deeqdeel);
00680 subv (epsn, epso, deps);
00681 subv (deps, epsft, deps);
00682 subv (deps, epstm, deps);
00683 scprd(deps, deeqdeel, depseq);
00684
00685 domega = domegadt(ipp, newt, chi, kappa)*depseq;
00686 domega += domegadkmd(ipp, newt, chi, kappa)*(newt-oldt);
00687 }
00688
00689
00690
00691
00692
00693
00694
00695
00696
00697
00698
00699 for (i=0;i<ncompstr;i++){
00700 epsirr[i]=epsft[i]+epsc[i]+epstm[i];
00701 }
00702
00703
00704
00705
00706 for (i=0;i<ncompstr;i++){
00707 Mm->ip[ipp].eqother[3*ncompstr+i]=epsirr[i];
00708 }
00709 Mm->ip[ipp].eqother[4*ncompstr+1]=kappa[0];
00710 Mm->ip[ipp].eqother[4*ncompstr+2]=omega;
00711 Mm->ip[ipp].eqother[4*ncompstr+3]=tkappa[0];
00712 Mm->ip[ipp].eqother[4*ncompstr+4]=chi;
00713 Mm->ip[ipp].eqother[4*ncompstr+5]=domega;
00714
00715
00716 Mm->elmatstiff(d,ipp);
00717
00718
00719 mxv (d,epsirr,sig);
00720
00721 Mm->storestress (0,ipp,sig);
00722 }
00723
00724
00725
00726 if (i==2){
00727
00728 for (i=0;i<ncompstr;i++){
00729
00730 epsn[i] = Mm->ip[ipp].strain[i];
00731
00732 epso[i] = Mm->ip[ipp].eqother[ido+ncompstr+i];
00733
00734 epstirr[i] = Mm->ip[ipp].eqother[ido+2*ncompstr+i];
00735
00736 epsirr[i] = Mm->ip[ipp].eqother[ido+3*ncompstr+i];
00737
00738
00739 epse[i]=epsn[i]-epsirr[i];
00740 }
00741 omega =Mm->ip[ipp].eqother[4*ncompstr+2];
00742 chi =Mm->ip[ipp].eqother[4*ncompstr+4];
00743 domega=Mm->ip[ipp].eqother[4*ncompstr+5];
00744
00745
00746 newt = 20.0;
00747
00748 oldt = Mm->ip[ipp].eqother[4*ncompstr];
00749
00750 dt = newt-oldt;
00751
00752 tkappa[0]=Mm->ip[ipp].eqother[4*ncompstr+3];
00753
00754 htd = dtdfdt (ipp,newt,tkappa);
00755
00756
00757
00758 Mm->matstiff(d,ipp);
00759
00760 mxv (d,epse,sigtrial);
00761
00762
00763 subv (epsn,epso,epso);
00764
00765 subv (epso,epsirr,epse);
00766
00767 cmulm ((1.0-omega)*(1.0-chi),d);
00768
00769 mxv (d,epse,sig);
00770
00771 cmulv (htd*dt*(1.0-omega)+domega*(1.0-chi),sigtrial);
00772
00773 subv (sig,sigtrial,sig);
00774
00775
00776
00777 Mm->ip[ipp].eqother[4*ncompstr+0]=newt;
00778 for (i=0;i<ncompstr;i++){
00779
00780 Mm->ip[ipp].eqother[ido+i] += sig[i];
00781
00782 Mm->ip[ipp].eqother[ido+ncompstr+i] = epsn[i];
00783
00784 Mm->ip[ipp].eqother[ido+2*ncompstr+i] += epsirr[i];
00785
00786 Mm->ip[ipp].stress[i]=sig[i];
00787 }
00788
00789 }
00790 }