Index: ../trunk-jpl/externalpackages/neopz/install.sh =================================================================== --- ../trunk-jpl/externalpackages/neopz/install.sh (revision 23064) +++ ../trunk-jpl/externalpackages/neopz/install.sh (revision 23065) @@ -19,7 +19,7 @@ #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 make Index: ../trunk-jpl/src/c/classes/AmrBamg.h =================================================================== --- ../trunk-jpl/src/c/classes/AmrBamg.h (revision 23064) +++ ../trunk-jpl/src/c/classes/AmrBamg.h (revision 23065) @@ -34,7 +34,7 @@ /*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); private: Index: ../trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h =================================================================== --- ../trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h (revision 23064) +++ ../trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h (revision 23065) @@ -66,11 +66,9 @@ /*General methods*/ 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: /*Private attributes{{{*/ @@ -85,7 +83,7 @@ void RefineMeshWithSmoothing(bool &verbose,TPZGeoMesh* gmesh); 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;} inline int GetNumberOfNodes(){return 3;} Index: ../trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp =================================================================== --- ../trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp (revision 23064) +++ ../trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp (revision 23065) @@ -99,7 +99,7 @@ /*}}}*/ /*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*/ /*NEOPZ works only in C indexing*/ @@ -122,10 +122,10 @@ this->RefineMeshOneLevel(verbose,gl_distance,if_distance,deviatoricerror,thicknesserror); /*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); } /*}}}*/ void AdaptiveMeshRefinement::RefineMeshOneLevel(bool &verbose,double* gl_distance,double* if_distance,double* deviatoricerror,double* thicknesserror){/*{{{*/ @@ -448,7 +448,7 @@ gmesh->BuildConnectivity(); } /*}}}*/ -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 NEOPZ works only in C indexing. @@ -461,8 +461,8 @@ int nconformelements,nconformvertices; 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); this->index2sid.clear(); this->index2sid.resize(gmesh->NElements()); @@ -493,11 +493,14 @@ } } + /* 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 } for(int i=0;isid2index.size();i++){//over the sid (fill the ISSM elements) @@ -524,9 +527,9 @@ b=newelements[i*this->GetNumberOfNodes()+1]-1; 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); @@ -538,10 +541,8 @@ } /*Setting outputs*/ - *nvertices = nconformvertices; - *nelements = nconformelements; - *px = newmeshX; - *py = newmeshY; + *pdata = newdata; + *pxy = newmeshXY; *pelements = newelements; /*Cleanup*/ @@ -692,37 +693,14 @@ return newgmesh; } /*}}}*/ -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"); - } - } - } /*}}}*/ void AdaptiveMeshRefinement::PrintGMeshVTK(TPZGeoMesh* gmesh,std::ofstream &file,bool matColor){/*{{{*/ Index: ../trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- ../trunk-jpl/src/c/classes/FemModel.cpp (revision 23064) +++ ../trunk-jpl/src/c/classes/FemModel.cpp (revision 23065) @@ -4876,7 +4876,9 @@ IssmDouble *newx = NULL; IssmDouble *newy = NULL; IssmDouble *newz = NULL; + IssmDouble *newxylist = NULL; int *newelementslist = NULL; + int* newdatalist = NULL; int newnumberofvertices = -1; int newnumberofelements = -1; @@ -4910,36 +4912,45 @@ } 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."); + this->amrbamg->ExecuteRefinementBamg(vector_serial,hmaxvertices_serial,&newdatalist,&newxylist,&newelementslist); + if(newdatalist[0]<=0 || newdatalist[1]<=0) _error_("Error in the refinement process."); } - /*Cleanup*/ - xDelete(vector_serial); - xDelete(hmaxvertices_serial); - delete vector; + /*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()); - /*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()); + /*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; } /*}}}*/ void FemModel::InitializeAdaptiveRefinementBamg(void){/*{{{*/ @@ -5189,7 +5200,6 @@ /*pnewelementslist keep vertices in Matlab indexing*/ int my_rank = IssmComm::GetRank(); - int numberofelements = this->elements->NumberOfElements(); IssmDouble* gl_distance = NULL; IssmDouble* if_distance = NULL; IssmDouble* deviatoricerror= NULL; @@ -5197,8 +5207,10 @@ IssmDouble* newx = NULL; IssmDouble* newy = NULL; IssmDouble* newz = NULL; + IssmDouble* newxylist = NULL; int* newelementslist = NULL; - int newnumberofvertices = -1; + int* newdatalist = NULL; + int newnumberofvertices = -1; int newnumberofelements = -1; /*Get fields, if requested*/ @@ -5209,38 +5221,46 @@ 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."); + &newdatalist,&newxylist,&newelementslist); + if(newdatalist[0]<=0 || newdatalist[1]<=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()); + /*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); } /*}}}*/ void FemModel::InitializeAdaptiveRefinementNeopz(void){/*{{{*/ Index: ../trunk-jpl/src/c/classes/AmrBamg.cpp =================================================================== --- ../trunk-jpl/src/c/classes/AmrBamg.cpp (revision 23064) +++ ../trunk-jpl/src/c/classes/AmrBamg.cpp (revision 23065) @@ -97,7 +97,7 @@ /*Cleanup and return*/ 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*/ BamgGeom* geomout=new BamgGeom(); @@ -149,17 +149,20 @@ this->previousmesh = meshout; /*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]); elementslist[3*i+1] = reCast(meshout->Triangles[4*i+1]); @@ -166,14 +169,13 @@ elementslist[3*i+2] = reCast(meshout->Triangles[4*i+2]); } + /*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){/*{{{*/