Index: /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 20607) +++ /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 20608) @@ -105,5 +105,5 @@ element->StrainRateHO(&epsilon[0],xyz_list,gauss,vx_input,vy_input); - element->ViscosityFSDerivativeEpsSquare(&mu_prime,&epsilon[0]); + element->material->ViscosityFSDerivativeEpsSquare(&mu_prime,&epsilon[0]); eps1[0]=epsilon[0]; eps2[0]=epsilon[2]; eps3[0]=epsilon[3]; eps1[1]=epsilon[2]; eps2[1]=epsilon[1]; eps3[1]=epsilon[4]; @@ -179,5 +179,5 @@ element->StrainRateHO(&epsilon[0],xyz_list,gauss,vx_input,vy_input); - element->ViscosityHODerivativeEpsSquare(&mu_prime,&epsilon[0]); + element->material->ViscosityHODerivativeEpsSquare(&mu_prime,&epsilon[0]); eps1[0]=2.*epsilon[0]+epsilon[1]; eps2[0]=epsilon[2]; eps1[1]=epsilon[2]; eps2[1]=epsilon[0]+2.*epsilon[1]; @@ -283,5 +283,5 @@ thickness_input->GetInputValue(&thickness, gauss); basalelement->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input); - basalelement->ViscositySSADerivativeEpsSquare(&mu_prime,&epsilon[0]); + basalelement->material->ViscositySSADerivativeEpsSquare(&mu_prime,&epsilon[0]); eps1[0]=2.*epsilon[0]+epsilon[1]; eps2[0]=epsilon[2]; eps1[1]=epsilon[2]; eps2[1]=epsilon[0]+2*epsilon[1]; Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 20607) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 20608) @@ -1200,5 +1200,5 @@ thickness_input->GetInputValue(&thickness, gauss); basalelement->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input); - basalelement->ViscositySSADerivativeEpsSquare(&mu_prime,&epsilon[0]); + basalelement->material->ViscositySSADerivativeEpsSquare(&mu_prime,&epsilon[0]); eps1[0]=2*epsilon[0]+epsilon[1]; eps2[0]=epsilon[2]; eps1[1]=epsilon[2]; eps2[1]=epsilon[0]+2*epsilon[1]; @@ -1474,6 +1474,6 @@ this->GetBSSAprime(Bprime,element,dim,xyz_list,gauss); - element->ViscositySSA(&viscosity,dim,xyz_list,gauss,vx_input,vy_input); - element->ViscositySSA(&oldviscosity,dim,xyz_list,gauss,vxold_input,vyold_input); + element->material->ViscositySSA(&viscosity,dim,xyz_list,gauss,vx_input,vy_input); + element->material->ViscositySSA(&oldviscosity,dim,xyz_list,gauss,vxold_input,vyold_input); thickness_input->GetInputValue(&thickness, gauss); @@ -1961,5 +1961,5 @@ this->GetBSSAprime(Bprime,basalelement,2,xyz_list,gauss_base); - element->ViscosityL1L2(&viscosity,xyz_list,gauss,vx_input,vy_input,surface_input); + element->material->ViscosityL1L2(&viscosity,xyz_list,gauss,vx_input,vy_input,surface_input); for(int i=0;i<3;i++) D[i*3+i]=2*viscosity*gauss->weight*Jdet; @@ -2257,5 +2257,5 @@ element->StrainRateHO(&epsilon[0],xyz_list,gauss,vx_input,vy_input); - element->ViscosityHODerivativeEpsSquare(&mu_prime,&epsilon[0]); + element->material->ViscosityHODerivativeEpsSquare(&mu_prime,&epsilon[0]); eps1[0]=2*epsilon[0]+epsilon[1]; eps2[0]=epsilon[2]; eps1[1]=epsilon[2]; eps2[1]=epsilon[0]+2*epsilon[1]; @@ -2428,6 +2428,6 @@ this->GetBHOprime(Bprime,element,dim,xyz_list,gauss); - element->ViscosityHO(&viscosity,dim,xyz_list,gauss,vx_input,vy_input); - element->ViscosityHO(&oldviscosity,dim,xyz_list,gauss,vxold_input,vyold_input); + element->material->ViscosityHO(&viscosity,dim,xyz_list,gauss,vx_input,vy_input); + element->material->ViscosityHO(&oldviscosity,dim,xyz_list,gauss,vxold_input,vyold_input); newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); @@ -2937,5 +2937,5 @@ eps1[1]=epsilon[2]; eps2[1]=epsilon[1]; eps3[1]=epsilon[4]; eps1[2]=epsilon[3]; eps2[2]=epsilon[4]; eps3[2]= -epsilon[0] -epsilon[1]; - element->ViscosityFSDerivativeEpsSquare(&mu_prime,&epsilon[0]); + element->material->ViscosityFSDerivativeEpsSquare(&mu_prime,&epsilon[0]); for(i=0;iGetBFSprime(Bprime,element,dim,xyz_list,gauss); - element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input); + element->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input); for(i=0;iweight*Jdet; for(i=epssize;iweight*Jdet; @@ -3191,5 +3191,5 @@ this->GetBFSprimevel(Bprime,element,dim,xyz_list,gauss); - element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input); + element->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input); for(i=0;iweight*Jdet; @@ -4881,5 +4881,5 @@ } - element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input); + element->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input); epsxx[i] = dvx[0]; sigmapxx[i] = 2.*viscosity*epsxx[i]; epsyy[i] = dvy[1]; sigmapyy[i] = 2.*viscosity*epsyy[i]; @@ -5508,5 +5508,5 @@ this->GetLprimeFSSSA(&LprimeFSSSA[0][0], element,xyz_list, gauss); - element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input); + element->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input); element->NormalBase(&bed_normal[0],xyz_list_tria); @@ -5625,5 +5625,5 @@ this->GetBprimeSSAFS(&Bprime2[0][0], element,xyz_list, gauss); - element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input); + element->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input); D_scalar=2*viscosity*gauss->weight*Jdet; @@ -5805,6 +5805,6 @@ this->GetBSSAHO(B, element,xyz_list, gauss); this->GetBSSAprime(Bprime,basaltria,2,xyz_list, gauss_tria); - element->ViscosityHO(&viscosity,3,xyz_list,gauss,vx_input,vy_input); - element->ViscosityHO(&oldviscosity,3,xyz_list,gauss,vxold_input,vyold_input); + element->material->ViscosityHO(&viscosity,3,xyz_list,gauss,vx_input,vy_input); + element->material->ViscosityHO(&oldviscosity,3,xyz_list,gauss,vxold_input,vyold_input); newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); @@ -5964,10 +5964,10 @@ if(approximation==SSAHOApproximationEnum){ - element->ViscosityHO(&viscosity,dim,xyz_list,gauss,vx_input,vy_input); - element->ViscosityHO(&oldviscosity,dim,xyz_list,gauss,vxold_input,vyold_input); + element->material->ViscosityHO(&viscosity,dim,xyz_list,gauss,vx_input,vy_input); + element->material->ViscosityHO(&oldviscosity,dim,xyz_list,gauss,vxold_input,vyold_input); newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); } else if (approximation==SSAFSApproximationEnum){ - element->ViscosityFS(&newviscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input); + element->material->ViscosityFS(&newviscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input); } else _error_("approximation " << approximation << " (" << EnumToStringx(approximation) << ") not supported yet"); @@ -6069,5 +6069,5 @@ element->NormalBase(&bed_normal[0],xyz_list_tria); - element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input); + element->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input); friction->GetAlpha2(&alpha2_gauss,gauss); @@ -6139,5 +6139,5 @@ element->NodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss); - element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input); + element->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input); vzHO_input->GetInputDerivativeValue(&dw[0],xyz_list,gauss); @@ -6230,5 +6230,5 @@ element->NormalBase(&bed_normal[0],xyz_list_tria); - element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input); + element->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input); friction->GetAlpha2(&alpha2_gauss,gauss); @@ -6299,5 +6299,5 @@ vzSSA_input->GetInputDerivativeValue(&dw[0],xyz_list,gauss); - element->ViscosityFS(&viscosity,3,xyz_list,gauss,vx_input,vy_input,vz_input); + element->material->ViscosityFS(&viscosity,3,xyz_list,gauss,vx_input,vy_input,vz_input); for(i=0;i<6;i++){ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 20607) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 20608) @@ -2879,161 +2879,4 @@ xDelete(values); }/*}}}*/ -void Element::ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/ - /*The effective strain rate is defined in Paterson 3d Ed p 91 eq 9, - * and Cuffey p 303 eq 8.18: - * - * 2 eps_eff^2 = eps_xx^2 + eps_yy^2 + eps_zz^2 + 2(eps_xy^2 + eps_xz^2 + eps_yz^2) - * - * or - * - * eps_eff = 1/sqrt(2) sqrt( \sum_ij eps_ij^2 ) - * - * = 1/sqrt(2) ||eps||_F - * - * where ||.||_F is the Frobenius norm */ - - /*Intermediaries*/ - IssmDouble viscosity; - IssmDouble epsilon3d[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ - IssmDouble epsilon2d[3]; /* epsilon=[exx,eyy,exy]; */ - IssmDouble eps_eff; - IssmDouble eps0=1.e-27; - - if(dim==3){ - /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */ - this->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 )*/ - this->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]); - } - - /*Get viscosity*/ - material->GetViscosity(&viscosity,eps_eff); - - /*Assign output pointer*/ - *pviscosity=viscosity; -} -/*}}}*/ -void Element::ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/ - this->material->GetViscosityDerivativeEpsSquare(pmu_prime,epsilon); -}/*}}}*/ -void Element::ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ - - /*Intermediaries*/ - IssmDouble viscosity; - IssmDouble epsilon3d[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/ - IssmDouble epsilon2d[2];/* epsilon=[exx,exy]; */ - IssmDouble eps_eff; - - if(dim==3){ - /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */ - this->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 (exx^2 + 2*exy^2 )*/ - this->StrainRateHO2dvertical(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input); - eps_eff = 1./sqrt(2.)*sqrt(epsilon2d[0]*epsilon2d[0] + 2.*epsilon2d[1]*epsilon2d[1]); - } - - /*Get viscosity*/ - material->GetViscosity(&viscosity,eps_eff); - - /*Assign output pointer*/ - *pviscosity=viscosity; -}/*}}}*/ -void Element::ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/ - this->material->GetViscosityDerivativeEpsSquare(pmu_prime,epsilon); -}/*}}}*/ -void Element::ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surface_input){/*{{{*/ - /*Compute the L1L2 viscosity - * - * 1 - * mu = - A^-1 (sigma'_e)^(1-n) - * 2 - * - * sigma'_e^2 = |sigma'_//|^2 + |sigma'_perp|^2 (see Perego 2012 eq. 17,18) - * - * L1L2 assumptions: - * - * (1) |eps_b|_// = A (|sigma'_//|^2 + |sigma'_perp|^2)^((n-1)/2) |sigma'_//| - * (2) |sigma'_perp|^2 = |rho g (s-z) grad(s)|^2 - * - * Assuming that n = 3, we have a polynom of degree 3 to solve (the only unkown is X=|sigma'_//|) - * - * A X^3 + A |rho g (s-z) grad(s)|^2 X - |eps_b|_// = 0 */ - - IssmDouble z,s,viscosity,p,q,delta; - IssmDouble tau_perp,tau_par,eps_b,A; - IssmDouble epsilon[5]; /*exx eyy exy exz eyz*/ - IssmDouble slope[3]; - - /*Check that both inputs have been found*/ - if (!vx_input || !vy_input || !surface_input) _error_("Input missing"); - - /*Get tau_perp*/ - surface_input->GetInputValue(&s,gauss); - surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss); - z=this->GetZcoord(xyz_list,gauss); - tau_perp = matpar->GetMaterialParameter(MaterialsRhoIceEnum) * matpar->GetMaterialParameter(ConstantsGEnum) * fabs(s-z)*sqrt(slope[0]*slope[0]+slope[1]*slope[1]); - - /* Get eps_b*/ - this->StrainRateHO(&epsilon[0],xyz_list,gauss,vx_input,vy_input); - eps_b = sqrt(epsilon[0]*epsilon[0] + epsilon[1]*epsilon[1] + epsilon[0]*epsilon[1] + epsilon[2]*epsilon[2]); - if(eps_b==0.){ - *pviscosity = 2.5e+17; - return; - } - - /*Get A*/ - _assert_(material->GetN()==3.0); - A=material->GetA(); - - /*Solve for tau_perp (http://fr.wikipedia.org/wiki/Méthode_de_Cardan)*/ - p = tau_perp *tau_perp; - q = - eps_b/A; - delta = q *q + p*p*p*4./27.; - _assert_(delta>0); - tau_par = pow(0.5*(-q+sqrt(delta)),1./3.) - pow(0.5*(q+sqrt(delta)),1./3.); - - /*Viscosity*/ - viscosity = 1./(2.*A*(tau_par*tau_par + tau_perp*tau_perp)); - _assert_(!xIsNan(viscosity)); - _assert_(viscosity > 0.); - - /*Assign output pointer*/ - *pviscosity = viscosity; -}/*}}}*/ -void Element::ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ - - /*Intermediaries*/ - IssmDouble viscosity; - IssmDouble epsilon2d[3];/* epsilon=[exx,eyy,exy]; */ - IssmDouble epsilon1d; /* epsilon=[exx]; */ - IssmDouble eps_eff; - - if(dim==2){ - /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exx*eyy*/ - this->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 = 1/2 exx^2*/ - this->StrainRateSSA1d(&epsilon1d,xyz_list,gauss,vx_input); - eps_eff = sqrt(epsilon1d*epsilon1d/2.); - } - - /*Get viscosity*/ - material->GetViscosityBar(&viscosity,eps_eff); - - /*Assign output pointer*/ - *pviscosity=viscosity; -}/*}}}*/ -void Element::ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/ - this->material->GetViscosity2dDerivativeEpsSquare(pmu_prime,epsilon); -}/*}}}*/ void Element::ViscousHeatingCreateInput(void){/*{{{*/ @@ -3064,5 +2907,5 @@ this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input); - this->ViscosityFS(&viscosity,3,xyz_list,gauss,vx_input,vy_input,vz_input); + this->material->ViscosityFS(&viscosity,3,xyz_list,gauss,vx_input,vy_input,vz_input); this->GetPhi(&phi,&epsilon[0],viscosity); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 20607) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 20608) @@ -167,11 +167,4 @@ void TransformStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes,int numnodes,int* cs_array); void TransformStiffnessMatrixCoord(ElementMatrix* Ke,int numnodes,int* transformenum_list){_error_("not implemented yet");};/*Tiling only*/ - void ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input); - void ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon); - void ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); - void ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon); - void ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf); - 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 ViscousHeatingCreateInput(void); Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 20607) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 20608) @@ -289,5 +289,5 @@ /*Compute strain rate viscosity and pressure: */ this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); - this->ViscosityFS(&viscosity,3,&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); + this->material->ViscosityFS(&viscosity,3,&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); pressure_input->GetInputValue(&pressure,gauss); @@ -346,5 +346,5 @@ /*Compute strain rate viscosity and pressure: */ this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); - this->ViscosityFS(&viscosity,3,&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); + this->material->ViscosityFS(&viscosity,3,&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); /*Compute Stress*/ @@ -398,5 +398,5 @@ /*Compute strain rate viscosity and pressure: */ this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); - this->ViscosityFS(&viscosity,3,&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); + this->material->ViscosityFS(&viscosity,3,&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); pressure_input->GetInputValue(&pressure,gauss); @@ -640,5 +640,5 @@ /*Check for basal force only if grounded and touching GL*/ this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input); - this->ViscosityFS(&viscosity,3,xyz_list,gauss,vx_input,vy_input,vz_input); + this->material->ViscosityFS(&viscosity,3,xyz_list,gauss,vx_input,vy_input,vz_input); pressure_input->GetInputValue(&pressure, gauss); base_input->GetInputValue(&base, gauss); _assert_(base<0.); @@ -3174,5 +3174,5 @@ _assert_(gauss->Enum()==GaussPentaEnum); this->StrainRateFS(&epsilon[0],xyz_list,(GaussPenta*)gauss,vx_input,vy_input,vz_input); - this->ViscosityFS(&viscosity,3,xyz_list,(GaussPenta*)gauss,vx_input,vy_input,vz_input); + this->material->ViscosityFS(&viscosity,3,xyz_list,(GaussPenta*)gauss,vx_input,vy_input,vz_input); GetPhi(&phi,&epsilon[0],viscosity); Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.cpp (revision 20607) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.cpp (revision 20608) @@ -997,5 +997,5 @@ _assert_(gauss->Enum()==GaussTetraEnum); this->StrainRateFS(&epsilon[0],xyz_list,(GaussTetra*)gauss,vx_input,vy_input,vz_input); - this->ViscosityFS(&viscosity,3,xyz_list,(GaussTetra*)gauss,vx_input,vy_input,vz_input); + this->material->ViscosityFS(&viscosity,3,xyz_list,(GaussTetra*)gauss,vx_input,vy_input,vz_input); GetPhi(&phi,&epsilon[0],viscosity); Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 20607) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 20608) @@ -498,5 +498,5 @@ /*Compute strain rate and viscosity: */ this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); - this->ViscositySSA(&viscosity,dim,&xyz_list[0][0],gauss,vx_input,vy_input); + this->material->ViscositySSA(&viscosity,dim,&xyz_list[0][0],gauss,vx_input,vy_input); /*Compute Stress*/ @@ -555,5 +555,5 @@ /*Compute strain rate viscosity and pressure: */ this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input); - this->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,NULL); + this->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,NULL); pressure_input->GetInputValue(&pressure,gauss); @@ -611,5 +611,5 @@ /*Compute strain rate viscosity and pressure: */ this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); - this->ViscositySSA(&viscosity,dim,&xyz_list[0][0],gauss,vx_input,vy_input); + this->material->ViscositySSA(&viscosity,dim,&xyz_list[0][0],gauss,vx_input,vy_input); pressure_input->GetInputValue(&pressure,gauss); @@ -906,5 +906,5 @@ // if(this->inputs->Min(MaskGroundediceLevelsetEnum)==0.){ this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input); - this->ViscosityFS(&viscosity,2,xyz_list,gauss,vx_input,vy_input,NULL); + this->material->ViscosityFS(&viscosity,2,xyz_list,gauss,vx_input,vy_input,NULL); pressure_input->GetInputValue(&pressure, gauss); base_input->GetInputValue(&base, gauss); Index: /issm/trunk-jpl/src/c/classes/Materials/Material.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Material.h (revision 20607) +++ /issm/trunk-jpl/src/c/classes/Materials/Material.h (revision 20608) @@ -14,4 +14,6 @@ class Element; class Elements; +class Gauss; +class Input; /*}}}*/ @@ -43,4 +45,12 @@ virtual void ResetHooks()=0; + virtual void ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input)=0; + virtual void ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon)=0; + virtual void ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0; + virtual void ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon)=0; + virtual void ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf)=0; + 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; + }; #endif Index: /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 20607) +++ /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 20608) @@ -337,6 +337,5 @@ 2 eps_eff ^[(n-1)/n] - where viscosity is the viscotiy, B the flow law parameter , eps_eff is the effective strain rate - and n the flow law exponent. + where B the flow law parameter, eps_eff is the effective strain rate and n the flow law exponent. If eps_eff = 0 , it means this is the first time SystemMatrices is being run, and we @@ -580,4 +579,161 @@ } /*}}}*/ +void Matice::ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/ + /*The effective strain rate is defined in Paterson 3d Ed p 91 eq 9, + * and Cuffey p 303 eq 8.18: + * + * 2 eps_eff^2 = eps_xx^2 + eps_yy^2 + eps_zz^2 + 2(eps_xy^2 + eps_xz^2 + eps_yz^2) + * + * or + * + * eps_eff = 1/sqrt(2) sqrt( \sum_ij eps_ij^2 ) + * + * = 1/sqrt(2) ||eps||_F + * + * where ||.||_F is the Frobenius norm */ + + /*Intermediaries*/ + IssmDouble viscosity; + IssmDouble epsilon3d[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ + IssmDouble epsilon2d[3]; /* epsilon=[exx,eyy,exy]; */ + IssmDouble eps_eff; + IssmDouble eps0=1.e-27; + + 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]); + } + + /*Get viscosity*/ + this->GetViscosity(&viscosity,eps_eff); + + /*Assign output pointer*/ + *pviscosity=viscosity; +} +/*}}}*/ +void Matice::ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/ + this->GetViscosityDerivativeEpsSquare(pmu_prime,epsilon); +}/*}}}*/ +void Matice::ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ + + /*Intermediaries*/ + IssmDouble viscosity; + IssmDouble epsilon3d[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/ + IssmDouble epsilon2d[2];/* epsilon=[exx,exy]; */ + IssmDouble eps_eff; + + 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 (exx^2 + 2*exy^2 )*/ + element->StrainRateHO2dvertical(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input); + eps_eff = 1./sqrt(2.)*sqrt(epsilon2d[0]*epsilon2d[0] + 2.*epsilon2d[1]*epsilon2d[1]); + } + + /*Get viscosity*/ + this->GetViscosity(&viscosity,eps_eff); + + /*Assign output pointer*/ + *pviscosity=viscosity; +}/*}}}*/ +void Matice::ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/ + this->GetViscosityDerivativeEpsSquare(pmu_prime,epsilon); +}/*}}}*/ +void Matice::ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surface_input){/*{{{*/ + /*Compute the L1L2 viscosity + * + * 1 + * mu = - A^-1 (sigma'_e)^(1-n) + * 2 + * + * sigma'_e^2 = |sigma'_//|^2 + |sigma'_perp|^2 (see Perego 2012 eq. 17,18) + * + * L1L2 assumptions: + * + * (1) |eps_b|_// = A (|sigma'_//|^2 + |sigma'_perp|^2)^((n-1)/2) |sigma'_//| + * (2) |sigma'_perp|^2 = |rho g (s-z) grad(s)|^2 + * + * Assuming that n = 3, we have a polynom of degree 3 to solve (the only unkown is X=|sigma'_//|) + * + * A X^3 + A |rho g (s-z) grad(s)|^2 X - |eps_b|_// = 0 */ + + IssmDouble z,s,viscosity,p,q,delta; + IssmDouble tau_perp,tau_par,eps_b,A; + IssmDouble epsilon[5]; /*exx eyy exy exz eyz*/ + IssmDouble slope[3]; + + /*Check that both inputs have been found*/ + if (!vx_input || !vy_input || !surface_input) _error_("Input missing"); + + /*Get tau_perp*/ + surface_input->GetInputValue(&s,gauss); + surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss); + z=this->element->GetZcoord(xyz_list,gauss); + tau_perp = element->matpar->GetMaterialParameter(MaterialsRhoIceEnum) * element->matpar->GetMaterialParameter(ConstantsGEnum) * fabs(s-z)*sqrt(slope[0]*slope[0]+slope[1]*slope[1]); + + /* Get eps_b*/ + element->StrainRateHO(&epsilon[0],xyz_list,gauss,vx_input,vy_input); + eps_b = sqrt(epsilon[0]*epsilon[0] + epsilon[1]*epsilon[1] + epsilon[0]*epsilon[1] + epsilon[2]*epsilon[2]); + if(eps_b==0.){ + *pviscosity = 2.5e+17; + return; + } + + /*Get A*/ + _assert_(this->GetN()==3.0); + A=this->GetA(); + + /*Solve for tau_perp (http://fr.wikipedia.org/wiki/Méthode_de_Cardan)*/ + p = tau_perp *tau_perp; + q = - eps_b/A; + delta = q *q + p*p*p*4./27.; + _assert_(delta>0); + tau_par = pow(0.5*(-q+sqrt(delta)),1./3.) - pow(0.5*(q+sqrt(delta)),1./3.); + + /*Viscosity*/ + viscosity = 1./(2.*A*(tau_par*tau_par + tau_perp*tau_perp)); + _assert_(!xIsNan(viscosity)); + _assert_(viscosity > 0.); + + /*Assign output pointer*/ + *pviscosity = viscosity; +}/*}}}*/ +void Matice::ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ + + /*Intermediaries*/ + IssmDouble viscosity; + IssmDouble epsilon2d[3];/* epsilon=[exx,eyy,exy]; */ + IssmDouble epsilon1d; /* epsilon=[exx]; */ + IssmDouble eps_eff; + + 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 = 1/2 exx^2*/ + element->StrainRateSSA1d(&epsilon1d,xyz_list,gauss,vx_input); + eps_eff = sqrt(epsilon1d*epsilon1d/2.); + } + + /*Get viscosity*/ + this->GetViscosityBar(&viscosity,eps_eff); + + /*Assign output pointer*/ + *pviscosity=viscosity; +}/*}}}*/ +void Matice::ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/ + this->GetViscosity2dDerivativeEpsSquare(pmu_prime,epsilon); +}/*}}}*/ void Matice::ResetHooks(){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Materials/Matice.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matice.h (revision 20607) +++ /issm/trunk-jpl/src/c/classes/Materials/Matice.h (revision 20608) @@ -18,4 +18,6 @@ class Materials; class Parameters; +class Gauss; +class Input; /*}}}*/ @@ -72,4 +74,12 @@ bool IsDamage(); void ResetHooks(); + + void ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input); + void ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon); + void ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); + void ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon); + void ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf); + void ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); + void ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon); /*}}}*/ }; Index: /issm/trunk-jpl/src/c/classes/Materials/Matpar.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matpar.h (revision 20607) +++ /issm/trunk-jpl/src/c/classes/Materials/Matpar.h (revision 20608) @@ -115,4 +115,12 @@ bool IsDamage(){_error_("not supported");}; void ResetHooks(); + + void ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not supported");}; + void ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){_error_("not supported");}; + void ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");}; + void ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){_error_("not supported");}; + void ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf){_error_("not supported");}; + 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");}; /*}}}*/ /*Numerics: {{{*/