Changeset 22241

Timestamp:

11/08/17 05:07:49 (7 years ago)

Author:

tsantos

Message:

CHG: Updated AMR methods for NeoPZ and Bamg, insert option to restart a simulation using error estimators and Bamg, deleted some non-used methods.

Location:

issm/trunk-jpl/src

Files:

• c/classes/AdaptiveMeshRefinement.cpp (modified) (15 diffs)
• c/classes/AdaptiveMeshRefinement.h (modified) (5 diffs)
• c/classes/AmrBamg.cpp (modified) (4 diffs)
• c/classes/AmrBamg.h (modified) (2 diffs)
• c/classes/FemModel.cpp (modified) (19 diffs)
• c/classes/FemModel.h (modified) (1 diff)
• c/cores/transient_core.cpp (modified) (2 diffs)
• c/modules/ModelProcessorx/CreateParameters.cpp (modified) (2 diffs)
• c/shared/Enum/EnumDefinitions.h (modified) (2 diffs)
• c/shared/Enum/EnumToStringx.cpp (modified) (2 diffs)
• m/classes/amr.js (modified) (5 diffs)
• m/classes/amr.m (modified) (5 diffs)
• m/classes/amr.py (modified) (5 diffs)
• m/contrib/tsantos/AMRexportVTK.m (modified) (2 diffs)
• m/contrib/tsantos/mismip/gl_position.m (modified) (1 diff)
• m/contrib/tsantos/mismip/ice_evolution.m (modified) (1 diff)
• m/contrib/tsantos/remesh.m (modified) (2 diffs)

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

@@ -30 +30 @@
         this->fathermesh                                                = new TPZGeoMesh(*cp.fathermesh);
         this->previousmesh                                      = new TPZGeoMesh(*cp.previousmesh);
+        this->refinement_type                           = cp.refinement_type;
         this->level_max                                         = cp.level_max;
         this->radius_level_max                          = cp.radius_level_max;
@@ -38 +39 @@
         this->thicknesserror_threshold  = cp.thicknesserror_threshold;
         this->deviatoricerror_threshold = cp.deviatoricerror_threshold;
+        this->max_deviatoricerror                       = cp.max_deviatoricerror;
+        this->max_thicknesserror                        = cp.max_thicknesserror;
         this->sid2index.clear();
         this->sid2index.resize(cp.sid2index.size());
@@ -61 +64 @@
         if(this->fathermesh)    delete this->fathermesh;
         if(this->previousmesh)  delete this->previousmesh;
+        this->refinement_type                           = -1;
         this->level_max                                         = -1;
         this->radius_level_max                          = -1;
@@ -69 +73 @@
         this->thicknesserror_threshold  = -1;
         this->deviatoricerror_threshold = -1;
+        this->max_deviatoricerror                       = -1;
+        this->max_thicknesserror                        = -1;
         this->sid2index.clear();
         this->index2sid.clear();
@@ -79 +85 @@
         this->fathermesh                                                = NULL;
         this->previousmesh                                      = NULL;
+        this->refinement_type                           = -1;
         this->level_max                                         = -1;
         this->radius_level_max                          = -1;
@@ -87 +94 @@
         this->thicknesserror_threshold  = -1;
         this->deviatoricerror_threshold = -1;
+        this->max_deviatoricerror                       = -1;
+        this->max_thicknesserror                        = -1;
         this->sid2index.clear();
         this->index2sid.clear();
@@ -94 +103 @@
 
 /*Mesh refinement methods*/
-void AdaptiveMeshRefinement::ExecuteRefinement(int numberofpoints,double* xylist,int* pnewnumberofvertices,int* pnewnumberofelements,double** px,double** py,int** pelementslist){/*{{{*/
-
+void AdaptiveMeshRefinement::ExecuteRefinement(double* gl_distance,double* if_distance,double* deviatoricerror,double* thicknesserror,int* pnewnumberofvertices,int* pnewnumberofelements,double** px,double** py,int** pelementslist){/*{{{*/
+        
         /*IMPORTANT! pelementslist are in Matlab indexing*/
         /*NEOPZ works only in C indexing*/
         if(!this->fathermesh || !this->previousmesh) _error_("Impossible to execute refinement: fathermesh or previousmesh is NULL!\n");
+        if(this->refinement_type!=0 && this->refinement_type!=1) _error_("Impossible to execute refinement: refinement type is not defined!\n");
+
+        /*Input verifications*/
+        if(this->deviatoricerror_threshold>0    && !deviatoricerror) _error_("deviatoricerror is NULL!\n");
+        if(this->thicknesserror_threshold>0             && !thicknesserror)     _error_("thicknesserror is NULL!\n");
+        if(this->groundingline_distance>0               && !gl_distance)                _error_("gl_distance is NULL!\n");
+        if(this->icefront_distance>0                            && !if_distance)                _error_("if_distance is NULL!\n");
 
         /*Intermediaries*/
         bool verbose=VerboseSolution();
-
-        /*Refine the mesh using level max*/
-        this->RefineMesh(verbose,this->previousmesh,numberofpoints,xylist);
-        
-        /*Get new geometric mesh in ISSM data structure*/
-        this->GetMesh(this->previousmesh,pnewnumberofvertices,pnewnumberofelements,px,py,pelementslist);
-
-        /*Verify the new geometry*/
-        this->CheckMesh(pnewnumberofvertices,pnewnumberofelements,px,py,pelementslist);
-        
-        if(verbose) _printf_("\trefinement process done!\n");
-}
-/*}}}*/
-void AdaptiveMeshRefinement::ExecuteRefinement(double* gl_elementdistance,double* if_elementdistance,double* deviatoricerror,double* thicknesserror,int* pnewnumberofvertices,int* pnewnumberofelements,double** px,double** py,int** pelementslist){/*{{{*/
-
-        /*IMPORTANT! pelementslist are in Matlab indexing*/
-        /*NEOPZ works only in C indexing*/
-        if(!this->fathermesh || !this->previousmesh) _error_("Impossible to execute refinement: fathermesh or previousmesh is NULL!\n");
-
-        /*Intermediaries*/
-        bool verbose=true;
         
         /*Execute refinement*/
-        this->RefinementProcess(verbose,gl_elementdistance,if_elementdistance,deviatoricerror,thicknesserror);
+        this->RefineMeshOneLevel(verbose,gl_distance,if_distance,deviatoricerror,thicknesserror);
    
         /*Get new geometric mesh in ISSM data structure*/
@@ -134 +129 @@
 }
 /*}}}*/
-void AdaptiveMeshRefinement::RefinementProcess(bool &verbose,double* gl_elementdistance,double* if_elementdistance,double* deviatoricerror,double* thicknesserror){/*{{{*/
-   
-        if(verbose) _printf_("\n\trefinement process started (level max = " << this->level_max << ")\n");
+void AdaptiveMeshRefinement::RefineMeshOneLevel(bool &verbose,double* gl_distance,double* if_distance,double* deviatoricerror,double* thicknesserror){/*{{{*/
         
         /*Intermediaries*/
-        TPZGeoMesh* nohangingnodesmesh=NULL;
-
-        double mean_mask                = 0;
-        double mean_tauerror = 0;
-        double mean_Herror      = 0;
-        double group_Herror     = 0;
-        int index,sid;
-        std::vector<int> specialfatherindex; specialfatherindex.clear();
-
-        /*Calculate mean values*/
-        /*
-        for(int i=0;i<this->sid2index.size();i++){
-                mean_mask               += masklevelset[i]; 
-                mean_tauerror   += deviatorictensorerror[i]; 
-                mean_Herror             += thicknesserror[i];
-        }
-        mean_mask               /= this->sid2index.size();
-        mean_tauerror   /= this->sid2index.size();
-        mean_Herror             /= this->sid2index.size();
-        */
-        if(verbose) _printf_("\t\tdeal with special elements...\n");
-        /*Deal with special elements*/
-        for(int i=0;i<this->specialelementsindex.size();i++){
-                if(this->specialelementsindex[i]==-1) continue;
-                /*Get special element and verify*/
-                TPZGeoEl* geoel=this->previousmesh->Element(this->specialelementsindex[i]);
-                if(!geoel)_error_("special element (sid) "<<i<<" is null!\n");
-                if(geoel->HasSubElement())_error_("special element (sid) "<<i<<" has "<<geoel->NSubElements()<<" subelements!\n");
-                if(geoel->MaterialId()!=this->GetElemMaterialID()) _error_("geoel->MaterialId is not GetElemMaterialID!\n");
-                if(!geoel->Father())_error_("father of special element (sid) "<<i<<" is null!\n");
-                
-                /*Get element's siblings and verify*/
-                TPZGeoEl* father=geoel->Father();
-                TPZVec<TPZGeoEl *> siblings;
-                father->GetHigherSubElements(siblings);
-                std::vector<int> sidvec; sidvec.resize(siblings.size());
-                for (int j=0;j<siblings.size();j++){
-                        if(!siblings[j]) _error_("special element (sid) "<<i<<" has a null siblings null!\n"); 
-                        sidvec[j]=this->index2sid[siblings[j]->Index()];
-                }
-                
-                /*Now, reset the data strucure and verify if the siblings should be deleted*/   
-                if(siblings.size()<4){
-                        /*Reset subelements in the father*/
-                        father->ResetSubElements();
-                }else{
-                        if(siblings.size()!=4) _error_("element (index) "<<father->Index()<<" has "<<father->NSubElements()<<" subelements!\n");
-                }
-                for (int j=0;j<siblings.size();j++){
-                        for(int k=0;k<this->specialelementsindex.size();k++){
-                                if(this->specialelementsindex[k]==siblings[j]->Index()){ 
-                                        index                                                                   = siblings[j]->Index();
-                                        if(index<0) _error_("index is null!\n");
-                                        sid                                                                     = this->index2sid[index];
-                                        if(sid<0) _error_("sid is null!\n");
-                                        this->specialelementsindex[k]   = -1;
-                                        this->index2sid[index]                  = -1;
-                                        this->sid2index[sid]                            = -1;
-                                }
-                        }
-                        if(siblings.size()<4){
-                                /*Ok, the special element can be deleted*/
-                                siblings[j]->ResetSubElements();
-                                this->previousmesh->DeleteElement(siblings[j],siblings[j]->Index());
-                        }
-                }
-                
-                /*Now, verify if the father should be refined with uniform pattern (smoother)*/
-                if(siblings.size()==3){//it keeps the mesh with uniform elements
-                        /*Father has uniform subelements now*/
-                        TPZVec<TPZGeoEl *> sons;
-                        father->Divide(sons);
-                }else{
-                        specialfatherindex.push_back(father->Index());
-                }
-                if(this->specialelementsindex[i]!=-1) _error_("special element "<<i<<" was not deleted!\n");    
-        }
-        this->previousmesh->BuildConnectivity();
-        
-        if(verbose) _printf_("\t\tuniform refinement...\n");
-        /*Deal with uniform elemnts*/
-        for(int i=0;i<this->sid2index.size();i++){
-                if(this->sid2index[i]==-1) continue;
-                /*Get element and verify*/
-                TPZGeoEl* geoel=this->previousmesh->Element(this->sid2index[i]);
-                if(geoel->HasSubElement()) _error_("element (sid) "<<i<<" has "<<geoel->NSubElements()<<" subelements!\n");
-                if(geoel->MaterialId()!=this->GetElemMaterialID()) _error_("geoel->MaterialId is not GetElemMaterialID!\n");
-
-                /*Refine process*/
-                if(thicknesserror[i]>mean_Herror)
-                {       
-                        int count=0;
-                        TPZGeoEl* father=geoel->Father();
-                        if(father){
-                                for(int j=3;j<6;j++){
-                                        index=father->Neighbour(j).Element()->Index();
-                                        for(int k=0;k<specialfatherindex.size();k++) if(specialfatherindex[k]==index) count++;
-                                }
-                        }
-                        TPZVec<TPZGeoEl *> sons;
-                        if(geoel->Level()<this->level_max && count==0) geoel->Divide(sons);
-                } 
-                else if(geoel->Level()>0)
-                {/*Unrefine process*/
-                        
-                        /*Get siblings and verify*/
-                        TPZVec<TPZGeoEl *> siblings;
-                        geoel->Father()->GetHigherSubElements(siblings);
-                        //if(siblings.size()<4) _error_("Impossible to refine: geoel (index) "<<this->sid2index[i]<<" has less than 3 siblings!\n");    
-                        if(siblings.size()>4) continue;//Father has more then 4 sons, this group should not be unrefined.
-                        
-                        /*Compute the error of the group*/
-                        group_Herror=0;
-                        for(int j=0;j<siblings.size();j++){
-                                index           = siblings[j]->Index();
-                                sid             = this->index2sid[index];
-                                if(sid==-1) continue;
-                                group_Herror+=thicknesserror[sid];
-                        }
-                        /*Verify if this group should be unrefined*/
-                        if(group_Herror>0 && group_Herror<0*mean_Herror){ //itapopo
-                                /*Reset subelements in the father*/
-                                this->previousmesh->Element(geoel->Father()->Index())->ResetSubElements();
-                                /*Delete the elements and set their indexes in the index2sid and sid2index*/
-                                for (int j=0;j<siblings.size();j++){
-                                        index   = siblings[j]->Index();
-                                        sid     = this->index2sid[index];
-                                        this->index2sid[index]=-1;
-                                        if(sid!=-1) this->sid2index[sid]=-1;
-                                        this->previousmesh->DeleteElement(siblings[j],siblings[j]->Index());
-                                }//for j
-                        }//if
-                }/*Unrefine process*/
-        }//for i
-        this->previousmesh->BuildConnectivity();
-        
-        if(verbose) _printf_("\t\trefine to avoid hanging nodes...\n");
-        this->RefineMeshToAvoidHangingNodes(verbose,this->previousmesh);
-        
-                //nohangingnodesmesh = this->CreateRefPatternMesh(newmesh); itapopo tentar otimizar
-
-        if(verbose) _printf_("\trefinement process done!\n");
-}
-/*}}}*/
-void AdaptiveMeshRefinement::RefineMesh(bool &verbose,TPZGeoMesh* gmesh,int numberofpoints,double* xylist){/*{{{*/
-        
-        /*Verify if there are points*/
-        if(numberofpoints==0) return;   
-
-        /*Intermediaries*/
-        int nelem                       =-1;
-        int side2D                      = 6;
-        double radius_h =-1;
-        double radius_hmax=std::max(this->radius_level_max,std::max(this->groundingline_distance,this->icefront_distance));
-        int count                       =-1;
-        double mindistance=0.;
-        double distance =0.;;
+        int nelem                                                       =-1;
+        int side2D                                                      = 6;
+        int sid                                                         =-1;
+        int count                                                       =-1;
+        int criteria                                            =-1;
+        int numberofcriteria                            =-1;
+        int nconformelements                            = this->sid2index.size();
+        double gl_radius_h                              =-1;
+        double gl_radius_hmax                   = std::max(this->radius_level_max,this->groundingline_distance);
+        double if_radius_h                              =-1;
+        double if_radius_hmax                   = std::max(this->radius_level_max,this->icefront_distance);
+        double gl_groupdistance                 =-1;
+        double if_groupdistance                 =-1;
+        double d_maxerror                                       =-1;
+        double t_maxerror                                       =-1;
+        double deviatoric_grouperror    =-1;
+        double thickness_grouperror     =-1;
+        TPZGeoMesh* gmesh                                       = NULL; 
         TPZVec<REAL> qsi(2,0.),cp(3,0.);
         TPZVec<TPZGeoEl*> sons;
- 
-        /*First, delete the special elements*/
-        this->DeleteSpecialElements(verbose,gmesh);
+        std::vector<int> index;
+
+        /*Calculate the number of criteria{{{*/
+        numberofcriteria=0;
+        if(this->deviatoricerror_threshold>0)   numberofcriteria++;
+        if(this->thicknesserror_threshold>0)    numberofcriteria++;
+        if(this->groundingline_distance>0)              numberofcriteria++;
+        if(this->icefront_distance>0)                           numberofcriteria++;
+        /*}}}*/
+
+        /*Calculate the maximum of the estimators, if requested{{{*/
+        if(this->deviatoricerror_threshold>0 && this->max_deviatoricerror<0){ 
+                for(int i=0;i<nconformelements;i++) this->max_deviatoricerror=max(this->max_deviatoricerror,deviatoricerror[i]);
+        }
+        if(this->thicknesserror_threshold>0 && this->max_thicknesserror<0){
+                for(int i=0;i<nconformelements;i++) this->max_thicknesserror=max(this->max_thicknesserror,thicknesserror[i]);
+        }
+        /*}}}*/
+
+        /*First, verify if special elements have min distance or high errors{{{*/
+        gmesh=this->previousmesh;
+        for(int i=0;i<this->specialelementsindex.size();i++){
+                if(this->specialelementsindex[i]==-1) _error_("index is -1!\n");
+                if(!gmesh->Element(this->specialelementsindex[i])) _error_("element is null!\n");
+                if(!gmesh->Element(this->specialelementsindex[i])->Father()) _error_("father is null!\n");
+                if(gmesh->Element(this->specialelementsindex[i])->HasSubElement()) _error_("special element has sub elements!\n");
+                sons.clear();
+                gmesh->Element(this->specialelementsindex[i])->Father()->GetHigherSubElements(sons);
+                /*Limits*/
+                gl_radius_h = gl_radius_hmax*std::pow(this->gradation,this->level_max-gmesh->Element(this->specialelementsindex[i])->Level());
+                if_radius_h = if_radius_hmax*std::pow(this->gradation,this->level_max-gmesh->Element(this->specialelementsindex[i])->Level());
+                d_maxerror      = this->deviatoricerror_threshold*this->max_deviatoricerror;
+                t_maxerror      = this->thicknesserror_threshold*this->max_thicknesserror;
+                /*Calculate the distance and error of the group (sons)*/
+                gl_groupdistance=INFINITY;if_groupdistance=INFINITY;deviatoric_grouperror=0;thickness_grouperror=0;
+                for(int s=0;s<sons.size();s++){
+                        sid=this->index2sid[sons[s]->Index()];
+                        if(sid<0) continue;
+                        if(this->groundingline_distance>0)              gl_groupdistance=std::min(gl_groupdistance,gl_distance[sid]); 
+                        if(this->icefront_distance>0)                           if_groupdistance=std::min(if_groupdistance,if_distance[sid]); 
+                        if(this->deviatoricerror_threshold>0)   deviatoric_grouperror+=deviatoricerror[sid];
+                        if(this->thicknesserror_threshold>0)    thickness_grouperror+=thicknesserror[sid];
+                }       
+                criteria=0;
+                if(this->groundingline_distance>0               && gl_groupdistance<gl_radius_h)                        criteria++;
+                if(this->icefront_distance>0                            && if_groupdistance<if_radius_h)                        criteria++;
+                if(this->deviatoricerror_threshold>0    && deviatoric_grouperror>d_maxerror)    criteria++;
+                if(this->thicknesserror_threshold>0             && thickness_grouperror>t_maxerror)             criteria++;
+                /*Finally, it keeps the father index if it must be refine*/
+                if(criteria) index.push_back(gmesh->Element(this->specialelementsindex[i])->FatherIndex());
+        }
+        /*}}}*/
+        
+        /*Now, detele the special elements{{{*/
+        if(this->refinement_type==1) this->DeleteSpecialElements(verbose,gmesh);
+        else this->specialelementsindex.clear();
+        /*}}}*/
+
+        /*Set the mesh and delete previousmesh if refinement type is 0{{{*/
+        if(this->refinement_type==0){
+                delete this->previousmesh;      
+                gmesh=this->fathermesh;
+        }
+        /*}}}*/
+        
+        /*Unrefinement process: loop over level of refinements{{{*/
+        if(verbose) _printf_("\tunrefinement process...\n");
+        if(verbose) _printf_("\ttotal: ");
+        count=0;
+         
+        nelem=gmesh->NElements();//must keep here
+        for(int i=0;i<nelem;i++){
+                if(!gmesh->Element(i)) continue;
+                if(gmesh->Element(i)->MaterialId()!=this->GetElemMaterialID()) continue;
+                if(gmesh->Element(i)->HasSubElement()) continue;
+                if(gmesh->Element(i)->Level()==0) continue;
+                if(!gmesh->Element(i)->Father()) _error_("father is NULL!\n");
+                /*Limits with lag*/
+                gl_radius_h = this->lag*gl_radius_hmax*std::pow(this->gradation,this->level_max-gmesh->Element(i)->Level());
+                if_radius_h = this->lag*if_radius_hmax*std::pow(this->gradation,this->level_max-gmesh->Element(i)->Level());
+                d_maxerror      = 0.05*this->max_deviatoricerror;//itapopo definir melhor
+                t_maxerror      = 0.05*this->max_thicknesserror;//itapopo definir melhor
+                /*Get the sons of the father (sibilings)*/      
+                sons.clear();
+                gmesh->Element(i)->Father()->GetHigherSubElements(sons);
+                if(sons.size()!=4) continue;//delete just group of 4 elements. This avoids big holes in the mesh
+                /*Find the minimal distance and the error of the group*/        
+                gl_groupdistance=INFINITY;if_groupdistance=INFINITY;deviatoric_grouperror=0;thickness_grouperror=0;
+                for(int s=0;s<sons.size();s++){
+                        sid=this->index2sid[sons[s]->Index()];
+                        /*Verify if this group have solutions*/
+                        if(sid<0){gl_groupdistance=INFINITY;if_groupdistance=INFINITY;deviatoric_grouperror=INFINITY;thickness_grouperror=INFINITY;continue;} 
+                        /*Distance and error of the group*/
+                        if(this->groundingline_distance>0)              gl_groupdistance=std::min(gl_groupdistance,gl_distance[sid]); 
+                        if(this->icefront_distance>0)                           if_groupdistance=std::min(if_groupdistance,if_distance[sid]); 
+                        if(this->deviatoricerror_threshold>0)   deviatoric_grouperror+=deviatoricerror[sid]; 
+                        if(this->thicknesserror_threshold>0)    thickness_grouperror+=thicknesserror[sid]; 
+                }
+                /*Verify the criteria*/
+                criteria=0;
+                if(this->groundingline_distance>0               && gl_groupdistance>gl_radius_h)                        criteria++;
+                if(this->icefront_distance>0                            && if_groupdistance>if_radius_h)                        criteria++;
+                if(this->deviatoricerror_threshold>0    && deviatoric_grouperror<d_maxerror)    criteria++;
+                if(this->thicknesserror_threshold>0             && thickness_grouperror<t_maxerror)             criteria++;
+                /*Now, if the group attends the criteria, unrefine it*/
+                if(criteria==numberofcriteria){ 
+                        gmesh->Element(i)->Father()->ResetSubElements(); count++;
+                        for(int s=0;s<sons.size();s++){this->index2sid[sons[s]->Index()]=-1;gmesh->DeleteElement(sons[s],sons[s]->Index());}
+                }
+        }
+        if(verbose) _printf_(""<<count<<"\n");
+        /*Adjust the connectivities before continue*/
+        gmesh->BuildConnectivity();
+        /*}}}*/
         
         /*Refinement process: loop over level of refinements{{{*/
         if(verbose) _printf_("\trefinement process (level max = "<<this->level_max<<")\n");
-        for(int h=1;h<=this->level_max;h++){
-                if(verbose) _printf_("\tlevel "<<h<<" (total: ");
+        if(verbose) _printf_("\ttotal: ");
+        count=0;
+        nelem=gmesh->NElements();//must keep here
+        for(int i=0;i<nelem;i++){
+                if(!gmesh->Element(i)) continue;
+                if(gmesh->Element(i)->MaterialId()!=this->GetElemMaterialID()) continue;
+                if(gmesh->Element(i)->HasSubElement()) continue;
+                if(gmesh->Element(i)->Level()==this->level_max) continue;
+                /*Verify if this element has solutions*/
+                sid=this->index2sid[gmesh->Element(i)->Index()];
+                if(sid<0) continue;
+                /*Filter: set the region/radius for level h*/
+                gl_radius_h     = gl_radius_hmax*std::pow(this->gradation,this->level_max-(gmesh->Element(i)->Level()+1));//+1, current element level is <level_max
+                if_radius_h     = if_radius_hmax*std::pow(this->gradation,this->level_max-(gmesh->Element(i)->Level()+1));//+1, current element level is <level_max
+                d_maxerror      = this->deviatoricerror_threshold*this->max_deviatoricerror;
+                t_maxerror      = this->thicknesserror_threshold*this->max_thicknesserror;
+                /*Verify distance and error of the element, if requested*/
+                criteria=0;
+                if(this->groundingline_distance>0               && gl_distance[sid]<gl_radius_h)                criteria++; 
+                if(this->icefront_distance>0                            && if_distance[sid]<if_radius_h)                criteria++; 
+                if(this->deviatoricerror_threshold>0    && deviatoricerror[sid]>d_maxerror)     criteria++; 
+                if(this->thicknesserror_threshold>0             && thicknesserror[sid]>t_maxerror)      criteria++; 
+                /*Now, if it attends any criterion, keep the element index to refine in next step*/
+                if(criteria)index.push_back(i);
+        }
+        /*Now, refine the elements*/
+        for(int i=0;i<index.size();i++){ 
+                if(!gmesh->Element(index[i])) DebugStop();
+                if(!gmesh->Element(index[i])->HasSubElement()){gmesh->Element(index[i])->Divide(sons);count++;}
+        }
+        if(verbose) _printf_(""<<count<<"\n");
+        /*Adjust the connectivities before continue*/
+        gmesh->BuildConnectivity();
+        /*}}}*/
+
+        /*Now, apply smoothing and insert special elements to avoid hanging nodes{{{*/
+        this->RefineMeshWithSmoothing(verbose,gmesh);
+        if(this->refinement_type==0) this->previousmesh=this->CreateRefPatternMesh(gmesh);//in this case, gmesh==this->fathermesh
+        gmesh=this->previousmesh;//previous mesh is always refined to avoid hanging nodes
+        this->RefineMeshToAvoidHangingNodes(verbose,gmesh);
+        /*}}}*/
+}
+/*}}}*/
+int AdaptiveMeshRefinement::VerifyRefinementType(TPZGeoEl* geoel){/*{{{*/
+
+        /*
+         * 0 : no refinement
+         * 1 : special refinement (to avoid hanging nodes)
+         * 2 : uniform refinment
+         * */
+        if(!geoel) _error_("geoel is NULL!\n");
+        
+        /*Output*/
+        int type=0;
+        int count=0;
+        
+        /*Intermediaries*/
+        TPZVec<TPZGeoEl*> sons;
+        
+        /*Loop over neighboors (sides 3, 4 and 5)*/
+        for(int j=3;j<6;j++){
+                sons.clear();
+                geoel->Neighbour(j).Element()->GetHigherSubElements(sons);
+                if(sons.size()) count++; //if neighbour was refined
+                if(sons.size()>4) count++; //if neighbour's level is > element level+1
+        }
+        
+        /*Verify and return*/
+        if(count>1) type=2;
+        else type=count;
+        
+        return type;
+}
+/*}}}*/
+void AdaptiveMeshRefinement::RefineMeshWithSmoothing(bool &verbose,TPZGeoMesh* gmesh){/*{{{*/
+
+        /*Intermediaries*/
+        int nelem               =-1;
+        int count               =-1;
+        int type                        =-1;
+        int typecount   =-1;
+
+        TPZVec<TPZGeoEl*> sons;
+
+        /*Refinement process: loop over level of refinements*/
+        if(verbose) _printf_("\tsmoothing process (level max = "<<this->level_max<<")\n");
+        if(verbose) _printf_("\ttotal: ");
+
+        count=1;
+        while(count>0){
                 count=0;
-
-                /*Filter: set the region/radius for level h*/
-                radius_h=radius_hmax*std::pow(this->gradation,this->level_max-h);
-
-                /*Find the minimal distance of the elements (center point) to the points */ 
                 nelem=gmesh->NElements();//must keep here
-                for(int i=0;i<nelem;i++){//itapopo pode-se reduzir o espaço de busca aqui
+                for(int i=0;i<nelem;i++){
                         if(!gmesh->Element(i)) continue;
                         if(gmesh->Element(i)->MaterialId()!=this->GetElemMaterialID()) continue;
                         if(gmesh->Element(i)->HasSubElement()) continue;
-                        if(gmesh->Element(i)->Level()>=h) continue;
-                        gmesh->Element(i)->CenterPoint(side2D,qsi);
-                        gmesh->Element(i)->X(qsi,cp);
-                        mindistance=radius_h;
-                        for (int j=0;j<numberofpoints;j++){
-                                distance                = std::sqrt( (xylist[2*j]-cp[0])*(xylist[2*j]-cp[0])+(xylist[2*j+1]-cp[1] )*(xylist[2*j+1]-cp[1]) );
-                                mindistance = std::min(mindistance,distance);//min distance to the point
+                        if(gmesh->Element(i)->Level()==this->level_max) continue;
+                        /*loop over neighboors (sides 3, 4 and 5). Important: neighbours has the same dimension of the element*/
+                        type=this->VerifyRefinementType(gmesh->Element(i));
+                        if(type<2){
+                                typecount=0;
+                                for(int j=3;j<6;j++){
+                                        if(gmesh->Element(i)->Neighbour(j).Element()->HasSubElement()) continue;
+                                        if(gmesh->Element(i)->Neighbour(j).Element()->Index()==i) typecount++;//neighbour==this element, element at the border
+                                        if(this->VerifyRefinementType(gmesh->Element(i)->Neighbour(j).Element())==1) typecount++;
+                                }
+                                if(typecount>1 && type==1) type=2;
+                                else if(typecount>2 && type==0) type=2;
                         }
-                        /*If the element is inside the region, refine it*/
-                        if(mindistance<radius_h){ 
-                                gmesh->Element(i)->Divide(sons);
-                                count++;
-                        }
-                }
-                if(verbose) _printf_(""<<count<<")\n");
-        }
-        /*Adjust the connectivities before continue*/
-        gmesh->BuildConnectivity();
-        /*}}}*/
-        
-        /*Unrefinement process: loop over level of refinements{{{*/
-        if(verbose) _printf_("\tunrefinement process...\n");
-        for(int h=this->level_max;h>=1;h--){
-                if(verbose) _printf_("\tlevel "<<h<<" (total: ");
-                count=0;
-                
-                /*Filter with lag: set the region/radius for level h*/
-                radius_h=this->lag*radius_hmax*std::pow(this->gradation,this->level_max-h);
-                
-                /*Find the minimal distance of the elements (center point) to the points */ 
-                nelem=gmesh->NElements();//must keep here
-                for(int i=0;i<nelem;i++){//itapopo pode-se reduzir o espaço de busca aqui
-                        if(!gmesh->Element(i)) continue;
-                        if(gmesh->Element(i)->MaterialId()!=this->GetElemMaterialID()) continue;
-                        if(gmesh->Element(i)->HasSubElement()) continue;
-                        if(gmesh->Element(i)->Level()!=h) continue;
-                        if(!gmesh->Element(i)->Father()) _error_("father is NULL!\n");
-                        /*Get the sons of the father (sibilings)*/      
-                        sons.clear();
-                        gmesh->Element(i)->Father()->GetHigherSubElements(sons);
-                        if(sons.size()!=4) continue;//delete just group of 4 elements. This avoids big holes in the mesh
-                        /*Find the minimal distance of the group*/      
-                        mindistance=INFINITY;
-                        for(int s=0;s<sons.size();s++){
-                                sons[s]->CenterPoint(side2D,qsi);
-                                sons[s]->X(qsi,cp);
-                                for (int j=0;j<numberofpoints;j++){
-                                        distance                = std::sqrt( (xylist[2*j]-cp[0])*(xylist[2*j]-cp[0])+(xylist[2*j+1]-cp[1] )*(xylist[2*j+1]-cp[1]) );
-                                        mindistance = std::min(mindistance,distance);//min distance to the point
-                                }
-                        }
-                        /*If the group is outside the region, unrefine the group*/
-                        if(mindistance>radius_h){ 
-                                gmesh->Element(i)->Father()->ResetSubElements();
-                                count++;
-                                for(int s=0;s<sons.size();s++){
-                                        gmesh->DeleteElement(sons[s],sons[s]->Index());
-                                }
-                        }
-                }
-                if(verbose) _printf_(""<<count<<")\n");
-        }
-        /*Adjust the connectivities before continue*/
-        gmesh->BuildConnectivity();
-        /*}}}*/
-        
-        /*Now, insert special elements to avoid hanging nodes*/
-        this->RefineMeshToAvoidHangingNodes(verbose,gmesh);
-}
-/*}}}*/
-void AdaptiveMeshRefinement::RefineMesh(TPZGeoMesh* gmesh,std::vector<int> &elements){/*{{{*/
-
-        /*Refine elements: uniform pattern refinement*/
-        for(int i=0;i<elements.size();i++){
-                /*Get geometric element and verify if it has already been refined*/
-                int index = elements[i];
-                TPZGeoEl * geoel = gmesh->Element(index);
-                if(geoel->HasSubElement()) _error_("Impossible to refine: geoel (index) " << index << " has subelements!\n");
-                if(geoel->MaterialId()!=this->GetElemMaterialID()) _error_("Impossible to refine: geoel->MaterialId is not GetElemMaterialID!\n");
-                /*Divide geoel*/
-                TPZVec<TPZGeoEl *> Sons;
-                geoel->Divide(Sons);
-        }
-        gmesh->BuildConnectivity();
+                        /*refine the element if requested*/
+                        if(type==2){gmesh->Element(i)->Divide(sons);    count++;}
+                }
+                if(verbose){
+                        if(count==0) _printf_(""<<count<<"\n");
+                        else _printf_(""<<count<<", ");
+                }
+                /*Adjust the connectivities before continue*/
+                gmesh->BuildConnectivity();
+        }
 }
 /*}}}*/
@@ -452 +439 @@
                 if(gmesh->Element(this->specialelementsindex[i])->Father()) gmesh->Element(this->specialelementsindex[i])->Father()->ResetSubElements();
                 gmesh->DeleteElement(gmesh->Element(this->specialelementsindex[i]),this->specialelementsindex[i]);
+      this->index2sid[this->specialelementsindex[i]]=-1;
                 count++;
         }
@@ -568 +556 @@
         if(nvertices<=0) _error_("Impossible to create initial mesh: nvertices is <= 0!\n");
    if(nelements<=0) _error_("Impossible to create initial mesh: nelements is <= 0!\n");
+        if(this->refinement_type!=0 && this->refinement_type!=1) _error_("Impossible to create initial mesh: refinement type is not defined!\n");
 
     /*Verify and creating initial mesh*/
@@ -591 +580 @@
    const int mat = this->GetElemMaterialID();
    TPZManVector<long> elem(this->GetNumberOfNodes(),0);
+        this->index2sid.clear(); this->index2sid.resize(nelements);
    this->sid2index.clear();
 
         for(int i=0;i<nelements;i++){
                 for(int j=0;j<this->GetNumberOfNodes();j++) elem[j]=elements[i*this->GetNumberOfNodes()+j]-1;//Convert Matlab to C indexing
-      /*reftype = 0: uniform, fast / reftype = 1: uniform and non-uniform (avoid hanging nodes), it is not too fast */
-      const int reftype = 1;
       switch(this->GetNumberOfNodes()){
-                        case 3: this->fathermesh->CreateGeoElement(ETriangle,elem,mat,index,reftype);   break;
+                        case 3: this->fathermesh->CreateGeoElement(ETriangle,elem,mat,index,this->refinement_type);     break;
          default:       _error_("mesh not supported yet");
                 }
       /*Define the element ID*/        
       this->fathermesh->ElementVec()[index]->SetId(i);
-                /*Initialize sid2index*/
+                /*Initialize sid2index and index2sid*/
                 this->sid2index.push_back((int)index);
+                this->index2sid[(int)index]=this->sid2index.size()-1;//keep the element sid
         }
    /*Build element and node connectivities*/
    this->fathermesh->BuildConnectivity();
         /*Set previous mesh*/
-        this->previousmesh=new TPZGeoMesh(*this->fathermesh);
+        if(this->refinement_type==1) this->previousmesh=new TPZGeoMesh(*this->fathermesh);
+        else this->previousmesh=this->CreateRefPatternMesh(this->fathermesh); 
 }
 /*}}}*/
@@ -622 +612 @@
    int reftype = 1;
    long index; 
-   
+
         //nodes
         newgmesh->NodeVec().Resize(nnodes);
@@ -630 +620 @@
    for(int i=0;i<nelem;i++){
         TPZGeoEl * geoel = gmesh->Element(i);
-      TPZManVector<long> elem(3,0);
+                
+                if(!geoel){
+                        index=newgmesh->ElementVec().AllocateNewElement();
+                        newgmesh->ElementVec()[index] = NULL;
+                        continue;
+                }
+      
+                TPZManVector<long> elem(3,0);
       for(int j=0;j<3;j++) elem[j] = geoel->NodeIndex(j);
      
       newgmesh->CreateGeoElement(ETriangle,elem,mat,index,reftype);
-      newgmesh->ElementVec()[index]->SetId(geoel->Id());
+                newgmesh->ElementVec()[index]->SetId(geoel->Id());
         
       TPZGeoElRefPattern<TPZGeoTriangle>* newgeoel = dynamic_cast<TPZGeoElRefPattern<TPZGeoTriangle>*>(newgmesh->ElementVec()[index]);
@@ -689 +686 @@
       }
    }
