Changeset 21516

Timestamp:

02/03/17 15:42:25 (8 years ago)

Author:

tsantos

Message:

CHG: in FemModel.h and FemModel.cpp, changes and improvements in some AMR methods (adaptative mesh refinement). NOTE: not working yet. In Makefile.am, InterpFromMeshToMesh2dx.cpp was moved to issm_cores.

Location:

issm/trunk-jpl/src/c

Files:

• Makefile.am (modified) (2 diffs)
• classes/FemModel.cpp (modified) (12 diffs)
• classes/FemModel.h (modified) (1 diff)

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

@@ -264 +264 @@
                                         ./modules/ConstraintsStatex/RiftConstraintsState.cpp\
                                         ./modules/ModelProcessorx/CreateOutputDefinitions.cpp\
-                                        ./modules/OutputDefinitionsResponsex/OutputDefinitionsResponsex.cpp\
+                                        ./modules/OutputDefinitionsResponsex/OutputDefinitionsResponsex.cpp\    
+                                        ./modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp\
                                         ./classes/Inputs/PentaInput.cpp\
                                         ./classes/Inputs/TetraInput.cpp
@@ -549 +550 @@
                         ./modules/InterpFromMesh2dx/InterpFromMesh2dx.cpp\
                         ./modules/InterpFromMesh2dx/InterpFromMesh2dxt.cpp\
-                        ./modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp\
                         ./modules/InterpFromMeshToMesh3dx/InterpFromMeshToMesh3dx.cpp\
                         ./modules/InterpFromMeshToGridx/InterpFromMeshToGridx.cpp\

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

@@ -84 +84 @@
         this->InitializeAdaptiveRefinement();
         FemModel *Test=this->ReMesh();//itapopo: just to test!;
-        printf("   AFTER REMESH!!!\n");
+        //Test->elements->DeepEcho();
         //Test->CleanUp();
-        printf("   AFTER CLEANUP!!!\n");
         //delete Test;
-        printf("   AFTER DELETE!!!\n");
         #endif
 
@@ -3012 +3010 @@
 FemModel* FemModel::ReMesh(void){/*{{{*/
         
-        /*All indexing here is in C type: 0..n-1*/
-        int my_rank=IssmComm::GetRank();
-   const int elementswidth=3;//just 2D mesh, tria elements
-        
-        /*Solutions which will be used to refine the elements*/
-        IssmDouble *vx=NULL; //This will be used in constraints
-        IssmDouble *vy=NULL; //This will be used in constraints
-        IssmDouble *masklevelset=NULL;
- 
-        this->GetMaskLevelSet(&masklevelset);
-        // #1: TEST MASK LEVEL SET
-        _printf_("   PRINTING MASKLEVELSET... \n");
-        int numberofvertices=this->vertices->NumberOfVertices();            
-        if(my_rank==0) for(int i=0;i<numberofvertices;i++) _printf_("vertex: " << i << "\t" << masklevelset[i] << "\n");
-
-        /*Refine the mesh and get the new mesh*/
-        IssmDouble *newx;
-        IssmDouble *newy;
-        IssmDouble *newz;
-        int **newelements;
-        int **newsegments;
-        int newnumberofvertices, newnumberofelements, newnumberofsegments;
-        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.");
-        
-        /*Cleanup masklevetset*/
-        xDelete<IssmDouble>(masklevelset);
-
-        // #2: TEST NEW MESH
-        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. Maybe ElementsAndVerticesPartitioning.cpp could be modified to set this without iomodel.*/
-        /*Fill the element list to partitioning*/       
-        int* newelementslist=NULL;
-        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
-                }
-        }       
+        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
+
+        /*Execute refinement and get the new mesh*/
+        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=NULL;
-        output=new FemModel(*this);
+        FemModel* output=new FemModel(*this);
 
         /*Copy basic attributes:*/
-        output->nummodels=this->nummodels;
-        output->solution_type=this->solution_type;
-        output->analysis_counter=this->analysis_counter;
+        output->nummodels        = this->nummodels;
+        output->solution_type    = this->solution_type;
+        output->analysis_counter = this->analysis_counter;
 
         /*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());      
-
+        if(this->loads->Size()!=0){
+                _error_("not supported yet");
+        }
+        else{
+                output->loads=new Loads();
+        }
+        
+        //itapopo
+        // output->results=static_cast<Results*>(this->results->Copy());
+        
         /*Create vertices*/
         output->vertices=new Vertices();
         this->CreateVertices(newnumberofvertices,newnumberofelements,elementswidth,newelementslist,my_vertices,newx,newy,newz,output->vertices);
