Index: ../trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp =================================================================== --- ../trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp (revision 21671) +++ ../trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp (revision 21672) @@ -26,15 +26,30 @@ this->CleanUp(); /*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; } /*}}}*/ 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(); } @@ -44,17 +59,21 @@ /*Verify and delete all data*/ 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; } /*}}}*/ void AdaptiveMeshRefinement::Initialize(){/*{{{*/ /*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; } /*}}}*/ int AdaptiveMeshRefinement::ClassId() const{/*{{{*/ @@ -81,10 +100,12 @@ } /* Read simple attributes */ - buf.Read(&this->hmax,1); + buf.Read(&this->levelmax,1); buf.Read(&this->elementswidth,1); + buf.Read(&this->regionlevel1,1); + buf.Read(&this->regionlevelmax,1); - /* Read geometric mesh*/ + /* Read geometric mesh*/ TPZSaveable *sv1 = TPZSaveable::Restore(buf,0); this->fathermesh = dynamic_cast(sv1); @@ -113,9 +134,11 @@ buf.Write(&classid,1); /* 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()); this->previousmesh->Write(buf, this->ClassId()); @@ -137,15 +160,15 @@ /*IMPORTANT! pelements (and psegments) are in Matlab indexing*/ /*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*/ std::vector > GLvec; 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*/ TPZGeoMesh *newmesh; @@ -157,8 +180,8 @@ 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*/ if(type_process==1){ @@ -167,15 +190,16 @@ } /*Refine elements to avoid hanging nodes*/ - TPZGeoMesh *nohangingnodesmesh = new TPZGeoMesh(*newmesh); + //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("/ronne_1/home/santos/mesh2.vtk"); - TPZVTKGeoMesh::PrintGMeshVTK(newmesh,file3); + //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*/ this->GetMesh(nohangingnodesmesh,nvertices,nelements,nsegments,px,py,pz,pelements,psegments); @@ -190,10 +214,10 @@ /*}}}*/ 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*/ std::vector ElemVec; //elements without children @@ -364,6 +388,7 @@ int vertex1 = this->previousmesh->Element(i)->NodeIndex(1); 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]; double mls2 = masklevelset[vertex2]; @@ -418,10 +443,12 @@ /*}}}*/ 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++){ @@ -435,20 +462,20 @@ gmesh->Element(i)->CenterPoint(side2D, qsi); gmesh->Element(i)->X(qsi, centerPoint); - REAL distance = MaxDistance; + REAL distance = Di; for (int j = 0; j < GLvec.size(); j++) { 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 if(value < distance) distance = value; } - if(distance < MaxDistance) ElemVec.push_back(i); + if(distance < Di) ElemVec.push_back(i); } } @@ -458,12 +485,12 @@ /*IMPORTANT! elements come in Matlab indexing*/ /*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(); this->fathermesh->NodeVec().Resize( nvertices ); @@ -490,7 +517,7 @@ for(int jel=0;jelelementswidth;jel++) elem[jel]=elements[iel*this->elementswidth+jel]-1;//Convert Matlab to C indexing /*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; case 4: this->fathermesh->CreateGeoElement(ETetraedro, elem, mat, index, reftype); DebugStop(); break; @@ -520,54 +547,107 @@ /*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); -} -/*}}}*/ void AdaptiveMeshRefinement::CheckMesh(int &nvertices,int &nelements,int &nsegments,int &width,double** px,double** py,double** pz,int** pelements, int** psegments){/*{{{*/ /*Basic verification*/ Index: ../trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- ../trunk-jpl/src/c/classes/FemModel.cpp (revision 21671) +++ ../trunk-jpl/src/c/classes/FemModel.cpp (revision 21672) @@ -2920,28 +2920,53 @@ void FemModel::InitializeAdaptiveRefinement(void){/*{{{*/ /*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)*/ /*elements comes in Matlab indexing*/ 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); } /*Free the vectors*/ @@ -2952,6 +2977,40 @@ } /*}}}*/ +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); +} +/*}}}*/ void FemModel::ReMesh(void){/*{{{*/ /*Variables*/ @@ -3066,6 +3125,12 @@ 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]; SetCurrentConfiguration(analysis_type); @@ -3082,6 +3147,138 @@ } /*}}}*/ +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; +} +/*}}}*/ void FemModel::InterpolateInputs(Vertices* newfemmodel_vertices,Elements* newfemmodel_elements){/*{{{*/ int maxinputs = MaximumNumberOfDefinitionsEnum; @@ -3215,7 +3412,7 @@ InterpFromMeshToMesh2dx(&P1inputsnew,Indexold,Xold,Yold,numverticesold,numelementsold, P1inputsold,numverticesold,numP1inputs, Xnew,Ynew,numverticesnew,NULL); - + /*Insert P0 inputs into the new elements.*/ vector=NULL; for(int i=0;i(numverticesnew); for(int j=0;j