+
+        /*Now, build connectivities*/
         newgmesh->BuildConnectivity();
     
@@ -715 +714 @@
         //Verify if there are inf or NaN in coords
         for(int i=0;i<*nvertices;i++){
-                if(isnan((*px)[i]) || isinf((*px)[i])) _error_("Impossible to continue: px i=" << i <<" is NaN or Inf!\n"); 
-                if(isnan((*py)[i]) || isinf((*py)[i])) _error_("Impossible to continue: py i=" << i <<" is NaN or Inf!\n");
+                if(std::isnan((*px)[i]) || std::isinf((*px)[i])) _error_("Impossible to continue: px i=" << i <<" is NaN or Inf!\n"); 
+                if(std::isnan((*py)[i]) || std::isinf((*py)[i])) _error_("Impossible to continue: py i=" << i <<" is NaN or Inf!\n");
         }
         for(int i=0;i<*nelements;i++){
                 for(int j=0;j<this->GetNumberOfNodes();j++){
-                        if(isnan((*pelements)[i*GetNumberOfNodes()+j])) _error_("Impossible to continue: px i=" << i <<" is NaN!\n");
-                        if(isinf((*pelements)[i*GetNumberOfNodes()+j])) _error_("Impossible to continue: px i=" << i <<" is Inf!\n");
+                        if(std::isnan((*pelements)[i*GetNumberOfNodes()+j])) _error_("Impossible to continue: px i=" << i <<" is NaN!\n");
+                        if(std::isinf((*pelements)[i*GetNumberOfNodes()+j])) _error_("Impossible to continue: px i=" << i <<" is Inf!\n");
                 }
         }

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

