Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/solutions/controltao_core.cpp
===================================================================
--- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/solutions/controltao_core.cpp	(revision 12909)
+++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/solutions/controltao_core.cpp	(revision 12910)
@@ -20,6 +20,7 @@
 int IssmMonitor(TaoSolver,void*);
 typedef struct {
 	FemModel* femmodel;
+	double*   J;
 } AppCtx;
 
 void controltao_core(FemModel* femmodel){
@@ -68,11 +69,12 @@
 	GetVectorFromControlInputsx(&X, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value");
 	GetVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
 	GetVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
-	TaoSetInitialVector(tao,X->vector);
-	TaoSetVariableBounds(tao,XL->vector,XU->vector);
+	TaoSetInitialVector(tao,X->pvector->vector);
+	TaoSetVariableBounds(tao,XL->pvector->vector,XU->pvector->vector);
 	xdelete(&XL);
 	xdelete(&XU);
 
+	user.J=(double*)xcalloc((maxiter+5),sizeof(double));
 	user.femmodel=femmodel;
 	TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void*)&user); 
 
@@ -80,11 +82,12 @@
 	if(VerboseControl()) _pprintLine_("   Starting optimization");
 	TaoSolve(tao);
 	TaoView(tao,PETSC_VIEWER_STDOUT_WORLD);
-	TaoGetSolutionVector(tao,&X->vector);
+	TaoGetSolutionVector(tao,&X->pvector->vector);
 	SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X);
 	for(int i=0;i<num_controls;i++){
 		InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_list[i]);
 	}
+	femmodel->results->AddObject(new DoubleVecExternalResult(femmodel->results->Size()+1,JEnum,user.J,maxiter+3,1,0));
 
 	/*Finalize*/
 	if(VerboseControl()) _pprintLine_("   preparing final solution");
@@ -141,7 +144,7 @@
 
 	/*Compute gradient*/
 	Gradjx(&gradient,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
-	VecCopy(gradient->vector,G); xdelete(&gradient);
+	VecCopy(gradient->pvector->vector,G); xdelete(&gradient);
 	VecScale(G,-1.);
 
 	/*Clean-up and return*/
@@ -165,6 +168,7 @@
 	if(its==0) _pprintLine_("Iter       Function      Residual  |  List of contributions");
 	if(its==0) _pprintLine_("-----------------------------------+-----------------------");
 	_pprintString_(setw(4)<<its<<"   "<<setw(12)<<setprecision(7)<<f<<"  "<<setw(12)<<setprecision(7)<<gnorm<<"  | ");
+	user->J[its]=f;
 
 	/*Retrieve objective functions independently*/
 	for(i=0;i<num_responses;i++){
Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/classes/matrix/Vector.cpp
===================================================================
--- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/classes/matrix/Vector.cpp	(revision 12909)
+++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/classes/matrix/Vector.cpp	(revision 12910)
@@ -60,6 +60,17 @@
 
 }
 /*}}}*/
+#ifdef _HAVE_PETSC_
+/*FUNCTION Vector::Vector(Vec petsc_vector){{{*/
+Vector::Vector(Vec petsc_vector){
+
+	this->type=PetscVecType;
+	this->svector=NULL;
+	this->pvector=new PetscVec(petsc_vector);
+
+}
+/*}}}*/
+#endif
 /*FUNCTION Vector::Vector(IssmDouble* serial_vec,int M,int type){{{*/
 Vector::Vector(IssmDouble* serial_vec,int M,int in_type){
 
Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/classes/matrix/Vector.h
===================================================================
--- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/classes/matrix/Vector.h	(revision 12909)
+++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/classes/matrix/Vector.h	(revision 12910)
@@ -35,6 +35,9 @@
 		Vector();
 		Vector(int M,bool fromlocalsize=false,int type=PetscVecType);
 		Vector(IssmDouble* serial_vec,int pM,int type=PetscVecType);
+		#ifdef _HAVE_PETSC_
+		Vector(Vec petsc_vector);
+		#endif
 		~Vector();
 		/*}}}*/
 		/*Vector specific routines {{{*/