Index: /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp (revision 23468) +++ /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp (revision 23469) @@ -264,5 +264,5 @@ if(verbose) _printf_(""<BuildConnectivity(); + //gmesh->BuildConnectivity(); this is not necessary /*}}}*/ @@ -301,5 +301,5 @@ if(verbose) _printf_(""<BuildConnectivity(); + //gmesh->BuildConnectivity();//this is not necessary /*}}}*/ @@ -312,5 +312,5 @@ } /*}}}*/ -int AdaptiveMeshRefinement::VerifyRefinementType(TPZGeoEl* geoel){/*{{{*/ +int AdaptiveMeshRefinement::VerifyRefinementType(TPZGeoEl* geoel,TPZGeoMesh* gmesh){/*{{{*/ /* @@ -323,5 +323,4 @@ /*Output*/ int type=0; - int count=0; /*Intermediaries*/ @@ -330,13 +329,14 @@ /*Loop over neighboors (sides 3, 4 and 5)*/ for(int j=3;j<6;j++){ + if(!gmesh->Element(geoel->NeighbourIndex(j))->HasSubElement()) continue; sons.clear(); - geoel->Neighbour(j).Element()->GetHigherSubElements(sons); - if(sons.size()) count++; //if neighbour was refined - if(sons.size()>4) count++; //if neighbour's level is > element level+1 - } - + gmesh->Element(geoel->NeighbourIndex(j))->GetHigherSubElements(sons); + if(sons.size()) type++; //if neighbour was refined + if(sons.size()>4) type++; //if neighbour's level is > element level+1 + if(type>1) break; + } + /*Verify and return*/ - if(count>1) type=2; - else type=count; + if(type>1) type=2; return type; @@ -358,4 +358,5 @@ count=1; + while(count>0){ count=0; @@ -367,15 +368,17 @@ if(gmesh->Element(i)->Level()==this->level_max) continue; /*loop over neighboors (sides 3, 4 and 5). Important: neighbours has the same dimension of the element*/ - type=this->VerifyRefinementType(gmesh->Element(i)); + type=this->VerifyRefinementType(gmesh->Element(i),gmesh); if(type<2){ typecount=0; for(int j=3;j<6;j++){ - if(gmesh->Element(i)->Neighbour(j).Element()->HasSubElement()) continue; - if(gmesh->Element(i)->Neighbour(j).Element()->Index()==i) typecount++;//neighbour==this element, element at the border - if(this->VerifyRefinementType(gmesh->Element(i)->Neighbour(j).Element())==1) typecount++; + if(gmesh->Element(gmesh->Element(i)->NeighbourIndex(j))->HasSubElement()) continue; + if(gmesh->Element(i)->NeighbourIndex(j)==i) typecount++;//neighbour==this element, element at the border + if(this->VerifyRefinementType(gmesh->Element(gmesh->Element(i)->NeighbourIndex(j)),gmesh)==1) typecount++; + if(typecount>1 && type==1) type=2; + else if(typecount>2 && type==0) type=2; + if(type==2) break; } - if(typecount>1 && type==1) type=2; - else if(typecount>2 && type==0) type=2; } + /*refine the element if requested*/ if(type==2){gmesh->Element(i)->Divide(sons); count++;} @@ -386,5 +389,5 @@ } /*Adjust the connectivities before continue*/ - gmesh->BuildConnectivity(); + //gmesh->BuildConnectivity();//this is not necessary } } @@ -424,5 +427,5 @@ else _printf_(""<BuildConnectivity(); + //gmesh->BuildConnectivity();//this is not necessary } } @@ -446,5 +449,5 @@ this->specialelementsindex.clear(); /*Adjust connectivities*/ - gmesh->BuildConnectivity(); + //gmesh->BuildConnectivity();//this is not necessary } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h =================================================================== --- /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h (revision 23468) +++ /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h (revision 23469) @@ -90,5 +90,5 @@ void PrintGMeshVTK(TPZGeoMesh *gmesh,std::ofstream &file,bool matColor=true); int GetVTK_ElType(TPZGeoEl* gel); - int VerifyRefinementType(TPZGeoEl* geoel); + int VerifyRefinementType(TPZGeoEl* geoel,TPZGeoMesh* gmesh); /*}}}*/ }; Index: /issm/trunk-jpl/src/c/classes/AmrBamg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/AmrBamg.h (revision 23468) +++ /issm/trunk-jpl/src/c/classes/AmrBamg.h (revision 23469) @@ -37,5 +37,8 @@ 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); - + + /*Access Method*/ + BamgOpts* GetBamgOpts(){return this->options;} + private: BamgGeom* geometry; Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 23468) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 23469) @@ -5106,4 +5106,5 @@ int numberofelements = this->elements->NumberOfElements(); int numberofvertices = this->vertices->NumberOfVertices(); + IssmDouble hmax = this->amrbamg->GetBamgOpts()->hmax; IssmDouble* maxlength = xNew(numberofelements); IssmDouble* error_vertices = xNewZeroInit(numberofvertices); @@ -5168,21 +5169,25 @@ /*Fill hmaxvertices with the criteria*/ for(int i=0;i phi*maxerror*/ - if(error_elements[i]>threshold*maxerror) refine=true; - /*If the element size is closer to the resolution, verify the sum of error in the vertices*/ - if(resolution/maxlength[i]>0.85){ + if(error_elements[i]>threshold*maxerror){ + /*Now, fill the hmaxvertices if requested*/ for(int j=0;jgroupthreshold*maxerror) refine=true; - } - } - /*Now, fill the hmaxvertices if requested*/ - if(refine){ - for(int j=0;j(hmaxvertices[vid])) hmaxvertices[vid]=resolution; - else hmaxvertices[vid]=min(resolution,hmaxvertices[vid]); - } + if(xIsNan(hmaxvertices[vid])) hmaxvertices[vid]=max(maxlength[i]/2.,resolution);//Try first dividing the element + else hmaxvertices[vid]=min(max(maxlength[i]/2.,resolution),hmaxvertices[vid]);//Try first dividing the element + } + } + else { + /*Try unrefine the element*/ + if(maxlength[i] < 1.1*hmax/2.){ + for(int j=0;jgroupthreshold*maxerror) hmaxvertices[vid]=maxlength[i]; //keep the current resolution + else{ + if(xIsNan(hmaxvertices[vid])) hmaxvertices[vid]=min(maxlength[i]*2.,hmax); + else hmaxvertices[vid]=min(min(maxlength[i]*2.,hmax),hmaxvertices[vid]);//Try first to duplicate the element + } + } + } } }