Changeset 1486 for inundation/ga/storm_surge/pyvolution/util_ext.h

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

File:

• inundation/ga/storm_surge/pyvolution/util_ext.h (modified) (1 diff)

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;
 }