Index: /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp	(revision 18544)
+++ /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp	(revision 18545)
@@ -30,4 +30,5 @@
 	iomodel->FetchDataToInput(elements,BaseEnum);
 	iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
+	iomodel->FetchDataToInput(elements,FrictionCoefficientEnum);
 	iomodel->FetchDataToInput(elements,BasalforcingsGroundediceMeltingRateEnum);
 	iomodel->FetchDataToInput(elements,SurfaceforcingsMassBalanceEnum);
@@ -105,5 +106,5 @@
 
 	/*Intermediaries */
-	IssmDouble  dhdt,mb,ms,dD[2],db[2],Jdet;
+	IssmDouble  dhdt,mb,ms,D,db[2],Jdet;
 	IssmDouble* xyz_list = NULL;
 
@@ -114,4 +115,5 @@
 	ElementVector* pe    = element->NewElementVector();
 	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+	IssmDouble*   dbasis = xNew<IssmDouble>(2*numnodes);
 
 	/*Retrieve all inputs and parameters*/
@@ -130,15 +132,18 @@
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
 		element->NodalFunctions(basis,gauss);
+		element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 
 		ms_input->GetInputValue(&ms,gauss);
 		mb_input->GetInputValue(&mb,gauss);
 		dhdt_input->GetInputValue(&dhdt,gauss);
+		D_input->GetInputValue(&D,gauss);
 		bed_input->GetInputDerivativeValue(&db[0],xyz_list,gauss);
-		D_input->GetInputDerivativeValue(&dD[0],xyz_list,gauss);
-		db[0]=0.;
-		db[1]=0.;
-
-		/*Since grad(b) is constant div(D grad(b) ) = grad(D).grad(b)*/
-		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(ms-mb-dhdt +dD[0]*db[0]+dD[1]*db[1])*basis[i];
+		//db[0]=0.;
+		//db[1]=0.;
+
+		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(
+					(ms-mb-dhdt)*basis[i]
+					-D*db[0]*dbasis[0*numnodes+i] - D*db[1]*dbasis[1*numnodes+i]
+					);
 	}
 
@@ -167,6 +172,6 @@
 
 	/*Intermediaries */
-	IssmDouble       omega,h,mu0,ds[2],Cmu,B;
-	const int        n = 3.;
+	IssmDouble       Gamma,h,mu0,ds[2],Cmu,B,k;
+	const int        n = 3;
 	const IssmDouble Hstar = 500.;
 	const IssmDouble Lstar = 500.e+3;
@@ -183,4 +188,5 @@
 	Input* surface_input    = element->GetInput(SurfaceEnum);            _assert_(surface_input);
 	Input* B_input          = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input);
+	Input* k_input          = element->GetInput(FrictionCoefficientEnum);_assert_(k_input);
 	IssmDouble rhog = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum);
 
@@ -191,12 +197,13 @@
 		
 		B_input->GetInputValue(&B,gauss);
+		k_input->GetInputValue(&k,gauss);
 		thickness_input->GetInputValue(&h,gauss);
 		surface_input->GetInputDerivativeValue(&ds[0],xyz_list,gauss);
 
-		mu0   = pow(2.,(1.-3*n)/(2.*n)) * B;
-		omega = pow(rhog,3) * pow(Hstar,2*(n+1)) * pow(Hstar/Lstar,2*(n+1)) * 1./pow(mu0,n);
-		Cmu   = 0.;
-
-		D[i] = omega*(ds[0]*ds[0]+ds[1]*ds[1])*pow(h,4)*(1./5.*h - Cmu);
+		mu0   = pow(2.,(1-3*n)/(2.*n)) * B;
+		Gamma = pow(rhog,n) * pow(Hstar,2*(n+1)) * pow(Hstar/Lstar,2*(n+1)) * 1./pow(mu0,n);
+		Cmu   = mu0/Hstar/(k*k);
+
+		D[i] = Gamma*pow(ds[0]*ds[0]+ds[1]*ds[1],(n-1)/2)*pow(h,n+1)*(1./(n+2.)*h - Cmu);
 	}