Changeset 5682


Ignore:
Timestamp:
09/07/10 08:04:36 (15 years ago)
Author:
seroussi
Message:

fixed friction stokes

Location:
issm/trunk/src/c
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk/src/c/Container/Inputs.cpp

    r5660 r5682  
    9999}
    100100/*}}}*/
     101/*FUNCTION Inputs::GetParameterValue(double* pvalue,GaussTria* gauss,int enum_type){{{1*/
     102void Inputs::GetParameterValue(double* pvalue,GaussPenta* gauss, int enum_type){
     103
     104        vector<Object*>::iterator object;
     105        Input* input=NULL;
     106        bool   found=false;
     107
     108        /*Go through inputs and check whether any input with the same name is already in: */
     109        for ( object=objects.begin() ; object < objects.end(); object++ ){
     110
     111                input=(Input*)(*object);
     112                if (input->EnumType()==enum_type){
     113                        found=true;
     114                        break;
     115                }
     116        }
     117
     118        if (!found){
     119                /*we could not find an input with the correct enum type. No defaults values were provided,
     120                 * error out: */
     121                ISSMERROR("could not find input with enum type %i (%s)",enum_type,EnumToString(enum_type));
     122        }
     123
     124        /*Ok, we have an input if we made it here, request the input to return the values: */
     125        input->GetParameterValue(pvalue,gauss);
     126
     127}
     128/*}}}*/
    101129/*FUNCTION Inputs::GetParameterValue(double* pvalue,double* gauss,int enum_type,double defaultvalue){{{1*/
    102130void Inputs::GetParameterValue(double* pvalue,double* gauss, int enum_type,double defaultvalue){
  • issm/trunk/src/c/Container/Inputs.h

    r5660 r5682  
    1717class Node;
    1818class GaussTria;
     19class GaussPenta;
    1920
    2021class Inputs: public DataSet{
     
    4950                void GetParameterValue(double* pvalue,double* gauss,int enum_type);
    5051                void GetParameterValue(double* pvalue,GaussTria* gauss,int enum_type);
     52                void GetParameterValue(double* pvalue,GaussPenta* gauss,int enum_type);
    5153                void GetParameterValue(double* pvalue,double* gauss,int enum_type,double defaultvalue);
    5254                void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss,int enum_type);
  • issm/trunk/src/c/objects/Elements/Penta.cpp

    r5678 r5682  
    28792879
    28802880                        /*Calculate DL on gauss point */
    2881                         friction->GetAlpha2(&alpha2_gauss, gauss_tria,VxEnum,VyEnum,VzEnum);
     2881                        friction->GetAlpha2(&alpha2_gauss, gauss,VxEnum,VyEnum,VzEnum);
    28822882
    28832883                        DLStokes[0][0]=alpha2_gauss*gauss->weight*Jdet2d;
  • issm/trunk/src/c/objects/Loads/Friction.cpp

    r5640 r5682  
    177177}
    178178/*}}}*/
     179/*FUNCTION Friction::GetAlpha2{{{1*/
     180void Friction::GetAlpha2(double* palpha2, GaussPenta* gauss,int vxenum,int vyenum,int vzenum){
     181
     182        /*This routine calculates the basal friction coefficient
     183          alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p**/
     184
     185        /*diverse: */
     186        double  r,s;
     187        double  drag_p, drag_q;
     188        double  gravity,rho_ice,rho_water;
     189        double  Neff;
     190        double  thickness,bed;
     191        double  vx,vy,vz,vmag;
     192        double  drag_coefficient;
     193        double  alpha2;
     194
     195
     196        /*Recover parameters: */
     197        inputs->GetParameterValue(&drag_p,DragPEnum);
     198        inputs->GetParameterValue(&drag_q,DragQEnum);
     199        inputs->GetParameterValue(&thickness, gauss,ThicknessEnum);
     200        inputs->GetParameterValue(&bed, gauss,BedEnum);
     201        inputs->GetParameterValue(&drag_coefficient, gauss,DragCoefficientEnum);
     202
     203        /*Get material parameters: */
     204        gravity=matpar->GetG();
     205        rho_ice=matpar->GetRhoIce();
     206        rho_water=matpar->GetRhoWater();
     207
     208        //compute r and q coefficients: */
     209        r=drag_q/drag_p;
     210        s=1./drag_p;
     211
     212        //From bed and thickness, compute effective pressure when drag is viscous:
     213        Neff=gravity*(rho_ice*thickness+rho_water*bed);
     214
     215        /*If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because
     216          the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour
     217          for friction should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should
     218          replace it by Neff=0 (ie, equival it to an ice shelf)*/
     219        if (Neff<0)Neff=0;
     220
     221        if(strcmp(element_type,"2d")==0){
     222                inputs->GetParameterValue(&vx, gauss,vxenum);
     223                inputs->GetParameterValue(&vy, gauss,vyenum);
     224                vmag=sqrt(pow(vx,2)+pow(vy,2));
     225        }
     226        else if (strcmp(element_type,"3d")==0){
     227                inputs->GetParameterValue(&vx, gauss,vxenum);
     228                inputs->GetParameterValue(&vy, gauss,vyenum);
     229                inputs->GetParameterValue(&vz, gauss,vzenum);
     230                vmag=sqrt(pow(vx,2)+pow(vy,2)+pow(vz,2));
     231        }
     232        else ISSMERROR("element_type %s not supported yet",element_type);
     233
     234        alpha2=pow(drag_coefficient,2)*pow(Neff,r)*pow(vmag,(s-1));
     235
     236        /*Assign output pointers:*/
     237        *palpha2=alpha2;
     238}
     239/*}}}*/
    179240/*FUNCTION Friction::GetAlphaComplement {{{1*/
    180241void Friction::GetAlphaComplement(double* palpha_complement, double* gauss,int vxenum,int vyenum){
  • issm/trunk/src/c/objects/Loads/Friction.h

    r5640 r5682  
    2929                void  GetAlpha2(double* palpha2, double* gauss,int vxenum,int vyenum,int vzenum);
    3030                void  GetAlpha2(double* palpha2, GaussTria* gauss,int vxenum,int vyenum,int vzenum);
     31                void  GetAlpha2(double* palpha2, GaussPenta* gauss,int vxenum,int vyenum,int vzenum);
    3132                void  GetAlphaComplement(double* alpha_complement, double* gauss,int vxenum,int vyenum);
    3233                void  GetAlphaComplement(double* alpha_complement, GaussTria* gauss,int vxenum,int vyenum);
Note: See TracChangeset for help on using the changeset viewer.