00001 #include <math.h>
00002 #include <stdio.h>
00003 #include <stdlib.h>
00004 #include "drprag.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
00013
00014 #define nijac 20
00015 #define limit 1.0e-8
00016
00017
00018
00019
00020
00021
00022
00023 drprag::drprag (void)
00024 {
00025 phi=0.0; c=0.0; psi=0.0;
00026 theta = 0.0; clim = 0.0;
00027 }
00028
00029
00030
00031
00032
00033
00034 drprag::~drprag (void)
00035 {
00036
00037 }
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050 void drprag::read (XFILE *in)
00051 {
00052 xfscanf (in,"%k%lf%k%lf%k%lf%k%lf%k%lf","phi",&phi,"coh",&c,"psi",&psi,"theta",&theta,"clim",&clim);
00053 sra.read (in);
00054
00055 alpha=6.0*sin(phi)/(sqrt(3.0)*(3.0-sin(phi)));
00056
00057 alpha1=6.0*sin(psi)/(sqrt(3.0)*(3.0-sin(psi)));
00058
00059 beta=6.0*cos(phi)/(sqrt(3.0)*(3.0-sin(phi)));
00060
00061 }
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073 double drprag::yieldfunction (matrix &sig, vector &q)
00074 {
00075 double i1s,j2s,sigm_max,h,f;
00076 matrix dev(3,3);
00077
00078 i1s = first_invar (sig)/3.0;
00079
00080 deviator (sig,dev);
00081
00082 j2s = second_invar(dev);
00083 j2s = sqrt(0.5)*tensornorm(dev);
00084
00085 sigm_max = cohesion(q)/tan(phi);
00086
00087 if (i1s > sigm_max)
00088 h = i1s/sqrt(3.0)/alpha - sigm_max/sqrt(3.0)/alpha;
00089 else
00090 h = 0.0;
00091 if (j2s >= h)
00092 {
00093 f = j2s + alpha*i1s - beta*cohesion(q);
00094 }
00095 else
00096 {
00097 f = i1s - sigm_max;
00098 }
00099
00100 return f;
00101 }
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118 void drprag::deryieldfsigma (matrix &sig, matrix &dfds, vector &q)
00119 {
00120 double h, i1s, j2s, sigm_max;
00121 matrix dev(3,3);
00122
00123 i1s = first_invar(sig)/3.0;
00124 deviator(sig,dev);
00125 j2s = sqrt(0.5)*tensornorm(dev);
00126 sigm_max = cohesion(q)/tan(phi);
00127
00128 if (i1s > sigm_max)
00129 h = i1s/sqrt(3.0)/alpha - sigm_max/sqrt(3.0)/alpha;
00130 else
00131 h = 0.0;
00132 if (j2s >= h)
00133 {
00134 cmulm (0.5/j2s,dev,dfds);
00135
00136 dfds[0][0]+=alpha/3.0;
00137 dfds[1][1]+=alpha/3.0;
00138 dfds[2][2]+=alpha/3.0;
00139 }
00140 else
00141 {
00142 fillm(0.0, dfds);
00143 dfds[0][0]=1.3;
00144 dfds[1][1]=1.3;
00145 dfds[2][2]=1.3;
00146 }
00147 }
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164 void drprag::derpotsigma (matrix &sig,matrix &dgds, vector &q)
00165 {
00166 double h, i1s, j2s, sigm_max;
00167 matrix dev(3,3);
00168
00169
00170 i1s = first_invar(sig)/3.0;
00171 deviator(sig,dev);
00172 j2s = sqrt(0.5)*tensornorm(dev);
00173 sigm_max = c/tan(psi);
00174
00175 if (i1s > sigm_max)
00176 h = i1s/sqrt(3.0)/alpha1 - sigm_max/sqrt(3.0)/alpha1;
00177 else
00178 h = 0.0;
00179
00180 if (j2s >= h)
00181 {
00182 cmulm (0.5/j2s,dev,dgds);
00183 dgds[0][0]+=alpha1/3.0;
00184 dgds[1][1]+=alpha1/3.0;
00185 dgds[2][2]+=alpha1/3.0;
00186 }
00187 else
00188 {
00189 fillm(0.0, dgds);
00190 dgds[0][0]=1.3;
00191 dgds[1][1]=1.3;
00192 dgds[2][2]=1.3;
00193 }
00194 }
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209 double drprag::cohesion(vector &qtr)
00210 {
00211 double tgt = tan(theta);
00212 double cq = c + tgt*qtr[0];
00213 if (((cq < clim) && (c < clim)) ||
00214 ((cq > clim) && (c > clim)))
00215 return cq;
00216 return clim;
00217 }
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233 void drprag::deryieldfq(vector &qtr, vector &dfq)
00234 {
00235 double tgt = tan(theta);
00236 double cq = c + tgt*qtr[0];
00237 if (((cq < clim) && (c < clim)) ||
00238 ((cq > clim) && (c > clim)))
00239 dfq[0] = -beta*tgt;
00240 else
00241 dfq[0] = 0.0;
00242 return;
00243 }
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258 void drprag::der_q_gamma(vector &dqdg)
00259 {
00260 dqdg[0] = sqrt(1 + 2/9*alpha1*alpha1);
00261 }
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274 double drprag::plasmodscalar(vector &qtr)
00275 {
00276 double ret;
00277 vector dfq(qtr.n);
00278 vector dqg(qtr.n);
00279
00280 deryieldfq(qtr, dfq);
00281 der_q_gamma(dqg);
00282 scprd(dfq, dqg, ret);
00283
00284
00285 return -ret;
00286 }
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302 void drprag::updateq(long ipp, vector &epsp, vector &q)
00303 {
00304 matrix epst(3,3);
00305 strastrestate ssst=Mm->ip[ipp].ssst;
00306
00307 vector_tensor (epsp, epst, ssst, strain);
00308 cumulstrain(epst, q[0]);
00309 }
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324 void drprag::matstiff (matrix &d, long ipp,long ido)
00325 {
00326 if (Mp->nlman->stmat==initial_stiff)
00327 {
00328
00329 Mm->elmatstiff (d,ipp);
00330 }
00331 if (Mp->nlman->stmat==tangent_stiff)
00332 {
00333
00334 matrix ad(d.m,d.n);
00335 Mm->elmatstiff (ad,ipp);
00336 tangentstiff (ad,d,ipp,ido);
00337 }
00338 }
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354 void drprag::tangentstiff (matrix &d,matrix &td,long ipp,long ido)
00355 {
00356 long ncomp=Mm->ip[ipp].ncompstr;
00357 double denom,gamma;
00358 vector str,av(d.m),q(1);
00359 matrix sig(3,3),am(d.m,d.n);
00360
00361 gamma=Mm->ip[ipp].eqother[ido+ncomp];
00362 if (gamma<1.0e-10){
00363 copym (d,td);
00364 }
00365 else{
00366
00367 allocv (ncomp,str);
00368
00369 Mm->givestress (0,ipp,str);
00370 vector_tensor (str,sig,Mm->ip[ipp].ssst,stress);
00371 deryieldfsigma (sig,sig,q);
00372 tensor_vector (str,sig,Mm->ip[ipp].ssst,stress);
00373
00374 if ((Mm->ip[ipp].ssst==planestress) || (Mm->ip[ipp].ssst==planestrain)){
00375 vector auxstr(3);
00376 auxstr[0]=str[0];auxstr[1]=str[1];auxstr[2]=str[2];
00377 destrv (str);
00378 allocv (d.m,str);
00379 str[0]=auxstr[0];str[1]=auxstr[1];str[2]=auxstr[2];
00380 }
00381
00382 mxv (d,str,av);
00383 scprd (av,str,denom);
00384
00385
00386 q[0] = Mm->ip[ipp].eqother[ido+ncomp+1];
00387 denom+= plasmodscalar(q);
00388
00389 if (fabs(denom)<1.0e-10){
00390 copym (d,td);
00391 }
00392 else{
00393 vxv (str,str,am);
00394 mxm (d,am,td);
00395 mxm (td,d,am);
00396
00397 cmulm (1.0/denom,am);
00398
00399 subm (d,am,td);
00400 }
00401 }
00402
00403 }
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421 void drprag::nlstresses (long ipp, long im, long ido)
00422 {
00423 long i,ni,n=Mm->ip[ipp].ncompstr;
00424 double gamma,err;
00425 vector epsn(n),epsp(n),q(1);
00426
00427
00428 for (i=0; i<n; i++)
00429 {
00430 epsn[i]=Mm->ip[ipp].strain[i];
00431 epsp[i]=Mm->ip[ipp].eqother[ido+i];
00432 }
00433 gamma=Mm->ip[ipp].eqother[ido+n];
00434 q[0] = Mm->ip[ipp].eqother[ido+n+1];
00435
00436
00437 if (sra.give_tsra () == cp){
00438 ni=sra.give_ni ();
00439 err=sra.give_err ();
00440
00441 Mm->cutting_plane (ipp,im,ido,gamma,epsn,epsp,q,ni,err);
00442 }
00443 else{
00444 print_err("wrong type of stress return algorithm is required", __FILE__, __LINE__, __func__);
00445 abort ();
00446 }
00447
00448 double i1s,j2s;
00449 matrix dev(3,3);
00450 matrix sigt(3,3);
00451 vector sig(n);
00452 for (i=0;i<n;i++)
00453 sig[i] = Mm->ip[ipp].stress[i];
00454 vector_tensor (sig,sigt,Mm->ip[ipp].ssst,stress);
00455 i1s = first_invar (sigt)/3.0;
00456 deviator (sigt,dev);
00457 j2s = second_invar (dev);
00458
00459
00460 for (i=0;i<n;i++){
00461 Mm->ip[ipp].other[ido+i]=epsp[i];
00462 }
00463 Mm->ip[ipp].other[ido+n]=gamma;
00464 Mm->ip[ipp].other[ido+n+1]=q[0];
00465 Mm->ip[ipp].other[ido+n+2]=i1s;
00466 Mm->ip[ipp].other[ido+n+3]=j2s;
00467 }
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485 void drprag::nonloc_nlstresses (long ipp, long im, long ido)
00486 {
00487 long i,ni,n=Mm->ip[ipp].ncompstr;
00488 double gamma,err;
00489 vector epsn(n),epsp(n),q(1);
00490
00491
00492 for (i=0; i<n; i++)
00493 {
00494 epsn[i]=Mm->ip[ipp].strain[i];
00495 epsp[i]=Mm->ip[ipp].nonloc[i];
00496 }
00497 gamma=Mm->ip[ipp].eqother[ido+n];
00498 q[0] = Mm->ip[ipp].eqother[ido+n+1];
00499
00500
00501 if (sra.give_tsra () == cp){
00502 ni=sra.give_ni ();
00503 err=sra.give_err ();
00504
00505 Mm->cutting_plane (ipp,im,ido,gamma,epsn,epsp,q,ni,err);
00506 }
00507 else{
00508 print_err("wrong type of stress return algorithm is required", __FILE__, __LINE__, __func__);
00509 abort ();
00510 }
00511
00512 double i1s,j2s;
00513 matrix dev(3,3);
00514 matrix sigt(3,3);
00515 vector sig(n);
00516 for (i=0;i<n;i++)
00517 sig[i] = Mm->ip[ipp].stress[i];
00518 vector_tensor (sig,sigt,Mm->ip[ipp].ssst,stress);
00519 i1s = first_invar (sigt)/3.0;
00520 deviator (sigt,dev);
00521 j2s = second_invar (dev);
00522
00523 for (i=0;i<n;i++){
00524 Mm->ip[ipp].other[ido+i]=epsp[i];
00525 }
00526 Mm->ip[ipp].other[ido+n]=gamma;
00527 Mm->ip[ipp].other[ido+n+1]=q[0];
00528 Mm->ip[ipp].other[ido+n+2]=i1s;
00529 Mm->ip[ipp].other[ido+n+3]=j2s;
00530 }
00531
00532
00533
00534
00535
00536
00537
00538
00539
00540
00541
00542
00543
00544
00545
00546 void drprag::updateval (long ipp, long im, long ido)
00547 {
00548 long i,n = Mm->givencompeqother(ipp,im);
00549
00550 for (i=0;i<n;i++){
00551 Mm->ip[ipp].eqother[ido+i]=Mm->ip[ipp].other[ido+i];
00552 }
00553 }
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568 void drprag::giveirrstrains (long ipp, long ido, vector &epsp)
00569 {
00570 long i;
00571 for (i=0;i<epsp.n;i++)
00572 epsp[i] = Mm->ip[ipp].eqother[ido+i];
00573 }
00574
00575
00576
00577
00578
00579
00580
00581
00582
00583
00584
00585
00586
00587
00588 double drprag::give_consparam (long ipp, long ido)
00589 {
00590 long ncompstr;
00591 double gamma;
00592
00593 ncompstr=Mm->ip[ipp].ncompstr;
00594 gamma = Mm->ip[ipp].eqother[ido+ncompstr];
00595
00596 return gamma;
00597 }
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613 void drprag::changeparam (atsel &atm,vector &val)
00614 {
00615 long i;
00616
00617 for (i=0;i<atm.num;i++){
00618 switch (atm.atrib[i]){
00619 case 0:{
00620 phi=val[i];
00621 break;
00622 }
00623 case 1:{
00624 c=val[i];
00625 break;
00626 }
00627 case 2:{
00628 psi=val[i];
00629 break;
00630 }
00631 case 3:{
00632 theta=val[i];
00633 break;
00634 }
00635 case 4:{
00636 clim=val[i];
00637 break;
00638 }
00639 default:{
00640 fprintf (stderr,"\n\n wrong number of atribute in function changeparam (file %s, line %d).\n",__FILE__,__LINE__);
00641 }
00642 }
00643 }
00644
00645 alpha=6.0*sin(phi)/(3.0-sin(phi));
00646 alpha1=6.0*sin(psi)/(3.0-sin(psi));
00647
00648 beta=6.0*cos(phi)/(3.0-sin(phi));
00649 }