FemModel.cpp

Go to the documentation of this file.

 
#ifdef HAVE_CONFIG_H
   #include <config.h>
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif
 
#include <stdio.h>
#include <math.h>
#include "../cores/cores.h"
#include "../shared/io/io.h"
#include "./classes.h"
#include "./modules/modules.h"
#include "../shared/Enum/Enum.h"
#include "../analyses/analyses.h"
#include "./Inputs2/DatasetInput2.h"
#include "./Inputs2/ElementInput2.h"
#include "./Inputs2/TransientInput2.h"
 
#if _HAVE_CODIPACK_
#include <sstream> // for output of the CoDiPack tape
extern CoDi_global codi_global;
#endif
 
#if defined(_HAVE_NEOPZ_) && !defined(_HAVE_AD_)
#include <TPZRefPatternDataBase.h>
#endif
 
/*module includes: {{{*/
#include "../modules/ModelProcessorx/ModelProcessorx.h"
#include "../modules/SpcNodesx/SpcNodesx.h"
#include "../modules/ConfigureObjectsx/ConfigureObjectsx.h"
#include "../modules/ParseToolkitsOptionsx/ParseToolkitsOptionsx.h"
#include "../modules/GetVectorFromInputsx/GetVectorFromInputsx.h"
#include "../modules/InputUpdateFromVectorx/InputUpdateFromVectorx.h"
#include "../modules/NodesDofx/NodesDofx.h"
#include "../modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.h"
#include "../modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.h"
#include "../modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.h"
#include "../modules/SurfaceLogVxVyMisfitx/SurfaceLogVxVyMisfitx.h"
#include "../modules/SurfaceAverageVelMisfitx/SurfaceAverageVelMisfitx.h"
#include "../modules/ThicknessAbsMisfitx/ThicknessAbsMisfitx.h"
#include "../modules/ThicknessAlongGradientx/ThicknessAlongGradientx.h"
#include "../modules/ThicknessAcrossGradientx/ThicknessAcrossGradientx.h"
#include "../modules/RheologyBbarAbsGradientx/RheologyBbarAbsGradientx.h"
#include "../modules/DragCoefficientAbsGradientx/DragCoefficientAbsGradientx.h"
#include "../modules/NodalValuex/NodalValuex.h"
#include "../modules/GetVectorFromInputsx/GetVectorFromInputsx.h"
#include "../modules/AverageOntoPartitionx/AverageOntoPartitionx.h"
/*}}}*/
 
/*Object constructors and destructor*/
FemModel::FemModel(int argc,char** argv,ISSM_MPI_Comm incomm,bool trace){/*{{{*/
 
   /*configuration: */
   int  solution_type,amrtype,amr_frequency;
   int  ierr;
 
   /*File names*/
   char *lockfilename   = NULL;
   char *binfilename    = NULL;
   char *outbinfilename = NULL;
   char *petscfilename  = NULL;
   char *restartfilename  = NULL;
   char *rootpath       = NULL;
 
   /*First things first, store the communicator, and set it as a global variable: */
   IssmComm::SetComm(incomm);
 
   /*Now, initialize PETSC: */
   #ifdef _HAVE_PETSC_
   PETSC_COMM_WORLD=incomm;
   ierr=PetscInitialize(&argc,&argv,(char*)0,"");  if(ierr) _error_("Could not initialize Petsc");
   #endif
 
   /*Start profiler: */
   this->profiler=new Profiler();
   profiler->Start(TOTAL);
 
   /*From command line arguments, retrieve different filenames needed to create the FemModel: */
   ProcessArguments(&solution_type,&binfilename,&outbinfilename,&petscfilename,&lockfilename,&restartfilename,&rootpath,argc,argv);
 
   /*Create femmodel from input files: */
   profiler->Start(MPROCESSOR);
   this->InitFromFiles(rootpath,binfilename,outbinfilename,petscfilename,lockfilename,restartfilename, solution_type,trace,NULL);
   profiler->Stop(MPROCESSOR);
 
   /*Save communicator in the parameters dataset: */
   this->parameters->AddObject(new GenericParam<ISSM_MPI_Comm>(incomm,FemModelCommEnum));
 
   /*AMR stuff*/
   this->parameters->FindParam(&amr_frequency,TransientAmrFrequencyEnum);
   this->parameters->FindParam(&amr_frequency,TransientAmrFrequencyEnum);
   #if defined(_HAVE_NEOPZ_) && !defined(_HAVE_AD_)
   this->amr = NULL;
   #endif
   #if defined(_HAVE_BAMG_) && !defined(_HAVE_AD_)
   this->amrbamg = NULL;
   #endif
   #if !defined(_HAVE_AD_)
   if(amr_frequency && solution_type==TransientSolutionEnum){
      /*Verifications. AMR supports SSA, P1 and horizontal 2D domain*/
      bool isSSA;
      int domaintype,element_type;
      this->analysis_counter=-1;
      this->parameters->FindParam(&isSSA,FlowequationIsSSAEnum);
      this->parameters->FindParam(&domaintype,DomainTypeEnum);
      for(int i=0;i<this->nummodels;i++) if(this->analysis_type_list[i]==StressbalanceAnalysisEnum){analysis_counter=i;break;}
      if(analysis_counter==-1) _error_("Could not find alias for analysis_type StressbalanceAnalysisEnum in list of FemModel analyses\n");
      for(int i=0;i<this->elements->Size();i++){
         Element* element  = xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
         element_type      = element->element_type_list[analysis_counter];
         if(element_type!=P1Enum) _error_("Element type "<<EnumToStringx(element_type)<<" not supported with AMR yet!\n");
      }
      if(!isSSA) _error_("Flow equation not supported with AMR yet!\n ");
      if(domaintype!=Domain2DhorizontalEnum) _error_("Domain "<<EnumToStringx(domaintype)<<" not supported with AMR yet!\n");
 
      this->parameters->FindParam(&amrtype,AmrTypeEnum);
      switch(amrtype){
 
         #if defined(_HAVE_NEOPZ_)
         case AmrNeopzEnum: this->InitializeAdaptiveRefinementNeopz(); break;
         #endif
 
         #if defined(_HAVE_BAMG_)
         case AmrBamgEnum: this->InitializeAdaptiveRefinementBamg(); break;
         #endif
 
         default: _error_("not implemented yet");
      }
   }
   #endif
 
   /*Free resources */
   xDelete<char>(lockfilename);
   xDelete<char>(binfilename);
   xDelete<char>(outbinfilename);
   xDelete<char>(petscfilename);
   xDelete<char>(restartfilename);
   xDelete<char>(rootpath);
 
}
/*}}}*/
FemModel::FemModel(char* rootpath, char* inputfilename, char* outputfilename, char* toolkitsfilename, char* lockfilename, char* restartfilename, ISSM_MPI_Comm incomm, int solution_type,IssmPDouble* X){ /*{{{*/
 
   bool traceon=true;
   this->profiler=NULL; /*avoid leak, as we are not using the profiler ever in ad control run. */
 
   /*Store the communicator, but do not set it as a global variable, as this has already
    * been done by the FemModel that called this copy constructor: */
   IssmComm::SetComm(incomm);
 
   /*Create femmodel from input files, with trace activated: */
   profiler->Start(MPROCESSOR);
   this->InitFromFiles(rootpath,inputfilename,outputfilename,toolkitsfilename,lockfilename,restartfilename, solution_type,traceon,X);
   profiler->Stop(MPROCESSOR);
 
   #if defined(_HAVE_NEOPZ_) && !defined(_HAVE_AD_)
   this->amr = NULL;
   #endif
   #if defined(_HAVE_BAMG_) && !defined(_HAVE_AD_)
   this->amrbamg = NULL;
   #endif
 
   /*Save communicator in the parameters dataset: */
   this->parameters->AddObject(new GenericParam<ISSM_MPI_Comm>(incomm,FemModelCommEnum));
}
/*}}}*/
FemModel::~FemModel(){/*{{{*/
 
   /*Intermediary*/
   FILE *output_fid;
   char *outbinfilename = NULL;
   char *lockfilename   = NULL;
 
   #ifndef _HAVE_JAVASCRIPT_
   if(this->parameters->Exist(OutputFileNameEnum)) this->parameters->FindParam(&outbinfilename,OutputFileNameEnum);
   if(this->parameters->Exist(LockFileNameEnum)) this->parameters->FindParam(&lockfilename,LockFileNameEnum);
   #endif
 
   /*Delete all the datasets: */
   if(analysis_type_list)xDelete<int>(analysis_type_list);
   if(outbinfilename)xDelete<char>(outbinfilename);
   if(lockfilename)xDelete<char>(lockfilename);
   if(elements)delete elements;
   if(vertices)delete vertices;
   if(this->constraints_list && this->nummodels){
      for(int i=0;i<this->nummodels;i++) delete this->constraints_list[i];
      xDelete<Constraints*>(constraints_list);
   }
   if(this->loads_list && this->nummodels){
      for(int i=0;i<this->nummodels;i++) delete this->loads_list[i];
      xDelete<Loads*>(loads_list);
   }
   if(this->nodes_list && this->nummodels){
      for(int i=0;i<this->nummodels;i++) delete this->nodes_list[i];
      xDelete<Nodes*>(nodes_list);
   }
   if(materials)delete materials;
   if(parameters)delete parameters;
   if(inputs2)delete inputs2;
   if(results)delete results;
 
   #if defined(_HAVE_NEOPZ_) && !defined(_HAVE_AD_)
   if(amr)delete amr;
   #endif
 
   #if defined(_HAVE_BAMG_) && !defined(_HAVE_AD_)
   if(amrbamg)delete amrbamg;
   #endif
 
   /*Now delete: */
   if(profiler)delete profiler;
}/*}}}*/
 
/*Object management*/
int FemModel::AnalysisIndex(int analysis_enum){/*{{{*/
 
   /*Checks in debugging mode*/
   _assert_(this->analysis_type_list);
 
   /*Find analysis in list*/
   for(int i=0;i<this->nummodels;i++){
      if(this->analysis_type_list[i]==analysis_enum){
         return i;
         break;
      }
   }
 
   /*If you reach this point, analysis has not been found*/
   _error_("Could not find index of analysis " << EnumToStringx(analysis_enum) << " in list of FemModel analyses");
 
 
}/*}}}*/
void FemModel::CheckPoint(void){/*{{{*/
 
   FILE* restartfid=NULL;
   char* restartfilename = NULL;
   int   femmodel_size;
   char* femmodel_buffer=NULL;
   char* femmodel_buffer_ini=NULL;
 
   /*First, recover the name of the restart file: */
   parameters->FindParam(&restartfilename,RestartFileNameEnum);
 
   /*Open file for writing: */
   restartfid=pfopen(restartfilename,"wb");
 
   /*Initialize: */
   femmodel_size=0;
 
   /*Create buffer to hold marshalled femmodel: */
   this->Marshall(NULL,&femmodel_size,MARSHALLING_SIZE);
   femmodel_buffer=xNew<char>(femmodel_size);
   /*Keep track of initial position of femmodel_buffer: */
   femmodel_buffer_ini=femmodel_buffer;
 
   /*Marshall:*/
   this->Marshall(&femmodel_buffer,NULL,MARSHALLING_FORWARD);
 
   /*Reset position of buffer: */
   femmodel_buffer=femmodel_buffer_ini;
 
   /*write buffer: */
   fwrite(femmodel_buffer,femmodel_size,sizeof(char),restartfid);
 
   /*Done, close file :*/
   pfclose(restartfid,restartfilename);
 
   /*Free ressources: */
   xDelete<char>(femmodel_buffer);
   xDelete<char>(restartfilename);
 
}
/*}}}*/
void FemModel::CleanUp(void){/*{{{*/
 
   /*Intermediary*/
   char *lockfilename   = NULL;
   bool  waitonlock     = false;
 
   /*Write lock file if requested: */
   this->parameters->FindParam(&waitonlock,SettingsWaitonlockEnum);
   this->parameters->FindParam(&lockfilename,LockFileNameEnum);
   if(waitonlock){
      _printf0_("write lock file:\n");
      WriteLockFile(lockfilename);
   }
 
   /*Before we delete the profiler, report statistics for this run: */
   profiler->Stop(TOTAL);  //final tagging
 
   _printf0_("\n");
   _printf0_("   "<<setw(40)<<left<<"FemModel initialization elapsed time:"<<setw(7)<<profiler->TotalTime(MPROCESSOR) << "\n");
   /*Total times*/
   _printf0_("   "<<setw(40)<<left<<"Total Core solution elapsed time:"<<setw(7)<<profiler->TotalTime(CORE) << "\n");
 
   /*Linear solver only*/
   _printf0_("   "<<setw(40)<<left<<"Linear solver elapsed time:"<<setw(7)<<profiler->TotalTime(SOLVER) << " ("<<setprecision(2)<<profiler->TotalTime(SOLVER)/profiler->TotalTime(CORE)*100.<<"%)\n");
   _printf0_("\n");
   _printf0_("   Total elapsed time: "
            <<profiler->TotalTimeModHour(TOTAL)<<" hrs "
            <<profiler->TotalTimeModMin(TOTAL)<<" min "
            <<profiler->TotalTimeModSec(TOTAL)<<" sec"
            );
   _printf0_("\n");
 
   /*Finalize PETSC for this model: */
   #ifdef _HAVE_PETSC_
   //_printf0_("closing PETSc\n");
   PetscFinalize();
   #endif
 
   /*Clean up*/
   xDelete<char>(lockfilename);
} /*}}}*/
FemModel* FemModel::copy(void){/*{{{*/
 
   FemModel* output=NULL;
   int       i;
   int       analysis_type;
 
   output=new FemModel(*this); //Use default copy constructor.
 
   output->nummodels = this->nummodels;
   output->solution_type = this->solution_type;
   output->analysis_counter = this->analysis_counter;
 
   /*Now, deep copy arrays: */
   output->analysis_type_list=xNew<int>(nummodels);
   xMemCpy<int>(output->analysis_type_list,this->analysis_type_list,this->nummodels);
 
   /*Analysis dependent arrays*/
   output->constraints_list=xNew<Constraints*>(this->nummodels);
   output->loads_list=xNew<Loads*>(this->nummodels);
   output->nodes_list=xNew<Nodes*>(this->nummodels);
 
   output->profiler=static_cast<Profiler*>(this->profiler->copy());
 
   output->materials=static_cast<Materials*>(this->materials->Copy());
   output->parameters=static_cast<Parameters*>(this->parameters->Copy());
   output->inputs2=static_cast<Inputs2*>(this->inputs2->Copy());
   output->results=static_cast<Results*>(this->results->Copy());
   output->vertices=static_cast<Vertices*>(this->vertices->Copy());
   output->elements=static_cast<Elements*>(this->elements->Copy());
 
   /*reset hooks for elements, loads and nodes: */
   output->elements->ResetHooks();
   output->loads->ResetHooks();
   output->materials->ResetHooks();
 
   /*do the post-processing of the datasets to get an FemModel that can actually run analyses: */
   for(i=0;i<nummodels;i++){
      output->constraints_list[i] = static_cast<Constraints*>(this->constraints_list[i]->Copy());
      output->loads_list[i] = static_cast<Loads*>(this->loads_list[i]->Copy());
      output->nodes_list[i] = static_cast<Nodes*>(this->nodes_list[i]->Copy());
      analysis_type=output->analysis_type_list[i];
      output->SetCurrentConfiguration(analysis_type);
      SpcNodesx(output->nodes_list[i],output->constraints_list[i],output->parameters);
      NodesDofx(output->nodes_list[i],output->parameters);
      ConfigureObjectsx(output->elements,output->loads_list[i],output->nodes_list[i],output->vertices,output->materials,output->parameters,output->inputs2);
   }
 
   /*AMR, no copy for now*/
   #if defined(_HAVE_NEOPZ_) && !defined(_HAVE_AD_)
   this->amr = NULL;
   #endif
   #if defined(_HAVE_BAMG_) && !defined(_HAVE_AD_)
   this->amrbamg = NULL;
   #endif
 
   /*Reset current configuration: */
   analysis_type=output->analysis_type_list[analysis_counter];
   output->SetCurrentConfiguration(analysis_type);
 
   return output;
}
/*}}}*/
void FemModel::Echo(void){/*{{{*/
 
   _printf_("FemModel echo: \n");
   _printf_("   number of fem models: " << nummodels << "\n");
   _printf_("   analysis_type_list: \n");
   for(int i=0;i<nummodels;i++)_printf_("     " << i << ": " << EnumToStringx(analysis_type_list[i]) << "\n");
   _printf_("   current analysis_type: \n");
   _printf_("     " << analysis_counter << ": " << EnumToStringx(analysis_type_list[analysis_counter]) << "\n");
 
}
/*}}}*/
void FemModel::InitFromFiles(char* rootpath, char* inputfilename, char* outputfilename, char* toolkitsfilename, char* lockfilename, char* restartfilename, const int in_solution_type,bool trace,IssmPDouble* X){/*{{{*/
 
   /*intermediary*/
   FILE *IOMODEL            = NULL;
   FILE *toolkitsoptionsfid = NULL;
 
   /*recover my_rank:*/
   int my_rank=IssmComm::GetRank();
 
   /*Open input file descriptor on cpu 0: */
   if(my_rank==0) IOMODEL = pfopen0(inputfilename ,"rb");
 
   /*Open toolkits file: */
   toolkitsoptionsfid=pfopen(toolkitsfilename,"r");
 
   /*Now, go create FemModel:*/
   this->InitFromFids(rootpath,IOMODEL,toolkitsoptionsfid,in_solution_type,trace,X);
 
   /*Close input file and toolkits file descriptors: */
   if(my_rank==0) pfclose(IOMODEL,inputfilename);
   pfclose(toolkitsoptionsfid,toolkitsfilename);
 
   /*Now save all of these file names into parameters, you never know when you might need them: */
   this->parameters->AddObject(new StringParam(ToolkitsFileNameEnum,toolkitsfilename));
   this->parameters->AddObject(new StringParam(RootPathEnum,rootpath));
   this->parameters->AddObject(new StringParam(InputFileNameEnum,inputfilename));
   this->parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename));
   this->parameters->AddObject(new StringParam(LockFileNameEnum,lockfilename));
   this->parameters->AddObject(new StringParam(RestartFileNameEnum,restartfilename));
 
}/*}}}*/
void FemModel::InitFromFids(char* rootpath, FILE* IOMODEL, FILE* toolkitsoptionsfid, int in_solution_type, bool trace, IssmPDouble* X){/*{{{*/
 
   /*Initialize internal data: */
   this->solution_type    = in_solution_type;
   this->analysis_counter = -1;
   this->results          = new Results(); //not initialized by CreateDataSets
 
   /*create IoModel */
   IoModel* iomodel = new IoModel(IOMODEL,in_solution_type,trace,X);
 
   /*Figure out what analyses are activated for this solution*/
   SolutionAnalysesList(&this->analysis_type_list,&this->nummodels,iomodel,this->solution_type);
 
   /*create datasets for all analyses*/
   ModelProcessorx(&this->elements,&this->nodes_list,&this->vertices,&this->materials,&this->constraints_list,&this->loads_list,&this->parameters,&this->inputs2,iomodel,toolkitsoptionsfid,rootpath,this->solution_type,this->nummodels,this->analysis_type_list);
 
   /*do the post-processing of the datasets to get an FemModel that can actually run analyses: */
   for(int i=0;i<nummodels;i++){
 
      if(VerboseMProcessor()) _printf0_("   Processing finite element model of analysis " << EnumToStringx(analysis_type_list[i]) << ":\n");
      this->SetCurrentConfiguration(analysis_type_list[i]);
 
      if(VerboseMProcessor()) _printf0_("      configuring element and loads\n");
      ConfigureObjectsx(this->elements,this->loads,this->nodes,this->vertices,this->materials,this->parameters,this->inputs2);
 
      if(i==0){
         if(VerboseMProcessor()) _printf0_("      detecting active vertices\n");
         GetMaskOfIceVerticesLSMx0(this);
      }
 
      if(VerboseMProcessor()) _printf0_("      resolving node constraints\n");
      SpcNodesx(nodes,this->constraints,parameters);
 
      if(VerboseMProcessor()) _printf0_("      creating nodal degrees of freedom\n");
      NodesDofx(nodes,parameters);
   }
 
   /*Clean up*/
   delete iomodel;
}/*}}}*/
void FemModel::Marshall(char** pmarshalled_data, int* pmarshalled_data_size, int marshall_direction){ /*{{{*/
 
   int       i;
   int       analysis_type;
 
   if(marshall_direction==MARSHALLING_BACKWARD){
      delete this->materials;
      delete this->parameters;
      delete this->inputs2;
      if(this->constraints_list && this->nummodels){
         for(i=0;i<this->nummodels;i++) delete this->constraints_list[i];
         xDelete<Constraints*>(constraints_list);
      }
      if(this->loads_list && this->nummodels){
         for(i=0;i<this->nummodels;i++) delete this->loads_list[i];
         xDelete<Loads*>(loads_list);
      }
      if(this->nodes_list && this->nummodels){
         for(i=0;i<this->nummodels;i++) delete this->nodes_list[i];
         xDelete<Nodes*>(nodes_list);
      }
      delete this->results;
      delete this->vertices;
      delete this->elements;
      xDelete<int>(this->analysis_type_list);
 
      this->materials   = new Materials();
      this->parameters  = new Parameters();
      this->inputs2     = new Inputs2();
      this->results     = new Results();
      this->nodes       = new Nodes();
      this->vertices    = new Vertices();
      this->elements    = new Elements();
   }
 
   MARSHALLING_ENUM(FemModelEnum);
 
   MARSHALLING(solution_type);
   MARSHALLING(analysis_counter);
   MARSHALLING(nummodels);
   MARSHALLING_DYNAMIC(analysis_type_list,int,nummodels);
 
   this->materials->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
   this->parameters->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
   this->inputs2->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
   this->results->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
   this->vertices->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
   this->elements->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
 
   if(marshall_direction==MARSHALLING_BACKWARD){
      this->constraints_list = xNew<Constraints*>(this->nummodels);
      for(i=0;i<nummodels;i++) this->constraints_list[i] = new Constraints();
      this->loads_list = xNew<Loads*>(this->nummodels);
      for(i=0;i<nummodels;i++) this->loads_list[i] = new Loads();
      this->nodes_list = xNew<Nodes*>(this->nummodels);
      for(i=0;i<nummodels;i++) this->nodes_list[i] = new Nodes();
   }
 
   for(i=0;i<nummodels;i++){
      this->constraints_list[i]->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
      this->loads_list[i]->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
      this->nodes_list[i]->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
   }
 
   if(marshall_direction==MARSHALLING_BACKWARD){
      /*reset hooks for elements, loads and nodes:*/
      this->elements->ResetHooks();
      this->materials->ResetHooks();
 
      /*do the post-processing of the datasets to get an FemModel that can actually run analyses:*/
      for(i=0;i<nummodels;i++){
         this->loads_list[i]->ResetHooks();
         analysis_type=this->analysis_type_list[i];
         SetCurrentConfiguration(analysis_type);
         SpcNodesx(this->nodes_list[i],this->constraints_list[i],this->parameters);
         NodesDofx(this->nodes_list[i],this->parameters);
         ConfigureObjectsx(this->elements,this->loads_list[i],this->nodes_list[i],this->vertices,this->materials,this->parameters,this->inputs2);
      }
 
      //Reset current configuration:
      analysis_type=this->analysis_type_list[analysis_counter];
      SetCurrentConfiguration(analysis_type);
   }
}
/*}}}*/
void FemModel::Restart(){ /*{{{*/
 
   FILE* restartfid=NULL;
   char* restartfilename = NULL;
   int   femmodel_size=0;
   int   fread_return=0;
   char* femmodel_buffer=NULL;
   char* femmodel_buffer_ini=NULL;
 
   /*First, recover the name of the restart file: */
   parameters->FindParam(&restartfilename,RestartFileNameEnum);
 
   /*Now, figure out whether this file actually exists!: */
   restartfid=pfopen(restartfilename,"r",false);
 
   if(restartfid==NULL){
      xDelete<char>(restartfilename);
      return; //could not find the file, so no restart possible.
   }
 
   /*Print banner*/
   _printf0_("                                                                    \n");
   _printf0_("====================================================================\n");
   _printf0_("                         RESTART DETECTED                           \n");
   _printf0_("                                                                    \n");
   _printf0_("  Restart file: "<<restartfilename<<"                               \n");
   _printf0_("====================================================================\n");
   _printf0_("                                                                    \n");
 
   /*Figure out size of buffer to be read: */
   fseek(restartfid, 0L, SEEK_END);
   femmodel_size = ftell(restartfid);
   fseek(restartfid, 0L, SEEK_SET);
 
   /*Allocate buffer: */
   femmodel_buffer=xNew<char>(femmodel_size);
 
   /*Read buffer from file: */
   fread_return=fread(femmodel_buffer,femmodel_size,sizeof(char),restartfid); if(fread_return!=1)_error_("error reading the buffer from marshalled file!");
   femmodel_buffer_ini=femmodel_buffer; //keep track of the initial position, so as to free later.
 
   /*Create new FemModel by demarshalling the buffer: */
   this->Marshall(&femmodel_buffer,NULL,MARSHALLING_BACKWARD);
 
   /*Reset position of buffer: */
   femmodel_buffer=femmodel_buffer_ini;
 
   /*Done, close file :*/
   pfclose(restartfid,restartfilename);
 
   /*Free ressources: */
   xDelete<char>(restartfilename);
   xDelete<char>(femmodel_buffer);
}/*}}}*/
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 index = AnalysisIndex(configuration_type);
 
   /*If we already have the right analysis, return*/
   //if(this->analysis_counter==index) return;
   this->analysis_counter=index;
 
   /*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->loads = this->loads_list[this->analysis_counter];
   this->constraints = this->constraints_list[this->analysis_counter];
   this->nodes = this->nodes_list[this->analysis_counter];
   this->loads->SetCurrentConfiguration(elements, loads, nodes,vertices, materials,parameters);
   this->elements->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: " << EnumToStringx(analysis_type) << "\n");
   }
 
}/*}}}*/
void FemModel::SetCurrentConfiguration(int configuration_type){/*{{{*/
   this->SetCurrentConfiguration(configuration_type,configuration_type);
}
/*}}}*/
int  FemModel::Size(){ /*{{{*/
   int   femmodel_size;
 
   this->Marshall(NULL,&femmodel_size,MARSHALLING_SIZE);
 
   return femmodel_size;
}
/*}}}*/
void FemModel::SolutionAnalysesList(int** panalyses,int* pnumanalyses,IoModel* iomodel,int solutiontype){/*{{{*/
 
   /*output: */
   int  numanalyses = 0;
   int *analyses    = NULL;
 
   /*Intermediaries*/
   const int MAXANALYSES = 30;
   int   analyses_temp[MAXANALYSES];
 
   /*Analyses lists*/
   switch(solutiontype){
 
      case StressbalanceSolutionEnum:{
         bool isSIA,isFS;
         int  fe_FS;
         iomodel->FindConstant(&fe_FS,"md.flowequation.fe_FS");
         iomodel->FindConstant(&isSIA,"md.flowequation.isSIA");
         iomodel->FindConstant(&isFS,"md.flowequation.isFS");
         analyses_temp[numanalyses++]=StressbalanceAnalysisEnum;
         analyses_temp[numanalyses++]=StressbalanceVerticalAnalysisEnum;
         if(isSIA){
            analyses_temp[numanalyses++]=StressbalanceSIAAnalysisEnum;
         }
         analyses_temp[numanalyses++]=L2ProjectionBaseAnalysisEnum;
         analyses_temp[numanalyses++]=ExtrudeFromBaseAnalysisEnum;
         analyses_temp[numanalyses++]=DepthAverageAnalysisEnum;
         if(fe_FS==LATaylorHoodEnum || fe_FS==LACrouzeixRaviartEnum){
            analyses_temp[numanalyses++]=UzawaPressureAnalysisEnum;
         }
         }
         break;
 
      case SteadystateSolutionEnum:{
         bool isSIA,isenthalpy;
         iomodel->FindConstant(&isSIA,"md.flowequation.isSIA");
         iomodel->FindConstant(&isenthalpy,"md.thermal.isenthalpy");
         analyses_temp[numanalyses++]=StressbalanceAnalysisEnum;
         analyses_temp[numanalyses++]=StressbalanceVerticalAnalysisEnum;
         if(isSIA){
            analyses_temp[numanalyses++]=StressbalanceSIAAnalysisEnum;
         }
         if(isenthalpy){
            analyses_temp[numanalyses++]=EnthalpyAnalysisEnum;
         }
         else{
            analyses_temp[numanalyses++]=ThermalAnalysisEnum;
            analyses_temp[numanalyses++]=MeltingAnalysisEnum;
         }
         analyses_temp[numanalyses++]=L2ProjectionBaseAnalysisEnum;
         }
         break;
 
      case ThermalSolutionEnum:{
         bool isenthalpy;
         iomodel->FindConstant(&isenthalpy,"md.thermal.isenthalpy");
         if(isenthalpy){
            analyses_temp[numanalyses++]=EnthalpyAnalysisEnum;
         }
         else{
            analyses_temp[numanalyses++]=ThermalAnalysisEnum;
            analyses_temp[numanalyses++]=MeltingAnalysisEnum;
         }
         }
         break;
 
      case HydrologySolutionEnum:{
         int hydrology_model;
         iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
         if(hydrology_model==HydrologyshreveEnum){
            analyses_temp[numanalyses++]=HydrologyShreveAnalysisEnum;
         }
         else if(hydrology_model==HydrologyGlaDSEnum){
            analyses_temp[numanalyses++]=HydrologyGlaDSAnalysisEnum;
         }
         if(hydrology_model==HydrologyshaktiEnum){
            analyses_temp[numanalyses++]=HydrologyShaktiAnalysisEnum;
         }
         if(hydrology_model==HydrologypismEnum){
            analyses_temp[numanalyses++]=HydrologyPismAnalysisEnum;
         }
         if(hydrology_model==HydrologydcEnum){
            analyses_temp[numanalyses++]=HydrologyDCInefficientAnalysisEnum;
            analyses_temp[numanalyses++]=HydrologyDCEfficientAnalysisEnum;
            analyses_temp[numanalyses++]=L2ProjectionEPLAnalysisEnum;
            analyses_temp[numanalyses++]=L2ProjectionBaseAnalysisEnum;
         }
      }
         break;
 
      case MasstransportSolutionEnum:
         analyses_temp[numanalyses++]=DepthAverageAnalysisEnum;
         analyses_temp[numanalyses++]=SmbAnalysisEnum;
         analyses_temp[numanalyses++]=MasstransportAnalysisEnum;
         analyses_temp[numanalyses++]=ExtrudeFromBaseAnalysisEnum;
 
         break;
 
      case BalancethicknessSolutionEnum:
         analyses_temp[numanalyses++]=BalancethicknessAnalysisEnum;
         break;
 
      case Balancethickness2SolutionEnum:
         analyses_temp[numanalyses++]=L2ProjectionBaseAnalysisEnum;
         analyses_temp[numanalyses++]=SmoothAnalysisEnum;
         analyses_temp[numanalyses++]=Balancethickness2AnalysisEnum;
         break;
 
      case BalancethicknessSoftSolutionEnum:
         analyses_temp[numanalyses++]=BalancethicknessAnalysisEnum;
         break;
 
      case BalancevelocitySolutionEnum:
         analyses_temp[numanalyses++]=BalancevelocityAnalysisEnum;
         analyses_temp[numanalyses++]=SmoothAnalysisEnum;
         break;
 
      case SurfaceSlopeSolutionEnum:
         analyses_temp[numanalyses++]=L2ProjectionBaseAnalysisEnum;
         break;
 
      case BedSlopeSolutionEnum:
         analyses_temp[numanalyses++]=L2ProjectionBaseAnalysisEnum;
         break;
 
      case GiaSolutionEnum:
         analyses_temp[numanalyses++]=GiaAnalysisEnum;
         break;
 
      case LoveSolutionEnum:
         analyses_temp[numanalyses++]=LoveAnalysisEnum;
         break;
 
      case EsaSolutionEnum:
         analyses_temp[numanalyses++]=EsaAnalysisEnum;
         break;
 
      case SealevelriseSolutionEnum:
         analyses_temp[numanalyses++]=SealevelriseAnalysisEnum;
         analyses_temp[numanalyses++]=GiaAnalysisEnum;
         break;
 
      case SmbSolutionEnum:
         analyses_temp[numanalyses++]=SmbAnalysisEnum;
         break;
 
      case DamageEvolutionSolutionEnum:
         analyses_temp[numanalyses++]=DamageEvolutionAnalysisEnum;
         break;
 
      case TransientSolutionEnum:{
         /*We have multiple analyses here, process one by one*/
         bool isSIA,isFS,isthermal,isenthalpy,ismasstransport,isgroundingline,isstressbalance,ismovingfront,ishydrology,isdamage,issmb,isslr,isesa,isgia;
         iomodel->FindConstant(&isthermal,"md.transient.isthermal");
         iomodel->FindConstant(&ismovingfront,"md.transient.ismovingfront");
         iomodel->FindConstant(&ismasstransport,"md.transient.ismasstransport");
         iomodel->FindConstant(&isstressbalance,"md.transient.isstressbalance");
         iomodel->FindConstant(&isgroundingline,"md.transient.isgroundingline");
         iomodel->FindConstant(&isdamage,"md.transient.isdamageevolution");
         iomodel->FindConstant(&ishydrology,"md.transient.ishydrology");
         iomodel->FindConstant(&issmb,"md.transient.issmb");
         iomodel->FindConstant(&isslr,"md.transient.isslr");
         iomodel->FindConstant(&isesa,"md.transient.isesa");
         iomodel->FindConstant(&isgia,"md.transient.isgia");
         int* analyses_iter     = NULL;
         int  num_analyses_iter = 0;
         if(isstressbalance){
            SolutionAnalysesList(&analyses_iter,&num_analyses_iter,iomodel,StressbalanceSolutionEnum);
            xMemCpy<int>(&analyses_temp[numanalyses],analyses_iter,num_analyses_iter);
            numanalyses+=num_analyses_iter; xDelete<int>(analyses_iter);
         }
         if(isthermal && iomodel->domaintype==Domain3DEnum){
            SolutionAnalysesList(&analyses_iter,&num_analyses_iter,iomodel,ThermalSolutionEnum);
            xMemCpy<int>(&analyses_temp[numanalyses],analyses_iter,num_analyses_iter);
            numanalyses+=num_analyses_iter; xDelete<int>(analyses_iter);
         }
         if(ismasstransport || isgroundingline){
            analyses_temp[numanalyses++]=MasstransportAnalysisEnum;
            int  basalforcing_model;
            iomodel->FindConstant(&basalforcing_model,"md.basalforcings.model");
            if(basalforcing_model==BasalforcingsPicoEnum){
               bool isplume;
               iomodel->FindConstant(&isplume,"md.basalforcings.isplume");
               if(isplume){
                  analyses_temp[numanalyses++]=GLheightadvectionAnalysisEnum;
               }
            }
         }
         if(issmb) analyses_temp[numanalyses++]=SmbAnalysisEnum;
         if(ismovingfront){
            analyses_temp[numanalyses++]=ExtrapolationAnalysisEnum;
            analyses_temp[numanalyses++]=LevelsetAnalysisEnum;
         }
         if(ishydrology){
            SolutionAnalysesList(&analyses_iter,&num_analyses_iter,iomodel,HydrologySolutionEnum);
            xMemCpy<int>(&analyses_temp[numanalyses],analyses_iter,num_analyses_iter);
            numanalyses+=num_analyses_iter; xDelete<int>(analyses_iter);
         }
         if(isdamage){
            analyses_temp[numanalyses++]=DamageEvolutionAnalysisEnum;
         }
         if(isslr){
            analyses_temp[numanalyses++]=GiaAnalysisEnum;
            analyses_temp[numanalyses++]=SealevelriseAnalysisEnum;
         }
         if(isesa){
            analyses_temp[numanalyses++]=EsaAnalysisEnum;
         }
         if(isgia){
            analyses_temp[numanalyses++]=GiaAnalysisEnum;
         }
 
         if(iomodel->domaintype==Domain2DverticalEnum || iomodel->domaintype==Domain3DEnum){
            analyses_temp[numanalyses++]=ExtrudeFromBaseAnalysisEnum;
            analyses_temp[numanalyses++]=ExtrudeFromTopAnalysisEnum;
            analyses_temp[numanalyses++]=FreeSurfaceBaseAnalysisEnum;
            analyses_temp[numanalyses++]=FreeSurfaceTopAnalysisEnum;
         }
         analyses_temp[numanalyses++]=L2ProjectionBaseAnalysisEnum;
         }
         break;
 
      default:
         _error_("solution type: " << EnumToStringx(solutiontype) << " not supported yet!");
         break;
   }
 
   /*Copy analyses from temp to output*/
   _assert_(numanalyses<MAXANALYSES);
   analyses=xNew<int>(numanalyses);
   for(int i=0;i<numanalyses;i++) analyses[i]=analyses_temp[i];
 
   /*Assign output pointers:*/
   if(pnumanalyses) *pnumanalyses=numanalyses;
   if(panalyses)    *panalyses=analyses;
   else              xDelete<int>(analyses);
}/*}}}*/
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 solution core: */
   profiler->Start(CORE);
   solutioncore(this);
   profiler->Stop(CORE);
 
   /*run AD core if needed: */
   profiler->Start(ADCORE);
   ad_core(this);
   profiler->Stop(ADCORE);
 
   /*some profiling results for the core: */
   parameters->FindParam(&profiling,DebugProfilingEnum);
   if(profiling){
 
      solution_time=profiler->TotalTime(CORE);
      solution_flops=profiler->TotalFlops(CORE);
      solution_memory=profiler->Memory(CORE);
 
      _printf0_("\nCore solution profiling\n");
      _printf0_("   elapsed time    : " << solution_time   << " Seconds\n");
      _printf0_("   number of flops : " << solution_flops  << " Flops\n");
      _printf0_("   memory used     : " << solution_memory << " Bytes\n");
 
      /*Individual cores*/
      _printf0_("Individual core profiling\n");
      if(profiler->Used(THERMALCORE)) _printf0_("   "<<setw(40)<<left<<"Thermal core elapsed time:"<<setw(7)<<setprecision(6)<<profiler->TotalTime(THERMALCORE) << " sec\n");
      if(profiler->Used(HYDROLOGYCORE)) _printf0_("   "<<setw(40)<<left<<"Hydrology core elapsed time:"<<setw(7)<<setprecision(6)<<profiler->TotalTime(HYDROLOGYCORE) << " sec\n");
      if(profiler->Used(STRESSBALANCECORE)) _printf0_("   "<<setw(40)<<left<<"Stress balance core elapsed time:"<<setw(7)<<setprecision(6)<<profiler->TotalTime(STRESSBALANCECORE) << " sec\n");
      if(profiler->Used(DAMAGECORE)) _printf0_("   "<<setw(40)<<left<<"Damage core elapsed time:"<<setw(7)<<setprecision(6)<<profiler->TotalTime(DAMAGECORE) << " sec\n");
      if(profiler->Used(MOVINGFRONTCORE)) _printf0_("   "<<setw(40)<<left<<"Moving front core elapsed time:"<<setw(7)<<setprecision(6)<<profiler->TotalTime(MOVINGFRONTCORE) << " sec\n");
      if(profiler->Used(MASSTRANSPORTCORE)) _printf0_("   "<<setw(40)<<left<<"Mass transport core elapsed time:"<<setw(7)<<setprecision(6)<<profiler->TotalTime(MASSTRANSPORTCORE) << " sec\n");
      if(profiler->Used(SMBCORE)) _printf0_("   "<<setw(40)<<left<<"SMB core elapsed time:"<<setw(7)<<setprecision(6)<<profiler->TotalTime(SMBCORE) << " sec\n");
      if(profiler->Used(GROUNDINGLINECORE)) _printf0_("   "<<setw(40)<<left<<"Groundingline migration core elapsed time:"<<setw(7)<<setprecision(6)<<profiler->TotalTime(GROUNDINGLINECORE) << " sec\n");
      if(profiler->Used(GIACORE)) _printf0_("   "<<setw(40)<<left<<"GIA core elapsed time:"<<setw(7)<<setprecision(6)<<profiler->TotalTime(GIACORE) << " sec\n");
      if(profiler->Used(ESACORE)) _printf0_("   "<<setw(40)<<left<<"ESA core elapsed time:"<<setw(7)<<setprecision(6)<<profiler->TotalTime(ESACORE) << " sec\n");
      if(profiler->Used(SLRCORE)) _printf0_("   "<<setw(40)<<left<<"SLR core elapsed time:"<<setw(7)<<setprecision(6)<<profiler->TotalTime(SLRCORE) << " sec\n");
      if(profiler->Used(MPISERIAL)) _printf0_("   "<<setw(40)<<left<<"MPISERIAL elapsed time:"<<setw(7)<<setprecision(6)<<profiler->TotalTime(MPISERIAL) << " sec\n");
 
      /*Add to results: */
      results->AddObject(new GenericExternalResult<IssmDouble>(results->Size()+1, ProfilingSolutionTimeEnum,  solution_time));
      results->AddObject(new GenericExternalResult<IssmDouble>(results->Size()+1, ProfilingCurrentMemEnum,  solution_memory));
      results->AddObject(new GenericExternalResult<IssmDouble>(results->Size()+1, ProfilingCurrentFlopsEnum, solution_flops));
 
      #ifdef _HAVE_AD_
      solution_time   = profiler->TotalTime(ADCORE);
      solution_flops  = profiler->TotalFlops(ADCORE);
      solution_memory = profiler->Memory(ADCORE);
 
      _printf0_("AD profiling\n");
      _printf0_("   elapsed time    : " << solution_time   << " Seconds\n");
      _printf0_("   number of flops : " << solution_flops  << " Flops\n");
      _printf0_("   memory used     : " << solution_memory << " Bytes\n");
      #endif
      _printf0_("\n");
   }
}
/*}}}*/
 
