Index: /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp	(revision 21805)
+++ /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp	(revision 21806)
@@ -10,4 +10,9 @@
 
 #include "./AdaptiveMeshRefinement.h"
+#include "TPZVTKGeoMesh.h"
+#include "../shared/shared.h"
+#include "pzgeotriangle.h"
+#include "pzreftriangle.h"
+using namespace pzgeom;
 
 /*Constructor, copy, clean up and destructor*/
@@ -17,5 +22,5 @@
 /*}}}*/
 AdaptiveMeshRefinement::AdaptiveMeshRefinement(const AdaptiveMeshRefinement &cp){/*{{{*/
-   this->Initialize(); 
+	this->Initialize(); 
 	this->operator =(cp);
 }
@@ -25,12 +30,15 @@
 	/*Clean all attributes*/
 	this->CleanUp();
-
 	/*Copy all data*/
-	this->fathermesh     = new TPZGeoMesh(*cp.fathermesh);
-	this->previousmesh   = new TPZGeoMesh(*cp.previousmesh); 
-	this->levelmax       = cp.levelmax;
-	this->elementswidth  = cp.elementswidth;
-	this->regionlevel1   = cp.regionlevel1;
-	this->regionlevelmax = cp.regionlevelmax;
+	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->sid2index.clear();
+	this->sid2index.resize(cp.sid2index.size());
+	for(int i=0;i<cp.sid2index.size();i++) this->sid2index[i]=cp.sid2index[i];
+	
 	return *this;
 
@@ -39,16 +47,16 @@
 AdaptiveMeshRefinement::~AdaptiveMeshRefinement(){/*{{{*/
 	
-	//bool ismismip = false;
-	//if(ismismip){//itapopo
-	//	TPZFileStream fstr;
-	//	std::stringstream ss;
+	bool ismismip = false;
+	if(ismismip){//itapopo
+		TPZFileStream fstr;
+		std::stringstream ss;
 	    
-	//	ss << this->levelmax;
-	//	std::string AMRfile	= "/home/santos/Misomip2/L" + ss.str() + "_tsai/amr.txt"; 
-	
-	//	fstr.OpenWrite(AMRfile.c_str());
-	//	int withclassid = 1;
-	//	this->Write(fstr,withclassid);
-	//}
+		ss << this->levelmax;
+		std::string AMRfile	= "/home/santos/L" + ss.str() + "_amr.txt"; 
+	
+		fstr.OpenWrite(AMRfile.c_str());
+		int withclassid = 1;
+		this->Write(fstr,withclassid);
+	}
 	this->CleanUp();
 	gRefDBase.clear();
@@ -57,11 +65,12 @@
 void AdaptiveMeshRefinement::CleanUp(){/*{{{*/
 
-    /*Verify and delete all data*/
+	/*Verify and delete all data*/
 	if(this->fathermesh)    delete this->fathermesh;
-   if(this->previousmesh)  delete this->previousmesh;
-	this->levelmax			= -1;
-	this->elementswidth  = -1;
-	this->regionlevel1	= -1;
-	this->regionlevelmax = -1;
+	if(this->currentmesh)   delete this->currentmesh;
+	this->levelmax				= -1;
+	this->elementswidth		= -1;
+	this->regionlevel1		= -1;
+	this->regionlevelmax		= -1;
+	this->sid2index.clear();
 }
 /*}}}*/
@@ -69,10 +78,11 @@
 
 	/*Set pointers to NULL*/
-	this->fathermesh		= NULL;
-	this->previousmesh	= NULL;
-	this->levelmax			= -1;
-	this->elementswidth	= -1;
-	this->regionlevel1	= -1;
-	this->regionlevelmax = -1;
+	this->fathermesh			= NULL;
+	this->currentmesh			= NULL;
+	this->levelmax				= -1;
+	this->elementswidth		= -1;
+	this->regionlevel1		= -1;
+	this->regionlevelmax		= -1;
+	this->sid2index.clear();
 }
 /*}}}*/
@@ -81,420 +91,339 @@
 }
 /*}}}*/
