Changeset 15020 for issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

Timestamp:

05/13/13 16:45:14 (12 years ago)

Author:

bdef

Message:

NEW: Introduction of the transfer term in the Pvectors

File:

• issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (modified) (9 diffs)

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

@@ -2078 +2078 @@
                                 name==SedimentHeadOldEnum ||
                                 name==SedimentHeadEnum ||
+                                name==WaterTransferEnum ||
                                 name==BasisIntegralEnum ||
                                 name==QmuVxEnum ||
@@ -6150 +6151 @@
         IssmDouble xyz_list[NUMVERTICES][3];
         IssmDouble dt,scalar,water_head;
-        IssmDouble water_load;
+        IssmDouble water_load,transfer;
         IssmDouble sediment_storing;
         IssmDouble basis[numdof];
@@ -6163 +6164 @@
         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
         Input* water_input=inputs->GetInput(BasalforcingsMeltingRateEnum);  _assert_(water_input);
+        Input* transfer_input=inputs->GetInput(WaterTransferEnum);  _assert_(transfer_input);
         Input* old_wh_input=NULL; 
 
@@ -6181 +6183 @@
                 /*Loading term*/
                 water_input->GetInputValue(&water_load,gauss);
-                scalar = Jdet*gauss->weight*water_load;
+                transfer_input->GetInputValue(&transfer,gauss);
+                scalar = Jdet*gauss->weight*(water_load-transfer);
                 if(reCast<bool,IssmDouble>(dt)) scalar = scalar*dt;
                 for(int i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
@@ -6207 +6210 @@
         IssmDouble xyz_list[NUMVERTICES][3];
         IssmDouble dt,scalar,water_head;
-        IssmDouble residual_load;
+        IssmDouble residual_load,transfer;
         IssmDouble epl_storing;
         IssmDouble basis[numdof];
@@ -6220 +6223 @@
         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
         Input* residual_input=inputs->GetInput(SedimentHeadResidualEnum);  _assert_(residual_input);
+        Input* transfer_input=inputs->GetInput(WaterTransferEnum);  _assert_(transfer_input);
         Input* old_wh_input=NULL; 
 
@@ -6237 +6241 @@
                 /*Loading term*/
                 residual_input->GetInputValue(&residual_load,gauss);
-                scalar = Jdet*gauss->weight*residual_load;
+                transfer_input->GetInputValue(&transfer,gauss);
+                scalar = Jdet*gauss->weight*(residual_load+transfer);
                 if(reCast<bool,IssmDouble>(dt)) scalar = scalar*dt;
                 for(int i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
@@ -6350 +6355 @@
                 for(int i=0;i<NUMVERTICES;i++){
                         this->GetHydrologyDCInefficientHmax(&h_max,nodes[i]);
-                        if(values[i]>h_max)
-                                residual[i]=kappa*(values[i]-h_max)/(dt*intbasis[i]);
+                        if(values[i]>h_max){
+                                residual[i]=kappa*(values[i]-h_max)/intbasis[i];
+                                if(reCast<bool,IssmDouble>(dt))residual[i]=residual[i]/dt;
+                        }
                         else
                                 residual[i]=0.0;
                 }
-        }       
-
+        }
         /*Add input to the element: */
 
@@ -6516 +6522 @@
                         }
                 }
+                break;
         default:
                 _error_("no case higher than 1 for the Transfer method");