Index: /issm/trunk/src/c/objects/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Penta.cpp	(revision 3159)
+++ /issm/trunk/src/c/objects/Penta.cpp	(revision 3160)
@@ -3049,6 +3049,4 @@
 	double z1,z2,z3,z4,z5,z6;
 
-	double sqrt3=SQRT3;
-
 	/*Figure out xi,eta and zi (parametric coordinates), for this gaussian point: */
 	A1=gauss_coord[0];
@@ -3057,5 +3055,5 @@
 
 	xi=A2-A1;
-	eta=sqrt3*A3;
+	eta=SQRT3*A3;
 	zi=gauss_coord[3];
 
@@ -3082,15 +3080,15 @@
 
 
-	*(J+NDOF3*0+0)=1.0/4.0*(x1-x2-x4+x5)*zi+1.0/4.0*(-x1+x2-x4+x5);
-	*(J+NDOF3*1+0)=sqrt3/12.0*(x1+x2-2*x3-x4-x5+2*x6)*zi+sqrt3/12.0*(-x1-x2+2*x3-x4-x5+2*x6);
-	*(J+NDOF3*2+0)=sqrt3/12.0*(x1+x2-2*x3-x4-x5+2*x6)*eta+1/4*(x1-x2-x4+x5)*xi +1.0/4.0*(-x1+x5-x2+x4);
-
-	*(J+NDOF3*0+1)=1.0/4.0*(y1-y2-y4+y5)*zi+1.0/4.0*(-y1+y2-y4+y5);
-	*(J+NDOF3*1+1)=sqrt3/12.0*(y1+y2-2*y3-y4-y5+2*y6)*zi+sqrt3/12.0*(-y1-y2+2*y3-y4-y5+2*y6);
-	*(J+NDOF3*2+1)=sqrt3/12.0*(y1+y2-2*y3-y4-y5+2*y6)*eta+1.0/4.0*(y1-y2-y4+y5)*xi+1.0/4.0*(y4-y1+y5-y2);
-
-	*(J+NDOF3*0+2)=1.0/4.0*(z1-z2-z4+z5)*zi+1.0/4.0*(-z1+z2-z4+z5);
-	*(J+NDOF3*1+2)=sqrt3/12.0*(z1+z2-2*z3-z4-z5+2*z6)*zi+sqrt3/12.0*(-z1-z2+2*z3-z4-z5+2*z6);
-	*(J+NDOF3*2+2)=sqrt3/12.0*(z1+z2-2*z3-z4-z5+2*z6)*eta+1.0/4.0*(z1-z2-z4+z5)*xi+1.0/4.0*(-z1+z5-z2+z4);
+	*(J+NDOF3*0+0)=0.25*(x1-x2-x4+x5)*zi+0.25*(-x1+x2-x4+x5);
+	*(J+NDOF3*1+0)=SQRT3/12.0*(x1+x2-2*x3-x4-x5+2*x6)*zi+SQRT3/12.0*(-x1-x2+2*x3-x4-x5+2*x6);
+	*(J+NDOF3*2+0)=SQRT3/12.0*(x1+x2-2*x3-x4-x5+2*x6)*eta+1/4*(x1-x2-x4+x5)*xi +0.25*(-x1+x5-x2+x4);
+
+	*(J+NDOF3*0+1)=0.25*(y1-y2-y4+y5)*zi+0.25*(-y1+y2-y4+y5);
+	*(J+NDOF3*1+1)=SQRT3/12.0*(y1+y2-2*y3-y4-y5+2*y6)*zi+SQRT3/12.0*(-y1-y2+2*y3-y4-y5+2*y6);
+	*(J+NDOF3*2+1)=SQRT3/12.0*(y1+y2-2*y3-y4-y5+2*y6)*eta+0.25*(y1-y2-y4+y5)*xi+0.25*(y4-y1+y5-y2);
+
+	*(J+NDOF3*0+2)=0.25*(z1-z2-z4+z5)*zi+0.25*(-z1+z2-z4+z5);
+	*(J+NDOF3*1+2)=SQRT3/12.0*(z1+z2-2*z3-z4-z5+2*z6)*zi+SQRT3/12.0*(-z1-z2+2*z3-z4-z5+2*z6);
+	*(J+NDOF3*2+2)=SQRT3/12.0*(z1+z2-2*z3-z4-z5+2*z6)*eta+0.25*(z1-z2-z4+z5)*xi+0.25*(-z1+z5-z2+z4);
 
 #ifdef _ISSM_DEBUG_