-void AdaptiveMeshRefinement::Read(TPZStream &buf, void *context){/*{{{*/
-
-    try
-    {
-        /* Read the id context*/
-        TPZSaveable::Read(buf,context);
-
-        /* Read class id*/
-        int classid;
-        buf.Read(&classid,1);
-        
-        /* Verify the class id*/
-        if (classid != this->ClassId() )
-        {
-            std::cout << "Error in restoring AdaptiveMeshRefinement!\n";
-            std::cout.flush();
-            DebugStop();
-        }
-        
-        /* Read simple attributes */
-        buf.Read(&this->levelmax,1);
-        buf.Read(&this->elementswidth,1);
-        buf.Read(&this->regionlevel1,1);
-        buf.Read(&this->regionlevelmax,1);
-        
-		/* Read geometric mesh*/
-        TPZSaveable *sv1 = TPZSaveable::Restore(buf,0);
-        this->fathermesh = dynamic_cast<TPZGeoMesh*>(sv1);
-        
-        TPZSaveable *sv2 = TPZSaveable::Restore(buf,0);
-        this->previousmesh = dynamic_cast<TPZGeoMesh*>(sv2);
+void AdaptiveMeshRefinement::Read(TPZStream &buf,void *context){/*{{{*/
+
+	try
+	{
+		/* Read the id context*/
+		TPZSaveable::Read(buf,context);
+		/* Read class id*/
+		int classid;
+		buf.Read(&classid,1); 
+		/* Verify the class id*/
+      if(classid!=this->ClassId()) _error_("AdaptiveMeshRefinement::Read: Error in restoring AdaptiveMeshRefinement!\n"); 
+		/* Read simple attributes */
+		buf.Read(&this->levelmax,1);
+		buf.Read(&this->elementswidth,1);
+		buf.Read(&this->regionlevel1,1);
+		buf.Read(&this->regionlevelmax,1);
+		/* Read vector attributes*/
+		int size;
+		buf.Read(&size,1);
+		int* psid2index=xNew<int>(size);
+		buf.Read(psid2index,size);
+		this->sid2index.clear();
+		this->sid2index.assign(psid2index,psid2index+size);
+		/* Read geometric mesh (father)*/
+		TPZSaveable *sv1 = TPZSaveable::Restore(buf,0);
+		this->fathermesh = dynamic_cast<TPZGeoMesh*>(sv1);
+		/* Read geometric mesh (current)*/
+		TPZSaveable *sv2 = TPZSaveable::Restore(buf,0);
+		this->currentmesh = dynamic_cast<TPZGeoMesh*>(sv2);
+		/* Cleanup*/
+		xDelete<int>(psid2index);
+	}
+	catch(const std::exception& e)
+	{
+		_error_("AdaptiveMeshRefinement::Read: Exception catched!\n");
+	}
+
+}
+/*}}}*/
+template class TPZRestoreClass<AdaptiveMeshRefinement,13829430>;/*{{{*/
+/*}}}*/
+void AdaptiveMeshRefinement::Write(TPZStream &buf,int withclassid){/*{{{*/
+    
+	try
+	{
+		/* Write context (this class) class ID*/
+		TPZSaveable::Write(buf,withclassid);
+		/* Write this class id*/
+		int classid = this->ClassId();
+		buf.Write(&classid,1);
+		/* Write simple attributes */
+		buf.Write(&this->levelmax,1);
+		buf.Write(&this->elementswidth,1);
+		buf.Write(&this->regionlevel1,1);
+		buf.Write(&this->regionlevelmax,1);
+		/* Write vector attributes*/
+		int size=this->sid2index.size();
+		int* psid2index=&this->sid2index[0];
+		buf.Write(&size,1);//vector size
+		buf.Write(psid2index,this->sid2index.size());
+		/* Write the geometric mesh*/
+		this->fathermesh->Write(buf,this->ClassId());
+		this->currentmesh->Write(buf,this->ClassId());
     }
     catch(const std::exception& e)
     {
-        std::cout << "Exception catched! " << e.what() << std::endl;
-        std::cout.flush();
-        DebugStop();
+		_error_("AdaptiveMeshRefinement::Write: Exception catched!\n");
     }
 }
 /*}}}*/
-template class TPZRestoreClass<AdaptiveMeshRefinement,13829430>;/*{{{*/
-/*}}}*/
-void AdaptiveMeshRefinement::Write(TPZStream &buf, int withclassid){/*{{{*/
-    
-    try
-    {
-        /* Write context (this class) class ID*/
-        TPZSaveable::Write(buf,withclassid);
-
-        /* Write this class id*/
-        int classid = ClassId();
-        buf.Write(&classid,1);
-
-        /* Write simple attributes */
-        buf.Write(&this->levelmax,1);
-        buf.Write(&this->elementswidth,1);
-        buf.Write(&this->regionlevel1,1);
-        buf.Write(&this->regionlevelmax,1);
-			
-        /* Write the geometric mesh*/
-        this->fathermesh->Write(buf, this->ClassId());
-        this->previousmesh->Write(buf, this->ClassId());
-    }
-    catch(const std::exception& e)
-    {
-        std::cout << "Exception catched! " << e.what() << std::endl;
-        std::cout.flush();
-        DebugStop();
-    }
-}
-/*}}}*/
 
 /*Mesh refinement methods*/
-#include "TPZVTKGeoMesh.h" //itapopo
-#include "../shared/shared.h" //itapopo
-void AdaptiveMeshRefinement::ExecuteRefinement(int &type_process,double *vx, double *vy, double *masklevelset, int &nvertices, int &nelements, int &nsegments, double** px, double** py, double** pz, int** pelements, int** psegments){/*{{{*/
-
-	/*IMPORTANT! pelements (and psegments) are in Matlab indexing*/
+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*/
 	/*NEOPZ works only in C indexing*/
-
-    _assert_(this->fathermesh);
-    _assert_(this->previousmesh);
-    
-    /*Calculate the position of the grounding line using previous mesh*/
-    std::vector<TPZVec<REAL> > GLvec;
-    this->CalcGroundingLinePosition(masklevelset, GLvec);
-    
-   // std::ofstream file1("/home/santos/mesh0.vtk");
-   // TPZVTKGeoMesh::PrintGMeshVTK(this->fathermesh,file1 );
-    
-    /*run refinement or unrefinement process*/
-    TPZGeoMesh *newmesh;
-    switch (type_process) {
-        case 0: newmesh = this->previousmesh; break;                    // refine previous mesh
-        case 1: newmesh = new TPZGeoMesh(*this->fathermesh); break;     // refine mesh 0 (unrefine process)
-        default: DebugStop(); break;//itapopo verificar se irá usar _assert_
-    }
-    
-    this->RefinementProcess(newmesh,GLvec);
-	
-    //std::ofstream file2("/home/santos/mesh1.vtk");
-    //TPZVTKGeoMesh::PrintGMeshVTK(this->previousmesh,file2 );
-    
-    /*Set new mesh pointer. Previous mesh just have uniform elements*/
-    if(type_process==1){
-        if(this->previousmesh) delete this->previousmesh;
-        this->previousmesh = newmesh;
-    }
-    
-    /*Refine elements to avoid hanging nodes*/
-	//TPZGeoMesh *nohangingnodesmesh = new TPZGeoMesh(*newmesh);//itapopo testando, este era o original
-   TPZGeoMesh *nohangingnodesmesh = this->CreateRefPatternMesh(newmesh);//itapopo testando, este eh novo metodo
-    
-    //std::ofstream file3("/home/santos/mesh2.vtk");
-    //TPZVTKGeoMesh::PrintGMeshVTK(this->previousmesh,file3);
-    
-    this->RefineMeshToAvoidHangingNodes(nohangingnodesmesh);
-    
-	 //std::ofstream file4("/home/santos/mesh3.vtk");
-    //TPZVTKGeoMesh::PrintGMeshVTK(nohangingnodesmesh,file4);
-    
-    /*Get new geometric mesh in ISSM data structure*/
-    this->GetMesh(nohangingnodesmesh,nvertices,nelements,nsegments,px,py,pz,pelements,psegments);
-	 
-    /*Verify the new geometry*/
-    this->CheckMesh(nvertices,nelements,nsegments,this->elementswidth,px,py,pz,pelements,psegments);
-
-	 _printf_("\trefinement process done!\n\n");
-
-    delete nohangingnodesmesh;
-}
-/*}}}*/
-void AdaptiveMeshRefinement::RefinementProcess(TPZGeoMesh *gmesh,std::vector<TPZVec<REAL> > &GLvec){/*{{{*/
-    
-    /*Refine mesh levelmax times*/
-   _printf_("\n\trefinement process (level max = " << this->levelmax << ")\n");
-	_printf_("\tprogress:  ");
-	for(int hlevel=1;hlevel<=this->levelmax;hlevel++){
-        
-        /*Set elements to be refined using some criteria*/
-        std::vector<int> ElemVec; //elements without children
-        this->SetElementsToRefine(gmesh,GLvec,hlevel,ElemVec);
-        
-        /*Refine the mesh*/
-        this->RefineMesh(gmesh, ElemVec);
-		  
-		  _printf_("*  ");
-    }
-    _printf_("\n");
-}
-/*}}}*/
-void AdaptiveMeshRefinement::RefineMesh(TPZGeoMesh *gmesh, std::vector<int> &ElemVec){/*{{{*/
-
-	/*Refine elements in ElemVec: uniform pattern refinement*/
-	for(int i = 0; i < ElemVec.size(); i++){
-		
-        /*Get geometric element and verify if it has already been refined*/
-        int index = ElemVec[i];
-        TPZGeoEl * geoel = gmesh->Element(index);
-        if(geoel->HasSubElement()) DebugStop();                              //itapopo _assert_(!geoel->HasSubElement());
-        if(geoel->MaterialId() != this->GetElemMaterialID()) DebugStop();   //itapopo verificar se usará _assert_
-        
-        /*Divide geoel*/
-        TPZVec<TPZGeoEl *> Sons;
-		  geoel->Divide(Sons);
-        
-        /*If a 1D segment is neighbor, it must be divided too*/
-        if(this->elementswidth != 3) DebugStop(); //itapopo verificar o segment para malha 3D
-        
-        std::vector<int> sides(3);
-        sides[0] = 3; sides[1] = 4; sides[2] = 5;
-        for(int j = 0; j < sides.size(); j++ ){
-            
-            TPZGeoElSide Neighbour = geoel->Neighbour(sides[j]);
-            
-            if( Neighbour.Element()->MaterialId() == this->GetBoundaryMaterialID() && !Neighbour.Element()->HasSubElement() ){
-                TPZVec<TPZGeoEl *> pv2;
-                Neighbour.Element()->Divide(pv2);
-            }
-        }
-	}
-    
-    gmesh->BuildConnectivity();
-
+	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");
+
+	/*Execute the refinement.*/
+	this->RefinementProcess(amr_verbose,partiallyfloatedelements,masklevelset,deviatorictensorerror,thicknesserror);
+    
+	/*Get new geometric mesh in ISSM data structure*/
+	this->GetMesh(newnumberofvertices,newnumberofelements,x,y,elementslist);
+	
+	/*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){/*{{{*/
+   
+	if(amr_verbose) _printf_("\n\trefinement process started (level max = " << this->levelmax << ")\n");
+	
+	/*Intermediaries*/
+	TPZGeoMesh* nohangingnodesmesh=NULL;
+	double mean_mask		= 0;
+	double mean_tauerror = 0;
+	double mean_Herror	= 0;
+	double group_error	= 0;
+   
+	/*Calculate mean values*/
+	for(int i=0;i<this->sid2index.size();i++){
+		mean_mask		+= masklevelset[i]; 
+		mean_tauerror	+= deviatorictensorerror[i]; 
+		mean_Herror		+= thicknesserror[i];
+	}
+	mean_mask		/= this->sid2index.size();
+	mean_tauerror	/= this->sid2index.size();
+	mean_Herror		/= this->sid2index.size();
+
+	if(amr_verbose) _printf_("\t\tuniform refinement...\n");
+	for(int i=0;i<this->sid2index.size();i++){
+		TPZGeoEl* geoel=this->currentmesh->Element(this->sid2index[i]);
+		if(geoel->HasSubElement()) _error_("Impossible to refine: geoel (index) "<<this->sid2index[i]<<" has subelements!\n");
+		if(geoel->MaterialId()!=this->GetElemMaterialID()) _error_("Impossible to refine: geoel->MaterialId is not GetElemMaterialID!\n");
+		/*Refine*/
+		if(thicknesserror[i]>mean_Herror){
+			TPZVec<TPZGeoEl *> sons;
+			if(geoel->Level()<this->levelmax) geoel->Divide(sons);
+		} 
+		else if(geoel->Level()>0){ /*try to unrefine*/
+			TPZVec<TPZGeoEl *> sons;
+			geoel->Father()->GetHigherSubElements(sons);
+			group_error=0;
+			for(int j=0;j<sons.size();j++){
+				sons[j]->Index();
+			}
+		}
+	}
+	this->currentmesh->BuildConnectivity();
+	
+	if(amr_verbose) _printf_("\t\trefine to avoid hanging nodes...\n");
+	this->RefineMeshToAvoidHangingNodes(this->currentmesh);
+	this->currentmesh->BuildConnectivity();
+	
+		//nohangingnodesmesh = this->CreateRefPatternMesh(newmesh); itapopo tentar otimizar
+	
+	if(amr_verbose) _printf_("\trefinement process done!\n");
+}
+/*}}}*/
+void AdaptiveMeshRefinement::RefineMesh(TPZGeoMesh *gmesh,std::vector<int> &elements){/*{{{*/
+
+	/*Refine elements: uniform pattern refinement*/
+	for(int i=0;i<elements.size();i++){
+		/*Get geometric element and verify if it has already been refined*/
+		int index = elements[i];
+		TPZGeoEl * geoel = gmesh->Element(index);
+		if(geoel->HasSubElement()) _error_("Impossible to refine: geoel (index) " << index << " has subelements!\n");
+		if(geoel->MaterialId()!=this->GetElemMaterialID()) _error_("Impossible to refine: geoel->MaterialId is not GetElemMaterialID!\n");
+		/*Divide geoel*/
+		TPZVec<TPZGeoEl *> Sons;
+		geoel->Divide(Sons);
+	}
+	gmesh->BuildConnectivity();
 }
 /*}}}*/
 void AdaptiveMeshRefinement::RefineMeshToAvoidHangingNodes(TPZGeoMesh *gmesh){/*{{{*/
    
-	 _printf_("\trefine to avoid hanging nodes...\n");
-    /*Refine elements to avoid hanging nodes: non-uniform refinement*/
-	 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){
-            TPZVec<TPZGeoEl *> Sons;
-            geoel->SetRefPattern(refp);
-            geoel->Divide(Sons);
-        }
-        
-    }
-    
-    gmesh->BuildConnectivity();
-    
-}
-/*}}}*/
-void AdaptiveMeshRefinement::GetMesh(TPZGeoMesh *gmesh,int &nvertices,int &nelements,int &nsegments,double** px,double** py,double** pz, int** pelements, int** psegments){/*{{{*/
-
-	/*IMPORTANT! pelements (and psegments) are in Matlab indexing*/
-	/*NEOPZ works only in C indexing*/
-
-	/* vertices */
-    int ntotalvertices = gmesh->NNodes();//total
-    
-    /* mesh coords */
-	 double* newmeshX = xNew<IssmDouble>(ntotalvertices);
-    double* newmeshY = xNew<IssmDouble>(ntotalvertices);
-    double* newmeshZ = xNew<IssmDouble>(ntotalvertices);
-   	
-   /* getting 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];
-        newmeshZ[i] = coords[2];
-    }
-    
-	/* elements */
-    std::vector<TPZGeoEl*> GeoVec; GeoVec.clear();
-    for(int i = 0; i < gmesh->NElements(); i++){ 
-        if( gmesh->ElementVec()[i]->HasSubElement() ) continue;
-        if( gmesh->ElementVec()[i]->MaterialId() != this->GetElemMaterialID() ) continue;
-        GeoVec.push_back( gmesh->ElementVec()[i]);
-    }
-    
-    int ntotalelements = (int)GeoVec.size();
-    int* newelements = xNew<int>(ntotalelements*this->elementswidth);
-
-    if ( !(this->elementswidth == 3) && !(this->elementswidth == 4) && !(this->elementswidth == 6) ) DebugStop();
-
-    for(int i=0;i<GeoVec.size();i++){
-        for(int j=0;j<this->elementswidth;j++) newelements[i*this->elementswidth+j]=(int)GeoVec[i]->NodeIndex(j)+1;//C to Matlab indexing
-	 }
-    
-    /* segments */
-    std::vector<TPZGeoEl*> SegVec; SegVec.clear();
-    for(int i = 0; i < gmesh->NElements(); i++){
-        if( gmesh->ElementVec()[i]->HasSubElement() ) continue;
-        if( gmesh->ElementVec()[i]->MaterialId() != this->GetBoundaryMaterialID() ) continue;
-        SegVec.push_back( gmesh->ElementVec()[i]);
-    }
-    
-    int ntotalsegments = (int)SegVec.size();
-    int *newsegments=NULL;
-	 if(ntotalsegments>0) newsegments=xNew<int>(ntotalsegments*3);
-    
-    for(int i=0;i<SegVec.size();i++){
-        
-        for(int j=0;j<2;j++) newsegments[i*3+j]=(int)SegVec[i]->NodeIndex(j)+1;//C to Matlab indexing
-        
-        int neighborindex = SegVec[i]->Neighbour(2).Element()->Index();
-        int neighbourid = -1;
-        
-        for(int j = 0; j < GeoVec.size(); j++){
-            if( GeoVec[j]->Index() == neighborindex || GeoVec[j]->FatherIndex() == neighborindex){
-                neighbourid = j;
-                break;
-            }
-        }
-        
-        if(neighbourid==-1) DebugStop(); //itapopo talvez passar para _assert_
-        newsegments[i*3+2] = neighbourid+1;//C to Matlab indexing
-    }
-    
-    //setting outputs
-    nvertices  = ntotalvertices;
-    nelements  = ntotalelements;
-    nsegments  = ntotalsegments;
-    *px		   = newmeshX;
-    *py		   = newmeshY;
-    *pz		   = newmeshZ;
-    *pelements = newelements;
-    *psegments = newsegments;
-    
-}
-/*}}}*/
-void AdaptiveMeshRefinement::CalcGroundingLinePosition(double *masklevelset,std::vector<TPZVec<REAL> > &GLvec){/*{{{*/
-    
-    /* Find grounding line using elments center point */
-    GLvec.clear();
-    for(int i=0;i<this->previousmesh->NElements();i++){
-        
-        if(this->previousmesh->Element(i)->MaterialId()!=this->GetElemMaterialID()) continue;
-        if(this->previousmesh->Element(i)->HasSubElement()) continue;
-        
-        //itapopo apenas malha 2D triangular!
-        int vertex0 = this->previousmesh->Element(i)->NodeIndex(0);
-        int vertex1 = this->previousmesh->Element(i)->NodeIndex(1);
-        int vertex2 = this->previousmesh->Element(i)->NodeIndex(2);
-        
-		  //itapopo inserir uma verificação para não acessar fora da memória
-        double mls0 = masklevelset[vertex0];
-        double mls1 = masklevelset[vertex1];
-        double mls2 = masklevelset[vertex2];
-        
-        if( mls0*mls1 < 0. || mls1*mls2 < 0. ){
-            const int side = 6;
-            TPZVec<double> qsi(2,0.);
-            TPZVec<double> X(3,0.);
-            this->previousmesh->Element(i)->CenterPoint(side, qsi);
-            this->previousmesh->Element(i)->X(qsi, X);
-            GLvec.push_back(X);
-        }
-    }
-    
-//    itapopo apenas para debugar
-//    std::ofstream fileGL("/Users/santos/Desktop/gl.nb");
-//    fileGL << "ListPlot[{";
-//    for(int i = 0; i < GLvec.size(); i++){
-//        fileGL << "{" << GLvec[i][0] << "," << GLvec[i][1] << /*"," << 0. << */"}";
-//        if(i != GLvec.size()-1) fileGL << ",";
-//    }
-//    fileGL << "}]";
-//    fileGL.flush();
-//    fileGL.close();
-    
-}
-/*}}}*/
-void AdaptiveMeshRefinement::SetElementsToRefine(TPZGeoMesh *gmesh,std::vector<TPZVec<REAL> > &GLvec,int &hlevel,std::vector<int> &ElemVec){/*{{{*/
-
-    if(!gmesh) DebugStop(); //itapopo verificar se usará _assert_
-
-    ElemVec.clear();
-    
-    // itapopo inserir modo de encontrar criterio TESTING!!!! Come to false!
-	 if(false) this->TagAllElements(gmesh,ElemVec); //uniform, refine all elements!
-
-    /* Adaptive refinement. This refines some elements following some criteria*/
-    this->TagElementsNearGroundingLine(gmesh, GLvec, hlevel, ElemVec);
-
-}
-/*}}}*/
-void AdaptiveMeshRefinement::TagAllElements(TPZGeoMesh *gmesh,std::vector<int> &ElemVec){/*{{{*/
-    
+	/*Refine elements to avoid hanging nodes: non-uniform refinement*/
+	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){
+			TPZVec<TPZGeoEl *> Sons;
+			geoel->SetRefPattern(refp);
+			geoel->Divide(Sons);
+      }
+	}
+   gmesh->BuildConnectivity();
+    
+}
+/*}}}*/
+void AdaptiveMeshRefinement::GetMesh(int &nvertices,int &nelements,double** px,double** py, int** pelements){/*{{{*/
+
+	/* IMPORTANT! pelements are in Matlab indexing
+	   NEOPZ works only in C indexing.
+		This method cleans up and updated the this->sid2index 
+		and fills in it with the new mesh.
+		Avoid to call this method before Refinement Process.*/
+
+	/*Intermediaries */
+	TPZGeoMesh* gmesh = this->currentmesh;//itapopo confirmar
+	
+	long sid,nodeindex;
+	int nconformelements,nconformvertices;	
+	int* newelements			= NULL;
+	double* newmeshX			= NULL;//xNew<double>(ntotalvertices);
+	double* newmeshY			= NULL;//xNew<double>(ntotalvertices);
+	TPZGeoEl* geoel			= NULL;
+	long* vertex_index2sid 	= xNew<long>(gmesh->NNodes());
+	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;
+	
+	/*Get elements without sons and fill in the vertex_index2sid with used vertices (indexes) */
+	sid=0;
+	for(int i=0;i<gmesh->NElements();i++){//over gmesh elements index 
+		geoel=gmesh->ElementVec()[i];
+		if(!geoel) continue;
+		if(geoel->HasSubElement()) continue;
+		if(geoel->MaterialId() != this->GetElemMaterialID()) continue;
+		this->sid2index.push_back(i);//keep the element index
+		for(int j=0;j<this->elementswidth;j++){
+      	nodeindex=geoel->NodeIndex(j);
+      	if(vertex_index2sid[nodeindex]==-1){
+      		vertex_index2sid[nodeindex]=sid; 
+				sid++;
+			}
+      }	
+	}
+
+	nconformelements	= (int)this->sid2index.size();
+	nconformvertices	= (int)sid;
+	newelements			= xNew<int>(nconformelements*this->elementswidth);
+	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<this->sid2index.size();i++){//over the sid (fill the ISSM elements)
+		for(int j=0;j<this->elementswidth;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
+		}
+	}
+ 
+	/*Setting outputs*/
+	nvertices	= nconformvertices;
+	nelements	= nconformelements;
+	*px			= newmeshX;
+	*py		   = newmeshY;
+	*pelements	= newelements;
+   
+	/*Cleanup*/
+	xDelete<long>(vertex_index2sid);
+
+}
+/*}}}*/
+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;
-        ElemVec.push_back(i);
+        elements.push_back(i);
     }
 }
 /*}}}*/
