Changeset 10141


Ignore:
Timestamp:
10/07/11 15:34:54 (14 years ago)
Author:
Mathieu Morlighem
Message:

removed many GetValuePtr in elements

Location:
issm/trunk/src
Files:
1 added
2 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/objects/Elements/Penta.cpp

    r10135 r10141  
    17451745        const int  numdof2d = NDOF1*NUMVERTICES2D;
    17461746
    1747         int    i,dummy,hydroadjustment;
     1747        int    i,hydroadjustment;
    17481748        int*   doflist = NULL;
    1749         double values[numdof];
    17501749        double rho_ice,rho_water;
    1751         double bed[numdof];
    1752         double surface[numdof];
    1753         double *bed_ptr = NULL;
    1754         double *thickness_ptr = NULL;
    1755         double *surface_ptr = NULL;
     1750        double newthickness[numdof];
     1751        double newbed[numdof];
     1752        double newsurface[numdof];
     1753        double oldbed[NUMVERTICES];
     1754        double oldsurface[NUMVERTICES];
     1755        double oldthickness[NUMVERTICES];
    17561756        Penta  *penta   = NULL;
    17571757
    … …  
    17641764        /*Use the dof list to index into the solution vector and extrude it */
    17651765        for(i=0;i<numdof2d;i++){
    1766                 values[i]         =solution[doflist[i]];
    1767                 if(isnan(values[i])) _error_("NaN found in solution vector");
     1766                newthickness[i]=solution[doflist[i]];
     1767                if(isnan(newthickness[i])) _error_("NaN found in solution vector");
    17681768                /*Constrain thickness to be at least 1m*/
    1769                 if(values[i]<1) values[i]=1;
    1770                 values[i+numdof2d]=values[i];
     1769                if(newthickness[i]<1) newthickness[i]=1;
     1770                newthickness[i+numdof2d]=newthickness[i];
    17711771        }
    17721772
    17731773        /*Get previous bed, thickness and surface*/
    1774         Input* bed_input=inputs->GetInput(BedEnum);             _assert_(bed_input);
    1775         Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
    1776         Input* surface_input=inputs->GetInput(SurfaceEnum);     _assert_(surface_input);
    1777         bed_input->GetValuesPtr(&bed_ptr,&dummy);
    1778         thickness_input->GetValuesPtr(&thickness_ptr,&dummy);
    1779         surface_input->GetValuesPtr(&surface_ptr,&dummy);
     1774        GetInputListOnVertices(&oldbed[0],BedEnum);
     1775        GetInputListOnVertices(&oldsurface[0],SurfaceEnum);
     1776        GetInputListOnVertices(&oldthickness[0],ThicknessEnum);
    17801777
    17811778        /*Fing PrognosticHydrostaticAdjustment to figure out how to update the geometry:*/
    … …  
    17891786                /*If shelf: hydrostatic equilibrium*/
    17901787                if (this->nodes[i]->IsOnSheet()){
    1791                         surface[i]=bed_ptr[i]+values[i]; //surface=oldbed+newthickness
    1792                         bed[i]=bed_ptr[i]; //bed does not change
     1788                        newsurface[i]=oldbed[i]+newthickness[i]; //surface = oldbed + newthickness
     1789                        newbed[i]=oldbed[i];               //same bed: do nothing
    17931790                }
    17941791                else{ //so it is an ice shelf
    17951792                        if(hydroadjustment==AbsoluteEnum){
    1796                                         surface[i]=values[i]*(1-rho_ice/rho_water);
    1797                                         bed[i]=values[i]*(-rho_ice/rho_water);
     1793                                newsurface[i]=newthickness[i]*(1-rho_ice/rho_water);
     1794                                newbed[i]=newthickness[i]*(-rho_ice/rho_water);
    17981795                        }
    17991796                        else if(hydroadjustment==IncrementalEnum){
    1800                                         surface[i]=surface_ptr[i]+(1.0-rho_ice/rho_water)*(values[i]-thickness_ptr[i]); //surface = oldsurface + (1-di) * dH
    1801                                         bed[i]=bed_ptr[i]-rho_ice/rho_water*(values[i]-thickness_ptr[i]); //bed = oldbed + di * dH
     1797                                newsurface[i]=oldsurface[i]+(1.0-rho_ice/rho_water)*(newthickness[i]-oldthickness[i]); //surface = oldsurface + (1-di) * dH
     1798                                newbed[i]=oldbed[i]-rho_ice/rho_water*(newthickness[i]-oldthickness[i]); //bed = oldbed + di * dH
    18021799                        }
    18031800                        else _error_("Hydrostatic adjustment %i (%s) not supported yet",hydroadjustment,EnumToStringx(hydroadjustment));
    … …  
    18091806        for(;;){
    18101807                /*Add input to the element: */
    1811                 penta->inputs->AddInput(new PentaVertexInput(ThicknessEnum,values));
    1812                 penta->inputs->AddInput(new PentaVertexInput(SurfaceEnum,surface));
    1813                 penta->inputs->AddInput(new PentaVertexInput(BedEnum,bed));
     1808                penta->inputs->AddInput(new PentaVertexInput(ThicknessEnum,newthickness));
     1809                penta->inputs->AddInput(new PentaVertexInput(SurfaceEnum,newsurface));
     1810                penta->inputs->AddInput(new PentaVertexInput(BedEnum,newbed));
    18141811
    18151812                /*Stop if we have reached the surface*/
    … …  
    70977094        const int    numdof=NDOF2*NUMVERTICES;
    70987095
    7099         int     i,dummy;
     7096        int     i;
    71007097        double  rho_ice,g;
    71017098        double  values[numdof];
    … …  
    71087105        double  xyz_list[NUMVERTICES][3];
    71097106        int    *doflist = NULL;
    7110         double *vz_ptr  = NULL;
    71117107        Penta  *penta   = NULL;
    71127108
    … …  
    71417137
    71427138                /*Now Compute vel*/
    7143                 Input* vz_input=inputs->GetInput(VzEnum);
    7144                 if (vz_input){
    7145                         if (vz_input->ObjectEnum()!=PentaVertexInputEnum) _error_("Cannot compute Vel as Vz is of type %s",EnumToStringx(vz_input->ObjectEnum()));
    7146                         vz_input->GetValuesPtr(&vz_ptr,&dummy);
    7147                         for(i=0;i<NUMVERTICES;i++) vz[i]=vz_ptr[i];
    7148                 }
    7149                 else{for(i=0;i<NUMVERTICES;i++) vz[i]=0.0;}
     7139                GetInputListOnVertices(&vz[0],VzEnum,0.0); //default is 0
    71507140                for(i=0;i<NUMVERTICES;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
    71517141
    … …  
    71837173        const int    numdof2d=NDOF2*NUMVERTICES2D;
    71847174
    7185         int     i,dummy;
     7175        int     i;
    71867176        double  rho_ice,g;
    71877177        double  macayeal_values[numdof];
    … …  
    71967186        int*    doflistp = NULL;
    71977187        int*    doflistm = NULL;
    7198         double  *vz_ptr  = NULL;
    71997188        Penta   *penta   = NULL;
    72007189
    … …  
    72307219        }
    72317220
    7232         /*Get Vz*/
    7233         Input* vz_input=inputs->GetInput(VzEnum);
    7234         if (vz_input){
    7235                 if (vz_input->ObjectEnum()!=PentaVertexInputEnum){
    7236                         _error_("Cannot compute Vel as Vz is of type %s",EnumToStringx(vz_input->ObjectEnum()));
    7237                 }
    7238                 vz_input->GetValuesPtr(&vz_ptr,&dummy);
    7239                 for(i=0;i<NUMVERTICES;i++) vz[i]=vz_ptr[i];
    7240         }
    7241         else{
    7242                 for(i=0;i<NUMVERTICES;i++) vz[i]=0.0;
    7243         }
    7244 
    72457221        /*Now Compute vel*/
     7222        GetInputListOnVertices(&vz[0],VzEnum,0.0); //default is 0
    72467223        for(i=0;i<NUMVERTICES;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
    72477224
    … …  
    72777254        const int    numdof2d=NDOF2*NUMVERTICES2D;
    72787255
    7279         int     i,dummy;
     7256        int     i;
    72807257        double  stokesreconditioning;
    72817258        double  macayeal_values[numdofm];
    … …  
    72917268        int*    doflistm        = NULL;
    72927269        int*    doflists        = NULL;
    7293         double  *vzmacayeal_ptr = NULL;
    72947270        Penta   *penta          = NULL;
    72957271
    … …  
    73357311                        _error_("Cannot compute Vel as VzMacAyeal is of type %s",EnumToStringx(vzmacayeal_input->ObjectEnum()));
    73367312                }
    7337                 vzmacayeal_input->GetValuesPtr(&vzmacayeal_ptr,&dummy);
    7338                 for(i=0;i<NUMVERTICES;i++) vzmacayeal[i]=vzmacayeal_ptr[i];
     7313                GetInputListOnVertices(&vzmacayeal[0],VzMacAyealEnum);
    73397314        }
    73407315        else{
    … …  
    73737348        const int    numdof=NDOF2*NUMVERTICES;
    73747349
    7375         int    i,dummy;
     7350        int    i;
    73767351        double rho_ice,g;
    73777352        double values[numdof];
    … …  
    74457420        const int    numdofs=NDOF4*NUMVERTICES;
    74467421
    7447         int    i,dummy;
     7422        int    i;
    74487423        double pattyn_values[numdofp];
    74497424        double stokes_values[numdofs];
    … …  
    74597434        int*   doflistp      = NULL;
    74607435        int*   doflists      = NULL;
    7461         double *vzpattyn_ptr = NULL;
    74627436        Penta  *penta        = NULL;
    74637437
    … …  
    74987472                        _error_("Cannot compute Vel as VzPattyn is of type %s",EnumToStringx(vzpattyn_input->ObjectEnum()));
    74997473                }
    7500                 vzpattyn_input->GetValuesPtr(&vzpattyn_ptr,&dummy);
    7501                 for(i=0;i<NUMVERTICES;i++) vzpattyn[i]=vzpattyn_ptr[i];
     7474                GetInputListOnVertices(&vzpattyn[0],VzPattynEnum);
    75027475        }
    75037476        else{
    … …  
    75367509        const int    numdof=NDOF2*NUMVERTICES;
    75377510
    7538         int     i,dummy;
     7511        int     i;
    75397512        double  rho_ice,g;
    75407513        double  values[numdof];
    … …  
    75477520        double  xyz_list[NUMVERTICES][3];
    75487521        int*    doflist = NULL;
    7549         double* vz_ptr  = NULL;
    75507522
    75517523        /*Get dof list: */
    … …  
    75687540        }
    75697541
    7570         /*Get Vz*/
    7571         Input* vz_input=inputs->GetInput(VzEnum);
    7572         if (vz_input){
    7573                 if (vz_input->ObjectEnum()!=PentaVertexInputEnum){
    7574                         _error_("Cannot compute Vel as Vz is of type %s",EnumToStringx(vz_input->ObjectEnum()));
    7575                 }
    7576                 vz_input->GetValuesPtr(&vz_ptr,&dummy);
    7577                 for(i=0;i<NUMVERTICES;i++) vz[i]=vz_ptr[i];
    7578         }
    7579         else{
    7580                 for(i=0;i<NUMVERTICES;i++) vz[i]=0.0;
    7581         }
    7582 
    75837542        /*Now Compute vel*/
     7543        GetInputListOnVertices(&vz[0],VzEnum,0.0); //default is 0
    75847544        for(i=0;i<NUMVERTICES;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
    75857545
    … …  
    76007560        this->inputs->AddInput(new PentaVertexInput(VxEnum,vx));
    76017561        this->inputs->AddInput(new PentaVertexInput(VyEnum,vy));
    7602         this->inputs->AddInput(new TriaVertexInput(VelEnum,vel));
     7562        this->inputs->AddInput(new PentaVertexInput(VelEnum,vel));
    76037563        this->inputs->AddInput(new PentaVertexInput(PressureEnum,pressure));
    76047564
    … …  
    76127572        const int numdof=NDOF1*NUMVERTICES;
    76137573       
    7614         int      i,dummy;
     7574        int      i;
    76157575        int      approximation;
    76167576        double   rho_ice,g;
    … …  
    76277587        double   xyz_list[NUMVERTICES][3];
    76287588        int*     doflist      = NULL;
    7629         double*  vx_ptr       = NULL;
    7630         double*  vy_ptr       = NULL;
    7631         double*  vzstokes_ptr = NULL;
    76327589
    76337590
    … …  
    76527609
    76537610        /*Get Vx and Vy*/
    7654         Input* vx_input=inputs->GetInput(VxEnum);
    7655         if (vx_input){
    7656                 if (vx_input->ObjectEnum()!=PentaVertexInputEnum) _error_("Cannot compute Vel as Vx is of type %s",EnumToStringx(vx_input->ObjectEnum()));
    7657                 vx_input->GetValuesPtr(&vx_ptr,&dummy);
    7658                 for(i=0;i<NUMVERTICES;i++) vx[i]=vx_ptr[i];
    7659         }
    7660         else for(i=0;i<NUMVERTICES;i++) vx[i]=0.0;
    7661 
    7662         Input* vy_input=inputs->GetInput(VyEnum);
    7663         if (vy_input){
    7664                 if (vy_input->ObjectEnum()!=PentaVertexInputEnum) _error_("Cannot compute Vel as Vy is of type %s",EnumToStringx(vy_input->ObjectEnum()));
    7665                 vy_input->GetValuesPtr(&vy_ptr,&dummy);
    7666                 for(i=0;i<NUMVERTICES;i++) vy[i]=vy_ptr[i];
    7667         }
    7668         else for(i=0;i<NUMVERTICES;i++) vy[i]=0.0;
     7611        GetInputListOnVertices(&vx[0],VxEnum,0.0); //default is 0
     7612        GetInputListOnVertices(&vy[0],VyEnum,0.0); //default is 0
    76697613
    76707614        /*Do some modifications if we actually have a PattynStokes or MacAyealStokes element*/
    … …  
    76727616                Input* vzstokes_input=inputs->GetInput(VzStokesEnum);
    76737617                if (vzstokes_input){
    7674                         if (vzstokes_input->ObjectEnum()!=PentaVertexInputEnum) _error_("Cannot compute Vel as VzStokes is of type %s",EnumToStringx(vy_input->ObjectEnum()));
    7675                         vzstokes_input->GetValuesPtr(&vzstokes_ptr,&dummy);
    7676                         for(i=0;i<NUMVERTICES;i++) vzstokes[i]=vzstokes_ptr[i];
     7618                        if (vzstokes_input->ObjectEnum()!=PentaVertexInputEnum) _error_("Cannot compute Vel as VzStokes is of type %s",EnumToStringx(vzstokes_input->ObjectEnum()));
     7619                        GetInputListOnVertices(&vzstokes[0],VzStokesEnum);
    76777620                }
    76787621                else _error_("Cannot compute Vz as VzStokes in not present in PattynStokes element");
    … …  
    76857628                Input* vzstokes_input=inputs->GetInput(VzStokesEnum);
    76867629                if (vzstokes_input){
    7687                         if (vzstokes_input->ObjectEnum()!=PentaVertexInputEnum) _error_("Cannot compute Vel as VzStokes is of type %s",EnumToStringx(vy_input->ObjectEnum()));
    7688                         vzstokes_input->GetValuesPtr(&vzstokes_ptr,&dummy);
    7689                         for(i=0;i<NUMVERTICES;i++) vzstokes[i]=vzstokes_ptr[i];
     7630                        if (vzstokes_input->ObjectEnum()!=PentaVertexInputEnum) _error_("Cannot compute Vel as VzStokes is of type %s",EnumToStringx(vzstokes_input->ObjectEnum()));
     7631                        GetInputListOnVertices(&vzstokes[0],VzStokesEnum);
    76907632                }
    76917633                else _error_("Cannot compute Vz as VzStokes in not present in MacAyealStokes element");
    … …  
    77227664                this->inputs->AddInput(new PentaVertexInput(VzMacAyealEnum,vzmacayeal));
    77237665        }
    7724 
    77257666        this->inputs->AddInput(new PentaVertexInput(VzEnum,vz));
    77267667        this->inputs->AddInput(new PentaVertexInput(VelEnum,vel));

  • issm/trunk/src/c/objects/Elements/Tria.cpp

    r10135 r10141  
    333333
    334334        /*Surface and bed are updated. Update inputs:*/
    335         surface_input->GetValuesPtr(&values,NULL); for(i=0;i<3;i++)values[i]=s[i];
    336         bed_input->GetValuesPtr(&values,NULL);     for(i=0;i<3;i++)values[i]=b[i];
    337 
     335        this->inputs->AddInput(new TriaVertexInput(SurfaceEnum,&s[0]));
     336        this->inputs->AddInput(new TriaVertexInput(BedEnum,&b[0]));
    338337        }
    339338/*}}}*/
    … …  
    16511650        const int numdof = NDOF1*NUMVERTICES;
    16521651
    1653         int       i,dummy,hydroadjustment;
     1652        int       i,hydroadjustment;
    16541653        int*      doflist=NULL;
    16551654        double    rho_ice,rho_water;
    1656         double    values[numdof];
    1657         double    bed[numdof];
    1658         double    surface[numdof];
    1659         double    *bed_ptr = NULL;
    1660         double    *thickness_ptr = NULL;
    1661         double    *surface_ptr = NULL;
     1655        double    newthickness[numdof];
     1656        double    newbed[numdof];
     1657        double    newsurface[numdof];
     1658        double    oldbed[NUMVERTICES];
     1659        double    oldsurface[NUMVERTICES];
     1660        double    oldthickness[NUMVERTICES];
    16621661
    16631662        /*Get dof list: */
    … …  
    16661665        /*Use the dof list to index into the solution vector: */
    16671666        for(i=0;i<numdof;i++){
    1668                 values[i]=solution[doflist[i]];
    1669                 if(isnan(values[i])) _error_("NaN found in solution vector");
     1667                newthickness[i]=solution[doflist[i]];
     1668                if(isnan(newthickness[i])) _error_("NaN found in solution vector");
    16701669                /*Constrain thickness to be at least 1m*/
    1671                 if(values[i]<1) values[i]=1;
     1670                if(newthickness[i]<1) newthickness[i]=1;
    16721671        }
    16731672
    16741673        /*Get previous bed, thickness and surface*/
    1675         Input* bed_input=inputs->GetInput(BedEnum);             _assert_(bed_input);
    1676         Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
    1677         Input* surface_input=inputs->GetInput(SurfaceEnum);     _assert_(surface_input);
    1678         bed_input->GetValuesPtr(&bed_ptr,&dummy);
    1679         thickness_input->GetValuesPtr(&thickness_ptr,&dummy);
    1680         surface_input->GetValuesPtr(&surface_ptr,&dummy);
     1674        GetInputListOnVertices(&oldbed[0],BedEnum);
     1675        GetInputListOnVertices(&oldsurface[0],SurfaceEnum);
     1676        GetInputListOnVertices(&oldthickness[0],ThicknessEnum);
    16811677
    16821678        /*Fing PrognosticHydrostaticAdjustment to figure out how to update the geometry:*/
    16831679        this->parameters->FindParam(&hydroadjustment,PrognosticHydrostaticAdjustmentEnum);
    1684 
    1685         /*recover material parameters: */
    16861680        rho_ice=matpar->GetRhoIce();
    16871681        rho_water=matpar->GetRhoWater();
    … …  
    16901684                /*If shelf: hydrostatic equilibrium*/
    16911685                if (this->nodes[i]->IsOnSheet()){
    1692                         surface[i]=bed_ptr[i]+values[i]; //surface=oldbed+newthickness
    1693                         bed[i]=bed_ptr[i]; //do nothing
     1686                        newsurface[i]=oldbed[i]+newthickness[i]; //surface = oldbed + newthickness
     1687                        newbed[i]=oldbed[i];               //same bed: do nothing
    16941688                }
    16951689                else{ //this is an ice shelf
    16961690
    16971691                        if(hydroadjustment==AbsoluteEnum){
    1698                                 surface[i]=values[i]*(1-rho_ice/rho_water);
    1699                                 bed[i]=values[i]*(-rho_ice/rho_water);
     1692                                newsurface[i]=newthickness[i]*(1-rho_ice/rho_water);
     1693                                newbed[i]=newthickness[i]*(-rho_ice/rho_water);
    17001694                        }
    17011695                        else if(hydroadjustment==IncrementalEnum){
    1702                                 surface[i]=surface_ptr[i]+(1.0-rho_ice/rho_water)*(values[i]-thickness_ptr[i]); //surface = oldsurface + (1-di) * dH
    1703                                 bed[i]=bed_ptr[i]-rho_ice/rho_water*(values[i]-thickness_ptr[i]); //bed = oldbed + di * dH
     1696                                newsurface[i]=oldsurface[i]+(1.0-rho_ice/rho_water)*(newthickness[i]-oldthickness[i]); //surface = oldsurface + (1-di) * dH
     1697                                newbed[i]=oldbed[i]-rho_ice/rho_water*(newthickness[i]-oldthickness[i]); //bed = oldbed + di * dH
    17041698                        }
    17051699                        else _error_("Hydrostatic adjustment %i (%s) not supported yet",hydroadjustment,EnumToStringx(hydroadjustment));
    … …  
    17081702
    17091703        /*Add input to the element: */
    1710         this->inputs->AddInput(new TriaVertexInput(ThicknessEnum,values));
    1711         this->inputs->AddInput(new TriaVertexInput(SurfaceEnum,surface));
    1712         this->inputs->AddInput(new TriaVertexInput(BedEnum,bed));
     1704        this->inputs->AddInput(new TriaVertexInput(ThicknessEnum,newthickness));
     1705        this->inputs->AddInput(new TriaVertexInput(SurfaceEnum,newsurface));
     1706        this->inputs->AddInput(new TriaVertexInput(BedEnum,newbed));
    17131707
    17141708        /*Free ressources:*/
    … …  
    19831977void  Tria::MigrateGroundingLine(void){
    19841978
    1985 
    19861979        double *values         = NULL;
    19871980        double  h[3],s[3],b[3],ba[3];
    … …  
    19921985        bool    elementonshelf = false;
    19931986
    1994 
    19951987        /*Recover info at the vertices: */
    1996         Input* surface_input =inputs->GetInput(SurfaceEnum); _assert_(surface_input);
    1997         Input* bed_input     =inputs->GetInput(BedEnum);     _assert_(bed_input);
    1998         if((surface_input->ObjectEnum()!=TriaVertexInputEnum) | (bed_input->ObjectEnum()!=TriaVertexInputEnum))_error_(" not supported yet for bed and surface interpolations not P1!");
    1999 
    20001988        GetInputListOnVertices(&h[0],ThicknessEnum);
    20011989        GetInputListOnVertices(&s[0],SurfaceEnum);
    … …  
    20202008                                b[i]=ba[i];
    20212009                                s[i]=b[i]+h[i];
    2022                                 if((nodes[i]->Sid()+1)==36)printf("MigrateGroundingLine: El %i Node %i grounding %g %g %g\n",this->Id(),nodes[i]->Sid()+1,b[i],ba[i],s[i]);
    20232010                        }
    20242011                        else{
    20252012                                /*do nothing, we are still floating.*/
    2026                                 if((nodes[i]->Sid()+1)==36)printf("MigrateGroundingLine: El %i Node %i still floating %g %g %g\n",this->Id(),nodes[i]->Sid()+1,b[i],ba[i],s[i]);
    20272013                        }
    20282014                }
    … …  
    20352021                                s[i]=(1-density)*h[i];
    20362022                                b[i]=-density*h[i];
    2037                                 printf("%i\n",nodes[i]->Sid()+1);
    2038                                 if((nodes[i]->Sid()+1)==36)printf("MigrateGroundingLine: El %i Node %i floating %g %g %g\n",this->Id(),nodes[i]->Sid()+1,b[i],ba[i],s[i]);
    20392023                        }
    20402024                        else{
    20412025                                /*do nothing, we are still grounded.*/
    2042                                 if((nodes[i]->Sid()+1)==36)printf("MigrateGroundingLine: El %i Node %i still grounded %g %g %g\n",this->Id(),nodes[i]->Sid()+1,b[i],ba[i],s[i]);
    20432026                        }
    20442027                }
    … …  
    20462029
    20472030        /*Surface and bed are updated. Update inputs:*/
    2048         surface_input->GetValuesPtr(&values,NULL); for(i=0;i<3;i++)values[i]=s[i];
    2049         bed_input->GetValuesPtr(&values,NULL);     for(i=0;i<3;i++)values[i]=b[i];
     2031        this->inputs->AddInput(new TriaVertexInput(SurfaceEnum,&s[0]));
     2032        this->inputs->AddInput(new TriaVertexInput(BedEnum,&b[0]));
    20502033       
    2051         for(i=0;i<3;i++){
    2052                 isonshelf[i]=nodes[i]->IsOnShelf();
    2053                 if((nodes[i]->Sid()+1)==36)printf("MigrateGroundingLine: El %i Node %i second shelf status %i\n",this->Id(),nodes[i]->Sid()+1,isonshelf[i]);
    2054         }
     2034        for(i=0;i<3;i++) isonshelf[i]=nodes[i]->IsOnShelf();
    20552035}
    20562036/*}}}*/
    … …  
    22072187void  Tria::ShelfSync(void){
    22082188
     2189        double  melting[NUMVERTICES];
     2190        double  h[NUMVERTICES],s[NUMVERTICES],b[NUMVERTICES],ba[NUMVERTICES];
     2191        double  bed_hydro;
     2192        double  rho_water,rho_ice,density;
     2193        int     i;
     2194        bool    elementonshelf = false;
     2195
     2196        /*melting rate at the grounding line: */
     2197        double  yts;
     2198        int     swap;
     2199        double  gl_melting_rate;
     2200
     2201        /*recover parameters: */
     2202        parameters->FindParam(&yts,ConstantsYtsEnum);
     2203        parameters->FindParam(&gl_melting_rate,GroundinglineMeltingRateEnum);
     2204
     2205        /*Recover info at the vertices: */
     2206        GetInputListOnVertices(&h[0],ThicknessEnum);
     2207        GetInputListOnVertices(&s[0],SurfaceEnum);
     2208        GetInputListOnVertices(&b[0],BedEnum);
     2209        GetInputListOnVertices(&ba[0],BathymetryEnum);
     2210
     2211        /*material parameters: */
     2212        rho_water=matpar->GetRhoWater();
     2213        rho_ice=matpar->GetRhoIce();
     2214        density=rho_ice/rho_water;
     2215
     2216        /*go through vertices, and update inputs, considering them to be TriaVertex type: */
     2217        for(i=0;i<NUMVERTICES;i++){
     2218                if(b[i]==ba[i]){
     2219                        nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,false));
     2220                        nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,true));
     2221                }
     2222                else{
     2223                        nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,true));
     2224                        nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,false));
     2225                }
     2226        }
     2227
     2228        /*Now, update  shelf status of element. An element can only be on shelf if all its nodes are on shelf: */
     2229        swap=0;
     2230        elementonshelf=false;
     2231        for(i=0;i<NUMVERTICES;i++){
     2232                if(nodes[i]->IsOnShelf()){
     2233                        elementonshelf=true;
     2234                        break;
     2235                }
     2236        }
     2237        if(!this->IsOnShelf() && elementonshelf==true)swap=1;
     2238    this->inputs->AddInput(new BoolInput(MaskElementonfloatingiceEnum,elementonshelf));
     2239       
     2240    /*If this element just  became ungrounded, set its basal melting rate at 50 m/yr:*/
     2241        if(swap){
     2242                for(i=0;i<NUMVERTICES;i++)melting[i]=gl_melting_rate/yts;
     2243                this->inputs->AddInput(new TriaVertexInput(BasalforcingsMeltingRateEnum,&melting[0]));
     2244        }
     2245}
     2246/*}}}*/
     2247/*FUNCTION Tria::SoftMigration{{{1*/
     2248void  Tria::SoftMigration(double* sheet_ungrounding){
    22092249
    22102250        double *values         = NULL;
    … …  
    22152255        bool    elementonshelf = false;
    22162256
    2217         /*melting rate at the grounding line: */
    2218         double  yts;
    2219         int     swap;
    2220         double  gl_melting_rate;
    2221 
    2222         /*recover parameters: */
    2223         parameters->FindParam(&yts,ConstantsYtsEnum);
    2224         parameters->FindParam(&gl_melting_rate,GroundinglineMeltingRateEnum);
    2225 
    22262257        /*Recover info at the vertices: */
    2227         Input* surface_input =inputs->GetInput(SurfaceEnum); _assert_(surface_input);
    2228         Input* bed_input     =inputs->GetInput(BedEnum);     _assert_(bed_input);
    2229         if((surface_input->ObjectEnum()!=TriaVertexInputEnum) | (bed_input->ObjectEnum()!=TriaVertexInputEnum))_error_(" not supported yet for bed and surface interpolations not P1!");
    2230 
    2231         GetInputListOnVertices(&h[0],ThicknessEnum);
    2232         GetInputListOnVertices(&s[0],SurfaceEnum);
    2233         GetInputListOnVertices(&b[0],BedEnum);
    2234         GetInputListOnVertices(&ba[0],BathymetryEnum);
    2235 
    2236         /*material parameters: */
    2237         rho_water=matpar->GetRhoWater();
    2238         rho_ice=matpar->GetRhoIce();
    2239         density=rho_ice/rho_water;
    2240 
    2241         /*go through vertices, and update inputs, considering them to be TriaVertex type: */
    2242         for(i=0;i<3;i++){
    2243                 if(b[i]==ba[i]){
    2244                                
    2245                         nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,false));
    2246                         nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,true));
    2247                 }
    2248                 else{
    2249                         nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,true));
    2250                         nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,false));
    2251 
    2252                 }
    2253         }
    2254 
    2255         /*Now, update  shelf status of element. An element can only be on shelf if all its nodes are on shelf: */
    2256         swap=0;
    2257         elementonshelf=false;
    2258         for(i=0;i<3;i++){
    2259                 if(nodes[i]->IsOnShelf()){
    2260                         elementonshelf=true;
    2261                         break;
    2262                 }
    2263         }
    2264         if(!this->IsOnShelf() && elementonshelf==true)swap=1;
    2265     this->inputs->AddInput(new BoolInput(MaskElementonfloatingiceEnum,elementonshelf));
    2266        
    2267     /*If this element just  became ungrounded, set its basal melting rate at 50 m/yr:*/
    2268         if(swap){
    2269                 Input* basal_melting_rate_input     =inputs->GetInput(BasalforcingsMeltingRateEnum);     _assert_(basal_melting_rate_input);
    2270                 basal_melting_rate_input->GetValuesPtr(&values,NULL); for(i=0;i<3;i++)values[i]=gl_melting_rate/yts;
    2271         }
    2272 
    2273 
    2274 }
    2275 /*}}}*/
    2276 /*FUNCTION Tria::SoftMigration{{{1*/
    2277 void  Tria::SoftMigration(double* sheet_ungrounding){
    2278 
    2279 
    2280         double *values         = NULL;
    2281         double  h[3],s[3],b[3],ba[3];
    2282         double  bed_hydro;
    2283         double  rho_water,rho_ice,density;
    2284         int     i;
    2285         bool    elementonshelf = false;
    2286 
    2287         /*Recover info at the vertices: */
    2288         Input* surface_input =inputs->GetInput(SurfaceEnum); _assert_(surface_input);
    2289         Input* bed_input     =inputs->GetInput(BedEnum);     _assert_(bed_input);
    2290         if((surface_input->ObjectEnum()!=TriaVertexInputEnum) | (bed_input->ObjectEnum()!=TriaVertexInputEnum))_error_(" not supported yet for bed and surface interpolations not P1!");
    2291 
    22922258        GetInputListOnVertices(&h[0],ThicknessEnum);
    22932259        GetInputListOnVertices(&s[0],SurfaceEnum);
    … …  
    23332299
    23342300        /*Surface and bed are updated. Update inputs:*/   
    2335         surface_input->GetValuesPtr(&values,NULL); for(i=0;i<3;i++)values[i]=s[i];
    2336         bed_input->GetValuesPtr(&values,NULL);     for(i=0;i<3;i++)values[i]=b[i];
    2337 
    2338         }
     2301        this->inputs->AddInput(new TriaVertexInput(SurfaceEnum,&s[0]));
     2302        this->inputs->AddInput(new TriaVertexInput(BedEnum,&b[0]));
     2303}
    23392304/*}}}*/
    23402305/*FUNCTION Tria::SurfaceArea {{{1*/
    … …  
    32443209        const int    numdof=NDOF2*NUMVERTICES;
    32453210
    3246         int i,dummy;
    3247         int*         doflist=NULL;
    3248         double       vx,vy;
    3249         double       values[numdof];
    3250         GaussTria*   gauss=NULL;
     3211        int        i;
     3212        double     vx,vy;
     3213        double     values[numdof];
     3214        int       *doflist = NULL;
     3215        GaussTria *gauss   = NULL;
    32513216
    32523217        /*Get dof list: */
    … …  
    32843249
    32853250        int       i;
    3286         int       dummy;
    32873251        int*      doflist=NULL;
    32883252        double    rho_ice,g;
    … …  
    33453309       
    33463310        int       i;
    3347         int       dummy;
    33483311        int*      doflist=NULL;
    33493312        double    rho_ice,g;
    … …  
    33553318        double    pressure[NUMVERTICES];
    33563319        double    thickness[NUMVERTICES];
    3357         double*   vz_ptr=NULL;
    3358         Input*    vz_input=NULL;
    33593320       
    33603321        /*Get dof list: */
    … …  
    33743335        }
    33753336
    3376         /*Get Vz*/
    3377         vz_input=inputs->GetInput(VzEnum);
    3378         if (vz_input){
    3379                 if (vz_input->ObjectEnum()!=TriaVertexInputEnum){
    3380                         _error_("Cannot compute Vel as Vz is of type %s",EnumToStringx(vz_input->ObjectEnum()));
    3381                 }
    3382                 vz_input->GetValuesPtr(&vz_ptr,&dummy);
    3383                 for(i=0;i<NUMVERTICES;i++) vz[i]=vz_ptr[i];
    3384         }
    3385         else{
    3386                 for(i=0;i<NUMVERTICES;i++) vz[i]=0.0;
    3387         }
    3388 
    33893337        /*Now Compute vel*/
     3338        GetInputListOnVertices(&vz[0],VzEnum,0.0); //default is 0
    33903339        for(i=0;i<NUMVERTICES;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
    33913340
    … …  
    34133362}
    34143363/*}}}*/
    3415 
    34163364#endif
    34173365
    … …  
    51035051        const int    numdof=NDOF1*NUMVERTICES;
    51045052
    5105         int i,dummy;
     5053        int i;
    51065054        int*         doflist=NULL;
    51075055        double       watercolumn;
    … …  
    55645512/*}}}*/
    55655513#endif
    5566 
Note: See TracChangeset for help on using the changeset viewer.