Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16995)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16996)
@@ -4520,109 +4520,4 @@
 }
 /*}}}*/
-/*FUNCTION Penta::KMatrixGLSstabilization{{{*/
-void Penta::KMatrixGLSstabilization(ElementMatrix* Ke){
-
-	int        numdof  = NUMVERTICES*NDOF4;
-
-	/*Intermediaries */
-	int        i,j;
-	IssmDouble Jdet,viscosity,FSreconditioning,diameter,rigidity;
-	IssmDouble xyz_list[NUMVERTICES][3];
-	IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
-	GaussPenta *gauss=NULL;
-
-	/*Stabilization*/
-	IssmDouble D_scalar;
-	IssmDouble dbasis[3][6];
-	IssmDouble dmu[3];
-	IssmDouble mu;
-	IssmDouble dnodalbasis[6][6][3];
-	IssmDouble SU[6][4][4];
-	IssmDouble SW[6][4][4];
-	int p,q,ii;
-	int c=3; //index of pressure
-
-	/*Retrieve all inputs and parameters*/
-	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-	parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
-	Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
-	Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
-	Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input);
-
-	/*Find minimal length and B*/
-	rigidity=material->GetB();
-	diameter=MinEdgeLength(&xyz_list[0][0]);
-
-		gauss=new GaussPenta();
-		for(int iv=0;iv<6;iv++){
-			gauss->GaussVertex(iv);
-			GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
-			for(i=0;i<6;i++){
-				for(j=0;j<3;j++){
-					dnodalbasis[i][iv][j] = dbasis[j][i];
-				}
-			}
-		}
-		delete gauss;
-
-	/* Start  looping on the number of gaussian points: */
-	gauss=new GaussPenta(5,5);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-
-		gauss->GaussPoint(ig);
-
-		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
-
-		this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
-		material->GetViscosity3dFS(&viscosity,&epsilon[0]);
-
-		D_scalar=gauss->weight*Jdet;
-
-		/*Add stabilization*/
-		GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
-		dmu[0]=0.; dmu[1]=0.; dmu[2]=0.;
-		mu = 2.*viscosity;
-		for(p=0;p<6;p++) for(i=0;i<4;i++) for(j=0;j<4;j++){
-			SU[p][i][j]=0.;
-			SW[p][i][j]=0.;
-		}
-		for(p=0;p<6;p++){
-			for(i=0;i<3;i++){
-				SU[p][i][c] += dbasis[i][p];
-				SW[p][c][i] += dbasis[i][p];
-				for(j=0;j<3;j++){
-					SU[p][i][i] += -dmu[j]*dbasis[j][p];
-					SU[p][i][j] += -dmu[i]*dbasis[i][p];
-					for(ii=0;ii<6;ii++){
-						SU[p][i][i] += -mu*dnodalbasis[p][ii][j]*dbasis[j][ii];
-						SU[p][i][j] += -mu*dnodalbasis[p][ii][i]*dbasis[j][ii];
-					}
-					SW[p][i][i] += -dmu[j]*dbasis[j][p];
-					SW[p][j][i] += -dmu[j]*dbasis[i][p];
-					for(ii=0;ii<6;ii++){
-						SW[p][i][i] += -mu*dnodalbasis[p][ii][j]*dbasis[j][ii];
-						SW[p][j][i] += -mu*dnodalbasis[p][ii][i]*dbasis[j][ii];
-					}
-				}
-			}
-		}
-		IssmDouble tau = 1./3.*pow(diameter,2)/(8.*2.*viscosity);
-		for(p=0;p<6;p++){
-			for(q=0;q<6;q++){
-				for(i=0;i<4;i++){
-					for(j=0;j<4;j++){
-						Ke->values[(p*4+i)*numdof+q*4+j] += gauss->weight*Jdet*tau*SW[p][i][0]*SU[q][0][j];
-						Ke->values[(p*4+i)*numdof+q*4+j] += gauss->weight*Jdet*tau*SW[p][i][1]*SU[q][1][j];
-						Ke->values[(p*4+i)*numdof+q*4+j] += gauss->weight*Jdet*tau*SW[p][i][2]*SU[q][2][j];
-					}
-				}
-			}
-		}
-	}
-
-	/*Clean up*/
-	delete gauss;
-}
-/*}}}*/
 #endif
 
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 16995)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 16996)
@@ -212,6 +212,4 @@
 		void           GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype);
 		Node*          GetNode(int node_number);
-		void           GetStrainRate3dHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss, Input* vx_input, Input* vy_input);
-		void           GetStrainRate3d(IssmDouble* epsilon,IssmDouble* xyz_list, Gauss* gauss, Input* vx_input, Input* vy_input, Input* vz_input);
 		Penta*         GetUpperElement(void);
 		Penta*         GetLowerElement(void);
@@ -257,5 +255,4 @@
 		ElementVector* CreateDVectorStressbalanceHoriz(void);
 		ElementVector* CreateDVectorStressbalanceFS(void);
