/*!\file FemModel.cpp
 * \brief: implementation of the FemModel object
 */

#ifdef HAVE_CONFIG_H
	#include <config.h>
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif

#include <stdio.h>
#include "../cores/cores.h"
#include "../shared/io/io.h"
#include "./classes.h"
#include "./modules/modules.h"
#include "../shared/Enum/Enum.h"

#include "../analyses/analyses.h"

/*module includes: {{{*/
#include "../modules/ModelProcessorx/ModelProcessorx.h"
#include "../modules/VerticesDofx/VerticesDofx.h"
#include "../modules/SpcNodesx/SpcNodesx.h"
#include "../modules/ConfigureObjectsx/ConfigureObjectsx.h"
#include "../modules/ParseToolkitsOptionsx/ParseToolkitsOptionsx.h"
#include "../modules/GetVectorFromInputsx/GetVectorFromInputsx.h"
#include "../modules/InputUpdateFromVectorx/InputUpdateFromVectorx.h"
#include "../modules/NodesDofx/NodesDofx.h"
#include "../modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.h"
#include "../modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.h"
#include "../modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.h"
#include "../modules/SurfaceLogVxVyMisfitx/SurfaceLogVxVyMisfitx.h"
#include "../modules/SurfaceAverageVelMisfitx/SurfaceAverageVelMisfitx.h"
#include "../modules/ThicknessAbsMisfitx/ThicknessAbsMisfitx.h"
#include "../modules/ThicknessAlongGradientx/ThicknessAlongGradientx.h"
#include "../modules/ThicknessAcrossGradientx/ThicknessAcrossGradientx.h"
#include "../modules/RheologyBbarAbsGradientx/RheologyBbarAbsGradientx.h"
#include "../modules/DragCoefficientAbsGradientx/DragCoefficientAbsGradientx.h"
#include "../modules/NodalValuex/NodalValuex.h"
#include "../modules/GetVectorFromInputsx/GetVectorFromInputsx.h"
#include "../modules/AverageOntoPartitionx/AverageOntoPartitionx.h"
/*}}}*/

/*Object constructors and destructor*/
FemModel::FemModel(int argc,char** argv,ISSM_MPI_Comm incomm){/*{{{*/

	/*configuration: */
	int* analyses=NULL;
	int  numanalyses;
	int  solution_type;
	int  ierr;

	/*File names*/
	char *lockfilename   = NULL;
	char *binfilename    = NULL;
	char *outbinfilename = NULL;
	char *petscfilename  = NULL;
	char *rootpath       = NULL;

	/*First things first, store the communicator, and set it as a global variable: */
	this->comm=incomm;
	this->SetStaticComm();

	/*Now, initialize PETSC: */
	#ifdef _HAVE_PETSC_
	PETSC_COMM_WORLD=this->comm;
	ierr=PetscInitialize(&argc,&argv,(char*)0,"");  if(ierr) _error_("Could not initialize Petsc");
	#endif

	/*Start profiler: */
	this->profiler=new Profiler();
	profiler->Tag(Start);

	/*From command line arguments, retrieve different filenames needed to create the FemModel: */
	ProcessArguments(&solution_type,&binfilename,&outbinfilename,&petscfilename,&lockfilename,&rootpath,argc,argv);

	/*out of solution_type, figure out types of analyses needed in the femmodel: */
	AnalysisConfiguration(&analyses,&numanalyses,solution_type);

	/*Create femmodel from input files: */
	profiler->Tag(StartInit);
	this->InitFromFiles(rootpath,binfilename,outbinfilename,petscfilename,lockfilename,solution_type,analyses,numanalyses);
	profiler->Tag(FinishInit);

	/*Free resources */
	xDelete<int>(analyses);
	xDelete<char>(lockfilename);
	xDelete<char>(binfilename);
	xDelete<char>(outbinfilename);
	xDelete<char>(petscfilename);
	xDelete<char>(rootpath);

}
/*}}}*/
FemModel::FemModel(char* rootpath, char* inputfilename, char* outputfilename, char* toolkitsfilename, char* lockfilename, const int in_solution_type,const int* analyses,const int nummodels){/*{{{*/

	/*Call InitFromFiles. This constructor is just a wrapper: */
	this->InitFromFiles(rootpath, inputfilename, outputfilename, toolkitsfilename, lockfilename, in_solution_type,analyses,nummodels);

}
/*}}}*/
FemModel::~FemModel(){/*{{{*/

	/*Intermediary*/
	FILE *output_fid;
	char *outbinfilename = NULL;
	char *lockfilename   = NULL;

	this->parameters->FindParam(&outbinfilename,OutputFileNameEnum);
	this->parameters->FindParam(&lockfilename,LockFileNameEnum);

	/*Delete all the datasets: */
	xDelete<int>(analysis_type_list);
	xDelete<char>(outbinfilename);
	xDelete<char>(lockfilename);
	delete elements;
	delete nodes;
	delete vertices;
	delete constraints;
	delete loads;
	delete materials;
	delete parameters;
	delete results;

	/*Now delete: */
	delete profiler;

}
/*}}}*/

/*Object management*/
void FemModel::Echo(void){/*{{{*/

	_printf_("FemModel echo: \n");
	_printf_("   number of fem models: " << nummodels << "\n");
	_printf_("   analysis_type_list: \n");
	for(int i=0;i<nummodels;i++)_printf_("     " << i << ": " << EnumToStringx(analysis_type_list[i]) << "\n");
	_printf_("   current analysis_type: \n");
	_printf_("     " << analysis_counter << ": " << EnumToStringx(analysis_type_list[analysis_counter]) << "\n");

}
/*}}}*/
void FemModel::InitFromFiles(char* rootpath, char* inputfilename, char* outputfilename, char* toolkitsfilename, char* lockfilename, const int in_solution_type,const int* analyses,const int nummodels){/*{{{*/

	/*intermediary*/
	int         i;
	int         analysis_type;
	FILE       *IOMODEL = NULL;
	FILE       *toolkitsoptionsfid = NULL;
	FILE       *output_fid = NULL;
	int         my_rank;

	/*recover my_rank:*/
	my_rank=IssmComm::GetRank();

	/*Open input file on cpu 0: */
	if(my_rank==0) IOMODEL = pfopen0(inputfilename ,"rb");

	/*Open toolkits file: */
	toolkitsoptionsfid=pfopen(toolkitsfilename,"r");

	/*Initialize internal data: */
	this->nummodels        = nummodels;
	this->solution_type    = in_solution_type;
	this->analysis_counter = nummodels-1;   //point to last analysis_type carried out.
	this->results          = new Results(); //not initialized by CreateDataSets

	/*Dynamically allocate whatever is a list of length nummodels: */
	analysis_type_list=xNew<int>(nummodels);

	/*Initialize: */
	for(i=0;i<nummodels;i++)analysis_type_list[i]=analyses[i];

	/*create datasets for all analyses*/
	ModelProcessorx(&this->elements,&this->nodes,&this->vertices,&this->materials,&this->constraints,&this->loads,&this->parameters,IOMODEL,toolkitsoptionsfid,rootpath,this->solution_type,nummodels,analyses);

	/*do the post-processing of the datasets to get an FemModel that can actually run analyses: */
	for(i=0;i<nummodels;i++){

		if(VerboseMProcessor()) _printf0_("   Processing finite element model of analysis " << EnumToStringx(analysis_type_list[i]) << ":\n");
		analysis_type=analysis_type_list[i];
		this->SetCurrentConfiguration(analysis_type);

		if(i==0){
			if(VerboseMProcessor()) _printf0_("      creating vertex PIDs\n");
			VerticesDofx(vertices,parameters); //only call once, we only have one set of vertices
		}

		if(VerboseMProcessor()) _printf0_("      resolving node constraints\n");
		SpcNodesx(nodes,constraints,parameters,analysis_type); 

		if(VerboseMProcessor()) _printf0_("      creating nodal degrees of freedom\n");
		NodesDofx(nodes,parameters,analysis_type);

		if(VerboseMProcessor()) _printf0_("      configuring element and loads\n");
		ConfigureObjectsx(elements, loads, nodes, vertices, materials,parameters);
	}

	/*Close input file and toolkits file descriptors: */
	if(my_rank==0) pfclose(IOMODEL,inputfilename);
	pfclose(toolkitsoptionsfid,toolkitsfilename);

	/*Open output file once for all and add output file name and file descriptor to parameters*/
	output_fid=pfopen(outputfilename,"wb");
	this->parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename));
	this->parameters->SetParam(output_fid,OutputFilePointerEnum);

	/*Save lock file name for later: */
	this->parameters->AddObject(new StringParam(LockFileNameEnum,lockfilename));

	}