@@ -3461,5 +3459,4 @@
 	const int numgrids=6;
 	double A1,A2,A3,z;
-	double sqrt3=SQRT3;
 
 	A1=gauss_coord[0]; //first area coordinate value. In term of xi and eta: A1=(1-xi)/2-eta/(2*SQRT3);
@@ -3470,32 +3467,32 @@
 
 	/*First nodal function derivatives. The corresponding nodal function is N=A1*(1-z)/2. Its derivatives follow*/
-	*(dl1dl6+numgrids*0+0)=-1.0/2.0*(1.0-z)/2.0;
-	*(dl1dl6+numgrids*1+0)=-1.0/2.0/sqrt3*(1.0-z)/2.0;
-	*(dl1dl6+numgrids*2+0)=-1.0/2.0*A1;
+	*(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;
 
 	/*Second nodal function: The corresponding nodal function is N=A2*(1-z)/2. Its derivatives follow*/
-	*(dl1dl6+numgrids*0+1)=1.0/2.0*(1.0-z)/2.0;
-	*(dl1dl6+numgrids*1+1)=-1.0/2.0/sqrt3*(1.0-z)/2.0;
-	*(dl1dl6+numgrids*2+1)=-1.0/2.0*A2;
+	*(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;
 
 	/*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)=-1.0/2.0*A3;
+	*(dl1dl6+numgrids*1+2)=1.0/SQRT3*(1.0-z)/2.0;
+	*(dl1dl6+numgrids*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)=-1.0/2.0*(1.0+z)/2.0;
-	*(dl1dl6+numgrids*1+3)=-1.0/2.0/sqrt3*(1.0+z)/2.0;
-	*(dl1dl6+numgrids*2+3)=1.0/2.0*A1;
+	*(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;
 
 	/*Fifth nodal function: The corresponding nodal function is N=A2*(1+z)/2. Its derivatives follow*/
-	*(dl1dl6+numgrids*0+4)=1.0/2.0*(1.0+z)/2.0;
-	*(dl1dl6+numgrids*1+4)=-1.0/2.0/sqrt3*(1.0+z)/2.0;
-	*(dl1dl6+numgrids*2+4)=1.0/2.0*A2;
+	*(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;
 
 	/*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)=1.0/2.0*A3;
+	*(dl1dl6+numgrids*1+5)=1.0/SQRT3*(1.0+z)/2.0;
+	*(dl1dl6+numgrids*2+5)=0.5*A3;
 }
 /*}}}*/
@@ -3508,44 +3505,43 @@
 	 * natural coordinate system) at the gaussian point. */
 
-	double sqrt3=SQRT3;
 	int    numgrids=7; //six plus bubble grids
 
 	double r=gauss_coord[1]-gauss_coord[0];
-	double s=-3.0/sqrt3*(gauss_coord[0]+gauss_coord[1]-2.0/3.0);
+	double s=-3.0/SQRT3*(gauss_coord[0]+gauss_coord[1]-2.0/3.0);
 	double zeta=gauss_coord[3];
 
 	/*First nodal function: */
-	*(dl1dl7+numgrids*0+0)=-1.0/2.0*(1.0-zeta)/2.0;
-	*(dl1dl7+numgrids*1+0)=-sqrt3/6.0*(1.0-zeta)/2.0;
-	*(dl1dl7+numgrids*2+0)=-1.0/2.0*(-1.0/2.0*r-sqrt3/6.0*s+1.0/3.0);
+	*(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);
 
 	/*Second nodal function: */
-	*(dl1dl7+numgrids*0+1)=1.0/2.0*(1.0-zeta)/2.0;
-	*(dl1dl7+numgrids*1+1)=-sqrt3/6.0*(1.0-zeta)/2.0;
-	*(dl1dl7+numgrids*2+1)=-1.0/2.0*(1.0/2.0*r-sqrt3/6.0*s+1.0/3.0);
+	*(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);
 
 	/*Third nodal function: */
 	*(dl1dl7+numgrids*0+2)=0;
-	*(dl1dl7+numgrids*1+2)=sqrt3/3.0*(1.0-zeta)/2.0;
-	*(dl1dl7+numgrids*2+2)=-1.0/2.0*(sqrt3/3.0*s+1.0/3.0);
+	*(dl1dl7+numgrids*1+2)=SQRT3/3.0*(1.0-zeta)/2.0;
+	*(dl1dl7+numgrids*2+2)=-0.5*(SQRT3/3.0*s+ONETHIRD);
 
 	/*Fourth nodal function: */
-	*(dl1dl7+numgrids*0+3)=-1.0/2.0*(1.0+zeta)/2.0;
-	*(dl1dl7+numgrids*1+3)=-sqrt3/6.0*(1.0+zeta)/2.0;
-	*(dl1dl7+numgrids*2+3)=1.0/2.0*(-1.0/2.0*r-sqrt3/6.0*s+1.0/3.0);
+	*(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);
 
 	/*Fith nodal function: */
-	*(dl1dl7+numgrids*0+4)=1.0/2.0*(1.0+zeta)/2.0;
-	*(dl1dl7+numgrids*1+4)=-sqrt3/6.0*(1.0+zeta)/2.0;
-	*(dl1dl7+numgrids*2+4)=1.0/2.0*(1.0/2.0*r-sqrt3/6.0*s+1.0/3.0);
+	*(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);
 
 	/*Sixth nodal function: */
 	*(dl1dl7+numgrids*0+5)=0;
-	*(dl1dl7+numgrids*1+5)=sqrt3/3.0*(1.0+zeta)/2.0;
-	*(dl1dl7+numgrids*2+5)=1.0/2.0*(sqrt3/3.0*s+1.0/3.0);
+	*(dl1dl7+numgrids*1+5)=SQRT3/3.0*(1.0+zeta)/2.0;
+	*(dl1dl7+numgrids*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*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_coord[0]*gauss_coord[1]*gauss_coord[2]*(-2.0*zeta);
 
@@ -3711,5 +3707,5 @@
 	epsilon_sqr[2][2]=pow(epsilon_matrix[2][2],2);
 
-	epsilon_eff=1/pow(2,1.0/2.0)*pow((epsilon_sqr[0][0]+epsilon_sqr[0][1]+ epsilon_sqr[0][2]+ epsilon_sqr[1][0]+ epsilon_sqr[1][1]+ epsilon_sqr[1][2]+ epsilon_sqr[2][0]+ epsilon_sqr[2][1]+ epsilon_sqr[2][2]),1.0/2.0);
+	epsilon_eff=1/pow(2,0.5)*pow((epsilon_sqr[0][0]+epsilon_sqr[0][1]+ epsilon_sqr[0][2]+ epsilon_sqr[1][0]+ epsilon_sqr[1][1]+ epsilon_sqr[1][2]+ epsilon_sqr[2][0]+ epsilon_sqr[2][1]+ epsilon_sqr[2][2]),0.5);
 	*phi=2*pow(epsilon_eff,2.0)*viscosity;
 }
Index: /issm/trunk/src/c/objects/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Tria.cpp	(revision 3159)
+++ /issm/trunk/src/c/objects/Tria.cpp	(revision 3160)
@@ -388,10 +388,10 @@
 	if(numpar->artdiff){
 		//Get the Jacobian determinant
-		gauss_l1l2l3[0]=1.0/3.0; gauss_l1l2l3[1]=1.0/3.0; gauss_l1l2l3[2]=1.0/3.0;
+		gauss_l1l2l3[0]=ONETHIRD; gauss_l1l2l3[1]=ONETHIRD; gauss_l1l2l3[2]=ONETHIRD;
 		GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
 
 		//Build K matrix (artificial diffusivity matrix)
-		v_gauss[0]=1.0/3.0*(vxvy_list[0][0]+vxvy_list[1][0]+vxvy_list[2][0]);
-		v_gauss[1]=1.0/3.0*(vxvy_list[0][1]+vxvy_list[1][1]+vxvy_list[2][1]);
+		v_gauss[0]=ONETHIRD*(vxvy_list[0][0]+vxvy_list[1][0]+vxvy_list[2][0]);
+		v_gauss[1]=ONETHIRD*(vxvy_list[0][1]+vxvy_list[1][1]+vxvy_list[2][1]);
 
 		K[0][0]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]);
@@ -611,10 +611,10 @@
 	if(numpar->artdiff){
 		//Get the Jacobian determinant
-		gauss_l1l2l3[0]=1.0/3.0; gauss_l1l2l3[1]=1.0/3.0; gauss_l1l2l3[2]=1.0/3.0;
+		gauss_l1l2l3[0]=ONETHIRD; gauss_l1l2l3[1]=ONETHIRD; gauss_l1l2l3[2]=ONETHIRD;
 		GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
 
 		//Build K matrix (artificial diffusivity matrix)
-		v_gauss[0]=1.0/3.0*(vxvy_list[0][0]+vxvy_list[1][0]+vxvy_list[2][0]);
-		v_gauss[1]=1.0/3.0*(vxvy_list[0][1]+vxvy_list[1][1]+vxvy_list[2][1]);
+		v_gauss[0]=ONETHIRD*(vxvy_list[0][0]+vxvy_list[1][0]+vxvy_list[2][0]);
+		v_gauss[1]=ONETHIRD*(vxvy_list[0][1]+vxvy_list[1][1]+vxvy_list[2][1]);
 
 		K[0][0]=pow(10,2)*pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]); //pow should be zero!!
