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

Last change on this file since 22758 was 22758, checked in by Mathieu Morlighem, 7 years ago

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