Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17538)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17539)
@@ -234,4 +234,5 @@
 		iomodel->FetchDataToInput(elements,PressureEnum,0.);
 		if(dakota_analysis)elements->InputDuplicate(PressureEnum,QmuPressureEnum);
+		if(finiteelement==XTaylorHoodEnum) this->InitializeXTH(elements);
 	}
 	if(islevelset){
@@ -820,4 +821,53 @@
 	}
 }/*}}}*/
+void StressbalanceAnalysis::InitializeXTH(Elements* elements){/*{{{*/
+
+	/*Intermediaries*/
+	int        meshtype,dim;
+	IssmDouble dvx[3],dvy[3],dvz[3];
+	IssmDouble viscosity;
+	IssmDouble *xyz_list = NULL;
+
+	for(int i=0;i<elements->Size();i++){
+		Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
+
+		/*Get problem dimension*/
+		element->FindParam(&meshtype,MeshTypeEnum);
+		switch(meshtype){
+			case Mesh2DverticalEnum: dim = 2; break;
+			case Mesh3DEnum:         dim = 3; break;
+			case Mesh3DtetrasEnum:   dim = 3; break;
+			default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+		}
+
+		int tnumnodes = element->GetNumberOfVertices();      //Tensors, P1 DG
+
+		/*Get inputs and parameters*/
+		element->GetVerticesCoordinates(&xyz_list);
+		Input* vx_input=element->GetInput(VxEnum); _assert_(vx_input);
+		Input* vy_input=element->GetInput(VyEnum); _assert_(vy_input);
+		Input* vz_input;
+		if(dim==3){vz_input=element->GetInput(VzEnum); _assert_(vz_input);}
+
+		/*Get d and tau*/
+		Gauss* gauss = element->NewGauss();
+		for(int i=0;i<tnumnodes;i++){
+			gauss->GaussNode(P1DGEnum,i);
+
+			vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
+			vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
+			if(dim==3){
+				vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss);
+			}
+
+			_error_("STOP");
+			element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
+		}
+
+		/*Clean up*/
+		xDelete<IssmDouble>(xyz_list);
+	}
+
+}/*}}}*/
 
 /*Finite Element Analysis*/
Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 17538)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 17539)
@@ -19,4 +19,5 @@
 		void CreateConstraints(Constraints* constraints,IoModel* iomodel);
 		void CreateLoads(Loads* loads, IoModel* iomodel);
+		void InitializeXTH(Elements* elements);
 
 		/*Finite element Analysis*/