Index: /issm/trunk/src/c/objects/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Tria.cpp	(revision 204)
+++ /issm/trunk/src/c/objects/Tria.cpp	(revision 205)
@@ -663,11 +663,5 @@
 		GetParameterValue(&alpha2, &alpha2_list[0],gauss_l1l2l3);
 
-		if (velocity_param){
-			DL_scalar=alpha2*gauss_weight*Jdet;
-		}
-		else{
-			DL_scalar=0;
-		}
-			
+		DL_scalar=alpha2*gauss_weight*Jdet;
 		for (i=0;i<2;i++){
 			DL[i][i]=DL_scalar;
@@ -922,4 +916,5 @@
 	B_param=ParameterInputsRecover(inputs,"B");
 	if(temperature_average_param){
+		temperature_average-0;
 		for(i=0;i<3;i++) temperature_average+=temperature_average_param[doflist[i*numberofdofspernode+0]];
 		temperature_average=temperature_average/3.0;
@@ -1374,4 +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];
 
 	/* gaussian points: */
@@ -1399,8 +1400,6 @@
 
 	/*relative and algorithmic fitting: */
-	double meanvel=0;
-	double epsvel=0;
-	double scalex=1;
-	double scaley=1;
+	double scalex=0;
+	double scaley=0;
 	double scale=0;
 	double* fit_param=NULL;
@@ -1427,7 +1426,39 @@
 		obs_vx_list[i]=obs_velocity[doflist[i*numberofdofspernode+0]];
 		obs_vy_list[i]=obs_velocity[doflist[i*numberofdofspernode+1]];
-
-	}
-
+	}
+
+	/*Get Du at the 3 nodes (integration of the linearized function)*/
+	if(fit==0){
+		/*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];
+		}
+	}
+	else if(fit==1){
+		/*We are using a relative misfit: */
+		for (i=0;i<numgrids;i++){
+			scalex=pow(meanvel/(obs_vx_list[i]+epsvel),2);
+			scaley=pow(meanvel/(obs_vy_list[i]+epsvel),2);
+			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]);
+		}
+	}
+	else if(fit==2){
+		/*We are using a logarithmic misfit: */
+		for (i=0;i<numgrids;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.
+			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];
+		}
+	}
+	else{
+		/*Not supported yet! : */
+		throw ErrorException(__FUNCT__,exprintf("%s%g","unsupported type of fit: ",fit));
+	}
 	
 	/* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
@@ -1479,6 +1510,6 @@
 			/*We are using an absolute misfit: */
 			for (i=0;i<numgrids;i++){
-				due_g_gaussian[i*NDOF2+0]=(velocity_x-obs_velocity_x)*Jdet*gauss_weight*l1l2l3[i]; 
-				due_g_gaussian[i*NDOF2+1]=(velocity_y-obs_velocity_y)*Jdet*gauss_weight*l1l2l3[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]; 
 			}
 		}
@@ -1490,10 +1521,9 @@
 			if(obs_velocity_y==0)scaley=0;
 			for (i=0;i<numgrids;i++){
-				due_g_gaussian[i*NDOF2+0]=scalex*(velocity_x-obs_velocity_x)*Jdet*gauss_weight*l1l2l3[i]; 
-				due_g_gaussian[i*NDOF2+1]=scaley*(velocity_y-obs_velocity_y)*Jdet*gauss_weight*l1l2l3[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]; 
 			}
 		}
 		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.
@@ -1501,6 +1531,6 @@
 			scale=8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag);
 			for (i=0;i<numgrids;i++){
-				due_g_gaussian[i*NDOF2+0]=scale*velocity_x*Jdet*gauss_weight*l1l2l3[i]; 
-				due_g_gaussian[i*NDOF2+1]=scale*velocity_y*Jdet*gauss_weight*l1l2l3[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]; 
 			}
 		}
@@ -1510,5 +1540,4 @@
 		}
 		
-		
 		/*Add due_g_gaussian vector to due_g: */
 		for( i=0; i<numdofs; i++){
@@ -1516,5 +1545,4 @@
 		}
 	}
-	
 	
 	/*Add due_g to global vector du_g: */
@@ -2127,4 +2155,5 @@
 					&Ke_gg_gaussian[0][0],0);
 
+		
 		/* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
 		for( i=0; i<numdofs; i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
@@ -2356,6 +2385,6 @@
 		/*Get bed slope: */
 		GetParameterDerivativeValue(&slope[0], &b[0],&xyz_list[0][0], gauss_l1l2l3);
-		dbdx=slope[0];
-		dbdy=slope[1];
+		dbdx=-slope[0];
+		dbdy=-slope[1];
 
 		/* Get Jacobian determinant: */