Index: /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp	(revision 24932)
+++ /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp	(revision 24933)
@@ -2502,6 +2502,6 @@
 	for(i=0;i<numnodes;i++){
 		if(domaintype!=Domain2DverticalEnum){
-			lambdax[i]=values[i*NDOF2+0];
-			lambday[i]=values[i*NDOF2+1];
+			lambdax[i]=values[i*2+0];
+			lambday[i]=values[i*2+1];
 		}
 		else {lambdax[i]=values[i];lambday[i]=0;}
Index: /issm/trunk-jpl/src/c/analyses/BalancethicknessAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/BalancethicknessAnalysis.cpp	(revision 24932)
+++ /issm/trunk-jpl/src/c/analyses/BalancethicknessAnalysis.cpp	(revision 24933)
@@ -424,5 +424,5 @@
 }/*}}}*/
 void           BalancethicknessAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
+	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*2. 
 	 * For node i, Bi can be expressed in the actual coordinate system
 	 * by: 
@@ -431,5 +431,5 @@
 	 * where N is the finiteelement function for node i.
 	 *
-	 * We assume B_prog has been allocated already, of size: 2x(NDOF1*numnodes)
+	 * We assume B_prog has been allocated already, of size: 2x(1*numnodes)
 	 */
 
@@ -451,5 +451,5 @@
 }/*}}}*/
 void           BalancethicknessAnalysis::GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
+	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*2. 
 	 * For node i, Bi' can be expressed in the actual coordinate system
 	 * by: 
@@ -458,5 +458,5 @@
 	 * where N is the finiteelement function for node i.
 	 *
-	 * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
+	 * We assume B' has been allocated already, of size: 3x(2*numnodes)
 	 */
 
Index: /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp	(revision 24932)
+++ /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp	(revision 24933)
@@ -667,5 +667,5 @@
 }/*}}}*/
 void           DamageEvolutionAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
+	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*2.
 	 * For node i, Bi can be expressed in the actual coordinate system
 	 * by:
@@ -674,5 +674,5 @@
 	 * where N is the finiteelement function for node i.
 	 *
-	 * We assume B_prog has been allocated already, of size: 2x(NDOF1*numnodes)
+	 * We assume B_prog has been allocated already, of size: 2x(1*numnodes)
 	 */
 
@@ -695,5 +695,5 @@
 }/*}}}*/
 void           DamageEvolutionAnalysis::GetBprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2.
+	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*2.
 	 * For node i, Bi' can be expressed in the actual coordinate system
 	 * by:
@@ -702,5 +702,5 @@
 	 * where N is the finiteelement function for node i.
 	 *
-	 * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
+	 * We assume B' has been allocated already, of size: 3x(2*numnodes)
 	 */
 
Index: /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp	(revision 24932)
+++ /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp	(revision 24933)
@@ -1288,5 +1288,5 @@
 }/*}}}*/
 void           EnthalpyAnalysis::GetBAdvec(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1.
+	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*1.
 	 * For node i, Bi' can be expressed in the actual coordinate system
 	 * by:
@@ -1296,5 +1296,5 @@
 	 * where h is the interpolation function for node i.
 	 *
-	 * We assume B has been allocated already, of size: 3x(NDOF1*NUMNODESP1)
+	 * We assume B has been allocated already, of size: 3x(1*NUMNODESP1)
 	 */
 
@@ -1317,5 +1317,5 @@
 }/*}}}*/
 void           EnthalpyAnalysis::GetBAdvecprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1.
+	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*1.
 	 * For node i, Bi' can be expressed in the actual coordinate system
 	 * by:
@@ -1325,5 +1325,5 @@
 	 * where h is the interpolation function for node i.
 	 *
-	 * We assume B has been allocated already, of size: 3x(NDOF1*numnodes)
+	 * We assume B has been allocated already, of size: 3x(1*numnodes)
 	 */
 