@@ -54 +54 @@
          * radius_h = lag * initial_radius * gradation ^ (level_max-h)
          */
+        int refinement_type;                                            //0 uniform (faster); 1 refpattern  
         int level_max;                                                          //max level of refinement
         double radius_level_max;                                //initial radius which in the elements will be refined with level max
@@ -63 +64 @@
         double thicknesserror_threshold;                //if ==0, it will not be used
         double deviatoricerror_threshold;       //if ==0, it will not be used
+        double max_deviatoricerror;                     // Max value of the error estimator; in general, it is defined in the first time step. Attention with restart
+        double max_thicknesserror;                              // Max value of the error estimator; in general, it is defined in the first time step. Attention with restart
         /*}}}*/
         /*Public methods{{{*/
@@ -74 +77 @@
         void Initialize();
         void ExecuteRefinement(int numberofpoints,double* xylist,int* pnewnumberofvertices,int* pnewnumberofelements,double** px,double** py,int** pelementslist);
-        void ExecuteRefinement(double* gl_elementdistance,double* if_elementdistance,double* deviatoricerror,double* thicknesserror,int* pnewnumberofvertices,int* pnewnumberofelements,double** px,double** py,int** pelementslist);
+        void ExecuteRefinement(int numberofpoints,double* xylist,double* deviatoricerror,double* thicknesserror,int* pnewnumberofvertices,int* pnewnumberofelements,double** px,double** py,int** pelementslist);
+   void ExecuteRefinement(double* gl_distance,double* if_distance,double* deviatoricerror,double* thicknesserror,int* pnewnumberofvertices,int* pnewnumberofelements,double** px,double** py,int** pelementslist);
         void CreateInitialMesh(int &nvertices,int &nelements,double* x,double* y,int* elements);
         void CheckMesh(int* nvertices,int* nelements,double** px,double** py,int** pelements);
@@ -80 +84 @@
 private:
         /*Private attributes{{{*/
-        std::vector<int> sid2index;                                                                     // Vector that keeps index of PZGeoMesh elements used in the ISSM mesh (sid) 
-        std::vector<int> index2sid;                                                                     // Vector that keeps sid of issm mesh elements used in the neopz mesh (index) 
-        std::vector<int> specialelementsindex;                                          // Vector that keeps index of the special elements (created to avoid haning nodes) 
-        TPZGeoMesh *fathermesh;                                                                                 // Father Mesh is the entire mesh without refinement
-        TPZGeoMesh *previousmesh;                                                                               // Previous Mesh is the refined mesh
+        std::vector<int> sid2index;                                     // Vector that keeps index of PZGeoMesh elements used in the ISSM mesh (sid) 
+        std::vector<int> index2sid;                                     // Vector that keeps sid of issm mesh elements used in the neopz mesh (index) 
+        std::vector<int> specialelementsindex;          // Vector that keeps index of the special elements (created to avoid haning nodes) 
+        TPZGeoMesh *fathermesh;                                                 // Entire mesh without refinement if refinement_type==1; refined with hanging nodes if efinement_type==0
+        TPZGeoMesh *previousmesh;                                               // Refined mesh without hanging nodes (it is always refpattern type), used to generate ISSM mesh
         /*}}}*/
         /*Private methods{{{*/
-   void RefinementProcess(bool &verbose,double* partiallyfloatedelements,double* masklevelset,double* deviatorictensorerror,double* thicknesserror);
-        void RefineMesh(bool &verbose,TPZGeoMesh* gmesh,int numberofpoints,double* xylist);
-        void RefineMesh(TPZGeoMesh *gmesh,std::vector<int> &elements); 
-   void RefineMeshToAvoidHangingNodes(bool &verbose,TPZGeoMesh* gmesh);
+   void RefineMeshOneLevel(bool &verbose,double* gl_distance,double* if_distance,double* deviatoricerror,double* thicknesserror);
+        void RefineMeshWithSmoothing(bool &verbose,TPZGeoMesh* gmesh);
+        void RefineMeshToAvoidHangingNodes(bool &verbose,TPZGeoMesh* gmesh);
         void DeleteSpecialElements(bool &verbose,TPZGeoMesh* gmesh);
         void GetMesh(TPZGeoMesh* gmesh,int* nvertices,int* nelements,double** px,double** py,int** pelements);
@@ -98 +101 @@
         void PrintGMeshVTK(TPZGeoMesh *gmesh,std::ofstream &file,bool matColor=true);
         int GetVTK_ElType(TPZGeoEl* gel);
+        int VerifyRefinementType(TPZGeoEl* geoel);
         /*}}}*/
 };

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

@@ -17 +17 @@
 
 /*Constructor, copy, clean up and destructor*/
-AmrBamg::AmrBamg(IssmDouble hmin,IssmDouble hmax,int fieldenum_in,IssmDouble err_in,int keepmetric_in,IssmDouble gradation_in,
-                                                IssmDouble groundingline_resolution_in,IssmDouble groundingline_distance_in,
-                                                IssmDouble icefront_resolution_in,IssmDouble icefront_distance_in,
-                                                IssmDouble thicknesserror_resolution_in,IssmDouble thicknesserror_threshold_in,
-                                                IssmDouble deviatoricerror_resolution_in,IssmDouble deviatoricerror_threshold_in){/*{{{*/
+AmrBamg::AmrBamg(){/*{{{*/
 
-        this->fieldenum                                         = fieldenum_in;
-        this->keepmetric                                        = keepmetric_in;
-        this->groundingline_resolution  = groundingline_resolution_in;
-        this->groundingline_distance            = groundingline_distance_in;
-        this->icefront_resolution                       = icefront_resolution_in;
-        this->icefront_distance                         = icefront_distance_in;
-        this->thicknesserror_resolution         = thicknesserror_resolution_in;
-        this->thicknesserror_threshold  = thicknesserror_threshold_in;
-        this->deviatoricerror_resolution = deviatoricerror_resolution_in;
-        this->deviatoricerror_threshold  = deviatoricerror_threshold_in;
+        /*These attributes MUST be setup by FemModel*/
+        this->fieldenum                                         = -1;//fieldenum_in;
+        this->keepmetric                                        = -1;//keepmetric_in;
+        this->groundingline_resolution  = -1;//groundingline_resolution_in;
+        this->groundingline_distance            = -1;//groundingline_distance_in;
+        this->icefront_resolution                       = -1;//icefront_resolution_in;
+        this->icefront_distance                         = -1;//icefront_distance_in;
+        this->thicknesserror_resolution         = -1;//thicknesserror_resolution_in;
+        this->thicknesserror_threshold  = -1;//thicknesserror_threshold_in;
+        this->thicknesserror_maximum            = -1;
+        this->deviatoricerror_resolution = -1;//deviatoricerror_resolution_in;
+        this->deviatoricerror_threshold  = -1;//deviatoricerror_threshold_in;
+        this->deviatoricerror_maximum           = -1;
+        
+        /*Geometry and mesh as NULL*/
         this->geometry                                          = NULL;
         this->fathermesh                                        = NULL;
@@ -43 +44 @@
         this->options->coeff             = 1;
         this->options->errg              = 0.1;
-        this->options->gradation         = gradation_in;
+        this->options->gradation         = -1; //MUST be setup by the FemModel 
         this->options->Hessiantype       = 0;
         this->options->maxnbv            = 1e6;
@@ -54 +55 @@
         this->options->verbose           = 0;
         this->options->Crack             = 0;
-        this->options->KeepVertices      = 1; /*!!!!! VERY IMPORTANT !!!!!*/
+        this->options->KeepVertices      = 1; /*!!!!! VERY IMPORTANT !!!!! This avoid numerical errors when remeshing*/
         this->options->splitcorners      = 1;
-        this->options->hmin              = hmin;
-        this->options->hmax              = hmax;
-
-        this->options->err=xNew<IssmDouble>(1);
-        this->options->errSize[0]=1;
-        this->options->errSize[1]=1;
-        this->options->err[0] = err_in;
+        this->options->hmin              = -1;/*MUST be setup by the FemModel*/
+        this->options->hmax              = -1;/*MUST be setup by the FemModel*/
+        this->options->err                                      = xNew<IssmDouble>(1);
+        this->options->err[0]                           = -1;/*MUST be setup by the FemModel*/
+        this->options->errSize[0]                       = 1;
+        this->options->errSize[1]                       = 1;
 }
 /*}}}*/
@@ -174 +174 @@
         *pelementslist = elementslist;
 }/*}}}*/
+void AmrBamg::SetBamgOpts(IssmDouble hmin_in,IssmDouble hmax_in,IssmDouble err_in,IssmDouble gradation_in){/*{{{*/
+
+        if(!this->options) _error_("AmrBamg->options is NULL!");
+        
+        this->options->hmin     = hmin_in; 
+        this->options->hmax     = hmax_in; 
+        this->options->err[0]   = err_in; 
+        this->options->gradation= gradation_in; 
+}/*}}}*/

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

@@ -21 +21 @@
                 IssmDouble thicknesserror_resolution;
                 IssmDouble thicknesserror_threshold;
+                IssmDouble thicknesserror_maximum;
                 IssmDouble deviatoricerror_resolution;
                 IssmDouble deviatoricerror_threshold;
+                IssmDouble deviatoricerror_maximum;
 
                 /* Constructor, destructor etc*/
-                AmrBamg(IssmDouble hmin,IssmDouble hmax,int fieldenum_in,IssmDouble err_in,int keepmetric_in,IssmDouble gradation_in,
-                        IssmDouble groundingline_resolution_in,IssmDouble groundingline_distance_in,
-                  IssmDouble icefront_resolution_in,IssmDouble icefront_distance_in,
-                  IssmDouble thicknesserror_resolution_in,IssmDouble thicknesserror_threshold_in,
-                  IssmDouble deviatoricerror_resolution_in,IssmDouble deviatoricerror_threshold_in);
+                AmrBamg();
+                
                 ~AmrBamg();
 
@@ -35 +34 @@
                 void Initialize(int* elements,IssmDouble* x,IssmDouble* y,int numberofvertices,int numberofelements);
                 void ExecuteRefinementBamg(IssmDouble* field,IssmDouble* hmaxVertices,int* pnewnumberofvertices,int *pnewnumberofelements,IssmDouble** px,IssmDouble** py,IssmDouble** pz,int** pelementslist);
+                void SetBamgOpts(IssmDouble hmin_in,IssmDouble hmax_in,IssmDouble err_in,IssmDouble gradation_in);
 
         private:

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

@@ -2648 +2648 @@
                         YC_new[i]+=Ynew[vid];
                 }
-                XC_new[i]=XC_new[i]/3;
-                YC_new[i]=YC_new[i]/3;
+                XC_new[i]=XC_new[i]/3.;
+                YC_new[i]=YC_new[i]/3.;
         }
 
@@ -3029 +3029 @@
                 elementslist[elementswidth*i+0] = reCast<int>(id1[i])+1; //InterpMesh wants Matlab indexing
                 elementslist[elementswidth*i+1] = reCast<int>(id2[i])+1; //InterpMesh wants Matlab indexing
