Index: /issm/trunk/src/c/objects/Elements/Element.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Element.h	(revision 5744)
+++ /issm/trunk/src/c/objects/Elements/Element.h	(revision 5745)
@@ -32,6 +32,6 @@
 		virtual void   GetSolutionFromInputs(Vec solution)=0;
 		virtual int    GetNodeIndex(Node* node)=0;
-		virtual bool   GetShelf()=0; 
-		virtual bool   GetOnBed()=0;
+		virtual bool   IsOnShelf()=0; 
+		virtual bool   IsOnBed()=0;
 		virtual void   GetParameterListOnVertices(double* pvalue,int enumtype)=0;
 		virtual void   GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue)=0;
Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 5744)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 5745)
@@ -897,24 +897,4 @@
 	ISSMERROR("Node provided not found among element nodes");
 
-}
-/*}}}*/
-/*FUNCTION Penta::GetOnBed {{{1*/
-bool Penta::GetOnBed(){
-	
-	bool onbed;
-
-	inputs->GetParameterValue(&onbed,ElementOnBedEnum);
-
-	return onbed;
-}
-/*}}}*/
-/*FUNCTION Penta::GetShelf {{{1*/
-bool   Penta::GetShelf(){
-	bool onshelf;
-
-	/*retrieve inputs :*/
-	inputs->GetParameterValue(&onshelf,ElementOnIceShelfEnum);
-
-	return onshelf;
 }
 /*}}}*/
@@ -2001,5 +1981,5 @@
 
 	/*If shelf: hydrostatic equilibrium*/
