Index: /issm/trunk-jpl/externalpackages/neopz/install.sh =================================================================== --- /issm/trunk-jpl/externalpackages/neopz/install.sh (revision 23064) +++ /issm/trunk-jpl/externalpackages/neopz/install.sh (revision 23065) @@ -20,5 +20,5 @@ #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 Index: /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp (revision 23064) +++ /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp (revision 23065) @@ -100,5 +100,5 @@ /*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,8 +123,8 @@ /*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,5 +449,5 @@ } /*}}}*/ -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,6 +462,6 @@ 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(ntotalvertices); @@ -494,9 +494,12 @@ } + /* Create new mesh structure and fill it */ nconformelements = (int)this->sid2index.size(); nconformvertices = (int)sid; newelements = xNew(nconformelements*this->GetNumberOfNodes()); - newmeshX = xNew(nconformvertices); - newmeshY = xNew(nconformvertices); + newmeshXY = xNew(nconformvertices*2); + newdata = xNew(2); + newdata[0] = nconformvertices; + newdata[1] = nconformelements; for(int i=0;i 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,7 +528,7 @@ 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,8 +542,6 @@ /*Setting outputs*/ - *nvertices = nconformvertices; - *nelements = nconformelements; - *px = newmeshX; - *py = newmeshY; + *pdata = newdata; + *pxy = newmeshXY; *pelements = newelements; @@ -693,35 +694,12 @@ } /*}}}*/ -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 elemvertices; - elemvertices.clear(); - for(int i=0;i<*nelements;i++){ - for(int j=0;jGetNumberOfNodes();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;jGetNumberOfNodes();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"); - } - } - } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h =================================================================== --- /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h (revision 23064) +++ /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h (revision 23065) @@ -67,9 +67,7 @@ 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,5 +84,5 @@ 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;} Index: /issm/trunk-jpl/src/c/classes/AmrBamg.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/AmrBamg.cpp (revision 23064) +++ /issm/trunk-jpl/src/c/classes/AmrBamg.cpp (revision 23065) @@ -98,5 +98,5 @@ 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,15 +150,18 @@ /*Prepare output*/ - int nbv = meshout->VerticesSize[0]; - int nbt = meshout->TrianglesSize[0]; - IssmDouble *x = xNew(nbv); - IssmDouble *y = xNew(nbv); - IssmDouble *z = xNew(nbv); + int nbv = meshout->VerticesSize[0]; + int nbt = meshout->TrianglesSize[0]; + int *datalist = xNew(2); + IssmDouble *xylist= xNew(nbv*2); + int* elementslist = xNew(nbt*3); + + datalist[0] = nbv; + datalist[1] = nbt; + for(int i=0;iVertices[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(nbt*3); + for(int i=0;i(meshout->Triangles[4*i+0]); @@ -167,12 +170,11 @@ } + /*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){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/AmrBamg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/AmrBamg.h (revision 23064) +++ /issm/trunk-jpl/src/c/classes/AmrBamg.h (revision 23065) @@ -35,5 +35,5 @@ /*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); Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 23064) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 23065) @@ -4877,5 +4877,7 @@ IssmDouble *newy = NULL; IssmDouble *newz = NULL; + IssmDouble *newxylist = NULL; int *newelementslist = NULL; + int* newdatalist = NULL; int newnumberofvertices = -1; int newnumberofelements = -1; @@ -4911,34 +4913,43 @@ 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(vector_serial); - xDelete(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(newnumberofvertices); - newy=xNew(newnumberofvertices); - newz=xNew(newnumberofvertices); - newelementslist=xNew(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(2); + ISSM_MPI_Bcast(newdatalist,2,ISSM_MPI_INT,0,IssmComm::GetComm()); + newnumberofvertices=newdatalist[0]; + newnumberofelements=newdatalist[1]; + if(my_rank){ + newxylist =xNew(newnumberofvertices*2); + newelementslist=xNew(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(newnumberofvertices); + newy=xNew(newnumberofvertices); + newz=xNewZeroInit(newnumberofvertices); + for(int i=0;i(newdatalist); + xDelete(newxylist); + xDelete(vector_serial); + xDelete(hmaxvertices_serial); + delete vector; } /*}}}*/ @@ -5190,5 +5201,4 @@ /*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,6 +5208,8 @@ 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,36 +5222,44 @@ if(my_rank==0){ this->amr->ExecuteRefinement(gl_distance,if_distance,deviatoricerror,thicknesserror, - &newnumberofvertices,&newnumberofelements,&newx,&newy,&newelementslist); - newz=xNewZeroInit(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(newnumberofvertices); - newy=xNew(newnumberofvertices); - newz=xNew(newnumberofvertices); - newelementslist=xNew(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(2); + ISSM_MPI_Bcast(newdatalist,2,ISSM_MPI_INT,0,IssmComm::GetComm()); + newnumberofvertices=newdatalist[0]; + newnumberofelements=newdatalist[1]; + if(my_rank){ + newxylist =xNew(newnumberofvertices*2); + newelementslist=xNew(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(newnumberofvertices); + newy=xNew(newnumberofvertices); + newz=xNewZeroInit(newnumberofvertices); + for(int i=0;i(deviatoricerror); - xDelete(thicknesserror); - xDelete(gl_distance); - xDelete(if_distance); + xDelete(newdatalist); + xDelete(newxylist); + xDelete(deviatoricerror); + xDelete(thicknesserror); + xDelete(gl_distance); + xDelete(if_distance); } /*}}}*/