Index: /issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 8595)
+++ /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 8596)
@@ -1655,12 +1655,6 @@
 	double     dux,duy,meanvel,epsvel;
 	double     scalex=0,scaley=0,scale=0,S=0;
+	double     vx,vy,vxobs,vyobs,weight;
 	double     xyz_list[NUMVERTICES][3];
-	double     vx_list[NUMVERTICES];
-	double     vy_list[NUMVERTICES];
-	double     obs_vx_list[NUMVERTICES];
-	double     obs_vy_list[NUMVERTICES];
-	double     dux_list[NUMVERTICES];
-	double     duy_list[NUMVERTICES];
-	double     weights_list[NUMVERTICES];
 	double     l1l2l3[3];
 	GaussTria* gauss=NULL;
@@ -1673,147 +1667,128 @@
 	this->parameters->FindParam(&meanvel,MeanVelEnum);
 	this->parameters->FindParam(&epsvel,EpsVelEnum);
-	GetParameterListOnVertices(&obs_vx_list[0],VxObsEnum);
-	GetParameterListOnVertices(&obs_vy_list[0],VyObsEnum);
-	GetParameterListOnVertices(&vx_list[0],VxEnum);
-	GetParameterListOnVertices(&vy_list[0],VyEnum);
-	GetParameterListOnVertices(&weights_list[0],WeightsEnum,0,0);
+	Input* weights_input=inputs->GetInput(WeightsEnum);   _assert_(weights_input);
+	Input* vx_input     =inputs->GetInput(VxEnum);        _assert_(vx_input);
+	Input* vy_input     =inputs->GetInput(VyEnum);        _assert_(vy_input);
+	Input* vxobs_input  =inputs->GetInput(VxObsEnum);     _assert_(vxobs_input);
+	Input* vyobs_input  =inputs->GetInput(VyObsEnum);     _assert_(vyobs_input);
+
 	inputs->GetParameterValue(&response,CmResponseEnum);
-	if(response==SurfaceAverageVelMisfitEnum){
-		inputs->GetParameterValue(&S,SurfaceAreaEnum);
-	}
-
-	/*Get Du at the 3 nodes (integration of the linearized function)
-	 * Here we integrate linearized functions:
-	 *               
-	 * J(E) = int_E   sum_{i=1}^3  J_i Phi_i
-	 *
-	 *       d J                  dJ_i
-	 * DU= - --- = sum_{i=1}^3  - ---  Phi_i = sum_{i=1}^3 DU_i Phi_i
-	 *       d u                  du_i
-	 *
-	 * where J_i are the misfits at the 3 nodes of the triangle
-	 *       Phi_i is the nodal function (P1) with respect to 
-	 *       the vertex i
-	 */
-	if(response==SurfaceAbsVelMisfitEnum){
-		/*We are using an absolute misfit:
-		 *
-		 *      1  [           2              2 ]
-		 * J = --- | (u - u   )  +  (v - v   )  |
-		 *      2  [       obs            obs   ]
-		 *
-		 *        dJ
-		 * DU = - -- = (u   - u )
-		 *        du     obs
-		 */
-		for (i=0;i<NUMVERTICES;i++){
-			dux_list[i]=obs_vx_list[i]-vx_list[i];
-			duy_list[i]=obs_vy_list[i]-vy_list[i];
-		}
-	}
-	else if(response==SurfaceRelVelMisfitEnum){
-		/*We are using a relative misfit: 
-		 *                        
-		 *      1  [     \bar{v}^2             2   \bar{v}^2              2 ]
-		 * J = --- | -------------  (u - u   ) + -------------  (v - v   )  |
-		 *      2  [  (u   + eps)^2       obs    (v   + eps)^2       obs    ]
-		 *              obs                        obs                      
-		 *
-		 *        dJ     \bar{v}^2
-		 * DU = - -- = ------------- (u   - u )
-		 *        du   (u   + eps)^2    obs
-		 *               obs
-		 */
-		for (i=0;i<NUMVERTICES;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;
-			dux_list[i]=scalex*(obs_vx_list[i]-vx_list[i]);
-			duy_list[i]=scaley*(obs_vy_list[i]-vy_list[i]);
-		}
-	}
-	else if(response==SurfaceLogVelMisfitEnum){
-		/*We are using a logarithmic misfit:
-		 *                        
-		 *                 [        vel + eps     ] 2
-		 * J = 4 \bar{v}^2 | log ( -----------  ) |  
-		 *                 [       vel   + eps    ]
-		 *                            obs
-		 *
-		 *        dJ                 2 * log(...)
-		 * DU = - -- = - 4 \bar{v}^2 -------------  u
-		 *        du                 vel^2 + eps
-		 *            
-		 */
-		for (i=0;i<NUMVERTICES;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);
-			dux_list[i]=scale*vx_list[i];
-			duy_list[i]=scale*vy_list[i];
-		}
-	}
-	else if(response==SurfaceAverageVelMisfitEnum){
-		/*We are using an spacially average absolute misfit:
-		 *
-		 *      1                    2              2
-		 * J = ---  sqrt(  (u - u   )  +  (v - v   )  )
-		 *      S                obs            obs
-		 *
-		 *        dJ      1       1 
-		 * DU = - -- = - --- ----------- * 2 (u - u   )
-		 *        du      S  2 sqrt(...)           obs
-		 */
-		for (i=0;i<NUMVERTICES;i++){
-			scale=1.0/(S*sqrt(pow(vx_list[i]-obs_vx_list[i],2)+pow(vy_list[i]-obs_vx_list[i],2))+epsvel);
-			dux_list[i]=scale*(obs_vx_list[i]-vx_list[i]);
-			duy_list[i]=scale*(obs_vy_list[i]-vy_list[i]);
-		}
-	}
-	else if(response==SurfaceLogVxVyMisfitEnum){
-		/*We are using an logarithmic 2 misfit:
-		 *
-		 *      1            [        |u| + eps     2          |v| + eps     2  ]
-		 * J = --- \bar{v}^2 | log ( -----------  )   +  log ( -----------  )   |  
-		 *      2            [       |u    |+ eps              |v    |+ eps     ]
-		 *                              obs                       obs
-		 *        dJ                              1      u                             1
-		 * DU = - -- = - \bar{v}^2 log(u...) --------- ----  ~ - \bar{v}^2 log(u...) ------
-		 *        du                         |u| + eps  |u|                           u + eps
-		 */
-		for (i=0;i<NUMVERTICES;i++){
-			dux_list[i] = - pow(meanvel,(double)2)*(
-						log((fabs(vx_list[i])+epsvel)/(fabs(obs_vx_list[i])+epsvel)) * 1/(vx_list[i]+epsvel));
-			duy_list[i] = - pow(meanvel,(double)2)*(
-						log((fabs(vy_list[i])+epsvel)/(fabs(obs_vy_list[i])+epsvel)) * 1/(vy_list[i]+epsvel));
-		}
-	}
-	else{
-		/*Not supported yet! : */
-		_error_("response %s not supported yet",EnumToStringx(response));
-	}
-
-	/*Apply weights to DU*/
-	for (i=0;i<NUMVERTICES;i++){
-		dux_list[i]=weights_list[i]*dux_list[i];
-		duy_list[i]=weights_list[i]*duy_list[i];
-	}
+	if(response==SurfaceAverageVelMisfitEnum) inputs->GetParameterValue(&S,SurfaceAreaEnum);
 
 	/* Start  looping on the number of gaussian points: */