-	if (this->GetShelf()){
+	if (this->IsOnShelf()){
 
 		/*recover material parameters: */
@@ -5483,4 +5463,13 @@
 }
 /*}}}*/
+/*FUNCTION Penta::IsOnShelf {{{1*/
+bool   Penta::IsOnShelf(){
+
+	bool onshelf;
+	inputs->GetParameterValue(&onshelf,ElementOnIceShelfEnum);
+	return onshelf;
+
+}
+/*}}}*/
 /*FUNCTION Penta::IsOnSurface{{{1*/
 bool Penta::IsOnSurface(void){
Index: /issm/trunk/src/c/objects/Elements/Penta.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.h	(revision 5744)
+++ /issm/trunk/src/c/objects/Elements/Penta.h	(revision 5745)
@@ -79,6 +79,4 @@
 		void   DeleteResults(void);
 		int    GetNodeIndex(Node* node);
-		bool   GetOnBed();
-		bool   GetShelf(); 
 		void   GetSolutionFromInputs(Vec solution);
 		void   GetVectorFromInputs(Vec vector,int NameEnum);
@@ -179,4 +177,5 @@
 		bool	  IsOnSurface(void);
 		bool	  IsOnBed(void);
+		bool    IsOnShelf(void); 
 		void	  ReduceMatrixStokes(double* Ke_reduced, double* Ke_temp);
 		void	  ReduceVectorStokes(double* Pe_reduced, double* Ke_temp, double* Pe_temp);
Index: /issm/trunk/src/c/objects/Elements/PentaRef.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/PentaRef.cpp	(revision 5744)
+++ /issm/trunk/src/c/objects/Elements/PentaRef.cpp	(revision 5745)
@@ -19,4 +19,13 @@
 #include "../../include/include.h"
 /*}}}*/
+
+/*Element macros*/
+#define NUMNODESP1    6
+#define NUMNODESP1_2d 3
+#define NUMNODESMINI  7
+#define NDOF1 1
+#define NDOF2 2
+#define NDOF3 3
+#define NDOF4 4
 
 /*Object constructors and destructor*/
@@ -62,13 +71,8 @@
 	 * where h is the interpolation function for grid i.
 	 *
-	 * We assume B has been allocated already, of size: 5x(NDOF2*numgrids)
-	 */
-
-	int i;
-	const int numgrids=6;
-	const int NDOF3=3;
-	const int NDOF2=2;
-
-	double dh1dh6[NDOF3][numgrids];
+	 * We assume B has been allocated already, of size: 5x(NDOF2*NUMNODESP1)
+	 */
+
+	double dh1dh6[3][NUMNODESP1];
 
 	/*Get dh1dh6 in actual coordinate system: */
@@ -76,13 +80,13 @@
 
 	/*Build B: */
-	for (i=0;i<numgrids;i++){
-		*(B+NDOF2*numgrids*0+NDOF2*i)=dh1dh6[0][i]; 
-		*(B+NDOF2*numgrids*0+NDOF2*i+1)=0.0;
-
-		*(B+NDOF2*numgrids*1+NDOF2*i)=0.0;
-		*(B+NDOF2*numgrids*1+NDOF2*i+1)=dh1dh6[1][i];
-
-		*(B+NDOF2*numgrids*2+NDOF2*i)=(float).5*dh1dh6[1][i]; 
-		*(B+NDOF2*numgrids*2+NDOF2*i+1)=(float).5*dh1dh6[0][i]; 
+	for (int i=0;i<NUMNODESP1;i++){
+		*(B+NDOF2*NUMNODESP1*0+NDOF2*i)=dh1dh6[0][i]; 
+		*(B+NDOF2*NUMNODESP1*0+NDOF2*i+1)=0.0;
+
+		*(B+NDOF2*NUMNODESP1*1+NDOF2*i)=0.0;
+		*(B+NDOF2*NUMNODESP1*1+NDOF2*i+1)=dh1dh6[1][i];
+
+		*(B+NDOF2*NUMNODESP1*2+NDOF2*i)=(float).5*dh1dh6[1][i]; 
+		*(B+NDOF2*NUMNODESP1*2+NDOF2*i+1)=(float).5*dh1dh6[0][i]; 
 
 	}
@@ -102,13 +106,8 @@
 	 * where h is the interpolation function for grid i.
 	 *
-	 * We assume B has been allocated already, of size: 5x(NDOF2*numgrids)
-	 */
-
-	int i;
-	const int numgrids=6;
-	const int NDOF3=3;
-	const int NDOF2=2;
-
-	double dh1dh6[NDOF3][numgrids];
+	 * We assume B has been allocated already, of size: 5x(NDOF2*NUMNODESP1)
+	 */
+
+	double dh1dh6[3][NUMNODESP1];
 
 	/*Get dh1dh6 in actual coordinate system: */
@@ -116,19 +115,19 @@
 
 	/*Build B: */
-	for (i=0;i<numgrids;i++){
-		*(B+NDOF2*numgrids*0+NDOF2*i)=dh1dh6[0][i]; 
-		*(B+NDOF2*numgrids*0+NDOF2*i+1)=0.0;
-
-		*(B+NDOF2*numgrids*1+NDOF2*i)=0.0;
-		*(B+NDOF2*numgrids*1+NDOF2*i+1)=dh1dh6[1][i];
-
-		*(B+NDOF2*numgrids*2+NDOF2*i)=(float).5*dh1dh6[1][i]; 
-		*(B+NDOF2*numgrids*2+NDOF2*i+1)=(float).5*dh1dh6[0][i]; 
-
-		*(B+NDOF2*numgrids*3+NDOF2*i)=(float).5*dh1dh6[2][i]; 
-		*(B+NDOF2*numgrids*3+NDOF2*i+1)=0.0;
-
-		*(B+NDOF2*numgrids*4+NDOF2*i)=0.0;
-		*(B+NDOF2*numgrids*4+NDOF2*i+1)=(float).5*dh1dh6[2][i]; 
+	for (int i=0;i<NUMNODESP1;i++){
+		*(B+NDOF2*NUMNODESP1*0+NDOF2*i)=dh1dh6[0][i]; 
+		*(B+NDOF2*NUMNODESP1*0+NDOF2*i+1)=0.0;
+
+		*(B+NDOF2*NUMNODESP1*1+NDOF2*i)=0.0;
+		*(B+NDOF2*NUMNODESP1*1+NDOF2*i+1)=dh1dh6[1][i];
+
+		*(B+NDOF2*NUMNODESP1*2+NDOF2*i)=(float).5*dh1dh6[1][i]; 
+		*(B+NDOF2*NUMNODESP1*2+NDOF2*i+1)=(float).5*dh1dh6[0][i]; 
+
+		*(B+NDOF2*NUMNODESP1*3+NDOF2*i)=(float).5*dh1dh6[2][i]; 
+		*(B+NDOF2*NUMNODESP1*3+NDOF2*i+1)=0.0;
+
+		*(B+NDOF2*NUMNODESP1*4+NDOF2*i)=0.0;
+		*(B+NDOF2*NUMNODESP1*4+NDOF2*i+1)=(float).5*dh1dh6[2][i]; 
 	}
 
@@ -147,13 +146,7 @@
 	 * where h is the interpolation function for grid i.
 	 *
-	 * We assume B has been allocated already, of size: 5x(NDOF2*numgrids)
-	 */
-
-	int i;
-	const int NDOF3=3;
-	const int NDOF2=2;
-	const int numgrids=6;
-
-	double dh1dh6[NDOF3][numgrids];
+	 * We assume B has been allocated already, of size: 5x(NDOF2*NUMNODESP1)
+	 */
+	double dh1dh6[3][NUMNODESP1];
 
 	/*Get dh1dh6 in actual coordinate system: */
@@ -161,19 +154,19 @@
 
 	/*Build BPrime: */
-	for (i=0;i<numgrids;i++){
-		*(B+NDOF2*numgrids*0+NDOF2*i)=2.0*dh1dh6[0][i]; 
-		*(B+NDOF2*numgrids*0+NDOF2*i+1)=dh1dh6[1][i];
-
-		*(B+NDOF2*numgrids*1+NDOF2*i)=dh1dh6[0][i];
-		*(B+NDOF2*numgrids*1+NDOF2*i+1)=2.0*dh1dh6[1][i];
-
-		*(B+NDOF2*numgrids*2+NDOF2*i)=dh1dh6[1][i]; 
-		*(B+NDOF2*numgrids*2+NDOF2*i+1)=dh1dh6[0][i]; 
-
-		*(B+NDOF2*numgrids*3+NDOF2*i)=dh1dh6[2][i]; 
-		*(B+NDOF2*numgrids*3+NDOF2*i+1)=0.0;
-
-		*(B+NDOF2*numgrids*4+NDOF2*i)=0.0;
-		*(B+NDOF2*numgrids*4+NDOF2*i+1)=dh1dh6[2][i]; 
+	for (int i=0;i<NUMNODESP1;i++){
+		*(B+NDOF2*NUMNODESP1*0+NDOF2*i)=2.0*dh1dh6[0][i]; 
+		*(B+NDOF2*NUMNODESP1*0+NDOF2*i+1)=dh1dh6[1][i];
+
+		*(B+NDOF2*NUMNODESP1*1+NDOF2*i)=dh1dh6[0][i];
+		*(B+NDOF2*NUMNODESP1*1+NDOF2*i+1)=2.0*dh1dh6[1][i];
+
+		*(B+NDOF2*NUMNODESP1*2+NDOF2*i)=dh1dh6[1][i]; 
+		*(B+NDOF2*NUMNODESP1*2+NDOF2*i+1)=dh1dh6[0][i]; 
+
+		*(B+NDOF2*NUMNODESP1*3+NDOF2*i)=dh1dh6[2][i]; 
+		*(B+NDOF2*NUMNODESP1*3+NDOF2*i+1)=0.0;
+
+		*(B+NDOF2*NUMNODESP1*4+NDOF2*i)=0.0;
+		*(B+NDOF2*NUMNODESP1*4+NDOF2*i+1)=dh1dh6[2][i]; 
 	}
 }
@@ -182,5 +175,5 @@
 void PentaRef::GetBStokes(double* B, double* xyz_list, GaussPenta* gauss){
 
-	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*DOFPERGRID. 
+	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*NDOF4. 
 	 * For grid i, Bi can be expressed in the actual coordinate system
 	 * by: 		Bi=[ dh/dx          0              0       0  ]
@@ -197,10 +190,7 @@
 
 	int i;
-	const int calculationdof=3;
-	const int numgrids=6;
-	int DOFPERGRID=4;
-
-	double dh1dh7[calculationdof][numgrids+1];
-	double l1l6[numgrids];
+
+	double dh1dh7[3][NUMNODESMINI];
+	double l1l6[NUMNODESP1];
 
 
@@ -210,40 +200,40 @@
 
 	/*Build B: */
-	for (i=0;i<numgrids+1;i++){
-		*(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i)=dh1dh7[0][i]; //B[0][DOFPERGRID*i]=dh1dh6[0][i];
-		*(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+1)=0;
-		*(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+2)=0;
-		*(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i)=0;
-		*(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+1)=dh1dh7[1][i];
-		*(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+2)=0;
-		*(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i)=0;
-		*(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+1)=0;
-		*(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+2)=dh1dh7[2][i];
-		*(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i)=(float).5*dh1dh7[1][i]; 
-		*(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+1)=(float).5*dh1dh7[0][i]; 
-		*(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+2)=0;
-		*(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i)=(float).5*dh1dh7[2][i];
-		*(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+1)=0;
-		*(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+2)=(float).5*dh1dh7[0][i];
-		*(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i)=0;
-		*(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+1)=(float).5*dh1dh7[2][i];
-		*(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+2)=(float).5*dh1dh7[1][i];
-		*(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i)=0;
-		*(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+1)=0;
-		*(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+2)=0;
-		*(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i)=dh1dh7[0][i];
-		*(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+1)=dh1dh7[1][i];
-		*(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+2)=dh1dh7[2][i];
-	}
-
-	for (i=0;i<numgrids;i++){ //last column not for the bubble function
-		*(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+3)=0;
-		*(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+3)=0;
-		*(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+3)=0;
-		*(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+3)=0;
-		*(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+3)=0;
-		*(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+3)=0;
-		*(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+3)=l1l6[i];
-		*(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+3)=0;
+	for (i=0;i<NUMNODESMINI;i++){
+		*(B+(NDOF4*NUMNODESP1+3)*0+NDOF4*i)=dh1dh7[0][i]; //B[0][NDOF4*i]=dh1dh6[0][i];
+		*(B+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+1)=0;
+		*(B+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+2)=0;
+		*(B+(NDOF4*NUMNODESP1+3)*1+NDOF4*i)=0;
+		*(B+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+1)=dh1dh7[1][i];
+		*(B+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+2)=0;
+		*(B+(NDOF4*NUMNODESP1+3)*2+NDOF4*i)=0;
+		*(B+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+1)=0;
+		*(B+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2)=dh1dh7[2][i];
+		*(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i)=(float).5*dh1dh7[1][i]; 
+		*(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+1)=(float).5*dh1dh7[0][i]; 
+		*(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+2)=0;
+		*(B+(NDOF4*NUMNODESP1+3)*4+NDOF4*i)=(float).5*dh1dh7[2][i];
+		*(B+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+1)=0;
+		*(B+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+2)=(float).5*dh1dh7[0][i];
+		*(B+(NDOF4*NUMNODESP1+3)*5+NDOF4*i)=0;
+		*(B+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+1)=(float).5*dh1dh7[2][i];
+		*(B+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+2)=(float).5*dh1dh7[1][i];
+		*(B+(NDOF4*NUMNODESP1+3)*6+NDOF4*i)=0;
+		*(B+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+1)=0;
+		*(B+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+2)=0;
+		*(B+(NDOF4*NUMNODESP1+3)*7+NDOF4*i)=dh1dh7[0][i];
+		*(B+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+1)=dh1dh7[1][i];
+		*(B+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+2)=dh1dh7[2][i];
+	}
+
+	for (i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
+		*(B+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+3)=0;
+		*(B+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+3)=0;
+		*(B+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+3)=0;
+		*(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+3)=0;
+		*(B+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+3)=0;
+		*(B+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+3)=0;
+		*(B+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+3)=l1l6[i];
+		*(B+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+3)=0;
 	}
 
@@ -269,53 +259,48 @@
 
 	int i;
-	const int calculationdof=3;
-	const int numgrids=6;
-	int DOFPERGRID=4;
-
-	double dh1dh7[calculationdof][numgrids+1];
-	double l1l6[numgrids];
+	double dh1dh7[3][NUMNODESMINI];
+	double l1l6[NUMNODESP1];
 
 	/*Get dh1dh7 in actual coordinate system: */
 	GetNodalFunctionsMINIDerivatives(&dh1dh7[0][0],xyz_list, gauss);
-
 	GetNodalFunctionsP1(l1l6, gauss);
 
 	/*B_primeuild B_prime: */
-	for (i=0;i<numgrids+1;i++){
-		*(B_prime+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i)=dh1dh7[0][i]; //B_prime[0][DOFPERGRID*i]=dh1dh6[0][i];
-		*(B_prime+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+1)=0;
-		*(B_prime+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+2)=0;
-		*(B_prime+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i)=0;
-		*(B_prime+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+1)=dh1dh7[1][i];
-		*(B_prime+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+2)=0;
-		*(B_prime+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i)=0;
-		*(B_prime+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+1)=0;
-		*(B_prime+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+2)=dh1dh7[2][i];
-		*(B_prime+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i)=dh1dh7[1][i]; 
-		*(B_prime+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+1)=dh1dh7[0][i]; 
-		*(B_prime+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+2)=0;
-		*(B_prime+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i)=dh1dh7[2][i];
-		*(B_prime+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+1)=0;
-		*(B_prime+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+2)=dh1dh7[0][i];
-		*(B_prime+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i)=0;
-		*(B_prime+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+1)=dh1dh7[2][i];
-		*(B_prime+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+2)=dh1dh7[1][i];
-		*(B_prime+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i)=dh1dh7[0][i];
-		*(B_prime+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+1)=dh1dh7[1][i];
-		*(B_prime+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+2)=dh1dh7[2][i];
-		*(B_prime+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i)=0;
-		*(B_prime+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+1)=0;
-		*(B_prime+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+2)=0;
-	}
-
-	for (i=0;i<numgrids;i++){ //last column not for the bubble function
-		*(B_prime+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+3)=0;
-		*(B_prime+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+3)=0;
-		*(B_prime+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+3)=0;
-		*(B_prime+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+3)=0;
-		*(B_prime+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+3)=0;
-		*(B_prime+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+3)=0;
-		*(B_prime+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+3)=0;
-		*(B_prime+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+3)=l1l6[i];
+	for (i=0;i<NUMNODESMINI;i++){
+		*(B_prime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i)=dh1dh7[0][i]; //B_prime[0][NDOF4*i]=dh1dh6[0][i];
+		*(B_prime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+1)=0;
+		*(B_prime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+2)=0;
+		*(B_prime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i)=0;
+		*(B_prime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+1)=dh1dh7[1][i];
+		*(B_prime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+2)=0;
+		*(B_prime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i)=0;
+		*(B_prime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+1)=0;
+		*(B_prime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2)=dh1dh7[2][i];
+		*(B_prime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i)=dh1dh7[1][i]; 
+		*(B_prime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+1)=dh1dh7[0][i]; 
+		*(B_prime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+2)=0;
+		*(B_prime+(NDOF4*NUMNODESP1+3)*4+NDOF4*i)=dh1dh7[2][i];
+		*(B_prime+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+1)=0;
+		*(B_prime+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+2)=dh1dh7[0][i];
+		*(B_prime+(NDOF4*NUMNODESP1+3)*5+NDOF4*i)=0;
+		*(B_prime+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+1)=dh1dh7[2][i];
+		*(B_prime+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+2)=dh1dh7[1][i];
+		*(B_prime+(NDOF4*NUMNODESP1+3)*6+NDOF4*i)=dh1dh7[0][i];
+		*(B_prime+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+1)=dh1dh7[1][i];
+		*(B_prime+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+2)=dh1dh7[2][i];
+		*(B_prime+(NDOF4*NUMNODESP1+3)*7+NDOF4*i)=0;
+		*(B_prime+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+1)=0;
+		*(B_prime+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+2)=0;
+	}
+
+	for (i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
+		*(B_prime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+3)=0;
+		*(B_prime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+3)=0;
+		*(B_prime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+3)=0;
+		*(B_prime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+3)=0;
+		*(B_prime+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+3)=0;
+		*(B_prime+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+3)=0;
+		*(B_prime+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+3)=0;
+		*(B_prime+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+3)=l1l6[i];
 	}
 
@@ -324,5 +309,5 @@
 /*FUNCTION PentaRef::GetBArtdiff {{{1*/
 void PentaRef::GetBArtdiff(double* B_artdiff, double* xyz_list, GaussPenta* gauss){
-	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID. 
+	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1. 
 	 * For grid i, Bi' can be expressed in the actual coordinate system
 	 * by: 
@@ -331,14 +316,9 @@
 	 * where h is the interpolation function for grid i.
 	 *
-	 * We assume B has been allocated already, of size: 2x(DOFPERGRID*numgrids)
-	 */
-
-	int i;
-	const int calculationdof=3;
-	const int numgrids=6;
-	int DOFPERGRID=1;
+	 * We assume B has been allocated already, of size: 2x(NDOF1*NUMNODESP1)
+	 */
 
 	/*Same thing in the actual coordinate system: */
-	double dh1dh6[calculationdof][numgrids];
+	double dh1dh6[3][NUMNODESP1];
 
 	/*Get dh1dh6 in actual coordinates system : */
@@ -346,7 +326,7 @@
 
 	/*Build B': */
-	for (i=0;i<numgrids;i++){
-		*(B_artdiff+DOFPERGRID*numgrids*0+DOFPERGRID*i)=dh1dh6[0][i]; 
-		*(B_artdiff+DOFPERGRID*numgrids*1+DOFPERGRID*i)=dh1dh6[1][i]; 
+	for (int i=0;i<NUMNODESP1;i++){
+		*(B_artdiff+NDOF1*NUMNODESP1*0+NDOF1*i)=dh1dh6[0][i]; 
+		*(B_artdiff+NDOF1*NUMNODESP1*1+NDOF1*i)=dh1dh6[1][i]; 
 	}
 }
@@ -354,5 +334,5 @@
 /*FUNCTION PentaRef::GetBAdvec{{{1*/
 void PentaRef::GetBAdvec(double* B_advec, double* xyz_list, GaussPenta* gauss){
-	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID. 
+	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1. 
 	 * For grid i, Bi' can be expressed in the actual coordinate system
 	 * by: 
@@ -362,11 +342,6 @@
 	 * where h is the interpolation function for grid i.
 	 *
-	 * We assume B has been allocated already, of size: 3x(DOFPERGRID*numgrids)
-	 */
-
-	int i;
-	int calculationdof=3;
-	int numgrids=6;
-	int DOFPERGRID=1;
+	 * We assume B has been allocated already, of size: 3x(NDOF1*NUMNODESP1)
+	 */
 
 	/*Same thing in the actual coordinate system: */
@@ -377,8 +352,8 @@
 
 	/*Build B': */
-	for (i=0;i<numgrids;i++){
-		*(B_advec+DOFPERGRID*numgrids*0+DOFPERGRID*i)=l1l6[i]; 
-		*(B_advec+DOFPERGRID*numgrids*1+DOFPERGRID*i)=l1l6[i]; 
-		*(B_advec+DOFPERGRID*numgrids*2+DOFPERGRID*i)=l1l6[i]; 
+	for (int i=0;i<NUMNODESP1;i++){
+		*(B_advec+NDOF1*NUMNODESP1*0+NDOF1*i)=l1l6[i]; 
+		*(B_advec+NDOF1*NUMNODESP1*1+NDOF1*i)=l1l6[i]; 
+		*(B_advec+NDOF1*NUMNODESP1*2+NDOF1*i)=l1l6[i]; 
 	}
 }
@@ -386,5 +361,5 @@
 /*FUNCTION PentaRef::GetBConduct{{{1*/
 void PentaRef::GetBConduct(double* B_conduct, double* xyz_list, GaussPenta* gauss){
-	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID. 
+	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1. 
 	 * For grid i, Bi' can be expressed in the actual coordinate system
 	 * by: 
@@ -394,14 +369,9 @@
 	 * where h is the interpolation function for grid i.
 	 *
-	 * We assume B has been allocated already, of size: 3x(DOFPERGRID*numgrids)
-	 */
-
-	int i;
-	const int calculationdof=3;
-	const int numgrids=6;
-	int DOFPERGRID=1;
+	 * We assume B has been allocated already, of size: 3x(NDOF1*NUMNODESP1)
+	 */
 
 	/*Same thing in the actual coordinate system: */
-	double dh1dh6[calculationdof][numgrids];
+	double dh1dh6[3][NUMNODESP1];
 
 	/*Get dh1dh2dh3 in actual coordinates system : */
@@ -409,8 +379,8 @@
 
 	/*Build B': */
-	for (i=0;i<numgrids;i++){
-		*(B_conduct+DOFPERGRID*numgrids*0+DOFPERGRID*i)=dh1dh6[0][i]; 
-		*(B_conduct+DOFPERGRID*numgrids*1+DOFPERGRID*i)=dh1dh6[1][i]; 
-		*(B_conduct+DOFPERGRID*numgrids*2+DOFPERGRID*i)=dh1dh6[2][i]; 
+	for (int i=0;i<NUMNODESP1;i++){
+		*(B_conduct+NDOF1*NUMNODESP1*0+NDOF1*i)=dh1dh6[0][i]; 
+		*(B_conduct+NDOF1*NUMNODESP1*1+NDOF1*i)=dh1dh6[1][i]; 
+		*(B_conduct+NDOF1*NUMNODESP1*2+NDOF1*i)=dh1dh6[2][i]; 
 	}
 }
@@ -422,7 +392,5 @@
 
 	int i;
-	const int NDOF3=3;
-	const int numgrids=6;
-	double dh1dh6[NDOF3][numgrids];
+	double dh1dh6[3][NUMNODESP1];
 
 	/*Get dh1dh6 in actual coordinate system: */
@@ -430,5 +398,5 @@
 
 	/*Build B: */
-	for (i=0;i<numgrids;i++){
+	for (i=0;i<NUMNODESP1;i++){
 		B[i]=dh1dh6[2][i];  
 	}
@@ -438,5 +406,5 @@
 /*FUNCTION PentaRef::GetBprimeAdvec{{{1*/
 void PentaRef::GetBprimeAdvec(double* Bprime_advec, double* xyz_list, GaussPenta* gauss){
-	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID. 
+	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1. 
 	 * For grid i, Bi' can be expressed in the actual coordinate system
 	 * by: 
@@ -446,14 +414,9 @@
 	 * where h is the interpolation function for grid i.
 	 *
-	 * We assume B has been allocated already, of size: 3x(DOFPERGRID*numgrids)
-	 */
-
-	int i;
-	const int calculationdof=3;
-	const int numgrids=6;
-	int DOFPERGRID=1;
+	 * We assume B has been allocated already, of size: 3x(NDOF1*NUMNODESP1)
+	 */
 
 	/*Same thing in the actual coordinate system: */
-	double dh1dh6[calculationdof][numgrids];
+	double dh1dh6[3][NUMNODESP1];
 
 	/*Get dh1dh2dh3 in actual coordinates system : */
@@ -461,8 +424,8 @@
 
 	/*Build B': */
-	for (i=0;i<numgrids;i++){
-		*(Bprime_advec+DOFPERGRID*numgrids*0+DOFPERGRID*i)=dh1dh6[0][i]; 
-		*(Bprime_advec+DOFPERGRID*numgrids*1+DOFPERGRID*i)=dh1dh6[1][i]; 
-		*(Bprime_advec+DOFPERGRID*numgrids*2+DOFPERGRID*i)=dh1dh6[2][i]; 
+	for (int i=0;i<NUMNODESP1;i++){
+		*(Bprime_advec+NDOF1*NUMNODESP1*0+NDOF1*i)=dh1dh6[0][i]; 
+		*(Bprime_advec+NDOF1*NUMNODESP1*1+NDOF1*i)=dh1dh6[1][i]; 
+		*(Bprime_advec+NDOF1*NUMNODESP1*2+NDOF1*i)=dh1dh6[2][i]; 
 	}
 }
@@ -500,8 +463,7 @@
 
 	int i;
-	const int numgrids2d=3;
 	int num_dof=4;
 
-	double l1l2l3[numgrids2d];
+	double l1l2l3[NUMNODESP1_2d];
 
 
@@ -513,60 +475,60 @@
 	/*Build LStokes: */
 	for (i=0;i<3;i++){
-		*(LStokes+num_dof*numgrids2d*0+num_dof*i)=l1l2l3[i]; //LStokes[0][NDOF2*i]=dh1dh3[0][i];
-		*(LStokes+num_dof*numgrids2d*0+num_dof*i+1)=0;
-		*(LStokes+num_dof*numgrids2d*0+num_dof*i+2)=0;
-		*(LStokes+num_dof*numgrids2d*0+num_dof*i+3)=0;
-		*(LStokes+num_dof*numgrids2d*1+num_dof*i)=0;
-		*(LStokes+num_dof*numgrids2d*1+num_dof*i+1)=l1l2l3[i];
-		*(LStokes+num_dof*numgrids2d*1+num_dof*i+2)=0;
-		*(LStokes+num_dof*numgrids2d*1+num_dof*i+3)=0;
-		*(LStokes+num_dof*numgrids2d*2+num_dof*i)=0;
-		*(LStokes+num_dof*numgrids2d*2+num_dof*i+1)=0;
-		*(LStokes+num_dof*numgrids2d*2+num_dof*i+2)=l1l2l3[i];
-		*(LStokes+num_dof*numgrids2d*2+num_dof*i+3)=0;
-		*(LStokes+num_dof*numgrids2d*3+num_dof*i)=0;
-		*(LStokes+num_dof*numgrids2d*3+num_dof*i+1)=0;
-		*(LStokes+num_dof*numgrids2d*3+num_dof*i+2)=l1l2l3[i];
-		*(LStokes+num_dof*numgrids2d*3+num_dof*i+3)=0;
-		*(LStokes+num_dof*numgrids2d*4+num_dof*i)=l1l2l3[i];
-		*(LStokes+num_dof*numgrids2d*4+num_dof*i+1)=0;
-		*(LStokes+num_dof*numgrids2d*4+num_dof*i+2)=0;
-		*(LStokes+num_dof*numgrids2d*4+num_dof*i+3)=0;
-		*(LStokes+num_dof*numgrids2d*5+num_dof*i)=0;
-		*(LStokes+num_dof*numgrids2d*5+num_dof*i+1)=l1l2l3[i];
-		*(LStokes+num_dof*numgrids2d*5+num_dof*i+2)=0;
-		*(LStokes+num_dof*numgrids2d*5+num_dof*i+3)=0;
-		*(LStokes+num_dof*numgrids2d*6+num_dof*i)=l1l2l3[i];
-		*(LStokes+num_dof*numgrids2d*6+num_dof*i+1)=0;
-		*(LStokes+num_dof*numgrids2d*6+num_dof*i+2)=0;
-		*(LStokes+num_dof*numgrids2d*6+num_dof*i+3)=0;
-		*(LStokes+num_dof*numgrids2d*7+num_dof*i)=0;
-		*(LStokes+num_dof*numgrids2d*7+num_dof*i+1)=l1l2l3[i];
-		*(LStokes+num_dof*numgrids2d*7+num_dof*i+2)=0;
-		*(LStokes+num_dof*numgrids2d*7+num_dof*i+3)=0;
-		*(LStokes+num_dof*numgrids2d*8+num_dof*i)=0;
-		*(LStokes+num_dof*numgrids2d*8+num_dof*i+1)=0;
-		*(LStokes+num_dof*numgrids2d*8+num_dof*i+2)=l1l2l3[i];
-		*(LStokes+num_dof*numgrids2d*8+num_dof*i+3)=0;
-		*(LStokes+num_dof*numgrids2d*9+num_dof*i)=0;
-		*(LStokes+num_dof*numgrids2d*9+num_dof*i+1)=0;
-		*(LStokes+num_dof*numgrids2d*9+num_dof*i+2)=l1l2l3[i];
-		*(LStokes+num_dof*numgrids2d*9+num_dof*i+3)=0;
-		*(LStokes+num_dof*numgrids2d*10+num_dof*i)=0;
-		*(LStokes+num_dof*numgrids2d*10+num_dof*i+1)=0;
-		*(LStokes+num_dof*numgrids2d*10+num_dof*i+2)=l1l2l3[i];
-		*(LStokes+num_dof*numgrids2d*10+num_dof*i+3)=0;
-		*(LStokes+num_dof*numgrids2d*11+num_dof*i)=l1l2l3[i];
-		*(LStokes+num_dof*numgrids2d*11+num_dof*i+1)=0;
-		*(LStokes+num_dof*numgrids2d*11+num_dof*i+2)=0;
-		*(LStokes+num_dof*numgrids2d*11+num_dof*i+3)=0;
-		*(LStokes+num_dof*numgrids2d*12+num_dof*i)=0;
-		*(LStokes+num_dof*numgrids2d*12+num_dof*i+1)=l1l2l3[i];
-		*(LStokes+num_dof*numgrids2d*12+num_dof*i+2)=0;
-		*(LStokes+num_dof*numgrids2d*12+num_dof*i+3)=0;
-		*(LStokes+num_dof*numgrids2d*13+num_dof*i)=0;
-		*(LStokes+num_dof*numgrids2d*13+num_dof*i+1)=0;
-		*(LStokes+num_dof*numgrids2d*13+num_dof*i+2)=l1l2l3[i];
-		*(LStokes+num_dof*numgrids2d*13+num_dof*i+3)=0;
+		*(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i)=l1l2l3[i]; //LStokes[0][NDOF2*i]=dh1dh3[0][i];
+		*(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0;
+		*(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+2)=0;
+		*(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+3)=0;
+		*(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i)=0;
+		*(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i];
+		*(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+2)=0;
+		*(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+3)=0;
+		*(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i)=0;
+		*(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+1)=0;
+		*(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+2)=l1l2l3[i];
+		*(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+3)=0;
+		*(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i)=0;
+		*(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+1)=0;
+		*(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+2)=l1l2l3[i];
+		*(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+3)=0;
+		*(LStokes+num_dof*NUMNODESP1_2d*4+num_dof*i)=l1l2l3[i];
+		*(LStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+1)=0;
+		*(LStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+2)=0;
+		*(LStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+3)=0;
+		*(LStokes+num_dof*NUMNODESP1_2d*5+num_dof*i)=0;
+		*(LStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+1)=l1l2l3[i];
+		*(LStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+2)=0;
+		*(LStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+3)=0;
+		*(LStokes+num_dof*NUMNODESP1_2d*6+num_dof*i)=l1l2l3[i];
+		*(LStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+1)=0;
+		*(LStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+2)=0;
+		*(LStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+3)=0;
+		*(LStokes+num_dof*NUMNODESP1_2d*7+num_dof*i)=0;
+		*(LStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+1)=l1l2l3[i];
+		*(LStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+2)=0;
+		*(LStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+3)=0;
+		*(LStokes+num_dof*NUMNODESP1_2d*8+num_dof*i)=0;
+		*(LStokes+num_dof*NUMNODESP1_2d*8+num_dof*i+1)=0;
+		*(LStokes+num_dof*NUMNODESP1_2d*8+num_dof*i+2)=l1l2l3[i];
+		*(LStokes+num_dof*NUMNODESP1_2d*8+num_dof*i+3)=0;
+		*(LStokes+num_dof*NUMNODESP1_2d*9+num_dof*i)=0;
+		*(LStokes+num_dof*NUMNODESP1_2d*9+num_dof*i+1)=0;
+		*(LStokes+num_dof*NUMNODESP1_2d*9+num_dof*i+2)=l1l2l3[i];
+		*(LStokes+num_dof*NUMNODESP1_2d*9+num_dof*i+3)=0;
+		*(LStokes+num_dof*NUMNODESP1_2d*10+num_dof*i)=0;
+		*(LStokes+num_dof*NUMNODESP1_2d*10+num_dof*i+1)=0;
+		*(LStokes+num_dof*NUMNODESP1_2d*10+num_dof*i+2)=l1l2l3[i];
+		*(LStokes+num_dof*NUMNODESP1_2d*10+num_dof*i+3)=0;
+		*(LStokes+num_dof*NUMNODESP1_2d*11+num_dof*i)=l1l2l3[i];
+		*(LStokes+num_dof*NUMNODESP1_2d*11+num_dof*i+1)=0;
+		*(LStokes+num_dof*NUMNODESP1_2d*11+num_dof*i+2)=0;
+		*(LStokes+num_dof*NUMNODESP1_2d*11+num_dof*i+3)=0;
+		*(LStokes+num_dof*NUMNODESP1_2d*12+num_dof*i)=0;
+		*(LStokes+num_dof*NUMNODESP1_2d*12+num_dof*i+1)=l1l2l3[i];
+		*(LStokes+num_dof*NUMNODESP1_2d*12+num_dof*i+2)=0;
+		*(LStokes+num_dof*NUMNODESP1_2d*12+num_dof*i+3)=0;
+		*(LStokes+num_dof*NUMNODESP1_2d*13+num_dof*i)=0;
+		*(LStokes+num_dof*NUMNODESP1_2d*13+num_dof*i+1)=0;
+		*(LStokes+num_dof*NUMNODESP1_2d*13+num_dof*i+2)=l1l2l3[i];
+		*(LStokes+num_dof*NUMNODESP1_2d*13+num_dof*i+3)=0;
 
 	}
@@ -597,9 +559,8 @@
 	 */
 	int i;
-	const int numgrids2d=3;
 	int num_dof=4;
 
-	double l1l2l3[numgrids2d];
-	double dh1dh6[3][6];
+	double l1l2l3[NUMNODESP1_2d];
+	double dh1dh6[3][NUMNODESP1];
 
 	/*Get l1l2l3 in actual coordinate system: */
@@ -612,61 +573,60 @@
 	/*Build LprimeStokes: */
 	for (i=0;i<3;i++){
-		*(LprimeStokes+num_dof*numgrids2d*0+num_dof*i)=l1l2l3[i]; //LprimeStokes[0][NDOF2*i]=dh1dh3[0][i];
-		*(LprimeStokes+num_dof*numgrids2d*0+num_dof*i+1)=0;
-		*(LprimeStokes+num_dof*numgrids2d*0+num_dof*i+2)=0;
-		*(LprimeStokes+num_dof*numgrids2d*0+num_dof*i+3)=0;
-		*(LprimeStokes+num_dof*numgrids2d*1+num_dof*i)=0;
-		*(LprimeStokes+num_dof*numgrids2d*1+num_dof*i+1)=l1l2l3[i];
-		*(LprimeStokes+num_dof*numgrids2d*1+num_dof*i+2)=0;
-		*(LprimeStokes+num_dof*numgrids2d*1+num_dof*i+3)=0;
-		*(LprimeStokes+num_dof*numgrids2d*2+num_dof*i)=l1l2l3[i];
-		*(LprimeStokes+num_dof*numgrids2d*2+num_dof*i+1)=0;
-		*(LprimeStokes+num_dof*numgrids2d*2+num_dof*i+2)=0;
-		*(LprimeStokes+num_dof*numgrids2d*2+num_dof*i+3)=0;
-		*(LprimeStokes+num_dof*numgrids2d*3+num_dof*i)=0;
-		*(LprimeStokes+num_dof*numgrids2d*3+num_dof*i+1)=l1l2l3[i];
-		*(LprimeStokes+num_dof*numgrids2d*3+num_dof*i+2)=0;
-		*(LprimeStokes+num_dof*numgrids2d*3+num_dof*i+3)=0;
-		*(LprimeStokes+num_dof*numgrids2d*4+num_dof*i)=0;
-		*(LprimeStokes+num_dof*numgrids2d*4+num_dof*i+1)=0;
-		*(LprimeStokes+num_dof*numgrids2d*4+num_dof*i+2)=l1l2l3[i];
-		*(LprimeStokes+num_dof*numgrids2d*4+num_dof*i+3)=0;
-		*(LprimeStokes+num_dof*numgrids2d*5+num_dof*i)=0;
-		*(LprimeStokes+num_dof*numgrids2d*5+num_dof*i+1)=0;
-		*(LprimeStokes+num_dof*numgrids2d*5+num_dof*i+2)=l1l2l3[i];
-		*(LprimeStokes+num_dof*numgrids2d*5+num_dof*i+3)=0;
-		*(LprimeStokes+num_dof*numgrids2d*6+num_dof*i)=0;
-		*(LprimeStokes+num_dof*numgrids2d*6+num_dof*i+1)=0;
-		*(LprimeStokes+num_dof*numgrids2d*6+num_dof*i+2)=dh1dh6[2][i];
-		*(LprimeStokes+num_dof*numgrids2d*6+num_dof*i+3)=0;
-		*(LprimeStokes+num_dof*numgrids2d*7+num_dof*i)=0;
-		*(LprimeStokes+num_dof*numgrids2d*7+num_dof*i+1)=0;
-		*(LprimeStokes+num_dof*numgrids2d*7+num_dof*i+2)=dh1dh6[2][i];
-		*(LprimeStokes+num_dof*numgrids2d*7+num_dof*i+3)=0;
-		*(LprimeStokes+num_dof*numgrids2d*8+num_dof*i)=0;
-		*(LprimeStokes+num_dof*numgrids2d*8+num_dof*i+1)=0;
-		*(LprimeStokes+num_dof*numgrids2d*8+num_dof*i+2)=dh1dh6[2][i];
-		*(LprimeStokes+num_dof*numgrids2d*8+num_dof*i+3)=0;
-		*(LprimeStokes+num_dof*numgrids2d*9+num_dof*i)=dh1dh6[2][i];
-		*(LprimeStokes+num_dof*numgrids2d*9+num_dof*i+1)=0;
-		*(LprimeStokes+num_dof*numgrids2d*9+num_dof*i+2)=dh1dh6[0][i];
-		*(LprimeStokes+num_dof*numgrids2d*9+num_dof*i+3)=0;
-		*(LprimeStokes+num_dof*numgrids2d*10+num_dof*i)=0;
-		*(LprimeStokes+num_dof*numgrids2d*10+num_dof*i+1)=dh1dh6[2][i];
-		*(LprimeStokes+num_dof*numgrids2d*10+num_dof*i+2)=dh1dh6[1][i];
-		*(LprimeStokes+num_dof*numgrids2d*10+num_dof*i+3)=0;
-		*(LprimeStokes+num_dof*numgrids2d*11+num_dof*i)=0;
-		*(LprimeStokes+num_dof*numgrids2d*11+num_dof*i+1)=0;
-		*(LprimeStokes+num_dof*numgrids2d*11+num_dof*i+2)=0;
-		*(LprimeStokes+num_dof*numgrids2d*11+num_dof*i+3)=l1l2l3[i];
-		*(LprimeStokes+num_dof*numgrids2d*12+num_dof*i)=0;
-		*(LprimeStokes+num_dof*numgrids2d*12+num_dof*i+1)=0;
-		*(LprimeStokes+num_dof*numgrids2d*12+num_dof*i+2)=0;
-		*(LprimeStokes+num_dof*numgrids2d*12+num_dof*i+3)=l1l2l3[i];
-		*(LprimeStokes+num_dof*numgrids2d*13+num_dof*i)=0;
-		*(LprimeStokes+num_dof*numgrids2d*13+num_dof*i+1)=0;
-		*(LprimeStokes+num_dof*numgrids2d*13+num_dof*i+2)=0;
-		*(LprimeStokes+num_dof*numgrids2d*13+num_dof*i+3)=l1l2l3[i];
-
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i)=l1l2l3[i]; //LprimeStokes[0][NDOF2*i]=dh1dh3[0][i];
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0;
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+2)=0;
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+3)=0;
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i)=0;
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i];
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+2)=0;
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+3)=0;
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i)=l1l2l3[i];
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+1)=0;
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+2)=0;
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+3)=0;
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i)=0;
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+1)=l1l2l3[i];
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+2)=0;
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+3)=0;
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i)=0;
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+1)=0;
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+2)=l1l2l3[i];
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+3)=0;
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i)=0;
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+1)=0;
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+2)=l1l2l3[i];
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+3)=0;
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i)=0;
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+1)=0;
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+2)=dh1dh6[2][i];
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+3)=0;
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i)=0;
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+1)=0;
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+2)=dh1dh6[2][i];
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+3)=0;
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*8+num_dof*i)=0;
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*8+num_dof*i+1)=0;
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*8+num_dof*i+2)=dh1dh6[2][i];
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*8+num_dof*i+3)=0;
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*9+num_dof*i)=dh1dh6[2][i];
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*9+num_dof*i+1)=0;
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*9+num_dof*i+2)=dh1dh6[0][i];
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*9+num_dof*i+3)=0;
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*10+num_dof*i)=0;
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*10+num_dof*i+1)=dh1dh6[2][i];
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*10+num_dof*i+2)=dh1dh6[1][i];
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*10+num_dof*i+3)=0;
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*11+num_dof*i)=0;
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*11+num_dof*i+1)=0;
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*11+num_dof*i+2)=0;
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*11+num_dof*i+3)=l1l2l3[i];
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*12+num_dof*i)=0;
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*12+num_dof*i+1)=0;
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*12+num_dof*i+2)=0;
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*12+num_dof*i+3)=l1l2l3[i];
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*13+num_dof*i)=0;
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*13+num_dof*i+1)=0;
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*13+num_dof*i+2)=0;
+		*(LprimeStokes+num_dof*NUMNODESP1_2d*13+num_dof*i+3)=l1l2l3[i];
 	}
 }
