/*!\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 "../analyses/analyses.h"
#include "../shared/io/io.h"
#include "./classes.h"
#include "../shared/Enum/Enum.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/OutputResultsx/OutputResultsx.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*/
/*FUNCTION FemModel::FemModel(int argc,char** argv,ISSM_MPI_Comm incomm){{{*/
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);

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

}
/*}}}*/
/*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){
		_printf0_("write lock file:\n");
		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
	_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");

	/*Now delete: */
	delete profiler;

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

}
/*}}}*/

/*Object management*/
/*FUNCTION FemModel::Echo {{{*/
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");

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

	/*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()) _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 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 toolkits options, which need to be put into the parameters dataset: */
	toolkitsoptionsfid=pfopen(toolkitsfilename,"r");
	ParseToolkitsOptionsx(this->parameters,toolkitsoptionsfid);
	pfclose(toolkitsoptionsfid,toolkitsfilename);
}
/*}}}*/
/*FUNCTION FemModel::OutputResults {{{*/
void FemModel::OutputResults(void){

	_printf0_("write results to disk:\n");

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

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

	_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

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

	Vector<IssmDouble> *vz        = NULL;
	IssmDouble         *thickness = NULL;
	IssmDouble         *bed       = NULL;

	/*get vertex vectors for bed and thickness: */
	GetVectorFromInputsx(&thickness,this, ThicknessEnum,VertexEnum);
	GetVectorFromInputsx(&bed      ,this, BedEnum,      VertexEnum);

	/*Allocate vector*/
	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(vz,parameters,thickness,bed);
	}

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

	/*Update element inputs*/
	InputUpdateFromVectorx(this,vz,VzMeshEnum,VertexPIdEnum);

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

	Element   *element = NULL;
	IssmDouble time;
	int        analysis_type;

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

	/*start module: */
	if(VerboseModule()) _printf0_("   Updating constraints for time: " << time << "\n");

	/*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 weight_index){/*{{{*/

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

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

	switch (response_descriptor_enum){

		#ifdef _HAVE_RESPONSES_
		case IceVolumeEnum:              this->IceVolumex(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;
		#ifdef _HAVE_CONTROL_
		case SurfaceAbsVelMisfitEnum:    SurfaceAbsVelMisfitx(responses, elements,nodes, vertices, loads, materials, parameters,weight_index); break;
		case SurfaceRelVelMisfitEnum:    SurfaceRelVelMisfitx(responses, elements,nodes, vertices, loads, materials, parameters,weight_index); break;
		case SurfaceLogVelMisfitEnum:    SurfaceLogVelMisfitx(responses, elements,nodes, vertices, loads, materials, parameters,weight_index); break;
		case SurfaceLogVxVyMisfitEnum:   SurfaceLogVxVyMisfitx(responses, elements,nodes, vertices, loads, materials, parameters,weight_index); break;
		case SurfaceAverageVelMisfitEnum:SurfaceAverageVelMisfitx(responses,this,weight_index); break;
		case ThicknessAbsMisfitEnum:     ThicknessAbsMisfitx(responses, elements,nodes, vertices, loads, materials, parameters,weight_index); break;
		case ThicknessAbsGradientEnum:   this->ThicknessAbsGradientx(responses,weight_index); break;
		case ThicknessAlongGradientEnum: ThicknessAlongGradientx(responses, elements,nodes, vertices, loads, materials, parameters,weight_index); break;
		case ThicknessAcrossGradientEnum:ThicknessAcrossGradientx(responses, elements,nodes, vertices, loads, materials, parameters,weight_index); break;
		case RheologyBbarAbsGradientEnum:RheologyBbarAbsGradientx(responses, elements,nodes, vertices, loads, materials, parameters,weight_index); break;
		case DragCoefficientAbsGradientEnum:DragCoefficientAbsGradientx(responses, elements,nodes, vertices, loads, materials, parameters,weight_index); break;
		case BalancethicknessMisfitEnum:BalancethicknessMisfitx(responses,weight_index); break;
		#endif
		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;
		#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",0);
					results->AddObject(new GenericExternalResult<double>(results->Size()+1,IceVolumeEnum,reCast<IssmPDouble>(output_value),step,time));
					break;
				case TotalSmbEnum:
					Responsex(&output_value,"TotalSmb",0);
					results->AddObject(new GenericExternalResult<double>(results->Size()+1,TotalSmbEnum,reCast<IssmPDouble>(output_value),step,time));
					break;
				case MaxVelEnum:
					Responsex(&output_value,"MaxVel",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;
	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;
		#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;
}
/*}}}*/
#ifdef _HAVE_RESPONSES_
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=(Element*)this->elements->GetObjectByOffset(i);
		element->MaxAbsVx(&element_maxabsvx);
		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=(Element*)this->elements->GetObjectByOffset(i);
		element->MaxAbsVy(&element_maxabsvy);
		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=(Element*)this->elements->GetObjectByOffset(i);
		element->MaxAbsVz(&element_maxabsvz);
		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=(Element*)this->elements->GetObjectByOffset(i);
		element->MaxVel(&element_maxvel);
		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=(Element*)this->elements->GetObjectByOffset(i);
		element->MaxVx(&element_maxvx);
		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=(Element*)this->elements->GetObjectByOffset(i);
		element->MaxVy(&element_maxvy);
		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=(Element*)this->elements->GetObjectByOffset(i);
		element->MaxVz(&element_maxvz);
		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=(Element*)this->elements->GetObjectByOffset(i);
		element->MinVel(&element_minvel);
		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=(Element*)this->elements->GetObjectByOffset(i);
		element->MinVx(&element_minvx);
		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=(Element*)this->elements->GetObjectByOffset(i);
		element->MinVy(&element_minvy);
		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=(Element*)this->elements->GetObjectByOffset(i);
		element->MinVz(&element_minvz);
		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::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::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;

}/*}}}*/
#endif
#ifdef _HAVE_CONTROL_
void FemModel::BalancethicknessMisfitx(IssmDouble* presponse,int weight_index){/*{{{*/

	IssmDouble J = 0;
	IssmDouble J_sum;

	for(int i=0;i<this->elements->Size();i++){
		Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
		J+=element->BalancethicknessMisfit(weight_index);
	}
	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, 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(weight_index);
	}

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

	/*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;
	//_printf_(true,"list of misfits: ");
	for(int i=0;i<num_responses;i++){
		this->Responsex(&Jplus,EnumToStringx(responses[i]),i);
		//_printf_(true," %12.7g ",Jplus);
		J+=Jplus;
	}
	//_printf_(true," \n");

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

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

  int      i;

	/*intermediary: */
	Element *element     = NULL;

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

	Vector<IssmDouble>* mask          = NULL;
	IssmDouble*         serial_mask   = NULL;
	Vector<IssmDouble>* active        = NULL;
	IssmDouble*         serial_active = NULL;

	/*Step 1: update maks, 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));
		element->HydrologyEPLGetMask(mask);
	}

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

	/*Update Mask*/
	InputUpdateFromVectorx(this,serial_mask,HydrologydcMaskEplactiveEnum,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));
		element->HydrologyEPLGetActive(active);
	}

	/*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);
	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::HydrologyTransferx(void){ /*{{{*/

	Vector<IssmDouble>* transferg=NULL; 

	/*Vector allocation*/
	transferg=new Vector<IssmDouble>(nodes->NumberOfNodes());

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

	/*Update Inputs*/
	InputUpdateFromVectorx(this,transferg,WaterTransferEnum,NodesEnum);
	delete transferg;
}
/*}}}*/
#endif
#ifdef _HAVE_THERMAL_
void FemModel::UpdateThermalBasalConstraintsx(void){ /*{{{*/

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

}
/*}}}*/
#endif
