Ignore:
Timestamp:
12/22/21 10:39:44 (3 years ago)
Author:
Mathieu Morlighem
Message:

merged trunk-jpl and trunk for revision 26742

Location:
issm/trunk
Files:
12 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk

  • issm/trunk/src

  • issm/trunk/src/c

    • Property svn:ignore

      •  

        old new  
        2020issm
        2121kriging
         22issm_slc
        2223issm_slr
        2324issm_ocean

  • issm/trunk/src/c/classes/Loads/Channel.cpp

    r25836 r26744  
    320320                        switch(set2_enum){
    321321                                case FsetEnum:
    322                                         if(nodes[i]->fsize){
     322                                        if(nodes[i]->FSize()){
    323323                                                if(this->nodes[i]->IsClone())
    324324                                                 o_nz += 1;
    … …  
    336336                                        break;
    337337                                case SsetEnum:
    338                                         if(nodes[i]->ssize){
     338                                        if(nodes[i]->SSize()){
    339339                                                if(this->nodes[i]->IsClone())
    340340                                                 o_nz += 1;

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

    r25836 r26744  
    1818Friction::Friction(){/*{{{*/
    1919        this->element=NULL;
    20         this->dim=0;
    2120        this->law=0;
    22 
     21        this->apply_dim = 1.;
     22        this->domaintype=-1;
     23        this->vx_input=NULL;
     24        this->vy_input=NULL;
     25        this->vz_input=NULL;
    2326}
    2427/*}}}*/
    25 Friction::Friction(Element* element_in,int dim_in){/*{{{*/
    26 
     28Friction::Friction(Element* element_in){/*{{{*/
     29        /* Determine the dimension according to the domain type automatically.
     30         * There are exceptions, e.g. HO, which needs the user to specify the dimension used in Friciton.*/
    2731        this->element=element_in;
    28         this->dim=dim_in;
     32
     33        /* Load necessary parameters */
    2934        element_in->FindParam(&this->law,FrictionLawEnum);
     35        element_in->FindParam(&this->domaintype,DomainTypeEnum);
     36
     37        /* Load VxBase and VyBase for this special case */
     38        switch(this->domaintype){
     39                case Domain2DhorizontalEnum:
     40                        this->apply_dim = 2.;
     41                        this->vx_input = element_in->GetInput(VxBaseEnum);                      _assert_(this->vx_input);
     42                        this->vy_input = element_in->GetInput(VyBaseEnum);                      _assert_(this->vy_input);
     43                        this->vz_input = NULL;
     44                        break;
     45      case Domain2DverticalEnum:
     46                        this->apply_dim = 2.;
     47                        this->vx_input = element_in->GetInput(VxEnum);                          _assert_(this->vx_input);
     48                        this->vy_input = element_in->GetInput(VyEnum);                          _assert_(this->vy_input);
     49                        this->vz_input = NULL;
     50                        break;
     51      case Domain3DEnum:           
     52                        this->apply_dim = 3.;
     53                        this->vx_input = element_in->GetInput(VxEnum);                          _assert_(this->vx_input);
     54                        this->vy_input = element_in->GetInput(VyEnum);                          _assert_(this->vy_input);
     55                        this->vz_input = element_in->GetInput(VzEnum);                          _assert_(this->vz_input);
     56                        break;
     57      default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
     58        }
     59}
     60/*}}}*/
     61Friction::Friction(Element* element_in,int dim) : Friction(element_in) {/*{{{*/
     62        this->apply_dim = reCast<IssmDouble>(dim);
     63}
     64/*}}}*/
     65Friction::Friction(Element* element_in,IssmDouble dim) : Friction(element_in) {/*{{{*/
     66        this->apply_dim = dim;
    3067}
    3168/*}}}*/
    … …  
    3471/*}}}*/
    3572
     73
    3674/*methods: */
    3775void Friction::Echo(void){/*{{{*/
    3876        _printf_("Friction:\n");
    39         _printf_("   dim: " << this->dim<< "\n");
     77        _printf_("   Domain type: " << this->domaintype<< "\n");
    4078}
    4179/*}}}*/
    … …  
    6199                        _error_("not supported");
    62100        }
    63 
    64101
    65102        /*Checks*/
    … …  
    160197        IssmDouble vmag = VelMag(gauss);
    161198
    162         /*Check to prevent dividing by zero if vmag==0*/
    163         if(vmag==0. && (s-1.)<0.) alpha_complement=0.;
    164         else alpha_complement=pow(Neff,r)*pow(vmag,(s-1));
     199        if(s==1.){
     200                /*This is to make AD happy and avoid 0^0*/
     201                alpha_complement=pow(Neff,r);
     202        }
     203        else{
     204                /*Check to prevent dividing by zero if vmag==0*/
     205                if(vmag==0. && (s-1.)<0.) alpha_complement=0.;
     206                else alpha_complement=pow(Neff,r)*pow(vmag,(s-1.));
     207        }
    165208
    166209        /*Assign output pointers:*/
    … …  
    219262        if(vmag==0. && (1./m-1.)<0.) alpha_complement=0.;
    220263        else alpha_complement= pow(vmag, 1.0/m-1.);
    221        
     264
    222265        /*Assign output pointers:*/
    223266        *palpha_complement=alpha_complement;
    … …  
    261304                case 12:
    262305                        GetAlpha2Tsai(palpha2,gauss);
     306                        break;
     307                case 13:
     308                        GetAlpha2Coulomb2(palpha2,gauss);
    263309                        break;
    264310          default:
    … …  
    304350        IssmDouble vmag = VelMag(gauss);
    305351
    306         /*Check to prevent dividing by zero if vmag==0*/
    307         if(vmag==0. && (s-1.)<0.){
    308                 alpha2=0.;
     352        if(s==1.){
     353                /*This is to make AD happy and avoid 0^0*/
     354                alpha2=drag_coefficient*drag_coefficient*pow(Neff,r);
    309355        }
    310356        else{
    311                 alpha2=drag_coefficient*drag_coefficient*pow(Neff,r)*pow(vmag,(s-1.));
     357                /*Check to prevent dividing by zero if vmag==0*/
     358                if(vmag==0. && (s-1.)<0.) alpha2=0.;
     359                else alpha2=drag_coefficient*drag_coefficient*pow(Neff,r)*pow(vmag,(s-1.));
    312360        }
    313361
    … …  
    461509        deltaT = T-Tpmp;
    462510
    463 
    464511        /*Compute gamma*/
    465512        element->parameters->FindParam(&time,TimeEnum);
    … …  
    498545
    499546        /*Get effective pressure and basal velocity*/
    500         IssmDouble Neff = EffectivePressure(gauss);
    501547        IssmDouble vmag = VelMag(gauss);
    502548
     549        bool ispwStochastic;
     550        IssmDouble Neff;
     551        element->parameters->FindParam(&ispwStochastic,StochasticForcingIsWaterPressureEnum);
     552        if(ispwStochastic){
     553                /*Retrieve stochastic water pressure and compute ice pressure*/
     554                IssmDouble p_ice,p_water,Neff_limit;
     555                element->GetInputValue(&p_water,gauss,FrictionWaterPressureEnum);
     556                element->parameters->FindParam(&Neff_limit,FrictionEffectivePressureLimitEnum);
     557                p_ice = IcePressure(gauss);
     558                Neff  = max(Neff_limit*p_ice, p_ice - p_water);
     559        }       
     560        else{
     561                /*Compute effective pressure directly*/
     562                Neff = EffectivePressure(gauss);
     563        }
     564
    503565        /*Check to prevent dividing by zero if vmag==0*/
    504         if(vmag==0. && (s-1.)<0.) alpha2=0.;
    505         else alpha2=drag_coefficient*drag_coefficient*pow(Neff,r)*pow(vmag,(s-1.));
     566        if(s==1.){
     567                /*This is to make AD happy and avoid 0^0*/
     568                alpha2=drag_coefficient*drag_coefficient*pow(Neff,r);
     569        }
     570        else{
     571                if(vmag==0. && (s-1.)<0.) alpha2=0.;
     572                else alpha2=drag_coefficient*drag_coefficient*pow(Neff,r)*pow(vmag,(s-1.));
     573        }
    506574
    507575        /*Assign output pointers:*/
    … …  
    547615        IssmDouble vmag = VelMag(gauss);
    548616
    549         /*Check to prevent dividing by zero if vmag==0*/
    550         if(vmag==0. && (s-1.)<0.) alpha2=0.;
    551         else alpha2=drag_coefficient*drag_coefficient*pow(Neff,r)*pow(vmag,(s-1.));
     617        if(s==1.){
     618                /*This is to make AD happy and avoid 0^0*/
     619                alpha2=drag_coefficient*drag_coefficient*pow(Neff,r);
     620        }
     621        else{
     622                /*Check to prevent dividing by zero if vmag==0*/
     623                if(vmag==0. && (s-1.)<0.) alpha2=0.;
     624                else alpha2=drag_coefficient*drag_coefficient*pow(Neff,r)*pow(vmag,(s-1.));
     625        }
    552626
    553627        /*Assign output pointers:*/
    … …  
    695769
    696770        /*Compute alpha^2*/
    697         if(ub<1e-10){
     771        if((ub<1e-10) ||(N==0.0)){
    698772                alpha2 = 0.;
    699773        }
    … …  
    735809
    736810                alpha2 = alpha2/ub;
     811        }
     812
     813        /*Assign output pointers:*/
     814        *palpha2=alpha2;
     815}/*}}}*/
     816void Friction::GetAlpha2Coulomb2(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
     817
     818        /*This routine calculates the basal friction coefficient
     819         *
     820         *               C^2 |u_b|^(m-1) * (.5*N)
     821         * alpha2= ___________________________________
     822         *          (C^(2/m) |u_b| + (0.5*N)^(1/m) )^m
     823         *
     824         * */
     825
     826        /*diverse: */
     827        IssmDouble  C,m,alpha2;
     828
     829        /*Recover parameters: */
     830        element->GetInputValue(&C,gauss,FrictionCEnum);
     831        element->GetInputValue(&m,gauss,FrictionMEnum);
     832
     833        /*Get effective pressure and velocity magnitude*/
     834        IssmDouble N  = EffectivePressure(gauss);
     835        IssmDouble ub = VelMag(gauss);
     836
     837        /*Compute alpha^2*/
     838        if(ub<1e-10){
     839                alpha2 = 0.;
     840        }
     841        else{
     842                alpha2= (pow(C,2)*pow(ub,m-1.)*(0.5*N)) / pow(pow(C,2./m)*ub + pow(0.5*N,1./m),m);
    737843        }
    738844
    … …  
    754860        element->parameters->FindParam(&Neff_limit,FrictionEffectivePressureLimitEnum);
    755861
     862        /*Compute ice pressure*/
     863        p_ice = IcePressure(gauss);
     864
    756865        /*From base and thickness, compute effective pressure when drag is viscous, or get Neff from forcing:*/
    757866        switch(coupled_flag){
    758867                case 0:{
    759                         element->GetInputValue(&thickness, gauss,ThicknessEnum);
    760868                        element->GetInputValue(&base, gauss,BaseEnum);
    761869                        element->GetInputValue(&sealevel, gauss,SealevelEnum);
    762870                        IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum);
    763                         IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
    764871                        IssmDouble gravity   = element->FindParam(ConstantsGEnum);
    765                         p_ice   = gravity*rho_ice*thickness;
    766872                        p_water = rho_water*gravity*(sealevel-base);
    767873                        Neff = p_ice - p_water;
    … …  
    769875                        break;
    770876                case 1:{
    771                         element->GetInputValue(&thickness, gauss,ThicknessEnum);
    772                         IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
    773                         IssmDouble gravity   = element->FindParam(ConstantsGEnum);
    774                         p_ice   = gravity*rho_ice*thickness;
    775877                        p_water = 0.;
    776878                        Neff = p_ice - p_water;
    … …  
    778880                        break;
    779881                case 2:{
    780                         element->GetInputValue(&thickness, gauss,ThicknessEnum);
    781882                        element->GetInputValue(&base, gauss,BaseEnum);
    782883                        element->GetInputValue(&sealevel, gauss,SealevelEnum);
    783884                        IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum);
    784                         IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
    785885                        IssmDouble gravity   = element->FindParam(ConstantsGEnum);
    786                         p_ice   = gravity*rho_ice*thickness;
    787886                        p_water = max(0.,rho_water*gravity*(sealevel-base));
    788887                        Neff = p_ice - p_water;
    … …  
    791890                case 3:{
    792891                        element->GetInputValue(&Neff,gauss,FrictionEffectivePressureEnum);
    793                         element->GetInputValue(&thickness, gauss,ThicknessEnum);
    794                         IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
    795                         IssmDouble gravity   = element->FindParam(ConstantsGEnum);
    796                         p_ice   = gravity*rho_ice*thickness;
    797892                }
    798893                        break;
    799894                case 4:{
    800895                        element->GetInputValue(&Neff,gauss,EffectivePressureEnum);
    801                         element->GetInputValue(&thickness, gauss,ThicknessEnum);
    802                         IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
    803                         IssmDouble gravity   = element->FindParam(ConstantsGEnum);
    804                         p_ice   = gravity*rho_ice*thickness;
    805896                }
    806897                        break;
    … …  
    816907
    817908}/*}}}*/
     909IssmDouble Friction::IcePressure(Gauss* gauss){/*{{{*/
     910        /*Get ice pressure*/
     911
     912        IssmDouble  thickness,p_ice;
     913        /*Recover Inputs and Parameters*/
     914        element->GetInputValue(&thickness, gauss,ThicknessEnum);
     915        IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
     916        IssmDouble gravity = element->FindParam(ConstantsGEnum);
     917
     918        /*Compute*/
     919        p_ice = gravity*rho_ice*thickness;
     920
     921        /*Return ice pressure*/
     922        return p_ice;
     923
     924}/*}}}*/
     925IssmDouble Friction::SubglacialWaterPressure(Gauss* gauss){/*{{{*/
     926        /*Get water pressure as a function of  flag */
     927
     928        int         coupled_flag;
     929        IssmDouble  base,sealevel,p_water;
     930
     931        /*Recover parameters: */
     932        element->parameters->FindParam(&coupled_flag,FrictionCouplingEnum);
     933
     934        switch(coupled_flag){
     935                case 0:{
     936                        element->GetInputValue(&base, gauss,BaseEnum);
     937                        element->GetInputValue(&sealevel, gauss,SealevelEnum);
     938                        IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum);
     939                        IssmDouble gravity   = element->FindParam(ConstantsGEnum);
     940                        p_water = rho_water*gravity*(sealevel-base);
     941                }
     942                        break;
     943                case 1:{
     944                        p_water = 0.;
     945                }
     946                        break;
     947                case 2:{
     948                        element->GetInputValue(&base, gauss,BaseEnum);
     949                        element->GetInputValue(&sealevel, gauss,SealevelEnum);
     950                        IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum);
     951                        IssmDouble gravity   = element->FindParam(ConstantsGEnum);
     952                        p_water = max(0.,rho_water*gravity*(sealevel-base));
     953                }
     954                        break;
     955                case 3:{
     956                        _error_("water pressure not computed for coupling==3 in friction law");
     957                }
     958                        break;
     959                case 4:{
     960                        _error_("water pressure not computed for coupling==4 in friction law");
     961                }
     962                        break;
     963                default:
     964                        _error_("not supported");
     965        }
     966
     967        /*Return water pressure*/
     968        return p_water;
     969
     970}/*}}}*/
    818971IssmDouble Friction::VelMag(Gauss* gauss){/*{{{*/
    819         /*Get effective pressure as a function of  flag */
    820 
     972        /*Get the velocity magnitude as a function of flag */
    821973
    822974        /*diverse*/
    823975        IssmDouble vx,vy,vz,vmag;
    824976
    825         switch(dim){
    826                 case 1:
    827                         element->GetInputValue(&vx,gauss,VxEnum);
    828                         vmag=sqrt(vx*vx);
    829                         break;
    830                 case 2:
    831                         element->GetInputValue(&vx,gauss,VxEnum);
    832                         element->GetInputValue(&vy,gauss,VyEnum);
    833                         vmag=sqrt(vx*vx+vy*vy);
    834                         break;
    835                 case 3:
    836                         element->GetInputValue(&vx,gauss,VxEnum);
    837                         element->GetInputValue(&vy,gauss,VyEnum);
    838                         element->GetInputValue(&vz,gauss,VzEnum);
    839                         vmag=sqrt(vx*vx+vy*vy+vz*vz);
    840                         break;
    841                 default:
    842                         _error_("not supported");
    843         }
    844 
     977        this->vx_input->GetInputValue(&vx, gauss);
     978        this->vy_input->GetInputValue(&vy, gauss);
     979
     980        if ((this->vz_input == NULL) || (this->apply_dim<3.)) vz = 0.0;
     981        else this->vz_input->GetInputValue(&vz, gauss);
     982
     983        if (this->apply_dim<2.) vy = 0.0;
     984
     985        vmag = sqrt(vx*vx+vy*vy+vz*vz);
    845986        return vmag;
    846 
    847 }/*}}}*/
     987}/*}}}*/
     988void Friction::GetBasalSlidingSpeeds(IssmDouble* pvx, Gauss* gauss){/*{{{*/
     989
     990        this->vx_input->GetInputValue(pvx, gauss);
     991        /*Checks*/
     992        _assert_(!xIsNan<IssmDouble>(*pvx));
     993        _assert_(!xIsInf<IssmDouble>(*pvx));
     994}/*}}}*/
     995void Friction::GetBasalSlidingSpeeds(IssmDouble* pvx, IssmDouble* pvy, Gauss* gauss){/*{{{*/
     996
     997        this->vx_input->GetInputValue(pvx, gauss);
     998        this->vy_input->GetInputValue(pvy, gauss);
     999        /*Checks*/
     1000        _assert_(!xIsNan<IssmDouble>(*pvx));
     1001        _assert_(!xIsInf<IssmDouble>(*pvx));
     1002        _assert_(!xIsNan<IssmDouble>(*pvy));
     1003        _assert_(!xIsInf<IssmDouble>(*pvy));
     1004}/*}}}*/
     1005void Friction::GetBasalSlidingSpeeds(IssmDouble* pvx, IssmDouble* pvy, IssmDouble* pvz, Gauss* gauss){/*{{{*/
     1006
     1007        this->vx_input->GetInputValue(pvx, gauss);
     1008        this->vy_input->GetInputValue(pvy, gauss);
     1009        this->vz_input->GetInputValue(pvz, gauss);
     1010        /*Checks*/
     1011        _assert_(!xIsNan<IssmDouble>(*pvx));
     1012        _assert_(!xIsInf<IssmDouble>(*pvx));
     1013        _assert_(!xIsNan<IssmDouble>(*pvy));
     1014        _assert_(!xIsInf<IssmDouble>(*pvy));
     1015        _assert_(!xIsNan<IssmDouble>(*pvz));
     1016        _assert_(!xIsInf<IssmDouble>(*pvz));
     1017}/*}}}*/

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

    r25836 r26744  
    1414
    1515        public:
    16                 Element* element;
    17                 int      dim;
    18                 int      law;
     16                Element    *element;
     17                int         law;
     18                int         domaintype;
     19                IssmDouble  apply_dim;
     20                Input      *vx_input;
     21                Input      *vy_input;
     22                Input      *vz_input;
    1923
    2024                /*methods: */
    2125                Friction();
    22                 Friction(Element* element_in,int dim_in);
     26                Friction(Element* element_in);
     27                Friction(Element* element_in, int dim);
     28                Friction(Element* element_in, IssmDouble dim);
    2329                ~Friction();
    2430
    … …  
    3238                void  GetAlpha2(IssmDouble* palpha2,Gauss* gauss);
    3339                void  GetAlpha2Coulomb(IssmDouble* palpha2,Gauss* gauss);
     40                void  GetAlpha2Coulomb2(IssmDouble* palpha2,Gauss* gauss);
    3441                void  GetAlpha2Hydro(IssmDouble* palpha2,Gauss* gauss);
    3542                void  GetAlpha2Josh(IssmDouble* palpha2,Gauss* gauss);
    … …  
    4552
    4653                IssmDouble EffectivePressure(Gauss* gauss);
     54                IssmDouble IcePressure(Gauss* gauss);
     55                IssmDouble SubglacialWaterPressure(Gauss* gauss);
    4756                IssmDouble VelMag(Gauss* gauss);
     57                void GetBasalSlidingSpeeds(IssmDouble* pvx, Gauss* gauss);
     58                void GetBasalSlidingSpeeds(IssmDouble* pvx, IssmDouble* pvy, Gauss* gauss);
     59                void GetBasalSlidingSpeeds(IssmDouble* pvx, IssmDouble* pvy, IssmDouble* pvz, Gauss* gauss);
    4860};
    4961

  • issm/trunk/src/c/classes/Loads/Moulin.cpp

    r25836 r26744  
    291291                switch(set2_enum){
    292292                        case FsetEnum:
    293                                 if(node->fsize){
     293                                if(node->FSize()){
    294294                                        if(this->node->IsClone())
    295295                                         o_nz += 1;
    … …  
    307307                                break;
    308308                        case SsetEnum:
    309                                 if(node->ssize){
     309                                if(node->SSize()){
    310310                                        if(this->node->IsClone())
    311311                                         o_nz += 1;

  • issm/trunk/src/c/classes/Loads/Neumannflux.cpp

    r25836 r26744  
    306306                        switch(set2_enum){
    307307                                case FsetEnum:
    308                                         if(nodes[i]->fsize){
     308                                        if(nodes[i]->FSize()){
    309309                                                if(this->nodes[i]->IsClone())
    310310                                                 o_nz += 1;
    … …  
    322322                                        break;
    323323                                case SsetEnum:
    324                                         if(nodes[i]->ssize){
     324                                        if(nodes[i]->SSize()){
    325325                                                if(this->nodes[i]->IsClone())
    326326                                                 o_nz += 1;
    … …  
    354354        Tria*  tria=(Tria*)element;
    355355        _assert_(tria->FiniteElement()==P1Enum);
    356         if(!tria->IsIceInElement() || tria->IsFloating()) return NULL;
     356        if(!tria->IsIceInElement() || tria->IsAllFloating()) return NULL;
    357357
    358358        /*Initialize Element vector and other vectors*/
    … …  
    396396        Tria*  tria=(Tria*)element;
    397397        _assert_(tria->FiniteElement()==P1Enum);
    398         if(!tria->IsIceInElement() || tria->IsFloating()) return NULL;
     398        if(!tria->IsIceInElement() || tria->IsAllFloating()) return NULL;
    399399
    400400        /*Initialize Element vector and other vectors*/

  • issm/trunk/src/c/classes/Loads/Numericalflux.cpp

    r25836 r26744  
    430430                        switch(set2_enum){
    431431                                case FsetEnum:
    432                                         if(nodes[i]->fsize){
     432                                        if(nodes[i]->FSize()){
    433433                                                if(this->nodes[i]->IsClone())
    434434                                                 o_nz += 1;
    … …  
    446446                                        break;
    447447                                case SsetEnum:
    448                                         if(nodes[i]->ssize){
     448                                        if(nodes[i]->SSize()){
    449449                                                if(this->nodes[i]->IsClone())
    450450                                                 o_nz += 1;

  • issm/trunk/src/c/classes/Loads/Pengrid.cpp

    r25836 r26744  
    294294                switch(set2_enum){
    295295                        case FsetEnum:
    296                                 if(node->fsize){
     296                                if(node->FSize()){
    297297                                        if(this->node->IsClone())
    298298                                         o_nz += 1;
    … …  
    310310                                break;
    311311                        case SsetEnum:
    312                                 if(node->ssize){
     312                                if(node->SSize()){
    313313                                        if(this->node->IsClone())
    314314                                         o_nz += 1;

  • issm/trunk/src/c/classes/Loads/Penpair.cpp

    r25836 r26744  
    235235                        switch(set2_enum){
    236236                                case FsetEnum:
    237                                         if(nodes[i]->fsize){
     237                                        if(nodes[i]->FSize()){
    238238                                                if(this->nodes[i]->IsClone())
    239239                                                 o_nz += 1;
    … …  
    251251                                        break;
    252252                                case SsetEnum:
    253                                         if(nodes[i]->ssize){
     253                                        if(nodes[i]->SSize()){
    254254                                                if(this->nodes[i]->IsClone())
    255255                                                 o_nz += 1;
    … …  
    355355                                default: _error_("Approximation "<<EnumToStringx(approximation1)<<" not supported yet");
    356356                        }
     357                case MLHOApproximationEnum:
     358                        switch(approximation1){
     359                                case MLHOApproximationEnum:   return PenaltyCreateKMatrixStressbalanceSSAHO(kmax);
     360                                default: _error_("Approximation "<<EnumToStringx(approximation1)<<" not supported yet");
     361                        }
     362                case L1L2ApproximationEnum:
     363                        switch(approximation1){
     364                                case L1L2ApproximationEnum:   return PenaltyCreateKMatrixStressbalanceSSAHO(kmax);
     365                                default: _error_("Approximation "<<EnumToStringx(approximation1)<<" not supported yet");
     366                        }
    357367                case FSvelocityEnum:
    358368                        switch(approximation1){

  • issm/trunk/src/c/classes/Loads/Riftfront.cpp

    r25836 r26744  
    392392                        switch(set2_enum){
    393393                                case FsetEnum:
    394                                         if(nodes[i]->fsize){
     394                                        if(nodes[i]->FSize()){
    395395                                                if(this->nodes[i]->IsClone())
    396396                                                 o_nz += 1;
    … …  
    408408                                        break;
    409409                                case SsetEnum:
    410                                         if(nodes[i]->ssize){
     410                                        if(nodes[i]->SSize()){
    411411                                                if(this->nodes[i]->IsClone())
    412412                                                 o_nz += 1;
Note: See TracChangeset for help on using the changeset viewer.