#include "./AdjointHorizAnalysis.h"
#include "../toolkits/toolkits.h"
#include "../classes/classes.h"
#include "../shared/shared.h"
#include "../modules/modules.h"

/*Model processing*/
int  AdjointHorizAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
	_error_("not implemented");
}/*}}}*/
void AdjointHorizAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
	   _error_("not implemented yet");
}/*}}}*/
void AdjointHorizAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
	   _error_("not implemented yet");
}/*}}}*/
void AdjointHorizAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
	   _error_("not implemented yet");
}/*}}}*/
void AdjointHorizAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
	   _error_("not implemented yet");
}/*}}}*/
void AdjointHorizAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
	   _error_("not implemented yet");
}/*}}}*/

/*Finite Element Analysis*/
void           AdjointHorizAnalysis::Core(FemModel* femmodel){/*{{{*/
	_error_("not implemented");
}/*}}}*/
ElementVector* AdjointHorizAnalysis::CreateDVector(Element* element){/*{{{*/
	/*Default, return NULL*/
	return NULL;
}/*}}}*/
ElementMatrix* AdjointHorizAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
_error_("Not implemented");
}/*}}}*/
ElementMatrix* AdjointHorizAnalysis::CreateKMatrix(Element* element){/*{{{*/
	int approximation;
	element->GetInputValue(&approximation,ApproximationEnum);
	switch(approximation){
		case SSAApproximationEnum: 
			return CreateKMatrixSSA(element);
		case HOApproximationEnum: 
			return CreateKMatrixHO(element);
		case FSApproximationEnum: 
			return CreateKMatrixFS(element);
		case NoneApproximationEnum:
			return NULL;
		default:
			_error_("Approximation "<<EnumToStringx(approximation)<<" not supported");
	}
}/*}}}*/
ElementMatrix* AdjointHorizAnalysis::CreateKMatrixSSA(Element* element){/*{{{*/

	/*Intermediaries*/
	int      domaintype;
	Element* basalelement;

	/*Get basal element*/
	element->FindParam(&domaintype,DomainTypeEnum);
	switch(domaintype){
		case Domain2DhorizontalEnum:
			basalelement = element;
			break;
		case Domain3DEnum:
			if(!element->IsOnBase()) return NULL;
			basalelement = element->SpawnBasalElement();
			break;
		default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
	}

	/*Intermediaries */
	bool        incomplete_adjoint;
	IssmDouble  Jdet,thickness,mu_prime;
	IssmDouble  eps1dotdphii,eps1dotdphij,eps2dotdphii,eps2dotdphij;
	IssmDouble  eps1[2],eps2[2],epsilon[3];
	IssmDouble *xyz_list = NULL;

	/*Fetch number of nodes and dof for this finite element*/
	int numnodes = basalelement->GetNumberOfNodes();

	/*Initialize Jacobian with regular SSA (first part of the Gateau derivative)*/
	basalelement->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum);
	StressbalanceAnalysis* analysis = new StressbalanceAnalysis();
	ElementMatrix* Ke=analysis->CreateKMatrix(element);
	delete analysis;
	if(incomplete_adjoint){
		if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
		return Ke;
	}

	/*Retrieve all inputs and parameters*/
	basalelement->GetVerticesCoordinates(&xyz_list);
	Input* vx_input        = basalelement->GetInput(VxEnum);       _assert_(vx_input);
	Input* vy_input        = basalelement->GetInput(VyEnum);       _assert_(vy_input);
	Input* thickness_input = basalelement->GetInput(ThicknessEnum); _assert_(thickness_input);

	/*Allocate dbasis*/
	IssmDouble* dbasis = xNew<IssmDouble>(2*numnodes);

	/* Start  looping on the number of gaussian points: */
	Gauss* gauss=basalelement->NewGauss(2);
	for(int ig=gauss->begin();ig<gauss->end();ig++){
		gauss->GaussPoint(ig);

		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
		basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);

		thickness_input->GetInputValue(&thickness, gauss);
		basalelement->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
		basalelement->ViscositySSADerivativeEpsSquare(&mu_prime,&epsilon[0]);
		eps1[0]=2.*epsilon[0]+epsilon[1]; eps2[0]=epsilon[2];
		eps1[1]=epsilon[2];               eps2[1]=epsilon[0]+2*epsilon[1];

		for(int i=0;i<numnodes;i++){
			for(int j=0;j<numnodes;j++){
				eps1dotdphii=eps1[0]*dbasis[0*numnodes+i]+eps1[1]*dbasis[1*numnodes+i];
				eps1dotdphij=eps1[0]*dbasis[0*numnodes+j]+eps1[1]*dbasis[1*numnodes+j];
				eps2dotdphii=eps2[0]*dbasis[0*numnodes+i]+eps2[1]*dbasis[1*numnodes+i];
				eps2dotdphij=eps2[0]*dbasis[0*numnodes+j]+eps2[1]*dbasis[1*numnodes+j];

				Ke->values[2*numnodes*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps1dotdphii;
				Ke->values[2*numnodes*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps1dotdphii;
				Ke->values[2*numnodes*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps2dotdphii;
				Ke->values[2*numnodes*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps2dotdphii;
			}
		}
	}

	/*Transform Coordinate System*/
	basalelement->TransformStiffnessMatrixCoord(Ke,XYEnum);

	/*Clean up and return*/
	delete gauss;
	xDelete<IssmDouble>(dbasis);
	xDelete<IssmDouble>(xyz_list);
	if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
	return Ke;
}/*}}}*/
ElementMatrix* AdjointHorizAnalysis::CreateKMatrixHO(Element* element){/*{{{*/

	/*Intermediaries */
	bool        incomplete_adjoint;
	IssmDouble  Jdet,mu_prime;
	IssmDouble  eps1dotdphii,eps1dotdphij,eps2dotdphii,eps2dotdphij;
	IssmDouble  eps1[3],eps2[3],epsilon[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/
	IssmDouble *xyz_list = NULL;

	/*Fetch number of nodes and dof for this finite element*/
	int numnodes = element->GetNumberOfNodes();

	/*Initialize Jacobian with regular HO (first part of the Gateau derivative)*/
	element->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum);
	StressbalanceAnalysis* analysis = new StressbalanceAnalysis();
	ElementMatrix* Ke=analysis->CreateKMatrix(element);
	delete analysis;
	if(incomplete_adjoint) return Ke;

	/*Retrieve all inputs and parameters*/
	element->GetVerticesCoordinates(&xyz_list);
	Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input);
	Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input);

	/*Allocate dbasis*/
	IssmDouble* dbasis = xNew<IssmDouble>(3*numnodes);

	/* Start  looping on the number of gaussian points: */
	Gauss* gauss=element->NewGauss(5);
	for(int ig=gauss->begin();ig<gauss->end();ig++){
		gauss->GaussPoint(ig);

		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
		element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);

		element->StrainRateHO(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
		element->ViscosityHODerivativeEpsSquare(&mu_prime,&epsilon[0]);
		eps1[0]=2.*epsilon[0]+epsilon[1];   eps2[0]=epsilon[2];
		eps1[1]=epsilon[2];                 eps2[1]=epsilon[0]+2.*epsilon[1];
		eps1[2]=epsilon[3];                 eps2[2]=epsilon[4];

		for(int i=0;i<numnodes;i++){
			for(int j=0;j<numnodes;j++){
				eps1dotdphii=eps1[0]*dbasis[0*numnodes+i]+eps1[1]*dbasis[1*numnodes+i]+eps1[2]*dbasis[2*numnodes+i];
				eps1dotdphij=eps1[0]*dbasis[0*numnodes+j]+eps1[1]*dbasis[1*numnodes+j]+eps1[2]*dbasis[2*numnodes+j];
				eps2dotdphii=eps2[0]*dbasis[0*numnodes+i]+eps2[1]*dbasis[1*numnodes+i]+eps2[2]*dbasis[2*numnodes+i];
				eps2dotdphij=eps2[0]*dbasis[0*numnodes+j]+eps2[1]*dbasis[1*numnodes+j]+eps2[2]*dbasis[2*numnodes+j];

				Ke->values[2*numnodes*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii;
				Ke->values[2*numnodes*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii;
				Ke->values[2*numnodes*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii;
				Ke->values[2*numnodes*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii;
			}
		}
	}

	/*Transform Coordinate System*/
	element->TransformStiffnessMatrixCoord(Ke,XYEnum);

	/*Clean up and return*/
	delete gauss;
	xDelete<IssmDouble>(dbasis);
	xDelete<IssmDouble>(xyz_list);
	return Ke;
}/*}}}*/
ElementMatrix* AdjointHorizAnalysis::CreateKMatrixFS(Element* element){/*{{{*/

	/*Intermediaries */
	bool        incomplete_adjoint;
	int         dim,epssize;
	IssmDouble  Jdet,mu_prime;
	IssmDouble  eps1dotdphii,eps1dotdphij,eps2dotdphii,eps2dotdphij,eps3dotdphii,eps3dotdphij;
	IssmDouble  eps1[3],eps2[3],eps3[3],epsilon[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/
	IssmDouble *xyz_list = NULL;

	/*Get problem dimension*/
	element->FindParam(&dim,DomainDimensionEnum);
	if(dim==2) epssize = 3;
	else       epssize = 6;

	/*Fetch number of nodes and dof for this finite element*/
	int vnumnodes = element->NumberofNodesVelocity();
	int pnumnodes = element->NumberofNodesPressure();
	int numdof    = vnumnodes*dim + pnumnodes;

	/*Initialize Jacobian with regular FS (first part of the Gateau derivative)*/
	element->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum);
	StressbalanceAnalysis* analysis = new StressbalanceAnalysis();
	ElementMatrix* Ke=analysis->CreateKMatrix(element);
	delete analysis;
	if(incomplete_adjoint) return Ke;

	/*Retrieve all inputs and parameters*/
	element->GetVerticesCoordinates(&xyz_list);
	Input* vx_input = element->GetInput(VxEnum);_assert_(vx_input);
	Input* vy_input = element->GetInput(VyEnum);_assert_(vy_input);
	Input* vz_input = NULL;
	if(dim==3){
		vz_input = element->GetInput(VzEnum);
	}
	else{
		_error_("Not implemented yet");
	}

	/*Allocate dbasis*/
	IssmDouble* dbasis = xNew<IssmDouble>(dim*vnumnodes);

	/* Start  looping on the number of gaussian points: */
	Gauss* gauss=element->NewGauss(5);
	for(int ig=gauss->begin();ig<gauss->end();ig++){
		gauss->GaussPoint(ig);

		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
		element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);

		element->StrainRateHO(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
		element->ViscosityFSDerivativeEpsSquare(&mu_prime,&epsilon[0]);
		eps1[0]=epsilon[0];   eps2[0]=epsilon[2];   eps3[0]=epsilon[3];
		eps1[1]=epsilon[2];   eps2[1]=epsilon[1];   eps3[1]=epsilon[4];
		eps1[2]=epsilon[3];   eps2[2]=epsilon[4];   eps3[2]= -epsilon[0] -epsilon[1];

		for(int i=0;i<vnumnodes;i++){
			for(int j=0;j<vnumnodes;j++){
				eps1dotdphii=eps1[0]*dbasis[0*vnumnodes+i]+eps1[1]*dbasis[1*vnumnodes+i]+eps1[2]*dbasis[2*vnumnodes+i];
				eps1dotdphij=eps1[0]*dbasis[0*vnumnodes+j]+eps1[1]*dbasis[1*vnumnodes+j]+eps1[2]*dbasis[2*vnumnodes+j];
				eps2dotdphii=eps2[0]*dbasis[0*vnumnodes+i]+eps2[1]*dbasis[1*vnumnodes+i]+eps2[2]*dbasis[2*vnumnodes+i];
				eps2dotdphij=eps2[0]*dbasis[0*vnumnodes+j]+eps2[1]*dbasis[1*vnumnodes+j]+eps2[2]*dbasis[2*vnumnodes+j];
				eps3dotdphii=eps3[0]*dbasis[0*vnumnodes+i]+eps3[1]*dbasis[1*vnumnodes+i]+eps3[2]*dbasis[2*vnumnodes+i];
				eps3dotdphij=eps3[0]*dbasis[0*vnumnodes+j]+eps3[1]*dbasis[1*vnumnodes+j]+eps3[2]*dbasis[2*vnumnodes+j];

				Ke->values[numdof*(4*i+0)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii;
				Ke->values[numdof*(4*i+0)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii;
				Ke->values[numdof*(4*i+0)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps1dotdphii;

				Ke->values[numdof*(4*i+1)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii;
				Ke->values[numdof*(4*i+1)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii;
				Ke->values[numdof*(4*i+1)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps2dotdphii;

				Ke->values[numdof*(4*i+2)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps3dotdphii;
				Ke->values[numdof*(4*i+2)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps3dotdphii;
				Ke->values[numdof*(4*i+2)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps3dotdphii;
			}
		}
	}

	/*Transform Coordinate System*/
	element->TransformStiffnessMatrixCoord(Ke,XYZEnum);

	/*Clean up and return*/
	delete gauss;
	xDelete<IssmDouble>(dbasis);
	xDelete<IssmDouble>(xyz_list);
	return Ke;
}/*}}}*/
ElementVector* AdjointHorizAnalysis::CreatePVector(Element* element){/*{{{*/

	int approximation;
	element->GetInputValue(&approximation,ApproximationEnum);
	switch(approximation){
		case SSAApproximationEnum: 
			return CreatePVectorSSA(element);
		case HOApproximationEnum: 
			return CreatePVectorHO(element);
		case FSApproximationEnum: 
			return CreatePVectorFS(element);
		case NoneApproximationEnum:
			return NULL;
		default:
			_error_("Approximation "<<EnumToStringx(approximation)<<" not supported");
	}
}/*}}}*/
ElementVector* AdjointHorizAnalysis::CreatePVectorFS(Element* element){/*{{{*/

	/*Nothing to be done if not on surface*/
	if(!element->IsOnSurface()) return NULL;

	/*Intermediaries */
	int        num_responses,i,dim;
	IssmDouble Jdet,obs_velocity_mag,velocity_mag;
	IssmDouble vx,vy,vxobs,vyobs,dux,duy,weight;
	IssmDouble scalex,scaley,scale,S;
	int        *responses    = NULL;
	IssmDouble *xyz_list_top = NULL;

	/*Get problem dimension*/
	element->FindParam(&dim,DomainDimensionEnum);

	/*Fetch number of nodes and dof for this finite element*/
	int vnumnodes = element->NumberofNodesVelocity();
	int pnumnodes = element->NumberofNodesPressure();

	/*Prepare coordinate system list*/
	int* cs_list = xNew<int>(vnumnodes+pnumnodes);
	if(dim==2) for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
	else       for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum;
	for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;

	/*Initialize Element vector and vectors*/
	ElementVector* pe     = element->NewElementVector(FSApproximationEnum);
	IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);

	/*Retrieve all inputs and parameters*/
	element->GetVerticesCoordinatesTop(&xyz_list_top);
	element->FindParam(&num_responses,InversionNumCostFunctionsEnum);
	element->FindParam(&responses,NULL,InversionCostFunctionsEnum);
	Input* weights_input = element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
	Input* vx_input      = element->GetInput(VxEnum);                                 _assert_(vx_input);
	Input* vy_input      = element->GetInput(VyEnum);                                 _assert_(vy_input);
	Input* vxobs_input   = element->GetInput(InversionVxObsEnum);                     _assert_(vxobs_input);
	Input* vyobs_input   = element->GetInput(InversionVyObsEnum);                     _assert_(vyobs_input);
	IssmDouble epsvel  = 2.220446049250313e-16;
	IssmDouble meanvel = 3.170979198376458e-05; /*1000 m/yr*/

	/*Get Surface if required by one response*/
	for(int resp=0;resp<num_responses;resp++){
		if(responses[resp]==SurfaceAverageVelMisfitEnum){
			element->GetInputValue(&S,SurfaceAreaEnum); break;
		}
	}

	/* Start  looping on the number of gaussian points: */
	Gauss* gauss=element->NewGaussTop(4);
	for(int ig=gauss->begin();ig<gauss->end();ig++){
		gauss->GaussPoint(ig);

		element->JacobianDeterminantTop(&Jdet,xyz_list_top,gauss);
		element->NodalFunctionsVelocity(vbasis,gauss);

		vx_input->GetInputValue(&vx,gauss);
		vy_input->GetInputValue(&vy,gauss);
		vxobs_input->GetInputValue(&vxobs,gauss);
		vyobs_input->GetInputValue(&vyobs,gauss);

		/*Loop over all requested responses*/
		for(int resp=0;resp<num_responses;resp++){
			weights_input->GetInputValue(&weight,gauss,responses[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<vnumnodes;i++){
						dux=vxobs-vx;
						pe->values[i*dim+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; 
						if(dim==3){
							duy=vyobs-vy;
							pe->values[i*dim+1]+=duy*weight*Jdet*gauss->weight*vbasis[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<vnumnodes;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*dim+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; 
						pe->values[i*dim+1]+=duy*weight*Jdet*gauss->weight*vbasis[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<vnumnodes;i++){
						if(dim==3){
							velocity_mag    =sqrt(vx*vx+vy*vy)+epsvel;
							obs_velocity_mag=sqrt(vxobs*vxobs+vyobs*vyobs)+epsvel;
							scale=-8.*meanvel*meanvel/(velocity_mag*velocity_mag)*log(velocity_mag/obs_velocity_mag);
							dux=scale*vx;
							duy=scale*vy;
							pe->values[i*dim+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; 
							pe->values[i*dim+1]+=duy*weight*Jdet*gauss->weight*vbasis[i]; 
						}
						else{
							velocity_mag    =fabs(vx)+epsvel;
							obs_velocity_mag=fabs(vxobs)+epsvel;
							scale=-8.*meanvel*meanvel/(velocity_mag*velocity_mag)*log(velocity_mag/obs_velocity_mag);
							dux=scale*vx;
							pe->values[i*dim+0]+=dux*weight*Jdet*gauss->weight*vbasis[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<vnumnodes;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*dim+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; 
						pe->values[i*dim+1]+=duy*weight*Jdet*gauss->weight*vbasis[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<vnumnodes;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*dim+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; 
						pe->values[i*dim+1]+=duy*weight*Jdet*gauss->weight*vbasis[i]; 
					}
					break;
				case DragCoefficientAbsGradientEnum:
					/*Nothing in P vector*/
					break;
				case ThicknessAbsGradientEnum:
					/*Nothing in P vector*/
					break;
				case ThicknessAlongGradientEnum:
					/*Nothing in P vector*/
					break;
				case ThicknessAcrossGradientEnum:
					/*Nothing in P vector*/
					break;
				case RheologyBbarAbsGradientEnum:
					/*Nothing in P vector*/
					break;
				default:
					_error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
			}
		}
	}

	/*Transform coordinate system*/
	element->TransformLoadVectorCoord(pe,cs_list);

	/*Clean up and return*/
	xDelete<int>(cs_list);
	xDelete<int>(responses);
	xDelete<IssmDouble>(xyz_list_top);
	xDelete<IssmDouble>(vbasis);
	delete gauss;
	return pe;
}/*}}}*/
ElementVector* AdjointHorizAnalysis::CreatePVectorHO(Element* element){/*{{{*/

	/*Nothing to be done if not on surface*/
	if(!element->IsOnSurface()) return NULL;

	/*Intermediaries */
	int        num_responses,i;
	IssmDouble Jdet,obs_velocity_mag,velocity_mag;
	IssmDouble vx,vy,vxobs,vyobs,dux,duy,weight;
	IssmDouble scalex,scaley,scale,S;
	int        *responses    = NULL;
	IssmDouble *xyz_list_top = NULL;

	/*Fetch number of nodes and dof for this finite element*/
	int numnodes = element->GetNumberOfNodes();

	/*Initialize Element vector and vectors*/
	ElementVector* pe    = element->NewElementVector(HOApproximationEnum);
	IssmDouble*    basis = xNew<IssmDouble>(numnodes);

	/*Retrieve all inputs and parameters*/
	element->GetVerticesCoordinatesTop(&xyz_list_top);
	element->FindParam(&num_responses,InversionNumCostFunctionsEnum);
	element->FindParam(&responses,NULL,InversionCostFunctionsEnum);
	Input* weights_input = element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
	Input* vx_input      = element->GetInput(VxEnum);                                 _assert_(vx_input);
	Input* vy_input      = element->GetInput(VyEnum);                                 _assert_(vy_input);
	Input* vxobs_input   = element->GetInput(InversionVxObsEnum);                     _assert_(vxobs_input);
	Input* vyobs_input   = element->GetInput(InversionVyObsEnum);                     _assert_(vyobs_input);
	IssmDouble epsvel  = 2.220446049250313e-16;
	IssmDouble meanvel = 3.170979198376458e-05; /*1000 m/yr*/

	/*Get Surface if required by one response*/
	for(int resp=0;resp<num_responses;resp++){
		if(responses[resp]==SurfaceAverageVelMisfitEnum){
			element->GetInputValue(&S,SurfaceAreaEnum); break;
		}
	}

	/* Start  looping on the number of gaussian points: */
	Gauss* gauss=element->NewGaussTop(4);
	for(int ig=gauss->begin();ig<gauss->end();ig++){
		gauss->GaussPoint(ig);

		element->JacobianDeterminantTop(&Jdet,xyz_list_top,gauss);
		element->NodalFunctions(basis, gauss);

		vx_input->GetInputValue(&vx,gauss);
		vy_input->GetInputValue(&vy,gauss);
		vxobs_input->GetInputValue(&vxobs,gauss);
		vyobs_input->GetInputValue(&vyobs,gauss);

		/*Loop over all requested responses*/
		for(int resp=0;resp<num_responses;resp++){
			weights_input->GetInputValue(&weight,gauss,responses[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<numnodes;i++){
						dux=vxobs-vx;
						duy=vyobs-vy;
						pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
						pe->values[i*2+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<numnodes;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*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
						pe->values[i*2+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<numnodes;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*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
						pe->values[i*2+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<numnodes;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*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
						pe->values[i*2+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<numnodes;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*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
						pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
					}
					break;
				case DragCoefficientAbsGradientEnum:
					/*Nothing in P vector*/
					break;
				case ThicknessAbsGradientEnum:
					/*Nothing in P vector*/
					break;
				case ThicknessAlongGradientEnum:
					/*Nothing in P vector*/
					break;
				case ThicknessAcrossGradientEnum:
					/*Nothing in P vector*/
					break;
				case RheologyBbarAbsGradientEnum:
					/*Nothing in P vector*/
					break;
				default:
					_error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
			}
		}
	}

	/*Transform coordinate system*/
	element->TransformLoadVectorCoord(pe,XYEnum);

	/*Clean up and return*/
	xDelete<int>(responses);
	xDelete<IssmDouble>(xyz_list_top);
	xDelete<IssmDouble>(basis);
	delete gauss;
	return pe;

}/*}}}*/
ElementVector* AdjointHorizAnalysis::CreatePVectorSSA(Element* element){/*{{{*/

	/*Intermediaries*/
	int      domaintype;
	Element* basalelement;

	/*Get basal element*/
	element->FindParam(&domaintype,DomainTypeEnum);
	switch(domaintype){
		case Domain2DhorizontalEnum:
			basalelement = element;
			break;
		case Domain3DEnum:
			if(!element->IsOnBase()) return NULL;
			basalelement = element->SpawnBasalElement();
			break;
		default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
	}

	/*Intermediaries */
	int         num_responses,i;
	IssmDouble  Jdet,obs_velocity_mag,velocity_mag;
	IssmDouble  vx,vy,vxobs,vyobs,dux,duy,weight;
	IssmDouble scalex,scaley,scale,S;
	int        *responses = NULL;
	IssmDouble *xyz_list  = NULL;

	/*Fetch number of nodes and dof for this finite element*/
	int numnodes = basalelement->GetNumberOfNodes();

	/*Initialize Element vector and vectors*/
	ElementVector* pe    = basalelement->NewElementVector(SSAApproximationEnum);
	IssmDouble*    basis = xNew<IssmDouble>(numnodes);

	/*Retrieve all inputs and parameters*/
	basalelement->GetVerticesCoordinates(&xyz_list);
	basalelement->FindParam(&num_responses,InversionNumCostFunctionsEnum);
	basalelement->FindParam(&responses,NULL,InversionCostFunctionsEnum);
	Input* weights_input = basalelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
	Input* vx_input      = basalelement->GetInput(VxEnum);                                 _assert_(vx_input);
	Input* vy_input      = basalelement->GetInput(VyEnum);                                 _assert_(vy_input);
	Input* vxobs_input   = basalelement->GetInput(InversionVxObsEnum);                     _assert_(vxobs_input);
	Input* vyobs_input   = basalelement->GetInput(InversionVyObsEnum);                     _assert_(vyobs_input);
	IssmDouble epsvel  = 2.220446049250313e-16;
	IssmDouble meanvel = 3.170979198376458e-05; /*1000 m/yr*/

	/*Get Surface if required by one response*/
	for(int resp=0;resp<num_responses;resp++){
		if(responses[resp]==SurfaceAverageVelMisfitEnum){
			basalelement->GetInputValue(&S,SurfaceAreaEnum); break;
		}
	}

	/* Start  looping on the number of gaussian points: */
	Gauss* gauss=basalelement->NewGauss(4);
	for(int ig=gauss->begin();ig<gauss->end();ig++){
		gauss->GaussPoint(ig);

		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
		basalelement->NodalFunctions(basis, gauss);

		vx_input->GetInputValue(&vx,gauss);
		vy_input->GetInputValue(&vy,gauss);
		vxobs_input->GetInputValue(&vxobs,gauss);
		vyobs_input->GetInputValue(&vyobs,gauss);

		/*Loop over all requested responses*/
		for(int resp=0;resp<num_responses;resp++){
			weights_input->GetInputValue(&weight,gauss,responses[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<numnodes;i++){
						dux=vxobs-vx;
						duy=vyobs-vy;
						pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
						pe->values[i*2+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<numnodes;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*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
						pe->values[i*2+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<numnodes;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*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
						pe->values[i*2+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<numnodes;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*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
						pe->values[i*2+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<numnodes;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*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
						pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
					}
					break;
				case DragCoefficientAbsGradientEnum:
					/*Nothing in P vector*/
					break;
				case ThicknessAbsGradientEnum:
					/*Nothing in P vector*/
					break;
				case ThicknessAlongGradientEnum:
					/*Nothing in P vector*/
					break;
				case ThicknessAcrossGradientEnum:
					/*Nothing in P vector*/
					break;
				case RheologyBbarAbsGradientEnum:
					/*Nothing in P vector*/
					break;
				default:
					_error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
			}
		}
	}

	/*Transform coordinate system*/
	basalelement->TransformLoadVectorCoord(pe,XYEnum);

	/*Clean up and return*/
	xDelete<int>(responses);
	xDelete<IssmDouble>(xyz_list);
	xDelete<IssmDouble>(basis);
	if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
	delete gauss;
	return pe;
}/*}}}*/
void AdjointHorizAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
	   _error_("not implemented yet");
}/*}}}*/
void AdjointHorizAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
	/*The gradient of the cost function is calculated in 2 parts.
	 *
	 * dJ    \partial J   \partial lambda^T(KU-F)
	 * --  = ---------- + ------------------------
	 * dk    \partial k   \parial k                  
	 *
	 * */

	/*If on water, grad = 0: */
	if(!element->IsIceInElement()) return;

	/*Get Approximation*/
	int approximation;
	element->GetInputValue(&approximation,ApproximationEnum);

	/*Get list of cost functions*/
	int *responses = NULL;
	int num_responses,resp;
	element->FindParam(&num_responses,InversionNumCostFunctionsEnum);
	element->FindParam(&responses,NULL,InversionCostFunctionsEnum);

	/*Check that control_type is supported*/
	if(control_type!=MaterialsRheologyBbarEnum && 
		control_type!=FrictionCoefficientEnum   && 
		control_type!=DamageDbarEnum){
		_error_("Control "<<EnumToStringx(control_type)<<" not supported");
	}

	/*Deal with first part (partial derivative a J with respect to k)*/
	for(resp=0;resp<num_responses;resp++) switch(responses[resp]){
		case SurfaceAbsVelMisfitEnum:     /*Nothing, \partial J/\partial k = 0*/ break;
		case SurfaceRelVelMisfitEnum:     /*Nothing, \partial J/\partial k = 0*/ break;
		case SurfaceLogVelMisfitEnum:     /*Nothing, \partial J/\partial k = 0*/ break;
		case SurfaceLogVxVyMisfitEnum:    /*Nothing, \partial J/\partial k = 0*/ break;
		case SurfaceAverageVelMisfitEnum: /*Nothing, \partial J/\partial k = 0*/ break;
		case DragCoefficientAbsGradientEnum: GradientJDragGradient(element,gradient,control_index); break;
		case RheologyBbarAbsGradientEnum:    GradientJBGradient(element,gradient,control_index);    break;
		default: _error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
	}

	/*Deal with second term*/
	switch(control_type){
		case FrictionCoefficientEnum:
			switch(approximation){
				case SSAApproximationEnum: GradientJDragSSA(element,gradient,control_index); break;
				case HOApproximationEnum:  GradientJDragHO( element,gradient,control_index); break;
				case FSApproximationEnum:  GradientJDragFS( element,gradient,control_index); break;
				case NoneApproximationEnum: /*Gradient is 0*/                    break;
				default: _error_("approximation " << EnumToStringx(approximation) << " not supported yet");
			}
			break;
		case MaterialsRheologyBbarEnum:
			switch(approximation){
				case SSAApproximationEnum: GradientJBbarSSA(element,gradient,control_index); break;
				case HOApproximationEnum:  GradientJBbarHO( element,gradient,control_index); break;
				case FSApproximationEnum:  GradientJBbarFS( element,gradient,control_index); break;
				case NoneApproximationEnum: /*Gradient is 0*/                    break;
				default: _error_("approximation " << EnumToStringx(approximation) << " not supported yet");
			}
			break;
		case DamageDbarEnum:
			switch(approximation){
				case SSAApproximationEnum: GradientJDSSA(element,gradient,control_index); break;
				case NoneApproximationEnum: /*Gradient is 0*/                 break;
				default: _error_("approximation " << EnumToStringx(approximation) << " not supported yet");
			}
			break;
		default: _error_("control type not supported yet: " << EnumToStringx(control_type));
	}

	/*Clean up and return*/
	xDelete<int>(responses);
                         
}/*}}}*/
void AdjointHorizAnalysis::GradientJDragGradient(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/

	/*return if floating (gradient is 0)*/
	if(element->IsFloating()) return;

	/*Intermediaries*/
	int      domaintype,dim;
	Element* basalelement;

	/*Get basal element*/
	element->FindParam(&domaintype,DomainTypeEnum);
	switch(domaintype){
		case Domain2DhorizontalEnum:
			basalelement = element;
			dim          = 2;
			break;
		case Domain2DverticalEnum:
			if(!element->IsOnBase()) return;
			basalelement = element->SpawnBasalElement();
			dim          = 1;
			break;
		case Domain3DEnum:
			if(!element->IsOnBase()) return;
			basalelement = element->SpawnBasalElement();
			dim          = 2;
			break;
		default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
	}

	/*Intermediaries*/
	IssmDouble Jdet,weight;
	IssmDouble dk[3]; 
	IssmDouble *xyz_list= NULL;

	/*Fetch number of vertices for this finite element*/
	int numvertices = basalelement->GetNumberOfVertices();

	/*Initialize some vectors*/
	IssmDouble* dbasis        = xNew<IssmDouble>(2*numvertices);
	IssmDouble* ge            = xNewZeroInit<IssmDouble>(numvertices);
	int*        vertexpidlist = xNew<int>(numvertices);

	/*Retrieve all inputs we will be needing: */
	basalelement->GetVerticesCoordinates(&xyz_list);
	basalelement->GradientIndexing(&vertexpidlist[0],control_index);
	Input* dragcoefficient_input = basalelement->GetInput(FrictionCoefficientEnum);                _assert_(dragcoefficient_input);
	Input* weights_input         = basalelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);

	/* Start  looping on the number of gaussian points: */
	Gauss* gauss=basalelement->NewGauss(2);
	for(int ig=gauss->begin();ig<gauss->end();ig++){
		gauss->GaussPoint(ig);

		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
		basalelement->NodalFunctionsP1Derivatives(dbasis,xyz_list,gauss);
		weights_input->GetInputValue(&weight,gauss,DragCoefficientAbsGradientEnum);

		/*Build alpha_complement_list: */
		dragcoefficient_input->GetInputDerivativeValue(&dk[0],xyz_list,gauss);

		/*Build gradient vector (actually -dJ/ddrag): */
		for(int i=0;i<numvertices;i++){
			if(dim==2){
				ge[i]+=-weight*Jdet*gauss->weight*(dbasis[0*numvertices+i]*dk[0]+dbasis[1*numvertices+i]*dk[1]);
			}
			else{
				ge[i]+=-weight*Jdet*gauss->weight*dbasis[0*numvertices+i]*dk[0];
			}
			_assert_(!xIsNan<IssmDouble>(ge[i]));
		}
	}
	gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);

	/*Clean up and return*/
	xDelete<IssmDouble>(xyz_list);
	xDelete<IssmDouble>(dbasis);
	xDelete<IssmDouble>(ge);
	xDelete<int>(vertexpidlist);
	delete gauss;
	if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};

}/*}}}*/
void AdjointHorizAnalysis::GradientJBGradient(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/

	/*Intermediaries*/
	int      domaintype,dim;
	Element* basalelement;

	/*Get basal element*/
	element->FindParam(&domaintype,DomainTypeEnum);
	switch(domaintype){
		case Domain2DhorizontalEnum:
			basalelement = element;
			dim          = 2;
			break;
		case Domain2DverticalEnum:
			if(!element->IsOnBase()) return;
			basalelement = element->SpawnBasalElement();
			dim          = 1;
			break;
		case Domain3DEnum:
			if(!element->IsOnBase()) return;
			basalelement = element->SpawnBasalElement();
			dim          = 2;
			break;
		default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
	}

	/*Intermediaries*/
	IssmDouble Jdet,weight;
	IssmDouble dk[3]; 
	IssmDouble *xyz_list= NULL;

	/*Fetch number of vertices for this finite element*/
	int numvertices = basalelement->GetNumberOfVertices();

	/*Initialize some vectors*/
	IssmDouble* dbasis        = xNew<IssmDouble>(2*numvertices);
	IssmDouble* ge            = xNewZeroInit<IssmDouble>(numvertices);
	int*        vertexpidlist = xNew<int>(numvertices);

	/*Retrieve all inputs we will be needing: */
	basalelement->GetVerticesCoordinates(&xyz_list);
	basalelement->GradientIndexing(&vertexpidlist[0],control_index);
	Input* rheologyb_input = basalelement->GetInput(MaterialsRheologyBbarEnum);              _assert_(rheologyb_input);
	Input* weights_input   = basalelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);

	/* Start  looping on the number of gaussian points: */
	Gauss* gauss=basalelement->NewGauss(2);
	for(int ig=gauss->begin();ig<gauss->end();ig++){
		gauss->GaussPoint(ig);

		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
		basalelement->NodalFunctionsP1Derivatives(dbasis,xyz_list,gauss);
		weights_input->GetInputValue(&weight,gauss,RheologyBbarAbsGradientEnum);

		/*Build alpha_complement_list: */
		rheologyb_input->GetInputDerivativeValue(&dk[0],xyz_list,gauss);

		/*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
		for(int i=0;i<numvertices;i++){
			if(dim==2){
				ge[i]+=-weight*Jdet*gauss->weight*(dbasis[0*numvertices+i]*dk[0]+dbasis[1*numvertices+i]*dk[1]);
			}
			else{
				ge[i]+=-weight*Jdet*gauss->weight*dbasis[0*numvertices+i]*dk[0];
			}
			_assert_(!xIsNan<IssmDouble>(ge[i]));
		}
	}
	gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);

	/*Clean up and return*/
	xDelete<IssmDouble>(xyz_list);
	xDelete<IssmDouble>(dbasis);
	xDelete<IssmDouble>(ge);
	xDelete<int>(vertexpidlist);
	delete gauss;
	if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
}/*}}}*/
void AdjointHorizAnalysis::GradientJDragSSA(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/

	/*return if floating (gradient is 0)*/
	if(element->IsFloating()) return;

	/*Intermediaries*/
	int      domaintype,dim;
	Element* basalelement;

	/*Get basal element*/
	element->FindParam(&domaintype,DomainTypeEnum);
	switch(domaintype){
		case Domain2DhorizontalEnum:
			basalelement = element;
			dim          = 2;
			break;
		case Domain2DverticalEnum:
			if(!element->IsOnBase()) return;
			basalelement = element->SpawnBasalElement();
			dim          = 1;
			break;
		case Domain3DEnum:
			if(!element->IsOnBase()) return;
			basalelement = element->SpawnBasalElement();
			dim          = 2;
			break;
		default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
	}

	/*Intermediaries*/
	IssmDouble Jdet,weight;
	IssmDouble drag,dalpha2dk;
	IssmDouble vx,vy,lambda,mu;
	IssmDouble *xyz_list= NULL;

	/*Fetch number of vertices for this finite element*/
	int numvertices = basalelement->GetNumberOfVertices();

	/*Initialize some vectors*/
	IssmDouble* basis         = xNew<IssmDouble>(numvertices);
	IssmDouble* ge            = xNewZeroInit<IssmDouble>(numvertices);
	int*        vertexpidlist = xNew<int>(numvertices);

	/*Build friction element, needed later: */
	Friction* friction=new Friction(basalelement,dim);

	/*Retrieve all inputs we will be needing: */
	basalelement->GetVerticesCoordinates(&xyz_list);
	basalelement->GradientIndexing(&vertexpidlist[0],control_index);
	Input* vx_input        = basalelement->GetInput(VxEnum);                   _assert_(vx_input);
	Input* vy_input        = basalelement->GetInput(VyEnum);                   _assert_(vy_input);
	Input* adjointx_input  = basalelement->GetInput(AdjointxEnum);             _assert_(adjointx_input);
	Input* adjointy_input  = basalelement->GetInput(AdjointyEnum);             _assert_(adjointy_input);
	Input* dragcoeff_input = basalelement->GetInput(FrictionCoefficientEnum);  _assert_(dragcoeff_input);

	/* Start  looping on the number of gaussian points: */
	Gauss* gauss=basalelement->NewGauss(4);
	for(int ig=gauss->begin();ig<gauss->end();ig++){
		gauss->GaussPoint(ig);

		adjointx_input->GetInputValue(&lambda, gauss);
		adjointy_input->GetInputValue(&mu, gauss);
		vx_input->GetInputValue(&vx,gauss);
		vy_input->GetInputValue(&vy,gauss);
		dragcoeff_input->GetInputValue(&drag, gauss);

		friction->GetAlphaComplement(&dalpha2dk,gauss);

		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
		basalelement->NodalFunctionsP1(basis,gauss);

		/*Build gradient vector (actually -dJ/dD): */
		for(int i=0;i<numvertices;i++){
			ge[i]+=-2.*drag*dalpha2dk*((lambda*vx+mu*vy))*Jdet*gauss->weight*basis[i];
			_assert_(!xIsNan<IssmDouble>(ge[i]));
		}
	}
	gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);

	/*Clean up and return*/
	xDelete<IssmDouble>(xyz_list);
	xDelete<IssmDouble>(basis);
	xDelete<IssmDouble>(ge);
	xDelete<int>(vertexpidlist);
	delete gauss;
	delete friction;
	if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
}/*}}}*/
void AdjointHorizAnalysis::GradientJDragHO(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/

	/*return if floating or not on bed (gradient is 0)*/
	if(element->IsFloating()) return;
	if(!element->IsOnBase()) return;

	/*Intermediaries*/
	int        dim=3;
	IssmDouble Jdet,weight;
	IssmDouble drag,dalpha2dk;
	IssmDouble vx,vy,lambda,mu;
	IssmDouble *xyz_list_base= NULL;

	/*Fetch number of vertices for this finite element*/
	int numvertices = element->GetNumberOfVertices();

	/*Initialize some vectors*/
	IssmDouble* basis         = xNew<IssmDouble>(numvertices);
	IssmDouble* ge            = xNewZeroInit<IssmDouble>(numvertices);
	int*        vertexpidlist = xNew<int>(numvertices);

	/*Build friction element, needed later: */
	Friction* friction=new Friction(element,dim);

	/*Retrieve all inputs we will be needing: */
	element->GetVerticesCoordinatesBase(&xyz_list_base);
	element->GradientIndexing(&vertexpidlist[0],control_index);
	Input* vx_input        = element->GetInput(VxEnum);                   _assert_(vx_input);
	Input* vy_input        = element->GetInput(VyEnum);                   _assert_(vy_input);
	Input* adjointx_input  = element->GetInput(AdjointxEnum);             _assert_(adjointx_input);
	Input* adjointy_input  = element->GetInput(AdjointyEnum);             _assert_(adjointy_input);
	Input* dragcoeff_input = element->GetInput(FrictionCoefficientEnum);  _assert_(dragcoeff_input);

	/* Start  looping on the number of gaussian points: */
	Gauss* gauss=element->NewGaussBase(4);
	for(int ig=gauss->begin();ig<gauss->end();ig++){
		gauss->GaussPoint(ig);

		adjointx_input->GetInputValue(&lambda, gauss);
		adjointy_input->GetInputValue(&mu, gauss);
		vx_input->GetInputValue(&vx,gauss);
		vy_input->GetInputValue(&vy,gauss);
		dragcoeff_input->GetInputValue(&drag, gauss);

		friction->GetAlphaComplement(&dalpha2dk,gauss);

		element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
		element->NodalFunctionsP1(basis,gauss);

		/*Build gradient vector (actually -dJ/dD): */
		for(int i=0;i<numvertices;i++){
			ge[i]+=-2.*drag*dalpha2dk*((lambda*vx+mu*vy))*Jdet*gauss->weight*basis[i];
			_assert_(!xIsNan<IssmDouble>(ge[i]));
		}
	}
	gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);

	/*Clean up and return*/
	xDelete<IssmDouble>(xyz_list_base);
	xDelete<IssmDouble>(basis);
	xDelete<IssmDouble>(ge);
	xDelete<int>(vertexpidlist);
	delete gauss;
	delete friction;
}/*}}}*/
void AdjointHorizAnalysis::GradientJDragFS(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/

	/*return if floating or not on bed (gradient is 0)*/
	if(element->IsFloating()) return;
	if(!element->IsOnBase()) return;

	/*Intermediaries*/
	int        dim=3;
	IssmDouble Jdet,weight;
	IssmDouble drag,dalpha2dk,normal[3];
	IssmDouble vx,vy,vz,lambda,mu,xi;
	IssmDouble *xyz_list_base= NULL;

	/*Fetch number of vertices for this finite element*/
	int numvertices = element->GetNumberOfVertices();

	/*Initialize some vectors*/
	IssmDouble* basis         = xNew<IssmDouble>(numvertices);
	IssmDouble* ge            = xNewZeroInit<IssmDouble>(numvertices);
	int*        vertexpidlist = xNew<int>(numvertices);

	/*Build friction element, needed later: */
	Friction* friction=new Friction(element,dim);

	/*Retrieve all inputs we will be needing: */
	element->GetVerticesCoordinatesBase(&xyz_list_base);
	element->GradientIndexing(&vertexpidlist[0],control_index);
	Input* vx_input        = element->GetInput(VxEnum);                   _assert_(vx_input);
	Input* vy_input        = element->GetInput(VyEnum);                   _assert_(vy_input);
	Input* vz_input        = element->GetInput(VzEnum);                   _assert_(vy_input);
	Input* adjointx_input  = element->GetInput(AdjointxEnum);             _assert_(adjointx_input);
	Input* adjointy_input  = element->GetInput(AdjointyEnum);             _assert_(adjointy_input);
	Input* adjointz_input  = element->GetInput(AdjointzEnum);             _assert_(adjointz_input);
	Input* dragcoeff_input = element->GetInput(FrictionCoefficientEnum);  _assert_(dragcoeff_input);

	/* Start  looping on the number of gaussian points: */
	Gauss* gauss=element->NewGaussBase(4);
	for(int ig=gauss->begin();ig<gauss->end();ig++){
		gauss->GaussPoint(ig);

		adjointx_input->GetInputValue(&lambda, gauss);
		adjointy_input->GetInputValue(&mu, gauss);
		adjointz_input->GetInputValue(&xi    ,gauss);
		vx_input->GetInputValue(&vx,gauss);
		vy_input->GetInputValue(&vy,gauss);
		vz_input->GetInputValue(&vz,gauss);
		dragcoeff_input->GetInputValue(&drag, gauss);

		friction->GetAlphaComplement(&dalpha2dk,gauss);
		element->NormalBase(&normal[0],xyz_list_base);

		element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
		element->NodalFunctionsP1(basis,gauss);

		/*Build gradient vector (actually -dJ/dk): */
		for(int i=0;i<numvertices;i++){
			ge[i]+=(
						-lambda*(2*drag*dalpha2dk*(vx - vz*normal[0]*normal[2]))
						-mu    *(2*drag*dalpha2dk*(vy - vz*normal[1]*normal[2]))
						-xi    *(2*drag*dalpha2dk*(-vx*normal[0]*normal[2]-vy*normal[1]*normal[2]))
						)*Jdet*gauss->weight*basis[i];
			_assert_(!xIsNan<IssmDouble>(ge[i]));
		}
	}
	gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);

	/*Clean up and return*/
	xDelete<IssmDouble>(xyz_list_base);
	xDelete<IssmDouble>(basis);
	xDelete<IssmDouble>(ge);
	xDelete<int>(vertexpidlist);
	delete gauss;
	delete friction;
}/*}}}*/
void AdjointHorizAnalysis::GradientJBbarSSA(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/

	/*Intermediaries*/
	int      domaintype,dim;
	Element* basalelement;

	/*Get basal element*/
	element->FindParam(&domaintype,DomainTypeEnum);
	switch(domaintype){
		case Domain2DhorizontalEnum:
			basalelement = element;
			dim          = 2;
			break;
		case Domain2DverticalEnum:
			if(!element->IsOnBase()) return;
			basalelement = element->SpawnBasalElement();
			dim          = 1;
			break;
		case Domain3DEnum:
			if(!element->IsOnBase()) return;
			basalelement = element->SpawnBasalElement();
			dim          = 2;
			break;
		default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
	}

	/*Intermediaries*/
	IssmDouble Jdet,weight;
	IssmDouble thickness,dmudB;
	IssmDouble dvx[3],dvy[3],dadjx[3],dadjy[3]; 
	IssmDouble *xyz_list= NULL;

	/*Fetch number of vertices for this finite element*/
	int numvertices = basalelement->GetNumberOfVertices();

	/*Initialize some vectors*/
	IssmDouble* basis         = xNew<IssmDouble>(numvertices);
	IssmDouble* ge            = xNewZeroInit<IssmDouble>(numvertices);
	int*        vertexpidlist = xNew<int>(numvertices);

	/*Retrieve all inputs we will be needing: */
	basalelement->GetVerticesCoordinates(&xyz_list);
	basalelement->GradientIndexing(&vertexpidlist[0],control_index);
	Input* thickness_input = basalelement->GetInput(ThicknessEnum);             _assert_(thickness_input);
	Input* vx_input        = basalelement->GetInput(VxEnum);                    _assert_(vx_input);
	Input* vy_input        = basalelement->GetInput(VyEnum);                    _assert_(vy_input);
	Input* adjointx_input  = basalelement->GetInput(AdjointxEnum);              _assert_(adjointx_input);
	Input* adjointy_input  = basalelement->GetInput(AdjointyEnum);              _assert_(adjointy_input);
	Input* rheologyb_input = basalelement->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input);

	/* Start  looping on the number of gaussian points: */
	Gauss* gauss=basalelement->NewGauss(4);
	for(int ig=gauss->begin();ig<gauss->end();ig++){
		gauss->GaussPoint(ig);

		thickness_input->GetInputValue(&thickness,gauss);
		vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
		vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
		adjointx_input->GetInputDerivativeValue(&dadjx[0],xyz_list,gauss);
		adjointy_input->GetInputDerivativeValue(&dadjy[0],xyz_list,gauss);

		basalelement->dViscositydBSSA(&dmudB,dim,xyz_list,gauss,vx_input,vy_input);

		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
		basalelement->NodalFunctionsP1(basis,gauss);

		/*Build gradient vector (actually -dJ/dB): */
		for(int i=0;i<numvertices;i++){
			ge[i]+=-dmudB*thickness*(
						(2*dvx[0]+dvy[1])*2*dadjx[0]+(dvx[1]+dvy[0])*(dadjx[1]+dadjy[0])+(2*dvy[1]+dvx[0])*2*dadjy[1]
						)*Jdet*gauss->weight*basis[i];
			_assert_(!xIsNan<IssmDouble>(ge[i]));
		}
	}
	gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);

	/*Clean up and return*/
	xDelete<IssmDouble>(xyz_list);
	xDelete<IssmDouble>(basis);
	xDelete<IssmDouble>(ge);
	xDelete<int>(vertexpidlist);
	delete gauss;
	if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
}/*}}}*/
void AdjointHorizAnalysis::GradientJBbarHO(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/

	/*WARNING: We use SSA as an estimate for now*/
	this->GradientJBbarSSA(element,gradient,control_index);
}/*}}}*/
void AdjointHorizAnalysis::GradientJBbarFS(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
	/*WARNING: We use SSA as an estimate for now*/
	this->GradientJBbarSSA(element,gradient,control_index);
}/*}}}*/
void AdjointHorizAnalysis::GradientJDSSA(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/

	/*Intermediaries*/
	int      domaintype,dim;
	Element* basalelement;

	/*Get basal element*/
	element->FindParam(&domaintype,DomainTypeEnum);
	switch(domaintype){
		case Domain2DhorizontalEnum:
			basalelement = element;
			dim          = 2;
			break;
		case Domain2DverticalEnum:
			if(!element->IsOnBase()) return;
			basalelement = element->SpawnBasalElement();
			dim          = 1;
			break;
		case Domain3DEnum:
			if(!element->IsOnBase()) return;
			basalelement = element->SpawnBasalElement();
			dim          = 2;
			break;
		default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
	}

	/*Intermediaries*/
	IssmDouble Jdet,weight;
	IssmDouble thickness,dmudD;
	IssmDouble dvx[3],dvy[3],dadjx[3],dadjy[3]; 
	IssmDouble *xyz_list= NULL;

	/*Fetch number of vertices for this finite element*/
	int numvertices = basalelement->GetNumberOfVertices();

	/*Initialize some vectors*/
	IssmDouble* basis         = xNew<IssmDouble>(numvertices);
	IssmDouble* ge            = xNewZeroInit<IssmDouble>(numvertices);
	int*        vertexpidlist = xNew<int>(numvertices);

	/*Retrieve all inputs we will be needing: */
	basalelement->GetVerticesCoordinates(&xyz_list);
	basalelement->GradientIndexing(&vertexpidlist[0],control_index);
	Input* thickness_input = basalelement->GetInput(ThicknessEnum);             _assert_(thickness_input);
	Input* vx_input        = basalelement->GetInput(VxEnum);                    _assert_(vx_input);
	Input* vy_input        = basalelement->GetInput(VyEnum);                    _assert_(vy_input);
	Input* adjointx_input  = basalelement->GetInput(AdjointxEnum);              _assert_(adjointx_input);
	Input* adjointy_input  = basalelement->GetInput(AdjointyEnum);              _assert_(adjointy_input);
	Input* rheologyb_input = basalelement->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input);

	/* Start  looping on the number of gaussian points: */
	Gauss* gauss=basalelement->NewGauss(4);
	for(int ig=gauss->begin();ig<gauss->end();ig++){
		gauss->GaussPoint(ig);

		thickness_input->GetInputValue(&thickness,gauss);
		vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
		vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
		adjointx_input->GetInputDerivativeValue(&dadjx[0],xyz_list,gauss);
		adjointy_input->GetInputDerivativeValue(&dadjy[0],xyz_list,gauss);

		basalelement->dViscositydDSSA(&dmudD,dim,xyz_list,gauss,vx_input,vy_input);

		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
		basalelement->NodalFunctionsP1(basis,gauss);

		for(int i=0;i<numvertices;i++){
			ge[i]+=-dmudD*thickness*(
						(2*dvx[0]+dvy[1])*2*dadjx[0]+(dvx[1]+dvy[0])*(dadjx[1]+dadjy[0])+(2*dvy[1]+dvx[0])*2*dadjy[1]
						)*Jdet*gauss->weight*basis[i];
			_assert_(!xIsNan<IssmDouble>(ge[i]));
		}
	}
	gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);

	/*Clean up and return*/
	xDelete<IssmDouble>(xyz_list);
	xDelete<IssmDouble>(basis);
	xDelete<IssmDouble>(ge);
	xDelete<int>(vertexpidlist);
	delete gauss;
	if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
}/*}}}*/
void AdjointHorizAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
	int approximation;
	element->GetInputValue(&approximation,ApproximationEnum);
	if(approximation==FSApproximationEnum || approximation==NoneApproximationEnum){
		InputUpdateFromSolutionFS(solution,element);
	}
	else{
		InputUpdateFromSolutionHoriz(solution,element);
	}
}/*}}}*/
void AdjointHorizAnalysis::InputUpdateFromSolutionHoriz(IssmDouble* solution,Element* element){/*{{{*/
	int  i;
	int* doflist=NULL;

	/*Fetch number of nodes and dof for this finite element*/
	int numnodes = element->GetNumberOfNodes();
	int numdof   = numnodes*2;

	/*Fetch dof list and allocate solution vectors*/
	element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
	IssmDouble* values  = xNew<IssmDouble>(numdof);
	IssmDouble* lambdax = xNew<IssmDouble>(numnodes);
	IssmDouble* lambday = xNew<IssmDouble>(numnodes);

	/*Use the dof list to index into the solution vector: */
	for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];

	/*Transform solution in Cartesian Space*/
	element->TransformSolutionCoord(&values[0],XYEnum);

	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
	for(i=0;i<numnodes;i++){
		lambdax[i]=values[i*NDOF2+0];
		lambday[i]=values[i*NDOF2+1];

		/*Check solution*/
		if(xIsNan<IssmDouble>(lambdax[i])) _error_("NaN found in solution vector");
		if(xIsNan<IssmDouble>(lambday[i])) _error_("NaN found in solution vector");
	}

	/*Add vx and vy as inputs to the tria element: */
	element->AddInput(AdjointxEnum,lambdax,element->GetElementType());
	element->AddInput(AdjointyEnum,lambday,element->GetElementType());

	/*Free ressources:*/
	xDelete<IssmDouble>(values);
	xDelete<IssmDouble>(lambdax);
	xDelete<IssmDouble>(lambday);
	xDelete<int>(doflist);
}/*}}}*/
void AdjointHorizAnalysis::InputUpdateFromSolutionFS(IssmDouble* solution,Element* element){/*{{{*/
	int          i,dim;
	int*         vdoflist=NULL;
	int*         pdoflist=NULL;
	IssmDouble   FSreconditioning;

	element->FindParam(&dim,DomainDimensionEnum);

	/*Fetch number of nodes and dof for this finite element*/
	int vnumnodes = element->NumberofNodesVelocity();
	int pnumnodes = element->NumberofNodesPressure();
	int vnumdof   = vnumnodes*dim;
	int pnumdof   = pnumnodes*1;

	/*Initialize values*/
	IssmDouble* values  = xNew<IssmDouble>(vnumdof+pnumdof);
	IssmDouble* lambdax = xNew<IssmDouble>(vnumnodes);
	IssmDouble* lambday = xNew<IssmDouble>(vnumnodes);
	IssmDouble* lambdaz = xNew<IssmDouble>(vnumnodes);
	IssmDouble* lambdap = xNew<IssmDouble>(pnumnodes);

	int* cs_list = xNew<int>(vnumnodes+pnumnodes);
	if(dim==2) for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
	else       for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum;
	for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;

	/*Get dof list: */
	element->GetDofListVelocity(&vdoflist,GsetEnum);
	element->GetDofListPressure(&pdoflist,GsetEnum);

	/*Use the dof list to index into the solution vector: */
	for(i=0;i<vnumdof;i++) values[i]        =solution[vdoflist[i]];
	for(i=0;i<pnumdof;i++) values[vnumdof+i]=solution[pdoflist[i]];

	/*Transform solution in Cartesian Space*/
	element->TransformSolutionCoord(values,cs_list);

	/*fill in all arrays: */
	for(i=0;i<vnumnodes;i++){
		lambdax[i] = values[i*NDOF3+0]; if(xIsNan<IssmDouble>(lambdax[i])) _error_("NaN found in solution vector");
		lambday[i] = values[i*NDOF3+1]; if(xIsNan<IssmDouble>(lambday[i])) _error_("NaN found in solution vector");
		lambdaz[i] = values[i*NDOF3+2]; if(xIsNan<IssmDouble>(lambdaz[i])) _error_("NaN found in solution vector");
	}
	for(i=0;i<pnumnodes;i++){
		lambdap[i] = values[vnumdof+i]; if(xIsNan<IssmDouble>(lambdap[i])) _error_("NaN found in solution vector");
	}

	/*Recondition pressure and compute vel: */
	element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
	for(i=0;i<pnumnodes;i++) lambdap[i]=lambdap[i]*FSreconditioning;

	/*Add vx and vy as inputs to the tria element: */
	element->AddInput(AdjointxEnum,lambdax,element->VelocityInterpolation());
	element->AddInput(AdjointyEnum,lambday,element->VelocityInterpolation());
	element->AddInput(AdjointzEnum,lambdaz,element->VelocityInterpolation());
	element->AddInput(AdjointpEnum,lambdap,element->PressureInterpolation());

	/*Free ressources:*/
	xDelete<int>(vdoflist);
	xDelete<int>(pdoflist);
	xDelete<int>(cs_list);
	xDelete<IssmDouble>(lambdap);
	xDelete<IssmDouble>(lambdaz);
	xDelete<IssmDouble>(lambday);
	xDelete<IssmDouble>(lambdax);
	xDelete<IssmDouble>(values);
}/*}}}*/
void AdjointHorizAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
	/*Default, do nothing*/
	return;
}/*}}}*/