@@ -675,5 +635,4 @@
 void PentaRef::GetJacobian(double* J, double* xyz_list,GaussPenta* gauss){
 
-	const int NDOF3=3;
 	int i,j;
 
@@ -822,6 +781,5 @@
 
 	int       i;
-	const int numgrids = 7;
-	double    dh1dh7_ref[3][numgrids];
+	double    dh1dh7_ref[3][NUMNODESMINI];
 	double    Jinv[3][3];
 
@@ -839,8 +797,8 @@
 	 */
 
-	for (i=0;i<numgrids;i++){
-		*(dh1dh7+numgrids*0+i)=Jinv[0][0]*dh1dh7_ref[0][i]+Jinv[0][1]*dh1dh7_ref[1][i]+Jinv[0][2]*dh1dh7_ref[2][i];
-		*(dh1dh7+numgrids*1+i)=Jinv[1][0]*dh1dh7_ref[0][i]+Jinv[1][1]*dh1dh7_ref[1][i]+Jinv[1][2]*dh1dh7_ref[2][i];
-		*(dh1dh7+numgrids*2+i)=Jinv[2][0]*dh1dh7_ref[0][i]+Jinv[2][1]*dh1dh7_ref[1][i]+Jinv[2][2]*dh1dh7_ref[2][i];
+	for (i=0;i<NUMNODESMINI;i++){
+		*(dh1dh7+NUMNODESMINI*0+i)=Jinv[0][0]*dh1dh7_ref[0][i]+Jinv[0][1]*dh1dh7_ref[1][i]+Jinv[0][2]*dh1dh7_ref[2][i];
+		*(dh1dh7+NUMNODESMINI*1+i)=Jinv[1][0]*dh1dh7_ref[0][i]+Jinv[1][1]*dh1dh7_ref[1][i]+Jinv[1][2]*dh1dh7_ref[2][i];
+		*(dh1dh7+NUMNODESMINI*2+i)=Jinv[2][0]*dh1dh7_ref[0][i]+Jinv[2][1]*dh1dh7_ref[1][i]+Jinv[2][2]*dh1dh7_ref[2][i];
 	}
 
@@ -852,7 +810,4 @@
 	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
 	 * natural coordinate system) at the gaussian point. */
-
-	int    numgrids=7; //six plus bubble grids
-
 	double r=gauss->coord2-gauss->coord1;
 	double s=-3.0/SQRT3*(gauss->coord1+gauss->coord2-2.0/3.0);
@@ -860,37 +815,37 @@
 
 	/*First nodal function: */
-	*(dl1dl7+numgrids*0+0)=-0.5*(1.0-zeta)/2.0;
-	*(dl1dl7+numgrids*1+0)=-SQRT3/6.0*(1.0-zeta)/2.0;
-	*(dl1dl7+numgrids*2+0)=-0.5*(-0.5*r-SQRT3/6.0*s+ONETHIRD);
+	*(dl1dl7+NUMNODESMINI*0+0)=-0.5*(1.0-zeta)/2.0;
+	*(dl1dl7+NUMNODESMINI*1+0)=-SQRT3/6.0*(1.0-zeta)/2.0;
+	*(dl1dl7+NUMNODESMINI*2+0)=-0.5*(-0.5*r-SQRT3/6.0*s+ONETHIRD);
 
 	/*Second nodal function: */
-	*(dl1dl7+numgrids*0+1)=0.5*(1.0-zeta)/2.0;
-	*(dl1dl7+numgrids*1+1)=-SQRT3/6.0*(1.0-zeta)/2.0;
-	*(dl1dl7+numgrids*2+1)=-0.5*(0.5*r-SQRT3/6.0*s+ONETHIRD);
+	*(dl1dl7+NUMNODESMINI*0+1)=0.5*(1.0-zeta)/2.0;
+	*(dl1dl7+NUMNODESMINI*1+1)=-SQRT3/6.0*(1.0-zeta)/2.0;
+	*(dl1dl7+NUMNODESMINI*2+1)=-0.5*(0.5*r-SQRT3/6.0*s+ONETHIRD);
 
 	/*Third nodal function: */
-	*(dl1dl7+numgrids*0+2)=0;
-	*(dl1dl7+numgrids*1+2)=SQRT3/3.0*(1.0-zeta)/2.0;
-	*(dl1dl7+numgrids*2+2)=-0.5*(SQRT3/3.0*s+ONETHIRD);
+	*(dl1dl7+NUMNODESMINI*0+2)=0;
+	*(dl1dl7+NUMNODESMINI*1+2)=SQRT3/3.0*(1.0-zeta)/2.0;
+	*(dl1dl7+NUMNODESMINI*2+2)=-0.5*(SQRT3/3.0*s+ONETHIRD);
 
 	/*Fourth nodal function: */
-	*(dl1dl7+numgrids*0+3)=-0.5*(1.0+zeta)/2.0;
-	*(dl1dl7+numgrids*1+3)=-SQRT3/6.0*(1.0+zeta)/2.0;
-	*(dl1dl7+numgrids*2+3)=0.5*(-0.5*r-SQRT3/6.0*s+ONETHIRD);
+	*(dl1dl7+NUMNODESMINI*0+3)=-0.5*(1.0+zeta)/2.0;
+	*(dl1dl7+NUMNODESMINI*1+3)=-SQRT3/6.0*(1.0+zeta)/2.0;
+	*(dl1dl7+NUMNODESMINI*2+3)=0.5*(-0.5*r-SQRT3/6.0*s+ONETHIRD);
 
 	/*Fith nodal function: */
-	*(dl1dl7+numgrids*0+4)=0.5*(1.0+zeta)/2.0;
-	*(dl1dl7+numgrids*1+4)=-SQRT3/6.0*(1.0+zeta)/2.0;
-	*(dl1dl7+numgrids*2+4)=0.5*(0.5*r-SQRT3/6.0*s+ONETHIRD);
+	*(dl1dl7+NUMNODESMINI*0+4)=0.5*(1.0+zeta)/2.0;
+	*(dl1dl7+NUMNODESMINI*1+4)=-SQRT3/6.0*(1.0+zeta)/2.0;
+	*(dl1dl7+NUMNODESMINI*2+4)=0.5*(0.5*r-SQRT3/6.0*s+ONETHIRD);
 
 	/*Sixth nodal function: */
-	*(dl1dl7+numgrids*0+5)=0;
-	*(dl1dl7+numgrids*1+5)=SQRT3/3.0*(1.0+zeta)/2.0;
-	*(dl1dl7+numgrids*2+5)=0.5*(SQRT3/3.0*s+ONETHIRD);
+	*(dl1dl7+NUMNODESMINI*0+5)=0;
+	*(dl1dl7+NUMNODESMINI*1+5)=SQRT3/3.0*(1.0+zeta)/2.0;
+	*(dl1dl7+NUMNODESMINI*2+5)=0.5*(SQRT3/3.0*s+ONETHIRD);
 
 	/*Seventh nodal function: */
-	*(dl1dl7+numgrids*0+6)=9.0/2.0*r*(1.0+zeta)*(zeta-1.0)*(SQRT3*s+1.0);
-	*(dl1dl7+numgrids*1+6)=9.0/4.0*(1+zeta)*(1-zeta)*(SQRT3*pow(s,2.0)-2.0*s-SQRT3*pow(r,2.0));
-	*(dl1dl7+numgrids*2+6)=27*gauss->coord1*gauss->coord2*gauss->coord3*(-2.0*zeta);
+	*(dl1dl7+NUMNODESMINI*0+6)=9.0/2.0*r*(1.0+zeta)*(zeta-1.0)*(SQRT3*s+1.0);
+	*(dl1dl7+NUMNODESMINI*1+6)=9.0/4.0*(1+zeta)*(1-zeta)*(SQRT3*pow(s,2.0)-2.0*s-SQRT3*pow(r,2.0));
+	*(dl1dl7+NUMNODESMINI*2+6)=27*gauss->coord1*gauss->coord2*gauss->coord3*(-2.0*zeta);
 
 }
