Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 16599) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 16600) @@ -14,4 +14,5 @@ int sedimentlimit_flag; int transfer_flag; + int penalty_lock; bool isefficientlayer; IssmDouble sedimentlimit; @@ -31,4 +32,5 @@ iomodel->FetchData(&penalty_factor,HydrologydcPenaltyFactorEnum); iomodel->FetchData(&rel_tol,HydrologydcRelTolEnum); + iomodel->FetchData(&penalty_lock,HydrologydcPenaltyLockEnum); if(sedimentlimit_flag==1){ @@ -48,4 +50,5 @@ parameters->AddObject(new IntParam(HydrologydcTransferFlagEnum,transfer_flag)); parameters->AddObject(new DoubleParam(HydrologydcRelTolEnum,rel_tol)); + parameters->AddObject(new IntParam(HydrologydcPenaltyLockEnum,penalty_lock)); }/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16599) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16600) @@ -136,4 +136,5 @@ virtual void HydrologyEPLGetMask(Vector* mask)=0; virtual void HydrologyEPLGetActive(Vector* active)=0; + virtual void ComputeEPLThickness(void)=0; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16599) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16600) @@ -10845,4 +10845,15 @@ } /*}}}*/ +/*FUNCTION Penta::ComputeEPLThickness{{{*/ +void Penta::ComputeEPLThickness(void){ + + if (!IsOnBed())return; + + Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1. + tria->ComputeEPLThickness(); + delete tria->material; delete tria; + +} +/*}}}*/ /*FUNCTION Penta::InputUpdateFromSolutionHydrologyDCInefficient{{{*/ void Penta::InputUpdateFromSolutionHydrologyDCInefficient(IssmDouble* solution){ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16599) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16600) @@ -323,4 +323,5 @@ void HydrologyEPLGetMask(Vector* vec_mask); void InputUpdateFromSolutionHydrologyDCInefficient(IssmDouble* solution); + void ComputeEPLThickness(void); #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16599) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16600) @@ -101,4 +101,5 @@ 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");}; #endif void GetSolutionFromInputs(Vector* solution){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16599) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16600) @@ -7224,4 +7224,10 @@ } /*}}}*/ +/*FUNCTION Tria::ComputeEPLThickness{{{*/ +void Tria::ComputeEPLThickness(void){ +} +/*}}}*/ + + #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16599) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16600) @@ -303,4 +303,5 @@ void HydrologyEPLGetActive(Vector* active_vec); void HydrologyEPLGetMask(Vector* vec_mask); + void ComputeEPLThickness(void); bool AllActive(void); bool AnyActive(void); Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 16599) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 16600) @@ -1412,3 +1412,11 @@ } /*}}}*/ +void FemModel::HydrologyEPLThicknessx(void){ /*{{{*/ + + for (int i=0;iSize();i++){ + Element* element=dynamic_cast(elements->GetObjectByOffset(i)); + element->ComputeEPLThickness(); + } +} +/*}}}*/ #endif Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 16599) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 16600) @@ -102,4 +102,5 @@ void HydrologyTransferx(void); void HydrologyEPLupdateDomainx(void); + void HydrologyEPLThicknessx(void); #endif }; Index: /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp (revision 16599) +++ /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp (revision 16600) @@ -640,6 +640,7 @@ // The penalty is stable if it doesn't change during two consecutive iterations. - int unstable = 0; + int unstable=0; int new_active; + int penalty_lock; IssmDouble pressure; IssmDouble h; @@ -656,4 +657,6 @@ element->GetInputValue(&h,node,SedimentHeadEnum); element->GetHydrologyDCInefficientHmax(&h_max,node); + parameters->FindParam(&penalty_lock,HydrologydcPenaltyLockEnum); + if (h>h_max) new_active=1; @@ -661,9 +664,19 @@ new_active=0; - if(this->active==new_active) - unstable=0; - else - unstable=1; - + if(this->active==new_active){ + unstable=0; + } + else{ + unstable=1; + if(penalty_lock)zigzag_counter++; + } + + /*If penalty keeps zigzagging more than penalty_lock times: */ + if(penalty_lock){ + if(zigzag_counter>penalty_lock){ + unstable=0; + active=1; + } + } /*Set penalty flag*/ this->active=new_active; Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 16599) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 16600) @@ -121,4 +121,5 @@ HydrologydcLeakageFactorEnum, HydrologydcPenaltyFactorEnum, + HydrologydcPenaltyLockEnum, HydrologyLayerEnum, HydrologySedimentEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 16599) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 16600) @@ -129,4 +129,5 @@ case HydrologydcLeakageFactorEnum : return "HydrologydcLeakageFactor"; case HydrologydcPenaltyFactorEnum : return "HydrologydcPenaltyFactor"; + case HydrologydcPenaltyLockEnum : return "HydrologydcPenaltyLock"; case HydrologyLayerEnum : return "HydrologyLayer"; case HydrologySedimentEnum : return "HydrologySediment"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 16599) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 16600) @@ -129,4 +129,5 @@ else if (strcmp(name,"HydrologydcLeakageFactor")==0) return HydrologydcLeakageFactorEnum; else if (strcmp(name,"HydrologydcPenaltyFactor")==0) return HydrologydcPenaltyFactorEnum; + else if (strcmp(name,"HydrologydcPenaltyLock")==0) return HydrologydcPenaltyLockEnum; else if (strcmp(name,"HydrologyLayer")==0) return HydrologyLayerEnum; else if (strcmp(name,"HydrologySediment")==0) return HydrologySedimentEnum; @@ -136,9 +137,9 @@ else if (strcmp(name,"IndependentObject")==0) return IndependentObjectEnum; else if (strcmp(name,"InversionControlParameters")==0) return InversionControlParametersEnum; - else if (strcmp(name,"InversionCostFunctionThreshold")==0) return InversionCostFunctionThresholdEnum; else stage=2; } if(stage==2){ - if (strcmp(name,"InversionCostFunctionsCoefficients")==0) return InversionCostFunctionsCoefficientsEnum; + if (strcmp(name,"InversionCostFunctionThreshold")==0) return InversionCostFunctionThresholdEnum; + else if (strcmp(name,"InversionCostFunctionsCoefficients")==0) return InversionCostFunctionsCoefficientsEnum; else if (strcmp(name,"InversionCostFunctions")==0) return InversionCostFunctionsEnum; else if (strcmp(name,"InversionGradientScaling")==0) return InversionGradientScalingEnum; @@ -259,9 +260,9 @@ else if (strcmp(name,"SteadystateRequestedOutputs")==0) return SteadystateRequestedOutputsEnum; else if (strcmp(name,"Surface")==0) return SurfaceEnum; - else if (strcmp(name,"SurfaceforcingsPrecipitation")==0) return SurfaceforcingsPrecipitationEnum; else stage=3; } if(stage==3){ - if (strcmp(name,"SurfaceforcingsMassBalance")==0) return SurfaceforcingsMassBalanceEnum; + if (strcmp(name,"SurfaceforcingsPrecipitation")==0) return SurfaceforcingsPrecipitationEnum; + else if (strcmp(name,"SurfaceforcingsMassBalance")==0) return SurfaceforcingsMassBalanceEnum; else if (strcmp(name,"SurfaceforcingsIspdd")==0) return SurfaceforcingsIspddEnum; else if (strcmp(name,"SurfaceforcingsDesfac")==0) return SurfaceforcingsDesfacEnum; @@ -382,9 +383,9 @@ else if (strcmp(name,"ExternalResult")==0) return ExternalResultEnum; else if (strcmp(name,"FileParam")==0) return FileParamEnum; - else if (strcmp(name,"Input")==0) return InputEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"IntInput")==0) return IntInputEnum; + if (strcmp(name,"Input")==0) return InputEnum; + else if (strcmp(name,"IntInput")==0) return IntInputEnum; else if (strcmp(name,"InputToExtrude")==0) return InputToExtrudeEnum; else if (strcmp(name,"InputToL2Project")==0) return InputToL2ProjectEnum; @@ -505,9 +506,9 @@ else if (strcmp(name,"HydrologyWaterVx")==0) return HydrologyWaterVxEnum; else if (strcmp(name,"HydrologyWaterVy")==0) return HydrologyWaterVyEnum; - else if (strcmp(name,"StressTensor")==0) return StressTensorEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"StressTensorxx")==0) return StressTensorxxEnum; + if (strcmp(name,"StressTensor")==0) return StressTensorEnum; + else if (strcmp(name,"StressTensorxx")==0) return StressTensorxxEnum; else if (strcmp(name,"StressTensorxy")==0) return StressTensorxyEnum; else if (strcmp(name,"StressTensorxz")==0) return StressTensorxzEnum; Index: /issm/trunk-jpl/src/m/classes/hydrologydc.m =================================================================== --- /issm/trunk-jpl/src/m/classes/hydrologydc.m (revision 16599) +++ /issm/trunk-jpl/src/m/classes/hydrologydc.m (revision 16600) @@ -9,4 +9,5 @@ isefficientlayer = 0; penalty_factor = 0; + penalty_lock = 0; rel_tol = 0; sedimentlimit_flag = 0; @@ -107,4 +108,5 @@ fielddisplay(obj,'isefficientlayer','do we use an efficient drainage system [1: true; 0: false]'); fielddisplay(obj,'penalty_factor','exponent of the value used in the penalisation method [dimensionless]'); + fielddisplay(obj,'penalty_lock','stabilize unstable constraints that keep zigzagging after n iteration (default is 0, no stabilization)'); fielddisplay(obj,'rel_tol','tolerance of the nonlinear iteration for the transfer between layers [dimensionless]'); fielddisplay(obj,'sedimentlimit_flag','what kind of upper limit is applied for the inefficient layer'); @@ -147,4 +149,5 @@ WriteData(fid,'object',obj,'fieldname','isefficientlayer','format','Boolean'); WriteData(fid,'object',obj,'fieldname','penalty_factor','format','Double'); + WriteData(fid,'object',obj,'fieldname','penalty_lock','format','Integer'); WriteData(fid,'object',obj,'fieldname','rel_tol','format','Double'); WriteData(fid,'object',obj,'fieldname','sedimentlimit_flag','format','Integer'); Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.py =================================================================== --- /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 16599) +++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 16600) @@ -121,4 +121,5 @@ def HydrologydcLeakageFactorEnum(): return StringToEnum("HydrologydcLeakageFactor")[0] def HydrologydcPenaltyFactorEnum(): return StringToEnum("HydrologydcPenaltyFactor")[0] +def HydrologydcPenaltyLockEnum(): return StringToEnum("HydrologydcPenaltyLock")[0] def HydrologyLayerEnum(): return StringToEnum("HydrologyLayer")[0] def HydrologySedimentEnum(): return StringToEnum("HydrologySediment")[0] Index: /issm/trunk-jpl/src/m/enum/HydrologydcPenaltyLockEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/HydrologydcPenaltyLockEnum.m (revision 16600) +++ /issm/trunk-jpl/src/m/enum/HydrologydcPenaltyLockEnum.m (revision 16600) @@ -0,0 +1,11 @@ +function macro=HydrologydcPenaltyLockEnum() +%HYDROLOGYDCPENALTYLOCKENUM - Enum of HydrologydcPenaltyLock +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=HydrologydcPenaltyLockEnum() + +macro=StringToEnum('HydrologydcPenaltyLock');