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

Last change on this file since 21729 was 21729, checked in by Mathieu Morlighem, 8 years ago

merged trunk-jpl and trunk for revision 21727
File size: 23.6 KB
Line 
1/*!\file AdaptiveMeshrefinement.cpp
2 * \brief: implementation of the adaptive mesh refinement tool based on NeoPZ library: github.com/labmec/neopz
3 */
4
5#ifdef HAVE_CONFIG_H
6 #include <config.h>
7#else
8#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
9#endif
10
11#include "./AdaptiveMeshRefinement.h"
12
13/*Constructor, copy, clean up and destructor*/
14AdaptiveMeshRefinement::AdaptiveMeshRefinement(){/*{{{*/
15 this->Initialize();
16}
17/*}}}*/
18AdaptiveMeshRefinement::AdaptiveMeshRefinement(const AdaptiveMeshRefinement &cp){/*{{{*/
19 this->Initialize();
20 this->operator =(cp);
21}
22/*}}}*/
23AdaptiveMeshRefinement & AdaptiveMeshRefinement::operator =(const AdaptiveMeshRefinement &cp){/*{{{*/
24
25 /*Clean all attributes*/
26 this->CleanUp();
27
28 /*Copy all data*/
29 this->fathermesh = new TPZGeoMesh(*cp.fathermesh);
30 this->previousmesh = new TPZGeoMesh(*cp.previousmesh);
31 this->levelmax = cp.levelmax;
32 this->elementswidth = cp.elementswidth;
33 this->regionlevel1 = cp.regionlevel1;
34 this->regionlevelmax = cp.regionlevelmax;
35 return *this;
36
37}
38/*}}}*/
39AdaptiveMeshRefinement::~AdaptiveMeshRefinement(){/*{{{*/
40
41 bool ismismip = true;
42 if(ismismip){//itapopo
43 TPZFileStream fstr;
44 std::stringstream ss;
45
46 ss << this->levelmax;
47 std::string AMRfile = "/home/santos/Misomip2/L" + ss.str() + "_tsai/amr.txt";
48
49 fstr.OpenWrite(AMRfile.c_str());
50 int withclassid = 1;
51 this->Write(fstr,withclassid);
52 }
53 this->CleanUp();
54 gRefDBase.clear();
55}
56/*}}}*/
57void AdaptiveMeshRefinement::CleanUp(){/*{{{*/
58
59 /*Verify and delete all data*/
60 if(this->fathermesh) delete this->fathermesh;
61 if(this->previousmesh) delete this->previousmesh;
62 this->levelmax = -1;
63 this->elementswidth = -1;
64 this->regionlevel1 = -1;
65 this->regionlevelmax = -1;
66}
67/*}}}*/
68void AdaptiveMeshRefinement::Initialize(){/*{{{*/
69
70 /*Set pointers to NULL*/
71 this->fathermesh = NULL;
72 this->previousmesh = NULL;
73 this->levelmax = -1;
74 this->elementswidth = -1;
75 this->regionlevel1 = -1;
76 this->regionlevelmax = -1;
77}
78/*}}}*/
79int AdaptiveMeshRefinement::ClassId() const{/*{{{*/
80 return 13829430; //Antartic area with ice shelves (km^2)
81}
82/*}}}*/
83void AdaptiveMeshRefinement::Read(TPZStream &buf, void *context){/*{{{*/
84
85 try
86 {
87 /* Read the id context*/
88 TPZSaveable::Read(buf,context);
89
90 /* Read class id*/
91 int classid;
92 buf.Read(&classid,1);
93
94 /* Verify the class id*/
95 if (classid != this->ClassId() )
96 {
97 std::cout << "Error in restoring AdaptiveMeshRefinement!\n";
98 std::cout.flush();
99 DebugStop();
100 }
101
102 /* Read simple attributes */
103 buf.Read(&this->levelmax,1);
104 buf.Read(&this->elementswidth,1);
105 buf.Read(&this->regionlevel1,1);
106 buf.Read(&this->regionlevelmax,1);
107
108 /* Read geometric mesh*/
109 TPZSaveable *sv1 = TPZSaveable::Restore(buf,0);
110 this->fathermesh = dynamic_cast<TPZGeoMesh*>(sv1);
111
112 TPZSaveable *sv2 = TPZSaveable::Restore(buf,0);
113 this->previousmesh = dynamic_cast<TPZGeoMesh*>(sv2);
114 }
115 catch(const std::exception& e)
116 {
117 std::cout << "Exception catched! " << e.what() << std::endl;
118 std::cout.flush();
119 DebugStop();
120 }
121}
122/*}}}*/
123template class TPZRestoreClass<AdaptiveMeshRefinement,13829430>;/*{{{*/
124/*}}}*/
125void AdaptiveMeshRefinement::Write(TPZStream &buf, int withclassid){/*{{{*/
126
127 try
128 {
129 /* Write context (this class) class ID*/
130 TPZSaveable::Write(buf,withclassid);
131
132 /* Write this class id*/
133 int classid = ClassId();
134 buf.Write(&classid,1);
135
136 /* Write simple attributes */
137 buf.Write(&this->levelmax,1);
138 buf.Write(&this->elementswidth,1);
139 buf.Write(&this->regionlevel1,1);
140 buf.Write(&this->regionlevelmax,1);
141
142 /* Write the geometric mesh*/
143 this->fathermesh->Write(buf, this->ClassId());
144 this->previousmesh->Write(buf, this->ClassId());
145 }
146 catch(const std::exception& e)
147 {
148 std::cout << "Exception catched! " << e.what() << std::endl;
149 std::cout.flush();
150 DebugStop();
151 }
152}
153/*}}}*/
154
155/*Mesh refinement methods*/
156#include "TPZVTKGeoMesh.h" //itapopo
157#include "../shared/shared.h" //itapopo
158void AdaptiveMeshRefinement::ExecuteRefinement(int &type_process,double *vx, double *vy, double *masklevelset, int &nvertices, int &nelements, int &nsegments, double** px, double** py, double** pz, int** pelements, int** psegments){/*{{{*/
159
160 /*IMPORTANT! pelements (and psegments) are in Matlab indexing*/
161 /*NEOPZ works only in C indexing*/
162
163 _assert_(this->fathermesh);
164 _assert_(this->previousmesh);
165
166 /*Calculate the position of the grounding line using previous mesh*/
167 std::vector<TPZVec<REAL> > GLvec;
168 this->CalcGroundingLinePosition(masklevelset, GLvec);
169
170 // std::ofstream file1("/home/santos/mesh0.vtk");
171 // TPZVTKGeoMesh::PrintGMeshVTK(this->fathermesh,file1 );
172
173 /*run refinement or unrefinement process*/
174 TPZGeoMesh *newmesh;
175 switch (type_process) {
176 case 0: newmesh = this->previousmesh; break; // refine previous mesh
177 case 1: newmesh = new TPZGeoMesh(*this->fathermesh); break; // refine mesh 0 (unrefine process)
178 default: DebugStop(); break;//itapopo verificar se irá usar _assert_
179 }
180
181 this->RefinementProcess(newmesh,GLvec);
182
183 //std::ofstream file2("/home/santos/mesh1.vtk");
184 //TPZVTKGeoMesh::PrintGMeshVTK(this->previousmesh,file2 );
185
186 /*Set new mesh pointer. Previous mesh just have uniform elements*/
187 if(type_process==1){
188 if(this->previousmesh) delete this->previousmesh;
189 this->previousmesh = newmesh;
190 }
191
192 /*Refine elements to avoid hanging nodes*/
193 //TPZGeoMesh *nohangingnodesmesh = new TPZGeoMesh(*newmesh);//itapopo testando, este era o original
194 TPZGeoMesh *nohangingnodesmesh = this->CreateRefPatternMesh(newmesh);//itapopo testando, este eh novo metodo
195
196 //std::ofstream file3("/home/santos/mesh2.vtk");
197 //TPZVTKGeoMesh::PrintGMeshVTK(this->previousmesh,file3);
198
199 this->RefineMeshToAvoidHangingNodes(nohangingnodesmesh);
200
201 //std::ofstream file4("/home/santos/mesh3.vtk");
202 //TPZVTKGeoMesh::PrintGMeshVTK(nohangingnodesmesh,file4);
203
204 /*Get new geometric mesh in ISSM data structure*/
205 this->GetMesh(nohangingnodesmesh,nvertices,nelements,nsegments,px,py,pz,pelements,psegments);
206
207 /*Verify the new geometry*/
208 this->CheckMesh(nvertices,nelements,nsegments,this->elementswidth,px,py,pz,pelements,psegments);
209
210 _printf_("\trefinement process done!\n\n");
211
212 delete nohangingnodesmesh;
213}
214/*}}}*/
215void AdaptiveMeshRefinement::RefinementProcess(TPZGeoMesh *gmesh,std::vector<TPZVec<REAL> > &GLvec){/*{{{*/
216
217 /*Refine mesh levelmax times*/
218 _printf_("\n\trefinement process (level max = " << this->levelmax << ")\n");
219 _printf_("\tprogress: ");
220 for(int hlevel=1;hlevel<=this->levelmax;hlevel++){
221
222 /*Set elements to be refined using some criteria*/
223 std::vector<int> ElemVec; //elements without children
224 this->SetElementsToRefine(gmesh,GLvec,hlevel,ElemVec);
225
226 /*Refine the mesh*/
227 this->RefineMesh(gmesh, ElemVec);
228
229 _printf_("* ");
230 }
231 _printf_("\n");
232}
233/*}}}*/
234void AdaptiveMeshRefinement::RefineMesh(TPZGeoMesh *gmesh, std::vector<int> &ElemVec){/*{{{*/
235
236 /*Refine elements in ElemVec: uniform pattern refinement*/
237 for(int i = 0; i < ElemVec.size(); i++){
238
239 /*Get geometric element and verify if it has already been refined*/
240 int index = ElemVec[i];
241 TPZGeoEl * geoel = gmesh->Element(index);
242 if(geoel->HasSubElement()) DebugStop(); //itapopo _assert_(!geoel->HasSubElement());
243 if(geoel->MaterialId() != this->GetElemMaterialID()) DebugStop(); //itapopo verificar se usará _assert_
244
245 /*Divide geoel*/
246 TPZVec<TPZGeoEl *> Sons;
247 geoel->Divide(Sons);
248
249 /*If a 1D segment is neighbor, it must be divided too*/
250 if(this->elementswidth != 3) DebugStop(); //itapopo verificar o segment para malha 3D
251
252 std::vector<int> sides(3);
253 sides[0] = 3; sides[1] = 4; sides[2] = 5;
254 for(int j = 0; j < sides.size(); j++ ){
255
256 TPZGeoElSide Neighbour = geoel->Neighbour(sides[j]);
257
258 if( Neighbour.Element()->MaterialId() == this->GetBoundaryMaterialID() && !Neighbour.Element()->HasSubElement() ){
259 TPZVec<TPZGeoEl *> pv2;
260 Neighbour.Element()->Divide(pv2);
261 }
262 }
263 }
264
265 gmesh->BuildConnectivity();
266
267}
268/*}}}*/
269void AdaptiveMeshRefinement::RefineMeshToAvoidHangingNodes(TPZGeoMesh *gmesh){/*{{{*/
270
271 _printf_("\trefine to avoid hanging nodes...\n");
272 /*Refine elements to avoid hanging nodes: non-uniform refinement*/
273 const int NElem = gmesh->NElements();
274 for(int i = 0; i < NElem; i++){
275
276 /*Get geometric element and verify if it has already been refined. Geoel may not have been previously refined*/
277 TPZGeoEl * geoel = gmesh->Element(i);
278 if(!geoel) continue;
279 if(geoel->HasSubElement()) continue;
280 if(geoel->MaterialId() != this->GetElemMaterialID()) continue;
281
282 /*Get the refinement pattern for this element and refine it*/
283 TPZAutoPointer<TPZRefPattern> refp = TPZRefPatternTools::PerfectMatchRefPattern(geoel);
284 if(refp){
285 TPZVec<TPZGeoEl *> Sons;
286 geoel->SetRefPattern(refp);
287 geoel->Divide(Sons);
288 }
289
290 }
291
292 gmesh->BuildConnectivity();
293
294}
295/*}}}*/
296void AdaptiveMeshRefinement::GetMesh(TPZGeoMesh *gmesh,int &nvertices,int &nelements,int &nsegments,double** px,double** py,double** pz, int** pelements, int** psegments){/*{{{*/
297
298 /*IMPORTANT! pelements (and psegments) are in Matlab indexing*/
299 /*NEOPZ works only in C indexing*/
300
301 /* vertices */
302 int ntotalvertices = gmesh->NNodes();//total
303
304 /* mesh coords */
305 double* newmeshX = xNew<IssmDouble>(ntotalvertices);
306 double* newmeshY = xNew<IssmDouble>(ntotalvertices);
307 double* newmeshZ = xNew<IssmDouble>(ntotalvertices);
308
309 /* getting mesh coords */
310 for(int i = 0; i < ntotalvertices; i++ ){
311 TPZVec<REAL> coords(3,0.);
312 gmesh->NodeVec()[i].GetCoordinates(coords);
313 newmeshX[i] = coords[0];
314 newmeshY[i] = coords[1];
315 newmeshZ[i] = coords[2];
316 }
317
318 /* elements */
319 std::vector<TPZGeoEl*> GeoVec; GeoVec.clear();
320 for(int i = 0; i < gmesh->NElements(); i++){
321 if( gmesh->ElementVec()[i]->HasSubElement() ) continue;
322 if( gmesh->ElementVec()[i]->MaterialId() != this->GetElemMaterialID() ) continue;
323 GeoVec.push_back( gmesh->ElementVec()[i]);
324 }
325
326 int ntotalelements = (int)GeoVec.size();
327 int* newelements = xNew<int>(ntotalelements*this->elementswidth);
328
329 if ( !(this->elementswidth == 3) && !(this->elementswidth == 4) && !(this->elementswidth == 6) ) DebugStop();
330
331 for(int i=0;i<GeoVec.size();i++){
332 for(int j=0;j<this->elementswidth;j++) newelements[i*this->elementswidth+j]=(int)GeoVec[i]->NodeIndex(j)+1;//C to Matlab indexing
333 }
334
335 /* segments */
336 std::vector<TPZGeoEl*> SegVec; SegVec.clear();
337 for(int i = 0; i < gmesh->NElements(); i++){
338 if( gmesh->ElementVec()[i]->HasSubElement() ) continue;
339 if( gmesh->ElementVec()[i]->MaterialId() != this->GetBoundaryMaterialID() ) continue;
340 SegVec.push_back( gmesh->ElementVec()[i]);
341 }
342
343 int ntotalsegments = (int)SegVec.size();
344 int *newsegments=NULL;
345 if(ntotalsegments>0) newsegments=xNew<int>(ntotalsegments*3);
346
347 for(int i=0;i<SegVec.size();i++){
348
349 for(int j=0;j<2;j++) newsegments[i*3+j]=(int)SegVec[i]->NodeIndex(j)+1;//C to Matlab indexing
350
351 int neighborindex = SegVec[i]->Neighbour(2).Element()->Index();
352 int neighbourid = -1;
353
354 for(int j = 0; j < GeoVec.size(); j++){
355 if( GeoVec[j]->Index() == neighborindex || GeoVec[j]->FatherIndex() == neighborindex){
356 neighbourid = j;
357 break;
358 }
359 }
360
361 if(neighbourid==-1) DebugStop(); //itapopo talvez passar para _assert_
362 newsegments[i*3+2] = neighbourid+1;//C to Matlab indexing
363 }
364
365 //setting outputs
366 nvertices = ntotalvertices;
367 nelements = ntotalelements;
368 nsegments = ntotalsegments;
369 *px = newmeshX;
370 *py = newmeshY;
371 *pz = newmeshZ;
372 *pelements = newelements;
373 *psegments = newsegments;
374
375}
376/*}}}*/
377void AdaptiveMeshRefinement::CalcGroundingLinePosition(double *masklevelset,std::vector<TPZVec<REAL> > &GLvec){/*{{{*/
378
379 /* Find grounding line using elments center point */
380 GLvec.clear();
381 for(int i=0;i<this->previousmesh->NElements();i++){
382
383 if(this->previousmesh->Element(i)->MaterialId()!=this->GetElemMaterialID()) continue;
384 if(this->previousmesh->Element(i)->HasSubElement()) continue;
385
386 //itapopo apenas malha 2D triangular!
387 int vertex0 = this->previousmesh->Element(i)->NodeIndex(0);
388 int vertex1 = this->previousmesh->Element(i)->NodeIndex(1);
389 int vertex2 = this->previousmesh->Element(i)->NodeIndex(2);
390
391 //itapopo inserir uma verificação para não acessar fora da memória
392 double mls0 = masklevelset[vertex0];
393 double mls1 = masklevelset[vertex1];
394 double mls2 = masklevelset[vertex2];
395
396 if( mls0*mls1 < 0. || mls1*mls2 < 0. ){
397 const int side = 6;
398 TPZVec<double> qsi(2,0.);
399 TPZVec<double> X(3,0.);
400 this->previousmesh->Element(i)->CenterPoint(side, qsi);
401 this->previousmesh->Element(i)->X(qsi, X);
402 GLvec.push_back(X);
403 }
404 }
405
406// itapopo apenas para debugar
407// std::ofstream fileGL("/Users/santos/Desktop/gl.nb");
408// fileGL << "ListPlot[{";
409// for(int i = 0; i < GLvec.size(); i++){
410// fileGL << "{" << GLvec[i][0] << "," << GLvec[i][1] << /*"," << 0. << */"}";
411// if(i != GLvec.size()-1) fileGL << ",";
412// }
413// fileGL << "}]";
414// fileGL.flush();
415// fileGL.close();
416
417}
418/*}}}*/
419void AdaptiveMeshRefinement::SetElementsToRefine(TPZGeoMesh *gmesh,std::vector<TPZVec<REAL> > &GLvec,int &hlevel,std::vector<int> &ElemVec){/*{{{*/
420
421 if(!gmesh) DebugStop(); //itapopo verificar se usará _assert_
422
423 ElemVec.clear();
424
425 // itapopo inserir modo de encontrar criterio TESTING!!!! Come to false!
426 if(false) this->TagAllElements(gmesh,ElemVec); //uniform, refine all elements!
427
428 /* Adaptive refinement. This refines some elements following some criteria*/
429 this->TagElementsNearGroundingLine(gmesh, GLvec, hlevel, ElemVec);
430
431}
432/*}}}*/
433void AdaptiveMeshRefinement::TagAllElements(TPZGeoMesh *gmesh,std::vector<int> &ElemVec){/*{{{*/
434
435 /* Uniform refinement. This refines the entire mesh */
436 int nelements = gmesh->NElements();
437 for(int i=0;i<nelements;i++){
438 if(gmesh->Element(i)->MaterialId()!=this->GetElemMaterialID()) continue;
439 if(gmesh->Element(i)->HasSubElement()) continue;
440 ElemVec.push_back(i);
441 }
442}
443/*}}}*/
444void AdaptiveMeshRefinement::TagElementsNearGroundingLine(TPZGeoMesh *gmesh,std::vector<TPZVec<REAL> > &GLvec,int &hlevel,std::vector<int> &ElemVec){/*{{{*/
445
446 /* Tag elements near grounding line */
447 double D1 = this->regionlevel1;
448 double Dhmax = this->regionlevelmax;
449 int hmax = this->levelmax;
450 double alpha = (hmax==1) ? 0. : log(D1/Dhmax)/(hmax-1.);
451 double Di = D1/exp(alpha*(hlevel-1));
452
453 for(int i=0;i<gmesh->NElements();i++){
454
455 if(gmesh->Element(i)->MaterialId()!=this->GetElemMaterialID()) continue;
456 if(gmesh->Element(i)->HasSubElement()) continue;
457 if(gmesh->Element(i)->Level()>=hlevel) continue;
458
459 const int side2D = 6;
460 TPZVec<REAL> qsi(2,0.);
461 TPZVec<REAL> centerPoint(3,0.);
462 gmesh->Element(i)->CenterPoint(side2D, qsi);
463 gmesh->Element(i)->X(qsi, centerPoint);
464
465 REAL distance = Di;
466
467 for (int j = 0; j < GLvec.size(); j++) {
468
469 REAL value = ( GLvec[j][0] - centerPoint[0] ) * ( GLvec[j][0] - centerPoint[0] ); // (x2-x1)^2
470 value += ( GLvec[j][1] - centerPoint[1] ) * ( GLvec[j][1] - centerPoint[1] );// (y2-y1)^2
471 value = std::sqrt(value); //Radius
472
473 //finding the min distance to the grounding line
474 if(value < distance) distance = value;
475
476 }
477
478 if(distance < Di) ElemVec.push_back(i);
479 }
480
481}
482/*}}}*/
483void AdaptiveMeshRefinement::CreateInitialMesh(int &nvertices,int &nelements,int &nsegments,int &width,double* x,double* y,double* z,int* elements,int* segments){/*{{{*/
484
485 /*IMPORTANT! elements come in Matlab indexing*/
486 /*NEOPZ works only in C indexing*/
487
488 _assert_(nvertices>0);
489 _assert_(nelements>0);
490 this->SetElementWidth(width);
491
492 /*Verify and creating initial mesh*/
493 if(this->fathermesh) _error_("Initial mesh already exists!");
494
495 this->fathermesh = new TPZGeoMesh();
496 this->fathermesh->NodeVec().Resize( nvertices );
497
498 /*Set the vertices (geometric nodes in NeoPZ context)*/
499 for(int i=0;i<nvertices;i++){
500 /*x,y,z coords*/
501 TPZManVector<REAL,3> coord(3,0.);
502 coord[0]= x[i];
503 coord[1]= y[i];
504 coord[2]= z[i];
505 /*Insert in the mesh*/
506 this->fathermesh->NodeVec()[i].SetCoord(coord);
507 this->fathermesh->NodeVec()[i].SetNodeId(i);
508 }
509
510 /*Generate the elements*/
511 long index;
512 const int mat = this->GetElemMaterialID();
513 TPZManVector<long> elem(this->elementswidth,0);
514
515 for(int iel=0;iel<nelements;iel++){
516
517 for(int jel=0;jel<this->elementswidth;jel++) elem[jel]=elements[iel*this->elementswidth+jel]-1;//Convert Matlab to C indexing
518
519 /*reftype = 0: uniform, fast / reftype = 1: uniform and non-uniform (avoid hanging nodes), it is not too fast */
520 const int reftype = 0;
521 switch(this->elementswidth){
522 case 3: this->fathermesh->CreateGeoElement(ETriangle, elem, mat, index, reftype); break;
523 case 4: this->fathermesh->CreateGeoElement(ETetraedro, elem, mat, index, reftype); DebugStop(); break;
524 case 6: this->fathermesh->CreateGeoElement(EPrisma, elem, mat, index, reftype); DebugStop(); break;
525 default: DebugStop();//itapopo _error_("mesh not supported yet");
526 }
527
528 /*Define the element ID*/
529 this->fathermesh->ElementVec()[index]->SetId(iel);
530 }
531
532 /*Generate the 1D segments elements (boundary)*/
533 const int matboundary = this->GetBoundaryMaterialID();
534 TPZManVector<long> boundary(2,0.);
535
536 for(int iel=nelements;iel<nelements+nsegments;iel++){
537 boundary[0] = segments[(iel-nelements)*2+0]-1;//Convert Matlab to C indexing
538 boundary[1] = segments[(iel-nelements)*2+1]-1;//Convert Matlab to C indexing
539 /*reftype = 0: uniform, fast / reftype = 1: uniform and non-uniform (avoid hanging nodes), it is not too fast */
540 const int reftype = 0;
541 this->fathermesh->CreateGeoElement(EOned, boundary, matboundary, index, reftype);//cria elemento unidimensional
542 this->fathermesh->ElementVec()[index]->SetId(iel);
543 }
544
545 /*Build element and node connectivities*/
546 this->fathermesh->BuildConnectivity();
547
548 /*Create previous mesh as a copy of father mesh*/
549 this->previousmesh = new TPZGeoMesh(*this->fathermesh);
550
551}
552/*}}}*/
553#include "pzgeotriangle.h" //itapopo
554#include "pzreftriangle.h" //itapopo
555using namespace pzgeom;
556TPZGeoMesh* AdaptiveMeshRefinement::CreateRefPatternMesh(TPZGeoMesh* gmesh){/*{{{*/
557
558 TPZGeoMesh *newgmesh = new TPZGeoMesh();
559 newgmesh->CleanUp();
560
561 int nnodes = gmesh->NNodes();
562 int nelem = gmesh->NElements();
563 int mat = this->GetElemMaterialID();;
564 int reftype = 1;
565 long index;
566
567 //nodes
568 newgmesh->NodeVec().Resize(nnodes);
569 for(int i=0;i<nnodes;i++) newgmesh->NodeVec()[i] = gmesh->NodeVec()[i];
570
571 //elements
572 for(int i=0;i<nelem;i++){
573 TPZGeoEl * geoel = gmesh->Element(i);
574 TPZManVector<long> elem(3,0);
575 for(int j=0;j<3;j++) elem[j] = geoel->NodeIndex(j);
576
577 newgmesh->CreateGeoElement(ETriangle,elem,mat,index,reftype);
578 newgmesh->ElementVec()[index]->SetId(geoel->Id());
579
580 TPZGeoElRefPattern<TPZGeoTriangle>* newgeoel = dynamic_cast<TPZGeoElRefPattern<TPZGeoTriangle>*>(newgmesh->ElementVec()[index]);
581
582 //old neighbourhood
583 const int nsides = TPZGeoTriangle::NSides;
584 TPZVec< std::vector<TPZGeoElSide> > neighbourhood(nsides);
585 TPZVec<long> NodesSequence(0);
586 for(int s = 0; s < nsides; s++){
587 neighbourhood[s].resize(0);
588 TPZGeoElSide mySide(geoel,s);
589 TPZGeoElSide neighS = mySide.Neighbour();
590 if(mySide.Dimension() == 0){
591 long oldSz = NodesSequence.NElements();
592 NodesSequence.resize(oldSz+1);
593 NodesSequence[oldSz] = geoel->NodeIndex(s);
594 }
595 while(mySide != neighS){
596 neighbourhood[s].push_back(neighS);
597 neighS = neighS.Neighbour();
598 }
599 }
600
601 //inserting in new element
602 for(int s = 0; s < nsides; s++){
603 TPZGeoEl * tempEl = newgeoel;
604 TPZGeoElSide tempSide(newgeoel,s);
605 int byside = s;
606 for(unsigned long n = 0; n < neighbourhood[s].size(); n++){
607 TPZGeoElSide neighS = neighbourhood[s][n];
608 tempEl->SetNeighbour(byside, neighS);
609 tempEl = neighS.Element();
610 byside = neighS.Side();
611 }
612 tempEl->SetNeighbour(byside, tempSide);
613 }
614
615 long fatherindex = geoel->FatherIndex();
616 if(fatherindex>-1) newgeoel->SetFather(fatherindex);
617
618 if(!geoel->HasSubElement()) continue;
619
620 int nsons = geoel->NSubElements();
621
622 TPZAutoPointer<TPZRefPattern> ref = gRefDBase.GetUniformRefPattern(ETriangle);
623 newgeoel->SetRefPattern(ref);
624
625 for(int j=0;j<nsons;j++){
626 TPZGeoEl* son = geoel->SubElement(j);
627 if(!son){
628 DebugStop();
629 }
630 newgeoel->SetSubElement(j,son);
631 }
632 }
633 newgmesh->BuildConnectivity();
634
635 return newgmesh;
636}
637/*}}}*/
638void AdaptiveMeshRefinement::SetLevelMax(int &h){/*{{{*/
639 this->levelmax = h;
640}
641/*}}}*/
642void AdaptiveMeshRefinement::SetRegions(double &D1,double Dhmax){/*{{{*/
643 this->regionlevel1 = D1;
644 this->regionlevelmax = Dhmax;
645}
646/*}}}*/
647void AdaptiveMeshRefinement::SetElementWidth(int &width){/*{{{*/
648 this->elementswidth = width;
649}
650/*}}}*/
651void AdaptiveMeshRefinement::CheckMesh(int &nvertices,int &nelements,int &nsegments,int &width,double** px,double** py,double** pz,int** pelements, int** psegments){/*{{{*/
652
653 /*Basic verification*/
654 if( !(nvertices > 0) || !(nelements > 0) ) DebugStop(); //itapopo verificar se irá usar o _assert_
655
656 if ( !(width == 3) && !(width == 4) && !(width == 6) ) DebugStop(); // itapopo verifcar se irá usar o _assert_
657
658 if( !px || !py || !pz || !pelements ) DebugStop(); // itapopo verifcar se irá usar o _assert_
659
660 /*Verify if there are orphan nodes*/
661 std::set<int> elemvertices;
662 elemvertices.clear();
663 for(int i = 0; i < nelements; i++){
664 for(int j = 0; j < width; j++) {
665 elemvertices.insert((*pelements)[i*width+j]);
666 }
667 }
668
669 if( elemvertices.size() != nvertices ) DebugStop();//itapopo verificar se irá usar o _assert_
670
671 //Verify if there are inf or NaN in coords
672 for(int i = 0; i < nvertices; i++) if(isnan((*px)[i]) || isinf((*px)[i])) DebugStop();
673 for(int i = 0; i < nvertices; i++) if(isnan((*py)[i]) || isinf((*py)[i])) DebugStop();
674 for(int i = 0; i < nvertices; i++) if(isnan((*pz)[i]) || isinf((*pz)[i])) DebugStop();
675
676 for(int i = 0; i < nelements; i++){
677 for(int j = 0; j < width; j++){
678 if( isnan((*pelements)[i*width+j]) || isinf((*pelements)[i*width+j]) ) DebugStop();
679 }
680 }
681
682}
683/*}}}*/
Note: See TracBrowser for help on using the repository browser.