Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17108)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17109)
@@ -5,4 +5,6 @@
 #include "../modules/modules.h"
 #include "../solutionsequences/solutionsequences.h"
+
+#define FSANALYTICAL 1
 
 /*Model processing*/
@@ -2719,4 +2721,8 @@
 ElementVector* StressbalanceAnalysis::CreatePVectorFS(Element* element){/*{{{*/
 
+	#ifdef FSANALYTICAL
+	ElementVector* pe=CreatePVectorFSAnalytical(element);
+
+	#else
 	/*compute all load vectors for this element*/
 	ElementVector* pe1=CreatePVectorFSViscous(element);
@@ -2729,4 +2735,68 @@
 	delete pe2;
 	delete pe3;
+	#endif
+
+	return pe;
+}/*}}}*/
+ElementVector* StressbalanceAnalysis::CreatePVectorFSAnalytical(Element* element){/*{{{*/
+
+	bool        analytical_solution=0;
+	int         i,meshtype,dim;
+	IssmDouble  x_coord,y_coord,z_coord;
+	IssmDouble  Jdet,forcex,forcey,forcez;
+	IssmDouble *xyz_list = NULL;
+
+	element->FindParam(&meshtype,MeshTypeEnum);
+	switch(meshtype){
+		case Mesh3DEnum:         dim = 3; break;
+		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+	}
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int vnumnodes = element->GetNumberOfNodesVelocity();
+	int pnumnodes = element->GetNumberOfNodesPressure();
+
+	/*Prepare coordinate system list*/
+	int* cs_list = xNew<int>(vnumnodes+pnumnodes);
+	if(dim==2) for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
+	else       for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum;
+	for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
+
+	/*Initialize vectors*/
+	ElementVector* pe     = element->NewElementVector(FSvelocityEnum);
+	IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinates(&xyz_list);
+	IssmDouble  rho_ice =element->GetMaterialParameter(MaterialsRhoIceEnum);
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=element->NewGauss(5);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		element->NodalFunctionsVelocity(vbasis,gauss);
+
+		forcex=1.0; forcey=1.0; forcez=1.0;
+//		forcex->fx(x_coord,y_coord,z_coord);
+//		forcex->fy(x_coord,y_coord,z_coord);
+//		forcex->fz(x_coord,y_coord,z_coord);
+
+		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];
+			pe->values[i*dim+2] += +rho_ice*forcez *Jdet*gauss->weight*vbasis[i];
+		}
+	}
+
+	/*Transform coordinate system*/
+	element->TransformLoadVectorCoord(pe,cs_list);
+
+	/*Clean up and return*/
+	delete gauss;
+	xDelete<int>(cs_list);
+	xDelete<IssmDouble>(vbasis);
+	xDelete<IssmDouble>(xyz_list);
 	return pe;
 }/*}}}*/
Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 17108)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 17109)
@@ -69,4 +69,5 @@
 		ElementMatrix* CreateKMatrixFSFriction(Element* element);
 		ElementVector* CreatePVectorFS(Element* element);
+		ElementVector* CreatePVectorFSAnalytical(Element* element);
 		ElementVector* CreatePVectorFSViscous(Element* element);
 		ElementVector* CreatePVectorFSShelf(Element* element);
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 17108)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 17109)
@@ -166,6 +166,7 @@
 		virtual void   GetVerticesCoordinatesBase(IssmDouble** xyz_list)=0;
 		virtual void   GetVerticesCoordinatesTop(IssmDouble** xyz_list)=0;
+		virtual IssmDouble GetXcoord(Gauss* gauss)=0;
+		virtual IssmDouble GetYcoord(Gauss* gauss)=0;
 		virtual IssmDouble GetZcoord(Gauss* gauss)=0;
-		virtual IssmDouble GetYcoord(Gauss* gauss)=0;
 		virtual void   GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype)=0;
 		virtual int    GetElementType(void)=0;
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 17108)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 17109)
@@ -1207,4 +1207,34 @@
 	/*We found the enum.  Use its values to fill into the vector, using the vertices ids: */
 	input->GetVectorFromInputs(vector,&vertexpidlist[0]);
+}
+/*}}}*/
+/*FUNCTION Penta::GetXcoord {{{*/
+IssmDouble Penta::GetXcoord(Gauss* gauss){
+
+	int    i;
+	IssmDouble x;
+	IssmDouble xyz_list[NUMVERTICES][3];
+	IssmDouble x_list[NUMVERTICES];
+
+	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+	for(i=0;i<NUMVERTICES;i++) x_list[i]=xyz_list[i][0];
+	PentaRef::GetInputValue(&x,x_list,gauss);
+
+	return x;
+}
+/*}}}*/
+/*FUNCTION Penta::GetYcoord {{{*/
+IssmDouble Penta::GetYcoord(Gauss* gauss){
+
+	int    i;
+	IssmDouble y;
+	IssmDouble xyz_list[NUMVERTICES][3];
+	IssmDouble y_list[NUMVERTICES];
+
+	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+	for(i=0;i<NUMVERTICES;i++) y_list[i]=xyz_list[i][1];
+	PentaRef::GetInputValue(&y,y_list,gauss);
+
+	return y;
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 17108)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 17109)
@@ -95,6 +95,7 @@
 		int    GetNumberOfVertices(void);
 		void   GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution,int enum_type);
+		IssmDouble GetXcoord(Gauss* gauss);
+		IssmDouble GetYcoord(Gauss* gauss);
 		IssmDouble GetZcoord(Gauss* gauss);
-		IssmDouble GetYcoord(Gauss* gauss){_error_("Not implemented");};
 		void   GetVectorFromInputs(Vector<IssmDouble>* vector,int name_enum);
 		void   GetVerticesCoordinates(IssmDouble** pxyz_list);
Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 17108)
+++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 17109)
@@ -139,4 +139,5 @@
 		void        GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype){_error_("not implemented yet");};
 		Node*       GetNode(int node_number){_error_("Not implemented");};
+		IssmDouble  GetXcoord(Gauss* gauss){_error_("Not implemented");};
 		IssmDouble  GetYcoord(Gauss* gauss){_error_("Not implemented");};
 		IssmDouble  GetZcoord(Gauss* gauss){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 17108)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 17109)
@@ -212,4 +212,5 @@
 
 		void	         GetVertexPidList(int* doflist);
+		IssmDouble     GetXcoord(Gauss* gauss){_error_("not implemented");};
 		IssmDouble     GetYcoord(Gauss* gauss);
 		IssmDouble     GetZcoord(Gauss* gauss){_error_("not implemented");};
Index: /issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscMat.cpp
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscMat.cpp	(revision 17108)
+++ /issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscMat.cpp	(revision 17109)
@@ -47,5 +47,5 @@
 	MatSetFromOptions(this->matrix);
 	MatMPIAIJSetPreallocation(this->matrix,0,d_nnz,0,o_nnz);
-	//MatSetOption(this->matrix,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
+	MatSetOption(this->matrix,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
 
 }