Changeset 23469

Timestamp:

11/17/18 11:53:24 (6 years ago)

Author:

tsantos

Message:

CHG: improving AMR (bamg) with error estimators, and optimizing AMR (neopz)

Location:

issm/trunk-jpl/src/c/classes

Files:

• AdaptiveMeshRefinement.cpp (modified) (10 diffs)
• AdaptiveMeshRefinement.h (modified) (1 diff)
• AmrBamg.h (modified) (1 diff)
• FemModel.cpp (modified) (2 diffs)

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

@@ -264 +264 @@
         if(verbose) _printf_(""<<count<<"\n");
         /*Adjust the connectivities before continue*/
-        gmesh->BuildConnectivity();
+        //gmesh->BuildConnectivity(); this is not necessary
         /*}}}*/
 
@@ -301 +301 @@
         if(verbose) _printf_(""<<count<<"\n");
         /*Adjust the connectivities before continue*/
-        gmesh->BuildConnectivity();
+        //gmesh->BuildConnectivity();//this is not necessary
         /*}}}*/
 
@@ -312 +312 @@
 }
 /*}}}*/
-int AdaptiveMeshRefinement::VerifyRefinementType(TPZGeoEl* geoel){/*{{{*/
+int AdaptiveMeshRefinement::VerifyRefinementType(TPZGeoEl* geoel,TPZGeoMesh* gmesh){/*{{{*/
 
         /*
@@ -323 +323 @@
         /*Output*/
         int type=0;
-        int count=0;
 
         /*Intermediaries*/
@@ -330 +329 @@
         /*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 +358 @@
 
         count=1;
+        
         while(count>0){
                 count=0;
@@ -367 +368 @@
                         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 +389 @@
                 }
                 /*Adjust the connectivities before continue*/
-                gmesh->BuildConnectivity();
+                //gmesh->BuildConnectivity();//this is not necessary
         }
 }
@@ -424 +427 @@
                         else _printf_(""<<count<<", ");
                 }
-                gmesh->BuildConnectivity();
+                //gmesh->BuildConnectivity();//this is not necessary
         }
 }
@@ -446 +449 @@
         this->specialelementsindex.clear();
         /*Adjust connectivities*/
-        gmesh->BuildConnectivity();
+        //gmesh->BuildConnectivity();//this is not necessary
 }
 /*}}}*/

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

@@ -90 +90 @@
         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);
         /*}}}*/
 };

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

@@ -37 +37 @@
                 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;

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

@@ -5106 +5106 @@
         int numberofelements                                            = this->elements->NumberOfElements();
         int numberofvertices                                            = this->vertices->NumberOfVertices();
+        IssmDouble      hmax                                                    = this->amrbamg->GetBamgOpts()->hmax;
         IssmDouble* maxlength                                   = xNew<IssmDouble>(numberofelements);
         IssmDouble* error_vertices                              = xNewZeroInit<IssmDouble>(numberofvertices);
@@ -5168 +5169 @@
         /*Fill hmaxvertices with the criteria*/
         for(int i=0;i<numberofelements;i++){
-                refine=false;
                 /*Refine any element if its error > 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;j<elementswidth;j++){
                                 vid=index[i*elementswidth+j]-1;//Matlab to C indexing
-                                if(error_vertices[vid]>groupthreshold*maxerror) refine=true;
-                        }
-                }
-                /*Now, fill the hmaxvertices if requested*/
-                if(refine){
-                        for(int j=0;j<elementswidth;j++){
-                                vid=index[i*elementswidth+j]-1;//Matlab to C indexing
-                                if(xIsNan<IssmDouble>(hmaxvertices[vid])) hmaxvertices[vid]=resolution;
-                                else hmaxvertices[vid]=min(resolution,hmaxvertices[vid]);
-                        }
+                                if(xIsNan<IssmDouble>(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;j<elementswidth;j++){
+                                        vid=index[i*elementswidth+j]-1;//Matlab to C indexing
+                                        if(error_vertices[vid]>groupthreshold*maxerror) hmaxvertices[vid]=maxlength[i]; //keep the current resolution
+                                        else{
+                                                if(xIsNan<IssmDouble>(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
+                                        }
+                                }
+                        } 
                 }
         }