-		void           KMatrixGLSstabilization(ElementMatrix* Ke);
 		#endif
 
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16995)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16996)
@@ -3062,21 +3062,4 @@
 #endif
 
-#ifdef _HAVE_STRESSBALANCE_
-/*FUNCTION Tria::GetYcoord {{{*/
-IssmDouble Tria::GetYcoord(GaussTria* gauss){
-
-	IssmDouble y;
-	IssmDouble xyz_list[NUMVERTICES][3];
-	IssmDouble y_list[NUMVERTICES];
-
-	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-	for(int i=0;i<NUMVERTICES;i++) y_list[i]=xyz_list[i][1];
-	TriaRef::GetInputValue(&y,&y_list[0],gauss,P1Enum);
-
-	return y;
-}
-/*}}}*/
-#endif
-
 #ifdef _HAVE_CONTROL_
 /*FUNCTION Tria::BalancethicknessMisfit{{{*/
@@ -4905,73 +4888,4 @@
 #endif
 
-#ifdef _HAVE_DAMAGE_
-/*FUNCTION Tria::DamageEvolutionF{{{*/
-void Tria::DamageEvolutionF(IssmDouble* f){
-
-	/*Intermediaries */
-	IssmDouble c1,c2,c3,healing,stress_threshold;
-	IssmDouble s_xx,s_xy,s_yy;
-	IssmDouble J2s;
-	IssmDouble Xis;
-	IssmDouble Psi;
-	IssmDouble PosPsi;
-	IssmDouble NegPsi;
-	Input* damage_input=NULL;
-	Input* sigma_xx_input  = NULL;
-	Input* sigma_xy_input  = NULL;
-	Input* sigma_yy_input  = NULL;
-	GaussTria* gauss=NULL;
-	IssmDouble damage,sigma_xx,sigma_xy,sigma_yy;
-
-	/*retrieve parameters:*/
-	this->parameters->FindParam(&c1,DamageC1Enum);
-	this->parameters->FindParam(&c2,DamageC2Enum);
-	this->parameters->FindParam(&c3,DamageC3Enum);
-	this->parameters->FindParam(&healing,DamageHealingEnum);
-	this->parameters->FindParam(&stress_threshold,DamageStressThresholdEnum);
-
-	/*Compute stress tensor: */
-	this->ComputeStressTensor();
-
-	/*retrieve what we need: */
-	sigma_xx_input  = inputs->GetInput(StressTensorxxEnum);  _assert_(sigma_xx_input);
-	sigma_xy_input  = inputs->GetInput(StressTensorxyEnum);  _assert_(sigma_xy_input);
-	sigma_yy_input  = inputs->GetInput(StressTensoryyEnum);  _assert_(sigma_yy_input);
-	damage_input  = this->material->inputs->GetInput(DamageDbarEnum);        _assert_(damage_input);
-
-	/*Damage evolution z mapping: */
-	gauss=new GaussTria();
-	J2s=0;
-	for (int iv=0;iv<NUMVERTICES;iv++){
-		gauss->GaussVertex(iv);
-		
-		damage_input->GetInputValue(&damage,gauss);
-		sigma_xx_input->GetInputValue(&sigma_xx,gauss);
-		sigma_xy_input->GetInputValue(&sigma_xy,gauss);
-		sigma_yy_input->GetInputValue(&sigma_yy,gauss);
-
-		s_xx=sigma_xx/(1-damage);
-		s_xy=sigma_xy/(1-damage);
-		s_yy=sigma_yy/(1-damage);
-
-		J2s=1.0/sqrt(2.0)*sqrt(pow(s_xx,2)+pow(s_yy,2)+2*pow(s_xy,2));
-		
-		Xis=sqrt(3.0)*J2s;
-
-		Psi=Xis-stress_threshold;
-
-		PosPsi=max(Psi,0.0);
-		NegPsi=max(-Psi,0.0);
-
-		f[iv]= c1* ( pow(PosPsi,c2)   -  healing * pow(NegPsi,c2) )  *   pow((1 - damage),-c3);
-
-	}
-	
-	/*Clean up and return*/
-	delete gauss;
-}
-/*}}}*/
-#endif
-
 #ifdef _HAVE_DAKOTA_
 /*FUNCTION Tria::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type);{{{*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 16995)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 16996)
@@ -261,9 +261,4 @@
 		void           ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not implemented yet");};
 
-		#ifdef _HAVE_STRESSBALANCE_
-		void           PVectorGLSstabilization(ElementVector* pe);
-		IssmDouble     GetYcoord(GaussTria* gauss);
-		#endif
-
 		void UpdateConstraintsExtrudeFromBase(void);
 		void UpdateConstraintsExtrudeFromTop(void);
@@ -283,12 +278,5 @@
 		void           HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask);
 		void           ComputeEPLThickness(void);
-
-		#endif
-
-		#ifdef _HAVE_DAMAGE_
-		void           DamageEvolutionF(IssmDouble* flist);
-		#endif
-
-
+		#endif
 		/*}}}*/