Index: /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp	(revision 21527)
+++ /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp	(revision 21528)
@@ -3,5 +3,5 @@
  */
 
- #ifdef HAVE_CONFIG_H
+#ifdef HAVE_CONFIG_H
     #include <config.h>
 #else
@@ -52,4 +52,5 @@
 	 this->hmax=-1;
 	 this->elementswidth=-1;
+	 gRefDBase.clear();
 }
 /*}}}*/
@@ -128,5 +129,9 @@
 /*Mesh refinement methods*/
 #include "TPZVTKGeoMesh.h" //itapopo
-void AdaptiveMeshRefinement::ExecuteRefinement(int &type_process,double *vx, double *vy, double *masklevelset, int &nvertices, int &nelements, int &nsegments, double** x, double** y, double** z, int*** elements, int*** segments){/*{{{*/
+#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*/
+	/*NEOPZ works only in C indexing*/
 
     //itapopo _assert_(this->fathermesh);
@@ -171,8 +176,8 @@
     
     /*Get new geometric mesh in ISSM data structure*/
-    this->GetMesh(nohangingnodesmesh,nvertices,nelements,nsegments,x,y,z,elements,segments);
+    this->GetMesh(nohangingnodesmesh,nvertices,nelements,nsegments,px,py,pz,pelements,psegments);
 	 
     /*Verify the new geometry*/
-    this->CheckMesh(nvertices,nelements,nsegments,this->elementswidth,x,y,z,elements,segments);
+    this->CheckMesh(nvertices,nelements,nsegments,this->elementswidth,px,py,pz,pelements,psegments);
      
     delete nohangingnodesmesh;
@@ -255,5 +260,8 @@
 }
 /*}}}*/
-void AdaptiveMeshRefinement::GetMesh(TPZGeoMesh *gmesh, int &nvertices, int &nelements, int &nsegments, double** meshX, double** meshY, double** meshZ, int*** elements, int*** segments){/*{{{*/
+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 */
@@ -261,7 +269,7 @@
     
     /* mesh coords */
-    double *newmeshX = new double[ntotalvertices];
-    double *newmeshY = new double[ntotalvertices];
-    double *newmeshZ = new double[ntotalvertices];
+	 double* newmeshX = xNew<IssmDouble>(ntotalvertices);
+    double* newmeshY = xNew<IssmDouble>(ntotalvertices);
+    double* newmeshZ = xNew<IssmDouble>(ntotalvertices);
    	
    /* getting mesh coords */
@@ -283,12 +291,11 @@
     
     int ntotalelements = (int)GeoVec.size();
-    int ** newelements = new int*[ntotalelements];
+    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++){
-        newelements[i] = new int[this->elementswidth];
-        for(int j = 0; j < this->elementswidth; j++) newelements[i][j] = (int)GeoVec[i]->NodeIndex(j);
-    }
+    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 */
@@ -301,10 +308,10 @@
     
     int ntotalsegments = (int)SegVec.size();
-    int ** newsegments = new int*[ntotalsegments];
-    
-    for(int i = 0; i < SegVec.size(); i++){
-        
-        newsegments[i] = new int[3];
-        for(int j = 0; j < 2; j++) newsegments[i][j] = SegVec[i]->NodeIndex(j);
+    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();
@@ -319,16 +326,16 @@
         
         if(neighbourid==-1) DebugStop(); //itapopo talvez passar para _assert_
-        newsegments[i][2] = neighbourid;
+        newsegments[i*3+2] = neighbourid+1;//C to Matlab indexing
     }
     
     //setting outputs
-    nvertices   = ntotalvertices;
-    nelements   = ntotalelements;
-    nsegments   = ntotalsegments;
-    *meshX      = newmeshX;
-    *meshY      = newmeshY;
-    *meshZ      = newmeshZ;
-    *elements   = newelements;
-    *segments   = newsegments;
+    nvertices  = ntotalvertices;
+    nelements  = ntotalelements;
+    nsegments  = ntotalsegments;
+    *px		   = newmeshX;
+    *py		   = newmeshY;
+    *pz		   = newmeshZ;
+    *pelements = newelements;
+    *psegments = newsegments;
     
 }
@@ -381,9 +388,9 @@
     ElemVec.clear();
     
-    // itapopo inserir modo de encontrar criterio
-    if(false) this->TagAllElements(gmesh,ElemVec); //uniform, refine all elements!
+    // itapopo inserir modo de encontrar criterio TESTING!!!! Come to false!
+	 if(true) this->TagAllElements(gmesh,ElemVec); //uniform, refine all elements!
 
     /* Adaptive refinement. This refines some elements following some criteria*/
-    this->TagElementsNearGroundingLine(gmesh, GLvec, hlevel, ElemVec);
+    //this->TagElementsNearGroundingLine(gmesh, GLvec, hlevel, ElemVec);
 
 }
@@ -437,9 +444,12 @@
 }
 /*}}}*/
-void AdaptiveMeshRefinement::CreateInitialMesh(int &nvertices, int &nelements, int &nsegments, int &width, double* x, double* y, double* z, int** elements, int** segments){/*{{{*/
-
+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*/
+	
 	// itapopo _assert_(nvertices>0);
     // itapoo _assert_(nelements>0);
-   this->SetElementWidth(width);
+	this->SetElementWidth(width);
 
     /*Verify and creating initial mesh*/
@@ -450,61 +460,59 @@
 
     /*Set the vertices (geometric nodes in NeoPZ context)*/
-	for(int 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);
+	for(int 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);
+	long index;
+   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][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");
+		for(int jel=0;jel<this->elementswidth;jel++) elem[jel]=elements[iel*this->elementswidth+jel]-1;//Convert Matlab to C indexing
+
+      /*reftype = 0: uniform, fast / reftype = 1: uniform and non-uniform (avoid hanging nodes), it is not too fast */
+      const int reftype = 1;
+      switch(this->elementswidth){
+			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");
 		}
         
-        /*Define the element ID*/        
-        this->fathermesh->ElementVec()[index]->SetId(iel);
-        
+      /*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][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);
-    
+   /*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);
+	}
+    
+   /*Build element and node connectivities*/
+   this->fathermesh->BuildConnectivity();
+    
+	/*Create previous mesh as a copy of father mesh*/
+   this->previousmesh = new TPZGeoMesh(*this->fathermesh);
+	 
+	/*Set refinement patterns*/
+	this->SetRefPatterns();
+
 }
 /*}}}*/