@@ -914,8 +869,5 @@
 	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
 	 * actual coordinate system): */
-	int       i;
-	const int NDOF3    = 3;
-	const int numgrids = 6;
-	double    dh1dh6_ref[NDOF3][numgrids];
+	double    dh1dh6_ref[NDOF3][NUMNODESP1];
 	double    Jinv[NDOF3][NDOF3];
 
@@ -933,8 +885,8 @@
 	 */
 
-	for (i=0;i<numgrids;i++){
-		*(dh1dh6+numgrids*0+i)=Jinv[0][0]*dh1dh6_ref[0][i]+Jinv[0][1]*dh1dh6_ref[1][i]+Jinv[0][2]*dh1dh6_ref[2][i];
-		*(dh1dh6+numgrids*1+i)=Jinv[1][0]*dh1dh6_ref[0][i]+Jinv[1][1]*dh1dh6_ref[1][i]+Jinv[1][2]*dh1dh6_ref[2][i];
-		*(dh1dh6+numgrids*2+i)=Jinv[2][0]*dh1dh6_ref[0][i]+Jinv[2][1]*dh1dh6_ref[1][i]+Jinv[2][2]*dh1dh6_ref[2][i];
+	for (int i=0;i<NUMNODESP1;i++){
+		*(dh1dh6+NUMNODESP1*0+i)=Jinv[0][0]*dh1dh6_ref[0][i]+Jinv[0][1]*dh1dh6_ref[1][i]+Jinv[0][2]*dh1dh6_ref[2][i];
+		*(dh1dh6+NUMNODESP1*1+i)=Jinv[1][0]*dh1dh6_ref[0][i]+Jinv[1][1]*dh1dh6_ref[1][i]+Jinv[1][2]*dh1dh6_ref[2][i];
+		*(dh1dh6+NUMNODESP1*2+i)=Jinv[2][0]*dh1dh6_ref[0][i]+Jinv[2][1]*dh1dh6_ref[1][i]+Jinv[2][2]*dh1dh6_ref[2][i];
 	}
 
@@ -947,5 +899,4 @@
 	 * natural coordinate system) at the gaussian point. Those values vary along xi,eta,z */
 
