Index: /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp (revision 21502) +++ /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp (revision 21503) @@ -50,5 +50,6 @@ if(this->fathermesh) delete this->fathermesh; if(this->previousmesh) delete this->previousmesh; - + this->hmax=-1; + this->elementswidth=-1; } /*}}}*/ @@ -127,16 +128,15 @@ /*Mesh refinement methods*/ #include "TPZVTKGeoMesh.h" //itapopo -void AdaptiveMeshRefinement::ExecuteRefinement(int &type_process,double *vx, double *vy, double *masklevelset, long &nvertices, long &nelements, long &nsegments, double** x, double** y, double** z, long*** elements, long*** segments){/*{{{*/ - - //ITAPOPO _assert_(this->fathermesh); - //ITAPOPO _assert_(this->previousmesh); +void AdaptiveMeshRefinement::ExecuteRefinement(int &type_process,double *vx, double *vy, double *masklevelset, int &nvertices, int &nelements, int &nsegments, double** x, double** y, double** z, int*** elements, int*** segments){/*{{{*/ + + //itapopo _assert_(this->fathermesh); + //itapopo _assert_(this->previousmesh); /*Calculate the position of the grounding line using previous mesh*/ std::vector > GLvec; this->CalcGroundingLinePosition(masklevelset, GLvec); - - - //std::ofstream file1("/Users/santos/Desktop/mesh0.vtk"); - //TPZVTKGeoMesh::PrintGMeshVTK(this->previousmesh,file1 ); + + std::ofstream file1("/ronne_1/home/santos/mesh0.vtk"); + TPZVTKGeoMesh::PrintGMeshVTK(this->previousmesh,file1 ); /*run refinement or unrefinement process*/ @@ -149,7 +149,7 @@ this->RefinementProcess(newmesh,GLvec); - - //std::ofstream file2("/Users/santos/Desktop/mesh1.vtk"); - //TPZVTKGeoMesh::PrintGMeshVTK(this->previousmesh,file2 ); + + std::ofstream file2("/ronne_1/home/santos/mesh1.vtk"); + TPZVTKGeoMesh::PrintGMeshVTK(this->previousmesh,file2 ); /*Set new mesh pointer. Previous mesh just have uniform elements*/ @@ -162,23 +162,19 @@ TPZGeoMesh *nohangingnodesmesh = new TPZGeoMesh(*newmesh); - - //std::ofstream file3("/Users/santos/Desktop/mesh2.vtk"); - //TPZVTKGeoMesh::PrintGMeshVTK(newmesh,file3); + std::ofstream file3("/ronne_1/home/santos/mesh2.vtk"); + TPZVTKGeoMesh::PrintGMeshVTK(newmesh,file3); this->RefineMeshToAvoidHangingNodes(nohangingnodesmesh); - - //std::ofstream file4("/Users/santos/Desktop/mesh3.vtk"); - //TPZVTKGeoMesh::PrintGMeshVTK(nohangingnodesmesh,file4); - + std::ofstream file4("/ronne_1/home/santos/mesh3.vtk"); + TPZVTKGeoMesh::PrintGMeshVTK(nohangingnodesmesh,file4); /*Get new geometric mesh in ISSM data structure*/ this->GetMesh(nohangingnodesmesh,nvertices,nelements,nsegments,x,y,z,elements,segments); - + /*Verify the new geometry*/ - this->CheckMesh(nvertices,nelements,nsegments,this->elementswidth,(*x),(*y),(*z),(*elements),(*segments)); - + this->CheckMesh(nvertices,nelements,nsegments,this->elementswidth,x,y,z,elements,segments); + delete nohangingnodesmesh; - } /*}}}*/ @@ -189,5 +185,5 @@ /*Set elements to be refined using some criteria*/ - std::vector ElemVec; //elements without children + std::vector ElemVec; //elements without children this->SetElementsToRefine(gmesh,GLvec,hlevel,ElemVec); @@ -198,11 +194,11 @@ } /*}}}*/ -void AdaptiveMeshRefinement::RefineMesh(TPZGeoMesh *gmesh, std::vector &ElemVec){/*{{{*/ +void AdaptiveMeshRefinement::RefineMesh(TPZGeoMesh *gmesh, std::vector &ElemVec){/*{{{*/ /*Refine elements in ElemVec: uniform pattern refinement*/ - for(long i = 0; i < ElemVec.size(); i++){ + for(int i = 0; i < ElemVec.size(); i++){ /*Get geometric element and verify if it has already been refined*/ - long index = ElemVec[i]; + int index = ElemVec[i]; TPZGeoEl * geoel = gmesh->Element(index); if(geoel->HasSubElement()) DebugStop(); //itapopo _assert_(!geoel->HasSubElement()); @@ -211,5 +207,5 @@ /*Divide geoel*/ TPZVec Sons; - geoel->Divide(Sons); + geoel->Divide(Sons); /*If a 1D segment is neighbor, it must be divided too*/ @@ -236,6 +232,6 @@ /*Refine elements to avoid hanging nodes: non-uniform refinement*/ - const long NElem = gmesh->NElements(); - for(long i = 0; i < NElem; i++){ + const int NElem = gmesh->NElements(); + for(int i = 0; i < NElem; i++){ /*Get geometric element and verify if it has already been refined. Geoel may not have been previously refined*/ @@ -253,5 +249,5 @@ } - } + } gmesh->BuildConnectivity(); @@ -259,8 +255,8 @@ } /*}}}*/ -void AdaptiveMeshRefinement::GetMesh(TPZGeoMesh *gmesh, long &nvertices, long &nelements, long &nsegments, double** meshX, double** meshY, double** meshZ, long*** elements, long*** segments){/*{{{*/ +void AdaptiveMeshRefinement::GetMesh(TPZGeoMesh *gmesh, int &nvertices, int &nelements, int &nsegments, double** meshX, double** meshY, double** meshZ, int*** elements, int*** segments){/*{{{*/ /* vertices */ - long ntotalvertices = gmesh->NNodes();//total + int ntotalvertices = gmesh->NNodes();//total /* mesh coords */ @@ -269,6 +265,6 @@ double *newmeshZ = new double[ntotalvertices]; - /* getting mesh coords */ - for(long i = 0; i < ntotalvertices; i++ ){ + /* getting mesh coords */ + for(int i = 0; i < ntotalvertices; i++ ){ TPZVec coords(3,0.); gmesh->NodeVec()[i].GetCoordinates(coords); @@ -280,45 +276,40 @@ /* elements */ std::vector GeoVec; GeoVec.clear(); - for(long i = 0; i < gmesh->NElements(); i++){ - + for(int i = 0; i < gmesh->NElements(); i++){ if( gmesh->ElementVec()[i]->HasSubElement() ) continue; if( gmesh->ElementVec()[i]->MaterialId() != this->GetElemMaterialID() ) continue; GeoVec.push_back( gmesh->ElementVec()[i]); - - } - - long ntotalelements = GeoVec.size(); - long ** newelements = new long*[ntotalelements]; + } + + int ntotalelements = (int)GeoVec.size(); + int ** newelements = new int*[ntotalelements]; if ( !(this->elementswidth == 3) && !(this->elementswidth == 4) && !(this->elementswidth == 6) ) DebugStop(); - for(long i = 0; i < GeoVec.size(); i++){ - - newelements[i] = new long[this->elementswidth]; - for(int j = 0; j < this->elementswidth; j++) newelements[i][j] = GeoVec[i]->NodeIndex(j); + for(int i = 0; i < GeoVec.size(); i++){ + newelements[i] = new int[this->elementswidth]; + for(int j = 0; j < this->elementswidth; j++) newelements[i][j] = (int)GeoVec[i]->NodeIndex(j); } /* segments */ std::vector SegVec; SegVec.clear(); - for(long i = 0; i < gmesh->NElements(); i++){ - + for(int i = 0; i < gmesh->NElements(); i++){ if( gmesh->ElementVec()[i]->HasSubElement() ) continue; if( gmesh->ElementVec()[i]->MaterialId() != this->GetBoundaryMaterialID() ) continue; SegVec.push_back( gmesh->ElementVec()[i]); - - } - - long ntotalsegments = SegVec.size(); - long ** newsegments = new long*[ntotalsegments]; - - for(long i = 0; i < SegVec.size(); i++){ - - newsegments[i] = new long[3]; + } + + int ntotalsegments = (int)SegVec.size(); + int ** newsegments = new int*[ntotalsegments]; + + for(int i = 0; i < SegVec.size(); i++){ + + newsegments[i] = new int[3]; for(int j = 0; j < 2; j++) newsegments[i][j] = SegVec[i]->NodeIndex(j); - long neighborindex = SegVec[i]->Neighbour(2).Element()->Index(); - long neighbourid = -1; - - for(long j = 0; j < GeoVec.size(); j++){ + int neighborindex = SegVec[i]->Neighbour(2).Element()->Index(); + int neighbourid = -1; + + for(int j = 0; j < GeoVec.size(); j++){ if( GeoVec[j]->Index() == neighborindex || GeoVec[j]->FatherIndex() == neighborindex){ neighbourid = j; @@ -347,5 +338,5 @@ /* Find grounding line using elments center point */ GLvec.clear(); - for(long i=0;ipreviousmesh->NElements();i++){ + for(int i=0;ipreviousmesh->NElements();i++){ if(this->previousmesh->Element(i)->MaterialId()!=this->GetElemMaterialID()) continue; @@ -353,7 +344,7 @@ //itapopo apenas malha 2D triangular! - long vertex0 = this->previousmesh->Element(i)->NodeIndex(0); - long vertex1 = this->previousmesh->Element(i)->NodeIndex(1); - long vertex2 = this->previousmesh->Element(i)->NodeIndex(2); + int vertex0 = this->previousmesh->Element(i)->NodeIndex(0); + int vertex1 = this->previousmesh->Element(i)->NodeIndex(1); + int vertex2 = this->previousmesh->Element(i)->NodeIndex(2); double mls0 = masklevelset[vertex0]; @@ -384,5 +375,5 @@ } /*}}}*/ -void AdaptiveMeshRefinement::SetElementsToRefine(TPZGeoMesh *gmesh,std::vector > &GLvec,int &hlevel,std::vector &ElemVec){/*{{{*/ +void AdaptiveMeshRefinement::SetElementsToRefine(TPZGeoMesh *gmesh,std::vector > &GLvec,int &hlevel,std::vector &ElemVec){/*{{{*/ if(!gmesh) DebugStop(); //itapopo verificar se usará _assert_ @@ -398,9 +389,9 @@ } /*}}}*/ -void AdaptiveMeshRefinement::TagAllElements(TPZGeoMesh *gmesh,std::vector &ElemVec){/*{{{*/ +void AdaptiveMeshRefinement::TagAllElements(TPZGeoMesh *gmesh,std::vector &ElemVec){/*{{{*/ /* Uniform refinement. This refines the entire mesh */ - long nelements = gmesh->NElements(); - for(long i=0;iNElements(); + for(int i=0;iElement(i)->MaterialId()!=this->GetElemMaterialID()) continue; if(gmesh->Element(i)->HasSubElement()) continue; @@ -409,5 +400,5 @@ } /*}}}*/ -void AdaptiveMeshRefinement::TagElementsNearGroundingLine(TPZGeoMesh *gmesh,std::vector > &GLvec,int &hlevel,std::vector &ElemVec){/*{{{*/ +void AdaptiveMeshRefinement::TagElementsNearGroundingLine(TPZGeoMesh *gmesh,std::vector > &GLvec,int &hlevel,std::vector &ElemVec){/*{{{*/ /* Tag elements near grounding line */ @@ -416,5 +407,5 @@ double MaxDistance = MaxRegion / std::exp(alpha*(hlevel-1)); - for(long i=0;iNElements();i++){ + for(int i=0;iNElements();i++){ if(gmesh->Element(i)->MaterialId()!=this->GetElemMaterialID()) continue; @@ -430,5 +421,5 @@ REAL distance = MaxDistance; - for (long j = 0; j < GLvec.size(); j++) { + for (int j = 0; j < GLvec.size(); j++) { REAL value = ( GLvec[j][0] - centerPoint[0] ) * ( GLvec[j][0] - centerPoint[0] ); // (x2-x1)^2 @@ -446,18 +437,18 @@ } /*}}}*/ -void AdaptiveMeshRefinement::CreateInitialMesh(long &nvertices, long &nelements, long &nsegments, int &width, double* x, double* y, double* z, long** elements, long** segments){/*{{{*/ +void AdaptiveMeshRefinement::CreateInitialMesh(int &nvertices, int &nelements, int &nsegments, int &width, double* x, double* y, double* z, int** elements, int** segments){/*{{{*/ // itapopo _assert_(nvertices>0); // itapoo _assert_(nelements>0); - this->SetElementWidth(width); + this->SetElementWidth(width); /*Verify and creating initial mesh*/ //itapopo if(this->fathermesh) _error_("Initial mesh already exists!"); - this->fathermesh = new TPZGeoMesh(); + this->fathermesh = new TPZGeoMesh(); this->fathermesh->NodeVec().Resize( nvertices ); /*Set the vertices (geometric nodes in NeoPZ context)*/ - for(long i=0;ifathermesh->NodeVec()[i].SetCoord(coord); - this->fathermesh->NodeVec()[i].SetNodeId(i); + this->fathermesh->NodeVec()[i].SetNodeId(i); } @@ -477,5 +468,5 @@ TPZManVector elem(this->elementswidth,0); - for(long iel=0;ielelementswidth;jel++) elem[jel]=elements[iel][jel]; @@ -499,5 +490,5 @@ TPZManVector boundary(2,0.); - for(long iel=nelements;iel elemvertices; - elemvertices.clear(); - - for(long i = 0; i < nelements; i++){ - for(long j = 0; j < width; j++) { - elemvertices.insert(elements[i][j]); - } - } + std::set elemvertices; + elemvertices.clear(); + for(int i = 0; i < nelements; i++){ + for(int j = 0; j < width; j++) { + elemvertices.insert((*elements)[i][j]); + } + } if( elemvertices.size() != nvertices ) DebugStop();//itapopo verificar se irá usar o _assert_ - + //Verify if there are inf or NaN in coords - for(long i = 0; i < nvertices; i++) if(isnan(x[i]) || isinf(x[i])) DebugStop(); - for(long i = 0; i < nvertices; i++) if(isnan(y[i]) || isinf(y[i])) DebugStop(); - for(long i = 0; i < nvertices; i++) if(isnan(z[i]) || isinf(z[i])) DebugStop(); - for(long i = 0; i < nvertices; i++){ - for(long j = 0; j < width; j++){ - if( isnan(elements[i][j]) || isinf(elements[i][j]) ) DebugStop(); - } - } - -} -/*}}}*/ + for(int i = 0; i < nvertices; i++) if(isnan((*x)[i]) || isinf((*x)[i])) DebugStop(); + for(int i = 0; i < nvertices; i++) if(isnan((*y)[i]) || isinf((*y)[i])) DebugStop(); + for(int i = 0; i < nvertices; i++) if(isnan((*z)[i]) || isinf((*z)[i])) DebugStop(); + + for(int i = 0; i < nelements; i++){ + for(int j = 0; j < width; j++){ + if( isnan((*elements)[i][j]) || isinf((*elements)[i][j]) ) DebugStop(); + } + } + +} +/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h =================================================================== --- /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h (revision 21502) +++ /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h (revision 21503) @@ -55,7 +55,7 @@ void SetHMax(int &h); // Define the max level of refinement void SetElementWidth(int &width); // Define elements width - void ExecuteRefinement(int &type_process,double *vx, double *vy, double *masklevelset, long &nvertices, long &nelements, long &nsegments, double** x, double** y, double** z, long*** elements, long*** segments=NULL); // A new mesh will be created and refined. This returns the new mesh - void CreateInitialMesh(long &nvertices, long &nelements, long &nsegments, int &width, double* x, double* y, double* z, long** elements, long** segments=NULL); // Create a NeoPZ geometric mesh by coords and elements - void CheckMesh(long &nvertices, long &nelements, long &nsegments, int &width, double* x, double* y, double* z, long** elements, long** segments=NULL); // Check the consistency of the mesh + void ExecuteRefinement(int &type_process,double *vx, double *vy, double *masklevelset, int &nvertices, int &nelements, int &nsegments, double** x, double** y, double** z, int*** elements, int*** segments=NULL); // A new mesh will be created and refined. This returns the new mesh + void CreateInitialMesh(int &nvertices, int &nelements, int &nsegments, int &width, double* x, double* y, double* z, int** elements, int** segments=NULL); // Create a NeoPZ geometric mesh by coords and elements + void CheckMesh(int &nvertices, int &nelements, int &nsegments, int &width, double** x, double** y, double** z, int*** elements, int*** segments=NULL); // Check the consistency of the mesh private: @@ -69,11 +69,11 @@ /*Private methods*/ void RefinementProcess(TPZGeoMesh *gmesh,std::vector > &GLvec); // Start the refinement process - void RefineMesh(TPZGeoMesh *gmesh, std::vector &ElemVec); // Refine the elements in ElemVec + void RefineMesh(TPZGeoMesh *gmesh, std::vector &ElemVec); // Refine the elements in ElemVec void RefineMeshToAvoidHangingNodes(TPZGeoMesh *gmesh); // Refine the elements to avoid hanging nodes - void SetElementsToRefine(TPZGeoMesh *gmesh,std::vector > &GLvec,int &hlevel, std::vector &ElemVec); //Define wich elements will be refined - void TagAllElements(TPZGeoMesh *gmesh,std::vector &ElemVec); // This tag all elements to be refined, that is, refine all elements - void TagElementsNearGroundingLine(TPZGeoMesh *gmesh,std::vector > &GLvec,int &hlevel,std::vector &ElemVec); // This tag elements near the grounding line + void SetElementsToRefine(TPZGeoMesh *gmesh,std::vector > &GLvec,int &hlevel, std::vector &ElemVec); //Define wich elements will be refined + void TagAllElements(TPZGeoMesh *gmesh,std::vector &ElemVec); // This tag all elements to be refined, that is, refine all elements + void TagElementsNearGroundingLine(TPZGeoMesh *gmesh,std::vector > &GLvec,int &hlevel,std::vector &ElemVec); // This tag elements near the grounding line void CalcGroundingLinePosition(double *masklevelset,std::vector > &GLvec); // Calculate the grounding line position using previous mesh - void GetMesh(TPZGeoMesh *gmesh, long &nvertices, long &nelements, long &nsegments, double** meshX, double** meshY, double** meshZ, long*** elements, long*** segments=NULL); // Return coords and elements in ISSM data structure + void GetMesh(TPZGeoMesh *gmesh, int &nvertices, int &nelements, int &nsegments, double** meshX, double** meshY, double** meshZ, int*** elements, int*** segments=NULL); // Return coords and elements in ISSM data structure inline int GetElemMaterialID(){return 1;} // Return element material ID inline int GetBoundaryMaterialID(){return 2;} // Return segment (2D boundary) material ID