@@ -517,5 +525,39 @@
 }
 /*}}}*/
-void AdaptiveMeshRefinement::CheckMesh(int &nvertices, int &nelements, int &nsegments,int &width, double** x, double** y, double** z, int*** elements, int*** segments){/*{{{*/
+void AdaptiveMeshRefinement::SetRefPatterns(){/*{{{*/
+
+	/*Initialize the global variable of refinement patterns*/ 
+	gRefDBase.InitializeUniformRefPattern(EOned);
+	gRefDBase.InitializeUniformRefPattern(ETriangle);
+    
+    //gRefDBase.InitializeRefPatterns();
+	/*Insert specifics patterns to ISSM core*/
+	std::string filepath	 = REFPATTERNDIR;
+   std::string filename1 = filepath + "/2D_Triang_Rib_3.rpt";
+   std::string filename2 = filepath + "/2D_Triang_Rib_4.rpt";
+   std::string filename3 = filepath + "/2D_Triang_Rib2_Side_3_4.rpt";
+   std::string filename4 = filepath + "/2D_Triang_Rib2_Side_3_4permuted.rpt";
+   std::string filename5 = filepath + "/2D_Triang_Rib2_Side_3_5.rpt";
+	std::string filename6 = filepath + "/2D_Triang_Rib2_Side_3_5permuted.rpt";
+   std::string filename7 = filepath + "/2D_Triang_Rib_5.rpt";
+
+   TPZAutoPointer<TPZRefPattern> refpat1 = new TPZRefPattern(filename1);
+   TPZAutoPointer<TPZRefPattern> refpat2 = new TPZRefPattern(filename2);
+   TPZAutoPointer<TPZRefPattern> refpat3 = new TPZRefPattern(filename3);
+	TPZAutoPointer<TPZRefPattern> refpat4 = new TPZRefPattern(filename4);
+   TPZAutoPointer<TPZRefPattern> refpat5 = new TPZRefPattern(filename5);
+   TPZAutoPointer<TPZRefPattern> refpat6 = new TPZRefPattern(filename6);
+   TPZAutoPointer<TPZRefPattern> refpat7 = new TPZRefPattern(filename7);
+   
+   if(!gRefDBase.FindRefPattern(refpat1)) gRefDBase.InsertRefPattern(refpat1);
+   if(!gRefDBase.FindRefPattern(refpat2)) gRefDBase.InsertRefPattern(refpat2);
+   if(!gRefDBase.FindRefPattern(refpat3)) gRefDBase.InsertRefPattern(refpat3);
+   if(!gRefDBase.FindRefPattern(refpat4)) gRefDBase.InsertRefPattern(refpat4);
+   if(!gRefDBase.FindRefPattern(refpat5)) gRefDBase.InsertRefPattern(refpat5);
+   if(!gRefDBase.FindRefPattern(refpat6)) gRefDBase.InsertRefPattern(refpat6);
+   if(!gRefDBase.FindRefPattern(refpat7)) gRefDBase.InsertRefPattern(refpat7);
+}
+/*}}}*/
+void AdaptiveMeshRefinement::CheckMesh(int &nvertices,int &nelements,int &nsegments,int &width,double** px,double** py,double** pz,int** pelements, int** psegments){/*{{{*/
 
     /*Basic verification*/
@@ -524,5 +566,5 @@
     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_
+    if( !px || !py || !pz || !pelements ) DebugStop(); // itapopo verifcar se irá usar o _assert_
     
     /*Verify if there are orphan nodes*/
@@ -531,5 +573,5 @@
     for(int i = 0; i < nelements; i++){
         for(int j = 0; j < width; j++) {
-            elemvertices.insert((*elements)[i][j]);
+            elemvertices.insert((*pelements)[i*width+j]);
 		  }
 	 }
@@ -538,14 +580,14 @@
 	
     //Verify if there are inf or NaN in coords
-    for(int i = 0; i < nvertices; i++) if(isnan((*x)[i]) || isinf((*x)[i])) DebugStop();
-    for(int i = 0; i < nvertices; i++) if(isnan((*y)[i]) || isinf((*y)[i])) DebugStop();
-    for(int i = 0; i < nvertices; i++) if(isnan((*z)[i]) || isinf((*z)[i])) DebugStop();
+    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((*elements)[i][j]) || isinf((*elements)[i][j]) ) DebugStop();
-        }
-    }
-    
-}
-/*}}}*/
+            if( isnan((*pelements)[i*width+j]) || isinf((*pelements)[i*width+j]) ) DebugStop();
+        }
+    }
+    
+}
+/*}}}*/
Index: /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h	(revision 21527)
+++ /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h	(revision 21528)
@@ -10,5 +10,9 @@
 
 /*NeoPZ includes*/
-/*REAL and STATE definitions, NeoPZ variables*/
+/*REAL and STATE definitions, NeoPZ variables itapopo should be read by NeoPZ's config.h*/
+#ifndef REFPATTERNDIR
+	# define REFPATTERNDIR "/home/santos/trunk-jpl/externalpackages/neopz/install/include/refpatterns"
+#endif
+
 #ifndef REALdouble
 	#define REALdouble
@@ -55,7 +59,7 @@
    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, int &nvertices, int &nelements, int &nsegments, double** x, double** y, double** z, int*** elements, int*** segments=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 CheckMesh(int &nvertices, int &nelements, int &nsegments, int &width, double** x, double** y, double** z, int*** elements, int*** segments=NULL); // Check the consistency of the mesh
+	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 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
 
 private:
@@ -75,8 +79,8 @@
    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** meshX, double** meshY, double** meshZ, int*** elements, int*** segments=NULL); // Return coords and elements in ISSM data structure
