![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 21919 for issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp
- Timestamp:
- 08/09/17 10:26:20 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp
r21806 r21919 10 10 11 11 #include "./AdaptiveMeshRefinement.h" 12 #include "TPZVTKGeoMesh.h" 13 #include "../shared/shared.h" 14 #include "pzgeotriangle.h" 15 #include "pzreftriangle.h" 12 16 13 using namespace pzgeom; 17 14 … … 40 37 this->sid2index.resize(cp.sid2index.size()); 41 38 for(int i=0;i<cp.sid2index.size();i++) this->sid2index[i]=cp.sid2index[i]; 42 39 this->index2sid.clear(); 40 this->index2sid.resize(cp.index2sid.size()); 41 for(int i=0;i<cp.index2sid.size();i++) this->index2sid[i]=cp.index2sid[i]; 42 this->specialelementsindex.clear(); 43 this->specialelementsindex.resize(cp.specialelementsindex.size()); 44 for(int i=0;i<cp.specialelementsindex.size();i++) this->specialelementsindex[i]=cp.specialelementsindex[i]; 45 43 46 return *this; 44 45 47 } 46 48 /*}}}*/ 47 49 AdaptiveMeshRefinement::~AdaptiveMeshRefinement(){/*{{{*/ 48 49 bool ismismip = false;50 if(ismismip){//itapopo51 TPZFileStream fstr;52 std::stringstream ss;53 54 ss << this->levelmax;55 std::string AMRfile = "/home/santos/L" + ss.str() + "_amr.txt";56 57 fstr.OpenWrite(AMRfile.c_str());58 int withclassid = 1;59 this->Write(fstr,withclassid);60 }61 50 this->CleanUp(); 62 51 gRefDBase.clear(); … … 73 62 this->regionlevelmax = -1; 74 63 this->sid2index.clear(); 64 this->index2sid.clear(); 65 this->specialelementsindex.clear(); 75 66 } 76 67 /*}}}*/ … … 84 75 this->regionlevel1 = -1; 85 76 this->regionlevelmax = -1; 77 this->step = 0; 78 this->smooth_frequency = 1; 86 79 this->sid2index.clear(); 87 } 88 /*}}}*/ 89 int AdaptiveMeshRefinement::ClassId() const{/*{{{*/ 90 return 13829430; //Antartic area with ice shelves (km^2) 91 } 92 /*}}}*/ 93 void AdaptiveMeshRefinement::Read(TPZStream &buf,void *context){/*{{{*/ 94 95 try 96 { 97 /* Read the id context*/ 98 TPZSaveable::Read(buf,context); 99 /* Read class id*/ 100 int classid; 101 buf.Read(&classid,1); 102 /* Verify the class id*/ 103 if(classid!=this->ClassId()) _error_("AdaptiveMeshRefinement::Read: Error in restoring AdaptiveMeshRefinement!\n"); 104 /* Read simple attributes */ 105 buf.Read(&this->levelmax,1); 106 buf.Read(&this->elementswidth,1); 107 buf.Read(&this->regionlevel1,1); 108 buf.Read(&this->regionlevelmax,1); 109 /* Read vector attributes*/ 110 int size; 111 buf.Read(&size,1); 112 int* psid2index=xNew<int>(size); 113 buf.Read(psid2index,size); 114 this->sid2index.clear(); 115 this->sid2index.assign(psid2index,psid2index+size); 116 /* Read geometric mesh (father)*/ 117 TPZSaveable *sv1 = TPZSaveable::Restore(buf,0); 118 this->fathermesh = dynamic_cast<TPZGeoMesh*>(sv1); 119 /* Read geometric mesh (current)*/ 120 TPZSaveable *sv2 = TPZSaveable::Restore(buf,0); 121 this->currentmesh = dynamic_cast<TPZGeoMesh*>(sv2); 122 /* Cleanup*/ 123 xDelete<int>(psid2index); 124 } 125 catch(const std::exception& e) 126 { 127 _error_("AdaptiveMeshRefinement::Read: Exception catched!\n"); 128 } 129 130 } 131 /*}}}*/ 132 template class TPZRestoreClass<AdaptiveMeshRefinement,13829430>;/*{{{*/ 133 /*}}}*/ 134 void AdaptiveMeshRefinement::Write(TPZStream &buf,int withclassid){/*{{{*/ 135 136 try 137 { 138 /* Write context (this class) class ID*/ 139 TPZSaveable::Write(buf,withclassid); 140 /* Write this class id*/ 141 int classid = this->ClassId(); 142 buf.Write(&classid,1); 143 /* Write simple attributes */ 144 buf.Write(&this->levelmax,1); 145 buf.Write(&this->elementswidth,1); 146 buf.Write(&this->regionlevel1,1); 147 buf.Write(&this->regionlevelmax,1); 148 /* Write vector attributes*/ 149 int size=this->sid2index.size(); 150 int* psid2index=&this->sid2index[0]; 151 buf.Write(&size,1);//vector size 152 buf.Write(psid2index,this->sid2index.size()); 153 /* Write the geometric mesh*/ 154 this->fathermesh->Write(buf,this->ClassId()); 155 this->currentmesh->Write(buf,this->ClassId()); 156 } 157 catch(const std::exception& e) 158 { 159 _error_("AdaptiveMeshRefinement::Write: Exception catched!\n"); 160 } 80 this->index2sid.clear(); 81 this->specialelementsindex.clear(); 161 82 } 162 83 /*}}}*/ … … 176 97 if(numberofelements!=this->sid2index.size()) _error_("Impossible to execute refinement: sid2index.size is not equal to numberofelements!\n"); 177 98 178 /*Execute the refinement.*/ 179 this->RefinementProcess(amr_verbose,partiallyfloatedelements,masklevelset,deviatorictensorerror,thicknesserror); 180 181 /*Get new geometric mesh in ISSM data structure*/ 182 this->GetMesh(newnumberofvertices,newnumberofelements,x,y,elementslist); 183 184 /*Verify the new geometry*/ 185 this->CheckMesh(newnumberofvertices,newnumberofelements,this->elementswidth,x,y,elementslist); 186 187 } 188 /*}}}*/ 189 void AdaptiveMeshRefinement::RefinementProcess(bool &amr_verbose,double* partiallyfloatedelements,double* masklevelset, 190 double* deviatorictensorerror,double* thicknesserror){/*{{{*/ 191 192 if(amr_verbose) _printf_("\n\trefinement process started (level max = " << this->levelmax << ")\n"); 193 194 /*Intermediaries*/ 195 TPZGeoMesh* nohangingnodesmesh=NULL; 99 /*Combine the fields*/ 196 100 double mean_mask = 0; 197 101 double mean_tauerror = 0; 198 102 double mean_Herror = 0; 199 double group_error = 0; 200 103 int* fmask = xNew<int>(numberofelements); 104 int* ftauerror = xNew<int>(numberofelements); 105 int* fHerror = xNew<int>(numberofelements); 106 int* phi = xNew<int>(numberofelements); 201 107 /*Calculate mean values*/ 202 108 for(int i=0;i<this->sid2index.size();i++){ … … 208 114 mean_tauerror /= this->sid2index.size(); 209 115 mean_Herror /= this->sid2index.size(); 210 116 /*Filter to thickness*/ 117 for(int i=0;i<this->sid2index.size();i++){ 118 fmask[i]=0; 119 if(thicknesserror[i]>mean_Herror) fmask[i]=1; 120 } 121 /*Filter to tau*/ 122 for(int i=0;i<this->sid2index.size();i++){ 123 ftauerror[i]=0; 124 if(deviatorictensorerror[i]>mean_tauerror) ftauerror[i]=1; 125 } 126 /*Sum*/ 127 for(int i=0;i<this->sid2index.size();i++){ 128 phi[i]=ftauerror[i]+fmask[i]; 129 } 130 131 /*Execute the refinement.*/ 132 this->RefinementProcess(amr_verbose,partiallyfloatedelements,masklevelset,deviatorictensorerror,thicknesserror); 133 134 /*Get new geometric mesh in ISSM data structure*/ 135 this->GetMesh(newnumberofvertices,newnumberofelements,x,y,elementslist); 136 137 std::ofstream file("/home/santos/mesh.vtk"); 138 TPZVTKGeoMesh::PrintGMeshVTK(this->currentmesh,file); 139 140 /*Verify the new geometry*/ 141 this->CheckMesh(newnumberofvertices,newnumberofelements,this->elementswidth,x,y,elementslist); 142 143 } 144 /*}}}*/ 145 void AdaptiveMeshRefinement::RefinementProcess(bool &amr_verbose,double* partiallyfloatedelements,double* masklevelset,double* deviatorictensorerror,double* thicknesserror){/*{{{*/ 146 147 if(amr_verbose) _printf_("\n\trefinement process started (level max = " << this->levelmax << ")\n"); 148 149 /*Intermediaries*/ 150 TPZGeoMesh* nohangingnodesmesh=NULL; 151 152 //itapopo 153 this->step++; 154 155 double mean_mask = 0; 156 double mean_tauerror = 0; 157 double mean_Herror = 0; 158 double group_Herror = 0; 159 int index,sid; 160 std::vector<int> specialfatherindex; specialfatherindex.clear(); 161 162 /*Calculate mean values*/ 163 for(int i=0;i<this->sid2index.size();i++){ 164 mean_mask += masklevelset[i]; 165 mean_tauerror += deviatorictensorerror[i]; 166 mean_Herror += thicknesserror[i]; 167 } 168 mean_mask /= this->sid2index.size(); 169 mean_tauerror /= this->sid2index.size(); 170 mean_Herror /= this->sid2index.size(); 171 172 if(amr_verbose) _printf_("\t\tdeal with special elements...\n"); 173 /*Deal with special elements*/ 174 for(int i=0;i<this->specialelementsindex.size();i++){ 175 if(this->specialelementsindex[i]==-1) continue; 176 /*Get special element and verify*/ 177 TPZGeoEl* geoel=this->currentmesh->Element(this->specialelementsindex[i]); 178 if(!geoel)_error_("special element (sid) "<<i<<" is null!\n"); 179 if(geoel->HasSubElement())_error_("special element (sid) "<<i<<" has "<<geoel->NSubElements()<<" subelements!\n"); 180 if(geoel->MaterialId()!=this->GetElemMaterialID()) _error_("geoel->MaterialId is not GetElemMaterialID!\n"); 181 if(!geoel->Father())_error_("father of special element (sid) "<<i<<" is null!\n"); 182 183 /*Get element's siblings and verify*/ 184 TPZGeoEl* father=geoel->Father(); 185 TPZVec<TPZGeoEl *> siblings; 186 father->GetHigherSubElements(siblings); 187 std::vector<int> sidvec; sidvec.resize(siblings.size()); 188 for (int j=0;j<siblings.size();j++){ 189 if(!siblings[j]) _error_("special element (sid) "<<i<<" has a null siblings null!\n"); 190 sidvec[j]=this->index2sid[siblings[j]->Index()]; 191 } 192 193 /*Now, reset the data strucure and verify if the siblings should be deleted*/ 194 if(siblings.size()<4){ 195 /*Reset subelements in the father*/ 196 father->ResetSubElements(); 197 }else{ 198 if(siblings.size()!=4) _error_("element (index) "<<father->Index()<<" has "<<father->NSubElements()<<" subelements!\n"); 199 } 200 for (int j=0;j<siblings.size();j++){ 201 for(int k=0;k<this->specialelementsindex.size();k++){ 202 if(this->specialelementsindex[k]==siblings[j]->Index()){ 203 index = siblings[j]->Index(); 204 if(index<0) _error_("index is null!\n"); 205 sid = this->index2sid[index]; 206 if(sid<0) _error_("sid is null!\n"); 207 this->specialelementsindex[k] = -1; 208 this->index2sid[index] = -1; 209 this->sid2index[sid] = -1; 210 } 211 } 212 if(siblings.size()<4){ 213 /*Ok, the special element can be deleted*/ 214 siblings[j]->ResetSubElements(); 215 this->currentmesh->DeleteElement(siblings[j],siblings[j]->Index()); 216 } 217 } 218 219 /*Now, verify if the father should be refined with uniform pattern (smoother)*/ 220 this->smooth_frequency=10000000; 221 if(siblings.size()==3 || this->step%this->smooth_frequency==0){//it keeps the mesh with uniform elements 222 /*Father has uniform subelements now*/ 223 TPZVec<TPZGeoEl *> sons; 224 father->Divide(sons); 225 this->smooth_frequency=this->step; 226 }else{ 227 specialfatherindex.push_back(father->Index()); 228 } 229 if(this->specialelementsindex[i]!=-1) _error_("special element "<<i<<" was not deleted!\n"); 230 } 231 this->currentmesh->BuildConnectivity(); 232 211 233 if(amr_verbose) _printf_("\t\tuniform refinement...\n"); 234 /*Deal with uniform elemnts*/ 212 235 for(int i=0;i<this->sid2index.size();i++){ 236 if(this->sid2index[i]==-1) continue; 237 /*Get element and verify*/ 213 238 TPZGeoEl* geoel=this->currentmesh->Element(this->sid2index[i]); 214 if(geoel->HasSubElement()) _error_("Impossible to refine: geoel (index) "<<this->sid2index[i]<<" has subelements!\n"); 215 if(geoel->MaterialId()!=this->GetElemMaterialID()) _error_("Impossible to refine: geoel->MaterialId is not GetElemMaterialID!\n"); 216 /*Refine*/ 217 if(thicknesserror[i]>mean_Herror){ 239 if(geoel->HasSubElement()) _error_("element (sid) "<<i<<" has "<<geoel->NSubElements()<<" subelements!\n"); 240 if(geoel->MaterialId()!=this->GetElemMaterialID()) _error_("geoel->MaterialId is not GetElemMaterialID!\n"); 241 242 /*Refine process*/ 243 if(thicknesserror[i]>mean_Herror) 244 { 245 int count=0; 246 TPZGeoEl* father=geoel->Father(); 247 if(father){ 248 for(int j=3;j<6;j++){ 249 index=father->Neighbour(j).Element()->Index(); 250 for(int k=0;k<specialfatherindex.size();k++) if(specialfatherindex[k]==index) count++; 251 } 252 } 218 253 TPZVec<TPZGeoEl *> sons; 219 if(geoel->Level()<this->levelmax ) geoel->Divide(sons);254 if(geoel->Level()<this->levelmax && count==0) geoel->Divide(sons); 220 255 } 221 else if(geoel->Level()>0){ /*try to unrefine*/ 222 TPZVec<TPZGeoEl *> sons; 223 geoel->Father()->GetHigherSubElements(sons); 224 group_error=0; 225 for(int j=0;j<sons.size();j++){ 226 sons[j]->Index(); 256 else if(geoel->Level()>0) 257 {/*Unrefine process*/ 258 259 /*Get siblings and verify*/ 260 TPZVec<TPZGeoEl *> siblings; 261 geoel->Father()->GetHigherSubElements(siblings); 262 //if(siblings.size()<4) _error_("Impossible to refine: geoel (index) "<<this->sid2index[i]<<" has less than 3 siblings!\n"); 263 if(siblings.size()>4) continue;//Father has more then 4 sons, this group should not be unrefined. 264 265 /*Compute the error of the group*/ 266 group_Herror=0; 267 for(int j=0;j<siblings.size();j++){ 268 index = siblings[j]->Index(); 269 sid = this->index2sid[index]; 270 if(sid==-1) continue; 271 group_Herror+=thicknesserror[sid]; 227 272 } 228 } 229 } 273 /*Verify if this group should be unrefined*/ 274 if(group_Herror>0 && group_Herror<0*mean_Herror){ //itapopo 275 /*Reset subelements in the father*/ 276 this->currentmesh->Element(geoel->Father()->Index())->ResetSubElements(); 277 /*Delete the elements and set their indexes in the index2sid and sid2index*/ 278 for (int j=0;j<siblings.size();j++){ 279 index = siblings[j]->Index(); 280 sid = this->index2sid[index]; 281 this->index2sid[index]=-1; 282 if(sid!=-1) this->sid2index[sid]=-1; 283 this->currentmesh->DeleteElement(siblings[j],siblings[j]->Index()); 284 }//for j 285 }//if 286 }/*Unrefine process*/ 287 }//for i 230 288 this->currentmesh->BuildConnectivity(); 231 289 232 290 if(amr_verbose) _printf_("\t\trefine to avoid hanging nodes...\n"); 233 291 this->RefineMeshToAvoidHangingNodes(this->currentmesh); 234 this->currentmesh->BuildConnectivity();235 292 236 293 //nohangingnodesmesh = this->CreateRefPatternMesh(newmesh); itapopo tentar otimizar 237 294 238 295 if(amr_verbose) _printf_("\trefinement process done!\n"); 239 296 } … … 257 314 void AdaptiveMeshRefinement::RefineMeshToAvoidHangingNodes(TPZGeoMesh *gmesh){/*{{{*/ 258 315 259 /*Refine elements to avoid hanging nodes: non-uniform refinement*/ 316 /*Now, insert special elements to avoid hanging nodes*/ 317 this->specialelementsindex.clear(); 260 318 const int NElem = gmesh->NElements(); 261 319 for(int i=0;i<NElem;i++){ … … 268 326 TPZAutoPointer<TPZRefPattern> refp=TPZRefPatternTools::PerfectMatchRefPattern(geoel); 269 327 if(refp){ 328 /*Non-uniform refinement*/ 270 329 TPZVec<TPZGeoEl *> Sons; 271 330 geoel->SetRefPattern(refp); 272 331 geoel->Divide(Sons); 273 } 274 } 275 gmesh->BuildConnectivity(); 276 332 /*Keep the index of the special elements*/ 333 for(int j=0;j<Sons.size();j++) this->specialelementsindex.push_back(Sons[j]->Index()); 334 } 335 } 336 gmesh->BuildConnectivity(); 277 337 } 278 338 /*}}}*/ … … 281 341 /* IMPORTANT! pelements are in Matlab indexing 282 342 NEOPZ works only in C indexing. 283 This method cleans up and updated the this->sid2index 284 and fills in it with the new mesh.343 This method cleans up and updated the this->sid2index 344 and this->index2sid and fills in it with the new mesh. 285 345 Avoid to call this method before Refinement Process.*/ 286 346 … … 295 355 TPZGeoEl* geoel = NULL; 296 356 long* vertex_index2sid = xNew<long>(gmesh->NNodes()); 357 this->index2sid.clear(); this->index2sid.resize(gmesh->NElements()); 297 358 this->sid2index.clear(); 298 359 … … 307 368 /*Fill in the vertex_index2sid vector with non usual index value*/ 308 369 for(int i=0;i<gmesh->NNodes();i++) vertex_index2sid[i]=-1; 370 371 /*Fill in the this->index2sid vector with non usual index value*/ 372 for(int i=0;i<gmesh->NElements();i++) this->index2sid[i]=-1; 309 373 310 374 /*Get elements without sons and fill in the vertex_index2sid with used vertices (indexes) */ … … 316 380 if(geoel->MaterialId() != this->GetElemMaterialID()) continue; 317 381 this->sid2index.push_back(i);//keep the element index 382 this->index2sid[i]=this->sid2index.size()-1;//keep the element sid 318 383 for(int j=0;j<this->elementswidth;j++){ 319 384 nodeindex=geoel->NodeIndex(j); … … 347 412 newelements[i*this->elementswidth+j]=(int)sid+1;//C to Matlab indexing 348 413 } 349 } 350 414 /*Verify the Jacobian determinant. If detJ<0, swap the 2 first postions: 415 a -> b 416 b -> a */ 417 double detJ,xa,xb,xc,ya,yb,yc; 418 int a,b,c; 419 420 a=newelements[i*this->elementswidth+0]-1; 421 b=newelements[i*this->elementswidth+1]-1; 422 c=newelements[i*this->elementswidth+2]-1; 423 424 xa=newmeshX[a]; ya=newmeshY[a]; 425 xb=newmeshX[b]; yb=newmeshY[b]; 426 xc=newmeshX[c]; yc=newmeshY[c]; 427 428 detJ=(xb-xa)*(yc-ya)-(xc-xa)*(yb-ya); 429 430 /*verify and swap, if necessary*/ 431 if(detJ<0) { 432 newelements[i*this->elementswidth+0]=b+1;//a->b 433 newelements[i*this->elementswidth+1]=a+1;//b->a 434 } 435 } 436 351 437 /*Setting outputs*/ 352 438 nvertices = nconformvertices; … … 358 444 /*Cleanup*/ 359 445 xDelete<long>(vertex_index2sid); 360 361 446 } 362 447 /*}}}*/
Note:
See TracChangeset
for help on using the changeset viewer.
