Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 16914) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 16915) @@ -108,2 +108,190 @@ return this->matpar->TMeltingPoint(pressure); }/*}}}*/ + +void Element::ViscosityFS(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/ + + /*Intermediaries*/ + IssmDouble viscosity; + IssmDouble epsilon3d[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ + IssmDouble epsilon2d[3]; /* epsilon=[exx,eyy,exy]; */ + + if(vz_input){ + /*3D*/ + this->StrainRateFS(&epsilon3d[0],xyz_list,gauss,vx_input,vy_input,vz_input); + material->GetViscosity3dFS(&viscosity, &epsilon3d[0]); + } + else{ + /*2D*/ + this->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input); + material->GetViscosity2dvertical(&viscosity,&epsilon2d[0]); + } + + /*Assign output pointer*/ + *pviscosity=viscosity; +} +/*}}}*/ +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 */ + + int i; + IssmDouble z,s,viscosity,p,q,delta; + IssmDouble tau_perp,tau_par,eps_b,A; + IssmDouble epsilonvx[5]; /*exx eyy exy exz eyz*/ + IssmDouble epsilonvy[5]; /*exx eyy exy exz eyz*/ + 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=GetZcoord(gauss); + tau_perp = matpar->GetRhoIce() * matpar->GetG() * fabs(s-z)*sqrt(slope[0]*slope[0]+slope[1]*slope[1]); + + /* Get eps_b*/ + vx_input->GetVxStrainRate3dHO(epsilonvx,xyz_list,gauss); + vy_input->GetVyStrainRate3dHO(epsilonvy,xyz_list,gauss); + for(i=0;i<5;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]; + 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::ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ + + /*Intermediaries*/ + IssmDouble viscosity; + IssmDouble epsilon[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/ + + this->StrainRateHO(&epsilon[0],xyz_list,gauss,vx_input,vy_input); + material->GetViscosity3d(&viscosity, &epsilon[0]); + + /*Assign output pointer*/ + *pviscosity=viscosity; +}/*}}}*/ +void Element::ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/ + this->material->GetViscosityDerivativeEpsSquare(pmu_prime,epsilon); +}/*}}}*/ +void Element::StrainRateHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ + /*Compute the 3d Blatter/HOStrain Rate (5 components): + * + * epsilon=[exx eyy exy exz eyz] + * + * with exz=1/2 du/dz + * eyz=1/2 dv/dz + * + * the contribution of vz is neglected + */ + + int i; + IssmDouble epsilonvx[5]; + IssmDouble epsilonvy[5]; + + /*Check that both inputs have been found*/ + if (!vx_input || !vy_input){ + _error_("Input missing. Here are the input pointers we have for vx: " << vx_input << ", vy: " << vy_input << "\n"); + } + + /*Get strain rate assuming that epsilon has been allocated*/ + vx_input->GetVxStrainRate3dHO(epsilonvx,xyz_list,gauss); + vy_input->GetVyStrainRate3dHO(epsilonvy,xyz_list,gauss); + + /*Sum all contributions*/ + for(i=0;i<5;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]; +}/*}}}*/ +void Element::ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/ + this->material->GetViscosityDerivativeEpsSquare(pmu_prime,epsilon); +}/*}}}*/ +void Element::StrainRateFS(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/ + /*Compute the 3d Strain Rate (6 components): + * + * epsilon=[exx eyy ezz exy exz eyz] + */ + + IssmDouble epsilonvx[6]; + IssmDouble epsilonvy[6]; + IssmDouble epsilonvz[6]; + + /*Check that both inputs have been found*/ + if (!vx_input || !vy_input || !vz_input){ + _error_("Input missing. Here are the input pointers we have for vx: " << vx_input << ", vy: " << vy_input << ", vz: " << vz_input << "\n"); + } + + /*Get strain rate assuming that epsilon has been allocated*/ + vx_input->GetVxStrainRate3d(epsilonvx,xyz_list,gauss); + vy_input->GetVyStrainRate3d(epsilonvy,xyz_list,gauss); + vz_input->GetVzStrainRate3d(epsilonvz,xyz_list,gauss); + + /*Sum all contributions*/ + for(int i=0;i<6;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]+epsilonvz[i]; +}/*}}}*/ +void Element::ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ + + /*Intermediaries*/ + IssmDouble viscosity; + IssmDouble epsilon[3];/* epsilon=[exx,eyy,exy]; */ + + this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input); + material->GetViscosity2d(&viscosity, &epsilon[0]); + + /*Assign output pointer*/ + *pviscosity=viscosity; +}/*}}}*/ +void Element::ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/ + this->material->GetViscosity2dDerivativeEpsSquare(pmu_prime,epsilon); +}/*}}}*/ +void Element::StrainRateSSA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ + + int i; + IssmDouble epsilonvx[3]; + IssmDouble epsilonvy[3]; + + /*Check that both inputs have been found*/ + if (!vx_input || !vy_input){ + _error_("Input missing. Here are the input pointers we have for vx: " << vx_input << ", vy: " << vy_input << "\n"); + } + + /*Get strain rate assuming that epsilon has been allocated*/ + vx_input->GetVxStrainRate2d(epsilonvx,xyz_list,gauss); + vy_input->GetVyStrainRate2d(epsilonvy,xyz_list,gauss); + + /*Sum all contributions*/ + for(i=0;i<3;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]; +}/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16914) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16915) @@ -62,5 +62,15 @@ IssmDouble GetMaterialParameter(int enum_in); bool IsFloating(); + void StrainRateSSA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); + void StrainRateHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); + void StrainRateFS(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input); IssmDouble TMeltingPoint(IssmDouble pressure); + void ViscosityFS(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input); + void ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); + void ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf); + void ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); + void ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon); + void ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon); + void ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon); /*Virtual functions*/ @@ -186,14 +196,5 @@ virtual void ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss)=0; virtual void ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input)=0; - virtual void ViscosityFS(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input)=0; - virtual void ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=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,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0; - virtual void ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon)=0; - virtual void ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon)=0; - virtual void ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon)=0; - virtual void StrainRateSSA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0; - virtual void StrainRateHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0; - virtual void StrainRateFS(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input)=0; + virtual int VelocityInterpolation()=0; virtual int PressureInterpolation()=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16914) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16915) @@ -3452,164 +3452,4 @@ /*Assign output pointer*/ *pphi = phi; -} -/*}}}*/ -/*FUNCTION Penta::ViscosityFS{{{*/ -void Penta::ViscosityFS(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){ - - /*Intermediaries*/ - IssmDouble viscosity; - IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ - - this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input); - material->GetViscosity3dFS(&viscosity, &epsilon[0]); - - /*Assign output pointer*/ - *pviscosity=viscosity; -} -/*}}}*/ -/*FUNCTION Penta::ViscosityL1L2{{{*/ -void Penta::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 */ - - int i; - IssmDouble z,s,viscosity,p,q,delta; - IssmDouble tau_perp,tau_par,eps_b,A; - IssmDouble epsilonvx[5]; /*exx eyy exy exz eyz*/ - IssmDouble epsilonvy[5]; /*exx eyy exy exz eyz*/ - IssmDouble epsilon[5]; /*exx eyy exy exz eyz*/ - IssmDouble z_list[NUMVERTICES]; - IssmDouble slope[3]; - - /*Check that both inputs have been found*/ - if (!vx_input || !vy_input || !surface_input) _error_("Input missing"); - - /*Get tau_perp*/ - for(i=0;iGetInputValue(&s,gauss); - surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss); - PentaRef::GetInputValue(&z,&z_list[0],gauss); - tau_perp = matpar->GetRhoIce() * matpar->GetG() * fabs(s-z)*sqrt(slope[0]*slope[0]+slope[1]*slope[1]); - - /* Get eps_b*/ - vx_input->GetVxStrainRate3dHO(epsilonvx,xyz_list,gauss); - vy_input->GetVyStrainRate3dHO(epsilonvy,xyz_list,gauss); - for(i=0;i<5;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]; - 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; -} -/*}}}*/ -/*FUNCTION Penta::ViscosityHO{{{*/ -void Penta::ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){ - - /*Intermediaries*/ - IssmDouble viscosity; - IssmDouble epsilon[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/ - - this->StrainRateHO(&epsilon[0],xyz_list,gauss,vx_input,vy_input); - material->GetViscosity3d(&viscosity, &epsilon[0]); - - /*Assign output pointer*/ - *pviscosity=viscosity; -} -/*}}}*/ -/*FUNCTION Penta::ViscosityHODerivativeEpsSquare{{{*/ -void Penta::ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){ - this->material->GetViscosityDerivativeEpsSquare(pmu_prime,epsilon); -} -/*}}}*/ -/*FUNCTION Penta::StrainRateHO{{{*/ -void Penta::StrainRateHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){ - /*Compute the 3d Blatter/HOStrain Rate (5 components): - * - * epsilon=[exx eyy exy exz eyz] - * - * with exz=1/2 du/dz - * eyz=1/2 dv/dz - * - * the contribution of vz is neglected - */ - - int i; - IssmDouble epsilonvx[5]; - IssmDouble epsilonvy[5]; - - /*Check that both inputs have been found*/ - if (!vx_input || !vy_input){ - _error_("Input missing. Here are the input pointers we have for vx: " << vx_input << ", vy: " << vy_input << "\n"); - } - - /*Get strain rate assuming that epsilon has been allocated*/ - vx_input->GetVxStrainRate3dHO(epsilonvx,xyz_list,gauss); - vy_input->GetVyStrainRate3dHO(epsilonvy,xyz_list,gauss); - - /*Sum all contributions*/ - for(i=0;i<5;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]; -} -/*}}}*/ -/*FUNCTION Penta::ViscosityFSDerivativeEpsSquare{{{*/ -void Penta::ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){ - this->material->GetViscosityDerivativeEpsSquare(pmu_prime,epsilon); -} -/*}}}*/ -/*FUNCTION Penta::StrainRateFS{{{*/ -void Penta::StrainRateFS(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){ - /*Compute the 3d Strain Rate (6 components): - * - * epsilon=[exx eyy ezz exy exz eyz] - */ - - IssmDouble epsilonvx[6]; - IssmDouble epsilonvy[6]; - IssmDouble epsilonvz[6]; - - /*Check that both inputs have been found*/ - if (!vx_input || !vy_input || !vz_input){ - _error_("Input missing. Here are the input pointers we have for vx: " << vx_input << ", vy: " << vy_input << ", vz: " << vz_input << "\n"); - } - - /*Get strain rate assuming that epsilon has been allocated*/ - vx_input->GetVxStrainRate3d(epsilonvx,xyz_list,gauss); - vy_input->GetVyStrainRate3d(epsilonvy,xyz_list,gauss); - vz_input->GetVzStrainRate3d(epsilonvz,xyz_list,gauss); - - /*Sum all contributions*/ - for(int i=0;i<6;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]+epsilonvz[i]; } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16914) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16915) @@ -278,14 +278,4 @@ void TransformStiffnessMatrixCoord(ElementMatrix* Ke,int numnodes,int* transformenum_list){_error_("not implemented yet");}; void TransformStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes_list,int numnodes,int* transformenum_list); - void ViscosityFS(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input); - void ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); - void ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf); - void ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented");}; - void ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon); - void ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon); - void ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){_error_("not implemented");}; - void StrainRateFS(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input); - void StrainRateHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); - void StrainRateSSA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented");}; #ifdef _HAVE_STRESSBALANCE_ Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16914) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16915) @@ -152,14 +152,4 @@ void VerticalSegmentIndices(int** pindices,int* pnumseg){_error_("not implemented yet");}; void ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not implemented yet");}; - void ViscosityFS(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not implemented");}; - void ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented");}; - void ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf){_error_("not implemented yet");}; - void ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented");}; - void ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){_error_("not implemented");}; - void ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){_error_("not implemented");}; - void ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){_error_("not implemented");}; - void StrainRateFS(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not implemented");}; - void StrainRateHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented");}; - void StrainRateSSA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented");}; bool IsZeroLevelset(int levelset_enum){_error_("not implemented");}; void ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum){_error_("not implemented");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16914) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16915) @@ -3000,57 +3000,4 @@ _assert_(this->vertices); return this->vertices[vertexindex]->Connectivity(); -} -/*}}}*/ -/*FUNCTION Tria::ViscosityFS{{{*/ -void Tria::ViscosityFS(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){ - - /*Intermediaries*/ - IssmDouble viscosity; - IssmDouble epsilon[3];/* epsilon=[exx,eyy,exy]; */ - - this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input); - material->GetViscosity2dvertical(&viscosity, &epsilon[0]); - - /*Assign output pointer*/ - *pviscosity=viscosity; -} -/*}}}*/ -/*FUNCTION Tria::ViscositySSA{{{*/ -void Tria::ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){ - - /*Intermediaries*/ - IssmDouble viscosity; - IssmDouble epsilon[3];/* epsilon=[exx,eyy,exy]; */ - - this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input); - material->GetViscosity2d(&viscosity, &epsilon[0]); - - /*Assign output pointer*/ - *pviscosity=viscosity; -} -/*}}}*/ -/*FUNCTION Tria::ViscositySSADerivativeEpsSquare{{{*/ -void Tria::ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){ - this->material->GetViscosity2dDerivativeEpsSquare(pmu_prime,epsilon); -} -/*}}}*/ -/*FUNCTION Tria::StrainRateSSA{{{*/ -void Tria::StrainRateSSA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){ - - int i; - IssmDouble epsilonvx[3]; - IssmDouble epsilonvy[3]; - - /*Check that both inputs have been found*/ - if (!vx_input || !vy_input){ - _error_("Input missing. Here are the input pointers we have for vx: " << vx_input << ", vy: " << vy_input << "\n"); - } - - /*Get strain rate assuming that epsilon has been allocated*/ - vx_input->GetVxStrainRate2d(epsilonvx,xyz_list,gauss); - vy_input->GetVyStrainRate2d(epsilonvy,xyz_list,gauss); - - /*Sum all contributions*/ - for(i=0;i<3;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]; } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16914) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16915) @@ -309,14 +309,4 @@ void TransformStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes_list,int numnodes,int* transformenum_list){_error_("not implemented yet");}; void ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not implemented yet");}; - void ViscosityFS(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input); - void ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented yet");}; - void ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf){_error_("not implemented yet");}; - void ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); - void ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){_error_("not implemented");}; - void ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){_error_("not implemented");}; - void ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon); - void StrainRateFS(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not implemented");}; - void StrainRateHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented");}; - void StrainRateSSA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); #ifdef _HAVE_STRESSBALANCE_