+	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 SetRefPatterns();																	// Initialize and insert the refinement patterns
 };
 
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 21527)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 21528)
@@ -1578,5 +1578,9 @@
 	         name==HydrologyHeadOldEnum ||		
 				name==StressbalanceConvergenceNumStepsEnum || 
-				name==MeshVertexonbaseEnum 
+				name==MeshVertexonbaseEnum || 
+				name==FrictionPEnum ||
+				name==FrictionQEnum ||
+				name==LoadingforceXEnum ||
+				name==LoadingforceYEnum
 
 				) {
Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 21527)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 21528)
@@ -82,8 +82,6 @@
 	#ifdef _HAVE_NEOPZ_
 	this->InitializeAdaptiveRefinement();
-	FemModel *Test=this->ReMesh();//itapopo: just to test!;
-	//Test->elements->DeepEcho();
-	//Test->CleanUp();
-	//delete Test;
+	FemModel *Test=this->ReMesh();//itapopo: just to test!;;
+	//if(IssmComm::GetRank()==0) Test->elements->DeepEcho();
 	#endif
 
@@ -2924,69 +2922,27 @@
 void FemModel::InitializeAdaptiveRefinement(void){/*{{{*/
 	
-	int my_rank=IssmComm::GetRank();
-	
-	/*Initialize AMR*/
-	if(my_rank==0){
-		/*Just CPU #0 should keep AMR object*/
+	/*Define variables*/
+	int my_rank				= IssmComm::GetRank();
+	this->amr				= NULL;//initialize amr as NULL
+	int numberofvertices = this->vertices->NumberOfVertices();
+	int numberofelements = this->elements->NumberOfElements();
+	int numberofsegments = 0; //used on matlab
+	IssmDouble* x			= NULL;
+	IssmDouble* y			= NULL;
+	IssmDouble* z			= NULL;
+	int* elements			= NULL;
+	int elementswidth		= this->GetElementsWidth(); //just tria elements in this version. Itapopo:
+
+	/*Get vertices coordinates of the coarse mesh (father mesh)*/
+	/*elements comes in Matlab indexing*/
+	this->GetMesh(this,&x,&y,&z,&elements);
+	
+	/*Create initial mesh (coarse mesh) in neopz data structure*/ 
+	/*Just CPU #0 should keep AMR object*/
+   if(my_rank==0){ 
 		this->amr = new AdaptiveMeshRefinement();
-		int hmax = 2; //itapopo: it must be defined by interface. Using 2 just to test
+		int hmax	 = 0; //itapopo: it must be defined by interface. Using 2 just to test
 		this->amr->SetHMax(hmax); //Set max level of refinement
-	} 
-	else{
-		this->amr=NULL;
-	}
-	
-	/*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, get vertices*/
-   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];
-	}
-
-	/*Create initial mesh (coarse mesh) in neopz data structure*/ 
-   if(my_rank==0){ 
-		int nsegments=0;
-		this->amr->CreateInitialMesh(numberofvertices, numberofelements, nsegments, elementswidth, x, y, z, elementsptr, NULL);
+		this->amr->CreateInitialMesh(numberofvertices,numberofelements,numberofsegments,elementswidth,x,y,z,elements,NULL);
 	}
 
@@ -2995,13 +2951,5 @@
 	xDelete<IssmDouble>(y);
 	xDelete<IssmDouble>(z);
-	delete vid1;
-	delete vid2;
-	delete vid3;
-	xDelete<IssmDouble>(id1);
-	xDelete<IssmDouble>(id2);
-	xDelete<IssmDouble>(id3);			
-	for(int i=0;i<numberofelements;i++) delete [] elementsptr[i];
-   if(elementsptr) delete [] elementsptr;
-	xDelete<int>(elem_vertices);
+	xDelete<int>(elements);
 
 }
@@ -3009,26 +2957,24 @@
 FemModel* FemModel::ReMesh(void){/*{{{*/
 	
-	int numprocs=IssmComm::GetSize();
-	if(numprocs>1) _error_("Multithreading not tested yet. Run with 1 cpu.");
-	
 	/*Variables*/
-	IssmDouble *newx=NULL;
-	IssmDouble *newy=NULL;
-	IssmDouble *newz=NULL;
-	int *newelementslist=NULL;
-	int newnumberofvertices=-1;
-	int newnumberofelements=-1;
-	bool* my_elements=NULL; 
-	int* my_vertices=NULL;
-	int elementswidth=3;//just tria elements in this version
+	IssmDouble *newx			= NULL;
+	IssmDouble *newy			= NULL;
+	IssmDouble *newz			= NULL;
+	int *newelementslist		= NULL;
+	int newnumberofvertices	= -1;
+	int newnumberofelements = -1;
+	bool* my_elements			= NULL; 
+	int* my_vertices			= NULL;
+	int elementswidth			= this->GetElementsWidth();//just tria elements in this version
 
 	/*Execute refinement and get the new mesh*/
+	/*newelementslist come in Matlab indexing*/
 	this->ExecuteRefinement(newnumberofvertices,newnumberofelements,&newx,&newy,&newz,&newelementslist);
 
 	/*Partitioning the new mesh. Maybe ElementsAndVerticesPartitioning.cpp could be modified to set this without iomodel.*/
 	this->ElementsAndVerticesPartitioning(newnumberofvertices,newnumberofelements,elementswidth,newelementslist,&my_elements,&my_vertices);
-	
+
 	/*Creating new femmodel with new mesh*/
-	FemModel* output=new FemModel(*this);
+	FemModel* output = new FemModel(*this);
 
 	/*Copy basic attributes:*/
@@ -3055,11 +3001,9 @@
 	output->vertices=new Vertices();
 	this->CreateVertices(newnumberofvertices,newnumberofelements,elementswidth,newelementslist,my_vertices,newx,newy,newz,output->vertices);
-
+ 
 	/*Creating elements*/
 	/*Just Tria in this version*/
 	output->elements=new Elements();
-	this->CreateElements(newnumberofelements,elementswidth,newelementslist,my_elements,output->elements);
-	/*Cleanup*/
-//	xDelete<int>(newelementslist); itapopo	    
+	this->CreateElements(newnumberofelements,elementswidth,newelementslist,my_elements,output->elements); 
 
 	/*Creating materials*/
@@ -3089,5 +3033,5 @@
 	/*Create constraints*/
 	output->constraints=new Constraints();
-	this->CreateConstraints(newnumberofvertices,newnumberofelements,newelementslist,newx,newy,my_vertices,output->constraints);
+	this->CreateConstraints(newnumberofvertices,newnumberofelements,newx,newy,my_vertices,output->constraints);
 	
 	output->elements->Presort();
@@ -3117,14 +3061,22 @@
 	}
 
