Changeset 1486

Timestamp:

May 31, 2005, 2:24:14 PM (20 years ago)

Author:

ole

Message:

Applyed new formula for computing gradients for triangles that have only one neighbour (thanks to Matt for pointing out the bug) and added unit tests as appropriate.
See util_ext.c for documentation of the derivation

Location:

inundation/ga/storm_surge/pyvolution

Files:

• quantity.py (modified) (1 diff)
• quantity_ext.c (modified) (2 diffs)
• test_quantity.py (modified) (4 diffs)
• test_util.py (modified) (2 diffs)
• util.py (modified) (2 diffs)
• util_ext.c (modified) (2 diffs)
• util_ext.h (modified) (1 diff)

inundation/ga/storm_surge/pyvolution/quantity.py

@@ -792 +792 @@
 
             #Gradient
-            det = x0*y1 - x1*y0
-            if det != 0.0:
-                a[k] = (y1*q0 - y0*q1)/det
-                b[k] = (x0*q1 - x1*q0)/det
-
+            a[k], b[k] = gradient2(x0, y0, x1, y1, q0, q1)          
         else:
             #No true neighbours -

inundation/ga/storm_surge/pyvolution/quantity_ext.c

@@ -29 +29 @@
 
   int i, k, k0, k1, k2, index3;
-  double x0, x1, x2, y0, y1, y2, q0, q1, q2, det;
+  double x0, x1, x2, y0, y1, y2, q0, q1, q2; //, det;
 
 
@@ -80 +80 @@
       x1 = centroids[k1*2]; y1 = centroids[k1*2+1];
 
-      //Gradient
-      det = x0*y1 - x1*y0;
-      if (det != 0.0) {
-        a[k] = (y1*q0 - y0*q1)/det;
-        b[k] = (x0*q1 - x1*q0)/det;
-      }
+      //Two point gradient
+      _gradient2(x0, y0, x1, y1, q0, q1, &a[k], &b[k]);      
+      
+      
+      //Old (wrong code)
+      //det = x0*y1 - x1*y0;
+      //if (det != 0.0) {
+      //        a[k] = (y1*q0 - y0*q1)/det;
+      //        b[k] = (x0*q1 - x1*q0)/det;
+      //}
     }
     //    else:

inundation/ga/storm_surge/pyvolution/test_quantity.py

@@ -84 +84 @@
         quantity.extrapolate_second_order()
 
-        assert allclose(quantity.vertex_values, [[2., 2., 2.],
-                                               [3.+2./3, 6.+2./3, 4.+2./3],
-                                               [7.5, 0.5, 1.],
-                                               [-5, -2.5, 7.5]])
-
+        #print quantity.vertex_values
+        #assert allclose(quantity.vertex_values, [[2., 2., 2.],
+        #                                         [3.+2./3, 6.+2./3, 4.+2./3],
+        #                                         [7.5, 0.5, 1.],
+        #                                         [-5, -2.5, 7.5]])
+
+        assert allclose(quantity.vertex_values[1,:],[3.+2./3, 6.+2./3, 4.+2./3])
+        #FIXME: Work out the others
+        
+                                         
         #print quantity.edge_values
         assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
@@ -267 +272 @@
         quantity = Conserved_quantity(self.mesh4)
 