-	gauss=new GaussTria(2);
-	for(ig=gauss->begin();ig<gauss->end();ig++){
+	gauss=new GaussTria(4);
+	for (ig=gauss->begin();ig<gauss->end();ig++){
 
 		gauss->GaussPoint(ig);
 
+		/* Get Jacobian determinant: */
 		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
+
+		/*Get all parameters at gaussian point*/
+		weights_input->GetParameterValue(&weight,gauss,0);
+		vx_input->GetParameterValue(&vx,gauss);
+		vy_input->GetParameterValue(&vy,gauss);
+		vxobs_input->GetParameterValue(&vxobs,gauss);
+		vyobs_input->GetParameterValue(&vyobs,gauss);
 		GetNodalFunctions(l1l2l3, gauss);
 
-		TriaRef::GetParameterValue(&dux, &dux_list[0],gauss);
-		TriaRef::GetParameterValue(&duy, &duy_list[0],gauss);
-
-		for (i=0;i<NUMVERTICES;i++){
-			pe->values[i*NDOF2+0]+=dux*Jdet*gauss->weight*l1l2l3[i]; 
-			pe->values[i*NDOF2+1]+=duy*Jdet*gauss->weight*l1l2l3[i]; 
+		if(response==SurfaceAbsVelMisfitEnum){
+			/*
+			 *      1  [           2              2 ]
+			 * J = --- | (u - u   )  +  (v - v   )  |
+			 *      2  [       obs            obs   ]
+			 *
+			 *        dJ
+			 * DU = - -- = (u   - u )
+			 *        du     obs
+			 */
+			for (i=0;i<NUMVERTICES;i++){
+				dux=vxobs-vx;
+				duy=vyobs-vy;
+				pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 
+				pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 
+			}
+		}
+		else if(response==SurfaceRelVelMisfitEnum){
+			/*
+			 *      1  [     \bar{v}^2             2   \bar{v}^2              2 ]
+			 * J = --- | -------------  (u - u   ) + -------------  (v - v   )  |
+			 *      2  [  (u   + eps)^2       obs    (v   + eps)^2       obs    ]
+			 *              obs                        obs                      
+			 *
+			 *        dJ     \bar{v}^2
+			 * DU = - -- = ------------- (u   - u )
+			 *        du   (u   + eps)^2    obs
+			 *               obs
+			 */
+			for (i=0;i<NUMVERTICES;i++){
+				scalex=pow(meanvel/(vxobs+epsvel),2.); if(vxobs==0)scalex=0;
+				scaley=pow(meanvel/(vyobs+epsvel),2.); if(vyobs==0)scaley=0;
+				dux=scalex*(vxobs-vx);
+				duy=scaley*(vyobs-vy);
+				pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 
+				pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 
+			}
+		}
+		else if(response==SurfaceLogVelMisfitEnum){
+			/*
+			 *                 [        vel + eps     ] 2
+			 * J = 4 \bar{v}^2 | log ( -----------  ) |  
+			 *                 [       vel   + eps    ]
+			 *                            obs
+			 *
+			 *        dJ                 2 * log(...)
+			 * DU = - -- = - 4 \bar{v}^2 -------------  u
+			 *        du                 vel^2 + eps
+			 *            
+			 */
+			for (i=0;i<NUMVERTICES;i++){
+				velocity_mag    =sqrt(pow(vx,   2.)+pow(vy,   2.))+epsvel;
+				obs_velocity_mag=sqrt(pow(vxobs,2.)+pow(vyobs,2.))+epsvel;
+				scale=-8*pow(meanvel,2.)/pow(velocity_mag,2.)*log(velocity_mag/obs_velocity_mag);
+				dux=scale*vx;
+				duy=scale*vy;
+				pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 
+				pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 
+			}
+		}
+		else if(response==SurfaceAverageVelMisfitEnum){
+			/*
+			 *      1                    2              2
+			 * J = ---  sqrt(  (u - u   )  +  (v - v   )  )
+			 *      S                obs            obs
+			 *
+			 *        dJ      1       1 
+			 * DU = - -- = - --- ----------- * 2 (u - u   )
+			 *        du      S  2 sqrt(...)           obs
+			 */
+			for (i=0;i<NUMVERTICES;i++){
+				scale=1./(S*2*sqrt(pow(vx-vxobs,2.)+pow(vy-vyobs,2.))+epsvel);
+				dux=scale*(vxobs-vx);
+				duy=scale*(vyobs-vy);
+				pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 
+				pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 
+			}
+		}
+		else if(response==SurfaceLogVxVyMisfitEnum){
+			/*
+			 *      1            [        |u| + eps     2          |v| + eps     2  ]
+			 * J = --- \bar{v}^2 | log ( -----------  )   +  log ( -----------  )   |  
+			 *      2            [       |u    |+ eps              |v    |+ eps     ]
+			 *                              obs                       obs
+			 *        dJ                              1      u                             1
+			 * DU = - -- = - \bar{v}^2 log(u...) --------- ----  ~ - \bar{v}^2 log(u...) ------
+			 *        du                         |u| + eps  |u|                           u + eps
+			 */
+			for (i=0;i<NUMVERTICES;i++){
+				dux = - pow(meanvel,2.) * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
+				duy = - pow(meanvel,2.) * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel);
+				pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 
+				pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 
+			}
+		}
+		else{
+			/*Not supported yet! : */
+			_error_("response %s not supported yet",EnumToStringx(response));
 		}
 	}
@@ -1828,19 +1803,11 @@
 
 	/*Intermediaries */
-	int        i,ig;
-	int        fit=-1;
-	int        response;
+	int        i,ig,response;
 	double     Jdet;
 	double     obs_velocity_mag,velocity_mag;
 	double     dux,duy,meanvel,epsvel;
 	double     scalex=0,scaley=0,scale=0,S=0;
+	double     vx,vy,vxobs,vyobs,weight;
 	double     xyz_list[NUMVERTICES][3];
-	double     vx_list[NUMVERTICES];
-	double     vy_list[NUMVERTICES];
-	double     obs_vx_list[NUMVERTICES];
-	double     obs_vy_list[NUMVERTICES];
-	double     dux_list[NUMVERTICES];
-	double     duy_list[NUMVERTICES];
-	double     weights_list[NUMVERTICES];
 	double     l1l2l3[3];
 	GaussTria* gauss=NULL;
@@ -1853,147 +1820,128 @@
 	this->parameters->FindParam(&meanvel,MeanVelEnum);
 	this->parameters->FindParam(&epsvel,EpsVelEnum);
-	GetParameterListOnVertices(&obs_vx_list[0],VxObsEnum);
-	GetParameterListOnVertices(&obs_vy_list[0],VyObsEnum);
-	GetParameterListOnVertices(&vx_list[0],VxEnum);
-	GetParameterListOnVertices(&vy_list[0],VyEnum);
-	GetParameterListOnVertices(&weights_list[0],WeightsEnum,0,0);
+	Input* weights_input=inputs->GetInput(WeightsEnum);   _assert_(weights_input);
+	Input* vx_input     =inputs->GetInput(VxEnum);        _assert_(vx_input);
+	Input* vy_input     =inputs->GetInput(VyEnum);        _assert_(vy_input);
+	Input* vxobs_input  =inputs->GetInput(VxObsEnum);     _assert_(vxobs_input);
+	Input* vyobs_input  =inputs->GetInput(VyObsEnum);     _assert_(vyobs_input);
+
 	inputs->GetParameterValue(&response,CmResponseEnum);
-	if(response==SurfaceAverageVelMisfitEnum){
-		inputs->GetParameterValue(&S,SurfaceAreaEnum);
-	}
-
-	/*Get Du at the 3 nodes (integration of the linearized function)
-	 * Here we integrate linearized functions:
-	 *               
-	 * J(E) = int_E   sum_{i=1}^3  J_i Phi_i
-	 *
-	 *       d J                  dJ_i
-	 * DU= - --- = sum_{i=1}^3  - ---  Phi_i = sum_{i=1}^3 DU_i Phi_i
-	 *       d u                  du_i
-	 *
-	 * where J_i are the misfits at the 3 nodes of the triangle
-	 *       Phi_i is the nodal function (P1) with respect to 
-	 *       the vertex i
-	 */
-	if(response==SurfaceAbsVelMisfitEnum){
-		/*We are using an absolute misfit:
-		 *
-		 *      1  [           2              2 ]
-		 * J = --- | (u - u   )  +  (v - v   )  |
-		 *      2  [       obs            obs   ]
-		 *
-		 *        dJ             2
-		 * DU = - -- = (u   - u )
-		 *        du     obs
-		 */
-		for (i=0;i<NUMVERTICES;i++){
-			dux_list[i]=obs_vx_list[i]-vx_list[i];
-			duy_list[i]=obs_vy_list[i]-vy_list[i];
-		}
-	}
-	else if(response==SurfaceRelVelMisfitEnum){
-		/*We are using a relative misfit: 
-		 *                        
-		 *      1  [     \bar{v}^2             2   \bar{v}^2              2 ]
-		 * J = --- | -------------  (u - u   ) + -------------  (v - v   )  |
-		 *      2  [  (u   + eps)^2       obs    (v   + eps)^2       obs    ]
-		 *              obs                        obs                      
-		 *
-		 *        dJ     \bar{v}^2
-		 * DU = - -- = ------------- (u   - u )
-		 *        du   (u   + eps)^2    obs
-		 *               obs
-		 */
-		for (i=0;i<NUMVERTICES;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;
-			dux_list[i]=scalex*(obs_vx_list[i]-vx_list[i]);
-			duy_list[i]=scaley*(obs_vy_list[i]-vy_list[i]);
-		}
-	}
-	else if(response==SurfaceLogVelMisfitEnum){
-		/*We are using a logarithmic misfit:
-		 *                        
-		 *                 [        vel + eps     ] 2
-		 * J = 4 \bar{v}^2 | log ( -----------  ) |  
-		 *                 [       vel   + eps    ]
-		 *                            obs
-		 *
-		 *        dJ                 2 * log(...)
-		 * DU = - -- = - 4 \bar{v}^2 -------------  u
-		 *        du                 vel^2 + eps
-		 *            
-		 */
-		for (i=0;i<NUMVERTICES;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);
-			dux_list[i]=scale*vx_list[i];
-			duy_list[i]=scale*vy_list[i];
-		}
-	}
-	else if(response==SurfaceAverageVelMisfitEnum){
-		/*We are using an spacially average absolute misfit:
-		 *
-		 *      1                    2              2
-		 * J = ---  sqrt(  (u - u   )  +  (v - v   )  )
-		 *      S                obs            obs
-		 *
-		 *        dJ      1       1 
-		 * DU = - -- = - --- ----------- * 2 (u - u   )
-		 *        du      S  2 sqrt(...)           obs
-		 */
-		for (i=0;i<NUMVERTICES;i++){
-			scale=1.0/(S*sqrt(pow(vx_list[i]-obs_vx_list[i],2)+pow(vy_list[i]-obs_vx_list[i],2))+epsvel);
-			dux_list[i]=scale*(obs_vx_list[i]-vx_list[i]);
-			duy_list[i]=scale*(obs_vy_list[i]-vy_list[i]);
-		}
-	}
-	else if(response==SurfaceLogVxVyMisfitEnum){
-		/*We are using an logarithmic 2 misfit:
-		 *
-		 *      1            [        |u| + eps     2          |v| + eps     2  ]
-		 * J = --- \bar{v}^2 | log ( -----------  )   +  log ( -----------  )   |  
-		 *      2            [       |u    |+ eps              |v    |+ eps     ]
-		 *                              obs                       obs
-		 *        dJ                              1      u                             1
-		 * DU = - -- = - \bar{v}^2 log(u...) --------- ----  ~ - \bar{v}^2 log(u...) ------
-		 *        du                         |u| + eps  |u|                           u + eps
-		 */
-		for (i=0;i<NUMVERTICES;i++){
-			dux_list[i] = - pow(meanvel,(double)2)*(
-						log((fabs(vx_list[i])+epsvel)/(fabs(obs_vx_list[i])+epsvel)) * 1/(vx_list[i]+epsvel));
-			duy_list[i] = - pow(meanvel,(double)2)*(
-						log((fabs(vy_list[i])+epsvel)/(fabs(obs_vy_list[i])+epsvel)) * 1/(vy_list[i]+epsvel));
-		}
-	}
-	else{
-		/*Not supported yet! : */
-		_error_("response %s not supported yet",EnumToStringx(response));
-	}
-
-	/*Apply weights to DU*/
-	for (i=0;i<NUMVERTICES;i++){
-		dux_list[i]=weights_list[i]*dux_list[i];
-		duy_list[i]=weights_list[i]*duy_list[i];
-	}
+	if(response==SurfaceAverageVelMisfitEnum) inputs->GetParameterValue(&S,SurfaceAreaEnum);
 
 	/* Start  looping on the number of gaussian points: */