/*}}}*/
void FemModel::CleanUp(void){/*{{{*/

	/*Intermediary*/
	FILE *output_fid;
	char *outbinfilename = NULL;
	char *lockfilename   = NULL;
	bool  waitonlock     = false;

	/*Close output file: */
	this->parameters->FindParam(&output_fid,OutputFilePointerEnum);
	this->parameters->FindParam(&outbinfilename,OutputFileNameEnum);
	pfclose(output_fid,outbinfilename);

	/*Write lock file if requested: */
	this->parameters->FindParam(&waitonlock,SettingsWaitonlockEnum);
	this->parameters->FindParam(&lockfilename,LockFileNameEnum);
	if(waitonlock){
		_printf0_("write lock file:\n");
		WriteLockFile(lockfilename);
	}

	/*Before we delete the profiler, report statistics for this run: */
	profiler->Tag(Finish);  //final tagging
	_printf0_("\n");
	_printf0_("   "<<setw(40)<<left<<"FemModel initialization elapsed time:"<<profiler->DeltaTime(StartInit,FinishInit) << "\n");
	_printf0_("   "<<setw(40)<<left<<"Core solution elapsed time:"<<profiler->DeltaTime(StartCore,FinishCore) << "\n");
	_printf0_("\n");
	_printf0_("   Total elapsed time:"
				<<profiler->DeltaTimeModHour(Start,Finish)<<" hrs "
				<<profiler->DeltaTimeModMin(Start,Finish)<<" min "
				<<profiler->DeltaTimeModSec(Start,Finish)<<" sec"
				);
	_printf0_("\n");

	/*Finalize PETSC for this model: */
	#ifdef _HAVE_PETSC_
	_printf0_("closing PETSc\n");
	PetscFinalize();
	#endif


	/*Clean up*/
	xDelete<char>(outbinfilename);
	xDelete<char>(lockfilename);

}
/*}}}*/
void FemModel::SetStaticComm(void){/*{{{*/

	/*This routine sets the global communicator variable hidden inside the IssmComm 
	 *class: */
	IssmComm::SetComm(this->comm);

}
/*}}}*/
void FemModel::SetCurrentConfiguration(int configuration_type,int analysis_type){/*{{{*/

	/*Use configuration_type to setup the analysis counter, the configurations of objects etc ... but use 
	 * analysis_type to drive the element numerics. This allows for use of 1 configuration_type for several 
	 * analyses. For example: do a SurfaceSlopeX, SurfaceSlopeY, BedSlopeX and BedSlopeY analysis using the 
	 * Slope configuration.*/
	int found=-1;
	for(int i=0;i<nummodels;i++){
	
		if (analysis_type_list[i]==configuration_type){
			found=i;
			break;
		}
	}
	if(found!=-1) analysis_counter=found;
	else _error_("Could not find alias for analysis_type " << EnumToStringx(configuration_type) << " in list of FemModel analyses");

	/*Now, plug analysis_counter and analysis_type inside the parameters: */
	this->parameters->SetParam(analysis_counter,AnalysisCounterEnum);
	this->parameters->SetParam(analysis_type,AnalysisTypeEnum);
	this->parameters->SetParam(configuration_type,ConfigurationTypeEnum);

	/*configure elements, loads and nodes, for this new analysis: */
	this->elements->SetCurrentConfiguration(elements,loads, nodes,vertices, materials,parameters);
	this->loads->SetCurrentConfiguration(elements, loads, nodes,vertices, materials,parameters);

	/*take care of toolkits options, that depend on this analysis type (present only after model processor)*/
	if(this->parameters->Exist(ToolkitsOptionsStringsEnum)){
		ToolkitsOptionsFromAnalysis(this->parameters,analysis_type);
		if(VerboseSolver()) _printf0_("      toolkits Options set for analysis type: " << EnumToStringx(analysis_type) << "\n");
	}

}
/*}}}*/
void FemModel::SetCurrentConfiguration(int configuration_type){/*{{{*/
	this->SetCurrentConfiguration(configuration_type,configuration_type);
}
/*}}}*/
void FemModel::Solve(void){/*{{{*/

	/*profiling: */
	bool profiling = false;
	IssmDouble solution_time;
	IssmDouble solution_flops;
	IssmDouble solution_memory;

	/*solution: */
	int solution_type;
	void (*solutioncore)(FemModel*)=NULL; //core solution function pointer

	_printf0_("call computational core:\n");

	/*Retrieve solution_type from parameters: */
	parameters->FindParam(&solution_type,SolutionTypeEnum);

	/*Figure out which solution core we are going to run with the current solution type: */
	WrapperCorePointerFromSolutionEnum(&solutioncore,this->parameters,solution_type);

	/*run solutoin core: */
	profiler->Tag(StartCore);   
	solutioncore(this); 
	profiler->Tag(FinishCore);

	/*run AD core if needed: */
	profiler->Tag(StartAdCore); 
	ad_core(this);      
	profiler->Tag(FinishAdCore);

	/*some profiling results for the core: */
	parameters->FindParam(&profiling,DebugProfilingEnum);
	if(profiling){

		solution_time=profiler->DeltaTime(StartCore,FinishCore);
		solution_flops=profiler->DeltaFlops(StartCore,FinishCore);
		solution_memory=profiler->Memory(FinishCore);

		_printf0_("Core solution elapsed time    : " << solution_time   << " Seconds\n");
		_printf0_("Core solution number of flops : " << solution_flops  << " Flops\n");
		_printf0_("Core solution memory used     : " << solution_memory << " Bytes\n");

		/*Add to results: */
		results->AddObject(new GenericExternalResult<IssmDouble>(results->Size()+1, ProfilingSolutionTimeEnum, solution_time, 1, 0));
		results->AddObject(new GenericExternalResult<IssmDouble>(results->Size()+1, ProfilingCurrentMemEnum, solution_memory, 1, 0));
		results->AddObject(new GenericExternalResult<IssmDouble>(results->Size()+1, ProfilingCurrentFlopsEnum, solution_flops, 1, 0));

		#ifdef _HAVE_ADOLC_
		solution_time=profiler->DeltaTime(StartAdCore,FinishAdCore);
		solution_flops=profiler->DeltaFlops(StartAdCore,FinishAdCore);
		solution_memory=profiler->Memory(FinishAdCore);

		_printf0_("AD Solution elapsed time    : " << solution_time   << " Seconds\n");
		_printf0_("AD Solution number of flops : " << solution_flops  << " Flops\n");
		_printf0_("AD Solution memory used     : " << solution_memory << " Bytes\n");

		#endif

			}
}
/*}}}*/
FemModel* FemModel::copy(void){/*{{{*/

	FemModel* output=NULL;
	int       i;
	int       analysis_type;

	output=new FemModel(*this); //Use default copy constructor.

	output->comm = this->comm;
	output->nummodels = this->nummodels;
	output->solution_type = this->solution_type;
	output->analysis_counter = this->analysis_counter;

	/*Now, deep copy arrays: */
	output->analysis_type_list=xNew<int>(nummodels);
	xMemCpy<int>(output->analysis_type_list,this->analysis_type_list,this->nummodels);

	output->profiler=static_cast<Profiler*>(this->profiler->copy());

	output->loads=static_cast<Loads*>(this->loads->Copy());
	output->materials=static_cast<Materials*>(this->materials->Copy());
	output->parameters=static_cast<Parameters*>(this->parameters->Copy());
	output->constraints=static_cast<Constraints*>(this->constraints->Copy());
	output->results=static_cast<Results*>(this->results->Copy());

	output->nodes=static_cast<Nodes*>(this->nodes->Copy());
	output->vertices=static_cast<Vertices*>(this->vertices->Copy());
	output->elements=static_cast<Elements*>(this->elements->Copy());

	/*reset hooks for elements, loads and nodes: */
	output->elements->ResetHooks();
	output->loads->ResetHooks();
	output->materials->ResetHooks();

	/*do the post-processing of the datasets to get an FemModel that can actually run analyses: */
	for(i=0;i<nummodels;i++){
		analysis_type=output->analysis_type_list[i];
		output->SetCurrentConfiguration(analysis_type);
		if(i==0) VerticesDofx(output->vertices,output->parameters); //only call once, we only have one set of vertices
		SpcNodesx(output->nodes,output->constraints,output->parameters,analysis_type);
		NodesDofx(output->nodes,output->parameters,analysis_type);
		ConfigureObjectsx(output->elements,output->loads,output->nodes,output->vertices,output->materials,output->parameters);
	}

	return output;
}
/*}}}*/

