Changeset 21504

Timestamp:

01/30/17 16:57:30 (8 years ago)

Author:

tsantos

Message:

CHG: adding methods in FemModel to generate new femmodel based on new mesh (adaptive mesh, not working yet).

Location:

issm/trunk-jpl/src/c/classes

Files:

• FemModel.cpp (modified) (7 diffs)
• FemModel.h (modified) (1 diff)

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

@@ -83 +83 @@
         #ifdef _HAVE_NEOPZ_
         this->InitializeAdaptiveRefinement();
+        FemModel *Test=this->ReMesh();//itapopo: just to test!;
+        printf("   AFTER REMESH!!!\n");
+        //Test->CleanUp();
+        printf("   AFTER CLEANUP!!!\n");
+        //delete Test;
+        printf("   AFTER DELETE!!!\n");
         #endif
 
@@ -2936 +2942 @@
         /*Get initial mesh*/
         /*Get total number of elements and total numbers of elements in the initial mesh*/
-        long numberofvertices, numberofelements, numberofsegments;
+        int numberofvertices, numberofelements, numberofsegments;
         numberofvertices = this->vertices->NumberOfVertices();
         numberofelements = this->elements->NumberOfElements();
@@ -2973 +2979 @@
         IssmDouble *id3 = vid3->ToMPISerial();
         
-        long **elementsptr;
-        elementsptr = new long*[numberofelements];
+        int **elementsptr;
+        elementsptr = new int*[numberofelements];
    for(int i=0;i<numberofelements;i++){
-                elementsptr[i] = new long[elementswidth];
-      elementsptr[i][0] = (long)id1[i];
-                elementsptr[i][1] = (long)id2[i];
-      elementsptr[i][2] = (long)id3[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){ 
-                long nsegments=0;
+                int nsegments=0;
                 this->amr->CreateInitialMesh(numberofvertices, numberofelements, nsegments, elementswidth, x, y, z, elementsptr, NULL);
         }
@@ -3007 +3013 @@
         
         /*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);
-    }
-                
+        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
+                }
+        }       
+        bool* my_elements=NULL; 
+        int* my_vertices=NULL;
+
+        this->ElementsAndVerticesPartitioning(newnumberofvertices,newnumberofelements,elementswidth,newelementslist,&my_elements,&my_vertices);
+        
+        /*Creating new femmodel with new mesh*/
+        FemModel* output=NULL;
+        output=new FemModel(*this);
+
+        /*Copy basic attributes:*/
+        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());      
+
+        /*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*/
+        /*Just Tria in this version*/
+        output->elements=new Elements();
+        this->CreateElements(newnumberofelements,elementswidth,newelementslist,my_elements,output->elements);
+
+        /*Cleanup*/
+        xDelete<int>(newelementslist);      
+
+        /*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++){
+                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*/
-    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);
+        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
+        return output;
+}
+/*}}}*/
+void FemModel::GetMaskLevelSet(IssmDouble **pmasklevelset){/*{{{*/
+
+        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++){
@@ -3083 +3318 @@
                 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->SetValue(elem_vertices[1],elementlevelset[1],INS_VAL);
+      vmasklevelset->SetValue(elem_vertices[2],elementlevelset[2],INS_VAL);
+        }
+
+   /*Assemble*/
         vmasklevelset->Assemble();
-
-        /*Serialize*/
-        masklevelset = vmasklevelset->ToMPISerial();
-
+        
+        /*Serialize and set output*/
+        (*pmasklevelset)=vmasklevelset->ToMPISerial();
+
+        /*Cleanup*/
         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");
-    //}
+}
+/*}}}*/
+void FemModel::CreateVertices(int newnumberofvertices,int newnumberofelements,int elementswidth,int* newelementslist,int* my_vertices,IssmDouble* newx,IssmDouble* newy,IssmDouble* newz,Vertices* vertices){/*{{{*/
+
+        /*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;
+                }
+        }       
+
+        /*Create vertex and insert in 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];
+                        vertices->AddObject(newvertex); 
+                } 
+        }
+
+        xDelete<int>(connectivity);
+}
+/*}}}*/
+void FemModel::CreateElements(int newnumberofelements,int elementswidth,int* newelementslist,bool* my_elements,Elements* 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]);  
+                        /*Setting the hooks*/
+                        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);
+                        elements->AddObject(newtria);   
+                } 
+        }
+}
+/*}}}*/
+void FemModel::CreateMaterials(int newnumberofelements,bool* my_elements,Materials* materials){/*{{{*/
+
+        /*Just Matice in this version*/
+        for(int i=0;i<newnumberofelements;i++){
+                if(my_elements[i]){
+                        materials->AddObject(new Matice(i+1,i,MaticeEnum));     
+                } 
+        }
+        
+        /*Add new constant material property to materials, at the end: */
+        Matpar *newmatpar=static_cast<Matpar*>(this->materials->GetObjectByOffset(this->materials->Size()-1)->copy());
+        materials->AddObject(newmatpar);//put it at the end of the materials        
+
+}
+/*}}}*/
+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;
+        
+        epart=xNew<int>(newnumberofelements);
+        npart=xNew<int>(newnumberofvertices);
 
         /*Partition using Metis:*/
@@ -3169 +3449 @@
         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");       
+                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");         
 
         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;
+        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]){ 
@@ -3197 +3474 @@
         }
 
+        /*Free ressources:*/
         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;
+        xDelete<int>(npart);        
+
+        /*Assign output pointers:*/
+        *pmy_elements=my_elements;
+        *pmy_vertices=my_vertices;
+
 }
 /*}}}*/

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

@@ -149 +149 @@
                 void InitializeAdaptiveRefinement(void);
                 FemModel* ReMesh(void);
+                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);
                 #endif
 };