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