source: issm/trunk/src/c/classes/AmrNeopz.cpp@ 26744

Last change on this file since 26744 was 26744, checked in by Mathieu Morlighem, 3 years ago

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