Index: /issm/trunk/src/c/Makefile.am
===================================================================
--- /issm/trunk/src/c/Makefile.am	(revision 4838)
+++ /issm/trunk/src/c/Makefile.am	(revision 4839)
@@ -409,6 +409,4 @@
 					./modules/OutputResultsx/ElementResultsToPatch.cpp\
 					./modules/OutputResultsx/MatlabWriteResults.cpp\
-					./modules/Dux/Dux.h\
-					./modules/Dux/Dux.cpp\
 					./modules/MinVelx/MinVelx.h\
 					./modules/MinVelx/MinVelx.cpp\
@@ -970,6 +968,4 @@
 					./modules/OutputResultsx/ElementResultsToPatch.cpp\
 					./modules/OutputResultsx/FileWriteResults.cpp\
-					./modules/Dux/Dux.h\
-					./modules/Dux/Dux.cpp\
 					./modules/MinVelx/MinVelx.h\
 					./modules/MinVelx/MinVelx.cpp\
@@ -1149,4 +1145,5 @@
 					./solutions/SolutionConfiguration.cpp\
 					./solvers/solver_linear.cpp\
+					./solvers/solver_adjoint_linear.cpp\
 					./solvers/solver_diagnostic_nonlinear.cpp\
 					./solvers/solver_thermal_nonlinear.cpp\
Index: /issm/trunk/src/c/modules/modules.h
===================================================================
--- /issm/trunk/src/c/modules/modules.h	(revision 4838)
+++ /issm/trunk/src/c/modules/modules.h	(revision 4839)
@@ -21,5 +21,4 @@
 #include "./CostFunctionx/CostFunctionx.h"
 #include "./DakotaResponsesx/DakotaResponsesx.h"
-#include "./Dux/Dux.h"
 #include "./ElementConnectivityx/ElementConnectivityx.h"
 #include "./FieldAverageOntoVerticesx/FieldAverageOntoVerticesx.h"
Index: /issm/trunk/src/c/objects/Elements/Beam.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Beam.cpp	(revision 4838)
+++ /issm/trunk/src/c/objects/Elements/Beam.cpp	(revision 4839)
@@ -285,9 +285,4 @@
 	}
 
-}
-/*}}}*/
-/*FUNCTION Beam::Du{{{1*/
-void  Beam::Du(Vec){
-	ISSMERROR(" not supported yet!");
 }
 /*}}}*/
Index: /issm/trunk/src/c/objects/Elements/Element.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Element.h	(revision 4838)
+++ /issm/trunk/src/c/objects/Elements/Element.h	(revision 4839)
@@ -37,5 +37,4 @@
 		virtual void   GetThicknessList(double* thickness_list)=0;
 		virtual void   GetBedList(double* bed_list)=0;
-		virtual void   Du(Vec du_g)=0;
 		virtual void   Gradj(Vec gradient,int control_type)=0;
 		virtual void   GradjDrag(Vec gradient)=0;
Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 4838)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 4839)
@@ -356,8 +356,5 @@
 
 	/*Just branch to the correct InputUpdateFromSolution generator, according to the type of analysis we are carrying out: */
-	if (analysis_type==ControlAnalysisEnum){
-		InputUpdateFromSolutionDiagnosticHoriz( solution);
-	}
-	else if (analysis_type==DiagnosticHorizAnalysisEnum){
+	if (analysis_type==DiagnosticHorizAnalysisEnum){
 		InputUpdateFromSolutionDiagnosticHoriz( solution);
 	}
@@ -730,8 +727,9 @@
 
 	/*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
-	if (analysis_type==ControlAnalysisEnum){
+	if (analysis_type==DiagnosticHorizAnalysisEnum){
 		CreateKMatrixDiagnosticHoriz( Kgg);
 	}
-	else if (analysis_type==DiagnosticHorizAnalysisEnum){
+	else if (analysis_type==AdjointHorizAnalysisEnum){
+		/*Same as diagnostic*/
 		CreateKMatrixDiagnosticHoriz( Kgg);
 	}
@@ -745,4 +743,8 @@
 		CreateKMatrixDiagnosticStokes( Kgg);
 	}
+	else if (analysis_type==AdjointStokesAnalysisEnum){
+		/*Same as diagnostic*/
+		CreateKMatrixDiagnosticStokes( Kgg);
+	}
 	else if (analysis_type==BedSlopeXAnalysisEnum || analysis_type==SurfaceSlopeXAnalysisEnum || analysis_type==BedSlopeYAnalysisEnum || analysis_type==SurfaceSlopeYAnalysisEnum){
 		CreateKMatrixSlope( Kgg);
@@ -779,9 +781,9 @@
 
 	/*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
-	if (analysis_type==ControlAnalysisEnum){
+	if (analysis_type==DiagnosticHorizAnalysisEnum){
 		CreatePVectorDiagnosticHoriz( pg);
 	}
-	else if (analysis_type==DiagnosticHorizAnalysisEnum){
-		CreatePVectorDiagnosticHoriz( pg);
+	else if (analysis_type==AdjointHorizAnalysisEnum){
+		CreatePVectorAdjointHoriz( pg);
 	}
 	else if (analysis_type==DiagnosticHutterAnalysisEnum){
@@ -794,4 +796,7 @@
 		CreatePVectorDiagnosticStokes( pg);
 	}
+	else if (analysis_type==AdjointStokesAnalysisEnum){
+		CreatePVectorAdjointStokes( pg);
+	}
 	else if (analysis_type==BedSlopeXAnalysisEnum || analysis_type==SurfaceSlopeXAnalysisEnum || analysis_type==BedSlopeYAnalysisEnum || analysis_type==SurfaceSlopeYAnalysisEnum){
 		CreatePVectorSlope( pg);
@@ -814,49 +819,4 @@
 	else ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumAsString(analysis_type));
 
-}
-/*}}}*/
-/*FUNCTION Penta::Du {{{1*/
-void  Penta::Du(Vec du_g){
-
-	int i;
-	Tria* tria=NULL;
-
-	/*inputs: */
-	bool onwater;
-	bool collapse;
-	bool onsurface;
-	bool onbed;
-
-	/*retrieve inputs :*/
-	inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
-	inputs->GetParameterValue(&collapse,CollapseEnum);
-	inputs->GetParameterValue(&onbed,ElementOnBedEnum);
-	inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum);
-
-	/*If on water, skip: */
-	if(onwater)return;
-
-	/*Bail out if this element if:
-	 * -> Non collapsed and not on the surface
-	 * -> collapsed (2d model) and not on bed) */
-	if ((!collapse && !onsurface) || (collapse && !onbed)){
-		return;
-	}
-	else if (collapse){
-
-		/*This element should be collapsed into a tria element at its base. Create this tria element, 
-		 * and compute Du*/
-		tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face).
-		tria->Du(du_g);
-		delete tria;
-		return;
-	}
-	else{
-
-		tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
-		tria->Du(du_g);
-		delete tria;
-		return;
-	}
 }
 /*}}}*/
@@ -4216,4 +4176,49 @@
 }
 /*}}}*/
+/*FUNCTION Penta::CreatePVectorAdjointHoriz{{{1*/
+void  Penta::CreatePVectorAdjointHoriz(Vec p_g){
+
+	int i;
+	Tria* tria=NULL;
+
+	/*inputs: */
+	bool onwater;
+	bool collapse;
+	bool onsurface;
+	bool onbed;
+
+	/*retrieve inputs :*/
+	inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
+	inputs->GetParameterValue(&collapse,CollapseEnum);
+	inputs->GetParameterValue(&onbed,ElementOnBedEnum);
+	inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum);
+
+	/*If on water, skip: */
+	if(onwater) return;
+
+	/*Bail out if this element if:
+	 * -> Non collapsed and not on the surface
+	 * -> collapsed (2d model) and not on bed) */
+	if ((!collapse && !onsurface) || (collapse && !onbed)){
+		return;
+	}
+	else if (collapse){
+
+		/*This element should be collapsed into a tria element at its base. Create this tria element, 
+		 * and compute pe_g*/
+		tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face).
+		tria->CreatePVectorAdjointHoriz(p_g);
+		delete tria;
+		return;
+	}
+	else{
+
+		tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
+		tria->CreatePVectorAdjointHoriz(p_g);
+		delete tria;
+		return;
+	}
+}
+/*}}}*/
 /*FUNCTION Penta::CreatePVectorDiagnosticHutter{{{1*/
 void  Penta::CreatePVectorDiagnosticHutter(Vec pg){
@@ -4497,4 +4502,28 @@
 	xfree((void**)&vert_gauss_weights);
 
+}
+/*}}}*/
+/*FUNCTION Penta::CreatePVectorAdjointStokes{{{1*/
+void  Penta::CreatePVectorAdjointStokes(Vec p_g){
+
+	int i;
+	Tria* tria=NULL;
+
+	/*inputs: */
+	bool onwater;
+	bool onsurface;
+
+	/*retrieve inputs :*/
+	inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
+	inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum);
+
+	/*If on water, skip: */
+	if(onwater || !onsurface)return;
+
+	/*Call Tria's function*/
+	tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
+	tria->CreatePVectorAdjointStokes(p_g);
+	delete tria;
+	return;
 }
 /*}}}*/
Index: /issm/trunk/src/c/objects/Elements/Penta.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.h	(revision 4838)
+++ /issm/trunk/src/c/objects/Elements/Penta.h	(revision 4839)
@@ -73,5 +73,4 @@
 		void   CreateKMatrix(Mat Kgg);
 		void   CreatePVector(Vec pg);
-		void   Du(Vec du_g);
 		void   GetBedList(double* bed_list);
 		void*  GetMatPar();
@@ -127,6 +126,8 @@
 		void	  CreatePVectorBalancedvelocities( Vec pg);
 		void	  CreatePVectorDiagnosticHoriz( Vec pg);
+		void	  CreatePVectorAdjointHoriz( Vec pg);
 		void	  CreatePVectorDiagnosticHutter( Vec pg);
 		void	  CreatePVectorDiagnosticStokes( Vec pg);
+		void	  CreatePVectorAdjointStokes( Vec pg);
 		void	  CreatePVectorDiagnosticVert( Vec pg);
 		void	  CreatePVectorMelting( Vec pg);
Index: /issm/trunk/src/c/objects/Elements/Sing.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Sing.cpp	(revision 4838)
+++ /issm/trunk/src/c/objects/Elements/Sing.cpp	(revision 4839)
@@ -256,9 +256,4 @@
 	}
 
-}
-/*}}}*/
-/*FUNCTION Sing::Du {{{1*/
-void  Sing::Du(Vec){
-	ISSMERROR(" not supported yet!");
 }
 /*}}}*/
Index: /issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 4838)
+++ /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 4839)
@@ -364,8 +364,5 @@
 
 	/*Just branch to the correct InputUpdateFromSolution generator, according to the type of analysis we are carrying out: */
-	if (analysis_type==ControlAnalysisEnum){
-		InputUpdateFromSolutionDiagnosticHoriz( solution);
-	}
-	else if (analysis_type==DiagnosticHorizAnalysisEnum){
+	if (analysis_type==DiagnosticHorizAnalysisEnum){
 		InputUpdateFromSolutionDiagnosticHoriz( solution);
 	}
@@ -664,10 +661,11 @@
 
 	/*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
-	if (analysis_type==ControlAnalysisEnum){
+	if (analysis_type==DiagnosticHorizAnalysisEnum){
 		CreateKMatrixDiagnosticHoriz( Kgg);
 	}
-	else if (analysis_type==DiagnosticHorizAnalysisEnum){
+	else if (analysis_type==AdjointHorizAnalysisEnum){
+		/*Same as diagnostic*/
 		CreateKMatrixDiagnosticHoriz( Kgg);
-		}
+	}
 	else if (analysis_type==DiagnosticHutterAnalysisEnum){
 		CreateKMatrixDiagnosticHutter( Kgg);
@@ -715,9 +713,9 @@
 
 	/*Just branch to the correct load generator, according to the type of analysis we are carrying out: */
-	if (analysis_type==ControlAnalysisEnum){
+	if (analysis_type==DiagnosticHorizAnalysisEnum){
 		CreatePVectorDiagnosticHoriz( pg);
 	}
-	else if (analysis_type==DiagnosticHorizAnalysisEnum){
-		CreatePVectorDiagnosticHoriz( pg);
+	else if (analysis_type==AdjointHorizAnalysisEnum){
+		CreatePVectorAdjointHoriz( pg);
 	}
 	else if (analysis_type==DiagnosticHutterAnalysisEnum){
@@ -749,243 +747,4 @@
 		ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumAsString(analysis_type));
 	}
-
-}
-/*}}}*/
-/*FUNCTION Tria::Du {{{1*/
-void Tria::Du(Vec du_g){
-
-	int i;
-
-	/* node data: */
-	const int    numgrids=3;
-	const int    numdof=2*numgrids;
-	const int    NDOF2=2;
-	double       xyz_list[numgrids][3];
-	int          doflist[numdof];
-	int          numberofdofspernode;
-	double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
-
-	/* grid data: */
-	double vx_list[numgrids];
-	double vy_list[numgrids];
-	double obs_vx_list[numgrids];
-	double obs_vy_list[numgrids];
-	double dux_list[numgrids];
-	double duy_list[numgrids];
-	double weights_list[numgrids];
-
-	/* gaussian points: */
-	int     num_gauss,ig;
-	double* first_gauss_area_coord  =  NULL;
-	double* second_gauss_area_coord =  NULL;
-	double* third_gauss_area_coord  =  NULL;
-	double* gauss_weights           =  NULL;
-	double  gauss_weight;
-	double  gauss_l1l2l3[3];
-
-	/* parameters: */
-	double  obs_velocity_mag,velocity_mag;
-	double  dux,duy;
-	double  meanvel, epsvel;
-
-	/*element vector : */
-	double  due_g[numdof]={0.0};
-	double  due_g_gaussian[numdof];
-
-	/* Jacobian: */
-	double Jdet;
-
-	/*nodal functions: */
-	double l1l2l3[3];
-
-	/*relative and algorithmic fitting: */
-	double scalex=0;
-	double scaley=0;
-	double scale=0;
-	double S=0;
-	int    fit=-1;
-
-	/*retrieve some parameters: */
-	this->parameters->FindParam(&meanvel,MeanVelEnum);
-	this->parameters->FindParam(&epsvel,EpsVelEnum);
-
-	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
-	GetDofList(&doflist[0],&numberofdofspernode);
-
-	/* Recover input data: */
-	inputs->GetParameterValues(&obs_vx_list[0],&gaussgrids[0][0],3,VxObsEnum);
-	inputs->GetParameterValues(&obs_vy_list[0],&gaussgrids[0][0],3,VyObsEnum);
-
-	inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3,VxEnum);
-	inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyEnum);
-	
-	inputs->GetParameterValues(&weights_list[0],&gaussgrids[0][0],3,WeightsEnum);
-	
-	inputs->GetParameterValue(&fit,FitEnum);
-	if(fit==3){
-		inputs->GetParameterValue(&S,SurfaceAreaEnum);
-	}
-
-	/*Get Du at the 3 nodes (integration of the linearized function)
-	 * Here we integrate linearized functions:
-	 *               
-	 * J(E) = int_E   sum_{i=1}^3  J_i Phi_i
-	 *
-	 *       d J                  dJ_i
-	 * DU= - --- = sum_{i=1}^3  - ---  Phi_i = sum_{i=1}^3 DU_i Phi_i
-	 *       d u                  du_i
-	 *
-	 * where J_i are the misfits at the 3 nodes of the triangle
-	 *       Phi_i is the nodal function (P1) with respect to 
-	 *       the vertex i
-	 */
-	if(fit==0){
-		/*We are using an absolute misfit:
-		 *
-		 *      1  [           2              2 ]
-		 * J = --- | (u - u   )  +  (v - v   )  |
-		 *      2  [       obs            obs   ]
-		 *
-		 *        dJ             2
-		 * DU = - -- = (u   - u )
-		 *        du     obs
-		 */
-		for (i=0;i<numgrids;i++){
-			dux_list[i]=obs_vx_list[i]-vx_list[i];
-			duy_list[i]=obs_vy_list[i]-vy_list[i];
-		}
-	}
-	else if(fit==1){
-		/*We are using a relative misfit: 
-		 *                        
-		 *      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<numgrids;i++){
-			scalex=pow(meanvel/(obs_vx_list[i]+epsvel),2);
-			scaley=pow(meanvel/(obs_vy_list[i]+epsvel),2);
-			if(obs_vx_list[i]==0)scalex=0;
-			if(obs_vy_list[i]==0)scaley=0;
-			dux_list[i]=scalex*(obs_vx_list[i]-vx_list[i]);
-			duy_list[i]=scaley*(obs_vy_list[i]-vy_list[i]);
-		}
-	}
-	else if(fit==2){
-		/*We are using a logarithmic misfit:
-		 *                        
-		 *                 [        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<numgrids;i++){
-			velocity_mag=sqrt(pow(vx_list[i],2)+pow(vy_list[i],2))+epsvel; //epsvel to avoid velocity being nil.
-			obs_velocity_mag=sqrt(pow(obs_vx_list[i],2)+pow(obs_vy_list[i],2))+epsvel; //epsvel to avoid observed velocity being nil.
-			scale=-8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag);
-			dux_list[i]=scale*vx_list[i];
-			duy_list[i]=scale*vy_list[i];
-		}
-	}
-	else if(fit==3){
-		/*We are using an spacially average absolute misfit:
-		 *
-		 *      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<numgrids;i++){
-			scale=1.0/(S*sqrt(pow(vx_list[i]-obs_vx_list[i],2)+pow(vy_list[i]-obs_vx_list[i],2))+epsvel);
-			dux_list[i]=scale*(obs_vx_list[i]-vx_list[i]);
-			duy_list[i]=scale*(obs_vy_list[i]-vy_list[i]);
-		}
-	}
-	else if(fit==4){
-		/*We are using an logarithmic 2 misfit:
-		 *
-		 *      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<numgrids;i++){
-			dux_list[i] = - pow(meanvel,(double)2)*(
-						log((fabs(vx_list[i])+epsvel)/(fabs(obs_vx_list[i])+epsvel)) * 1/(vx_list[i]+epsvel));
-			duy_list[i] = - pow(meanvel,(double)2)*(
-						log((fabs(vy_list[i])+epsvel)/(fabs(obs_vy_list[i])+epsvel)) * 1/(vy_list[i]+epsvel));
-		}
-	}
-	else{
-		/*Not supported yet! : */
-		ISSMERROR("%s%g","unsupported type of fit: ",fit);
-	}
-
-	/*Apply weights to DU*/
-	for (i=0;i<numgrids;i++){
-		dux_list[i]=weights_list[i]*dux_list[i];
-		duy_list[i]=weights_list[i]*duy_list[i];
-	}
-
-	/* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
-	GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
-
-	/* Start  looping on the number of gaussian points: */
-	for (ig=0; ig<num_gauss; ig++){
-		/*Pick up the gaussian point: */
-		gauss_weight=*(gauss_weights+ig);
-		gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 
-		gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
-		gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
-
-		/* Get Jacobian determinant: */
-		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
-
-		/* Get nodal functions value at gaussian point:*/
-		GetNodalFunctions(l1l2l3, gauss_l1l2l3);
-
-		/*Build due_g_gaussian vector: we have three cases here, according to which type of misfit we are using. */
-
-		/*Compute absolute(x/y) at gaussian point: */
-		GetParameterValue(&dux, &dux_list[0],gauss_l1l2l3);
-		GetParameterValue(&duy, &duy_list[0],gauss_l1l2l3);
-
-		/*compute Du*/
-		for (i=0;i<numgrids;i++){
-			due_g_gaussian[i*NDOF2+0]=dux*Jdet*gauss_weight*l1l2l3[i]; 
-			due_g_gaussian[i*NDOF2+1]=duy*Jdet*gauss_weight*l1l2l3[i]; 
-		}
-
-		/*Add due_g_gaussian vector to due_g: */
-		for( i=0; i<numdof; i++){
-			due_g[i]+=due_g_gaussian[i];
-		}
-	}
-
-	/*Add due_g to global vector du_g: */
-	VecSetValues(du_g,numdof,doflist,(const double*)due_g,ADD_VALUES);
-
-	xfree((void**)&first_gauss_area_coord);
-	xfree((void**)&second_gauss_area_coord);
-	xfree((void**)&third_gauss_area_coord);
-	xfree((void**)&gauss_weights);
 
 }
@@ -4389,4 +4148,484 @@
 }
 /*}}}*/
+/*FUNCTION Tria::CreatePVectorAdjointHoriz{{{1*/
+void Tria::CreatePVectorAdjointHoriz(Vec p_g){
+
+	int i;
+
+	/* node data: */
+	const int    numgrids=3;
+	const int    numdof=2*numgrids;
+	const int    NDOF2=2;
+	double       xyz_list[numgrids][3];
+	int          doflist[numdof];
+	int          numberofdofspernode;
+	double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
+
+	/* grid data: */
+	double vx_list[numgrids];
+	double vy_list[numgrids];
+	double obs_vx_list[numgrids];
+	double obs_vy_list[numgrids];
+	double dux_list[numgrids];
+	double duy_list[numgrids];
+	double weights_list[numgrids];
+
+	/* gaussian points: */
+	int     num_gauss,ig;
+	double* first_gauss_area_coord  =  NULL;
+	double* second_gauss_area_coord =  NULL;
+	double* third_gauss_area_coord  =  NULL;
+	double* gauss_weights           =  NULL;
+	double  gauss_weight;
+	double  gauss_l1l2l3[3];
+
+	/* parameters: */
+	double  obs_velocity_mag,velocity_mag;
+	double  dux,duy;
+	double  meanvel, epsvel;
+
+	/*element vector : */
+	double  pe_g[numdof]={0.0};
+	double  pe_g_gaussian[numdof];
+
+	/* Jacobian: */
+	double Jdet;
+
+	/*nodal functions: */
+	double l1l2l3[3];
+
+	/*relative and algorithmic fitting: */
+	double scalex=0;
+	double scaley=0;
+	double scale=0;
+	double S=0;
+	int    fit=-1;
+
+	/*retrieve some parameters: */
+	this->parameters->FindParam(&meanvel,MeanVelEnum);
+	this->parameters->FindParam(&epsvel,EpsVelEnum);
+
+	/* Get node coordinates and dof list: */
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
+	GetDofList(&doflist[0],&numberofdofspernode);
+
+	/* Recover input data: */
+	inputs->GetParameterValues(&obs_vx_list[0],&gaussgrids[0][0],3,VxObsEnum);
+	inputs->GetParameterValues(&obs_vy_list[0],&gaussgrids[0][0],3,VyObsEnum);
+
+	inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3,VxEnum);
+	inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyEnum);
+
+	inputs->GetParameterValues(&weights_list[0],&gaussgrids[0][0],3,WeightsEnum);
+
+	inputs->GetParameterValue(&fit,FitEnum);
+	if(fit==3){
+		inputs->GetParameterValue(&S,SurfaceAreaEnum);
+	}
+
+	/*Get Du at the 3 nodes (integration of the linearized function)
+	 * Here we integrate linearized functions:
+	 *               
+	 * J(E) = int_E   sum_{i=1}^3  J_i Phi_i
+	 *
+	 *       d J                  dJ_i
+	 * DU= - --- = sum_{i=1}^3  - ---  Phi_i = sum_{i=1}^3 DU_i Phi_i
+	 *       d u                  du_i
+	 *
+	 * where J_i are the misfits at the 3 nodes of the triangle
+	 *       Phi_i is the nodal function (P1) with respect to 
+	 *       the vertex i
+	 */
+	if(fit==0){
+		/*We are using an absolute misfit:
+		 *
+		 *      1  [           2              2 ]
+		 * J = --- | (u - u   )  +  (v - v   )  |
+		 *      2  [       obs            obs   ]
+		 *
+		 *        dJ             2
+		 * DU = - -- = (u   - u )
+		 *        du     obs
+		 */
+		for (i=0;i<numgrids;i++){
+			dux_list[i]=obs_vx_list[i]-vx_list[i];
+			duy_list[i]=obs_vy_list[i]-vy_list[i];
+		}
+	}
+	else if(fit==1){
+		/*We are using a relative misfit: 
+		 *                        
+		 *      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<numgrids;i++){
+			scalex=pow(meanvel/(obs_vx_list[i]+epsvel),2);
+			scaley=pow(meanvel/(obs_vy_list[i]+epsvel),2);
+			if(obs_vx_list[i]==0)scalex=0;
+			if(obs_vy_list[i]==0)scaley=0;
+			dux_list[i]=scalex*(obs_vx_list[i]-vx_list[i]);
+			duy_list[i]=scaley*(obs_vy_list[i]-vy_list[i]);
+		}
+	}
+	else if(fit==2){
+		/*We are using a logarithmic misfit:
+		 *                        
+		 *                 [        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<numgrids;i++){
+			velocity_mag=sqrt(pow(vx_list[i],2)+pow(vy_list[i],2))+epsvel; //epsvel to avoid velocity being nil.
+			obs_velocity_mag=sqrt(pow(obs_vx_list[i],2)+pow(obs_vy_list[i],2))+epsvel; //epsvel to avoid observed velocity being nil.
+			scale=-8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag);
+			dux_list[i]=scale*vx_list[i];
+			duy_list[i]=scale*vy_list[i];
+		}
+	}
+	else if(fit==3){
+		/*We are using an spacially average absolute misfit:
+		 *
+		 *      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<numgrids;i++){
+			scale=1.0/(S*sqrt(pow(vx_list[i]-obs_vx_list[i],2)+pow(vy_list[i]-obs_vx_list[i],2))+epsvel);
+			dux_list[i]=scale*(obs_vx_list[i]-vx_list[i]);
+			duy_list[i]=scale*(obs_vy_list[i]-vy_list[i]);
+		}
+	}
+	else if(fit==4){
+		/*We are using an logarithmic 2 misfit:
+		 *
+		 *      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<numgrids;i++){
+			dux_list[i] = - pow(meanvel,(double)2)*(
+						log((fabs(vx_list[i])+epsvel)/(fabs(obs_vx_list[i])+epsvel)) * 1/(vx_list[i]+epsvel));
+			duy_list[i] = - pow(meanvel,(double)2)*(
+						log((fabs(vy_list[i])+epsvel)/(fabs(obs_vy_list[i])+epsvel)) * 1/(vy_list[i]+epsvel));
+		}
+	}
+	else{
+		/*Not supported yet! : */
+		ISSMERROR("%s%g","unsupported type of fit: ",fit);
+	}
+
+	/*Apply weights to DU*/
+	for (i=0;i<numgrids;i++){
+		dux_list[i]=weights_list[i]*dux_list[i];
+		duy_list[i]=weights_list[i]*duy_list[i];
+	}
+
+	/* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
+	GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
+
+	/* Start  looping on the number of gaussian points: */
+	for (ig=0; ig<num_gauss; ig++){
+		/*Pick up the gaussian point: */
+		gauss_weight=*(gauss_weights+ig);
+		gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 
+		gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
+		gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
+
+		/* Get Jacobian determinant: */
+		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
+
+		/* Get nodal functions value at gaussian point:*/
+		GetNodalFunctions(l1l2l3, gauss_l1l2l3);
+
+		/*Build due_g_gaussian vector: we have three cases here, according to which type of misfit we are using. */
+
+		/*Compute absolute(x/y) at gaussian point: */
+		GetParameterValue(&dux, &dux_list[0],gauss_l1l2l3);
+		GetParameterValue(&duy, &duy_list[0],gauss_l1l2l3);
+
+		/*compute Du*/
+		for (i=0;i<numgrids;i++){
+			pe_g_gaussian[i*NDOF2+0]=dux*Jdet*gauss_weight*l1l2l3[i]; 
+			pe_g_gaussian[i*NDOF2+1]=duy*Jdet*gauss_weight*l1l2l3[i]; 
+		}
+
+		/*Add pe_g_gaussian vector to pe_g: */
+		for( i=0; i<numdof; i++){
+			pe_g[i]+=pe_g_gaussian[i];
+		}
+	}
+
+	/*Add pe_g to global vector p_g: */
+	VecSetValues(p_g,numdof,doflist,(const double*)pe_g,ADD_VALUES);
+
+	/*Clean up*/
+	xfree((void**)&first_gauss_area_coord);
+	xfree((void**)&second_gauss_area_coord);
+	xfree((void**)&third_gauss_area_coord);
+	xfree((void**)&gauss_weights);
+
+}
+/*}}}*/
+/*FUNCTION Tria::CreatePVectorAdjointStokes{{{1*/
+void Tria::CreatePVectorAdjointStokes(Vec p_g){
+
+	int i;
+
+	/* node data: */
+	const int    numgrids=3;
+	const int    numdof=2*numgrids;
+	const int    NDOF2=2;
+	double       xyz_list[numgrids][3];
+	int          doflist[numdof];
+	int          numberofdofspernode;
+	double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
+
+	/* grid data: */
+	double vx_list[numgrids];
+	double vy_list[numgrids];
+	double obs_vx_list[numgrids];
+	double obs_vy_list[numgrids];
+	double dux_list[numgrids];
+	double duy_list[numgrids];
+	double weights_list[numgrids];
+
+	/* gaussian points: */
+	int     num_gauss,ig;
+	double* first_gauss_area_coord  =  NULL;
+	double* second_gauss_area_coord =  NULL;
+	double* third_gauss_area_coord  =  NULL;
+	double* gauss_weights           =  NULL;
+	double  gauss_weight;
+	double  gauss_l1l2l3[3];
+
+	/* parameters: */
+	double  obs_velocity_mag,velocity_mag;
+	double  dux,duy;
+	double  meanvel, epsvel;
+
+	/*element vector : */
+	double  pe_g[numdof]={0.0};
+	double  pe_g_gaussian[numdof];
+
+	/* Jacobian: */
+	double Jdet;
+
+	/*nodal functions: */
+	double l1l2l3[3];
+
+	/*relative and algorithmic fitting: */
+	double scalex=0;
+	double scaley=0;
+	double scale=0;
+	double S=0;
+	int    fit=-1;
+
+	/*retrieve some parameters: */
+	this->parameters->FindParam(&meanvel,MeanVelEnum);
+	this->parameters->FindParam(&epsvel,EpsVelEnum);
+
+	/* Get node coordinates and dof list: */
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
+	GetDofList(&doflist[0],&numberofdofspernode);
+
+	/* Recover input data: */
+	inputs->GetParameterValues(&obs_vx_list[0],&gaussgrids[0][0],3,VxObsEnum);
+	inputs->GetParameterValues(&obs_vy_list[0],&gaussgrids[0][0],3,VyObsEnum);
+
+	inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3,VxEnum);
+	inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyEnum);
+
+	inputs->GetParameterValues(&weights_list[0],&gaussgrids[0][0],3,WeightsEnum);
+
+	inputs->GetParameterValue(&fit,FitEnum);
+	if(fit==3){
+		inputs->GetParameterValue(&S,SurfaceAreaEnum);
+	}
+
+	/*Get Du at the 3 nodes (integration of the linearized function)
+	 * Here we integrate linearized functions:
+	 *               
+	 * J(E) = int_E   sum_{i=1}^3  J_i Phi_i
+	 *
+	 *       d J                  dJ_i
+	 * DU= - --- = sum_{i=1}^3  - ---  Phi_i = sum_{i=1}^3 DU_i Phi_i
+	 *       d u                  du_i
+	 *
+	 * where J_i are the misfits at the 3 nodes of the triangle
+	 *       Phi_i is the nodal function (P1) with respect to 
+	 *       the vertex i
+	 */
+	if(fit==0){
+		/*We are using an absolute misfit:
+		 *
+		 *      1  [           2              2 ]
+		 * J = --- | (u - u   )  +  (v - v   )  |
+		 *      2  [       obs            obs   ]
+		 *
+		 *        dJ             2
+		 * DU = - -- = (u   - u )
+		 *        du     obs
+		 */
+		for (i=0;i<numgrids;i++){
+			dux_list[i]=obs_vx_list[i]-vx_list[i];
+			duy_list[i]=obs_vy_list[i]-vy_list[i];
+		}
+	}
+	else if(fit==1){
+		/*We are using a relative misfit: 
+		 *                        
+		 *      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<numgrids;i++){
+			scalex=pow(meanvel/(obs_vx_list[i]+epsvel),2);
+			scaley=pow(meanvel/(obs_vy_list[i]+epsvel),2);
+			if(obs_vx_list[i]==0)scalex=0;
+			if(obs_vy_list[i]==0)scaley=0;
+			dux_list[i]=scalex*(obs_vx_list[i]-vx_list[i]);
+			duy_list[i]=scaley*(obs_vy_list[i]-vy_list[i]);
+		}
+	}
+	else if(fit==2){
+		/*We are using a logarithmic misfit:
+		 *                        
+		 *                 [        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<numgrids;i++){
+			velocity_mag=sqrt(pow(vx_list[i],2)+pow(vy_list[i],2))+epsvel; //epsvel to avoid velocity being nil.
+			obs_velocity_mag=sqrt(pow(obs_vx_list[i],2)+pow(obs_vy_list[i],2))+epsvel; //epsvel to avoid observed velocity being nil.
+			scale=-8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag);
+			dux_list[i]=scale*vx_list[i];
+			duy_list[i]=scale*vy_list[i];
+		}
+	}
+	else if(fit==3){
+		/*We are using an spacially average absolute misfit:
+		 *
+		 *      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<numgrids;i++){
+			scale=1.0/(S*sqrt(pow(vx_list[i]-obs_vx_list[i],2)+pow(vy_list[i]-obs_vx_list[i],2))+epsvel);
+			dux_list[i]=scale*(obs_vx_list[i]-vx_list[i]);
+			duy_list[i]=scale*(obs_vy_list[i]-vy_list[i]);
+		}
+	}
+	else if(fit==4){
+		/*We are using an logarithmic 2 misfit:
+		 *
+		 *      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<numgrids;i++){
+			dux_list[i] = - pow(meanvel,(double)2)*(
+						log((fabs(vx_list[i])+epsvel)/(fabs(obs_vx_list[i])+epsvel)) * 1/(vx_list[i]+epsvel));
+			duy_list[i] = - pow(meanvel,(double)2)*(
+						log((fabs(vy_list[i])+epsvel)/(fabs(obs_vy_list[i])+epsvel)) * 1/(vy_list[i]+epsvel));
+		}
+	}
+	else{
+		/*Not supported yet! : */
+		ISSMERROR("%s%g","unsupported type of fit: ",fit);
+	}
+
+	/*Apply weights to DU*/
+	for (i=0;i<numgrids;i++){
+		dux_list[i]=weights_list[i]*dux_list[i];
+		duy_list[i]=weights_list[i]*duy_list[i];
+	}
+
+	/* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
+	GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
+
+	/* Start  looping on the number of gaussian points: */
+	for (ig=0; ig<num_gauss; ig++){
+		/*Pick up the gaussian point: */
+		gauss_weight=*(gauss_weights+ig);
+		gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 
+		gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
+		gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
+
+		/* Get Jacobian determinant: */
+		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
+
+		/* Get nodal functions value at gaussian point:*/
+		GetNodalFunctions(l1l2l3, gauss_l1l2l3);
+
+		/*Build due_g_gaussian vector: we have three cases here, according to which type of misfit we are using. */
+
+		/*Compute absolute(x/y) at gaussian point: */
+		GetParameterValue(&dux, &dux_list[0],gauss_l1l2l3);
+		GetParameterValue(&duy, &duy_list[0],gauss_l1l2l3);
+
+		/*compute Du*/
+		for (i=0;i<numgrids;i++){
+			pe_g_gaussian[i*NDOF2+0]=dux*Jdet*gauss_weight*l1l2l3[i]; 
+			pe_g_gaussian[i*NDOF2+1]=duy*Jdet*gauss_weight*l1l2l3[i]; 
+		}
+
+		/*Add pe_g_gaussian vector to pe_g: */
+		for( i=0; i<numdof; i++){
+			pe_g[i]+=pe_g_gaussian[i];
+		}
+	}
+
+	/*Add pe_g to global vector p_g: */
+	VecSetValues(p_g,numdof,doflist,(const double*)pe_g,ADD_VALUES);
+
+	/*Clean up*/
+	xfree((void**)&first_gauss_area_coord);
+	xfree((void**)&second_gauss_area_coord);
+	xfree((void**)&third_gauss_area_coord);
+	xfree((void**)&gauss_weights);
+
+}
+/*}}}*/
 /*FUNCTION Tria::CreatePVectorDiagnosticHutter{{{1*/
 void  Tria::CreatePVectorDiagnosticHutter(Vec pg){
Index: /issm/trunk/src/c/objects/Elements/Tria.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.h	(revision 4838)
+++ /issm/trunk/src/c/objects/Elements/Tria.h	(revision 4839)
@@ -70,5 +70,4 @@
 		void   CreateKMatrix(Mat Kgg);
 		void   CreatePVector(Vec pg);
-		void   Du(Vec du_g);
 		void   GetBedList(double* bed_list);
 		void*  GetMatPar();
@@ -128,4 +127,6 @@
 		void	  CreatePVectorDiagnosticBaseVert(Vec pg);
 		void	  CreatePVectorDiagnosticHoriz(Vec pg);
+		void	  CreatePVectorAdjointHoriz(Vec pg);
+		void	  CreatePVectorAdjointStokes(Vec pg);
 		void	  CreatePVectorDiagnosticHutter(Vec pg);
 		void	  CreatePVectorPrognostic_CG(Vec pg);
Index: /issm/trunk/src/c/objects/Loads/Icefront.cpp
===================================================================
--- /issm/trunk/src/c/objects/Loads/Icefront.cpp	(revision 4838)
+++ /issm/trunk/src/c/objects/Loads/Icefront.cpp	(revision 4839)
@@ -312,8 +312,5 @@
 
 	/*Just branch to the correct element icefront vector generator, according to the type of analysis we are carrying out: */
-	if (analysis_type==ControlAnalysisEnum){
-		CreatePVectorDiagnosticHoriz( pg);
-	}
-	else if (analysis_type==DiagnosticHorizAnalysisEnum){
+	if (analysis_type==DiagnosticHorizAnalysisEnum){
 		CreatePVectorDiagnosticHoriz( pg);
 	}
Index: /issm/trunk/src/c/solutions/adjoint_core.cpp
===================================================================
--- /issm/trunk/src/c/solutions/adjoint_core.cpp	(revision 4838)
+++ /issm/trunk/src/c/solutions/adjoint_core.cpp	(revision 4839)
@@ -14,81 +14,35 @@
 void adjoint_core(FemModel* femmodel){
 	
-	
 	/*parameters: */
-	int   verbose=0;
-	char* solverstring=NULL;
-	bool  isstokes=false;
-	bool  conserve_loads=true;
-	int   dim;
-	int   solution_type;
-	
-	/*intermediary: */
-	Vec u_g=NULL;
-	Mat K_ff0=NULL;
-	Mat K_fs0=NULL;
-
-	Vec du_g=NULL;
-	Vec du_f=NULL;
-
-	Vec adjoint_f=NULL;
-	Vec adjoint_g=NULL;
+	int  verbose        = 0;
+	bool isstokes;
+	bool conserve_loads = true;
+	int  solution_type;
 
 	/*retrieve parameters:*/
-	femmodel->parameters->FindParam(&solverstring,SolverStringEnum);
 	femmodel->parameters->FindParam(&verbose,VerboseEnum);
 	femmodel->parameters->FindParam(&isstokes,IsStokesEnum);
-	femmodel->parameters->FindParam(&dim,DimEnum);
 	femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
 
 	/*set analysis type to compute velocity: */
+	_printf_("%s\n","      computing velocities");
 	if(isstokes)femmodel->SetCurrentConfiguration(DiagnosticStokesAnalysisEnum);
 	else femmodel->SetCurrentConfiguration(DiagnosticHorizAnalysisEnum);
-	
-	_printf_("%s\n","      recover solution for this stiffness and right hand side:");
-	solver_diagnostic_nonlinear(NULL,&K_ff0,&K_fs0, femmodel,conserve_loads); 
-
-	_printf_("%s\n","      buid Du, difference between observed velocity and model velocity:");
-	Dux(&du_g, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
-
-	_printf_("%s\n","      reduce adjoint load from g-set to f-set:");
-	Reduceloadfromgtofx(&du_f, du_g, femmodel->Gmn, K_fs0, femmodel->ys, femmodel->nodesets,true); //true means that ys0 flag is activated: all spcs show 0 displacement
-	
-	/*free some ressources: */
-	VecFree(&du_g);MatFree(&K_fs0);
-
-	_printf_("%s\n","      solve for adjoint vector:");
-	Solverx(&adjoint_f, K_ff0, du_f, NULL, solverstring);
-	
-	/*free some ressources: */
-	VecFree(&du_f); MatFree(&K_ff0);
-	
-	_printf_("%s\n","      merge back to g set:");
-	Mergesolutionfromftogx(&adjoint_g, adjoint_f,femmodel->Gmn,femmodel->ys,femmodel->nodesets,true);//true means that ys0 flag is activated: all spc are 0
-	
-	/*free some ressources: */
-	VecFree(&adjoint_f);
+	solver_diagnostic_nonlinear(femmodel,conserve_loads); 
 
 	/*Update inputs using adjoint solution, and same type of setup as diagnostic solution: */
+	_printf_("%s\n","      computing adjoint");
 	if(isstokes)femmodel->SetCurrentConfiguration(DiagnosticStokesAnalysisEnum,AdjointStokesAnalysisEnum);
 	else femmodel->SetCurrentConfiguration(DiagnosticHorizAnalysisEnum,AdjointHorizAnalysisEnum);
+	solver_adjoint_linear(femmodel);
 	
-	InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,adjoint_g);
-
 	if(verbose)_printf_("saving results:\n");
 	if(solution_type==AdjointSolutionEnum){
 		InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,AdjointxEnum);
 		InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,AdjointyEnum);
-		if(dim==3) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,AdjointzEnum);
+		if (isstokes){
+			InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,AdjointzEnum);
+			InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,AdjointpEnum);
+		}
 	}
-	
-	/*Free ressources:*/
-	xfree((void**)&solverstring);
-	VecFree(&u_g);
-	MatFree(&K_ff0);
-	MatFree(&K_fs0);
-	VecFree(&du_g);
-	VecFree(&du_f);
-	VecFree(&adjoint_f);
-	VecFree(&adjoint_g);
-
 }
Index: /issm/trunk/src/c/solutions/diagnostic_core.cpp
===================================================================
--- /issm/trunk/src/c/solutions/diagnostic_core.cpp	(revision 4838)
+++ /issm/trunk/src/c/solutions/diagnostic_core.cpp	(revision 4839)
@@ -15,14 +15,14 @@
 
 	/*parameters: */
-	int  verbose=0;
-	bool qmu_analysis=false;
-	int  dim=-1;
-	bool ishutter=false;
-	bool ismacayealpattyn=false;
-	bool isstokes=false;
+	int    verbose              = 0;
+	bool   qmu_analysis         = false;
+	int    dim                  = -1;
+	bool   ishutter             = false;
+	bool   ismacayealpattyn     = false;
+	bool   isstokes             = false;
 	double stokesreconditioning;
-	bool conserve_loads=true;
-	bool modify_loads=true;
-	int solution_type;
+	bool   conserve_loads       = true;
+	bool   modify_loads         = true;
+	int    solution_type;
 
 	/* recover parameters:*/
@@ -65,5 +65,5 @@
 		if(verbose)_printf_("%s\n"," computing horizontal velocities...");
 		femmodel->SetCurrentConfiguration(DiagnosticHorizAnalysisEnum);
-		solver_diagnostic_nonlinear(NULL,NULL,NULL,femmodel,modify_loads); 
+		solver_diagnostic_nonlinear(femmodel,modify_loads); 
 	}
 	
@@ -82,5 +82,5 @@
 			if(verbose)_printf_("%s\n"," computing stokes velocities and pressure ...");
 			femmodel->SetCurrentConfiguration(DiagnosticStokesAnalysisEnum);
-			solver_diagnostic_nonlinear(NULL,NULL,NULL,femmodel,conserve_loads);
+			solver_diagnostic_nonlinear(femmodel,conserve_loads);
 		}
 	}
Index: /issm/trunk/src/c/solutions/objectivefunctionC.cpp
===================================================================
--- /issm/trunk/src/c/solutions/objectivefunctionC.cpp	(revision 4838)
+++ /issm/trunk/src/c/solutions/objectivefunctionC.cpp	(revision 4839)
@@ -69,5 +69,5 @@
 
 	/*Run diagnostic with updated inputs: */
-	if(!control_steady) solver_diagnostic_nonlinear(NULL,NULL,NULL,femmodel,conserve_loads);  //true means we conserve loads at each diagnostic run
+	if(!control_steady) solver_diagnostic_nonlinear(femmodel,conserve_loads);  //true means we conserve loads at each diagnostic run
 	else                diagnostic_core(femmodel);	//We need a 3D velocity!! (vz is required for the next thermal run)
 
Index: /issm/trunk/src/c/solutions/stokescontrolinit.cpp
===================================================================
--- /issm/trunk/src/c/solutions/stokescontrolinit.cpp	(revision 4838)
+++ /issm/trunk/src/c/solutions/stokescontrolinit.cpp	(revision 4839)
@@ -41,5 +41,5 @@
 	/*Run a complete diagnostic to update the Stokes spcs: */
 	femmodel->SetCurrentConfiguration(DiagnosticHorizAnalysisEnum);
-	solver_diagnostic_nonlinear(NULL,NULL,NULL,femmodel,conserve_loads);
+	solver_diagnostic_nonlinear(femmodel,conserve_loads);
 	femmodel->SetCurrentConfiguration(DiagnosticVertAnalysisEnum);
 	solver_linear(femmodel);
@@ -50,4 +50,4 @@
 	if(verbose)_printf_("%s\n"," computing stokes velocities and pressure ...");
 	femmodel->SetCurrentConfiguration(DiagnosticStokesAnalysisEnum);
