Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17140)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17141)
@@ -6,5 +6,5 @@
 #include "../solutionsequences/solutionsequences.h"
 
-//#define FSANALYTICAL 1
+//#define FSANALYTICAL 2
 
 /*Model processing*/
@@ -2762,7 +2762,15 @@
 		z_coord=element->GetZcoord(gauss);
 
-		forcex=fx1(x_coord,y_coord,z_coord);
-		forcey=fy1(x_coord,y_coord,z_coord);
-		forcez=fz1(x_coord,y_coord,z_coord);
+		#if FSANALYTICAL == 1
+			forcex=fx1(x_coord,y_coord,z_coord);
+			forcey=fy1(x_coord,y_coord,z_coord);
+			forcez=fz1(x_coord,y_coord,z_coord);
+		#elseif FSANALYTICAL == 2
+			forcex=fx2(x_coord,y_coord,z_coord);
+			forcey=fy2(x_coord,y_coord,z_coord);
+			forcez=fz2(x_coord,y_coord,z_coord);
+		#else 
+			_error_("FS analytical not implemented yet");
+		#endif
 
 		for(i=0;i<vnumnodes;i++){
@@ -3155,12 +3163,12 @@
 	 *
 	 *	In 3d
-	 *     	   Bvi=[ dh/dx          0             0      ]
-	 *					[   0           dh/dy           0      ]
-	 *					[   0             0           dh/dz    ]
-	 *					[ 1/2*dh/dy    1/2*dh/dx        0      ]
-	 *					[ 1/2*dh/dz       0         1/2*dh/dx  ]
-	 *					[   0          1/2*dh/dz    1/2*dh/dy  ]
-	 *					[   0             0             0      ]
-	 *					[ dh/dx         dh/dy         dh/dz    ]
+	 *     	   Bvi=[ dh/dx     0        0    ]
+	 *					[   0      dh/dy      0    ]
+	 *					[   0        0      dh/dz  ]
+	 *					[ dh/dy    dh/dx      0    ]
+	 *					[ dh/dz      0      dh/dx  ]
+	 *					[   0      dh/dz    dh/dy  ]
+	 *					[ dh/dx    dh/dy    dh/dz  ]
+	 *					[   0        0        0    ]
 	 *
 	 *         Bpi=[ 0 ]
@@ -3170,6 +3178,6 @@
 	 *					[ 0 ]
 	 *					[ 0 ]
+	 *					[ 0 ]
 	 *					[ h ]
-	 *					[ 0 ]
 	 *	where phi is the finiteelement function for node i.
 	 *	In 3d: