source: issm/trunk/src/c/classes/AdaptiveMeshRefinement.cpp@ 24313

Last change on this file since 24313 was 24313, checked in by Mathieu Morlighem, 5 years ago

merged trunk-jpl and trunk for revision 24310
File size: 34.6 KB
Line 
1/*!\file AdaptiveMeshrefinement.cpp
2 * \brief: implementation of the adaptive mesh refinement tool based on NeoPZ library: github.com/labmec/neopz
3 */
4
5#ifdef HAVE_CONFIG_H
6 #include <config.h>
7#else
8#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
9#endif
10
11#include "./AdaptiveMeshRefinement.h"
12
13using namespace pzgeom;
14
15/*Constructor, copy, clean up and destructor*/
16AdaptiveMeshRefinement::AdaptiveMeshRefinement(){/*{{{*/
17
18 /*Set pointers to NULL*/
19 this->fathermesh = NULL;
20 this->previousmesh = NULL;
21 this->refinement_type = -1;
22 this->level_max = -1;
23 this->gradation = -1;
24 this->lag = -1;
25 this->groundingline_distance = -1;
26 this->icefront_distance = -1;
27 this->thicknesserror_threshold = -1;
28 this->deviatoricerror_threshold = -1;
29 this->deviatoricerror_maximum = -1;
30 this->thicknesserror_maximum = -1;
31 this->sid2index.clear();
32 this->index2sid.clear();
33 this->specialelementsindex.clear();
34 this->x = NULL;
35 this->y = NULL;
36 this->elementslist = NULL;
37 this->numberofvertices = -1;
38 this->numberofelements = -1;
39}
40/*}}}*/
41AdaptiveMeshRefinement::AdaptiveMeshRefinement(const AdaptiveMeshRefinement &cp){/*{{{*/
42 this->Initialize();
43 this->operator =(cp);
44}
45/*}}}*/
46AdaptiveMeshRefinement & AdaptiveMeshRefinement::operator =(const AdaptiveMeshRefinement &cp){/*{{{*/
47
48 /*Clean all attributes*/
49 this->CleanUp();
50 /*Copy all data*/
51 this->fathermesh = new TPZGeoMesh(*cp.fathermesh);
52 this->previousmesh = new TPZGeoMesh(*cp.previousmesh);
53 this->refinement_type = cp.refinement_type;
54 this->level_max = cp.level_max;
55 this->gradation = cp.gradation;
56 this->lag = cp.lag;
57 this->groundingline_distance = cp.groundingline_distance;
58 this->icefront_distance = cp.icefront_distance;
59 this->thicknesserror_threshold = cp.thicknesserror_threshold;
60 this->deviatoricerror_threshold = cp.deviatoricerror_threshold;
61 this->deviatoricerror_maximum = cp.deviatoricerror_maximum;
62 this->thicknesserror_maximum = cp.thicknesserror_maximum;
63 this->sid2index.clear();
64 this->sid2index.resize(cp.sid2index.size());
65 for(int i=0;i<cp.sid2index.size();i++) this->sid2index[i]=cp.sid2index[i];
66 this->index2sid.clear();
67 this->index2sid.resize(cp.index2sid.size());
68 for(int i=0;i<cp.index2sid.size();i++) this->index2sid[i]=cp.index2sid[i];
69 this->specialelementsindex.clear();
70 this->specialelementsindex.resize(cp.specialelementsindex.size());
71 for(int i=0;i<cp.specialelementsindex.size();i++) this->specialelementsindex[i]=cp.specialelementsindex[i];
72
73 return *this;
74}
75/*}}}*/
76AdaptiveMeshRefinement::~AdaptiveMeshRefinement(){/*{{{*/
77 int writemesh = 0;//only to restart
78 if(writemesh) this->WriteMesh();
79 this->CleanUp();
80 gRefDBase.clear();
81}
82/*}}}*/
83void AdaptiveMeshRefinement::CleanUp(){/*{{{*/
84
85 /*Verify and delete all data*/
86 if(this->fathermesh) delete this->fathermesh;
87 if(this->previousmesh) delete this->previousmesh;
88 if(this->x) delete this->x;
89 if(this->y) delete this->y;
90 if(this->elementslist) delete this->elementslist;
91 this->refinement_type = -1;
92 this->level_max = -1;
93 this->gradation = -1;
94 this->lag = -1;
95 this->groundingline_distance = -1;
96 this->icefront_distance = -1;
97 this->thicknesserror_threshold = -1;
98 this->deviatoricerror_threshold = -1;
99 this->deviatoricerror_maximum = -1;
100 this->thicknesserror_maximum = -1;
101 this->numberofvertices = -1;
102 this->numberofelements = -1;
103 this->sid2index.clear();
104 this->index2sid.clear();
105 this->specialelementsindex.clear();
106}
107/*}}}*/
108
109/*Mesh refinement methods*/
110void AdaptiveMeshRefinement::SetMesh(int** elementslist_in,IssmDouble** x_in,IssmDouble** y_in,int* numberofvertices_in,int* numberofelements_in){/*{{{*/
111
112 /*Delete previous mesh and keep the entire mesh*/
113 if(this->elementslist) xDelete<int>(this->elementslist);
114 if(this->x) xDelete<IssmDouble>(this->x);
115 if(this->y) xDelete<IssmDouble>(this->y);
116
117 this->elementslist = *elementslist_in;
118 this->x = *x_in;
119 this->y = *y_in;
120 this->numberofvertices = *numberofvertices_in;
121 this->numberofelements = *numberofelements_in;
122}/*}}}*/
123void AdaptiveMeshRefinement::GetMesh(int** elementslist_out,IssmDouble** x_out,IssmDouble** y_out,int* numberofvertices_out,int* numberofelements_out){/*{{{*/
124
125 /*Get the entire mesh*/
126 *elementslist_out = this->elementslist;
127 *x_out = this->x;
128 *y_out = this->y;
129 *numberofvertices_out= this->numberofvertices;
130 *numberofelements_out= this->numberofelements;
131}/*}}}*/
132void AdaptiveMeshRefinement::ExecuteRefinement(double* gl_distance,double* if_distance,double* deviatoricerror,double* thicknesserror,int** pdatalist,double** pxylist,int** pelementslist){/*{{{*/
133
134 /*IMPORTANT! pelementslist are in Matlab indexing*/
135 /*NEOPZ works only in C indexing*/
136 if(!this->fathermesh || !this->previousmesh) _error_("Impossible to execute refinement: fathermesh or previousmesh is NULL!\n");
137 if(this->refinement_type!=0 && this->refinement_type!=1) _error_("Impossible to execute refinement: refinement type is not defined!\n");
138
139 /*Input verifications*/
140 if(this->deviatoricerror_threshold>0 && !deviatoricerror) _error_("deviatoricerror is NULL!\n");
141 if(this->thicknesserror_threshold>0 && !thicknesserror) _error_("thicknesserror is NULL!\n");
142 if(this->groundingline_distance>0 && !gl_distance) _error_("gl_distance is NULL!\n");
143 if(this->icefront_distance>0 && !if_distance) _error_("if_distance is NULL!\n");
144 /*Attributes verifications*/
145 if(this->deviatoricerror_threshold>0 && this->deviatoricerror_groupthreshold<DBL_EPSILON) _error_("group threshold is too small!");
146 if(this->thicknesserror_threshold>0 && this->thicknesserror_groupthreshold<DBL_EPSILON) _error_("group threshold is too small!");
147
148 /*Intermediaries*/
149 bool verbose=VerboseSolution();
150
151 /*Execute refinement*/
152 this->RefineMeshOneLevel(verbose,gl_distance,if_distance,deviatoricerror,thicknesserror);
153
154 /*Get new geometric mesh in ISSM data structure*/
155 this->GetMesh(this->previousmesh,pdatalist,pxylist,pelementslist);
156
157 /*Verify the new geometry*/
158 this->CheckMesh(pdatalist,pxylist,pelementslist);
159}
160/*}}}*/
161void AdaptiveMeshRefinement::RefineMeshOneLevel(bool &verbose,double* gl_distance,double* if_distance,double* deviatoricerror,double* thicknesserror){/*{{{*/
162
163 /*Intermediaries*/
164 int nelem =-1;
165 int side2D = 6;
166 int sid =-1;
167 int count =-1;
168 int criteria =-1;
169 int numberofcriteria =-1;
170 int nconformelements = this->sid2index.size();
171 double gl_distance_h =-1;
172 double gl_distance_hmax = this->groundingline_distance;
173 double if_distance_h =-1;
174 double if_distance_hmax = this->icefront_distance;
175 double gl_groupdistance =-1;
176 double if_groupdistance =-1;
177 double d_maxerror =-1;
178 double t_maxerror =-1;
179 double deviatoric_grouperror =-1;
180 double thickness_grouperror =-1;
181 TPZGeoMesh* gmesh = NULL;
182 TPZVec<REAL> qsi(2,0.),cp(3,0.);
183 TPZVec<TPZGeoEl*> sons;
184 std::vector<int> index;
185
186 /*Calculate the number of criteria{{{*/
187 numberofcriteria=0;
188 if(this->deviatoricerror_threshold>0) numberofcriteria++;
189 if(this->thicknesserror_threshold>0) numberofcriteria++;
190 if(this->groundingline_distance>0) numberofcriteria++;
191 if(this->icefront_distance>0) numberofcriteria++;
192 /*}}}*/
193
194 /*Calculate the maximum of the estimators, if requested{{{*/
195 if(this->deviatoricerror_threshold>0 && this->deviatoricerror_maximum<DBL_EPSILON){
196 for(int i=0;i<nconformelements;i++) this->deviatoricerror_maximum=max(this->deviatoricerror_maximum,deviatoricerror[i]);
197 }
198 if(this->thicknesserror_threshold>0 && this->thicknesserror_maximum<DBL_EPSILON){
199 for(int i=0;i<nconformelements;i++) this->thicknesserror_maximum=max(this->thicknesserror_maximum,thicknesserror[i]);
200 }
201 /*}}}*/
202
203 /*First, verify if special elements have min distance or high errors{{{*/
204 gmesh=this->previousmesh;
205 for(int i=0;i<this->specialelementsindex.size();i++){
206 if(this->specialelementsindex[i]==-1) _error_("index is -1!\n");
207 if(!gmesh->Element(this->specialelementsindex[i])) _error_("element is null!\n");
208 if(!gmesh->Element(this->specialelementsindex[i])->Father()) _error_("father is null!\n");
209 if(gmesh->Element(this->specialelementsindex[i])->HasSubElement()) _error_("special element has sub elements!\n");
210 sons.clear();
211 gmesh->Element(this->specialelementsindex[i])->Father()->GetHigherSubElements(sons);
212 /*Limits*/
213 gl_distance_h = gl_distance_hmax*std::pow(this->gradation,this->level_max-gmesh->Element(this->specialelementsindex[i])->Level());
214 if_distance_h = if_distance_hmax*std::pow(this->gradation,this->level_max-gmesh->Element(this->specialelementsindex[i])->Level());
215 d_maxerror = this->deviatoricerror_threshold*this->deviatoricerror_maximum;
216 t_maxerror = this->thicknesserror_threshold*this->thicknesserror_maximum;
217 /*Calculate the distance and error of the group (sons)*/
218 gl_groupdistance=INFINITY;if_groupdistance=INFINITY;deviatoric_grouperror=0;thickness_grouperror=0;
219 for(int s=0;s<sons.size();s++){
220 sid=this->index2sid[sons[s]->Index()];
221 if(sid<0) continue;
222 if(this->groundingline_distance>0) gl_groupdistance=std::min(gl_groupdistance,gl_distance[sid]);
223 if(this->icefront_distance>0) if_groupdistance=std::min(if_groupdistance,if_distance[sid]);
224 if(this->deviatoricerror_threshold>0) deviatoric_grouperror+=deviatoricerror[sid];
225 if(this->thicknesserror_threshold>0) thickness_grouperror+=thicknesserror[sid];
226 }
227 criteria=0;
228 if(this->groundingline_distance>0 && gl_groupdistance<gl_distance_h+DBL_EPSILON) criteria++;
229 if(this->icefront_distance>0 && if_groupdistance<if_distance_h+DBL_EPSILON) criteria++;
230 if(this->deviatoricerror_threshold>0 && deviatoric_grouperror>d_maxerror-DBL_EPSILON) criteria++;
231 if(this->thicknesserror_threshold>0 && thickness_grouperror>t_maxerror-DBL_EPSILON) criteria++;
232 /*Finally, it keeps the father index if it must be refine*/
233 if(criteria) index.push_back(gmesh->Element(this->specialelementsindex[i])->FatherIndex());
234 }
235 /*}}}*/
236
237 /*Now, detele the special elements{{{*/
238 if(this->refinement_type==1) this->DeleteSpecialElements(verbose,gmesh);
239 else this->specialelementsindex.clear();
240 /*}}}*/
241
242 /*Set the mesh and delete previousmesh if refinement type is 0{{{*/
243 if(this->refinement_type==0){
244 delete this->previousmesh;
245 gmesh=this->fathermesh;
246 }
247 /*}}}*/
248
249 /*Unrefinement process: loop over level of refinements{{{*/
250 if(verbose) _printf_("\tunrefinement process...\n");
251 if(verbose) _printf_("\ttotal: ");
252 count=0;
253
254 nelem=gmesh->NElements();//must keep here
255 for(int i=0;i<nelem;i++){
256 if(!gmesh->Element(i)) continue;
257 if(gmesh->Element(i)->MaterialId()!=this->GetElemMaterialID()) continue;
258 if(gmesh->Element(i)->HasSubElement()) continue;
259 if(gmesh->Element(i)->Level()==0) continue;
260 if(!gmesh->Element(i)->Father()) _error_("father is NULL!\n");
261 /*Limits with lag*/
262 gl_distance_h = this->lag*gl_distance_hmax*std::pow(this->gradation,this->level_max-gmesh->Element(i)->Level());
263 if_distance_h = this->lag*if_distance_hmax*std::pow(this->gradation,this->level_max-gmesh->Element(i)->Level());
264 d_maxerror = this->deviatoricerror_groupthreshold*this->deviatoricerror_maximum;
265 t_maxerror = this->thicknesserror_groupthreshold*this->thicknesserror_maximum;
266 /*Get the sons of the father (sibilings)*/
267 sons.clear();
268 gmesh->Element(i)->Father()->GetHigherSubElements(sons);
269 if(sons.size()!=4) continue;//delete just group of 4 elements. This avoids big holes in the mesh
270 /*Find the minimal distance and the error of the group*/
271 gl_groupdistance=INFINITY;if_groupdistance=INFINITY;deviatoric_grouperror=0;thickness_grouperror=0;
272 for(int s=0;s<sons.size();s++){
273 sid=this->index2sid[sons[s]->Index()];
274 /*Verify if this group have solutions*/
275 if(sid<0){gl_groupdistance=INFINITY;if_groupdistance=INFINITY;deviatoric_grouperror=INFINITY;thickness_grouperror=INFINITY;continue;}
276 /*Distance and error of the group*/
277 if(this->groundingline_distance>0) gl_groupdistance=std::min(gl_groupdistance,gl_distance[sid]);
278 if(this->icefront_distance>0) if_groupdistance=std::min(if_groupdistance,if_distance[sid]);
279 if(this->deviatoricerror_threshold>0) deviatoric_grouperror+=deviatoricerror[sid];
280 if(this->thicknesserror_threshold>0) thickness_grouperror+=thicknesserror[sid];
281 }
282 /*Verify the criteria*/
283 criteria=0;
284 if(this->groundingline_distance>0 && gl_groupdistance>gl_distance_h-DBL_EPSILON) criteria++;
285 if(this->icefront_distance>0 && if_groupdistance>if_distance_h-DBL_EPSILON) criteria++;
286 if(this->deviatoricerror_threshold>0 && deviatoric_grouperror<d_maxerror+DBL_EPSILON) criteria++;
287 if(this->thicknesserror_threshold>0 && thickness_grouperror<t_maxerror+DBL_EPSILON) criteria++;
288 /*Now, if the group attends the criteria, unrefine it*/
289 if(criteria==numberofcriteria){
290 gmesh->Element(i)->Father()->ResetSubElements(); count++;
291 for(int s=0;s<sons.size();s++){this->index2sid[sons[s]->Index()]=-1;gmesh->DeleteElement(sons[s],sons[s]->Index());}
292 }
293 }
294 if(verbose) _printf_(""<<count<<"\n");
295 /*Adjust the connectivities before continue*/
296 //gmesh->BuildConnectivity(); this is not necessary
297 /*}}}*/
298
299 /*Refinement process: loop over level of refinements{{{*/
300 if(verbose) _printf_("\trefinement process (level max = "<<this->level_max<<")\n");
301 if(verbose) _printf_("\ttotal: ");
302 count=0;
303 nelem=gmesh->NElements();//must keep here
304 for(int i=0;i<nelem;i++){
305 if(!gmesh->Element(i)) continue;
306 if(gmesh->Element(i)->MaterialId()!=this->GetElemMaterialID()) continue;
307 if(gmesh->Element(i)->HasSubElement()) continue;
308 if(gmesh->Element(i)->Level()==this->level_max) continue;
309 /*Verify if this element has solutions*/
310 sid=this->index2sid[gmesh->Element(i)->Index()];
311 if(sid<0) continue;
312 /*Set the distance for level h*/
313 gl_distance_h = gl_distance_hmax*std::pow(this->gradation,this->level_max-(gmesh->Element(i)->Level()+1));//+1: current element level is <level_max
314 if_distance_h = if_distance_hmax*std::pow(this->gradation,this->level_max-(gmesh->Element(i)->Level()+1));//+1: current element level is <level_max
315 d_maxerror = this->deviatoricerror_threshold*this->deviatoricerror_maximum;
316 t_maxerror = this->thicknesserror_threshold*this->thicknesserror_maximum;
317 /*Verify distance and error of the element, if requested*/
318 criteria=0;
319 if(this->groundingline_distance>0 && gl_distance[sid]<gl_distance_h+DBL_EPSILON) criteria++;
320 if(this->icefront_distance>0 && if_distance[sid]<if_distance_h+DBL_EPSILON) criteria++;
321 if(this->deviatoricerror_threshold>0 && deviatoricerror[sid]>d_maxerror-DBL_EPSILON) criteria++;
322 if(this->thicknesserror_threshold>0 && thicknesserror[sid]>t_maxerror-DBL_EPSILON) criteria++;
323 /*Now, if it attends any criterion, keep the element index to refine in next step*/
324 if(criteria)index.push_back(i);
325 }
326 /*Now, refine the elements*/
327 for(int i=0;i<index.size();i++){
328 if(!gmesh->Element(index[i])) DebugStop();
329 if(!gmesh->Element(index[i])->HasSubElement()){gmesh->Element(index[i])->Divide(sons);count++;}
330 }
331 if(verbose) _printf_(""<<count<<"\n");
332 /*Adjust the connectivities before continue*/
333 //gmesh->BuildConnectivity();//this is not necessary
334 /*}}}*/
335
336 /*Now, apply smoothing and insert special elements to avoid hanging nodes{{{*/
337 this->RefineMeshWithSmoothing(verbose,gmesh);
338 if(this->refinement_type==0) this->previousmesh=this->CreateRefPatternMesh(gmesh);//in this case, gmesh==this->fathermesh
339 gmesh=this->previousmesh;//previous mesh is always refined to avoid hanging nodes
340 this->RefineMeshToAvoidHangingNodes(verbose,gmesh);
341 /*}}}*/
342}
343/*}}}*/
344int AdaptiveMeshRefinement::VerifyRefinementType(TPZGeoEl* geoel,TPZGeoMesh* gmesh){/*{{{*/
345
346 /*
347 * 0 : no refinement
348 * 1 : special refinement (to avoid hanging nodes)
349 * 2 : uniform refinment
350 * */
351 if(!geoel) _error_("geoel is NULL!\n");
352
353 /*Output*/
354 int type=0;
355
356 /*Intermediaries*/
357 TPZVec<TPZGeoEl*> sons;
358
359 /*Loop over neighboors (sides 3, 4 and 5)*/
360 for(int j=3;j<6;j++){
361 if(!gmesh->Element(geoel->NeighbourIndex(j))->HasSubElement()) continue;
362 sons.clear();
363 gmesh->Element(geoel->NeighbourIndex(j))->GetHigherSubElements(sons);
364 if(sons.size()) type++; //if neighbour was refined
365 if(sons.size()>4) type++; //if neighbour's level is > element level+1
366 if(type>1) break;
367 }
368
369 /*Verify and return*/
370 if(type>1) type=2;
371
372 return type;
373}
374/*}}}*/
375void AdaptiveMeshRefinement::RefineMeshWithSmoothing(bool &verbose,TPZGeoMesh* gmesh){/*{{{*/
376
377 /*Intermediaries*/
378 int nelem =-1;
379 int count =-1;
380 int type =-1;
381 int typecount =-1;
382
383 TPZVec<TPZGeoEl*> sons;
384
385 /*Refinement process: loop over level of refinements*/
386 if(verbose) _printf_("\tsmoothing process (level max = "<<this->level_max<<")\n");
387 if(verbose) _printf_("\ttotal: ");
388
389 count=1;
390
391 while(count>0){
392 count=0;
393 nelem=gmesh->NElements();//must keep here
394 for(int i=0;i<nelem;i++){
395 if(!gmesh->Element(i)) continue;
396 if(gmesh->Element(i)->MaterialId()!=this->GetElemMaterialID()) continue;
397 if(gmesh->Element(i)->HasSubElement()) continue;
398 if(gmesh->Element(i)->Level()==this->level_max) continue;
399 /*loop over neighboors (sides 3, 4 and 5). Important: neighbours has the same dimension of the element*/
400 type=this->VerifyRefinementType(gmesh->Element(i),gmesh);
401 if(type<2){
402 typecount=0;
403 for(int j=3;j<6;j++){
404 if(gmesh->Element(gmesh->Element(i)->NeighbourIndex(j))->HasSubElement()) continue;
405 if(gmesh->Element(i)->NeighbourIndex(j)==i) typecount++;//neighbour==this element, element at the border
406 if(this->VerifyRefinementType(gmesh->Element(gmesh->Element(i)->NeighbourIndex(j)),gmesh)==1) typecount++;
407 if(typecount>1 && type==1) type=2;
408 else if(typecount>2 && type==0) type=2;
409 if(type==2) break;
410 }
411 }
412
413 /*refine the element if requested*/
414 if(type==2){gmesh->Element(i)->Divide(sons); count++;}
415 }
416 if(verbose){
417 if(count==0) _printf_(""<<count<<"\n");
418 else _printf_(""<<count<<", ");
419 }
420 /*Adjust the connectivities before continue*/
421 //gmesh->BuildConnectivity();//this is not necessary
422 }
423}
424/*}}}*/
425void AdaptiveMeshRefinement::RefineMeshToAvoidHangingNodes(bool &verbose,TPZGeoMesh* gmesh){/*{{{*/
426
427 /*Now, insert special elements to avoid hanging nodes*/
428 if(verbose) _printf_("\trefining to avoid hanging nodes (total: ");
429
430 /*Intermediaries*/
431 int nelem=-1;
432 int count= 1;
433
434 while(count>0){
435 nelem=gmesh->NElements();//must keep here
436 count=0;
437 for(int i=0;i<nelem;i++){
438 /*Get geometric element and verify if it has already been refined. Geoel may not have been previously refined*/
439 TPZGeoEl * geoel=gmesh->Element(i);
440 if(!geoel) continue;
441 if(geoel->HasSubElement()) continue;
442 if(geoel->MaterialId() != this->GetElemMaterialID()) continue;
443 /*Get the refinement pattern for this element and refine it*/
444 TPZAutoPointer<TPZRefPattern> refp=TPZRefPatternTools::PerfectMatchRefPattern(geoel);
445 if(refp){
446 /*Non-uniform refinement*/
447 TPZVec<TPZGeoEl *> sons;
448 geoel->SetRefPattern(refp);
449 geoel->Divide(sons);
450 count++;
451 /*Keep the index of the special elements*/
452 for(int j=0;j<sons.size();j++) this->specialelementsindex.push_back(sons[j]->Index());
453 }
454 }
455 if(verbose){
456 if(count==0) _printf_(""<<count<<")\n");
457 else _printf_(""<<count<<", ");
458 }
459 //gmesh->BuildConnectivity();//this is not necessary
460 }
461}
462/*}}}*/
463void AdaptiveMeshRefinement::DeleteSpecialElements(bool &verbose,TPZGeoMesh* gmesh){/*{{{*/
464
465 /*Intermediaries*/
466 int count=0;
467
468 if(verbose) _printf_("\tdelete "<<this->specialelementsindex.size()<<" special elements (total: ");
469 for(int i=0;i<this->specialelementsindex.size();i++){
470 if(this->specialelementsindex[i]==-1) continue;
471 if(!gmesh->Element(this->specialelementsindex[i])) continue;
472 if(gmesh->Element(this->specialelementsindex[i])->Father()) gmesh->Element(this->specialelementsindex[i])->Father()->ResetSubElements();
473 gmesh->DeleteElement(gmesh->Element(this->specialelementsindex[i]),this->specialelementsindex[i]);
474 this->index2sid[this->specialelementsindex[i]]=-1;
475 count++;
476 }
477 if(verbose) _printf_(""<<count<<")\n");
478 /*Cleanup*/
479 this->specialelementsindex.clear();
480 /*Adjust connectivities*/
481 //gmesh->BuildConnectivity();//this is not necessary
482}
483/*}}}*/
484void AdaptiveMeshRefinement::GetMesh(TPZGeoMesh* gmesh,int** pdata,double** pxy, int** pelements){/*{{{*/
485
486 /* IMPORTANT! pelements are in Matlab indexing
487 NEOPZ works only in C indexing.
488 This method cleans up and updated the this->sid2index
489 and this->index2sid and fills in it with the new mesh.
490 Avoid to call this method before Refinement Process.*/
491
492 /*Intermediaries */
493 long sid,nodeindex;
494 int nconformelements,nconformvertices;
495 int ntotalvertices = gmesh->NNodes();
496 int* newelements = NULL;
497 int* newdata = NULL;
498 double* newmeshXY = NULL;
499 TPZGeoEl* geoel = NULL;
500 long* vertex_index2sid = xNew<long>(ntotalvertices);
501 this->index2sid.clear(); this->index2sid.resize(gmesh->NElements());
502 this->sid2index.clear();
503
504 /*Fill in the vertex_index2sid vector with non usual index value*/
505 for(int i=0;i<gmesh->NNodes();i++) vertex_index2sid[i]=-1;
506
507 /*Fill in the this->index2sid vector with non usual index value*/
508 for(int i=0;i<gmesh->NElements();i++) this->index2sid[i]=-1;
509
510 /*Get elements without sons and fill in the vertex_index2sid with used vertices (indexes) */
511 sid=0;
512 for(int i=0;i<gmesh->NElements();i++){//over gmesh elements index
513 geoel=gmesh->ElementVec()[i];
514 if(!geoel) continue;
515 if(geoel->HasSubElement()) continue;
516 if(geoel->MaterialId() != this->GetElemMaterialID()) continue;
517 this->sid2index.push_back(geoel->Index());//keep the element index
518 this->index2sid[geoel->Index()]=this->sid2index.size()-1;//keep the element sid
519 for(int j=0;j<this->GetNumberOfNodes();j++){
520 nodeindex=geoel->NodeIndex(j);
521 if(nodeindex<0) _error_("nodeindex is <0\n");
522 if(vertex_index2sid[nodeindex]==-1){
523 vertex_index2sid[nodeindex]=sid;
524 sid++;
525 }
526 }
527 }
528
529 /* Create new mesh structure and fill it */
530 nconformelements = (int)this->sid2index.size();
531 nconformvertices = (int)sid;
532 newelements = xNew<int>(nconformelements*this->GetNumberOfNodes());
533 newmeshXY = xNew<double>(nconformvertices*2);
534 newdata = xNew<int>(2);
535 newdata[0] = nconformvertices;
536 newdata[1] = nconformelements;
537
538 for(int i=0;i<ntotalvertices;i++){//over the TPZNode index (fill in the ISSM vertices coords)
539 sid=vertex_index2sid[i];
540 if(sid==-1) continue;//skip this index (node no used)
541 TPZVec<REAL> coords(3,0.);
542 gmesh->NodeVec()[i].GetCoordinates(coords);
543 newmeshXY[2*sid] = coords[0]; // X
544 newmeshXY[2*sid+1] = coords[1]; // Y
545 }
546
547 for(int i=0;i<this->sid2index.size();i++){//over the sid (fill the ISSM elements)
548 for(int j=0;j<this->GetNumberOfNodes();j++) {
549 geoel = gmesh->ElementVec()[this->sid2index[i]];
550 sid = vertex_index2sid[geoel->NodeIndex(j)];
551 newelements[i*this->GetNumberOfNodes()+j]=(int)sid+1;//C to Matlab indexing
552 }
553 /*Verify the Jacobian determinant. If detJ<0, swap the 2 first postions:
554 a -> b
555 b -> a */
556 double detJ,xa,xb,xc,ya,yb,yc;
557 int a,b,c;
558
559 a=newelements[i*this->GetNumberOfNodes()+0]-1;
560 b=newelements[i*this->GetNumberOfNodes()+1]-1;
561 c=newelements[i*this->GetNumberOfNodes()+2]-1;
562
563 xa=newmeshXY[2*a]; ya=newmeshXY[2*a+1];
564 xb=newmeshXY[2*b]; yb=newmeshXY[2*b+1];
565 xc=newmeshXY[2*c]; yc=newmeshXY[2*c+1];
566
567 detJ=(xb-xa)*(yc-ya)-(xc-xa)*(yb-ya);
568
569 /*verify and swap, if necessary*/
570 if(detJ<0) {
571 newelements[i*this->GetNumberOfNodes()+0]=b+1;//a->b
572 newelements[i*this->GetNumberOfNodes()+1]=a+1;//b->a
573 }
574 }
575
576 /*Setting outputs*/
577 *pdata = newdata;
578 *pxy = newmeshXY;
579 *pelements = newelements;
580
581 /*Cleanup*/
582 xDelete<long>(vertex_index2sid);
583}
584/*}}}*/
585void AdaptiveMeshRefinement::Initialize(){/*{{{*/
586
587 /* IMPORTANT! elements come in Matlab indexing
588 NEOPZ works only in C indexing*/
589
590 if(this->numberofvertices<=0) _error_("Impossible to create initial mesh: nvertices is <= 0!\n");
591 if(this->numberofelements<=0) _error_("Impossible to create initial mesh: nelements is <= 0!\n");
592 if(this->refinement_type!=0 && this->refinement_type!=1) _error_("Impossible to create initial mesh: refinement type is not defined!\n");
593
594 /*Verify and creating initial mesh*/
595 if(!this->x || !this->y || !this->elementslist) _error_("Mesh data structure is NULL!\n");
596 if(this->fathermesh || this->previousmesh) _error_("Initial mesh already exists!\n");
597
598 this->fathermesh = new TPZGeoMesh();
599 this->fathermesh->NodeVec().Resize(this->numberofvertices);
600
601 /*Set the vertices (geometric nodes in NeoPZ context)*/
602 for(int i=0;i<this->numberofvertices;i++){
603 /*x,y,z coords*/
604 TPZManVector<REAL,3> coord(3,0.);
605 coord[0]= this->x[i];
606 coord[1]= this->y[i];
607 coord[2]= 0.;
608 /*Insert in the mesh*/
609 this->fathermesh->NodeVec()[i].SetCoord(coord);
610 this->fathermesh->NodeVec()[i].SetNodeId(i);
611 }
612
613 /*Generate the elements*/
614 long index;
615 const int mat = this->GetElemMaterialID();
616 TPZManVector<long> elem(this->GetNumberOfNodes(),0);
617 this->index2sid.clear(); this->index2sid.resize(this->numberofelements);
618 this->sid2index.clear();
619
620 for(int i=0;i<this->numberofelements;i++){
621 for(int j=0;j<this->GetNumberOfNodes();j++) elem[j]=this->elementslist[i*this->GetNumberOfNodes()+j]-1;//Convert Matlab to C indexing
622 switch(this->GetNumberOfNodes()){
623 case 3: this->fathermesh->CreateGeoElement(ETriangle,elem,mat,index,this->refinement_type); break;
624 default: _error_("mesh not supported yet");
625 }
626 /*Define the element ID*/
627 this->fathermesh->ElementVec()[index]->SetId(i);
628 /*Initialize sid2index and index2sid*/
629 this->sid2index.push_back((int)index);
630 this->index2sid[(int)index]=this->sid2index.size()-1;//keep the element sid
631 }
632 /*Build element and node connectivities*/
633 this->fathermesh->BuildConnectivity();
634 /*Set previous mesh*/
635 if(this->refinement_type==1) this->previousmesh=new TPZGeoMesh(*this->fathermesh);
636 else this->previousmesh=this->CreateRefPatternMesh(this->fathermesh);
637}
638/*}}}*/
639TPZGeoMesh* AdaptiveMeshRefinement::CreateRefPatternMesh(TPZGeoMesh* gmesh){/*{{{*/
640
641 TPZGeoMesh *newgmesh = new TPZGeoMesh();
642 newgmesh->CleanUp();
643
644 int nnodes = gmesh->NNodes();
645 int nelem = gmesh->NElements();
646 int mat = this->GetElemMaterialID();;
647 int reftype = 1;
648 long index;
649
650 //nodes
651 newgmesh->NodeVec().Resize(nnodes);
652 for(int i=0;i<nnodes;i++) newgmesh->NodeVec()[i] = gmesh->NodeVec()[i];
653
654 //elements
655 for(int i=0;i<nelem;i++){
656 TPZGeoEl * geoel = gmesh->Element(i);
657
658 if(!geoel){
659 index=newgmesh->ElementVec().AllocateNewElement();
660 newgmesh->ElementVec()[index] = NULL;
661 continue;
662 }
663
664 TPZManVector<long> elem(3,0);
665 for(int j=0;j<3;j++) elem[j] = geoel->NodeIndex(j);
666
667 newgmesh->CreateGeoElement(ETriangle,elem,mat,index,reftype);
668 newgmesh->ElementVec()[index]->SetId(geoel->Id());
669
670 TPZGeoElRefPattern<TPZGeoTriangle>* newgeoel = dynamic_cast<TPZGeoElRefPattern<TPZGeoTriangle>*>(newgmesh->ElementVec()[index]);
671
672 //old neighbourhood
673 const int nsides = TPZGeoTriangle::NSides;
674 TPZVec< std::vector<TPZGeoElSide> > neighbourhood(nsides);
675 TPZVec<long> NodesSequence(0);
676 for(int s = 0; s < nsides; s++){
677 neighbourhood[s].resize(0);
678 TPZGeoElSide mySide(geoel,s);
679 TPZGeoElSide neighS = mySide.Neighbour();
680 if(mySide.Dimension() == 0){
681 long oldSz = NodesSequence.NElements();
682 NodesSequence.resize(oldSz+1);
683 NodesSequence[oldSz] = geoel->NodeIndex(s);
684 }
685 while(mySide != neighS){
686 neighbourhood[s].push_back(neighS);
687 neighS = neighS.Neighbour();
688 }
689 }
690
691 //inserting in new element
692 for(int s = 0; s < nsides; s++){
693 TPZGeoEl * tempEl = newgeoel;
694 TPZGeoElSide tempSide(newgeoel,s);
695 int byside = s;
696 for(unsigned long n = 0; n < neighbourhood[s].size(); n++){
697 TPZGeoElSide neighS = neighbourhood[s][n];
698 tempEl->SetNeighbour(byside, neighS);
699 tempEl = neighS.Element();
700 byside = neighS.Side();
701 }
702 tempEl->SetNeighbour(byside, tempSide);
703 }
704
705 long fatherindex = geoel->FatherIndex();
706 if(fatherindex>-1) newgeoel->SetFather(fatherindex);
707
708 if(!geoel->HasSubElement()) continue;
709
710 int nsons = geoel->NSubElements();
711
712 TPZAutoPointer<TPZRefPattern> ref = gRefDBase.GetUniformRefPattern(ETriangle);
713 newgeoel->SetRefPattern(ref);
714
715 for(int j=0;j<nsons;j++){
716 TPZGeoEl* son = geoel->SubElement(j);
717 if(!son){
718 DebugStop();
719 }
720 newgeoel->SetSubElement(j,son);
721 }
722 }
723
724 /*Now, build connectivities*/
725 newgmesh->BuildConnectivity();
726
727 return newgmesh;
728}
729/*}}}*/
730void AdaptiveMeshRefinement::CheckMesh(int** pdata,double** pxy,int** pelements){/*{{{*/
731
732 /*Basic verification*/
733 if(!pdata) _error_("Impossible to continue: pdata is NULL!\n");
734 if(pdata[0]<=0) _error_("Impossible to continue: nvertices <=0!\n");
735 if(pdata[1]<=0) _error_("Impossible to continue: nelements <=0!\n");
736 if(!pxy) _error_("Impossible to continue: pxy is NULL!\n");
737 if(!pelements) _error_("Impossible to continue: pelements is NULL!\n");
738}
739/*}}}*/
740void AdaptiveMeshRefinement::PrintGMeshVTK(TPZGeoMesh* gmesh,std::ofstream &file,bool matColor){/*{{{*/
741
742 file.clear();
743 long nelements = gmesh->NElements();
744 TPZGeoEl *gel;
745 std::stringstream node, connectivity, type, material;
746
747 //Header
748 file << "# vtk DataFile Version 3.0" << std::endl;
749 file << "TPZGeoMesh VTK Visualization" << std::endl;
750 file << "ASCII" << std::endl << std::endl;
751 file << "DATASET UNSTRUCTURED_GRID" << std::endl;
752 file << "POINTS ";
753
754 long actualNode = -1, size = 0, nVALIDelements = 0;
755 for(long el = 0; el < nelements; el++){
756 gel = gmesh->ElementVec()[el];
757 if(!gel )//|| (gel->Type() == EOned && !gel->IsLinearMapping()))//Exclude Arc3D and Ellipse3D
758 {
759 continue;
760 }
761 if(gel->HasSubElement())
762 {
763 continue;
764 }
765 MElementType elt = gel->Type();
766 int elNnodes = MElementType_NNodes(elt);
767
768 size += (1+elNnodes);
769 connectivity << elNnodes;
770
771 for(int t = 0; t < elNnodes; t++)
772 {
773 for(int c = 0; c < 3; c++)
774 {
775 double coord = gmesh->NodeVec()[gel->NodeIndex(t)].Coord(c);
776 node << coord << " ";
777 }
778 node << std::endl;
779
780 actualNode++;
781 connectivity << " " << actualNode;
782 }
783 connectivity << std::endl;
784
785 int elType = this->GetVTK_ElType(gel);
786 type << elType << std::endl;
787
788 if(matColor == true)
789 {
790 material << gel->MaterialId() << std::endl;
791 }
792 else
793 {
794 material << gel->Index() << std::endl;
795 }
796
797 nVALIDelements++;
798 }
799 node << std::endl;
800 actualNode++;
801 file << actualNode << " float" << std::endl << node.str();
802
803 file << "CELLS " << nVALIDelements << " ";
804
805 file << size << std::endl;
806 file << connectivity.str() << std::endl;
807
808 file << "CELL_TYPES " << nVALIDelements << std::endl;
809 file << type.str() << std::endl;
810
811 file << "CELL_DATA" << " " << nVALIDelements << std::endl;
812 file << "FIELD FieldData 1" << std::endl;
813 if(matColor == true)
814 {
815 file << "material 1 " << nVALIDelements << " int" << std::endl;
816 }
817 else
818 {
819 file << "ElementIndex 1 " << nVALIDelements << " int" << std::endl;
820 }
821 file << material.str();
822 file.close();
823}
824/*}}}*/
825int AdaptiveMeshRefinement::GetVTK_ElType(TPZGeoEl * gel){/*{{{*/
826
827 MElementType pzElType = gel->Type();
828
829 int elType = -1;
830 switch (pzElType)
831 {
832 case(EPoint):
833 {
834 elType = 1;
835 break;
836 }
837 case(EOned):
838 {
839 elType = 3;
840 break;
841 }
842 case (ETriangle):
843 {
844 elType = 5;
845 break;
846 }
847 case (EQuadrilateral):
848 {
849 elType = 9;
850 break;
851 }
852 case (ETetraedro):
853 {
854 elType = 10;
855 break;
856 }
857 case (EPiramide):
858 {
859 elType = 14;
860 break;
861 }
862 case (EPrisma):
863 {
864 elType = 13;
865 break;
866 }
867 case (ECube):
868 {
869 elType = 12;
870 break;
871 }
872 default:
873 {
874 std::cout << "Element type not found on " << __PRETTY_FUNCTION__ << std::endl;
875 DebugStop();
876 break;
877 }
878 }
879 if(elType == -1)
880 {
881 std::cout << "Element type not found on " << __PRETTY_FUNCTION__ << std::endl;
882 std::cout << "MIGHT BE CURVED ELEMENT (quadratic or quarter point)" << std::endl;
883 DebugStop();
884 }
885
886 return elType;
887}
888/*}}}*/
889void AdaptiveMeshRefinement::ReadMesh(){/*{{{*/
890
891 std::string fathermeshfile = "/home/santos/issm_fathermesh.txt";
892 std::string previousmeshfile = "/home/santos/issm_previousmesh.txt";
893 std::string amrfile = "/home/santos/issm_amr.txt";
894 std::ifstream amrifstream(amrfile.c_str());
895 int size,value;
896
897 TPZPersistenceManager::OpenRead(fathermeshfile);
898 this->fathermesh = dynamic_cast<TPZGeoMesh *>(TPZPersistenceManager::ReadFromFile());
899 TPZPersistenceManager::CloseRead();
900
901 TPZPersistenceManager::OpenRead(previousmeshfile);
902 this->previousmesh = dynamic_cast<TPZGeoMesh *>(TPZPersistenceManager::ReadFromFile());
903 TPZPersistenceManager::CloseRead();
904
905 if(!amrifstream.is_open()) _error_("amr ifstream is not open!");
906 amrifstream.seekg(0);
907
908 this->sid2index.clear();
909 this->index2sid.clear();
910 this->specialelementsindex.clear();
911
912 amrifstream>>size;
913 for(int i=0;i<size;i++){
914 amrifstream>>value;
915 this->sid2index.push_back(value);
916 }
917
918 amrifstream>>size;
919 for(int i=0;i<size;i++){
920 amrifstream>>value;
921 this->index2sid.push_back(value);
922 }
923
924 amrifstream>>size;
925 for(int i=0;i<size;i++){
926 amrifstream>>value;
927 this->specialelementsindex.push_back(value);
928 }
929}
930/*}}}*/
931void AdaptiveMeshRefinement::WriteMesh(){/*{{{*/
932
933 std::string fathermeshfile = "/home/santos/issm_fathermesh.txt";
934 std::string previousmeshfile = "/home/santos/issm_previousmesh.txt";
935 std::string amrfile = "/home/santos/issm_amr.txt";
936 std::ofstream amrofstream(amrfile.c_str());
937
938 if(this->fathermesh){
939 TPZPersistenceManager::OpenWrite(fathermeshfile);
940 TPZPersistenceManager::WriteToFile(this->fathermesh);
941 TPZPersistenceManager::CloseWrite();
942 }
943
944 if(this->previousmesh){
945 TPZPersistenceManager::OpenWrite(previousmeshfile);
946 TPZPersistenceManager::WriteToFile(this->previousmesh);
947 TPZPersistenceManager::CloseWrite();
948 }
949
950 if(this->sid2index.size()>0){
951 amrofstream << this->sid2index.size() << std::endl;
952 for(int i=0;i<this->sid2index.size();i++) {
953 amrofstream << this->sid2index[i] << std::endl;
954 }
955 }
956
957 if(this->index2sid.size()>0){
958 amrofstream << this->index2sid.size() << std::endl;
959 for(int i=0;i<this->index2sid.size();i++) {
960 amrofstream << this->index2sid[i] << std::endl;
961 }
962 }
963
964 if(this->specialelementsindex.size()){
965 amrofstream << this->specialelementsindex.size() << std::endl;
966 for(int i=0;i<this->specialelementsindex.size();i++) {
967 amrofstream << this->specialelementsindex[i] << std::endl;
968 }
969 }
970 amrofstream.flush();
971 amrofstream.close();
972}
973/*}}}*/
Note: See TracBrowser for help on using the repository browser.