-        #Set up for a gradient of (3,0) at mid triangle
-        quantity.set_values([2.0, 4.0, 8.0, 2.0],
+        #Set up for a gradient of (2,0) at mid triangle
+        quantity.set_values([2.0, 4.0, 6.0, 2.0],
                             location = 'centroids')
-
 
 
@@ -276 +280 @@
         a, b = quantity.compute_gradients()
 
-        #gradient bewteen t0 and t1 is undefined as det==0
-        assert a[0] == 0.0
-        assert b[0] == 0.0
-        #The others are OK
-        for i in range(1,4):
-            assert a[i] == 3.0
-            assert b[i] == 0.0
-
-
+        #print self.mesh4.centroid_coordinates
+        #print a, b
+
+        #The central triangle (1)
+        #(using standard gradient based on neigbours controid values)
+        assert a[1] == 2.0
+        assert b[1] == 0.0
+
+
+        #Left triangle (0) using two point gradient
+        #q0 = q1 + a*(x0-x1) + b*(y0-y1)  <=>
+        #2  = 4  + a*(-2/3)  + b*(-2/3)
+        assert a[0] + b[0] == 3
+        #From orthogonality (a*(y0-y1) + b*(x0-x1) == 0)        
+        assert a[0] - b[0] == 0 
+
+
+        #Right triangle (2) using two point gradient
+        #q2 = q1 + a*(x2-x1) + b*(y2-y1)  <=>
+        #6  = 4  + a*(4/3)  + b*(-2/3)        
+        assert 2*a[2] - b[2] == 3
+        #From orthogonality (a*(y1-y2) + b*(x2-x1) == 0)
+        assert a[2] + 2*b[2] == 0 
+
+
+        #Top triangle (3) using two point gradient
+        #q3 = q1 + a*(x3-x1) + b*(y3-y1)  <=>
+        #2  = 4  + a*(-2/3)  + b*(4/3)        
+        assert a[3] - 2*b[3] == 3
+        #From orthogonality (a*(y1-y3) + b*(x3-x1) == 0)        
+        assert 2*a[3] + b[3] == 0 
+
+
+
+        #print a, b
         quantity.extrapolate_second_order()
 
-        #print quantity.vertex_values
-        assert allclose(quantity.vertex_values, [[2., 2.,  2.],
-                                                 [0., 6.,  6.],
-                                                 [6., 6., 12.],
-                                                 [0., 0.,  6.]])
+        #Apply q(x,y) = qc + a*(x-xc) + b*(y-yc)
+        assert allclose(quantity.vertex_values[0,:], [3., 0.,  3.])
+        assert allclose(quantity.vertex_values[1,:], [4./3, 16./3,  16./3])
+
+
+        #a = 1.2, b=-0.6
+        #q(4,0) = 6 + a*(4 - 8/3) + b*(-2/3)
+        assert allclose(quantity.vertex_values[2,2], 8)       
 
 
@@ -304 +337 @@
         a, b = quantity.compute_gradients()
 
-        #gradient bewteen t0 and t1 is undefined as det==0
-        assert a[0] == 0.0
-        assert b[0] == 0.0
-        #The others are OK
-        for i in range(1,4):
-            assert allclose(a[i], 3.0)
-            assert allclose(b[i], 1.0)
-
+        #print a, b
+
+        assert a[1] == 3.0
+        assert b[1] == 1.0
+        
+        #Work out the others
+        
         quantity.extrapolate_second_order()
 

inundation/ga/storm_surge/pyvolution/test_util.py

@@ -31 +31 @@
         assert zy == 0.0
 
-
-
-
-    def test_gradient2(self):
-
-        from util import gradient
+    def test_gradient_more(self):
         x0 = 2.0/3; y0 = 2.0/3
         x1=  8.0/3; y1 = 2.0/3
@@ -51 +46 @@
         assert abs(b - 1.0) < epsilon
 
+
+    def test_gradient2(self):
+        """Test two-point gradient
+        """
+        
+        x0 = 5.0; y0 = 5.0; z0 = 10.0
+        x1 = 8.0; y1 = 2.0; z1 = 1.0
+        x2 = 8.0; y2 = 8.0; z2 = 10.0
+
+        #Reference
+        zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)
+        a, b = gradient2(x0, y0, x1, y1, z0, z1)
+
+        assert zx == a
+        assert zy == b
+
+        z2_computed = z0 + a*(x2-x0) + b*(y2-y0)
+        assert z2_computed == z2
+        
+    def test_gradient2_more(self):
+        """Test two-point gradient more
+        """
+        x0 = 2.0; y0 = 2.0
+        x1 = 8.0; y1 = 3.0
+        x2 = 1.0; y2 = 8.0
+
+        q0 = 2.0
+        q1 = 8.0
+        q2 = q0
+
+        #Gradient of fitted pwl surface
+        a_ref, b_ref = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
+        a, b = gradient2(x0, y0, x1, y1, q0, q1)        
+
+        assert a == a_ref
+        assert b == b_ref
 
 

inundation/ga/storm_surge/pyvolution/util.py

