Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17358)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17359)
@@ -7,5 +7,5 @@
 #include "../cores/cores.h"
 
-//#define FSANALYTICAL 3
+//#define FSANALYTICAL 101
 
 /*Model processing*/
@@ -2273,4 +2273,63 @@
 	return Ke;
 }/*}}}*/
+#ifdef FSANALYTICAL
+ElementVector* StressbalanceAnalysis::CreatePVectorHO(Element* element){/*{{{*/
+
+	/*Intermediaries */
+	int         dim,meshtype;
+	IssmDouble  x_coord,y_coord,z_coord;
+	IssmDouble  Jdet,forcex,forcey,forcez;
+	IssmDouble* xyz_list = NULL;
+
+	/*Get problem dimension*/
+	element->FindParam(&meshtype,MeshTypeEnum);
+	switch(meshtype){
+		case Mesh2DverticalEnum: dim = 2; break;
+		case Mesh3DEnum:         dim = 3; break;
+		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+	}
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Initialize Element vector and vectors*/
+	ElementVector* pe=element->NewElementVector(HOApproximationEnum);
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinates(&xyz_list);
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=element->NewGauss(3);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		element->NodalFunctions(basis, gauss);
+
+		x_coord=element->GetXcoord(gauss);
+		y_coord=element->GetYcoord(gauss);
+		if(dim==3) z_coord=element->GetZcoord(gauss);
+		else z_coord=0.;
+
+		forcex=fx(x_coord,y_coord,z_coord,FSANALYTICAL);
+		forcey=fy(x_coord,y_coord,z_coord,FSANALYTICAL);
+
+		for(int i=0;i<numnodes;i++){
+			pe->values[i*(dim-1)+0]+=forcex*Jdet*gauss->weight*basis[i];
+			pe->values[i*(dim-1)+1]+=forcey*Jdet*gauss->weight*basis[i];
+		}
+	}
+
+	/*Transform coordinate system*/
+	if(dim==3) element->TransformLoadVectorCoord(pe,XYEnum);
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(basis);
+	xDelete<IssmDouble>(xyz_list);
+	delete gauss;
+	return pe;
+}/*}}}*/
+#else
 ElementVector* StressbalanceAnalysis::CreatePVectorHO(Element* element){/*{{{*/
 
@@ -2285,4 +2344,5 @@
 	return pe;
 }/*}}}*/
+#endif
 ElementVector* StressbalanceAnalysis::CreatePVectorHODrivingStress(Element* element){/*{{{*/
 
Index: /issm/trunk-jpl/src/c/shared/FSanalyticals/fsanalyticals.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/FSanalyticals/fsanalyticals.cpp	(revision 17358)
+++ /issm/trunk-jpl/src/c/shared/FSanalyticals/fsanalyticals.cpp	(revision 17359)
@@ -36,4 +36,7 @@
 		case 8: 
 			return fx=1.0;
+
+		case 101: 
+			return fx=fx101(x,y,z);
 		default:
 			_error_("FS analytical solution"<<testid<<" not implemented yet");
@@ -61,4 +64,7 @@
 		case 8: 
 			return fy=1.0;
+
+		case 101: 
+			return fy=fy101(x,y,z);
 		default:
 			_error_("FS analytical solution"<<testid<<" not implemented yet");
@@ -253,2 +259,13 @@
 }
 /*}}}*/
+
+IssmDouble fx101(IssmDouble x,IssmDouble y,IssmDouble z){   /*{{{*/
+
+   return 4*pow(x, 2)*y*z*pow(x - 1, 2)*(z - 1) + 8*pow(x, 2)*y*z*(y - 1)*(2*y - 1)*(z - 1) + 2*pow(x, 2)*y*pow(x - 1, 2)*(y - 1)*(2*y - 1) + 4*pow(x, 2)*z*pow(x - 1, 2)*(y - 1)*(z - 1) + 2*pow(x, 2)*z*pow(x - 1, 2)*(2*y - 1)*(z - 1) - 12*x*pow(y, 2)*z*(x - 1)*(y - 1)*(z - 1) - 6*x*pow(y, 2)*z*(2*x - 1)*(y - 1)*(z - 1) - 12*x*y*z*(x - 1)*pow(y - 1, 2)*(z - 1) + 32*x*y*z*(x - 1)*(y - 1)*(2*y - 1)*(z - 1) - 6*x*y*z*(2*x - 1)*pow(y - 1, 2)*(z - 1) - 6*pow(y, 2)*z*(x - 1)*(2*x - 1)*(y - 1)*(z - 1) + 8*y*z*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1) - 6*y*z*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1);
+}
+/*}}}*/
+IssmDouble fy101(IssmDouble x,IssmDouble y,IssmDouble z){   /*{{{*/
+
+	return 12*pow(x, 2)*y*z*(x - 1)*(y - 1)*(z - 1) + 6*pow(x, 2)*y*z*(x - 1)*(2*y - 1)*(z - 1) + 6*pow(x, 2)*z*(x - 1)*(y - 1)*(2*y - 1)*(z - 1) - 8*x*pow(y, 2)*z*(x - 1)*(2*x - 1)*(z - 1) - 4*x*pow(y, 2)*z*pow(y - 1, 2)*(z - 1) - 2*x*pow(y, 2)*(x - 1)*(2*x - 1)*pow(y - 1, 2) + 12*x*y*z*pow(x - 1, 2)*(y - 1)*(z - 1) + 6*x*y*z*pow(x - 1, 2)*(2*y - 1)*(z - 1) - 32*x*y*z*(x - 1)*(2*x - 1)*(y - 1)*(z - 1) + 6*x*z*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1) - 8*x*z*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1) - 4*pow(y, 2)*z*(x - 1)*pow(y - 1, 2)*(z - 1) - 2*pow(y, 2)*z*(2*x - 1)*pow(y - 1, 2)*(z - 1);
+}
+/*}}}*/
Index: /issm/trunk-jpl/src/c/shared/FSanalyticals/fsanalyticals.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/FSanalyticals/fsanalyticals.h	(revision 17358)
+++ /issm/trunk-jpl/src/c/shared/FSanalyticals/fsanalyticals.h	(revision 17359)
@@ -30,3 +30,5 @@
 IssmDouble fy7(IssmDouble x_coord, IssmDouble y_coord);
 
+IssmDouble fx101(IssmDouble x_coord, IssmDouble y_coord, IssmDouble z_coord);
+IssmDouble fy101(IssmDouble x_coord, IssmDouble y_coord, IssmDouble z_coord);
 #endif //ifndef _SHARED_ANALYTICALS_H_