Index: /issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 8602)
+++ /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 8603)
@@ -1599,7 +1599,9 @@
 
 	/*Intermediaries */
-	int         i,ig;
+	int         i,ig,resp;
 	double      Jdet;
 	double      thickness,thicknessobs,weight;
+	int        *responses = NULL;
+	int         num_responses;
 	double      xyz_list[NUMVERTICES][3];
 	double      l1l2l3[3];
@@ -1613,7 +1615,9 @@
 	/*Retrieve all inputs and parameters*/
 	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
-	Input* thickness_input   =inputs->GetInput(ThicknessEnum);   _assert_(thickness_input);
-	Input* thicknessobs_input=inputs->GetInput(ThicknessObsEnum);_assert_(thicknessobs_input);
-	Input* weights_input     =inputs->GetInput(WeightsEnum);     _assert_(weights_input);
+	this->parameters->FindParam(&num_responses,NumResponsesEnum);
+	this->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum);
+	Input* thickness_input    = inputs->GetInput(ThicknessEnum);   _assert_(thickness_input);
+	Input* thicknessobs_input = inputs->GetInput(ThicknessObsEnum);_assert_(thicknessobs_input);
+	Input* weights_input      = inputs->GetInput(WeightsEnum);     _assert_(weights_input);
 
 	/* Start  looping on the number of gaussian points: */
@@ -1632,8 +1636,16 @@
 		weights_input->GetParameterValue(&weight, gauss,0);
 
-		for(i=0;i<numdof;i++) pe->values[i]+=(thicknessobs-thickness)*weight*Jdet*gauss->weight*l1l2l3[i];
-		/*Regularization of the constraint: 2000000 79 N*/
-		//for(i=0;i<numdof;i++) pe->values[i]+= - 1*100000*dH[0]*dbasis[0][i]*Jdet*gauss->weight;
-		//for(i=0;i<numdof;i++) pe->values[i]+= - 1*100000*dH[1]*dbasis[1][i]*Jdet*gauss->weight;
+		/*Loop over all requested responses*/
+		for(resp=0;resp<num_responses;resp++) switch(responses[resp]){
+
+			case ThicknessAbsMisfitEnum:
+				for(i=0;i<numdof;i++) pe->values[i]+=(thicknessobs-thickness)*weight*Jdet*gauss->weight*l1l2l3[i];
+				/*Regularization of the constraint: 2000000 79 N*/
+				//for(i=0;i<numdof;i++) pe->values[i]+= - 1*100000*dH[0]*dbasis[0][i]*Jdet*gauss->weight;
+				//for(i=0;i<numdof;i++) pe->values[i]+= - 1*100000*dH[1]*dbasis[1][i]*Jdet*gauss->weight;
+				break;
+			default:
+				_error_("response %s not supported yet",EnumToStringx(responses[resp]));
+		}
 	}
 
@@ -1650,7 +1662,7 @@
 
 	/*Intermediaries */
-	int        i,ig;
+	int        i,resp,ig;
 	int       *responses=NULL;
-	int        response;
+	int        num_responses;
 	double     Jdet;
 	double     obs_velocity_mag,velocity_mag;
@@ -1667,7 +1679,8 @@
 	/*Retrieve all inputs and parameters*/
 	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
+	this->parameters->FindParam(&num_responses,NumResponsesEnum);
+	this->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum);
 	this->parameters->FindParam(&meanvel,MeanVelEnum);
 	this->parameters->FindParam(&epsvel,EpsVelEnum);
-	this->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum); response=responses[0];
 	Input* weights_input=inputs->GetInput(WeightsEnum);   _assert_(weights_input);
 	Input* vx_input     =inputs->GetInput(VxEnum);        _assert_(vx_input);
@@ -1676,5 +1689,10 @@
 	Input* vyobs_input  =inputs->GetInput(VyObsEnum);     _assert_(vyobs_input);
 
-	if(response==SurfaceAverageVelMisfitEnum) inputs->GetParameterValue(&S,SurfaceAreaEnum);
+	/*Get Surface if required by one response*/
+	for(resp=0;resp<num_responses;resp++){
+		if(responses[resp]==SurfaceAverageVelMisfitEnum){
+			inputs->GetParameterValue(&S,SurfaceAreaEnum); break;
+		}
+	}
 
 	/* Start  looping on the number of gaussian points: */
@@ -1695,102 +1713,104 @@
 		GetNodalFunctions(l1l2l3, gauss);
 
-		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));
+		/*Loop over all requested responses*/
+		for(resp=0;resp<num_responses;resp++) switch(responses[resp]){
+
+			case 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]; 
+				}
+				break;
+			case 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]; 
+				}
+				break;
+			case 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]; 
+				}
+				break;
+			case 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]; 
+				}
+				break;
+			case 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]; 
+				}
+				break;
+			default:
+				_error_("response %s not supported yet",EnumToStringx(responses[resp]));
 		}
 	}
@@ -1806,7 +1826,7 @@
 
 	/*Intermediaries */
-	int        i,ig;
+	int        i,resp,ig;
 	int       *responses=NULL;
-	int        response;
+	int        num_responses;
 	double     Jdet;
 	double     obs_velocity_mag,velocity_mag;
@@ -1823,14 +1843,20 @@
 	/*Retrieve all inputs and parameters*/
 	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
+	this->parameters->FindParam(&num_responses,NumResponsesEnum);
+	this->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum);
 	this->parameters->FindParam(&meanvel,MeanVelEnum);
 	this->parameters->FindParam(&epsvel,EpsVelEnum);
-	this->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum); response=responses[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);
-
-	if(response==SurfaceAverageVelMisfitEnum) inputs->GetParameterValue(&S,SurfaceAreaEnum);
+	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);
+
+	/*Get Surface if required by one response*/
+	for(resp=0;resp<num_responses;resp++){
+		if(responses[resp]==SurfaceAverageVelMisfitEnum){
+			inputs->GetParameterValue(&S,SurfaceAreaEnum); break;
+		}
+	}
 
 	/* Start  looping on the number of gaussian points: */
@@ -1851,102 +1877,104 @@
 		GetNodalFunctions(l1l2l3, gauss);
 
-		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));
+		/*Loop over all requested responses*/
+		for(resp=0;resp<num_responses;resp++) switch(responses[resp]){
+
+			case 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]; 
+				}
+				break;
+			case 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]; 
+				}
+				break;
+			case 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]; 
+				}
+				break;
+			case 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]; 
+				}
+				break;
+			case 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]; 
+				}
+				break;
+			default:
+				_error_("response %s not supported yet",EnumToStringx(responses[resp]));
 		}
 	}