Changeset 21528

Timestamp:

02/08/17 17:13:30 (8 years ago)

Author:

tsantos

Message:

CHG: changes in AMR methods (in FemModel) and in AdaptiveMeshRefinement class. Cleanup done and parallel tested (ReMesh method, simple example). AMR not working yet.

Location:

issm/trunk-jpl/src/c/classes

Files:

• AdaptiveMeshRefinement.cpp (modified) (16 diffs)
• AdaptiveMeshRefinement.h (modified) (3 diffs)
• Elements/Element.cpp (modified) (1 diff)
• FemModel.cpp (modified) (30 diffs)
• FemModel.h (modified) (2 diffs)

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

@@ -3 +3 @@
  */
 
- #ifdef HAVE_CONFIG_H
+#ifdef HAVE_CONFIG_H
     #include <config.h>
 #else
@@ -52 +52 @@
          this->hmax=-1;
          this->elementswidth=-1;
+         gRefDBase.clear();
 }
 /*}}}*/
@@ -128 +129 @@
 /*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 +176 @@
     
     /*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 +260 @@
 }
 /*}}}*/
-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 +269 @@
     
     /* 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 +291 @@
     
     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 +308 @@
     
     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 +326 @@
         
         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 +388 @@
     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 +444 @@
 }
 /*}}}*/
-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 +460 @@
 
     /*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 +525 @@
 }
 /*}}}*/
-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 +566 @@
     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 +573 @@
     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 +580 @@
         
     //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();
+        }
+    }
+    
+}
+/*}}}*/

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

@@ -10 +10 @@
 
 /*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 +59 @@
    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 +79 @@
    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
 };
 

issm/trunk-jpl/src/c/classes/Elements/Element.cpp

@@ -1578 +1578 @@
                  name==HydrologyHeadOldEnum ||          
                                 name==StressbalanceConvergenceNumStepsEnum || 
-                                name==MeshVertexonbaseEnum 
+                                name==MeshVertexonbaseEnum || 
+                                name==FrictionPEnum ||
+                                name==FrictionQEnum ||
+                                name==LoadingforceXEnum ||
+                                name==LoadingforceYEnum
 
                                 ) {

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

@@ -82 +82 @@
         #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 +2922 @@
 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 +2951 @@
         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 +2957 @@
 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 +3001 @@
         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 +3033 @@
         /*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 +3061 @@
         }
 
-        /*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 +3142 @@
         }
 
-        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 +3152 @@
 
         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 +3158 @@
                 xDelete<IssmDouble>(vector);
         }
-        printf("Printing array to make sure it looks alright\n");
-        printarray(P0inputsold,numelementsold,numP0inputs);
 
         /*========== Deal with P1 inputs ==========*/
@@ -3219 +3165 @@
 
         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 +3171 @@
                 xDelete<IssmDouble>(vector);
         }
-        printf("Printing array to make sure it looks alright\n");
-        printarray(P1inputsold,numverticesold,numP1inputs);
         
         /*Old mesh coordinates*/
@@ -3234 +3177 @@
         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 +3215 @@
                                 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 +3368 @@
 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 +3376 @@
         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 +3404 @@
 /*}}}*/
 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 +3430 @@
                         /*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 +3446 @@
         }
 
-        //itapopo there is this line in CreateElementsVerticesAndMaterials.
-        //elements->InputDuplicate(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum);
 }
 /*}}}*/
@@ -3441 +3459 @@
         /*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 +3467 @@
 
         int lid=0;
-//itapopo use the GetMesh!!
         for(int j=0;j<newnumberofvertices;j++){
                 if(my_vertices[j]){                             
@@ -3497 +3515 @@
         
         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 +3556 @@
         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 +3565 @@
         *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 +3637 @@
 
         /*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 +3682 @@
                         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 +3718 @@
                         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 +3731 @@
 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 +3748 @@
         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 +3786 @@
                          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);
 
 }

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

@@ -150 +150 @@
                 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 +157 @@
                 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);