-	const int numgrids=6;
 	double A1,A2,A3,z;
 
@@ -956,32 +907,32 @@
 
 	/*First nodal function derivatives. The corresponding nodal function is N=A1*(1-z)/2. Its derivatives follow*/
-	*(dl1dl6+numgrids*0+0)=-0.5*(1.0-z)/2.0;
-	*(dl1dl6+numgrids*1+0)=-0.5/SQRT3*(1.0-z)/2.0;
-	*(dl1dl6+numgrids*2+0)=-0.5*A1;
+	*(dl1dl6+NUMNODESP1*0+0)=-0.5*(1.0-z)/2.0;
+	*(dl1dl6+NUMNODESP1*1+0)=-0.5/SQRT3*(1.0-z)/2.0;
+	*(dl1dl6+NUMNODESP1*2+0)=-0.5*A1;
 
 	/*Second nodal function: The corresponding nodal function is N=A2*(1-z)/2. Its derivatives follow*/
-	*(dl1dl6+numgrids*0+1)=0.5*(1.0-z)/2.0;
-	*(dl1dl6+numgrids*1+1)=-0.5/SQRT3*(1.0-z)/2.0;
-	*(dl1dl6+numgrids*2+1)=-0.5*A2;
+	*(dl1dl6+NUMNODESP1*0+1)=0.5*(1.0-z)/2.0;
+	*(dl1dl6+NUMNODESP1*1+1)=-0.5/SQRT3*(1.0-z)/2.0;
+	*(dl1dl6+NUMNODESP1*2+1)=-0.5*A2;
 
 	/*Third nodal function: The corresponding nodal function is N=A3*(1-z)/2. Its derivatives follow*/