@@ -1401 +1401 @@
 
 
+def gradient2_python(x0, y0, x1, y1, q0, q1):
+    """Compute radient based on two points and enforce zero gradient
+    in the direction orthogonal to (x1-x0), (y1-y0) 
+    """
+
+    #Old code
+    #det = x0*y1 - x1*y0
+    #if det != 0.0:
+    #    a = (y1*q0 - y0*q1)/det
+    #    b = (x0*q1 - x1*q0)/det
+
+    #Correct code (ON)
+    det = (x1-x0)**2 + (y1-y0)**2
+    if det != 0.0:
+        a = (x1-x0)*(q1-q0)/det
+        b = (y1-y0)*(q1-q0)/det
+        
+    return a, b        
+
 
 ##############################################
@@ -1407 +1426 @@
 import compile
 if compile.can_use_C_extension('util_ext.c'):
-    from util_ext import gradient, point_on_line
+    from util_ext import gradient, gradient2, point_on_line
     separate_points_by_polygon = separate_points_by_polygon_c
 else:
     gradient = gradient_python
+    gradient2 = gradient2_python    
 
 

inundation/ga/storm_surge/pyvolution/util_ext.c

@@ -272 +272 @@
 }
 
+PyObject *gradient2(PyObject *self, PyObject *args) {
+  //
+  // a,b = gradient2(x0, y0, x1, y1, q0, q1)
+  //
+
+  double x0, y0, x1, y1, q0, q1, a, b;
+  PyObject *result;
+
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "dddddd", &x0, &y0, &x1, &y1, &q0, &q1))
+    return NULL;
+
+
+  // Call underlying routine
+  _gradient2(x0, y0, x1, y1, q0, q1, &a, &b);
+
+  // Return values a and b
+  result = Py_BuildValue("dd", a, b);
+  return result;
+}
 
 // Method table for python module
@@ -282 +302 @@
   //{"rotate", (PyCFunction)rotate, METH_VARARGS | METH_KEYWORDS, "Print out"},
   {"gradient", gradient, METH_VARARGS, "Print out"},
+  {"gradient2", gradient2, METH_VARARGS, "Print out"},  
   {"point_on_line", point_on_line, METH_VARARGS, "Print out"},
   //{"inside_polygon", inside_polygon, METH_VARARGS, "Print out"},

inundation/ga/storm_surge/pyvolution/util_ext.h

@@ -47 +47 @@
   *b /= det;
 
+  return 0;
+}
+
+
+int _gradient2(double x0, double y0, 
+               double x1, double y1, 
+               double q0, double q1, 
+               double *a, double *b) {
+  /*Compute gradient (a,b) between two points (x0,y0) and (x1,y1) 
+  with values q0 and q1 such that the plane is constant in the direction 
+  orthogonal to (x1-x0, y1-y0).
+  
+  Extrapolation formula
+    q(x,y) = q0 + a*(x-x0) + b*(y-y0)                    (1)
+  
+  Substituting the known values for q1 into (1) yields an 
+  under determined  equation for a and b 
+      q1-q0 = a*(x1-x0) + b*(y1-y0)                      (2)
+      
+      
+  Now add the additional requirement that the gradient in the direction 
+  orthogonal to (x1-x0, y1-y0) should be zero. The orthogonal direction 
+  is given by the vector (y0-y1, x1-x0).
+  
+  Define the point (x2, y2) = (x0 + y0-y1, y0 + x1-x0) on the orthognal line. 
+  Then we know that the corresponding value q2 should be equal to q0 in order 
+  to obtain the zero gradient, hence applying (1) again    
+      q0 = q2 = q(x2, y2) = q0 + a*(x2-x0) + b*(y2-y0)
+                          = q0 + a*(x0 + y0-y1-x0) + b*(y0 + x1-x0 - y0)
+                          = q0 + a*(y0-y1) + b*(x1-x0)
+                          
+  leads to the orthogonality constraint
+     a*(y0-y1) + b*(x1-x0) = 0                           (3) 
+     
+  which closes the system and yields
+  
+  /               \  /   \   /       \  
+  |  x1-x0  y1-y0 |  | a |   | q1-q0 |
+  |               |  |   | = |       | 
+  |  y0-y1  x1-x0 |  | b |   |   0   |
+  \               /  \   /   \       /
+   
+  which is solved using the standard determinant technique    
+      
+  */
+
+  double det, xx, yy, qq;
+  
+  xx = x1-x0;
+  yy = y1-y0;
+  qq = q1-q0;
+    
+  det = xx*xx + yy*yy;  //FIXME  catch det == 0
+  *a = xx*qq/det;
+  *b = yy*qq/det;
+        
   return 0;
 }