Changeset 18736


Ignore:
Timestamp:
11/05/14 15:40:11 (10 years ago)
Author:
srebuffi
Message:

CHG: added CalvingRate computation

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

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp

    r18470 r18736  
    4848
    4949        int finiteelement;
     50        bool   islevelset;
    5051
    5152        iomodel->Constant(&finiteelement,DamageElementinterpEnum);
     53        iomodel->Constant(&islevelset,TransientIslevelsetEnum);
    5254
    5355        /*Update elements: */
    … …  
    6870        iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
    6971        iomodel->FetchDataToInput(elements,PressureEnum);
     72
     73        if(islevelset){
     74                iomodel->FetchDataToInput(elements,IceMaskNodeActivationEnum);
     75        }
    7076
    7177}/*}}}*/
    … …  
    441447void DamageEvolutionAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
    442448        /*Default, do nothing*/
     449        bool islevelset;
     450        femmodel->parameters->FindParam(&islevelset,TransientIslevelsetEnum);
     451        if(islevelset){
     452                SetActiveNodesLSMx(femmodel);
     453        }
    443454        return;
    444455}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp

    r18733 r18736  
    105105        IssmDouble Jdet, dt, D_scalar;
    106106        IssmDouble h,hx,hy,hz;
    107         IssmDouble vel, calvingrate;
    108         IssmDouble norm_dlsf;
     107        IssmDouble vel;
     108//      IssmDouble norm_dlsf;
    109109        IssmDouble* xyz_list = NULL;
    110110
    … …  
    130130        IssmDouble*    w        = xNew<IssmDouble>(dim);
    131131        IssmDouble*    c        = xNew<IssmDouble>(dim);
    132         IssmDouble*    dlsf     = xNew<IssmDouble>(dim);
     132        //IssmDouble*    dlsf     = xNew<IssmDouble>(dim);
    133133
    134134        /*Retrieve all inputs and parameters*/
    … …  
    137137        Input* vx_input=NULL;
    138138        Input* vy_input=NULL;
     139        Input* calvingratex_input=NULL;
     140        Input* calvingratey_input=NULL;
    139141        if(domaintype==Domain2DhorizontalEnum){
    140142                vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input);
    141143                vy_input=basalelement->GetInput(VyEnum); _assert_(vy_input);
     144                calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
     145                calvingratey_input=basalelement->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
    142146        }
    143147        else{
    144148                if(dim==1){
    145149                        vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input);
     150                        calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
    146151                }
    147152                if(dim==2){
    148153                        vx_input=basalelement->GetInput(VxAverageEnum); _assert_(vx_input);
    149154                        vy_input=basalelement->GetInput(VyAverageEnum); _assert_(vy_input);
    150                 }
    151         }
    152 
    153         Input* lsf_slopex_input  = basalelement->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
    154         Input* lsf_slopey_input  = basalelement->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
    155         Input* calvingrate_input  = basalelement->GetInput(MasstransportCalvingrateEnum);     _assert_(calvingrate_input);
     155                        calvingratex_input=basalelement->GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input);
     156                        calvingratey_input=basalelement->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input);
     157                }
     158        }
     159
     160//      Input* lsf_slopex_input  = basalelement->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
     161//      Input* lsf_slopey_input  = basalelement->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
     162        //Input* calvingrate_input  = basalelement->GetInput(MasstransportCalvingrateEnum);     _assert_(calvingrate_input);
    156163       
    157164        /* Start  looping on the number of gaussian points: */
    … …  
    178185                vx_input->GetInputValue(&v[0],gauss); // in 3D case, add mesh velocity
    179186                vy_input->GetInputValue(&v[1],gauss);
    180                 lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
    181                 lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
    182                 calvingrate_input->GetInputValue(&calvingrate,gauss);
    183 
    184                 norm_dlsf=0.;
    185                 for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
    186                 norm_dlsf=sqrt(norm_dlsf);
    187 
    188                 if(norm_dlsf>1.e-10)
    189                         for(i=0;i<dim;i++) c[i]=calvingrate*dlsf[i]/norm_dlsf;
    190                 else
    191                         for(i=0;i<dim;i++) c[i]=0.;
     187                calvingratex_input->GetInputValue(&c[0],gauss); // in 3D case, add mesh velocity
     188                calvingratey_input->GetInputValue(&c[1],gauss);
     189                //lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
     190                //lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
     191                //calvingrate_input->GetInputValue(&calvingrate,gauss);
     192
     193                //norm_dlsf=0.;
     194                //for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
     195                //norm_dlsf=sqrt(norm_dlsf);
     196
     197                //if(norm_dlsf>1.e-10)
     198                //      for(i=0;i<dim;i++) c[i]=calvingrate*dlsf[i]/norm_dlsf;
     199                //else
     200                //for(i=0;i<dim;i++) c[i]=0.;
    192201               
    193202                for(i=0;i<dim;i++) w[i]=v[i]-c[i];
    … …  
    206215
    207216                /* Stabilization */
    208                 int stabilization=2;
     217                int stabilization=1;
    209218                vel=0.;
    210                 for(i=0;i<dim;i++) vel+=v[i]*v[i];
     219                for(i=0;i<dim;i++) vel+=w[i]*w[i];
    211220                vel=sqrt(vel)+1.e-14;
    212221                switch(stabilization){
    … …  
    217226                                /* Artificial Diffusion */
    218227                                basalelement->ElementSizes(&hx,&hy,&hz);
    219                                 h=sqrt( pow(hx*v[0]/vel,2) + pow(hy*v[1]/vel,2) );
     228                                h=sqrt( pow(hx*w[0]/vel,2) + pow(hy*w[1]/vel,2) );
    220229                                kappa=h*vel/2.;
    221230                                for(row=0;row<dim;row++)
    … …  
    234243                                /* Streamline Upwinding */
    235244                                basalelement->ElementSizes(&hx,&hy,&hz);
    236                                 h=sqrt( pow(hx*v[0]/vel,2) + pow(hy*v[1]/vel,2) );
     245                                h=sqrt( pow(hx*w[0]/vel,2) + pow(hy*w[1]/vel,2) );
    237246                                for(row=0;row<dim;row++)
    238247                                        for(col=0;col<dim;col++)
    239                                                 D[row*dim+col] = D_scalar*h/(2.*vel)*v[row]*v[col];
     248                                                D[row*dim+col] = D_scalar*h/(2.*vel)*w[row]*w[col];
    240249
    241250                                TripleMultiply(Bprime,dim,numnodes,1,
    … …  
    258267        xDelete<IssmDouble>(w);
    259268        xDelete<IssmDouble>(c);
    260         xDelete<IssmDouble>(dlsf);
     269        //xDelete<IssmDouble>(dlsf);
    261270        delete gauss;
    262271        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};

  • issm/trunk-jpl/src/c/classes/Elements/Element.cpp

    r18699 r18736  
    12101210                                input=this->inputs->GetInput(output_enum);
    12111211                                break;
     1212                        case CalvingratexEnum:
     1213                        case CalvingrateyEnum:
     1214                        case MasstransportCalvingrateEnum:
     1215                                this->StrainRateparallel();
     1216                                this->StrainRateperpendicular();
     1217                                this->CalvingRateLevermann();
     1218                                input=this->inputs->GetInput(output_enum);
     1219                                break;
    12121220                        case StrainRateparallelEnum:
    12131221                                this->StrainRateparallel();

  • issm/trunk-jpl/src/c/classes/Elements/Element.h

    r18699 r18736  
    221221                virtual void    StressIntensityFactor(void)=0;
    222222                virtual void    StrainRateparallel(void)=0;
    223                 virtual void   StrainRateperpendicular(void)=0;
     223                virtual void    StrainRateperpendicular(void)=0;
     224                virtual void    CalvingRateLevermann(void)=0;
    224225
    225226                virtual void   Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finite_element)=0;

  • issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

    r18717 r18736  
    345345        /*Clean up and return*/
    346346        delete gauss;
     347}
     348/*}}}*/
     349void       Penta::CalvingRateLevermann(){/*{{{*/
     350
     351        IssmDouble  xyz_list[NUMVERTICES][3];
     352        GaussPenta* gauss=NULL;
     353        IssmDouble  vx,vy,vel;
     354        IssmDouble  strainparallel;
     355        IssmDouble  propcoeff;
     356        IssmDouble  strainperpendicular;
     357        IssmDouble  calvingratex[NUMVERTICES];
     358        IssmDouble  calvingratey[NUMVERTICES];
     359        IssmDouble  calvingrate[NUMVERTICES];
     360
     361
     362        /* Get node coordinates and dof list: */
     363        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     364
     365        /*Retrieve all inputs and parameters we will need*/
     366        Input* vx_input=inputs->GetInput(VxEnum);                                                                                                                                               _assert_(vx_input);
     367        Input* vy_input=inputs->GetInput(VyEnum);                                                                                                                                               _assert_(vy_input);
     368        Input* strainparallel_input=inputs->GetInput(StrainRateparallelEnum);                                                           _assert_(strainparallel_input);
     369        Input* strainperpendicular_input=inputs->GetInput(StrainRateperpendicularEnum);             _assert_(strainperpendicular_input);
     370        this->parameters->FindParam(&propcoeff,MasstransportLevermannCalvingCoeffEnum);
     371
     372        /* Start looping on the number of vertices: */
     373        gauss=new GaussPenta();
     374        for (int iv=0;iv<NUMVERTICES;iv++){
     375                gauss->GaussVertex(iv);
     376
     377                /* Get the value we need*/
     378                vx_input->GetInputValue(&vx,gauss);
     379                vy_input->GetInputValue(&vy,gauss);
     380                vel=vx*vx+vy*vy;
     381                strainparallel_input->GetInputValue(&strainparallel,gauss);
     382                strainperpendicular_input->GetInputValue(&strainperpendicular,gauss);
     383
     384                /*Calving rate proportionnal to the positive product of the strain rate along the ice flow direction and the strain rate perpendicular to the ice flow */
     385                calvingrate[iv]=propcoeff*strainparallel*strainperpendicular;   
     386                if(calvingrate[iv]<0){
     387                        calvingrate[iv]=0;
     388                }
     389                calvingratex[iv]=calvingrate[iv]*vx/(vel+1.e-6);
     390                calvingratey[iv]=calvingrate[iv]*vy/(vel+1.e-6);
     391        }
     392
     393        /*Add input*/
     394        this->inputs->AddInput(new PentaInput(CalvingratexEnum,&calvingratex[0],P1Enum));
     395        this->inputs->AddInput(new PentaInput(CalvingrateyEnum,&calvingratey[0],P1Enum));
     396        this->inputs->AddInput(new PentaInput(MasstransportCalvingrateEnum,&calvingrate[0],P1Enum));
     397
     398        /*Clean up and return*/
     399        delete gauss;
     400
    347401}
    348402/*}}}*/
    … …  
    426480void       Penta::StrainRateparallel(){/*{{{*/
    427481
    428         IssmDouble  xyz_list[NUMVERTICES][3];
     482        IssmDouble *xyz_list = NULL;
     483        IssmDouble  epsilon[6];
    429484        GaussPenta* gauss=NULL;
    430485        IssmDouble  vx,vy,vel;
    … …  
    435490
    436491        /* Get node coordinates and dof list: */
    437         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    438        
     492        this->GetVerticesCoordinates(&xyz_list);
     493
    439494        /*Retrieve all inputs we will need*/
    440495        Input* vx_input=inputs->GetInput(VxEnum);                                  _assert_(vx_input);
    441496        Input* vy_input=inputs->GetInput(VyEnum);                                  _assert_(vy_input);
    442         Input* strainxx_input=inputs->GetInput(StrainRatexxEnum);             _assert_(strainxx_input);
    443         Input* strainxy_input=inputs->GetInput(StrainRatexyEnum);             _assert_(strainxy_input);
    444         Input* strainyy_input=inputs->GetInput(StrainRateyyEnum);             _assert_(strainyy_input);
     497        Input* vz_input=inputs->GetInput(VzEnum);                                                                                               _assert_(vz_input);
    445498
    446499        /* Start looping on the number of vertices: */
    … …  
    453506                vy_input->GetInputValue(&vy,gauss);
    454507                vel=vx*vx+vy*vy;
    455                 strainxx_input->GetInputValue(&strainxx,gauss);
    456                 strainxy_input->GetInputValue(&strainxy,gauss);
    457                 strainyy_input->GetInputValue(&strainyy,gauss);
     508
     509                /*Compute strain rate and viscosity: */
     510                this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input);
     511                strainxx=epsilon[0];
     512                strainyy=epsilon[1];
     513                strainxy=epsilon[3];
    458514
    459515                /*strainparallel= Strain rate along the ice flow direction */
    … …  
    466522        /*Clean up and return*/
    467523        delete gauss;
     524        xDelete<IssmDouble>(xyz_list);
    468525}
    469526/*}}}*/
    470527void       Penta::StrainRateperpendicular(){/*{{{*/
    471528
    472         IssmDouble  xyz_list[NUMVERTICES][3];
     529        IssmDouble *xyz_list = NULL;
     530        IssmDouble  epsilon[6];
    473531        GaussPenta* gauss=NULL;
    474532        IssmDouble  vx,vy,vel;
    … …  
    479537
    480538        /* Get node coordinates and dof list: */
    481         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     539        this->GetVerticesCoordinates(&xyz_list);
    482540
    483541        /*Retrieve all inputs we will need*/
    484542        Input* vx_input=inputs->GetInput(VxEnum);                                  _assert_(vx_input);
    485543        Input* vy_input=inputs->GetInput(VyEnum);                                  _assert_(vy_input);
    486         Input* strainxx_input=inputs->GetInput(StrainRatexxEnum);             _assert_(strainxx_input);
    487         Input* strainxy_input=inputs->GetInput(StrainRatexyEnum);             _assert_(strainxy_input);
    488         Input* strainyy_input=inputs->GetInput(StrainRateyyEnum);             _assert_(strainyy_input);
    489 
     544        Input* vz_input=inputs->GetInput(VzEnum);                                                                                               _assert_(vz_input);
    490545
    491546        /* Start looping on the number of vertices: */
    … …  
    498553                vy_input->GetInputValue(&vy,gauss);
    499554                vel=vx*vx+vy*vy;
    500                 strainxx_input->GetInputValue(&strainxx,gauss);
    501                 strainxy_input->GetInputValue(&strainxy,gauss);
    502                 strainyy_input->GetInputValue(&strainyy,gauss);
     555
     556                /*Compute strain rate and viscosity: */
     557                this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input);
     558                strainxx=epsilon[0];
     559                strainyy=epsilon[1];
     560                strainxy=epsilon[3];
    503561
    504562                /*strainperpendicular= Strain rate perpendicular to the ice flow direction */
    … …  
    511569        /*Clean up and return*/
    512570        delete gauss;
     571        xDelete<IssmDouble>(xyz_list);
    513572}
    514573/*}}}*/
    … …  
    21762235        if(this->inputs->GetInput(VxEnum)) this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
    21772236        if(this->inputs->GetInput(VyEnum)) this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
     2237        if(this->inputs->GetInput(CalvingratexEnum)) this->InputDepthAverageAtBase(CalvingratexEnum,CalvingratexAverageEnum);
     2238        if(this->inputs->GetInput(CalvingrateyEnum)) this->InputDepthAverageAtBase(CalvingrateyEnum,CalvingrateyAverageEnum);
    21782239        Tria* tria=(Tria*)SpawnTria(0,1,2);
    21792240        this->inputs->DeleteInput(MaterialsRheologyBbarEnum);
    … …  
    21812242        this->inputs->DeleteInput(VxAverageEnum);
    21822243        this->inputs->DeleteInput(VyAverageEnum);
     2244        this->inputs->DeleteInput(CalvingratexAverageEnum);
     2245        this->inputs->DeleteInput(CalvingrateyAverageEnum);
    21832246
    21842247        return tria;

  • issm/trunk-jpl/src/c/classes/Elements/Penta.h

    r18699 r18736  
    6060                void   StrainRateparallel();
    6161                void   StrainRateperpendicular();
     62                void   CalvingRateLevermann();
    6263                void   Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
    6364                void   ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz);

  • issm/trunk-jpl/src/c/classes/Elements/Seg.h

    r18699 r18736  
    5858                void        ComputeDeviatoricStressTensor(){_error_("not implemented yet");};
    5959                void        StressIntensityFactor(void){_error_("not implemented yet");};
    60                 void                    StrainRateparallel(void){_error_("not implemented yet");};
    61                 void                    StrainRateperpendicular(void){_error_("not implemented yet");};
     60                void        StrainRateparallel(void){_error_("not implemented yet");};
     61                void        StrainRateperpendicular(void){_error_("not implemented yet");};
     62                void        CalvingRateLevermann(void){_error_("not implemented yet");};
    6263                void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters){_error_("not implemented yet");};
    6364                void        SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters){_error_("not implemented yet");};

  • issm/trunk-jpl/src/c/classes/Elements/Tetra.h

    r18699 r18736  
    6060                void        StrainRateparallel(void){_error_("not implemented yet");};
    6161                void        StrainRateperpendicular(void){_error_("not implemented yet");};
     62                void        CalvingRateLevermann(void){_error_("not implemented yet");};
    6263                void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters);
    6364                void        SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters);

  • issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

    r18717 r18736  
    276276}
    277277/*}}}*/
    278 void       Tria::StrainRateperpendicular(){/*{{{*/
    279 
    280         IssmDouble  xyz_list[NUMVERTICES][3];
    281         GaussPenta* gauss=NULL;
    282         IssmDouble  vx,vy,vel;
    283         IssmDouble  strainxx;
    284         IssmDouble  strainxy;
    285         IssmDouble  strainyy;
    286         IssmDouble  strainperpendicular[NUMVERTICES];
    287 
    288         /* Get node coordinates and dof list: */
    289         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    290 
    291         /*Retrieve all inputs we will need*/
    292         Input* vx_input=inputs->GetInput(VxEnum);                                  _assert_(vx_input);
    293         Input* vy_input=inputs->GetInput(VyEnum);                                  _assert_(vy_input);
    294         Input* strainxx_input=inputs->GetInput(StrainRatexxEnum);             _assert_(strainxx_input);
    295         Input* strainxy_input=inputs->GetInput(StrainRatexyEnum);             _assert_(strainxy_input);
    296         Input* strainyy_input=inputs->GetInput(StrainRateyyEnum);             _assert_(strainyy_input);
    297 
    298 
    299         /* Start looping on the number of vertices: */
    300         gauss=new GaussPenta();
    301         for (int iv=0;iv<NUMVERTICES;iv++){
    302                 gauss->GaussVertex(iv);
    303 
    304                 /* Get the value we need*/
    305                 vx_input->GetInputValue(&vx,gauss);
    306                 vy_input->GetInputValue(&vy,gauss);
    307                 vel=vx*vx+vy*vy;
    308                 strainxx_input->GetInputValue(&strainxx,gauss);
    309                 strainxy_input->GetInputValue(&strainxy,gauss);
    310                 strainyy_input->GetInputValue(&strainyy,gauss);
    311 
    312                 /*strainperpendicular= Strain rate perpendicular to the ice flow direction */
    313                 strainperpendicular[iv]=(vx*vx*(strainyy)+vy*vy*(strainxx)-2*vy*vx*strainxy)/(vel+1.e-6);
    314         }
    315 
    316         /*Add input*/
    317         this->inputs->AddInput(new PentaInput(StrainRateperpendicularEnum,&strainperpendicular[0],P1Enum));
    318 
    319         /*Clean up and return*/
    320         delete gauss;
    321 }
    322 /*}}}*/
    323 void       Tria::StrainRateparallel(){/*{{{*/
    324 
    325         IssmDouble  xyz_list[NUMVERTICES][3];
    326         GaussPenta* gauss=NULL;
    327         IssmDouble  vx,vy,vel;
    328         IssmDouble  strainxx;
    329         IssmDouble  strainxy;
    330         IssmDouble  strainyy;
    331         IssmDouble  strainparallel[NUMVERTICES];
    332 
    333         /* Get node coordinates and dof list: */
    334         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    335        
    336         /*Retrieve all inputs we will need*/
    337         Input* vx_input=inputs->GetInput(VxEnum);                                  _assert_(vx_input);
    338         Input* vy_input=inputs->GetInput(VyEnum);                                  _assert_(vy_input);
    339         Input* strainxx_input=inputs->GetInput(StrainRatexxEnum);             _assert_(strainxx_input);
    340         Input* strainxy_input=inputs->GetInput(StrainRatexyEnum);             _assert_(strainxy_input);
    341         Input* strainyy_input=inputs->GetInput(StrainRateyyEnum);             _assert_(strainyy_input);
    342 
    343         /* Start looping on the number of vertices: */
    344         gauss=new GaussPenta();
    345         for (int iv=0;iv<NUMVERTICES;iv++){
    346                 gauss->GaussVertex(iv);
    347 
    348                 /* Get the value we need*/
    349                 vx_input->GetInputValue(&vx,gauss);
    350                 vy_input->GetInputValue(&vy,gauss);
    351                 vel=vx*vx+vy*vy;
    352                 strainxx_input->GetInputValue(&strainxx,gauss);
    353                 strainxy_input->GetInputValue(&strainxy,gauss);
    354                 strainyy_input->GetInputValue(&strainyy,gauss);
    355 
    356                 /*strainparallel= Strain rate along the ice flow direction */
    357                 strainparallel[iv]=(vx*vx*(strainxx)+vy*vy*(strainyy)+2*vy*vx*strainxy)/(vel+1.e-6);
    358         }
    359 
    360         /*Add input*/
    361         this->inputs->AddInput(new PentaInput(StrainRateparallelEnum,&strainparallel[0],P1Enum));
    362 
    363         /*Clean up and return*/
    364         delete gauss;
    365 }
    366 /*}}}*/
    367278void       Tria::ComputeDeviatoricStressTensor(){/*{{{*/
    368279
    … …  
    413324        /*Clean up and return*/
    414325        delete gauss;
     326}
     327/*}}}*/
     328void       Tria::StrainRateparallel(){/*{{{*/
     329
     330        IssmDouble *xyz_list = NULL;
     331        IssmDouble  epsilon[3];
     332        GaussTria* gauss=NULL;
     333        IssmDouble  vx,vy,vel;
     334        IssmDouble  strainxx;
     335        IssmDouble  strainxy;
     336        IssmDouble  strainyy;
     337        IssmDouble  strainparallel[NUMVERTICES];
     338
     339        /* Get node coordinates and dof list: */
     340        this->GetVerticesCoordinates(&xyz_list);
     341
     342        /*Retrieve all inputs we will need*/
     343        Input* vx_input=inputs->GetInput(VxEnum);                                  _assert_(vx_input);
     344        Input* vy_input=inputs->GetInput(VyEnum);                                  _assert_(vy_input);
     345
     346        /* Start looping on the number of vertices: */
     347        gauss=new GaussTria();
     348        for (int iv=0;iv<NUMVERTICES;iv++){
     349                gauss->GaussVertex(iv);
     350
     351                /* Get the value we need*/
     352                vx_input->GetInputValue(&vx,gauss);
     353                vy_input->GetInputValue(&vy,gauss);
     354                vel=vx*vx+vy*vy;
     355
     356                /*Compute strain rate viscosity and pressure: */
     357                this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
     358                strainxx=epsilon[0];
     359                strainyy=epsilon[1];
     360                strainxy=epsilon[2];
     361
     362                /*strainparallel= Strain rate along the ice flow direction */
     363                strainparallel[iv]=(vx*vx*(strainxx)+vy*vy*(strainyy)+2*vy*vx*strainxy)/(vel+1.e-6);
     364        }
     365
     366        /*Add input*/
     367        this->inputs->AddInput(new TriaInput(StrainRateparallelEnum,&strainparallel[0],P1Enum));
     368
     369        /*Clean up and return*/
     370        delete gauss;
     371        xDelete<IssmDouble>(xyz_list);
     372}
     373/*}}}*/
     374void       Tria::StrainRateperpendicular(){/*{{{*/
     375
     376        IssmDouble *xyz_list = NULL;
     377        GaussTria* gauss=NULL;
     378        IssmDouble  epsilon[3];
     379        IssmDouble  vx,vy,vel;
     380        IssmDouble  strainxx;
     381        IssmDouble  strainxy;
     382        IssmDouble  strainyy;
     383        IssmDouble  strainperpendicular[NUMVERTICES];
     384
     385        /* Get node coordinates and dof list: */
     386        this->GetVerticesCoordinates(&xyz_list);
     387
     388        /*Retrieve all inputs we will need*/
     389        Input* vx_input=inputs->GetInput(VxEnum);                                  _assert_(vx_input);
     390        Input* vy_input=inputs->GetInput(VyEnum);                                  _assert_(vy_input);
     391
     392        /* Start looping on the number of vertices: */
     393        gauss=new GaussTria();
     394        for (int iv=0;iv<NUMVERTICES;iv++){
     395                gauss->GaussVertex(iv);
     396
     397                /* Get the value we need*/
     398                vx_input->GetInputValue(&vx,gauss);
     399                vy_input->GetInputValue(&vy,gauss);
     400                vel=vx*vx+vy*vy;
     401
     402                /*Compute strain rate viscosity and pressure: */
     403                this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
     404                strainxx=epsilon[0];
     405                strainyy=epsilon[1];
     406                strainxy=epsilon[2];
     407
     408                /*strainperpendicular= Strain rate perpendicular to the ice flow direction */
     409                strainperpendicular[iv]=(vx*vx*(strainyy)+vy*vy*(strainxx)-2*vy*vx*strainxy)/(vel+1.e-6);
     410        }
     411
     412        /*Add input*/
     413        this->inputs->AddInput(new TriaInput(StrainRateperpendicularEnum,&strainperpendicular[0],P1Enum));
     414
     415        /*Clean up and return*/
     416        delete gauss;
     417        xDelete<IssmDouble>(xyz_list);
     418}
     419/*}}}*/
     420void       Tria::CalvingRateLevermann(){/*{{{*/
     421
     422        IssmDouble  xyz_list[NUMVERTICES][3];
     423        GaussTria* gauss=NULL;
     424        IssmDouble  vx,vy,vel;
     425        IssmDouble  strainparallel;
     426        IssmDouble  propcoeff;
     427        IssmDouble  strainperpendicular;
     428        IssmDouble  calvingratex[NUMVERTICES];
     429        IssmDouble  calvingratey[NUMVERTICES];
     430        IssmDouble  calvingrate[NUMVERTICES];
     431
     432
     433        /* Get node coordinates and dof list: */
     434        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     435
     436        /*Retrieve all inputs and parameters we will need*/
     437        Input* vx_input=inputs->GetInput(VxEnum);                                                                                                                                               _assert_(vx_input);
     438        Input* vy_input=inputs->GetInput(VyEnum);                                                                                                                                               _assert_(vy_input);
     439        Input* strainparallel_input=inputs->GetInput(StrainRateparallelEnum);                                                           _assert_(strainparallel_input);
     440        Input* strainperpendicular_input=inputs->GetInput(StrainRateperpendicularEnum);                                 _assert_(strainperpendicular_input);
     441        this->parameters->FindParam(&propcoeff,MasstransportLevermannCalvingCoeffEnum);
     442
     443        /* Start looping on the number of vertices: */
     444        gauss=new GaussTria();
     445        for (int iv=0;iv<NUMVERTICES;iv++){
     446                gauss->GaussVertex(iv);
     447
     448                /* Get the value we need*/
     449                vx_input->GetInputValue(&vx,gauss);
     450                vy_input->GetInputValue(&vy,gauss);
     451                vel=vx*vx+vy*vy;
     452                strainparallel_input->GetInputValue(&strainparallel,gauss);
     453                strainperpendicular_input->GetInputValue(&strainperpendicular,gauss);
     454
     455                /*Calving rate proportionnal to the positive product of the strain rate along the ice flow direction and the strain rate perpendicular to the ice flow */
     456                calvingrate[iv]=propcoeff*strainparallel*strainperpendicular;
     457                if(calvingrate[iv]<0){
     458                        calvingrate[iv]=0;
     459                }
     460                calvingratex[iv]=calvingrate[iv]*vx/(vel+1.e-6);
     461                calvingratey[iv]=calvingrate[iv]*vy/(vel+1.e-6);
     462        }
     463
     464        /*Add input*/
     465        this->inputs->AddInput(new TriaInput(CalvingratexEnum,&calvingratex[0],P1Enum));
     466        this->inputs->AddInput(new TriaInput(CalvingrateyEnum,&calvingratey[0],P1Enum));
     467        this->inputs->AddInput(new TriaInput(MasstransportCalvingrateEnum,&calvingrate[0],P1Enum));
     468
     469        /*Clean up and return*/
     470        delete gauss;
     471
    415472}
    416473/*}}}*/

  • issm/trunk-jpl/src/c/classes/Elements/Tria.h

    r18699 r18736  
    5555                void        ComputeStressTensor();
    5656                void        ComputeDeviatoricStressTensor();
     57                void                    StrainRateparallel();
     58                void                    StrainRateperpendicular();
    5759                void        ComputeSurfaceNormalVelocity();
    5860                void        StressIntensityFactor(void){_error_("not implemented yet");};
    59                 void                    StrainRateparallel();
    60                 void                    StrainRateperpendicular();
     61                void                    CalvingRateLevermann();
    6162                void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters);
    6263                void        SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters);

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

    r18730 r18736  
    16951695}
    16961696        /*}}}*/
     1697void FemModel::CalvingRateLevermannx(){/*{{{*/
     1698
     1699           for(int i=0;i<elements->Size();i++){
     1700                              Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
     1701                                              element->CalvingRateLevermann();
     1702                                                           }
     1703}
     1704/*}}}*/
     1705void FemModel::StrainRateparallelx(){/*{{{*/
     1706
     1707           for(int i=0;i<elements->Size();i++){
     1708                              Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
     1709                                              element->StrainRateparallel();
     1710                                                           }
     1711}
     1712/*}}}*/
     1713void FemModel::StrainRateperpendicularx(){/*{{{*/
     1714
     1715           for(int i=0;i<elements->Size();i++){
     1716                              Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
     1717                                              element->StrainRateperpendicular();
     1718                                                           }
     1719}
     1720/*}}}*/
    16971721#ifdef  _HAVE_DAKOTA_
    16981722void FemModel::DakotaResponsesx(double* d_responses,char** responses_descriptors,int numresponsedescriptors,int d_numresponses){/*{{{*/

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

    r18613 r18736  
    8282                void BalancethicknessMisfitx(IssmDouble* pV);
    8383                void StressIntensityFactorx();
     84                void StrainRateparallelx();
     85                void StrainRateperpendicularx();
     86                void CalvingRateLevermannx();
    8487                #ifdef  _HAVE_DAKOTA_
    8588                void DakotaResponsesx(double* d_responses,char** responses_descriptors,int numresponsedescriptors,int d_numresponses);

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

    r18237 r18736  
    2121        /*parameters: */
    2222        IssmDouble starttime,finaltime,dt,yts;
    23         bool       isstressbalance,ismasstransport,isFS,isthermal,isgroundingline,isgia,islevelset,isdamageevolution,ishydrology;
     23        bool       isstressbalance,ismasstransport,isFS,isthermal,isgroundingline,isgia,islevelset,isdamageevolution,ishydrology,iscalvingrate;
    2424        bool       save_results,dakota_analysis;
    2525        bool       time_adapt=false;
    … …  
    5353        femmodel->parameters->FindParam(&ishydrology,TransientIshydrologyEnum);
    5454        femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum);
     55        femmodel->parameters->FindParam(&iscalvingrate,MasstransportIscalvingrateEnum);
    5556        if(isgroundingline) femmodel->parameters->FindParam(&groundingline_migration,GroundinglineMigrationEnum);
    5657        femmodel->parameters->FindParam(&numoutputs,TransientNumRequestedOutputsEnum);
    … …  
    101102
    102103                if(islevelset){
     104                        if(iscalvingrate){
     105                                if(VerboseSolution()) _printf0_("   computing calving rate\n");
     106                                femmodel->StrainRateparallelx();
     107                                femmodel->StrainRateperpendicularx();
     108                                femmodel->CalvingRateLevermannx();
     109                        }
    103110                        if(VerboseSolution()) _printf0_("   computing movement of ice boundaries\n");
    104111                        /* smoothen slope of lsf for computation of normal on ice domain*/

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp

    r18717 r18736  
    8383                parameters->AddObject(iomodel->CopyConstantObject(TransientIsdamageevolutionEnum));
    8484                parameters->AddObject(iomodel->CopyConstantObject(TransientIshydrologyEnum));
     85                parameters->AddObject(iomodel->CopyConstantObject(MasstransportIscalvingrateEnum));
     86                parameters->AddObject(iomodel->CopyConstantObject(MasstransportLevermannCalvingCoeffEnum));
    8587                parameters->AddObject(iomodel->CopyConstantObject(MaterialsRheologyLawEnum));
    8688                parameters->AddObject(iomodel->CopyConstantObject(GiaCrossSectionShapeEnum));

  • issm/trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.cpp

    r18521 r18736  
    5959                                        analysis_type==BalancethicknessAnalysisEnum ||
    6060                                        analysis_type==HydrologyDCInefficientAnalysisEnum ||
    61                                         analysis_type==DamageEvolutionAnalysisEnum ||
     61                                        //analysis_type==DamageEvolutionAnalysisEnum ||
    6262                                        analysis_type==HydrologyDCEfficientAnalysisEnum ||
    6363                                        analysis_type==LevelsetAnalysisEnum ||

  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r18732 r18736  
    208208        NewDamageEnum,
    209209        StressIntensityFactorEnum,
     210        CalvingratexEnum,
     211        CalvingrateyEnum,
     212        CalvingratexAverageEnum,
     213        CalvingrateyAverageEnum,
    210214        StrainRateparallelEnum,
    211215        StrainRateperpendicularEnum,
    … …  
    245249        MiscellaneousNameEnum, //FIXME: only used by qmu, should not be marshalled (already in queueing script)
    246250        MasstransportHydrostaticAdjustmentEnum,
     251        MasstransportIscalvingrateEnum,
     252        MasstransportLevermannCalvingCoeffEnum,
    247253        MasstransportIsfreesurfaceEnum,
    248254        MasstransportMinThicknessEnum,

  • issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r18732 r18736  
    216216                case NewDamageEnum : return "NewDamage";
    217217                case StressIntensityFactorEnum : return "StressIntensityFactor";
     218                case CalvingratexEnum : return "Calvingratex";
     219                case CalvingrateyEnum : return "Calvingratey";
     220                case CalvingratexAverageEnum : return "CalvingratexAverage";
     221                case CalvingrateyAverageEnum : return "CalvingrateyAverage";
    218222                case StrainRateparallelEnum : return "StrainRateparallel";
    219223                case StrainRateperpendicularEnum : return "StrainRateperpendicular";
    … …  
    253257                case MiscellaneousNameEnum : return "MiscellaneousName";
    254258                case MasstransportHydrostaticAdjustmentEnum : return "MasstransportHydrostaticAdjustment";
     259                case MasstransportIscalvingrateEnum : return "MasstransportIscalvingrate";
     260                case MasstransportLevermannCalvingCoeffEnum : return "MasstransportLevermannCalvingCoeff";
    255261                case MasstransportIsfreesurfaceEnum : return "MasstransportIsfreesurface";
    256262                case MasstransportMinThicknessEnum : return "MasstransportMinThickness";

  • issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

    r18732 r18736  
    219219              else if (strcmp(name,"NewDamage")==0) return NewDamageEnum;
    220220              else if (strcmp(name,"StressIntensityFactor")==0) return StressIntensityFactorEnum;
     221              else if (strcmp(name,"Calvingratex")==0) return CalvingratexEnum;
     222              else if (strcmp(name,"Calvingratey")==0) return CalvingrateyEnum;
     223              else if (strcmp(name,"CalvingratexAverage")==0) return CalvingratexAverageEnum;
     224              else if (strcmp(name,"CalvingrateyAverage")==0) return CalvingrateyAverageEnum;
    221225              else if (strcmp(name,"StrainRateparallel")==0) return StrainRateparallelEnum;
    222226              else if (strcmp(name,"StrainRateperpendicular")==0) return StrainRateperpendicularEnum;
    … …  
    256260              else if (strcmp(name,"MiscellaneousName")==0) return MiscellaneousNameEnum;
    257261              else if (strcmp(name,"MasstransportHydrostaticAdjustment")==0) return MasstransportHydrostaticAdjustmentEnum;
     262         else stage=3;
     263   }
     264   if(stage==3){
     265              if (strcmp(name,"MasstransportIscalvingrate")==0) return MasstransportIscalvingrateEnum;
     266              else if (strcmp(name,"MasstransportLevermannCalvingCoeff")==0) return MasstransportLevermannCalvingCoeffEnum;
    258267              else if (strcmp(name,"MasstransportIsfreesurface")==0) return MasstransportIsfreesurfaceEnum;
    259268              else if (strcmp(name,"MasstransportMinThickness")==0) return MasstransportMinThicknessEnum;
    260269              else if (strcmp(name,"MasstransportPenaltyFactor")==0) return MasstransportPenaltyFactorEnum;
    261270              else if (strcmp(name,"MasstransportSpcthickness")==0) return MasstransportSpcthicknessEnum;
    262          else stage=3;
    263    }
    264    if(stage==3){
    265               if (strcmp(name,"MasstransportCalvingrate")==0) return MasstransportCalvingrateEnum;
     271              else if (strcmp(name,"MasstransportCalvingrate")==0) return MasstransportCalvingrateEnum;
    266272              else if (strcmp(name,"MasstransportStabilization")==0) return MasstransportStabilizationEnum;
    267273              else if (strcmp(name,"MasstransportVertexPairing")==0) return MasstransportVertexPairingEnum;
    … …  
    377383              else if (strcmp(name,"BalancethicknessSolution")==0) return BalancethicknessSolutionEnum;
    378384              else if (strcmp(name,"Balancethickness2Analysis")==0) return Balancethickness2AnalysisEnum;
    379               else if (strcmp(name,"Balancethickness2Solution")==0) return Balancethickness2SolutionEnum;
     385         else stage=4;
     386   }
     387   if(stage==4){
     388              if (strcmp(name,"Balancethickness2Solution")==0) return Balancethickness2SolutionEnum;
    380389              else if (strcmp(name,"BalancethicknessSoftAnalysis")==0) return BalancethicknessSoftAnalysisEnum;
    381390              else if (strcmp(name,"BalancethicknessSoftSolution")==0) return BalancethicknessSoftSolutionEnum;
    … …  
    383392              else if (strcmp(name,"BalancevelocitySolution")==0) return BalancevelocitySolutionEnum;
    384393              else if (strcmp(name,"L2ProjectionEPLAnalysis")==0) return L2ProjectionEPLAnalysisEnum;
    385          else stage=4;
    386    }
    387    if(stage==4){
    388               if (strcmp(name,"L2ProjectionBaseAnalysis")==0) return L2ProjectionBaseAnalysisEnum;
     394              else if (strcmp(name,"L2ProjectionBaseAnalysis")==0) return L2ProjectionBaseAnalysisEnum;
    389395              else if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum;
    390396              else if (strcmp(name,"DamageEvolutionSolution")==0) return DamageEvolutionSolutionEnum;
    … …  
    500506              else if (strcmp(name,"TriaInput")==0) return TriaInputEnum;
    501507              else if (strcmp(name,"Tetra")==0) return TetraEnum;
    502               else if (strcmp(name,"TetraInput")==0) return TetraInputEnum;
     508         else stage=5;
     509   }
     510   if(stage==5){
     511              if (strcmp(name,"TetraInput")==0) return TetraInputEnum;
    503512              else if (strcmp(name,"Penta")==0) return PentaEnum;
    504513              else if (strcmp(name,"PentaInput")==0) return PentaInputEnum;
    … …  
    506515              else if (strcmp(name,"VertexPId")==0) return VertexPIdEnum;
    507516              else if (strcmp(name,"VertexSId")==0) return VertexSIdEnum;
    508          else stage=5;
    509    }
    510    if(stage==5){
    511               if (strcmp(name,"Air")==0) return AirEnum;
     517              else if (strcmp(name,"Air")==0) return AirEnum;
    512518              else if (strcmp(name,"Ice")==0) return IceEnum;
    513519              else if (strcmp(name,"Melange")==0) return MelangeEnum;
    … …  
    623629              else if (strcmp(name,"P2xP4")==0) return P2xP4Enum;
    624630              else if (strcmp(name,"P1P1")==0) return P1P1Enum;
    625               else if (strcmp(name,"P1P1GLS")==0) return P1P1GLSEnum;
     631         else stage=6;
     632   }
     633   if(stage==6){
     634              if (strcmp(name,"P1P1GLS")==0) return P1P1GLSEnum;
    626635              else if (strcmp(name,"MINI")==0) return MINIEnum;
    627636              else if (strcmp(name,"MINIcondensed")==0) return MINIcondensedEnum;
    … …  
    629638              else if (strcmp(name,"LATaylorHood")==0) return LATaylorHoodEnum;
    630639              else if (strcmp(name,"XTaylorHood")==0) return XTaylorHoodEnum;
    631          else stage=6;
    632    }
    633    if(stage==6){
    634               if (strcmp(name,"OneLayerP4z")==0) return OneLayerP4zEnum;
     640              else if (strcmp(name,"OneLayerP4z")==0) return OneLayerP4zEnum;
    635641              else if (strcmp(name,"CrouzeixRaviart")==0) return CrouzeixRaviartEnum;
    636642              else if (strcmp(name,"LACrouzeixRaviart")==0) return LACrouzeixRaviartEnum;
    … …  
    746752              else if (strcmp(name,"LevelsetfunctionSlopeY")==0) return LevelsetfunctionSlopeYEnum;
    747753              else if (strcmp(name,"LevelsetfunctionPicard")==0) return LevelsetfunctionPicardEnum;
    748               else if (strcmp(name,"Seaiceatm")==0) return SeaiceatmEnum;
     754         else stage=7;
     755   }
     756   if(stage==7){
     757              if (strcmp(name,"Seaiceatm")==0) return SeaiceatmEnum;
    749758              else if (strcmp(name,"Seaiceocean")==0) return SeaiceoceanEnum;
    750759              else if (strcmp(name,"SeaiceThickness")==0) return SeaiceThicknessEnum;
    … …  
    752761              else if (strcmp(name,"SeaiceMinConcentration")==0) return SeaiceMinConcentrationEnum;
    753762              else if (strcmp(name,"SeaiceMinThickness")==0) return SeaiceMinThicknessEnum;
    754          else stage=7;
    755    }
    756    if(stage==7){
    757               if (strcmp(name,"SeaiceMaxThickness")==0) return SeaiceMaxThicknessEnum;
     763              else if (strcmp(name,"SeaiceMaxThickness")==0) return SeaiceMaxThicknessEnum;
    758764              else if (strcmp(name,"SeaiceSpcvx")==0) return SeaiceSpcvxEnum;
    759765              else if (strcmp(name,"SeaiceSpcvy")==0) return SeaiceSpcvyEnum;
Note: See TracChangeset for help on using the changeset viewer.