Index: /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp (revision 23500) +++ /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp (revision 23501) @@ -15,5 +15,26 @@ /*Constructor, copy, clean up and destructor*/ AdaptiveMeshRefinement::AdaptiveMeshRefinement(){/*{{{*/ - this->Initialize(); + + /*Set pointers to NULL*/ + this->fathermesh = NULL; + this->previousmesh = NULL; + this->refinement_type = -1; + this->level_max = -1; + this->gradation = -1; + this->lag = -1; + this->groundingline_distance = -1; + this->icefront_distance = -1; + this->thicknesserror_threshold = -1; + this->deviatoricerror_threshold = -1; + this->deviatoricerror_maximum = -1; + this->thicknesserror_maximum = -1; + this->sid2index.clear(); + this->index2sid.clear(); + this->specialelementsindex.clear(); + this->x = NULL; + this->y = NULL; + this->elementslist = NULL; + this->numberofvertices = -1; + this->numberofelements = -1; } /*}}}*/ @@ -54,5 +75,5 @@ /*}}}*/ AdaptiveMeshRefinement::~AdaptiveMeshRefinement(){/*{{{*/ - int writemesh = 1;//itapopo keep 0 + int writemesh = 0;//only to restart if(writemesh) this->WriteMesh(); this->CleanUp(); @@ -65,4 +86,7 @@ if(this->fathermesh) delete this->fathermesh; if(this->previousmesh) delete this->previousmesh; + if(this->x) delete this->x; + if(this->y) delete this->y; + if(this->elementslist) delete this->elementslist; this->refinement_type = -1; this->level_max = -1; @@ -75,4 +99,6 @@ this->deviatoricerror_maximum = -1; this->thicknesserror_maximum = -1; + this->numberofvertices = -1; + this->numberofelements = -1; this->sid2index.clear(); this->index2sid.clear(); @@ -80,26 +106,28 @@ } /*}}}*/ -void AdaptiveMeshRefinement::Initialize(){/*{{{*/ - - /*Set pointers to NULL*/ - this->fathermesh = NULL; - this->previousmesh = NULL; - this->refinement_type = -1; - this->level_max = -1; - this->gradation = -1; - this->lag = -1; - this->groundingline_distance = -1; - this->icefront_distance = -1; - this->thicknesserror_threshold = -1; - this->deviatoricerror_threshold = -1; - this->deviatoricerror_maximum = -1; - this->thicknesserror_maximum = -1; - this->sid2index.clear(); - this->index2sid.clear(); - this->specialelementsindex.clear(); -} -/*}}}*/ /*Mesh refinement methods*/ +void AdaptiveMeshRefinement::SetMesh(int** elementslist_in,IssmDouble** x_in,IssmDouble** y_in,int* numberofvertices_in,int* numberofelements_in){/*{{{*/ + + /*Delete previous mesh and keep the entire mesh*/ + if(this->elementslist) xDelete(this->elementslist); + if(this->x) xDelete(this->x); + if(this->y) xDelete(this->y); + + this->elementslist = *elementslist_in; + this->x = *x_in; + this->y = *y_in; + this->numberofvertices = *numberofvertices_in; + this->numberofelements = *numberofelements_in; +}/*}}}*/ +void AdaptiveMeshRefinement::GetMesh(int** elementslist_out,IssmDouble** x_out,IssmDouble** y_out,int* numberofvertices_out,int* numberofelements_out){/*{{{*/ + + /*Get the entire mesh*/ + *elementslist_out = this->elementslist; + *x_out = this->x; + *y_out = this->y; + *numberofvertices_out= this->numberofvertices; + *numberofelements_out= this->numberofelements; +}/*}}}*/ void AdaptiveMeshRefinement::ExecuteRefinement(double* gl_distance,double* if_distance,double* deviatoricerror,double* thicknesserror,int** pdatalist,double** pxylist,int** pelementslist){/*{{{*/ @@ -555,25 +583,26 @@ } /*}}}*/ -void AdaptiveMeshRefinement::CreateInitialMesh(int &nvertices,int &nelements,double* x,double* y,int* elements){/*{{{*/ +void AdaptiveMeshRefinement::Initialize(){/*{{{*/ /* IMPORTANT! elements come in Matlab indexing NEOPZ works only in C indexing*/ - if(nvertices<=0) _error_("Impossible to create initial mesh: nvertices is <= 0!\n"); - if(nelements<=0) _error_("Impossible to create initial mesh: nelements is <= 0!\n"); + if(this->numberofvertices<=0) _error_("Impossible to create initial mesh: nvertices is <= 0!\n"); + if(this->numberofelements<=0) _error_("Impossible to create initial mesh: nelements is <= 0!\n"); if(this->refinement_type!=0 && this->refinement_type!=1) _error_("Impossible to create initial mesh: refinement type is not defined!\n"); /*Verify and creating initial mesh*/ - if(this->fathermesh || this->previousmesh) _error_("Initial mesh already exists!"); + if(!this->x || !this->y || !this->elementslist) _error_("Mesh data structure is NULL!\n"); + if(this->fathermesh || this->previousmesh) _error_("Initial mesh already exists!\n"); this->fathermesh = new TPZGeoMesh(); - this->fathermesh->NodeVec().Resize(nvertices); + this->fathermesh->NodeVec().Resize(this->numberofvertices); /*Set the vertices (geometric nodes in NeoPZ context)*/ - for(int i=0;inumberofvertices;i++){ /*x,y,z coords*/ TPZManVector coord(3,0.); - coord[0]= x[i]; - coord[1]= y[i]; + coord[0]= this->x[i]; + coord[1]= this->y[i]; coord[2]= 0.; /*Insert in the mesh*/ @@ -586,9 +615,9 @@ const int mat = this->GetElemMaterialID(); TPZManVector elem(this->GetNumberOfNodes(),0); - this->index2sid.clear(); this->index2sid.resize(nelements); + this->index2sid.clear(); this->index2sid.resize(this->numberofelements); this->sid2index.clear(); - for(int i=0;iGetNumberOfNodes();j++) elem[j]=elements[i*this->GetNumberOfNodes()+j]-1;//Convert Matlab to C indexing + for(int i=0;inumberofelements;i++){ + for(int j=0;jGetNumberOfNodes();j++) elem[j]=this->elementslist[i*this->GetNumberOfNodes()+j]-1;//Convert Matlab to C indexing switch(this->GetNumberOfNodes()){ case 3: this->fathermesh->CreateGeoElement(ETriangle,elem,mat,index,this->refinement_type); break; @@ -865,5 +894,5 @@ std::ifstream amrifstream(amrfile.c_str()); int size,value; - + TPZPersistenceManager::OpenRead(fathermeshfile); this->fathermesh = dynamic_cast(TPZPersistenceManager::ReadFromFile()); Index: /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h =================================================================== --- /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h (revision 23500) +++ /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h (revision 23501) @@ -68,5 +68,6 @@ void Initialize(); void ExecuteRefinement(double* gl_distance,double* if_distance,double* deviatoricerror,double* thicknesserror,int** pdatalist,double** pxy,int** pelementslist); - void CreateInitialMesh(int &nvertices,int &nelements,double* x,double* y,int* elements); + void SetMesh(int** elementslist_in,IssmDouble** x_in,IssmDouble** y_in,int* numberofvertices,int* numberofelements); + void GetMesh(int** elementslist_out,IssmDouble** x_out,IssmDouble** y_out,int* numberofvertices,int* numberofelements); void CheckMesh(int** pdata,double** pxy,int** pelements); void ReadMesh(); @@ -80,4 +81,9 @@ TPZGeoMesh *fathermesh; // Entire mesh without refinement if refinement_type==1; refined with hanging nodes if efinement_type==0 TPZGeoMesh *previousmesh; // Refined mesh without hanging nodes (it is always refpattern type), used to generate ISSM mesh + IssmDouble* x; // Entire mesh + IssmDouble* y; + int* elementslist; + int numberofvertices; + int numberofelements; /*}}}*/ /*Private methods{{{*/ Index: /issm/trunk-jpl/src/c/classes/AmrBamg.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/AmrBamg.cpp (revision 23500) +++ /issm/trunk-jpl/src/c/classes/AmrBamg.cpp (revision 23501) @@ -39,4 +39,9 @@ this->fathermesh = NULL; this->previousmesh = NULL; + this->elementslist = NULL; + this->x = NULL; + this->y = NULL; + this->numberofvertices = -1; + this->numberofelements = -1; /*Only initialize options for now (same as bamg.m)*/ @@ -73,10 +78,34 @@ if(this->previousmesh) delete this->previousmesh; if(this->options) delete this->options; - + if(this->x) xDelete(this->x); + if(this->y) xDelete(this->y); + if(this->elementslist) xDelete(this->elementslist); } /*}}}*/ /*Methods*/ -void AmrBamg::Initialize(int* elements,IssmDouble* x,IssmDouble* y,int numberofvertices,int numberofelements){/*{{{*/ +void AmrBamg::SetMesh(int** elementslist_in,IssmDouble** x_in,IssmDouble** y_in,int* numberofvertices_in,int* numberofelements_in){/*{{{*/ + + /*Delete previous mesh and keep the entire mesh*/ + if(this->elementslist) xDelete(this->elementslist); + if(this->x) xDelete(this->x); + if(this->y) xDelete(this->y); + + this->elementslist = *elementslist_in; + this->x = *x_in; + this->y = *y_in; + this->numberofvertices = *numberofvertices_in; + this->numberofelements = *numberofelements_in; +}/*}}}*/ +void AmrBamg::GetMesh(int** elementslist_out,IssmDouble** x_out,IssmDouble** y_out,int* numberofvertices_out,int* numberofelements_out){/*{{{*/ + + /*Get the entire mesh*/ + *elementslist_out = this->elementslist; + *x_out = this->x; + *y_out = this->y; + *numberofvertices_out= this->numberofvertices; + *numberofelements_out= this->numberofelements; +}/*}}}*/ +void AmrBamg::Initialize(){/*{{{*/ /*Check options*/ @@ -85,5 +114,5 @@ /*Read father mesh and create geometry*/ - Mesh* Th=new Mesh(elements,x,y,numberofvertices,numberofelements,this->options); + Mesh* Th=new Mesh(this->elementslist,this->x,this->y,this->numberofvertices,this->numberofelements,this->options); /*Write geometry*/ @@ -99,5 +128,5 @@ }/*}}}*/ void AmrBamg::ExecuteRefinementBamg(IssmDouble* field,IssmDouble* hmaxVertices,int** pdatalist,IssmDouble** pxylist,int** pelementslist){/*{{{*/ - + /*Intermediaries*/ BamgGeom* geomout=new BamgGeom(); @@ -113,5 +142,5 @@ this->options->field = field; this->options->hmaxVertices = hmaxVertices; - + /*remesh*/ if(this->previousmesh){ Index: /issm/trunk-jpl/src/c/classes/AmrBamg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/AmrBamg.h (revision 23500) +++ /issm/trunk-jpl/src/c/classes/AmrBamg.h (revision 23501) @@ -28,4 +28,5 @@ IssmDouble deviatoricerror_maximum; + /* Constructor, destructor etc*/ AmrBamg(); @@ -34,5 +35,7 @@ /*General methods*/ - void Initialize(int* elements,IssmDouble* x,IssmDouble* y,int numberofvertices,int numberofelements); + void Initialize(); + void SetMesh(int** elementslist_in,IssmDouble** x_in,IssmDouble** y_in,int* numberofvertices,int* numberofelements); + void GetMesh(int** elementslist_out,IssmDouble** x_out,IssmDouble** y_out,int* numberofvertices,int* numberofelements); void ExecuteRefinementBamg(IssmDouble* field,IssmDouble* hmaxVertices,int** pdatalist,IssmDouble** pxylist,int** pelementslist); void SetBamgOpts(IssmDouble hmin_in,IssmDouble hmax_in,IssmDouble err_in,IssmDouble gradation_in); @@ -46,4 +49,10 @@ BamgMesh* previousmesh; BamgOpts* options; + /*entire mesh*/ + IssmDouble* x; + IssmDouble* y; + int* elementslist; + int numberofvertices; + int numberofelements; }; Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 23500) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 23501) @@ -2719,5 +2719,5 @@ /*Finally: interpolate all inputs and insert them into the new elements.*/ this->InterpolateInputs(new_vertices,new_elements); - + /*Delete old structure and set new pointers*/ delete this->vertices; this->vertices = new_vertices; @@ -2749,9 +2749,9 @@ SetCurrentConfiguration(analysis_type); + /*Set the new mesh*/ + this->SetMesh(&newelementslist,&newx,&newy,&newnumberofvertices,&newnumberofelements); + /*Cleanup*/ - xDelete(newx); - xDelete(newy); xDelete(newz); - xDelete(newelementslist); xDelete(my_vertices); xDelete(my_elements); @@ -2967,7 +2967,7 @@ void FemModel::InterpolateInputs(Vertices* newfemmodel_vertices,Elements* newfemmodel_elements){/*{{{*/ - int numberofelements = this->elements->NumberOfElements(); //global, entire old mesh + int numberofelements = -1; //global, entire old mesh int newnumberofelements = newfemmodel_elements->Size(); //just on the new partition - int numberofvertices = this->vertices->NumberOfVertices(); //global, entire old mesh + int numberofvertices = -1; //global, entire old mesh int newnumberofvertices = newfemmodel_vertices->Size(); //just on the new partition int elementswidth = this->GetElementsWidth(); //just tria in this version @@ -2986,5 +2986,4 @@ IssmDouble* x = NULL;//global, entire old mesh IssmDouble* y = NULL;//global, entire old mesh - IssmDouble* z = NULL;//global, entire old mesh int* elementslist = NULL;//global, entire old mesh IssmDouble* newx = NULL;//just on the new partition @@ -3000,5 +2999,5 @@ /*Get the old mesh (global, entire mesh)*/ - this->GetMesh(this->vertices,this->elements,&x,&y,&z,&elementslist); + this->GetMesh(&elementslist,&x,&y,&numberofvertices,&numberofelements); /*Get the new mesh (just on the new partition)*/ @@ -3065,8 +3064,4 @@ xDelete(P1input_enums); xDelete(P1input_interp); - xDelete(x); - xDelete(y); - xDelete(z); - xDelete(elementslist); xDelete(newx); xDelete(newy); @@ -3090,5 +3085,4 @@ IssmDouble* x = NULL; IssmDouble* y = NULL; - IssmDouble* z = NULL; int* elementslist = NULL; @@ -3097,9 +3091,7 @@ 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); + this->GetMesh(&elementslist,&x,&y,&numberofvertices,&numberofelements); /*Write mesh in Results*/ @@ -3112,10 +3104,4 @@ this->results->AddResult(new GenericExternalResult(this->results->Size()+1,MeshYEnum, y,numberofvertices,1,step,time)); - - /*Cleanup*/ - xDelete(x); - xDelete(y); - xDelete(z); - xDelete(elementslist); } /*}}}*/ @@ -3326,5 +3312,5 @@ } /*}}}*/ -void FemModel::GetMesh(Vertices* femmodel_vertices, Elements* femmodel_elements,IssmDouble** px, IssmDouble** py, IssmDouble** pz, int** pelementslist){/*{{{*/ +void FemModel::GetMesh(Vertices* femmodel_vertices, Elements* femmodel_elements,IssmDouble** px, IssmDouble** py, int** pelementslist){/*{{{*/ if(!femmodel_vertices) _error_("GetMesh: vertices are NULL."); @@ -3383,5 +3369,4 @@ *px = x; *py = y; - *pz = z; *pelementslist = elementslist; //Matlab indexing. InterMesh uses this type. @@ -3391,4 +3376,5 @@ xDelete(id2); xDelete(id3); + xDelete(z); delete vid1; delete vid2; @@ -3396,4 +3382,40 @@ } /*}}}*/ +void FemModel::GetMesh(int** elementslist, IssmDouble** x, IssmDouble** y, int* numberofvertices, int* numberofelements){/*{{{*/ + + int amrtype; + this->parameters->FindParam(&amrtype,AmrTypeEnum); + + switch(amrtype){ + + #if defined(_HAVE_NEOPZ_) + case AmrNeopzEnum: this->amr->GetMesh(elementslist,x,y,numberofvertices,numberofelements); break; + #endif + + #if defined(_HAVE_BAMG_) + case AmrBamgEnum: this->amrbamg->GetMesh(elementslist,x,y,numberofvertices,numberofelements); break; + #endif + + default: _error_("not implemented yet"); + } +}/*}}}*/ +void FemModel::SetMesh(int** elementslist, IssmDouble** x, IssmDouble** y, int* numberofvertices, int* numberofelements){/*{{{*/ + + int amrtype; + this->parameters->FindParam(&amrtype,AmrTypeEnum); + + switch(amrtype){ + + #if defined(_HAVE_NEOPZ_) + case AmrNeopzEnum: this->amr->SetMesh(elementslist,x,y,numberofvertices,numberofelements); break; + #endif + + #if defined(_HAVE_BAMG_) + case AmrBamgEnum: this->amrbamg->SetMesh(elementslist,x,y,numberofvertices,numberofelements); break; + #endif + + default: _error_("not implemented yet"); + } +}/*}}}*/ void FemModel::GetMeshOnPartition(Vertices* femmodel_vertices,Elements* femmodel_elements,IssmDouble** px,IssmDouble** py,IssmDouble** pz,int** pelementslist,int** psidtoindex){/*{{{*/ @@ -3460,14 +3482,13 @@ int analysis_index = AnalysisIndex(analysis_enum); - int numberofnodes_analysistype= this->nodes_list[analysis_index]->NumberOfNodes(); + int numberofnodes_analysistype= this->nodes_list[analysis_index]->NumberOfNodes(); int dofpernode = 2; //vx and vy int numberofcols = dofpernode*2; //to keep dofs and flags in the vspc vector - int numberofvertices = this->vertices->NumberOfVertices(); //global, entire old mesh - int numberofelements = this->elements->NumberOfElements(); //global, entire old mesh + int numberofvertices = -1; //global, entire old mesh + int numberofelements = -1; //global, entire old mesh int newnumberofvertices = newfemmodel_vertices->Size(); //local, just the new partition int count = 0; IssmDouble* x = NULL; //global, entire old mesh IssmDouble* y = NULL; //global, entire old mesh - IssmDouble* z = NULL; //global, entire old mesh int* elementslist = NULL; //global, entire old mesh IssmDouble* spc = NULL; //global, entire old mesh @@ -3479,5 +3500,5 @@ /*Get old mesh (global, entire mesh). Elementslist comes in Matlab indexing*/ - this->GetMesh(this->vertices,this->elements,&x,&y,&z,&elementslist); + this->GetMesh(&elementslist,&x,&y,&numberofvertices,&numberofelements); /*Get vertices coordinates of the new partition*/ @@ -3545,10 +3566,6 @@ } } - + /*Cleanup*/ - xDelete(x); - xDelete(y); - xDelete(z); - xDelete(elementslist); xDelete(spc); xDelete(newspc); @@ -5077,5 +5094,4 @@ IssmDouble* x = NULL; IssmDouble* y = NULL; - IssmDouble* z = NULL; int* elements = NULL; IssmDouble hmin,hmax,err,gradation; @@ -5088,5 +5104,5 @@ /*Get vertices coordinates of the coarse mesh (father mesh)*/ - this->GetMesh(this->vertices,this->elements,&x,&y,&z,&elements); + this->GetMesh(this->vertices,this->elements,&x,&y,&elements); /*Create bamg data structures for bamg*/ @@ -5116,13 +5132,8 @@ /*Re-create original mesh and put it in bamg structure (only cpu 0)*/ + this->amrbamg->SetMesh(&elements,&x,&y,&numberofvertices,&numberofelements); if(my_rank==0){ - this->amrbamg->Initialize(elements,x,y,numberofvertices,numberofelements); - } - - /*Free the vectors*/ - xDelete(x); - xDelete(y); - xDelete(z); - xDelete(elements); + this->amrbamg->Initialize(); + } } /*}}}*/ @@ -5132,7 +5143,12 @@ /*Intermediaries*/ - int numberofvertices = this->vertices->NumberOfVertices(); - IssmDouble* verticedistance = NULL; + int numberofvertices = this->vertices->NumberOfVertices(); + Vector* vminvertexdistance = new Vector(numberofvertices); + IssmDouble* pminvertexdistance = NULL; + IssmDouble* levelset_points = NULL; + IssmDouble x,y; IssmDouble threshold,resolution; + IssmDouble minvertexdistance,distance; + int sid,numberofpoints; switch(levelset_type){ @@ -5148,11 +5164,32 @@ } - /*Get vertice distance to zero level set points*/ - this->GetVerticeDistanceToZeroLevelSet(&verticedistance,levelset_type); - if(!verticedistance) _error_("verticedistance is NULL!\n"); + /*Get points which level set is zero (center of elements with zero level set)*/ + this->GetZeroLevelSetPoints(&levelset_points,numberofpoints,levelset_type);//levelset_points is serial (global) + + for(int i=0;ivertices->Size();i++){//only on this partition + Vertex* vertex=(Vertex*)this->vertices->GetObjectByOffset(i); + /*Attention: no spherical coordinates*/ + x = vertex->GetX(); + y = vertex->GetY(); + sid= vertex->Sid(); + minvertexdistance=INFINITY; + + /*Find the minimum vertex distance*/ + for(int j=0;jSetValue(sid,minvertexdistance,INS_VAL); + } + /*Assemble*/ + vminvertexdistance->Assemble(); + + /*Assign the pointer*/ + pminvertexdistance=vminvertexdistance->ToMPISerial(); /*Fill hmaxVertices*/ for(int i=0;i(hmaxvertices[i])) hmaxvertices[i]=resolution; else hmaxvertices[i]=min(resolution,hmaxvertices[i]); @@ -5161,5 +5198,7 @@ /*Cleanup*/ - xDelete(verticedistance); + xDelete(pminvertexdistance); + xDelete(levelset_points); + delete vminvertexdistance; } /*}}}*/ @@ -5170,13 +5209,12 @@ /*Intermediaries*/ int elementswidth = this->GetElementsWidth(); - int numberofelements = this->elements->NumberOfElements(); - int numberofvertices = this->vertices->NumberOfVertices(); + int numberofelements = -1; + int numberofvertices = -1; IssmDouble hmax = this->amrbamg->GetBamgOpts()->hmax; - IssmDouble* maxlength = xNew(numberofelements); - IssmDouble* error_vertices = xNewZeroInit(numberofvertices); + IssmDouble* maxlength = NULL; + IssmDouble* error_vertices = NULL; IssmDouble* error_elements = NULL; IssmDouble* x = NULL; IssmDouble* y = NULL; - IssmDouble* z = NULL; int* index = NULL; IssmDouble maxerror,threshold,groupthreshold,resolution,length; @@ -5216,5 +5254,7 @@ /*Get mesh*/ - this->GetMesh(this->vertices,this->elements,&x,&y,&z,&index); + this->GetMesh(&index,&x,&y,&numberofvertices,&numberofelements); + maxlength = xNew(numberofelements); + error_vertices = xNewZeroInit(numberofvertices); /*Fill error_vertices (this is the sum of all elements connected to the vertex)*/ @@ -5260,8 +5300,4 @@ /*Cleanup*/ - xDelete(x); - xDelete(y); - xDelete(z); - xDelete(index); xDelete(error_elements); xDelete(error_vertices); @@ -5270,4 +5306,7 @@ /*}}}*/ void FemModel::GetVerticeDistanceToZeroLevelSet(IssmDouble** pverticedistance,int levelset_type){/*{{{*/ + + //itapopo esse metodo pode ser deletado + /*Here, "zero level set" means grounding line or ice front, depending on the level set type*/ @@ -5281,14 +5320,16 @@ /*Intermediaries*/ - int numberofvertices = this->vertices->NumberOfVertices(); + int numberofvertices = -1; + int numberofelements = -1; IssmDouble* levelset_points= NULL; IssmDouble* x = NULL; IssmDouble* y = NULL; - IssmDouble* z = NULL; + int* elementslist = NULL; int numberofpoints; IssmDouble distance; /*Get vertices coordinates*/ - VertexCoordinatesx(&x,&y,&z,this->vertices,false) ; + this->GetMesh(&elementslist,&x,&y,&numberofvertices,&numberofelements); + //this->GetMeshOnPartition(this->vertices,this->elements,&x,&y,&z,&elementslist,&sidtoindex); /*Get points which level set is zero (center of elements with zero level set)*/ @@ -5310,7 +5351,4 @@ /*Cleanup*/ xDelete(levelset_points); - xDelete(x); - xDelete(y); - xDelete(z); } /*}}}*/ @@ -5393,5 +5431,4 @@ IssmDouble* x = NULL; IssmDouble* y = NULL; - IssmDouble* z = NULL; int* elements = NULL; int amr_restart; @@ -5402,5 +5439,5 @@ /*Get vertices coordinates of the coarse mesh (father mesh)*/ /*elements comes in Matlab indexing*/ - this->GetMesh(this->vertices,this->elements,&x,&y,&z,&elements); + this->GetMesh(this->vertices,this->elements,&x,&y,&elements); /*Create initial mesh (coarse mesh) in neopz data structure*/ @@ -5422,4 +5459,7 @@ this->parameters->FindParam(&this->amr->deviatoricerror_groupthreshold,AmrDeviatoricErrorGroupThresholdEnum); this->parameters->FindParam(&this->amr->deviatoricerror_maximum,AmrDeviatoricErrorMaximumEnum); + + /*Initialize NeoPZ data structure*/ + this->amr->SetMesh(&elements,&x,&y,&numberofvertices,&numberofelements); if(my_rank==0){ this->parameters->FindParam(&amr_restart,AmrRestartEnum); @@ -5427,13 +5467,7 @@ this->amr->ReadMesh(); } else {//this is the default method - this->amr->CreateInitialMesh(numberofvertices,numberofelements,x,y,elements); - } - } - - /*Free the vectors*/ - xDelete(x); - xDelete(y); - xDelete(z); - xDelete(elements); + this->amr->Initialize(); + } + } } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 23500) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 23501) @@ -174,5 +174,7 @@ void BedrockFromMismipPlus(void); void AdjustBaseThicknessAndMask(void); - void GetMesh(Vertices* femmodel_vertices,Elements* femmodel_elements,IssmDouble** px, IssmDouble** py, IssmDouble** pz, int** pelementslist); + void GetMesh(Vertices* femmodel_vertices,Elements* femmodel_elements,IssmDouble** px, IssmDouble** py, int** pelementslist); + void GetMesh(int** elementslist, IssmDouble** x, IssmDouble** y, int* numberofvertices, int* numberofelements); + void SetMesh(int** elementslist, IssmDouble** x, IssmDouble** y, int* numberofvertices, int* numberofelements); void GetMeshOnPartition(Vertices* femmodel_vertices,Elements* femmodel_elements,IssmDouble** px, IssmDouble** py, IssmDouble** pz, int** pelementslist,int** psidtoindex); void GetGroundediceLevelSet(IssmDouble** pmasklevelset); Index: /issm/trunk-jpl/src/c/cores/transient_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 23500) +++ /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 23501) @@ -101,5 +101,4 @@ IssmDouble* lat_ice = NULL; IssmDouble* lon_ice = NULL; - IssmDouble* z_ice = NULL; IssmDouble* icebase = NULL; IssmDouble* icemask = NULL; @@ -143,5 +142,5 @@ /*Interpolate ice base and mask onto ocean grid*/ - femmodel->GetMesh(femmodel->vertices,femmodel->elements,&x_ice,&y_ice,&z_ice,&index_ice); + femmodel->GetMesh(femmodel->vertices,femmodel->elements,&x_ice,&y_ice,&index_ice); BamgTriangulatex(&index_ocean,&nels_ocean,oceangridx,oceangridy,ngrids_ocean); femmodel->vertices->LatLonList(&lat_ice,&lon_ice); @@ -194,5 +193,4 @@ xDelete(x_ice); xDelete(y_ice); - xDelete(z_ice); xDelete(icebase_oceangrid); xDelete(oceangridx); @@ -271,5 +269,4 @@ IssmDouble *lat_ice = NULL; IssmDouble *lon_ice = NULL; - IssmDouble *z_ice = NULL; IssmDouble *icebase = NULL; IssmDouble *icemask = NULL; @@ -282,5 +279,5 @@ /*Recover mesh and inputs needed*/ femmodel->parameters->FindParam(&rho_ice,MaterialsRhoIceEnum); - femmodel->GetMesh(femmodel->vertices,femmodel->elements,&x_ice,&y_ice,&z_ice,&index_ice); + femmodel->GetMesh(femmodel->vertices,femmodel->elements,&x_ice,&y_ice,&index_ice); femmodel->parameters->FindParam(&oceangridx,&ngrids_ocean,OceanGridXEnum); femmodel->parameters->FindParam(&oceangridy,&ngrids_ocean,OceanGridYEnum); @@ -354,5 +351,4 @@ xDelete(x_ice); xDelete(y_ice); - xDelete(z_ice); xDelete(melt_mesh); xDelete(oceangridx);