Index: /issm/trunk/src/c/Makefile.am
===================================================================
--- /issm/trunk/src/c/Makefile.am	(revision 3158)
+++ /issm/trunk/src/c/Makefile.am	(revision 3159)
@@ -129,5 +129,5 @@
 					./shared/Elements/ResolvePointers.cpp\
 					./shared/Elements/Paterson.cpp\
-					./shared/Elements/GetElementNodeData.cpp\
+					./shared/Elements/GetVerticesCoordinates.cpp\
 					./shared/String/isdistributed.cpp\
 					./shared/String/sharedstring.h\
@@ -456,5 +456,5 @@
 					./shared/Elements/ResolvePointers.cpp\
 					./shared/Elements/Paterson.cpp\
-					./shared/Elements/GetElementNodeData.cpp\
+					./shared/Elements/GetVerticesCoordinates.cpp\
 					./shared/String/isdistributed.cpp\
 					./shared/String/sharedstring.h\
Index: /issm/trunk/src/c/include/typedefs.h
===================================================================
--- /issm/trunk/src/c/include/typedefs.h	(revision 3158)
+++ /issm/trunk/src/c/include/typedefs.h	(revision 3159)
@@ -7,5 +7,8 @@
 
 #define UNDEF -9999
+#define ONETHIRD 0.333333333333333333333333333333333333333333333333333333333333
+#define SQRT2 1.414213562373095048801688724209698078569671875376948073176679738
+#define SQRT3 1.732050807568877293527446341505872366942805253810380628055806979
+#define PI 3.141592653589793238462643383279502884197169399375105820974944592308
 
 #endif //ifndef _ISSMTYPEDEFS_H_
-
Index: /issm/trunk/src/c/objects/Beam.cpp
===================================================================
--- /issm/trunk/src/c/objects/Beam.cpp	(revision 3158)
+++ /issm/trunk/src/c/objects/Beam.cpp	(revision 3159)
@@ -179,5 +179,5 @@
 
 	/*Get node data: */
-	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
         
 	/*Get dof list on which we will plug the pressure values: */
@@ -369,5 +369,5 @@
 
 	//Get all element grid data:
-	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
 	for(i=0;i<numgrids;i++)z_list[i]=xyz_list[i][2];
 	
Index: /issm/trunk/src/c/objects/Icefront.cpp
===================================================================
--- /issm/trunk/src/c/objects/Icefront.cpp	(revision 3158)
+++ /issm/trunk/src/c/objects/Icefront.cpp	(revision 3159)
@@ -315,5 +315,5 @@
 	/* Get dof list and node coordinates: */
 	GetDofList(&doflist[0],&numberofdofspernode);
-	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
 	
 	/*Now build thickness_list_element and bed_list_element: */
@@ -423,5 +423,5 @@
 	/* Get dof list and node coordinates: */
 	GetDofList(&doflist[0],&numberofdofspernode);
-	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
 	
 	//Identify which grids are comprised in the quad: 
@@ -582,5 +582,5 @@
 	/* Get dof list and node coordinates: */
 	GetDofList(&doflist[0],&numberofdofspernode);
-	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
 	
 	//Identify which grids are comprised in the quad: 
Index: /issm/trunk/src/c/objects/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Penta.cpp	(revision 3158)
+++ /issm/trunk/src/c/objects/Penta.cpp	(revision 3159)
@@ -229,5 +229,5 @@
 
 	/*Get node data: */
-	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
 
 	/*pressure is lithostatic: */
@@ -510,5 +510,5 @@
 
 		/* Get node coordinates and dof list: */
-		GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+		GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
 		GetDofList(&doflist[0],&numberofdofspernode);
 
@@ -743,5 +743,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
 	GetDofList(&doflist[0],&numberofdofspernode);
 
@@ -995,5 +995,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
 	GetDofList(&doflist[0],&numberofdofspernode);
 
@@ -1228,5 +1228,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
 	GetDofList(&doflist[0],&numberofdofspernode);
 