- 
-        /* #4: TEST THE VERTICES!*/
-        //_printf_("     Old vertices deep echo: \n");
-   //this->vertices->DeepEcho();
-   //_printf_("     New vertices deep echo: \n");
-   //output->vertices->DeepEcho();
 
         /*Creating elements*/
@@ -3094 +3061 @@
         output->elements=new Elements();
         this->CreateElements(newnumberofelements,elementswidth,newelementslist,my_elements,output->elements);
-
         /*Cleanup*/
-        xDelete<int>(newelementslist);      
+//      xDelete<int>(newelementslist); itapopo      
 
         /*Creating materials*/
         output->materials=new Materials();
         this->CreateMaterials(newnumberofelements,my_elements,output->materials);
+        
         /*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++){//create nodes for each analysis in analysis_type_list
+        
+                int analysis_enum = output->analysis_type_list[i];
+                
+                /*As the domain is 2D, it is not necessary to create nodes for this analysis*/
+                /*itapopo must verify if domain is not 3D. Only 2D in this version!*/
+                if(analysis_enum==StressbalanceVerticalAnalysisEnum) continue;      
+                
+                this->CreateNodes(newnumberofvertices,my_vertices,nodecounter,analysis_enum,output->nodes);
+                this->UpdateElements(newnumberofelements,newelementslist,my_elements,nodecounter,i,output->elements);
+
+                if(output->nodes->Size()) nodecounter = output->nodes->MaximumId();
+        }
+
+        /*Create constraints*/
+        output->constraints=new Constraints();
+        this->CreateConstraints(newnumberofvertices,newnumberofelements,newelementslist,newx,newy,my_vertices,output->constraints);
+        
+        output->elements->Presort();
+        output->nodes->Presort();
+        output->vertices->Presort();
+        output->loads->Presort();
+        output->materials->Presort();
+        output->constraints->Presort();
+
+        /*reset hooks for elements, loads and nodes: */
+        output->elements->ResetHooks();
+        output->loads->ResetHooks();
+        output->materials->ResetHooks();
+
+        printf("-------------- file: FemModel.cpp line: %i\n",__LINE__); 
+        /*Finally: interpolate all inputs*/
+//      this->InterpolateInputs(output);
+
+        printf("-------------- file: FemModel.cpp line: %i\n",__LINE__); 
+        /*do the post-processing of the datasets to get an FemModel that can actually run analyses: */
+        int analysis_type;
         for(int i=0;i<output->nummodels;i++){
-                int analysis_enum = output->analysis_type_list[i];
-                //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();
-        }
-        
-        /* #4: TEST THE NODES!*/
-        _printf_("              Old nodes deep echo: \n");
-        this->nodes->DeepEcho();
-        _printf_("              New nodes deep echo: \n");
-        output->nodes->DeepEcho();
-
-        printf(" I arrived here!!!!!!!\n");
-#ifdef _CONTINUE_NEOPZ_
-
-        /*Create constraints*/
-        IssmDouble *spcvx=NULL;
-        IssmDouble *spcvy=NULL;
-        int numberofnodes_analysistype=this->nodes->NumberOfNodes(StressbalanceAnalysisEnum);
-        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;
-        
-        IssmDouble *newspcvx=NULL;
-        IssmDouble *newspcvy=NULL;
-        int *oldelements=newelementslist; //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;
-
-        InterpFromMeshToMesh2dx(&newspcvx,oldelements,vx,vy,nods_data,nels_data,spcvx,M_data,N_data,newx,newy,N_interp,options);
-        InterpFromMeshToMesh2dx(&newspcvx,oldelements,vx,vy,nods_data,nels_data,spcvx,M_data,N_data,newx,newy,N_interp,options);
-// TEST: New SPCVX and VY must be tested!
-
-        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++;
-                }
-
-        }
-
-        /* #5: TEST CONSTRAINTS*/           
-        this->constraints->DeepEcho();
-        output->constraints->DeepEcho();
-
-        // 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
-        
-        */
-#endif
+                analysis_type=output->analysis_type_list[i];
+                output->SetCurrentConfiguration(analysis_type);
+                ConfigureObjectsx(output->elements,output->loads,output->nodes,output->vertices,output->materials,output->parameters);
+                if(i==0){ 
+                        VerticesDofx(output->vertices,output->parameters); //only call once, we only have one set of vertices
+                        //GetMaskOfIceVerticesLSMx(output); //itapopo it needs element->inputs
+                }
+                SpcNodesx(output->nodes,output->constraints,output->parameters,analysis_type);
+                NodesDofx(output->nodes,output->parameters,analysis_type);
+        }
+
+        /*Reset current configuration: */
+        analysis_type=output->analysis_type_list[this->analysis_counter];
+        output->SetCurrentConfiguration(analysis_type);
+
+
         return output;
