Index: sm/trunk/src/c/parallel/GradJCompute.cpp
===================================================================
--- /issm/trunk/src/c/parallel/GradJCompute.cpp	(revision 1999)
+++ 	(revision )
@@ -1,101 +1,0 @@
-/*!\file:  GradJCompute.cpp
- * \brief compute inverse method gradient direction.
- */ 
-
-#include "../issm.h"
-#include "./parallel.h"
-
-#undef __FUNCT__ 
-#define __FUNCT__ "GradJCompute"
-
-#ifdef HAVE_CONFIG_H
-	#include "config.h"
-#else
-#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
-#endif
-
-Vec GradJCompute(ParameterInputs* inputs,FemModel* femmodel){
-	
-	
-	/*intermediary: */
-	int analysis_type;
-	int sub_analysis_type;
-	int numberofnodes;
-	int numberofdofspernode;
-	char* solverstring=NULL;
-	char* control_type=NULL;
-	
-	Vec u_g=NULL;
-	double* u_g_double=NULL;
-
-	Vec du_g=NULL;
-	Vec du_f=NULL;
-
-	Vec lambda_f=NULL;
-	Vec lambda_g=NULL;
-	double* lambda_g_double=NULL;
-
-	Mat K_ff0=NULL;
-	Mat K_fs0=NULL;
-	
-	/*output: */
-	Vec grad_g=NULL;
-
-	/*flags: */
-	int debug=0;
-	int dim=-1;
-	int extrude_param=0;
-	
-	/*some parameters:*/
-	femmodel->parameters->FindParam((void*)&analysis_type,"analysis_type");
-	femmodel->parameters->FindParam((void*)&sub_analysis_type,"sub_analysis_type");
-	femmodel->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
-	femmodel->parameters->FindParam((void*)&numberofdofspernode,"numberofdofspernode");
-	femmodel->parameters->FindParam((void*)&solverstring,"solverstring");
-	femmodel->parameters->FindParam((void*)&control_type,"control_type");
-	femmodel->parameters->FindParam((void*)&extrude_param,"extrude_param");
-	femmodel->parameters->FindParam((void*)&debug,"debug");
-	femmodel->parameters->FindParam((void*)&dim,"dim");
-
-	_printf_("%s\n","      recover solution for this stiffness and right hand side:");
-	diagnostic_core_nonlinear(&u_g,&K_ff0,&K_fs0,NULL, femmodel,inputs,DiagnosticAnalysisEnum(),sub_analysis_type);
-	VecToMPISerial(&u_g_double,u_g); VecFree(&u_g);
-	inputs->Add("velocity",u_g_double,numberofdofspernode,numberofnodes);
-
-	_printf_("%s\n","      buid Du, difference between observed velocity and model velocity:");
-	Dux( &du_g, femmodel->elements,femmodel->nodes,femmodel->loads,femmodel->materials,inputs,analysis_type,sub_analysis_type);
-
-	_printf_("%s\n","      reduce adjoint load from g-set to f-set:");
-	Reduceloadfromgtofx(&du_f, du_g, femmodel->Gmn, K_fs0, femmodel->ys0, femmodel->nodesets);
-	VecFree(&du_g);MatFree(&K_fs0);
-
-	_printf_("%s\n","      solve for adjoint vector:");
-	Solverx(&lambda_f, K_ff0, du_f, NULL, solverstring);
-	VecFree(&du_f); MatFree(&K_ff0);
-	
-	_printf_("%s\n","      merge back to g set:");
-	Mergesolutionfromftogx(&lambda_g, lambda_f,femmodel->Gmn,femmodel->ys0,femmodel->nodesets);
-	VecFree(&lambda_f);
-
-	_printf_("%s\n","      compute gradJ:");
-	VecToMPISerial(&lambda_g_double,lambda_g);VecFree(&lambda_g);
-	inputs->Add("adjoint",lambda_g_double,numberofdofspernode,numberofnodes);
-
-	Gradjx( &grad_g, numberofnodes,femmodel->elements,femmodel->nodes, femmodel->loads, femmodel->materials, 
-				inputs,analysis_type,sub_analysis_type,control_type);
-
-	if (dim==3 && extrude_param){
-
-		_printf_("%s\n","      extruding gradient...");
-		FieldExtrudex( grad_g, femmodel->elements,femmodel->nodes,femmodel->loads,femmodel->materials,"gradj",0);
-	}
-	
-	/*Free ressources:*/
-	xfree((void**)&solverstring);
-	xfree((void**)&control_type);
-	xfree((void**)&u_g_double);
-	xfree((void**)&lambda_g_double);
-
-	return grad_g;
-
-}
Index: /issm/trunk/src/c/parallel/control_core.cpp
===================================================================
--- /issm/trunk/src/c/parallel/control_core.cpp	(revision 1999)
+++ /issm/trunk/src/c/parallel/control_core.cpp	(revision 2000)
@@ -21,4 +21,5 @@
 	/*Intermediary: */
 	DataSet* diagnostic_results=NULL;
