/*!\file FemModel.c
 * \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 "../Container/Container.h"
#include "../modules/ModelProcessorx/ModelProcessorx.h"
#include "../solutions/solutions.h"
#include "../io/io.h"
#include "./classes.h"
#include "../include/include.h"
#include "../EnumDefinitions/EnumDefinitions.h"
#include "../modules/modules.h"

/*Object constructors and destructor*/
/*FUNCTION FemModel::FemModel(int argc,char** argv){{{*/
FemModel::FemModel(int argc,char** argv,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);

}
/*}}}*/
/*FUNCTION FemModel::FemModel(char* rootpath, char* inputfilename, char* outputfilename, char* petscfilename, char* lockfilename, const int in_solution_type,const int* analyses,const int nummodels){{{*/
FemModel::FemModel(char* rootpath, char* inputfilename, char* outputfilename, char* petscfilename, 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, petscfilename, lockfilename, in_solution_type,analyses,nummodels);

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

	/*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>0){
		_pprintLine_("write lock file:");
		WriteLockFile(lockfilename);
	}

	/*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;

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

	/*Now delete: */
	delete profiler;
	
	/*Finalize PETSC for this model: */
	#ifdef _HAVE_PETSC_
	_pprintLine_("closing Petsc");
	PetscFinalize(); 
	#endif

}
/*}}}*/

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

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

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

	/*intermediary*/
	int         i;
	int         analysis_type;
	FILE       *IOMODEL = NULL;
	FILE       *petscoptionsfid = 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 = pfopen(inputfilename ,"rb");

	/*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,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()) _pprintLine_("   Processing finite element model of analysis " << EnumToStringx(analysis_type_list[i]) << ":");
		analysis_type=analysis_type_list[i];
		this->SetCurrentConfiguration(analysis_type);

		if(i==0){
			if(VerboseMProcessor()) _pprintLine_("      creating vertex degrees of freedom");
			VerticesDofx(vertices,parameters); //only call once, we only have one set of vertices
		}

		if(VerboseMProcessor()) _pprintLine_("      resolving node constraints");
		SpcNodesx(nodes,constraints,parameters,analysis_type); 

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

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

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

	/*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));

	/*Now, deal with petsc options, which need to be put into the parameters dataset: */
	petscoptionsfid=pfopen(petscfilename,"r");
	ParsePetscOptionsx(this->parameters,petscoptionsfid);
	pfclose(petscoptionsfid,petscfilename);
}
/*}}}*/
/*FUNCTION FemModel::OutputResults {{{*/
void FemModel::OutputResults(void){

	_pprintLine_("write results to disk:");

	/*Just call the OutputResultsx module: */
	OutputResultsx(this->elements, this->nodes, this->vertices, this->loads, this->materials, this->parameters,this->results);

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

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

}
/*}}}*/
/*FUNCTION FemModel::SetCurrentConfiguration(int configuration_type,int analysis_type){{{*/
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->nodes->SetCurrentConfiguration(elements,loads, nodes,vertices, materials,parameters);
	this->loads->SetCurrentConfiguration(elements, loads, nodes,vertices, materials,parameters);

	#ifdef _HAVE_PETSC_
	/*take care of petsc options, that depend on this analysis type (present only after model processor)*/
	if(this->parameters->Exist(PetscOptionsStringsEnum)){
		PetscOptionsFromAnalysis(this->parameters,analysis_type);
		if(VerboseSolver()) _pprintLine_("      petsc Options set for analysis type: " << EnumToStringx(analysis_type));
	}
	#endif

}
/*}}}*/
/*FUNCTION FemModel::SetCurrentConfiguration(int configuration_type){{{*/
void FemModel::SetCurrentConfiguration(int configuration_type){
	this->SetCurrentConfiguration(configuration_type,configuration_type);
}
/*}}}*/
/*FUNCTION FemModel::Solve {{{*/
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

	_pprintLine_("call computational core:");

	/*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);

		_pprintLine_("Solution elapsed time  : " << solution_time << "  Seconds");
		_pprintLine_("Solution elapsed flops : " << solution_flops << "  Flops");
		_pprintLine_("Solution memory used   : " << solution_memory << "  Bytes");

		/*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));
	}
}
/*}}}*/