+}
+/*}}}*/
+void FemModel::InterpolateInputs(FemModel* newfemmodel){/*{{{*/
+
+        int maxinputs = MaximumNumberOfDefinitionsEnum;
+
+        /*Figure out how many inputs we have and their respective interpolation*/
+        Vector<IssmDouble>* input_interpolations=new Vector<IssmDouble>(maxinputs);
+        if(this->elements->Size()){
+                Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(0));
+                element->GetInputsInterpolations(input_interpolations);
+        }
+        input_interpolations->Assemble();
+        
+        /*Serialize and set output*/
+        IssmDouble* input_interpolations_serial = input_interpolations->ToMPISerial();
+        delete input_interpolations;
+
+        /*Count and get enums of all inputs in old mesh*/
+        int  numP0inputs;
+        int  numP1inputs;
+        int *P0input_enums  = NULL;
+        int *P1input_enums  = NULL;
+        int *P0input_interp = NULL;
+        int *P1input_interp = NULL;
+        for(int step=0;step<2;step++){
+                if(step){
+                        P0input_enums = xNew<int>(numP0inputs);
+                        P1input_enums = xNew<int>(numP1inputs);
+                }
+                numP0inputs = 0;
+                numP1inputs = 0;
+                for(int i=0;i<maxinputs;i++){
+                        int inputinterp = reCast<int>(input_interpolations_serial[i]);
+                        switch(inputinterp){
+                                case 0:
+                                        /*Input not found, go to the next*/
+                                        break;
+                                case P1Enum:
+                                        if(step){
+                                                P1input_enums[numP1inputs]  = i;
+                                                P1input_interp[numP1inputs] = inputinterp;
+                                        }
+                                         numP1inputs++;
+                                        break;
+                                case P0Enum:
+                                case DoubleInputEnum:
+                                case IntInputEnum:
+                                case BoolInputEnum:
+                                        if(step){
+                                                P0input_enums[numP0inputs]  = i;
+                                                P0input_interp[numP0inputs] = inputinterp;
+                                        }
+                                        numP0inputs++;
+                                        break;
+                                default:
+                                        _error_(EnumToStringx(inputinterp)<<" Not supported yet");
+                        }
+                }
+        }
+
+        printf("Found %i %i inputs\n",numP0inputs,numP1inputs);
+
+        /*========== Deal with P0 inputs ==========*/
+        int numelementsold = this->elements->NumberOfElements();
+        int numelementsnew = newfemmodel->elements->NumberOfElements();
+        int numverticesold = this->vertices->NumberOfVertices();
+        int numverticesnew = newfemmodel->vertices->NumberOfVertices();
+        IssmDouble* P0inputsold = xNew<IssmDouble>(numelementsold*numP0inputs);
+        IssmDouble* P0inputsnew = NULL;
+
+        for(int i=0;i<numP0inputs;i++){
+                printf("Dealing with %s (type: %s)\n",EnumToStringx(P0input_enums[i]),EnumToStringx(P0input_interp[i]));
+        }
+
+        _error_("stop");
+        /*Old mesh coordinates*/
+        IssmDouble *Xold     = NULL;
+        IssmDouble *Yold     = NULL;
+        int        *Indexold = NULL;
+        IssmDouble *Xnew     = NULL;
+        IssmDouble *Ynew     = NULL;
+
+
+        InterpFromMeshToMesh2dx(&P0inputsnew,Indexold,Xold,Yold,numverticesold,numelementsold,
+                                P0inputsold,numelementsold,numP0inputs,
+                                Xnew,Ynew,numelementsnew,NULL);
+
+        _error_("STOP");
+
+        xDelete<IssmDouble>(P0inputsold);
+        xDelete<IssmDouble>(P0inputsnew);
+
+
+        _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
+        
+        /*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*/
+        xDelete<IssmDouble>(masklevelset);
 }
 /*}}}*/