-                elementslist[elementswidth*i+2] = reCast<int>(id3[i])+1; //InterpMesh wants Matlab indexinf
+                elementslist[elementswidth*i+2] = reCast<int>(id3[i])+1; //InterpMesh wants Matlab indexing
         }
         
@@ -3608 +3608 @@
    Vector<IssmDouble>* vxc              = new Vector<IssmDouble>(numberofelements);
    Vector<IssmDouble>* vyc              = new Vector<IssmDouble>(numberofelements);
-        IssmDouble *x                                   = NULL;
-        IssmDouble *y                                   = NULL;
-        IssmDouble *z                                   = NULL; 
-
-        /*Get vertices coordinates*/
-        VertexCoordinatesx(&x,&y,&z,this->vertices,false) ;
+        IssmDouble* x                                   = NULL;
+        IssmDouble* y                                   = NULL;
+        IssmDouble* z                                   = NULL; 
+        IssmDouble* xyz_list                    = NULL;
+        IssmDouble x1,y1,x2,y2,x3,y3;
 
         /*Insert the element center coordinates*/
    for(int i=0;i<this->elements->Size();i++){
       Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
-      element->GetVerticesSidList(elem_vertices);
+      //element->GetVerticesSidList(elem_vertices);
       int sid = element->Sid();
-      vxc->SetValue(sid,(x[elem_vertices[0]]+x[elem_vertices[1]]+x[elem_vertices[2]])/3.,INS_VAL);
-      vyc->SetValue(sid,(y[elem_vertices[0]]+y[elem_vertices[1]]+y[elem_vertices[2]])/3.,INS_VAL);
+                element->GetVerticesCoordinates(&xyz_list); 
+                x1 = xyz_list[3*0+0];y1 = xyz_list[3*0+1];
+                x2 = xyz_list[3*1+0];y2 = xyz_list[3*1+1];
+                x3 = xyz_list[3*2+0];y3 = xyz_list[3*2+1];
+                vxc->SetValue(sid,(x1+x2+x3)/3.,INS_VAL);
+      vyc->SetValue(sid,(y1+y2+y3)/3.,INS_VAL);
    }
 
@@ -3636 +3639 @@
         xDelete<IssmDouble>(y);
         xDelete<IssmDouble>(z);
+        xDelete<IssmDouble>(xyz_list);
    xDelete<int>(elem_vertices);
    delete vxc;
@@ -3658 +3662 @@
         int* elem_vertices                                      = xNew<int>(elementswidth);
    IssmDouble* levelset                                         = xNew<IssmDouble>(elementswidth);
-   IssmDouble* x                                                                        = NULL;
-   IssmDouble* y                                                                        = NULL;
-   IssmDouble* z                                                                        = NULL;
+   IssmDouble* xyz_list                                                 = NULL;
         Vector<IssmDouble>* vx_zerolevelset             = new Vector<IssmDouble>(numberofelements);
         Vector<IssmDouble>* vy_zerolevelset             = new Vector<IssmDouble>(numberofelements);
@@ -3666 +3668 @@
         IssmDouble* y_zerolevelset                                      = NULL;
         int count,sid;
-        IssmDouble xcenter,ycenter;
-        
-        /*Get vertices coordinates*/
-        VertexCoordinatesx(&x,&y,&z,this->vertices,false) ;
+        IssmDouble xc,yc,x1,y1,x2,y2,x3,y3;
         
         /*Use the element center coordinate if level set is zero (grounding line or ice front), otherwise set NAN*/
@@ -3676 +3675 @@
       element->GetInputListOnVertices(levelset,levelset_type);
                 element->GetVerticesSidList(elem_vertices);
-                sid                     = element->Sid();
-                xcenter         = NAN;
-                ycenter         = NAN;  
+                sid= element->Sid();
+                element->GetVerticesCoordinates(&xyz_list); 
+                x1 = xyz_list[3*0+0];y1 = xyz_list[3*0+1];
+                x2 = xyz_list[3*1+0];y2 = xyz_list[3*1+1];
+                x3 = xyz_list[3*2+0];y3 = xyz_list[3*2+1];
+                xc      = NAN;
+                yc      = NAN;  
         Tria* tria      = xDynamicCast<Tria*>(element);
                 if(tria->IsIceInElement()){/*verify if there is ice in the element*/
                         if(levelset[0]*levelset[1]<0. || levelset[0]*levelset[2]<0. ||  
                                 abs(levelset[0]*levelset[1])<DBL_EPSILON || abs(levelset[0]*levelset[2])<DBL_EPSILON) {
-                                xcenter=(x[elem_vertices[0]]+x[elem_vertices[1]]+x[elem_vertices[2]])/3.;
-                                ycenter=(y[elem_vertices[0]]+y[elem_vertices[1]]+y[elem_vertices[2]])/3.;
+                                xc=(x1+x2+x3)/3.;
+                                yc=(y1+y2+y3)/3.;
                         }
                 }
-                vx_zerolevelset->SetValue(sid,xcenter,INS_VAL);
-                vy_zerolevelset->SetValue(sid,ycenter,INS_VAL);
+                vx_zerolevelset->SetValue(sid,xc,INS_VAL);
+                vy_zerolevelset->SetValue(sid,yc,INS_VAL);
+                xDelete<IssmDouble>(xyz_list);
         }
    /*Assemble and serialize*/
@@ -3720 +3724 @@
         xDelete<IssmDouble>(x_zerolevelset);
         xDelete<IssmDouble>(y_zerolevelset);
-   xDelete<IssmDouble>(x);
-   xDelete<IssmDouble>(y);
-   xDelete<IssmDouble>(z);
+   xDelete<IssmDouble>(xyz_list);
         delete vx_zerolevelset;
         delete vy_zerolevelset;
@@ -4474 +4476 @@
                 if(this->amrbamg->deviatoricerror_threshold>0)  this->GethmaxVerticesFromEstimators(hmaxvertices_serial,DeviatoricStressErrorEstimatorEnum);
         }