-void AdaptiveMeshRefinement::TagElementsNearGroundingLine(TPZGeoMesh *gmesh,std::vector<TPZVec<REAL> > &GLvec,int &hlevel,std::vector<int> &ElemVec){/*{{{*/
-    
-    /* Tag elements near grounding line */ 
-	 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));
-    
-    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;
-        
-        const int side2D = 6;
-        TPZVec<REAL> qsi(2,0.);
-        TPZVec<REAL> centerPoint(3,0.);
-        gmesh->Element(i)->CenterPoint(side2D, qsi);
-        gmesh->Element(i)->X(qsi, centerPoint);
-        
-        REAL distance = Di;
-        
-        for (int j = 0; j < GLvec.size(); j++) {
-            
-            REAL value = ( GLvec[j][0] - centerPoint[0] ) * ( GLvec[j][0] - centerPoint[0] ); // (x2-x1)^2
-            value += ( GLvec[j][1] - centerPoint[1] ) * ( GLvec[j][1] - centerPoint[1] );// (y2-y1)^2
-            value = std::sqrt(value); //Radius
-            
-            //finding the min distance to the grounding line
-            if(value < distance) distance = value;
-            
-        }
-        
-        if(distance < Di) ElemVec.push_back(i);
-    }
-    
-}
-/*}}}*/
-void AdaptiveMeshRefinement::CreateInitialMesh(int &nvertices,int &nelements,int &nsegments,int &width,double* x,double* y,double* z,int* elements,int* segments){/*{{{*/
-
-	/*IMPORTANT! elements come in Matlab indexing*/
-	/*NEOPZ works only in C indexing*/
-	
-	_assert_(nvertices>0);
-   _assert_(nelements>0);
+void AdaptiveMeshRefinement::CreateInitialMesh(int &nvertices,int &nelements,int &width,double* x,double* y,int* elements){/*{{{*/
+
+	/* IMPORTANT! elements come in Matlab indexing
+		NEOPZ works only in C 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) _error_("Initial mesh already exists!");
+   if(this->fathermesh || this->currentmesh) _error_("Initial mesh already exists!");
     
    this->fathermesh = new TPZGeoMesh();
-	this->fathermesh->NodeVec().Resize( nvertices );
-
-    /*Set the vertices (geometric nodes in NeoPZ context)*/
+	this->fathermesh->NodeVec().Resize(nvertices);
+
+	/*Set the vertices (geometric nodes in NeoPZ context)*/
 	for(int i=0;i<nvertices;i++){  
       /*x,y,z coords*/
@@ -502,5 +431,5 @@
       coord[0]= x[i];
       coord[1]= y[i];
-      coord[2]= z[i];
+      coord[2]= 0.;
       /*Insert in the mesh*/
       this->fathermesh->NodeVec()[i].SetCoord(coord);
@@ -512,46 +441,25 @@
    const int mat = this->GetElemMaterialID();
    TPZManVector<long> elem(this->elementswidth,0);
-    
-	for(int iel=0;iel<nelements;iel++){
-
-		for(int jel=0;jel<this->elementswidth;jel++) elem[jel]=elements[iel*this->elementswidth+jel]-1;//Convert Matlab to C indexing
-
+   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
       /*reftype = 0: uniform, fast / reftype = 1: uniform and non-uniform (avoid hanging nodes), it is not too fast */
-      const int reftype = 0;
+      const int reftype = 1;
       switch(this->elementswidth){
-			case 3: this->fathermesh->CreateGeoElement(ETriangle, elem, mat, index, reftype);	break;
-         case 4: this->fathermesh->CreateGeoElement(ETetraedro, elem, mat, index, reftype); DebugStop(); break;
-			case 6: this->fathermesh->CreateGeoElement(EPrisma, elem, mat, index, reftype); DebugStop(); break;
-         default:	DebugStop();//itapopo _error_("mesh not supported yet");
+			case 3: this->fathermesh->CreateGeoElement(ETriangle,elem,mat,index,reftype);	break;
+         default:	_error_("mesh not supported yet");
 		}
-        
       /*Define the element ID*/        
-      this->fathermesh->ElementVec()[index]->SetId(iel); 
-	}
-    
-   /*Generate the 1D segments elements (boundary)*/
-   const int matboundary = this->GetBoundaryMaterialID();
-   TPZManVector<long> boundary(2,0.);
-    
-   for(int iel=nelements;iel<nelements+nsegments;iel++){     
-		boundary[0] = segments[(iel-nelements)*2+0]-1;//Convert Matlab to C indexing
-      boundary[1] = segments[(iel-nelements)*2+1]-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 = 0;
-      this->fathermesh->CreateGeoElement(EOned, boundary, matboundary, index, reftype);//cria elemento unidimensional
-      this->fathermesh->ElementVec()[index]->SetId(iel);
-	}
-    
+      this->fathermesh->ElementVec()[index]->SetId(i);
+		/*Initialize sid2index*/
+		this->sid2index.push_back((int)index);
+	}
    /*Build element and node connectivities*/
    this->fathermesh->BuildConnectivity();
-    
-	/*Create previous mesh as a copy of father mesh*/
-   this->previousmesh = new TPZGeoMesh(*this->fathermesh);
-
-}
-/*}}}*/
-#include "pzgeotriangle.h" //itapopo
-#include "pzreftriangle.h" //itapopo
-using namespace pzgeom;
+	/*Set current mesh*/
+	this->currentmesh=new TPZGeoMesh(*this->fathermesh);
+}
+/*}}}*/
 TPZGeoMesh* AdaptiveMeshRefinement::CreateRefPatternMesh(TPZGeoMesh* gmesh){/*{{{*/
 	
@@ -646,38 +554,39 @@
 /*}}}*/
 void AdaptiveMeshRefinement::SetElementWidth(int &width){/*{{{*/
-    this->elementswidth = width;
-}
-/*}}}*/
-void AdaptiveMeshRefinement::CheckMesh(int &nvertices,int &nelements,int &nsegments,int &width,double** px,double** py,double** pz,int** pelements, int** psegments){/*{{{*/
-
-    /*Basic verification*/
-    if( !(nvertices > 0) || !(nelements > 0) ) DebugStop(); //itapopo verificar se irá usar o _assert_
-    
-    if ( !(width == 3) && !(width == 4) && !(width == 6) ) DebugStop(); // itapopo verifcar se irá usar o _assert_
-    
-    if( !px || !py || !pz || !pelements ) DebugStop(); // itapopo verifcar se irá usar o _assert_
-    
-    /*Verify if there are orphan nodes*/
-    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 ) DebugStop();//itapopo verificar se irá usar o _assert_
-	
-    //Verify if there are inf or NaN in coords
-    for(int i = 0; i < nvertices; i++) if(isnan((*px)[i]) || isinf((*px)[i])) DebugStop();
-    for(int i = 0; i < nvertices; i++) if(isnan((*py)[i]) || isinf((*py)[i])) DebugStop();
-    for(int i = 0; i < nvertices; i++) if(isnan((*pz)[i]) || isinf((*pz)[i])) DebugStop();
-   
-	 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]) ) DebugStop();
-        }
-    }
-    
-}
-/*}}}*/
+	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){/*{{{*/
+
+	/*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");
+	if(!pelements) _error_("Impossible to continue: pelements is NULL!\n");
+
+	/*Verify if there are orphan nodes*/
+	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");
+	
+	//Verify if there are inf or NaN in coords
+	for(int i=0;i<nvertices;i++){
+		if(isnan((*px)[i]) || isinf((*px)[i])) _error_("Impossible to continue: px i=" << i <<" is NaN or Inf!\n"); 
+		if(isnan((*py)[i]) || isinf((*py)[i])) _error_("Impossible to continue: py i=" << i <<" is NaN or Inf!\n");
+	}
+	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");
+		}
+	}
+    
+}
+/*}}}*/
Index: /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h	(revision 21805)
+++ /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h	(revision 21806)
@@ -45,46 +45,45 @@
 	/*Public methods*/
 	/* Constructor, destructor etc*/
