Changeset 18920

Timestamp:

12/03/14 20:36:06 (10 years ago)

Author:

seroussi

Message:

CHG: alphabetical order in Element.cpp

File:

• issm/trunk-jpl/src/c/classes/Elements/Element.cpp (modified) (20 diffs)

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

@@ -306 +306 @@
 }
 /*}}}*/
-void       Element::DeleteMaterials(void){/*{{{*/
-        delete this->material;
-}/*}}}*/
 void       Element::DeepEcho(void){/*{{{*/
 
@@ -343 +340 @@
 }
 /*}}}*/
-void       Element::Echo(void){/*{{{*/
-        _printf_(EnumToStringx(this->ObjectEnum())<<" element:\n");
-        _printf_("   id : "<<this->id <<"\n");
-        _printf_("   sid: "<<this->sid<<"\n");
-        if(vertices){
-                int numvertices = this->GetNumberOfVertices();
-                for(int i=0;i<numvertices;i++) vertices[i]->Echo();
-        }
-        else _printf_("vertices = NULL\n");
-
-        if(nodes){
-                int numnodes = this->GetNumberOfNodes();
-                for(int i=0;i<numnodes;i++) {
-                        _printf_("nodes[" << i << "] = " << nodes[i]);  
-                        nodes[i]->Echo();
-                }
-        }
-        else _printf_("nodes = NULL\n");
-
-        if (material) material->Echo();
-        else _printf_("material = NULL\n");
-
-        if (matpar) matpar->Echo();
-        else _printf_("matpar = NULL\n");
-
-        _printf_("   parameters\n");
-        if (parameters) parameters->Echo();
-        else _printf_("parameters = NULL\n");
-
-        _printf_("   inputs\n");
-        if (inputs) inputs->Echo();
-        else _printf_("inputs=NULL\n");
-}
-/*}}}*/
+void       Element::DeleteInput(int input_enum){/*{{{*/
+
+        inputs->DeleteInput(input_enum);
+
+}
+/*}}}*/
+void       Element::DeleteMaterials(void){/*{{{*/
+        delete this->material;
+}/*}}}*/
 IssmDouble Element::Divergence(void){/*{{{*/
         /*Compute element divergence*/
@@ -419 +391 @@
         return divergence;
 }/*}}}*/
+void       Element::dViscositydBFS(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/
+
+        /*Intermediaries*/
+        IssmDouble dmudB;
+        IssmDouble epsilon3d[6];/* epsilon=[exx,eyy,exy,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_B(&dmudB,eps_eff);
+
+        /*Assign output pointer*/
+        *pdmudB=dmudB;
+
+}
+/*}}}*/
+void       Element::dViscositydBHO(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
+
+        /*Intermediaries*/
+        IssmDouble dmudB;
+        IssmDouble epsilon3d[5];/* epsilon=[exx,eyy,exy,exy,exz,eyz];    */
+        IssmDouble epsilon2d[2];/* 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->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]+eps0*eps0);
+        }
+        else{
+                /* eps_eff^2 = 1/2 ( exx^2 + eyy^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] + eps0*eps0);
+        }
+        /*Get viscosity*/
+        material->GetViscosity_B(&dmudB,eps_eff);
+
+        /*Assign output pointer*/
+        *pdmudB=dmudB;
+
+}
+/*}}}*/
 void       Element::dViscositydBSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
 
@@ -446 +472 @@
 }
 /*}}}*/
-void       Element::dViscositydBHO(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
-
-        /*Intermediaries*/
-        IssmDouble dmudB;
-        IssmDouble epsilon3d[5];/* epsilon=[exx,eyy,exy,exy,exz,eyz];    */
-        IssmDouble epsilon2d[2];/* 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->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]+eps0*eps0);
-        }
-        else{
-                /* eps_eff^2 = 1/2 ( exx^2 + eyy^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] + eps0*eps0);
-        }
-        /*Get viscosity*/
-        material->GetViscosity_B(&dmudB,eps_eff);
-
-        /*Assign output pointer*/
-        *pdmudB=dmudB;
-
-}
-/*}}}*/
-void       Element::dViscositydBFS(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/
-
-        /*Intermediaries*/
-        IssmDouble dmudB;
-        IssmDouble epsilon3d[6];/* epsilon=[exx,eyy,exy,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_B(&dmudB,eps_eff);
-
-        /*Assign output pointer*/
-        *pdmudB=dmudB;
-
-}
-/*}}}*/
 void       Element::dViscositydDSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
 
@@ -527 +499 @@
 }
 /*}}}*/
-void       Element::ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){/*{{{*/
-        matpar->ThermalToEnthalpy(penthalpy,temperature,waterfraction,pressure);
+void       Element::Echo(void){/*{{{*/
+        _printf_(EnumToStringx(this->ObjectEnum())<<" element:\n");
+        _printf_("   id : "<<this->id <<"\n");
+        _printf_("   sid: "<<this->sid<<"\n");
+        if(vertices){
+                int numvertices = this->GetNumberOfVertices();
+                for(int i=0;i<numvertices;i++) vertices[i]->Echo();
+        }
+        else _printf_("vertices = NULL\n");
+
+        if(nodes){
+                int numnodes = this->GetNumberOfNodes();
+                for(int i=0;i<numnodes;i++) {
+                        _printf_("nodes[" << i << "] = " << nodes[i]);  
+                        nodes[i]->Echo();
+                }
+        }
+        else _printf_("nodes = NULL\n");
+
+        if (material) material->Echo();
+        else _printf_("material = NULL\n");
+
+        if (matpar) matpar->Echo();
+        else _printf_("matpar = NULL\n");
+
+        _printf_("   parameters\n");
+        if (parameters) parameters->Echo();
+        else _printf_("parameters = NULL\n");
+
+        _printf_("   inputs\n");
+        if (inputs) inputs->Echo();
+        else _printf_("inputs=NULL\n");
+}
+/*}}}*/
+IssmDouble Element::EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){/*{{{*/
+        return matpar->GetEnthalpyDiffusionParameter(enthalpy,pressure);
+}/*}}}*/
+IssmDouble Element::EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){/*{{{*/
+        return matpar->GetEnthalpyDiffusionParameterVolume(numvertices,enthalpy,pressure);
 }/*}}}*/
 void       Element::EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){/*{{{*/
         matpar->EnthalpyToThermal(ptemperature,pwaterfraction,enthalpy,pressure);
-}/*}}}*/
-IssmDouble Element::EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){/*{{{*/
-        return matpar->GetEnthalpyDiffusionParameter(enthalpy,pressure);
-}/*}}}*/
-IssmDouble Element::EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){/*{{{*/
-        return matpar->GetEnthalpyDiffusionParameterVolume(numvertices,enthalpy,pressure);
 }/*}}}*/
 void       Element::FindParam(bool* pvalue,int paramenum){/*{{{*/
@@ -574 +577 @@
 }
 /*}}}*/