-	gauss=new GaussTria(2);
-	for(ig=gauss->begin();ig<gauss->end();ig++){
+	gauss=new GaussTria(4);
+	for (ig=gauss->begin();ig<gauss->end();ig++){
 
 		gauss->GaussPoint(ig);
 
+		/* Get Jacobian determinant: */
 		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
+
+		/*Get all parameters at gaussian point*/
+		weights_input->GetParameterValue(&weight,gauss,0);
+		vx_input->GetParameterValue(&vx,gauss);
+		vy_input->GetParameterValue(&vy,gauss);
+		vxobs_input->GetParameterValue(&vxobs,gauss);
+		vyobs_input->GetParameterValue(&vyobs,gauss);
 		GetNodalFunctions(l1l2l3, gauss);
 
-		TriaRef::GetParameterValue(&dux, &dux_list[0],gauss);
-		TriaRef::GetParameterValue(&duy, &duy_list[0],gauss);
-
-		for (i=0;i<NUMVERTICES;i++){
-			pe->values[i*NDOF4+0]+=dux*Jdet*gauss->weight*l1l2l3[i];
-			pe->values[i*NDOF4+1]+=duy*Jdet*gauss->weight*l1l2l3[i]; 
+		if(response==SurfaceAbsVelMisfitEnum){
+			/*
+			 *      1  [           2              2 ]
+			 * J = --- | (u - u   )  +  (v - v   )  |
+			 *      2  [       obs            obs   ]
+			 *
+			 *        dJ
+			 * DU = - -- = (u   - u )
+			 *        du     obs
+			 */
+			for (i=0;i<NUMVERTICES;i++){
+				dux=vxobs-vx;
+				duy=vyobs-vy;
+				pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 
+				pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 
+			}
+		}
+		else if(response==SurfaceRelVelMisfitEnum){
+			/*
+			 *      1  [     \bar{v}^2             2   \bar{v}^2              2 ]
+			 * J = --- | -------------  (u - u   ) + -------------  (v - v   )  |
+			 *      2  [  (u   + eps)^2       obs    (v   + eps)^2       obs    ]
+			 *              obs                        obs                      
+			 *
+			 *        dJ     \bar{v}^2
+			 * DU = - -- = ------------- (u   - u )
+			 *        du   (u   + eps)^2    obs
+			 *               obs
+			 */
+			for (i=0;i<NUMVERTICES;i++){
+				scalex=pow(meanvel/(vxobs+epsvel),2.); if(vxobs==0)scalex=0;
+				scaley=pow(meanvel/(vyobs+epsvel),2.); if(vyobs==0)scaley=0;
+				dux=scalex*(vxobs-vx);
+				duy=scaley*(vyobs-vy);
+				pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 
+				pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 
+			}
+		}
+		else if(response==SurfaceLogVelMisfitEnum){
+			/*
+			 *                 [        vel + eps     ] 2
+			 * J = 4 \bar{v}^2 | log ( -----------  ) |  
+			 *                 [       vel   + eps    ]
+			 *                            obs
+			 *
+			 *        dJ                 2 * log(...)
+			 * DU = - -- = - 4 \bar{v}^2 -------------  u
+			 *        du                 vel^2 + eps
+			 *            
+			 */
+			for (i=0;i<NUMVERTICES;i++){
+				velocity_mag    =sqrt(pow(vx,   2.)+pow(vy,   2.))+epsvel;
+				obs_velocity_mag=sqrt(pow(vxobs,2.)+pow(vyobs,2.))+epsvel;
+				scale=-8*pow(meanvel,2.)/pow(velocity_mag,2.)*log(velocity_mag/obs_velocity_mag);
+				dux=scale*vx;
+				duy=scale*vy;
+				pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 
+				pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 
+			}
+		}
+		else if(response==SurfaceAverageVelMisfitEnum){
+			/*
+			 *      1                    2              2
+			 * J = ---  sqrt(  (u - u   )  +  (v - v   )  )
+			 *      S                obs            obs
+			 *
+			 *        dJ      1       1 
+			 * DU = - -- = - --- ----------- * 2 (u - u   )
+			 *        du      S  2 sqrt(...)           obs
+			 */
+			for (i=0;i<NUMVERTICES;i++){
+				scale=1./(S*2*sqrt(pow(vx-vxobs,2.)+pow(vy-vyobs,2.))+epsvel);
+				dux=scale*(vxobs-vx);
+				duy=scale*(vyobs-vy);
+				pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 
+				pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 
+			}
+		}
+		else if(response==SurfaceLogVxVyMisfitEnum){
+			/*
+			 *      1            [        |u| + eps     2          |v| + eps     2  ]
+			 * J = --- \bar{v}^2 | log ( -----------  )   +  log ( -----------  )   |  
+			 *      2            [       |u    |+ eps              |v    |+ eps     ]
+			 *                              obs                       obs
+			 *        dJ                              1      u                             1
+			 * DU = - -- = - \bar{v}^2 log(u...) --------- ----  ~ - \bar{v}^2 log(u...) ------
+			 *        du                         |u| + eps  |u|                           u + eps
+			 */
+			for (i=0;i<NUMVERTICES;i++){
+				dux = - pow(meanvel,2.) * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
+				duy = - pow(meanvel,2.) * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel);
+				pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 
+				pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 
+			}
+		}
+		else{
+			/*Not supported yet! : */
+			_error_("response %s not supported yet",EnumToStringx(response));
 		}
 	}
@@ -4845,14 +4793,7 @@
 	int        i,ig;
 	double     Jelem=0;
-	double     velocity_mag,obs_velocity_mag;
-	double     meanvel, epsvel,misfit,Jdet;
+	double     misfit,Jdet;
+	double     vx,vy,vxobs,vyobs,weight;
 	double     xyz_list[NUMVERTICES][3];
-	double     vx_list[NUMVERTICES];
-	double     vy_list[NUMVERTICES];
-	double     obs_vx_list[NUMVERTICES];
-	double     obs_vy_list[NUMVERTICES];
-	double     misfit_square_list[NUMVERTICES];
-	double     misfit_list[NUMVERTICES];
-	double     weights_list[NUMVERTICES];
 	GaussTria *gauss=NULL;
 
@@ -4863,50 +4804,40 @@
 	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 
-	/* Recover input data: */
-	GetParameterListOnVertices(&obs_vx_list[0],VxObsEnum);
-	GetParameterListOnVertices(&obs_vy_list[0],VyObsEnum);
-	GetParameterListOnVertices(&vx_list[0],VxEnum);
-	GetParameterListOnVertices(&vy_list[0],VyEnum);
-	GetParameterListOnVertices(&weights_list[0],WeightsEnum,0,0);
-
-	/*retrieve some parameters: */
-	this->parameters->FindParam(&meanvel,MeanVelEnum);
-	this->parameters->FindParam(&epsvel,EpsVelEnum);
-	
-	/* Compute SurfaceAbsVelMisfit at the 3 nodes
-	 * Here we integrate linearized functions:
-	 *               
-	 * J(E) = int_E   sum_{i=1}^3  J_i Phi_i
-	 *
-	 * where J_i are the misfits at the 3 nodes of the triangle
-	 *       Phi_i is the nodal function (P1) with respect to 
-	 *       the vertex i
-	 */
-
-	/*We are using an absolute misfit:
-	 *
-	 *      1  [           2              2 ]
-	 * J = --- | (u - u   )  +  (v - v   )  |
-	 *      2  [       obs            obs   ]
-	 *
-	 */
-	for (i=0;i<NUMVERTICES;i++){
-		misfit_list[i]=0.5*(pow((vx_list[i]-obs_vx_list[i]),(double)2)+pow((vy_list[i]-obs_vy_list[i]),(double)2));
-	}
-	/*Process units: */
-	if(process_units)UnitConversion(&misfit_list[0],NUMVERTICES,IuToExtEnum,SurfaceAbsVelMisfitEnum,this->parameters);
-
-	/*Apply weights to misfits*/
-	for (i=0;i<NUMVERTICES;i++) misfit_list[i]=weights_list[i]*misfit_list[i];
+	/*Retrieve all inputs we will be needing: */
+	Input* weights_input=inputs->GetInput(WeightsEnum);   _assert_(weights_input);
+	Input* vx_input     =inputs->GetInput(VxEnum);        _assert_(vx_input);
+	Input* vy_input     =inputs->GetInput(VyEnum);        _assert_(vy_input);
+	Input* vxobs_input  =inputs->GetInput(VxObsEnum);     _assert_(vxobs_input);
+	Input* vyobs_input  =inputs->GetInput(VyObsEnum);     _assert_(vyobs_input);
 
 	/* Start  looping on the number of gaussian points: */
 	gauss=new GaussTria(2);
-	for (ig=gauss->begin();ig<gauss->end(); ig++){
+	for (ig=gauss->begin();ig<gauss->end();ig++){
 
 		gauss->GaussPoint(ig);
 
+		/* Get Jacobian determinant: */
 		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
-		TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss);
-		Jelem+=misfit*Jdet*gauss->weight;
+
+		/*Get all parameters at gaussian point*/
+		weights_input->GetParameterValue(&weight,gauss,0);
+		vx_input->GetParameterValue(&vx,gauss);
+		vy_input->GetParameterValue(&vy,gauss);
+		vxobs_input->GetParameterValue(&vxobs,gauss);
+		vyobs_input->GetParameterValue(&vyobs,gauss);
+
+		/*Compute SurfaceAbsVelMisfitEnum:
+		 *
+		 *      1  [           2              2 ]
+		 * J = --- | (u - u   )  +  (v - v   )  |
+		 *      2  [       obs            obs   ]
+		 *
+		 */
+		misfit=0.5*( pow(vx-vxobs,2.) + pow(vy-vyobs,2.) );
+
+		if(process_units)UnitConversion(misfit,IuToExtEnum,SurfaceAverageVelMisfitEnum,this->parameters);
+
+		/*Add to cost function*/
+		Jelem+=misfit*weight*Jdet*gauss->weight;
 	}
 
@@ -4951,17 +4882,8 @@
 
 	int        i,ig;
-	int        fit=-1;
-	double     Jelem=0,S=0;
-	double     scalex=1, scaley=1;
-	double     meanvel, epsvel,Jdet;
-	double     velocity_mag,obs_velocity_mag,misfit;
+	double     Jelem=0,S,Jdet;
+	double     misfit;
+	double     vx,vy,vxobs,vyobs,weight;
 	double     xyz_list[NUMVERTICES][3];
-	double     vx_list[NUMVERTICES];
-	double     vy_list[NUMVERTICES];
-	double     obs_vx_list[NUMVERTICES];
-	double     obs_vy_list[NUMVERTICES];
-	double     misfit_square_list[NUMVERTICES];
-	double     misfit_list[NUMVERTICES];
-	double     weights_list[NUMVERTICES];
 	GaussTria *gauss=NULL;
 
@@ -4972,52 +4894,40 @@
 	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 
-	/* Recover input data: */
+	/*Retrieve all inputs we will be needing: */
 	inputs->GetParameterValue(&S,SurfaceAreaEnum);
-	GetParameterListOnVertices(&obs_vx_list[0],VxObsEnum);
-	GetParameterListOnVertices(&obs_vy_list[0],VyObsEnum);
-	GetParameterListOnVertices(&vx_list[0],VxEnum);
-	GetParameterListOnVertices(&vy_list[0],VyEnum);
-	GetParameterListOnVertices(&weights_list[0],WeightsEnum,0,0);
-
-	/*retrieve some parameters: */
-	this->parameters->FindParam(&meanvel,MeanVelEnum);
-	this->parameters->FindParam(&epsvel,EpsVelEnum);
-
-	/* Compute SurfaceAverageVelMisfit at the 3 nodes
-	 * Here we integrate linearized functions:
-	 *               
-	 * J(E) = int_E   sum_{i=1}^3  J_i Phi_i
-	 *
-	 * where J_i are the misfits at the 3 nodes of the triangle
-	 *       Phi_i is the nodal function (P1) with respect to 
-	 *       the vertex i
-	 */
-
-	/*We are using a spacially average absolute misfit:
-	 *
-	 *      1                    2              2
-	 * J = ---  sqrt(  (u - u   )  +  (v - v   )  )
-	 *      S                obs            obs
-	 */
-	for (i=0;i<NUMVERTICES;i++) misfit_square_list[i]=pow(vx_list[i]-obs_vx_list[i],2)+pow(vy_list[i]-obs_vx_list[i],2);
-
-	/*Process units: */
-	if(process_units)UnitConversion(&misfit_square_list[0],NUMVERTICES,IuToExtEnum,SurfaceAverageVelMisfitEnum,this->parameters);
-
-	/*Take the square root, and scale by surface: */
-	for (i=0;i<NUMVERTICES;i++)misfit_list[i]=pow(misfit_square_list[i],2)/S;
-
-	/*Apply weights to misfits*/
-	for (i=0;i<NUMVERTICES;i++) misfit_list[i]=weights_list[i]*misfit_list[i];
+	Input* weights_input=inputs->GetInput(WeightsEnum);   _assert_(weights_input);
+	Input* vx_input     =inputs->GetInput(VxEnum);        _assert_(vx_input);
+	Input* vy_input     =inputs->GetInput(VyEnum);        _assert_(vy_input);
+	Input* vxobs_input  =inputs->GetInput(VxObsEnum);     _assert_(vxobs_input);
+	Input* vyobs_input  =inputs->GetInput(VyObsEnum);     _assert_(vyobs_input);
 
 	/* Start  looping on the number of gaussian points: */
-	gauss=new GaussTria(2);
+	gauss=new GaussTria(3);
 	for (ig=gauss->begin();ig<gauss->end();ig++){
 
 		gauss->GaussPoint(ig);
 
+		/* Get Jacobian determinant: */
 		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
-		TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss);
-		Jelem+=misfit*Jdet*gauss->weight;
+
+		/*Get all parameters at gaussian point*/
+		weights_input->GetParameterValue(&weight,gauss,0);
+		vx_input->GetParameterValue(&vx,gauss);
+		vy_input->GetParameterValue(&vy,gauss);
+		vxobs_input->GetParameterValue(&vxobs,gauss);
+		vyobs_input->GetParameterValue(&vyobs,gauss);
+
+		/*Compute SurfaceAverageVelMisfitEnum:
+		 *
+		 *      1                    2              2
+		 * J = ---  sqrt(  (u - u   )  +  (v - v   )  )
+		 *      S                obs            obs
+		 */
+		misfit=1/S*pow( pow(vx-vxobs,2.) + pow(vy-vyobs,2.) ,0.5);
+
+		if(process_units)UnitConversion(misfit,IuToExtEnum,SurfaceAverageVelMisfitEnum,this->parameters);
+
+		/*Add to cost function*/
+		Jelem+=misfit*weight*Jdet*gauss->weight;
 	}
 
