Index: /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp (revision 21671) +++ /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp (revision 21672) @@ -27,8 +27,10 @@ /*Copy all data*/ - this->fathermesh = new TPZGeoMesh(*cp.fathermesh); - this->previousmesh = new TPZGeoMesh(*cp.previousmesh); - this->hmax = cp.hmax; - this->elementswidth = cp.elementswidth; + this->fathermesh = new TPZGeoMesh(*cp.fathermesh); + this->previousmesh = new TPZGeoMesh(*cp.previousmesh); + this->levelmax = cp.levelmax; + this->elementswidth = cp.elementswidth; + this->regionlevel1 = cp.regionlevel1; + this->regionlevelmax = cp.regionlevelmax; return *this; @@ -36,4 +38,17 @@ /*}}}*/ AdaptiveMeshRefinement::~AdaptiveMeshRefinement(){/*{{{*/ + + bool ismismip = true; + if(ismismip){//itapopo + TPZFileStream fstr; + std::stringstream ss; + + ss << this->levelmax; + std::string AMRfile = "/home/santos/Misomip2/L" + ss.str() + "_tsai/amr.txt"; + + fstr.OpenWrite(AMRfile.c_str()); + int withclassid = 1; + this->Write(fstr,withclassid); + } this->CleanUp(); gRefDBase.clear(); @@ -45,6 +60,8 @@ if(this->fathermesh) delete this->fathermesh; if(this->previousmesh) delete this->previousmesh; - this->hmax=-1; - this->elementswidth=-1; + this->levelmax = -1; + this->elementswidth = -1; + this->regionlevel1 = -1; + this->regionlevelmax = -1; } /*}}}*/ @@ -52,8 +69,10 @@ /*Set pointers to NULL*/ - this->fathermesh = NULL; - this->previousmesh = NULL; - this->hmax = -1; - this->elementswidth = -1; + this->fathermesh = NULL; + this->previousmesh = NULL; + this->levelmax = -1; + this->elementswidth = -1; + this->regionlevel1 = -1; + this->regionlevelmax = -1; } /*}}}*/ @@ -82,8 +101,10 @@ /* Read simple attributes */ - buf.Read(&this->hmax,1); + buf.Read(&this->levelmax,1); buf.Read(&this->elementswidth,1); - - /* Read geometric mesh*/ + buf.Read(&this->regionlevel1,1); + buf.Read(&this->regionlevelmax,1); + + /* Read geometric mesh*/ TPZSaveable *sv1 = TPZSaveable::Restore(buf,0); this->fathermesh = dynamic_cast(sv1); @@ -114,7 +135,9 @@ /* Write simple attributes */ - buf.Write(&this->hmax,1); + buf.Write(&this->levelmax,1); buf.Write(&this->elementswidth,1); - + buf.Write(&this->regionlevel1,1); + buf.Write(&this->regionlevelmax,1); + /* Write the geometric mesh*/ this->fathermesh->Write(buf, this->ClassId()); @@ -138,6 +161,6 @@ /*NEOPZ works only in C indexing*/ - //itapopo _assert_(this->fathermesh); - //itapopo _assert_(this->previousmesh); + _assert_(this->fathermesh); + _assert_(this->previousmesh); /*Calculate the position of the grounding line using previous mesh*/ @@ -145,6 +168,6 @@ this->CalcGroundingLinePosition(masklevelset, GLvec); - std::ofstream file1("/ronne_1/home/santos/mesh0.vtk"); - TPZVTKGeoMesh::PrintGMeshVTK(this->previousmesh,file1 ); + // std::ofstream file1("/home/santos/mesh0.vtk"); + // TPZVTKGeoMesh::PrintGMeshVTK(this->fathermesh,file1 ); /*run refinement or unrefinement process*/ @@ -158,6 +181,6 @@ this->RefinementProcess(newmesh,GLvec); - std::ofstream file2("/ronne_1/home/santos/mesh1.vtk"); - TPZVTKGeoMesh::PrintGMeshVTK(this->previousmesh,file2 ); + //std::ofstream file2("/home/santos/mesh1.vtk"); + //TPZVTKGeoMesh::PrintGMeshVTK(this->previousmesh,file2 ); /*Set new mesh pointer. Previous mesh just have uniform elements*/ @@ -168,13 +191,14 @@ /*Refine elements to avoid hanging nodes*/ - TPZGeoMesh *nohangingnodesmesh = new TPZGeoMesh(*newmesh); - - std::ofstream file3("/ronne_1/home/santos/mesh2.vtk"); - TPZVTKGeoMesh::PrintGMeshVTK(newmesh,file3); + //TPZGeoMesh *nohangingnodesmesh = new TPZGeoMesh(*newmesh);//itapopo testando, este era o original + TPZGeoMesh *nohangingnodesmesh = this->CreateRefPatternMesh(newmesh);//itapopo testando, este eh novo metodo + + //std::ofstream file3("/home/santos/mesh2.vtk"); + //TPZVTKGeoMesh::PrintGMeshVTK(this->previousmesh,file3); this->RefineMeshToAvoidHangingNodes(nohangingnodesmesh); - std::ofstream file4("/ronne_1/home/santos/mesh3.vtk"); - TPZVTKGeoMesh::PrintGMeshVTK(nohangingnodesmesh,file4); + //std::ofstream file4("/home/santos/mesh3.vtk"); + //TPZVTKGeoMesh::PrintGMeshVTK(nohangingnodesmesh,file4); /*Get new geometric mesh in ISSM data structure*/ @@ -191,8 +215,8 @@ void AdaptiveMeshRefinement::RefinementProcess(TPZGeoMesh *gmesh,std::vector > &GLvec){/*{{{*/ - /*Refine mesh hmax times*/ - _printf_("\n\trefinement process (hmax = " << this->hmax << ")\n"); + /*Refine mesh levelmax times*/ + _printf_("\n\trefinement process (level max = " << this->levelmax << ")\n"); _printf_("\tprogress: "); - for(int hlevel=1;hlevel<=this->hmax;hlevel++){ + for(int hlevel=1;hlevel<=this->levelmax;hlevel++){ /*Set elements to be refined using some criteria*/ @@ -365,4 +389,5 @@ int vertex2 = this->previousmesh->Element(i)->NodeIndex(2); + //itapopo inserir uma verificação para não acessar fora da memória double mls0 = masklevelset[vertex0]; double mls1 = masklevelset[vertex1]; @@ -419,8 +444,10 @@ void AdaptiveMeshRefinement::TagElementsNearGroundingLine(TPZGeoMesh *gmesh,std::vector > &GLvec,int &hlevel,std::vector &ElemVec){/*{{{*/ - /* Tag elements near grounding line */ - double MaxRegion = 40000.; //itapopo - double alpha = log(1.5); //itapopo - double MaxDistance = MaxRegion / std::exp(alpha*(hlevel-1)); + /* Tag elements near grounding line */ + double D1 = this->regionlevel1; + double Dhmax = this->regionlevelmax; + int hmax = this->levelmax; + double alpha = (hmax==1) ? 0. : log(D1/Dhmax)/(hmax-1.); + double Di = D1/exp(alpha*(hlevel-1)); for(int i=0;iNElements();i++){ @@ -436,5 +463,5 @@ gmesh->Element(i)->X(qsi, centerPoint); - REAL distance = MaxDistance; + REAL distance = Di; for (int j = 0; j < GLvec.size(); j++) { @@ -442,5 +469,5 @@ REAL value = ( GLvec[j][0] - centerPoint[0] ) * ( GLvec[j][0] - centerPoint[0] ); // (x2-x1)^2 value += ( GLvec[j][1] - centerPoint[1] ) * ( GLvec[j][1] - centerPoint[1] );// (y2-y1)^2 - value = std::sqrt(value); ///Radius + value = std::sqrt(value); //Radius //finding the min distance to the grounding line @@ -449,5 +476,5 @@ } - if(distance < MaxDistance) ElemVec.push_back(i); + if(distance < Di) ElemVec.push_back(i); } @@ -459,10 +486,10 @@ /*NEOPZ works only in C indexing*/ - // itapopo _assert_(nvertices>0); - // itapoo _assert_(nelements>0); + _assert_(nvertices>0); + _assert_(nelements>0); this->SetElementWidth(width); /*Verify and creating initial mesh*/ - //itapopo if(this->fathermesh) _error_("Initial mesh already exists!"); + if(this->fathermesh) _error_("Initial mesh already exists!"); this->fathermesh = new TPZGeoMesh(); @@ -491,5 +518,5 @@ /*reftype = 0: uniform, fast / reftype = 1: uniform and non-uniform (avoid hanging nodes), it is not too fast */ - const int reftype = 1; + const int reftype = 0; switch(this->elementswidth){ case 3: this->fathermesh->CreateGeoElement(ETriangle, elem, mat, index, reftype); break; @@ -521,50 +548,103 @@ /*Create previous mesh as a copy of father mesh*/ this->previousmesh = new TPZGeoMesh(*this->fathermesh); - - /*Set refinement patterns*/ - this->SetRefPatterns(); - -} -/*}}}*/ -void AdaptiveMeshRefinement::SetHMax(int &h){/*{{{*/ - this->hmax = h; + +} +/*}}}*/ +#include "pzgeotriangle.h" //itapopo +#include "pzreftriangle.h" //itapopo +using namespace pzgeom; +TPZGeoMesh* AdaptiveMeshRefinement::CreateRefPatternMesh(TPZGeoMesh* gmesh){/*{{{*/ + + TPZGeoMesh *newgmesh = new TPZGeoMesh(); + newgmesh->CleanUp(); + + int nnodes = gmesh->NNodes(); + int nelem = gmesh->NElements(); + int mat = this->GetElemMaterialID();; + int reftype = 1; + long index; + + //nodes + newgmesh->NodeVec().Resize(nnodes); + for(int i=0;iNodeVec()[i] = gmesh->NodeVec()[i]; + + //elements + for(int i=0;iElement(i); + TPZManVector elem(3,0); + for(int j=0;j<3;j++) elem[j] = geoel->NodeIndex(j); + + newgmesh->CreateGeoElement(ETriangle,elem,mat,index,reftype); + newgmesh->ElementVec()[index]->SetId(geoel->Id()); + + TPZGeoElRefPattern* newgeoel = dynamic_cast*>(newgmesh->ElementVec()[index]); + + //old neighbourhood + const int nsides = TPZGeoTriangle::NSides; + TPZVec< std::vector > neighbourhood(nsides); + TPZVec NodesSequence(0); + for(int s = 0; s < nsides; s++){ + neighbourhood[s].resize(0); + TPZGeoElSide mySide(geoel,s); + TPZGeoElSide neighS = mySide.Neighbour(); + if(mySide.Dimension() == 0){ + long oldSz = NodesSequence.NElements(); + NodesSequence.resize(oldSz+1); + NodesSequence[oldSz] = geoel->NodeIndex(s); + } + while(mySide != neighS){ + neighbourhood[s].push_back(neighS); + neighS = neighS.Neighbour(); + } + } + + //inserting in new element + for(int s = 0; s < nsides; s++){ + TPZGeoEl * tempEl = newgeoel; + TPZGeoElSide tempSide(newgeoel,s); + int byside = s; + for(unsigned long n = 0; n < neighbourhood[s].size(); n++){ + TPZGeoElSide neighS = neighbourhood[s][n]; + tempEl->SetNeighbour(byside, neighS); + tempEl = neighS.Element(); + byside = neighS.Side(); + } + tempEl->SetNeighbour(byside, tempSide); + } + + long fatherindex = geoel->FatherIndex(); + if(fatherindex>-1) newgeoel->SetFather(fatherindex); + + if(!geoel->HasSubElement()) continue; + + int nsons = geoel->NSubElements(); + + TPZAutoPointer ref = gRefDBase.GetUniformRefPattern(ETriangle); + newgeoel->SetRefPattern(ref); + + for(int j=0;jSubElement(j); + if(!son){ + DebugStop(); + } + newgeoel->SetSubElement(j,son); + } + } + newgmesh->BuildConnectivity(); + + return newgmesh; +} +/*}}}*/ +void AdaptiveMeshRefinement::SetLevelMax(int &h){/*{{{*/ + this->levelmax = h; +} +/*}}}*/ +void AdaptiveMeshRefinement::SetRegions(double &D1,double Dhmax){/*{{{*/ + this->regionlevel1 = D1; + this->regionlevelmax = Dhmax; } /*}}}*/ void AdaptiveMeshRefinement::SetElementWidth(int &width){/*{{{*/ this->elementswidth = width; -} -/*}}}*/ -void AdaptiveMeshRefinement::SetRefPatterns(){/*{{{*/ - - /*Initialize the global variable of refinement patterns*/ - gRefDBase.InitializeUniformRefPattern(EOned); - gRefDBase.InitializeUniformRefPattern(ETriangle); - - //gRefDBase.InitializeRefPatterns(); - /*Insert specifics patterns to ISSM core*/ - std::string filepath = REFPATTERNDIR; - std::string filename1 = filepath + "/2D_Triang_Rib_3.rpt"; - std::string filename2 = filepath + "/2D_Triang_Rib_4.rpt"; - std::string filename3 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_4.rpt"; - std::string filename4 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_4_permuted.rpt"; - std::string filename5 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_5.rpt"; - std::string filename6 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_5_permuted.rpt"; - std::string filename7 = filepath + "/2D_Triang_Rib_5.rpt"; - - TPZAutoPointer refpat1 = new TPZRefPattern(filename1); - TPZAutoPointer refpat2 = new TPZRefPattern(filename2); - TPZAutoPointer refpat3 = new TPZRefPattern(filename3); - TPZAutoPointer refpat4 = new TPZRefPattern(filename4); - TPZAutoPointer refpat5 = new TPZRefPattern(filename5); - TPZAutoPointer refpat6 = new TPZRefPattern(filename6); - TPZAutoPointer refpat7 = new TPZRefPattern(filename7); - - if(!gRefDBase.FindRefPattern(refpat1)) gRefDBase.InsertRefPattern(refpat1); - if(!gRefDBase.FindRefPattern(refpat2)) gRefDBase.InsertRefPattern(refpat2); - if(!gRefDBase.FindRefPattern(refpat3)) gRefDBase.InsertRefPattern(refpat3); - if(!gRefDBase.FindRefPattern(refpat4)) gRefDBase.InsertRefPattern(refpat4); - if(!gRefDBase.FindRefPattern(refpat5)) gRefDBase.InsertRefPattern(refpat5); - if(!gRefDBase.FindRefPattern(refpat6)) gRefDBase.InsertRefPattern(refpat6); - if(!gRefDBase.FindRefPattern(refpat7)) gRefDBase.InsertRefPattern(refpat7); } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h =================================================================== --- /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h (revision 21671) +++ /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h (revision 21672) @@ -58,9 +58,11 @@ void CleanUp(); // Clean all attributes void Initialize(); // Initialize the attributes with NULL and values out of usually range - void SetHMax(int &h); // Define the max level of refinement - void SetElementWidth(int &width); // Define elements width + void SetLevelMax(int &h); // Define the max level of refinement + void SetRegions(double &D1,double Dhmax); // Define the regions which will be refined + void SetElementWidth(int &width); // Define elements width 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 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 - 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 + TPZGeoMesh* CreateRefPatternMesh(TPZGeoMesh* gmesh); + 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 private: @@ -68,5 +70,7 @@ /*Private attributes*/ int elementswidth; // Geometric nodes for element: 3 == Tria, 4 == Tetra, 6 == Penta - int hmax; // Max level of refinement + int levelmax; // Max level of refinement + double regionlevel1; // Region which will be refined with level 1 + double regionlevelmax; // Region which will be refined with level max TPZGeoMesh *fathermesh; // Father Mesh is the entire mesh without refinement TPZGeoMesh *previousmesh; // Previous mesh is a refined mesh of last step @@ -83,5 +87,4 @@ inline int GetElemMaterialID(){return 1;} // Return element material ID inline int GetBoundaryMaterialID(){return 2;} // Return segment (2D boundary) material ID - void SetRefPatterns(); // Initialize and insert the refinement patterns }; Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 21671) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 21672) @@ -2921,14 +2921,17 @@ /*Define variables*/ - int my_rank = IssmComm::GetRank(); - this->amr = NULL;//initialize amr as NULL - int numberofvertices = this->vertices->NumberOfVertices(); - int numberofelements = this->elements->NumberOfElements(); - int numberofsegments = 0; //used on matlab - IssmDouble* x = NULL; - IssmDouble* y = NULL; - IssmDouble* z = NULL; - int* elements = NULL; - int elementswidth = this->GetElementsWidth(); //just tria elements in this version. Itapopo: + int my_rank = IssmComm::GetRank(); + this->amr = NULL;//initialize amr as NULL + int numberofvertices = this->vertices->NumberOfVertices(); + int numberofelements = this->elements->NumberOfElements(); + int numberofsegments = 0; //used on matlab + IssmDouble* x = NULL; + IssmDouble* y = NULL; + IssmDouble* z = NULL; + int* elements = NULL; + int elementswidth = this->GetElementsWidth(); //just tria elements in this version. Itapopo: + int levelmax = 0; + IssmDouble regionlevel1 = 0.; + IssmDouble regionlevelmax = 0.; /*Get vertices coordinates of the coarse mesh (father mesh)*/ @@ -2936,11 +2939,33 @@ this->GetMesh(this->vertices,this->elements,&x,&y,&z,&elements); + /*Get amr parameters*/ + this->parameters->FindParam(&levelmax,AmrLevelMaxEnum); + this->parameters->FindParam(®ionlevel1,AmrRegionLevel1Enum); + this->parameters->FindParam(®ionlevelmax,AmrRegionLevelMaxEnum); + /*Create initial mesh (coarse mesh) in neopz data structure*/ /*Just CPU #0 should keep AMR object*/ - if(my_rank==0){ - this->amr = new AdaptiveMeshRefinement(); - int hmax = 0; //itapopo: it must be defined by interface. Using 2 just to test - this->amr->SetHMax(hmax); //Set max level of refinement - this->amr->CreateInitialMesh(numberofvertices,numberofelements,numberofsegments,elementswidth,x,y,z,elements,NULL); + this->SetRefPatterns(); + if(my_rank==0){ + bool ismisomip = true; + if(ismisomip){//itapopo + TPZFileStream fstr; + std::stringstream ss; + + ss << levelmax; + std::string AMRfile = "/home/santos/Misomip2/L" + ss.str() + "_tsai/amr.txt"; + fstr.OpenRead(AMRfile.c_str()); + + TPZSaveable *sv = TPZSaveable::Restore(fstr,0); + this->amr = dynamic_cast(sv); + } + else{ + this->amr = new AdaptiveMeshRefinement(); + //this->amr->SetLevelMax(levelmax); //Set max level of refinement + //this->amr->SetRegions(regionlevel1,regionlevelmax); + this->amr->CreateInitialMesh(numberofvertices,numberofelements,numberofsegments,elementswidth,x,y,z,elements,NULL); + } + this->amr->SetLevelMax(levelmax); //Set max level of refinement + this->amr->SetRegions(regionlevel1,regionlevelmax); } @@ -2951,4 +2976,38 @@ xDelete(elements); +} +/*}}}*/ +void FemModel::SetRefPatterns(){/*{{{*/ + + /*Initialize the global variable of refinement patterns*/ + gRefDBase.InitializeUniformRefPattern(EOned); + gRefDBase.InitializeUniformRefPattern(ETriangle); + + //gRefDBase.InitializeRefPatterns(); + /*Insert specifics patterns to ISSM core*/ + std::string filepath = REFPATTERNDIR; + std::string filename1 = filepath + "/2D_Triang_Rib_3.rpt"; + std::string filename2 = filepath + "/2D_Triang_Rib_4.rpt"; + std::string filename3 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_4.rpt"; + std::string filename4 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_4_permuted.rpt"; + std::string filename5 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_5.rpt"; + std::string filename6 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_5_permuted.rpt"; + std::string filename7 = filepath + "/2D_Triang_Rib_5.rpt"; + + TPZAutoPointer refpat1 = new TPZRefPattern(filename1); + TPZAutoPointer refpat2 = new TPZRefPattern(filename2); + TPZAutoPointer refpat3 = new TPZRefPattern(filename3); + TPZAutoPointer refpat4 = new TPZRefPattern(filename4); + TPZAutoPointer refpat5 = new TPZRefPattern(filename5); + TPZAutoPointer refpat6 = new TPZRefPattern(filename6); + TPZAutoPointer refpat7 = new TPZRefPattern(filename7); + + if(!gRefDBase.FindRefPattern(refpat1)) gRefDBase.InsertRefPattern(refpat1); + if(!gRefDBase.FindRefPattern(refpat2)) gRefDBase.InsertRefPattern(refpat2); + if(!gRefDBase.FindRefPattern(refpat3)) gRefDBase.InsertRefPattern(refpat3); + if(!gRefDBase.FindRefPattern(refpat4)) gRefDBase.InsertRefPattern(refpat4); + if(!gRefDBase.FindRefPattern(refpat5)) gRefDBase.InsertRefPattern(refpat5); + if(!gRefDBase.FindRefPattern(refpat6)) gRefDBase.InsertRefPattern(refpat6); + if(!gRefDBase.FindRefPattern(refpat7)) gRefDBase.InsertRefPattern(refpat7); } /*}}}*/ @@ -3067,4 +3126,10 @@ GetMaskOfIceVerticesLSMx(this); + /*Insert MISMIP+ bed topography*/ + if(true) this->BedrockFromMismipPlus(); + + /*Adjust base, thickness and mask grounded ice leve set*/ + if(true) this->AdjustBaseThicknessAndMask(); + /*Reset current configuration: */ analysis_type=this->analysis_type_list[this->analysis_counter]; @@ -3081,4 +3146,136 @@ return; +} +/*}}}*/ +void FemModel::BedrockFromMismipPlus(void){/*{{{*/ + + /*Insert bedrock from mismip+ setup*/ + /*This was used to Misomip project/simulations*/ + + IssmDouble x,y,bx,by; + int numvertices = this->GetElementsWidth(); + IssmDouble* xyz_list = NULL; + IssmDouble* r = xNew(numvertices); + + for(int el=0;elelements->Size();el++){ + Element* element=xDynamicCast(this->elements->GetObjectByOffset(el)); + + element->GetVerticesCoordinates(&xyz_list); + for(int i=0;iAddInput(BedEnum,&r[0],P1Enum); + + /*Cleanup*/ + xDelete(xyz_list); + } + + /*Delete*/ + xDelete(r); + + return; +} +/*}}}*/ +void FemModel::AdjustBaseThicknessAndMask(void){/*{{{*/ + + int numvertices = this->GetElementsWidth(); + IssmDouble rho_water,rho_ice,density,base_float; + IssmDouble* phi = xNew(numvertices); + IssmDouble* h = xNew(numvertices); + IssmDouble* s = xNew(numvertices); + IssmDouble* b = xNew(numvertices); + IssmDouble* r = xNew(numvertices); + IssmDouble* sl = xNew(numvertices); + + for(int el=0;elelements->Size();el++){ + Element* element=xDynamicCast(this->elements->GetObjectByOffset(el)); + + element->GetInputListOnVertices(&s[0],SurfaceEnum); + element->GetInputListOnVertices(&r[0],BedEnum); + element->GetInputListOnVertices(&sl[0],SealevelEnum); + rho_water = element->matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum); + rho_ice = element->matpar->GetMaterialParameter(MaterialsRhoIceEnum); + density = rho_ice/rho_water; + + for(int i=0;ibase_float){ + b[i] = r[i]; + } + else { + b[i] = base_float; + } + + if(abs(sl[i])>0) _error_("Sea level value not supported!"); + /*update thickness and mask grounded ice level set*/ + h[i] = s[i]-b[i]; + phi[i] = h[i]+r[i]/density; + } + + /*Update inputs*/ + element->AddInput(MaskGroundediceLevelsetEnum,&phi[0],P1Enum); + element->AddInput(ThicknessEnum,&h[0],P1Enum); + element->AddInput(BaseEnum,&b[0],P1Enum); + + } + + /*Delete*/ + xDelete(phi); + xDelete(h); + xDelete(s); + xDelete(b); + xDelete(r); + xDelete(sl); + + return; +} +/*}}}*/ +void FemModel::WriteMeshInResults(void){/*{{{*/ + + int step = -1; + int numberofelements = -1; + int numberofvertices = -1; + IssmDouble time = -1; + IssmDouble* x = NULL; + IssmDouble* y = NULL; + IssmDouble* z = NULL; + int* elementslist = NULL; + + if(!this->elements || !this->vertices || !this->results || !this->parameters) return; + + parameters->FindParam(&step,StepEnum); + parameters->FindParam(&time,TimeEnum); + numberofelements=this->elements->NumberOfElements(); + numberofvertices=this->vertices->NumberOfVertices(); + + /*Get mesh. Elementslist comes in Matlab indexing*/ + this->GetMesh(this->vertices,this->elements,&x,&y,&z,&elementslist); + + /*Write mesh in Results*/ + this->results->AddResult(new GenericExternalResult(this->results->Size()+1,MeshElementsEnum, + elementslist,numberofelements,this->GetElementsWidth(),step,time)); + + this->results->AddResult(new GenericExternalResult(this->results->Size()+1,MeshXEnum, + x,numberofvertices,1,step,time)); + + this->results->AddResult(new GenericExternalResult(this->results->Size()+1,MeshYEnum, + y,numberofvertices,1,step,time)); + + /*Cleanup*/ + xDelete(x); + xDelete(y); + xDelete(z); + xDelete(elementslist); + + return; } /*}}}*/ @@ -3216,5 +3413,5 @@ P1inputsold,numverticesold,numP1inputs, Xnew,Ynew,numverticesnew,NULL); - + /*Insert P0 inputs into the new elements.*/ vector=NULL; @@ -3251,5 +3448,5 @@ for(int i=0;i(numverticesnew); for(int j=0;j