/*Modules:*/
void FemModel::AllocateSystemMatrices(Matrix<IssmDouble>** pKff,Matrix<IssmDouble>** pKfs,Vector<IssmDouble>** pdf,Vector<IssmDouble>** ppf){ /*{{{*/

	/*Intermediary*/
	int  fsize,ssize,flocalsize,slocalsize;
	int  connectivity, numberofdofspernode;
	int  configuration_type;
	int  m,n,M,N;
	int *d_nnz = NULL;
	int *o_nnz = NULL;

	/*output*/
	Matrix<IssmDouble> *Kff  = NULL;
	Matrix<IssmDouble> *Kfs  = NULL;
	Vector<IssmDouble> *pf   = NULL;
	Vector<IssmDouble> *df   = NULL;

	bool oldalloc=false;

	/*retrieve parameters: */
	this->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
	this->parameters->FindParam(&connectivity,MeshAverageVertexConnectivityEnum);

	/*retrieve node info*/
	fsize      = this->nodes->NumberOfDofs(configuration_type,FsetEnum);
	ssize      = this->nodes->NumberOfDofs(configuration_type,SsetEnum);
	flocalsize = this->nodes->NumberOfDofsLocal(configuration_type,FsetEnum);
	slocalsize = this->nodes->NumberOfDofsLocal(configuration_type,SsetEnum);

	numberofdofspernode=this->nodes->MaxNumDofs(configuration_type,GsetEnum);

	if(oldalloc){
		if(pKff) Kff=new Matrix<IssmDouble>(fsize,fsize,connectivity,numberofdofspernode);
		if(pKfs) Kfs=new Matrix<IssmDouble>(fsize,ssize,connectivity,numberofdofspernode);
		if(pdf)  df =new Vector<IssmDouble>(fsize);
		if(ppf)  pf =new Vector<IssmDouble>(fsize);
	}
	else{
		if(pKff){
			m=flocalsize; n=flocalsize; /*local  sizes*/
			M=fsize;      N=fsize;      /*global sizes*/
			this->MatrixNonzeros(&d_nnz,&o_nnz,FsetEnum,FsetEnum);
			Kff=new Matrix<IssmDouble>(m,n,M,N,d_nnz,o_nnz);
			xDelete<int>(d_nnz);
			xDelete<int>(o_nnz);
		}
		if(pKfs){
			m=flocalsize; n=slocalsize; /*local  sizes*/
			M=fsize;      N=ssize;      /*global sizes*/
			this->MatrixNonzeros(&d_nnz,&o_nnz,FsetEnum,SsetEnum);
			Kfs=new Matrix<IssmDouble>(m,n,M,N,d_nnz,o_nnz);
			xDelete<int>(d_nnz);
			xDelete<int>(o_nnz);
		}
		if(pdf) df =new Vector<IssmDouble>(flocalsize,fsize);
		if(ppf) pf =new Vector<IssmDouble>(flocalsize,fsize);
	}

	/*Allocate output pointers*/
	if(pKff) *pKff = Kff;
	if(pKfs) *pKfs = Kfs;
	if(pdf)  *pdf  = df;
	if(ppf)  *ppf  = pf;
}
/*}}}*/
void FemModel::MatrixNonzeros(int** pd_nnz,int** po_nnz,int set1enum,int set2enum){/*{{{*/

	/*Intermediary*/
	int      i,j,k,index,offset,count;
	int      configuration_type;
	int      d_nz,o_nz;
	Element *element            = NULL;
	Load    *load               = NULL;
	int     *head_e             = NULL;
	int     *next_e             = NULL;
	int     *count2offset_e     = NULL;
	int     *head_l             = NULL;
	int     *next_l             = NULL;
	int     *count2offset_l     = NULL;
	int     *sidlist            = NULL;

	/*output*/
	int *d_nnz = NULL;
	int *o_nnz = NULL;

	/*retrive parameters: */
	this->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);

	/*Get vector size and number of nodes*/
	int numnodes            = nodes->NumberOfNodes(configuration_type);
	int numberofdofspernode = nodes->MaxNumDofs(configuration_type,GsetEnum);
	int M                   = nodes->NumberOfDofs(configuration_type,set1enum);
	int N                   = nodes->NumberOfDofs(configuration_type,set2enum);
	int m                   = nodes->NumberOfDofsLocal(configuration_type,set1enum);
	int n                   = nodes->NumberOfDofsLocal(configuration_type,set2enum);
	int numnodesperelement  = elements->MaxNumNodes();
	int numnodesperload     = loads->MaxNumNodes(configuration_type);

	/*First, we are building chaining vectors so that we know what nodes are
	 * connected to what elements. These vectors are such that:
	 *   for(int i=head[id];i!=-1;i=next[i])
	 * will loop over all the elements that are connected to the node number
	 * id*/
	head_e         = xNew<int>(numnodes); for(i=0;i<numnodes;i++) head_e[i]=-1;
	next_e         = xNew<int>(elements->Size()*numnodesperelement);
	count2offset_e = xNew<int>(elements->Size()*numnodesperelement);

	k=0;
	for(i=0;i<elements->Size();i++){
		element = dynamic_cast<Element*>(elements->GetObjectByOffset(i));
		sidlist = xNew<int>(element->GetNumberOfNodes());
		element->GetNodesSidList(sidlist);

		for(j=0;j<element->GetNumberOfNodes();j++){
			index = sidlist[j];
			_assert_(index>=0 && index<numnodes);

			count2offset_e[k]=i;
			next_e[k]=head_e[index];
			head_e[index]=k++;
		}
		for(j=0;j<numnodesperelement-element->GetNumberOfNodes();j++) k++;

		xDelete<int>(sidlist);
	}

	/*Chain for loads*/
	head_l         = xNew<int>(numnodes); for(i=0;i<numnodes;i++) head_l[i]=-1;
	next_l         = xNew<int>(loads->Size(configuration_type)*numnodesperload);
	count2offset_l = xNew<int>(loads->Size(configuration_type)*numnodesperload);
	k=0;
	for(i=0;i<loads->Size();i++){
		load = dynamic_cast<Load*>(loads->GetObjectByOffset(i));
		if(!load->InAnalysis(configuration_type)) continue;
		sidlist = xNew<int>(load->GetNumberOfNodes());
		load->GetNodesSidList(sidlist);

		for(j=0;j<load->GetNumberOfNodes();j++){
			index = sidlist[j];
			_assert_(index>=0 && index<numnodes);

			count2offset_l[k]=i;
			next_l[k]=head_l[index];
			head_l[index]=k++;
		}
		for(j=0;j<numnodesperload-load->GetNumberOfNodes();j++) k++;

		xDelete<int>(sidlist);
	}

	/*OK now count number of dofs and flag each nodes for each node i*/
	bool *flags                  = xNew<bool>(numnodes);
	int  *d_connectivity         = xNewZeroInit<int>(numnodes);
	int  *o_connectivity         = xNewZeroInit<int>(numnodes);
	int  *connectivity_clone     = xNewZeroInit<int>(numnodes);
	int  *all_connectivity_clone = xNewZeroInit<int>(numnodes);

	/*Create connectivity vector*/
	for(i=0;i<nodes->Size();i++){
		Node* node=dynamic_cast<Node*>(nodes->GetObjectByOffset(i));
		if(node->InAnalysis(configuration_type)){

			/*Reinitialize flags to 0*/
			for(j=0;j<numnodes;j++) flags[j]=false;

			/*Loop over elements that hold node number i*/
			//if(head_e[node->Sid()]==-1 && head_l[node->Sid()]==-1){
			//	printf("[%i] vertex %i\n",IssmComm::GetRank(),node->Sid()+1);
			//}
			for(j=head_e[node->Sid()];j!=-1;j=next_e[j]){
				offset=count2offset_e[j];
				element=dynamic_cast<Element*>(elements->GetObjectByOffset(offset));
				element->SetwiseNodeConnectivity(&d_nz,&o_nz,node,flags,set1enum,set2enum);
				if(node->IsClone()){
					connectivity_clone[node->Sid()]+=d_nz+o_nz;
				}
				else{
					d_connectivity[node->Sid()]+=d_nz;
					o_connectivity[node->Sid()]+=o_nz;
				}
			}
			for(j=head_l[node->Sid()];j!=-1;j=next_l[j]){
				offset=count2offset_l[j];
				load=dynamic_cast<Load*>(loads->GetObjectByOffset(offset));
				load->SetwiseNodeConnectivity(&d_nz,&o_nz,node,flags,set1enum,set2enum);
				if(node->IsClone()){
					connectivity_clone[node->Sid()]+=d_nz+o_nz;
				}
				else{
					d_connectivity[node->Sid()]+=d_nz;
					o_connectivity[node->Sid()]+=o_nz;
				}
			}
		}
	}
	xDelete<bool>(flags);
	xDelete<int>(count2offset_e);
	xDelete<int>(head_e);
	xDelete<int>(next_e);
	xDelete<int>(count2offset_l);
	xDelete<int>(head_l);
	xDelete<int>(next_l);

	/*sum over all cpus*/
