![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 21672
- Timestamp:
- 04/18/17 07:33:24 (8 years ago)
- Location:
- issm/trunk-jpl/src/c/classes
- Files:
-
- 4 edited
-
AdaptiveMeshRefinement.cpp (modified) (19 diffs)
-
AdaptiveMeshRefinement.h (modified) (3 diffs)
-
FemModel.cpp (modified) (7 diffs)
-
FemModel.h (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp
r21573 r21672 27 27 28 28 /*Copy all data*/ 29 this->fathermesh = new TPZGeoMesh(*cp.fathermesh); 30 this->previousmesh = new TPZGeoMesh(*cp.previousmesh); 31 this->hmax = cp.hmax; 32 this->elementswidth = cp.elementswidth; 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; 33 35 return *this; 34 36 … … 36 38 /*}}}*/ 37 39 AdaptiveMeshRefinement::~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 } 38 53 this->CleanUp(); 39 54 gRefDBase.clear(); … … 45 60 if(this->fathermesh) delete this->fathermesh; 46 61 if(this->previousmesh) delete this->previousmesh; 47 this->hmax=-1; 48 this->elementswidth=-1; 62 this->levelmax = -1; 63 this->elementswidth = -1; 64 this->regionlevel1 = -1; 65 this->regionlevelmax = -1; 49 66 } 50 67 /*}}}*/ … … 52 69 53 70 /*Set pointers to NULL*/ 54 this->fathermesh = NULL; 55 this->previousmesh = NULL; 56 this->hmax = -1; 57 this->elementswidth = -1; 71 this->fathermesh = NULL; 72 this->previousmesh = NULL; 73 this->levelmax = -1; 74 this->elementswidth = -1; 75 this->regionlevel1 = -1; 76 this->regionlevelmax = -1; 58 77 } 59 78 /*}}}*/ … … 82 101 83 102 /* Read simple attributes */ 84 buf.Read(&this-> hmax,1);103 buf.Read(&this->levelmax,1); 85 104 buf.Read(&this->elementswidth,1); 86 87 /* Read geometric mesh*/ 105 buf.Read(&this->regionlevel1,1); 106 buf.Read(&this->regionlevelmax,1); 107 108 /* Read geometric mesh*/ 88 109 TPZSaveable *sv1 = TPZSaveable::Restore(buf,0); 89 110 this->fathermesh = dynamic_cast<TPZGeoMesh*>(sv1); … … 114 135 115 136 /* Write simple attributes */ 116 buf.Write(&this-> hmax,1);137 buf.Write(&this->levelmax,1); 117 138 buf.Write(&this->elementswidth,1); 118 139 buf.Write(&this->regionlevel1,1); 140 buf.Write(&this->regionlevelmax,1); 141 119 142 /* Write the geometric mesh*/ 120 143 this->fathermesh->Write(buf, this->ClassId()); … … 138 161 /*NEOPZ works only in C indexing*/ 139 162 140 //itapopo_assert_(this->fathermesh);141 //itapopo_assert_(this->previousmesh);163 _assert_(this->fathermesh); 164 _assert_(this->previousmesh); 142 165 143 166 /*Calculate the position of the grounding line using previous mesh*/ … … 145 168 this->CalcGroundingLinePosition(masklevelset, GLvec); 146 169 147 std::ofstream file1("/ronne_1/home/santos/mesh0.vtk");148 TPZVTKGeoMesh::PrintGMeshVTK(this->previousmesh,file1 );170 // std::ofstream file1("/home/santos/mesh0.vtk"); 171 // TPZVTKGeoMesh::PrintGMeshVTK(this->fathermesh,file1 ); 149 172 150 173 /*run refinement or unrefinement process*/ … … 158 181 this->RefinementProcess(newmesh,GLvec); 159 182 160 std::ofstream file2("/ronne_1/home/santos/mesh1.vtk");161 TPZVTKGeoMesh::PrintGMeshVTK(this->previousmesh,file2 );183 //std::ofstream file2("/home/santos/mesh1.vtk"); 184 //TPZVTKGeoMesh::PrintGMeshVTK(this->previousmesh,file2 ); 162 185 163 186 /*Set new mesh pointer. Previous mesh just have uniform elements*/ … … 168 191 169 192 /*Refine elements to avoid hanging nodes*/ 170 TPZGeoMesh *nohangingnodesmesh = new TPZGeoMesh(*newmesh); 171 172 std::ofstream file3("/ronne_1/home/santos/mesh2.vtk"); 173 TPZVTKGeoMesh::PrintGMeshVTK(newmesh,file3); 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); 174 198 175 199 this->RefineMeshToAvoidHangingNodes(nohangingnodesmesh); 176 200 177 std::ofstream file4("/ronne_1/home/santos/mesh3.vtk");178 TPZVTKGeoMesh::PrintGMeshVTK(nohangingnodesmesh,file4);201 //std::ofstream file4("/home/santos/mesh3.vtk"); 202 //TPZVTKGeoMesh::PrintGMeshVTK(nohangingnodesmesh,file4); 179 203 180 204 /*Get new geometric mesh in ISSM data structure*/ … … 191 215 void AdaptiveMeshRefinement::RefinementProcess(TPZGeoMesh *gmesh,std::vector<TPZVec<REAL> > &GLvec){/*{{{*/ 192 216 193 /*Refine mesh hmax times*/194 _printf_("\n\trefinement process ( hmax = " << this->hmax << ")\n");217 /*Refine mesh levelmax times*/ 218 _printf_("\n\trefinement process (level max = " << this->levelmax << ")\n"); 195 219 _printf_("\tprogress: "); 196 for(int hlevel=1;hlevel<=this-> hmax;hlevel++){220 for(int hlevel=1;hlevel<=this->levelmax;hlevel++){ 197 221 198 222 /*Set elements to be refined using some criteria*/ … … 365 389 int vertex2 = this->previousmesh->Element(i)->NodeIndex(2); 366 390 391 //itapopo inserir uma verificação para não acessar fora da memória 367 392 double mls0 = masklevelset[vertex0]; 368 393 double mls1 = masklevelset[vertex1]; … … 419 444 void AdaptiveMeshRefinement::TagElementsNearGroundingLine(TPZGeoMesh *gmesh,std::vector<TPZVec<REAL> > &GLvec,int &hlevel,std::vector<int> &ElemVec){/*{{{*/ 420 445 421 /* Tag elements near grounding line */ 422 double MaxRegion = 40000.; //itapopo 423 double alpha = log(1.5); //itapopo 424 double MaxDistance = MaxRegion / std::exp(alpha*(hlevel-1)); 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)); 425 452 426 453 for(int i=0;i<gmesh->NElements();i++){ … … 436 463 gmesh->Element(i)->X(qsi, centerPoint); 437 464 438 REAL distance = MaxDistance;465 REAL distance = Di; 439 466 440 467 for (int j = 0; j < GLvec.size(); j++) { … … 442 469 REAL value = ( GLvec[j][0] - centerPoint[0] ) * ( GLvec[j][0] - centerPoint[0] ); // (x2-x1)^2 443 470 value += ( GLvec[j][1] - centerPoint[1] ) * ( GLvec[j][1] - centerPoint[1] );// (y2-y1)^2 444 value = std::sqrt(value); // /Radius471 value = std::sqrt(value); //Radius 445 472 446 473 //finding the min distance to the grounding line … … 449 476 } 450 477 451 if(distance < MaxDistance) ElemVec.push_back(i);478 if(distance < Di) ElemVec.push_back(i); 452 479 } 453 480 … … 459 486 /*NEOPZ works only in C indexing*/ 460 487 461 // itapopo_assert_(nvertices>0);462 // itapoo_assert_(nelements>0);488 _assert_(nvertices>0); 489 _assert_(nelements>0); 463 490 this->SetElementWidth(width); 464 491 465 492 /*Verify and creating initial mesh*/ 466 //itapopoif(this->fathermesh) _error_("Initial mesh already exists!");493 if(this->fathermesh) _error_("Initial mesh already exists!"); 467 494 468 495 this->fathermesh = new TPZGeoMesh(); … … 491 518 492 519 /*reftype = 0: uniform, fast / reftype = 1: uniform and non-uniform (avoid hanging nodes), it is not too fast */ 493 const int reftype = 1;520 const int reftype = 0; 494 521 switch(this->elementswidth){ 495 522 case 3: this->fathermesh->CreateGeoElement(ETriangle, elem, mat, index, reftype); break; … … 521 548 /*Create previous mesh as a copy of father mesh*/ 522 549 this->previousmesh = new TPZGeoMesh(*this->fathermesh); 523 524 /*Set refinement patterns*/ 525 this->SetRefPatterns(); 526 527 } 528 /*}}}*/ 529 void AdaptiveMeshRefinement::SetHMax(int &h){/*{{{*/ 530 this->hmax = h; 550 551 } 552 /*}}}*/ 553 #include "pzgeotriangle.h" //itapopo 554 #include "pzreftriangle.h" //itapopo 555 using namespace pzgeom; 556 TPZGeoMesh* 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 /*}}}*/ 638 void AdaptiveMeshRefinement::SetLevelMax(int &h){/*{{{*/ 639 this->levelmax = h; 640 } 641 /*}}}*/ 642 void AdaptiveMeshRefinement::SetRegions(double &D1,double Dhmax){/*{{{*/ 643 this->regionlevel1 = D1; 644 this->regionlevelmax = Dhmax; 531 645 } 532 646 /*}}}*/ 533 647 void AdaptiveMeshRefinement::SetElementWidth(int &width){/*{{{*/ 534 648 this->elementswidth = width; 535 }536 /*}}}*/537 void AdaptiveMeshRefinement::SetRefPatterns(){/*{{{*/538 539 /*Initialize the global variable of refinement patterns*/540 gRefDBase.InitializeUniformRefPattern(EOned);541 gRefDBase.InitializeUniformRefPattern(ETriangle);542 543 //gRefDBase.InitializeRefPatterns();544 /*Insert specifics patterns to ISSM core*/545 std::string filepath = REFPATTERNDIR;546 std::string filename1 = filepath + "/2D_Triang_Rib_3.rpt";547 std::string filename2 = filepath + "/2D_Triang_Rib_4.rpt";548 std::string filename3 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_4.rpt";549 std::string filename4 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_4_permuted.rpt";550 std::string filename5 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_5.rpt";551 std::string filename6 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_5_permuted.rpt";552 std::string filename7 = filepath + "/2D_Triang_Rib_5.rpt";553 554 TPZAutoPointer<TPZRefPattern> refpat1 = new TPZRefPattern(filename1);555 TPZAutoPointer<TPZRefPattern> refpat2 = new TPZRefPattern(filename2);556 TPZAutoPointer<TPZRefPattern> refpat3 = new TPZRefPattern(filename3);557 TPZAutoPointer<TPZRefPattern> refpat4 = new TPZRefPattern(filename4);558 TPZAutoPointer<TPZRefPattern> refpat5 = new TPZRefPattern(filename5);559 TPZAutoPointer<TPZRefPattern> refpat6 = new TPZRefPattern(filename6);560 TPZAutoPointer<TPZRefPattern> refpat7 = new TPZRefPattern(filename7);561 562 if(!gRefDBase.FindRefPattern(refpat1)) gRefDBase.InsertRefPattern(refpat1);563 if(!gRefDBase.FindRefPattern(refpat2)) gRefDBase.InsertRefPattern(refpat2);564 if(!gRefDBase.FindRefPattern(refpat3)) gRefDBase.InsertRefPattern(refpat3);565 if(!gRefDBase.FindRefPattern(refpat4)) gRefDBase.InsertRefPattern(refpat4);566 if(!gRefDBase.FindRefPattern(refpat5)) gRefDBase.InsertRefPattern(refpat5);567 if(!gRefDBase.FindRefPattern(refpat6)) gRefDBase.InsertRefPattern(refpat6);568 if(!gRefDBase.FindRefPattern(refpat7)) gRefDBase.InsertRefPattern(refpat7);569 649 } 570 650 /*}}}*/ -
issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h
r21562 r21672 58 58 void CleanUp(); // Clean all attributes 59 59 void Initialize(); // Initialize the attributes with NULL and values out of usually range 60 void SetHMax(int &h); // Define the max level of refinement 61 void SetElementWidth(int &width); // Define elements width 60 void SetLevelMax(int &h); // Define the max level of refinement 61 void SetRegions(double &D1,double Dhmax); // Define the regions which will be refined 62 void SetElementWidth(int &width); // Define elements width 62 63 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 63 64 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 64 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 65 TPZGeoMesh* CreateRefPatternMesh(TPZGeoMesh* gmesh); 66 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 65 67 66 68 private: … … 68 70 /*Private attributes*/ 69 71 int elementswidth; // Geometric nodes for element: 3 == Tria, 4 == Tetra, 6 == Penta 70 int hmax; // Max level of refinement 72 int levelmax; // Max level of refinement 73 double regionlevel1; // Region which will be refined with level 1 74 double regionlevelmax; // Region which will be refined with level max 71 75 TPZGeoMesh *fathermesh; // Father Mesh is the entire mesh without refinement 72 76 TPZGeoMesh *previousmesh; // Previous mesh is a refined mesh of last step … … 83 87 inline int GetElemMaterialID(){return 1;} // Return element material ID 84 88 inline int GetBoundaryMaterialID(){return 2;} // Return segment (2D boundary) material ID 85 void SetRefPatterns(); // Initialize and insert the refinement patterns86 89 }; 87 90 -
issm/trunk-jpl/src/c/classes/FemModel.cpp
r21583 r21672 2921 2921 2922 2922 /*Define variables*/ 2923 int my_rank = IssmComm::GetRank(); 2924 this->amr = NULL;//initialize amr as NULL 2925 int numberofvertices = this->vertices->NumberOfVertices(); 2926 int numberofelements = this->elements->NumberOfElements(); 2927 int numberofsegments = 0; //used on matlab 2928 IssmDouble* x = NULL; 2929 IssmDouble* y = NULL; 2930 IssmDouble* z = NULL; 2931 int* elements = NULL; 2932 int elementswidth = this->GetElementsWidth(); //just tria elements in this version. Itapopo: 2923 int my_rank = IssmComm::GetRank(); 2924 this->amr = NULL;//initialize amr as NULL 2925 int numberofvertices = this->vertices->NumberOfVertices(); 2926 int numberofelements = this->elements->NumberOfElements(); 2927 int numberofsegments = 0; //used on matlab 2928 IssmDouble* x = NULL; 2929 IssmDouble* y = NULL; 2930 IssmDouble* z = NULL; 2931 int* elements = NULL; 2932 int elementswidth = this->GetElementsWidth(); //just tria elements in this version. Itapopo: 2933 int levelmax = 0; 2934 IssmDouble regionlevel1 = 0.; 2935 IssmDouble regionlevelmax = 0.; 2933 2936 2934 2937 /*Get vertices coordinates of the coarse mesh (father mesh)*/ … … 2936 2939 this->GetMesh(this->vertices,this->elements,&x,&y,&z,&elements); 2937 2940 2941 /*Get amr parameters*/ 2942 this->parameters->FindParam(&levelmax,AmrLevelMaxEnum); 2943 this->parameters->FindParam(®ionlevel1,AmrRegionLevel1Enum); 2944 this->parameters->FindParam(®ionlevelmax,AmrRegionLevelMaxEnum); 2945 2938 2946 /*Create initial mesh (coarse mesh) in neopz data structure*/ 2939 2947 /*Just CPU #0 should keep AMR object*/ 2940 if(my_rank==0){ 2941 this->amr = new AdaptiveMeshRefinement(); 2942 int hmax = 0; //itapopo: it must be defined by interface. Using 2 just to test 2943 this->amr->SetHMax(hmax); //Set max level of refinement 2944 this->amr->CreateInitialMesh(numberofvertices,numberofelements,numberofsegments,elementswidth,x,y,z,elements,NULL); 2948 this->SetRefPatterns(); 2949 if(my_rank==0){ 2950 bool ismisomip = true; 2951 if(ismisomip){//itapopo 2952 TPZFileStream fstr; 2953 std::stringstream ss; 2954 2955 ss << levelmax; 2956 std::string AMRfile = "/home/santos/Misomip2/L" + ss.str() + "_tsai/amr.txt"; 2957 fstr.OpenRead(AMRfile.c_str()); 2958 2959 TPZSaveable *sv = TPZSaveable::Restore(fstr,0); 2960 this->amr = dynamic_cast<AdaptiveMeshRefinement*>(sv); 2961 } 2962 else{ 2963 this->amr = new AdaptiveMeshRefinement(); 2964 //this->amr->SetLevelMax(levelmax); //Set max level of refinement 2965 //this->amr->SetRegions(regionlevel1,regionlevelmax); 2966 this->amr->CreateInitialMesh(numberofvertices,numberofelements,numberofsegments,elementswidth,x,y,z,elements,NULL); 2967 } 2968 this->amr->SetLevelMax(levelmax); //Set max level of refinement 2969 this->amr->SetRegions(regionlevel1,regionlevelmax); 2945 2970 } 2946 2971 … … 2951 2976 xDelete<int>(elements); 2952 2977 2978 } 2979 /*}}}*/ 2980 void FemModel::SetRefPatterns(){/*{{{*/ 2981 2982 /*Initialize the global variable of refinement patterns*/ 2983 gRefDBase.InitializeUniformRefPattern(EOned); 2984 gRefDBase.InitializeUniformRefPattern(ETriangle); 2985 2986 //gRefDBase.InitializeRefPatterns(); 2987 /*Insert specifics patterns to ISSM core*/ 2988 std::string filepath = REFPATTERNDIR; 2989 std::string filename1 = filepath + "/2D_Triang_Rib_3.rpt"; 2990 std::string filename2 = filepath + "/2D_Triang_Rib_4.rpt"; 2991 std::string filename3 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_4.rpt"; 2992 std::string filename4 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_4_permuted.rpt"; 2993 std::string filename5 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_5.rpt"; 2994 std::string filename6 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_5_permuted.rpt"; 2995 std::string filename7 = filepath + "/2D_Triang_Rib_5.rpt"; 2996 2997 TPZAutoPointer<TPZRefPattern> refpat1 = new TPZRefPattern(filename1); 2998 TPZAutoPointer<TPZRefPattern> refpat2 = new TPZRefPattern(filename2); 2999 TPZAutoPointer<TPZRefPattern> refpat3 = new TPZRefPattern(filename3); 3000 TPZAutoPointer<TPZRefPattern> refpat4 = new TPZRefPattern(filename4); 3001 TPZAutoPointer<TPZRefPattern> refpat5 = new TPZRefPattern(filename5); 3002 TPZAutoPointer<TPZRefPattern> refpat6 = new TPZRefPattern(filename6); 3003 TPZAutoPointer<TPZRefPattern> refpat7 = new TPZRefPattern(filename7); 3004 3005 if(!gRefDBase.FindRefPattern(refpat1)) gRefDBase.InsertRefPattern(refpat1); 3006 if(!gRefDBase.FindRefPattern(refpat2)) gRefDBase.InsertRefPattern(refpat2); 3007 if(!gRefDBase.FindRefPattern(refpat3)) gRefDBase.InsertRefPattern(refpat3); 3008 if(!gRefDBase.FindRefPattern(refpat4)) gRefDBase.InsertRefPattern(refpat4); 3009 if(!gRefDBase.FindRefPattern(refpat5)) gRefDBase.InsertRefPattern(refpat5); 3010 if(!gRefDBase.FindRefPattern(refpat6)) gRefDBase.InsertRefPattern(refpat6); 3011 if(!gRefDBase.FindRefPattern(refpat7)) gRefDBase.InsertRefPattern(refpat7); 2953 3012 } 2954 3013 /*}}}*/ … … 3067 3126 GetMaskOfIceVerticesLSMx(this); 3068 3127 3128 /*Insert MISMIP+ bed topography*/ 3129 if(true) this->BedrockFromMismipPlus(); 3130 3131 /*Adjust base, thickness and mask grounded ice leve set*/ 3132 if(true) this->AdjustBaseThicknessAndMask(); 3133 3069 3134 /*Reset current configuration: */ 3070 3135 analysis_type=this->analysis_type_list[this->analysis_counter]; … … 3081 3146 return; 3082 3147 3148 } 3149 /*}}}*/ 3150 void FemModel::BedrockFromMismipPlus(void){/*{{{*/ 3151 3152 /*Insert bedrock from mismip+ setup*/ 3153 /*This was used to Misomip project/simulations*/ 3154 3155 IssmDouble x,y,bx,by; 3156 int numvertices = this->GetElementsWidth(); 3157 IssmDouble* xyz_list = NULL; 3158 IssmDouble* r = xNew<IssmDouble>(numvertices); 3159 3160 for(int el=0;el<this->elements->Size();el++){ 3161 Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(el)); 3162 3163 element->GetVerticesCoordinates(&xyz_list); 3164 for(int i=0;i<numvertices;i++){ 3165 x = *(xyz_list+3*i+0); 3166 y = *(xyz_list+3*i+1); 3167 3168 bx=-150.-728.8*std::pow(x/300000.,2)+343.91*std::pow(x/300000.,4)-50.57*std::pow(x/300000.,6); 3169 by=500./(1.+std::exp((-2./4000.)*(y-80000./2.-24000.)))+500./(1.+std::exp((2./4000.)*(y-80000./2.+24000.))); 3170 3171 r[i]=std::max(bx+by,-720.); 3172 } 3173 3174 /*insert new bedrock*/ 3175 element->AddInput(BedEnum,&r[0],P1Enum); 3176 3177 /*Cleanup*/ 3178 xDelete<IssmDouble>(xyz_list); 3179 } 3180 3181 /*Delete*/ 3182 xDelete<IssmDouble>(r); 3183 3184 return; 3185 } 3186 /*}}}*/ 3187 void FemModel::AdjustBaseThicknessAndMask(void){/*{{{*/ 3188 3189 int numvertices = this->GetElementsWidth(); 3190 IssmDouble rho_water,rho_ice,density,base_float; 3191 IssmDouble* phi = xNew<IssmDouble>(numvertices); 3192 IssmDouble* h = xNew<IssmDouble>(numvertices); 3193 IssmDouble* s = xNew<IssmDouble>(numvertices); 3194 IssmDouble* b = xNew<IssmDouble>(numvertices); 3195 IssmDouble* r = xNew<IssmDouble>(numvertices); 3196 IssmDouble* sl = xNew<IssmDouble>(numvertices); 3197 3198 for(int el=0;el<this->elements->Size();el++){ 3199 Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(el)); 3200 3201 element->GetInputListOnVertices(&s[0],SurfaceEnum); 3202 element->GetInputListOnVertices(&r[0],BedEnum); 3203 element->GetInputListOnVertices(&sl[0],SealevelEnum); 3204 rho_water = element->matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum); 3205 rho_ice = element->matpar->GetMaterialParameter(MaterialsRhoIceEnum); 3206 density = rho_ice/rho_water; 3207 3208 for(int i=0;i<numvertices;i++){ 3209 /*calculate base floatation (which supports given surface*/ 3210 base_float = rho_ice*s[i]/(rho_ice-rho_water); 3211 if(r[i]>base_float){ 3212 b[i] = r[i]; 3213 } 3214 else { 3215 b[i] = base_float; 3216 } 3217 3218 if(abs(sl[i])>0) _error_("Sea level value not supported!"); 3219 /*update thickness and mask grounded ice level set*/ 3220 h[i] = s[i]-b[i]; 3221 phi[i] = h[i]+r[i]/density; 3222 } 3223 3224 /*Update inputs*/ 3225 element->AddInput(MaskGroundediceLevelsetEnum,&phi[0],P1Enum); 3226 element->AddInput(ThicknessEnum,&h[0],P1Enum); 3227 element->AddInput(BaseEnum,&b[0],P1Enum); 3228 3229 } 3230 3231 /*Delete*/ 3232 xDelete<IssmDouble>(phi); 3233 xDelete<IssmDouble>(h); 3234 xDelete<IssmDouble>(s); 3235 xDelete<IssmDouble>(b); 3236 xDelete<IssmDouble>(r); 3237 xDelete<IssmDouble>(sl); 3238 3239 return; 3240 } 3241 /*}}}*/ 3242 void FemModel::WriteMeshInResults(void){/*{{{*/ 3243 3244 int step = -1; 3245 int numberofelements = -1; 3246 int numberofvertices = -1; 3247 IssmDouble time = -1; 3248 IssmDouble* x = NULL; 3249 IssmDouble* y = NULL; 3250 IssmDouble* z = NULL; 3251 int* elementslist = NULL; 3252 3253 if(!this->elements || !this->vertices || !this->results || !this->parameters) return; 3254 3255 parameters->FindParam(&step,StepEnum); 3256 parameters->FindParam(&time,TimeEnum); 3257 numberofelements=this->elements->NumberOfElements(); 3258 numberofvertices=this->vertices->NumberOfVertices(); 3259 3260 /*Get mesh. Elementslist comes in Matlab indexing*/ 3261 this->GetMesh(this->vertices,this->elements,&x,&y,&z,&elementslist); 3262 3263 /*Write mesh in Results*/ 3264 this->results->AddResult(new GenericExternalResult<int*>(this->results->Size()+1,MeshElementsEnum, 3265 elementslist,numberofelements,this->GetElementsWidth(),step,time)); 3266 3267 this->results->AddResult(new GenericExternalResult<IssmDouble*>(this->results->Size()+1,MeshXEnum, 3268 x,numberofvertices,1,step,time)); 3269 3270 this->results->AddResult(new GenericExternalResult<IssmDouble*>(this->results->Size()+1,MeshYEnum, 3271 y,numberofvertices,1,step,time)); 3272 3273 /*Cleanup*/ 3274 xDelete<IssmDouble>(x); 3275 xDelete<IssmDouble>(y); 3276 xDelete<IssmDouble>(z); 3277 xDelete<int>(elementslist); 3278 3279 return; 3083 3280 } 3084 3281 /*}}}*/ … … 3216 3413 P1inputsold,numverticesold,numP1inputs, 3217 3414 Xnew,Ynew,numverticesnew,NULL); 3218 3415 3219 3416 /*Insert P0 inputs into the new elements.*/ 3220 3417 vector=NULL; … … 3251 3448 for(int i=0;i<numP1inputs;i++){ 3252 3449 3253 /*Get P 0input vector from the interpolated matrix*/3450 /*Get P1 input vector from the interpolated matrix*/ 3254 3451 vector=xNew<IssmDouble>(numverticesnew); 3255 3452 for(int j=0;j<numverticesnew;j++) vector[j]=P1inputsnew[j*numP1inputs+i];//vector has all vertices (serial) -
issm/trunk-jpl/src/c/classes/FemModel.h
r21583 r21672 149 149 void InitializeAdaptiveRefinement(void); 150 150 void ReMesh(void); 151 void BedrockFromMismipPlus(void); 152 void AdjustBaseThicknessAndMask(void); 151 153 void GetMesh(Vertices* femmodel_vertices,Elements* femmodel_elements,IssmDouble** px, IssmDouble** py, IssmDouble** pz, int** pelementslist); 152 154 int GetElementsWidth(){return 3;};//just tria elements in this first version … … 161 163 void UpdateElements(int newnumberofelements,int* newelementslist,bool* my_elements,int nodecounter,int analysis_counter,Elements* newelements); 162 164 void ElementsAndVerticesPartitioning(int& newnumberofvertices,int& newnumberofelements,int& elementswidth,int* newelementslist,bool** pmy_elements,int** pmy_vertices); 165 void WriteMeshInResults(void); 166 void SetRefPatterns(void); 163 167 #endif 164 168 };
Note:
See TracChangeset
for help on using the changeset viewer.
