Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 21825) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 21826) @@ -48,5 +48,5 @@ IssmDouble vx,vy,vz,vmag; IssmDouble dvx[3],dvy[3],dvz[3],dvmag[3]; - IssmDouble epseff,epsprime; + IssmDouble eps[3][3],epseff,epsprime; int dim; IssmDouble *xyz_list = NULL; @@ -101,5 +101,16 @@ dvmag[2]=1./(2*sqrt(vmag))*(2*vx*dvx[2]+2*vy*dvy[2]+2*vz*dvz[2]); } - EstarStrainrateQuantities(&epseff,&epsprime,vx,vy,vz,vmag,&dvx[0],&dvy[0],&dvz[0],&dvmag[0]); + + /*Build strain rate tensor*/ + eps[0][0] = dvx[0]; eps[0][1] = .5*(dvx[1]+dvy[0]); eps[0][2] = .5*(dvx[2]+dvz[0]); + eps[1][0] = .5*(dvx[1]+dvy[0]); eps[1][1] = dvy[1]; eps[1][2] = .5*(dvy[2]+dvz[1]); + eps[2][0] = .5*(dvx[2]+dvz[0]); eps[2][1] = .5*(dvy[2]+dvz[1]); eps[2][2] = dvz[2]; + + /*effective strain rate*/ + epseff = 0.; + /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */ + epseff = sqrt(eps[0][0]*eps[0][0] + eps[1][1]*eps[1][1] + eps[0][1]*eps[0][1] + eps[0][2]*eps[0][2] + eps[1][2]*eps[1][2] + eps[0][0]*eps[1][1]); + + EstarStrainrateQuantities(&epsprime,vx,vy,vz,vmag,&dvx[0],&dvy[0],&dvz[0],&dvmag[0]); lambdas[iv]=EstarLambdaS(epseff,epsprime); } @@ -677,5 +688,5 @@ break; case MatestarEnum: - material->ViscosityBFS(&dmudB,dim,xyz_list,gauss,vx_input,vy_input,vz_input); + material->ViscosityBFS(&dmudB,dim,xyz_list,gauss,vx_input,vy_input,vz_input,eps_eff); break; default: _error_("not supported"); @@ -714,5 +725,5 @@ break; case MatestarEnum: - material->ViscosityBHO(&dmudB,dim,xyz_list,gauss,vx_input,vy_input); + material->ViscosityBHO(&dmudB,dim,xyz_list,gauss,vx_input,vy_input,eps_eff); break; default: _error_("not supported"); @@ -751,5 +762,5 @@ break; case MatestarEnum: - material->ViscosityBSSA(&dmudB,dim,xyz_list,gauss,vx_input,vy_input); + material->ViscosityBSSA(&dmudB,dim,xyz_list,gauss,vx_input,vy_input,eps_eff); break; default: _error_("not supported"); Index: /issm/trunk-jpl/src/c/classes/Materials/Material.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Material.h (revision 21825) +++ /issm/trunk-jpl/src/c/classes/Materials/Material.h (revision 21826) @@ -53,7 +53,7 @@ virtual void ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0; virtual void ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon)=0; - virtual void ViscosityBFS(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input)=0; - virtual void ViscosityBHO(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0; - virtual void ViscosityBSSA(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0; + virtual void ViscosityBFS(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input,IssmDouble epseff)=0; + virtual void ViscosityBHO(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,IssmDouble epseff)=0; + virtual void ViscosityBSSA(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,IssmDouble epseff)=0; }; Index: /issm/trunk-jpl/src/c/classes/Materials/Matestar.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matestar.cpp (revision 21825) +++ /issm/trunk-jpl/src/c/classes/Materials/Matestar.cpp (revision 21826) @@ -249,5 +249,5 @@ } /*}}}*/ -IssmDouble Matestar::GetViscosityGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,bool isdepthaveraged){/*{{{*/ +IssmDouble Matestar::GetViscosityGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,IssmDouble eps_eff,bool isdepthaveraged){/*{{{*/ /*output: */ @@ -255,5 +255,5 @@ /*Intermediaries*/ - IssmDouble epseff,epsprime_norm; + IssmDouble epsprime_norm; IssmDouble lambdas; IssmDouble vmag,dvmag[3]; @@ -273,6 +273,6 @@ } - EstarStrainrateQuantities(&epseff,&epsprime_norm,vx,vy,vz,vmag,dvx,dvy,dvz,&dvmag[0]); - lambdas=EstarLambdaS(epseff,epsprime_norm); + EstarStrainrateQuantities(&epsprime_norm,vx,vy,vz,vmag,dvx,dvy,dvz,&dvmag[0]); + lambdas=EstarLambdaS(eps_eff,epsprime_norm); /*Get B and enhancement*/ @@ -294,5 +294,5 @@ /*Compute viscosity*/ /*if no strain rate, return maximum viscosity*/ - if(epseff==0.){ + if(eps_eff==0.){ viscosity = 1.e+14/2.; //viscosity = B; @@ -300,6 +300,9 @@ } else{ - viscosity = B/(2.*pow(E,1./n)*pow(epseff,2./n)); - } + viscosity = B/(2.*pow(E,1./n)*pow(eps_eff,2./n)); + } + + /*Checks in debugging mode*/ + if(viscosity<=0) _error_("Negative viscosity"); /*Assign output pointer*/ @@ -307,9 +310,9 @@ } /*}}}*/ -IssmDouble Matestar::GetViscosity_BGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,bool isdepthaveraged){/*{{{*/ +IssmDouble Matestar::GetViscosity_BGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,IssmDouble eps_eff,bool isdepthaveraged){/*{{{*/ /*Intermediaries*/ IssmDouble dmudB; - IssmDouble epseff,epsprime_norm; + IssmDouble epsprime_norm; IssmDouble lambdas; IssmDouble vmag,dvmag[3]; @@ -329,6 +332,6 @@ } - EstarStrainrateQuantities(&epseff,&epsprime_norm,vx,vy,vz,vmag,dvx,dvy,dvz,&dvmag[0]); - lambdas=EstarLambdaS(epseff,epsprime_norm); + EstarStrainrateQuantities(&epsprime_norm,vx,vy,vz,vmag,dvx,dvy,dvz,&dvmag[0]); + lambdas=EstarLambdaS(eps_eff,epsprime_norm); /*Get enhancement*/ @@ -346,6 +349,6 @@ /*Compute dmudB*/ - if(epseff==0.) dmudB = 0.; - else dmudB = 1./(2.*pow(E,1./3.)*pow(epseff,2./3.)); + if(eps_eff==0.) dmudB = 0.; + else dmudB = 1./(2.*pow(E,1./3.)*pow(eps_eff,2./3.)); /*Assign output*/ @@ -408,5 +411,5 @@ } /*}}}*/ -void Matestar::ViscosityBFS(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/ +void Matestar::ViscosityBFS(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input,IssmDouble eps_eff){/*{{{*/ /*Intermediaries*/ @@ -434,8 +437,8 @@ /*Compute dmudB*/ - *pdmudB=GetViscosity_BGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],isdepthaveraged); -} -/*}}}*/ -void Matestar::ViscosityBHO(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ + *pdmudB=GetViscosity_BGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],eps_eff,isdepthaveraged); +} +/*}}}*/ +void Matestar::ViscosityBHO(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,IssmDouble eps_eff){/*{{{*/ /*Intermediaries*/ @@ -463,7 +466,7 @@ /*Compute viscosity*/ - *pdmudB=GetViscosity_BGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],isdepthaveraged); -}/*}}}*/ -void Matestar::ViscosityBSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ + *pdmudB=GetViscosity_BGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],eps_eff,isdepthaveraged); +}/*}}}*/ +void Matestar::ViscosityBSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,IssmDouble eps_eff){/*{{{*/ /*Intermediaries*/ IssmDouble vx,vy,vz; @@ -493,5 +496,5 @@ /*Compute viscosity*/ - *pdmudB=GetViscosity_BGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],isdepthaveraged); + *pdmudB=GetViscosity_BGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],eps_eff,isdepthaveraged); }/*}}}*/ void Matestar::ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/ @@ -500,4 +503,7 @@ IssmDouble vx,vy,vz; IssmDouble dvx[3],dvy[3],dvz[3]; + IssmDouble epsilon3d[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ + IssmDouble epsilon2d[3]; /* epsilon=[exx,eyy,exy];*/ + IssmDouble eps_eff,eps0=1.e-27; bool isdepthaveraged=0.; @@ -520,6 +526,17 @@ } + if(dim==3){ + /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */ + element->StrainRateFS(&epsilon3d[0],xyz_list,gauss,vx_input,vy_input,vz_input); + eps_eff = sqrt(epsilon3d[0]*epsilon3d[0] + epsilon3d[1]*epsilon3d[1] + epsilon3d[3]*epsilon3d[3] + epsilon3d[4]*epsilon3d[4] + epsilon3d[5]*epsilon3d[5] + epsilon3d[0]*epsilon3d[1]+eps0*eps0); + } + else{ + /* eps_eff^2 = 1/2 ( exx^2 + eyy^2 + 2*exy^2 )*/ + element->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input); + eps_eff = 1./sqrt(2.)*sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + 2.*epsilon2d[2]*epsilon2d[2]); + } + /*Compute viscosity*/ - *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],isdepthaveraged); + *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],eps_eff,isdepthaveraged); } /*}}}*/ @@ -532,5 +549,19 @@ IssmDouble vx,vy,vz; IssmDouble dvx[3],dvy[3],dvz[3]; + IssmDouble epsilon3d[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ + IssmDouble epsilon2d[5]; /* epsilon=[exx,exy];*/ + IssmDouble eps_eff; bool isdepthaveraged=0.; + + if(dim==3){ + /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */ + element->StrainRateHO(&epsilon3d[0],xyz_list,gauss,vx_input,vy_input); + eps_eff = sqrt(epsilon3d[0]*epsilon3d[0] + epsilon3d[1]*epsilon3d[1] + epsilon3d[2]*epsilon3d[2] + epsilon3d[3]*epsilon3d[3] + epsilon3d[4]*epsilon3d[4] + epsilon3d[0]*epsilon3d[1]); + } + else{ + /* eps_eff^2 = 1/2 (2*exx^2 + 2*exy^2 ) (since eps_zz = - eps_xx)*/ + element->StrainRateHO2dvertical(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input); + eps_eff = 1./sqrt(2.)*sqrt(2*epsilon2d[0]*epsilon2d[0] + 2*epsilon2d[1]*epsilon2d[1]); + } /*Get velocity derivatives in all directions*/ @@ -553,5 +584,5 @@ /*Compute viscosity*/ - *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],isdepthaveraged); + *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],eps_eff,isdepthaveraged); }/*}}}*/ void Matestar::ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/ @@ -566,4 +597,7 @@ IssmDouble vx,vy,vz; IssmDouble dvx[3],dvy[3],dvz[3]; + IssmDouble epsilon2d[3];/* epsilon=[exx,eyy,exy]; */ + IssmDouble epsilon1d; /* epsilon=[exx]; */ + IssmDouble eps_eff; bool isdepthaveraged=1.; @@ -589,6 +623,17 @@ dvz[0] = 0.; dvz[1] = 0.; dvz[2] = -dvx[0]-dvy[1]; + if(dim==2){ + /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exx*eyy*/ + element->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input); + eps_eff = sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + epsilon2d[2]*epsilon2d[2] + epsilon2d[0]*epsilon2d[1]); + } + else{ + /* eps_eff^2 = exx^2*/ + element->StrainRateSSA1d(&epsilon1d,xyz_list,gauss,vx_input); + eps_eff = fabs(epsilon1d); + } + /*Compute viscosity*/ - *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],isdepthaveraged); + *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],eps_eff,isdepthaveraged); }/*}}}*/ void Matestar::ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Materials/Matestar.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matestar.h (revision 21825) +++ /issm/trunk-jpl/src/c/classes/Materials/Matestar.h (revision 21826) @@ -86,10 +86,10 @@ void ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); void ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon); - void ViscosityBFS(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input); - void ViscosityBHO(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); - void ViscosityBSSA(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); + void ViscosityBFS(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input,IssmDouble eps_eff); + void ViscosityBHO(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,IssmDouble eps_eff); + void ViscosityBSSA(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,IssmDouble eps_eff); /*}}}*/ - IssmDouble GetViscosityGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,bool isdepthaveraged); - IssmDouble GetViscosity_BGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,bool isdepthaveraged); + IssmDouble GetViscosityGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,IssmDouble eps_eff,bool isdepthaveraged); + IssmDouble GetViscosity_BGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,IssmDouble eps_eff,bool isdepthaveraged); }; Index: /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 21825) +++ /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 21826) @@ -724,5 +724,5 @@ IssmDouble epsilon3d[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/ IssmDouble epsilon2d[2];/* epsilon=[exx,exy]; */ - IssmDouble eps_eff,E=1.0; + IssmDouble eps_eff; if(dim==3){ Index: /issm/trunk-jpl/src/c/classes/Materials/Matice.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matice.h (revision 21825) +++ /issm/trunk-jpl/src/c/classes/Materials/Matice.h (revision 21826) @@ -88,7 +88,7 @@ void ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); void ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon); - void ViscosityBFS(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not supported");}; - void ViscosityBHO(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");}; - void ViscosityBSSA(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");}; + void ViscosityBFS(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input,IssmDouble eps_eff){_error_("not supported");}; + void ViscosityBHO(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,IssmDouble eps_eff){_error_("not supported");}; + void ViscosityBSSA(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,IssmDouble eps_eff){_error_("not supported");}; /*}}}*/ }; Index: /issm/trunk-jpl/src/c/classes/Materials/Matpar.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matpar.h (revision 21825) +++ /issm/trunk-jpl/src/c/classes/Materials/Matpar.h (revision 21826) @@ -125,7 +125,7 @@ void ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");}; void ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){_error_("not supported");}; - void ViscosityBFS(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not supported");}; - void ViscosityBHO(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");}; - void ViscosityBSSA(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");}; + void ViscosityBFS(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input,IssmDouble eps_eff){_error_("not supported");}; + void ViscosityBHO(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,IssmDouble eps_eff){_error_("not supported");}; + void ViscosityBSSA(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,IssmDouble eps_eff){_error_("not supported");}; /*}}}*/ /*Numerics: {{{*/ Index: /issm/trunk-jpl/src/c/shared/Elements/EstarComponents.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/EstarComponents.cpp (revision 21825) +++ /issm/trunk-jpl/src/c/shared/Elements/EstarComponents.cpp (revision 21826) @@ -3,10 +3,10 @@ #include "../Exceptions/exceptions.h" #include "./elements.h" -void EstarStrainrateQuantities(IssmDouble *pepseff, IssmDouble *pepsprime_norm, IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble vmag,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,IssmDouble* dvmag){/*{{{*/ +void EstarStrainrateQuantities(IssmDouble *pepsprime_norm, IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble vmag,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,IssmDouble* dvmag){/*{{{*/ /*Intermediaries*/ IssmDouble omega[3]; IssmDouble nrsp[3],nrsp_norm; - IssmDouble eps[3][3],epseff; + IssmDouble eps[3][3]; IssmDouble epsprime[3],epsprime_norm; @@ -36,9 +36,4 @@ eps[1][0] = .5*(dvx[1]+dvy[0]); eps[1][1] = dvy[1]; eps[1][2] = .5*(dvy[2]+dvz[1]); eps[2][0] = .5*(dvx[2]+dvz[0]); eps[2][1] = .5*(dvy[2]+dvz[1]); eps[2][2] = dvz[2]; - - /*effective strain rate*/ - epseff = 0.; - /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */ - epseff = sqrt(eps[0][0]*eps[0][0] + eps[1][1]*eps[1][1] + eps[0][1]*eps[0][1] + eps[0][2]*eps[0][2] + eps[1][2]*eps[1][2] + eps[0][0]*eps[1][1]); /*Compute the shear strain rate on the non rotating shear plane*/ @@ -73,5 +68,4 @@ /*Assign output pointers*/ - *pepseff=epseff; *pepsprime_norm=epsprime_norm; }/*}}}*/ Index: /issm/trunk-jpl/src/c/shared/Elements/elements.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/elements.h (revision 21825) +++ /issm/trunk-jpl/src/c/shared/Elements/elements.h (revision 21826) @@ -17,5 +17,5 @@ IssmDouble EstarLambdaS(IssmDouble epseff, IssmDouble epsprime_norm); void EstarOmega(IssmDouble* omega,IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble vmag,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz, IssmDouble* dvmag); -void EstarStrainrateQuantities(IssmDouble *pepseff, IssmDouble *pepsprime_norm, IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble vmag,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,IssmDouble* dvmag); +void EstarStrainrateQuantities(IssmDouble *pepsprime_norm, IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble vmag,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,IssmDouble* dvmag); IssmDouble PddSurfaceMassBalance(IssmDouble* monthlytemperatures, IssmDouble* monthlyprec, IssmDouble* pdds, IssmDouble* pds, IssmDouble* melt, IssmDouble* accu, IssmDouble signorm,