Changeset 21484

Timestamp:

01/11/17 14:56:54 (8 years ago)

Author:

tsantos

Message:

NEW: adding AMR functionality, not working yet but structure in place

Location:

issm/trunk-jpl

Files:

• m4/issm_options.m4 (modified) (1 diff)
• src/c/Makefile.am (modified) (2 diffs)
• src/c/classes/AdaptiveMeshRefinement.cpp (added)
• src/c/classes/AdaptiveMeshRefinement.h (added)
• src/c/classes/FemModel.cpp (modified) (2 diffs)
• src/c/classes/FemModel.h (modified) (1 diff)
• src/c/classes/Materials/Matice.cpp (modified) (2 diffs)
• src/c/classes/Materials/Matice.h (modified) (1 diff)
• src/c/classes/Node.cpp (modified) (1 diff)
• src/c/classes/Node.h (modified) (1 diff)

issm/trunk-jpl/m4/issm_options.m4

@@ -1896 +1896 @@
         AC_MSG_RESULT($HAVE_KRIGING)
         dnl }}}
+        dnl with-amr{{{
+        AC_ARG_WITH([amr],
+                AS_HELP_STRING([--with-amr = YES],[compile with Adaptive Mesh Refinment capability (default is no)]),
+                [AMR=$withval],[AMR=no]) 
+        AC_MSG_CHECKING(for AMR capability compilation)
+
+        HAVE_AMR=no
+        if test "x$AMR" = "xyes"; then
+                HAVE_AMR=yes
+                AC_DEFINE([_HAVE_AMR_],[1],[with amr capability])
+        fi
+        AM_CONDITIONAL([AMR], [test x$HAVE_AMR = xyes])
+        AC_MSG_RESULT($HAVE_AMR)
+        dnl }}}
         AX_ANALYSES_SELECTION
 

issm/trunk-jpl/src/c/Makefile.am

@@ -530 +530 @@
                                   ./kml/KMLFileReadUtils.cpp
 #}}}
+#AMR sources  {{{
+amr_sources = ./classes/AdaptiveMeshRefinement.cpp
+#}}}
 #Modules sources{{{
 modules_sources= ./shared/Threads/LaunchThread.cpp\
@@ -613 +616 @@
 libISSMModules_la_SOURCES += $(kml_sources)
 endif
+if AMR
+libISSMModules_la_SOURCES += $(amr_sources)
+endif
 libISSMModules_la_CXXFLAGS = $(ALLCXXFLAGS)
 if !WINDOWS

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

@@ -80 +80 @@
         /*Save communicator in the parameters dataset: */
         this->parameters->AddObject(new GenericParam<ISSM_MPI_Comm>(incomm,FemModelCommEnum));
+
+        #ifdef _HAVE_AMR_
+        this->InitializeAdaptiveRefinement();
+        #endif
 
         /*Free resources */
@@ -2911 +2915 @@
 }/*}}}*/
 #endif