@@ -1509,5 +1509,5 @@
 }/*}}}*/
 void           EnthalpyAnalysis::GetBConduct(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1.
+	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*1.
 	 * For node i, Bi' can be expressed in the actual coordinate system
 	 * by:
@@ -1517,5 +1517,5 @@
 	 * where h is the interpolation function for node i.
 	 *
-	 * We assume B has been allocated already, of size: 3x(NDOF1*numnodes)
+	 * We assume B has been allocated already, of size: 3x(1*numnodes)
 	 */
 
Index: /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp	(revision 24932)
+++ /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp	(revision 24933)
@@ -335,5 +335,5 @@
 }/*}}}*/
 void           FreeSurfaceBaseAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
+	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*2. 
 	 * For node i, Bi can be expressed in the actual coordinate system
 	 * by: 
@@ -342,5 +342,5 @@
 	 * where N is the finiteelement function for node i.
 	 *
-	 * We assume B_prog has been allocated already, of size: 2x(NDOF1*numnodes)
+	 * We assume B_prog has been allocated already, of size: 2x(1*numnodes)
 	 */
 
@@ -363,5 +363,5 @@
 }/*}}}*/
 void           FreeSurfaceBaseAnalysis::GetBprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
+	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*2. 
 	 * For node i, Bi' can be expressed in the actual coordinate system
 	 * by: 
@@ -370,5 +370,5 @@
 	 * where N is the finiteelement function for node i.
 	 *
-	 * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
+	 * We assume B' has been allocated already, of size: 3x(2*numnodes)
 	 */
 
Index: /issm/trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.cpp	(revision 24932)
+++ /issm/trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.cpp	(revision 24933)
@@ -307,5 +307,5 @@
 }/*}}}*/
 void           FreeSurfaceTopAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
+	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*2. 
 	 * For node i, Bi can be expressed in the actual coordinate system
 	 * by: 
@@ -314,5 +314,5 @@
 	 * where N is the finiteelement function for node i.
 	 *
-	 * We assume B_prog has been allocated already, of size: 2x(NDOF1*numnodes)
+	 * We assume B_prog has been allocated already, of size: 2x(1*numnodes)
 	 */
 
@@ -335,5 +335,5 @@
 }/*}}}*/
 void           FreeSurfaceTopAnalysis::GetBprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
+	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*2. 
 	 * For node i, Bi' can be expressed in the actual coordinate system
 	 * by: 
@@ -342,5 +342,5 @@
 	 * where N is the finiteelement function for node i.
 	 *
-	 * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
+	 * We assume B' has been allocated already, of size: 3x(2*numnodes)
 	 */
 
Index: /issm/trunk-jpl/src/c/analyses/GLheightadvectionAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/GLheightadvectionAnalysis.cpp	(revision 24932)
+++ /issm/trunk-jpl/src/c/analyses/GLheightadvectionAnalysis.cpp	(revision 24933)
@@ -246,5 +246,5 @@
 
 void           GLheightadvectionAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
+	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*2. 
 	 * For node i, Bi can be expressed in the actual coordinate system
 	 * by: 
@@ -253,5 +253,5 @@
 	 * where N is the finiteelement function for node i.
 	 *
-	 * We assume B_prog has been allocated already, of size: 2x(NDOF1*numnodes)
+	 * We assume B_prog has been allocated already, of size: 2x(1*numnodes)
 	 */
 
@@ -274,5 +274,5 @@
 }/*}}}*/
 void           GLheightadvectionAnalysis::GetBprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
+	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*2. 
 	 * For node i, Bi' can be expressed in the actual coordinate system
 	 * by: 
@@ -281,5 +281,5 @@
 	 * where N is the finiteelement function for node i.
 	 *
-	 * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
+	 * We assume B' has been allocated already, of size: 3x(2*numnodes)
 	 */
 
Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp	(revision 24932)
+++ /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp	(revision 24933)
@@ -384,5 +384,5 @@
 }/*}}}*/
 void HydrologyDCEfficientAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