@@ -1408,10 +1408,10 @@
 	if(numpar->artdiff==1){
 		//Get the Jacobian determinant
-		gauss_l1l2l3[0]=1.0/3.0; gauss_l1l2l3[1]=1.0/3.0; gauss_l1l2l3[2]=1.0/3.0;
+		gauss_l1l2l3[0]=ONETHIRD; gauss_l1l2l3[1]=ONETHIRD; gauss_l1l2l3[2]=ONETHIRD;
 		GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
 
 		//Build K matrix (artificial diffusivity matrix)
-		v_gauss[0]=1.0/3.0*(vxvy_list[0][0]+vxvy_list[1][0]+vxvy_list[2][0]);
-		v_gauss[1]=1.0/3.0*(vxvy_list[0][1]+vxvy_list[1][1]+vxvy_list[2][1]);
+		v_gauss[0]=ONETHIRD*(vxvy_list[0][0]+vxvy_list[1][0]+vxvy_list[2][0]);
+		v_gauss[1]=ONETHIRD*(vxvy_list[0][1]+vxvy_list[1][1]+vxvy_list[2][1]);
 
 		K[0][0]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]);
@@ -3262,7 +3262,7 @@
 
 
-	*(J+NDOF2*0+0)=1.0/2.0*(x2-x1);
+	*(J+NDOF2*0+0)=0.5*(x2-x1);
 	*(J+NDOF2*1+0)=SQRT3/6.0*(2*x3-x1-x2);