-
+        
         if(my_rank==0){
                 this->amrbamg->ExecuteRefinementBamg(vector_serial,hmaxvertices_serial,&newnumberofvertices,&newnumberofelements,&newx,&newy,&newz,&newelementslist);
@@ -4518 +4520 @@
         int* elements             = NULL;
         IssmDouble hmin,hmax,err,gradation;
-        IssmDouble groundingline_resolution,groundingline_distance,icefront_resolution,icefront_distance;
-        IssmDouble thicknesserror_resolution,thicknesserror_threshold,deviatoricerror_resolution,deviatoricerror_threshold;
-        int        fieldenum,keepmetric;
 
    /*Get rank*/
@@ -4531 +4530 @@
         this->GetMesh(this->vertices,this->elements,&x,&y,&z,&elements);
 
+        /*Create bamg data structures for bamg*/
+        this->amrbamg = new AmrBamg();
+        
         /*Get amr parameters*/
         this->parameters->FindParam(&hmin,AmrHminEnum);
         this->parameters->FindParam(&hmax,AmrHmaxEnum);
-        this->parameters->FindParam(&fieldenum,AmrFieldEnum);
         this->parameters->FindParam(&err,AmrErrEnum);
-        this->parameters->FindParam(&keepmetric,AmrKeepMetricEnum);
         this->parameters->FindParam(&gradation,AmrGradationEnum);
-        this->parameters->FindParam(&groundingline_resolution,AmrGroundingLineResolutionEnum);
-        this->parameters->FindParam(&groundingline_distance,AmrGroundingLineDistanceEnum);
-        this->parameters->FindParam(&icefront_resolution,AmrIceFrontResolutionEnum);
-        this->parameters->FindParam(&icefront_distance,AmrIceFrontDistanceEnum);
-        this->parameters->FindParam(&thicknesserror_resolution,AmrThicknessErrorResolutionEnum);
-        this->parameters->FindParam(&thicknesserror_threshold,AmrThicknessErrorThresholdEnum);
-        this->parameters->FindParam(&deviatoricerror_resolution,AmrDeviatoricErrorResolutionEnum);
-        this->parameters->FindParam(&deviatoricerror_threshold,AmrDeviatoricErrorThresholdEnum);
-
-        /*Create bamg data structures for bamg*/
-        this->amrbamg = new AmrBamg(hmin,hmax,fieldenum,err,keepmetric,gradation,
-                                                                                 groundingline_resolution,groundingline_distance,
-                                                                                 icefront_resolution,icefront_distance,
-                                                                                 thicknesserror_resolution,thicknesserror_threshold,
-                                                                                 deviatoricerror_resolution,deviatoricerror_threshold);
+        this->parameters->FindParam(&this->amrbamg->fieldenum,AmrFieldEnum);
+        this->parameters->FindParam(&this->amrbamg->keepmetric,AmrKeepMetricEnum);
+        this->parameters->FindParam(&this->amrbamg->groundingline_resolution,AmrGroundingLineResolutionEnum);
+        this->parameters->FindParam(&this->amrbamg->groundingline_distance,AmrGroundingLineDistanceEnum);
+        this->parameters->FindParam(&this->amrbamg->icefront_resolution,AmrIceFrontResolutionEnum);
+        this->parameters->FindParam(&this->amrbamg->icefront_distance,AmrIceFrontDistanceEnum);
+        this->parameters->FindParam(&this->amrbamg->thicknesserror_resolution,AmrThicknessErrorResolutionEnum);
+        this->parameters->FindParam(&this->amrbamg->thicknesserror_threshold,AmrThicknessErrorThresholdEnum);
+        this->parameters->FindParam(&this->amrbamg->thicknesserror_maximum,AmrThicknessErrorMaximumEnum);
+        this->parameters->FindParam(&this->amrbamg->deviatoricerror_resolution,AmrDeviatoricErrorResolutionEnum);
+        this->parameters->FindParam(&this->amrbamg->deviatoricerror_threshold,AmrDeviatoricErrorThresholdEnum);
+        this->parameters->FindParam(&this->amrbamg->deviatoricerror_maximum,AmrDeviatoricErrorMaximumEnum);
+        /*Set BamgOpts*/
+        this->amrbamg->SetBamgOpts(hmin,hmax,err,gradation);
 
         /*Re-create original mesh and put it in bamg structure (only cpu 0)*/
@@ -4608 +4607 @@
 
         /*Intermediaries*/
-        int numberofelements                            = this->elements->NumberOfElements();
-        IssmDouble* error_elements              = NULL;
-        IssmDouble *x                                           = NULL;
-        IssmDouble *y                                           = NULL;
-        IssmDouble *z                                           = NULL;
-        int *index                                                      = NULL;
-        IssmDouble maxerror,threshold,resolution;
-        int vid;
-        
+        int elementswidth                                                       = this->GetElementsWidth();
+        int numberofelements                                            = this->elements->NumberOfElements();
+        int numberofvertices                                            = this->vertices->NumberOfVertices();
+        IssmDouble* maxlength                                   = xNew<IssmDouble>(numberofelements);
+        IssmDouble* error_vertices                              = xNewZeroInit<IssmDouble>(numberofvertices);   
+        IssmDouble* error_elements                              = NULL;
+        IssmDouble* x                                                           = NULL;
+        IssmDouble* y                                                           = NULL;
+        IssmDouble* z                                                           = NULL;
+        int* index                                                                      = NULL;
+        IssmDouble maxerror,threshold,resolution,length;
+        IssmDouble L1,L2,L3;
+        int vid,v1,v2,v3;
+        bool refine;
+
+        /*Fill variables*/
         switch(errorestimator_type){
                 case ThicknessErrorEstimatorEnum: 
                         threshold       = this->amrbamg->thicknesserror_threshold;
                         resolution      = this->amrbamg->thicknesserror_resolution;
-                        this->ThicknessZZErrorEstimator(&error_elements);
+                        maxerror                = this->amrbamg->thicknesserror_maximum;
+                        this->ThicknessZZErrorEstimator(&error_elements);//error is serial, but the calculation is parallel
                         break;
                 case DeviatoricStressErrorEstimatorEnum:
                         threshold       = this->amrbamg->deviatoricerror_threshold;
                         resolution      = this->amrbamg->deviatoricerror_resolution;
-                        this->ZZErrorEstimator(&error_elements);
+                        maxerror                = this->amrbamg->deviatoricerror_maximum;
+                        this->ZZErrorEstimator(&error_elements);//error is serial, but the calculation is parallel
                         break;
                 default: _error_("not implemented yet");
         }
-
         if(!error_elements) _error_("error_elements is NULL!\n");
+
+        /*Find the max of the estimators if it was not provided*/
+        if(maxerror<DBL_EPSILON){
+                for(int i=0;i<numberofelements;i++) maxerror=max(maxerror,error_elements[i]);
+                switch(errorestimator_type){
+        case ThicknessErrorEstimatorEnum:                       this->amrbamg->thicknesserror_maximum   = maxerror;break;
+        case DeviatoricStressErrorEstimatorEnum:        this->amrbamg->deviatoricerror_maximum = maxerror;break;
+        }       
+        }
 
         /*Get mesh*/
         this->GetMesh(this->vertices,this->elements,&x,&y,&z,&index);
-
-        /*Find the max of the estimators (use error_elements)*/
-        maxerror=error_elements[0];
-        for(int i=0;i<numberofelements;i++) maxerror=max(maxerror,error_elements[i]);
-        
-        /*Fill hmaxvertices*/
+        
+        /*Fill error_vertices (this is the sum of all elements connected to the vertex)*/
         for(int i=0;i<numberofelements;i++){
-                if(error_elements[i]>threshold*maxerror){
-                        /*ok, fill the hmaxvertices using the element vertices*/
-                        for(int j=0;j<this->GetElementsWidth();j++){
-                                vid=index[i*this->GetElementsWidth()+j]-1;//Matlab to C indexing
+                v1=index[i*elementswidth+0]-1;//Matlab to C indexing
+                v2=index[i*elementswidth+1]-1;//Matlab to C indexing
+                v3=index[i*elementswidth+2]-1;//Matlab to C indexing
+                L1=sqrt(pow(x[v2]-x[v1],2)+pow(y[v2]-y[v1],2));
+                L2=sqrt(pow(x[v3]-x[v2],2)+pow(y[v3]-y[v2],2));
+                L3=sqrt(pow(x[v1]-x[v3],2)+pow(y[v1]-y[v3],2));
+                /*Fill the vectors*/
+                maxlength[i]            =max(L1,max(L2,L3));
+                error_vertices[v1]+=error_elements[i];
+                error_vertices[v2]+=error_elements[i];
+                error_vertices[v3]+=error_elements[i];
+        }       
+
+        /*Fill hmaxvertices with the criteria*/
+        for(int i=0;i<numberofelements;i++){
+                refine=false;
+                /*Refine any element if its error > phi*maxerror*/
+                if(error_elements[i]>threshold*maxerror) refine=true;
+                /*If the element size is closer to the resolution, verify the sum of error in the vertices*/
+                if(resolution/maxlength[i]>0.85){
+                        for(int j=0;j<elementswidth;j++){
+                                vid=index[i*elementswidth+j]-1;//Matlab to C indexing
+                                if(error_vertices[vid]>0.005*maxerror) refine=true;//itapopo must be better defined
+                        }
+                }
+                /*Now, fill the hmaxvertices if requested*/       
+                if(refine){
+                        for(int j=0;j<elementswidth;j++){
+                                vid=index[i*elementswidth+j]-1;//Matlab to C indexing
                                 if(xIsNan<IssmDouble>(hmaxvertices[vid])) hmaxvertices[vid]=resolution;
                                 else hmaxvertices[vid]=min(resolution,hmaxvertices[vid]);
@@ -4658 +4695 @@
         xDelete<int>(index);
    xDelete<IssmDouble>(error_elements);
+   xDelete<IssmDouble>(error_vertices);
+   xDelete<IssmDouble>(maxlength);
 }
 /*}}}*/
@@ -4714 +4753 @@
    int my_rank                                          = IssmComm::GetRank();
    int numberofelements                 = this->elements->NumberOfElements();
-        IssmDouble* element_label       = xNewZeroInit<IssmDouble>(numberofelements);
-        int numberofpoints                      = -1;
-        IssmDouble* xylist                      = NULL;
+        IssmDouble* gl_distance         = NULL;
+        IssmDouble* if_distance         = NULL;
+        IssmDouble* deviatoricerror= NULL;
+        IssmDouble* thicknesserror      = NULL;
         IssmDouble* newx                                = NULL;
    IssmDouble* newy                             = NULL;
@@ -4724 +4764 @@
         int newnumberofelements         = -1;
 
-        /*Get element_label, if requested*/
-        if(this->amr->groundingline_distance>0)         this->GetElementLabelFromZeroLevelSet(element_label,MaskGroundediceLevelsetEnum);
-   if(this->amr->icefront_distance>0)                           this->GetElementLabelFromZeroLevelSet(element_label,MaskIceLevelsetEnum);
-   if(this->amr->thicknesserror_threshold>0)            this->GetElementLabelFromEstimators(element_label,ThicknessErrorEstimatorEnum);
-   if(this->amr->deviatoricerror_threshold>0)   this->GetElementLabelFromEstimators(element_label,DeviatoricStressErrorEstimatorEnum);
-        this->GetPointsFromElementLabel(element_label,&numberofpoints,&xylist);
-
+        /*Get fields, if requested*/
+        if(this->amr->groundingline_distance>0)         this->GetElementDistanceToZeroLevelSet(&gl_distance,MaskGroundediceLevelsetEnum);
+   if(this->amr->icefront_distance>0)                           this->GetElementDistanceToZeroLevelSet(&if_distance,MaskIceLevelsetEnum);
+   if(this->amr->thicknesserror_threshold>0)            this->ThicknessZZErrorEstimator(&thicknesserror);       
+        if(this->amr->deviatoricerror_threshold>0)      this->ZZErrorEstimator(&deviatoricerror);       
+        
         if(my_rank==0){
-                this->amr->ExecuteRefinement(numberofpoints,xylist,&newnumberofvertices,&newnumberofelements,&newx,&newy,&newelementslist);
-      newz=xNewZeroInit<IssmDouble>(newnumberofvertices);
+                this->amr->ExecuteRefinement(gl_distance,if_distance,deviatoricerror,thicknesserror,
+                                                                                                &newnumberofvertices,&newnumberofelements,&newx,&newy,&newelementslist); 
+                newz=xNewZeroInit<IssmDouble>(newnumberofvertices);
                 if(newnumberofvertices<=0 || newnumberofelements<=0) _error_("Error in the ReMeshNeopz.");
         }
@@ -4760 +4800 @@
 
         /*Cleanup*/
-        xDelete<IssmDouble>(element_label);
-        xDelete<IssmDouble>(xylist);
+        xDelete<IssmDouble>(deviatoricerror);
+        xDelete<IssmDouble>(thicknesserror);
+        xDelete<IssmDouble>(gl_distance);
+        xDelete<IssmDouble>(if_distance);
 }
 /*}}}*/
@@ -4805 +4847 @@
         this->SetRefPatterns();
         this->amr = new AdaptiveMeshRefinement();
+        this->amr->refinement_type                                      = 1;//1 is refpattern; 0 is uniform (faster) 
         this->amr->level_max                                                    = level_max;
         this->amr->radius_level_max                             = radius_level_max;
@@ -4824 +4867 @@
 }
 /*}}}*/
-void FemModel::GetPointsFromElementLabel(IssmDouble* element_label,int* pnumberofpoints,IssmDouble** pxylist){/*{{{*/
-
-        if(!element_label) _error_("element_label is NULL!\n");
-
-        /*Outputs*/
-        int numberofpoints      = -1;
-        IssmDouble* xylist      = NULL;
-
-        /*Intermediaries*/
-   int numberofelements = this->elements->NumberOfElements();
-        int count                               = -1;
-   IssmDouble* xc                       = NULL;
-   IssmDouble* yc                       = NULL;
-
-        /*First, find the number of labeled elements (points)*/
-        count=0;
-        for(int i=0;i<numberofelements;i++){ 
-                if(element_label[i]>DBL_EPSILON) count++;
-        }
-         
-        /*Set number of points*/
-        numberofpoints=count;
-        if(count>0) xylist=xNew<IssmDouble>(2*numberofpoints);
-        
-        /*Get element center coordinates*/
-        this->GetElementCenterCoordinates(&xc,&yc);
-
-        /*Now, fill xylist data*/
-        count=0;
-        for(int i=0;i<numberofelements;i++){
-                if(element_label[i]>DBL_EPSILON){
-                        xylist[2*count] = xc[i];                        
-                        xylist[2*count+1]       = yc[i];
-                        count++;
-                }
-        }
-
-        /*Assign pointers*/
-        (*pxylist)=xylist;
-        *pnumberofpoints=numberofpoints;
-
-        /*Cleanup*/
-        xDelete<IssmDouble>(xc);
-        xDelete<IssmDouble>(yc);
-}
-/*}}}*/
-void FemModel::GetElementLabelFromZeroLevelSet(IssmDouble* element_label,int levelset_type){/*{{{*/
-
-        /*Here, "zero level set" means grounding line or ice front, depending on the level set type*/
-        /*element_label is 1 if the element zero level set, NAN otherwise*/
-        if(levelset_type!=MaskGroundediceLevelsetEnum && levelset_type!=MaskIceLevelsetEnum) _error_("level set type not implemented yet!");
-        if(!element_label) _error_("element_label is NULL!\n");
-        
-        /*Intermediaries*/
-        int elementswidth                       = this->GetElementsWidth();
-   int numberofelements                 = this->elements->NumberOfElements();
-        int* elem_vertices                                      = xNew<int>(elementswidth);
-   IssmDouble* levelset                                         = xNew<IssmDouble>(elementswidth);
-        Vector<IssmDouble>* velement_label              = new Vector<IssmDouble>(numberofelements);
-        IssmDouble* element_label_serial                        = NULL;
-        int sid                                                                                 = -1;
-        IssmDouble label                                                                = -1.;
-
-        /*Use the element center coordinate if level set is zero (grounding line or ice front), otherwise set NAN*/
-   for(int i=0;i<this->elements->Size();i++){
-      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
-      element->GetInputListOnVertices(levelset,levelset_type);
-                element->GetVerticesSidList(elem_vertices);
-                sid     = element->Sid();
-                label = NAN;    
-        Tria* tria      = xDynamicCast<Tria*>(element);
-                if(tria->IsIceInElement()){/*verify if there is ice in the element*/
-                        if(levelset[0]*levelset[1]<0. || levelset[0]*levelset[2]<0. ||  
-                                abs(levelset[0]*levelset[1])<DBL_EPSILON || abs(levelset[0]*levelset[2])<DBL_EPSILON) {
-                                label=1.;
-                        }
-                }
-                velement_label->SetValue(sid,label,INS_VAL);
-        }
-   
-        /*Assemble and serialize*/
-   velement_label->Assemble();
-   element_label_serial=velement_label->ToMPISerial();
-
-        /*Merge with the output*/
-   for(int i=0;i<numberofelements;i++){
-                if(!xIsNan<IssmDouble>(element_label_serial[i])) element_label[i]=element_label_serial[i];
-                else; //do nothing
-        }
-
-        /*Cleanup*/
-        xDelete<int>(elem_vertices);
-   xDelete<IssmDouble>(levelset);
-   xDelete<IssmDouble>(element_label_serial);
-        delete velement_label;
-}
-/*}}}*/
-void FemModel::GetElementLabelFromEstimators(IssmDouble* element_label,int estimator_type){/*{{{*/
-
-        /*element_label is 1 if the element zero level set, NAN otherwise*/
-        if(!element_label) _error_("element_label is NULL!\n");
-        
-        /*Intermediaries*/
-   int numberofelements                 = this->elements->NumberOfElements();
-   IssmDouble* elementerror     = NULL;
-        IssmDouble threshold                    = -1.;
-        IssmDouble maxerror                     = -1.;
-
-        switch(estimator_type){
-                case ThicknessErrorEstimatorEnum: 
-                        threshold=this->amr->thicknesserror_threshold;
-                        this->ThicknessZZErrorEstimator(&elementerror);
-                        break;
-                case DeviatoricStressErrorEstimatorEnum:
-                        threshold=this->amr->deviatoricerror_threshold;
-                        this->ZZErrorEstimator(&elementerror);
-                        break;
-                default: _error_("not implemented yet");
-        }
-
-        /*Find the max of the estimators*/
-        maxerror=elementerror[0];
-        for(int i=0;i<numberofelements;i++) maxerror=max(maxerror,elementerror[i]);
-        
-        /*Merge with the output*/
-   for(int i=0;i<numberofelements;i++){
-                if(elementerror[i]>threshold*maxerror) element_label[i]=1.;
-                else; //do nothing
-        }
-
-        /*Cleanup*/
-        xDelete<IssmDouble>(elementerror);
-}
-/*}}}*/
 void FemModel::GetElementDistanceToZeroLevelSet(IssmDouble** pelementdistance,int levelset_type){/*{{{*/
 
@@ -4970 +4879 @@
         
         /*Intermediaries*/
-   int numberofelements       = this->elements->NumberOfElements();
-   IssmDouble* levelset_points= NULL;
-   IssmDouble* xc                                       = NULL;
-   IssmDouble* yc                                       = NULL;
+   int numberofelements                                                 = this->elements->NumberOfElements();
+   Vector<IssmDouble>* velementdistance = new Vector<IssmDouble>(numberofelements);
+   IssmDouble* levelset_points                          = NULL;
+        IssmDouble* xyz_list                                                    = NULL;
+        IssmDouble mindistance,distance;
+        IssmDouble xc,yc,x1,y1,x2,y2,x3,y3;
         int numberofpoints;
-        IssmDouble distance;
-
-        /*Get element center coordinates*/
-        this->GetElementCenterCoordinates(&xc,&yc);
-        
-        /*Get points which level set is zero (center of elements with zero level set)*/ 
+        
+        /*Get points which level set is zero (center of elements with zero level set, levelset_points is serial)*/      
         this->GetZeroLevelSetPoints(&levelset_points,numberofpoints,levelset_type);
 
-        /*Find the minimal element distance to the zero levelset (grounding line or ice front)*/
-        elementdistance=xNew<IssmDouble>(numberofelements);
-        for(int i=0;i<numberofelements;i++){
-                elementdistance[i]=INFINITY;
+        for(int i=0;i<this->elements->Size();i++){//parallel
+      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
+      int sid = element->Sid();
+                element->GetVerticesCoordinates(&xyz_list); 
+                x1 = xyz_list[3*0+0];y1 = xyz_list[3*0+1];
+                x2 = xyz_list[3*1+0];y2 = xyz_list[3*1+1];
+                x3 = xyz_list[3*2+0];y3 = xyz_list[3*2+1];
+                xc = (x1+x2+x3)/3.;
+                yc = (y1+y2+y3)/3.;
+                mindistance=INFINITY;
+                /*Loop over each point (where level set is zero)*/
                 for(int j=0;j<numberofpoints;j++){
-                        distance=sqrt((xc[i]-levelset_points[2*j])*(xc[i]-levelset_points[2*j])+(yc[i]-levelset_points[2*j+1])*(yc[i]-levelset_points[2*j+1]));
-                        elementdistance[i]=min(distance,elementdistance[i]);            
-                }
+                        distance =sqrt((xc-levelset_points[2*j])*(xc-levelset_points[2*j])+(yc-levelset_points[2*j+1])*(yc-levelset_points[2*j+1]));
+                        mindistance=min(distance,mindistance);          
+                }
+                velementdistance->SetValue(sid,mindistance,INS_VAL);
+                xDelete<IssmDouble>(xyz_list);
         }       
 
+   /*Assemble*/
+   velementdistance->Assemble();
+
         /*Assign the pointer*/
-        (*pelementdistance)=elementdistance;
+        (*pelementdistance)=velementdistance->ToMPISerial();
 
         /*Cleanup*/
    xDelete<IssmDouble>(levelset_points);
-   xDelete<IssmDouble>(xc);
-   xDelete<IssmDouble>(yc);
+   xDelete<IssmDouble>(xyz_list);
+        delete velementdistance;
 }
 /*}}}*/

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

@@ -193 +193 @@
                 void ReMeshNeopz(int* pnewnumberofvertices,int* pnewnumberofelements,IssmDouble** pnewx,IssmDouble** pnewy,IssmDouble** pnewz,int** pnewelementslist);
                 void InitializeAdaptiveRefinementNeopz(void);
-                void GetPointsFromElementLabel(IssmDouble* element_label,int* numberofpoints,IssmDouble** xylist);
-                void GetElementLabelFromZeroLevelSet(IssmDouble* element_label,int levelset_type);
-                void GetElementLabelFromEstimators(IssmDouble* element_label,int estimator_type);       
                 void GetElementDistanceToZeroLevelSet(IssmDouble** pelementdistance,int levelset_type);
                 void SetRefPatterns(void);

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

@@ -26 +26 @@
         int        output_frequency;
         int        recording_frequency;
-        int        domaintype,groundingline_migration,smb_model,amr_frequency;
+        int        domaintype,groundingline_migration,smb_model,amr_frequency,amr_restart;
         int        numoutputs;
         Analysis  *analysis          = NULL;
@@ -65 +65 @@
         if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,TransientRequestedOutputsEnum);
 
-        #ifdef _HAVE_NEOPZ_
-        bool ismismip = false;//itapopo testing restart 
-        if(ismismip) femmodel->ReMesh();
+        #ifdef _HAVE_BAMG_ //#ifdef _HAVE_NEOPZ_ itapopo
+        if(amr_frequency){
+                femmodel->parameters->FindParam(&amr_restart,AmrRestartEnum);
+                if(amr_restart) femmodel->ReMesh();
+        }
         #endif
 

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

@@ -136 +136 @@
                                 parameters->AddObject(iomodel->CopyConstantObject("md.amr.icefront_distance",AmrIceFrontDistanceEnum));
                                 parameters->AddObject(iomodel->CopyConstantObject("md.amr.thicknesserror_threshold",AmrThicknessErrorThresholdEnum));
+                                parameters->AddObject(iomodel->CopyConstantObject("md.amr.thicknesserror_maximum",AmrThicknessErrorMaximumEnum));
                                 parameters->AddObject(iomodel->CopyConstantObject("md.amr.deviatoricerror_threshold",AmrDeviatoricErrorThresholdEnum));
+                                parameters->AddObject(iomodel->CopyConstantObject("md.amr.deviatoricerror_maximum",AmrDeviatoricErrorMaximumEnum));
+                                parameters->AddObject(iomodel->CopyConstantObject("md.amr.restart",AmrRestartEnum));
                                 break;
                         #endif
@@ -153 +156 @@
                                 parameters->AddObject(iomodel->CopyConstantObject("md.amr.thicknesserror_resolution",AmrThicknessErrorResolutionEnum));
                                 parameters->AddObject(iomodel->CopyConstantObject("md.amr.thicknesserror_threshold",AmrThicknessErrorThresholdEnum));
+                                parameters->AddObject(iomodel->CopyConstantObject("md.amr.thicknesserror_maximum",AmrThicknessErrorMaximumEnum));
                                 parameters->AddObject(iomodel->CopyConstantObject("md.amr.deviatoricerror_resolution",AmrDeviatoricErrorResolutionEnum));
                                 parameters->AddObject(iomodel->CopyConstantObject("md.amr.deviatoricerror_threshold",AmrDeviatoricErrorThresholdEnum));
+                                parameters->AddObject(iomodel->CopyConstantObject("md.amr.deviatoricerror_maximum",AmrDeviatoricErrorMaximumEnum));
+                                parameters->AddObject(iomodel->CopyConstantObject("md.amr.restart",AmrRestartEnum));
                                 /*Convert fieldname to enum and put it in params*/
                                 iomodel->FindConstant(&fieldname,"md.amr.fieldname");

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

@@ -870 +870 @@
         TransientAmrFrequencyEnum,
         AmrTypeEnum,
+        AmrRestartEnum,
         AmrNeopzEnum,
         AmrLevelMaxEnum,
@@ -887 +888 @@
         AmrThicknessErrorResolutionEnum,
         AmrThicknessErrorThresholdEnum,
+        AmrThicknessErrorMaximumEnum,
         AmrDeviatoricErrorResolutionEnum,
         AmrDeviatoricErrorThresholdEnum,
+        AmrDeviatoricErrorMaximumEnum,
         DeviatoricStressErrorEstimatorEnum,
         ThicknessErrorEstimatorEnum,

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

@@ -843 +843 @@
                 case TransientAmrFrequencyEnum : return "TransientAmrFrequency";
                 case AmrTypeEnum : return "AmrType";
+                case AmrRestartEnum : return "AmrRestart";
                 case AmrNeopzEnum : return "AmrNeopz";
                 case AmrLevelMaxEnum : return "AmrLevelMax";
@@ -860 +861 @@
                 case AmrThicknessErrorResolutionEnum : return "AmrThicknessErrorResolution";
                 case AmrThicknessErrorThresholdEnum : return "AmrThicknessErrorThreshold";
+                case AmrThicknessErrorMaximumEnum : return "AmrThicknessErrorMaximum";
                 case AmrDeviatoricErrorResolutionEnum : return "AmrDeviatoricErrorResolution";
                 case AmrDeviatoricErrorThresholdEnum : return "AmrDeviatoricErrorThreshold";
+                case AmrDeviatoricErrorMaximumEnum : return "AmrDeviatoricErrorMaximum";
                 case DeviatoricStressErrorEstimatorEnum : return "DeviatoricStressErrorEstimator";
                 case ThicknessErrorEstimatorEnum : return "ThicknessErrorEstimator";

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

@@ -19 +19 @@
       this.thicknesserror_resolution    = 500;
       this.thicknesserror_threshold     = 0;
+      this.thicknesserror_maximum               = 0;
       this.deviatoricerror_resolution   = 500;  
       this.deviatoricerror_threshold    = 0;    
+      this.deviatoricerror_maximum              = 0;    
         }// }}}
         this.disp= function(){// {{{
@@ -35 +37 @@
                 fielddisplay(this,'thicknesserror_resolution','element length when thickness error estimator is used');
                 fielddisplay(this,'thicknesserror_threshold','maximum threshold thickness error permitted');
+                fielddisplay(this,'thicknesserror_maximum','maximum thickness error permitted');
                 fielddisplay(this,'deviatoricerror_resolution','element length when deviatoric stress error estimator is used');
                 fielddisplay(this,'deviatoricerror_threshold','maximum threshold deviatoricstress error permitted');
+                fielddisplay(this,'deviatoricerror_maximum','maximum deviatoricstress error permitted');
         }// }}}
         this.classname= function(){// {{{
@@ -53 +57 @@
          checkfield(md,'fieldname','amr.thicknesserror_resolution','numel',[1],'>',0,'<',this.hmax,'NaN',1);
          checkfield(md,'fieldname','amr.thicknesserror_threshold','numel',[1],'>=',0,'<=',1,'NaN',1);
+         checkfield(md,'fieldname','amr.thicknesserror_maximum','numel',[1],'>=',0,'NaN',1,'Inf',1);
          checkfield(md,'fieldname','amr.deviatoricerror_resolution','numel',[1],'>',0,'<',this.hmax,'NaN',1);
          checkfield(md,'fieldname','amr.deviatoricerror_threshold','numel',[1],'>=',0,'<=',1,'NaN',1);
+         checkfield(md,'fieldname','amr.deviatoricerror_maximum','numel',[1],'>=',0,'NaN',1,'Inf',1);
                 } // }}}
                 this.marshall=function(md,prefix,fid) { //{{{
@@ -70 +76 @@
          WriteData(fid,prefix,'object',this,'fieldname','thicknesserror_resolution','format','Double');
          WriteData(fid,prefix,'object',this,'fieldname','thicknesserror_threshold','format','Double');
+         WriteData(fid,prefix,'object',this,'fieldname','thicknesserror_maximum','format','Double');
          WriteData(fid,prefix,'object',this,'fieldname','deviatoricerror_resolution','format','Double');
          WriteData(fid,prefix,'object',this,'fieldname','deviatoricerror_threshold','format','Double');
+         WriteData(fid,prefix,'object',this,'fieldname','deviatoricerror_maximum','format','Double');
                 }//}}}
                 this.fix=function() { //{{{
@@ -89 +97 @@
         this.thicknesserror_resolution  = 0.;
         this.thicknesserror_threshold           = 0.;
+        this.thicknesserror_maximum             = 0.;
         this.deviatoricerror_resolution = 0.;
         this.deviatoricerror_threshold  = 0.;
+        this.deviatoricerror_maximum            = 0.;
 
         this.setdefaultparameters();

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

@@ -18 +18 @@
                 thicknesserror_resolution = 0.;
                 thicknesserror_threshold = 0.;
+                thicknesserror_maximum = 0.;
                 deviatoricerror_resolution = 0.;
                 deviatoricerror_threshold = 0.;
+                deviatoricerror_maximum = 0.;
+                restart=0.;
         end
         methods (Static)
@@ -86 +89 @@
                         self.thicknesserror_resolution=500.;
                         self.thicknesserror_threshold=0.;
+                        self.thicknesserror_maximum=0.;
                         self.deviatoricerror_resolution=500.;
                         self.deviatoricerror_threshold=0.;
+                        self.deviatoricerror_maximum=0.;
+                        
+                        %is restart? This calls femmodel->ReMesh() before first time step. 
+                        self.restart=0;
 
                 end % }}}