+void       Element::GetDofListPressure(int** pdoflist,int setenum){/*{{{*/
+
+        /*Fetch number of nodes and dof for this finite element*/
+        int vnumnodes = this->NumberofNodesVelocity();
+        int pnumnodes = this->NumberofNodesPressure();
+
+        /*First, figure out size of doflist and create it: */
+        int numberofdofs=0;
+        for(int i=vnumnodes;i<vnumnodes+pnumnodes;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(FSApproximationEnum,setenum);
+
+        /*Allocate output*/
+        int* doflist=xNew<int>(numberofdofs);
+
+        /*Populate: */
+        int count=0;
+        for(int i=vnumnodes;i<vnumnodes+pnumnodes;i++){
+                nodes[i]->GetDofList(doflist+count,FSApproximationEnum,setenum);
+                count+=nodes[i]->GetNumberOfDofs(FSApproximationEnum,setenum);
+        }
+
+        /*Assign output pointers:*/
+        *pdoflist=doflist;
+}
+/*}}}*/
 void       Element::GetDofListVelocity(int** pdoflist,int setenum){/*{{{*/
 
@@ -597 +624 @@
 }
 /*}}}*/
-void       Element::GetDofListPressure(int** pdoflist,int setenum){/*{{{*/
-
-        /*Fetch number of nodes and dof for this finite element*/
-        int vnumnodes = this->NumberofNodesVelocity();
-        int pnumnodes = this->NumberofNodesPressure();
-
-        /*First, figure out size of doflist and create it: */
-        int numberofdofs=0;
-        for(int i=vnumnodes;i<vnumnodes+pnumnodes;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(FSApproximationEnum,setenum);
-
-        /*Allocate output*/
-        int* doflist=xNew<int>(numberofdofs);
-
-        /*Populate: */
-        int count=0;
-        for(int i=vnumnodes;i<vnumnodes+pnumnodes;i++){
-                nodes[i]->GetDofList(doflist+count,FSApproximationEnum,setenum);
-                count+=nodes[i]->GetNumberOfDofs(FSApproximationEnum,setenum);
-        }
-
-        /*Assign output pointers:*/
-        *pdoflist=doflist;
-}
-/*}}}*/
+Input*     Element::GetInput(int inputenum){/*{{{*/
+        return inputs->GetInput(inputenum);
+}/*}}}*/
+void       Element::GetInputListOnNodes(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue){/*{{{*/
+
+        _assert_(pvalue);
+
+        Input *input    = this->GetInput(enumtype);
+        int    numnodes = this->GetNumberOfNodes();
+
+        /* Start looping on the number of vertices: */
+        if(input){
+                Gauss* gauss=this->NewGauss();
+                for(int iv=0;iv<numnodes;iv++){
+                        gauss->GaussNode(this->FiniteElement(),iv);
+                        input->GetInputValue(&pvalue[iv],gauss);
+                }
+                delete gauss;
+        }
+        else{
+                for(int iv=0;iv<numnodes;iv++) pvalue[iv]=defaultvalue;
+        }
+}
+/*}}}*/
+void       Element::GetInputListOnNodes(IssmDouble* pvalue,int enumtype){/*{{{*/
+
+        _assert_(pvalue);
+
+        int    numnodes = this->GetNumberOfNodes();
+        Input *input    = this->GetInput(enumtype);
+        if(!input) _error_("Input " << EnumToStringx(enumtype) << " not found in element");
+
+        /* Start looping on the number of vertices: */
+        Gauss* gauss=this->NewGauss();
+        for(int iv=0;iv<numnodes;iv++){
+                gauss->GaussNode(this->FiniteElement(),iv);
+                input->GetInputValue(&pvalue[iv],gauss);
+        }
+        delete gauss;
+}
+/*}}}*/
+void       Element::GetInputListOnNodesVelocity(IssmDouble* pvalue,int enumtype){/*{{{*/
+
+        _assert_(pvalue);
+
+        int    numnodes = this->NumberofNodesVelocity();
+        Input *input    = this->GetInput(enumtype);
+        if(!input) _error_("Input " << EnumToStringx(enumtype) << " not found in element");
+
+        /* Start looping on the number of vertices: */
+        Gauss* gauss=this->NewGauss();
+        for(int iv=0;iv<numnodes;iv++){
+                gauss->GaussNode(this->VelocityInterpolation(),iv);
+                input->GetInputValue(&pvalue[iv],gauss);
+        }
+        delete gauss;
+}
+/*}}}*/
+void       Element::GetInputListOnVertices(IssmDouble* pvalue,int enumtype){/*{{{*/
+
+        /*Recover input*/
+        Input* input=this->GetInput(enumtype);
+        if (!input) _error_("Input " << EnumToStringx(enumtype) << " not found in element");
+
+        /*Fetch number vertices for this element*/
+        int numvertices = this->GetNumberOfVertices();
+
+        /*Checks in debugging mode*/
+        _assert_(pvalue);
+
+        /* Start looping on the number of vertices: */
+        Gauss*gauss=this->NewGauss();
+        for(int iv=0;iv<numvertices;iv++){
+                gauss->GaussVertex(iv);
+                input->GetInputValue(&pvalue[iv],gauss);
+        }
+
+        /*clean-up*/
+        delete gauss;
+}
+/*}}}*/
+void       Element::GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue){/*{{{*/
+
+        /*Recover input*/
+        Input* input=this->GetInput(enumtype);
+
+        /*Checks in debugging mode*/
+        _assert_(pvalue);
+
+        /*Fetch number vertices for this element*/
+        int numvertices = this->GetNumberOfVertices();
+
+        /* Start looping on the number of vertices: */
+        if (input){
+                Gauss* gauss=this->NewGauss();
+                for (int iv=0;iv<numvertices;iv++){
+                        gauss->GaussVertex(iv);
+                        input->GetInputValue(&pvalue[iv],gauss);
+                }
+                delete gauss;
+        }
+        else{
+                for(int iv=0;iv<numvertices;iv++) pvalue[iv]=defaultvalue;
+        }
+}
+/*}}}*/
+void       Element::GetInputLocalMinMaxOnNodes(IssmDouble* min,IssmDouble* max,IssmDouble* ug){/*{{{*/
+
+
+        /*Get number of nodes for this element*/
+        int numnodes = this->GetNumberOfNodes();
+
+        /*Some checks to avoid segmentation faults*/
+        _assert_(ug);
+        _assert_(numnodes>0);
+        _assert_(nodes);
+
+        /*Get element minimum/maximum*/
+        IssmDouble input_min = ug[nodes[0]->GetDof(0,GsetEnum)];
+        IssmDouble input_max = input_min;
+        for(int i=1;i<numnodes;i++){
+                if(ug[nodes[i]->GetDof(0,GsetEnum)] < input_min) input_min = ug[nodes[i]->GetDof(0,GsetEnum)];
+                if(ug[nodes[i]->GetDof(0,GsetEnum)] > input_max) input_max = ug[nodes[i]->GetDof(0,GsetEnum)];
+        }
+
+
+        /*Second loop to reassign min and max with local extrema*/
+        for(int i=0;i<numnodes;i++){
+                if(min[nodes[i]->GetDof(0,GsetEnum)]>input_min) min[nodes[i]->GetDof(0,GsetEnum)] = input_min;
+                if(max[nodes[i]->GetDof(0,GsetEnum)]<input_max) max[nodes[i]->GetDof(0,GsetEnum)] = input_max;
+        }
+}
+/*}}}*/
+void       Element::GetInputValue(bool* pvalue,int inputenum){/*{{{*/
+
+        Input* input=inputs->GetInput(inputenum);
+        if(!input) _error_("Input " << EnumToStringx(inputenum) << " not found in element");
+        input->GetInputValue(pvalue);
+
+}/*}}}*/
+void       Element::GetInputValue(int* pvalue,int inputenum){/*{{{*/
+
+        Input* input=inputs->GetInput(inputenum);
+        if(!input) _error_("Input " << EnumToStringx(inputenum) << " not found in element");
+        input->GetInputValue(pvalue);
+
+}/*}}}*/
+void       Element::GetInputValue(IssmDouble* pvalue,int inputenum){/*{{{*/
+
+        Input* input=inputs->GetInput(inputenum);
+        if(!input) _error_("Input " << EnumToStringx(inputenum) << " not found in element");
+        input->GetInputValue(pvalue);
+
+}/*}}}*/
+void       Element::GetInputValue(IssmDouble* pvalue,Gauss* gauss,int inputenum){/*{{{*/
+
+        Input* input=inputs->GetInput(inputenum);
+        if(!input) _error_("Input " << EnumToStringx(inputenum) << " not found in element");
+        input->GetInputValue(pvalue,gauss);
+
+}/*}}}*/
 IssmDouble Element::GetMaterialParameter(int enum_in){/*{{{*/
 
@@ -635 +799 @@
         }
 }/*}}}*/
