Index: /issm/trunk-jpl/src/c/classes/DofIndexing.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/DofIndexing.cpp (revision 15222) +++ /issm/trunk-jpl/src/c/classes/DofIndexing.cpp (revision 15223) @@ -160,17 +160,4 @@ } /*}}}*/ -/*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: */ @@ -235,2 +222,28 @@ } /*}}}*/ +/*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: /issm/trunk-jpl/src/c/classes/DofIndexing.h =================================================================== --- /issm/trunk-jpl/src/c/classes/DofIndexing.h (revision 15222) +++ /issm/trunk-jpl/src/c/classes/DofIndexing.h (revision 15223) @@ -49,4 +49,5 @@ /*DofIndexing management: {{{*/ DofIndexing* Spawn(int* indices, int numindices); + void Activate(void); void Deactivate(void); /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 15222) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 15223) @@ -135,4 +135,5 @@ virtual void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode)=0; virtual void GetHydrologyTransfer(Vector* transfer)=0; + virtual void UpdateConstraints(void)=0; #endif }; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15222) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15223) @@ -3189,4 +3189,11 @@ } /*}}}*/ +/*FUNCTION Penta::UpdateConstraints{{{*/ +void Penta::UpdateConstraints(void){ + + /*Do nothing for now*/ + +} +/*}}}*/ /*FUNCTION Penta::ViscousHeatingCreateInput {{{*/ void Penta::ViscousHeatingCreateInput(void){ @@ -9231,5 +9238,4 @@ /*}}}*/ #endif - #ifdef _HAVE_BALANCED_ /*FUNCTION Penta::CreateKMatrixBalancethickness {{{*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 15222) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 15223) @@ -304,5 +304,5 @@ #ifdef _HAVE_HYDROLOGY_ - void CreateHydrologyWaterVelocityInput(void); + void UpdateConstraints(void); void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode); void GetHydrologyTransfer(Vector* transfer); Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15222) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15223) @@ -2818,4 +2818,24 @@ //this->parameters: we still can't point to it, it may not even exist. Configure will handle this. 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 + } + } /*}}}*/ @@ -6407,5 +6427,5 @@ /*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: */ @@ -6591,4 +6611,28 @@ /*}}}*/ +/*FUNCTION Tria::HydrologyUpdateEplDomain{{{*/ +void Tria::HydrologyUpdateEplDomain(void){ + + /*Intermediaries*/ + const int numdof = NDOF1 *NUMVERTICES; + bool isefficientlayer; + IssmDouble activeEpl[numdof]; + + /*Get parameter and Inout list*/ + this->parameters->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 Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 15222) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 15223) @@ -250,4 +250,6 @@ void GetSolutionFromInputsHydrologyDCEfficient(Vector* solution); void CreateHydrologyWaterVelocityInput(void); + void HydrologyUpdateEplDomain(void); + void UpdateConstraints(void); void InputUpdateFromSolutionHydrology(IssmDouble* solution); void InputUpdateFromSolutionHydrologyShreve(IssmDouble* solution); Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 15222) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 15223) @@ -668,6 +668,7 @@ void FemModel::UpdateConstraintsx(void){ /*{{{*/ + Element *element = NULL; IssmDouble time; - int analysis_type; + int analysis_type; /*retrieve parameters: */ @@ -678,5 +679,11 @@ 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); Index: /issm/trunk-jpl/src/c/classes/Node.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Node.cpp (revision 15222) +++ /issm/trunk-jpl/src/c/classes/Node.cpp (revision 15223) @@ -496,4 +496,11 @@ indexing.Deactivate(); + +} +/*}}}*/ +/*FUNCTION Node::Activate{{{*/ +void Node::Activate(void){ + + indexing.Activate(); } Index: /issm/trunk-jpl/src/c/classes/Node.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Node.h (revision 15222) +++ /issm/trunk-jpl/src/c/classes/Node.h (revision 15223) @@ -87,4 +87,5 @@ void GetLocalDofList(int* poutdoflist,int approximation_enum,int setenum); void FreezeDof(int dof); + void Activate(void); void Deactivate(void); int IsFloating(); Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 15222) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 15223) @@ -11,12 +11,12 @@ /*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; @@ -50,11 +50,9 @@ /*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(;;){ @@ -62,8 +60,9 @@ 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; } @@ -73,6 +72,5 @@ femmodel->UpdateConstraintsx(); femmodel->parameters->SetParam(HydrologySedimentEnum,HydrologyLayerEnum); - /* femmodel->HydrologyTransferx(); - */ + /*Iteration on the sediment layer*/ for(;;){ @@ -81,11 +79,12 @@ 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; - - 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); + 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_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); @@ -102,5 +101,5 @@ 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; @@ -122,11 +121,12 @@ 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); @@ -143,5 +143,5 @@ 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; } @@ -151,9 +151,8 @@ 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)){ @@ -169,6 +168,10 @@ } 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*/ @@ -193,14 +196,15 @@ } - 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); + 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; }