@@ -103 +111 @@
                         md = checkfield(md,'fieldname','amr.thicknesserror_resolution','numel',[1],'>',0,'<',self.hmax,'NaN',1);
                         md = checkfield(md,'fieldname','amr.thicknesserror_threshold','numel',[1],'>=',0,'<=',1,'NaN',1);
+                        md = checkfield(md,'fieldname','amr.thicknesserror_maximum','numel',[1],'>=',0,'NaN',1,'Inf',1);
                         md = checkfield(md,'fieldname','amr.deviatoricerror_resolution','numel',[1],'>',0,'<',self.hmax,'NaN',1);
                         md = checkfield(md,'fieldname','amr.deviatoricerror_threshold','numel',[1],'>=',0,'<=',1,'NaN',1);
+                        md = checkfield(md,'fieldname','amr.deviatoricerror_maximum','numel',[1],'>=',0,'NaN',1,'Inf',1);
+                        md = checkfield(md,'fieldname','amr.restart','numel',[1],'>=',0,'<=',1,'NaN',1);
                 end % }}}
                 function disp(self) % {{{
@@ -120 +131 @@
                         fielddisplay(self,'thicknesserror_resolution',['element length when thickness error estimator is used']);
                         fielddisplay(self,'thicknesserror_threshold',['maximum threshold thickness error permitted']);
+                        fielddisplay(self,'thicknesserror_maximum',['maximum thickness error permitted']);
                         fielddisplay(self,'deviatoricerror_resolution',['element length when deviatoric stress error estimator is used']);
                         fielddisplay(self,'deviatoricerror_threshold',['maximum threshold deviatoricstress error permitted']);
-
+                        fielddisplay(self,'deviatoricerror_maximum',['maximum deviatoricstress error permitted']);
+                        fielddisplay(self,'restart',['indicates if ReMesh() will call before first time step']);
                 end % }}}
                 function marshall(self,prefix,md,fid) % {{{
@@ -139 +152 @@
                         WriteData(fid,prefix,'object',self,'class','amr','fieldname','thicknesserror_resolution','format','Double');
                         WriteData(fid,prefix,'object',self,'class','amr','fieldname','thicknesserror_threshold','format','Double');
+                        WriteData(fid,prefix,'object',self,'class','amr','fieldname','thicknesserror_maximum','format','Double');
                         WriteData(fid,prefix,'object',self,'class','amr','fieldname','deviatoricerror_resolution','format','Double');
                         WriteData(fid,prefix,'object',self,'class','amr','fieldname','deviatoricerror_threshold','format','Double');
+                        WriteData(fid,prefix,'object',self,'class','amr','fieldname','deviatoricerror_maximum','format','Double');
+                        WriteData(fid,prefix,'object',self,'class','amr','fieldname','restart','format','Integer');
 
                 end % }}}

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