-	/*Finally: interpolate all inputs*/
-	printf("-------------- file: FemModel.cpp line: %i\n",__LINE__); 
+	/*Finally: interpolate all inputs and insert them into the new elements.*/
 	this->InterpolateInputs(output);
-	
+	GetMaskOfIceVerticesLSMx(output);
+
 	/*Reset current configuration: */
 	analysis_type=output->analysis_type_list[this->analysis_counter];
 	output->SetCurrentConfiguration(analysis_type);
 
+	/*Cleanup*/
+	xDelete<IssmDouble>(newx);
+	xDelete<IssmDouble>(newy);
+	xDelete<IssmDouble>(newz);
+	xDelete<int>(newelementslist);
+	xDelete<int>(my_vertices);
+	xDelete<bool>(my_elements);
 
 	return output;
+
 }
 /*}}}*/
@@ -3190,7 +3142,4 @@
 	}
 
-	printf("-------------- file: FemModel.cpp line: %i\n",__LINE__); 
-	printf("Found %i %i inputs\n",numP0inputs,numP1inputs);
-
 	/*========== Deal with P0 inputs ==========*/
 	int numelementsold = this->elements->NumberOfElements();
@@ -3203,5 +3152,4 @@
 
 	for(int i=0;i<numP0inputs;i++){
-		printf("Dealing with %s (type: %s)\n",EnumToStringx(P0input_enums[i]),EnumToStringx(P0input_interp[i]));
 		GetVectorFromInputsx(&vector,this,P0input_enums[i],ElementSIdEnum);
 
@@ -3210,6 +3158,4 @@
 		xDelete<IssmDouble>(vector);
 	}
-	printf("Printing array to make sure it looks alright\n");
-	printarray(P0inputsold,numelementsold,numP0inputs);
 
 	/*========== Deal with P1 inputs ==========*/
@@ -3219,5 +3165,4 @@
 
 	for(int i=0;i<numP1inputs;i++){
-		printf("Dealing with %s (type: %s)\n",EnumToStringx(P1input_enums[i]),EnumToStringx(P1input_interp[i]));
 		GetVectorFromInputsx(&vector,this,P1input_enums[i],VertexSIdEnum);
 
@@ -3226,6 +3171,4 @@
 		xDelete<IssmDouble>(vector);
 	}
-	printf("Printing array to make sure it looks alright\n");
-	printarray(P1inputsold,numverticesold,numP1inputs);
 	
 	/*Old mesh coordinates*/
@@ -3234,25 +3177,36 @@
 	IssmDouble *Zold		= NULL;
 	int        *Indexold = NULL;
-	int        *Indexnew = NULL;
 	IssmDouble *Xnew     = NULL;
 	IssmDouble *Ynew     = NULL;
 	IssmDouble *Znew		= NULL;
+	IssmDouble* XC_new   = NULL;
+	IssmDouble* YC_new   = NULL;
+	int        *Indexnew = NULL;
 	
 	/*Get the old mesh*/
-	printf("-------------- file: FemModel.cpp line: %i\n",__LINE__); 
 	this->GetMesh(this,&Xold,&Yold,&Zold,&Indexold);
 
 	/*Get the new mesh*/
-	printf("-------------- file: FemModel.cpp line: %i\n",__LINE__); 
 	this->GetMesh(newfemmodel,&Xnew,&Ynew,&Znew,&Indexnew);
-	
-	printf("-------------- file: FemModel.cpp line: %i\n",__LINE__); 
+
+	/*Calculate the center points xc and xy*/
+	int elementswidth = this->GetElementsWidth(); //just tria in this version
+	/*new mesh*/
+	XC_new=xNewZeroInit<IssmDouble>(numelementsnew);
+	YC_new=xNewZeroInit<IssmDouble>(numelementsnew);
+	for(int i=0;i<numelementsnew;i++){
+		for(int j=0;j<elementswidth;j++){
+			int vid = Indexnew[i*elementswidth+j]-1;//Transform to C indexing
+			XC_new[i]+=Xnew[vid];
+			YC_new[i]+=Ynew[vid];
+		}
+		XC_new[i]=XC_new[i]/3;
+		YC_new[i]=YC_new[i]/3;
+	}
+
 	/*Interplate P0 inputs in the new mesh*/
 	InterpFromMeshToMesh2dx(&P0inputsnew,Indexold,Xold,Yold,numverticesold,numelementsold,
 				P0inputsold,numelementsold,numP0inputs,
-				Xnew,Ynew,numelementsnew,NULL);
-	
-	printf("Printing array AFTER INTERP - P0 INPUTS - to make sure it looks alright\n");
-	printarray(P0inputsnew,numelementsnew,numP0inputs);
+				XC_new,YC_new,numelementsnew,NULL);
 	
 	/*Interpolate P1 inputs in the new mesh*/
@@ -3261,66 +3215,128 @@
 				Xnew,Ynew,numverticesnew,NULL);
 