#ifdef _HAVE_MPI_
	MPI_Allreduce((void*)connectivity_clone,(void*)all_connectivity_clone,numnodes,MPI_INT,MPI_SUM,IssmComm::GetComm());
#endif
	xDelete<int>(connectivity_clone);

	if(set1enum==FsetEnum){
		count=0;
		d_nnz=xNew<int>(m);
		o_nnz=xNew<int>(m);
		for(i=0;i<nodes->Size();i++){
			Node* node=dynamic_cast<Node*>(nodes->GetObjectByOffset(i));
			if(node->InAnalysis(configuration_type) && !node->IsClone()){
				for(j=0;j<node->indexing.fsize;j++){
					_assert_(count<m);
					d_nnz[count]=numberofdofspernode*(d_connectivity[node->Sid()] + all_connectivity_clone[node->Sid()]);
					o_nnz[count]=numberofdofspernode*(o_connectivity[node->Sid()] + all_connectivity_clone[node->Sid()]);
					if(d_nnz[count]>n)   d_nnz[count]=n;
					if(o_nnz[count]>N-n) o_nnz[count]=N-n;
					count++;
				}
			}
		}
		_assert_(m==count);
	}
	else{
		_error_("STOP not implemented");
	}
	xDelete<int>(d_connectivity);
	xDelete<int>(o_connectivity);
	xDelete<int>(all_connectivity_clone);

	/*Allocate ouptput pointer*/
	*pd_nnz=d_nnz;
	*po_nnz=o_nnz;

}/*}}}*/
void FemModel::CreateJacobianMatrixx(Matrix<IssmDouble>** pJff,IssmDouble kmax){/*{{{*/

	int      i,connectivity;
	int      numberofdofspernode;
	int      fsize,configuration_type;
	Element *element = NULL;
	Load    *load    = NULL;
	Matrix<IssmDouble>* Jff = NULL;

	/*Checks*/
	_assert_(nodes && elements);

	/*Recover some parameters*/
	parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
	parameters->FindParam(&connectivity,MeshAverageVertexConnectivityEnum);
	fsize=nodes->NumberOfDofs(configuration_type,FsetEnum);
	numberofdofspernode=nodes->MaxNumDofs(configuration_type,GsetEnum);

	/*Initialize Jacobian Matrix*/
	this->AllocateSystemMatrices(&Jff,NULL,NULL,NULL);

	/*Create and assemble matrix*/
	for(i=0;i<elements->Size();i++){
		element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
		element->CreateJacobianMatrix(Jff);
	}
	for (i=0;i<loads->Size();i++){
		load=(Load*)loads->GetObjectByOffset(i);
		if(load->InAnalysis(configuration_type)) load->CreateJacobianMatrix(Jff);
		if(load->InAnalysis(configuration_type)) load->PenaltyCreateJacobianMatrix(Jff,kmax);
	}
	Jff->Assemble();

	/*Assign output pointer*/
	*pJff=Jff;

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

	int     i;
	Vector<IssmDouble>*     vz        = NULL;
	Vertex *vertex    = NULL;
	IssmDouble *thickness = NULL;
	IssmDouble *bed       = NULL;

	/*get vertex vectors for bed and thickness: */
	GetVectorFromInputsx(&thickness,elements,nodes, vertices, loads, materials, parameters, ThicknessEnum,VertexEnum);
	GetVectorFromInputsx(&bed      ,elements,nodes, vertices, loads, materials, parameters, BedEnum,      VertexEnum);

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

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

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

	/*Update element inputs*/
	InputUpdateFromVectorx(elements,nodes,vertices,loads,materials,parameters,vz,VzMeshEnum,VertexEnum);

	/*Free ressources:*/
	xDelete<IssmDouble>(thickness);
	xDelete<IssmDouble>(bed);
	xdelete(&vz);
	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()) _pprintLine_("   Updating constraints for time: " << time);

	/*First, update dof constraints in nodes, using constraints: */
	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,bool process_units,int weight_index){/*{{{*/

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

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


	switch (response_descriptor_enum){

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

}
/*}}}*/
void FemModel::RequestedOutputsx(int* requested_outputs, int numoutputs){/*{{{*/


	int      output_enum;
	int      step;
	IssmDouble   time;
	IssmDouble   output_value;
	Element *element      = NULL;

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

	/*retrieve Inputs*/
	if(numoutputs){
		for(int i=0;i<numoutputs;i++){
			output_enum=requested_outputs[i];

			switch(output_enum){

				case IceVolumeEnum:
					Responsex(&output_value,"IceVolume",false,0);
					results->AddObject(new GenericExternalResult<double>(results->Size()+1,IceVolumeEnum,reCast<IssmPDouble>(output_value),step,time));
					break;
				case TotalSmbEnum:
					Responsex(&output_value,"TotalSmb",false,0);
					results->AddObject(new GenericExternalResult<double>(results->Size()+1,TotalSmbEnum,reCast<IssmPDouble>(output_value),step,time));
					break;
				case MaxVelEnum:
					Responsex(&output_value,"MaxVel",false,0);
					results->AddObject(new GenericExternalResult<double>(results->Size()+1,MaxVelEnum,reCast<IssmPDouble>(output_value),step,time));
					break;
				default:
					/*create this output in the element inputs, and then transfer to results:*/
					for(int j=0;j<elements->Size();j++){
						element=(Element*)elements->GetObjectByOffset(j);
						element->RequestedOutput(output_enum,step,time);
					}
					break;
			}
		}
	}
}
/*}}}*/
void FemModel::RequestedDependentsx(void){/*{{{*/


	bool        isautodiff      = false;
	IssmDouble  output_value;

	int         num_dependents;
	IssmPDouble *dependents;
	DataSet*    dependent_objects=NULL;

	/*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);
				output_value>>=dependents[i];
			}
		}
		delete dependent_objects;
		#else
		_error_("Should not be requesting dependents when an AD library is not available!");
		#endif
	}
}
/*}}}*/
void FemModel::SystemMatricesx(Matrix<IssmDouble>** pKff, Matrix<IssmDouble>** pKfs, Vector<IssmDouble>** ppf, Vector<IssmDouble>** pdf, IssmDouble* pkmax){/*{{{*/

	/*intermediary: */
	int      i;
	int      configuration_type;
	Element *element = NULL;
	Load    *load    = NULL;

	/*output: */
	Matrix<IssmDouble> *Kff  = NULL;
	Matrix<IssmDouble> *Kfs  = NULL;
	Vector<IssmDouble> *pf   = NULL;
	Vector<IssmDouble> *df   = NULL;
	IssmDouble          kmax = 0;

	/*Display message*/
	if(VerboseModule()) _pprintLine_("   Generating matrices");

	/*retrive parameters: */
	this->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);

	/*First, we might need to do a dry run to get kmax if penalties are employed*/
	if(loads->IsPenalty(configuration_type)){

		/*Allocate Kff_temp*/
		Matrix<IssmDouble> *Kff_temp = NULL;
		this->AllocateSystemMatrices(&Kff_temp,NULL,NULL,NULL);

		/*Get complete stiffness matrix without penalties*/
		for (i=0;i<this->elements->Size();i++){
			element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
			element->CreateKMatrix(Kff_temp,NULL,NULL);
		}
		for (i=0;i<this->loads->Size();i++){
			load=dynamic_cast<Load*>(this->loads->GetObjectByOffset(i));
			if(load->InAnalysis(configuration_type)) load->CreateKMatrix(Kff_temp,NULL);
		}
		Kff_temp->Assemble();

		/*Now, figure out maximum value of stiffness matrix, so that we can penalize it correctly: */
		kmax=Kff_temp->Norm(NORM_INF);
		delete Kff_temp;
	}

	/*Allocate stiffness matrices and load vector*/
	this->AllocateSystemMatrices(&Kff,&Kfs,&df,&pf);

	/*Fill stiffness matrix from elements and loads */
	for (i=0;i<this->elements->Size();i++){
		element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
		element->CreateKMatrix(Kff,Kfs,df);
	}
	for (i=0;i<this->loads->Size();i++){
		load=dynamic_cast<Load*>(this->loads->GetObjectByOffset(i));
		if(load->InAnalysis(configuration_type)) load->CreateKMatrix(Kff,Kfs);
	}

	/*Fill right hand side vector, from elements and loads */
	for (i=0;i<this->elements->Size();i++){
		element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
		element->CreatePVector(pf);
	}
	for (i=0;i<this->loads->Size();i++){
		load=dynamic_cast<Load*>(this->loads->GetObjectByOffset(i));
		if(load->InAnalysis(configuration_type)) load->CreatePVector(pf);
	}

	/*Now deal with penalties (only in loads)*/
	if(loads->IsPenalty(configuration_type)){
		for (i=0;i<this->loads->Size();i++){
			load=dynamic_cast<Load*>(this->loads->GetObjectByOffset(i));
			if(load->InAnalysis(configuration_type)) load->PenaltyCreateKMatrix(Kff,Kfs,kmax);
		}
		for (i=0;i<this->loads->Size();i++){
			load=dynamic_cast<Load*>(this->loads->GetObjectByOffset(i));
			if(load->InAnalysis(configuration_type)) load->PenaltyCreatePVector(pf,kmax);
		}
	}

	/*Assemble matrices and vector*/
	Kff->Assemble();
	Kfs->Assemble();
	pf->Assemble();
	df->Assemble();
	//Kff->AllocationInfo();
	//Kfs->AllocationInfo();

	/*Assign output pointers: */
	if(pKff) *pKff=Kff;
	else      xdelete(&Kff);
	if(pKfs) *pKfs=Kfs;
	else      xdelete(&Kfs);
	if(ppf)  *ppf=pf;
	else      xdelete(&pf);
	if(pdf)  *pdf=df;
	else      xdelete(&df);
	if(pkmax) *pkmax=kmax;
}
/*}}}*/
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: */
#ifdef _HAVE_MPI_
	MPI_Reduce (&min_dt,&node_min_dt,1,MPI_DOUBLE,MPI_MIN,0,IssmComm::GetComm() );
	MPI_Bcast(&node_min_dt,1,MPI_DOUBLE,0,IssmComm::GetComm());
	min_dt=node_min_dt;
#endif

	/*Assign output pointers:*/
	*pdt=min_dt;
}
/*}}}*/
#ifdef _HAVE_RESPONSES_
void FemModel::MassFluxx(IssmDouble* pmass_flux,bool process_units){/*{{{*/

	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,process_units);
				break;
			}
		}
	}