+
+#ifdef _HAVE_AMR_
+void FemModel::InitializeAdaptiveRefinement(void){/*{{{*/
+        
+        int my_rank=IssmComm::GetRank();
+        
+        /*Initialize AMR*/
+        //AdaptiveMeshRefinement *AMR = new AdaptiveMeshRefinement();
+        
+        /*Set max level of refinement*/
+        int hmax = 2;
+        //AMR->SetHMax(hmax);
+        
+        /*Get initial mesh*/
+        /*Get total number of elements and total numbers of elements in the initial mesh*/
+        int numberofvertices, numberofelements, numberofsegments;
+        numberofvertices = this->vertices->NumberOfVertices();
+        numberofelements = this->elements->NumberOfElements();
+        numberofsegments = -1; //used on matlab
+
+        /*Get vertices coordinates*/
+        IssmDouble *x = NULL;
+        IssmDouble *y = NULL;
+        IssmDouble *z = NULL;
+        VertexCoordinatesx(&x, &y, &z, this->vertices,false) ;
+
+        /*Get element vertices*/
+        int elementswidth = 3; //just 2D mesh in this version (just tria elements)
+        int* elem_vertices=xNew<int>(elementswidth);
+        Vector<IssmDouble>* vid1= new Vector<IssmDouble>(numberofelements);
+        Vector<IssmDouble>* vid2= new Vector<IssmDouble>(numberofelements);
+        Vector<IssmDouble>* vid3= new Vector<IssmDouble>(numberofelements);
+
+        /*Go through elements, and for each element, object*/
+    for(int i=0;i<this->elements->Size();i++){
+        Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
+        element->GetVerticesSidList(elem_vertices);
+        vid1->SetValue(element->sid,elem_vertices[0],INS_VAL);
+        vid2->SetValue(element->sid,elem_vertices[1],INS_VAL);
+        vid3->SetValue(element->sid,elem_vertices[2],INS_VAL);
+    }
+                
+        /*Assemble*/
+    vid1->Assemble();
+    vid2->Assemble();
+    vid3->Assemble();
+
+    /*Serialize*/
+        IssmDouble *id1 = vid1->ToMPISerial();
+   IssmDouble *id2 = vid2->ToMPISerial();
+        IssmDouble *id3 = vid3->ToMPISerial();
+        
+        int **elementsptr;
+        elementsptr = new int*[numberofelements];
+    for(int i=0;i<numberofelements;i++){
+        elementsptr[i] = new int[elementswidth];
+        elementsptr[i][0] = (int)id1[i];
+        elementsptr[i][1] = (int)id2[i];
+        elementsptr[i][2] = (int)id3[i];
+    }
+
+        //if(my_rank==0){
+        //      _printf_("   PRINTING COORDINATES... \n");          
+   //           for(long i=0;i<numberofvertices;i++) _printf_("ID: " << i << "\t" << "X: " << x[i] << "\t" << "Y: " << y[i] << "\n");
+   //   _printf_("   PRINTING ELEMENTS... \n");    
+   //   for(int i=0;i<numberofelements;i++) _printf_("El: " << i << "\t" << elementsptr[i][0] << "\t" << elementsptr[i][1] << "\t" << elementsptr[i][2] << "\n");
+        //}
+
+    //AMR->CreateInitialMesh(numberofvertices, numberofelements, 0, elementswidth, x, y, z, elementsptr, NULL);
+
+        /*Free the vectors*/
+        delete x;
+        delete y;
+        delete z;
+        delete vid1;
+        delete vid2;
+        delete vid3;
+        delete id1;
+        delete id2;
+        delete id3;
+        for(int i=0;i<numberofelements;i++) delete elementsptr[i];
+   if(elementsptr) delete elementsptr;
+        xDelete<int>(elem_vertices); 
+}
+/*}}}*/
+FemModel* FemModel::ReMesh(void){/*{{{*/
+        
+        /*All indexing here is in C type: 0..n-1*/
+
+        //////// to test! using the current femmodel
+        int numberofvertices2, numberofelements2;
+        numberofvertices2 = this->vertices->NumberOfVertices();
+        numberofelements2 = this->elements->NumberOfElements();
+
+        /*Get vertices coordinates*/
+        IssmDouble *x = NULL;
+        IssmDouble *y = NULL;
+        IssmDouble *z = NULL;
+        VertexCoordinatesx(&x, &y, &z, this->vertices,false) ;
+
+        /*Get element vertices*/
+        int elementswidth2 = 3; //just 2D mesh in this version (just tria elements)
+        int* elem_vertices2=xNew<int>(elementswidth2);
+        Vector<IssmDouble>* vid1= new Vector<IssmDouble>(numberofelements2);
+        Vector<IssmDouble>* vid2= new Vector<IssmDouble>(numberofelements2);
+        Vector<IssmDouble>* vid3= new Vector<IssmDouble>(numberofelements2);
+
+        /*Go through elements, and for each element, object, report it cpu:*/
+    for(int i=0;i<this->elements->Size();i++){
+        Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
+        element->GetVerticesSidList(elem_vertices2);
+        vid1->SetValue(element->sid,elem_vertices2[0],INS_VAL);
+        vid2->SetValue(element->sid,elem_vertices2[1],INS_VAL);
+        vid3->SetValue(element->sid,elem_vertices2[2],INS_VAL);
+    }
+                
+        /*Assemble*/
+    vid1->Assemble();
+    vid2->Assemble();
+    vid3->Assemble();
+
+    /*Serialize*/
+        IssmDouble *id1 = vid1->ToMPISerial();
+    IssmDouble *id2 = vid2->ToMPISerial();
+        IssmDouble *id3 = vid3->ToMPISerial();
+
+        delete vid1;
+        delete vid2;
+        delete vid3;
+        
+        int **elementsptr;
+        elementsptr = new int*[numberofelements2];
+    for(int i=0;i<numberofelements2;i++){
+        elementsptr[i] = new int[elementswidth2];
+        elementsptr[i][0] = (int)id1[i];
+        elementsptr[i][1] = (int)id2[i];
+        elementsptr[i][2] = (int)id3[i];
+    }
+
+   delete id1;
+        delete id2;
+        delete id3;
+        //////////////////////////////////////////////
+
+        int my_rank=IssmComm::GetRank();
+
+        /*Solutions which will be used to refine the elements*/
+        double *vx = NULL;
+        double *vy = NULL;
+        double *masklevelset = NULL;
+
+        int elementswidth = 3;//just 2D mesh, tria elements
+        int numberofelements = this->elements->NumberOfElements();
+        int numberofvertices = this->vertices->NumberOfVertices();
+
+        IssmDouble* elementlevelset = xNew<IssmDouble>(elementswidth);
+        int* elem_vertices = xNew<int>(elementswidth);
+        Vector<IssmDouble>* vmasklevelset = new Vector<IssmDouble>(numberofvertices);
+
+        for(int i=0;i<this->elements->Size();i++){
+                Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
+                element->GetInputListOnVertices(elementlevelset,MaskIceLevelsetEnum);
+                element->GetVerticesSidList(elem_vertices);
+                vmasklevelset->SetValue(elem_vertices[0],elementlevelset[0],INS_VAL);
+        vmasklevelset->SetValue(elem_vertices[1],elementlevelset[1],INS_VAL);
+        vmasklevelset->SetValue(elem_vertices[2],elementlevelset[2],INS_VAL);
+        }
+
+        /*Assemble*/
+        vmasklevelset->Assemble();
+
+        /*Serialize*/
+        masklevelset = vmasklevelset->ToMPISerial();
+
+        xDelete<IssmDouble>(elementlevelset);
+        xDelete<int>(elem_vertices);
+        delete vmasklevelset;
+
+        //_printf_("   PRINTING MASKLEVELSET... \n");       
+        //if(my_rank==0){
+        //      for(int i=0;i<numberofvertices;i++) _printf_("vertex: " << i << "\t" << masklevelset[i] << "\n");
+        //}
+
+        /*Refine the mesh*/
+        double *newx;
+        double *newy;
+        double *newz;
+        int **newelements;
+        int **newsegments;
+        int newnumberofvertices, newnumberofelements, newnumberofsegments;
+        int type_process=1;
+        //AMR->ExecuteRefinement(type_process,vx,vy,masklevelset,newnumberofvertices,newnumberofelements,newnumberofsegments,&newx,&newy,&newz,&newelements,&newsegments);
+        
+        //if(newnumberofvertices<=0 || newnumberofelements<=0 || newnumberofsegments=<0) _error_("Error in the mesh refinement process.");
+
+        delete masklevelset;
+
+        ///////////////////////// to test using the current femmodel
+        newx = x;
+        newy = y;
+        newz = z;
+        newelements = elementsptr;
+        newnumberofvertices = numberofvertices2;
+        newnumberofelements = numberofelements2;
+
+        //if(my_rank==0){
+        //      _printf_("   PRINTING COORDINATES... \n");          
+   //           for(int i=0;i<newnumberofvertices;i++) _printf_("ID: " << i << "\t" << "X: " << newx[i] << "\t" << "Y: " << newy[i] << "\n");
+   //   _printf_("   PRINTING ELEMENTS... \n");    
+   //   for(int i=0;i<newnumberofelements;i++) _printf_("El: " << i << "\t" << newelements[i][0] << "\t" << newelements[i][1] << "\t" << newelements[i][2] << "\n");
+        //}
+        /////////////////////////
+
+        /*Partitioning the new mesh*/
+        int        *epart                       = NULL; //element partitioning.
+        int        *npart                       = NULL; //node partitioning.
+        int        *newelementslist     = NULL;
+        int             edgecut                         = 1;
+        int             numflag                         = 0;
+        int             etype                           = 1;
+        int             numprocs                        = IssmComm::GetSize();
+        bool            *my_elements            = NULL;
+        int             *my_vertices            = NULL;
+        
+        epart                                                   = xNew<int>(newnumberofelements);
+        npart                                                   = xNew<int>(newnumberofvertices);
+        newelementslist                                 = xNew<int>(newnumberofelements*elementswidth);
+
+        /*Fill the element list to partitioning*/
+        for(int i=0;i<newnumberofelements;i++){ //itapopo os elementos poderão sair do ExecuteRefinement nesse formato
+                for(int j=0;j<elementswidth;j++){
+                        newelementslist[elementswidth*i+j] = newelements[i][j]; //C indexing
+                }
+        }
+
+        //if(my_rank==0){
+        //      _printf_("   PRINTING ELEMENTS in ELEMENTLIST... \n");     
+   //   for(int i=0;i<newnumberofelements;i++) _printf_("El: " << i << "\t" << newelementslist[i*elementswidth+0] << "\t" << newelementslist[i*elementswidth+1] << "\t" << newelementslist[i*elementswidth+2] << "\n");
+    //}
+
+        /*Partition using Metis:*/
+        if (numprocs>1){
+#ifdef _HAVE_METIS_
+                METIS_PartMeshNodalPatch(&newnumberofelements,&newnumberofvertices, newelementslist, &etype, &numflag, &numprocs, &edgecut, epart, npart);
+#else
+                _error_("metis has not beed installed. Cannot run with more than 1 cpu");
+#endif
+        }
+        else if (numprocs==1){
+                /*METIS does not know how to deal with one cpu only!*/
+                for (int i=0;i<newnumberofelements;i++) epart[i]=0;
+                for (int i=0;i<newnumberofvertices;i++) npart[i]=0;
+        }
+        else _error_("At least one processor is required");
+
+        //_printf_("   PRINTING AFTER MESTIS... \n");       
+
+        my_vertices=xNew<int>(newnumberofvertices);
+        my_elements=xNew<bool>(newnumberofelements);
+        for(int i=0;i<newnumberofvertices;i++) my_vertices[i] = 0;
+        for(int i=0;i<newnumberofelements;i++) my_elements[i] = false;
+
+        /*Start figuring out, out of the partition, which elements belong to this cpu: */
+        for(int i=0;i<newnumberofelements;i++){
+
+                /*!All elements have been partitioned above, only deal with elements for this cpu: */
+                if(my_rank==epart[i]){ 
+                        my_elements[i]=true;
+                        /*Now that we are here, we can also start building the list of vertices belonging to this cpu partition: we use 
+                         *the  element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing) 
+                         into the vertices coordinates. If we start plugging 1 into my_vertices for each index[n][i] (i=0:2), then my_vertices 
+                         will hold which vertices belong to this partition*/
+                        for(int j=0;j<elementswidth;j++){
+                                my_vertices[newelementslist[elementswidth*i+j]]=1;
+                        }
+                }
+        }
+
+        xDelete<int>(epart);
+        xDelete<int>(npart);
+
+        //_printf_("   PRINTING AFTER MY_VERTICES... \n");          
+
+        /*Creating new femmodel with new mesh*/
+        FemModel* output = NULL;
+        int       analysis_type;
+
+        output = new FemModel(*this);
+
+        //_printf_("   PRINTING AFTER FEMMODEL... \n");     
+
+        /*Copy basic attributes:*/
+        output->nummodels = this->nummodels;
+        output->solution_type = this->solution_type;
+        output->analysis_counter = this->analysis_counter;
+
+        //_printf_("   PRINTING AFTER BASIC COPIES... \n");         
+
+        /*Now, deep copy arrays:*/
+        output->analysis_type_list=xNew<int>(this->nummodels);
+        xMemCpy<int>(output->analysis_type_list,this->analysis_type_list,this->nummodels);
+
+        output->profiler=static_cast<Profiler*>(this->profiler->copy());
+        output->parameters=static_cast<Parameters*>(this->parameters->Copy());
+
+        //_printf_("   PRINTING AFTER PARAMETERS... \n");           
+
+        /*Creating connectivity table*/
+        int* connectivity = xNew<int>(newnumberofvertices);
+        for(int i=0;i<newnumberofvertices;i++) connectivity[i] = 0;
+
+        for (int i=0;i<newnumberofelements;i++){
+                for (int j=0;j<elementswidth;j++){
+                        int vertexid = newelementslist[elementswidth*i+j];
+                        _assert_(vertexid>-1 && vertexid<newnumberofvertices);
+                        connectivity[vertexid]+=1;
+                }
+        }       
+
+        //_printf_("   PRINTING AFTER CONNECTIVITY... \n");         
+
+        /*Creating vertices*/
+        output->vertices = new Vertices();
+
+        for(int i=0;i<newnumberofvertices;i++){
+                if(my_vertices[i]){
+                        Vertex *newvertex = new Vertex();
+                        
+                        newvertex->id           = i+1;
+                        newvertex->sid          = i;
+                        newvertex->pid          = UNDEF;
+
+                        newvertex->x            = newx[i];
+                        newvertex->y            = newy[i];
+                        newvertex->z            = newz[i];
+                        newvertex->domaintype   = Domain2DhorizontalEnum;
+                        newvertex->sigma                = 0.;
+
+                        newvertex->connectivity = connectivity[i];
+
+                        output->vertices->AddObject(newvertex); 
+                } 
+        }
+
+        xDelete<int>(connectivity);
+
+        //_printf_("   PRINTING AFTER VERTICES... \n");     
+
+        /*Creating elements*/
+        /*Just Tria in this version*/
+        output->elements = new Elements();
+
+        for(int i=0;i<newnumberofelements;i++){
+                if(my_elements[i]){
+
+                        Tria *newtria = new Tria();
+
+                        newtria->id  = i+1;
+                        newtria->sid = i;
+                        newtria->parameters = NULL;
+
+                        newtria->inputs  = new Inputs();
+
+                        newtria->nodes    = NULL;
+                        newtria->vertices = NULL;
+                        newtria->material = NULL;
+                        newtria->matpar   = NULL;
+
+                        if(this->nummodels>0){
+                                newtria->element_type_list=xNew<int>(this->nummodels);
+                                for(int j=0;j<nummodels;j++) newtria->element_type_list[j] = 0;
+                        }
+                        else newtria->element_type_list = NULL;
+
+                        /*Element hook*/
+                        int matpar_id           = newnumberofelements+1; //retrieve material parameter id (last pointer in femodel->materials)
+                        int material_id         = i+1; // retrieve material_id = i+1;
+
+                        /*retrieve vertices ids*/
+                        int* vertex_ids = xNew<int>(elementswidth);
+                        
+                        for(int j=0;j<elementswidth;j++){ 
+                                vertex_ids[j]=reCast<int>(newelementslist[elementswidth*i+j]);
+                        }
+
+                        newtria->numanalyses = this->nummodels;
+                        newtria->hnodes      = new Hook*[this->nummodels];
+                        newtria->hvertices   = new Hook(&vertex_ids[0],elementswidth);
+                        newtria->hmaterial   = new Hook(&material_id,1);
+                        newtria->hmatpar     = new Hook(&matpar_id,1);
+                        newtria->hneighbors  = NULL;
+
+                        /*Initialize hnodes as NULL*/
+                        for(int j=0;j<this->nummodels;j++){
+                                newtria->hnodes[j]=NULL;
+                        }
+
+                        /*Clean up*/
+                        xDelete<int>(vertex_ids);
+                        output->elements->AddObject(newtria);   
+                } 
+        }
+
+        xDelete<int>(newelementslist);
+        //_printf_("   PRINTING AFTER ELEMENTS... \n");     
+
+        /*Creating materials*/
+        /*Just Matice in this version*/
+        output->materials = new Materials();
+
+        for(int i=0;i<newnumberofelements;i++){
+                if(my_elements[i]){
+                        output->materials->AddObject(new Matice(i+1,i,MaticeEnum));     
+                } 
+        }
+        /*This is done to follow CreateElementsVerticesAndMaterials, line 57*/
+        //output->elements->InputDuplicate(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum);
+
+        /*Add new constant material property to materials, at the end: */
+        Matpar *newmatpar = static_cast<Matpar*>(this->materials->GetObjectByOffset(this->materials->Size()-1)->copy());
+        output->materials->AddObject(newmatpar);//put it at the end of the materials
+
+        //_printf_("   PRINTING AFTER MATERIALS... \n");            
+
+        delete x;
+        delete y;
+        delete z;
+        for(int i=0;i<numberofelements;i++) delete elementsptr[i];
+    if(elementsptr) delete elementsptr;
+
+        //itapopo to test and print the elements and coordinates
+        //output->InitializeAdaptiveRefinement();
+
+        /*Creating nodes*/
+        /*Just SSA (2D) and P1 in this version*/
+        output->nodes = new Nodes();
+
+        int nodecounter = 0;
+        int lid                 = 0;
+        for(int i=0;i<output->nummodels;i++){
+                
+                int analysis_enum = output->analysis_type_list[i];
+
+                //_printf_("    Analysis type:" << EnumToStringx(analysis_enum) << "\n");
+
+                //itapopo as the domain is 2D, it is not necessary to create nodes for this analysis
+                if(analysis_enum==StressbalanceVerticalAnalysisEnum) continue;      
+
+                for(int j=0;j<newnumberofvertices;j++){
+                        if(my_vertices[j]){
+                                
+                                Node* newnode = new Node();     
+
+                                /*id: */
+                                newnode->id            = nodecounter+j+1;
+                                newnode->sid           = j;
+                                newnode->lid           = lid++;
+                                newnode->analysis_enum = analysis_enum;
+
+                                /*Initialize coord_system: Identity matrix by default*/
+                                for(int k=0;k<3;k++) for(int l=0;l<3;l++) newnode->coord_system[k][l]=0.0;
+                                for(int k=0;k<3;k++) newnode->coord_system[k][k]=1.0;
+
+                                /*indexing:*/
+                                newnode->indexingupdate = true;
+
+                                Analysis* analysis = EnumToAnalysis(analysis_enum);
+                                int *doftypes = NULL;
+                                int numdofs        = analysis->DofsPerNode(&doftypes,Domain2DhorizontalEnum,SSAApproximationEnum);
+                                newnode->indexing.Init(numdofs,doftypes);
+                                xDelete<int>(doftypes);
+                                delete analysis;
+
+                                if(analysis_enum==StressbalanceAnalysisEnum)
+                                 newnode->SetApproximation(SSAApproximationEnum);
+                                else
+                                 newnode->SetApproximation(0);
+
+                                /*Stressbalance Horiz*/
+                                if(analysis_enum==StressbalanceAnalysisEnum){
+                                        // itapopo
+                                        /*Coordinate system provided, convert to coord_system matrix*/
+                                        //XZvectorsToCoordinateSystem(&this->coord_system[0][0],&iomodel->Data(StressbalanceReferentialEnum)[j*6]);
+                                        //_assert_(sqrt( coord_system[0][0]*coord_system[0][0] + coord_system[1][0]*coord_system[1][0]) >1.e-4);
+
+                                }
+
+                                output->nodes->AddObject(newnode);
+                        }
+                }
+
+                if(output->nodes->Size()) nodecounter = output->nodes->MaximumId();
+        }
+        
+        //_printf_("   Old node size: " << this->nodes->Size() << " \n");           
+
+        //for(int i=0;i<this->nummodels;i++){
+        //      int analysis_type = this->analysis_type_list[i];
+        //      _printf_("      Analysis type:" << EnumToStringx(analysis_type) << " size: "<< this->nodes->NumberOfNodes(analysis_type) << " NDofs: " <<  this->nodes->NumberOfDofs(analysis_type, GsetEnum) << " \n");
+        //}
+        
+        //_printf_("   New number of nodes: " << output->nodes->Size() << " \n");           
+        //for(int i=0;i<output->nummodels;i++){
+        //      int analysis_type = output->analysis_type_list[i];
+        //      _printf_("      Analysis type:" << EnumToStringx(analysis_type) << " size: "<< output->nodes->NumberOfNodes(analysis_type) << " NDofs: " <<  this->nodes->NumberOfDofs(analysis_type, GsetEnum) << " \n");
+        //}
+        
+        //_printf_("            Old nodes deep echo: \n");
+        //this->nodes->DeepEcho();
+
+        //_printf_("            New nodes deep echo: \n");
+        //output->nodes->DeepEcho();
+
+        //_printf_("   PRINTING AFTER NODES... \n");        
+
+        /*Create constraints*/
+        IssmDouble *spcvx = NULL;
+        IssmDouble *spcvy = NULL;
+
+        int numberofnodes_analysistype  = this->nodes->NumberOfNodes(StressbalanceAnalysisEnum);
+        
+        //_printf_("   Number of nodes: " << numberofnodes_analysistype << " \n");          
+
+        Vector<IssmDouble>* vspcvx              = new Vector<IssmDouble>(numberofnodes_analysistype);
+        Vector<IssmDouble>* vspcvy              = new Vector<IssmDouble>(numberofnodes_analysistype);
+
+        IssmDouble BigNumber                    = 1.e8;
+
+        for(int i=0;i<numberofnodes_analysistype;i++){
+                vspcvx->SetValue(i,BigNumber,INS_VAL);
+                vspcvy->SetValue(i,BigNumber,INS_VAL);
+        }
+
+        IssmDouble absmaxspcvx = 0;
+        IssmDouble absmaxspcvy = 0;
+
+        for(int i=0;i<this->constraints->Size();i++){
+                SpcStatic* spc          = xDynamicCast<SpcStatic*>(this->constraints->GetObjectByOffset(i));
+                int dof                         = spc->GetDof();
+                int node                        = spc->GetNodeId();
+                IssmDouble spcvalue     = spc->GetValue(); 
+                int nodeindex           = node-1;
+
+                if(dof==0) {
+                        vspcvx->SetValue(nodeindex,spcvalue,INS_VAL);
+                        if(fabs(spcvalue)>absmaxspcvx) absmaxspcvx = fabs(spcvalue);
+                }
+                else {
+                        vspcvy->SetValue(nodeindex,spcvalue,INS_VAL);
+                        if(fabs(spcvalue)>absmaxspcvy) absmaxspcvy = fabs(spcvalue);
+                }
+        }
+
+        /*Assemble*/
+        vspcvx->Assemble();
+        vspcvy->Assemble();
+
+        /*Serialize*/
+        spcvx = vspcvx->ToMPISerial();
+        spcvy = vspcvy->ToMPISerial();
+
+        /*Free the data*/
+        delete vspcvx;
+        delete vspcvy;
+        
+        double *newspcvx        = NULL;
+        double *newspcvy        = NULL;
+        int *oldelements        = newelementslist; //itapopo
+        double *oldx            = x; //itapopo
+        double *oldy            = y; //itapopo
+        int nods_data           = numberofnodes_analysistype;
+        int nels_data           = newnumberofelements;
+        int M_data                      = numberofnodes_analysistype;
+        int N_data              = 1;
+        int N_interp            = newnumberofvertices;//itapopo
+        Options *options        = NULL;
+
+        //itapopo voltar aqui
+        InterpFromMeshToMesh2dx(&newspcvx,oldelements,oldx,oldy,nods_data,nels_data,spcvx,M_data,N_data,newx,newy,N_interp,options);
+        InterpFromMeshToMesh2dx(&newspcvx,oldelements,oldx,oldy,nods_data,nels_data,spcvx,M_data,N_data,newx,newy,N_interp,options);
+
+        output->constraints = new Constraints();
+
+        nodecounter                     = 0; //itapopo deve começar pelo primeiro nó do StressbalanceAnalysis
+        int count                               = 0;
+        int constraintcounter   = 0; //itapopo
+        IssmDouble eps                  = 1.e-2;
+
+        for(int i=0;i<newnumberofvertices;i++){
+                if(my_vertices[i])
+                /*spcvx*/
+                if(fabs(spcvx[i]) < absmaxspcvx+eps){//itapopo
+                        output->constraints->AddObject(new SpcStatic(constraintcounter+count+1,nodecounter+i+1,0,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
+                        count++;
+                }
+        }
+        count=0;
+        for(int i=0;i<newnumberofvertices;i++){
+                if(my_vertices[i])
+                /*spcvy*/
+                if(fabs(spcvy[i]) < absmaxspcvy+eps){//itapopo
+                        output->constraints->AddObject(new SpcStatic(constraintcounter+count+1,nodecounter+i+1,1,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
+                        count++;
+                }
+
+        }
+
+        //_printf_("   Old constraints size: " << this->constraints->Size() << " \n");      
+        //this->constraints->DeepEcho();
+
+        //_printf_("\n");
+
+        //_printf_("   New constraints size: " << output->constraints->Size() << " \n");
+        
+        //_printf_("SPCVX\n");
+        //for(int i = 0;i<newnumberofvertices;i++){
+        //      _printf_("value: " << spcvx[i] << "\n");
+        //}         
+        //_printf_("SPCVY\n");
+        //for(int i = 0;i<newnumberofvertices;i++){
+        //      _printf_("value: " << spcvy[i] << "\n");
+        //}
+
+        //output->constraints->DeepEcho();
+
+        //_printf_("   PRINTING AFTER CONSTRAINTS... \n");          
+
+        //_printf_("   Old loads size: " << this->loads->Size() << " \n");          
+        //this->loads->DeepEcho();
+
+        //_printf_("   PRINTING AFTER LOADS... \n");        
+
+        // output->loads=static_cast<Loads*>(this->loads->Copy());
+        // output->constraints=static_cast<Constraints*>(this->constraints->Copy());
+        // output->results=static_cast<Results*>(this->results->Copy());
+
+        // /*reset hooks for elements, loads and nodes: */
+        // output->elements->ResetHooks();
+        // output->loads->ResetHooks();
+        // output->materials->ResetHooks();
+
+        // /*do the post-processing of the datasets to get an FemModel that can actually run analyses: */
+        // for(i=0;i<nummodels;i++){
+        //      analysis_type=output->analysis_type_list[i];
+        //      output->SetCurrentConfiguration(analysis_type);
+        //      if(i==0) VerticesDofx(output->vertices,output->parameters); //only call once, we only have one set of vertices
+        //      SpcNodesx(output->nodes,output->constraints,output->parameters,analysis_type);
+        //      NodesDofx(output->nodes,output->parameters,analysis_type);
+        //      ConfigureObjectsx(output->elements,output->loads,output->nodes,output->vertices,output->materials,output->parameters);
+        // }
+
+        // /*Reset current configuration: */
+        // analysis_type=output->analysis_type_list[analysis_counter];
+        // output->SetCurrentConfiguration(analysis_type);
+
+        /** TODO
+
+        - generate the required input objects for NeoPZ
+                AMR has fathermesh and previousmesh. On first call, these meshes are generated. 
+
+        - call NeoPZ
+                This creates a newmesh which will be refined. 
+
+        - get the new mesh (index,x,y) from NeoPZ
+                Is is doing by GetNewMesh method.
+
+        - Create a new FemModel* that will be consistent with the new mesh
+                It can be done by a method. (attention with CPU #)
+
+        - Initialize new FemModel based on new mesh (and copy the old FemModel parameters)
+                It can be done by a method. (attention with CPU #)
+
+        - The last and most difficult thing to do is to update the inputs of the elements of the new mesh:
+                It can be done in just one method, which calls GatherInputs, InterpFromMeshToMesh and InputUpdateFromVector
+
+                - CPU #0 will gather all inputs from FemModel
+                - call InterpFromMeshToMesh to interpolate them onto the new Mesh
+                - broadcast the new fields to all cpus
+                - call InputUpdateFromVector so that all inputs are updated
+        
+        - return new FemModel
+        
+        */
+
+        return output;
+}
+/*}}}*/
+#endif

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

@@ -137 +137 @@
                 void InitFromBuffers(char* buffer, int buffersize, char* toolkits, int solution_type,bool trace,IssmPDouble* X=NULL);
                 #endif
+
+                #ifdef _HAVE_AMR_
+                /*Adaptive mesh refinement methods*/
+                void InitializeAdaptiveRefinement(void);
+                FemModel* ReMesh(void);
+                #endif
 };
                 

issm/trunk-jpl/src/c/classes/Materials/Matice.cpp

@@ -37 +37 @@
 Matice::Matice(int matice_mid,int index, IoModel* iomodel){/*{{{*/
 
-        /*Intermediaries:*/
-        int    matice_eid;
-
-        /*Initialize id*/
-        this->mid=matice_mid;
-
-        /*Hooks: */
-        matice_eid=index+1;
-        this->helement=new Hook(&matice_eid,1);
-        this->element=NULL;
-
-         /*Other perporties*/
-   int    materialtype;
+         /*Get material type and initialize object*/
+   int materialtype;
    iomodel->FindConstant(&materialtype,"md.materials.type");
-   if(materialtype==MatdamageiceEnum){
-                this->isdamaged = true;
-                this->isenhanced = false;
-        }
-        else if(materialtype==MaticeEnum){
-                this->isdamaged = false;
-                this->isenhanced = false;
-        }
-        else if(materialtype==MatenhancediceEnum){
-                this->isdamaged = false;
-                this->isenhanced = true;
-        }
-   else _error_("Material type not recognized");
+        this->Init(matice_mid,index,materialtype);
+
+}
+/*}}}*/
+Matice::Matice(int matice_mid,int index,int materialtype){/*{{{*/
+
+        this->Init(matice_mid,index,materialtype);
         return;
-
-}
-/*}}}*/
+} /*}}}*/
 Matice::~Matice(){/*{{{*/
         delete helement;
@@ -73 +54 @@
 }
 /*}}}*/
+void Matice::Init(int matice_mid,int index,int materialtype){/*{{{*/
+
+        /*Initialize id*/
+        this->mid=matice_mid;
+
+        /*Hooks: */
+        int matice_eid=index+1;
+        this->helement=new Hook(&matice_eid,1);
+        this->element=NULL;
+
+        /*Material specific properties*/
+        switch(materialtype){
+                case MatdamageiceEnum:
+                        this->isdamaged = true;
+                        this->isenhanced = false;
+                        break;
+                case MaticeEnum:
+                        this->isdamaged = false;
+                        this->isenhanced = false;
+                        break;
+                case MatenhancediceEnum:
+                        this->isdamaged = false;
+                        this->isenhanced = true;
+                        break;
+                default:
+                        _error_("Material type not recognized");
+        }
+
+        return;
+} /*}}}*/
 
 /*Object virtual functions definitions:*/

issm/trunk-jpl/src/c/classes/Materials/Matice.h

@@ -35 +35 @@
                 Matice();
                 Matice(int mid,int i, IoModel* iomodel);
+                Matice(int mid,int i, int materialtype);
                 ~Matice();
+                void Init(int mid,int i, int materialtype);
                 /*}}}*/
                 /*Object virtual functions definitions:{{{ */

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

@@ -497 +497 @@
 }
 /*}}}*/
+void  Node::SetApproximation(int in_approximation){/*{{{*/
+        
+        this->approximation = in_approximation;
+}
+/*}}}*/
 int  Node::GetNumberOfDofs(int approximation_enum,int setenum){/*{{{*/
 

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

@@ -85 +85 @@
                 void  VecMerge(Vector<IssmDouble>* ug, IssmDouble* vector_serial,int setenum);
                 void  VecReduce(Vector<IssmDouble>* vector, IssmDouble* ug_serial,int setnum);
+                void  SetApproximation(int in_approximation);
 };