Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 15717)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 15718)
@@ -5187,187 +5187,4 @@
 }
 /*}}}*/
-/*FUNCTION Tria::CreatePVectorAdjointFS{{{*/
-ElementVector* Tria::CreatePVectorAdjointFS(void){
-
-	/*Intermediaries */
-	int        i,resp;
-	int       *responses=NULL;
-	int        num_responses;
-	IssmDouble     Jdet;
-	IssmDouble     obs_velocity_mag,velocity_mag;
-	IssmDouble     dux,duy;
-	IssmDouble     epsvel=2.220446049250313e-16;
-	IssmDouble     meanvel=3.170979198376458e-05; /*1000 m/yr*/
-	IssmDouble     scalex=0,scaley=0,scale=0,S=0;
-	IssmDouble     vx,vy,vxobs,vyobs,weight;
-	IssmDouble     xyz_list[NUMVERTICES][3];
-	IssmDouble     basis[3];
-	GaussTria* gauss=NULL;
-
-	/*Initialize Element vector*/
-	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,FSApproximationEnum);
-
-	/*Retrieve all inputs and parameters*/
-	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-	this->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
-	this->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum);
-	Input* weights_input = inputs->GetInput(InversionCostFunctionsCoefficientsEnum);   _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(InversionVxObsEnum);     _assert_(vxobs_input);
-	Input* vyobs_input   = inputs->GetInput(InversionVyObsEnum);     _assert_(vyobs_input);
-
-	/*Get Surface if required by one response*/
-	for(resp=0;resp<num_responses;resp++){
-		if(responses[resp]==SurfaceAverageVelMisfitEnum){
-			inputs->GetInputValue(&S,SurfaceAreaEnum); break;
-		}
-	}
-
-	/* Start  looping on the number of gaussian points: */
-	gauss=new GaussTria(4);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-
-		gauss->GaussPoint(ig);
-
-		/* Get Jacobian determinant: */
-		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
-
-		/*Get all parameters at gaussian point*/
-		vx_input->GetInputValue(&vx,gauss);
-		vy_input->GetInputValue(&vy,gauss);
-		vxobs_input->GetInputValue(&vxobs,gauss);
-		vyobs_input->GetInputValue(&vyobs,gauss);
-		GetNodalFunctions(basis, gauss);
-
-		/*Loop over all requested responses*/
-		for(resp=0;resp<num_responses;resp++){
-
-			weights_input->GetInputValue(&weight,gauss,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*basis[i]; 
-						pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*basis[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*basis[i]; 
-						pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*basis[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*basis[i]; 
-						pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*basis[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*basis[i]; 
-						pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*basis[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 = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
-						duy = - meanvel*meanvel * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel);
-						pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
-						pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
-					}
-					break;
-				case DragCoefficientAbsGradientEnum:
-					/*Nothing in P vector*/
-					break;
-				case ThicknessAbsGradientEnum:
-					/*Nothing in P vector*/
-					break;
-				case ThicknessAcrossGradientEnum:
-					/*Nothing in P vector*/
-					break;
-				case ThicknessAlongGradientEnum:
-					/*Nothing in P vector*/
-					break;
-				case RheologyBbarAbsGradientEnum:
-					/*Nothing in P vector*/
-					break;
-				default:
-					_error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
-			}
-		}
-	}
-
-	/*Clean up and return*/
-	delete gauss;
-	xDelete<int>(responses);
-	return pe;
-}
-/*}}}*/
 /*FUNCTION Tria::DragCoefficientAbsGradient{{{*/
 IssmDouble Tria::DragCoefficientAbsGradient(int weight_index){
@@ -6050,7 +5867,4 @@
 ElementVector* Tria::CreatePVectorHydrologyDCInefficient(void){
 
-	/*Constants*/
-	const int    numdof=NDOF1*NUMVERTICES;
-
 	/*Intermediaries */
 	IssmDouble Jdet;
@@ -6059,9 +5873,11 @@
 	IssmDouble water_load,transfer;
 	IssmDouble sediment_storing;
-	IssmDouble basis[numdof];
-	GaussTria* gauss=NULL;
-
-	/*Initialize Element vector*/
-	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = this->NumberofNodes();
+
+	/*Initialize Element vector and other vectors*/
+	ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters);
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
 
 	/*Retrieve all inputs and parameters*/
@@ -6078,5 +5894,5 @@
 
 	/* Start  looping on the number of gaussian points: */
-	gauss=new GaussTria(2);
+	GaussTria* gauss=new GaussTria(2);
 	for(int ig=gauss->begin();ig<gauss->end();ig++){
 
@@ -6090,5 +5906,5 @@
 		scalar = Jdet*gauss->weight*(water_load+transfer);
 		if(reCast<bool,IssmDouble>(dt)) scalar = scalar*dt;
-		for(int i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
+		for(int i=0;i<numnodes;i++) pe->values[i]+=scalar*basis[i];
 
 		/*Transient term*/
@@ -6096,8 +5912,10 @@
 			old_wh_input->GetInputValue(&water_head,gauss);
 			scalar = Jdet*gauss->weight*water_head*sediment_storing;
-			for(int i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
-		}
-	}
+			for(int i=0;i<numnodes;i++) pe->values[i]+=scalar*basis[i];
+		}
+	}
+
 	/*Clean up and return*/