-	*(dl1dl6+numgrids*0+2)=0.0;
-	*(dl1dl6+numgrids*1+2)=1.0/SQRT3*(1.0-z)/2.0;
-	*(dl1dl6+numgrids*2+2)=-0.5*A3;
+	*(dl1dl6+NUMNODESP1*0+2)=0.0;
+	*(dl1dl6+NUMNODESP1*1+2)=1.0/SQRT3*(1.0-z)/2.0;
+	*(dl1dl6+NUMNODESP1*2+2)=-0.5*A3;
 
 	/*Fourth nodal function: The corresponding nodal function is N=A1*(1+z)/2. Its derivatives follow*/
-	*(dl1dl6+numgrids*0+3)=-0.5*(1.0+z)/2.0;
-	*(dl1dl6+numgrids*1+3)=-0.5/SQRT3*(1.0+z)/2.0;
-	*(dl1dl6+numgrids*2+3)=0.5*A1;
+	*(dl1dl6+NUMNODESP1*0+3)=-0.5*(1.0+z)/2.0;
+	*(dl1dl6+NUMNODESP1*1+3)=-0.5/SQRT3*(1.0+z)/2.0;
+	*(dl1dl6+NUMNODESP1*2+3)=0.5*A1;
 
 	/*Fifth nodal function: The corresponding nodal function is N=A2*(1+z)/2. Its derivatives follow*/
