Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17525)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17526)
@@ -3387,8 +3387,7 @@
 }/*}}}*/
 ElementVector* StressbalanceAnalysis::CreatePVectorFSViscousXTH(Element* element){/*{{{*/
-	_error_("STOP");
-
-	int         i,meshtype,dim;
-	IssmDouble  Jdet,forcex,forcey,forcez;
+
+	int         i,tausize,meshtype,dim;
+	IssmDouble  Jdet,r;
 	IssmDouble *xyz_list = NULL;
 
@@ -3396,13 +3395,17 @@
 	element->FindParam(&meshtype,MeshTypeEnum);
 	switch(meshtype){
-		case Mesh2DverticalEnum: dim = 2; break;
-		case Mesh3DEnum:         dim = 3; break;
-		case Mesh3DtetrasEnum:   dim = 3; break;
+		case Mesh2DverticalEnum: dim = 2; tausize = 3; break;
+		case Mesh3DEnum:         dim = 3; tausize = 6; break;
+		case Mesh3DtetrasEnum:   dim = 3; tausize = 6; break;
 		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
 	}
+
+	/*FIXME*/
+	r = 1.;
 
 	/*Fetch number of nodes and dof for this finite element*/
 	int vnumnodes = element->NumberofNodesVelocity();
 	int pnumnodes = element->NumberofNodesPressure();
+	int tnumnodes = element->GetNumberOfVertices();      //Tensors, P1 DG
 
 	/*Prepare coordinate system list*/
@@ -3413,17 +3416,10 @@
 
 	/*Initialize vectors*/
-	ElementVector* pe     = element->NewElementVector(FSvelocityEnum);
-	IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
+	ElementVector* pe    = element->NewElementVector(FSvelocityEnum);
+	IssmDouble*    Dstar = xNew<IssmDouble>(tausize*tnumnodes);
+	IssmDouble*    tau   = xNew<IssmDouble>(tausize);
 
 	/*Retrieve all inputs and parameters*/
 	element->GetVerticesCoordinates(&xyz_list);
-	IssmDouble  rho_ice =element->GetMaterialParameter(MaterialsRhoIceEnum);
-	IssmDouble  gravity =element->GetMaterialParameter(ConstantsGEnum);
-	Input*      loadingforcex_input=element->GetInput(LoadingforceXEnum);  _assert_(loadingforcex_input);
-	Input*      loadingforcey_input=element->GetInput(LoadingforceYEnum);  _assert_(loadingforcey_input);
-	Input*      loadingforcez_input=NULL;
-	if(dim==3){
-		loadingforcez_input=element->GetInput(LoadingforceZEnum);  _assert_(loadingforcez_input);
-	}
 
 	/* Start  looping on the number of gaussian points: */
@@ -3431,23 +3427,6 @@
 	for(int ig=gauss->begin();ig<gauss->end();ig++){
 		gauss->GaussPoint(ig);
-
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-		element->NodalFunctionsVelocity(vbasis,gauss);
-
-		loadingforcex_input->GetInputValue(&forcex,gauss);
-		loadingforcey_input->GetInputValue(&forcey,gauss);
-		if(dim==3) loadingforcez_input->GetInputValue(&forcez,gauss);
-
-		for(i=0;i<vnumnodes;i++){
-			pe->values[i*dim+0] += +rho_ice*forcex *Jdet*gauss->weight*vbasis[i];
-			pe->values[i*dim+1] += +rho_ice*forcey *Jdet*gauss->weight*vbasis[i];
-			if(dim==3){
-				pe->values[i*dim+2] += +rho_ice*forcez*Jdet*gauss->weight*vbasis[i];
-				pe->values[i*dim+2] += -rho_ice*gravity*Jdet*gauss->weight*vbasis[i];
-			}
-			else{
-				pe->values[i*dim+1] += -rho_ice*gravity*Jdet*gauss->weight*vbasis[i];
-			}
-		}
+		_error_("STOP");
 	}
 
@@ -3458,5 +3437,4 @@
 	delete gauss;
 	xDelete<int>(cs_list);
-	xDelete<IssmDouble>(vbasis);
 	xDelete<IssmDouble>(xyz_list);
 	return pe;