![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 5314
- Timestamp:
- 08/17/10 12:05:52 (15 years ago)
- Location:
- issm/trunk/src/c/objects
- Files:
-
- 3 edited
-
Elements/Penta.cpp (modified) (7 diffs)
-
Elements/Tria.cpp (modified) (1 diff)
-
Materials/Matice.cpp (modified) (9 diffs)
Legend:
- Unmodified
- Added
- Removed
-
TabularUnified issm/trunk/src/c/objects/Elements/Penta.cpp ¶
r5311 r5314 593 593 594 594 /*Compute Stress*/ 595 sigma_xx= viscosity*epsilon[0]-pressure*stokesreconditioning; // sigma = nu eps - pressure596 sigma_yy= viscosity*epsilon[1]-pressure*stokesreconditioning;597 sigma_zz= viscosity*epsilon[2]-pressure*stokesreconditioning;598 sigma_xy= viscosity*epsilon[3];599 sigma_xz= viscosity*epsilon[4];600 sigma_yz= viscosity*epsilon[5];595 sigma_xx=2*viscosity*epsilon[0]-pressure*stokesreconditioning; // sigma = nu eps - pressure 596 sigma_yy=2*viscosity*epsilon[1]-pressure*stokesreconditioning; 597 sigma_zz=2*viscosity*epsilon[2]-pressure*stokesreconditioning; 598 sigma_xy=2*viscosity*epsilon[3]; 599 sigma_xz=2*viscosity*epsilon[4]; 600 sigma_yz=2*viscosity*epsilon[5]; 601 601 602 602 /*Get normal vector to the bed */ … … 2491 2491 2492 2492 newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); 2493 D_scalar= newviscosity*gauss_weight*Jdet;2493 D_scalar=2*newviscosity*gauss_weight*Jdet; 2494 2494 for (i=0;i<3;i++){ 2495 2495 D[i][i]=D_scalar; … … 2825 2825 2826 2826 newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); 2827 D_scalar= newviscosity*gauss_weight*Jdet;2827 D_scalar=2*newviscosity*gauss_weight*Jdet; 2828 2828 for (i=0;i<5;i++){ 2829 2829 D[i][i]=D_scalar; … … 3029 3029 D_scalar=gauss_weight*Jdet; 3030 3030 for (i=0;i<6;i++){ 3031 D[i][i]=D_scalar* viscosity;3031 D[i][i]=D_scalar*2*viscosity; 3032 3032 } 3033 3033 for (i=6;i<8;i++){ … … 3103 3103 DLStokes[4][4]=-alpha2_gauss*gauss_weight*Jdet2d*bed_normal[0]*bed_normal[2]; 3104 3104 DLStokes[5][5]=-alpha2_gauss*gauss_weight*Jdet2d*bed_normal[1]*bed_normal[2]; 3105 DLStokes[6][6]=- viscosity*gauss_weight*Jdet2d*bed_normal[0];3106 DLStokes[7][7]=- viscosity*gauss_weight*Jdet2d*bed_normal[1];3107 DLStokes[8][8]=- viscosity*gauss_weight*Jdet2d*bed_normal[2];3108 DLStokes[9][8]=- viscosity*gauss_weight*Jdet2d*bed_normal[0]/2.0;3109 DLStokes[10][10]=- viscosity*gauss_weight*Jdet2d*bed_normal[1]/2.0;3105 DLStokes[6][6]=-2*viscosity*gauss_weight*Jdet2d*bed_normal[0]; 3106 DLStokes[7][7]=-2*viscosity*gauss_weight*Jdet2d*bed_normal[1]; 3107 DLStokes[8][8]=-2*viscosity*gauss_weight*Jdet2d*bed_normal[2]; 3108 DLStokes[9][8]=-2*viscosity*gauss_weight*Jdet2d*bed_normal[0]/2.0; 3109 DLStokes[10][10]=-2*viscosity*gauss_weight*Jdet2d*bed_normal[1]/2.0; 3110 3110 DLStokes[11][11]=stokesreconditioning*gauss_weight*Jdet2d*bed_normal[0]; 3111 3111 DLStokes[12][12]=stokesreconditioning*gauss_weight*Jdet2d*bed_normal[1]; … … 4230 4230 D_scalar=gauss_weight*Jdet; 4231 4231 for (i=0;i<6;i++){ 4232 D[i][i]=D_scalar* viscosity;4232 D[i][i]=D_scalar*2*viscosity; 4233 4233 } 4234 4234 for (i=6;i<8;i++){ … … 4972 4972 epsilon_sqr[1][2]=pow(epsilon_matrix[1][2],2); 4973 4973 epsilon_sqr[2][2]=pow(epsilon_matrix[2][2],2); 4974 4975 4974 epsilon_eff=1/pow(2,0.5)*pow((epsilon_sqr[0][0]+epsilon_sqr[0][1]+ epsilon_sqr[0][2]+ epsilon_sqr[1][0]+ epsilon_sqr[1][1]+ epsilon_sqr[1][2]+ epsilon_sqr[2][0]+ epsilon_sqr[2][1]+ epsilon_sqr[2][2]),0.5); 4976 *phi=2*pow(epsilon_eff,2.0)*viscosity; 4975 4976 /*Phi = Tr(sigma * eps) 4977 * = Tr(sigma'* eps) 4978 * = 2 * eps_eff * sigma'_eff 4979 * = 4 * eps_eff ^2*/ 4980 *phi=4*pow(epsilon_eff,2.0)*viscosity; 4977 4981 } 4978 4982 /*}}}*/ -
TabularUnified issm/trunk/src/c/objects/Elements/Tria.cpp ¶
r5311 r5314 3524 3524 onto this scalar matrix, so that we win some computational time: */ 3525 3525 newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); 3526 D_scalar= newviscosity*thickness*gauss_weight*Jdet;3526 D_scalar=2*newviscosity*thickness*gauss_weight*Jdet; 3527 3527 3528 3528 for (i=0;i<3;i++){ -
TabularUnified issm/trunk/src/c/objects/Materials/Matice.cpp ¶
r5311 r5314 280 280 void Matice::GetViscosity2d(double* pviscosity, double* epsilon){ 281 281 /*From a string tensor and a material object, return viscosity, using Glen's flow law. 282 2*B282 B 283 283 viscosity= ------------------------------------------------------------------- 284 284 2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n] … … 311 311 else{ 312 312 if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0)){ 313 viscosity= pow((double)10,(double)14);313 viscosity=0.5*pow((double)10,(double)14); 314 314 } 315 315 else{ … … 319 319 exy=*(epsilon+2); 320 320 321 /*Build viscosity: viscosity= 2*B/(2*A^e) */321 /*Build viscosity: viscosity=B/(2*A^e) */ 322 322 A=pow(exx,2)+pow(eyy,2)+pow(exy,2)+exx*eyy; 323 323 if(A==0){ 324 324 /*Maxiviscositym viscosity for 0 shear areas: */ 325 viscosity= 4.5*pow((double)10,(double)17);325 viscosity=2.5*pow((double)10,(double)17); 326 326 } 327 327 else{ 328 e=(n-1)/ 2/n;329 viscosity= 2*B/(2*pow(A,e));328 e=(n-1)/(2*n); 329 viscosity=B/(2*pow(A,e)); 330 330 } 331 331 } … … 346 346 /*Return viscosity accounting for steady state power law creep [Thomas and MacAyeal, 1982]: 347 347 * 348 * 2*B348 * B 349 349 * viscosity3d= ------------------------------------------------------------------- 350 * 2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n]350 * 2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n] 351 351 * 352 352 * where mu is the viscotiy, B the flow law parameter , (u,v) the velocity … … 393 393 if(A==0){ 394 394 /*Maxiviscosity3dm viscosity for 0 shear areas: */ 395 viscosity3d= 4.5*pow((double)10,(double)17);395 viscosity3d=2.25*pow((double)10,(double)17); 396 396 } 397 397 else{ 398 398 e=(n-1)/2/n; 399 399 400 viscosity3d= 2*B/(2*pow(A,e));400 viscosity3d=B/(2*pow(A,e)); 401 401 } 402 402 } … … 416 416 /*Return viscosity accounting for steady state power law creep [Thomas and MacAyeal, 1982]: 417 417 * 418 * 2*B418 * B 419 419 * viscosity3d= ------------------------------------------------------------------- 420 420 * 2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n] … … 462 462 eyz=*(epsilon+5); 463 463 464 /*Build viscosity: viscosity3d= 2*B/(2*A^e) */464 /*Build viscosity: viscosity3d=B/(2*A^e) */ 465 465 A=pow(exx,2)+pow(eyy,2)+pow(exy,2)+pow(exz,2)+pow(eyz,2)+exx*eyy+pow(eps0,2); 466 466 if(A==0){ 467 467 /*Maxiviscosity3dm viscosity for 0 shear areas: */ 468 viscosity3d= 4.5*pow((double)10,(double)17);468 viscosity3d=2.25*pow((double)10,(double)17); 469 469 } 470 470 else{ 471 471 e=(n-1)/2/n; 472 viscosity3d= 2*B/(2*pow(A,e));472 viscosity3d=B/(2*pow(A,e)); 473 473 } 474 474 } … … 488 488 /*Return viscosity accounting for steady state power law creep [Thomas and MacAyeal, 1982]: 489 489 * 490 * 2* (1-n)/2n 491 * mu2= ------------------------------------------------------------------- 492 * 2[ (du/dx)^2+(dv/dy)^2+1/4*(du/dy+dv/dx)^2+du/dx*dv/dy ]^[(3n-1)/2n] 493 * 494 * where mu2 is the second viscosity, (u,v) the velocity 495 * vector, and n the flow law exponent. 490 * 1 491 * viscosity= ------------------------------------------------------------------- 492 * 2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n] 496 493 * 497 494 * If epsilon is NULL, it means this is the first time Gradjb is being run, and we … … 522 519 if(A==0){ 523 520 /*Maximum viscosity_complement for 0 shear areas: */ 524 viscosity_complement= 4.5*pow((double)10,(double)17);521 viscosity_complement=2.25*pow((double)10,(double)17); 525 522 } 526 523 else{ 527 e=(n-1)/ 2/n;524 e=(n-1)/(2*n); 528 525 529 526 viscosity_complement=1/(2*pow(A,e));
Note:
See TracChangeset
for help on using the changeset viewer.
