source: issm/oecreview/Archive/18296-19100/ISSM-18910-18911.diff@ 19102

Last change on this file since 19102 was 19102, checked in by Mathieu Morlighem, 10 years ago

NEW: added 18296-19100
File size: 109.9 KB

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

     
    154154        this->inputs->AddInput(new TriaInput(input_enum,values,interpolation_enum));
    155155}
    156156/*}}}*/
     157void       Tria::AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part){/*{{{*/
     158
     159        bool       already = false;
     160        int        i,j;
     161        int        partition[NUMVERTICES];
     162        int        offsetsid[NUMVERTICES];
     163        int        offsetdof[NUMVERTICES];
     164        IssmDouble area;
     165        IssmDouble mean;
     166
     167        /*First, get the area: */
     168        area=this->GetArea();
     169
     170        /*Figure out the average for this element: */
     171        this->GetVerticesSidList(&offsetsid[0]);
     172        this->GetVertexPidList(&offsetdof[0]);
     173        mean=0;
     174        for(i=0;i<NUMVERTICES;i++){
     175                partition[i]=reCast<int>(qmu_part[offsetsid[i]]);
     176                mean=mean+1.0/NUMVERTICES*vertex_response[offsetdof[i]];
     177        }
     178
     179        /*Add contribution: */
     180        for(i=0;i<NUMVERTICES;i++){
     181                already=false;
     182                for(j=0;j<i;j++){
     183                        if (partition[i]==partition[j]){
     184                                already=true;
     185                                break;
     186                        }
     187                }
     188                if(!already){
     189                        partition_contributions->SetValue(partition[i],mean*area,ADD_VAL);
     190                        partition_areas->SetValue(partition[i],area,ADD_VAL);
     191                };
     192        }
     193}
     194/*}}}*/
     195void       Tria::CalvingRateLevermann(){/*{{{*/
     196
     197        IssmDouble  xyz_list[NUMVERTICES][3];
     198        GaussTria* gauss=NULL;
     199        IssmDouble  vx,vy,vel;
     200        IssmDouble  strainparallel;
     201        IssmDouble  propcoeff;
     202        IssmDouble  strainperpendicular;
     203        IssmDouble  calvingratex[NUMVERTICES];
     204        IssmDouble  calvingratey[NUMVERTICES];
     205        IssmDouble  calvingrate[NUMVERTICES];
     206
     207
     208        /* Get node coordinates and dof list: */
     209        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     210
     211        /*Retrieve all inputs and parameters we will need*/
     212        Input* vx_input=inputs->GetInput(VxEnum);                                                                                                                                               _assert_(vx_input);
     213        Input* vy_input=inputs->GetInput(VyEnum);                                                                                                                                               _assert_(vy_input);
     214        Input* strainparallel_input=inputs->GetInput(StrainRateparallelEnum);                                                           _assert_(strainparallel_input);
     215        Input* strainperpendicular_input=inputs->GetInput(StrainRateperpendicularEnum);                                 _assert_(strainperpendicular_input);
     216        Input* levermanncoeff_input=inputs->GetInput(CalvinglevermannCoeffEnum);                     _assert_(levermanncoeff_input);
     217
     218        /* Start looping on the number of vertices: */
     219        gauss=new GaussTria();
     220        for (int iv=0;iv<NUMVERTICES;iv++){
     221                gauss->GaussVertex(iv);
     222
     223                /* Get the value we need*/
     224                vx_input->GetInputValue(&vx,gauss);
     225                vy_input->GetInputValue(&vy,gauss);
     226                vel=vx*vx+vy*vy;
     227                strainparallel_input->GetInputValue(&strainparallel,gauss);
     228                strainperpendicular_input->GetInputValue(&strainperpendicular,gauss);
     229                levermanncoeff_input->GetInputValue(&propcoeff,gauss);
     230
     231                /*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 */
     232                calvingrate[iv]=propcoeff*strainparallel*strainperpendicular;
     233                if(calvingrate[iv]<0){
     234                        calvingrate[iv]=0;
     235                }
     236                calvingratex[iv]=calvingrate[iv]*vx/(vel+1.e-6);
     237                calvingratey[iv]=calvingrate[iv]*vy/(vel+1.e-6);
     238        }
     239
     240        /*Add input*/
     241        this->inputs->AddInput(new TriaInput(CalvingratexEnum,&calvingratex[0],P1Enum));
     242        this->inputs->AddInput(new TriaInput(CalvingrateyEnum,&calvingratey[0],P1Enum));
     243        this->inputs->AddInput(new TriaInput(CalvingCalvingrateEnum,&calvingrate[0],P1Enum));
     244
     245        /*Clean up and return*/
     246        delete gauss;
     247
     248}
     249/*}}}*/
    157250IssmDouble Tria::CharacteristicLength(void){/*{{{*/
    158251
    159252        return sqrt(2*this->GetArea());
    … …  
    163256        _error_("Not Implemented yet");
    164257}
    165258/*}}}*/
     259void       Tria::ComputeDeviatoricStressTensor(){/*{{{*/
     260
     261        IssmDouble  xyz_list[NUMVERTICES][3];
     262        IssmDouble  viscosity;
     263        IssmDouble  epsilon[3]; /* epsilon=[exx,eyy,exy];*/
     264        IssmDouble  tau_xx[NUMVERTICES];
     265        IssmDouble      tau_yy[NUMVERTICES];
     266        IssmDouble      tau_zz[NUMVERTICES]={0,0,0};
     267        IssmDouble  tau_xy[NUMVERTICES];
     268        IssmDouble      tau_xz[NUMVERTICES]={0,0,0};
     269        IssmDouble      tau_yz[NUMVERTICES]={0,0,0};
     270        GaussTria*  gauss=NULL;
     271        int domaintype,dim=2;
     272
     273        /* Get node coordinates and dof list: */
     274        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     275
     276        /*Retrieve all inputs we will be needing: */
     277        this->FindParam(&domaintype,DomainTypeEnum);
     278        if(domaintype!=Domain2DhorizontalEnum) _error_("deviatoric stress tensor calculation not implemented for mesh of type " <<EnumToStringx(domaintype));
     279        Input* vx_input=inputs->GetInput(VxEnum);             _assert_(vx_input);
     280        Input* vy_input=inputs->GetInput(VyEnum);             _assert_(vy_input);
     281
     282        /* Start looping on the number of vertices: */
     283        gauss=new GaussTria();
     284        for (int iv=0;iv<NUMVERTICES;iv++){
     285                gauss->GaussVertex(iv);
     286
     287                /*Compute strain rate and viscosity: */
     288                this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
     289                this->ViscositySSA(&viscosity,dim,&xyz_list[0][0],gauss,vx_input,vy_input);
     290
     291                /*Compute Stress*/
     292                tau_xx[iv]=2*viscosity*epsilon[0]; // tau = nu eps
     293                tau_yy[iv]=2*viscosity*epsilon[1];
     294                tau_xy[iv]=2*viscosity*epsilon[2];
     295        }
     296
     297        /*Add Stress tensor components into inputs*/
     298        this->inputs->AddInput(new TriaInput(DeviatoricStressxxEnum,&tau_xx[0],P1Enum));
     299        this->inputs->AddInput(new TriaInput(DeviatoricStressxyEnum,&tau_xy[0],P1Enum));
     300        this->inputs->AddInput(new TriaInput(DeviatoricStressxzEnum,&tau_xz[0],P1Enum));
     301        this->inputs->AddInput(new TriaInput(DeviatoricStressyyEnum,&tau_yy[0],P1Enum));
     302        this->inputs->AddInput(new TriaInput(DeviatoricStressyzEnum,&tau_yz[0],P1Enum));
     303        this->inputs->AddInput(new TriaInput(DeviatoricStresszzEnum,&tau_zz[0],P1Enum));
     304
     305        /*Clean up and return*/
     306        delete gauss;
     307}
     308/*}}}*/
    166309void       Tria::ComputeSigmaNN(){/*{{{*/
    167310
    168311        if(!IsOnBase()){
    … …  
    275418        delete gauss;
    276419}
    277420/*}}}*/
    278 void       Tria::ComputeDeviatoricStressTensor(){/*{{{*/
    279 
    280         IssmDouble  xyz_list[NUMVERTICES][3];
    281         IssmDouble  viscosity;
    282         IssmDouble  epsilon[3]; /* epsilon=[exx,eyy,exy];*/
    283         IssmDouble  tau_xx[NUMVERTICES];
    284         IssmDouble      tau_yy[NUMVERTICES];
    285         IssmDouble      tau_zz[NUMVERTICES]={0,0,0};
    286         IssmDouble  tau_xy[NUMVERTICES];
    287         IssmDouble      tau_xz[NUMVERTICES]={0,0,0};
    288         IssmDouble      tau_yz[NUMVERTICES]={0,0,0};
    289         GaussTria*  gauss=NULL;
    290         int domaintype,dim=2;
    291 
    292         /* Get node coordinates and dof list: */
    293         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    294 
    295         /*Retrieve all inputs we will be needing: */
    296         this->FindParam(&domaintype,DomainTypeEnum);
    297         if(domaintype!=Domain2DhorizontalEnum) _error_("deviatoric stress tensor calculation not implemented for mesh of type " <<EnumToStringx(domaintype));
    298         Input* vx_input=inputs->GetInput(VxEnum);             _assert_(vx_input);
    299         Input* vy_input=inputs->GetInput(VyEnum);             _assert_(vy_input);
    300 
    301         /* Start looping on the number of vertices: */
    302         gauss=new GaussTria();
    303         for (int iv=0;iv<NUMVERTICES;iv++){
    304                 gauss->GaussVertex(iv);
    305 
    306                 /*Compute strain rate and viscosity: */
    307                 this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    308                 this->ViscositySSA(&viscosity,dim,&xyz_list[0][0],gauss,vx_input,vy_input);
    309 
    310                 /*Compute Stress*/
    311                 tau_xx[iv]=2*viscosity*epsilon[0]; // tau = nu eps
    312                 tau_yy[iv]=2*viscosity*epsilon[1];
    313                 tau_xy[iv]=2*viscosity*epsilon[2];
    314         }
    315 
    316         /*Add Stress tensor components into inputs*/
    317         this->inputs->AddInput(new TriaInput(DeviatoricStressxxEnum,&tau_xx[0],P1Enum));
    318         this->inputs->AddInput(new TriaInput(DeviatoricStressxyEnum,&tau_xy[0],P1Enum));
    319         this->inputs->AddInput(new TriaInput(DeviatoricStressxzEnum,&tau_xz[0],P1Enum));
    320         this->inputs->AddInput(new TriaInput(DeviatoricStressyyEnum,&tau_yy[0],P1Enum));
    321         this->inputs->AddInput(new TriaInput(DeviatoricStressyzEnum,&tau_yz[0],P1Enum));
    322         this->inputs->AddInput(new TriaInput(DeviatoricStresszzEnum,&tau_zz[0],P1Enum));
    323 
    324         /*Clean up and return*/
    325         delete gauss;
    326 }
    327 /*}}}*/
    328 void       Tria::StrainRateparallel(){/*{{{*/
    329 
    330         IssmDouble *xyz_list = NULL;
    331         IssmDouble  epsilon[3];
    332         GaussTria* gauss=NULL;
    333         IssmDouble  vx,vy,vel;
    334         IssmDouble  strainxx;
    335         IssmDouble  strainxy;
    336         IssmDouble  strainyy;
    337         IssmDouble  strainparallel[NUMVERTICES];
    338 
    339         /* Get node coordinates and dof list: */
    340         this->GetVerticesCoordinates(&xyz_list);
    341 
    342         /*Retrieve all inputs we will need*/
    343         Input* vx_input=inputs->GetInput(VxEnum);                                  _assert_(vx_input);
    344         Input* vy_input=inputs->GetInput(VyEnum);                                  _assert_(vy_input);
    345 
    346         /* Start looping on the number of vertices: */
    347         gauss=new GaussTria();
    348         for (int iv=0;iv<NUMVERTICES;iv++){
    349                 gauss->GaussVertex(iv);
    350 
    351                 /* Get the value we need*/
    352                 vx_input->GetInputValue(&vx,gauss);
    353                 vy_input->GetInputValue(&vy,gauss);
    354                 vel=vx*vx+vy*vy;
    355 
    356                 /*Compute strain rate viscosity and pressure: */
    357                 this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
    358                 strainxx=epsilon[0];
    359                 strainyy=epsilon[1];
    360                 strainxy=epsilon[2];
    361 
    362                 /*strainparallel= Strain rate along the ice flow direction */
    363                 strainparallel[iv]=(vx*vx*(strainxx)+vy*vy*(strainyy)+2*vy*vx*strainxy)/(vel+1.e-6);
    364         }
    365 
    366         /*Add input*/
    367         this->inputs->AddInput(new TriaInput(StrainRateparallelEnum,&strainparallel[0],P1Enum));
    368 
    369         /*Clean up and return*/
    370         delete gauss;
    371         xDelete<IssmDouble>(xyz_list);
    372 }
    373 /*}}}*/
    374 void       Tria::StrainRateperpendicular(){/*{{{*/
    375 
    376         IssmDouble *xyz_list = NULL;
    377         GaussTria* gauss=NULL;
    378         IssmDouble  epsilon[3];
    379         IssmDouble  vx,vy,vel;
    380         IssmDouble  strainxx;
    381         IssmDouble  strainxy;
    382         IssmDouble  strainyy;
    383         IssmDouble  strainperpendicular[NUMVERTICES];
    384 
    385         /* Get node coordinates and dof list: */
    386         this->GetVerticesCoordinates(&xyz_list);
    387 
    388         /*Retrieve all inputs we will need*/
    389         Input* vx_input=inputs->GetInput(VxEnum);                                  _assert_(vx_input);
    390         Input* vy_input=inputs->GetInput(VyEnum);                                  _assert_(vy_input);
    391 
    392         /* Start looping on the number of vertices: */
    393         gauss=new GaussTria();
    394         for (int iv=0;iv<NUMVERTICES;iv++){
    395                 gauss->GaussVertex(iv);
    396 
    397                 /* Get the value we need*/
    398                 vx_input->GetInputValue(&vx,gauss);
    399                 vy_input->GetInputValue(&vy,gauss);
    400                 vel=vx*vx+vy*vy;
    401 
    402                 /*Compute strain rate viscosity and pressure: */
    403                 this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
    404                 strainxx=epsilon[0];
    405                 strainyy=epsilon[1];
    406                 strainxy=epsilon[2];
    407 
    408                 /*strainperpendicular= Strain rate perpendicular to the ice flow direction */
    409                 strainperpendicular[iv]=(vx*vx*(strainyy)+vy*vy*(strainxx)-2*vy*vx*strainxy)/(vel+1.e-6);
    410         }
    411 
    412         /*Add input*/
    413         this->inputs->AddInput(new TriaInput(StrainRateperpendicularEnum,&strainperpendicular[0],P1Enum));
    414 
    415         /*Clean up and return*/
    416         delete gauss;
    417         xDelete<IssmDouble>(xyz_list);
    418 }
    419 /*}}}*/
    420 void       Tria::CalvingRateLevermann(){/*{{{*/
    421 
    422         IssmDouble  xyz_list[NUMVERTICES][3];
    423         GaussTria* gauss=NULL;
    424         IssmDouble  vx,vy,vel;
    425         IssmDouble  strainparallel;
    426         IssmDouble  propcoeff;
    427         IssmDouble  strainperpendicular;
    428         IssmDouble  calvingratex[NUMVERTICES];
    429         IssmDouble  calvingratey[NUMVERTICES];
    430         IssmDouble  calvingrate[NUMVERTICES];
    431 
    432 
    433         /* Get node coordinates and dof list: */
    434         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    435 
    436         /*Retrieve all inputs and parameters we will need*/
    437         Input* vx_input=inputs->GetInput(VxEnum);                                                                                                                                               _assert_(vx_input);
    438         Input* vy_input=inputs->GetInput(VyEnum);                                                                                                                                               _assert_(vy_input);
    439         Input* strainparallel_input=inputs->GetInput(StrainRateparallelEnum);                                                           _assert_(strainparallel_input);
    440         Input* strainperpendicular_input=inputs->GetInput(StrainRateperpendicularEnum);                                 _assert_(strainperpendicular_input);
    441         Input* levermanncoeff_input=inputs->GetInput(CalvinglevermannCoeffEnum);                     _assert_(levermanncoeff_input);
    442 
    443         /* Start looping on the number of vertices: */
    444         gauss=new GaussTria();
    445         for (int iv=0;iv<NUMVERTICES;iv++){
    446                 gauss->GaussVertex(iv);
    447 
    448                 /* Get the value we need*/
    449                 vx_input->GetInputValue(&vx,gauss);
    450                 vy_input->GetInputValue(&vy,gauss);
    451                 vel=vx*vx+vy*vy;
    452                 strainparallel_input->GetInputValue(&strainparallel,gauss);
    453                 strainperpendicular_input->GetInputValue(&strainperpendicular,gauss);
    454                 levermanncoeff_input->GetInputValue(&propcoeff,gauss);
    455 
    456                 /*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 */
    457                 calvingrate[iv]=propcoeff*strainparallel*strainperpendicular;
    458                 if(calvingrate[iv]<0){
    459                         calvingrate[iv]=0;
    460                 }
    461                 calvingratex[iv]=calvingrate[iv]*vx/(vel+1.e-6);
    462                 calvingratey[iv]=calvingrate[iv]*vy/(vel+1.e-6);
    463         }
    464 
    465         /*Add input*/
    466         this->inputs->AddInput(new TriaInput(CalvingratexEnum,&calvingratex[0],P1Enum));
    467         this->inputs->AddInput(new TriaInput(CalvingrateyEnum,&calvingratey[0],P1Enum));
    468         this->inputs->AddInput(new TriaInput(CalvingCalvingrateEnum,&calvingrate[0],P1Enum));
    469 
    470         /*Clean up and return*/
    471         delete gauss;
    472 
    473 }
    474 /*}}}*/
    475421void       Tria::Configure(Elements* elementsin, Loads* loadsin,Nodes* nodesin,Vertices *verticesin,Materials* materialsin, Parameters* parametersin){/*{{{*/
    476422
    477423        /*go into parameters and get the analysis_counter: */
    … …  
    507453
    508454}
    509455/*}}}*/
     456void       Tria::ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){/*{{{*/
     457
     458        int    vertexpidlist[NUMVERTICES];
     459        IssmDouble grad_list[NUMVERTICES];
     460        Input* grad_input=NULL;
     461
     462        Input* input=inputs->GetInput(enum_type);
     463        if (!input) _error_("Input " << EnumToStringx(enum_type) << " not found");
     464        if (input->ObjectEnum()!=ControlInputEnum) _error_("Input " << EnumToStringx(enum_type) << " is not a ControlInput");
     465
     466        GradientIndexing(&vertexpidlist[0],control_index);
     467        for(int i=0;i<NUMVERTICES;i++) grad_list[i]=gradient[vertexpidlist[i]];
     468        grad_input=new TriaInput(GradientEnum,grad_list,P1Enum);
     469
     470        ((ControlInput*)input)->SetGradient(grad_input);
     471
     472}/*}}}*/
     473void       Tria::ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum){/*{{{*/
     474
     475        Input* input=inputs->GetInput(control_enum);
     476        if (!input) _error_("Input " << EnumToStringx(control_enum) << " not found");
     477        if (input->ObjectEnum()!=ControlInputEnum) _error_("Input " << EnumToStringx(control_enum) << " is not a ControlInput");
     478
     479        int         sidlist[NUMVERTICES];
     480        int         connectivity[NUMVERTICES];
     481        IssmPDouble values[NUMVERTICES];
     482        IssmPDouble gradients[NUMVERTICES];
     483        IssmDouble  value,gradient;
     484
     485        this->GetVerticesConnectivityList(&connectivity[0]);
     486        this->GetVerticesSidList(&sidlist[0]);
     487
     488        GaussTria* gauss=new GaussTria();
     489        for (int iv=0;iv<NUMVERTICES;iv++){
     490                gauss->GaussVertex(iv);
     491
     492                ((ControlInput*)input)->GetInputValue(&value,gauss);
     493                ((ControlInput*)input)->GetGradientValue(&gradient,gauss);
     494
     495                values[iv]    = reCast<IssmPDouble>(value)/reCast<IssmPDouble>(connectivity[iv]);
     496                gradients[iv] = reCast<IssmPDouble>(gradient)/reCast<IssmPDouble>(connectivity[iv]);
     497        }
     498        delete gauss;
     499
     500        vector_control->SetValues(NUMVERTICES,&sidlist[0],&values[0],ADD_VAL);
     501        vector_gradient->SetValues(NUMVERTICES,&sidlist[0],&gradients[0],ADD_VAL);
     502
     503}/*}}}*/
    510504void       Tria::Delta18oParameterization(void){/*{{{*/
    511505
    512506        int        i;
    … …  
    577571        delete gauss;
    578572}
    579573/*}}}*/
     574int        Tria::EdgeOnBaseIndex(void){/*{{{*/
     575
     576        IssmDouble values[NUMVERTICES];
     577        int        indices[3][2] = {{1,2},{2,0},{0,1}};
     578
     579        /*Retrieve all inputs and parameters*/
     580        GetInputListOnVertices(&values[0],MeshVertexonbaseEnum);
     581
     582        for(int i=0;i<3;i++){
     583                if(values[indices[i][0]] == 1. && values[indices[i][1]] == 1.){
     584                        return i;
     585                }
     586        }
     587
     588        _printf_("list of vertices on bed: "<<values[0]<<" "<<values[1]<<" "<<values[2]);
     589        _error_("Could not find 2 vertices on bed");
     590}
     591/*}}}*/
     592void       Tria::EdgeOnBaseIndices(int* pindex1,int* pindex2){/*{{{*/
     593
     594        IssmDouble values[NUMVERTICES];
     595        int        indices[3][2] = {{1,2},{2,0},{0,1}};
     596
     597        /*Retrieve all inputs and parameters*/
     598        GetInputListOnVertices(&values[0],MeshVertexonbaseEnum);
     599
     600        for(int i=0;i<3;i++){
     601                if(values[indices[i][0]] == 1. && values[indices[i][1]] == 1.){
     602                        *pindex1 = indices[i][0];
     603                        *pindex2 = indices[i][1];
     604                        return;
     605                }
     606        }
     607
     608        _printf_("list of vertices on bed: "<<values[0]<<" "<<values[1]<<" "<<values[2]);
     609        _error_("Could not find 2 vertices on bed");
     610}
     611/*}}}*/
     612int        Tria::EdgeOnSurfaceIndex(void){/*{{{*/
     613
     614        IssmDouble values[NUMVERTICES];
     615        int        indices[3][2] = {{1,2},{2,0},{0,1}};
     616
     617        /*Retrieve all inputs and parameters*/
     618        GetInputListOnVertices(&values[0],MeshVertexonsurfaceEnum);
     619
     620        for(int i=0;i<3;i++){
     621                if(values[indices[i][0]] == 1. && values[indices[i][1]] == 1.){
     622                        return i;
     623                }
     624        }
     625
     626        _printf_("list of vertices on surface: "<<values[0]<<" "<<values[1]<<" "<<values[2]);
     627        _error_("Could not find 2 vertices on surface");
     628}
     629/*}}}*/
     630void       Tria::EdgeOnSurfaceIndices(int* pindex1,int* pindex2){/*{{{*/
     631
     632        IssmDouble values[NUMVERTICES];
     633        int        indices[3][2] = {{1,2},{2,0},{0,1}};
     634
     635        /*Retrieve all inputs and parameters*/
     636        GetInputListOnVertices(&values[0],MeshVertexonsurfaceEnum);
     637
     638        for(int i=0;i<3;i++){
     639                if(values[indices[i][0]] == 1. && values[indices[i][1]] == 1.){
     640                        *pindex1 = indices[i][0];
     641                        *pindex2 = indices[i][1];
     642                        return;
     643                }
     644        }
     645
     646        _printf_("list of vertices on surface: "<<values[0]<<" "<<values[1]<<" "<<values[2]);
     647        _error_("Could not find 2 vertices on surface");
     648}
     649/*}}}*/
     650void       Tria::ElementResponse(IssmDouble* presponse,int response_enum){/*{{{*/
     651
     652        switch(response_enum){
     653                case MaterialsRheologyBbarEnum:
     654                        *presponse=this->material->GetBbar();
     655                        break;
     656
     657                case VelEnum:{
     658
     659                        /*Get input:*/
     660                        IssmDouble vel;
     661                        Input* vel_input;
     662
     663                        vel_input=this->inputs->GetInput(VelEnum); _assert_(vel_input);
     664                        vel_input->GetInputAverage(&vel);
     665
     666                        /*Assign output pointers:*/
     667                        *presponse=vel;}
     668                        break;
     669                default: 
     670                        _error_("Response type " << EnumToStringx(response_enum) << " not supported yet!");
     671        }
     672
     673}
     674/*}}}*/
    580675void       Tria::ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){/*{{{*/
    581676
    582677        IssmDouble xyz_list[NUMVERTICES][3];
    … …  
    604699        return this->element_type;
    605700}
    606701/*}}}*/
    607 int        Tria::ObjectEnum(void){/*{{{*/
     702void       Tria::FSContactMigration(Vector<IssmDouble>* vertexgrounded,Vector<IssmDouble>* vertexfloating){/*{{{*/
    608703
    609         return TriaEnum;
     704        if(!IsOnBase()) return;
    610705
     706        int approximation;
     707        inputs->GetInputValue(&approximation,ApproximationEnum);
     708
     709        if(approximation==HOApproximationEnum || approximation==SSAApproximationEnum || approximation==SSAHOApproximationEnum){
     710                for(int i=0;i<NUMVERTICES;i++){
     711                        vertexgrounded->SetValue(vertices[i]->Pid(),+9999.,INS_VAL);
     712                        vertexfloating->SetValue(vertices[i]->Pid(),+9999.,INS_VAL);
     713                }
     714        }
     715        else{
     716                /*Intermediaries*/
     717                IssmDouble* xyz_list = NULL;
     718                IssmDouble* xyz_list_base = NULL;
     719                IssmDouble  pressure,water_pressure,sigma_nn,viscosity,bed,base;
     720                IssmDouble  bed_normal[2];
     721                IssmDouble  epsilon[3]; /* epsilon=[exx,eyy,exy];*/
     722                IssmDouble  surface=0,value=0;
     723                bool grounded;
     724
     725                /* Get node coordinates and dof list: */
     726                GetVerticesCoordinates(&xyz_list);
     727                GetVerticesCoordinatesBase(&xyz_list_base);
     728
     729                /*Retrieve all inputs we will be needing: */
     730                Input* pressure_input = inputs->GetInput(PressureEnum); _assert_(pressure_input);
     731                Input* base_input     = inputs->GetInput(BaseEnum);     _assert_(base_input);
     732                Input* bed_input      = inputs->GetInput(BedEnum);      _assert_(bed_input);
     733                Input* vx_input       = inputs->GetInput(VxEnum);       _assert_(vx_input);
     734                Input* vy_input       = inputs->GetInput(VyEnum);       _assert_(vy_input);
     735
     736                /*Create gauss point in the middle of the basal edge*/
     737                Gauss* gauss=NewGaussBase(1);
     738                gauss->GaussPoint(0);
     739
     740                if(!IsFloating()){
     741                        /*Check for basal force only if grounded and touching GL*/
     742                        //              if(this->inputs->Min(MaskGroundediceLevelsetEnum)==0.){
     743                        this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
     744                        this->ViscosityFS(&viscosity,2,xyz_list,gauss,vx_input,vy_input,NULL);
     745                        pressure_input->GetInputValue(&pressure, gauss);
     746                        base_input->GetInputValue(&base, gauss);
     747
     748                        /*Compute Stress*/
     749                        IssmDouble sigma_xx=2.*viscosity*epsilon[0]-pressure;
     750                        IssmDouble sigma_yy=2.*viscosity*epsilon[1]-pressure;
     751                        IssmDouble sigma_xy=2.*viscosity*epsilon[2];
     752
     753                        /*Get normal vector to the bed */
     754                        NormalBase(&bed_normal[0],xyz_list_base);
     755
     756                        /*basalforce*/
     757                        sigma_nn = sigma_xx*bed_normal[0]*bed_normal[0] + sigma_yy*bed_normal[1]*bed_normal[1] + 2.*sigma_xy*bed_normal[0]*bed_normal[1];
     758
     759                        /*Compute water pressure*/
     760                        IssmDouble rho_ice   = matpar->GetRhoIce();
     761                        IssmDouble rho_water = matpar->GetRhoWater();
     762                        IssmDouble gravity   = matpar->GetG();
     763                        water_pressure=gravity*rho_water*base;
     764
     765                        /*Compare basal stress to water pressure and determine whether it should ground*/
     766                        if (sigma_nn<water_pressure) grounded=true;
     767                        else                         grounded=false;
     768                }
     769                else{
     770                        /*Check for basal elevation if floating*/
     771                        base_input->GetInputValue(&base, gauss);
     772                        bed_input->GetInputValue(&bed, gauss);
     773                        if(base<bed) grounded=true;
     774                        else         grounded=false;
     775                }
     776                for(int i=0;i<NUMVERTICES;i++){
     777                        if(grounded) vertexgrounded->SetValue(vertices[i]->Pid(),+1.,INS_VAL);
     778                        else         vertexfloating->SetValue(vertices[i]->Pid(),+1.,INS_VAL);
     779                }
     780
     781                /*clean up*/
     782                delete gauss;
     783                xDelete<IssmDouble>(xyz_list);
     784                xDelete<IssmDouble>(xyz_list_base);
     785        }
    611786}
    612787/*}}}*/
    613788IssmDouble Tria::GetArea(void){/*{{{*/
    … …  
    667842
    668843}
    669844/*}}}*/
    670 void        Tria::GetLevelsetPositivePart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlynegative,IssmDouble* gl){/*{{{*/
    671        
    672         /*Computeportion of the element that has a positive levelset*/
    673 
    674         bool               negative=true;
    675         int                point;
    676         const IssmPDouble  epsilon= 1.e-15;
    677         IssmDouble         f1,f2;
    678 
    679         /*Be sure that values are not zero*/
    680         if(gl[0]==0.) gl[0]=gl[0]+epsilon;
    681         if(gl[1]==0.) gl[1]=gl[1]+epsilon;
    682         if(gl[2]==0.) gl[2]=gl[2]+epsilon;
    683 
    684         /*Check that not all nodes are positive or negative*/
    685         if(gl[0]>0 && gl[1]>0 && gl[2]>0){ // All positive
    686                 point=0;
    687                 f1=1.;
    688                 f2=1.;
    689         }
    690         else if(gl[0]<0 && gl[1]<0 && gl[2]<0){ //All negative
    691                 point=0;
    692                 f1=0.;
    693                 f2=0.;
    694         }
    695         else{
    696                 if(gl[0]*gl[1]*gl[2]<0) negative=false;
    697 
    698                 if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
    699                         point=2;
    700                         f1=gl[2]/(gl[2]-gl[0]);
    701                         f2=gl[2]/(gl[2]-gl[1]);
    702                 }
    703                 else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
    704                         point=0;
    705                         f1=gl[0]/(gl[0]-gl[1]);
    706                         f2=gl[0]/(gl[0]-gl[2]);
    707                 }
    708                 else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
    709                         point=1;
    710                         f1=gl[1]/(gl[1]-gl[2]);
    711                         f2=gl[1]/(gl[1]-gl[0]);
    712                 }
    713         }
    714         *point1=point;
    715         *fraction1=f1;
    716         *fraction2=f2;
    717         *mainlynegative=negative;
    718 }
    719 /*}}}*/
    720845void       Tria::GetGroundedPart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlyfloating){/*{{{*/
    721846        /*Computeportion of the element that is grounded*/
    722847
    … …  
    8881013        return phi;
    8891014}
    8901015/*}}}*/
    891 void       Tria::GetVerticesCoordinatesBase(IssmDouble** pxyz_list){/*{{{*/
    892 
    893         int        indices[2];
    894         IssmDouble xyz_list[NUMVERTICES][3];
    895 
    896         /*Element XYZ list*/
    897         ::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES);
    898 
    899         /*Allocate Output*/
    900         IssmDouble* xyz_list_edge = xNew<IssmDouble>(2*3);
    901         this->EdgeOnBaseIndices(&indices[0],&indices[1]);
    902         for(int i=0;i<2;i++) for(int j=0;j<2;j++) xyz_list_edge[i*3+j]=xyz_list[indices[i]][j];
    903 
    904         /*Assign output pointer*/
    905         *pxyz_list = xyz_list_edge;
    906 
    907 }/*}}}*/
    908 void       Tria::GetVerticesCoordinatesTop(IssmDouble** pxyz_list){/*{{{*/
    909 
    910         int        indices[2];
    911         IssmDouble xyz_list[NUMVERTICES][3];
    912 
    913         /*Element XYZ list*/
    914         ::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES);
    915 
    916         /*Allocate Output*/
    917         IssmDouble* xyz_list_edge = xNew<IssmDouble>(2*3);
    918         this->EdgeOnSurfaceIndices(&indices[0],&indices[1]);
    919         for(int i=0;i<2;i++) for(int j=0;j<2;j++) xyz_list_edge[i*3+j]=xyz_list[indices[i]][j];
    920 
    921         /*Assign output pointer*/
    922         *pxyz_list = xyz_list_edge;
    923 
    924 }/*}}}*/
    925 void       Tria::NormalSection(IssmDouble* normal,IssmDouble* xyz_list){/*{{{*/
    926 
    927         /*Build unit outward pointing vector*/
    928         IssmDouble vector[2];
    929         IssmDouble norm;
    930 
    931         vector[0]=xyz_list[1*3+0] - xyz_list[0*3+0];
    932         vector[1]=xyz_list[1*3+1] - xyz_list[0*3+1];
    933 
    934         norm=sqrt(vector[0]*vector[0] + vector[1]*vector[1]);
    935 
    936         normal[0]= + vector[1]/norm;
    937         normal[1]= - vector[0]/norm;
    938 }
    939 /*}}}*/
    940 void       Tria::ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum){/*{{{*/
    941 
    942         int         normal_orientation=0;
    943         IssmDouble  s1,s2;
    944         IssmDouble  levelset[NUMVERTICES];
    945 
    946         /*Recover parameters and values*/
    947         IssmDouble* xyz_zero = xNew<IssmDouble>(2*3);
    948         GetInputListOnVertices(&levelset[0],levelsetenum);
    949 
    950         if(levelset[0]*levelset[1]>0.){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
    951                 /*Portion of the segments*/
    952                 s1=levelset[2]/(levelset[2]-levelset[1]);
    953                 s2=levelset[2]/(levelset[2]-levelset[0]);
    954 
    955                 if(levelset[2]<0.) normal_orientation=1; //orientation of quadrangle depending on distribution of levelsetfunction
    956                 /*New point 1*/
    957                 xyz_zero[3*normal_orientation+0]=xyz_list[2*3+0]+s1*(xyz_list[1*3+0]-xyz_list[2*3+0]);
    958                 xyz_zero[3*normal_orientation+1]=xyz_list[2*3+1]+s1*(xyz_list[1*3+1]-xyz_list[2*3+1]);
    959                 xyz_zero[3*normal_orientation+2]=xyz_list[2*3+2]+s1*(xyz_list[1*3+2]-xyz_list[2*3+2]);
    960 
    961                 /*New point 0*/
    962                 xyz_zero[3*(1-normal_orientation)+0]=xyz_list[2*3+0]+s2*(xyz_list[0*3+0]-xyz_list[2*3+0]);
    963                 xyz_zero[3*(1-normal_orientation)+1]=xyz_list[2*3+1]+s2*(xyz_list[0*3+1]-xyz_list[2*3+1]);
    964                 xyz_zero[3*(1-normal_orientation)+2]=xyz_list[2*3+2]+s2*(xyz_list[0*3+2]-xyz_list[2*3+2]);
    965         }
    966         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
    967                 /*Portion of the segments*/
    968                 s1=levelset[0]/(levelset[0]-levelset[2]);
    969                 s2=levelset[0]/(levelset[0]-levelset[1]);
    970 
    971                 if(levelset[0]<0.) normal_orientation=1;
    972                 /*New point 1*/
    973                 xyz_zero[3*normal_orientation+0]=xyz_list[0*3+0]+s1*(xyz_list[2*3+0]-xyz_list[0*3+0]);
    974                 xyz_zero[3*normal_orientation+1]=xyz_list[0*3+1]+s1*(xyz_list[2*3+1]-xyz_list[0*3+1]);
    975                 xyz_zero[3*normal_orientation+2]=xyz_list[0*3+2]+s1*(xyz_list[2*3+2]-xyz_list[0*3+2]);
    976 
    977                 /*New point 2*/
    978                 xyz_zero[3*(1-normal_orientation)+0]=xyz_list[0*3+0]+s2*(xyz_list[1*3+0]-xyz_list[0*3+0]);
    979                 xyz_zero[3*(1-normal_orientation)+1]=xyz_list[0*3+1]+s2*(xyz_list[1*3+1]-xyz_list[0*3+1]);
    980                 xyz_zero[3*(1-normal_orientation)+2]=xyz_list[0*3+2]+s2*(xyz_list[1*3+2]-xyz_list[0*3+2]);
    981         }
    982         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
    983                 /*Portion of the segments*/
    984                 s1=levelset[1]/(levelset[1]-levelset[0]);
    985                 s2=levelset[1]/(levelset[1]-levelset[2]);
    986 
    987                 if(levelset[1]<0.) normal_orientation=1;
    988                 /*New point 0*/
    989                 xyz_zero[3*normal_orientation+0]=xyz_list[1*3+0]+s1*(xyz_list[0*3+0]-xyz_list[1*3+0]);
    990                 xyz_zero[3*normal_orientation+1]=xyz_list[1*3+1]+s1*(xyz_list[0*3+1]-xyz_list[1*3+1]);
    991                 xyz_zero[3*normal_orientation+2]=xyz_list[1*3+2]+s1*(xyz_list[0*3+2]-xyz_list[1*3+2]);
    992 
    993                 /*New point 2*/
    994                 xyz_zero[3*(1-normal_orientation)+0]=xyz_list[1*3+0]+s2*(xyz_list[2*3+0]-xyz_list[1*3+0]);
    995                 xyz_zero[3*(1-normal_orientation)+1]=xyz_list[1*3+1]+s2*(xyz_list[2*3+1]-xyz_list[1*3+1]);
    996                 xyz_zero[3*(1-normal_orientation)+2]=xyz_list[1*3+2]+s2*(xyz_list[2*3+2]-xyz_list[1*3+2]);
    997         }
    998         else if(levelset[0]==0. && levelset[1]==0.){ //front is on point 0 and 1
    999                 xyz_zero[3*0+0]=xyz_list[0*3+0];
    1000                 xyz_zero[3*0+1]=xyz_list[0*3+1];
    1001                 xyz_zero[3*0+2]=xyz_list[0*3+2];
    1002 
    1003                 /*New point 2*/
    1004                 xyz_zero[3*1+0]=xyz_list[1*3+0];
    1005                 xyz_zero[3*1+1]=xyz_list[1*3+1];
    1006                 xyz_zero[3*1+2]=xyz_list[1*3+2];
    1007         }
    1008         else if(levelset[0]==0. && levelset[2]==0.){ //front is on point 0 and 1
    1009                 xyz_zero[3*0+0]=xyz_list[2*3+0];
    1010                 xyz_zero[3*0+1]=xyz_list[2*3+1];
    1011                 xyz_zero[3*0+2]=xyz_list[2*3+2];
    1012 
    1013                 /*New point 2*/
    1014                 xyz_zero[3*1+0]=xyz_list[0*3+0];
    1015                 xyz_zero[3*1+1]=xyz_list[0*3+1];
    1016                 xyz_zero[3*1+2]=xyz_list[0*3+2];
    1017         }
    1018         else if(levelset[1]==0. && levelset[2]==0.){ //front is on point 0 and 1
    1019                 xyz_zero[3*0+0]=xyz_list[1*3+0];
    1020                 xyz_zero[3*0+1]=xyz_list[1*3+1];
    1021                 xyz_zero[3*0+2]=xyz_list[1*3+2];
    1022 
    1023                 /*New point 2*/
    1024                 xyz_zero[3*1+0]=xyz_list[2*3+0];
    1025                 xyz_zero[3*1+1]=xyz_list[2*3+1];
    1026                 xyz_zero[3*1+2]=xyz_list[2*3+2];
    1027         }
    1028         else _error_("Case not covered");
    1029 
    1030         /*Assign output pointer*/
    1031         *pxyz_zero= xyz_zero;
    1032 }
    1033 /*}}}*/
    10341016void       Tria::GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){/*{{{*/
    10351017       
    10361018        /* Intermediaries */
    … …  
    10711053
    10721054        xDelete<int>(indicesfront);
    10731055}/*}}}*/
     1056void       Tria::GetInputValue(IssmDouble* pvalue,Node* node,int enumtype){/*{{{*/
     1057
     1058        Input* input=inputs->GetInput(enumtype);
     1059        if(!input) _error_("No input of type " << EnumToStringx(enumtype) << " found in tria");
     1060
     1061        GaussTria* gauss=new GaussTria();
     1062        gauss->GaussVertex(this->GetNodeIndex(node));
     1063
     1064        input->GetInputValue(pvalue,gauss);
     1065        delete gauss;
     1066}
     1067/*}}}*/
    10741068void       Tria::GetLevelCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum,IssmDouble level){/*{{{*/
    10751069
    10761070        /* Intermediaries */
    … …  
    11111105
    11121106        xDelete<int>(indicesfront);
    11131107}/*}}}*/
     1108void       Tria::GetLevelsetPositivePart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlynegative,IssmDouble* gl){/*{{{*/
     1109       
     1110        /*Computeportion of the element that has a positive levelset*/
     1111
     1112        bool               negative=true;
     1113        int                point;
     1114        const IssmPDouble  epsilon= 1.e-15;
     1115        IssmDouble         f1,f2;
     1116
     1117        /*Be sure that values are not zero*/
     1118        if(gl[0]==0.) gl[0]=gl[0]+epsilon;
     1119        if(gl[1]==0.) gl[1]=gl[1]+epsilon;
     1120        if(gl[2]==0.) gl[2]=gl[2]+epsilon;
     1121
     1122        /*Check that not all nodes are positive or negative*/
     1123        if(gl[0]>0 && gl[1]>0 && gl[2]>0){ // All positive
     1124                point=0;
     1125                f1=1.;
     1126                f2=1.;
     1127        }
     1128        else if(gl[0]<0 && gl[1]<0 && gl[2]<0){ //All negative
     1129                point=0;
     1130                f1=0.;
     1131                f2=0.;
     1132        }
     1133        else{
     1134                if(gl[0]*gl[1]*gl[2]<0) negative=false;
     1135
     1136                if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
     1137                        point=2;
     1138                        f1=gl[2]/(gl[2]-gl[0]);
     1139                        f2=gl[2]/(gl[2]-gl[1]);
     1140                }
     1141                else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
     1142                        point=0;
     1143                        f1=gl[0]/(gl[0]-gl[1]);
     1144                        f2=gl[0]/(gl[0]-gl[2]);
     1145                }
     1146                else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
     1147                        point=1;
     1148                        f1=gl[1]/(gl[1]-gl[2]);
     1149                        f2=gl[1]/(gl[1]-gl[0]);
     1150                }
     1151        }
     1152        *point1=point;
     1153        *fraction1=f1;
     1154        *fraction2=f2;
     1155        *mainlynegative=negative;
     1156}
     1157/*}}}*/
     1158Node*      Tria::GetNode(int node_number){/*{{{*/
     1159        _assert_(node_number>=0);
     1160        _assert_(node_number<this->NumberofNodes(this->element_type));
     1161        return this->nodes[node_number];
     1162
     1163}/*}}}*/
    11141164int        Tria::GetNodeIndex(Node* node){/*{{{*/
    11151165
    11161166        _assert_(nodes);
    … …  
    11341184        return NUMVERTICES;
    11351185}
    11361186/*}}}*/
    1137 void       Tria::GetInputValue(IssmDouble* pvalue,Node* node,int enumtype){/*{{{*/
     1187void       Tria::GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution, int enum_type){/*{{{*/
    11381188
    1139         Input* input=inputs->GetInput(enumtype);
    1140         if(!input) _error_("No input of type " << EnumToStringx(enumtype) << " found in tria");
     1189        int        *doflist = NULL;
     1190        IssmDouble  value;
    11411191
     1192        /*Fetch number of nodes for this finite element*/
     1193        int numnodes = this->NumberofNodes(this->element_type);
     1194
     1195        /*Fetch dof list and allocate solution vector*/
     1196        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     1197        IssmDouble* values = xNew<IssmDouble>(numnodes);
     1198
     1199        /*Get inputs*/
     1200        Input* enum_input=inputs->GetInput(enum_type); _assert_(enum_input);
     1201
     1202        /*Ok, we have the values, fill in the array: */
    11421203        GaussTria* gauss=new GaussTria();
    1143         gauss->GaussVertex(this->GetNodeIndex(node));
     1204        for(int i=0;i<numnodes;i++){
     1205                gauss->GaussNode(this->element_type,i);
    11441206
    1145         input->GetInputValue(pvalue,gauss);
     1207                enum_input->GetInputValue(&value,gauss);
     1208                values[i]=value;
     1209        }
     1210
     1211        solution->SetValues(numnodes,doflist,values,INS_VAL);
     1212
     1213        /*Free ressources:*/
     1214        xDelete<int>(doflist);
     1215        xDelete<IssmDouble>(values);
    11461216        delete gauss;
    11471217}
    11481218/*}}}*/
    1149 Node*      Tria::GetNode(int node_number){/*{{{*/
    1150         _assert_(node_number>=0);
    1151         _assert_(node_number<this->NumberofNodes(this->element_type));
    1152         return this->nodes[node_number];
     1219void       Tria::GetVectorFromControlInputs(Vector<IssmDouble>* vector,int control_enum,int control_index,const char* data,bool onsid){/*{{{*/
    11531220
     1221        int vertexidlist[NUMVERTICES];
     1222        Input *input=NULL;
     1223
     1224        /*Get out if this is not an element input*/
     1225        if(!IsInput(control_enum)) return;
     1226
     1227        /*Prepare index list*/
     1228        GradientIndexing(&vertexidlist[0],control_index,onsid);
     1229
     1230        /*Get input (either in element or material)*/
     1231        input=(Input*)this->inputs->GetInput(control_enum);   _assert_(input);
     1232
     1233        /*Check that it is a ControlInput*/
     1234        if (input->ObjectEnum()!=ControlInputEnum){
     1235                _error_("input " << EnumToStringx(control_enum) << " is not a ControlInput");
     1236        }
     1237
     1238        ((ControlInput*)input)->GetVectorFromInputs(vector,&vertexidlist[0],data);
     1239}
     1240/*}}}*/
     1241void       Tria::GetVerticesCoordinatesBase(IssmDouble** pxyz_list){/*{{{*/
     1242
     1243        int        indices[2];
     1244        IssmDouble xyz_list[NUMVERTICES][3];
     1245
     1246        /*Element XYZ list*/
     1247        ::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES);
     1248
     1249        /*Allocate Output*/
     1250        IssmDouble* xyz_list_edge = xNew<IssmDouble>(2*3);
     1251        this->EdgeOnBaseIndices(&indices[0],&indices[1]);
     1252        for(int i=0;i<2;i++) for(int j=0;j<2;j++) xyz_list_edge[i*3+j]=xyz_list[indices[i]][j];
     1253
     1254        /*Assign output pointer*/
     1255        *pxyz_list = xyz_list_edge;
     1256
    11541257}/*}}}*/
     1258void       Tria::GetVerticesCoordinatesTop(IssmDouble** pxyz_list){/*{{{*/
     1259
     1260        int        indices[2];
     1261        IssmDouble xyz_list[NUMVERTICES][3];
     1262
     1263        /*Element XYZ list*/
     1264        ::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES);
     1265
     1266        /*Allocate Output*/
     1267        IssmDouble* xyz_list_edge = xNew<IssmDouble>(2*3);
     1268        this->EdgeOnSurfaceIndices(&indices[0],&indices[1]);
     1269        for(int i=0;i<2;i++) for(int j=0;j<2;j++) xyz_list_edge[i*3+j]=xyz_list[indices[i]][j];
     1270
     1271        /*Assign output pointer*/
     1272        *pxyz_list = xyz_list_edge;
     1273
     1274}/*}}}*/
     1275bool       Tria::HasEdgeOnBase(){/*{{{*/
     1276
     1277        IssmDouble values[NUMVERTICES];
     1278        IssmDouble sum;
     1279
     1280        /*Retrieve all inputs and parameters*/
     1281        GetInputListOnVertices(&values[0],MeshVertexonbaseEnum);
     1282        sum = values[0]+values[1]+values[2];
     1283
     1284        _assert_(sum==0. || sum==1. || sum==2.);
     1285
     1286        if(sum==3.)  _error_("Two edges on bed not supported yet...");
     1287
     1288        if(sum>1.){
     1289                return true;
     1290        }
     1291        else{
     1292                return false;
     1293        }
     1294}
     1295/*}}}*/
     1296bool       Tria::HasEdgeOnSurface(){/*{{{*/
     1297
     1298        IssmDouble values[NUMVERTICES];
     1299        IssmDouble sum;
     1300
     1301        /*Retrieve all inputs and parameters*/
     1302        GetInputListOnVertices(&values[0],MeshVertexonsurfaceEnum);
     1303        sum = values[0]+values[1]+values[2];
     1304
     1305        _assert_(sum==0. || sum==1. || sum==2.);
     1306
     1307        if(sum==3.)  _error_("Two edges on surface not supported yet...");
     1308
     1309        if(sum>1.){
     1310                return true;
     1311        }
     1312        else{
     1313                return false;
     1314        }
     1315}
     1316/*}}}*/
     1317IssmDouble Tria::IceVolume(void){/*{{{*/
     1318
     1319        /*The volume of a troncated prism is base * 1/3 sum(length of edges)*/
     1320        IssmDouble base,surface,bed;
     1321        IssmDouble xyz_list[NUMVERTICES][3];
     1322
     1323        if(!IsIceInElement())return 0;
     1324
     1325        /*First get back the area of the base*/
     1326        base=this->GetArea();
     1327
     1328        /*Now get the average height*/
     1329        Input* surface_input = inputs->GetInput(SurfaceEnum); _assert_(surface_input);
     1330        Input* base_input     = inputs->GetInput(BaseEnum);     _assert_(base_input);
     1331        surface_input->GetInputAverage(&surface);
     1332        base_input->GetInputAverage(&bed);
     1333
     1334        /*Return: */
     1335        int domaintype;
     1336        parameters->FindParam(&domaintype,DomainTypeEnum);
     1337        if(domaintype==Domain2DverticalEnum){
     1338          return base;
     1339        }
     1340        else{
     1341          return base*(surface-bed);
     1342        }
     1343}
     1344/*}}}*/
     1345IssmDouble Tria::IceVolumeAboveFloatation(void){/*{{{*/
     1346
     1347        /*The volume above floatation: H + rho_water/rho_ice * bathymetry */
     1348        IssmDouble rho_ice,rho_water;
     1349        IssmDouble base,surface,bed,bathymetry;
     1350        IssmDouble xyz_list[NUMVERTICES][3];
     1351
     1352        if(!IsIceInElement() || IsFloating())return 0;
     1353
     1354        rho_ice=matpar->GetRhoIce();
     1355        rho_water=matpar->GetRhoWater();
     1356        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     1357
     1358        /*First calculate the area of the base (cross section triangle)
     1359         * http://en.wikipedia.org/wiki/Triangle
     1360         * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
     1361        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]));
     1362
     1363        /*Now get the average height and bathymetry*/
     1364        Input* surface_input    = inputs->GetInput(SurfaceEnum);    _assert_(surface_input);
     1365        Input* base_input        = inputs->GetInput(BaseEnum);        _assert_(base_input);
     1366        Input* bed_input = inputs->GetInput(BedEnum); _assert_(bed_input);
     1367        surface_input->GetInputAverage(&surface);
     1368        base_input->GetInputAverage(&bed);
     1369        bed_input->GetInputAverage(&bathymetry);
     1370       
     1371        /*Return: */
     1372        return base*(surface-bed+min(rho_water/rho_ice*bathymetry,0.));
     1373}
     1374/*}}}*/
     1375void       Tria::InputControlUpdate(IssmDouble scalar,bool save_parameter){/*{{{*/
     1376
     1377        /*Intermediary*/
     1378        int    num_controls;
     1379        int*   control_type=NULL;
     1380        Input* input=NULL;
     1381
     1382        /*retrieve some parameters: */
     1383        this->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
     1384        this->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
     1385
     1386        for(int i=0;i<num_controls;i++){
     1387                input=(Input*)this->inputs->GetInput(control_type[i]);   _assert_(input);
     1388                if (input->ObjectEnum()!=ControlInputEnum){
     1389                        _error_("input " << EnumToStringx(control_type[i]) << " is not a ControlInput");
     1390                }
     1391
     1392                ((ControlInput*)input)->UpdateValue(scalar);
     1393                ((ControlInput*)input)->Constrain();
     1394                if (save_parameter) ((ControlInput*)input)->SaveValue();
     1395
     1396        }
     1397
     1398        /*Clean up and return*/
     1399        xDelete<int>(control_type);
     1400}
     1401/*}}}*/
    11551402void       Tria::InputDepthAverageAtBase(int enum_type,int average_enum_type){/*{{{*/
    11561403
    11571404        /*New input*/
    … …  
    13701617
    13711618}
    13721619/*}}}*/
     1620bool       Tria::IsFaceOnBoundary(void){/*{{{*/
     1621
     1622        IssmDouble values[NUMVERTICES];
     1623        IssmDouble sum;
     1624
     1625        /*Retrieve all inputs and parameters*/
     1626        GetInputListOnVertices(&values[0],MeshVertexonboundaryEnum);
     1627        sum = values[0]+values[1]+values[2];
     1628
     1629        _assert_(sum==0. || sum==1. || sum==2.);
     1630
     1631        if(sum==3.)  _error_("Two edges on boundary not supported yet...");
     1632
     1633        if(sum>1.){
     1634                return true;
     1635        }
     1636        else{
     1637                return false;
     1638        }
     1639}/*}}}*/
     1640bool       Tria::IsIcefront(void){/*{{{*/
     1641
     1642        bool isicefront;
     1643        int i,nrice;
     1644   IssmDouble ls[NUMVERTICES];
     1645
     1646        /*Retrieve all inputs and parameters*/
     1647        GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum);
     1648
     1649        /* If only one vertex has ice, there is an ice front here */
     1650        isicefront=false;
     1651        if(IsIceInElement()){
     1652                nrice=0;       
     1653                for(i=0;i<NUMVERTICES;i++)
     1654                        if(ls[i]<0.) nrice++;
     1655                if(nrice==1) isicefront= true;
     1656        }
     1657        return isicefront;
     1658}/*}}}*/
     1659bool       Tria::IsNodeOnShelfFromFlags(IssmDouble* flags){/*{{{*/
     1660
     1661        int  i;
     1662        bool shelf=false;
     1663
     1664        for(i=0;i<NUMVERTICES;i++){
     1665                if (flags[vertices[i]->Pid()]<0.){
     1666                        shelf=true;
     1667                        break;
     1668                }
     1669        }
     1670        return shelf;
     1671}
     1672/*}}}*/
    13731673bool       Tria::IsOnBase(){/*{{{*/
    13741674
    13751675        int domaintype;
    … …  
    13961696        }
    13971697}
    13981698/*}}}*/
     1699bool       Tria::IsZeroLevelset(int levelset_enum){/*{{{*/
     1700
     1701        bool iszerols;
     1702        IssmDouble ls[NUMVERTICES];
     1703
     1704        /*Retrieve all inputs and parameters*/
     1705        GetInputListOnVertices(&ls[0],levelset_enum);
     1706
     1707        /*If the level set is awlays <0, there is no ice front here*/
     1708        iszerols= false;
     1709        if(IsIceInElement()){
     1710                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.)){
     1711                        iszerols = true;
     1712                }
     1713        }
     1714
     1715        return iszerols;
     1716}
     1717/*}}}*/
    13991718void       Tria::JacobianDeterminant(IssmDouble* pJdet,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    14001719
    14011720        _assert_(gauss->Enum()==GaussTriaEnum);
    … …  
    14241743
    14251744}
    14261745/*}}}*/
    1427 bool       Tria::HasEdgeOnBase(){/*{{{*/
     1746IssmDouble Tria::Masscon(IssmDouble* levelset){ /*{{{*/
    14281747
    1429         IssmDouble values[NUMVERTICES];
    1430         IssmDouble sum;
    14311748
    1432         /*Retrieve all inputs and parameters*/
    1433         GetInputListOnVertices(&values[0],MeshVertexonbaseEnum);
    1434         sum = values[0]+values[1]+values[2];
     1749        /*intermediary: */
     1750        IssmDouble* values=NULL;
     1751        Input*      thickness_input=NULL;
     1752        IssmDouble  thickness;
     1753        IssmDouble  weight;
     1754        IssmDouble  Jdet;
     1755        IssmDouble  volume;
     1756        IssmDouble  rho_ice;
     1757        IssmDouble* xyz_list=NULL;
     1758        int         point1;
     1759        IssmDouble  fraction1,fraction2;
     1760        bool        mainlynegative=true;
     1761       
     1762        /*Output:*/
     1763        volume=0;
    14351764
    1436         _assert_(sum==0. || sum==1. || sum==2.);
     1765        /* Get node coordinates and dof list: */
     1766        GetVerticesCoordinates(&xyz_list);
    14371767
    1438         if(sum==3.)  _error_("Two edges on bed not supported yet...");
     1768        /*Retrieve inputs required:*/
     1769        thickness_input=this->GetInput(ThicknessEnum); _assert_(thickness_input);
     1770       
     1771        /*Retrieve material parameters: */
     1772        rho_ice=matpar->GetRhoIce();
    14391773
    1440         if(sum>1.){
    1441                 return true;
     1774        /*Retrieve values of the levelset defining the masscon: */
     1775        values = xNew<IssmDouble>(NUMVERTICES);
     1776        for(int i=0;i<NUMVERTICES;i++){
     1777                values[i]=levelset[this->vertices[i]->Sid()];
    14421778        }
    1443         else{
    1444                 return false;
     1779               
     1780        /*Ok, use the level set values to figure out where we put our gaussian points:*/
     1781        this->GetLevelsetPositivePart(&point1,&fraction1,&fraction2,&mainlynegative,values);
     1782        Gauss* gauss = this->NewGauss(point1,fraction1,fraction2,mainlynegative,4);
     1783
     1784        volume=0;
     1785
     1786        for(int ig=gauss->begin();ig<gauss->end();ig++){
     1787                gauss->GaussPoint(ig);
     1788
     1789                this->JacobianDeterminant(&Jdet,xyz_list,gauss);
     1790                thickness_input->GetInputValue(&thickness, gauss);
     1791
     1792                volume+=thickness*gauss->weight*Jdet;
    14451793        }
     1794
     1795        /* clean up and Return: */
     1796        xDelete<IssmDouble>(xyz_list);
     1797        xDelete<IssmDouble>(values);
     1798        delete gauss;
     1799        return rho_ice*volume;
    14461800}
    14471801/*}}}*/
    1448 bool       Tria::HasEdgeOnSurface(){/*{{{*/
     1802IssmDouble Tria::MassFlux(IssmDouble x1, IssmDouble y1, IssmDouble x2, IssmDouble y2,int segment_id){/*{{{*/
    14491803
    1450         IssmDouble values[NUMVERTICES];
    1451         IssmDouble sum;
     1804        int        domaintype;
     1805        IssmDouble mass_flux=0.;
     1806        IssmDouble xyz_list[NUMVERTICES][3];
     1807        IssmDouble normal[2];
     1808        IssmDouble length,rho_ice;
     1809        IssmDouble h1,h2;
     1810        IssmDouble vx1,vx2,vy1,vy2;
     1811        GaussTria* gauss_1=NULL;
     1812        GaussTria* gauss_2=NULL;
    14521813
    1453         /*Retrieve all inputs and parameters*/
    1454         GetInputListOnVertices(&values[0],MeshVertexonsurfaceEnum);
    1455         sum = values[0]+values[1]+values[2];
     1814        /*Get material parameters :*/
     1815        rho_ice=matpar->GetRhoIce();
    14561816
    1457         _assert_(sum==0. || sum==1. || sum==2.);
     1817        /*First off, check that this segment belongs to this element: */
     1818        if (segment_id!=this->id)_error_("error message: segment with id " << segment_id << " does not belong to element with id:" << this->id);
    14581819
    1459         if(sum==3.)  _error_("Two edges on surface not supported yet...");
     1820        /*Get xyz list: */
     1821        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    14601822
    1461         if(sum>1.){
    1462                 return true;
     1823        /*get area coordinates of 0 and 1 locations: */
     1824        gauss_1=new GaussTria();
     1825        gauss_1->GaussFromCoords(x1,y1,&xyz_list[0][0]);
     1826        gauss_2=new GaussTria();
     1827        gauss_2->GaussFromCoords(x2,y2,&xyz_list[0][0]);
     1828
     1829        normal[0]=cos(atan2(x1-x2,y2-y1));
     1830        normal[1]=sin(atan2(x1-x2,y2-y1));
     1831
     1832        length=sqrt(pow(x2-x1,2)+pow(y2-y1,2));
     1833
     1834        Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
     1835        this->parameters->FindParam(&domaintype,DomainTypeEnum);
     1836        Input* vx_input=NULL;
     1837        Input* vy_input=NULL;
     1838        if(domaintype==Domain2DhorizontalEnum){
     1839                vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
     1840                vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
    14631841        }
    14641842        else{
    1465                 return false;
     1843                vx_input=inputs->GetInput(VxAverageEnum); _assert_(vx_input);
     1844                vy_input=inputs->GetInput(VyAverageEnum); _assert_(vy_input);
    14661845        }
    1467 }
    1468 /*}}}*/
    1469 void       Tria::EdgeOnBaseIndices(int* pindex1,int* pindex2){/*{{{*/
    14701846
    1471         IssmDouble values[NUMVERTICES];
    1472         int        indices[3][2] = {{1,2},{2,0},{0,1}};
     1847        thickness_input->GetInputValue(&h1, gauss_1);
     1848        thickness_input->GetInputValue(&h2, gauss_2);
     1849        vx_input->GetInputValue(&vx1,gauss_1);
     1850        vx_input->GetInputValue(&vx2,gauss_2);
     1851        vy_input->GetInputValue(&vy1,gauss_1);
     1852        vy_input->GetInputValue(&vy2,gauss_2);
    14731853
    1474         /*Retrieve all inputs and parameters*/
    1475         GetInputListOnVertices(&values[0],MeshVertexonbaseEnum);
     1854        mass_flux= rho_ice*length*( 
     1855                                (ONETHIRD*(h1-h2)*(vx1-vx2)+0.5*h2*(vx1-vx2)+0.5*(h1-h2)*vx2+h2*vx2)*normal[0]+
     1856                                (ONETHIRD*(h1-h2)*(vy1-vy2)+0.5*h2*(vy1-vy2)+0.5*(h1-h2)*vy2+h2*vy2)*normal[1]
     1857                                );
    14761858
    1477         for(int i=0;i<3;i++){
    1478                 if(values[indices[i][0]] == 1. && values[indices[i][1]] == 1.){
    1479                         *pindex1 = indices[i][0];
    1480                         *pindex2 = indices[i][1];
    1481                         return;
    1482                 }
    1483         }
    1484 
    1485         _printf_("list of vertices on bed: "<<values[0]<<" "<<values[1]<<" "<<values[2]);
    1486         _error_("Could not find 2 vertices on bed");
     1859        /*clean up and return:*/
     1860        delete gauss_1;
     1861        delete gauss_2;
     1862        return mass_flux;
    14871863}
    14881864/*}}}*/
    1489 void       Tria::EdgeOnSurfaceIndices(int* pindex1,int* pindex2){/*{{{*/
     1865IssmDouble Tria::MassFlux(IssmDouble* segment){/*{{{*/
    14901866
    1491         IssmDouble values[NUMVERTICES];
    1492         int        indices[3][2] = {{1,2},{2,0},{0,1}};
     1867        int        domaintype;
     1868        IssmDouble mass_flux=0.;
     1869        IssmDouble xyz_list[NUMVERTICES][3];
     1870        IssmDouble normal[2];
     1871        IssmDouble length,rho_ice;
     1872        IssmDouble x1,y1,x2,y2,h1,h2;
     1873        IssmDouble vx1,vx2,vy1,vy2;
     1874        GaussTria* gauss_1=NULL;
     1875        GaussTria* gauss_2=NULL;
    14931876
    1494         /*Retrieve all inputs and parameters*/
    1495         GetInputListOnVertices(&values[0],MeshVertexonsurfaceEnum);
     1877        /*Get material parameters :*/
     1878        rho_ice=matpar->GetRhoIce();
    14961879
    1497         for(int i=0;i<3;i++){
    1498                 if(values[indices[i][0]] == 1. && values[indices[i][1]] == 1.){
    1499                         *pindex1 = indices[i][0];
    1500                         *pindex2 = indices[i][1];
    1501                         return;
    1502                 }
    1503         }
     1880        /*First off, check that this segment belongs to this element: */
     1881        if (reCast<int>(*(segment+4))!=this->id)_error_("error message: segment with id " << reCast<int>(*(segment+4)) << " does not belong to element with id:" << this->id);
    15041882
    1505         _printf_("list of vertices on surface: "<<values[0]<<" "<<values[1]<<" "<<values[2]);
    1506         _error_("Could not find 2 vertices on surface");
    1507 }
    1508 /*}}}*/
    1509 int        Tria::EdgeOnBaseIndex(void){/*{{{*/
     1883        /*Recover segment node locations: */
     1884        x1=*(segment+0); y1=*(segment+1); x2=*(segment+2); y2=*(segment+3);
    15101885
    1511         IssmDouble values[NUMVERTICES];
    1512         int        indices[3][2] = {{1,2},{2,0},{0,1}};
     1886        /*Get xyz list: */
     1887        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    15131888
    1514         /*Retrieve all inputs and parameters*/
    1515         GetInputListOnVertices(&values[0],MeshVertexonbaseEnum);
     1889        /*get area coordinates of 0 and 1 locations: */
     1890        gauss_1=new GaussTria();
     1891        gauss_1->GaussFromCoords(x1,y1,&xyz_list[0][0]);
     1892        gauss_2=new GaussTria();
     1893        gauss_2->GaussFromCoords(x2,y2,&xyz_list[0][0]);
    15161894
    1517         for(int i=0;i<3;i++){
    1518                 if(values[indices[i][0]] == 1. && values[indices[i][1]] == 1.){
    1519                         return i;
    1520                 }
     1895        normal[0]=cos(atan2(x1-x2,y2-y1));
     1896        normal[1]=sin(atan2(x1-x2,y2-y1));
     1897
     1898        length=sqrt(pow(x2-x1,2)+pow(y2-y1,2));
     1899
     1900        Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
     1901        this->parameters->FindParam(&domaintype,DomainTypeEnum);
     1902        Input* vx_input=NULL;
     1903        Input* vy_input=NULL;
     1904        if(domaintype==Domain2DhorizontalEnum){
     1905                vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
     1906                vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
    15211907        }
     1908        else{
     1909                vx_input=inputs->GetInput(VxAverageEnum); _assert_(vx_input);
     1910                vy_input=inputs->GetInput(VyAverageEnum); _assert_(vy_input);
     1911        }
    15221912
    1523         _printf_("list of vertices on bed: "<<values[0]<<" "<<values[1]<<" "<<values[2]);
    1524         _error_("Could not find 2 vertices on bed");
     1913        thickness_input->GetInputValue(&h1, gauss_1);
     1914        thickness_input->GetInputValue(&h2, gauss_2);
     1915        vx_input->GetInputValue(&vx1,gauss_1);
     1916        vx_input->GetInputValue(&vx2,gauss_2);
     1917        vy_input->GetInputValue(&vy1,gauss_1);
     1918        vy_input->GetInputValue(&vy2,gauss_2);
     1919
     1920        mass_flux= rho_ice*length*( 
     1921                                (ONETHIRD*(h1-h2)*(vx1-vx2)+0.5*h2*(vx1-vx2)+0.5*(h1-h2)*vx2+h2*vx2)*normal[0]+
     1922                                (ONETHIRD*(h1-h2)*(vy1-vy2)+0.5*h2*(vy1-vy2)+0.5*(h1-h2)*vy2+h2*vy2)*normal[1]
     1923                                );
     1924
     1925        /*clean up and return:*/
     1926        delete gauss_1;
     1927        delete gauss_2;
     1928        return mass_flux;
    15251929}
    15261930/*}}}*/
    1527 int        Tria::EdgeOnSurfaceIndex(void){/*{{{*/
     1931IssmDouble Tria::Misfit(int modelenum,int observationenum,int weightsenum){/*{{{*/
    15281932
    1529         IssmDouble values[NUMVERTICES];
    1530         int        indices[3][2] = {{1,2},{2,0},{0,1}};
     1933        /*Intermediaries*/
     1934        IssmDouble model,observation,weight;
     1935        IssmDouble Jdet;
     1936        IssmDouble Jelem = 0;
     1937        IssmDouble xyz_list[NUMVERTICES][3];
     1938        GaussTria *gauss = NULL;
    15311939
    1532         /*Retrieve all inputs and parameters*/
    1533         GetInputListOnVertices(&values[0],MeshVertexonsurfaceEnum);
     1940        /*If on water, return 0: */
     1941        if(!IsIceInElement())return 0;
    15341942
    1535         for(int i=0;i<3;i++){
    1536                 if(values[indices[i][0]] == 1. && values[indices[i][1]] == 1.){
    1537                         return i;
    1538                 }
    1539         }
     1943        /*Retrieve all inputs we will be needing: */
     1944        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     1945        Input* model_input=inputs->GetInput(modelenum);   _assert_(model_input);
     1946        Input* observation_input=inputs->GetInput(observationenum);_assert_(observation_input);
     1947        Input* weights_input     =inputs->GetInput(weightsenum);     _assert_(weights_input);
    15401948
    1541         _printf_("list of vertices on surface: "<<values[0]<<" "<<values[1]<<" "<<values[2]);
    1542         _error_("Could not find 2 vertices on surface");
    1543 }
    1544 /*}}}*/
    1545 void       Tria::FSContactMigration(Vector<IssmDouble>* vertexgrounded,Vector<IssmDouble>* vertexfloating){/*{{{*/
     1949        /* Start  looping on the number of gaussian points: */
     1950        gauss=new GaussTria(2);
     1951        for(int ig=gauss->begin();ig<gauss->end();ig++){
    15461952
    1547         if(!IsOnBase()) return;
     1953                gauss->GaussPoint(ig);
    15481954
    1549         int approximation;
    1550         inputs->GetInputValue(&approximation,ApproximationEnum);
     1955                /* Get Jacobian determinant: */
     1956                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    15511957
    1552         if(approximation==HOApproximationEnum || approximation==SSAApproximationEnum || approximation==SSAHOApproximationEnum){
    1553                 for(int i=0;i<NUMVERTICES;i++){
    1554                         vertexgrounded->SetValue(vertices[i]->Pid(),+9999.,INS_VAL);
    1555                         vertexfloating->SetValue(vertices[i]->Pid(),+9999.,INS_VAL);
    1556                 }
     1958                /*Get parameters at gauss point*/
     1959                model_input->GetInputValue(&model,gauss);
     1960                observation_input->GetInputValue(&observation,gauss);
     1961                weights_input->GetInputValue(&weight,gauss);
     1962
     1963                /*compute misfit between model and observation */
     1964                Jelem+=0.5*(model-observation)*(model-observation)*Jdet*weight*gauss->weight;
    15571965        }
    1558         else{
    1559                 /*Intermediaries*/
    1560                 IssmDouble* xyz_list = NULL;
    1561                 IssmDouble* xyz_list_base = NULL;
    1562                 IssmDouble  pressure,water_pressure,sigma_nn,viscosity,bed,base;
    1563                 IssmDouble  bed_normal[2];
    1564                 IssmDouble  epsilon[3]; /* epsilon=[exx,eyy,exy];*/
    1565                 IssmDouble  surface=0,value=0;
    1566                 bool grounded;
    15671966
    1568                 /* Get node coordinates and dof list: */
    1569                 GetVerticesCoordinates(&xyz_list);
    1570                 GetVerticesCoordinatesBase(&xyz_list_base);
     1967        /* clean up and Return: */
     1968        delete gauss;
     1969        return Jelem;
     1970}
     1971/*}}}*/
     1972IssmDouble Tria::MisfitArea(int weightsenum){/*{{{*/
    15711973
    1572                 /*Retrieve all inputs we will be needing: */
    1573                 Input* pressure_input = inputs->GetInput(PressureEnum); _assert_(pressure_input);
    1574                 Input* base_input     = inputs->GetInput(BaseEnum);     _assert_(base_input);
    1575                 Input* bed_input      = inputs->GetInput(BedEnum);      _assert_(bed_input);
    1576                 Input* vx_input       = inputs->GetInput(VxEnum);       _assert_(vx_input);
    1577                 Input* vy_input       = inputs->GetInput(VyEnum);       _assert_(vy_input);
     1974        /*Intermediaries*/
     1975        IssmDouble weight;
     1976        IssmDouble Jdet;
     1977        IssmDouble Jelem = 0;
     1978        IssmDouble xyz_list[NUMVERTICES][3];
     1979        GaussTria *gauss = NULL;
    15781980
    1579                 /*Create gauss point in the middle of the basal edge*/
    1580                 Gauss* gauss=NewGaussBase(1);
    1581                 gauss->GaussPoint(0);
     1981        /*If on water, return 0: */
     1982        if(!IsIceInElement())return 0;
    15821983
    1583                 if(!IsFloating()){
    1584                         /*Check for basal force only if grounded and touching GL*/
    1585                         //              if(this->inputs->Min(MaskGroundediceLevelsetEnum)==0.){
    1586                         this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
    1587                         this->ViscosityFS(&viscosity,2,xyz_list,gauss,vx_input,vy_input,NULL);
    1588                         pressure_input->GetInputValue(&pressure, gauss);
    1589                         base_input->GetInputValue(&base, gauss);
     1984        /*Retrieve all inputs we will be needing: */
     1985        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     1986        Input* weights_input     =inputs->GetInput(weightsenum);     _assert_(weights_input);
    15901987
    1591                         /*Compute Stress*/
    1592                         IssmDouble sigma_xx=2.*viscosity*epsilon[0]-pressure;
    1593                         IssmDouble sigma_yy=2.*viscosity*epsilon[1]-pressure;
    1594                         IssmDouble sigma_xy=2.*viscosity*epsilon[2];
     1988        /* Start  looping on the number of gaussian points: */
     1989        gauss=new GaussTria(2);
     1990        for(int ig=gauss->begin();ig<gauss->end();ig++){
    15951991
    1596                         /*Get normal vector to the bed */
    1597                         NormalBase(&bed_normal[0],xyz_list_base);
     1992                gauss->GaussPoint(ig);
    15981993
    1599                         /*basalforce*/
    1600                         sigma_nn = sigma_xx*bed_normal[0]*bed_normal[0] + sigma_yy*bed_normal[1]*bed_normal[1] + 2.*sigma_xy*bed_normal[0]*bed_normal[1];
     1994                /* Get Jacobian determinant: */
     1995                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    16011996
    1602                         /*Compute water pressure*/
    1603                         IssmDouble rho_ice   = matpar->GetRhoIce();
    1604                         IssmDouble rho_water = matpar->GetRhoWater();
    1605                         IssmDouble gravity   = matpar->GetG();
    1606                         water_pressure=gravity*rho_water*base;
     1997                /*Get parameters at gauss point*/
     1998                weights_input->GetInputValue(&weight,gauss);
    16071999
    1608                         /*Compare basal stress to water pressure and determine whether it should ground*/
    1609                         if (sigma_nn<water_pressure) grounded=true;
    1610                         else                         grounded=false;
    1611                 }
    1612                 else{
    1613                         /*Check for basal elevation if floating*/
    1614                         base_input->GetInputValue(&base, gauss);
    1615                         bed_input->GetInputValue(&bed, gauss);
    1616                         if(base<bed) grounded=true;
    1617                         else         grounded=false;
    1618                 }
    1619                 for(int i=0;i<NUMVERTICES;i++){
    1620                         if(grounded) vertexgrounded->SetValue(vertices[i]->Pid(),+1.,INS_VAL);
    1621                         else         vertexfloating->SetValue(vertices[i]->Pid(),+1.,INS_VAL);
    1622                 }
    1623 
    1624                 /*clean up*/
    1625                 delete gauss;
    1626                 xDelete<IssmDouble>(xyz_list);
    1627                 xDelete<IssmDouble>(xyz_list_base);
     2000                /*compute misfit between model and observation */
     2001                Jelem+=Jdet*weight*gauss->weight;
    16282002        }
    1629 }
    1630 /*}}}*/
    1631 bool       Tria::IsNodeOnShelfFromFlags(IssmDouble* flags){/*{{{*/
    16322003
    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;
     2004        /* clean up and Return: */
     2005        delete gauss;
     2006        return Jelem;
    16432007}
    16442008/*}}}*/
    16452009Gauss*     Tria::NewGauss(void){/*{{{*/
    … …  
    16902054
    16912055}
    16922056/*}}}*/
    1693 void       Tria::NodalFunctionsP1(IssmDouble* basis, Gauss* gauss){/*{{{*/
     2057void       Tria::NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    16942058
    16952059        _assert_(gauss->Enum()==GaussTriaEnum);
    1696         this->GetNodalFunctions(basis,(GaussTria*)gauss,P1Enum);
     2060        this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussTria*)gauss,this->element_type);
    16972061
    16982062}
    16992063/*}}}*/
    1700 void       Tria::NodalFunctionsP2(IssmDouble* basis, Gauss* gauss){/*{{{*/
     2064void       Tria::NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    17012065
    17022066        _assert_(gauss->Enum()==GaussTriaEnum);
    1703         this->GetNodalFunctions(basis,(GaussTria*)gauss,P2Enum);
     2067        this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussTria*)gauss,this->VelocityInterpolation());
    17042068
    17052069}
    17062070/*}}}*/
    1707 void       Tria::NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     2071void       Tria::NodalFunctionsPressure(IssmDouble* basis, Gauss* gauss){/*{{{*/
    17082072
    17092073        _assert_(gauss->Enum()==GaussTriaEnum);
    1710         this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussTria*)gauss,this->element_type);
     2074        this->GetNodalFunctions(basis,(GaussTria*)gauss,this->PressureInterpolation());
    17112075
    17122076}
    17132077/*}}}*/
    1714 void       Tria::NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     2078void       Tria::NodalFunctionsP1(IssmDouble* basis, Gauss* gauss){/*{{{*/
    17152079
    17162080        _assert_(gauss->Enum()==GaussTriaEnum);
    1717         this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussTria*)gauss,P1Enum);
     2081        this->GetNodalFunctions(basis,(GaussTria*)gauss,P1Enum);
    17182082
    17192083}
    17202084/*}}}*/
    1721 void       Tria::NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     2085void       Tria::NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    17222086
    17232087        _assert_(gauss->Enum()==GaussTriaEnum);
    1724         this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussTria*)gauss,this->VelocityInterpolation());
     2088        this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussTria*)gauss,P1Enum);
    17252089
    17262090}
    17272091/*}}}*/
    1728 void       Tria::NodalFunctionsVelocity(IssmDouble* basis, Gauss* gauss){/*{{{*/
     2092void       Tria::NodalFunctionsP2(IssmDouble* basis, Gauss* gauss){/*{{{*/
    17292093
    17302094        _assert_(gauss->Enum()==GaussTriaEnum);
    1731         this->GetNodalFunctions(basis,(GaussTria*)gauss,this->VelocityInterpolation());
     2095        this->GetNodalFunctions(basis,(GaussTria*)gauss,P2Enum);
    17322096
    17332097}
    17342098/*}}}*/
    1735 void       Tria::NodalFunctionsPressure(IssmDouble* basis, Gauss* gauss){/*{{{*/
     2099void       Tria::NodalFunctionsTensor(IssmDouble* basis, Gauss* gauss){/*{{{*/
    17362100
    17372101        _assert_(gauss->Enum()==GaussTriaEnum);
    1738         this->GetNodalFunctions(basis,(GaussTria*)gauss,this->PressureInterpolation());
     2102        this->GetNodalFunctions(basis,(GaussTria*)gauss,this->TensorInterpolation());
    17392103
    17402104}
    17412105/*}}}*/
    1742 void       Tria::NodalFunctionsTensor(IssmDouble* basis, Gauss* gauss){/*{{{*/
     2106void       Tria::NodalFunctionsVelocity(IssmDouble* basis, Gauss* gauss){/*{{{*/
    17432107
    17442108        _assert_(gauss->Enum()==GaussTriaEnum);
    1745         this->GetNodalFunctions(basis,(GaussTria*)gauss,this->TensorInterpolation());
     2109        this->GetNodalFunctions(basis,(GaussTria*)gauss,this->VelocityInterpolation());
    17462110
    17472111}
    17482112/*}}}*/
    … …  
    17942158        _assert_(bed_normal[1]<0);
    17952159}
    17962160/*}}}*/
     2161void       Tria::NormalSection(IssmDouble* normal,IssmDouble* xyz_list){/*{{{*/
     2162
     2163        /*Build unit outward pointing vector*/
     2164        IssmDouble vector[2];
     2165        IssmDouble norm;
     2166
     2167        vector[0]=xyz_list[1*3+0] - xyz_list[0*3+0];
     2168        vector[1]=xyz_list[1*3+1] - xyz_list[0*3+1];
     2169
     2170        norm=sqrt(vector[0]*vector[0] + vector[1]*vector[1]);
     2171
     2172        normal[0]= + vector[1]/norm;
     2173        normal[1]= - vector[0]/norm;
     2174}
     2175/*}}}*/
    17972176void       Tria::NormalTop(IssmDouble* top_normal,IssmDouble* xyz_list){/*{{{*/
    17982177
    17992178        /*Build unit outward pointing vector*/
    … …  
    18122191        _assert_(top_normal[1]>0);
    18132192}
    18142193/*}}}*/
    1815 int        Tria::VelocityInterpolation(void){/*{{{*/
    1816         return TriaRef::VelocityInterpolation(this->element_type);
     2194int        Tria::ObjectEnum(void){/*{{{*/
     2195
     2196        return TriaEnum;
     2197
    18172198}
    18182199/*}}}*/
    18192200int        Tria::PressureInterpolation(void){/*{{{*/
    18202201        return TriaRef::PressureInterpolation(this->element_type);
    18212202}
    18222203/*}}}*/
    1823 int        Tria::TensorInterpolation(void){/*{{{*/
    1824         return TriaRef::TensorInterpolation(this->element_type);
    1825 }
    1826 /*}}}*/
    18272204int        Tria::NumberofNodesPressure(void){/*{{{*/
    18282205        return TriaRef::NumberofNodes(this->PressureInterpolation());
    18292206}
    … …  
    19002277        delete gauss;
    19012278}
    19022279/*}}}*/
     2280void       Tria::PotentialUngrounding(Vector<IssmDouble>* potential_ungrounding){/*{{{*/
     2281
     2282        IssmDouble  h[NUMVERTICES],r[NUMVERTICES],gl[NUMVERTICES];
     2283        IssmDouble  bed_hydro;
     2284        IssmDouble  rho_water,rho_ice,density;
     2285
     2286        /*material parameters: */
     2287        rho_water=matpar->GetRhoWater();
     2288        rho_ice=matpar->GetRhoIce();
     2289        density=rho_ice/rho_water;
     2290        GetInputListOnVertices(&h[0],ThicknessEnum);
     2291        GetInputListOnVertices(&r[0],BedEnum);
     2292        GetInputListOnVertices(&gl[0],MaskGroundediceLevelsetEnum);
     2293
     2294        /*go through vertices, and figure out which ones are grounded and want to unground: */
     2295        for(int i=0;i<NUMVERTICES;i++){
     2296                /*Find if grounded vertices want to start floating*/
     2297                if (gl[i]>0.){
     2298                        bed_hydro=-density*h[i];
     2299                        if(bed_hydro>r[i]){
     2300                                /*Vertex that could potentially unground, flag it*/
     2301                                potential_ungrounding->SetValue(vertices[i]->Pid(),1,INS_VAL);
     2302                        }
     2303                }
     2304        }
     2305}
     2306/*}}}*/
    19032307void       Tria::ReduceMatrices(ElementMatrix* Ke,ElementVector* pe){/*{{{*/
    19042308
    19052309        /*Static condensation if requested*/
    … …  
    20102414        _error_("not implemented yet");
    20112415}
    20122416/*}}}*/
     2417void       Tria::SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index){/*{{{*/
     2418
     2419        IssmDouble  values[NUMVERTICES];
     2420        int         vertexpidlist[NUMVERTICES],control_init;
     2421
     2422
     2423        /*Get Domain type*/
     2424        int domaintype;
     2425        parameters->FindParam(&domaintype,DomainTypeEnum);
     2426
     2427        /*Specific case for depth averaged quantities*/
     2428        control_init=control_enum;
     2429        if(domaintype==Domain2DverticalEnum){
     2430                if(control_enum==MaterialsRheologyBbarEnum){
     2431                        control_enum=MaterialsRheologyBEnum;
     2432                        if(!IsOnBase()) return;
     2433                }
     2434                if(control_enum==DamageDbarEnum){
     2435                        control_enum=DamageDEnum;
     2436                        if(!IsOnBase()) return;
     2437                }
     2438        }
     2439
     2440        /*Get out if this is not an element input*/
     2441        if(!IsInput(control_enum)) return;
     2442
     2443        /*Prepare index list*/
     2444        GradientIndexing(&vertexpidlist[0],control_index);
     2445
     2446        /*Get values on vertices*/
     2447        for(int i=0;i<NUMVERTICES;i++){
     2448                values[i]=vector[vertexpidlist[i]];
     2449        }
     2450        Input* new_input = new TriaInput(control_enum,values,P1Enum);
     2451        Input* input     = (Input*)this->inputs->GetInput(control_enum);   _assert_(input);
     2452        if(input->ObjectEnum()!=ControlInputEnum){
     2453                _error_("input " << EnumToStringx(control_enum) << " is not a ControlInput");
     2454        }
     2455
     2456        ((ControlInput*)input)->SetInput(new_input);
     2457}
     2458/*}}}*/
     2459void       Tria::SetCurrentConfiguration(Elements* elementsin, Loads* loadsin, Nodes* nodesin, Materials* materialsin, Parameters* parametersin){/*{{{*/
     2460
     2461        /*go into parameters and get the analysis_counter: */
     2462        int analysis_counter;
     2463        parametersin->FindParam(&analysis_counter,AnalysisCounterEnum);
     2464
     2465        /*Get Element type*/
     2466        if(this->element_type_list) this->element_type=this->element_type_list[analysis_counter];
     2467
     2468        /*Pick up nodes*/
     2469        if(this->hnodes && this->hnodes[analysis_counter]){
     2470                this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp();
     2471        }
     2472
     2473}
     2474/*}}}*/
     2475Element*   Tria::SpawnBasalElement(void){/*{{{*/
     2476
     2477        int index1,index2;
     2478        int domaintype;
     2479
     2480        this->parameters->FindParam(&domaintype,DomainTypeEnum);
     2481        switch(domaintype){
     2482                case Domain2DhorizontalEnum:
     2483                        return this;
     2484                case Domain2DverticalEnum:
     2485                        _assert_(HasEdgeOnBase());
     2486                        this->EdgeOnBaseIndices(&index1,&index2);
     2487                        return SpawnSeg(index1,index2);
     2488                default:
     2489                        _error_("not implemented yet");
     2490        }
     2491}
     2492/*}}}*/
    20132493Seg*       Tria::SpawnSeg(int index1,int index2){/*{{{*/
    20142494
    20152495        int analysis_counter;
    … …  
    20372517        return seg;
    20382518}
    20392519/*}}}*/
    2040 Element*   Tria::SpawnBasalElement(void){/*{{{*/
     2520Element*   Tria::SpawnTopElement(void){/*{{{*/
    20412521
    20422522        int index1,index2;
    20432523        int domaintype;
    … …  
    20472527                case Domain2DhorizontalEnum:
    20482528                        return this;
    20492529                case Domain2DverticalEnum:
    2050                         _assert_(HasEdgeOnBase());
    2051                         this->EdgeOnBaseIndices(&index1,&index2);
    2052                         return SpawnSeg(index1,index2);
     2530                        _assert_(HasEdgeOnSurface());
     2531                        this->EdgeOnSurfaceIndices(&index1,&index2);
     2532                        return SpawnSeg(index2,index1); //reverse order
    20532533                default:
    20542534                        _error_("not implemented yet");
    20552535        }
    20562536}
    20572537/*}}}*/
    2058 Element*   Tria::SpawnTopElement(void){/*{{{*/
     2538void       Tria::StrainRateparallel(){/*{{{*/
    20592539
    2060         int index1,index2;
    2061         int domaintype;
     2540        IssmDouble *xyz_list = NULL;
     2541        IssmDouble  epsilon[3];
     2542        GaussTria* gauss=NULL;
     2543        IssmDouble  vx,vy,vel;
     2544        IssmDouble  strainxx;
     2545        IssmDouble  strainxy;
     2546        IssmDouble  strainyy;
     2547        IssmDouble  strainparallel[NUMVERTICES];
    20622548
    2063         this->parameters->FindParam(&domaintype,DomainTypeEnum);
    2064         switch(domaintype){
    2065                 case Domain2DhorizontalEnum:
    2066                         return this;
    2067                 case Domain2DverticalEnum:
    2068                         _assert_(HasEdgeOnSurface());
    2069                         this->EdgeOnSurfaceIndices(&index1,&index2);
    2070                         return SpawnSeg(index2,index1); //reverse order
    2071                 default:
    2072                         _error_("not implemented yet");
     2549        /* Get node coordinates and dof list: */
     2550        this->GetVerticesCoordinates(&xyz_list);
     2551
     2552        /*Retrieve all inputs we will need*/
     2553        Input* vx_input=inputs->GetInput(VxEnum);                                  _assert_(vx_input);
     2554        Input* vy_input=inputs->GetInput(VyEnum);                                  _assert_(vy_input);
     2555
     2556        /* Start looping on the number of vertices: */
     2557        gauss=new GaussTria();
     2558        for (int iv=0;iv<NUMVERTICES;iv++){
     2559                gauss->GaussVertex(iv);
     2560
     2561                /* Get the value we need*/
     2562                vx_input->GetInputValue(&vx,gauss);
     2563                vy_input->GetInputValue(&vy,gauss);
     2564                vel=vx*vx+vy*vy;
     2565
     2566                /*Compute strain rate viscosity and pressure: */
     2567                this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
     2568                strainxx=epsilon[0];
     2569                strainyy=epsilon[1];
     2570                strainxy=epsilon[2];
     2571
     2572                /*strainparallel= Strain rate along the ice flow direction */
     2573                strainparallel[iv]=(vx*vx*(strainxx)+vy*vy*(strainyy)+2*vy*vx*strainxy)/(vel+1.e-6);
    20732574        }
     2575
     2576        /*Add input*/
     2577        this->inputs->AddInput(new TriaInput(StrainRateparallelEnum,&strainparallel[0],P1Enum));
     2578
     2579        /*Clean up and return*/
     2580        delete gauss;
     2581        xDelete<IssmDouble>(xyz_list);
    20742582}
    20752583/*}}}*/
    2076 void       Tria::SetCurrentConfiguration(Elements* elementsin, Loads* loadsin, Nodes* nodesin, Materials* materialsin, Parameters* parametersin){/*{{{*/
     2584void       Tria::StrainRateperpendicular(){/*{{{*/
    20772585
    2078         /*go into parameters and get the analysis_counter: */
    2079         int analysis_counter;
    2080         parametersin->FindParam(&analysis_counter,AnalysisCounterEnum);
     2586        IssmDouble *xyz_list = NULL;
     2587        GaussTria* gauss=NULL;
     2588        IssmDouble  epsilon[3];
     2589        IssmDouble  vx,vy,vel;
     2590        IssmDouble  strainxx;
     2591        IssmDouble  strainxy;
     2592        IssmDouble  strainyy;
     2593        IssmDouble  strainperpendicular[NUMVERTICES];
    20812594
    2082         /*Get Element type*/
    2083         if(this->element_type_list) this->element_type=this->element_type_list[analysis_counter];
     2595        /* Get node coordinates and dof list: */
     2596        this->GetVerticesCoordinates(&xyz_list);
    20842597
    2085         /*Pick up nodes*/
    2086         if(this->hnodes && this->hnodes[analysis_counter]){
    2087                 this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp();
     2598        /*Retrieve all inputs we will need*/
     2599        Input* vx_input=inputs->GetInput(VxEnum);                                  _assert_(vx_input);
     2600        Input* vy_input=inputs->GetInput(VyEnum);                                  _assert_(vy_input);
     2601
     2602        /* Start looping on the number of vertices: */
     2603        gauss=new GaussTria();
     2604        for (int iv=0;iv<NUMVERTICES;iv++){
     2605                gauss->GaussVertex(iv);
     2606
     2607                /* Get the value we need*/
     2608                vx_input->GetInputValue(&vx,gauss);
     2609                vy_input->GetInputValue(&vy,gauss);
     2610                vel=vx*vx+vy*vy;
     2611
     2612                /*Compute strain rate viscosity and pressure: */
     2613                this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
     2614                strainxx=epsilon[0];
     2615                strainyy=epsilon[1];
     2616                strainxy=epsilon[2];
     2617
     2618                /*strainperpendicular= Strain rate perpendicular to the ice flow direction */
     2619                strainperpendicular[iv]=(vx*vx*(strainyy)+vy*vy*(strainxx)-2*vy*vx*strainxy)/(vel+1.e-6);
    20882620        }
    20892621
     2622        /*Add input*/
     2623        this->inputs->AddInput(new TriaInput(StrainRateperpendicularEnum,&strainperpendicular[0],P1Enum));
     2624
     2625        /*Clean up and return*/
     2626        delete gauss;
     2627        xDelete<IssmDouble>(xyz_list);
    20902628}
    20912629/*}}}*/
    20922630IssmDouble Tria::SurfaceArea(void){/*{{{*/
    … …  
    21162654        return S;
    21172655}
    21182656/*}}}*/
     2657int        Tria::TensorInterpolation(void){/*{{{*/
     2658        return TriaRef::TensorInterpolation(this->element_type);
     2659}
     2660/*}}}*/
    21192661IssmDouble Tria::TimeAdapt(void){/*{{{*/
    21202662
    21212663        /*intermediary: */
    … …  
    21572699        return dt;
    21582700}
    21592701/*}}}*/
     2702IssmDouble Tria::TotalSmb(void){/*{{{*/
     2703
     2704        /*The smb[kg yr-1] of one element is area[m2] * smb [kg m^-2 yr^-1]*/
     2705        IssmDouble base,smb,rho_ice;
     2706        IssmDouble Total_Smb=0;
     2707        IssmDouble xyz_list[NUMVERTICES][3];
     2708
     2709        /*Get material parameters :*/
     2710        rho_ice=matpar->GetRhoIce();
     2711
     2712   if(!IsIceInElement())return 0;
     2713
     2714        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     2715
     2716        /*First calculate the area of the base (cross section triangle)
     2717         * http://en.wikipedia.org/wiki/Triangle
     2718         * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
     2719        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])); // area of element in m2
     2720
     2721        /*Now get the average SMB over the element*/
     2722        Input* smb_input = inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(smb_input);
     2723        smb_input->GetInputAverage(&smb);                                                                                                                                                                                               // average smb on element in m ice s-1
     2724   Total_Smb=rho_ice*base*smb;                                                                                                                                                                                                                  // smb on element in kg s-1
     2725
     2726        /*Return: */
     2727        return Total_Smb;
     2728}
     2729/*}}}*/
    21602730void       Tria::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finiteelement_type){/*{{{*/
    21612731
    21622732        /*Intermediaries*/
    … …  
    23462916
    23472917}
    23482918/*}}}*/
    2349 void       Tria::ValueP1OnGauss(IssmDouble* pvalue,IssmDouble* values,Gauss* gauss){/*{{{*/
    2350         TriaRef::GetInputValue(pvalue,values,gauss,P1Enum);
    2351 }
    2352 /*}}}*/
    2353 void       Tria::ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    2354         TriaRef::GetInputDerivativeValue(dvalue,values,xyz_list,gauss,P1Enum);
    2355 }
    2356 /*}}}*/
    2357 int        Tria::VertexConnectivity(int vertexindex){/*{{{*/
    2358         _assert_(this->vertices);
    2359         return this->vertices[vertexindex]->Connectivity();
    2360 }
    2361 /*}}}*/
    2362 bool       Tria::IsZeroLevelset(int levelset_enum){/*{{{*/
     2919int        Tria::UpdatePotentialUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf){/*{{{*/
    23632920
    2364         bool iszerols;
    2365         IssmDouble ls[NUMVERTICES];
     2921        int i;
     2922        int nflipped=0;
    23662923
    2367         /*Retrieve all inputs and parameters*/
    2368         GetInputListOnVertices(&ls[0],levelset_enum);
     2924        /*Go through nodes, and whoever is on the potential_ungrounding, ends up in nodes_on_iceshelf: */
     2925        for(i=0;i<3;i++){
     2926                if (reCast<bool>(vertices_potentially_ungrounding[vertices[i]->Pid()])){
     2927                        vec_nodes_on_iceshelf->SetValue(vertices[i]->Pid(),-1.,INS_VAL);
    23692928
    2370         /*If the level set is awlays <0, there is no ice front here*/
    2371         iszerols= false;
    2372         if(IsIceInElement()){
    2373                 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.)){
    2374                         iszerols = true;
     2929                        /*If node was not on ice shelf, we flipped*/
     2930                        if(nodes_on_iceshelf[vertices[i]->Pid()]>=0.){
     2931                                nflipped++;
     2932                        }
    23752933                }
    23762934        }
    2377 
    2378         return iszerols;
     2935        return nflipped;
    23792936}
    23802937/*}}}*/
    2381 bool       Tria::IsIcefront(void){/*{{{*/
    2382 
    2383         bool isicefront;
    2384         int i,nrice;
    2385    IssmDouble ls[NUMVERTICES];
    2386 
    2387         /*Retrieve all inputs and parameters*/
    2388         GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum);
    2389 
    2390         /* If only one vertex has ice, there is an ice front here */
    2391         isicefront=false;
    2392         if(IsIceInElement()){
    2393                 nrice=0;       
    2394                 for(i=0;i<NUMVERTICES;i++)
    2395                         if(ls[i]<0.) nrice++;
    2396                 if(nrice==1) isicefront= true;
    2397         }
    2398         return isicefront;
    2399 }/*}}}*/
    2400 bool       Tria::IsFaceOnBoundary(void){/*{{{*/
    2401 
    2402         IssmDouble values[NUMVERTICES];
    2403         IssmDouble sum;
    2404 
    2405         /*Retrieve all inputs and parameters*/
    2406         GetInputListOnVertices(&values[0],MeshVertexonboundaryEnum);
    2407         sum = values[0]+values[1]+values[2];
    2408 
    2409         _assert_(sum==0. || sum==1. || sum==2.);
    2410 
    2411         if(sum==3.)  _error_("Two edges on boundary not supported yet...");
    2412 
    2413         if(sum>1.){
    2414                 return true;
    2415         }
    2416         else{
    2417                 return false;
    2418         }
    2419 }/*}}}*/
    2420 void       Tria::AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part){/*{{{*/
    2421 
    2422         bool       already = false;
    2423         int        i,j;
    2424         int        partition[NUMVERTICES];
    2425         int        offsetsid[NUMVERTICES];
    2426         int        offsetdof[NUMVERTICES];
    2427         IssmDouble area;
    2428         IssmDouble mean;
    2429 
    2430         /*First, get the area: */
    2431         area=this->GetArea();
    2432 
    2433         /*Figure out the average for this element: */
    2434         this->GetVerticesSidList(&offsetsid[0]);
    2435         this->GetVertexPidList(&offsetdof[0]);
    2436         mean=0;
    2437         for(i=0;i<NUMVERTICES;i++){
    2438                 partition[i]=reCast<int>(qmu_part[offsetsid[i]]);
    2439                 mean=mean+1.0/NUMVERTICES*vertex_response[offsetdof[i]];
    2440         }
    2441 
    2442         /*Add contribution: */
    2443         for(i=0;i<NUMVERTICES;i++){
    2444                 already=false;
    2445                 for(j=0;j<i;j++){
    2446                         if (partition[i]==partition[j]){
    2447                                 already=true;
    2448                                 break;
    2449                         }
    2450                 }
    2451                 if(!already){
    2452                         partition_contributions->SetValue(partition[i],mean*area,ADD_VAL);
    2453                         partition_areas->SetValue(partition[i],area,ADD_VAL);
    2454                 };
    2455         }
     2938void       Tria::ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     2939        TriaRef::GetInputDerivativeValue(dvalue,values,xyz_list,gauss,P1Enum);
    24562940}
    24572941/*}}}*/
    2458 IssmDouble Tria::IceVolume(void){/*{{{*/
    2459 
    2460         /*The volume of a troncated prism is base * 1/3 sum(length of edges)*/
    2461         IssmDouble base,surface,bed;
    2462         IssmDouble xyz_list[NUMVERTICES][3];
    2463 
    2464         if(!IsIceInElement())return 0;
    2465 
    2466         /*First get back the area of the base*/
    2467         base=this->GetArea();
    2468 
    2469         /*Now get the average height*/
    2470         Input* surface_input = inputs->GetInput(SurfaceEnum); _assert_(surface_input);
    2471         Input* base_input     = inputs->GetInput(BaseEnum);     _assert_(base_input);
    2472         surface_input->GetInputAverage(&surface);
    2473         base_input->GetInputAverage(&bed);
    2474 
    2475         /*Return: */
    2476         int domaintype;
    2477         parameters->FindParam(&domaintype,DomainTypeEnum);
    2478         if(domaintype==Domain2DverticalEnum){
    2479           return base;
    2480         }
    2481         else{
    2482           return base*(surface-bed);
    2483         }
     2942void       Tria::ValueP1OnGauss(IssmDouble* pvalue,IssmDouble* values,Gauss* gauss){/*{{{*/
     2943        TriaRef::GetInputValue(pvalue,values,gauss,P1Enum);
    24842944}
    24852945/*}}}*/
    2486 IssmDouble Tria::IceVolumeAboveFloatation(void){/*{{{*/
    2487 
    2488         /*The volume above floatation: H + rho_water/rho_ice * bathymetry */
    2489         IssmDouble rho_ice,rho_water;
    2490         IssmDouble base,surface,bed,bathymetry;
    2491         IssmDouble xyz_list[NUMVERTICES][3];
    2492 
    2493         if(!IsIceInElement() || IsFloating())return 0;
    2494 
    2495         rho_ice=matpar->GetRhoIce();
    2496         rho_water=matpar->GetRhoWater();
    2497         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    2498 
    2499         /*First calculate the area of the base (cross section triangle)
    2500          * http://en.wikipedia.org/wiki/Triangle
    2501          * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
    2502         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]));
    2503 
    2504         /*Now get the average height and bathymetry*/
    2505         Input* surface_input    = inputs->GetInput(SurfaceEnum);    _assert_(surface_input);
    2506         Input* base_input        = inputs->GetInput(BaseEnum);        _assert_(base_input);
    2507         Input* bed_input = inputs->GetInput(BedEnum); _assert_(bed_input);
    2508         surface_input->GetInputAverage(&surface);
    2509         base_input->GetInputAverage(&bed);
    2510         bed_input->GetInputAverage(&bathymetry);
    2511        
    2512         /*Return: */
    2513         return base*(surface-bed+min(rho_water/rho_ice*bathymetry,0.));
     2946int        Tria::VelocityInterpolation(void){/*{{{*/
     2947        return TriaRef::VelocityInterpolation(this->element_type);
    25142948}
    25152949/*}}}*/
    2516 IssmDouble Tria::MassFlux( IssmDouble x1, IssmDouble y1, IssmDouble x2, IssmDouble y2,int segment_id){/*{{{*/
    2517 
    2518         int        domaintype;
    2519         IssmDouble mass_flux=0.;
    2520         IssmDouble xyz_list[NUMVERTICES][3];
    2521         IssmDouble normal[2];
    2522         IssmDouble length,rho_ice;
    2523         IssmDouble h1,h2;
    2524         IssmDouble vx1,vx2,vy1,vy2;
    2525         GaussTria* gauss_1=NULL;
    2526         GaussTria* gauss_2=NULL;
    2527 
    2528         /*Get material parameters :*/
    2529         rho_ice=matpar->GetRhoIce();
    2530 
    2531         /*First off, check that this segment belongs to this element: */
    2532         if (segment_id!=this->id)_error_("error message: segment with id " << segment_id << " does not belong to element with id:" << this->id);
    2533 
    2534         /*Get xyz list: */
    2535         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    2536 
    2537         /*get area coordinates of 0 and 1 locations: */
    2538         gauss_1=new GaussTria();
    2539         gauss_1->GaussFromCoords(x1,y1,&xyz_list[0][0]);
    2540         gauss_2=new GaussTria();
    2541         gauss_2->GaussFromCoords(x2,y2,&xyz_list[0][0]);
    2542 
    2543         normal[0]=cos(atan2(x1-x2,y2-y1));
    2544         normal[1]=sin(atan2(x1-x2,y2-y1));
    2545 
    2546         length=sqrt(pow(x2-x1,2)+pow(y2-y1,2));
    2547 
    2548         Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
    2549         this->parameters->FindParam(&domaintype,DomainTypeEnum);
    2550         Input* vx_input=NULL;
    2551         Input* vy_input=NULL;
    2552         if(domaintype==Domain2DhorizontalEnum){
    2553                 vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
    2554                 vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
    2555         }
    2556         else{
    2557                 vx_input=inputs->GetInput(VxAverageEnum); _assert_(vx_input);
    2558                 vy_input=inputs->GetInput(VyAverageEnum); _assert_(vy_input);
    2559         }
    2560 
    2561         thickness_input->GetInputValue(&h1, gauss_1);
    2562         thickness_input->GetInputValue(&h2, gauss_2);
    2563         vx_input->GetInputValue(&vx1,gauss_1);
    2564         vx_input->GetInputValue(&vx2,gauss_2);
    2565         vy_input->GetInputValue(&vy1,gauss_1);
    2566         vy_input->GetInputValue(&vy2,gauss_2);
    2567 
    2568         mass_flux= rho_ice*length*( 
    2569                                 (ONETHIRD*(h1-h2)*(vx1-vx2)+0.5*h2*(vx1-vx2)+0.5*(h1-h2)*vx2+h2*vx2)*normal[0]+
    2570                                 (ONETHIRD*(h1-h2)*(vy1-vy2)+0.5*h2*(vy1-vy2)+0.5*(h1-h2)*vy2+h2*vy2)*normal[1]
    2571                                 );
    2572 
    2573         /*clean up and return:*/
    2574         delete gauss_1;
    2575         delete gauss_2;
    2576         return mass_flux;
     2950int        Tria::VertexConnectivity(int vertexindex){/*{{{*/
     2951        _assert_(this->vertices);
     2952        return this->vertices[vertexindex]->Connectivity();
    25772953}
    25782954/*}}}*/
    2579 IssmDouble Tria::MassFlux( IssmDouble* segment){/*{{{*/
     2955void       Tria::ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum){/*{{{*/
    25802956
    2581         int        domaintype;
    2582         IssmDouble mass_flux=0.;
    2583         IssmDouble xyz_list[NUMVERTICES][3];
    2584         IssmDouble normal[2];
    2585         IssmDouble length,rho_ice;
    2586         IssmDouble x1,y1,x2,y2,h1,h2;
    2587         IssmDouble vx1,vx2,vy1,vy2;
    2588         GaussTria* gauss_1=NULL;
    2589         GaussTria* gauss_2=NULL;
     2957        int         normal_orientation=0;
     2958        IssmDouble  s1,s2;
     2959        IssmDouble  levelset[NUMVERTICES];
    25902960
    2591         /*Get material parameters :*/
    2592         rho_ice=matpar->GetRhoIce();
     2961        /*Recover parameters and values*/
     2962        IssmDouble* xyz_zero = xNew<IssmDouble>(2*3);
     2963        GetInputListOnVertices(&levelset[0],levelsetenum);
    25932964
    2594         /*First off, check that this segment belongs to this element: */
    2595         if (reCast<int>(*(segment+4))!=this->id)_error_("error message: segment with id " << reCast<int>(*(segment+4)) << " does not belong to element with id:" << this->id);
     2965        if(levelset[0]*levelset[1]>0.){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
     2966                /*Portion of the segments*/
     2967                s1=levelset[2]/(levelset[2]-levelset[1]);
     2968                s2=levelset[2]/(levelset[2]-levelset[0]);
    25962969
    2597         /*Recover segment node locations: */
    2598         x1=*(segment+0); y1=*(segment+1); x2=*(segment+2); y2=*(segment+3);
     2970                if(levelset[2]<0.) normal_orientation=1; //orientation of quadrangle depending on distribution of levelsetfunction
     2971                /*New point 1*/
     2972                xyz_zero[3*normal_orientation+0]=xyz_list[2*3+0]+s1*(xyz_list[1*3+0]-xyz_list[2*3+0]);
     2973                xyz_zero[3*normal_orientation+1]=xyz_list[2*3+1]+s1*(xyz_list[1*3+1]-xyz_list[2*3+1]);
     2974                xyz_zero[3*normal_orientation+2]=xyz_list[2*3+2]+s1*(xyz_list[1*3+2]-xyz_list[2*3+2]);
    25992975
    2600         /*Get xyz list: */
    2601         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    2602 
    2603         /*get area coordinates of 0 and 1 locations: */
    2604         gauss_1=new GaussTria();
    2605         gauss_1->GaussFromCoords(x1,y1,&xyz_list[0][0]);
    2606         gauss_2=new GaussTria();
    2607         gauss_2->GaussFromCoords(x2,y2,&xyz_list[0][0]);
    2608 
    2609         normal[0]=cos(atan2(x1-x2,y2-y1));
    2610         normal[1]=sin(atan2(x1-x2,y2-y1));
    2611 
    2612         length=sqrt(pow(x2-x1,2)+pow(y2-y1,2));
    2613 
    2614         Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
    2615         this->parameters->FindParam(&domaintype,DomainTypeEnum);
    2616         Input* vx_input=NULL;
    2617         Input* vy_input=NULL;
    2618         if(domaintype==Domain2DhorizontalEnum){
    2619                 vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
    2620                 vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
     2976                /*New point 0*/
     2977                xyz_zero[3*(1-normal_orientation)+0]=xyz_list[2*3+0]+s2*(xyz_list[0*3+0]-xyz_list[2*3+0]);
     2978                xyz_zero[3*(1-normal_orientation)+1]=xyz_list[2*3+1]+s2*(xyz_list[0*3+1]-xyz_list[2*3+1]);
     2979                xyz_zero[3*(1-normal_orientation)+2]=xyz_list[2*3+2]+s2*(xyz_list[0*3+2]-xyz_list[2*3+2]);
    26212980        }
    2622         else{
    2623                 vx_input=inputs->GetInput(VxAverageEnum); _assert_(vx_input);
    2624                 vy_input=inputs->GetInput(VyAverageEnum); _assert_(vy_input);
    2625         }
     2981        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
     2982                /*Portion of the segments*/
     2983                s1=levelset[0]/(levelset[0]-levelset[2]);
     2984                s2=levelset[0]/(levelset[0]-levelset[1]);
    26262985
    2627         thickness_input->GetInputValue(&h1, gauss_1);
    2628         thickness_input->GetInputValue(&h2, gauss_2);
    2629         vx_input->GetInputValue(&vx1,gauss_1);
    2630         vx_input->GetInputValue(&vx2,gauss_2);
    2631         vy_input->GetInputValue(&vy1,gauss_1);
    2632         vy_input->GetInputValue(&vy2,gauss_2);
     2986                if(levelset[0]<0.) normal_orientation=1;
     2987                /*New point 1*/
     2988                xyz_zero[3*normal_orientation+0]=xyz_list[0*3+0]+s1*(xyz_list[2*3+0]-xyz_list[0*3+0]);
     2989                xyz_zero[3*normal_orientation+1]=xyz_list[0*3+1]+s1*(xyz_list[2*3+1]-xyz_list[0*3+1]);
     2990                xyz_zero[3*normal_orientation+2]=xyz_list[0*3+2]+s1*(xyz_list[2*3+2]-xyz_list[0*3+2]);
    26332991
    2634         mass_flux= rho_ice*length*( 
    2635                                 (ONETHIRD*(h1-h2)*(vx1-vx2)+0.5*h2*(vx1-vx2)+0.5*(h1-h2)*vx2+h2*vx2)*normal[0]+
    2636                                 (ONETHIRD*(h1-h2)*(vy1-vy2)+0.5*h2*(vy1-vy2)+0.5*(h1-h2)*vy2+h2*vy2)*normal[1]
    2637                                 );
    2638 
    2639         /*clean up and return:*/
    2640         delete gauss_1;
    2641         delete gauss_2;
    2642         return mass_flux;
    2643 }
    2644 /*}}}*/
    2645 void       Tria::ElementResponse(IssmDouble* presponse,int response_enum){/*{{{*/
    2646 
    2647         switch(response_enum){
    2648                 case MaterialsRheologyBbarEnum:
    2649                         *presponse=this->material->GetBbar();
    2650                         break;
    2651 
    2652                 case VelEnum:{
    2653 
    2654                         /*Get input:*/
    2655                         IssmDouble vel;
    2656                         Input* vel_input;
    2657 
    2658                         vel_input=this->inputs->GetInput(VelEnum); _assert_(vel_input);
    2659                         vel_input->GetInputAverage(&vel);
    2660 
    2661                         /*Assign output pointers:*/
    2662                         *presponse=vel;}
    2663                         break;
    2664                 default: 
    2665                         _error_("Response type " << EnumToStringx(response_enum) << " not supported yet!");
     2992                /*New point 2*/
     2993                xyz_zero[3*(1-normal_orientation)+0]=xyz_list[0*3+0]+s2*(xyz_list[1*3+0]-xyz_list[0*3+0]);
     2994                xyz_zero[3*(1-normal_orientation)+1]=xyz_list[0*3+1]+s2*(xyz_list[1*3+1]-xyz_list[0*3+1]);
     2995                xyz_zero[3*(1-normal_orientation)+2]=xyz_list[0*3+2]+s2*(xyz_list[1*3+2]-xyz_list[0*3+2]);
    26662996        }
     2997        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
     2998                /*Portion of the segments*/
     2999                s1=levelset[1]/(levelset[1]-levelset[0]);
     3000                s2=levelset[1]/(levelset[1]-levelset[2]);
    26673001
    2668 }
    2669 /*}}}*/
    2670 IssmDouble Tria::TotalSmb(void){/*{{{*/
     3002                if(levelset[1]<0.) normal_orientation=1;
     3003                /*New point 0*/
     3004                xyz_zero[3*normal_orientation+0]=xyz_list[1*3+0]+s1*(xyz_list[0*3+0]-xyz_list[1*3+0]);
     3005                xyz_zero[3*normal_orientation+1]=xyz_list[1*3+1]+s1*(xyz_list[0*3+1]-xyz_list[1*3+1]);
     3006                xyz_zero[3*normal_orientation+2]=xyz_list[1*3+2]+s1*(xyz_list[0*3+2]-xyz_list[1*3+2]);
    26713007
    2672         /*The smb[kg yr-1] of one element is area[m2] * smb [kg m^-2 yr^-1]*/
    2673         IssmDouble base,smb,rho_ice;
    2674         IssmDouble Total_Smb=0;
    2675         IssmDouble xyz_list[NUMVERTICES][3];
    2676 
    2677         /*Get material parameters :*/
    2678         rho_ice=matpar->GetRhoIce();
    2679 
    2680    if(!IsIceInElement())return 0;
    2681 
    2682         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    2683 
    2684         /*First calculate the area of the base (cross section triangle)
    2685          * http://en.wikipedia.org/wiki/Triangle
    2686          * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
    2687         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])); // area of element in m2
    2688 
    2689         /*Now get the average SMB over the element*/
    2690         Input* smb_input = inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(smb_input);
    2691         smb_input->GetInputAverage(&smb);                                                                                                                                                                                               // average smb on element in m ice s-1
    2692    Total_Smb=rho_ice*base*smb;                                                                                                                                                                                                                  // smb on element in kg s-1
    2693 
    2694         /*Return: */
    2695         return Total_Smb;
    2696 }
    2697 /*}}}*/
    2698 IssmDouble Tria::MisfitArea(int weightsenum){/*{{{*/
    2699 
    2700         /*Intermediaries*/
    2701         IssmDouble weight;
    2702         IssmDouble Jdet;
    2703         IssmDouble Jelem = 0;
    2704         IssmDouble xyz_list[NUMVERTICES][3];
    2705         GaussTria *gauss = NULL;
    2706 
    2707         /*If on water, return 0: */
    2708         if(!IsIceInElement())return 0;
    2709 
    2710         /*Retrieve all inputs we will be needing: */
    2711         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    2712         Input* weights_input     =inputs->GetInput(weightsenum);     _assert_(weights_input);
    2713 
    2714         /* Start  looping on the number of gaussian points: */
    2715         gauss=new GaussTria(2);
    2716         for(int ig=gauss->begin();ig<gauss->end();ig++){
    2717 
    2718                 gauss->GaussPoint(ig);
    2719 
    2720                 /* Get Jacobian determinant: */
    2721                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    2722 
    2723                 /*Get parameters at gauss point*/
    2724                 weights_input->GetInputValue(&weight,gauss);
    2725 
    2726                 /*compute misfit between model and observation */
    2727                 Jelem+=Jdet*weight*gauss->weight;
     3008                /*New point 2*/
     3009                xyz_zero[3*(1-normal_orientation)+0]=xyz_list[1*3+0]+s2*(xyz_list[2*3+0]-xyz_list[1*3+0]);
     3010                xyz_zero[3*(1-normal_orientation)+1]=xyz_list[1*3+1]+s2*(xyz_list[2*3+1]-xyz_list[1*3+1]);
     3011                xyz_zero[3*(1-normal_orientation)+2]=xyz_list[1*3+2]+s2*(xyz_list[2*3+2]-xyz_list[1*3+2]);
    27283012        }
     3013        else if(levelset[0]==0. && levelset[1]==0.){ //front is on point 0 and 1
     3014                xyz_zero[3*0+0]=xyz_list[0*3+0];
     3015                xyz_zero[3*0+1]=xyz_list[0*3+1];
     3016                xyz_zero[3*0+2]=xyz_list[0*3+2];
    27293017
    2730         /* clean up and Return: */
    2731         delete gauss;
    2732         return Jelem;
    2733 }
    2734 /*}}}*/
    2735 IssmDouble Tria::Misfit(int modelenum,int observationenum,int weightsenum){/*{{{*/
    2736 
    2737         /*Intermediaries*/
    2738         IssmDouble model,observation,weight;
    2739         IssmDouble Jdet;
    2740         IssmDouble Jelem = 0;
    2741         IssmDouble xyz_list[NUMVERTICES][3];
    2742         GaussTria *gauss = NULL;
    2743 
    2744         /*If on water, return 0: */
    2745         if(!IsIceInElement())return 0;
    2746 
    2747         /*Retrieve all inputs we will be needing: */
    2748         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    2749         Input* model_input=inputs->GetInput(modelenum);   _assert_(model_input);
    2750         Input* observation_input=inputs->GetInput(observationenum);_assert_(observation_input);
    2751         Input* weights_input     =inputs->GetInput(weightsenum);     _assert_(weights_input);
    2752 
    2753         /* Start  looping on the number of gaussian points: */
    2754         gauss=new GaussTria(2);
    2755         for(int ig=gauss->begin();ig<gauss->end();ig++){
    2756 
    2757                 gauss->GaussPoint(ig);
    2758 
    2759                 /* Get Jacobian determinant: */
    2760                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    2761 
    2762                 /*Get parameters at gauss point*/
    2763                 model_input->GetInputValue(&model,gauss);
    2764                 observation_input->GetInputValue(&observation,gauss);
    2765                 weights_input->GetInputValue(&weight,gauss);
    2766 
    2767                 /*compute misfit between model and observation */
    2768                 Jelem+=0.5*(model-observation)*(model-observation)*Jdet*weight*gauss->weight;
     3018                /*New point 2*/
     3019                xyz_zero[3*1+0]=xyz_list[1*3+0];
     3020                xyz_zero[3*1+1]=xyz_list[1*3+1];
     3021                xyz_zero[3*1+2]=xyz_list[1*3+2];
    27693022        }
     3023        else if(levelset[0]==0. && levelset[2]==0.){ //front is on point 0 and 1
     3024                xyz_zero[3*0+0]=xyz_list[2*3+0];
     3025                xyz_zero[3*0+1]=xyz_list[2*3+1];
     3026                xyz_zero[3*0+2]=xyz_list[2*3+2];
    27703027
    2771         /* clean up and Return: */
    2772         delete gauss;
    2773         return Jelem;
    2774 }
    2775 /*}}}*/
    2776 IssmDouble Tria::Masscon(IssmDouble* levelset){ /*{{{*/
    2777 
    2778 
    2779         /*intermediary: */
    2780         IssmDouble* values=NULL;
    2781         Input*      thickness_input=NULL;
    2782         IssmDouble  thickness;
    2783         IssmDouble  weight;
    2784         IssmDouble  Jdet;
    2785         IssmDouble  volume;
    2786         IssmDouble  rho_ice;
    2787         IssmDouble* xyz_list=NULL;
    2788         int         point1;
    2789         IssmDouble  fraction1,fraction2;
    2790         bool        mainlynegative=true;
    2791        
    2792         /*Output:*/
    2793         volume=0;
    2794 
    2795         /* Get node coordinates and dof list: */
    2796         GetVerticesCoordinates(&xyz_list);
    2797 
    2798         /*Retrieve inputs required:*/
    2799         thickness_input=this->GetInput(ThicknessEnum); _assert_(thickness_input);
    2800        
    2801         /*Retrieve material parameters: */
    2802         rho_ice=matpar->GetRhoIce();
    2803 
    2804         /*Retrieve values of the levelset defining the masscon: */
    2805         values = xNew<IssmDouble>(NUMVERTICES);
    2806         for(int i=0;i<NUMVERTICES;i++){
    2807                 values[i]=levelset[this->vertices[i]->Sid()];
     3028                /*New point 2*/
     3029                xyz_zero[3*1+0]=xyz_list[0*3+0];
     3030                xyz_zero[3*1+1]=xyz_list[0*3+1];
     3031                xyz_zero[3*1+2]=xyz_list[0*3+2];
    28083032        }
    2809                
    2810         /*Ok, use the level set values to figure out where we put our gaussian points:*/
    2811         this->GetLevelsetPositivePart(&point1,&fraction1,&fraction2,&mainlynegative,values);
    2812         Gauss* gauss = this->NewGauss(point1,fraction1,fraction2,mainlynegative,4);
     3033        else if(levelset[1]==0. && levelset[2]==0.){ //front is on point 0 and 1
     3034                xyz_zero[3*0+0]=xyz_list[1*3+0];
     3035                xyz_zero[3*0+1]=xyz_list[1*3+1];
     3036                xyz_zero[3*0+2]=xyz_list[1*3+2];
    28133037
    2814         volume=0;
    2815 
    2816         for(int ig=gauss->begin();ig<gauss->end();ig++){
    2817                 gauss->GaussPoint(ig);
    2818 
    2819                 this->JacobianDeterminant(&Jdet,xyz_list,gauss);
    2820                 thickness_input->GetInputValue(&thickness, gauss);
    2821 
    2822                 volume+=thickness*gauss->weight*Jdet;
     3038                /*New point 2*/
     3039                xyz_zero[3*1+0]=xyz_list[2*3+0];
     3040                xyz_zero[3*1+1]=xyz_list[2*3+1];
     3041                xyz_zero[3*1+2]=xyz_list[2*3+2];
    28233042        }
     3043        else _error_("Case not covered");
    28243044
    2825         /* clean up and Return: */
    2826         xDelete<IssmDouble>(xyz_list);
    2827         xDelete<IssmDouble>(values);
    2828         delete gauss;
    2829         return rho_ice*volume;
     3045        /*Assign output pointer*/
     3046        *pxyz_zero= xyz_zero;
    28303047}
    28313048/*}}}*/
    28323049
    … …  
    29583175/*}}}*/
    29593176#endif
    29603177
    2961 void       Tria::InputControlUpdate(IssmDouble scalar,bool save_parameter){/*{{{*/
     3178#ifdef _HAVE_DAKOTA_
     3179void       Tria::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type){/*{{{*/
    29623180
    2963         /*Intermediary*/
    2964         int    num_controls;
    2965         int*   control_type=NULL;
    2966         Input* input=NULL;
     3181        int             i,t,row;
     3182        IssmDouble      time;
     3183        TransientInput *transientinput = NULL;
     3184        IssmDouble      values[3];
    29673185
    2968         /*retrieve some parameters: */
    2969         this->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
    2970         this->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
     3186        /*Check that name is an element input*/
     3187        if (!IsInput(name)) return;
    29713188
    2972         for(int i=0;i<num_controls;i++){
    2973                 input=(Input*)this->inputs->GetInput(control_type[i]);   _assert_(input);
    2974                 if (input->ObjectEnum()!=ControlInputEnum){
    2975                         _error_("input " << EnumToStringx(control_type[i]) << " is not a ControlInput");
    2976                 }
     3189        switch(type){
    29773190
    2978                 ((ControlInput*)input)->UpdateValue(scalar);
    2979                 ((ControlInput*)input)->Constrain();
    2980                 if (save_parameter) ((ControlInput*)input)->SaveValue();
     3191                case VertexEnum:
     3192                        /*Create transient input: */
     3193                        for(t=0;t<ncols;t++){ //ncols is the number of times
    29813194
    2982         }
     3195                                /*create input values: */
     3196                                for(i=0;i<3;i++){
     3197                                        row=this->vertices[i]->Sid();
     3198                                        values[i]=matrix[ncols*row+t];
     3199                                }
    29833200
    2984         /*Clean up and return*/
    2985         xDelete<int>(control_type);
    2986 }
    2987 /*}}}*/
    2988 void       Tria::ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){/*{{{*/
     3201                                /*time:*/
     3202                                time=matrix[(nrows-1)*ncols+t];
    29893203
    2990         int    vertexpidlist[NUMVERTICES];
    2991         IssmDouble grad_list[NUMVERTICES];
    2992         Input* grad_input=NULL;
     3204                                if(t==0) transientinput=new TransientInput(name);
     3205                                transientinput->AddTimeInput(new TriaInput(name,values,P1Enum),time);
     3206                                transientinput->Configure(parameters);
     3207                        }
     3208                        this->inputs->AddInput(transientinput);
     3209                        break;
    29933210
    2994         Input* input=inputs->GetInput(enum_type);
    2995         if (!input) _error_("Input " << EnumToStringx(enum_type) << " not found");
    2996         if (input->ObjectEnum()!=ControlInputEnum) _error_("Input " << EnumToStringx(enum_type) << " is not a ControlInput");
    2997 
    2998         GradientIndexing(&vertexpidlist[0],control_index);
    2999         for(int i=0;i<NUMVERTICES;i++) grad_list[i]=gradient[vertexpidlist[i]];
    3000         grad_input=new TriaInput(GradientEnum,grad_list,P1Enum);
    3001 
    3002         ((ControlInput*)input)->SetGradient(grad_input);
    3003 
    3004 }/*}}}*/
    3005 void       Tria::ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum){/*{{{*/
    3006 
    3007         Input* input=inputs->GetInput(control_enum);
    3008         if (!input) _error_("Input " << EnumToStringx(control_enum) << " not found");
    3009         if (input->ObjectEnum()!=ControlInputEnum) _error_("Input " << EnumToStringx(control_enum) << " is not a ControlInput");
    3010 
    3011         int         sidlist[NUMVERTICES];
    3012         int         connectivity[NUMVERTICES];
    3013         IssmPDouble values[NUMVERTICES];
    3014         IssmPDouble gradients[NUMVERTICES];
    3015         IssmDouble  value,gradient;
    3016 
    3017         this->GetVerticesConnectivityList(&connectivity[0]);
    3018         this->GetVerticesSidList(&sidlist[0]);
    3019 
    3020         GaussTria* gauss=new GaussTria();
    3021         for (int iv=0;iv<NUMVERTICES;iv++){
    3022                 gauss->GaussVertex(iv);
    3023 
    3024                 ((ControlInput*)input)->GetInputValue(&value,gauss);
    3025                 ((ControlInput*)input)->GetGradientValue(&gradient,gauss);
    3026 
    3027                 values[iv]    = reCast<IssmPDouble>(value)/reCast<IssmPDouble>(connectivity[iv]);
    3028                 gradients[iv] = reCast<IssmPDouble>(gradient)/reCast<IssmPDouble>(connectivity[iv]);
     3211                default:
     3212                        _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet");
    30293213        }
    3030         delete gauss;
    30313214
    3032         vector_control->SetValues(NUMVERTICES,&sidlist[0],&values[0],ADD_VAL);
    3033         vector_gradient->SetValues(NUMVERTICES,&sidlist[0],&gradients[0],ADD_VAL);
    3034 
    3035 }/*}}}*/
    3036 void       Tria::GetVectorFromControlInputs(Vector<IssmDouble>* vector,int control_enum,int control_index,const char* data,bool onsid){/*{{{*/
    3037 
    3038         int vertexidlist[NUMVERTICES];
    3039         Input *input=NULL;
    3040 
    3041         /*Get out if this is not an element input*/
    3042         if(!IsInput(control_enum)) return;
    3043 
    3044         /*Prepare index list*/
    3045         GradientIndexing(&vertexidlist[0],control_index,onsid);
    3046 
    3047         /*Get input (either in element or material)*/
    3048         input=(Input*)this->inputs->GetInput(control_enum);   _assert_(input);
    3049 
    3050         /*Check that it is a ControlInput*/
    3051         if (input->ObjectEnum()!=ControlInputEnum){
    3052                 _error_("input " << EnumToStringx(control_enum) << " is not a ControlInput");
    3053         }
    3054 
    3055         ((ControlInput*)input)->GetVectorFromInputs(vector,&vertexidlist[0],data);
    30563215}
    30573216/*}}}*/
    3058 void       Tria::SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index){/*{{{*/
    3059 
    3060         IssmDouble  values[NUMVERTICES];
    3061         int         vertexpidlist[NUMVERTICES],control_init;
    3062 
    3063 
    3064         /*Get Domain type*/
    3065         int domaintype;
    3066         parameters->FindParam(&domaintype,DomainTypeEnum);
    3067 
    3068         /*Specific case for depth averaged quantities*/
    3069         control_init=control_enum;
    3070         if(domaintype==Domain2DverticalEnum){
    3071                 if(control_enum==MaterialsRheologyBbarEnum){
    3072                         control_enum=MaterialsRheologyBEnum;
    3073                         if(!IsOnBase()) return;
    3074                 }
    3075                 if(control_enum==DamageDbarEnum){
    3076                         control_enum=DamageDEnum;
    3077                         if(!IsOnBase()) return;
    3078                 }
    3079         }
    3080 
    3081         /*Get out if this is not an element input*/
    3082         if(!IsInput(control_enum)) return;
    3083 
    3084         /*Prepare index list*/
    3085         GradientIndexing(&vertexpidlist[0],control_index);
    3086 
    3087         /*Get values on vertices*/
    3088         for(int i=0;i<NUMVERTICES;i++){
    3089                 values[i]=vector[vertexpidlist[i]];
    3090         }
    3091         Input* new_input = new TriaInput(control_enum,values,P1Enum);
    3092         Input* input     = (Input*)this->inputs->GetInput(control_enum);   _assert_(input);
    3093         if(input->ObjectEnum()!=ControlInputEnum){
    3094                 _error_("input " << EnumToStringx(control_enum) << " is not a ControlInput");
    3095         }
    3096 
    3097         ((ControlInput*)input)->SetInput(new_input);
    3098 }
    3099 /*}}}*/
    3100 void       Tria::GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution, int enum_type){/*{{{*/
    3101 
    3102         int        *doflist = NULL;
    3103         IssmDouble  value;
    3104 
    3105         /*Fetch number of nodes for this finite element*/
    3106         int numnodes = this->NumberofNodes(this->element_type);
    3107 
    3108         /*Fetch dof list and allocate solution vector*/
    3109         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    3110         IssmDouble* values = xNew<IssmDouble>(numnodes);
    3111 
    3112         /*Get inputs*/
    3113         Input* enum_input=inputs->GetInput(enum_type); _assert_(enum_input);
    3114 
    3115         /*Ok, we have the values, fill in the array: */
    3116         GaussTria* gauss=new GaussTria();
    3117         for(int i=0;i<numnodes;i++){
    3118                 gauss->GaussNode(this->element_type,i);
    3119 
    3120                 enum_input->GetInputValue(&value,gauss);
    3121                 values[i]=value;
    3122         }
    3123 
    3124         solution->SetValues(numnodes,doflist,values,INS_VAL);
    3125 
    3126         /*Free ressources:*/
    3127         xDelete<int>(doflist);
    3128         xDelete<IssmDouble>(values);
    3129         delete gauss;
    3130 }
    3131 /*}}}*/
    3132 
    3133 #ifdef _HAVE_DAKOTA_
    31343217void       Tria::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){/*{{{*/
    31353218
    31363219        int i,j;
    … …  
    32223305
    32233306}
    32243307/*}}}*/
    3225 void       Tria::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type){/*{{{*/
    3226 
    3227         int             i,t,row;
    3228         IssmDouble      time;
    3229         TransientInput *transientinput = NULL;
    3230         IssmDouble      values[3];
    3231 
    3232         /*Check that name is an element input*/
    3233         if (!IsInput(name)) return;
    3234 
    3235         switch(type){
    3236 
    3237                 case VertexEnum:
    3238                         /*Create transient input: */
    3239                         for(t=0;t<ncols;t++){ //ncols is the number of times
    3240 
    3241                                 /*create input values: */
    3242                                 for(i=0;i<3;i++){
    3243                                         row=this->vertices[i]->Sid();
    3244                                         values[i]=matrix[ncols*row+t];
    3245                                 }
    3246 
    3247                                 /*time:*/
    3248                                 time=matrix[(nrows-1)*ncols+t];
    3249 
    3250                                 if(t==0) transientinput=new TransientInput(name);
    3251                                 transientinput->AddTimeInput(new TriaInput(name,values,P1Enum),time);
    3252                                 transientinput->Configure(parameters);
    3253                         }
    3254                         this->inputs->AddInput(transientinput);
    3255                         break;
    3256 
    3257                 default:
    3258                         _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet");
    3259         }
    3260 
    3261 }
    3262 /*}}}*/
    32633308#endif
    3264 
    3265 void       Tria::PotentialUngrounding(Vector<IssmDouble>* potential_ungrounding){/*{{{*/
    3266 
    3267         IssmDouble  h[NUMVERTICES],r[NUMVERTICES],gl[NUMVERTICES];
    3268         IssmDouble  bed_hydro;
    3269         IssmDouble  rho_water,rho_ice,density;
    3270 
    3271         /*material parameters: */
    3272         rho_water=matpar->GetRhoWater();
    3273         rho_ice=matpar->GetRhoIce();
    3274         density=rho_ice/rho_water;
    3275         GetInputListOnVertices(&h[0],ThicknessEnum);
    3276         GetInputListOnVertices(&r[0],BedEnum);
    3277         GetInputListOnVertices(&gl[0],MaskGroundediceLevelsetEnum);
    3278 
    3279         /*go through vertices, and figure out which ones are grounded and want to unground: */
    3280         for(int i=0;i<NUMVERTICES;i++){
    3281                 /*Find if grounded vertices want to start floating*/
    3282                 if (gl[i]>0.){
    3283                         bed_hydro=-density*h[i];
    3284                         if(bed_hydro>r[i]){
    3285                                 /*Vertex that could potentially unground, flag it*/
    3286                                 potential_ungrounding->SetValue(vertices[i]->Pid(),1,INS_VAL);
    3287                         }
    3288                 }
    3289         }
    3290 }
    3291 /*}}}*/
    3292 int        Tria::UpdatePotentialUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf){/*{{{*/
    3293 
    3294         int i;
    3295         int nflipped=0;
    3296 
    3297         /*Go through nodes, and whoever is on the potential_ungrounding, ends up in nodes_on_iceshelf: */
    3298         for(i=0;i<3;i++){
    3299                 if (reCast<bool>(vertices_potentially_ungrounding[vertices[i]->Pid()])){
    3300                         vec_nodes_on_iceshelf->SetValue(vertices[i]->Pid(),-1.,INS_VAL);
    3301 
    3302                         /*If node was not on ice shelf, we flipped*/
    3303                         if(nodes_on_iceshelf[vertices[i]->Pid()]>=0.){
    3304                                 nflipped++;
    3305                         }
    3306                 }
    3307         }
    3308         return nflipped;
    3309 }
    3310 /*}}}*/

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

     
    4141                Object *copy();
    4242                /*}}}*/
    4343                /*Update virtual functions resolution: {{{*/
    44                 void  InputUpdateFromVector(IssmDouble* vector, int name, int type);
    4544                #ifdef _HAVE_DAKOTA_
     45                void  InputUpdateFromMatrixDakota(IssmDouble* matrix, int nows, int ncols, int name, int type);
    4646                void  InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type);
    47                 void  InputUpdateFromMatrixDakota(IssmDouble* matrix, int nows, int ncols, int name, int type);
    4847                #endif
    4948                void  InputUpdateFromIoModel(int index, IoModel* iomodel);
     49                void  InputUpdateFromVector(IssmDouble* vector, int name, int type);
    5050                /*}}}*/
    5151                /*Element virtual functions definitions: {{{*/
     52                void        AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part);
     53                void                    CalvingRateLevermann();
    5254                IssmDouble  CharacteristicLength(void);
    5355                void        ComputeBasalStress(Vector<IssmDouble>* sigma_b);
     56                void        ComputeDeviatoricStressTensor();
    5457                void        ComputeSigmaNN();
    5558                void        ComputeStressTensor();
    56                 void        ComputeDeviatoricStressTensor();
    57                 void                    StrainRateparallel();
    58                 void                    StrainRateperpendicular();
    5959                void        ComputeSurfaceNormalVelocity();
    60                 void        StressIntensityFactor(void){_error_("not implemented yet");};
    61                 void                    CalvingRateLevermann();
    6260                void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters);
    63                 void        SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters);
    64                 void        ResetHooks();
     61                void        ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index);
     62                void        ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum);
    6563                void        Delta18oParameterization(void);
     64                int         EdgeOnBaseIndex();
     65                void        EdgeOnBaseIndices(int* pindex1,int* pindex);
     66                int         EdgeOnSurfaceIndex();
     67                void        EdgeOnSurfaceIndices(int* pindex1,int* pindex);
     68                void        ElementResponse(IssmDouble* presponse,int response_enum);
    6669                void        ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz);
     70                int         FiniteElement(void);
    6771                void        FSContactMigration(Vector<IssmDouble>* vertexgrounded,Vector<IssmDouble>* vertexfloating);
    68                 int         FiniteElement(void);
    69                 Element*    GetUpperElement(void){_error_("not implemented yet");};
    7072                Element*    GetBasalElement(void){_error_("not implemented yet");};
    7173                void        GetLevelsetPositivePart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlynegative,IssmDouble* levelsetvalues);
    7274                void        GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating);
    7375                IssmDouble  GetGroundedPortion(IssmDouble* xyz_list);
     76                void          GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum);
     77                void          GetLevelCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum,IssmDouble level);
    7478                int         GetNodeIndex(Node* node);
    7579                int         GetNumberOfNodes(void);
    7680                int         GetNumberOfNodes(int enum_type);
    7781                int         GetNumberOfVertices(void);
    78                 bool        IsOnBase();
    79                 bool        IsOnSurface();
    80                 bool        HasEdgeOnBase();
    81                 bool        HasEdgeOnSurface();
    82                 void        EdgeOnSurfaceIndices(int* pindex1,int* pindex);
    83                 void        EdgeOnBaseIndices(int* pindex1,int* pindex);
    84                 int         EdgeOnBaseIndex();
    85                 int         EdgeOnSurfaceIndex();
    86                 bool        IsNodeOnShelfFromFlags(IssmDouble* flags);
    87                 int         NumberofNodesVelocity(void);
    88                 int         NumberofNodesPressure(void);
    8982                void        GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution,int enum_type);
     83                Element*    GetUpperElement(void){_error_("not implemented yet");};
     84                void        GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,const char* data,bool onsid);
    9085                void        GetVerticesCoordinatesBase(IssmDouble** pxyz_list);
    9186                void        GetVerticesCoordinatesTop(IssmDouble** pxyz_list);
     87                bool        HasEdgeOnBase();
     88                bool        HasEdgeOnSurface();
     89                IssmDouble  IceVolume(void);
     90                IssmDouble  IceVolumeAboveFloatation(void);
     91                void        InputControlUpdate(IssmDouble scalar,bool save_parameter);
    9292                void        InputDepthAverageAtBase(int enum_type,int average_enum_type);
    9393                void        InputExtrude(int enum_type,int start){_error_("not implemented"); /*For penta only*/};
    9494                void        InputScale(int enum_type,IssmDouble scale_factor);
     95                bool            IsFaceOnBoundary(void);
     96                bool            IsIcefront(void);
     97                bool        IsNodeOnShelfFromFlags(IssmDouble* flags);
     98                bool        IsOnBase();
     99                bool        IsOnSurface();
     100                bool        IsZeroLevelset(int levelset_enum);
     101                IssmDouble  Masscon(IssmDouble* levelset);
     102                IssmDouble  MassFlux(IssmDouble* segment);
     103                IssmDouble  MassFlux(IssmDouble x1,IssmDouble y1, IssmDouble x2, IssmDouble y2,int segment_id);
    95104                void        MaterialUpdateFromTemperature(void){_error_("not implemented yet");};
     105                IssmDouble  Misfit(int modelenum,int observationenum,int weightsenum);
     106                IssmDouble  MisfitArea(int weightsenum);
    96107                int         NodalValue(IssmDouble* pvalue, int index, int natureofdataenum);
     108                int         NumberofNodesPressure(void);
     109                int         NumberofNodesVelocity(void);
    97110                void        PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm);
     111                void        PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding);
     112                int         PressureInterpolation();
    98113                void        ReduceMatrices(ElementMatrix* Ke,ElementVector* pe);
    99114                void        ResetFSBasalBoundaryCondition(void);
     115                void        ResetHooks();
     116                void        SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index);
     117                void        SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters);
    100118           Element*    SpawnBasalElement(void);
    101119                Element*    SpawnTopElement(void);
    102                 int         VelocityInterpolation();
    103                 int         PressureInterpolation();
     120                void                    StrainRateparallel();
     121                void                    StrainRateperpendicular();
     122                void        StressIntensityFactor(void){_error_("not implemented yet");};
     123                IssmDouble  SurfaceArea(void);
    104124                int         TensorInterpolation();
    105                 IssmDouble  SurfaceArea(void);
     125                IssmDouble  TimeAdapt();
     126                IssmDouble  TotalSmb(void);
    106127                void        Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finitelement);
    107                 IssmDouble  TimeAdapt();
     128                int         UpdatePotentialUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf);
     129                void        ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss);
    108130                void        ValueP1OnGauss(IssmDouble* pvalue,IssmDouble* values,Gauss* gauss);
    109                 void        ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss);
     131                int         VelocityInterpolation();
    110132                int         VertexConnectivity(int vertexindex);
    111133                void        VerticalSegmentIndices(int** pindices,int* pnumseg){_error_("not implemented yet");};
    112134                void        ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum);
    113                 void          GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum);
    114                 void          GetLevelCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum,IssmDouble level);
    115                 bool        IsZeroLevelset(int levelset_enum);
    116                 bool            IsIcefront(void);
    117                 bool            IsFaceOnBoundary(void);
    118135
    119                 void       AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part);
    120                 IssmDouble IceVolume(void);
    121                 IssmDouble IceVolumeAboveFloatation(void);
    122                 IssmDouble TotalSmb(void);
    123                 IssmDouble MassFlux(IssmDouble* segment);
    124                 IssmDouble MassFlux(IssmDouble x1,IssmDouble y1, IssmDouble x2, IssmDouble y2,int segment_id);
    125                 void       ElementResponse(IssmDouble* presponse,int response_enum);
    126                 IssmDouble Masscon(IssmDouble* levelset);
    127                 IssmDouble Misfit(int modelenum,int observationenum,int weightsenum);
    128                 IssmDouble MisfitArea(int weightsenum);
    129 
    130136                #ifdef _HAVE_GIA_
    131137                void   GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,IssmDouble* x,IssmDouble* y);
    132138                #endif
    133 
    134                 void       GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,const char* data,bool onsid);
    135                 void       SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index);
    136                 void       ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index);
    137                 void       ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum);
    138                 void       InputControlUpdate(IssmDouble scalar,bool save_parameter);
    139 
    140                 void   PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding);
    141                 int    UpdatePotentialUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf);
    142 
    143139                /*}}}*/
    144140                /*Tria specific routines:{{{*/
    145141                void           AddBasalInput(int input_enum, IssmDouble* values, int interpolation_enum);
    … …  
    147143                IssmDouble     GetArea(void);
    148144                void           GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints);
    149145                int            GetElementType(void);
    150                 void           NormalSection(IssmDouble* normal,IssmDouble* xyz_list);
    151                 void           NormalTop(IssmDouble* normal,IssmDouble* xyz_list);
    152                 void           NormalBase(IssmDouble* normal,IssmDouble* xyz_list);
    153146                void           GetInputValue(IssmDouble* pvalue,Node* node,int enumtype);
    154147                void           GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype);
    155148                Node*          GetNode(int node_number);
    156149                void             InputUpdateFromSolutionOneDof(IssmDouble* solution,int enum_type);
    157150                void             InputUpdateFromSolutionOneDofCollapsed(IssmDouble* solution,int enum_type){_error_("not implemented yet");};
    158151                void           JacobianDeterminant(IssmDouble*  pJdet, IssmDouble* xyz_list,Gauss* gauss);
     152                void           JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss);
    159153                void           JacobianDeterminantLine(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
    160154                void           JacobianDeterminantSurface(IssmDouble*  pJdet, IssmDouble* xyz_list,Gauss* gauss);
    161                 void           JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss);
    162155                void           JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss);
    163156                IssmDouble     MinEdgeLength(IssmDouble* xyz_list){_error_("not implemented yet");};
    164157                Gauss*         NewGauss(void);
    … …  
    170163                Gauss*         NewGaussLine(int vertex1,int vertex2,int order){_error_("not implemented yet");};
    171164                Gauss*         NewGaussTop(int order);
    172165                void           NodalFunctions(IssmDouble* basis,Gauss* gauss);
    173                 void           NodalFunctionsP1(IssmDouble* basis,Gauss* gauss);
    174                 void           NodalFunctionsP2(IssmDouble* basis,Gauss* gauss);
    175166                void           NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
    176                 void           NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
     167                void           NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
    177168                void           NodalFunctionsMINIDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
    178                 void           NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
    179                 void           NodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss);
    180169                void           NodalFunctionsPressure(IssmDouble* basis,Gauss* gauss);
     170                void           NodalFunctionsP1(IssmDouble* basis,Gauss* gauss);
     171                void           NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
     172                void           NodalFunctionsP2(IssmDouble* basis,Gauss* gauss);
    181173                void           NodalFunctionsTensor(IssmDouble* basis,Gauss* gauss);
     174                void           NodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss);
     175                void           NormalBase(IssmDouble* normal,IssmDouble* xyz_list);
     176                void           NormalSection(IssmDouble* normal,IssmDouble* xyz_list);
     177                void           NormalTop(IssmDouble* normal,IssmDouble* xyz_list);
    182178                void             SetClone(int* minranks);
    183179                void           SetTemporaryElementType(int element_type_in){_error_("not implemented yet");};
    184180                Seg*             SpawnSeg(int index1,int index2);
    185181                IssmDouble     StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa){_error_("not implemented yet");};
     182                void           UpdateConstraintsExtrudeFromBase(void);
     183                void           UpdateConstraintsExtrudeFromTop(void);
    186184                void           ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not implemented yet");};
    187 
    188                 void UpdateConstraintsExtrudeFromBase(void);
    189                 void UpdateConstraintsExtrudeFromTop(void);
    190185                /*}}}*/
    191186
    192187};
Note: See TracBrowser for help on using the repository browser.