/*Modules:*/
int  FemModel::UpdateVertexPositionsx(void){ /*{{{*/

	IssmDouble         *surface = NULL;
	IssmDouble         *bed     = NULL;

	/*get vertex vectors for bed and thickness: */
	GetVectorFromInputsx(&surface  ,this, SurfaceEnum,VertexEnum);
	GetVectorFromInputsx(&bed      ,this, BaseEnum,    VertexEnum);

	/*Allocate vector*/
	Vector<IssmDouble> *vx=new Vector<IssmDouble>(vertices->NumberOfVertices());
	Vector<IssmDouble> *vy=new Vector<IssmDouble>(vertices->NumberOfVertices());
	Vector<IssmDouble> *vz=new Vector<IssmDouble>(vertices->NumberOfVertices());

	/*Update verices new geometry: */
	for(int i=0;i<vertices->Size();i++){
		Vertex* vertex=(Vertex*)vertices->GetObjectByOffset(i);
		vertex->UpdatePosition(vx,vy,vz,parameters,surface,bed);
	}

	/*Assemble mesh velocity*/
	vx->Assemble();
	vy->Assemble();
	vz->Assemble();

	/*Update element inputs*/
	InputUpdateFromVectorx(this,vx,VxMeshEnum,VertexPIdEnum);
	InputUpdateFromVectorx(this,vy,VyMeshEnum,VertexPIdEnum);
	InputUpdateFromVectorx(this,vz,VzMeshEnum,VertexPIdEnum);

	/*Free ressources:*/
	delete vx;
	delete vy;
	delete vz;
	xDelete<IssmDouble>(bed);
	xDelete<IssmDouble>(surface);
	return 1;
}
/*}}}*/
void FemModel::UpdateConstraintsx(void){ /*{{{*/

	IssmDouble time;
	int        analysis_type;

	/*retrieve parameters: */
	parameters->FindParam(&analysis_type,AnalysisTypeEnum);
	parameters->FindParam(&time,TimeEnum);

	/*start module: */
	if(VerboseModule()) _printf0_("   Updating constraints and active domain of analysis " << EnumToStringx(analysis_type)  << " for time: " << time << "\n");

	Analysis* analysis= EnumToAnalysis(analysis_type);
	analysis->UpdateConstraints(this);
	delete analysis;
	
	/*Second, constraints might be time dependent: */
	SpcNodesx(nodes,constraints,parameters,analysis_type); 

	/*Now, update degrees of freedoms: */
	NodesDofx(nodes,parameters,analysis_type);

}
/*}}}*/
void FemModel::Responsex(IssmDouble* responses,const char* response_descriptor){/*{{{*/

	int response_descriptor_enum=StringToEnumx(response_descriptor);
	this->Responsex(responses, response_descriptor_enum);

}
/*}}}*/
void FemModel::Responsex(IssmDouble* responses,int response_descriptor_enum){/*{{{*/

	switch (response_descriptor_enum){

		case DivergenceEnum:               this->Divergencex(responses); break;
		case IceVolumeEnum:                this->IceVolumex(responses); break;
		case IceVolumeAboveFloatationEnum: this->IceVolumeAboveFloatationx(responses); break;
		case MinVelEnum:                   this->MinVelx(responses); break;
		case MaxVelEnum:                   this->MaxVelx(responses); break;
		case MinVxEnum:                    this->MinVxx(responses); break;
		case MaxVxEnum:                    this->MaxVxx(responses); break;
		case MaxAbsVxEnum:                 this->MaxAbsVxx(responses); break;
		case MinVyEnum:                    this->MinVyx(responses); break;
		case MaxVyEnum:                    this->MaxVyx(responses); break;
		case MaxAbsVyEnum:                 this->MaxAbsVyx(responses); break;
		case MinVzEnum:                    this->MinVzx(responses); break;
		case MaxVzEnum:                    this->MaxVzx(responses); break;
		case MaxAbsVzEnum:                 this->MaxAbsVzx(responses); break;
		case MassFluxEnum:                 this->MassFluxx(responses); break;
		case SurfaceAbsVelMisfitEnum:      SurfaceAbsVelMisfitx(responses, elements,nodes, vertices, loads, materials,parameters); break;
		case SurfaceRelVelMisfitEnum:      SurfaceRelVelMisfitx(responses, elements,nodes, vertices, loads, materials,parameters); break;
		case SurfaceLogVelMisfitEnum:      SurfaceLogVelMisfitx(responses, elements,nodes, vertices, loads, materials,parameters); break;
		case SurfaceLogVxVyMisfitEnum:     SurfaceLogVxVyMisfitx(responses, elements,nodes, vertices, loads, materials,parameters); break;
		case SurfaceAverageVelMisfitEnum:  SurfaceAverageVelMisfitx(responses,this); break;
		case ThicknessAbsMisfitEnum:       ThicknessAbsMisfitx(responses, elements,nodes, vertices, loads, materials, parameters); break;
		case ThicknessAbsGradientEnum:     this->ThicknessAbsGradientx(responses); break;
		case ThicknessAlongGradientEnum:   ThicknessAlongGradientx(responses, elements,nodes, vertices, loads, materials, parameters); break;
		case ThicknessAcrossGradientEnum:  ThicknessAcrossGradientx(responses, elements,nodes, vertices, loads, materials, parameters); break;
		case RheologyBbarAbsGradientEnum:  RheologyBbarAbsGradientx(responses, elements,nodes, vertices, loads, materials, parameters); break;
		case DragCoefficientAbsGradientEnum:DragCoefficientAbsGradientx(responses, elements,nodes, vertices, loads, materials, parameters); break;
		case BalancethicknessMisfitEnum:   BalancethicknessMisfitx(responses); break;
		case Balancethickness2MisfitEnum:  Balancethickness2Misfitx(responses); break;
		case TotalSmbEnum:					  this->TotalSmbx(responses); break;
		case MaterialsRheologyBbarEnum:    this->ElementResponsex(responses,MaterialsRheologyBbarEnum); break;
		case VelEnum:                      this->ElementResponsex(responses,VelEnum); break;
		case FrictionCoefficientEnum:      NodalValuex(responses, FrictionCoefficientEnum,elements,nodes, vertices, loads, materials, parameters); break;
		default: _error_("response descriptor \"" << EnumToStringx(response_descriptor_enum) << "\" not supported yet!"); break; 
	}

}
/*}}}*/
void FemModel::RequestedOutputsx(Results **presults,char** requested_outputs, int numoutputs, bool save_results){/*{{{*/

	/*Intermediaries*/
	bool        isvec,results_on_nodes;
	int         step,output_enum;
	IssmDouble  time;
	IssmDouble  double_result;
	const char *output_string = NULL;

	/*recover results*/
	Results* results = *presults;
	if(!results) results = new Results();

	/*Get time and step*/
	parameters->FindParam(&step,StepEnum);
	parameters->FindParam(&time,TimeEnum);
	parameters->FindParam(&results_on_nodes,SettingsResultsOnNodesEnum);

	/*Go through all requested output*/
	for(int i=0;i<numoutputs;i++){

		output_string = requested_outputs[i];
		output_enum   = StringToEnumx(output_string,false);
		isvec         = false;

		/*If string is not an enum, it is defined in output definitions*/
		if(output_enum<0){
			double_result = OutputDefinitionsResponsex(this,output_string);
			if(save_results){
				results->AddResult(new GenericExternalResult<IssmPDouble>(results->Size()+1,output_string,reCast<IssmPDouble>(double_result),step,time));
				continue;
			}
		}
		else{
			switch(output_enum){

				/*Scalar output*/
				case DivergenceEnum:               this->Divergencex(&double_result);               break;
				case IceVolumeEnum:                this->IceVolumex(&double_result);                break;
				case IceVolumeAboveFloatationEnum: this->IceVolumeAboveFloatationx(&double_result); break;
				case MinVelEnum:                   this->MinVelx(&double_result);                   break;
				case MaxVelEnum:                   this->MaxVelx(&double_result);                   break;
				case MinVxEnum:                    this->MinVxx(&double_result);                    break;
				case MaxVxEnum:                    this->MaxVxx(&double_result);                    break;
				case MaxAbsVxEnum:                 this->MaxAbsVxx(&double_result);                 break;
				case MinVyEnum:                    this->MinVyx(&double_result);                    break;
				case MaxVyEnum:                    this->MaxVyx(&double_result);                    break;
				case MaxAbsVyEnum:                 this->MaxAbsVyx(&double_result);                 break;
				case MinVzEnum:                    this->MinVzx(&double_result);                    break;
				case MaxVzEnum:                    this->MaxVzx(&double_result);                    break;
				case MaxAbsVzEnum:                 this->MaxAbsVzx(&double_result);                 break;
				case MassFluxEnum:                 this->MassFluxx(&double_result);                 break;
				case TotalSmbEnum:                 this->TotalSmbx(&double_result);                 break;

			   /*Scalar control output*/
				case SurfaceAbsVelMisfitEnum:       SurfaceAbsVelMisfitx(&double_result,elements,nodes,vertices,loads,materials,parameters);        break;
				case SurfaceRelVelMisfitEnum:       SurfaceRelVelMisfitx(&double_result,elements,nodes,vertices,loads,materials,parameters);        break;
				case SurfaceLogVelMisfitEnum:       SurfaceLogVelMisfitx(&double_result,elements,nodes,vertices,loads,materials,parameters);        break;
				case SurfaceLogVxVyMisfitEnum:      SurfaceLogVxVyMisfitx(&double_result,elements,nodes,vertices,loads,materials,parameters);       break;
				case SurfaceAverageVelMisfitEnum:   SurfaceAverageVelMisfitx(&double_result,this);                                                  break;
				case ThicknessAbsMisfitEnum:        ThicknessAbsMisfitx(&double_result,elements,nodes,vertices,loads,materials,parameters);         break;
				case ThicknessAbsGradientEnum:      this->ThicknessAbsGradientx(&double_result);                                                    break;
				case ThicknessAlongGradientEnum:    ThicknessAlongGradientx(&double_result,elements,nodes,vertices,loads,materials,parameters);     break;
				case ThicknessAcrossGradientEnum:   ThicknessAcrossGradientx(&double_result,elements,nodes,vertices,loads,materials,parameters);    break;
				case RheologyBbarAbsGradientEnum:   RheologyBbarAbsGradientx(&double_result,elements,nodes,vertices,loads,materials,parameters);    break;
				case DragCoefficientAbsGradientEnum:DragCoefficientAbsGradientx(&double_result,elements,nodes,vertices,loads,materials,parameters); break;
				case BalancethicknessMisfitEnum:    BalancethicknessMisfitx(&double_result);                                                        break;
				case Balancethickness2MisfitEnum:  Balancethickness2Misfitx(&double_result); break;

			   /*Vector */
				default:

					/*Vector layout*/
					int interpolation,nodesperelement,size;
					int rank_interpolation=-1,rank_nodesperelement=-1;

					/*Get interpolation (and compute input if necessary)*/
					for(int j=0;j<elements->Size();j++){
						Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(j));
						element->ResultInterpolation(&rank_interpolation,&rank_nodesperelement,output_enum);
					}

					/*Broadcast for cpus that do not have any elements*/
					ISSM_MPI_Reduce(&rank_interpolation,&interpolation,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm());
					ISSM_MPI_Reduce(&rank_nodesperelement,&nodesperelement,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm());
					ISSM_MPI_Bcast(&interpolation,1,ISSM_MPI_INT,0,IssmComm::GetComm());
					ISSM_MPI_Bcast(&nodesperelement,1,ISSM_MPI_INT,0,IssmComm::GetComm());

					if(results_on_nodes){

						/*Allocate matrices*/
						int         nbe       = this->elements->NumberOfElements();
						IssmDouble* values    = xNewZeroInit<IssmDouble>(nbe*nodesperelement);
						IssmDouble* allvalues = xNew<IssmDouble>(nbe*nodesperelement);

						/*Fill-in matrix*/
						for(int j=0;j<elements->Size();j++){
							Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(j));
							element->ResultToPatch(values,nodesperelement,output_enum);
						}

						/*Gather from all cpus*/
						ISSM_MPI_Allreduce((void*)values,(void*)allvalues,nbe*nodesperelement,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
						xDelete<IssmDouble>(values);

						if(save_results)results->AddResult(new GenericExternalResult<IssmDouble*>(results->Size()+1,output_enum,allvalues,nbe,nodesperelement,step,time));
						xDelete<IssmDouble>(allvalues);

					}
					else{

						/*Allocate vector depending on interpolation*/
						switch(interpolation){
							case P0Enum: size = this->elements->NumberOfElements(); break;
							case P1Enum: size = this->vertices->NumberOfVertices(); break;
							default:     _error_("Interpolation "<<EnumToStringx(interpolation)<<" not supported yet");

						}
						Vector<IssmDouble> *vector_result = new Vector<IssmDouble>(size);

						/*Fill in vector*/
						for(int j=0;j<elements->Size();j++){
							Element* element=(Element*)elements->GetObjectByOffset(j);
							element->ResultToVector(vector_result,output_enum);
						}
						vector_result->Assemble();

						if(save_results)results->AddResult(new GenericExternalResult<Vector<IssmDouble>*>(results->Size()+1,output_enum,vector_result,step,time));
					}
					isvec = true;
					break;
			}
		}

		/*Add result to Results*/
		if(!isvec && save_results){
			results->AddResult(new GenericExternalResult<IssmPDouble>(results->Size()+1,output_string,reCast<IssmPDouble>(double_result),step,time));
		}
	}

	/*Assign pointer and clean up*/
	*presults = results;
}
/*}}}*/
void FemModel::RequestedOutputsx(Results **presults,int* requested_outputs, int numoutputs,bool save_results){/*{{{*/

	/*Convert list of enums to list of string*/
	char** enumlist = xNew<char*>(numoutputs);
	for(int i=0;i<numoutputs;i++){
		EnumToStringx(&enumlist[i],requested_outputs[i]);
	}

	/*Call main module*/
	this->RequestedOutputsx(presults,enumlist,numoutputs,save_results);

	/*clean up and return*/
	for(int i=0;i<numoutputs;i++) xDelete<char>(enumlist[i]);
	xDelete<char*>(enumlist);
	return;
}
/*}}}*/
void FemModel::RequestedDependentsx(void){/*{{{*/

	bool        isautodiff      = false;
	IssmDouble  output_value;

	int         num_dependents;
	IssmPDouble *dependents;
	DataSet*    dependent_objects=NULL;
	int my_rank=IssmComm::GetRank();

	/*AD mode on?: */
	parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum);

	if(isautodiff){
		#ifdef _HAVE_ADOLC_
		parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum);
		parameters->FindParam(&dependent_objects,AutodiffDependentObjectsEnum);
		if(num_dependents){
			dependents=xNew<IssmPDouble>(num_dependents);

			/*Go through our dependent variables, and compute the response:*/
			for(int i=0;i<dependent_objects->Size();i++){
				DependentObject* dep=(DependentObject*)dependent_objects->GetObjectByOffset(i);
				dep->Responsex(&output_value,this);
				if (my_rank==0) {
					output_value>>=dependents[i];
				}
			}
		}
		delete dependent_objects;
		if(num_dependents)xDelete<IssmPDouble>(dependents);
		#else
		_error_("Should not be requesting dependents when an AD library is not available!");
		#endif
	}
}
/*}}}*/
void FemModel::TimeAdaptx(IssmDouble* pdt){/*{{{*/

	int      i;

	/*output: */
	IssmDouble   dt;

	/*intermediary: */
	Element *element     = NULL;
	IssmDouble   min_dt      = 0;
	IssmDouble   node_min_dt = 0;

	/*Go through elements, and figure out the minimum of the time steps for each element (using CFL criterion): */
	element=(Element*)elements->GetObjectByOffset(0); min_dt=element->TimeAdapt();

	for (i=1;i<elements->Size();i++){
		element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
		dt=element->TimeAdapt();
		if(dt<min_dt)min_dt=dt;
	}

	/*Figure out minimum across the cluster: */
	ISSM_MPI_Reduce (&min_dt,&node_min_dt,1,ISSM_MPI_DOUBLE,ISSM_MPI_MIN,0,IssmComm::GetComm() );
	ISSM_MPI_Bcast(&node_min_dt,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
	min_dt=node_min_dt;

	/*Assign output pointers:*/
	*pdt=min_dt;
}
/*}}}*/
void FemModel::UpdateConstraintsExtrudeFromBasex(void){ /*{{{*/

	for(int i=0;i<elements->Size();i++){
		Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
		element->UpdateConstraintsExtrudeFromBase();
	}

}
/*}}}*/
void FemModel::UpdateConstraintsExtrudeFromTopx(void){ /*{{{*/

	for(int i=0;i<elements->Size();i++){
		Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
		element->UpdateConstraintsExtrudeFromTop();
	}

}
/*}}}*/
void FemModel::MassFluxx(IssmDouble* pmass_flux){/*{{{*/

	int          i,j;
	Element     *element       = NULL;
	int          element_id;
	bool         ispresent     = false;
	IssmDouble   mass_flux     = 0;
	IssmDouble   all_mass_flux = 0;
	int          counter;
	IssmDouble **array         = NULL;
	int          M;
	int         *mdims_array   = NULL;
	int         *ndims_array   = NULL;
	IssmDouble  *segments      = NULL;
	int          num_segments;

	/*First, figure out which segment to compute our mass flux on. Start with retrieving qmu_mass_flux_segments: */
	this->parameters->FindParam(&ispresent,MassFluxSegmentsPresentEnum);
	if(!ispresent)_error_("no mass flux segments available!");
	this->parameters->FindParam(&array,&M,&mdims_array,&ndims_array,MassFluxSegmentsEnum);

	/*Retrieve index of segments being used for MassFlux computation: */
	parameters->FindParam(&counter,IndexEnum);

	/*retrieve segments from array: */
	segments     = array[counter-1]; //matlab to "C" indexing
	num_segments = mdims_array[counter-1];

	/*Go through segments, and then elements, and figure out which elements belong to a segment. 
	 * When we find one, use the element to compute the mass flux on the segment: */
	for(i=0;i<num_segments;i++){
		element_id=reCast<int,IssmDouble>(*(segments+5*i+4));
		for(j=0;j<elements->Size();j++){
			element=(Element*)this->elements->GetObjectByOffset(j);
			if (element->Id()==element_id){
				/*We found the element which owns this segment, use it to compute the mass flux: */
				mass_flux+=element->MassFlux(segments+5*i+0);
				break;
			}
		}
	}

	ISSM_MPI_Allreduce ( (void*)&mass_flux,(void*)&all_mass_flux,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
	mass_flux=all_mass_flux;

	/*Free ressources:*/
	for(i=0;i<M;i++){
		IssmDouble* matrix=array[i];
		xDelete<IssmDouble>(matrix);
	}
	xDelete<int>(mdims_array);
	xDelete<int>(ndims_array);
	xDelete<IssmDouble*>(array);

	/*Assign output pointers: */
	*pmass_flux=mass_flux;

}/*}}}*/
void FemModel::MaxAbsVxx(IssmDouble* pmaxabsvx){/*{{{*/

	int i;
	IssmDouble maxabsvx;
	IssmDouble node_maxabsvx;
	IssmDouble element_maxabsvx;

	/*Go through elements, and request velocity: */
	maxabsvx=-INFINITY;
	for(i=0;i<this->elements->Size();i++){
		Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
		element_maxabsvx=element->inputs->MaxAbs(VxEnum);
		if(element_maxabsvx>maxabsvx) maxabsvx=element_maxabsvx;
	}

	/*Figure out maximum across the cluster: */
	ISSM_MPI_Reduce(&maxabsvx,&node_maxabsvx,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
	ISSM_MPI_Bcast(&node_maxabsvx,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());   
	maxabsvx=node_maxabsvx;

	/*Assign output pointers:*/
	*pmaxabsvx=maxabsvx;

}/*}}}*/
void FemModel::MaxAbsVyx(IssmDouble* pmaxabsvy){/*{{{*/

	int i;
	IssmDouble maxabsvy;
	IssmDouble node_maxabsvy;
	IssmDouble element_maxabsvy;

	/*Go through elements, and request velocity: */
	maxabsvy=-INFINITY;
	for(i=0;i<this->elements->Size();i++){
		Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
		element_maxabsvy=element->inputs->MaxAbs(VyEnum);
		if(element_maxabsvy>maxabsvy) maxabsvy=element_maxabsvy;
	}

	/*Figure out maximum across the cluster: */
	ISSM_MPI_Reduce(&maxabsvy,&node_maxabsvy,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
	ISSM_MPI_Bcast(&node_maxabsvy,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());   
	maxabsvy=node_maxabsvy;

	/*Assign output pointers:*/
	*pmaxabsvy=maxabsvy;

}/*}}}*/
void FemModel::MaxAbsVzx(IssmDouble* pmaxabsvz){/*{{{*/

	int i;
	IssmDouble maxabsvz;
	IssmDouble node_maxabsvz;
	IssmDouble element_maxabsvz;

	/*Go through elements, and request velocity: */
	maxabsvz=-INFINITY;
	for(i=0;i<this->elements->Size();i++){
		Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
		element_maxabsvz=element->inputs->MaxAbs(VzEnum);
		if(element_maxabsvz>maxabsvz) maxabsvz=element_maxabsvz;
	}

	/*Figure out maximum across the cluster: */
	ISSM_MPI_Reduce(&maxabsvz,&node_maxabsvz,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
	ISSM_MPI_Bcast(&node_maxabsvz,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());   
	maxabsvz=node_maxabsvz;

	/*Assign output pointers:*/
	*pmaxabsvz=maxabsvz;

}/*}}}*/
void FemModel::MaxVelx(IssmDouble* pmaxvel){/*{{{*/

	int i;
	IssmDouble maxvel;
	IssmDouble node_maxvel;
	IssmDouble element_maxvel;

	/*Go through elements, and request velocity: */
	maxvel=-INFINITY;
	for(i=0;i<this->elements->Size();i++){
		Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
		element_maxvel = element->inputs->Max(VelEnum);
		if(element_maxvel>maxvel) maxvel=element_maxvel;
	}

	/*Figure out maximum across the cluster: */
	ISSM_MPI_Reduce(&maxvel,&node_maxvel,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
	ISSM_MPI_Bcast(&node_maxvel,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());   
	maxvel=node_maxvel;

	/*Assign output pointers:*/
	*pmaxvel=maxvel;

}/*}}}*/
void FemModel::MaxVxx(IssmDouble* pmaxvx){/*{{{*/

	int i;
	IssmDouble maxvx;
	IssmDouble node_maxvx;
	IssmDouble element_maxvx;

	/*Go through elements, and request velocity: */
	maxvx=-INFINITY;
	for(i=0;i<this->elements->Size();i++){
		Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
		element_maxvx = element->inputs->Max(VxEnum);
		if(element_maxvx>maxvx) maxvx=element_maxvx;
	}

	/*Figure out maximum across the cluster: */
	ISSM_MPI_Reduce(&maxvx,&node_maxvx,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
	ISSM_MPI_Bcast(&node_maxvx,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());   
	maxvx=node_maxvx;

	/*Assign output pointers:*/
	*pmaxvx=maxvx;

}/*}}}*/
void FemModel::MaxVyx(IssmDouble* pmaxvy){/*{{{*/

	int i;
	IssmDouble maxvy;
	IssmDouble node_maxvy;
	IssmDouble element_maxvy;

	/*Go through elements, and request velocity: */
	maxvy=-INFINITY;
	for(i=0;i<this->elements->Size();i++){
		Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
		element_maxvy = element->inputs->Max(VyEnum);
		if(element_maxvy>maxvy) maxvy=element_maxvy;
	}

	/*Figure out maximum across the cluster: */
	ISSM_MPI_Reduce(&maxvy,&node_maxvy,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
	ISSM_MPI_Bcast(&node_maxvy,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());   
	maxvy=node_maxvy;

	/*Assign output pointers:*/
	*pmaxvy=maxvy;

}/*}}}*/
void FemModel::MaxVzx(IssmDouble* pmaxvz){/*{{{*/

	int i;
	IssmDouble maxvz;
	IssmDouble node_maxvz;
	IssmDouble element_maxvz;

	/*Go through elements, and request velocity: */
	maxvz=-INFINITY;
	for(i=0;i<this->elements->Size();i++){
		Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
		element_maxvz = element->inputs->Max(VzEnum);
		if(element_maxvz>maxvz) maxvz=element_maxvz;
	}

	/*Figure out maximum across the cluster: */
	ISSM_MPI_Reduce(&maxvz,&node_maxvz,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
	ISSM_MPI_Bcast(&node_maxvz,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());   
	maxvz=node_maxvz;

	/*Assign output pointers:*/
	*pmaxvz=maxvz;

}/*}}}*/
void FemModel::MinVelx(IssmDouble* pminvel){/*{{{*/

	int i;
	IssmDouble minvel;
	IssmDouble node_minvel;
	IssmDouble element_minvel;

	/*Go through elements, and request velocity: */
	minvel=INFINITY;
	for(i=0;i<this->elements->Size();i++){
		Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
		element_minvel = element->inputs->Min(VelEnum);
		if(element_minvel<minvel) minvel=element_minvel;
	}

	/*Figure out minimum across the cluster: */
	ISSM_MPI_Reduce(&minvel,&node_minvel,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
	ISSM_MPI_Bcast(&node_minvel,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());   
	minvel=node_minvel;

	/*Assign output pointers:*/
	*pminvel=minvel;

}/*}}}*/
void FemModel::MinVxx(IssmDouble* pminvx){/*{{{*/

	int i;
	IssmDouble minvx;
	IssmDouble node_minvx;
	IssmDouble element_minvx;

	/*Go through elements, and request velocity: */
	minvx=INFINITY;
	for(i=0;i<this->elements->Size();i++){
		Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
		element_minvx = element->inputs->Min(VxEnum);
		if(element_minvx<minvx) minvx=element_minvx;
	}

	/*Figure out minimum across the cluster: */
	ISSM_MPI_Reduce(&minvx,&node_minvx,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
	ISSM_MPI_Bcast(&node_minvx,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());   
	minvx=node_minvx;

	/*Assign output pointers:*/
	*pminvx=minvx;

}/*}}}*/
void FemModel::MinVyx(IssmDouble* pminvy){/*{{{*/

	int i;
	IssmDouble minvy;
	IssmDouble node_minvy;
	IssmDouble element_minvy;

	/*Go through elements, and request velocity: */
	minvy=INFINITY;
	for(i=0;i<this->elements->Size();i++){
		Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
		element_minvy = element->inputs->Min(VyEnum);
		if(element_minvy<minvy) minvy=element_minvy;
	}

	/*Figure out minimum across the cluster: */
	ISSM_MPI_Reduce(&minvy,&node_minvy,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
	ISSM_MPI_Bcast(&node_minvy,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());   
	minvy=node_minvy;

	/*Assign output pointers:*/
	*pminvy=minvy;

}/*}}}*/
void FemModel::MinVzx(IssmDouble* pminvz){/*{{{*/

	int i;
	IssmDouble minvz;
	IssmDouble node_minvz;
	IssmDouble element_minvz;

	/*Go through elements, and request velocity: */
	minvz=INFINITY;
	for(i=0;i<this->elements->Size();i++){
		Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
		element_minvz = element->inputs->Min(VzEnum);
		if(element_minvz<minvz) minvz=element_minvz;
	}

	/*Figure out minimum across the cluster: */
	ISSM_MPI_Reduce(&minvz,&node_minvz,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
	ISSM_MPI_Bcast(&node_minvz,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());   
	minvz=node_minvz;

	/*Assign output pointers:*/
	*pminvz=minvz;

}/*}}}*/
void FemModel::TotalSmbx(IssmDouble* pSmb){/*{{{*/

	IssmDouble local_smb = 0;
	IssmDouble total_smb;

	for(int i=0;i<this->elements->Size();i++){
		Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
		local_smb+=element->TotalSmb();
	}
	ISSM_MPI_Reduce(&local_smb,&total_smb,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
	ISSM_MPI_Bcast(&total_smb,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());

	/*Assign output pointers: */
	*pSmb=total_smb;

}/*}}}*/
void FemModel::Divergencex(IssmDouble* pdiv){/*{{{*/

	IssmDouble local_divergence=0;
	IssmDouble total_divergence;

	for(int i=0;i<this->elements->Size();i++){
		Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
		local_divergence+=element->Divergence();
	}
	ISSM_MPI_Reduce(&local_divergence,&total_divergence,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
	ISSM_MPI_Bcast(&total_divergence,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());

	/*Assign output pointers: */
	*pdiv=total_divergence;

}/*}}}*/
void FemModel::IceVolumex(IssmDouble* pV){/*{{{*/

	IssmDouble local_ice_volume = 0;
	IssmDouble total_ice_volume;

	for(int i=0;i<this->elements->Size();i++){
		Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
		local_ice_volume+=element->IceVolume();
	}
	ISSM_MPI_Reduce(&local_ice_volume,&total_ice_volume,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
	ISSM_MPI_Bcast(&total_ice_volume,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());

	/*Assign output pointers: */
	*pV=total_ice_volume;

}/*}}}*/
void FemModel::IceVolumeAboveFloatationx(IssmDouble* pV){/*{{{*/

	IssmDouble local_ice_volume_af = 0;
	IssmDouble total_ice_volume_af;

	for(int i=0;i<this->elements->Size();i++){
		Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
		local_ice_volume_af+=element->IceVolumeAboveFloatation();
	}
	ISSM_MPI_Reduce(&local_ice_volume_af,&total_ice_volume_af,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
	ISSM_MPI_Bcast(&total_ice_volume_af,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());

	/*Assign output pointers: */
	*pV=total_ice_volume_af;

}/*}}}*/
void FemModel::ElementResponsex(IssmDouble* presponse,int response_enum){/*{{{*/

	int found=0;
	int sumfound=0;
	int cpu_found=-1;
	int index;
	IssmDouble response;
	Element* element=NULL;

	/*retrieve element we are interested in: */
	this->parameters->FindParam(&index,IndexEnum);
	int my_rank=IssmComm::GetRank();

	/*now, go through our elements, and retrieve the one with this id: index: */
	for(int i=0;i<this->elements->Size();i++){
		element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
		if (element->Id()==index){
			found=1;
			cpu_found=my_rank;
			break;
		}
	}

	/*Broadcast whether we found the element: */
	ISSM_MPI_Allreduce ( &found,&sumfound,1,ISSM_MPI_INT,ISSM_MPI_SUM,IssmComm::GetComm());
	if(!sumfound)_error_("could not find material with id" << index << " to compute ElementResponse");

	/*Ok, we found the element, compute responseocity: */
	if(my_rank==cpu_found){
		element->ElementResponse(&response,response_enum);
	}

	/*Broadcast and plug into response: */
	ISSM_MPI_Allreduce ( &cpu_found,&cpu_found,1,ISSM_MPI_INT,ISSM_MPI_MAX,IssmComm::GetComm());
	ISSM_MPI_Bcast(&response,1,ISSM_MPI_DOUBLE,cpu_found,IssmComm::GetComm()); 

	/*Assign output pointers: */
	*presponse=response;

}/*}}}*/
void FemModel::BalancethicknessMisfitx(IssmDouble* presponse){/*{{{*/

	/*output: */
	IssmDouble J=0.;
	IssmDouble J_sum;

	IssmDouble  weight,vx,vy,H,dvx[2],dvy[2],dH[2];
	IssmDouble  temp,Jdet,dhdt,groundedice_melting,surface_mass_balance;
	IssmDouble* xyz_list = NULL;
	IssmDouble  dp[3];

	/*Compute Misfit: */
	for(int i=0;i<elements->Size();i++){
		Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));

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

		/* Get node coordinates*/
		element->GetVerticesCoordinates(&xyz_list);
		Input* weights_input                   = element->GetInput(InversionCostFunctionsCoefficientsEnum);   _assert_(weights_input);
		Input* thickness_input                 = element->GetInput(ThicknessEnum); _assert_(thickness_input);
		Input* vx_input                        = element->GetInput(VxEnum);                                  _assert_(vx_input);
		Input* vy_input                        = element->GetInput(VyEnum);                                  _assert_(vy_input);
		Input* surface_mass_balance_input      = element->GetInput(SurfaceforcingsMassBalanceEnum);          _assert_(surface_mass_balance_input);
		Input* groundedice_melting_input       = element->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(groundedice_melting_input);
		Input* dhdt_input                      = element->GetInput(BalancethicknessThickeningRateEnum);      _assert_(dhdt_input);

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

			gauss->GaussPoint(ig);

			/* Get Jacobian determinant: */
			element->JacobianDeterminant(&Jdet,xyz_list,gauss);

			/*Get all parameters at gaussian point*/
			weights_input->GetInputValue(&weight,gauss,BalancethicknessMisfitEnum);
			thickness_input->GetInputValue(&H, gauss);
			thickness_input->GetInputDerivativeValue(&dH[0],xyz_list,gauss);
			surface_mass_balance_input->GetInputValue(&surface_mass_balance,gauss);
			groundedice_melting_input->GetInputValue(&groundedice_melting,gauss);
			dhdt_input->GetInputValue(&dhdt,gauss);
			vx_input->GetInputValue(&vx,gauss);
			vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
			vy_input->GetInputValue(&vy,gauss);
			vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);

			/*Balance thickness soft constraint J = 1/2 (div(Hv)-a)^2*/
			temp  = vx*dH[0]+vy*dH[1]+H*(dvx[0]+dvy[1]) - (surface_mass_balance-groundedice_melting-dhdt);
			J    +=weight*1/2*temp*temp*Jdet*gauss->weight;
		}

		/*clean up and Return: */
		xDelete<IssmDouble>(xyz_list);
		delete gauss;
	}

	/*Sum all J from all cpus of the cluster:*/
	ISSM_MPI_Reduce (&J,&J_sum,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
	ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
	J=J_sum;

	/*Assign output pointers: */
	*presponse=J;

}/*}}}*/
void FemModel::Balancethickness2Misfitx(IssmDouble* presponse){/*{{{*/

	/*output: */
	IssmDouble J=0.;
	IssmDouble J_sum;

	IssmDouble  weight,thicknessobs,thickness,potential;
	IssmDouble  vx,vy,vxbar,vybar,vxobs,vyobs,vxbarobs,vybarobs,nux,nuy;
	IssmDouble  Jdet;
	IssmDouble* xyz_list = NULL;

	/*Compute Misfit: */
	for(int i=0;i<elements->Size();i++){
		Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));

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

		/* Get node coordinates*/
		element->GetVerticesCoordinates(&xyz_list);
		Input* weights_input     =element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
		Input* thickness_input   =element->GetInput(ThicknessEnum);                          _assert_(thickness_input);
		Input* thicknessobs_input=element->GetInput(InversionThicknessObsEnum);              _assert_(thicknessobs_input);
		Input* potential_input   =element->GetInput(PotentialEnum);                          _assert_(potential_input);
		Input* vxobs_input       =element->GetInput(BalancethicknessVxObsEnum); _assert_(vxobs_input);
		Input* vyobs_input       =element->GetInput(BalancethicknessVyObsEnum); _assert_(vyobs_input);
		Input* vx_input          =element->GetInput(VxEnum); _assert_(vx_input);
		Input* vy_input          =element->GetInput(VyEnum); _assert_(vy_input);
		Input* nux_input       = element->GetInput(BalancethicknessNuxEnum);   _assert_(nux_input);
		Input* nuy_input       = element->GetInput(BalancethicknessNuyEnum);   _assert_(nuy_input);

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

			gauss->GaussPoint(ig);

			/* Get Jacobian determinant: */
			element->JacobianDeterminant(&Jdet,xyz_list,gauss);

			/*Get all parameters at gaussian point*/
			weights_input->GetInputValue(&weight,gauss,Balancethickness2MisfitEnum);
			thickness_input->GetInputValue(&thickness,gauss);
			thicknessobs_input->GetInputValue(&thicknessobs,gauss);
			potential_input->GetInputValue(&potential,gauss);
			vxobs_input->GetInputValue(&vxobs,gauss);
			vyobs_input->GetInputValue(&vyobs,gauss);
			vx_input->GetInputValue(&vxbar,gauss);
			vy_input->GetInputValue(&vybar,gauss);
			nux_input->GetInputValue(&nux,gauss);
			nuy_input->GetInputValue(&nuy,gauss);

			vxbarobs = nux*vxobs;
			vybarobs = nuy*vyobs;

			/*J = (H^2 - Hobs^2)^2*/
			//J +=0.5*(thickness*thickness - thicknessobs*thicknessobs)*(thickness*thickness - thicknessobs*thicknessobs)*weight*Jdet*gauss->weight;
			/*J = phi^2*/
			//J +=.5*potential*potential*weight*Jdet*gauss->weight; OK
			/*J = (ubar - nux*uobs)^2*/
			J +=0.5*((vxbarobs - vxbar)*(vxbarobs - vxbar) + (vybarobs - vybar)*(vybarobs - vybar))*weight*Jdet*gauss->weight;
		}

		/*clean up and Return: */
		xDelete<IssmDouble>(xyz_list);
		delete gauss;
	}

	/*Sum all J from all cpus of the cluster:*/
	ISSM_MPI_Reduce (&J,&J_sum,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
	ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
	J=J_sum;

	/*Assign output pointers: */
	*presponse=J;

}/*}}}*/
void FemModel::ThicknessAbsGradientx( IssmDouble* pJ){/*{{{*/

	/*output: */
	IssmDouble J=0.;
	IssmDouble J_sum;

	IssmDouble  thickness,thicknessobs,weight;
	IssmDouble  Jdet;
	IssmDouble* xyz_list = NULL;
	IssmDouble  dp[3];

	/*Compute Misfit: */
	for(int i=0;i<elements->Size();i++){
		Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));

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

		/* Get node coordinates*/
		element->GetVerticesCoordinates(&xyz_list);

		/*Retrieve all inputs we will be needing: */
		Input* weights_input   =element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
		Input* thickness_input =element->GetInput(ThicknessEnum);                          _assert_(thickness_input);

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

			gauss->GaussPoint(ig);

			/* Get Jacobian determinant: */
			element->JacobianDeterminant(&Jdet,xyz_list,gauss);

			/*Get all parameters at gaussian point*/
			weights_input->GetInputValue(&weight,gauss,ThicknessAbsMisfitEnum);
			thickness_input->GetInputDerivativeValue(&dp[0],xyz_list,gauss);

			/*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ 
			J+=weight*1/2*(dp[0]*dp[0]+dp[1]*dp[1])*Jdet*gauss->weight;
		}

		/*clean up and Return: */
		xDelete<IssmDouble>(xyz_list);
		delete gauss;
	}

	/*Sum all J from all cpus of the cluster:*/
	ISSM_MPI_Reduce (&J,&J_sum,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
	ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
	J=J_sum;

	/*Assign output pointers: */
	*pJ=J;
}
/*}}}*/
void FemModel::CostFunctionx(IssmDouble* pJ,IssmDouble** pJlist,int* pn){/*{{{*/

	/*Intermediary*/
	int      num_responses;
	int     *responses      = NULL;
	Results *cost_functions = NULL;

	/*Recover parameters*/
	parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
	parameters->FindParam(&responses,NULL,InversionCostFunctionsEnum);

	/*Get the value of all cost functions*/
	this->RequestedOutputsx(&cost_functions,responses,num_responses);

	/*Get and add all contributions one by one*/
	IssmDouble  J=0;
	IssmDouble* Jlist = xNew<IssmDouble>(num_responses);
	for(int i=0;i<num_responses;i++){
		ExternalResult* result=(ExternalResult*)cost_functions->GetObjectByOffset(i);
		Jlist[i] = reCast<IssmDouble>(result->GetValue());
		J       += Jlist[i];
	}
	_assert_(cost_functions->Size()==num_responses);

	/*Assign output pointers: */
	delete cost_functions;
	xDelete<int>(responses);
	if(pJ)     *pJ     = J;
	if(pJlist) *pJlist = Jlist;
	else        xDelete<IssmDouble>(Jlist);
	if(pn)     *pn     = num_responses;
}
/*}}}*/
void FemModel::OutputControlsx(Results **presults){/*{{{*/

	/*parameters: */
	int         num_controls,step;
	IssmDouble  time;
	int        *control_type = NULL;

	/*recover results*/
	Results* results = *presults;
	if(!results) results = new Results();

	/*Get list of Controls*/
	this->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
	this->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
	this->parameters->FindParam(&step,StepEnum);
	this->parameters->FindParam(&time,TimeEnum);

	for(int i=0;i<num_controls;i++){

		int control_enum = control_type[i];
		int gradient_enum;

		switch(i){
			case 0: gradient_enum = Gradient1Enum; break;
			case 1: gradient_enum = Gradient2Enum; break;
			case 2: gradient_enum = Gradient3Enum; break;
			default: _error_("more than 3 controls not implemented yet");
		}

		/*Allocate vector*/
		Vector<IssmPDouble> *vector_control  = new Vector<IssmPDouble>(this->vertices->NumberOfVertices());
		Vector<IssmPDouble> *vector_gradient = new Vector<IssmPDouble>(this->vertices->NumberOfVertices());

		/*Fill in vector*/
		for(int j=0;j<elements->Size();j++){
			Element* element=(Element*)elements->GetObjectByOffset(j);
			element->ControlToVectors(vector_control,vector_gradient,control_enum);
		}
		vector_control->Assemble();
		vector_gradient->Assemble();

		results->AddResult(new GenericExternalResult<Vector<IssmPDouble>*>(results->Size()+1,control_enum,vector_control ,step,time));
		results->AddResult(new GenericExternalResult<Vector<IssmPDouble>*>(results->Size()+1,gradient_enum,vector_gradient,step,time));
	}

	/*Clean up and return*/
	xDelete<int>(control_type);
}
/*}}}*/
#ifdef  _HAVE_DAKOTA_
void FemModel::DakotaResponsesx(double* d_responses,char** responses_descriptors,int numresponsedescriptors,int d_numresponses){/*{{{*/

	int        i,j;
	int        my_rank;

	/*intermediary: */
	char   root[50];
	int    index;
	int    npart;
	double femmodel_response;
	int    flag;
	double *vertex_response   = NULL;
	double *qmu_response      = NULL;
	double *responses_pointer = NULL;

	/*retrieve npart: */
	parameters->FindParam(&npart,QmuNumberofpartitionsEnum);
	my_rank=IssmComm::GetRank();

	/*save the d_responses pointer: */
	responses_pointer=d_responses;

	//watch out, we have more d_numresponses than numresponsedescriptors, because the responses have been expanded if they were scaled. 
	//because we don't know the d_responses descriptors (the scaled ones) we can't key off them, so we will key off the responses_descriptors: */

	for(i=0;i<numresponsedescriptors;i++){

		flag=DescriptorIndex(root,&index,responses_descriptors[i]);

		if(flag==ScaledEnum){

			/*this response was scaled. pick up the response from the inputs: */
			GetVectorFromInputsx(&vertex_response,this, StringToEnumx(root),VertexEnum);

			/*Now, average it onto the partition nodes: */
			AverageOntoPartitionx(&qmu_response,elements,nodes,vertices,loads,materials,parameters,vertex_response);

			/*Copy onto our dakota responses: */
			if(my_rank==0){
				/*plug response: */
				for(j=0;j<npart;j++)responses_pointer[j]=qmu_response[j];

				/*increment response_pointer :*/
				responses_pointer+=npart;
			}

			/*Free ressources:*/
			xDelete<double>(vertex_response);
			xDelete<double>(qmu_response);

		}
		else if (flag==IndexedEnum){

			/*indexed response: plug index into parameters and call response module: */
			parameters->SetParam(index,IndexEnum);

			this->Responsex(&femmodel_response,root);

			if(my_rank==0){
				/*plug response: */
				responses_pointer[0]=femmodel_response;

				/*increment response_pointer :*/
				responses_pointer++;
			}
		}
		else if (flag==NodalEnum){
			_error_("nodal response functions not supported yet!");

			/*increment response_pointer :*/
			responses_pointer++;
		}
		else if (flag==RegularEnum){

			/*perfectly normal response function: */
			this->Responsex(&femmodel_response,root);

			if(my_rank==0){
				/*plug response: */
				responses_pointer[0]=femmodel_response;

				/*increment response_pointer :*/
				responses_pointer++;
			}
		}
		else _error_("flag type " << flag << " not supported yet for response analysis");
	}

	/*Synthesize echo: {{{*/
	if(my_rank==0){
		_printf_("   responses: " << d_numresponses << ": ");
		for(i=0;i<d_numresponses-1;i++)_printf_(d_responses[i] << "|");
		_printf_(d_responses[d_numresponses-1]);
		_printf_("\n");
	}
	/*}}}*/

}
/*}}}*/
#endif
#ifdef _HAVE_GIA_
void FemModel::Deflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt, IssmDouble* x, IssmDouble* y){ /*{{{*/

	/*Go through elements, and add contribution from each element to the deflection vector wg:*/
	for(int i=0;i<elements->Size();i++){
		Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
		element->GiaDeflection(wg,dwgdt, x,y);
	}
}
/*}}}*/
#endif
void FemModel::HydrologyEPLupdateDomainx(void){ /*{{{*/

	Vector<IssmDouble>* mask          = NULL;
	IssmDouble*         serial_mask   = NULL;
	Vector<IssmDouble>* active        = NULL;
	IssmDouble*         serial_active = NULL;
	
	HydrologyDCEfficientAnalysis* effanalysis =  new HydrologyDCEfficientAnalysis();
	/*Step 1: update mask, the mask might be extended by residual and/or using downstream sediment head*/
	mask=new Vector<IssmDouble>(nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum));

	for (int i=0;i<elements->Size();i++){
		Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
		effanalysis->HydrologyEPLGetMask(mask,element);
	}

	/*Assemble and serialize*/
	mask->Assemble();
	serial_mask=mask->ToMPISerial();
	delete mask;

	/*Update Mask*/
	InputUpdateFromVectorx(this,serial_mask,HydrologydcMaskEplactiveNodeEnum,NodeSIdEnum);
	xDelete<IssmDouble>(serial_mask);

	/*Step 2: update node activity. If one element is connected to mask=1, all nodes are active*/
	active=new Vector<IssmDouble>(nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum));
	for (int i=0;i<elements->Size();i++){
		Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
		effanalysis->HydrologyEPLGetActive(active,element);
	}

	/*Assemble and serialize*/
	active->Assemble();
	serial_active=active->ToMPISerial();
	delete active;

	/*Update node activation accordingly*/
	int counter =0;
	for (int i=0;i<nodes->Size();i++){
		Node* node=dynamic_cast<Node*>(nodes->GetObjectByOffset(i));
		if(node->InAnalysis(HydrologyDCEfficientAnalysisEnum)){
			if(serial_active[node->Sid()]==1.){
				node->Activate();
				counter++;
			}
			else{
				node->Deactivate();
			}
		}
	}
	xDelete<IssmDouble>(serial_active);
	delete effanalysis;
	int sum_counter;
	ISSM_MPI_Reduce(&counter,&sum_counter,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,IssmComm::GetComm() );
	ISSM_MPI_Bcast(&sum_counter,1,ISSM_MPI_INT,0,IssmComm::GetComm());                
	counter=sum_counter;
	if(VerboseSolution()) _printf0_("   Number of active nodes in EPL layer: "<< counter <<"\n");

	/*Update dof indexings*/
	this->UpdateConstraintsx();

}
/*}}}*/
void FemModel::UpdateConstraintsL2ProjectionEPLx(void){ /*{{{*/

	Vector<IssmDouble>* active        = NULL;
	IssmDouble*         serial_active = NULL;
	HydrologyDCEfficientAnalysis* effanalysis = new HydrologyDCEfficientAnalysis();

	/*update node activity. If one element is connected to mask=1, all nodes are active*/
	active=new Vector<IssmDouble>(nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum));
	for (int i=0;i<elements->Size();i++){
		Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
		effanalysis->HydrologyEPLGetActive(active,element);
	}

	/*Assemble and serialize*/
	active->Assemble();
	serial_active=active->ToMPISerial();
	delete active;
	delete effanalysis;

	/*Update node activation accordingly*/
	int counter =0;
	for (int i=0;i<nodes->Size();i++){
		Node* node=dynamic_cast<Node*>(nodes->GetObjectByOffset(i));
		if(node->InAnalysis(L2ProjectionEPLAnalysisEnum)){
			if(serial_active[node->Sid()]==1.){
				node->Activate();
				counter++;
			}
			else{
				node->Deactivate();
			}
		}
	}
	xDelete<IssmDouble>(serial_active);
	int sum_counter;
	ISSM_MPI_Reduce(&counter,&sum_counter,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,IssmComm::GetComm() );
	ISSM_MPI_Bcast(&sum_counter,1,ISSM_MPI_INT,0,IssmComm::GetComm());                
	counter=sum_counter;
	if(VerboseSolution()) _printf0_("   Number of active nodes L2 Projection: "<< counter <<"\n");

	/*Update dof indexings*/
	this->UpdateConstraintsx();
}
/*}}}*/
