Index: /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 16875) +++ /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 16876) @@ -27,5 +27,108 @@ /*Finite Element Analysis*/ ElementMatrix* AdjointHorizAnalysis::CreateKMatrix(Element* element){/*{{{*/ - _error_("not implemented yet"); + int approximation; + element->GetInputValue(&approximation,ApproximationEnum); + switch(approximation){ + case SSAApproximationEnum: + return CreateKMatrixSSA(element); + case HOApproximationEnum: + return CreateKMatrixHO(element); + case FSApproximationEnum: + return CreateKMatrixFS(element); + case NoneApproximationEnum: + return NULL; + default: + _error_("Approximation "<FindParam(&meshtype,MeshTypeEnum); + switch(meshtype){ + case Mesh2DhorizontalEnum: + basalelement = element; + break; + case Mesh3DEnum: + if(!element->IsOnBed()) return NULL; + basalelement = element->SpawnBasalElement(); + break; + default: _error_("mesh "<GetNumberOfNodes(); + + /*Initialize Jacobian with regular SSA (first part of the Gateau derivative)*/ + basalelement->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum); + StressbalanceAnalysis* analysis = new StressbalanceAnalysis(); + ElementMatrix* Ke=analysis->CreateKMatrix(element); + delete analysis; + if(incomplete_adjoint) return Ke; + + /*Retrieve all inputs and parameters*/ + basalelement->GetVerticesCoordinates(&xyz_list); + Input* vx_input = basalelement->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = basalelement->GetInput(VyEnum); _assert_(vy_input); + Input* thickness_input = basalelement->GetInput(ThicknessEnum); _assert_(thickness_input); + + /*Allocate dbasis*/ + IssmDouble* dbasis = xNew(2*numnodes); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=basalelement->NewGauss(2); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss); + basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); + + thickness_input->GetInputValue(&thickness, gauss); + basalelement->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input); + basalelement->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]; + + for(int i=0;ivalues[2*numnodes*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps1dotdphii; + Ke->values[2*numnodes*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps1dotdphii; + Ke->values[2*numnodes*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps2dotdphii; + Ke->values[2*numnodes*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps2dotdphii; + } + } + } + + /*Transform Coordinate System*/ + basalelement->TransformStiffnessMatrixCoord(Ke,XYEnum); + + /*Clean up and return*/ + delete gauss; + xDelete(dbasis); + xDelete(xyz_list); + if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; + return Ke; +}/*}}}*/ +ElementMatrix* AdjointHorizAnalysis::CreateKMatrixHO(Element* element){/*{{{*/ + _error_("S"); +}/*}}}*/ +ElementMatrix* AdjointHorizAnalysis::CreateKMatrixFS(Element* element){/*{{{*/ + _error_("S"); }/*}}}*/ ElementVector* AdjointHorizAnalysis::CreatePVector(Element* element){/*{{{*/ Index: /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.h (revision 16875) +++ /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.h (revision 16876) @@ -22,4 +22,7 @@ /*Finite element Analysis*/ ElementMatrix* CreateKMatrix(Element* element); + ElementMatrix* CreateKMatrixSSA(Element* element); + ElementMatrix* CreateKMatrixHO(Element* element); + ElementMatrix* CreateKMatrixFS(Element* element); ElementVector* CreatePVector(Element* element); ElementVector* CreatePVectorSSA(Element* element); Index: /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp (revision 16875) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp (revision 16876) @@ -136,5 +136,95 @@ /*Finite Element Analysis*/ ElementMatrix* StressbalanceSIAAnalysis::CreateKMatrix(Element* element){/*{{{*/ - _error_("not implemented yet"); + int meshtype; + element->FindParam(&meshtype,MeshTypeEnum); + switch(meshtype){ + case Mesh2DhorizontalEnum: + return CreateKMatrix2D(element); + case Mesh3DEnum: + return CreateKMatrix3D(element); + default: _error_("mesh "<GetNumberOfVertices(); + + /*Initialize Element vector*/ + ElementMatrix* Ke=element->NewElementMatrix(); + for(int iv=0;ivVertexConnectivity(iv); + Ke->values[(2*iv+0)*2*numvertices+(2*iv+0)]=1./connectivity; + Ke->values[(2*iv+1)*2*numvertices+(2*iv+1)]=1./connectivity; + } + + /*Clean up and return*/ + return Ke; +}/*}}}*/ +ElementMatrix* StressbalanceSIAAnalysis::CreateKMatrix3D(Element* element){/*{{{*/ + + /*Intermediaries */ + int i0,i1,j0,j1,nodeup,nodedown,numsegments; + IssmDouble slope[2],connectivity[2],one0,one1; + int *pairindices = NULL; + + /*Fetch number vertices for this element*/ + int numvertices = element->GetNumberOfVertices(); + int numdof = 2*numvertices; + + /*Initialize Element vector*/ + ElementMatrix* Ke=element->NewElementMatrix(); + + element->VerticalSegmentIndices(&pairindices,&numsegments); + for(int is=0;isVertexConnectivity(nodedown); + connectivity[1]=(IssmDouble)element->VertexConnectivity(nodeup); + one0=1./connectivity[0]; + one1=1./connectivity[1]; + + /*2 dofs of first node*/ + i0=2*nodedown; i1=2*nodedown+1; + /*2 dofs of second node*/ + j0=2*nodeup; j1=2*nodeup+1; + + /*Create matrix for these two nodes*/ + if(element->IsOnBed() && element->IsOnSurface()){ + Ke->values[i0*numdof+i0] = +one0; + Ke->values[i1*numdof+i1] = +one0; + Ke->values[j0*numdof+i0] = -one1; + Ke->values[j0*numdof+j0] = +one1; + Ke->values[j1*numdof+i1] = -one1; + Ke->values[j1*numdof+j1] = +one1; + } + else if(element->IsOnBed()){ + Ke->values[i0*numdof+i0] = one0; + Ke->values[i1*numdof+i1] = one0; + Ke->values[j0*numdof+i0] = -2.*one1; + Ke->values[j0*numdof+j0] = +2.*one1; + Ke->values[j1*numdof+i1] = -2.*one1; + Ke->values[j1*numdof+j1] = +2.*one1; + } + else if(element->IsOnSurface()){ + Ke->values[j0*numdof+i0] = -one1; + Ke->values[j0*numdof+j0] = +one1; + Ke->values[j1*numdof+i1] = -one1; + Ke->values[j1*numdof+j1] = +one1; + } + else{ //node is on two horizontal layers and beams include the values only once, so the have to use half of the connectivity + Ke->values[j0*numdof+i0] = -2.*one1; + Ke->values[j0*numdof+j0] = +2.*one1; + Ke->values[j1*numdof+i1] = -2.*one1; + Ke->values[j1*numdof+j1] = +2.*one1; + } + } + + /*Clean up and return*/ + xDelete(pairindices); + return Ke; }/*}}}*/ ElementVector* StressbalanceSIAAnalysis::CreatePVector(Element* element){/*{{{*/ Index: /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.h (revision 16875) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.h (revision 16876) @@ -22,4 +22,6 @@ /*Finite element Analysis*/ ElementMatrix* CreateKMatrix(Element* element); + ElementMatrix* CreateKMatrix2D(Element* element); + ElementMatrix* CreateKMatrix3D(Element* element); ElementVector* CreatePVector(Element* element); ElementVector* CreatePVector2D(Element* element); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16875) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16876) @@ -56,4 +56,5 @@ virtual IssmDouble EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure)=0; virtual IssmDouble EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure)=0; + virtual void FindParam(bool* pvalue,int paramenum)=0; virtual void FindParam(int* pvalue,int paramenum)=0; virtual void FindParam(IssmDouble* pvalue,int paramenum)=0; @@ -170,4 +171,6 @@ virtual void ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=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 StrainRateSSA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_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 16875) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16876) @@ -867,4 +867,9 @@ IssmDouble Penta::EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){ return matpar->GetEnthalpyDiffusionParameterVolume(numvertices,enthalpy,pressure); +} +/*}}}*/ +/*FUNCTION Penta::FindParam(bool* pvalue,int paramenum){{{*/ +void Penta::FindParam(bool* pvalue,int paramenum){ + this->parameters->FindParam(pvalue,paramenum); } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16875) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16876) @@ -78,4 +78,5 @@ void DeleteMaterials(void){_error_("not implemented yet");}; void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz); + void FindParam(bool* pvalue,int paramenum); void FindParam(int* pvalue,int paramenum); void FindParam(IssmDouble* pvalue,int paramenum); @@ -297,4 +298,6 @@ void ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); void ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented");}; + void ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){_error_("not implemented");}; + 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.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.cpp (revision 16875) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.cpp (revision 16876) @@ -88,4 +88,9 @@ } /*}}}*/ +/*FUNCTION Seg::FindParam(bool* pvalue,int paramenum){{{*/ +void Seg::FindParam(bool* pvalue,int paramenum){ + this->parameters->FindParam(pvalue,paramenum); +} +/*}}}*/ /*FUNCTION Seg::FindParam(int* pvalue,int paramenum){{{*/ void Seg::FindParam(int* pvalue,int paramenum){ Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16875) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16876) @@ -91,4 +91,5 @@ IssmDouble EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented");}; IssmDouble EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){_error_("not implemented");}; + void FindParam(bool* pvalue,int paramenum); void FindParam(int* pvalue,int paramenum); void FindParam(IssmDouble* pvalue,int paramenum); @@ -172,4 +173,6 @@ void ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented");}; void ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented");}; + void ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){_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 16875) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16876) @@ -663,5 +663,5 @@ /*Compute strain rate viscosity and pressure: */ - this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); + this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); material->GetViscosity2d(&viscosity,&epsilon[0]); pressure_input->GetInputValue(&pressure,gauss); @@ -848,4 +848,9 @@ } /*}}}*/ +/*FUNCTION Tria::FindParam(bool* pvalue,int paramenum){{{*/ +void Tria::FindParam(bool* pvalue,int paramenum){ + this->parameters->FindParam(pvalue,paramenum); +} +/*}}}*/ /*FUNCTION Tria::FindParam(int* pvalue,int paramenum){{{*/ void Tria::FindParam(int* pvalue,int paramenum){ @@ -1562,26 +1567,4 @@ void Tria::GetConnectivityList(int* connectivity){ for(int i=0;iConnectivity(); -} -/*}}}*/ -/*FUNCTION Tria::GetStrainRate2d(IssmDouble* epsilon,IssmDouble* xyz_list, Gauss* gauss, Input* vx_input, Input* vy_input){{{*/ -void Tria::GetStrainRate2d(IssmDouble* epsilon,IssmDouble* xyz_list, Gauss* gauss, Input* vx_input, Input* vy_input){ - /*Compute the 2d Strain Rate (3 components): - * epsilon=[exx eyy exy] */ - - 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]; } /*}}}*/ @@ -3073,5 +3056,5 @@ IssmDouble epsilon[3];/* epsilon=[exx,eyy,exy]; */ - this->GetStrainRate2d(&epsilon[0],xyz_list,gauss,vx_input,vy_input); + this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input); material->GetViscosity2dvertical(&viscosity, &epsilon[0]); @@ -3087,9 +3070,34 @@ IssmDouble epsilon[3];/* epsilon=[exx,eyy,exy]; */ - this->GetStrainRate2d(&epsilon[0],xyz_list,gauss,vx_input,vy_input); + 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]; } /*}}}*/ @@ -3747,5 +3755,5 @@ GetBprimeFS(Bprime,&xyz_list[0][0],gauss); - this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); + this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); material->GetViscosity2dvertical(&viscosity,&epsilon[0]); @@ -3897,6 +3905,6 @@ GetBprimeSSA(&Bprime[0], &xyz_list[0][0], gauss); - this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); - this->GetStrainRate2d(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input); + this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); + this->StrainRateSSA(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input); material->GetViscosity2d(&viscosity, &epsilon[0]); material->GetViscosity2d(&oldviscosity, &oldepsilon[0]); @@ -4486,5 +4494,5 @@ thickness_input->GetInputValue(&thickness, gauss); - this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); + this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); material->GetViscosity2dDerivativeEpsSquare(&mu_prime,&epsilon[0]); eps1[0]=2*epsilon[0]+epsilon[1]; eps2[0]=epsilon[2]; @@ -4869,5 +4877,5 @@ adjointy_input->GetInputDerivativeValue(&dadjy[0],&xyz_list[0][0],gauss); - this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); + this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); material->GetViscosityComplement(&viscosity_complement,&epsilon[0]); @@ -4925,5 +4933,5 @@ adjointy_input->GetInputDerivativeValue(&dadjy[0],&xyz_list[0][0],gauss); - this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); + this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); material->GetViscosityDComplement(&viscosity_complement,&epsilon[0]); @@ -6203,5 +6211,5 @@ thickness_input->GetInputValue(&thickness, gauss); - this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); + this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); material->GetViscosity2dDerivativeEpsSquare(&mu_prime,&epsilon[0]); eps1[0]=2*epsilon[0]+epsilon[1]; eps2[0]=epsilon[2]; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16875) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16876) @@ -83,4 +83,5 @@ void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){_error_("not implemented yet");}; void EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented yet");}; + void FindParam(bool* pvalue,int paramenum); void FindParam(int* pvalue,int paramenum); void FindParam(IssmDouble* pvalue,int paramenum); @@ -285,5 +286,4 @@ void GetInputValue(IssmDouble* pvalue,Gauss* gauss,int enum_type); void GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype); - void GetStrainRate2d(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss, Input* vx_input, Input* vy_input); void InputChangeName(int enum_type,int enum_type_old); void InputUpdateFromSolutionOneDof(IssmDouble* solution,int enum_type); @@ -330,4 +330,6 @@ void ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented yet");}; void ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); + void ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon); + void StrainRateSSA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); #ifdef _HAVE_STRESSBALANCE_