Index: /issm/trunk/src/c/objects/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Tria.cpp	(revision 231)
+++ /issm/trunk/src/c/objects/Tria.cpp	(revision 232)
@@ -1369,10 +1369,10 @@
 	double obs_vx_list[numgrids];
 	double obs_vy_list[numgrids];
-	double absolutex[numgrids];
-	double absolutey[numgrids];
-	double relativex[numgrids];
-	double relativey[numgrids];
-	double logarithmicx[numgrids];
-	double logarithmicy[numgrids];
+	double absolutex_list[numgrids];
+	double absolutey_list[numgrids];
+	double relativex_list[numgrids];
+	double relativey_list[numgrids];
+	double logarithmicx_list[numgrids];
+	double logarithmicy_list[numgrids];
 
 	/* gaussian points: */
@@ -1388,4 +1388,5 @@
 	double  velocity_x,velocity_y,velocity_mag;
 	double  obs_velocity_x,obs_velocity_y,obs_velocity_mag;
+	double  absolutex,absolutey,relativex,relativey,logarithmicx,logarithmicy;
 
 	/*element vector : */
@@ -1432,6 +1433,6 @@
 		/*We are using an absolute misfit: */
 		for (i=0;i<numgrids;i++){
-			absolutex[i]=vx_list[i]-obs_vx_list[i];
-			absolutey[i]=vy_list[i]-obs_vy_list[i];
+			absolutex_list[i]=vx_list[i]-obs_vx_list[i];
+			absolutey_list[i]=vy_list[i]-obs_vy_list[i];
 		}
 	}
@@ -1443,6 +1444,6 @@
 			if(obs_vx_list[i]==0)scalex=0;
 			if(obs_vy_list[i]==0)scaley=0;
-			relativex[i]=scalex*(vx_list[i]-obs_vx_list[i]);
-			relativey[i]=scaley*(vy_list[i]-obs_vy_list[i]);
+			relativex_list[i]=scalex*(vx_list[i]-obs_vx_list[i]);
+			relativey_list[i]=scaley*(vy_list[i]-obs_vy_list[i]);
 		}
 	}
@@ -1453,6 +1454,6 @@
 			obs_velocity_mag=sqrt(pow(obs_vx_list[i],2)+pow(obs_vy_list[i],2))+epsvel; //epsvel to avoid observed velocity being nil.
 			scale=8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag);
-			logarithmicx[i]=scale*vx_list[i];
-			logarithmicy[i]=scale*vy_list[i];
+			logarithmicx_list[i]=scale*vx_list[i];
+			logarithmicy_list[i]=scale*vy_list[i];
 		}
 	}
@@ -1509,28 +1510,39 @@
 		if(fit==0){
 			/*We are using an absolute misfit: */
+
+			/*Compute absolute(x/y) at gaussian point: */
+			GetParameterValue(&absolutex, &absolutex_list[0],gauss_l1l2l3);
+			GetParameterValue(&absolutey, &absolutey_list[0],gauss_l1l2l3);
+
+			/*compute Du*/
 			for (i=0;i<numgrids;i++){
-				due_g_gaussian[i*NDOF2+0]=absolutex[i]*Jdet*gauss_weight*l1l2l3[i]; 
-				due_g_gaussian[i*NDOF2+1]=absolutey[i]*Jdet*gauss_weight*l1l2l3[i]; 
+				due_g_gaussian[i*NDOF2+0]=absolutex*Jdet*gauss_weight*l1l2l3[i]; 
+				due_g_gaussian[i*NDOF2+1]=absolutey*Jdet*gauss_weight*l1l2l3[i]; 
 			}
 		}
 		else if(fit==1){
 			/*We are using a relative misfit: */
-			scalex=pow(meanvel/(obs_velocity_x+epsvel),2);
-			scaley=pow(meanvel/(obs_velocity_y+epsvel),2);
-			if(obs_velocity_x==0)scalex=0;
-			if(obs_velocity_y==0)scaley=0;
+
+			/*Compute relative(x/y) at gaussian point: */
+			GetParameterValue(&relativex, &relativex_list[0],gauss_l1l2l3);
+			GetParameterValue(&relativey, &relativey_list[0],gauss_l1l2l3);
+
+			/*compute Du*/
 			for (i=0;i<numgrids;i++){
-				due_g_gaussian[i*NDOF2+0]=relativex[i]*Jdet*gauss_weight*l1l2l3[i]; 
-				due_g_gaussian[i*NDOF2+1]=relativey[i]*Jdet*gauss_weight*l1l2l3[i]; 
+				due_g_gaussian[i*NDOF2+0]=relativex*Jdet*gauss_weight*l1l2l3[i]; 
+				due_g_gaussian[i*NDOF2+1]=relativey*Jdet*gauss_weight*l1l2l3[i]; 
 			}
 		}
 		else if(fit==2){
 			/*We are using a logarithmic misfit: */
-			velocity_mag=sqrt(pow(velocity_x,2)+pow(velocity_y,2))+epsvel; //epsvel to avoid velocity being nil.
-			obs_velocity_mag=sqrt(pow(obs_velocity_x,2)+pow(obs_velocity_y,2))+epsvel; //epsvel to avoid observed velocity being nil.
-			scale=8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag);
+
+			/*Compute logarithmic(x/y) at gaussian point: */
+			GetParameterValue(&logarithmicx, &logarithmicx_list[0],gauss_l1l2l3);
+			GetParameterValue(&logarithmicy, &logarithmicy_list[0],gauss_l1l2l3);
+
+			/*compute Du*/
 			for (i=0;i<numgrids;i++){
-				due_g_gaussian[i*NDOF2+0]=logarithmicx[i]*Jdet*gauss_weight*l1l2l3[i]; 
-				due_g_gaussian[i*NDOF2+1]=logarithmicy[i]*Jdet*gauss_weight*l1l2l3[i]; 
+				due_g_gaussian[i*NDOF2+0]=logarithmicx*Jdet*gauss_weight*l1l2l3[i]; 
+				due_g_gaussian[i*NDOF2+1]=logarithmicy*Jdet*gauss_weight*l1l2l3[i]; 
 			}
 		}
@@ -1816,5 +1828,5 @@
 
 	/*element vector at the gaussian points: */
-	double  grade_g[numgrids];
+	double  grade_g[numdofs];
 	double  grade_g_gaussian[numgrids];
 
@@ -1912,5 +1924,5 @@
 
 		/*Add grade_g_gaussian to grade_g: */
-		for( i=0; i<numdofs;i++)grade_g[i]+=grade_g_gaussian[i];
+		for( i=0; i<numgrids;i++) grade_g[i*numberofdofspernode]+=grade_g_gaussian[i];
 	}