00001 #include "tensor.h"
00002 #include <math.h>
00003 #include <stdlib.h>
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019 void deviator (matrix &tens,matrix &dev)
00020 {
00021 double isopart;
00022
00023
00024 isopart = first_invar (tens)/3.0;
00025
00026
00027 if (tens.a != dev.a)
00028 copym (tens,dev);
00029 dev[0][0]=tens[0][0]-isopart;
00030 dev[1][1]=tens[1][1]-isopart;
00031 dev[2][2]=tens[2][2]-isopart;
00032 }
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049 void normedtensor (matrix &tens,matrix &ntens)
00050 {
00051 double norm,zero;
00052 zero=1.0e-10;
00053 norm = tensornorm (tens);
00054 copym (tens,ntens);
00055 if (norm<zero){}
00056 else{
00057 cmulm (1.0/norm,ntens);
00058 }
00059 }
00060
00061
00062
00063
00064
00065
00066 double tensornorm (matrix &tens)
00067 {
00068 long i,j;
00069 double nor;
00070
00071 nor=0.0;
00072 for (i=0;i<3;i++){
00073 for (j=0;j<3;j++){
00074 nor+=tens[i][j]*tens[i][j];
00075 }
00076 }
00077 nor=sqrt(nor);
00078 return nor;
00079 }
00080
00081
00082
00083
00084
00085
00086
00087
00088 void cumulstrain (matrix &eps,double &cs)
00089 {
00090 cs+=sqrt(2.0/3.0)*tensornorm (eps);
00091 }
00092
00093
00094
00095
00096
00097 double first_invar (matrix &tens)
00098 {
00099 double inv;
00100 inv=tens[0][0]+tens[1][1]+tens[2][2];
00101 return inv;
00102 }
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113 double second_invar (matrix &tens)
00114 {
00115 double inv;
00116 inv=tens[0][0]*tens[1][1]-tens[0][1]*tens[1][0];
00117 inv+=tens[0][0]*tens[2][2]-tens[0][2]*tens[2][0];
00118 inv+=tens[1][1]*tens[2][2]-tens[1][2]*tens[2][1];
00119 return (inv*(-1.0));
00120 }
00121
00122
00123
00124
00125
00126 double third_invar (matrix &tens)
00127 {
00128 double inv;
00129 inv=tens[0][0]*tens[1][1]*tens[2][2]+tens[0][1]*tens[1][2]*tens[2][0]+tens[0][2]*tens[1][0]*tens[2][1];
00130 inv-=tens[0][2]*tens[1][1]*tens[2][0]+tens[1][2]*tens[2][1]*tens[0][0]+tens[2][2]*tens[0][1]*tens[1][0];
00131 return inv;
00132 }
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160 void glob_loc_tens_trans (vector &tens, matrix &tmat)
00161 {
00162 int n;
00163
00164 switch ( tens.n )
00165 {
00166 case 3:
00167 if ( tmat.n != 2 ) fprintf(stderr,"\nError multiplying matrix by vector - incompatible size (tens_trans in file tensor.*)\n");
00168 if ( tmat.m != 2 ) fprintf(stderr,"\nError multiplying matrix by vector - incompatible size (tens_trans in file tensor.*)\n");
00169 n=2;
00170 break;
00171
00172 case 4:
00173 if ( tmat.n != 3 ) fprintf(stderr,"\nError multiplying matrix by vector - incompatible size (tens_trans in file tensor.*)\n");
00174 if ( tmat.m != 3 ) fprintf(stderr,"\nError multiplying matrix by vector - incompatible size (tens_trans in file tensor.*)\n");
00175
00176 n=3;
00177 break;
00178
00179 case 6:
00180 if ( tmat.n != 3 ) fprintf(stderr,"\nError multiplying matrix by vector - incompatible size (tens_trans in file tensor.*)\n");
00181 if ( tmat.m != 3 ) fprintf(stderr,"\nError multiplying matrix by vector - incompatible size (tens_trans in file tensor.*)\n");
00182
00183 n=3;
00184 break;
00185
00186 default :
00187 fprintf(stderr,"\nError multiplying matrix by vector - incompatible size (tens_trans in file tensor.*)\n");
00188 }
00189
00190 matrix a(n,n);
00191 matrix b(n,n);
00192
00193
00194 switch ( tens.n )
00195 {
00196 case 3 :
00197
00198 a[0][0]=tens[0];
00199 a[1][1]=tens[1];
00200 a[0][1]=tens[2];
00201 a[1][0]=tens[2];
00202 mtxm(tmat,a,b);
00203 mxm(b,tmat,a);
00204 tens[0]=a[0][0];
00205 tens[1]=a[1][1];
00206 tens[2]=a[0][1];
00207 break;
00208
00209 case 4 :
00210 a[0][0]=tens[0];
00211 a[1][1]=tens[1];
00212 a[2][2]=tens[2];
00213 a[1][2]=tens[3];
00214 a[2][1]=tens[3];
00215 mtxm(tmat,a,b);
00216 mxm(b,tmat,a);
00217 tens[0]=a[0][0];
00218 tens[1]=a[1][1];
00219 tens[2]=a[2][2];
00220 tens[3]=a[1][2];
00221 break;
00222
00223 case 6 :
00224 a[0][0]=tens[0];
00225 a[1][1]=tens[1];
00226 a[2][2]=tens[2];
00227 a[1][2]=tens[3];
00228 a[2][1]=tens[3];
00229 a[0][2]=tens[4];
00230 a[2][0]=tens[4];
00231 a[0][1]=tens[5];
00232 a[1][0]=tens[5];
00233 mtxm(tmat,a,b);
00234 mxm(b,tmat,a);
00235 tens[0]=a[0][0];
00236 tens[1]=a[1][1];
00237 tens[2]=a[2][2];
00238 tens[3]=a[1][2];
00239 tens[4]=a[0][2];
00240 tens[5]=a[0][1];
00241 break;
00242 }
00243 }
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271 void loc_glob_tens_trans (vector &tens, matrix &tmat)
00272 {
00273 int n;
00274
00275 switch ( tens.n )
00276 {
00277 case 3:
00278 if ( tmat.n != 2 ) fprintf(stderr,"\nError multiplying matrix by vector - incompatible size (tens_trans in file tensor.*)\n");
00279 if ( tmat.m != 2 ) fprintf(stderr,"\nError multiplying matrix by vector - incompatible size (tens_trans in file tensor.*)\n");
00280 n=2;
00281 break;
00282
00283 case 4:
00284 if ( tmat.n != 3 ) fprintf(stderr,"\nError multiplying matrix by vector - incompatible size (tens_trans in file tensor.*)\n");
00285 if ( tmat.m != 3 ) fprintf(stderr,"\nError multiplying matrix by vector - incompatible size (tens_trans in file tensor.*)\n");
00286
00287 n=3;
00288 break;
00289
00290 case 6:
00291 if ( tmat.n != 3 ) fprintf(stderr,"\nError multiplying matrix by vector - incompatible size (tens_trans in file tensor.*)\n");
00292 if ( tmat.m != 3 ) fprintf(stderr,"\nError multiplying matrix by vector - incompatible size (tens_trans in file tensor.*)\n");
00293
00294 n=3;
00295 break;
00296
00297 default :
00298 fprintf(stderr,"\nError multiplying matrix by vector - incompatible size (tens_trans in file tensor.*)\n");
00299 }
00300
00301 matrix a(n,n);
00302 matrix b(n,n);
00303
00304
00305 switch ( tens.n )
00306 {
00307 case 3 :
00308
00309 a[0][0]=tens[0];
00310 a[1][1]=tens[1];
00311 a[0][1]=tens[2];
00312 a[1][0]=tens[2];
00313 mxm(tmat,a,b);
00314 mxmt(b,tmat,a);
00315 tens[0]=a[0][0];
00316 tens[1]=a[1][1];
00317 tens[2]=a[0][1];
00318 break;
00319
00320 case 4 :
00321 a[0][0]=tens[0];
00322 a[1][1]=tens[1];
00323 a[2][2]=tens[2];
00324 a[1][2]=tens[3];
00325 a[2][1]=tens[3];
00326 mxm(tmat,a,b);
00327 mxmt(b,tmat,a);
00328 tens[0]=a[0][0];
00329 tens[1]=a[1][1];
00330 tens[2]=a[2][2];
00331 tens[3]=a[1][2];
00332 break;
00333
00334 case 6 :
00335 a[0][0]=tens[0];
00336 a[1][1]=tens[1];
00337 a[2][2]=tens[2];
00338 a[1][2]=tens[3];
00339 a[2][1]=tens[3];
00340 a[0][2]=tens[4];
00341 a[2][0]=tens[4];
00342 a[0][1]=tens[5];
00343 a[1][0]=tens[5];
00344 mxm(tmat,a,b);
00345 mxmt(b,tmat,a);
00346 tens[0]=a[0][0];
00347 tens[1]=a[1][1];
00348 tens[2]=a[2][2];
00349 tens[3]=a[1][2];
00350 tens[4]=a[0][2];
00351 tens[5]=a[0][1];
00352 break;
00353 }
00354 }
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364 double tensor2_inner_prod (matrix &a,matrix &b)
00365 {
00366 long i,j,n;
00367 double ip;
00368
00369 n=a.n;
00370 if (a.m!=3 || a.n!=3 || b.m!=3 || b.n!=3){
00371 fprintf (stderr,"\n\n incompatible tensor(s) in function tensor2_inner_prod (file %s, line %d),\n",__FILE__,__LINE__);
00372 }
00373
00374 ip=0.0;
00375 for (i=0;i<n;i++){
00376 for (j=0;j<n;j++){
00377 ip+=a[i][j]*b[i][j];
00378 }
00379 }
00380
00381 return ip;
00382 }
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397 void f_tensor(matrix &a, double (*f)(double), long nijac, double error, double zero, matrix &af)
00398 {
00399 vector pa(3);
00400 matrix t(3,3);
00401 long i;
00402 princ_val (a, pa, t, nijac, error, zero, 3, 1);
00403 fillm(0.0, af);
00404 for(i=0; i<3; i++)
00405 af[i][i] = (*f)(pa[i]);
00406 lgmatrixtransf(af, t);
00407 }
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420 void f_tensor(matrix &a, double (*f)(double), matrix &t, matrix &af)
00421 {
00422 long i;
00423 copym(a, af);
00424 glmatrixtransf(af, t);
00425 fillm(0.0, af);
00426 for(i=0; i<3; i++)
00427 af[i][i] = (*f)(a[i][i]);
00428 lgmatrixtransf(af, t);
00429 }
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444 long dpdir_da(matrix &t, vector &pa, long alpha, long i, matrix &dpd_da)
00445 {
00446 long k,l,beta,ret;
00447
00448 ret = 0;
00449 fillm(0.0, dpd_da);
00450 for(k=0; k<3; k++)
00451 {
00452 for(l=0; l<3; l++)
00453 {
00454 for(beta=0; beta<3; beta++)
00455 {
00456 if (alpha == beta)
00457 continue;
00458 if (pa[alpha] == pa[beta])
00459 {
00460 ret++;
00461 continue;
00462 }
00463
00464
00465
00466
00467
00468
00469 dpd_da[k][l] += (t[k][alpha]*t[l][beta]+t[k][beta]*t[l][alpha])*t[i][beta]/(pa[alpha]-pa[beta]);
00470 }
00471 dpd_da[k][l] /= 2.0;
00472 }
00473 }
00474 return ret;
00475 }