-	AdaptiveMeshRefinement();																// Default constructor
-	AdaptiveMeshRefinement(const AdaptiveMeshRefinement &cp); 					// Copy constructor
-	AdaptiveMeshRefinement & operator= (const AdaptiveMeshRefinement &cp);	// Operator of copy
-	virtual ~AdaptiveMeshRefinement();													// Destructor
+	AdaptiveMeshRefinement();															
+	AdaptiveMeshRefinement(const AdaptiveMeshRefinement &cp); 					
+	AdaptiveMeshRefinement & operator= (const AdaptiveMeshRefinement &cp);	
+	virtual ~AdaptiveMeshRefinement();												
 
     /*Savable methods*/
-	virtual int ClassId() const;                                            // ClassId to save the class
-   virtual void Read(TPZStream &buf, void *context);								// Read this class
-   virtual void Write(TPZStream &buf, int withclassid);                    // Write this class, using ClassId to identify
+	virtual int ClassId() const;                                   
+   virtual void Read(TPZStream &buf,void *context);							
+   virtual void Write(TPZStream &buf,int withclassid);
     
 	/*General methods*/
-	void CleanUp();																			// Clean all attributes
-	void Initialize();																		// Initialize the attributes with NULL and values out of usually range
-	void SetLevelMax(int &h);                                               // Define the max level of refinement
-   void SetRegions(double &D1,double Dhmax);										// Define the regions which will be refined
-	void SetElementWidth(int &width);                                       // Define elements width
-	void ExecuteRefinement(int &type_process,double *vx,double *vy,double *masklevelset,int &nvertices,int &nelements,int &nsegments,double** px,double** py,double** pz,int** pelements,int** psegments=NULL);					// A new mesh will be created and refined. This returns the new mesh
-	void CreateInitialMesh(int &nvertices,int &nelements,int &nsegments,int &width,double* x,double* y,double* z,int* elements,int* segments=NULL); // Create a NeoPZ geometric mesh by coords and elements
+	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 &nsegments,int &width,double** px,double** py,double** pz,int** pelements,int** psegments=NULL); // Check the consistency of the mesh
+	void CheckMesh(int &nvertices,int &nelements,int &width,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
-	TPZGeoMesh *fathermesh;																	// Father Mesh is the entire mesh without refinement
-	TPZGeoMesh *previousmesh;																// Previous mesh is a refined mesh of last step
+   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
+	std::vector<int> sid2index;									// Vector that keeps index of PZGeoMesh elements used in the ISSM mesh (sid) 
+	TPZGeoMesh *fathermesh;											// Father Mesh is the entire mesh without refinement
+	TPZGeoMesh *currentmesh;										// Current Mesh is the refined mesh
 
 	/*Private methods*/
-   void RefinementProcess(TPZGeoMesh *gmesh,std::vector<TPZVec<REAL> > &GLvec);  // Start the refinement process
-	void RefineMesh(TPZGeoMesh *gmesh, std::vector<int> &ElemVec); 					// Refine the elements in ElemVec
-   void RefineMeshToAvoidHangingNodes(TPZGeoMesh *gmesh);                        // Refine the elements to avoid hanging nodes
-	void SetElementsToRefine(TPZGeoMesh *gmesh,std::vector<TPZVec<REAL> > &GLvec,int &hlevel, std::vector<int> &ElemVec); 	//Define wich elements will be refined
-   void TagAllElements(TPZGeoMesh *gmesh,std::vector<int> &ElemVec);				 // This tag all elements to be refined, that is, refine all elements
-   void TagElementsNearGroundingLine(TPZGeoMesh *gmesh,std::vector<TPZVec<REAL> > &GLvec,int &hlevel,std::vector<int> &ElemVec);    // This tag elements near the grounding line
-   void CalcGroundingLinePosition(double *masklevelset,std::vector<TPZVec<REAL> > &GLvec);	// Calculate the grounding line position using previous mesh
-	void GetMesh(TPZGeoMesh *gmesh,int &nvertices,int &nelements,int &nsegments,double** px,double** py,double** pz,int** pelements,int** psegments=NULL); // Return coords and elements in ISSM data structure
-   inline int GetElemMaterialID(){return 1;}                               // Return element material ID
-   inline int GetBoundaryMaterialID(){return 2;}                           // Return segment (2D boundary) material ID
+   void RefinementProcess(bool &amr_verbose,double* partiallyfloatedelements,double* masklevelset,double* deviatorictensorerror,double* thicknesserror);
+	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;}          
 };
 
Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 21805)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 21806)
@@ -2693,4 +2693,10 @@
 void FemModel::WriteMeshInResults(void){/*{{{*/
 
+	//itapopo
+	#ifdef _HAVE_NEOPZ_
+	this->WriteErrorEstimatorsInResults();
+	#endif
+	//itapopo
+
 	int step					= -1;
 	int numberofelements = -1;
@@ -2751,4 +2757,37 @@
 	xDelete<IssmDouble>(z);
 	xDelete<int>(elementslist);
+}
+/*}}}*/
+void FemModel::WriteErrorEstimatorsInResults(void){/*{{{*/
+
+   int step                   = -1;
+   int numberofelements       = -1;
+   IssmDouble time            = -1;
+   IssmDouble* stresserror    = NULL;
+   IssmDouble* thicknesserror = NULL;
+
+   if(!this->elements || !this->vertices || !this->results || !this->parameters) return;
+
+   parameters->FindParam(&step,StepEnum);
+   parameters->FindParam(&time,TimeEnum);
+   numberofelements=this->elements->NumberOfElements();
+
+   /*Compute the deviatoric stress tensor*/
+   this->ZZErrorEstimator(&stresserror);
+
+   /*Compute the thickness error*/
+   this->ThicknessZZErrorEstimator(&thicknesserror);
+
+   /*Write error estimators in Results*/
+   this->results->AddResult(new GenericExternalResult<IssmDouble*>(this->results->Size()+1,DeviatoricStressErrorEstimatorEnum,
+                                                                  stresserror,numberofelements,1,step,time));
+
+   this->results->AddResult(new GenericExternalResult<IssmDouble*>(this->results->Size()+1,ThicknessErrorEstimatorEnum,
+                                                                  thicknesserror,numberofelements,1,step,time));
+   /*Cleanup*/
+   xDelete<IssmDouble>(stresserror);
+   xDelete<IssmDouble>(thicknesserror);
+
+   return;
 }
 /*}}}*/