#ifdef _HAVE_MPI_
	MPI_Allreduce ( (void*)&mass_flux,(void*)&all_mass_flux,1,MPI_DOUBLE,MPI_SUM,IssmComm::GetComm());
	mass_flux=all_mass_flux;
#endif

	/*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,bool process_units){/*{{{*/

	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=(Element*)this->elements->GetObjectByOffset(i);
		element->MaxAbsVx(&element_maxabsvx,process_units);
		if(element_maxabsvx>maxabsvx) maxabsvx=element_maxabsvx;
	}

	/*Figure out maximum across the cluster: */
#ifdef _HAVE_MPI_
	MPI_Reduce(&maxabsvx,&node_maxabsvx,1,MPI_DOUBLE,MPI_MAX,0,IssmComm::GetComm() );
	MPI_Bcast(&node_maxabsvx,1,MPI_DOUBLE,0,IssmComm::GetComm());   
	maxabsvx=node_maxabsvx;
#endif

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

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

	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=(Element*)this->elements->GetObjectByOffset(i);
		element->MaxAbsVy(&element_maxabsvy,process_units);
		if(element_maxabsvy>maxabsvy) maxabsvy=element_maxabsvy;
	}

	/*Figure out maximum across the cluster: */
#ifdef _HAVE_MPI_
	MPI_Reduce(&maxabsvy,&node_maxabsvy,1,MPI_DOUBLE,MPI_MAX,0,IssmComm::GetComm() );
	MPI_Bcast(&node_maxabsvy,1,MPI_DOUBLE,0,IssmComm::GetComm());   
	maxabsvy=node_maxabsvy;
