Index: /issm/trunk-jpl/src/c/classes/ElementResults/TriaP1ElementResult.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/ElementResults/TriaP1ElementResult.cpp (revision 15376) +++ /issm/trunk-jpl/src/c/classes/ElementResults/TriaP1ElementResult.cpp (revision 15377) @@ -118,4 +118,4 @@ /*FUNCTION TriaP1ElementResult::GetElementVectorFromResults{{{*/ void TriaP1ElementResult::GetElementVectorFromResults(Vector* vector,int dof){ - _error_("Result " << EnumToStringx(enum_type) << " is a TriaP1ElementResult and should not write vector of size numberofelemenrs"); + _error_("Result " << EnumToStringx(enum_type) << " is a TriaP1ElementResult and should not write vector of size numberofelements"); } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15376) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15377) @@ -455,4 +455,12 @@ case MeltingAnalysisEnum: Ke=CreateKMatrixMelting(); + break; + #endif + #ifdef _HAVE_HYDROLOGY_ + case HydrologyDCInefficientAnalysisEnum: + Ke=CreateKMatrixHydrologyDCInefficient(); + break; + case HydrologyDCEfficientAnalysisEnum: + Ke=CreateKMatrixHydrologyDCEfficient(); break; #endif @@ -585,4 +593,12 @@ break; #endif + #ifdef _HAVE_HYDROLOGY_ + case HydrologyDCInefficientAnalysisEnum: + pe=CreatePVectorHydrologyDCInefficient(); + break; + case HydrologyDCEfficientAnalysisEnum: + pe=CreatePVectorHydrologyDCEfficient(); + break; + #endif default: _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet"); @@ -1162,15 +1178,26 @@ break; case DiagnosticVertAnalysisEnum: - GetSolutionFromInputsDiagnosticVert(solution); + //GetSolutionFromInputsDiagnosticVert(solution); + GetSolutionFromInputsOneDof(solution, VzEnum); break; #endif #ifdef _HAVE_THERMAL_ case ThermalAnalysisEnum: - GetSolutionFromInputsThermal(solution); + //GetSolutionFromInputsThermal(solution); + GetSolutionFromInputsOneDof(solution, TemperatureEnum); break; case EnthalpyAnalysisEnum: - GetSolutionFromInputsEnthalpy(solution); + //GetSolutionFromInputsEnthalpy(solution); + GetSolutionFromInputsOneDof(solution, EnthalpyEnum); break; #endif + #ifdef _HAVE_HYDROLOGY_ + case HydrologyDCInefficientAnalysisEnum: + GetSolutionFromInputsOneDof(solution, SedimentHeadEnum); + break; + case HydrologyDCEfficientAnalysisEnum: + GetSolutionFromInputsOneDof(solution, EplHeadEnum); + break; + #endif default: _error_("analysis: " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet"); @@ -1880,4 +1907,12 @@ break; #endif + #ifdef _HAVE_HYDROLOGY_ + case HydrologyDCInefficientAnalysisEnum: + InputUpdateFromSolutionHydrologyDCInefficient(solution); + break; + case HydrologyDCEfficientAnalysisEnum: + InputUpdateFromSolutionOneDofCollapsed(solution,EplHeadEnum); + break; + #endif default: _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet"); @@ -2085,5 +2120,17 @@ return; - default: + case NodeSIdEnum: + for(int i=0;iSid()]; + if(xIsNan(values[i])) _error_("NaN found in solution vector"); + } + /*Add input to the element: */ + this->inputs->AddInput(new PentaP1Input(name,values)); + + /*Free ressources:*/ + xDelete(doflist); + return; + + default: _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet"); } @@ -2152,6 +2199,12 @@ name==QmuMeltingEnum || name==GiaWEnum || - name==GiadWdtEnum - + name==GiadWdtEnum || + name==SedimentHeadEnum || + name==EplHeadEnum || + name==SedimentHeadOldEnum || + name==EplHeadOldEnum || + name==HydrologydcMaskEplactiveEnum || + name==WaterTransferEnum + ) { return true; @@ -9267,23 +9320,197 @@ #endif #ifdef _HAVE_HYDROLOGY_ + +/*FUNCTION Penta::CreateKMatrixHydrologyDCInefficient {{{*/ +ElementMatrix* Penta::CreateKMatrixHydrologyDCInefficient(void){ + + if (!IsOnBed()) return NULL; + + Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. + ElementMatrix* Ke=tria->CreateKMatrixHydrologyDCInefficient(); + delete tria->material; delete tria; + + /*clean up and return*/ + return Ke; +} +/*}}}*/ +/*FUNCTION Penta::CreateKMatrixHydrologyDCEfficient {{{*/ +ElementMatrix* Penta::CreateKMatrixHydrologyDCEfficient(void){ + + if (!IsOnBed()) return NULL; + + Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. + ElementMatrix* Ke=tria->CreateKMatrixHydrologyDCEfficient(); + delete tria->material; delete tria; + + /*clean up and return*/ + return Ke; +} +/*}}}*/ +/*FUNCTION Penta::CreatePVectorHydrologyDCInefficient {{{*/ +ElementVector* Penta::CreatePVectorHydrologyDCInefficient(void){ + + if (!IsOnBed()) return NULL; + + /*Call Tria function*/ + Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. + ElementVector* pe=tria->CreatePVectorHydrologyDCInefficient(); + delete tria->material; delete tria; + + /*Clean up and return*/ + return pe; +} +/*}}}*/ +/*FUNCTION Penta::CreatePVectorHydrologyDCEfficient {{{*/ +ElementVector* Penta::CreatePVectorHydrologyDCEfficient(void){ + + if (!IsOnBed()) return NULL; + + /*Call Tria function*/ + Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. + ElementVector* pe=tria->CreatePVectorHydrologyDCEfficient(); + delete tria->material; delete tria; + + /*Clean up and return*/ + return pe; +} +/*}}}*/ /*FUNCTION Penta::GetHydrologyDCInefficientHmax{{{*/ void Penta::GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode){ - _error_("Hydrological stuff not suported in Penta"); + + if (!IsOnBed()) return; + + Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. + tria->GetHydrologyDCInefficientHmax(ph_max,innode); + delete tria->material; delete tria; } /*}}}*/ /*FUNCTION Penta::GetHydrologyTransfer{{{*/ void Penta::GetHydrologyTransfer(Vector* transfer){ - _error_("Hydrological stuff not suported in Penta"); -} -/*}}}*/ - + + if (!IsOnBed()) return; + + Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. + tria->GetHydrologyTransfer(transfer); + delete tria->material; delete tria; +} +/*}}}*/ +/*FUNCTION Penta::GetSolutionFromInputsOneDof {{{*/ +void Penta::GetSolutionFromInputsOneDof(Vector* solution, int enum_type){ + + const int numdof=NDOF1*NUMVERTICES; + + int i; + int* doflist=NULL; + IssmDouble values[numdof]; + IssmDouble enum_value; + GaussPenta *gauss=NULL; + + /*Get dof list: */ + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); + Input* enum_input=inputs->GetInput(enum_type); _assert_(enum_input); + + gauss=new GaussPenta(); + for(i=0;iGaussVertex(i); + enum_input->GetInputValue(&enum_value,gauss); + values[i]=enum_value; + } + + /*Add value to global vector*/ + solution->SetValues(numdof,doflist,values,INS_VAL); + + /*Free ressources:*/ + delete gauss; + xDelete(doflist); +} +/*{{{*/ /*FUNCTION Penta::HydrologyEPLGetActive {{{*/ void Penta::HydrologyEPLGetActive(Vector* active_vec){ - _error_("Hydrological stuff not suported in Penta"); + + if (!IsOnBed())return; + + Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. + tria->HydrologyEPLGetActive(active_vec); + delete tria->material; delete tria; + } /*}}}*/ /*FUNCTION Penta::HydrologyEPLGetMask{{{*/ void Penta::HydrologyEPLGetMask(Vector* vec_mask){ - _error_("Hydrological stuff not suported in Penta"); + + if (!IsOnBed())return; + + Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. + tria->HydrologyEPLGetMask(vec_mask); + delete tria->material; delete tria; + +} +/*}}}*/ +/*FUNCTION Penta::InputUpdateFromSolutionHydrologyDCInefficient{{{*/ +void Penta::InputUpdateFromSolutionHydrologyDCInefficient(IssmDouble* solution){ + + const int numdof = NDOF1*NUMVERTICES; + const int numdof2d = NDOF1*NUMVERTICES2D; + int* doflist = NULL; + bool converged; + IssmDouble values[numdof]; + IssmDouble residual[numdof]; + IssmDouble penalty_factor; + IssmDouble kmax, kappa, h_max; + Penta *penta = NULL; + + /*If not on bed, return*/ + if (!IsOnBed()) return; + + /*Get dof list: */ + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); + + /*Use the dof list to index into the solution vector and extrude it */ + for(int i=0;i(values[i])) _error_("NaN found in solution vector"); + } + + /*If converged keep the residual in mind*/ + this->inputs->GetInputValue(&converged,ConvergedEnum); + + /*Get inputs*/ + if(converged){ + this->parameters->FindParam(&kmax,HydrologySedimentKmaxEnum); + this->parameters->FindParam(&penalty_factor,HydrologydcPenaltyFactorEnum); + + kappa=kmax*pow(10.,penalty_factor); + + Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. + for(int i=0;iGetHydrologyDCInefficientHmax(&h_max,nodes[i]); + if(values[i]>h_max){ + residual[i]=kappa*(values[i]-h_max); + } + else{ + residual[i]=0.0; + } + } + delete tria->material; delete tria; + } + + /*Start looping over all elements above current element and update all inputs*/ + penta=this; + for(;;){ + /*Add input to the element: */ + penta->inputs->AddInput(new PentaP1Input(SedimentHeadEnum,values)); + penta->inputs->AddInput(new PentaP1Input(SedimentHeadResidualEnum,residual)); + + /*Stop if we have reached the surface*/ + if (penta->IsOnSurface()) break; + + /* get upper Penta*/ + penta=penta->GetUpperElement(); _assert_(penta->Id()!=this->id); + } + + /*Free ressources:*/ + xDelete(doflist); } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 15376) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 15377) @@ -302,9 +302,17 @@ #ifdef _HAVE_HYDROLOGY_ - void UpdateConstraints(void); + + ElementMatrix* CreateKMatrixHydrologyDCInefficient(void); + ElementMatrix* CreateKMatrixHydrologyDCEfficient(void); + ElementVector* CreatePVectorHydrologyDCInefficient(void); + ElementVector* CreatePVectorHydrologyDCEfficient(void); + void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode); void GetHydrologyTransfer(Vector* transfer); + void GetSolutionFromInputsOneDof(Vector* solution, int enum_type); void HydrologyEPLGetActive(Vector* active_vec); void HydrologyEPLGetMask(Vector* vec_mask); + void InputUpdateFromSolutionHydrologyDCInefficient(IssmDouble* solution); + void UpdateConstraints(void); #endif #ifdef _HAVE_THERMAL_ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15376) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15377) @@ -6406,5 +6406,4 @@ int *doflist = NULL; bool converged; - bool isefficientlayer; IssmDouble values[numdof]; IssmDouble residual[numdof]; @@ -6415,7 +6414,4 @@ GetDofList(&doflist,NoneApproximationEnum,GsetEnum); - /*Get the flag to know if the efficient layer is present*/ - this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); - /*Use the dof list to index into the solution vector: */ for(int i=0;iparameters->FindParam(&dt,TimesteppingTimeStepEnum); this->parameters->FindParam(&kmax,HydrologySedimentKmaxEnum); this->parameters->FindParam(&penalty_factor,HydrologydcPenaltyFactorEnum); @@ -6461,5 +6456,4 @@ IssmDouble rho_ice,rho_water; IssmDouble thickness,bed; - /*Get the flag to the limitation method*/ this->parameters->FindParam(&hmax_flag,HydrologydcSedimentlimitFlagEnum); Index: /issm/trunk-jpl/src/c/classes/Node.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Node.cpp (revision 15376) +++ /issm/trunk-jpl/src/c/classes/Node.cpp (revision 15377) @@ -122,5 +122,7 @@ analysis_type==BedSlopeAnalysisEnum || analysis_type==SurfaceSlopeAnalysisEnum || - analysis_type==BalancethicknessAnalysisEnum + analysis_type==BalancethicknessAnalysisEnum || + analysis_type==HydrologyDCInefficientAnalysisEnum || + analysis_type==HydrologyDCEfficientAnalysisEnum ){ if (dim==3){