@@ -3212,9 +3251,6 @@
 void FemModel::SmoothedDeviatoricStressTensor(IssmDouble** ptauxx,IssmDouble** ptauyy,IssmDouble** ptauxy){/*{{{*/
 	
-	this->DeviatoricStressx();//itapopo
-
 	int elementswidth							= this->GetElementsWidth();//just 2D mesh, tria elements
    int numberofvertices						= this->vertices->NumberOfVertices();
-
    IssmDouble weight 						= 0.;
 	IssmDouble*	tauxx							= NULL; 
@@ -3231,5 +3267,9 @@
    Vector<IssmDouble>* vectotalweight	= new Vector<IssmDouble>(numberofvertices);
 	
-   for(int i=0;i<this->elements->Size();i++){
+	/*Update the Deviatoric Stress tensor over the elements*/
+	this->DeviatoricStressx();
+	
+   /*Calculate the Smoothed Deviatoric Stress tensor*/
+	for(int i=0;i<this->elements->Size();i++){
       Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
       element->GetInputListOnVertices(deviatoricstressxx,DeviatoricStressxxEnum);
@@ -3299,5 +3339,5 @@
 void FemModel::ZZErrorEstimator(IssmDouble** pelementerror){/*{{{*/
 
-	/*Compute the Zienkiewicz and Zhu (ZZ) error estimator. 
+	/*Compute the Zienkiewicz and Zhu (ZZ) error estimator for the deviatoric stress tensor. 
 	 * Ref.: Zienkiewicz and Zhu, A Simple Error Estimator and Adaptive Procedure for Practical Engineering Analysis, Int. J. Numer. Meth. Eng, 1987*/
 
@@ -3342,9 +3382,9 @@
 				ftxy+=(tauxy[n]-smoothedtauxy[elem_vertices[n]])*basis[n];
 			}
-			error+=Jdet*gauss->weight*( std::pow(ftxx,2)+std::pow(ftyy,2)+std::pow(ftxy,2) ); 
+			error+=Jdet*gauss->weight*( std::pow(ftxx,2)+std::pow(ftyy,2)+std::pow(ftxy,2) ); //e^2 
 		}
 		/*Set the error in the global vector*/	
       sid=element->Sid();
-		velementerror->SetValue(sid,error,INS_VAL);	
+		velementerror->SetValue(sid,std::sqrt(error),INS_VAL);//sqrt( e^2 )
 		/*Cleanup intermediaries*/
 		xDelete<IssmDouble>(xyz_list);
@@ -3368,4 +3408,179 @@
 	xDelete<int>(elem_vertices);
 	delete velementerror;
+}
+/*}}}*/
+void FemModel::SmoothedGradThickness(IssmDouble** pdHdx,IssmDouble** pdHdy){/*{{{*/
+
+   int elementswidth                   = this->GetElementsWidth();//just 2D mesh, tria elements
+   int numberofvertices                = this->vertices->NumberOfVertices();
+
+   IssmDouble weight                   = 0.;
+   IssmDouble* dHdx                    = NULL;
+   IssmDouble* dHdy                    = NULL;
+   IssmDouble* totalweight             = NULL;
+   IssmDouble* xyz_list                = NULL;
+   IssmDouble* H                       = xNew<IssmDouble>(elementswidth);
+   IssmDouble* GradH                   = xNew<IssmDouble>(2);
+   int* elem_vertices                  = xNew<int>(elementswidth);
+   Vector<IssmDouble>* vecdHdx         = new Vector<IssmDouble>(numberofvertices);
+   Vector<IssmDouble>* vecdHdy         = new Vector<IssmDouble>(numberofvertices);
+   Vector<IssmDouble>* vectotalweight  = new Vector<IssmDouble>(numberofvertices);
+
+   for(int i=0;i<this->elements->Size();i++){
+      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
+      element->GetInputListOnVertices(H,ThicknessEnum);
+      element->GetVerticesSidList(elem_vertices);
+      element->GetVerticesCoordinates(&xyz_list);
+
+      /*Get the gradient of thickness at the center point (in fact, GradH is constante over the element)*/
+      Gauss* gauss=element->NewGauss(1);
+      gauss->GaussPoint(gauss->begin());
+      element->ValueP1DerivativesOnGauss(GradH,H,xyz_list,gauss);
+
+      /*weight to calculate the smoothed grad H*/
+      Tria* triaelement = xDynamicCast<Tria*>(element);
+      weight            = triaelement->GetArea();//the tria area is a choice for the weight
+      
+		/*dH/dx*/
+      vecdHdx->SetValue(elem_vertices[0],weight*GradH[0],ADD_VAL);
+      vecdHdx->SetValue(elem_vertices[1],weight*GradH[0],ADD_VAL);
+      vecdHdx->SetValue(elem_vertices[2],weight*GradH[0],ADD_VAL);
+      /*dH/dy*/
+      vecdHdy->SetValue(elem_vertices[0],weight*GradH[1],ADD_VAL);
+      vecdHdy->SetValue(elem_vertices[1],weight*GradH[1],ADD_VAL);
+      vecdHdy->SetValue(elem_vertices[2],weight*GradH[1],ADD_VAL);
+      /*total weight*/
+      vectotalweight->SetValue(elem_vertices[0],weight,ADD_VAL);
+      vectotalweight->SetValue(elem_vertices[1],weight,ADD_VAL);
+      vectotalweight->SetValue(elem_vertices[2],weight,ADD_VAL);
+      /*Cleanup intermediaries*/
+      xDelete<IssmDouble>(xyz_list);
+      delete gauss;
+   }
+
+   /*Assemble*/
+   vecdHdx->Assemble();
+   vecdHdy->Assemble();
+   vectotalweight->Assemble();
+
+   /*Serialize*/
+   dHdx        = vecdHdx->ToMPISerial();
+   dHdy        = vecdHdy->ToMPISerial();
+   totalweight = vectotalweight->ToMPISerial();
+
+   /*Divide for the total weight*/
+   for(int i=0;i<numberofvertices;i++){
+      _assert_(totalweight[i]>0);
+      dHdx[i] = dHdx[i]/totalweight[i];
+      dHdy[i] = dHdy[i]/totalweight[i];
+   }
+
+   /*Set output*/
+   (*pdHdx) = dHdx;
+   (*pdHdy) = dHdy;
+
+ 	/*Cleanup*/
+   delete vecdHdx;
+   delete vecdHdy;
+   delete vectotalweight;
+   xDelete<IssmDouble>(H);
+   xDelete<IssmDouble>(GradH);
+   xDelete<IssmDouble>(totalweight);
+   xDelete<int>(elem_vertices);
+}
+/*}}}*/
+void FemModel::ThicknessZZErrorEstimator(IssmDouble** pelementerror){/*{{{*/
+   /*Compute the Zienkiewicz and Zhu (ZZ) error estimator for the thickness
+    * Ref.: Zienkiewicz and Zhu, A Simple Error Estimator and Adaptive Procedure for Practical Engineering Analysis, Int. J. Numer. Meth. Eng, 1987*/
+
+   IssmDouble Jdet,error,fdHdx,fdHdy;
+   int sid;
+   int numnodes                     = this->GetElementsWidth();//just 2D mesh, tria elements, P1
+   int numberofelements             = this->elements->NumberOfElements();
+   IssmDouble* xyz_list             = NULL;
+   IssmDouble* smoothed_dHdx        = NULL;
+   IssmDouble* smoothed_dHdy        = NULL;
+   IssmDouble* H                    = xNew<IssmDouble>(numnodes);
+   IssmDouble* GradH                = xNew<IssmDouble>(2);
+   IssmDouble* basis                = xNew<IssmDouble>(numnodes);
+   int* elem_vertices               = xNew<int>(numnodes);
+   Vector<IssmDouble>* velementerror= new Vector<IssmDouble>(numberofelements);
+
+   /*Get smoothed deviatoric stress tensor*/
+   this->SmoothedGradThickness(&smoothed_dHdx,&smoothed_dHdy);
+   
+	/*Integrate the error over elements*/
+   for(int i=0;i<this->elements->Size();i++){
+      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
+      element->GetInputListOnVertices(H,ThicknessEnum);
+      element->GetVerticesSidList(elem_vertices);
+      element->GetVerticesCoordinates(&xyz_list);
+      /*Get the gradient of thickness*/
+      Gauss* gaussH=element->NewGauss(1);
+      gaussH->GaussPoint(gaussH->begin());
+      element->ValueP1DerivativesOnGauss(GradH,H,xyz_list,gaussH);
+      /*Integrate*/
+      Gauss* gauss=element->NewGauss(2);
+      error=0.;
+      for(int ig=gauss->begin();ig<gauss->end();ig++){
+         gauss->GaussPoint(ig);
+         element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+         element->NodalFunctions(basis,gauss);
+         fdHdx=0;fdHdy=0;
+         for(int n=0;n<numnodes;n++) {
+            fdHdx+=(GradH[0]-smoothed_dHdx[elem_vertices[n]])*basis[n];
+            fdHdy+=(GradH[1]-smoothed_dHdy[elem_vertices[n]])*basis[n];
+         }
+         error+=Jdet*gauss->weight*( std::pow(fdHdx,2)+std::pow(fdHdy,2) ); //e^2
+      }
+      /*Set the error in the global vector*/
+      sid=element->Sid();
+      velementerror->SetValue(sid,std::sqrt(error),INS_VAL);//sqrt( e^2 )
+      /*Cleanup intermediaries*/
+      xDelete<IssmDouble>(xyz_list);
+      delete gaussH;
+      delete gauss;
+   }
+
+   /*Assemble*/
+   velementerror->Assemble();
+
+   /*Serialize and set output*/
+   (*pelementerror)=velementerror->ToMPISerial();
+
+   /*Cleanup*/
+   xDelete<IssmDouble>(smoothed_dHdx);
+   xDelete<IssmDouble>(smoothed_dHdy);
+   xDelete<IssmDouble>(H);
+   xDelete<IssmDouble>(GradH);
+   xDelete<IssmDouble>(basis);
+   xDelete<int>(elem_vertices);
+   delete velementerror;
+}
+/*}}}*/
+void FemModel::MeanGroundedIceLevelSet(IssmDouble** pmasklevelset){/*{{{*/
+
+   int elementswidth                   = this->GetElementsWidth();
+   int numberofelements                = this->elements->NumberOfElements();
+   IssmDouble* elementlevelset         = xNew<IssmDouble>(elementswidth);
+   Vector<IssmDouble>* vmasklevelset   = new Vector<IssmDouble>(numberofelements);
+
+   for(int i=0;i<this->elements->Size();i++){
+      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
+      element->GetInputListOnVertices(elementlevelset,MaskGroundediceLevelsetEnum);
+      int sid = element->Sid();
+      vmasklevelset->SetValue(sid,(elementlevelset[0]+elementlevelset[1]+elementlevelset[2])/3.,INS_VAL);
+   }
+
+   /*Assemble*/
+   vmasklevelset->Assemble();
+
+   /*Serialize and set output*/
+   (*pmasklevelset)=vmasklevelset->ToMPISerial();
+
+   /*Cleanup*/
+   xDelete<IssmDouble>(elementlevelset);
+   delete vmasklevelset;
+
 }
 /*}}}*/
@@ -4160,64 +4375,63 @@
 
 #ifdef _HAVE_NEOPZ_
-void FemModel::ReMeshNeopz(int* pnumberofvertices,int* pnumberofelements,IssmDouble** px,IssmDouble** py,IssmDouble** pz,int** pelementslist){/*{{{*/
-
-	/*elements is in Matlab indexing*/
-	int my_rank						= IssmComm::GetRank();
-	int numberofsegments			= -1;
-	int numberofvertices,numberofelements;
-	IssmDouble* vx					= NULL; //itapopo this is not being used
-	IssmDouble* vy					= NULL; //itapopo this is not being used
-	IssmDouble* x					= NULL;
-	IssmDouble* y					= NULL;
-	IssmDouble* z					= NULL;
-	int* elementslist				= NULL;
-	int* segments					= NULL;
-	IssmDouble* masklevelset	= NULL;
-	IssmDouble* pelementerror	= NULL;
-	const int elementswidth		= this->GetElementsWidth();//just 2D mesh, tria elements
-
-	/*Solutions which will be used to refine the elements*/
-	this->GetGroundediceLevelSet(&masklevelset);//itapopo verificar se já existe um método igual a esse
-
-	//Compute the ZZ error estimator per element
-	this->ZZErrorEstimator(&pelementerror);
-
-	_printf0_("P Element error\n");
-	for(int i=0;i<this->elements->NumberOfElements();i++)	_printf0_(""<<pelementerror[i]<< "\n");
-	_printf0_("\n");
+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){
-		int type_process=1; //1: it refines father mesh. See AdaptiveMeshRefinement.h (.cpp)
-		this->amr->ExecuteRefinement(type_process,vx,vy,masklevelset,
-					numberofvertices,numberofelements,numberofsegments,&x,&y,&z,&elementslist,&segments);
-		if(numberofvertices<=0 || numberofelements<=0 /*|| newnumberofsegments<=0*/) _error_("Error in the refinement process.");
+		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>(numberofvertices);
-		y=xNew<IssmDouble>(numberofvertices);
-		z=xNew<IssmDouble>(numberofvertices);
-		elementslist=xNew<int>(numberofelements*this->GetElementsWidth());
+		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(&numberofvertices,1,ISSM_MPI_INT,0,IssmComm::GetComm());
-	ISSM_MPI_Bcast(&numberofelements,1,ISSM_MPI_INT,0,IssmComm::GetComm());
-	ISSM_MPI_Bcast(x,numberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());	
-	ISSM_MPI_Bcast(y,numberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());	
-	ISSM_MPI_Bcast(z,numberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());	
-	ISSM_MPI_Bcast(elementslist,numberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm());	
+	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*/	
-	(*pelementslist) = elementslist; //Matlab indexing
-	(*px)				  = x;
-	(*py)				  = y;
-	(*pz)				  = z;
-	*pnumberofelements = numberofelements;
-	*pnumberofvertices = numberofvertices;
+	(*pnewelementslist) 	= elementslist; //Matlab indexing
+	(*pnewx)				  	= x;
+	(*pnewy)				  	= y;
+	(*pnewz)				  	= z;
+	*pnewnumberofvertices= newnumberofvertices;
+	*pnewnumberofelements= newnumberofelements;
 
 	/*Cleanup*/
-	if(segments) xDelete<int>(segments);
 	xDelete<IssmDouble>(masklevelset);
-	xDelete<IssmDouble>(pelementerror);
+	xDelete<IssmDouble>(deviatorictensorerror);
+   xDelete<IssmDouble>(thicknesserror);
 
 }
@@ -4229,5 +4443,4 @@
 	int numberofvertices			= this->vertices->NumberOfVertices();
 	int numberofelements			= this->elements->NumberOfElements();
-	int numberofsegments			= 0; //used on matlab
 	IssmDouble* x					= NULL;
 	IssmDouble* y					= NULL;
@@ -4282,5 +4495,5 @@
 			//this->amr->SetLevelMax(levelmax); //Set max level of refinement
 			//this->amr->SetRegions(regionlevel1,regionlevelmax);
-			this->amr->CreateInitialMesh(numberofvertices,numberofelements,numberofsegments,elementswidth,x,y,z,elements,NULL);
+			this->amr->CreateInitialMesh(numberofvertices,numberofelements,elementswidth,x,y,elements);
 		}
 		this->amr->SetLevelMax(levelmax); //Set max level of refinement
