Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 19090) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 19091) @@ -646,5 +646,5 @@ }/*}}}*/ -void HydrologyDCEfficientAnalysis::HydrologyEPLGetMask(Vector* vec_mask, int* eplzigzag_counter, Element* element){ +void HydrologyDCEfficientAnalysis::HydrologyEPLGetMask(Vector* vec_mask, Vector* recurence, int* eplzigzag_counter, Element* element){ bool active_element; @@ -672,4 +672,6 @@ int eplflip_lock; int numnodes =basalelement->GetNumberOfNodes(); + IssmDouble new_active; + IssmDouble recurent; IssmDouble* epl_thickness =xNew(numnodes); IssmDouble* old_active =xNew(numnodes); @@ -697,23 +699,29 @@ for(i=0;i0 */ + + /*Dealing with the initial value to define zigzaging, preceding iteration the first time we see the node and current evolving vector after*/ + recurence->GetValue(&recurent,basalelement->nodes[i]->Sid()); + if (recurent>0)vec_mask->GetValue(&old_active[i],basalelement->nodes[i]->Sid()); + new_active=old_active[i]; + recurence->SetValue(basalelement->nodes[i]->Sid(),1.,ADD_VAL); + + /*Now starting to look at the activations*/ if(residual[i]>0.){ - if(old_active[i]==0) eplzigzag_counter[basalelement->nodes[i]->Lid()] ++; vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL); - } - + new_active=1.0; + } /*If mask was already one, keep one or colapse*/ else if(old_active[i]>0.){ - vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL); + new_active=1.0; /*If epl thickness gets under colapse thickness, close the layer*/ if(epl_thickness[i]nodes[i]->Lid()] ++; - /*Avoid flipfloping between open and closed states*/ - if(eplzigzag_counter[basalelement->nodes[i]->Lid()]SetValue(basalelement->nodes[i]->Sid(),0.,INS_VAL); - epl_thickness[i]=init_thick; - } + vec_mask->SetValue(basalelement->nodes[i]->Sid(),0.,INS_VAL); + new_active=0.0; } + } + + if(basalelement->nodes[i]->Sid()==42){ + printf("old active is %e, and new active is %e, with counter= %i\n",old_active[i],new_active,eplzigzag_counter[basalelement->nodes[i]->Lid()]); } /*Increase of the efficient system is needed if the epl head reach the maximum value (sediment max value for now)*/ @@ -723,15 +731,20 @@ /*Increase of the domain is on the downstream node in term of sediment head*/ if(sedhead[j] == sedheadmin){ - if(old_active[j]==0) { - eplzigzag_counter[basalelement->nodes[j]->Lid()] ++; - if(eplzigzag_counter[basalelement->nodes[i]->Lid()]SetValue(basalelement->nodes[j]->Sid(),1.,INS_VAL); - } - else{ - vec_mask->SetValue(basalelement->nodes[j]->Sid(),0.,INS_VAL); - } - } + vec_mask->SetValue(basalelement->nodes[j]->Sid(),1.,INS_VAL); } } + } + /*Now checking for nodes zigzaging between open and close*/ + if (old_active[i] != new_active){ + eplzigzag_counter[basalelement->nodes[i]->Lid()] ++; + /*If node changed too much of state, fix it to it's last known state*/ + if(eplzigzag_counter[basalelement->nodes[i]->Lid()]>eplflip_lock & eplflip_lock!=0){ + vec_mask->SetValue(basalelement->nodes[i]->Sid(),old_active[i],INS_VAL); + new_active=old_active[i]; + } + } + /*If node is now closed bring its thickness back to initial*/ + if (new_active==0.){ + epl_thickness[i]=init_thick; } } Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h (revision 19090) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h (revision 19091) @@ -40,5 +40,5 @@ IssmDouble GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* thick_input, Input* base_input); IssmDouble GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* thick_input, Input* base_input); - void HydrologyEPLGetMask(Vector* vec_mask, int* eplzigzag_counter, Element* element); + void HydrologyEPLGetMask(Vector* vec_mask, Vector* recurence, int* eplzigzag_counter, Element* element); void HydrologyEPLGetActive(Vector* active_vec, Element* element); void GetHydrologyDCInefficientHmax(IssmDouble* ph_max,Element* element, Node* innode); Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 19090) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 19091) @@ -1913,4 +1913,5 @@ Vector* mask = NULL; + Vector* recurence = NULL; IssmDouble* serial_mask = NULL; Vector* active = NULL; @@ -1922,9 +1923,9 @@ /*Step 1: update mask, the mask might be extended by residual and/or using downstream sediment head*/ mask=new Vector(this->nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum)); + recurence=new Vector(this->nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum)); this->parameters->FindParam(&eplzigzag_counter,NULL,EplZigZagCounterEnum); - for (int i=0;iSize();i++){ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - effanalysis->HydrologyEPLGetMask(mask,eplzigzag_counter,element); + effanalysis->HydrologyEPLGetMask(mask,recurence,eplzigzag_counter,element); } this->parameters->SetParam(eplzigzag_counter,this->nodes->Size(),EplZigZagCounterEnum); @@ -1934,4 +1935,5 @@ xDelete(eplzigzag_counter); delete mask; + delete recurence; /*Update Mask*/ Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 19090) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 19091) @@ -92,5 +92,5 @@ femmodel->parameters->SetParam(HydrologySedimentEnum,HydrologyLayerEnum); - /*Reset constraint on the ZigZag Lock, this thing doesn't work, it have to disapear*/ + /*Reset constraint on the ZigZag Lock*/ ResetConstraintsx(femmodel);