Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 18718) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 18719) @@ -30,4 +30,5 @@ bool isefficientlayer; int hydrology_model; + /*Now, do we really want DC?*/ @@ -61,6 +62,6 @@ iomodel->FetchDataToInput(elements,MeshVertexonsurfaceEnum); } - -}/*}}}*/ +}/*}}}*/ + void HydrologyDCEfficientAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/ @@ -98,16 +99,41 @@ }/*}}}*/ +void HydrologyDCEfficientAnalysis::InitZigZagCounter(FemModel* femmodel){/*{{{*/ + int* eplzigzag_counter =NULL; + + eplzigzag_counter=xNewZeroInit(femmodel->nodes->Size()); + femmodel->parameters->AddObject(new IntVecParam(EplZigZagCounterEnum,eplzigzag_counter,femmodel->nodes->Size())); + xDelete(eplzigzag_counter); +}/*}}}*/ + +void HydrologyDCEfficientAnalysis::ResetCounter(FemModel* femmodel){/*{{{*/ + + int* eplzigzag_counter=NULL; + Element* element=NULL; + + femmodel->parameters->FindParam(&eplzigzag_counter,NULL,EplZigZagCounterEnum); + for(int i=0;inodes->Size();i++){ + + eplzigzag_counter[i]=0; + } + femmodel->parameters->SetParam(eplzigzag_counter,femmodel->nodes->Size(),EplZigZagCounterEnum); + xDelete(eplzigzag_counter); +}/*}}}*/ + /*Finite Element Analysis*/ void HydrologyDCEfficientAnalysis::Core(FemModel* femmodel){/*{{{*/ _error_("not implemented"); }/*}}}*/ -ElementVector* HydrologyDCEfficientAnalysis::CreateDVector(Element* element){/*{{{*/ + + ElementVector* HydrologyDCEfficientAnalysis::CreateDVector(Element* element){/*{{{*/ /*Default, return NULL*/ return NULL; }/*}}}*/ -ElementMatrix* HydrologyDCEfficientAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/ + + ElementMatrix* HydrologyDCEfficientAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/ _error_("Not implemented"); }/*}}}*/ -ElementMatrix* HydrologyDCEfficientAnalysis::CreateKMatrix(Element* element){/*{{{*/ + + ElementMatrix* HydrologyDCEfficientAnalysis::CreateKMatrix(Element* element){/*{{{*/ /*Intermediaries*/ @@ -612,5 +638,5 @@ xDelete(dbasis); }/*}}}*/ -void HydrologyDCEfficientAnalysis::HydrologyEPLGetMask(Vector* vec_mask,Element* element){ +void HydrologyDCEfficientAnalysis::HydrologyEPLGetMask(Vector* vec_mask, int* eplzigzag_counter, Element* element){ bool active_element; @@ -619,5 +645,5 @@ IssmDouble h_max; IssmDouble sedheadmin; - Element* basalelement=NULL; + Element* basalelement=NULL; /*Get basal element*/ @@ -636,4 +662,5 @@ /*Intermediaries*/ + int penalty_lock; int numnodes =basalelement->GetNumberOfNodes(); IssmDouble* epl_thickness =xNew(numnodes); @@ -642,5 +669,5 @@ IssmDouble* eplhead =xNew(numnodes); IssmDouble* residual =xNew(numnodes); - + IssmDouble init_thick =basalelement->GetMaterialParameter(HydrologydcEplInitialThicknessEnum); IssmDouble colapse_thick =basalelement->GetMaterialParameter(HydrologydcEplColapseThicknessEnum); @@ -648,4 +675,6 @@ Input* active_element_input=basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input); active_element_input->GetInputValue(&active_element); + + basalelement->parameters->FindParam(&penalty_lock,HydrologydcPenaltyLockEnum); basalelement-> GetInputListOnVertices(&old_active[0],HydrologydcMaskEplactiveNodeEnum); @@ -662,4 +691,5 @@ /*Activate EPL if residual is >0 */ 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); } @@ -670,6 +700,12 @@ /*If epl thickness gets under 10-3 initial thickness, close the layer*/ if(epl_thickness[i]SetValue(basalelement->nodes[i]->Sid(),0.,INS_VAL); - epl_thickness[i]=init_thick; + eplzigzag_counter[basalelement->nodes[i]->Lid()] ++; + if(eplzigzag_counter[basalelement->nodes[i]->Sid()]SetValue(basalelement->nodes[i]->Sid(),0.,INS_VAL); + epl_thickness[i]=init_thick; + } + else{ + vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL); + } } } @@ -678,5 +714,4 @@ if(eplhead[i]>=h_max && active_element){ - for(j=0;j0.){ @@ -685,10 +720,12 @@ /*Increase of the domain is on the downstream node in term of sediment head*/ if(sedhead[j] == sedheadmin){ + if(old_active[i]==0) eplzigzag_counter[basalelement->nodes[j]->Lid()] ++; vec_mask->SetValue(basalelement->nodes[j]->Sid(),1.,INS_VAL); } } } - - } + } + basalelement->AddInput(HydrologydcEplThicknessEnum,epl_thickness,basalelement->GetElementType()); + if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; xDelete(epl_thickness); @@ -723,4 +760,10 @@ basalelement->GetInputListOnVertices(&active[0],HydrologydcMaskEplactiveNodeEnum); + + bool active_element; + + Input* active_element_input=basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input); + active_element_input->GetInputValue(&active_element); + for(int i=0;iSetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL); } + } + else if(active_element){ + /*Checking Stuff*/ + for(int i=0;iSetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL); + } } else{ Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h (revision 18718) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h (revision 18719) @@ -20,5 +20,7 @@ void CreateConstraints(Constraints* constraints,IoModel* iomodel); void CreateLoads(Loads* loads, IoModel* iomodel); - + void InitZigZagCounter(FemModel* femmodel); + void ResetCounter(FemModel* femmodel); + /*Finite element Analysis*/ void Core(FemModel* femmodel); @@ -38,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,Element* element); + void HydrologyEPLGetMask(Vector* vec_mask, 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/analyses/HydrologyDCInefficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 18718) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 18719) @@ -17,5 +17,5 @@ int penalty_lock; int hydro_maxiter; - int* elementactive_counter =NULL; + /* int* elementactive_counter =NULL; */ bool isefficientlayer; IssmDouble sedimentlimit; @@ -56,10 +56,4 @@ parameters->AddObject(new IntParam(HydrologydcPenaltyLockEnum,penalty_lock)); parameters->AddObject(new IntParam(HydrologydcMaxIterEnum,hydro_maxiter)); - - if(isefficientlayer==1){ - elementactive_counter=xNewZeroInit(iomodel->numberofelements); - parameters->AddObject(new IntVecParam(ElementActiveCounterEnum,elementactive_counter,iomodel->numberofelements)); - xDelete(elementactive_counter); - } }/*}}}*/ @@ -104,16 +98,16 @@ iomodel->FetchDataToInput(elements,HydrologydcMaskEplactiveNodeEnum); - for(int i=0;iSize();i++){ - Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - Input* node_mask_input = element->GetInput(HydrologydcMaskEplactiveNodeEnum); _assert_(node_mask_input); + /* for(int i=0;iSize();i++){ */ + /* Element* element=xDynamicCast(elements->GetObjectByOffset(i)); */ + /* Input* node_mask_input = element->GetInput(HydrologydcMaskEplactiveNodeEnum); _assert_(node_mask_input); */ - if(node_mask_input->Max()>0.) { - element_active = true; - } - else{ - element_active = false; - } - element->AddInput(new BoolInput(HydrologydcMaskEplactiveEltEnum,element_active)); - } + /* if(node_mask_input->Max()>0.) { */ + /* element_active = true; */ + /* } */ + /* else{ */ + /* element_active = false; */ + /* } */ + /* element->AddInput(new BoolInput(HydrologydcMaskEplactiveEltEnum,element_active)); */ + /* } */ } }/*}}}*/ @@ -684,17 +678,9 @@ bool element_active; - bool old_active; - int penalty_lock; - int* elementactive_counter=NULL; Element* element=NULL; - femmodel->parameters->FindParam(&elementactive_counter,NULL,ElementActiveCounterEnum); for(int i=0;ielements->Size();i++){ element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); - element->parameters->FindParam(&penalty_lock,HydrologydcPenaltyLockEnum); - - Input* node_mask_input = element->GetInput(HydrologydcMaskEplactiveNodeEnum); _assert_(node_mask_input); - Input* old_active_input = element->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(old_active_input); - old_active_input->GetInputValue(&old_active); + Input* node_mask_input = element->GetInput(HydrologydcMaskEplactiveNodeEnum); _assert_(node_mask_input); if(node_mask_input->Max()>0.){ @@ -704,29 +690,5 @@ element_active = false; } - if(old_active!=element_active){ - elementactive_counter[i]=elementactive_counter[i]+1; - } - if(elementactive_counter[i]>penalty_lock){ - element_active = true; - } element->AddInput(new BoolInput(HydrologydcMaskEplactiveEltEnum,element_active)); } - femmodel->parameters->SetParam(elementactive_counter,femmodel->elements->Size(),ElementActiveCounterEnum); - xDelete(elementactive_counter); -}/*}}}*/ - -void HydrologyDCInefficientAnalysis::ResetCounter(FemModel* femmodel){/*{{{*/ - - int* elementactive_counter=NULL; - Element* element=NULL; - int N; - - femmodel->parameters->FindParam(&elementactive_counter,&N,ElementActiveCounterEnum); - for(int i=0;ielements->Size();i++){ - element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); - - elementactive_counter[i]=0; - } - femmodel->parameters->SetParam(elementactive_counter,femmodel->elements->Size(),ElementActiveCounterEnum); - xDelete(elementactive_counter); -}/*}}}*/ +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.h (revision 18718) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.h (revision 18719) @@ -41,5 +41,4 @@ IssmDouble GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* thick_input, Input* base_input); void ElementizeEplMask(FemModel* femmodel); - void ResetCounter(FemModel* femmodel); }; #endif Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 18718) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 18719) @@ -1809,21 +1809,25 @@ void FemModel::HydrologyEPLupdateDomainx(void){ /*{{{*/ - Vector* mask = NULL; - IssmDouble* serial_mask = NULL; - Vector* active = NULL; - IssmDouble* serial_active = NULL; - + Vector* mask = NULL; + IssmDouble* serial_mask = NULL; + Vector* active = NULL; + IssmDouble* serial_active = NULL; + int* eplzigzag_counter = NULL; + 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)); + mask=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,element); - } - + effanalysis->HydrologyEPLGetMask(mask,eplzigzag_counter,element); + } + printf("POINTER = %p\n",eplzigzag_counter); + this->parameters->SetParam(eplzigzag_counter,this->nodes->Size(),EplZigZagCounterEnum); /*Assemble and serialize*/ mask->Assemble(); serial_mask=mask->ToMPISerial(); + xDelete(eplzigzag_counter); delete mask; Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 18718) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 18719) @@ -109,5 +109,5 @@ EplHeadSlopeXEnum, EplHeadSlopeYEnum, - ElementActiveCounterEnum, + EplZigZagCounterEnum, HydrologydcMaxIterEnum, HydrologydcRelTolEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 18718) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 18719) @@ -117,5 +117,5 @@ case EplHeadSlopeXEnum : return "EplHeadSlopeX"; case EplHeadSlopeYEnum : return "EplHeadSlopeY"; - case ElementActiveCounterEnum : return "ElementActiveCounter"; + case EplZigZagCounterEnum : return "EplZigZagCounter"; case HydrologydcMaxIterEnum : return "HydrologydcMaxIter"; case HydrologydcRelTolEnum : return "HydrologydcRelTol"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 18718) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 18719) @@ -117,5 +117,5 @@ else if (strcmp(name,"EplHeadSlopeX")==0) return EplHeadSlopeXEnum; else if (strcmp(name,"EplHeadSlopeY")==0) return EplHeadSlopeYEnum; - else if (strcmp(name,"ElementActiveCounter")==0) return ElementActiveCounterEnum; + else if (strcmp(name,"EplZigZagCounter")==0) return EplZigZagCounterEnum; else if (strcmp(name,"HydrologydcMaxIter")==0) return HydrologydcMaxIterEnum; else if (strcmp(name,"HydrologydcRelTol")==0) return HydrologydcRelTolEnum; Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 18718) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 18719) @@ -69,4 +69,5 @@ /*Initialize the element mask*/ inefanalysis->ElementizeEplMask(femmodel); + effanalysis->InitZigZagCounter(femmodel); } /*The real computation starts here, outermost loop is on the two layer system*/ @@ -222,5 +223,5 @@ InputUpdateFromConstantx(femmodel,eplconverged,ConvergedEnum); InputUpdateFromSolutionx(femmodel,ug_epl); - inefanalysis->ResetCounter(femmodel); + effanalysis->ResetCounter(femmodel); break; }