Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 17029) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 17030) @@ -58,4 +58,6 @@ iomodel->FetchDataToInput(elements,EplHeadEnum); iomodel->FetchDataToInput(elements,HydrologydcEplInitialThicknessEnum); + + // iomodel->FetchDataToInput(elements,HydrologydcMaskEplactiveNodeEnum); elements->InputDuplicate(HydrologydcEplInitialThicknessEnum,HydrologydcEplThicknessEnum); @@ -98,19 +100,24 @@ /*Finite Element Analysis*/ -void HydrologyDCEfficientAnalysis::Core(FemModel* femmodel){/*{{{*/ +void HydrologyDCEfficientAnalysis::Core(FemModel* femmodel){/*{{{*/ _error_("not implemented"); }/*}}}*/ + ElementVector* HydrologyDCEfficientAnalysis::CreateDVector(Element* element){/*{{{*/ /*Default, return NULL*/ return NULL; }/*}}}*/ + ElementMatrix* HydrologyDCEfficientAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/ _error_("Not implemented"); }/*}}}*/ + ElementMatrix* HydrologyDCEfficientAnalysis::CreateKMatrix(Element* element){/*{{{*/ /*Intermediaries*/ int meshtype; + bool active_element; Element* basalelement; + Input* active_element_input=NULL; /*Get basal element*/ @@ -127,7 +134,11 @@ } + /* active_element_input=basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input); */ + /* active_element_input->GetInputValue(&active_element); */ + /*Check that all nodes are active, else return empty matrix*/ if(!basalelement->AllActive()) { - if(meshtype!=Mesh2DhorizontalEnum){ + //if(!active_element){ + if(meshtype!=Mesh2DhorizontalEnum){ basalelement->DeleteMaterials(); delete basalelement; @@ -199,5 +210,8 @@ /*Intermediaries*/ int meshtype; + bool active_element; Element* basalelement; + Input* active_element_input=NULL; + /*Get basal element*/ @@ -213,8 +227,11 @@ default: _error_("mesh "<GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input); */ + /* active_element_input->GetInputValue(&active_element); */ /*Check that all nodes are active, else return empty matrix*/ - if(!basalelement->AllActive()){ - if(meshtype!=Mesh2DhorizontalEnum){ + if(!basalelement->AllActive()) { + //if(!active_element){ + if(meshtype!=Mesh2DhorizontalEnum){ basalelement->DeleteMaterials(); delete basalelement; @@ -289,7 +306,9 @@ return pe; }/*}}}*/ + void HydrologyDCEfficientAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ element->GetSolutionFromInputsOneDof(solution,EplHeadEnum); }/*}}}*/ + void HydrologyDCEfficientAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ int meshtype; @@ -315,4 +334,5 @@ return rho_freshwater*g*epl_porosity*(water_compressibility+(epl_compressibility/epl_porosity)); }/*}}}*/ + void HydrologyDCEfficientAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 17029) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 17030) @@ -296,4 +296,5 @@ return pe; }/*}}}*/ + void HydrologyDCInefficientAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. @@ -323,7 +324,9 @@ xDelete(dbasis); }/*}}}*/ + void HydrologyDCInefficientAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ element->GetSolutionFromInputsOneDof(solution,SedimentHeadEnum); }/*}}}*/ + void HydrologyDCInefficientAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ @@ -395,4 +398,5 @@ return rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity)); }/*}}}*/ + void HydrologyDCInefficientAnalysis::ElementizeEplMask(FemModel* femmodel){/*{{{*/ Index: /issm/trunk-jpl/src/c/analyses/L2ProjectionEPLAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/L2ProjectionEPLAnalysis.cpp (revision 17029) +++ /issm/trunk-jpl/src/c/analyses/L2ProjectionEPLAnalysis.cpp (revision 17030) @@ -73,16 +73,21 @@ _error_("not implemented"); }/*}}}*/ + ElementVector* L2ProjectionEPLAnalysis::CreateDVector(Element* element){/*{{{*/ /*Default, return NULL*/ return NULL; }/*}}}*/ + ElementMatrix* L2ProjectionEPLAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/ _error_("Not implemented"); }/*}}}*/ + ElementMatrix* L2ProjectionEPLAnalysis::CreateKMatrix(Element* element){/*{{{*/ /*Intermediaries*/ int meshtype; + bool active_element; Element* basalelement; + Input* active_element_input=NULL; /*Get basal element*/ @@ -101,4 +106,17 @@ break; default: _error_("mesh "<GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input); */ + /* active_element_input->GetInputValue(&active_element); */ + + /* Check that all nodes are active, else return empty matrix */ + if(!basalelement->AllActive()) { + // if(!active_element){ + if(meshtype!=Mesh2DhorizontalEnum){ + basalelement->DeleteMaterials(); + delete basalelement; + } + return NULL; } @@ -143,5 +161,7 @@ /*Intermediaries*/ int meshtype; + bool active_element; Element* basalelement; + Input* active_element_input=NULL; /*Get basal element*/ @@ -158,4 +178,17 @@ } + /* active_element_input=basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input); */ + /* active_element_input->GetInputValue(&active_element); */ + + /*Check that all nodes are active, else return empty matrix*/ + if(!basalelement->AllActive()) { + /* if(!active_element){ */ + if(meshtype!=Mesh2DhorizontalEnum){ + basalelement->DeleteMaterials(); + delete basalelement; + } + return NULL; + } + /*Intermediaries */ int input_enum,index; @@ -199,7 +232,9 @@ return pe; }/*}}}*/ + void L2ProjectionEPLAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ _error_("not implemented yet"); }/*}}}*/ + void L2ProjectionEPLAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ int inputenum,meshtype; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17029) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17030) @@ -4729,5 +4729,9 @@ /*Introduce relaxation*/ - wh_trans=old_transfer[i]+0.66*(wh_trans-old_transfer[i]); + // wh_trans=old_transfer[i]+0.66*(wh_trans-old_transfer[i]); + + if(this->Id()==27){ + printf("old %e, transfer is %e ,eplhead %e, sedhead %e, factor %e \n",old_transfer[i],wh_trans,epl_head[i],sed_head[i],epl_specificstoring[i]*epl_thickness[i]*sed_trans[i]); + } /*Assign output pointer*/ @@ -4782,4 +4786,10 @@ IssmDouble residual[numdof]; + bool active_element; + Input* active_element_input=NULL; + + active_element_input=inputs->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input); + active_element_input->GetInputValue(&active_element); + GetInputListOnVertices(&old_active[0],HydrologydcMaskEplactiveNodeEnum); GetInputListOnVertices(&epl_thickness[0],HydrologydcEplThicknessEnum); @@ -4808,5 +4818,6 @@ /*Increase of the efficient system is needed if the epl head reach the maximum value (sediment max value for now)*/ this->GetHydrologyDCInefficientHmax(&h_max,i); - if(eplhead[i]>=h_max && this->AnyActive()){ + //if(eplhead[i]>=h_max && this->AnyActive()){ + if(eplhead[i]>=h_max && active_element){ for(j=0;j0.){ @@ -4838,4 +4849,6 @@ IssmDouble ice_thickness[numdof],bed[numdof]; + bool active_element; + Input* active_element_input=NULL; /*Get the flag to know if the efficient layer is present*/ @@ -4854,5 +4867,9 @@ A = material->GetAbar(); - GetInputListOnVertices(&activeEpl[0],HydrologydcMaskEplactiveNodeEnum); + // GetInputListOnVertices(&activeEpl[0],HydrologydcMaskEplactiveNodeEnum); + + active_element_input=inputs->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input); + active_element_input->GetInputValue(&active_element); + GetInputListOnVertices(&eplhead[0],EplHeadEnum); GetInputListOnVertices(&epl_slopeX[0],EplHeadSlopeXEnum); @@ -4862,5 +4879,6 @@ GetInputListOnVertices(&bed[0],BedEnum); - if(!this->AnyActive()){ + //if(!this->AnyActive()){ + if(!active_element){ /*Keeping thickness to initial value if EPL is not active*/ for(int i=0;iId()==27){ + printf("old %e, thickness is %e ,eplhead %e, slopeX %e, slopeY %e\n",old_thickness[i],thickness[i],eplhead[i],epl_slopeX[i],epl_slopeY[i]); + } } } Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 17029) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 17030) @@ -169,10 +169,25 @@ femmodel->parameters->SetParam(HydrologyEfficientEnum,HydrologyLayerEnum); - eplconverged = false; - + eplconverged=false; for(;;){ ug_epl_sub_iter=ug_epl->Duplicate();_assert_(ug_epl_sub_iter); ug_epl->Copy(ug_epl_sub_iter); - + + /*Loop on the epl layer to deal with the penalization, have been removed TO CHECK*/ + + SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel); + CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCEfficientAnalysisEnum); + Reduceloadx(pf,Kfs,ys); delete Kfs; + delete uf_epl; + Solverx(&uf_epl,Kff,pf,uf_epl_sub_iter,df,femmodel->parameters); + delete Kff; delete pf; delete df; + delete uf_epl_sub_iter; + uf_epl_sub_iter=uf_epl->Duplicate(); + uf_epl->Copy(uf_epl_sub_iter); + delete ug_epl; + Mergesolutionfromftogx(&ug_epl,uf_epl,ys,femmodel->nodes,femmodel->parameters); delete ys; + InputUpdateFromSolutionx(femmodel,ug_epl); + ConstraintsStatex(&constraints_converged,&num_unstable_constraints,femmodel); + /* {{{ *//*Start by retrieving the EPL head slopes and compute EPL Thickness*/ if(VerboseSolution()) _printf0_("computing EPL Head slope...\n"); @@ -183,40 +198,11 @@ femmodel->parameters->SetParam(EplHeadSlopeYEnum,InputToL2ProjectEnum); solutionsequence_linear(femmodel); + femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum); femmodel->HydrologyEPLThicknessx(); + //updating mask after the computation of the epl thickness (Allow to close too thin EPL) femmodel->HydrologyEPLupdateDomainx(); /* }}} */ - femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum); - - /* {{{ *//*Loop on the epl layer to deal with the penalization*/ - for(;;){ - SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel); - CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCEfficientAnalysisEnum); - Reduceloadx(pf,Kfs,ys); delete Kfs; - delete uf_epl; - Solverx(&uf_epl,Kff,pf,uf_epl_sub_iter,df,femmodel->parameters); - delete Kff; delete pf; delete df; - delete uf_epl_sub_iter; - uf_epl_sub_iter=uf_epl->Duplicate(); - uf_epl->Copy(uf_epl_sub_iter); - delete ug_epl; - Mergesolutionfromftogx(&ug_epl,uf_epl,ys,femmodel->nodes,femmodel->parameters); delete ys; - InputUpdateFromSolutionx(femmodel,ug_epl); - ConstraintsStatex(&constraints_converged,&num_unstable_constraints,femmodel); - - if (!eplconverged){ - if(VerboseConvergence()) _printf0_(" # EPL unstable constraints = " << num_unstable_constraints << "\n"); - if(num_unstable_constraints==0) eplconverged = true; - if (eplcount>=hydro_maxiter){ - _error_(" maximum number of EPL iterations (" << hydro_maxiter << ") exceeded"); - } - } - /*Add an iteration and get out of the loop if the penalisation is converged*/ - eplcount++; - if(eplconverged) break; - } - /* }}} */ /*End of the EPL penalization loop*/ - - eplconverged=false; + /*Update Elemental Mask and compute transfer*/ HydrologyDCInefficientAnalysis* analysis = new HydrologyDCInefficientAnalysis(); @@ -224,5 +210,5 @@ delete analysis; femmodel->HydrologyTransferx(); - + dug=ug_epl_sub_iter->Duplicate();_assert_(dug); ug_epl_sub_iter->Copy(dug); @@ -235,5 +221,9 @@ if(VerboseConvergence()) _printf0_(setw(50) << left << " Inner EPL Convergence criterion:" << ndu_epl/nu_epl*100 << " aiming lower than " << eps_hyd*100 << " %\n"); if((ndu_epl/nu_epl)=hydro_maxiter){ + _error_(" maximum number of EPL iterations (" << hydro_maxiter << ") exceeded"); + } + eplcount++; + delete ug_epl_sub_iter; if(eplconverged){