#endif

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

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

	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=(Element*)this->elements->GetObjectByOffset(i);
		element->MaxAbsVz(&element_maxabsvz,process_units);
		if(element_maxabsvz>maxabsvz) maxabsvz=element_maxabsvz;
	}

	/*Figure out maximum across the cluster: */
#ifdef _HAVE_MPI_
	MPI_Reduce(&maxabsvz,&node_maxabsvz,1,MPI_DOUBLE,MPI_MAX,0,IssmComm::GetComm() );
	MPI_Bcast(&node_maxabsvz,1,MPI_DOUBLE,0,IssmComm::GetComm());   
	maxabsvz=node_maxabsvz;
#endif

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

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

	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=(Element*)this->elements->GetObjectByOffset(i);
		element->MaxVel(&element_maxvel,process_units);
		if(element_maxvel>maxvel) maxvel=element_maxvel;
	}

	/*Figure out maximum across the cluster: */
#ifdef _HAVE_MPI_
	MPI_Reduce(&maxvel,&node_maxvel,1,MPI_DOUBLE,MPI_MAX,0,IssmComm::GetComm() );
	MPI_Bcast(&node_maxvel,1,MPI_DOUBLE,0,IssmComm::GetComm());   
	maxvel=node_maxvel;
#endif

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

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

	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=(Element*)this->elements->GetObjectByOffset(i);
		element->MaxVx(&element_maxvx,process_units);
		if(element_maxvx>maxvx) maxvx=element_maxvx;
	}

	/*Figure out maximum across the cluster: */
