Ignore:
Timestamp:
09/12/11 12:42:33 (14 years ago)
Author:
Eric.Larour
Message:

Reduced compilation to the maximum when adic flag is on, so that we do not have to forward differentiate the whole code.
Added 3D, Rifts, Transient, Prognostic, Steadystate conditional compilation.
Added eraselinks in src/ad/ directory to get rid of symlinks.
From mods in ad/Makefile.am, applied to c/Makefile.am -> better compilation.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r9761 r9775  
    376376}
    377377/*}}}*/
    378 /*FUNCTION Tria::CreateHydrologyWaterVelocityInput {{{1*/
    379 void Tria::CreateHydrologyWaterVelocityInput(void){
    380 
    381         /*material parameters: */
    382         double mu_water=0.001787;  //unit= [N s/m2]
    383         double w;
    384         double rho_ice, rho_water, g;
    385         double dsdx,dsdy,dbdx,dbdy;
    386         double vx[NUMVERTICES];
    387         double vy[NUMVERTICES];
    388         GaussTria *gauss = NULL;
    389 
    390         /*Retrieve all inputs and parameters*/
    391         rho_ice=matpar->GetRhoIce();
    392         rho_water=matpar->GetRhoWater();
    393         g=matpar->GetG();
    394         Input* surfaceslopex_input=inputs->GetInput(SurfaceSlopeXEnum); _assert_(surfaceslopex_input);
    395         Input* surfaceslopey_input=inputs->GetInput(SurfaceSlopeYEnum); _assert_(surfaceslopey_input);
    396         Input* bedslopex_input=inputs->GetInput(BedSlopeXEnum);         _assert_(bedslopex_input);
    397         Input* bedslopey_input=inputs->GetInput(BedSlopeYEnum);         _assert_(bedslopey_input);
    398         Input* watercolumn_input=inputs->GetInput(WatercolumnEnum);     _assert_(watercolumn_input);
    399 
    400         gauss=new GaussTria();
    401         for (int iv=0;iv<NUMVERTICES;iv++){
    402                 gauss->GaussVertex(iv);
    403                 surfaceslopex_input->GetParameterValue(&dsdx,gauss);
    404                 surfaceslopey_input->GetParameterValue(&dsdy,gauss);
    405                 bedslopex_input->GetParameterValue(&dbdx,gauss);
    406                 bedslopey_input->GetParameterValue(&dbdy,gauss);
    407                 watercolumn_input->GetParameterValue(&w,gauss);
    408 
    409                 /* Water velocity x and y components */
    410                 vx[iv]= - pow(w,2)/(12 * mu_water)*(rho_ice*g*dsdx+(rho_water-rho_ice)*g*dbdx);
    411                 vy[iv]= - pow(w,2)/(12 * mu_water)*(rho_ice*g*dsdy+(rho_water-rho_ice)*g*dbdy);
    412                 //vx[iv]=0.001;
    413                 //vy[iv]=0;
    414         }
    415 
    416         /*clean-up*/
    417         delete gauss;
    418 
    419         /*Add to inputs*/
    420         this->inputs->AddInput(new TriaVertexInput(HydrologyWaterVxEnum,vx));
    421         this->inputs->AddInput(new TriaVertexInput(HydrologyWaterVyEnum,vy));
    422 }
    423 /*}}}*/
    424378/*FUNCTION Tria::CreateKMatrix {{{1*/
    425379void  Tria::CreateKMatrix(Mat Kff, Mat Kfs,Vec df){
    … …  
    439393        /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
    440394        switch(analysis_type){
     395                #ifdef _HAVE_DIAGNOSTIC_
    441396                case DiagnosticHorizAnalysisEnum: case AdjointHorizAnalysisEnum:
    442397                        Ke=CreateKMatrixDiagnosticMacAyeal();
    … …  
    445400                        Ke=CreateKMatrixDiagnosticHutter();
    446401                        break;
     402                 #endif
    447403                case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum:
    448404                        Ke=CreateKMatrixSlope();
    … …  
    451407                        Ke=CreateKMatrixPrognostic();
    452408                        break;
     409                #ifdef _HAVE_HYDROLOGY_
    453410                case HydrologyAnalysisEnum:
    454411                        Ke=CreateKMatrixHydrology();
    455412                        break;
     413                #endif
     414                #ifdef _HAVE_BALANCED_
    456415                case BalancethicknessAnalysisEnum:
    457416                        Ke=CreateKMatrixBalancethickness();
    458417                        break;
     418                #endif
     419                #ifdef _HAVE_CONTROL_
    459420                case AdjointBalancethicknessAnalysisEnum:
    460421                        Ke=CreateKMatrixAdjointBalancethickness();
    461422                        break;
     423                #endif
    462424                default:
    463425                        _error_("analysis %i (%s) not supported yet",analysis_type,EnumToStringx(analysis_type));
    … …  
    469431                delete Ke;
    470432        }
    471 }
    472 /*}}}*/
    473 /*FUNCTION Tria::CreateKMatrixAdjointBalancethickness {{{1*/
    474 ElementMatrix* Tria::CreateKMatrixAdjointBalancethickness(void){
    475 
    476         ElementMatrix* Ke=NULL;
    477 
    478         /*Get Element Matrix of the forward model*/
    479         switch(GetElementType()){
    480                 case P1Enum:
    481                         Ke=CreateKMatrixBalancethickness_CG();
    482                         break;
    483                 case P1DGEnum:
    484                         Ke=CreateKMatrixBalancethickness_DG();
    485                         break;
    486                 default:
    487                         _error_("Element type %s not supported yet",EnumToStringx(GetElementType()));
    488         }
    489 
    490         /*Transpose and return Ke*/
    491         Ke->Transpose();
    492         return Ke;
    493 }
    494 /*}}}*/
    495 /*FUNCTION Tria::CreateKMatrixBalancethickness {{{1*/
    496 ElementMatrix* Tria::CreateKMatrixBalancethickness(void){
    497 
    498         switch(GetElementType()){
    499                 case P1Enum:
    500                         return CreateKMatrixBalancethickness_CG();
    501                 case P1DGEnum:
    502                         return CreateKMatrixBalancethickness_DG();
    503                 default:
    504                         _error_("Element type %s not supported yet",EnumToStringx(GetElementType()));
    505         }
    506 
    507 }
    508 /*}}}*/
    509 /*FUNCTION Tria::CreateKMatrixBalancethickness_CG {{{1*/
    510 ElementMatrix* Tria::CreateKMatrixBalancethickness_CG(void){
    511 
    512         /*Constants*/
    513         const int    numdof=NDOF1*NUMVERTICES;
    514 
    515         /*Intermediaries */
    516         int        stabilization;
    517         int        i,j,ig,dim;
    518         double     Jdettria,vx,vy,dvxdx,dvydy,vel,h;
    519         double     dvx[2],dvy[2];
    520         double     xyz_list[NUMVERTICES][3];
    521         double     L[NUMVERTICES];
    522         double     B[2][NUMVERTICES];
    523         double     Bprime[2][NUMVERTICES];
    524         double     K[2][2]                          = {0.0};
    525         double     KDL[2][2]                        = {0.0};
    526         double     DL[2][2]                         = {0.0};
    527         double     DLprime[2][2]                    = {0.0};
    528         double     DL_scalar;
    529         GaussTria *gauss                            = NULL;
    530 
    531         /*Initialize Element matrix*/
    532         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
    533 
    534         /*Retrieve all Inputs and parameters: */
    535         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    536         this->parameters->FindParam(&stabilization,BalancethicknessStabilizationEnum);
    537         this->parameters->FindParam(&dim,MeshDimensionEnum);
    538         Input* vxaverage_input=NULL;
    539         Input* vyaverage_input=NULL;
    540         if(dim==2){
    541                 vxaverage_input=inputs->GetInput(VxEnum); _assert_(vxaverage_input);
    542                 vyaverage_input=inputs->GetInput(VyEnum); _assert_(vyaverage_input);
    543         }
    544         else{
    545                 vxaverage_input=inputs->GetInput(VxAverageEnum); _assert_(vxaverage_input);
    546                 vyaverage_input=inputs->GetInput(VyAverageEnum); _assert_(vyaverage_input);
    547         }
    548         h=sqrt(2*this->GetArea());
    549 
    550         ///*Create Artificial diffusivity once for all if requested*/
    551         //if(stabilization){
    552         //      gauss=new GaussTria();
    553         //      gauss->GaussCenter();
    554         //      GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
    555         //      delete gauss;
    556 
    557         //      vxaverage_input->GetParameterAverage(&vx);
    558         //      vyaverage_input->GetParameterAverage(&vy);
    559         //      K[0][0]=pow(Jdettria,(double).5)/2.0*fabs(vx);
    560         //      K[1][1]=pow(Jdettria,(double).5)/2.0*fabs(vy);
    561         //}
    562 
    563         /*Start looping on the number of gaussian points:*/
    564         gauss=new GaussTria(2);
    565         for (ig=gauss->begin();ig<gauss->end();ig++){
    566 
    567                 gauss->GaussPoint(ig);
    568 
    569                 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
    570                 GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss);
    571                 GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
    572 
    573                 vxaverage_input->GetParameterValue(&vx,gauss);
    574                 vyaverage_input->GetParameterValue(&vy,gauss);
    575                 vxaverage_input->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
    576                 vyaverage_input->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
    577 
    578                 dvxdx=dvx[0];
    579                 dvydy=dvy[1];
    580                 DL_scalar=gauss->weight*Jdettria;
    581 
    582                 DL[0][0]=DL_scalar*dvxdx;
    583                 DL[1][1]=DL_scalar*dvydy;
    584 
    585                 DLprime[0][0]=DL_scalar*vx;
    586                 DLprime[1][1]=DL_scalar*vy;
    587 
    588                 TripleMultiply( &B[0][0],2,numdof,1,
    589                                         &DL[0][0],2,2,0,
    590                                         &B[0][0],2,numdof,0,
    591                                         &Ke->values[0],1);
    592 
    593                 TripleMultiply( &B[0][0],2,numdof,1,
    594                                         &DLprime[0][0],2,2,0,
    595                                         &Bprime[0][0],2,numdof,0,
    596                                         &Ke->values[0],1);
    597 
    598                 if(stabilization==1){
    599                         vel=sqrt(pow(vx,2.)+pow(vy,2.));
    600                         K[0][0]=h/(2*vel)*vx*vx;
    601                         K[1][0]=h/(2*vel)*vy*vx;
    602                         K[0][1]=h/(2*vel)*vx*vy;
    603                         K[1][1]=h/(2*vel)*vy*vy;
    604                         KDL[0][0]=DL_scalar*K[0][0];
    605                         KDL[1][0]=DL_scalar*K[1][0];
    606                         KDL[0][1]=DL_scalar*K[0][1];
    607                         KDL[1][1]=DL_scalar*K[1][1];
    608 
    609                         //KDL[0][0]=DL_scalar*K[0][0];
    610                         //KDL[1][1]=DL_scalar*K[1][1];
    611 
    612                         TripleMultiply( &Bprime[0][0],2,numdof,1,
    613                                                 &KDL[0][0],2,2,0,
    614                                                 &Bprime[0][0],2,numdof,0,
    615                                                 &Ke->values[0],1);
    616                 }
    617         }
    618 
    619         /*Clean up and return*/
    620         delete gauss;
    621         return Ke;
    622 }
    623 /*}}}*/
    624 /*FUNCTION Tria::CreateKMatrixBalancethickness_DG {{{1*/
    625 ElementMatrix* Tria::CreateKMatrixBalancethickness_DG(void){
    626 
    627         /*Constants*/
    628         const int  numdof=NDOF1*NUMVERTICES;
    629 
    630         /*Intermediaries*/
    631         int        i,j,ig,dim;
    632         double     vx,vy,Jdettria;
    633         double     xyz_list[NUMVERTICES][3];
    634         double     B[2][NUMVERTICES];
    635         double     Bprime[2][NUMVERTICES];
    636         double     DL[2][2]={0.0};
    637         double     DL_scalar;
    638         GaussTria  *gauss=NULL;
    639 
    640         /*Initialize Element matrix*/
    641         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
    642 
    643         /*Retrieve all inputs and parameters*/
    644         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    645         this->parameters->FindParam(&dim,MeshDimensionEnum);
    646         Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
    647         Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
    648 
    649         /*Start looping on the number of gaussian points:*/
    650         gauss=new GaussTria(2);
    651         for (ig=gauss->begin();ig<gauss->end();ig++){
    652 
    653                 gauss->GaussPoint(ig);
    654 
    655                 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
    656                 /*WARNING: B and Bprime are inverted compared to usual prognostic!!!!*/
    657                 GetBPrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
    658                 GetBprimePrognostic(&B[0][0], &xyz_list[0][0], gauss);
    659 
    660                 vx_input->GetParameterValue(&vx,gauss);
    661                 vy_input->GetParameterValue(&vy,gauss);
    662 
    663                 DL_scalar=-gauss->weight*Jdettria;
    664                 DL[0][0]=DL_scalar*vx;
    665                 DL[1][1]=DL_scalar*vy;
    666 
    667                 TripleMultiply( &B[0][0],2,numdof,1,
    668                                         &DL[0][0],2,2,0,
    669                                         &Bprime[0][0],2,numdof,0,
    670                                         &Ke->values[0],1);
    671         }
    672 
    673         /*Clean up and return*/
    674         delete gauss;
    675         return Ke;
    676 }
    677 /*}}}*/
    678 /*FUNCTION Tria::CreateKMatrixDiagnosticMacAyeal {{{1*/
    679 ElementMatrix* Tria::CreateKMatrixDiagnosticMacAyeal(void){
    680 
    681         /*compute all stiffness matrices for this element*/
    682         ElementMatrix* Ke1=CreateKMatrixDiagnosticMacAyealViscous();
    683         ElementMatrix* Ke2=CreateKMatrixDiagnosticMacAyealFriction();
    684         ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
    685        
    686         /*clean-up and return*/
    687         delete Ke1;
    688         delete Ke2;
    689         return Ke;
    690 }
    691 /*}}}*/
    692 /*FUNCTION Tria::CreateKMatrixDiagnosticMacAyealViscous{{{1*/
    693 ElementMatrix* Tria::CreateKMatrixDiagnosticMacAyealViscous(void){
    694 
    695         /*Constants*/
    696         const int  numdof=NDOF2*NUMVERTICES;
    697 
    698         /*Intermediaries*/
    699         int        i,j,ig;
    700         double     xyz_list[NUMVERTICES][3];
    701         double     viscosity,newviscosity,oldviscosity;
    702         double     viscosity_overshoot,thickness,Jdet;
    703         double     epsilon[3],oldepsilon[3];    /* epsilon=[exx,eyy,exy];    */
    704         double     B[3][numdof];
    705         double     Bprime[3][numdof];
    706         double     D[3][3]   = {0.0};
    707         double     D_scalar;
    708         GaussTria *gauss = NULL;
    709 
    710         /*Initialize Element matrix*/
    711         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum);
    712 
    713         /*Retrieve all inputs and parameters*/
    714         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    715         Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
    716         Input* vx_input=inputs->GetInput(VxEnum);               _assert_(vx_input);
    717         Input* vy_input=inputs->GetInput(VyEnum);               _assert_(vy_input);
    718         Input* vxold_input=inputs->GetInput(VxPicardEnum);      _assert_(vxold_input);
    719         Input* vyold_input=inputs->GetInput(VyPicardEnum);      _assert_(vyold_input);
    720         this->parameters->FindParam(&viscosity_overshoot,DiagnosticViscosityOvershootEnum);
    721 
    722         /* Start  looping on the number of gaussian points: */
    723         gauss=new GaussTria(2);
    724         for (ig=gauss->begin();ig<gauss->end();ig++){
    725 
    726                 gauss->GaussPoint(ig);
    727 
    728                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    729                 GetBMacAyeal(&B[0][0], &xyz_list[0][0], gauss);
    730                 GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss);
    731 
    732                 this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    733                 this->GetStrainRate2d(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input);
    734                 matice->GetViscosity2d(&viscosity, &epsilon[0]);
    735                 matice->GetViscosity2d(&oldviscosity, &oldepsilon[0]);
    736                 thickness_input->GetParameterValue(&thickness, gauss);
    737 
    738                 newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
    739                 D_scalar=2*newviscosity*thickness*gauss->weight*Jdet;
    740                 for (i=0;i<3;i++) D[i][i]=D_scalar;
    741 
    742                 TripleMultiply(&B[0][0],3,numdof,1,
    743                                         &D[0][0],3,3,0,
    744                                         &Bprime[0][0],3,numdof,0,
    745                                         &Ke->values[0],1);
    746         }
    747 
    748         /*Clean up and return*/
    749         delete gauss;
    750         return Ke;
    751 }
    752 /*}}}*/
    753 /*FUNCTION Tria::CreateKMatrixDiagnosticMacAyealFriction {{{1*/
    754 ElementMatrix* Tria::CreateKMatrixDiagnosticMacAyealFriction(void){
    755 
    756         /*Constants*/
    757         const int  numdof=NDOF2*NUMVERTICES;
    758 
    759         /*Intermediaries*/
    760         int        i,j,ig;
    761         int        analysis_type;
    762         double     MAXSLOPE  = .06; // 6 %
    763         double     MOUNTAINKEXPONENT = 10;
    764         double     slope_magnitude,alpha2;
    765         double     Jdet;
    766         double     L[2][numdof];
    767         double     DL[2][2]  = {{ 0,0 },{0,0}};
    768         double     DL_scalar;
    769         double     slope[2]  = {0.0,0.0};
    770         double     xyz_list[NUMVERTICES][3];
    771         Friction  *friction = NULL;
    772         GaussTria *gauss    = NULL;
    773 
    774         /*Initialize Element matrix and return if necessary*/
    775         if(IsOnShelf()) return NULL;
    776         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum);
    777 
    778         /*Retrieve all inputs and parameters*/
    779         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    780         Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input);
    781         Input* vx_input=inputs->GetInput(VxEnum);           _assert_(vx_input);
    782         Input* vy_input=inputs->GetInput(VyEnum);           _assert_(vy_input);
    783         Input* vz_input=inputs->GetInput(VzEnum);           _assert_(vz_input);
    784         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    785 
    786         /*build friction object, used later on: */
    787         friction=new Friction("2d",inputs,matpar,analysis_type);
    788 
    789         /* Start  looping on the number of gaussian points: */
    790         gauss=new GaussTria(2);
    791         for (ig=gauss->begin();ig<gauss->end();ig++){
    792 
    793                 gauss->GaussPoint(ig);
    794 
    795                 // If we have a slope > 6% for this element,  it means  we are on a mountain. In this particular case,
    796                 //velocity should be = 0. To achieve this result, we set alpha2_list to a very high value: */
    797                 surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
    798                 slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2));
    799                 if(slope_magnitude>MAXSLOPE) alpha2=pow((double)10,MOUNTAINKEXPONENT);
    800                 else friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum);
    801 
    802                 GetL(&L[0][0], &xyz_list[0][0], gauss,NDOF2);
    803                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    804                 DL_scalar=alpha2*gauss->weight*Jdet;
    805                 for (i=0;i<2;i++) DL[i][i]=DL_scalar;
    806                
    807                 TripleMultiply( &L[0][0],2,numdof,1,
    808                                         &DL[0][0],2,2,0,
    809                                         &L[0][0],2,numdof,0,
    810                                         &Ke->values[0],1);
    811         }
    812 
    813         /*Clean up and return*/
    814         delete gauss;
    815         delete friction;
    816         return Ke;
    817 }
    818 /*}}}*/
    819 /*FUNCTION Tria::CreateKMatrixDiagnosticHutter{{{1*/
    820 ElementMatrix* Tria::CreateKMatrixDiagnosticHutter(void){
    821 
    822         /*Intermediaries*/
    823         const int numdof=NUMVERTICES*NDOF2;
    824         int    i,connectivity;
    825 
    826         /*Initialize Element matrix*/
    827         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
    828 
    829         /*Create Element matrix*/
    830         for(i=0;i<NUMVERTICES;i++){
    831                 connectivity=nodes[i]->GetConnectivity();
    832                 Ke->values[(2*i)*numdof  +(2*i)  ]=1/(double)connectivity;
    833                 Ke->values[(2*i+1)*numdof+(2*i+1)]=1/(double)connectivity;
    834         }
    835 
    836         /*Clean up and return*/
    837         return Ke;
    838 }
    839 /*}}}*/
    840 /*FUNCTION Tria::CreateKMatrixHydrology{{{1*/
    841 ElementMatrix* Tria::CreateKMatrixHydrology(void){
    842 
    843         /*Constants*/
    844         const int    numdof=NDOF1*NUMVERTICES;
    845 
    846         /*Intermediaries */
    847         int        artdiff;
    848         int        i,j,ig;
    849         double     Jdettria,DL_scalar,dt,h;
    850         double     vx,vy,vel,dvxdx,dvydy;
    851         double     dvx[2],dvy[2];
    852         double     v_gauss[2]={0.0};
    853         double     xyz_list[NUMVERTICES][3];
    854         double     L[NUMVERTICES];
    855         double     B[2][NUMVERTICES];
    856         double     Bprime[2][NUMVERTICES];
    857         double     K[2][2]                        ={0.0};
    858         double     KDL[2][2]                      ={0.0};
    859         double     DL[2][2]                        ={0.0};
    860         double     DLprime[2][2]                   ={0.0};
    861         GaussTria *gauss=NULL;
    862 
    863         /*Initialize Element matrix*/
    864         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
    865 
    866         /*Create water velocity vx and vy from current inputs*/
    867         CreateHydrologyWaterVelocityInput();
    868 
    869         /*Retrieve all inputs and parameters*/
    870         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    871         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    872         this->parameters->FindParam(&artdiff,HydrologyStabilizationEnum);
    873         Input* vx_input=inputs->GetInput(HydrologyWaterVxEnum); _assert_(vx_input);
    874         Input* vy_input=inputs->GetInput(HydrologyWaterVyEnum); _assert_(vy_input);
    875         h=this->GetArea();
    876 
    877         /* Start  looping on the number of gaussian points: */
    878         gauss=new GaussTria(2);
    879         for (ig=gauss->begin();ig<gauss->end();ig++){
    880 
    881                 gauss->GaussPoint(ig);
    882 
    883                 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
    884                 GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
    885 
    886                 vx_input->GetParameterValue(&vx,gauss);
    887                 vy_input->GetParameterValue(&vy,gauss);
    888                 vx_input->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
    889                 vy_input->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
    890 
    891                 DL_scalar=gauss->weight*Jdettria;
    892 
    893                 TripleMultiply( &L[0],1,numdof,1,
    894                                         &DL_scalar,1,1,0,
    895                                         &L[0],1,numdof,0,
    896                                         &Ke->values[0],1);
    897 
    898                 GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss);
    899                 GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
    900 
    901                 dvxdx=dvx[0];
    902                 dvydy=dvy[1];
    903                 DL_scalar=dt*gauss->weight*Jdettria;
    904 
    905                 DL[0][0]=DL_scalar*dvxdx;
    906                 DL[1][1]=DL_scalar*dvydy;
    907                 DLprime[0][0]=DL_scalar*vx;
    908                 DLprime[1][1]=DL_scalar*vy;
    909 
    910                 TripleMultiply( &B[0][0],2,numdof,1,
    911                                         &DL[0][0],2,2,0,
    912                                         &B[0][0],2,numdof,0,
    913                                         &Ke->values[0],1);
    914 
    915                 TripleMultiply( &B[0][0],2,numdof,1,
    916                                         &DLprime[0][0],2,2,0,
    917                                         &Bprime[0][0],2,numdof,0,
    918                                         &Ke->values[0],1);
    919 
    920                 if(artdiff){
    921                         vel=sqrt(pow(vx,2.)+pow(vy,2.));
    922                         K[0][0]=h/(2*vel)*vx*vx;
    923                         K[1][0]=h/(2*vel)*vy*vx;
    924                         K[0][1]=h/(2*vel)*vx*vy;
    925                         K[1][1]=h/(2*vel)*vy*vy;
    926                         KDL[0][0]=DL_scalar*K[0][0];
    927                         KDL[1][0]=DL_scalar*K[1][0];
    928                         KDL[0][1]=DL_scalar*K[0][1];
    929                         KDL[1][1]=DL_scalar*K[1][1];
    930 
    931                         TripleMultiply( &Bprime[0][0],2,numdof,1,
    932                                                 &KDL[0][0],2,2,0,
    933                                                 &Bprime[0][0],2,numdof,0,
    934                                                 &Ke->values[0],1);
    935                 }
    936         }
    937 
    938         /*Clean up and return*/
    939         delete gauss;
    940         return Ke;
    941433}
    942434/*}}}*/
    … …  
    1247739        /*Just branch to the correct load generator, according to the type of analysis we are carrying out: */
    1248740        switch(analysis_type){
     741                #ifdef _HAVE_DIAGNOSTIC_
    1249742                case DiagnosticHorizAnalysisEnum:
    1250743                        pe=CreatePVectorDiagnosticMacAyeal();
    … …  
    1253746                        pe=CreatePVectorDiagnosticHutter();
    1254747                        break;
     748                #endif
    1255749                case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum:
    1256750                        pe=CreatePVectorSlope();
    … …  
    1259753                        pe=CreatePVectorPrognostic();
    1260754                        break;
     755                #ifdef _HAVE_HYDROLOGY_
    1261756                case HydrologyAnalysisEnum:
    1262757                        pe=CreatePVectorHydrology();
    1263758                        break;
     759                #endif
     760                #ifdef _HAVE_BALANCED_
    1264761                case BalancethicknessAnalysisEnum:
    1265762                        pe=CreatePVectorBalancethickness();
    1266763                        break;
     764                #endif
    1267765                #ifdef _HAVE_CONTROL_
    1268766                case AdjointBalancethicknessAnalysisEnum:
    … …  
    1284782}
    1285783/*}}}*/
    1286 /*FUNCTION Tria::CreatePVectorBalancethickness{{{1*/
    1287 ElementVector* Tria::CreatePVectorBalancethickness(void){
    1288 
    1289         switch(GetElementType()){
    1290                 case P1Enum:
    1291                         return CreatePVectorBalancethickness_CG();
    1292                         break;
    1293                 case P1DGEnum:
    1294                         return CreatePVectorBalancethickness_DG();
    1295                 default:
    1296                         _error_("Element type %s not supported yet",EnumToStringx(GetElementType()));
    1297         }
    1298 }
    1299 /*}}}*/
    1300 /*FUNCTION Tria::CreatePVectorBalancethickness_CG{{{1*/
    1301 ElementVector* Tria::CreatePVectorBalancethickness_CG(void){
    1302 
    1303         /*Constants*/
    1304         const int    numdof=NDOF1*NUMVERTICES;
    1305        
    1306         /*Intermediaries */
    1307         int        i,j,ig;
    1308         double     xyz_list[NUMVERTICES][3];
    1309         double     dhdt_g,basal_melting_g,surface_mass_balance_g,Jdettria;
    1310         double     L[NUMVERTICES];
    1311         GaussTria* gauss=NULL;
    1312 
    1313         /*Initialize Element vector*/
    1314         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
    1315 
    1316         /*Retrieve all inputs and parameters*/
    1317         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    1318         Input* surface_mass_balance_input=inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(surface_mass_balance_input);
    1319         Input* basal_melting_input=inputs->GetInput(BasalforcingsMeltingRateEnum);          _assert_(basal_melting_input);
    1320         Input* dhdt_input=inputs->GetInput(BalancethicknessThickeningRateEnum);             _assert_(dhdt_input);
    1321        
    1322         /* Start  looping on the number of gaussian points: */
    1323         gauss=new GaussTria(2);
    1324         for(ig=gauss->begin();ig<gauss->end();ig++){
    1325 
    1326                 gauss->GaussPoint(ig);
    1327 
    1328                 surface_mass_balance_input->GetParameterValue(&surface_mass_balance_g,gauss);
    1329                 basal_melting_input->GetParameterValue(&basal_melting_g,gauss);
    1330                 dhdt_input->GetParameterValue(&dhdt_g,gauss);
    1331 
    1332                 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
    1333                 GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
    1334 
    1335                 for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(surface_mass_balance_g-basal_melting_g-dhdt_g)*L[i];
    1336         }
    1337 
    1338         /*Clean up and return*/
    1339         delete gauss;
    1340         return pe;
    1341 }
    1342 /*}}}*/
    1343 /*FUNCTION Tria::CreatePVectorBalancethickness_DG {{{1*/
    1344 ElementVector* Tria::CreatePVectorBalancethickness_DG(void){
    1345 
    1346         /*Constants*/
    1347         const int    numdof=NDOF1*NUMVERTICES;
    1348 
    1349         /*Intermediaries */
    1350         int        i,j,ig;
    1351         double     xyz_list[NUMVERTICES][3];
    1352         double     basal_melting_g,surface_mass_balance_g,dhdt_g,Jdettria;
    1353         double     L[NUMVERTICES];
    1354         GaussTria* gauss=NULL;
    1355 
    1356         /*Initialize Element vector*/
    1357         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
    1358 
    1359         /*Retrieve all inputs and parameters*/
    1360         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    1361         Input* surface_mass_balance_input=inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(surface_mass_balance_input);
    1362         Input* basal_melting_input=inputs->GetInput(BasalforcingsMeltingRateEnum);          _assert_(basal_melting_input);
    1363         Input* dhdt_input=inputs->GetInput(BalancethicknessThickeningRateEnum);                                       _assert_(dhdt_input);
    1364 
    1365         /* Start  looping on the number of gaussian points: */
    1366         gauss=new GaussTria(2);
    1367         for(ig=gauss->begin();ig<gauss->end();ig++){
    1368 
    1369                 gauss->GaussPoint(ig);
    1370 
    1371                 surface_mass_balance_input->GetParameterValue(&surface_mass_balance_g,gauss);
    1372                 basal_melting_input->GetParameterValue(&basal_melting_g,gauss);
    1373                 dhdt_input->GetParameterValue(&dhdt_g,gauss);
    1374 
    1375                 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
    1376                 GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
    1377 
    1378                 for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(surface_mass_balance_g-basal_melting_g-dhdt_g)*L[i];
    1379         }
    1380 
    1381         /*Clean up and return*/
    1382         delete gauss;
    1383         return pe;
    1384 }
    1385 /*}}}*/
    1386 /*FUNCTION Tria::CreatePVectorDiagnosticMacAyeal {{{1*/
    1387 ElementVector* Tria::CreatePVectorDiagnosticMacAyeal(){
    1388 
    1389         /*Constants*/
    1390         const int    numdof=NDOF2*NUMVERTICES;
    1391 
    1392         /*Intermediaries */
    1393         int            i,j,ig;
    1394         double         driving_stress_baseline,thickness;
    1395         double         Jdet;
    1396         double         xyz_list[NUMVERTICES][3];
    1397         double         slope[2];
    1398         double         basis[3];
    1399         double         pe_g_gaussian[numdof];
    1400         GaussTria*     gauss=NULL;
    1401 
    1402         /*Initialize Element vector*/
    1403         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum);
    1404 
    1405         /*Retrieve all inputs and parameters*/
    1406         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    1407         Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
    1408         Input* surface_input=inputs->GetInput(SurfaceEnum);     _assert_(surface_input);
    1409         Input* drag_input=inputs->GetInput(FrictionCoefficientEnum);_assert_(drag_input);
    1410 
    1411         /* Start  looping on the number of gaussian points: */
    1412         gauss=new GaussTria(2);
    1413         for(ig=gauss->begin();ig<gauss->end();ig++){
    1414 
    1415                 gauss->GaussPoint(ig);
    1416 
    1417                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    1418                 GetNodalFunctions(basis, gauss);
    1419 
    1420                 thickness_input->GetParameterValue(&thickness,gauss);
    1421                 surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
    1422                 driving_stress_baseline=matpar->GetRhoIce()*matpar->GetG()*thickness;
    1423 
    1424                 /*Build pe_g_gaussian vector: */
    1425                 for (i=0;i<NUMVERTICES;i++){
    1426                         for (j=0;j<NDOF2;j++){
    1427                                 pe->values[i*NDOF2+j]+=-driving_stress_baseline*slope[j]*Jdet*gauss->weight*basis[i];
    1428                         }
    1429                 }
    1430         }
    1431 
    1432         /*Clean up and return*/
    1433         delete gauss;
    1434         return pe;
    1435 }
    1436 /*}}}*/
    1437 /*FUNCTION Tria::CreatePVectorDiagnosticHutter{{{1*/
    1438 ElementVector* Tria::CreatePVectorDiagnosticHutter(void){
    1439 
    1440         /*Intermediaries */
    1441         int        i,connectivity;
    1442         double     constant_part,ub,vb;
    1443         double     rho_ice,gravity,n,B;
    1444         double     slope2,thickness;
    1445         double     slope[2];
    1446         GaussTria* gauss=NULL;
    1447 
    1448         /*Initialize Element vector*/
    1449         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
    1450 
    1451         /*Retrieve all inputs and parameters*/
    1452         rho_ice=matpar->GetRhoIce();
    1453         gravity=matpar->GetG();
    1454         n=matice->GetN();
    1455         B=matice->GetBbar();
    1456         Input* slopex_input=inputs->GetInput(SurfaceSlopeXEnum); _assert_(slopex_input);
    1457         Input* slopey_input=inputs->GetInput(SurfaceSlopeYEnum); _assert_(slopey_input);
    1458         Input* thickness_input=inputs->GetInput(ThicknessEnum);  _assert_(thickness_input);
    1459 
    1460         /*Spawn 3 sing elements: */
    1461         gauss=new GaussTria();
    1462         for(i=0;i<NUMVERTICES;i++){
    1463 
    1464                 gauss->GaussVertex(i);
    1465 
    1466                 connectivity=nodes[i]->GetConnectivity();
    1467 
    1468                 thickness_input->GetParameterValue(&thickness,gauss);
    1469                 slopex_input->GetParameterValue(&slope[0],gauss);
    1470                 slopey_input->GetParameterValue(&slope[1],gauss);
    1471                 slope2=pow(slope[0],2)+pow(slope[1],2);
    1472 
    1473                 constant_part=-2*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2));
    1474 
    1475                 ub=-1.58*pow((double)10.0,(double)-10.0)*rho_ice*gravity*thickness*slope[0];
    1476                 vb=-1.58*pow((double)10.0,(double)-10.0)*rho_ice*gravity*thickness*slope[1];
    1477 
    1478                 pe->values[2*i]  =(ub-2.0*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2.0))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[0])/(double)connectivity;
    1479                 pe->values[2*i+1]=(vb-2.0*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2.0))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[1])/(double)connectivity;
    1480         }
    1481 
    1482         /*Clean up and return*/
    1483         delete gauss;
    1484         return pe;
    1485 }
    1486 /*}}}*/
    1487784/*FUNCTION Tria::CreatePVectorPrognostic{{{1*/
    1488785ElementVector* Tria::CreatePVectorPrognostic(void){
    … …  
    1496793                        _error_("Element type %s not supported yet",EnumToStringx(GetElementType()));
    1497794        }
    1498 }
    1499 /*}}}*/
    1500 /*FUNCTION Tria::CreatePVectorHydrology {{{1*/
    1501 ElementVector* Tria::CreatePVectorHydrology(void){
    1502 
    1503         /*Constants*/
    1504         const int    numdof=NDOF1*NUMVERTICES;
    1505 
    1506         /*Intermediaries */
    1507         int        i,j,ig;
    1508         double     Jdettria,dt;
    1509         double     basal_melting_g;
    1510         double     old_watercolumn_g;
    1511         double     xyz_list[NUMVERTICES][3];
    1512         double     basis[numdof];
    1513         GaussTria* gauss=NULL;
    1514 
    1515         /*Initialize Element vector*/
    1516         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
    1517 
    1518         /*Retrieve all inputs and parameters*/
    1519         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    1520         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    1521         Input* basal_melting_input=inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(basal_melting_input);
    1522         Input* old_watercolumn_input=inputs->GetInput(WaterColumnOldEnum);         _assert_(old_watercolumn_input);
    1523 
    1524         /*Initialize basal_melting_correction_g to 0, do not forget!:*/
    1525         /* Start  looping on the number of gaussian points: */
    1526         gauss=new GaussTria(2);
    1527         for(ig=gauss->begin();ig<gauss->end();ig++){
    1528 
    1529                 gauss->GaussPoint(ig);
    1530 
    1531                 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
    1532                 GetNodalFunctions(basis, gauss);
    1533 
    1534                 basal_melting_input->GetParameterValue(&basal_melting_g,gauss);
    1535                 old_watercolumn_input->GetParameterValue(&old_watercolumn_g,gauss);
    1536 
    1537                 if(dt)for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(old_watercolumn_g+dt*basal_melting_g)*basis[i];
    1538                 else  for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*basal_melting_g*basis[i];
    1539         }
    1540                
    1541         /*Clean up and return*/
    1542         delete gauss;
    1543         return pe;
    1544795}
    1545796/*}}}*/
    … …  
    19931244
    19941245        /*Just branch to the correct InputUpdateFromSolution generator, according to the type of analysis we are carrying out: */
    1995         if (analysis_type==DiagnosticHorizAnalysisEnum)
    1996          GetSolutionFromInputsDiagnosticHoriz(solution);
    1997         else if (analysis_type==DiagnosticHutterAnalysisEnum)
    1998          GetSolutionFromInputsDiagnosticHutter(solution);
    1999         else if (analysis_type==HydrologyAnalysisEnum)
    2000          GetSolutionFromInputsHydrology(solution);
     1246        switch(analysis_type){
     1247        #ifdef _HAVE_DIAGNOSTIC_
     1248        case DiagnosticHorizAnalysisEnum:
     1249                GetSolutionFromInputsDiagnosticHoriz(solution);
     1250                break;
     1251        case DiagnosticHutterAnalysisEnum:
     1252                GetSolutionFromInputsDiagnosticHutter(solution);
     1253                break;
     1254        #endif
     1255        #ifdef _HAVE_HYDROLOGY_
     1256        case HydrologyAnalysisEnum:
     1257                GetSolutionFromInputsHydrology(solution);
     1258                break;
     1259        #endif
     1260        default:
     1261                _error_("analysis: %s not supported yet",EnumToStringx(analysis_type));
     1262        }
     1263
     1264}
     1265/*}}}*/
     1266/*FUNCTION Tria::GetStrainRate2d(double* epsilon,double* xyz_list, GaussTria* gauss, Input* vx_input, Input* vy_input){{{1*/
     1267void Tria::GetStrainRate2d(double* epsilon,double* xyz_list, GaussTria* gauss, Input* vx_input, Input* vy_input){
     1268        /*Compute the 2d Strain Rate (3 components):
     1269         * epsilon=[exx eyy exy] */
     1270
     1271        int i;
     1272        double epsilonvx[3];
     1273        double epsilonvy[3];
     1274
     1275        /*Check that both inputs have been found*/
     1276        if (!vx_input || !vy_input){
     1277                _error_("Input missing. Here are the input pointers we have for vx: %p, vy: %p\n",vx_input,vy_input);
     1278        }
     1279
     1280        /*Get strain rate assuming that epsilon has been allocated*/
     1281        vx_input->GetVxStrainRate2d(epsilonvx,xyz_list,gauss);
     1282        vy_input->GetVyStrainRate2d(epsilonvy,xyz_list,gauss);
     1283
     1284        /*Sum all contributions*/
     1285        for(i=0;i<3;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i];
     1286}
     1287/*}}}*/
     1288/*FUNCTION Tria::GetVectorFromInputs{{{1*/
     1289void  Tria::GetVectorFromInputs(Vec vector,int input_enum){
     1290
     1291        int doflist1[NUMVERTICES];
     1292
     1293        /*Get out if this is not an element input*/
     1294        if (!IsInput(input_enum)) return;
     1295
     1296        /*Prepare index list*/
     1297        this->GetDofList1(&doflist1[0]);
     1298
     1299        /*Get input (either in element or material)*/
     1300        Input* input=inputs->GetInput(input_enum);
     1301        if(!input) _error_("Input %s not found in element",EnumToStringx(input_enum));
     1302
     1303        /*We found the enum.  Use its values to fill into the vector, using the vertices ids: */
     1304        input->GetVectorFromInputs(vector,&doflist1[0]);
     1305}
     1306/*}}}*/
     1307/*FUNCTION Tria::Id {{{1*/
     1308int    Tria::Id(){
     1309       
     1310        return id;
     1311
     1312}
     1313/*}}}*/
     1314/*FUNCTION Tria::Sid {{{1*/
     1315int    Tria::Sid(){
     1316       
     1317        return sid;
     1318
     1319}
     1320/*}}}*/
     1321/*FUNCTION Tria::InputArtificialNoise{{{1*/
     1322void  Tria::InputArtificialNoise(int enum_type,double min,double max){
     1323
     1324        Input* input=NULL;
     1325
     1326        /*Make a copy of the original input: */
     1327        input=(Input*)this->inputs->GetInput(enum_type);
     1328        if(!input)_error_(" could not find old input with enum: %s",EnumToStringx(enum_type));
     1329
     1330        /*ArtificialNoise: */
     1331        input->ArtificialNoise(min,max);
     1332}
     1333/*}}}*/
     1334/*FUNCTION Tria::InputConvergence{{{1*/
     1335bool Tria::InputConvergence(double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){
     1336
     1337        bool    converged=true;
     1338        int     i;
     1339        Input** new_inputs=NULL;
     1340        Input** old_inputs=NULL;
     1341
     1342        new_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the new inputs
     1343        old_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the old inputs
     1344
     1345        for(i=0;i<num_enums/2;i++){
     1346                new_inputs[i]=(Input*)this->inputs->GetInput(enums[2*i+0]);
     1347                old_inputs[i]=(Input*)this->inputs->GetInput(enums[2*i+1]);
     1348                if(!new_inputs[i])_error_("%s%s"," could not find input with enum ",EnumToStringx(enums[2*i+0]));
     1349                if(!old_inputs[i])_error_("%s%s"," could not find input with enum ",EnumToStringx(enums[2*i+0]));
     1350        }
     1351
     1352        /*ok, we've got the inputs (new and old), now loop throught the number of criterions and fill the eps array:*/
     1353        for(i=0;i<num_criterionenums;i++){
     1354                IsInputConverged(eps+i,new_inputs,old_inputs,num_enums/2,criterionenums[i]);
     1355                if(eps[i]>criterionvalues[i]) converged=false;
     1356        }
     1357
     1358        /*clean up and return*/
     1359        xfree((void**)&new_inputs);
     1360        xfree((void**)&old_inputs);
     1361        return converged;
     1362}
     1363/*}}}*/
     1364/*FUNCTION Tria::InputDepthAverageAtBase {{{1*/
     1365void  Tria::InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum){
     1366
     1367        /*New input*/
     1368        Input* oldinput=NULL;
     1369        Input* newinput=NULL;
     1370
     1371        /*copy input of enum_type*/
     1372        if (object_enum==MeshElementsEnum)
     1373         oldinput=(Input*)this->inputs->GetInput(enum_type);
     1374        else if (object_enum==MaterialsEnum)
     1375         oldinput=(Input*)this->matice->inputs->GetInput(enum_type);
    20011376        else
    2002          _error_("analysis: %s not supported yet",EnumToStringx(analysis_type));
    2003 
     1377         _error_("object %s not supported yet",EnumToStringx(object_enum));
     1378        if(!oldinput)_error_("%s%s"," could not find old input with enum: ",EnumToStringx(enum_type));
     1379        newinput=(Input*)oldinput->copy();
     1380
     1381        /*Assign new name (average)*/
     1382        newinput->ChangeEnum(average_enum_type);
     1383
     1384        /*Add new input to current element*/
     1385        if (object_enum==MeshElementsEnum)
     1386         this->inputs->AddInput((Input*)newinput);
     1387        else if (object_enum==MaterialsEnum)
     1388         this->matice->inputs->AddInput((Input*)newinput);
     1389        else
     1390         _error_("object %s not supported yet",EnumToStringx(object_enum));
     1391}
     1392/*}}}*/
     1393/*FUNCTION Tria::InputDuplicate{{{1*/
     1394void  Tria::InputDuplicate(int original_enum,int new_enum){
     1395
     1396        /*Call inputs method*/
     1397        if (IsInput(original_enum)) inputs->DuplicateInput(original_enum,new_enum);
     1398
     1399}
     1400/*}}}*/
     1401/*FUNCTION Tria::InputScale{{{1*/
     1402void  Tria::InputScale(int enum_type,double scale_factor){
     1403
     1404        Input* input=NULL;
     1405
     1406        /*Make a copy of the original input: */
     1407        input=(Input*)this->inputs->GetInput(enum_type);
     1408        if(!input)_error_(" could not find old input with enum: %s",EnumToStringx(enum_type));
     1409
     1410        /*Scale: */
     1411        input->Scale(scale_factor);
     1412}
     1413/*}}}*/
     1414/*FUNCTION Tria::InputToResult{{{1*/
     1415void  Tria::InputToResult(int enum_type,int step,double time){
     1416
     1417        int    i;
     1418        Input *input = NULL;
     1419
     1420        /*Go through all the input objects, and find the one corresponding to enum_type, if it exists: */
     1421        if (enum_type==MaterialsRheologyBbarEnum) input=this->matice->inputs->GetInput(enum_type);
     1422        else input=this->inputs->GetInput(enum_type);
     1423        if (!input) _error_("Input %s not found in tria->inputs",EnumToStringx(enum_type));
     1424
     1425        /*If we don't find it, no big deal, just don't do the transfer. Otherwise, build a new Result
     1426         * object out of the input, with the additional step and time information: */
     1427        this->results->AddObject((Object*)input->SpawnResult(step,time));
     1428        #ifdef _HAVE_CONTROL_
     1429        if(input->Enum()==ControlInputEnum) this->results->AddObject((Object*)((ControlInput*)input)->SpawnGradient(step,time));
     1430        #endif
     1431}
     1432/*}}}*/
     1433/*FUNCTION Tria::InputUpdateFromConstant(int value, int name);{{{1*/
     1434void  Tria::InputUpdateFromConstant(int constant, int name){
     1435        /*Check that name is an element input*/
     1436        if (!IsInput(name)) return;
     1437
     1438        /*update input*/
     1439        this->inputs->AddInput(new IntInput(name,constant));
     1440}
     1441/*}}}*/
     1442/*FUNCTION Tria::InputUpdateFromConstant(double value, int name);{{{1*/
     1443void  Tria::InputUpdateFromConstant(double constant, int name){
     1444        /*Check that name is an element input*/
     1445        if (!IsInput(name)) return;
     1446
     1447        /*update input*/
     1448        this->inputs->AddInput(new DoubleInput(name,constant));
     1449}
     1450/*}}}*/
     1451/*FUNCTION Tria::InputUpdateFromConstant(bool value, int name);{{{1*/
     1452void  Tria::InputUpdateFromConstant(bool constant, int name){
     1453        /*Check that name is an element input*/
     1454        if (!IsInput(name)) return;
     1455
     1456        /*update input*/
     1457        this->inputs->AddInput(new BoolInput(name,constant));
     1458}
     1459/*}}}*/
     1460/*FUNCTION Tria::InputUpdateFromIoModel{{{1*/
     1461void Tria::InputUpdateFromIoModel(int index, IoModel* iomodel){ //i is the element index
     1462
     1463        /*Intermediaries*/
     1464        int    i,j;
     1465        int    tria_vertex_ids[3];
     1466        double nodeinputs[3];
     1467        double cmmininputs[3];
     1468        double cmmaxinputs[3];
     1469        bool   control_analysis=false;
     1470        int    num_control_type;
     1471        double yts;
     1472        int    num_cm_responses;
     1473   
     1474        /*Get parameters: */
     1475        iomodel->Constant(&control_analysis,InversionIscontrolEnum);
     1476        iomodel->Constant(&num_control_type,InversionNumControlParametersEnum);
     1477        iomodel->Constant(&yts,ConstantsYtsEnum);
     1478        iomodel->Constant(&num_cm_responses,InversionNumCostFunctionsEnum);
     1479
     1480        /*Recover vertices ids needed to initialize inputs*/
     1481        for(i=0;i<3;i++){
     1482                tria_vertex_ids[i]=(int)iomodel->Data(MeshElementsEnum)[3*index+i]; //ids for vertices are in the elements array from Matlab
     1483        }
     1484
     1485        /*Control Inputs*/
     1486        #ifdef _HAVE_CONTROL_
     1487        if (control_analysis && iomodel->Data(InversionControlParametersEnum)){
     1488                for(i=0;i<num_control_type;i++){
     1489                        switch((int)iomodel->Data(InversionControlParametersEnum)[i]){
     1490                                case BalancethicknessThickeningRateEnum:
     1491                                        if (iomodel->Data(BalancethicknessThickeningRateEnum)){
     1492                                                for(j=0;j<3;j++)nodeinputs[j]=iomodel->Data(BalancethicknessThickeningRateEnum)[tria_vertex_ids[j]-1]/yts;
     1493                                                for(j=0;j<3;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
     1494                                                for(j=0;j<3;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
     1495                                                this->inputs->AddInput(new ControlInput(BalancethicknessThickeningRateEnum,TriaVertexInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     1496                                        }
     1497                                        break;
     1498                                case VxEnum:
     1499                                        if (iomodel->Data(VxEnum)){
     1500                                                for(j=0;j<3;j++)nodeinputs[j]=iomodel->Data(VxEnum)[tria_vertex_ids[j]-1]/yts;
     1501                                                for(j=0;j<3;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
     1502                                                for(j=0;j<3;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
     1503                                                this->inputs->AddInput(new ControlInput(VxEnum,TriaVertexInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     1504                                        }
     1505                                        break;
     1506                                case VyEnum:
     1507                                        if (iomodel->Data(VyEnum)){
     1508                                                for(j=0;j<3;j++)nodeinputs[j]=iomodel->Data(VyEnum)[tria_vertex_ids[j]-1]/yts;
     1509                                                for(j=0;j<3;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
     1510                                                for(j=0;j<3;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
     1511                                                this->inputs->AddInput(new ControlInput(VyEnum,TriaVertexInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     1512                                        }
     1513                                        break;
     1514                                case FrictionCoefficientEnum:
     1515                                        if (iomodel->Data(FrictionCoefficientEnum)){
     1516                                                for(j=0;j<3;j++)nodeinputs[j]=iomodel->Data(FrictionCoefficientEnum)[tria_vertex_ids[j]-1];
     1517                                                for(j=0;j<3;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i];
     1518                                                for(j=0;j<3;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i];
     1519                                                this->inputs->AddInput(new ControlInput(FrictionCoefficientEnum,TriaVertexInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     1520                                        }
     1521                                        break;
     1522                                case MaterialsRheologyBbarEnum:
     1523                                        /*Matice will take care of it*/ break;
     1524                                default:
     1525                                        _error_("Control %s not implemented yet",EnumToStringx((int)iomodel->Data(InversionControlParametersEnum)[i]));
     1526                        }
     1527                }
     1528        }
     1529        #endif
     1530
     1531        /*DatasetInputs*/
     1532        if (control_analysis && iomodel->Data(InversionCostFunctionsCoefficientsEnum)){
     1533
     1534                /*Create inputs and add to DataSetInput*/
     1535                DatasetInput* datasetinput=new DatasetInput(InversionCostFunctionsCoefficientsEnum);
     1536                for(i=0;i<num_cm_responses;i++){
     1537                        for(j=0;j<3;j++)nodeinputs[j]=iomodel->Data(InversionCostFunctionsCoefficientsEnum)[(tria_vertex_ids[j]-1)*num_cm_responses+i];
     1538                        datasetinput->inputs->AddObject(new TriaVertexInput(InversionCostFunctionsCoefficientsEnum,nodeinputs));
     1539                }
     1540
     1541                /*Add datasetinput to element inputs*/
     1542                this->inputs->AddInput(datasetinput);
     1543        }
     1544}
     1545/*}}}*/
     1546/*FUNCTION Tria::InputUpdateFromSolution {{{1*/
     1547void  Tria::InputUpdateFromSolution(double* solution){
     1548
     1549        /*retrive parameters: */
     1550        int analysis_type;
     1551        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
     1552
     1553        /*Just branch to the correct InputUpdateFromSolution generator, according to the type of analysis we are carrying out: */
     1554        switch(analysis_type){
     1555                #ifdef _HAVE_DIAGNOSTIC_
     1556                case DiagnosticHorizAnalysisEnum:
     1557                        InputUpdateFromSolutionDiagnosticHoriz( solution);
     1558                        break;
     1559                case DiagnosticHutterAnalysisEnum:
     1560                        InputUpdateFromSolutionDiagnosticHoriz( solution);
     1561                        break;
     1562                #endif
     1563                #ifdef _HAVE_CONTROL_
     1564                case AdjointHorizAnalysisEnum:
     1565                        InputUpdateFromSolutionAdjointHoriz( solution);
     1566                        break;
     1567                case AdjointBalancethicknessAnalysisEnum:
     1568                        InputUpdateFromSolutionAdjointBalancethickness( solution);
     1569                        break;
     1570                #endif
     1571                #ifdef _HAVE_HYDROLOGY_
     1572                case HydrologyAnalysisEnum:
     1573                        InputUpdateFromSolutionHydrology(solution);
     1574                        break ;
     1575                #endif
     1576                #ifdef _HAVE_BALANCED_
     1577                case BalancethicknessAnalysisEnum:
     1578                        InputUpdateFromSolutionOneDof(solution,ThicknessEnum);
     1579                        break;
     1580                #endif
     1581                case BedSlopeXAnalysisEnum:
     1582                        InputUpdateFromSolutionOneDof(solution,BedSlopeXEnum);
     1583                        break;
     1584                case BedSlopeYAnalysisEnum:
     1585                        InputUpdateFromSolutionOneDof(solution,BedSlopeYEnum);
     1586                        break;
     1587                case SurfaceSlopeXAnalysisEnum:
     1588                        InputUpdateFromSolutionOneDof(solution,SurfaceSlopeXEnum);
     1589                        break;
     1590                case SurfaceSlopeYAnalysisEnum:
     1591                        InputUpdateFromSolutionOneDof(solution,SurfaceSlopeYEnum);
     1592                        break;
     1593                case PrognosticAnalysisEnum:
     1594                        InputUpdateFromSolutionPrognostic(solution);
     1595                        break;
     1596                default:
     1597                        _error_("analysis %i (%s) not supported yet",analysis_type,EnumToStringx(analysis_type));
     1598        }
     1599}
     1600/*}}}*/
     1601/*FUNCTION Tria::InputUpdateFromSolutionOneDof{{{1*/
     1602void  Tria::InputUpdateFromSolutionOneDof(double* solution,int enum_type){
     1603
     1604        const int numdof          = NDOF1*NUMVERTICES;
     1605
     1606        int*      doflist=NULL;
     1607        double    values[numdof];
     1608
     1609        /*Get dof list: */
     1610        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     1611
     1612        /*Use the dof list to index into the solution vector: */
     1613        for(int i=0;i<numdof;i++){
     1614                values[i]=solution[doflist[i]];
     1615                if(isnan(values[i])) _error_("NaN found in solution vector");
     1616        }
     1617
     1618        /*Add input to the element: */
     1619        this->inputs->AddInput(new TriaVertexInput(enum_type,values));
     1620
     1621        /*Free ressources:*/
     1622        xfree((void**)&doflist);
     1623}
     1624/*}}}*/
     1625/*FUNCTION Tria::InputUpdateFromSolutionPrognostic{{{1*/
     1626void  Tria::InputUpdateFromSolutionPrognostic(double* solution){
     1627
     1628        /*Intermediaries*/
     1629        const int numdof = NDOF1*NUMVERTICES;
     1630
     1631        int       i,dummy,hydroadjustment;
     1632        int*      doflist=NULL;
     1633        double    rho_ice,rho_water;
     1634        double    values[numdof];
     1635        double    bed[numdof];
     1636        double    surface[numdof];
     1637        double    *bed_ptr = NULL;
     1638        double    *thickness_ptr = NULL;
     1639        double    *surface_ptr = NULL;
     1640
     1641        /*Get dof list: */
     1642        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     1643
     1644        /*Use the dof list to index into the solution vector: */
     1645        for(i=0;i<numdof;i++){
     1646                values[i]=solution[doflist[i]];
     1647                if(isnan(values[i])) _error_("NaN found in solution vector");
     1648                /*Constrain thickness to be at least 1m*/
     1649                if(values[i]<1) values[i]=1;
     1650        }
     1651
     1652        /*Get previous bed, thickness and surface*/
     1653        Input* bed_input=inputs->GetInput(BedEnum);             _assert_(bed_input);
     1654        Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
     1655        Input* surface_input=inputs->GetInput(SurfaceEnum);     _assert_(surface_input);
     1656        bed_input->GetValuesPtr(&bed_ptr,&dummy);
     1657        thickness_input->GetValuesPtr(&thickness_ptr,&dummy);
     1658        surface_input->GetValuesPtr(&surface_ptr,&dummy);
     1659
     1660        /*Fing PrognosticHydrostaticAdjustment to figure out how to update the geometry:*/
     1661        this->parameters->FindParam(&hydroadjustment,PrognosticHydrostaticAdjustmentEnum);
     1662
     1663        /*recover material parameters: */
     1664        rho_ice=matpar->GetRhoIce();
     1665        rho_water=matpar->GetRhoWater();
     1666
     1667        for(i=0;i<numdof;i++) {
     1668                /*If shelf: hydrostatic equilibrium*/
     1669                if (this->nodes[i]->IsOnSheet()){
     1670                        surface[i]=bed_ptr[i]+values[i]; //surface=oldbed+newthickness
     1671                        bed[i]=bed_ptr[i]; //do nothing
     1672                }
     1673                else{ //this is an ice shelf
     1674
     1675                        if(hydroadjustment==AbsoluteEnum){
     1676                                surface[i]=values[i]*(1-rho_ice/rho_water);
     1677                                bed[i]=values[i]*(-rho_ice/rho_water);
     1678                        }
     1679                        else if(hydroadjustment==IncrementalEnum){
     1680                                surface[i]=surface_ptr[i]+(1.0-rho_ice/rho_water)*(values[i]-thickness_ptr[i]); //surface = oldsurface + (1-di) * dH
     1681                                bed[i]=bed_ptr[i]-rho_ice/rho_water*(values[i]-thickness_ptr[i]); //bed = oldbed + di * dH
     1682                        }
     1683                        else _error_("Hydrostatic adjustment %i (%s) not supported yet",hydroadjustment,EnumToStringx(hydroadjustment));
     1684                }
     1685        }
     1686
     1687        /*Add input to the element: */
     1688        this->inputs->AddInput(new TriaVertexInput(ThicknessEnum,values));
     1689        this->inputs->AddInput(new TriaVertexInput(SurfaceEnum,surface));
     1690        this->inputs->AddInput(new TriaVertexInput(BedEnum,bed));
     1691
     1692        /*Free ressources:*/
     1693        xfree((void**)&doflist);
     1694}
     1695/*}}}*/
     1696/*FUNCTION Tria::InputUpdateFromVector(double* vector, int name, int type);{{{1*/
     1697void  Tria::InputUpdateFromVector(double* vector, int name, int type){
     1698
     1699        /*Check that name is an element input*/
     1700        if (!IsInput(name)) return;
     1701
     1702        switch(type){
     1703
     1704                case VertexEnum:
     1705
     1706                        /*New TriaVertexInput*/
     1707                        double values[3];
     1708
     1709                        /*Get values on the 3 vertices*/
     1710                        for (int i=0;i<3;i++){
     1711                                values[i]=vector[this->nodes[i]->GetVertexDof()];
     1712                        }
     1713
     1714                        /*update input*/
     1715                        if (name==MaterialsRheologyBbarEnum || name==MaterialsRheologyBEnum){
     1716                                matice->inputs->AddInput(new TriaVertexInput(name,values));
     1717                        }
     1718                        else{
     1719                                this->inputs->AddInput(new TriaVertexInput(name,values));
     1720                        }
     1721                        return;
     1722
     1723                default:
     1724                        _error_("type %i (%s) not implemented yet",type,EnumToStringx(type));
     1725        }
     1726}
     1727/*}}}*/
     1728/*FUNCTION Tria::InputUpdateFromVector(int* vector, int name, int type);{{{1*/
     1729void  Tria::InputUpdateFromVector(int* vector, int name, int type){
     1730        _error_(" not supported yet!");
     1731}
     1732/*}}}*/
     1733/*FUNCTION Tria::InputUpdateFromVector(bool* vector, int name, int type);{{{1*/
     1734void  Tria::InputUpdateFromVector(bool* vector, int name, int type){
     1735        _error_(" not supported yet!");
     1736}
     1737/*}}}*/
     1738/*FUNCTION Tria::InputCreate(double scalar,int enum,int code);{{{1*/
     1739void Tria::InputCreate(double scalar,int name,int code){
     1740
     1741        /*Check that name is an element input*/
     1742        if (!IsInput(name)) return;
     1743       
     1744        if ((code==5) || (code==1)){ //boolean
     1745                this->inputs->AddInput(new BoolInput(name,(bool)scalar));
     1746        }
     1747        else if ((code==6) || (code==2)){ //integer
     1748                this->inputs->AddInput(new IntInput(name,(int)scalar));
     1749        }
     1750        else if ((code==7) || (code==3)){ //double
     1751                this->inputs->AddInput(new DoubleInput(name,(int)scalar));
     1752        }
     1753        else _error_("%s%i"," could not recognize nature of vector from code ",code);
     1754
     1755}
     1756/*}}}*/
     1757/*FUNCTION Tria::InputCreate(double* vector,int index,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){{{1*/
     1758void Tria::InputCreate(double* vector, int index,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){//index into elements
     1759
     1760        /*Intermediaries*/
     1761        int    i,j,t;
     1762        int    tria_vertex_ids[3];
     1763        int    row;
     1764        double nodeinputs[3];
     1765        double time;
     1766        TransientInput* transientinput=NULL;
     1767        int    numberofvertices;
     1768        int    numberofelements;
     1769        double yts;
     1770
     1771
     1772        /*Fetch parameters: */
     1773        iomodel->Constant(&numberofvertices,MeshNumberofverticesEnum);
     1774        iomodel->Constant(&numberofelements,MeshNumberofelementsEnum);
     1775        iomodel->Constant(&yts,ConstantsYtsEnum);
     1776
     1777        /*Branch on type of vector: nodal or elementary: */
     1778        if(vector_type==1){ //nodal vector
     1779
     1780                /*Recover vertices ids needed to initialize inputs*/
     1781                for(i=0;i<3;i++){
     1782                        tria_vertex_ids[i]=(int)iomodel->Data(MeshElementsEnum)[3*index+i]; //ids for vertices are in the elements array from Matlab
     1783                }
     1784
     1785                /*Are we in transient or static? */
     1786                if(M==numberofvertices){
     1787
     1788                        /*create input values: */
     1789                        for(i=0;i<3;i++)nodeinputs[i]=(double)vector[tria_vertex_ids[i]-1];
     1790
     1791                        /*process units: */
     1792                        UnitConversion(&nodeinputs[0], 3 ,ExtToIuEnum, vector_enum);
     1793
     1794                        /*create static input: */
     1795                        this->inputs->AddInput(new TriaVertexInput(vector_enum,nodeinputs));
     1796                }
     1797                else if(M==numberofvertices+1){
     1798                        /*create transient input: */
     1799                        for(t=0;t<N;t++){ //N is the number of times
     1800
     1801                                /*create input values: */
     1802                                for(i=0;i<3;i++){
     1803                                        row=tria_vertex_ids[i]-1;
     1804                                        nodeinputs[i]=(double)vector[N*row+t];
     1805                                }
     1806
     1807                                /*process units: */
     1808                                UnitConversion(&nodeinputs[0], 3 ,ExtToIuEnum, vector_enum);
     1809
     1810                                /*time? :*/
     1811                                time=(double)vector[(M-1)*N+t]*yts;
     1812
     1813                                if(t==0) transientinput=new TransientInput(vector_enum);
     1814                                transientinput->AddTimeInput(new TriaVertexInput(vector_enum,nodeinputs),time);
     1815                        }
     1816                        this->inputs->AddInput(transientinput);
     1817                }
     1818                else _error_("nodal vector is either numberofnodes (%i), or numberofnodes+1 long. Field provided is %i long",numberofvertices,M);
     1819        }
     1820        else if(vector_type==2){ //element vector
     1821                /*Are we in transient or static? */
     1822                if(M==numberofelements){
     1823
     1824                        /*static mode: create an input out of the element value: */
     1825
     1826                        if (code==5){ //boolean
     1827                                this->inputs->AddInput(new BoolInput(vector_enum,(bool)vector[index]));
     1828                        }
     1829                        else if (code==6){ //integer
     1830                                this->inputs->AddInput(new IntInput(vector_enum,(int)vector[index]));
     1831                        }
     1832                        else if (code==7){ //double
     1833                                this->inputs->AddInput(new DoubleInput(vector_enum,(double)vector[index]));
     1834                        }
     1835                        else _error_("%s%i"," could not recognize nature of vector from code ",code);
     1836                }
     1837                else {
     1838                        _error_("transient elementary inputs not supported yet!");
     1839                }
     1840        }
     1841        else{
     1842                _error_("Cannot add input for vector type %i (not supported)",vector_type);
     1843        }
     1844
     1845}
     1846/*}}}*/
     1847/*FUNCTION Tria::IsInput{{{1*/
     1848bool Tria::IsInput(int name){
     1849        if (
     1850                                name==ThicknessEnum ||
     1851                                name==SurfaceEnum ||
     1852                                name==BedEnum ||
     1853                                name==SurfaceSlopeXEnum ||
     1854                                name==SurfaceSlopeYEnum ||
     1855                                name==BasalforcingsMeltingRateEnum ||
     1856                                name==WatercolumnEnum ||
     1857                                name==SurfaceforcingsMassBalanceEnum ||
     1858                                name==SurfaceAreaEnum||
     1859                                name==VxEnum ||
     1860                                name==VyEnum ||
     1861                                name==InversionVxObsEnum ||
     1862                                name==InversionVyObsEnum ||
     1863                                name==FrictionCoefficientEnum ||
     1864                                name==GradientEnum ||
     1865                                name==OldGradientEnum
     1866                ){
     1867                return true;
     1868        }
     1869        else return false;
     1870}
     1871/*}}}*/
     1872/*FUNCTION Tria::IsOnBed {{{1*/
     1873bool Tria::IsOnBed(){
     1874       
     1875        bool onbed;
     1876        inputs->GetParameterValue(&onbed,MeshElementonbedEnum);
     1877        return onbed;
     1878}
     1879/*}}}*/
     1880/*FUNCTION Tria::IsOnShelf {{{1*/
     1881bool   Tria::IsOnShelf(){
     1882
     1883        bool shelf;
     1884        inputs->GetParameterValue(&shelf,MaskElementonfloatingiceEnum);
     1885        return shelf;
     1886}
     1887/*}}}*/
     1888/*FUNCTION Tria::IsNodeOnShelf {{{1*/
     1889bool   Tria::IsNodeOnShelf(){
     1890
     1891        int  i;
     1892        bool shelf=false;
     1893
     1894        for(i=0;i<3;i++){
     1895                if (nodes[i]->IsOnShelf()){
     1896                        shelf=true;
     1897                        break;
     1898                }
     1899        }
     1900        return shelf;
     1901}
     1902/*}}}*/
     1903/*FUNCTION Tria::IsNodeOnShelfFromFlags {{{1*/
     1904bool   Tria::IsNodeOnShelfFromFlags(double* flags){
     1905
     1906        int  i;
     1907        bool shelf=false;
     1908
     1909        for(i=0;i<3;i++){
     1910                if (flags[nodes[i]->Sid()]){
     1911                        shelf=true;
     1912                        break;
     1913                }
     1914        }
     1915        return shelf;
     1916}
     1917/*}}}*/
     1918/*FUNCTION Tria::IsOnWater {{{1*/
     1919bool   Tria::IsOnWater(){
     1920
     1921        bool water;
     1922        inputs->GetParameterValue(&water,MaskElementonwaterEnum);
     1923        return water;
     1924}
     1925/*}}}*/
     1926/*FUNCTION Tria::MassFlux {{{1*/
     1927double Tria::MassFlux( double* segment,bool process_units){
     1928
     1929        const int    numdofs=2;
     1930
     1931        int        i;
     1932        double     mass_flux=0;
     1933        double     xyz_list[NUMVERTICES][3];
     1934        double     normal[2];
     1935        double     length,rho_ice;
     1936        double     x1,y1,x2,y2,h1,h2;
     1937        double     vx1,vx2,vy1,vy2;
     1938        GaussTria* gauss_1=NULL;
     1939        GaussTria* gauss_2=NULL;
     1940
     1941        /*Get material parameters :*/
     1942        rho_ice=matpar->GetRhoIce();
     1943
     1944        /*First off, check that this segment belongs to this element: */
     1945        if ((int)*(segment+4)!=this->id)_error_("%s%i%s%i","error message: segment with id ",(int)*(segment+4)," does not belong to element with id:",this->id);
     1946
     1947        /*Recover segment node locations: */
     1948        x1=*(segment+0); y1=*(segment+1); x2=*(segment+2); y2=*(segment+3);
     1949
     1950        /*Get xyz list: */
     1951        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     1952
     1953        /*get area coordinates of 0 and 1 locations: */
     1954        gauss_1=new GaussTria();
     1955        gauss_1->GaussFromCoords(x1,y1,&xyz_list[0][0]);
     1956        gauss_2=new GaussTria();
     1957        gauss_2->GaussFromCoords(x2,y2,&xyz_list[0][0]);
     1958
     1959        normal[0]=cos(atan2(x1-x2,y2-y1));
     1960        normal[1]=sin(atan2(x1-x2,y2-y1));
     1961
     1962        length=sqrt(pow(x2-x1,2.0)+pow(y2-y1,2));
     1963
     1964        Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
     1965        Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
     1966        Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
     1967
     1968        thickness_input->GetParameterValue(&h1, gauss_1);
     1969        thickness_input->GetParameterValue(&h2, gauss_2);
     1970        vx_input->GetParameterValue(&vx1,gauss_1);
     1971        vx_input->GetParameterValue(&vx2,gauss_2);
     1972        vy_input->GetParameterValue(&vy1,gauss_1);
     1973        vy_input->GetParameterValue(&vy2,gauss_2);
     1974
     1975        mass_flux= rho_ice*length*( 
     1976                                (ONETHIRD*(h1-h2)*(vx1-vx2)+0.5*h2*(vx1-vx2)+0.5*(h1-h2)*vx2+h2*vx2)*normal[0]+
     1977                                (ONETHIRD*(h1-h2)*(vy1-vy2)+0.5*h2*(vy1-vy2)+0.5*(h1-h2)*vy2+h2*vy2)*normal[1]
     1978                                );
     1979
     1980        /*Process units: */
     1981        mass_flux=UnitConversion(mass_flux,IuToExtEnum,MassFluxEnum);
     1982
     1983        /*clean up and return:*/
     1984        delete gauss_1;
     1985        delete gauss_2;
     1986        return mass_flux;
     1987}
     1988/*}}}*/
     1989/*FUNCTION Tria::MaxAbsVx{{{1*/
     1990void  Tria::MaxAbsVx(double* pmaxabsvx, bool process_units){
     1991
     1992        /*Get maximum:*/
     1993        double maxabsvx=this->inputs->MaxAbs(VxEnum);
     1994
     1995        /*process units if requested: */
     1996        if(process_units) maxabsvx=UnitConversion(maxabsvx,IuToExtEnum,VxEnum);
     1997
     1998        /*Assign output pointers:*/
     1999        *pmaxabsvx=maxabsvx;
     2000}
     2001/*}}}*/
     2002/*FUNCTION Tria::MaxAbsVy{{{1*/
     2003void  Tria::MaxAbsVy(double* pmaxabsvy, bool process_units){
     2004
     2005        /*Get maximum:*/
     2006        double maxabsvy=this->inputs->MaxAbs(VyEnum);
     2007
     2008        /*process units if requested: */
     2009        if(process_units) maxabsvy=UnitConversion(maxabsvy,IuToExtEnum,VyEnum);
     2010
     2011        /*Assign output pointers:*/
     2012        *pmaxabsvy=maxabsvy;
     2013}
     2014/*}}}*/
     2015/*FUNCTION Tria::MaxAbsVz{{{1*/
     2016void  Tria::MaxAbsVz(double* pmaxabsvz, bool process_units){
     2017
     2018        /*Get maximum:*/
     2019        double maxabsvz=this->inputs->MaxAbs(VzEnum);
     2020
     2021        /*process units if requested: */
     2022        if(process_units) maxabsvz=UnitConversion(maxabsvz,IuToExtEnum,VyEnum);
     2023
     2024        /*Assign output pointers:*/
     2025        *pmaxabsvz=maxabsvz;
     2026}
     2027/*}}}*/
     2028/*FUNCTION Tria::MaxVel{{{1*/
     2029void  Tria::MaxVel(double* pmaxvel, bool process_units){
     2030
     2031        /*Get maximum:*/
     2032        double maxvel=this->inputs->Max(VelEnum);
     2033
     2034        /*process units if requested: */
     2035        if(process_units) maxvel=UnitConversion(maxvel,IuToExtEnum,VelEnum);
     2036
     2037        /*Assign output pointers:*/
     2038        *pmaxvel=maxvel;
     2039}
     2040/*}}}*/
     2041/*FUNCTION Tria::MaxVx{{{1*/
     2042void  Tria::MaxVx(double* pmaxvx, bool process_units){
     2043
     2044        /*Get maximum:*/
     2045        double maxvx=this->inputs->Max(VxEnum);
     2046
     2047        /*process units if requested: */
     2048        if(process_units) maxvx=UnitConversion(maxvx,IuToExtEnum,VxEnum);
     2049
     2050        /*Assign output pointers:*/
     2051        *pmaxvx=maxvx;
     2052}
     2053/*}}}*/
     2054/*FUNCTION Tria::MaxVy{{{1*/
     2055void  Tria::MaxVy(double* pmaxvy, bool process_units){
     2056
     2057        /*Get maximum:*/
     2058        double maxvy=this->inputs->Max(VyEnum);
     2059
     2060        /*process units if requested: */
     2061        if(process_units) maxvy=UnitConversion(maxvy,IuToExtEnum,VyEnum);
     2062
     2063        /*Assign output pointers:*/
     2064        *pmaxvy=maxvy;
     2065
     2066}
     2067/*}}}*/
     2068/*FUNCTION Tria::MaxVz{{{1*/
     2069void  Tria::MaxVz(double* pmaxvz, bool process_units){
     2070
     2071        /*Get maximum:*/
     2072        double maxvz=this->inputs->Max(VzEnum);
     2073
     2074        /*process units if requested: */
     2075        if(process_units) maxvz=UnitConversion(maxvz,IuToExtEnum,VzEnum);
     2076
     2077        /*Assign output pointers:*/
     2078        *pmaxvz=maxvz;
     2079}
     2080/*}}}*/
     2081/*FUNCTION Tria::MigrateGroundingLine{{{1*/
     2082void  Tria::MigrateGroundingLine(void){
     2083
     2084
     2085        double *values         = NULL;
     2086        double  h[3],s[3],b[3],ba[3];
     2087        double  bed_hydro;
     2088        double  rho_water,rho_ice,density;
     2089        int     isonshelf[3];
     2090        int     i;
     2091        bool    elementonshelf = false;
     2092
     2093
     2094        /*Recover info at the vertices: */
     2095        Input* surface_input =inputs->GetInput(SurfaceEnum); _assert_(surface_input);
     2096        Input* bed_input     =inputs->GetInput(BedEnum);     _assert_(bed_input);
     2097        if((surface_input->Enum()!=TriaVertexInputEnum) | (bed_input->Enum()!=TriaVertexInputEnum))_error_(" not supported yet for bed and surface interpolations not P1!");
     2098
     2099        GetParameterListOnVertices(&h[0],ThicknessEnum);
     2100        GetParameterListOnVertices(&s[0],SurfaceEnum);
     2101        GetParameterListOnVertices(&b[0],BedEnum);
     2102        GetParameterListOnVertices(&ba[0],BathymetryEnum);
     2103        for(i=0;i<3;i++){
     2104                isonshelf[i]=nodes[i]->IsOnShelf();
     2105                if((nodes[i]->Sid()+1)==36)printf("MigrateGroundingLine: El %i Node %i shelf status %i\n",this->Id(),nodes[i]->Sid()+1,isonshelf[i]);
     2106        }
     2107
     2108        /*material parameters: */
     2109        rho_water=matpar->GetRhoWater();
     2110        rho_ice=matpar->GetRhoIce();
     2111        density=rho_ice/rho_water;
     2112
     2113        /*go through vertices, and update inputs, considering them to be TriaVertex type: */
     2114        for(i=0;i<3;i++){
     2115                if (isonshelf[i]){
     2116                        /*This node is on the shelf. See if its bed is going under the bathymetry: */
     2117                        if(b[i]<=ba[i]){ //<= because Neff being 0 when b=ba, drag will be 0 anyway.
     2118                                /*The ice shelf is getting grounded, the thickness is the same, so just update the bed to stick to the bathymetry and elevate the surface accordingly: */
     2119                                b[i]=ba[i];
     2120                                s[i]=b[i]+h[i];
     2121                                if((nodes[i]->Sid()+1)==36)printf("MigrateGroundingLine: El %i Node %i grounding %g %g %g\n",this->Id(),nodes[i]->Sid()+1,b[i],ba[i],s[i]);
     2122                        }
     2123                        else{
     2124                                /*do nothing, we are still floating.*/
     2125                                if((nodes[i]->Sid()+1)==36)printf("MigrateGroundingLine: El %i Node %i still floating %g %g %g\n",this->Id(),nodes[i]->Sid()+1,b[i],ba[i],s[i]);
     2126                        }
     2127                }
     2128                else{
     2129                        /*This node is on the sheet, near the grounding line. See if wants to unground. To
     2130                         * do so, we compute the hydrostatic bed, and if it is > bathymetry, then we unground: */
     2131                        bed_hydro=-density*h[i];
     2132                        if (bed_hydro>ba[i]){
     2133                                /*We are now floating, bed and surface are determined from hydrostatic equilibrium: */
     2134                                s[i]=(1-density)*h[i];
     2135                                b[i]=-density*h[i];
     2136                                printf("%i\n",nodes[i]->Sid()+1);
     2137                                if((nodes[i]->Sid()+1)==36)printf("MigrateGroundingLine: El %i Node %i floating %g %g %g\n",this->Id(),nodes[i]->Sid()+1,b[i],ba[i],s[i]);
     2138                        }
     2139                        else{
     2140                                /*do nothing, we are still grounded.*/
     2141                                if((nodes[i]->Sid()+1)==36)printf("MigrateGroundingLine: El %i Node %i still grounded %g %g %g\n",this->Id(),nodes[i]->Sid()+1,b[i],ba[i],s[i]);
     2142                        }
     2143                }
     2144        }
     2145
     2146        /*Surface and bed are updated. Update inputs:*/
     2147        surface_input->GetValuesPtr(&values,NULL); for(i=0;i<3;i++)values[i]=s[i];
     2148        bed_input->GetValuesPtr(&values,NULL);     for(i=0;i<3;i++)values[i]=b[i];
     2149       
     2150        for(i=0;i<3;i++){
     2151                isonshelf[i]=nodes[i]->IsOnShelf();
     2152                if((nodes[i]->Sid()+1)==36)printf("MigrateGroundingLine: El %i Node %i second shelf status %i\n",this->Id(),nodes[i]->Sid()+1,isonshelf[i]);
     2153        }
     2154}
     2155/*}}}*/
     2156/*FUNCTION Tria::MinVel{{{1*/
     2157void  Tria::MinVel(double* pminvel, bool process_units){
     2158
     2159        /*Get minimum:*/
     2160        double minvel=this->inputs->Min(VelEnum);
     2161
     2162        /*process units if requested: */
     2163        if(process_units) minvel=UnitConversion(minvel,IuToExtEnum,VelEnum);
     2164
     2165        /*Assign output pointers:*/
     2166        *pminvel=minvel;
     2167}
     2168/*}}}*/
     2169/*FUNCTION Tria::MinVx{{{1*/
     2170void  Tria::MinVx(double* pminvx, bool process_units){
     2171
     2172        /*Get minimum:*/
     2173        double minvx=this->inputs->Min(VxEnum);
     2174
     2175        /*process units if requested: */
     2176        if(process_units) minvx=UnitConversion(minvx,IuToExtEnum,VxEnum);
     2177
     2178        /*Assign output pointers:*/
     2179        *pminvx=minvx;
     2180}
     2181/*}}}*/
     2182/*FUNCTION Tria::MinVy{{{1*/
     2183void  Tria::MinVy(double* pminvy, bool process_units){
     2184
     2185        /*Get minimum:*/
     2186        double minvy=this->inputs->Min(VyEnum);
     2187
     2188        /*process units if requested: */
     2189        if(process_units) minvy=UnitConversion(minvy,IuToExtEnum,VyEnum);
     2190
     2191        /*Assign output pointers:*/
     2192        *pminvy=minvy;
     2193}
     2194/*}}}*/
     2195/*FUNCTION Tria::MinVz{{{1*/
     2196void  Tria::MinVz(double* pminvz, bool process_units){
     2197
     2198        /*Get minimum:*/
     2199        double minvz=this->inputs->Min(VzEnum);
     2200
     2201        /*process units if requested: */
     2202        if(process_units) minvz=UnitConversion(minvz,IuToExtEnum,VzEnum);
     2203
     2204        /*Assign output pointers:*/
     2205        *pminvz=minvz;
     2206}
     2207/*}}}*/
     2208/*FUNCTION Tria::MyRank {{{1*/
     2209int    Tria::MyRank(void){
     2210        extern int my_rank;
     2211        return my_rank;
     2212}
     2213/*}}}*/
     2214/*FUNCTION Tria::NodalValue {{{1*/
     2215int    Tria::NodalValue(double* pvalue, int index, int natureofdataenum,bool process_units){
     2216
     2217        int i;
     2218        int found=0;
     2219        double value;
     2220        Input* data=NULL;
     2221        GaussTria *gauss                            = NULL;
     2222
     2223        /*First, serarch the input: */
     2224        data=inputs->GetInput(natureofdataenum);
     2225
     2226        /*figure out if we have the vertex id: */
     2227        found=0;
     2228        for(i=0;i<NUMVERTICES;i++){
     2229                if(index==nodes[i]->GetVertexId()){
     2230                        /*Do we have natureofdataenum in our inputs? :*/
     2231                        if(data){
     2232                                /*ok, we are good. retrieve value of input at vertex :*/
     2233                                gauss=new GaussTria(); gauss->GaussVertex(i);
     2234                                data->GetParameterValue(&value,gauss);
     2235                                found=1;
     2236                                break;
     2237                        }
     2238                }
     2239        }
     2240
     2241        if(found)*pvalue=value;
     2242        return found;
     2243}
     2244/*}}}*/
     2245/*FUNCTION Tria::PatchFill{{{1*/
     2246void  Tria::PatchFill(int* prow, Patch* patch){
     2247
     2248        int i,row;
     2249        int vertices_ids[3];
     2250
     2251        /*recover pointer: */
     2252        row=*prow;
     2253               
     2254        for(i=0;i<3;i++) vertices_ids[i]=nodes[i]->GetVertexId(); //vertices id start at column 3 of the patch.
     2255
     2256        for(i=0;i<this->results->Size();i++){
     2257                ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i);
     2258
     2259                /*For this result,fill the information in the Patch object (element id + vertices ids), and then hand
     2260                 *it to the result object, to fill the rest: */
     2261                patch->fillelementinfo(row,this->sid+1,vertices_ids,3);
     2262                elementresult->PatchFill(row,patch);
     2263
     2264                /*increment rower: */
     2265                row++;
     2266        }
     2267
     2268        /*Assign output pointers:*/
     2269        *prow=row;
     2270}
     2271/*}}}*/
     2272/*FUNCTION Tria::PatchSize{{{1*/
     2273void  Tria::PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes){
     2274
     2275        int     i;
     2276        int     numrows     = 0;
     2277        int     numnodes    = 0;
     2278        int     temp_numnodes    = 0;
     2279
     2280        /*Go through all the results objects, and update the counters: */
     2281        for (i=0;i<this->results->Size();i++){
     2282                ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i);
     2283                /*first, we have one more result: */
     2284                numrows++;
     2285                /*now, how many vertices and how many nodal values for this result? :*/
     2286                temp_numnodes=elementresult->NumberOfNodalValues(); //ask result object.
     2287                if(temp_numnodes>numnodes)numnodes=temp_numnodes;
     2288        }
     2289
     2290        /*Assign output pointers:*/
     2291        *pnumrows=numrows;
     2292        *pnumvertices=NUMVERTICES;
     2293        *pnumnodes=numnodes;
     2294}
     2295/*}}}*/
     2296/*FUNCTION Tria::PotentialSheetUngrounding{{{1*/
     2297void  Tria::PotentialSheetUngrounding(Vec potential_sheet_ungrounding){
     2298
     2299
     2300        double *values         = NULL;
     2301        double  h[3],s[3],b[3],ba[3];
     2302        double  bed_hydro;
     2303        double  rho_water,rho_ice,density;
     2304        int     i;
     2305        bool    elementonshelf = false;
     2306
     2307        /*Recover info at the vertices: */
     2308        Input* surface_input =inputs->GetInput(SurfaceEnum); _assert_(surface_input);
     2309        Input* bed_input     =inputs->GetInput(BedEnum);     _assert_(bed_input);
     2310        if((surface_input->Enum()!=TriaVertexInputEnum) | (bed_input->Enum()!=TriaVertexInputEnum))_error_(" not supported yet for bed and surface interpolations not P1!");
     2311
     2312        GetParameterListOnVertices(&h[0],ThicknessEnum);
     2313        GetParameterListOnVertices(&s[0],SurfaceEnum);
     2314        GetParameterListOnVertices(&b[0],BedEnum);
     2315        GetParameterListOnVertices(&ba[0],BathymetryEnum);
     2316
     2317        /*material parameters: */
     2318        rho_water=matpar->GetRhoWater();
     2319        rho_ice=matpar->GetRhoIce();
     2320        density=rho_ice/rho_water;
     2321
     2322        /*go through vertices, and figure out which ones are on the ice sheet, and want to unground: */
     2323        for(i=0;i<3;i++){
     2324                if (!nodes[i]->IsOnShelf()){
     2325                       
     2326                        /*This node is on the sheet, near the grounding line. See if wants to unground. To
     2327                         * do so, we compute the hydrostatic bed, and if it is > bathymetry, then we unground: */
     2328                        bed_hydro=-density*h[i];
     2329                        if (bed_hydro>ba[i]){
     2330                                /*ok, this node wants to unground. flag it: */
     2331                                VecSetValue(potential_sheet_ungrounding,nodes[i]->Sid(),1,INSERT_VALUES);
     2332                        }
     2333                }
     2334        }
     2335}
     2336/*}}}*/
     2337/*FUNCTION Tria::ProcessResultsUnits{{{1*/
     2338void  Tria::ProcessResultsUnits(void){
     2339
     2340        int i;
     2341
     2342        for(i=0;i<this->results->Size();i++){
     2343                ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i);
     2344                elementresult->ProcessUnits(this->parameters);
     2345        }
     2346}
     2347/*}}}*/
     2348/*FUNCTION Tria::RheologyBbarx{{{1*/
     2349double Tria::RheologyBbarx(void){
     2350
     2351        return this->matice->GetBbar();
     2352
     2353}
     2354/*}}}*/
     2355/*FUNCTION Tria::RequestedOutput{{{1*/
     2356void Tria::RequestedOutput(int output_enum,int step,double time){
     2357
     2358        if(IsInput(output_enum)){
     2359                /*just transfer this input to results, and we are done: */
     2360                InputToResult(output_enum,step,time);
     2361        }
     2362        else{
     2363                /*this input does not exist, compute it, and then transfer to results: */
     2364                switch(output_enum){
     2365                        default:
     2366                                /*do nothing, no need to derail the computation because one of the outputs requested cannot be found: */
     2367                                break;
     2368                }
     2369        }
     2370
     2371}
     2372/*}}}*/
     2373/*FUNCTION Tria::SetClone {{{1*/
     2374void  Tria::SetClone(int* minranks){
     2375
     2376        _error_("not implemented yet");
     2377}
     2378/*}}}1*/
     2379/*FUNCTION Tria::SetCurrentConfiguration {{{1*/
     2380void  Tria::SetCurrentConfiguration(Elements* elementsin, Loads* loadsin, DataSet* nodesin, Materials* materialsin, Parameters* parametersin){
     2381       
     2382        /*go into parameters and get the analysis_counter: */
     2383        int analysis_counter;
     2384        parametersin->FindParam(&analysis_counter,AnalysisCounterEnum);
     2385
     2386        /*Get Element type*/
     2387        this->element_type=this->element_type_list[analysis_counter];
     2388
     2389        /*Pick up nodes*/
     2390        if(this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp();
     2391        else this->nodes=NULL;
     2392
     2393}
     2394/*}}}*/
     2395/*FUNCTION Tria::ShelfSync{{{1*/
     2396void  Tria::ShelfSync(void){
     2397
     2398
     2399        double *values         = NULL;
     2400        double  h[3],s[3],b[3],ba[3];
     2401        double  bed_hydro;
     2402        double  rho_water,rho_ice,density;
     2403        int     i;
     2404        bool    elementonshelf = false;
     2405
     2406        /*melting rate at the grounding line: */
     2407        double  yts;
     2408        int     swap;
     2409        double  gl_melting_rate;
     2410
     2411        /*recover parameters: */
     2412        parameters->FindParam(&yts,ConstantsYtsEnum);
     2413        parameters->FindParam(&gl_melting_rate,GroundinglineMeltingRateEnum);
     2414
     2415        /*Recover info at the vertices: */
     2416        Input* surface_input =inputs->GetInput(SurfaceEnum); _assert_(surface_input);
     2417        Input* bed_input     =inputs->GetInput(BedEnum);     _assert_(bed_input);
     2418        if((surface_input->Enum()!=TriaVertexInputEnum) | (bed_input->Enum()!=TriaVertexInputEnum))_error_(" not supported yet for bed and surface interpolations not P1!");
     2419
     2420        GetParameterListOnVertices(&h[0],ThicknessEnum);
     2421        GetParameterListOnVertices(&s[0],SurfaceEnum);
     2422        GetParameterListOnVertices(&b[0],BedEnum);
     2423        GetParameterListOnVertices(&ba[0],BathymetryEnum);
     2424
     2425        /*material parameters: */
     2426        rho_water=matpar->GetRhoWater();
     2427        rho_ice=matpar->GetRhoIce();
     2428        density=rho_ice/rho_water;
     2429
     2430        /*go through vertices, and update inputs, considering them to be TriaVertex type: */
     2431        for(i=0;i<3;i++){
     2432                if(b[i]==ba[i]){
     2433                               
     2434                        nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,false));
     2435                        nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,true));
     2436                }
     2437                else{
     2438                        nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,true));
     2439                        nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,false));
     2440
     2441                }
     2442        }
     2443
     2444        /*Now, update  shelf status of element. An element can only be on shelf if all its nodes are on shelf: */
     2445        swap=0;
     2446        elementonshelf=false;
     2447        for(i=0;i<3;i++){
     2448                if(nodes[i]->IsOnShelf()){
     2449                        elementonshelf=true;
     2450                        break;
     2451                }
     2452        }
     2453        if(!this->IsOnShelf() && elementonshelf==true)swap=1;
     2454    this->inputs->AddInput(new BoolInput(MaskElementonfloatingiceEnum,elementonshelf));
     2455       
     2456    /*If this element just  became ungrounded, set its basal melting rate at 50 m/yr:*/
     2457        if(swap){
     2458                Input* basal_melting_rate_input     =inputs->GetInput(BasalforcingsMeltingRateEnum);     _assert_(basal_melting_rate_input);
     2459                basal_melting_rate_input->GetValuesPtr(&values,NULL); for(i=0;i<3;i++)values[i]=gl_melting_rate/yts;
     2460        }
     2461
     2462
     2463}
     2464/*}}}*/
     2465/*FUNCTION Tria::SoftMigration{{{1*/
     2466void  Tria::SoftMigration(double* sheet_ungrounding){
     2467
     2468
     2469        double *values         = NULL;
     2470        double  h[3],s[3],b[3],ba[3];
     2471        double  bed_hydro;
     2472        double  rho_water,rho_ice,density;
     2473        int     i;
     2474        bool    elementonshelf = false;
     2475
     2476        /*Recover info at the vertices: */
     2477        Input* surface_input =inputs->GetInput(SurfaceEnum); _assert_(surface_input);
     2478        Input* bed_input     =inputs->GetInput(BedEnum);     _assert_(bed_input);
     2479        if((surface_input->Enum()!=TriaVertexInputEnum) | (bed_input->Enum()!=TriaVertexInputEnum))_error_(" not supported yet for bed and surface interpolations not P1!");
     2480
     2481        GetParameterListOnVertices(&h[0],ThicknessEnum);
     2482        GetParameterListOnVertices(&s[0],SurfaceEnum);
     2483        GetParameterListOnVertices(&b[0],BedEnum);
     2484        GetParameterListOnVertices(&ba[0],BathymetryEnum);
     2485
     2486        /*material parameters: */
     2487        rho_water=matpar->GetRhoWater();
     2488        rho_ice=matpar->GetRhoIce();
     2489        density=rho_ice/rho_water;
     2490       
     2491        /*go through vertices, and update inputs, considering them to be TriaVertex type: */
     2492        for(i=0;i<3;i++){
     2493                if (nodes[i]->IsOnShelf()){
     2494                        /*This node is on the shelf. See if its bed is going under the bathymetry: */
     2495                        if(b[i]<=ba[i]){ //<= because Neff being 0 when b=ba, drag will be 0 anyway.
     2496                                /*The ice shelf is getting grounded, the thickness is the same, so just update the bed to stick to the bathymetry and elevate the surface accordingly: */
     2497                                b[i]=ba[i];
     2498                                s[i]=b[i]+h[i];
     2499                        }
     2500                        else{
     2501                                /*do nothing, we are still floating.*/
     2502                        }
     2503                }
     2504                else{
     2505                        /*This node is on the sheet, near the grounding line. See if wants to unground. To
     2506                         * do so, we compute the hydrostatic bed, and if it is > bathymetry, then we unground: */
     2507                        bed_hydro=-density*h[i];
     2508                        if (bed_hydro>ba[i]){
     2509
     2510                                /*Now, are we connected to the grounding line, if so, go forward, otherwise, bail out: */
     2511                                if(sheet_ungrounding[nodes[i]->Sid()]){
     2512                                        /*We are now floating, bed and surface are determined from hydrostatic equilibrium: */
     2513                                        s[i]=(1-density)*h[i];
     2514                                        b[i]=-density*h[i];
     2515                                }
     2516                        }
     2517                        else{
     2518                                /*do nothing, we are still grounded.*/
     2519                        }
     2520                }
     2521        }
     2522
     2523        /*Surface and bed are updated. Update inputs:*/   
     2524        surface_input->GetValuesPtr(&values,NULL); for(i=0;i<3;i++)values[i]=s[i];
     2525        bed_input->GetValuesPtr(&values,NULL);     for(i=0;i<3;i++)values[i]=b[i];
     2526
     2527        }
     2528/*}}}*/
     2529/*FUNCTION Tria::SurfaceArea {{{1*/
     2530double Tria::SurfaceArea(void){
     2531
     2532        int    i;
     2533        double S;
     2534        double normal[3];
     2535        double v13[3],v23[3];
     2536        double xyz_list[NUMVERTICES][3];
     2537
     2538        /*If on water, return 0: */
     2539        if(IsOnWater())return 0;
     2540
     2541        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     2542
     2543        for (i=0;i<3;i++){
     2544                v13[i]=xyz_list[0][i]-xyz_list[2][i];
     2545                v23[i]=xyz_list[1][i]-xyz_list[2][i];
     2546        }
     2547
     2548        normal[0]=v13[1]*v23[2]-v13[2]*v23[1];
     2549        normal[1]=v13[2]*v23[0]-v13[0]*v23[2];
     2550        normal[2]=v13[0]*v23[1]-v13[1]*v23[0];
     2551
     2552        S = 0.5 * sqrt(pow(normal[0],(double)2)+pow(normal[1],(double)2)+pow(normal[2],(double)2));
     2553
     2554        /*Return: */
     2555        return S;
     2556}
     2557/*}}}*/
     2558/*FUNCTION Tria::SurfaceNormal{{{1*/
     2559void Tria::SurfaceNormal(double* surface_normal, double xyz_list[3][3]){
     2560
     2561        int i;
     2562        double v13[3],v23[3];
     2563        double normal[3];
     2564        double normal_norm;
     2565
     2566        for (i=0;i<3;i++){
     2567                v13[i]=xyz_list[0][i]-xyz_list[2][i];
     2568                v23[i]=xyz_list[1][i]-xyz_list[2][i];
     2569        }
     2570
     2571        normal[0]=v13[1]*v23[2]-v13[2]*v23[1];
     2572        normal[1]=v13[2]*v23[0]-v13[0]*v23[2];
     2573        normal[2]=v13[0]*v23[1]-v13[1]*v23[0];
     2574
     2575        normal_norm=sqrt( pow(normal[0],(double)2)+pow(normal[1],(double)2)+pow(normal[2],(double)2) );
     2576
     2577        *(surface_normal)=normal[0]/normal_norm;
     2578        *(surface_normal+1)=normal[1]/normal_norm;
     2579        *(surface_normal+2)=normal[2]/normal_norm;
     2580}
     2581/*}}}*/
     2582/*FUNCTION Tria::TimeAdapt{{{1*/
     2583double  Tria::TimeAdapt(void){
     2584
     2585        /*intermediary: */
     2586        int    i;
     2587        double C,dt;
     2588        double dx,dy;
     2589        double maxx,minx;
     2590        double maxy,miny;
     2591        double maxabsvx,maxabsvy;
     2592        double xyz_list[NUMVERTICES][3];
     2593
     2594        /*get CFL coefficient:*/
     2595        this->parameters->FindParam(&C,TimesteppingCflCoefficientEnum);
     2596
     2597        /*Get for Vx and Vy, the max of abs value: */
     2598        this->MaxAbsVx(&maxabsvx,false);
     2599        this->MaxAbsVy(&maxabsvy,false);
     2600
     2601        /* Get node coordinates and dof list: */
     2602        GetVerticesCoordinates(&xyz_list[0][0], this->nodes, NUMVERTICES);
     2603
     2604        minx=xyz_list[0][0];
     2605        maxx=xyz_list[0][0];
     2606        miny=xyz_list[0][1];
     2607        maxy=xyz_list[0][1];
     2608
     2609        for(i=1;i<NUMVERTICES;i++){
     2610                if (xyz_list[i][0]<minx)minx=xyz_list[i][0];
     2611                if (xyz_list[i][0]>maxx)maxx=xyz_list[i][0];
     2612                if (xyz_list[i][1]<miny)miny=xyz_list[i][1];
     2613                if (xyz_list[i][1]>maxy)maxy=xyz_list[i][1];
     2614        }
     2615        dx=maxx-minx;
     2616        dy=maxy-miny;
     2617
     2618        /*CFL criterion: */
     2619        dt=C/(maxabsvy/dx+maxabsvy/dy);
     2620
     2621        return dt;
     2622}
     2623/*}}}*/
     2624/*FUNCTION Tria::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type){{{1*/
     2625void Tria::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type){ //i is the element index
     2626
     2627        /*Intermediaries*/
     2628        int    i,j;
     2629        int    tria_node_ids[3];
     2630        int    tria_vertex_ids[3];
     2631        int    tria_type;
     2632        double nodeinputs[3];
     2633        double yts;
     2634        int    progstabilization,balancestabilization;
     2635        bool   dakota_analysis;
     2636
     2637        /*Checks if debuging*/
     2638        /*{{{2*/
     2639        _assert_(iomodel->Data(MeshElementsEnum));
     2640        /*}}}*/
     2641
     2642        /*Fetch parameters: */
     2643        iomodel->Constant(&yts,ConstantsYtsEnum);
     2644        iomodel->Constant(&progstabilization,PrognosticStabilizationEnum);
     2645        iomodel->Constant(&balancestabilization,BalancethicknessStabilizationEnum);
     2646        iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum);
     2647
     2648        /*Recover element type*/
     2649        if ((analysis_type==PrognosticAnalysisEnum && progstabilization==3) || (analysis_type==BalancethicknessAnalysisEnum && balancestabilization==3)){
     2650                /*P1 Discontinuous Galerkin*/
     2651                tria_type=P1DGEnum;
     2652        }
     2653        else{
     2654                /*P1 Continuous Galerkin*/
     2655                tria_type=P1Enum;
     2656        }
     2657        this->SetElementType(tria_type,analysis_counter);
     2658
     2659        /*Recover vertices ids needed to initialize inputs*/
     2660        for(i=0;i<3;i++){
     2661                tria_vertex_ids[i]=(int)iomodel->Data(MeshElementsEnum)[3*index+i]; //ids for vertices are in the elements array from Matlab
     2662        }
     2663
     2664        /*Recover nodes ids needed to initialize the node hook.*/
     2665        if (tria_type==P1DGEnum){
     2666                /*Discontinuous Galerkin*/
     2667                tria_node_ids[0]=iomodel->nodecounter+3*index+1;
     2668                tria_node_ids[1]=iomodel->nodecounter+3*index+2;
     2669                tria_node_ids[2]=iomodel->nodecounter+3*index+3;
     2670        }
     2671        else{
     2672                /*Continuous Galerkin*/
     2673                for(i=0;i<3;i++){
     2674                        tria_node_ids[i]=iomodel->nodecounter+(int)*(iomodel->Data(MeshElementsEnum)+3*index+i); //ids for vertices are in the elements array from Matlab
     2675                }
     2676        }
     2677
     2678        /*hooks: */
     2679        this->SetHookNodes(tria_node_ids,analysis_counter); this->nodes=NULL; //set hook to nodes, for this analysis type
     2680
     2681        /*Fill with IoModel*/
     2682        this->InputUpdateFromIoModel(index,iomodel);
     2683
     2684        /*Defaults if not provided in iomodel*/
     2685        switch(analysis_type){
     2686
     2687                case DiagnosticHorizAnalysisEnum:
     2688
     2689                        /*default vx,vy and vz: either observation or 0 */
     2690                        if(!iomodel->Data(VxEnum)){
     2691                                if (iomodel->Data(InversionVxObsEnum)) for(i=0;i<3;i++)nodeinputs[i]=iomodel->Data(InversionVxObsEnum)[tria_vertex_ids[i]-1]/yts;
     2692                                else                 for(i=0;i<3;i++)nodeinputs[i]=0;
     2693                                this->inputs->AddInput(new TriaVertexInput(VxEnum,nodeinputs));
     2694                                this->inputs->AddInput(new TriaVertexInput(VxPicardEnum,nodeinputs));
     2695                                if(dakota_analysis) this->inputs->AddInput(new TriaVertexInput(QmuVxEnum,nodeinputs));
     2696                        }
     2697                        if(!iomodel->Data(VyEnum)){
     2698                                if (iomodel->Data(InversionVyObsEnum)) for(i=0;i<3;i++)nodeinputs[i]=iomodel->Data(InversionVyObsEnum)[tria_vertex_ids[i]-1]/yts;
     2699                                else                 for(i=0;i<3;i++)nodeinputs[i]=0;
     2700                                this->inputs->AddInput(new TriaVertexInput(VyEnum,nodeinputs));
     2701                                this->inputs->AddInput(new TriaVertexInput(VyPicardEnum,nodeinputs));
     2702                                if(dakota_analysis) this->inputs->AddInput(new TriaVertexInput(QmuVyEnum,nodeinputs));
     2703                        }
     2704                        if(!iomodel->Data(VzEnum)){
     2705                                if (iomodel->Data(InversionVzObsEnum)) for(i=0;i<3;i++)nodeinputs[i]=iomodel->Data(InversionVzObsEnum)[tria_vertex_ids[i]-1]/yts;
     2706                                else                 for(i=0;i<3;i++)nodeinputs[i]=0;
     2707                                this->inputs->AddInput(new TriaVertexInput(VzEnum,nodeinputs));
     2708                                this->inputs->AddInput(new TriaVertexInput(VzPicardEnum,nodeinputs));
     2709                                if(dakota_analysis) this->inputs->AddInput(new TriaVertexInput(QmuVzEnum,nodeinputs));
     2710                        }
     2711                        if(!iomodel->Data(PressureEnum)){
     2712                                for(i=0;i<3;i++)nodeinputs[i]=0;
     2713                                if(dakota_analysis){
     2714                                        this->inputs->AddInput(new TriaVertexInput(PressureEnum,nodeinputs));
     2715                                        this->inputs->AddInput(new TriaVertexInput(QmuPressureEnum,nodeinputs));
     2716                                }
     2717                        }
     2718                        break;
     2719
     2720                default:
     2721                        /*No update for other solution types*/
     2722                        break;
     2723
     2724        }
     2725
     2726        //this->parameters: we still can't point to it, it may not even exist. Configure will handle this.
     2727        this->parameters=NULL;
     2728}
     2729/*}}}*/
     2730/*FUNCTION Tria::UpdateShelfStatus{{{1*/
     2731int  Tria::UpdateShelfStatus(Vec new_shelf_nodes){
     2732
     2733        int     i;
     2734
     2735        double  h[3];
     2736        double  s[3];
     2737        double  b[3];
     2738        double  ba[3];
     2739        Input  *surface_input  = NULL;
     2740        Input  *bed_input      = NULL;
     2741        double  rho_water;
     2742        double  rho_ice;
     2743        double  density;
     2744        bool    elementonshelf = false;
     2745        int     flipped=0;
     2746        int     shelfstatus[3];
     2747
     2748
     2749        /*Recover info at the vertices: */
     2750        surface_input   =inputs->GetInput(SurfaceEnum);   _assert_(surface_input);
     2751        bed_input   =inputs->GetInput(BedEnum);   _assert_(bed_input);
     2752        if((surface_input->Enum()!=TriaVertexInputEnum) | (bed_input->Enum()!=TriaVertexInputEnum))_error_(" not supported yet for bed and surface interpolations not P1!");
     2753
     2754        GetParameterListOnVertices(&h[0],ThicknessEnum);
     2755        GetParameterListOnVertices(&s[0],SurfaceEnum);
     2756        GetParameterListOnVertices(&b[0],BedEnum);
     2757        GetParameterListOnVertices(&ba[0],BathymetryEnum);
     2758
     2759        /*material parameters: */
     2760        rho_water=matpar->GetRhoWater();
     2761        rho_ice=matpar->GetRhoIce();
     2762        density=rho_ice/rho_water;
     2763
     2764        /*Initialize current status of nodes: */
     2765        for(i=0;i<3;i++){
     2766                shelfstatus[i]=nodes[i]->IsOnShelf();
     2767                if((nodes[i]->Sid()+1)==36) printf("UpdateShelfStatus: El %i Node %i shelf status %i\n",this->Id(),nodes[i]->Sid()+1,shelfstatus[i]);
     2768        }
     2769       
     2770        /*go through vertices, and figure out if they are grounded or not, then update their status: */
     2771        flipped=0;
     2772        for(i=0;i<3;i++){
     2773                if(b[i]<=ba[i]){ //the = will lead to oscillations.
     2774                        nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,false));
     2775                        nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,true));
     2776                        if(shelfstatus[i]){
     2777                                flipped++;
     2778                                if((nodes[i]->Sid()+1)==36)printf("UpdateShelfStatus: El %i Node %i grounding %g %g %g\n",this->Id(),nodes[i]->Sid()+1,b[i],ba[i],s[i]);
     2779                        }
     2780                        VecSetValue(new_shelf_nodes,nodes[i]->Sid(),0,INSERT_VALUES);
     2781                }
     2782                else{
     2783                        nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,true));
     2784                        nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,false));
     2785                        if(!shelfstatus[i]){
     2786                                flipped++;
     2787                                if((nodes[i]->Sid()+1)==36)printf("UpdateShelfStatus: El %i Node %i floating %g %g %g\n",this->Id(),nodes[i]->Sid()+1,b[i],ba[i],s[i]);
     2788                        }
     2789                        VecSetValue(new_shelf_nodes,nodes[i]->Sid(),1,INSERT_VALUES);
     2790                }
     2791        }
     2792
     2793        /*Now, update  shelf status of element. An element can only be on shelf if all its nodes are on shelf: */
     2794        elementonshelf=false;
     2795        for(i=0;i<3;i++){
     2796                if(nodes[i]->IsOnShelf()){
     2797                        elementonshelf=true;
     2798                        break;
     2799                }
     2800        }
     2801    this->inputs->AddInput(new BoolInput(MaskElementonfloatingiceEnum,elementonshelf));
     2802
     2803         return flipped;
     2804}
     2805/*}}}*/
     2806/*FUNCTION Tria::UpdateShelfFlags{{{1*/
     2807void Tria::UpdateShelfFlags(double* new_shelf_nodes){
     2808
     2809        /*go through vertices, and update the status of MaskVertexonfloatingiceEnum and MaskVertexongroundediceEnum flags: */
     2810        bool flag;
     2811        int  i;
     2812        double  h[3];
     2813        double  s[3];
     2814        double  b[3];
     2815        double  ba[3];
     2816
     2817        GetParameterListOnVertices(&h[0],ThicknessEnum);
     2818        GetParameterListOnVertices(&s[0],SurfaceEnum);
     2819        GetParameterListOnVertices(&b[0],BedEnum);
     2820        GetParameterListOnVertices(&ba[0],BathymetryEnum);
     2821
     2822        for(i=0;i<3;i++){
     2823                flag=(bool)new_shelf_nodes[nodes[i]->Sid()];
     2824                nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,flag));
     2825                nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,!flag));
     2826        }
     2827}
     2828/*}}}*/
     2829/*FUNCTION Tria::UpdatePotentialSheetUngrounding{{{1*/
     2830int Tria::UpdatePotentialSheetUngrounding(double* potential_sheet_ungrounding,Vec vec_nodes_on_iceshelf,double* nodes_on_iceshelf){
     2831
     2832        /*intermediary: */
     2833        int i;
     2834        int nflipped=0;
     2835
     2836        /*Ok, go through our 3 nodes, and whoever is on the potential_sheet_ungrounding, ends up in nodes_on_iceshelf: */
     2837        for(i=0;i<3;i++){
     2838                if (potential_sheet_ungrounding[nodes[i]->Sid()]){
     2839                        VecSetValue(vec_nodes_on_iceshelf,nodes[i]->Sid(),1,INSERT_VALUES);
     2840               
     2841                        /*Figure out if we flipped: */
     2842                        if (potential_sheet_ungrounding[nodes[i]->Sid()] != nodes_on_iceshelf[nodes[i]->Sid()])nflipped++;
     2843                }
     2844        }
     2845        return nflipped;
     2846}
     2847/*}}}*/
     2848
     2849#ifdef _HAVE_DIAGNOSTIC_
     2850/*FUNCTION Tria::CreateKMatrixDiagnosticMacAyeal {{{1*/
     2851ElementMatrix* Tria::CreateKMatrixDiagnosticMacAyeal(void){
     2852
     2853        /*compute all stiffness matrices for this element*/
     2854        ElementMatrix* Ke1=CreateKMatrixDiagnosticMacAyealViscous();
     2855        ElementMatrix* Ke2=CreateKMatrixDiagnosticMacAyealFriction();
     2856        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
     2857       
     2858        /*clean-up and return*/
     2859        delete Ke1;
     2860        delete Ke2;
     2861        return Ke;
     2862}
     2863/*}}}*/
     2864/*FUNCTION Tria::CreateKMatrixDiagnosticMacAyealViscous{{{1*/
     2865ElementMatrix* Tria::CreateKMatrixDiagnosticMacAyealViscous(void){
     2866
     2867        /*Constants*/
     2868        const int  numdof=NDOF2*NUMVERTICES;
     2869
     2870        /*Intermediaries*/
     2871        int        i,j,ig;
     2872        double     xyz_list[NUMVERTICES][3];
     2873        double     viscosity,newviscosity,oldviscosity;
     2874        double     viscosity_overshoot,thickness,Jdet;
     2875        double     epsilon[3],oldepsilon[3];    /* epsilon=[exx,eyy,exy];    */
     2876        double     B[3][numdof];
     2877        double     Bprime[3][numdof];
     2878        double     D[3][3]   = {0.0};
     2879        double     D_scalar;
     2880        GaussTria *gauss = NULL;
     2881
     2882        /*Initialize Element matrix*/
     2883        ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum);
     2884
     2885        /*Retrieve all inputs and parameters*/
     2886        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     2887        Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
     2888        Input* vx_input=inputs->GetInput(VxEnum);               _assert_(vx_input);
     2889        Input* vy_input=inputs->GetInput(VyEnum);               _assert_(vy_input);
     2890        Input* vxold_input=inputs->GetInput(VxPicardEnum);      _assert_(vxold_input);
     2891        Input* vyold_input=inputs->GetInput(VyPicardEnum);      _assert_(vyold_input);
     2892        this->parameters->FindParam(&viscosity_overshoot,DiagnosticViscosityOvershootEnum);
     2893
     2894        /* Start  looping on the number of gaussian points: */
     2895        gauss=new GaussTria(2);
     2896        for (ig=gauss->begin();ig<gauss->end();ig++){
     2897
     2898                gauss->GaussPoint(ig);
     2899
     2900                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
     2901                GetBMacAyeal(&B[0][0], &xyz_list[0][0], gauss);
     2902                GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss);
     2903
     2904                this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
     2905                this->GetStrainRate2d(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input);
     2906                matice->GetViscosity2d(&viscosity, &epsilon[0]);
     2907                matice->GetViscosity2d(&oldviscosity, &oldepsilon[0]);
     2908                thickness_input->GetParameterValue(&thickness, gauss);
     2909
     2910                newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
     2911                D_scalar=2*newviscosity*thickness*gauss->weight*Jdet;
     2912                for (i=0;i<3;i++) D[i][i]=D_scalar;
     2913
     2914                TripleMultiply(&B[0][0],3,numdof,1,
     2915                                        &D[0][0],3,3,0,
     2916                                        &Bprime[0][0],3,numdof,0,
     2917                                        &Ke->values[0],1);
     2918        }
     2919
     2920        /*Clean up and return*/
     2921        delete gauss;
     2922        return Ke;
     2923}
     2924/*}}}*/
     2925/*FUNCTION Tria::CreateKMatrixDiagnosticMacAyealFriction {{{1*/
     2926ElementMatrix* Tria::CreateKMatrixDiagnosticMacAyealFriction(void){
     2927
     2928        /*Constants*/
     2929        const int  numdof=NDOF2*NUMVERTICES;
     2930
     2931        /*Intermediaries*/
     2932        int        i,j,ig;
     2933        int        analysis_type;
     2934        double     MAXSLOPE  = .06; // 6 %
     2935        double     MOUNTAINKEXPONENT = 10;
     2936        double     slope_magnitude,alpha2;
     2937        double     Jdet;
     2938        double     L[2][numdof];
     2939        double     DL[2][2]  = {{ 0,0 },{0,0}};
     2940        double     DL_scalar;
     2941        double     slope[2]  = {0.0,0.0};
     2942        double     xyz_list[NUMVERTICES][3];
     2943        Friction  *friction = NULL;
     2944        GaussTria *gauss    = NULL;
     2945
     2946        /*Initialize Element matrix and return if necessary*/
     2947        if(IsOnShelf()) return NULL;
     2948        ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum);
     2949
     2950        /*Retrieve all inputs and parameters*/
     2951        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     2952        Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input);
     2953        Input* vx_input=inputs->GetInput(VxEnum);           _assert_(vx_input);
     2954        Input* vy_input=inputs->GetInput(VyEnum);           _assert_(vy_input);
     2955        Input* vz_input=inputs->GetInput(VzEnum);           _assert_(vz_input);
     2956        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
     2957
     2958        /*build friction object, used later on: */
     2959        friction=new Friction("2d",inputs,matpar,analysis_type);
     2960
     2961        /* Start  looping on the number of gaussian points: */
     2962        gauss=new GaussTria(2);
     2963        for (ig=gauss->begin();ig<gauss->end();ig++){
     2964
     2965                gauss->GaussPoint(ig);
     2966
     2967                // If we have a slope > 6% for this element,  it means  we are on a mountain. In this particular case,
     2968                //velocity should be = 0. To achieve this result, we set alpha2_list to a very high value: */
     2969                surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
     2970                slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2));
     2971                if(slope_magnitude>MAXSLOPE) alpha2=pow((double)10,MOUNTAINKEXPONENT);
     2972                else friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum);
     2973
     2974                GetL(&L[0][0], &xyz_list[0][0], gauss,NDOF2);
     2975                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
     2976                DL_scalar=alpha2*gauss->weight*Jdet;
     2977                for (i=0;i<2;i++) DL[i][i]=DL_scalar;
     2978               
     2979                TripleMultiply( &L[0][0],2,numdof,1,
     2980                                        &DL[0][0],2,2,0,
     2981                                        &L[0][0],2,numdof,0,
     2982                                        &Ke->values[0],1);
     2983        }
     2984
     2985        /*Clean up and return*/
     2986        delete gauss;
     2987        delete friction;
     2988        return Ke;
     2989}
     2990/*}}}*/
     2991/*FUNCTION Tria::CreateKMatrixDiagnosticHutter{{{1*/
     2992ElementMatrix* Tria::CreateKMatrixDiagnosticHutter(void){
     2993
     2994        /*Intermediaries*/
     2995        const int numdof=NUMVERTICES*NDOF2;
     2996        int    i,connectivity;
     2997
     2998        /*Initialize Element matrix*/
     2999        ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
     3000
     3001        /*Create Element matrix*/
     3002        for(i=0;i<NUMVERTICES;i++){
     3003                connectivity=nodes[i]->GetConnectivity();
     3004                Ke->values[(2*i)*numdof  +(2*i)  ]=1/(double)connectivity;
     3005                Ke->values[(2*i+1)*numdof+(2*i+1)]=1/(double)connectivity;
     3006        }
     3007
     3008        /*Clean up and return*/
     3009        return Ke;
     3010}
     3011/*}}}*/
     3012/*FUNCTION Tria::CreatePVectorDiagnosticMacAyeal {{{1*/
     3013ElementVector* Tria::CreatePVectorDiagnosticMacAyeal(){
     3014
     3015        /*Constants*/
     3016        const int    numdof=NDOF2*NUMVERTICES;
     3017
     3018        /*Intermediaries */
     3019        int            i,j,ig;
     3020        double         driving_stress_baseline,thickness;
     3021        double         Jdet;
     3022        double         xyz_list[NUMVERTICES][3];
     3023        double         slope[2];
     3024        double         basis[3];
     3025        double         pe_g_gaussian[numdof];
     3026        GaussTria*     gauss=NULL;
     3027
     3028        /*Initialize Element vector*/
     3029        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum);
     3030
     3031        /*Retrieve all inputs and parameters*/
     3032        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     3033        Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
     3034        Input* surface_input=inputs->GetInput(SurfaceEnum);     _assert_(surface_input);
     3035        Input* drag_input=inputs->GetInput(FrictionCoefficientEnum);_assert_(drag_input);
     3036
     3037        /* Start  looping on the number of gaussian points: */
     3038        gauss=new GaussTria(2);
     3039        for(ig=gauss->begin();ig<gauss->end();ig++){
     3040
     3041                gauss->GaussPoint(ig);
     3042
     3043                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
     3044                GetNodalFunctions(basis, gauss);
     3045
     3046                thickness_input->GetParameterValue(&thickness,gauss);
     3047                surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
     3048                driving_stress_baseline=matpar->GetRhoIce()*matpar->GetG()*thickness;
     3049
     3050                /*Build pe_g_gaussian vector: */
     3051                for (i=0;i<NUMVERTICES;i++){
     3052                        for (j=0;j<NDOF2;j++){
     3053                                pe->values[i*NDOF2+j]+=-driving_stress_baseline*slope[j]*Jdet*gauss->weight*basis[i];
     3054                        }
     3055                }
     3056        }
     3057
     3058        /*Clean up and return*/
     3059        delete gauss;
     3060        return pe;
     3061}
     3062/*}}}*/
     3063/*FUNCTION Tria::CreatePVectorDiagnosticHutter{{{1*/
     3064ElementVector* Tria::CreatePVectorDiagnosticHutter(void){
     3065
     3066        /*Intermediaries */
     3067        int        i,connectivity;
     3068        double     constant_part,ub,vb;
     3069        double     rho_ice,gravity,n,B;
     3070        double     slope2,thickness;
     3071        double     slope[2];
     3072        GaussTria* gauss=NULL;
     3073
     3074        /*Initialize Element vector*/
     3075        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
     3076
     3077        /*Retrieve all inputs and parameters*/
     3078        rho_ice=matpar->GetRhoIce();
     3079        gravity=matpar->GetG();
     3080        n=matice->GetN();
     3081        B=matice->GetBbar();
     3082        Input* slopex_input=inputs->GetInput(SurfaceSlopeXEnum); _assert_(slopex_input);
     3083        Input* slopey_input=inputs->GetInput(SurfaceSlopeYEnum); _assert_(slopey_input);
     3084        Input* thickness_input=inputs->GetInput(ThicknessEnum);  _assert_(thickness_input);
     3085
     3086        /*Spawn 3 sing elements: */
     3087        gauss=new GaussTria();
     3088        for(i=0;i<NUMVERTICES;i++){
     3089
     3090                gauss->GaussVertex(i);
     3091
     3092                connectivity=nodes[i]->GetConnectivity();
     3093
     3094                thickness_input->GetParameterValue(&thickness,gauss);
     3095                slopex_input->GetParameterValue(&slope[0],gauss);
     3096                slopey_input->GetParameterValue(&slope[1],gauss);
     3097                slope2=pow(slope[0],2)+pow(slope[1],2);
     3098
     3099                constant_part=-2*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2));
     3100
     3101                ub=-1.58*pow((double)10.0,(double)-10.0)*rho_ice*gravity*thickness*slope[0];
     3102                vb=-1.58*pow((double)10.0,(double)-10.0)*rho_ice*gravity*thickness*slope[1];
     3103
     3104                pe->values[2*i]  =(ub-2.0*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2.0))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[0])/(double)connectivity;
     3105                pe->values[2*i+1]=(vb-2.0*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2.0))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[1])/(double)connectivity;
     3106        }
     3107
     3108        /*Clean up and return*/
     3109        delete gauss;
     3110        return pe;
    20043111}
    20053112/*}}}*/
    … …  
    20823189}
    20833190/*}}}*/
    2084 /*FUNCTION Tria::GetSolutionFromInputsHydrology{{{1*/
    2085 void  Tria::GetSolutionFromInputsHydrology(Vec solution){
    2086 
    2087         const int    numdof=NDOF1*NUMVERTICES;
    2088 
    2089         int i,dummy;
    2090         int*         doflist=NULL;
    2091         double       watercolumn;
    2092         double       values[numdof];
    2093         GaussTria*   gauss=NULL;
    2094 
     3191/*FUNCTION Tria::InputUpdateFromSolutionDiagnosticHoriz {{{1*/
     3192void  Tria::InputUpdateFromSolutionDiagnosticHoriz(double* solution){
     3193       
     3194        const int numdof=NDOF2*NUMVERTICES;
     3195
     3196        int       i;
     3197        int       dummy;
     3198        int*      doflist=NULL;
     3199        double    rho_ice,g;
     3200        double    values[numdof];
     3201        double    vx[NUMVERTICES];
     3202        double    vy[NUMVERTICES];
     3203        double    vz[NUMVERTICES];
     3204        double    vel[NUMVERTICES];
     3205        double    pressure[NUMVERTICES];
     3206        double    thickness[NUMVERTICES];
     3207       
    20953208        /*Get dof list: */
    20963209        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    20973210
    2098         /*Get inputs*/
    2099         Input* watercolumn_input=inputs->GetInput(WatercolumnEnum); _assert_(watercolumn_input);
    2100 
    2101         /*Ok, we have watercolumn values, fill in watercolumn array: */
    2102         /*P1 element only for now*/
    2103         gauss=new GaussTria();
     3211        /*Use the dof list to index into the solution vector: */
     3212        for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
     3213
     3214        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    21043215        for(i=0;i<NUMVERTICES;i++){
    2105 
    2106                 gauss->GaussVertex(i);
    2107 
    2108                 /*Recover watercolumn*/
    2109                 watercolumn_input->GetParameterValue(&watercolumn,gauss);
    2110                 values[i]=watercolumn;
    2111         }
    2112 
    2113         VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES);
     3216                vx[i]=values[i*NDOF2+0];
     3217                vy[i]=values[i*NDOF2+1];
     3218
     3219                /*Check solution*/
     3220                if(isnan(vx[i])) _error_("NaN found in solution vector");
     3221                if(isnan(vy[i])) _error_("NaN found in solution vector");
     3222        }
     3223
     3224        /*Get Vz and compute vel*/
     3225        GetParameterListOnVertices(&vz[0],VzEnum,0);
     3226        for(i=0;i<NUMVERTICES;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
     3227
     3228        /*For pressure: we have not computed pressure in this analysis, for this element. We are in 2D,
     3229         *so the pressure is just the pressure at the bedrock: */
     3230        rho_ice=matpar->GetRhoIce();
     3231        g=matpar->GetG();
     3232        GetParameterListOnVertices(&thickness[0],ThicknessEnum);
     3233        for(i=0;i<NUMVERTICES;i++) pressure[i]=rho_ice*g*thickness[i];
     3234
     3235        /*Now, we have to move the previous Vx and Vy inputs  to old
     3236         * status, otherwise, we'll wipe them off: */
     3237        this->inputs->ChangeEnum(VxEnum,VxPicardEnum);
     3238        this->inputs->ChangeEnum(VyEnum,VyPicardEnum);
     3239        this->inputs->ChangeEnum(PressureEnum,PressurePicardEnum);
     3240
     3241        /*Add vx and vy as inputs to the tria element: */
     3242        this->inputs->AddInput(new TriaVertexInput(VxEnum,vx));
     3243        this->inputs->AddInput(new TriaVertexInput(VyEnum,vy));
     3244        this->inputs->AddInput(new TriaVertexInput(VelEnum,vel));
     3245        this->inputs->AddInput(new TriaVertexInput(PressureEnum,pressure));
    21143246
    21153247        /*Free ressources:*/
    2116         delete gauss;
    21173248        xfree((void**)&doflist);
    2118 }
    2119 /*}}}*/
    2120 /*FUNCTION Tria::GetStrainRate2d(double* epsilon,double* xyz_list, GaussTria* gauss, Input* vx_input, Input* vy_input){{{1*/
    2121 void Tria::GetStrainRate2d(double* epsilon,double* xyz_list, GaussTria* gauss, Input* vx_input, Input* vy_input){
    2122         /*Compute the 2d Strain Rate (3 components):
    2123          * epsilon=[exx eyy exy] */
    2124 
    2125         int i;
    2126         double epsilonvx[3];
    2127         double epsilonvy[3];
    2128 
    2129         /*Check that both inputs have been found*/
    2130         if (!vx_input || !vy_input){
    2131                 _error_("Input missing. Here are the input pointers we have for vx: %p, vy: %p\n",vx_input,vy_input);
    2132         }
    2133 
    2134         /*Get strain rate assuming that epsilon has been allocated*/
    2135         vx_input->GetVxStrainRate2d(epsilonvx,xyz_list,gauss);
    2136         vy_input->GetVyStrainRate2d(epsilonvy,xyz_list,gauss);
    2137 
    2138         /*Sum all contributions*/
    2139         for(i=0;i<3;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i];
    2140 }
    2141 /*}}}*/
    2142 /*FUNCTION Tria::GetVectorFromInputs{{{1*/
    2143 void  Tria::GetVectorFromInputs(Vec vector,int input_enum){
    2144 
    2145         int doflist1[NUMVERTICES];
    2146 
    2147         /*Get out if this is not an element input*/
    2148         if (!IsInput(input_enum)) return;
    2149 
    2150         /*Prepare index list*/
    2151         this->GetDofList1(&doflist1[0]);
    2152 
    2153         /*Get input (either in element or material)*/
    2154         Input* input=inputs->GetInput(input_enum);
    2155         if(!input) _error_("Input %s not found in element",EnumToStringx(input_enum));
    2156 
    2157         /*We found the enum.  Use its values to fill into the vector, using the vertices ids: */
    2158         input->GetVectorFromInputs(vector,&doflist1[0]);
    2159 }
    2160 /*}}}*/
    2161 /*FUNCTION Tria::Id {{{1*/
    2162 int    Tria::Id(){
     3249
     3250}
     3251/*}}}*/
     3252/*FUNCTION Tria::InputUpdateFromSolutionDiagnosticHutter {{{1*/
     3253void  Tria::InputUpdateFromSolutionDiagnosticHutter(double* solution){
    21633254       
    2164         return id;
    2165 
    2166 }
    2167 /*}}}*/
    2168 /*FUNCTION Tria::Sid {{{1*/
    2169 int    Tria::Sid(){
     3255        const int numdof=NDOF2*NUMVERTICES;
    21703256       
    2171         return sid;
    2172 
    2173 }
    2174 /*}}}*/
    2175 /*FUNCTION Tria::InputArtificialNoise{{{1*/
    2176 void  Tria::InputArtificialNoise(int enum_type,double min,double max){
    2177 
    2178         Input* input=NULL;
    2179 
    2180         /*Make a copy of the original input: */
    2181         input=(Input*)this->inputs->GetInput(enum_type);
    2182         if(!input)_error_(" could not find old input with enum: %s",EnumToStringx(enum_type));
    2183 
    2184         /*ArtificialNoise: */
    2185         input->ArtificialNoise(min,max);
    2186 }
    2187 /*}}}*/
     3257        int       i;
     3258        int       dummy;
     3259        int*      doflist=NULL;
     3260        double    rho_ice,g;
     3261        double    values[numdof];
     3262        double    vx[NUMVERTICES];
     3263        double    vy[NUMVERTICES];
     3264        double    vz[NUMVERTICES];
     3265        double    vel[NUMVERTICES];
     3266        double    pressure[NUMVERTICES];
     3267        double    thickness[NUMVERTICES];
     3268        double*   vz_ptr=NULL;
     3269        Input*    vz_input=NULL;
     3270       
     3271        /*Get dof list: */
     3272        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     3273
     3274        /*Use the dof list to index into the solution vector: */
     3275        for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
     3276
     3277        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
     3278        for(i=0;i<NUMVERTICES;i++){
     3279                vx[i]=values[i*NDOF2+0];
     3280                vy[i]=values[i*NDOF2+1];
     3281
     3282                /*Check solution*/
     3283                if(isnan(vx[i])) _error_("NaN found in solution vector");
     3284                if(isnan(vy[i])) _error_("NaN found in solution vector");
     3285        }
     3286
     3287        /*Get Vz*/
     3288        vz_input=inputs->GetInput(VzEnum);
     3289        if (vz_input){
     3290                if (vz_input->Enum()!=TriaVertexInputEnum){
     3291                        _error_("Cannot compute Vel as Vz is of type %s",EnumToStringx(vz_input->Enum()));
     3292                }
     3293                vz_input->GetValuesPtr(&vz_ptr,&dummy);
     3294                for(i=0;i<NUMVERTICES;i++) vz[i]=vz_ptr[i];
     3295        }
     3296        else{
     3297                for(i=0;i<NUMVERTICES;i++) vz[i]=0.0;
     3298        }
     3299
     3300        /*Now Compute vel*/
     3301        for(i=0;i<NUMVERTICES;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
     3302
     3303        /*For pressure: we have not computed pressure in this analysis, for this element. We are in 2D,
     3304         *so the pressure is just the pressure at the bedrock: */
     3305        rho_ice=matpar->GetRhoIce();
     3306        g=matpar->GetG();
     3307        GetParameterListOnVertices(&thickness[0],ThicknessEnum);
     3308        for(i=0;i<NUMVERTICES;i++) pressure[i]=rho_ice*g*thickness[i];
     3309
     3310        /*Now, we have to move the previous Vx and Vy inputs  to old
     3311         * status, otherwise, we'll wipe them off: */
     3312        this->inputs->ChangeEnum(VxEnum,VxPicardEnum);
     3313        this->inputs->ChangeEnum(VyEnum,VyPicardEnum);
     3314        this->inputs->ChangeEnum(PressureEnum,PressurePicardEnum);
     3315
     3316        /*Add vx and vy as inputs to the tria element: */
     3317        this->inputs->AddInput(new TriaVertexInput(VxEnum,vx));
     3318        this->inputs->AddInput(new TriaVertexInput(VyEnum,vy));
     3319        this->inputs->AddInput(new TriaVertexInput(VelEnum,vel));
     3320        this->inputs->AddInput(new TriaVertexInput(PressureEnum,pressure));
     3321
     3322        /*Free ressources:*/
     3323        xfree((void**)&doflist);
     3324}
     3325/*}}}*/
     3326
     3327#endif
    21883328
    21893329#ifdef _HAVE_CONTROL_
    … …  
    35864726}
    35874727/*}}}*/
    3588 #endif
    3589 
    3590 /*FUNCTION Tria::InputConvergence{{{1*/
    3591 bool Tria::InputConvergence(double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){
    3592 
    3593         bool    converged=true;
    3594         int     i;
    3595         Input** new_inputs=NULL;
    3596         Input** old_inputs=NULL;
    3597 
    3598         new_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the new inputs
    3599         old_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the old inputs
    3600 
    3601         for(i=0;i<num_enums/2;i++){
    3602                 new_inputs[i]=(Input*)this->inputs->GetInput(enums[2*i+0]);
    3603                 old_inputs[i]=(Input*)this->inputs->GetInput(enums[2*i+1]);
    3604                 if(!new_inputs[i])_error_("%s%s"," could not find input with enum ",EnumToStringx(enums[2*i+0]));
    3605                 if(!old_inputs[i])_error_("%s%s"," could not find input with enum ",EnumToStringx(enums[2*i+0]));
    3606         }
    3607 
    3608         /*ok, we've got the inputs (new and old), now loop throught the number of criterions and fill the eps array:*/
    3609         for(i=0;i<num_criterionenums;i++){
    3610                 IsInputConverged(eps+i,new_inputs,old_inputs,num_enums/2,criterionenums[i]);
    3611                 if(eps[i]>criterionvalues[i]) converged=false;
    3612         }
    3613 
    3614         /*clean up and return*/
    3615         xfree((void**)&new_inputs);
    3616         xfree((void**)&old_inputs);
    3617         return converged;
    3618 }
    3619 /*}}}*/
    3620 /*FUNCTION Tria::InputDepthAverageAtBase {{{1*/
    3621 void  Tria::InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum){
    3622 
    3623         /*New input*/
    3624         Input* oldinput=NULL;
    3625         Input* newinput=NULL;
    3626 
    3627         /*copy input of enum_type*/
    3628         if (object_enum==MeshElementsEnum)
    3629          oldinput=(Input*)this->inputs->GetInput(enum_type);
    3630         else if (object_enum==MaterialsEnum)
    3631          oldinput=(Input*)this->matice->inputs->GetInput(enum_type);
    3632         else
    3633          _error_("object %s not supported yet",EnumToStringx(object_enum));
    3634         if(!oldinput)_error_("%s%s"," could not find old input with enum: ",EnumToStringx(enum_type));
    3635         newinput=(Input*)oldinput->copy();
    3636 
    3637         /*Assign new name (average)*/
    3638         newinput->ChangeEnum(average_enum_type);
    3639 
    3640         /*Add new input to current element*/
    3641         if (object_enum==MeshElementsEnum)
    3642          this->inputs->AddInput((Input*)newinput);
    3643         else if (object_enum==MaterialsEnum)
    3644          this->matice->inputs->AddInput((Input*)newinput);
    3645         else
    3646          _error_("object %s not supported yet",EnumToStringx(object_enum));
    3647 }
    3648 /*}}}*/
    3649 /*FUNCTION Tria::InputDuplicate{{{1*/
    3650 void  Tria::InputDuplicate(int original_enum,int new_enum){
    3651 
    3652         /*Call inputs method*/
    3653         if (IsInput(original_enum)) inputs->DuplicateInput(original_enum,new_enum);
    3654 
    3655 }
    3656 /*}}}*/
    3657 /*FUNCTION Tria::InputScale{{{1*/
    3658 void  Tria::InputScale(int enum_type,double scale_factor){
    3659 
    3660         Input* input=NULL;
    3661 
    3662         /*Make a copy of the original input: */
    3663         input=(Input*)this->inputs->GetInput(enum_type);
    3664         if(!input)_error_(" could not find old input with enum: %s",EnumToStringx(enum_type));
    3665 
    3666         /*Scale: */
    3667         input->Scale(scale_factor);
    3668 }
    3669 /*}}}*/
    3670 /*FUNCTION Tria::InputToResult{{{1*/
    3671 void  Tria::InputToResult(int enum_type,int step,double time){
    3672 
    3673         int    i;
    3674         Input *input = NULL;
    3675 
    3676         /*Go through all the input objects, and find the one corresponding to enum_type, if it exists: */
    3677         if (enum_type==MaterialsRheologyBbarEnum) input=this->matice->inputs->GetInput(enum_type);
    3678         else input=this->inputs->GetInput(enum_type);
    3679         if (!input) _error_("Input %s not found in tria->inputs",EnumToStringx(enum_type));
    3680 
    3681         /*If we don't find it, no big deal, just don't do the transfer. Otherwise, build a new Result
    3682          * object out of the input, with the additional step and time information: */
    3683         this->results->AddObject((Object*)input->SpawnResult(step,time));
    3684         #ifdef _HAVE_CONTROL_
    3685         if(input->Enum()==ControlInputEnum) this->results->AddObject((Object*)((ControlInput*)input)->SpawnGradient(step,time));
    3686         #endif
    3687 }
    3688 /*}}}*/
    3689 /*FUNCTION Tria::InputUpdateFromConstant(int value, int name);{{{1*/
    3690 void  Tria::InputUpdateFromConstant(int constant, int name){
    3691         /*Check that name is an element input*/
    3692         if (!IsInput(name)) return;
    3693 
    3694         /*update input*/
    3695         this->inputs->AddInput(new IntInput(name,constant));
    3696 }
    3697 /*}}}*/
    3698 /*FUNCTION Tria::InputUpdateFromConstant(double value, int name);{{{1*/
    3699 void  Tria::InputUpdateFromConstant(double constant, int name){
    3700         /*Check that name is an element input*/
    3701         if (!IsInput(name)) return;
    3702 
    3703         /*update input*/
    3704         this->inputs->AddInput(new DoubleInput(name,constant));
    3705 }
    3706 /*}}}*/
    3707 /*FUNCTION Tria::InputUpdateFromConstant(bool value, int name);{{{1*/
    3708 void  Tria::InputUpdateFromConstant(bool constant, int name){
    3709         /*Check that name is an element input*/
    3710         if (!IsInput(name)) return;
    3711 
    3712         /*update input*/
    3713         this->inputs->AddInput(new BoolInput(name,constant));
    3714 }
    3715 /*}}}*/
    3716 /*FUNCTION Tria::InputUpdateFromIoModel{{{1*/
    3717 void Tria::InputUpdateFromIoModel(int index, IoModel* iomodel){ //i is the element index
    3718 
    3719         /*Intermediaries*/
    3720         int    i,j;
    3721         int    tria_vertex_ids[3];
    3722         double nodeinputs[3];
    3723         double cmmininputs[3];
    3724         double cmmaxinputs[3];
    3725         bool   control_analysis=false;
    3726         int    num_control_type;
    3727         double yts;
    3728         int    num_cm_responses;
    3729    
    3730         /*Get parameters: */
    3731         iomodel->Constant(&control_analysis,InversionIscontrolEnum);
    3732         iomodel->Constant(&num_control_type,InversionNumControlParametersEnum);
    3733         iomodel->Constant(&yts,ConstantsYtsEnum);
    3734         iomodel->Constant(&num_cm_responses,InversionNumCostFunctionsEnum);
    3735 
    3736         /*Recover vertices ids needed to initialize inputs*/
    3737         for(i=0;i<3;i++){
    3738                 tria_vertex_ids[i]=(int)iomodel->Data(MeshElementsEnum)[3*index+i]; //ids for vertices are in the elements array from Matlab
    3739         }
    3740 
    3741         /*Control Inputs*/
    3742         #ifdef _HAVE_CONTROL_
    3743         if (control_analysis && iomodel->Data(InversionControlParametersEnum)){
    3744                 for(i=0;i<num_control_type;i++){
    3745                         switch((int)iomodel->Data(InversionControlParametersEnum)[i]){
    3746                                 case BalancethicknessThickeningRateEnum:
    3747                                         if (iomodel->Data(BalancethicknessThickeningRateEnum)){
    3748                                                 for(j=0;j<3;j++)nodeinputs[j]=iomodel->Data(BalancethicknessThickeningRateEnum)[tria_vertex_ids[j]-1]/yts;
    3749                                                 for(j=0;j<3;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
    3750                                                 for(j=0;j<3;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
    3751                                                 this->inputs->AddInput(new ControlInput(BalancethicknessThickeningRateEnum,TriaVertexInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
    3752                                         }
    3753                                         break;
    3754                                 case VxEnum:
    3755                                         if (iomodel->Data(VxEnum)){
    3756                                                 for(j=0;j<3;j++)nodeinputs[j]=iomodel->Data(VxEnum)[tria_vertex_ids[j]-1]/yts;
    3757                                                 for(j=0;j<3;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
    3758                                                 for(j=0;j<3;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
    3759                                                 this->inputs->AddInput(new ControlInput(VxEnum,TriaVertexInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
    3760                                         }
    3761                                         break;
    3762                                 case VyEnum:
    3763                                         if (iomodel->Data(VyEnum)){
    3764                                                 for(j=0;j<3;j++)nodeinputs[j]=iomodel->Data(VyEnum)[tria_vertex_ids[j]-1]/yts;
    3765                                                 for(j=0;j<3;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
    3766                                                 for(j=0;j<3;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
    3767                                                 this->inputs->AddInput(new ControlInput(VyEnum,TriaVertexInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
    3768                                         }
    3769                                         break;
    3770                                 case FrictionCoefficientEnum:
    3771                                         if (iomodel->Data(FrictionCoefficientEnum)){
    3772                                                 for(j=0;j<3;j++)nodeinputs[j]=iomodel->Data(FrictionCoefficientEnum)[tria_vertex_ids[j]-1];
    3773                                                 for(j=0;j<3;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i];
    3774                                                 for(j=0;j<3;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i];
    3775                                                 this->inputs->AddInput(new ControlInput(FrictionCoefficientEnum,TriaVertexInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
    3776                                         }
    3777                                         break;
    3778                                 case MaterialsRheologyBbarEnum:
    3779                                         /*Matice will take care of it*/ break;
    3780                                 default:
    3781                                         _error_("Control %s not implemented yet",EnumToStringx((int)iomodel->Data(InversionControlParametersEnum)[i]));
    3782                         }
    3783                 }
    3784         }
    3785         #endif
    3786 
    3787         /*DatasetInputs*/
    3788         if (control_analysis && iomodel->Data(InversionCostFunctionsCoefficientsEnum)){
    3789 
    3790                 /*Create inputs and add to DataSetInput*/
    3791                 DatasetInput* datasetinput=new DatasetInput(InversionCostFunctionsCoefficientsEnum);
    3792                 for(i=0;i<num_cm_responses;i++){
    3793                         for(j=0;j<3;j++)nodeinputs[j]=iomodel->Data(InversionCostFunctionsCoefficientsEnum)[(tria_vertex_ids[j]-1)*num_cm_responses+i];
    3794                         datasetinput->inputs->AddObject(new TriaVertexInput(InversionCostFunctionsCoefficientsEnum,nodeinputs));
    3795                 }
    3796 
    3797                 /*Add datasetinput to element inputs*/
    3798                 this->inputs->AddInput(datasetinput);
    3799         }
    3800 }
    3801 /*}}}*/
    3802 /*FUNCTION Tria::InputUpdateFromSolution {{{1*/
    3803 void  Tria::InputUpdateFromSolution(double* solution){
    3804 
    3805         /*retrive parameters: */
    3806         int analysis_type;
    3807         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    3808 
    3809         /*Just branch to the correct InputUpdateFromSolution generator, according to the type of analysis we are carrying out: */
    3810         switch(analysis_type){
    3811                 case DiagnosticHorizAnalysisEnum:
    3812                         InputUpdateFromSolutionDiagnosticHoriz( solution);
     4728/*FUNCTION Tria::CreateKMatrixAdjointBalancethickness {{{1*/
     4729ElementMatrix* Tria::CreateKMatrixAdjointBalancethickness(void){
     4730
     4731        ElementMatrix* Ke=NULL;
     4732
     4733        /*Get Element Matrix of the forward model*/
     4734        switch(GetElementType()){
     4735                case P1Enum:
     4736                        Ke=CreateKMatrixBalancethickness_CG();
    38134737                        break;
    3814                 case DiagnosticHutterAnalysisEnum:
    3815                         InputUpdateFromSolutionDiagnosticHoriz( solution);
    3816                         break;
    3817                 #ifdef _HAVE_CONTROL_
    3818                 case AdjointHorizAnalysisEnum:
    3819                         InputUpdateFromSolutionAdjointHoriz( solution);
    3820                         break;
    3821                 case AdjointBalancethicknessAnalysisEnum:
    3822                         InputUpdateFromSolutionAdjointBalancethickness( solution);
    3823                         break;
    3824                 #endif
    3825                 case BedSlopeXAnalysisEnum:
    3826                         InputUpdateFromSolutionOneDof(solution,BedSlopeXEnum);
    3827                         break;
    3828                 case BedSlopeYAnalysisEnum:
    3829                         InputUpdateFromSolutionOneDof(solution,BedSlopeYEnum);
    3830                         break;
    3831                 case SurfaceSlopeXAnalysisEnum:
    3832                         InputUpdateFromSolutionOneDof(solution,SurfaceSlopeXEnum);
    3833                         break;
    3834                 case SurfaceSlopeYAnalysisEnum:
    3835                         InputUpdateFromSolutionOneDof(solution,SurfaceSlopeYEnum);
    3836                         break;
    3837                 case PrognosticAnalysisEnum:
    3838                         InputUpdateFromSolutionPrognostic(solution);
    3839                         break;
    3840                 case HydrologyAnalysisEnum:
    3841                         InputUpdateFromSolutionHydrology(solution);
    3842                         break;
    3843                 case BalancethicknessAnalysisEnum:
    3844                         InputUpdateFromSolutionOneDof(solution,ThicknessEnum);
     4738                case P1DGEnum:
     4739                        Ke=CreateKMatrixBalancethickness_DG();
    38454740                        break;
    38464741                default:
    3847                         _error_("analysis %i (%s) not supported yet",analysis_type,EnumToStringx(analysis_type));
    3848         }
    3849 }
    3850 /*}}}*/
    3851 /*FUNCTION Tria::InputUpdateFromSolutionDiagnosticHoriz {{{1*/
    3852 void  Tria::InputUpdateFromSolutionDiagnosticHoriz(double* solution){
    3853        
     4742                        _error_("Element type %s not supported yet",EnumToStringx(GetElementType()));
     4743        }
     4744
     4745        /*Transpose and return Ke*/
     4746        Ke->Transpose();
     4747        return Ke;
     4748}
     4749/*}}}*/
     4750/*FUNCTION Tria::InputUpdateFromSolutionAdjointHoriz {{{1*/
     4751void  Tria::InputUpdateFromSolutionAdjointHoriz(double* solution){
     4752
    38544753        const int numdof=NDOF2*NUMVERTICES;
    38554754
    38564755        int       i;
    3857         int       dummy;
    38584756        int*      doflist=NULL;
    3859         double    rho_ice,g;
    38604757        double    values[numdof];
    3861         double    vx[NUMVERTICES];
    3862         double    vy[NUMVERTICES];
    3863         double    vz[NUMVERTICES];
    3864         double    vel[NUMVERTICES];
    3865         double    pressure[NUMVERTICES];
    3866         double    thickness[NUMVERTICES];
    3867        
     4758        double    lambdax[NUMVERTICES];
     4759        double    lambday[NUMVERTICES];
     4760
    38684761        /*Get dof list: */
    38694762        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    … …  
    38744767        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    38754768        for(i=0;i<NUMVERTICES;i++){
    3876                 vx[i]=values[i*NDOF2+0];
    3877                 vy[i]=values[i*NDOF2+1];
     4769                lambdax[i]=values[i*NDOF2+0];
     4770                lambday[i]=values[i*NDOF2+1];
    38784771
    38794772                /*Check solution*/
    3880                 if(isnan(vx[i])) _error_("NaN found in solution vector");
    3881                 if(isnan(vy[i])) _error_("NaN found in solution vector");
    3882         }
    3883 
    3884         /*Get Vz and compute vel*/
    3885         GetParameterListOnVertices(&vz[0],VzEnum,0);
    3886         for(i=0;i<NUMVERTICES;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
    3887 
    3888         /*For pressure: we have not computed pressure in this analysis, for this element. We are in 2D,
    3889          *so the pressure is just the pressure at the bedrock: */
    3890         rho_ice=matpar->GetRhoIce();
    3891         g=matpar->GetG();
    3892         GetParameterListOnVertices(&thickness[0],ThicknessEnum);
    3893         for(i=0;i<NUMVERTICES;i++) pressure[i]=rho_ice*g*thickness[i];
    3894 
    3895         /*Now, we have to move the previous Vx and Vy inputs  to old
    3896          * status, otherwise, we'll wipe them off: */
    3897         this->inputs->ChangeEnum(VxEnum,VxPicardEnum);
    3898         this->inputs->ChangeEnum(VyEnum,VyPicardEnum);
    3899         this->inputs->ChangeEnum(PressureEnum,PressurePicardEnum);
     4773                if(isnan(lambdax[i])) _error_("NaN found in solution vector");
     4774                if(isnan(lambday[i])) _error_("NaN found in solution vector");
     4775        }
    39004776
    39014777        /*Add vx and vy as inputs to the tria element: */
    3902         this->inputs->AddInput(new TriaVertexInput(VxEnum,vx));
    3903         this->inputs->AddInput(new TriaVertexInput(VyEnum,vy));
    3904         this->inputs->AddInput(new TriaVertexInput(VelEnum,vel));
    3905         this->inputs->AddInput(new TriaVertexInput(PressureEnum,pressure));
     4778        this->inputs->AddInput(new TriaVertexInput(AdjointxEnum,lambdax));
     4779        this->inputs->AddInput(new TriaVertexInput(AdjointyEnum,lambday));
    39064780
    39074781        /*Free ressources:*/
    39084782        xfree((void**)&doflist);
    3909 
    3910 }
    3911 /*}}}*/
    3912 /*FUNCTION Tria::InputUpdateFromSolutionDiagnosticHutter {{{1*/
    3913 void  Tria::InputUpdateFromSolutionDiagnosticHutter(double* solution){
    3914        
    3915         const int numdof=NDOF2*NUMVERTICES;
    3916        
     4783}
     4784/*}}}*/
     4785/*FUNCTION Tria::InputUpdateFromSolutionAdjointBalancethickness {{{1*/
     4786void  Tria::InputUpdateFromSolutionAdjointBalancethickness(double* solution){
     4787
     4788        const int numdof=NDOF1*NUMVERTICES;
     4789
    39174790        int       i;
    3918         int       dummy;
    39194791        int*      doflist=NULL;
    3920         double    rho_ice,g;
    39214792        double    values[numdof];
    3922         double    vx[NUMVERTICES];
    3923         double    vy[NUMVERTICES];
    3924         double    vz[NUMVERTICES];
    3925         double    vel[NUMVERTICES];
    3926         double    pressure[NUMVERTICES];
    3927         double    thickness[NUMVERTICES];
    3928         double*   vz_ptr=NULL;
    3929         Input*    vz_input=NULL;
    3930        
     4793        double    lambda[NUMVERTICES];
     4794
    39314795        /*Get dof list: */
    39324796        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    … …  
    39364800
    39374801        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    3938         for(i=0;i<NUMVERTICES;i++){
    3939                 vx[i]=values[i*NDOF2+0];
    3940                 vy[i]=values[i*NDOF2+1];
    3941 
    3942                 /*Check solution*/
    3943                 if(isnan(vx[i])) _error_("NaN found in solution vector");
    3944                 if(isnan(vy[i])) _error_("NaN found in solution vector");
    3945         }
    3946 
    3947         /*Get Vz*/
    3948         vz_input=inputs->GetInput(VzEnum);
    3949         if (vz_input){
    3950                 if (vz_input->Enum()!=TriaVertexInputEnum){
    3951                         _error_("Cannot compute Vel as Vz is of type %s",EnumToStringx(vz_input->Enum()));
    3952                 }
    3953                 vz_input->GetValuesPtr(&vz_ptr,&dummy);
    3954                 for(i=0;i<NUMVERTICES;i++) vz[i]=vz_ptr[i];
    3955         }
    3956         else{
    3957                 for(i=0;i<NUMVERTICES;i++) vz[i]=0.0;
    3958         }
    3959 
    3960         /*Now Compute vel*/
    3961         for(i=0;i<NUMVERTICES;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
    3962 
    3963         /*For pressure: we have not computed pressure in this analysis, for this element. We are in 2D,
    3964          *so the pressure is just the pressure at the bedrock: */
    3965         rho_ice=matpar->GetRhoIce();
    3966         g=matpar->GetG();
    3967         GetParameterListOnVertices(&thickness[0],ThicknessEnum);
    3968         for(i=0;i<NUMVERTICES;i++) pressure[i]=rho_ice*g*thickness[i];
    3969 
    3970         /*Now, we have to move the previous Vx and Vy inputs  to old
    3971          * status, otherwise, we'll wipe them off: */
    3972         this->inputs->ChangeEnum(VxEnum,VxPicardEnum);
    3973         this->inputs->ChangeEnum(VyEnum,VyPicardEnum);
    3974         this->inputs->ChangeEnum(PressureEnum,PressurePicardEnum);
     4802        for(i=0;i<numdof;i++){
     4803                lambda[i]=values[i];
     4804                if(isnan(lambda[i])) _error_("NaN found in solution vector");
     4805        }
    39754806
    39764807        /*Add vx and vy as inputs to the tria element: */
    3977         this->inputs->AddInput(new TriaVertexInput(VxEnum,vx));
    3978         this->inputs->AddInput(new TriaVertexInput(VyEnum,vy));
    3979         this->inputs->AddInput(new TriaVertexInput(VelEnum,vel));
    3980         this->inputs->AddInput(new TriaVertexInput(PressureEnum,pressure));
     4808        this->inputs->AddInput(new TriaVertexInput(AdjointEnum,lambda));
    39814809
    39824810        /*Free ressources:*/
    … …  
    39844812}
    39854813/*}}}*/
    3986 /*FUNCTION Tria::InputUpdateFromSolutionOneDof{{{1*/
    3987 void  Tria::InputUpdateFromSolutionOneDof(double* solution,int enum_type){
    3988 
    3989         const int numdof          = NDOF1*NUMVERTICES;
    3990 
    3991         int*      doflist=NULL;
    3992         double    values[numdof];
     4814#endif
     4815
     4816#ifdef _HAVE_HYDROLOGY_
     4817/*FUNCTION Tria::CreateHydrologyWaterVelocityInput {{{1*/
     4818void Tria::CreateHydrologyWaterVelocityInput(void){
     4819
     4820        /*material parameters: */
     4821        double mu_water=0.001787;  //unit= [N s/m2]
     4822        double w;
     4823        double rho_ice, rho_water, g;
     4824        double dsdx,dsdy,dbdx,dbdy;
     4825        double vx[NUMVERTICES];
     4826        double vy[NUMVERTICES];
     4827        GaussTria *gauss = NULL;
     4828
     4829        /*Retrieve all inputs and parameters*/
     4830        rho_ice=matpar->GetRhoIce();
     4831        rho_water=matpar->GetRhoWater();
     4832        g=matpar->GetG();
     4833        Input* surfaceslopex_input=inputs->GetInput(SurfaceSlopeXEnum); _assert_(surfaceslopex_input);
     4834        Input* surfaceslopey_input=inputs->GetInput(SurfaceSlopeYEnum); _assert_(surfaceslopey_input);
     4835        Input* bedslopex_input=inputs->GetInput(BedSlopeXEnum);         _assert_(bedslopex_input);
     4836        Input* bedslopey_input=inputs->GetInput(BedSlopeYEnum);         _assert_(bedslopey_input);
     4837        Input* watercolumn_input=inputs->GetInput(WatercolumnEnum);     _assert_(watercolumn_input);
     4838
     4839        gauss=new GaussTria();
     4840        for (int iv=0;iv<NUMVERTICES;iv++){
     4841                gauss->GaussVertex(iv);
     4842                surfaceslopex_input->GetParameterValue(&dsdx,gauss);
     4843                surfaceslopey_input->GetParameterValue(&dsdy,gauss);
     4844                bedslopex_input->GetParameterValue(&dbdx,gauss);
     4845                bedslopey_input->GetParameterValue(&dbdy,gauss);
     4846                watercolumn_input->GetParameterValue(&w,gauss);
     4847
     4848                /* Water velocity x and y components */
     4849                vx[iv]= - pow(w,2)/(12 * mu_water)*(rho_ice*g*dsdx+(rho_water-rho_ice)*g*dbdx);
     4850                vy[iv]= - pow(w,2)/(12 * mu_water)*(rho_ice*g*dsdy+(rho_water-rho_ice)*g*dbdy);
     4851                //vx[iv]=0.001;
     4852                //vy[iv]=0;
     4853        }
     4854
     4855        /*clean-up*/
     4856        delete gauss;
     4857
     4858        /*Add to inputs*/
     4859        this->inputs->AddInput(new TriaVertexInput(HydrologyWaterVxEnum,vx));
     4860        this->inputs->AddInput(new TriaVertexInput(HydrologyWaterVyEnum,vy));
     4861}
     4862/*}}}*/
     4863/*FUNCTION Tria::CreateKMatrixHydrology{{{1*/
     4864ElementMatrix* Tria::CreateKMatrixHydrology(void){
     4865
     4866        /*Constants*/
     4867        const int    numdof=NDOF1*NUMVERTICES;
     4868
     4869        /*Intermediaries */
     4870        int        artdiff;
     4871        int        i,j,ig;
     4872        double     Jdettria,DL_scalar,dt,h;
     4873        double     vx,vy,vel,dvxdx,dvydy;
     4874        double     dvx[2],dvy[2];
     4875        double     v_gauss[2]={0.0};
     4876        double     xyz_list[NUMVERTICES][3];
     4877        double     L[NUMVERTICES];
     4878        double     B[2][NUMVERTICES];
     4879        double     Bprime[2][NUMVERTICES];
     4880        double     K[2][2]                        ={0.0};
     4881        double     KDL[2][2]                      ={0.0};
     4882        double     DL[2][2]                        ={0.0};
     4883        double     DLprime[2][2]                   ={0.0};
     4884        GaussTria *gauss=NULL;
     4885
     4886        /*Initialize Element matrix*/
     4887        ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
     4888
     4889        /*Create water velocity vx and vy from current inputs*/
     4890        CreateHydrologyWaterVelocityInput();
     4891
     4892        /*Retrieve all inputs and parameters*/
     4893        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     4894        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     4895        this->parameters->FindParam(&artdiff,HydrologyStabilizationEnum);
     4896        Input* vx_input=inputs->GetInput(HydrologyWaterVxEnum); _assert_(vx_input);
     4897        Input* vy_input=inputs->GetInput(HydrologyWaterVyEnum); _assert_(vy_input);
     4898        h=this->GetArea();
     4899
     4900        /* Start  looping on the number of gaussian points: */
     4901        gauss=new GaussTria(2);
     4902        for (ig=gauss->begin();ig<gauss->end();ig++){
     4903
     4904                gauss->GaussPoint(ig);
     4905
     4906                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
     4907                GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
     4908
     4909                vx_input->GetParameterValue(&vx,gauss);
     4910                vy_input->GetParameterValue(&vy,gauss);
     4911                vx_input->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
     4912                vy_input->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
     4913
     4914                DL_scalar=gauss->weight*Jdettria;
     4915
     4916                TripleMultiply( &L[0],1,numdof,1,
     4917                                        &DL_scalar,1,1,0,
     4918                                        &L[0],1,numdof,0,
     4919                                        &Ke->values[0],1);
     4920
     4921                GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss);
     4922                GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
     4923
     4924                dvxdx=dvx[0];
     4925                dvydy=dvy[1];
     4926                DL_scalar=dt*gauss->weight*Jdettria;
     4927
     4928                DL[0][0]=DL_scalar*dvxdx;
     4929                DL[1][1]=DL_scalar*dvydy;
     4930                DLprime[0][0]=DL_scalar*vx;
     4931                DLprime[1][1]=DL_scalar*vy;
     4932
     4933                TripleMultiply( &B[0][0],2,numdof,1,
     4934                                        &DL[0][0],2,2,0,
     4935                                        &B[0][0],2,numdof,0,
     4936                                        &Ke->values[0],1);
     4937
     4938                TripleMultiply( &B[0][0],2,numdof,1,
     4939                                        &DLprime[0][0],2,2,0,
     4940                                        &Bprime[0][0],2,numdof,0,
     4941                                        &Ke->values[0],1);
     4942
     4943                if(artdiff){
     4944                        vel=sqrt(pow(vx,2.)+pow(vy,2.));
     4945                        K[0][0]=h/(2*vel)*vx*vx;
     4946                        K[1][0]=h/(2*vel)*vy*vx;
     4947                        K[0][1]=h/(2*vel)*vx*vy;
     4948                        K[1][1]=h/(2*vel)*vy*vy;
     4949                        KDL[0][0]=DL_scalar*K[0][0];
     4950                        KDL[1][0]=DL_scalar*K[1][0];
     4951                        KDL[0][1]=DL_scalar*K[0][1];
     4952                        KDL[1][1]=DL_scalar*K[1][1];
     4953
     4954                        TripleMultiply( &Bprime[0][0],2,numdof,1,
     4955                                                &KDL[0][0],2,2,0,
     4956                                                &Bprime[0][0],2,numdof,0,
     4957                                                &Ke->values[0],1);
     4958                }
     4959        }
     4960
     4961        /*Clean up and return*/
     4962        delete gauss;
     4963        return Ke;
     4964}
     4965/*}}}*/
     4966/*FUNCTION Tria::CreatePVectorHydrology {{{1*/
     4967ElementVector* Tria::CreatePVectorHydrology(void){
     4968
     4969        /*Constants*/
     4970        const int    numdof=NDOF1*NUMVERTICES;
     4971
     4972        /*Intermediaries */
     4973        int        i,j,ig;
     4974        double     Jdettria,dt;
     4975        double     basal_melting_g;
     4976        double     old_watercolumn_g;
     4977        double     xyz_list[NUMVERTICES][3];
     4978        double     basis[numdof];
     4979        GaussTria* gauss=NULL;
     4980
     4981        /*Initialize Element vector*/
     4982        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
     4983
     4984        /*Retrieve all inputs and parameters*/
     4985        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     4986        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     4987        Input* basal_melting_input=inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(basal_melting_input);
     4988        Input* old_watercolumn_input=inputs->GetInput(WaterColumnOldEnum);         _assert_(old_watercolumn_input);
     4989
     4990        /*Initialize basal_melting_correction_g to 0, do not forget!:*/
     4991        /* Start  looping on the number of gaussian points: */
     4992        gauss=new GaussTria(2);
     4993        for(ig=gauss->begin();ig<gauss->end();ig++){
     4994
     4995                gauss->GaussPoint(ig);
     4996
     4997                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
     4998                GetNodalFunctions(basis, gauss);
     4999
     5000                basal_melting_input->GetParameterValue(&basal_melting_g,gauss);
     5001                old_watercolumn_input->GetParameterValue(&old_watercolumn_g,gauss);
     5002
     5003                if(dt)for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(old_watercolumn_g+dt*basal_melting_g)*basis[i];
     5004                else  for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*basal_melting_g*basis[i];
     5005        }
     5006               
     5007        /*Clean up and return*/
     5008        delete gauss;
     5009        return pe;
     5010}
     5011/*}}}*/
     5012/*FUNCTION Tria::GetSolutionFromInputsHydrology{{{1*/
     5013void  Tria::GetSolutionFromInputsHydrology(Vec solution){
     5014
     5015        const int    numdof=NDOF1*NUMVERTICES;
     5016
     5017        int i,dummy;
     5018        int*         doflist=NULL;
     5019        double       watercolumn;
     5020        double       values[numdof];
     5021        GaussTria*   gauss=NULL;
    39935022
    39945023        /*Get dof list: */
    39955024        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    39965025
    3997         /*Use the dof list to index into the solution vector: */
    3998         for(int i=0;i<numdof;i++){
    3999                 values[i]=solution[doflist[i]];
    4000                 if(isnan(values[i])) _error_("NaN found in solution vector");
    4001         }
    4002 
    4003         /*Add input to the element: */
    4004         this->inputs->AddInput(new TriaVertexInput(enum_type,values));
     5026        /*Get inputs*/
     5027        Input* watercolumn_input=inputs->GetInput(WatercolumnEnum); _assert_(watercolumn_input);
     5028
     5029        /*Ok, we have watercolumn values, fill in watercolumn array: */
     5030        /*P1 element only for now*/
     5031        gauss=new GaussTria();
     5032        for(i=0;i<NUMVERTICES;i++){
     5033
     5034                gauss->GaussVertex(i);
     5035
     5036                /*Recover watercolumn*/
     5037                watercolumn_input->GetParameterValue(&watercolumn,gauss);
     5038                values[i]=watercolumn;
     5039        }
     5040
     5041        VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES);
    40055042
    40065043        /*Free ressources:*/
    4007         xfree((void**)&doflist);
    4008 }
    4009 /*}}}*/
    4010 /*FUNCTION Tria::InputUpdateFromSolutionPrognostic{{{1*/
    4011 void  Tria::InputUpdateFromSolutionPrognostic(double* solution){
    4012 
    4013         /*Intermediaries*/
    4014         const int numdof = NDOF1*NUMVERTICES;
    4015 
    4016         int       i,dummy,hydroadjustment;
    4017         int*      doflist=NULL;
    4018         double    rho_ice,rho_water;
    4019         double    values[numdof];
    4020         double    bed[numdof];
    4021         double    surface[numdof];
    4022         double    *bed_ptr = NULL;
    4023         double    *thickness_ptr = NULL;
    4024         double    *surface_ptr = NULL;
    4025 
    4026         /*Get dof list: */
    4027         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    4028 
    4029         /*Use the dof list to index into the solution vector: */
    4030         for(i=0;i<numdof;i++){
    4031                 values[i]=solution[doflist[i]];
    4032                 if(isnan(values[i])) _error_("NaN found in solution vector");
    4033                 /*Constrain thickness to be at least 1m*/
    4034                 if(values[i]<1) values[i]=1;
    4035         }
    4036 
    4037         /*Get previous bed, thickness and surface*/
    4038         Input* bed_input=inputs->GetInput(BedEnum);             _assert_(bed_input);
    4039         Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
    4040         Input* surface_input=inputs->GetInput(SurfaceEnum);     _assert_(surface_input);
    4041         bed_input->GetValuesPtr(&bed_ptr,&dummy);
    4042         thickness_input->GetValuesPtr(&thickness_ptr,&dummy);
    4043         surface_input->GetValuesPtr(&surface_ptr,&dummy);
    4044 
    4045         /*Fing PrognosticHydrostaticAdjustment to figure out how to update the geometry:*/
    4046         this->parameters->FindParam(&hydroadjustment,PrognosticHydrostaticAdjustmentEnum);
    4047 
    4048         /*recover material parameters: */
    4049         rho_ice=matpar->GetRhoIce();
    4050         rho_water=matpar->GetRhoWater();
    4051 
    4052         for(i=0;i<numdof;i++) {
    4053                 /*If shelf: hydrostatic equilibrium*/
    4054                 if (this->nodes[i]->IsOnSheet()){
    4055                         surface[i]=bed_ptr[i]+values[i]; //surface=oldbed+newthickness
    4056                         bed[i]=bed_ptr[i]; //do nothing
    4057                 }
    4058                 else{ //this is an ice shelf
    4059 
    4060                         if(hydroadjustment==AbsoluteEnum){
    4061                                 surface[i]=values[i]*(1-rho_ice/rho_water);
    4062                                 bed[i]=values[i]*(-rho_ice/rho_water);
    4063                         }
    4064                         else if(hydroadjustment==IncrementalEnum){
    4065                                 surface[i]=surface_ptr[i]+(1.0-rho_ice/rho_water)*(values[i]-thickness_ptr[i]); //surface = oldsurface + (1-di) * dH
    4066                                 bed[i]=bed_ptr[i]-rho_ice/rho_water*(values[i]-thickness_ptr[i]); //bed = oldbed + di * dH
    4067                         }
    4068                         else _error_("Hydrostatic adjustment %i (%s) not supported yet",hydroadjustment,EnumToStringx(hydroadjustment));
    4069                 }
    4070         }
    4071 
    4072         /*Add input to the element: */
    4073         this->inputs->AddInput(new TriaVertexInput(ThicknessEnum,values));
    4074         this->inputs->AddInput(new TriaVertexInput(SurfaceEnum,surface));
    4075         this->inputs->AddInput(new TriaVertexInput(BedEnum,bed));
    4076 
    4077         /*Free ressources:*/
     5044        delete gauss;
    40785045        xfree((void**)&doflist);
    40795046}
    … …  
    41075074}
    41085075/*}}}*/
    4109 /*FUNCTION Tria::InputUpdateFromVector(double* vector, int name, int type);{{{1*/
    4110 void  Tria::InputUpdateFromVector(double* vector, int name, int type){
    4111 
    4112         /*Check that name is an element input*/
    4113         if (!IsInput(name)) return;
    4114 
    4115         switch(type){
    4116 
    4117                 case VertexEnum:
    4118 
    4119                         /*New TriaVertexInput*/
    4120                         double values[3];
    4121 
    4122                         /*Get values on the 3 vertices*/
    4123                         for (int i=0;i<3;i++){
    4124                                 values[i]=vector[this->nodes[i]->GetVertexDof()];
    4125                         }
    4126 
    4127                         /*update input*/
    4128                         if (name==MaterialsRheologyBbarEnum || name==MaterialsRheologyBEnum){
    4129                                 matice->inputs->AddInput(new TriaVertexInput(name,values));
    4130                         }
    4131                         else{
    4132                                 this->inputs->AddInput(new TriaVertexInput(name,values));
    4133                         }
    4134                         return;
    4135 
    4136                 default:
    4137                         _error_("type %i (%s) not implemented yet",type,EnumToStringx(type));
    4138         }
    4139 }
    4140 /*}}}*/
    4141 /*FUNCTION Tria::InputUpdateFromVector(int* vector, int name, int type);{{{1*/
    4142 void  Tria::InputUpdateFromVector(int* vector, int name, int type){
    4143         _error_(" not supported yet!");
    4144 }
    4145 /*}}}*/
    4146 /*FUNCTION Tria::InputUpdateFromVector(bool* vector, int name, int type);{{{1*/
    4147 void  Tria::InputUpdateFromVector(bool* vector, int name, int type){
    4148         _error_(" not supported yet!");
    4149 }
    4150 /*}}}*/
    4151 
     5076#endif
    41525077
    41535078#ifdef _HAVE_DAKOTA_
    … …  
    42655190/*}}}*/
    42665191#endif
    4267 /*FUNCTION Tria::InputCreate(double scalar,int enum,int code);{{{1*/
    4268 void Tria::InputCreate(double scalar,int name,int code){
    4269 
    4270         /*Check that name is an element input*/
    4271         if (!IsInput(name)) return;
    4272        
    4273         if ((code==5) || (code==1)){ //boolean
    4274                 this->inputs->AddInput(new BoolInput(name,(bool)scalar));
    4275         }
    4276         else if ((code==6) || (code==2)){ //integer
    4277                 this->inputs->AddInput(new IntInput(name,(int)scalar));
    4278         }
    4279         else if ((code==7) || (code==3)){ //double
    4280                 this->inputs->AddInput(new DoubleInput(name,(int)scalar));
    4281         }
    4282         else _error_("%s%i"," could not recognize nature of vector from code ",code);
    4283 
    4284 }
    4285 /*}}}*/
    4286 /*FUNCTION Tria::InputCreate(double* vector,int index,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){{{1*/
    4287 void Tria::InputCreate(double* vector, int index,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){//index into elements
     5192
     5193#ifdef _HAVE_BALANCED_
     5194/*FUNCTION Tria::CreateKMatrixBalancethickness {{{1*/
     5195ElementMatrix* Tria::CreateKMatrixBalancethickness(void){
     5196
     5197        switch(GetElementType()){
     5198                case P1Enum:
     5199                        return CreateKMatrixBalancethickness_CG();
     5200                case P1DGEnum:
     5201                        return CreateKMatrixBalancethickness_DG();
     5202                default:
     5203                        _error_("Element type %s not supported yet",EnumToStringx(GetElementType()));
     5204        }
     5205
     5206}
     5207/*}}}*/
     5208/*FUNCTION Tria::CreateKMatrixBalancethickness_CG {{{1*/
     5209ElementMatrix* Tria::CreateKMatrixBalancethickness_CG(void){
     5210
     5211        /*Constants*/
     5212        const int    numdof=NDOF1*NUMVERTICES;
     5213
     5214        /*Intermediaries */
     5215        int        stabilization;
     5216        int        i,j,ig,dim;
     5217        double     Jdettria,vx,vy,dvxdx,dvydy,vel,h;
     5218        double     dvx[2],dvy[2];
     5219        double     xyz_list[NUMVERTICES][3];
     5220        double     L[NUMVERTICES];
     5221        double     B[2][NUMVERTICES];
     5222        double     Bprime[2][NUMVERTICES];
     5223        double     K[2][2]                          = {0.0};
     5224        double     KDL[2][2]                        = {0.0};
     5225        double     DL[2][2]                         = {0.0};
     5226        double     DLprime[2][2]                    = {0.0};
     5227        double     DL_scalar;
     5228        GaussTria *gauss                            = NULL;
     5229
     5230        /*Initialize Element matrix*/
     5231        ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
     5232
     5233        /*Retrieve all Inputs and parameters: */
     5234        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     5235        this->parameters->FindParam(&stabilization,BalancethicknessStabilizationEnum);
     5236        this->parameters->FindParam(&dim,MeshDimensionEnum);
     5237        Input* vxaverage_input=NULL;
     5238        Input* vyaverage_input=NULL;
     5239        if(dim==2){
     5240                vxaverage_input=inputs->GetInput(VxEnum); _assert_(vxaverage_input);
     5241                vyaverage_input=inputs->GetInput(VyEnum); _assert_(vyaverage_input);
     5242        }
     5243        else{
     5244                vxaverage_input=inputs->GetInput(VxAverageEnum); _assert_(vxaverage_input);
     5245                vyaverage_input=inputs->GetInput(VyAverageEnum); _assert_(vyaverage_input);
     5246        }
     5247        h=sqrt(2*this->GetArea());
     5248
     5249        ///*Create Artificial diffusivity once for all if requested*/
     5250        //if(stabilization){
     5251        //      gauss=new GaussTria();
     5252        //      gauss->GaussCenter();
     5253        //      GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
     5254        //      delete gauss;
     5255
     5256        //      vxaverage_input->GetParameterAverage(&vx);
     5257        //      vyaverage_input->GetParameterAverage(&vy);
     5258        //      K[0][0]=pow(Jdettria,(double).5)/2.0*fabs(vx);
     5259        //      K[1][1]=pow(Jdettria,(double).5)/2.0*fabs(vy);
     5260        //}
     5261
     5262        /*Start looping on the number of gaussian points:*/
     5263        gauss=new GaussTria(2);
     5264        for (ig=gauss->begin();ig<gauss->end();ig++){
     5265
     5266                gauss->GaussPoint(ig);
     5267
     5268                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
     5269                GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss);
     5270                GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
     5271
     5272                vxaverage_input->GetParameterValue(&vx,gauss);
     5273                vyaverage_input->GetParameterValue(&vy,gauss);
     5274                vxaverage_input->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
     5275                vyaverage_input->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
     5276
     5277                dvxdx=dvx[0];
     5278                dvydy=dvy[1];
     5279                DL_scalar=gauss->weight*Jdettria;
     5280
     5281                DL[0][0]=DL_scalar*dvxdx;
     5282                DL[1][1]=DL_scalar*dvydy;
     5283
     5284                DLprime[0][0]=DL_scalar*vx;
     5285                DLprime[1][1]=DL_scalar*vy;
     5286
     5287                TripleMultiply( &B[0][0],2,numdof,1,
     5288                                        &DL[0][0],2,2,0,
     5289                                        &B[0][0],2,numdof,0,
     5290                                        &Ke->values[0],1);
     5291
     5292                TripleMultiply( &B[0][0],2,numdof,1,
     5293                                        &DLprime[0][0],2,2,0,
     5294                                        &Bprime[0][0],2,numdof,0,
     5295                                        &Ke->values[0],1);
     5296
     5297                if(stabilization==1){
     5298                        vel=sqrt(pow(vx,2.)+pow(vy,2.));
     5299                        K[0][0]=h/(2*vel)*vx*vx;
     5300                        K[1][0]=h/(2*vel)*vy*vx;
     5301                        K[0][1]=h/(2*vel)*vx*vy;
     5302                        K[1][1]=h/(2*vel)*vy*vy;
     5303                        KDL[0][0]=DL_scalar*K[0][0];
     5304                        KDL[1][0]=DL_scalar*K[1][0];
     5305                        KDL[0][1]=DL_scalar*K[0][1];
     5306                        KDL[1][1]=DL_scalar*K[1][1];
     5307
     5308                        //KDL[0][0]=DL_scalar*K[0][0];
     5309                        //KDL[1][1]=DL_scalar*K[1][1];
     5310
     5311                        TripleMultiply( &Bprime[0][0],2,numdof,1,
     5312                                                &KDL[0][0],2,2,0,
     5313                                                &Bprime[0][0],2,numdof,0,
     5314                                                &Ke->values[0],1);
     5315                }
     5316        }
     5317
     5318        /*Clean up and return*/
     5319        delete gauss;
     5320        return Ke;
     5321}
     5322/*}}}*/
     5323/*FUNCTION Tria::CreateKMatrixBalancethickness_DG {{{1*/
     5324ElementMatrix* Tria::CreateKMatrixBalancethickness_DG(void){
     5325
     5326        /*Constants*/
     5327        const int  numdof=NDOF1*NUMVERTICES;
    42885328
    42895329        /*Intermediaries*/
    4290         int    i,j,t;
    4291         int    tria_vertex_ids[3];
    4292         int    row;
    4293         double nodeinputs[3];
    4294         double time;
    4295         TransientInput* transientinput=NULL;
    4296         int    numberofvertices;
    4297         int    numberofelements;
    4298         double yts;
    4299 
    4300 
    4301         /*Fetch parameters: */
    4302         iomodel->Constant(&numberofvertices,MeshNumberofverticesEnum);
    4303         iomodel->Constant(&numberofelements,MeshNumberofelementsEnum);
    4304         iomodel->Constant(&yts,ConstantsYtsEnum);
    4305 
    4306         /*Branch on type of vector: nodal or elementary: */
    4307         if(vector_type==1){ //nodal vector
    4308 
    4309                 /*Recover vertices ids needed to initialize inputs*/
    4310                 for(i=0;i<3;i++){
    4311                         tria_vertex_ids[i]=(int)iomodel->Data(MeshElementsEnum)[3*index+i]; //ids for vertices are in the elements array from Matlab
    4312                 }
    4313 
    4314                 /*Are we in transient or static? */
    4315                 if(M==numberofvertices){
    4316 
    4317                         /*create input values: */
    4318                         for(i=0;i<3;i++)nodeinputs[i]=(double)vector[tria_vertex_ids[i]-1];
    4319 
    4320                         /*process units: */
    4321                         UnitConversion(&nodeinputs[0], 3 ,ExtToIuEnum, vector_enum);
    4322 
    4323                         /*create static input: */
    4324                         this->inputs->AddInput(new TriaVertexInput(vector_enum,nodeinputs));
    4325                 }
    4326                 else if(M==numberofvertices+1){
    4327                         /*create transient input: */
    4328                         for(t=0;t<N;t++){ //N is the number of times
    4329 
    4330                                 /*create input values: */
    4331                                 for(i=0;i<3;i++){
    4332                                         row=tria_vertex_ids[i]-1;
    4333                                         nodeinputs[i]=(double)vector[N*row+t];
    4334                                 }
    4335 
    4336                                 /*process units: */
    4337                                 UnitConversion(&nodeinputs[0], 3 ,ExtToIuEnum, vector_enum);
    4338 
    4339                                 /*time? :*/
    4340                                 time=(double)vector[(M-1)*N+t]*yts;
    4341 
    4342                                 if(t==0) transientinput=new TransientInput(vector_enum);
    4343                                 transientinput->AddTimeInput(new TriaVertexInput(vector_enum,nodeinputs),time);
    4344                         }
    4345                         this->inputs->AddInput(transientinput);
    4346                 }
    4347                 else _error_("nodal vector is either numberofnodes (%i), or numberofnodes+1 long. Field provided is %i long",numberofvertices,M);
    4348         }
    4349         else if(vector_type==2){ //element vector
    4350                 /*Are we in transient or static? */
    4351                 if(M==numberofelements){
    4352 
    4353                         /*static mode: create an input out of the element value: */
    4354 
    4355                         if (code==5){ //boolean
    4356                                 this->inputs->AddInput(new BoolInput(vector_enum,(bool)vector[index]));
    4357                         }
    4358                         else if (code==6){ //integer
    4359                                 this->inputs->AddInput(new IntInput(vector_enum,(int)vector[index]));
    4360                         }
    4361                         else if (code==7){ //double
    4362                                 this->inputs->AddInput(new DoubleInput(vector_enum,(double)vector[index]));
    4363                         }
    4364                         else _error_("%s%i"," could not recognize nature of vector from code ",code);
    4365                 }
    4366                 else {
    4367                         _error_("transient elementary inputs not supported yet!");
    4368                 }
    4369         }
    4370         else{
    4371                 _error_("Cannot add input for vector type %i (not supported)",vector_type);
    4372         }
    4373 
    4374 }
    4375 /*}}}*/
    4376 /*FUNCTION Tria::IsInput{{{1*/
    4377 bool Tria::IsInput(int name){
    4378         if (
    4379                                 name==ThicknessEnum ||
    4380                                 name==SurfaceEnum ||
    4381                                 name==BedEnum ||
    4382                                 name==SurfaceSlopeXEnum ||
    4383                                 name==SurfaceSlopeYEnum ||
    4384                                 name==BasalforcingsMeltingRateEnum ||
    4385                                 name==WatercolumnEnum ||
    4386                                 name==SurfaceforcingsMassBalanceEnum ||
    4387                                 name==SurfaceAreaEnum||
    4388                                 name==VxEnum ||
    4389                                 name==VyEnum ||
    4390                                 name==InversionVxObsEnum ||
    4391                                 name==InversionVyObsEnum ||
    4392                                 name==FrictionCoefficientEnum ||
    4393                                 name==GradientEnum ||
    4394                                 name==OldGradientEnum
    4395                 ){
    4396                 return true;
    4397         }
    4398         else return false;
    4399 }
    4400 /*}}}*/
    4401 /*FUNCTION Tria::IsOnBed {{{1*/
    4402 bool Tria::IsOnBed(){
    4403        
    4404         bool onbed;
    4405         inputs->GetParameterValue(&onbed,MeshElementonbedEnum);
    4406         return onbed;
    4407 }
    4408 /*}}}*/
    4409 /*FUNCTION Tria::IsOnShelf {{{1*/
    4410 bool   Tria::IsOnShelf(){
    4411 
    4412         bool shelf;
    4413         inputs->GetParameterValue(&shelf,MaskElementonfloatingiceEnum);
    4414         return shelf;
    4415 }
    4416 /*}}}*/
    4417 /*FUNCTION Tria::IsNodeOnShelf {{{1*/
    4418 bool   Tria::IsNodeOnShelf(){
    4419 
    4420         int  i;
    4421         bool shelf=false;
    4422 
    4423         for(i=0;i<3;i++){
    4424                 if (nodes[i]->IsOnShelf()){
    4425                         shelf=true;
    4426                         break;
    4427                 }
    4428         }
    4429         return shelf;
    4430 }
    4431 /*}}}*/
    4432 /*FUNCTION Tria::IsNodeOnShelfFromFlags {{{1*/
    4433 bool   Tria::IsNodeOnShelfFromFlags(double* flags){
    4434 
    4435         int  i;
    4436         bool shelf=false;
    4437 
    4438         for(i=0;i<3;i++){
    4439                 if (flags[nodes[i]->Sid()]){
    4440                         shelf=true;
    4441                         break;
    4442                 }
    4443         }
    4444         return shelf;
    4445 }
    4446 /*}}}*/
    4447 /*FUNCTION Tria::IsOnWater {{{1*/
    4448 bool   Tria::IsOnWater(){
    4449 
    4450         bool water;
    4451         inputs->GetParameterValue(&water,MaskElementonwaterEnum);
    4452         return water;
    4453 }
    4454 /*}}}*/
    4455 /*FUNCTION Tria::MassFlux {{{1*/
    4456 double Tria::MassFlux( double* segment,bool process_units){
    4457 
    4458         const int    numdofs=2;
    4459 
    4460         int        i;
    4461         double     mass_flux=0;
     5330        int        i,j,ig,dim;
     5331        double     vx,vy,Jdettria;
    44625332        double     xyz_list[NUMVERTICES][3];
    4463         double     normal[2];
    4464         double     length,rho_ice;
    4465         double     x1,y1,x2,y2,h1,h2;
    4466         double     vx1,vx2,vy1,vy2;
    4467         GaussTria* gauss_1=NULL;
    4468         GaussTria* gauss_2=NULL;
    4469 
    4470         /*Get material parameters :*/
    4471         rho_ice=matpar->GetRhoIce();
    4472 
    4473         /*First off, check that this segment belongs to this element: */
    4474         if ((int)*(segment+4)!=this->id)_error_("%s%i%s%i","error message: segment with id ",(int)*(segment+4)," does not belong to element with id:",this->id);
    4475 
    4476         /*Recover segment node locations: */
    4477         x1=*(segment+0); y1=*(segment+1); x2=*(segment+2); y2=*(segment+3);
    4478 
    4479         /*Get xyz list: */
     5333        double     B[2][NUMVERTICES];
     5334        double     Bprime[2][NUMVERTICES];
     5335        double     DL[2][2]={0.0};
     5336        double     DL_scalar;
     5337        GaussTria  *gauss=NULL;
     5338
     5339        /*Initialize Element matrix*/
     5340        ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
     5341
     5342        /*Retrieve all inputs and parameters*/
    44805343        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    4481 
    4482         /*get area coordinates of 0 and 1 locations: */
    4483         gauss_1=new GaussTria();
    4484         gauss_1->GaussFromCoords(x1,y1,&xyz_list[0][0]);
    4485         gauss_2=new GaussTria();
    4486         gauss_2->GaussFromCoords(x2,y2,&xyz_list[0][0]);
    4487 
    4488         normal[0]=cos(atan2(x1-x2,y2-y1));
    4489         normal[1]=sin(atan2(x1-x2,y2-y1));
    4490 
    4491         length=sqrt(pow(x2-x1,2.0)+pow(y2-y1,2));
    4492 
    4493         Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
     5344        this->parameters->FindParam(&dim,MeshDimensionEnum);
    44945345        Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
    44955346        Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
    44965347
    4497         thickness_input->GetParameterValue(&h1, gauss_1);
    4498         thickness_input->GetParameterValue(&h2, gauss_2);
    4499         vx_input->GetParameterValue(&vx1,gauss_1);
    4500         vx_input->GetParameterValue(&vx2,gauss_2);
    4501         vy_input->GetParameterValue(&vy1,gauss_1);
    4502         vy_input->GetParameterValue(&vy2,gauss_2);
    4503 
    4504         mass_flux= rho_ice*length*( 
    4505                                 (ONETHIRD*(h1-h2)*(vx1-vx2)+0.5*h2*(vx1-vx2)+0.5*(h1-h2)*vx2+h2*vx2)*normal[0]+
    4506                                 (ONETHIRD*(h1-h2)*(vy1-vy2)+0.5*h2*(vy1-vy2)+0.5*(h1-h2)*vy2+h2*vy2)*normal[1]
    4507                                 );
    4508 
    4509         /*Process units: */
    4510         mass_flux=UnitConversion(mass_flux,IuToExtEnum,MassFluxEnum);
    4511 
    4512         /*clean up and return:*/
    4513         delete gauss_1;
    4514         delete gauss_2;
    4515         return mass_flux;
    4516 }
    4517 /*}}}*/
    4518 /*FUNCTION Tria::MaxAbsVx{{{1*/
    4519 void  Tria::MaxAbsVx(double* pmaxabsvx, bool process_units){
    4520 
    4521         /*Get maximum:*/
    4522         double maxabsvx=this->inputs->MaxAbs(VxEnum);
    4523 
    4524         /*process units if requested: */
    4525         if(process_units) maxabsvx=UnitConversion(maxabsvx,IuToExtEnum,VxEnum);
    4526 
    4527         /*Assign output pointers:*/
    4528         *pmaxabsvx=maxabsvx;
    4529 }
    4530 /*}}}*/
    4531 /*FUNCTION Tria::MaxAbsVy{{{1*/
    4532 void  Tria::MaxAbsVy(double* pmaxabsvy, bool process_units){
    4533 
    4534         /*Get maximum:*/
    4535         double maxabsvy=this->inputs->MaxAbs(VyEnum);
    4536 
    4537         /*process units if requested: */
    4538         if(process_units) maxabsvy=UnitConversion(maxabsvy,IuToExtEnum,VyEnum);
    4539 
    4540         /*Assign output pointers:*/
    4541         *pmaxabsvy=maxabsvy;
    4542 }
    4543 /*}}}*/
    4544 /*FUNCTION Tria::MaxAbsVz{{{1*/
    4545 void  Tria::MaxAbsVz(double* pmaxabsvz, bool process_units){
    4546 
    4547         /*Get maximum:*/
    4548         double maxabsvz=this->inputs->MaxAbs(VzEnum);
    4549 
    4550         /*process units if requested: */
    4551         if(process_units) maxabsvz=UnitConversion(maxabsvz,IuToExtEnum,VyEnum);
    4552 
    4553         /*Assign output pointers:*/
    4554         *pmaxabsvz=maxabsvz;
    4555 }
    4556 /*}}}*/
    4557 /*FUNCTION Tria::MaxVel{{{1*/
    4558 void  Tria::MaxVel(double* pmaxvel, bool process_units){
    4559 
    4560         /*Get maximum:*/
    4561         double maxvel=this->inputs->Max(VelEnum);
    4562 
    4563         /*process units if requested: */
    4564         if(process_units) maxvel=UnitConversion(maxvel,IuToExtEnum,VelEnum);
    4565 
    4566         /*Assign output pointers:*/
    4567         *pmaxvel=maxvel;
    4568 }
    4569 /*}}}*/
    4570 /*FUNCTION Tria::MaxVx{{{1*/
    4571 void  Tria::MaxVx(double* pmaxvx, bool process_units){
    4572 
    4573         /*Get maximum:*/
    4574         double maxvx=this->inputs->Max(VxEnum);
    4575 
    4576         /*process units if requested: */
    4577         if(process_units) maxvx=UnitConversion(maxvx,IuToExtEnum,VxEnum);
    4578 
    4579         /*Assign output pointers:*/
    4580         *pmaxvx=maxvx;
    4581 }
    4582 /*}}}*/
    4583 /*FUNCTION Tria::MaxVy{{{1*/
    4584 void  Tria::MaxVy(double* pmaxvy, bool process_units){
    4585 
    4586         /*Get maximum:*/
    4587         double maxvy=this->inputs->Max(VyEnum);
    4588 
    4589         /*process units if requested: */
    4590         if(process_units) maxvy=UnitConversion(maxvy,IuToExtEnum,VyEnum);
    4591 
    4592         /*Assign output pointers:*/
    4593         *pmaxvy=maxvy;
    4594 
    4595 }
    4596 /*}}}*/
    4597 /*FUNCTION Tria::MaxVz{{{1*/
    4598 void  Tria::MaxVz(double* pmaxvz, bool process_units){
    4599 
    4600         /*Get maximum:*/
    4601         double maxvz=this->inputs->Max(VzEnum);
    4602 
    4603         /*process units if requested: */
    4604         if(process_units) maxvz=UnitConversion(maxvz,IuToExtEnum,VzEnum);
    4605 
    4606         /*Assign output pointers:*/
    4607         *pmaxvz=maxvz;
    4608 }
    4609 /*}}}*/
    4610 /*FUNCTION Tria::MigrateGroundingLine{{{1*/
    4611 void  Tria::MigrateGroundingLine(void){
    4612 
    4613 
    4614         double *values         = NULL;
    4615         double  h[3],s[3],b[3],ba[3];
    4616         double  bed_hydro;
    4617         double  rho_water,rho_ice,density;
    4618         int     isonshelf[3];
    4619         int     i;
    4620         bool    elementonshelf = false;
    4621 
    4622 
    4623         /*Recover info at the vertices: */
    4624         Input* surface_input =inputs->GetInput(SurfaceEnum); _assert_(surface_input);
    4625         Input* bed_input     =inputs->GetInput(BedEnum);     _assert_(bed_input);
    4626         if((surface_input->Enum()!=TriaVertexInputEnum) | (bed_input->Enum()!=TriaVertexInputEnum))_error_(" not supported yet for bed and surface interpolations not P1!");
    4627 
    4628         GetParameterListOnVertices(&h[0],ThicknessEnum);
    4629         GetParameterListOnVertices(&s[0],SurfaceEnum);
    4630         GetParameterListOnVertices(&b[0],BedEnum);
    4631         GetParameterListOnVertices(&ba[0],BathymetryEnum);
    4632         for(i=0;i<3;i++){
    4633                 isonshelf[i]=nodes[i]->IsOnShelf();
    4634                 if((nodes[i]->Sid()+1)==36)printf("MigrateGroundingLine: El %i Node %i shelf status %i\n",this->Id(),nodes[i]->Sid()+1,isonshelf[i]);
    4635         }
    4636 
    4637         /*material parameters: */
    4638         rho_water=matpar->GetRhoWater();
    4639         rho_ice=matpar->GetRhoIce();
    4640         density=rho_ice/rho_water;
    4641 
    4642         /*go through vertices, and update inputs, considering them to be TriaVertex type: */
    4643         for(i=0;i<3;i++){
    4644                 if (isonshelf[i]){
    4645                         /*This node is on the shelf. See if its bed is going under the bathymetry: */
    4646                         if(b[i]<=ba[i]){ //<= because Neff being 0 when b=ba, drag will be 0 anyway.
    4647                                 /*The ice shelf is getting grounded, the thickness is the same, so just update the bed to stick to the bathymetry and elevate the surface accordingly: */
    4648                                 b[i]=ba[i];
    4649                                 s[i]=b[i]+h[i];
    4650                                 if((nodes[i]->Sid()+1)==36)printf("MigrateGroundingLine: El %i Node %i grounding %g %g %g\n",this->Id(),nodes[i]->Sid()+1,b[i],ba[i],s[i]);
    4651                         }
    4652                         else{
    4653                                 /*do nothing, we are still floating.*/
    4654                                 if((nodes[i]->Sid()+1)==36)printf("MigrateGroundingLine: El %i Node %i still floating %g %g %g\n",this->Id(),nodes[i]->Sid()+1,b[i],ba[i],s[i]);
    4655                         }
    4656                 }
    4657                 else{
    4658                         /*This node is on the sheet, near the grounding line. See if wants to unground. To
    4659                          * do so, we compute the hydrostatic bed, and if it is > bathymetry, then we unground: */
    4660                         bed_hydro=-density*h[i];
    4661                         if (bed_hydro>ba[i]){
    4662                                 /*We are now floating, bed and surface are determined from hydrostatic equilibrium: */
    4663                                 s[i]=(1-density)*h[i];
    4664                                 b[i]=-density*h[i];
    4665                                 printf("%i\n",nodes[i]->Sid()+1);
    4666                                 if((nodes[i]->Sid()+1)==36)printf("MigrateGroundingLine: El %i Node %i floating %g %g %g\n",this->Id(),nodes[i]->Sid()+1,b[i],ba[i],s[i]);
    4667                         }
    4668                         else{
    4669                                 /*do nothing, we are still grounded.*/
    4670                                 if((nodes[i]->Sid()+1)==36)printf("MigrateGroundingLine: El %i Node %i still grounded %g %g %g\n",this->Id(),nodes[i]->Sid()+1,b[i],ba[i],s[i]);
    4671                         }
    4672                 }
    4673         }
    4674 
    4675         /*Surface and bed are updated. Update inputs:*/
    4676         surface_input->GetValuesPtr(&values,NULL); for(i=0;i<3;i++)values[i]=s[i];
    4677         bed_input->GetValuesPtr(&values,NULL);     for(i=0;i<3;i++)values[i]=b[i];
     5348        /*Start looping on the number of gaussian points:*/
     5349        gauss=new GaussTria(2);
     5350        for (ig=gauss->begin();ig<gauss->end();ig++){
     5351
     5352                gauss->GaussPoint(ig);
     5353
     5354                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
     5355                /*WARNING: B and Bprime are inverted compared to usual prognostic!!!!*/
     5356                GetBPrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
     5357                GetBprimePrognostic(&B[0][0], &xyz_list[0][0], gauss);
     5358
     5359                vx_input->GetParameterValue(&vx,gauss);
     5360                vy_input->GetParameterValue(&vy,gauss);
     5361
     5362                DL_scalar=-gauss->weight*Jdettria;
     5363                DL[0][0]=DL_scalar*vx;
     5364                DL[1][1]=DL_scalar*vy;
     5365
     5366                TripleMultiply( &B[0][0],2,numdof,1,
     5367                                        &DL[0][0],2,2,0,
     5368                                        &Bprime[0][0],2,numdof,0,
     5369                                        &Ke->values[0],1);
     5370        }
     5371
     5372        /*Clean up and return*/
     5373        delete gauss;
     5374        return Ke;
     5375}
     5376/*}}}*/
     5377/*FUNCTION Tria::CreatePVectorBalancethickness{{{1*/
     5378ElementVector* Tria::CreatePVectorBalancethickness(void){
     5379
     5380        switch(GetElementType()){
     5381                case P1Enum:
     5382                        return CreatePVectorBalancethickness_CG();
     5383                        break;
     5384                case P1DGEnum:
     5385                        return CreatePVectorBalancethickness_DG();
     5386                default:
     5387                        _error_("Element type %s not supported yet",EnumToStringx(GetElementType()));
     5388        }
     5389}
     5390/*}}}*/
     5391/*FUNCTION Tria::CreatePVectorBalancethickness_CG{{{1*/
     5392ElementVector* Tria::CreatePVectorBalancethickness_CG(void){
     5393
     5394        /*Constants*/
     5395        const int    numdof=NDOF1*NUMVERTICES;
    46785396       
    4679         for(i=0;i<3;i++){
    4680                 isonshelf[i]=nodes[i]->IsOnShelf();
    4681                 if((nodes[i]->Sid()+1)==36)printf("MigrateGroundingLine: El %i Node %i second shelf status %i\n",this->Id(),nodes[i]->Sid()+1,isonshelf[i]);
    4682         }
    4683 }
    4684 /*}}}*/
    4685 /*FUNCTION Tria::MinVel{{{1*/
    4686 void  Tria::MinVel(double* pminvel, bool process_units){
    4687 
    4688         /*Get minimum:*/
    4689         double minvel=this->inputs->Min(VelEnum);
    4690 
    4691         /*process units if requested: */
    4692         if(process_units) minvel=UnitConversion(minvel,IuToExtEnum,VelEnum);
    4693 
    4694         /*Assign output pointers:*/
    4695         *pminvel=minvel;
    4696 }
    4697 /*}}}*/
    4698 /*FUNCTION Tria::MinVx{{{1*/
    4699 void  Tria::MinVx(double* pminvx, bool process_units){
    4700 
    4701         /*Get minimum:*/
    4702         double minvx=this->inputs->Min(VxEnum);
    4703 
    4704         /*process units if requested: */
    4705         if(process_units) minvx=UnitConversion(minvx,IuToExtEnum,VxEnum);
    4706 
    4707         /*Assign output pointers:*/
    4708         *pminvx=minvx;
    4709 }
    4710 /*}}}*/
    4711 /*FUNCTION Tria::MinVy{{{1*/
    4712 void  Tria::MinVy(double* pminvy, bool process_units){
    4713 
    4714         /*Get minimum:*/
    4715         double minvy=this->inputs->Min(VyEnum);
    4716 
    4717         /*process units if requested: */
    4718         if(process_units) minvy=UnitConversion(minvy,IuToExtEnum,VyEnum);
    4719 
    4720         /*Assign output pointers:*/
    4721         *pminvy=minvy;
    4722 }
    4723 /*}}}*/
    4724 /*FUNCTION Tria::MinVz{{{1*/
    4725 void  Tria::MinVz(double* pminvz, bool process_units){
    4726 
    4727         /*Get minimum:*/
    4728         double minvz=this->inputs->Min(VzEnum);
    4729 
    4730         /*process units if requested: */
    4731         if(process_units) minvz=UnitConversion(minvz,IuToExtEnum,VzEnum);
    4732 
    4733         /*Assign output pointers:*/
    4734         *pminvz=minvz;
    4735 }
    4736 /*}}}*/
    4737 /*FUNCTION Tria::MyRank {{{1*/
    4738 int    Tria::MyRank(void){
    4739         extern int my_rank;
    4740         return my_rank;
    4741 }
    4742 /*}}}*/
    4743 /*FUNCTION Tria::NodalValue {{{1*/
    4744 int    Tria::NodalValue(double* pvalue, int index, int natureofdataenum,bool process_units){
    4745 
    4746         int i;
    4747         int found=0;
    4748         double value;
    4749         Input* data=NULL;
    4750         GaussTria *gauss                            = NULL;
    4751 
    4752         /*First, serarch the input: */
    4753         data=inputs->GetInput(natureofdataenum);
    4754 
    4755         /*figure out if we have the vertex id: */
    4756         found=0;
    4757         for(i=0;i<NUMVERTICES;i++){
    4758                 if(index==nodes[i]->GetVertexId()){
    4759                         /*Do we have natureofdataenum in our inputs? :*/
    4760                         if(data){
    4761                                 /*ok, we are good. retrieve value of input at vertex :*/
    4762                                 gauss=new GaussTria(); gauss->GaussVertex(i);
    4763                                 data->GetParameterValue(&value,gauss);
    4764                                 found=1;
    4765                                 break;
    4766                         }
    4767                 }
    4768         }
    4769 
    4770         if(found)*pvalue=value;
    4771         return found;
    4772 }
    4773 /*}}}*/
    4774 /*FUNCTION Tria::PatchFill{{{1*/
    4775 void  Tria::PatchFill(int* prow, Patch* patch){
    4776 
    4777         int i,row;
    4778         int vertices_ids[3];
    4779 
    4780         /*recover pointer: */
    4781         row=*prow;
    4782                
    4783         for(i=0;i<3;i++) vertices_ids[i]=nodes[i]->GetVertexId(); //vertices id start at column 3 of the patch.
    4784 
    4785         for(i=0;i<this->results->Size();i++){
    4786                 ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i);
    4787 
    4788                 /*For this result,fill the information in the Patch object (element id + vertices ids), and then hand
    4789                  *it to the result object, to fill the rest: */
    4790                 patch->fillelementinfo(row,this->sid+1,vertices_ids,3);
    4791                 elementresult->PatchFill(row,patch);
    4792 
    4793                 /*increment rower: */
    4794                 row++;
    4795         }
    4796 
    4797         /*Assign output pointers:*/
    4798         *prow=row;
    4799 }
    4800 /*}}}*/
    4801 /*FUNCTION Tria::PatchSize{{{1*/
    4802 void  Tria::PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes){
    4803 
    4804         int     i;
    4805         int     numrows     = 0;
    4806         int     numnodes    = 0;
    4807         int     temp_numnodes    = 0;
    4808 
    4809         /*Go through all the results objects, and update the counters: */
    4810         for (i=0;i<this->results->Size();i++){
    4811                 ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i);
    4812                 /*first, we have one more result: */
    4813                 numrows++;
    4814                 /*now, how many vertices and how many nodal values for this result? :*/
    4815                 temp_numnodes=elementresult->NumberOfNodalValues(); //ask result object.
    4816                 if(temp_numnodes>numnodes)numnodes=temp_numnodes;
    4817         }
    4818 
    4819         /*Assign output pointers:*/
    4820         *pnumrows=numrows;
    4821         *pnumvertices=NUMVERTICES;
    4822         *pnumnodes=numnodes;
    4823 }
    4824 /*}}}*/
    4825 /*FUNCTION Tria::PotentialSheetUngrounding{{{1*/
    4826 void  Tria::PotentialSheetUngrounding(Vec potential_sheet_ungrounding){
    4827 
    4828 
    4829         double *values         = NULL;
    4830         double  h[3],s[3],b[3],ba[3];
    4831         double  bed_hydro;
    4832         double  rho_water,rho_ice,density;
    4833         int     i;
    4834         bool    elementonshelf = false;
    4835 
    4836         /*Recover info at the vertices: */
    4837         Input* surface_input =inputs->GetInput(SurfaceEnum); _assert_(surface_input);
    4838         Input* bed_input     =inputs->GetInput(BedEnum);     _assert_(bed_input);
    4839         if((surface_input->Enum()!=TriaVertexInputEnum) | (bed_input->Enum()!=TriaVertexInputEnum))_error_(" not supported yet for bed and surface interpolations not P1!");
    4840 
    4841         GetParameterListOnVertices(&h[0],ThicknessEnum);
    4842         GetParameterListOnVertices(&s[0],SurfaceEnum);
    4843         GetParameterListOnVertices(&b[0],BedEnum);
    4844         GetParameterListOnVertices(&ba[0],BathymetryEnum);
    4845 
    4846         /*material parameters: */
    4847         rho_water=matpar->GetRhoWater();
    4848         rho_ice=matpar->GetRhoIce();
    4849         density=rho_ice/rho_water;
    4850 
    4851         /*go through vertices, and figure out which ones are on the ice sheet, and want to unground: */
    4852         for(i=0;i<3;i++){
    4853                 if (!nodes[i]->IsOnShelf()){
    4854                        
    4855                         /*This node is on the sheet, near the grounding line. See if wants to unground. To
    4856                          * do so, we compute the hydrostatic bed, and if it is > bathymetry, then we unground: */
    4857                         bed_hydro=-density*h[i];
    4858                         if (bed_hydro>ba[i]){
    4859                                 /*ok, this node wants to unground. flag it: */
    4860                                 VecSetValue(potential_sheet_ungrounding,nodes[i]->Sid(),1,INSERT_VALUES);
    4861                         }
    4862                 }
    4863         }
    4864 }
    4865 /*}}}*/
    4866 /*FUNCTION Tria::ProcessResultsUnits{{{1*/
    4867 void  Tria::ProcessResultsUnits(void){
    4868 
    4869         int i;
    4870 
    4871         for(i=0;i<this->results->Size();i++){
    4872                 ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i);
    4873                 elementresult->ProcessUnits(this->parameters);
    4874         }
    4875 }
    4876 /*}}}*/
    4877 /*FUNCTION Tria::RheologyBbarx{{{1*/
    4878 double Tria::RheologyBbarx(void){
    4879 
    4880         return this->matice->GetBbar();
    4881 
    4882 }
    4883 /*}}}*/
    4884 /*FUNCTION Tria::RequestedOutput{{{1*/
    4885 void Tria::RequestedOutput(int output_enum,int step,double time){
    4886 
    4887         if(IsInput(output_enum)){
    4888                 /*just transfer this input to results, and we are done: */
    4889                 InputToResult(output_enum,step,time);
    4890         }
    4891         else{
    4892                 /*this input does not exist, compute it, and then transfer to results: */
    4893                 switch(output_enum){
    4894                         default:
    4895                                 /*do nothing, no need to derail the computation because one of the outputs requested cannot be found: */
    4896                                 break;
    4897                 }
    4898         }
    4899 
    4900 }
    4901 /*}}}*/
    4902 /*FUNCTION Tria::SetClone {{{1*/
    4903 void  Tria::SetClone(int* minranks){
    4904 
    4905         _error_("not implemented yet");
    4906 }
    4907 /*}}}1*/
    4908 /*FUNCTION Tria::SetCurrentConfiguration {{{1*/
    4909 void  Tria::SetCurrentConfiguration(Elements* elementsin, Loads* loadsin, DataSet* nodesin, Materials* materialsin, Parameters* parametersin){
     5397        /*Intermediaries */
     5398        int        i,j,ig;
     5399        double     xyz_list[NUMVERTICES][3];
     5400        double     dhdt_g,basal_melting_g,surface_mass_balance_g,Jdettria;
     5401        double     L[NUMVERTICES];
     5402        GaussTria* gauss=NULL;
     5403
     5404        /*Initialize Element vector*/
     5405        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
     5406
     5407        /*Retrieve all inputs and parameters*/
     5408        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     5409        Input* surface_mass_balance_input=inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(surface_mass_balance_input);
     5410        Input* basal_melting_input=inputs->GetInput(BasalforcingsMeltingRateEnum);          _assert_(basal_melting_input);
     5411        Input* dhdt_input=inputs->GetInput(BalancethicknessThickeningRateEnum);             _assert_(dhdt_input);
    49105412       
    4911         /*go into parameters and get the analysis_counter: */
    4912         int analysis_counter;
    4913         parametersin->FindParam(&analysis_counter,AnalysisCounterEnum);
    4914 
    4915         /*Get Element type*/
    4916         this->element_type=this->element_type_list[analysis_counter];
    4917 
    4918         /*Pick up nodes*/
    4919         if(this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp();
    4920         else this->nodes=NULL;
    4921 
    4922 }
    4923 /*}}}*/
    4924 /*FUNCTION Tria::ShelfSync{{{1*/
    4925 void  Tria::ShelfSync(void){
    4926 
    4927 
    4928         double *values         = NULL;
    4929         double  h[3],s[3],b[3],ba[3];
    4930         double  bed_hydro;
    4931         double  rho_water,rho_ice,density;
    4932         int     i;
    4933         bool    elementonshelf = false;
    4934 
    4935         /*melting rate at the grounding line: */
    4936         double  yts;
    4937         int     swap;
    4938         double  gl_melting_rate;
    4939 
    4940         /*recover parameters: */
    4941         parameters->FindParam(&yts,ConstantsYtsEnum);
    4942         parameters->FindParam(&gl_melting_rate,GroundinglineMeltingRateEnum);
    4943 
    4944         /*Recover info at the vertices: */
    4945         Input* surface_input =inputs->GetInput(SurfaceEnum); _assert_(surface_input);
    4946         Input* bed_input     =inputs->GetInput(BedEnum);     _assert_(bed_input);
    4947         if((surface_input->Enum()!=TriaVertexInputEnum) | (bed_input->Enum()!=TriaVertexInputEnum))_error_(" not supported yet for bed and surface interpolations not P1!");
    4948 
    4949         GetParameterListOnVertices(&h[0],ThicknessEnum);
    4950         GetParameterListOnVertices(&s[0],SurfaceEnum);
    4951         GetParameterListOnVertices(&b[0],BedEnum);
    4952         GetParameterListOnVertices(&ba[0],BathymetryEnum);
    4953 
    4954         /*material parameters: */
    4955         rho_water=matpar->GetRhoWater();
    4956         rho_ice=matpar->GetRhoIce();
    4957         density=rho_ice/rho_water;
    4958 
    4959         /*go through vertices, and update inputs, considering them to be TriaVertex type: */
    4960         for(i=0;i<3;i++){
    4961                 if(b[i]==ba[i]){
    4962                                
    4963                         nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,false));
    4964                         nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,true));
    4965                 }
    4966                 else{
    4967                         nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,true));
    4968                         nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,false));
    4969 
    4970                 }
    4971         }
    4972 
    4973         /*Now, update  shelf status of element. An element can only be on shelf if all its nodes are on shelf: */
    4974         swap=0;
    4975         elementonshelf=false;
    4976         for(i=0;i<3;i++){
    4977                 if(nodes[i]->IsOnShelf()){
    4978                         elementonshelf=true;
    4979                         break;
    4980                 }
    4981         }
    4982         if(!this->IsOnShelf() && elementonshelf==true)swap=1;
    4983     this->inputs->AddInput(new BoolInput(MaskElementonfloatingiceEnum,elementonshelf));
    4984        
    4985     /*If this element just  became ungrounded, set its basal melting rate at 50 m/yr:*/
    4986         if(swap){
    4987                 Input* basal_melting_rate_input     =inputs->GetInput(BasalforcingsMeltingRateEnum);     _assert_(basal_melting_rate_input);
    4988                 basal_melting_rate_input->GetValuesPtr(&values,NULL); for(i=0;i<3;i++)values[i]=gl_melting_rate/yts;
    4989         }
    4990 
    4991 
    4992 }
    4993 /*}}}*/
    4994 /*FUNCTION Tria::SoftMigration{{{1*/
    4995 void  Tria::SoftMigration(double* sheet_ungrounding){
    4996 
    4997 
    4998         double *values         = NULL;
    4999         double  h[3],s[3],b[3],ba[3];
    5000         double  bed_hydro;
    5001         double  rho_water,rho_ice,density;
    5002         int     i;
    5003         bool    elementonshelf = false;
    5004 
    5005         /*Recover info at the vertices: */
    5006         Input* surface_input =inputs->GetInput(SurfaceEnum); _assert_(surface_input);
    5007         Input* bed_input     =inputs->GetInput(BedEnum);     _assert_(bed_input);
    5008         if((surface_input->Enum()!=TriaVertexInputEnum) | (bed_input->Enum()!=TriaVertexInputEnum))_error_(" not supported yet for bed and surface interpolations not P1!");
    5009 
    5010         GetParameterListOnVertices(&h[0],ThicknessEnum);
    5011         GetParameterListOnVertices(&s[0],SurfaceEnum);
    5012         GetParameterListOnVertices(&b[0],BedEnum);
    5013         GetParameterListOnVertices(&ba[0],BathymetryEnum);
    5014 
    5015         /*material parameters: */
    5016         rho_water=matpar->GetRhoWater();
    5017         rho_ice=matpar->GetRhoIce();
    5018         density=rho_ice/rho_water;
    5019        
    5020         /*go through vertices, and update inputs, considering them to be TriaVertex type: */
    5021         for(i=0;i<3;i++){
    5022                 if (nodes[i]->IsOnShelf()){
    5023                         /*This node is on the shelf. See if its bed is going under the bathymetry: */
    5024                         if(b[i]<=ba[i]){ //<= because Neff being 0 when b=ba, drag will be 0 anyway.
    5025                                 /*The ice shelf is getting grounded, the thickness is the same, so just update the bed to stick to the bathymetry and elevate the surface accordingly: */
    5026                                 b[i]=ba[i];
    5027                                 s[i]=b[i]+h[i];
    5028                         }
    5029                         else{
    5030                                 /*do nothing, we are still floating.*/
    5031                         }
    5032                 }
    5033                 else{
    5034                         /*This node is on the sheet, near the grounding line. See if wants to unground. To
    5035                          * do so, we compute the hydrostatic bed, and if it is > bathymetry, then we unground: */
    5036                         bed_hydro=-density*h[i];
    5037                         if (bed_hydro>ba[i]){
    5038 
    5039                                 /*Now, are we connected to the grounding line, if so, go forward, otherwise, bail out: */
    5040                                 if(sheet_ungrounding[nodes[i]->Sid()]){
    5041                                         /*We are now floating, bed and surface are determined from hydrostatic equilibrium: */
    5042                                         s[i]=(1-density)*h[i];
    5043                                         b[i]=-density*h[i];
    5044                                 }
    5045                         }
    5046                         else{
    5047                                 /*do nothing, we are still grounded.*/
    5048                         }
    5049                 }
    5050         }
    5051 
    5052         /*Surface and bed are updated. Update inputs:*/   
    5053         surface_input->GetValuesPtr(&values,NULL); for(i=0;i<3;i++)values[i]=s[i];
    5054         bed_input->GetValuesPtr(&values,NULL);     for(i=0;i<3;i++)values[i]=b[i];
    5055 
    5056         }
    5057 /*}}}*/
    5058 /*FUNCTION Tria::SurfaceArea {{{1*/
    5059 double Tria::SurfaceArea(void){
    5060 
    5061         int    i;
    5062         double S;
    5063         double normal[3];
    5064         double v13[3],v23[3];
    5065         double xyz_list[NUMVERTICES][3];
    5066 
    5067         /*If on water, return 0: */
    5068         if(IsOnWater())return 0;
    5069 
     5413        /* Start  looping on the number of gaussian points: */
     5414        gauss=new GaussTria(2);
     5415        for(ig=gauss->begin();ig<gauss->end();ig++){
     5416
     5417                gauss->GaussPoint(ig);
     5418
     5419                surface_mass_balance_input->GetParameterValue(&surface_mass_balance_g,gauss);
     5420                basal_melting_input->GetParameterValue(&basal_melting_g,gauss);
     5421                dhdt_input->GetParameterValue(&dhdt_g,gauss);
     5422
     5423                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
     5424                GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
     5425
     5426                for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(surface_mass_balance_g-basal_melting_g-dhdt_g)*L[i];
     5427        }
     5428
     5429        /*Clean up and return*/
     5430        delete gauss;
     5431        return pe;
     5432}
     5433/*}}}*/
     5434/*FUNCTION Tria::CreatePVectorBalancethickness_DG {{{1*/
     5435ElementVector* Tria::CreatePVectorBalancethickness_DG(void){
     5436
     5437        /*Constants*/
     5438        const int    numdof=NDOF1*NUMVERTICES;
     5439
     5440        /*Intermediaries */
     5441        int        i,j,ig;
     5442        double     xyz_list[NUMVERTICES][3];
     5443        double     basal_melting_g,surface_mass_balance_g,dhdt_g,Jdettria;
     5444        double     L[NUMVERTICES];
     5445        GaussTria* gauss=NULL;
     5446
     5447        /*Initialize Element vector*/
     5448        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
     5449
     5450        /*Retrieve all inputs and parameters*/
    50705451        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    5071 
    5072         for (i=0;i<3;i++){
    5073                 v13[i]=xyz_list[0][i]-xyz_list[2][i];
    5074                 v23[i]=xyz_list[1][i]-xyz_list[2][i];
    5075         }
    5076 
    5077         normal[0]=v13[1]*v23[2]-v13[2]*v23[1];
    5078         normal[1]=v13[2]*v23[0]-v13[0]*v23[2];
    5079         normal[2]=v13[0]*v23[1]-v13[1]*v23[0];
    5080 
    5081         S = 0.5 * sqrt(pow(normal[0],(double)2)+pow(normal[1],(double)2)+pow(normal[2],(double)2));
    5082 
    5083         /*Return: */
    5084         return S;
    5085 }
    5086 /*}}}*/
    5087 /*FUNCTION Tria::SurfaceNormal{{{1*/
    5088 void Tria::SurfaceNormal(double* surface_normal, double xyz_list[3][3]){
    5089 
    5090         int i;
    5091         double v13[3],v23[3];
    5092         double normal[3];
    5093         double normal_norm;
    5094 
    5095         for (i=0;i<3;i++){
    5096                 v13[i]=xyz_list[0][i]-xyz_list[2][i];
    5097                 v23[i]=xyz_list[1][i]-xyz_list[2][i];
    5098         }
    5099 
    5100         normal[0]=v13[1]*v23[2]-v13[2]*v23[1];
    5101         normal[1]=v13[2]*v23[0]-v13[0]*v23[2];
    5102         normal[2]=v13[0]*v23[1]-v13[1]*v23[0];
    5103 
    5104         normal_norm=sqrt( pow(normal[0],(double)2)+pow(normal[1],(double)2)+pow(normal[2],(double)2) );
    5105 
    5106         *(surface_normal)=normal[0]/normal_norm;
    5107         *(surface_normal+1)=normal[1]/normal_norm;
    5108         *(surface_normal+2)=normal[2]/normal_norm;
    5109 }
    5110 /*}}}*/
    5111 /*FUNCTION Tria::TimeAdapt{{{1*/
    5112 double  Tria::TimeAdapt(void){
    5113 
    5114         /*intermediary: */
    5115         int    i;
    5116         double C,dt;
    5117         double dx,dy;
    5118         double maxx,minx;
    5119         double maxy,miny;
    5120         double maxabsvx,maxabsvy;
    5121         double xyz_list[NUMVERTICES][3];
    5122 
    5123         /*get CFL coefficient:*/
    5124         this->parameters->FindParam(&C,TimesteppingCflCoefficientEnum);
    5125 
    5126         /*Get for Vx and Vy, the max of abs value: */
    5127         this->MaxAbsVx(&maxabsvx,false);
    5128         this->MaxAbsVy(&maxabsvy,false);
    5129 
    5130         /* Get node coordinates and dof list: */
    5131         GetVerticesCoordinates(&xyz_list[0][0], this->nodes, NUMVERTICES);
    5132 
    5133         minx=xyz_list[0][0];
    5134         maxx=xyz_list[0][0];
    5135         miny=xyz_list[0][1];
    5136         maxy=xyz_list[0][1];
    5137 
    5138         for(i=1;i<NUMVERTICES;i++){
    5139                 if (xyz_list[i][0]<minx)minx=xyz_list[i][0];
    5140                 if (xyz_list[i][0]>maxx)maxx=xyz_list[i][0];
    5141                 if (xyz_list[i][1]<miny)miny=xyz_list[i][1];
    5142                 if (xyz_list[i][1]>maxy)maxy=xyz_list[i][1];
    5143         }
    5144         dx=maxx-minx;
    5145         dy=maxy-miny;
    5146 
    5147         /*CFL criterion: */
    5148         dt=C/(maxabsvy/dx+maxabsvy/dy);
    5149 
    5150         return dt;
    5151 }
    5152 /*}}}*/
    5153 /*FUNCTION Tria::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type){{{1*/
    5154 void Tria::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type){ //i is the element index
    5155 
    5156         /*Intermediaries*/
    5157         int    i,j;
    5158         int    tria_node_ids[3];
    5159         int    tria_vertex_ids[3];
    5160         int    tria_type;
    5161         double nodeinputs[3];
    5162         double yts;
    5163         int    progstabilization,balancestabilization;
    5164         bool   dakota_analysis;
    5165 
    5166         /*Checks if debuging*/
    5167         /*{{{2*/
    5168         _assert_(iomodel->Data(MeshElementsEnum));
    5169         /*}}}*/
    5170 
    5171         /*Fetch parameters: */
    5172         iomodel->Constant(&yts,ConstantsYtsEnum);
    5173         iomodel->Constant(&progstabilization,PrognosticStabilizationEnum);
    5174         iomodel->Constant(&balancestabilization,BalancethicknessStabilizationEnum);
    5175         iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum);
    5176 
    5177         /*Recover element type*/
    5178         if ((analysis_type==PrognosticAnalysisEnum && progstabilization==3) || (analysis_type==BalancethicknessAnalysisEnum && balancestabilization==3)){
    5179                 /*P1 Discontinuous Galerkin*/
    5180                 tria_type=P1DGEnum;
    5181         }
    5182         else{
    5183                 /*P1 Continuous Galerkin*/
    5184                 tria_type=P1Enum;
    5185         }
    5186         this->SetElementType(tria_type,analysis_counter);
    5187 
    5188         /*Recover vertices ids needed to initialize inputs*/
    5189         for(i=0;i<3;i++){
    5190                 tria_vertex_ids[i]=(int)iomodel->Data(MeshElementsEnum)[3*index+i]; //ids for vertices are in the elements array from Matlab
    5191         }
    5192 
    5193         /*Recover nodes ids needed to initialize the node hook.*/
    5194         if (tria_type==P1DGEnum){
    5195                 /*Discontinuous Galerkin*/
    5196                 tria_node_ids[0]=iomodel->nodecounter+3*index+1;
    5197                 tria_node_ids[1]=iomodel->nodecounter+3*index+2;
    5198                 tria_node_ids[2]=iomodel->nodecounter+3*index+3;
    5199         }
    5200         else{
    5201                 /*Continuous Galerkin*/
    5202                 for(i=0;i<3;i++){
    5203                         tria_node_ids[i]=iomodel->nodecounter+(int)*(iomodel->Data(MeshElementsEnum)+3*index+i); //ids for vertices are in the elements array from Matlab
    5204                 }
    5205         }
    5206 
    5207         /*hooks: */
    5208         this->SetHookNodes(tria_node_ids,analysis_counter); this->nodes=NULL; //set hook to nodes, for this analysis type
    5209 
    5210         /*Fill with IoModel*/
    5211         this->InputUpdateFromIoModel(index,iomodel);
    5212 
    5213         /*Defaults if not provided in iomodel*/
    5214         switch(analysis_type){
    5215 
    5216                 case DiagnosticHorizAnalysisEnum:
    5217 
    5218                         /*default vx,vy and vz: either observation or 0 */
    5219                         if(!iomodel->Data(VxEnum)){
    5220                                 if (iomodel->Data(InversionVxObsEnum)) for(i=0;i<3;i++)nodeinputs[i]=iomodel->Data(InversionVxObsEnum)[tria_vertex_ids[i]-1]/yts;
    5221                                 else                 for(i=0;i<3;i++)nodeinputs[i]=0;
    5222                                 this->inputs->AddInput(new TriaVertexInput(VxEnum,nodeinputs));
    5223                                 this->inputs->AddInput(new TriaVertexInput(VxPicardEnum,nodeinputs));
    5224                                 if(dakota_analysis) this->inputs->AddInput(new TriaVertexInput(QmuVxEnum,nodeinputs));
    5225                         }
    5226                         if(!iomodel->Data(VyEnum)){
    5227                                 if (iomodel->Data(InversionVyObsEnum)) for(i=0;i<3;i++)nodeinputs[i]=iomodel->Data(InversionVyObsEnum)[tria_vertex_ids[i]-1]/yts;
    5228                                 else                 for(i=0;i<3;i++)nodeinputs[i]=0;
    5229                                 this->inputs->AddInput(new TriaVertexInput(VyEnum,nodeinputs));
    5230                                 this->inputs->AddInput(new TriaVertexInput(VyPicardEnum,nodeinputs));
    5231                                 if(dakota_analysis) this->inputs->AddInput(new TriaVertexInput(QmuVyEnum,nodeinputs));
    5232                         }
    5233                         if(!iomodel->Data(VzEnum)){
    5234                                 if (iomodel->Data(InversionVzObsEnum)) for(i=0;i<3;i++)nodeinputs[i]=iomodel->Data(InversionVzObsEnum)[tria_vertex_ids[i]-1]/yts;
    5235                                 else                 for(i=0;i<3;i++)nodeinputs[i]=0;
    5236                                 this->inputs->AddInput(new TriaVertexInput(VzEnum,nodeinputs));
    5237                                 this->inputs->AddInput(new TriaVertexInput(VzPicardEnum,nodeinputs));
    5238                                 if(dakota_analysis) this->inputs->AddInput(new TriaVertexInput(QmuVzEnum,nodeinputs));
    5239                         }
    5240                         if(!iomodel->Data(PressureEnum)){
    5241                                 for(i=0;i<3;i++)nodeinputs[i]=0;
    5242                                 if(dakota_analysis){
    5243                                         this->inputs->AddInput(new TriaVertexInput(PressureEnum,nodeinputs));
    5244                                         this->inputs->AddInput(new TriaVertexInput(QmuPressureEnum,nodeinputs));
    5245                                 }
    5246                         }
    5247                         break;
    5248 
    5249                 default:
    5250                         /*No update for other solution types*/
    5251                         break;
    5252 
    5253         }
    5254 
    5255         //this->parameters: we still can't point to it, it may not even exist. Configure will handle this.
    5256         this->parameters=NULL;
    5257 }
    5258 /*}}}*/
    5259 /*FUNCTION Tria::UpdateShelfStatus{{{1*/
    5260 int  Tria::UpdateShelfStatus(Vec new_shelf_nodes){
    5261 
    5262         int     i;
    5263 
    5264         double  h[3];
    5265         double  s[3];
    5266         double  b[3];
    5267         double  ba[3];
    5268         Input  *surface_input  = NULL;
    5269         Input  *bed_input      = NULL;
    5270         double  rho_water;
    5271         double  rho_ice;
    5272         double  density;
    5273         bool    elementonshelf = false;
    5274         int     flipped=0;
    5275         int     shelfstatus[3];
    5276 
    5277 
    5278         /*Recover info at the vertices: */
    5279         surface_input   =inputs->GetInput(SurfaceEnum);   _assert_(surface_input);
    5280         bed_input   =inputs->GetInput(BedEnum);   _assert_(bed_input);
    5281         if((surface_input->Enum()!=TriaVertexInputEnum) | (bed_input->Enum()!=TriaVertexInputEnum))_error_(" not supported yet for bed and surface interpolations not P1!");
    5282 
    5283         GetParameterListOnVertices(&h[0],ThicknessEnum);
    5284         GetParameterListOnVertices(&s[0],SurfaceEnum);
    5285         GetParameterListOnVertices(&b[0],BedEnum);
    5286         GetParameterListOnVertices(&ba[0],BathymetryEnum);
    5287 
    5288         /*material parameters: */
    5289         rho_water=matpar->GetRhoWater();
    5290         rho_ice=matpar->GetRhoIce();
    5291         density=rho_ice/rho_water;
    5292 
    5293         /*Initialize current status of nodes: */
    5294         for(i=0;i<3;i++){
    5295                 shelfstatus[i]=nodes[i]->IsOnShelf();
    5296                 if((nodes[i]->Sid()+1)==36) printf("UpdateShelfStatus: El %i Node %i shelf status %i\n",this->Id(),nodes[i]->Sid()+1,shelfstatus[i]);
    5297         }
    5298        
    5299         /*go through vertices, and figure out if they are grounded or not, then update their status: */
    5300         flipped=0;
    5301         for(i=0;i<3;i++){
    5302                 if(b[i]<=ba[i]){ //the = will lead to oscillations.
    5303                         nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,false));
    5304                         nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,true));
    5305                         if(shelfstatus[i]){
    5306                                 flipped++;
    5307                                 if((nodes[i]->Sid()+1)==36)printf("UpdateShelfStatus: El %i Node %i grounding %g %g %g\n",this->Id(),nodes[i]->Sid()+1,b[i],ba[i],s[i]);
    5308                         }
    5309                         VecSetValue(new_shelf_nodes,nodes[i]->Sid(),0,INSERT_VALUES);
    5310                 }
    5311                 else{
    5312                         nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,true));
    5313                         nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,false));
    5314                         if(!shelfstatus[i]){
    5315                                 flipped++;
    5316                                 if((nodes[i]->Sid()+1)==36)printf("UpdateShelfStatus: El %i Node %i floating %g %g %g\n",this->Id(),nodes[i]->Sid()+1,b[i],ba[i],s[i]);
    5317                         }
    5318                         VecSetValue(new_shelf_nodes,nodes[i]->Sid(),1,INSERT_VALUES);
    5319                 }
    5320         }
    5321 
    5322         /*Now, update  shelf status of element. An element can only be on shelf if all its nodes are on shelf: */
    5323         elementonshelf=false;
    5324         for(i=0;i<3;i++){
    5325                 if(nodes[i]->IsOnShelf()){
    5326                         elementonshelf=true;
    5327                         break;
    5328                 }
    5329         }
    5330     this->inputs->AddInput(new BoolInput(MaskElementonfloatingiceEnum,elementonshelf));
    5331 
    5332          return flipped;
    5333 }
    5334 /*}}}*/
    5335 /*FUNCTION Tria::UpdateShelfFlags{{{1*/
    5336 void Tria::UpdateShelfFlags(double* new_shelf_nodes){
    5337 
    5338         /*go through vertices, and update the status of MaskVertexonfloatingiceEnum and MaskVertexongroundediceEnum flags: */
    5339         bool flag;
    5340         int  i;
    5341         double  h[3];
    5342         double  s[3];
    5343         double  b[3];
    5344         double  ba[3];
    5345 
    5346         GetParameterListOnVertices(&h[0],ThicknessEnum);
    5347         GetParameterListOnVertices(&s[0],SurfaceEnum);
    5348         GetParameterListOnVertices(&b[0],BedEnum);
    5349         GetParameterListOnVertices(&ba[0],BathymetryEnum);
    5350 
    5351         for(i=0;i<3;i++){
    5352                 flag=(bool)new_shelf_nodes[nodes[i]->Sid()];
    5353                 nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,flag));
    5354                 nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,!flag));
    5355         }
    5356 }
    5357 /*}}}*/
    5358 /*FUNCTION Tria::UpdatePotentialSheetUngrounding{{{1*/
    5359 int Tria::UpdatePotentialSheetUngrounding(double* potential_sheet_ungrounding,Vec vec_nodes_on_iceshelf,double* nodes_on_iceshelf){
    5360 
    5361         /*intermediary: */
    5362         int i;
    5363         int nflipped=0;
    5364 
    5365         /*Ok, go through our 3 nodes, and whoever is on the potential_sheet_ungrounding, ends up in nodes_on_iceshelf: */
    5366         for(i=0;i<3;i++){
    5367                 if (potential_sheet_ungrounding[nodes[i]->Sid()]){
    5368                         VecSetValue(vec_nodes_on_iceshelf,nodes[i]->Sid(),1,INSERT_VALUES);
    5369                
    5370                         /*Figure out if we flipped: */
    5371                         if (potential_sheet_ungrounding[nodes[i]->Sid()] != nodes_on_iceshelf[nodes[i]->Sid()])nflipped++;
    5372                 }
    5373         }
    5374         return nflipped;
    5375 }
    5376 /*}}}*/
     5452        Input* surface_mass_balance_input=inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(surface_mass_balance_input);
     5453        Input* basal_melting_input=inputs->GetInput(BasalforcingsMeltingRateEnum);          _assert_(basal_melting_input);
     5454        Input* dhdt_input=inputs->GetInput(BalancethicknessThickeningRateEnum);                                       _assert_(dhdt_input);
     5455
     5456        /* Start  looping on the number of gaussian points: */
     5457        gauss=new GaussTria(2);
     5458        for(ig=gauss->begin();ig<gauss->end();ig++){
     5459
     5460                gauss->GaussPoint(ig);
     5461
     5462                surface_mass_balance_input->GetParameterValue(&surface_mass_balance_g,gauss);
     5463                basal_melting_input->GetParameterValue(&basal_melting_g,gauss);
     5464                dhdt_input->GetParameterValue(&dhdt_g,gauss);
     5465
     5466                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
     5467                GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
     5468
     5469                for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(surface_mass_balance_g-basal_melting_g-dhdt_g)*L[i];
     5470        }
     5471
     5472        /*Clean up and return*/
     5473        delete gauss;
     5474        return pe;
     5475}
     5476/*}}}*/
     5477#endif
     5478
Note: See TracChangeset for help on using the changeset viewer.