+	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*2.
 	 * For node i, Bi can be expressed in the actual coordinate system
 	 * by:
@@ -391,5 +391,5 @@
 	 * where N is the finiteelement function for node i.
 	 *
-	 * We assume B has been allocated already, of size: 3x(NDOF2*numnodes)
+	 * We assume B has been allocated already, of size: 3x(2*numnodes)
 	 */
 
Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp	(revision 24932)
+++ /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp	(revision 24933)
@@ -473,5 +473,5 @@
 }/*}}}*/
 void HydrologyDCInefficientAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
+	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*2.
 	 * For node i, Bi can be expressed in the actual coordinate system
 	 * by:
@@ -480,5 +480,5 @@
 	 * where N is the finiteelement function for node i.
 	 *
-	 * We assume B has been allocated already, of size: 3x(NDOF2*numnodes)
+	 * We assume B has been allocated already, of size: 3x(2*numnodes)
 	 */
 
Index: /issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp	(revision 24932)
+++ /issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp	(revision 24933)
@@ -287,5 +287,5 @@
 }/*}}}*/
 void           HydrologyShreveAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
+	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*2. 
 	 * For node i, Bi can be expressed in the actual coordinate system
 	 * by: 
@@ -294,5 +294,5 @@
 	 * where N is the finiteelement function for node i.
 	 *
-	 * We assume B_prog has been allocated already, of size: 2x(NDOF1*numnodes)
+	 * We assume B_prog has been allocated already, of size: 2x(1*numnodes)
 	 */
 
@@ -314,5 +314,5 @@
 }/*}}}*/
 void           HydrologyShreveAnalysis::GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
+	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*2. 
 	 * For node i, Bi' can be expressed in the actual coordinate system
 	 * by: 
@@ -321,5 +321,5 @@
 	 * where N is the finiteelement function for node i.
 	 *
-	 * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
+	 * We assume B' has been allocated already, of size: 3x(2*numnodes)
 	 */
 
Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 24932)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 24933)
@@ -1815,5 +1815,5 @@
 }/*}}}*/
 void           StressbalanceAnalysis::GetBSSA(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
+	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*2.
 	 * For node i, Bi can be expressed in the actual coordinate system
 	 * by:
@@ -1824,5 +1824,5 @@
 	 * where N is the finiteelement function for node i.
 	 *
-	 * We assume B has been allocated already, of size: 3x(NDOF2*numnodes)
+	 * We assume B has been allocated already, of size: 3x(2*numnodes)
 	 */
 
@@ -1892,5 +1892,5 @@
 }/*}}}*/
 void           StressbalanceAnalysis::GetBSSAprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2.
+	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*2.
 	 * For node i, Bi' can be expressed in the actual coordinate system
 	 * by:
@@ -1901,5 +1901,5 @@
 	 * where hNis the finiteelement function for node i.
 	 *
-	 * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
+	 * We assume B' has been allocated already, of size: 3x(2*numnodes)
 	 */
 
@@ -2807,5 +2807,5 @@
 }/*}}}*/
 void           StressbalanceAnalysis::GetBHO(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2.
+	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*2.
 	 * For node i, Bi can be expressed in the actual coordinate system
 	 * by:
@@ -2819,5 +2819,5 @@
 	 * where h is the interpolation function for node i.
 	 *
-	 * We assume B has been allocated already, of size: 5x(NDOF2*numnodes)
+	 * We assume B has been allocated already, of size: 5x(2*numnodes)
 	 */
 
@@ -2892,5 +2892,5 @@
 }/*}}}*/
 void           StressbalanceAnalysis::GetBHOprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2.
+	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*2.
 	 * For node i, Bi' can be expressed in the actual coordinate system
 	 * by:
@@ -2901,5 +2901,5 @@
 	 * where hNis the finiteelement function for node i.
 	 *
