source: issm/oecreview/Archive/16554-17801/ISSM-17443-17444.diff@ 17802

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

@@ -123 +123 @@
 
         iomodel->FetchData(1,MeshVertexonbedEnum);
 
-        //create penalties for nodes: no node can have a temperature over the melting point
+        //create penalties for nodes: no node can have water above the max
         CreateSingleNodeToElementConnectivity(iomodel);
         for(int i=0;i<iomodel->numberofvertices;i++){
-                if (iomodel->meshtype==Mesh3DEnum){
+                if (iomodel->meshtype!=Mesh3DEnum){
                         /*keep only this partition's nodes:*/
                         if(iomodel->my_vertices[i]){
                                 loads->AddObject(new Pengrid(iomodel->loadcounter+i+1,i,iomodel,HydrologyDCInefficientAnalysisEnum));
@@ -279 +279 @@
         /*Intermediaries */
         bool       active_element,isefficientlayer;
         IssmDouble dt,scalar;
-        IssmDouble moulin_load,water_head;
+        IssmDouble water_head;
         IssmDouble water_load,transfer;
         IssmDouble Jdet;
 
@@ -306 +306 @@
         Input* thickness_input   = basalelement->GetInput(ThicknessEnum);
         Input* bed_input         = basalelement->GetInput(BedEnum);
         Input* water_input       = basalelement->GetInput(BasalforcingsMeltingRateEnum);    _assert_(water_input);
-        Input* moulin_input      = basalelement->GetInput(HydrologydcBasalMoulinInputEnum); _assert_(moulin_input);
         if(dt!= 0.){old_wh_input = basalelement->GetInput(SedimentHeadOldEnum);             _assert_(old_wh_input);}
 
         IssmDouble sediment_storing    = SedimentStoring(basalelement);
@@ -326 +325 @@
 
                 /*Loading term*/
                 water_input->GetInputValue(&water_load,gauss);
-                moulin_input->GetInputValue(&moulin_load,gauss);
-        
+
                 scalar = Jdet*gauss->weight*(water_load);
-                scalar = scalar + Jdet* moulin_load;
                 if(dt!=0.) scalar = scalar*dt;
                 for(int i=0;i<numnodes;i++){
                         pe->values[i]+=scalar*basis[i];

../trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp

@@ -75 +75 @@
                 ElementVector* pe = analysis->CreatePVector(element);
                 element->ReduceMatrices(Ke,pe);
                 if(Ke) Ke->AddToGlobal(Kff,Kfs);
-                if(pe) pe->AddToGlobal(pf);
+                if(pe){
+                        pe->AddToGlobal(pf);
+                        /* printf("-------------------USUAL \n"); */
+                        /* pf->Echo(); */
+                }
                 delete Ke;
                 delete pe;
         }
@@ -86 +90 @@
                 if(load->InAnalysis(configuration_type)){
                         load->CreateKMatrix(Kff,Kfs);
                         load->CreatePVector(pf);
+                        /* printf("-------------------LOADING \n"); */
+                        /* pf->Echo(); */
                 }
         }
 
@@ -96 +102 @@
                         if(load->InAnalysis(configuration_type)){
                                 load->PenaltyCreateKMatrix(Kff,Kfs,kmax);
                                 load->PenaltyCreatePVector(pf,kmax);
+                                /* printf("-------------------PENALTY \n"); */
+                                /* pf->Echo(); */
                         }
                 }
         }

../trunk-jpl/src/c/classes/Loads/Pengrid.cpp

@@ -180 +180 @@
 /*FUNCTION Pengrid::CreatePVector {{{*/
 void  Pengrid::CreatePVector(Vector<IssmDouble>* pf){
 
-        /*No loads applied, do nothing: */
-        return;
+        ElementVector* pe=NULL;
+        int analysis_type;
+        this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
 
+        switch(analysis_type){
+                
+        case HydrologyDCInefficientAnalysisEnum:
+                pe = CreatePVectorHydrologyDCInefficient();
+                break;
+        default:
+                /*No loads applied, do nothing: */
+                return;
+        }
+        if(pe){
+                pe->AddToGlobal(pf);
+                delete pe;
+        }
+
 }
 /*}}}*/
 /*FUNCTION Pengrid::GetNodesSidList{{{*/
@@ -792 +807 @@
         return pe;
 }
 /*}}}*/
+/*FUNCTION Pengrid::CreatePVectorHydrologyDCInefficient {{{*/
+ElementVector* Pengrid::CreatePVectorHydrologyDCInefficient(void){
+
+        IssmDouble moulin_load,dt;
+
+        /*Initialize Element matrix*/
+        ElementVector* pe=new ElementVector(&node,1,this->parameters);
+
+        this->element->GetInputValue(&moulin_load,node,HydrologydcBasalMoulinInputEnum);
+        parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+
+        if(dt!=0.0) pe->values[0]=moulin_load*dt;
+
+        /*Clean up and return*/
+        return pe;
+ }
+/*}}}*/
 /*FUNCTION Pengrid::ResetConstraint {{{*/
 void  Pengrid::ResetConstraint(void){
         active         = 0;

../trunk-jpl/src/c/classes/Loads/Pengrid.h

@@ -93 +93 @@
                 ElementVector* PenaltyCreatePVectorHydrologyDCInefficient(IssmDouble kmax);
                 void           ConstraintActivateHydrologyDCInefficient(int* punstable);
                 void  ConstraintActivate(int* punstable);
+                ElementVector* CreatePVectorHydrologyDCInefficient(void);
                 void  ResetConstraint(void);
                 /*}}}*/