+void       Element::GetNodesLidList(int* lidlist){/*{{{*/
+
+        _assert_(lidlist);
+        _assert_(nodes);
+        int numnodes = this->GetNumberOfNodes();
+        for(int i=0;i<numnodes;i++){
+                lidlist[i]=nodes[i]->Lid();
+        }
+}
+/*}}}*/
+void       Element::GetNodesSidList(int* sidlist){/*{{{*/
+
+        _assert_(sidlist);
+        _assert_(nodes);
+        int numnodes = this->GetNumberOfNodes();
+        for(int i=0;i<numnodes;i++){
+                sidlist[i]=nodes[i]->Sid();
+        }
+}
+/*}}}*/
 void       Element::GetPhi(IssmDouble* phi, IssmDouble*  epsilon, IssmDouble viscosity){/*{{{*/
         /*Compute deformational heating from epsilon and viscosity */
@@ -672 +856 @@
 }
 /*}}}*/
-Input*     Element::GetInput(int inputenum){/*{{{*/
-        return inputs->GetInput(inputenum);
-}/*}}}*/
-void       Element::GetInputListOnVertices(IssmDouble* pvalue,int enumtype){/*{{{*/
-
-        /*Recover input*/
-        Input* input=this->GetInput(enumtype);
-        if (!input) _error_("Input " << EnumToStringx(enumtype) << " not found in element");
-
-        /*Fetch number vertices for this element*/
-        int numvertices = this->GetNumberOfVertices();
-
-        /*Checks in debugging mode*/
-        _assert_(pvalue);
-
-        /* Start looping on the number of vertices: */
-        Gauss*gauss=this->NewGauss();
-        for(int iv=0;iv<numvertices;iv++){
-                gauss->GaussVertex(iv);
-                input->GetInputValue(&pvalue[iv],gauss);
-        }
-
-        /*clean-up*/
-        delete gauss;
-}
-/*}}}*/
-void       Element::GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue){/*{{{*/
-
-        /*Recover input*/
-        Input* input=this->GetInput(enumtype);
-
-        /*Checks in debugging mode*/
-        _assert_(pvalue);
-
-        /*Fetch number vertices for this element*/
-        int numvertices = this->GetNumberOfVertices();
-
-        /* Start looping on the number of vertices: */
-        if (input){
-                Gauss* gauss=this->NewGauss();
-                for (int iv=0;iv<numvertices;iv++){
-                        gauss->GaussVertex(iv);
-                        input->GetInputValue(&pvalue[iv],gauss);
-                }
-                delete gauss;
-        }
-        else{
-                for(int iv=0;iv<numvertices;iv++) pvalue[iv]=defaultvalue;
-        }
-}
-/*}}}*/
-void       Element::GetInputListOnNodes(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue){/*{{{*/
-
-        _assert_(pvalue);
-
-        Input *input    = this->GetInput(enumtype);
-        int    numnodes = this->GetNumberOfNodes();
-
-        /* Start looping on the number of vertices: */
-        if(input){
-                Gauss* gauss=this->NewGauss();
-                for(int iv=0;iv<numnodes;iv++){
-                        gauss->GaussNode(this->FiniteElement(),iv);
-                        input->GetInputValue(&pvalue[iv],gauss);
-                }
-                delete gauss;
-        }
-        else{
-                for(int iv=0;iv<numnodes;iv++) pvalue[iv]=defaultvalue;
-        }
-}
-/*}}}*/
-void       Element::GetInputListOnNodes(IssmDouble* pvalue,int enumtype){/*{{{*/
-
-        _assert_(pvalue);
-
-        int    numnodes = this->GetNumberOfNodes();
-        Input *input    = this->GetInput(enumtype);
-        if(!input) _error_("Input " << EnumToStringx(enumtype) << " not found in element");
-
-        /* Start looping on the number of vertices: */
-        Gauss* gauss=this->NewGauss();
-        for(int iv=0;iv<numnodes;iv++){
-                gauss->GaussNode(this->FiniteElement(),iv);
-                input->GetInputValue(&pvalue[iv],gauss);
-        }
-        delete gauss;
-}
-/*}}}*/
-void       Element::GetInputListOnNodesVelocity(IssmDouble* pvalue,int enumtype){/*{{{*/
-
-        _assert_(pvalue);
-
-        int    numnodes = this->NumberofNodesVelocity();
-        Input *input    = this->GetInput(enumtype);
-        if(!input) _error_("Input " << EnumToStringx(enumtype) << " not found in element");
-
-        /* Start looping on the number of vertices: */
-        Gauss* gauss=this->NewGauss();
-        for(int iv=0;iv<numnodes;iv++){
-                gauss->GaussNode(this->VelocityInterpolation(),iv);
-                input->GetInputValue(&pvalue[iv],gauss);
-        }
-        delete gauss;
-}
-/*}}}*/
-void       Element::GetInputLocalMinMaxOnNodes(IssmDouble* min,IssmDouble* max,IssmDouble* ug){/*{{{*/
-
-
-        /*Get number of nodes for this element*/
-        int numnodes = this->GetNumberOfNodes();
-
-        /*Some checks to avoid segmentation faults*/
-        _assert_(ug);
-        _assert_(numnodes>0);
-        _assert_(nodes);
-
-        /*Get element minimum/maximum*/
-        IssmDouble input_min = ug[nodes[0]->GetDof(0,GsetEnum)];
-        IssmDouble input_max = input_min;
-        for(int i=1;i<numnodes;i++){
-                if(ug[nodes[i]->GetDof(0,GsetEnum)] < input_min) input_min = ug[nodes[i]->GetDof(0,GsetEnum)];
-                if(ug[nodes[i]->GetDof(0,GsetEnum)] > input_max) input_max = ug[nodes[i]->GetDof(0,GsetEnum)];
-        }
-
-
-        /*Second loop to reassign min and max with local extrema*/
-        for(int i=0;i<numnodes;i++){
-                if(min[nodes[i]->GetDof(0,GsetEnum)]>input_min) min[nodes[i]->GetDof(0,GsetEnum)] = input_min;
-                if(max[nodes[i]->GetDof(0,GsetEnum)]<input_max) max[nodes[i]->GetDof(0,GsetEnum)] = input_max;
-        }
-}
-/*}}}*/
-void       Element::GetInputValue(bool* pvalue,int inputenum){/*{{{*/
-
-        Input* input=inputs->GetInput(inputenum);
-        if(!input) _error_("Input " << EnumToStringx(inputenum) << " not found in element");
-        input->GetInputValue(pvalue);
-
-}/*}}}*/
-void       Element::GetInputValue(int* pvalue,int inputenum){/*{{{*/
-
-        Input* input=inputs->GetInput(inputenum);
-        if(!input) _error_("Input " << EnumToStringx(inputenum) << " not found in element");
-        input->GetInputValue(pvalue);
-
-}/*}}}*/
-void       Element::GetInputValue(IssmDouble* pvalue,int inputenum){/*{{{*/
-
-        Input* input=inputs->GetInput(inputenum);
-        if(!input) _error_("Input " << EnumToStringx(inputenum) << " not found in element");
-        input->GetInputValue(pvalue);
-
-}/*}}}*/
-void       Element::GetInputValue(IssmDouble* pvalue,Gauss* gauss,int inputenum){/*{{{*/
-
-        Input* input=inputs->GetInput(inputenum);
-        if(!input) _error_("Input " << EnumToStringx(inputenum) << " not found in element");
-        input->GetInputValue(pvalue,gauss);
-
-}/*}}}*/
-void       Element::GetNodesSidList(int* sidlist){/*{{{*/
-
-        _assert_(sidlist);
-        _assert_(nodes);
-        int numnodes = this->GetNumberOfNodes();
-        for(int i=0;i<numnodes;i++){
-                sidlist[i]=nodes[i]->Sid();
-        }
-}
-/*}}}*/
-void       Element::GetNodesLidList(int* lidlist){/*{{{*/
-
-        _assert_(lidlist);
-        _assert_(nodes);
-        int numnodes = this->GetNumberOfNodes();
-        for(int i=0;i<numnodes;i++){
-                lidlist[i]=nodes[i]->Lid();
-        }
-}
-/*}}}*/
 void       Element::GetVectorFromInputs(Vector<IssmDouble>* vector,int input_enum){/*{{{*/
 
@@ -878 +881 @@
 }
 /*}}}*/
