Changeset 16876

Timestamp:

11/21/13 20:55:55 (11 years ago)

Author:

Mathieu Morlighem

Message:

NEW: Adjoint SSA and SIA

Location:

issm/trunk-jpl/src/c

Files:

• analyses/AdjointHorizAnalysis.cpp (modified) (1 diff)
• analyses/AdjointHorizAnalysis.h (modified) (1 diff)
• analyses/StressbalanceSIAAnalysis.cpp (modified) (1 diff)
• analyses/StressbalanceSIAAnalysis.h (modified) (1 diff)
• classes/Elements/Element.h (modified) (2 diffs)
• classes/Elements/Penta.cpp (modified) (1 diff)
• classes/Elements/Penta.h (modified) (2 diffs)
• classes/Elements/Seg.cpp (modified) (1 diff)
• classes/Elements/Seg.h (modified) (2 diffs)
• classes/Elements/Tria.cpp (modified) (11 diffs)
• classes/Elements/Tria.h (modified) (3 diffs)

issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp

@@ -27 +27 @@
 /*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 "<<EnumToStringx(approximation)<<" not supported");
+        }
+}/*}}}*/
+ElementMatrix* AdjointHorizAnalysis::CreateKMatrixSSA(Element* element){/*{{{*/
+
+        /*Intermediaries*/
+        int      meshtype;
+        Element* basalelement;
+
+        /*Get basal element*/
+        element->FindParam(&meshtype,MeshTypeEnum);
+        switch(meshtype){
+                case Mesh2DhorizontalEnum:
+                        basalelement = element;
+                        break;
+                case Mesh3DEnum:
+                        if(!element->IsOnBed()) return NULL;
+                        basalelement = element->SpawnBasalElement();
+                        break;
+                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+        }
+
+        /*Intermediaries */
+        bool        incomplete_adjoint;
+        IssmDouble  Jdet,thickness,mu_prime;
+        IssmDouble  eps1dotdphii,eps1dotdphij,eps2dotdphii,eps2dotdphij;
+        IssmDouble  eps1[2],eps2[2],epsilon[3];
+        IssmDouble *xyz_list = NULL;
+
+        /*Fetch number of nodes and dof for this finite element*/
+        int numnodes = basalelement->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<IssmDouble>(2*numnodes);
+
+        /* Start  looping on the number of gaussian points: */
+        Gauss* gauss=basalelement->NewGauss(2);
+        for(int ig=gauss->begin();ig<gauss->end();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;i<numnodes;i++){
+                        for(int j=0;j<numnodes;j++){
+                                eps1dotdphii=eps1[0]*dbasis[0*numnodes+i]+eps1[1]*dbasis[1*numnodes+i];
+                                eps1dotdphij=eps1[0]*dbasis[0*numnodes+j]+eps1[1]*dbasis[1*numnodes+j];
+                                eps2dotdphii=eps2[0]*dbasis[0*numnodes+i]+eps2[1]*dbasis[1*numnodes+i];
+                                eps2dotdphij=eps2[0]*dbasis[0*numnodes+j]+eps2[1]*dbasis[1*numnodes+j];
+
+                                Ke->values[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<IssmDouble>(dbasis);
+        xDelete<IssmDouble>(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){/*{{{*/

issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.h

@@ -22 +22 @@
                 /*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);

issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp

@@ -136 +136 @@
 /*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 "<<EnumToStringx(meshtype)<<" not supported yet");
+        }
+}/*}}}*/
+ElementMatrix* StressbalanceSIAAnalysis::CreateKMatrix2D(Element* element){/*{{{*/
+
+        /*Intermediaries */
+        IssmDouble connectivity;
+
+        /*Fetch number vertices for this element*/
+        int numvertices = element->GetNumberOfVertices();
+
+        /*Initialize Element vector*/
+        ElementMatrix* Ke=element->NewElementMatrix();
+        for(int iv=0;iv<numvertices;iv++){
+                connectivity=(IssmDouble)element->VertexConnectivity(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;is<numsegments;is++){
+                nodedown = pairindices[is*2+0];
+                nodeup   = pairindices[is*2+1];
+                connectivity[0]=(IssmDouble)element->VertexConnectivity(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<int>(pairindices);
+        return Ke;
 }/*}}}*/
 ElementVector* StressbalanceSIAAnalysis::CreatePVector(Element* element){/*{{{*/

issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.h

@@ -22 +22 @@
                 /*Finite element Analysis*/
                 ElementMatrix* CreateKMatrix(Element* element);
+                ElementMatrix* CreateKMatrix2D(Element* element);
+                ElementMatrix* CreateKMatrix3D(Element* element);
                 ElementVector* CreatePVector(Element* element);
                 ElementVector* CreatePVector2D(Element* element);

issm/trunk-jpl/src/c/classes/Elements/Element.h

@@ -56 +56 @@
                 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 +171 @@
                 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;

issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

@@ -867 +867 @@
 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);
 }
 /*}}}*/

issm/trunk-jpl/src/c/classes/Elements/Penta.h

@@ -78 +78 @@
                 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 +298 @@
                 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_

issm/trunk-jpl/src/c/classes/Elements/Seg.cpp

@@ -88 +88 @@
 }
 /*}}}*/
+/*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){

issm/trunk-jpl/src/c/classes/Elements/Seg.h

@@ -91 +91 @@
                 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 +173 @@
                 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");};

issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -663 +663 @@
 
                 /*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 +848 @@
 }
 /*}}}*/
+/*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 +1567 @@
 void  Tria::GetConnectivityList(int* connectivity){
         for(int i=0;i<NUMVERTICES;i++) connectivity[i]=vertices[i]->Connectivity();
-}
-/*}}}*/
-/*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 +3056 @@
         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 +3070 @@
         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 +3755 @@
                 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 +3905 @@
                 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 +4494 @@
 
                 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 +4877 @@
                 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 +4933 @@
                 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 +6211 @@
 
                 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];

issm/trunk-jpl/src/c/classes/Elements/Tria.h

@@ -83 +83 @@
                 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 +286 @@
                 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 +330 @@
                 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_