Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 23660) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 23661) @@ -410,5 +410,4 @@ void HydrologyDCEfficientAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ /*Intermediaries*/ - bool active_element; int domaintype; Element* basalelement=NULL; @@ -433,29 +432,20 @@ /*Fetch dof list and allocate solution vector*/ - IssmDouble* eplHeads = xNew(numnodes); - IssmDouble* basevalue = xNew(numnodes); basalelement->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum); - - Input* active_element_input=basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input); - active_element_input->GetInputValue(&active_element); - + IssmDouble* eplHeads = xNew(numnodes); + IssmDouble* base = xNew(numnodes); + + basalelement->GetInputListOnVertices(&base[0],BaseEnum); /*Use the dof list to index into the solution vector: */ /*If the EPL is not active we revert to the bedrock elevation*/ - - if(active_element){ - for(int i=0;inodes[i]->IsActive()){ eplHeads[i]=solution[doflist[i]]; - if(xIsNan(eplHeads[i])) _error_("NaN found in solution vector"); - if(xIsInf(eplHeads[i])) _error_("Inf found in solution vector"); - } - } - else{ - //Fixing epl head at bedrock elevation when deactivating - basalelement->GetInputListOnVertices(&basevalue[0],BaseEnum); - for(int i=0;i(eplHeads[i])) _error_("NaN found in solution vector"); - if(xIsInf(eplHeads[i])) _error_("Inf found in solution vector"); - } + } + else{ + eplHeads[i]=base[i]; + } + if(xIsNan(eplHeads[i])) _error_("NaN found in solution vector"); + if(xIsInf(eplHeads[i])) _error_("Inf found in solution vector"); } /*Add input to the element: */ @@ -463,5 +453,5 @@ /*Free ressources:*/ xDelete(eplHeads); - xDelete(basevalue); + xDelete(base); xDelete(doflist); if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 23660) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 23661) @@ -488,9 +488,7 @@ /*Intermediaries*/ - bool thawed_element; int domaintype; Element* basalelement=NULL; bool converged; - int* doflist = NULL; /*Get basal element*/ @@ -506,4 +504,6 @@ default: _error_("mesh "<(numnodes); IssmDouble* residual = xNew(numnodes); - - /*Use the dof list to index into the solution vector: */ - + IssmDouble* base = xNew(numnodes); + + basalelement->GetInputListOnVertices(&base[0],BaseEnum); + /*Use the dof list to index into the solution vector, frozen nodes are set to initial value: */ for(int i=0;inodes[i]->IsActive()){ + values[i] =solution[doflist[i]]; + } + else{ + values[i] = base[i]; + } if(xIsNan(values[i])) _error_("NaN found in solution vector"); if(xIsInf(values[i])) _error_("Inf found in solution vector"); @@ -531,5 +537,4 @@ IssmDouble penalty_factor,kmax,kappa,h_max; IssmDouble* thickness = xNew(numnodes); - IssmDouble* base = xNew(numnodes); IssmDouble* transmitivity = xNew(numnodes); @@ -540,18 +545,10 @@ IssmDouble g = basalelement->FindParam(ConstantsGEnum); - basalelement->GetInputListOnVertices(thickness,ThicknessEnum); - basalelement->GetInputListOnVertices(base,BaseEnum); - basalelement->GetInputListOnVertices(transmitivity,HydrologydcSedimentTransmitivityEnum); + basalelement->GetInputListOnVertices(&thickness[0],ThicknessEnum); + basalelement->GetInputListOnVertices(&transmitivity[0],HydrologydcSedimentTransmitivityEnum); + kappa=kmax*pow(10.,penalty_factor); - - Input* thawed_element_input = basalelement->GetInput(HydrologydcMaskThawedEltEnum); _assert_(thawed_element_input); - thawed_element_input->GetInputValue(&thawed_element); - for(int i=0;iGetNode(i)); if(values[i]>h_max) { Index: /issm/trunk-jpl/src/c/cores/hydrology_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/hydrology_core.cpp (revision 23660) +++ /issm/trunk-jpl/src/c/cores/hydrology_core.cpp (revision 23661) @@ -64,5 +64,5 @@ femmodel->HydrologyIDSupdateDomainx(&ThawedNodes); if(ThawedNodes>0){ - init_time = time-dt; //getting the time back to the start of the timestep + init_time=time-dt; //getting the time back to the start of the timestep hydrotime=init_time; hydrodt=dt/hydroslices; //computing hydro dt from dt and a divider @@ -120,5 +120,5 @@ } } - else if (hydrology_model==HydrologyshaktiEnum){ + else if (hydrology_model==HydrologyshaktiEnum){ femmodel->SetCurrentConfiguration(HydrologyShaktiAnalysisEnum); InputDuplicatex(femmodel,HydrologyHeadEnum,HydrologyHeadOldEnum); Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 23660) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 23661) @@ -57,5 +57,5 @@ if(!isefficientlayer) hydroconverged=true; - /*Retrieve inputs as the initial state for the non linear iteration*/ + /*{{{*//*Retrieve inputs as the initial state for the non linear iteration*/ GetSolutionFromInputsx(&ug_sed,femmodel); Reducevectorgtofx(&uf_sed, ug_sed, femmodel->nodes,femmodel->parameters); @@ -65,8 +65,12 @@ femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum); GetSolutionFromInputsx(&ug_epl,femmodel); - /*Initialize the element mask*/ + /*Initialize the EPL element mask*/ inefanalysis->ElementizeEplMask(femmodel); effanalysis->InitZigZagCounter(femmodel); + /*Initialize the IDS element mask*/ + femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum); + inefanalysis->ElementizeIdsMask(femmodel); } + /*}}}*/ /*The real computation starts here, outermost loop is on the two layer system*/ for(;;){ @@ -98,5 +102,5 @@ /*{{{*//*Loop on the sediment layer to deal with the penalization*/ for(;;){ - /*{{{*/ /*Core of the computation*/ + /*{{{*//*Core of the computation*/ if(VerboseSolution()) _printf0_("Building Sediment Matrix...\n"); SystemMatricesx(&Kff,&Kfs,&pf,&df,&sediment_kmax,femmodel); @@ -130,5 +134,5 @@ sedconverged=false; - /*Checking convegence on the value of the sediment head*/ + /*Checking convergence on the value of the sediment head*/ duf=uf_sed_sub_iter->Duplicate();_assert_(duf); uf_sed_sub_iter->Copy(duf); @@ -160,5 +164,4 @@ if(VerboseSolution()) _printf0_("==updating mask...\n"); femmodel->HydrologyEPLupdateDomainx(&ThickCount); - inefanalysis->ElementizeEplMask(femmodel); ResetZigzagCounterx(femmodel); InputUpdateFromConstantx(femmodel,false,ConvergedEnum); @@ -172,5 +175,4 @@ femmodel->SetCurrentConfiguration(L2ProjectionEPLAnalysisEnum); femmodel->UpdateConstraintsL2ProjectionEPLx(&L2Count); - inefanalysis->ElementizeEplMask(femmodel); femmodel->parameters->SetParam(EplHeadSlopeXEnum,InputToL2ProjectEnum); solutionsequence_linear(femmodel); @@ -182,5 +184,4 @@ //updating mask after the computation of the epl thickness (Allow to close too thin EPL) femmodel->HydrologyEPLupdateDomainx(&ThickCount); - inefanalysis->ElementizeEplMask(femmodel); /*}}}*/ @@ -229,7 +230,6 @@ } } - /*}}}*/ /*End of the global EPL loop*/ - - /*{{{*/ /*Now dealing with the convergence of the whole system*/ + /*}}}*//*End of the global EPL loop*/ + /*{{{*//*Now dealing with the convergence of the whole system*/ if(!hydroconverged){ //compute norm(du)/norm(u)