Changeset 8657


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

Fixed alphacomplement for 3d full-Stokes

Location:
issm/trunk/src/c/objects
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • TabularUnified issm/trunk/src/c/objects/Elements/Penta.cpp

    r8655 r8657  
    45964596
    45974597                /*Build alpha_complement_list: */
    4598                 if (drag_type==2) friction->GetAlphaComplement(&alpha_complement, gauss,VxEnum,VyEnum);
     4598                if (drag_type==2) friction->GetAlphaComplement(&alpha_complement, gauss,VxEnum,VyEnum,VzEnum);
    45994599                else alpha_complement=0;
    46004600
    … …  
    46624662        /*Build frictoin element, needed later: */
    46634663        inputs->GetParameterValue(&drag_type,DragTypeEnum);
    4664         friction=new Friction("2d",inputs,matpar,analysis_type);
     4664        friction=new Friction("3d",inputs,matpar,analysis_type);
    46654665
    46664666        /* Start  looping on the number of gaussian points: */
    … …  
    46714671
    46724672                /*Recover alpha_complement and drag: */
    4673                 if (drag_type==2) friction->GetAlphaComplement(&alpha_complement, gauss,VxEnum,VyEnum);
     4673                if (drag_type==2) friction->GetAlphaComplement(&alpha_complement, gauss,VxEnum,VyEnum,VzEnum);
    46744674                else alpha_complement=0;
    46754675                drag_input->GetParameterValue(&drag,gauss);

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

    r8654 r8657  
    30843084
    30853085                /*Build alpha_complement_list: */
    3086                 if (drag_type==2) friction->GetAlphaComplement(&alpha_complement, gauss,VxEnum,VyEnum);
     3086                if (drag_type==2) friction->GetAlphaComplement(&alpha_complement, gauss,VxEnum,VyEnum,VzEnum);
    30873087                else alpha_complement=0;
    30883088       

  • TabularUnified 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*/

  • TabularUnified issm/trunk/src/c/objects/Loads/Friction.h

    r8654 r8657  
    2929                void  GetAlpha2(double* palpha2, GaussTria* gauss,int vxenum,int vyenum,int vzenum);
    3030                void  GetAlpha2(double* palpha2, GaussPenta* gauss,int vxenum,int vyenum,int vzenum);
    31                 void  GetAlphaComplement(double* alpha_complement, GaussTria* gauss,int vxenum,int vyenum);
    32                 void  GetAlphaComplement(double* alpha_complement, GaussPenta* gauss,int vxenum,int vyenum);
     31                void  GetAlphaComplement(double* alpha_complement, GaussTria* gauss,int vxenum,int vyenum,int vzenum);
     32                void  GetAlphaComplement(double* alpha_complement, GaussPenta* gauss,int vxenum,int vyenum,int vzenum);
    3333                void  GetParameterValue(double* pvalue,GaussTria* gauss,int enum_type);
    3434                void  GetParameterValue(double* pvalue,GaussPenta* gauss,int enum_type);
Note: See TracChangeset for help on using the changeset viewer.