Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 17371) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 17372) @@ -367,36 +367,4 @@ return rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity)); }/*}}}*/ -IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyDCInefficientHmax(Element* element, Gauss* gauss){/*{{{*/ - int hmax_flag; - IssmDouble h_max; - IssmDouble rho_ice,rho_water; - IssmDouble thickness,bed; - - /*Get the flag to the limitation method*/ - element->FindParam(&hmax_flag,HydrologydcSedimentlimitFlagEnum); - - /*Switch between the different cases*/ - switch(hmax_flag){ - case 0: - h_max=1.0e+10; - break; - case 1: - element->FindParam(&h_max,HydrologydcSedimentlimitEnum); - break; - case 2: - rho_water= element->GetMaterialParameter(MaterialsRhoFreshwaterEnum); - rho_ice= element->GetMaterialParameter(MaterialsRhoIceEnum); - element->GetInputValue(&thickness,gauss,ThicknessEnum); - element->GetInputValue(&bed,gauss,BedEnum); - h_max=((rho_ice*thickness)/rho_water)+bed; - break; - case 3: - _error_("Using normal stress not supported yet"); - break; - default: - _error_("no case higher than 3 for SedimentlimitFlag"); - } - return h_max; -}/*}}}*/ IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss){/*{{{*/ @@ -422,5 +390,5 @@ /* get input */ - element->GetInputValue(&epl_thickness,gauss,HydrologydcEplThicknessEnum); + element->GetInputValue(&epl_thickness,gauss,HydrologydcEplThicknessEnum); element->GetInputValue(&sed_head,gauss,SedimentHeadEnum); element->GetInputValue(&epl_head,gauss,EplHeadEnum); @@ -519,10 +487,12 @@ } - int numnodes = element->GetNumberOfNodes(); - IssmDouble thickness[numnodes]; - IssmDouble eplhead[numnodes]; - IssmDouble epl_slopeX[numnodes],epl_slopeY[numnodes]; - IssmDouble old_thickness[numnodes]; - IssmDouble ice_thickness[numnodes],bed[numnodes]; + int numnodes = element->GetNumberOfNodes(); + IssmDouble* thickness = xNew(numnodes); + IssmDouble* eplhead = xNew(numnodes); + IssmDouble* epl_slopeX = xNew(numnodes); + IssmDouble* epl_slopeY = xNew(numnodes); + IssmDouble* old_thickness = xNew(numnodes); + IssmDouble* ice_thickness = xNew(numnodes); + IssmDouble* bed = xNew(numnodes); Input* active_element_input=element->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input); @@ -574,4 +544,11 @@ } element->AddInput(HydrologydcEplThicknessEnum,thickness,P1Enum); + xDelete(thickness); + xDelete(eplhead); + xDelete(epl_slopeX); + xDelete(epl_slopeY); + xDelete(old_thickness); + xDelete(ice_thickness); + xDelete(bed); } } @@ -604,2 +581,159 @@ xDelete(dbasis); }/*}}}*/ +void HydrologyDCEfficientAnalysis::HydrologyEPLGetMask(Vector* vec_mask,Element* element){ + + bool active_element; + int i,j; + int meshtype; + IssmDouble h_max; + IssmDouble sedheadmin; + Element* basalelement=NULL; + + /*Get basal element*/ + element->FindParam(&meshtype,MeshTypeEnum); + switch(meshtype){ + case Mesh2DhorizontalEnum: + basalelement = element; + break; + case Mesh3DEnum: + if(!element->IsOnBed()) return; + basalelement = element->SpawnBasalElement(); + break; + default: _error_("mesh "<GetNumberOfNodes(); + IssmDouble* epl_thickness =xNew(numnodes); + IssmDouble* old_active =xNew(numnodes); + IssmDouble* sedhead =xNew(numnodes); + IssmDouble* eplhead =xNew(numnodes); + IssmDouble* residual =xNew(numnodes); + + IssmDouble init_thick = basalelement->GetMaterialParameter(HydrologydcEplInitialThicknessEnum); + + Input* active_element_input=basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input); + active_element_input->GetInputValue(&active_element); + + basalelement-> GetInputListOnVertices(&old_active[0],HydrologydcMaskEplactiveNodeEnum); + basalelement-> GetInputListOnVertices(&epl_thickness[0],HydrologydcEplThicknessEnum); + basalelement-> GetInputListOnVertices(&sedhead[0],SedimentHeadEnum); + basalelement-> GetInputListOnVertices(&eplhead[0],EplHeadEnum); + basalelement-> GetInputListOnVertices(&residual[0],SedimentHeadResidualEnum); + + /*Get minimum sediment head of the element*/ + sedheadmin=sedhead[0]; + for(i=1;i0 */ + if(residual[i]>0.){ + vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL); + } + + /*If mask was already one, keep one*/ + else if(old_active[i]>0.){ + vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL); + /*If epl thickness gets under , close the layer*/ + if(epl_thickness[i]<0.001*init_thick){ + vec_mask->SetValue(basalelement->nodes[i]->Sid(),0.,INS_VAL); + epl_thickness[i]=init_thick; + } + } + /*Increase of the efficient system is needed if the epl head reach the maximum value (sediment max value for now)*/ + GetHydrologyDCInefficientHmax(&h_max,basalelement,basalelement->nodes[i]); + if(eplhead[i]>=h_max && active_element){ + for(j=0;j0.){ + vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL); + } + /*Increase of the domain is on the downstream node in term of sediment head*/ + if(sedhead[j] == sedheadmin){ + vec_mask->SetValue(basalelement->nodes[j]->Sid(),1.,INS_VAL); + } + } + } + } + if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; + xDelete(epl_thickness); + xDelete(old_active); + xDelete(sedhead); + xDelete(eplhead); + xDelete(residual); +} +/*}}}*/ +void HydrologyDCEfficientAnalysis::HydrologyEPLGetActive(Vector* active_vec, Element* element){ + /*Constants*/ + + int meshtype; + Element* basalelement=NULL; + + /*Get basal element*/ + element->FindParam(&meshtype,MeshTypeEnum); + switch(meshtype){ + case Mesh2DhorizontalEnum: + basalelement = element; + break; + case Mesh3DEnum: + if(!element->IsOnBed()) return; + basalelement = element->SpawnBasalElement(); + break; + default: _error_("mesh "<GetNumberOfNodes(); + IssmDouble flag = 0.; + IssmDouble* active = xNew(numnodes); + + basalelement->GetInputListOnVertices(&active[0],HydrologydcMaskEplactiveNodeEnum); + + for(int i=0;i0.){ + for(int i=0;iSetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL); + } + } + else{ + /*Do not do anything: at least one node is active for this element but this element is not solved for*/ + } + if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; + xDelete(active); +} + +void HydrologyDCEfficientAnalysis::GetHydrologyDCInefficientHmax(IssmDouble* ph_max,Element* element, Node* innode){/*{{{*/ + + int hmax_flag; + IssmDouble h_max; + IssmDouble rho_ice,rho_water; + IssmDouble thickness,bed; + /*Get the flag to the limitation method*/ + element->FindParam(&hmax_flag,HydrologydcSedimentlimitFlagEnum); + + /*Switch between the different cases*/ + switch(hmax_flag){ + case 0: + h_max=1.0e+10; + break; + case 1: + element->FindParam(&h_max,HydrologydcSedimentlimitEnum); + break; + case 2: + /*Compute max*/ + rho_water = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum); + rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum); + element->GetInputValue(&thickness,innode,ThicknessEnum); + element->GetInputValue(&bed,innode,BedEnum); + h_max=((rho_ice*thickness)/rho_water)+bed; + break; + case 3: + _error_("Using normal stress not supported yet"); + break; + default: + _error_("no case higher than 3 for SedimentlimitFlag"); + } + /*Assign output pointer*/ + *ph_max=h_max; +} +/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h (revision 17371) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h (revision 17372) @@ -8,5 +8,5 @@ /*Headers*/ #include "./Analysis.h" - +class Node; class HydrologyDCEfficientAnalysis: public Analysis{ @@ -34,7 +34,9 @@ IssmDouble EplSpecificStoring(Element* element); IssmDouble SedimentStoring(Element* element); - IssmDouble GetHydrologyDCInefficientHmax(Element* element, Gauss* gauss); IssmDouble GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss); IssmDouble GetHydrologyPVectorTransfer(Element* element, Gauss* gauss); + void HydrologyEPLGetMask(Vector* vec_mask,Element* element); + void HydrologyEPLGetActive(Vector* active_vec, Element* element); + void GetHydrologyDCInefficientHmax(IssmDouble* ph_max,Element* element, Node* innode); void ComputeEPLThickness(FemModel* femmodel); }; Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 17371) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 17372) @@ -2,4 +2,5 @@ #include "../toolkits/toolkits.h" #include "../classes/classes.h" +#include "../classes/Node.h" #include "../shared/shared.h" #include "../modules/modules.h" @@ -419,5 +420,5 @@ for(int i=0;iGetHydrologyDCInefficientHmax(&h_max,basalelement->GetNode(i)); + GetHydrologyDCInefficientHmax(&h_max,basalelement,basalelement->GetNode(i)); if(values[i]>h_max) residual[i] = kappa*(values[i]-h_max); else residual[i] = 0.; @@ -496,5 +497,39 @@ return h_max; }/*}}}*/ - +void HydrologyDCInefficientAnalysis::GetHydrologyDCInefficientHmax(IssmDouble* ph_max,Element* element, Node* innode){/*{{{*/ + + int hmax_flag; + IssmDouble h_max; + IssmDouble rho_ice,rho_water; + IssmDouble thickness,bed; + /*Get the flag to the limitation method*/ + element->FindParam(&hmax_flag,HydrologydcSedimentlimitFlagEnum); + + /*Switch between the different cases*/ + switch(hmax_flag){ + case 0: + h_max=1.0e+10; + break; + case 1: + element->FindParam(&h_max,HydrologydcSedimentlimitEnum); + break; + case 2: + /*Compute max*/ + rho_water = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum); + rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum); + element->GetInputValue(&thickness,innode,ThicknessEnum); + element->GetInputValue(&bed,innode,BedEnum); + h_max=((rho_ice*thickness)/rho_water)+bed; + break; + case 3: + _error_("Using normal stress not supported yet"); + break; + default: + _error_("no case higher than 3 for SedimentlimitFlag"); + } + /*Assign output pointer*/ + *ph_max=h_max; +} +/*}}}*/ IssmDouble HydrologyDCInefficientAnalysis::GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss){/*{{{*/ Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.h (revision 17371) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.h (revision 17372) @@ -8,5 +8,5 @@ /*Headers*/ #include "./Analysis.h" - +class Node; class HydrologyDCInefficientAnalysis: public Analysis{ @@ -35,4 +35,5 @@ IssmDouble EplSpecificStoring(Element* element); IssmDouble GetHydrologyDCInefficientHmax(Element* element, Gauss* gauss); + void GetHydrologyDCInefficientHmax(IssmDouble* ph_max,Element* element, Node* innode); IssmDouble GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss); IssmDouble GetHydrologyPVectorTransfer(Element* element, Gauss* gauss); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 17371) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 17372) @@ -283,9 +283,4 @@ virtual void UpdateConstraintsExtrudeFromTop(void)=0; - virtual void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode)=0; - virtual void HydrologyEPLGetMask(Vector* mask)=0; - virtual void HydrologyEPLGetActive(Vector* active)=0; - // virtual void ComputeEPLThickness(void)=0; - virtual void MigrateGroundingLine(IssmDouble* sheet_ungrounding)=0; virtual void PotentialUngrounding(Vector* potential_sheet_ungrounding)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 17371) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 17372) @@ -4155,15 +4155,4 @@ /*}}}*/ #endif - -/*FUNCTION Penta::GetHydrologyDCInefficientHmax{{{*/ -void Penta::GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode){ - - if (!IsOnBed()) return; - - Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1. - tria->GetHydrologyDCInefficientHmax(ph_max,innode); - delete tria->material; delete tria; -} -/*}}}*/ /*FUNCTION Penta::GetSolutionFromInputsOneDof {{{*/ void Penta::GetSolutionFromInputsOneDof(Vector* solution, int enum_type){ @@ -4197,104 +4186,4 @@ } /*}}}*/ -/*FUNCTION Penta::HydrologyEPLGetActive {{{*/ -void Penta::HydrologyEPLGetActive(Vector* active_vec){ - - if (!IsOnBed()){ - return; - } - Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1. - tria->HydrologyEPLGetActive(active_vec); - delete tria->material; delete tria; - -} -/*}}}*/ -/*FUNCTION Penta::HydrologyEPLGetMask{{{*/ -void Penta::HydrologyEPLGetMask(Vector* vec_mask){ - - if (!IsOnBed())return; - - Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1. - tria->HydrologyEPLGetMask(vec_mask); - delete tria->material; delete tria; - -} -/*}}}*/ -/* /\*FUNCTION Penta::ComputeEPLThickness{{{*\/ */ -/* void Penta::ComputeEPLThickness(void){ */ - -/* int i; */ -/* const int numdof = NDOF1 *NUMVERTICES; */ -/* const int numdof2d = NDOF1 *NUMVERTICES2D; */ -/* bool isefficientlayer; */ -/* IssmDouble n,A,dt,init_thick; */ -/* IssmDouble rho_water,rho_ice; */ -/* IssmDouble gravity,latentheat,EPLgrad; */ -/* IssmDouble EPL_N,epl_conductivity; */ -/* IssmDouble activeEpl[numdof],thickness[numdof]; */ -/* IssmDouble eplhead[numdof], old_thickness[numdof]; */ -/* IssmDouble epl_slopeX[numdof],epl_slopeY[numdof]; */ -/* IssmDouble ice_thickness[numdof],bed[numdof]; */ -/* Penta *penta = NULL; */ -/* /\*If not on bed, return*\/ */ -/* if (!IsOnBed())return; */ - -/* /\*Get the flag to know if the efficient layer is present*\/ */ -/* this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); */ -/* this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); */ - -/* if(isefficientlayer){ */ -/* /\*For now, assuming just one way to compute EPL thickness*\/ */ -/* rho_water = matpar->GetRhoWater(); */ -/* rho_ice = matpar->GetRhoIce(); */ -/* gravity = matpar->GetG(); */ -/* latentheat = matpar->GetLatentHeat(); */ -/* epl_conductivity = matpar->GetEplConductivity(); */ -/* init_thick = matpar->GetEplInitialThickness(); */ -/* n = material->GetN(); */ -/* A = material->GetA(); */ - -/* GetInputListOnVertices(&activeEpl[0],HydrologydcMaskEplactiveNodeEnum); */ -/* GetInputListOnVertices(&eplhead[0],EplHeadEnum); */ -/* GetInputListOnVertices(&epl_slopeX[0],EplHeadSlopeXEnum); */ -/* GetInputListOnVertices(&epl_slopeY[0],EplHeadSlopeYEnum); */ -/* GetInputListOnVertices(&old_thickness[0],HydrologydcEplThicknessOldEnum); */ -/* GetInputListOnVertices(&ice_thickness[0],ThicknessEnum); */ -/* GetInputListOnVertices(&bed[0],BedEnum); */ - -/* for(int i=0;iinputs->AddInput(new PentaInput(HydrologydcEplThicknessEnum,thickness,P1Enum)); */ - -/* /\*Stop if we have reached the surface*\/ */ -/* if (penta->IsOnSurface()) break; */ - -/* /\* get upper Penta*\/ */ -/* penta=penta->GetUpperPenta(); _assert_(penta->Id()!=this->id); */ -/* } */ -/* } */ -/* } */ -/* /\*}}}*\/ */ - /*FUNCTION Penta::MigrateGroundingLine{{{*/ void Penta::MigrateGroundingLine(IssmDouble* phi_ungrounding){ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 17371) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 17372) @@ -241,9 +241,4 @@ IssmDouble StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa); - void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode); - void HydrologyEPLGetActive(Vector* active_vec); - void HydrologyEPLGetMask(Vector* vec_mask); - // void ComputeEPLThickness(void); - void UpdateConstraintsExtrudeFromBase(void); void UpdateConstraintsExtrudeFromTop(void); Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 17371) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 17372) @@ -150,7 +150,7 @@ void GetNormalFromLSF(IssmDouble *pnormal){_error_("not implemented yet");}; - void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode){_error_("not implemented yet");}; - void HydrologyEPLGetActive(Vector* active_vec){_error_("not implemented yet");}; - void HydrologyEPLGetMask(Vector* vec_mask){_error_("not implemented yet");}; + //void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode){_error_("not implemented yet");}; + //void HydrologyEPLGetActive(Vector* active_vec){_error_("not implemented yet");}; + //void HydrologyEPLGetMask(Vector* vec_mask){_error_("not implemented yet");}; // void ComputeEPLThickness(void){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17371) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17372) @@ -4425,197 +4425,4 @@ } /*}}}*/ -/*FUNCTION Tria::GetHydrologyDCInefficientHmax{{{*/ -void Tria::GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode){ - - int hmax_flag; - IssmDouble h_max; - IssmDouble rho_ice,rho_water; - IssmDouble thickness,bed; - /*Get the flag to the limitation method*/ - this->parameters->FindParam(&hmax_flag,HydrologydcSedimentlimitFlagEnum); - - /*Switch between the different cases*/ - switch(hmax_flag){ - case 0: - h_max=1.0e+10; - break; - case 1: - parameters->FindParam(&h_max,HydrologydcSedimentlimitEnum); - break; - case 2: - rho_ice=matpar->GetRhoIce(); - rho_water=matpar->GetRhoFreshwater(); - this->GetInputValue(&thickness,innode,ThicknessEnum); - this->GetInputValue(&bed,innode,BedEnum); - h_max=((rho_ice*thickness)/rho_water)+bed; - break; - case 3: - _error_("Using normal stress not supported yet"); - break; - default: - _error_("no case higher than 3 for SedimentlimitFlag"); - } - /*Assign output pointer*/ - *ph_max=h_max; -} -/*}}}*/ -/*FUNCTION Tria::HydrologyEPLGetActive {{{*/ -void Tria::HydrologyEPLGetActive(Vector* active_vec){ - - /*Constants*/ - const int numnodes = NUMVERTICES; - IssmDouble flag = 0.; - IssmDouble active[numnodes]; - - GetInputListOnVertices(&active[0],HydrologydcMaskEplactiveNodeEnum); - - for(int i=0;i0.){ - for(int i=0;iSetValue(nodes[i]->Sid(),1.,INS_VAL); - } - } - else{ - /*Do not do anything: at least one node is active for this element but this element is not solved for*/ - } -} -/*}}}*/ -/*FUNCTION Tria::HydrologyEPLGetMask{{{*/ -void Tria::HydrologyEPLGetMask(Vector* vec_mask){ - - /*Intermediaries*/ - int i,j; - const int numdof = NDOF1 *NUMVERTICES; - IssmDouble init_thick; - IssmDouble h_max; - IssmDouble sedheadmin; - IssmDouble epl_thickness[numdof]; - IssmDouble old_active[numdof]; - IssmDouble sedhead[numdof]; - IssmDouble eplhead[numdof]; - IssmDouble residual[numdof]; - - bool active_element; - Input* active_element_input=NULL; - - init_thick = matpar->GetEplInitialThickness(); - - active_element_input=inputs->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input); - active_element_input->GetInputValue(&active_element); - - GetInputListOnVertices(&old_active[0],HydrologydcMaskEplactiveNodeEnum); - GetInputListOnVertices(&epl_thickness[0],HydrologydcEplThicknessEnum); - GetInputListOnVertices(&sedhead[0],SedimentHeadEnum); - GetInputListOnVertices(&eplhead[0],EplHeadEnum); - GetInputListOnVertices(&residual[0],SedimentHeadResidualEnum); - - /*Get minimum sediment head of the element*/ - sedheadmin=sedhead[0]; - for(i=1;i0 */ - if(residual[i]>0.){ - vec_mask->SetValue(nodes[i]->Sid(),1.,INS_VAL); - } - - /*If mask was already one, keep one*/ - else if(old_active[i]>0.){ - vec_mask->SetValue(nodes[i]->Sid(),1.,INS_VAL); - /*If epl thickness gets under , close the layer*/ - if(epl_thickness[i]<0.001*init_thick){ - vec_mask->SetValue(nodes[i]->Sid(),0.,INS_VAL); - epl_thickness[i]=init_thick; - } - } - /*Increase of the efficient system is needed if the epl head reach the maximum value (sediment max value for now)*/ - this->GetHydrologyDCInefficientHmax(&h_max,this->nodes[i]); - if(eplhead[i]>=h_max && active_element){ - for(j=0;j0.){ - vec_mask->SetValue(nodes[i]->Sid(),1.,INS_VAL); - } - /*Increase of the domain is on the downstream node in term of sediment head*/ - if(sedhead[j] == sedheadmin){ - vec_mask->SetValue(nodes[j]->Sid(),1.,INS_VAL); - } - } - } - } -} -/*}}}*/ -//* *FUNCTION Tria::ComputeEPLThickness{{{*\/ */ -/* void Tria::ComputeEPLThickness(void){ */ - -/* const int numdof = NDOF1 *NUMVERTICES; */ -/* bool isefficientlayer; */ -/* bool active_element; */ -/* IssmDouble n,A,dt,init_thick; */ -/* IssmDouble rho_water,rho_ice; */ -/* IssmDouble gravity,latentheat,EPLgrad2; */ -/* IssmDouble EPL_N,epl_conductivity; */ -/* IssmDouble thickness[numdof]; */ -/* IssmDouble eplhead[numdof]; */ -/* IssmDouble epl_slopeX[numdof],epl_slopeY[numdof]; */ -/* IssmDouble old_thickness[numdof]; */ -/* IssmDouble ice_thickness[numdof],bed[numdof]; */ - -/* Input* active_element_input=NULL; */ - -/* /\*Get the flag to know if the efficient layer is present*\/ */ -/* this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); */ -/* this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); */ - -/* if(isefficientlayer){ */ - -/* /\*For now, assuming just one way to compute EPL thickness*\/ */ -/* rho_water = matpar->GetRhoWater(); */ -/* rho_ice = matpar->GetRhoIce(); */ -/* gravity = matpar->GetG(); */ -/* latentheat = matpar->GetLatentHeat(); */ -/* epl_conductivity = matpar->GetEplConductivity(); */ -/* init_thick = matpar->GetEplInitialThickness(); */ -/* n = material->GetN(); */ -/* A = material->GetAbar(); */ - -/* active_element_input=inputs->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input); */ -/* active_element_input->GetInputValue(&active_element); */ - -/* GetInputListOnVertices(&eplhead[0],EplHeadEnum); */ -/* GetInputListOnVertices(&epl_slopeX[0],EplHeadSlopeXEnum); */ -/* GetInputListOnVertices(&epl_slopeY[0],EplHeadSlopeYEnum); */ -/* GetInputListOnVertices(&old_thickness[0],HydrologydcEplThicknessOldEnum); */ -/* GetInputListOnVertices(&ice_thickness[0],ThicknessEnum); */ -/* GetInputListOnVertices(&bed[0],BedEnum); */ - -/* if(!active_element){ */ -/* /\*Keeping thickness to initial value if EPL is not active*\/ */ -/* for(int i=0;i10.0*init_thick){ */ -/* thickness[i] = 10.0*init_thick; */ -/* } */ -/* } */ -/* } */ -/* this->inputs->AddInput(new TriaInput(HydrologydcEplThicknessEnum,thickness,P1Enum)); */ -/* } */ -/* } */ -/*}}}*/ #ifdef _HAVE_DAKOTA_ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 17371) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 17372) @@ -253,9 +253,4 @@ void CreateHydrologyWaterVelocityInput(void); - - void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode); - void HydrologyEPLGetActive(Vector* active_vec); - void HydrologyEPLGetMask(Vector* vec_mask); - // void ComputeEPLThickness(void); /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 17371) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 17372) @@ -15,4 +15,6 @@ #include "./modules/modules.h" #include "../shared/Enum/Enum.h" + +#include "../analyses/analyses.h" /*module includes: {{{*/ @@ -1352,10 +1354,12 @@ Vector* active = NULL; IssmDouble* serial_active = NULL; - - /*Step 1: update maks, the mask might be extended by residual and/or using downstream sediment head*/ + + HydrologyDCEfficientAnalysis* effanalysis = new HydrologyDCEfficientAnalysis(); + /*Step 1: update mask, the mask might be extended by residual and/or using downstream sediment head*/ mask=new Vector(nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum)); + for (int i=0;iSize();i++){ Element* element=dynamic_cast(elements->GetObjectByOffset(i)); - element->HydrologyEPLGetMask(mask); + effanalysis->HydrologyEPLGetMask(mask,element); } @@ -1373,5 +1377,5 @@ for (int i=0;iSize();i++){ Element* element=dynamic_cast(elements->GetObjectByOffset(i)); - element->HydrologyEPLGetActive(active); + effanalysis->HydrologyEPLGetActive(active,element); } @@ -1396,4 +1400,5 @@ } xDelete(serial_active); + delete effanalysis; int sum_counter; ISSM_MPI_Reduce(&counter,&sum_counter,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,IssmComm::GetComm() ); @@ -1407,16 +1412,9 @@ } /*}}}*/ -/* void FemModel::HydrologyEPLThicknessx(void){ /\*{{{*\/ */ - -/* for (int i=0;iSize();i++){ */ -/* Element* element=dynamic_cast(elements->GetObjectByOffset(i)); */ -/* element->ComputeEPLThickness(); */ -/* } */ -/* } */ -/*}}}*/ void FemModel::UpdateConstraintsL2ProjectionEPLx(void){ /*{{{*/ Vector* active = NULL; IssmDouble* serial_active = NULL; + HydrologyDCEfficientAnalysis* effanalysis = new HydrologyDCEfficientAnalysis(); /*update node activity. If one element is connected to mask=1, all nodes are active*/ @@ -1424,5 +1422,5 @@ for (int i=0;iSize();i++){ Element* element=dynamic_cast(elements->GetObjectByOffset(i)); - element->HydrologyEPLGetActive(active); + effanalysis->HydrologyEPLGetActive(active,element); } @@ -1431,4 +1429,5 @@ serial_active=active->ToMPISerial(); delete active; + delete effanalysis; /*Update node activation accordingly*/ Index: /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp (revision 17371) +++ /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp (revision 17372) @@ -13,4 +13,5 @@ #include "../classes.h" #include "shared/shared.h" +#include "../../analyses/analyses.h" /*}}}*/ @@ -601,4 +602,5 @@ IssmDouble h; IssmDouble h_max; + HydrologyDCInefficientAnalysis* inefanalysis = NULL; /*check that pengrid is not a clone (penalty to be added only once)*/ @@ -616,6 +618,7 @@ /*Get sediment water head h*/ + inefanalysis = new HydrologyDCInefficientAnalysis(); element->GetInputValue(&h,node,SedimentHeadEnum); - element->GetHydrologyDCInefficientHmax(&h_max,node); + inefanalysis->GetHydrologyDCInefficientHmax(&h_max,element,node); parameters->FindParam(&penalty_lock,HydrologydcPenaltyLockEnum); @@ -646,4 +649,5 @@ /*Assign output pointers:*/ + delete inefanalysis; *punstable=unstable; } @@ -651,5 +655,4 @@ /*FUNCTION Pengrid::PenaltyCreateKMatrixHydrologyDCInefficient {{{*/ ElementMatrix* Pengrid::PenaltyCreateKMatrixHydrologyDCInefficient(IssmDouble kmax){ - IssmDouble penalty_factor; @@ -672,8 +675,10 @@ IssmDouble h_max; IssmDouble penalty_factor; + HydrologyDCInefficientAnalysis* inefanalysis = NULL; /*Initialize Element matrix and return if necessary*/ if(!this->active) return NULL; ElementVector* pe=new ElementVector(&node,1,this->parameters); + inefanalysis = new HydrologyDCInefficientAnalysis(); /*Retrieve parameters*/ @@ -681,8 +686,10 @@ /*Get h_max and compute penalty*/ - element->GetHydrologyDCInefficientHmax(&h_max,node); + inefanalysis->GetHydrologyDCInefficientHmax(&h_max,element,node); + pe->values[0]=kmax*pow(10.,penalty_factor)*h_max; /*Clean up and return*/ + delete inefanalysis; return pe; } Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 17371) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 17372) @@ -156,5 +156,4 @@ } /* }}} *//*End of the global sediment loop*/ - /* {{{ *//*Now dealing with the EPL in the same way*/ if(isefficientlayer){ @@ -170,5 +169,4 @@ ug_epl_sub_iter=ug_epl->Duplicate();_assert_(ug_epl_sub_iter); ug_epl->Copy(ug_epl_sub_iter); - /* {{{ *//*Retrieve the EPL head slopes and compute EPL Thickness*/ if(VerboseSolution()) _printf0_("computing EPL Head slope...\n"); Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_linear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_linear.cpp (revision 17371) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_linear.cpp (revision 17372) @@ -24,5 +24,4 @@ femmodel->UpdateConstraintsx(); SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel); - CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type); Reduceloadx(pf, Kfs, ys); delete Kfs;