source: issm/oecreview/Archive/21337-21723/ISSM-21671-21672.diff@ 21726

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

CHG added Archive/21337-21723
File size: 29.1 KB
RevLine 
[21726]1Index: ../trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp
2===================================================================
3--- ../trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp (revision 21671)
4+++ ../trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp (revision 21672)
5@@ -26,15 +26,30 @@
6 this->CleanUp();
7
8 /*Copy all data*/
9- this->fathermesh = new TPZGeoMesh(*cp.fathermesh);
10- this->previousmesh = new TPZGeoMesh(*cp.previousmesh);
11- this->hmax = cp.hmax;
12- this->elementswidth = cp.elementswidth;
13+ this->fathermesh = new TPZGeoMesh(*cp.fathermesh);
14+ this->previousmesh = new TPZGeoMesh(*cp.previousmesh);
15+ this->levelmax = cp.levelmax;
16+ this->elementswidth = cp.elementswidth;
17+ this->regionlevel1 = cp.regionlevel1;
18+ this->regionlevelmax = cp.regionlevelmax;
19 return *this;
20
21 }
22 /*}}}*/
23 AdaptiveMeshRefinement::~AdaptiveMeshRefinement(){/*{{{*/
24+
25+ bool ismismip = true;
26+ if(ismismip){//itapopo
27+ TPZFileStream fstr;
28+ std::stringstream ss;
29+
30+ ss << this->levelmax;
31+ std::string AMRfile = "/home/santos/Misomip2/L" + ss.str() + "_tsai/amr.txt";
32+
33+ fstr.OpenWrite(AMRfile.c_str());
34+ int withclassid = 1;
35+ this->Write(fstr,withclassid);
36+ }
37 this->CleanUp();
38 gRefDBase.clear();
39 }
40@@ -44,17 +59,21 @@
41 /*Verify and delete all data*/
42 if(this->fathermesh) delete this->fathermesh;
43 if(this->previousmesh) delete this->previousmesh;
44- this->hmax=-1;
45- this->elementswidth=-1;
46+ this->levelmax = -1;
47+ this->elementswidth = -1;
48+ this->regionlevel1 = -1;
49+ this->regionlevelmax = -1;
50 }
51 /*}}}*/
52 void AdaptiveMeshRefinement::Initialize(){/*{{{*/
53
54 /*Set pointers to NULL*/
55- this->fathermesh = NULL;
56- this->previousmesh = NULL;
57- this->hmax = -1;
58- this->elementswidth = -1;
59+ this->fathermesh = NULL;
60+ this->previousmesh = NULL;
61+ this->levelmax = -1;
62+ this->elementswidth = -1;
63+ this->regionlevel1 = -1;
64+ this->regionlevelmax = -1;
65 }
66 /*}}}*/
67 int AdaptiveMeshRefinement::ClassId() const{/*{{{*/
68@@ -81,10 +100,12 @@
69 }
70
71 /* Read simple attributes */
72- buf.Read(&this->hmax,1);
73+ buf.Read(&this->levelmax,1);
74 buf.Read(&this->elementswidth,1);
75+ buf.Read(&this->regionlevel1,1);
76+ buf.Read(&this->regionlevelmax,1);
77
78- /* Read geometric mesh*/
79+ /* Read geometric mesh*/
80 TPZSaveable *sv1 = TPZSaveable::Restore(buf,0);
81 this->fathermesh = dynamic_cast<TPZGeoMesh*>(sv1);
82
83@@ -113,9 +134,11 @@
84 buf.Write(&classid,1);
85
86 /* Write simple attributes */
87- buf.Write(&this->hmax,1);
88+ buf.Write(&this->levelmax,1);
89 buf.Write(&this->elementswidth,1);
90-
91+ buf.Write(&this->regionlevel1,1);
92+ buf.Write(&this->regionlevelmax,1);
93+
94 /* Write the geometric mesh*/
95 this->fathermesh->Write(buf, this->ClassId());
96 this->previousmesh->Write(buf, this->ClassId());
97@@ -137,15 +160,15 @@
98 /*IMPORTANT! pelements (and psegments) are in Matlab indexing*/
99 /*NEOPZ works only in C indexing*/
100
101- //itapopo _assert_(this->fathermesh);
102- //itapopo _assert_(this->previousmesh);
103+ _assert_(this->fathermesh);
104+ _assert_(this->previousmesh);
105
106 /*Calculate the position of the grounding line using previous mesh*/
107 std::vector<TPZVec<REAL> > GLvec;
108 this->CalcGroundingLinePosition(masklevelset, GLvec);
109
110- std::ofstream file1("/ronne_1/home/santos/mesh0.vtk");
111- TPZVTKGeoMesh::PrintGMeshVTK(this->previousmesh,file1 );
112+ // std::ofstream file1("/home/santos/mesh0.vtk");
113+ // TPZVTKGeoMesh::PrintGMeshVTK(this->fathermesh,file1 );
114
115 /*run refinement or unrefinement process*/
116 TPZGeoMesh *newmesh;
117@@ -157,8 +180,8 @@
118
119 this->RefinementProcess(newmesh,GLvec);
120
121- std::ofstream file2("/ronne_1/home/santos/mesh1.vtk");
122- TPZVTKGeoMesh::PrintGMeshVTK(this->previousmesh,file2 );
123+ //std::ofstream file2("/home/santos/mesh1.vtk");
124+ //TPZVTKGeoMesh::PrintGMeshVTK(this->previousmesh,file2 );
125
126 /*Set new mesh pointer. Previous mesh just have uniform elements*/
127 if(type_process==1){
128@@ -167,15 +190,16 @@
129 }
130
131 /*Refine elements to avoid hanging nodes*/
132- TPZGeoMesh *nohangingnodesmesh = new TPZGeoMesh(*newmesh);
133+ //TPZGeoMesh *nohangingnodesmesh = new TPZGeoMesh(*newmesh);//itapopo testando, este era o original
134+ TPZGeoMesh *nohangingnodesmesh = this->CreateRefPatternMesh(newmesh);//itapopo testando, este eh novo metodo
135
136- std::ofstream file3("/ronne_1/home/santos/mesh2.vtk");
137- TPZVTKGeoMesh::PrintGMeshVTK(newmesh,file3);
138+ //std::ofstream file3("/home/santos/mesh2.vtk");
139+ //TPZVTKGeoMesh::PrintGMeshVTK(this->previousmesh,file3);
140
141 this->RefineMeshToAvoidHangingNodes(nohangingnodesmesh);
142
143- std::ofstream file4("/ronne_1/home/santos/mesh3.vtk");
144- TPZVTKGeoMesh::PrintGMeshVTK(nohangingnodesmesh,file4);
145+ //std::ofstream file4("/home/santos/mesh3.vtk");
146+ //TPZVTKGeoMesh::PrintGMeshVTK(nohangingnodesmesh,file4);
147
148 /*Get new geometric mesh in ISSM data structure*/
149 this->GetMesh(nohangingnodesmesh,nvertices,nelements,nsegments,px,py,pz,pelements,psegments);
150@@ -190,10 +214,10 @@
151 /*}}}*/
152 void AdaptiveMeshRefinement::RefinementProcess(TPZGeoMesh *gmesh,std::vector<TPZVec<REAL> > &GLvec){/*{{{*/
153
154- /*Refine mesh hmax times*/
155- _printf_("\n\trefinement process (hmax = " << this->hmax << ")\n");
156+ /*Refine mesh levelmax times*/
157+ _printf_("\n\trefinement process (level max = " << this->levelmax << ")\n");
158 _printf_("\tprogress: ");
159- for(int hlevel=1;hlevel<=this->hmax;hlevel++){
160+ for(int hlevel=1;hlevel<=this->levelmax;hlevel++){
161
162 /*Set elements to be refined using some criteria*/
163 std::vector<int> ElemVec; //elements without children
164@@ -364,6 +388,7 @@
165 int vertex1 = this->previousmesh->Element(i)->NodeIndex(1);
166 int vertex2 = this->previousmesh->Element(i)->NodeIndex(2);
167
168+ //itapopo inserir uma verificação para não acessar fora da memória
169 double mls0 = masklevelset[vertex0];
170 double mls1 = masklevelset[vertex1];
171 double mls2 = masklevelset[vertex2];
172@@ -418,10 +443,12 @@
173 /*}}}*/
174 void AdaptiveMeshRefinement::TagElementsNearGroundingLine(TPZGeoMesh *gmesh,std::vector<TPZVec<REAL> > &GLvec,int &hlevel,std::vector<int> &ElemVec){/*{{{*/
175
176- /* Tag elements near grounding line */
177- double MaxRegion = 40000.; //itapopo
178- double alpha = log(1.5); //itapopo
179- double MaxDistance = MaxRegion / std::exp(alpha*(hlevel-1));
180+ /* Tag elements near grounding line */
181+ double D1 = this->regionlevel1;
182+ double Dhmax = this->regionlevelmax;
183+ int hmax = this->levelmax;
184+ double alpha = (hmax==1) ? 0. : log(D1/Dhmax)/(hmax-1.);
185+ double Di = D1/exp(alpha*(hlevel-1));
186
187 for(int i=0;i<gmesh->NElements();i++){
188
189@@ -435,20 +462,20 @@
190 gmesh->Element(i)->CenterPoint(side2D, qsi);
191 gmesh->Element(i)->X(qsi, centerPoint);
192
193- REAL distance = MaxDistance;
194+ REAL distance = Di;
195
196 for (int j = 0; j < GLvec.size(); j++) {
197
198 REAL value = ( GLvec[j][0] - centerPoint[0] ) * ( GLvec[j][0] - centerPoint[0] ); // (x2-x1)^2
199 value += ( GLvec[j][1] - centerPoint[1] ) * ( GLvec[j][1] - centerPoint[1] );// (y2-y1)^2
200- value = std::sqrt(value); ///Radius
201+ value = std::sqrt(value); //Radius
202
203 //finding the min distance to the grounding line
204 if(value < distance) distance = value;
205
206 }
207
208- if(distance < MaxDistance) ElemVec.push_back(i);
209+ if(distance < Di) ElemVec.push_back(i);
210 }
211
212 }
213@@ -458,12 +485,12 @@
214 /*IMPORTANT! elements come in Matlab indexing*/
215 /*NEOPZ works only in C indexing*/
216
217- // itapopo _assert_(nvertices>0);
218- // itapoo _assert_(nelements>0);
219+ _assert_(nvertices>0);
220+ _assert_(nelements>0);
221 this->SetElementWidth(width);
222
223 /*Verify and creating initial mesh*/
224- //itapopo if(this->fathermesh) _error_("Initial mesh already exists!");
225+ if(this->fathermesh) _error_("Initial mesh already exists!");
226
227 this->fathermesh = new TPZGeoMesh();
228 this->fathermesh->NodeVec().Resize( nvertices );
229@@ -490,7 +517,7 @@
230 for(int jel=0;jel<this->elementswidth;jel++) elem[jel]=elements[iel*this->elementswidth+jel]-1;//Convert Matlab to C indexing
231
232 /*reftype = 0: uniform, fast / reftype = 1: uniform and non-uniform (avoid hanging nodes), it is not too fast */
233- const int reftype = 1;
234+ const int reftype = 0;
235 switch(this->elementswidth){
236 case 3: this->fathermesh->CreateGeoElement(ETriangle, elem, mat, index, reftype); break;
237 case 4: this->fathermesh->CreateGeoElement(ETetraedro, elem, mat, index, reftype); DebugStop(); break;
238@@ -520,54 +547,107 @@
239
240 /*Create previous mesh as a copy of father mesh*/
241 this->previousmesh = new TPZGeoMesh(*this->fathermesh);
242-
243- /*Set refinement patterns*/
244- this->SetRefPatterns();
245
246 }
247 /*}}}*/
248-void AdaptiveMeshRefinement::SetHMax(int &h){/*{{{*/
249- this->hmax = h;
250+#include "pzgeotriangle.h" //itapopo
251+#include "pzreftriangle.h" //itapopo
252+using namespace pzgeom;
253+TPZGeoMesh* AdaptiveMeshRefinement::CreateRefPatternMesh(TPZGeoMesh* gmesh){/*{{{*/
254+
255+ TPZGeoMesh *newgmesh = new TPZGeoMesh();
256+ newgmesh->CleanUp();
257+
258+ int nnodes = gmesh->NNodes();
259+ int nelem = gmesh->NElements();
260+ int mat = this->GetElemMaterialID();;
261+ int reftype = 1;
262+ long index;
263+
264+ //nodes
265+ newgmesh->NodeVec().Resize(nnodes);
266+ for(int i=0;i<nnodes;i++) newgmesh->NodeVec()[i] = gmesh->NodeVec()[i];
267+
268+ //elements
269+ for(int i=0;i<nelem;i++){
270+ TPZGeoEl * geoel = gmesh->Element(i);
271+ TPZManVector<long> elem(3,0);
272+ for(int j=0;j<3;j++) elem[j] = geoel->NodeIndex(j);
273+
274+ newgmesh->CreateGeoElement(ETriangle,elem,mat,index,reftype);
275+ newgmesh->ElementVec()[index]->SetId(geoel->Id());
276+
277+ TPZGeoElRefPattern<TPZGeoTriangle>* newgeoel = dynamic_cast<TPZGeoElRefPattern<TPZGeoTriangle>*>(newgmesh->ElementVec()[index]);
278+
279+ //old neighbourhood
280+ const int nsides = TPZGeoTriangle::NSides;
281+ TPZVec< std::vector<TPZGeoElSide> > neighbourhood(nsides);
282+ TPZVec<long> NodesSequence(0);
283+ for(int s = 0; s < nsides; s++){
284+ neighbourhood[s].resize(0);
285+ TPZGeoElSide mySide(geoel,s);
286+ TPZGeoElSide neighS = mySide.Neighbour();
287+ if(mySide.Dimension() == 0){
288+ long oldSz = NodesSequence.NElements();
289+ NodesSequence.resize(oldSz+1);
290+ NodesSequence[oldSz] = geoel->NodeIndex(s);
291+ }
292+ while(mySide != neighS){
293+ neighbourhood[s].push_back(neighS);
294+ neighS = neighS.Neighbour();
295+ }
296+ }
297+
298+ //inserting in new element
299+ for(int s = 0; s < nsides; s++){
300+ TPZGeoEl * tempEl = newgeoel;
301+ TPZGeoElSide tempSide(newgeoel,s);
302+ int byside = s;
303+ for(unsigned long n = 0; n < neighbourhood[s].size(); n++){
304+ TPZGeoElSide neighS = neighbourhood[s][n];
305+ tempEl->SetNeighbour(byside, neighS);
306+ tempEl = neighS.Element();
307+ byside = neighS.Side();
308+ }
309+ tempEl->SetNeighbour(byside, tempSide);
310+ }
311+
312+ long fatherindex = geoel->FatherIndex();
313+ if(fatherindex>-1) newgeoel->SetFather(fatherindex);
314+
315+ if(!geoel->HasSubElement()) continue;
316+
317+ int nsons = geoel->NSubElements();
318+
319+ TPZAutoPointer<TPZRefPattern> ref = gRefDBase.GetUniformRefPattern(ETriangle);
320+ newgeoel->SetRefPattern(ref);
321+
322+ for(int j=0;j<nsons;j++){
323+ TPZGeoEl* son = geoel->SubElement(j);
324+ if(!son){
325+ DebugStop();
326+ }
327+ newgeoel->SetSubElement(j,son);
328+ }
329+ }
330+ newgmesh->BuildConnectivity();
331+
332+ return newgmesh;
333 }
334 /*}}}*/
335+void AdaptiveMeshRefinement::SetLevelMax(int &h){/*{{{*/
336+ this->levelmax = h;
337+}
338+/*}}}*/
339+void AdaptiveMeshRefinement::SetRegions(double &D1,double Dhmax){/*{{{*/
340+ this->regionlevel1 = D1;
341+ this->regionlevelmax = Dhmax;
342+}
343+/*}}}*/
344 void AdaptiveMeshRefinement::SetElementWidth(int &width){/*{{{*/
345 this->elementswidth = width;
346 }
347 /*}}}*/
348-void AdaptiveMeshRefinement::SetRefPatterns(){/*{{{*/
349-
350- /*Initialize the global variable of refinement patterns*/
351- gRefDBase.InitializeUniformRefPattern(EOned);
352- gRefDBase.InitializeUniformRefPattern(ETriangle);
353-
354- //gRefDBase.InitializeRefPatterns();
355- /*Insert specifics patterns to ISSM core*/
356- std::string filepath = REFPATTERNDIR;
357- std::string filename1 = filepath + "/2D_Triang_Rib_3.rpt";
358- std::string filename2 = filepath + "/2D_Triang_Rib_4.rpt";
359- std::string filename3 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_4.rpt";
360- std::string filename4 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_4_permuted.rpt";
361- std::string filename5 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_5.rpt";
362- std::string filename6 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_5_permuted.rpt";
363- std::string filename7 = filepath + "/2D_Triang_Rib_5.rpt";
364-
365- TPZAutoPointer<TPZRefPattern> refpat1 = new TPZRefPattern(filename1);
366- TPZAutoPointer<TPZRefPattern> refpat2 = new TPZRefPattern(filename2);
367- TPZAutoPointer<TPZRefPattern> refpat3 = new TPZRefPattern(filename3);
368- TPZAutoPointer<TPZRefPattern> refpat4 = new TPZRefPattern(filename4);
369- TPZAutoPointer<TPZRefPattern> refpat5 = new TPZRefPattern(filename5);
370- TPZAutoPointer<TPZRefPattern> refpat6 = new TPZRefPattern(filename6);
371- TPZAutoPointer<TPZRefPattern> refpat7 = new TPZRefPattern(filename7);
372-
373- if(!gRefDBase.FindRefPattern(refpat1)) gRefDBase.InsertRefPattern(refpat1);
374- if(!gRefDBase.FindRefPattern(refpat2)) gRefDBase.InsertRefPattern(refpat2);
375- if(!gRefDBase.FindRefPattern(refpat3)) gRefDBase.InsertRefPattern(refpat3);
376- if(!gRefDBase.FindRefPattern(refpat4)) gRefDBase.InsertRefPattern(refpat4);
377- if(!gRefDBase.FindRefPattern(refpat5)) gRefDBase.InsertRefPattern(refpat5);
378- if(!gRefDBase.FindRefPattern(refpat6)) gRefDBase.InsertRefPattern(refpat6);
379- if(!gRefDBase.FindRefPattern(refpat7)) gRefDBase.InsertRefPattern(refpat7);
380-}
381-/*}}}*/
382 void AdaptiveMeshRefinement::CheckMesh(int &nvertices,int &nelements,int &nsegments,int &width,double** px,double** py,double** pz,int** pelements, int** psegments){/*{{{*/
383
384 /*Basic verification*/
385Index: ../trunk-jpl/src/c/classes/FemModel.cpp
386===================================================================
387--- ../trunk-jpl/src/c/classes/FemModel.cpp (revision 21671)
388+++ ../trunk-jpl/src/c/classes/FemModel.cpp (revision 21672)
389@@ -2920,28 +2920,53 @@
390 void FemModel::InitializeAdaptiveRefinement(void){/*{{{*/
391
392 /*Define variables*/
393- int my_rank = IssmComm::GetRank();
394- this->amr = NULL;//initialize amr as NULL
395- int numberofvertices = this->vertices->NumberOfVertices();
396- int numberofelements = this->elements->NumberOfElements();
397- int numberofsegments = 0; //used on matlab
398- IssmDouble* x = NULL;
399- IssmDouble* y = NULL;
400- IssmDouble* z = NULL;
401- int* elements = NULL;
402- int elementswidth = this->GetElementsWidth(); //just tria elements in this version. Itapopo:
403+ int my_rank = IssmComm::GetRank();
404+ this->amr = NULL;//initialize amr as NULL
405+ int numberofvertices = this->vertices->NumberOfVertices();
406+ int numberofelements = this->elements->NumberOfElements();
407+ int numberofsegments = 0; //used on matlab
408+ IssmDouble* x = NULL;
409+ IssmDouble* y = NULL;
410+ IssmDouble* z = NULL;
411+ int* elements = NULL;
412+ int elementswidth = this->GetElementsWidth(); //just tria elements in this version. Itapopo:
413+ int levelmax = 0;
414+ IssmDouble regionlevel1 = 0.;
415+ IssmDouble regionlevelmax = 0.;
416
417 /*Get vertices coordinates of the coarse mesh (father mesh)*/
418 /*elements comes in Matlab indexing*/
419 this->GetMesh(this->vertices,this->elements,&x,&y,&z,&elements);
420
421+ /*Get amr parameters*/
422+ this->parameters->FindParam(&levelmax,AmrLevelMaxEnum);
423+ this->parameters->FindParam(&regionlevel1,AmrRegionLevel1Enum);
424+ this->parameters->FindParam(&regionlevelmax,AmrRegionLevelMaxEnum);
425+
426 /*Create initial mesh (coarse mesh) in neopz data structure*/
427 /*Just CPU #0 should keep AMR object*/
428- if(my_rank==0){
429- this->amr = new AdaptiveMeshRefinement();
430- int hmax = 0; //itapopo: it must be defined by interface. Using 2 just to test
431- this->amr->SetHMax(hmax); //Set max level of refinement
432- this->amr->CreateInitialMesh(numberofvertices,numberofelements,numberofsegments,elementswidth,x,y,z,elements,NULL);
433+ this->SetRefPatterns();
434+ if(my_rank==0){
435+ bool ismisomip = true;
436+ if(ismisomip){//itapopo
437+ TPZFileStream fstr;
438+ std::stringstream ss;
439+
440+ ss << levelmax;
441+ std::string AMRfile = "/home/santos/Misomip2/L" + ss.str() + "_tsai/amr.txt";
442+ fstr.OpenRead(AMRfile.c_str());
443+
444+ TPZSaveable *sv = TPZSaveable::Restore(fstr,0);
445+ this->amr = dynamic_cast<AdaptiveMeshRefinement*>(sv);
446+ }
447+ else{
448+ this->amr = new AdaptiveMeshRefinement();
449+ //this->amr->SetLevelMax(levelmax); //Set max level of refinement
450+ //this->amr->SetRegions(regionlevel1,regionlevelmax);
451+ this->amr->CreateInitialMesh(numberofvertices,numberofelements,numberofsegments,elementswidth,x,y,z,elements,NULL);
452+ }
453+ this->amr->SetLevelMax(levelmax); //Set max level of refinement
454+ this->amr->SetRegions(regionlevel1,regionlevelmax);
455 }
456
457 /*Free the vectors*/
458@@ -2952,6 +2977,40 @@
459
460 }
461 /*}}}*/
462+void FemModel::SetRefPatterns(){/*{{{*/
463+
464+ /*Initialize the global variable of refinement patterns*/
465+ gRefDBase.InitializeUniformRefPattern(EOned);
466+ gRefDBase.InitializeUniformRefPattern(ETriangle);
467+
468+ //gRefDBase.InitializeRefPatterns();
469+ /*Insert specifics patterns to ISSM core*/
470+ std::string filepath = REFPATTERNDIR;
471+ std::string filename1 = filepath + "/2D_Triang_Rib_3.rpt";
472+ std::string filename2 = filepath + "/2D_Triang_Rib_4.rpt";
473+ std::string filename3 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_4.rpt";
474+ std::string filename4 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_4_permuted.rpt";
475+ std::string filename5 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_5.rpt";
476+ std::string filename6 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_5_permuted.rpt";
477+ std::string filename7 = filepath + "/2D_Triang_Rib_5.rpt";
478+
479+ TPZAutoPointer<TPZRefPattern> refpat1 = new TPZRefPattern(filename1);
480+ TPZAutoPointer<TPZRefPattern> refpat2 = new TPZRefPattern(filename2);
481+ TPZAutoPointer<TPZRefPattern> refpat3 = new TPZRefPattern(filename3);
482+ TPZAutoPointer<TPZRefPattern> refpat4 = new TPZRefPattern(filename4);
483+ TPZAutoPointer<TPZRefPattern> refpat5 = new TPZRefPattern(filename5);
484+ TPZAutoPointer<TPZRefPattern> refpat6 = new TPZRefPattern(filename6);
485+ TPZAutoPointer<TPZRefPattern> refpat7 = new TPZRefPattern(filename7);
486+
487+ if(!gRefDBase.FindRefPattern(refpat1)) gRefDBase.InsertRefPattern(refpat1);
488+ if(!gRefDBase.FindRefPattern(refpat2)) gRefDBase.InsertRefPattern(refpat2);
489+ if(!gRefDBase.FindRefPattern(refpat3)) gRefDBase.InsertRefPattern(refpat3);
490+ if(!gRefDBase.FindRefPattern(refpat4)) gRefDBase.InsertRefPattern(refpat4);
491+ if(!gRefDBase.FindRefPattern(refpat5)) gRefDBase.InsertRefPattern(refpat5);
492+ if(!gRefDBase.FindRefPattern(refpat6)) gRefDBase.InsertRefPattern(refpat6);
493+ if(!gRefDBase.FindRefPattern(refpat7)) gRefDBase.InsertRefPattern(refpat7);
494+}
495+/*}}}*/
496 void FemModel::ReMesh(void){/*{{{*/
497
498 /*Variables*/
499@@ -3066,6 +3125,12 @@
500
501 GetMaskOfIceVerticesLSMx(this);
502
503+ /*Insert MISMIP+ bed topography*/
504+ if(true) this->BedrockFromMismipPlus();
505+
506+ /*Adjust base, thickness and mask grounded ice leve set*/
507+ if(true) this->AdjustBaseThicknessAndMask();
508+
509 /*Reset current configuration: */
510 analysis_type=this->analysis_type_list[this->analysis_counter];
511 SetCurrentConfiguration(analysis_type);
512@@ -3082,6 +3147,138 @@
513
514 }
515 /*}}}*/
516+void FemModel::BedrockFromMismipPlus(void){/*{{{*/
517+
518+ /*Insert bedrock from mismip+ setup*/
519+ /*This was used to Misomip project/simulations*/
520+
521+ IssmDouble x,y,bx,by;
522+ int numvertices = this->GetElementsWidth();
523+ IssmDouble* xyz_list = NULL;
524+ IssmDouble* r = xNew<IssmDouble>(numvertices);
525+
526+ for(int el=0;el<this->elements->Size();el++){
527+ Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(el));
528+
529+ element->GetVerticesCoordinates(&xyz_list);
530+ for(int i=0;i<numvertices;i++){
531+ x = *(xyz_list+3*i+0);
532+ y = *(xyz_list+3*i+1);
533+
534+ bx=-150.-728.8*std::pow(x/300000.,2)+343.91*std::pow(x/300000.,4)-50.57*std::pow(x/300000.,6);
535+ by=500./(1.+std::exp((-2./4000.)*(y-80000./2.-24000.)))+500./(1.+std::exp((2./4000.)*(y-80000./2.+24000.)));
536+
537+ r[i]=std::max(bx+by,-720.);
538+ }
539+
540+ /*insert new bedrock*/
541+ element->AddInput(BedEnum,&r[0],P1Enum);
542+
543+ /*Cleanup*/
544+ xDelete<IssmDouble>(xyz_list);
545+ }
546+
547+ /*Delete*/
548+ xDelete<IssmDouble>(r);
549+
550+ return;
551+}
552+/*}}}*/
553+void FemModel::AdjustBaseThicknessAndMask(void){/*{{{*/
554+
555+ int numvertices = this->GetElementsWidth();
556+ IssmDouble rho_water,rho_ice,density,base_float;
557+ IssmDouble* phi = xNew<IssmDouble>(numvertices);
558+ IssmDouble* h = xNew<IssmDouble>(numvertices);
559+ IssmDouble* s = xNew<IssmDouble>(numvertices);
560+ IssmDouble* b = xNew<IssmDouble>(numvertices);
561+ IssmDouble* r = xNew<IssmDouble>(numvertices);
562+ IssmDouble* sl = xNew<IssmDouble>(numvertices);
563+
564+ for(int el=0;el<this->elements->Size();el++){
565+ Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(el));
566+
567+ element->GetInputListOnVertices(&s[0],SurfaceEnum);
568+ element->GetInputListOnVertices(&r[0],BedEnum);
569+ element->GetInputListOnVertices(&sl[0],SealevelEnum);
570+ rho_water = element->matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum);
571+ rho_ice = element->matpar->GetMaterialParameter(MaterialsRhoIceEnum);
572+ density = rho_ice/rho_water;
573+
574+ for(int i=0;i<numvertices;i++){
575+ /*calculate base floatation (which supports given surface*/
576+ base_float = rho_ice*s[i]/(rho_ice-rho_water);
577+ if(r[i]>base_float){
578+ b[i] = r[i];
579+ }
580+ else {
581+ b[i] = base_float;
582+ }
583+
584+ if(abs(sl[i])>0) _error_("Sea level value not supported!");
585+ /*update thickness and mask grounded ice level set*/
586+ h[i] = s[i]-b[i];
587+ phi[i] = h[i]+r[i]/density;
588+ }
589+
590+ /*Update inputs*/
591+ element->AddInput(MaskGroundediceLevelsetEnum,&phi[0],P1Enum);
592+ element->AddInput(ThicknessEnum,&h[0],P1Enum);
593+ element->AddInput(BaseEnum,&b[0],P1Enum);
594+
595+ }
596+
597+ /*Delete*/
598+ xDelete<IssmDouble>(phi);
599+ xDelete<IssmDouble>(h);
600+ xDelete<IssmDouble>(s);
601+ xDelete<IssmDouble>(b);
602+ xDelete<IssmDouble>(r);
603+ xDelete<IssmDouble>(sl);
604+
605+ return;
606+}
607+/*}}}*/
608+void FemModel::WriteMeshInResults(void){/*{{{*/
609+
610+ int step = -1;
611+ int numberofelements = -1;
612+ int numberofvertices = -1;
613+ IssmDouble time = -1;
614+ IssmDouble* x = NULL;
615+ IssmDouble* y = NULL;
616+ IssmDouble* z = NULL;
617+ int* elementslist = NULL;
618+
619+ if(!this->elements || !this->vertices || !this->results || !this->parameters) return;
620+
621+ parameters->FindParam(&step,StepEnum);
622+ parameters->FindParam(&time,TimeEnum);
623+ numberofelements=this->elements->NumberOfElements();
624+ numberofvertices=this->vertices->NumberOfVertices();
625+
626+ /*Get mesh. Elementslist comes in Matlab indexing*/
627+ this->GetMesh(this->vertices,this->elements,&x,&y,&z,&elementslist);
628+
629+ /*Write mesh in Results*/
630+ this->results->AddResult(new GenericExternalResult<int*>(this->results->Size()+1,MeshElementsEnum,
631+ elementslist,numberofelements,this->GetElementsWidth(),step,time));
632+
633+ this->results->AddResult(new GenericExternalResult<IssmDouble*>(this->results->Size()+1,MeshXEnum,
634+ x,numberofvertices,1,step,time));
635+
636+ this->results->AddResult(new GenericExternalResult<IssmDouble*>(this->results->Size()+1,MeshYEnum,
637+ y,numberofvertices,1,step,time));
638+
639+ /*Cleanup*/
640+ xDelete<IssmDouble>(x);
641+ xDelete<IssmDouble>(y);
642+ xDelete<IssmDouble>(z);
643+ xDelete<int>(elementslist);
644+
645+ return;
646+}
647+/*}}}*/
648 void FemModel::InterpolateInputs(Vertices* newfemmodel_vertices,Elements* newfemmodel_elements){/*{{{*/
649
650 int maxinputs = MaximumNumberOfDefinitionsEnum;
651@@ -3215,7 +3412,7 @@
652 InterpFromMeshToMesh2dx(&P1inputsnew,Indexold,Xold,Yold,numverticesold,numelementsold,
653 P1inputsold,numverticesold,numP1inputs,
654 Xnew,Ynew,numverticesnew,NULL);
655-
656+
657 /*Insert P0 inputs into the new elements.*/
658 vector=NULL;
659 for(int i=0;i<numP0inputs;i++){
660@@ -3250,7 +3447,7 @@
661 vector=NULL;
662 for(int i=0;i<numP1inputs;i++){
663
664- /*Get P0 input vector from the interpolated matrix*/
665+ /*Get P1 input vector from the interpolated matrix*/
666 vector=xNew<IssmDouble>(numverticesnew);
667 for(int j=0;j<numverticesnew;j++) vector[j]=P1inputsnew[j*numP1inputs+i];//vector has all vertices (serial)
668
669Index: ../trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h
670===================================================================
671--- ../trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h (revision 21671)
672+++ ../trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h (revision 21672)
673@@ -57,17 +57,21 @@
674 /*General methods*/
675 void CleanUp(); // Clean all attributes
676 void Initialize(); // Initialize the attributes with NULL and values out of usually range
677- void SetHMax(int &h); // Define the max level of refinement
678- void SetElementWidth(int &width); // Define elements width
679+ void SetLevelMax(int &h); // Define the max level of refinement
680+ void SetRegions(double &D1,double Dhmax); // Define the regions which will be refined
681+ void SetElementWidth(int &width); // Define elements width
682 void 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=NULL); // A new mesh will be created and refined. This returns the new mesh
683 void CreateInitialMesh(int &nvertices,int &nelements,int &nsegments,int &width,double* x,double* y,double* z,int* elements,int* segments=NULL); // Create a NeoPZ geometric mesh by coords and elements
684- void CheckMesh(int &nvertices,int &nelements,int &nsegments,int &width,double** px,double** py,double** pz,int** pelements,int** psegments=NULL); // Check the consistency of the mesh
685+ TPZGeoMesh* CreateRefPatternMesh(TPZGeoMesh* gmesh);
686+ void CheckMesh(int &nvertices,int &nelements,int &nsegments,int &width,double** px,double** py,double** pz,int** pelements,int** psegments=NULL); // Check the consistency of the mesh
687
688 private:
689
690 /*Private attributes*/
691 int elementswidth; // Geometric nodes for element: 3 == Tria, 4 == Tetra, 6 == Penta
692- int hmax; // Max level of refinement
693+ int levelmax; // Max level of refinement
694+ double regionlevel1; // Region which will be refined with level 1
695+ double regionlevelmax; // Region which will be refined with level max
696 TPZGeoMesh *fathermesh; // Father Mesh is the entire mesh without refinement
697 TPZGeoMesh *previousmesh; // Previous mesh is a refined mesh of last step
698
699@@ -82,7 +86,6 @@
700 void GetMesh(TPZGeoMesh *gmesh,int &nvertices,int &nelements,int &nsegments,double** px,double** py,double** pz,int** pelements,int** psegments=NULL); // Return coords and elements in ISSM data structure
701 inline int GetElemMaterialID(){return 1;} // Return element material ID
702 inline int GetBoundaryMaterialID(){return 2;} // Return segment (2D boundary) material ID
703- void SetRefPatterns(); // Initialize and insert the refinement patterns
704 };
705
706 #endif
707Index: ../trunk-jpl/src/c/classes/FemModel.h
708===================================================================
709--- ../trunk-jpl/src/c/classes/FemModel.h (revision 21671)
710+++ ../trunk-jpl/src/c/classes/FemModel.h (revision 21672)
711@@ -148,6 +148,8 @@
712 /*Adaptive mesh refinement methods*/
713 void InitializeAdaptiveRefinement(void);
714 void ReMesh(void);
715+ void BedrockFromMismipPlus(void);
716+ void AdjustBaseThicknessAndMask(void);
717 void GetMesh(Vertices* femmodel_vertices,Elements* femmodel_elements,IssmDouble** px, IssmDouble** py, IssmDouble** pz, int** pelementslist);
718 int GetElementsWidth(){return 3;};//just tria elements in this first version
719 void ExecuteRefinement(int &numberofvertices,int &numberofelements,IssmDouble** px,IssmDouble** py,IssmDouble** pz,int** pelementslist);
720@@ -160,6 +162,8 @@
721 void InterpolateInputs(Vertices* newfemmodel_vertices,Elements* newfemmodel_elements);
722 void UpdateElements(int newnumberofelements,int* newelementslist,bool* my_elements,int nodecounter,int analysis_counter,Elements* newelements);
723 void ElementsAndVerticesPartitioning(int& newnumberofvertices,int& newnumberofelements,int& elementswidth,int* newelementslist,bool** pmy_elements,int** pmy_vertices);
724+ void WriteMeshInResults(void);
725+ void SetRefPatterns(void);
726 #endif
727 };
728
Note: See TracBrowser for help on using the repository browser.