Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 8645)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 8646)
@@ -3876,4 +3876,5 @@
 /*FUNCTION Penta::Gradj {{{1*/
 void  Penta::Gradj(Vec gradient,int control_type){
+	/*dJ/dalpha = ∂L/∂alpha = ∂J/∂alpha + ∂/∂alpha(KU-F)*/
 
 	int              i,approximation;
@@ -3883,4 +3884,5 @@
 	if(IsOnWater())return;
 
+	/*First deal with ∂/∂alpha(KU-F)*/
 	switch(control_type){
 
@@ -3962,4 +3964,30 @@
 		default:
 			_error_("control type %s not supported yet: ",EnumToStringx(control_type));
+	}
+
+	/*Now deal with ∂J/∂alpha*/
+	int        *responses = NULL;
+	int         num_responses,resp;
+	this->parameters->FindParam(&num_responses,NumResponsesEnum);
+	this->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum);
+
+	for(resp=0;resp<num_responses;resp++) switch(responses[resp]){
+
+		case ThicknessAbsMisfitEnum:
+		case SurfaceAbsVelMisfitEnum:
+		case SurfaceRelVelMisfitEnum:
+		case SurfaceLogVelMisfitEnum:
+		case SurfaceLogVxVyMisfitEnum:
+		case SurfaceAverageVelMisfitEnum:
+			/*Nothing, J does not depends on the parameter being inverted for*/
+			break;
+		case DragCoefficientAbsGradientEnum:
+			if (!IsOnBed()) return;
+			tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
+			tria->GradjDragGradient(gradient,resp);
+			delete tria->matice; delete tria;
+			break;
+		default:
+			_error_("response %s not supported yet",EnumToStringx(responses[resp]));
 	}
 }
@@ -5046,6 +5074,6 @@
 		/*Create inputs and add to DataSetInput*/
 		DatasetInput* datasetinput=new DatasetInput(WeightsEnum);
-		for(j=0;j<iomodel->num_cm_responses;j++){
-			for(i=0;i<6;i++)nodeinputs[i]=iomodel->weights[penta_vertex_ids[i]-1];
+		for(i=0;i<iomodel->num_cm_responses;i++){
+			for(j=0;j<6;j++)nodeinputs[j]=iomodel->weights[(penta_vertex_ids[j]-1)*iomodel->num_cm_responses+i];
 			datasetinput->inputs->AddObject(new PentaVertexInput(WeightsEnum,nodeinputs));
 		}
Index: /issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 8645)
+++ /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 8646)
@@ -1990,4 +1990,13 @@
 					}
 					break;
+				case DragCoefficientAbsGradientEnum:
+					/*Nothing in P vector*/
+					break;
+				case ThicknessAbsGradientEnum:
+					/*Nothing in P vector*/
+					break;
+				case RheologyBbarAbsGradientEnum:
+					/*Nothing in P vector*/
+					break;
 				default:
 					_error_("response %s not supported yet",EnumToStringx(responses[resp]));
@@ -2876,8 +2885,10 @@
 /*FUNCTION Tria::Gradj {{{1*/
 void  Tria::Gradj(Vec gradient,int control_type){
+	/*dJ/dalpha = ∂L/∂alpha = ∂J/∂alpha + ∂/∂alpha(KU-F)*/
 
 	/*If on water, grad = 0: */
 	if(IsOnWater()) return;
 
+	/*First deal with ∂/∂alpha(KU-F)*/
 	switch(control_type){
 		case DragCoefficientEnum:
@@ -2899,4 +2910,27 @@
 			_error_("%s%i","control type not supported yet: ",control_type);
 	}
+
+	/*Now deal with ∂J/∂alpha*/
+	int        *responses = NULL;
+	int         num_responses,resp;
+	this->parameters->FindParam(&num_responses,NumResponsesEnum);
+	this->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum);
+
+	for(resp=0;resp<num_responses;resp++) switch(responses[resp]){
+
+		case ThicknessAbsMisfitEnum:
+		case SurfaceAbsVelMisfitEnum:
+		case SurfaceRelVelMisfitEnum:
+		case SurfaceLogVelMisfitEnum:
+		case SurfaceLogVxVyMisfitEnum:
+		case SurfaceAverageVelMisfitEnum:
+			/*Nothing, J does not depends on the parameter being inverted for*/
+			break;
+		case DragCoefficientAbsGradientEnum:
+			GradjDragGradient(gradient,resp);
+			break;
+		default:
+			_error_("response %s not supported yet",EnumToStringx(responses[resp]));
+	}
 }
 /*}}}*/
@@ -3157,4 +3191,48 @@
 
 	delete friction;
+	delete gauss;
+}
+/*}}}*/
+/*FUNCTION Tria::GradjDragGradient{{{1*/
+void  Tria::GradjDragGradient(Vec gradient, int weight_index){
+
+	int        i,ig;
+	int        doflist1[NUMVERTICES];
+	double     Jdet,weight;
+	double     xyz_list[NUMVERTICES][3];
+	double     dh1dh3[NDOF2][NUMVERTICES];
+	double     dk[NDOF2]; 
+	double     grade_g[NUMVERTICES]={0.0};
+	GaussTria  *gauss=NULL;
+
+	/*Retrieve all inputs we will be needing: */
+	if(IsOnShelf())return;
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
+	GetDofList1(&doflist1[0]);
+	Input* dragcoefficient_input=inputs->GetInput(DragCoefficientEnum); _assert_(dragcoefficient_input);
+	Input* weights_input=inputs->GetInput(WeightsEnum);                 _assert_(weights_input);
+
+	/* 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);
+		GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss);
+		weights_input->GetParameterValue(&weight,gauss,weight_index);
+
+		/*Build alpha_complement_list: */
+		dragcoefficient_input->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],gauss);
+
+		/*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
+		for (i=0;i<NUMVERTICES;i++){
+			grade_g[i]+=-weight*Jdet*gauss->weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]);
+			_assert_(!isnan(grade_g[i]));
+		}
+	}
+	VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES);
+
+	/*Clean up and return*/
 	delete gauss;
 }
