Changeset 8017 for trunk/anuga_core/source/anuga/utilities

Timestamp:

Sep 16, 2010, 10:38:49 AM (14 years ago)

Author:

steve

Message:

Added code to pickup if elevation is discontinuous. compute_fluxes now produces an error in this case

Location:

trunk/anuga_core/source/anuga/utilities

Files:

• polygon_ext.c (modified) (2 diffs)
• util_ext.h (modified) (1 diff)

trunk/anuga_core/source/anuga/utilities/polygon_ext.c

@@ -20 +20 @@
 #include "util_ext.h"
 
+#define YES 1
+#define NO 0
+
 
 double dist(double x,
@@ -89 +92 @@
 
 
+//  public domain function by Darel Rex Finley, 2006
+//  http://www.alienryderflex.com/intersect/
+//
+//  Determines the intersection point of the line segment defined by points A and B
+//  with the line segment defined by points C and D.
+//
+//  Returns YES if the intersection point was found, and stores that point in X,Y.
+//  Returns NO if there is no determinable intersection point, in which case X,Y will
+//  be unmodified.
+
+int __lineSegmentIntersection(
+        double Ax, double Ay,
+        double Bx, double By,
+        double Cx, double Cy,
+        double Dx, double Dy,
+        double *X, double *Y) {
+
+    double distAB, theCos, theSin, newX, ABpos;
+
+    //  Fail if either line segment is zero-length.
+    if ( (Ax == Bx && Ay == By) || (Cx == Dx && Cy == Dy) ) return NO ;
+
+    //  Fail if the segments share an end-point.
+    if ( (Ax == Cx && Ay == Cy) || (Bx == Cx && By == Cy)
+            || (Ax == Dx && Ay == Dy) || (Bx == Dx && By == Dy) ) {
+        return NO;
+    }
+
+    //  (1) Translate the system so that point A is on the origin.
+    Bx -= Ax;
+    By -= Ay;
+    Cx -= Ax;
+    Cy -= Ay;
+    Dx -= Ax;
+    Dy -= Ay;
+
+    //  Discover the length of segment A-B.
+    distAB = sqrt(Bx * Bx + By * By);
+
+    //  (2) Rotate the system so that point B is on the positive X axis.
+    theCos = Bx / distAB;
+    theSin = By / distAB;
+    newX = Cx * theCos + Cy*theSin;
+    Cy = Cy * theCos - Cx*theSin;
+    Cx = newX;
+    newX = Dx * theCos + Dy*theSin;
+    Dy = Dy * theCos - Dx*theSin;
+    Dx = newX;
+
+    //  Fail if segment C-D doesn't cross line A-B.
+    if ( (Cy < 0. && Dy < 0.) || (Cy >= 0. && Dy >= 0.) ) return NO;
+
+    //  (3) Discover the position of the intersection point along line A-B.
+    ABpos = Dx + (Cx - Dx) * Dy / (Dy - Cy);
+
+    //  Fail if segment C-D crosses line A-B outside of segment A-B.
+    if (ABpos < 0. || ABpos > distAB) return NO;
+
+    //  (4) Apply the discovered position to line A-B in the original coordinate system.
+    *X = Ax + ABpos*theCos;
+    *Y = Ay + ABpos*theSin;
+
+    //  Success.
+    return YES;
+}
 
 /*

trunk/anuga_core/source/anuga/utilities/util_ext.h

@@ -14 +14 @@
 #include "math.h"
 
+#include <stdio.h>
+#define STRINGIFY(x) #x
+#define TOSTRING(x) STRINGIFY(x)
+#define AT __FILE__ ":" TOSTRING(__LINE__)
+#define P_ERROR_BUFFER_SIZE 65
+void report_python_error(const char *location, const char *msg)
+{
+
+    char buf[P_ERROR_BUFFER_SIZE];
+    int n;
+    n = snprintf(buf, P_ERROR_BUFFER_SIZE, "Error at %s: %s\n", location, msg);
+
+    PyErr_SetString(PyExc_RuntimeError, buf);
+}
+
+
 
 double max(double x, double y) {