Ignore:
Timestamp:
06/17/11 13:20:05 (14 years ago)
Author:
Mathieu Morlighem
Message:

Fixed alphacomplement for 3d full-Stokes

File:
1 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/objects/Loads/Friction.cpp

    r8654 r8657  
    7171        double  alpha2;
    7272
    73 
    7473        /*Recover parameters: */
    7574        inputs->GetParameterValue(&drag_p,DragPEnum);
    … …  
    136135        double  alpha2;
    137136
    138 
    139137        /*Recover parameters: */
    140138        inputs->GetParameterValue(&drag_p,DragPEnum);
    … …  
    185183}
    186184/*}}}*/
    187 /*FUNCTION Friction::GetAlphaComplement(double* palpha_complement, GaussTria* gauss,int vxenum,int vyenum) {{{1*/
    188 void Friction::GetAlphaComplement(double* palpha_complement, GaussTria* gauss,int vxenum,int vyenum){
     185/*FUNCTION Friction::GetAlphaComplement(double* palpha_complement, GaussTria* gauss,int vxenum,int vyenum,int vzenum) {{{1*/
     186void Friction::GetAlphaComplement(double* palpha_complement, GaussTria* gauss,int vxenum,int vyenum,int vzenum){
    189187
    190188        /* 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.
    … …  
    194192        /*diverse: */
    195193        int     i;
     194        double  r,s;
     195        double  vx,vy,vz,vmag;
     196        double  drag_p,drag_q;
    196197        double  Neff;
    197         double  r,s;
    198         double  vx;
    199         double  vy;
    200         double  vmag;
    201         double  drag_p,drag_q;
    202198        double  drag_coefficient;
    203199        double  bed,thickness;
    … …  
    217213        rho_water=matpar->GetRhoWater();
    218214
    219 
    220215        //compute r and q coefficients: */
    221216        r=drag_q/drag_p;
    … …  
    232227
    233228        //We need the velocity magnitude to evaluate the basal stress:
    234         this->GetParameterValue(&vx, gauss,vxenum);
    235         this->GetParameterValue(&vy, gauss,vyenum);
    236         vmag=sqrt(pow(vx,2)+pow(vy,2));
     229        if(strcmp(element_type,"2d")==0){
     230                this->GetParameterValue(&vx, gauss,vxenum);
     231                this->GetParameterValue(&vy, gauss,vyenum);
     232                vmag=sqrt(pow(vx,2)+pow(vy,2));
     233        }
     234        else if (strcmp(element_type,"3d")==0){
     235                this->GetParameterValue(&vx, gauss,vxenum);
     236                this->GetParameterValue(&vy, gauss,vyenum);
     237                this->GetParameterValue(&vz, gauss,vzenum);
     238                vmag=sqrt(pow(vx,2)+pow(vy,2)+pow(vz,2));
     239        }
     240        else _error_("element_type %s not supported yet",element_type);
    237241
    238242        /*Checks that s-1>0 if v=0*/
    … …  
    245249}
    246250/*}}}*/
    247 /*FUNCTION Friction::GetAlphaComplement(double* palpha_complement, GaussPenta* gauss,int vxenum,int vyenum) {{{1*/
    248 void Friction::GetAlphaComplement(double* palpha_complement, GaussPenta* gauss,int vxenum,int vyenum){
     251/*FUNCTION Friction::GetAlphaComplement(double* palpha_complement, GaussPenta* gauss,int vxenum,int vyenum,int vzenum) {{{1*/
     252void Friction::GetAlphaComplement(double* palpha_complement, GaussPenta* gauss,int vxenum,int vyenum,int vzenum){
    249253
    250254        /* 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.
    … …  
    254258        /*diverse: */
    255259        int     i;
     260        double  r,s;
     261        double  vx,vy,vz,vmag;
     262        double  drag_p,drag_q;
    256263        double  Neff;
    257         double  r,s;
    258         double  vx;
    259         double  vy;
    260         double  vmag;
    261         double  drag_p,drag_q;
    262264        double  drag_coefficient;
    263265        double  bed,thickness;
    … …  
    277279        rho_water=matpar->GetRhoWater();
    278280
    279 
    280281        //compute r and q coefficients: */
    281282        r=drag_q/drag_p;
    … …  
    292293
    293294        //We need the velocity magnitude to evaluate the basal stress:
    294         this->GetParameterValue(&vx, gauss,vxenum);
    295         this->GetParameterValue(&vy, gauss,vyenum);
    296         vmag=sqrt(pow(vx,2)+pow(vy,2));
     295        if(strcmp(element_type,"2d")==0){
     296                this->GetParameterValue(&vx, gauss,vxenum);
     297                this->GetParameterValue(&vy, gauss,vyenum);
     298                vmag=sqrt(pow(vx,2)+pow(vy,2));
     299        }
     300        else if (strcmp(element_type,"3d")==0){
     301                this->GetParameterValue(&vx, gauss,vxenum);
     302                this->GetParameterValue(&vy, gauss,vyenum);
     303                this->GetParameterValue(&vz, gauss,vzenum);
     304                vmag=sqrt(pow(vx,2)+pow(vy,2)+pow(vz,2));
     305        }
     306        else _error_("element_type %s not supported yet",element_type);
    297307
    298308        /*Checks that s-1>0 if v=0*/
Note: See TracChangeset for help on using the changeset viewer.