00001 #include <stdio.h>
00002 #include <stdlib.h>
00003 #include <math.h>
00004 #include "splas1d.h"
00005 #include "global.h"
00006 #include "vecttens.h"
00007 #include "intpoints.h"
00008 #include "matrix.h"
00009 #include "vector.h"
00010 #include "stochdriver.h"
00011
00012
00013 splas1d::splas1d (void)
00014 {
00015
00016 fs=0.0;
00017 }
00018
00019 splas1d::~splas1d (void)
00020 {
00021
00022 }
00023
00024
00025
00026
00027
00028
00029
00030
00031 void splas1d::read (XFILE *in)
00032 {
00033
00034 xfscanf (in,"%lf",&fs);
00035
00036 sra.read (in);
00037
00038 hs.read (in);
00039 }
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053 double splas1d::yieldfunction (matrix &sig,vector &q)
00054 {
00055 double f;
00056 f = fabs(sig[0][0])-(fs+q[0]);
00057 return f;
00058 }
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073 void splas1d::dfdsigma (matrix &sig,matrix &dfds)
00074 {
00075 if (sig[0][0]<0.0) dfds[0][0] = -1.0;
00076 else dfds[0][0] = 1.0;
00077 }
00078
00079
00080
00081
00082
00083
00084
00085 void splas1d::dfdsigmadsigma(matrix &dfdsdst)
00086 {
00087 fillm (0.0,dfdsdst);
00088 }
00089
00090
00091
00092
00093
00094
00095
00096
00097 void splas1d::dfdqpar (vector &dq)
00098 {
00099 dq[0]=-1.0;
00100 }
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111 void splas1d::dfdsigmadq (matrix &dfdsdq)
00112 {
00113 fillm (0.0,dfdsdq);
00114 }
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127 void splas1d::dhdgamma (long ipp,vector &epsp,matrix &sig,vector &dhdc)
00128 {
00129 matrix dgds(3,3);
00130
00131
00132 dfdsigma (sig,dgds);
00133
00134
00135
00136
00137 dhdc[0] *= -1.0;
00138 }
00139
00140
00141
00142
00143
00144
00145
00146 void splas1d::plasmod (matrix &h)
00147 {
00148 h[0][0]=0.0;
00149 }
00150
00151 double splas1d::plasmodscalar(matrix &sig,vector &epsp,vector &qtr,double gamma)
00152 {
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189 double ret,a,b,c,e;
00190 vector dfq(qtr.n),res(1);
00191 matrix dfds(1,1),h(qtr.n, qtr.n),hp(1,1);
00192
00193 dfdsigma (sig,dfds);
00194 dfdqpar(dfq);
00195 plasmod (h);
00196
00197 a=dfds[0][0];
00198 b=dfq[0];
00199 c=h[0][0];
00200 e=epsp[0];
00201
00202 if (e<1.0e-6)
00203 e=1.0e-6;
00204
00205 ret = b*0.5*c*a/sqrt(e);
00206
00207 ret=-1.0*ret;
00208 return ret;
00209 }
00210
00211 void splas1d::updateq(long ipp,double dgamma, vector &epsp,vector &q)
00212 {
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229 }
00230
00231 void splas1d::nlstresses (long ipp, long im, long ido)
00232 {
00233 long ni;
00234 double err;
00235 vector epsn(1),epsp(1),q(1);
00236 double gamma;
00237
00238
00239 epsn[0] = Mm->ip[ipp].strain[0];
00240 epsp[0] = Mm->ip[ipp].eqother[ido+0];
00241 gamma = Mm->ip[ipp].eqother[ido+1];
00242 q[0] = Mm->ip[ipp].eqother[ido+2];
00243
00244
00245 if (sra.give_tsra () == cp){
00246 ni=sra.give_ni ();
00247 err=sra.give_err ();
00248
00249
00250
00251 Mm->newton_stress_return (ipp,im,ido,gamma,epsn,epsp,q,ni,err);
00252
00253 }
00254 else{
00255 print_err("wrong type of stress return algorithm is required", __FILE__, __LINE__, __func__);
00256 abort ();
00257 }
00258
00259
00260 Mm->ip[ipp].other[ido+0]=epsp[0];
00261 Mm->ip[ipp].other[ido+1]=gamma;
00262 Mm->ip[ipp].other[ido+2]=q[0];
00263 }
00264
00265 void splas1d::nonloc_nlstresses (long ipp, long im, long ido)
00266 {
00267 long ni;
00268 double err;
00269 vector epsn(1),epsp(1),q(1);
00270 double gamma;
00271
00272
00273 epsn[0] = Mm->ip[ipp].strain[0];
00274 gamma = Mm->ip[ipp].eqother[ido+1];
00275 epsp[0] = Mm->ip[ipp].nonloc[ido+0];
00276 q[0] = Mm->ip[ipp].eqother[ido+2];
00277
00278
00279 if (sra.give_tsra () == cp){
00280 ni=sra.give_ni ();
00281 err=sra.give_err ();
00282
00283 Mm->cutting_plane (ipp,im,ido,gamma,epsn,epsp,q,ni,err);
00284 }
00285 else{
00286 fprintf (stderr,"\n\n wrong type of stress return algorithm is required in nlstresses (file %s, line %d).\n",__FILE__,__LINE__);
00287 abort ();
00288 }
00289
00290
00291 Mm->ip[ipp].other[ido+0]=epsp[0];
00292 Mm->ip[ipp].other[ido+1]=gamma;
00293 Mm->ip[ipp].other[ido+2]=q[0];
00294 }
00295
00296 void splas1d::matstiff (matrix &d,long ipp,long ido)
00297 {
00298 double denom1,denom2;
00299 vector sig(1),dq(1),u(1);
00300 matrix de(1,1),h(1,1),g(1,1),sigt(3,3);
00301
00302 if (Mp->nlman->stmat==0){
00303 Mm->elmatstiff (d,ipp);
00304 }
00305 else{
00306
00307
00308 Mm->elmatstiff (de,ipp);
00309
00310
00311 sig[0]=Mm->ip[ipp].stress[0];
00312 if (Mm->ip[ipp].eqother[ido+1]<Mp->zero){
00313 Mm->elmatstiff (d,ipp);
00314 }
00315 else{
00316
00317
00318 dfdsigma (sigt,sigt);
00319 tensor_vector (sig,sigt,Mm->ip[ipp].ssst,stress);
00320
00321 vxv (sig,sig,h);
00322 mxm (de,h,g);
00323 mxm (g,de,h);
00324
00325 mxv (de,sig,u);
00326 scprd (u,sig,denom1);
00327
00328 dfdqpar (dq);
00329 scprd (dq,dq,denom2);
00330
00331 cmulm (1.0/(denom1+denom2),h);
00332
00333 subm (de,h,d);
00334 }
00335 }
00336 }
00337
00338 void splas1d::updateval (long ipp, long ido)
00339 {
00340 Mm->ip[ipp].eqother[ido+0]=Mm->ip[ipp].other[ido+0];
00341 Mm->ip[ipp].eqother[ido+1]=Mm->ip[ipp].other[ido+1];
00342 Mm->ip[ipp].eqother[ido+2]=Mm->ip[ipp].other[ido+2];
00343 }
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356 void splas1d::giveirrstrains (long ipp, long ido, vector &epsp)
00357 {
00358 epsp[0] = Mm->ip[ipp].eqother[ido+0];
00359 }
00360
00361
00362
00363 void splas1d::changeparam (atsel &atm,vector &val)
00364 {
00365 long i;
00366
00367 for (i=0;i<atm.num;i++){
00368 switch (atm.atrib[i]){
00369 case 0:{
00370 fs=val[i];
00371 break;
00372 }
00373 case 1:{
00374
00375 break;
00376 }
00377 default:{
00378 print_err("wrong attribute number",__FILE__,__LINE__,__func__);
00379 }
00380 }
00381 }
00382 }
00383
00384 double splas1d::give_consparam (long ipp, long ido)
00385 {
00386 double gamma;
00387 gamma = Mm->ip[ipp].other[ido+1];
00388 return gamma;
00389 }
00390
00391 long splas1d::give_num_interparam ()
00392 {
00393 return 1;
00394 }
00395
00396 void splas1d::give_interparam (long ipp,long ido,vector &q)
00397 {
00398 long ncompstr=Mm->ip[ipp].ncompstr;
00399
00400 q[0]=Mm->ip[ipp].eqother[ido+ncompstr+1];
00401 }