@@ -5031,4 +4941,158 @@
 
 	const int    numdof=NDOF2*NUMVERTICES;
+
+	int        i,ig;
+	double     Jelem=0;
+	double     meanvel, epsvel,misfit,Jdet;
+	double     velocity_mag,obs_velocity_mag;
+	double     xyz_list[NUMVERTICES][3];
+	double     vx,vy,vxobs,vyobs,weight;
+	GaussTria *gauss=NULL;
+
+	/*If on water, return 0: */
+	if(IsOnWater())return 0;
+
+	/* Get node coordinates and dof list: */
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
+
+	/*Retrieve all inputs we will be needing: */
+	this->parameters->FindParam(&meanvel,MeanVelEnum);
+	this->parameters->FindParam(&epsvel,EpsVelEnum);
+	Input* weights_input=inputs->GetInput(WeightsEnum);   _assert_(weights_input);
+	Input* vx_input     =inputs->GetInput(VxEnum);        _assert_(vx_input);
+	Input* vy_input     =inputs->GetInput(VyEnum);        _assert_(vy_input);
+	Input* vxobs_input  =inputs->GetInput(VxObsEnum);     _assert_(vxobs_input);
+	Input* vyobs_input  =inputs->GetInput(VyObsEnum);     _assert_(vyobs_input);
+
+	/* Start  looping on the number of gaussian points: */
+	gauss=new GaussTria(4);
+	for (ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+
+		/* Get Jacobian determinant: */
+		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
+
+		/*Get all parameters at gaussian point*/
+		weights_input->GetParameterValue(&weight,gauss,0);
+		vx_input->GetParameterValue(&vx,gauss);
+		vy_input->GetParameterValue(&vy,gauss);
+		vxobs_input->GetParameterValue(&vxobs,gauss);
+		vyobs_input->GetParameterValue(&vyobs,gauss);
+
+		/*Compute SurfaceLogVelMisfit:
+		 *                 [        vel + eps     ] 2
+		 * J = 4 \bar{v}^2 | log ( -----------  ) |  
+		 *                 [       vel   + eps    ]
+		 *                            obs
+		 */
+		velocity_mag    =sqrt(pow(vx,   2.)+pow(vy,   2.))+epsvel;
+		obs_velocity_mag=sqrt(pow(vxobs,2.)+pow(vyobs,2.))+epsvel;
+		misfit=4*pow(meanvel,2.)*pow(log(velocity_mag/obs_velocity_mag),2.);
+
+		if(process_units)UnitConversion(misfit,IuToExtEnum,SurfaceLogVelMisfitEnum,this->parameters);
+
+		/*Add to cost function*/
+		Jelem+=misfit*weight*Jdet*gauss->weight;
+	}
+
+	/*clean-up and Return: */
+	delete gauss;
+	return Jelem;
+}
+/*}}}*/
+/*FUNCTION Tria::SurfaceLogVxVyMisfit {{{1*/
+double Tria::SurfaceLogVxVyMisfit(bool process_units){
+
+	const int    numdof=NDOF2*NUMVERTICES;
+
+	int        i,ig;
+	int        fit=-1;
+	double     Jelem=0, S=0;
+	double     meanvel, epsvel, misfit, Jdet;
+	double     vx,vy,vxobs,vyobs,weight;
+	double     xyz_list[NUMVERTICES][3];
+	GaussTria *gauss=NULL;
+
+	/*If on water, return 0: */
+	if(IsOnWater())return 0;
+
+	/* Get node coordinates and dof list: */
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
+
+	/*Retrieve all inputs we will be needing: */
+	this->parameters->FindParam(&meanvel,MeanVelEnum);
+	this->parameters->FindParam(&epsvel,EpsVelEnum);
+	Input* weights_input=inputs->GetInput(WeightsEnum);   _assert_(weights_input);
+	Input* vx_input     =inputs->GetInput(VxEnum);        _assert_(vx_input);
+	Input* vy_input     =inputs->GetInput(VyEnum);        _assert_(vy_input);
+	Input* vxobs_input  =inputs->GetInput(VxObsEnum);     _assert_(vxobs_input);
+	Input* vyobs_input  =inputs->GetInput(VyObsEnum);     _assert_(vyobs_input);
+	
+	/* Start  looping on the number of gaussian points: */
+	gauss=new GaussTria(4);
+	for (ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+
+		/* Get Jacobian determinant: */
+		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
+
+		/*Get all parameters at gaussian point*/
+		weights_input->GetParameterValue(&weight,gauss,0);
+		vx_input->GetParameterValue(&vx,gauss);
+		vy_input->GetParameterValue(&vy,gauss);
+		vxobs_input->GetParameterValue(&vxobs,gauss);
+		vyobs_input->GetParameterValue(&vyobs,gauss);
+
+		/*Compute SurfaceRelVelMisfit:
+		 *
+		 *      1            [        |u| + eps     2          |v| + eps     2  ]
+		 * J = --- \bar{v}^2 | log ( -----------  )   +  log ( -----------  )   |  
+		 *      2            [       |u    |+ eps              |v    |+ eps     ]
+		 *                              obs                       obs
+		 */
+		misfit=0.5*pow(meanvel,2.)*(
+					pow(log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)),2.) +
+					pow(log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)),2.) );
+
+		if(process_units)UnitConversion(misfit,IuToExtEnum,SurfaceLogVxVyMisfitEnum,this->parameters);
+
+		/*Add to cost function*/
+		Jelem+=misfit*weight*Jdet*gauss->weight;
+	}
+
+	/*clean-up and Return: */
+	delete gauss;
+	return Jelem;
+}
+/*}}}*/
+/*FUNCTION Tria::SurfaceNormal{{{1*/
+void Tria::SurfaceNormal(double* surface_normal, double xyz_list[3][3]){
+
+	int i;
+	double v13[3],v23[3];
+	double normal[3];
+	double normal_norm;
+
+	for (i=0;i<3;i++){
+		v13[i]=xyz_list[0][i]-xyz_list[2][i];
+		v23[i]=xyz_list[1][i]-xyz_list[2][i];
+	}
+
+	normal[0]=v13[1]*v23[2]-v13[2]*v23[1];
+	normal[1]=v13[2]*v23[0]-v13[0]*v23[2];
+	normal[2]=v13[0]*v23[1]-v13[1]*v23[0];
+
+	normal_norm=sqrt( pow(normal[0],(double)2)+pow(normal[1],(double)2)+pow(normal[2],(double)2) );
+
+	*(surface_normal)=normal[0]/normal_norm;
+	*(surface_normal+1)=normal[1]/normal_norm;
+	*(surface_normal+2)=normal[2]/normal_norm;
+}
+/*}}}*/
+/*FUNCTION Tria::SurfaceRelVelMisfit {{{1*/
+double Tria::SurfaceRelVelMisfit(bool process_units){
+	const int  numdof=2*NUMVERTICES;
 
 	int        i,ig;
@@ -5036,13 +5100,6 @@
 	double     scalex=1,scaley=1;
 	double     meanvel, epsvel,misfit,Jdet;
-	double     velocity_mag,obs_velocity_mag;
+	double     vx,vy,vxobs,vyobs,weight;
 	double     xyz_list[NUMVERTICES][3];
-	double     vx_list[NUMVERTICES];
-	double     vy_list[NUMVERTICES];
-	double     obs_vx_list[NUMVERTICES];
-	double     obs_vy_list[NUMVERTICES];
-	double     misfit_square_list[NUMVERTICES];
-	double     misfit_list[NUMVERTICES];
-	double     weights_list[NUMVERTICES];
 	GaussTria *gauss=NULL;
 
@@ -5053,244 +5110,43 @@
 	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 
-	/* Recover input data: */
-	GetParameterListOnVertices(&obs_vx_list[0],VxObsEnum);
-	GetParameterListOnVertices(&obs_vy_list[0],VyObsEnum);
-	GetParameterListOnVertices(&vx_list[0],VxEnum);
-	GetParameterListOnVertices(&vy_list[0],VyEnum);
-	GetParameterListOnVertices(&weights_list[0],WeightsEnum,0,0);
-
-	/*retrieve some parameters: */
+	/*Retrieve all inputs we will be needing: */
 	this->parameters->FindParam(&meanvel,MeanVelEnum);
 	this->parameters->FindParam(&epsvel,EpsVelEnum);
-	
-	/* Compute SurfaceLogVelMisfit at the 3 nodes
-	 * Here we integrate linearized functions:
-	 *               
-	 * J(E) = int_E   sum_{i=1}^3  J_i Phi_i
-	 *
-	 * where J_i are the misfits at the 3 nodes of the triangle
-	 *       Phi_i is the nodal function (P1) with respect to 
-	 *       the vertex i
-	 */
-
-	/*We are using a logarithmic misfit:
-	 *                        
-	 *                 [        vel + eps     ] 2
-	 * J = 4 \bar{v}^2 | log ( -----------  ) |  
-	 *                 [       vel   + eps    ]
-	 *                            obs
-	 */
-	for (i=0;i<NUMVERTICES;i++){
-		velocity_mag=sqrt(pow(vx_list[i],(double)2)+pow(vy_list[i],(double)2))+epsvel; //epsvel to avoid velocity being nil.
-		obs_velocity_mag=sqrt(pow(obs_vx_list[i],(double)2)+pow(obs_vy_list[i],(double)2))+epsvel; //epsvel to avoid observed velocity being nil.
-		misfit_list[i]=4*pow(meanvel,(double)2)*pow(log(velocity_mag/obs_velocity_mag),(double)2);
-	}
-
-	/*Process units: */
-	if(process_units)UnitConversion(&misfit_list[0],NUMVERTICES,IuToExtEnum,SurfaceLogVelMisfitEnum,this->parameters);
-
-	/*Apply weights to misfits*/
-	for (i=0;i<NUMVERTICES;i++) misfit_list[i]=weights_list[i]*misfit_list[i];
+	Input* weights_input=inputs->GetInput(WeightsEnum);   _assert_(weights_input);
+	Input* vx_input     =inputs->GetInput(VxEnum);        _assert_(vx_input);
+	Input* vy_input     =inputs->GetInput(VyEnum);        _assert_(vy_input);
+	Input* vxobs_input  =inputs->GetInput(VxObsEnum);     _assert_(vxobs_input);
+	Input* vyobs_input  =inputs->GetInput(VyObsEnum);     _assert_(vyobs_input);
 
 	/* Start  looping on the number of gaussian points: */
-	gauss=new GaussTria(2);
+	gauss=new GaussTria(4);
 	for (ig=gauss->begin();ig<gauss->end();ig++){
 
 		gauss->GaussPoint(ig);
 
+		/* Get Jacobian determinant: */
 		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
-		TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss);
-		Jelem+=misfit*Jdet*gauss->weight;
-	}
-
-	/*clean-up and Return: */
-	delete gauss;
-	return Jelem;
-}
-/*}}}*/
-/*FUNCTION Tria::SurfaceLogVxVyMisfit {{{1*/
-double Tria::SurfaceLogVxVyMisfit(bool process_units){
-
-	const int    numdof=NDOF2*NUMVERTICES;
-
-	int        i,ig;
-	int        fit=-1;
-	double     Jelem=0, S=0;
-	double     scalex=1,scaley=1;
-	double     meanvel, epsvel, misfit, Jdet;
-	double     velocity_mag,obs_velocity_mag;
-	double     xyz_list[NUMVERTICES][3];
-	double     vx_list[NUMVERTICES];
-	double     vy_list[NUMVERTICES];
-	double     obs_vx_list[NUMVERTICES];
-	double     obs_vy_list[NUMVERTICES];
-	double     misfit_square_list[NUMVERTICES];
-	double     misfit_list[NUMVERTICES];
-	double     weights_list[NUMVERTICES];
-	GaussTria *gauss=NULL;
-
-	/*If on water, return 0: */
-	if(IsOnWater())return 0;
-
-	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
-
-	/* Recover input data: */
-	GetParameterListOnVertices(&obs_vx_list[0],VxObsEnum);
-	GetParameterListOnVertices(&obs_vy_list[0],VyObsEnum);
-	GetParameterListOnVertices(&vx_list[0],VxEnum);
-	GetParameterListOnVertices(&vy_list[0],VyEnum);
-	GetParameterListOnVertices(&weights_list[0],WeightsEnum,0,0);
-
-	/*retrieve some parameters: */
-	this->parameters->FindParam(&meanvel,MeanVelEnum);
-	this->parameters->FindParam(&epsvel,EpsVelEnum);
-	
-	/* Compute SurfaceLogVxVyMisfit at the 3 nodes
-	 * Here we integrate linearized functions:
-	 *               
-	 * J(E) = int_E   sum_{i=1}^3  J_i Phi_i
-	 *
-	 * where J_i are the misfits at the 3 nodes of the triangle
-	 *       Phi_i is the nodal function (P1) with respect to 
-	 *       the vertex i
-	 */
-
-	/*We are using an logarithmic 2 misfit:
-	 *
-	 *      1            [        |u| + eps     2          |v| + eps     2  ]
-	 * J = --- \bar{v}^2 | log ( -----------  )   +  log ( -----------  )   |  
-	 *      2            [       |u    |+ eps              |v    |+ eps     ]
-	 *                              obs                       obs
-	 */
-	for (i=0;i<NUMVERTICES;i++){
-		misfit_list[i]=0.5*pow(meanvel,(double)2)*(
-					pow(log((fabs(vx_list[i])+epsvel)/(fabs(obs_vx_list[i])+epsvel)),(double)2) +
-					pow(log((fabs(vy_list[i])+epsvel)/(fabs(obs_vy_list[i])+epsvel)),(double)2) );
-	}
-
-	/*Process units: */
-	if(process_units)UnitConversion(&misfit_list[0],NUMVERTICES,IuToExtEnum,SurfaceLogVxVyMisfitEnum,this->parameters);
-
-	/*Apply weights to misfits*/
-	for (i=0;i<NUMVERTICES;i++) misfit_list[i]=weights_list[i]*misfit_list[i];
-
-	/* Start  looping on the number of gaussian points: */
-	gauss=new GaussTria(2);
-	for (ig=gauss->begin();ig<gauss->end();ig++){
-
-		gauss->GaussPoint(ig);
-
-		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
-		TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss);
-		Jelem+=misfit*Jdet*gauss->weight;
-	}
-
-	/*clean-up and Return: */
-	delete gauss;
-	return Jelem;
-}
-/*}}}*/
-/*FUNCTION Tria::SurfaceNormal{{{1*/
-void Tria::SurfaceNormal(double* surface_normal, double xyz_list[3][3]){
-
-	int i;
-	double v13[3],v23[3];
-	double normal[3];
-	double normal_norm;
-
-	for (i=0;i<3;i++){
-		v13[i]=xyz_list[0][i]-xyz_list[2][i];
-		v23[i]=xyz_list[1][i]-xyz_list[2][i];
-	}
-
-	normal[0]=v13[1]*v23[2]-v13[2]*v23[1];
-	normal[1]=v13[2]*v23[0]-v13[0]*v23[2];
-	normal[2]=v13[0]*v23[1]-v13[1]*v23[0];
-
-	normal_norm=sqrt( pow(normal[0],(double)2)+pow(normal[1],(double)2)+pow(normal[2],(double)2) );
-
-	*(surface_normal)=normal[0]/normal_norm;
-	*(surface_normal+1)=normal[1]/normal_norm;
-	*(surface_normal+2)=normal[2]/normal_norm;
-}
-/*}}}*/
-/*FUNCTION Tria::SurfaceRelVelMisfit {{{1*/
-double Tria::SurfaceRelVelMisfit(bool process_units){
-
-	const int    numdof=2*NUMVERTICES;
-
-	int        i,ig;
-	double     Jelem=0;
-	double     scalex=1,scaley=1;
-	double     meanvel, epsvel,misfit,Jdet;
-	double     velocity_mag,obs_velocity_mag;
-	double     xyz_list[NUMVERTICES][3];
-	double     vx_list[NUMVERTICES];
-	double     vy_list[NUMVERTICES];
-	double     obs_vx_list[NUMVERTICES];
-	double     obs_vy_list[NUMVERTICES];
-	double     misfit_square_list[NUMVERTICES];
-	double     misfit_list[NUMVERTICES];
-	double     weights_list[NUMVERTICES];
-	GaussTria *gauss=NULL;
-
-	/*If on water, return 0: */
-	if(IsOnWater())return 0;
-
-	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
-
-	/* Recover input data: */
-	GetParameterListOnVertices(&weights_list[0],WeightsEnum,0,0);
-	GetParameterListOnVertices(&obs_vx_list[0],VxObsEnum);
-	GetParameterListOnVertices(&obs_vy_list[0],VyObsEnum);
-	GetParameterListOnVertices(&vx_list[0],VxEnum);
-	GetParameterListOnVertices(&vy_list[0],VyEnum);
-
-	/*retrieve some parameters: */
-	this->parameters->FindParam(&meanvel,MeanVelEnum);
-	this->parameters->FindParam(&epsvel,EpsVelEnum);
-
-	/* Compute SurfaceRelVelMisfit at the 3 nodes
-	 * Here we integrate linearized functions:
-	 *               
-	 * J(E) = int_E   sum_{i=1}^3  J_i Phi_i
-	 *
-	 * where J_i are the misfits at the 3 nodes of the triangle
-	 *       Phi_i is the nodal function (P1) with respect to 
-	 *       the vertex i
-	 */
-
-	/*We are using a relative misfit: 
-	 *                        
-	 *      1  [     \bar{v}^2             2   \bar{v}^2              2 ]
-	 * J = --- | -------------  (u - u   ) + -------------  (v - v   )  |
-	 *      2  [  (u   + eps)^2       obs    (v   + eps)^2       obs    ]
-	 *              obs                        obs                      
-	 */
-	for (i=0;i<NUMVERTICES;i++){
-		scalex=pow(meanvel/(obs_vx_list[i]+epsvel),(double)2);
-		scaley=pow(meanvel/(obs_vy_list[i]+epsvel),(double)2);
-		if(obs_vx_list[i]==0)scalex=0;
-		if(obs_vy_list[i]==0)scaley=0;
-		misfit_list[i]=0.5*(scalex*pow((vx_list[i]-obs_vx_list[i]),2)+scaley*pow((vy_list[i]-obs_vy_list[i]),2));
-	}
-
-	/*Process units: */
-	if(process_units)UnitConversion(&misfit_list[0],NUMVERTICES,IuToExtEnum,SurfaceRelVelMisfitEnum,this->parameters);
-
-	/*Apply weights to misfits*/
-	for (i=0;i<NUMVERTICES;i++) misfit_list[i]=weights_list[i]*misfit_list[i];
-
-	/* Start  looping on the number of gaussian points: */
-	gauss=new GaussTria(2);
-	for (ig=gauss->begin();ig<gauss->end();ig++){
-
-		gauss->GaussPoint(ig);
-
-		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
-		TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss);
-		Jelem+=misfit*Jdet*gauss->weight;
+
+		/*Get all parameters at gaussian point*/
+		weights_input->GetParameterValue(&weight,gauss,0);
+		vx_input->GetParameterValue(&vx,gauss);
+		vy_input->GetParameterValue(&vy,gauss);
+		vxobs_input->GetParameterValue(&vxobs,gauss);
+		vyobs_input->GetParameterValue(&vyobs,gauss);
+
+		/*Compute SurfaceRelVelMisfit:
+		 *                        
+		 *      1  [     \bar{v}^2             2   \bar{v}^2              2 ]
+		 * J = --- | -------------  (u - u   ) + -------------  (v - v   )  |
+		 *      2  [  (u   + eps)^2       obs    (v   + eps)^2       obs    ]
+		 *              obs                        obs                      
+		 */
+		scalex=pow(meanvel/(vxobs+epsvel),2.); if(vxobs==0)scalex=0;
+		scaley=pow(meanvel/(vyobs+epsvel),2.); if(vyobs==0)scaley=0;
+		misfit=0.5*(scalex*pow((vx-vxobs),2.)+scaley*pow((vy-vyobs),2.));
+		if(process_units)UnitConversion(misfit,IuToExtEnum,SurfaceRelVelMisfitEnum,this->parameters);
+
+		/*Add to cost function*/
+		Jelem+=misfit*weight*Jdet*gauss->weight;
 	}