Changeset 17248


Ignore:
Timestamp:
02/10/14 14:38:22 (11 years ago)
Author:
Mathieu Morlighem
Message:

CHG: simplifying getviscosity calls, now only provide effective strain rate

Location:
issm/trunk-jpl
Files:
8 edited

Legend:

Unmodified
Added
Removed

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

    r17242 r17248  
    283283
    284284                this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input);
    285                 this->material->GetViscosity3dFS(&viscosity,&epsilon[0]);
     285                this->ViscosityFS(&viscosity,3,xyz_list,gauss,vx_input,vy_input,vz_input);
    286286                this->GetPhi(&phi,&epsilon[0],viscosity);
    287287
    … …  
    303303        IssmDouble epsilon3d[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
    304304        IssmDouble epsilon2d[3]; /* epsilon=[exx,eyy,exy];            */
     305        IssmDouble eps_eff;
     306        IssmDouble eps0=1.e-27;
    305307
    306308        if(dim==3){
    307                 /*3D*/
     309                /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */
    308310                this->StrainRateFS(&epsilon3d[0],xyz_list,gauss,vx_input,vy_input,vz_input);
    309                 material->GetViscosity3dFS(&viscosity, &epsilon3d[0]);
     311                eps_eff = sqrt(epsilon3d[0]*epsilon3d[0] + epsilon3d[1]*epsilon3d[1] + epsilon3d[3]*epsilon3d[3] +  epsilon3d[4]*epsilon3d[4] + epsilon3d[5]*epsilon3d[5] + epsilon3d[0]*epsilon3d[1]+eps0*eps0);
    310312        }
    311313        else{
    312                 /*2D*/
     314                /* eps_eff^2 = exx^2 + eyy^2 + 2*exy^2 */
    313315                this->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
    314                 material->GetViscosity2dvertical(&viscosity,&epsilon2d[0]);
    315         }
     316                eps_eff = sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + 2.*epsilon2d[2]*epsilon2d[2]);
     317        }
     318
     319        /*Get viscosity*/
     320        material->GetViscosity(&viscosity,eps_eff);
    316321
    317322        /*Assign output pointer*/
    … …  
    384389        IssmDouble epsilon3d[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/
    385390        IssmDouble epsilon2d[2]; /* epsilon=[exx,exy];           */
     391        IssmDouble eps_eff;
    386392
    387393        if(dim==3){
     394                /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */
    388395                this->StrainRateHO(&epsilon3d[0],xyz_list,gauss,vx_input,vy_input);
    389                 material->GetViscosity3d(&viscosity, &epsilon3d[0]);
     396                eps_eff = sqrt(epsilon3d[0]*epsilon3d[0] + epsilon3d[1]*epsilon3d[1] + epsilon3d[2]*epsilon3d[2] +  epsilon3d[3]*epsilon3d[3] + epsilon3d[4]*epsilon3d[4] + epsilon3d[0]*epsilon3d[1]);
    390397        }
    391398        else{
     399                /* eps_eff^2 = exx^2 + exy^2 */
    392400                this->StrainRateHO2dvertical(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
    393                 material->GetViscosity2dverticalHO(&viscosity, &epsilon2d[0]);
    394         }
     401                eps_eff = sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1]);
     402        }
     403
     404        /*Get viscosity*/
     405        material->GetViscosity(&viscosity,eps_eff);
    395406
    396407        /*Assign output pointer*/
    … …  
    402413        IssmDouble viscosity;
    403414        IssmDouble epsilon[3];/* epsilon=[exx,eyy,exy];    */
    404 
     415        IssmDouble eps_eff;
     416
     417        /*Get effective strain rate
     418         * eps_eff^2 = exx^2 + eyy^2 + exy^2 + exx*eyy+*/
    405419        this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
    406         material->GetViscosity2d(&viscosity, &epsilon[0]);
     420        eps_eff = sqrt(epsilon[0]*epsilon[0] + epsilon[1]*epsilon[1] + epsilon[2]*epsilon[2] + epsilon[0]*epsilon[1]);
     421
     422        /*Get viscosity*/
     423        material->GetViscosityBar(&viscosity,eps_eff);
    407424
    408425        /*Assign output pointer*/

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

    r17236 r17248  
    254254                /*Compute strain rate viscosity and pressure: */
    255255                this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    256                 material->GetViscosity3dFS(&viscosity,&epsilon[0]);
     256                this->ViscosityFS(&viscosity,3,&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    257257                pressure_input->GetInputValue(&pressure,gauss);
    258258
    … …  
    320320                /*Compute strain rate viscosity and pressure: */
    321321                this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    322                 material->GetViscosity3d(&viscosity,&epsilon[0]);
     322                this->ViscosityFS(&viscosity,3,&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    323323                pressure_input->GetInputValue(&pressure,gauss);
    324324
    … …  
    31383138        _assert_(gauss->Enum()==GaussPentaEnum);
    31393139        this->StrainRateFS(&epsilon[0],xyz_list,(GaussPenta*)gauss,vx_input,vy_input,vz_input);
    3140         material->GetViscosity3dFS(&viscosity,&epsilon[0]);
     3140        this->ViscosityFS(&viscosity,3,xyz_list,(GaussPenta*)gauss,vx_input,vy_input,vz_input);
    31413141        GetPhi(&phi,&epsilon[0],viscosity);
    31423142

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

    r17236 r17248  
    204204void  Tria::ComputeStressTensor(){
    205205
    206         IssmDouble      xyz_list[NUMVERTICES][3];
    207         IssmDouble      pressure,viscosity;
    208         IssmDouble      epsilon[3]; /* epsilon=[exx,eyy,exy];*/
    209         IssmDouble      sigma_xx[NUMVERTICES];
    210         IssmDouble              sigma_yy[NUMVERTICES];
    211         IssmDouble              sigma_zz[NUMVERTICES]={0,0,0};
    212         IssmDouble      sigma_xy[NUMVERTICES];
    213         IssmDouble              sigma_xz[NUMVERTICES]={0,0,0};
    214         IssmDouble              sigma_yz[NUMVERTICES]={0,0,0};
    215         GaussTria* gauss=NULL;
     206        IssmDouble  xyz_list[NUMVERTICES][3];
     207        IssmDouble  pressure,viscosity;
     208        IssmDouble  epsilon[3]; /* epsilon=[exx,eyy,exy];*/
     209        IssmDouble  sigma_xx[NUMVERTICES];
     210        IssmDouble      sigma_yy[NUMVERTICES];
     211        IssmDouble      sigma_zz[NUMVERTICES]={0,0,0};
     212        IssmDouble  sigma_xy[NUMVERTICES];
     213        IssmDouble      sigma_xz[NUMVERTICES]={0,0,0};
     214        IssmDouble      sigma_yz[NUMVERTICES]={0,0,0};
     215        GaussTria*  gauss=NULL;
    216216
    217217        /* Get node coordinates and dof list: */
    … …  
    230230                /*Compute strain rate viscosity and pressure: */
    231231                this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    232                 material->GetViscosity2d(&viscosity,&epsilon[0]);
     232                this->ViscositySSA(&viscosity,&xyz_list[0][0],gauss,vx_input,vy_input);
    233233                pressure_input->GetInputValue(&pressure,gauss);
    234234

  • issm/trunk-jpl/src/c/classes/Materials/Material.h

    r17242 r17248  
    2626                virtual void       Configure(Elements* elements)=0;
    2727                virtual void       GetVectorFromInputs(Vector<IssmDouble>* vector,int input_enum)=0;
    28                 virtual void       GetViscosity2d(IssmDouble* pviscosity, IssmDouble* pepsilon)=0;
    29                 virtual void       GetViscosity2dvertical(IssmDouble* pviscosity, IssmDouble* pepsilon)=0;
    30                 virtual void       GetViscosity2dverticalHO(IssmDouble* pviscosity, IssmDouble* pepsilon)=0;
    31                 virtual void       GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* pepsilon)=0;
    32                 virtual void       GetViscosity3dFS(IssmDouble* pviscosity3d, IssmDouble* epsilon)=0;
     28                virtual void       GetViscosity(IssmDouble* pviscosity,IssmDouble epseff)=0;
     29                virtual void       GetViscosityBar(IssmDouble* pviscosity,IssmDouble epseff)=0;
    3330                virtual void       GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon)=0;
    3431                virtual void       GetViscosityDComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon)=0;

  • issm/trunk-jpl/src/c/classes/Materials/Matice.cpp

    r17244 r17248  
    240240}
    241241/*}}}*/
    242 /*FUNCTION Matice::GetViscosity2d {{{*/
    243 void  Matice::GetViscosity2d(IssmDouble* pviscosity, IssmDouble* epsilon){
     242/*FUNCTION Matice::GetViscosity {{{*/
     243void  Matice::GetViscosity(IssmDouble* pviscosity,IssmDouble eps_eff){
    244244        /*From a string tensor and a material object, return viscosity, using Glen's flow law.
    245                                                                                                     (1-D) B
    246           viscosity= -------------------------------------------------------------------
    247                                                   2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n]
    248 
    249           where viscosity is the viscotiy, B the flow law parameter , (u,v) the velocity
    250           vector, and n the flow law exponent.
    251 
    252           If epsilon is NULL, it means this is the first time SystemMatrices is being run, and we
     245                                                                (1-D) B
     246          viscosity= -------------------------
     247                                                  2 eps_eff ^[(n-1)/n]
     248
     249          where viscosity is the viscotiy, B the flow law parameter , eps_eff is the effective strain rate
     250          and n the flow law exponent.
     251
     252          If eps_eff = 0 , it means this is the first time SystemMatrices is being run, and we
    253253          return 10^14, initial viscosity.
    254254          */
    … …  
    257257        IssmDouble viscosity;
    258258
    259         /*input strain rate: */
    260         IssmDouble exx,eyy,exy;
    261 
    262259        /*Intermediary: */
    263         IssmDouble A,e;
    264260        IssmDouble B,D,n;
    265261
    266262        /*Get B and n*/
    267         B=GetBbar();
    268         n=GetN();
    269         D=GetDbar();
    270 
    271         if (n==1){
    272                 /*Viscous behaviour! viscosity=B: */
    273                 viscosity=(1-D)*B/2;
     263        B=GetB(); _assert_(B>0.);
     264        n=GetN(); _assert_(n>0.);
     265        D=GetD(); _assert_(D>=0. && D<1.);
     266
     267        if (n==1.){
     268                /*Linear Viscous behavior (Newtonian fluid) viscosity=B/2: */
     269                viscosity=(1.-D)*B/2.;
    274270        }
    275271        else{
    276                 if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0)){
    277                         viscosity=0.5*pow(10.,14);
    278                 }
     272
     273                /*if no strain rate, return maximum viscosity*/
     274                if(eps_eff==0.){
     275                        viscosity = 1.e+14/2.;
     276                        //viscosity=2.5*pow(10.,17);
     277                }
     278
    279279                else{
    280                         /*Retrive strain rate components: */
    281                         exx=*(epsilon+0);
    282                         eyy=*(epsilon+1);
    283                         exy=*(epsilon+2);
    284 
    285                         /*Build viscosity: viscosity=B/(2*A^e) */
    286                         A=pow(exx,2)+pow(eyy,2)+pow(exy,2)+exx*eyy;
    287                         if(A==0){
    288                                 /*Maxiviscositym viscosity for 0 shear areas: */
    289                                 viscosity=2.5*pow(10.,17);
    290                         }
    291                         else{
    292                                 e=(n-1)/(2*n);
    293                                 viscosity=(1-D)*B/(2*pow(A,e));
    294                         }
     280                        viscosity=(1.-D)*B/(2.*pow(eps_eff,(n-1.)/n));
    295281                }
    296282        }
    … …  
    298284        /*Checks in debugging mode*/
    299285        if(viscosity<=0) _error_("Negative viscosity");
    300         _assert_(B>0);
    301         _assert_(n>0);
    302         _assert_(D>=0 && D<1);
    303286
    304287        /*Return: */
    … …  
    306289}
    307290/*}}}*/
    308 /*FUNCTION Matice::GetViscosity2dvertical {{{*/
    309 void  Matice::GetViscosity2dvertical(IssmDouble* pviscosity, IssmDouble* epsilon){
     291/*FUNCTION Matice::GetViscosityBar {{{*/
     292void  Matice::GetViscosityBar(IssmDouble* pviscosity,IssmDouble eps_eff){
    310293        /*From a string tensor and a material object, return viscosity, using Glen's flow law.
    311                                                                            (1-D) B
    312           viscosity= --------------------------------------
    313                                                   2[ exx^2+eyy^2+ 2exy^2]^[(n-1)/2n]
    314 
    315           where viscosity is the viscotiy, B the flow law parameter , (u,v) the velocity
    316           vector, and n the flow law exponent.
    317 
    318           If epsilon is NULL, it means this is the first time SystemMatrices is being run, and we
     294                                                                (1-D) B
     295          viscosity= -------------------------
     296                                                  2 eps_eff ^[(n-1)/n]
     297
     298          where viscosity is the viscotiy, B the flow law parameter , eps_eff is the effective strain rate
     299          and n the flow law exponent.
     300
     301          If eps_eff = 0 , it means this is the first time SystemMatrices is being run, and we
    319302          return 10^14, initial viscosity.
    320303          */
    … …  
    323306        IssmDouble viscosity;
    324307
    325         /*input strain rate: */
    326         IssmDouble exx,eyy,exy;
    327 
    328308        /*Intermediary: */
    329         IssmDouble A,e;
    330309        IssmDouble B,D,n;
    331310
    332311        /*Get B and n*/
    333         B=GetB();
    334         n=GetN();
    335         D=GetD();
    336 
    337         if (n==1){
    338                 /*Viscous behaviour! viscosity=B: */
    339                 viscosity=(1-D)*B/2;
     312        B=GetBbar(); _assert_(B>0.);
     313        n=GetN();    _assert_(n>0.);
     314        D=GetDbar(); _assert_(D>=0. && D<1.);
     315
     316        if (n==1.){
     317                /*Linear Viscous behavior (Newtonian fluid) viscosity=B/2: */
     318                viscosity=(1.-D)*B/2.;
    340319        }
    341320        else{
    342                 if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0)){
    343                         viscosity=0.5*pow(10.,14);
    344                 }
     321
     322                /*if no strain rate, return maximum viscosity*/
     323                if(eps_eff==0.){
     324                        viscosity = 1.e+14/2.;
     325                        //viscosity=2.5*pow(10.,17);
     326                }
     327
    345328                else{
    346                         /*Retrive strain rate components: */
    347                         exx=epsilon[0];
    348                         eyy=epsilon[1];
    349                         exy=epsilon[2];
    350 
    351                         /*Build viscosity: viscosity=B/(2*A^e) */
    352                         A=exx*exx+eyy*eyy+2.*exy*exy;
    353                         if(A==0.){
    354                                 /*Maxiviscositym viscosity for 0 shear areas: */
    355                                 viscosity=2.5*2.e+17;
    356                         }
    357                         else{
    358                                 e=(n-1.)/(2.*n);
    359                                 viscosity=(1.-D)*B/(2.*pow(A,e));
    360                         }
     329                        viscosity=(1.-D)*B/(2.*pow(eps_eff,(n-1.)/n));
    361330                }
    362331        }
    … …  
    364333        /*Checks in debugging mode*/
    365334        if(viscosity<=0) _error_("Negative viscosity");
    366         _assert_(B>0);
    367         _assert_(n>0);
    368         _assert_(D>=0 && D<1);
    369335
    370336        /*Return: */
    371337        *pviscosity=viscosity;
    372 }
    373 /*}}}*/
    374 /*FUNCTION Matice::GetViscosity2dverticalHO {{{*/
    375 void  Matice::GetViscosity2dverticalHO(IssmDouble* pviscosity, IssmDouble* epsilon){
    376         /*From a string tensor and a material object, return viscosity, using Glen's flow law.
    377                                                                            (1-D) B
    378           viscosity= --------------------------------------
    379                                                   2[ 2exx^2 + 2exy^2]^[(n-1)/2n]
    380 
    381           where viscosity is the viscotiy, B the flow law parameter , (u,v) the velocity
    382           vector, and n the flow law exponent.
    383 
    384           If epsilon is NULL, it means this is the first time SystemMatrices is being run, and we
    385           return 10^14, initial viscosity.
    386           */
    387 
    388         /*output: */
    389         IssmDouble viscosity;
    390 
    391         /*input strain rate: */
    392         IssmDouble exx,exy;
    393 
    394         /*Intermediary: */
    395         IssmDouble A,e;
    396         IssmDouble B,D,n;
    397 
    398         /*Get B and n*/
    399         B=GetB();
    400         n=GetN();
    401         D=GetD();
    402 
    403         if (n==1){
    404                 /*Viscous behaviour! viscosity=B: */
    405                 viscosity=(1-D)*B/2;
    406         }
    407         else{
    408                 if((epsilon[0]==0) && (epsilon[1]==0)){
    409                         viscosity=0.5*pow(10.,14);
    410                 }
    411                 else{
    412                         /*Retrive strain rate components: */
    413                         exx=epsilon[0];
    414                         exy=epsilon[1];
    415 
    416                         /*Build viscosity: viscosity=B/(2*A^e) */
    417                         A=2*exx*exx+2.*exy*exy;
    418                         if(A==0.){
    419                                 /*Maxiviscositym viscosity for 0 shear areas: */
    420                                 viscosity=2.5*2.e+17;
    421                         }
    422                         else{
    423                                 e=(n-1.)/(2.*n);
    424                                 viscosity=(1.-D)*B/(2.*pow(A,e));
    425                         }
    426                 }
    427         }
    428 
    429         /*Checks in debugging mode*/
    430         if(viscosity<=0) _error_("Negative viscosity");
    431         _assert_(B>0);
    432         _assert_(n>0);
    433         _assert_(D>=0 && D<1);
    434 
    435         /*Return: */
    436         *pviscosity=viscosity;
    437 }
    438 /*}}}*/
    439 /*FUNCTION Matice::GetViscosity3d {{{*/
    440 void  Matice::GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* epsilon){
    441 
    442         /*Return viscosity accounting for steady state power law creep [Thomas and SSA, 1982]:
    443          *
    444          *                                               (1-D)*B
    445          * viscosity3d= -------------------------------------------------------------------
    446          *                      2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n]
    447          *
    448          *     where mu is the viscotiy, B the flow law parameter , (u,v) the velocity
    449          *     vector, and n the flow law exponent.
    450          *
    451          * If epsilon is NULL, it means this is the first time Emg is being run, and we
    452          * return g, initial viscosity.
    453          */
    454 
    455         /*output: */
    456         IssmDouble viscosity3d;
    457 
    458         /*input strain rate: */
    459         IssmDouble exx,eyy,exy,exz,eyz;
    460 
    461         /*Intermediaries: */
    462         IssmDouble A,e;
    463         IssmDouble B,D,n;
    464 
    465         /*Get B and n*/
    466         B=GetB();
    467         D=GetD();
    468         n=GetN();
    469 
    470         if (n==1){
    471                 /*Viscous behaviour! viscosity3d=B: */
    472                 viscosity3d=(1-D)*B/2;
    473         }
    474         else{
    475                 if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) &&
    476                                 (epsilon[3]==0) && (epsilon[4]==0)){
    477                         viscosity3d=0.5*pow(10.,14);
    478                 }
    479                 else{
    480 
    481                         /*Retrive strain rate components: */
    482                         exx=*(epsilon+0);
    483                         eyy=*(epsilon+1);
    484                         exy=*(epsilon+2);
    485                         exz=*(epsilon+3);
    486                         eyz=*(epsilon+4);
    487 
    488                         /*Build viscosity: viscosity3d=2*B/(2*A^e) */
    489                         A=pow(exx,2)+pow(eyy,2)+pow(exy,2)+pow(exz,2)+pow(eyz,2)+exx*eyy;
    490                         if(A==0){
    491                                 /*Maxiviscosity3dm viscosity for 0 shear areas: */
    492                                 viscosity3d=2.25*pow(10.,17);
    493                         }
    494                         else{
    495                                 e=(n-1)/2/n;
    496 
    497                                 viscosity3d=(1-D)*B/(2*pow(A,e));
    498                         }
    499                 }
    500         }
    501 
    502         /*Checks in debugging mode*/
    503         if(viscosity3d<=0) _error_("Negative viscosity");
    504         _assert_(B>0);
    505         _assert_(n>0);
    506         _assert_(D>=0 && D<1);
    507 
    508         /*Assign output pointers:*/
    509         *pviscosity3d=viscosity3d;
    510 }
    511 /*}}}*/
    512 /*FUNCTION Matice::GetViscosity3dFS {{{*/
    513 void  Matice::GetViscosity3dFS(IssmDouble* pviscosity3d, IssmDouble* epsilon){
    514         /*Return viscosity accounting for steady state power law creep [Thomas and SSA, 1982]:
    515          *
    516          *                                          (1-D)*B
    517          * viscosity3d= -------------------------------------------------------------------
    518          *                   2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n]
    519          *
    520          *     where mu is the viscotiy, B the flow law parameter , (u,v) the velocity
    521          *     vector, and n the flow law exponent.
    522          *
    523          * If epsilon is NULL, it means this is the first time Emg is being run, and we
    524          * return g, initial viscosity.
    525          */
    526 
    527         /*output: */
    528         IssmDouble viscosity3d;
    529 
    530         /*input strain rate: */
    531         IssmDouble exx,eyy,exy,exz,eyz,ezz;
    532 
    533         /*Intermediaries: */
    534         IssmDouble A,e;
    535         IssmDouble B,D,n;
    536         IssmDouble eps0;
    537 
    538         /*Get B and n*/
    539         eps0=pow(10.,-27);
    540         B=GetB();
    541         D=GetD();
    542         n=GetN();
    543 
    544         if (n==1){
    545                 /*Viscous behaviour! viscosity3d=B: */
    546                 viscosity3d=(1-D)*B/2;
    547         }
    548         else{
    549                 if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) &&
    550                                 (epsilon[3]==0) && (epsilon[4]==0) && (epsilon[5]==0)){
    551                         viscosity3d=0.5*pow(10.,14);
    552                 }
    553                 else{
    554 
    555                         /*Retrive strain rate components: */
    556                         exx=*(epsilon+0);
    557                         eyy=*(epsilon+1);
    558                         ezz=*(epsilon+2); //not used
    559                         exy=*(epsilon+3);
    560                         exz=*(epsilon+4);
    561                         eyz=*(epsilon+5);
    562 
    563                         /*Build viscosity: viscosity3d=B/(2*A^e) */
    564                         A=pow(exx,2)+pow(eyy,2)+pow(exy,2)+pow(exz,2)+pow(eyz,2)+exx*eyy+pow(eps0,2);
    565                         if(A==0){
    566                                 /*Maxiviscosity3dm viscosity for 0 shear areas: */
    567                                 viscosity3d=2.25*pow(10.,17);
    568                         }
    569                         else{
    570                                 e=(n-1)/2/n;
    571                                 viscosity3d=(1-D)*B/(2*pow(A,e));
    572                         }
    573                 }
    574         }
    575 
    576         /*Checks in debugging mode*/
    577         if(viscosity3d<=0) _error_("Negative viscosity");
    578         _assert_(B>0);
    579         _assert_(n>0);
    580         _assert_(D>=0 && D<1);
    581 
    582         /*Assign output pointers:*/
    583         *pviscosity3d=viscosity3d;
    584338}
    585339/*}}}*/
    … …  
    706460        IssmDouble exx,eyy,exy,exz,eyz;
    707461
    708         /*Get visocisty and n*/
    709         GetViscosity3d(&mu,epsilon);
    710         n=GetN();
    711462
    712463        if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) &&
    … …  
    715466        }
    716467        else{
     468
    717469                /*Retrive strain rate components: */
    718470                exx=epsilon[0];
    … …  
    723475                eff2 = exx*exx + eyy*eyy + exx*eyy + exy*exy + exz*exz + eyz*eyz;
    724476
    725                 mu_prime=(1-n)/(2*n) * mu/eff2;
     477                GetViscosity(&mu,sqrt(eff2));
     478                n=GetN();
     479                mu_prime=(1.-n)/(2.*n) * mu/eff2;
    726480        }
    727481
    … …  
    739493        /*input strain rate: */
    740494        IssmDouble exx,eyy,exy,exz,eyz,ezz;
    741 
    742         /*Get visocisty and n*/
    743         GetViscosity3d(&mu,epsilon);
    744         n=GetN();
    745495
    746496        if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) &&
    … …  
    758508                eff2 = exx*exx + eyy*eyy + exx*eyy + exy*exy + exz*exz + eyz*eyz;
    759509
     510                GetViscosity(&mu,sqrt(eff2));
     511                n=GetN();
    760512                mu_prime=(1-n)/(2*n) * mu/eff2;
    761513        }
    … …  
    774526        /*input strain rate: */
    775527        IssmDouble exx,eyy,exy;
    776 
    777         /*Get visocisty and n*/
    778         GetViscosity2d(&mu,epsilon);
    779         n=GetN();
    780528
    781529        if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0)){
    … …  
    789537                eff2 = exx*exx + eyy*eyy + exx*eyy + exy*exy ;
    790538
    791                 mu_prime=(1-n)/(2*n) * mu/eff2;
     539                GetViscosityBar(&mu,sqrt(eff2));
     540                n=GetN();
     541                mu_prime=(1.-n)/(2.*n)*mu/eff2;
    792542        }
    793543

  • issm/trunk-jpl/src/c/classes/Materials/Matice.h

    r17242 r17248  
    5252                void   GetVectorFromInputs(Vector<IssmDouble>* vector,int input_enum);
    5353                void       SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
    54                 void       GetViscosity2d(IssmDouble* pviscosity, IssmDouble* pepsilon);
    55                 void       GetViscosity2dvertical(IssmDouble* pviscosity, IssmDouble* pepsilon);
    56                 void       GetViscosity2dverticalHO(IssmDouble* pviscosity, IssmDouble* pepsilon);
    57                 void       GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* pepsilon);
    58                 void       GetViscosity3dFS(IssmDouble* pviscosity3d, IssmDouble* epsilon);
     54                void       GetViscosity(IssmDouble* pviscosity, IssmDouble eps_eff);
     55                void       GetViscosityBar(IssmDouble* pviscosity, IssmDouble eps_eff);
    5956                void       GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon);
    6057                void       GetViscosityDComplement(IssmDouble*, IssmDouble*);

  • issm/trunk-jpl/src/c/classes/Materials/Matpar.h

    r17242 r17248  
    8282                void       Configure(Elements* elements);
    8383                void       GetVectorFromInputs(Vector<IssmDouble>* vector,int input_enum){return;}
    84                 void       GetViscosity2d(IssmDouble* pviscosity, IssmDouble* pepsilon){_error_("not supported");};
    85                 void       GetViscosity2dvertical(IssmDouble* pviscosity, IssmDouble* pepsilon){_error_("not supported");};
    86                 void       GetViscosity2dverticalHO(IssmDouble* pviscosity, IssmDouble* pepsilon){_error_("not supported");};
    87                 void       GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* pepsilon){_error_("not supported");};
    88                 void       GetViscosity3dFS(IssmDouble* pviscosity3d, IssmDouble* epsilon){_error_("not supported");};
     84                void       GetViscosity(IssmDouble* pviscosity,IssmDouble eps_eff){_error_("not supported");};
     85                void       GetViscosityBar(IssmDouble* pviscosity,IssmDouble eps_eff){_error_("not supported");};
    8986                void       GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon){_error_("not supported");};
    9087                void       GetViscosityDComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon){_error_("not supported");};
Note: See TracChangeset for help on using the changeset viewer.