@@ -4298,8 +4511,6 @@
 
    /*Initialize the global variable of refinement patterns*/
-   gRefDBase.InitializeUniformRefPattern(EOned);
    gRefDBase.InitializeUniformRefPattern(ETriangle);
 
-    //gRefDBase.InitializeRefPatterns();
    /*Insert specifics patterns to ISSM core*/
    std::string filepath  = REFPATTERNDIR;
Index: /issm/trunk-jpl/src/c/classes/FemModel.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.h	(revision 21805)
+++ /issm/trunk-jpl/src/c/classes/FemModel.h	(revision 21806)
@@ -170,7 +170,10 @@
 		void ElementsAndVerticesPartitioning(int& newnumberofvertices,int& newnumberofelements,int& elementswidth,int* newelementslist,bool** pmy_elements,int** pmy_vertices);
 		void WriteMeshInResults(void);
-		void SetRefPatterns(void);
+		void WriteErrorEstimatorsInResults(void);
 		void SmoothedDeviatoricStressTensor(IssmDouble** ptauxx,IssmDouble** ptauyy,IssmDouble** ptauxy); //nodal values, just for SSA-P1: TauXX, TauYY, TauXY
 		void ZZErrorEstimator(IssmDouble** pelementerror);
+		void SmoothedGradThickness(IssmDouble** pdHdx,IssmDouble** pdHdy);
+		void ThicknessZZErrorEstimator(IssmDouble** pelementerror);
+		void MeanGroundedIceLevelSet(IssmDouble** pmasklevelset);
 
 		#ifdef _HAVE_BAMG_
