Changeset 23501


Ignore:
Timestamp:
12/02/18 15:34:43 (6 years ago)
Author:
tsantos
Message:

CHG: optimizing AMR (distance to GL parallelized, global mesh data structure is kept in AMR object, GetMesh is changed such that no z-coordinate is necessary)

Location:
issm/trunk-jpl/src/c
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp

    r23471 r23501  
    1515/*Constructor, copy, clean up and destructor*/
    1616AdaptiveMeshRefinement::AdaptiveMeshRefinement(){/*{{{*/
    17         this->Initialize();
     17
     18        /*Set pointers to NULL*/
     19        this->fathermesh                                                = NULL;
     20        this->previousmesh                                      = NULL;
     21        this->refinement_type                           = -1;
     22        this->level_max                                         = -1;
     23        this->gradation                                         = -1;
     24        this->lag                                                               = -1;
     25   this->groundingline_distance         = -1;
     26        this->icefront_distance                         = -1;
     27        this->thicknesserror_threshold  = -1;
     28        this->deviatoricerror_threshold = -1;
     29        this->deviatoricerror_maximum           = -1;
     30        this->thicknesserror_maximum            = -1;
     31        this->sid2index.clear();
     32        this->index2sid.clear();
     33        this->specialelementsindex.clear();
     34        this->x                                                                 = NULL;
     35        this->y                                                                 = NULL;
     36        this->elementslist                                      = NULL;
     37        this->numberofvertices                          = -1;
     38        this->numberofelements                          = -1;
    1839}
    1940/*}}}*/
    … …  
    5475/*}}}*/
    5576AdaptiveMeshRefinement::~AdaptiveMeshRefinement(){/*{{{*/
    56         int writemesh = 1;//itapopo keep 0
     77        int writemesh = 0;//only to restart
    5778        if(writemesh) this->WriteMesh();
    5879        this->CleanUp();
    … …  
    6586        if(this->fathermesh)    delete this->fathermesh;
    6687        if(this->previousmesh)  delete this->previousmesh;
     88        if(this->x)                                     delete this->x;
     89        if(this->y)                                     delete this->y;
     90        if(this->elementslist)  delete this->elementslist;
    6791        this->refinement_type                           = -1;
    6892        this->level_max                                         = -1;
    … …  
    7599        this->deviatoricerror_maximum           = -1;
    76100        this->thicknesserror_maximum            = -1;
     101        this->numberofvertices                          = -1;
     102        this->numberofelements                          = -1;
    77103        this->sid2index.clear();
    78104        this->index2sid.clear();
    … …  
    80106}
    81107/*}}}*/
    82 void AdaptiveMeshRefinement::Initialize(){/*{{{*/
    83 
    84         /*Set pointers to NULL*/
    85         this->fathermesh                                                = NULL;
    86         this->previousmesh                                      = NULL;
    87         this->refinement_type                           = -1;
    88         this->level_max                                         = -1;
    89         this->gradation                                         = -1;
    90         this->lag                                                               = -1;
    91    this->groundingline_distance         = -1;
    92         this->icefront_distance                         = -1;
    93         this->thicknesserror_threshold  = -1;
    94         this->deviatoricerror_threshold = -1;
    95         this->deviatoricerror_maximum           = -1;
    96         this->thicknesserror_maximum            = -1;
    97         this->sid2index.clear();
    98         this->index2sid.clear();
    99         this->specialelementsindex.clear();
    100 }
    101 /*}}}*/
    102108
    103109/*Mesh refinement methods*/
     110void AdaptiveMeshRefinement::SetMesh(int** elementslist_in,IssmDouble** x_in,IssmDouble** y_in,int* numberofvertices_in,int* numberofelements_in){/*{{{*/
     111   
     112   /*Delete previous mesh and keep the entire mesh*/
     113   if(this->elementslist) xDelete<int>(this->elementslist);
     114   if(this->x) xDelete<IssmDouble>(this->x);
     115   if(this->y) xDelete<IssmDouble>(this->y);
     116
     117   this->elementslist      = *elementslist_in;
     118   this->x                 = *x_in;
     119   this->y                 = *y_in;
     120   this->numberofvertices  = *numberofvertices_in;
     121   this->numberofelements  = *numberofelements_in;
     122}/*}}}*/
     123void AdaptiveMeshRefinement::GetMesh(int** elementslist_out,IssmDouble** x_out,IssmDouble** y_out,int* numberofvertices_out,int* numberofelements_out){/*{{{*/
     124   
     125   /*Get the entire mesh*/
     126   *elementslist_out    = this->elementslist;
     127   *x_out               = this->x;
     128   *y_out               = this->y;
     129   *numberofvertices_out= this->numberofvertices;
     130   *numberofelements_out= this->numberofelements;
     131}/*}}}*/
    104132void AdaptiveMeshRefinement::ExecuteRefinement(double* gl_distance,double* if_distance,double* deviatoricerror,double* thicknesserror,int** pdatalist,double** pxylist,int** pelementslist){/*{{{*/
    105133
    … …  
    555583}
    556584/*}}}*/
    557 void AdaptiveMeshRefinement::CreateInitialMesh(int &nvertices,int &nelements,double* x,double* y,int* elements){/*{{{*/
     585void AdaptiveMeshRefinement::Initialize(){/*{{{*/
    558586
    559587        /* IMPORTANT! elements come in Matlab indexing
    560588                NEOPZ works only in C indexing*/
    561589
    562         if(nvertices<=0) _error_("Impossible to create initial mesh: nvertices is <= 0!\n");
    563    if(nelements<=0) _error_("Impossible to create initial mesh: nelements is <= 0!\n");
     590        if(this->numberofvertices<=0) _error_("Impossible to create initial mesh: nvertices is <= 0!\n");
     591   if(this->numberofelements<=0) _error_("Impossible to create initial mesh: nelements is <= 0!\n");
    564592        if(this->refinement_type!=0 && this->refinement_type!=1) _error_("Impossible to create initial mesh: refinement type is not defined!\n");
    565593
    566594    /*Verify and creating initial mesh*/
    567    if(this->fathermesh || this->previousmesh) _error_("Initial mesh already exists!");
     595   if(!this->x || !this->y || !this->elementslist) _error_("Mesh data structure is NULL!\n");
     596        if(this->fathermesh || this->previousmesh) _error_("Initial mesh already exists!\n");
    568597
    569598   this->fathermesh = new TPZGeoMesh();
    570         this->fathermesh->NodeVec().Resize(nvertices);
     599        this->fathermesh->NodeVec().Resize(this->numberofvertices);
    571600
    572601        /*Set the vertices (geometric nodes in NeoPZ context)*/
    573         for(int i=0;i<nvertices;i++){ 
     602        for(int i=0;i<this->numberofvertices;i++){ 
    574603      /*x,y,z coords*/
    575604                TPZManVector<REAL,3> coord(3,0.);
    576       coord[0]= x[i];
    577       coord[1]= y[i];
     605      coord[0]= this->x[i];
     606      coord[1]= this->y[i];
    578607      coord[2]= 0.;
    579608      /*Insert in the mesh*/
    … …  
    586615   const int mat = this->GetElemMaterialID();
    587616   TPZManVector<long> elem(this->GetNumberOfNodes(),0);
    588         this->index2sid.clear(); this->index2sid.resize(nelements);
     617        this->index2sid.clear(); this->index2sid.resize(this->numberofelements);
    589618   this->sid2index.clear();
    590619
    591         for(int i=0;i<nelements;i++){
    592                 for(int j=0;j<this->GetNumberOfNodes();j++) elem[j]=elements[i*this->GetNumberOfNodes()+j]-1;//Convert Matlab to C indexing
     620        for(int i=0;i<this->numberofelements;i++){
     621                for(int j=0;j<this->GetNumberOfNodes();j++) elem[j]=this->elementslist[i*this->GetNumberOfNodes()+j]-1;//Convert Matlab to C indexing
    593622      switch(this->GetNumberOfNodes()){
    594623                        case 3: this->fathermesh->CreateGeoElement(ETriangle,elem,mat,index,this->refinement_type);     break;
    … …  
    865894        std::ifstream amrifstream(amrfile.c_str());                             
    866895        int size,value;
    867        
     896
    868897        TPZPersistenceManager::OpenRead(fathermeshfile);
    869898        this->fathermesh = dynamic_cast<TPZGeoMesh *>(TPZPersistenceManager::ReadFromFile());

  • issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h

    r23471 r23501  
    6868        void Initialize();
    6969   void ExecuteRefinement(double* gl_distance,double* if_distance,double* deviatoricerror,double* thicknesserror,int** pdatalist,double** pxy,int** pelementslist);
    70         void CreateInitialMesh(int &nvertices,int &nelements,double* x,double* y,int* elements);
     70        void SetMesh(int** elementslist_in,IssmDouble** x_in,IssmDouble** y_in,int* numberofvertices,int* numberofelements);
     71        void GetMesh(int** elementslist_out,IssmDouble** x_out,IssmDouble** y_out,int* numberofvertices,int* numberofelements);
    7172        void CheckMesh(int** pdata,double** pxy,int** pelements);
    7273        void ReadMesh();
    … …  
    8081        TPZGeoMesh *fathermesh;                                                 // Entire mesh without refinement if refinement_type==1; refined with hanging nodes if efinement_type==0
    8182        TPZGeoMesh *previousmesh;                                               // Refined mesh without hanging nodes (it is always refpattern type), used to generate ISSM mesh
     83        IssmDouble* x;                                                                  // Entire mesh
     84   IssmDouble* y;
     85   int* elementslist;
     86   int numberofvertices;
     87   int numberofelements;       
    8288        /*}}}*/
    8389        /*Private methods{{{*/

  • issm/trunk-jpl/src/c/classes/AmrBamg.cpp

    r23066 r23501  
    3939        this->fathermesh                                                                = NULL;
    4040        this->previousmesh                                                      = NULL;
     41        this->elementslist                                                      = NULL;
     42        this->x                                                                                 = NULL;
     43        this->y                                                                                 = NULL;
     44        this->numberofvertices                                          = -1;
     45        this->numberofelements                                          = -1;
    4146
    4247        /*Only initialize options for now (same as bamg.m)*/
    … …  
    7378        if(this->previousmesh) delete this->previousmesh;
    7479        if(this->options) delete this->options;
    75 
     80        if(this->x) xDelete<IssmDouble>(this->x);
     81        if(this->y) xDelete<IssmDouble>(this->y);
     82        if(this->elementslist) xDelete<int>(this->elementslist);
    7683}
    7784/*}}}*/
    7885
    7986/*Methods*/
    80 void AmrBamg::Initialize(int* elements,IssmDouble* x,IssmDouble* y,int numberofvertices,int numberofelements){/*{{{*/
     87void AmrBamg::SetMesh(int** elementslist_in,IssmDouble** x_in,IssmDouble** y_in,int* numberofvertices_in,int* numberofelements_in){/*{{{*/
     88       
     89        /*Delete previous mesh and keep the entire mesh*/
     90        if(this->elementslist) xDelete<int>(this->elementslist);
     91        if(this->x) xDelete<IssmDouble>(this->x);
     92        if(this->y) xDelete<IssmDouble>(this->y);
     93
     94        this->elementslist              = *elementslist_in;
     95        this->x                                         = *x_in;
     96        this->y                                         = *y_in;
     97        this->numberofvertices  = *numberofvertices_in;
     98        this->numberofelements  = *numberofelements_in;
     99}/*}}}*/
     100void AmrBamg::GetMesh(int** elementslist_out,IssmDouble** x_out,IssmDouble** y_out,int* numberofvertices_out,int* numberofelements_out){/*{{{*/
     101       
     102        /*Get the entire mesh*/
     103        *elementslist_out               = this->elementslist;
     104        *x_out                                  = this->x;
     105        *y_out                                  = this->y;
     106        *numberofvertices_out= this->numberofvertices;
     107        *numberofelements_out= this->numberofelements;
     108}/*}}}*/
     109void AmrBamg::Initialize(){/*{{{*/
    81110
    82111        /*Check options*/
    … …  
    85114
    86115        /*Read father mesh and create geometry*/
    87         Mesh* Th=new Mesh(elements,x,y,numberofvertices,numberofelements,this->options);
     116        Mesh* Th=new Mesh(this->elementslist,this->x,this->y,this->numberofvertices,this->numberofelements,this->options);
    88117
    89118        /*Write geometry*/
    … …  
    99128}/*}}}*/
    100129void AmrBamg::ExecuteRefinementBamg(IssmDouble* field,IssmDouble* hmaxVertices,int** pdatalist,IssmDouble** pxylist,int** pelementslist){/*{{{*/
    101 
     130       
    102131        /*Intermediaries*/
    103132        BamgGeom* geomout=new BamgGeom();
    … …  
    113142        this->options->field                     = field;
    114143        this->options->hmaxVertices = hmaxVertices;
    115 
     144       
    116145        /*remesh*/
    117146        if(this->previousmesh){

  • issm/trunk-jpl/src/c/classes/AmrBamg.h

    r23469 r23501  
    2828                IssmDouble deviatoricerror_maximum;
    2929
     30
    3031                /* Constructor, destructor etc*/
    3132                AmrBamg();
    … …  
    3435
    3536                /*General methods*/
    36                 void Initialize(int* elements,IssmDouble* x,IssmDouble* y,int numberofvertices,int numberofelements);
     37                void Initialize();
     38                void SetMesh(int** elementslist_in,IssmDouble** x_in,IssmDouble** y_in,int* numberofvertices,int* numberofelements);
     39                void GetMesh(int** elementslist_out,IssmDouble** x_out,IssmDouble** y_out,int* numberofvertices,int* numberofelements);
    3740                void ExecuteRefinementBamg(IssmDouble* field,IssmDouble* hmaxVertices,int** pdatalist,IssmDouble** pxylist,int** pelementslist);
    3841                void SetBamgOpts(IssmDouble hmin_in,IssmDouble hmax_in,IssmDouble err_in,IssmDouble gradation_in);
    … …  
    4649                BamgMesh* previousmesh;
    4750                BamgOpts* options;
     51                /*entire mesh*/
     52                IssmDouble* x;
     53                IssmDouble* y;
     54                int* elementslist;
     55                int numberofvertices;
     56                int numberofelements;
    4857};
    4958

  • issm/trunk-jpl/src/c/classes/FemModel.cpp

    r23500 r23501  
    27192719        /*Finally: interpolate all inputs and insert them into the new elements.*/
    27202720        this->InterpolateInputs(new_vertices,new_elements);
    2721 
     2721       
    27222722        /*Delete old structure and set new pointers*/
    27232723        delete this->vertices;          this->vertices          = new_vertices;
    … …  
    27492749        SetCurrentConfiguration(analysis_type);
    27502750
     2751        /*Set the new mesh*/
     2752        this->SetMesh(&newelementslist,&newx,&newy,&newnumberofvertices,&newnumberofelements);
     2753
    27512754        /*Cleanup*/
    2752         xDelete<IssmDouble>(newx);
    2753         xDelete<IssmDouble>(newy);
    27542755        xDelete<IssmDouble>(newz);
    2755         xDelete<int>(newelementslist);
    27562756        xDelete<int>(my_vertices);
    27572757        xDelete<bool>(my_elements);
    … …  
    29672967void FemModel::InterpolateInputs(Vertices* newfemmodel_vertices,Elements* newfemmodel_elements){/*{{{*/
    29682968
    2969         int numberofelements                    = this->elements->NumberOfElements();   //global, entire old mesh
     2969        int numberofelements                    = -1;                                                                                           //global, entire old mesh
    29702970        int newnumberofelements         = newfemmodel_elements->Size();                 //just on the new partition
    2971         int numberofvertices                    = this->vertices->NumberOfVertices();   //global, entire old mesh
     2971        int numberofvertices                    = -1;                                                                                           //global, entire old mesh
    29722972        int newnumberofvertices         = newfemmodel_vertices->Size();                 //just on the new partition
    29732973        int elementswidth                               = this->GetElementsWidth(); //just tria in this version
    … …  
    29862986        IssmDouble* x                                   = NULL;//global, entire old mesh
    29872987        IssmDouble* y                                   = NULL;//global, entire old mesh
    2988         IssmDouble* z                                   = NULL;//global, entire old mesh
    29892988        int* elementslist                               = NULL;//global, entire old mesh
    29902989        IssmDouble* newx                                = NULL;//just on the new partition
    … …  
    30002999
    30013000        /*Get the old mesh (global, entire mesh)*/
    3002         this->GetMesh(this->vertices,this->elements,&x,&y,&z,&elementslist);
     3001        this->GetMesh(&elementslist,&x,&y,&numberofvertices,&numberofelements);
    30033002
    30043003        /*Get the new mesh (just on the new partition)*/
    … …  
    30653064        xDelete<int>(P1input_enums);
    30663065        xDelete<int>(P1input_interp);
    3067         xDelete<IssmDouble>(x);
    3068         xDelete<IssmDouble>(y);
    3069         xDelete<IssmDouble>(z);
    3070         xDelete<int>(elementslist);
    30713066        xDelete<IssmDouble>(newx);
    30723067        xDelete<IssmDouble>(newy);
    … …  
    30903085        IssmDouble* x                   = NULL;
    30913086        IssmDouble* y                   = NULL;
    3092         IssmDouble* z                   = NULL;
    30933087        int* elementslist               = NULL;
    30943088
    … …  
    30973091        parameters->FindParam(&step,StepEnum);
    30983092        parameters->FindParam(&time,TimeEnum);
    3099         numberofelements=this->elements->NumberOfElements();
    3100         numberofvertices=this->vertices->NumberOfVertices();
    31013093
    31023094        /*Get mesh. Elementslist comes in Matlab indexing*/
    3103         this->GetMesh(this->vertices,this->elements,&x,&y,&z,&elementslist);
     3095        this->GetMesh(&elementslist,&x,&y,&numberofvertices,&numberofelements);
    31043096
    31053097        /*Write mesh in Results*/
    … …  
    31123104        this->results->AddResult(new GenericExternalResult<IssmDouble*>(this->results->Size()+1,MeshYEnum,
    31133105                                        y,numberofvertices,1,step,time));
    3114 
    3115         /*Cleanup*/
    3116         xDelete<IssmDouble>(x);
    3117         xDelete<IssmDouble>(y);
    3118         xDelete<IssmDouble>(z);
    3119         xDelete<int>(elementslist);
    31203106}
    31213107/*}}}*/
    … …  
    33263312}
    33273313/*}}}*/
    3328 void FemModel::GetMesh(Vertices* femmodel_vertices, Elements* femmodel_elements,IssmDouble** px, IssmDouble** py, IssmDouble** pz, int** pelementslist){/*{{{*/
     3314void FemModel::GetMesh(Vertices* femmodel_vertices, Elements* femmodel_elements,IssmDouble** px, IssmDouble** py, int** pelementslist){/*{{{*/
    33293315
    33303316        if(!femmodel_vertices) _error_("GetMesh: vertices are NULL.");
    … …  
    33833369        *px                             = x;
    33843370        *py                             = y;
    3385         *pz                             = z;
    33863371        *pelementslist = elementslist; //Matlab indexing. InterMesh uses this type.
    33873372
    … …  
    33913376        xDelete<IssmDouble>(id2);
    33923377        xDelete<IssmDouble>(id3);
     3378        xDelete<IssmDouble>(z);
    33933379        delete vid1;
    33943380        delete vid2;
    … …  
    33963382}
    33973383/*}}}*/
     3384void FemModel::GetMesh(int** elementslist, IssmDouble** x, IssmDouble** y, int* numberofvertices, int* numberofelements){/*{{{*/
     3385       
     3386        int amrtype;   
     3387        this->parameters->FindParam(&amrtype,AmrTypeEnum);
     3388   
     3389        switch(amrtype){
     3390
     3391      #if defined(_HAVE_NEOPZ_)
     3392      case AmrNeopzEnum: this->amr->GetMesh(elementslist,x,y,numberofvertices,numberofelements); break;
     3393      #endif
     3394
     3395      #if defined(_HAVE_BAMG_)
     3396      case AmrBamgEnum: this->amrbamg->GetMesh(elementslist,x,y,numberofvertices,numberofelements); break;
     3397      #endif
     3398
     3399      default: _error_("not implemented yet");
     3400   }
     3401}/*}}}*/
     3402void FemModel::SetMesh(int** elementslist, IssmDouble** x, IssmDouble** y, int* numberofvertices, int* numberofelements){/*{{{*/
     3403       
     3404        int amrtype;   
     3405        this->parameters->FindParam(&amrtype,AmrTypeEnum);
     3406   
     3407        switch(amrtype){
     3408
     3409      #if defined(_HAVE_NEOPZ_)
     3410      case AmrNeopzEnum: this->amr->SetMesh(elementslist,x,y,numberofvertices,numberofelements); break;
     3411      #endif
     3412
     3413      #if defined(_HAVE_BAMG_)
     3414      case AmrBamgEnum: this->amrbamg->SetMesh(elementslist,x,y,numberofvertices,numberofelements); break;
     3415      #endif
     3416
     3417      default: _error_("not implemented yet");
     3418   }
     3419}/*}}}*/
    33983420void FemModel::GetMeshOnPartition(Vertices* femmodel_vertices,Elements* femmodel_elements,IssmDouble** px,IssmDouble** py,IssmDouble** pz,int** pelementslist,int** psidtoindex){/*{{{*/
    33993421
    … …  
    34603482        int analysis_index = AnalysisIndex(analysis_enum);
    34613483
    3462         int numberofnodes_analysistype=  this->nodes_list[analysis_index]->NumberOfNodes();
     3484        int numberofnodes_analysistype= this->nodes_list[analysis_index]->NumberOfNodes();
    34633485        int dofpernode                                          = 2;                                                                                                            //vx and vy
    34643486        int numberofcols                                        = dofpernode*2;                                                                         //to keep dofs and flags in the vspc vector
    3465         int numberofvertices                            = this->vertices->NumberOfVertices();                   //global, entire old mesh
    3466         int numberofelements                            = this->elements->NumberOfElements();                   //global, entire old mesh
     3487        int numberofvertices                            = -1;                                                                                                           //global, entire old mesh
     3488        int numberofelements                            = -1;                                                                                                           //global, entire old mesh
    34673489        int newnumberofvertices                 = newfemmodel_vertices->Size();                                 //local, just the new partition
    34683490        int count                                                       = 0;
    34693491        IssmDouble* x                                           = NULL;                                                                                                 //global, entire old mesh
    34703492        IssmDouble* y                                           = NULL;                                                                                                 //global, entire old mesh
    3471         IssmDouble* z                                           = NULL;                                                                                                 //global, entire old mesh
    34723493        int*                    elementslist            = NULL;                                                                                                 //global, entire old mesh
    34733494        IssmDouble* spc                                 = NULL;                                                                                                 //global, entire old mesh
    … …  
    34793500
    34803501        /*Get old mesh (global, entire mesh). Elementslist comes in Matlab indexing*/
    3481         this->GetMesh(this->vertices,this->elements,&x,&y,&z,&elementslist);
     3502        this->GetMesh(&elementslist,&x,&y,&numberofvertices,&numberofelements);
    34823503
    34833504        /*Get vertices coordinates of the new partition*/
    … …  
    35453566                }
    35463567        }
    3547 
     3568       
    35483569        /*Cleanup*/
    3549         xDelete<IssmDouble>(x);
    3550         xDelete<IssmDouble>(y);
    3551         xDelete<IssmDouble>(z);
    3552         xDelete<int>(elementslist);
    35533570        xDelete<IssmDouble>(spc);
    35543571        xDelete<IssmDouble>(newspc);
    … …  
    50775094        IssmDouble* x             = NULL;
    50785095        IssmDouble* y             = NULL;
    5079         IssmDouble* z             = NULL;
    50805096        int* elements             = NULL;
    50815097        IssmDouble hmin,hmax,err,gradation;
    … …  
    50885104
    50895105        /*Get vertices coordinates of the coarse mesh (father mesh)*/
    5090         this->GetMesh(this->vertices,this->elements,&x,&y,&z,&elements);
     5106        this->GetMesh(this->vertices,this->elements,&x,&y,&elements);
    50915107
    50925108        /*Create bamg data structures for bamg*/
    … …  
    51165132
    51175133        /*Re-create original mesh and put it in bamg structure (only cpu 0)*/
     5134        this->amrbamg->SetMesh(&elements,&x,&y,&numberofvertices,&numberofelements);
    51185135        if(my_rank==0){
    5119                 this->amrbamg->Initialize(elements,x,y,numberofvertices,numberofelements);
    5120         }
    5121 
    5122         /*Free the vectors*/
    5123         xDelete<IssmDouble>(x);
    5124         xDelete<IssmDouble>(y);
    5125         xDelete<IssmDouble>(z);
    5126         xDelete<int>(elements);
     5136                this->amrbamg->Initialize();
     5137        }
    51275138}
    51285139/*}}}*/
    … …  
    51325143
    51335144        /*Intermediaries*/
    5134         int numberofvertices                     = this->vertices->NumberOfVertices();
    5135         IssmDouble* verticedistance = NULL;
     5145        int numberofvertices                                                    = this->vertices->NumberOfVertices();
     5146   Vector<IssmDouble>* vminvertexdistance = new Vector<IssmDouble>(numberofvertices);
     5147        IssmDouble* pminvertexdistance                  = NULL;
     5148        IssmDouble* levelset_points                             = NULL;
     5149        IssmDouble x,y;
    51365150        IssmDouble threshold,resolution;
     5151        IssmDouble minvertexdistance,distance;
     5152        int sid,numberofpoints;
    51375153
    51385154        switch(levelset_type){
    … …  
    51485164        }
    51495165
    5150         /*Get vertice distance to zero level set points*/
    5151         this->GetVerticeDistanceToZeroLevelSet(&verticedistance,levelset_type);
    5152         if(!verticedistance) _error_("verticedistance is NULL!\n");
     5166        /*Get points which level set is zero (center of elements with zero level set)*/
     5167   this->GetZeroLevelSetPoints(&levelset_points,numberofpoints,levelset_type);//levelset_points is serial (global)
     5168
     5169        for(int i=0;i<this->vertices->Size();i++){//only on this partition
     5170                Vertex* vertex=(Vertex*)this->vertices->GetObjectByOffset(i);
     5171      /*Attention: no spherical coordinates*/
     5172      x = vertex->GetX();
     5173      y         = vertex->GetY();
     5174      sid= vertex->Sid();
     5175                minvertexdistance=INFINITY;
     5176     
     5177                /*Find the minimum vertex distance*/
     5178                for(int j=0;j<numberofpoints;j++){
     5179         distance=sqrt((x-levelset_points[2*j])*(x-levelset_points[2*j])+(y-levelset_points[2*j+1])*(y-levelset_points[2*j+1]));
     5180         minvertexdistance=min(distance,minvertexdistance);
     5181        }
     5182                /*Now, insert in the vector*/
     5183                vminvertexdistance->SetValue(sid,minvertexdistance,INS_VAL);   
     5184   }
     5185        /*Assemble*/
     5186   vminvertexdistance->Assemble();
     5187
     5188   /*Assign the pointer*/
     5189        pminvertexdistance=vminvertexdistance->ToMPISerial();
    51535190
    51545191        /*Fill hmaxVertices*/
    51555192        for(int i=0;i<numberofvertices;i++){
    5156                 if(verticedistance[i]<threshold){
     5193                if(pminvertexdistance[i]<threshold){ //hmaxvertices is serial (global)
    51575194                        if(xIsNan<IssmDouble>(hmaxvertices[i])) hmaxvertices[i]=resolution;
    51585195                        else hmaxvertices[i]=min(resolution,hmaxvertices[i]);
    … …  
    51615198
    51625199        /*Cleanup*/
    5163         xDelete<IssmDouble>(verticedistance);
     5200        xDelete<IssmDouble>(pminvertexdistance);
     5201   xDelete<IssmDouble>(levelset_points);
     5202   delete vminvertexdistance;
    51645203}
    51655204/*}}}*/
    … …  
    51705209        /*Intermediaries*/
    51715210        int elementswidth                                                       = this->GetElementsWidth();
    5172         int numberofelements                                            = this->elements->NumberOfElements();
    5173         int numberofvertices                                            = this->vertices->NumberOfVertices();
     5211        int numberofelements                                            = -1;
     5212        int numberofvertices                                            = -1;
    51745213        IssmDouble      hmax                                                    = this->amrbamg->GetBamgOpts()->hmax;
    5175         IssmDouble* maxlength                                   = xNew<IssmDouble>(numberofelements);
    5176         IssmDouble* error_vertices                              = xNewZeroInit<IssmDouble>(numberofvertices);
     5214        IssmDouble* maxlength                                   = NULL;
     5215        IssmDouble* error_vertices                              = NULL;
    51775216        IssmDouble* error_elements                              = NULL;
    51785217        IssmDouble* x                                                           = NULL;
    51795218        IssmDouble* y                                                           = NULL;
    5180         IssmDouble* z                                                           = NULL;
    51815219        int* index                                                                      = NULL;
    51825220        IssmDouble maxerror,threshold,groupthreshold,resolution,length;
    … …  
    52165254
    52175255        /*Get mesh*/
    5218         this->GetMesh(this->vertices,this->elements,&x,&y,&z,&index);
     5256        this->GetMesh(&index,&x,&y,&numberofvertices,&numberofelements);
     5257        maxlength               = xNew<IssmDouble>(numberofelements);
     5258        error_vertices  = xNewZeroInit<IssmDouble>(numberofvertices);
    52195259
    52205260        /*Fill error_vertices (this is the sum of all elements connected to the vertex)*/
    … …  
    52605300
    52615301        /*Cleanup*/
    5262         xDelete<IssmDouble>(x);
    5263         xDelete<IssmDouble>(y);
    5264         xDelete<IssmDouble>(z);
    5265         xDelete<int>(index);
    52665302   xDelete<IssmDouble>(error_elements);
    52675303   xDelete<IssmDouble>(error_vertices);
    … …  
    52705306/*}}}*/
    52715307void FemModel::GetVerticeDistanceToZeroLevelSet(IssmDouble** pverticedistance,int levelset_type){/*{{{*/
     5308       
     5309        //itapopo esse metodo pode ser deletado
     5310
    52725311
    52735312        /*Here, "zero level set" means grounding line or ice front, depending on the level set type*/
    … …  
    52815320
    52825321        /*Intermediaries*/
    5283    int numberofvertices       = this->vertices->NumberOfVertices();
     5322   int numberofvertices       = -1;
     5323        int numberofelements                    = -1;
    52845324   IssmDouble* levelset_points= NULL;
    52855325   IssmDouble* x                                        = NULL;
    52865326   IssmDouble* y                                        = NULL;
    5287    IssmDouble* z                                        = NULL;
     5327        int* elementslist                               = NULL;
    52885328        int numberofpoints;
    52895329        IssmDouble distance;
    52905330
    52915331        /*Get vertices coordinates*/
    5292         VertexCoordinatesx(&x,&y,&z,this->vertices,false) ;
     5332        this->GetMesh(&elementslist,&x,&y,&numberofvertices,&numberofelements);
     5333        //this->GetMeshOnPartition(this->vertices,this->elements,&x,&y,&z,&elementslist,&sidtoindex);
    52935334
    52945335        /*Get points which level set is zero (center of elements with zero level set)*/
    … …  
    53105351        /*Cleanup*/
    53115352   xDelete<IssmDouble>(levelset_points);
    5312    xDelete<IssmDouble>(x);
    5313    xDelete<IssmDouble>(y);
    5314    xDelete<IssmDouble>(z);
    53155353}
    53165354/*}}}*/
    … …  
    53935431        IssmDouble* x                                                                   = NULL;
    53945432        IssmDouble* y                                                                   = NULL;
    5395         IssmDouble* z                                                                   = NULL;
    53965433        int* elements                                                                   = NULL;
    53975434        int amr_restart;
    … …  
    54025439        /*Get vertices coordinates of the coarse mesh (father mesh)*/
    54035440        /*elements comes in Matlab indexing*/
    5404         this->GetMesh(this->vertices,this->elements,&x,&y,&z,&elements);
     5441        this->GetMesh(this->vertices,this->elements,&x,&y,&elements);
    54055442
    54065443        /*Create initial mesh (coarse mesh) in neopz data structure*/
    … …  
    54225459        this->parameters->FindParam(&this->amr->deviatoricerror_groupthreshold,AmrDeviatoricErrorGroupThresholdEnum);
    54235460        this->parameters->FindParam(&this->amr->deviatoricerror_maximum,AmrDeviatoricErrorMaximumEnum);
     5461       
     5462        /*Initialize NeoPZ data structure*/
     5463        this->amr->SetMesh(&elements,&x,&y,&numberofvertices,&numberofelements);
    54245464        if(my_rank==0){
    54255465                this->parameters->FindParam(&amr_restart,AmrRestartEnum);
    … …  
    54275467                        this->amr->ReadMesh();
    54285468                } else {//this is the default method
    5429                         this->amr->CreateInitialMesh(numberofvertices,numberofelements,x,y,elements);
    5430                 }
    5431         }
    5432 
    5433         /*Free the vectors*/
    5434         xDelete<IssmDouble>(x);
    5435         xDelete<IssmDouble>(y);
    5436         xDelete<IssmDouble>(z);
    5437         xDelete<int>(elements);
     5469                        this->amr->Initialize();
     5470                }
     5471        }
    54385472}
    54395473/*}}}*/

  • issm/trunk-jpl/src/c/classes/FemModel.h

    r23493 r23501  
    174174                void BedrockFromMismipPlus(void);
    175175                void AdjustBaseThicknessAndMask(void);
    176                 void GetMesh(Vertices* femmodel_vertices,Elements* femmodel_elements,IssmDouble** px, IssmDouble** py, IssmDouble** pz, int** pelementslist);
     176                void GetMesh(Vertices* femmodel_vertices,Elements* femmodel_elements,IssmDouble** px, IssmDouble** py, int** pelementslist);
     177                void GetMesh(int** elementslist, IssmDouble** x, IssmDouble** y, int* numberofvertices, int* numberofelements);
     178                void SetMesh(int** elementslist, IssmDouble** x, IssmDouble** y, int* numberofvertices, int* numberofelements);
    177179                void GetMeshOnPartition(Vertices* femmodel_vertices,Elements* femmodel_elements,IssmDouble** px, IssmDouble** py, IssmDouble** pz, int** pelementslist,int** psidtoindex);
    178180                void GetGroundediceLevelSet(IssmDouble** pmasklevelset);

  • issm/trunk-jpl/src/c/cores/transient_core.cpp

    r23485 r23501  
    101101                IssmDouble* lat_ice           = NULL;
    102102                IssmDouble* lon_ice           = NULL;
    103                 IssmDouble* z_ice             = NULL;
    104103                IssmDouble* icebase           = NULL;
    105104                IssmDouble* icemask           = NULL;
    … …  
    143142
    144143                /*Interpolate ice base and mask onto ocean grid*/
    145                 femmodel->GetMesh(femmodel->vertices,femmodel->elements,&x_ice,&y_ice,&z_ice,&index_ice);
     144                femmodel->GetMesh(femmodel->vertices,femmodel->elements,&x_ice,&y_ice,&index_ice);
    146145                BamgTriangulatex(&index_ocean,&nels_ocean,oceangridx,oceangridy,ngrids_ocean);
    147146                femmodel->vertices->LatLonList(&lat_ice,&lon_ice);
    … …  
    194193                xDelete<IssmDouble>(x_ice);
    195194                xDelete<IssmDouble>(y_ice);
    196                 xDelete<IssmDouble>(z_ice);
    197195                xDelete<IssmDouble>(icebase_oceangrid);
    198196                xDelete<IssmDouble>(oceangridx);
    … …  
    271269                        IssmDouble *lat_ice           = NULL;
    272270                        IssmDouble *lon_ice           = NULL;
    273                         IssmDouble *z_ice             = NULL;
    274271                        IssmDouble *icebase           = NULL;
    275272                        IssmDouble *icemask           = NULL;
    … …  
    282279                        /*Recover mesh and inputs needed*/
    283280                        femmodel->parameters->FindParam(&rho_ice,MaterialsRhoIceEnum);
    284                         femmodel->GetMesh(femmodel->vertices,femmodel->elements,&x_ice,&y_ice,&z_ice,&index_ice);
     281                        femmodel->GetMesh(femmodel->vertices,femmodel->elements,&x_ice,&y_ice,&index_ice);
    285282                        femmodel->parameters->FindParam(&oceangridx,&ngrids_ocean,OceanGridXEnum);
    286283                        femmodel->parameters->FindParam(&oceangridy,&ngrids_ocean,OceanGridYEnum);
    … …  
    354351                        xDelete<IssmDouble>(x_ice);
    355352                        xDelete<IssmDouble>(y_ice);
    356                         xDelete<IssmDouble>(z_ice);
    357353                        xDelete<IssmDouble>(melt_mesh);
    358354                        xDelete<IssmDouble>(oceangridx);
Note: See TracChangeset for help on using the changeset viewer.