Index: /issm/trunk-jpl/m4/issm_options.m4
===================================================================
--- /issm/trunk-jpl/m4/issm_options.m4	(revision 21483)
+++ /issm/trunk-jpl/m4/issm_options.m4	(revision 21484)
@@ -1896,4 +1896,18 @@
 	AC_MSG_RESULT($HAVE_KRIGING)
 	dnl }}}
+	dnl with-amr{{{
+	AC_ARG_WITH([amr],
+		AS_HELP_STRING([--with-amr = YES],[compile with Adaptive Mesh Refinment capability (default is no)]),
+		[AMR=$withval],[AMR=no]) 
+	AC_MSG_CHECKING(for AMR capability compilation)
+
+	HAVE_AMR=no
+	if test "x$AMR" = "xyes"; then
+		HAVE_AMR=yes
+		AC_DEFINE([_HAVE_AMR_],[1],[with amr capability])
+	fi
+	AM_CONDITIONAL([AMR], [test x$HAVE_AMR = xyes])
+	AC_MSG_RESULT($HAVE_AMR)
+	dnl }}}
 	AX_ANALYSES_SELECTION
 
Index: /issm/trunk-jpl/src/c/Makefile.am
===================================================================
--- /issm/trunk-jpl/src/c/Makefile.am	(revision 21483)
+++ /issm/trunk-jpl/src/c/Makefile.am	(revision 21484)
@@ -530,4 +530,7 @@
 				  ./kml/KMLFileReadUtils.cpp
 #}}}
+#AMR sources  {{{
+amr_sources = ./classes/AdaptiveMeshRefinement.cpp
+#}}}
 #Modules sources{{{
 modules_sources= ./shared/Threads/LaunchThread.cpp\
@@ -613,4 +616,7 @@
 libISSMModules_la_SOURCES += $(kml_sources)
 endif
+if AMR
+libISSMModules_la_SOURCES += $(amr_sources)
+endif
 libISSMModules_la_CXXFLAGS = $(ALLCXXFLAGS)
 if !WINDOWS
