source: issm/oecreview/Archive/21337-21723/ISSM-21429-21430.diff@ 21726

Last change on this file since 21726 was 21726, checked in by Mathieu Morlighem, 8 years ago

CHG added Archive/21337-21723
File size: 23.2 KB

  • ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp

     
    1818        int         penalty_lock;
    1919        int         hydro_maxiter;
    2020        bool        isefficientlayer;
    21         IssmDouble  sedimentlimit;
    2221        IssmDouble  penalty_factor;
    23         IssmDouble  leakagefactor;
    2422        IssmDouble  rel_tol;
     23        IssmDouble  leakagefactor;
     24        IssmDouble  sedimentlimit;
    2525
    2626        /*retrieve some parameters: */
    2727        iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
    … …  
    2929        /*Now, do we really want DC?*/
    3030        if(hydrology_model!=HydrologydcEnum) return;
    3131
    32         iomodel->FetchData(&isefficientlayer,   "md.hydrology.isefficientlayer");
    3332        iomodel->FetchData(&sedimentlimit_flag, "md.hydrology.sedimentlimit_flag" );
    3433        iomodel->FetchData(&transfer_flag,      "md.hydrology.transfer_flag" );
    3534        iomodel->FetchData(&penalty_factor,     "md.hydrology.penalty_factor" );
    36         iomodel->FetchData(&rel_tol,            "md.hydrology.rel_tol" );
    37         iomodel->FetchData(&penalty_lock,       "md.hydrology.penalty_lock" );
     35        iomodel->FetchData(&isefficientlayer,   "md.hydrology.isefficientlayer");
    3836        iomodel->FetchData(&hydro_maxiter,      "md.hydrology.max_iter" );
     37        iomodel->FetchData(&penalty_lock,       "md.hydrology.penalty_lock" );
     38        iomodel->FetchData(&rel_tol,            "md.hydrology.rel_tol" );
    3939
    40         if(sedimentlimit_flag==1){
    41                 iomodel->FetchData(&sedimentlimit,"md.hydrology.sedimentlimit");
    42                 parameters->AddObject(new DoubleParam(HydrologydcSedimentlimitEnum,sedimentlimit));
    43         }
    44 
    45         if(transfer_flag==1){
    46                 iomodel->FetchData(&leakagefactor,"md.hydrology.leakage_factor");
    47                 parameters->AddObject(new DoubleParam(HydrologydcLeakageFactorEnum,leakagefactor));
    48         }
    49 
    50         parameters->AddObject(new DoubleParam(HydrologydcPenaltyFactorEnum,penalty_factor));
    5140        parameters->AddObject(new IntParam(HydrologyModelEnum,hydrology_model));
    52         parameters->AddObject(new BoolParam(HydrologydcIsefficientlayerEnum,isefficientlayer));
    5341        parameters->AddObject(new IntParam(HydrologydcSedimentlimitFlagEnum,sedimentlimit_flag));
    5442        parameters->AddObject(new IntParam(HydrologydcTransferFlagEnum,transfer_flag));
    55         parameters->AddObject(new DoubleParam(HydrologydcRelTolEnum,rel_tol));
    5643        parameters->AddObject(new IntParam(HydrologydcPenaltyLockEnum,penalty_lock));
    5744        parameters->AddObject(new IntParam(HydrologydcMaxIterEnum,hydro_maxiter));
     45        parameters->AddObject(new BoolParam(HydrologydcIsefficientlayerEnum,isefficientlayer));
     46        parameters->AddObject(new DoubleParam(HydrologydcPenaltyFactorEnum,penalty_factor));
     47        parameters->AddObject(new DoubleParam(HydrologydcRelTolEnum,rel_tol));
     48        if(transfer_flag==1){
     49                iomodel->FetchData(&leakagefactor,"md.hydrology.leakage_factor");
     50                parameters->AddObject(new DoubleParam(HydrologydcLeakageFactorEnum,leakagefactor));
     51        }
     52        if(sedimentlimit_flag==1){
     53                iomodel->FetchData(&sedimentlimit,"md.hydrology.sedimentlimit");
     54                parameters->AddObject(new DoubleParam(HydrologydcSedimentlimitEnum,sedimentlimit));
     55        }
    5856}/*}}}*/
    5957
    6058void HydrologyDCInefficientAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
    … …  
    108106        /*Now, do we really want DC?*/
    109107        if(hydrology_model!=HydrologydcEnum) return;
    110108
    111         if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
     109        if(iomodel->domaintype!=Domain2DhorizontalEnum){
     110                iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
     111        }
    112112        ::CreateNodes(nodes,iomodel,HydrologyDCInefficientAnalysisEnum,P1Enum);
    113113        iomodel->DeleteData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
    114114}/*}}}*/
    … …  
    130130        iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
    131131        if(hydrology_model!=HydrologydcEnum) return;
    132132
    133         if(iomodel->domaintype==Domain3DEnum) iomodel->FetchData(1,"md.mesh.vertexonbase");
    134 
     133        if(iomodel->domaintype==Domain3DEnum){
     134                iomodel->FetchData(1,"md.mesh.vertexonbase");
     135        }
    135136        //create penalties for nodes: no node can have water above the max
    136137        CreateSingleNodeToElementConnectivity(iomodel);
    137138        for(int i=0;i<iomodel->numberofvertices;i++){
    … …  
    189190        bool        active_element,isefficientlayer;
    190191        IssmDouble  D_scalar,Jdet,dt;
    191192        IssmDouble  sediment_transmitivity;
    192         IssmDouble  prestep_head,base_elev;
     193        IssmDouble  base_elev;
    193194        IssmDouble  transfer,sediment_storing;
    194195        IssmDouble *xyz_list  = NULL;
    195196
    … …  
    300301        IssmDouble Jdet;
    301302
    302303        IssmDouble *xyz_list             = NULL;
    303         Input*      old_wh_input         = NULL;
    304304        Input*      active_element_input = NULL;
    305305
    306306        /*Fetch number of nodes and dof for this finite element*/
    … …  
    320320        Input* thick_input        = basalelement->GetInput(ThicknessEnum);
    321321        Input* base_input                 = basalelement->GetInput(BaseEnum);
    322322        Input* water_input        = basalelement->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(water_input);
    323         if(dt!= 0.){old_wh_input = basalelement->GetInput(SedimentHeadOldEnum);                     _assert_(old_wh_input);}
    324323
    325324        /*Transfer related Inputs*/
    326325        if(isefficientlayer){
    327326                active_element_input = basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input);
    328327        }
    329328
    330 
    331329        /* Start  looping on the number of gaussian points: */
    332330        Gauss* gauss=basalelement->NewGauss(2);
    333331        for(int ig=gauss->begin();ig<gauss->end();ig++){
    334332                gauss->GaussPoint(ig);
    335        
    336333                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
    337334                basalelement->NodalFunctions(basis,gauss);
    338335
    … …  
    358355                        }
    359356                }
    360357
    361 
    362358                /*Transient and transfer terms*/
    363359                if(dt!=0.){
    364                         old_wh_input    ->GetInputValue(&water_head,gauss);
    365                         sediment_storing = SedimentStoring(basalelement,gauss,sed_head_input,base_input);//sed_head_input
     360                        sediment_storing = SedimentStoring(basalelement,gauss,sed_head_input,base_input);
    366361                        if(isefficientlayer){
    367362                                /*Dealing with the sediment part of the transfer term*/
    368363                                active_element_input->GetInputValue(&active_element);
    … …  
    427422
    428423void HydrologyDCInefficientAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
    429424
    430         int        domaintype;
    431         bool       converged;
    432         int*       doflist=NULL;
    433         Element*   basalelement=NULL;
     425        /*Intermediaries*/
     426        int                                             domaintype;
     427        Element*                        basalelement;
     428        bool                                    converged;
     429        int*                                    doflist =       NULL;
    434430
     431        /*Get basal element*/
    435432        element->FindParam(&domaintype,DomainTypeEnum);
    436         if(domaintype!=Domain2DhorizontalEnum){
    437                 if(!element->IsOnBase()) return;
    438                 basalelement=element->SpawnBasalElement();
     433        switch(domaintype){
     434                case Domain2DhorizontalEnum:
     435                        basalelement = element;
     436                        break;
     437                case Domain3DEnum:
     438                        if(!element->IsOnBase()) return NULL;
     439                        basalelement = element->SpawnBasalElement();
     440                        break;
     441                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
    439442        }
    440         else{
    441                 basalelement = element;
    442         }
    443443
     444
    444445        /*Fetch number of nodes for this finite element*/
    445446        int numnodes = basalelement->GetNumberOfNodes();
    446447
    … …  
    478479                kappa=kmax*pow(10.,penalty_factor);
    479480
    480481                for(int i=0;i<numnodes;i++){
    481 
    482482                        GetHydrologyDCInefficientHmax(&h_max,basalelement,basalelement->GetNode(i));
    483483                        if(values[i]>h_max) {
    484484                                residual[i] = kappa*(values[i]-h_max);
    … …  
    526526        sed_head_input->GetInputValue(&prestep_head,gauss);
    527527        water_sheet=max(0.0,(prestep_head-(base_elev-sediment_thickness)));
    528528        specific_porosity=sediment_porosity-rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity));
    529         if (water_sheet<=0.9*sediment_thickness){//porosity for unconfined region
     529        //porosity for unconfined region
     530        if (water_sheet<=0.9*sediment_thickness){
    530531                sediment_storing=sediment_porosity;
    531532        }
    532         else if((water_sheet<sediment_thickness) && (water_sheet>0.9*sediment_thickness)){//continuity ramp
     533        //continuity ramp
     534        else if((water_sheet<sediment_thickness) && (water_sheet>0.9*sediment_thickness)){
    533535                sediment_storing=(sediment_thickness-water_sheet)*specific_porosity/(0.1*sediment_thickness)+
    534536                        rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity));
    535537        }
    536         else{//storing coefficient for confined
     538        //storing coefficient for confined
     539        else{
    537540                sediment_storing=rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity));
    538541        }
    539542        return sediment_storing;
    540543}/*}}}*/
     544
    541545IssmDouble HydrologyDCInefficientAnalysis::SedimentTransmitivity(Element* element,Gauss* gauss,Input* sed_head_input, Input* base_input,Input* SedTrans_input){/*{{{*/
    542546        IssmDouble sediment_transmitivity;
    543547        IssmDouble FullLayer_transmitivity;
    … …  
    573577                element->FindParam(&h_max,HydrologydcSedimentlimitEnum);
    574578                break;
    575579        case 2:
    576        
    577580                rho_water = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
    578581                rho_ice   = element->GetMaterialParameter(MaterialsRhoIceEnum);
    579 
    580582                _assert_(thick_input);
    581583                _assert_(base_input);
    582 
    583584                /*Compute max*/
    584585                thick_input->GetInputValue(&thickness,gauss);
    585586                base_input->GetInputValue(&bed,gauss);
    … …  
    633634IssmDouble HydrologyDCInefficientAnalysis::GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* thick_input, Input* base_input){/*{{{*/
    634635
    635636        int transfermethod;
    636         IssmDouble hmax;
    637         IssmDouble epl_head,sediment_head;
    638637        IssmDouble leakage,transfer;
    639         IssmDouble continuum, factor;
    640        
    641638        element->FindParam(&transfermethod,HydrologydcTransferFlagEnum);
    642639
    643640        /*Switch between the different transfer methods cases*/
    … …  
    647644                transfer=0.0;
    648645                break;
    649646        case 1:
    650                 _assert_(sed_head_input);
    651                 _assert_(epl_head_input);
    652                
    653                 sed_head_input->GetInputValue(&sediment_head,gauss);
    654                 epl_head_input->GetInputValue(&epl_head,gauss);
    655647                element->FindParam(&leakage,HydrologydcLeakageFactorEnum);
    656 
    657                 hmax=GetHydrologyDCInefficientHmax(element, gauss, thick_input, base_input);
    658        
    659                 //Computing continuum function to apply to transfer term, transfer is null only if
    660                 //epl_head>sediment_head AND sediment_head>h_max
    661                 //let's try without that for a while
    662                 /* continuum=((1.0/(1.0+exp(-20.0*(sediment_head-epl_head)))))+(1.0/(1.0+exp(-20.0*(hmax-sediment_head)))); */
    663                 /* factor=min(continuum,1.0); */
    664                 /* transfer=leakage*factor; */
    665648                transfer=leakage;
    666 
    667649                break;
    668650        default:
    669651                _error_("no case higher than 1 for the Transfer method");
    … …  
    674656IssmDouble HydrologyDCInefficientAnalysis::GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* thick_input, Input* base_input){/*{{{*/
    675657
    676658        int transfermethod;
    677         IssmDouble hmax;
    678         IssmDouble epl_head,sediment_head;
     659        IssmDouble epl_head;
    679660        IssmDouble leakage,transfer;
    680         IssmDouble continuum, factor;
    681        
    682661        element->FindParam(&transfermethod,HydrologydcTransferFlagEnum);
    683662
    684663        /*Switch between the different transfer methods cases*/
    … …  
    688667                transfer=0.0;
    689668                break;
    690669        case 1:
    691                 _assert_(sed_head_input);
    692670                _assert_(epl_head_input);
    693                
    694                 sed_head_input->GetInputValue(&sediment_head,gauss);
    695671                epl_head_input->GetInputValue(&epl_head,gauss);
    696672                element->FindParam(&leakage,HydrologydcLeakageFactorEnum);
    697 
    698                 hmax=GetHydrologyDCInefficientHmax(element, gauss, thick_input, base_input);
    699                
    700                 //Computing continuum function to apply to transfer term, transfer is null only if
    701                 //epl_head>sediment_head AND sediment_head>h_max
    702                 //let's try without that for a while
    703                 /* continuum=((1.0/(1.0+exp(-20.0*(sediment_head-epl_head)))))+(1.0/(1.0+exp(-20.0*(hmax-sediment_head)))); */
    704                 /* factor=min(continuum,1.0); */
    705                 /* transfer=epl_head*leakage*factor; */
    706 
    707673                transfer=epl_head*leakage;
    708 
    709674                break;
    710675        default:
    711676                _error_("no case higher than 1 for the Transfer method");
    … …  
    730695                }
    731696                element->AddInput(new BoolInput(HydrologydcMaskEplactiveEltEnum,element_active));
    732697        }
     698        delete element;
    733699}/*}}}*/

  • ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp

     
    88int  HydrologyDCEfficientAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
    99        return 1;
    1010}/*}}}*/
     11
    1112void HydrologyDCEfficientAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
    1213
    1314        int         hydrology_model;
    … …  
    3233
    3334        iomodel->FetchData(&eplthickcomp,"md.hydrology.epl_thick_comp");
    3435        parameters->AddObject(new IntParam(HydrologydcEplThickCompEnum,eplthickcomp));
    35 
    36        
    3736}/*}}}*/
     37
    3838void HydrologyDCEfficientAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
    3939
    4040        bool   isefficientlayer;
    … …  
    8282        iomodel->FindConstant(&isefficientlayer,"md.hydrology.isefficientlayer");
    8383        if(!isefficientlayer) return;
    8484
    85         if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
     85        if(iomodel->domaintype!=Domain2DhorizontalEnum){
     86                iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
     87        }
    8688        ::CreateNodes(nodes,iomodel,HydrologyDCEfficientAnalysisEnum,P1Enum);
    8789        iomodel->DeleteData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
    8890}/*}}}*/
    … …  
    105107void HydrologyDCEfficientAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
    106108        /*Nothing for now*/
    107109
    108         /*Fetch parameters: */
     110        /*Do we really want DC?*/
    109111        int hydrology_model;
    110112        iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
    111113        if(hydrology_model!=HydrologydcEnum) return;
    112114
    113         if(iomodel->domaintype==Domain3DEnum) iomodel->FetchData(1,"md.mesh.vertexonbase");
     115        /*Do we want an efficient layer*/
     116        bool isefficientlayer;
     117        iomodel->FindConstant(&isefficientlayer,"md.hydrology.isefficientlayer");
     118        if(!isefficientlayer) return;
    114119
    115         //create penalties for nodes: no node can have water above the max
     120        /*Fetch parameters: */
     121        if(iomodel->domaintype==Domain3DEnum){
     122                iomodel->FetchData(1,"md.mesh.vertexonbase");
     123        }
     124
     125        //Add moulin inputs as loads
    116126        CreateSingleNodeToElementConnectivity(iomodel);
    117127        for(int i=0;i<iomodel->numberofvertices;i++){
    118128                if (iomodel->domaintype!=Domain3DEnum){
    … …  
    131141}/*}}}*/
    132142
    133143void HydrologyDCEfficientAnalysis::InitZigZagCounter(FemModel* femmodel){/*{{{*/
     144
    134145        int*   eplzigzag_counter =NULL;
    135        
    136146        eplzigzag_counter=xNewZeroInit<int>(femmodel->nodes->Size());
    137147        femmodel->parameters->AddObject(new IntVecParam(EplZigZagCounterEnum,eplzigzag_counter,femmodel->nodes->Size()));
    138148        xDelete<int>(eplzigzag_counter);
    … …  
    141151void HydrologyDCEfficientAnalysis::ResetCounter(FemModel* femmodel){/*{{{*/
    142152
    143153        int*     eplzigzag_counter=NULL;
    144         Element* element=NULL;
    145 
    146154        femmodel->parameters->FindParam(&eplzigzag_counter,NULL,EplZigZagCounterEnum);
    147155        for(int i=0;i<femmodel->nodes->Size();i++){
    148 
    149156                eplzigzag_counter[i]=0;
    150157        }
    151158        femmodel->parameters->SetParam(eplzigzag_counter,femmodel->nodes->Size(),EplZigZagCounterEnum);
    … …  
    298305
    299306        /*Check that all nodes are active, else return empty matrix*/
    300307        if(!active_element) {
    301         if(domaintype!=Domain2DhorizontalEnum){
     308                if(domaintype!=Domain2DhorizontalEnum){
    302309                        basalelement->DeleteMaterials();
    303310                        delete basalelement;
    304311                }
    … …  
    341348        Gauss* gauss=basalelement->NewGauss(2);
    342349        for(int ig=gauss->begin();ig<gauss->end();ig++){
    343350                gauss->GaussPoint(ig);
    344 
    345351                basalelement ->JacobianDeterminant(&Jdet,xyz_list,gauss);
    346352                basalelement ->NodalFunctions(basis,gauss);
    347353
    … …  
    395401
    396402void HydrologyDCEfficientAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
    397403
    398         bool active_element;
    399         int domaintype,i;
    400         Element*   basalelement=NULL;
     404        /*Intermediaries*/
     405        bool     active_element;
     406        int      domaintype;
     407        Element* basalelement;
    401408
     409        /*Get basal element*/
    402410        element->FindParam(&domaintype,DomainTypeEnum);
    403 
    404         if(domaintype!=Domain2DhorizontalEnum){
    405                 if(!element->IsOnBase()) return;
    406                 basalelement=element->SpawnBasalElement();
     411        switch(domaintype){
     412                case Domain2DhorizontalEnum:
     413                        basalelement = element;
     414                        break;
     415                case Domain3DEnum:
     416                        if(!element->IsOnBase()) return NULL;
     417                        basalelement = element->SpawnBasalElement();
     418                        break;
     419                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
    407420        }
    408         else{
    409                 basalelement = element;
    410         }
    411        
     421
    412422        /*Intermediary*/
    413423        int* doflist = NULL;
    414424
    … …  
    469479IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* thick_input, Input* base_input){/*{{{*/
    470480
    471481        int transfermethod;
    472         IssmDouble hmax;
    473         IssmDouble epl_head,sediment_head;
    474482        IssmDouble leakage,transfer;
    475         IssmDouble continuum, factor;
    476         HydrologyDCInefficientAnalysis* inefanalysis = NULL;
    477483
    478484        element->FindParam(&transfermethod,HydrologydcTransferFlagEnum);
    479485
    … …  
    484490                transfer=0.0;
    485491                break;
    486492        case 1:
    487                 _assert_(sed_head_input);
    488                 _assert_(epl_head_input);
    489                
    490                 inefanalysis = new HydrologyDCInefficientAnalysis();
    491                 hmax = inefanalysis->GetHydrologyDCInefficientHmax(element,gauss,thick_input,base_input);
    492                 delete inefanalysis;
    493 
    494                 sed_head_input->GetInputValue(&sediment_head,gauss);
    495                 epl_head_input->GetInputValue(&epl_head,gauss);
    496493                element->FindParam(&leakage,HydrologydcLeakageFactorEnum);
    497                
    498                 //Computing continuum function to apply to transfer term, transfer is null only if
    499                 // epl_head>sediment_head AND sediment_head>h_max
    500                 //let's try without that for a while
    501                 /* continuum=((1.0/(1.0+exp(-20.0*(sediment_head-epl_head)))))+(1.0/(1.0+exp(-20.0*(hmax-sediment_head)))); */
    502                 /* factor=min(continuum,1.0); */
    503                 /* transfer=leakage*factor; */
    504494                transfer=leakage;
    505495                break;
    506496        default:
    … …  
    512502IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* thick_input, Input* base_input){/*{{{*/
    513503
    514504        int transfermethod;
    515         IssmDouble hmax;
    516         IssmDouble epl_head,sediment_head;
     505        IssmDouble sediment_head;
    517506        IssmDouble leakage,transfer;
    518         IssmDouble continuum, factor;
    519507
    520         HydrologyDCInefficientAnalysis* inefanalysis = NULL;
    521 
    522508        element->FindParam(&transfermethod,HydrologydcTransferFlagEnum);
    523509
    524510        /*Switch between the different transfer methods cases*/
    … …  
    529515                break;
    530516        case 1:
    531517                _assert_(sed_head_input);
    532                 _assert_(epl_head_input);
    533 
    534                 inefanalysis = new HydrologyDCInefficientAnalysis();
    535                 hmax = inefanalysis->GetHydrologyDCInefficientHmax(element,gauss,thick_input,base_input);
    536                 delete inefanalysis;
    537                
    538518                sed_head_input->GetInputValue(&sediment_head,gauss);
    539                 epl_head_input->GetInputValue(&epl_head,gauss);
    540519                element->FindParam(&leakage,HydrologydcLeakageFactorEnum);
    541 
    542                 //Computing continuum function to apply to transfer term, transfer is null only if
    543                 // epl_head>sediment_head AND sediment_head>h_max
    544                 //let's try without that for a while
    545                 /* continuum=((1.0/(1.0+exp(-20.0*(sediment_head-epl_head)))))+(1.0/(1.0+exp(-20.0*(hmax-sediment_head)))); */
    546                 /* factor=min(continuum,1.0); */
    547                 /* transfer=sediment_head*leakage*factor; */
    548520                transfer=sediment_head*leakage;
    549 
    550521                break;
    551522        default:
    552523                _error_("no case higher than 1 for the Transfer method");
    … …  
    563534        IssmDouble  dt,A,B;
    564535        IssmDouble  EPLgrad2;
    565536        IssmDouble  EPL_N;
    566 
    567537       
    568538        femmodel->parameters->FindParam(&domaintype,DomainTypeEnum);
    569539
    … …  
    618588                element->GetInputListOnVertices(&bed[0],BaseEnum);
    619589               
    620590                if(!active_element){
    621                        
    622591                        /*Keeping thickness to initial value if EPL is not active*/
    623592                        for(int i=0;i<numnodes;i++){
    624593                                thickness[i]=init_thick;
    … …  
    626595                }
    627596                else{
    628597                        for(int i=0;i<numnodes;i++){
    629                                
    630598                                /*Compute first the effective pressure in the EPL*/
    631                                 //EPL_N=gravity*((rho_ice*ice_thickness[i])-(rho_water*max(0.0,(eplhead[i]-bed[i]))));
    632                                 //allowing EPLN larger than Pi
    633599                                EPL_N=gravity*((rho_ice*ice_thickness[i])-(rho_water*(eplhead[i]-bed[i])));
    634600                                if(EPL_N<0.0)EPL_N=0.0;
    635601                                /*Get then the square of the gradient of EPL heads*/
    636602                                EPLgrad2 = (epl_slopeX[i]*epl_slopeX[i])+(epl_slopeY[i]*epl_slopeY[i]);
    637                                
    638603                                /*And proceed to the real thing*/
    639                                 /* thickness[i] = old_thickness[i]*(1.0+((rho_water*gravity*dt)/(rho_ice*latentheat))*epl_conductivity*EPLgrad2- */
    640                                 /*                                                                                                                               (2.0*(A*dt/(pow(n,n)))*(pow(EPL_N,n)))); */
    641604                                thickness[i] = old_thickness[i]*
    642605                                        (2.0+((rho_water*gravity*dt)/(rho_ice*latentheat))*epl_conductivity*EPLgrad2-(2.0*(A*dt/(pow(n,n)))*(pow(EPL_N,n))))/
    643606                                        (2.0-((rho_water*gravity*dt)/(rho_ice*latentheat))*epl_conductivity*EPLgrad2+(2.0*(A*dt/(pow(n,n)))*(pow(EPL_N,n))));
    644                                 //thickness[i] = 0.5*(thickness[i]+old_thickness[i]);
    645607                                /*Take care of otherthikening*/
    646608                                if(thickness[i]>max_thick){
    647609                                        thickness[i] = max_thick;
    … …  
    694656        int         domaintype;
    695657        IssmDouble  h_max;
    696658        IssmDouble  sedheadmin;
    697         Element*    basalelement=NULL;
     659        Element*    basalelement;
    698660
    699661        /*Get basal element*/
    700662        element->FindParam(&domaintype,DomainTypeEnum);
    … …  
    703665                        basalelement = element;
    704666                        break;
    705667                case Domain3DEnum:
    706                         if(!element->IsOnBase()) return;
     668                        if(!element->IsOnBase()) return NULL;
    707669                        basalelement = element->SpawnBasalElement();
    708670                        break;
    709671                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
    … …  
    736698                if (old_active[i]==0.){
    737699                        epl_thickness[i]=init_thick;
    738700                }
    739 
    740701                /*Now starting to look at the activations*/
    741702                if(residual[i]>0.){
    742703                        vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
    743                         if(old_active[i]==0.)   recurence->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
     704                        if(old_active[i]==0.){
     705                                recurence->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
     706                        }
    744707                }
    745708                /*If mask was already one, keep one or colapse*/
    746709                else if(old_active[i]>0.){
    … …  
    758721                                /*Increase of the domain is on the downstream node in term of sediment head*/
    759722                                if(sedhead[j] == sedheadmin){
    760723                                        vec_mask->SetValue(basalelement->nodes[j]->Sid(),1.,INS_VAL);
    761                                         if(old_active[i]==0.)   recurence->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
     724                                        if(old_active[i]==0.){
     725                                                recurence->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
     726                                        }
    762727                                }
    763728                        }
    764729                }
    … …  
    777742void HydrologyDCEfficientAnalysis::HydrologyEPLGetActive(Vector<IssmDouble>* active_vec, Element* element){
    778743        /*Constants*/
    779744        int      domaintype;
    780         Element*   basalelement=NULL;
     745        Element*   basalelement;
    781746
    782747        /*Get basal element*/
    783748        element->FindParam(&domaintype,DomainTypeEnum);
    … …  
    786751                        basalelement = element;
    787752                        break;
    788753                case Domain3DEnum:
    789                         if(!element->IsOnBase()) return;
     754                        if(!element->IsOnBase()) return NULL;
    790755                        basalelement = element->SpawnBasalElement();
    791756                        break;
    792757                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
    … …  
    795760        const int   numnodes = basalelement->GetNumberOfNodes();
    796761        IssmDouble  flag     = 0.;
    797762        IssmDouble* active   = xNew<IssmDouble>(numnodes);
     763        bool active_element;
    798764
    799765        /*Pass the activity mask from elements to nodes*/
    800766        basalelement->GetInputListOnVertices(&active[0],HydrologydcMaskEplactiveNodeEnum);
    801         bool active_element;
    802767        Input*  active_element_input=basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input);           
    803768        active_element_input->GetInputValue(&active_element);
    804769       
Note: See TracBrowser for help on using the repository browser.