-	*(dl1dl6+numgrids*0+4)=0.5*(1.0+z)/2.0;
-	*(dl1dl6+numgrids*1+4)=-0.5/SQRT3*(1.0+z)/2.0;
-	*(dl1dl6+numgrids*2+4)=0.5*A2;
+	*(dl1dl6+NUMNODESP1*0+4)=0.5*(1.0+z)/2.0;
+	*(dl1dl6+NUMNODESP1*1+4)=-0.5/SQRT3*(1.0+z)/2.0;
+	*(dl1dl6+NUMNODESP1*2+4)=0.5*A2;
 
 	/*Sixth nodal function: The corresponding nodal function is N=A3*(1+z)/2. Its derivatives follow*/
-	*(dl1dl6+numgrids*0+5)=0.0;
-	*(dl1dl6+numgrids*1+5)=1.0/SQRT3*(1.0+z)/2.0;
-	*(dl1dl6+numgrids*2+5)=0.5*A3;
+	*(dl1dl6+NUMNODESP1*0+5)=0.0;
+	*(dl1dl6+NUMNODESP1*1+5)=1.0/SQRT3*(1.0+z)/2.0;
+	*(dl1dl6+NUMNODESP1*2+5)=0.5*A3;
 }
 /*}}}*/
@@ -1010,5 +961,5 @@
 	 *   p is a vector of size 3x1 already allocated.
 	 */
-	double dh1dh6[3][6];
+	double dh1dh6[3][NUMNODESP1];
 
 	/*Get nodal funnctions derivatives in actual coordinate system: */
Index: /issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 5744)
+++ /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 5745)
@@ -489,5 +489,5 @@
 
 					/*build new bed and surface: */