+void       Element::GetVerticesConnectivityList(int* connectivity){/*{{{*/
+
+        int numvertices = this->GetNumberOfVertices();
+        for(int i=0;i<numvertices;i++) connectivity[i]=this->vertices[i]->Connectivity();
+}
+/*}}}*/
 void       Element::GetVerticesCoordinates(IssmDouble** pxyz_list){/*{{{*/
 
@@ -891 +900 @@
         int numvertices = this->GetNumberOfVertices();
         for(int i=0;i<numvertices;i++) sidlist[i]=this->vertices[i]->Sid();
-}
-/*}}}*/
-void       Element::GetVerticesConnectivityList(int* connectivity){/*{{{*/
-
-        int numvertices = this->GetNumberOfVertices();
-        for(int i=0;i<numvertices;i++) connectivity[i]=this->vertices[i]->Connectivity();
 }
 /*}}}*/
@@ -1055 +1058 @@
         /*Call inputs method*/
         if (IsInput(original_enum)) inputs->DuplicateInput(original_enum,new_enum);
-
-}
-/*}}}*/
-void       Element::DeleteInput(int input_enum){/*{{{*/
-
-        inputs->DeleteInput(input_enum);
 
 }
@@ -1298 +1295 @@
 }
 /*}}}*/
+ElementMatrix* Element::NewElementMatrix(int approximation_enum){/*{{{*/
+        return new ElementMatrix(nodes,this->GetNumberOfNodes(),this->parameters,approximation_enum);
+}
+/*}}}*/
+ElementMatrix* Element::NewElementMatrixCoupling(int number_nodes,int approximation_enum){/*{{{*/
+        return new ElementMatrix(nodes,number_nodes,this->parameters,approximation_enum);
+}
+/*}}}*/
 ElementVector* Element::NewElementVector(int approximation_enum){/*{{{*/
         return new ElementVector(nodes,this->GetNumberOfNodes(),this->parameters,approximation_enum);
-}
-/*}}}*/
-ElementMatrix* Element::NewElementMatrix(int approximation_enum){/*{{{*/
-        return new ElementMatrix(nodes,this->GetNumberOfNodes(),this->parameters,approximation_enum);
-}
-/*}}}*/
-ElementMatrix* Element::NewElementMatrixCoupling(int number_nodes,int approximation_enum){/*{{{*/
-        return new ElementMatrix(nodes,number_nodes,this->parameters,approximation_enum);
 }
 /*}}}*/
@@ -1393 +1390 @@
         *pnodesperelement = input->GetResultNumberOfNodes();
 }/*}}}*/
+void       Element::ResultToPatch(IssmDouble* values,int nodesperelement,int output_enum){/*{{{*/
+
+        Input* input=this->inputs->GetInput(output_enum);
+        if(!input) _error_("input "<<EnumToStringx(output_enum)<<" not found in element");
+
+        input->ResultToPatch(values,nodesperelement,this->Sid());
+
+} /*}}}*/
 void       Element::ResultToVector(Vector<IssmDouble>* vector,int output_enum){/*{{{*/
 
@@ -1440 +1445 @@
                                          _error_("interpolation "<<EnumToStringx(input->GetResultInterpolation())<<" not supported yet");
         }
-} /*}}}*/
-void       Element::ResultToPatch(IssmDouble* values,int nodesperelement,int output_enum){/*{{{*/
-
-        Input* input=this->inputs->GetInput(output_enum);
-        if(!input) _error_("input "<<EnumToStringx(output_enum)<<" not found in element");
-
-        input->ResultToPatch(values,nodesperelement,this->Sid());
-
 } /*}}}*/
 void       Element::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){/*{{{*/
@@ -1525 +1522 @@
 }
 /*}}}*/
