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