-					if (this->GetShelf()){
+					if (this->IsOnShelf()){
 						/*hydrostatic equilibrium: */
 						double rho_ice,rho_water,di;
@@ -819,6 +819,6 @@
 }
 /*}}}*/
-/*FUNCTION Tria::GetOnBed {{{1*/
-bool Tria::GetOnBed(){
+/*FUNCTION Tria::IsOnBed {{{1*/
+bool Tria::IsOnBed(){
 	
 	bool onbed;
@@ -828,6 +828,6 @@
 }
 /*}}}*/
-/*FUNCTION Tria::GetShelf {{{1*/
-bool   Tria::GetShelf(){
+/*FUNCTION Tria::IsOnShelf {{{1*/
+bool   Tria::IsOnShelf(){
 	/*inputs: */
 	bool shelf;
@@ -2421,5 +2421,5 @@
 
 	/*If shelf: hydrostatic equilibrium*/
-	if (this->GetShelf()){
+	if (this->IsOnShelf()){
 
 		/*recover material parameters: */
Index: /issm/trunk/src/c/objects/Elements/Tria.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.h	(revision 5744)
+++ /issm/trunk/src/c/objects/Elements/Tria.h	(revision 5745)
@@ -75,6 +75,6 @@
 		void   CreatePVector(Vec pg);
 		int    GetNodeIndex(Node* node);
-		bool   GetOnBed();
-		bool   GetShelf(); 
+		bool   IsOnBed();
+		bool   IsOnShelf(); 
 		void   GetSolutionFromInputs(Vec solution);
 		void   GetVectorFromInputs(Vec vector,int NameEnum);
Index: /issm/trunk/src/c/objects/Elements/TriaRef.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/TriaRef.cpp	(revision 5744)
+++ /issm/trunk/src/c/objects/Elements/TriaRef.cpp	(revision 5745)
@@ -20,4 +20,11 @@
 /*}}}*/
 
+/*Element macros*/
+#define NUMNODES 3
+#define NDOF1 1
+#define NDOF2 2
+#define NDOF3 3
+#define NDOF4 4
+
 /*Object constructors and destructor*/
 /*FUNCTION TriaRef::TriaRef(){{{1*/
@@ -27,4 +34,5 @@
 /*}}}*/
 /*FUNCTION TriaRef::TriaRef(int* types,int nummodels){{{1*/
+
 TriaRef::TriaRef(const int nummodels){
 
@@ -62,13 +70,9 @@
 	 * where h is the interpolation function for grid i.
 	 *
-	 * We assume B has been allocated already, of size: 3x(NDOF2*numgrids)
+	 * We assume B has been allocated already, of size: 3x(NDOF2*NUMNODES)
 	 */
 
 	int i;
-	const int NDOF2=2;
-	const int numgrids=3;
-
-	double dh1dh3[NDOF2][numgrids];
-
+	double dh1dh3[NDOF2][NUMNODES];
 
 	/*Get dh1dh2dh3 in actual coordinate system: */
@@ -76,11 +80,11 @@
 
 	/*Build B: */
-	for (i=0;i<numgrids;i++){
-		*(B+NDOF2*numgrids*0+NDOF2*i)=dh1dh3[0][i]; //B[0][NDOF2*i]=dh1dh3[0][i];
-		*(B+NDOF2*numgrids*0+NDOF2*i+1)=0;
-		*(B+NDOF2*numgrids*1+NDOF2*i)=0;
-		*(B+NDOF2*numgrids*1+NDOF2*i+1)=dh1dh3[1][i];
-		*(B+NDOF2*numgrids*2+NDOF2*i)=(float).5*dh1dh3[1][i]; 
-		*(B+NDOF2*numgrids*2+NDOF2*i+1)=(float).5*dh1dh3[0][i]; 
+	for (i=0;i<NUMNODES;i++){
+		*(B+NDOF2*NUMNODES*0+NDOF2*i)=dh1dh3[0][i]; //B[0][NDOF2*i]=dh1dh3[0][i];
+		*(B+NDOF2*NUMNODES*0+NDOF2*i+1)=0;
+		*(B+NDOF2*NUMNODES*1+NDOF2*i)=0;
+		*(B+NDOF2*NUMNODES*1+NDOF2*i+1)=dh1dh3[1][i];
+		*(B+NDOF2*NUMNODES*2+NDOF2*i)=(float).5*dh1dh3[1][i]; 
+		*(B+NDOF2*NUMNODES*2+NDOF2*i+1)=(float).5*dh1dh3[0][i]; 
 	}
 }
@@ -95,6 +99,5 @@
 	 */
 
-	const int numgrids=3;
-	double l1l3[numgrids];
+	double l1l3[NUMNODES];
 
 	GetNodalFunctions(&l1l3[0],gauss);
@@ -115,6 +118,5 @@
 	 */
 
-	const int numgrids=3;
-	double l1l3[numgrids];
+	double l1l3[NUMNODES];
 
 	GetNodalFunctions(&l1l3[0],gauss);
@@ -135,13 +137,8 @@
 	 * where h is the interpolation function for grid i.
 	 *
-	 * We assume B_prog has been allocated already, of size: 2x(NDOF1*numgrids)
-	 */
-
-	int i;
-	const int NDOF1=1;
-	const int numgrids=3;
-
-	double l1l2l3[numgrids];
-
+	 * We assume B_prog has been allocated already, of size: 2x(NDOF1*NUMNODES)
+	 */
+
+	double l1l2l3[NUMNODES];
 
 	/*Get dh1dh2dh3 in actual coordinate system: */
@@ -149,7 +146,7 @@
 
 	/*Build B_prog: */
-	for (i=0;i<numgrids;i++){
-		*(B_prog+NDOF1*numgrids*0+NDOF1*i)=l1l2l3[i];
-		*(B_prog+NDOF1*numgrids*1+NDOF1*i)=l1l2l3[i];
+	for (int i=0;i<NUMNODES;i++){
+		*(B_prog+NDOF1*NUMNODES*0+NDOF1*i)=l1l2l3[i];
+		*(B_prog+NDOF1*NUMNODES*1+NDOF1*i)=l1l2l3[i];
 	}
 }
@@ -166,13 +163,9 @@
 	 * where h is the interpolation function for grid i.
 	 *
-	 * We assume B' has been allocated already, of size: 3x(NDOF2*numgrids)
-	 */
-
-	int i;
-	const int NDOF2=2;
-	const int numgrids=3;
+	 * We assume B' has been allocated already, of size: 3x(NDOF2*NUMNODES)
+	 */
 
 	/*Same thing in the actual coordinate system: */
-	double dh1dh3[NDOF2][numgrids];
+	double dh1dh3[NDOF2][NUMNODES];
 
 	/*Get dh1dh2dh3 in actual coordinates system : */
@@ -180,11 +173,11 @@
 
 	/*Build B': */
-	for (i=0;i<numgrids;i++){
-		*(Bprime+NDOF2*numgrids*0+NDOF2*i)=2*dh1dh3[0][i]; 
-		*(Bprime+NDOF2*numgrids*0+NDOF2*i+1)=dh1dh3[1][i]; 
-		*(Bprime+NDOF2*numgrids*1+NDOF2*i)=dh1dh3[0][i]; 
-		*(Bprime+NDOF2*numgrids*1+NDOF2*i+1)=2*dh1dh3[1][i]; 
-		*(Bprime+NDOF2*numgrids*2+NDOF2*i)=dh1dh3[1][i]; 
-		*(Bprime+NDOF2*numgrids*2+NDOF2*i+1)=dh1dh3[0][i]; 
+	for (int i=0;i<NUMNODES;i++){
+		*(Bprime+NDOF2*NUMNODES*0+NDOF2*i)=2*dh1dh3[0][i]; 
+		*(Bprime+NDOF2*NUMNODES*0+NDOF2*i+1)=dh1dh3[1][i]; 
+		*(Bprime+NDOF2*NUMNODES*1+NDOF2*i)=dh1dh3[0][i]; 
+		*(Bprime+NDOF2*NUMNODES*1+NDOF2*i+1)=2*dh1dh3[1][i]; 
+		*(Bprime+NDOF2*NUMNODES*2+NDOF2*i)=dh1dh3[1][i]; 
+		*(Bprime+NDOF2*NUMNODES*2+NDOF2*i+1)=dh1dh3[0][i]; 
 	}
 }
@@ -199,14 +192,9 @@
 	 * where h is the interpolation function for grid i.
 	 *
-	 * We assume B' has been allocated already, of size: 3x(NDOF2*numgrids)
-	 */
-
-	int i;
-	const int NDOF1=1;
-	const int NDOF2=2;
-	const int numgrids=3;
+	 * We assume B' has been allocated already, of size: 3x(NDOF2*NUMNODES)
+	 */
 
 	/*Same thing in the actual coordinate system: */
-	double dh1dh3[NDOF2][numgrids];
+	double dh1dh3[NDOF2][NUMNODES];
 
 	/*Get dh1dh2dh3 in actual coordinates system : */
@@ -214,7 +202,7 @@
 
 	/*Build B': */
-	for (i=0;i<numgrids;i++){
-		*(Bprime_prog+NDOF1*numgrids*0+NDOF1*i)=dh1dh3[0][i]; 
-		*(Bprime_prog+NDOF1*numgrids*1+NDOF1*i)=dh1dh3[1][i]; 
+	for (int i=0;i<NUMNODES;i++){
+		*(Bprime_prog+NDOF1*NUMNODES*0+NDOF1*i)=dh1dh3[0][i]; 
+		*(Bprime_prog+NDOF1*NUMNODES*1+NDOF1*i)=dh1dh3[1][i]; 
 	}
 }
@@ -232,10 +220,8 @@
 	 * where h is the interpolation function for grid i.
 	 *
-	 * We assume L has been allocated already, of size: numgrids (numdof=1), or numdofx(numdof*numgrids) (numdof=2)
+	 * We assume L has been allocated already, of size: NUMNODES (numdof=1), or numdofx(numdof*NUMNODES) (numdof=2)
 	 */
 
 	int i;
-	const int NDOF2=2;
-	const int numgrids=3;
 	double l1l2l3[3];
 
@@ -245,14 +231,14 @@
 	/*Build L: */
 	if(numdof==1){
-		for (i=0;i<numgrids;i++){
+		for (i=0;i<NUMNODES;i++){
 			L[i]=l1l2l3[i]; 
 		}
 	}
 	else{
-		for (i=0;i<numgrids;i++){
-			*(L+numdof*numgrids*0+numdof*i)=l1l2l3[i]; //L[0][NDOF2*i]=dh1dh3[0][i];
-			*(L+numdof*numgrids*0+numdof*i+1)=0;
-			*(L+numdof*numgrids*1+numdof*i)=0;
-			*(L+numdof*numgrids*1+numdof*i+1)=l1l2l3[i];
+		for (i=0;i<NUMNODES;i++){
+			*(L+numdof*NUMNODES*0+numdof*i)=l1l2l3[i]; //L[0][NDOF2*i]=dh1dh3[0][i];
+			*(L+numdof*NUMNODES*0+numdof*i+1)=0;
+			*(L+numdof*NUMNODES*1+numdof*i)=0;
+			*(L+numdof*NUMNODES*1+numdof*i+1)=l1l2l3[i];
 		}
 	}
@@ -263,15 +249,12 @@
 	/*The Jacobian is constant over the element, discard the gaussian points. 
 	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
-
-	const int NDOF2=2;
-	const int numgrids=3;
 	double x1,y1,x2,y2,x3,y3;
 
-	x1=*(xyz_list+numgrids*0+0);
-	y1=*(xyz_list+numgrids*0+1);
-	x2=*(xyz_list+numgrids*1+0);
-	y2=*(xyz_list+numgrids*1+1);
-	x3=*(xyz_list+numgrids*2+0);
-	y3=*(xyz_list+numgrids*2+1);
+	x1=*(xyz_list+NUMNODES*0+0);
+	y1=*(xyz_list+NUMNODES*0+1);
+	x2=*(xyz_list+NUMNODES*1+0);
+	y2=*(xyz_list+NUMNODES*1+1);
+	x3=*(xyz_list+NUMNODES*2+0);
+	y3=*(xyz_list+NUMNODES*2+1);
 
 
@@ -388,7 +371,5 @@
 	 * actual coordinate system): */
 	int       i;
-	const int NDOF2    = 2;
-	const int numgrids = 3;
-	double    dh1dh3_ref[NDOF2][numgrids];
+	double    dh1dh3_ref[NDOF2][NUMNODES];
 	double    Jinv[NDOF2][NDOF2];
 
@@ -404,7 +385,7 @@
 	 * [dhi/dy]       [dhi/ds]
 	 */
-	for (i=0;i<numgrids;i++){
-		dh1dh3[numgrids*0+i]=Jinv[0][0]*dh1dh3_ref[0][i]+Jinv[0][1]*dh1dh3_ref[1][i];
-		dh1dh3[numgrids*1+i]=Jinv[1][0]*dh1dh3_ref[0][i]+Jinv[1][1]*dh1dh3_ref[1][i];
+	for (i=0;i<NUMNODES;i++){
+		dh1dh3[NUMNODES*0+i]=Jinv[0][0]*dh1dh3_ref[0][i]+Jinv[0][1]*dh1dh3_ref[1][i];
+		dh1dh3[NUMNODES*1+i]=Jinv[1][0]*dh1dh3_ref[0][i]+Jinv[1][1]*dh1dh3_ref[1][i];
 	}
 
@@ -413,22 +394,18 @@
 /*FUNCTION TriaRef::GetNodalFunctionsDerivativesReference{{{1*/
 void TriaRef::GetNodalFunctionsDerivativesReference(double* dl1dl3,GaussTria* gauss){
-
 	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
 	 * natural coordinate system) at the gaussian point. */
 
-	const int NDOF2=2;
-	const int numgrids=3;
-
 	/*First nodal function: */
-	*(dl1dl3+numgrids*0+0)=-0.5; 
-	*(dl1dl3+numgrids*1+0)=-1.0/(2.0*SQRT3);
+	*(dl1dl3+NUMNODES*0+0)=-0.5; 
+	*(dl1dl3+NUMNODES*1+0)=-1.0/(2.0*SQRT3);
 
 	/*Second nodal function: */
-	*(dl1dl3+numgrids*0+1)=0.5;
-	*(dl1dl3+numgrids*1+1)=-1.0/(2.0*SQRT3);
+	*(dl1dl3+NUMNODES*0+1)=0.5;
+	*(dl1dl3+NUMNODES*1+1)=-1.0/(2.0*SQRT3);
 
 	/*Third nodal function: */
-	*(dl1dl3+numgrids*0+2)=0;
-	*(dl1dl3+numgrids*1+2)=1.0/SQRT3;
+	*(dl1dl3+NUMNODES*0+2)=0;
+	*(dl1dl3+NUMNODES*1+2)=1.0/SQRT3;
 
 }