Index: /issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 3830)
+++ /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 3831)
@@ -972,4 +972,5 @@
 	int          doflist[numdof];
 	int          numberofdofspernode;
+	double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
 
 	/* gaussian points: */
@@ -981,5 +982,4 @@
 	double  gauss_weight;
 	double  gauss_l1l2l3[3];
-	double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
 
 	/* matrices: */
@@ -1000,6 +1000,4 @@
 	double  surface_list[3];
 	double  nx,ny,norm;
-	double  vx_list[numgrids]={0.0};
-	double  vy_list[numgrids]={0.0};
 	double  dvx[2];
 	double  dvy[2];
@@ -1015,8 +1013,4 @@
 	double artdiff;
 
-	/*Recover velocity: */
-	inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3,VxAverageEnum);
-	inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyAverageEnum);
-
 	/*retrieve some parameters: */
 	this->parameters->FindParam(&artdiff,ArtDiffEnum);
@@ -1031,6 +1025,6 @@
 
 	/*Get normal vector to the surface*/
-	nx=(vx_list[0]+vx_list[1]+vx_list[2])/3;
-	ny=(vy_list[0]+vy_list[1]+vy_list[2])/3;
+	inputs->GetParameterAverage(&nx,VxAverageEnum);
+	inputs->GetParameterAverage(&ny,VyAverageEnum);
 	if(nx==0 && ny==0){
 		SurfaceNormal(&surface_normal[0],xyz_list);
@@ -1053,6 +1047,6 @@
 
 		//Build K matrix (artificial diffusivity matrix)
-		v_gauss[0]=ONETHIRD*(vx_list[0]+vx_list[1]+vx_list[2]);
-		v_gauss[1]=ONETHIRD*(vy_list[0]+vy_list[1]+vy_list[2]);
+		inputs->GetParameterAverage(&v_gauss[0],VxAverageEnum);
+		inputs->GetParameterAverage(&v_gauss[1],VyAverageEnum);
 
 		K[0][0]=pow(10,2)*pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]); //pow should be zero!!
@@ -1079,9 +1073,9 @@
 
 		//Get vx, vy and their derivatives at gauss point
-		GetParameterValue(&vx, &vx_list[0],gauss_l1l2l3);
-		GetParameterValue(&vy, &vy_list[0],gauss_l1l2l3);
-
-		GetParameterDerivativeValue(&dvx[0], &vx_list[0],&xyz_list[0][0], gauss_l1l2l3);
-		GetParameterDerivativeValue(&dvy[0], &vy_list[0],&xyz_list[0][0], gauss_l1l2l3);
+		inputs->GetParameterValue(&vx,&gauss_l1l2l3[0],VxAverageEnum);
+		inputs->GetParameterValue(&vy,&gauss_l1l2l3[0],VyAverageEnum);
+
+		inputs->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],&gauss_l1l2l3[0],VxAverageEnum);
+		inputs->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],&gauss_l1l2l3[0],VyAverageEnum);
 
 		dvxdx=dvx[0];
@@ -1654,5 +1648,4 @@
 	double  gauss_weight;
 	double  gauss_l1l2l3[3];
-	double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
 
 	/* matrices: */
@@ -1670,6 +1663,4 @@
 
 	/*input parameters for structural analysis (diagnostic): */
-	double  vx_list[numgrids]={0.0};
-	double  vy_list[numgrids]={0.0};
 	double  dvx[2];
 	double  dvy[2];
@@ -1694,8 +1685,4 @@
 	GetDofList(&doflist[0],&numberofdofspernode);
 
-	/*Recover velocity: */
-	inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3,VxAverageEnum);
-	inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyAverageEnum);
-
 	//Create Artificial diffusivity once for all if requested
 	if(artdiff==1){
@@ -1705,6 +1692,6 @@
 
 		//Build K matrix (artificial diffusivity matrix)
-		v_gauss[0]=ONETHIRD*(vx_list[0]+vx_list[1]+vx_list[2]);
-		v_gauss[1]=ONETHIRD*(vy_list[0]+vy_list[1]+vy_list[2]);
+		inputs->GetParameterAverage(&v_gauss[0],VxAverageEnum);
+		inputs->GetParameterAverage(&v_gauss[1],VyAverageEnum);
 
 		K[0][0]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]);
@@ -1743,9 +1730,9 @@
 
 		//Get vx, vy and their derivatives at gauss point
-		GetParameterValue(&vx, &vx_list[0],gauss_l1l2l3);
-		GetParameterValue(&vy, &vy_list[0],gauss_l1l2l3);
-
-		GetParameterDerivativeValue(&dvx[0], &vx_list[0],&xyz_list[0][0], gauss_l1l2l3);
-		GetParameterDerivativeValue(&dvy[0], &vy_list[0],&xyz_list[0][0], gauss_l1l2l3);
+		inputs->GetParameterValue(&vx,&gauss_l1l2l3[0],VxAverageEnum);
+		inputs->GetParameterValue(&vy,&gauss_l1l2l3[0],VyAverageEnum);
+
+		inputs->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],&gauss_l1l2l3[0],VxAverageEnum);
+		inputs->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],&gauss_l1l2l3[0],VyAverageEnum);
 
 		dvxdx=dvx[0];