Index: ../trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp (revision 16760) +++ ../trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp (revision 16761) @@ -85,5 +85,28 @@ element->GetSolutionFromInputsOneDof(solution,WatercolumnEnum); }/*}}}*/ void HydrologyShreveAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ - _error_("not implemented yet"); + + /*Intermediary*/ + int* doflist = NULL; + + /*Fetch number of nodes for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Fetch dof list and allocate solution vector*/ + element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum); + IssmDouble* values = xNew(numnodes); + + /*Use the dof list to index into the solution vector: */ + for(int i=0;i(values[i])) _error_("NaN found in solution vector"); + if (values[i]<10e-10) values[i]=10e-10; //correcting the water column to positive values + } + + /*Add input to the element: */ + element->AddInput(WatercolumnEnum,values,P1Enum); + + /*Free ressources:*/ + xDelete(values); + xDelete(doflist); }/*}}}*/ Index: ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 16760) +++ ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 16761) @@ -99,5 +99,15 @@ element->GetSolutionFromInputsOneDof(solution,EplHeadEnum); }/*}}}*/ void HydrologyDCEfficientAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ - _error_("not implemented yet"); + int meshtype; + element->FindParam(&meshtype,MeshTypeEnum); + switch(meshtype){ + case Mesh2DhorizontalEnum: + element->InputUpdateFromSolutionOneDof(solution,EplHeadEnum); + break; + case Mesh3DEnum: + element->InputUpdateFromSolutionOneDofCollapsed(solution,EplHeadEnum); + break; + default: _error_("mesh "<GetSolutionFromInputsOneDof(solution,SedimentHeadEnum); }/*}}}*/ void HydrologyDCInefficientAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ - _error_("not implemented yet"); + + int meshtype; + bool converged; + int* doflist=NULL; + Element* basalelement=NULL; + + element->FindParam(&meshtype,MeshTypeEnum); + if(meshtype!=Mesh2DhorizontalEnum){ + if(!element->IsOnBed()) return; + basalelement=element->SpawnBasalElement(); + } + else{ + basalelement = element; + } + + /*Fetch number of nodes for this finite element*/ + int numnodes = basalelement->GetNumberOfNodes(); + + /*Fetch dof list and allocate solution vector*/ + basalelement->GetDofList(&doflist,NoneApproximationEnum,GsetEnum); + IssmDouble* values = xNew(numnodes); + IssmDouble* residual = xNew(numnodes); + + /*Use the dof list to index into the solution vector: */ + for(int i=0;i(values[i])) _error_("NaN found in solution vector"); + } + + /*If converged keep the residual in mind*/ + element->GetInputValue(&converged,ConvergedEnum); + + /*Get inputs*/ + if(converged){ + IssmDouble penalty_factor,kmax,kappa,h_max; + element->FindParam(&kmax,HydrologySedimentKmaxEnum); + element->FindParam(&penalty_factor,HydrologydcPenaltyFactorEnum); + + kappa=kmax*pow(10.,penalty_factor); + + for(int i=0;iGetHydrologyDCInefficientHmax(&h_max,i); + if(values[i]>h_max) residual[i] = kappa*(values[i]-h_max); + else residual[i] = 0.; + } + } + + /*Add input to the element: */ + element->AddBasalInput(SedimentHeadEnum,values,P1Enum); + element->AddBasalInput(SedimentHeadResidualEnum,residual,P1Enum); + + /*Free ressources:*/ + xDelete(values); + xDelete(residual); + xDelete(doflist); + if(meshtype!=Mesh2DhorizontalEnum) delete basalelement; }/*}}}*/ Index: ../trunk-jpl/src/c/analyses/L2ProjectionEPLAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/L2ProjectionEPLAnalysis.cpp (revision 16760) +++ ../trunk-jpl/src/c/analyses/L2ProjectionEPLAnalysis.cpp (revision 16761) @@ -73,6 +73,21 @@ _error_("not implemented yet"); }/*}}}*/ void L2ProjectionEPLAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ - _error_("not implemented yet"); + int inputenum,meshtype; + + element->FindParam(&inputenum,InputToL2ProjectEnum); + element->FindParam(&meshtype,MeshTypeEnum); + switch(meshtype){ + case Mesh2DhorizontalEnum: + element->InputUpdateFromSolutionOneDof(solution,inputenum); + break; + case Mesh2DverticalEnum: + element->InputUpdateFromSolutionOneDof(solution,inputenum); + break; + case Mesh3DEnum: + element->InputUpdateFromSolutionOneDofCollapsed(solution,inputenum); + break; + default: _error_("mesh "<* transfer)=0; virtual void HydrologyEPLGetMask(Vector* mask)=0; virtual void HydrologyEPLGetActive(Vector* active)=0; Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16760) +++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16761) @@ -7226,6 +7226,41 @@ *ph_max=h_max; } /*}}}*/ +/*FUNCTION Tria::GetHydrologyDCInefficientHmax{{{*/ +void Tria::GetHydrologyDCInefficientHmax(IssmDouble* ph_max,int index){ + + int hmax_flag; + IssmDouble h_max; + IssmDouble rho_ice,rho_water; + IssmDouble thickness,bed; + /*Get the flag to the limitation method*/ + this->parameters->FindParam(&hmax_flag,HydrologydcSedimentlimitFlagEnum); + + /*Switch between the different cases*/ + switch(hmax_flag){ + case 0: + h_max=1.0e+10; + break; + case 1: + parameters->FindParam(&h_max,HydrologydcSedimentlimitEnum); + break; + case 2: + rho_ice=matpar->GetRhoIce(); + rho_water=matpar->GetRhoFreshwater(); + this->GetInputValue(&thickness,this->nodes[index],ThicknessEnum); + this->GetInputValue(&bed,this->nodes[index],BedEnum); + h_max=((rho_ice*thickness)/rho_water)+bed; + break; + case 3: + _error_("Using normal stress not supported yet"); + break; + default: + _error_("no case higher than 3 for SedimentlimitFlag"); + } + /*Assign output pointer*/ + *ph_max=h_max; +} +/*}}}*/ /*FUNCTION Tria::GetHydrologyTransfer{{{*/ void Tria::GetHydrologyTransfer(Vector* transfer){ Index: ../trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 16760) +++ ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 16761) @@ -327,6 +327,7 @@ void InputUpdateFromSolutionHydrologyDCInefficient(IssmDouble* solution); void InputUpdateFromSolutionHydrologyDCEfficient(IssmDouble* solution); void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode); + void GetHydrologyDCInefficientHmax(IssmDouble* ph_max,int index); void GetHydrologyTransfer(Vector* transfer); void HydrologyEPLGetActive(Vector* active_vec); void HydrologyEPLGetMask(Vector* vec_mask); Index: ../trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Penta.h (revision 16760) +++ ../trunk-jpl/src/c/classes/Elements/Penta.h (revision 16761) @@ -343,6 +343,7 @@ ElementVector* CreatePVectorHydrologyDCEfficient(void); ElementVector* CreatePVectorL2ProjectionEPL(void); void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode); + void GetHydrologyDCInefficientHmax(IssmDouble* ph_max,int index){_error_("not implemented yet");}; void GetHydrologyTransfer(Vector* transfer); void HydrologyEPLGetActive(Vector* active_vec); void HydrologyEPLGetMask(Vector* vec_mask); Index: ../trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Seg.h (revision 16760) +++ ../trunk-jpl/src/c/classes/Elements/Seg.h (revision 16761) @@ -127,6 +127,7 @@ #endif #ifdef _HAVE_HYDROLOGY_ void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode){_error_("not implemented yet");}; + void GetHydrologyDCInefficientHmax(IssmDouble* ph_max,int index){_error_("not implemented yet");}; void GetHydrologyTransfer(Vector* transfer){_error_("not implemented yet");}; void HydrologyEPLGetActive(Vector* active_vec){_error_("not implemented yet");}; void HydrologyEPLGetMask(Vector* vec_mask){_error_("not implemented yet");};