Index: /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp	(revision 21484)
+++ /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp	(revision 21484)
@@ -0,0 +1,560 @@
+/*!\file AdaptiveMeshrefinement.cpp
+ * \brief: implementation of the adaptive mesh refinement tool based on NeoPZ library: github.com/labmec/neopz
+ */
+
+ #ifdef HAVE_CONFIG_H
+    #include <config.h>
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include "./AdaptiveMeshRefinement.h"
+
+/*Constructor, copy, clean up and destructor*/
+AdaptiveMeshRefinement::AdaptiveMeshRefinement(){/*{{{*/
+
+    /*Set pointers to NULL*/
+    this->fathermesh    = NULL;
+    this->previousmesh  = NULL;
+    this->hmax          = -1;
+    this->elementswidth = -1;
+
+}
+/*}}}*/
+AdaptiveMeshRefinement::AdaptiveMeshRefinement(const AdaptiveMeshRefinement &cp){/*{{{*/
+    this->operator =(cp);
+}
+/*}}}*/
+AdaptiveMeshRefinement & AdaptiveMeshRefinement::operator =(const AdaptiveMeshRefinement &cp){/*{{{*/
+
+    /*Clean all attributes*/
+    this->CleanUp();
+
+    /*Copy all data*/
+    this->fathermesh    = cp.fathermesh;
+    this->previousmesh  = cp.previousmesh;
+    this->hmax          = cp.hmax;
+    this->elementswidth = cp.elementswidth;
+
+    return *this;
+
+}
+/*}}}*/
+AdaptiveMeshRefinement::~AdaptiveMeshRefinement(){/*{{{*/
+    this->CleanUp();
+}
+/*}}}*/
+void AdaptiveMeshRefinement::CleanUp(){/*{{{*/
+
+    /*Verify and delete all data*/
+    if(this->fathermesh)    delete this->fathermesh;
+    if(this->previousmesh)  delete this->previousmesh;
+
+}
+/*}}}*/
+int AdaptiveMeshRefinement::ClassId() const{/*{{{*/
+    return 13829430; //Antartic area with ice shelves (km^2)
+}
+/*}}}*/
+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->hmax,1);
+        buf.Read(&this->elementswidth,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);
+    }
+    catch(const std::exception& e)
+    {
+        std::cout << "Exception catched! " << e.what() << std::endl;
+        std::cout.flush();
+        DebugStop();
+    }
+}
+/*}}}*/
+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->hmax,1);
+        buf.Write(&this->elementswidth,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
+void AdaptiveMeshRefinement::ExecuteRefinement(int &type_process,double *vx, double *vy, double *masklevelset, long &nvertices, long &nelements, long &nsegments, double** x, double** y, double** z, long*** elements, long*** segments){/*{{{*/
+
+    //ITAPOPO _assert_(this->fathermesh);
+    //ITAPOPO _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("/Users/santos/Desktop/mesh0.vtk");
+    //TPZVTKGeoMesh::PrintGMeshVTK(this->previousmesh,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("/Users/santos/Desktop/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);
+    
+    
+    //std::ofstream file3("/Users/santos/Desktop/mesh2.vtk");
+    //TPZVTKGeoMesh::PrintGMeshVTK(newmesh,file3);
+    
+    this->RefineMeshToAvoidHangingNodes(nohangingnodesmesh);
+    
+    
+    //std::ofstream file4("/Users/santos/Desktop/mesh3.vtk");
+    //TPZVTKGeoMesh::PrintGMeshVTK(nohangingnodesmesh,file4);
+    
+    
+    /*Get new geometric mesh in ISSM data structure*/
+    this->GetMesh(nohangingnodesmesh,nvertices,nelements,nsegments,x,y,z,elements,segments);
+
+    /*Verify the new geometry*/
+    this->CheckMesh(nvertices,nelements,nsegments,this->elementswidth,(*x),(*y),(*z),(*elements),(*segments));
+    
+    delete nohangingnodesmesh;
+
+}
+/*}}}*/
+void AdaptiveMeshRefinement::RefinementProcess(TPZGeoMesh *gmesh,std::vector<TPZVec<REAL> > &GLvec){/*{{{*/
+    
+    /*Refine mesh hmax times*/
+    for(int hlevel=1;hlevel<=this->hmax;hlevel++){
+        
+        /*Set elements to be refined using some criteria*/
+        std::vector<long> ElemVec; //elements without children
+        this->SetElementsToRefine(gmesh,GLvec,hlevel,ElemVec);
+        
+        /*Refine the mesh*/
+        this->RefineMesh(gmesh, ElemVec);
+    }
+    
+}
+/*}}}*/
+void AdaptiveMeshRefinement::RefineMesh(TPZGeoMesh *gmesh, std::vector<long> &ElemVec){/*{{{*/
+
+	/*Refine elements in ElemVec: uniform pattern refinement*/
+	for(long i = 0; i < ElemVec.size(); i++){
+		
+        /*Get geometric element and verify if it has already been refined*/
+        long 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();
+
+}
+/*}}}*/
+void AdaptiveMeshRefinement::RefineMeshToAvoidHangingNodes(TPZGeoMesh *gmesh){/*{{{*/
+    
+    /*Refine elements to avoid hanging nodes: non-uniform refinement*/
+   	const long NElem = gmesh->NElements();
+    for(long 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, long &nvertices, long &nelements, long &nsegments, double** meshX, double** meshY, double** meshZ, long*** elements, long*** segments){/*{{{*/
+
+	/* vertices */
+    long ntotalvertices = gmesh->NNodes();//total
+    
+    /* mesh coords */
+    double *newmeshX = new double[ntotalvertices];
+    double *newmeshY = new double[ntotalvertices];
+    double *newmeshZ = new double[ntotalvertices];
+   	
+   	/* getting mesh coords */
+    for(long 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(long 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]);
+               
+    }
+    
+    long ntotalelements = GeoVec.size();
+    long ** newelements = new long*[ntotalelements];
+
+    if ( !(this->elementswidth == 3) && !(this->elementswidth == 4) && !(this->elementswidth == 6) ) DebugStop();
+
+    for(long i = 0; i < GeoVec.size(); i++){
+
+        newelements[i] = new long[this->elementswidth];
+        for(int j = 0; j < this->elementswidth; j++) newelements[i][j] = GeoVec[i]->NodeIndex(j);
+    }
+    
+    /* segments */
+    std::vector<TPZGeoEl*> SegVec; SegVec.clear();
+    for(long 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]);
+        
+    }
+    
+    long ntotalsegments = SegVec.size();
+    long ** newsegments = new long*[ntotalsegments];
+    
+    for(long i = 0; i < SegVec.size(); i++){
+        
+        newsegments[i] = new long[3];
+        for(int j = 0; j < 2; j++) newsegments[i][j] = SegVec[i]->NodeIndex(j);
+        
+        long neighborindex = SegVec[i]->Neighbour(2).Element()->Index();
+        long neighbourid = -1;
+        
+        for(long 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][2] = neighbourid;
+    }
+    
+    //setting outputs
+    nvertices   = ntotalvertices;
+    nelements   = ntotalelements;
+    nsegments   = ntotalsegments;
+    *meshX      = newmeshX;
+    *meshY      = newmeshY;
+    *meshZ      = newmeshZ;
+    *elements   = newelements;
+    *segments   = newsegments;
+    
+}
+/*}}}*/
+void AdaptiveMeshRefinement::CalcGroundingLinePosition(double *masklevelset,std::vector<TPZVec<REAL> > &GLvec){/*{{{*/
+    
+    /* Find grounding line using elments center point */
+    GLvec.clear();
+    for(long 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!
+        long vertex0 = this->previousmesh->Element(i)->NodeIndex(0);
+        long vertex1 = this->previousmesh->Element(i)->NodeIndex(1);
+        long vertex2 = this->previousmesh->Element(i)->NodeIndex(2);
+        
+        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<long> &ElemVec){/*{{{*/
+
+    if(!gmesh) DebugStop(); //itapopo verificar se usará _assert_
+
+    ElemVec.clear();
+    
+    // itapopo inserir modo de encontrar criterio
+    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<long> &ElemVec){/*{{{*/
+    
+    /* Uniform refinement. This refines the entire mesh */
+    long nelements = gmesh->NElements();
+    for(long i=0;i<nelements;i++){
+        if(gmesh->Element(i)->MaterialId()!=this->GetElemMaterialID()) continue;
+        if(gmesh->Element(i)->HasSubElement()) continue;
+        ElemVec.push_back(i);
+    }
+}
+/*}}}*/
+void AdaptiveMeshRefinement::TagElementsNearGroundingLine(TPZGeoMesh *gmesh,std::vector<TPZVec<REAL> > &GLvec,int &hlevel,std::vector<long> &ElemVec){/*{{{*/
+    
+    /* Tag elements near grounding line */
+    double MaxRegion = 20000.; //itapopo
+    double alpha = 1.0;         //itapopo
+    double MaxDistance = MaxRegion / std::exp(alpha*(hlevel-1));
+    
+    for(long 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 = MaxDistance;
+        
+        for (long 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 < MaxDistance) ElemVec.push_back(i);
+    }
+    
+}
+/*}}}*/
+void AdaptiveMeshRefinement::CreateInitialMesh(long &nvertices, long &nelements, long &nsegments, int &width, double* x, double* y, double* z, long** elements, long** segments){/*{{{*/
+
+	// itapopo _assert_(nvertices>0);
+    // itapoo _assert_(nelements>0);
+    this->SetElementWidth(width);
+
+    /*Verify and creating initial mesh*/
+    //itapopo if(this->fathermesh) _error_("Initial mesh already exists!");
+    
+    this->fathermesh = new TPZGeoMesh();
+	this->fathermesh->NodeVec().Resize( nvertices );
+
+    /*Set the vertices (geometric nodes in NeoPZ context)*/
+	for(long i=0;i<nvertices;i++){
+        
+        /*x,y,z coords*/
+        TPZManVector<REAL,3> coord(3,0.);
+        coord[0]= x[i];
+        coord[1]= y[i];
+        coord[2]= z[i];
+		
+        /*Insert in the mesh*/
+        this->fathermesh->NodeVec()[i].SetCoord(coord);
+		this->fathermesh->NodeVec()[i].SetNodeId(i);
+	}
+	
+	/*Generate the elements*/
+    long index;
+    const int mat = this->GetElemMaterialID();
+    TPZManVector<long> elem(this->elementswidth,0);
+    
+	for(long iel=0;iel<nelements;iel++){
+
+		for(int jel=0;jel<this->elementswidth;jel++) elem[jel]=elements[iel][jel];
+
+        /*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){
+			case 3: this->fathermesh->CreateGeoElement(ETriangle, elem, mat, index, reftype);			break;
+            case 4:	this->fathermesh->CreateGeoElement(ETetraedro, elem, mat, index, reftype); DebugStop(); break;// itapopo _error_("tetra elements must be verified!")	break;
+			case 6:	this->fathermesh->CreateGeoElement(EPrisma, elem, mat, index, reftype);		break;
+            default:	DebugStop();//itapopo _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(long iel=nelements;iel<nelements+nsegments;iel++){
+        
+        boundary[0] = segments[iel-nelements][0];
+        boundary[1] = segments[iel-nelements][1];
+        
+        /*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);
+        
+    }
+    
+    /*Build element and node connectivities*/
+    this->fathermesh->BuildConnectivity();
+    /*Create previous mesh as a copy of father mesh*/
+    this->previousmesh = new TPZGeoMesh(*this->fathermesh);
+    
+}
+/*}}}*/
+void AdaptiveMeshRefinement::SetHMax(int &h){/*{{{*/
+    this->hmax = h;
+}
+/*}}}*/
+void AdaptiveMeshRefinement::SetElementWidth(int &width){/*{{{*/
+    this->elementswidth = width;
+}
+/*}}}*/
+void AdaptiveMeshRefinement::CheckMesh(long &nvertices, long &nelements, long &nsegments,int &width, double* x, double* y, double* z, long** elements, long** segments){/*{{{*/
+
+    /*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( !x || !y || !z || !elements ) DebugStop(); // itapopo verifcar se irá usar o _assert_
+    
+    /*Verify if there are orphan nodes*/
+    std::set<long> elemvertices;
+    elemvertices.clear();
+    
+    for(long i = 0; i < nelements; i++){
+        for(long j = 0; j < width; j++) {
+            elemvertices.insert(elements[i][j]);
+        }
+    }
+    
+    if( elemvertices.size() != nvertices ) DebugStop();//itapopo verificar se irá usar o _assert_
+    
+    //Verify if there are inf or NaN in coords
+    for(long i = 0; i < nvertices; i++) if(isnan(x[i]) || isinf(x[i])) DebugStop();
+    for(long i = 0; i < nvertices; i++) if(isnan(y[i]) || isinf(y[i])) DebugStop();
+    for(long i = 0; i < nvertices; i++) if(isnan(z[i]) || isinf(z[i])) DebugStop();
+    for(long i = 0; i < nvertices; i++){
+        for(long j = 0; j < width; j++){
+            if( isnan(elements[i][j]) || isinf(elements[i][j]) ) DebugStop();
+        }
+    }
+    
+}
+/*}}}*/
Index: /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h	(revision 21484)
+++ /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h	(revision 21484)
@@ -0,0 +1,70 @@
+#ifndef ADAPTIVEMESHREFINEMENT
+#define ADAPTIVEMESHREFINEMENT
+
+#include <iostream>
+#include <fstream>
+#include <string>
+
+#include "pzsave.h"
+#include "pzgmesh.h"
+#include "pzreal.h"
+#include "pzvec.h"
+#include "pzeltype.h"
+
+#include "TPZRefPatternTools.h"
+#include "TPZRefPatternDataBase.h"
+#include "TPZRefPattern.h"
+
+#include "tpzchangeel.h"
+#include "TPZGeoElement.h"
+#include "pzreftriangle.h"
+#include "tpzgeoelrefpattern.h"
+#include "TPZRefPattern.h"
+
+class AdaptiveMeshRefinement : public TPZSaveable {
+
+public:
+
+	/*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
+
+    /*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
+    
+	/*General methods*/
+	void CleanUp();																		//Clean all attributes
+    void SetHMax(int &h);                                                               //Define the max level of refinement
+    void SetElementWidth(int &width);                                                   //Define elements width
+	void ExecuteRefinement(int &type_process,double *vx, double *vy, double *masklevelset, long &nvertices, long &nelements, long &nsegments, double** x, double** y, double** z, long*** elements, long*** segments=NULL);           //A new mesh will be created and refined. This returns the new mesh
+	void CreateInitialMesh(long &nvertices, long &nelements, long &nsegments, int &width, double* x, double* y, double* z, long** elements, long** segments=NULL); //Create a NeoPZ geometric mesh by coords and elements
+    void CheckMesh(long &nvertices, long &nelements, long &nsegments, int &width, double* x, double* y, double* z, long** elements, long** segments=NULL); //Check the consistency of the mesh
+
+private:
+
+	/*Private attributes*/
+    int elementswidth;                                                                 	// geometric nodes for element: 3 == Tria, 4 == Tetra, 6 == Penta
+    int hmax;                                                                          	// max level of refinement
+	TPZGeoMesh *fathermesh;																// Father Mesh is the entire mesh without refinement
+	TPZGeoMesh *previousmesh;															// Previous mesh is a refined mesh of last step
+
+	/*Private methods*/
+    void RefinementProcess(TPZGeoMesh *gmesh,std::vector<TPZVec<REAL> > &GLvec);        // Start the refinement process
+	void RefineMesh(TPZGeoMesh *gmesh, std::vector<long> &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<long> &ElemVec); 	//Define wich elements will be refined
+    void TagAllElements(TPZGeoMesh *gmesh,std::vector<long> &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<long> &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, long &nvertices, long &nelements, long &nsegments, double** meshX, double** meshY, double** meshZ, long*** elements, long*** segments=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
+
+};
+
+#endif
Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 21483)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 21484)
@@ -80,4 +80,8 @@
 	/*Save communicator in the parameters dataset: */
 	this->parameters->AddObject(new GenericParam<ISSM_MPI_Comm>(incomm,FemModelCommEnum));
+
+	#ifdef _HAVE_AMR_
+	this->InitializeAdaptiveRefinement();
+	#endif
 
 	/*Free resources */
@@ -2911,2 +2915,690 @@
 }/*}}}*/
 #endif