+	DataSet* gradjcompute_results=NULL;
 	Vec     u_g=NULL;
 	double  search_scalar;
@@ -88,5 +89,8 @@
 
 		_printf_("%s\n","      computing gradJ...");
-		grad_g=GradJCompute(inputs,fem_model);
+		gradjcompute_results=new DataSet(ResultsEnum()); 
+		gradjcompute_core(gradjcompute_results,model, inputs);
+		gradjcompute_results->FindResult(&grad_g,"grad_g");
+		delete gradjcompute_results;
 		_printf_("%s\n","      done.");
 
Index: /issm/trunk/src/c/parallel/gradjcompute_core.cpp
===================================================================
--- /issm/trunk/src/c/parallel/gradjcompute_core.cpp	(revision 2000)
+++ /issm/trunk/src/c/parallel/gradjcompute_core.cpp	(revision 2000)
@@ -0,0 +1,118 @@
+/*!\file:  gradjcompute_core.cpp
+ * \brief compute inverse method gradient direction.
+ */ 
+
+#include "../issm.h"
+#include "./parallel.h"
+
+#undef __FUNCT__ 
+#define __FUNCT__ "gradjcompute_core"
+
+#ifdef HAVE_CONFIG_H
+	#include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+void gradjcompute_core(DataSet* results,Model* model, ParameterInputs* inputs){
+	
+	
+	/*intermediary: */
+	FemModel* femmodel=NULL;
+	DataSet* diagnostic_results=NULL;
+	int analysis_type;
+	int sub_analysis_type;
+	int numberofnodes;
+	int numberofdofspernode;
+	char* solverstring=NULL;
+	char* control_type=NULL;
+	
+	Vec u_g=NULL;
+	double* u_g_double=NULL;
+
+	Vec du_g=NULL;
+	Vec du_f=NULL;
+	Vec grad_g=NULL;
+
+	Vec lambda_f=NULL;
+	Vec lambda_g=NULL;
+	double* lambda_g_double=NULL;
+
+	Mat K_ff0=NULL;
+	Mat K_fs0=NULL;
+	
+	/*output: */
+	Result* result=NULL;
+
+	/*flags: */
+	int debug=0;
+	int dim=-1;
+	int extrude_param=0;
+	
+	/*some parameters:*/
+	femmodel=model->GetActiveFormulation();
+	femmodel->parameters->FindParam((void*)&analysis_type,"analysis_type");
+	femmodel->parameters->FindParam((void*)&sub_analysis_type,"sub_analysis_type");
+	femmodel->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
+	femmodel->parameters->FindParam((void*)&numberofdofspernode,"numberofdofspernode");
+	femmodel->parameters->FindParam((void*)&solverstring,"solverstring");
+	femmodel->parameters->FindParam((void*)&control_type,"control_type");
+	femmodel->parameters->FindParam((void*)&extrude_param,"extrude_param");
+	femmodel->parameters->FindParam((void*)&debug,"debug");
+	femmodel->parameters->FindParam((void*)&dim,"dim");
+
+	_printf_("%s\n","      recover solution for this stiffness and right hand side:");
+	diagnostic_core_nonlinear(&u_g,&K_ff0,&K_fs0,NULL, femmodel,inputs,DiagnosticAnalysisEnum(),sub_analysis_type);
+	VecToMPISerial(&u_g_double,u_g); VecFree(&u_g);
+	inputs->Add("velocity",u_g_double,numberofdofspernode,numberofnodes);
+
+	_printf_("%s\n","      buid Du, difference between observed velocity and model velocity:");
+	Dux( &du_g, femmodel->elements,femmodel->nodes,femmodel->loads,femmodel->materials,inputs,analysis_type,sub_analysis_type);
+
+	_printf_("%s\n","      reduce adjoint load from g-set to f-set:");
+	Reduceloadfromgtofx(&du_f, du_g, femmodel->Gmn, K_fs0, femmodel->ys0, femmodel->nodesets);
+	VecFree(&du_g);MatFree(&K_fs0);
+
+	_printf_("%s\n","      solve for adjoint vector:");
+	Solverx(&lambda_f, K_ff0, du_f, NULL, solverstring);
+	VecFree(&du_f); MatFree(&K_ff0);
+	
+	_printf_("%s\n","      merge back to g set:");
+	Mergesolutionfromftogx(&lambda_g, lambda_f,femmodel->Gmn,femmodel->ys0,femmodel->nodesets);
+	VecFree(&lambda_f);
+
+	_printf_("%s\n","      compute gradJ:");
+	VecToMPISerial(&lambda_g_double,lambda_g);VecFree(&lambda_g);
+	inputs->Add("adjoint",lambda_g_double,numberofdofspernode,numberofnodes);
+
+	Gradjx( &grad_g, numberofnodes,femmodel->elements,femmodel->nodes, femmodel->loads, femmodel->materials, 
+				inputs,analysis_type,sub_analysis_type,control_type);
+
+	if (dim==3 && extrude_param){
+
+		_printf_("%s\n","      extruding gradient...");
+		FieldExtrudex( grad_g, femmodel->elements,femmodel->nodes,femmodel->loads,femmodel->materials,"gradj",0);
+	}
+
+	if(analysis_type==SteadystateAnalysisEnum()){
+		diagnostic_results=new DataSet(ResultsEnum()); 
+		diagnostic_core(diagnostic_results,model, inputs);
+
+		//extract u_g and add it to input (3d velocity needed by thermal_core)
+		diagnostic_results->FindResult(&u_g,"u_g");
+		inputs->Add("velocity",u_g,3,numberofnodes);
+		delete diagnostic_results;
+	}
+
+	/*Plug results into output dataset: */
+	result=new Result(results->Size()+1,0,1,"grad_g",grad_g);
+	results->AddObject(result);
+	
+	/*Free ressources:*/
+	xfree((void**)&solverstring);
+	xfree((void**)&control_type);
+	xfree((void**)&u_g_double);
+	xfree((void**)&grad_g);
+	xfree((void**)&lambda_g_double);
+
+}
Index: /issm/trunk/src/c/parallel/parallel.h
===================================================================
--- /issm/trunk/src/c/parallel/parallel.h	(revision 1999)
+++ /issm/trunk/src/c/parallel/parallel.h	(revision 2000)
@@ -13,5 +13,5 @@
 class FemModel;
 
-Vec GradJCompute(ParameterInputs* inputs,FemModel* femmodel);
+void gradjcompute_core(DataSet* results,Model* model, ParameterInputs* inputs);
 
 void diagnostic_core(DataSet* results,Model* model, ParameterInputs* inputs);