#ifdef _HAVE_MPI_
	MPI_Reduce(&maxvx,&node_maxvx,1,MPI_DOUBLE,MPI_MAX,0,IssmComm::GetComm() );
	MPI_Bcast(&node_maxvx,1,MPI_DOUBLE,0,IssmComm::GetComm());   
	maxvx=node_maxvx;
#endif

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

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

	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=(Element*)this->elements->GetObjectByOffset(i);
		element->MaxVy(&element_maxvy,process_units);
		if(element_maxvy>maxvy) maxvy=element_maxvy;
	}

	/*Figure out maximum across the cluster: */
#ifdef _HAVE_MPI_
	MPI_Reduce(&maxvy,&node_maxvy,1,MPI_DOUBLE,MPI_MAX,0,IssmComm::GetComm() );
	MPI_Bcast(&node_maxvy,1,MPI_DOUBLE,0,IssmComm::GetComm());   
	maxvy=node_maxvy;
#endif

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

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

	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=(Element*)this->elements->GetObjectByOffset(i);
		element->MaxVz(&element_maxvz,process_units);
		if(element_maxvz>maxvz) maxvz=element_maxvz;
	}

	/*Figure out maximum across the cluster: */
#ifdef _HAVE_MPI_
	MPI_Reduce(&maxvz,&node_maxvz,1,MPI_DOUBLE,MPI_MAX,0,IssmComm::GetComm() );
	MPI_Bcast(&node_maxvz,1,MPI_DOUBLE,0,IssmComm::GetComm());   
	maxvz=node_maxvz;
