Changeset 21213

Timestamp:

09/20/16 09:30:18 (9 years ago)

Author:

bdef

Message:

NEW:Moved the treatment of DC moulins inut to Moulin.cpp and define input in EPL when it exists

Location:

issm/trunk-jpl/src/c

Files:

• analyses/HydrologyDCEfficientAnalysis.cpp (modified) (1 diff)
• analyses/HydrologyDCInefficientAnalysis.cpp (modified) (2 diffs)
• classes/Loads/Moulin.cpp (modified) (3 diffs)
• classes/Loads/Moulin.h (modified) (1 diff)
• classes/Loads/Pengrid.cpp (modified) (3 diffs)
• classes/Loads/Pengrid.h (modified) (1 diff)

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

@@ -105 +105 @@
 void HydrologyDCEfficientAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
         /*Nothing for now*/
+
+        /*Fetch parameters: */
+        int hydrology_model;
+        iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
+        if(hydrology_model!=HydrologydcEnum) return;
+
+        if(iomodel->domaintype==Domain3DEnum) iomodel->FetchData(1,"md.mesh.vertexonbase");
+
+        //create penalties for nodes: no node can have water above the max
+        CreateSingleNodeToElementConnectivity(iomodel);
+        for(int i=0;i<iomodel->numberofvertices;i++){
+                if (iomodel->domaintype!=Domain3DEnum){
+                        /*keep only this partition's nodes:*/
+                        if(iomodel->my_vertices[i]){
+                                loads->AddObject(new Moulin(iomodel->loadcounter+i+1,i,iomodel,HydrologyDCEfficientAnalysisEnum));
+                        }
+                }
+                else if(reCast<int>(iomodel->Data("md.mesh.vertexonbase")[i])){
+                        if(iomodel->my_vertices[i]){
+                                loads->AddObject(new Moulin(iomodel->loadcounter+i+1,i,iomodel,HydrologyDCEfficientAnalysisEnum));
+                        }       
+                }
+        }
+        iomodel->DeleteData(1,"md.mesh.vertexonbase");
 }/*}}}*/
 

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

@@ -140 +140 @@
                         if(iomodel->my_vertices[i]){
                                 loads->AddObject(new Pengrid(iomodel->loadcounter+i+1,i,iomodel,HydrologyDCInefficientAnalysisEnum));
+                                loads->AddObject(new Moulin(iomodel->loadcounter+i+1,i,iomodel,HydrologyDCInefficientAnalysisEnum));
                         }
                 }
@@ -145 +146 @@
                         if(iomodel->my_vertices[i]){
                                 loads->AddObject(new Pengrid(iomodel->loadcounter+i+1,i,iomodel,HydrologyDCInefficientAnalysisEnum));
+                                loads->AddObject(new Moulin(iomodel->loadcounter+i+1,i,iomodel,HydrologyDCInefficientAnalysisEnum));
                         }       
                 }

issm/trunk-jpl/src/c/classes/Loads/Moulin.cpp

@@ -180 +180 @@
                 
         case HydrologySommersAnalysisEnum:
-                pe = this->CreatePVectorMoulin();
+                pe = this->CreatePVectorHydrologySommers();
+                break;
+        case HydrologyDCInefficientAnalysisEnum:
+                pe = CreatePVectorHydrologyDCInefficient();
+                break;
+        case HydrologyDCEfficientAnalysisEnum:
+                pe = CreatePVectorHydrologyDCEfficient();
                 break;
         default:
@@ -331 +337 @@
 /*}}}*/
 