-	solver_diagnostic_nonlinear(NULL,NULL,NULL,femmodel,conserve_loads);
+	solver_diagnostic_nonlinear(femmodel,conserve_loads);
 }
Index: /issm/trunk/src/c/solutions/thermal_core_step.cpp
===================================================================
--- /issm/trunk/src/c/solutions/thermal_core_step.cpp	(revision 4838)
+++ /issm/trunk/src/c/solutions/thermal_core_step.cpp	(revision 4839)
@@ -21,5 +21,5 @@
 	if(verbose)_printf_("computing temperatures:\n");
 	femmodel->SetCurrentConfiguration(ThermalAnalysisEnum);
-	solver_thermal_nonlinear(NULL,NULL,femmodel);
+	solver_thermal_nonlinear(femmodel);
 
 	if(verbose)_printf_("computing melting:\n");
Index: /issm/trunk/src/c/solvers/solver_adjoint_linear.cpp
===================================================================
--- /issm/trunk/src/c/solvers/solver_adjoint_linear.cpp	(revision 4839)
+++ /issm/trunk/src/c/solvers/solver_adjoint_linear.cpp	(revision 4839)
@@ -0,0 +1,66 @@
+/*!\file: solver_adjoint_linear.cpp
+ * \brief: numerical core of linear solutions
+ */ 
+
+#include "../toolkits/toolkits.h"
+#include "../objects/objects.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+#include "../modules/modules.h"
+
+void solver_adjoint_linear(FemModel* fem){
+	/*This is axactly the same solver as solver_linear except that Reduceloadfromgtofx and Mergesolutionfromftogx
+	 * use the flag "true" so that all spc are taken as 0*/
+
+	/*parameters:*/
+	int kflag,pflag;
+	int verbose=0;
+	char* solver_string=NULL;
+
+	/*output: */
+	Vec ug=NULL;
+	Vec uf=NULL; 
+	
+	/*intermediary: */
+	Mat Kgg=NULL;
+	Mat Kff=NULL;
+	Mat Kfs=NULL;
+	Vec pg=NULL;
+	Vec pf=NULL;
+
+	/*Recover parameters: */
+	kflag=1; pflag=1;
+	fem->parameters->FindParam(&verbose,VerboseEnum);
+	fem->parameters->FindParam(&solver_string,SolverStringEnum);
+
+	//*Generate system matrices
+	if(verbose)_printf_("   Generating matrices\n");
+	SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag);
+
+	if(verbose)_printf_("   Generating penalty matrices\n");
+	//*Generate penalty system matrices
+	PenaltySystemMatricesx(Kgg, pg,NULL,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag);
+
+	/*!Reduce matrix from g to f size:*/
+	if(verbose)_printf_("   reducing matrix from g to f set\n");
+	Reducematrixfromgtofx(&Kff,&Kfs,Kgg,fem->Gmn,fem->nodesets); MatFree(&Kgg);
+
+	/*!Reduce load from g to f size: */
+	if(verbose)_printf_("   reducing load from g to f set\n"); //true means spc = 0
+	Reduceloadfromgtofx(&pf, pg, fem->Gmn, Kfs, fem->ys, fem->nodesets,true);VecFree(&pg); MatFree(&Kfs);
+
+	/*Solve: */
+	if(verbose)_printf_("   solving\n");
+	Solverx(&uf, Kff, pf, NULL, solver_string); MatFree(&Kff); VecFree(&pf);
+
+	//Merge back to g set
+	if(verbose)_printf_("   merging solution from f to g set\n");//true means spc=0
+	Mergesolutionfromftogx(&ug, uf,fem->Gmn,fem->ys,fem->nodesets,true);VecFree(&uf);
+
+	//Update inputs using new solution:
+	InputUpdateFromSolutionx( fem->elements,fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters,ug);
+
+	/*free ressources: */
+	xfree((void**)&solver_string);
+	VecFree(&ug);
+	VecFree(&uf);
+}
Index: /issm/trunk/src/c/solvers/solver_diagnostic_nonlinear.cpp
===================================================================
--- /issm/trunk/src/c/solvers/solver_diagnostic_nonlinear.cpp	(revision 4838)
+++ /issm/trunk/src/c/solvers/solver_diagnostic_nonlinear.cpp	(revision 4839)
@@ -10,5 +10,5 @@
 #include "./solvers.h"
 
-void solver_diagnostic_nonlinear(Vec* pug,Mat* pKff0,Mat* pKfs0, FemModel* femmodel,bool conserve_loads){
+void solver_diagnostic_nonlinear(FemModel* femmodel,bool conserve_loads){
 
 
@@ -127,15 +127,4 @@
 		}
 	}
-	
-	//more output might be needed, when running in control
-	if(pKff0){
-
-		kflag=1; pflag=0; //stiffness generation only
-	
-		SystemMatricesx(&Kgg, &pg,femmodel->elements,femmodel->nodes,femmodel->vertices,loads,femmodel->materials,femmodel->parameters,kflag,pflag);
-		Reducematrixfromgtofx(&Kff,&Kfs,Kgg,femmodel->Gmn,femmodel->nodesets);
-		MatFree(&Kgg);VecFree(&pg);
-
-	}
 
 	/*Delete loads only if no ouput was requested: */
@@ -144,14 +133,8 @@
 	/*clean up*/
 	VecFree(&uf);
+	VecFree(&ug);
 	VecFree(&old_uf);
 	VecFree(&old_ug);
 	xfree((void**)&solver_string);
 	
-	/*Assign output pointers: */
-	if(pug)*pug=ug;
-	else VecFree(&ug);
-
-	if(pKff0)*pKff0=Kff;
-	if(pKfs0)*pKfs0=Kfs;
-
 }
Index: /issm/trunk/src/c/solvers/solver_thermal_nonlinear.cpp
===================================================================
--- /issm/trunk/src/c/solvers/solver_thermal_nonlinear.cpp	(revision 4838)
+++ /issm/trunk/src/c/solvers/solver_thermal_nonlinear.cpp	(revision 4839)
@@ -8,5 +8,5 @@
 #include "../modules/modules.h"
 