-	*(J+NDOF2*0+1)=1.0/2.0*(y2-y1);
+	*(J+NDOF2*0+1)=0.5*(y2-y1);
 	*(J+NDOF2*1+1)=SQRT3/6.0*(2*y3-y1-y2);
 }
@@ -3468,9 +3468,9 @@
 
 	/*First nodal function: */
-	*(dl1dl3+numgrids*0+0)=-1.0/2.0; 
+	*(dl1dl3+numgrids*0+0)=-0.5; 
 	*(dl1dl3+numgrids*1+0)=-1.0/(2.0*SQRT3);
 
 	/*Second nodal function: */
-	*(dl1dl3+numgrids*0+1)=1.0/2.0;
+	*(dl1dl3+numgrids*0+1)=0.5;
 	*(dl1dl3+numgrids*1+1)=-1.0/(2.0*SQRT3);
 
@@ -4291,6 +4291,6 @@
 
 	mass_flux= rho_ice*length*(  
-				(1.0/3.0*(h1-h2)*(vx1-vx2)+1.0/2.0*h2*(vx1-vx2)+1.0/2.0*(h1-h2)*vx2+h2*vx2)*normal[0]+
-				(1.0/3.0*(h1-h2)*(vy1-vy2)+1.0/2.0*h2*(vy1-vy2)+1.0/2.0*(h1-h2)*vy2+h2*vy2)*normal[1]
+				(ONETHIRD*(h1-h2)*(vx1-vx2)+0.5*h2*(vx1-vx2)+0.5*(h1-h2)*vx2+h2*vx2)*normal[0]+
+				(ONETHIRD*(h1-h2)*(vy1-vy2)+0.5*h2*(vy1-vy2)+0.5*(h1-h2)*vy2+h2*vy2)*normal[1]
 				);
 	return mass_flux;