@@ -3338 +3305 @@
 
         /*Creating connectivity table*/
-        int* connectivity=xNew<int>(newnumberofvertices);
-        for(int i=0;i<newnumberofvertices;i++) connectivity[i]=0;
+        int* connectivity=NULL;
+        connectivity=xNewZeroInit<int>(newnumberofvertices);
 
         for (int i=0;i<newnumberofelements;i++){
                 for (int j=0;j<elementswidth;j++){
-                        int vertexid = newelementslist[elementswidth*i+j];
+                        int vertexid = newelementslist[elementswidth*i+j];//C indexing
                         _assert_(vertexid>-1 && vertexid<newnumberofvertices);
                         connectivity[vertexid]+=1;
@@ -3362 +3329 @@
                         newvertex->sigma=0.;
                         newvertex->connectivity=connectivity[i];
+                        newvertex->clone=false;//itapopo check this
                         vertices->AddObject(newvertex); 
                 } 
@@ -3373 +3341 @@
         for(int i=0;i<newnumberofelements;i++){
                 if(my_elements[i]){
+                        /*Create element - just tria in this version*/
                         Tria *newtria=new Tria();
                         newtria->id=i+1;
@@ -3387 +3356 @@
                         }
                         else newtria->element_type_list=NULL;
+                        
                         /*Element hook*/
                         int matpar_id=newnumberofelements+1; //retrieve material parameter id (last pointer in femodel->materials)
@@ -3392 +3362 @@
                         /*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]);  
+                        for(int j=0;j<elementswidth;j++)        vertex_ids[j]=reCast<int>(newelementslist[elementswidth*i+j])+1;//this Hook wants Matlab indexing       
                         /*Setting the hooks*/
                         newtria->numanalyses =this->nummodels;
@@ -3407 +3377 @@
                 } 
         }
+
+        //itapopo there is this line in CreateElementsVerticesAndMaterials.
+        //elements->InputDuplicate(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum);
 }
 /*}}}*/
@@ -3424 +3397 @@
 }
 /*}}}*/
