Changeset 21802

Timestamp:

07/17/17 21:18:32 (8 years ago)

Author:

Mathieu Morlighem

Message:

CHG: working on bamg amr

Location:

issm/trunk-jpl/src

Files:

• c/Makefile.am (modified) (1 diff)
• c/bamg/BamgOpts.cpp (modified) (1 diff)
• c/bamg/Mesh.cpp (modified) (3 diffs)
• c/classes/AmrBamg.cpp (added)
• c/classes/AmrBamg.h (added)
• c/classes/FemModel.cpp (modified) (20 diffs)
• c/classes/FemModel.h (modified) (5 diffs)
• c/classes/IoModel.cpp (modified) (1 diff)
• c/cores/transient_core.cpp (modified) (2 diffs)
• c/modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp (modified) (4 diffs)
• c/modules/ModelProcessorx/CreateParameters.cpp (modified) (5 diffs)
• c/shared/Enum/EnumDefinitions.h (modified) (1 diff)
• c/shared/Enum/EnumToStringx.cpp (modified) (1 diff)
• c/shared/Enum/StringToEnumx.cpp (modified) (5 diffs)
• c/shared/io/Marshalling/IoCodeConversions.cpp (modified) (1 diff)
• c/shared/io/Marshalling/IoCodeConversions.h (modified) (1 diff)
• m/classes/amr.m (modified) (1 diff)
• m/classes/amrbamg.m (added)
• m/classes/trans.js (modified) (1 diff)
• m/classes/transient.m (modified) (1 diff)
• m/classes/transient.py (modified) (1 diff)

issm/trunk-jpl/src/c/Makefile.am

@@ -47 +47 @@
                                          ./bamg/Mesh.cpp\
                                          ./shared/Bamg/BigPrimeNumber.cpp\
+                                         ./classes/AmrBamg.cpp\
                                          ./modules/Bamgx/Bamgx.cpp\
                                          ./modules/BamgConvertMeshx/BamgConvertMeshx.cpp\

issm/trunk-jpl/src/c/bamg/BamgOpts.cpp

@@ -5 +5 @@
 BamgOpts::BamgOpts(){/*{{{*/
 
-        this->anisomax=0;
-        this->cutoff=0;
-        this->coeff=0;
-        this->errg=0;
-        this->gradation=0;
-        this->Hessiantype=0;
-        this->MaxCornerAngle=0;
-        this->maxnbv=0;
-        this->maxsubdiv=0;
-        this->Metrictype=0;
-        this->nbjacobi=0;
-        this->nbsmooth=0;
-        this->omega=0;
-        this->power=0;
-        this->random=0;
-        this->verbose=0;
+        this->anisomax          = 0;
+        this->cutoff            = 0;
+        this->coeff             = 0;
+        this->errg              = 0;
+        this->gradation         = 0;
+        this->Hessiantype       = 0;
+        this->MaxCornerAngle    = 0;
+        this->maxnbv            = 0;
+        this->maxsubdiv         = 0;
+        this->Metrictype        = 0;
+        this->nbjacobi          = 0;
+        this->nbsmooth          = 0;
+        this->omega             = 0;
+        this->power             = 0;
+        this->random            = 0;
+        this->verbose           = 0;
 
-        this->Crack=0;
-        this->geometricalmetric=0;
-        this->KeepVertices=0;
-        this->splitcorners=0;
+        this->Crack             = 0;
+        this->geometricalmetric = 0;
+        this->KeepVertices      = 0;
+        this->splitcorners      = 0;
 
-        this->hmin=0;
-        this->hmax=0;
+        this->hmin              = 0;
+        this->hmax              = 0;
         this->hminVertices=NULL; this->hminVerticesSize[0]=this->hminVerticesSize[1]=0;
         this->hmaxVertices=NULL; this->hmaxVerticesSize[0]=this->hmaxVerticesSize[1]=0;

issm/trunk-jpl/src/c/bamg/Mesh.cpp

@@ -698 +698 @@
                         for (i=0;i<nbsubdomains;i++){
                                 bamgmesh->SubDomains[i*4+0]=3;
-                                bamgmesh->SubDomains[i*4+1]=reft[GetId(subdomains[i].head)];
+                                bamgmesh->SubDomains[i*4+1]=reft[GetId(subdomains[i].head)]+1;//MATLAB indexing
                                 bamgmesh->SubDomains[i*4+2]=1;
                                 bamgmesh->SubDomains[i*4+3]=subdomains[i].ReferenceNumber;
@@ -3076 +3076 @@
                                 IssmDouble area_3 =((bv.r.x -(*tcvj)(1)->r.x)*((*tcvj)(0)->r.y-(*tcvj)(1)->r.y)
                                                 - (bv.r.y -(*tcvj)(1)->r.y)*((*tcvj)(0)->r.x-(*tcvj)(1)->r.x));
-                                if(area_1<0 || area_2<0 || area_3<0){
-                                        pointsoutside = true;
-                                        continue;
-                                }
+                        //      if(area_1<0 || area_2<0 || area_3<0){
+                        //              pointsoutside = true;
+                        //              continue;
+                        //      }
                                 if(!bv.GeomEdgeHook){
                                         vertices[nbv].r              = bv.r;
@@ -3086 +3086 @@
                                 }
                         }
-                        if(pointsoutside) _printf_("WARNING: One or more points of the initial mesh fall outside of the geometric boundary\n");
+                        //if(pointsoutside) _printf_("WARNING: One or more points of the initial mesh fall outside of the geometric boundary\n");
                         Bh.CreateSingleVertexToTriangleConnectivity();     
                         InsertNewPoints(nbvold,NbTSwap,bamgopts->random);

issm/trunk-jpl/src/c/classes/FemModel.cpp

@@ -45 +45 @@
 
         /*configuration: */
-        int  solution_type;
+        int  solution_type,amrtype,amr_frequency;
         int  ierr;
 
@@ -80 +80 @@
         this->parameters->AddObject(new GenericParam<ISSM_MPI_Comm>(incomm,FemModelCommEnum));
 
-        #ifdef _HAVE_NEOPZ_
-        this->InitializeAdaptiveRefinement();
-        #endif
+   /*AMR stuff*/
+        this->parameters->FindParam(&amr_frequency,TransientAmrFrequencyEnum);
+        if(amr_frequency){
+                this->parameters->FindParam(&amrtype,AmrTypeEnum);
+                switch(amrtype){
+
+                        #ifdef _HAVE_NEOPZ_
+                        case AmrNeopzEnum: this->InitializeAdaptiveRefinementNeopz(); break;
+                        #endif
+
+                        #ifdef _HAVE_BAMG_
+                        case AmrBamgEnum: this->InitializeAdaptiveRefinementBamg(); break;
+                        #endif
+
+                        default: _error_("not implemented yet");
+                }
+        }
 
         /*Free resources */
@@ -140 +154 @@
         #ifdef _HAVE_NEOPZ_
         if(amr)delete amr;
+        #endif
+
+        #ifdef _HAVE_BAMG_
+        if(amrbamg)delete amrbamg;
         #endif
 
@@ -1502 +1520 @@
 }
 /*}}}*/
+void FemModel::ReMesh(void){/*{{{*/
+
+        /*Intermediaries*/
+        IssmDouble *newx                        = NULL;
+        IssmDouble *newy                        = NULL;
+        IssmDouble *newz                        = NULL;
+        int *newelementslist            = NULL;
+        int newnumberofvertices = -1;
+        int newnumberofelements = -1;
+        bool* my_elements                       = NULL; 
+        int* my_vertices                        = NULL;
+        int elementswidth       = this->GetElementsWidth();//just tria elements in this version
+        int amrtype;
+
+        /*Branch to specific amr depending on requested method*/        
+        parameters->FindParam(&amrtype,AmrTypeEnum);
+        switch(amrtype){
+                #ifdef _HAVE_NEOPZ_
+                case AmrNeopzEnum: this->ReMeshNeopz(&newnumberofvertices,&newnumberofelements,&newx,&newy,&newz,&newelementslist); break;
+                #endif
+
+                #ifdef _HAVE_BAMG_
+                case AmrBamgEnum: this->ReMeshBamg(&newnumberofvertices,&newnumberofelements,&newx,&newy,&newz,&newelementslist); break;
+                #endif
+
+                default: _error_("not implemented yet");
+        }
+
+        /*Partitioning the new mesh. Maybe ElementsAndVerticesPartitioning.cpp could be modified to set this without iomodel.*/
+        this->ElementsAndVerticesPartitioning(newnumberofvertices,newnumberofelements,elementswidth,newelementslist,&my_elements,&my_vertices);
+
+        if(this->loads->Size()!=0) _error_("not supported yet");
+
+        /*Create vertices*/
+        Vertices* new_vertices=new Vertices();
+        this->CreateVertices(newnumberofvertices,newnumberofelements,elementswidth,newelementslist,my_vertices,newx,newy,newz,new_vertices);
+
+        /*Creating elements*/
+        /*Just Tria in this version*/
+        Elements* new_elements=new Elements();
+        this->CreateElements(newnumberofelements,elementswidth,newelementslist,my_elements,new_elements);
+
+        /*Creating materials*/
+        Materials* new_materials=new Materials();
+        this->CreateMaterials(newnumberofelements,my_elements,new_materials);
+
+        /*Creating nodes and constraints*/
+        /*Just SSA (2D) and P1 in this version*/
+        Nodes* new_nodes=new Nodes();
+        Constraints* new_constraints=new Constraints();
+
+        int nodecounter=0;
+        int constraintcounter=0;
+        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];
+
+                /*As the domain is 2D, it is not necessary to create nodes for this analysis*/
+                /*itapopo must verify if domain is not 3D. Only 2D in this version!*/
+                if(analysis_enum==StressbalanceVerticalAnalysisEnum) continue;      
+
+                this->CreateNodes(newnumberofvertices,my_vertices,nodecounter,analysis_enum,new_nodes);
+                if(analysis_enum==StressbalanceAnalysisEnum) this->CreateConstraints(newnumberofvertices,newnumberofelements,nodecounter,constraintcounter,newx,newy,my_vertices,new_constraints);
+                this->UpdateElements(newnumberofelements,newelementslist,my_elements,nodecounter,i,new_elements);
+
+                if(new_nodes->Size()) nodecounter=new_nodes->MaximumId();
+                constraintcounter = new_constraints->NumberOfConstraints();
+                /*Make sure nodecounter is at least 0 (if no node exists, maxid will be -1*/
+                _assert_(nodecounter>=0);
+        }
+
+        new_elements->Presort();
+        new_nodes->Presort();
+        new_vertices->Presort();
+        this->loads->Presort();
+        new_materials->Presort();
+        new_constraints->Presort();
+
+        /*reset hooks for elements, loads and nodes: */
+        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,new_vertices,new_materials,this->parameters);
+                this->loads->SetCurrentConfiguration(new_elements,this->loads,new_nodes,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,new_vertices,new_materials,this->parameters);
+                if(i==0){ 
+                        VerticesDofx(new_vertices,this->parameters); //only call once, we only have one set of vertices
+                }
+                SpcNodesx(new_nodes,new_constraints,this->parameters,analysis_type);
+                NodesDofx(new_nodes,this->parameters,analysis_type);
+        }
+
+        /*Finally: interpolate all inputs and insert them into the new elements.*/
+        this->InterpolateInputs(new_vertices,new_elements);
+
+        /*Delete old structure and set new pointers*/
+        delete this->vertices;          this->vertices          = new_vertices;
+        delete this->elements;          this->elements          = new_elements;
+        delete this->nodes;                     this->nodes                     = new_nodes;
+        delete this->constraints;       this->constraints       = new_constraints;
+        delete this->materials;         this->materials = new_materials;
+
+        GetMaskOfIceVerticesLSMx(this);
+
+        /*Insert MISMIP+ bed topography*/
+        if(false) this->BedrockFromMismipPlus();
+
+        /*Adjust base, thickness and mask grounded ice leve set*/
+        if(true) this->AdjustBaseThicknessAndMask();
+
+        /*Reset current configuration: */
+        analysis_type=this->analysis_type_list[this->analysis_counter];
+        SetCurrentConfiguration(analysis_type);
+
+        /*Cleanup*/
+        xDelete<IssmDouble>(newx);
+        xDelete<IssmDouble>(newy);
+        xDelete<IssmDouble>(newz);
+        xDelete<int>(newelementslist);
+        xDelete<int>(my_vertices);
+        xDelete<bool>(my_elements);
+}
+/*}}}*/
 void FemModel::RequestedDependentsx(void){/*{{{*/
 
@@ -2234 +2396 @@
 }
 /*}}}*/
