Changeset 16873


Ignore:
Timestamp:
11/21/13 18:55:31 (11 years ago)
Author:
seroussi
Message:

CHG: changed friction's interface

Location:
issm/trunk-jpl/src/c/classes
Files:
8 edited

Legend:

Unmodified
Added
Removed

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

    r16872 r16873  
    118118                virtual void   GetInputValue(int* pvalue,int enum_type)=0;
    119119                virtual void   GetInputValue(IssmDouble* pvalue,int enum_type)=0;
     120                virtual void   GetInputValue(IssmDouble* pvalue,Gauss* gauss,int enum_type)=0;
    120121                virtual void   GetVerticesCoordinates(IssmDouble** xyz_list)=0;
    121122                virtual void   GetVerticesCoordinatesBase(IssmDouble** xyz_list)=0;

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

    r16864 r16873  
    193193
    194194        /*Build friction element, needed later: */
    195         friction=new Friction("3d",inputs,matpar,StressbalanceAnalysisEnum);
     195        friction=new Friction(this,3);
    196196
    197197        /* Start looping on the number of gauss 2d (nodes on the bedrock) */
    … …  
    202202                gauss->GaussPoint(ig);
    203203
    204                 friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
     204                friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,vz_input);
    205205                vx_input->GetInputValue(&vx,gauss);
    206206                vy_input->GetInputValue(&vy,gauss);
    … …  
    14451445        if(!input) _error_("Input " << EnumToStringx(inputenum) << " not found in element");
    14461446        input->GetInputValue(pvalue);
     1447
     1448}/*}}}*/
     1449/*FUNCTION Penta::GetInputValue(IssmDouble* pvalue,Gauss* gauss,int inputenum) {{{*/
     1450void Penta::GetInputValue(IssmDouble* pvalue,Gauss* gauss,int inputenum){
     1451
     1452        Input* input=inputs->GetInput(inputenum);
     1453        if(!input) _error_("Input " << EnumToStringx(inputenum) << " not found in element");
     1454        input->GetInputValue(pvalue,gauss);
    14471455
    14481456}/*}}}*/
    … …  
    45804588
    45814589        /*Build friction element, needed later: */
    4582         friction=new Friction("3d",inputs,matpar,analysis_type);
     4590        friction=new Friction(this,3);
    45834591
    45844592        /* Start looping on the number of gauss 2d (nodes on the bedrock) */
    … …  
    46034611                        geothermalflux_input->GetInputValue(&geothermalflux_value,gauss);
    46044612
    4605                         friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
     4613                        friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,vz_input);
    46064614                        vx_input->GetInputValue(&vx,gauss);
    46074615                        vy_input->GetInputValue(&vy,gauss);
    … …  
    48394847
    48404848        /*Build frictoin element, needed later: */
    4841         friction=new Friction("3d",inputs,matpar,analysis_type);
     4849        friction=new Friction(this,3);
    48424850
    48434851        /* Start looping on the number of gauss 2d (nodes on the bedrock) */
    … …  
    48514859
    48524860                        geothermalflux_input->GetInputValue(&geothermalflux_value,gauss);
    4853                         friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
     4861                        friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,vz_input);
    48544862                        vx_input->GetInputValue(&vx,gauss);
    48554863                        vy_input->GetInputValue(&vy,gauss);
    … …  
    49854993        GetInputListOnVertices(&geothermalflux[0],BasalforcingsGeothermalfluxEnum);
    49864994        Input* enthalpy_input=inputs->GetInput(EnthalpyEnum); _assert_(enthalpy_input);
     4995        Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
     4996        Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
     4997        Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input);
    49874998
    49884999        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    … …  
    49915002        /*Build friction element, needed later: */
    49925003        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    4993         friction=new Friction("3d",inputs,matpar,analysis_type);
     5004        friction=new Friction(this,3);
    49945005
    49955006        /******** MELTING RATES  ************************************/
    … …  
    50325043
    50335044                        /*basal friction*/
    5034                         friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
     5045                        friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,vz_input);
    50355046                        basalfriction=alpha2*(pow(vx[iv],2.0)+pow(vy[iv],2.0)+pow(vz[iv],2.0));
    50365047
    … …  
    60286039
    60296040        /*Build frictoin element, needed later: */
    6030         friction=new Friction("2d",inputs,matpar,analysis_type);
     6041        friction=new Friction(this,2);
    60316042
    60326043        /* Start  looping on the number of gaussian points: */
    … …  
    60406051
    60416052                /*Build alpha_complement_list: */
    6042                 friction->GetAlphaComplement(&alpha_complement, gauss,VxEnum,VyEnum,VzEnum);
     6053                friction->GetAlphaComplement(&alpha_complement,gauss,vx_input,vy_input,NULL);
    60436054
    60446055                dragcoefficient_input->GetInputValue(&drag, gauss);
    … …  
    61036114
    61046115        /*Build frictoin element, needed later: */
    6105         friction=new Friction("3d",inputs,matpar,analysis_type);
     6116        friction=new Friction(this,3);
    61066117
    61076118        /* Start  looping on the number of gaussian points: */
    … …  
    61126123
    61136124                /*Recover alpha_complement and drag: */
    6114                 friction->GetAlphaComplement(&alpha_complement, gauss,VxEnum,VyEnum,VzEnum);
     6125                friction->GetAlphaComplement(&alpha_complement,gauss,vx_input,vy_input,vz_input);
    61156126                drag_input->GetInputValue(&drag,gauss);
    61166127
    … …  
    69056916
    69066917        /*build friction object, used later on: */
    6907         friction=new Friction("2d",inputs,matpar,analysis_type);
     6918        friction=new Friction(this,2);
    69086919
    69096920        /* Start  looping on the number of gaussian points: */
    … …  
    69146925
    69156926                /*Friction: */
    6916                 friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum);
     6927                friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,vz_input);
    69176928
    69186929                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss);
    … …  
    71397150
    71407151        /*build friction object, used later on: */
    7141         friction=new Friction("3d",inputs,matpar,analysis_type);
     7152        friction=new Friction(this,3);
    71427153
    71437154        /* Start  looping on the number of gaussian points: */
    … …  
    71577168
    71587169                NormalBase(&bed_normal[0],&xyz_list_tria[0][0]);
    7159                 friction->GetAlpha2(&alpha2_gauss, gauss,VxEnum,VyEnum,VzEnum);
     7170                friction->GetAlpha2(&alpha2_gauss, gauss,vx_input,vy_input,vz_input);
    71607171
    71617172                DLSSAFS[0][0]=alpha2_gauss*gauss->weight*Jdet2d;
    … …  
    77627773
    77637774        /*build friction object, used later on: */
    7764         friction=new Friction("2d",inputs,matpar,analysis_type);
     7775        friction=new Friction(this,2);
    77657776
    77667777        /*Recover portion of element that is grounded*/
    … …  
    77847795                GetBHOFriction(&B[0],gauss);
    77857796
    7786                 friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum);
     7797                friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,vz_input);
    77877798                if(migration_style==SubelementMigrationEnum) alpha2=phi*alpha2;
    77887799                if(migration_style==SubelementMigration2Enum){
    … …  
    80668077
    80678078        /*build friction object, used later on: */
    8068         friction=new Friction("3d",inputs,matpar,analysis_type);
     8079        friction=new Friction(this,3);
    80698080
    80708081        /* Start  looping on the number of gaussian points: */
    … …  
    80778088                GetLFS(BFriction,gauss);
    80788089
    8079                 friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
     8090                friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,vz_input);
    80808091
    80818092                D[0*2+0] = +alpha2*gauss->weight*Jdet2d; //taub_x = -alpha2 vx
    … …  
    83248335
    83258336        /*build friction object, used later on: */
    8326         friction=new Friction("3d",inputs,matpar,analysis_type);
     8337        friction=new Friction(this,3);
    83278338
    83288339        /* Start looping on the number of gauss 2d (nodes on the bedrock) */
    … …  
    83418352                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    83428353                material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    8343                 friction->GetAlpha2(&alpha2_gauss, gauss,VxEnum,VyEnum,VzEnum);
     8354                friction->GetAlpha2(&alpha2_gauss, gauss,vx_input,vy_input,vz_input);
    83448355
    83458356                for(i=0;i<NUMVERTICES2D;i++){
    … …  
    84898500
    84908501        /*build friction object, used later on: */
    8491         friction=new Friction("3d",inputs,matpar,analysis_type);
     8502        friction=new Friction(this,3);
    84928503
    84938504        /* Start looping on the number of gauss 2d (nodes on the bedrock) */
    … …  
    85068517                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    85078518                material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    8508                 friction->GetAlpha2(&alpha2_gauss, gauss,VxEnum,VyEnum,VzEnum);
     8519                friction->GetAlpha2(&alpha2_gauss, gauss,vx_input,vy_input,vz_input);
    85098520
    85108521                for(i=0;i<NUMVERTICES2D;i++){

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

    r16872 r16873  
    241241                void           GetInputValue(int* pvalue,int enum_type);
    242242                void           GetInputValue(IssmDouble* pvalue,int enum_type);
     243                void           GetInputValue(IssmDouble* pvalue,Gauss* gauss,int enum_type);
    243244                void           GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype);
    244245                void             GetPhi(IssmDouble* phi, IssmDouble*  epsilon, IssmDouble viscosity);

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

    r16872 r16873  
    152152                void        GetInputValue(int* pvalue,int enum_type){_error_("not implemented yet");};
    153153                void        GetInputValue(IssmDouble* pvalue,int enum_type){_error_("not implemented yet");};
     154                void        GetInputValue(IssmDouble* pvalue,Gauss* gauss,int enum_type){_error_("not implemented yet");};
    154155                void        GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype){_error_("not implemented yet");};
    155156                IssmDouble  GetYcoord(Gauss* gauss){_error_("Not implemented");};

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

    r16872 r16873  
    15261526
    15271527}/*}}}*/
     1528/*FUNCTION Tria::GetInputValue(IssmDouble* pvalue,Gauss* gauss,int inputenum) {{{*/
     1529void Tria::GetInputValue(IssmDouble* pvalue,Gauss* gauss,int inputenum){
     1530
     1531        Input* input=inputs->GetInput(inputenum);
     1532        if(!input) _error_("Input " << EnumToStringx(inputenum) << " not found in element");
     1533        input->GetInputValue(pvalue,gauss);
     1534
     1535}/*}}}*/
    15281536/*FUNCTION Tria::GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype) {{{*/
    15291537void Tria::GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype){
    … …  
    38053813
    38063814        /*build friction object, used later on: */
    3807         friction=new Friction("1d",inputs,matpar,analysis_type);
     3815        friction=new Friction(this,1);
    38083816
    38093817        /* Start looping on the number of gauss 1d (nodes on the bedrock) */
    … …  
    38163824                GetLFS(BFriction,gauss);
    38173825
    3818                 friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
     3826                friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,NULL);
    38193827                D = +alpha2*gauss->weight*Jdet; //taub_x = -alpha2 vx
    38203828
    … …  
    39543962
    39553963        /*build friction object, used later on: */
    3956         friction=new Friction("2d",inputs,matpar,analysis_type);
     3964        friction=new Friction(this,2);
    39573965
    39583966        /*Recover portion of element that is grounded*/
    … …  
    39723980                gauss->GaussPoint(ig);
    39733981
    3974                 friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum);
     3982                friction->GetAlpha2(&alpha2, gauss,vx_input,vy_input,NULL);
    39753983                if(migration_style==SubelementMigrationEnum) alpha2=phi*alpha2;
    39763984                if(migration_style==SubelementMigration2Enum){
    … …  
    49614969
    49624970        /*Build frictoin element, needed later: */
    4963         friction=new Friction("2d",inputs,matpar,analysis_type);
     4971        friction=new Friction(this,2);
    49644972
    49654973        /*Retrieve all inputs we will be needing: */
    … …  
    49804988
    49814989                /*Build alpha_complement_list: */
    4982                 friction->GetAlphaComplement(&alpha_complement, gauss,VxEnum,VyEnum,VzEnum);
     4990                friction->GetAlphaComplement(&alpha_complement,gauss,vx_input,vy_input,NULL);
    49834991
    49844992                dragcoefficient_input->GetInputValue(&drag, gauss);
    … …  
    50055013        //for (int iv=0;iv<NUMVERTICES;iv++){
    50065014        //      gauss->GaussVertex(iv);
    5007         //      friction->GetAlphaComplement(&alpha_complement, gauss,VxEnum,VyEnum,VzEnum);
     5015        //      friction->GetAlphaComplement(&alpha_complement,gauss,vx_input,vy_input,vz_input);
    50085016        //      dragcoefficient_input->GetInputValue(&drag, gauss);
    50095017        //      adjointx_input->GetInputValue(&lambda, gauss);

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

    r16872 r16873  
    283283                void           GetInputValue(int* pvalue,int enum_type);
    284284                void           GetInputValue(IssmDouble* pvalue,int enum_type);
     285                void           GetInputValue(IssmDouble* pvalue,Gauss* gauss,int enum_type);
    285286                void           GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype);
    286287                void           GetStrainRate2d(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss, Input* vx_input, Input* vy_input);

  • issm/trunk-jpl/src/c/classes/Loads/Friction.cpp

    r16373 r16873  
    1818/*FUNCTION Friction::Friction() {{{*/
    1919Friction::Friction(){
    20         this->element_type=NULL;
    21         this->inputs=NULL;
    22         this->matpar=NULL;
     20        this->element=NULL;
     21        this->dim=0;
    2322}
    2423/*}}}*/
    25 /*FUNCTION Friction::Friction(const char* element_type, Inputs* inputs,Matpar* matpar,int analysis_type){{{*/
    26 Friction::Friction(const char* element_type_in,Inputs* inputs_in,Matpar* matpar_in, int in_analysis_type){
     24/*FUNCTION Friction::Friction(Element* element,int dim){{{*/
     25Friction::Friction(Element* element_in,int dim_in){
    2726
    28         this->analysis_type=in_analysis_type;
    29         this->inputs=inputs_in;
    30         this->element_type=xNew<char>(strlen(element_type_in)+1);
    31         xMemCpy<char>(this->element_type,element_type_in,(strlen(element_type_in)+1));
    32 
    33         this->matpar=matpar_in;
     27        this->element=element_in;
     28        this->dim=dim_in;
    3429}
    3530/*}}}*/
    3631/*FUNCTION Friction::~Friction() {{{*/
    3732Friction::~Friction(){
    38         xDelete<char>(element_type);
    3933}
    4034/*}}}*/
    … …  
    4438void Friction::Echo(void){
    4539        _printf_("Friction:\n");
    46         _printf_("   analysis_type: " << EnumToStringx(analysis_type) << "\n");
    47         _printf_("   element_type: " << this->element_type << "\n");
    48         inputs->Echo();
    49         matpar->Echo();
     40        _printf_("   dim: " << this->dim<< "\n");
    5041}
    5142/*}}}*/
    52 /*FUNCTION Friction::GetAlpha2(IssmDouble* palpha2, GaussTria* gauss,int vxenum,int vyenum,int vzenum){{{*/
    53 void Friction::GetAlpha2(IssmDouble* palpha2, GaussTria* gauss,int vxenum,int vyenum,int vzenum){
     43/*FUNCTION Friction::GetAlpha2{{{*/
     44void Friction::GetAlpha2(IssmDouble* palpha2, Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){
    5445
    5546        /*This routine calculates the basal friction coefficient
    … …  
    5950        IssmDouble  r,s;
    6051        IssmDouble  drag_p, drag_q;
    61         IssmDouble  gravity,rho_ice,rho_water;
    6252        IssmDouble  Neff;
    6353        IssmDouble  thickness,bed;
    … …  
    6757
    6858        /*Recover parameters: */
    69         inputs->GetInputValue(&drag_p,FrictionPEnum);
    70         inputs->GetInputValue(&drag_q,FrictionQEnum);
    71         this->GetInputValue(&thickness, gauss,ThicknessEnum);
    72         this->GetInputValue(&bed, gauss,BedEnum);
    73         this->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
    74 
    75         /*Get material parameters: */
    76         gravity=matpar->GetG();
    77         rho_ice=matpar->GetRhoIce();
    78         rho_water=matpar->GetRhoWater();
     59        element->GetInputValue(&drag_p,FrictionPEnum);
     60        element->GetInputValue(&drag_q,FrictionQEnum);
     61        element->GetInputValue(&thickness, gauss,ThicknessEnum);
     62        element->GetInputValue(&bed, gauss,BedEnum);
     63        element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
     64        IssmDouble rho_water   = element->GetMaterialParameter(MaterialsRhoWaterEnum);
     65        IssmDouble rho_ice     = element->GetMaterialParameter(MaterialsRhoIceEnum);
     66        IssmDouble gravity     = element->GetMaterialParameter(ConstantsGEnum);
    7967
    8068        //compute r and q coefficients: */
    … …  
    8472        //From bed and thickness, compute effective pressure when drag is viscous:
    8573        Neff=gravity*(rho_ice*thickness+rho_water*bed);
     74        if(Neff<0)Neff=0;
    8675
    87         /*If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because
    88           the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour
    89           for friction should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should
    90           replace it by Neff=0 (ie, equival it to an ice shelf)*/
    91         if (Neff<0)Neff=0;
    92 
    93         if(strcmp(element_type,"1d")==0){
    94                 this->GetInputValue(&vx, gauss,vxenum);
    95                 vmag=sqrt(vx*vx);
     76        switch(dim){
     77                case 1:
     78                        vx_input->GetInputValue(&vx,gauss);
     79                        vmag=sqrt(vx*vx);
     80                        break;
     81                case 2:
     82                        vx_input->GetInputValue(&vx,gauss);
     83                        vy_input->GetInputValue(&vy,gauss);
     84                        vmag=sqrt(vx*vx+vy*vy);
     85                        break;
     86                case 3:
     87                        vx_input->GetInputValue(&vx,gauss);
     88                        vy_input->GetInputValue(&vy,gauss);
     89                        vz_input->GetInputValue(&vz,gauss);
     90                        vmag=sqrt(vx*vx+vy*vy+vz*vz);
     91                        break;
     92                default:
     93                        _error_("not supported");
    9694        }
    97         else if(strcmp(element_type,"2d")==0){
    98                 this->GetInputValue(&vx, gauss,vxenum);
    99                 this->GetInputValue(&vy, gauss,vyenum);
    100                 vmag=sqrt(pow(vx,2)+pow(vy,2));
    101         }
    102         else if (strcmp(element_type,"3d")==0){
    103                 this->GetInputValue(&vx, gauss,vxenum);
    104                 this->GetInputValue(&vy, gauss,vyenum);
    105                 this->GetInputValue(&vz, gauss,vzenum);
    106                 vmag=sqrt(pow(vx,2)+pow(vy,2)+pow(vz,2));
    107         }
    108         else _error_("element_type "<< element_type << " not supported yet");
    10995
    11096        /*Checks that s-1>0 if v=0*/
    111         if(vmag==0 && (s-1)<0) _error_("velocity is 0 and (s-1)=" << (s-1) << "<0, alpha2 is Inf");
     97        if(vmag==0. && (s-1.)<0.) _error_("velocity is 0 and (s-1)=" << (s-1.) << "<0, alpha2 is Inf");
    11298
    113         alpha2=pow(drag_coefficient,2)*pow(Neff,r)*pow(vmag,(s-1));
     99        alpha2=drag_coefficient*drag_coefficient*pow(Neff,r)*pow(vmag,(s-1.));
    114100        _assert_(!xIsNan<IssmDouble>(alpha2));
    115101
    … …  
    118104}
    119105/*}}}*/
    120 /*FUNCTION Friction::GetAlpha2(IssmDouble* palpha2, GaussPenta* gauss,int vxenum,int vyenum,int vzenum){{{*/
    121 void Friction::GetAlpha2(IssmDouble* palpha2, GaussPenta* gauss,int vxenum,int vyenum,int vzenum){
    122 
    123         /*This routine calculates the basal friction coefficient
    124           alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p**/
    125 
    126         /*diverse: */
    127         IssmDouble  r,s;
    128         IssmDouble  drag_p, drag_q;
    129         IssmDouble  gravity,rho_ice,rho_water;
    130         IssmDouble  Neff;
    131         IssmDouble  thickness,bed;
    132         IssmDouble  vx,vy,vz,vmag;
    133         IssmDouble  drag_coefficient;
    134         IssmDouble  alpha2;
    135 
    136         /*Recover parameters: */
    137         inputs->GetInputValue(&drag_p,FrictionPEnum);
    138         inputs->GetInputValue(&drag_q,FrictionQEnum);
    139         this->GetInputValue(&thickness, gauss,ThicknessEnum);
    140         this->GetInputValue(&bed, gauss,BedEnum);
    141         this->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
    142 
    143         /*Get material parameters: */
    144         gravity=matpar->GetG();
    145         rho_ice=matpar->GetRhoIce();
    146         rho_water=matpar->GetRhoWater();
    147 
    148         //compute r and q coefficients: */
    149         r=drag_q/drag_p;
    150         s=1./drag_p;
    151 
    152         //From bed and thickness, compute effective pressure when drag is viscous:
    153         Neff=gravity*(rho_ice*thickness+rho_water*bed);
    154 
    155         /*If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because
    156           the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour
    157           for friction should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should
    158           replace it by Neff=0 (ie, equival it to an ice shelf)*/
    159         if (Neff<0)Neff=0;
    160 
    161         if(strcmp(element_type,"2d")==0){
    162                 this->GetInputValue(&vx, gauss,vxenum);
    163                 this->GetInputValue(&vy, gauss,vyenum);
    164                 vmag=sqrt(pow(vx,2)+pow(vy,2));
    165         }
    166         else if (strcmp(element_type,"3d")==0){
    167                 this->GetInputValue(&vx, gauss,vxenum);
    168                 this->GetInputValue(&vy, gauss,vyenum);
    169                 this->GetInputValue(&vz, gauss,vzenum);
    170                 vmag=sqrt(pow(vx,2)+pow(vy,2)+pow(vz,2));
    171         }
    172         else _error_("element_type "<< element_type << " not supported yet");
    173 
    174         /*Checks that s-1>0 if v=0*/
    175         if(vmag==0 && (s-1)<0) _error_("velocity is 0 and (s-1)=" << (s-1) << "<0, alpha_complement is Inf");
    176 
    177         alpha2=pow(drag_coefficient,2)*pow(Neff,r)*pow(vmag,(s-1));
    178         _assert_(!xIsNan<IssmDouble>(alpha2));
    179 
    180         /*Assign output pointers:*/
    181         *palpha2=alpha2;
    182 }
    183 /*}}}*/
    184 /*FUNCTION Friction::GetAlphaComplement(IssmDouble* palpha_complement, GaussTria* gauss,int vxenum,int vyenum,int vzenum) {{{*/
    185 void Friction::GetAlphaComplement(IssmDouble* palpha_complement, GaussTria* gauss,int vxenum,int vyenum,int vzenum){
     106/*FUNCTION Friction::GetAlphaComplement(IssmDouble* palpha_complement, Gauss* gauss,int vxenum,int vyenum,int vzenum) {{{*/
     107void Friction::GetAlphaComplement(IssmDouble* palpha_complement, Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){
    186108
    187109        /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p.
    … …  
    196118        IssmDouble  drag_coefficient;
    197119        IssmDouble  bed,thickness;
    198         IssmDouble  gravity,rho_ice,rho_water;
    199120        IssmDouble  alpha_complement;
    200121
    201122        /*Recover parameters: */
    202         inputs->GetInputValue(&drag_p,FrictionPEnum);
    203         inputs->GetInputValue(&drag_q,FrictionQEnum);
    204         this->GetInputValue(&thickness, gauss,ThicknessEnum);
    205         this->GetInputValue(&bed, gauss,BedEnum);
    206         this->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
    207 
    208         /*Get material parameters: */
    209         gravity=matpar->GetG();
    210         rho_ice=matpar->GetRhoIce();
    211         rho_water=matpar->GetRhoWater();
     123        element->GetInputValue(&drag_p,FrictionPEnum);
     124        element->GetInputValue(&drag_q,FrictionQEnum);
     125        element->GetInputValue(&thickness, gauss,ThicknessEnum);
     126        element->GetInputValue(&bed, gauss,BedEnum);
     127        element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
     128        IssmDouble rho_water   = element->GetMaterialParameter(MaterialsRhoWaterEnum);
     129        IssmDouble rho_ice     = element->GetMaterialParameter(MaterialsRhoIceEnum);
     130        IssmDouble gravity     = element->GetMaterialParameter(ConstantsGEnum);
    212131
    213132        //compute r and q coefficients: */
    … …  
    217136        //From bed and thickness, compute effective pressure when drag is viscous:
    218137        Neff=gravity*(rho_ice*thickness+rho_water*bed);
    219 
    220         /*If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because
    221           the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour
    222           for friction should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should
    223           replace it by Neff=0 (ie, equival it to an ice shelf)*/
    224         if (Neff<0)Neff=0;
     138        if(Neff<0)Neff=0;
    225139
    226140        //We need the velocity magnitude to evaluate the basal stress:
    227         if(strcmp(element_type,"2d")==0){
    228                 this->GetInputValue(&vx, gauss,vxenum);
    229                 this->GetInputValue(&vy, gauss,vyenum);
    230                 vmag=sqrt(pow(vx,2)+pow(vy,2));
     141        switch(dim){
     142                case 1:
     143                        vx_input->GetInputValue(&vx,gauss);
     144                        vmag=sqrt(vx*vx);
     145                        break;
     146                case 2:
     147                        vx_input->GetInputValue(&vx,gauss);
     148                        vy_input->GetInputValue(&vy,gauss);
     149                        vmag=sqrt(vx*vx+vy*vy);
     150                        break;
     151                case 3:
     152                        vx_input->GetInputValue(&vx,gauss);
     153                        vy_input->GetInputValue(&vy,gauss);
     154                        vz_input->GetInputValue(&vz,gauss);
     155                        vmag=sqrt(vx*vx+vy*vy+vz*vz);
     156                        break;
     157                default:
     158                        _error_("not supported");
    231159        }
    232         else if (strcmp(element_type,"3d")==0){
    233                 this->GetInputValue(&vx, gauss,vxenum);
    234                 this->GetInputValue(&vy, gauss,vyenum);
    235                 this->GetInputValue(&vz, gauss,vzenum);
    236                 vmag=sqrt(pow(vx,2)+pow(vy,2)+pow(vz,2));
    237         }
    238         else _error_("element_type "<< element_type << " not supported yet");
    239160
    240161        /*Checks that s-1>0 if v=0*/
    241         if(vmag==0 && (s-1)<0) _error_("velocity is 0 and (s-1)=" << (s-1) << "<0, alpha_complement is Inf");
     162        if(vmag==0. && (s-1.)<0.) _error_("velocity is 0 and (s-1)=" << (s-1.) << "<0, alpha2 is Inf");
    242163
    243         alpha_complement=pow(Neff,r)*pow(vmag,(s-1));            _assert_(!xIsNan<IssmDouble>(alpha_complement));
     164        alpha_complement=pow(Neff,r)*pow(vmag,(s-1));_assert_(!xIsNan<IssmDouble>(alpha_complement));
    244165
    245166        /*Assign output pointers:*/
    … …  
    247168}
    248169/*}}}*/
    249 /*FUNCTION Friction::GetAlphaComplement(IssmDouble* palpha_complement, GaussPenta* gauss,int vxenum,int vyenum,int vzenum) {{{*/
    250 void Friction::GetAlphaComplement(IssmDouble* palpha_complement, GaussPenta* gauss,int vxenum,int vyenum,int vzenum){
    251 
    252         /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p.
    253          * FrictionGetAlphaComplement is used in control methods on drag, and it computes:
    254          * alpha_complement= Neff ^r * vel ^s*/
    255 
    256         /*diverse: */
    257         IssmDouble  r,s;
    258         IssmDouble  vx,vy,vz,vmag;
    259         IssmDouble  drag_p,drag_q;
    260         IssmDouble  Neff;
    261         IssmDouble  drag_coefficient;
    262         IssmDouble  bed,thickness;
    263         IssmDouble  gravity,rho_ice,rho_water;
    264         IssmDouble  alpha_complement;
    265 
    266         /*Recover parameters: */
    267         inputs->GetInputValue(&drag_p,FrictionPEnum);
    268         inputs->GetInputValue(&drag_q,FrictionQEnum);
    269         this->GetInputValue(&thickness, gauss,ThicknessEnum);
    270         this->GetInputValue(&bed, gauss,BedEnum);
    271         this->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
    272 
    273         /*Get material parameters: */
    274         gravity=matpar->GetG();
    275         rho_ice=matpar->GetRhoIce();
    276         rho_water=matpar->GetRhoWater();
    277 
    278         //compute r and q coefficients: */
    279         r=drag_q/drag_p;
    280         s=1./drag_p;
    281 
    282         //From bed and thickness, compute effective pressure when drag is viscous:
    283         Neff=gravity*(rho_ice*thickness+rho_water*bed);
    284 
    285         /*If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because
    286           the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour
    287           for friction should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should
    288           replace it by Neff=0 (ie, equival it to an ice shelf)*/
    289         if (Neff<0)Neff=0;
    290 
    291         //We need the velocity magnitude to evaluate the basal stress:
    292         if(strcmp(element_type,"2d")==0){
    293                 this->GetInputValue(&vx, gauss,vxenum);
    294                 this->GetInputValue(&vy, gauss,vyenum);
    295                 vmag=sqrt(pow(vx,2)+pow(vy,2));
    296         }
    297         else if (strcmp(element_type,"3d")==0){
    298                 this->GetInputValue(&vx, gauss,vxenum);
    299                 this->GetInputValue(&vy, gauss,vyenum);
    300                 this->GetInputValue(&vz, gauss,vzenum);
    301                 vmag=sqrt(pow(vx,2)+pow(vy,2)+pow(vz,2));
    302         }
    303         else _error_("element_type "<< element_type << " not supported yet");
    304 
    305         /*Checks that s-1>0 if v=0*/
    306         if(vmag==0 && (s-1)<0) _error_("velocity is 0 and (s-1)=" << (s-1) << "<0, alpha_complement is Inf");
    307 
    308         alpha_complement=pow(Neff,r)*pow(vmag,(s-1));            _assert_(!xIsNan<IssmDouble>(alpha_complement));
    309 
    310         /*Assign output pointers:*/
    311         *palpha_complement=alpha_complement;
    312 }
    313 /*}}}*/
    314 /*FUNCTION Friction::GetInputValue{{{*/
    315 void Friction::GetInputValue(IssmDouble* pvalue,GaussTria* gauss,int enum_type){
    316 
    317         Input* input=inputs->GetInput(enum_type);
    318         if(!input) _error_("input " << EnumToStringx(enum_type) << " not found");
    319         input->GetInputValue(pvalue,gauss);
    320 
    321 }
    322 /*}}}*/
    323 /*FUNCTION Friction::GetInputValue{{{*/
    324 void Friction::GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,int enum_type){
    325 
    326         Input* input=inputs->GetInput(enum_type);
    327         if(!input) _error_("input " << EnumToStringx(enum_type) << " not found");
    328         input->GetInputValue(pvalue,gauss);
    329 
    330 }
    331 /*}}}*/

  • issm/trunk-jpl/src/c/classes/Loads/Friction.h

    r14951 r16873  
    1919                int analysis_type;
    2020
    21                 char* element_type;
    22                 Inputs* inputs;
    23                 Matpar* matpar;
     21                Element* element;
     22                int      dim;
    2423
    2524                /*methods: */
    2625                Friction();
    27                 Friction(const char* element_type, Inputs* inputs,Matpar* matpar, int analysis_type);
     26                Friction(Element* element_in,int dim_in);
    2827                ~Friction();
    2928
    3029                void  Echo(void);
    31                 void  GetAlpha2(IssmDouble* palpha2, GaussTria* gauss,int vxenum,int vyenum,int vzenum);
    32                 void  GetAlpha2(IssmDouble* palpha2, GaussPenta* gauss,int vxenum,int vyenum,int vzenum);
    33                 void  GetAlphaComplement(IssmDouble* alpha_complement, GaussTria* gauss,int vxenum,int vyenum,int vzenum);
    34                 void  GetAlphaComplement(IssmDouble* alpha_complement, GaussPenta* gauss,int vxenum,int vyenum,int vzenum);
    35                 void  GetInputValue(IssmDouble* pvalue,GaussTria* gauss,int enum_type);
    36                 void  GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,int enum_type);
    37 
     30                void  GetAlpha2(IssmDouble* palpha2, Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input);
     31                void  GetAlphaComplement(IssmDouble* alpha_complement,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input);
    3832};
    3933
Note: See TracChangeset for help on using the changeset viewer.