-void FemModel::ElementsAndVerticesPartitioning(int newnumberofvertices,int newnumberofelements,int elementswidth,int* newelementslist,bool** pmy_elements,int** pmy_vertices){/*{{{*/
+void FemModel::CreateNodes(int newnumberofvertices,int* my_vertices,int nodecounter,int analysis_enum,Nodes* nodes){/*{{{*/
+
+        int lid=0;
+
+        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 this code is rarely used. 
+                                /*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);
+
+                        }
+                        nodes->AddObject(newnode);
+                }
+        }
+        return;
+}
+/*}}}*/
+void FemModel::CreateConstraints(int newnumberofvertices,int newnumberofelements,int* newelementslist,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
+        }       
+
+        /*itapopo ATTENTION: JUST SPCVX AND SPCVY TO TEST!!!*/
+        /*OTHERS CONSTRAINTS MUST BE IMPLEMENTED!!!*/
+
+        /*Get spcvx and spcvy for mesh i-1*/
+        IssmDouble *spcvx=NULL;
+        IssmDouble *spcvy=NULL;
+        IssmDouble *spcvxflag=NULL;
+        IssmDouble *spcvyflag=NULL;
+        int numberofnodes_analysistype=this->nodes->NumberOfNodes(StressbalanceAnalysisEnum);
+        Vector<IssmDouble>* vspcvx=new Vector<IssmDouble>(numberofnodes_analysistype);
+        Vector<IssmDouble>* vspcvy=new Vector<IssmDouble>(numberofnodes_analysistype);
+        Vector<IssmDouble>* vspcvxflag=new Vector<IssmDouble>(numberofnodes_analysistype);
+        Vector<IssmDouble>* vspcvyflag=new Vector<IssmDouble>(numberofnodes_analysistype);
+        for(int i=0;i<numberofnodes_analysistype;i++){ 
+                vspcvx->SetValue(i,0.,INS_VAL);
+                vspcvy->SetValue(i,0.,INS_VAL);
+                vspcvxflag->SetValue(i,0.,INS_VAL);
+                vspcvyflag->SetValue(i,0.,INS_VAL);
+        }
+
+        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) {//vx
+                        vspcvx->SetValue(nodeindex,spcvalue,INS_VAL);
+                        vspcvxflag->SetValue(nodeindex,1,INS_VAL);
+                } 
+                else if(dof==1){//vy
+                        vspcvy->SetValue(nodeindex,spcvalue,INS_VAL);
+                        vspcvyflag->SetValue(nodeindex,1,INS_VAL);
+                }
+                else{
+                        /*nothing here*/
+                }
+        }
+
+        /*Assemble*/
+        vspcvx->Assemble();
+        vspcvy->Assemble();
+        vspcvxflag->Assemble();
+        vspcvyflag->Assemble();
+
+        /*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();
+
+  /*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);
+        
+        int nodecounter                 = 0; //itapopo deve começar pelo primeiro nó do StressbalanceAnalysis
+        int count                                       = 0;
+        int constraintcounter   = 0; //itapopo
+        IssmDouble eps                          = 1.e-8;
+
+        /*Now, insert the interpolated constraints in the data set (constraints)*/
+        for(int i=0;i<newnumberofvertices;i++){
+                if(my_vertices[i])
+                /*spcvx*/
+                if(!xIsNan<IssmDouble>(newspcvx[i]) && newspcvxflag[i]>(1-eps)){
+                        constraints->AddObject(new SpcStatic(constraintcounter+count+1,nodecounter+i+1,0,newspcvx[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(!xIsNan<IssmDouble>(newspcvy[i]) && newspcvyflag[i]>(1-eps) ){
+                        constraints->AddObject(new SpcStatic(constraintcounter+count+1,nodecounter+i+1,1,newspcvy[i],StressbalanceAnalysisEnum)); 
+                        //add count'th spc, on node i+1, setting dof 1 to vx.
+                        count++;
+                }
+
+        }
+}
+/*}}}*/
+void FemModel::UpdateElements(int newnumberofelements,int* newelementslist,bool* my_elements,int nodecounter,int analysis_counter,Elements* newelements){/*{{{*/
+
+        /*Update elements, set hnode.
+        This code is in all analysis */
+        int elemcounter=0;
+        for(int iel=0;iel<newnumberofelements;iel++){
+                if(my_elements[iel]){
+                        Tria* tria=(Tria*)newelements->GetObjectByOffset(elemcounter);
+                        //element update
+                        tria->element_type_list[analysis_counter]=P1Enum;
+                        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->SetHookNodes(tria_node_ids,numnodes,analysis_counter); tria->nodes=NULL;
+                xDelete<int>(tria_node_ids);
+                        elemcounter++;
+                }
+        }
+        return;
+}
+/*}}}*/
+void FemModel::ElementsAndVerticesPartitioning(int& newnumberofvertices,int& newnumberofelements,int& elementswidth,int* newelementslist,bool** pmy_elements,int** pmy_vertices){/*{{{*/
 
         int *epart=NULL; //element partitioning.
@@ -3436 +3628 @@
         int *my_vertices=NULL;
         
+        _assert_(newnumberofvertices>0); 
+        _assert_(newnumberofelements>0); 
         epart=xNew<int>(newnumberofelements);
         npart=xNew<int>(newnumberofvertices);
@@ -3469 +3663 @@
                          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;
                         }

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

@@ -149 +149 @@
                 void InitializeAdaptiveRefinement(void);
                 FemModel* ReMesh(void);
+                void ExecuteRefinement(int &newnumberofvertices,int &newnumberofelements,IssmDouble** newx,IssmDouble** newy,IssmDouble** newz,int** newelementslist);
                 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);
                 void CreateElements(int newnumberofelements,int elementswidth,int* newelementslist,bool* my_elements,Elements* elements);
                 void CreateMaterials(int newnumberofelements,bool* my_elements,Materials* materials);
-                void ElementsAndVerticesPartitioning(int newnumberofvertices,int newnumberofelements,int elementswidth,int* newelementslist,bool** pmy_elements,int** pmy_vertices);
+                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 InterpolateInputs(FemModel* femmodel);
+                void UpdateElements(int newnumberofelements,int* newelementslist,bool* my_elements,int nodecounter,int analysis_counter,Elements* newelements);
+                void ElementsAndVerticesPartitioning(int& newnumberofvertices,int& newnumberofelements,int& elementswidth,int* newelementslist,bool** pmy_elements,int** pmy_vertices);
                 #endif
 };