Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 21212) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 21213) @@ -105,4 +105,28 @@ 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;inumberofvertices;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(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"); }/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 21212) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 21213) @@ -140,4 +140,5 @@ 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,4 +146,5 @@ 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)); } } Index: /issm/trunk-jpl/src/c/classes/Loads/Moulin.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Moulin.cpp (revision 21212) +++ /issm/trunk-jpl/src/c/classes/Loads/Moulin.cpp (revision 21213) @@ -180,5 +180,11 @@ case HydrologySommersAnalysisEnum: - pe = this->CreatePVectorMoulin(); + pe = this->CreatePVectorHydrologySommers(); + break; + case HydrologyDCInefficientAnalysisEnum: + pe = CreatePVectorHydrologyDCInefficient(); + break; + case HydrologyDCEfficientAnalysisEnum: + pe = CreatePVectorHydrologyDCEfficient(); break; default: @@ -331,5 +337,5 @@ /*}}}*/ -ElementVector* Moulin::CreatePVectorMoulin(void){/*{{{*/ +ElementVector* Moulin::CreatePVectorHydrologySommers(void){/*{{{*/ /*If this node is not the master node (belongs to another partition of the @@ -349,2 +355,51 @@ } /*}}}*/ +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; +} +/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Loads/Moulin.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Moulin.h (revision 21212) +++ /issm/trunk-jpl/src/c/classes/Loads/Moulin.h (revision 21213) @@ -80,5 +80,7 @@ /*}}}*/ - ElementVector* CreatePVectorMoulin(void); + ElementVector* CreatePVectorHydrologySommers(void); + ElementVector* CreatePVectorHydrologyDCInefficient(void); + ElementVector* CreatePVectorHydrologyDCEfficient(void); }; Index: /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp (revision 21212) +++ /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp (revision 21213) @@ -193,26 +193,7 @@ /*}}}*/ void Pengrid::CreatePVector(Vector* 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,23 +396,20 @@ 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,56 +537,4 @@ } /*}}}*/ -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; Index: /issm/trunk-jpl/src/c/classes/Loads/Pengrid.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Pengrid.h (revision 21212) +++ /issm/trunk-jpl/src/c/classes/Loads/Pengrid.h (revision 21213) @@ -87,6 +87,4 @@ void ConstraintActivateHydrologyDCInefficient(int* punstable); void ConstraintActivateThermal(int* punstable); - ElementVector* CreatePVectorHydrologyDCInefficient(void); - ElementVector* CreatePVectorHydrologyDCEfficient(void); ElementMatrix* PenaltyCreateKMatrixHydrologyDCInefficient(IssmDouble kmax); ElementMatrix* PenaltyCreateKMatrixMelting(IssmDouble kmax);