#endif

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

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

	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=(Element*)this->elements->GetObjectByOffset(i);
		element->MinVel(&element_minvel,process_units);
		if(element_minvel<minvel) minvel=element_minvel;
	}

	/*Figure out minimum across the cluster: */
#ifdef _HAVE_MPI_
	MPI_Reduce(&minvel,&node_minvel,1,MPI_DOUBLE,MPI_MAX,0,IssmComm::GetComm() );
	MPI_Bcast(&node_minvel,1,MPI_DOUBLE,0,IssmComm::GetComm());   
	minvel=node_minvel;
#endif

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

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

	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=(Element*)this->elements->GetObjectByOffset(i);
		element->MinVx(&element_minvx,process_units);
		if(element_minvx<minvx) minvx=element_minvx;
	}

	/*Figure out minimum across the cluster: */
#ifdef _HAVE_MPI_
	MPI_Reduce(&minvx,&node_minvx,1,MPI_DOUBLE,MPI_MAX,0,IssmComm::GetComm() );
	MPI_Bcast(&node_minvx,1,MPI_DOUBLE,0,IssmComm::GetComm());   
	minvx=node_minvx;
#endif

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

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

	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=(Element*)this->elements->GetObjectByOffset(i);
		element->MinVy(&element_minvy,process_units);
		if(element_minvy<minvy) minvy=element_minvy;
	}

	/*Figure out minimum across the cluster: */