@@ -3641,6 +3719,6 @@
 		/*Create inputs and add to DataSetInput*/
 		DatasetInput* datasetinput=new DatasetInput(WeightsEnum);
-		for(j=0;j<iomodel->num_cm_responses;j++){
-			for(i=0;i<3;i++)nodeinputs[i]=iomodel->weights[tria_vertex_ids[i]-1];
+		for(i=0;i<iomodel->num_cm_responses;i++){
+			for(j=0;j<3;j++)nodeinputs[j]=iomodel->weights[(tria_vertex_ids[j]-1)*iomodel->num_cm_responses+i];
 			datasetinput->inputs->AddObject(new TriaVertexInput(WeightsEnum,nodeinputs));
 		}
Index: /issm/trunk/src/c/objects/Elements/Tria.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.h	(revision 8645)
+++ /issm/trunk/src/c/objects/Elements/Tria.h	(revision 8646)
@@ -93,4 +93,5 @@
 		void   GradjDragMacAyeal(Vec gradient);
 		void   GradjDragStokes(Vec gradient);
+		void   GradjDragGradient(Vec gradient,int weight_index);
 		void   GradjDhDtBalancedthickness(Vec gradient);
 		void   GradjVxBalancedthickness(Vec gradient);
Index: /issm/trunk/src/c/objects/Params/DoubleMatParam.cpp
===================================================================
--- /issm/trunk/src/c/objects/Params/DoubleMatParam.cpp	(revision 8645)
+++ /issm/trunk/src/c/objects/Params/DoubleMatParam.cpp	(revision 8646)
@@ -152,5 +152,5 @@
 
 /*DoubleMatParam virtual functions definitions: */
-/*FUNCTION DoubleMatParam::GetParameterValue{{{1*/
+/*FUNCTION DoubleMatParam::GetParameterValue(double** pdoublearray,int* pM,int* pN){{{1*/
 void  DoubleMatParam::GetParameterValue(double** pdoublearray,int* pM,int* pN){
 	double* output=NULL;
@@ -163,4 +163,22 @@
 	if(pN) *pN=N;
 	*pdoublearray=output;
+}
+/*}}}*/
+/*FUNCTION DoubleMatParam::GetParameterValue(int** pintarray,int* pM,int* pN){{{1*/
+void  DoubleMatParam::GetParameterValue(int** pintarray,int* pM,int* pN){
+#ifdef _SERIAL_
+	int* output=NULL;
+	int  i;
+
+	output=(int*)xmalloc((int)(M*N*sizeof(int)));
+	for(i=0;i<M*N;i++) output[i]=(int)value[i];
+
+	/*Assign output pointers:*/
+	if(pM) *pM=M;
+	if(pN) *pN=N;
+	*pintarray=output;
+#else
+	_error_("DoubleMat of enum %i (%s) cannot return an array of int",enum_type,EnumToStringx(enum_type));
+#endif
 }
 /*}}}*/
Index: /issm/trunk/src/c/objects/Params/DoubleMatParam.h
===================================================================
--- /issm/trunk/src/c/objects/Params/DoubleMatParam.h	(revision 8645)
+++ /issm/trunk/src/c/objects/Params/DoubleMatParam.h	(revision 8646)
@@ -54,5 +54,5 @@
 		void  GetParameterValue(int* pinteger){_error_("DoubleMat param of enum %i (%s) cannot return an integer",enum_type,EnumToStringx(enum_type));}
 		void  GetParameterValue(int** pintarray,int* pM){_error_("DoubleMat param of enum %i (%s) cannot return an array of integers",enum_type,EnumToStringx(enum_type));}
-		void  GetParameterValue(int** pintarray,int* pM,int* pN){_error_("DoubleMat param of enum %i (%s) cannot return an array of integers",enum_type,EnumToStringx(enum_type));}
+		void  GetParameterValue(int** pintarray,int* pM,int* pN);
 		void  GetParameterValue(double* pdouble){_error_("DoubleMat param of enum %i (%s) cannot return a double",enum_type,EnumToStringx(enum_type));}
 		void  GetParameterValue(char** pstring){_error_("DoubleMat param of enum %i (%s) cannot return a string",enum_type,EnumToStringx(enum_type));}
Index: /issm/trunk/src/c/objects/Params/DoubleVecParam.cpp
===================================================================
--- /issm/trunk/src/c/objects/Params/DoubleVecParam.cpp	(revision 8645)
+++ /issm/trunk/src/c/objects/Params/DoubleVecParam.cpp	(revision 8646)
@@ -174,5 +174,5 @@
 	*pintarray=output;
 #else
-	_error_("Double param of enum %i (%s) cannot return an array of double",enum_type,EnumToStringx(enum_type));
+	_error_("DoubleVec param of enum %i (%s) cannot return an array of int",enum_type,EnumToStringx(enum_type));
 #endif
 }