-#ifdef  _HAVE_DAKOTA_
-void FemModel::DakotaResponsesx(double* d_responses,char** responses_descriptors,int numresponsedescriptors,int d_numresponses){/*{{{*/
-
-        int        i,j;
-        int        my_rank;
-
-        /*intermediary: */
-        char   root[50];
-        int    index;
-        int    npart;
-        double femmodel_response;
-        int    flag;
-        double *vertex_response   = NULL;
-        double *qmu_response      = NULL;
-        double *responses_pointer = NULL;
-
-        /*retrieve npart: */
-        parameters->FindParam(&npart,QmuNumberofpartitionsEnum);
-        my_rank=IssmComm::GetRank();
-
-        /*save the d_responses pointer: */
-        responses_pointer=d_responses;
-
-        //watch out, we have more d_numresponses than numresponsedescriptors, because the responses have been expanded if they were scaled. 
-        //because we don't know the d_responses descriptors (the scaled ones) we can't key off them, so we will key off the responses_descriptors: */
-
-        for(i=0;i<numresponsedescriptors;i++){
-
-                flag=DescriptorIndex(root,&index,responses_descriptors[i]);
-
-                if(flag==ScaledEnum){
-
-                        /*this response was scaled. pick up the response from the inputs: */
-                        GetVectorFromInputsx(&vertex_response,this, StringToEnumx(root),VertexPIdEnum);
-
-                        /*Now, average it onto the partition nodes: */
-                        AverageOntoPartitionx(&qmu_response,elements,nodes,vertices,loads,materials,parameters,vertex_response);
-
-                        /*Copy onto our dakota responses: */
-                        if(my_rank==0){
-                                /*plug response: */
-                                for(j=0;j<npart;j++)responses_pointer[j]=qmu_response[j];
-
-                                /*increment response_pointer :*/
-                                responses_pointer+=npart;
-                        }
-
-                        /*Free ressources:*/
-                        xDelete<double>(vertex_response);
-                        xDelete<double>(qmu_response);
-
-                }
-                else if (flag==IndexedEnum){
-
-                        /*indexed response: plug index into parameters and call response module: */
-                        parameters->SetParam(index,IndexEnum);
-
-                        this->Responsex(&femmodel_response,root);
-
-                        if(my_rank==0){
-                                /*plug response: */
-                                responses_pointer[0]=femmodel_response;
-
-                                /*increment response_pointer :*/
-                                responses_pointer++;
-                        }
-                }
-                else if (flag==NodalEnum){
-                        _error_("nodal response functions not supported yet!");
-
-                        /*increment response_pointer :*/
-                        responses_pointer++;
-                }
-                else if (flag==RegularEnum){
-
-                        /*perfectly normal response function: */
-                        this->Responsex(&femmodel_response,root);
-
-                        if(my_rank==0){
-                                /*plug response: */
-                                responses_pointer[0]=femmodel_response;
-
-                                /*increment response_pointer :*/
-                                responses_pointer++;
-                        }
-                }
-                else _error_("flag type " << flag << " not supported yet for response analysis");
-        }
-
-        /*Synthesize echo: {{{*/
-        if(my_rank==0){
-                _printf_("   responses: " << d_numresponses << ": ");
-                for(i=0;i<d_numresponses-1;i++)_printf_(d_responses[i] << "|");
-                _printf_(d_responses[d_numresponses-1]);
-                _printf_("\n");
-        }
-        /*}}}*/
-
-}
-/*}}}*/
-#endif
-#ifdef _HAVE_GIAIVINS_
-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, 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,xx,yy);
-                }
-                if(i%100==0){
-                        pUp->Assemble();
-                        pNorth->Assemble();
-                        pEast->Assemble();
-                }
-        }
-        
-        /*One last time: */
-        pUp->Assemble();
-        pNorth->Assemble();
-        pEast->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){/*{{{*/
-
-        IssmDouble  eartharea=0;
-        IssmDouble  eartharea_cpu=0;
-
-        int         ns,nsmax;
-        
-        /*Go through elements, and add contribution from each element to the deflection vector wg:*/
-        ns = elements->Size();
-        
-        /*First, figure out the surface area of Earth: */ 
-        for(int i=0;i<ns;i++){
-                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
-                eartharea_cpu += element->GetAreaSpherical();
-        }
-        ISSM_MPI_Reduce (&eartharea_cpu,&eartharea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
-        ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-
-        /*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,eartharea);
-                }
-                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>* pSgi, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius) { /*{{{*/
-
-        /*serialized vectors:*/
-        IssmDouble  eustatic       = 0.;
-        IssmDouble  eustatic_cpu   = 0.;
-        IssmDouble  eustatic_cpu_e = 0.;
-        IssmDouble  oceanarea      = 0.;
-        IssmDouble  oceanarea_cpu  = 0.;
-        IssmDouble  eartharea      = 0.;
-        IssmDouble  eartharea_cpu  = 0.;
-        int         ns,nsmax;
-        
-        /*Go through elements, and add contribution from each element to the deflection vector wg:*/
-        ns = elements->Size();
-
-        /*First, figure out the area of the ocean, which is needed to compute the eustatic component: */
-        for(int i=0;i<ns;i++){
-                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
-                oceanarea_cpu += element->OceanArea();
-                eartharea_cpu += element->GetAreaSpherical();
-        }
-        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.);
-
-        ISSM_MPI_Reduce (&eartharea_cpu,&eartharea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
-        ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-
-        /*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 sea level rise core: */
-        for(int i=0;i<nsmax;i++){
-                if(i<ns){
-                
-                        if(VerboseConvergence())if(i%100==0)_printf0_("\r" << "      convolution progress: " << (double)i/(double)ns*100 << "%  ");
-                
-                        Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
-                        element->SealevelriseEustatic(pSgi,&eustatic_cpu_e,latitude,longitude,radius,oceanarea,eartharea);
-                        eustatic_cpu+=eustatic_cpu_e;
-                }
-                if(i%100==0)pSgi->Assemble();
-        }
-        if(VerboseConvergence())_printf0_("\n");
-                
-        /*One last time: */
-        pSgi->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));
-
-        /*Assign output pointers:*/
-        *peustatic=eustatic;
-
-}
-/*}}}*/
-void FemModel::SealevelriseNonEustatic(Vector<IssmDouble>* pSgo, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, bool verboseconvolution){/*{{{*/
-
-        /*serialized vectors:*/
-        IssmDouble* Sg_old=NULL;
-        
-        IssmDouble  eartharea=0;
-        IssmDouble  eartharea_cpu=0;
-
-        int         ns,nsmax;
-        
-        /*Serialize vectors from previous iteration:*/
-        Sg_old=pSg_old->ToMPISerial();
-
-        /*Go through elements, and add contribution from each element to the deflection vector wg:*/
-        ns = elements->Size();
-
-        /*First, figure out the area of the ocean, which is needed to compute the eustatic component: */
-        for(int i=0;i<ns;i++){
-                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
-                eartharea_cpu += element->GetAreaSpherical();
-        }
-        
-        ISSM_MPI_Reduce (&eartharea_cpu,&eartharea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
-        ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-
-        /*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 sea level rise core: */
-        for(int i=0;i<nsmax;i++){
-                if(i<ns){
-                        if(verboseconvolution)if(VerboseConvergence())if(i%100==0)_printf_("\r" << "      convolution progress: " << (double)i/(double)ns*100 << "%   ");
-                        Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
-                        element->SealevelriseNonEustatic(pSgo,Sg_old,latitude,longitude,radius,eartharea);
-                }
-                if(i%100==0)pSgo->Assemble();
-        }
-        if(verboseconvolution)if(VerboseConvergence())_printf_("\n");
-        
-        /*Free ressources:*/
-        xDelete<IssmDouble>(Sg_old);
-}
-/*}}}*/
-void FemModel::SealevelriseRotationalFeedback(Vector<IssmDouble>* pSgo_rot, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius){/*{{{*/
-
-        /*serialized vectors:*/
-        IssmDouble* Sg_old=NULL;
-        IssmDouble  eartharea=0;
-        IssmDouble  eartharea_cpu=0;
-        IssmDouble      tide_love_h, tide_love_k, fluid_love, moi_e, moi_p, omega, g;
-        IssmDouble      load_love_k2 = -0.30922675; //degree 2 load Love number 
-        IssmDouble      m1, m2, m3; 
-        IssmDouble      lati, longi, radi, value; 
-
-        /*Serialize vectors from previous iteration:*/
-        Sg_old=pSg_old->ToMPISerial();
-
-        /*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));
-                eartharea_cpu += element->GetAreaSpherical();
-        }
-        ISSM_MPI_Reduce (&eartharea_cpu,&eartharea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
-        ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-
-        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],Sg_old,eartharea);
-                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(&tide_love_h,SealevelriseTidalLoveHEnum);
-        parameters->FindParam(&tide_love_k,SealevelriseTidalLoveKEnum);
-        parameters->FindParam(&fluid_love,SealevelriseFluidLoveEnum);
-        parameters->FindParam(&moi_e,SealevelriseEquatorialMoiEnum);
-        parameters->FindParam(&moi_p,SealevelrisePolarMoiEnum);
-        parameters->FindParam(&omega,SealevelriseAngularVelocityEnum);
-
-        /*compute perturbation terms for angular velocity vector: */
-        m1 = 1/(1-tide_love_k/fluid_love) * (1+load_love_k2)/(moi_p-moi_e) * moi_list[0]; 
-        m2 = 1/(1-tide_love_k/fluid_love) * (1+load_love_k2)/(moi_p-moi_e) * moi_list[1]; 
-        m3 = -(1+load_love_k2)/moi_p * moi_list[2];     // term associated with fluid number (3-order-of-magnitude smaller) is negelected  
-
-        /* 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-tide_love_h)/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))); 
-        
-                pSgo_rot->SetValue(sid,value,INS_VAL); //INS_VAL ensures that you don't add several times
-        }
-
-        /*Assemble mesh velocity*/
-        pSgo_rot->Assemble();
-        
-        /*Free ressources:*/
-        xDelete<IssmDouble>(Sg_old);
-        
-}
-/*}}}*/
-void FemModel::SealevelriseGeodetic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pSg, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){/*{{{*/
-
-        /*serialized vectors:*/
-        IssmDouble* Sg=NULL;
-        
-        IssmDouble  eartharea=0;
-        IssmDouble  eartharea_cpu=0;
-
-        int         ns,nsmax;
-        
-        /*Serialize vectors from previous iteration:*/
-        Sg=pSg->ToMPISerial();
-
-        /*Go through elements, and add contribution from each element to the deflection vector wg:*/
-        ns = elements->Size();
-        
-        /*First, figure out the area of the ocean, which is needed to compute the eustatic component: */
-        for(int i=0;i<ns;i++){
-                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
-                eartharea_cpu += element->GetAreaSpherical();
-        }
-        ISSM_MPI_Reduce (&eartharea_cpu,&eartharea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
-        ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-
-        /*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 sea level rise core: */
-        for(int i=0;i<nsmax;i++){
-                if(i<ns){
-                        Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
-                        element->SealevelriseGeodetic(pUp,pNorth,pEast,Sg,latitude,longitude,radius,xx,yy,zz,eartharea);
-                }
-                if(i%100==0){
-                        pUp->Assemble();
-                        pNorth->Assemble();
-                        pEast->Assemble();
-                }
-        }
-        
-        /*One last time: */
-        pUp->Assemble();
-        pNorth->Assemble();
-        pEast->Assemble();
-
-        /*Free ressources:*/
-        xDelete<IssmDouble>(Sg);
-        xDelete<IssmDouble>(latitude);
-        xDelete<IssmDouble>(longitude);
-        xDelete<IssmDouble>(radius);
-        xDelete<IssmDouble>(xx);
-        xDelete<IssmDouble>(yy);
-        xDelete<IssmDouble>(zz);
-}
-/*}}}*/
-IssmDouble FemModel::SealevelriseOceanAverage(Vector<IssmDouble>* Sg) { /*{{{*/
-
-        IssmDouble* Sg_serial=NULL;
-        IssmDouble  oceanvalue,oceanvalue_cpu;
-        IssmDouble  oceanarea,oceanarea_cpu;
-
-        /*Serialize vectors from previous iteration:*/
-        Sg_serial=Sg->ToMPISerial();
-
-        /*Initialize:*/
-        oceanvalue_cpu=0;
-        oceanarea_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));
-                oceanarea_cpu += element->OceanArea();
-                oceanvalue_cpu += element->OceanAverage(Sg_serial);
-        }
-        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());
-        
-        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>(Sg_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*/
-        mask=new Vector<IssmDouble>(this->nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum));
-        recurence=new Vector<IssmDouble>(this->nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum));
-        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,eplzigzag_counter,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(HydrologyDCEfficientAnalysisEnum));
-        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;
-        for (int i=0;i<nodes->Size();i++){
-                Node* node=xDynamicCast<Node*>(nodes->GetObjectByOffset(i));
-                if(node->InAnalysis(HydrologyDCEfficientAnalysisEnum)){
-                        if(serial_active[node->Sid()]==1.){
-                                node->Activate();
-                                if(!node->IsClone()) counter++;
-                        }
-                        else{
-                                node->Deactivate();
-                        }
-                }
-        }
-        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()) _printf0_("   Number of active nodes in EPL layer: "<< counter <<"\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*/
-        active=new Vector<IssmDouble>(nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum));
-        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;
-        for (int i=0;i<nodes->Size();i++){
-                Node* node=xDynamicCast<Node*>(nodes->GetObjectByOffset(i));
-                if(node->InAnalysis(L2ProjectionEPLAnalysisEnum)){
-                        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");
-}
-/*}}}*/
-#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->Tag(START);
-
-        /*From command line arguments, retrieve different filenames needed to create the FemModel: */
-        solution_type=StringToEnumx(solution);
-        
-        /*Create femmodel from input files: */
-        profiler->Tag(STARTINIT);
-        this->InitFromBuffers((char*)buffer,buffersize,toolkits, solution_type,trace,NULL);
-        profiler->Tag(FINISHINIT);
-        
-        /*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->Tag(FINISH);  //final tagging
-        _printf0_("\n");
-        _printf0_("   "<<setw(40)<<left<<"FemModel initialization elapsed time:"<<profiler->DeltaTime(STARTINIT,FINISHINIT) << "\n");
-        _printf0_("   "<<setw(40)<<left<<"Core solution elapsed time:"<<profiler->DeltaTime(STARTCORE,FINISHCORE) << "\n");
-        _printf0_("\n");
-        _printf0_("   Total elapsed time: "
-                                <<profiler->DeltaTimeModHour(START,FINISH)<<" hrs "
-                                <<profiler->DeltaTimeModMin(START,FINISH)<<" min "
-                                <<profiler->DeltaTimeModSec(START,FINISH)<<" sec"
-                                );
-        _printf0_("\n");
-        
-        /*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
-
-#ifdef _HAVE_NEOPZ_
-void FemModel::InitializeAdaptiveRefinement(void){/*{{{*/
-        
-        /*Define variables*/
-        int my_rank                                             = IssmComm::GetRank();
-        this->amr                                               = NULL;//initialize amr as NULL
-        int numberofvertices                    = this->vertices->NumberOfVertices();
-        int numberofelements                    = this->elements->NumberOfElements();
-        int numberofsegments                    = 0; //used on matlab
-        IssmDouble* x                                   = NULL;
-        IssmDouble* y                                   = NULL;
-        IssmDouble* z                                   = NULL;
-        int* elements                                   = NULL;
-        int elementswidth                               = this->GetElementsWidth(); //just tria elements in this version. Itapopo:
-        int levelmax                                    = 0;
-        IssmDouble regionlevel1         = 0.;
-        IssmDouble regionlevelmax       = 0.;
-
-        /*Get vertices coordinates of the coarse mesh (father mesh)*/
-        /*elements comes in Matlab indexing*/
-        this->GetMesh(this->vertices,this->elements,&x,&y,&z,&elements);
-        
-        /*Get amr parameters*/
-        this->parameters->FindParam(&levelmax,AmrLevelMaxEnum);
-        this->parameters->FindParam(&regionlevel1,AmrRegionLevel1Enum);
-        this->parameters->FindParam(&regionlevelmax,AmrRegionLevelMaxEnum);
-
-        /*Create initial mesh (coarse mesh) in neopz data structure*/ 
-        /*Just CPU #0 should keep AMR object*/
-   this->SetRefPatterns();
-        if(my_rank==0){ 
-           bool ismisomip       = false;
-                if(ismisomip){//itapopo
-                        TPZFileStream fstr;
-                        std::stringstream ss;
-                        int frictionlaw;
-                        
-                        this->parameters->FindParam(&frictionlaw,FrictionLawEnum);
-                
-                        ss      << levelmax;
-                        if(frictionlaw==1){
-                                ss << "_viscous/amr.txt";
-                        }else if(frictionlaw==7){
-                                ss << "_tsai/amr.txt";
-                        }else{
-                                _error_("friction law not supported here.");
-                        }
-                        
-                        std::string AMRfile  = "/home/santos/Misomip2/L" + ss.str();
-                        fstr.OpenRead(AMRfile.c_str());
-                        
-                        TPZSaveable *sv         = TPZSaveable::Restore(fstr,0);
-                        this->amr                               = dynamic_cast<AdaptiveMeshRefinement*>(sv);
-                }
-                else{
-                        this->amr = new AdaptiveMeshRefinement();
-                        //this->amr->SetLevelMax(levelmax); //Set max level of refinement
-                        //this->amr->SetRegions(regionlevel1,regionlevelmax);
-                        this->amr->CreateInitialMesh(numberofvertices,numberofelements,numberofsegments,elementswidth,x,y,z,elements,NULL);
-                }
-                this->amr->SetLevelMax(levelmax); //Set max level of refinement
-                this->amr->SetRegions(regionlevel1,regionlevelmax);
-        }
-
-        /*Free the vectors*/
-        xDelete<IssmDouble>(x);
-        xDelete<IssmDouble>(y);
-        xDelete<IssmDouble>(z);
-        xDelete<int>(elements);
-
-}
-/*}}}*/
-void FemModel::SetRefPatterns(){/*{{{*/
-
-   /*Initialize the global variable of refinement patterns*/
-   gRefDBase.InitializeUniformRefPattern(EOned);
-   gRefDBase.InitializeUniformRefPattern(ETriangle);
-
-    //gRefDBase.InitializeRefPatterns();
-   /*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_OnlyTriang_Side_3_4.rpt";
-   std::string filename4 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_4_permuted.rpt";
-   std::string filename5 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_5.rpt";
-   std::string filename6 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_5_permuted.rpt";
-   std::string filename7 = filepath + "/2D_Triang_Rib_5.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);
-
-   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);
-}
-/*}}}*/
-void FemModel::ReMesh(void){/*{{{*/
-
-        /*Variables*/
-        IssmDouble *newx                        = NULL;
-        IssmDouble *newy                        = NULL;
-        IssmDouble *newz                        = NULL;
-        int *newelementslist            = NULL;
-        int newnumberofvertices = -1;
-        int newnumberofelements = -1;
-        bool* my_elements                       = NULL; 
-        int* my_vertices                        = NULL;
-        int elementswidth                       = this->GetElementsWidth();//just tria elements in this version
-
-        /*Execute refinement and get the new mesh*/
-        /*newelementslist come in Matlab indexing*/
-        this->ExecuteRefinement(newnumberofvertices,newnumberofelements,&newx,&newy,&newz,&newelementslist);
-
-        /*Partitioning the new mesh. Maybe ElementsAndVerticesPartitioning.cpp could be modified to set this without iomodel.*/
-        this->ElementsAndVerticesPartitioning(newnumberofvertices,newnumberofelements,elementswidth,newelementslist,&my_elements,&my_vertices);
-
-        if(this->loads->Size()!=0) _error_("not supported yet");
-
-        /*Create vertices*/
-        Vertices* new_vertices=new Vertices();
-        this->CreateVertices(newnumberofvertices,newnumberofelements,elementswidth,newelementslist,my_vertices,newx,newy,newz,new_vertices);
- 
-        /*Creating elements*/
-        /*Just Tria in this version*/
-        Elements* new_elements=new Elements();
-        this->CreateElements(newnumberofelements,elementswidth,newelementslist,my_elements,new_elements);
-
-        /*Creating materials*/
-        Materials* new_materials=new Materials();
-        this->CreateMaterials(newnumberofelements,my_elements,new_materials);
-        
-        /*Creating nodes and constraints*/
-        /*Just SSA (2D) and P1 in this version*/
-        Nodes* new_nodes=new Nodes();
-        Constraints* new_constraints=new Constraints();
-
-        int nodecounter=0;
-        int constraintcounter=0;
-        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];
-                
-                /*As the domain is 2D, it is not necessary to create nodes for this analysis*/
-                /*itapopo must verify if domain is not 3D. Only 2D in this version!*/
-                if(analysis_enum==StressbalanceVerticalAnalysisEnum) continue;      
-                
-                this->CreateNodes(newnumberofvertices,my_vertices,nodecounter,analysis_enum,new_nodes);
-                if(analysis_enum==StressbalanceAnalysisEnum) this->CreateConstraints(newnumberofvertices,newnumberofelements,nodecounter,constraintcounter,newx,newy,my_vertices,new_constraints);
-                this->UpdateElements(newnumberofelements,newelementslist,my_elements,nodecounter,i,new_elements);
-
-                if(new_nodes->Size()) nodecounter=new_nodes->MaximumId();
-                constraintcounter = new_constraints->NumberOfConstraints();
-                /*Make sure nodecounter is at least 0 (if no node exists, maxid will be -1*/
-                _assert_(nodecounter>=0);
-        }
-
-        new_elements->Presort();
-        new_nodes->Presort();
-        new_vertices->Presort();
-        this->loads->Presort();
-        new_materials->Presort();
-        new_constraints->Presort();
-
-        /*reset hooks for elements, loads and nodes: */
-        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,new_vertices,new_materials,this->parameters);
-                this->loads->SetCurrentConfiguration(new_elements,this->loads,new_nodes,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,new_vertices,new_materials,this->parameters);
-                if(i==0){ 
-                        VerticesDofx(new_vertices,this->parameters); //only call once, we only have one set of vertices
-                }
-                SpcNodesx(new_nodes,new_constraints,this->parameters,analysis_type);
-                NodesDofx(new_nodes,this->parameters,analysis_type);
-        }
-
-        /*Finally: interpolate all inputs and insert them into the new elements.*/
-        this->InterpolateInputs(new_vertices,new_elements);
-
-        /*Delete old structure and set new pointers*/
-        delete this->vertices;          this->vertices          = new_vertices;
-        delete this->elements;          this->elements          = new_elements;
-        delete this->nodes;                     this->nodes                     = new_nodes;
-        delete this->constraints;       this->constraints       = new_constraints;
-        delete this->materials;         this->materials = new_materials;
-        
-        GetMaskOfIceVerticesLSMx(this);
-
-        /*Insert MISMIP+ bed topography*/
-        if(false) this->BedrockFromMismipPlus();
-        
-        /*Adjust base, thickness and mask grounded ice leve set*/
-        if(true) this->AdjustBaseThicknessAndMask();
-
-        /*Reset current configuration: */
-        analysis_type=this->analysis_type_list[this->analysis_counter];
-        SetCurrentConfiguration(analysis_type);
-
-        /*Cleanup*/
-        xDelete<IssmDouble>(newx);
-        xDelete<IssmDouble>(newy);
-        xDelete<IssmDouble>(newz);
-        xDelete<int>(newelementslist);
-        xDelete<int>(my_vertices);
-        xDelete<bool>(my_elements);
-
-        return;
-
-}
-/*}}}*/
+
+/*AMR*/
 void FemModel::BedrockFromMismipPlus(void){/*{{{*/
 
@@ -3245 +2466 @@
                         } 
 