@@ -183,4 +186,5 @@
 		void ReMeshNeopz(int* pnewnumberofvertices,int* pnewnumberofelements,IssmDouble** pnewx,IssmDouble** pnewy,IssmDouble** pnewz,int** pnewelementslist);
 		void InitializeAdaptiveRefinementNeopz(void);
+		void SetRefPatterns(void);
 		#endif
 };
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 21805)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 21806)
@@ -844,4 +844,6 @@
 	AmrFieldEnum,
 	AmrErrEnum,
+	DeviatoricStressErrorEstimatorEnum,
+	ThicknessErrorEstimatorEnum,
 	/*}}}*/
 	ParametersENDEnum,
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 21805)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 21806)
@@ -820,4 +820,6 @@
 		case AmrFieldEnum : return "AmrField";
 		case AmrErrEnum : return "AmrErr";
+		case DeviatoricStressErrorEstimatorEnum : return "DeviatoricStressErrorEstimator";
+		case ThicknessErrorEstimatorEnum : return "ThicknessErrorEstimator";
 		case ParametersENDEnum : return "ParametersEND";
 		case XYEnum : return "XY";
Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 21805)
+++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 21806)
@@ -838,4 +838,6 @@
 	      else if (strcmp(name,"AmrField")==0) return AmrFieldEnum;
 	      else if (strcmp(name,"AmrErr")==0) return AmrErrEnum;