/*Modules:*/
void FemModel::BalancethicknessMisfitx(IssmDouble* presponse){/*{{{*/
 
   /*output: */
   IssmDouble J=0.;
   IssmDouble J_sum;
 
   IssmDouble  weight,vx,vy,H,dvx[2],dvy[2],dH[2];
   IssmDouble  temp,Jdet,dhdt,groundedice_melting,surface_mass_balance;
   IssmDouble* xyz_list = NULL;
   IssmDouble  dp[3];
 
   /*Compute Misfit: */
   for(int i=0;i<elements->Size();i++){
      Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
 
      /*If on water, return 0: */
      if(!element->IsIceInElement()) continue;
 
      /* Get node coordinates*/
      element->GetVerticesCoordinates(&xyz_list);
      DatasetInput2* weights_input                   = element->GetDatasetInput2(InversionCostFunctionsCoefficientsEnum);   _assert_(weights_input);
      Input2* thickness_input                 = element->GetInput2(ThicknessEnum); _assert_(thickness_input);
      Input2* vx_input                        = element->GetInput2(VxEnum);                                  _assert_(vx_input);
      Input2* vy_input                        = element->GetInput2(VyEnum);                                  _assert_(vy_input);
      Input2* surface_mass_balance_input      = element->GetInput2(SmbMassBalanceEnum);          _assert_(surface_mass_balance_input);
      Input2* groundedice_melting_input       = element->GetInput2(BasalforcingsGroundediceMeltingRateEnum); _assert_(groundedice_melting_input);
      Input2* dhdt_input                      = element->GetInput2(BalancethicknessThickeningRateEnum);      _assert_(dhdt_input);
 
      /* Start  looping on the number of gaussian points: */
      Gauss* gauss=element->NewGauss(2);
      for(int ig=gauss->begin();ig<gauss->end();ig++){
 
         gauss->GaussPoint(ig);
 
         /* Get Jacobian determinant: */
         element->JacobianDeterminant(&Jdet,xyz_list,gauss);
 
         /*Get all parameters at gaussian point*/
         weights_input->GetInputValue(&weight,gauss,BalancethicknessMisfitEnum);
         thickness_input->GetInputValue(&H, gauss);
         thickness_input->GetInputDerivativeValue(&dH[0],xyz_list,gauss);
         surface_mass_balance_input->GetInputValue(&surface_mass_balance,gauss);
         groundedice_melting_input->GetInputValue(&groundedice_melting,gauss);
         dhdt_input->GetInputValue(&dhdt,gauss);
         vx_input->GetInputValue(&vx,gauss);
         vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
         vy_input->GetInputValue(&vy,gauss);
         vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
 
         /*Balance thickness soft constraint J = 1/2 (div(Hv)-a)^2*/
         temp  = vx*dH[0]+vy*dH[1]+H*(dvx[0]+dvy[1]) - (surface_mass_balance-groundedice_melting-dhdt);
         J    +=weight*1/2*temp*temp*Jdet*gauss->weight;
      }
 
      /*clean up and Return: */
      xDelete<IssmDouble>(xyz_list);
      delete gauss;
   }
 
   /*Sum all J from all cpus of the cluster:*/
   ISSM_MPI_Reduce (&J,&J_sum,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
   ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
   J=J_sum;
 
   /*Assign output pointers: */
   *presponse=J;
 
}/*}}}*/
void FemModel::CalvingRateVonmisesx(){/*{{{*/
 
   for(int i=0;i<elements->Size();i++){
      Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
      element->CalvingRateVonmises();
   }
}
/*}}}*/
void FemModel::CalvingRateLevermannx(){/*{{{*/
 
   for(int i=0;i<elements->Size();i++){
      Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
      element->CalvingRateLevermann();
   }
}
/*}}}*/
void FemModel::CalvingFluxLevelsetx(){/*{{{*/
 
   for(int i=0;i<elements->Size();i++){
      Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
      element->CalvingFluxLevelset();
   }
}
/*}}}*/
void FemModel::CalvingMeltingFluxLevelsetx(){/*{{{*/
 
   for(int i=0;i<elements->Size();i++){
      Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
      element->CalvingMeltingFluxLevelset();
   }
}
/*}}}*/
void FemModel::CostFunctionx(IssmDouble* pJ,IssmDouble** pJlist,int* pn){/*{{{*/
 
   /*Intermediary*/
   int      num_responses;
   int     *responses      = NULL;
   Results *cost_functions = NULL;
 
   /*Recover parameters*/
   parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
   parameters->FindParam(&responses,NULL,InversionCostFunctionsEnum);
 
   /*Get the value of all cost functions*/
   this->RequestedOutputsx(&cost_functions,responses,num_responses);
 
   /*Get and add all contributions one by one*/
   IssmDouble  J=0;
   IssmDouble* Jlist = xNew<IssmDouble>(num_responses);
   for(int i=0;i<num_responses;i++){
      ExternalResult* result=(ExternalResult*)cost_functions->GetObjectByOffset(i);
      Jlist[i] = reCast<IssmDouble>(result->GetValue());
      J       += Jlist[i];
   }
   _assert_(cost_functions->Size()==num_responses);
 
   /*Assign output pointers: */
   delete cost_functions;
   xDelete<int>(responses);
   if(pJ)     *pJ     = J;
   if(pJlist) *pJlist = Jlist;
   else        xDelete<IssmDouble>(Jlist);
   if(pn)     *pn     = num_responses;
}
/*}}}*/
void FemModel::DeviatoricStressx(){/*{{{*/
 
   for(int i=0;i<elements->Size();i++){
      Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
      element->ComputeDeviatoricStressTensor();
   }
}
/*}}}*/
void FemModel::DistanceToFieldValue(int fieldenum,IssmDouble fieldvalue,int distanceenum){/*{{{*/
 
   /*recover my_rank:*/
   int my_rank   = IssmComm::GetRank();
   int num_procs = IssmComm::GetSize();
 
   /*Get domain type (2d or 3d)*/
   int domaintype;
   this->parameters->FindParam(&domaintype,DomainTypeEnum);
 
   /*1: go throug all elements of this partition and figure out how many
    * segments we have (corresopnding to field = value)*/
   DataSet* segments=new DataSet();
   for(int i=0;i<elements->Size();i++){
      Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
      if(!element->IsOnBase()) continue;
      Element* basalelement = element->SpawnBasalElement();
      basalelement->WriteFieldIsovalueSegment(segments,fieldenum,fieldvalue);
      if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
   }
 
   /*2: now get the segments from all partitions*/
   int  segcount=segments->Size();
   int* allsegcount=xNew<int>(num_procs);
   ISSM_MPI_Gather(&segcount,1,ISSM_MPI_INT,allsegcount,1,ISSM_MPI_INT,0,IssmComm::GetComm());
   ISSM_MPI_Bcast(allsegcount,num_procs,ISSM_MPI_INT,0,IssmComm::GetComm());
 
   /* Every cpu should start its own dof count at the end of the dofcount from cpu-1*/
   int numseg_offset=0;
   int numseg=0;
   for(int i=0;i<my_rank;  i++) numseg_offset+=allsegcount[i];
   for(int i=0;i<num_procs;i++) numseg+=allsegcount[i];
   IssmDouble* segmentlist    = xNewZeroInit<IssmDouble>(4*numseg);
   IssmDouble* allsegmentlist = xNewZeroInit<IssmDouble>(4*numseg);
   for(int i=0;i<segments->Size();i++){
      Contour<IssmDouble>* segment=(Contour<IssmDouble>*)segments->GetObjectByOffset(i);
      _assert_(segment->nods == 2);
      segmentlist[(numseg_offset+i)*4 + 0] = segment->x[0];
      segmentlist[(numseg_offset+i)*4 + 1] = segment->y[0];
      segmentlist[(numseg_offset+i)*4 + 2] = segment->x[1];
      segmentlist[(numseg_offset+i)*4 + 3] = segment->y[1];
   }
 
   ISSM_MPI_Allreduce((void*)segmentlist,(void*)allsegmentlist,4*numseg,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
   delete segments;
   xDelete<IssmDouble>(segmentlist);
   xDelete<int>(allsegcount);
 
   /*3: Add distance input to all elements*/
   IssmDouble* distances = xNew<IssmDouble>(vertices->Size());
   IssmDouble  d,xn,yn,dmin;
   int         last = -1;
   for(int v=0;v<vertices->Size();v++){
 
      Vertex* vertex=dynamic_cast<Vertex*>(this->vertices->GetObjectByOffset(v));
      IssmDouble x = vertex->x;
      IssmDouble y = vertex->y;
 
      /*Most of the time the last checked segment is the closest so start with that one*/
      if(last>0){
         dmin = pow(allsegmentlist[4*last+0] - x,2) + pow(y-allsegmentlist[4*last+1],2);
      }
      else{
         dmin = 1.e+50;
      }
 
      for(int i=0;i<numseg;i++){
 
         /*Skip if tip is more than 10xdmin away*/
         if(pow(allsegmentlist[4*i+0] - x,2)>10*dmin) continue;
         if(pow(allsegmentlist[4*i+1] - y,2)>10*dmin) continue;
 
         IssmDouble l2 = (allsegmentlist[4*i+2]-allsegmentlist[4*i+0])*(allsegmentlist[4*i+2]-allsegmentlist[4*i+0]) + (allsegmentlist[4*i+3]-allsegmentlist[4*i+1])*(allsegmentlist[4*i+3]-allsegmentlist[4*i+1]);
 
         /*Segment has a length of 0*/
         if(l2==0.){
            d = (x-allsegmentlist[4*i+0])*(x-allsegmentlist[4*i+0])+(y-allsegmentlist[4*i+1])*(y-allsegmentlist[4*i+1]);
            if(d<dmin){
               dmin = d;
               last = i;
            }
            continue;
         }
 
         /*Consider the line extending the segment, parameterized as v + t (w - v).
          *We find projection of point p onto the line.
          *It falls where t = [(p-v) . (w-v)] / |w-v|^2*/
         IssmDouble t = ((x-allsegmentlist[4*i+0])*(allsegmentlist[4*i+2]-allsegmentlist[4*i+0]) + (y-allsegmentlist[4*i+1])*(allsegmentlist[4*i+3]-allsegmentlist[4*i+1]))/l2;
         if(t < 0.0){
            // Beyond the 'v' end of the segment
            d = (x-allsegmentlist[4*i+0])*(x-allsegmentlist[4*i+0])+(y-allsegmentlist[4*i+1])*(y-allsegmentlist[4*i+1]);
         }
         else if (t > 1.0){
            // Beyond the 'w' end of the segment
            d = (x-allsegmentlist[4*i+2])*(x-allsegmentlist[4*i+2])+(y-allsegmentlist[4*i+3])*(y-allsegmentlist[4*i+3]);
         }
         else{
            // Projection falls on the segment
            xn = allsegmentlist[4*i+0] + t * (allsegmentlist[4*i+2] - allsegmentlist[4*i+0]);
            yn = allsegmentlist[4*i+1] + t * (allsegmentlist[4*i+3] - allsegmentlist[4*i+1]);
            d = (x-xn)*(x-xn)+(y-yn)*(y-yn);
         }
 
         if(d<dmin){
            dmin = d;
            last = i;
         }
      }
 
      /*Update signed distance*/
      _assert_(vertex->lid<vertices->Size());
      distances[vertex->lid] = sqrt(dmin);
   }
 
   for(int i=0;i<this->elements->Size();i++){
      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
      element->CreateDistanceInputFromSegmentlist(distances,distanceenum);
   }
   //InputUpdateFromVectorx(this,distances,distanceenum,VertexLIdEnum);
 
   /*Clean up and return*/
   xDelete<IssmDouble>(distances);
   xDelete<IssmDouble>(allsegmentlist);
}/*}}}*/
void FemModel::Divergencex(IssmDouble* pdiv){/*{{{*/
 
   IssmDouble local_divergence=0;
   IssmDouble total_divergence;
 
   for(int i=0;i<this->elements->Size();i++){
      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
      local_divergence+=element->Divergence();
   }
   ISSM_MPI_Reduce(&local_divergence,&total_divergence,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
   ISSM_MPI_Bcast(&total_divergence,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
 
   /*Assign output pointers: */
   *pdiv=total_divergence;
 
}/*}}}*/
void FemModel::ElementOperationx(void (Element::*function)(void)){ /*{{{*/
 
   for(int i=0;i<elements->Size();i++){
      Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
      (element->*function)();
   }
 
}
/*}}}*/
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=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
      if (element->Id()==index){
         found=1;
         cpu_found=my_rank;
         break;
      }
   }
 
   /*Broadcast whether we found the element: */
   ISSM_MPI_Allreduce ( &found,&sumfound,1,ISSM_MPI_INT,ISSM_MPI_SUM,IssmComm::GetComm());
   if(!sumfound)_error_("could not find material with id" << index << " to compute ElementResponse");
 
   /*Ok, we found the element, compute responseocity: */
   if(my_rank==cpu_found){
      element->ElementResponse(&response,response_enum);
   }
 
   /*Broadcast and plug into response: */
   ISSM_MPI_Allreduce ( &cpu_found,&cpu_found,1,ISSM_MPI_INT,ISSM_MPI_MAX,IssmComm::GetComm());
   ISSM_MPI_Bcast(&response,1,ISSM_MPI_DOUBLE,cpu_found,IssmComm::GetComm());
 
   /*Assign output pointers: */
   *presponse=response;
 
}/*}}}*/
void FemModel::FloatingAreax(IssmDouble* pV, bool scaled){/*{{{*/
 
   IssmDouble local_floating_area= 0;
   IssmDouble total_floating_area;
 
   for(int i=0;i<this->elements->Size();i++){
      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
      local_floating_area+=element->FloatingArea(scaled);
   }
   ISSM_MPI_Reduce(&local_floating_area,&total_floating_area,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
   ISSM_MPI_Bcast(&total_floating_area,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
 
   /*Assign output pointers: */
   *pV=total_floating_area;
 
}/*}}}*/
void FemModel::GetInputLocalMinMaxOnNodesx(IssmDouble** pmin,IssmDouble** pmax,IssmDouble* ug){/*{{{*/
 
   /*Initialize output vectors*/
   int numnodes = this->nodes->NumberOfNodes();
   IssmDouble* uLmin_local = xNew<IssmDouble>(numnodes);
   IssmDouble* uLmax_local = xNew<IssmDouble>(numnodes);
   IssmDouble* uLmin = xNew<IssmDouble>(numnodes);
   IssmDouble* uLmax = xNew<IssmDouble>(numnodes);
   for(int i=0;i<numnodes;i++){
      uLmin_local[i] = +1.e+50;
      uLmax_local[i] = -1.e+50;
   }
 
   for(int i=0;i<this->elements->Size();i++){
      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
      element->GetInputLocalMinMaxOnNodes(uLmin_local,uLmax_local,ug);
   }
 
   /*Synchronize all CPUs*/
   ISSM_MPI_Allreduce((void*)uLmin_local,(void*)uLmin,numnodes,ISSM_MPI_DOUBLE,ISSM_MPI_MIN,IssmComm::GetComm());
   ISSM_MPI_Allreduce((void*)uLmax_local,(void*)uLmax,numnodes,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
   xDelete<IssmDouble>(uLmin_local);
   xDelete<IssmDouble>(uLmax_local);
 
   /*Assign output pointers: */
   *pmin=uLmin;
   *pmax=uLmax;
 
}/*}}}*/
void FemModel::GetLocalVectorWithClonesGset(IssmDouble** plocal_ug,Vector<IssmDouble> *ug){/*{{{*/
 
   this->nodes->GetLocalVectorWithClonesGset(plocal_ug,ug);
 
}/*}}}*/
void FemModel::GetLocalVectorWithClonesVertices(IssmDouble** plocal_vector,Vector<IssmDouble> *vector){/*{{{*/
 
   /*retrieve vertex info*/
   int localsize         = this->vertices->NumberOfVerticesLocalAll();
   int localsize_masters = this->vertices->NumberOfVerticesLocal();
 
   /*Get local vector of vector*/
   int        *indices_vector_masters = NULL;
   IssmDouble *local_vector_masters   = NULL;
   vector->GetLocalVector(&local_vector_masters,&indices_vector_masters);
   _assert_(localsize_masters==indices_vector_masters[localsize_masters-1] - indices_vector_masters[0]+1);
   xDelete<int>(indices_vector_masters);
 
   /*Now, extend vectors to account for clones (make vectors longer, for clones at the end)*/
   IssmDouble *local_vector  = xNew<IssmDouble>(localsize);
   xMemCpy<IssmDouble>(local_vector,local_vector_masters,localsize_masters);
   xDelete<IssmDouble>(local_vector_masters);
 
   /*Now send and receive vector for vertices on partition edge*/
   SyncLocalVectorWithClonesVertices(local_vector);
 
   /*Assign output pointer*/
   *plocal_vector = local_vector;
}/*}}}*/
void FemModel::SyncLocalVectorWithClonesVertices(IssmDouble* local_vector){/*{{{*/
 
   /*recover my_rank:*/
   ISSM_MPI_Status status;
   int my_rank   = IssmComm::GetRank();
   int num_procs = IssmComm::GetSize();
 
   /*Now send and receive vector for vertices on partition edge*/
   #ifdef _HAVE_AD_
   IssmDouble* buffer = xNew<IssmDouble>(this->vertices->Size(),"t"); //only one alloc, "t" is required by adolc
   #else
   IssmDouble* buffer = xNew<IssmDouble>(this->vertices->Size());
   #endif
   for(int rank=0;rank<num_procs;rank++){
      if(this->vertices->common_send[rank]){
         int  numids = this->vertices->common_send[rank];
         for(int i=0;i<numids;i++){
            int   master_lid = this->vertices->common_send_ids[rank][i];
            Vertex* vertex=xDynamicCast<Vertex*>(this->vertices->GetObjectByOffset(master_lid));
            _assert_(!vertex->clone);
            buffer[i] = local_vector[vertex->lid];
         }
         ISSM_MPI_Send(buffer,numids,ISSM_MPI_DOUBLE,rank,0,IssmComm::GetComm());
      }
   }
   for(int rank=0;rank<num_procs;rank++){
      if(this->vertices->common_recv[rank]){
         int  numids = this->vertices->common_recv[rank];
         ISSM_MPI_Recv(buffer,numids,ISSM_MPI_DOUBLE,rank,0,IssmComm::GetComm(),&status);
         for(int i=0;i<numids;i++){
            int   master_lid = this->vertices->common_recv_ids[rank][i];
            Vertex* vertex=xDynamicCast<Vertex*>(this->vertices->GetObjectByOffset(master_lid));
            _assert_(vertex->clone);
            local_vector[vertex->lid] = buffer[i];
         }
      }
   }
   xDelete<IssmDouble>(buffer);
}/*}}}*/
void FemModel::SyncLocalVectorWithClonesVerticesAdd(IssmDouble* local_vector){/*{{{*/
 
   /*recover my_rank:*/
   ISSM_MPI_Status status;
   int my_rank   = IssmComm::GetRank();
   int num_procs = IssmComm::GetSize();
 
   /*Now send and receive vector for vertices on partition edge*/
   #ifdef _HAVE_AD_
   IssmDouble* buffer = xNew<IssmDouble>(this->vertices->Size(),"t"); //only one alloc, "t" is required by adolc
   #else
   IssmDouble* buffer = xNew<IssmDouble>(this->vertices->Size());
   #endif
 
   /*1st: add slaves to master values (reverse of what we usually do)*/
   for(int rank=0;rank<num_procs;rank++){
      if(this->vertices->common_recv[rank]){
         int  numids = this->vertices->common_recv[rank];
         for(int i=0;i<numids;i++){
            int   master_lid = this->vertices->common_recv_ids[rank][i];
            Vertex* vertex=xDynamicCast<Vertex*>(this->vertices->GetObjectByOffset(master_lid));
            _assert_(vertex->clone);
            buffer[i] = local_vector[vertex->lid];
         }
         ISSM_MPI_Send(buffer,numids,ISSM_MPI_DOUBLE,rank,0,IssmComm::GetComm());
      }
   }
   for(int rank=0;rank<num_procs;rank++){
      if(this->vertices->common_send[rank]){
         int  numids = this->vertices->common_send[rank];
         ISSM_MPI_Recv(buffer,numids,ISSM_MPI_DOUBLE,rank,0,IssmComm::GetComm(),&status);
         for(int i=0;i<numids;i++){
            int   master_lid = this->vertices->common_send_ids[rank][i];
            Vertex* vertex=xDynamicCast<Vertex*>(this->vertices->GetObjectByOffset(master_lid));
            _assert_(!vertex->clone);
            local_vector[vertex->lid] += buffer[i];
         }
      }
   }
 
   /*Now sync masters across partitions*/
   for(int rank=0;rank<num_procs;rank++){
      if(this->vertices->common_send[rank]){
         int  numids = this->vertices->common_send[rank];
         for(int i=0;i<numids;i++){
            int   master_lid = this->vertices->common_send_ids[rank][i];
            Vertex* vertex=xDynamicCast<Vertex*>(this->vertices->GetObjectByOffset(master_lid));
            _assert_(!vertex->clone);
            buffer[i] = local_vector[vertex->lid];
         }
         ISSM_MPI_Send(buffer,numids,ISSM_MPI_DOUBLE,rank,0,IssmComm::GetComm());
      }
   }
   for(int rank=0;rank<num_procs;rank++){
      if(this->vertices->common_recv[rank]){
         int  numids = this->vertices->common_recv[rank];
         ISSM_MPI_Recv(buffer,numids,ISSM_MPI_DOUBLE,rank,0,IssmComm::GetComm(),&status);
         for(int i=0;i<numids;i++){
            int   master_lid = this->vertices->common_recv_ids[rank][i];
            Vertex* vertex=xDynamicCast<Vertex*>(this->vertices->GetObjectByOffset(master_lid));
            _assert_(vertex->clone);
            local_vector[vertex->lid] = buffer[i];
         }
      }
   }
   xDelete<IssmDouble>(buffer);
}/*}}}*/
void FemModel::GetLocalVectorWithClonesNodes(IssmDouble** plocal_vector,Vector<IssmDouble> *vector){/*{{{*/
 
   /*recover my_rank:*/
   ISSM_MPI_Status status;
   int my_rank   = IssmComm::GetRank();
   int num_procs = IssmComm::GetSize();
 
   /*retrieve vertex info*/
   int localsize         = this->nodes->NumberOfNodesLocalAll();
   int localsize_masters = this->nodes->NumberOfNodesLocal();
 
   /*Get local vector of vector*/
   int        *indices_vector_masters = NULL;
   IssmDouble *local_vector_masters   = NULL;
   vector->GetLocalVector(&local_vector_masters,&indices_vector_masters);
   _assert_(localsize_masters==indices_vector_masters[localsize_masters-1] - indices_vector_masters[0]+1);
   xDelete<int>(indices_vector_masters);
 
   /*Now, extend vectors to account for clones (make vectors longer, for clones at the end)*/
   IssmDouble *local_vector  = xNew<IssmDouble>(localsize);
   xMemCpy<IssmDouble>(local_vector,local_vector_masters,localsize_masters);
   xDelete<IssmDouble>(local_vector_masters);
 
   /*Now send and receive vector for nodes on partition edge*/
   #ifdef _HAVE_AD_
   IssmDouble* buffer = xNew<IssmDouble>(this->nodes->Size(),"t"); //only one alloc, "t" is required by adolc
   #else
   IssmDouble* buffer = xNew<IssmDouble>(this->nodes->Size());
   #endif
   for(int rank=0;rank<num_procs;rank++){
      if(this->nodes->common_send[rank]){
         int  numids = this->nodes->common_send[rank];
         for(int i=0;i<numids;i++){
            int   master_lid = this->nodes->common_send_ids[rank][i];
            Node* vertex=xDynamicCast<Node*>(this->nodes->GetObjectByOffset(master_lid));
            _assert_(!vertex->clone);
            buffer[i] = local_vector[vertex->lid];
         }
         ISSM_MPI_Send(buffer,numids,ISSM_MPI_DOUBLE,rank,0,IssmComm::GetComm());
      }
   }
   for(int rank=0;rank<num_procs;rank++){
      if(this->nodes->common_recv[rank]){
         int  numids = this->nodes->common_recv[rank];
         ISSM_MPI_Recv(buffer,numids,ISSM_MPI_DOUBLE,rank,0,IssmComm::GetComm(),&status);
         for(int i=0;i<numids;i++){
            int   master_lid = this->nodes->common_recv_ids[rank][i];
            Node* vertex=xDynamicCast<Node*>(this->nodes->GetObjectByOffset(master_lid));
            _assert_(vertex->clone);
            local_vector[vertex->lid] = buffer[i];
         }
      }
   }
   xDelete<IssmDouble>(buffer);
 
   /*Assign output pointer*/
   *plocal_vector = local_vector;
}/*}}}*/
void FemModel::GroundedAreax(IssmDouble* pV, bool scaled){/*{{{*/
 
   IssmDouble local_grounded_area= 0;
   IssmDouble total_grounded_area;
 
   for(int i=0;i<this->elements->Size();i++){
      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
      local_grounded_area+=element->GroundedArea(scaled);
   }
   ISSM_MPI_Reduce(&local_grounded_area,&total_grounded_area,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
   ISSM_MPI_Bcast(&total_grounded_area,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
 
   /*Assign output pointers: */
   *pV=total_grounded_area;
 
}/*}}}*/
void FemModel::IcefrontAreax(){/*{{{*/
 
   int numvertices      = 6;
   int numbasins;
   IssmDouble* BasinId   = xNew<IssmDouble>(numvertices);
   this->parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum);
   IssmDouble* basin_icefront_area           = xNewZeroInit<IssmDouble>(numbasins);
 
   for(int basin=1;basin<numbasins+1;basin++){
      IssmDouble local_icefront_area = 0;
      IssmDouble total_icefront_area;
 
      for(int i=0;i<this->elements->Size();i++){
         Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
         element->GetInputListOnVertices(BasinId,FrontalForcingsBasinIdEnum);
         for(int j=0;j<numvertices;j++){
            if(BasinId[j]==basin){
               local_icefront_area+=element->GetIcefrontArea();
               break;
            }
         }
      }
      ISSM_MPI_Reduce(&local_icefront_area,&total_icefront_area,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
      ISSM_MPI_Bcast(&total_icefront_area,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
 
      basin_icefront_area[basin-1]=total_icefront_area;
   }
 
   this->parameters->AddObject(new DoubleVecParam(FrontalForcingsBasinIcefrontAreaEnum,basin_icefront_area,numbasins));
 
   xDelete<IssmDouble>(basin_icefront_area);
   xDelete<IssmDouble>(BasinId);
}/*}}}*/
void FemModel::IcefrontMassFluxx(IssmDouble* pM, bool scaled){/*{{{*/
 
   IssmDouble local_mass_flux = 0;
   IssmDouble total_mass_flux;
 
   for(int i=0;i<this->elements->Size();i++){
      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
      local_mass_flux+=element->IcefrontMassFlux(scaled);
   }
   ISSM_MPI_Reduce(&local_mass_flux,&total_mass_flux,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
   ISSM_MPI_Bcast(&total_mass_flux,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
 
   /*Assign output pointers: */
   *pM=total_mass_flux;
 
}/*}}}*/
void FemModel::IcefrontMassFluxLevelsetx(IssmDouble* pM, bool scaled){/*{{{*/
 
   IssmDouble local_mass_flux = 0;
   IssmDouble total_mass_flux;
 
   for(int i=0;i<this->elements->Size();i++){
      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
      local_mass_flux+=element->IcefrontMassFluxLevelset(scaled);
   }
   ISSM_MPI_Reduce(&local_mass_flux,&total_mass_flux,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
   ISSM_MPI_Bcast(&total_mass_flux,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
 
   /*Assign output pointers: */
   *pM=total_mass_flux;
 
}/*}}}*/
void FemModel::InputMakeDiscontinuous(int enum_in){/*{{{*/
 
   int numvertices  = 6;
   IssmDouble* P1DGlist = xNew<IssmDouble>(numvertices);
 
   for(int i=0;i<this->elements->Size();i++){
      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
      element->GetInputListOnVertices(P1DGlist,enum_in);
      element->AddInput2(DummyEnum,P1DGlist,P1DGEnum);
   }
   xDelete<IssmDouble>(P1DGlist);
 
   this->inputs2->ChangeEnum(DummyEnum,enum_in);
   this->inputs2->DeleteInput(DummyEnum);
 
}/*}}}*/
void FemModel::GroundinglineMassFluxx(IssmDouble* pM, bool scaled){/*{{{*/
 
   /*First we need to depth average the velocities*/
   int domaintype;
   this->parameters->FindParam(&domaintype,DomainTypeEnum);
   this->parameters->SetParam(VxEnum,InputToDepthaverageInEnum);
   this->parameters->SetParam(VxAverageEnum,InputToDepthaverageOutEnum);
   depthaverage_core(this);
   if(domaintype!=Domain2DverticalEnum){
      this->parameters->SetParam(VyEnum,InputToDepthaverageInEnum);
      this->parameters->SetParam(VyAverageEnum,InputToDepthaverageOutEnum);
      depthaverage_core(this);
   }
 
   IssmDouble local_mass_flux = 0;
   IssmDouble total_mass_flux;
 
   for(int i=0;i<this->elements->Size();i++){
      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
      local_mass_flux+=element->GroundinglineMassFlux(scaled);
   }
   ISSM_MPI_Reduce(&local_mass_flux,&total_mass_flux,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
   ISSM_MPI_Bcast(&total_mass_flux,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
 
   /*Assign output pointers: */
   *pM=total_mass_flux;
 
}/*}}}*/
void FemModel::IceMassx(IssmDouble* pM, bool scaled){/*{{{*/
 
   IssmDouble local_ice_mass = 0;
   IssmDouble total_ice_mass;
 
   for(int i=0;i<this->elements->Size();i++){
      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
      local_ice_mass+=element->IceMass(scaled);
   }
   ISSM_MPI_Reduce(&local_ice_mass,&total_ice_mass,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
   ISSM_MPI_Bcast(&total_ice_mass,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
 
   /*Assign output pointers: */
   *pM=total_ice_mass;
 
}/*}}}*/
void FemModel::IceVolumeAboveFloatationx(IssmDouble* pV, bool scaled){/*{{{*/
 
   IssmDouble local_ice_volume_af = 0;
   IssmDouble total_ice_volume_af;
 
   for(int i=0;i<this->elements->Size();i++){
      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
      local_ice_volume_af+=element->IceVolumeAboveFloatation(scaled);
   }
   ISSM_MPI_Reduce(&local_ice_volume_af,&total_ice_volume_af,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
   ISSM_MPI_Bcast(&total_ice_volume_af,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
 
   /*Assign output pointers: */
   *pV=total_ice_volume_af;
 
}/*}}}*/
void FemModel::IceVolumex(IssmDouble* pV, bool scaled){/*{{{*/
 
   IssmDouble local_ice_volume = 0;
   IssmDouble total_ice_volume;
 
   for(int i=0;i<this->elements->Size();i++){
      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
      local_ice_volume+=element->IceVolume(scaled);
   }
   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::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=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
      Input2*  input = element->GetInput2(VxEnum);
      element_maxabsvx=input->GetInputMaxAbs();
      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=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
      Input2*  input = element->GetInput2(VyEnum);
      element_maxabsvy=input->GetInputMaxAbs();
      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=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
      Input2*  input = element->GetInput2(VzEnum);
      element_maxabsvz=input->GetInputMaxAbs();
      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::MaxDivergencex(IssmDouble* pdiv){/*{{{*/
 
   IssmDouble local_divergence;
   IssmDouble node_max_divergence;
   IssmDouble max_divergence = -INFINITY;
 
   for(int i=0;i<this->elements->Size();i++){
      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
      local_divergence=element->Divergence();
      if(fabs(local_divergence)>max_divergence) max_divergence=fabs(local_divergence);
   }
   ISSM_MPI_Reduce(&max_divergence,&node_max_divergence,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
   ISSM_MPI_Bcast(&node_max_divergence,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
   max_divergence=node_max_divergence;
 
   /*Assign output pointers: */
   *pdiv=max_divergence;
 
}/*}}}*/
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=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
      Input2* vel_input = element->GetInput2(VelEnum); _assert_(vel_input);
      element_maxvel = vel_input->GetInputMax();
      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=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
      Input2* vx_input = element->GetInput2(VxEnum); _assert_(vx_input);
      element_maxvx = vx_input->GetInputMax();
      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=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
      Input2* vy_input = element->GetInput2(VyEnum); _assert_(vy_input);
      element_maxvy = vy_input->GetInputMax();
      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=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
      Input2* vz_input = element->GetInput2(VzEnum); _assert_(vz_input);
      element_maxvz = vz_input->GetInputMax();
      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=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
      Input2*  input = element->GetInput2(VelEnum);
      element_minvel =input->GetInputMin();
      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=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
      Input2*  input = element->GetInput2(VxEnum);
      element_minvx =input->GetInputMin();
      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=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
      Input2*  input = element->GetInput2(VyEnum);
      element_minvy =input->GetInputMin();
      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=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
      Input2*  input = element->GetInput2(VzEnum);
      element_minvz =input->GetInputMin();
      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::OmegaAbsGradientx( IssmDouble* pJ){/*{{{*/
 
   /*output: */
   IssmDouble J=0.;
   IssmDouble J_sum;
 
   IssmDouble  omega,weight;
   IssmDouble  Jdet;
   IssmDouble* xyz_list = NULL;
   IssmDouble  dp[3];
 
   /*Compute Misfit: */
   for(int i=0;i<elements->Size();i++){
      Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
 
      /*If on water, return 0: */
      if(!element->IsIceInElement()) continue;
 
      /* Get node coordinates*/
      element->GetVerticesCoordinates(&xyz_list);
 
      /*Retrieve all inputs we will be needing: */
      DatasetInput2* weights_input = element->GetDatasetInput2(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
      Input2* omega_input   = element->GetInput2(BalancethicknessOmegaEnum);              _assert_(omega_input);
 
      /* Start  looping on the number of gaussian points: */
      Gauss* gauss=element->NewGauss(2);
      for(int ig=gauss->begin();ig<gauss->end();ig++){
 
         gauss->GaussPoint(ig);
 
         /* Get Jacobian determinant: */
         element->JacobianDeterminant(&Jdet,xyz_list,gauss);
 
         /*Get all parameters at gaussian point*/
         weights_input->GetInputValue(&weight,gauss,OmegaAbsGradientEnum);
         omega_input->GetInputDerivativeValue(&dp[0],xyz_list,gauss);
 
         /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */
         //J+=weight*1/2*(dp[0]*dp[0]+dp[1]*dp[1])*Jdet*gauss->weight;
         J+=weight*1/2*pow(dp[0]*dp[0]+dp[1]*dp[1],2)*Jdet*gauss->weight;
      }
 
      /*clean up and Return: */
      xDelete<IssmDouble>(xyz_list);
      delete gauss;
   }
 
   /*Sum all J from all cpus of the cluster:*/
   ISSM_MPI_Reduce (&J,&J_sum,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
   ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
   J=J_sum;
 
   /*Assign output pointers: */
   *pJ=J;
}
/*}}}*/
void FemModel::EtaDiffx( IssmDouble* pJ){/*{{{*/
 
   /*output: */
   IssmDouble J=0.;
   IssmDouble J_sum;
 
   IssmDouble  omega,weight;
   IssmDouble  Jdet;
   IssmDouble* xyz_list = NULL;
   IssmDouble  p,p0;
 
   /*Compute Misfit: */
   for(int i=0;i<elements->Size();i++){
      Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
 
      /*If on water, return 0: */
      if(!element->IsIceInElement()) continue;
 
      /* Get node coordinates*/
      element->GetVerticesCoordinates(&xyz_list);
 
      /*Retrieve all inputs we will be needing: */
      DatasetInput2* weights_input =element->GetDatasetInput2(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
      Input2* omega_input   =element->GetInput2(BalancethicknessOmegaEnum);              _assert_(omega_input);
      Input2* omega0_input  =element->GetInput2(BalancethicknessOmega0Enum);             _assert_(omega0_input);
 
      /* Start  looping on the number of gaussian points: */
      Gauss* gauss=element->NewGauss(2);
      for(int ig=gauss->begin();ig<gauss->end();ig++){
 
         gauss->GaussPoint(ig);
 
         /* Get Jacobian determinant: */
         element->JacobianDeterminant(&Jdet,xyz_list,gauss);
 
         /*Get all parameters at gaussian point*/
         weights_input->GetInputValue(&weight,gauss,EtaDiffEnum);
         omega_input->GetInputValue(&p,gauss);
         omega0_input->GetInputValue(&p0,gauss);
 
         /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */
         //J+=weight*1/2*(dp[0]*dp[0]+dp[1]*dp[1])*Jdet*gauss->weight;
         J+=weight*1/2*pow(p - p0,2)*Jdet*gauss->weight;
      }
 
      /*clean up and Return: */
      xDelete<IssmDouble>(xyz_list);
      delete gauss;
   }
 
   /*Sum all J from all cpus of the cluster:*/
   ISSM_MPI_Reduce (&J,&J_sum,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
   ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
   J=J_sum;
 
   /*Assign output pointers: */
   *pJ=J;
}
/*}}}*/
void FemModel::OutputControlsx(Results **presults){/*{{{*/
 
   /*parameters: */
   int         num_controls,step;
   IssmDouble  time;
   int        *control_type = NULL;
 
   /*recover results*/
   Results* results = *presults;
   if(!results) results = new Results();
 
   /*Get list of Controls*/
   this->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
   this->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
   this->parameters->FindParam(&step,StepEnum);
   this->parameters->FindParam(&time,TimeEnum);
 
   for(int i=0;i<num_controls;i++){
 
      int control_enum = control_type[i];
      int gradient_enum;
 
      switch(i){
         case 0: gradient_enum = Gradient1Enum; break;
         case 1: gradient_enum = Gradient2Enum; break;
         case 2: gradient_enum = Gradient3Enum; break;
         default: _error_("more than 3 controls not implemented yet");
      }
 
      /*Allocate vector*/
      Vector<IssmPDouble> *vector_control  = new Vector<IssmPDouble>(this->vertices->NumberOfVertices());
      Vector<IssmPDouble> *vector_gradient = new Vector<IssmPDouble>(this->vertices->NumberOfVertices());
 
      /*Fill in vector*/
      for(int j=0;j<elements->Size();j++){
         Element* element=(Element*)elements->GetObjectByOffset(j);
         element->ControlToVectors(vector_control,vector_gradient,control_enum);
      }
      vector_control->Assemble();
      vector_gradient->Assemble();
 
      results->AddResult(new GenericExternalResult<Vector<IssmPDouble>*>(results->Size()+1,control_enum,vector_control ,step,time));
      results->AddResult(new GenericExternalResult<Vector<IssmPDouble>*>(results->Size()+1,gradient_enum,vector_gradient,step,time));
   }
 
   /*Clean up and return*/
   xDelete<int>(control_type);
}
/*}}}*/
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_AD_
      parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum);
      parameters->FindParam(&dependent_objects,AutodiffDependentObjectsEnum);
      if(num_dependents){
         dependents=xNew<IssmPDouble>(num_dependents);
 
         #if defined(_HAVE_CODIPACK_)
         auto& tape_codi = IssmDouble::getGlobalTape();
         #endif
 
         /*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) {
               #if defined(_HAVE_CODIPACK_)
                  tape_codi.registerOutput(output_value);
                  dependents[i] = output_value.getValue();
                  codi_global.output_indices.push_back(output_value.getGradientData());
               #else
                  output_value>>=dependents[i];
               #endif
            }
         }
      }
      delete dependent_objects;
      if(num_dependents)xDelete<IssmPDouble>(dependents);
      #else
      _error_("Should not be requesting dependents when an AD library is not available!");
      #endif
   }
}
/*}}}*/
void FemModel::RequestedOutputsx(Results **presults,char** requested_outputs, int numoutputs, bool save_results){/*{{{*/
 
   /*Intermediaries*/
   bool        isvec,results_on_nodes;
   int         step,output_enum,numonnodes;
   IssmDouble  time;
   IssmDouble  double_result;
   const char *output_string = NULL;
   char**      resultsonnodes = NULL;
 
   /*recover results*/
   Results* results = *presults;
   if(!results) results = new Results();
 
   /*Get time and step*/
   parameters->FindParam(&step,StepEnum);
   parameters->FindParam(&time,TimeEnum);
   parameters->FindParam(&numonnodes,SettingsNumResultsOnNodesEnum);
   if(numonnodes) parameters->FindParam(&resultsonnodes,&numonnodes,SettingsResultsOnNodesEnum);
 
 
   /*Go through all requested output*/
   for(int i=0;i<numoutputs;i++){
      output_string = requested_outputs[i];
      output_enum   = StringToEnumx(output_string,false);
      isvec         = false;
 
 
      /*If string is not an enum, it is defined in output definitions*/
      if(output_enum<0){
         double_result = OutputDefinitionsResponsex(this,output_string);
         if(save_results){
            results->AddResult(new GenericExternalResult<IssmPDouble>(results->Size()+1,output_string,reCast<IssmPDouble>(double_result),step,time));
            continue;
         }
      }
      else{
         /*last chance for the output definition, if the enum is one of Outputdefinition[1-10]Enum:*/
         if(output_enum>=Outputdefinition1Enum && output_enum <=Outputdefinition100Enum){
            double_result = OutputDefinitionsResponsex(this,output_enum);
            if(save_results){
               results->AddResult(new GenericExternalResult<IssmPDouble>(results->Size()+1,output_string,reCast<IssmPDouble>(double_result),step,time));
               continue;
            }
         }
         else{
            switch(output_enum){
 
               /*Scalar output*/
               case DivergenceEnum:                     this->Divergencex(&double_result);                     break;
               case MaxDivergenceEnum:                  this->MaxDivergencex(&double_result);                  break;
               case IceMassEnum:                        this->IceMassx(&double_result,false);                  break;
               case IcefrontMassFluxEnum:               this->IcefrontMassFluxx(&double_result,false);         break;
               case IcefrontMassFluxLevelsetEnum:       this->IcefrontMassFluxLevelsetx(&double_result,false);         break;
               case IceMassScaledEnum:                  this->IceMassx(&double_result,true);                   break;
               case IceVolumeEnum:                      this->IceVolumex(&double_result,false);                break;
               case IceVolumeScaledEnum:                this->IceVolumex(&double_result,true);                 break;
               case IceVolumeAboveFloatationEnum:       this->IceVolumeAboveFloatationx(&double_result,false); break;
               case IceVolumeAboveFloatationScaledEnum: this->IceVolumeAboveFloatationx(&double_result,true);  break;
               case GroundedAreaEnum:                   this->GroundedAreax(&double_result,false);             break;
               case GroundedAreaScaledEnum:             this->GroundedAreax(&double_result,true);              break;
               case GroundinglineMassFluxEnum:          this->GroundinglineMassFluxx(&double_result,false);    break;
               case FloatingAreaEnum:                   this->FloatingAreax(&double_result,false);             break;
               case FloatingAreaScaledEnum:             this->FloatingAreax(&double_result,true);              break;
               case MinVelEnum:                         this->MinVelx(&double_result);                         break;
               case MaxVelEnum:                         this->MaxVelx(&double_result);                         break;
               case MinVxEnum:                          this->MinVxx(&double_result);                          break;
               case MaxVxEnum:                          this->MaxVxx(&double_result);                          break;
               case MaxAbsVxEnum:                       this->MaxAbsVxx(&double_result);                       break;
               case MinVyEnum:                          this->MinVyx(&double_result);                          break;
               case MaxVyEnum:                          this->MaxVyx(&double_result);                          break;
               case MaxAbsVyEnum:                       this->MaxAbsVyx(&double_result);                       break;
               case MinVzEnum:                          this->MinVzx(&double_result);                          break;
               case MaxVzEnum:                          this->MaxVzx(&double_result);                          break;
               case MaxAbsVzEnum:                       this->MaxAbsVzx(&double_result);                       break;
               case MassFluxEnum:                       this->MassFluxx(&double_result);                       break;
               case TotalCalvingFluxLevelsetEnum:       this->TotalCalvingFluxLevelsetx(&double_result,false);          break;
               case TotalCalvingMeltingFluxLevelsetEnum:this->TotalCalvingMeltingFluxLevelsetx(&double_result,false);          break;
               case TotalFloatingBmbEnum:               this->TotalFloatingBmbx(&double_result,false);         break;
               case TotalFloatingBmbScaledEnum:         this->TotalFloatingBmbx(&double_result,true);          break;
               case TotalGroundedBmbEnum:               this->TotalGroundedBmbx(&double_result,false);         break;
               case TotalGroundedBmbScaledEnum:         this->TotalGroundedBmbx(&double_result,true);          break;
               case TotalSmbEnum:                       this->TotalSmbx(&double_result,false);                 break;
               case TotalSmbScaledEnum:                 this->TotalSmbx(&double_result,true);                  break;
 
               /*Scalar control output*/
               case SurfaceAbsVelMisfitEnum:       SurfaceAbsVelMisfitx(&double_result,elements,nodes,vertices,loads,materials,parameters);        break;
               case SurfaceRelVelMisfitEnum:       SurfaceRelVelMisfitx(&double_result,elements,nodes,vertices,loads,materials,parameters);        break;
               case SurfaceLogVelMisfitEnum:       SurfaceLogVelMisfitx(&double_result,elements,nodes,vertices,loads,materials,parameters);        break;
               case SurfaceLogVxVyMisfitEnum:      SurfaceLogVxVyMisfitx(&double_result,elements,nodes,vertices,loads,materials,parameters);       break;
               case SurfaceAverageVelMisfitEnum:   SurfaceAverageVelMisfitx(&double_result,this);                                                  break;
               case ThicknessAbsMisfitEnum:        ThicknessAbsMisfitx(&double_result,elements,nodes,vertices,loads,materials,parameters);         break;
               case ThicknessAbsGradientEnum:      this->ThicknessAbsGradientx(&double_result);                                                    break;
               case ThicknessAlongGradientEnum:    ThicknessAlongGradientx(&double_result,elements,nodes,vertices,loads,materials,parameters);     break;
               case ThicknessAcrossGradientEnum:   ThicknessAcrossGradientx(&double_result,elements,nodes,vertices,loads,materials,parameters);    break;
               case ThicknessPositiveEnum:         this->ThicknessPositivex(&double_result);                                                       break;
               case RheologyBbarAbsGradientEnum:   RheologyBbarAbsGradientx(&double_result,elements,nodes,vertices,loads,materials,parameters);    break;
               case RheologyBAbsGradientEnum:      RheologyBAbsGradientx(&double_result,elements,nodes,vertices,loads,materials,parameters);       break;
               case RheologyBInitialguessMisfitEnum:  RheologyBInitialguessMisfitx(&double_result,elements,nodes,vertices,loads,materials,parameters);       break;
               case DragCoefficientAbsGradientEnum:DragCoefficientAbsGradientx(&double_result,elements,nodes,vertices,loads,materials,parameters); break;
               case BalancethicknessMisfitEnum:    BalancethicknessMisfitx(&double_result);                                                        break;
               case SurfaceAbsMisfitEnum:          SurfaceAbsMisfitx(&double_result); break;
               case OmegaAbsGradientEnum:          OmegaAbsGradientx(&double_result); break;
               case EtaDiffEnum:                   EtaDiffx(&double_result); break;
 
               /*Vector special case (maybe should go to specific analysis?)*/
               case ChannelAreaEnum:{
 
                     /*Get Number of Channels*/
                     int numchannels_local=0,numchannels;
                     for(int j=0;j<this->loads->Size();j++){
                        if(this->loads->GetEnum(j)==ChannelEnum) numchannels_local++;
                     }
                     ISSM_MPI_Reduce(&numchannels_local,&numchannels,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,IssmComm::GetComm() );
                     ISSM_MPI_Bcast(&numchannels,1,ISSM_MPI_INT,0,IssmComm::GetComm());
 
                     IssmPDouble* values    = xNewZeroInit<IssmPDouble>(numchannels);
                     IssmPDouble* allvalues = xNew<IssmPDouble>(numchannels);
 
                     /*Fill-in vector*/
                     for(int j=0;j<this->loads->Size();j++){
                        if(this->loads->GetEnum(j)==ChannelEnum){
                           Channel* channel=(Channel*)this->loads->GetObjectByOffset(j);
                           channel->WriteChannelCrossSection(values);
                        }
                     }
 
                     /*Gather from all cpus*/
                     ISSM_MPI_Allreduce((void*)values,(void*)allvalues,numchannels,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
                     xDelete<IssmPDouble>(values);
 
                     if(save_results)results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,output_enum,allvalues,numchannels,1,step,time));
                     xDelete<IssmPDouble>(allvalues);
 
                     isvec = true;
               }
               break;
 
 
               /*Default is always Vector */
               default:
 
               /*Some preliminary calculation may be required (use similar syntax for other inputs)*/
                  if(output_enum==NewDamageEnum){
                     InputDuplicatex(this,DamageDEnum,DamageDOldEnum);
                     InputDuplicatex(this,DamageDbarEnum,DamageDbarOldEnum);
                     this->ElementOperationx(&Element::ComputeNewDamage);
                  }
 
                  /*Vector layout*/
                  if(!IsInputEnum(output_enum)) _error_("Cannot output \""<<EnumToStringx(output_enum)<<"\" because it is not an input");
                  int interpolation,nodesperelement,size,nlines,ncols,array_size;
                  int rank_interpolation=-1,rank_nodesperelement=-1,rank_arraysize=-1,max_rank_arraysize=0;
                  bool isarray=false;
 
                  /*Get interpolation (and compute input if necessary)*/
                  for(int j=0;j<elements->Size();j++){
                     Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(j));
                     element->ResultInterpolation(&rank_interpolation,&rank_nodesperelement,&rank_arraysize,output_enum);
                     if(rank_arraysize>max_rank_arraysize)max_rank_arraysize=rank_arraysize;
                  }
                  rank_arraysize=max_rank_arraysize;
 
                  /*Broadcast for cpus that do not have any elements*/
                  ISSM_MPI_Reduce(&rank_interpolation,&interpolation,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm());
                  ISSM_MPI_Reduce(&rank_nodesperelement,&nodesperelement,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm());
                  ISSM_MPI_Reduce(&rank_arraysize,&array_size,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm());
                  ISSM_MPI_Bcast(&interpolation,1,ISSM_MPI_INT,0,IssmComm::GetComm());
                  ISSM_MPI_Bcast(&nodesperelement,1,ISSM_MPI_INT,0,IssmComm::GetComm());
                  ISSM_MPI_Bcast(&array_size,1,ISSM_MPI_INT,0,IssmComm::GetComm());
 
                  results_on_nodes=false;
                  /*Loop to see if this output was requested on nodes*/
                  for(int j=0;j<numonnodes & results_on_nodes==false;j++){
                     if(strcmp(resultsonnodes[j],output_string) == 0 || strcmp(resultsonnodes[j],"all") == 0) results_on_nodes=true;
                  }
 
                  if(results_on_nodes){
 
                     /*Allocate matrices*/
                     int         nbe       = this->elements->NumberOfElements();
                     IssmDouble* values    = xNewZeroInit<IssmDouble>(nbe*nodesperelement);
                     IssmDouble* allvalues = xNew<IssmDouble>(nbe*nodesperelement);
 
                     /*Fill-in matrix*/
                     for(int j=0;j<elements->Size();j++){
                        Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(j));
                        element->ResultToPatch(values,nodesperelement,output_enum);
                     }
 
                     /*Gather from all cpus*/
                     ISSM_MPI_Allreduce((void*)values,(void*)allvalues,nbe*nodesperelement,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
                     xDelete<IssmDouble>(values);
 
                     if(save_results)results->AddResult(new GenericExternalResult<IssmDouble*>(results->Size()+1,output_enum,allvalues,nbe,nodesperelement,step,time));
                     xDelete<IssmDouble>(allvalues);
 
                  }
                  else{
 
                     /*Allocate vector depending on interpolation*/
                     switch(interpolation){
                        case P0Enum: isarray = false; size = this->elements->NumberOfElements(); break;
                        case P1Enum: isarray = false; size = this->vertices->NumberOfVertices(); break;
                        case P0ArrayEnum: isarray = true; nlines = this->elements->NumberOfElements(); ncols= array_size; break;
                        default:     _error_("Interpolation "<<EnumToStringx(interpolation)<<" not supported yet");
 
                     }
                     if(!isarray){
                        Vector<IssmDouble> *vector_result = new Vector<IssmDouble>(size);
 
                        /*Fill in vector*/
                        for(int j=0;j<elements->Size();j++){
                           Element* element=(Element*)elements->GetObjectByOffset(j);
                           element->ResultToVector(vector_result,output_enum);
                        }
                        vector_result->Assemble();
 
                        if(save_results){
                           results->AddResult(new GenericExternalResult<Vector<IssmDouble>*>(results->Size()+1,output_enum,vector_result,step,time));
                           /*We do not do a copy for Vectors, so don't delete*/
                        }
                        else{
                           delete vector_result;
                        }
                     }
                     else{
                        IssmDouble* values    = xNewZeroInit<IssmDouble>(nlines*ncols);
                        IssmDouble* allvalues = xNew<IssmDouble>(nlines*ncols);
 
                        /*Fill-in matrix*/
                        for(int j=0;j<elements->Size();j++){
                           Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(j));
                           element->ResultToMatrix(values,ncols,output_enum);
                        }
                        /*Gather from all cpus*/
                        ISSM_MPI_Allreduce((void*)values,(void*)allvalues,ncols*nlines,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
                        xDelete<IssmDouble>(values);
 
                        if(save_results)results->AddResult(new GenericExternalResult<IssmDouble*>(results->Size()+1,output_enum,allvalues,nlines,ncols,step,time));
                        xDelete<IssmDouble>(allvalues);
                     }
                  }
                  isvec = true;
                  break;
            }
         }
 
      }
 
      /*Add result to Results*/
      if(!isvec && save_results){
         results->AddResult(new GenericExternalResult<IssmPDouble>(results->Size()+1,output_string,reCast<IssmPDouble>(double_result),step,time));
      }
   }
 
   /*Clean up*/
   for(int i=0;i<numonnodes;i++) xDelete<char>(resultsonnodes[i]);
   xDelete<char*>(resultsonnodes);
 
   /*Assign pointer and clean up*/
   *presults = results;
}
/*}}}*/
void FemModel::RequestedOutputsx(Results **presults,int* requested_outputs, int numoutputs,bool save_results){/*{{{*/
 
   /*Convert list of enums to list of string*/
   char** enumlist = xNew<char*>(numoutputs);
   for(int i=0;i<numoutputs;i++) EnumToStringx(&enumlist[i],requested_outputs[i]);
 
   /*Call main module*/
   this->RequestedOutputsx(presults,enumlist,numoutputs,save_results);
 
   /*clean up and return*/
   for(int i=0;i<numoutputs;i++) xDelete<char>(enumlist[i]);
   xDelete<char*>(enumlist);
}
/*}}}*/
void FemModel::ResetLevelset(void){/*{{{*/
 
   this->DistanceToFieldValue(MaskIceLevelsetEnum,0.,MaskIceLevelsetEnum);
 
}
/*}}}*/
void FemModel::Responsex(IssmDouble* responses,const char* response_descriptor){/*{{{*/
 
   int response_descriptor_enum=StringToEnumx(response_descriptor);
   this->Responsex(responses, response_descriptor_enum);
 
}
/*}}}*/
void FemModel::Responsex(IssmDouble* responses,int response_descriptor_enum){/*{{{*/
 
   switch (response_descriptor_enum){
 
      case DivergenceEnum:                     this->Divergencex(responses); break;
      case MaxDivergenceEnum:                  this->MaxDivergencex(responses); break;
      case IceMassEnum:                        this->IceMassx(responses, false); break;
      case IceMassScaledEnum:                  this->IceMassx(responses, true); break;
      case IceVolumeEnum:                      this->IceVolumex(responses, false); break;
      case IceVolumeScaledEnum:                this->IceVolumex(responses, true); break;
      case IceVolumeAboveFloatationEnum:       this->IceVolumeAboveFloatationx(responses, false); break;
      case IceVolumeAboveFloatationScaledEnum: this->IceVolumeAboveFloatationx(responses, true); break;
      case IcefrontMassFluxEnum:               this->IcefrontMassFluxx(responses, false); break;
      case IcefrontMassFluxLevelsetEnum:       this->IcefrontMassFluxLevelsetx(responses, false); break;
      case GroundedAreaEnum:                   this->GroundedAreax(responses, false); break;
      case GroundedAreaScaledEnum:             this->GroundedAreax(responses, true); break;
      case FloatingAreaEnum:                   this->FloatingAreax(responses, false); break;
      case FloatingAreaScaledEnum:             this->FloatingAreax(responses, true); break;
      case MinVelEnum:                         this->MinVelx(responses); break;
      case MaxVelEnum:                         this->MaxVelx(responses); break;
      case MinVxEnum:                          this->MinVxx(responses); break;
      case MaxVxEnum:                          this->MaxVxx(responses); break;
      case MaxAbsVxEnum:                       this->MaxAbsVxx(responses); break;
      case MinVyEnum:                          this->MinVyx(responses); break;
      case MaxVyEnum:                          this->MaxVyx(responses); break;
      case MaxAbsVyEnum:                       this->MaxAbsVyx(responses); break;
      case MinVzEnum:                          this->MinVzx(responses); break;
      case MaxVzEnum:                          this->MaxVzx(responses); break;
      case MaxAbsVzEnum:                       this->MaxAbsVzx(responses); break;
      case MassFluxEnum:                       this->MassFluxx(responses); break;
      case SurfaceAbsVelMisfitEnum:            SurfaceAbsVelMisfitx(responses, elements,nodes, vertices, loads, materials,parameters); break;
      case SurfaceRelVelMisfitEnum:            SurfaceRelVelMisfitx(responses, elements,nodes, vertices, loads, materials,parameters); break;
      case SurfaceLogVelMisfitEnum:            SurfaceLogVelMisfitx(responses, elements,nodes, vertices, loads, materials,parameters); break;
      case SurfaceLogVxVyMisfitEnum:           SurfaceLogVxVyMisfitx(responses, elements,nodes, vertices, loads, materials,parameters); break;
      case SurfaceAverageVelMisfitEnum:        SurfaceAverageVelMisfitx(responses,this); break;
      case ThicknessAbsMisfitEnum:             ThicknessAbsMisfitx(responses, elements,nodes, vertices, loads, materials, parameters); break;
      case ThicknessAbsGradientEnum:           this->ThicknessAbsGradientx(responses); break;
      case ThicknessAlongGradientEnum:         ThicknessAlongGradientx(responses, elements,nodes, vertices, loads, materials, parameters); break;
      case ThicknessAcrossGradientEnum:        ThicknessAcrossGradientx(responses, elements,nodes, vertices, loads, materials, parameters); break;
      case RheologyBbarAbsGradientEnum:        RheologyBbarAbsGradientx(responses, elements,nodes, vertices, loads, materials, parameters); break;
      case DragCoefficientAbsGradientEnum:     DragCoefficientAbsGradientx(responses, elements,nodes, vertices, loads, materials, parameters); break;
      case BalancethicknessMisfitEnum:         BalancethicknessMisfitx(responses); break;
      case TotalCalvingFluxLevelsetEnum:       this->TotalCalvingFluxLevelsetx(responses, false); break;
      case TotalCalvingMeltingFluxLevelsetEnum:this->TotalCalvingMeltingFluxLevelsetx(responses, false); break;
      case TotalFloatingBmbEnum:               this->TotalFloatingBmbx(responses, false); break;
      case TotalFloatingBmbScaledEnum:         this->TotalFloatingBmbx(responses, true); break;
      case TotalGroundedBmbEnum:               this->TotalGroundedBmbx(responses, false); break;
      case TotalGroundedBmbScaledEnum:         this->TotalGroundedBmbx(responses, true); break;
      case TotalSmbEnum:                       this->TotalSmbx(responses, false); break;
      case TotalSmbScaledEnum:                 this->TotalSmbx(responses, true); 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:
         if(response_descriptor_enum>=Outputdefinition1Enum && response_descriptor_enum <=Outputdefinition100Enum){
            IssmDouble double_result = OutputDefinitionsResponsex(this,response_descriptor_enum);
            *responses=double_result;
         }
         else _error_("response descriptor \"" << EnumToStringx(response_descriptor_enum) << "\" not supported yet!");
         break;
   }
 
}
/*}}}*/
void FemModel::RignotMeltParameterizationx(){/*{{{*/
 
   for(int i=0;i<elements->Size();i++){
      Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
      element->RignotMeltParameterization();
   }
}
/*}}}*/
void FemModel::StrainRateparallelx(){/*{{{*/
 
   for(int i=0;i<elements->Size();i++){
      Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
      element->StrainRateparallel();
   }
}
/*}}}*/
void FemModel::StrainRateperpendicularx(){/*{{{*/
 
   for(int i=0;i<elements->Size();i++){
      Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
      element->StrainRateperpendicular();
   }
}
/*}}}*/
void FemModel::StrainRateeffectivex(){/*{{{*/
 
   for(int i=0;i<elements->Size();i++){
      Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
      element->ComputeStrainRate();
   }
}
/*}}}*/
void FemModel::StressIntensityFactorx(){/*{{{*/
 
   /*Update input for basal element only*/
   for(int i=0;i<elements->Size();i++){
      Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
      element->StressIntensityFactor();
   }
}
   /*}}}*/
void FemModel::SurfaceAbsMisfitx(IssmDouble* presponse){/*{{{*/
 
   /*output: */
   IssmDouble J=0.;
   IssmDouble J_sum;
 
   IssmDouble  surface,surfaceobs,weight;
   IssmDouble  Jdet;
   IssmDouble* xyz_list = NULL;
 
   /*Compute Misfit: */
   for(int i=0;i<elements->Size();i++){
      Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
 
      /*If on water, return 0: */
      if(!element->IsIceInElement()) continue;
 
      /* Get node coordinates*/
      element->GetVerticesCoordinates(&xyz_list);
 
       /*Retrieve all inputs we will be needing: */
       DatasetInput2* weights_input   =element->GetDatasetInput2(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
       Input2* surface_input   =element->GetInput2(SurfaceEnum);                            _assert_(surface_input);
       Input2* surfaceobs_input=element->GetInput2(InversionSurfaceObsEnum);                _assert_(surfaceobs_input);
 
       /* Start  looping on the number of gaussian points: */
       Gauss* gauss=element->NewGauss(2);
       for(int ig=gauss->begin();ig<gauss->end();ig++){
 
          gauss->GaussPoint(ig);
 
          /* Get Jacobian determinant: */
          element->JacobianDeterminant(&Jdet,xyz_list,gauss);
 
          /*Get all parameters at gaussian point*/
          weights_input->GetInputValue(&weight,gauss,SurfaceAbsMisfitEnum);
          surface_input->GetInputValue(&surface,gauss);
          surfaceobs_input->GetInputValue(&surfaceobs,gauss);
 
          /*Compute SurfaceAbsMisfitEnum*/
          J+=0.5*(surface-surfaceobs)*(surface-surfaceobs)*weight*Jdet*gauss->weight;
       }
       delete gauss;
       xDelete<IssmDouble>(xyz_list);
   }
 
   /*Sum all J from all cpus of the cluster:*/
   ISSM_MPI_Reduce (&J,&J_sum,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
   ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
   J=J_sum;
 
   /*Assign output pointers: */
   *presponse=J;
 
}/*}}}*/
void FemModel::ThicknessAbsGradientx( IssmDouble* pJ){/*{{{*/
 
   /*output: */
   IssmDouble J=0.;
   IssmDouble J_sum;
 
   IssmDouble  thickness,weight;
   IssmDouble  Jdet;
   IssmDouble* xyz_list = NULL;
   IssmDouble  dp[3];
 
   /*Compute Misfit: */
   for(int i=0;i<elements->Size();i++){
      Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
 
      /*If on water, return 0: */
      if(!element->IsIceInElement()) continue;
 
      /* Get node coordinates*/
      element->GetVerticesCoordinates(&xyz_list);
 
      /*Retrieve all inputs we will be needing: */
      DatasetInput2* weights_input   =element->GetDatasetInput2(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
      Input2* thickness_input =element->GetInput2(ThicknessEnum);                          _assert_(thickness_input);
 
      /* Start  looping on the number of gaussian points: */
      Gauss* gauss=element->NewGauss(2);
      for(int ig=gauss->begin();ig<gauss->end();ig++){
 
         gauss->GaussPoint(ig);
 
         /* Get Jacobian determinant: */
         element->JacobianDeterminant(&Jdet,xyz_list,gauss);
 
         /*Get all parameters at gaussian point*/
         weights_input->GetInputValue(&weight,gauss,ThicknessAbsGradientEnum);
         thickness_input->GetInputDerivativeValue(&dp[0],xyz_list,gauss);
 
         /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */
         J+=weight*1/2*(dp[0]*dp[0]+dp[1]*dp[1])*Jdet*gauss->weight;
      }
 
      /*clean up and Return: */
      xDelete<IssmDouble>(xyz_list);
      delete gauss;
   }
 
   /*Sum all J from all cpus of the cluster:*/
   ISSM_MPI_Reduce (&J,&J_sum,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
   ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
   J=J_sum;
 
   /*Assign output pointers: */
   *pJ=J;
}
/*}}}*/
void FemModel::ThicknessAverage(){/*{{{*/
 
   int elementswidth                   = this->GetElementsWidth();//just 2D mesh, tria elements
   int numberofvertices                = this->vertices->NumberOfVertices();//total number of vertices
 
   IssmDouble weight                   = 0.;
   IssmDouble* totalweight             = NULL;
   IssmDouble* Hserial                 = NULL;
   IssmDouble* H                       = xNew<IssmDouble>(elementswidth);
   Vector<IssmDouble>* vecH            = new Vector<IssmDouble>(numberofvertices);
   Vector<IssmDouble>* vectotalweight  = new Vector<IssmDouble>(numberofvertices);
 
   for(int i=0;i<this->elements->Size();i++){
      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
 
      /*check if there is ice in this element*/
      if(!element->IsIceInElement()) continue;
 
      /*get H on the vertices*/
      element->GetInputListOnVertices(H,ThicknessEnum);
 
      /*weight to calculate the smoothed H*/
      weight=1.;//simple average
 
      /*add in the serial vector*/
      vecH->SetValue(element->vertices[0]->Sid(),weight*H[0],ADD_VAL);
      vecH->SetValue(element->vertices[1]->Sid(),weight*H[1],ADD_VAL);
      vecH->SetValue(element->vertices[2]->Sid(),weight*H[2],ADD_VAL);
      /*total weight*/
      vectotalweight->SetValue(element->vertices[0]->Sid(),weight,ADD_VAL);
      vectotalweight->SetValue(element->vertices[1]->Sid(),weight,ADD_VAL);
      vectotalweight->SetValue(element->vertices[2]->Sid(),weight,ADD_VAL);
   }
 
   /*Assemble and serialize*/
   vecH->Assemble();
   vectotalweight->Assemble();
   Hserial=vecH->ToMPISerial();
   totalweight=vectotalweight->ToMPISerial();
 
   /*Divide for the total weight*/
   for(int i=0;i<numberofvertices;i++){
      _assert_(totalweight[i]>0);
      Hserial[i]=Hserial[i]/totalweight[i];
   }
 
   /*Set element inputs*/
   for(int i=0;i<this->elements->Size();i++){
      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
      H[0]=Hserial[element->vertices[0]->Sid()];
      H[1]=Hserial[element->vertices[1]->Sid()];
      H[2]=Hserial[element->vertices[2]->Sid()];
      element->AddInput2(ThicknessEnum,H,P1Enum);
   }
 
   /*Cleanup*/
   delete vecH;
   delete vectotalweight;
   xDelete<IssmDouble>(H);
   xDelete<IssmDouble>(Hserial);
   xDelete<IssmDouble>(totalweight);
}
/*}}}*/
void FemModel::ThicknessPositivex(IssmDouble* pJ){/*{{{*/
 
   /*output: */
   IssmDouble J=0.;
   IssmDouble J_sum;
 
   IssmDouble  thickness,weight;
   IssmDouble  Jdet;
   IssmDouble* xyz_list = NULL;
   IssmDouble  H;
 
   /*Compute Misfit: */
   for(int i=0;i<elements->Size();i++){
      Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
 
      /*If on water, return 0: */
      if(!element->IsIceInElement()) continue;
 
      /* Get node coordinates*/
      element->GetVerticesCoordinates(&xyz_list);
 
      /*Retrieve all inputs we will be needing: */
      DatasetInput2* weights_input   =element->GetDatasetInput2(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
      Input2* thickness_input =element->GetInput2(ThicknessEnum);                          _assert_(thickness_input);
 
      /* Start  looping on the number of gaussian points: */
      Gauss* gauss=element->NewGauss(2);
      for(int ig=gauss->begin();ig<gauss->end();ig++){
 
         gauss->GaussPoint(ig);
 
         /* Get Jacobian determinant: */
         element->JacobianDeterminant(&Jdet,xyz_list,gauss);
 
         /*Get all parameters at gaussian point*/
         weights_input->GetInputValue(&weight,gauss,ThicknessPositiveEnum);
         thickness_input->GetInputValue(&H,gauss);
 
         /*int min(H,0)^2 */
         if(H<=0){
            J+=weight*H*H*Jdet*gauss->weight;
         }
      }
 
      /*clean up and Return: */
      xDelete<IssmDouble>(xyz_list);
      delete gauss;
   }
 
   /*Sum all J from all cpus of the cluster:*/
   ISSM_MPI_Reduce (&J,&J_sum,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
   ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
   J=J_sum;
 
   /*Assign output pointers: */
   *pJ=J;
}
/*}}}*/
void FemModel::TimeAdaptx(IssmDouble* pdt){/*{{{*/
 
   /*output: */
   IssmDouble   dt;
 
   /*intermediary: */
   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=(Element*)elements->GetObjectByOffset(0); min_dt=element->TimeAdapt();
 
   for(int i=1;i<elements->Size();i++){
      element=xDynamicCast<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;
 
   /*Constrain dt */
   IssmDouble dt_low,dt_high;
   parameters->FindParam(&dt_low,TimesteppingTimeStepMinEnum);
   parameters->FindParam(&dt_high,TimesteppingTimeStepMaxEnum);
   if(min_dt<dt_low)  min_dt = dt_low;
   if(min_dt>dt_high) min_dt = dt_high;
 
   /*Assign output pointers:*/
   *pdt=min_dt;
}
/*}}}*/
void FemModel::TotalCalvingFluxLevelsetx(IssmDouble* pM, bool scaled){/*{{{*/
 
   IssmDouble local_calving_flux = 0;
   IssmDouble total_calving_flux;
 
   for(int i=0;i<this->elements->Size();i++){
      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
      local_calving_flux+=element->TotalCalvingFluxLevelset(scaled);
   }
   ISSM_MPI_Reduce(&local_calving_flux,&total_calving_flux,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
   ISSM_MPI_Bcast(&total_calving_flux,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
 
   /*Assign output pointers: */
   *pM=total_calving_flux;
 
}/*}}}*/
void FemModel::TotalCalvingMeltingFluxLevelsetx(IssmDouble* pM, bool scaled){/*{{{*/
 
   IssmDouble local_calving_flux = 0;
   IssmDouble total_calving_flux;
 
   for(int i=0;i<this->elements->Size();i++){
      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
      local_calving_flux+=element->TotalCalvingMeltingFluxLevelset(scaled);
   }
   ISSM_MPI_Reduce(&local_calving_flux,&total_calving_flux,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
   ISSM_MPI_Bcast(&total_calving_flux,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
 
   /*Assign output pointers: */
   *pM=total_calving_flux;
 
}/*}}}*/
void FemModel::TotalFloatingBmbx(IssmDouble* pFbmb, bool scaled){/*{{{*/
 
   IssmDouble local_fbmb = 0;
   IssmDouble total_fbmb;
 
   for(int i=0;i<this->elements->Size();i++){
      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
      local_fbmb+=element->TotalFloatingBmb(scaled);
   }
   ISSM_MPI_Reduce(&local_fbmb,&total_fbmb,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
   ISSM_MPI_Bcast(&total_fbmb,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
 
   /*Assign output pointers: */
   *pFbmb=total_fbmb;
 
}/*}}}*/
void FemModel::TotalGroundedBmbx(IssmDouble* pGbmb, bool scaled){/*{{{*/
 
   IssmDouble local_gbmb = 0;
   IssmDouble total_gbmb;
 
   for(int i=0;i<this->elements->Size();i++){
      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
      local_gbmb+=element->TotalGroundedBmb(scaled);
   }
   ISSM_MPI_Reduce(&local_gbmb,&total_gbmb,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
   ISSM_MPI_Bcast(&total_gbmb,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
 
   /*Assign output pointers: */
   *pGbmb=total_gbmb;
 
}/*}}}*/
void FemModel::TotalSmbx(IssmDouble* pSmb, bool scaled){/*{{{*/
 
   IssmDouble local_smb = 0;
   IssmDouble total_smb;
 
   for(int i=0;i<this->elements->Size();i++){
      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
      local_smb+=element->TotalSmb(scaled);
   }
   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::UpdateConstraintsExtrudeFromBasex(void){ /*{{{*/
 
   for(int i=0;i<elements->Size();i++){
      Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
      element->UpdateConstraintsExtrudeFromBase();
   }
 
}
/*}}}*/
void FemModel::UpdateConstraintsExtrudeFromTopx(void){ /*{{{*/
 
   for(int i=0;i<elements->Size();i++){
      Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
      element->UpdateConstraintsExtrudeFromTop();
   }
 
}
/*}}}*/
void FemModel::UpdateConstraintsx(void){ /*{{{*/
 
   IssmDouble time;
   int        analysis_type,config_type;
 
   /*retrieve parameters: */
   parameters->FindParam(&analysis_type,AnalysisTypeEnum);
   parameters->FindParam(&config_type,ConfigurationTypeEnum);
   parameters->FindParam(&time,TimeEnum);
 
   int index=AnalysisIndex(config_type);
   _assert_(this->analysis_type_list[index]==config_type);
 
   /*start module: */
   if(VerboseModule()) _printf0_("   Updating constraints and active domain of analysis " << EnumToStringx(analysis_type)  << " for time: " << time << "\n");
 
   Analysis* analysis= EnumToAnalysis(analysis_type);
   analysis->UpdateConstraints(this);
   delete analysis;
 
   /*Second, constraints might be time dependent: */
   SpcNodesx(nodes,constraints,parameters);
 
   /*Now, update degrees of freedoms: */
   NodesDofx(nodes,parameters);
 
}
/*}}}*/
int  FemModel::UpdateVertexPositionsx(void){ /*{{{*/
 
   IssmDouble         *surface = NULL;
   IssmDouble         *bed     = NULL;
 
   if(VerboseSolution()) _printf0_("   updating vertices positions\n");
 
   /*get vertex vectors for bed and thickness: */
   GetVectorFromInputsx(&surface  ,this, SurfaceEnum,VertexPIdEnum);
   GetVectorFromInputsx(&bed      ,this, BaseEnum,   VertexPIdEnum);
 
   /*Allocate vector*/
   int numvert       = vertices->NumberOfVertices();
   int numvert_local = vertices->NumberOfVerticesLocal();
   Vector<IssmDouble> *vx=new Vector<IssmDouble>(numvert_local,numvert);
   Vector<IssmDouble> *vy=new Vector<IssmDouble>(numvert_local,numvert);
   Vector<IssmDouble> *vz=new Vector<IssmDouble>(numvert_local,numvert);
 
   /*Update verices new geometry: */
   for(int i=0;i<vertices->Size();i++){
      Vertex* vertex=(Vertex*)vertices->GetObjectByOffset(i);
      vertex->UpdatePosition(vx,vy,vz,parameters,surface,bed);
   }
 
   /*Assemble mesh velocity*/
   vx->Assemble();
   vy->Assemble();
   vz->Assemble();
 
   /*Update element inputs*/
   InputUpdateFromVectorx(this,vx,VxMeshEnum,VertexPIdEnum);
   InputUpdateFromVectorx(this,vy,VyMeshEnum,VertexPIdEnum);
   InputUpdateFromVectorx(this,vz,VzMeshEnum,VertexPIdEnum);
 
   /*Free ressources:*/
   delete vx;
   delete vy;
   delete vz;
   xDelete<IssmDouble>(bed);
   xDelete<IssmDouble>(surface);
   return 1;
}
/*}}}*/
 
/*AMR*/
#ifndef _HAVE_AD_
void FemModel::ReMesh(void){/*{{{*/
 
   /*Intermediaries*/
   IssmDouble *newx        = NULL;
   IssmDouble *newy        = NULL;
   IssmDouble *newz        = NULL;
   int *newelementslist    = NULL;
   int newnumberofvertices = -1;
   int newnumberofelements = -1;
 
   int elementswidth       = this->GetElementsWidth();//just tria elements in this version
   int amrtype,basalforcing_model;
   bool isgroundingline;
 
   /*Branch to specific amr depending on requested method*/
   this->parameters->FindParam(&amrtype,AmrTypeEnum);
   switch(amrtype){
      #if defined(_HAVE_NEOPZ_) && !defined(_HAVE_AD_)
      case AmrNeopzEnum: this->ReMeshNeopz(&newnumberofvertices,&newnumberofelements,&newx,&newy,&newz,&newelementslist); break;
      #endif
 
      #if defined(_HAVE_BAMG_) && !defined(_HAVE_AD_)
      case AmrBamgEnum: this->ReMeshBamg(&newnumberofvertices,&newnumberofelements,&newx,&newy,&newz,&newelementslist); break;
      #endif
 
      default: _error_("not implemented yet");
   }
 
   /*Create iomodel for model processing*/
   IoModel* iomodel = new IoModel();
   this->parameters->FindParam(&iomodel->domaintype,DomainTypeEnum);
   this->parameters->FindParam(&iomodel->domaindim ,DomainDimensionEnum);
   this->parameters->FindParam(&iomodel->meshelementtype,MeshElementtypeEnum);
   iomodel->numberofvertices = newnumberofvertices;
   iomodel->numberofelements = newnumberofelements;
   iomodel->elements         = newelementslist;
   iomodel->AddConstant(new IoConstant(0,"md.rifts.numrifts"));
   iomodel->AddConstant(new IoConstant(false,"md.transient.isoceancoupling"));
   bool temp; int tempint;
   this->parameters->FindParam(&temp,FlowequationIsSIAEnum); iomodel->AddConstant(new IoConstant(temp,"md.flowequation.isSIA"));
   this->parameters->FindParam(&temp,FlowequationIsSSAEnum); iomodel->AddConstant(new IoConstant(temp,"md.flowequation.isSSA"));
   this->parameters->FindParam(&temp,FlowequationIsL1L2Enum); iomodel->AddConstant(new IoConstant(temp,"md.flowequation.isL1L2"));
   this->parameters->FindParam(&temp,FlowequationIsHOEnum); iomodel->AddConstant(new IoConstant(temp,"md.flowequation.isHO"));
   this->parameters->FindParam(&temp,FlowequationIsFSEnum); iomodel->AddConstant(new IoConstant(temp,"md.flowequation.isFS"));
   this->parameters->FindParam(&tempint,MasstransportStabilizationEnum); iomodel->AddConstant(new IoConstant(tempint,"md.masstransport.stabilization"));
   iomodel->AddConstant(new IoConstant(P1Enum,"md.flowequation.fe_SSA"));
 
   /*Partitioning the new mesh. Maybe ElementsAndVerticesPartitioning.cpp could be modified to set this without iomodel.*/
   ::ElementsAndVerticesPartitioning(iomodel);
 
   /*Creating elements*/
   /*Just Tria in this version*/
   Elements* new_elements=new Elements();
   this->CreateElements(newnumberofelements,elementswidth,newelementslist,iomodel->my_elements,new_elements);
 
   /*Create vertices*/
   Vertices* new_vertices=new Vertices();
   CreateNumberNodeToElementConnectivity(iomodel);
   ::CreateVertices(new_elements,new_vertices,iomodel,TransientSolutionEnum,true);
   for(int i=0;i<new_vertices->Size();i++){
      Vertex *vertex=(Vertex*)new_vertices->GetObjectByOffset(i);
      int     sid = vertex->Sid();
      vertex->x=newx[sid];
      vertex->y=newy[sid];
      vertex->z=newz[sid];
   }
 
   /*Creating inputs*/
   Inputs2* new_inputs2=new Inputs2(newnumberofelements,newnumberofvertices);
 
   /*Creating materials*/
   Materials* new_materials=new Materials();
   this->CreateMaterials(newnumberofelements,iomodel->my_elements,new_materials);
 
   /*Creating nodes and constraints*/
   /*Just SSA (2D) and P1 in this version*/
   Constraints **new_constraints_list = xNew<Constraints*>(this->nummodels);
   Nodes       **new_nodes_list       = xNew<Nodes*>(this->nummodels);
 
   this->analysis_counter=-1;
   for(int i=0;i<this->nummodels;i++){//create nodes for each analysis in analysis_type_list
 
      int analysis_enum = this->analysis_type_list[i];
      if(VerboseMProcessor()) _printf0_("   creating datasets for analysis " << EnumToStringx(analysis_enum) << "\n");
 
      if(this->loads_list[i]->Size()!=0) _error_("not supported yet");
      new_constraints_list[i] = new Constraints();
      new_nodes_list[i] = new Nodes();
 
      /*As the domain is 2D, it is not necessary to create nodes for this analysis*/
      if(analysis_enum==StressbalanceVerticalAnalysisEnum) continue;
      Analysis* analysis = EnumToAnalysis(analysis_enum);
      analysis->CreateNodes(new_nodes_list[i],iomodel,true);
      delete analysis;
      this->UpdateElements(newnumberofelements,newelementslist,iomodel->my_elements,i,new_elements);
      this->CreateConstraints(new_vertices,analysis_enum,new_constraints_list[i]);
 
      new_constraints_list[i]->Presort();
      new_nodes_list[i]->Presort();
   }
 
   new_elements->Presort();
   new_vertices->Presort();
   //this->loads->Presort();
   new_materials->Presort();
 
   /*reset hooks*/
   new_elements->ResetHooks();
   //this->loads->ResetHooks();
   new_materials->ResetHooks();
 
   /*do the post-processing of the datasets to get an FemModel that can actually run analyses: */
   int analysis_type;
   for(int i=0;i<this->nummodels;i++){
      analysis_type=this->analysis_type_list[i];
      SetCurrentConfiguration(analysis_type);
 
      this->analysis_counter=i;
      /*Now, plug analysis_counter and analysis_type inside the parameters: */
      this->parameters->SetParam(this->analysis_counter,AnalysisCounterEnum);
      this->parameters->SetParam(analysis_type,AnalysisTypeEnum);
      this->parameters->SetParam(analysis_type,ConfigurationTypeEnum);
 
      /*configure elements, loads and nodes, for this new analysis: */
      new_elements->SetCurrentConfiguration(new_elements,this->loads,new_nodes_list[i],new_vertices,new_materials,this->parameters);
      this->loads->SetCurrentConfiguration(new_elements,this->loads,new_nodes_list[i],new_vertices,new_materials,this->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");
      }
 
      ConfigureObjectsx(new_elements,this->loads,new_nodes_list[i],new_vertices,new_materials,this->parameters,new_inputs2);
      SpcNodesx(new_nodes_list[i],new_constraints_list[i],this->parameters);
      NodesDofx(new_nodes_list[i],this->parameters);
   }
 
   /*Interpolate all inputs and insert them into the new elements.*/
   this->InterpolateInputs(new_vertices,new_elements,new_inputs2);
 
   /*Delete old structure and set new pointers*/
   delete this->inputs2;   this->inputs2 = new_inputs2;
   delete this->vertices;  this->vertices = new_vertices;
   delete this->elements;  this->elements = new_elements;
   delete this->materials; this->materials = new_materials;
   if(this->constraints_list && this->nummodels){
      for(int i=0;i<this->nummodels;i++) delete this->constraints_list[i];
      xDelete<Constraints*>(this->constraints_list);
   }
   this->constraints_list= new_constraints_list;
   if(this->nodes_list && this->nummodels){
      for(int i=0;i<this->nummodels;i++) delete this->nodes_list[i];
      xDelete<Nodes*>(this->nodes_list);
   }
   this->nodes_list = new_nodes_list;
 
   /*Reset mask*/
   GetMaskOfIceVerticesLSMx0(this);
 
   /*Insert MISMIP+ bed topography FIXME it could be stay in another place*/
   this->parameters->FindParam(&basalforcing_model,BasalforcingsEnum);
   if(basalforcing_model==MismipFloatingMeltRateEnum) this->BedrockFromMismipPlus();
 
   /*Adjust base, thickness and mask grounded ice leve set*/
   this->parameters->FindParam(&isgroundingline,TransientIsgroundinglineEnum);
   if(isgroundingline) this->AdjustBaseThicknessAndMask();
 
   /*Reset current configuration: */
   analysis_type=this->analysis_type_list[this->analysis_counter];
   SetCurrentConfiguration(analysis_type);
 
   /*Set the new mesh*/
   this->SetMesh(&newelementslist,&newx,&newy,&newnumberofvertices,&newnumberofelements);
 
   /*Cleanup*/
   xDelete<IssmDouble>(newz);
   /*Delete iomodel, but make sure to not erase some pointers*/
   iomodel->elements = NULL;
   delete iomodel;
}
/*}}}*/
void FemModel::BedrockFromMismipPlus(void){/*{{{*/
 
   /*Insert bedrock from mismip+ setup*/
   /*This was used to Misomip project/simulations*/
 
   if(VerboseSolution())_printf0_(" call Mismip bedrock adjust module\n");
 
   IssmDouble x,y,bx,by;
   int numvertices      = this->GetElementsWidth();
   IssmDouble* xyz_list = NULL;
   IssmDouble* r        = xNew<IssmDouble>(numvertices);
 
   for(int el=0;el<this->elements->Size();el++){
      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(el));
      element->GetVerticesCoordinates(&xyz_list);
      for(int i=0;i<numvertices;i++){
         x     = xyz_list[3*i+0];
         y     = xyz_list[3*i+1];
         bx    = -150.-728.8*pow(x/300000.,2)+343.91*pow(x/300000.,4)-50.57*pow(x/300000.,6);
         by    = 500./(1.+exp((-2./4000.)*(y-80000./2.-24000.)))+500./(1.+exp((2./4000.)*(y-80000./2.+24000.)));
         r[i]  = max(bx+by,-720.);
      }
      /*insert new bedrock*/
      element->AddInput2(BedEnum,&r[0],P1Enum);
      /*Cleanup*/
      xDelete<IssmDouble>(xyz_list);
   }
   /*Delete*/
   xDelete<IssmDouble>(r);
}
/*}}}*/
void FemModel::AdjustBaseThicknessAndMask(void){/*{{{*/
 
   if(VerboseSolution())_printf0_(" call adjust base and thickness module\n");
 
   int     numvertices = this->GetElementsWidth();
   IssmDouble rho_water,rho_ice,density,base_float;
   IssmDouble* phi     = xNew<IssmDouble>(numvertices);
   IssmDouble* h       = xNew<IssmDouble>(numvertices);
   IssmDouble* s       = xNew<IssmDouble>(numvertices);
   IssmDouble* b       = xNew<IssmDouble>(numvertices);
   IssmDouble* r       = xNew<IssmDouble>(numvertices);
   IssmDouble* sl      = xNew<IssmDouble>(numvertices);
 
   for(int el=0;el<this->elements->Size();el++){
      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(el));
 
      element->GetInputListOnVertices(&s[0],SurfaceEnum);
      element->GetInputListOnVertices(&r[0],BedEnum);
      element->GetInputListOnVertices(&sl[0],SealevelEnum);
      rho_water   = element->FindParam(MaterialsRhoSeawaterEnum);
      rho_ice     = element->FindParam(MaterialsRhoIceEnum);
      density     = rho_ice/rho_water;
 
      for(int i=0;i<numvertices;i++){
         /*calculate base floatation (which supports given surface*/
         base_float = rho_ice*s[i]/(rho_ice-rho_water);
         if(r[i]>base_float){
            b[i] = r[i];
         }
         else {
            b[i] = base_float;
         }
 
         if(fabs(sl[i])>0) _error_("Sea level value "<<sl[i]<<" not supported!");
         /*update thickness and mask grounded ice level set*/
         h[i]    = s[i]-b[i];
         phi[i]  = h[i]+r[i]/density;
      }
 
      /*Update inputs*/
      element->AddInput2(MaskOceanLevelsetEnum,&phi[0],P1Enum);
      element->AddInput2(ThicknessEnum,&h[0],P1Enum);
      element->AddInput2(BaseEnum,&b[0],P1Enum);
   }
 
   /*Delete*/
   xDelete<IssmDouble>(phi);
   xDelete<IssmDouble>(h);
   xDelete<IssmDouble>(s);
   xDelete<IssmDouble>(b);
   xDelete<IssmDouble>(r);
   xDelete<IssmDouble>(sl);
}
/*}}}*/
void FemModel::GetInputs(int* pnumP0inputs,IssmDouble** pP0inputs,int** pP0input_enums,int** pP0input_interp,int* pnumP1inputs,IssmDouble** pP1inputs,int** pP1input_enums,int** pP1input_interp){/*{{{*/
 
   int numberofvertices = this->vertices->NumberOfVertices();
   int numberofelements = this->elements->NumberOfElements();
   int elementswidth    = this->GetElementsWidth();
   int numinputs,numP0inputs,numP1inputs;
   IssmDouble* P0inputs                      = NULL;
   Vector<IssmDouble>* vP0inputs             = NULL;
   int* P0input_enums                        = NULL;
   int* P0input_interp                       = NULL;
   IssmDouble* P1inputs                      = NULL;
   Vector<IssmDouble>* vP1inputs             = NULL;
   int* P1input_enums                        = NULL;
   int* P1input_interp                       = NULL;
   int* input_interpolations                 = NULL;
   int* input_enums                          = NULL;
   int* pos                                  = NULL;
   IssmDouble value                          = 0;
 
   /*Figure out how many inputs we have and their respective interpolation*/
   this->inputs2->GetInputsInterpolations(&numinputs,&input_interpolations,&input_enums);
 
   /*Count and get enums of all inputs in old mesh*/
   for(int step=0;step<2;step++){
      if(step){
         P0input_enums  = xNew<int>(numP0inputs);
         P0input_interp = xNew<int>(numP0inputs);
         P1input_enums  = xNew<int>(numP1inputs);
         P1input_interp = xNew<int>(numP1inputs);
      }
      numP0inputs = 0;
      numP1inputs = 0;
      for(int i=0;i<numinputs;i++){
         int inputinterp = input_interpolations[i];
         switch(inputinterp){
            case 0:
               /*Input not found, go to the next*/
               break;
            case P1Enum:
               if(step){
                  P1input_enums[numP1inputs]  = input_enums[i];
                  P1input_interp[numP1inputs] = inputinterp;
               }
               numP1inputs++;
               break;
            case P0Enum:
            case IntInput2Enum:
            case BoolInput2Enum:
               if(step){
                  P0input_enums[numP0inputs]  = input_enums[i];
                  P0input_interp[numP0inputs] = inputinterp;
               }
               numP0inputs++;
               break;
            default:
               _error_(EnumToStringx(inputinterp)<<" ("<<inputinterp<<") Not supported yet");
         }
      }
   }
 
   /*Get P0 and P1 inputs over the elements*/
   pos      = xNew<int>(elementswidth);
   vP0inputs= new Vector<IssmDouble>(numberofelements*numP0inputs);
   vP1inputs= new Vector<IssmDouble>(numberofvertices*numP1inputs);
   for(int i=0;i<this->elements->Size();i++){
      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
 
      /*Get P0 inputs*/
      for(int j=0;j<numP0inputs;j++){
         switch(P0input_interp[j]){
            case P0Enum:{
               Input2* input=element->GetInput2(P0input_enums[j]);
               input->GetInputAverage(&value);
            }
                      break;
            case IntInput2Enum:{
               int valueint;
               element->GetInput2Value(&valueint,P0input_enums[j]);
               value = reCast<IssmDouble>(valueint);
            }
                            break;
            case BoolInput2Enum:{
               bool valuebool;
               element->GetInput2Value(&valuebool,P0input_enums[j]);
               value = reCast<IssmDouble>(valuebool);
            }
               break;
            default:
               _error_(EnumToStringx(P0input_interp[j])<<" ("<<P0input_interp[j]<<") Not supported yet");
         }
         pos[0]=element->Sid()*numP0inputs+j;
         /*Insert input in the vector*/
         vP0inputs->SetValues(1,pos,&value,INS_VAL);
      }
 
      /*Get P1 inputs*/
      for(int j=0;j<numP1inputs;j++){
         Input2* temp = element->GetInput2(P1input_enums[j]); _assert_(temp);
         ElementInput2* input=xDynamicCast<ElementInput2*>(temp);
         pos[0]=element->vertices[0]->Sid()*numP1inputs+j;
         pos[1]=element->vertices[1]->Sid()*numP1inputs+j;
         pos[2]=element->vertices[2]->Sid()*numP1inputs+j;
         /*Insert input in the vector*/
         vP1inputs->SetValues(elementswidth,pos,input->element_values,INS_VAL);
      }
   }
 
   /*Assemble and serialize*/
   vP0inputs->Assemble();
   vP1inputs->Assemble();
   P0inputs=vP0inputs->ToMPISerial();
   P1inputs=vP1inputs->ToMPISerial();
 
   /*Assign pointers*/
   *pnumP0inputs     = numP0inputs;
   *pP0inputs        = P0inputs;
   *pP0input_enums   = P0input_enums;
   *pP0input_interp  = P0input_interp;
   *pnumP1inputs     = numP1inputs;
   *pP1inputs        = P1inputs;
   *pP1input_enums   = P1input_enums;
   *pP1input_interp  = P1input_interp;
 
   /*Cleanup*/
   delete vP0inputs;
   delete vP1inputs;
   xDelete<int>(input_interpolations);
   xDelete<int>(input_enums);
   xDelete<int>(pos);
}
/*}}}*/
void FemModel::InterpolateInputs(Vertices* newfemmodel_vertices,Elements* newfemmodel_elements,Inputs2* newinputs2){/*{{{*/
 
   int numberofelements       = -1;                                  //global, entire old mesh
   int newnumberofelements    = newfemmodel_elements->Size();        //just on the new partition
   int numberofvertices       = -1;                                  //global, entire old mesh
   int newnumberofvertices    = newfemmodel_vertices->Size();        //just on the new partition
   int elementswidth          = this->GetElementsWidth(); //just tria in this version
   int numP0inputs            = -1;
   IssmDouble* P0inputs       = NULL; //global, entire old mesh
   IssmDouble* newP0inputs    = NULL; //just on the new partition
   int* P0input_enums         = NULL;
   int* P0input_interp        = NULL;
   int numP1inputs            = -1;
   IssmDouble* P1inputs       = NULL; //global, entire old mesh
   IssmDouble* newP1inputs    = NULL; //just on the new partition
   int* P1input_enums         = NULL;
   int* P1input_interp        = NULL;
   IssmDouble* values         = NULL;
   IssmDouble* vector         = NULL;
   IssmDouble* x              = NULL;//global, entire old mesh
   IssmDouble* y              = NULL;//global, entire old mesh
   int* elementslist          = NULL;//global, entire old mesh
   IssmDouble* newx           = NULL;//just on the new partition
   IssmDouble* newy           = NULL;//just on the new partition
   IssmDouble* newz           = NULL;//just on the new partition
   IssmDouble* newxc          = NULL;//just on the new partition
   IssmDouble* newyc          = NULL;//just on the new partition
   int* newelementslist       = NULL;//just on the new partition
   int* sidtoindex            = NULL;//global vertices sid to partition index
 
   /*Get old P0 and P1  inputs (entire mesh)*/
   this->GetInputs(&numP0inputs,&P0inputs,&P0input_enums,&P0input_interp,&numP1inputs,&P1inputs,&P1input_enums,&P1input_interp);
 
   /*Get the old mesh (global, entire mesh)*/
   this->GetMesh(&elementslist,&x,&y,&numberofvertices,&numberofelements);
 
   /*Get the new mesh (just on the new partition)*/
   this->GetMeshOnPartition(newfemmodel_vertices,newfemmodel_elements,&newx,&newy,&newz,&newelementslist,&sidtoindex);
 
   /*Calculate the center points xc and xy (new mesh, new partition)*/
   newxc=xNewZeroInit<IssmDouble>(newnumberofelements);
   newyc=xNewZeroInit<IssmDouble>(newnumberofelements);
   for(int i=0;i<newnumberofelements;i++){
      for(int j=0;j<elementswidth;j++){
         int vid = newelementslist[i*elementswidth+j]-1;//Transform to C indexing
         newxc[i]+=newx[vid]/elementswidth;
         newyc[i]+=newy[vid]/elementswidth;
      }
   }
 
   /*Interplate P0 inputs in the new mesh (just on the new partition)*/
   InterpFromMeshToMesh2dx(&newP0inputs,elementslist,x,y,numberofvertices,numberofelements,
            P0inputs,numberofelements,numP0inputs,
            newxc,newyc,newnumberofelements,NULL);
 
   /*Interpolate P1 inputs in the new mesh (just on the new partition)*/
   InterpFromMeshToMesh2dx(&newP1inputs,elementslist,x,y,numberofvertices,numberofelements,
            P1inputs,numberofvertices,numP1inputs,
            newx,newy,newnumberofvertices,NULL);
 
   /*Insert P0 and P1 inputs into the new elements (just on the new partition)*/
   int vertexlids[3];
   values=xNew<IssmDouble>(elementswidth);
   for(int i=0;i<newfemmodel_elements->Size();i++){//just on the new partition
      Element* element=xDynamicCast<Element*>(newfemmodel_elements->GetObjectByOffset(i));
      /*newP0inputs is just on the new partition*/
      for(int j=0;j<numP0inputs;j++){
         switch(P0input_interp[j]){
            case P0Enum:
               element->SetElementInput(newinputs2,P0input_enums[j],newP0inputs[i*numP0inputs+j]);
               break;
            case IntInput2Enum:
               element->SetIntInput(newinputs2,P0input_enums[j],reCast<int>(newP0inputs[i*numP0inputs+j]));
               break;
            case BoolInput2Enum:
               element->SetBoolInput(newinputs2,P0input_enums[j],reCast<bool>(newP0inputs[i*numP0inputs+j]));
               break;
            default:
               _error_(EnumToStringx(P0input_interp[j])<<" Not supported yet");
         }
      }
      /*newP1inputs is just on the new partition*/
      for(int i=0;i<3;i++) vertexlids[i]=element->vertices[i]->lid;
      for(int j=0;j<numP1inputs;j++){
         values[0]=newP1inputs[sidtoindex[element->vertices[0]->Sid()]*numP1inputs+j];
         values[1]=newP1inputs[sidtoindex[element->vertices[1]->Sid()]*numP1inputs+j];
         values[2]=newP1inputs[sidtoindex[element->vertices[2]->Sid()]*numP1inputs+j];
         newinputs2->SetTriaInput(P1input_enums[j],P1Enum,3,vertexlids,values);
      }
   }
 
 
   /*Cleanup*/
   xDelete<IssmDouble>(P0inputs);
   xDelete<IssmDouble>(newP0inputs);
   xDelete<int>(P0input_enums);
   xDelete<int>(P0input_interp);
   xDelete<IssmDouble>(P1inputs);
   xDelete<IssmDouble>(newP1inputs);
   xDelete<int>(P1input_enums);
   xDelete<int>(P1input_interp);
   xDelete<IssmDouble>(newx);
   xDelete<IssmDouble>(newy);
   xDelete<IssmDouble>(newz);
   xDelete<IssmDouble>(newxc);
   xDelete<IssmDouble>(newyc);
   xDelete<int>(newelementslist);
   xDelete<int>(sidtoindex);
   xDelete<IssmDouble>(values);
}
/*}}}*/
void FemModel::WriteMeshInResults(void){/*{{{*/
 
   /*Write the erros estimators*/
   this->WriteErrorEstimatorsInResults();
 
   int step             = -1;
   int numberofelements = -1;
   int numberofvertices = -1;
   IssmDouble time      = -1;
   IssmDouble* x        = NULL;
   IssmDouble* y        = NULL;
   int* elementslist    = NULL;
 
   if(!this->elements || !this->vertices || !this->results || !this->parameters) return;
 
   parameters->FindParam(&step,StepEnum);
   parameters->FindParam(&time,TimeEnum);
 
   /*Get mesh. Elementslist comes in Matlab indexing*/
   this->GetMesh(&elementslist,&x,&y,&numberofvertices,&numberofelements);
 
   /*Write mesh in Results*/
   this->results->AddResult(new GenericExternalResult<int*>(this->results->Size()+1,MeshElementsEnum,
               elementslist,numberofelements,this->GetElementsWidth(),step,time));
 
   this->results->AddResult(new GenericExternalResult<IssmDouble*>(this->results->Size()+1,MeshXEnum,
               x,numberofvertices,1,step,time));
 
   this->results->AddResult(new GenericExternalResult<IssmDouble*>(this->results->Size()+1,MeshYEnum,
               y,numberofvertices,1,step,time));
}
/*}}}*/
void FemModel::WriteErrorEstimatorsInResults(void){/*{{{*/
 
   int step                   = -1;
   int numberofelements       = -1;
   IssmDouble time            = -1;
   IssmDouble* stresserror    = NULL;
   IssmDouble* thicknesserror = NULL;
 
   if(!this->elements || !this->vertices || !this->results || !this->parameters) return;
 
   parameters->FindParam(&step,StepEnum);
   parameters->FindParam(&time,TimeEnum);
   numberofelements=this->elements->NumberOfElements();
 
   /*Compute the deviatoric stress tensor*/
   this->ZZErrorEstimator(&stresserror);
 
   /*Compute the thickness error*/
   this->ThicknessZZErrorEstimator(&thicknesserror);
 
   /*Write error estimators in Results*/
   this->results->AddResult(new GenericExternalResult<IssmDouble*>(this->results->Size()+1,DeviatoricStressErrorEstimatorEnum,
                                                                  stresserror,numberofelements,1,step,time));
 
   this->results->AddResult(new GenericExternalResult<IssmDouble*>(this->results->Size()+1,ThicknessErrorEstimatorEnum,
                                                                  thicknesserror,numberofelements,1,step,time));
   /*Cleanup*/
   xDelete<IssmDouble>(stresserror);
   xDelete<IssmDouble>(thicknesserror);
 
   return;
}
/*}}}*/
void FemModel::CreateElements(int newnumberofelements,int elementswidth,int* newelementslist,bool* my_elements,Elements* elements){/*{{{*/
 
   /*newlementslist is in Matlab indexing*/
   int lid=0;
   for(int i=0;i<newnumberofelements;i++){
      if(my_elements[i]){
         /*Create element - just tria in this version*/
         Tria *newtria=new Tria();
         newtria->id=i+1;
         newtria->sid=i;
         newtria->lid=lid++;
         newtria->iscollapsed=0;
         newtria->isonsurface = true;
         newtria->isonbase = true;
         newtria->parameters=NULL;
         newtria->inputs2=NULL;
         newtria->nodes=NULL;
         newtria->vertices=NULL;
         newtria->material=NULL;
         if(this->nummodels>0){
            newtria->element_type_list=xNew<int>(this->nummodels);
            for(int j=0;j<nummodels;j++) newtria->element_type_list[j]=0;
         }
         else newtria->element_type_list=NULL;
 
         /*Element hook*/
         int material_id=i+1; // retrieve material_id = i+1;
         /*retrieve vertices ids*/
         int* vertex_ids=xNew<int>(elementswidth);
         for(int j=0;j<elementswidth;j++) vertex_ids[j]=reCast<int>(newelementslist[elementswidth*i+j]);//this Hook wants Matlab indexing
         /*Setting the hooks*/
         newtria->numanalyses =this->nummodels;
         newtria->hnodes      =new Hook*[this->nummodels];
         newtria->hvertices   =new Hook(&vertex_ids[0],elementswidth);
         newtria->hmaterial   =new Hook(&material_id,1);
         newtria->hneighbors  =NULL;
         /*Initialize hnodes as NULL*/
         for(int j=0;j<this->nummodels;j++) newtria->hnodes[j]=NULL;
         /*Clean up*/
         xDelete<int>(vertex_ids);
         elements->AddObject(newtria);
      }
   }
 
}
/*}}}*/
void FemModel::CreateMaterials(int newnumberofelements,bool* my_elements,Materials* materials){/*{{{*/
 
   /*Just Matice in this version*/
   for(int i=0;i<newnumberofelements;i++){
      if(my_elements[i]){
         materials->AddObject(new Matice(i+1,i,MaticeEnum));
      }
   }
}
/*}}}*/
void FemModel::GetMesh(Vertices* femmodel_vertices, Elements* femmodel_elements,IssmDouble** px, IssmDouble** py, int** pelementslist){/*{{{*/
 
   if(!femmodel_vertices) _error_("GetMesh: vertices are NULL.");
   if(!femmodel_elements) _error_("GetMesh: elements are NULL.");
 
   int numberofvertices = femmodel_vertices->NumberOfVertices();
   int numberofelements = femmodel_elements->NumberOfElements();
   int elementswidth    = this->GetElementsWidth(); // just 2D mesh in this version (just tria elements)
   IssmDouble* x        = NULL;
   IssmDouble* y        = NULL;
   IssmDouble* z        = NULL;
   int* elementslist    = NULL;
   int* elem_vertices   = NULL;
   IssmDouble *id1      = NULL;
   IssmDouble *id2      = NULL;
   IssmDouble *id3      = NULL;
 
   /*Get vertices coordinates*/
   VertexCoordinatesx(&x,&y,&z,femmodel_vertices,false) ;
 
   /*Get element vertices*/
   elem_vertices           = xNew<int>(elementswidth);
   Vector<IssmDouble>* vid1= new Vector<IssmDouble>(numberofelements);
   Vector<IssmDouble>* vid2= new Vector<IssmDouble>(numberofelements);
   Vector<IssmDouble>* vid3= new Vector<IssmDouble>(numberofelements);
 
   /*Go through elements, and for each element, get vertices*/
   for(int i=0;i<femmodel_elements->Size();i++){
      Element* element=xDynamicCast<Element*>(femmodel_elements->GetObjectByOffset(i));
      element->GetVerticesSidList(elem_vertices);
      vid1->SetValue(element->sid,elem_vertices[0],INS_VAL);
      vid2->SetValue(element->sid,elem_vertices[1],INS_VAL);
      vid3->SetValue(element->sid,elem_vertices[2],INS_VAL);
   }
 
   /*Assemble*/
   vid1->Assemble();
   vid2->Assemble();
   vid3->Assemble();
 
   /*Serialize*/
   id1 = vid1->ToMPISerial();
   id2 = vid2->ToMPISerial();
   id3 = vid3->ToMPISerial();
 
   /*Construct elements list*/
   elementslist=xNew<int>(numberofelements*elementswidth);
   if(numberofelements*elementswidth<0) _error_("numberofelements negative.");
   for(int i=0;i<numberofelements;i++){
      elementslist[elementswidth*i+0] = reCast<int>(id1[i])+1; //InterpMesh wants Matlab indexing
      elementslist[elementswidth*i+1] = reCast<int>(id2[i])+1; //InterpMesh wants Matlab indexing
      elementslist[elementswidth*i+2] = reCast<int>(id3[i])+1; //InterpMesh wants Matlab indexing
   }
 
   /*Assign pointers*/
   *px            = x;
   *py            = y;
   *pelementslist = elementslist; //Matlab indexing. InterMesh uses this type.
 
   /*Cleanup*/
   xDelete<int>(elem_vertices);
   xDelete<IssmDouble>(id1);
   xDelete<IssmDouble>(id2);
   xDelete<IssmDouble>(id3);
   xDelete<IssmDouble>(z);
   delete vid1;
   delete vid2;
   delete vid3;
}
/*}}}*/
void FemModel::GetMesh(int** elementslist, IssmDouble** x, IssmDouble** y, int* numberofvertices, int* numberofelements){/*{{{*/
 
   int amrtype;
   this->parameters->FindParam(&amrtype,AmrTypeEnum);
 
   switch(amrtype){
 
      #if defined(_HAVE_NEOPZ_)
      case AmrNeopzEnum: this->amr->GetMesh(elementslist,x,y,numberofvertices,numberofelements); break;
      #endif
 
      #if defined(_HAVE_BAMG_)
      case AmrBamgEnum: this->amrbamg->GetMesh(elementslist,x,y,numberofvertices,numberofelements); break;
      #endif
 
      default: _error_("not implemented yet");
   }
}/*}}}*/
void FemModel::SetMesh(int** elementslist, IssmDouble** x, IssmDouble** y, int* numberofvertices, int* numberofelements){/*{{{*/
 
   int amrtype;
   this->parameters->FindParam(&amrtype,AmrTypeEnum);
 
   switch(amrtype){
 
      #if defined(_HAVE_NEOPZ_)
      case AmrNeopzEnum: this->amr->SetMesh(elementslist,x,y,numberofvertices,numberofelements); break;
      #endif
 
      #if defined(_HAVE_BAMG_)
      case AmrBamgEnum: this->amrbamg->SetMesh(elementslist,x,y,numberofvertices,numberofelements); break;
      #endif
 
      default: _error_("not implemented yet");
   }
}/*}}}*/
void FemModel::GetMeshOnPartition(Vertices* femmodel_vertices,Elements* femmodel_elements,IssmDouble** px,IssmDouble** py,IssmDouble** pz,int** pelementslist,int** psidtoindex){/*{{{*/
 
   if(!femmodel_vertices) _error_("GetMesh: vertices are NULL.");
   if(!femmodel_elements) _error_("GetMesh: elements are NULL.");
 
   int numberofvertices       = femmodel_vertices->Size();  //number of vertices of this partition
   int numbertotalofvertices  = femmodel_vertices->NumberOfVertices();  //number total of vertices (entire mesh)
   int numberofelements       = femmodel_elements->Size();  //number of elements of this partition
   int elementswidth          = this->GetElementsWidth();   //just 2D mesh in this version (just tria elements)
   IssmDouble* x              = NULL;
   IssmDouble* y              = NULL;
   IssmDouble* z              = NULL;
   int* elementslist          = NULL;
   int* sidtoindex            = NULL;
   int* elem_vertices         = NULL;
 
   /*Get vertices coordinates of this partition*/
   sidtoindex  = xNewZeroInit<int>(numbertotalofvertices);//entire mesh, all vertices
   x           = xNew<IssmDouble>(numberofvertices);//just this partition
   y           = xNew<IssmDouble>(numberofvertices);//just this partitio;
   z           = xNew<IssmDouble>(numberofvertices);//just this partitio;
 
   /*Go through in this partition (vertices)*/
   for(int i=0;i<numberofvertices;i++){//just this partition
      Vertex* vertex=(Vertex*)femmodel_vertices->GetObjectByOffset(i);
      /*Attention: no spherical coordinates*/
      x[i]=vertex->GetX();
      y[i]=vertex->GetY();
      z[i]=vertex->GetZ();
      /*Keep the index and sid pair*/
      sidtoindex[vertex->Sid()]=i;
   }
 
   /*Go through in this partition (elements) and build the element list*/
   elem_vertices= xNew<int>(elementswidth);
   elementslist = xNew<int>(numberofelements*elementswidth);
   if(numberofelements*elementswidth<0) _error_("numberofelements negative.");
 
   for(int i=0;i<numberofelements;i++){//just this partition
      Element* element=xDynamicCast<Element*>(femmodel_elements->GetObjectByOffset(i));
      element->GetVerticesSidList(elem_vertices);
      elementslist[elementswidth*i+0] = sidtoindex[elem_vertices[0]]+1; //InterpMesh wants Matlab indexing
      elementslist[elementswidth*i+1] = sidtoindex[elem_vertices[1]]+1; //InterpMesh wants Matlab indexing
      elementslist[elementswidth*i+2] = sidtoindex[elem_vertices[2]]+1; //InterpMesh wants Matlab indexing
   }
 
   /*Assign pointers*/
   *px            = x;
   *py            = y;
   *pz            = z;
   *pelementslist = elementslist; //Matlab indexing. InterMesh uses this type.
   *psidtoindex   = sidtoindex;  //it is ncessary to insert inputs
 
   /*Cleanup*/
   xDelete<int>(elem_vertices);
}
/*}}}*/
void FemModel::CreateConstraints(Vertices* newfemmodel_vertices,int analysis_enum,Constraints* newfemmodel_constraints){/*{{{*/
 
   /*ATTENTION: JUST SPCVX AND SPCVY*/
   /*OTHERS CONSTRAINTS MUST BE IMPLEMENTED*/
   if(analysis_enum!=StressbalanceAnalysisEnum) return;
   int analysis_index = AnalysisIndex(analysis_enum);
 
   int numberofnodes_analysistype= this->nodes_list[analysis_index]->NumberOfNodes();
   int dofpernode                = 2;                                         //vx and vy
   int numberofcols              = dofpernode*2;                              //to keep dofs and flags in the vspc vector
   int numberofvertices          = -1;                                        //global, entire old mesh
   int numberofelements          = -1;                                        //global, entire old mesh
   int newnumberofvertices       = newfemmodel_vertices->Size();              //local, just the new partition
   int count                     = 0;
   IssmDouble* x                 = NULL;                                      //global, entire old mesh
   IssmDouble* y                 = NULL;                                      //global, entire old mesh
   int*        elementslist      = NULL;                                      //global, entire old mesh
   IssmDouble* spc               = NULL;                                      //global, entire old mesh
   IssmDouble* newx              = NULL;                                      //local, just new partition
   IssmDouble* newy              = NULL;                                      //local, just new partition
   IssmDouble* newspc            = NULL;                                      //local, just new partition
   IssmDouble eps                = 1.e-8;
   Vector<IssmDouble>* vspc      = new Vector<IssmDouble>(numberofnodes_analysistype*numberofcols);
 
   /*Get old mesh (global, entire mesh). Elementslist comes in Matlab indexing*/
   this->GetMesh(&elementslist,&x,&y,&numberofvertices,&numberofelements);
 
   /*Get vertices coordinates of the new partition*/
   newx=xNew<IssmDouble>(newnumberofvertices);//just the new partition
   newy=xNew<IssmDouble>(newnumberofvertices);//just the new partition
   for(int i=0;i<newnumberofvertices;i++){//just the new partition
      Vertex* vertex=(Vertex*)newfemmodel_vertices->GetObjectByOffset(i);
      /*Attention: no spherical coordinates*/
      newx[i]=vertex->GetX();
      newy[i]=vertex->GetY();
   }
 
   /*Get spcvx and spcvy of old mesh*/
   for(int i=0;i<this->constraints_list[analysis_index]->Size();i++){
 
      Constraint* constraint=(Constraint*)this->constraints_list[analysis_index]->GetObjectByOffset(i);
 
      SpcStatic* spcstatic = xDynamicCast<SpcStatic*>(constraint);
      int dof              = spcstatic->GetDof();
      int node             = spcstatic->GetNodeId();
      IssmDouble spcvalue  = spcstatic->GetValue();
      int nodeindex        = node-1;
 
      /*vx and vx flag insertion*/
      if(dof==0) {//vx
         vspc->SetValue(nodeindex*numberofcols,spcvalue,INS_VAL);    //vx
         vspc->SetValue(nodeindex*numberofcols+dofpernode,1,INS_VAL);//vxflag
      }
      /*vy and vy flag insertion*/
      if(dof==1){//vy
         vspc->SetValue(nodeindex*numberofcols+1,spcvalue,INS_VAL);  //vy
         vspc->SetValue(nodeindex*numberofcols+dofpernode+1,1,INS_VAL);//vyflag
      }
   }
 
   /*Assemble and serialize*/
   vspc->Assemble();
   spc=vspc->ToMPISerial();
 
   /*Interpolate spc values and flags in the new partition*/
   InterpFromMeshToMesh2dx(&newspc,elementslist,x,y,numberofvertices,numberofelements,
                        spc,numberofvertices,numberofcols,
                        newx,newy,newnumberofvertices,NULL);
 
   /*Now, insert the interpolated constraints in the data set (constraints)*/
   count=0;
   for(int i=0;i<newnumberofvertices;i++){//just in the new partition
      Vertex* vertex=(Vertex*)newfemmodel_vertices->GetObjectByOffset(i);
      /*spcvx*/
      if(!xIsNan<IssmDouble>(newspc[i*numberofcols]) && newspc[i*numberofcols+dofpernode]>(1-eps)){
         newfemmodel_constraints->AddObject(new SpcStatic(count+1,vertex->Sid()+1,0,newspc[i*numberofcols],analysis_enum));
         //add count'th spc, on node i+1, setting dof 1 to vx.
         count++;
      }
   }
   count=0;
   for(int i=0;i<newnumberofvertices;i++){//just in the new partition
      Vertex* vertex=(Vertex*)newfemmodel_vertices->GetObjectByOffset(i);
      /*spcvy*/
      if(!xIsNan<IssmDouble>(newspc[i*numberofcols+1]) && newspc[i*numberofcols+dofpernode+1]>(1-eps) ){
         newfemmodel_constraints->AddObject(new SpcStatic(count+1,vertex->Sid()+1,1,newspc[i*numberofcols+1],analysis_enum));
         //add count'th spc, on node i+1, setting dof 1 to vx.
         count++;
      }
   }
 
   /*Cleanup*/
   xDelete<IssmDouble>(spc);
   xDelete<IssmDouble>(newspc);
   xDelete<IssmDouble>(newx);
   xDelete<IssmDouble>(newy);
   delete vspc;
}
/*}}}*/
void FemModel::UpdateElements(int newnumberofelements,int* newelementslist,bool* my_elements,int analysis_counter,Elements* newelements){/*{{{*/
 
   /*newelementslist is in Matlab indexing*/
 
   /*Update elements, set hnode.
   This code is in all analysis */
   int elemcounter=0;
   for(int iel=0;iel<newnumberofelements;iel++){
      if(my_elements[iel]){
         Tria* tria=(Tria*)newelements->GetObjectByOffset(elemcounter);
         //element update
         tria->element_type_list[analysis_counter]=P1Enum;
         int numnodes=3;
         int* tria_node_ids=xNew<int>(numnodes);
         tria_node_ids[0]=newelementslist[3*iel+0]; //matlab indexing
         tria_node_ids[1]=newelementslist[3*iel+1]; //matlab indexing
         tria_node_ids[2]=newelementslist[3*iel+2]; //matlab indexing
         tria->SetHookNodes(tria_node_ids,numnodes,analysis_counter); tria->nodes=NULL;
         xDelete<int>(tria_node_ids);
         elemcounter++;
      }
   }
   return;
}
/*}}}*/
void FemModel::SmoothedDeviatoricStressTensor(IssmDouble** ptauxx,IssmDouble** ptauyy,IssmDouble** ptauxy){/*{{{*/
 
   int elementswidth                   = this->GetElementsWidth();//just 2D mesh, tria elements
   int numberofvertices                = this->vertices->NumberOfVertices();
   IssmDouble weight                   = 0.;
   IssmDouble* tauxx                   = NULL;
   IssmDouble* tauyy                   = NULL;
   IssmDouble* tauxy                   = NULL;
   IssmDouble* totalweight             = NULL;
   IssmDouble* deviatoricstressxx      = xNew<IssmDouble>(elementswidth);
   IssmDouble* deviatoricstressyy      = xNew<IssmDouble>(elementswidth);
   IssmDouble* deviatoricstressxy      = xNew<IssmDouble>(elementswidth);
   int* elem_vertices                  = xNew<int>(elementswidth);
   Vector<IssmDouble>* vectauxx        = new Vector<IssmDouble>(numberofvertices);
   Vector<IssmDouble>* vectauyy        = new Vector<IssmDouble>(numberofvertices);
   Vector<IssmDouble>* vectauxy        = new Vector<IssmDouble>(numberofvertices);
   Vector<IssmDouble>* vectotalweight  = new Vector<IssmDouble>(numberofvertices);
 
   /*Update the Deviatoric Stress tensor over the elements*/
   this->DeviatoricStressx();
 
   /*Calculate the Smoothed Deviatoric Stress tensor*/
   for(int i=0;i<this->elements->Size();i++){
      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
      element->GetInputListOnVertices(deviatoricstressxx,DeviatoricStressxxEnum);
      element->GetInputListOnVertices(deviatoricstressyy,DeviatoricStressyyEnum);
      element->GetInputListOnVertices(deviatoricstressxy,DeviatoricStressxyEnum);
      element->GetVerticesSidList(elem_vertices);
 
      /*weight to calculate the smoothed deviatoric stress*/
      Tria* triaelement = xDynamicCast<Tria*>(element);
      weight            = triaelement->GetArea();//the tria area is a choice for the weight
 
      /*taux xx*/
      vectauxx->SetValue(elem_vertices[0],weight*deviatoricstressxx[0],ADD_VAL);
      vectauxx->SetValue(elem_vertices[1],weight*deviatoricstressxx[1],ADD_VAL);
      vectauxx->SetValue(elem_vertices[2],weight*deviatoricstressxx[2],ADD_VAL);
      /*tau yy*/
      vectauyy->SetValue(elem_vertices[0],weight*deviatoricstressyy[0],ADD_VAL);
      vectauyy->SetValue(elem_vertices[1],weight*deviatoricstressyy[1],ADD_VAL);
      vectauyy->SetValue(elem_vertices[2],weight*deviatoricstressyy[2],ADD_VAL);
      /*tau xy*/
      vectauxy->SetValue(elem_vertices[0],weight*deviatoricstressxy[0],ADD_VAL);
      vectauxy->SetValue(elem_vertices[1],weight*deviatoricstressxy[1],ADD_VAL);
      vectauxy->SetValue(elem_vertices[2],weight*deviatoricstressxy[2],ADD_VAL);
      /*total weight*/
      vectotalweight->SetValue(elem_vertices[0],weight,ADD_VAL);
      vectotalweight->SetValue(elem_vertices[1],weight,ADD_VAL);
      vectotalweight->SetValue(elem_vertices[2],weight,ADD_VAL);
   }
 
   /*Assemble*/
   vectauxx->Assemble();
   vectauyy->Assemble();
   vectauxy->Assemble();
   vectotalweight->Assemble();
 
   /*Serialize*/
   tauxx       = vectauxx->ToMPISerial();
   tauyy       = vectauyy->ToMPISerial();
   tauxy       = vectauxy->ToMPISerial();
   totalweight = vectotalweight->ToMPISerial();
 
   /*Divide for the total weight*/
   for(int i=0;i<numberofvertices;i++){
      _assert_(totalweight[i]>0);
      tauxx[i] = tauxx[i]/totalweight[i];
      tauyy[i] = tauyy[i]/totalweight[i];
      tauxy[i] = tauxy[i]/totalweight[i];
   }
 
   /*Set output*/
   (*ptauxx) = tauxx;
   (*ptauyy) = tauyy;
   (*ptauxy) = tauxy;
 
   /*Cleanup*/
   delete vectauxx;
   delete vectauyy;
   delete vectauxy;
   delete vectotalweight;
   xDelete<IssmDouble>(deviatoricstressxx);
   xDelete<IssmDouble>(deviatoricstressyy);
   xDelete<IssmDouble>(deviatoricstressxy);
   xDelete<IssmDouble>(totalweight);
   xDelete<int>(elem_vertices);
}
/*}}}*/
void FemModel::ZZErrorEstimator(IssmDouble** pelementerror){/*{{{*/
 
   /*Compute the Zienkiewicz and Zhu (ZZ) error estimator for the deviatoric stress tensor.
    * Ref.: Zienkiewicz and Zhu, A Simple Error Estimator and Adaptive Procedure for Practical Engineering Analysis, Int. J. Numer. Meth. Eng, 1987*/
 
   IssmDouble Jdet,error,ftxx,ftyy,ftxy;
   int sid;
   int numnodes                     = this->GetElementsWidth();//just 2D mesh, tria elements, P1
   int numberofelements             = this->elements->NumberOfElements();
   IssmDouble* xyz_list             = NULL;
   IssmDouble* smoothedtauxx        = NULL;
   IssmDouble* smoothedtauyy        = NULL;
   IssmDouble* smoothedtauxy        = NULL;
   IssmDouble* tauxx                = xNew<IssmDouble>(numnodes);
   IssmDouble* tauyy                = xNew<IssmDouble>(numnodes);
   IssmDouble* tauxy                = xNew<IssmDouble>(numnodes);
   IssmDouble* basis                = xNew<IssmDouble>(numnodes);
   int* elem_vertices               = xNew<int>(numnodes);
   Vector<IssmDouble>* velementerror= new Vector<IssmDouble>(numberofelements);
 
   /*Get smoothed deviatoric stress tensor*/
   this->SmoothedDeviatoricStressTensor(&smoothedtauxx,&smoothedtauyy,&smoothedtauxy);
 
   /*Integrate the error over elements*/
   for(int i=0;i<this->elements->Size();i++){
      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
      element->GetInputListOnVertices(tauxx,DeviatoricStressxxEnum);
      element->GetInputListOnVertices(tauyy,DeviatoricStressyyEnum);
      element->GetInputListOnVertices(tauxy,DeviatoricStressxyEnum);
      element->GetVerticesSidList(elem_vertices);
 
      /*Integrate*/
      element->GetVerticesCoordinates(&xyz_list);
      Gauss* gauss=element->NewGauss(2);
      error=0.;
      for(int ig=gauss->begin();ig<gauss->end();ig++){
         gauss->GaussPoint(ig);
         element->JacobianDeterminant(&Jdet,xyz_list,gauss);
         element->NodalFunctions(basis,gauss);
         ftxx=0;ftyy=0;ftxy=0;
         for(int n=0;n<numnodes;n++) {
            ftxx+=(tauxx[n]-smoothedtauxx[elem_vertices[n]])*basis[n];
            ftyy+=(tauyy[n]-smoothedtauyy[elem_vertices[n]])*basis[n];
            ftxy+=(tauxy[n]-smoothedtauxy[elem_vertices[n]])*basis[n];
         }
         error+=Jdet*gauss->weight*( pow(ftxx,2)+pow(ftyy,2)+pow(ftxy,2) ); //e^2
      }
      /*Set the error in the global vector*/
      sid=element->Sid();
      error = sqrt(error);//sqrt(e^2)
      velementerror->SetValue(sid,error,INS_VAL);
      /*Cleanup intermediaries*/
      xDelete<IssmDouble>(xyz_list);
      delete gauss;
   }
 
   /*Assemble*/
   velementerror->Assemble();
 
   /*Serialize and set output*/
   (*pelementerror)=velementerror->ToMPISerial();
 
   /*Cleanup*/
   xDelete<IssmDouble>(smoothedtauxx);
   xDelete<IssmDouble>(smoothedtauyy);
   xDelete<IssmDouble>(smoothedtauxy);
   xDelete<IssmDouble>(tauxx);
   xDelete<IssmDouble>(tauyy);
   xDelete<IssmDouble>(tauxy);
   xDelete<IssmDouble>(basis);
   xDelete<int>(elem_vertices);
   delete velementerror;
}
/*}}}*/
void FemModel::SmoothedGradThickness(IssmDouble** pdHdx,IssmDouble** pdHdy){/*{{{*/
 
   int elementswidth                   = this->GetElementsWidth();//just 2D mesh, tria elements
   int numberofvertices                = this->vertices->NumberOfVertices();
 
   IssmDouble weight                   = 0.;
   IssmDouble* dHdx                    = NULL;
   IssmDouble* dHdy                    = NULL;
   IssmDouble* totalweight             = NULL;
   IssmDouble* xyz_list                = NULL;
   IssmDouble* H                       = xNew<IssmDouble>(elementswidth);
   IssmDouble* GradH                   = xNew<IssmDouble>(2);
   int* elem_vertices                  = xNew<int>(elementswidth);
   Vector<IssmDouble>* vecdHdx         = new Vector<IssmDouble>(numberofvertices);
   Vector<IssmDouble>* vecdHdy         = new Vector<IssmDouble>(numberofvertices);
   Vector<IssmDouble>* vectotalweight  = new Vector<IssmDouble>(numberofvertices);
 
   for(int i=0;i<this->elements->Size();i++){
      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
      element->GetInputListOnVertices(H,ThicknessEnum);
      element->GetVerticesSidList(elem_vertices);
      element->GetVerticesCoordinates(&xyz_list);
 
      /*Get the gradient of thickness at the center point (in fact, GradH is constante over the element)*/
      Gauss* gauss=element->NewGauss(1);
      gauss->GaussPoint(gauss->begin());
      element->ValueP1DerivativesOnGauss(GradH,H,xyz_list,gauss);
 
      /*weight to calculate the smoothed grad H*/
      Tria* triaelement = xDynamicCast<Tria*>(element);
      weight            = triaelement->GetArea();//the tria area is a choice for the weight
 
      /*dH/dx*/
      vecdHdx->SetValue(elem_vertices[0],weight*GradH[0],ADD_VAL);
      vecdHdx->SetValue(elem_vertices[1],weight*GradH[0],ADD_VAL);
      vecdHdx->SetValue(elem_vertices[2],weight*GradH[0],ADD_VAL);
      /*dH/dy*/
      vecdHdy->SetValue(elem_vertices[0],weight*GradH[1],ADD_VAL);
      vecdHdy->SetValue(elem_vertices[1],weight*GradH[1],ADD_VAL);
      vecdHdy->SetValue(elem_vertices[2],weight*GradH[1],ADD_VAL);
      /*total weight*/
      vectotalweight->SetValue(elem_vertices[0],weight,ADD_VAL);
      vectotalweight->SetValue(elem_vertices[1],weight,ADD_VAL);
      vectotalweight->SetValue(elem_vertices[2],weight,ADD_VAL);
      /*Cleanup intermediaries*/
      xDelete<IssmDouble>(xyz_list);
      delete gauss;
   }
 
   /*Assemble*/
   vecdHdx->Assemble();
   vecdHdy->Assemble();
   vectotalweight->Assemble();
 
   /*Serialize*/
   dHdx        = vecdHdx->ToMPISerial();
   dHdy        = vecdHdy->ToMPISerial();
   totalweight = vectotalweight->ToMPISerial();
 
   /*Divide for the total weight*/
   for(int i=0;i<numberofvertices;i++){
      _assert_(totalweight[i]>0);
      dHdx[i] = dHdx[i]/totalweight[i];
      dHdy[i] = dHdy[i]/totalweight[i];
   }
 
   /*Set output*/
   (*pdHdx) = dHdx;
   (*pdHdy) = dHdy;
 
   /*Cleanup*/
   delete vecdHdx;
   delete vecdHdy;
   delete vectotalweight;
   xDelete<IssmDouble>(H);
   xDelete<IssmDouble>(GradH);
   xDelete<IssmDouble>(totalweight);
   xDelete<int>(elem_vertices);
}
/*}}}*/
void FemModel::ThicknessZZErrorEstimator(IssmDouble** pelementerror){/*{{{*/
   /*Compute the Zienkiewicz and Zhu (ZZ) error estimator for the thickness
    * Ref.: Zienkiewicz and Zhu, A Simple Error Estimator and Adaptive Procedure for Practical Engineering Analysis, Int. J. Numer. Meth. Eng, 1987*/
 
   IssmDouble Jdet,error,fdHdx,fdHdy;
   int sid;
   int numnodes                     = this->GetElementsWidth();//just 2D mesh, tria elements, P1
   int numberofelements             = this->elements->NumberOfElements();
   IssmDouble* xyz_list             = NULL;
   IssmDouble* smoothed_dHdx        = NULL;
   IssmDouble* smoothed_dHdy        = NULL;
   IssmDouble* H                    = xNew<IssmDouble>(numnodes);
   IssmDouble* GradH                = xNew<IssmDouble>(2);
   IssmDouble* basis                = xNew<IssmDouble>(numnodes);
   int* elem_vertices               = xNew<int>(numnodes);
   Vector<IssmDouble>* velementerror= new Vector<IssmDouble>(numberofelements);
 
   /*Get smoothed deviatoric stress tensor*/
   this->SmoothedGradThickness(&smoothed_dHdx,&smoothed_dHdy);
 
   /*Integrate the error over elements*/
   for(int i=0;i<this->elements->Size();i++){
      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
      element->GetInputListOnVertices(H,ThicknessEnum);
      element->GetVerticesSidList(elem_vertices);
      element->GetVerticesCoordinates(&xyz_list);
      /*Get the gradient of thickness*/
      Gauss* gaussH=element->NewGauss(1);
      gaussH->GaussPoint(gaussH->begin());
      element->ValueP1DerivativesOnGauss(GradH,H,xyz_list,gaussH);
      /*Integrate*/
      Gauss* gauss=element->NewGauss(2);
      error=0.;
      for(int ig=gauss->begin();ig<gauss->end();ig++){
         gauss->GaussPoint(ig);
         element->JacobianDeterminant(&Jdet,xyz_list,gauss);
         element->NodalFunctions(basis,gauss);
         fdHdx=0;fdHdy=0;
         for(int n=0;n<numnodes;n++) {
            fdHdx+=(GradH[0]-smoothed_dHdx[elem_vertices[n]])*basis[n];
            fdHdy+=(GradH[1]-smoothed_dHdy[elem_vertices[n]])*basis[n];
         }
         error+=Jdet*gauss->weight*(pow(fdHdx,2)+pow(fdHdy,2) ); //e^2
      }
      /*Set the error in the global vector*/
      sid=element->Sid();
      error = sqrt(error); //sqrt( e^2 )
      velementerror->SetValue(sid,error,INS_VAL);
      /*Cleanup intermediaries*/
      xDelete<IssmDouble>(xyz_list);
      delete gaussH;
      delete gauss;
   }
 
   /*Assemble*/
   velementerror->Assemble();
 
   /*Serialize and set output*/
   (*pelementerror)=velementerror->ToMPISerial();
 
   /*Cleanup*/
   xDelete<IssmDouble>(smoothed_dHdx);
   xDelete<IssmDouble>(smoothed_dHdy);
   xDelete<IssmDouble>(H);
   xDelete<IssmDouble>(GradH);
   xDelete<IssmDouble>(basis);
   xDelete<int>(elem_vertices);
   delete velementerror;
}
/*}}}*/
void FemModel::MeanGroundedIceLevelSet(IssmDouble** pmasklevelset){/*{{{*/
 
   int elementswidth                   = this->GetElementsWidth();
   int numberofelements                = this->elements->NumberOfElements();
   IssmDouble* elementlevelset         = xNew<IssmDouble>(elementswidth);
   Vector<IssmDouble>* vmasklevelset   = new Vector<IssmDouble>(numberofelements);
 
   for(int i=0;i<this->elements->Size();i++){
      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
      element->GetInputListOnVertices(elementlevelset,MaskOceanLevelsetEnum);
      int sid = element->Sid();
      vmasklevelset->SetValue(sid,(elementlevelset[0]+elementlevelset[1]+elementlevelset[2])/3.,INS_VAL);
   }
 
   /*Assemble*/
   vmasklevelset->Assemble();
 
   /*Serialize and set output*/
   (*pmasklevelset)=vmasklevelset->ToMPISerial();
 
   /*Cleanup*/
   xDelete<IssmDouble>(elementlevelset);
   delete vmasklevelset;
}
/*}}}*/
void FemModel::GetElementCenterCoordinates(IssmDouble** pxc,IssmDouble** pyc){/*{{{*/
 
   /*Intermediaries*/
   int elementswidth          = this->GetElementsWidth();
   int numberofelements       = this->elements->NumberOfElements();
   int* elem_vertices         = xNew<int>(elementswidth);
   Vector<IssmDouble>* vxc    = new Vector<IssmDouble>(numberofelements);
   Vector<IssmDouble>* vyc    = new Vector<IssmDouble>(numberofelements);
   IssmDouble* x              = NULL;
   IssmDouble* y              = NULL;
   IssmDouble* z              = NULL;
   IssmDouble* xyz_list       = NULL;
   IssmDouble x1,y1,x2,y2,x3,y3;
 
   /*Insert the element center coordinates*/
   for(int i=0;i<this->elements->Size();i++){
      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
      //element->GetVerticesSidList(elem_vertices);
      int sid = element->Sid();
      element->GetVerticesCoordinates(&xyz_list);
      x1 = xyz_list[3*0+0];y1 = xyz_list[3*0+1];
      x2 = xyz_list[3*1+0];y2 = xyz_list[3*1+1];
      x3 = xyz_list[3*2+0];y3 = xyz_list[3*2+1];
      vxc->SetValue(sid,(x1+x2+x3)/3.,INS_VAL);
      vyc->SetValue(sid,(y1+y2+y3)/3.,INS_VAL);
   }
 
   /*Assemble*/
   vxc->Assemble();
   vyc->Assemble();
 
   /*Serialize and set output*/
   (*pxc)=vxc->ToMPISerial();
   (*pyc)=vyc->ToMPISerial();
 
   /*Cleanup*/
   xDelete<IssmDouble>(x);
   xDelete<IssmDouble>(y);
   xDelete<IssmDouble>(z);
   xDelete<IssmDouble>(xyz_list);
   xDelete<int>(elem_vertices);
   delete vxc;
   delete vyc;
}
/*}}}*/
void FemModel::GetZeroLevelSetPoints(IssmDouble** pzerolevelset_points,int &numberofpoints,int levelset_type){/*{{{*/
 
   /*Here, "zero level set" means grounding line or ice front, depending on the level set type*/
   /*pzerolevelset_points are the element center points with zero level set. X and Y coords*/
   if(levelset_type!=MaskOceanLevelsetEnum && levelset_type!=MaskIceLevelsetEnum){
      _error_("level set type not implemented yet!");
   }
 
   /*Outputs*/
   IssmDouble* zerolevelset_points        = NULL;
   int npoints                            = 0;
 
   /*Intermediaries*/
   int elementswidth                      = this->GetElementsWidth();
   int numberofelements                   = this->elements->NumberOfElements();
   int* elem_vertices                     = xNew<int>(elementswidth);
   IssmDouble* levelset                   = xNew<IssmDouble>(elementswidth);
   IssmDouble* xyz_list                   = NULL;
   Vector<IssmDouble>* vx_zerolevelset    = new Vector<IssmDouble>(numberofelements);
   Vector<IssmDouble>* vy_zerolevelset    = new Vector<IssmDouble>(numberofelements);
   IssmDouble* x_zerolevelset             = NULL;
   IssmDouble* y_zerolevelset             = NULL;
   int count,sid;
   IssmDouble xc,yc,x1,y1,x2,y2,x3,y3;
 
   /*Use the element center coordinate if level set is zero (grounding line or ice front), otherwise set NAN*/
   for(int i=0;i<this->elements->Size();i++){
      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
      element->GetInputListOnVertices(levelset,levelset_type);
      element->GetVerticesSidList(elem_vertices);
      sid= element->Sid();
      element->GetVerticesCoordinates(&xyz_list);
      x1 = xyz_list[3*0+0];y1 = xyz_list[3*0+1];
      x2 = xyz_list[3*1+0];y2 = xyz_list[3*1+1];
      x3 = xyz_list[3*2+0];y3 = xyz_list[3*2+1];
      xc = NAN;
      yc = NAN;
      Tria* tria  = xDynamicCast<Tria*>(element);
      if(tria->IsIceInElement()){/*verify if there is ice in the element*/
         if(levelset[0]*levelset[1]<0. || levelset[0]*levelset[2]<0. ||
            abs(levelset[0]*levelset[1])<DBL_EPSILON || abs(levelset[0]*levelset[2])<DBL_EPSILON) {
            xc=(x1+x2+x3)/3.;
            yc=(y1+y2+y3)/3.;
         }
      }
      vx_zerolevelset->SetValue(sid,xc,INS_VAL);
      vy_zerolevelset->SetValue(sid,yc,INS_VAL);
      xDelete<IssmDouble>(xyz_list);
   }
   /*Assemble and serialize*/
   vx_zerolevelset->Assemble();
   vy_zerolevelset->Assemble();
   x_zerolevelset=vx_zerolevelset->ToMPISerial();
   y_zerolevelset=vy_zerolevelset->ToMPISerial();
 
   /*Find the number of points*/
   npoints=0;
   for(int i=0;i<numberofelements;i++) if(!xIsNan<IssmDouble>(x_zerolevelset[i])) npoints++;
 
   /*Keep just the element center coordinates with zero level set (compact the structure)*/
   zerolevelset_points=xNew<IssmDouble>(2*npoints);//x and y
   count=0;
   for(int i=0;i<numberofelements;i++){
      if(!xIsNan<IssmDouble>(x_zerolevelset[i])){
         zerolevelset_points[2*count]   = x_zerolevelset[i];
         zerolevelset_points[2*count+1] = y_zerolevelset[i];
         count++;
      }
   }
 
   /*Assign outputs*/
   numberofpoints          = npoints;
   (*pzerolevelset_points) = zerolevelset_points;
 
   /*Cleanup*/
   xDelete<int>(elem_vertices);
   xDelete<IssmDouble>(levelset);
   xDelete<IssmDouble>(x_zerolevelset);
   xDelete<IssmDouble>(y_zerolevelset);
   xDelete<IssmDouble>(xyz_list);
   delete vx_zerolevelset;
   delete vy_zerolevelset;
}
/*}}}*/
#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;
   double femmodel_response;
   int    flag;
   double *vertex_response   = NULL;
   double *qmu_response      = NULL;
   double *responses_pointer = NULL;
 
   IssmDouble **response_partitions         = NULL;
   IssmDouble * response_partition         = NULL;
   int * response_partitions_npart         = NULL;
   int          response_partitions_num;
   int          npart;
 
   /*retrieve partition vectors for responses that are scaled:*/
   this->parameters->FindParam(&response_partitions,&response_partitions_num,NULL,NULL,QmuResponsePartitionsEnum); 
   this->parameters->FindParam(&response_partitions_npart,NULL,NULL,QmuResponsePartitionsNpartEnum); 
   
   /*retrieve my_rank: */
   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),VertexPIdEnum);
 
         /*recover partition vector: */
         response_partition=response_partitions[i];
         npart=response_partitions_npart[i];
 
         /*Now, average it onto the partition nodes: */
         AverageOntoPartitionx(&qmu_response,elements,nodes,vertices,loads,materials,parameters,vertex_response,response_partition,npart);
 
         /*Copy onto our dakota responses: */
         if(my_rank==0){
            /*plug response: */
            for(j=0;j<npart;j++)responses_pointer[j]=qmu_response[j];
 
            /*increment response_pointer :*/
            responses_pointer+=npart;
         }
 
         /*Free ressources:*/
         xDelete<double>(vertex_response);
         xDelete<double>(qmu_response);
 
      }
      else if (flag==IndexedEnum){
 
         /*indexed response: plug index into parameters and call response module: */
         parameters->SetParam(index,IndexEnum);
 
         this->Responsex(&femmodel_response,root);
 
         if(my_rank==0){
            /*plug response: */
            responses_pointer[0]=femmodel_response;
 
            /*increment response_pointer :*/
            responses_pointer++;
         }
      }
      else if (flag==NodalEnum){
         _error_("nodal response functions not supported yet!");
 
         /*increment response_pointer :*/
         responses_pointer++;
      }
      else if (flag==RegularEnum){
 
         /*perfectly normal response function: */
         this->Responsex(&femmodel_response,root);
 
         if(my_rank==0){
            /*plug response: */
            responses_pointer[0]=femmodel_response;
 
            /*increment response_pointer :*/
            responses_pointer++;
         }
      }
      else _error_("flag type " << flag << " not supported yet for response analysis");
   }
 
   /*Synthesize echo: {{{*/
   if(my_rank==0){
      _printf_("   responses: " << d_numresponses << ": ");
      for(i=0;i<d_numresponses-1;i++)_printf_(d_responses[i] << "|");
      _printf_(d_responses[d_numresponses-1]);
      _printf_("\n");
   }
   /*}}}*/
 
   /*Free ressources:*/
   for(i=0;i<response_partitions_num;i++){
      IssmDouble* matrix=response_partitions[i];
      xDelete<IssmDouble>(matrix);
   }
   xDelete<IssmDouble*>(response_partitions);
 
}
/*}}}*/
#endif
#ifdef _HAVE_GIA_
void FemModel::Deflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt, IssmDouble* x, IssmDouble* y){ /*{{{*/
 
   /*Go through elements, and add contribution from each element to the deflection vector wg:*/
   for(int i=0;i<elements->Size();i++){
      Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
      element->GiaDeflection(wg,dwgdt, x,y);
   }
}
/*}}}*/
#endif
#ifdef _HAVE_ESA_
void FemModel::EsaGeodetic2D(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pX, Vector<IssmDouble>* pY, IssmDouble* xx, IssmDouble* yy){/*{{{*/
 
   int         ns,nsmax;
 
   /*Go through elements, and add contribution from each element to the deflection vector wg:*/
   ns = elements->Size();
 
   /*Figure out max of ns: */
   ISSM_MPI_Reduce(&ns,&nsmax,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm());
   ISSM_MPI_Bcast(&nsmax,1,ISSM_MPI_INT,0,IssmComm::GetComm());
 
   /*Call the esa geodetic core: */
   for(int i=0;i<nsmax;i++){
      if(i<ns){
         Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
         element->EsaGeodetic2D(pUp,pNorth,pEast,pX,pY,xx,yy);
      }
      if(i%100==0){
         pUp->Assemble();
         pNorth->Assemble();
         pEast->Assemble();
         pX->Assemble();
         pY->Assemble();
      }
   }
 
   /*One last time: */
   pUp->Assemble();
   pNorth->Assemble();
   pEast->Assemble();
   pX->Assemble();
   pY->Assemble();
 
   /*Free ressources:*/
   xDelete<IssmDouble>(xx);
   xDelete<IssmDouble>(yy);
}
/*}}}*/
void FemModel::EsaGeodetic3D(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){/*{{{*/
 
   int         ns,nsmax;
 
   /*Go through elements, and add contribution from each element to the deflection vector wg:*/
   ns = elements->Size();
 
   /*Figure out max of ns: */
   ISSM_MPI_Reduce(&ns,&nsmax,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm());
   ISSM_MPI_Bcast(&nsmax,1,ISSM_MPI_INT,0,IssmComm::GetComm());
 
   /*Call the esa geodetic core: */
   for(int i=0;i<nsmax;i++){
      if(i<ns){
         Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
         element->EsaGeodetic3D(pUp,pNorth,pEast,latitude,longitude,radius,xx,yy,zz);
      }
      if(i%100==0){
         pUp->Assemble();
         pNorth->Assemble();
         pEast->Assemble();
      }
   }
 
   /*One last time: */
   pUp->Assemble();
   pNorth->Assemble();
   pEast->Assemble();
 
   /*Free ressources:*/
   xDelete<IssmDouble>(latitude);
   xDelete<IssmDouble>(longitude);
   xDelete<IssmDouble>(radius);
   xDelete<IssmDouble>(xx);
   xDelete<IssmDouble>(yy);
   xDelete<IssmDouble>(zz);
}
/*}}}*/
#endif
#ifdef _HAVE_SEALEVELRISE_
void FemModel::SealevelriseEustatic(Vector<IssmDouble>* pRSLgi, IssmDouble* poceanarea, IssmDouble* peustatic, SealevelMasks* masks) { /*{{{*/
 
   /*serialized vectors:*/
   IssmDouble  eustatic       = 0.;
   IssmDouble  eustatic_cpu   = 0.;
   IssmDouble  eustatic_cpu_e = 0.;
   IssmDouble  area      = 0.;
   IssmDouble  oceanarea      = 0.;
   IssmDouble  oceanarea_cpu  = 0.;
 
   /*Initialize temporary vector that will be used to sum eustatic components
    * on all local elements, prior to assembly:*/
   int gsize = this->nodes->NumberOfDofs(GsetEnum);
   IssmDouble* RSLgi=xNewZeroInit<IssmDouble>(gsize);
   int* indices=xNew<int>(gsize);
   for(int i=0;i<gsize;i++) indices[i]=i;
 
   /*First, figure out the area of the ocean, which is needed to compute the eustatic component: */
   for(int i=0;i<elements->Size();i++){
      Element*   element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
      element->GetInput2Value(&area,AreaEnum);
      if (masks->isoceanin[i]) oceanarea_cpu += area;
   }
   ISSM_MPI_Reduce (&oceanarea_cpu,&oceanarea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
   ISSM_MPI_Bcast(&oceanarea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
   _assert_(oceanarea>0.);
 
   /*Call the sea level rise core: */
   for(int i=0;i<elements->Size();i++){
      Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
      element->SealevelriseEustatic(RSLgi,&eustatic_cpu_e,masks, oceanarea);
      eustatic_cpu+=eustatic_cpu_e;
   }
 
   /*Plug values once and assemble: */
   pRSLgi->SetValues(gsize,indices,RSLgi,ADD_VAL);
   pRSLgi->Assemble();
 
   /*Sum all eustatic components from all cpus:*/
   ISSM_MPI_Reduce (&eustatic_cpu,&eustatic,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
   ISSM_MPI_Bcast(&eustatic,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
   _assert_(!xIsNan<IssmDouble>(eustatic));
 
   /*Free ressources:*/
   xDelete<int>(indices);
   xDelete<IssmDouble>(RSLgi);
 
   /*Assign output pointers:*/
   *poceanarea = oceanarea;
   *peustatic  = eustatic;
 
}
/*}}}*/
void FemModel::SealevelriseNonEustatic(Vector<IssmDouble>* pRSLgo, Vector<IssmDouble>* pRSLg_old,  SealevelMasks* masks, bool verboseconvolution){/*{{{*/
 
   /*serialized vectors:*/
   IssmDouble* RSLg_old=NULL;
 
   IssmDouble* RSLgo  = NULL;
   int* indices = NULL;
   int         gsize;
 
   bool computerigid = true;
   bool computeelastic= true;
 
   /*recover computational flags: */
   this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum);
   this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum);
 
   /*Initialize temporary vector that will be used to sum eustatic components on all local elements, prior
    * to assembly:*/
   gsize = this->nodes->NumberOfDofs(GsetEnum);
   RSLgo=xNewZeroInit<IssmDouble>(gsize);
   indices=xNew<int>(gsize); for (int i=0;i<gsize;i++)indices[i]=i;
 
   /*Serialize vectors from previous iteration:*/
   RSLg_old=pRSLg_old->ToMPISerial();
 
   /*Call the sea level rise non-eustatic core only if required: */
   if(computerigid | computeelastic){
      for(int i=0;i<elements->Size();i++){
         Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
         element->SealevelriseNonEustatic(RSLgo,RSLg_old,masks);
      }
   }
   pRSLgo->SetValues(gsize,indices,RSLgo,ADD_VAL);
   pRSLgo->Assemble();
 
   /*Free ressources:*/
   xDelete<int>(indices);
   xDelete<IssmDouble>(RSLgo);
   xDelete<IssmDouble>(RSLg_old);
 
}
/*}}}*/
void FemModel::SealevelriseRotationalFeedback(Vector<IssmDouble>* pRSLgo_rot, Vector<IssmDouble>* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz,  SealevelMasks* masks){/*{{{*/
 
   /*serialized vectors:*/
   bool spherical=true;
   IssmDouble* RSLg_old=NULL;
   IssmDouble* tide_love_h  = NULL;
   IssmDouble* tide_love_k  = NULL;
   IssmDouble* load_love_k  = NULL;
   IssmDouble  tide_love_k2secular;
   IssmDouble  moi_e, moi_p, omega, g;
   IssmDouble  m1, m2, m3;
   IssmDouble  lati, longi, radi, value;
   IssmDouble          *latitude    = NULL;
   IssmDouble          *longitude    = NULL;
   IssmDouble          *radius    = NULL;
 
   /*Serialize vectors from previous iteration:*/
   RSLg_old=pRSLg_old->ToMPISerial();
 
   IssmDouble moi_list[3]={0,0,0};
   IssmDouble moi_list_cpu[3]={0,0,0};
   for(int i=0;i<elements->Size();i++){
      Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
      element->SealevelriseMomentOfInertia(&moi_list[0],RSLg_old,masks );
      moi_list_cpu[0] += moi_list[0];
      moi_list_cpu[1] += moi_list[1];
      moi_list_cpu[2] += moi_list[2];
   }
   ISSM_MPI_Reduce (&moi_list_cpu[0],&moi_list[0],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
   ISSM_MPI_Bcast(&moi_list[0],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
   //
   ISSM_MPI_Reduce (&moi_list_cpu[1],&moi_list[1],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
   ISSM_MPI_Bcast(&moi_list[1],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
   //
   ISSM_MPI_Reduce (&moi_list_cpu[2],&moi_list[2],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
   ISSM_MPI_Bcast(&moi_list[2],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
 
   /*pull out some useful parameters: */
   parameters->FindParam(&load_love_k,NULL,NULL,LoadLoveKEnum);
   parameters->FindParam(&tide_love_h,NULL,NULL,TidalLoveHEnum);
   parameters->FindParam(&tide_love_k,NULL,NULL,TidalLoveKEnum);
   parameters->FindParam(&tide_love_k2secular,TidalLoveK2SecularEnum);
   parameters->FindParam(&moi_e,RotationalEquatorialMoiEnum);
   parameters->FindParam(&moi_p,RotationalPolarMoiEnum);
   parameters->FindParam(&omega,RotationalAngularVelocityEnum);
 
   /*compute perturbation terms for angular velocity vector: */
   m1 = 1/(1-tide_love_k[2]/tide_love_k2secular) * (1+load_love_k[2])/(moi_p-moi_e) * moi_list[0];
   m2 = 1/(1-tide_love_k[2]/tide_love_k2secular) * (1+load_love_k[2])/(moi_p-moi_e) * moi_list[1];
   m3 = -(1+load_love_k[2])/moi_p * moi_list[2];   // term associated with fluid number (3-order-of-magnitude smaller) is negelected
 
   /*recover lat,long and radius vectors from vertices: */
   VertexCoordinatesx(&latitude,&longitude,&radius,this->vertices,spherical);
 
   /* Green's function (1+k_2-h_2/g): checked against Glenn Milne's thesis Chapter 3 (eqs: 3.3-4, 3.10-11)
    * Perturbation terms for angular velocity vector (m1, m2, m3): checked against Mitrovica (2005 Appendix) & Jensen et al (2013 Appendix A3)
    * Sea level rotational feedback: checked against GMD eqs 8-9 (only first order terms, i.e., degree 2 order 0 & 1 considered)
    * all DONE in Geographic coordinates: theta \in [-90,90], lambda \in [-180 180]
    */
   for(int i=0;i<vertices->Size();i++){
      int sid;
      //Vertex* vertex=(Vertex*)vertices->GetObjectByOffset(i);
      Vertex* vertex=xDynamicCast<Vertex*>(vertices->GetObjectByOffset(i));
      sid=vertex->Sid();
 
      lati=latitude[sid]/180*PI; longi=longitude[sid]/180*PI; radi=radius[sid];
 
      /*only first order terms are considered now: */
      value=((1.0+tide_love_k[2]-tide_love_h[2])/9.81)*pow(omega*radi,2.0)*
                  (-m3/6.0 + 0.5*m3*cos(2.0*lati) - 0.5*sin(2.*lati)*(m1*cos(longi)+m2*sin(longi)));
 
      pRSLgo_rot->SetValue(sid,value,INS_VAL); //INS_VAL ensures that you don't add several times
   }
 
   /*Assemble mesh velocity*/
   pRSLgo_rot->Assemble();
 
   /*Assign output pointers:*/
   if(pIxz)*pIxz=moi_list[0];
   if(pIyz)*pIyz=moi_list[1];
   if(pIzz)*pIzz=moi_list[2];
   xDelete<IssmDouble>(latitude);
   xDelete<IssmDouble>(longitude);
   xDelete<IssmDouble>(tide_love_h);
   xDelete<IssmDouble>(tide_love_k);
   xDelete<IssmDouble>(load_love_k);
 
   xDelete<IssmDouble>(radius);
 
   /*Free ressources:*/
   xDelete<IssmDouble>(RSLg_old);
 
}
/*}}}*/
void FemModel::SealevelriseElastic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pRSLg, SealevelMasks* masks){/*{{{*/
 
   /*serialized vectors:*/
   IssmDouble* RSLg=NULL;
 
   IssmDouble* Up  = NULL;
   IssmDouble* North  = NULL;
   IssmDouble* East  = NULL;
   int* indices = NULL;
   int  gsize;
   int  horiz;
 
   /*retrieve parameters:*/
   this->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
 
   /*Serialize vectors from previous iteration:*/
   RSLg=pRSLg->ToMPISerial();
 
   /*Initialize temporary vector that will be used to sum eustatic components on all local elements, prior
    * to assembly:*/
   gsize = this->nodes->NumberOfDofs(GsetEnum);
   Up=xNewZeroInit<IssmDouble>(gsize);
   if(horiz){
      North=xNewZeroInit<IssmDouble>(gsize);
      East=xNewZeroInit<IssmDouble>(gsize);
   }
   indices=xNew<int>(gsize); for (int i=0;i<gsize;i++)indices[i]=i;
 
   /*Call the sea level rise core: */
   for(int i=0;i<elements->Size();i++){
      Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
      element->SealevelriseGeodetic(Up,North,East,RSLg,masks);
   }
 
   pUp->SetValues(gsize,indices,Up,ADD_VAL);
   pUp->Assemble();
   if (horiz){
      pNorth->SetValues(gsize,indices,North,ADD_VAL);
      pNorth->Assemble();
      pEast->SetValues(gsize,indices,East,ADD_VAL);
      pEast->Assemble();
   }
 
   /*Free ressources:*/
   xDelete<IssmDouble>(Up);
   if(horiz){
      xDelete<IssmDouble>(North);
      xDelete<IssmDouble>(East);
   }
   xDelete<int>(indices);
   xDelete<IssmDouble>(RSLg);
}
/*}}}*/
IssmDouble FemModel::SealevelriseOceanAverage(Vector<IssmDouble>* RSLg,SealevelMasks* masks, IssmDouble oceanarea) { /*{{{*/
 
   IssmDouble* RSLg_serial=NULL;
   IssmDouble  oceanvalue,oceanvalue_cpu;
 
   /*Serialize vectors from previous iteration:*/
   RSLg_serial=RSLg->ToMPISerial();
 
   /*Initialize:*/
   oceanvalue_cpu=0;
 
   /*Go through elements, and add contribution from each element and divide by overall ocean area:*/
   for(int i=0;i<elements->Size();i++){
      Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
      oceanvalue_cpu += element->OceanAverage(RSLg_serial,masks);
   }
 
   ISSM_MPI_Reduce (&oceanvalue_cpu,&oceanvalue,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
   ISSM_MPI_Bcast(&oceanvalue,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
 
   /*Free ressources:*/
   xDelete<IssmDouble>(RSLg_serial);
 
   return oceanvalue/oceanarea;
}
/*}}}*/
#endif
void FemModel::HydrologyEPLupdateDomainx(IssmDouble* pEplcount){ /*{{{*/
 
   Vector<IssmDouble> *mask = NULL;
   Vector<IssmDouble> *recurence = NULL;
   Vector<IssmDouble> *active = NULL;
   IssmDouble         *serial_mask       = NULL;
   IssmDouble         *serial_rec        = NULL;
   IssmDouble         *serial_active     = NULL;
   IssmDouble         *old_active        = NULL;
   int                *eplzigzag_counter = NULL;
   int                 eplflip_lock;
 
   HydrologyDCEfficientAnalysis   *effanalysis  = new HydrologyDCEfficientAnalysis();
   HydrologyDCInefficientAnalysis *inefanalysis = new HydrologyDCInefficientAnalysis();
 
   /*Step 1: update mask, the mask might be extended by residual and/or using downstream sediment head*/
   int numnodes = this->nodes->NumberOfNodes();
   mask=new Vector<IssmDouble>(numnodes);
   recurence=new Vector<IssmDouble>(numnodes);
   this->parameters->FindParam(&eplzigzag_counter,NULL,EplZigZagCounterEnum);
   this->parameters->FindParam(&eplflip_lock,HydrologydcEplflipLockEnum);
   GetVectorFromInputsx(&old_active,this,HydrologydcMaskEplactiveNodeEnum,NodeSIdEnum);
 
   for (int i=0;i<elements->Size();i++){
      Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
      effanalysis->HydrologyEPLGetMask(mask,recurence,element);
   }
 
   /*check for changes and increment zigzag counter, change the mask if necessary*/
   recurence->Assemble();
   serial_rec=recurence->ToMPISerial();
   for (int i=0;i<nodes->Size();i++){
      Node* node=xDynamicCast<Node*>(nodes->GetObjectByOffset(i));
      if(serial_rec[node->Sid()]==1.)eplzigzag_counter[node->Lid()] ++;
      if(eplzigzag_counter[node->Lid()]>eplflip_lock && eplflip_lock!=0){
         mask->SetValue(node->Sid(),old_active[node->Sid()],INS_VAL);
      }
   }
   this->parameters->SetParam(eplzigzag_counter,this->nodes->Size(),EplZigZagCounterEnum);
   /*Assemble and serialize*/
   mask->Assemble();
   serial_mask=mask->ToMPISerial();
 
   xDelete<int>(eplzigzag_counter);
   xDelete<IssmDouble>(serial_rec);
   xDelete<IssmDouble>(old_active);
   delete mask;
   delete recurence;
 
   /*Update Mask*/
   InputUpdateFromVectorx(this,serial_mask,HydrologydcMaskEplactiveNodeEnum,NodeSIdEnum);
   xDelete<IssmDouble>(serial_mask);
   inefanalysis->ElementizeEplMask(this);
   /*Step 2: update node activity. If one element is connected to mask=1, all nodes are active*/
   active=new Vector<IssmDouble>(nodes->NumberOfNodes());
   for (int i=0;i<elements->Size();i++){
      Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
      effanalysis->HydrologyEPLGetActive(active,element);
   }
 
   /*Assemble and serialize*/
   active->Assemble();
   serial_active=active->ToMPISerial();
   delete active;
 
 
   /*Update node activation accordingly*/
   int         counter  = 0; //this is probably not acurate but we are only interested in positivity
   for(int i=0;i<elements->Size();i++){
      Element    *element  = xDynamicCast<Element*>(elements->GetObjectByOffset(i));
      int         numnodes = element->GetNumberOfNodes();
      IssmDouble *base     = xNew<IssmDouble>(numnodes);
      element->GetInputListOnNodes(&base[0],BaseEnum);
      for(int in=0;in<numnodes;in++){
         Node* node=element->GetNode(in);
         if(serial_active[node->Sid()]==1.){
            node->Activate();
            if(!node->IsClone()) counter++;
         }
         else{
            node->Deactivate();
            node->ApplyConstraint(0,base[in]);
         }
      }
      xDelete<IssmDouble>(base);
   }
   xDelete<IssmDouble>(serial_active);
   delete effanalysis;
   delete inefanalysis;
   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;
   *pEplcount = counter;
   if(VerboseSolution()) {
      if(counter==0){
         _printf0_("   No nodes are active in EPL layer \n");
      }
      else {
         _printf0_("   Some active nodes in EPL layer \n");
      }
   }
 
   /*Update dof indexings*/
   this->UpdateConstraintsx();
}
/*}}}*/
void FemModel::HydrologyIDSupdateDomainx(IssmDouble* pIDScount){ /*{{{*/
 
   bool                isthermal;
   Vector<IssmDouble>* mask            = NULL;
   Vector<IssmDouble>* active          = NULL;
   IssmDouble*         serial_mask  = NULL;
   IssmDouble*         serial_active   = NULL;
 
   HydrologyDCInefficientAnalysis* inefanalysis =  new HydrologyDCInefficientAnalysis();
   parameters->FindParam(&isthermal,TransientIsthermalEnum);
 
   /*When solving a thermal model we update the thawed nodes*/
   if(isthermal){
      /*Step 1: update mask, the mask correspond to thawed nodes (that have a meltingrate)*/
      mask=new Vector<IssmDouble>(this->nodes->NumberOfNodes());
 
      for (int i=0;i<elements->Size();i++){
         Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
         inefanalysis->HydrologyIDSGetMask(mask,element);
      }
      /*Assemble and serialize*/
      mask->Assemble();
      serial_mask=mask->ToMPISerial();
      delete mask;
   }
   /*for other cases we just grab the mask from the initialisation value*/
   else{
      GetVectorFromInputsx(&serial_mask,this,HydrologydcMaskThawedNodeEnum,NodeSIdEnum);
   }
   /*Update Mask and elementize*/
   InputUpdateFromVectorx(this,serial_mask,HydrologydcMaskThawedNodeEnum,NodeSIdEnum);
   xDelete<IssmDouble>(serial_mask);
   inefanalysis->ElementizeIdsMask(this);
 
   /*get node mask coherent with element mask*/
   active=new Vector<IssmDouble>(nodes->NumberOfNodes());
   for (int i=0;i<elements->Size();i++){
      Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
      inefanalysis->HydrologyIdsGetActive(active,element);
   }
 
   /*Assemble and serialize*/
   active->Assemble();
   serial_active=active->ToMPISerial();
   delete active;
 
   /*Update node activation accordingly*/
   int         counter  = 0; //this is probably not acurate but we are only interested in positivity
   for(int i=0;i<elements->Size();i++){
      Element    *element  = xDynamicCast<Element*>(elements->GetObjectByOffset(i));
      int         numnodes = element->GetNumberOfNodes();
      IssmDouble *base     = xNew<IssmDouble>(numnodes);
 
      element->GetInputListOnNodes(&base[0],BaseEnum);
 
      for(int in=0;in<numnodes;in++){
         Node* node=element->GetNode(in);
         if(serial_active[node->Sid()]==1.){
            node->Activate();
            if(!node->IsClone()) counter++;
         }
         else{
            node->Deactivate();
            node->ApplyConstraint(0,base[in]);
         }
      }
      xDelete<IssmDouble>(base);
   }
   xDelete<IssmDouble>(serial_active);
   delete inefanalysis;
   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;
   *pIDScount = counter;
   if(VerboseSolution()) {
      if(counter==0){
         _printf0_("   No nodes are active in IDS layer \n");
      }
      else {
         _printf0_("   Some active nodes in IDS layer \n");
      }
   }
   /*Update dof indexings*/
   this->UpdateConstraintsx();
 
}
/*}}}*/
void FemModel::UpdateConstraintsL2ProjectionEPLx(IssmDouble* pL2count){ /*{{{*/
 
   Vector<IssmDouble>* active        = NULL;
   IssmDouble*         serial_active = NULL;
   HydrologyDCEfficientAnalysis* effanalysis = new HydrologyDCEfficientAnalysis();
 
   /*update node activity. If one element is connected to mask=1, all nodes are active*/
   this->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum);
   active=new Vector<IssmDouble>(nodes->NumberOfNodes());
   for (int i=0;i<elements->Size();i++){
      Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
      effanalysis->HydrologyEPLGetActive(active,element);
   }
 
   /*Assemble and serialize*/
   active->Assemble();
   serial_active=active->ToMPISerial();
   delete active;
   delete effanalysis;
 
   /*Update node activation accordingly*/
   int counter =0;
   this->SetCurrentConfiguration(L2ProjectionEPLAnalysisEnum);
   for (int i=0;i<nodes->Size();i++){
      Node* node=xDynamicCast<Node*>(nodes->GetObjectByOffset(i));
      if(serial_active[node->Sid()]==1.){
         node->Activate();
         if(!node->IsClone()) 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;
   *pL2count = counter;
   if(VerboseSolution()) _printf0_("   Number of active nodes L2 Projection: "<< counter <<"\n");
}
/*}}}*/
void FemModel::InitTransientInputx(int* transientinput_enum,int numoutputs){ /*{{{*/
 
   for(int i=0;i<numoutputs;i++){
      this->inputs2->DeleteInput(transientinput_enum[i]);
      this->inputs2->SetTransientInput(transientinput_enum[i],NULL,0);
      /*We need to configure this input!*/
      TransientInput2* transientinput = this->inputs2->GetTransientInput(transientinput_enum[i]); _assert_(transientinput);
      transientinput->Configure(this->parameters);
   }
}
/*}}}*/
void FemModel::StackTransientInputx(int* input_enum,int* transientinput_enum,IssmDouble subtime,int numoutputs){ /*{{{*/
 
  for(int i=0;i<numoutputs;i++){
      if(input_enum[i]<0){
         _error_("Can't deal with non enum fields for result Stack");
      }
      else{
         for(int j=0;j<elements->Size();j++){
 
            /*Get the right transient input*/
            Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
            TransientInput2* transientinput = this->inputs2->GetTransientInput(transientinput_enum[i]);
 
            /*Get values and lid list*/
            const int   numvertices = element->GetNumberOfVertices();
            IssmDouble* values=xNew<IssmDouble>(numvertices);
            int        *vertexlids = xNew<int>(numvertices);
            element->GetInputListOnVertices(&values[0],input_enum[i]);   //this is the enum to stack
 
            element->GetVerticesLidList(vertexlids);
 
            switch(element->ObjectEnum()){
               case TriaEnum:  transientinput->AddTriaTimeInput(subtime,numvertices,vertexlids,values,P1Enum); break;
               case PentaEnum: transientinput->AddPentaTimeInput(subtime,numvertices,vertexlids,values,P1Enum); break;
               default: _error_("Not implemented yet");
            }
            xDelete<IssmDouble>(values);
            xDelete<int>(vertexlids);
         }
      }
   }
}
/*}}}*/
void FemModel::AverageTransientInputx(int* transientinput_enum,int* averagedinput_enum,IssmDouble init_time,IssmDouble end_time,int numoutputs, int averaging_method){ /*{{{*/
 
   for(int i=0;i<numoutputs;i++){
      for(int j=0;j<this->elements->Size();j++){
         Element* element = xDynamicCast<Element*>(elements->GetObjectByOffset(j));
         element->CreateInputTimeAverage(transientinput_enum[i],averagedinput_enum[i],init_time,end_time,averaging_method);
      }
   }
}/*}}}*/
#ifdef _HAVE_JAVASCRIPT_
FemModel::FemModel(IssmDouble* buffer, int buffersize, char* toolkits, char* solution, char* modelname,ISSM_MPI_Comm incomm, bool trace){ /*{{{*/
   /*configuration: */
   int  solution_type;
   int  ierr;
 
   /*First things first, store the communicator, and set it as a global variable: */
   IssmComm::SetComm(incomm);
 
   /*Start profiler: */
   this->profiler=new Profiler();
   profiler->Start(TOTAL);
 
   /*From command line arguments, retrieve different filenames needed to create the FemModel: */
   solution_type=StringToEnumx(solution);
 
   /*Create femmodel from input files: */
   profiler->Start(MPROCESSOR);
   this->InitFromBuffers((char*)buffer,buffersize,toolkits, solution_type,trace,NULL);
   profiler->Stop(MPROCESSOR);
 
   /*Save communicator in the parameters dataset: */
   this->parameters->AddObject(new GenericParam<ISSM_MPI_Comm>(incomm,FemModelCommEnum));
 
}
/*}}}*/
void FemModel::CleanUpJs(char** poutput, size_t* psize){/*{{{*/
 
   /*Intermediary*/
   FILE *output_fid;
   GenericParam<char**>* outputbufferparam=NULL;
   GenericParam<size_t*>* outputbuffersizeparam=NULL;
   char** poutputbuffer;
   size_t* poutputbuffersize;
 
   /*Before we delete the profiler, report statistics for this run: */
   profiler->Stop(TOTAL);  //final tagging
   _printf0_("\n");
   _printf0_("   "<<setw(40)<<left<<"FemModel initialization elapsed time:"<<profiler->TotalTime(MPROCESSOR) << "\n");
   _printf0_("   "<<setw(40)<<left<<"Core solution elapsed time:"<<profiler->TotalTime(CORE) << "\n");
   _printf0_("\n");
   _printf0_("   Total elapsed time: "
            <<profiler->TotalTimeModHour(TOTAL)<<" hrs "
            <<profiler->TotalTimeModMin(TOTAL)<<" min "
            <<profiler->TotalTimeModSec(TOTAL)<<" sec"
            );
   _printf0_("\n");
 
   /*Before we close the output file, recover the buffer and size:*/
   outputbufferparam = xDynamicCast<GenericParam<char**>*>(this->parameters->FindParamObject(OutputBufferPointerEnum));
   poutputbuffer=outputbufferparam->GetParameterValue();
   outputbuffersizeparam = xDynamicCast<GenericParam<size_t*>*>(this->parameters->FindParamObject(OutputBufferSizePointerEnum));
   poutputbuffersize=outputbuffersizeparam->GetParameterValue();
 
   /*Assign output values: */
   *poutput=*poutputbuffer;
   *psize=*poutputbuffersize;
}
/*}}}*/
void FemModel::InitFromBuffers(char* buffer, int buffersize, char* toolkits, int in_solution_type, bool trace, IssmPDouble* X){/*{{{*/
 
   /*intermediary*/
   FILE       *IOMODEL = NULL;
   FILE       *toolkitsoptionsfid = NULL;
   FILE       *output_fid = NULL;
   int         my_rank;
   size_t      outputsize;
   char       *outputbuffer;
   const char *rootpath = "";   //needed for Dakota runs only, which we won't do here.
 
   /*recover my_rank:*/
   my_rank=IssmComm::GetRank();
 
   /*Open input file descriptor on cpu 0: */
   if(my_rank==0) IOMODEL = fmemopen((void*)buffer, buffersize, "rb");
 
   /*Open toolkits file descriptor: */
   toolkitsoptionsfid=fmemopen((void*)toolkits, strlen(toolkits)+1, "r");
 
   /*Now, go create FemModel:*/
   this->InitFromFids((char*)rootpath,IOMODEL,toolkitsoptionsfid,in_solution_type,trace,X);
 
   /*Close input file and toolkits file descriptors: */
   if(my_rank==0) fclose(IOMODEL);
   fclose(toolkitsoptionsfid);
 
   /*Open output file once for all and add output file descriptor to parameters*/
   output_fid=open_memstream(&outputbuffer,&outputsize);
   if(output_fid==NULL)_error_("could not initialize output stream");
   this->parameters->SetParam(output_fid,OutputFilePointerEnum);
   this->parameters->AddObject(new GenericParam<char**>(&outputbuffer,OutputBufferPointerEnum));
   this->parameters->AddObject(new GenericParam<size_t*>(&outputsize,OutputBufferSizePointerEnum));
 
}/*}}}*/
#endif
 
#if defined(_HAVE_BAMG_) && !defined(_HAVE_AD_)
void FemModel::ReMeshBamg(int* pnewnumberofvertices,int* pnewnumberofelements,IssmDouble** pnewx,IssmDouble** pnewy,IssmDouble** pnewz,int** pnewelementslist){/*{{{*/
 
   /*Output*/
   IssmDouble *newx        = NULL;
   IssmDouble *newy        = NULL;
   IssmDouble *newz        = NULL;
   IssmDouble *newxylist   = NULL;
   int *newelementslist    = NULL;
   int* newdatalist        = NULL;
   int newnumberofvertices = -1;
   int newnumberofelements = -1;
 
   /*Get Rank*/
   int my_rank = IssmComm::GetRank();
 
   /*Intermediaries*/
   int numberofvertices             = this->vertices->NumberOfVertices();
   IssmDouble* vector_serial        = NULL;
   IssmDouble* hmaxvertices_serial  = NULL;
   Vector<IssmDouble> *vector       = NULL;
 
   /*Get vector to create metric*/
   if(this->amrbamg->fieldenum!=NoneEnum){
      GetVectorFromInputsx(&vector,this,this->amrbamg->fieldenum,VertexSIdEnum);
      vector->Assemble();
      vector_serial = vector->ToMPISerial();
   }
 
   /*Get hmaxVertices to create metric*/
   if(this->amrbamg->groundingline_distance>0||this->amrbamg->icefront_distance>0||
      this->amrbamg->thicknesserror_threshold>0||this->amrbamg->deviatoricerror_threshold>0){
      /*Initialize hmaxvertices with NAN*/
      hmaxvertices_serial=xNew<IssmDouble>(numberofvertices);
      for(int i=0;i<numberofvertices;i++) hmaxvertices_serial[i]=NAN;
      /*Fill hmaxvertices*/
      if(this->amrbamg->thicknesserror_threshold>0)   this->GethmaxVerticesFromEstimators(hmaxvertices_serial,ThicknessErrorEstimatorEnum);
      if(this->amrbamg->deviatoricerror_threshold>0)  this->GethmaxVerticesFromEstimators(hmaxvertices_serial,DeviatoricStressErrorEstimatorEnum);
      if(this->amrbamg->groundingline_distance>0)     this->GethmaxVerticesFromZeroLevelSetDistance(hmaxvertices_serial,MaskOceanLevelsetEnum);
      if(this->amrbamg->icefront_distance>0)          this->GethmaxVerticesFromZeroLevelSetDistance(hmaxvertices_serial,MaskIceLevelsetEnum);
   }
 
   if(my_rank==0){
      this->amrbamg->ExecuteRefinementBamg(vector_serial,hmaxvertices_serial,&newdatalist,&newxylist,&newelementslist);
      if(newdatalist[0]<=0 || newdatalist[1]<=0) _error_("Error in the refinement process.");
   }
 
   /*Send new mesh to others CPU's*/
   if(my_rank) newdatalist=xNew<int>(2);
   ISSM_MPI_Bcast(newdatalist,2,ISSM_MPI_INT,0,IssmComm::GetComm());
   newnumberofvertices=newdatalist[0];
   newnumberofelements=newdatalist[1];
   if(my_rank){
      newxylist      =xNew<IssmDouble>(newnumberofvertices*2);
      newelementslist=xNew<int>(newnumberofelements*this->GetElementsWidth());
   }
   ISSM_MPI_Bcast(newxylist,newnumberofvertices*2,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
   ISSM_MPI_Bcast(newelementslist,newnumberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm());
 
   /*Reorganize the data*/
   newx=xNew<IssmDouble>(newnumberofvertices);
   newy=xNew<IssmDouble>(newnumberofvertices);
   newz=xNewZeroInit<IssmDouble>(newnumberofvertices);
   for(int i=0;i<newnumberofvertices;i++){
      newx[i] = newxylist[2*i];
      newy[i] = newxylist[2*i+1];
   }
 
   /*Assign output pointers*/
   *pnewnumberofvertices = newnumberofvertices;
   *pnewnumberofelements = newnumberofelements;
   *pnewx = newx;
   *pnewy = newy;
   *pnewz = newz;
   *pnewelementslist = newelementslist;
 
   /*Cleanup*/
   xDelete<int>(newdatalist);
   xDelete<IssmDouble>(newxylist);
   xDelete<IssmDouble>(vector_serial);
   xDelete<IssmDouble>(hmaxvertices_serial);
   delete vector;
}
/*}}}*/
void FemModel::InitializeAdaptiveRefinementBamg(void){/*{{{*/
 
   /*Define variables*/
   int numberofvertices      = this->vertices->NumberOfVertices();
   int numberofelements      = this->elements->NumberOfElements();
   IssmDouble* x             = NULL;
   IssmDouble* y             = NULL;
   int* elements             = NULL;
   IssmDouble hmin,hmax,err,gradation;
 
   /*Get rank*/
   int my_rank = IssmComm::GetRank();
 
   /*Initialize field as NULL for now*/
   this->amrbamg = NULL;
 
   /*Get vertices coordinates of the coarse mesh (father mesh)*/
   this->GetMesh(this->vertices,this->elements,&x,&y,&elements);
 
   /*Create bamg data structures for bamg*/
   this->amrbamg = new AmrBamg();
 
   /*Get amr parameters*/
   this->parameters->FindParam(&hmin,AmrHminEnum);
   this->parameters->FindParam(&hmax,AmrHmaxEnum);
   this->parameters->FindParam(&err,AmrErrEnum);
   this->parameters->FindParam(&gradation,AmrGradationEnum);
   this->parameters->FindParam(&this->amrbamg->fieldenum,AmrFieldEnum);
   this->parameters->FindParam(&this->amrbamg->keepmetric,AmrKeepMetricEnum);
   this->parameters->FindParam(&this->amrbamg->groundingline_resolution,AmrGroundingLineResolutionEnum);
   this->parameters->FindParam(&this->amrbamg->groundingline_distance,AmrGroundingLineDistanceEnum);
   this->parameters->FindParam(&this->amrbamg->icefront_resolution,AmrIceFrontResolutionEnum);
   this->parameters->FindParam(&this->amrbamg->icefront_distance,AmrIceFrontDistanceEnum);
   this->parameters->FindParam(&this->amrbamg->thicknesserror_resolution,AmrThicknessErrorResolutionEnum);
   this->parameters->FindParam(&this->amrbamg->thicknesserror_threshold,AmrThicknessErrorThresholdEnum);
   this->parameters->FindParam(&this->amrbamg->thicknesserror_groupthreshold,AmrThicknessErrorGroupThresholdEnum);
   this->parameters->FindParam(&this->amrbamg->thicknesserror_maximum,AmrThicknessErrorMaximumEnum);
   this->parameters->FindParam(&this->amrbamg->deviatoricerror_resolution,AmrDeviatoricErrorResolutionEnum);
   this->parameters->FindParam(&this->amrbamg->deviatoricerror_threshold,AmrDeviatoricErrorThresholdEnum);
   this->parameters->FindParam(&this->amrbamg->deviatoricerror_groupthreshold,AmrDeviatoricErrorGroupThresholdEnum);
   this->parameters->FindParam(&this->amrbamg->deviatoricerror_maximum,AmrDeviatoricErrorMaximumEnum);
   /*Set BamgOpts*/
   this->amrbamg->SetBamgOpts(hmin,hmax,err,gradation);
 
   /*Re-create original mesh and put it in bamg structure (only cpu 0)*/
   this->amrbamg->SetMesh(&elements,&x,&y,&numberofvertices,&numberofelements);
   if(my_rank==0){
      this->amrbamg->Initialize();
   }
}
/*}}}*/
void FemModel::GethmaxVerticesFromZeroLevelSetDistance(IssmDouble* hmaxvertices,int levelset_type){/*{{{*/
 
   if(!hmaxvertices) _error_("hmaxvertices is NULL!\n");
 
   /*Intermediaries*/
   int numberofvertices                   = this->vertices->NumberOfVertices();
   Vector<IssmDouble>* vminvertexdistance = new Vector<IssmDouble>(numberofvertices);
   IssmDouble* pminvertexdistance         = NULL;
   IssmDouble* levelset_points            = NULL;
   IssmDouble x,y;
   IssmDouble threshold,resolution;
   IssmDouble minvertexdistance,distance;
   int sid,numberofpoints;
 
   switch(levelset_type){
      case MaskOceanLevelsetEnum:
         threshold   = this->amrbamg->groundingline_distance;
         resolution  = this->amrbamg->groundingline_resolution;
         break;
      case MaskIceLevelsetEnum:
         threshold   = this->amrbamg->icefront_distance;
         resolution  = this->amrbamg->icefront_resolution;
         break;
      default: _error_("not implemented yet");
   }
 
   /*Get points which level set is zero (center of elements with zero level set)*/
   this->GetZeroLevelSetPoints(&levelset_points,numberofpoints,levelset_type);//levelset_points is serial (global)
 
   for(int i=0;i<this->vertices->Size();i++){//only on this partition
      Vertex* vertex=(Vertex*)this->vertices->GetObjectByOffset(i);
      /*Attention: no spherical coordinates*/
      x  = vertex->GetX();
      y  = vertex->GetY();
      sid= vertex->Sid();
      minvertexdistance=INFINITY;
 
      /*Find the minimum vertex distance*/
      for(int j=0;j<numberofpoints;j++){
         distance=sqrt((x-levelset_points[2*j])*(x-levelset_points[2*j])+(y-levelset_points[2*j+1])*(y-levelset_points[2*j+1]));
         minvertexdistance=min(distance,minvertexdistance);
      }
      /*Now, insert in the vector*/
      vminvertexdistance->SetValue(sid,minvertexdistance,INS_VAL);
   }
   /*Assemble*/
   vminvertexdistance->Assemble();
 
   /*Assign the pointer*/
   pminvertexdistance=vminvertexdistance->ToMPISerial();
 
   /*Fill hmaxVertices*/
   for(int i=0;i<numberofvertices;i++){
      if(pminvertexdistance[i]<threshold){ //hmaxvertices is serial (global)
         if(xIsNan<IssmDouble>(hmaxvertices[i])) hmaxvertices[i]=resolution;
         else hmaxvertices[i]=min(resolution,hmaxvertices[i]);
      }
   }
 
   /*Cleanup*/
   xDelete<IssmDouble>(pminvertexdistance);
   xDelete<IssmDouble>(levelset_points);
   delete vminvertexdistance;
}
/*}}}*/
void FemModel::GethmaxVerticesFromEstimators(IssmDouble* hmaxvertices,int errorestimator_type){/*{{{*/
 
   if(!hmaxvertices) _error_("hmaxvertices is NULL!\n");
 
   /*Intermediaries*/
   int elementswidth                   = this->GetElementsWidth();
   int numberofelements                = -1;
   int numberofvertices                = -1;
   IssmDouble  hmax                    = this->amrbamg->GetBamgOpts()->hmax;
   IssmDouble* maxlength               = NULL;
   IssmDouble* error_vertices          = NULL;
   IssmDouble* error_elements          = NULL;
   IssmDouble* x                       = NULL;
   IssmDouble* y                       = NULL;
   int* index                          = NULL;
   IssmDouble maxerror,threshold,groupthreshold,resolution,length;
   IssmDouble L1,L2,L3;
   int vid,v1,v2,v3;
   bool refine;
 
 
   /*Fill variables*/
   switch(errorestimator_type){
      case ThicknessErrorEstimatorEnum:
         threshold      = this->amrbamg->thicknesserror_threshold;
         groupthreshold = this->amrbamg->thicknesserror_groupthreshold;
         resolution     = this->amrbamg->thicknesserror_resolution;
         maxerror       = this->amrbamg->thicknesserror_maximum;
         this->ThicknessZZErrorEstimator(&error_elements);//error is serial, but the calculation is parallel
         break;
      case DeviatoricStressErrorEstimatorEnum:
         threshold      = this->amrbamg->deviatoricerror_threshold;
         groupthreshold = this->amrbamg->deviatoricerror_groupthreshold;
         resolution     = this->amrbamg->deviatoricerror_resolution;
         maxerror       = this->amrbamg->deviatoricerror_maximum;
         this->ZZErrorEstimator(&error_elements);//error is serial, but the calculation is parallel
         break;
      default: _error_("not implemented yet");
   }
   if(!error_elements) _error_("error_elements is NULL!\n");
   if(groupthreshold<DBL_EPSILON) _error_("group threshold is too small!");
 
   /*Get mesh*/
   this->GetMesh(&index,&x,&y,&numberofvertices,&numberofelements);
   if(numberofelements<0) _error_("number of elements is negative!\n");
   if(numberofvertices<0) _error_("number of vertices is negative!\n");
   maxlength      = xNew<IssmDouble>(numberofelements);
   error_vertices = xNewZeroInit<IssmDouble>(numberofvertices);
 
   /*Find the max of the estimators if it was not provided*/
   if(maxerror<DBL_EPSILON){
      for(int i=0;i<numberofelements;i++) maxerror=max(maxerror,error_elements[i]);
      switch(errorestimator_type){
         case ThicknessErrorEstimatorEnum:         this->amrbamg->thicknesserror_maximum  = maxerror;break;
         case DeviatoricStressErrorEstimatorEnum:  this->amrbamg->deviatoricerror_maximum = maxerror;break;
      }
   }
 
   /*Fill error_vertices (this is the sum of all elements connected to the vertex)*/
   for(int i=0;i<numberofelements;i++){
      v1=index[i*elementswidth+0]-1;//Matlab to C indexing
      v2=index[i*elementswidth+1]-1;//Matlab to C indexing
      v3=index[i*elementswidth+2]-1;//Matlab to C indexing
      L1=sqrt(pow(x[v2]-x[v1],2)+pow(y[v2]-y[v1],2));
      L2=sqrt(pow(x[v3]-x[v2],2)+pow(y[v3]-y[v2],2));
      L3=sqrt(pow(x[v1]-x[v3],2)+pow(y[v1]-y[v3],2));
      /*Fill the vectors*/
      maxlength[i]      =max(L1,max(L2,L3));
      error_vertices[v1]+=error_elements[i];
      error_vertices[v2]+=error_elements[i];
      error_vertices[v3]+=error_elements[i];
   }
 
   /*Fill hmaxvertices with the criteria*/
   for(int i=0;i<numberofelements;i++){
      /*Refine any element if its error > phi*maxerror*/
      if(error_elements[i]>threshold*maxerror){
         /*Now, fill the hmaxvertices if requested*/
         for(int j=0;j<elementswidth;j++){
            vid=index[i*elementswidth+j]-1;//Matlab to C indexing
            if(xIsNan<IssmDouble>(hmaxvertices[vid])) hmaxvertices[vid]=max(maxlength[i]/2.,resolution);//Try first dividing the element
            else hmaxvertices[vid]=min(max(maxlength[i]/2.,resolution),hmaxvertices[vid]);//Try first dividing the element
         }
      }
      else {
         /*Try unrefine the element*/
         if(maxlength[i] < 1.1*hmax/2.){
            for(int j=0;j<elementswidth;j++){
               vid=index[i*elementswidth+j]-1;//Matlab to C indexing
               if(error_vertices[vid]>groupthreshold*maxerror) hmaxvertices[vid]=maxlength[i]; //keep the current resolution
               else{
                  if(xIsNan<IssmDouble>(hmaxvertices[vid])) hmaxvertices[vid]=min(maxlength[i]*2.,hmax);
                  else hmaxvertices[vid]=min(min(maxlength[i]*2.,hmax),hmaxvertices[vid]);//Try first to duplicate the element
               }
            }
         }
      }
   }
 
   /*Cleanup*/
   xDelete<IssmDouble>(error_elements);
   xDelete<IssmDouble>(error_vertices);
   xDelete<IssmDouble>(maxlength);
}
/*}}}*/
void FemModel::GetVerticeDistanceToZeroLevelSet(IssmDouble** pverticedistance,int levelset_type){/*{{{*/
 
   //itapopo esse metodo pode ser deletado
 
 
   /*Here, "zero level set" means grounding line or ice front, depending on the level set type*/
   /*pverticedistance is the minimal vertice distance to the grounding line or ice front*/
   if(levelset_type!=MaskOceanLevelsetEnum && levelset_type!=MaskIceLevelsetEnum){
      _error_("level set type not implemented yet!");
   }
 
   /*Output*/
   IssmDouble* verticedistance;
 
   /*Intermediaries*/
   int numberofvertices       = -1;
   int numberofelements       = -1;
   IssmDouble* levelset_points= NULL;
   IssmDouble* x              = NULL;
   IssmDouble* y              = NULL;
   int* elementslist          = NULL;
   int numberofpoints;
   IssmDouble distance;
 
   /*Get vertices coordinates*/
   this->GetMesh(&elementslist,&x,&y,&numberofvertices,&numberofelements);
   //this->GetMeshOnPartition(this->vertices,this->elements,&x,&y,&z,&elementslist,&sidtoindex);
 
   /*Get points which level set is zero (center of elements with zero level set)*/
   this->GetZeroLevelSetPoints(&levelset_points,numberofpoints,levelset_type);
 
   /*Find the minimal vertice distance to the zero levelset (grounding line or ice front)*/
   verticedistance=xNew<IssmDouble>(numberofvertices);
   for(int i=0;i<numberofvertices;i++){
      verticedistance[i]=INFINITY;
      for(int j=0;j<numberofpoints;j++){
         distance=sqrt((x[i]-levelset_points[2*j])*(x[i]-levelset_points[2*j])+(y[i]-levelset_points[2*j+1])*(y[i]-levelset_points[2*j+1]));
         verticedistance[i]=min(distance,verticedistance[i]);
      }
   }
 
   /*Assign the pointer*/
   (*pverticedistance)=verticedistance;
 
   /*Cleanup*/
   xDelete<IssmDouble>(levelset_points);
}
/*}}}*/
#endif
 
#if defined(_HAVE_NEOPZ_) && !defined(_HAVE_AD_)
void FemModel::ReMeshNeopz(int* pnewnumberofvertices,int* pnewnumberofelements,IssmDouble** pnewx,IssmDouble** pnewy,IssmDouble** pnewz,int** pnewelementslist){/*{{{*/
 
   /*pnewelementslist keep vertices in Matlab indexing*/
   int my_rank                = IssmComm::GetRank();
   IssmDouble* gl_distance    = NULL;
   IssmDouble* if_distance    = NULL;
   IssmDouble* deviatoricerror= NULL;
   IssmDouble* thicknesserror = NULL;
   IssmDouble* newx           = NULL;
   IssmDouble* newy           = NULL;
   IssmDouble* newz           = NULL;
   IssmDouble* newxylist      = NULL;
   int* newelementslist       = NULL;
   int* newdatalist           = NULL;
   int newnumberofvertices    = -1;
   int newnumberofelements    = -1;
 
   /*Get fields, if requested*/
   if(this->amr->groundingline_distance>0)      this->GetElementDistanceToZeroLevelSet(&gl_distance,MaskOceanLevelsetEnum);
   if(this->amr->icefront_distance>0)           this->GetElementDistanceToZeroLevelSet(&if_distance,MaskIceLevelsetEnum);
   if(this->amr->thicknesserror_threshold>0)    this->ThicknessZZErrorEstimator(&thicknesserror);
   if(this->amr->deviatoricerror_threshold>0)   this->ZZErrorEstimator(&deviatoricerror);
 
   if(my_rank==0){
      this->amr->ExecuteRefinement(gl_distance,if_distance,deviatoricerror,thicknesserror,
                                    &newdatalist,&newxylist,&newelementslist);
      if(newdatalist[0]<=0 || newdatalist[1]<=0) _error_("Error in the ReMeshNeopz.");
   }
 
   /*Send new mesh to others CPU's*/
   if(my_rank) newdatalist=xNew<int>(2);
   ISSM_MPI_Bcast(newdatalist,2,ISSM_MPI_INT,0,IssmComm::GetComm());
   newnumberofvertices=newdatalist[0];
   newnumberofelements=newdatalist[1];
   if(my_rank){
      newxylist      =xNew<IssmDouble>(newnumberofvertices*2);
      newelementslist=xNew<int>(newnumberofelements*this->GetElementsWidth());
   }
   ISSM_MPI_Bcast(newxylist,newnumberofvertices*2,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
   ISSM_MPI_Bcast(newelementslist,newnumberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm());
 
   /*Reorganize the data*/
   newx=xNew<IssmDouble>(newnumberofvertices);
   newy=xNew<IssmDouble>(newnumberofvertices);
   newz=xNewZeroInit<IssmDouble>(newnumberofvertices);
   for(int i=0;i<newnumberofvertices;i++){
      newx[i] = newxylist[2*i];
      newy[i] = newxylist[2*i+1];
   }
 
   /*Assign the pointers*/
   (*pnewelementslist)  = newelementslist; //Matlab indexing
   (*pnewx)             = newx;
   (*pnewy)             = newy;
   (*pnewz)             = newz;
   *pnewnumberofvertices= newnumberofvertices;
   *pnewnumberofelements= newnumberofelements;
 
   /*Cleanup*/
   xDelete<int>(newdatalist);
   xDelete<IssmDouble>(newxylist);
   xDelete<IssmDouble>(deviatoricerror);
   xDelete<IssmDouble>(thicknesserror);
   xDelete<IssmDouble>(gl_distance);
   xDelete<IssmDouble>(if_distance);
}
/*}}}*/
void FemModel::InitializeAdaptiveRefinementNeopz(void){/*{{{*/
 
   /*Define variables*/
   int my_rank                            = IssmComm::GetRank();
   int numberofvertices                   = this->vertices->NumberOfVertices();
   int numberofelements                   = this->elements->NumberOfElements();
   IssmDouble* x                          = NULL;
   IssmDouble* y                          = NULL;
   int* elements                          = NULL;
   int amr_restart;
 
   /*Initialize field as NULL for now*/
   this->amr = NULL;
 
   /*Get vertices coordinates of the coarse mesh (father mesh)*/
   /*elements comes in Matlab indexing*/
   this->GetMesh(this->vertices,this->elements,&x,&y,&elements);
 
   /*Create initial mesh (coarse mesh) in neopz data structure*/
   /*Just CPU #0 should keep AMR object*/
   /*Initialize refinement pattern*/
   this->SetRefPatterns();
   this->amr = new AdaptiveMeshRefinement();
   this->amr->refinement_type=1;//1 is refpattern; 0 is uniform (faster)
   /*Get amr parameters*/
   this->parameters->FindParam(&this->amr->level_max,AmrLevelMaxEnum);
   this->parameters->FindParam(&this->amr->gradation,AmrGradationEnum);
   this->parameters->FindParam(&this->amr->lag,AmrLagEnum);
   this->parameters->FindParam(&this->amr->groundingline_distance,AmrGroundingLineDistanceEnum);
   this->parameters->FindParam(&this->amr->icefront_distance,AmrIceFrontDistanceEnum);
   this->parameters->FindParam(&this->amr->thicknesserror_threshold,AmrThicknessErrorThresholdEnum);
   this->parameters->FindParam(&this->amr->thicknesserror_groupthreshold,AmrThicknessErrorGroupThresholdEnum);
   this->parameters->FindParam(&this->amr->thicknesserror_maximum,AmrThicknessErrorMaximumEnum);
   this->parameters->FindParam(&this->amr->deviatoricerror_threshold,AmrDeviatoricErrorThresholdEnum);
   this->parameters->FindParam(&this->amr->deviatoricerror_groupthreshold,AmrDeviatoricErrorGroupThresholdEnum);
   this->parameters->FindParam(&this->amr->deviatoricerror_maximum,AmrDeviatoricErrorMaximumEnum);
 
   /*Initialize NeoPZ data structure*/
   this->amr->SetMesh(&elements,&x,&y,&numberofvertices,&numberofelements);
   if(my_rank==0){
      this->parameters->FindParam(&amr_restart,AmrRestartEnum);
      if(amr_restart){//experimental
         this->amr->ReadMesh();
      } else {//this is the default method
         this->amr->Initialize();
      }
   }
}
/*}}}*/
void FemModel::GetElementDistanceToZeroLevelSet(IssmDouble** pelementdistance,int levelset_type){/*{{{*/
 
   /*Here, "zero level set" means grounding line or ice front, depending on the level set type*/
   /*pverticedistance is the minimal vertice distance to the grounding line or ice front*/
   if(levelset_type!=MaskOceanLevelsetEnum && levelset_type!=MaskIceLevelsetEnum){
      _error_("level set type not implemented yet!");
   }
 
   /*Output*/
   IssmDouble* elementdistance;
 
   /*Intermediaries*/
   int numberofelements                   = this->elements->NumberOfElements();
   Vector<IssmDouble>* velementdistance   = new Vector<IssmDouble>(numberofelements);
   IssmDouble* levelset_points            = NULL;
   IssmDouble* xyz_list                   = NULL;
   IssmDouble mindistance,distance;
   IssmDouble xc,yc,x1,y1,x2,y2,x3,y3;
   int numberofpoints;
 
   /*Get points which level set is zero (center of elements with zero level set, levelset_points is serial)*/
   this->GetZeroLevelSetPoints(&levelset_points,numberofpoints,levelset_type);
 
   for(int i=0;i<this->elements->Size();i++){//parallel
      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
      int sid = element->Sid();
      element->GetVerticesCoordinates(&xyz_list);
      x1 = xyz_list[3*0+0];y1 = xyz_list[3*0+1];
      x2 = xyz_list[3*1+0];y2 = xyz_list[3*1+1];
      x3 = xyz_list[3*2+0];y3 = xyz_list[3*2+1];
      xc = (x1+x2+x3)/3.;
      yc = (y1+y2+y3)/3.;
      mindistance=INFINITY;
      /*Loop over each point (where level set is zero)*/
      for(int j=0;j<numberofpoints;j++){
         distance =sqrt((xc-levelset_points[2*j])*(xc-levelset_points[2*j])+(yc-levelset_points[2*j+1])*(yc-levelset_points[2*j+1]));
         mindistance=min(distance,mindistance);
      }
      velementdistance->SetValue(sid,mindistance,INS_VAL);
      xDelete<IssmDouble>(xyz_list);
   }
 
   /*Assemble*/
   velementdistance->Assemble();
 
   /*Assign the pointer*/
   (*pelementdistance)=velementdistance->ToMPISerial();
 
   /*Cleanup*/
   xDelete<IssmDouble>(levelset_points);
   xDelete<IssmDouble>(xyz_list);
   delete velementdistance;
}
/*}}}*/
void FemModel::SetRefPatterns(){/*{{{*/
 
   /*Initialize the global variable of refinement patterns*/
   gRefDBase.InitializeUniformRefPattern(ETriangle);
 
   /*Insert specifics patterns to ISSM core*/
   std::string filepath  = REFPATTERNDIR;
   std::string filename1 = filepath + "/2D_Triang_Rib_3.rpt";
   std::string filename2 = filepath + "/2D_Triang_Rib_4.rpt";
   std::string filename3 = filepath + "/2D_Triang_Rib_5.rpt";
   std::string filename4 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_4.rpt";
   std::string filename5 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_4_permuted.rpt";
   std::string filename6 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_5.rpt";
   std::string filename7 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_5_permuted.rpt";
   std::string filename8 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_4_5.rpt";
   std::string filename9 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_4_5_permuted.rpt";
 
   TPZAutoPointer<TPZRefPattern> refpat1 = new TPZRefPattern(filename1);
   TPZAutoPointer<TPZRefPattern> refpat2 = new TPZRefPattern(filename2);
   TPZAutoPointer<TPZRefPattern> refpat3 = new TPZRefPattern(filename3);
   TPZAutoPointer<TPZRefPattern> refpat4 = new TPZRefPattern(filename4);
   TPZAutoPointer<TPZRefPattern> refpat5 = new TPZRefPattern(filename5);
   TPZAutoPointer<TPZRefPattern> refpat6 = new TPZRefPattern(filename6);
   TPZAutoPointer<TPZRefPattern> refpat7 = new TPZRefPattern(filename7);
   TPZAutoPointer<TPZRefPattern> refpat8 = new TPZRefPattern(filename8);
   TPZAutoPointer<TPZRefPattern> refpat9 = new TPZRefPattern(filename9);
 
   if(!gRefDBase.FindRefPattern(refpat1)) gRefDBase.InsertRefPattern(refpat1);
   if(!gRefDBase.FindRefPattern(refpat2)) gRefDBase.InsertRefPattern(refpat2);
   if(!gRefDBase.FindRefPattern(refpat3)) gRefDBase.InsertRefPattern(refpat3);
   if(!gRefDBase.FindRefPattern(refpat4)) gRefDBase.InsertRefPattern(refpat4);
   if(!gRefDBase.FindRefPattern(refpat5)) gRefDBase.InsertRefPattern(refpat5);
   if(!gRefDBase.FindRefPattern(refpat6)) gRefDBase.InsertRefPattern(refpat6);
   if(!gRefDBase.FindRefPattern(refpat7)) gRefDBase.InsertRefPattern(refpat7);
   if(!gRefDBase.FindRefPattern(refpat8)) gRefDBase.InsertRefPattern(refpat8);
   if(!gRefDBase.FindRefPattern(refpat9)) gRefDBase.InsertRefPattern(refpat9);
}
/*}}}*/
#endif