Index: /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp	(revision 22263)
+++ /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp	(revision 22264)
@@ -71,6 +71,6 @@
 
 	/*Intermediaries */
-	int         n=3;
-	IssmDouble  Jdet,omega,ds[2],slope;
+	IssmDouble  yts = 365*24*3600.;
+	IssmDouble  Jdet,vx,vy,vel;
 	IssmDouble* xyz_list = NULL;
 
@@ -84,8 +84,18 @@
 	/*Retrieve all inputs and parameters*/
 	element->GetVerticesCoordinates(&xyz_list);
-	Input* omega_input         = element->GetInput(BalancethicknessOmegaEnum); _assert_(omega_input);
-	Input* surfaceslopex_input = element->GetInput(SurfaceSlopeXEnum); _assert_(surfaceslopex_input); 
-	Input* surfaceslopey_input = element->GetInput(SurfaceSlopeYEnum); _assert_(surfaceslopey_input); 
-	IssmDouble rhog    = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum);
+	Input* vx_input = element->GetInput(SurfaceSlopeXEnum); _assert_(vx_input); 
+	Input* vy_input = element->GetInput(SurfaceSlopeYEnum); _assert_(vy_input);
+
+	/*make sure are diffusivisty is large enough*/
+	vel = sqrt(vx*vx+vy*vy);
+	if(sqrt(vx*vx+vy*vy)==0.){
+		vx = 0.1/yts;
+		vy = 0.1/yts;
+	}
+	else if(vel<0.1/yts){
+		vx = vx/vel*0.1;
+		vy = vy/vel*0.1;
+
+	}
 
 	/* Start  looping on the number of gaussian points: */
@@ -95,13 +105,10 @@
 		element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-		omega_input->GetInputValue(&omega,gauss);
-		surfaceslopex_input->GetInputValue(&ds[0],gauss);
-		surfaceslopey_input->GetInputValue(&ds[1],gauss);
-
-		slope = sqrt(ds[0]*ds[0] + ds[1]*ds[1]);
+		vx_input->GetInputValue(&vx,gauss);
+		vy_input->GetInputValue(&vy,gauss);
 
 		for(int i=0;i<numnodes;i++){
 			for(int j=0;j<numnodes;j++){
-				Ke->values[i*numnodes+j] += pow(rhog,n)*exp(omega)*pow(slope,n-1)*gauss->weight*Jdet*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j]);
+				Ke->values[i*numnodes+j] += gauss->weight*Jdet*(vx*dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + vy*dbasis[1*numnodes+i]*dbasis[1*numnodes+j]);
 			}
 		}
@@ -117,5 +124,5 @@
 
 	/*Intermediaries */
-	IssmDouble  dhdt,mb,ms,Jdet;
+	IssmDouble  dhdt[2],mb[2],ms[2],Jdet;
 	IssmDouble* xyz_list = NULL;
 
@@ -141,10 +148,12 @@
 		element->NodalFunctions(basis,gauss);
 
-		ms_input->GetInputValue(&ms,gauss);
-		mb_input->GetInputValue(&mb,gauss);
-		dhdt_input->GetInputValue(&dhdt,gauss);
+		ms_input->GetInputDerivativeValue(&ms[0],xyz_list,gauss);
+
+		ms_input->GetInputDerivativeValue(&ms[0],xyz_list,gauss);
+		mb_input->GetInputDerivativeValue(&mb[0],xyz_list,gauss);
+		dhdt_input->GetInputDerivativeValue(&dhdt[0],xyz_list,gauss);
 
 		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(
-					(ms-mb-dhdt)*basis[i]
+					(ms[0]+ms[1]-mb[0]-mb[1]-dhdt[0]-dhdt[1])*basis[i]
 					);
 	}
@@ -157,5 +166,5 @@
 }/*}}}*/
 void           Balancethickness2Analysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
-		element->GetSolutionFromInputsOneDof(solution,SurfaceEnum);
+		element->GetSolutionFromInputsOneDof(solution,ThicknessEnum);
 }/*}}}*/
 void           Balancethickness2Analysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
@@ -164,55 +173,6 @@
 void           Balancethickness2Analysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
 
-	/*Intermediaries*/
-	IssmDouble  ds[2],s,b,D;
-	IssmDouble* xyz_list = NULL;
+			element->InputUpdateFromSolutionOneDof(solution,ThicknessEnum);
 
-	//element->InputUpdateFromSolutionOneDof(solution,ThicknessEnum);
-	element->InputUpdateFromSolutionOneDof(solution,SurfaceEnum);
-
-	/*Fetch number of vertices and allocate velocity vectors*/
-	int numvertices = element->GetNumberOfVertices();
-	IssmDouble* vel_list = xNew<IssmDouble>(numvertices);
-	IssmDouble* vx_list  = xNew<IssmDouble>(numvertices);
-	IssmDouble* vy_list  = xNew<IssmDouble>(numvertices);
-
-	/*Retrieve all inputs and parameters*/
-	element->GetVerticesCoordinates(&xyz_list);
-	Input* D_input   = element->GetInput(BalancethicknessDiffusionCoefficientEnum);
-	Input* H_input   = element->GetInput(ThicknessEnum);                            _assert_(H_input);
-	Input* s_input   = element->GetInput(SurfaceEnum);                              _assert_(s_input);
-	Input* b_input   = element->GetInput(BaseEnum);                                 _assert_(b_input);
-
-	/*Calculate velocities*/
-	Gauss* gauss=element->NewGauss();
-	for(int iv=0;iv<numvertices;iv++){
-		gauss->GaussVertex(iv);
-
-		if(D_input){
-			D_input->GetInputValue(&D,gauss);
-		}
-		else{
-			D = 0.;
-		}
-		b_input->GetInputValue(&b,gauss);
-		s_input->GetInputValue(&s,gauss);
-		s_input->GetInputDerivativeValue(&ds[0],xyz_list,gauss);
-
-		vx_list[iv] = -1./(s-b)*D*ds[0];
-		vy_list[iv] = -1./(s-b)*D*ds[1];
-		vel_list[iv] = sqrt(pow(vx_list[iv],2) + pow(vy_list[iv],2));
-	}
-
-	/*Add vx and vy as inputs to the tria element: */
-	element->AddInput(VxEnum,vx_list,P1Enum);
-	element->AddInput(VyEnum,vy_list,P1Enum);
-	element->AddInput(VelEnum,vel_list,P1Enum);
-
-	/*Free ressources:*/
-	delete gauss;
-	xDelete<IssmDouble>(vy_list);
-	xDelete<IssmDouble>(vx_list);
-	xDelete<IssmDouble>(vel_list);
-	xDelete<IssmDouble>(xyz_list);
 }/*}}}*/
 void           Balancethickness2Analysis::UpdateConstraints(FemModel* femmodel){/*{{{*/