-                        if(abs(sl[i])>0) _error_("Sea level value not supported!");
+                        if(abs(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];
@@ -3255 +2476 @@
                 element->AddInput(ThicknessEnum,&h[0],P1Enum);
                 element->AddInput(BaseEnum,&b[0],P1Enum);
-
         }
         
@@ -3265 +2485 @@
    xDelete<IssmDouble>(r);
    xDelete<IssmDouble>(sl);
-
-        return;
-}
-/*}}}*/
-void FemModel::WriteMeshInResults(void){/*{{{*/
-
-        int step                                        = -1;
-        int numberofelements = -1;
-        int numberofvertices = -1;
-        IssmDouble time         = -1;
-        IssmDouble* x                   = NULL;
-        IssmDouble* y                   = NULL;
-        IssmDouble* z                   = NULL;
-        int* elementslist               = NULL;
-
-        if(!this->elements || !this->vertices || !this->results || !this->parameters) return;
-         
-        parameters->FindParam(&step,StepEnum);
-        parameters->FindParam(&time,TimeEnum);
-        numberofelements=this->elements->NumberOfElements();
-        numberofvertices=this->vertices->NumberOfVertices();
-
-        /*Get mesh. Elementslist comes in Matlab indexing*/
-        this->GetMesh(this->vertices,this->elements,&x,&y,&z,&elementslist);
-
-        /*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));
-        
-        //itapopo
-        if(IssmComm::GetRank()==0){
-                TPZFileStream fstr;
-                std::stringstream ss;
-                int frictionlaw,levelmax;
-                this->parameters->FindParam(&frictionlaw,FrictionLawEnum);
-                this->parameters->FindParam(&levelmax,AmrLevelMaxEnum);
-                ss << levelmax;
-                if(frictionlaw==1){
-                        ss << "_viscous/amr.txt";
-                }else if(frictionlaw==7){
-                        ss << "_tsai/amr.txt";
-                }else{
-                        _error_("friction law not supported here.");
-                }
-                std::string AMRfile  = "/home/santos/Misomip2/L" + ss.str();
-                fstr.OpenWrite(AMRfile.c_str());
-                int withclassid = 1;
-                this->amr->Write(fstr,withclassid);
-        }
-        //itapopo
-        
-        /*Cleanup*/
-        xDelete<IssmDouble>(x);
-        xDelete<IssmDouble>(y);
-        xDelete<IssmDouble>(z);
-        xDelete<int>(elementslist);
-
-        return;
 }
 /*}}}*/