+
+#ifdef _HAVE_AMR_
+void FemModel::InitializeAdaptiveRefinement(void){/*{{{*/
+	
+	int my_rank=IssmComm::GetRank();
+	
+	/*Initialize AMR*/
+	//AdaptiveMeshRefinement *AMR = new AdaptiveMeshRefinement();
+	
+	/*Set max level of refinement*/
+	int hmax = 2;
+	//AMR->SetHMax(hmax);
+	
+	/*Get initial mesh*/
+	/*Get total number of elements and total numbers of elements in the initial mesh*/
+	int numberofvertices, numberofelements, numberofsegments;
+	numberofvertices = this->vertices->NumberOfVertices();
+	numberofelements = this->elements->NumberOfElements();
+	numberofsegments = -1; //used on matlab
+
+	/*Get vertices coordinates*/
+	IssmDouble *x = NULL;
+	IssmDouble *y = NULL;
+	IssmDouble *z = NULL;
+	VertexCoordinatesx(&x, &y, &z, this->vertices,false) ;
+
+	/*Get element vertices*/
+	int elementswidth = 3; //just 2D mesh in this version (just tria elements)
+	int* elem_vertices=xNew<int>(elementswidth);
+	Vector<IssmDouble>* vid1= new Vector<IssmDouble>(numberofelements);
+	Vector<IssmDouble>* vid2= new Vector<IssmDouble>(numberofelements);
+	Vector<IssmDouble>* vid3= new Vector<IssmDouble>(numberofelements);
+
+	/*Go through elements, and for each element, object*/
+    for(int i=0;i<this->elements->Size();i++){
+    	Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
+    	element->GetVerticesSidList(elem_vertices);
+    	vid1->SetValue(element->sid,elem_vertices[0],INS_VAL);
+    	vid2->SetValue(element->sid,elem_vertices[1],INS_VAL);
+    	vid3->SetValue(element->sid,elem_vertices[2],INS_VAL);
+    }
+		
+	/*Assemble*/
+    vid1->Assemble();
+    vid2->Assemble();
+    vid3->Assemble();
+
+    /*Serialize*/
+	IssmDouble *id1 = vid1->ToMPISerial();
+   IssmDouble *id2 = vid2->ToMPISerial();
+	IssmDouble *id3 = vid3->ToMPISerial();
+	
+	int **elementsptr;
+	elementsptr = new int*[numberofelements];
+    for(int i=0;i<numberofelements;i++){
+        elementsptr[i] = new int[elementswidth];
+        elementsptr[i][0] = (int)id1[i];
+        elementsptr[i][1] = (int)id2[i];
+        elementsptr[i][2] = (int)id3[i];
+    }
+
+	//if(my_rank==0){
+	//	_printf_("   PRINTING COORDINATES... \n");	    
+   // 		for(long i=0;i<numberofvertices;i++) _printf_("ID: " << i << "\t" << "X: " << x[i] << "\t" << "Y: " << y[i] << "\n");
+   // 	_printf_("   PRINTING ELEMENTS... \n");	   
+   // 	for(int i=0;i<numberofelements;i++) _printf_("El: " << i << "\t" << elementsptr[i][0] << "\t" << elementsptr[i][1] << "\t" << elementsptr[i][2] << "\n");
+	//}
+
+    //AMR->CreateInitialMesh(numberofvertices, numberofelements, 0, elementswidth, x, y, z, elementsptr, NULL);
+
+	/*Free the vectors*/
+	delete x;
+	delete y;
+	delete z;
+	delete vid1;
+	delete vid2;
+	delete vid3;
+	delete id1;
+	delete id2;
+	delete id3;
+	for(int i=0;i<numberofelements;i++) delete elementsptr[i];
+   if(elementsptr) delete elementsptr;
+	xDelete<int>(elem_vertices); 
+}
+/*}}}*/
+FemModel* FemModel::ReMesh(void){/*{{{*/
+	
+	/*All indexing here is in C type: 0..n-1*/
+
+	//////// to test! using the current femmodel
+	int numberofvertices2, numberofelements2;
+	numberofvertices2 = this->vertices->NumberOfVertices();
+	numberofelements2 = this->elements->NumberOfElements();
+
+	/*Get vertices coordinates*/
+	IssmDouble *x = NULL;
+	IssmDouble *y = NULL;
+	IssmDouble *z = NULL;
+	VertexCoordinatesx(&x, &y, &z, this->vertices,false) ;
+
+	/*Get element vertices*/
+	int elementswidth2 = 3; //just 2D mesh in this version (just tria elements)
+	int* elem_vertices2=xNew<int>(elementswidth2);
+	Vector<IssmDouble>* vid1= new Vector<IssmDouble>(numberofelements2);
+	Vector<IssmDouble>* vid2= new Vector<IssmDouble>(numberofelements2);
+	Vector<IssmDouble>* vid3= new Vector<IssmDouble>(numberofelements2);
+
+	/*Go through elements, and for each element, object, report it cpu:*/
+    for(int i=0;i<this->elements->Size();i++){
+    	Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
+    	element->GetVerticesSidList(elem_vertices2);
+    	vid1->SetValue(element->sid,elem_vertices2[0],INS_VAL);
+    	vid2->SetValue(element->sid,elem_vertices2[1],INS_VAL);
+    	vid3->SetValue(element->sid,elem_vertices2[2],INS_VAL);
+    }
+		
+	/*Assemble*/
+    vid1->Assemble();
+    vid2->Assemble();
+    vid3->Assemble();
+
+    /*Serialize*/
+	IssmDouble *id1 = vid1->ToMPISerial();
+    IssmDouble *id2 = vid2->ToMPISerial();
+	IssmDouble *id3 = vid3->ToMPISerial();
+
+	delete vid1;
+	delete vid2;
+	delete vid3;
+	
+	int **elementsptr;
+	elementsptr = new int*[numberofelements2];
+    for(int i=0;i<numberofelements2;i++){
+        elementsptr[i] = new int[elementswidth2];
+        elementsptr[i][0] = (int)id1[i];
+        elementsptr[i][1] = (int)id2[i];
+        elementsptr[i][2] = (int)id3[i];
+    }
+
+   delete id1;
+	delete id2;
+	delete id3;
+	//////////////////////////////////////////////
+
+	int my_rank=IssmComm::GetRank();
+
+	/*Solutions which will be used to refine the elements*/
+	double *vx = NULL;
+	double *vy = NULL;
+	double *masklevelset = NULL;
+
+	int elementswidth = 3;//just 2D mesh, tria elements
+	int numberofelements = this->elements->NumberOfElements();
+	int numberofvertices = this->vertices->NumberOfVertices();
+
+	IssmDouble* elementlevelset = xNew<IssmDouble>(elementswidth);
+	int* elem_vertices = xNew<int>(elementswidth);
+	Vector<IssmDouble>* vmasklevelset = new Vector<IssmDouble>(numberofvertices);
+
+	for(int i=0;i<this->elements->Size();i++){
+		Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
+		element->GetInputListOnVertices(elementlevelset,MaskIceLevelsetEnum);
+		element->GetVerticesSidList(elem_vertices);
+		vmasklevelset->SetValue(elem_vertices[0],elementlevelset[0],INS_VAL);
+    	vmasklevelset->SetValue(elem_vertices[1],elementlevelset[1],INS_VAL);
+    	vmasklevelset->SetValue(elem_vertices[2],elementlevelset[2],INS_VAL);
+	}
+
+	/*Assemble*/
+	vmasklevelset->Assemble();
+
+	/*Serialize*/
+	masklevelset = vmasklevelset->ToMPISerial();
+
+	xDelete<IssmDouble>(elementlevelset);
+	xDelete<int>(elem_vertices);
+	delete vmasklevelset;
+
+	//_printf_("   PRINTING MASKLEVELSET... \n");	    
+	//if(my_rank==0){
+	//	for(int i=0;i<numberofvertices;i++) _printf_("vertex: " << i << "\t" << masklevelset[i] << "\n");
+	//}
+
+	/*Refine the mesh*/
+	double *newx;
+	double *newy;
+	double *newz;
+	int **newelements;
+	int **newsegments;
+	int newnumberofvertices, newnumberofelements, newnumberofsegments;
+	int type_process=1;
+	//AMR->ExecuteRefinement(type_process,vx,vy,masklevelset,newnumberofvertices,newnumberofelements,newnumberofsegments,&newx,&newy,&newz,&newelements,&newsegments);
+	
+	//if(newnumberofvertices<=0 || newnumberofelements<=0 || newnumberofsegments=<0) _error_("Error in the mesh refinement process.");
+
+	delete masklevelset;
+
+	///////////////////////// to test using the current femmodel
+	newx = x;
+	newy = y;
+	newz = z;
+	newelements = elementsptr;
+	newnumberofvertices = numberofvertices2;
+	newnumberofelements = numberofelements2;
+
+	//if(my_rank==0){
+	//	_printf_("   PRINTING COORDINATES... \n");	    
+   // 		for(int i=0;i<newnumberofvertices;i++) _printf_("ID: " << i << "\t" << "X: " << newx[i] << "\t" << "Y: " << newy[i] << "\n");
+   // 	_printf_("   PRINTING ELEMENTS... \n");	   
+   // 	for(int i=0;i<newnumberofelements;i++) _printf_("El: " << i << "\t" << newelements[i][0] << "\t" << newelements[i][1] << "\t" << newelements[i][2] << "\n");
+	//}
+	/////////////////////////
+
+	/*Partitioning the new mesh*/
+	int        *epart          		= NULL; //element partitioning.
+	int        *npart          		= NULL; //node partitioning.
+	int        *newelementslist     = NULL;
+	int        	edgecut				= 1;
+	int 		numflag				= 0;
+	int 		etype 	 			= 1;
+	int 		numprocs 			= IssmComm::GetSize();
+	bool 		*my_elements 		= NULL;
+	int  		*my_vertices 		= NULL;
+	
+	epart 							= xNew<int>(newnumberofelements);
+	npart 							= xNew<int>(newnumberofvertices);
+	newelementslist 				= xNew<int>(newnumberofelements*elementswidth);
+
+	/*Fill the element list to partitioning*/
+	for(int i=0;i<newnumberofelements;i++){ //itapopo os elementos poderão sair do ExecuteRefinement nesse formato
+		for(int j=0;j<elementswidth;j++){
+			newelementslist[elementswidth*i+j] = newelements[i][j]; //C indexing
+		}
+	}
+
+	//if(my_rank==0){
+	//	_printf_("   PRINTING ELEMENTS in ELEMENTLIST... \n");	   
+   // 	for(int i=0;i<newnumberofelements;i++) _printf_("El: " << i << "\t" << newelementslist[i*elementswidth+0] << "\t" << newelementslist[i*elementswidth+1] << "\t" << newelementslist[i*elementswidth+2] << "\n");
+    //}
+
+	/*Partition using Metis:*/
+	if (numprocs>1){
+#ifdef _HAVE_METIS_
+		METIS_PartMeshNodalPatch(&newnumberofelements,&newnumberofvertices, newelementslist, &etype, &numflag, &numprocs, &edgecut, epart, npart);
+#else
+		_error_("metis has not beed installed. Cannot run with more than 1 cpu");
+#endif
+	}
+	else if (numprocs==1){
+		/*METIS does not know how to deal with one cpu only!*/
+		for (int i=0;i<newnumberofelements;i++)	epart[i]=0;
+		for (int i=0;i<newnumberofvertices;i++)	npart[i]=0;
+	}
+	else _error_("At least one processor is required");
+
+	//_printf_("   PRINTING AFTER MESTIS... \n");	    
+
+	my_vertices=xNew<int>(newnumberofvertices);
+	my_elements=xNew<bool>(newnumberofelements);
+	for(int i=0;i<newnumberofvertices;i++) my_vertices[i] = 0;
+	for(int i=0;i<newnumberofelements;i++) my_elements[i] = false;
+
+	/*Start figuring out, out of the partition, which elements belong to this cpu: */
+	for(int i=0;i<newnumberofelements;i++){
+
+		/*!All elements have been partitioned above, only deal with elements for this cpu: */
+		if(my_rank==epart[i]){ 
+			my_elements[i]=true;
+			/*Now that we are here, we can also start building the list of vertices belonging to this cpu partition: we use 
+			 *the  element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing) 
+			 into the vertices coordinates. If we start plugging 1 into my_vertices for each index[n][i] (i=0:2), then my_vertices 
+			 will hold which vertices belong to this partition*/
+			for(int j=0;j<elementswidth;j++){
+				my_vertices[newelementslist[elementswidth*i+j]]=1;
+			}
+		}
+	}
+
+	xDelete<int>(epart);
+	xDelete<int>(npart);
+
+	//_printf_("   PRINTING AFTER MY_VERTICES... \n");	    
+
+	/*Creating new femmodel with new mesh*/
+	FemModel* output = NULL;
+	int       analysis_type;
+
+	output = new FemModel(*this);
+
+	//_printf_("   PRINTING AFTER FEMMODEL... \n");	    
+
+	/*Copy basic attributes:*/
+	output->nummodels = this->nummodels;
+	output->solution_type = this->solution_type;
+	output->analysis_counter = this->analysis_counter;
+
+	//_printf_("   PRINTING AFTER BASIC COPIES... \n");	    
+
+	/*Now, deep copy arrays:*/
+	output->analysis_type_list=xNew<int>(this->nummodels);
+	xMemCpy<int>(output->analysis_type_list,this->analysis_type_list,this->nummodels);
+
+	output->profiler=static_cast<Profiler*>(this->profiler->copy());
+	output->parameters=static_cast<Parameters*>(this->parameters->Copy());
+
+	//_printf_("   PRINTING AFTER PARAMETERS... \n");	    
+
+	/*Creating connectivity table*/
+	int* connectivity = xNew<int>(newnumberofvertices);
+	for(int i=0;i<newnumberofvertices;i++) connectivity[i] = 0;
+
+	for (int i=0;i<newnumberofelements;i++){
+		for (int j=0;j<elementswidth;j++){
+			int vertexid = newelementslist[elementswidth*i+j];
+			_assert_(vertexid>-1 && vertexid<newnumberofvertices);
+			connectivity[vertexid]+=1;
+		}
+	}	
+
+	//_printf_("   PRINTING AFTER CONNECTIVITY... \n");	    
+
+	/*Creating vertices*/
+	output->vertices = new Vertices();
+
+	for(int i=0;i<newnumberofvertices;i++){
+		if(my_vertices[i]){
+			Vertex *newvertex = new Vertex();
+			
+			newvertex->id 		= i+1;
+			newvertex->sid 		= i;
+			newvertex->pid 		= UNDEF;
+
+			newvertex->x         	= newx[i];
+			newvertex->y         	= newy[i];
+			newvertex->z         	= newz[i];
+			newvertex->domaintype	= Domain2DhorizontalEnum;
+			newvertex->sigma		= 0.;
+
+			newvertex->connectivity = connectivity[i];
+
+			output->vertices->AddObject(newvertex);	
+		} 
+	}
+
+	xDelete<int>(connectivity);
+
+	//_printf_("   PRINTING AFTER VERTICES... \n");	    
+
+	/*Creating elements*/
+	/*Just Tria in this version*/
+	output->elements = new Elements();
+
+	for(int i=0;i<newnumberofelements;i++){
+		if(my_elements[i]){
+
+			Tria *newtria = new Tria();
+
+			newtria->id  = i+1;
+			newtria->sid = i;
+			newtria->parameters = NULL;
+
+			newtria->inputs  = new Inputs();
+
+			newtria->nodes    = NULL;
+			newtria->vertices = NULL;
+			newtria->material = NULL;
+			newtria->matpar   = NULL;
+
+			if(this->nummodels>0){
+				newtria->element_type_list=xNew<int>(this->nummodels);
+				for(int j=0;j<nummodels;j++) newtria->element_type_list[j] = 0;
+			}
+			else newtria->element_type_list = NULL;
+
+			/*Element hook*/
+			int matpar_id		= newnumberofelements+1; //retrieve material parameter id (last pointer in femodel->materials)
+			int material_id 	= i+1; // retrieve material_id = i+1;
+
+			/*retrieve vertices ids*/
+			int* vertex_ids = xNew<int>(elementswidth);
+			
+			for(int j=0;j<elementswidth;j++){ 
+				vertex_ids[j]=reCast<int>(newelementslist[elementswidth*i+j]);
+			}
+
+			newtria->numanalyses = this->nummodels;
+			newtria->hnodes      = new Hook*[this->nummodels];
+			newtria->hvertices   = new Hook(&vertex_ids[0],elementswidth);
+			newtria->hmaterial   = new Hook(&material_id,1);
+			newtria->hmatpar     = new Hook(&matpar_id,1);
+			newtria->hneighbors  = NULL;
+
+			/*Initialize hnodes as NULL*/
+			for(int j=0;j<this->nummodels;j++){
+				newtria->hnodes[j]=NULL;
+			}
+
+			/*Clean up*/
+			xDelete<int>(vertex_ids);
+			output->elements->AddObject(newtria);	
+		} 
+	}
+
+	xDelete<int>(newelementslist);
+	//_printf_("   PRINTING AFTER ELEMENTS... \n");	    
+
+	/*Creating materials*/
+	/*Just Matice in this version*/
+	output->materials = new Materials();
+
+	for(int i=0;i<newnumberofelements;i++){
+		if(my_elements[i]){
+			output->materials->AddObject(new Matice(i+1,i,MaticeEnum));	
+		} 
+	}
+	/*This is done to follow CreateElementsVerticesAndMaterials, line 57*/
+	//output->elements->InputDuplicate(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum);
+
+	/*Add new constant material property to materials, at the end: */
+	Matpar *newmatpar = static_cast<Matpar*>(this->materials->GetObjectByOffset(this->materials->Size()-1)->copy());
+	output->materials->AddObject(newmatpar);//put it at the end of the materials
+
+	//_printf_("   PRINTING AFTER MATERIALS... \n");	    
+
+	delete x;
+	delete y;
+	delete z;
+	for(int i=0;i<numberofelements;i++) delete elementsptr[i];
+    if(elementsptr) delete elementsptr;
+
+	//itapopo to test and print the elements and coordinates
+	//output->InitializeAdaptiveRefinement();
+
+	/*Creating nodes*/
+	/*Just SSA (2D) and P1 in this version*/
+	output->nodes = new Nodes();
+
+	int nodecounter = 0;
+	int lid 		= 0;
+	for(int i=0;i<output->nummodels;i++){
+		
+		int analysis_enum = output->analysis_type_list[i];
+
+		//_printf_("   	Analysis type:" << EnumToStringx(analysis_enum) << "\n");
+
+		//itapopo as the domain is 2D, it is not necessary to create nodes for this analysis
+		if(analysis_enum==StressbalanceVerticalAnalysisEnum) continue;	    
+
+		for(int j=0;j<newnumberofvertices;j++){
+			if(my_vertices[j]){
+				
+				Node* newnode = new Node();	
+
+				/*id: */
+				newnode->id            = nodecounter+j+1;
+				newnode->sid           = j;
+				newnode->lid           = lid++;
+				newnode->analysis_enum = analysis_enum;
+
+				/*Initialize coord_system: Identity matrix by default*/
+				for(int k=0;k<3;k++) for(int l=0;l<3;l++) newnode->coord_system[k][l]=0.0;
+				for(int k=0;k<3;k++) newnode->coord_system[k][k]=1.0;
+
+				/*indexing:*/
+				newnode->indexingupdate = true;
+
+				Analysis* analysis = EnumToAnalysis(analysis_enum);
+				int *doftypes = NULL;
+				int numdofs        = analysis->DofsPerNode(&doftypes,Domain2DhorizontalEnum,SSAApproximationEnum);
+				newnode->indexing.Init(numdofs,doftypes);
+				xDelete<int>(doftypes);
+				delete analysis;
+
+				if(analysis_enum==StressbalanceAnalysisEnum)
+				 newnode->SetApproximation(SSAApproximationEnum);
+				else
+				 newnode->SetApproximation(0);
+
+				/*Stressbalance Horiz*/
+				if(analysis_enum==StressbalanceAnalysisEnum){
+					// itapopo
+					/*Coordinate system provided, convert to coord_system matrix*/
+					//XZvectorsToCoordinateSystem(&this->coord_system[0][0],&iomodel->Data(StressbalanceReferentialEnum)[j*6]);
+					//_assert_(sqrt( coord_system[0][0]*coord_system[0][0] + coord_system[1][0]*coord_system[1][0]) >1.e-4);
+
+				}
+
+				output->nodes->AddObject(newnode);
+			}
+		}
+
+		if(output->nodes->Size()) nodecounter = output->nodes->MaximumId();
+	}
+	
+	//_printf_("   Old node size: " << this->nodes->Size() << " \n");	    
+
+	//for(int i=0;i<this->nummodels;i++){
+	//	int analysis_type = this->analysis_type_list[i];
+	//	_printf_("   	Analysis type:" << EnumToStringx(analysis_type) << " size: "<< this->nodes->NumberOfNodes(analysis_type) << " NDofs: " <<  this->nodes->NumberOfDofs(analysis_type, GsetEnum) << " \n");
+	//}
+	
+	//_printf_("   New number of nodes: " << output->nodes->Size() << " \n");	    
+	//for(int i=0;i<output->nummodels;i++){
+	//	int analysis_type = output->analysis_type_list[i];
+	//	_printf_("   	Analysis type:" << EnumToStringx(analysis_type) << " size: "<< output->nodes->NumberOfNodes(analysis_type) << " NDofs: " <<  this->nodes->NumberOfDofs(analysis_type, GsetEnum) << " \n");
+	//}
+	
+	//_printf_("		Old nodes deep echo: \n");
+	//this->nodes->DeepEcho();
+
+	//_printf_("		New nodes deep echo: \n");
+	//output->nodes->DeepEcho();
+
+	//_printf_("   PRINTING AFTER NODES... \n");	    
+
+	/*Create constraints*/
+	IssmDouble *spcvx = NULL;
+	IssmDouble *spcvy = NULL;
+
+	int numberofnodes_analysistype 	= this->nodes->NumberOfNodes(StressbalanceAnalysisEnum);
+	
+	//_printf_("   Number of nodes: " << numberofnodes_analysistype << " \n");	    
+
+	Vector<IssmDouble>* vspcvx 		= new Vector<IssmDouble>(numberofnodes_analysistype);
+	Vector<IssmDouble>* vspcvy		= new Vector<IssmDouble>(numberofnodes_analysistype);
+
+	IssmDouble BigNumber 			= 1.e8;
+
+	for(int i=0;i<numberofnodes_analysistype;i++){
+		vspcvx->SetValue(i,BigNumber,INS_VAL);
+		vspcvy->SetValue(i,BigNumber,INS_VAL);
+	}
+
+	IssmDouble absmaxspcvx = 0;
+	IssmDouble absmaxspcvy = 0;
+
+	for(int i=0;i<this->constraints->Size();i++){
+		SpcStatic* spc 		= xDynamicCast<SpcStatic*>(this->constraints->GetObjectByOffset(i));
+		int dof 			= spc->GetDof();
+		int node 			= spc->GetNodeId();
+		IssmDouble spcvalue	= spc->GetValue(); 
+		int nodeindex		= node-1;
+
+		if(dof==0) {
+			vspcvx->SetValue(nodeindex,spcvalue,INS_VAL);
+			if(fabs(spcvalue)>absmaxspcvx) absmaxspcvx = fabs(spcvalue);
+		}
+		else {
+			vspcvy->SetValue(nodeindex,spcvalue,INS_VAL);
+			if(fabs(spcvalue)>absmaxspcvy) absmaxspcvy = fabs(spcvalue);
+		}
+	}
+
+	/*Assemble*/
+	vspcvx->Assemble();
+	vspcvy->Assemble();
+
+	/*Serialize*/
+	spcvx = vspcvx->ToMPISerial();
+	spcvy = vspcvy->ToMPISerial();
+
+	/*Free the data*/
+	delete vspcvx;
+	delete vspcvy;
+	
+	double *newspcvx 	= NULL;
+	double *newspcvy 	= NULL;
+	int *oldelements 	= newelementslist; //itapopo
+	double *oldx		= x; //itapopo
+	double *oldy		= y; //itapopo
+	int nods_data		= numberofnodes_analysistype;
+	int nels_data 		= newnumberofelements;
+	int M_data 			= numberofnodes_analysistype;
+	int N_data  		= 1;
+	int N_interp 		= newnumberofvertices;//itapopo
+	Options *options   	= NULL;
+
+	//itapopo voltar aqui
+	InterpFromMeshToMesh2dx(&newspcvx,oldelements,oldx,oldy,nods_data,nels_data,spcvx,M_data,N_data,newx,newy,N_interp,options);
+	InterpFromMeshToMesh2dx(&newspcvx,oldelements,oldx,oldy,nods_data,nels_data,spcvx,M_data,N_data,newx,newy,N_interp,options);
+
+	output->constraints = new Constraints();
+
+	nodecounter 			= 0; //itapopo deve começar pelo primeiro nó do StressbalanceAnalysis
+	int count 				= 0;
+	int constraintcounter 	= 0; //itapopo
+	IssmDouble eps			= 1.e-2;
+
+	for(int i=0;i<newnumberofvertices;i++){
+		if(my_vertices[i])
+		/*spcvx*/
+		if(fabs(spcvx[i]) < absmaxspcvx+eps){//itapopo
+			output->constraints->AddObject(new SpcStatic(constraintcounter+count+1,nodecounter+i+1,0,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
+			count++;
+		}
+	}
+	count=0;
+	for(int i=0;i<newnumberofvertices;i++){
+		if(my_vertices[i])
+		/*spcvy*/
+		if(fabs(spcvy[i]) < absmaxspcvy+eps){//itapopo
+			output->constraints->AddObject(new SpcStatic(constraintcounter+count+1,nodecounter+i+1,1,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
+			count++;
+		}
+
+	}
+
+	//_printf_("   Old constraints size: " << this->constraints->Size() << " \n");	    
+	//this->constraints->DeepEcho();
+
+	//_printf_("\n");
+
+	//_printf_("   New constraints size: " << output->constraints->Size() << " \n");
+	
+	//_printf_("SPCVX\n");
+	//for(int i = 0;i<newnumberofvertices;i++){
+	//	_printf_("value: " << spcvx[i] << "\n");
+	//}	    
+	//_printf_("SPCVY\n");
+	//for(int i = 0;i<newnumberofvertices;i++){
+	//	_printf_("value: " << spcvy[i] << "\n");
+	//}
+
+	//output->constraints->DeepEcho();
+
+	//_printf_("   PRINTING AFTER CONSTRAINTS... \n");	    
+
+	//_printf_("   Old loads size: " << this->loads->Size() << " \n");	    
+	//this->loads->DeepEcho();
+
+	//_printf_("   PRINTING AFTER LOADS... \n");	    
+
+	// output->loads=static_cast<Loads*>(this->loads->Copy());
+	// output->constraints=static_cast<Constraints*>(this->constraints->Copy());
+	// output->results=static_cast<Results*>(this->results->Copy());
+
+	// /*reset hooks for elements, loads and nodes: */
+	// output->elements->ResetHooks();
+	// output->loads->ResetHooks();
+	// output->materials->ResetHooks();
+
+	// /*do the post-processing of the datasets to get an FemModel that can actually run analyses: */
+	// for(i=0;i<nummodels;i++){
+	// 	analysis_type=output->analysis_type_list[i];
+	// 	output->SetCurrentConfiguration(analysis_type);
+	// 	if(i==0) VerticesDofx(output->vertices,output->parameters); //only call once, we only have one set of vertices
+	// 	SpcNodesx(output->nodes,output->constraints,output->parameters,analysis_type);
+	// 	NodesDofx(output->nodes,output->parameters,analysis_type);
+	// 	ConfigureObjectsx(output->elements,output->loads,output->nodes,output->vertices,output->materials,output->parameters);
+	// }
+
+	// /*Reset current configuration: */
+	// analysis_type=output->analysis_type_list[analysis_counter];
+	// output->SetCurrentConfiguration(analysis_type);
+
+	/** TODO
+
+	- generate the required input objects for NeoPZ
+		AMR has fathermesh and previousmesh. On first call, these meshes are generated. 
+
+	- call NeoPZ
+		This creates a newmesh which will be refined. 
+
+	- get the new mesh (index,x,y) from NeoPZ
+		Is is doing by GetNewMesh method.
+
+	- Create a new FemModel* that will be consistent with the new mesh
+		It can be done by a method. (attention with CPU #)
+
+	- Initialize new FemModel based on new mesh (and copy the old FemModel parameters)
+		It can be done by a method. (attention with CPU #)
+
+	- The last and most difficult thing to do is to update the inputs of the elements of the new mesh:
+		It can be done in just one method, which calls GatherInputs, InterpFromMeshToMesh and InputUpdateFromVector
+
+		- CPU #0 will gather all inputs from FemModel
+        	- call InterpFromMeshToMesh to interpolate them onto the new Mesh
+        	- broadcast the new fields to all cpus
+		- call InputUpdateFromVector so that all inputs are updated
+	
+	- return new FemModel
+	
+	*/
+
+	return output;
+}
+/*}}}*/
+#endif
Index: /issm/trunk-jpl/src/c/classes/FemModel.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.h	(revision 21483)
+++ /issm/trunk-jpl/src/c/classes/FemModel.h	(revision 21484)
@@ -137,4 +137,10 @@
 		void InitFromBuffers(char* buffer, int buffersize, char* toolkits, int solution_type,bool trace,IssmPDouble* X=NULL);
 		#endif
+
+		#ifdef _HAVE_AMR_
+		/*Adaptive mesh refinement methods*/
+		void InitializeAdaptiveRefinement(void);
+		FemModel* ReMesh(void);
+		#endif
 };
 		
Index: /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp	(revision 21483)
+++ /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp	(revision 21484)
@@ -37,35 +37,16 @@
 Matice::Matice(int matice_mid,int index, IoModel* iomodel){/*{{{*/
 
-	/*Intermediaries:*/
-	int    matice_eid;
-
-	/*Initialize id*/
-	this->mid=matice_mid;
-
-	/*Hooks: */
-	matice_eid=index+1;
-	this->helement=new Hook(&matice_eid,1);
-	this->element=NULL;
-
-	 /*Other perporties*/
-   int    materialtype;
+	 /*Get material type and initialize object*/
+   int materialtype;
    iomodel->FindConstant(&materialtype,"md.materials.type");
-   if(materialtype==MatdamageiceEnum){
-		this->isdamaged = true;
-		this->isenhanced = false;
-	}
-	else if(materialtype==MaticeEnum){
-		this->isdamaged = false;
-		this->isenhanced = false;
-	}
-	else if(materialtype==MatenhancediceEnum){
-		this->isdamaged = false;
-		this->isenhanced = true;
-	}
-   else _error_("Material type not recognized");
+	this->Init(matice_mid,index,materialtype);
+
+}
+/*}}}*/
+Matice::Matice(int matice_mid,int index,int materialtype){/*{{{*/
+
+	this->Init(matice_mid,index,materialtype);
 	return;
-
-}
-/*}}}*/
+} /*}}}*/
 Matice::~Matice(){/*{{{*/
 	delete helement;
@@ -73,4 +54,34 @@
 }
 /*}}}*/
+void Matice::Init(int matice_mid,int index,int materialtype){/*{{{*/
+
+	/*Initialize id*/
+	this->mid=matice_mid;
+
+	/*Hooks: */
+	int matice_eid=index+1;
+	this->helement=new Hook(&matice_eid,1);
+	this->element=NULL;
+
+	/*Material specific properties*/
+	switch(materialtype){
+		case MatdamageiceEnum:
+			this->isdamaged = true;
+			this->isenhanced = false;
+			break;
+		case MaticeEnum:
+			this->isdamaged = false;
+			this->isenhanced = false;
+			break;
+		case MatenhancediceEnum:
+			this->isdamaged = false;
+			this->isenhanced = true;
+			break;
+		default:
+			_error_("Material type not recognized");
+	}
+
+	return;
+} /*}}}*/
 
 /*Object virtual functions definitions:*/
Index: /issm/trunk-jpl/src/c/classes/Materials/Matice.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Materials/Matice.h	(revision 21483)
+++ /issm/trunk-jpl/src/c/classes/Materials/Matice.h	(revision 21484)
@@ -35,5 +35,7 @@
 		Matice();
 		Matice(int mid,int i, IoModel* iomodel);
+		Matice(int mid,int i, int materialtype);
 		~Matice();
+		void Init(int mid,int i, int materialtype);
 		/*}}}*/
 		/*Object virtual functions definitions:{{{ */
Index: /issm/trunk-jpl/src/c/classes/Node.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Node.cpp	(revision 21483)
+++ /issm/trunk-jpl/src/c/classes/Node.cpp	(revision 21484)
@@ -497,4 +497,9 @@
 }
 /*}}}*/
+void  Node::SetApproximation(int in_approximation){/*{{{*/
+	
+	this->approximation = in_approximation;
+}
+/*}}}*/
 int  Node::GetNumberOfDofs(int approximation_enum,int setenum){/*{{{*/
 
Index: /issm/trunk-jpl/src/c/classes/Node.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Node.h	(revision 21483)
+++ /issm/trunk-jpl/src/c/classes/Node.h	(revision 21484)
@@ -85,4 +85,5 @@
 		void  VecMerge(Vector<IssmDouble>* ug, IssmDouble* vector_serial,int setenum);
 		void  VecReduce(Vector<IssmDouble>* vector, IssmDouble* ug_serial,int setnum);
+		void  SetApproximation(int in_approximation);
 };