-	 * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
+	 * We assume B' has been allocated already, of size: 3x(2*numnodes)
 	 */
 
@@ -3065,5 +3065,5 @@
 	int vnumnodes = element->NumberofNodesVelocity();
 	int pnumnodes = element->NumberofNodesPressure();
-	int numdof    = vnumnodes*NDOF3 + pnumnodes*NDOF1;
+	int numdof    = vnumnodes*3 + pnumnodes*1;
 
 	/*Prepare coordinate system list*/
@@ -4544,5 +4544,5 @@
 }/*}}}*/
 void           StressbalanceAnalysis::GetBFS(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B  matrix. B=[Bv1 Bv2 ... Bp1 Bp2 ...] where Bvi is of size 3*NDOF3.
+	/*Compute B  matrix. B=[Bv1 Bv2 ... Bp1 Bp2 ...] where Bvi is of size 3*3.
 	 * For node i, Bvi can be expressed in the actual coordinate system
 	 * by: 	   Bvi=[ dphi/dx          0        ]
@@ -4716,5 +4716,5 @@
 }/*}}}*/
 void           StressbalanceAnalysis::GetBFSprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*	Compute B'  matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2.
+	/*	Compute B'  matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*2.
 	 *	For node i, Bi' can be expressed in the actual coordinate system
 	 *	by:
@@ -4829,5 +4829,5 @@
 }/*}}}*/
 void           StressbalanceAnalysis::GetBFSprimeUzawa(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*	Compute B'  matrix. B'=[B1' B2' B3' B4' B5' B6'] where Bi' is of size 3*NDOF2.
+	/*	Compute B'  matrix. B'=[B1' B2' B3' B4' B5' B6'] where Bi' is of size 3*2.
 	 *	For node i, Bi' can be expressed in the actual coordinate system
 	 *	by:
@@ -4865,5 +4865,5 @@
 }/*}}}*/
 void           StressbalanceAnalysis::GetBFSprimevel(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*	Compute B'  matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2.
+	/*	Compute B'  matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*2.
 	 *	For node i, Bi' can be expressed in the actual coordinate system
 	 *	by:
@@ -4950,5 +4950,5 @@
 }/*}}}*/
 void           StressbalanceAnalysis::GetBFSvel(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B  matrix. B=[Bv1 Bv2 ... Bp1 Bp2 ...] where Bvi is of size 3*NDOF3.
+	/*Compute B  matrix. B=[Bv1 Bv2 ... Bp1 Bp2 ...] where Bvi is of size 3*3.
 	 * For node i, Bvi can be expressed in the actual coordinate system
 	 * by: 	   Bvi=[ dphi/dx          0        ]
@@ -6717,5 +6717,5 @@
 }/*}}}*/
 void           StressbalanceAnalysis::GetBprimeSSAFS(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute Bprime  matrix. Bprime=[Bprime1 Bprime2 Bprime3 Bprime4 Bprime5 Bprime6] where Bprimei is of size 5*NDOF2.
+	/*Compute Bprime  matrix. Bprime=[Bprime1 Bprime2 Bprime3 Bprime4 Bprime5 Bprime6] where Bprimei is of size 5*2.
 	 * For node i, Bprimei can be expressed in the actual coordinate system
 	 * by:
@@ -6725,5 +6725,5 @@
 	 * where h is the interpolation function for node i.
 	 *
-	 * We assume Bprime has been allocated already, of size: 5x(NDOF2*NUMNODESP1)
+	 * We assume Bprime has been allocated already, of size: 5x(2*NUMNODESP1)
 	 */
 
@@ -6766,5 +6766,5 @@
 }/*}}}*/
 void           StressbalanceAnalysis::GetBprimeSSAFSTria(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute Bprime  matrix. Bprime=[Bprime1 Bprime2 Bprime3] where Bprimei is of size 3*NDOF2.
+	/*Compute Bprime  matrix. Bprime=[Bprime1 Bprime2 Bprime3] where Bprimei is of size 3*2.
 	 * For node i, Bprimei can be expressed in the actual coordinate system
 	 * by:
@@ -6775,5 +6775,5 @@
 	 * where N is the finiteelement function for node i.
 	 *
-	 * We assume Bprime has been allocated already, of size: 3x(NDOF2*numnodes)
+	 * We assume Bprime has been allocated already, of size: 3x(2*numnodes)
 	 */
 
@@ -6801,5 +6801,5 @@
 }/*}}}*/
 void           StressbalanceAnalysis::GetBSSAFS(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2.
+	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*2.
 	 * For node i, Bi can be expressed in the actual coordinate system
 	 * by:
@@ -6810,5 +6810,5 @@
 	 * where h is the interpolation function for node i.
 	 *
-	 * We assume B has been allocated already, of size: 5x(NDOF2*NUMNODESP1)
+	 * We assume B has been allocated already, of size: 5x(2*NUMNODESP1)
 	 */
 
@@ -6859,5 +6859,5 @@
 }/*}}}*/
 void           StressbalanceAnalysis::GetBSSAFSTria(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
+	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*2.
 	 * For node i, Bi can be expressed in the actual coordinate system
 	 * by:
@@ -6867,5 +6867,5 @@
 	 * where N is the finiteelement function for node i.
 	 *
-	 * We assume B has been allocated already, of size: 3x(NDOF2*numnodes)
+	 * We assume B has been allocated already, of size: 3x(2*numnodes)
 	 */
 
@@ -6891,5 +6891,5 @@
 }/*}}}*/
 void           StressbalanceAnalysis::GetBSSAHO(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*NDOF2.
+	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*2.
 	 * For node i, Bi can be expressed in the actual coordinate system
 	 * by:
@@ -6899,5 +6899,5 @@
 	 * where h is the interpolation function for node i.
 	 *
-	 * We assume B has been allocated already, of size: 5x(NDOF2*NUMNODESP1)
+	 * We assume B has been allocated already, of size: 5x(2*NUMNODESP1)
 	 */
 
Index: /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp	(revision 24932)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp	(revision 24933)
@@ -589,6 +589,6 @@
 	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
 	for(i=0;i<numnodes;i++){
-		vx[i]=values[i*NDOF2+0];
-		vy[i]=values[i*NDOF2+1];
+		vx[i]=values[i*2+0];
+		vy[i]=values[i*2+1];
 
 		/*Check solution*/
Index: /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp	(revision 24932)
+++ /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp	(revision 24933)
@@ -746,5 +746,5 @@
 }/*}}}*/
 void           ThermalAnalysis::GetBAdvec(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1.
+	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*1.
 	 * For node i, Bi' can be expressed in the actual coordinate system
 	 * by:
@@ -754,5 +754,5 @@
 	 * where h is the interpolation function for node i.
 	 *
-	 * We assume B has been allocated already, of size: 3x(NDOF1*NUMNODESP1)
+	 * We assume B has been allocated already, of size: 3x(1*NUMNODESP1)
 	 */
 
@@ -775,5 +775,5 @@
 }/*}}}*/
 void           ThermalAnalysis::GetBAdvecprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1.
+	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*1.
 	 * For node i, Bi' can be expressed in the actual coordinate system
 	 * by:
@@ -783,5 +783,5 @@
 	 * where h is the interpolation function for node i.
 	 *
-	 * We assume B has been allocated already, of size: 3x(NDOF1*numnodes)
+	 * We assume B has been allocated already, of size: 3x(1*numnodes)
 	 */
 
@@ -804,5 +804,5 @@
 }/*}}}*/
 void           ThermalAnalysis::GetBConduct(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1.
+	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*1.
 	 * For node i, Bi' can be expressed in the actual coordinate system
 	 * by:
@@ -812,5 +812,5 @@
 	 * where h is the interpolation function for node i.
 	 *
-	 * We assume B has been allocated already, of size: 3x(NDOF1*numnodes)
+	 * We assume B has been allocated already, of size: 3x(1*numnodes)
 	 */
 
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 24932)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 24933)
@@ -2686,6 +2686,6 @@
 void       Penta::InputUpdateFromSolutionOneDofCollapsed(IssmDouble* solution,int enum_type){/*{{{*/
 
-	const int  numdof   = NDOF1*NUMVERTICES;
-	const int  numdof2d = NDOF1*NUMVERTICES2D;
+	const int  numdof   = NUMVERTICES;
+	const int  numdof2d = NUMVERTICES2D;
 
 	IssmDouble  values[numdof];
@@ -2726,5 +2726,5 @@
 void       Penta::InputUpdateFromVector(IssmDouble* vector, int name, int type){/*{{{*/
 
-	const int   numdof         = NDOF1 *NUMVERTICES;
+	const int   numdof         = NUMVERTICES;
 	int        *doflist        = NULL;
 	IssmDouble  values[numdof];
Index: /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 24932)
+++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 24933)
@@ -191,5 +191,5 @@
 void PentaRef::GetJacobian(IssmDouble* J, IssmDouble* xyz_list,Gauss* gauss_in){/*{{{*/
 	/*The Jacobian is constant over the element, discard the gaussian points. 
-	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
+	 * J is assumed to have been allocated of size 2x2.*/
 
 	IssmDouble A1,A2,A3;  // area coordinates
@@ -235,15 +235,15 @@
 	j_const_reciprocal=SQRT3/12;
 
-	J[NDOF3*0+0] = 0.25*(x1-x2-x4+x5)*zi+0.25*(-x1+x2-x4+x5);
-	J[NDOF3*1+0] = j_const_reciprocal*(x1+x2-2*x3-x4-x5+2*x6)*zi+j_const_reciprocal*(-x1-x2+2*x3-x4-x5+2*x6);
-	J[NDOF3*2+0] = j_const_reciprocal*(x1+x2-2*x3-x4-x5+2*x6)*eta+0.25*(x1-x2-x4+x5)*xi +0.25*(-x1+x5-x2+x4);
-
-	J[NDOF3*0+1] = 0.25*(y1-y2-y4+y5)*zi+0.25*(-y1+y2-y4+y5);
-	J[NDOF3*1+1] = j_const_reciprocal*(y1+y2-2*y3-y4-y5+2*y6)*zi+j_const_reciprocal*(-y1-y2+2*y3-y4-y5+2*y6);
-	J[NDOF3*2+1] = j_const_reciprocal*(y1+y2-2*y3-y4-y5+2*y6)*eta+0.25*(y1-y2-y4+y5)*xi+0.25*(y4-y1+y5-y2);
-
-	J[NDOF3*0+2] = 0.25*(z1-z2-z4+z5)*zi+0.25*(-z1+z2-z4+z5);
-	J[NDOF3*1+2] = j_const_reciprocal*(z1+z2-2*z3-z4-z5+2*z6)*zi+j_const_reciprocal*(-z1-z2+2*z3-z4-z5+2*z6);
-	J[NDOF3*2+2] = j_const_reciprocal*(z1+z2-2*z3-z4-z5+2*z6)*eta+0.25*(z1-z2-z4+z5)*xi+0.25*(-z1+z5-z2+z4);
+	J[3*0+0] = 0.25*(x1-x2-x4+x5)*zi+0.25*(-x1+x2-x4+x5);
+	J[3*1+0] = j_const_reciprocal*(x1+x2-2*x3-x4-x5+2*x6)*zi+j_const_reciprocal*(-x1-x2+2*x3-x4-x5+2*x6);
+	J[3*2+0] = j_const_reciprocal*(x1+x2-2*x3-x4-x5+2*x6)*eta+0.25*(x1-x2-x4+x5)*xi +0.25*(-x1+x5-x2+x4);
+
+	J[3*0+1] = 0.25*(y1-y2-y4+y5)*zi+0.25*(-y1+y2-y4+y5);
+	J[3*1+1] = j_const_reciprocal*(y1+y2-2*y3-y4-y5+2*y6)*zi+j_const_reciprocal*(-y1-y2+2*y3-y4-y5+2*y6);
+	J[3*2+1] = j_const_reciprocal*(y1+y2-2*y3-y4-y5+2*y6)*eta+0.25*(y1-y2-y4+y5)*xi+0.25*(y4-y1+y5-y2);
+
+	J[3*0+2] = 0.25*(z1-z2-z4+z5)*zi+0.25*(-z1+z2-z4+z5);
+	J[3*1+2] = j_const_reciprocal*(z1+z2-2*z3-z4-z5+2*z6)*zi+j_const_reciprocal*(-z1-z2+2*z3-z4-z5+2*z6);
+	J[3*2+2] = j_const_reciprocal*(z1+z2-2*z3-z4-z5+2*z6)*eta+0.25*(z1-z2-z4+z5)*xi+0.25*(-z1+z5-z2+z4);
 }
 /*}}}*/
@@ -1090,5 +1090,5 @@
 void PentaRef::GetSegmentJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
 	/*The Jacobian determinant is constant over the element, discard the gaussian points. 
-	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
+	 * J is assumed to have been allocated of size 2x2.*/
 
 	IssmDouble x1=xyz_list[3*0+0];
@@ -1106,5 +1106,5 @@
 void PentaRef::GetTriaJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
 	/*The Jacobian determinant is constant over the element, discard the gaussian points. 
-	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
+	 * J is assumed to have been allocated of size 2x2.*/
 
 	IssmDouble x1=xyz_list[3*0+0];
Index: /issm/trunk-jpl/src/c/classes/Elements/SegRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/SegRef.cpp	(revision 24932)
+++ /issm/trunk-jpl/src/c/classes/Elements/SegRef.cpp	(revision 24933)
@@ -93,5 +93,5 @@
 void SegRef::GetJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,GaussSeg* gauss){/*{{{*/
 	/*The Jacobian determinant is constant over the element, discard the gaussian points. 
-	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
+	 * J is assumed to have been allocated of size 1.*/
 
 	/*Call Jacobian routine to get the jacobian:*/
Index: /issm/trunk-jpl/src/c/classes/Elements/TetraRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/TetraRef.cpp	(revision 24932)
+++ /issm/trunk-jpl/src/c/classes/Elements/TetraRef.cpp	(revision 24933)
@@ -110,20 +110,20 @@
 	IssmDouble z4=xyz_list[3*3+2];
 
-	J[NDOF3*0+0] = x2-x1;
-	J[NDOF3*0+1] = y2-y1;
-	J[NDOF3*0+2] = z2-z1;
-
-	J[NDOF3*1+0] = x3-x1;
-	J[NDOF3*1+1] = y3-y1;
-	J[NDOF3*1+2] = z3-z1;
-
-	J[NDOF3*2+0] = x4-x1;
-	J[NDOF3*2+1] = y4-y1;
-	J[NDOF3*2+2] = z4-z1;
+	J[3*0+0] = x2-x1;
+	J[3*0+1] = y2-y1;
+	J[3*0+2] = z2-z1;
+
+	J[3*1+0] = x3-x1;
+	J[3*1+1] = y3-y1;
+	J[3*1+2] = z3-z1;
+
+	J[3*2+0] = x4-x1;
+	J[3*2+1] = y4-y1;
+	J[3*2+2] = z4-z1;
 }
 /*}}}*/
 void TetraRef::GetJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,GaussTetra* gauss){/*{{{*/
 	/*The Jacobian determinant is constant over the element, discard the gaussian points. 
-	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
+	 * J is assumed to have been allocated of size 2x2.*/
 	IssmDouble J[3][3];
 
@@ -139,5 +139,5 @@
 void TetraRef::GetJacobianDeterminantFace(IssmDouble*  Jdet, IssmDouble* xyz_list,GaussTetra* gauss){/*{{{*/
 	/*The Jacobian determinant is constant over the element, discard the gaussian points. 
-	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
+	 * J is assumed to have been allocated of size 2x2.*/
 
 	IssmDouble x1=xyz_list[3*0+0];
Index: /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp	(revision 24932)
+++ /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp	(revision 24933)
@@ -90,5 +90,5 @@
 void TriaRef::GetJacobian(IssmDouble* J, IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
 	/*The Jacobian is constant over the element, discard the gaussian points.
-	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
+	 * J is assumed to have been allocated of size 2x2.*/
 
 	IssmDouble x1 = xyz_list[3*0+0];
@@ -107,5 +107,5 @@
 void TriaRef::GetJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
 	/*The Jacobian determinant is constant over the element, discard the gaussian points.
-	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
+	 * J is assumed to have been allocated of size 2x2.*/
 	IssmDouble J[2][2];
 
@@ -121,5 +121,5 @@
 void TriaRef::GetJacobianDeterminant3D(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
 	/*The Jacobian determinant is constant over the element, discard the gaussian points.
-	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
+	 * J is assumed to have been allocated of size 2x2.*/
 	IssmDouble J[2][2];
 
Index: /issm/trunk-jpl/src/c/classes/Loads/Penpair.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Penpair.cpp	(revision 24932)
+++ /issm/trunk-jpl/src/c/classes/Loads/Penpair.cpp	(revision 24933)
@@ -290,5 +290,5 @@
 ElementMatrix* Penpair::PenaltyCreateKMatrixMasstransport(IssmDouble kmax){/*{{{*/
 
-	const int numdof=NUMVERTICES*NDOF1;
+	const int numdof=NUMVERTICES*1;
 	IssmDouble penalty_factor;
 
Index: /issm/trunk-jpl/src/c/classes/Loads/Riftfront.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Riftfront.cpp	(revision 24932)
+++ /issm/trunk-jpl/src/c/classes/Loads/Riftfront.cpp	(revision 24933)
@@ -431,5 +431,5 @@
 ElementMatrix* Riftfront::PenaltyCreateKMatrixStressbalanceHoriz(IssmDouble kmax){/*{{{*/
 
-	const int   numdof = NDOF2*NUMVERTICES;
+	const int   numdof = 2*NUMVERTICES;
 	IssmDouble  thickness;
 	IssmDouble  h[2];
Index: /issm/trunk-jpl/src/c/shared/Numerics/constants.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Numerics/constants.h	(revision 24932)
+++ /issm/trunk-jpl/src/c/shared/Numerics/constants.h	(revision 24933)
@@ -11,25 +11,3 @@
 const double PI=3.141592653589793238462643383279502884197169399375105820974944592308; // Macro definition conflicts with Dakota's declaration of PI
 
-#define NDOF1 1
-#define NDOF2 2
-#define NDOF3 3
-#define NDOF4 4
-
-// /*Windows specific typefefs: */
-// #ifdef _INTEL_WIN_
-//
-// #ifndef NAN
-// //For reference, for Intel compile on win64
-// //#define NAN 0.0/0.0
-// #define NAN (INFINITY-INFINITY)
-// #endif
-//
-// #ifndef INFINITY
-// //For reference, for Intel compile on win64
-// //#define INFINITY 1.0/0.0
-// #define INFINITY (DBL_MAX+DBL_MAX)
-// #endif
-//
-// #endif /*_INTEL_WIN_*/
-
 #endif /*_ISSM_CONSTANTS_H_*/