@@ -3342 +2498 @@
         }
         input_interpolations->Assemble();
-        
+
         /*Serialize and set output*/
         IssmDouble* input_interpolations_serial = input_interpolations->ToMPISerial();
@@ -3374 +2530 @@
                                                 P1input_interp[numP1inputs] = inputinterp;
                                         }
-                                         numP1inputs++;
+                                        numP1inputs++;
                                         break;
                                 case P0Enum:
@@ -3421 +2577 @@
                 xDelete<IssmDouble>(vector);
         }
-        
+
         /*Old mesh coordinates*/
         IssmDouble *Xold     = NULL;
@@ -3433 +2589 @@
         IssmDouble* YC_new   = NULL;
         int        *Indexnew = NULL;
-        
+
         /*Get the old mesh*/
         this->GetMesh(this->vertices,this->elements,&Xold,&Yold,&Zold,&Indexold);
@@ -3459 +2615 @@
                                 P0inputsold,numelementsold,numP0inputs,
                                 XC_new,YC_new,numelementsnew,NULL);
-        
+
         /*Interpolate P1 inputs in the new mesh*/
         InterpFromMeshToMesh2dx(&P1inputsnew,Indexold,Xold,Yold,numverticesold,numelementsold,
                                 P1inputsold,numverticesold,numP1inputs,
                                 Xnew,Ynew,numverticesnew,NULL);
-        
+
         /*Insert P0 inputs into the new elements.*/
         vector=NULL;
         for(int i=0;i<numP0inputs;i++){
-                
+
                 /*Get P0 input vector from the interpolated matrix*/
                 vector=xNew<IssmDouble>(numelementsnew);
@@ -3533 +2689 @@
         xDelete<IssmDouble>(YC_new);
         xDelete<int>(Indexnew);
-
-}
-/*}}}*/
-void FemModel::ExecuteRefinement(int &numberofvertices,int &numberofelements,IssmDouble** px,IssmDouble** py,IssmDouble** pz,int** pelementslist){/*{{{*/
-        
-        /*elements is in Matlab indexing*/
-        int my_rank                                             = IssmComm::GetRank();
-        int numberofsegments                    = -1;
-        IssmDouble* vx                                  = NULL; //itapopo this is not being used
-        IssmDouble* vy                                  = NULL; //itapopo this is not being used
-        IssmDouble* x                                   = NULL;
-        IssmDouble* y                                   = NULL;
-        IssmDouble* z                                   = NULL;
-        int* elementslist                               = NULL;
-        int* segments                                   = NULL;
-        IssmDouble* masklevelset        = NULL;
-        IssmDouble* pelementerror       = NULL;
-   const int elementswidth              = this->GetElementsWidth();//just 2D mesh, tria elements
-        
-        /*Solutions which will be used to refine the elements*/
-        this->GetGroundediceLevelSet(&masklevelset);//itapopo verificar se já existe um método igual a esse
-
-        //Compute the ZZ error estimator per element
-        this->ZZErrorEstimator(&pelementerror);
-        
-        _printf0_("P Element error\n");
-        for(int i=0;i<this->elements->NumberOfElements();i++)   _printf0_(""<<pelementerror[i]<< "\n");
-        _printf0_("\n");
-
-        if(my_rank==0){
-                int type_process=1; //1: it refines father mesh. See AdaptiveMeshRefinement.h (.cpp)
-                this->amr->ExecuteRefinement(type_process,vx,vy,masklevelset,
-                                                                                                numberofvertices,numberofelements,numberofsegments,&x,&y,&z,&elementslist,&segments);
-                if(numberofvertices<=0 || numberofelements<=0 /*|| newnumberofsegments<=0*/) _error_("Error in the refinement process.");
-        }
-        else{
-                x=xNew<IssmDouble>(numberofvertices);
-                y=xNew<IssmDouble>(numberofvertices);
-                z=xNew<IssmDouble>(numberofvertices);
-                elementslist=xNew<int>(numberofelements*this->GetElementsWidth());
-        }
-
-        /*Send new mesh to others CPU*/
-        ISSM_MPI_Bcast(&numberofvertices,1,ISSM_MPI_INT,0,IssmComm::GetComm());
-        ISSM_MPI_Bcast(&numberofelements,1,ISSM_MPI_INT,0,IssmComm::GetComm());
-        ISSM_MPI_Bcast(x,numberofvertices,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());      
-        ISSM_MPI_Bcast(y,numberofvertices,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());      
-        ISSM_MPI_Bcast(z,numberofvertices,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());      
-        ISSM_MPI_Bcast(elementslist,numberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm());      
-
-        /*Assign the pointers*/ 
-        (*pelementslist) = elementslist; //Matlab indexing
-        (*px)                             = x;
-        (*py)                             = y;
-        (*pz)                             = z;
-
+}
+/*}}}*/
+void FemModel::WriteMeshInResults(void){/*{{{*/
+
+        int step                                        = -1;
+        int numberofelements = -1;
+        int numberofvertices = -1;
+        IssmDouble time         = -1;
+        IssmDouble* x                   = NULL;
+        IssmDouble* y                   = NULL;
+        IssmDouble* z                   = NULL;
+        int* elementslist               = NULL;
+
+        if(!this->elements || !this->vertices || !this->results || !this->parameters) return;
+         
+        parameters->FindParam(&step,StepEnum);
+        parameters->FindParam(&time,TimeEnum);
+        numberofelements=this->elements->NumberOfElements();
+        numberofvertices=this->vertices->NumberOfVertices();
+
+        /*Get mesh. Elementslist comes in Matlab indexing*/
+        this->GetMesh(this->vertices,this->elements,&x,&y,&z,&elementslist);
+
+        /*Write mesh in Results*/
+        printf("-------------- file: FemModel.cpp line: %i\n",__LINE__); 
+        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));
+        
+        //itapopo
+        #ifdef _HAVE_NEOPZ_
+        if(IssmComm::GetRank()==0){
+                TPZFileStream fstr;
+                std::stringstream ss;
+                int frictionlaw,levelmax;
+                this->parameters->FindParam(&frictionlaw,FrictionLawEnum);
+                this->parameters->FindParam(&levelmax,AmrLevelMaxEnum);
+                ss << levelmax;
+                if(frictionlaw==1){
+                        ss << "_viscous/amr.txt";
+                }else if(frictionlaw==7){
+                        ss << "_tsai/amr.txt";
+                }else{
+                        _error_("friction law not supported here.");
+                }
+                std::string AMRfile  = "/home/santos/Misomip2/L" + ss.str();
+                fstr.OpenWrite(AMRfile.c_str());
+                int withclassid = 1;
+                this->amr->Write(fstr,withclassid);
+        }
+        #endif
+        //itapopo
+        
         /*Cleanup*/
-        if(segments) xDelete<int>(segments);
-        xDelete<IssmDouble>(masklevelset);
-        xDelete<IssmDouble>(pelementerror);
-
+        xDelete<IssmDouble>(x);
+        xDelete<IssmDouble>(y);
+        xDelete<IssmDouble>(z);
+        xDelete<int>(elementslist);
 }
 /*}}}*/
@@ -3635 +2793 @@
         connectivity=xNewZeroInit<int>(newnumberofvertices);
 
-        for (int i=0;i<newnumberofelements;i++){
-                for (int j=0;j<elementswidth;j++){
+        for(int i=0;i<newnumberofelements;i++){
+                for(int j=0;j<elementswidth;j++){
                         int vertexid = newelementslist[elementswidth*i+j];
                         _assert_(vertexid>0 && vertexid-1<newnumberofvertices);//Matlab indexing
@@ -3722 +2880 @@
         newmatpar->SetMid(newnumberofelements+1);
         materials->AddObject(newmatpar);//put it at the end of the materials        
-
 }
 /*}}}*/
@@ -3954 +3111 @@
         delete vspcvxflag;
         delete vspcvyflag;
-
 }
 /*}}}*/
@@ -4053 +3209 @@
         xDelete<int>(npart);        
         xDelete<int>(index);
-
 }
 /*}}}*/
@@ -4141 +3296 @@
    xDelete<IssmDouble>(totalweight);
    xDelete<int>(elem_vertices);
-
 }
 /*}}}*/
@@ -4215 +3369 @@
         xDelete<int>(elem_vertices);
         delete velementerror;
+}
+/*}}}*/
+
+#ifdef  _HAVE_DAKOTA_
+void FemModel::DakotaResponsesx(double* d_responses,char** responses_descriptors,int numresponsedescriptors,int d_numresponses){/*{{{*/
+
+        int        i,j;
+        int        my_rank;
+
+        /*intermediary: */
+        char   root[50];
+        int    index;
+        int    npart;
+        double femmodel_response;
+        int    flag;
+        double *vertex_response   = NULL;
+        double *qmu_response      = NULL;
+        double *responses_pointer = NULL;
+
+        /*retrieve npart: */
+        parameters->FindParam(&npart,QmuNumberofpartitionsEnum);
+        my_rank=IssmComm::GetRank();
+
+        /*save the d_responses pointer: */
+        responses_pointer=d_responses;
+
+        //watch out, we have more d_numresponses than numresponsedescriptors, because the responses have been expanded if they were scaled. 
+        //because we don't know the d_responses descriptors (the scaled ones) we can't key off them, so we will key off the responses_descriptors: */
+
+        for(i=0;i<numresponsedescriptors;i++){
+
+                flag=DescriptorIndex(root,&index,responses_descriptors[i]);
+
+                if(flag==ScaledEnum){
+
+                        /*this response was scaled. pick up the response from the inputs: */
+                        GetVectorFromInputsx(&vertex_response,this, StringToEnumx(root),VertexPIdEnum);
+
+                        /*Now, average it onto the partition nodes: */
+                        AverageOntoPartitionx(&qmu_response,elements,nodes,vertices,loads,materials,parameters,vertex_response);
+
+                        /*Copy onto our dakota responses: */
+                        if(my_rank==0){
+                                /*plug response: */
+                                for(j=0;j<npart;j++)responses_pointer[j]=qmu_response[j];
+
+                                /*increment response_pointer :*/
+                                responses_pointer+=npart;
+                        }
+
+                        /*Free ressources:*/
+                        xDelete<double>(vertex_response);
+                        xDelete<double>(qmu_response);
+
+                }
+                else if (flag==IndexedEnum){
+
+                        /*indexed response: plug index into parameters and call response module: */
+                        parameters->SetParam(index,IndexEnum);
+
+                        this->Responsex(&femmodel_response,root);
+
+                        if(my_rank==0){
+                                /*plug response: */
+                                responses_pointer[0]=femmodel_response;
+
+                                /*increment response_pointer :*/
+                                responses_pointer++;
+                        }
+                }
+                else if (flag==NodalEnum){
+                        _error_("nodal response functions not supported yet!");
+
+                        /*increment response_pointer :*/
+                        responses_pointer++;
+                }
+                else if (flag==RegularEnum){
+
+                        /*perfectly normal response function: */
+                        this->Responsex(&femmodel_response,root);
+
+                        if(my_rank==0){
+                                /*plug response: */
+                                responses_pointer[0]=femmodel_response;
+
+                                /*increment response_pointer :*/
+                                responses_pointer++;
+                        }
+                }
+                else _error_("flag type " << flag << " not supported yet for response analysis");
+        }
+
+        /*Synthesize echo: {{{*/
+        if(my_rank==0){
+                _printf_("   responses: " << d_numresponses << ": ");
+                for(i=0;i<d_numresponses-1;i++)_printf_(d_responses[i] << "|");
+                _printf_(d_responses[d_numresponses-1]);
+                _printf_("\n");
+        }
+        /*}}}*/
 
 }
 /*}}}*/
 #endif
+#ifdef _HAVE_GIAIVINS_
+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, 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,xx,yy);
+                }
+                if(i%100==0){
+                        pUp->Assemble();
+                        pNorth->Assemble();
+                        pEast->Assemble();
+                }
+        }
+        
+        /*One last time: */
+        pUp->Assemble();
+        pNorth->Assemble();
+        pEast->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){/*{{{*/
+
+        IssmDouble  eartharea=0;
+        IssmDouble  eartharea_cpu=0;
+
+        int         ns,nsmax;
+        
+        /*Go through elements, and add contribution from each element to the deflection vector wg:*/
+        ns = elements->Size();
+        
+        /*First, figure out the surface area of Earth: */ 
+        for(int i=0;i<ns;i++){
+                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
+                eartharea_cpu += element->GetAreaSpherical();
+        }
+        ISSM_MPI_Reduce (&eartharea_cpu,&eartharea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
+        ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+
+        /*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,eartharea);
+                }
+                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>* pSgi, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius) { /*{{{*/
+
+        /*serialized vectors:*/
+        IssmDouble  eustatic       = 0.;
+        IssmDouble  eustatic_cpu   = 0.;
+        IssmDouble  eustatic_cpu_e = 0.;
+        IssmDouble  oceanarea      = 0.;
+        IssmDouble  oceanarea_cpu  = 0.;
+        IssmDouble  eartharea      = 0.;
+        IssmDouble  eartharea_cpu  = 0.;
+        int         ns,nsmax;
+        
+        /*Go through elements, and add contribution from each element to the deflection vector wg:*/
+        ns = elements->Size();
+
+        /*First, figure out the area of the ocean, which is needed to compute the eustatic component: */
+        for(int i=0;i<ns;i++){
+                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
+                oceanarea_cpu += element->OceanArea();
+                eartharea_cpu += element->GetAreaSpherical();
+        }
+        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.);
+
+        ISSM_MPI_Reduce (&eartharea_cpu,&eartharea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
+        ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+
+        /*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 sea level rise core: */
+        for(int i=0;i<nsmax;i++){
+                if(i<ns){
+                
+                        if(VerboseConvergence())if(i%100==0)_printf0_("\r" << "      convolution progress: " << (double)i/(double)ns*100 << "%  ");
+                
+                        Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
+                        element->SealevelriseEustatic(pSgi,&eustatic_cpu_e,latitude,longitude,radius,oceanarea,eartharea);
+                        eustatic_cpu+=eustatic_cpu_e;
+                }
+                if(i%100==0)pSgi->Assemble();
+        }
+        if(VerboseConvergence())_printf0_("\n");
+                
+        /*One last time: */
+        pSgi->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));
+
+        /*Assign output pointers:*/
+        *peustatic=eustatic;
+
+}
+/*}}}*/
+void FemModel::SealevelriseNonEustatic(Vector<IssmDouble>* pSgo, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, bool verboseconvolution){/*{{{*/
+
+        /*serialized vectors:*/
+        IssmDouble* Sg_old=NULL;
+        
+        IssmDouble  eartharea=0;
+        IssmDouble  eartharea_cpu=0;
+
+        int         ns,nsmax;
+        
+        /*Serialize vectors from previous iteration:*/
+        Sg_old=pSg_old->ToMPISerial();
+
+        /*Go through elements, and add contribution from each element to the deflection vector wg:*/
+        ns = elements->Size();
+
+        /*First, figure out the area of the ocean, which is needed to compute the eustatic component: */
+        for(int i=0;i<ns;i++){
+                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
+                eartharea_cpu += element->GetAreaSpherical();
+        }
+        
+        ISSM_MPI_Reduce (&eartharea_cpu,&eartharea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
+        ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+
+        /*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 sea level rise core: */
+        for(int i=0;i<nsmax;i++){
+                if(i<ns){
+                        if(verboseconvolution)if(VerboseConvergence())if(i%100==0)_printf_("\r" << "      convolution progress: " << (double)i/(double)ns*100 << "%   ");
+                        Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
+                        element->SealevelriseNonEustatic(pSgo,Sg_old,latitude,longitude,radius,eartharea);
+                }
+                if(i%100==0)pSgo->Assemble();
+        }
+        if(verboseconvolution)if(VerboseConvergence())_printf_("\n");
+        
+        /*Free ressources:*/
+        xDelete<IssmDouble>(Sg_old);
+}
+/*}}}*/
+void FemModel::SealevelriseRotationalFeedback(Vector<IssmDouble>* pSgo_rot, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius){/*{{{*/
+
+        /*serialized vectors:*/
+        IssmDouble* Sg_old=NULL;
+        IssmDouble  eartharea=0;
+        IssmDouble  eartharea_cpu=0;
+        IssmDouble      tide_love_h, tide_love_k, fluid_love, moi_e, moi_p, omega, g;
+        IssmDouble      load_love_k2 = -0.30922675; //degree 2 load Love number 
+        IssmDouble      m1, m2, m3; 
+        IssmDouble      lati, longi, radi, value; 
+
+        /*Serialize vectors from previous iteration:*/
+        Sg_old=pSg_old->ToMPISerial();
+
+        /*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));
+                eartharea_cpu += element->GetAreaSpherical();
+        }
+        ISSM_MPI_Reduce (&eartharea_cpu,&eartharea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
+        ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+
+        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],Sg_old,eartharea);
+                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(&tide_love_h,SealevelriseTidalLoveHEnum);
+        parameters->FindParam(&tide_love_k,SealevelriseTidalLoveKEnum);
+        parameters->FindParam(&fluid_love,SealevelriseFluidLoveEnum);
+        parameters->FindParam(&moi_e,SealevelriseEquatorialMoiEnum);
+        parameters->FindParam(&moi_p,SealevelrisePolarMoiEnum);
+        parameters->FindParam(&omega,SealevelriseAngularVelocityEnum);
+
+        /*compute perturbation terms for angular velocity vector: */
+        m1 = 1/(1-tide_love_k/fluid_love) * (1+load_love_k2)/(moi_p-moi_e) * moi_list[0]; 
+        m2 = 1/(1-tide_love_k/fluid_love) * (1+load_love_k2)/(moi_p-moi_e) * moi_list[1]; 
+        m3 = -(1+load_love_k2)/moi_p * moi_list[2];     // term associated with fluid number (3-order-of-magnitude smaller) is negelected  
+
+        /* 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-tide_love_h)/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))); 
+        
+                pSgo_rot->SetValue(sid,value,INS_VAL); //INS_VAL ensures that you don't add several times
+        }
+
+        /*Assemble mesh velocity*/
+        pSgo_rot->Assemble();
+        
+        /*Free ressources:*/
+        xDelete<IssmDouble>(Sg_old);
+        
+}
+/*}}}*/
+void FemModel::SealevelriseGeodetic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pSg, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){/*{{{*/
+
+        /*serialized vectors:*/
+        IssmDouble* Sg=NULL;
+        
+        IssmDouble  eartharea=0;
+        IssmDouble  eartharea_cpu=0;
+
+        int         ns,nsmax;
+        
+        /*Serialize vectors from previous iteration:*/
+        Sg=pSg->ToMPISerial();
+
+        /*Go through elements, and add contribution from each element to the deflection vector wg:*/
+        ns = elements->Size();
+        
+        /*First, figure out the area of the ocean, which is needed to compute the eustatic component: */
+        for(int i=0;i<ns;i++){
+                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
+                eartharea_cpu += element->GetAreaSpherical();
+        }
+        ISSM_MPI_Reduce (&eartharea_cpu,&eartharea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
+        ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+
+        /*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 sea level rise core: */
+        for(int i=0;i<nsmax;i++){
+                if(i<ns){
+                        Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
+                        element->SealevelriseGeodetic(pUp,pNorth,pEast,Sg,latitude,longitude,radius,xx,yy,zz,eartharea);
+                }
+                if(i%100==0){
+                        pUp->Assemble();
+                        pNorth->Assemble();
+                        pEast->Assemble();
+                }
+        }
+        
+        /*One last time: */
+        pUp->Assemble();
+        pNorth->Assemble();
+        pEast->Assemble();
+
+        /*Free ressources:*/
+        xDelete<IssmDouble>(Sg);
+        xDelete<IssmDouble>(latitude);
+        xDelete<IssmDouble>(longitude);
+        xDelete<IssmDouble>(radius);
+        xDelete<IssmDouble>(xx);
+        xDelete<IssmDouble>(yy);
+        xDelete<IssmDouble>(zz);
+}
+/*}}}*/
+IssmDouble FemModel::SealevelriseOceanAverage(Vector<IssmDouble>* Sg) { /*{{{*/
+
+        IssmDouble* Sg_serial=NULL;
+        IssmDouble  oceanvalue,oceanvalue_cpu;
+        IssmDouble  oceanarea,oceanarea_cpu;
+
+        /*Serialize vectors from previous iteration:*/
+        Sg_serial=Sg->ToMPISerial();
+
+        /*Initialize:*/
+        oceanvalue_cpu=0;
+        oceanarea_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));
+                oceanarea_cpu += element->OceanArea();
+                oceanvalue_cpu += element->OceanAverage(Sg_serial);
+        }
+        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());
+        
+        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>(Sg_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*/
+        mask=new Vector<IssmDouble>(this->nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum));
+        recurence=new Vector<IssmDouble>(this->nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum));
+        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,eplzigzag_counter,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(HydrologyDCEfficientAnalysisEnum));
+        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;
+        for (int i=0;i<nodes->Size();i++){
+                Node* node=xDynamicCast<Node*>(nodes->GetObjectByOffset(i));
+                if(node->InAnalysis(HydrologyDCEfficientAnalysisEnum)){
+                        if(serial_active[node->Sid()]==1.){
+                                node->Activate();
+                                if(!node->IsClone()) counter++;
+                        }
+                        else{
+                                node->Deactivate();
+                        }
+                }
+        }
+        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()) _printf0_("   Number of active nodes in EPL layer: "<< counter <<"\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*/
+        active=new Vector<IssmDouble>(nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum));
+        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;
+        for (int i=0;i<nodes->Size();i++){
+                Node* node=xDynamicCast<Node*>(nodes->GetObjectByOffset(i));
+                if(node->InAnalysis(L2ProjectionEPLAnalysisEnum)){
+                        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");
+}
+/*}}}*/
+#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->Tag(START);
+
+        /*From command line arguments, retrieve different filenames needed to create the FemModel: */
+        solution_type=StringToEnumx(solution);
+        
+        /*Create femmodel from input files: */
+        profiler->Tag(STARTINIT);
+        this->InitFromBuffers((char*)buffer,buffersize,toolkits, solution_type,trace,NULL);
+        profiler->Tag(FINISHINIT);
+        
+        /*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->Tag(FINISH);  //final tagging
+        _printf0_("\n");
+        _printf0_("   "<<setw(40)<<left<<"FemModel initialization elapsed time:"<<profiler->DeltaTime(STARTINIT,FINISHINIT) << "\n");
+        _printf0_("   "<<setw(40)<<left<<"Core solution elapsed time:"<<profiler->DeltaTime(STARTCORE,FINISHCORE) << "\n");
+        _printf0_("\n");
+        _printf0_("   Total elapsed time: "
+                                <<profiler->DeltaTimeModHour(START,FINISH)<<" hrs "
+                                <<profiler->DeltaTimeModMin(START,FINISH)<<" min "
+                                <<profiler->DeltaTimeModSec(START,FINISH)<<" sec"
+                                );
+        _printf0_("\n");
+        
+        /*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
+
+#ifdef _HAVE_BAMG_
+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;
+        int *newelementslist            = NULL;
+        int newnumberofvertices = -1;
+        int newnumberofelements = -1;
+
+        /*Get Rank*/
+        int my_rank     = IssmComm::GetRank();
+
+        if(my_rank==0){
+                this->amrbamg->ExecuteRefinementBamg(&newnumberofvertices,&newnumberofelements,&newx,&newy,&newz,&newelementslist);
+                if(newnumberofvertices<=0 || newnumberofelements<=0) _error_("Error in the refinement process.");
+        }
+        else{
+                newx=xNew<IssmDouble>(newnumberofvertices);
+                newy=xNew<IssmDouble>(newnumberofvertices);
+                newz=xNew<IssmDouble>(newnumberofvertices);
+                newelementslist=xNew<int>(newnumberofelements*this->GetElementsWidth());
+        }
+
+        /*Send new mesh to others CPU*/
+        ISSM_MPI_Bcast(&newnumberofvertices,1,ISSM_MPI_INT,0,IssmComm::GetComm());
+        ISSM_MPI_Bcast(&newnumberofelements,1,ISSM_MPI_INT,0,IssmComm::GetComm());
+        ISSM_MPI_Bcast(newx,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 
+        ISSM_MPI_Bcast(newy,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 
+        ISSM_MPI_Bcast(newz,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 
+        ISSM_MPI_Bcast(newelementslist,newnumberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm());        
+
+        /*Assign output pointers*/
+        *pnewnumberofvertices = newnumberofvertices;
+        *pnewnumberofelements = newnumberofelements;
+        *pnewx = newx;
+        *pnewy = newy;
+        *pnewz = newz;
+        *pnewelementslist = newelementslist;
+}
+/*}}}*/
+void FemModel::InitializeAdaptiveRefinementBamg(void){/*{{{*/
+
+        /*Define variables*/
+        int numberofvertices      = this->vertices->NumberOfVertices();
+        int numberofelements      = this->elements->NumberOfElements();
+        int numberofsegments      = 0; //used on matlab
+        IssmDouble* x             = NULL;
+        IssmDouble* y             = NULL;
+        IssmDouble* z             = NULL;
+        int* elements             = NULL;
+        int elementswidth         = this->GetElementsWidth(); //just tria elements in this version. Itapopo:
+        IssmDouble hmin,hmax,err;
+        int        fieldenum;
+
+   /*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,&z,&elements);
+
+        /*Get amr parameters*/
+        this->parameters->FindParam(&hmin,AmrHminEnum);
+        this->parameters->FindParam(&hmax,AmrHmaxEnum);
+        this->parameters->FindParam(&fieldenum,AmrFieldEnum);
+        this->parameters->FindParam(&err,AmrErrEnum);
+
+        /*Create bamg data structures for bamg (only cpu 0)*/
+        if(my_rank==0){ 
+                /*Re-create original mesh and put it in bamg structure*/
+                this->amrbamg = new AmrBamg(hmin,hmax);
+                this->amrbamg->Initialize(elements,x,y,numberofvertices,numberofelements);
+        }
+
+        /*Free the vectors*/
+        xDelete<IssmDouble>(x);
+        xDelete<IssmDouble>(y);
+        xDelete<IssmDouble>(z);
+        xDelete<int>(elements);
+}
+/*}}}*/
+#endif
+
+#ifdef _HAVE_NEOPZ_
+void FemModel::ReMeshNeopz(int* pnumberofvertices,int* pnumberofelements,IssmDouble** px,IssmDouble** py,IssmDouble** pz,int** pelementslist){/*{{{*/
+
+        /*elements is in Matlab indexing*/
+        int my_rank                                             = IssmComm::GetRank();
+        int numberofsegments                    = -1;
+        int numberofvertices,numberofelements;
+        IssmDouble* vx                                  = NULL; //itapopo this is not being used
+        IssmDouble* vy                                  = NULL; //itapopo this is not being used
+        IssmDouble* x                                   = NULL;
+        IssmDouble* y                                   = NULL;
+        IssmDouble* z                                   = NULL;
+        int* elementslist                               = NULL;
+        int* segments                                   = NULL;
+        IssmDouble* masklevelset        = NULL;
+        IssmDouble* pelementerror       = NULL;
+        const int elementswidth         = this->GetElementsWidth();//just 2D mesh, tria elements
+
+        /*Solutions which will be used to refine the elements*/
+        this->GetGroundediceLevelSet(&masklevelset);//itapopo verificar se já existe um método igual a esse
+
+        //Compute the ZZ error estimator per element
+        this->ZZErrorEstimator(&pelementerror);
+
+        _printf0_("P Element error\n");
+        for(int i=0;i<this->elements->NumberOfElements();i++)   _printf0_(""<<pelementerror[i]<< "\n");
+        _printf0_("\n");
+
+        if(my_rank==0){
+                int type_process=1; //1: it refines father mesh. See AdaptiveMeshRefinement.h (.cpp)
+                this->amr->ExecuteRefinement(type_process,vx,vy,masklevelset,
+                                        numberofvertices,numberofelements,numberofsegments,&x,&y,&z,&elementslist,&segments);
+                if(numberofvertices<=0 || numberofelements<=0 /*|| newnumberofsegments<=0*/) _error_("Error in the refinement process.");
+        }
+        else{
+                x=xNew<IssmDouble>(numberofvertices);
+                y=xNew<IssmDouble>(numberofvertices);
+                z=xNew<IssmDouble>(numberofvertices);
+                elementslist=xNew<int>(numberofelements*this->GetElementsWidth());
+        }
+
+        /*Send new mesh to others CPU*/
+        ISSM_MPI_Bcast(&numberofvertices,1,ISSM_MPI_INT,0,IssmComm::GetComm());
+        ISSM_MPI_Bcast(&numberofelements,1,ISSM_MPI_INT,0,IssmComm::GetComm());
+        ISSM_MPI_Bcast(x,numberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());       
+        ISSM_MPI_Bcast(y,numberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());       
+        ISSM_MPI_Bcast(z,numberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());       
+        ISSM_MPI_Bcast(elementslist,numberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm());      
+
+        /*Assign the pointers*/ 
+        (*pelementslist) = elementslist; //Matlab indexing
+        (*px)                             = x;
+        (*py)                             = y;
+        (*pz)                             = z;
+        *pnumberofelements = numberofelements;
+        *pnumberofvertices = numberofvertices;
+
+        /*Cleanup*/
+        if(segments) xDelete<int>(segments);
+        xDelete<IssmDouble>(masklevelset);
+        xDelete<IssmDouble>(pelementerror);
+
+}
+/*}}}*/
+void FemModel::InitializeAdaptiveRefinementNeopz(void){/*{{{*/
+        
+        /*Define variables*/
+        int my_rank                                             = IssmComm::GetRank();
+        int numberofvertices                    = this->vertices->NumberOfVertices();
+        int numberofelements                    = this->elements->NumberOfElements();
+        int numberofsegments                    = 0; //used on matlab
+        IssmDouble* x                                   = NULL;
+        IssmDouble* y                                   = NULL;
+        IssmDouble* z                                   = NULL;
+        int* elements                                   = NULL;
+        int elementswidth                               = this->GetElementsWidth(); //just tria elements in this version. Itapopo:
+        int levelmax                                    = 0;
+        IssmDouble regionlevel1         = 0.;
+        IssmDouble regionlevelmax       = 0.;
+
+        /*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,&z,&elements);
+        
+        /*Get amr parameters*/
+        this->parameters->FindParam(&levelmax,AmrLevelMaxEnum);
+        this->parameters->FindParam(&regionlevel1,AmrRegionLevel1Enum);
+        this->parameters->FindParam(&regionlevelmax,AmrRegionLevelMaxEnum);
+
+        /*Create initial mesh (coarse mesh) in neopz data structure*/ 
+        /*Just CPU #0 should keep AMR object*/
+   this->SetRefPatterns();
+        if(my_rank==0){ 
+           bool ismisomip       = false;
+                if(ismisomip){//itapopo
+                        TPZFileStream fstr;
+                        std::stringstream ss;
+                        int frictionlaw;
+                        
+                        this->parameters->FindParam(&frictionlaw,FrictionLawEnum);
+                
+                        ss      << levelmax;
+                        if(frictionlaw==1){
+                                ss << "_viscous/amr.txt";
+                        }else if(frictionlaw==7){
+                                ss << "_tsai/amr.txt";
+                        }else{
+                                _error_("friction law not supported here.");
+                        }
+                        
+                        std::string AMRfile  = "/home/santos/Misomip2/L" + ss.str();
+                        fstr.OpenRead(AMRfile.c_str());
+                        
+                        TPZSaveable *sv         = TPZSaveable::Restore(fstr,0);
+                        this->amr                               = dynamic_cast<AdaptiveMeshRefinement*>(sv);
+                }
+                else{
+                        this->amr = new AdaptiveMeshRefinement();
+                        //this->amr->SetLevelMax(levelmax); //Set max level of refinement
+                        //this->amr->SetRegions(regionlevel1,regionlevelmax);
+                        this->amr->CreateInitialMesh(numberofvertices,numberofelements,numberofsegments,elementswidth,x,y,z,elements,NULL);
+                }
+                this->amr->SetLevelMax(levelmax); //Set max level of refinement
+                this->amr->SetRegions(regionlevel1,regionlevelmax);
+        }
+
+        /*Free the vectors*/
+        xDelete<IssmDouble>(x);
+        xDelete<IssmDouble>(y);
+        xDelete<IssmDouble>(z);
+        xDelete<int>(elements);
+}
+/*}}}*/
+void FemModel::SetRefPatterns(){/*{{{*/
+
+   /*Initialize the global variable of refinement patterns*/
+   gRefDBase.InitializeUniformRefPattern(EOned);
+   gRefDBase.InitializeUniformRefPattern(ETriangle);
+
+    //gRefDBase.InitializeRefPatterns();
+   /*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_OnlyTriang_Side_3_4.rpt";
+   std::string filename4 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_4_permuted.rpt";
+   std::string filename5 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_5.rpt";
+   std::string filename6 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_5_permuted.rpt";
+   std::string filename7 = filepath + "/2D_Triang_Rib_5.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);
+
+   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);
+}
+/*}}}*/
+#endif

issm/trunk-jpl/src/c/classes/FemModel.h

@@ -22 +22 @@
 #ifdef _HAVE_NEOPZ_
 #include "./AdaptiveMeshRefinement.h"
+#endif
+#ifdef _HAVE_BAMG_
+#include "./AmrBamg.h"
 #endif
 /*}}}*/
@@ -47 +50 @@
                 Vertices    *vertices;             //one set of vertices
 
+                //FIXME: do we want only one class and have virtual functions? or keep 2 classes, at least rename AdaptiveMeshRefinement -> AmrNeopz
                 #ifdef _HAVE_NEOPZ_
                 AdaptiveMeshRefinement *amr;              //adaptive mesh refinement object. It keeps coarse mesh and execute refinement process
+                #endif
+
+                #ifdef _HAVE_BAMG_
+                AmrBamg *amrbamg; //adaptive mesh refinement object. It keeps coarse mesh and execute refinement process
                 #endif
 
@@ -61 +69 @@
                 FemModel* copy();
                 void Echo();
+                int  GetElementsWidth(){return 3;};//just tria elements in this first version
                 void InitFromFiles(char* rootpath, char* inputfilename, char* outputfilename, char* petscfilename, char* lockfilename, char* restartfilename, const int solution_type,bool trace,IssmPDouble* X=NULL);
                 void InitFromFids(char* rootpath, FILE* IOMODEL, FILE* toolkitsoptionsfid, int in_solution_type, bool trace, IssmPDouble* X=NULL);
                 void Marshall(char** pmarshalled_data, int* pmarshalled_data_size, int marshall_direction);
+                void ReMesh(void);
                 void Restart(void);
                 void SetCurrentConfiguration(int configuration_type);
@@ -146 +156 @@
                 #endif
 
-                #ifdef _HAVE_NEOPZ_
-                /*Adaptive mesh refinement methods*/
-                void InitializeAdaptiveRefinement(void);
-                void ReMesh(void);
+                /*AMR*/
                 void BedrockFromMismipPlus(void);
                 void AdjustBaseThicknessAndMask(void);
                 void GetMesh(Vertices* femmodel_vertices,Elements* femmodel_elements,IssmDouble** px, IssmDouble** py, IssmDouble** pz, int** pelementslist);
-                int GetElementsWidth(){return 3;};//just tria elements in this first version
-                void ExecuteRefinement(int &numberofvertices,int &numberofelements,IssmDouble** px,IssmDouble** py,IssmDouble** pz,int** pelementslist);
                 void GetGroundediceLevelSet(IssmDouble** pmasklevelset);
                 void CreateVertices(int newnumberofvertices,int newnumberofelements,int elementswidth,int* newelementslist,int* my_vertices,IssmDouble* newx,IssmDouble* newy,IssmDouble* newz,Vertices* vertices);
@@ -168 +173 @@
                 void SmoothedDeviatoricStressTensor(IssmDouble** ptauxx,IssmDouble** ptauyy,IssmDouble** ptauxy); //nodal values, just for SSA-P1: TauXX, TauYY, TauXY
                 void ZZErrorEstimator(IssmDouble** pelementerror);
+
+                #ifdef _HAVE_BAMG_
+                void ReMeshBamg(int* pnewnumberofvertices,int* pnewnumberofelements,IssmDouble** pnewx,IssmDouble** pnewy,IssmDouble** pnewz,int** pnewelementslist);
+                void InitializeAdaptiveRefinementBamg(void);
+                #endif
+
+                #ifdef _HAVE_NEOPZ_
+                /*Adaptive mesh refinement methods*/
+                void ReMeshNeopz(int* pnewnumberofvertices,int* pnewnumberofelements,IssmDouble** pnewx,IssmDouble** pnewy,IssmDouble** pnewz,int** pnewelementslist);
+                void InitializeAdaptiveRefinementNeopz(void);
                 #endif
 };
-                
 
 #endif

issm/trunk-jpl/src/c/classes/IoModel.cpp

@@ -605 +605 @@
                                         if(strcmp(record_name,"md.materials.type")==0) integer = IoCodeToEnumMaterials(integer);
                                         if(strcmp(record_name,"md.timestepping.type")==0) integer = IoCodeToEnumTimestepping(integer);
+                                        if(strcmp(record_name,"md.amr.type")==0) integer = IoCodeToEnumAmr(integer);
 
                                         /*Broadcast to other cpus*/

issm/trunk-jpl/src/c/cores/transient_core.cpp

@@ -111 +111 @@
                                 }
                         }
-                        if(VerboseSolution()) _printf0_("   computing thermal regime\n");
                         thermal_core(femmodel);
                 }
