Changeset 18921


Ignore:
Timestamp:
12/03/14 20:48:06 (10 years ago)
Author:
seroussi
Message:

CHG: minor reordering Penta.cpp

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r18884 r18921  
    125125
    126126/*Other*/
    127 void       Penta::AddInput(int input_enum,IssmDouble* values, int interpolation_enum){/*{{{*/
    128 
    129         _assert_(this->inputs);
    130         this->inputs->AddInput(new PentaInput(input_enum,values,interpolation_enum));
    131 }
    132 /*}}}*/
    133127void       Penta::AddBasalInput(int input_enum,IssmDouble* values, int interpolation_enum){/*{{{*/
    134128
    … …  
    154148                else _error_("not implemented yet");
    155149        }
     150}
     151/*}}}*/
     152void       Penta::AddInput(int input_enum,IssmDouble* values, int interpolation_enum){/*{{{*/
     153
     154        _assert_(this->inputs);
     155        this->inputs->AddInput(new PentaInput(input_enum,values,interpolation_enum));
     156}
     157/*}}}*/
     158void       Penta::AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part){/*{{{*/
     159        _error_("Not supported yet!");
     160}
     161/*}}}*/
     162void       Penta::CalvingRateLevermann(){/*{{{*/
     163
     164        IssmDouble  xyz_list[NUMVERTICES][3];
     165        GaussPenta* gauss=NULL;
     166        IssmDouble  vx,vy,vel;
     167        IssmDouble  strainparallel;
     168        IssmDouble  propcoeff;
     169        IssmDouble  strainperpendicular;
     170        IssmDouble  calvingratex[NUMVERTICES];
     171        IssmDouble  calvingratey[NUMVERTICES];
     172        IssmDouble  calvingrate[NUMVERTICES];
     173
     174
     175        /* Get node coordinates and dof list: */
     176        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     177
     178        /*Retrieve all inputs and parameters we will need*/
     179        Input* vx_input=inputs->GetInput(VxEnum);                                                                                                                                               _assert_(vx_input);
     180        Input* vy_input=inputs->GetInput(VyEnum);                                                                                                                                               _assert_(vy_input);
     181        Input* strainparallel_input=inputs->GetInput(StrainRateparallelEnum);                                                           _assert_(strainparallel_input);
     182        Input* strainperpendicular_input=inputs->GetInput(StrainRateperpendicularEnum);              _assert_(strainperpendicular_input);
     183        Input* levermanncoeff_input=inputs->GetInput(CalvinglevermannCoeffEnum);                     _assert_(levermanncoeff_input);
     184
     185        /* Start looping on the number of vertices: */
     186        gauss=new GaussPenta();
     187        for (int iv=0;iv<NUMVERTICES;iv++){
     188                gauss->GaussVertex(iv);
     189
     190                /* Get the value we need*/
     191                vx_input->GetInputValue(&vx,gauss);
     192                vy_input->GetInputValue(&vy,gauss);
     193                vel=vx*vx+vy*vy;
     194                strainparallel_input->GetInputValue(&strainparallel,gauss);
     195                strainperpendicular_input->GetInputValue(&strainperpendicular,gauss);
     196                levermanncoeff_input->GetInputValue(&propcoeff,gauss);
     197
     198                /*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 */
     199                calvingrate[iv]=propcoeff*strainparallel*strainperpendicular;   
     200                if(calvingrate[iv]<0){
     201                        calvingrate[iv]=0;
     202                }
     203                calvingratex[iv]=calvingrate[iv]*vx/(vel+1.e-6);
     204                calvingratey[iv]=calvingrate[iv]*vy/(vel+1.e-6);
     205        }
     206
     207        /*Add input*/
     208        this->inputs->AddInput(new PentaInput(CalvingratexEnum,&calvingratex[0],P1Enum));
     209        this->inputs->AddInput(new PentaInput(CalvingrateyEnum,&calvingratey[0],P1Enum));
     210        this->inputs->AddInput(new PentaInput(CalvingCalvingrateEnum,&calvingrate[0],P1Enum));
     211
     212        /*Clean up and return*/
     213        delete gauss;
     214
    156215}
    157216/*}}}*/
    … …  
    243302}
    244303/*}}}*/
     304void       Penta::ComputeDeviatoricStressTensor(){/*{{{*/
     305
     306        IssmDouble  xyz_list[NUMVERTICES][3];
     307        IssmDouble  viscosity;
     308        IssmDouble  epsilon[6]; /* epsilon=[exx,eyy,exy];*/
     309        IssmDouble  tau_xx[NUMVERTICES];
     310        IssmDouble      tau_yy[NUMVERTICES];
     311        IssmDouble      tau_zz[NUMVERTICES];
     312        IssmDouble  tau_xy[NUMVERTICES];
     313        IssmDouble      tau_xz[NUMVERTICES];
     314        IssmDouble      tau_yz[NUMVERTICES];
     315        GaussPenta* gauss=NULL;
     316
     317        /* Get node coordinates and dof list: */
     318        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     319
     320        /*Retrieve all inputs we will be needing: */
     321        Input* vx_input=inputs->GetInput(VxEnum);             _assert_(vx_input);
     322        Input* vy_input=inputs->GetInput(VyEnum);             _assert_(vy_input);
     323        Input* vz_input=inputs->GetInput(VzEnum);             _assert_(vz_input);
     324
     325        /* Start looping on the number of vertices: */
     326        gauss=new GaussPenta();
     327        for (int iv=0;iv<NUMVERTICES;iv++){
     328                gauss->GaussVertex(iv);
     329
     330                /*Compute strain rate viscosity and pressure: */
     331                this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     332                this->ViscosityFS(&viscosity,3,&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     333
     334                /*Compute Stress*/
     335                tau_xx[iv]=2*viscosity*epsilon[0]; // tau = nu eps
     336                tau_yy[iv]=2*viscosity*epsilon[1];
     337                tau_zz[iv]=2*viscosity*epsilon[2];
     338                tau_xy[iv]=2*viscosity*epsilon[3];
     339                tau_xz[iv]=2*viscosity*epsilon[4];
     340                tau_yz[iv]=2*viscosity*epsilon[5];
     341        }
     342
     343        /*Add Stress tensor components into inputs*/
     344        this->inputs->AddInput(new PentaInput(DeviatoricStressxxEnum,&tau_xx[0],P1Enum));
     345        this->inputs->AddInput(new PentaInput(DeviatoricStressxyEnum,&tau_xy[0],P1Enum));
     346        this->inputs->AddInput(new PentaInput(DeviatoricStressxzEnum,&tau_xz[0],P1Enum));
     347        this->inputs->AddInput(new PentaInput(DeviatoricStressyyEnum,&tau_yy[0],P1Enum));
     348        this->inputs->AddInput(new PentaInput(DeviatoricStressyzEnum,&tau_yz[0],P1Enum));
     349        this->inputs->AddInput(new PentaInput(DeviatoricStresszzEnum,&tau_zz[0],P1Enum));
     350
     351        /*Clean up and return*/
     352        delete gauss;
     353}
     354/*}}}*/
    245355void       Penta::ComputeStressTensor(){/*{{{*/
    246356
    … …  
    296406}
    297407/*}}}*/
    298 void       Penta::ComputeDeviatoricStressTensor(){/*{{{*/
    299 
    300         IssmDouble  xyz_list[NUMVERTICES][3];
    301         IssmDouble  viscosity;
    302         IssmDouble  epsilon[6]; /* epsilon=[exx,eyy,exy];*/
    303         IssmDouble  tau_xx[NUMVERTICES];
    304         IssmDouble      tau_yy[NUMVERTICES];
    305         IssmDouble      tau_zz[NUMVERTICES];
    306         IssmDouble  tau_xy[NUMVERTICES];
    307         IssmDouble      tau_xz[NUMVERTICES];
    308         IssmDouble      tau_yz[NUMVERTICES];
    309         GaussPenta* gauss=NULL;
    310 
    311         /* Get node coordinates and dof list: */
    312         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    313 
    314         /*Retrieve all inputs we will be needing: */
    315         Input* vx_input=inputs->GetInput(VxEnum);             _assert_(vx_input);
    316         Input* vy_input=inputs->GetInput(VyEnum);             _assert_(vy_input);
    317         Input* vz_input=inputs->GetInput(VzEnum);             _assert_(vz_input);
    318 
    319         /* Start looping on the number of vertices: */
    320         gauss=new GaussPenta();
    321         for (int iv=0;iv<NUMVERTICES;iv++){
    322                 gauss->GaussVertex(iv);
    323 
    324                 /*Compute strain rate viscosity and pressure: */
    325                 this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    326                 this->ViscosityFS(&viscosity,3,&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    327 
    328                 /*Compute Stress*/
    329                 tau_xx[iv]=2*viscosity*epsilon[0]; // tau = nu eps
    330                 tau_yy[iv]=2*viscosity*epsilon[1];
    331                 tau_zz[iv]=2*viscosity*epsilon[2];
    332                 tau_xy[iv]=2*viscosity*epsilon[3];
    333                 tau_xz[iv]=2*viscosity*epsilon[4];
    334                 tau_yz[iv]=2*viscosity*epsilon[5];
    335         }
    336 
    337         /*Add Stress tensor components into inputs*/
    338         this->inputs->AddInput(new PentaInput(DeviatoricStressxxEnum,&tau_xx[0],P1Enum));
    339         this->inputs->AddInput(new PentaInput(DeviatoricStressxyEnum,&tau_xy[0],P1Enum));
    340         this->inputs->AddInput(new PentaInput(DeviatoricStressxzEnum,&tau_xz[0],P1Enum));
    341         this->inputs->AddInput(new PentaInput(DeviatoricStressyyEnum,&tau_yy[0],P1Enum));
    342         this->inputs->AddInput(new PentaInput(DeviatoricStressyzEnum,&tau_yz[0],P1Enum));
    343         this->inputs->AddInput(new PentaInput(DeviatoricStresszzEnum,&tau_zz[0],P1Enum));
    344 
    345         /*Clean up and return*/
    346         delete gauss;
    347 }
    348 /*}}}*/
    349 void       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         Input* levermanncoeff_input=inputs->GetInput(CalvinglevermannCoeffEnum);                     _assert_(levermanncoeff_input);
    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                 levermanncoeff_input->GetInputValue(&propcoeff,gauss);
    384 
    385                 /*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 */
    386                 calvingrate[iv]=propcoeff*strainparallel*strainperpendicular;   
    387                 if(calvingrate[iv]<0){
    388                         calvingrate[iv]=0;
    389                 }
    390                 calvingratex[iv]=calvingrate[iv]*vx/(vel+1.e-6);
    391                 calvingratey[iv]=calvingrate[iv]*vy/(vel+1.e-6);
    392         }
    393 
    394         /*Add input*/
    395         this->inputs->AddInput(new PentaInput(CalvingratexEnum,&calvingratex[0],P1Enum));
    396         this->inputs->AddInput(new PentaInput(CalvingrateyEnum,&calvingratey[0],P1Enum));
    397         this->inputs->AddInput(new PentaInput(CalvingCalvingrateEnum,&calvingrate[0],P1Enum));
    398 
    399         /*Clean up and return*/
    400         delete gauss;
    401 
    402 }
    403 /*}}}*/
    404 void       Penta::StressIntensityFactor(){/*{{{*/
    405 
    406         /* Check if we are on the base */
    407         if(!IsOnBase()) return;
    408 
    409         IssmDouble  ki[6]={0.};
    410         IssmDouble  const_grav=9.81;
    411         IssmDouble  rho_ice=900;
    412         IssmDouble  rho_water=1000;
    413         IssmDouble  Jdet[3];
    414         IssmDouble  pressure,vx,vy,vel,deviaxx,deviaxy,deviayy,water_depth,prof,stress_xx,thickness;
    415 
    416         Penta* penta=this;
    417         for(;;){
    418        
    419                 IssmDouble  xyz_list[NUMVERTICES][3];
    420                 /* Get node coordinates and dof list: */
    421                 ::GetVerticesCoordinates(&xyz_list[0][0],penta->vertices,NUMVERTICES);
    422 
    423                 ///*Compute the Jacobian for the vertical integration*/
    424                 Jdet[0]=(xyz_list[3][2]-xyz_list[0][2])*0.5;
    425                 Jdet[1]=(xyz_list[4][2]-xyz_list[1][2])*0.5;
    426                 Jdet[2]=(xyz_list[5][2]-xyz_list[2][2])*0.5;
    427        
    428                 /*Retrieve all inputs we will need*/
    429                 Input* vx_input=inputs->GetInput(VxEnum);                                  _assert_(vx_input);
    430                 Input* vy_input=inputs->GetInput(VyEnum);                                  _assert_(vy_input);
    431                 Input* vel_input=inputs->GetInput(VelEnum);                                _assert_(vel_input);
    432                 Input* pressure_input=inputs->GetInput(PressureEnum);                      _assert_(pressure_input);
    433                 Input* deviaxx_input=inputs->GetInput(DeviatoricStressxxEnum);             _assert_(deviaxx_input);
    434                 Input* deviaxy_input=inputs->GetInput(DeviatoricStressxyEnum);             _assert_(deviaxy_input);
    435                 Input* deviayy_input=inputs->GetInput(DeviatoricStressyyEnum);             _assert_(deviayy_input);
    436                 Input* surface_input=inputs->GetInput(SurfaceEnum);                                                             _assert_(surface_input);
    437                 Input* thickness_input=inputs->GetInput(ThicknessEnum);                                                 _assert_(thickness_input);
    438                
    439                 /* Start looping on the number of 2D vertices: */
    440                 for(int ig=0;ig<3;ig++){
    441                         GaussPenta* gauss=new GaussPenta(ig,3+ig,11);
    442                         for (int iv=gauss->begin();iv<gauss->end();iv++){
    443                                 gauss->GaussPoint(iv);
    444 
    445                                 /* Get the value we need*/
    446                                 pressure_input->GetInputValue(&pressure,gauss);
    447                                 vx_input->GetInputValue(&vx,gauss);
    448                                 vy_input->GetInputValue(&vy,gauss);
    449                                 vel_input->GetInputValue(&vel,gauss);
    450                                 deviaxx_input->GetInputValue(&deviaxx,gauss);
    451                                 deviaxy_input->GetInputValue(&deviaxy,gauss);
    452                                 deviayy_input->GetInputValue(&deviayy,gauss);
    453                                 surface_input->GetInputValue(&water_depth,gauss);
    454                                 thickness_input->GetInputValue(&thickness,gauss);
    455                                 prof=water_depth-penta->GetZcoord(&xyz_list[0][0],gauss);
    456 
    457                                 /*stress_xx= Deviatoric stress along the ice flow direction plus cryostatic pressure */
    458                                 stress_xx=(vx*vx*(deviaxx)+vy*vy*(deviayy)+2*vy*vx*deviaxy)/(vel*vel+1.e-6);
    459 
    460                                 if(prof<water_depth&prof<thickness){
    461                                         /* Compute the local stress intensity factor*/
    462                                         ki[ig]+=Jdet[ig]*gauss->weight*stress_xx*StressIntensityIntegralWeight(prof,min(water_depth,thickness),thickness);
    463                                 }
    464                         }
    465                         delete gauss;
    466                 }
    467                        
    468                 /*Stop if we have reached the surface/base*/
    469                 if(penta->IsOnSurface()) break;
    470                
    471                 /*get upper Penta*/
    472                 penta=penta->GetUpperPenta();
    473                 _assert_(penta->Id()!=this->id);
    474         }
    475 
    476         /*Add input*/
    477         this->inputs->AddInput(new PentaInput(StressIntensityFactorEnum,&ki[0],P1Enum));
    478         this->InputExtrude(StressIntensityFactorEnum,-1);
    479 }
    480 /*}}}*/
    481 void       Penta::StrainRateparallel(){/*{{{*/
    482 
    483         IssmDouble *xyz_list = NULL;
    484         IssmDouble  epsilon[6];
    485         GaussPenta* gauss=NULL;
    486         IssmDouble  vx,vy,vel;
    487         IssmDouble  strainxx;
    488         IssmDouble  strainxy;
    489         IssmDouble  strainyy;
    490         IssmDouble  strainparallel[NUMVERTICES];
    491 
    492         /* Get node coordinates and dof list: */
    493         this->GetVerticesCoordinates(&xyz_list);
    494 
    495         /*Retrieve all inputs we will need*/
    496         Input* vx_input=inputs->GetInput(VxEnum);                                  _assert_(vx_input);
    497         Input* vy_input=inputs->GetInput(VyEnum);                                  _assert_(vy_input);
    498         Input* vz_input=inputs->GetInput(VzEnum);                                                                                               _assert_(vz_input);
    499 
    500         /* Start looping on the number of vertices: */
    501         gauss=new GaussPenta();
    502         for (int iv=0;iv<NUMVERTICES;iv++){
    503                 gauss->GaussVertex(iv);
    504 
    505                 /* Get the value we need*/
    506                 vx_input->GetInputValue(&vx,gauss);
    507                 vy_input->GetInputValue(&vy,gauss);
    508                 vel=vx*vx+vy*vy;
    509 
    510                 /*Compute strain rate and viscosity: */
    511                 this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input);
    512                 strainxx=epsilon[0];
    513                 strainyy=epsilon[1];
    514                 strainxy=epsilon[3];
    515 
    516                 /*strainparallel= Strain rate along the ice flow direction */
    517                 strainparallel[iv]=(vx*vx*(strainxx)+vy*vy*(strainyy)+2*vy*vx*strainxy)/(vel+1.e-6);
    518         }
    519 
    520         /*Add input*/
    521         this->inputs->AddInput(new PentaInput(StrainRateparallelEnum,&strainparallel[0],P1Enum));
    522 
    523         /*Clean up and return*/
    524         delete gauss;
    525         xDelete<IssmDouble>(xyz_list);
    526 }
    527 /*}}}*/
    528 void       Penta::StrainRateperpendicular(){/*{{{*/
    529 
    530         IssmDouble *xyz_list = NULL;
    531         IssmDouble  epsilon[6];
    532         GaussPenta* gauss=NULL;
    533         IssmDouble  vx,vy,vel;
    534         IssmDouble  strainxx;
    535         IssmDouble  strainxy;
    536         IssmDouble  strainyy;
    537         IssmDouble  strainperpendicular[NUMVERTICES];
    538 
    539         /* Get node coordinates and dof list: */
    540         this->GetVerticesCoordinates(&xyz_list);
    541 
    542         /*Retrieve all inputs we will need*/
    543         Input* vx_input=inputs->GetInput(VxEnum);                                  _assert_(vx_input);
    544         Input* vy_input=inputs->GetInput(VyEnum);                                  _assert_(vy_input);
    545         Input* vz_input=inputs->GetInput(VzEnum);                                                                                               _assert_(vz_input);
    546 
    547         /* Start looping on the number of vertices: */
    548         gauss=new GaussPenta();
    549         for (int iv=0;iv<NUMVERTICES;iv++){
    550                 gauss->GaussVertex(iv);
    551 
    552                 /* Get the value we need*/
    553                 vx_input->GetInputValue(&vx,gauss);
    554                 vy_input->GetInputValue(&vy,gauss);
    555                 vel=vx*vx+vy*vy;
    556 
    557                 /*Compute strain rate and viscosity: */
    558                 this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input);
    559                 strainxx=epsilon[0];
    560                 strainyy=epsilon[1];
    561                 strainxy=epsilon[3];
    562 
    563                 /*strainperpendicular= Strain rate perpendicular to the ice flow direction */
    564                 strainperpendicular[iv]=(vx*vx*(strainyy)+vy*vy*(strainxx)-2*vy*vx*strainxy)/(vel+1.e-6);
    565         }
    566 
    567         /*Add input*/
    568         this->inputs->AddInput(new PentaInput(StrainRateperpendicularEnum,&strainperpendicular[0],P1Enum));
    569 
    570         /*Clean up and return*/
    571         delete gauss;
    572         xDelete<IssmDouble>(xyz_list);
    573 }
    574 /*}}}*/
    575408void       Penta::Configure(Elements* elementsin, Loads* loadsin, Nodes* nodesin,Vertices* verticesin, Materials* materialsin, Parameters* parametersin){/*{{{*/
    576409
    … …  
    606439}
    607440/*}}}*/
     441void       Penta::ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){/*{{{*/
     442
     443        int    vertexpidlist[NUMVERTICES];
     444        IssmDouble grad_list[NUMVERTICES];
     445        Input* grad_input=NULL;
     446        Input* input=NULL;
     447
     448        if(enum_type==MaterialsRheologyBbarEnum){
     449                input=(Input*)inputs->GetInput(MaterialsRheologyBEnum);
     450        }
     451        else if(enum_type==DamageDbarEnum){
     452                input=(Input*)inputs->GetInput(DamageDEnum);
     453        }
     454        else{
     455                input=inputs->GetInput(enum_type);
     456        }
     457        if (!input) _error_("Input " << EnumToStringx(enum_type) << " not found");
     458        if (input->ObjectEnum()!=ControlInputEnum) _error_("Input " << EnumToStringx(enum_type) << " is not a ControlInput");
     459
     460        GradientIndexing(&vertexpidlist[0],control_index);
     461        for(int i=0;i<NUMVERTICES;i++) grad_list[i]=gradient[vertexpidlist[i]];
     462        grad_input=new PentaInput(GradientEnum,grad_list,P1Enum);
     463        ((ControlInput*)input)->SetGradient(grad_input);
     464
     465}/*}}}*/
     466void       Penta::ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum){/*{{{*/
     467
     468        Input* input=NULL;
     469
     470        if(control_enum==MaterialsRheologyBbarEnum){
     471                input=(Input*)inputs->GetInput(MaterialsRheologyBEnum);
     472        }
     473        else if(control_enum==DamageDbarEnum){
     474                input=(Input*)inputs->GetInput(DamageDEnum);
     475        }
     476        else{
     477                input=inputs->GetInput(control_enum);
     478        }
     479        if (!input) _error_("Input " << EnumToStringx(control_enum) << " not found");
     480        if (input->ObjectEnum()!=ControlInputEnum) _error_("Input " << EnumToStringx(control_enum) << " is not a ControlInput");
     481
     482        int         sidlist[NUMVERTICES];
     483        int         connectivity[NUMVERTICES];
     484        IssmPDouble values[NUMVERTICES];
     485        IssmPDouble gradients[NUMVERTICES];
     486        IssmDouble  value,gradient;
     487
     488        this->GetVerticesConnectivityList(&connectivity[0]);
     489        this->GetVerticesSidList(&sidlist[0]);
     490
     491        GaussPenta* gauss=new GaussPenta();
     492        for (int iv=0;iv<NUMVERTICES;iv++){
     493                gauss->GaussVertex(iv);
     494
     495                ((ControlInput*)input)->GetInputValue(&value,gauss);
     496                ((ControlInput*)input)->GetGradientValue(&gradient,gauss);
     497
     498                values[iv]    = reCast<IssmPDouble>(value)/reCast<IssmPDouble>(connectivity[iv]);
     499                gradients[iv] = reCast<IssmPDouble>(gradient)/reCast<IssmPDouble>(connectivity[iv]);
     500        }
     501        delete gauss;
     502
     503        vector_control->SetValues(NUMVERTICES,&sidlist[0],&values[0],ADD_VAL);
     504        vector_gradient->SetValues(NUMVERTICES,&sidlist[0],&gradients[0],ADD_VAL);
     505
     506}/*}}}*/
    608507void       Penta::Delta18oParameterization(void){/*{{{*/
    609508
    … …  
    680579}
    681580/*}}}*/
     581void       Penta::ElementResponse(IssmDouble* presponse,int response_enum){/*{{{*/
     582
     583        switch(response_enum){
     584                case MaterialsRheologyBbarEnum:
     585                        *presponse=this->material->GetBbar();
     586                        break;
     587                case DamageDbarEnum:
     588                        *presponse=this->material->GetDbar();
     589                        break;
     590                case VelEnum:
     591                        {
     592
     593                                /*Get input:*/
     594                                IssmDouble vel;
     595                                Input* vel_input;
     596
     597                                vel_input=this->inputs->GetInput(VelEnum); _assert_(vel_input);
     598                                vel_input->GetInputAverage(&vel);
     599
     600                                /*Assign output pointers:*/
     601                                *presponse=vel;
     602                        }
     603                        break;
     604                default: 
     605                        _error_("Response type " << EnumToStringx(response_enum) << " not supported yet!");
     606        }
     607
     608}
     609/*}}}*/
     610void       Penta::ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){/*{{{*/
     611
     612        IssmDouble xyz_list[NUMVERTICES][3];
     613        IssmDouble xmin,ymin,zmin;
     614        IssmDouble xmax,ymax,zmax;
     615
     616        /*Get xyz list: */
     617        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     618        xmin=xyz_list[0][0]; xmax=xyz_list[0][0];
     619        ymin=xyz_list[0][1]; ymax=xyz_list[0][1];
     620        zmin=xyz_list[0][2]; zmax=xyz_list[0][2];
     621
     622        for(int i=1;i<NUMVERTICES;i++){
     623                if(xyz_list[i][0]<xmin) xmin=xyz_list[i][0];
     624                if(xyz_list[i][0]>xmax) xmax=xyz_list[i][0];
     625                if(xyz_list[i][1]<ymin) ymin=xyz_list[i][1];
     626                if(xyz_list[i][1]>ymax) ymax=xyz_list[i][1];
     627                if(xyz_list[i][2]<zmin) zmin=xyz_list[i][2];
     628                if(xyz_list[i][2]>zmax) zmax=xyz_list[i][2];
     629        }
     630
     631        *hx=xmax-xmin;
     632        *hy=ymax-ymin;
     633        *hz=zmax-zmin;
     634}
     635/*}}}*/
    682636int        Penta::FiniteElement(void){/*{{{*/
    683637        return this->element_type;
    … …  
    770724}
    771725/*}}}*/
    772 int        Penta::ObjectEnum(void){/*{{{*/
    773 
    774         return PentaEnum;
    775 
    776 }
    777 /*}}}*/
    778726void       Penta::GetAreaCoordinates(IssmDouble* area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints){/*{{{*/
    779727        /*Computeportion of the element that is grounded*/
    … …  
    811759}
    812760/*}}}*/
    813 Penta*     Penta::GetUpperPenta(void){/*{{{*/
    814 
    815         Penta* upper_penta=NULL;
    816 
    817         upper_penta=(Penta*)verticalneighbors[1]; //first one (0) under, second one (1) above
    818 
    819         return upper_penta;
    820 }
    821 /*}}}*/
    822 Penta*     Penta::GetLowerPenta(void){/*{{{*/
    823 
    824         Penta* lower_penta=NULL;
    825 
    826         lower_penta=(Penta*)verticalneighbors[0]; //first one (0) under, second one (1) above
    827 
    828         return lower_penta;
    829 }
    830 /*}}}*/
    831 Penta*     Penta::GetSurfacePenta(void){/*{{{*/
     761Element*   Penta::GetBasalElement(void){/*{{{*/
    832762
    833763        /*Output*/
    834         Penta* penta=NULL;
    835 
    836         /*Go through all pentas till the surface is reached*/
    837         penta=this;
    838         for(;;){
    839                 /*Stop if we have reached the surface, else, take upper penta*/
    840                 if (penta->IsOnSurface()) break;
    841 
    842                 /* get upper Penta*/
    843                 penta=penta->GetUpperPenta();
    844                 _assert_(penta->Id()!=this->id);
    845         }
    846 
    847         /*return output*/
    848         return penta;
     764        Element* element=this->GetBasalPenta();
     765        return element;
    849766}
    850767/*}}}*/
    … …  
    869786}
    870787/*}}}*/
    871 Element*   Penta::GetUpperElement(void){/*{{{*/
    872 
    873         /*Output*/
    874         Element* upper_element=this->GetUpperPenta();
    875         return upper_element;
    876 }
    877 /*}}}*/
    878 Element*   Penta::GetBasalElement(void){/*{{{*/
    879 
    880         /*Output*/
    881         Element* element=this->GetBasalPenta();
    882         return element;
     788int        Penta::GetElementType(){/*{{{*/
     789
     790        /*return PentaRef field*/
     791        return this->element_type;
    883792}
    884793/*}}}*/
    … …  
    1035944
    1036945        return phi;
    1037 }
    1038 /*}}}*/
    1039 int        Penta::GetElementType(){/*{{{*/
    1040 
    1041         /*return PentaRef field*/
    1042         return this->element_type;
    1043 }
    1044 /*}}}*/
    1045 void       Penta::ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){/*{{{*/
    1046 
    1047         IssmDouble xyz_list[NUMVERTICES][3];
    1048         IssmDouble xmin,ymin,zmin;
    1049         IssmDouble xmax,ymax,zmax;
    1050 
    1051         /*Get xyz list: */
    1052         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    1053         xmin=xyz_list[0][0]; xmax=xyz_list[0][0];
    1054         ymin=xyz_list[0][1]; ymax=xyz_list[0][1];
    1055         zmin=xyz_list[0][2]; zmax=xyz_list[0][2];
    1056 
    1057         for(int i=1;i<NUMVERTICES;i++){
    1058                 if(xyz_list[i][0]<xmin) xmin=xyz_list[i][0];
    1059                 if(xyz_list[i][0]>xmax) xmax=xyz_list[i][0];
    1060                 if(xyz_list[i][1]<ymin) ymin=xyz_list[i][1];
    1061                 if(xyz_list[i][1]>ymax) ymax=xyz_list[i][1];
    1062                 if(xyz_list[i][2]<zmin) zmin=xyz_list[i][2];
    1063                 if(xyz_list[i][2]>zmax) zmax=xyz_list[i][2];
    1064         }
    1065 
    1066         *hx=xmax-xmin;
    1067         *hy=ymax-ymin;
    1068         *hz=zmax-zmin;
    1069 }
    1070 /*}}}*/
    1071 int        Penta::GetNodeIndex(Node* node){/*{{{*/
    1072 
    1073         _assert_(nodes);
    1074         int numnodes = this->NumberofNodes(this->element_type);
    1075 
    1076         for(int i=0;i<numnodes;i++){
    1077                 if(node==nodes[i]) return i;
    1078         }
    1079         _error_("Node provided not found among element nodes");
    1080 
    1081 }
    1082 /*}}}*/
    1083 int        Penta::GetNumberOfNodes(void){/*{{{*/
    1084         return this->NumberofNodes(this->element_type);
    1085 }
    1086 /*}}}*/
    1087 int        Penta::GetNumberOfNodes(int enum_type){/*{{{*/
    1088         return this->NumberofNodes(enum_type);
    1089 }
    1090 /*}}}*/
    1091 int        Penta::GetNumberOfVertices(void){/*{{{*/
    1092         return NUMVERTICES;
    1093 }
    1094 /*}}}*/
    1095 Node*      Penta::GetNode(int node_number){/*{{{*/
    1096         _assert_(node_number>=0);
    1097         _assert_(node_number<this->NumberofNodes(this->element_type));
    1098         return this->nodes[node_number];
    1099 }
    1100 /*}}}*/
    1101 void       Penta::GetInputValue(IssmDouble* pvalue,Node* node,int enumtype){/*{{{*/
    1102 
    1103         Input* input=inputs->GetInput(enumtype);
    1104         if(!input) _error_("No input of type " << EnumToStringx(enumtype) << " found in tria");
    1105 
    1106         GaussPenta* gauss=new GaussPenta();
    1107         gauss->GaussVertex(this->GetNodeIndex(node));
    1108 
    1109         input->GetInputValue(pvalue,gauss);
    1110         delete gauss;
    1111 }
    1112 /*}}}*/
    1113 void       Penta::GetVerticesCoordinatesBase(IssmDouble** pxyz_list){/*{{{*/
    1114 
    1115         IssmDouble* xyz_list = xNew<IssmDouble>(NUMVERTICES2D*3);
    1116         ::GetVerticesCoordinates(xyz_list,this->vertices,NUMVERTICES2D);
    1117 
    1118         /*Assign output pointer*/
    1119         *pxyz_list = xyz_list;
    1120 
    1121 }/*}}}*/
    1122 void       Penta::GetVerticesCoordinatesTop(IssmDouble** pxyz_list){/*{{{*/
    1123 
    1124         IssmDouble* xyz_list = xNew<IssmDouble>(NUMVERTICES2D*3);
    1125         ::GetVerticesCoordinates(xyz_list,&this->vertices[3],NUMVERTICES2D);
    1126 
    1127         /*Assign output pointer*/
    1128         *pxyz_list = xyz_list;
    1129 
    1130 }/*}}}*/
    1131 void       Penta::NormalSection(IssmDouble* normal,IssmDouble* xyz_list){/*{{{*/
    1132 
    1133         /*Build unit outward pointing vector*/
    1134         IssmDouble AB[3];
    1135         IssmDouble AC[3];
    1136         IssmDouble norm;
    1137 
    1138         AB[0]=xyz_list[1*3+0] - xyz_list[0*3+0];
    1139         AB[1]=xyz_list[1*3+1] - xyz_list[0*3+1];
    1140         AB[2]=xyz_list[1*3+2] - xyz_list[0*3+2];
    1141         AC[0]=xyz_list[2*3+0] - xyz_list[0*3+0];
    1142         AC[1]=xyz_list[2*3+1] - xyz_list[0*3+1];
    1143         AC[2]=xyz_list[2*3+2] - xyz_list[0*3+2];
    1144 
    1145         cross(normal,AB,AC);
    1146         norm=sqrt(normal[0]*normal[0]+normal[1]*normal[1]+normal[2]*normal[2]);
    1147 
    1148         for(int i=0;i<3;i++) normal[i]=normal[i]/norm;
    1149 }
    1150 /*}}}*/
    1151 IssmDouble Penta::StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa){/*{{{*/
    1152         /*Compute stabilization parameter*/
    1153         /*kappa=thermalconductivity/(rho_ice*hearcapacity) for thermal model*/
    1154         /*kappa=enthalpydiffusionparameter for enthalpy model*/
    1155 
    1156         IssmDouble normu;
    1157         IssmDouble tau_parameter;
    1158 
    1159         normu=pow(pow(u,2)+pow(v,2)+pow(w,2),0.5);
    1160         if(normu*diameter/(3*2*kappa)<1){
    1161                 tau_parameter=pow(diameter,2)/(3*2*2*kappa);
    1162         }
    1163         else tau_parameter=diameter/(2*normu);
    1164 
    1165         return tau_parameter;
    1166 }
    1167 /*}}}*/
    1168 void       Penta::ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum){/*{{{*/
    1169         /*Compute portion of the element that is grounded*/
    1170 
    1171         int         normal_orientation=0;
    1172         IssmDouble  s1,s2;
    1173         IssmDouble  levelset[NUMVERTICES];
    1174         IssmDouble* xyz_zero = xNew<IssmDouble>(4*3);
    1175 
    1176         /*Recover parameters and values*/
    1177         GetInputListOnVertices(&levelset[0],levelsetenum);
    1178 
    1179         if(levelset[0]*levelset[1]>0.){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
    1180                 /*Portion of the segments*/
    1181                 s1=levelset[2]/(levelset[2]-levelset[1]);
    1182                 s2=levelset[2]/(levelset[2]-levelset[0]);
    1183 
    1184                 if(levelset[2]<0.) normal_orientation=1; //orientation of quadrangle at base and top, depending on distribution of levelsetfunction
    1185                 /*New point 1*/
    1186                 xyz_zero[3*normal_orientation+0]=xyz_list[2*3+0]+s1*(xyz_list[1*3+0]-xyz_list[2*3+0]);
    1187                 xyz_zero[3*normal_orientation+1]=xyz_list[2*3+1]+s1*(xyz_list[1*3+1]-xyz_list[2*3+1]);
    1188                 xyz_zero[3*normal_orientation+2]=xyz_list[2*3+2]+s1*(xyz_list[1*3+2]-xyz_list[2*3+2]);
    1189 
    1190                 /*New point 0*/
    1191                 xyz_zero[3*(1-normal_orientation)+0]=xyz_list[2*3+0]+s2*(xyz_list[0*3+0]-xyz_list[2*3+0]);
    1192                 xyz_zero[3*(1-normal_orientation)+1]=xyz_list[2*3+1]+s2*(xyz_list[0*3+1]-xyz_list[2*3+1]);
    1193                 xyz_zero[3*(1-normal_orientation)+2]=xyz_list[2*3+2]+s2*(xyz_list[0*3+2]-xyz_list[2*3+2]);
    1194 
    1195                 /*New point 3*/
    1196                 xyz_zero[3*(2+1-normal_orientation)+0]=xyz_list[5*3+0]+s1*(xyz_list[4*3+0]-xyz_list[5*3+0]);
    1197                 xyz_zero[3*(2+1-normal_orientation)+1]=xyz_list[5*3+1]+s1*(xyz_list[4*3+1]-xyz_list[5*3+1]);
    1198                 xyz_zero[3*(2+1-normal_orientation)+2]=xyz_list[5*3+2]+s1*(xyz_list[4*3+2]-xyz_list[5*3+2]);
    1199 
    1200                 /*New point 4*/
    1201                 xyz_zero[3*(2+normal_orientation)+0]=xyz_list[5*3+0]+s2*(xyz_list[3*3+0]-xyz_list[5*3+0]);
    1202                 xyz_zero[3*(2+normal_orientation)+1]=xyz_list[5*3+1]+s2*(xyz_list[3*3+1]-xyz_list[5*3+1]);
    1203                 xyz_zero[3*(2+normal_orientation)+2]=xyz_list[5*3+2]+s2*(xyz_list[3*3+2]-xyz_list[5*3+2]);
    1204         }
    1205         else if(levelset[1]*levelset[2]>0.){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
    1206                 /*Portion of the segments*/
    1207                 s1=levelset[0]/(levelset[0]-levelset[2]);
    1208                 s2=levelset[0]/(levelset[0]-levelset[1]);
    1209 
    1210                 if(levelset[0]<0.) normal_orientation=1;
    1211                 /*New point 1*/
    1212                 xyz_zero[3*normal_orientation+0]=xyz_list[0*3+0]+s1*(xyz_list[2*3+0]-xyz_list[0*3+0]);
    1213                 xyz_zero[3*normal_orientation+1]=xyz_list[0*3+1]+s1*(xyz_list[2*3+1]-xyz_list[0*3+1]);
    1214                 xyz_zero[3*normal_orientation+2]=xyz_list[0*3+2]+s1*(xyz_list[2*3+2]-xyz_list[0*3+2]);
    1215 
    1216                 /*New point 2*/
    1217                 xyz_zero[3*(1-normal_orientation)+0]=xyz_list[0*3+0]+s2*(xyz_list[1*3+0]-xyz_list[0*3+0]);
    1218                 xyz_zero[3*(1-normal_orientation)+1]=xyz_list[0*3+1]+s2*(xyz_list[1*3+1]-xyz_list[0*3+1]);
    1219                 xyz_zero[3*(1-normal_orientation)+2]=xyz_list[0*3+2]+s2*(xyz_list[1*3+2]-xyz_list[0*3+2]);
    1220 
    1221                 /*New point 3*/
    1222                 xyz_zero[3*(2+1-normal_orientation)+0]=xyz_list[3*3+0]+s1*(xyz_list[5*3+0]-xyz_list[3*3+0]);
    1223                 xyz_zero[3*(2+1-normal_orientation)+1]=xyz_list[3*3+1]+s1*(xyz_list[5*3+1]-xyz_list[3*3+1]);
    1224                 xyz_zero[3*(2+1-normal_orientation)+2]=xyz_list[3*3+2]+s1*(xyz_list[5*3+2]-xyz_list[3*3+2]);
    1225 
    1226                 /*New point 4*/
    1227                 xyz_zero[3*(2+normal_orientation)+0]=xyz_list[3*3+0]+s2*(xyz_list[4*3+0]-xyz_list[3*3+0]);
    1228                 xyz_zero[3*(2+normal_orientation)+1]=xyz_list[3*3+1]+s2*(xyz_list[4*3+1]-xyz_list[3*3+1]);
    1229                 xyz_zero[3*(2+normal_orientation)+2]=xyz_list[3*3+2]+s2*(xyz_list[4*3+2]-xyz_list[3*3+2]);
    1230         }
    1231         else if(levelset[0]*levelset[2]>0.){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
    1232                 /*Portion of the segments*/
    1233                 s1=levelset[1]/(levelset[1]-levelset[0]);
    1234                 s2=levelset[1]/(levelset[1]-levelset[2]);
    1235 
    1236                 if(levelset[1]<0.) normal_orientation=1;
    1237                 /*New point 0*/
    1238                 xyz_zero[3*normal_orientation+0]=xyz_list[1*3+0]+s1*(xyz_list[0*3+0]-xyz_list[1*3+0]);
    1239                 xyz_zero[3*normal_orientation+1]=xyz_list[1*3+1]+s1*(xyz_list[0*3+1]-xyz_list[1*3+1]);
    1240                 xyz_zero[3*normal_orientation+2]=xyz_list[1*3+2]+s1*(xyz_list[0*3+2]-xyz_list[1*3+2]);
    1241 
    1242                 /*New point 2*/
    1243                 xyz_zero[3*(1-normal_orientation)+0]=xyz_list[1*3+0]+s2*(xyz_list[2*3+0]-xyz_list[1*3+0]);
    1244                 xyz_zero[3*(1-normal_orientation)+1]=xyz_list[1*3+1]+s2*(xyz_list[2*3+1]-xyz_list[1*3+1]);
    1245                 xyz_zero[3*(1-normal_orientation)+2]=xyz_list[1*3+2]+s2*(xyz_list[2*3+2]-xyz_list[1*3+2]);
    1246 
    1247                 /*New point 3*/
    1248                 xyz_zero[3*(2+1-normal_orientation)+0]=xyz_list[4*3+0]+s1*(xyz_list[3*3+0]-xyz_list[4*3+0]);
    1249                 xyz_zero[3*(2+1-normal_orientation)+1]=xyz_list[4*3+1]+s1*(xyz_list[3*3+1]-xyz_list[4*3+1]);
    1250                 xyz_zero[3*(2+1-normal_orientation)+2]=xyz_list[4*3+2]+s1*(xyz_list[3*3+2]-xyz_list[4*3+2]);
    1251 
    1252                 /*New point 4*/
    1253                 xyz_zero[3*(2+normal_orientation)+0]=xyz_list[4*3+0]+s2*(xyz_list[5*3+0]-xyz_list[4*3+0]);
    1254                 xyz_zero[3*(2+normal_orientation)+1]=xyz_list[4*3+1]+s2*(xyz_list[5*3+1]-xyz_list[4*3+1]);
    1255                 xyz_zero[3*(2+normal_orientation)+2]=xyz_list[4*3+2]+s2*(xyz_list[5*3+2]-xyz_list[4*3+2]);
    1256         }
    1257         else if(levelset[0]==0. && levelset[1]==0.){ //front is on point 0 and 1
    1258                 xyz_zero[3*0+0]=xyz_list[0*3+0];
    1259                 xyz_zero[3*0+1]=xyz_list[0*3+1];
    1260                 xyz_zero[3*0+2]=xyz_list[0*3+2];
    1261 
    1262                 /*New point 2*/
    1263                 xyz_zero[3*1+0]=xyz_list[1*3+0];
    1264                 xyz_zero[3*1+1]=xyz_list[1*3+1];
    1265                 xyz_zero[3*1+2]=xyz_list[1*3+2];
    1266 
    1267                 /*New point 3*/
    1268                 xyz_zero[3*2+0]=xyz_list[4*3+0];
    1269                 xyz_zero[3*2+1]=xyz_list[4*3+1];
    1270                 xyz_zero[3*2+2]=xyz_list[4*3+2];
    1271 
    1272                 /*New point 4*/
    1273                 xyz_zero[3*3+0]=xyz_list[3*3+0];
    1274                 xyz_zero[3*3+1]=xyz_list[3*3+1];
    1275                 xyz_zero[3*3+2]=xyz_list[3*3+2];
    1276         }
    1277         else if(levelset[0]==0. && levelset[2]==0.){ //front is on point 0 and 1
    1278                 xyz_zero[3*0+0]=xyz_list[2*3+0];
    1279                 xyz_zero[3*0+1]=xyz_list[2*3+1];
    1280                 xyz_zero[3*0+2]=xyz_list[2*3+2];
    1281 
    1282                 /*New point 2*/
    1283                 xyz_zero[3*1+0]=xyz_list[0*3+0];
    1284                 xyz_zero[3*1+1]=xyz_list[0*3+1];
    1285                 xyz_zero[3*1+2]=xyz_list[0*3+2];
    1286 
    1287                 /*New point 3*/
    1288                 xyz_zero[3*2+0]=xyz_list[3*3+0];
    1289                 xyz_zero[3*2+1]=xyz_list[3*3+1];
    1290                 xyz_zero[3*2+2]=xyz_list[3*3+2];
    1291 
    1292                 /*New point 4*/
    1293                 xyz_zero[3*3+0]=xyz_list[5*3+0];
    1294                 xyz_zero[3*3+1]=xyz_list[5*3+1];
    1295                 xyz_zero[3*3+2]=xyz_list[5*3+2];
    1296         }
    1297         else if(levelset[1]==0. && levelset[2]==0.){ //front is on point 0 and 1
    1298                 xyz_zero[3*0+0]=xyz_list[1*3+0];
    1299                 xyz_zero[3*0+1]=xyz_list[1*3+1];
    1300                 xyz_zero[3*0+2]=xyz_list[1*3+2];
    1301 
    1302                 /*New point 2*/
    1303                 xyz_zero[3*1+0]=xyz_list[2*3+0];
    1304                 xyz_zero[3*1+1]=xyz_list[2*3+1];
    1305                 xyz_zero[3*1+2]=xyz_list[2*3+2];
    1306 
    1307                 /*New point 3*/
    1308                 xyz_zero[3*2+0]=xyz_list[5*3+0];
    1309                 xyz_zero[3*2+1]=xyz_list[5*3+1];
    1310                 xyz_zero[3*2+2]=xyz_list[5*3+2];
    1311 
    1312                 /*New point 4*/
    1313                 xyz_zero[3*3+0]=xyz_list[4*3+0];
    1314                 xyz_zero[3*3+1]=xyz_list[4*3+1];
    1315                 xyz_zero[3*3+2]=xyz_list[4*3+2];
    1316         }
    1317         else _error_("Case not covered");
    1318 
    1319         /*Assign output pointer*/
    1320         *pxyz_zero= xyz_zero;
    1321946}
    1322947/*}}}*/
    … …  
    1363988        xDelete<int>(indicesfront);
    1364989}/*}}}*/
     990void       Penta::GetInputValue(IssmDouble* pvalue,Node* node,int enumtype){/*{{{*/
     991
     992        Input* input=inputs->GetInput(enumtype);
     993        if(!input) _error_("No input of type " << EnumToStringx(enumtype) << " found in tria");
     994
     995        GaussPenta* gauss=new GaussPenta();
     996        gauss->GaussVertex(this->GetNodeIndex(node));
     997
     998        input->GetInputValue(pvalue,gauss);
     999        delete gauss;
     1000}
     1001/*}}}*/
     1002Penta*     Penta::GetLowerPenta(void){/*{{{*/
     1003
     1004        Penta* lower_penta=NULL;
     1005
     1006        lower_penta=(Penta*)verticalneighbors[0]; //first one (0) under, second one (1) above
     1007
     1008        return lower_penta;
     1009}
     1010/*}}}*/
     1011Node*      Penta::GetNode(int node_number){/*{{{*/
     1012        _assert_(node_number>=0);
     1013        _assert_(node_number<this->NumberofNodes(this->element_type));
     1014        return this->nodes[node_number];
     1015}
     1016/*}}}*/
     1017int        Penta::GetNodeIndex(Node* node){/*{{{*/
     1018
     1019        _assert_(nodes);
     1020        int numnodes = this->NumberofNodes(this->element_type);
     1021
     1022        for(int i=0;i<numnodes;i++){
     1023                if(node==nodes[i]) return i;
     1024        }
     1025        _error_("Node provided not found among element nodes");
     1026
     1027}
     1028/*}}}*/
     1029int        Penta::GetNumberOfNodes(void){/*{{{*/
     1030        return this->NumberofNodes(this->element_type);
     1031}
     1032/*}}}*/
     1033int        Penta::GetNumberOfNodes(int enum_type){/*{{{*/
     1034        return this->NumberofNodes(enum_type);
     1035}
     1036/*}}}*/
     1037int        Penta::GetNumberOfVertices(void){/*{{{*/
     1038        return NUMVERTICES;
     1039}
     1040/*}}}*/
     1041void       Penta::GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution, int enum_type){/*{{{*/
     1042
     1043        const int    numdof=NDOF1*NUMVERTICES;
     1044
     1045        int          i;
     1046        int*         doflist=NULL;
     1047        IssmDouble   values[numdof];
     1048        IssmDouble   enum_value;
     1049        GaussPenta   *gauss=NULL;
     1050
     1051        /*Get dof list: */
     1052        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     1053        Input* enum_input=inputs->GetInput(enum_type); _assert_(enum_input);
     1054
     1055        gauss=new GaussPenta();
     1056        for(i=0;i<NUMVERTICES;i++){
     1057                /*Recover temperature*/
     1058                gauss->GaussVertex(i);
     1059                enum_input->GetInputValue(&enum_value,gauss);
     1060                values[i]=enum_value;
     1061        }
     1062
     1063        /*Add value to global vector*/
     1064        solution->SetValues(numdof,doflist,values,INS_VAL);
     1065
     1066        /*Free ressources:*/
     1067        delete gauss;
     1068        xDelete<int>(doflist);
     1069}
     1070/*}}}*/
     1071Penta*     Penta::GetSurfacePenta(void){/*{{{*/
     1072
     1073        /*Output*/
     1074        Penta* penta=NULL;
     1075
     1076        /*Go through all pentas till the surface is reached*/
     1077        penta=this;
     1078        for(;;){
     1079                /*Stop if we have reached the surface, else, take upper penta*/
     1080                if (penta->IsOnSurface()) break;
     1081
     1082                /* get upper Penta*/
     1083                penta=penta->GetUpperPenta();
     1084                _assert_(penta->Id()!=this->id);
     1085        }
     1086
     1087        /*return output*/
     1088        return penta;
     1089}
     1090/*}}}*/
     1091Element*   Penta::GetUpperElement(void){/*{{{*/
     1092
     1093        /*Output*/
     1094        Element* upper_element=this->GetUpperPenta();
     1095        return upper_element;
     1096}
     1097/*}}}*/
     1098Penta*     Penta::GetUpperPenta(void){/*{{{*/
     1099
     1100        Penta* upper_penta=NULL;
     1101
     1102        upper_penta=(Penta*)verticalneighbors[1]; //first one (0) under, second one (1) above
     1103
     1104        return upper_penta;
     1105}
     1106/*}}}*/
     1107void       Penta::GetVectorFromControlInputs(Vector<IssmDouble>* vector,int control_enum,int control_index,const char* data,bool onsid){/*{{{*/
     1108
     1109        int vertexidlist[NUMVERTICES];
     1110
     1111        /*Get out if this is not an element input*/
     1112        if(!IsInput(control_enum)) return;
     1113
     1114        /*Prepare index list*/
     1115        GradientIndexing(&vertexidlist[0],control_index,onsid);
     1116
     1117        /*Get input (either in element or material)*/
     1118        if(control_enum==MaterialsRheologyBbarEnum) control_enum=MaterialsRheologyBEnum;
     1119        Input* input=inputs->GetInput(control_enum);
     1120        if(!input) _error_("Input " << EnumToStringx(control_enum) << " not found in element");
     1121
     1122        /*Check that it is a ControlInput*/
     1123        if (input->ObjectEnum()!=ControlInputEnum){
     1124                _error_("input " << EnumToStringx(control_enum) << " is not a ControlInput");
     1125        }
     1126
     1127        ((ControlInput*)input)->GetVectorFromInputs(vector,&vertexidlist[0],data);
     1128}
     1129/*}}}*/
     1130void       Penta::GetVerticesCoordinatesBase(IssmDouble** pxyz_list){/*{{{*/
     1131
     1132        IssmDouble* xyz_list = xNew<IssmDouble>(NUMVERTICES2D*3);
     1133        ::GetVerticesCoordinates(xyz_list,this->vertices,NUMVERTICES2D);
     1134
     1135        /*Assign output pointer*/
     1136        *pxyz_list = xyz_list;
     1137
     1138}/*}}}*/
     1139void       Penta::GetVerticesCoordinatesTop(IssmDouble** pxyz_list){/*{{{*/
     1140
     1141        IssmDouble* xyz_list = xNew<IssmDouble>(NUMVERTICES2D*3);
     1142        ::GetVerticesCoordinates(xyz_list,&this->vertices[3],NUMVERTICES2D);
     1143
     1144        /*Assign output pointer*/
     1145        *pxyz_list = xyz_list;
     1146
     1147}/*}}}*/
     1148IssmDouble Penta::IceVolume(void){/*{{{*/
     1149
     1150        /*The volume of a troncated prism is base * 1/3 sum(length of edges)*/
     1151        IssmDouble base,height;
     1152        IssmDouble xyz_list[NUMVERTICES][3];
     1153
     1154        if(!IsIceInElement())return 0;
     1155
     1156        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     1157
     1158        /*First calculate the area of the base (cross section triangle)
     1159         * http://en.wikipedia.org/wiki/Pentangle
     1160         * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
     1161        base = 1./2.*fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1]));
     1162
     1163        /*Now get the average height*/
     1164        height = 1./3.*((xyz_list[3][2]-xyz_list[0][2])+(xyz_list[4][2]-xyz_list[1][2])+(xyz_list[5][2]-xyz_list[2][2]));
     1165
     1166        /*Return: */
     1167        return base*height;
     1168}
     1169/*}}}*/
     1170IssmDouble Penta::IceVolumeAboveFloatation(void){/*{{{*/
     1171
     1172        /*Volume above floatation: H + rho_water/rho_ice*bathymetry for nodes on the bed*/
     1173        IssmDouble rho_ice,rho_water;
     1174        IssmDouble base,bed,surface,bathymetry;
     1175        IssmDouble xyz_list[NUMVERTICES][3];
     1176
     1177        if(!IsIceInElement() || IsFloating() || !IsOnBase())return 0;
     1178
     1179        rho_ice=matpar->GetRhoIce();
     1180        rho_water=matpar->GetRhoWater();
     1181        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     1182
     1183        /*First calculate the area of the base (cross section triangle)
     1184         * http://en.wikipedia.org/wiki/Pentangle
     1185         * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
     1186        base = 1./2.*fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1]));
     1187
     1188        /*Now get the average height above floatation*/
     1189        Input* surface_input    = inputs->GetInput(SurfaceEnum);    _assert_(surface_input);
     1190        Input* base_input        = inputs->GetInput(BaseEnum);        _assert_(base_input);
     1191        Input* bed_input = inputs->GetInput(BedEnum); _assert_(bed_input);
     1192        surface_input->GetInputAverage(&surface);
     1193        base_input->GetInputAverage(&bed);
     1194        bed_input->GetInputAverage(&bathymetry);
     1195
     1196        /*Return: */
     1197        return base*(surface - bed + min( rho_water/rho_ice * bathymetry, 0.) );
     1198}
     1199/*}}}*/
     1200void       Penta::InputControlUpdate(IssmDouble scalar,bool save_parameter){/*{{{*/
     1201
     1202        /*Intermediary*/
     1203        int    num_controls;
     1204        int*   control_type=NULL;
     1205        Input* input=NULL;
     1206
     1207        /*retrieve some parameters: */
     1208        this->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
     1209        this->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
     1210
     1211        for(int i=0;i<num_controls;i++){
     1212
     1213                if(control_type[i]==MaterialsRheologyBbarEnum){
     1214                        if (!IsOnBase()) goto cleanup_and_return;
     1215                        input=(Input*)this->inputs->GetInput(MaterialsRheologyBEnum); _assert_(input);
     1216                }
     1217                else if(control_type[i]==DamageDbarEnum){
     1218                        if (!IsOnBase()) goto cleanup_and_return;
     1219                        input=(Input*)this->inputs->GetInput(DamageDEnum); _assert_(input);
     1220                }
     1221                else{
     1222                        input=(Input*)this->inputs->GetInput(control_type[i]); _assert_(input);
     1223                }
     1224                if(input->ObjectEnum()!=ControlInputEnum) _error_("input " << EnumToStringx(control_type[i]) << " is not a ControlInput");
     1225
     1226                ((ControlInput*)input)->UpdateValue(scalar);
     1227                ((ControlInput*)input)->Constrain();
     1228                if (save_parameter) ((ControlInput*)input)->SaveValue();
     1229
     1230                if(control_type[i]==MaterialsRheologyBbarEnum){
     1231                        this->InputExtrude(MaterialsRheologyBEnum,-1);
     1232                }
     1233                else if(control_type[i]==DamageDbarEnum){
     1234                        this->InputExtrude(DamageDEnum,-1);
     1235                }
     1236        }
     1237
     1238        /*Clean up and return*/
     1239cleanup_and_return:
     1240        xDelete<int>(control_type);
     1241}
     1242/*}}}*/
    13651243void       Penta::InputDepthAverageAtBase(int enum_type,int average_enum_type){/*{{{*/
    13661244
    … …  
    17321610}
    17331611/*}}}*/
     1612bool       Penta::IsIcefront(void){/*{{{*/
     1613
     1614        bool isicefront;
     1615        int i,nrice;
     1616   IssmDouble ls[NUMVERTICES];
     1617
     1618        /*Retrieve all inputs and parameters*/
     1619        GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum);
     1620
     1621        /* If only one vertex has ice, there is an ice front here */
     1622        isicefront=false;
     1623        if(IsIceInElement()){
     1624                nrice=0;       
     1625                for(i=0;i<NUMVERTICES2D;i++)
     1626                        if(ls[i]<0.) nrice++;
     1627                if(nrice==1) isicefront= true;
     1628        }
     1629        return isicefront;
     1630}/*}}}*/
     1631bool       Penta::IsNodeOnShelfFromFlags(IssmDouble* flags){/*{{{*/
     1632
     1633        int  i;
     1634        bool shelf=false;
     1635
     1636        for(i=0;i<NUMVERTICES;i++){
     1637                if (flags[vertices[i]->Pid()]<0.){
     1638                        shelf=true;
     1639                        break;
     1640                }
     1641        }
     1642        return shelf;
     1643}
     1644/*}}}*/
    17341645bool       Penta::IsOnBase(void){/*{{{*/
    17351646
    … …  
    17681679}
    17691680/*}}}*/
    1770 bool       Penta::IsNodeOnShelfFromFlags(IssmDouble* flags){/*{{{*/
    1771 
    1772         int  i;
    1773         bool shelf=false;
    1774 
    1775         for(i=0;i<NUMVERTICES;i++){
    1776                 if (flags[vertices[i]->Pid()]<0.){
    1777                         shelf=true;
    1778                         break;
    1779                 }
    1780         }
    1781         return shelf;
     1681bool       Penta::IsZeroLevelset(int levelset_enum){/*{{{*/
     1682
     1683        bool        iszerols;
     1684        IssmDouble  ls[NUMVERTICES];
     1685
     1686        /*Retrieve all inputs and parameters*/
     1687        GetInputListOnVertices(&ls[0],levelset_enum);
     1688
     1689        /*If the level set has always same sign, there is no ice front here*/
     1690        iszerols = false;
     1691        if(IsIceInElement()){
     1692                if(ls[0]*ls[1]<0. || ls[0]*ls[2]<0. || (ls[0]*ls[1]*ls[2]==0. && ls[0]*ls[1]+ls[0]*ls[2]+ls[1]*ls[2]<=0.)){
     1693                        iszerols = true;
     1694                }
     1695        }
     1696        return iszerols;
    17821697}
    17831698/*}}}*/
    … …  
    18031718}
    18041719/*}}}*/
     1720void       Penta::JacobianDeterminantSurface(IssmDouble* pJdet,IssmDouble* xyz_list_quad,Gauss* gauss){/*{{{*/
     1721
     1722        _assert_(gauss->Enum()==GaussPentaEnum);
     1723        this->GetQuadJacobianDeterminant(pJdet,xyz_list_quad,(GaussPenta*)gauss);
     1724
     1725}
     1726/*}}}*/
    18051727void       Penta::JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_top,Gauss* gauss){/*{{{*/
    18061728
    … …  
    18101732}
    18111733/*}}}*/
    1812 void       Penta::JacobianDeterminantSurface(IssmDouble* pJdet,IssmDouble* xyz_list_quad,Gauss* gauss){/*{{{*/
    1813 
    1814         _assert_(gauss->Enum()==GaussPentaEnum);
    1815         this->GetQuadJacobianDeterminant(pJdet,xyz_list_quad,(GaussPenta*)gauss);
    1816 
     1734IssmDouble Penta::MassFlux(IssmDouble* segment){/*{{{*/
     1735
     1736        IssmDouble mass_flux=0;
     1737
     1738        if(!IsOnBase()) return mass_flux;
     1739
     1740        /*Depth Averaging Vx and Vy*/
     1741        this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
     1742        this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
     1743
     1744        /*Spawn Tria element from the base of the Penta: */
     1745        Tria* tria=(Tria*)SpawnTria(0,1,2);
     1746        mass_flux=tria->MassFlux(segment);
     1747        delete tria->material; delete tria;
     1748
     1749        /*Delete Vx and Vy averaged*/
     1750        this->inputs->DeleteInput(VxAverageEnum);
     1751        this->inputs->DeleteInput(VyAverageEnum);
     1752
     1753        /*clean up and return*/
     1754        return mass_flux;
     1755}
     1756/*}}}*/
     1757IssmDouble Penta::MassFlux(IssmDouble x1, IssmDouble y1, IssmDouble x2, IssmDouble y2,int segment_id){/*{{{*/
     1758
     1759        IssmDouble mass_flux=0;
     1760
     1761        if(!IsOnBase()) return mass_flux;
     1762
     1763        /*Depth Averaging Vx and Vy*/
     1764        this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
     1765        this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
     1766
     1767        /*Spawn Tria element from the base of the Penta: */
     1768        Tria* tria=(Tria*)SpawnTria(0,1,2);
     1769        mass_flux=tria->MassFlux(x1,y1,x2,y2,segment_id);
     1770        delete tria->material; delete tria;
     1771
     1772        /*Delete Vx and Vy averaged*/
     1773        this->inputs->DeleteInput(VxAverageEnum);
     1774        this->inputs->DeleteInput(VyAverageEnum);
     1775
     1776        /*clean up and return*/
     1777        return mass_flux;
    18171778}
    18181779/*}}}*/
    … …  
    18361797
    18371798        return minlength;
     1799}
     1800/*}}}*/
     1801Gauss*     Penta::NewGauss(void){/*{{{*/
     1802        return new GaussPenta();
     1803}
     1804/*}}}*/
     1805Gauss*     Penta::NewGauss(int order){/*{{{*/
     1806        return new GaussPenta(order,order);
     1807}
     1808/*}}}*/
     1809Gauss*     Penta::NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert){/*{{{*/
     1810
     1811        IssmDouble  area_coordinates[4][3];
     1812
     1813        GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,xyz_list,4);
     1814
     1815        return new GaussPenta(area_coordinates,order_horiz,order_vert);
     1816}
     1817/*}}}*/
     1818Gauss*     Penta::NewGaussBase(int order){/*{{{*/
     1819        return new GaussPenta(0,1,2,order);
     1820}
     1821/*}}}*/
     1822Gauss*     Penta::NewGaussLine(int vertex1,int vertex2,int order){/*{{{*/
     1823        return new GaussPenta(vertex1,vertex2,order);
     1824}
     1825/*}}}*/
     1826Gauss*     Penta::NewGaussTop(int order){/*{{{*/
     1827        return new GaussPenta(3,4,5,order);
     1828}
     1829/*}}}*/
     1830void       Penta::NodalFunctions(IssmDouble* basis, Gauss* gauss){/*{{{*/
     1831
     1832        _assert_(gauss->Enum()==GaussPentaEnum);
     1833        this->GetNodalFunctions(basis,(GaussPenta*)gauss,this->element_type);
     1834
     1835}
     1836/*}}}*/
     1837void       Penta::NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     1838
     1839        _assert_(gauss->Enum()==GaussPentaEnum);
     1840        this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussPenta*)gauss,this->element_type);
     1841
     1842}
     1843/*}}}*/
     1844void       Penta::NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     1845
     1846        _assert_(gauss->Enum()==GaussPentaEnum);
     1847        this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussPenta*)gauss,this->VelocityInterpolation());
     1848
     1849}
     1850/*}}}*/
     1851void       Penta::NodalFunctionsMINIDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     1852
     1853        _assert_(gauss->Enum()==GaussPentaEnum);
     1854        this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussPenta*)gauss,P1bubbleEnum);
     1855
     1856}
     1857/*}}}*/
     1858void       Penta::NodalFunctionsPressure(IssmDouble* basis, Gauss* gauss){/*{{{*/
     1859
     1860        _assert_(gauss->Enum()==GaussPentaEnum);
     1861        this->GetNodalFunctions(basis,(GaussPenta*)gauss,this->PressureInterpolation());
     1862
     1863}
     1864/*}}}*/
     1865void       Penta::NodalFunctionsP1(IssmDouble* basis, Gauss* gauss){/*{{{*/
     1866
     1867        _assert_(gauss->Enum()==GaussPentaEnum);
     1868        this->GetNodalFunctions(basis,(GaussPenta*)gauss,P1Enum);
     1869
     1870}
     1871/*}}}*/
     1872void       Penta::NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     1873
     1874        _assert_(gauss->Enum()==GaussPentaEnum);
     1875        this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussPenta*)gauss,P1Enum);
     1876
     1877}
     1878/*}}}*/
     1879void       Penta::NodalFunctionsP2(IssmDouble* basis, Gauss* gauss){/*{{{*/
     1880
     1881        _assert_(gauss->Enum()==GaussPentaEnum);
     1882        this->GetNodalFunctions(basis,(GaussPenta*)gauss,P2Enum);
     1883
     1884}
     1885/*}}}*/
     1886void       Penta::NodalFunctionsVelocity(IssmDouble* basis, Gauss* gauss){/*{{{*/
     1887
     1888        _assert_(gauss->Enum()==GaussPentaEnum);
     1889        this->GetNodalFunctions(basis,(GaussPenta*)gauss,this->VelocityInterpolation());
     1890
     1891}
     1892/*}}}*/
     1893void       Penta::NodalFunctionsTensor(IssmDouble* basis, Gauss* gauss){/*{{{*/
     1894
     1895        _assert_(gauss->Enum()==GaussPentaEnum);
     1896        this->GetNodalFunctions(basis,(GaussPenta*)gauss,this->TensorInterpolation());
     1897
    18381898}
    18391899/*}}}*/
    … …  
    18691929}
    18701930/*}}}*/
    1871 Gauss*     Penta::NewGauss(void){/*{{{*/
    1872         return new GaussPenta();
    1873 }
    1874 /*}}}*/
    1875 Gauss*     Penta::NewGauss(int order){/*{{{*/
    1876         return new GaussPenta(order,order);
    1877 }
    1878 /*}}}*/
    1879 Gauss*     Penta::NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert){/*{{{*/
    1880 
    1881         IssmDouble  area_coordinates[4][3];
    1882 
    1883         GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,xyz_list,4);
    1884 
    1885         return new GaussPenta(area_coordinates,order_horiz,order_vert);
    1886 }
    1887 /*}}}*/
    1888 Gauss*     Penta::NewGaussBase(int order){/*{{{*/
    1889         return new GaussPenta(0,1,2,order);
    1890 }
    1891 /*}}}*/
    1892 Gauss*     Penta::NewGaussLine(int vertex1,int vertex2,int order){/*{{{*/
    1893         return new GaussPenta(vertex1,vertex2,order);
    1894 }
    1895 /*}}}*/
    1896 Gauss*     Penta::NewGaussTop(int order){/*{{{*/
    1897         return new GaussPenta(3,4,5,order);
    1898 }
    1899 /*}}}*/
    1900 void       Penta::NodalFunctions(IssmDouble* basis, Gauss* gauss){/*{{{*/
    1901 
    1902         _assert_(gauss->Enum()==GaussPentaEnum);
    1903         this->GetNodalFunctions(basis,(GaussPenta*)gauss,this->element_type);
    1904 
    1905 }
    1906 /*}}}*/
    1907 void       Penta::NodalFunctionsP1(IssmDouble* basis, Gauss* gauss){/*{{{*/
    1908 
    1909         _assert_(gauss->Enum()==GaussPentaEnum);
    1910         this->GetNodalFunctions(basis,(GaussPenta*)gauss,P1Enum);
    1911 
    1912 }
    1913 /*}}}*/
    1914 void       Penta::NodalFunctionsP2(IssmDouble* basis, Gauss* gauss){/*{{{*/
    1915 
    1916         _assert_(gauss->Enum()==GaussPentaEnum);
    1917         this->GetNodalFunctions(basis,(GaussPenta*)gauss,P2Enum);
    1918 
    1919 }
    1920 /*}}}*/
    1921 void       Penta::NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    1922 
    1923         _assert_(gauss->Enum()==GaussPentaEnum);
    1924         this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussPenta*)gauss,this->element_type);
    1925 
    1926 }
    1927 /*}}}*/
    1928 void       Penta::NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    1929 
    1930         _assert_(gauss->Enum()==GaussPentaEnum);
    1931         this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussPenta*)gauss,P1Enum);
    1932 
    1933 }
    1934 /*}}}*/
    1935 void       Penta::NodalFunctionsMINIDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    1936 
    1937         _assert_(gauss->Enum()==GaussPentaEnum);
    1938         this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussPenta*)gauss,P1bubbleEnum);
    1939 
    1940 }
    1941 /*}}}*/
    1942 void       Penta::NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    1943 
    1944         _assert_(gauss->Enum()==GaussPentaEnum);
    1945         this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussPenta*)gauss,this->VelocityInterpolation());
    1946 
    1947 }
    1948 /*}}}*/
    1949 void       Penta::NodalFunctionsVelocity(IssmDouble* basis, Gauss* gauss){/*{{{*/
    1950 
    1951         _assert_(gauss->Enum()==GaussPentaEnum);
    1952         this->GetNodalFunctions(basis,(GaussPenta*)gauss,this->VelocityInterpolation());
    1953 
    1954 }
    1955 /*}}}*/
    1956 void       Penta::NodalFunctionsPressure(IssmDouble* basis, Gauss* gauss){/*{{{*/
    1957 
    1958         _assert_(gauss->Enum()==GaussPentaEnum);
    1959         this->GetNodalFunctions(basis,(GaussPenta*)gauss,this->PressureInterpolation());
    1960 
    1961 }
    1962 /*}}}*/
    1963 void       Penta::NodalFunctionsTensor(IssmDouble* basis, Gauss* gauss){/*{{{*/
    1964 
    1965         _assert_(gauss->Enum()==GaussPentaEnum);
    1966         this->GetNodalFunctions(basis,(GaussPenta*)gauss,this->TensorInterpolation());
    1967 
    1968 }
    1969 /*}}}*/
    19701931void       Penta::NormalBase(IssmDouble* bed_normal,IssmDouble* xyz_list){/*{{{*/
    19711932
    … …  
    19901951}
    19911952/*}}}*/
     1953void       Penta::NormalSection(IssmDouble* normal,IssmDouble* xyz_list){/*{{{*/
     1954
     1955        /*Build unit outward pointing vector*/
     1956        IssmDouble AB[3];
     1957        IssmDouble AC[3];
     1958        IssmDouble norm;
     1959
     1960        AB[0]=xyz_list[1*3+0] - xyz_list[0*3+0];
     1961        AB[1]=xyz_list[1*3+1] - xyz_list[0*3+1];
     1962        AB[2]=xyz_list[1*3+2] - xyz_list[0*3+2];
     1963        AC[0]=xyz_list[2*3+0] - xyz_list[0*3+0];
     1964        AC[1]=xyz_list[2*3+1] - xyz_list[0*3+1];
     1965        AC[2]=xyz_list[2*3+2] - xyz_list[0*3+2];
     1966
     1967        cross(normal,AB,AC);
     1968        norm=sqrt(normal[0]*normal[0]+normal[1]*normal[1]+normal[2]*normal[2]);
     1969
     1970        for(int i=0;i<3;i++) normal[i]=normal[i]/norm;
     1971}
     1972/*}}}*/
    19921973void       Penta::NormalTop(IssmDouble* top_normal,IssmDouble* xyz_list){/*{{{*/
    19931974
    … …  
    20181999int        Penta::NumberofNodesVelocity(void){/*{{{*/
    20192000        return PentaRef::NumberofNodes(this->VelocityInterpolation());
     2001}
     2002/*}}}*/
     2003int        Penta::ObjectEnum(void){/*{{{*/
     2004
     2005        return PentaEnum;
     2006
    20202007}
    20212008/*}}}*/
    … …  
    20892076}
    20902077/*}}}*/
     2078void       Penta::PotentialUngrounding(Vector<IssmDouble>* potential_ungrounding){/*{{{*/
     2079
     2080        IssmDouble  h[NUMVERTICES],r[NUMVERTICES],gl[NUMVERTICES];
     2081        IssmDouble  bed_hydro;
     2082        IssmDouble  rho_water,rho_ice,density;
     2083
     2084        /*material parameters: */
     2085        rho_water=matpar->GetRhoWater();
     2086        rho_ice=matpar->GetRhoIce();
     2087        density=rho_ice/rho_water;
     2088        GetInputListOnVertices(&h[0],ThicknessEnum);
     2089        GetInputListOnVertices(&r[0],BedEnum);
     2090        GetInputListOnVertices(&gl[0],MaskGroundediceLevelsetEnum);
     2091
     2092        /*go through vertices, and figure out which ones are on the ice sheet, and want to unground: */
     2093        for(int i=0;i<NUMVERTICES;i++){
     2094                /*Find if grounded vertices want to start floating*/
     2095                if (gl[i]>0.){
     2096                        bed_hydro=-density*h[i];
     2097                        if(bed_hydro>r[i]){
     2098                                /*Vertex that could potentially unground, flag it*/
     2099                                potential_ungrounding->SetValue(vertices[i]->Pid(),1,INS_VAL);
     2100                        }
     2101                }
     2102        }
     2103}
     2104/*}}}*/
     2105int        Penta::PressureInterpolation(void){/*{{{*/
     2106        return PentaRef::PressureInterpolation(this->element_type);
     2107}
     2108/*}}}*/
    20912109void       Penta::ReduceMatrices(ElementMatrix* Ke,ElementVector* pe){/*{{{*/
    20922110
    … …  
    22072225}
    22082226/*}}}*/
    2209 void       Penta::SetClone(int* minranks){/*{{{*/
    2210 
    2211         _error_("not implemented yet");
    2212 }
    2213 /*}}}*/
    2214 void       Penta::SetCurrentConfiguration(Elements* elementsin, Loads* loadsin, Nodes* nodesin, Materials* materialsin, Parameters* parametersin){/*{{{*/
    2215 
    2216         /*go into parameters and get the analysis_counter: */
    2217         int analysis_counter;
    2218         parametersin->FindParam(&analysis_counter,AnalysisCounterEnum);
    2219 
    2220         /*Get Element type*/
    2221         this->element_type=this->element_type_list[analysis_counter];
    2222 
    2223         /*Pick up nodes*/
    2224         if(this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp();
    2225         else this->nodes=NULL;
    2226 
    2227 }
    2228 /*}}}*/
    22292227void       Penta::ResetHooks(){/*{{{*/
    22302228
    … …  
    22452243}
    22462244/*}}}*/
     2245void       Penta::SetClone(int* minranks){/*{{{*/
     2246
     2247        _error_("not implemented yet");
     2248}
     2249/*}}}*/
     2250void       Penta::SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index){/*{{{*/
     2251
     2252        IssmDouble  values[NUMVERTICES];
     2253        int         vertexpidlist[NUMVERTICES],control_init;
     2254
     2255        /*Specific case for depth averaged quantities*/
     2256        control_init=control_enum;
     2257        if(control_enum==MaterialsRheologyBbarEnum){
     2258                control_enum=MaterialsRheologyBEnum;
     2259                if(!IsOnBase()) return;
     2260        }
     2261        if(control_enum==DamageDbarEnum){
     2262                control_enum=DamageDEnum;
     2263                if(!IsOnBase()) return;
     2264        }
     2265
     2266        /*Get out if this is not an element input*/
     2267        if(!IsInput(control_enum)) return;
     2268
     2269        /*Prepare index list*/
     2270        GradientIndexing(&vertexpidlist[0],control_index);
     2271
     2272        /*Get values on vertices*/
     2273        for(int i=0;i<NUMVERTICES;i++){
     2274                values[i]=vector[vertexpidlist[i]];
     2275        }
     2276        Input* new_input = new PentaInput(control_enum,values,P1Enum);
     2277        Input* input=(Input*)this->inputs->GetInput(control_enum);   _assert_(input);
     2278        if(input->ObjectEnum()!=ControlInputEnum){
     2279                _error_("input " << EnumToStringx(control_enum) << " is not a ControlInput");
     2280        }
     2281
     2282        ((ControlInput*)input)->SetInput(new_input);
     2283
     2284        if(control_init==MaterialsRheologyBbarEnum){
     2285                this->InputExtrude(control_enum,-1);
     2286        }
     2287        if(control_init==DamageDbarEnum){
     2288                this->InputExtrude(control_enum,-1);
     2289        }
     2290}
     2291/*}}}*/
     2292void       Penta::SetCurrentConfiguration(Elements* elementsin, Loads* loadsin, Nodes* nodesin, Materials* materialsin, Parameters* parametersin){/*{{{*/
     2293
     2294        /*go into parameters and get the analysis_counter: */
     2295        int analysis_counter;
     2296        parametersin->FindParam(&analysis_counter,AnalysisCounterEnum);
     2297
     2298        /*Get Element type*/
     2299        this->element_type=this->element_type_list[analysis_counter];
     2300
     2301        /*Pick up nodes*/
     2302        if(this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp();
     2303        else this->nodes=NULL;
     2304
     2305}
     2306/*}}}*/
    22472307void       Penta::SetTemporaryElementType(int element_type_in){/*{{{*/
    22482308        this->element_type=element_type_in;
    2249 }
    2250 /*}}}*/
    2251 Tria*      Penta::SpawnTria(int index1,int index2,int index3){/*{{{*/
    2252 
    2253         int analysis_counter;
    2254 
    2255         /*go into parameters and get the analysis_counter: */
    2256         this->parameters->FindParam(&analysis_counter,AnalysisCounterEnum);
    2257 
    2258         /*Create Tria*/
    2259         Tria* tria=new Tria();
    2260         tria->id=this->id;
    2261         tria->inputs=(Inputs*)this->inputs->SpawnTriaInputs(index1,index2,index3);
    2262         tria->parameters=this->parameters;
    2263         tria->element_type=P1Enum; //Only P1 CG for now (TO BE CHANGED)
    2264         this->SpawnTriaHook(xDynamicCast<ElementHook*>(tria),index1,index2,index3);
    2265 
    2266         /*Spawn material*/
    2267         tria->material=(Material*)this->material->copy2(tria);
    2268 
    2269         /*recover nodes, material and matpar: */
    2270         tria->nodes=(Node**)tria->hnodes[analysis_counter]->deliverp();
    2271         tria->vertices=(Vertex**)tria->hvertices->deliverp();
    2272         tria->matpar=(Matpar*)tria->hmatpar->delivers();
    2273 
    2274         /*Return new Tria*/
    2275         return tria;
    22762309}
    22772310/*}}}*/
    … …  
    23062339}
    23072340/*}}}*/
     2341Tria*      Penta::SpawnTria(int index1,int index2,int index3){/*{{{*/
     2342
     2343        int analysis_counter;
     2344
     2345        /*go into parameters and get the analysis_counter: */
     2346        this->parameters->FindParam(&analysis_counter,AnalysisCounterEnum);
     2347
     2348        /*Create Tria*/
     2349        Tria* tria=new Tria();
     2350        tria->id=this->id;
     2351        tria->inputs=(Inputs*)this->inputs->SpawnTriaInputs(index1,index2,index3);
     2352        tria->parameters=this->parameters;
     2353        tria->element_type=P1Enum; //Only P1 CG for now (TO BE CHANGED)
     2354        this->SpawnTriaHook(xDynamicCast<ElementHook*>(tria),index1,index2,index3);
     2355
     2356        /*Spawn material*/
     2357        tria->material=(Material*)this->material->copy2(tria);
     2358
     2359        /*recover nodes, material and matpar: */
     2360        tria->nodes=(Node**)tria->hnodes[analysis_counter]->deliverp();
     2361        tria->vertices=(Vertex**)tria->hvertices->deliverp();
     2362        tria->matpar=(Matpar*)tria->hmatpar->delivers();
     2363
     2364        /*Return new Tria*/
     2365        return tria;
     2366}
     2367/*}}}*/
     2368IssmDouble Penta::StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa){/*{{{*/
     2369        /*Compute stabilization parameter*/
     2370        /*kappa=thermalconductivity/(rho_ice*hearcapacity) for thermal model*/
     2371        /*kappa=enthalpydiffusionparameter for enthalpy model*/
     2372
     2373        IssmDouble normu;
     2374        IssmDouble tau_parameter;
     2375
     2376        normu=pow(pow(u,2)+pow(v,2)+pow(w,2),0.5);
     2377        if(normu*diameter/(3*2*kappa)<1){
     2378                tau_parameter=pow(diameter,2)/(3*2*2*kappa);
     2379        }
     2380        else tau_parameter=diameter/(2*normu);
     2381
     2382        return tau_parameter;
     2383}
     2384/*}}}*/
     2385void       Penta::StrainRateparallel(){/*{{{*/
     2386
     2387        IssmDouble *xyz_list = NULL;
     2388        IssmDouble  epsilon[6];
     2389        GaussPenta* gauss=NULL;
     2390        IssmDouble  vx,vy,vel;
     2391        IssmDouble  strainxx;
     2392        IssmDouble  strainxy;
     2393        IssmDouble  strainyy;
     2394        IssmDouble  strainparallel[NUMVERTICES];
     2395
     2396        /* Get node coordinates and dof list: */
     2397        this->GetVerticesCoordinates(&xyz_list);
     2398
     2399        /*Retrieve all inputs we will need*/
     2400        Input* vx_input=inputs->GetInput(VxEnum);                                  _assert_(vx_input);
     2401        Input* vy_input=inputs->GetInput(VyEnum);                                  _assert_(vy_input);
     2402        Input* vz_input=inputs->GetInput(VzEnum);                                                                                               _assert_(vz_input);
     2403
     2404        /* Start looping on the number of vertices: */
     2405        gauss=new GaussPenta();
     2406        for (int iv=0;iv<NUMVERTICES;iv++){
     2407                gauss->GaussVertex(iv);
     2408
     2409                /* Get the value we need*/
     2410                vx_input->GetInputValue(&vx,gauss);
     2411                vy_input->GetInputValue(&vy,gauss);
     2412                vel=vx*vx+vy*vy;
     2413
     2414                /*Compute strain rate and viscosity: */
     2415                this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input);
     2416                strainxx=epsilon[0];
     2417                strainyy=epsilon[1];
     2418                strainxy=epsilon[3];
     2419
     2420                /*strainparallel= Strain rate along the ice flow direction */
     2421                strainparallel[iv]=(vx*vx*(strainxx)+vy*vy*(strainyy)+2*vy*vx*strainxy)/(vel+1.e-6);
     2422        }
     2423
     2424        /*Add input*/
     2425        this->inputs->AddInput(new PentaInput(StrainRateparallelEnum,&strainparallel[0],P1Enum));
     2426
     2427        /*Clean up and return*/
     2428        delete gauss;
     2429        xDelete<IssmDouble>(xyz_list);
     2430}
     2431/*}}}*/
     2432void       Penta::StrainRateperpendicular(){/*{{{*/
     2433
     2434        IssmDouble *xyz_list = NULL;
     2435        IssmDouble  epsilon[6];
     2436        GaussPenta* gauss=NULL;
     2437        IssmDouble  vx,vy,vel;
     2438        IssmDouble  strainxx;
     2439        IssmDouble  strainxy;
     2440        IssmDouble  strainyy;
     2441        IssmDouble  strainperpendicular[NUMVERTICES];
     2442
     2443        /* Get node coordinates and dof list: */
     2444        this->GetVerticesCoordinates(&xyz_list);
     2445
     2446        /*Retrieve all inputs we will need*/
     2447        Input* vx_input=inputs->GetInput(VxEnum);                                  _assert_(vx_input);
     2448        Input* vy_input=inputs->GetInput(VyEnum);                                  _assert_(vy_input);
     2449        Input* vz_input=inputs->GetInput(VzEnum);                                                                                               _assert_(vz_input);
     2450
     2451        /* Start looping on the number of vertices: */
     2452        gauss=new GaussPenta();
     2453        for (int iv=0;iv<NUMVERTICES;iv++){
     2454                gauss->GaussVertex(iv);
     2455
     2456                /* Get the value we need*/
     2457                vx_input->GetInputValue(&vx,gauss);
     2458                vy_input->GetInputValue(&vy,gauss);
     2459                vel=vx*vx+vy*vy;
     2460
     2461                /*Compute strain rate and viscosity: */
     2462                this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input);
     2463                strainxx=epsilon[0];
     2464                strainyy=epsilon[1];
     2465                strainxy=epsilon[3];
     2466
     2467                /*strainperpendicular= Strain rate perpendicular to the ice flow direction */
     2468                strainperpendicular[iv]=(vx*vx*(strainyy)+vy*vy*(strainxx)-2*vy*vx*strainxy)/(vel+1.e-6);
     2469        }
     2470
     2471        /*Add input*/
     2472        this->inputs->AddInput(new PentaInput(StrainRateperpendicularEnum,&strainperpendicular[0],P1Enum));
     2473
     2474        /*Clean up and return*/
     2475        delete gauss;
     2476        xDelete<IssmDouble>(xyz_list);
     2477}
     2478/*}}}*/
     2479void       Penta::StressIntensityFactor(){/*{{{*/
     2480
     2481        /* Check if we are on the base */
     2482        if(!IsOnBase()) return;
     2483
     2484        IssmDouble  ki[6]={0.};
     2485        IssmDouble  const_grav=9.81;
     2486        IssmDouble  rho_ice=900;
     2487        IssmDouble  rho_water=1000;
     2488        IssmDouble  Jdet[3];
     2489        IssmDouble  pressure,vx,vy,vel,deviaxx,deviaxy,deviayy,water_depth,prof,stress_xx,thickness;
     2490
     2491        Penta* penta=this;
     2492        for(;;){
     2493       
     2494                IssmDouble  xyz_list[NUMVERTICES][3];
     2495                /* Get node coordinates and dof list: */
     2496                ::GetVerticesCoordinates(&xyz_list[0][0],penta->vertices,NUMVERTICES);
     2497
     2498                ///*Compute the Jacobian for the vertical integration*/
     2499                Jdet[0]=(xyz_list[3][2]-xyz_list[0][2])*0.5;
     2500                Jdet[1]=(xyz_list[4][2]-xyz_list[1][2])*0.5;
     2501                Jdet[2]=(xyz_list[5][2]-xyz_list[2][2])*0.5;
     2502       
     2503                /*Retrieve all inputs we will need*/
     2504                Input* vx_input=inputs->GetInput(VxEnum);                                  _assert_(vx_input);
     2505                Input* vy_input=inputs->GetInput(VyEnum);                                  _assert_(vy_input);
     2506                Input* vel_input=inputs->GetInput(VelEnum);                                _assert_(vel_input);
     2507                Input* pressure_input=inputs->GetInput(PressureEnum);                      _assert_(pressure_input);
     2508                Input* deviaxx_input=inputs->GetInput(DeviatoricStressxxEnum);             _assert_(deviaxx_input);
     2509                Input* deviaxy_input=inputs->GetInput(DeviatoricStressxyEnum);             _assert_(deviaxy_input);
     2510                Input* deviayy_input=inputs->GetInput(DeviatoricStressyyEnum);             _assert_(deviayy_input);
     2511                Input* surface_input=inputs->GetInput(SurfaceEnum);                                                             _assert_(surface_input);
     2512                Input* thickness_input=inputs->GetInput(ThicknessEnum);                                                 _assert_(thickness_input);
     2513               
     2514                /* Start looping on the number of 2D vertices: */
     2515                for(int ig=0;ig<3;ig++){
     2516                        GaussPenta* gauss=new GaussPenta(ig,3+ig,11);
     2517                        for (int iv=gauss->begin();iv<gauss->end();iv++){
     2518                                gauss->GaussPoint(iv);
     2519
     2520                                /* Get the value we need*/
     2521                                pressure_input->GetInputValue(&pressure,gauss);
     2522                                vx_input->GetInputValue(&vx,gauss);
     2523                                vy_input->GetInputValue(&vy,gauss);
     2524                                vel_input->GetInputValue(&vel,gauss);
     2525                                deviaxx_input->GetInputValue(&deviaxx,gauss);
     2526                                deviaxy_input->GetInputValue(&deviaxy,gauss);
     2527                                deviayy_input->GetInputValue(&deviayy,gauss);
     2528                                surface_input->GetInputValue(&water_depth,gauss);
     2529                                thickness_input->GetInputValue(&thickness,gauss);
     2530                                prof=water_depth-penta->GetZcoord(&xyz_list[0][0],gauss);
     2531
     2532                                /*stress_xx= Deviatoric stress along the ice flow direction plus cryostatic pressure */
     2533                                stress_xx=(vx*vx*(deviaxx)+vy*vy*(deviayy)+2*vy*vx*deviaxy)/(vel*vel+1.e-6);
     2534
     2535                                if(prof<water_depth&prof<thickness){
     2536                                        /* Compute the local stress intensity factor*/
     2537                                        ki[ig]+=Jdet[ig]*gauss->weight*stress_xx*StressIntensityIntegralWeight(prof,min(water_depth,thickness),thickness);
     2538                                }
     2539                        }
     2540                        delete gauss;
     2541                }
     2542                       
     2543                /*Stop if we have reached the surface/base*/
     2544                if(penta->IsOnSurface()) break;
     2545               
     2546                /*get upper Penta*/
     2547                penta=penta->GetUpperPenta();
     2548                _assert_(penta->Id()!=this->id);
     2549        }
     2550
     2551        /*Add input*/
     2552        this->inputs->AddInput(new PentaInput(StressIntensityFactorEnum,&ki[0],P1Enum));
     2553        this->InputExtrude(StressIntensityFactorEnum,-1);
     2554}
     2555/*}}}*/
    23082556IssmDouble Penta::SurfaceArea(void){/*{{{*/
    23092557
    … …  
    23852633        return dt;
    23862634}/*}}}*/
     2635IssmDouble Penta::TotalSmb(void){/*{{{*/
     2636
     2637        /*The smb[Gt yr-1] of one element is area[m2] * smb [ m ice yr^-1] * rho_ice [kg m-3] / 1e+10^12 */
     2638        IssmDouble base,smb,rho_ice;
     2639        IssmDouble Total_Smb=0;
     2640        IssmDouble xyz_list[NUMVERTICES][3];
     2641
     2642        /*Get material parameters :*/
     2643        rho_ice=matpar->GetRhoIce();
     2644
     2645        if(!IsIceInElement() || !IsOnSurface()) return 0.;
     2646
     2647        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     2648
     2649        /*First calculate the area of the base (cross section triangle)
     2650         * http://en.wikipedia.org/wiki/Triangle
     2651         * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
     2652        base = 1./2. * fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1]));
     2653
     2654        /*Now get the average SMB over the element*/
     2655        Input* smb_input = inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(smb_input);
     2656
     2657        smb_input->GetInputAverage(&smb);
     2658        Total_Smb=rho_ice*base*smb;// smb on element in kg s-1
     2659
     2660        /*Return: */
     2661        return Total_Smb;
     2662}
     2663/*}}}*/
    23872664void       Penta::Update(int index,IoModel* iomodel,int analysis_counter,int analysis_type,int finiteelement_type){ /*{{{*/
    23882665
    … …  
    28333110}
    28343111/*}}}*/
     3112int        Penta::UpdatePotentialUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf){/*{{{*/
     3113
     3114        int i;
     3115        int nflipped=0;
     3116
     3117        /*Go through nodes, and whoever is on the potential_ungrounding, ends up in nodes_on_iceshelf: */
     3118        for(i=0;i<NUMVERTICES;i++){
     3119                if (reCast<bool,IssmDouble>(vertices_potentially_ungrounding[vertices[i]->Pid()])){
     3120                        vec_nodes_on_iceshelf->SetValue(vertices[i]->Pid(),-1.,INS_VAL);
     3121
     3122                        /*If node was not on ice shelf, we flipped*/
     3123                        if(nodes_on_iceshelf[vertices[i]->Pid()]>=0.){
     3124                                nflipped++;
     3125                        }
     3126                }
     3127        }
     3128        return nflipped;
     3129}
     3130/*}}}*/
     3131void       Penta::ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     3132        PentaRef::GetInputDerivativeValue(dvalue,values,xyz_list,gauss,P1Enum);
     3133}
     3134/*}}}*/
    28353135void       Penta::ValueP1OnGauss(IssmDouble* pvalue,IssmDouble* values,Gauss* gauss){/*{{{*/
    28363136        PentaRef::GetInputValue(pvalue,values,gauss,P1Enum);
    28373137}
    28383138/*}}}*/
    2839 void       Penta::ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    2840         PentaRef::GetInputDerivativeValue(dvalue,values,xyz_list,gauss,P1Enum);
     3139int        Penta::VelocityInterpolation(void){/*{{{*/
     3140        return PentaRef::VelocityInterpolation(this->element_type);
    28413141}
    28423142/*}}}*/
    … …  
    28753175}
    28763176/*}}}*/
    2877 int        Penta::VelocityInterpolation(void){/*{{{*/
    2878         return PentaRef::VelocityInterpolation(this->element_type);
    2879 }
    2880 /*}}}*/
    2881 int        Penta::PressureInterpolation(void){/*{{{*/
    2882         return PentaRef::PressureInterpolation(this->element_type);
    2883 }
    2884 /*}}}*/
    2885 bool       Penta::IsZeroLevelset(int levelset_enum){/*{{{*/
    2886 
    2887         bool        iszerols;
    2888         IssmDouble  ls[NUMVERTICES];
    2889 
    2890         /*Retrieve all inputs and parameters*/
    2891         GetInputListOnVertices(&ls[0],levelset_enum);
    2892 
    2893         /*If the level set has always same sign, there is no ice front here*/
    2894         iszerols = false;
    2895         if(IsIceInElement()){
    2896                 if(ls[0]*ls[1]<0. || ls[0]*ls[2]<0. || (ls[0]*ls[1]*ls[2]==0. && ls[0]*ls[1]+ls[0]*ls[2]+ls[1]*ls[2]<=0.)){
    2897                         iszerols = true;
    2898                 }
    2899         }
    2900         return iszerols;
    2901 }
    2902 /*}}}*/
    2903 bool       Penta::IsIcefront(void){/*{{{*/
    2904 
    2905         bool isicefront;
    2906         int i,nrice;
    2907    IssmDouble ls[NUMVERTICES];
    2908 
    2909         /*Retrieve all inputs and parameters*/
    2910         GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum);
    2911 
    2912         /* If only one vertex has ice, there is an ice front here */
    2913         isicefront=false;
    2914         if(IsIceInElement()){
    2915                 nrice=0;       
    2916                 for(i=0;i<NUMVERTICES2D;i++)
    2917                         if(ls[i]<0.) nrice++;
    2918                 if(nrice==1) isicefront= true;
    2919         }
    2920         return isicefront;
    2921 }/*}}}*/
    2922 void       Penta::AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part){/*{{{*/
    2923         _error_("Not supported yet!");
    2924 }
    2925 /*}}}*/
    2926 IssmDouble Penta::IceVolume(void){/*{{{*/
    2927 
    2928         /*The volume of a troncated prism is base * 1/3 sum(length of edges)*/
    2929         IssmDouble base,height;
    2930         IssmDouble xyz_list[NUMVERTICES][3];
    2931 
    2932         if(!IsIceInElement())return 0;
    2933 
    2934         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    2935 
    2936         /*First calculate the area of the base (cross section triangle)
    2937          * http://en.wikipedia.org/wiki/Pentangle
    2938          * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
    2939         base = 1./2.*fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1]));
    2940 
    2941         /*Now get the average height*/
    2942         height = 1./3.*((xyz_list[3][2]-xyz_list[0][2])+(xyz_list[4][2]-xyz_list[1][2])+(xyz_list[5][2]-xyz_list[2][2]));
    2943 
    2944         /*Return: */
    2945         return base*height;
    2946 }
    2947 /*}}}*/
    2948 IssmDouble Penta::IceVolumeAboveFloatation(void){/*{{{*/
    2949 
    2950         /*Volume above floatation: H + rho_water/rho_ice*bathymetry for nodes on the bed*/
    2951         IssmDouble rho_ice,rho_water;
    2952         IssmDouble base,bed,surface,bathymetry;
    2953         IssmDouble xyz_list[NUMVERTICES][3];
    2954 
    2955         if(!IsIceInElement() || IsFloating() || !IsOnBase())return 0;
    2956 
    2957         rho_ice=matpar->GetRhoIce();
    2958         rho_water=matpar->GetRhoWater();
    2959         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    2960 
    2961         /*First calculate the area of the base (cross section triangle)
    2962          * http://en.wikipedia.org/wiki/Pentangle
    2963          * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
    2964         base = 1./2.*fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1]));
    2965 
    2966         /*Now get the average height above floatation*/
    2967         Input* surface_input    = inputs->GetInput(SurfaceEnum);    _assert_(surface_input);
    2968         Input* base_input        = inputs->GetInput(BaseEnum);        _assert_(base_input);
    2969         Input* bed_input = inputs->GetInput(BedEnum); _assert_(bed_input);
    2970         surface_input->GetInputAverage(&surface);
    2971         base_input->GetInputAverage(&bed);
    2972         bed_input->GetInputAverage(&bathymetry);
    2973 
    2974         /*Return: */
    2975         return base*(surface - bed + min( rho_water/rho_ice * bathymetry, 0.) );
    2976 }
    2977 /*}}}*/
    2978 IssmDouble Penta::MassFlux( IssmDouble* segment){/*{{{*/
    2979 
    2980         IssmDouble mass_flux=0;
    2981 
    2982         if(!IsOnBase()) return mass_flux;
    2983 
    2984         /*Depth Averaging Vx and Vy*/
    2985         this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
    2986         this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
    2987 
    2988         /*Spawn Tria element from the base of the Penta: */
    2989         Tria* tria=(Tria*)SpawnTria(0,1,2);
    2990         mass_flux=tria->MassFlux(segment);
    2991         delete tria->material; delete tria;
    2992 
    2993         /*Delete Vx and Vy averaged*/
    2994         this->inputs->DeleteInput(VxAverageEnum);
    2995         this->inputs->DeleteInput(VyAverageEnum);
    2996 
    2997         /*clean up and return*/
    2998         return mass_flux;
    2999 }
    3000 /*}}}*/
    3001 IssmDouble Penta::MassFlux( IssmDouble x1, IssmDouble y1, IssmDouble x2, IssmDouble y2,int segment_id){/*{{{*/
    3002 
    3003         IssmDouble mass_flux=0;
    3004 
    3005         if(!IsOnBase()) return mass_flux;
    3006 
    3007         /*Depth Averaging Vx and Vy*/
    3008         this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
    3009         this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
    3010 
    3011         /*Spawn Tria element from the base of the Penta: */
    3012         Tria* tria=(Tria*)SpawnTria(0,1,2);
    3013         mass_flux=tria->MassFlux(x1,y1,x2,y2,segment_id);
    3014         delete tria->material; delete tria;
    3015 
    3016         /*Delete Vx and Vy averaged*/
    3017         this->inputs->DeleteInput(VxAverageEnum);
    3018         this->inputs->DeleteInput(VyAverageEnum);
    3019 
    3020         /*clean up and return*/
    3021         return mass_flux;
    3022 }
    3023 /*}}}*/
    3024 void       Penta::ElementResponse(IssmDouble* presponse,int response_enum){/*{{{*/
    3025 
    3026         switch(response_enum){
    3027                 case MaterialsRheologyBbarEnum:
    3028                         *presponse=this->material->GetBbar();
    3029                         break;
    3030                 case DamageDbarEnum:
    3031                         *presponse=this->material->GetDbar();
    3032                         break;
    3033                 case VelEnum:
    3034                         {
    3035 
    3036                                 /*Get input:*/
    3037                                 IssmDouble vel;
    3038                                 Input* vel_input;
    3039 
    3040                                 vel_input=this->inputs->GetInput(VelEnum); _assert_(vel_input);
    3041                                 vel_input->GetInputAverage(&vel);
    3042 
    3043                                 /*Assign output pointers:*/
    3044                                 *presponse=vel;
    3045                         }
    3046                         break;
    3047                 default: 
    3048                         _error_("Response type " << EnumToStringx(response_enum) << " not supported yet!");
    3049         }
    3050 
    3051 }
    3052 /*}}}*/
    3053 IssmDouble Penta::TotalSmb(void){/*{{{*/
    3054 
    3055         /*The smb[Gt yr-1] of one element is area[m2] * smb [ m ice yr^-1] * rho_ice [kg m-3] / 1e+10^12 */
    3056         IssmDouble base,smb,rho_ice;
    3057         IssmDouble Total_Smb=0;
    3058         IssmDouble xyz_list[NUMVERTICES][3];
    3059 
    3060         /*Get material parameters :*/
    3061         rho_ice=matpar->GetRhoIce();
    3062 
    3063         if(!IsIceInElement() || !IsOnSurface()) return 0.;
    3064 
    3065         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    3066 
    3067         /*First calculate the area of the base (cross section triangle)
    3068          * http://en.wikipedia.org/wiki/Triangle
    3069          * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
    3070         base = 1./2. * fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1]));
    3071 
    3072         /*Now get the average SMB over the element*/
    3073         Input* smb_input = inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(smb_input);
    3074 
    3075         smb_input->GetInputAverage(&smb);
    3076         Total_Smb=rho_ice*base*smb;// smb on element in kg s-1
    3077 
    3078         /*Return: */
    3079         return Total_Smb;
     3177void       Penta::ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum){/*{{{*/
     3178        /*Compute portion of the element that is grounded*/
     3179
     3180        int         normal_orientation=0;
     3181        IssmDouble  s1,s2;
     3182        IssmDouble  levelset[NUMVERTICES];
     3183        IssmDouble* xyz_zero = xNew<IssmDouble>(4*3);
     3184
     3185        /*Recover parameters and values*/
     3186        GetInputListOnVertices(&levelset[0],levelsetenum);
     3187
     3188        if(levelset[0]*levelset[1]>0.){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
     3189                /*Portion of the segments*/
     3190                s1=levelset[2]/(levelset[2]-levelset[1]);
     3191                s2=levelset[2]/(levelset[2]-levelset[0]);
     3192
     3193                if(levelset[2]<0.) normal_orientation=1; //orientation of quadrangle at base and top, depending on distribution of levelsetfunction
     3194                /*New point 1*/
     3195                xyz_zero[3*normal_orientation+0]=xyz_list[2*3+0]+s1*(xyz_list[1*3+0]-xyz_list[2*3+0]);
     3196                xyz_zero[3*normal_orientation+1]=xyz_list[2*3+1]+s1*(xyz_list[1*3+1]-xyz_list[2*3+1]);
     3197                xyz_zero[3*normal_orientation+2]=xyz_list[2*3+2]+s1*(xyz_list[1*3+2]-xyz_list[2*3+2]);
     3198
     3199                /*New point 0*/
     3200                xyz_zero[3*(1-normal_orientation)+0]=xyz_list[2*3+0]+s2*(xyz_list[0*3+0]-xyz_list[2*3+0]);
     3201                xyz_zero[3*(1-normal_orientation)+1]=xyz_list[2*3+1]+s2*(xyz_list[0*3+1]-xyz_list[2*3+1]);
     3202                xyz_zero[3*(1-normal_orientation)+2]=xyz_list[2*3+2]+s2*(xyz_list[0*3+2]-xyz_list[2*3+2]);
     3203
     3204                /*New point 3*/
     3205                xyz_zero[3*(2+1-normal_orientation)+0]=xyz_list[5*3+0]+s1*(xyz_list[4*3+0]-xyz_list[5*3+0]);
     3206                xyz_zero[3*(2+1-normal_orientation)+1]=xyz_list[5*3+1]+s1*(xyz_list[4*3+1]-xyz_list[5*3+1]);
     3207                xyz_zero[3*(2+1-normal_orientation)+2]=xyz_list[5*3+2]+s1*(xyz_list[4*3+2]-xyz_list[5*3+2]);
     3208
     3209                /*New point 4*/
     3210                xyz_zero[3*(2+normal_orientation)+0]=xyz_list[5*3+0]+s2*(xyz_list[3*3+0]-xyz_list[5*3+0]);
     3211                xyz_zero[3*(2+normal_orientation)+1]=xyz_list[5*3+1]+s2*(xyz_list[3*3+1]-xyz_list[5*3+1]);
     3212                xyz_zero[3*(2+normal_orientation)+2]=xyz_list[5*3+2]+s2*(xyz_list[3*3+2]-xyz_list[5*3+2]);
     3213        }
     3214        else if(levelset[1]*levelset[2]>0.){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
     3215                /*Portion of the segments*/
     3216                s1=levelset[0]/(levelset[0]-levelset[2]);
     3217                s2=levelset[0]/(levelset[0]-levelset[1]);
     3218
     3219                if(levelset[0]<0.) normal_orientation=1;
     3220                /*New point 1*/
     3221                xyz_zero[3*normal_orientation+0]=xyz_list[0*3+0]+s1*(xyz_list[2*3+0]-xyz_list[0*3+0]);
     3222                xyz_zero[3*normal_orientation+1]=xyz_list[0*3+1]+s1*(xyz_list[2*3+1]-xyz_list[0*3+1]);
     3223                xyz_zero[3*normal_orientation+2]=xyz_list[0*3+2]+s1*(xyz_list[2*3+2]-xyz_list[0*3+2]);
     3224
     3225                /*New point 2*/
     3226                xyz_zero[3*(1-normal_orientation)+0]=xyz_list[0*3+0]+s2*(xyz_list[1*3+0]-xyz_list[0*3+0]);
     3227                xyz_zero[3*(1-normal_orientation)+1]=xyz_list[0*3+1]+s2*(xyz_list[1*3+1]-xyz_list[0*3+1]);
     3228                xyz_zero[3*(1-normal_orientation)+2]=xyz_list[0*3+2]+s2*(xyz_list[1*3+2]-xyz_list[0*3+2]);
     3229
     3230                /*New point 3*/
     3231                xyz_zero[3*(2+1-normal_orientation)+0]=xyz_list[3*3+0]+s1*(xyz_list[5*3+0]-xyz_list[3*3+0]);
     3232                xyz_zero[3*(2+1-normal_orientation)+1]=xyz_list[3*3+1]+s1*(xyz_list[5*3+1]-xyz_list[3*3+1]);
     3233                xyz_zero[3*(2+1-normal_orientation)+2]=xyz_list[3*3+2]+s1*(xyz_list[5*3+2]-xyz_list[3*3+2]);
     3234
     3235                /*New point 4*/
     3236                xyz_zero[3*(2+normal_orientation)+0]=xyz_list[3*3+0]+s2*(xyz_list[4*3+0]-xyz_list[3*3+0]);
     3237                xyz_zero[3*(2+normal_orientation)+1]=xyz_list[3*3+1]+s2*(xyz_list[4*3+1]-xyz_list[3*3+1]);
     3238                xyz_zero[3*(2+normal_orientation)+2]=xyz_list[3*3+2]+s2*(xyz_list[4*3+2]-xyz_list[3*3+2]);
     3239        }
     3240        else if(levelset[0]*levelset[2]>0.){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
     3241                /*Portion of the segments*/
     3242                s1=levelset[1]/(levelset[1]-levelset[0]);
     3243                s2=levelset[1]/(levelset[1]-levelset[2]);
     3244
     3245                if(levelset[1]<0.) normal_orientation=1;
     3246                /*New point 0*/
     3247                xyz_zero[3*normal_orientation+0]=xyz_list[1*3+0]+s1*(xyz_list[0*3+0]-xyz_list[1*3+0]);
     3248                xyz_zero[3*normal_orientation+1]=xyz_list[1*3+1]+s1*(xyz_list[0*3+1]-xyz_list[1*3+1]);
     3249                xyz_zero[3*normal_orientation+2]=xyz_list[1*3+2]+s1*(xyz_list[0*3+2]-xyz_list[1*3+2]);
     3250
     3251                /*New point 2*/
     3252                xyz_zero[3*(1-normal_orientation)+0]=xyz_list[1*3+0]+s2*(xyz_list[2*3+0]-xyz_list[1*3+0]);
     3253                xyz_zero[3*(1-normal_orientation)+1]=xyz_list[1*3+1]+s2*(xyz_list[2*3+1]-xyz_list[1*3+1]);
     3254                xyz_zero[3*(1-normal_orientation)+2]=xyz_list[1*3+2]+s2*(xyz_list[2*3+2]-xyz_list[1*3+2]);
     3255
     3256                /*New point 3*/
     3257                xyz_zero[3*(2+1-normal_orientation)+0]=xyz_list[4*3+0]+s1*(xyz_list[3*3+0]-xyz_list[4*3+0]);
     3258                xyz_zero[3*(2+1-normal_orientation)+1]=xyz_list[4*3+1]+s1*(xyz_list[3*3+1]-xyz_list[4*3+1]);
     3259                xyz_zero[3*(2+1-normal_orientation)+2]=xyz_list[4*3+2]+s1*(xyz_list[3*3+2]-xyz_list[4*3+2]);
     3260
     3261                /*New point 4*/
     3262                xyz_zero[3*(2+normal_orientation)+0]=xyz_list[4*3+0]+s2*(xyz_list[5*3+0]-xyz_list[4*3+0]);
     3263                xyz_zero[3*(2+normal_orientation)+1]=xyz_list[4*3+1]+s2*(xyz_list[5*3+1]-xyz_list[4*3+1]);
     3264                xyz_zero[3*(2+normal_orientation)+2]=xyz_list[4*3+2]+s2*(xyz_list[5*3+2]-xyz_list[4*3+2]);
     3265        }
     3266        else if(levelset[0]==0. && levelset[1]==0.){ //front is on point 0 and 1
     3267                xyz_zero[3*0+0]=xyz_list[0*3+0];
     3268                xyz_zero[3*0+1]=xyz_list[0*3+1];
     3269                xyz_zero[3*0+2]=xyz_list[0*3+2];
     3270
     3271                /*New point 2*/
     3272                xyz_zero[3*1+0]=xyz_list[1*3+0];
     3273                xyz_zero[3*1+1]=xyz_list[1*3+1];
     3274                xyz_zero[3*1+2]=xyz_list[1*3+2];
     3275
     3276                /*New point 3*/
     3277                xyz_zero[3*2+0]=xyz_list[4*3+0];
     3278                xyz_zero[3*2+1]=xyz_list[4*3+1];
     3279                xyz_zero[3*2+2]=xyz_list[4*3+2];
     3280
     3281                /*New point 4*/
     3282                xyz_zero[3*3+0]=xyz_list[3*3+0];
     3283                xyz_zero[3*3+1]=xyz_list[3*3+1];
     3284                xyz_zero[3*3+2]=xyz_list[3*3+2];
     3285        }
     3286        else if(levelset[0]==0. && levelset[2]==0.){ //front is on point 0 and 1
     3287                xyz_zero[3*0+0]=xyz_list[2*3+0];
     3288                xyz_zero[3*0+1]=xyz_list[2*3+1];
     3289                xyz_zero[3*0+2]=xyz_list[2*3+2];
     3290
     3291                /*New point 2*/
     3292                xyz_zero[3*1+0]=xyz_list[0*3+0];
     3293                xyz_zero[3*1+1]=xyz_list[0*3+1];
     3294                xyz_zero[3*1+2]=xyz_list[0*3+2];
     3295
     3296                /*New point 3*/
     3297                xyz_zero[3*2+0]=xyz_list[3*3+0];
     3298                xyz_zero[3*2+1]=xyz_list[3*3+1];
     3299                xyz_zero[3*2+2]=xyz_list[3*3+2];
     3300
     3301                /*New point 4*/
     3302                xyz_zero[3*3+0]=xyz_list[5*3+0];
     3303                xyz_zero[3*3+1]=xyz_list[5*3+1];
     3304                xyz_zero[3*3+2]=xyz_list[5*3+2];
     3305        }
     3306        else if(levelset[1]==0. && levelset[2]==0.){ //front is on point 0 and 1
     3307                xyz_zero[3*0+0]=xyz_list[1*3+0];
     3308                xyz_zero[3*0+1]=xyz_list[1*3+1];
     3309                xyz_zero[3*0+2]=xyz_list[1*3+2];
     3310
     3311                /*New point 2*/
     3312                xyz_zero[3*1+0]=xyz_list[2*3+0];
     3313                xyz_zero[3*1+1]=xyz_list[2*3+1];
     3314                xyz_zero[3*1+2]=xyz_list[2*3+2];
     3315
     3316                /*New point 3*/
     3317                xyz_zero[3*2+0]=xyz_list[5*3+0];
     3318                xyz_zero[3*2+1]=xyz_list[5*3+1];
     3319                xyz_zero[3*2+2]=xyz_list[5*3+2];
     3320
     3321                /*New point 4*/
     3322                xyz_zero[3*3+0]=xyz_list[4*3+0];
     3323                xyz_zero[3*3+1]=xyz_list[4*3+1];
     3324                xyz_zero[3*3+2]=xyz_list[4*3+2];
     3325        }
     3326        else _error_("Case not covered");
     3327
     3328        /*Assign output pointer*/
     3329        *pxyz_zero= xyz_zero;
    30803330}
    30813331/*}}}*/
    … …  
    30873337/*}}}*/
    30883338#endif
    3089 
    3090 void       Penta::ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){/*{{{*/
    3091 
    3092         int    vertexpidlist[NUMVERTICES];
    3093         IssmDouble grad_list[NUMVERTICES];
    3094         Input* grad_input=NULL;
    3095         Input* input=NULL;
    3096 
    3097         if(enum_type==MaterialsRheologyBbarEnum){
    3098                 input=(Input*)inputs->GetInput(MaterialsRheologyBEnum);
    3099         }
    3100         else if(enum_type==DamageDbarEnum){
    3101                 input=(Input*)inputs->GetInput(DamageDEnum);
    3102         }
    3103         else{
    3104                 input=inputs->GetInput(enum_type);
    3105         }
    3106         if (!input) _error_("Input " << EnumToStringx(enum_type) << " not found");
    3107         if (input->ObjectEnum()!=ControlInputEnum) _error_("Input " << EnumToStringx(enum_type) << " is not a ControlInput");
    3108 
    3109         GradientIndexing(&vertexpidlist[0],control_index);
    3110         for(int i=0;i<NUMVERTICES;i++) grad_list[i]=gradient[vertexpidlist[i]];
    3111         grad_input=new PentaInput(GradientEnum,grad_list,P1Enum);
    3112         ((ControlInput*)input)->SetGradient(grad_input);
    3113 
    3114 }/*}}}*/
    3115 void       Penta::ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum){/*{{{*/
    3116 
    3117         Input* input=NULL;
    3118 
    3119         if(control_enum==MaterialsRheologyBbarEnum){
    3120                 input=(Input*)inputs->GetInput(MaterialsRheologyBEnum);
    3121         }
    3122         else if(control_enum==DamageDbarEnum){
    3123                 input=(Input*)inputs->GetInput(DamageDEnum);
    3124         }
    3125         else{
    3126                 input=inputs->GetInput(control_enum);
    3127         }
    3128         if (!input) _error_("Input " << EnumToStringx(control_enum) << " not found");
    3129         if (input->ObjectEnum()!=ControlInputEnum) _error_("Input " << EnumToStringx(control_enum) << " is not a ControlInput");
    3130 
    3131         int         sidlist[NUMVERTICES];
    3132         int         connectivity[NUMVERTICES];
    3133         IssmPDouble values[NUMVERTICES];
    3134         IssmPDouble gradients[NUMVERTICES];
    3135         IssmDouble  value,gradient;
    3136 
    3137         this->GetVerticesConnectivityList(&connectivity[0]);
    3138         this->GetVerticesSidList(&sidlist[0]);
    3139 
    3140         GaussPenta* gauss=new GaussPenta();
    3141         for (int iv=0;iv<NUMVERTICES;iv++){
    3142                 gauss->GaussVertex(iv);
    3143 
    3144                 ((ControlInput*)input)->GetInputValue(&value,gauss);
    3145                 ((ControlInput*)input)->GetGradientValue(&gradient,gauss);
    3146 
    3147                 values[iv]    = reCast<IssmPDouble>(value)/reCast<IssmPDouble>(connectivity[iv]);
    3148                 gradients[iv] = reCast<IssmPDouble>(gradient)/reCast<IssmPDouble>(connectivity[iv]);
    3149         }
    3150         delete gauss;
    3151 
    3152         vector_control->SetValues(NUMVERTICES,&sidlist[0],&values[0],ADD_VAL);
    3153         vector_gradient->SetValues(NUMVERTICES,&sidlist[0],&gradients[0],ADD_VAL);
    3154 
    3155 }/*}}}*/
    3156 void       Penta::InputControlUpdate(IssmDouble scalar,bool save_parameter){/*{{{*/
    3157 
    3158         /*Intermediary*/
    3159         int    num_controls;
    3160         int*   control_type=NULL;
    3161         Input* input=NULL;
    3162 
    3163         /*retrieve some parameters: */
    3164         this->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
    3165         this->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
    3166 
    3167         for(int i=0;i<num_controls;i++){
    3168 
    3169                 if(control_type[i]==MaterialsRheologyBbarEnum){
    3170                         if (!IsOnBase()) goto cleanup_and_return;
    3171                         input=(Input*)this->inputs->GetInput(MaterialsRheologyBEnum); _assert_(input);
    3172                 }
    3173                 else if(control_type[i]==DamageDbarEnum){
    3174                         if (!IsOnBase()) goto cleanup_and_return;
    3175                         input=(Input*)this->inputs->GetInput(DamageDEnum); _assert_(input);
    3176                 }
    3177                 else{
    3178                         input=(Input*)this->inputs->GetInput(control_type[i]); _assert_(input);
    3179                 }
    3180                 if(input->ObjectEnum()!=ControlInputEnum) _error_("input " << EnumToStringx(control_type[i]) << " is not a ControlInput");
    3181 
    3182                 ((ControlInput*)input)->UpdateValue(scalar);
    3183                 ((ControlInput*)input)->Constrain();
    3184                 if (save_parameter) ((ControlInput*)input)->SaveValue();
    3185 
    3186                 if(control_type[i]==MaterialsRheologyBbarEnum){
    3187                         this->InputExtrude(MaterialsRheologyBEnum,-1);
    3188                 }
    3189                 else if(control_type[i]==DamageDbarEnum){
    3190                         this->InputExtrude(DamageDEnum,-1);
    3191                 }
    3192         }
    3193 
    3194         /*Clean up and return*/
    3195 cleanup_and_return:
    3196         xDelete<int>(control_type);
    3197 }
    3198 /*}}}*/
    3199 void       Penta::GetVectorFromControlInputs(Vector<IssmDouble>* vector,int control_enum,int control_index,const char* data,bool onsid){/*{{{*/
    3200 
    3201         int vertexidlist[NUMVERTICES];
    3202 
    3203         /*Get out if this is not an element input*/
    3204         if(!IsInput(control_enum)) return;
    3205 
    3206         /*Prepare index list*/
    3207         GradientIndexing(&vertexidlist[0],control_index,onsid);
    3208 
    3209         /*Get input (either in element or material)*/
    3210         if(control_enum==MaterialsRheologyBbarEnum) control_enum=MaterialsRheologyBEnum;
    3211         Input* input=inputs->GetInput(control_enum);
    3212         if(!input) _error_("Input " << EnumToStringx(control_enum) << " not found in element");
    3213 
    3214         /*Check that it is a ControlInput*/
    3215         if (input->ObjectEnum()!=ControlInputEnum){
    3216                 _error_("input " << EnumToStringx(control_enum) << " is not a ControlInput");
    3217         }
    3218 
    3219         ((ControlInput*)input)->GetVectorFromInputs(vector,&vertexidlist[0],data);
    3220 }
    3221 /*}}}*/
    3222 void       Penta::SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index){/*{{{*/
    3223 
    3224         IssmDouble  values[NUMVERTICES];
    3225         int         vertexpidlist[NUMVERTICES],control_init;
    3226 
    3227         /*Specific case for depth averaged quantities*/
    3228         control_init=control_enum;
    3229         if(control_enum==MaterialsRheologyBbarEnum){
    3230                 control_enum=MaterialsRheologyBEnum;
    3231                 if(!IsOnBase()) return;
    3232         }
    3233         if(control_enum==DamageDbarEnum){
    3234                 control_enum=DamageDEnum;
    3235                 if(!IsOnBase()) return;
    3236         }
    3237 
    3238         /*Get out if this is not an element input*/
    3239         if(!IsInput(control_enum)) return;
    3240 
    3241         /*Prepare index list*/
    3242         GradientIndexing(&vertexpidlist[0],control_index);
    3243 
    3244         /*Get values on vertices*/
    3245         for(int i=0;i<NUMVERTICES;i++){
    3246                 values[i]=vector[vertexpidlist[i]];
    3247         }
    3248         Input* new_input = new PentaInput(control_enum,values,P1Enum);
    3249         Input* input=(Input*)this->inputs->GetInput(control_enum);   _assert_(input);
    3250         if(input->ObjectEnum()!=ControlInputEnum){
    3251                 _error_("input " << EnumToStringx(control_enum) << " is not a ControlInput");
    3252         }
    3253 
    3254         ((ControlInput*)input)->SetInput(new_input);
    3255 
    3256         if(control_init==MaterialsRheologyBbarEnum){
    3257                 this->InputExtrude(control_enum,-1);
    3258         }
    3259         if(control_init==DamageDbarEnum){
    3260                 this->InputExtrude(control_enum,-1);
    3261         }
    3262 }
    3263 /*}}}*/
    32643339
    32653340#ifdef _HAVE_DAKOTA_
    … …  
    34033478/*}}}*/
    34043479#endif
    3405 void       Penta::GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution, int enum_type){/*{{{*/
    3406 
    3407         const int    numdof=NDOF1*NUMVERTICES;
    3408 
    3409         int          i;
    3410         int*         doflist=NULL;
    3411         IssmDouble   values[numdof];
    3412         IssmDouble   enum_value;
    3413         GaussPenta   *gauss=NULL;
    3414 
    3415         /*Get dof list: */
    3416         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    3417         Input* enum_input=inputs->GetInput(enum_type); _assert_(enum_input);
    3418 
    3419         gauss=new GaussPenta();
    3420         for(i=0;i<NUMVERTICES;i++){
    3421                 /*Recover temperature*/
    3422                 gauss->GaussVertex(i);
    3423                 enum_input->GetInputValue(&enum_value,gauss);
    3424                 values[i]=enum_value;
    3425         }
    3426 
    3427         /*Add value to global vector*/
    3428         solution->SetValues(numdof,doflist,values,INS_VAL);
    3429 
    3430         /*Free ressources:*/
    3431         delete gauss;
    3432         xDelete<int>(doflist);
    3433 }
    3434 /*}}}*/
    3435 void       Penta::PotentialUngrounding(Vector<IssmDouble>* potential_ungrounding){/*{{{*/
    3436 
    3437         IssmDouble  h[NUMVERTICES],r[NUMVERTICES],gl[NUMVERTICES];
    3438         IssmDouble  bed_hydro;
    3439         IssmDouble  rho_water,rho_ice,density;
    3440 
    3441         /*material parameters: */
    3442         rho_water=matpar->GetRhoWater();
    3443         rho_ice=matpar->GetRhoIce();
    3444         density=rho_ice/rho_water;
    3445         GetInputListOnVertices(&h[0],ThicknessEnum);
    3446         GetInputListOnVertices(&r[0],BedEnum);
    3447         GetInputListOnVertices(&gl[0],MaskGroundediceLevelsetEnum);
    3448 
    3449         /*go through vertices, and figure out which ones are on the ice sheet, and want to unground: */
    3450         for(int i=0;i<NUMVERTICES;i++){
    3451                 /*Find if grounded vertices want to start floating*/
    3452                 if (gl[i]>0.){
    3453                         bed_hydro=-density*h[i];
    3454                         if(bed_hydro>r[i]){
    3455                                 /*Vertex that could potentially unground, flag it*/
    3456                                 potential_ungrounding->SetValue(vertices[i]->Pid(),1,INS_VAL);
    3457                         }
    3458                 }
    3459         }
    3460 }
    3461 /*}}}*/
    3462 int        Penta::UpdatePotentialUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf){/*{{{*/
    3463 
    3464         int i;
    3465         int nflipped=0;
    3466 
    3467         /*Go through nodes, and whoever is on the potential_ungrounding, ends up in nodes_on_iceshelf: */
    3468         for(i=0;i<NUMVERTICES;i++){
    3469                 if (reCast<bool,IssmDouble>(vertices_potentially_ungrounding[vertices[i]->Pid()])){
    3470                         vec_nodes_on_iceshelf->SetValue(vertices[i]->Pid(),-1.,INS_VAL);
    3471 
    3472                         /*If node was not on ice shelf, we flipped*/
    3473                         if(nodes_on_iceshelf[vertices[i]->Pid()]>=0.){
    3474                                 nflipped++;
    3475                         }
    3476                 }
    3477         }
    3478         return nflipped;
    3479 }
    3480 /*}}}*/
Note: See TracChangeset for help on using the changeset viewer.