-void solver_thermal_nonlinear(Vec* ptg,double* pmelting_offset,FemModel* fem){
+void solver_thermal_nonlinear(FemModel* fem){
 
 	/*solution : */
@@ -125,11 +125,7 @@
 	MatFree(&Kgg_nopenalty);
 	VecFree(&pg_nopenalty);
+	VecFree(&tg);
 	VecFree(&tf);
 	VecFree(&tf_old);
 	xfree((void**)&solver_string);
-
-	/*Assign output pointers: */
-	if(ptg)*ptg=tg;
-	else VecFree(&tg);
-	if(pmelting_offset)*pmelting_offset=melting_offset;
 }
Index: /issm/trunk/src/c/solvers/solvers.h
===================================================================
--- /issm/trunk/src/c/solvers/solvers.h	(revision 4838)
+++ /issm/trunk/src/c/solvers/solvers.h	(revision 4839)
@@ -12,7 +12,8 @@
 class FemModel;
 
-void solver_thermal_nonlinear(Vec* ptg,double* pmelting_offset,FemModel* femmodel);
-void solver_diagnostic_nonlinear(Vec* pug,Mat* pK_ff0,Mat* pK_fs0, FemModel* femmodel,bool conserve_loads);
+void solver_thermal_nonlinear(FemModel* femmodel);
+void solver_diagnostic_nonlinear(FemModel* femmodel,bool conserve_loads);
 void solver_linear(FemModel* femmodel);
+void solver_adjoint_linear(FemModel* femmodel);
 
 #endif
Index: /issm/trunk/src/m/solutions/adjoint_core.m
===================================================================
--- /issm/trunk/src/m/solutions/adjoint_core.m	(revision 4838)
+++ /issm/trunk/src/m/solutions/adjoint_core.m	(revision 4839)
@@ -10,4 +10,5 @@
 
 	%set analysis type to compute velocity:
+	displaystring('\n%s',['      computing velocities']);
 	if(isstokes),
 		femmodel=SetCurrentConfiguration(femmodel,DiagnosticStokesAnalysisEnum);
@@ -15,21 +16,7 @@
 		femmodel=SetCurrentConfiguration(femmodel,DiagnosticHorizAnalysisEnum);
 	end
+	femmodel=solver_diagnostic_nonlinear(femmodel,conserve_loads);
 
-	displaystring('\n%s',['      recover solution for this stiffness and right hand side:']);
-	[femmodel,ug,K_ff0,K_fs0]=solver_diagnostic_nonlinear(femmodel,conserve_loads);
-
-	displaystring('\n%s',['      Build Du, difference between observed and modeled velocities:']);
-	du_g=Du(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
-
-	displaystring('\n%s',['      reduce adjoint load from g_set to f_set:']);
-	du_f=Reduceloadfromgtof(du_g,femmodel.Gmn,K_fs0,femmodel.ys,femmodel.nodesets,true); %true means that ys0 flag is activated: all spcs show 0 displacement
-
-	displaystring('\n%s',['      solve for adjoint vector:']);
-	adjoint_f=Solver(K_ff0,du_f,[],femmodel.parameters);
-
-	displaystring('\n%s',['      merge back to g set:']);
-	adjoint_g=Mergesolutionfromftog(adjoint_f,femmodel.Gmn,femmodel.ys,femmodel.nodesets,true);%true means that ys0 flag is activated: all spc are 0
-
-	%Update inputs using adjoint solution, and same type of setup as diagnostic solution:
+	displaystring('\n%s',['      computing asjoint']);
 	if(isstokes),
 		femmodel=SetCurrentConfiguration(femmodel,DiagnosticStokesAnalysisEnum,AdjointStokesAnalysisEnum);
@@ -37,6 +24,5 @@
 		femmodel=SetCurrentConfiguration(femmodel,DiagnosticHorizAnalysisEnum,AdjointHorizAnalysisEnum);
 	end
-
-	[femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,adjoint_g);
+	femmodel=solver_adjoint_linear(femmodel);
 
 	displaystring(verbose,'\n%s',['      saving results...']);
@@ -44,6 +30,7 @@
 		femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,AdjointxEnum);
 		femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,AdjointyEnum);
-		if(dim==3),
+		if(isstokes),
 			femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,AdjointzEnum);
+			femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,AdjointpEnum);
 		end
 	end
Index: /issm/trunk/src/m/solvers/solver_adjoint_linear.m
===================================================================
--- /issm/trunk/src/m/solvers/solver_adjoint_linear.m	(revision 4839)
+++ /issm/trunk/src/m/solvers/solver_adjoint_linear.m	(revision 4839)
@@ -0,0 +1,29 @@
+function femmodel=solver_adjoint_linear(femmodel)
+%SOLVER_LINEAR - core solver of any linear solution sequence
+%
+%   Usage:
+%      femmodel =solver_adjoint_linear(femmodel)
+
+	%stiffness and load generation only:
+	femmodel.parameters.Kflag=1; femmodel.parameters.Pflag=1;
+
+	%system matrices
+	[K_gg, p_g]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
+	[K_gg, p_g,kmax]=PenaltySystemMatrices(K_gg,p_g,femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
+	
+	%Reduce tangent matrix from g size to f size
+	[K_ff, K_fs] = Reducematrixfromgtof( K_gg, femmodel.Gmn, femmodel.nodesets); 
+	displaystring(femmodel.parameters.Verbose>1,'%s%g','      condition number of stiffness matrix: ',condest(K_ff));
+	
+	%Reduce load from g size to f size
+	p_f = Reduceloadfromgtof( p_g, femmodel.Gmn, K_fs, femmodel.ys, femmodel.nodesets,true);
+
+	%Solve	
+	u_f=Solver(K_ff,p_f,[],femmodel.parameters);
+	
+	%Merge back to g set
+	u_g= Mergesolutionfromftog( u_f, femmodel.Gmn, femmodel.ys, femmodel.nodesets,true); 
+
+	%Update inputs using new solution
+	[femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,u_g);
+
Index: /issm/trunk/src/m/solvers/solver_diagnostic_nonlinear.m
===================================================================
--- /issm/trunk/src/m/solvers/solver_diagnostic_nonlinear.m	(revision 4838)
+++ /issm/trunk/src/m/solvers/solver_diagnostic_nonlinear.m	(revision 4839)
@@ -1,7 +1,7 @@
-function [femmodel ug varargout]=solver_diagnostic_nonlinear(femmodel,conserve_loads)
+function femmodel=solver_diagnostic_nonlinear(femmodel,conserve_loads)
 %SOLVER_DIAGNOSTIC_NONLINEAR - core solver of diagnostic run
 %
 %   Usage:
-%      [femmodel ug varargout]=solver_diagnostic_nonlinear(femmodel,conserve_loads)
+%      [femmodel]=solver_diagnostic_nonlinear(femmodel,conserve_loads)
 	
 
@@ -82,15 +82,4 @@
 	end
 
-	%more output might be needed, when running in control_core.m
-	nout=max(nargout,1)-2;
-	if nout==2,
-		femmodel.parameters.Kflag=1; femmodel.parameters.Pflag=0;
-		[K_gg, p_g]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters);
-		[K_ff, K_fs] = Reducematrixfromgtof( K_gg, femmodel.Gmn, femmodel.nodesets); 
-		varargout(1)={K_ff};
-		varargout(2)={K_fs};
-	end
-
-
 	%deal with loads:
 	if conserve_loads==false,
Index: /issm/trunk/src/m/solvers/solver_linear.m
===================================================================
--- /issm/trunk/src/m/solvers/solver_linear.m	(revision 4838)
+++ /issm/trunk/src/m/solvers/solver_linear.m	(revision 4839)
@@ -1,7 +1,7 @@
-function [femmodel u_g]=solver_linear(femmodel)
+function femmodel=solver_linear(femmodel)
 %SOLVER_LINEAR - core solver of any linear solution sequence
 %
 %   Usage:
-%      [femmodel, u_g]=solver_linear(femmodel)
+%      femmodel =solver_linear(femmodel)
 
 	%stiffness and load generation only:
Index: /issm/trunk/src/m/solvers/solver_thermal_nonlinear.m
===================================================================
--- /issm/trunk/src/m/solvers/solver_thermal_nonlinear.m	(revision 4838)
+++ /issm/trunk/src/m/solvers/solver_thermal_nonlinear.m	(revision 4839)
@@ -1,8 +1,8 @@
-function [femmodel t_g ,melting_offset]=solver_thermal_nonlinear(femmodel)
+function femmodel=solver_thermal_nonlinear(femmodel)
 %SOLVER_THERMAL_NONLINEAR - core of thermal solution sequence.
 %   femmodel is returned together with temperature and melting_offset, in case loads have been modified
 %
 %   Usage:
-%      [femmodel t_g melting_offset]=solver_thermal_nonlinear(femmodel)
+%      [femmodel]=solver_thermal_nonlinear(femmodel)
 
 	count=1;
Index: /issm/trunk/src/mex/Makefile.am
===================================================================
--- /issm/trunk/src/mex/Makefile.am	(revision 4838)
+++ /issm/trunk/src/mex/Makefile.am	(revision 4839)
@@ -19,5 +19,4 @@
 				CostFunction \
 				DakotaResponses\
-				Du\
 				Echo\
 				ElementConnectivity\
@@ -157,7 +156,4 @@
 			  DakotaResponses/DakotaResponses.h
 
-Du_SOURCES = Du/Du.cpp\
-			  Du/Du.h
-
 Echo_SOURCES = Echo/Echo.cpp\
 			  Echo/Echo.h