-IssmDouble Element::TMeltingPoint(IssmDouble pressure){/*{{{*/
-        _assert_(matpar);
-        return this->matpar->TMeltingPoint(pressure);
+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]
+         */
+
+        /*Intermediaries*/
+        IssmDouble dvx[3];
+        IssmDouble dvy[3];
+        IssmDouble dvz[3];
+
+        /*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->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
+        vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
+        vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss);
+        epsilon[0] = dvx[0];
+        epsilon[1] = dvy[1];
+        epsilon[2] = dvz[2];
+        epsilon[3] = 0.5*(dvx[1] + dvy[0]);
+        epsilon[4] = 0.5*(dvx[2] + dvz[0]);
+        epsilon[5] = 0.5*(dvy[2] + dvz[1]);
+
+}/*}}}*/
+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
+         */
+
+        /*Intermediaries*/
+        IssmDouble dvx[3];
+        IssmDouble dvy[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->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
+        vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
+        epsilon[0] = dvx[0];
+        epsilon[1] = dvy[1];
+        epsilon[2] = 0.5*(dvx[1] + dvy[0]);
+        epsilon[3] = 0.5*dvx[2];
+        epsilon[4] = 0.5*dvy[2];
+
+}/*}}}*/
+void       Element::StrainRateHO2dvertical(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
+        /*Compute the 2d Blatter/HOStrain Rate (2 components):
+         *
+         * epsilon=[exx exz]
+         *
+         * with exz=1/2 du/dz
+         *
+         * the contribution of vz is neglected
+         */
+
+        /*Intermediaries*/
+        IssmDouble dvx[3];
+
+        /*Check that both inputs have been found*/
+        if (!vx_input){
+                _error_("Input missing. Here are the input pointers we have for vx: " << vx_input <<"\n");
+        }
+
+        /*Get strain rate assuming that epsilon has been allocated*/
+        vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
+        epsilon[0] = dvx[0];
+        epsilon[1] = 0.5*dvx[1];
+
+}/*}}}*/
+void       Element::StrainRateSSA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
+
+        /*Intermediaries*/
+        IssmDouble dvx[3];
+        IssmDouble dvy[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->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
+        vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
+        epsilon[0] = dvx[0];
+        epsilon[1] = dvy[1];
+        epsilon[2] = 0.5*(dvx[1] + dvy[0]);
+
+}/*}}}*/
+void       Element::StrainRateSSA1d(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input){/*{{{*/
+
+        /*Intermediaries*/
+        IssmDouble dvx[3];
+
+        /*Check that both inputs have been found*/
+        if (!vx_input){
+                _error_("Input missing. Here are the input pointers we have for vx: " << vx_input << "\n");
+        }
+
+        /*Get strain rate assuming that epsilon has been allocated*/
+        vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
+        *epsilon = dvx[0];
+
 }/*}}}*/
 void       Element::StressMaxPrincipalCreateInput(void){/*{{{*/
@@ -1615 +1724 @@
 }
 /*}}}*/