-ElementVector* Moulin::CreatePVectorMoulin(void){/*{{{*/
+ElementVector* Moulin::CreatePVectorHydrologySommers(void){/*{{{*/
 
         /*If this node is not the master node (belongs to another partition of the
@@ -349 +355 @@
 }
 /*}}}*/
+ElementVector* Moulin::CreatePVectorHydrologyDCInefficient(void){/*{{{*/
+
+        /*If this node is not the master node (belongs to another partition of the
+         * mesh), don't add the moulin input a second time*/
+        if(node->IsClone()) return NULL;
+        bool isefficientlayer;
+        IssmDouble moulin_load,dt;
+        IssmDouble epl_active;
+
+        /*Initialize Element matrix*/
+        ElementVector* pe=new ElementVector(&node,1,this->parameters);
+
+        this->element->GetInputValue(&moulin_load,node,HydrologydcBasalMoulinInputEnum);
+        parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+        parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
+        // Test version input in EPL when active
+        if(isefficientlayer){
+                this->element->GetInputValue(&epl_active,node,HydrologydcMaskEplactiveNodeEnum);
+                if(epl_active){
+                        pe->values[0]=moulin_load*0.0;
+                }
+                else{
+                        pe->values[0]=moulin_load*dt;
+                }
+        }
+        else{
+                pe->values[0]=moulin_load*dt;
+        }
+        /*Clean up and return*/
+        return pe;
+ }
+/*}}}*/
+ElementVector* Moulin::CreatePVectorHydrologyDCEfficient(void){/*{{{*/
+
+        /*If this node is not the master node (belongs to another partition of the
+         * mesh), don't add the moulin input a second time*/
+        if(node->IsClone()) return NULL;
+        if(!this->node->IsActive()) return NULL;
+        IssmDouble moulin_load,dt;
+        ElementVector* pe=new ElementVector(&node,1,this->parameters);
+        
+        this->element->GetInputValue(&moulin_load,node,HydrologydcBasalMoulinInputEnum);
+        parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+        
+        pe->values[0]=moulin_load*dt;
+        /*Clean up and return*/
+        return pe;
+}
+/*}}}*/

issm/trunk-jpl/src/c/classes/Loads/Moulin.h

@@ -80 +80 @@
                 /*}}}*/
 
-                ElementVector* CreatePVectorMoulin(void);
+                ElementVector* CreatePVectorHydrologySommers(void);
+                ElementVector* CreatePVectorHydrologyDCInefficient(void);
+                ElementVector* CreatePVectorHydrologyDCEfficient(void);
 };
 

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

@@ -193 +193 @@
 /*}}}*/
 void  Pengrid::CreatePVector(Vector<IssmDouble>* pf){/*{{{*/
-
-        ElementVector* pe=NULL;
-        int analysis_type;
-        this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
-
-        switch(analysis_type){
-                
-        case HydrologyDCInefficientAnalysisEnum:
-                pe = CreatePVectorHydrologyDCInefficient();
-                break;
-        case HydrologyDCEfficientAnalysisEnum:
-                pe = CreatePVectorHydrologyDCEfficient();
-                break;
-        default:
-                /*No loads applied, do nothing: */
-                return;
-        }
-        if(pe){
-                pe->AddToGlobal(pf);
-                delete pe;
-        }
-
+        /*No loads applied, do nothing, originaly used for moulin input: */
+        return;
+        
 }
 /*}}}*/
@@ -415 +396 @@
         this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
 
-        if (analysis_type==StressbalanceAnalysisEnum){
+        switch(analysis_type){
+        case StressbalanceAnalysisEnum:
                 /*No penalty to check*/
                 return;
-        }
-        else if (analysis_type==ThermalAnalysisEnum){
+        case ThermalAnalysisEnum:
                 ConstraintActivateThermal(punstable);
-        }
-        else if (analysis_type==MeltingAnalysisEnum){
+                break;
+        case MeltingAnalysisEnum:
                 /*No penalty to check*/
                 return;
-        }
-        else if (analysis_type==HydrologyDCInefficientAnalysisEnum){
+        case HydrologyDCInefficientAnalysisEnum:
                 ConstraintActivateHydrologyDCInefficient(punstable);
-                return;
-        }
-        else{
+                break;
+        default:
                 _error_("analysis: " << EnumToStringx(analysis_type) << " not supported yet");
         }
-
 }
 /*}}}*/
@@ -559 +537 @@
 }
 /*}}}*/
-ElementVector* Pengrid::CreatePVectorHydrologyDCInefficient(void){/*{{{*/
-
-        bool active_element,isefficientlayer;
-        IssmDouble moulin_load,dt;
-        IssmDouble epl_active;
-
-        /*Initialize Element matrix*/
-        ElementVector* pe=new ElementVector(&node,1,this->parameters);
-
-        this->element->GetInputValue(&moulin_load,node,HydrologydcBasalMoulinInputEnum);
-        parameters->FindParam(&dt,TimesteppingTimeStepEnum);
-        parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
-
-
-        // Test version input in EPL when active
-        if(isefficientlayer){
-                this->element->GetInputValue(&epl_active,node,HydrologydcMaskEplactiveNodeEnum);
-                if(epl_active){
-                        pe->values[0]=moulin_load*0.0;
-                }
-                else{
-                        if(dt!=0.0) pe->values[0]=moulin_load*dt;
-                }
-        }
-        else{
-                        if(dt!=0.0) pe->values[0]=moulin_load*dt;
-        }
-
-        //Initial treatment, allways input in sediment
-        /* if(dt!=0.0) pe->values[0]=moulin_load*dt; */
-
-        /*Clean up and return*/
-        return pe;
- }
-/*}}}*/
-ElementVector* Pengrid::CreatePVectorHydrologyDCEfficient(void){/*{{{*/
-
-        IssmDouble moulin_load,dt;
-
-        /*Initialize Element matrix*/
-        if(!this->active) return NULL;
-        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;
- }
-/*}}}*/
 ElementMatrix* Pengrid::PenaltyCreateKMatrixHydrologyDCInefficient(IssmDouble kmax){/*{{{*/
         IssmDouble    penalty_factor;

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

@@ -87 +87 @@
                 void           ConstraintActivateHydrologyDCInefficient(int* punstable);
                 void           ConstraintActivateThermal(int* punstable);
-                ElementVector* CreatePVectorHydrologyDCInefficient(void);
-                ElementVector* CreatePVectorHydrologyDCEfficient(void);
                 ElementMatrix* PenaltyCreateKMatrixHydrologyDCInefficient(IssmDouble kmax);
                 ElementMatrix* PenaltyCreateKMatrixMelting(IssmDouble kmax);