#ifdef _HAVE_MPI_
	MPI_Reduce(&minvy,&node_minvy,1,MPI_DOUBLE,MPI_MAX,0,IssmComm::GetComm() );
	MPI_Bcast(&node_minvy,1,MPI_DOUBLE,0,IssmComm::GetComm());   
	minvy=node_minvy;
#endif

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

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

	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=(Element*)this->elements->GetObjectByOffset(i);
		element->MinVz(&element_minvz,process_units);
		if(element_minvz<minvz) minvz=element_minvz;
	}

	/*Figure out minimum across the cluster: */
#ifdef _HAVE_MPI_
	MPI_Reduce(&minvz,&node_minvz,1,MPI_DOUBLE,MPI_MAX,0,IssmComm::GetComm() );
	MPI_Bcast(&node_minvz,1,MPI_DOUBLE,0,IssmComm::GetComm());   
	minvz=node_minvz;
#endif

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

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

	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();
	}
#ifdef _HAVE_MPI_
	MPI_Reduce(&local_smb,&total_smb,1,MPI_DOUBLE,MPI_SUM,0,IssmComm::GetComm() );
	MPI_Bcast(&total_smb,1,MPI_DOUBLE,0,IssmComm::GetComm());
#else
	total_smb=local_smb;
#endif

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

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

	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();
	}
	#ifdef _HAVE_MPI_
	MPI_Reduce(&local_ice_volume,&total_ice_volume,1,MPI_DOUBLE,MPI_SUM,0,IssmComm::GetComm() );
	MPI_Bcast(&total_ice_volume,1,MPI_DOUBLE,0,IssmComm::GetComm());
	#else
	total_ice_volume=local_ice_volume;
	#endif

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

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

	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: */
#ifdef _HAVE_MPI_
	MPI_Allreduce ( &found,&sumfound,1,MPI_INT,MPI_SUM,IssmComm::GetComm());
	if(!sumfound)_error_("could not find material with id" << index << " to compute ElementResponse");
#endif

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

	/*Broadcast and plug into response: */
#ifdef _HAVE_MPI_
	MPI_Allreduce ( &cpu_found,&cpu_found,1,MPI_INT,MPI_MAX,IssmComm::GetComm());
	MPI_Bcast(&response,1,MPI_DOUBLE,cpu_found,IssmComm::GetComm()); 
#endif

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

}/*}}}*/
#endif
#ifdef _HAVE_CONTROL_
void FemModel::ThicknessAbsGradientx( IssmDouble* pJ, bool process_units, int weight_index){/*{{{*/


	/*Intermediary*/
	int i;
	Element* element=NULL;

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

	/*Compute Misfit: */
	for (i=0;i<elements->Size();i++){
		element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
		J+=element->ThicknessAbsGradient(process_units,weight_index);
	}

	/*Sum all J from all cpus of the cluster:*/
	#ifdef _HAVE_MPI_
	MPI_Reduce (&J,&J_sum,1,MPI_DOUBLE,MPI_SUM,0,IssmComm::GetComm() );
	MPI_Bcast(&J_sum,1,MPI_DOUBLE,0,IssmComm::GetComm());
	J=J_sum;
	#endif

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


	/*Intermediary*/
	int      num_responses;
	Element *element       = NULL;
	int     *responses     = NULL;

	/*output: */
	IssmDouble J,Jplus;

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

	/*Get response*/
	J=0;
	for(int i=0;i<num_responses;i++){
		this->Responsex(&Jplus,EnumToStringx(responses[i]),false,i); //False means DO NOT process units
		J+=Jplus;
	}

	/*Assign output pointers: */
	xDelete<int>(responses);
	*pJ=J;
}
/*}}}*/
#endif
#ifdef  _HAVE_DAKOTA_
void FemModel::DakotaResponsesx(double* d_responses,char** responses_descriptors,int numresponsedescriptors,int d_numresponses){/*{{{*/


	int        i,j;
	int        my_rank;
	bool       process_units = true;

	/*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,elements,nodes, vertices, loads, materials, parameters, 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,process_units,0);//0 is the index for weights

			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,process_units,0);//0 is the weight index

			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){
		_printString_("   responses: " << d_numresponses << ": ");
		for(i=0;i<d_numresponses-1;i++)_printString_(d_responses[i] << "|");
		_printString_(d_responses[d_numresponses-1]);
		_printLine_("");
	}
	/*}}}*/

}
/*}}}*/
#endif
