Ignore:
Timestamp:
04/01/20 21:54:40 (5 years ago)
Author:
Mathieu Morlighem
Message:

merged trunk-jpl and trunk for revision 24684

Location:
issm/trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk

  • issm/trunk/src

  • issm/trunk/src/c/analyses/HydrologyDCEfficientAnalysis.cpp

    r24313 r24686  
    99        return 1;
    1010}/*}}}*/
    11 
    1211void HydrologyDCEfficientAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
    1312
    … …  
    3534        parameters->AddObject(new IntParam(HydrologydcEplThickCompEnum,eplthickcomp));
    3635}/*}}}*/
    37 
    38 void HydrologyDCEfficientAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
     36void HydrologyDCEfficientAnalysis::UpdateElements(Elements* elements,Inputs2* inputs2,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
    3937
    4038        bool   isefficientlayer;
    … …  
    5452                if(iomodel->my_elements[i]){
    5553                        Element* element=(Element*)elements->GetObjectByOffset(counter);
    56                         element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum);
     54                        element->Update(inputs2,i,iomodel,analysis_counter,analysis_type,P1Enum);
    5755                        counter++;
    5856                }
    5957        }
    60         iomodel->FetchDataToInput(elements,"md.geometry.thickness",ThicknessEnum);
    61         iomodel->FetchDataToInput(elements,"md.geometry.base",BaseEnum);
    62         iomodel->FetchDataToInput(elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
    63         iomodel->FetchDataToInput(elements,"md.initialization.epl_head",EplHeadSubstepEnum);
    64         iomodel->FetchDataToInput(elements,"md.initialization.sediment_head",SedimentHeadSubstepEnum);
    65         iomodel->FetchDataToInput(elements,"md.initialization.epl_thickness",HydrologydcEplThicknessSubstepEnum);
    66         iomodel->FetchDataToInput(elements,"md.hydrology.basal_moulin_input",HydrologydcBasalMoulinInputEnum);
     58        iomodel->FetchDataToInput(inputs2,elements,"md.geometry.thickness",ThicknessEnum);
     59        iomodel->FetchDataToInput(inputs2,elements,"md.geometry.base",BaseEnum);
     60        iomodel->FetchDataToInput(inputs2,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
     61        iomodel->FetchDataToInput(inputs2,elements,"md.initialization.epl_head",EplHeadSubstepEnum);
     62        iomodel->FetchDataToInput(inputs2,elements,"md.initialization.sediment_head",SedimentHeadSubstepEnum);
     63        iomodel->FetchDataToInput(inputs2,elements,"md.initialization.epl_thickness",HydrologydcEplThicknessSubstepEnum);
     64        iomodel->FetchDataToInput(inputs2,elements,"md.hydrology.basal_moulin_input",HydrologydcBasalMoulinInputEnum);
    6765        if(iomodel->domaintype!=Domain2DhorizontalEnum){
    68                 iomodel->FetchDataToInput(elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum);
    69                 iomodel->FetchDataToInput(elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum);
    70         }
    71 }/*}}}*/
    72 
     66                iomodel->FetchDataToInput(inputs2,elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum);
     67                iomodel->FetchDataToInput(inputs2,elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum);
     68        }
     69}/*}}}*/
    7370void HydrologyDCEfficientAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/
    7471
    … …  
    8986        iomodel->DeleteData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
    9087}/*}}}*/
    91 
    9288void HydrologyDCEfficientAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
    9389
    … …  
    104100        IoModelToConstraintsx(constraints,iomodel,"md.hydrology.spcepl_head",HydrologyDCEfficientAnalysisEnum,P1Enum);
    105101}/*}}}*/
    106 
    107102void HydrologyDCEfficientAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
    108103
    … …  
    139134        iomodel->DeleteData(1,"md.mesh.vertexonbase");
    140135}/*}}}*/
    141 
    142136void HydrologyDCEfficientAnalysis::InitZigZagCounter(FemModel* femmodel){/*{{{*/
    143137
    … …  
    147141        xDelete<int>(eplzigzag_counter);
    148142}/*}}}*/
    149 
    150143void HydrologyDCEfficientAnalysis::ResetCounter(FemModel* femmodel){/*{{{*/
    151144
    … …  
    163156        _error_("not implemented");
    164157}/*}}}*/
    165 
    166158ElementVector* HydrologyDCEfficientAnalysis::CreateDVector(Element* element){/*{{{*/
    167159        /*Default, return NULL*/
    168160        return NULL;
    169161}/*}}}*/
    170 
    171162ElementMatrix* HydrologyDCEfficientAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
    172163_error_("Not implemented");
    173164}/*}}}*/
    174 
    175165ElementMatrix* HydrologyDCEfficientAnalysis::CreateKMatrix(Element* element){/*{{{*/
    176166
    … …  
    193183        }
    194184
    195         Input* active_element_input = basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input);
    196         active_element_input->GetInputValue(&active_element);
     185        basalelement->GetInputValue(&active_element,HydrologydcMaskEplactiveEltEnum);
    197186
    198187        /*Check that all nodes are active, else return empty matrix*/
    … …  
    224213        basalelement->GetVerticesCoordinates(&xyz_list);
    225214        basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
    226 
    227         Input* epl_thick_input = basalelement->GetInput(HydrologydcEplThicknessSubstepEnum); _assert_(epl_thick_input);
    228         Input* epl_head_input   = basalelement->GetInput(EplHeadSubstepEnum);  _assert_(epl_head_input);
    229         Input* base_input                       = basalelement->GetInput(BaseEnum); _assert_(base_input);
     215        Input2* epl_thick_input = basalelement->GetInput2(HydrologydcEplThicknessSubstepEnum); _assert_(epl_thick_input);
     216        Input2* epl_head_input  = basalelement->GetInput2(EplHeadSubstepEnum);  _assert_(epl_head_input);
     217        Input2* base_input      = basalelement->GetInput2(BaseEnum); _assert_(base_input);
    230218
    231219        /* Start  looping on the number of gaussian points: */
    … …  
    279267        return Ke;
    280268}/*}}}*/
    281 
    282269ElementVector* HydrologyDCEfficientAnalysis::CreatePVector(Element* element){/*{{{*/
    283270
    … …  
    300287        }
    301288
    302         Input* active_element_input = basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input);
    303         active_element_input->GetInputValue(&active_element);
     289        basalelement->GetInputValue(&active_element,HydrologydcMaskEplactiveEltEnum);
    304290
    305291        /*Check that all nodes are active, else return empty matrix*/
    … …  
    319305        IssmDouble residual,connectivity;
    320306
    321         IssmDouble              *xyz_list                                = NULL;
    322         Input*                   old_wh_input                    = NULL;
    323         Input*                   surface_runoff_input = NULL;
     307        IssmDouble *xyz_list            = NULL;
     308        Input2     *old_wh_input        = NULL;
     309        Input2     *surface_runoff_input = NULL;
    324310
    325311        /*Fetch number of nodes and dof for this finite element*/
    … …  
    336322        basalelement ->FindParam(&smb_model,SmbEnum);
    337323
    338         Input*  epl_thick_input                  = basalelement->GetInput(HydrologydcEplThicknessSubstepEnum); _assert_(epl_thick_input);
    339         Input*  sed_head_input                   = basalelement->GetInput(SedimentHeadSubstepEnum); _assert_(sed_head_input);
    340         Input*  epl_head_input                   = basalelement->GetInput(EplHeadSubstepEnum); _assert_(epl_head_input);
    341         Input*  basal_melt_input                 = basalelement->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(basal_melt_input);
    342         Input*  residual_input                   = basalelement->GetInput(SedimentHeadResidualEnum); _assert_(residual_input);
    343         Input*  base_input                                       = basalelement->GetInput(BaseEnum); _assert_(base_input);
     324        Input2* epl_thick_input  = basalelement->GetInput2(HydrologydcEplThicknessSubstepEnum); _assert_(epl_thick_input);
     325        Input2* sed_head_input   = basalelement->GetInput2(SedimentHeadSubstepEnum); _assert_(sed_head_input);
     326        Input2* epl_head_input   = basalelement->GetInput2(EplHeadSubstepEnum); _assert_(epl_head_input);
     327        Input2* basal_melt_input = basalelement->GetInput2(BasalforcingsGroundediceMeltingRateEnum); _assert_(basal_melt_input);
     328        Input2* residual_input   = basalelement->GetInput2(SedimentHeadResidualEnum); _assert_(residual_input);
     329        Input2* base_input       = basalelement->GetInput2(BaseEnum); _assert_(base_input);
    344330
    345331        if(dt!= 0.){
    346                 old_wh_input = basalelement->GetInput(EplHeadOldEnum);            _assert_(old_wh_input);
     332                old_wh_input = basalelement->GetInput2(EplHeadOldEnum);            _assert_(old_wh_input);
    347333        }
    348334        if(smb_model==SMBgradientscomponentsEnum){
    349                 surface_runoff_input = basalelement->GetInput(SmbRunoffEnum); _assert_(surface_runoff_input);
     335                surface_runoff_input = basalelement->GetInput2(SmbRunoffEnum); _assert_(surface_runoff_input);
    350336        }
    351337
    … …  
    397383        return pe;
    398384}/*}}}*/
    399 
     385void HydrologyDCEfficientAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     386        /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
     387         * For node i, Bi can be expressed in the actual coordinate system
     388         * by:
     389         *       Bi=[ dN/dx ]
     390         *          [ dN/dy ]
     391         * where N is the finiteelement function for node i.
     392         *
     393         * We assume B has been allocated already, of size: 3x(NDOF2*numnodes)
     394         */
     395
     396        /*Fetch number of nodes for this finite element*/
     397        int numnodes = element->GetNumberOfNodes();
     398
     399        /*Get nodal functions derivatives*/
     400        IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
     401        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     402
     403        /*Build B: */
     404        for(int i=0;i<numnodes;i++){
     405                B[numnodes*0+i] = dbasis[0*numnodes+i];
     406                B[numnodes*1+i] = dbasis[1*numnodes+i];
     407        }
     408
     409        /*Clean-up*/
     410        xDelete<IssmDouble>(dbasis);
     411}/*}}}*/
    400412void HydrologyDCEfficientAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
    401413        element->GetSolutionFromInputsOneDof(solution,EplHeadSubstepEnum);
    402414}/*}}}*/
    403 
    404415void HydrologyDCEfficientAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
    405416        _error_("Not implemented yet");
    406417}/*}}}*/
    407 
    408418void HydrologyDCEfficientAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
    409419        /*Intermediaries*/
    … …  
    441451        }
    442452        /*Add input to the element: */
    443         element->AddBasalInput(EplHeadSubstepEnum,eplHeads,P1Enum);
     453        element->AddBasalInput2(EplHeadSubstepEnum,eplHeads,P1Enum);
     454
    444455        /*Free ressources:*/
    445456        xDelete<IssmDouble>(eplHeads);
    … …  
    447458        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    448459} /*}}}*/
    449 
    450460void HydrologyDCEfficientAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
    451461        /*Default, do nothing*/
    … …  
    454464
    455465/*Intermediaries*/
     466IssmDouble HydrologyDCEfficientAnalysis::EplStoring(Element* element,Gauss* gauss, Input2* epl_thick_input, Input2* epl_head_input, Input2* base_input){/*{{{*/
     467        IssmDouble epl_storing;
     468        IssmDouble water_sheet,storing;
     469        IssmDouble epl_thickness,prestep_head,base_elev;
     470        IssmDouble rho_freshwater        = element->FindParam(MaterialsRhoFreshwaterEnum);
     471        IssmDouble g                     = element->FindParam(ConstantsGEnum);
     472        IssmDouble epl_porosity                                  = element->FindParam(HydrologydcEplPorosityEnum);
     473        IssmDouble epl_compressibility   = element->FindParam(HydrologydcEplCompressibilityEnum);
     474        IssmDouble water_compressibility = element->FindParam(HydrologydcWaterCompressibilityEnum);
     475
     476        epl_thick_input->GetInputValue(&epl_thickness,gauss);
     477        epl_head_input->GetInputValue(&prestep_head,gauss);
     478        base_input->GetInputValue(&base_elev,gauss);
     479        water_sheet=max(0.0,(prestep_head-base_elev));
     480        storing=rho_freshwater*g*epl_porosity*epl_thickness*(water_compressibility+(epl_compressibility/epl_porosity));
     481
     482        /* //porosity for unconfined region */
     483        /* if (water_sheet<=0.9*epl_thickness){ */
     484        /*      epl_storing=epl_porosity; */
     485        /* } */
     486        /* //continuity ramp */
     487        /* else if((water_sheet<epl_thickness) && (water_sheet>0.9*epl_thickness)){ */
     488        /*      epl_storing=(epl_thickness-water_sheet)*(epl_porosity-storing)/(0.1*epl_thickness)+storing; */
     489        /* } */
     490        /* //storing coefficient for confined */
     491        /* else{ */
     492        /*      epl_storing=storing; */
     493        /* } */
     494        /* return epl_storing; */
     495        return storing;
     496}/*}}}*/
     497IssmDouble HydrologyDCEfficientAnalysis::EplTransmitivity(Element* element,Gauss* gauss, Input2* epl_thick_input, Input2* epl_head_input, Input2* base_input){/*{{{*/
     498        IssmDouble epl_transmitivity;
     499        IssmDouble water_sheet;
     500        IssmDouble epl_thickness,base_elev,prestep_head;
     501        IssmDouble epl_conductivity      = element->FindParam(HydrologydcEplConductivityEnum);
     502        epl_thick_input->GetInputValue(&epl_thickness,gauss);
     503        epl_head_input->GetInputValue(&prestep_head,gauss);
     504        base_input->GetInputValue(&base_elev,gauss);
     505
     506        water_sheet=max(0.0,(prestep_head-base_elev));
     507        epl_transmitivity=epl_conductivity*epl_thickness;
     508        //epl_transmitivity=max(1.0e-6,(epl_conductivity*min(water_sheet,epl_thickness)));
     509        return epl_transmitivity;
     510}/*}}}*/
     511void HydrologyDCEfficientAnalysis::GetHydrologyDCInefficientHmax(IssmDouble* ph_max,Element* element, Node* innode){/*{{{*/
     512
     513        int        hmax_flag;
     514        IssmDouble h_max;
     515        IssmDouble rho_ice,rho_water;
     516        IssmDouble thickness,bed;
     517        /*Get the flag to the limitation method*/
     518        element->FindParam(&hmax_flag,HydrologydcSedimentlimitFlagEnum);
     519
     520        /*Switch between the different cases*/
     521        switch(hmax_flag){
     522        case 0:
     523                h_max=1.0e+10;
     524                break;
     525        case 1:
     526                element->FindParam(&h_max,HydrologydcSedimentlimitEnum);
     527                break;
     528        case 2:
     529                /*Compute max*/
     530                rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
     531                rho_ice   = element->FindParam(MaterialsRhoIceEnum);
     532                element-> GetInputValue(&thickness,innode,ThicknessEnum);
     533                element-> GetInputValue(&bed,innode,BaseEnum);
     534                h_max=((rho_ice*thickness)/rho_water)+bed;
     535                break;
     536        case 3:
     537                _error_("Using normal stress  not supported yet");
     538                break;
     539        default:
     540                _error_("no case higher than 3 for SedimentlimitFlag");
     541        }
     542        /*Assign output pointer*/
     543        *ph_max=h_max;
     544}
     545/*}}}*/
    456546IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyKMatrixTransfer(Element* element){/*{{{*/
    457547
    … …  
    475565        return transfer;
    476566}/*}}}*/
    477 
    478 IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* sed_head_input){/*{{{*/
     567IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input2* sed_head_input){/*{{{*/
    479568
    480569        int transfermethod;
    … …  
    501590        return transfer;
    502591}/*}}}*/
    503 
    504592void HydrologyDCEfficientAnalysis::ComputeEPLThickness(FemModel* femmodel){/*{{{*/
    505593
    … …  
    534622                IssmDouble* bed           = xNew<IssmDouble>(numnodes);
    535623
    536                 Input*  active_element_input=element->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input);
    537                 active_element_input->GetInputValue(&active_element);
     624                element->GetInputValue(&active_element,HydrologydcMaskEplactiveEltEnum);
    538625                element->FindParam(&dt,TimesteppingTimeStepEnum);
    539626
    … …  
    590677                        }
    591678                }
    592                 element->AddInput(HydrologydcEplThicknessSubstepEnum,thickness,element->GetElementType());
     679                element->AddInput2(HydrologydcEplThicknessSubstepEnum,thickness,element->GetElementType());
    593680                xDelete<IssmDouble>(thickness);
    594681                xDelete<IssmDouble>(eplhead);
    … …  
    602689        }
    603690}/*}}}*/
    604 
    605 void HydrologyDCEfficientAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    606         /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
    607          * For node i, Bi can be expressed in the actual coordinate system
    608          * by:
    609          *       Bi=[ dN/dx ]
    610          *          [ dN/dy ]
    611          * where N is the finiteelement function for node i.
    612          *
    613          * We assume B has been allocated already, of size: 3x(NDOF2*numnodes)
    614          */
    615 
    616         /*Fetch number of nodes for this finite element*/
    617         int numnodes = element->GetNumberOfNodes();
    618 
    619         /*Get nodal functions derivatives*/
    620         IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
    621         element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    622 
    623         /*Build B: */
    624         for(int i=0;i<numnodes;i++){
    625                 B[numnodes*0+i] = dbasis[0*numnodes+i];
    626                 B[numnodes*1+i] = dbasis[1*numnodes+i];
    627         }
    628 
    629         /*Clean-up*/
    630         xDelete<IssmDouble>(dbasis);
    631 }/*}}}*/
    632 
    633 void  HydrologyDCEfficientAnalysis::HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask, Vector<IssmDouble>* recurence, Element* element){
     691void  HydrologyDCEfficientAnalysis::HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask, Vector<IssmDouble>* recurence, Element* element){/*{{{*/
    634692
    635693        bool        active_element;
    … …  
    663721        IssmDouble colapse_thick =basalelement->FindParam(HydrologydcEplColapseThicknessEnum);
    664722
    665         Input* active_element_input=basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input);
    666         active_element_input->GetInputValue(&active_element);
     723        basalelement->GetInputValue(&active_element,HydrologydcMaskEplactiveEltEnum);
    667724
    668725        basalelement-> GetInputListOnVertices(&old_active[0],HydrologydcMaskEplactiveNodeEnum);
    … …  
    712769                }
    713770        }
    714         basalelement->AddInput(HydrologydcEplThicknessSubstepEnum,epl_thickness,basalelement->GetElementType());
     771        element->AddBasalInput2(HydrologydcEplThicknessSubstepEnum,epl_thickness,basalelement->GetElementType());
    715772
    716773        if(domaintype!=Domain2DhorizontalEnum){
    … …  
    726783}
    727784/*}}}*/
    728 IssmDouble HydrologyDCEfficientAnalysis::EplStoring(Element* element,Gauss* gauss, Input* epl_thick_input, Input* epl_head_input, Input* base_input){/*{{{*/
    729         IssmDouble epl_storing;
    730         IssmDouble water_sheet,storing;
    731         IssmDouble epl_thickness,prestep_head,base_elev;
    732         IssmDouble rho_freshwater        = element->FindParam(MaterialsRhoFreshwaterEnum);
    733         IssmDouble g                     = element->FindParam(ConstantsGEnum);
    734         IssmDouble epl_porosity                                  = element->FindParam(HydrologydcEplPorosityEnum);
    735         IssmDouble epl_compressibility   = element->FindParam(HydrologydcEplCompressibilityEnum);
    736         IssmDouble water_compressibility = element->FindParam(HydrologydcWaterCompressibilityEnum);
    737 
    738         epl_thick_input->GetInputValue(&epl_thickness,gauss);
    739         epl_head_input->GetInputValue(&prestep_head,gauss);
    740         base_input->GetInputValue(&base_elev,gauss);
    741         water_sheet=max(0.0,(prestep_head-base_elev));
    742         storing=rho_freshwater*g*epl_porosity*epl_thickness*(water_compressibility+(epl_compressibility/epl_porosity));
    743 
    744         /* //porosity for unconfined region */
    745         /* if (water_sheet<=0.9*epl_thickness){ */
    746         /*      epl_storing=epl_porosity; */
    747         /* } */
    748         /* //continuity ramp */
    749         /* else if((water_sheet<epl_thickness) && (water_sheet>0.9*epl_thickness)){ */
    750         /*      epl_storing=(epl_thickness-water_sheet)*(epl_porosity-storing)/(0.1*epl_thickness)+storing; */
    751         /* } */
    752         /* //storing coefficient for confined */
    753         /* else{ */
    754         /*      epl_storing=storing; */
    755         /* } */
    756         /* return epl_storing; */
    757         return storing;
    758 }/*}}}*/
    759 
    760 IssmDouble HydrologyDCEfficientAnalysis::EplTransmitivity(Element* element,Gauss* gauss, Input* epl_thick_input, Input* epl_head_input, Input* base_input){/*{{{*/
    761         IssmDouble epl_transmitivity;
    762         IssmDouble water_sheet;
    763         IssmDouble epl_thickness,base_elev,prestep_head;
    764         IssmDouble epl_conductivity      = element->FindParam(HydrologydcEplConductivityEnum);
    765         epl_thick_input->GetInputValue(&epl_thickness,gauss);
    766         epl_head_input->GetInputValue(&prestep_head,gauss);
    767         base_input->GetInputValue(&base_elev,gauss);
    768 
    769         water_sheet=max(0.0,(prestep_head-base_elev));
    770         epl_transmitivity=epl_conductivity*epl_thickness;
    771         //epl_transmitivity=max(1.0e-6,(epl_conductivity*min(water_sheet,epl_thickness)));
    772         return epl_transmitivity;
    773 }/*}}}*/
    774 
    775785void HydrologyDCEfficientAnalysis::HydrologyEPLGetActive(Vector<IssmDouble>* active_vec, Element* element){/*{{{*/
    776786        /*Constants*/
    … …  
    798808        /*Pass the activity mask from elements to nodes*/
    799809        basalelement->GetInputListOnVertices(&active[0],HydrologydcMaskEplactiveNodeEnum);
    800         Input*  active_element_input=basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input);
    801         active_element_input->GetInputValue(&active_element);
     810        basalelement->GetInputValue(&active_element,HydrologydcMaskEplactiveEltEnum);
    802811
    803812        for(int i=0;i<numnodes;i++) flag+=active[i];
    … …  
    822831}
    823832/*}}}*/
    824 
    825 void HydrologyDCEfficientAnalysis::GetHydrologyDCInefficientHmax(IssmDouble* ph_max,Element* element, Node* innode){/*{{{*/
    826 
    827         int        hmax_flag;
    828         IssmDouble h_max;
    829         IssmDouble rho_ice,rho_water;
    830         IssmDouble thickness,bed;
    831         /*Get the flag to the limitation method*/
    832         element->FindParam(&hmax_flag,HydrologydcSedimentlimitFlagEnum);
    833 
    834         /*Switch between the different cases*/
    835         switch(hmax_flag){
    836         case 0:
    837                 h_max=1.0e+10;
    838                 break;
    839         case 1:
    840                 element->FindParam(&h_max,HydrologydcSedimentlimitEnum);
    841                 break;
    842         case 2:
    843                 /*Compute max*/
    844                 rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
    845                 rho_ice   = element->FindParam(MaterialsRhoIceEnum);
    846                 element-> GetInputValue(&thickness,innode,ThicknessEnum);
    847                 element-> GetInputValue(&bed,innode,BaseEnum);
    848                 h_max=((rho_ice*thickness)/rho_water)+bed;
    849                 break;
    850         case 3:
    851                 _error_("Using normal stress  not supported yet");
    852                 break;
    853         default:
    854                 _error_("no case higher than 3 for SedimentlimitFlag");
    855         }
    856         /*Assign output pointer*/
    857         *ph_max=h_max;
    858 }
    859 /*}}}*/
Note: See TracChangeset for help on using the changeset viewer.