+	      else if (strcmp(name,"DeviatoricStressErrorEstimator")==0) return DeviatoricStressErrorEstimatorEnum;
+	      else if (strcmp(name,"ThicknessErrorEstimator")==0) return ThicknessErrorEstimatorEnum;
 	      else if (strcmp(name,"ParametersEND")==0) return ParametersENDEnum;
 	      else if (strcmp(name,"XY")==0) return XYEnum;
@@ -873,10 +875,10 @@
 	      else if (strcmp(name,"Neumannflux")==0) return NeumannfluxEnum;
 	      else if (strcmp(name,"Param")==0) return ParamEnum;
-	      else if (strcmp(name,"Moulin")==0) return MoulinEnum;
-	      else if (strcmp(name,"Pengrid")==0) return PengridEnum;
          else stage=8;
    }
    if(stage==8){
-	      if (strcmp(name,"Penpair")==0) return PenpairEnum;
+	      if (strcmp(name,"Moulin")==0) return MoulinEnum;
+	      else if (strcmp(name,"Pengrid")==0) return PengridEnum;
+	      else if (strcmp(name,"Penpair")==0) return PenpairEnum;
 	      else if (strcmp(name,"Profiler")==0) return ProfilerEnum;
 	      else if (strcmp(name,"MatrixParam")==0) return MatrixParamEnum;
@@ -996,10 +998,10 @@
 	      else if (strcmp(name,"IceVolumeAboveFloatation")==0) return IceVolumeAboveFloatationEnum;
 	      else if (strcmp(name,"TotalFloatingBmb")==0) return TotalFloatingBmbEnum;
-	      else if (strcmp(name,"TotalGroundedBmb")==0) return TotalGroundedBmbEnum;
-	      else if (strcmp(name,"TotalSmb")==0) return TotalSmbEnum;
          else stage=9;
    }
    if(stage==9){
-	      if (strcmp(name,"P0")==0) return P0Enum;
+	      if (strcmp(name,"TotalGroundedBmb")==0) return TotalGroundedBmbEnum;
+	      else if (strcmp(name,"TotalSmb")==0) return TotalSmbEnum;
+	      else if (strcmp(name,"P0")==0) return P0Enum;
 	      else if (strcmp(name,"P0Array")==0) return P0ArrayEnum;
 	      else if (strcmp(name,"P1")==0) return P1Enum;
Index: /issm/trunk-jpl/src/m/contrib/tsantos/AMRexportVTK.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/tsantos/AMRexportVTK.m	(revision 21805)
+++ /issm/trunk-jpl/src/m/contrib/tsantos/AMRexportVTK.m	(revision 21806)
@@ -37,5 +37,5 @@
 end
 
-%this is the result structure
+%this is the result structure (just the Transient solution)
 res_struct=model.results;
 %checking for results
@@ -51,5 +51,5 @@
 		if(size(sol_struct{i},2)>num_of_timesteps);
 			num_of_timesteps=size(sol_struct{i},2);
-      outstep=model.timestepping.time_step*model.settings.output_frequency;
+			outstep=model.timestepping.time_step*model.settings.output_frequency;
 	  end
   end
@@ -101,5 +101,5 @@
 			else
 				timestep = size(sol_struct{j},2);
-	    end
+			end
 			
 			%getting the number of fields in the solution
@@ -120,6 +120,21 @@
 					s='%e\n';
 					fprintf(fid,s,sol_struct{j}(timestep).(fieldnames{k}));
-		    end		
-	    end 
+				end
+			end
+			fprintf(fid,'CELL_DATA %s \n',num2str(num_of_elt));
+			for k=1:num_of_fields
+				if ((numel(sol_struct{j}(timestep).(fieldnames{k})))==num_of_elt);
+					%paraview does not like NaN, replacing
+					nanval=find(isnan(sol_struct{j}(timestep).(fieldnames{k})));
+					sol_struct{j}(timestep).(fieldnames{k})(nanval)=-9999;
+					%also checking for verry small value that mess up
+					smallval=(abs(sol_struct{j}(timestep).(fieldnames{k}))<1.0e-20);
+					sol_struct{j}(timestep).(fieldnames{k})(smallval)=0.0;
+					fprintf(fid,'SCALARS %s float 1 \n',fieldnames{k});
+					fprintf(fid,'LOOKUP_TABLE default\n');
+					s='%e\n';
+					fprintf(fid,s,sol_struct{j}(timestep).(fieldnames{k}));
+				end		
+			end
 	  end
   end