Changeset 23065

Timestamp:

08/07/18 06:15:18 (7 years ago)

Author:

tsantos

Message:

CHG: only one MPI call after Remesh core

Location:

issm/trunk-jpl

Files:

• externalpackages/neopz/install.sh (modified) (1 diff)
• src/c/classes/AdaptiveMeshRefinement.cpp (modified) (9 diffs)
• src/c/classes/AdaptiveMeshRefinement.h (modified) (2 diffs)
• src/c/classes/AmrBamg.cpp (modified) (3 diffs)
• src/c/classes/AmrBamg.h (modified) (1 diff)
• src/c/classes/FemModel.cpp (modified) (5 diffs)

issm/trunk-jpl/externalpackages/neopz/install.sh

@@ -20 +20 @@
 #Configure neopz using cmake
 cd $PROJECT_SOURCE_DIR
-cmake -DCMAKE_INSTALL_PREFIX:PATH=$PROJECT_BINARY_DIR
+cmake -DCMAKE_INSTALL_PREFIX:PATH=$PROJECT_BINARY_DIR -DCMAKE_BUILD_TYPE=Release -DCMAKE_CXX_FLAGS="-g -O3"
 
 cd $PROJECT_SOURCE_DIR

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

@@ -100 +100 @@
 
 /*Mesh refinement methods*/
-void AdaptiveMeshRefinement::ExecuteRefinement(double* gl_distance,double* if_distance,double* deviatoricerror,double* thicknesserror,int* pnewnumberofvertices,int* pnewnumberofelements,double** px,double** py,int** pelementslist){/*{{{*/
+void AdaptiveMeshRefinement::ExecuteRefinement(double* gl_distance,double* if_distance,double* deviatoricerror,double* thicknesserror,int** pdatalist,double** pxylist,int** pelementslist){/*{{{*/
         
         /*IMPORTANT! pelementslist are in Matlab indexing*/
@@ -123 +123 @@
    
         /*Get new geometric mesh in ISSM data structure*/
-        this->GetMesh(this->previousmesh,pnewnumberofvertices,pnewnumberofelements,px,py,pelementslist);
+        this->GetMesh(this->previousmesh,pdatalist,pxylist,pelementslist);
 
         /*Verify the new geometry*/
-        this->CheckMesh(pnewnumberofvertices,pnewnumberofelements,px,py,pelementslist);
+        this->CheckMesh(pdatalist,pxylist,pelementslist);
 }
 /*}}}*/
@@ -449 +449 @@
 }
 /*}}}*/
