Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 5661)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 5662)
@@ -18,4 +18,7 @@
 #include "../../Container/Container.h"
 /*}}}*/
+
+/*Element macros*/
+#define NUMVERTICES 6
 
 /*Penta constructors and destructor*/
@@ -482,7 +485,6 @@
 
 	int i,j;
-	const int numgrids=6;
-	int    doflist[numgrids];
-	double xyz_list[numgrids][3];
+	int    doflist[NUMVERTICES];
+	double xyz_list[NUMVERTICES][3];
 	double xyz_list_tria[3][3];
 
@@ -555,5 +557,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	for(i=0;i<3;i++){
 		for(j=0;j<3;j++){
@@ -968,6 +970,5 @@
 
 	int i;
-	const int numvertices=6;
-	int doflist1[numvertices];
+	int doflist1[NUMVERTICES];
 
 	/*Find NameEnum input in the inputs dataset, and get it to fill in the vector: */
@@ -1107,5 +1108,4 @@
 
 	/*Intermediaries*/
-	const int numvertices=6;
 	bool onbed;
 	int  step,i;
@@ -1119,7 +1119,7 @@
 	Input* depth_averaged_input=NULL;
 
-	double  xyz_list[numvertices][3];
-	double  Helem_list[numvertices];
-	double  zeros_list[numvertices]={0.0};
+	double  xyz_list[NUMVERTICES][3];
+	double  Helem_list[NUMVERTICES];
+	double  zeros_list[NUMVERTICES]={0.0};
 
 	/*recover parameters: */
@@ -1157,5 +1157,5 @@
 
 		/*Step2: Create element thickness input*/
-		GetVerticesCoordinates(&xyz_list[0][0],penta->nodes,numvertices);
+		GetVerticesCoordinates(&xyz_list[0][0],penta->nodes,NUMVERTICES);
 		for(i=0;i<3;i++){
 			Helem_list[i]=xyz_list[i+3][2]-xyz_list[i][2];
@@ -1705,5 +1705,4 @@
 	/*output: */
 	int     numrows     = 0;
-	int     numvertices = 0;
 	int     numnodes    = 0;
 
@@ -1714,5 +1713,4 @@
 		numrows++;
 		/*now, how many vertices and how many nodal values for this result? :*/
-		numvertices=6; //this is a penta element, with 6 vertices
 		numnodes=elementresult->NumberOfNodalValues(); //ask result object.
 	}
@@ -1720,5 +1718,5 @@
 	/*Assign output pointers:*/
 	*pnumrows=numrows;
-	*pnumvertices=numvertices;
+	*pnumvertices=NUMVERTICES;
 	*pnumnodes=numnodes;
 	
@@ -2153,9 +2151,9 @@
 
 	/* node data: */
-	const int    numgridsm=3;  //MacAyealnumgrids
-	const int    numdofm=2*numgridsm;
-	const int    numgridsp=6; //Pattyn numgrids
-	const int    numdofp=2*numgridsp;
-	double       xyz_list[numgridsp][3];
+	const int    NUMVERTICESm=3;  //MacAyealNUMVERTICES
+	const int    numdofm=2*NUMVERTICESm;
+	const int    NUMVERTICESp=6; //Pattyn NUMVERTICES
+	const int    numdofp=2*NUMVERTICESp;
+	double       xyz_list[NUMVERTICESp][3];
 	int*         doflistm=NULL;
 	int*         doflistp=NULL;
@@ -2226,5 +2224,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgridsp);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICESp);
 	tria->GetDofList(&doflistm,MacAyealApproximationEnum);  //Pattyn dof list
 	GetDofList(&doflistp,PattynApproximationEnum); //MacAyeal dof list
@@ -2304,7 +2302,6 @@
 	/*Collapsed formulation: */
 	int       i;
-	const int numgrids=6;
 	const int NDOF2=2;
-	const int numdofs=NDOF2*numgrids;
+	const int numdofs=NDOF2*NUMVERTICES;
 	int*      doflist=NULL;
 	double    Ke_gg[numdofs][numdofs]={0.0};
@@ -2389,8 +2386,8 @@
 
 	/* node data: */
-	const int    numgrids3d=6;
-	const int    numgrids2d=3;
-	const int    numdof2d=2*numgrids2d;
-	double       xyz_list[numgrids3d][3];
+	const int    NUMVERTICES3d=6;
+	const int    NUMVERTICES2d=3;
+	const int    numdof2d=2*NUMVERTICES2d;
+	double       xyz_list[NUMVERTICES3d][3];
 	int*         doflist=NULL;
 
@@ -2500,5 +2497,5 @@
 
 		/* Get node coordinates and dof list: */
-		GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids3d);
+		GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES3d);
 		tria->GetDofList(&doflist,MacAyealApproximationEnum);  //Pattyn dof list
 
@@ -2576,7 +2573,6 @@
 
 	/* node data: */
-	const int    numgrids=6;
-	const int    numdof=2*numgrids;
-	double       xyz_list[numgrids][3];
+	const int    numdof=2*NUMVERTICES;
+	double       xyz_list[NUMVERTICES][3];
 	int*         doflist=NULL;
 
@@ -2648,5 +2644,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	GetDofList(&doflist,PattynApproximationEnum);
 
@@ -2719,11 +2715,10 @@
 
 	int i,j;
-	const int numgrids=6;
 	const int DOFPERGRID=4;
-	const int numdof=numgrids*DOFPERGRID;
+	const int numdof=NUMVERTICES*DOFPERGRID;
 	int*      doflist=NULL;
 
-	const int numgrids2d=3;
-	const int numdof2d=numgrids2d*DOFPERGRID;
+	const int NUMVERTICES2d=3;
+	const int numdof2d=NUMVERTICES2d*DOFPERGRID;
 
 	/*Collapsed formulation: */
@@ -2731,6 +2726,6 @@
 
 	/*Grid data: */
-	double     xyz_list[numgrids][3];
-	double	  xyz_list_tria[numgrids2d][3];
+	double     xyz_list[NUMVERTICES][3];
+	double	  xyz_list_tria[NUMVERTICES2d][3];
 	double	  bed_normal[3];
 
@@ -2759,5 +2754,5 @@
 	double* area_gauss_weights  =  NULL;
 	double* vert_gauss_weights  =  NULL;