@@ -1576,5 +1576,5 @@
 
 		/* Get node coordinates and dof list: */
-		GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+		GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
 		GetDofList(&doflist[0],&numberofdofspernode);
 
@@ -1743,5 +1743,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
 	GetDofList(&doflist[0],&numberofdofspernode);
 
@@ -1962,5 +1962,5 @@
 	/*Now, handle the standard penta element: */
 	/* Get node coordinates and dof list: */
-	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
 	GetDofList(&doflist[0],&numberofdofspernode);
 
@@ -2172,5 +2172,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
 	GetDofList(&doflist[0],&numberofdofspernode);
 
@@ -3049,5 +3049,5 @@
 	double z1,z2,z3,z4,z5,z6;
 
-	double sqrt3=sqrt(3.0);
+	double sqrt3=SQRT3;
 
 	/*Figure out xi,eta and zi (parametric coordinates), for this gaussian point: */
@@ -3461,9 +3461,9 @@
 	const int numgrids=6;
 	double A1,A2,A3,z;
-	double sqrt3=sqrt(3.0);
-
-	A1=gauss_coord[0]; //first area coordinate value. In term of xi and eta: A1=(1-xi)/2-eta/(2*sqrt(3));
-	A2=gauss_coord[1]; //second area coordinate value In term of xi and eta: A2=(1+xi)/2-eta/(2*sqrt(3));
-	A3=gauss_coord[2]; //third area coordinate value  In term of xi and eta: A3=y/sqrt(3);
+	double sqrt3=SQRT3;
+
+	A1=gauss_coord[0]; //first area coordinate value. In term of xi and eta: A1=(1-xi)/2-eta/(2*SQRT3);
+	A2=gauss_coord[1]; //second area coordinate value In term of xi and eta: A2=(1+xi)/2-eta/(2*SQRT3);
+	A3=gauss_coord[2]; //third area coordinate value  In term of xi and eta: A3=y/SQRT3;
 	z=gauss_coord[3]; //fourth vertical coordinate value. Corresponding nodal function: (1-z)/2 and (1+z)/2
 
@@ -3508,5 +3508,5 @@
 	 * natural coordinate system) at the gaussian point. */
 
-	double sqrt3=sqrt(3.0);
+	double sqrt3=SQRT3;
 	int    numgrids=7; //six plus bubble grids
 
Index: /issm/trunk/src/c/objects/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Tria.cpp	(revision 3158)
+++ /issm/trunk/src/c/objects/Tria.cpp	(revision 3159)
@@ -382,5 +382,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
 	GetDofList(&doflist[0],&numberofdofspernode);
 
@@ -586,5 +586,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
 	GetDofList(&doflist[0],&numberofdofspernode);
 
@@ -795,5 +795,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
 	GetDofList(&doflist[0],&numberofdofspernode);
 
@@ -1005,5 +1005,5 @@
 	
 	/* Get node coordinates and dof list: */
-	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
 	GetDofList(&doflist[0],&numberofdofspernode);
 
@@ -1167,5 +1167,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
 	GetDofList(&doflist[0],&numberofdofspernode);
 
@@ -1287,5 +1287,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
 	GetDofList(&doflist[0],&numberofdofspernode);
 
