Index: /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp	(revision 16896)
+++ /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp	(revision 16897)
@@ -231,5 +231,5 @@
 				for(int j=0;j<numnodes;j++){
 					Ke->values[i*numnodes+j]+=tau_parameter*D_scalar*
-					  ((u-um)*dbasis[0*3+i]+(v-vm)*dbasis[1*3+i]+(w-wm)*dbasis[2*3+i])*((u-um)*dbasis[0*3+j]+(v-vm)*dbasis[1*3+j]+(w-wm)*dbasis[2*3+j]);
+					  ((u-um)*dbasis[0*numnodes+i]+(v-vm)*dbasis[1*numnodes+i]+(w-wm)*dbasis[2*numnodes+i])*((u-um)*dbasis[0*numnodes+j]+(v-vm)*dbasis[1*numnodes+j]+(w-wm)*dbasis[2*numnodes+j]);
 				}
 			}
@@ -237,5 +237,5 @@
 				for(int i=0;i<numnodes;i++){
 					for(int j=0;j<numnodes;j++){
-						Ke->values[i*numnodes+j]+=tau_parameter*D_scalar*basis[j]*((u-um)*dbasis[0*3+i]+(v-vm)*dbasis[1*3+i]+(w-wm)*dbasis[2*3+i]);
+						Ke->values[i*numnodes+j]+=tau_parameter*D_scalar*basis[j]*((u-um)*dbasis[0*numnodes+i]+(v-vm)*dbasis[1*numnodes+i]+(w-wm)*dbasis[2*numnodes+i]);
 					}
 				}
@@ -377,7 +377,7 @@
 			tau_parameter=element->StabilizationParameter(u,v,w,diameter,kappa);
 
-			for(int i=0;i<numnodes;i++) pe->values[i]+=tau_parameter*scalar_def*(u*dbasis[0*3+i]+v*dbasis[1*3+i]+w*dbasis[2*3+i]);
+			for(int i=0;i<numnodes;i++) pe->values[i]+=tau_parameter*scalar_def*(u*dbasis[0*numnodes+i]+v*dbasis[1*numnodes+i]+w*dbasis[2*numnodes+i]);
 			if(reCast<bool,IssmDouble>(dt)){
-				for(int i=0;i<numnodes;i++) pe->values[i]+=tau_parameter*scalar_transient*(u*dbasis[0*3+i]+v*dbasis[1*3+i]+w*dbasis[2*3+i]);
+				for(int i=0;i<numnodes;i++) pe->values[i]+=tau_parameter*scalar_transient*(u*dbasis[0*numnodes+i]+v*dbasis[1*numnodes+i]+w*dbasis[2*numnodes+i]);
 			}
 		}
@@ -393,5 +393,61 @@
 }/*}}}*/
 ElementVector* ThermalAnalysis::CreatePVectorSheet(Element* element){/*{{{*/
-	return NULL;
+
+	/* Geothermal flux on ice sheet base and basal friction */
+	if(!element->IsOnBed() || element->IsFloating()) return NULL;
+
+	IssmDouble  dt,Jdet,geothermalflux,vx,vy,vz;
+	IssmDouble  alpha2,scalar,basalfriction,heatflux;
+	IssmDouble *xyz_list_base = NULL;
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Initialize vectors*/
+	ElementVector* pe    = element->NewElementVector();
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinatesBase(&xyz_list_base);
+	element->FindParam(&dt,TimesteppingTimeStepEnum);
+	Input* vx_input             = element->GetInput(VxEnum);                          _assert_(vx_input);
+	Input* vy_input             = element->GetInput(VyEnum);                          _assert_(vy_input);
+	Input* vz_input             = element->GetInput(VzEnum);                          _assert_(vz_input);
+	Input* geothermalflux_input = element->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input);
+	IssmDouble  rho_ice             = element->GetMaterialParameter(MaterialsRhoIceEnum);
+	IssmDouble  heatcapacity        = element->GetMaterialParameter(MaterialsHeatcapacityEnum);
+
+	/*Build friction element, needed later: */
+	Friction* friction=new Friction(element,3);
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss   = element->NewGaussBase(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
+		element->NodalFunctions(basis,gauss);
+
+		geothermalflux_input->GetInputValue(&geothermalflux,gauss);
+		friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,vz_input);
+		vx_input->GetInputValue(&vx,gauss);
+		vy_input->GetInputValue(&vy,gauss);
+		vz_input->GetInputValue(&vz,gauss);
+		vz = 0.;//FIXME
+		basalfriction = alpha2*(vx*vx + vy*vy + vz*vz);
+		heatflux      = (basalfriction+geothermalflux)/(rho_ice*heatcapacity);
+
+		scalar = gauss->weight*Jdet*heatflux;
+		if(dt!=0.) scalar=dt*scalar;
+
+		for(int i=0;i<numnodes;i++) pe->values[i]+=scalar*basis[i];
+	}
+
+	/*Clean up and return*/
+	delete gauss;
+	delete friction;
+	xDelete<IssmDouble>(basis);
+	xDelete<IssmDouble>(xyz_list_base);
+	return pe;
 }/*}}}*/
 ElementVector* ThermalAnalysis::CreatePVectorShelf(Element* element){/*{{{*/