Index: ../trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp =================================================================== --- ../trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 15222) +++ ../trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 15223) @@ -10,14 +10,14 @@ void solutionsequence_hydro_nonlinear(FemModel* femmodel){ /*solution : */ - Vector* ug=NULL; - Vector* uf_sed=NULL; - Vector* uf_epl=NULL; + Vector* ug_sed=NULL; + Vector* ug_epl=NULL; + Vector* uf=NULL; Vector* old_uf=NULL; - Vector* aged_uf_sed=NULL; - Vector* aged_uf_epl=NULL; + Vector* aged_ug_sed=NULL; + Vector* aged_ug_epl=NULL; Vector* ys=NULL; - Vector* duf=NULL; + Vector* dug=NULL; Matrix* Kff=NULL; Matrix* Kfs=NULL; @@ -49,22 +49,21 @@ /*Iteration on the two layers + transfer*/ femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum); - GetSolutionFromInputsx(&ug, femmodel->elements, femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters); - Reducevectorgtofx(&uf_sed, ug, femmodel->nodes,femmodel->parameters); _assert_(uf_sed); + GetSolutionFromInputsx(&ug_sed, femmodel->elements, femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters); if(isefficientlayer) { femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum); - GetSolutionFromInputsx(&ug, femmodel->elements, femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters); - Reducevectorgtofx(&uf_epl, ug, femmodel->nodes,femmodel->parameters); - } + GetSolutionFromInputsx(&ug_epl, femmodel->elements, femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters); +} for(;;){ sedcount=1; eplcount=1; //save pointer to old velocity - delete aged_uf_sed;aged_uf_sed=uf_sed; + delete aged_ug_sed; + aged_ug_sed=ug_sed; if(isefficientlayer){ - delete aged_uf_epl; - aged_uf_epl=uf_epl; + delete aged_ug_epl; + aged_ug_epl=ug_epl; } femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum); @@ -72,21 +71,21 @@ InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,false,ConvergedEnum); femmodel->UpdateConstraintsx(); femmodel->parameters->SetParam(HydrologySedimentEnum,HydrologyLayerEnum); - /* femmodel->HydrologyTransferx(); - */ + /*Iteration on the sediment layer*/ for(;;){ femmodel->HydrologyTransferx(); femmodel->SystemMatricesx(&Kff, &Kfs, &pf,&df, &sediment_kmax); CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCInefficientAnalysisEnum); Reduceloadx(pf,Kfs,ys); delete Kfs; - if(sedcount>1)delete uf_sed; - Solverx(&uf_sed, Kff, pf,old_uf, df, femmodel->parameters); - delete old_uf; old_uf=uf_sed->Duplicate(); - delete Kff; delete pf; delete ug; delete df; + 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*/ + delete Kff; delete pf; delete df; - Mergesolutionfromftogx(&ug,uf_sed,ys,femmodel->nodes,femmodel->parameters); delete ys; - InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug); + Mergesolutionfromftogx(&ug_sed,uf,ys,femmodel->nodes,femmodel->parameters); delete ys; + InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug_sed); ConstraintsStatex(&constraints_converged,&num_unstable_constraints, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); if (!sedconverged){ @@ -101,7 +100,7 @@ if(sedconverged){ femmodel->parameters->SetParam(sediment_kmax,HydrologySedimentKmaxEnum); InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,sedconverged,ConvergedEnum); - InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug); + InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug_sed); InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,sediment_kmax,HydrologySedimentKmaxEnum); break; } @@ -121,13 +120,14 @@ femmodel->SystemMatricesx(&Kff, &Kfs, &pf,&df,NULL); CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCEfficientAnalysisEnum); Reduceloadx(pf,Kfs,ys); delete Kfs; - if(eplcount>1) delete uf_epl; - Solverx(&uf_epl, Kff, pf,old_uf, df, femmodel->parameters); - delete old_uf; old_uf=uf_epl->Duplicate(); - delete Kff;delete pf; delete ug; delete df; - Mergesolutionfromftogx(&ug,uf_epl,ys,femmodel->nodes,femmodel->parameters); delete ys; - InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug); - + delete uf; + Solverx(&uf, Kff, pf,old_uf, df, femmodel->parameters); + delete old_uf; old_uf=uf->Duplicate(); + if(eplcount>1) delete ug_epl; + delete Kff;delete pf; + delete df; + 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); ConstraintsStatex(&constraints_converged,&num_unstable_constraints, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); if (!eplconverged){ @@ -142,7 +142,7 @@ 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); + InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug_epl); break; } } @@ -150,11 +150,10 @@ /*System convergence check*/ if(!hydroconverged){ //compute norm(du)/norm(u) - _assert_(aged_uf_sed); _assert_(uf_sed); - duf=uf_sed->Duplicate(); _assert_(duf); - aged_uf_sed->Copy(duf); - duf->AYPX(uf_sed,-1.0); - ndu_sed=duf->Norm(NORM_TWO); nu_sed=aged_uf_sed->Norm(NORM_TWO); + dug=ug_sed->Duplicate(); _assert_(dug); + aged_ug_sed->Copy(dug); + dug->AYPX(ug_sed,-1.0); + ndu_sed=dug->Norm(NORM_TWO); nu_sed=aged_ug_sed->Norm(NORM_TWO); if (xIsNan(ndu_sed) || xIsNan(nu_sed)) _error_("Sed convergence criterion is NaN!"); if (!xIsNan(eps_hyd)){ if (!isefficientlayer){ @@ -168,8 +167,12 @@ } } else{ - duf=aged_uf_epl->Duplicate(); aged_uf_epl->Copy(duf); duf->AYPX(uf_epl,-1.0); - ndu_epl=duf->Norm(NORM_TWO); nu_epl=aged_uf_epl->Norm(NORM_TWO); + dug=ug_epl->Duplicate();_assert_(dug); + aged_ug_epl->Copy(dug);_assert_(aged_ug_epl); + dug->AYPX(ug_epl,-1.0); + ndu_epl=dug->Norm(NORM_TWO); + nu_epl=aged_ug_epl->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)elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug); + InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug_sed); + InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug_epl); /*Free ressources: */ - delete ug; - delete uf_sed; - delete uf_epl; + delete ug_epl; + delete ug_sed; + delete uf; delete old_uf; - delete aged_uf_sed; - delete aged_uf_epl; - delete duf; + delete aged_ug_sed; + delete aged_ug_epl; + delete dug; } Index: ../trunk-jpl/src/c/classes/DofIndexing.cpp =================================================================== --- ../trunk-jpl/src/c/classes/DofIndexing.cpp (revision 15222) +++ ../trunk-jpl/src/c/classes/DofIndexing.cpp (revision 15223) @@ -159,20 +159,7 @@ else _error_("set of enum type " << EnumToStringx(setenum) << " not supported yet!"); } /*}}}*/ -/*FUNCTION DofIndexing::Deactivate{{{*/ -void DofIndexing::Deactivate(void){ - this->active = false; - /*Constrain to 0. at this point*/ - for(int i=0;igsize;i++){ - this->f_set[i] = false; - this->s_set[i] = true; - this->svalues[i] = 0.; - } - return; -} -/*}}}*/ - /*Some of the Object functionality: */ /*FUNCTION DofIndexing::Echo{{{*/ void DofIndexing::Echo(void){ @@ -234,3 +221,29 @@ _printf_("\n"); } /*}}}*/ +/*FUNCTION DofIndexing::Deactivate{{{*/ +void DofIndexing::Deactivate(void){ + this->active = false; + + /*Constrain to 0. at this point*/ + for(int i=0;igsize;i++){ + this->f_set[i] = false; + this->s_set[i] = true; + this->svalues[i] = 0.; + } + return; +} +/*}}}*/ +/*FUNCTION DofIndexing::Activate{{{*/ +void DofIndexing::Activate(void){ + this->active = true; + + /*Constrain to 0. at this point*/ + for(int i=0;igsize;i++){ + this->f_set[i] = true; + this->s_set[i] = false; + this->svalues[i] = 0.; + } + return; +} +/*}}}*/ Index: ../trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- ../trunk-jpl/src/c/classes/FemModel.cpp (revision 15222) +++ ../trunk-jpl/src/c/classes/FemModel.cpp (revision 15223) @@ -667,8 +667,9 @@ /*}}}*/ void FemModel::UpdateConstraintsx(void){ /*{{{*/ + Element *element = NULL; IssmDouble time; - int analysis_type; + int analysis_type; /*retrieve parameters: */ parameters->FindParam(&analysis_type,AnalysisTypeEnum); @@ -677,7 +678,13 @@ /*start module: */ if(VerboseModule()) _printf0_(" Updating constraints for time: " << time << "\n"); - /*First, update dof constraints in nodes, using constraints: */ + /*First, Nodes might be activated/deactivated by element*/ + for(int i=0;iSize();i++){ + element=dynamic_cast(elements->GetObjectByOffset(i)); + element->UpdateConstraints(); + } + + /*Second, constraints might be time dependent: */ SpcNodesx(nodes,constraints,parameters,analysis_type); /*Now, update degrees of freedoms: */ Index: ../trunk-jpl/src/c/classes/Node.h =================================================================== --- ../trunk-jpl/src/c/classes/Node.h (revision 15222) +++ ../trunk-jpl/src/c/classes/Node.h (revision 15223) @@ -86,6 +86,7 @@ void GetDofList(int* poutdoflist,int approximation_enum,int setenum); void GetLocalDofList(int* poutdoflist,int approximation_enum,int setenum); void FreezeDof(int dof); + void Activate(void); void Deactivate(void); int IsFloating(); int IsGrounded(); Index: ../trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Element.h (revision 15222) +++ ../trunk-jpl/src/c/classes/Elements/Element.h (revision 15223) @@ -134,6 +134,7 @@ #ifdef _HAVE_HYDROLOGY_ virtual void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode)=0; virtual void GetHydrologyTransfer(Vector* transfer)=0; + virtual void UpdateConstraints(void)=0; #endif }; #endif Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15222) +++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15223) @@ -2819,6 +2819,26 @@ this->parameters=NULL; } /*}}}*/ +/*FUNCTION Tria::UpdateConstraints{{{*/ +void Tria::UpdateConstraints(void){ + + int analysis_type; + parameters->FindParam(&analysis_type,AnalysisTypeEnum); + + /*Skip if water element*/ + if(IsOnWater()) return; + + /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ + switch(analysis_type){ + #ifdef _HAVE_HYDROLOGY_ + case HydrologyDCEfficientAnalysisEnum: + this->HydrologyUpdateEplDomain(); + break; + #endif + } + +} +/*}}}*/ /*FUNCTION Tria::UpdatePotentialUngrounding{{{*/ int Tria::UpdatePotentialUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf){ @@ -6406,7 +6426,7 @@ GetDofList(&doflist,NoneApproximationEnum,GsetEnum); /*Get the flag to know if the efficient layer is present*/ - this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); + this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); /*Use the dof list to index into the solution vector: */ for(int i=0;iparameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); + GetInputListOnVertices(&activeEpl[0],HydrologydcMaskEplactiveEnum); + + if(isefficientlayer){ + for(int i=0;iActivate(); + else if (activeEpl[i]==0.0) + nodes[i]->Deactivate(); + else + _error_("activation value "<< activeEpl[i] <<" not supported"); + } + } +} +/*}}}*/ #endif #ifdef _HAVE_DAKOTA_ Index: ../trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 15222) +++ ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 15223) @@ -249,6 +249,8 @@ void GetSolutionFromInputsHydrologyDCInefficient(Vector* solution); void GetSolutionFromInputsHydrologyDCEfficient(Vector* solution); void CreateHydrologyWaterVelocityInput(void); + void HydrologyUpdateEplDomain(void); + void UpdateConstraints(void); void InputUpdateFromSolutionHydrology(IssmDouble* solution); void InputUpdateFromSolutionHydrologyShreve(IssmDouble* solution); void InputUpdateFromSolutionHydrologyDC(IssmDouble* solution); Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15222) +++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15223) @@ -3188,6 +3188,13 @@ return nflipped; } /*}}}*/ +/*FUNCTION Penta::UpdateConstraints{{{*/ +void Penta::UpdateConstraints(void){ + + /*Do nothing for now*/ + +} +/*}}}*/ /*FUNCTION Penta::ViscousHeatingCreateInput {{{*/ void Penta::ViscousHeatingCreateInput(void){ @@ -9230,7 +9237,6 @@ } /*}}}*/ #endif - #ifdef _HAVE_BALANCED_ /*FUNCTION Penta::CreateKMatrixBalancethickness {{{*/ ElementMatrix* Penta::CreateKMatrixBalancethickness(void){ Index: ../trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Penta.h (revision 15222) +++ ../trunk-jpl/src/c/classes/Elements/Penta.h (revision 15223) @@ -303,7 +303,7 @@ #endif #ifdef _HAVE_HYDROLOGY_ - void CreateHydrologyWaterVelocityInput(void); + void UpdateConstraints(void); void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode); void GetHydrologyTransfer(Vector* transfer); #endif Index: ../trunk-jpl/src/c/classes/DofIndexing.h =================================================================== --- ../trunk-jpl/src/c/classes/DofIndexing.h (revision 15222) +++ ../trunk-jpl/src/c/classes/DofIndexing.h (revision 15223) @@ -48,6 +48,7 @@ /*}}}*/ /*DofIndexing management: {{{*/ DofIndexing* Spawn(int* indices, int numindices); + void Activate(void); void Deactivate(void); /*}}}*/ Index: ../trunk-jpl/src/c/classes/Node.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Node.cpp (revision 15222) +++ ../trunk-jpl/src/c/classes/Node.cpp (revision 15223) @@ -498,6 +498,13 @@ } /*}}}*/ +/*FUNCTION Node::Activate{{{*/ +void Node::Activate(void){ + + indexing.Activate(); + +} +/*}}}*/ /*FUNCTION Node::GetApproximation {{{*/ int Node::GetApproximation(){