-	printf("Printing array AFTER INTERP - P1 INPUTS - to make sure it looks alright\n");
-	printarray(P1inputsnew,numverticesnew,numP1inputs);
-	
-	_error_("STOP");
-
+	/*Insert P0 inputs into the new elements.*/
+	vector=NULL;
+	for(int i=0;i<numP0inputs;i++){
+		
+		/*Get P0 input vector from the interpolated matrix*/
+		vector=xNew<IssmDouble>(numelementsnew);
+		for(int j=0;j<numelementsnew;j++) vector[j]=P0inputsnew[j*numP0inputs+i];//vector has values for all elements (serial)
+
+		/*Update elements from inputs: */
+		for(int j=0;j<newfemmodel->elements->Size();j++){
+			Element* element=xDynamicCast<Element*>(newfemmodel->elements->GetObjectByOffset(j));
+			switch(P0input_interp[i]){	
+				case P0Enum:
+				case DoubleInputEnum:
+					element->AddInput(new DoubleInput(P0input_enums[i],vector[element->sid]));//sid because newfemmodel has just a partitioning 
+					break;
+				case IntInputEnum: 
+					element->AddInput(new IntInput(P0input_enums[i],vector[element->sid]));//sid because newfemmodel has just a partitioning
+					break;
+				case BoolInputEnum:
+					element->AddInput(new BoolInput(P0input_enums[i],vector[element->sid]));//sid because newfemmodel has just a partitioning
+					break;
+				default:
+					_error_(EnumToStringx(P0input_enums[i])<<" Not supported yet");
+			}
+		}
+
+		xDelete<IssmDouble>(vector);
+	}
+
+	/*Insert P1 inputs into the new elements.*/
+	vector=NULL;
+	for(int i=0;i<numP1inputs;i++){
+
+		/*Get P0 input vector from the interpolated matrix*/
+		vector=xNew<IssmDouble>(numverticesnew);
+		for(int j=0;j<numverticesnew;j++) vector[j]=P1inputsnew[j*numP1inputs+i];//vector has all vertices	(serial)
+
+		/*Update elements from inputs: */
+		InputUpdateFromVectorx(newfemmodel,vector,P1input_enums[i],VertexSIdEnum);//VertexSId because vector is serial in SId indexing
+		
+		xDelete<IssmDouble>(vector);
+	}
+
+	/*Cleanup*/
+	xDelete<IssmDouble>(input_interpolations_serial);
 	xDelete<IssmDouble>(P0inputsold);
 	xDelete<IssmDouble>(P0inputsnew);
+	xDelete<int>(P0input_enums);
+	xDelete<int>(P0input_interp);
 	xDelete<IssmDouble>(P1inputsold);
 	xDelete<IssmDouble>(P1inputsnew);
-
-	_error_("STOP");
-}
-/*}}}*/
-void FemModel::ExecuteRefinement(int &newnumberofvertices,int &newnumberofelements,IssmDouble** newx,IssmDouble** newy,IssmDouble** newz,int** pnewelementslist){/*{{{*/
-
-	int **newelements=NULL;
-	int **newsegments=NULL;	
-	int newnumberofsegments=-1;
-	
-	/*All indexing here is in C type: 0..n-1*/
-	int my_rank=IssmComm::GetRank();
-   const int elementswidth=3;//just 2D mesh, tria elements
+	xDelete<int>(P1input_enums);
+	xDelete<int>(P1input_interp);
+	xDelete<IssmDouble>(Xold);
+	xDelete<IssmDouble>(Yold);
+	xDelete<IssmDouble>(Zold);
+	xDelete<int>(Indexold);
+	xDelete<IssmDouble>(Xnew);
+	xDelete<IssmDouble>(Ynew);
+	xDelete<IssmDouble>(Znew);
+	xDelete<IssmDouble>(XC_new);
+	xDelete<IssmDouble>(YC_new);
+	xDelete<int>(Indexnew);
+
+}
+/*}}}*/
+void FemModel::ExecuteRefinement(int &numberofvertices,int &numberofelements,IssmDouble** px,IssmDouble** py,IssmDouble** pz,int** pelementslist){/*{{{*/
+	
+	/*elements is in Matlab indexing*/
+	
+	int my_rank					 = IssmComm::GetRank();
+	int numberofsegments		 = -1;
+	IssmDouble* vx				 = NULL; //itapopo this is not being used
+	IssmDouble* vy				 = NULL; //itapopo this is not being used
+	IssmDouble* x				 = NULL;
+	IssmDouble* y				 = NULL;
+	IssmDouble* z				 = NULL;
+	int* elementslist			 = NULL;
+	int* segments				 = NULL;
+	IssmDouble* masklevelset = NULL;
+   const int elementswidth  = this->GetElementsWidth();//just 2D mesh, tria elements
 	
 	/*Solutions which will be used to refine the elements*/
-	IssmDouble* vx=NULL; //itapopo this is not being used
-	IssmDouble* vy=NULL; //itapopo this is not being used
-	IssmDouble* masklevelset=NULL;
-	int* newelementslist=NULL;
-
 	this->GetMaskLevelSet(&masklevelset);//itapopo verificar se já existe um método igual a esse
 	
-	int type_process=1; //1: it refines father mesh. See AdaptiveMeshRefinement.h (.cpp)
 	if(my_rank==0){
-		this->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 refinement process.");
-		
-//		printf("   PRINTING COORDINATES... \n");	    
- //		for(int i=0;i<newnumberofvertices;i++) std::cout << "ID: " << i << "\t" << "X: " << newx[i] << "\t" << "Y: " << newy[i] << "\n";
-  //  	printf("   PRINTING ELEMENTS... \n");	   
-  // 	for(int i=0;i<newnumberofelements;i++) std::cout << "El: " << i << "\t" << newelements[i][0] << "\t" << newelements[i][1] << "\t" << newelements[i][2] << "\n";
-	}
-	
-	/*Fill the element list to partitioning*/	
-	if(my_rank==0){
-		newelementslist=xNew<int>(newnumberofelements*elementswidth);
-		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
-			}
-		}	
-		(*pnewelementslist)=newelementslist;
-	}
-
-	/*Cleanup masklevetset*/
+		int type_process=1; //1: it refines father mesh. See AdaptiveMeshRefinement.h (.cpp)
+		this->amr->ExecuteRefinement(type_process,vx,vy,masklevelset,
+												numberofvertices,numberofelements,numberofsegments,&x,&y,&z,&elementslist,&segments);
+		if(numberofvertices<=0 || numberofelements<=0 /*|| newnumberofsegments<=0*/) _error_("Error in the refinement process.");
+	}
+	else{
+		x=xNew<IssmDouble>(numberofvertices);
+		y=xNew<IssmDouble>(numberofvertices);
+		z=xNew<IssmDouble>(numberofvertices);
+		elementslist=xNew<int>(numberofelements*this->GetElementsWidth());
+	}
+
+	/*Send new mesh to others CPU*/
+	ISSM_MPI_Bcast(&numberofvertices,1,ISSM_MPI_INT,0,IssmComm::GetComm());
+	ISSM_MPI_Bcast(&numberofelements,1,ISSM_MPI_INT,0,IssmComm::GetComm());
+	ISSM_MPI_Bcast(x,numberofvertices,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());	
+	ISSM_MPI_Bcast(y,numberofvertices,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());	
+	ISSM_MPI_Bcast(z,numberofvertices,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());	
+	ISSM_MPI_Bcast(elementslist,numberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm());	
+
+	/*Assign the pointers*/	
+	(*pelementslist) = elementslist; //Matlab indexing
+	(*px)				  = x;
+	(*py)				  = y;
+	(*pz)				  = z;
+
+	/*Cleanup*/
+	if(segments) xDelete<int>(segments);
 	xDelete<IssmDouble>(masklevelset);
+
 }
 /*}}}*/
 void FemModel::GetMaskLevelSet(IssmDouble **pmasklevelset){/*{{{*/
 
-	int elementswidth=3;//just 2D mesh, tria elements
-	int numberofelements=this->elements->NumberOfElements();
-	int numberofvertices=this->vertices->NumberOfVertices();
+	int elementswidth		= this->GetElementsWidth();//just 2D mesh, tria elements
+	int numberofelements = this->elements->NumberOfElements();
+	int numberofvertices = this->vertices->NumberOfVertices();
 
 	IssmDouble* elementlevelset=xNew<IssmDouble>(elementswidth);
@@ -3352,4 +3368,6 @@
 void FemModel::CreateVertices(int newnumberofvertices,int newnumberofelements,int elementswidth,int* newelementslist,int* my_vertices,IssmDouble* newx,IssmDouble* newy,IssmDouble* newz,Vertices* vertices){/*{{{*/
 
+	/*newelementslist is in Matlab indexing*/
+	
 	/*Creating connectivity table*/
 	int* connectivity=NULL;
@@ -3358,7 +3376,7 @@
 	for (int i=0;i<newnumberofelements;i++){
 		for (int j=0;j<elementswidth;j++){
-			int vertexid = newelementslist[elementswidth*i+j];//C indexing
-			_assert_(vertexid>-1 && vertexid<newnumberofvertices);
-			connectivity[vertexid]+=1;
+			int vertexid = newelementslist[elementswidth*i+j];
+			_assert_(vertexid>0 && vertexid-1<newnumberofvertices);//Matlab indexing
+			connectivity[vertexid-1]+=1;//Matlab to C indexing
 		}
 	}	
@@ -3386,4 +3404,6 @@
 /*}}}*/
 void FemModel::CreateElements(int newnumberofelements,int elementswidth,int* newelementslist,bool* my_elements,Elements* elements){/*{{{*/
+
+	/*newlementslist is in Matlab indexing*/
 
 	for(int i=0;i<newnumberofelements;i++){
@@ -3410,5 +3430,5 @@
 			/*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])+1;//this Hook wants Matlab indexing	
+			for(int j=0;j<elementswidth;j++)	vertex_ids[j]=reCast<int>(newelementslist[elementswidth*i+j]);//this Hook wants Matlab indexing	
 			/*Setting the hooks*/
 			newtria->numanalyses =this->nummodels;
@@ -3426,6 +3446,4 @@
 	}
 
-	//itapopo there is this line in CreateElementsVerticesAndMaterials.
-	//elements->InputDuplicate(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum);
 }
 /*}}}*/
@@ -3441,4 +3459,5 @@
 	/*Add new constant material property to materials, at the end: */
 	Matpar *newmatpar=static_cast<Matpar*>(this->materials->GetObjectByOffset(this->materials->Size()-1)->copy());
+	newmatpar->mid = newnumberofelements+1;//itapopo testando
 	materials->AddObject(newmatpar);//put it at the end of the materials	    
 
@@ -3448,5 +3467,4 @@
 
 	int lid=0;
-//itapopo use the GetMesh!!
 	for(int j=0;j<newnumberofvertices;j++){
 		if(my_vertices[j]){				
@@ -3497,5 +3515,5 @@
 	
 	int numberofvertices, numberofelements;
-	int elementswidth = 3; //itapopo just 2D mesh in this version (just tria elements)
+	int elementswidth = this->GetElementsWidth(); // just 2D mesh in this version (just tria elements)
 	IssmDouble *x		= NULL;
 	IssmDouble *y		= NULL;
@@ -3538,5 +3556,4 @@
 	if(numberofelements*elementswidth<0) _error_("numberofelements negative.");
 	for(int i=0;i<numberofelements;i++){
-		std::cout << (int)id1[i]+1 << " " << (int)id2[i]+1 <<" " << (int)id3[i]+1 <<  std::endl;
 		elementslist[elementswidth*i+0] = (int)id1[i]+1; //InterpMesh wants Matlab indexing
 		elementslist[elementswidth*i+1] = (int)id2[i]+1; //InterpMesh wants Matlab indexing
@@ -3548,55 +3565,30 @@
 	*py				= y;
 	*pz				= z;
-	*pelementslist = elementslist;
-}
-/*}}}*/
-void FemModel::CreateConstraints(int newnumberofvertices,int newnumberofelements,int* newelementslist,IssmDouble* newx,IssmDouble* newy,int* my_vertices,Constraints* constraints){/*{{{*/
+	*pelementslist = elementslist; //Matlab indexing. InterMesh uses this type.
+
+	/*Cleanup*/
+	xDelete<int>(elem_vertices);
+	xDelete<IssmDouble>(id1);
+	xDelete<IssmDouble>(id2);
+	xDelete<IssmDouble>(id3);
+	delete vid1;
+	delete vid2;
+	delete vid3;
+}
+/*}}}*/
+void FemModel::CreateConstraints(int newnumberofvertices,int newnumberofelements,IssmDouble* newx,IssmDouble* newy,int* my_vertices,Constraints* constraints){/*{{{*/
 
 	/*Get x and y of the mesh i-1*/
-	int numberofvertices, numberofelements;
-	numberofvertices = this->vertices->NumberOfVertices();
-	numberofelements = 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 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, get vertices*/
-   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();
-	
-	/*Construct elements list (mesh i-1)*/
-	int* elementslist=NULL;
-	elementslist=xNew<int>(numberofelements*elementswidth);
-	for(int i=0;i<numberofelements;i++){ //itapopo os elementos poderão sair do ExecuteRefinement nesse formato
-		elementslist[elementswidth*i+0] = (int)id1[i]+1; //InterpMesh wants Matlab indexing
-		elementslist[elementswidth*i+1] = (int)id2[i]+1; //InterpMesh wants Matlab indexing
-		elementslist[elementswidth*i+2] = (int)id3[i]+1; //InterpMesh wants Matlab indexinf
-	}	
-
+	int my_rank						= IssmComm::GetRank();
+	int numberofvertices			= this->vertices->NumberOfVertices();
+	int numberofelements			= this->elements->NumberOfElements();
+	IssmDouble *x					= NULL;
+	IssmDouble *y					= NULL;
+	IssmDouble *z					= NULL;
+	int        *elementslist	= NULL;
+
+	/*elementslist is in Matlab indexing*/
+	this->GetMesh(this,&x,&y,&z,&elementslist);
+	
 	/*itapopo ATTENTION: JUST SPCVX AND SPCVY TO TEST!!!*/
 	/*OTHERS CONSTRAINTS MUST BE IMPLEMENTED!!!*/
@@ -3645,30 +3637,24 @@
 
 	/*Serialize*/
-	spcvx=vspcvx->ToMPISerial();
-	spcvy=vspcvy->ToMPISerial();
-	spcvxflag=vspcvxflag->ToMPISerial();
-	spcvyflag=vspcvyflag->ToMPISerial();
-	/*Free the data*/
-	delete vspcvx;
-	delete vspcvy;	
-	delete vspcvxflag;
-	delete vspcvyflag;
-	
-	IssmDouble *newspcvx=NULL;
-	IssmDouble *newspcvy=NULL;
-	IssmDouble *newspcvxflag=NULL;
-	IssmDouble *newspcvyflag=NULL;
-	int nods_data=numberofvertices;
-	int nels_data=numberofelements;
-	int M_data=numberofvertices;
-	int N_data=1;
-	int N_interp=newnumberofvertices;
-	Options *options=new Options();
+	spcvx		 = vspcvx->ToMPISerial();
+	spcvy		 = vspcvy->ToMPISerial();
+	spcvxflag = vspcvxflag->ToMPISerial();
+	spcvyflag = vspcvyflag->ToMPISerial();
+	
+	IssmDouble *newspcvx		 = NULL;
+	IssmDouble *newspcvy		 = NULL;
+	IssmDouble *newspcvxflag = NULL;
+	IssmDouble *newspcvyflag = NULL;
+	int nods_data				 = numberofvertices;
+	int nels_data				 = numberofelements;
+	int M_data					 = numberofvertices;
+	int N_data					 = 1;
+	int N_interp				 = newnumberofvertices;
 
   /*Interpolate spcvx and spcvy in the new mesh*/
-	InterpFromMeshToMesh2dx(&newspcvx,elementslist,x,y,nods_data,nels_data,spcvx,M_data,N_data,newx,newy,N_interp,options);
-	InterpFromMeshToMesh2dx(&newspcvy,elementslist,x,y,nods_data,nels_data,spcvy,M_data,N_data,newx,newy,N_interp,options);
-	InterpFromMeshToMesh2dx(&newspcvxflag,elementslist,x,y,nods_data,nels_data,spcvxflag,M_data,N_data,newx,newy,N_interp,options);
-	InterpFromMeshToMesh2dx(&newspcvyflag,elementslist,x,y,nods_data,nels_data,spcvyflag,M_data,N_data,newx,newy,N_interp,options);
+	InterpFromMeshToMesh2dx(&newspcvx,elementslist,x,y,nods_data,nels_data,spcvx,M_data,N_data,newx,newy,N_interp,NULL);
+	InterpFromMeshToMesh2dx(&newspcvy,elementslist,x,y,nods_data,nels_data,spcvy,M_data,N_data,newx,newy,N_interp,NULL);
+	InterpFromMeshToMesh2dx(&newspcvxflag,elementslist,x,y,nods_data,nels_data,spcvxflag,M_data,N_data,newx,newy,N_interp,NULL);
+	InterpFromMeshToMesh2dx(&newspcvyflag,elementslist,x,y,nods_data,nels_data,spcvyflag,M_data,N_data,newx,newy,N_interp,NULL);
 	
 	int nodecounter			= 0; //itapopo deve começar pelo primeiro nó do StressbalanceAnalysis
@@ -3696,9 +3682,29 @@
 			count++;
 		}
-
-	}
+	}
+
+	/*Cleanup*/
+	xDelete<IssmDouble>(x);
+	xDelete<IssmDouble>(y);
+	xDelete<IssmDouble>(z);
+	xDelete<int>(elementslist);
+	xDelete<IssmDouble>(spcvx);
+	xDelete<IssmDouble>(spcvy);
+	xDelete<IssmDouble>(spcvxflag);
+	xDelete<IssmDouble>(spcvyflag);
+	xDelete<IssmDouble>(newspcvx);
+	xDelete<IssmDouble>(newspcvy);
+	xDelete<IssmDouble>(newspcvxflag);
+	xDelete<IssmDouble>(newspcvyflag);
+	delete vspcvx;
+	delete vspcvy;	
+	delete vspcvxflag;
+	delete vspcvyflag;
+
 }
 /*}}}*/
 void FemModel::UpdateElements(int newnumberofelements,int* newelementslist,bool* my_elements,int nodecounter,int analysis_counter,Elements* newelements){/*{{{*/
+
+	/*newelementslist is in Matlab indexing*/
 
 	/*Update elements, set hnode.
@@ -3712,7 +3718,7 @@
 			int numnodes=3;
          int* tria_node_ids=xNew<int>(numnodes);
-         tria_node_ids[0]=nodecounter+newelementslist[3*iel+0]+1; //matlab indexing
-         tria_node_ids[1]=nodecounter+newelementslist[3*iel+1]+1; //matlab indexing
-         tria_node_ids[2]=nodecounter+newelementslist[3*iel+2]+1; //matlab indexing
+         tria_node_ids[0]=nodecounter+newelementslist[3*iel+0]; //matlab indexing
+         tria_node_ids[1]=nodecounter+newelementslist[3*iel+1]; //matlab indexing
+         tria_node_ids[2]=nodecounter+newelementslist[3*iel+2]; //matlab indexing
 			tria->SetHookNodes(tria_node_ids,numnodes,analysis_counter); tria->nodes=NULL;
    		xDelete<int>(tria_node_ids);
@@ -3725,13 +3731,16 @@
 void FemModel::ElementsAndVerticesPartitioning(int& newnumberofvertices,int& newnumberofelements,int& elementswidth,int* newelementslist,bool** pmy_elements,int** pmy_vertices){/*{{{*/
 
-	int *epart=NULL; //element partitioning.
-	int *npart=NULL; //node partitioning.
-	int edgecut=1;
-	int numflag=0;
-	int etype=1;
-	int my_rank = IssmComm::GetRank();
-	int numprocs=IssmComm::GetSize();
-	bool *my_elements=NULL;
-	int *my_vertices=NULL;
+	/*newelementslist come in Matlab indexing*/
+
+	int *epart			= NULL; //element partitioning.
+	int *npart			= NULL; //node partitioning.
+	int *index 			= NULL; //elements in C indexing
+	int edgecut			= 1;
+	int numflag			= 0;
+	int etype			= 1;
+	int my_rank			= IssmComm::GetRank();
+	int numprocs		= IssmComm::GetSize();
+	bool *my_elements = NULL;
+	int *my_vertices  = NULL;
 	
 	_assert_(newnumberofvertices>0); 
@@ -3739,9 +3748,16 @@
 	epart=xNew<int>(newnumberofelements);
 	npart=xNew<int>(newnumberofvertices);
+   index=xNew<int>(elementswidth*newnumberofelements);
+   
+	for (int i=0;i<newnumberofelements;i++){
+   	for (int j=0;j<elementswidth;j++){
+      	*(index+elementswidth*i+j)=(*(newelementslist+elementswidth*i+j))-1; //-1 for C indexing in Metis
+      }
+   }
 
 	/*Partition using Metis:*/
 	if (numprocs>1){
 #ifdef _HAVE_METIS_
-		METIS_PartMeshNodalPatch(&newnumberofelements,&newnumberofvertices, newelementslist, &etype, &numflag, &numprocs, &edgecut, epart, npart);
+		METIS_PartMeshNodalPatch(&newnumberofelements,&newnumberofvertices, index, &etype, &numflag, &numprocs, &edgecut, epart, npart);
 #else
 		_error_("metis has not beed installed. Cannot run with more than 1 cpu");
@@ -3770,17 +3786,18 @@
 			 will hold which vertices belong to this partition*/
 			for(int j=0;j<elementswidth;j++){
-				_assert_(newelementslist[elementswidth*i+j]<newnumberofvertices); 
-				my_vertices[newelementslist[elementswidth*i+j]]=1;
+				_assert_(newelementslist[elementswidth*i+j]-1<newnumberofvertices);//newelementslist is in Matlab indexing 
+				my_vertices[newelementslist[elementswidth*i+j]-1]=1;//newelementslist is in Matlab indexing
 			}
 		}
 	}
+
+	/*Assign output pointers:*/
+	*pmy_elements=my_elements;
+	*pmy_vertices=my_vertices;
 
 	/*Free ressources:*/
 	xDelete<int>(epart);
 	xDelete<int>(npart);	    
-
-	/*Assign output pointers:*/
-	*pmy_elements=my_elements;
-	*pmy_vertices=my_vertices;
+	xDelete<int>(index);
 
 }
Index: /issm/trunk-jpl/src/c/classes/FemModel.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.h	(revision 21527)
+++ /issm/trunk-jpl/src/c/classes/FemModel.h	(revision 21528)
@@ -150,5 +150,6 @@
 		FemModel* ReMesh(void);
 		void GetMesh(FemModel* femmodel,IssmDouble** px, IssmDouble** py, IssmDouble** pz, int** pelementslist);
-		void ExecuteRefinement(int &newnumberofvertices,int &newnumberofelements,IssmDouble** newx,IssmDouble** newy,IssmDouble** newz,int** newelementslist);
+		int GetElementsWidth(){return 3;};//just tria elements in this first version
+		void ExecuteRefinement(int &numberofvertices,int &numberofelements,IssmDouble** px,IssmDouble** py,IssmDouble** pz,int** pelementslist);
 		void GetMaskLevelSet(IssmDouble** pmasklevelset);
 		void CreateVertices(int newnumberofvertices,int newnumberofelements,int elementswidth,int* newelementslist,int* my_vertices,IssmDouble* newx,IssmDouble* newy,IssmDouble* newz,Vertices* vertices);
@@ -156,5 +157,5 @@
 		void CreateMaterials(int newnumberofelements,bool* my_elements,Materials* materials);
 		void CreateNodes(int newnumberofvertices,int* my_vertices,int nodecounter,int analysis_enum,Nodes* nodes);
-		void CreateConstraints(int newnumberofvertices,int newnumberofelements,int* newelementslist,IssmDouble* newx,IssmDouble* newy,int* my_vertices,Constraints* constraints);
+		void CreateConstraints(int newnumberofvertices,int newnumberofelements,IssmDouble* newx,IssmDouble* newy,int* my_vertices,Constraints* constraints);
 		void InterpolateInputs(FemModel* femmodel);
 		void UpdateElements(int newnumberofelements,int* newelementslist,bool* my_elements,int nodecounter,int analysis_counter,Elements* newelements);