+	xDelete<IssmDouble>(basis);
 	delete gauss;
 	return pe;
@@ -6106,7 +5924,4 @@
 /*FUNCTION Tria::CreatePVectorHydrologyDCEfficient {{{*/
 ElementVector* Tria::CreatePVectorHydrologyDCEfficient(void){
-
-	/*Constants*/
-	const int    numdof=NDOF1*NUMVERTICES;
 
 	/*Intermediaries */
@@ -6117,6 +5932,5 @@
 	IssmDouble transfer,residual;
 	IssmDouble epl_storing;
-	IssmDouble basis[numdof];
-	GaussTria* gauss=NULL;
+	GaussTria* gauss = NULL;
 
 	/*Check that all nodes are active, else return empty matrix*/
@@ -6124,6 +5938,11 @@
 		return NULL;
 	}
-	/*Initialize Element vector*/
-	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = this->NumberofNodes();
+
+	/*Initialize Element vector and other vectors*/
+	ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters);
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
 
 	/*Retrieve all inputs and parameters*/
@@ -6145,5 +5964,5 @@
 		gauss->GaussPoint(ig);
 		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
-		GetNodalFunctions(basis, gauss);
+		GetNodalFunctions(basis,gauss);
 
 		/*Loading term*/
@@ -6151,5 +5970,5 @@
 		scalar = Jdet*gauss->weight*(-transfer);
 		if(reCast<bool,IssmDouble>(dt)) scalar = scalar*dt;
-		for(int i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
+		for(int i=0;i<numnodes;i++) pe->values[i]+=scalar*basis[i];
 
 		/*Transient term*/
@@ -6157,5 +5976,5 @@
 			old_wh_input->GetInputValue(&water_head,gauss);
 			scalar = Jdet*gauss->weight*water_head*epl_storing;
-			for(int i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
+			for(int i=0;i<numnodes;i++) pe->values[i]+=scalar*basis[i];
 		}
 	}
@@ -6171,5 +5990,7 @@
 		pe->values[iv]+=residual/connectivity;
 	}
+
 	/*Clean up and return*/
+	xDelete<IssmDouble>(basis);
 	delete gauss;
 	return pe;
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 15717)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 15718)
@@ -241,5 +241,4 @@
 		ElementMatrix* CreateKMatrixAdjointSSA(void);
 		ElementVector* CreatePVectorAdjointHoriz(void);
-		ElementVector* CreatePVectorAdjointFS(void);
 		ElementVector* CreatePVectorAdjointBalancethickness(void);
 		void	  InputUpdateFromSolutionAdjointBalancethickness( IssmDouble* solution);