Changeset 4251

Timestamp:

Feb 9, 2007, 4:23:15 PM (18 years ago)

Author:

ole

Message:

Pulled speed computation out into a separate function that can be used by all flux functions.

File:

• anuga_core/source/anuga/shallow_water/shallow_water_ext.c (modified) (3 diffs)

anuga_core/source/anuga/shallow_water/shallow_water_ext.c

@@ -125 +125 @@
 }
 
+
+// Function to obtain speed from momentum and depth.
+// This is used by flux functions
+// Input parameters uh and h may be modified by this function.
+double _compute_speed(double *uh, 
+                      double *h, 
+                      double epsilon, 
+                      double h0) {
+  
+  double u;
+  
+  if (*h < epsilon) {
+    *h = 0.0;  //Could have been negative
+    u = 0.0;
+  } else {
+    u = *uh/(*h + h0/ *h);
+  }
+  
+  return u;
+}
+
 // Computational function for flux computation (using stage w=z+h)
 int flux_function_central(double *q_left, double *q_right,
@@ -167 +188 @@
 
   //Compute speeds in x-direction
-  w_left = q_left_copy[0];              // h+z
+  w_left = q_left_copy[0];          
   h_left = w_left-z;
   uh_left = q_left_copy[1];
-
-  if (h_left < epsilon) {
-    h_left = 0.0;  //Could have been negative
-    u_left = 0.0;
-  } else {
-    u_left = uh_left/(h_left + h0/h_left);
-  }
+  u_left =_compute_speed(&uh_left, &h_left, epsilon, h0);
 
   w_right = q_right_copy[0];
   h_right = w_right-z;
   uh_right = q_right_copy[1];
-
-  if (h_right < epsilon) {
-    h_right = 0.0; //Could have been negative
-    u_right = 0.0;
-  } else {
-    u_right = uh_right/(h_right + h0/h_right);
-  }
+  u_right =_compute_speed(&uh_right, &h_right, epsilon, h0);  
 
   //Momentum in y-direction
@@ -288 +297 @@
 
   //Compute speeds in x-direction
-  w_left = q_left_copy[0];              // h+z
+  w_left = q_left_copy[0];          
   h_left = w_left-z;
   uh_left = q_left_copy[1];
-
-  if (h_left < epsilon) {
-    h_left = 0.0;  //Could have been negative
-    u_left = 0.0;
-  } else {
-    u_left = uh_left/(h_left + h0/h_left);
-  }
+  u_left =_compute_speed(&uh_left, &h_left, epsilon, h0);
 
   w_right = q_right_copy[0];
   h_right = w_right-z;
   uh_right = q_right_copy[1];
-
-  if (h_right < epsilon) {
-    h_right = 0.0; //Could have been negative
-    u_right = 0.0;
-  } else {
-    u_right = uh_right/(h_right + h0/h_right);
-  }
+  u_right =_compute_speed(&uh_right, &h_right, epsilon, h0);  
+
 
   //Momentum in y-direction