@@ -182 +181 @@
 
                 /*Adaptive mesh refinement*/
-                #ifdef _HAVE_NEOPZ_
                 if(amr_frequency){
                         if(save_results) femmodel->WriteMeshInResults();
-                        if(step%amr_frequency==0 && time<finaltime) femmodel->ReMesh();//Do not refine the last step
-                }
-                #endif
+                        if(step%amr_frequency==0 && time<finaltime){
+                                if(VerboseSolution()) _printf0_("   refining mesh\n");
+                                femmodel->ReMesh();//Do not refine the last step
+                        }
+                        //if(save_results){
+                        //      printf("-------------- file: transient_core.cpp line: %i\n",__LINE__); 
+                        //      femmodel->WriteMeshInResults();
+                        //      OutputResultsx(femmodel);
+                        //}
+                }
         }
         

issm/trunk-jpl/src/c/modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp

@@ -24 +24 @@
         R2     r;
         I2     I;
-        int    i,j;
+        int    i,j,k;
         int    it;
         int    i0,i1,i2;
@@ -64 +64 @@
                 if(y_data[i]<ymin) ymin=y_data[i];
                 if(y_data[i]>ymax) ymax=y_data[i];
+        }
+
+        /*Create Single vertex to element connectivity*/
+        int* connectivity = xNew<int>(nods_data);
+        for(i=0;i<nels_data;i++){
+                for(j=0;j<3;j++){
+                        k = index_data[i*3+j]-1;
+                        _assert_(k>=0 & k<nods_data);
+                        connectivity[k]=i;
+                }
         }
 
@@ -137 +147 @@
                         /*If we fall outside of the convex or outside of the mesh, return NaN*/
                         if(tb.det<0 || reft[it]<0){
-                                for (j=0;j<N_data;j++){
-                                        data_interp[i*N_data+j]=NAN;
-                                }
+                                _assert_(i0>=0 & i0<nods_data);
+                                it=connectivity[i0]; //or i1 or i2
+                                _assert_(it>=0 && it<nels_data);
+                                for(j=0;j<N_data;j++) data_interp[i*N_data+j]=data[N_data*it+j];
                         }
                         else{
+                                /*Inside the mesh!*/
                                 if(it<0 || it>=nels_data){
                                         _error_("Triangle number " << it << " not in [0 " << nels_data
                                                                 << "], report bug to developers (interpolation point: " <<x_interp[i]<<" "<<y_interp[i]<<")");
                                 }
-                                for (j=0;j<N_data;j++){
-
-                                        data_interp[i*N_data+j]=data[N_data*it+j];
-                                }
+                                for (j=0;j<N_data;j++) data_interp[i*N_data+j]=data[N_data*it+j];
                         }
                 }
@@ -158 +167 @@
         delete Th;
         xDelete<long>(reft);
+        xDelete<int>(connectivity);
         *pdata_interp=data_interp;
         return 1;

issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp

@@ -21 +21 @@
         int         numoutputs,materialtype,smb_model,basalforcing_model,timestepping_type;
         char**      requestedoutputs = NULL;
+        char*       fieldname = NULL;
         IssmDouble  time;
 
@@ -60 +61 @@
         parameters->AddObject(iomodel->CopyConstantObject("md.inversion.type",InversionTypeEnum));
         parameters->AddObject(iomodel->CopyConstantObject("md.calving.law",CalvingLawEnum));