-void       Element::ViscousHeatingCreateInput(void){/*{{{*/
-
-        /*Intermediaries*/
-        IssmDouble phi;
-        IssmDouble viscosity;
-        IssmDouble epsilon[6];
-        IssmDouble thickness;
-        IssmDouble *xyz_list = NULL;
-
-        /*Fetch number vertices and allocate memory*/
-        int         numvertices    = this->GetNumberOfVertices();
-        IssmDouble* viscousheating = xNew<IssmDouble>(numvertices);
-
-        /*Retrieve all inputs and parameters*/
-        this->GetVerticesCoordinatesBase(&xyz_list);
-        Input* vx_input        = this->GetInput(VxEnum); _assert_(vx_input);
-        Input* vy_input        = this->GetInput(VyEnum); _assert_(vy_input);
-        Input* vz_input        = this->GetInput(VzEnum); _assert_(vz_input);
-        Input* thickness_input = this->GetInput(ThicknessEnum); _assert_(thickness_input);
-
-        /*loop over vertices: */
-        Gauss* gauss=this->NewGauss();
-        for (int iv=0;iv<numvertices;iv++){
-                gauss->GaussVertex(iv);
-
-                thickness_input->GetInputValue(&thickness,gauss);
-
-                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->GetPhi(&phi,&epsilon[0],viscosity);
-
-                viscousheating[iv]=phi*thickness;
-        }
-
-        /*Create PentaVertex input, which will hold the basal friction:*/
-        this->AddInput(ViscousHeatingEnum,viscousheating,P1Enum);
-
-        /*Clean up and return*/
-        xDelete<IssmDouble>(viscousheating);
-        delete gauss;
-}
-/*}}}*/
+void       Element::ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){/*{{{*/
+        matpar->ThermalToEnthalpy(penthalpy,temperature,waterfraction,pressure);
+}/*}}}*/
+IssmDouble Element::TMeltingPoint(IssmDouble pressure){/*{{{*/
+        _assert_(matpar);
+        return this->matpar->TMeltingPoint(pressure);
+}/*}}}*/
+void       Element::TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,int transformenum){/*{{{*/
+
+        /*All nodes have the same Coordinate System*/
+        int numnodes  = this->GetNumberOfNodes();
+        int* cs_array = xNew<int>(numnodes);
+        for(int i=0;i<numnodes;i++) cs_array[i]=transformenum;
+
+        /*Call core*/
+        TransformInvStiffnessMatrixCoord(Ke,this->nodes,numnodes,cs_array);
+
+        /*Clean-up*/
+        xDelete<int>(cs_array);
+}/*}}}*/
+void       Element::TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes_list,int numnodes,int* cs_array){/*{{{*/
+
+        int         i,j;
+        int         numdofs   = 0;
+        IssmDouble *transform = NULL;
+        IssmDouble *values    = NULL;
+
+        /*Get total number of dofs*/
+        for(i=0;i<numnodes;i++){
+                switch(cs_array[i]){
+                        case PressureEnum: numdofs+=1; break;
+                        case XYEnum:       numdofs+=2; break;
+                        case XYZEnum:      numdofs+=3; break;
+                        default: _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
+                }
+        }
+
+        /*Copy current stiffness matrix*/
+        values=xNew<IssmDouble>(Ke->nrows*Ke->ncols);
+        for(i=0;i<Ke->nrows;i++) for(j=0;j<Ke->ncols;j++) values[i*Ke->ncols+j]=Ke->values[i*Ke->ncols+j];
+
+        /*Get Coordinate Systems transform matrix*/
+        CoordinateSystemTransform(&transform,nodes_list,numnodes,cs_array);
+
+        /*Transform matrix: R*Ke*R^T */
+        TripleMultiply(transform,numdofs,numdofs,0,
+                                values,Ke->nrows,Ke->ncols,0,
+                                transform,numdofs,numdofs,1,
+                                &Ke->values[0],0);
+
+        /*Free Matrix*/
+        xDelete<IssmDouble>(transform);
+        xDelete<IssmDouble>(values);
+}/*}}}*/
+void       Element::TransformLoadVectorCoord(ElementVector* pe,int transformenum){/*{{{*/
+
+        /*All nodes have the same Coordinate System*/
+        int  numnodes = this->GetNumberOfNodes();
+        int* cs_array = xNew<int>(numnodes);
+        for(int i=0;i<numnodes;i++) cs_array[i]=transformenum;
+
+        /*Call core*/
+        this->TransformLoadVectorCoord(pe,this->nodes,numnodes,cs_array);
+
+        /*Clean-up*/
+        xDelete<int>(cs_array);
+}/*}}}*/
+void       Element::TransformLoadVectorCoord(ElementVector* pe,int* cs_array){/*{{{*/
+
+        this->TransformLoadVectorCoord(pe,this->nodes,this->GetNumberOfNodes(),cs_array);
+
+}/*}}}*/
+void       Element::TransformLoadVectorCoord(ElementVector* pe,Node** nodes_list,int numnodes,int* cs_array){/*{{{*/
+
+        int         i;
+        int         numdofs   = 0;
+        IssmDouble *transform = NULL;
+        IssmDouble *values    = NULL;
+
+        /*Get total number of dofs*/
+        for(i=0;i<numnodes;i++){
+                switch(cs_array[i]){
+                        case PressureEnum: numdofs+=1; break;
+                        case XYEnum:       numdofs+=2; break;
+                        case XYZEnum:      numdofs+=3; break;
+                        default: _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
+                }
+        }
+
+        /*Copy current load vector*/
+        values=xNew<IssmDouble>(pe->nrows);
+        for(i=0;i<pe->nrows;i++) values[i]=pe->values[i];
+
+        /*Get Coordinate Systems transform matrix*/
+        CoordinateSystemTransform(&transform,nodes_list,numnodes,cs_array);
+
+        /*Transform matrix: R^T*pe */
+        MatrixMultiply(transform,numdofs,numdofs,1,
+                                values,pe->nrows,1,0,
+                                &pe->values[0],0);
+
+        /*Free Matrices*/
+        xDelete<IssmDouble>(transform);
+        xDelete<IssmDouble>(values);
+}/*}}}*/
+void       Element::TransformSolutionCoord(IssmDouble* values,int transformenum){/*{{{*/
+
+        /*All nodes have the same Coordinate System*/
+        int  numnodes = this->GetNumberOfNodes();
+        int* cs_array = xNew<int>(numnodes);
+        for(int i=0;i<numnodes;i++) cs_array[i]=transformenum;
+
+        /*Call core*/
+        this->TransformSolutionCoord(values,this->nodes,numnodes,cs_array);
+
+        /*Clean-up*/
+        xDelete<int>(cs_array);
+}/*}}}*/
+void       Element::TransformSolutionCoord(IssmDouble* values,int* transformenum_list){/*{{{*/
+        this->TransformSolutionCoord(values,this->nodes,this->GetNumberOfNodes(),transformenum_list);
+}/*}}}*/
+void       Element::TransformSolutionCoord(IssmDouble* values,int numnodes,int transformenum){/*{{{*/
+
+        /*All nodes have the same Coordinate System*/
+        int* cs_array = xNew<int>(numnodes);
+        for(int i=0;i<numnodes;i++) cs_array[i]=transformenum;
+
+        /*Call core*/
+        this->TransformSolutionCoord(values,this->nodes,numnodes,cs_array);
+
+        /*Clean-up*/
+        xDelete<int>(cs_array);
+}/*}}}*/
+void       Element::TransformSolutionCoord(IssmDouble* solution,int numnodes,int* cs_array){/*{{{*/
+        this->TransformSolutionCoord(solution,this->nodes,numnodes,cs_array);
+}/*}}}*/
+void       Element::TransformSolutionCoord(IssmDouble* values,Node** nodes_list,int numnodes,int transformenum){/*{{{*/
+        /*NOT NEEDED*/
+        /*All nodes have the same Coordinate System*/
+        int* cs_array = xNew<int>(numnodes);
+        for(int i=0;i<numnodes;i++) cs_array[i]=transformenum;
+
+        /*Call core*/
+        this->TransformSolutionCoord(values,nodes_list,numnodes,cs_array);
+
+        /*Clean-up*/
+        xDelete<int>(cs_array);
+}/*}}}*/
+void       Element::TransformSolutionCoord(IssmDouble* solution,Node** nodes_list,int numnodes,int* cs_array){/*{{{*/
+
+        int         i;
+        int         numdofs   = 0;
+        IssmDouble *transform = NULL;
+        IssmDouble *values    = NULL;
+
+        /*Get total number of dofs*/
+        for(i=0;i<numnodes;i++){
+                switch(cs_array[i]){
+                        case PressureEnum: numdofs+=1; break;
+                        case XYEnum:       numdofs+=2; break;
+                        case XYZEnum:      numdofs+=3; break;
+                        default: _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
+                }
+        }
+
+        /*Copy current solution vector*/
+        values=xNew<IssmDouble>(numdofs);
+        for(i=0;i<numdofs;i++) values[i]=solution[i];
+
+        /*Get Coordinate Systems transform matrix*/
+        CoordinateSystemTransform(&transform,nodes_list,numnodes,cs_array);
+
+        /*Transform matrix: R*U */
+        MatrixMultiply(transform,numdofs,numdofs,0,
+                                values,numdofs,1,0,
+                                &solution[0],0);
+
+        /*Free Matrices*/
+        xDelete<IssmDouble>(transform);
+        xDelete<IssmDouble>(values);
+}/*}}}*/
+void       Element::TransformStiffnessMatrixCoord(ElementMatrix* Ke,int transformenum){/*{{{*/
+
+        /*All nodes have the same Coordinate System*/
+        int  numnodes = this->GetNumberOfNodes();
+        int* cs_array = xNew<int>(numnodes);
+        for(int i=0;i<numnodes;i++) cs_array[i]=transformenum;
+
+        /*Call core*/
+        this->TransformStiffnessMatrixCoord(Ke,this->nodes,numnodes,cs_array);
+
+        /*Clean-up*/
+        xDelete<int>(cs_array);
+}/*}}}*/
+void       Element::TransformStiffnessMatrixCoord(ElementMatrix* Ke,int* transformenum_list){/*{{{*/
+        this->TransformStiffnessMatrixCoord(Ke,this->nodes,this->GetNumberOfNodes(),transformenum_list);
+}/*}}}*/
+void       Element::TransformStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes_list,int numnodes,int* cs_array){/*{{{*/
+
+        int         numdofs = 0;
+        IssmDouble *transform = NULL;
+        IssmDouble *values    = NULL;
+
+        /*Get total number of dofs*/
+        for(int i=0;i<numnodes;i++){
+                switch(cs_array[i]){
+                        case PressureEnum: numdofs+=1; break;
+                        case XYEnum:       numdofs+=2; break;
+                        case XYZEnum:      numdofs+=3; break;
+                        default: _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
+                }
+        }
+
+        /*Copy current stiffness matrix*/
+        values=xNew<IssmDouble>(Ke->nrows*Ke->ncols);
+        for(int i=0;i<Ke->nrows*Ke->ncols;i++) values[i]=Ke->values[i];
+
+        /*Get Coordinate Systems transform matrix*/
+        CoordinateSystemTransform(&transform,nodes_list,numnodes,cs_array);
+
+        /*Transform matrix: R^T*Ke*R */
+        TripleMultiply(transform,numdofs,numdofs,1,
+                                values,Ke->nrows,Ke->ncols,0,
+                                transform,numdofs,numdofs,0,
+                                &Ke->values[0],0);
+
+        /*Free Matrix*/
+        xDelete<IssmDouble>(transform);
+        xDelete<IssmDouble>(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,
@@ -1696 +1993 @@
 }
 /*}}}*/
+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
@@ -1755 +2083 @@
         *pviscosity = viscosity;
 }/*}}}*/
-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::ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
 
