Changeset 22100

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

-              r21919
+              r22100
         this->CleanUp();
         /*Copy all data*/
+        this->fathermesh                        = new TPZGeoMesh(*cp.fathermesh);
+        this->currentmesh                       = new TPZGeoMesh(*cp.currentmesh);
+        this->levelmax                          = cp.levelmax;
+        this->elementswidth             = cp.elementswidth;
+        this->regionlevel1              = cp.regionlevel1;
+        this->regionlevelmax            = cp.regionlevelmax;
+        this->fathermesh                                                = new TPZGeoMesh(*cp.fathermesh);
+        this->previousmesh                                      = new TPZGeoMesh(*cp.previousmesh);
+        this->level_max                                         = cp.level_max;
+        this->radius_level_max                          = cp.radius_level_max;
+        this->gradation                                         = cp.gradation;
+        this->lag                                                               = cp.lag;
+   this->groundingline_distance         = cp.groundingline_distance;
+        this->icefront_distance                         = cp.icefront_distance;
+        this->thicknesserror_threshold  = cp.thicknesserror_threshold;
+        this->deviatoricerror_threshold = cp.deviatoricerror_threshold;
         this->sid2index.clear();
         this->sid2index.resize(cp.sid2index.size());
 …
         /*Verify and delete all data*/
         if(this->fathermesh)    delete this->fathermesh;
+        if(this->currentmesh)   delete this->currentmesh;
+        this->levelmax                          = -1;
+        this->elementswidth             = -1;
+        this->regionlevel1              = -1;
+        this->regionlevelmax            = -1;
+        if(this->previousmesh)  delete this->previousmesh;
+        this->level_max                                         = -1;
+        this->radius_level_max                          = -1;
+        this->gradation                                         = -1;
+        this->lag                                                               = -1;
+   this->groundingline_distance         = -1;
+        this->icefront_distance                         = -1;
+        this->thicknesserror_threshold  = -1;
+        this->deviatoricerror_threshold = -1;
         this->sid2index.clear();
         this->index2sid.clear();
 …
         /*Set pointers to NULL*/
+        this->fathermesh                        = NULL;
+        this->currentmesh                       = NULL;
+        this->levelmax                          = -1;
+        this->elementswidth             = -1;
+        this->regionlevel1              = -1;
+        this->regionlevelmax            = -1;
+        this->step                                      = 0;
+        this->smooth_frequency  = 1;
+        this->fathermesh                                                = NULL;
+        this->previousmesh                                      = NULL;
+        this->level_max                                         = -1;
+        this->radius_level_max                          = -1;
+        this->gradation                                         = -1;
+        this->lag                                                               = -1;
+   this->groundingline_distance         = -1;
+        this->icefront_distance                         = -1;
+        this->thicknesserror_threshold  = -1;
+        this->deviatoricerror_threshold = -1;
         this->sid2index.clear();
         this->index2sid.clear();
 …
 /*Mesh refinement methods*/
+void AdaptiveMeshRefinement::Execute(bool &amr_verbose,
+                                                                                                int &numberofelements,
+                                                                                                double* partiallyfloatedelements,
+                                                                                                double *masklevelset,
+                                                                                                double* deviatorictensorerror,
+                                                                                                double* thicknesserror,
+                                                                                                int &newnumberofvertices,int &newnumberofelements,double** x,double** y,int** elementslist){/*{{{*/
+        /*IMPORTANT! newelements are in Matlab indexing*/
+void AdaptiveMeshRefinement::ExecuteRefinement(int numberofpoints,double* xylist,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->currentmesh) _error_("Impossible to execute refinement: fathermesh or currentmesh is NULL!\n");
+        if(numberofelements!=this->sid2index.size()) _error_("Impossible to execute refinement: sid2index.size is not equal to numberofelements!\n");
+        /*Combine the fields*/
+        double mean_mask                = 0;
+        double mean_tauerror = 0;
+        double mean_Herror      = 0;
+        int* fmask                              = xNew<int>(numberofelements);
+        int* ftauerror                  = xNew<int>(numberofelements);
+        int* fHerror                    = xNew<int>(numberofelements);
+        int* phi                                        = xNew<int>(numberofelements);
+        /*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();
+        /*Filter to thickness*/
+        for(int i=0;i<this->sid2index.size();i++){
+                fmask[i]=0;
+                if(thicknesserror[i]>mean_Herror) fmask[i]=1;
+        }
+        /*Filter to tau*/
+        for(int i=0;i<this->sid2index.size();i++){
+                ftauerror[i]=0;
+                if(deviatorictensorerror[i]>mean_tauerror) ftauerror[i]=1;
+        }
+        /*Sum*/
+        for(int i=0;i<this->sid2index.size();i++){
+                phi[i]=ftauerror[i]+fmask[i];
+        }
+        /*Execute the refinement.*/
+        this->RefinementProcess(amr_verbose,partiallyfloatedelements,masklevelset,deviatorictensorerror,thicknesserror);
+        if(!this->fathermesh || !this->previousmesh) _error_("Impossible to execute refinement: fathermesh or previousmesh 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);
         /*Get new geometric mesh in ISSM data structure*/
+        this->GetMesh(newnumberofvertices,newnumberofelements,x,y,elementslist);
+        std::ofstream file("/home/santos/mesh.vtk");
+        TPZVTKGeoMesh::PrintGMeshVTK(this->currentmesh,file);
+        this->GetMesh(this->previousmesh,pnewnumberofvertices,pnewnumberofelements,px,py,pelementslist);
         /*Verify the new geometry*/
+        this->CheckMesh(newnumberofvertices,newnumberofelements,this->elementswidth,x,y,elementslist);
+}
+/*}}}*/
+void AdaptiveMeshRefinement::RefinementProcess(bool &amr_verbose,double* partiallyfloatedelements,double* masklevelset,double* deviatorictensorerror,double* thicknesserror){/*{{{*/
+        this->CheckMesh(pnewnumberofvertices,pnewnumberofelements,px,py,pelementslist);
+}
+/*}}}*/
+void AdaptiveMeshRefinement::RefinementProcess(bool &verbose,double* gl_elementdistance,double* if_elementdistance,double* deviatoricerror,double* thicknesserror){/*{{{*/
         if(amr_verbose) _printf_("\n\trefinement process started (level max = " << this->levelmax << ")\n");
+        if(verbose) _printf_("\n\trefinement process started (level max = " << this->level_max << ")\n");
         /*Intermediaries*/
         TPZGeoMesh* nohangingnodesmesh=NULL;
-        //itapopo
-        this->step++;
         double mean_mask                = 0;
 …
         /*Calculate mean values*/
+        /*
         for(int i=0;i<this->sid2index.size();i++){
                 mean_mask               += masklevelset[i];
 …
         mean_tauerror   /= this->sid2index.size();
         mean_Herror             /= this->sid2index.size();
         if(amr_verbose) _printf_("\t\tdeal with special elements...\n");
+        */
+        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->currentmesh->Element(this->specialelementsindex[i]);
+                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");
 …
                                 /*Ok, the special element can be deleted*/
                                 siblings[j]->ResetSubElements();
                                 this->currentmesh->DeleteElement(siblings[j],siblings[j]->Index());
+                                this->previousmesh->DeleteElement(siblings[j],siblings[j]->Index());
+                        }
+                }
                 /*Now, verify if the father should be refined with uniform pattern (smoother)*/
+                this->smooth_frequency=10000000;
+                if(siblings.size()==3 || this->step%this->smooth_frequency==0){//it keeps the mesh with uniform elements
+                if(siblings.size()==3){//it keeps the mesh with uniform elements
                         /*Father has uniform subelements now*/
                         TPZVec<TPZGeoEl *> sons;
                         father->Divide(sons);
-                        this->smooth_frequency=this->step;
                 }else{
                         specialfatherindex.push_back(father->Index());
 …
                 if(this->specialelementsindex[i]!=-1) _error_("special element "<<i<<" was not deleted!\n");
+        }
         this->currentmesh->BuildConnectivity();
         if(amr_verbose) _printf_("\t\tuniform refinement...\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->currentmesh->Element(this->sid2index[i]);
+                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");
 …
+                        }
                         TPZVec<TPZGeoEl *> sons;
                         if(geoel->Level()<this->levelmax && count==0) geoel->Divide(sons);
+                        if(geoel->Level()<this->level_max && count==0) geoel->Divide(sons);
+                }
                 else if(geoel->Level()>0)
 …
                         if(group_Herror>0 && group_Herror<0*mean_Herror){ //itapopo
                                 /*Reset subelements in the father*/
                                 this->currentmesh->Element(geoel->Father()->Index())->ResetSubElements();
+                                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++){
 …
                                         this->index2sid[index]=-1;
                                         if(sid!=-1) this->sid2index[sid]=-1;
                                         this->currentmesh->DeleteElement(siblings[j],siblings[j]->Index());
+                                        this->previousmesh->DeleteElement(siblings[j],siblings[j]->Index());
                                 }//for j
                         }//if
                 }/*Unrefine process*/
         }//for i
         this->currentmesh->BuildConnectivity();
         if(amr_verbose) _printf_("\t\trefine to avoid hanging nodes...\n");
         this->RefineMeshToAvoidHangingNodes(this->currentmesh);
