Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5313) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5314) @@ -593,10 +593,10 @@ /*Compute Stress*/ - sigma_xx=viscosity*epsilon[0]-pressure*stokesreconditioning; // sigma = nu eps - pressure - sigma_yy=viscosity*epsilon[1]-pressure*stokesreconditioning; - sigma_zz=viscosity*epsilon[2]-pressure*stokesreconditioning; - sigma_xy=viscosity*epsilon[3]; - sigma_xz=viscosity*epsilon[4]; - sigma_yz=viscosity*epsilon[5]; + sigma_xx=2*viscosity*epsilon[0]-pressure*stokesreconditioning; // sigma = nu eps - pressure + sigma_yy=2*viscosity*epsilon[1]-pressure*stokesreconditioning; + sigma_zz=2*viscosity*epsilon[2]-pressure*stokesreconditioning; + sigma_xy=2*viscosity*epsilon[3]; + sigma_xz=2*viscosity*epsilon[4]; + sigma_yz=2*viscosity*epsilon[5]; /*Get normal vector to the bed */ @@ -2491,5 +2491,5 @@ newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); - D_scalar=newviscosity*gauss_weight*Jdet; + D_scalar=2*newviscosity*gauss_weight*Jdet; for (i=0;i<3;i++){ D[i][i]=D_scalar; @@ -2825,5 +2825,5 @@ newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); - D_scalar=newviscosity*gauss_weight*Jdet; + D_scalar=2*newviscosity*gauss_weight*Jdet; for (i=0;i<5;i++){ D[i][i]=D_scalar; @@ -3029,5 +3029,5 @@ D_scalar=gauss_weight*Jdet; for (i=0;i<6;i++){ - D[i][i]=D_scalar*viscosity; + D[i][i]=D_scalar*2*viscosity; } for (i=6;i<8;i++){ @@ -3103,9 +3103,9 @@ DLStokes[4][4]=-alpha2_gauss*gauss_weight*Jdet2d*bed_normal[0]*bed_normal[2]; DLStokes[5][5]=-alpha2_gauss*gauss_weight*Jdet2d*bed_normal[1]*bed_normal[2]; - DLStokes[6][6]=-viscosity*gauss_weight*Jdet2d*bed_normal[0]; - DLStokes[7][7]=-viscosity*gauss_weight*Jdet2d*bed_normal[1]; - DLStokes[8][8]=-viscosity*gauss_weight*Jdet2d*bed_normal[2]; - DLStokes[9][8]=-viscosity*gauss_weight*Jdet2d*bed_normal[0]/2.0; - DLStokes[10][10]=-viscosity*gauss_weight*Jdet2d*bed_normal[1]/2.0; + DLStokes[6][6]=-2*viscosity*gauss_weight*Jdet2d*bed_normal[0]; + DLStokes[7][7]=-2*viscosity*gauss_weight*Jdet2d*bed_normal[1]; + DLStokes[8][8]=-2*viscosity*gauss_weight*Jdet2d*bed_normal[2]; + DLStokes[9][8]=-2*viscosity*gauss_weight*Jdet2d*bed_normal[0]/2.0; + DLStokes[10][10]=-2*viscosity*gauss_weight*Jdet2d*bed_normal[1]/2.0; DLStokes[11][11]=stokesreconditioning*gauss_weight*Jdet2d*bed_normal[0]; DLStokes[12][12]=stokesreconditioning*gauss_weight*Jdet2d*bed_normal[1]; @@ -4230,5 +4230,5 @@ D_scalar=gauss_weight*Jdet; for (i=0;i<6;i++){ - D[i][i]=D_scalar*viscosity; + D[i][i]=D_scalar*2*viscosity; } for (i=6;i<8;i++){ @@ -4972,7 +4972,11 @@ epsilon_sqr[1][2]=pow(epsilon_matrix[1][2],2); epsilon_sqr[2][2]=pow(epsilon_matrix[2][2],2); - 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); - *phi=2*pow(epsilon_eff,2.0)*viscosity; + + /*Phi = Tr(sigma * eps) + * = Tr(sigma'* eps) + * = 2 * eps_eff * sigma'_eff + * = 4 * eps_eff ^2*/ + *phi=4*pow(epsilon_eff,2.0)*viscosity; } /*}}}*/ Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5313) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5314) @@ -3524,5 +3524,5 @@ onto this scalar matrix, so that we win some computational time: */ newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); - D_scalar=newviscosity*thickness*gauss_weight*Jdet; + D_scalar=2*newviscosity*thickness*gauss_weight*Jdet; for (i=0;i<3;i++){ Index: /issm/trunk/src/c/objects/Materials/Matice.cpp =================================================================== --- /issm/trunk/src/c/objects/Materials/Matice.cpp (revision 5313) +++ /issm/trunk/src/c/objects/Materials/Matice.cpp (revision 5314) @@ -280,5 +280,5 @@ void Matice::GetViscosity2d(double* pviscosity, double* epsilon){ /*From a string tensor and a material object, return viscosity, using Glen's flow law. - 2*B + B viscosity= ------------------------------------------------------------------- 2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n] @@ -311,5 +311,5 @@ else{ if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0)){ - viscosity=pow((double)10,(double)14); + viscosity=0.5*pow((double)10,(double)14); } else{ @@ -319,13 +319,13 @@ exy=*(epsilon+2); - /*Build viscosity: viscosity=2*B/(2*A^e) */ + /*Build viscosity: viscosity=B/(2*A^e) */ A=pow(exx,2)+pow(eyy,2)+pow(exy,2)+exx*eyy; if(A==0){ /*Maxiviscositym viscosity for 0 shear areas: */ - viscosity=4.5*pow((double)10,(double)17); + viscosity=2.5*pow((double)10,(double)17); } else{ - e=(n-1)/2/n; - viscosity=2*B/(2*pow(A,e)); + e=(n-1)/(2*n); + viscosity=B/(2*pow(A,e)); } } @@ -346,7 +346,7 @@ /*Return viscosity accounting for steady state power law creep [Thomas and MacAyeal, 1982]: * - * 2*B + * B * viscosity3d= ------------------------------------------------------------------- - * 2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n] + * 2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n] * * where mu is the viscotiy, B the flow law parameter , (u,v) the velocity @@ -393,10 +393,10 @@ if(A==0){ /*Maxiviscosity3dm viscosity for 0 shear areas: */ - viscosity3d=4.5*pow((double)10,(double)17); + viscosity3d=2.25*pow((double)10,(double)17); } else{ e=(n-1)/2/n; - viscosity3d=2*B/(2*pow(A,e)); + viscosity3d=B/(2*pow(A,e)); } } @@ -416,5 +416,5 @@ /*Return viscosity accounting for steady state power law creep [Thomas and MacAyeal, 1982]: * - * 2*B + * B * viscosity3d= ------------------------------------------------------------------- * 2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n] @@ -462,13 +462,13 @@ eyz=*(epsilon+5); - /*Build viscosity: viscosity3d=2*B/(2*A^e) */ + /*Build viscosity: viscosity3d=B/(2*A^e) */ A=pow(exx,2)+pow(eyy,2)+pow(exy,2)+pow(exz,2)+pow(eyz,2)+exx*eyy+pow(eps0,2); if(A==0){ /*Maxiviscosity3dm viscosity for 0 shear areas: */ - viscosity3d=4.5*pow((double)10,(double)17); + viscosity3d=2.25*pow((double)10,(double)17); } else{ e=(n-1)/2/n; - viscosity3d=2*B/(2*pow(A,e)); + viscosity3d=B/(2*pow(A,e)); } } @@ -488,10 +488,7 @@ /*Return viscosity accounting for steady state power law creep [Thomas and MacAyeal, 1982]: * - * 2* (1-n)/2n - * mu2= ------------------------------------------------------------------- - * 2[ (du/dx)^2+(dv/dy)^2+1/4*(du/dy+dv/dx)^2+du/dx*dv/dy ]^[(3n-1)/2n] - * - * where mu2 is the second viscosity, (u,v) the velocity - * vector, and n the flow law exponent. + * 1 + * viscosity= ------------------------------------------------------------------- + * 2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n] * * If epsilon is NULL, it means this is the first time Gradjb is being run, and we @@ -522,8 +519,8 @@ if(A==0){ /*Maximum viscosity_complement for 0 shear areas: */ - viscosity_complement=4.5*pow((double)10,(double)17); + viscosity_complement=2.25*pow((double)10,(double)17); } else{ - e=(n-1)/2/n; + e=(n-1)/(2*n); viscosity_complement=1/(2*pow(A,e));