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

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

merged trunk-jpl and trunk for revision 23187
File size: 30.5 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** pdatalist,double** pxylist,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,pdatalist,pxylist,pelementslist);
126
127 /*Verify the new geometry*/
128 this->CheckMesh(pdatalist,pxylist,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** pdata,double** pxy, 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 int* newdata = NULL;
465 double* newmeshXY = 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 /* Create new mesh structure and fill it */
497 nconformelements = (int)this->sid2index.size();
498 nconformvertices = (int)sid;
499 newelements = xNew<int>(nconformelements*this->GetNumberOfNodes());
500 newmeshXY = xNew<double>(nconformvertices*2);
501 newdata = xNew<int>(2);
502 newdata[0] = nconformvertices;
503 newdata[1] = nconformelements;
504
505 for(int i=0;i<ntotalvertices;i++){//over the TPZNode index (fill in the ISSM vertices coords)
506 sid=vertex_index2sid[i];
507 if(sid==-1) continue;//skip this index (node no used)
508 TPZVec<REAL> coords(3,0.);
509 gmesh->NodeVec()[i].GetCoordinates(coords);
510 newmeshXY[2*sid] = coords[0]; // X
511 newmeshXY[2*sid+1] = coords[1]; // Y
512 }
513
514 for(int i=0;i<this->sid2index.size();i++){//over the sid (fill the ISSM elements)
515 for(int j=0;j<this->GetNumberOfNodes();j++) {
516 geoel = gmesh->ElementVec()[this->sid2index[i]];
517 sid = vertex_index2sid[geoel->NodeIndex(j)];
518 newelements[i*this->GetNumberOfNodes()+j]=(int)sid+1;//C to Matlab indexing
519 }
520 /*Verify the Jacobian determinant. If detJ<0, swap the 2 first postions:
521 a -> b
522 b -> a */
523 double detJ,xa,xb,xc,ya,yb,yc;
524 int a,b,c;
525
526 a=newelements[i*this->GetNumberOfNodes()+0]-1;
527 b=newelements[i*this->GetNumberOfNodes()+1]-1;
528 c=newelements[i*this->GetNumberOfNodes()+2]-1;
529
530 xa=newmeshXY[2*a]; ya=newmeshXY[2*a+1];
531 xb=newmeshXY[2*b]; yb=newmeshXY[2*b+1];
532 xc=newmeshXY[2*c]; yc=newmeshXY[2*c+1];
533
534 detJ=(xb-xa)*(yc-ya)-(xc-xa)*(yb-ya);
535
536 /*verify and swap, if necessary*/
537 if(detJ<0) {
538 newelements[i*this->GetNumberOfNodes()+0]=b+1;//a->b
539 newelements[i*this->GetNumberOfNodes()+1]=a+1;//b->a
540 }
541 }
542
543 /*Setting outputs*/
544 *pdata = newdata;
545 *pxy = newmeshXY;
546 *pelements = newelements;
547
548 /*Cleanup*/
549 xDelete<long>(vertex_index2sid);
550}
551/*}}}*/
552void AdaptiveMeshRefinement::CreateInitialMesh(int &nvertices,int &nelements,double* x,double* y,int* elements){/*{{{*/
553
554 /* IMPORTANT! elements come in Matlab indexing
555 NEOPZ works only in C indexing*/
556
557 if(nvertices<=0) _error_("Impossible to create initial mesh: nvertices is <= 0!\n");
558 if(nelements<=0) _error_("Impossible to create initial mesh: nelements is <= 0!\n");
559 if(this->refinement_type!=0 && this->refinement_type!=1) _error_("Impossible to create initial mesh: refinement type is not defined!\n");
560
561 /*Verify and creating initial mesh*/
562 if(this->fathermesh || this->previousmesh) _error_("Initial mesh already exists!");
563
564 this->fathermesh = new TPZGeoMesh();
565 this->fathermesh->NodeVec().Resize(nvertices);
566
567 /*Set the vertices (geometric nodes in NeoPZ context)*/
568 for(int i=0;i<nvertices;i++){
569 /*x,y,z coords*/
570 TPZManVector<REAL,3> coord(3,0.);
571 coord[0]= x[i];
572 coord[1]= y[i];
573 coord[2]= 0.;
574 /*Insert in the mesh*/
575 this->fathermesh->NodeVec()[i].SetCoord(coord);
576 this->fathermesh->NodeVec()[i].SetNodeId(i);
577 }
578
579 /*Generate the elements*/
580 long index;
581 const int mat = this->GetElemMaterialID();
582 TPZManVector<long> elem(this->GetNumberOfNodes(),0);
583 this->index2sid.clear(); this->index2sid.resize(nelements);
584 this->sid2index.clear();
585
586 for(int i=0;i<nelements;i++){
587 for(int j=0;j<this->GetNumberOfNodes();j++) elem[j]=elements[i*this->GetNumberOfNodes()+j]-1;//Convert Matlab to C indexing
588 switch(this->GetNumberOfNodes()){
589 case 3: this->fathermesh->CreateGeoElement(ETriangle,elem,mat,index,this->refinement_type); break;
590 default: _error_("mesh not supported yet");
591 }
592 /*Define the element ID*/
593 this->fathermesh->ElementVec()[index]->SetId(i);
594 /*Initialize sid2index and index2sid*/
595 this->sid2index.push_back((int)index);
596 this->index2sid[(int)index]=this->sid2index.size()-1;//keep the element sid
597 }
598 /*Build element and node connectivities*/
599 this->fathermesh->BuildConnectivity();
600 /*Set previous mesh*/
601 if(this->refinement_type==1) this->previousmesh=new TPZGeoMesh(*this->fathermesh);
602 else this->previousmesh=this->CreateRefPatternMesh(this->fathermesh);
603}
604/*}}}*/
605TPZGeoMesh* AdaptiveMeshRefinement::CreateRefPatternMesh(TPZGeoMesh* gmesh){/*{{{*/
606
607 TPZGeoMesh *newgmesh = new TPZGeoMesh();
608 newgmesh->CleanUp();
609
610 int nnodes = gmesh->NNodes();
611 int nelem = gmesh->NElements();
612 int mat = this->GetElemMaterialID();;
613 int reftype = 1;
614 long index;
615
616 //nodes
617 newgmesh->NodeVec().Resize(nnodes);
618 for(int i=0;i<nnodes;i++) newgmesh->NodeVec()[i] = gmesh->NodeVec()[i];
619
620 //elements
621 for(int i=0;i<nelem;i++){
622 TPZGeoEl * geoel = gmesh->Element(i);
623
624 if(!geoel){
625 index=newgmesh->ElementVec().AllocateNewElement();
626 newgmesh->ElementVec()[index] = NULL;
627 continue;
628 }
629
630 TPZManVector<long> elem(3,0);
631 for(int j=0;j<3;j++) elem[j] = geoel->NodeIndex(j);
632
633 newgmesh->CreateGeoElement(ETriangle,elem,mat,index,reftype);
634 newgmesh->ElementVec()[index]->SetId(geoel->Id());
635
636 TPZGeoElRefPattern<TPZGeoTriangle>* newgeoel = dynamic_cast<TPZGeoElRefPattern<TPZGeoTriangle>*>(newgmesh->ElementVec()[index]);
637
638 //old neighbourhood
639 const int nsides = TPZGeoTriangle::NSides;
640 TPZVec< std::vector<TPZGeoElSide> > neighbourhood(nsides);
641 TPZVec<long> NodesSequence(0);
642 for(int s = 0; s < nsides; s++){
643 neighbourhood[s].resize(0);
644 TPZGeoElSide mySide(geoel,s);
645 TPZGeoElSide neighS = mySide.Neighbour();
646 if(mySide.Dimension() == 0){
647 long oldSz = NodesSequence.NElements();
648 NodesSequence.resize(oldSz+1);
649 NodesSequence[oldSz] = geoel->NodeIndex(s);
650 }
651 while(mySide != neighS){
652 neighbourhood[s].push_back(neighS);
653 neighS = neighS.Neighbour();
654 }
655 }
656
657 //inserting in new element
658 for(int s = 0; s < nsides; s++){
659 TPZGeoEl * tempEl = newgeoel;
660 TPZGeoElSide tempSide(newgeoel,s);
661 int byside = s;
662 for(unsigned long n = 0; n < neighbourhood[s].size(); n++){
663 TPZGeoElSide neighS = neighbourhood[s][n];
664 tempEl->SetNeighbour(byside, neighS);
665 tempEl = neighS.Element();
666 byside = neighS.Side();
667 }
668 tempEl->SetNeighbour(byside, tempSide);
669 }
670
671 long fatherindex = geoel->FatherIndex();
672 if(fatherindex>-1) newgeoel->SetFather(fatherindex);
673
674 if(!geoel->HasSubElement()) continue;
675
676 int nsons = geoel->NSubElements();
677
678 TPZAutoPointer<TPZRefPattern> ref = gRefDBase.GetUniformRefPattern(ETriangle);
679 newgeoel->SetRefPattern(ref);
680
681 for(int j=0;j<nsons;j++){
682 TPZGeoEl* son = geoel->SubElement(j);
683 if(!son){
684 DebugStop();
685 }
686 newgeoel->SetSubElement(j,son);
687 }
688 }
689
690 /*Now, build connectivities*/
691 newgmesh->BuildConnectivity();
692
693 return newgmesh;
694}
695/*}}}*/
696void AdaptiveMeshRefinement::CheckMesh(int** pdata,double** pxy,int** pelements){/*{{{*/
697
698 /*Basic verification*/
699 if(!pdata) _error_("Impossible to continue: pdata is NULL!\n");
700 if(pdata[0]<=0) _error_("Impossible to continue: nvertices <=0!\n");
701 if(pdata[1]<=0) _error_("Impossible to continue: nelements <=0!\n");
702 if(!pxy) _error_("Impossible to continue: pxy is NULL!\n");
703 if(!pelements) _error_("Impossible to continue: pelements is NULL!\n");
704}
705/*}}}*/
706void AdaptiveMeshRefinement::PrintGMeshVTK(TPZGeoMesh* gmesh,std::ofstream &file,bool matColor){/*{{{*/
707
708 file.clear();
709 long nelements = gmesh->NElements();
710 TPZGeoEl *gel;
711 std::stringstream node, connectivity, type, material;
712
713 //Header
714 file << "# vtk DataFile Version 3.0" << std::endl;
715 file << "TPZGeoMesh VTK Visualization" << std::endl;
716 file << "ASCII" << std::endl << std::endl;
717 file << "DATASET UNSTRUCTURED_GRID" << std::endl;
718 file << "POINTS ";
719
720 long actualNode = -1, size = 0, nVALIDelements = 0;
721 for(long el = 0; el < nelements; el++){
722 gel = gmesh->ElementVec()[el];
723 if(!gel )//|| (gel->Type() == EOned && !gel->IsLinearMapping()))//Exclude Arc3D and Ellipse3D
724 {
725 continue;
726 }
727 if(gel->HasSubElement())
728 {
729 continue;
730 }
731 MElementType elt = gel->Type();
732 int elNnodes = MElementType_NNodes(elt);
733
734 size += (1+elNnodes);
735 connectivity << elNnodes;
736
737 for(int t = 0; t < elNnodes; t++)
738 {
739 for(int c = 0; c < 3; c++)
740 {
741 double coord = gmesh->NodeVec()[gel->NodeIndex(t)].Coord(c);
742 node << coord << " ";
743 }
744 node << std::endl;
745
746 actualNode++;
747 connectivity << " " << actualNode;
748 }
749 connectivity << std::endl;
750
751 int elType = this->GetVTK_ElType(gel);
752 type << elType << std::endl;
753
754 if(matColor == true)
755 {
756 material << gel->MaterialId() << std::endl;
757 }
758 else
759 {
760 material << gel->Index() << std::endl;
761 }
762
763 nVALIDelements++;
764 }
765 node << std::endl;
766 actualNode++;
767 file << actualNode << " float" << std::endl << node.str();
768
769 file << "CELLS " << nVALIDelements << " ";
770
771 file << size << std::endl;
772 file << connectivity.str() << std::endl;
773
774 file << "CELL_TYPES " << nVALIDelements << std::endl;
775 file << type.str() << std::endl;
776
777 file << "CELL_DATA" << " " << nVALIDelements << std::endl;
778 file << "FIELD FieldData 1" << std::endl;
779 if(matColor == true)
780 {
781 file << "material 1 " << nVALIDelements << " int" << std::endl;
782 }
783 else
784 {
785 file << "ElementIndex 1 " << nVALIDelements << " int" << std::endl;
786 }
787 file << material.str();
788 file.close();
789}
790/*}}}*/
791int AdaptiveMeshRefinement::GetVTK_ElType(TPZGeoEl * gel){/*{{{*/
792
793 MElementType pzElType = gel->Type();
794
795 int elType = -1;
796 switch (pzElType)
797 {
798 case(EPoint):
799 {
800 elType = 1;
801 break;
802 }
803 case(EOned):
804 {
805 elType = 3;
806 break;
807 }
808 case (ETriangle):
809 {
810 elType = 5;
811 break;
812 }
813 case (EQuadrilateral):
814 {
815 elType = 9;
816 break;
817 }
818 case (ETetraedro):
819 {
820 elType = 10;
821 break;
822 }
823 case (EPiramide):
824 {
825 elType = 14;
826 break;
827 }
828 case (EPrisma):
829 {
830 elType = 13;
831 break;
832 }
833 case (ECube):
834 {
835 elType = 12;
836 break;
837 }
838 default:
839 {
840 std::cout << "Element type not found on " << __PRETTY_FUNCTION__ << std::endl;
841 DebugStop();
842 break;
843 }
844 }
845 if(elType == -1)
846 {
847 std::cout << "Element type not found on " << __PRETTY_FUNCTION__ << std::endl;
848 std::cout << "MIGHT BE CURVED ELEMENT (quadratic or quarter point)" << std::endl;
849 DebugStop();
850 }
851
852 return elType;
853}
854/*}}}*/
Note: See TracBrowser for help on using the repository browser.