Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15351) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15352) @@ -6132,5 +6132,4 @@ gauss->GaussPoint(ig); - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); @@ -6157,5 +6156,4 @@ } } - /*Clean up and return*/ delete gauss; @@ -6244,9 +6242,7 @@ /* Start looping on the number of gaussian points: */ gauss=new GaussTria(2); - for(int ig=gauss->begin();igend();ig++){ gauss->GaussPoint(ig); - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); GetNodalFunctions(basis, gauss); @@ -6311,5 +6307,4 @@ gauss->GaussPoint(ig); - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); GetNodalFunctions(basis, gauss); @@ -6691,5 +6686,5 @@ /*Increase of the efficient system is needed if the epl head reach the maximum value (sediment max value for now)*/ this->GetHydrologyDCInefficientHmax(&h_max,nodes[i]); - if(eplhead[i]>=h_max && this->AnyActive()){ + if(eplhead[i]>=h_max && this->AnyActive()){ for(j=0;j0.){ Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 15351) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 15352) @@ -853,4 +853,5 @@ element->CreateKMatrix(Kff_temp,NULL); } + for (i=0;iloads->Size();i++){ load=dynamic_cast(this->loads->GetObjectByOffset(i)); @@ -872,4 +873,5 @@ element->CreateKMatrix(Kff,Kfs); } + for (i=0;iloads->Size();i++){ load=dynamic_cast(this->loads->GetObjectByOffset(i)); Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 15351) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 15352) @@ -9,15 +9,15 @@ void solutionsequence_hydro_nonlinear(FemModel* femmodel){ - /*solution : */ Vector* ug_sed=NULL; Vector* ug_epl=NULL; Vector* uf=NULL; - Vector* old_uf=NULL; - Vector* aged_ug_sed=NULL; - Vector* aged_ug_epl=NULL; + Vector* uf_int_iter=NULL; + Vector* ug_sed_main_iter=NULL; + Vector* ug_epl_main_iter=NULL; Vector* ys=NULL; Vector* dug=NULL; - + Vector* old_ug=NULL; + Matrix* Kff=NULL; Matrix* Kfs=NULL; @@ -34,5 +34,5 @@ IssmDouble eps_hyd; IssmDouble ndu_sed,nu_sed; - IssmDouble ndu_epl,nu_epl; + IssmDouble ndu_epl,nu_epl; /*Recover parameters: */ @@ -41,5 +41,5 @@ femmodel->parameters->FindParam(&eps_hyd,HydrologydcRelTolEnum); femmodel->parameters->FindParam(&time,TimeEnum); - hydro_maxiter=100; + hydro_maxiter=150; hydrocount=1; @@ -57,9 +57,9 @@ eplcount=1; //save pointer to old velocity - delete aged_ug_sed; - aged_ug_sed=ug_sed; + delete ug_sed_main_iter; + ug_sed_main_iter=ug_sed; if(isefficientlayer){ - delete aged_ug_epl; - aged_ug_epl=ug_epl; + delete ug_epl_main_iter; + ug_epl_main_iter=ug_epl; } @@ -78,7 +78,7 @@ Reduceloadx(pf,Kfs,ys); delete Kfs; delete uf; - Solverx(&uf, Kff, pf,old_uf, df, femmodel->parameters); - delete old_uf; old_uf=uf->Duplicate(); - if(sedcount>1)delete ug_sed; /*Not on first time to avoid deleting aged_ug_sed*/ + Solverx(&uf, Kff, pf,uf_int_iter, df, femmodel->parameters); + delete uf_int_iter; uf_int_iter=uf->Duplicate(); + if(sedcount>1)delete ug_sed; /*Not on first time to avoid deleting ug_sed_main_iter*/ delete Kff; delete pf; delete df; @@ -120,6 +120,6 @@ Reduceloadx(pf,Kfs,ys); delete Kfs; delete uf; - Solverx(&uf, Kff, pf,old_uf, df, femmodel->parameters); - delete old_uf; old_uf=uf->Duplicate(); + Solverx(&uf, Kff, pf,uf_int_iter, df, femmodel->parameters); + delete uf_int_iter; uf_int_iter=uf->Duplicate(); if(eplcount>1) delete ug_epl; delete Kff;delete pf; @@ -127,6 +127,6 @@ Mergesolutionfromftogx(&ug_epl,uf,ys,femmodel->nodes,femmodel->parameters); delete ys; InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug_epl); - //femmodel->HydrologyEPLupdateDomainx(); ConstraintsStatex(&constraints_converged,&num_unstable_constraints, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + femmodel->HydrologyEPLupdateDomainx(); if (!eplconverged){ @@ -138,10 +138,9 @@ } eplcount++; - + if(eplconverged){ InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,eplconverged,ConvergedEnum); InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,sediment_kmax,MeltingOffsetEnum); InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug_epl); - //femmodel->HydrologyEPLupdateDomainx(); break; } @@ -152,12 +151,12 @@ //compute norm(du)/norm(u) dug=ug_sed->Duplicate(); _assert_(dug); - aged_ug_sed->Copy(dug); + ug_sed_main_iter->Copy(dug); dug->AYPX(ug_sed,-1.0); - ndu_sed=dug->Norm(NORM_TWO); nu_sed=aged_ug_sed->Norm(NORM_TWO); + ndu_sed=dug->Norm(NORM_TWO); nu_sed=ug_sed_main_iter->Norm(NORM_TWO); if (xIsNan(ndu_sed) || xIsNan(nu_sed)) _error_("Sed convergence criterion is NaN!"); if (!xIsNan(eps_hyd)){ if (!isefficientlayer){ if ((ndu_sed/nu_sed)Duplicate();_assert_(dug); - aged_ug_epl->Copy(dug);_assert_(aged_ug_epl); + ug_epl_main_iter->Copy(dug);_assert_(ug_epl_main_iter); dug->AYPX(ug_epl,-1.0); ndu_epl=dug->Norm(NORM_TWO); - nu_epl=aged_ug_epl->Norm(NORM_TWO); + nu_epl=ug_epl_main_iter->Norm(NORM_TWO); if (xIsNan(ndu_epl) || xIsNan(nu_epl)) _error_("EPL convergence criterion is NaN!"); if (ndu_epl==0.0 && nu_epl==0.0) nu_epl=1.0e-6; /*Hacking the case where the EPL is used but empty*/ - if ((ndu_epl/nu_epl)