-void AdaptiveMeshRefinement::GetMesh(TPZGeoMesh* gmesh,int* nvertices,int* nelements,double** px,double** py, int** pelements){/*{{{*/
+void AdaptiveMeshRefinement::GetMesh(TPZGeoMesh* gmesh,int** pdata,double** pxy, int** pelements){/*{{{*/
 
         /* IMPORTANT! pelements are in Matlab indexing
@@ -462 +462 @@
         int ntotalvertices              = gmesh->NNodes();
         int* newelements                        = NULL;
-        double* newmeshX                        = NULL;
-        double* newmeshY                        = NULL;
+        int* newdata                            = NULL;
+        double* newmeshXY                       = NULL;
         TPZGeoEl* geoel                 = NULL;
         long* vertex_index2sid  = xNew<long>(ntotalvertices);
@@ -494 +494 @@
         }
 
+        /* Create new mesh structure and fill it */
         nconformelements        = (int)this->sid2index.size();
         nconformvertices        = (int)sid;
         newelements                     = xNew<int>(nconformelements*this->GetNumberOfNodes());
-        newmeshX                                = xNew<double>(nconformvertices);
-   newmeshY                             = xNew<double>(nconformvertices);
+        newmeshXY                       = xNew<double>(nconformvertices*2);
+        newdata                         = xNew<int>(2);
+        newdata[0]                      = nconformvertices;
+        newdata[1]                      = nconformelements;
 
         for(int i=0;i<ntotalvertices;i++){//over the TPZNode index (fill in the ISSM vertices coords)
@@ -505 +508 @@
                 TPZVec<REAL> coords(3,0.);
                 gmesh->NodeVec()[i].GetCoordinates(coords);
-                newmeshX[sid] = coords[0];
-                newmeshY[sid] = coords[1];
+                newmeshXY[2*sid]                = coords[0]; // X
+                newmeshXY[2*sid+1]      = coords[1]; // Y
         }
                 
@@ -525 +528 @@
                 c=newelements[i*this->GetNumberOfNodes()+2]-1;
 
-                xa=newmeshX[a]; ya=newmeshY[a];
-                xb=newmeshX[b]; yb=newmeshY[b];
-                xc=newmeshX[c]; yc=newmeshY[c];
+                xa=newmeshXY[2*a]; ya=newmeshXY[2*a+1];
+                xb=newmeshXY[2*b]; yb=newmeshXY[2*b+1];
+                xc=newmeshXY[2*c]; yc=newmeshXY[2*c+1];
 
                 detJ=(xb-xa)*(yc-ya)-(xc-xa)*(yb-ya);
@@ -539 +542 @@
 
         /*Setting outputs*/
-        *nvertices      = nconformvertices;
-        *nelements      = nconformelements;
-        *px                     = newmeshX;
-        *py                = newmeshY;
+        *pdata          = newdata;
+        *pxy               = newmeshXY;
         *pelements      = newelements;
    
@@ -693 +694 @@
 }
 /*}}}*/
-void AdaptiveMeshRefinement::CheckMesh(int* nvertices,int* nelements,double** px,double** py,int** pelements){/*{{{*/
+void AdaptiveMeshRefinement::CheckMesh(int** pdata,double** pxy,int** pelements){/*{{{*/
 
         /*Basic verification*/
-        if(nvertices<=0) _error_("Impossible to continue: nvertices <=0!\n");
-        if(nelements<=0) _error_("Impossible to continue: nelements <=0!\n");
-        if(!px) _error_("Impossible to continue: px is NULL!\n");
-        if(!py) _error_("Impossible to continue: py is NULL!\n");
+        if(!pdata) _error_("Impossible to continue: pdata is NULL!\n");
+        if(pdata[0]<=0) _error_("Impossible to continue: nvertices <=0!\n");
+        if(pdata[1]<=0) _error_("Impossible to continue: nelements <=0!\n");
+        if(!pxy) _error_("Impossible to continue: pxy is NULL!\n");
         if(!pelements) _error_("Impossible to continue: pelements is NULL!\n");
-
-        /*Verify if there are orphan nodes*/
-        std::set<int> elemvertices;
-        elemvertices.clear(); 
-        for(int i=0;i<*nelements;i++){
-                for(int j=0;j<this->GetNumberOfNodes();j++) {
-                        elemvertices.insert((*pelements)[i*this->GetNumberOfNodes()+j]);
-                }
-        }
-        if(elemvertices.size()!=*nvertices) _error_("Impossible to continue: elemvertices.size() != nvertices!\n");
-        
-        //Verify if there are inf or NaN in coords
-        for(int i=0;i<*nvertices;i++){
-                if(std::isnan((*px)[i]) || std::isinf((*px)[i])) _error_("Impossible to continue: px i=" << i <<" is NaN or Inf!\n"); 
-                if(std::isnan((*py)[i]) || std::isinf((*py)[i])) _error_("Impossible to continue: py i=" << i <<" is NaN or Inf!\n");
-        }
-        for(int i=0;i<*nelements;i++){
-                for(int j=0;j<this->GetNumberOfNodes();j++){
-                        if(std::isnan((*pelements)[i*GetNumberOfNodes()+j])) _error_("Impossible to continue: px i=" << i <<" is NaN!\n");
-                        if(std::isinf((*pelements)[i*GetNumberOfNodes()+j])) _error_("Impossible to continue: px i=" << i <<" is Inf!\n");
-                }
-        }
-    
 }
 /*}}}*/

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

@@ -67 +67 @@
         void CleanUp();
         void Initialize();
-        void ExecuteRefinement(int numberofpoints,double* xylist,int* pnewnumberofvertices,int* pnewnumberofelements,double** px,double** py,int** pelementslist);
-        void ExecuteRefinement(int numberofpoints,double* xylist,double* deviatoricerror,double* thicknesserror,int* pnewnumberofvertices,int* pnewnumberofelements,double** px,double** py,int** pelementslist);
-   void ExecuteRefinement(double* gl_distance,double* if_distance,double* deviatoricerror,double* thicknesserror,int* pnewnumberofvertices,int* pnewnumberofelements,double** px,double** py,int** pelementslist);
+   void ExecuteRefinement(double* gl_distance,double* if_distance,double* deviatoricerror,double* thicknesserror,int** pdatalist,double** pxy,int** pelementslist);
         void CreateInitialMesh(int &nvertices,int &nelements,double* x,double* y,int* elements);
-        void CheckMesh(int* nvertices,int* nelements,double** px,double** py,int** pelements);
+        void CheckMesh(int** pdata,double** pxy,int** pelements);
         /*}}}*/
 private:
@@ -86 +84 @@
         void RefineMeshToAvoidHangingNodes(bool &verbose,TPZGeoMesh* gmesh);
         void DeleteSpecialElements(bool &verbose,TPZGeoMesh* gmesh);
-        void GetMesh(TPZGeoMesh* gmesh,int* nvertices,int* nelements,double** px,double** py,int** pelements);
+        void GetMesh(TPZGeoMesh* gmesh,int** pdata,double** pxy,int** pelements);
         TPZGeoMesh* CreateRefPatternMesh(TPZGeoMesh* gmesh);
    inline int GetElemMaterialID(){return 1;} 

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

@@ -98 +98 @@
         delete Th;
 }/*}}}*/
-void AmrBamg::ExecuteRefinementBamg(IssmDouble* field,IssmDouble* hmaxVertices,int* pnewnumberofvertices,int *pnewnumberofelements,IssmDouble** px,IssmDouble** py,IssmDouble** pz,int** pelementslist){/*{{{*/
+void AmrBamg::ExecuteRefinementBamg(IssmDouble* field,IssmDouble* hmaxVertices,int** pdatalist,IssmDouble** pxylist,int** pelementslist){/*{{{*/
 
         /*Intermediaries*/
@@ -150 +150 @@
 
         /*Prepare output*/
-        int nbv = meshout->VerticesSize[0];
-        int nbt = meshout->TrianglesSize[0];
-        IssmDouble *x = xNew<IssmDouble>(nbv);
-        IssmDouble *y = xNew<IssmDouble>(nbv);
-        IssmDouble *z = xNew<IssmDouble>(nbv);
+        int nbv                         = meshout->VerticesSize[0];
+        int nbt                         = meshout->TrianglesSize[0];
+        int *datalist           = xNew<int>(2);
+        IssmDouble *xylist= xNew<IssmDouble>(nbv*2);
+        int* elementslist = xNew<int>(nbt*3);
+        
+        datalist[0] = nbv;
+        datalist[1] = nbt;
+        
         for(int i=0;i<nbv;i++){
-                x[i] = meshout->Vertices[i*3+0];
-                y[i] = meshout->Vertices[i*3+1];
-                z[i] = 0.;
+                xylist[2*i]             = meshout->Vertices[i*3+0];
+                xylist[2*i+1]   = meshout->Vertices[i*3+1];
         }
-        int* elementslist= xNew<int>(nbt*3);
+        
         for(int i=0;i<nbt;i++){
                 elementslist[3*i+0] = reCast<int>(meshout->Triangles[4*i+0]);
@@ -167 +170 @@
         }
 
+        /*Assign pointers*/
+        *pdatalist              = datalist;
+        *pxylist                        = xylist;
+        *pelementslist = elementslist;
+        
         /*Cleanup and return*/
         delete geomout;
-        *pnewnumberofvertices = nbv;
-        *pnewnumberofelements = nbt;
-        *px = x;
-        *py = y;
-        *pz = z;
-        *pelementslist = elementslist;
 }/*}}}*/
 void AmrBamg::SetBamgOpts(IssmDouble hmin_in,IssmDouble hmax_in,IssmDouble err_in,IssmDouble gradation_in){/*{{{*/

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

@@ -35 +35 @@
                 /*General methods*/
                 void Initialize(int* elements,IssmDouble* x,IssmDouble* y,int numberofvertices,int numberofelements);
-                void ExecuteRefinementBamg(IssmDouble* field,IssmDouble* hmaxVertices,int* pnewnumberofvertices,int *pnewnumberofelements,IssmDouble** px,IssmDouble** py,IssmDouble** pz,int** pelementslist);
+                void ExecuteRefinementBamg(IssmDouble* field,IssmDouble* hmaxVertices,int** pdatalist,IssmDouble** pxylist,int** pelementslist);
                 void SetBamgOpts(IssmDouble hmin_in,IssmDouble hmax_in,IssmDouble err_in,IssmDouble gradation_in);
 

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

@@ -4877 +4877 @@
         IssmDouble *newy                        = NULL;
         IssmDouble *newz                        = NULL;
+        IssmDouble *newxylist   = NULL;
         int *newelementslist            = NULL;
+        int* newdatalist        = NULL;
         int newnumberofvertices = -1;
         int newnumberofelements = -1;
@@ -4911 +4913 @@
 
         if(my_rank==0){
-                this->amrbamg->ExecuteRefinementBamg(vector_serial,hmaxvertices_serial,&newnumberofvertices,&newnumberofelements,&newx,&newy,&newz,&newelementslist);
-                if(newnumberofvertices<=0 || newnumberofelements<=0) _error_("Error in the refinement process.");
-        }
-
-        /*Cleanup*/
-        xDelete<IssmDouble>(vector_serial);
-        xDelete<IssmDouble>(hmaxvertices_serial);
-        delete vector;
-
-        /*Send new mesh to others CPU*/
-        ISSM_MPI_Bcast(&newnumberofvertices,1,ISSM_MPI_INT,0,IssmComm::GetComm());
-        ISSM_MPI_Bcast(&newnumberofelements,1,ISSM_MPI_INT,0,IssmComm::GetComm());
-        if(my_rank){
-                newx=xNew<IssmDouble>(newnumberofvertices);
-                newy=xNew<IssmDouble>(newnumberofvertices);
-                newz=xNew<IssmDouble>(newnumberofvertices);
-                newelementslist=xNew<int>(newnumberofelements*this->GetElementsWidth());
-        }
-        ISSM_MPI_Bcast(newx,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-        ISSM_MPI_Bcast(newy,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-        ISSM_MPI_Bcast(newz,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-        ISSM_MPI_Bcast(newelementslist,newnumberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm());
-
-        /*Assign output pointers*/
-        *pnewnumberofvertices = newnumberofvertices;
-        *pnewnumberofelements = newnumberofelements;
-        *pnewx = newx;
-        *pnewy = newy;
-        *pnewz = newz;
-        *pnewelementslist = newelementslist;
+                this->amrbamg->ExecuteRefinementBamg(vector_serial,hmaxvertices_serial,&newdatalist,&newxylist,&newelementslist);
+                if(newdatalist[0]<=0 || newdatalist[1]<=0) _error_("Error in the refinement process.");
+        }
+
+   /*Send new mesh to others CPU's*/
+   if(my_rank) newdatalist=xNew<int>(2);
+   ISSM_MPI_Bcast(newdatalist,2,ISSM_MPI_INT,0,IssmComm::GetComm());
+   newnumberofvertices=newdatalist[0];
+   newnumberofelements=newdatalist[1];
+   if(my_rank){
+      newxylist      =xNew<IssmDouble>(newnumberofvertices*2);
+      newelementslist=xNew<int>(newnumberofelements*this->GetElementsWidth());
+   }
+   ISSM_MPI_Bcast(newxylist,newnumberofvertices*2,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+   ISSM_MPI_Bcast(newelementslist,newnumberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm());
+
+        /*Reorganize the data*/
+   newx=xNew<IssmDouble>(newnumberofvertices);
+   newy=xNew<IssmDouble>(newnumberofvertices);
+   newz=xNewZeroInit<IssmDouble>(newnumberofvertices);
+   for(int i=0;i<newnumberofvertices;i++){
+      newx[i] = newxylist[2*i];
+      newy[i] = newxylist[2*i+1];
+   }
+
+   /*Assign output pointers*/
+   *pnewnumberofvertices = newnumberofvertices;
+   *pnewnumberofelements = newnumberofelements;
+   *pnewx = newx;
+   *pnewy = newy;
+   *pnewz = newz;
+   *pnewelementslist = newelementslist;
+
+   /*Cleanup*/
+   xDelete<int>(newdatalist);
+   xDelete<IssmDouble>(newxylist);
+   xDelete<IssmDouble>(vector_serial);
+   xDelete<IssmDouble>(hmaxvertices_serial);
+   delete vector;
 }
 /*}}}*/
@@ -5190 +5201 @@
         /*pnewelementslist keep vertices in Matlab indexing*/
    int my_rank                                          = IssmComm::GetRank();
-   int numberofelements                 = this->elements->NumberOfElements();
         IssmDouble* gl_distance         = NULL;
         IssmDouble* if_distance         = NULL;
@@ -5198 +5208 @@
    IssmDouble* newy                             = NULL;
    IssmDouble* newz                             = NULL;
+        IssmDouble* newxylist           = NULL;
    int* newelementslist                 = NULL;
-   int newnumberofvertices              = -1;
+        int* newdatalist                                = NULL;
+        int newnumberofvertices         = -1;
         int newnumberofelements         = -1;
 
@@ -5210 +5222 @@
         if(my_rank==0){
                 this->amr->ExecuteRefinement(gl_distance,if_distance,deviatoricerror,thicknesserror,
-                                                                                                &newnumberofvertices,&newnumberofelements,&newx,&newy,&newelementslist);
-                newz=xNewZeroInit<IssmDouble>(newnumberofvertices);
-                if(newnumberofvertices<=0 || newnumberofelements<=0) _error_("Error in the ReMeshNeopz.");
-        }
-
-        /*Send new mesh to others CPU*/
-        ISSM_MPI_Bcast(&newnumberofvertices,1,ISSM_MPI_INT,0,IssmComm::GetComm());
-        ISSM_MPI_Bcast(&newnumberofelements,1,ISSM_MPI_INT,0,IssmComm::GetComm());
-        if(my_rank){
-                newx=xNew<IssmDouble>(newnumberofvertices);
-                newy=xNew<IssmDouble>(newnumberofvertices);
-                newz=xNew<IssmDouble>(newnumberofvertices);
-                newelementslist=xNew<int>(newnumberofelements*this->GetElementsWidth());
-        }
-        ISSM_MPI_Bcast(newx,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-        ISSM_MPI_Bcast(newy,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-        ISSM_MPI_Bcast(newz,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-        ISSM_MPI_Bcast(newelementslist,newnumberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm());
-
+                                                                                                &newdatalist,&newxylist,&newelementslist);
+                if(newdatalist[0]<=0 || newdatalist[1]<=0) _error_("Error in the ReMeshNeopz.");        
+        }
+
+        /*Send new mesh to others CPU's*/
+        if(my_rank) newdatalist=xNew<int>(2);
+   ISSM_MPI_Bcast(newdatalist,2,ISSM_MPI_INT,0,IssmComm::GetComm());
+   newnumberofvertices=newdatalist[0];
+   newnumberofelements=newdatalist[1];
+   if(my_rank){
+      newxylist      =xNew<IssmDouble>(newnumberofvertices*2);
+      newelementslist=xNew<int>(newnumberofelements*this->GetElementsWidth());
+   }
+   ISSM_MPI_Bcast(newxylist,newnumberofvertices*2,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+   ISSM_MPI_Bcast(newelementslist,newnumberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm());
+
+   /*Reorganize the data*/
+   newx=xNew<IssmDouble>(newnumberofvertices);
+   newy=xNew<IssmDouble>(newnumberofvertices);
+   newz=xNewZeroInit<IssmDouble>(newnumberofvertices);
+   for(int i=0;i<newnumberofvertices;i++){
+      newx[i] = newxylist[2*i];
+      newy[i] = newxylist[2*i+1];
+   }
+        
         /*Assign the pointers*/
-        (*pnewelementslist)     = newelementslist; //Matlab indexing
-        (*pnewx)                                        = newx;
-        (*pnewy)                                        = newy;
-        (*pnewz)                                        = newz;
-        *pnewnumberofvertices= newnumberofvertices;
-        *pnewnumberofelements= newnumberofelements;
+   (*pnewelementslist)  = newelementslist; //Matlab indexing
+   (*pnewx)             = newx;
+   (*pnewy)             = newy;
+   (*pnewz)             = newz;
+   *pnewnumberofvertices= newnumberofvertices;
+   *pnewnumberofelements= newnumberofelements;
 
         /*Cleanup*/
-        xDelete<IssmDouble>(deviatoricerror);
-        xDelete<IssmDouble>(thicknesserror);
-        xDelete<IssmDouble>(gl_distance);
-        xDelete<IssmDouble>(if_distance);
+   xDelete<int>(newdatalist);
+   xDelete<IssmDouble>(newxylist);
+   xDelete<IssmDouble>(deviatoricerror);
+   xDelete<IssmDouble>(thicknesserror);
+   xDelete<IssmDouble>(gl_distance);
+   xDelete<IssmDouble>(if_distance);
 }
 /*}}}*/