@@ -1335,5 +1335,4 @@
 #define __FUNCT__ "Tria::CreateKMatrixPrognostic"
 void  Tria::CreateKMatrixPrognostic(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
-
 
 	/* local declarations */
@@ -1364,9 +1363,8 @@
 	double DLprime[2][2]={0.0};
 	double DL_scalar;
-	double Ke_gg[numdof][numdof]={0.0};//local element stiffness matrix 
-	double Ke_gg_gaussian[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
-	double Ke_gg_thickness1[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
-	double Ke_gg_thickness2[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
-
+	double Ke_gg[numdof][numdof]={0.0};
+	double Ke_gg_gaussian[numdof][numdof]={0.0};
+	double Ke_gg_thickness1[numdof][numdof]={0.0};
+	double Ke_gg_thickness2[numdof][numdof]={0.0};
 	double Jdettria;
 
@@ -1384,5 +1382,5 @@
 	double  dt;
 	int     dofs[2]={0,1};
-	int     found=0;
+	int     found;
 
 	ParameterInputs* inputs=NULL;
@@ -1404,9 +1402,9 @@
 
 	/* Get node coordinates and dof list: */
-	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
 	GetDofList(&doflist[0],&numberofdofspernode);
 
 	//Create Artificial diffusivity once for all if requested
-	if(numpar->artdiff){
+	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;
@@ -1426,4 +1424,5 @@
 	/* Start  looping on the number of gaussian points: */
 	for (ig=0; ig<num_gauss; ig++){
+
 		/*Pick up the gaussian point: */
 		gauss_weight=*(gauss_weights+ig);
@@ -1487,5 +1486,5 @@
 		for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_thickness2[i][j];
 
-		if(numpar->artdiff){
+		if(numpar->artdiff==1){
 
 			/* Compute artificial diffusivity */
@@ -1588,5 +1587,5 @@
 	
 	/* Get node coordinates and dof list: */
-	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
 	GetDofList(&doflist[0],&numberofdofspernode);
 
@@ -1682,5 +1681,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
 	GetDofList(&doflist[0],&numberofdofspernode);
 
@@ -1830,5 +1829,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
 	GetDofList(&doflist[0],&numberofdofspernode);
 
@@ -1921,5 +1920,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
 	GetDofList(&doflist[0],&numberofdofspernode);
 
@@ -2019,5 +2018,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
 	GetDofList(&doflist[0],&numberofdofspernode);
 
@@ -2129,5 +2128,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
 	GetDofList(&doflist[0],&numberofdofspernode);
 
@@ -2307,5 +2306,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
 	GetDofList(&doflist[0],&numberofdofspernode);
 
@@ -2386,5 +2385,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
 	GetDofList(&doflist[0],&numberofdofspernode);
 
@@ -2494,5 +2493,5 @@
 	
 	/* Get node coordinates and dof list: */
-	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
 	GetDofList(&doflist[0],&numberofdofspernode);
 
@@ -2609,5 +2608,5 @@
 	
 	/* Get node coordinates and dof list: */
-	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
 	GetDofList(&doflist[0],&numberofdofspernode);
 
@@ -2803,5 +2802,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
 	GetDofList(&doflist[0],&numberofdofspernode);
 
@@ -3005,5 +3004,5 @@
 
 	/*Get xyz list: */
-	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
 	x1=xyz_list[0][0]; y1=xyz_list[0][1];
 	x2=xyz_list[1][0]; y2=xyz_list[1][1];
@@ -3027,5 +3026,5 @@
 
 	/*Get xyz list: */
-	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
 	x1=xyz_list[0][0]; y1=xyz_list[0][1];
 	x2=xyz_list[1][0]; y2=xyz_list[1][1];
@@ -3254,5 +3253,4 @@
 	const int numgrids=3;
 	double x1,y1,x2,y2,x3,y3;
-	double sqrt3=sqrt(3.0);
 	
 	x1=*(xyz_list+numgrids*0+0);
@@ -3265,7 +3263,7 @@
 
 	*(J+NDOF2*0+0)=1.0/2.0*(x2-x1);
-	*(J+NDOF2*1+0)=sqrt3/6.0*(2*x3-x1-x2);
+	*(J+NDOF2*1+0)=SQRT3/6.0*(2*x3-x1-x2);
 	*(J+NDOF2*0+1)=1.0/2.0*(y2-y1);
-	*(J+NDOF2*1+1)=sqrt3/6.0*(2*y3-y1-y2);
+	*(J+NDOF2*1+1)=SQRT3/6.0*(2*y3-y1-y2);
 }
 /*}}}*/
@@ -3288,5 +3286,5 @@
 
 
-	*Jdet=sqrt(3.0)/6.0*((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1));
+	*Jdet=SQRT3/6.0*((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1));
 
 
@@ -3318,5 +3316,5 @@
 
 
-	*Jdet=sqrt(3.0)/6.0*pow(pow(((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1)),2.0)+pow(((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1)),2.0)+pow(((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)),2.0),0.5);
+	*Jdet=SQRT3/6.0*pow(pow(((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1)),2.0)+pow(((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1)),2.0)+pow(((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)),2.0),0.5);
 
 
@@ -3469,17 +3467,15 @@
 	const int numgrids=3;
 
-	double sqrt3=sqrt(3.0);
-
 	/*First nodal function: */
 	*(dl1dl3+numgrids*0+0)=-1.0/2.0; 
-	*(dl1dl3+numgrids*1+0)=-1.0/(2.0*sqrt3);
+	*(dl1dl3+numgrids*1+0)=-1.0/(2.0*SQRT3);
 
 	/*Second nodal function: */
 	*(dl1dl3+numgrids*0+1)=1.0/2.0;
-	*(dl1dl3+numgrids*1+1)=-1.0/(2.0*sqrt3);
+	*(dl1dl3+numgrids*1+1)=-1.0/(2.0*SQRT3);
 
 	/*Third nodal function: */
 	*(dl1dl3+numgrids*0+2)=0;
-	*(dl1dl3+numgrids*1+2)=1.0/sqrt3;
+	*(dl1dl3+numgrids*1+2)=1.0/SQRT3;
 
 }
@@ -3666,5 +3662,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
 	GetDofList1(&doflist1[0]);
 
@@ -3843,5 +3839,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
 	GetDofList1(&doflist1[0]);
 
@@ -4065,5 +4061,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
 	GetDofList1(&doflist1[0]);
 
@@ -4263,5 +4259,5 @@
 
 	/*Get xyz list: */
-	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
 
 	/*recover velocity at three element nodes: */
@@ -4383,5 +4379,5 @@
 
 	/* Get node coordinates and dof list: */
-	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
 
 	/* Recover input data: */
Index: sm/trunk/src/c/shared/Elements/GetElementNodeData.cpp
===================================================================
--- /issm/trunk/src/c/shared/Elements/GetElementNodeData.cpp	(revision 3158)
+++ 	(revision )
@@ -1,18 +1,0 @@
-/*!\file:  GetElementNodeData.cpp
- * \brief get node coordinates
- */ 
-
-#include "../../objects/Node.h"
-
-int GetElementNodeData( double* xyz,  Node** nodes, int numgrids){
-
-	int i;
-
-	for ( i=0; i<numgrids; i++) {
-		xyz[i*3+0]=nodes[i]->GetX();
-		xyz[i*3+1]=nodes[i]->GetY();
-		xyz[i*3+2]=nodes[i]->GetZ();
-	}
-
-}
-
Index: /issm/trunk/src/c/shared/Elements/ResolvePointers.cpp
===================================================================
--- /issm/trunk/src/c/shared/Elements/ResolvePointers.cpp	(revision 3158)
+++ /issm/trunk/src/c/shared/Elements/ResolvePointers.cpp	(revision 3159)
@@ -25,5 +25,4 @@
 	Object* object=NULL;
 	int i;
-	
 
 	for(i=0;i<num_objects;i++){
@@ -57,5 +56,3 @@
 		}
 	}
-	
 }
-
Index: /issm/trunk/src/c/shared/Elements/elements.h
===================================================================
--- /issm/trunk/src/c/shared/Elements/elements.h	(revision 3158)
+++ /issm/trunk/src/c/shared/Elements/elements.h	(revision 3159)
@@ -11,5 +11,5 @@
 void ResolvePointers(Object** objects,int* object_ids,int* object_offsets,int num_objects,DataSet* dataset);
 double Paterson(double temperature);
-int GetElementNodeData( double* xyz,  Node** nodes, int numgrids);
+int GetVerticesCoordinates(double* xyz,  Node** nodes, int numgrids);
 
 #endif //ifndef _SHARED_ELEMENTS_H_
Index: /issm/trunk/src/c/shared/Numerics/GaussPoints.cpp
===================================================================
--- /issm/trunk/src/c/shared/Numerics/GaussPoints.cpp	(revision 3158)
+++ /issm/trunk/src/c/shared/Numerics/GaussPoints.cpp	(revision 3159)
@@ -6,8 +6,8 @@
 #include "../Alloc/alloc.h"
 #include "../../include/macros.h"
+#include "../../include/typedefs.h"
 #include "../Exceptions/exceptions.h"
 #include <math.h>
 #include <float.h>
-
 
 /*=======1=========2=========3=========4=========5=========6=========7=========8
@@ -1605,10 +1605,10 @@
 			for (i=0; i<nigaus; i++) {
 				xi        =      1./2.*(1.-egaus[j])*xgaus[i];
-				eta       =sqrt(3.)/2.*(1.+egaus[j]);
-				(*pwgt)[ipt]=xwgt[i]*ewgt[j]*(sqrt(3.)/4.*(1.-egaus[j]));
-
-				(*pl1 )[ipt]=(1.-xi-eta/sqrt(3.))/2.;
-				(*pl2 )[ipt]=(1.+xi-eta/sqrt(3.))/2.;
-				(*pl3 )[ipt]=       eta/sqrt(3.);
+				eta       =SQRT3/2.*(1.+egaus[j]);
+				(*pwgt)[ipt]=xwgt[i]*ewgt[j]*(SQRT3/4.*(1.-egaus[j]));
+
+				(*pl1 )[ipt]=(1.-xi-eta/SQRT3)/2.;
+				(*pl2 )[ipt]=(1.+xi-eta/SQRT3)/2.;
+				(*pl3 )[ipt]=       eta/SQRT3;
 
 				ipt++;
@@ -1955,5 +1955,5 @@
 			(*pl3 )[i]=l3p [iord-1][i];
 			(*pl4 )[i]=l4p [iord-1][i];
-			(*pwgt)[i]=wgtp[iord-1][i]*2./3.*sqrt(2.);
+			(*pwgt)[i]=wgtp[iord-1][i]*2./3.*SQRT2;
 		}
 	}
@@ -1989,15 +1989,15 @@
 				for (i=0; i<nigaus; i++) {
 					xi        =1./4.*(1.-egaus[j])*(1.-zgaus[k])*xgaus[i];
-					eta       =1./4./sqrt(3.)
+					eta       =1./4./SQRT3
 							   *(5.+3.*egaus[j]-zgaus[k]-3.*egaus[j]*zgaus[k]);
 					zeta      =sqrt(2./3.)*(1.+zgaus[k]);
 					(*pwgt)[ipt]=xwgt[i]*ewgt[j]*zwgt[k]
-								 *(sqrt(2.)/16.
+								 *(SQRT2/16.
 							 	   *(1.-egaus[j])*pow(1.-zgaus[k],2.));
 
-					(*pl1 )[ipt]=(1.-xi-eta/sqrt(3.)-zeta/sqrt(6.))/2.;
-					(*pl2 )[ipt]=(1.+xi-eta/sqrt(3.)-zeta/sqrt(6.))/2.;
-					(*pl3 )[ipt]=(      eta         -zeta/sqrt(8.))/sqrt(3.);
-					(*pl4 )[ipt]=(                   zeta/sqrt(8.))*sqrt(3.);
+					(*pl1 )[ipt]=(1.-xi-eta/SQRT3-zeta/sqrt(6.))/2.;
+					(*pl2 )[ipt]=(1.+xi-eta/SQRT3-zeta/sqrt(6.))/2.;
+					(*pl3 )[ipt]=(      eta         -zeta/sqrt(8.))/SQRT3;
+					(*pl4 )[ipt]=(                   zeta/sqrt(8.))*SQRT3;
 
 					ipt++;