-	double  gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
+	double  gaussgrids[NUMVERTICES][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
 
 	/* specific gaussian point: */
@@ -2802,5 +2797,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	GetDofList(&doflist,StokesApproximationEnum);
 
@@ -2877,5 +2872,5 @@
 		friction=new Friction("3d",inputs,matpar,analysis_type);
 
-		for(i=0;i<numgrids2d;i++){
+		for(i=0;i<NUMVERTICES2d;i++){
 			for(j=0;j<3;j++){
 				xyz_list_tria[i][j]=xyz_list[i][j];
@@ -2975,8 +2970,7 @@
 
 	/* node data: */
-	const int    numgrids=6;
 	const int    NDOF1=1;
-	const int    numdof=NDOF1*numgrids;
-	double       xyz_list[numgrids][3];
+	const int    numdof=NDOF1*NUMVERTICES;
+	double       xyz_list[NUMVERTICES][3];
 	int*         doflist=NULL;
 
@@ -2989,6 +2983,6 @@
 	double  Ke_gg_gaussian[numdof][numdof];
 	double  Jdet;
-	double  B[NDOF1][numgrids];
-	double  Bprime[NDOF1][numgrids];
+	double  B[NDOF1][NUMVERTICES];
+	double  Bprime[NDOF1][NUMVERTICES];
 	double  DL_scalar;
 
@@ -3018,5 +3012,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	GetDofList(&doflist);
 
@@ -3036,7 +3030,7 @@
 
 		/*  Do the triple product tB*D*Bprime: */
-		TripleMultiply( &B[0][0],1,numgrids,1,
+		TripleMultiply( &B[0][0],1,NUMVERTICES,1,
 					&DL_scalar,1,1,0,
-					&Bprime[0][0],1,numgrids,0,
+					&Bprime[0][0],1,NUMVERTICES,0,
 					&Ke_gg_gaussian[0][0],0);
 
@@ -3154,8 +3148,7 @@
 
 	/* node data: */
-	const int    numgrids=6;
 	const int    NDOF1=1;
-	const int    numdof=NDOF1*numgrids;
-	double       xyz_list[numgrids][3];
+	const int    numdof=NDOF1*NUMVERTICES;
+	double       xyz_list[NUMVERTICES][3];
 	int*         doflist=NULL;
 
@@ -3222,5 +3215,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	GetDofList(&doflist);
 
@@ -3487,11 +3480,10 @@
 	int i,j,k;
 	
-	const int numgrids=6;
 	const int NDOF2=2;
-	const int numdofs=NDOF2*numgrids;
+	const int numdofs=NDOF2*NUMVERTICES;
 	int*      doflist=NULL;
 	double    pe_g[numdofs]={0.0};
-	double    xyz_list[numgrids][3];
-	double    z_list[numgrids];
+	double    xyz_list[NUMVERTICES][3];
+	double    z_list[NUMVERTICES];
 	double    z_segment[2];
 	double    slope2,constant_part;
@@ -3555,6 +3547,6 @@
 	slopey_input=inputs->GetInput(SurfaceSlopeYEnum);
 
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
-	for(i=0;i<numgrids;i++)z_list[i]=xyz_list[i][2];
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
+	for(i=0;i<NUMVERTICES;i++)z_list[i]=xyz_list[i][2];
 
 	//Get gaussian points and weights. order 2 since we have a product of 2 nodal functions
@@ -3675,8 +3667,7 @@
 
 	/* node data: */
-	const int    numgrids=6;
 	const int    NDOF2=2;
-	const int    numdof=NDOF2*numgrids;
-	double       xyz_list[numgrids][3];
+	const int    numdof=NDOF2*NUMVERTICES;
+	double       xyz_list[NUMVERTICES][3];
 	int*         doflist=NULL;
 
@@ -3714,5 +3705,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	GetDofList(&doflist,PattynApproximationEnum);
 
@@ -3743,5 +3734,5 @@
 
 		/*Build pe_g_gaussian vector: */
-		for (i=0;i<numgrids;i++){
+		for (i=0;i<NUMVERTICES;i++){
 			for (j=0;j<NDOF2;j++){
 				pe_g_gaussian[i*NDOF2+j]=-driving_stress_baseline*slope[j]*Jdet*gauss->weight*l1l6[i];
@@ -3766,9 +3757,8 @@
 	int i,j;
 
-	const int numgrids=6;
 	const int DOFPERGRID=4;
-	const int numdof=numgrids*DOFPERGRID;
-	const int numgrids2d=3;
-	int       numdof2d=numgrids2d*DOFPERGRID;
+	const int numdof=NUMVERTICES*DOFPERGRID;
+	const int NUMVERTICES2d=3;
+	int       numdof2d=NUMVERTICES2d*DOFPERGRID;
 	int*      doflist=NULL;
 
@@ -3777,6 +3767,6 @@
 
 	/*parameters: */
-	double		   xyz_list_tria[numgrids2d][3];
-	double         xyz_list[numgrids][3];
+	double		   xyz_list_tria[NUMVERTICES2d][3];
+	double         xyz_list[NUMVERTICES][3];
 	double		   bed_normal[3];
 	double         bed;
@@ -3857,5 +3847,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	GetDofList(&doflist,StokesApproximationEnum);
 
@@ -3898,5 +3888,5 @@
 
 			/* Build gaussian vector */
-			for(i=0;i<numgrids+1;i++){
+			for(i=0;i<NUMVERTICES+1;i++){
 				Pe_gaussian[i*DOFPERGRID+2]=-rho_ice*gravity*Jdet*gauss_weight*l1l7[i];
 			}
@@ -3944,5 +3934,5 @@
 	if ( (onbed==1) && (shelf==1)){
 
-		for(i=0;i<numgrids2d;i++){
+		for(i=0;i<NUMVERTICES2d;i++){
 			for(j=0;j<3;j++){
 				xyz_list_tria[i][j]=xyz_list[i][j];
@@ -3980,5 +3970,5 @@
 			BedNormal(&bed_normal[0],xyz_list_tria);
 
-			for(i=0;i<numgrids2d;i++){
+			for(i=0;i<NUMVERTICES2d;i++){
 				for(j=0;j<3;j++){
 					Pe_temp[i*DOFPERGRID+j]+=water_pressure*gauss_weight*Jdet2d*L[i]*bed_normal[j];
@@ -4035,8 +4025,7 @@
 
 	/* node data: */
-	const int    numgrids=6;
 	const int    NDOF1=1;
-	const int    numdof=NDOF1*numgrids;
-	double       xyz_list[numgrids][3];
+	const int    numdof=NDOF1*NUMVERTICES;
+	double       xyz_list[NUMVERTICES][3];
 	int*         doflist=NULL;
 
@@ -4084,5 +4073,5 @@
 	/*Now, handle the standard penta element: */
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	GetDofList(&doflist);
 
@@ -4110,5 +4099,5 @@
 
 		/*Build pe_g_gaussian vector: */
-		for (i=0;i<numgrids;i++){
+		for (i=0;i<NUMVERTICES;i++){
 			pe_g_gaussian[i]=(dudx+dvdy)*Jdet*gauss->weight*l1l6[i];
 		}
@@ -4200,11 +4189,10 @@
 	int found=0;
 
-	const int  numgrids=6;
 	const int  NDOF1=1;
-	const int  numdof=numgrids*NDOF1;
+	const int  numdof=NUMVERTICES*NDOF1;
 	int*       doflist=NULL;
 
 	/*Grid data: */
-	double        xyz_list[numgrids][3];
+	double        xyz_list[NUMVERTICES][3];
 
 	/* gaussian points: */
@@ -4212,5 +4200,5 @@
 	GaussPenta *gauss=NULL;
 
-	double temperature_list[numgrids];
+	double temperature_list[NUMVERTICES];
 	double temperature;
 
@@ -4269,5 +4257,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	GetDofList(&doflist);
 
@@ -4309,5 +4297,5 @@
 		if(dt) scalar_def=scalar_def*dt;
 
-		for(i=0;i<numgrids;i++) P_terms[i]+=scalar_def*L[i];
+		for(i=0;i<NUMVERTICES;i++) P_terms[i]+=scalar_def*L[i];
 
 		/* Build transient now */
@@ -4315,5 +4303,5 @@
 			temperature_input->GetParameterValue(&temperature, gauss);
 			scalar_transient=temperature*Jdet*gauss->weight;
-			for(i=0;i<numgrids;i++) P_terms[i]+=scalar_transient*L[i];
+			for(i=0;i<NUMVERTICES;i++) P_terms[i]+=scalar_transient*L[i];
 		}
 	}
@@ -4435,6 +4423,5 @@
 
 	/*Intermediaries*/
-	const int  numvertices = 6;
-	double     value[numvertices];
+	double     value[NUMVERTICES];
 	GaussPenta *gauss              = NULL;
 
@@ -4448,5 +4435,5 @@
 	/* Start looping on the number of vertices: */
 	gauss=new GaussPenta();
-	for (int iv=0;iv<numvertices;iv++){
+	for (int iv=0;iv<NUMVERTICES;iv++){
 		gauss->GaussVertex(iv);
 		input->GetParameterValue(&pvalue[iv],gauss);
@@ -4461,6 +4448,5 @@
 
 	/*Intermediaries*/
-	const int  numvertices = 6;
-	double     value[numvertices];
+	double     value[NUMVERTICES];
 	GaussPenta *gauss              = NULL;
 
@@ -4474,5 +4460,5 @@
 	if (input){
 		gauss=new GaussPenta();
-		for (int iv=0;iv<numvertices;iv++){
+		for (int iv=0;iv<NUMVERTICES;iv++){
 			gauss->GaussVertex(iv);
 			input->GetParameterValue(&pvalue[iv],gauss);
@@ -4480,5 +4466,5 @@
 	}
 	else{
-		for (int iv=0;iv<numvertices;iv++) pvalue[iv]=defaultvalue;
+		for (int iv=0;iv<NUMVERTICES;iv++) pvalue[iv]=defaultvalue;
 	}
 
@@ -4656,8 +4642,7 @@
 	int i;
 
-	const int    numvertices=6;
 	const int    numdofpervertex=2;
-	const int    numdof=numdofpervertex*numvertices;
-	double       gauss[numvertices][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
+	const int    numdof=numdofpervertex*NUMVERTICES;
+	double       gauss[NUMVERTICES][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
 	int*         doflist=NULL;
 	double       values[numdof];
@@ -4680,5 +4665,5 @@
 	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
 	/*P1 element only for now*/
-	for(i=0;i<numvertices;i++){
+	for(i=0;i<NUMVERTICES;i++){
 
 		/*Recover vx and vy*/
@@ -4701,8 +4686,7 @@
 	int i;
 
-	const int    numvertices=6;
 	const int    numdofpervertex=2;
-	const int    numdof=numdofpervertex*numvertices;
-	double       gauss[numvertices][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
+	const int    numdof=numdofpervertex*NUMVERTICES;
+	double       gauss[NUMVERTICES][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
 	int*         doflist=NULL;
 	double       values[numdof];
@@ -4715,5 +4699,5 @@
 	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
 	/*P1 element only for now*/
-	for(i=0;i<numvertices;i++){
+	for(i=0;i<NUMVERTICES;i++){
 
 		/*Recover vx and vy*/
@@ -4736,8 +4720,7 @@
 	int i;
 
-	const int    numvertices=6;
 	const int    numdofpervertex=1;
-	const int    numdof=numdofpervertex*numvertices;
-	double       gauss[numvertices][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
+	const int    numdof=numdofpervertex*NUMVERTICES;
+	double       gauss[NUMVERTICES][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
 	int*         doflist=NULL;
 	double       values[numdof];
@@ -4749,5 +4732,5 @@
 	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
 	/*P1 element only for now*/
-	for(i=0;i<numvertices;i++){
+	for(i=0;i<NUMVERTICES;i++){
 
 		/*Recover vz */
@@ -4768,8 +4751,7 @@
 	int i;
 
-	const int    numvertices=6;
 	const int    numdofpervertex=4;
-	const int    numdof=numdofpervertex*numvertices;
-	double       gauss[numvertices][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
+	const int    numdof=numdofpervertex*NUMVERTICES;
+	double       gauss[NUMVERTICES][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
 	int*         doflist=NULL;
 	double       values[numdof];
@@ -4785,5 +4767,5 @@
 	/*Ok, we have vx vy vz and P in values, fill in vx vy vz P arrays: */
 	/*P1 element only for now*/
-	for(i=0;i<numvertices;i++){
+	for(i=0;i<NUMVERTICES;i++){
 
 		/*Recover vx and vy*/
@@ -4810,8 +4792,7 @@
 	int i;
 
-	const int    numvertices=6;
 	const int    numdofpervertex=1;
-	const int    numdof=numdofpervertex*numvertices;
-	double       gauss[numvertices][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
+	const int    numdof=numdofpervertex*NUMVERTICES;
+	double       gauss[NUMVERTICES][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
 	int*         doflist=NULL;
 	double       values[numdof];
@@ -4824,5 +4805,5 @@
 	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
 	/*P1 element only for now*/
-	for(i=0;i<numvertices;i++){
+	for(i=0;i<NUMVERTICES;i++){
 
 		/*Recover vz */
@@ -5196,19 +5177,18 @@
 	int i;
 
-	const int    numvertices=6;
 	const int    numdofpervertex=2;
-	const int    numdof=numdofpervertex*numvertices;
+	const int    numdof=numdofpervertex*NUMVERTICES;
 	int*         doflist=NULL;
 	double       values[numdof];
-	double       vx[numvertices];
-	double       vy[numvertices];
-	double       vz[numvertices];
-	double       vel[numvertices];
+	double       vx[NUMVERTICES];
+	double       vy[NUMVERTICES];
+	double       vz[NUMVERTICES];
+	double       vel[NUMVERTICES];
 	int          dummy;
-	double       pressure[numvertices];
-	double       surface[numvertices];
+	double       pressure[NUMVERTICES];
+	double       surface[NUMVERTICES];
 	double       rho_ice,g;
-	double       xyz_list[numvertices][3];
-	double       gauss[numvertices][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
+	double       xyz_list[NUMVERTICES][3];
+	double       gauss[NUMVERTICES][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
 
 	Input  *vz_input        = NULL;
@@ -5243,5 +5223,5 @@
 
 		/*Get node data: */
-		GetVerticesCoordinates(&xyz_list[0][0],penta->nodes,numvertices);
+		GetVerticesCoordinates(&xyz_list[0][0],penta->nodes,NUMVERTICES);
 
 		/*Now Compute vel*/
@@ -5250,12 +5230,12 @@
 			if (vz_input->Enum()!=PentaVertexInputEnum) ISSMERROR("Cannot compute Vel as Vz is of type %s",EnumToString(vz_input->Enum()));
 			vz_input->GetValuesPtr(&vz_ptr,&dummy);
-			for(i=0;i<numvertices;i++) vz[i]=vz_ptr[i];
-		}
-		else{for(i=0;i<numvertices;i++) vz[i]=0.0;}
-		for(i=0;i<numvertices;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
+			for(i=0;i<NUMVERTICES;i++) vz[i]=vz_ptr[i];
+		}
+		else{for(i=0;i<NUMVERTICES;i++) vz[i]=0.0;}
+		for(i=0;i<NUMVERTICES;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
 
 		/*Now compute pressure*/
 		GetParameterListOnVertices(&surface[0],SurfaceEnum);
-		for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]);
+		for(i=0;i<NUMVERTICES;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]);
 
 		/*Now, we have to move the previous Vx and Vy inputs  to old 
@@ -5287,23 +5267,22 @@
 	int i;
 
-	const int    numvertices=6;
-	const int    numvertices2d=3;
+	const int    NUMVERTICES2d=3;
 	const int    numdofpervertex=2;
-	const int    numdof=numdofpervertex*numvertices;
-	const int    numdof2d=numdofpervertex*numvertices2d;
+	const int    numdof=numdofpervertex*NUMVERTICES;
+	const int    numdof2d=numdofpervertex*NUMVERTICES2d;
 	int*         doflistp=NULL;
 	int*         doflistm=NULL;
 	double       macayeal_values[numdof];
 	double       pattyn_values[numdof];
-	double       vx[numvertices];
-	double       vy[numvertices];
-	double       vz[numvertices];
-	double       vel[numvertices];
+	double       vx[NUMVERTICES];
+	double       vy[NUMVERTICES];
+	double       vz[NUMVERTICES];
+	double       vel[NUMVERTICES];
 	int          dummy;
-	double       pressure[numvertices];
-	double       surface[numvertices];
+	double       pressure[NUMVERTICES];
+	double       surface[NUMVERTICES];
 	double       rho_ice,g;
-	double       xyz_list[numvertices][3];
-	double       gauss[numvertices][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
+	double       xyz_list[NUMVERTICES][3];
+	double       gauss[NUMVERTICES][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
 
 	Input  *vz_input       = NULL;
@@ -5320,5 +5299,5 @@
 
 	/*Get node data: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 
 	/*Use the dof list to index into the solution vector: */
@@ -5333,5 +5312,5 @@
 
 	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
-	for(i=0;i<numvertices;i++){
+	for(i=0;i<NUMVERTICES;i++){
 		vx[i]=macayeal_values[i*numdofpervertex+0]+pattyn_values[i*numdofpervertex+0];
 		vy[i]=macayeal_values[i*numdofpervertex+1]+pattyn_values[i*numdofpervertex+1];
@@ -5345,12 +5324,12 @@
 		}
 		vz_input->GetValuesPtr(&vz_ptr,&dummy);
-		for(i=0;i<numvertices;i++) vz[i]=vz_ptr[i];
+		for(i=0;i<NUMVERTICES;i++) vz[i]=vz_ptr[i];
 	}
 	else{
-		for(i=0;i<numvertices;i++) vz[i]=0.0;
+		for(i=0;i<NUMVERTICES;i++) vz[i]=0.0;
 	}
 
 	/*Now Compute vel*/
-	for(i=0;i<numvertices;i++) {
+	for(i=0;i<NUMVERTICES;i++) {
 		vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
 	}
@@ -5361,5 +5340,5 @@
 	g=matpar->GetG();
 	GetParameterListOnVertices(&surface[0],SurfaceEnum);
-	for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]);
+	for(i=0;i<NUMVERTICES;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]);
 
 	/*Now, we have to move the previous Vx and Vy inputs  to old 
@@ -5385,19 +5364,18 @@
 	int i;
 
-	const int    numvertices=6;
 	const int    numdofpervertex=2;
-	const int    numdof=numdofpervertex*numvertices;
+	const int    numdof=numdofpervertex*NUMVERTICES;
 	int*         doflist=NULL;
 	double       values[numdof];
-	double       vx[numvertices];
-	double       vy[numvertices];
-	double       vz[numvertices];
-	double       vel[numvertices];
+	double       vx[NUMVERTICES];
+	double       vy[NUMVERTICES];
+	double       vz[NUMVERTICES];
+	double       vel[NUMVERTICES];
 	int          dummy;
-	double       pressure[numvertices];
-	double       surface[numvertices];
+	double       pressure[NUMVERTICES];
+	double       surface[NUMVERTICES];
 	double       rho_ice,g;
-	double       xyz_list[numvertices][3];
-	double       gauss[numvertices][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
+	double       xyz_list[NUMVERTICES][3];
+	double       gauss[NUMVERTICES][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
 	
 	Input  *vz_input        = NULL;
@@ -5408,5 +5386,5 @@
 
 	/*Get node data: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 
 	/*Use the dof list to index into the solution vector: */
@@ -5416,5 +5394,5 @@
 
 	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
-	for(i=0;i<numvertices;i++){
+	for(i=0;i<NUMVERTICES;i++){
 		vx[i]=values[i*numdofpervertex+0];
 		vy[i]=values[i*numdofpervertex+1];
@@ -5428,12 +5406,12 @@
 		}
 		vz_input->GetValuesPtr(&vz_ptr,&dummy);
-		for(i=0;i<numvertices;i++) vz[i]=vz_ptr[i];
+		for(i=0;i<NUMVERTICES;i++) vz[i]=vz_ptr[i];
 	}
 	else{
-		for(i=0;i<numvertices;i++) vz[i]=0.0;
+		for(i=0;i<NUMVERTICES;i++) vz[i]=0.0;
 	}
 
 	/*Now Compute vel*/
-	for(i=0;i<numvertices;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
+	for(i=0;i<NUMVERTICES;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
 
 	/*For pressure: we have not computed pressure in this analysis, for this element. We are in 3D, 
@@ -5442,5 +5420,5 @@
 	g=matpar->GetG();
 	GetParameterListOnVertices(&surface[0],SurfaceEnum);
-	for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]);
+	for(i=0;i<NUMVERTICES;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]);
 
 	/*Now, we have to move the previous Vx and Vy inputs  to old 
@@ -5466,19 +5444,18 @@
 	int i;
 
-	const int    numvertices=6;
 	const int    numdofpervertex=2;
-	const int    numdof=numdofpervertex*numvertices;
+	const int    numdof=numdofpervertex*NUMVERTICES;
 	int*         doflist=NULL;
 	double       values[numdof];
-	double       vx[numvertices];
-	double       vy[numvertices];
-	double       vz[numvertices];
-	double       vel[numvertices];
+	double       vx[NUMVERTICES];
+	double       vy[NUMVERTICES];
+	double       vz[NUMVERTICES];
+	double       vel[NUMVERTICES];
 	int          dummy;
-	double       pressure[numvertices];
-	double       surface[numvertices];
+	double       pressure[NUMVERTICES];
+	double       surface[NUMVERTICES];
 	double       rho_ice,g;
-	double       xyz_list[numvertices][3];
-	double       gauss[numvertices][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
+	double       xyz_list[NUMVERTICES][3];
+	double       gauss[NUMVERTICES][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
 	
 	Input*       vz_input=NULL;
@@ -5489,5 +5466,5 @@
 
 	/*Get node data: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 
 	/*Use the dof list to index into the solution vector: */
@@ -5497,5 +5474,5 @@
 
 	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
-	for(i=0;i<numvertices;i++){
+	for(i=0;i<NUMVERTICES;i++){
 		vx[i]=values[i*numdofpervertex+0];
 		vy[i]=values[i*numdofpervertex+1];
@@ -5509,12 +5486,12 @@
 		}
 		vz_input->GetValuesPtr(&vz_ptr,&dummy);
-		for(i=0;i<numvertices;i++) vz[i]=vz_ptr[i];
+		for(i=0;i<NUMVERTICES;i++) vz[i]=vz_ptr[i];
 	}
 	else{
-		for(i=0;i<numvertices;i++) vz[i]=0.0;
+		for(i=0;i<NUMVERTICES;i++) vz[i]=0.0;
 	}
 
 	/*Now Compute vel*/
-	for(i=0;i<numvertices;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
+	for(i=0;i<NUMVERTICES;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
 
 	/*For pressure: we have not computed pressure in this analysis, for this element. We are in 3D, 
@@ -5523,5 +5500,5 @@
 	g=matpar->GetG();
 	GetParameterListOnVertices(&surface[0],SurfaceEnum);
-	for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]);
+	for(i=0;i<NUMVERTICES;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]);
 
 	/*Now, we have to move the previous Vx and Vy inputs  to old 
@@ -5547,19 +5524,18 @@
 	int i;
 
-	const int    numvertices=6;
 	const int    numdofpervertex=1;
-	const int    numdof=numdofpervertex*numvertices;
+	const int    numdof=numdofpervertex*NUMVERTICES;
 	int*         doflist=NULL;
 	double       values[numdof];
-	double       vx[numvertices];
-	double       vy[numvertices];
-	double       vz[numvertices];
-	double       vel[numvertices];
+	double       vx[NUMVERTICES];
+	double       vy[NUMVERTICES];
+	double       vz[NUMVERTICES];
+	double       vel[NUMVERTICES];
 	int          dummy;
-	double       pressure[numvertices];
-	double       surface[numvertices];
+	double       pressure[NUMVERTICES];
+	double       surface[NUMVERTICES];
 	double       rho_ice,g;
-	double       xyz_list[numvertices][3];
-	double       gauss[numvertices][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
+	double       xyz_list[NUMVERTICES][3];
+	double       gauss[NUMVERTICES][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
 
 	Input*       vx_input=NULL;
@@ -5572,5 +5548,5 @@
 
 	/*Get node data: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 
 	/*Use the dof list to index into the solution vector: */
@@ -5580,5 +5556,5 @@
 
 	/*Ok, we have vz in values, fill in vz array: */
-	for(i=0;i<numvertices;i++){
+	for(i=0;i<NUMVERTICES;i++){
 		vz[i]=values[i*numdofpervertex+0];
 	}
@@ -5589,7 +5565,7 @@
 		if (vx_input->Enum()!=PentaVertexInputEnum) ISSMERROR("Cannot compute Vel as Vx is of type %s",EnumToString(vx_input->Enum()));
 		vx_input->GetValuesPtr(&vx_ptr,&dummy);
-		for(i=0;i<numvertices;i++) vx[i]=vx_ptr[i];
-	}
-	else for(i=0;i<numvertices;i++) vx[i]=0.0;
+		for(i=0;i<NUMVERTICES;i++) vx[i]=vx_ptr[i];
+	}
+	else for(i=0;i<NUMVERTICES;i++) vx[i]=0.0;
 
 	vy_input=inputs->GetInput(VyEnum);
@@ -5597,10 +5573,10 @@
 		if (vy_input->Enum()!=PentaVertexInputEnum) ISSMERROR("Cannot compute Vel as Vy is of type %s",EnumToString(vy_input->Enum()));
 		vy_input->GetValuesPtr(&vy_ptr,&dummy);
-		for(i=0;i<numvertices;i++) vy[i]=vy_ptr[i];
-	}
-	else for(i=0;i<numvertices;i++) vy[i]=0.0;
+		for(i=0;i<NUMVERTICES;i++) vy[i]=vy_ptr[i];
+	}
+	else for(i=0;i<NUMVERTICES;i++) vy[i]=0.0;
 
 	/*Now Compute vel*/
-	for(i=0;i<numvertices;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
+	for(i=0;i<NUMVERTICES;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
 
 	/*For pressure: we have not computed pressure in this analysis, for this element. We are in 3D, 
@@ -5609,5 +5585,5 @@
 	g=matpar->GetG();
 	GetParameterListOnVertices(&surface[0],SurfaceEnum);
-	for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]);
+	for(i=0;i<NUMVERTICES;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]);
 
 	/*Now, we have to move the previous Vz inputs to old 
@@ -5630,14 +5606,13 @@
 	int i;
 
-	const int    numvertices=6;
 	const int    numdofpervertex=4;
-	const int    numdof=numdofpervertex*numvertices;
+	const int    numdof=numdofpervertex*NUMVERTICES;
 	int*         doflist=NULL;
 	double       values[numdof];
-	double       vx[numvertices];
-	double       vy[numvertices];
-	double       vz[numvertices];
-	double       vel[numvertices];
-	double       pressure[numvertices];
+	double       vx[NUMVERTICES];
+	double       vy[NUMVERTICES];
+	double       vz[NUMVERTICES];
+	double       vel[NUMVERTICES];
+	double       pressure[NUMVERTICES];
 	double       stokesreconditioning;
 
@@ -5651,5 +5626,5 @@
 
 	/*Ok, we have vx and vy in values, fill in all arrays: */
-	for(i=0;i<numvertices;i++){
+	for(i=0;i<NUMVERTICES;i++){
 		vx[i]=values[i*numdofpervertex+0];
 		vy[i]=values[i*numdofpervertex+1];
@@ -5660,10 +5635,10 @@
 	/*Recondition pressure: */
 	this->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum);
-	for(i=0;i<numvertices;i++){
+	for(i=0;i<NUMVERTICES;i++){
 		pressure[i]=pressure[i]*stokesreconditioning;
 	}
 
 	/*Now Compute vel*/
-	for(i=0;i<numvertices;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
+	for(i=0;i<NUMVERTICES;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
 	
 	/*Now, we have to move the previous inputs  to old 
@@ -5691,13 +5666,12 @@
 	int i;
 
-	const int    numvertices=6;
 	const int    numdofpervertex=4;
-	const int    numdof=numdofpervertex*numvertices;
+	const int    numdof=numdofpervertex*NUMVERTICES;
 	int*         doflist=NULL;
 	double       values[numdof];
-	double       lambdax[numvertices];
-	double       lambday[numvertices];
-	double       lambdaz[numvertices];
-	double       lambdap[numvertices];
+	double       lambdax[NUMVERTICES];
+	double       lambday[NUMVERTICES];
+	double       lambdaz[NUMVERTICES];
+	double       lambdap[NUMVERTICES];
 
 	/*Get dof list: */
@@ -5710,5 +5684,5 @@
 
 	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
-	for(i=0;i<numvertices;i++){
+	for(i=0;i<NUMVERTICES;i++){
 		lambdax[i]=values[i*numdofpervertex+0];
 		lambday[i]=values[i*numdofpervertex+1];
@@ -5732,11 +5706,10 @@
 	int i;
 
-	const int    numvertices=6;
 	const int    numdofpervertex=2;
-	const int    numdof=numdofpervertex*numvertices;
+	const int    numdof=numdofpervertex*NUMVERTICES;
 	int*         doflist=NULL;
 	double       values[numdof];
-	double       lambdax[numvertices];
-	double       lambday[numvertices];
+	double       lambdax[NUMVERTICES];
+	double       lambday[NUMVERTICES];
 
 	/*Get dof list: */
@@ -5749,5 +5722,5 @@
 
 	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
-	for(i=0;i<numvertices;i++){
+	for(i=0;i<NUMVERTICES;i++){
 		lambdax[i]=values[i*numdofpervertex+0];
 		lambday[i]=values[i*numdofpervertex+1];
@@ -5767,7 +5740,6 @@
 	int i;
 
-	const int    numvertices=6;
 	const int    numdofpervertex=1;
-	const int    numdof=numdofpervertex*numvertices;
+	const int    numdof=numdofpervertex*NUMVERTICES;
 	int*         doflist=NULL;
 	double       values[numdof];
@@ -5801,7 +5773,6 @@
 void  Penta::InputUpdateFromSolutionOneDof(double* solution,int enum_type){
 
-	const int numvertices     = 6;
 	const int numdofpervertex = 1;
-	const int numdof          = numdofpervertex *numvertices;
+	const int numdof          = numdofpervertex *NUMVERTICES;
 	int*         doflist=NULL;
 	double    values[numdof];
@@ -5825,7 +5796,6 @@
 void  Penta::InputUpdateFromSolutionOneDofCollapsed(double* solution,int enum_type){
 
-	const int  numvertices     = 6;
 	const int  numdofpervertex = 1;
-	const int  numdof          = numdofpervertex *numvertices;
+	const int  numdof          = numdofpervertex *NUMVERTICES;
 	const int  numdof2d        = numdof/2;
 	int*         doflist=NULL;
Index: /issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 5661)
+++ /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 5662)
@@ -18,4 +18,7 @@
 #include "../../include/include.h"
 /*}}}*/
+
+/*Element macros*/
+#define NUMVERTICES 3
 
 /*Tria constructors and destructor*/
@@ -449,5 +452,4 @@
 void  Tria::InputUpdateFromVectorDakota(double* vector, int name, int type){
 	
-	const int numvertices=3;
 	int i,j;
 
@@ -535,9 +537,8 @@
 
 	int       i,j;
-	const int numvertices=3;
 	double    area;
 	double    mean;
-	int       partition[numvertices];
-	int       offset[numvertices];
+	int       partition[NUMVERTICES];
+	int       offset[NUMVERTICES];
 	double    values[3];
 	bool      already=false;
@@ -549,11 +550,11 @@
 	this->GetDofList1(&offset[0]);
 	mean=0;
-	for(i=0;i<numvertices;i++){
+	for(i=0;i<NUMVERTICES;i++){
 		partition[i]=(int)qmu_part[offset[i]];
-		mean=mean+1.0/numvertices*vertex_response[offset[i]];
+		mean=mean+1.0/NUMVERTICES*vertex_response[offset[i]];
 	}
 
 	/*Add contribution: */
-	for(i=0;i<numvertices;i++){
+	for(i=0;i<NUMVERTICES;i++){
 		already=false;
 		for(j=0;j<i;j++){
@@ -627,7 +628,6 @@
 
 	/*constants: */
-	const int    numvertices=3;
 	const int    NDOF2=2;
-	const int    numdof=NDOF2*numvertices;
+	const int    numdof=NDOF2*NUMVERTICES;
 
 	/* Intermediaries */
@@ -638,5 +638,5 @@
 	double     cm_noisedmp;
 	double     Jdet;
-	double     xyz_list[numvertices][3];
+	double     xyz_list[NUMVERTICES][3];
 	double     dk[NDOF2],dB[NDOF2];
 	GaussTria *gauss = NULL;
@@ -654,5 +654,5 @@
 	if(onwater) return 0;
 
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 
 	/* Start looping on the number of gaussian points: */
@@ -857,6 +857,5 @@
 void  Tria::GetVectorFromInputs(Vec vector,int NameEnum){
 
-	const int numvertices=3;
-	int doflist1[numvertices];
+	int doflist1[NUMVERTICES];
 
 	/*Find NameEnum input in the inputs dataset, and get it to fill in the vector: */
@@ -901,5 +900,4 @@
 
 	/* constants*/
-	const int    numvertices=3;
 	const int    NDOF2=2;
 
@@ -910,11 +908,11 @@
 	double     vx,vy,lambda,mu,thickness;
 	double     cm_noisedmp;
-	int        doflist[numvertices];
+	int        doflist[NUMVERTICES];
 	double     dvx[NDOF2],dvy[NDOF2],dadjx[NDOF2],dadjy[NDOF2],dB[NDOF2]; 
-	double     xyz_list[numvertices][3];
+	double     xyz_list[NUMVERTICES][3];
 	double     basis[3];
-	double     dbasis[NDOF2][numvertices];
-	double     grad[numvertices]={0.0};
-	double     grad_g[numvertices];
+	double     dbasis[NDOF2][NUMVERTICES];
+	double     grad[NUMVERTICES]={0.0};
+	double     grad_g[NUMVERTICES];
 	double     epsilon[3]; /* epsilon=[exx,eyy,exy];*/
 	GaussTria *gauss = NULL;
@@ -924,5 +922,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	GetDofList1(&doflist[0]);
 
@@ -956,17 +954,17 @@
 
 		/*standard gradient dJ/dki*/
-		for (i=0;i<numvertices;i++){
+		for (i=0;i<NUMVERTICES;i++){
 			grad_g[i]=-viscosity_complement*thickness*( (2*dvx[0]+dvy[1])*2*dadjx[0]+(dvx[1]+dvy[0])*(dadjx[1]+dadjy[0])+(2*dvy[1]+dvx[0])*2*dadjy[1])*Jdet*gauss->weight*basis[i]; 
 		}
 		/*Add regularization term*/
-		for (i=0;i<numvertices;i++){
+		for (i=0;i<NUMVERTICES;i++){
 			grad_g[i]-=cm_noisedmp*Jdet*gauss->weight*(dbasis[0][i]*dB[0]+dbasis[1][i]*dB[1]);
 		}
 
-		for(i=0;i<numvertices;i++) grad[i]+=grad_g[i];
+		for(i=0;i<NUMVERTICES;i++) grad[i]+=grad_g[i];
 	}
 
 	/*Add grade_g to global vector gradient: */
-	VecSetValues(gradient,numvertices,doflist,(const double*)grad,ADD_VALUES);
+	VecSetValues(gradient,NUMVERTICES,doflist,(const double*)grad,ADD_VALUES);
 
 	/*clean-up*/
@@ -980,9 +978,8 @@
 
 	/* node data: */
-	const int    numvertices=3;
 	const int    NDOF2=2;
-	double       xyz_list[numvertices][3];
-	int          doflist1[numvertices];
-	double       dh1dh3[NDOF2][numvertices];
+	double       xyz_list[NUMVERTICES][3];
+	int          doflist1[NUMVERTICES];
+	double       dh1dh3[NDOF2][NUMVERTICES];
 
 	/* gaussian points: */
@@ -1001,6 +998,6 @@
 
 	/*element vector at the gaussian points: */
-	double  grade_g[numvertices]={0.0};
-	double  grade_g_gaussian[numvertices];
+	double  grade_g[NUMVERTICES]={0.0};
+	double  grade_g_gaussian[NUMVERTICES];
 
 	/* Jacobian: */
@@ -1039,5 +1036,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	GetDofList1(&doflist1[0]);
 
@@ -1087,5 +1084,5 @@
 
 		/*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
-		for (i=0;i<numvertices;i++){
+		for (i=0;i<NUMVERTICES;i++){
 
 			//standard term dJ/dki
@@ -1097,9 +1094,9 @@
 		
 		/*Add gradje_g_gaussian vector to gradje_g: */
-		for( i=0; i<numvertices; i++)grade_g[i]+=grade_g_gaussian[i];
+		for( i=0; i<NUMVERTICES; i++)grade_g[i]+=grade_g_gaussian[i];
 	}
 
 	/*Add grade_g to global vector gradient: */
-	VecSetValues(gradient,numvertices,doflist1,(const double*)grade_g,ADD_VALUES);
+	VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES);
 
 	/*Clean up and return*/
@@ -1115,9 +1112,8 @@
 
 	/* node data: */
-	const int    numvertices=3;
-	int          doflist1[numvertices];
+	int          doflist1[NUMVERTICES];
 
 	/*Gauss*/
-	double  gauss[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}};
+	double  gauss[NUMVERTICES][NUMVERTICES]={{1,0,0},{0,1,0},{0,0,1}};
 
 	/* grid data: */
@@ -1125,5 +1121,5 @@
 
 	/*element vector at the gaussian points: */
-	double  gradient_g[numvertices];
+	double  gradient_g[NUMVERTICES];
 
 	/*Inputs*/
@@ -1137,5 +1133,5 @@
 
 	/* Start  looping on the vertices: */
-	for(i=0; i<numvertices;i++){
+	for(i=0; i<NUMVERTICES;i++){
 		adjoint_input->GetParameterValue(&lambda,&gauss[i][0]);
 		gradient_g[i]=-lambda;
@@ -1143,5 +1139,5 @@
 
 	/*Add grade_g to global vector gradient: */
-	VecSetValues(gradient,numvertices,doflist1,(const double*)gradient_g,INSERT_VALUES);
+	VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)gradient_g,INSERT_VALUES);
 }
 /*}}}*/
@@ -1340,8 +1336,7 @@
 	int i;
 
-	const int    numvertices=3;
 	const int    numdofs=2;
 	double mass_flux=0;
-	double xyz_list[numvertices][3];
+	double xyz_list[NUMVERTICES][3];
 	double gauss_1[3];
 	double gauss_2[3];
@@ -1363,5 +1358,5 @@
 
 	/*Get xyz list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 
 	/*get area coordinates of 0 and 1 locations: */
@@ -1552,6 +1547,5 @@
 
 	/* Constants */
-	const int    numvertices=3;
-	const int    numdof=1*numvertices;
+	const int    numdof=1*NUMVERTICES;
 
 	/*Intermediaries*/
@@ -1561,5 +1555,5 @@
 	double     Jdet;
 	double     Jelem = 0;
-	double     xyz_list[numvertices][3];
+	double     xyz_list[NUMVERTICES][3];
 	GaussTria *gauss = NULL;
 
@@ -1568,5 +1562,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 
 	/*Retrieve all inputs we will be needing: */
@@ -1608,17 +1602,16 @@
 
 	/* node data: */
-	const int    numvertices=3;
-	const int    numdof=2*numvertices;
+	const int    numdof=2*NUMVERTICES;
 	const int    NDOF2=2;
-	double       xyz_list[numvertices][3];
+	double       xyz_list[NUMVERTICES][3];
 
 	/* grid data: */
-	double vx_list[numvertices];
-	double vy_list[numvertices];
-	double obs_vx_list[numvertices];
-	double obs_vy_list[numvertices];
-	double misfit_square_list[numvertices];
-	double misfit_list[numvertices];
-	double weights_list[numvertices];
+	double vx_list[NUMVERTICES];
+	double vy_list[NUMVERTICES];
+	double obs_vx_list[NUMVERTICES];
+	double obs_vy_list[NUMVERTICES];
+	double misfit_square_list[NUMVERTICES];
+	double misfit_list[NUMVERTICES];
+	double weights_list[NUMVERTICES];
 
 	/* gaussian points: */
@@ -1645,5 +1638,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 
 	/* Recover input data: */
@@ -1675,12 +1668,12 @@
 	 *
 	 */
-	for (i=0;i<numvertices;i++){
+	for (i=0;i<NUMVERTICES;i++){
 		misfit_list[i]=0.5*(pow((vx_list[i]-obs_vx_list[i]),(double)2)+pow((vy_list[i]-obs_vy_list[i]),(double)2));
 	}
 	/*Process units: */
-	if(process_units)UnitConversion(&misfit_list[0],numvertices,IuToExtEnum,SurfaceAbsVelMisfitEnum,this->parameters);
+	if(process_units)UnitConversion(&misfit_list[0],NUMVERTICES,IuToExtEnum,SurfaceAbsVelMisfitEnum,this->parameters);
 
 	/*Apply weights to misfits*/
-	for (i=0;i<numvertices;i++){
+	for (i=0;i<NUMVERTICES;i++){
 		misfit_list[i]=weights_list[i]*misfit_list[i];
 	}
@@ -1711,17 +1704,16 @@
 
 	/* node data: */
-	const int    numvertices=3;
-	const int    numdof=2*numvertices;
+	const int    numdof=2*NUMVERTICES;
 	const int    NDOF2=2;
-	double       xyz_list[numvertices][3];
+	double       xyz_list[NUMVERTICES][3];
 
 	/* grid data: */
-	double vx_list[numvertices];
-	double vy_list[numvertices];
-	double obs_vx_list[numvertices];
-	double obs_vy_list[numvertices];
-	double misfit_square_list[numvertices];
-	double misfit_list[numvertices];
-	double weights_list[numvertices];
+	double vx_list[NUMVERTICES];
+	double vy_list[NUMVERTICES];
+	double obs_vx_list[NUMVERTICES];
+	double obs_vy_list[NUMVERTICES];
+	double misfit_square_list[NUMVERTICES];
+	double misfit_list[NUMVERTICES];
+	double weights_list[NUMVERTICES];
 
 	/* gaussian points: */
@@ -1752,5 +1744,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 
 	/* Recover input data: */
@@ -1782,5 +1774,5 @@
 	 *              obs                        obs                      
 	 */
-	for (i=0;i<numvertices;i++){
+	for (i=0;i<NUMVERTICES;i++){
 		scalex=pow(meanvel/(obs_vx_list[i]+epsvel),(double)2);
 		scaley=pow(meanvel/(obs_vy_list[i]+epsvel),(double)2);
@@ -1791,8 +1783,8 @@
 
 	/*Process units: */
-	if(process_units)UnitConversion(&misfit_list[0],numvertices,IuToExtEnum,SurfaceRelVelMisfitEnum,this->parameters);
+	if(process_units)UnitConversion(&misfit_list[0],NUMVERTICES,IuToExtEnum,SurfaceRelVelMisfitEnum,this->parameters);
 
 	/*Apply weights to misfits*/
-	for (i=0;i<numvertices;i++){
+	for (i=0;i<NUMVERTICES;i++){
 		misfit_list[i]=weights_list[i]*misfit_list[i];
 	}
@@ -1823,17 +1815,16 @@
 
 	/* node data: */
-	const int    numvertices=3;
-	const int    numdof=2*numvertices;
+	const int    numdof=2*NUMVERTICES;
 	const int    NDOF2=2;
-	double       xyz_list[numvertices][3];
+	double       xyz_list[NUMVERTICES][3];
 
 	/* grid data: */
-	double vx_list[numvertices];
-	double vy_list[numvertices];
-	double obs_vx_list[numvertices];
-	double obs_vy_list[numvertices];
-	double misfit_square_list[numvertices];
-	double misfit_list[numvertices];
-	double weights_list[numvertices];
+	double vx_list[NUMVERTICES];
+	double vy_list[NUMVERTICES];
+	double obs_vx_list[NUMVERTICES];
+	double obs_vy_list[NUMVERTICES];
+	double misfit_square_list[NUMVERTICES];
+	double misfit_list[NUMVERTICES];
+	double weights_list[NUMVERTICES];
 
 	/* gaussian points: */
@@ -1863,5 +1854,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 
 	/* Recover input data: */
@@ -1893,5 +1884,5 @@
 	 *                            obs
 	 */
-	for (i=0;i<numvertices;i++){
+	for (i=0;i<NUMVERTICES;i++){
 		velocity_mag=sqrt(pow(vx_list[i],(double)2)+pow(vy_list[i],(double)2))+epsvel; //epsvel to avoid velocity being nil.
 		obs_velocity_mag=sqrt(pow(obs_vx_list[i],(double)2)+pow(obs_vy_list[i],(double)2))+epsvel; //epsvel to avoid observed velocity being nil.
@@ -1900,8 +1891,8 @@
 
 	/*Process units: */
-	if(process_units)UnitConversion(&misfit_list[0],numvertices,IuToExtEnum,SurfaceLogVelMisfitEnum,this->parameters);
+	if(process_units)UnitConversion(&misfit_list[0],NUMVERTICES,IuToExtEnum,SurfaceLogVelMisfitEnum,this->parameters);
 
 	/*Apply weights to misfits*/
-	for (i=0;i<numvertices;i++){
+	for (i=0;i<NUMVERTICES;i++){
 		misfit_list[i]=weights_list[i]*misfit_list[i];
 	}
@@ -1932,17 +1923,16 @@
 
 	/* node data: */
-	const int    numvertices=3;
-	const int    numdof=2*numvertices;
+	const int    numdof=2*NUMVERTICES;
 	const int    NDOF2=2;
-	double       xyz_list[numvertices][3];
+	double       xyz_list[NUMVERTICES][3];
 
 	/* grid data: */
-	double vx_list[numvertices];
-	double vy_list[numvertices];
-	double obs_vx_list[numvertices];
-	double obs_vy_list[numvertices];
-	double misfit_square_list[numvertices];
-	double misfit_list[numvertices];
-	double weights_list[numvertices];
+	double vx_list[NUMVERTICES];
+	double vy_list[NUMVERTICES];
+	double obs_vx_list[NUMVERTICES];
+	double obs_vy_list[NUMVERTICES];
+	double misfit_square_list[NUMVERTICES];
+	double misfit_list[NUMVERTICES];
+	double weights_list[NUMVERTICES];
 
 	/* gaussian points: */
@@ -1975,5 +1965,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 
 	/* Recover input data: */
@@ -2005,5 +1995,5 @@
 	 *                              obs                       obs
 	 */
-	for (i=0;i<numvertices;i++){
+	for (i=0;i<NUMVERTICES;i++){
 		misfit_list[i]=0.5*pow(meanvel,(double)2)*(
 					pow(log((fabs(vx_list[i])+epsvel)/(fabs(obs_vx_list[i])+epsvel)),(double)2) +
@@ -2012,8 +2002,8 @@
 
 	/*Process units: */
-	if(process_units)UnitConversion(&misfit_list[0],numvertices,IuToExtEnum,SurfaceLogVxVyMisfitEnum,this->parameters);
+	if(process_units)UnitConversion(&misfit_list[0],NUMVERTICES,IuToExtEnum,SurfaceLogVxVyMisfitEnum,this->parameters);
 
 	/*Apply weights to misfits*/
-	for (i=0;i<numvertices;i++){
+	for (i=0;i<NUMVERTICES;i++){
 		misfit_list[i]=weights_list[i]*misfit_list[i];
 	}
@@ -2044,17 +2034,16 @@
 
 	/* node data: */
-	const int    numvertices=3;
-	const int    numdof=2*numvertices;
+	const int    numdof=2*NUMVERTICES;
 	const int    NDOF2=2;
-	double       xyz_list[numvertices][3];
+	double       xyz_list[NUMVERTICES][3];
 
 	/* grid data: */
-	double vx_list[numvertices];
-	double vy_list[numvertices];
-	double obs_vx_list[numvertices];
-	double obs_vy_list[numvertices];
-	double misfit_square_list[numvertices];
-	double misfit_list[numvertices];
-	double weights_list[numvertices];
+	double vx_list[NUMVERTICES];
+	double vy_list[NUMVERTICES];
+	double obs_vx_list[NUMVERTICES];
+	double obs_vy_list[NUMVERTICES];
+	double misfit_square_list[NUMVERTICES];
+	double misfit_list[NUMVERTICES];
+	double weights_list[NUMVERTICES];
 
 	/* gaussian points: */
@@ -2087,5 +2076,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 
 	/* Recover input data: */
@@ -2117,14 +2106,14 @@
 	 *      S                obs            obs
 	 */
-	for (i=0;i<numvertices;i++) misfit_square_list[i]=pow(vx_list[i]-obs_vx_list[i],2)+pow(vy_list[i]-obs_vx_list[i],2);
+	for (i=0;i<NUMVERTICES;i++) misfit_square_list[i]=pow(vx_list[i]-obs_vx_list[i],2)+pow(vy_list[i]-obs_vx_list[i],2);
 
 	/*Process units: */
-	if(process_units)UnitConversion(&misfit_square_list[0],numvertices,IuToExtEnum,SurfaceAverageVelMisfitEnum,this->parameters);
+	if(process_units)UnitConversion(&misfit_square_list[0],NUMVERTICES,IuToExtEnum,SurfaceAverageVelMisfitEnum,this->parameters);
 
 	/*Take the square root, and scale by surface: */
-	for (i=0;i<numvertices;i++)misfit_list[i]=pow(misfit_square_list[i],2)/S;
+	for (i=0;i<NUMVERTICES;i++)misfit_list[i]=pow(misfit_square_list[i],2)/S;
 
 	/*Apply weights to misfits*/
-	for (i=0;i<numvertices;i++){
+	for (i=0;i<NUMVERTICES;i++){
 		misfit_list[i]=weights_list[i]*misfit_list[i];
 	}
@@ -2182,5 +2171,4 @@
 	/*output: */
 	int     numrows     = 0;
-	int     numvertices = 0;
 	int     numnodes    = 0;
 
@@ -2191,5 +2179,4 @@
 		numrows++;
 		/*now, how many vertices and how many nodal values for this result? :*/
-		numvertices=3; //this is a tria element, with 3 vertices
 		numnodes=elementresult->NumberOfNodalValues(); //ask result object.
 	}
@@ -2197,5 +2184,5 @@
 	/*Assign output pointers:*/
 	*pnumrows=numrows;
-	*pnumvertices=numvertices;
+	*pnumvertices=NUMVERTICES;
 	*pnumnodes=numnodes;
 	
@@ -2223,6 +2210,5 @@
 
 	/* node data: */
-	int numvertices=3;
-	double xyz_list[numvertices][3];
+	double xyz_list[NUMVERTICES][3];
 	double v13[3];
 	double v23[3];
@@ -2239,5 +2225,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 
 	for (i=0;i<3;i++){
@@ -2489,8 +2475,7 @@
 
 	/* node data: */
-	const int    numvertices=3;
 	const int    NDOF1=1;
-	const int    numdof=NDOF1*numvertices;
-	double       xyz_list[numvertices][3];
+	const int    numdof=NDOF1*NUMVERTICES;
+	double       xyz_list[NUMVERTICES][3];
 	int*         doflist=NULL;
 
@@ -2500,7 +2485,7 @@
 
 	/* matrices: */
-	double L[numvertices];
-	double B[2][numvertices];
-	double Bprime[2][numvertices];
+	double L[NUMVERTICES];
+	double B[2][NUMVERTICES];
+	double Bprime[2][NUMVERTICES];
 	double DL[2][2]={0.0};
 	double DLprime[2][2]={0.0};
@@ -2531,5 +2516,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	GetDofList(&doflist);
 
@@ -2645,8 +2630,7 @@
 
 	/* node data: */
-	const int    numvertices=3;
 	const int    NDOF1=1;
-	const int    numdof=NDOF1*numvertices;
-	double       xyz_list[numvertices][3];
+	const int    numdof=NDOF1*NUMVERTICES;
+	double       xyz_list[NUMVERTICES][3];
 	int*         doflist=NULL;
 
@@ -2656,6 +2640,6 @@
 
 	/* matrices: */
-	double B[2][numvertices];
-	double Bprime[2][numvertices];
+	double B[2][NUMVERTICES];
+	double Bprime[2][NUMVERTICES];
 	double DL[2][2]={0.0};
 	double DLprime[2][2]={0.0};
@@ -2675,5 +2659,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	GetDofList(&doflist);
 	this->parameters->FindParam(&dim,DimEnum);
@@ -2734,8 +2718,7 @@
 
 	/* node data: */
-	const int    numvertices=3;
 	const int    NDOF1=1;
-	const int    numdof=NDOF1*numvertices;
-	double       xyz_list[numvertices][3];
+	const int    numdof=NDOF1*NUMVERTICES;
+	double       xyz_list[NUMVERTICES][3];
 	int*         doflist=NULL;
 
@@ -2745,7 +2728,7 @@
 
 	/* matrices: */
-	double L[numvertices];
-	double B[2][numvertices];
-	double Bprime[2][numvertices];
+	double L[NUMVERTICES];
+	double B[2][NUMVERTICES];
+	double Bprime[2][NUMVERTICES];
 	double DL[2][2]={0.0};
 	double DLprime[2][2]={0.0};
@@ -2782,5 +2765,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	GetDofList(&doflist);
 
@@ -2798,5 +2781,5 @@
 	/*Modify z so that it reflects the surface*/
 	GetParameterListOnVertices(&surface_list[0],SurfaceEnum);
-	for(i=0;i<numvertices;i++) xyz_list[i][2]=surface_list[i];
+	for(i=0;i<NUMVERTICES;i++) xyz_list[i][2]=surface_list[i];
 
 	/*Get normal vector to the surface*/
@@ -2903,7 +2886,6 @@
 
 	/* node data: */
-	const int    numvertices=3;
-	const int    numdof=2*numvertices;
-	double       xyz_list[numvertices][3];
+	const int    numdof=2*NUMVERTICES;
+	double       xyz_list[NUMVERTICES][3];
 	int*         doflist=NULL;
 
@@ -2958,5 +2940,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	GetDofList(&doflist,MacAyealApproximationEnum);
 
@@ -3034,7 +3016,6 @@
 
 	/* node data: */
-	const int numvertices = 3;
-	const int numdof   = 2 *numvertices;
-	double    xyz_list[numvertices][3];
+	const int numdof   = 2 *NUMVERTICES;
+	double    xyz_list[NUMVERTICES][3];
 	int*      doflistm=NULL;
 	int*      doflistp=NULL;
@@ -3087,5 +3068,5 @@
 	
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	GetDofList(&doflistm,MacAyealApproximationEnum);
 	GetDofList(&doflistp,PattynApproximationEnum);
@@ -3159,7 +3140,6 @@
 
 	/* node data: */
-	const int numvertices = 3;
-	const int numdof   = 2 *numvertices;
-	double    xyz_list[numvertices][3];
+	const int numdof   = 2 *NUMVERTICES;
+	double    xyz_list[NUMVERTICES][3];
 	int*      doflist=NULL;
 	int       numberofdofspernode=2;
@@ -3211,5 +3191,5 @@
 	
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	GetDofList(&doflist,MacAyealApproximationEnum);
 
@@ -3280,7 +3260,6 @@
 
 	/* node data: */
-	const int numvertices = 3;
-	const int numdof   = 2 *numvertices;
-	double    xyz_list[numvertices][3];
+	const int numdof   = 2 *NUMVERTICES;
+	double    xyz_list[NUMVERTICES][3];
 	int*      doflist=NULL;
 	int       numberofdofspernode=2;
@@ -3332,5 +3311,5 @@
 	
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	GetDofList(&doflist,PattynApproximationEnum);
 
@@ -3399,7 +3378,6 @@
 	int    i;
 	int    connectivity;
-	const int numvertices=3;
 	const int NDOF2=2;
-	const int numdofs=numvertices*NDOF2;
+	const int numdofs=NUMVERTICES*NDOF2;
 	int*         doflist=NULL;
 
@@ -3436,8 +3414,7 @@
 
 	/* node data: */
-	const int    numvertices=3;
 	const int    NDOF1=1;
-	const int    numdof=NDOF1*numvertices;
-	double       xyz_list[numvertices][3];
+	const int    numdof=NDOF1*NUMVERTICES;
+	double       xyz_list[NUMVERTICES][3];
 	int*         doflist=NULL;
 
@@ -3466,5 +3443,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	GetDofList(&doflist);
 
@@ -3539,11 +3516,10 @@
 	int i,j;
 
-	const int  numvertices=3;
 	const int  NDOF1=1;
-	const int  numdof=numvertices*NDOF1;
+	const int  numdof=NUMVERTICES*NDOF1;
 	int*       doflist=NULL;
 
 	/*Grid data: */
-	double     xyz_list[numvertices][3];
+	double     xyz_list[NUMVERTICES][3];
 
 	/*Material constants */
@@ -3567,5 +3543,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	GetDofList(&doflist);
 
@@ -3589,6 +3565,6 @@
 		MatrixMultiply(&tLD[0],numdof,1,0,&L[0],1,numdof,0,&Ke_gaussian[0][0],0);
 
-		for(i=0;i<numvertices;i++){
-			for(j=0;j<numvertices;j++){
+		for(i=0;i<NUMVERTICES;i++){
+			for(j=0;j<NUMVERTICES;j++){
 				K_terms[i][j]+=Ke_gaussian[i][j];
 			}
@@ -3611,8 +3587,7 @@
 
 	/* node data: */
-	const int    numvertices=3;
 	const int    NDOF1=1;
-	const int    numdof=NDOF1*numvertices;
-	double       xyz_list[numvertices][3];
+	const int    numdof=NDOF1*NUMVERTICES;
+	double       xyz_list[NUMVERTICES][3];
 	int*         doflist=NULL;
 	int          numberofdofspernode=1;
@@ -3623,7 +3598,7 @@
 
 	/* matrices: */
-	double L[numvertices];
-	double B[2][numvertices];
-	double Bprime[2][numvertices];
+	double L[NUMVERTICES];
+	double B[2][NUMVERTICES];
+	double Bprime[2][NUMVERTICES];
 	double DL[2][2]={0.0};
 	double DLprime[2][2]={0.0};
@@ -3658,5 +3633,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	GetDofList(&doflist);
 
@@ -3781,8 +3756,7 @@
 
 	/* node data: */
-	const int    numvertices=3;
 	const int    NDOF1=1;
-	const int    numdof=NDOF1*numvertices;
-	double       xyz_list[numvertices][3];
+	const int    numdof=NDOF1*NUMVERTICES;
+	double       xyz_list[NUMVERTICES][3];
 	int*         doflist=NULL;
 	int          numberofdofspernode=1;
@@ -3793,7 +3767,7 @@
 
 	/* matrices: */
-	double L[numvertices];
-	double B[2][numvertices];
-	double Bprime[2][numvertices];
+	double L[NUMVERTICES];
+	double B[2][NUMVERTICES];
+	double Bprime[2][NUMVERTICES];
 	double DL[2][2]={0.0};
 	double DLprime[2][2]={0.0};
@@ -3819,5 +3793,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	GetDofList(&doflist);
 
@@ -3941,8 +3915,7 @@
 	
 	/* node data: */
-	const int    numvertices=3;
 	const int    NDOF1=1;
-	const int    numdof=NDOF1*numvertices;
-	double       xyz_list[numvertices][3];
+	const int    numdof=NDOF1*NUMVERTICES;
+	double       xyz_list[NUMVERTICES][3];
 	int*         doflist=NULL;
 
@@ -3960,5 +3933,5 @@
 	double  K_terms[numdof][numdof]={0.0};
 	double  Ke_gaussian[numdof][numdof]={0.0};
-	double  l1l2l3[numvertices];
+	double  l1l2l3[NUMVERTICES];
 	double     tl1l2l3D[3];
 	double  D_scalar;
@@ -3971,5 +3944,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	GetDofList(&doflist);
 
@@ -4026,8 +3999,7 @@
 
 	/* node data: */
-	const int    numvertices=3;
 	const int    NDOF1=1;
-	const int    numdof=NDOF1*numvertices;
-	double       xyz_list[numvertices][3];
+	const int    numdof=NDOF1*NUMVERTICES;
+	double       xyz_list[NUMVERTICES][3];
 	int*         doflist=NULL;
 	int          numberofdofspernode=1;
@@ -4038,6 +4010,6 @@
 
 	/* matrix */
-	double pe_g[numvertices]={0.0};
-	double L[numvertices];
+	double pe_g[NUMVERTICES]={0.0};
+	double L[NUMVERTICES];
 	double Jdettria;
 
@@ -4053,5 +4025,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	GetDofList(&doflist);
 
@@ -4098,8 +4070,7 @@
 
 	/* node data: */
-	const int    numvertices=3;
 	const int    NDOF1=1;
-	const int    numdof=NDOF1*numvertices;
-	double       xyz_list[numvertices][3];
+	const int    numdof=NDOF1*NUMVERTICES;
+	double       xyz_list[NUMVERTICES][3];
 	int*         doflist=NULL;
 	int          numberofdofspernode=1;
@@ -4110,6 +4081,6 @@
 
 	/* matrix */
-	double pe_g[numvertices]={0.0};
-	double L[numvertices];
+	double pe_g[NUMVERTICES]={0.0};
+	double L[NUMVERTICES];
 	double Jdettria;
 
@@ -4125,5 +4096,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	GetDofList(&doflist);
 
@@ -4170,8 +4141,7 @@
 
 	/* node data: */
-	const int    numvertices=3;
 	const int    NDOF1=1;
-	const int    numdof=NDOF1*numvertices;
-	double       xyz_list[numvertices][3];
+	const int    numdof=NDOF1*NUMVERTICES;
+	double       xyz_list[NUMVERTICES][3];
 	int*         doflist=NULL;
 	int          numberofdofspernode=1;
@@ -4182,6 +4152,6 @@
 
 	/* matrix */
-	double pe_g[numvertices]={0.0};
-	double L[numvertices];
+	double pe_g[NUMVERTICES]={0.0};
+	double L[NUMVERTICES];
 	double Jdettria;
 
@@ -4195,5 +4165,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	GetDofList(&doflist);
 
@@ -4237,8 +4207,7 @@
 
 	/* node data: */
-	const int    numvertices=3;
 	const int    NDOF1=1;
-	const int    numdof=NDOF1*numvertices;
-	double       xyz_list[numvertices][3];
+	const int    numdof=NDOF1*NUMVERTICES;
+	double       xyz_list[NUMVERTICES][3];
 	int*         doflist=NULL;
 
@@ -4258,5 +4227,5 @@
 
 	/* matrices: */
-	double L[numvertices];
+	double L[NUMVERTICES];
 
 	/*input parameters for structural analysis (diagnostic): */
@@ -4273,5 +4242,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	GetDofList(&doflist);
 
@@ -4308,5 +4277,5 @@
 
 		/*Build gaussian vector: */
-		for(i=0;i<numvertices;i++){
+		for(i=0;i<NUMVERTICES;i++){
 			pe_g_gaussian[i]=-Jdet*gauss->weight*(vx*dbdx+vy*dbdy-meltingvalue)*L[i];
 		}
@@ -4331,8 +4300,7 @@
 
 	/* node data: */
-	const int    numvertices=3;
-	const int    numdof=2*numvertices;
+	const int    numdof=2*NUMVERTICES;
 	const int    NDOF2=2;
-	double       xyz_list[numvertices][3];
+	double       xyz_list[NUMVERTICES][3];
 	int*         doflist=NULL;
 	
@@ -4382,5 +4350,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	GetDofList(&doflist,MacAyealApproximationEnum);
 
@@ -4418,5 +4386,5 @@
 		/*Build pe_g_gaussian vector: */
 		if(drag_type==1){
-			for (i=0;i<numvertices;i++){
+			for (i=0;i<NUMVERTICES;i++){
 				for (j=0;j<NDOF2;j++){
 					pe_g_gaussian[i*NDOF2+j]=(-driving_stress_baseline*slope[j]-plastic_stress)*Jdet*gauss_weight*l1l2l3[i]; 
@@ -4425,5 +4393,5 @@
 		}
 		else {
-			for (i=0;i<numvertices;i++){
+			for (i=0;i<NUMVERTICES;i++){
 				for (j=0;j<NDOF2;j++){
 					pe_g_gaussian[i*NDOF2+j]=-driving_stress_baseline*slope[j]*Jdet*gauss_weight*l1l2l3[i];
@@ -4454,7 +4422,6 @@
 
 	/* node data: */
-	const int    numvertices=3;
-	const int    numdof=1*numvertices;
-	double       xyz_list[numvertices][3];
+	const int    numdof=1*NUMVERTICES;
+	double       xyz_list[NUMVERTICES][3];
 	int*         doflist=NULL;
 
@@ -4482,5 +4449,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	GetDofList(&doflist);
 
@@ -4507,5 +4474,5 @@
 		weights_input->GetParameterValue(&weight, gauss);
 
-		for (i=0;i<numvertices;i++){
+		for (i=0;i<NUMVERTICES;i++){
 			pe_g_gaussian[i]=(thicknessobs-thickness)*weight*Jdet*gauss->weight*l1l2l3[i]; 
 		}
@@ -4531,18 +4498,17 @@
 
 	/* node data: */
-	const int    numvertices=3;
-	const int    numdof=2*numvertices;
+	const int    numdof=2*NUMVERTICES;
 	const int    NDOF2=2;
-	double       xyz_list[numvertices][3];
+	double       xyz_list[NUMVERTICES][3];
 	int*         doflist=NULL;
 
 	/* grid data: */
-	double vx_list[numvertices];
-	double vy_list[numvertices];
-	double obs_vx_list[numvertices];
-	double obs_vy_list[numvertices];
-	double dux_list[numvertices];
-	double duy_list[numvertices];
-	double weights_list[numvertices];
+	double vx_list[NUMVERTICES];
+	double vy_list[NUMVERTICES];
+	double obs_vx_list[NUMVERTICES];
+	double obs_vy_list[NUMVERTICES];
+	double dux_list[NUMVERTICES];
+	double duy_list[NUMVERTICES];
+	double weights_list[NUMVERTICES];
 
 	/* gaussian points: */
@@ -4577,5 +4543,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	GetDofList(&doflist);
 
@@ -4616,5 +4582,5 @@
 		 *        du     obs
 		 */
-		for (i=0;i<numvertices;i++){
+		for (i=0;i<NUMVERTICES;i++){
 			dux_list[i]=obs_vx_list[i]-vx_list[i];
 			duy_list[i]=obs_vy_list[i]-vy_list[i];
@@ -4634,5 +4600,5 @@
 		 *               obs
 		 */
-		for (i=0;i<numvertices;i++){
+		for (i=0;i<NUMVERTICES;i++){
 			scalex=pow(meanvel/(obs_vx_list[i]+epsvel),2);
 			scaley=pow(meanvel/(obs_vy_list[i]+epsvel),2);
@@ -4656,5 +4622,5 @@
 		 *            
 		 */
-		for (i=0;i<numvertices;i++){
+		for (i=0;i<NUMVERTICES;i++){
 			velocity_mag=sqrt(pow(vx_list[i],2)+pow(vy_list[i],2))+epsvel; //epsvel to avoid velocity being nil.
 			obs_velocity_mag=sqrt(pow(obs_vx_list[i],2)+pow(obs_vy_list[i],2))+epsvel; //epsvel to avoid observed velocity being nil.
@@ -4675,5 +4641,5 @@
 		 *        du      S  2 sqrt(...)           obs
 		 */
-		for (i=0;i<numvertices;i++){
+		for (i=0;i<NUMVERTICES;i++){
 			scale=1.0/(S*sqrt(pow(vx_list[i]-obs_vx_list[i],2)+pow(vy_list[i]-obs_vx_list[i],2))+epsvel);
 			dux_list[i]=scale*(obs_vx_list[i]-vx_list[i]);
@@ -4692,5 +4658,5 @@
 		 *        du                         |u| + eps  |u|                           u + eps
 		 */
-		for (i=0;i<numvertices;i++){
+		for (i=0;i<NUMVERTICES;i++){
 			dux_list[i] = - pow(meanvel,(double)2)*(
 						log((fabs(vx_list[i])+epsvel)/(fabs(obs_vx_list[i])+epsvel)) * 1/(vx_list[i]+epsvel));
@@ -4705,5 +4671,5 @@
 
 	/*Apply weights to DU*/
-	for (i=0;i<numvertices;i++){
+	for (i=0;i<NUMVERTICES;i++){
 		dux_list[i]=weights_list[i]*dux_list[i];
 		duy_list[i]=weights_list[i]*duy_list[i];
@@ -4729,5 +4695,5 @@
 
 		/*compute Du*/
-		for (i=0;i<numvertices;i++){
+		for (i=0;i<NUMVERTICES;i++){
 			pe_g_gaussian[i*NDOF2+0]=dux*Jdet*gauss->weight*l1l2l3[i]; 
 			pe_g_gaussian[i*NDOF2+1]=duy*Jdet*gauss->weight*l1l2l3[i]; 
@@ -4754,18 +4720,17 @@
 
 	/* node data: */
-	const int    numvertices=3;
 	const int    NDOF4=4;
-	const int    numdof=NDOF4*numvertices;
-	double       xyz_list[numvertices][3];
+	const int    numdof=NDOF4*NUMVERTICES;
+	double       xyz_list[NUMVERTICES][3];
 	int*         doflist=NULL;
 
 	/* grid data: */
-	double vx_list[numvertices];
-	double vy_list[numvertices];
-	double obs_vx_list[numvertices];
-	double obs_vy_list[numvertices];
-	double dux_list[numvertices];
-	double duy_list[numvertices];
-	double weights_list[numvertices];
+	double vx_list[NUMVERTICES];
+	double vy_list[NUMVERTICES];
+	double obs_vx_list[NUMVERTICES];
+	double obs_vy_list[NUMVERTICES];
+	double dux_list[NUMVERTICES];
+	double duy_list[NUMVERTICES];
+	double weights_list[NUMVERTICES];
 
 	/* gaussian points: */
@@ -4800,5 +4765,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	GetDofList(&doflist);
 
@@ -4839,5 +4804,5 @@
 		 *        du     obs
 		 */
-		for (i=0;i<numvertices;i++){
+		for (i=0;i<NUMVERTICES;i++){
 			dux_list[i]=obs_vx_list[i]-vx_list[i];
 			duy_list[i]=obs_vy_list[i]-vy_list[i];
@@ -4857,5 +4822,5 @@
 		 *               obs
 		 */
-		for (i=0;i<numvertices;i++){
+		for (i=0;i<NUMVERTICES;i++){
 			scalex=pow(meanvel/(obs_vx_list[i]+epsvel),2);
 			scaley=pow(meanvel/(obs_vy_list[i]+epsvel),2);
@@ -4879,5 +4844,5 @@
 		 *            
 		 */
-		for (i=0;i<numvertices;i++){
+		for (i=0;i<NUMVERTICES;i++){
 			velocity_mag=sqrt(pow(vx_list[i],2)+pow(vy_list[i],2))+epsvel; //epsvel to avoid velocity being nil.
 			obs_velocity_mag=sqrt(pow(obs_vx_list[i],2)+pow(obs_vy_list[i],2))+epsvel; //epsvel to avoid observed velocity being nil.
@@ -4898,5 +4863,5 @@
 		 *        du      S  2 sqrt(...)           obs
 		 */
-		for (i=0;i<numvertices;i++){
+		for (i=0;i<NUMVERTICES;i++){
 			scale=1.0/(S*sqrt(pow(vx_list[i]-obs_vx_list[i],2)+pow(vy_list[i]-obs_vx_list[i],2))+epsvel);
 			dux_list[i]=scale*(obs_vx_list[i]-vx_list[i]);
@@ -4915,5 +4880,5 @@
 		 *        du                         |u| + eps  |u|                           u + eps
 		 */
-		for (i=0;i<numvertices;i++){
+		for (i=0;i<NUMVERTICES;i++){
 			dux_list[i] = - pow(meanvel,(double)2)*(
 						log((fabs(vx_list[i])+epsvel)/(fabs(obs_vx_list[i])+epsvel)) * 1/(vx_list[i]+epsvel));
@@ -4928,5 +4893,5 @@
 
 	/*Apply weights to DU*/
-	for (i=0;i<numvertices;i++){
+	for (i=0;i<NUMVERTICES;i++){
 		dux_list[i]=weights_list[i]*dux_list[i];
 		duy_list[i]=weights_list[i]*duy_list[i];
@@ -4952,5 +4917,5 @@
 
 		/*compute Du*/
-		for (i=0;i<numvertices;i++){
+		for (i=0;i<NUMVERTICES;i++){
 			pe_g[i*NDOF4+0]+=dux*Jdet*gauss->weight*l1l2l3[i];
 			pe_g[i*NDOF4+1]+=duy*Jdet*gauss->weight*l1l2l3[i]; 
@@ -4971,7 +4936,6 @@
 	/*Collapsed formulation: */
 	int       i;
-	const int numvertices=3;
 	const int NDOF2=2;
-	const int numdofs=NDOF2*numvertices;
+	const int numdofs=NDOF2*NUMVERTICES;
 	int*         doflist=NULL;
 	double    constant_part,ub,vb;
@@ -5011,5 +4975,5 @@
 	/*Spawn 3 sing elements: */
 	gauss=new GaussTria();
-	for(i=0;i<numvertices;i++){
+	for(i=0;i<NUMVERTICES;i++){
 
 		gauss->GaussVertex(i);
@@ -5048,8 +5012,7 @@
 
 	/* node data: */
-	const int    numvertices=3;
 	const int    NDOF1=1;
-	const int    numdof=NDOF1*numvertices;
-	double       xyz_list[numvertices][3];
+	const int    numdof=NDOF1*NUMVERTICES;
+	double       xyz_list[NUMVERTICES][3];
 	int*         doflist=NULL;
 	int          numberofdofspernode=1;
@@ -5060,6 +5023,6 @@
 
 	/* matrix */
-	double pe_g[numvertices]={0.0};
-	double L[numvertices];
+	double pe_g[NUMVERTICES]={0.0};
+	double L[NUMVERTICES];
 	double Jdettria;
 
@@ -5079,5 +5042,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	GetDofList(&doflist);
 
@@ -5125,8 +5088,7 @@
 
 	/* node data: */
-	const int    numvertices=3;
 	const int    NDOF1=1;
-	const int    numdof=NDOF1*numvertices;
-	double       xyz_list[numvertices][3];
+	const int    numdof=NDOF1*NUMVERTICES;
+	double       xyz_list[NUMVERTICES][3];
 	int*         doflist=NULL;
 	int          numberofdofspernode=1;
@@ -5137,6 +5099,6 @@
 
 	/* matrix */
-	double pe_g[numvertices]={0.0};
-	double L[numvertices];
+	double pe_g[NUMVERTICES]={0.0};
+	double L[NUMVERTICES];
 	double Jdettria;
 
@@ -5154,5 +5116,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	GetDofList(&doflist);
 
@@ -5198,8 +5160,7 @@
 
 	/* node data: */
-	const int    numvertices=3;
 	const int    NDOF1=1;
-	const int    numdof=NDOF1*numvertices;
-	double       xyz_list[numvertices][3];
+	const int    numdof=NDOF1*NUMVERTICES;
+	double       xyz_list[NUMVERTICES][3];
 	int*         doflist=NULL;
 	
@@ -5227,5 +5188,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	GetDofList(&doflist);
 
@@ -5278,9 +5239,8 @@
 	int i,found;
 	
-	const int  numvertices=3;
 	const int  NDOF1=1;
-	const int  numdof=numvertices*NDOF1;
+	const int  numdof=NUMVERTICES*NDOF1;
 	int*         doflist=NULL;
-	double       xyz_list[numvertices][3];
+	double       xyz_list[NUMVERTICES][3];
 
 	double mixed_layer_capacity;
@@ -5303,5 +5263,5 @@
 	double  Jdet;
 	double  P_terms[numdof]={0.0};
-	double  l1l2l3[numvertices];
+	double  l1l2l3[NUMVERTICES];
 
 	double  t_pmp;
@@ -5312,5 +5272,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	GetDofList(&doflist);
 
@@ -5372,9 +5332,8 @@
 	int i,found;
 	
-	const int  numvertices=3;
 	const int  NDOF1=1;
-	const int  numdof=numvertices*NDOF1;
+	const int  numdof=NUMVERTICES*NDOF1;
 	int*       doflist=NULL;
-	double     xyz_list[numvertices][3];
+	double     xyz_list[NUMVERTICES][3];
 
 	double rho_ice;
@@ -5398,5 +5357,5 @@
 	double  Jdet;
 	double  P_terms[numdof]={0.0};
-	double  l1l2l3[numvertices];
+	double  l1l2l3[NUMVERTICES];
 	double  scalar;
 
@@ -5413,5 +5372,5 @@
 	
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	GetDofList(&doflist);
 
@@ -5479,10 +5438,9 @@
 
 	double area=0;
-	const int    numvertices=3;
-	double xyz_list[numvertices][3];
+	double xyz_list[NUMVERTICES][3];
 	double x1,y1,x2,y2,x3,y3;
 
 	/*Get xyz list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	x1=xyz_list[0][0]; y1=xyz_list[0][1];
 	x2=xyz_list[1][0]; y2=xyz_list[1][1];
@@ -5497,6 +5455,5 @@
 	/*Intermediaries*/
 	double    area = 0;
-	const int numvertices = 3;
-	double    xyz_list[numvertices][3];
+	double    xyz_list[NUMVERTICES][3];
 	double    x1,y1,x2,y2,x3,y3;
 
@@ -5505,5 +5462,5 @@
 
 	/*Get xyz list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	x1=xyz_list[0][0]; y1=xyz_list[0][1];
 	x2=xyz_list[1][0]; y2=xyz_list[1][1];
@@ -5577,6 +5534,5 @@
 
 	/*Intermediaries*/
-	const int  numvertices        = 3;
-	double     value[numvertices];
+	double     value[NUMVERTICES];
 	GaussTria *gauss              = NULL;
 
@@ -5590,5 +5546,5 @@
 	/* Start looping on the number of vertices: */
 	gauss=new GaussTria();
-	for (int iv=0;iv<numvertices;iv++){
+	for (int iv=0;iv<NUMVERTICES;iv++){
 		gauss->GaussVertex(iv);
 		input->GetParameterValue(&pvalue[iv],gauss);
@@ -5603,6 +5559,5 @@
 
 	/*Intermediaries*/
-	const int  numvertices        = 3;
-	double     value[numvertices];
+	double     value[NUMVERTICES];
 	GaussTria *gauss              = NULL;
 
@@ -5616,5 +5571,5 @@
 	if (input){
 		gauss=new GaussTria();
-		for (int iv=0;iv<numvertices;iv++){
+		for (int iv=0;iv<NUMVERTICES;iv++){
 			gauss->GaussVertex(iv);
 			input->GetParameterValue(&pvalue[iv],gauss);
@@ -5622,5 +5577,5 @@
 	}
 	else{
-		for (int iv=0;iv<numvertices;iv++) pvalue[iv]=defaultvalue;
+		for (int iv=0;iv<NUMVERTICES;iv++) pvalue[iv]=defaultvalue;
 	}
 
@@ -5786,7 +5741,6 @@
 
 	int i;
-	const int    numvertices=3;
 	const int    numdofpervertex=2;
-	const int    numdof=numdofpervertex*numvertices;
+	const int    numdof=numdofpervertex*NUMVERTICES;
 	int*         doflist=NULL;
 	double       values[numdof];
@@ -5805,5 +5759,5 @@
 	/*P1 element only for now*/
 	gauss=new GaussTria();
-	for(i=0;i<numvertices;i++){
+	for(i=0;i<NUMVERTICES;i++){
 
 		gauss->GaussVertex(i);
@@ -5829,7 +5783,6 @@
 
 	int i;
-	const int    numvertices=3;
 	const int    numdofpervertex=2;
-	const int    numdof=numdofpervertex*numvertices;
+	const int    numdof=numdofpervertex*NUMVERTICES;
 	int*         doflist=NULL;
 	double       values[numdof];
@@ -5849,5 +5802,5 @@
 	/*P1 element only for now*/
 	gauss=new GaussTria();
-	for(i=0;i<numvertices;i++){
+	for(i=0;i<NUMVERTICES;i++){
 
 		gauss->GaussVertex(i);
@@ -5927,9 +5880,8 @@
 
 	/* node data: */
-	const int    numvertices=3;
 	const int    NDOF2=2;
-	double       xyz_list[numvertices][3];
-	int          doflist1[numvertices];
-	double       dh1dh3[NDOF2][numvertices];
+	double       xyz_list[NUMVERTICES][3];
+	int          doflist1[NUMVERTICES];
+	double       dh1dh3[NDOF2][NUMVERTICES];
 
 	/* grid data: */
@@ -5956,6 +5908,6 @@
 
 	/*element vector at the gaussian points: */
-	double  grade_g[numvertices]={0.0};
-	double  grade_g_gaussian[numvertices];
+	double  grade_g[NUMVERTICES]={0.0};
+	double  grade_g_gaussian[NUMVERTICES];
 
 	/* Jacobian: */
@@ -5991,5 +5943,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	GetDofList1(&doflist1[0]);
 
@@ -6044,5 +5996,5 @@
 
 		/*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
-		for (i=0;i<numvertices;i++){
+		for (i=0;i<NUMVERTICES;i++){
 			//standard gradient dJ/dki
 			grade_g_gaussian[i]=(
@@ -6057,9 +6009,9 @@
 
 		/*Add gradje_g_gaussian vector to gradje_g: */
-		for( i=0; i<numvertices; i++)grade_g[i]+=grade_g_gaussian[i];
+		for( i=0; i<NUMVERTICES; i++)grade_g[i]+=grade_g_gaussian[i];
 	}
 
 	/*Add grade_g to global vector gradient: */
-	VecSetValues(gradient,numvertices,doflist1,(const double*)grade_g,ADD_VALUES);
+	VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES);
 
 	/*Add grade_g to the inputs of this element: */
@@ -6079,10 +6031,9 @@
 	int i;
 
-	const int    numvertices=3;
 	const int    numdofpervertex=1;
-	const int    numdof=numdofpervertex*numvertices;
+	const int    numdof=numdofpervertex*NUMVERTICES;
 	int*         doflist=NULL;
 	double       values[numdof];
-	double       lambda[numvertices];
+	double       lambda[NUMVERTICES];
 
 	/*Get dof list: */
@@ -6112,11 +6063,10 @@
 	int i;
 
-	const int    numvertices=3;
 	const int    numdofpervertex=2;
-	const int    numdof=numdofpervertex*numvertices;
+	const int    numdof=numdofpervertex*NUMVERTICES;
 	int*         doflist=NULL;
 	double       values[numdof];
-	double       lambdax[numvertices];
-	double       lambday[numvertices];
+	double       lambdax[NUMVERTICES];
+	double       lambday[NUMVERTICES];
 
 	/*Get dof list: */
@@ -6129,5 +6079,5 @@
 
 	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
-	for(i=0;i<numvertices;i++){
+	for(i=0;i<NUMVERTICES;i++){
 		lambdax[i]=values[i*numdofpervertex+0];
 		lambday[i]=values[i*numdofpervertex+1];
@@ -6148,15 +6098,14 @@
 	int i;
 
-	const int    numvertices=3;
 	const int    numdofpervertex=2;
-	const int    numdof=numdofpervertex*numvertices;
+	const int    numdof=numdofpervertex*NUMVERTICES;
 	int*         doflist=NULL;
 	double       values[numdof];
-	double       vx[numvertices];
-	double       vy[numvertices];
-	double       vz[numvertices];
-	double       vel[numvertices];
-	double       pressure[numvertices];
-	double       thickness[numvertices];
+	double       vx[NUMVERTICES];
+	double       vy[NUMVERTICES];
+	double       vz[NUMVERTICES];
+	double       vel[NUMVERTICES];
+	double       pressure[NUMVERTICES];
+	double       thickness[NUMVERTICES];
 	double       rho_ice,g;
 	int          dummy;
@@ -6171,5 +6120,5 @@
 
 	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
-	for(i=0;i<numvertices;i++){
+	for(i=0;i<NUMVERTICES;i++){
 		vx[i]=values[i*numdofpervertex+0];
 		vy[i]=values[i*numdofpervertex+1];
@@ -6180,5 +6129,5 @@
 
 	/*Now Compute vel*/
-	for(i=0;i<numvertices;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
+	for(i=0;i<NUMVERTICES;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
 
 	/*For pressure: we have not computed pressure in this analysis, for this element. We are in 2D, 
@@ -6187,5 +6136,5 @@
 	g=matpar->GetG();
 	GetParameterListOnVertices(&thickness[0],ThicknessEnum);
-	for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*thickness[i];
+	for(i=0;i<NUMVERTICES;i++) pressure[i]=rho_ice*g*thickness[i];
 
 	/*Now, we have to move the previous Vx and Vy inputs  to old 
@@ -6211,15 +6160,14 @@
 	int i;
 
-	const int    numvertices=3;
 	const int    numdofpervertex=2;
-	const int    numdof=numdofpervertex*numvertices;
+	const int    numdof=numdofpervertex*NUMVERTICES;
 	int*         doflist=NULL;
 	double       values[numdof];
-	double       vx[numvertices];
-	double       vy[numvertices];
-	double       vz[numvertices];
-	double       vel[numvertices];
-	double       pressure[numvertices];
-	double       thickness[numvertices];
+	double       vx[NUMVERTICES];
+	double       vy[NUMVERTICES];
+	double       vz[NUMVERTICES];
+	double       vel[NUMVERTICES];
+	double       pressure[NUMVERTICES];
+	double       thickness[NUMVERTICES];
 	double       rho_ice,g;
 	Input*       vz_input=NULL;
@@ -6236,5 +6184,5 @@
 
 	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
-	for(i=0;i<numvertices;i++){
+	for(i=0;i<NUMVERTICES;i++){
 		vx[i]=values[i*numdofpervertex+0];
 		vy[i]=values[i*numdofpervertex+1];
@@ -6248,12 +6196,12 @@
 		}
 		vz_input->GetValuesPtr(&vz_ptr,&dummy);
-		for(i=0;i<numvertices;i++) vz[i]=vz_ptr[i];
+		for(i=0;i<NUMVERTICES;i++) vz[i]=vz_ptr[i];
 	}
 	else{
-		for(i=0;i<numvertices;i++) vz[i]=0.0;
+		for(i=0;i<NUMVERTICES;i++) vz[i]=0.0;
 	}
 
 	/*Now Compute vel*/
-	for(i=0;i<numvertices;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
+	for(i=0;i<NUMVERTICES;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
 
 	/*For pressure: we have not computed pressure in this analysis, for this element. We are in 2D, 
@@ -6263,5 +6211,5 @@
 	GetParameterListOnVertices(&thickness[0],ThicknessEnum);
 	
-	for(i=0;i<numvertices;i++){
+	for(i=0;i<NUMVERTICES;i++){
 		pressure[i]=rho_ice*g*thickness[i];
 	}
@@ -6287,7 +6235,6 @@
 void  Tria::InputUpdateFromSolutionOneDof(double* solution,int enum_type){
 
-	const int numvertices     = 3;
 	const int numdofpervertex = 1;
-	const int numdof          = numdofpervertex *numvertices;
+	const int numdof          = numdofpervertex *NUMVERTICES;
 	int*         doflist=NULL;
 	double    values[numdof];
Index: /issm/trunk/src/c/objects/Elements/Tria.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.h	(revision 5661)
+++ /issm/trunk/src/c/objects/Elements/Tria.h	(revision 5662)
@@ -22,7 +22,4 @@
 #include "../../EnumDefinitions/EnumDefinitions.h"
 /*}}}*/
-
-/*Element macros*/
-#define NUMVERTICES 3
 
 class Tria: public Element,public TriaHook,public TriaRef{