@@ -1808 +2111 @@
         this->material->GetViscosity2dDerivativeEpsSquare(pmu_prime,epsilon);
 }/*}}}*/
-void       Element::ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
-        this->material->GetViscosityDerivativeEpsSquare(pmu_prime,epsilon);
-}/*}}}*/
-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]
-         */
+void       Element::ViscousHeatingCreateInput(void){/*{{{*/
 
         /*Intermediaries*/
-        IssmDouble dvx[3];
-        IssmDouble dvy[3];
-        IssmDouble dvz[3];
-
-        /*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->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
-        vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
-        vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss);
-        epsilon[0] = dvx[0];
-        epsilon[1] = dvy[1];
-        epsilon[2] = dvz[2];
-        epsilon[3] = 0.5*(dvx[1] + dvy[0]);
-        epsilon[4] = 0.5*(dvx[2] + dvz[0]);
-        epsilon[5] = 0.5*(dvy[2] + dvz[1]);
-
-}/*}}}*/
-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
-         */
-
-        /*Intermediaries*/
-        IssmDouble dvx[3];
-        IssmDouble dvy[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->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
-        vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
-        epsilon[0] = dvx[0];
-        epsilon[1] = dvy[1];
-        epsilon[2] = 0.5*(dvx[1] + dvy[0]);
-        epsilon[3] = 0.5*dvx[2];
-        epsilon[4] = 0.5*dvy[2];
-
-}/*}}}*/
-void       Element::StrainRateHO2dvertical(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
-        /*Compute the 2d Blatter/HOStrain Rate (2 components):
-         *
-         * epsilon=[exx exz]
-         *
-         * with exz=1/2 du/dz
-         *
-         * the contribution of vz is neglected
-         */
-
-        /*Intermediaries*/
-        IssmDouble dvx[3];
-
-        /*Check that both inputs have been found*/
-        if (!vx_input){
-                _error_("Input missing. Here are the input pointers we have for vx: " << vx_input <<"\n");
-        }
-
-        /*Get strain rate assuming that epsilon has been allocated*/
-        vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
-        epsilon[0] = dvx[0];
-        epsilon[1] = 0.5*dvx[1];
-
-}/*}}}*/
-void       Element::StrainRateSSA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
-
-        /*Intermediaries*/
-        IssmDouble dvx[3];
-        IssmDouble dvy[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->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
-        vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
-        epsilon[0] = dvx[0];
-        epsilon[1] = dvy[1];
-        epsilon[2] = 0.5*(dvx[1] + dvy[0]);
-
-}/*}}}*/
-void       Element::StrainRateSSA1d(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input){/*{{{*/
-
-        /*Intermediaries*/
-        IssmDouble dvx[3];
-
-        /*Check that both inputs have been found*/
-        if (!vx_input){
-                _error_("Input missing. Here are the input pointers we have for vx: " << vx_input << "\n");
-        }
-
-        /*Get strain rate assuming that epsilon has been allocated*/
-        vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
-        *epsilon = dvx[0];
-
-}/*}}}*/
-void       Element::TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,int transformenum){/*{{{*/
-
-        /*All nodes have the same Coordinate System*/
-        int numnodes  = this->GetNumberOfNodes();
-        int* cs_array = xNew<int>(numnodes);
-        for(int i=0;i<numnodes;i++) cs_array[i]=transformenum;
-
-        /*Call core*/
-        TransformInvStiffnessMatrixCoord(Ke,this->nodes,numnodes,cs_array);
-
-        /*Clean-up*/
-        xDelete<int>(cs_array);
-}/*}}}*/
-void       Element::TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes_list,int numnodes,int* cs_array){/*{{{*/
-
-        int         i,j;
-        int         numdofs   = 0;
-        IssmDouble *transform = NULL;
-        IssmDouble *values    = NULL;
-
-        /*Get total number of dofs*/
-        for(i=0;i<numnodes;i++){
-                switch(cs_array[i]){
-                        case PressureEnum: numdofs+=1; break;
-                        case XYEnum:       numdofs+=2; break;
-                        case XYZEnum:      numdofs+=3; break;
-                        default: _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
-                }
-        }
-
-        /*Copy current stiffness matrix*/
-        values=xNew<IssmDouble>(Ke->nrows*Ke->ncols);
-        for(i=0;i<Ke->nrows;i++) for(j=0;j<Ke->ncols;j++) values[i*Ke->ncols+j]=Ke->values[i*Ke->ncols+j];
-
-        /*Get Coordinate Systems transform matrix*/
-        CoordinateSystemTransform(&transform,nodes_list,numnodes,cs_array);
-
-        /*Transform matrix: R*Ke*R^T */
-        TripleMultiply(transform,numdofs,numdofs,0,
-                                values,Ke->nrows,Ke->ncols,0,
-                                transform,numdofs,numdofs,1,
-                                &Ke->values[0],0);
-
-        /*Free Matrix*/
-        xDelete<IssmDouble>(transform);
-        xDelete<IssmDouble>(values);
-}/*}}}*/
-void       Element::TransformLoadVectorCoord(ElementVector* pe,int transformenum){/*{{{*/
-
-        /*All nodes have the same Coordinate System*/
-        int  numnodes = this->GetNumberOfNodes();
-        int* cs_array = xNew<int>(numnodes);
-        for(int i=0;i<numnodes;i++) cs_array[i]=transformenum;
-
-        /*Call core*/
-        this->TransformLoadVectorCoord(pe,this->nodes,numnodes,cs_array);
-
-        /*Clean-up*/
-        xDelete<int>(cs_array);
-}/*}}}*/
-void       Element::TransformLoadVectorCoord(ElementVector* pe,int* cs_array){/*{{{*/
-
-        this->TransformLoadVectorCoord(pe,this->nodes,this->GetNumberOfNodes(),cs_array);
-
-}/*}}}*/
-void       Element::TransformLoadVectorCoord(ElementVector* pe,Node** nodes_list,int numnodes,int* cs_array){/*{{{*/
-
-        int         i;
-        int         numdofs   = 0;
-        IssmDouble *transform = NULL;
-        IssmDouble *values    = NULL;
-
-        /*Get total number of dofs*/
-        for(i=0;i<numnodes;i++){
-                switch(cs_array[i]){
-                        case PressureEnum: numdofs+=1; break;
-                        case XYEnum:       numdofs+=2; break;
-                        case XYZEnum:      numdofs+=3; break;
-                        default: _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
-                }
-        }
-
-        /*Copy current load vector*/
-        values=xNew<IssmDouble>(pe->nrows);
-        for(i=0;i<pe->nrows;i++) values[i]=pe->values[i];
-
-        /*Get Coordinate Systems transform matrix*/
-        CoordinateSystemTransform(&transform,nodes_list,numnodes,cs_array);
-
-        /*Transform matrix: R^T*pe */
-        MatrixMultiply(transform,numdofs,numdofs,1,
-                                values,pe->nrows,1,0,
-                                &pe->values[0],0);
-
-        /*Free Matrices*/
-        xDelete<IssmDouble>(transform);
-        xDelete<IssmDouble>(values);
-}/*}}}*/
-void       Element::TransformSolutionCoord(IssmDouble* values,int transformenum){/*{{{*/
-
-        /*All nodes have the same Coordinate System*/
-        int  numnodes = this->GetNumberOfNodes();
-        int* cs_array = xNew<int>(numnodes);
-        for(int i=0;i<numnodes;i++) cs_array[i]=transformenum;
-
-        /*Call core*/
-        this->TransformSolutionCoord(values,this->nodes,numnodes,cs_array);
-
-        /*Clean-up*/
-        xDelete<int>(cs_array);
-}/*}}}*/
-void       Element::TransformSolutionCoord(IssmDouble* values,int* transformenum_list){/*{{{*/
-        this->TransformSolutionCoord(values,this->nodes,this->GetNumberOfNodes(),transformenum_list);
-}/*}}}*/
-void       Element::TransformSolutionCoord(IssmDouble* values,int numnodes,int transformenum){/*{{{*/
-
-        /*All nodes have the same Coordinate System*/
-        int* cs_array = xNew<int>(numnodes);
-        for(int i=0;i<numnodes;i++) cs_array[i]=transformenum;
-
-        /*Call core*/
-        this->TransformSolutionCoord(values,this->nodes,numnodes,cs_array);
-
-        /*Clean-up*/
-        xDelete<int>(cs_array);
-}/*}}}*/
-void       Element::TransformSolutionCoord(IssmDouble* solution,int numnodes,int* cs_array){/*{{{*/
-        this->TransformSolutionCoord(solution,this->nodes,numnodes,cs_array);
-}/*}}}*/
-void       Element::TransformSolutionCoord(IssmDouble* values,Node** nodes_list,int numnodes,int transformenum){/*{{{*/
-        /*NOT NEEDED*/
-        /*All nodes have the same Coordinate System*/
-        int* cs_array = xNew<int>(numnodes);
-        for(int i=0;i<numnodes;i++) cs_array[i]=transformenum;
-
-        /*Call core*/
-        this->TransformSolutionCoord(values,nodes_list,numnodes,cs_array);
-
-        /*Clean-up*/
-        xDelete<int>(cs_array);
-}/*}}}*/
-void       Element::TransformSolutionCoord(IssmDouble* solution,Node** nodes_list,int numnodes,int* cs_array){/*{{{*/
-
-        int         i;
-        int         numdofs   = 0;
-        IssmDouble *transform = NULL;
-        IssmDouble *values    = NULL;
-
-        /*Get total number of dofs*/
-        for(i=0;i<numnodes;i++){
-                switch(cs_array[i]){
-                        case PressureEnum: numdofs+=1; break;
-                        case XYEnum:       numdofs+=2; break;
-                        case XYZEnum:      numdofs+=3; break;
-                        default: _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
-                }
-        }
-
-        /*Copy current solution vector*/
-        values=xNew<IssmDouble>(numdofs);
-        for(i=0;i<numdofs;i++) values[i]=solution[i];
-
-        /*Get Coordinate Systems transform matrix*/
-        CoordinateSystemTransform(&transform,nodes_list,numnodes,cs_array);
-
-        /*Transform matrix: R*U */
-        MatrixMultiply(transform,numdofs,numdofs,0,
-                                values,numdofs,1,0,
-                                &solution[0],0);
-
-        /*Free Matrices*/
-        xDelete<IssmDouble>(transform);
-        xDelete<IssmDouble>(values);
-}/*}}}*/
-void       Element::TransformStiffnessMatrixCoord(ElementMatrix* Ke,int transformenum){/*{{{*/
-
-        /*All nodes have the same Coordinate System*/
-        int  numnodes = this->GetNumberOfNodes();
-        int* cs_array = xNew<int>(numnodes);
-        for(int i=0;i<numnodes;i++) cs_array[i]=transformenum;
-
-        /*Call core*/
-        this->TransformStiffnessMatrixCoord(Ke,this->nodes,numnodes,cs_array);
-
-        /*Clean-up*/
-        xDelete<int>(cs_array);
-}/*}}}*/
-void       Element::TransformStiffnessMatrixCoord(ElementMatrix* Ke,int* transformenum_list){/*{{{*/
-        this->TransformStiffnessMatrixCoord(Ke,this->nodes,this->GetNumberOfNodes(),transformenum_list);
-}/*}}}*/
-void       Element::TransformStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes_list,int numnodes,int* cs_array){/*{{{*/
-
-        int         numdofs = 0;
-        IssmDouble *transform = NULL;
-        IssmDouble *values    = NULL;
-
-        /*Get total number of dofs*/
-        for(int i=0;i<numnodes;i++){
-                switch(cs_array[i]){
-                        case PressureEnum: numdofs+=1; break;
-                        case XYEnum:       numdofs+=2; break;
-                        case XYZEnum:      numdofs+=3; break;
-                        default: _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
-                }
-        }
-
-        /*Copy current stiffness matrix*/
-        values=xNew<IssmDouble>(Ke->nrows*Ke->ncols);
-        for(int i=0;i<Ke->nrows*Ke->ncols;i++) values[i]=Ke->values[i];
-
-        /*Get Coordinate Systems transform matrix*/
-        CoordinateSystemTransform(&transform,nodes_list,numnodes,cs_array);
-
-        /*Transform matrix: R^T*Ke*R */
-        TripleMultiply(transform,numdofs,numdofs,1,
-                                values,Ke->nrows,Ke->ncols,0,
-                                transform,numdofs,numdofs,0,
-                                &Ke->values[0],0);
-
-        /*Free Matrix*/
-        xDelete<IssmDouble>(transform);
-        xDelete<IssmDouble>(values);
-}/*}}}*/
+        IssmDouble phi;
+        IssmDouble viscosity;
+        IssmDouble epsilon[6];
+        IssmDouble thickness;
+        IssmDouble *xyz_list = NULL;
+
+        /*Fetch number vertices and allocate memory*/
+        int         numvertices    = this->GetNumberOfVertices();
+        IssmDouble* viscousheating = xNew<IssmDouble>(numvertices);
+
+        /*Retrieve all inputs and parameters*/
+        this->GetVerticesCoordinatesBase(&xyz_list);
+        Input* vx_input        = this->GetInput(VxEnum); _assert_(vx_input);
+        Input* vy_input        = this->GetInput(VyEnum); _assert_(vy_input);
+        Input* vz_input        = this->GetInput(VzEnum); _assert_(vz_input);
+        Input* thickness_input = this->GetInput(ThicknessEnum); _assert_(thickness_input);
+
+        /*loop over vertices: */
+        Gauss* gauss=this->NewGauss();
+        for (int iv=0;iv<numvertices;iv++){
+                gauss->GaussVertex(iv);
+
+                thickness_input->GetInputValue(&thickness,gauss);
+
+                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->GetPhi(&phi,&epsilon[0],viscosity);
+
+                viscousheating[iv]=phi*thickness;
+        }
+
+        /*Create PentaVertex input, which will hold the basal friction:*/
+        this->AddInput(ViscousHeatingEnum,viscousheating,P1Enum);
+
+        /*Clean up and return*/
+        xDelete<IssmDouble>(viscousheating);
+        delete gauss;
+}
+/*}}}*/