source: issm/oecreview/Archive/16554-17801/ISSM-16760-16761.diff@ 17802

../trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp

@@ -85 +85 @@
         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<IssmDouble>(numnodes);
+
+        /*Use the dof list to index into the solution vector: */
+        for(int i=0;i<numnodes;i++){
+                values[i]=solution[doflist[i]];
+                if(xIsNan<IssmDouble>(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<IssmDouble>(values);
+        xDelete<int>(doflist);
 }/*}}}*/

../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp

@@ -99 +99 @@
         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 "<<EnumToStringx(meshtype)<<" not supported yet");
+        }
 }/*}}}*/

../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp

@@ -140 +140 @@
         element->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<IssmDouble>(numnodes);
+        IssmDouble* residual = xNew<IssmDouble>(numnodes);
+
+        /*Use the dof list to index into the solution vector: */
+        for(int i=0;i<numnodes;i++){
+                values[i] =solution[doflist[i]];
+                if(xIsNan<IssmDouble>(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;i<numnodes;i++){
+                        basalelement->GetHydrologyDCInefficientHmax(&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<IssmDouble>(values);
+        xDelete<IssmDouble>(residual);
+        xDelete<int>(doflist);
+        if(meshtype!=Mesh2DhorizontalEnum) delete basalelement;
 }/*}}}*/

../trunk-jpl/src/c/analyses/L2ProjectionEPLAnalysis.cpp

@@ -73 +73 @@
            _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 "<<EnumToStringx(meshtype)<<" not supported yet");
+        }
 }/*}}}*/
 

../trunk-jpl/src/c/classes/Elements/Element.h

@@ -168 +168 @@
 
                 #ifdef _HAVE_HYDROLOGY_
                 virtual void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode)=0;
+                virtual void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, int index)=0;
                 virtual void GetHydrologyTransfer(Vector<IssmDouble>* transfer)=0; 
                 virtual void HydrologyEPLGetMask(Vector<IssmDouble>* mask)=0;
                 virtual void HydrologyEPLGetActive(Vector<IssmDouble>* active)=0;

../trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -7226 +7226 @@
                 *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<IssmDouble>* transfer){
 

../trunk-jpl/src/c/classes/Elements/Tria.h

@@ -327 +327 @@
                 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<IssmDouble>* transfer);
                 void           HydrologyEPLGetActive(Vector<IssmDouble>* active_vec);
                 void           HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask);

../trunk-jpl/src/c/classes/Elements/Penta.h

@@ -343 +343 @@
                 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<IssmDouble>* transfer);
                 void           HydrologyEPLGetActive(Vector<IssmDouble>* active_vec);
                 void           HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask);

../trunk-jpl/src/c/classes/Elements/Seg.h

@@ -127 +127 @@
                 #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<IssmDouble>* transfer){_error_("not implemented yet");};
                 void    HydrologyEPLGetActive(Vector<IssmDouble>* active_vec){_error_("not implemented yet");};
                 void    HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask){_error_("not implemented yet");};