-        {/*This is specific to ice...*/
+          {/*This is specific to ice...*/
                 parameters->AddObject(iomodel->CopyConstantObject("md.mesh.elementtype",MeshElementtypeEnum));
                 parameters->AddObject(iomodel->CopyConstantObject("md.steadystate.reltol",SteadystateReltolEnum));
@@ -81 +82 @@
                 parameters->AddObject(iomodel->CopyConstantObject("md.materials.rheology_law",MaterialsRheologyLawEnum));
                 parameters->AddObject(iomodel->CopyConstantObject("md.gia.cross_section_shape",GiaCrossSectionShapeEnum));
-                /*amr properties*/      
-                parameters->AddObject(iomodel->CopyConstantObject("md.amr.level_max",AmrLevelMaxEnum));
-                parameters->AddObject(iomodel->CopyConstantObject("md.amr.region_level_1",AmrRegionLevel1Enum));
-                parameters->AddObject(iomodel->CopyConstantObject("md.amr.region_level_max",AmrRegionLevelMaxEnum));
 
                 /*For stress balance only*/
@@ -95 +92 @@
                 if(iomodel->domaintype==Domain3DEnum)
                  parameters->AddObject(iomodel->CopyConstantObject("md.mesh.numberoflayers",MeshNumberoflayersEnum));
-        }
-
-        
+          }
+
+        /*amr properties*/      
+        int amrtype,amr_frequency;
+        iomodel->FindConstant(&amr_frequency,"md.transient.amr_frequency");
+        if(solution_type==TransientSolutionEnum && amr_frequency){
+                parameters->AddObject(iomodel->CopyConstantObject("md.amr.type",AmrTypeEnum));
+                iomodel->FindConstant(&amrtype,"md.amr.type");
+                switch(amrtype){
+                        #ifdef _HAVE_NEOPZ_
+                        case AmrNeopzEnum:
+                                parameters->AddObject(iomodel->CopyConstantObject("md.amr.level_max",AmrLevelMaxEnum));
+                                parameters->AddObject(iomodel->CopyConstantObject("md.amr.region_level_1",AmrRegionLevel1Enum));
+                                parameters->AddObject(iomodel->CopyConstantObject("md.amr.region_level_max",AmrRegionLevelMaxEnum));
+                                break;
+                        #endif
+
+                        #ifdef _HAVE_BAMG_
+                        case AmrBamgEnum:
+                                parameters->AddObject(iomodel->CopyConstantObject("md.amr.hmin",AmrHminEnum));
+                                parameters->AddObject(iomodel->CopyConstantObject("md.amr.hmax",AmrHmaxEnum));
+                                parameters->AddObject(iomodel->CopyConstantObject("md.amr.err",AmrErrEnum));
+                                /*Convert fieldname to enum and put it in params*/
+                                iomodel->FindConstant(&fieldname,"md.amr.fieldname");
+                                parameters->AddObject(new IntParam(AmrFieldEnum,StringToEnumx(fieldname)));
+                                xDelete<char>(fieldname);
+                                break;
+                        #endif
+
+                        default:
+                                _error_("Adaptive mesh refinement "<<EnumToStringx(amrtype)<<" not implemented yet");
+                }
+        }
+
         /*Basal forcing parameters*/
         parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.model",BasalforcingsEnum));
@@ -255 +283 @@
         ParseToolkitsOptionsx(parameters,toolkitsoptionsfid);
 
-        #ifdef _HAVE_ADOLC_
+        #ifdef _HAVE_ADOLC_
         if(VerboseMProcessor()) _printf0_("   starting autodiff parameters \n");
         CreateParametersAutodiff(parameters,iomodel);

issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

@@ -833 +833 @@
         /*}}}*/
         /*Adaptive mesh refinement (AMR){{{*/
+        TransientAmrFrequencyEnum,
+        AmrTypeEnum,
+        AmrNeopzEnum,
         AmrLevelMaxEnum,
         AmrRegionLevel1Enum,
         AmrRegionLevelMaxEnum,
-        TransientAmrFrequencyEnum,
+        AmrBamgEnum,
+        AmrHminEnum,
+        AmrHmaxEnum,
+        AmrFieldEnum,
+        AmrErrEnum,
         /*}}}*/
         ParametersENDEnum,

issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

@@ -809 +809 @@
                 case EsaRequestedOutputsEnum : return "EsaRequestedOutputs";
                 case EsaNumRequestedOutputsEnum : return "EsaNumRequestedOutputs";
+                case TransientAmrFrequencyEnum : return "TransientAmrFrequency";
+                case AmrTypeEnum : return "AmrType";
+                case AmrNeopzEnum : return "AmrNeopz";
                 case AmrLevelMaxEnum : return "AmrLevelMax";
                 case AmrRegionLevel1Enum : return "AmrRegionLevel1";
                 case AmrRegionLevelMaxEnum : return "AmrRegionLevelMax";
-                case TransientAmrFrequencyEnum : return "TransientAmrFrequency";
+                case AmrBamgEnum : return "AmrBamg";
+                case AmrHminEnum : return "AmrHmin";
+                case AmrHmaxEnum : return "AmrHmax";
+                case AmrFieldEnum : return "AmrField";
+                case AmrErrEnum : return "AmrErr";
                 case ParametersENDEnum : return "ParametersEND";
                 case XYEnum : return "XY";

issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

@@ -827 +827 @@
               else if (strcmp(name,"EsaRequestedOutputs")==0) return EsaRequestedOutputsEnum;
               else if (strcmp(name,"EsaNumRequestedOutputs")==0) return EsaNumRequestedOutputsEnum;
+              else if (strcmp(name,"TransientAmrFrequency")==0) return TransientAmrFrequencyEnum;
+              else if (strcmp(name,"AmrType")==0) return AmrTypeEnum;
+              else if (strcmp(name,"AmrNeopz")==0) return AmrNeopzEnum;
               else if (strcmp(name,"AmrLevelMax")==0) return AmrLevelMaxEnum;
               else if (strcmp(name,"AmrRegionLevel1")==0) return AmrRegionLevel1Enum;
               else if (strcmp(name,"AmrRegionLevelMax")==0) return AmrRegionLevelMaxEnum;
-              else if (strcmp(name,"TransientAmrFrequency")==0) return TransientAmrFrequencyEnum;
+              else if (strcmp(name,"AmrBamg")==0) return AmrBamgEnum;
+              else if (strcmp(name,"AmrHmin")==0) return AmrHminEnum;
+              else if (strcmp(name,"AmrHmax")==0) return AmrHmaxEnum;
+              else if (strcmp(name,"AmrField")==0) return AmrFieldEnum;
+              else if (strcmp(name,"AmrErr")==0) return AmrErrEnum;
               else if (strcmp(name,"ParametersEND")==0) return ParametersENDEnum;
               else if (strcmp(name,"XY")==0) return XYEnum;
@@ -868 +875 @@
               else if (strcmp(name,"Moulin")==0) return MoulinEnum;
               else if (strcmp(name,"Pengrid")==0) return PengridEnum;
-              else if (strcmp(name,"Penpair")==0) return PenpairEnum;
+         else stage=8;
+   }
+   if(stage==8){
+              if (strcmp(name,"Penpair")==0) return PenpairEnum;
               else if (strcmp(name,"Profiler")==0) return ProfilerEnum;
               else if (strcmp(name,"MatrixParam")==0) return MatrixParamEnum;
@@ -875 +885 @@
               else if (strcmp(name,"NodeSId")==0) return NodeSIdEnum;
               else if (strcmp(name,"ElementSId")==0) return ElementSIdEnum;
-         else stage=8;
-   }
-   if(stage==8){
-              if (strcmp(name,"VectorParam")==0) return VectorParamEnum;
+              else if (strcmp(name,"VectorParam")==0) return VectorParamEnum;
               else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum;
               else if (strcmp(name,"Segment")==0) return SegmentEnum;
@@ -991 +998 @@
               else if (strcmp(name,"TotalGroundedBmb")==0) return TotalGroundedBmbEnum;
               else if (strcmp(name,"TotalSmb")==0) return TotalSmbEnum;
-              else if (strcmp(name,"P0")==0) return P0Enum;
+         else stage=9;
+   }
+   if(stage==9){
+              if (strcmp(name,"P0")==0) return P0Enum;
               else if (strcmp(name,"P0Array")==0) return P0ArrayEnum;
               else if (strcmp(name,"P1")==0) return P1Enum;
@@ -998 +1008 @@
               else if (strcmp(name,"P1bubblecondensed")==0) return P1bubblecondensedEnum;
               else if (strcmp(name,"P2")==0) return P2Enum;
-         else stage=9;
-   }
-   if(stage==9){
-              if (strcmp(name,"P2bubble")==0) return P2bubbleEnum;
+              else if (strcmp(name,"P2bubble")==0) return P2bubbleEnum;
               else if (strcmp(name,"P2bubblecondensed")==0) return P2bubblecondensedEnum;
               else if (strcmp(name,"P2xP1")==0) return P2xP1Enum;

issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp

@@ -128 +128 @@
         }
 }/*}}}*/
+int IoCodeToEnumAmr(int enum_in){/*{{{*/
+        switch(enum_in){
+                case 1: return AmrBamgEnum;
+                case 2: return AmrNeopzEnum; 
+                default: _error_("Marshalled AMR code \""<<enum_in<<"\" not supported yet"); 
+        }
+}/*}}}*/
 
 int IoCodeToEnumVertexEquation(int enum_in){/*{{{*/

issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.h

@@ -10 +10 @@
 int IoCodeToEnumMaterials(int enum_in);
 int IoCodeToEnumTimestepping(int enum_in);
+int IoCodeToEnumAmr(int enum_in);
 
 int IoCodeToEnumVertexEquation(int enum_in);

issm/trunk-jpl/src/m/classes/amr.m

@@ -50 +50 @@
                 function marshall(self,prefix,md,fid) % {{{
 
-                        scale = md.constants.yts;
+                        WriteData(fid,prefix,'name','md.amr.type','data',2,'format','Integer');
                         WriteData(fid,prefix,'object',self,'fieldname','level_max','format','Integer');
                         WriteData(fid,prefix,'object',self,'fieldname','region_level_1','format','Double');

issm/trunk-jpl/src/m/classes/trans.js

@@ -9 +9 @@
 
                 //full analysis: Stressbalance, Masstransport and Thermal but no groundingline migration for now
-                this.issmb           = 1;
-                this.ismasstransport = 1;
-                this.isstressbalance = 1;
-                this.isthermal       = 1;
-                this.isgroundingline = 0;
-                this.isgia           = 0;
+                this.issmb             = 1;
+                this.ismasstransport   = 1;
+                this.isstressbalance   = 1;
+                this.isthermal         = 1;
+                this.isgroundingline   = 0;
+                this.isgia             = 0;
                 this.isdamageevolution = 0;
-                this.ismovingfront   = 0;
-                this.ishydrology     = 0;
-                this.isslr           = 0;
-                this.iscoupler       = 0;
-                this.amr_frequency   = 1;
+                this.ismovingfront     = 0;
+                this.ishydrology       = 0;
+                this.isslr             = 0;
+                this.iscoupler         = 0;
+                this.amr_frequency     = 0;
 
                 //default output

issm/trunk-jpl/src/m/classes/transient.m

@@ -68 +68 @@
                         self.isoceancoupling = 0;
                         self.iscoupler       = 0;
-                        self.amr_frequency      = 1;
+                        self.amr_frequency      = 0;
 
                         %default output

issm/trunk-jpl/src/m/classes/transient.py

@@ -84 +84 @@
                 
                 #full analysis: Stressbalance, Masstransport and Thermal but no groundingline migration for now
-                self.issmb = True
-                self.ismasstransport = True
-                self.isstressbalance = True
-                self.isthermal       = True
-                self.isgroundingline = False
-                self.isgia           = False
-                self.isesa           = False
+                self.issmb             = True
+                self.ismasstransport   = True
+                self.isstressbalance   = True
+                self.isthermal         = True
+                self.isgroundingline   = False
+                self.isgia             = False
+                self.isesa             = False
                 self.isdamageevolution = False
-                self.ismovingfront   = False
-                self.ishydrology     = False
-                self.isslr           = False
-                self.isoceancoupling = False
-                self.iscoupler       = False
-                self.amr_frequency      = 1
+                self.ismovingfront     = False
+                self.ishydrology       = False
+                self.isslr             = False
+                self.isoceancoupling   = False
+                self.iscoupler         = False
+                self.amr_frequency     = 0
 
                 #default output