00001 #include <math.h>
00002 #include <stdio.h>
00003 #include <stdlib.h>
00004 #include "mohrc.h"
00005 #include "global.h"
00006 #include "matrix.h"
00007 #include "vector.h"
00008 #include "intpoints.h"
00009 #include "tensor.h"
00010 #include "vecttens.h"
00011 #include "elastisomat.h"
00012
00013
00014
00015 #define nijac 20
00016 #define limit 1.0e-8
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 mohrcoulomb::mohrcoulomb (void)
00027 {
00028 phi=0.0; c=0.0; psi=0.0;
00029 }
00030
00031
00032
00033
00034
00035
00036
00037
00038 mohrcoulomb::~mohrcoulomb (void)
00039 {
00040
00041 }
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055 void mohrcoulomb::read (XFILE *in)
00056 {
00057 xfscanf (in,"%lf %lf %lf",&phi,&c,&psi);
00058 sra.read (in);
00059 }
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072 double mohrcoulomb::yieldfunction (vector &psig)
00073 {
00074 double f,k1,k2, s1, s2, sigma, tau, maxsigt;
00075
00076 s1 = psig[0];
00077 s2 = psig[2];
00078 sigma = (s1 + s2) / 2.0;
00079 tau = -(s1 - s2) / (2.0 * cos(phi));
00080 maxsigt = c / tan(phi);
00081
00082 k1 = -1.0 + sin(phi); k2 = 1.0 + sin(phi);
00083 f = k1*s1 + k2*s2 - 2.0*c*cos(phi);
00084 return f;
00085 }
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099 void mohrcoulomb::deryieldfsigma (vector &dfds)
00100 {
00101 dfds[0] = -1.0 + sin(phi);
00102 dfds[1] = 0.0;
00103 dfds[2] = 1.0 + sin(phi);
00104 return;
00105 }
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119 void mohrcoulomb::derplaspotsigma (vector &dgds)
00120 {
00121 dgds[0] = -1.0 + sin(psi);
00122 dgds[1] = 0.0;
00123 dgds[2] = 1.0 + sin(psi);
00124 return;
00125 }
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141 void mohrcoulomb::matstiff (matrix &d, long ipp,long ido)
00142 {
00143 matrix ad(d.m,d.n);
00144 switch (Mp->nlman->stmat)
00145 {
00146 case initial_stiff:
00147
00148 Mm->elmatstiff (d,ipp);
00149 break;
00150 case tangent_stiff:
00151 Mm->elmatstiff (ad,ipp);
00152 tangentstiff (ad,d,ipp,ido);
00153 break;
00154 default:
00155 print_err("unknown type of the stiffness matrix is required", __FILE__, __LINE__, __func__);
00156 break;
00157 }
00158 }
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175 void mohrcoulomb::tangentstiff (matrix &d,matrix &td,long ipp,long ido)
00176 {
00177 long ncomp=Mm->ip[ipp].ncompstr;
00178 double denom,gamma;
00179 vector str,av(d.m),q(1),psig(3), b(3);
00180 matrix sig(3,3),am(d.m,d.n),pvect(3,3);
00181
00182 gamma=Mm->ip[ipp].eqother[ido+ncomp];
00183 if (gamma<1.0e-10){
00184 copym (d,td);
00185 }
00186 else{
00187
00188 allocv (ncomp,str);
00189
00190 Mm->givestress (0,ipp,str);
00191 vector_tensor (str,sig,Mm->ip[ipp].ssst,stress);
00192 princ_val (sig, psig, pvect, nijac, limit, Mp->zero, 3,1);
00193 deryieldfsigma (psig);
00194 mtxv(pvect, psig, b);
00195 if ((Mm->ip[ipp].ssst==planestress) || (Mm->ip[ipp].ssst==planestrain)){
00196 vector auxstr(3);
00197 auxstr[0]=str[0];auxstr[1]=str[1];auxstr[2]=str[2];
00198 destrv (str);
00199 allocv (d.m,str);
00200 str[0]=auxstr[0];str[1]=auxstr[1];str[2]=auxstr[2];
00201 }
00202
00203 mxv (d,b,av);
00204 scprd (av,b,denom);
00205
00206
00207
00208
00209
00210 if (fabs(denom)<1.0e-10){
00211 copym (d,td);
00212 }
00213 else{
00214 vxv (b,b,am);
00215 mxm (d,am,td);
00216 mxm (td,d,am);
00217
00218 cmulm (1.0/denom,am);
00219
00220 subm (d,am,td);
00221 }
00222 }
00223
00224 }
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241 void mohrcoulomb::pelmatstiff (long ipp, matrix &d, double e, double nu)
00242 {
00243 double b;
00244
00245 switch (Mm->ip[ipp].ssst)
00246 {
00247 case planestress:
00248 b = e/(1.0-nu*nu);
00249 d[0][0] = b; d[0][1] = b*nu; d[0][2] = b*nu;
00250 d[1][0] = b*nu; d[1][1] = b; d[1][2] = b*nu;
00251 d[2][0] = b*nu; d[2][1] = b*nu; d[2][2] = b;
00252 break;
00253 case planestrain:
00254 case spacestress:
00255 b = e /((1.0 + nu)*(1.0 - 2.0*nu));
00256 d[0][0] = d[1][1] = d[2][2] = 1.0 - nu;
00257 d[0][1] = d[0][2] = d[1][0] = d[1][2] = d[2][0] = d[2][1] = nu;
00258 cmulm(b, d);
00259 break;
00260 default:
00261 print_err("unknown type of stress/strain state is required", __FILE__, __LINE__, __func__);
00262 break;
00263 }
00264 }
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280 void mohrcoulomb::nlstresses (long ipp, long ido)
00281
00282 {
00283 long i,ni,n=Mm->ip[ipp].ncompstr,mu;
00284 double gamma,err;
00285 vector epsn(n),epsp(n),q(0);
00286
00287
00288 for (i=0; i<n; i++)
00289 {
00290 epsn[i]=Mm->ip[ipp].strain[i];
00291 epsp[i]=Mm->ip[ipp].eqother[ido+i];
00292 }
00293 gamma=Mm->ip[ipp].eqother[ido+n];
00294
00295
00296 if (sra.give_tsra () == cp){
00297 ni=sra.give_ni ();
00298 err=sra.give_err ();
00299
00300 mu = cutting_plane (ipp,gamma,epsn,epsp,q,ni,err);
00301 }
00302 else{
00303 print_err("wrong type of stress return algorithm is required", __FILE__, __LINE__, __func__);
00304 abort ();
00305 }
00306
00307
00308 if (mu > 0)
00309 {
00310
00311 gamma=Mm->ip[ipp].eqother[ido+n];
00312
00313 ni=sra.give_ni ();
00314 err=sra.give_err ();
00315 mc_msurf_cp(ipp, gamma, epsn, epsp, q, mu,ni,err);
00316 }
00317
00318
00319 for (i=0;i<n;i++){
00320 Mm->ip[ipp].other[ido+i]=epsp[i];
00321 }
00322 Mm->ip[ipp].other[ido+n]=gamma;
00323 }
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340 void mohrcoulomb::nonloc_nlstresses (long ipp,long ido)
00341 {
00342 long i,ni,n=Mm->ip[ipp].ncompstr,mu;
00343 double gamma,err;
00344 vector epsn(n),epsp(n),q(0);
00345
00346
00347 for (i=0; i<n; i++)
00348 {
00349 epsn[i]=Mm->ip[ipp].strain[i];
00350 epsp[i]=Mm->ip[ipp].nonloc[i];
00351 }
00352 gamma=Mm->ip[ipp].eqother[ido+n];
00353
00354
00355 if (sra.give_tsra () == cp){
00356 ni=sra.give_ni ();
00357 err=sra.give_err ();
00358
00359 mu = cutting_plane (ipp,gamma,epsn,epsp,q,ni,err);
00360 }
00361 else{
00362 print_err("wrong type of stress return algorithm is required", __FILE__, __LINE__, __func__);
00363 abort ();
00364 }
00365
00366
00367 if (mu > 0)
00368
00369
00370 ni=sra.give_ni ();
00371 err=sra.give_err ();
00372 mc_msurf_cp(ipp, gamma, epsn, epsp, q, mu,ni,err);
00373
00374
00375 for (i=0;i<n;i++){
00376 Mm->ip[ipp].other[ido+i]=epsp[i];
00377 }
00378 Mm->ip[ipp].other[ido+n]=gamma;
00379
00380 }
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405 long mohrcoulomb::cutting_plane(long ipp, double &gamma, vector &epsn, vector &epsp, vector &q, long ni, double err)
00406 {
00407 long i,j,ncomp=epsn.n,idem,mu,zps;
00408 double f,denom,dgamma,e,nu;
00409 vector epsa(ncomp), sig(ncomp), epsptr(ncomp);
00410 vector psig(3), peps(3), pdepsp(3), pepsa(3), dfds(3), dgds(3), depsp(ncomp);
00411 matrix d(ncomp,ncomp), sigt(3,3), epsat(3,3), epspt(3,3);
00412 matrix pd(3,3), pvect(3,3);
00413
00414
00415
00416 idem = Mm->ip[ipp].gemid();
00417 if (Mm->ip[ipp].tm[idem] != elisomat)
00418 {
00419 print_err("Mohr-Coulomb material is combined with\n"
00420 " the unsupported elastic material.\n"
00421 " only elisomat type is supported", __FILE__, __LINE__, __func__);
00422 abort ();
00423 }
00424
00425 dgamma=0.0;
00426 e = Mm->give_actual_ym(ipp);
00427 nu = Mm->give_actual_nu(ipp);
00428
00429 Mm->elmatstiff(d, ipp);
00430 pelmatstiff (ipp, pd, e, nu);
00431
00432 if (Mm->ip[ipp].ssst == planestrain)
00433 {
00434 d[0][3] = d[0][1]; d[1][3] = d[1][0];
00435 d[3][0] = d[1][0]; d[3][1] = d[1][0]; d[3][3] = d[1][1];
00436 }
00437 if (Mm->ip[ipp].ssst == planestress)
00438 epsn[3] = -nu / (1.0 - nu) * (epsn[0]+epsn[1]);
00439 copyv(epsp, epsptr);
00440
00441
00442
00443 subv (epsn, epsp, epsa);
00444
00445 mxv (d, epsa, sig);
00446
00447 vector_tensor (sig, sigt, Mm->ip[ipp].ssst, stress);
00448 princ_val (sigt, psig, pvect, nijac, limit, Mp->zero, 3,1);
00449
00450 vector_tensor (epsa, epsat, Mm->ip[ipp].ssst, strain);
00451 glmatrixtransf(epsat, pvect);
00452 peps[0] = epsat[0][0];
00453 peps[1] = epsat[1][1];
00454 peps[2] = epsat[2][2];
00455 if (Mm->ip[ipp].ssst == planestress)
00456 {
00457 zps = checkzeropsig(psig);
00458 fillrow(0.0, zps, pd);
00459 fillcol(0.0, zps, pd);
00460 }
00461
00462 for (i=0;i<ni;i++){
00463
00464 subv (peps, pdepsp, pepsa);
00465
00466 mxv (pd, pepsa, psig);
00467
00468 mu = checkpsig(psig);
00469 if (mu > 0)
00470 return mu;
00471
00472 f = yieldfunction (psig);
00473 if (f<err) break;
00474 if (i==ni-1 && f>err)
00475 {
00476 print_err("yield surface was not reached for ipp = %ld", __FILE__, __LINE__, __func__, ipp);
00477 exit(0);
00478 }
00479
00480 deryieldfsigma (dfds);
00481 derplaspotsigma (dgds);
00482
00483 mxv (pd, dgds, psig);
00484 scprd (dfds, psig, denom);
00485 denom += plasmodscalar(ipp, q);
00486
00487 dgamma = f/denom;
00488
00489 gamma+=dgamma;
00490
00491 for (j=0; j<3; j++)
00492 {
00493 epspt[j][j]+=dgamma*dgds[j];
00494 pdepsp[j]+=dgamma*dgds[j];
00495 }
00496 lgmatrixtransf(epspt, pvect);
00497 tensor_vector(depsp, epspt, Mm->ip[ipp].ssst, strain);
00498 addv(epsp, depsp, epsptr);
00499
00500 updateq(ipp, epsptr, q);
00501 }
00502
00503 copyv(epsptr, epsp);
00504
00505 subv (epsn, epsp, epsa);
00506
00507 mxv (d, epsa, sig);
00508
00509 Mm->storestress (0, ipp, sig);
00510 return 0;
00511 }
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525 double mohrcoulomb::plasmodscalar(long ipp, vector &q)
00526 {
00527 return 0.0;
00528 }
00529
00530
00531
00532
00533
00534
00535
00536
00537
00538
00539
00540
00541
00542
00543 void mohrcoulomb::updateq(long ipp, vector &epsp, vector &q)
00544 {
00545 return;
00546 }
00547
00548
00549
00550
00551
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562 long mohrcoulomb::checkpsig(vector &psig)
00563 {
00564 double s1, s2, sigma, tau, maxsigt;
00565
00566
00567 s1 = psig[0];
00568 s2 = psig[2];
00569 sigma = (s1 + s2) / 2.0;
00570 tau = -(s1 - s2) / (2.0 * cos(phi));
00571 maxsigt = c / tan(phi);
00572 if ((phi != 0.0) && (tau < ((sigma - maxsigt)/tan(phi))))
00573 return 1;
00574
00575
00576 if (psig[0] > psig[1])
00577 return 12;
00578 if (psig[1] > psig[2])
00579 return 23;
00580
00581
00582 return 0;
00583 }
00584
00585
00586
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596 long mohrcoulomb::checkzeropsig(vector &psig)
00597 {
00598 long i;
00599
00600 for (i=0; i<3; i++)
00601 {
00602 if (psig[i] == 0.0)
00603 return i;
00604 }
00605 return -1;
00606 }
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620
00621
00622
00623
00624
00625
00626
00627
00628
00629 void mohrcoulomb::mc_msurf_cp(long ipp, double &gamma, vector &epsn, vector &epsp, vector &q, long mu,long ni,double err)
00630 {
00631 long i,j,ncomp=epsn.n,idem,zps;
00632 double e,nu;
00633 long stat[2], nas;
00634 vector epsa(ncomp), sig(ncomp), epsptr(ncomp), depsp(ncomp);
00635 vector psig(3),peps(3),pdepsp(3),pepsa(3),dpepsp(3), f(2);
00636 matrix d(ncomp,ncomp),sigt(3,3),epsat(3,3),epspt(3,3);
00637 matrix pd(3,3),pvect(3,3);
00638 matrix dfds, dgds, am, cpm, hcpm;
00639 vector dgamma, ff;
00640
00641
00642
00643 idem = Mm->ip[ipp].gemid();
00644 if (Mm->ip[ipp].tm[idem] != elisomat)
00645 {
00646 print_err("Mohr-Coulomb material is combined with\n"
00647 " the nonsupported elastic material.\n"
00648 " only elisomat type is supported", __FILE__, __LINE__, __func__);
00649 exit(0);
00650 }
00651
00652 e = Mm->give_actual_ym(ipp);
00653 nu = Mm->give_actual_nu(ipp);
00654
00655 Mm->elmatstiff(d, ipp);
00656 pelmatstiff (ipp, pd, e, nu);
00657
00658 if (Mm->ip[ipp].ssst == planestrain)
00659 {
00660 d[0][3] = d[0][1]; d[1][3] = d[1][0];
00661 d[3][0] = d[1][0]; d[3][1] = d[1][0]; d[3][3] = d[1][1];
00662 }
00663
00664
00665 copyv(epsp, epsptr);
00666
00667
00668
00669 subv (epsn, epsp, epsa);
00670
00671 mxv (d, epsa, sig);
00672
00673 vector_tensor (sig, sigt, Mm->ip[ipp].ssst, stress);
00674 princ_val (sigt, psig, pvect, nijac, limit, Mp->zero, 3,1);
00675
00676 vector_tensor (epsa, epsat, Mm->ip[ipp].ssst, strain);
00677 glmatrixtransf(epsat, pvect);
00678 peps[0] = epsat[0][0];
00679 peps[1] = epsat[1][1];
00680 peps[2] = epsat[2][2];
00681 if (Mm->ip[ipp].ssst == planestress)
00682 {
00683 zps = checkzeropsig(psig);
00684 fillrow(0.0, zps, pd);
00685 fillcol(0.0, zps, pd);
00686 }
00687
00688
00689
00690 for (i=0; i<ni; i++)
00691 {
00692
00693 subv (peps, pdepsp, pepsa);
00694
00695 mxv (pd, pepsa, psig);
00696
00697 yieldfunction(psig, mu, f);
00698
00699 nas=0;
00700 for (j=0; j<2; j++)
00701 {
00702 if (f[j] >= err)
00703 {
00704 stat[j] = 1;
00705 nas++;
00706 }
00707 else
00708 stat[j] = 0;
00709 }
00710 if (nas==0)
00711 break;
00712
00713
00714 allocm (nas, 3, dfds); allocm (nas, 3, dgds);
00715 allocm (nas, 3, am); allocm (nas, nas, cpm);
00716 allocm (nas, nas, hcpm);
00717 allocv (nas, dgamma);
00718 allocv (nas, ff);
00719 yieldfunction(psig, mu, stat, ff);
00720
00721 dfdsigma(stat, mu, dfds);
00722 dgdsigma(stat, mu, dgds);
00723 plasmodscalar(ipp, q, mu, hcpm);
00724
00725 mxm(dfds, pd, am);
00726 mxmt(am, dgds, cpm);
00727 addm(cpm, hcpm, cpm);
00728
00729 gemp(cpm.a, dgamma.a, ff.a, nas, 1, Mp->zero, 1);
00730
00731 mtxv(dgds, dgamma, dpepsp);
00732
00733 for (j=0; j<3; j++)
00734 {
00735 epspt[j][j]+=dpepsp[j];
00736 pdepsp[j]+=dpepsp[j];
00737 }
00738 for (j=0; j<nas; j++)
00739 gamma += dgamma[j];
00740 lgmatrixtransf(epspt, pvect);
00741 tensor_vector(depsp, epspt, Mm->ip[ipp].ssst, strain);
00742 addv(epsp, depsp, epsptr);
00743
00744 updateq(ipp,epsptr, q);
00745
00746 destrm (dfds); destrm (dgds);
00747 destrm (am); destrm (cpm);
00748 destrm (hcpm);
00749 destrv (dgamma);
00750 destrv(ff);
00751 }
00752
00753 copyv(epsptr, epsp);
00754
00755 subv (epsn, epsp, epsa);
00756
00757 mxv (d, epsa, sig);
00758
00759 Mm->storestress (0, ipp, sig);
00760 return;
00761 }
00762
00763
00764
00765
00766
00767
00768
00769
00770
00771
00772
00773
00774
00775
00776
00777
00778 void mohrcoulomb::plasmodscalar(long ipp, vector &q, long mu, matrix &hcpm)
00779 {
00780 fillm(0.0, hcpm);
00781 }
00782
00783
00784
00785
00786
00787
00788
00789
00790
00791
00792
00793
00794
00795
00796
00797
00798
00799
00800
00801
00802 void mohrcoulomb::yieldfunction(vector &psig, long mu, long *stat, vector &f)
00803 {
00804 long i = 0;
00805 double k1, k2, s1, s2;
00806
00807 s1 = psig[0];
00808 s2 = psig[2];
00809 k1 = -1.0 + sin(phi); k2 = 1.0 + sin(phi);
00810
00811 if (stat[0])
00812 {
00813 s1 = psig[0];
00814 s2 = psig[2];
00815 f[i] = k1*s1 + k2*s2 - 2.0*c*cos(phi);
00816 i++;
00817 }
00818 if (stat[1])
00819 {
00820 switch (mu)
00821 {
00822 case 12:
00823 s1 = psig[1];
00824 s2 = psig[2];
00825 f[i] = k1*s1 + k2*s2 - 2.0*c*cos(phi);
00826 break;
00827 case 23:
00828 s1 = psig[0];
00829 s2 = psig[1];
00830 f[i] = k1*s1 + k2*s2 - 2.0*c*cos(phi);
00831 break;
00832 case 1:
00833 f[i] = s1 + s2 - 2.0*c/tan(phi);
00834 break;
00835 default:
00836 break;
00837 }
00838 }
00839 return;
00840 }
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858
00859
00860
00861 void mohrcoulomb::yieldfunction(vector &psig, long mu, vector &f)
00862 {
00863 double k1, k2, s1, s2;
00864
00865 k1 = -1.0 + sin(phi); k2 = 1.0 + sin(phi);
00866 s1 = psig[0];
00867 s2 = psig[2];
00868
00869 f[0] = k1*s1 + k2*s2 - 2.0*c*cos(phi);
00870 switch (mu)
00871 {
00872 case 12:
00873 s1 = psig[1];
00874 s2 = psig[2];
00875 f[1] = k1*s1 + k2*s2 - 2.0*c*cos(phi);
00876 break;
00877 case 23:
00878 s1 = psig[0];
00879 s2 = psig[1];
00880 f[1] = k1*s1 + k2*s2 - 2.0*c*cos(phi);
00881 break;
00882 case 1:
00883 f[1] = s1 + s2 - 2.0*c/tan(phi);
00884 break;
00885 default:
00886 break;
00887 }
00888 return;
00889 }
00890
00891
00892
00893
00894
00895
00896
00897
00898
00899
00900
00901
00902
00903
00904
00905
00906
00907
00908
00909 void mohrcoulomb::dfdsigma(long *stat, long mu, matrix &dfds)
00910 {
00911 long i = 0;
00912 if (stat[0])
00913 {
00914 dfds[i][0] = -1.0 + sin(phi);
00915 dfds[i][1] = 0.0;
00916 dfds[i][2] = 1.0 + sin(phi);
00917 i++;
00918 }
00919 if (stat[1])
00920 {
00921 switch (mu)
00922 {
00923 case 12:
00924 dfds[i][0] = 0.0;
00925 dfds[i][1] = -1.0 + sin(phi);
00926 dfds[i][2] = 1.0 + sin(phi);
00927 break;
00928 case 23:
00929 dfds[i][0] = -1.0 + sin(phi);
00930 dfds[i][1] = 1.0 + sin(phi);
00931 dfds[i][2] = 0.0;
00932 break;
00933 case 1:
00934 dfds[i][0] = 1.0;
00935 dfds[i][1] = 0.0;
00936 dfds[i][2] = 1.0;
00937 break;
00938 default:
00939 break;
00940 }
00941 }
00942 return;
00943 }
00944
00945
00946
00947
00948
00949
00950
00951
00952
00953
00954
00955
00956
00957
00958
00959
00960
00961
00962
00963 void mohrcoulomb::dgdsigma(long *stat, long mu, matrix &dgds)
00964 {
00965 long i = 0;
00966 if (stat[0])
00967 {
00968 dgds[i][0] = -1.0 + sin(psi);
00969 dgds[i][1] = 0.0;
00970 dgds[i][2] = 1.0 + sin(psi);
00971 i++;
00972 }
00973 if (stat[1])
00974 {
00975 switch (mu)
00976 {
00977 case 12:
00978 dgds[i][0] = 0.0;
00979 dgds[i][1] = -1.0 + sin(psi);
00980 dgds[i][2] = 1.0 + sin(psi);
00981 break;
00982 case 23:
00983 dgds[i][0] = -1.0 + sin(psi);
00984 dgds[i][1] = 1.0 + sin(psi);
00985 dgds[i][2] = 0.0;
00986 break;
00987 case 1:
00988 dgds[i][0] = 1.0;
00989 dgds[i][1] = 0.0;
00990 dgds[i][2] = 1.0;
00991 break;
00992 default:
00993 break;
00994 }
00995 }
00996 return;
00997 }
00998
00999
01000
01001
01002
01003
01004
01005
01006
01007
01008
01009
01010
01011
01012
01013 void mohrcoulomb::updateval (long ipp,long im, long ido)
01014 {
01015 long i,n = Mm->givencompeqother(ipp, im);
01016
01017 for (i=0; i<n; i++)
01018 Mm->ip[ipp].eqother[ido+i]=Mm->ip[ipp].other[ido+i];
01019 }
01020
01021
01022
01023
01024
01025
01026
01027
01028
01029
01030
01031
01032
01033
01034 void mohrcoulomb::giveirrstrains (long ipp, long ido, vector &epsp)
01035 {
01036 long i;
01037 for (i=0;i<epsp.n;i++)
01038 epsp[i] = Mm->ip[ipp].eqother[ido+i];
01039 }
01040
01041
01042
01043
01044
01045
01046
01047
01048
01049
01050
01051
01052
01053 double mohrcoulomb::give_consparam (long ipp,long ido)
01054 {
01055 long ncompstr;
01056 double gamma;
01057
01058 ncompstr=Mm->ip[ipp].ncompstr;
01059 gamma = Mm->ip[ipp].eqother[ido+ncompstr];
01060
01061 return gamma;
01062 }
01063
01064
01065
01066
01067
01068
01069
01070
01071
01072
01073
01074
01075
01076 void mohrcoulomb::changeparam (atsel &atm,vector &val)
01077 {
01078 long i;
01079
01080 for (i=0;i<atm.num;i++){
01081 switch (atm.atrib[i]){
01082 case 0:{
01083 phi=val[i];
01084 break;
01085 }
01086 case 1:{
01087 c=val[i];
01088 break;
01089 }
01090 case 2:{
01091 psi=val[i];
01092 break;
01093 }
01094 default:
01095 print_err("wrong number of atribute", __FILE__, __LINE__, __func__);
01096 }
01097 }
01098 }