@@ -24 +24 @@
         self.thicknesserror_resolution = 0.
         self.thicknesserror_threshold   = 0.
+        self.thicknesserror_maximum     = 0.
         self.deviatoricerror_resolution= 0.
         self.deviatoricerror_threshold = 0.
+        self.deviatoricerror_maximum    = 0.
         #set defaults
         self.setdefaultparameters()
@@ -42 +44 @@
         string="%s\n%s"%(string,fielddisplay(self,"thicknesserror_resolution","element length when thickness error estimator is used"))
         string="%s\n%s"%(string,fielddisplay(self,"thicknesserror_threshold","maximum threshold thickness error permitted"))
+        string="%s\n%s"%(string,fielddisplay(self,"thicknesserror_maximum","maximum thickness error permitted"))
         string="%s\n%s"%(string,fielddisplay(self,"deviatoricerror_resolution","element length when deviatoric stress error estimator is used"))
         string="%s\n%s"%(string,fielddisplay(self,"deviatoricerror_threshold","maximum threshold deviatoricstress error permitted"))
+        string="%s\n%s"%(string,fielddisplay(self,"deviatoricerror_maximum","maximum deviatoricstress error permitted"))
         return string
     #}}}
@@ -59 +63 @@
         self.thicknesserror_resolution = 500.
         self.thicknesserror_threshold   = 0
+        self.thicknesserror_maximum     = 0
         self.deviatoricerror_resolution= 500.
         self.deviatoricerror_threshold = 0
+        self.deviatoricerror_maximum    = 0
         return self
     #}}}
@@ -74 +80 @@
         md = checkfield(md,'fieldname','amr.thicknesserror_resolution','numel',[1],'>',0,'<',self.hmax,'NaN',1);
         md = checkfield(md,'fieldname','amr.thicknesserror_threshold','numel',[1],'>=',0,'<=',1,'NaN',1);
+        md = checkfield(md,'fieldname','amr.thicknesserror_maximum','numel',[1],'>=',0,'NaN',1,'Inf',1);
         md = checkfield(md,'fieldname','amr.deviatoricerror_resolution','numel',[1],'>',0,'<',self.hmax,'NaN',1);
         md = checkfield(md,'fieldname','amr.deviatoricerror_threshold','numel',[1],'>=',0,'<=',1,'NaN',1);        
+        md = checkfield(md,'fieldname','amr.deviatoricerror_maximum','numel',[1],'>=',0,'NaN',1,'Inf',1);
         return md
     # }}}
@@ -92 +100 @@
         WriteData(fid,prefix,'object',self,'fieldname','thicknesserror_resolution','format','Double');
         WriteData(fid,prefix,'object',self,'fieldname','thicknesserror_threshold','format','Double');
+        WriteData(fid,prefix,'object',self,'fieldname','thicknesserror_maximum','format','Double');
         WriteData(fid,prefix,'object',self,'fieldname','deviatoricerror_resolution','format','Double');
         WriteData(fid,prefix,'object',self,'fieldname','deviatoricerror_threshold','format','Double'); 
+        WriteData(fid,prefix,'object',self,'fieldname','deviatoricerror_maximum','format','Double'); 
     # }}}

issm/trunk-jpl/src/m/contrib/tsantos/AMRexportVTK.m

@@ -60 +60 @@
         
         timestep=step;
-
-        points=[model.results.TransientSolution(step).MeshX model.results.TransientSolution(step).MeshY zeros(size(model.results.TransientSolution(step).MeshX))];
+        if(isfield(model.results.TransientSolution,'MeshElements'))
+        index = model.results.TransientSolution(step).MeshElements;
+        x     = model.results.TransientSolution(step).MeshX;
+        y     = model.results.TransientSolution(step).MeshY;
+   else
+        index = model.mesh.elements;
+        x     = model.mesh.x;
+        y     = model.mesh.y;
+   end
+        
+        points=[x y zeros(size(x))];
         [num_of_points,dim]=size(points);
-        [num_of_elt]=size(model.results.TransientSolution(step).MeshElements,1);
+        [num_of_elt]=size(index,1);
 
         fid = fopen(strcat(path,filesep,name,filesep,'timestep.vtk',int2str(timestep),'.vtk'),'w+');
@@ -86 +95 @@
   end
         s=cell2mat(horzcat(s,{'\n'}));
-                fprintf(fid,s,[(point_per_elt)*ones(num_of_elt,1)       model.results.TransientSolution(step).MeshElements-1]');
+                fprintf(fid,s,[(point_per_elt)*ones(num_of_elt,1) index-1]');
         
         fprintf(fid,'CELL_TYPES %d\n',num_of_elt);

issm/trunk-jpl/src/m/contrib/tsantos/mismip/gl_position.m

@@ -3 +3 @@
                 %initialization of some variables
                 data                                    = md.results.TransientSolution(step).MaskGroundediceLevelset;
-                index                                   = md.results.TransientSolution(step).MeshElements;
-                x                                               = md.results.TransientSolution(step).MeshX;
-                y                                               = md.results.TransientSolution(step).MeshY;
+                if(isfield(md.results.TransientSolution,'MeshElements'))
+                        index                                   = md.results.TransientSolution(step).MeshElements;
+                        x                                               = md.results.TransientSolution(step).MeshX;
+                        y                                               = md.results.TransientSolution(step).MeshY;
+                else
+                        index                                   = md.mesh.elements;
+                        x                                               = md.mesh.x;
+                        y                                               = md.mesh.y;
+                end
                 numberofelements        = size(index,1);
                 elementslist            = 1:numberofelements;

issm/trunk-jpl/src/m/contrib/tsantos/mismip/ice_evolution.m

@@ -2 +2 @@
 %iv: ice volume
 %ivaf: ice volume above floatation
-%GLy40 : grounding line position @ y=40km
+%GL_y : grounding line position @ y (y comes in m)
 %nelem : number of elements
+% usage:
+%
+% Default: y=40km, i0=1
+% [ga iv ivaf GL_y nelem t] = ice_evolution(md);
+%
+% Default: y=40km
+% [ga iv ivaf GL_y nelem t] = ice_evolution(md,i0);
+%
+% Use this for y the borders (y=0 or y=ymax)
+% [ga iv ivaf GL_y nelem t] = ice_evolution(md,i0,y);
+%
+%
 
-function [ga iv ivaf GLy40 nelem t] = ice_evolution(md),
+function [ga iv ivaf GL_y nelem t] = ice_evolution(varargin),
 
         ga                      = [];
         iv                      = [];
         ivaf            = [];
-        GLy40           = [];
+        GL_y            = [];
         nelem           = [];
         t                       = [];
+
+        if(nargin==0)
+                error('it is necessary the model!')
+        elseif(nargin==1)
+                % Defoult is y=40km
+                i0      = 1;
+                y0      = 35000;
+                y1 = 45000;
+                dy = 100;
+                y  = 40000;
+        elseif(nargin==2)
+                i0=varargin{2};
+                % Defoult is y=40km
+                y0      = 35000;
+                y1 = 45000;
+                dy = 100;
+                y  = 40000;
+        elseif(nargin==3)
+                i0 = varargin{2};
+                y  = varargin{3};
+                dy = 10;
+                y0      = y-dy;
+                y1 = y+dy;
+        else
+                error('number of inputs is more than 3');
+        end
+        %set the model
+        md                      =varargin{1};
         nsteps  = length(md.results.TransientSolution);
 
-        for i=1:nsteps,
+
+        for i=i0:nsteps,
                 ga(i)                   = md.results.TransientSolution(i).GroundedArea;
                 iv(i)                   = md.results.TransientSolution(i).IceVolume;
                 ivaf(i)         = md.results.TransientSolution(i).IceVolumeAboveFloatation;
-                nelem(i)                = size(md.results.TransientSolution(i).MeshElements,1);
+                if(isfield(md.results.TransientSolution,'MeshElements'))
+                        nelem(i)                = size(md.results.TransientSolution(i).MeshElements,1);
+                else
+                        nelem(i)                = md.mesh.numberofelements;
+                end
                 t(i)                    = md.results.TransientSolution(i).time; 
-                %find GL position at y=40km
+                %find GL position between y0 and y1 
                 [glx gly]       = gl_position(md,i,0);
-                pos=find(gly<45000 & gly > 35000);
-                x=gly(pos);
-                v=glx(pos);
-                xq=[38000:100:42000];
-                vq = interp1(x,v,xq,'linear');
-                pos=find(xq==40000);
-                GLy40(i)=vq(pos);
+                pos                     = find(gly<y1 & gly>y0);
+                x                               = gly(pos);
+                v                               = glx(pos);
+                if(length(pos)==0)
+                        error('pos is null')
+                elseif(length(pos)==1)
+                        %this should be used for y=0 or y=ymax
+                        GL_y(i) = v;
+                else
+                        %this should be used when y is inside the domain; so, use linear interpolation
+                        xq                      = [y0:dy:y1];
+                        vq                      = interp1(x,v,xq,'linear');
+                        pos             = find(xq==y);
+                        if(pos)
+                                GL_y(i) = vq(pos);
+                        else
+                                error('pos is null')
+                        end
+                end
+
         end
 

issm/trunk-jpl/src/m/contrib/tsantos/remesh.m

@@ -58 +58 @@
 NewModel.geometry.surface                               = md.results.TransientSolution(end).Surface;
 NewModel.geometry.base                                  = md.results.TransientSolution(end).Base;
-%NewModel.geometry.bed                                  = md.geometry.bed; %use from parameterize
+NewModel.geometry.bed                                   = md.results.TransientSolution(end).Bed;%md.geometry.bed; %use from parameterize
 NewModel.geometry.thickness                     = md.results.TransientSolution(end).Thickness;
 NewModel.mask.groundedice_levelset  = md.results.TransientSolution(end).MaskGroundediceLevelset;
@@ -72 +72 @@
 NewModel.cluster                = md.cluster;
 NewModel.transient              = md.transient;
+NewModel.amr                    = md.amr;
 
 mdOut = NewModel;