Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 19093) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 19094) @@ -672,6 +672,4 @@ int eplflip_lock; int numnodes =basalelement->GetNumberOfNodes(); - IssmDouble new_active; - IssmDouble recurent; IssmDouble* epl_thickness =xNew(numnodes); IssmDouble* old_active =xNew(numnodes); @@ -700,23 +698,21 @@ for(i=0;iGetValue(&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); + /*If node is now closed bring its thickness back to initial*/ + if (old_active[i]==0.){ + epl_thickness[i]=init_thick; + } /*Now starting to look at the activations*/ if(residual[i]>0.){ vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL); - new_active=1.0; + if(old_active[i]==0.) recurence->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL); } /*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]SetValue(basalelement->nodes[i]->Sid(),0.,INS_VAL); - new_active=0.0; + recurence->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL); } } @@ -728,19 +724,7 @@ if(sedhead[j] == sedheadmin){ vec_mask->SetValue(basalelement->nodes[j]->Sid(),1.,INS_VAL); + if(old_active[i]==0.) recurence->SetValue(basalelement->nodes[i]->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/LevelsetAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 19093) +++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 19094) @@ -556,5 +556,5 @@ Vector* vec_dist_zerolevelset = NULL; - GetVectorFromInputsx(&vec_dist_zerolevelset, femmodel, MaskIceLevelsetEnum, VertexEnum); + GetVectorFromInputsx(&vec_dist_zerolevelset, femmodel, MaskIceLevelsetEnum, VertexPIdEnum); /* set NaN on elements intersected by zero levelset */ Index: /issm/trunk-jpl/src/c/analyses/LsfReinitializationAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/LsfReinitializationAnalysis.cpp (revision 19093) +++ /issm/trunk-jpl/src/c/analyses/LsfReinitializationAnalysis.cpp (revision 19094) @@ -417,5 +417,5 @@ Vector* vec_dist_zerolevelset = NULL; - GetVectorFromInputsx(&vec_dist_zerolevelset, femmodel, MaskIceLevelsetEnum, VertexEnum); + GetVectorFromInputsx(&vec_dist_zerolevelset, femmodel, MaskIceLevelsetEnum, VertexPIdEnum); /* set distance on elements intersected by zero levelset */ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 19093) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 19094) @@ -862,18 +862,69 @@ } /*}}}*/ -void Element::GetVectorFromInputs(Vector* vector,int input_enum){/*{{{*/ +/* void Element::GetVectorFromInputs(Vector* vector,int input_enum){/\*{{{*\/ */ + +/* /\*Fetch number vertices for this element and allocate arrays*\/ */ +/* int numvertices = this->GetNumberOfVertices(); */ +/* int* vertexpidlist = xNew(numvertices); */ +/* IssmDouble* values = xNew(numvertices); */ + +/* /\*Fill in values*\/ */ +/* this->GetVertexPidList(vertexpidlist); */ +/* this->GetInputListOnVertices(values,input_enum); */ +/* vector->SetValues(numvertices,vertexpidlist,values,INS_VAL); */ + +/* /\*Clean up*\/ */ +/* xDelete(vertexpidlist); */ +/* xDelete(values); */ + +/* } */ +/* /\*}}}*\/ */ +void Element::GetVectorFromInputs(Vector* vector,int input_enum,int type){/*{{{*/ /*Fetch number vertices for this element and allocate arrays*/ - int numvertices = this->GetNumberOfVertices(); - int* vertexpidlist = xNew(numvertices); - IssmDouble* values = xNew(numvertices); - - /*Fill in values*/ - this->GetVertexPidList(vertexpidlist); - this->GetInputListOnVertices(values,input_enum); - vector->SetValues(numvertices,vertexpidlist,values,INS_VAL); - + int numvertices = this->GetNumberOfVertices(); + int numnodes = this->GetNumberOfNodes(); + int* doflist = NULL; + IssmDouble* values = NULL; + + switch(type){ + case VertexPIdEnum: + doflist = xNew(numvertices); + values = xNew(numvertices); + /*Fill in values*/ + this->GetVertexPidList(doflist); + this->GetInputListOnVertices(values,input_enum); + vector->SetValues(numvertices,doflist,values,INS_VAL); + break; + case VertexSIdEnum: + doflist = xNew(numvertices); + values = xNew(numvertices); + /*Fill in values*/ + this->GetVerticesSidList(doflist); + this->GetInputListOnVertices(values,input_enum); + vector->SetValues(numvertices,doflist,values,INS_VAL); + break; + case NodesEnum: + doflist = xNew(numnodes); + values = xNew(numnodes); + /*Fill in values*/ + this->GetInputListOnNodes(values,input_enum); + this->GetDofList(&doflist,NoneApproximationEnum,GsetEnum); + vector->SetValues(numnodes,doflist,values,INS_VAL); + break; + case NodeSIdEnum: + doflist = xNew(numnodes); + values = xNew(numnodes); + /*Fill in values*/ + this->GetNodesSidList(doflist); + this->GetInputListOnNodes(values,input_enum); + vector->SetValues(numnodes,doflist,values,INS_VAL); + break; + default: + _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet"); + } + /*Clean up*/ - xDelete(vertexpidlist); + xDelete(doflist); xDelete(values); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 19093) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 19094) @@ -96,5 +96,5 @@ void GetNodesSidList(int* sidlist); void GetPhi(IssmDouble* phi, IssmDouble* epsilon, IssmDouble viscosity); - void GetVectorFromInputs(Vector* vector, int name_enum); + void GetVectorFromInputs(Vector* vector, int name_enum, int type); void GetVertexPidList(int* pidlist); void GetVerticesConnectivityList(int* connectivitylist); Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 19093) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 19094) @@ -1943,5 +1943,4 @@ case NodesEnum: - /*Get number of nodes and dof list: */ numnodes = this->NumberofNodes(this->element_type); @@ -1957,5 +1956,4 @@ case NodeSIdEnum: - /*Get number of nodes and dof list: */ numnodes = this->NumberofNodes(this->element_type); Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 19093) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 19094) @@ -611,6 +611,6 @@ /*get vertex vectors for bed and thickness: */ - GetVectorFromInputsx(&surface ,this, SurfaceEnum,VertexEnum); - GetVectorFromInputsx(&bed ,this, BaseEnum, VertexEnum); + GetVectorFromInputsx(&surface ,this, SurfaceEnum,VertexPIdEnum); + GetVectorFromInputsx(&bed ,this, BaseEnum, VertexPIdEnum); /*Allocate vector*/ @@ -1831,5 +1831,5 @@ /*this response was scaled. pick up the response from the inputs: */ - GetVectorFromInputsx(&vertex_response,this, StringToEnumx(root),VertexEnum); + GetVectorFromInputsx(&vertex_response,this, StringToEnumx(root),VertexPIdEnum); /*Now, average it onto the partition nodes: */ @@ -1914,23 +1914,43 @@ Vector* mask = NULL; Vector* recurence = NULL; + Vector* active = NULL; IssmDouble* serial_mask = NULL; - Vector* active = NULL; + IssmDouble* serial_rec = NULL; IssmDouble* serial_active = NULL; + IssmDouble* old_active = NULL; int* eplzigzag_counter = NULL; - + int eplflip_lock; + HydrologyDCEfficientAnalysis* effanalysis = new HydrologyDCEfficientAnalysis(); 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); + this->parameters->FindParam(&eplflip_lock,HydrologydcEplflipLockEnum); + GetVectorFromInputsx(&old_active,this,HydrologydcMaskEplactiveNodeEnum,NodeSIdEnum); + for (int i=0;iSize();i++){ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); effanalysis->HydrologyEPLGetMask(mask,recurence,eplzigzag_counter,element); } + /*check for changes and increment zigzag counter, change the mask if necessary*/ + recurence->Assemble(); + serial_rec=recurence->ToMPISerial(); + for (int i=0;iSize();i++){ + Node* node=xDynamicCast(nodes->GetObjectByOffset(i)); + if(serial_rec[node->Sid()]==1.)eplzigzag_counter[node->Lid()] ++; + if(eplzigzag_counter[node->Lid()]>eplflip_lock & eplflip_lock!=0){ + mask->SetValue(node->Sid(),old_active[node->Sid()],INS_VAL); + } + } + + this->parameters->SetParam(eplzigzag_counter,this->nodes->Size(),EplZigZagCounterEnum); /*Assemble and serialize*/ mask->Assemble(); - serial_mask=mask->ToMPISerial(); + serial_mask=mask->ToMPISerial(); + xDelete(eplzigzag_counter); delete mask; Index: /issm/trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp (revision 19093) +++ /issm/trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp (revision 19094) @@ -7,22 +7,23 @@ #include "../../toolkits/toolkits.h" -void GetVectorFromInputsx( Vector** pvector,FemModel* femmodel,int name,int type){ +void GetVectorFromInputsx(Vector** pvector,FemModel* femmodel,int name,int type){ int i; Vector* vector=NULL; - if(type==VertexEnum){ - - /*Allocate vector*/ + switch(type){ + case VertexPIdEnum: case VertexSIdEnum: vector=new Vector(femmodel->vertices->NumberOfVertices()); - - /*Look up in elements*/ - for(i=0;ielements->Size();i++){ - Element* element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); - element->GetVectorFromInputs(vector,name); - } + break; + case NodesEnum:case NodeSIdEnum: + vector=new Vector(femmodel->nodes->NumberOfNodes()); + break; + default: + _error_("vector type: " << EnumToStringx(type) << " not supported yet!"); } - else{ - _error_("vector type: " << EnumToStringx(type) << " not supported yet!"); + /*Look up in elements*/ + for(i=0;ielements->Size();i++){ + Element* element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); + element->GetVectorFromInputs(vector,name,type); } Index: /issm/trunk-jpl/src/c/modules/GroundinglineMigrationx/GroundinglineMigrationx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/GroundinglineMigrationx/GroundinglineMigrationx.cpp (revision 19093) +++ /issm/trunk-jpl/src/c/modules/GroundinglineMigrationx/GroundinglineMigrationx.cpp (revision 19094) @@ -142,5 +142,5 @@ for(i=0;iSize();i++){ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - element->GetVectorFromInputs(vec_phi,MaskGroundediceLevelsetEnum); + element->GetVectorFromInputs(vec_phi,MaskGroundediceLevelsetEnum,VertexPIdEnum); } vec_phi->Assemble(); Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_la.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_la.cpp (revision 19093) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_la.cpp (revision 19094) @@ -44,6 +44,6 @@ Vector* vel_old = NULL; Vector* vel_old_local = NULL; - GetVectorFromInputsx(&vel,femmodel,VxEnum,VertexEnum); - GetVectorFromInputsx(&pug,femmodel,PressureEnum,VertexEnum); + GetVectorFromInputsx(&vel,femmodel,VxEnum,VertexPIdEnum); + GetVectorFromInputsx(&pug,femmodel,PressureEnum,VertexPIdEnum); while(true){ @@ -70,5 +70,5 @@ /*Update solution*/ InputUpdateFromSolutionx(femmodel,ug); delete ug; - GetVectorFromInputsx(&vel,femmodel,VxEnum,VertexEnum); + GetVectorFromInputsx(&vel,femmodel,VxEnum,VertexPIdEnum); /*Check for convergence*/ Vector* dvel_local=vel_old_local->Duplicate(); vel_old_local->Copy(dvel_local); dvel_local->AYPX(vel,-1.0); @@ -101,5 +101,5 @@ /*Update solution*/ InputUpdateFromSolutionx(femmodel,pug); delete pug; - GetVectorFromInputsx(&pug,femmodel,PressureEnum,VertexEnum); + GetVectorFromInputsx(&pug,femmodel,PressureEnum,VertexPIdEnum); /*Check for convergence*/ Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_la_theta.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_la_theta.cpp (revision 19093) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_la_theta.cpp (revision 19094) @@ -41,5 +41,5 @@ Vector* vx = NULL; Vector* vx_old = NULL; - GetVectorFromInputsx(&vx,femmodel,VxEnum,VertexEnum); + GetVectorFromInputsx(&vx,femmodel,VxEnum,VertexPIdEnum); while(true){ @@ -69,5 +69,5 @@ analysis->InputUpdateFromSolutionFSXTH_d( femmodel->elements,femmodel->parameters); analysis->InputUpdateFromSolutionFSXTH_tau(femmodel->elements,femmodel->parameters); - GetVectorFromInputsx(&vx,femmodel,VxEnum,VertexEnum); + GetVectorFromInputsx(&vx,femmodel,VxEnum,VertexPIdEnum); /*Check for convergence*/