Index: ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 19090) +++ ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 19091) @@ -645,7 +645,7 @@ xDelete(dbasis); }/*}}}*/ -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; int i,j; @@ -671,6 +671,8 @@ int eplflip_lock; int numnodes =basalelement->GetNumberOfNodes(); + IssmDouble new_active; + IssmDouble recurent; IssmDouble* epl_thickness =xNew(numnodes); IssmDouble* old_active =xNew(numnodes); IssmDouble* sedhead =xNew(numnodes); @@ -696,44 +698,55 @@ for(i=1;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)*/ GetHydrologyDCInefficientHmax(&h_max,basalelement,basalelement->nodes[i]); if(eplhead[i]>=h_max && active_element){ for(j=0;jnodes[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; + } } basalelement->AddInput(HydrologydcEplThicknessEnum,epl_thickness,basalelement->GetElementType()); Index: ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h =================================================================== --- ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h (revision 19090) +++ ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h (revision 19091) @@ -39,7 +39,7 @@ IssmDouble SedimentStoring(Element* element); 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); void ComputeEPLThickness(FemModel* femmodel); Index: ../trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp =================================================================== --- ../trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 19090) +++ ../trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 19091) @@ -91,7 +91,7 @@ femmodel->UpdateConstraintsx(); 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); /*{{{*//*Treating the sediment*/ Index: ../trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- ../trunk-jpl/src/c/classes/FemModel.cpp (revision 19090) +++ ../trunk-jpl/src/c/classes/FemModel.cpp (revision 19091) @@ -1912,6 +1912,7 @@ void FemModel::HydrologyEPLupdateDomainx(IssmDouble* pEplcount){ /*{{{*/ Vector* mask = NULL; + Vector* recurence = NULL; IssmDouble* serial_mask = NULL; Vector* active = NULL; IssmDouble* serial_active = NULL; @@ -1921,11 +1922,11 @@ HydrologyDCInefficientAnalysis* inefanalysis = new HydrologyDCInefficientAnalysis(); /*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); /*Assemble and serialize*/ @@ -1933,6 +1934,7 @@ serial_mask=mask->ToMPISerial(); xDelete(eplzigzag_counter); delete mask; + delete recurence; /*Update Mask*/ InputUpdateFromVectorx(this,serial_mask,HydrologydcMaskEplactiveNodeEnum,NodeSIdEnum);