+        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(amr_verbose) _printf_("\trefinement process done!\n");
+}
+/*}}}*/
+void AdaptiveMeshRefinement::RefineMesh(TPZGeoMesh *gmesh,std::vector<int> &elements){/*{{{*/
+        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.;;
+        TPZVec<REAL> qsi(2,0.),cp(3,0.);
+        TPZVec<TPZGeoEl*> sons;
+        /*First, delete the special elements*/
+        this->DeleteSpecialElements(verbose,gmesh);
+        /*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: ");
+                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
+                        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 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*/
 …
+}
 /*}}}*/
 void AdaptiveMeshRefinement::RefineMeshToAvoidHangingNodes(TPZGeoMesh *gmesh){/*{{{*/
+void AdaptiveMeshRefinement::RefineMeshToAvoidHangingNodes(bool &verbose,TPZGeoMesh* gmesh){/*{{{*/
         /*Now, insert special elements to avoid hanging nodes*/
+        if(verbose) _printf_("\trefining to avoid hanging nodes (total: ");
+        /*Intermediaries*/
+        int nelem=-1;
+        int count= 1;
+        while(count>0){
+                nelem=gmesh->NElements();//must keep here
+                count=0;
+                for(int i=0;i<nelem;i++){
+                        /*Get geometric element and verify if it has already been refined. Geoel may not have been previously refined*/
+                        TPZGeoEl * geoel=gmesh->Element(i);
+                        if(!geoel) continue;
+                        if(geoel->HasSubElement()) continue;
+                        if(geoel->MaterialId() != this->GetElemMaterialID()) continue;
+                        /*Get the refinement pattern for this element and refine it*/
+                        TPZAutoPointer<TPZRefPattern> refp=TPZRefPatternTools::PerfectMatchRefPattern(geoel);
+                        if(refp){
+                                /*Non-uniform refinement*/
+                                TPZVec<TPZGeoEl *> sons;
+                                geoel->SetRefPattern(refp);
+                                geoel->Divide(sons);
+                                count++;
+                                /*Keep the index of the special elements*/
+                                for(int j=0;j<sons.size();j++) this->specialelementsindex.push_back(sons[j]->Index());
+                        }
+                }
+                if(verbose){
+                        if(count==0) _printf_(""<<count<<")\n");
+                        else _printf_(""<<count<<", ");
+                }
+                gmesh->BuildConnectivity();
+        }
+}
+/*}}}*/
+void AdaptiveMeshRefinement::DeleteSpecialElements(bool &verbose,TPZGeoMesh* gmesh){/*{{{*/
+        /*Intermediaries*/
+        int count=0;
+        if(verbose) _printf_("\tdelete "<<this->specialelementsindex.size()<<" special elements (total: ");
+        for(int i=0;i<this->specialelementsindex.size();i++){
+                if(this->specialelementsindex[i]==-1) continue;
+                if(!gmesh->Element(this->specialelementsindex[i])) continue;
+                if(gmesh->Element(this->specialelementsindex[i])->Father()) gmesh->Element(this->specialelementsindex[i])->Father()->ResetSubElements();
+                gmesh->DeleteElement(gmesh->Element(this->specialelementsindex[i]),this->specialelementsindex[i]);
+                count++;
+        }
+        if(verbose) _printf_(""<<count<<")\n");
+        /*Cleanup*/
         this->specialelementsindex.clear();
+        const int NElem = gmesh->NElements();
+        for(int i=0;i<NElem;i++){
+                /*Get geometric element and verify if it has already been refined. Geoel may not have been previously refined*/
+                TPZGeoEl * geoel=gmesh->Element(i);
+                if(!geoel) continue;
+                if(geoel->HasSubElement()) continue;
+                if(geoel->MaterialId() != this->GetElemMaterialID()) continue;
+                /*Get the refinement pattern for this element and refine it*/
+                TPZAutoPointer<TPZRefPattern> refp=TPZRefPatternTools::PerfectMatchRefPattern(geoel);
+                if(refp){
+                        /*Non-uniform refinement*/
+                        TPZVec<TPZGeoEl *> Sons;
+                        geoel->SetRefPattern(refp);
+                        geoel->Divide(Sons);
+                        /*Keep the index of the special elements*/
+                        for(int j=0;j<Sons.size();j++) this->specialelementsindex.push_back(Sons[j]->Index());
+                }
+        }
+        /*Adjust connectivities*/
         gmesh->BuildConnectivity();
+}
 /*}}}*/
 void AdaptiveMeshRefinement::GetMesh(int &nvertices,int &nelements,double** px,double** py, int** pelements){/*{{{*/
+void AdaptiveMeshRefinement::GetMesh(TPZGeoMesh* gmesh,int* nvertices,int* nelements,double** px,double** py, int** pelements){/*{{{*/
         /* IMPORTANT! pelements are in Matlab indexing
 …
         /*Intermediaries */
-        TPZGeoMesh* gmesh = this->currentmesh;//itapopo confirmar
         long sid,nodeindex;
+        int nconformelements,nconformvertices;
+        int nconformelements,nconformvertices;
+        int ntotalvertices              = gmesh->NNodes();
         int* newelements                        = NULL;
         double* newmeshX                        = NULL;//xNew<double>(ntotalvertices);
         double* newmeshY                        = NULL;//xNew<double>(ntotalvertices);
+        double* newmeshX                        = NULL;
+        double* newmeshY                        = NULL;
         TPZGeoEl* geoel                 = NULL;
         long* vertex_index2sid  = xNew<long>(gmesh->NNodes());
+        long* vertex_index2sid  = xNew<long>(ntotalvertices);
         this->index2sid.clear(); this->index2sid.resize(gmesh->NElements());
         this->sid2index.clear();
-        /*Get mesh coords */
-        //for(int i=0;i<ntotalvertices;i++ ){
-        //      TPZVec<REAL> coords(3,0.);
-        //      gmesh->NodeVec()[i].GetCoordinates(coords);
-        //      newmeshX[i] = coords[0];
-        //      newmeshY[i] = coords[1];
-        //}
         /*Fill in the vertex_index2sid vector with non usual index value*/
         for(int i=0;i<gmesh->NNodes();i++) vertex_index2sid[i]=-1;
 …
                 if(geoel->HasSubElement()) continue;
                 if(geoel->MaterialId() != this->GetElemMaterialID()) continue;
                 this->sid2index.push_back(i);//keep the element index
                 this->index2sid[i]=this->sid2index.size()-1;//keep the element sid
                 for(int j=0;j<this->elementswidth;j++){
+                this->sid2index.push_back(geoel->Index());//keep the element index
+                this->index2sid[geoel->Index()]=this->sid2index.size()-1;//keep the element sid
+                for(int j=0;j<this->GetNumberOfNodes();j++){
         nodeindex=geoel->NodeIndex(j);
+        if(vertex_index2sid[nodeindex]==-1){
+        if(nodeindex<0) _error_("nodeindex is <0\n");
+                        if(vertex_index2sid[nodeindex]==-1){
                 vertex_index2sid[nodeindex]=sid;
                                 sid++;
 …
         nconformelements        = (int)this->sid2index.size();
         nconformvertices        = (int)sid;
         newelements                     = xNew<int>(nconformelements*this->elementswidth);
+        newelements                     = xNew<int>(nconformelements*this->GetNumberOfNodes());
         newmeshX                                = xNew<double>(nconformvertices);
    newmeshY                             = xNew<double>(nconformvertices);
+        for(int i=0;i<nconformvertices;i++){//over the TPZNode index (fill in the ISSM vertices coords)
+                sid = vertex_index2sid[i];
+                if(sid!=-1){
+                        TPZVec<REAL> coords(3,0.);
+                        gmesh->NodeVec()[i].GetCoordinates(coords);
+                        newmeshX[sid] = coords[0];
+                        newmeshY[sid] = coords[1];
+                }
+        for(int i=0;i<ntotalvertices;i++){//over the TPZNode index (fill in the ISSM vertices coords)
+                sid=vertex_index2sid[i];
+                if(sid==-1) continue;//skip this index (node no used)
+                TPZVec<REAL> coords(3,0.);
+                gmesh->NodeVec()[i].GetCoordinates(coords);
+                newmeshX[sid] = coords[0];
+                newmeshY[sid] = coords[1];
+        }
         for(int i=0;i<this->sid2index.size();i++){//over the sid (fill the ISSM elements)
                 for(int j=0;j<this->elementswidth;j++) {
+                for(int j=0;j<this->GetNumberOfNodes();j++) {
                         geoel   = gmesh->ElementVec()[this->sid2index[i]];
                         sid     = vertex_index2sid[geoel->NodeIndex(j)];
                         newelements[i*this->elementswidth+j]=(int)sid+1;//C to Matlab indexing
+                        newelements[i*this->GetNumberOfNodes()+j]=(int)sid+1;//C to Matlab indexing
+                }
                 /*Verify the Jacobian determinant. If detJ<0, swap the 2 first postions:
 …
                 int a,b,c;
                 a=newelements[i*this->elementswidth+0]-1;
                 b=newelements[i*this->elementswidth+1]-1;
                 c=newelements[i*this->elementswidth+2]-1;
+                a=newelements[i*this->GetNumberOfNodes()+0]-1;
+                b=newelements[i*this->GetNumberOfNodes()+1]-1;
+                c=newelements[i*this->GetNumberOfNodes()+2]-1;
                 xa=newmeshX[a]; ya=newmeshY[a];
 …
                 /*verify and swap, if necessary*/
                 if(detJ<0) {
                         newelements[i*this->elementswidth+0]=b+1;//a->b
                         newelements[i*this->elementswidth+1]=a+1;//b->a
+                        newelements[i*this->GetNumberOfNodes()+0]=b+1;//a->b
+                        newelements[i*this->GetNumberOfNodes()+1]=a+1;//b->a
+                }
+        }
         /*Setting outputs*/
         nvertices       = nconformvertices;
         nelements       = nconformelements;
+        *nvertices      = nconformvertices;
+        *nelements      = nconformelements;
         *px                     = newmeshX;
         *py                = newmeshY;
 …
+}
 /*}}}*/
+void AdaptiveMeshRefinement::FindElements(int &numberofpoints,double* xp,double* yp,TPZGeoMesh *gmesh,int &hlevel,std::vector<int> &elements){/*{{{*/
+        if(!gmesh) _error_("Impossible to set elements: gmesh is NULL!\n");
+        if(false){
+                this->AllElements(gmesh,elements); //uniform, refine all elements!
+                return;
+        }
+        /*Intermediaries*/
+        elements.clear();
+        double D1               = this->regionlevel1;
+        double Dhmax    = this->regionlevelmax;
+        int hmax                        = this->levelmax;
+        double alpha    = (hmax==1) ? 0. : log(D1/Dhmax)/(hmax-1.);
+        double Di               = D1/exp(alpha*(hlevel-1));
+        int side2D              = 6;
+        double distance,value;
+        /*Find elements near the points */
+        for(int i=0;i<gmesh->NElements();i++){
+                if(gmesh->Element(i)->MaterialId()!=this->GetElemMaterialID()) continue;
+                if(gmesh->Element(i)->HasSubElement()) continue;
+                if(gmesh->Element(i)->Level()>=hlevel) continue;
+                TPZVec<REAL> qsi(2,0.);
+                TPZVec<REAL> centerPoint(3,0.);
+                gmesh->Element(i)->CenterPoint(side2D, qsi);
+                gmesh->Element(i)->X(qsi, centerPoint);
+                distance = Di;
+                for (int j=0;j<numberofpoints;j++){
+                        value = std::sqrt( (xp[j]-centerPoint[0])*(xp[j]-centerPoint[0])+(yp[j]-centerPoint[1] )*(yp[j]-centerPoint[1]) );//sqrt( (x2-x1)^2 + (y2-y1)^2 )
+                        if(value<distance) distance=value; //min distance to the point
+                }
+                if(distance<Di) elements.push_back(i);
+        }
+}
+/*}}}*/
+void AdaptiveMeshRefinement::AllElements(TPZGeoMesh *gmesh,std::vector<int> &elements){/*{{{*/
+    /* Uniform refinement. This refines the entire mesh */
+    int nelements = gmesh->NElements();
+         elements.clear();
+    for(int i=0;i<nelements;i++){
+        if(gmesh->Element(i)->MaterialId()!=this->GetElemMaterialID()) continue;
+        if(gmesh->Element(i)->HasSubElement()) continue;
+        elements.push_back(i);
+    }
+}
+/*}}}*/
+void AdaptiveMeshRefinement::CreateInitialMesh(int &nvertices,int &nelements,int &width,double* x,double* y,int* elements){/*{{{*/
+void AdaptiveMeshRefinement::CreateInitialMesh(int &nvertices,int &nelements,double* x,double* y,int* elements){/*{{{*/
         /* IMPORTANT! elements come in Matlab indexing
 …
         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");
-        this->SetElementWidth(width);
     /*Verify and creating initial mesh*/
    if(this->fathermesh || this->currentmesh) _error_("Initial mesh already exists!");
+   if(this->fathermesh || this->previousmesh) _error_("Initial mesh already exists!");
    this->fathermesh = new TPZGeoMesh();
 …
         long index;
    const int mat = this->GetElemMaterialID();
    TPZManVector<long> elem(this->elementswidth,0);
+   TPZManVector<long> elem(this->GetNumberOfNodes(),0);
    this->sid2index.clear();
         for(int i=0;i<nelements;i++){
                 for(int j=0;j<this->elementswidth;j++) elem[j]=elements[i*this->elementswidth+j]-1;//Convert Matlab to C indexing
+                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->elementswidth){
+      switch(this->GetNumberOfNodes()){
                         case 3: this->fathermesh->CreateGeoElement(ETriangle,elem,mat,index,reftype);   break;
          default:       _error_("mesh not supported yet");
 …
    /*Build element and node connectivities*/
    this->fathermesh->BuildConnectivity();
         /*Set current mesh*/
         this->currentmesh=new TPZGeoMesh(*this->fathermesh);
+        /*Set previous mesh*/
+        this->previousmesh=new TPZGeoMesh(*this->fathermesh);
+}
 /*}}}*/
 …
+}
 /*}}}*/
+void AdaptiveMeshRefinement::SetLevelMax(int &h){/*{{{*/
+    this->levelmax = h;
+}
+/*}}}*/
+void AdaptiveMeshRefinement::SetRegions(double &D1,double Dhmax){/*{{{*/
+    this->regionlevel1   = D1;
+    this->regionlevelmax = Dhmax;
+}
+/*}}}*/
+void AdaptiveMeshRefinement::SetElementWidth(int &width){/*{{{*/
+        if(width!=3) _error_("elementswidth not supported yet!");
+   this->elementswidth = width;
+}
+/*}}}*/
+void AdaptiveMeshRefinement::CheckMesh(int &nvertices,int &nelements,int &width,double** px,double** py,int** pelements){/*{{{*/
+void AdaptiveMeshRefinement::CheckMesh(int* nvertices,int* nelements,double** px,double** py,int** pelements){/*{{{*/
         /*Basic verification*/
         if(nvertices<=0) _error_("Impossible to continue: nvertices <=0!\n");
         if(nelements<=0) _error_("Impossible to continue: nelements <=0!\n");
-        if(width!=3) _error_("Impossible to continue: width !=3!\n");
         if(!px) _error_("Impossible to continue: px is NULL!\n");
         if(!py) _error_("Impossible to continue: py is NULL!\n");
 …
         std::set<int> elemvertices;
         elemvertices.clear();
         for(int i=0;i<nelements;i++){
                 for(int j=0;j<width;j++) {
                         elemvertices.insert((*pelements)[i*width+j]);
+                }
+        }
         if(elemvertices.size()!=nvertices) _error_("Impossible to continue: elemvertices.size() != nvertices!\n");
+        for(int i=0;i<*nelements;i++){
+                for(int j=0;j<this->GetNumberOfNodes();j++) {
+                        elemvertices.insert((*pelements)[i*this->GetNumberOfNodes()+j]);
+                }
+        }
+        if(elemvertices.size()!=*nvertices) _error_("Impossible to continue: elemvertices.size() != nvertices!\n");
         //Verify if there are inf or NaN in coords
         for(int i=0;i<nvertices;i++){
+        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");
+        }
+        for(int i=0;i<nelements;i++){
+                for(int j=0;j<width;j++){
+                        if(isnan((*pelements)[i*width+j]) || isinf((*pelements)[i*width+j]) ) _error_("Impossible to continue: px 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");
+                }
+        }
 …
+}
 /*}}}*/
+void AdaptiveMeshRefinement::PrintGMeshVTK(TPZGeoMesh* gmesh,std::ofstream &file,bool matColor){/*{{{*/
+        file.clear();
+        long nelements = gmesh->NElements();
+        TPZGeoEl *gel;
+        std::stringstream node, connectivity, type, material;
+        //Header
+        file << "# vtk DataFile Version 3.0" << std::endl;
+        file << "TPZGeoMesh VTK Visualization" << std::endl;
+        file << "ASCII" << std::endl << std::endl;
+        file << "DATASET UNSTRUCTURED_GRID" << std::endl;
+        file << "POINTS ";
+        long actualNode = -1, size = 0, nVALIDelements = 0;
+        for(long el = 0; el < nelements; el++){
+          gel = gmesh->ElementVec()[el];
+          if(!gel )//|| (gel->Type() == EOned && !gel->IsLinearMapping()))//Exclude Arc3D and Ellipse3D
+          {
+                        continue;
+          }
+          if(gel->HasSubElement())
+          {
+                        continue;
+          }
+          MElementType elt = gel->Type();
+          int elNnodes = MElementType_NNodes(elt);
+          size += (1+elNnodes);
+          connectivity << elNnodes;
+          for(int t = 0; t < elNnodes; t++)
+          {
+                        for(int c = 0; c < 3; c++)
+                        {
+                                 double coord = gmesh->NodeVec()[gel->NodeIndex(t)].Coord(c);
+                                 node << coord << " ";
+                        }
+                        node << std::endl;
+                        actualNode++;
+                        connectivity << " " << actualNode;
+          }
+          connectivity << std::endl;
+          int elType = this->GetVTK_ElType(gel);
+          type << elType << std::endl;
+          if(matColor == true)
+          {
+                        material << gel->MaterialId() << std::endl;
+          }
+          else
+          {
+                        material << gel->Index() << std::endl;
+          }
+          nVALIDelements++;
+        }
+        node << std::endl;
+        actualNode++;
+        file << actualNode << " float" << std::endl << node.str();
+        file << "CELLS " << nVALIDelements << " ";
+        file << size << std::endl;
+        file << connectivity.str() << std::endl;
+        file << "CELL_TYPES " << nVALIDelements << std::endl;
+        file << type.str() << std::endl;
+        file << "CELL_DATA" << " " << nVALIDelements << std::endl;
+        file << "FIELD FieldData 1" << std::endl;
+        if(matColor == true)
+        {
+          file << "material 1 " << nVALIDelements << " int" << std::endl;
+        }
+        else
+        {
+          file << "ElementIndex 1 " << nVALIDelements << " int" << std::endl;
+        }
+        file << material.str();
+        file.close();
+}
+/*}}}*/
+int AdaptiveMeshRefinement::GetVTK_ElType(TPZGeoEl * gel){/*{{{*/
+        MElementType pzElType = gel->Type();
+    int elType = -1;
+    switch (pzElType)
+    {
+        case(EPoint):
+        {
+            elType = 1;
+            break;
+        }
+        case(EOned):
+        {
+            elType = 3;
+            break;
+        }
+        case (ETriangle):
+        {
+            elType = 5;
+            break;
+        }
+        case (EQuadrilateral):
+        {
+            elType = 9;
+            break;
+        }
+        case (ETetraedro):
+        {
+            elType = 10;
+            break;
+        }
+        case (EPiramide):
+        {
+            elType = 14;
+            break;
+        }
+        case (EPrisma):
+        {
+            elType = 13;
+            break;
+        }
+        case (ECube):
+        {
+            elType = 12;
+            break;
+        }
+        default:
+        {
+            std::cout << "Element type not found on " << __PRETTY_FUNCTION__ << std::endl;
+            DebugStop();
+            break;
+        }
+    }
+    if(elType == -1)
+    {
+        std::cout << "Element type not found on " << __PRETTY_FUNCTION__ << std::endl;
+        std::cout << "MIGHT BE CURVED ELEMENT (quadratic or quarter point)" << std::endl;
+        DebugStop();
+    }
+    return elType;
+}
+/*}}}*/

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

-              r21919
+              r22100
 #include <pzgeotriangle.h>
 #include <tpzgeoelrefpattern.h>
+#include <pzgraphmesh.h>
 #include <TPZVTKGeoMesh.h>
 …
 public:
+        /*Public methods*/
+        /*Public attributes{{{*/
+        /* Filter (distance to the points)
+         * to refine:
+         * radius_h = initial_radius * gradation ^ (level_max-h)
+         * to unirefine
+         * radius_h = lag * initial_radius * gradation ^ (level_max-h)
+         */
+        int level_max;                                                          //max level of refinement
+        double radius_level_max;                                //initial radius which in the elements will be refined with level max
+        double gradation;                                                       //geometric progression ratio to calculate radius of level h
+        double lag;                                                                     //lag used in the unrefine process
+        /*Target and estimators*/
+        double groundingline_distance;          //if groundingline_distance>initial_radius, groundingline_distance will be used instead initial_radius
+        double icefront_distance;                               //if icefront_distance>initial_radius, icefront_distance will be used instead initial_radius
+        double thicknesserror_threshold;                //if ==0, it will not be used
+        double deviatoricerror_threshold;       //if ==0, it will not be used
+        /*}}}*/
+        /*Public methods{{{*/
         /* Constructor, destructor etc*/
         AdaptiveMeshRefinement();
 …
         AdaptiveMeshRefinement & operator= (const AdaptiveMeshRefinement &cp);
         virtual ~AdaptiveMeshRefinement();
         /*General methods*/
         void CleanUp();
         void Initialize();
+        void SetLevelMax(int &h);
+   void SetRegions(double &D1,double Dhmax);
+        void SetElementWidth(int &width);
+        void Execute(bool &amr_verbose,int &numberofelements,
+                                                double* partiallyfloatedelements,double *masklevelset,double* deviatorictensorerror,double* thicknesserror,
+                                                int &newnumberofvertices,int &newnumberofelements,double** x,double** y,int** elementslist);
+        void CreateInitialMesh(int &nvertices,int &nelements,int &width,double* x,double* y,int* elements);
+        TPZGeoMesh* CreateRefPatternMesh(TPZGeoMesh* gmesh);
+        void CheckMesh(int &nvertices,int &nelements,int &width,double** px,double** py,int** pelements);
+        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 CreateInitialMesh(int &nvertices,int &nelements,double* x,double* y,int* elements);
+        void CheckMesh(int* nvertices,int* nelements,double** px,double** py,int** pelements);
+        /*}}}*/
 private:
+        /*Private attributes*/
+        int elementswidth;                                    // Geometric nodes for element: 3 == Tria, 4 == Tetra, 6 == Penta
+   int levelmax;                                         // Max level of refinement
+        double regionlevel1;                                                                                            // Region which will be refined with level 1
+        double regionlevelmax;                                                                                  // Region which will be refined with level max
+        int step;       //itapopo testando
+        int smooth_frequency; //itapopo testando
+        /*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 *currentmesh;                                                                                // Current Mesh is the refined mesh
+        /*Private methods*/
+   void RefinementProcess(bool &amr_verbose,double* partiallyfloatedelements,double* masklevelset,double* deviatorictensorerror,double* thicknesserror);
+        TPZGeoMesh *previousmesh;                                                                               // Previous Mesh is the refined 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(TPZGeoMesh *gmesh);
+        void GetMesh(int &nvertices,int &nelements,double** px,double** py,int** pelements);
+        void FindElements(int &numberofpoints,double* xp,double* yp,TPZGeoMesh *gmesh,int &hlevel,std::vector<int> &elements);
+   void AllElements(TPZGeoMesh *gmesh,std::vector<int> &elements);
+   inline int GetElemMaterialID(){return 1;}
+   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);
+        TPZGeoMesh* CreateRefPatternMesh(TPZGeoMesh* gmesh);
+   inline int GetElemMaterialID(){return 1;}
+        inline int GetNumberOfNodes(){return 3;}
+        void PrintGMeshVTK(TPZGeoMesh *gmesh,std::ofstream &file,bool matColor=true);
+        int GetVTK_ElType(TPZGeoEl* gel);
+        /*}}}*/
 };

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

r21930	r22100
137	137
138	138	/verify if the metric will be reseted or not/
139		if(~~!this->keepmetric~~){
	139	if(this->keepmetric==0){
140	140	if(this->options->metric) xDelete<IssmDouble>(this->options->metric);
141	141	this->options->metricSize[0] = 0;

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

-              r22004
+              r22100
    xDelete<IssmDouble>(elementlevelset);
    delete vmasklevelset;
+}
 /*}}}*/
 …
    delete vxc;
    delete vyc;
+}
+/*}}}*/
+void FemModel::GetZeroLevelSetPoints(IssmDouble** pzerolevelset_points,int &numberofpoints,int levelset_type){/*{{{*/
+        /*Here, "zero level set" means grounding line or ice front, depending on the level set type*/
+        /*pzerolevelset_points are the element center points with zero level set. X and Y coords*/
+        if(levelset_type!=MaskGroundediceLevelsetEnum && levelset_type!=MaskIceLevelsetEnum){
+                _error_("level set type not implemented yet!");
+        }
+        /*Outputs*/
+        IssmDouble* zerolevelset_points                 = NULL;
+        int npoints                                                                             = 0;
+        /*Intermediaries*/
+        int elementswidth                       = this->GetElementsWidth();
+   int numberofelements                 = this->elements->NumberOfElements();
+        int* elem_vertices                                      = xNew<int>(elementswidth);
+   IssmDouble* levelset                                         = xNew<IssmDouble>(elementswidth);
+   IssmDouble* x                                                                        = NULL;
+   IssmDouble* y                                                                        = NULL;
+   IssmDouble* z                                                                        = NULL;
+        Vector<IssmDouble>* vx_zerolevelset             = new Vector<IssmDouble>(numberofelements);
+        Vector<IssmDouble>* vy_zerolevelset             = new Vector<IssmDouble>(numberofelements);
+        IssmDouble* x_zerolevelset                                      = NULL;
+        IssmDouble* y_zerolevelset                                      = NULL;
+        int count,sid;
+        IssmDouble xcenter,ycenter;
+        /*Get vertices coordinates*/
+        VertexCoordinatesx(&x,&y,&z,this->vertices,false) ;
+        /*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();
+                xcenter         = NAN;
+                ycenter         = 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.;
+                        }
+                }
+                vx_zerolevelset->SetValue(sid,xcenter,INS_VAL);
+                vy_zerolevelset->SetValue(sid,ycenter,INS_VAL);
+        }
+   /*Assemble and serialize*/
+   vx_zerolevelset->Assemble();
+   vy_zerolevelset->Assemble();
+   x_zerolevelset=vx_zerolevelset->ToMPISerial();
+   y_zerolevelset=vy_zerolevelset->ToMPISerial();
+        /*Find the number of points*/
+        npoints=0;
+        for(int i=0;i<numberofelements;i++) if(!xIsNan<IssmDouble>(x_zerolevelset[i])) npoints++;
+        /*Keep just the element center coordinates with zero level set (compact the structure)*/
+        zerolevelset_points=xNew<IssmDouble>(2*npoints);//x and y
+        count=0;
+        for(int i=0;i<numberofelements;i++){
+                if(!xIsNan<IssmDouble>(x_zerolevelset[i])){
+                        zerolevelset_points[2*count]     = x_zerolevelset[i];
+                        zerolevelset_points[2*count+1] = y_zerolevelset[i];
+                        count++;
+                }
+        }
+        /*Assign outputs*/
+        numberofpoints                          = npoints;
+        (*pzerolevelset_points) = zerolevelset_points;
+        /*Cleanup*/
+   xDelete<int>(elem_vertices);
+   xDelete<IssmDouble>(levelset);
+        xDelete<IssmDouble>(x_zerolevelset);
+        xDelete<IssmDouble>(y_zerolevelset);
+   xDelete<IssmDouble>(x);
+   xDelete<IssmDouble>(y);
+   xDelete<IssmDouble>(z);
+        delete vx_zerolevelset;
+        delete vy_zerolevelset;
+}
 /*}}}*/
 …
         /*Intermediaries*/
+   int numberofvertices       = this->vertices->NumberOfVertices();
+   IssmDouble* levelset_points= NULL;
+   IssmDouble* x                                        = NULL;
+   IssmDouble* y                                        = NULL;
+   IssmDouble* z                                        = NULL;
+        int numberofpoints;
+        IssmDouble distance;
+        /*Get vertices coordinates*/
+        VertexCoordinatesx(&x,&y,&z,this->vertices,false) ;
+        /*Get points which level set is zero (center of elements with zero level set)*/
+        this->GetZeroLevelSetPoints(&levelset_points,numberofpoints,levelset_type);
+        /*Find the minimal vertice distance to the zero levelset (grounding line or ice front)*/
+        verticedistance=xNew<IssmDouble>(numberofvertices);
+        for(int i=0;i<numberofvertices;i++){
+                verticedistance[i]=INFINITY;
+                for(int j=0;j<numberofpoints;j++){
+                        distance=sqrt((x[i]-levelset_points[2*j])*(x[i]-levelset_points[2*j])+(y[i]-levelset_points[2*j+1])*(y[i]-levelset_points[2*j+1]));
+                        verticedistance[i]=min(distance,verticedistance[i]);
+                }
+        }
+        /*Assign the pointer*/
+        (*pverticedistance)=verticedistance;
+        /*Cleanup*/
+   xDelete<IssmDouble>(levelset_points);
+   xDelete<IssmDouble>(x);
+   xDelete<IssmDouble>(y);
+   xDelete<IssmDouble>(z);
+}
+/*}}}*/
+#endif
+#if defined(_HAVE_NEOPZ_) && !defined(_HAVE_ADOLC_)
+void FemModel::ReMeshNeopz(int* pnewnumberofvertices,int* pnewnumberofelements,IssmDouble** pnewx,IssmDouble** pnewy,IssmDouble** pnewz,int** pnewelementslist){/*{{{*/
+        /*pnewelementslist keep vertices in Matlab indexing*/
+   int my_rank                                          = IssmComm::GetRank();
+   int numberofelements                 = this->elements->NumberOfElements();
+        IssmDouble* element_label       = xNewZeroInit<IssmDouble>(numberofelements);
+        int numberofpoints                      = -1;
+        IssmDouble* xylist                      = NULL;
+        IssmDouble* newx                                = NULL;
+   IssmDouble* newy                             = NULL;
+   IssmDouble* newz                             = NULL;
+   int* newelementslist                 = NULL;
+   int newnumberofvertices              = -1;
+        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);
+        if(my_rank==0){
+                this->amr->ExecuteRefinement(numberofpoints,xylist,&newnumberofvertices,&newnumberofelements,&newx,&newy,&newelementslist);
+      newz=xNewZeroInit<IssmDouble>(newnumberofvertices);
+                if(newnumberofvertices<=0 || newnumberofelements<=0) _error_("Error in the ReMeshNeopz.");
+        }
+        else{
+                newx=xNew<IssmDouble>(newnumberofvertices);
+                newy=xNew<IssmDouble>(newnumberofvertices);
+                newz=xNew<IssmDouble>(newnumberofvertices);
+                newelementslist=xNew<int>(newnumberofelements*this->GetElementsWidth());
+        }
+        /*Send new mesh to others CPU*/
+        ISSM_MPI_Bcast(&newnumberofvertices,1,ISSM_MPI_INT,0,IssmComm::GetComm());
+        ISSM_MPI_Bcast(&newnumberofelements,1,ISSM_MPI_INT,0,IssmComm::GetComm());
+        ISSM_MPI_Bcast(newx,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+        ISSM_MPI_Bcast(newy,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+        ISSM_MPI_Bcast(newz,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+        ISSM_MPI_Bcast(newelementslist,newnumberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm());
+        /*Assign the pointers*/
+        (*pnewelementslist)     = newelementslist; //Matlab indexing
+        (*pnewx)                                        = newx;
+        (*pnewy)                                        = newy;
+        (*pnewz)                                        = newz;
+        *pnewnumberofvertices= newnumberofvertices;
+        *pnewnumberofelements= newnumberofelements;
+        /*Cleanup*/
+        xDelete<IssmDouble>(element_label);
+        xDelete<IssmDouble>(xylist);
+}
+/*}}}*/
+void FemModel::InitializeAdaptiveRefinementNeopz(void){/*{{{*/
+        /*Define variables*/
+        int my_rank                                                                             = IssmComm::GetRank();
+        int numberofvertices                                                    = this->vertices->NumberOfVertices();
+        int numberofelements                                                    = this->elements->NumberOfElements();
+        IssmDouble* x                                                                   = NULL;
+        IssmDouble* y                                                                   = NULL;
+        IssmDouble* z                                                                   = NULL;
+        int* elements                                                                   = NULL;
+        int level_max                                                                   = -1;
+        IssmDouble radius_level_max                             = -1;
+        IssmDouble gradation                                                    = -1;
+        IssmDouble lag                                                                  = -1;
+        IssmDouble groundingline_distance               = -1;
+        IssmDouble icefront_distance                            = -1;
+   IssmDouble thicknesserror_threshold          = -1;
+        IssmDouble deviatoricerror_threshold    = -1;
+        /*Initialize field as NULL for now*/
+        this->amr = NULL;
+        /*Get vertices coordinates of the coarse mesh (father mesh)*/
+        /*elements comes in Matlab indexing*/
+        this->GetMesh(this->vertices,this->elements,&x,&y,&z,&elements);
+        /*Get amr parameters*/
+        this->parameters->FindParam(&level_max,AmrLevelMaxEnum);
+        this->parameters->FindParam(&radius_level_max,AmrRadiusLevelMaxEnum);
+        this->parameters->FindParam(&gradation,AmrGradationEnum);
+        this->parameters->FindParam(&lag,AmrLagEnum);
+        this->parameters->FindParam(&groundingline_distance,AmrGroundingLineDistanceEnum);
+        this->parameters->FindParam(&icefront_distance,AmrIceFrontDistanceEnum);
+        this->parameters->FindParam(&thicknesserror_threshold,AmrThicknessErrorThresholdEnum);
+        this->parameters->FindParam(&deviatoricerror_threshold,AmrDeviatoricErrorThresholdEnum);
+        /*Create initial mesh (coarse mesh) in neopz data structure*/
+        /*Just CPU #0 should keep AMR object*/
+   /*Initialize refinement pattern*/
+        this->SetRefPatterns();
+        this->amr = new AdaptiveMeshRefinement();
+        this->amr->level_max                                                    = level_max;
+        this->amr->radius_level_max                             = radius_level_max;
+        this->amr->gradation                                                    = gradation;
+        this->amr->lag                                                                  = lag;
+        this->amr->groundingline_distance               = groundingline_distance;
+        this->amr->icefront_distance                            = icefront_distance;
+        this->amr->thicknesserror_threshold             = thicknesserror_threshold;
+        this->amr->deviatoricerror_threshold    = deviatoricerror_threshold;
+        if(my_rank==0){
+                this->amr->CreateInitialMesh(numberofvertices,numberofelements,x,y,elements);
+        }
+        /*Free the vectors*/
+        xDelete<IssmDouble>(x);
+        xDelete<IssmDouble>(y);
+        xDelete<IssmDouble>(z);
+        xDelete<int>(elements);
+}
+/*}}}*/
+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 numberofvertices                 = this->vertices->NumberOfVertices();
    int numberofelements                 = this->elements->NumberOfElements();
         int* elem_vertices                                      = xNew<int>(elementswidth);
    IssmDouble* levelset                                         = xNew<IssmDouble>(elementswidth);
+   IssmDouble* xp                                                                       = NULL;
+   IssmDouble* yp                                                                       = NULL;
+   IssmDouble* x                                                                        = NULL;
+   IssmDouble* y                                                                        = NULL;
+   IssmDouble* z                                                                        = NULL;
+        Vector<IssmDouble>* vx_zerolevelset             = new Vector<IssmDouble>(numberofelements);
+        Vector<IssmDouble>* vy_zerolevelset             = new Vector<IssmDouble>(numberofelements);
+        IssmDouble* x_zerolevelset                                      = NULL;
+        IssmDouble* y_zerolevelset                                      = NULL;
+        int count,sid,npoints;
+        IssmDouble xcenter,ycenter,distance;
+        /*Get vertices coordinates*/
+        VertexCoordinatesx(&x,&y,&z,this->vertices,false) ;
+        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->GetInputListOnVertices(levelset,levelset_type);
                 element->GetVerticesSidList(elem_vertices);
+                sid                     = element->Sid();
+                xcenter         = NAN;
+                ycenter         = NAN;
+                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) {
+                                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.;
+                                label=1.;
+                        }
+                }
+                vx_zerolevelset->SetValue(sid,xcenter,INS_VAL);
+                vy_zerolevelset->SetValue(sid,ycenter,INS_VAL);
+        }
+   /*Assemble and serialize*/
+   vx_zerolevelset->Assemble();
+   vy_zerolevelset->Assemble();
+   x_zerolevelset=vx_zerolevelset->ToMPISerial();
+   y_zerolevelset=vy_zerolevelset->ToMPISerial();
+        /*keep just the element center coordinates with zero level set (compact the structure to save time in the next step)*/
+        count = 0;
+        xp      = xNewZeroInit<IssmDouble>(numberofelements);
+        yp      = xNewZeroInit<IssmDouble>(numberofelements);
+                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){/*{{{*/
+        /*Here, "zero level set" means grounding line or ice front, depending on the level set type*/
+        /*pverticedistance is the minimal vertice distance to the grounding line or ice front*/
+        if(levelset_type!=MaskGroundediceLevelsetEnum && levelset_type!=MaskIceLevelsetEnum){
+                _error_("level set type not implemented yet!");
+        }
+        /*Output*/
+        IssmDouble* elementdistance;
+        /*Intermediaries*/
+   int numberofelements       = this->elements->NumberOfElements();
+   IssmDouble* levelset_points= NULL;
+   IssmDouble* xc                                       = NULL;
+   IssmDouble* yc                                       = NULL;
+        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)*/
+        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++){
+                if(!xIsNan<IssmDouble>(x_zerolevelset[i])){
+                        xp[count]=x_zerolevelset[i];
+                        yp[count]=y_zerolevelset[i];
+                        count++;
+                }
+        }
+        npoints=count;
+        /*Find the minimal vertice distance to the zero levelset (grounding line or ice front)*/
+        verticedistance=xNew<IssmDouble>(numberofvertices);
+        for(int i=0;i<numberofvertices;i++){
+                verticedistance[i]=INFINITY;
+                for(int j=0;j<npoints;j++){
+                        distance=sqrt((x[i]-xp[j])*(x[i]-xp[j])+(y[i]-yp[j])*(y[i]-yp[j]));
+                        verticedistance[i]=min(distance,verticedistance[i]);
+                elementdistance[i]=INFINITY;
+                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]);
+                }
+        }
         /*Assign the pointer*/
         (*pverticedistance)=verticedistance;
+        (*pelementdistance)=elementdistance;
         /*Cleanup*/
+   xDelete<int>(elem_vertices);
+   xDelete<IssmDouble>(levelset);
+        xDelete<IssmDouble>(x_zerolevelset);
+        xDelete<IssmDouble>(y_zerolevelset);
+   xDelete<IssmDouble>(xp);
+   xDelete<IssmDouble>(yp);
+   xDelete<IssmDouble>(x);
+   xDelete<IssmDouble>(y);
+   xDelete<IssmDouble>(z);
+        delete vx_zerolevelset;
+        delete vy_zerolevelset;
+}
+/*}}}*/
+#endif
+#if defined(_HAVE_NEOPZ_) && !defined(_HAVE_ADOLC_)
+void FemModel::ReMeshNeopz(int* pnewnumberofvertices,int* pnewnumberofelements,IssmDouble** pnewx,IssmDouble** pnewy,IssmDouble** pnewz,int** pnewelementslist){/*{{{*/
+        /*pnewelementslist keep vertices in Matlab indexing*/
+   int my_rank                         = IssmComm::GetRank();
+   bool amr_verbose                    = true; //itapopo
+   IssmDouble* x                       = NULL;
+   IssmDouble* y                       = NULL;
+   IssmDouble* z                       = NULL;
+   int* elementslist                   = NULL;
+   int oldnumberofelements             = this->elements->NumberOfElements();
+   int newnumberofvertices                              = -1;
+        int newnumberofelements                                 = -1;
+        IssmDouble* partiallyfloatedelements= NULL;//itapopo verify if it will be used
+   IssmDouble* masklevelset            = NULL;
+   IssmDouble* deviatorictensorerror   = NULL;
+   IssmDouble* thicknesserror          = NULL;
+        /*Get the elements in which grounding line goes through*/
+   //itapopo verificar se irá usar esse this->GetPartiallyFloatedElements(&partiallyfloatedelements);
+   /*Get mean mask level set over each element*/
+   this->MeanGroundedIceLevelSet(&masklevelset);
+   /*Get the deviatoric error estimator*/
+   this->ZZErrorEstimator(&deviatorictensorerror);
+   /*Get the thickness error estimator*/
+   this->ThicknessZZErrorEstimator(&thicknesserror);
+        if(my_rank==0){
+                this->amr->Execute(amr_verbose,oldnumberofelements,partiallyfloatedelements,masklevelset,deviatorictensorerror,thicknesserror,
+                           newnumberofvertices,newnumberofelements,&x,&y,&elementslist);
+      z=xNewZeroInit<IssmDouble>(newnumberofvertices);
+                if(newnumberofvertices<=0 || newnumberofelements<=0) _error_("Error in the ReMeshNeopz.");
+        }
+        else{
+                x=xNew<IssmDouble>(newnumberofvertices);
+                y=xNew<IssmDouble>(newnumberofvertices);
+                z=xNew<IssmDouble>(newnumberofvertices);
+                elementslist=xNew<int>(newnumberofelements*this->GetElementsWidth());
+        }
+        /*Send new mesh to others CPU*/
+        ISSM_MPI_Bcast(&newnumberofvertices,1,ISSM_MPI_INT,0,IssmComm::GetComm());
+        ISSM_MPI_Bcast(&newnumberofelements,1,ISSM_MPI_INT,0,IssmComm::GetComm());
+        ISSM_MPI_Bcast(x,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+        ISSM_MPI_Bcast(y,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+        ISSM_MPI_Bcast(z,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+        ISSM_MPI_Bcast(elementslist,newnumberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm());
+        /*Assign the pointers*/
+        (*pnewelementslist)     = elementslist; //Matlab indexing
+        (*pnewx)                                        = x;
+        (*pnewy)                                        = y;
+        (*pnewz)                                        = z;
+        *pnewnumberofvertices= newnumberofvertices;
+        *pnewnumberofelements= newnumberofelements;
+        /*Cleanup*/
+        xDelete<IssmDouble>(masklevelset);
+        xDelete<IssmDouble>(deviatorictensorerror);
+   xDelete<IssmDouble>(thicknesserror);
+}
+/*}}}*/
+void FemModel::InitializeAdaptiveRefinementNeopz(void){/*{{{*/
+        /*Define variables*/
+        int my_rank                                             = IssmComm::GetRank();
+        int numberofvertices                    = this->vertices->NumberOfVertices();
+        int numberofelements                    = this->elements->NumberOfElements();
+        IssmDouble* x                                   = NULL;
+        IssmDouble* y                                   = NULL;
+        IssmDouble* z                                   = NULL;
+        int* elements                                   = NULL;
+        int elementswidth                               = this->GetElementsWidth(); //just tria elements in this version. Itapopo:
+        int levelmax                                    = 0;
+        IssmDouble regionlevel1         = 0.;
+        IssmDouble regionlevelmax       = 0.;
+        /*Initialize field as NULL for now*/
+        this->amr = NULL;
+        /*Get vertices coordinates of the coarse mesh (father mesh)*/
+        /*elements comes in Matlab indexing*/
+        this->GetMesh(this->vertices,this->elements,&x,&y,&z,&elements);
+        /*Get amr parameters*/
+        this->parameters->FindParam(&levelmax,AmrLevelMaxEnum);
+        this->parameters->FindParam(&regionlevel1,AmrRegionLevel1Enum);
+        this->parameters->FindParam(&regionlevelmax,AmrRegionLevelMaxEnum);
+        /*Create initial mesh (coarse mesh) in neopz data structure*/
+        /*Just CPU #0 should keep AMR object*/
+   /*Initialize refinement pattern*/
+        this->SetRefPatterns();
+        if(my_rank==0){
+                this->amr = new AdaptiveMeshRefinement();
+                this->amr->CreateInitialMesh(numberofvertices,numberofelements,elementswidth,x,y,elements);
+                this->amr->SetLevelMax(levelmax); //Set max level of refinement
+                this->amr->SetRegions(regionlevel1,regionlevelmax);
+        }
+        /*Free the vectors*/
+        xDelete<IssmDouble>(x);
+        xDelete<IssmDouble>(y);
+        xDelete<IssmDouble>(z);
+        xDelete<int>(elements);
+   xDelete<IssmDouble>(levelset_points);
+   xDelete<IssmDouble>(xc);
+   xDelete<IssmDouble>(yc);
+}
 /*}}}*/

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

-              r21919
+              r22100
                 void MeanGroundedIceLevelSet(IssmDouble** pmasklevelset);
                 void GetElementCenterCoordinates(IssmDouble** pxc,IssmDouble** pyc);
+                void GetZeroLevelSetPoints(IssmDouble** pzerolevelset_points,int &numberofpoints,int levelset_type);
                 #endif
 …
                 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);
                 #endif

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

-              r22004
+              r22100
                         case AmrNeopzEnum:
                                 parameters->AddObject(iomodel->CopyConstantObject("md.amr.level_max",AmrLevelMaxEnum));
+                                parameters->AddObject(iomodel->CopyConstantObject("md.amr.region_level_1",AmrRegionLevel1Enum));
+                                parameters->AddObject(iomodel->CopyConstantObject("md.amr.region_level_max",AmrRegionLevelMaxEnum));
+                                parameters->AddObject(iomodel->CopyConstantObject("md.amr.radius_level_max",AmrRadiusLevelMaxEnum));
+                                parameters->AddObject(iomodel->CopyConstantObject("md.amr.gradation",AmrGradationEnum));
+                                parameters->AddObject(iomodel->CopyConstantObject("md.amr.lag",AmrLagEnum));
+                                parameters->AddObject(iomodel->CopyConstantObject("md.amr.groundingline_distance",AmrGroundingLineDistanceEnum));
+                                parameters->AddObject(iomodel->CopyConstantObject("md.amr.icefront_distance",AmrIceFrontDistanceEnum));
+                                parameters->AddObject(iomodel->CopyConstantObject("md.amr.thicknesserror_threshold",AmrThicknessErrorThresholdEnum));
+                                parameters->AddObject(iomodel->CopyConstantObject("md.amr.deviatoricerror_threshold",AmrDeviatoricErrorThresholdEnum));
                                 break;
                         #endif

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

-              r22075
+              r22100
         AmrNeopzEnum,
         AmrLevelMaxEnum,
         AmrRegionLevel1Enum,
         AmrRegionLevelMaxEnum,
+        AmrLagEnum,
+        AmrRadiusLevelMaxEnum,
         AmrBamgEnum,
         AmrHminEnum,

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

-              r22075
+              r22100
                 case AmrNeopzEnum : return "AmrNeopz";
                 case AmrLevelMaxEnum : return "AmrLevelMax";
                 case AmrRegionLevel1Enum : return "AmrRegionLevel1";
                 case AmrRegionLevelMaxEnum : return "AmrRegionLevelMax";
+                case AmrLagEnum : return "AmrLag";
+                case AmrRadiusLevelMaxEnum : return "AmrRadiusLevelMax";
                 case AmrBamgEnum : return "AmrBamg";
                 case AmrHminEnum : return "AmrHmin";

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

-              r22075
+              r22100
               else if (strcmp(name,"AmrNeopz")==0) return AmrNeopzEnum;
               else if (strcmp(name,"AmrLevelMax")==0) return AmrLevelMaxEnum;
               else if (strcmp(name,"AmrRegionLevel1")==0) return AmrRegionLevel1Enum;
               else if (strcmp(name,"AmrRegionLevelMax")==0) return AmrRegionLevelMaxEnum;
+              else if (strcmp(name,"AmrLag")==0) return AmrLagEnum;
+              else if (strcmp(name,"AmrRadiusLevelMax")==0) return AmrRadiusLevelMaxEnum;
               else if (strcmp(name,"AmrBamg")==0) return AmrBamgEnum;
               else if (strcmp(name,"AmrHmin")==0) return AmrHminEnum;

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

-              r21674
+              r22100
         //methods
         this.setdefaultparameters = function(){// {{{
+                //level_max: 2 to 4
+                this.level_max=2;
+                //region_level_1: region around (m) the discontinuity (grounding line or ice front) where the mesh will be refined once (h=1).
+                this.region_level_1=20000.;
+                //region_level_max: region around (m) the discontinuity (grounding line or ice front) where the mesh will be refined with max level of refinement (h=level_max).
+                this.region_level_max=15000.;
+        this.hmin                                                       = 100.;
+      this.hmax                                                         = 100.e3;
+                this.fieldname                                          = "Vel";
+                this.err                                                                = 3.;
+        this.keepmetric                         = 1;
+        this.gradation                                  = 1.5;
+                this.groundingline_resolution           = 500.;
+        this.groundingline_distance     = 0;
+        this.icefront_resolution                = 500;
+        this.icefront_distance          = 0;
+      this.thicknesserror_resolution    = 500;
+      this.thicknesserror_threshold     = 0;
+      this.deviatoricerror_resolution   = 500;
+      this.deviatoricerror_threshold    = 0;
         }// }}}
         this.disp= function(){// {{{
                 console.log(sprintf('   amr parameters:'));
+                fielddisplay(this,'level_max','maximum refinement level (1, 2, 3 or 4)');
+                fielddisplay(this,'region_level_1','region which will be refined once (level 1) [ m ]');
+                fielddisplay(this,'region_level_max','region which will be refined with level_max [ m ]');
+                fielddisplay(this,'hmin','minimum element length');
+      fielddisplay(this,'hmax','maximum element length');
+      fielddisplay(this,'fieldname','name of input that will be used to compute the metric (should be an input of FemModel)');
+                fielddisplay(this,'keepmetric','indicates whether the metric should be kept every remeshing time');
+                fielddisplay(this,'gradation','maximum ratio between two adjacent edges');
+                fielddisplay(this,'groundingline_resolution','element length near the grounding line');
+                fielddisplay(this,'groundingline_distance','distance around the grounding line which elements will be refined');
+                fielddisplay(this,'icefront_resolution','element length near the ice front');
+                fielddisplay(this,'icefront_distance','distance around the ice front which elements will be refined');
+                fielddisplay(this,'thicknesserror_resolution','element length when thickness error estimator is used');
+                fielddisplay(this,'thicknesserror_threshold','maximum threshold 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');
         }// }}}
         this.classname= function(){// {{{
 …
         }// }}}
                 this.checkconsistency = function(md,solution,analyses) { //{{{
+                        checkfield(md,'fieldname','amr.level_max','numel',[1],'>=',0,'<=',4);
+                        checkfield(md,'fieldname','amr.region_level_1','numel',[1],'>',0,'NaN',1,'Inf',1);
+                        checkfield(md,'fieldname','amr.region_level_max','numel',[1],'>',0,'NaN',1,'Inf',1);
+                        if (this.region_level_1-this.region_level_max<0.2*this.region_level_1){
+                                md.checkmessage('region_level_max should be lower than 80% of region_level_1');
+                        }
+         checkfield(md,'fieldname','amr.hmax','numel',[1],'>',0,'NaN',1);
+         checkfield(md,'fieldname','amr.hmin','numel',[1],'>',0,'<',this.hmax,'NaN',1);
+         checkfield(md,'fieldname','amr.keepmetric','numel',[1],'>=',0,'<=',1,'NaN',1);
+         checkfield(md,'fieldname','amr.gradation','numel',[1],'>=',1.1,'<=',5,'NaN',1);
+         checkfield(md,'fieldname','amr.groundingline_resolution','numel',[1],'>',0,'<',this.hmax,'NaN',1);
+         checkfield(md,'fieldname','amr.groundingline_distance','numel',[1],'>=',0,'NaN',1,'Inf',1);
+         checkfield(md,'fieldname','amr.icefront_resolution','numel',[1],'>',0,'<',this.hmax,'NaN',1);
+         checkfield(md,'fieldname','amr.icefront_distance','numel',[1],'>=',0,'NaN',1,'Inf',1);
+         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.deviatoricerror_resolution','numel',[1],'>',0,'<',this.hmax,'NaN',1);
+         checkfield(md,'fieldname','amr.deviatoricerror_threshold','numel',[1],'>=',0,'<=',1,'NaN',1);
                 } // }}}
                 this.marshall=function(md,prefix,fid) { //{{{
+                        WriteData(fid,prefix,'object',this,'fieldname','level_max','format','Integer');
+                        WriteData(fid,prefix,'object',this,'fieldname','region_level_1','format','Double');
+                        WriteData(fid,prefix,'object',this,'fieldname','region_level_max','format','Double');
+         WriteData(fid,prefix,'name','md.amr.type','data',1,'format','Integer');
+         WriteData(fid,prefix,'object',this,'fieldname','hmin','format','Double');
+         WriteData(fid,prefix,'object',this,'fieldname','hmax','format','Double');
+         WriteData(fid,prefix,'object',this,'fieldname','fieldname','format','String');
+         WriteData(fid,prefix,'object',this,'fieldname','err','format','Double');
+         WriteData(fid,prefix,'object',this,'fieldname','keepmetric','format','Integer');
+         WriteData(fid,prefix,'object',this,'fieldname','gradation','format','Double');
+         WriteData(fid,prefix,'object',this,'fieldname','groundingline_resolution','format','Double');
+         WriteData(fid,prefix,'object',this,'fieldname','groundingline_distance','format','Double');
+         WriteData(fid,prefix,'object',this,'fieldname','icefront_resolution','format','Double');
+         WriteData(fid,prefix,'object',this,'fieldname','icefront_distance','format','Double');
+         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','deviatoricerror_resolution','format','Double');
+         WriteData(fid,prefix,'object',this,'fieldname','deviatoricerror_threshold','format','Double');
                 }//}}}
                 this.fix=function() { //{{{
 …
         //properties
         // {{{
+        this.level_max                          = 0;
+        this.region_level_1     = 0.;
+        this.region_level_max   = 0.;
+        this.hmin                                                               = 0.;
+        this.hmax                                                               = 0.;
+        this.fieldname                                                  = "";
+        this.err                                                                        = 0.;
+        this.keepmetric                                         = 0;
+        this.gradation                                                  = 0.;
+        this.groundingline_resolution           = 0.;
+        this.groundingline_distance             = 0.;
+        this.icefront_resolution                        = 0.;
+        this.icefront_distance                          = 0.;
+        this.thicknesserror_resolution  = 0.;
+        this.thicknesserror_threshold           = 0.;
+        this.deviatoricerror_resolution = 0.;
+        this.deviatoricerror_threshold  = 0.;
         this.setdefaultparameters();

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

-              r21802
+              r22100
+%
 %   Usage:
 %      amr=amr();
+%      md.amr=amr();
 classdef amr
         properties (SetAccess=public)
+                level_max                       = 0;
+                region_level_1          = 0;
+                region_level_max        = 0;
+                hmin = 0.;
+                hmax = 0.;
+                fieldname = '';
+                err = 0.;
+                keepmetric = 0;
+                gradation = 0.;
+                groundingline_resolution = 0.;
+                groundingline_distance = 0.;
+                icefront_resolution = 0.;
+                icefront_distance = 0.;
+                thicknesserror_resolution = 0.;
+                thicknesserror_threshold = 0.;
+                deviatoricerror_resolution = 0.;
+                deviatoricerror_threshold = 0.;
+        end
+        methods (Static)
+                function self = loadobj(self) % {{{
+         % This function is directly called by matlab when a model object is
+         % loaded. Update old properties here
+         if verLessThan('matlab','7.9'),
+            disp('Warning: your matlab version is old and there is a risk that load does not work correctly');
+            disp('         if the model is not loaded correctly, rename temporarily loadobj so that matlab does not use it');
+            % This is a Matlab bug: all the fields of md have their default value
+            % Example of error message:
+            % Warning: Error loading an object of class 'model':
+            % Undefined function or method 'exist' for input arguments of type 'cell'
+            %
+            % This has been fixed in MATLAB 7.9 (R2009b) and later versions
+         end
+         %2017 September 15th
+         if isstruct(self),
+            disp('WARNING: updating amr. Now the default is amr with bamg');
+            disp('         some old fields were not converted');
+            disp('         see the new fields typing md.amr and modify them properly');
+            obj2 = self;
+            self = amr();
+            %Converting region_level_max to groundingline_distance
+            if(obj2.region_level_max>0 && obj2.region_level_1>obj2.region_level_max)
+               self.groundingline_distance      = obj2.region_level_max;
+                                        self.keepmetric                                 = 0;
+                                        self.fieldname                                          = 'None';
+            end
+         end
+      end% }}}
         end
         methods
 …
                 function self = setdefaultparameters(self) % {{{
+                        %level_max: 2 to 4
+                        self.level_max=2;
+                        %hmin and hmax
+                        self.hmin=100.;
+                        self.hmax=100.e3;
+                        %region_level_1: region around (m) the discontinuity (grounding line or ice front) where the mesh will be refined once (h=1).
+                        self.region_level_1=20000.;
+                        %region_level_max: region around (m) the discontinuity (grounding line or ice front) where the mesh will be refined with max level of refinement (h=level_max).
+                        self.region_level_max=15000.;
+                        %fields
+                        self.fieldname ='Vel';
+                        self.err=3.;
+                        %keep metric?
+                        self.keepmetric=1;
+                        %control of element lengths
+                        self.gradation=1.5;
+                        %other criterias
+                        self.groundingline_resolution=500.;
+                        self.groundingline_distance=0.;
+                        self.icefront_resolution=500.;
+                        self.icefront_distance=0.;
+                        self.thicknesserror_resolution=500.;
+                        self.thicknesserror_threshold=0.;
+                        self.deviatoricerror_resolution=500.;
+                        self.deviatoricerror_threshold=0.;
                 end % }}}
                 function md = checkconsistency(self,md,solution,analyses) % {{{
+                        md = checkfield(md,'fieldname','amr.level_max','numel',[1],'>=',0,'<=',4);
+                        md = checkfield(md,'fieldname','amr.region_level_1','numel',[1],'>',0,'NaN',1,'Inf',1);
+                        md = checkfield(md,'fieldname','amr.region_level_max','numel',[1],'>',0,'NaN',1,'Inf',1);
+                        if self.region_level_1-self.region_level_max<0.2*self.region_level_1, %it was adopted 20% of the region_level_1
+                                md = checkmessage(md,'region_level_max should be lower than 80% of region_level_1');
+                        end
+                        md = checkfield(md,'fieldname','amr.hmax','numel',[1],'>',0,'NaN',1);
+                        md = checkfield(md,'fieldname','amr.hmin','numel',[1],'>',0,'<',self.hmax,'NaN',1);
+                        %md = checkfield(md,'fieldname','amr.fieldname','string',[1]);
+                        md = checkfield(md,'fieldname','amr.keepmetric','numel',[1],'>=',0,'<=',1,'NaN',1);
+                        md = checkfield(md,'fieldname','amr.gradation','numel',[1],'>=',1.1,'<=',5,'NaN',1);
+                        md = checkfield(md,'fieldname','amr.groundingline_resolution','numel',[1],'>',0,'<',self.hmax,'NaN',1);
+                        md = checkfield(md,'fieldname','amr.groundingline_distance','numel',[1],'>=',0,'NaN',1,'Inf',1);
+                        md = checkfield(md,'fieldname','amr.icefront_resolution','numel',[1],'>',0,'<',self.hmax,'NaN',1);
+                        md = checkfield(md,'fieldname','amr.icefront_distance','numel',[1],'>=',0,'NaN',1,'Inf',1);
+                        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.deviatoricerror_resolution','numel',[1],'>',0,'<',self.hmax,'NaN',1);
+                        md = checkfield(md,'fieldname','amr.deviatoricerror_threshold','numel',[1],'>=',0,'<=',1,'NaN',1);
                 end % }}}
                 function disp(self) % {{{
                         disp(sprintf('   amr parameters:'));
+                        fielddisplay(self,'level_max',['maximum refinement level (1, 2, 3 or 4)']);
+                        fielddisplay(self,'region_level_1',['region which will be refined once (level 1) [ m ]']);
+                        fielddisplay(self,'region_level_max',['region which will be refined with level_max [ m ]']);
+                        fielddisplay(self,'hmin',['minimum element length']);
+                        fielddisplay(self,'hmax',['maximum element length']);
+                        fielddisplay(self,'fieldname',['name of input that will be used to compute the metric (should be an input of FemModel)']);
+                        fielddisplay(self,'keepmetric',['indicates whether the metric should be kept every remeshing time']);
+                        fielddisplay(self,'gradation',['maximum ratio between two adjacent edges']);
+                        fielddisplay(self,'groundingline_resolution',['element length near the grounding line']);
+                        fielddisplay(self,'groundingline_distance',['distance around the grounding line which elements will be refined']);
+                        fielddisplay(self,'icefront_resolution',['element length near the ice front']);
+                        fielddisplay(self,'icefront_distance',['distance around the ice front which elements will be refined']);
+                        fielddisplay(self,'thicknesserror_resolution',['element length when thickness error estimator is used']);
+                        fielddisplay(self,'thicknesserror_threshold',['maximum threshold 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']);
                 end % }}}
                 function marshall(self,prefix,md,fid) % {{{
+                        WriteData(fid,prefix,'name','md.amr.type','data',2,'format','Integer');
+                        WriteData(fid,prefix,'object',self,'fieldname','level_max','format','Integer');
+                        WriteData(fid,prefix,'object',self,'fieldname','region_level_1','format','Double');
+                        WriteData(fid,prefix,'object',self,'fieldname','region_level_max','format','Double');
+                end % }}}
+                function savemodeljs(self,fid,modelname) % {{{
+                        writejsdouble(fid,[modelname '.amr.level_max'],self.level_max);
+                        writejsdouble(fid,[modelname '.amr.region_level_1'],self.region_level_1);
+                        writejsdouble(fid,[modelname '.amr.region_level_max'],self.region_level_max);
+                        WriteData(fid,prefix,'name','md.amr.type','data',1,'format','Integer');
+                        WriteData(fid,prefix,'object',self,'class','amr','fieldname','hmin','format','Double');
+                        WriteData(fid,prefix,'object',self,'class','amr','fieldname','hmax','format','Double');
+                        WriteData(fid,prefix,'object',self,'class','amr','fieldname','fieldname','format','String');
+                        WriteData(fid,prefix,'object',self,'class','amr','fieldname','err','format','Double');
+                        WriteData(fid,prefix,'object',self,'class','amr','fieldname','keepmetric','format','Integer');
+                        WriteData(fid,prefix,'object',self,'class','amr','fieldname','gradation','format','Double');
+                        WriteData(fid,prefix,'object',self,'class','amr','fieldname','groundingline_resolution','format','Double');
+                        WriteData(fid,prefix,'object',self,'class','amr','fieldname','groundingline_distance','format','Double');
+                        WriteData(fid,prefix,'object',self,'class','amr','fieldname','icefront_resolution','format','Double');
+                        WriteData(fid,prefix,'object',self,'class','amr','fieldname','icefront_distance','format','Double');
+                        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','deviatoricerror_resolution','format','Double');
+                        WriteData(fid,prefix,'object',self,'class','amr','fieldname','deviatoricerror_threshold','format','Double');
                 end % }}}

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

-              r21676
+              r22100
     def __init__(self): # {{{
+        self.level_max        = 0.
+        self.region_level_1   = 0.
+        self.region_level_max = 0.
+        self.hmin                                                               = 0.
+        self.hmax                                                               = 0.
+        self.fieldname                                          =''
+        self.err                                                                = 0.
+        self.keepmetric                                         = 0.
+        self.gradation                                          = 0.
+        self.groundingline_resolution   = 0.
+        self.groundingline_distance     = 0.
+        self.icefront_resolution                = 0.
+        self.icefront_distance                  = 0.
+        self.thicknesserror_resolution = 0.
+        self.thicknesserror_threshold   = 0.
+        self.deviatoricerror_resolution= 0.
+        self.deviatoricerror_threshold = 0.
         #set defaults
         self.setdefaultparameters()
 …
     def __repr__(self): # {{{
         string="   amr parameters:"
+        string="%s\n%s"%(string,fielddisplay(self,"level_max","maximum refinement level (1, 2, 3 or 4)"))
+        string="%s\n%s"%(string,fielddisplay(self,"region_level_1","region which will be refined once (level 1) [ m ]"))
+        string="%s\n%s"%(string,fielddisplay(self,"region_level_max","region which will be refined with level_max [ m ]"))
+        string="%s\n%s"%(string,fielddisplay(self,"hmin","minimum element length"))
+        string="%s\n%s"%(string,fielddisplay(self,"hmax","maximum element length"))
+        string="%s\n%s"%(string,fielddisplay(self,"fieldname","name of input that will be used to compute the metric (should be an input of FemModel)"))
+        string="%s\n%s"%(string,fielddisplay(self,"keepmetric","indicates whether the metric should be kept every remeshing time"))
+        string="%s\n%s"%(string,fielddisplay(self,"gradation","maximum ratio between two adjacent edges"))
+        string="%s\n%s"%(string,fielddisplay(self,"groundingline_resolution","element length near the grounding line"))
+        string="%s\n%s"%(string,fielddisplay(self,"groundingline_distance","distance around the grounding line which elements will be refined"))
+        string="%s\n%s"%(string,fielddisplay(self,"icefront_resolution","element length near the ice front"))
+        string="%s\n%s"%(string,fielddisplay(self,"icefront_distance","distance around the ice front which elements will be refined"))
+        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,"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"))
         return string
     #}}}
     def setdefaultparameters(self): # {{{
+        #level_max: 2 to 4
+        self.level_max=2
+        #region_level_1: region around (m) the discontinuity (grounding line or ice front) where the mesh will be refined once (h=1).
+        self.region_level_1=20000.
+        #region_level_max: region around (m) the discontinuity (grounding line or ice front) where the mesh will be refined with max level of refinement (h=level_max).
+        self.region_level_max=15000.
+        self.hmin                                                               = 100.
+        self.hmax                                                               = 100.e3
+        self.fieldname                                          = 'Vel'
+        self.err                                                                = 3.
+        self.keepmetric                                         = 1
+        self.gradation                                          = 1.5
+        self.groundingline_resolution   = 500.
+        self.groundingline_distance     = 0
+        self.icefront_resolution                = 500.
+        self.icefront_distance                  = 0
+        self.thicknesserror_resolution = 500.
+        self.thicknesserror_threshold   = 0
+        self.deviatoricerror_resolution= 500.
+        self.deviatoricerror_threshold = 0
         return self
         #}}}
+    #}}}
     def checkconsistency(self,md,solution,analyses):    # {{{
+        md = checkfield(md,'fieldname','amr.level_max','numel',[1],'>=',0,'<=',4)
+        md = checkfield(md,'fieldname','amr.region_level_1','numel',[1],'>',0,'NaN',1,'Inf',1)
+        md = checkfield(md,'fieldname','amr.region_level_max','numel',[1],'>',0,'NaN',1,'Inf',1)
+                #it was adopted 20% of the region_level_1
+        if self.region_level_1-self.region_level_max<0.2*self.region_level_1:
+            md.checkmessage("region_level_max should be lower than 80% of region_level_1")
+        md = checkfield(md,'fieldname','amr.hmax','numel',[1],'>',0,'NaN',1)
+        md = checkfield(md,'fieldname','amr.hmin','numel',[1],'>',0,'<',self.hmax,'NaN',1)
+        md = checkfield(md,'fieldname','amr.keepmetric','numel',[1],'>=',0,'<=',1,'NaN',1);
+        md = checkfield(md,'fieldname','amr.gradation','numel',[1],'>=',1.1,'<=',5,'NaN',1);
+        md = checkfield(md,'fieldname','amr.groundingline_resolution','numel',[1],'>',0,'<',self.hmax,'NaN',1);
+        md = checkfield(md,'fieldname','amr.groundingline_distance','numel',[1],'>=',0,'NaN',1,'Inf',1);
+        md = checkfield(md,'fieldname','amr.icefront_resolution','numel',[1],'>',0,'<',self.hmax,'NaN',1);
+        md = checkfield(md,'fieldname','amr.icefront_distance','numel',[1],'>=',0,'NaN',1,'Inf',1);
+        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.deviatoricerror_resolution','numel',[1],'>',0,'<',self.hmax,'NaN',1);
+        md = checkfield(md,'fieldname','amr.deviatoricerror_threshold','numel',[1],'>=',0,'<=',1,'NaN',1);
         return md
     # }}}
     def marshall(self,prefix,md,fid):    # {{{
+        WriteData(fid,prefix,'object',self,'fieldname','level_max','format','Integer')
+        WriteData(fid,prefix,'object',self,'fieldname','region_level_1','format','Double')
+        WriteData(fid,prefix,'object',self,'fieldname','region_level_max','format','Double')
+        WriteData(fid,prefix,'name','md.amr.type','data',1,'format','Integer')
+        WriteData(fid,prefix,'object',self,'fieldname','hmin','format','Double');
+        WriteData(fid,prefix,'object',self,'fieldname','hmax','format','Double');
+        WriteData(fid,prefix,'object',self,'fieldname','fieldname','format','String');
+        WriteData(fid,prefix,'object',self,'fieldname','err','format','Double');
+        WriteData(fid,prefix,'object',self,'fieldname','keepmetric','format','Integer');
+        WriteData(fid,prefix,'object',self,'fieldname','gradation','format','Double');
+        WriteData(fid,prefix,'object',self,'fieldname','groundingline_resolution','format','Double');
+        WriteData(fid,prefix,'object',self,'fieldname','groundingline_distance','format','Double');
+        WriteData(fid,prefix,'object',self,'fieldname','icefront_resolution','format','Double');
+        WriteData(fid,prefix,'object',self,'fieldname','icefront_distance','format','Double');
+        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','deviatoricerror_resolution','format','Double');
+        WriteData(fid,prefix,'object',self,'fieldname','deviatoricerror_threshold','format','Double');
     # }}}

issm/trunk-jpl/test/NightlyRun/test462.m

r21955	r22100
11	11	md.transient.isgroundingline=0;
12	12	%amr bamg settings, just field
13		~~md.amr=amrbamg();~~
14	13	md.amr.hmin=10000;
15	14	md.amr.hmax=100000;

issm/trunk-jpl/test/NightlyRun/test463.m

r21956	r22100
11	11	md.transient.isgroundingline=0;
12	12	%amr bamg settings, just grounding line
13		~~md.amr=amrbamg();~~
14	13	md.amr.hmin=10000;
15	14	md.amr.hmax=100000;

issm/trunk-jpl/test/NightlyRun/test464.m

r21955	r22100
11	11	md.transient.isgroundingline=0;
12	12	%amr bamg settings, just ice front
13		~~md.amr=amrbamg();~~
14	13	md.amr.hmin=10000;
15	14	md.amr.hmax=100000;

issm/trunk-jpl/test/NightlyRun/test465.m

r21956	r22100
11	11	md.transient.isgroundingline=0;
12	12	%amr bamg settings, field, grounding line and ice front
13		~~md.amr=amrbamg();~~
14	13	md.amr.hmin=20000;
15	14	md.amr.hmax=100000;

Context Navigation

Legend:

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

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

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

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

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

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

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

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

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

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

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

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

issm/trunk-jpl/test/NightlyRun/test462.m

issm/trunk-jpl/test/NightlyRun/test463.m

issm/trunk-jpl/test/NightlyRun/test464.m

issm/trunk-jpl/test/NightlyRun/test465.m

Download in other formats: