Changeset 8017
 Timestamp:
 Sep 16, 2010, 10:38:49 AM (9 years ago)
 Location:
 trunk/anuga_core/source/anuga
 Files:

 5 edited
Legend:
 Unmodified
 Added
 Removed

trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py
r7737 r8017 1309 1309 self.set_vertex_values(A, use_cache=use_cache, verbose=verbose) 1310 1310 1311 1312 1311 1313 ############################################################################ 1312 1314 # Methods for outputting model results 
trunk/anuga_core/source/anuga/shallow_water/shallow_water_domain.py
r7967 r8017 243 243 244 244 245 def set_ all_betas(self, beta):246 """Shorthand to assign one constant value [0, 1] to all limiters.245 def set_beta(self, beta): 246 """Shorthand to assign one constant value [0,2] to all limiters. 247 247 0 Corresponds to first order, where as larger values make use of 248 248 the second order scheme. 
trunk/anuga_core/source/anuga/shallow_water/shallow_water_ext.c
r8009 r8017 22 22 // Shared code snippets 23 23 #include "util_ext.h" 24 25 26 24 27 25 28 … … 1045 1048 &normal, &ql, &qr, &zl, &zr, &edgeflux, 1046 1049 &epsilon, &g, &H0)) { 1047 PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: flux_function_centralcould not parse input arguments");1048 return NULL;1050 report_python_error(AT, "could not parse input arguments"); 1051 return NULL; 1049 1052 } 1050 1053 … … 1090 1093 &xmom, &ymom)) { 1091 1094 //&epsilon)) { 1092 PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: gravitycould not parse input arguments");1093 return NULL;1095 report_python_error(AT, "could not parse input arguments"); 1096 return NULL; 1094 1097 } 1095 1098 … … 1154 1157 if (!PyArg_ParseTuple(args, "ddOOOOOOOO", 1155 1158 &g, &eps, &x, &w, &uh, &vh, &z, &eta, &xmom, &ymom)) { 1156 PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: manning_friction_newcould not parse input arguments");1157 return NULL;1159 report_python_error(AT, "could not parse input arguments"); 1160 return NULL; 1158 1161 } 1159 1162 … … 1196 1199 if (!PyArg_ParseTuple(args, "ddOOOOOOOO", 1197 1200 &g, &eps, &x, &w, &uh, &vh, &z, &chezy, &xmom, &ymom)) { 1198 PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: chezy_frictioncould not parse input arguments");1199 return NULL;1201 report_python_error(AT, "could not parse input arguments"); 1202 return NULL; 1200 1203 } 1201 1204 … … 1238 1241 if (!PyArg_ParseTuple(args, "ddOOOOOOO", 1239 1242 &g, &eps, &w, &uh, &vh, &z, &eta, &xmom, &ymom)) { 1240 PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: manning_friction_oldcould not parse input arguments");1241 return NULL;1243 report_python_error(AT, "could not parse input arguments"); 1244 return NULL; 1242 1245 } 1243 1246 … … 1438 1441 if (area2 <= 0) 1439 1442 { 1440 PyErr_SetString(PyExc_RuntimeError, 1441 "shallow_water_ext.c: Negative triangle area encountered. This might be due to a triangle being flat or it's vertices not order anticlockwise."); 1442 1443 return 1; 1443 report_python_error(AT, "Negative triangle area"); 1444 return 1; 1444 1445 } 1445 1446 … … 1630 1631 { 1631 1632 // If we didn't find an internal neighbour 1632 PyErr_SetString(PyExc_RuntimeError, 1633 "shallow_water_ext.c: Internal neighbour not found"); 1633 report_python_error(AT, "Internal neighbour not found"); 1634 1634 return 1; 1635 1635 } … … 1866 1866 &elevation_vertex_values, 1867 1867 &optimise_dry_cells)) { 1868 1869 PyErr_SetString(PyExc_RuntimeError, 1870 "Input arguments to extrapolate_second_order_sw failed"); 1868 1869 report_python_error(AT, "could not parse input arguments"); 1871 1870 return NULL; 1872 1871 } … … 1957 1956 if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OOi", argnames, 1958 1957 &Q, &Normal, &direction)) { 1959 PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: rotatecould not parse input arguments");1958 report_python_error(AT, "could not parse input arguments"); 1960 1959 return NULL; 1961 1960 } … … 1969 1968 1970 1969 if (normal > dimensions[0] != 2) { 1971 PyErr_SetString(PyExc_RuntimeError, "Normal vector must have 2 components");1970 report_python_error(AT, "Normal vector must have 2 components"); 1972 1971 return NULL; 1973 1972 } … … 2000 1999 2001 2000 // Computational function for flux computation 2001 2002 2002 double _compute_fluxes_central(int number_of_elements, 2003 double timestep, 2004 double epsilon, 2005 double H0, 2006 double g, 2007 long* neighbours, 2008 long* neighbour_edges, 2009 double* normals, 2010 double* edgelengths, 2011 double* radii, 2012 double* areas, 2013 long* tri_full_flag, 2014 double* stage_edge_values, 2015 double* xmom_edge_values, 2016 double* ymom_edge_values, 2017 double* bed_edge_values, 2018 double* stage_boundary_values, 2019 double* xmom_boundary_values, 2020 double* ymom_boundary_values, 2021 double* stage_explicit_update, 2022 double* xmom_explicit_update, 2023 double* ymom_explicit_update, 2024 long* already_computed_flux, 2025 double* max_speed_array, 2026 int optimise_dry_cells) 2027 { 2028 // Local variables 2029 double max_speed, length, inv_area, zl, zr; 2030 double h0 = H0*H0; // This ensures a good balance when h approaches H0. 2031 2032 double limiting_threshold = 10*H0; // Avoid applying limiter below this 2033 // threshold for performance reasons. 2034 // See ANUGA manual under flux limiting 2035 int k, i, m, n; 2036 int ki, nm=0, ki2; // Index shorthands 2037 2038 2039 // Workspace (making them static actually made function slightly slower (Ole)) 2040 double ql[3], qr[3], edgeflux[3]; // Work array for summing up fluxes 2041 2042 static long call = 1; // Static local variable flagging already computed flux 2043 2044 // Start computation 2045 call++; // Flag 'id' of flux calculation for this timestep 2046 2047 // Set explicit_update to zero for all conserved_quantities. 2048 // This assumes compute_fluxes called before forcing terms 2049 memset((char*) stage_explicit_update, 0, number_of_elements*sizeof(double)); 2050 memset((char*) xmom_explicit_update, 0, number_of_elements*sizeof(double)); 2051 memset((char*) ymom_explicit_update, 0, number_of_elements*sizeof(double)); 2052 2053 // For all triangles 2054 for (k = 0; k < number_of_elements; k++) 2055 { 2056 // Loop through neighbours and compute edge flux for each 2057 for (i = 0; i < 3; i++) 2058 { 2059 ki = k*3 + i; // Linear index to edge i of triangle k 2060 2061 if (already_computed_flux[ki] == call) 2062 { 2063 // We've already computed the flux across this edge 2064 continue; 2065 } 2066 2067 // Get left hand side values from triangle k, edge i 2068 ql[0] = stage_edge_values[ki]; 2069 ql[1] = xmom_edge_values[ki]; 2070 ql[2] = ymom_edge_values[ki]; 2071 zl = bed_edge_values[ki]; 2072 2073 // Get right hand side values either from neighbouring triangle 2074 // or from boundary array (Quantities at neighbour on nearest face). 2075 n = neighbours[ki]; 2076 if (n < 0) 2077 { 2078 // Neighbour is a boundary condition 2079 m = n  1; // Convert negative flag to boundary index 2080 2081 qr[0] = stage_boundary_values[m]; 2082 qr[1] = xmom_boundary_values[m]; 2083 qr[2] = ymom_boundary_values[m]; 2084 zr = zl; // Extend bed elevation to boundary 2085 } 2086 else 2087 { 2088 // Neighbour is a real triangle 2089 m = neighbour_edges[ki]; 2090 nm = n*3 + m; // Linear index (triangle n, edge m) 2091 2092 qr[0] = stage_edge_values[nm]; 2093 qr[1] = xmom_edge_values[nm]; 2094 qr[2] = ymom_edge_values[nm]; 2095 zr = bed_edge_values[nm]; 2096 } 2097 2098 // Now we have values for this edge  both from left and right side. 2099 2100 if (optimise_dry_cells) 2101 { 2102 // Check if flux calculation is necessary across this edge 2103 // This check will exclude dry cells. 2104 // This will also optimise cases where zl != zr as 2105 // long as both are dry 2106 2107 if (fabs(ql[0]  zl) < epsilon && 2108 fabs(qr[0]  zr) < epsilon) 2109 { 2110 // Cell boundary is dry 2111 2112 already_computed_flux[ki] = call; // #k Done 2113 if (n >= 0) 2114 { 2115 already_computed_flux[nm] = call; // #n Done 2116 } 2117 2118 max_speed = 0.0; 2119 continue; 2120 } 2121 } 2122 2123 // Outward pointing normal vector (domain.normals[k, 2*i:2*i+2]) 2124 ki2 = 2*ki; //k*6 + i*2 2125 2126 // Edge flux computation (triangle k, edge i) 2127 _flux_function_central(ql, qr, zl, zr, 2128 normals[ki2], normals[ki2+1], 2129 epsilon, h0, limiting_threshold, g, 2130 edgeflux, &max_speed); 2131 2132 2133 // Multiply edgeflux by edgelength 2134 length = edgelengths[ki]; 2135 edgeflux[0] *= length; 2136 edgeflux[1] *= length; 2137 edgeflux[2] *= length; 2138 2139 2140 // Update triangle k with flux from edge i 2141 stage_explicit_update[k] = edgeflux[0]; 2142 xmom_explicit_update[k] = edgeflux[1]; 2143 ymom_explicit_update[k] = edgeflux[2]; 2144 2145 already_computed_flux[ki] = call; // #k Done 2146 2147 2148 // Update neighbour n with same flux but reversed sign 2149 if (n >= 0) 2150 { 2151 stage_explicit_update[n] += edgeflux[0]; 2152 xmom_explicit_update[n] += edgeflux[1]; 2153 ymom_explicit_update[n] += edgeflux[2]; 2154 2155 already_computed_flux[nm] = call; // #n Done 2156 } 2157 2158 // Update timestep based on edge i and possibly neighbour n 2159 if (tri_full_flag[k] == 1) 2160 { 2161 if (max_speed > epsilon) 2162 { 2163 // Apply CFL condition for triangles joining this edge (triangle k and triangle n) 2164 2165 // CFL for triangle k 2166 timestep = min(timestep, radii[k]/max_speed); 2167 2168 if (n >= 0) 2169 { 2170 // Apply CFL condition for neigbour n (which is on the ith edge of triangle k) 2171 timestep = min(timestep, radii[n]/max_speed); 2172 } 2173 2174 // Ted Rigby's suggested less conservative version 2175 //if (n>=0) { 2176 // timestep = min(timestep, (radii[k]+radii[n])/max_speed); 2177 //} else { 2178 // timestep = min(timestep, radii[k]/max_speed); 2179 // } 2180 } 2181 } 2182 2183 } // End edge i (and neighbour n) 2184 2185 2186 // Normalise triangle k by area and store for when all conserved 2187 // quantities get updated 2188 inv_area = 1.0/areas[k]; 2189 stage_explicit_update[k] *= inv_area; 2190 xmom_explicit_update[k] *= inv_area; 2191 ymom_explicit_update[k] *= inv_area; 2192 2193 2194 // Keep track of maximal speeds 2195 max_speed_array[k] = max_speed; 2196 2197 } // End triangle k 2198 2199 return timestep; 2003 double timestep, 2004 double epsilon, 2005 double H0, 2006 double g, 2007 long* neighbours, 2008 long* neighbour_edges, 2009 double* normals, 2010 double* edgelengths, 2011 double* radii, 2012 double* areas, 2013 long* tri_full_flag, 2014 double* stage_edge_values, 2015 double* xmom_edge_values, 2016 double* ymom_edge_values, 2017 double* bed_edge_values, 2018 double* stage_boundary_values, 2019 double* xmom_boundary_values, 2020 double* ymom_boundary_values, 2021 double* stage_explicit_update, 2022 double* xmom_explicit_update, 2023 double* ymom_explicit_update, 2024 long* already_computed_flux, 2025 double* max_speed_array, 2026 int optimise_dry_cells) { 2027 // Local variables 2028 double max_speed, length, inv_area, zl, zr; 2029 double h0 = H0*H0; // This ensures a good balance when h approaches H0. 2030 2031 double limiting_threshold = 10 * H0; // Avoid applying limiter below this 2032 // threshold for performance reasons. 2033 // See ANUGA manual under flux limiting 2034 int k, i, m, n; 2035 int ki, nm = 0, ki2; // Index shorthands 2036 2037 2038 // Workspace (making them static actually made function slightly slower (Ole)) 2039 double ql[3], qr[3], edgeflux[3]; // Work array for summing up fluxes 2040 2041 static long call = 1; // Static local variable flagging already computed flux 2042 2043 // Start computation 2044 call++; // Flag 'id' of flux calculation for this timestep 2045 2046 // Set explicit_update to zero for all conserved_quantities. 2047 // This assumes compute_fluxes called before forcing terms 2048 memset((char*) stage_explicit_update, 0, number_of_elements * sizeof (double)); 2049 memset((char*) xmom_explicit_update, 0, number_of_elements * sizeof (double)); 2050 memset((char*) ymom_explicit_update, 0, number_of_elements * sizeof (double)); 2051 2052 // For all triangles 2053 for (k = 0; k < number_of_elements; k++) { 2054 // Loop through neighbours and compute edge flux for each 2055 for (i = 0; i < 3; i++) { 2056 ki = k * 3 + i; // Linear index to edge i of triangle k 2057 2058 if (already_computed_flux[ki] == call) { 2059 // We've already computed the flux across this edge 2060 continue; 2061 } 2062 2063 // Get left hand side values from triangle k, edge i 2064 ql[0] = stage_edge_values[ki]; 2065 ql[1] = xmom_edge_values[ki]; 2066 ql[2] = ymom_edge_values[ki]; 2067 zl = bed_edge_values[ki]; 2068 2069 // Get right hand side values either from neighbouring triangle 2070 // or from boundary array (Quantities at neighbour on nearest face). 2071 n = neighbours[ki]; 2072 if (n < 0) { 2073 // Neighbour is a boundary condition 2074 m = n  1; // Convert negative flag to boundary index 2075 2076 qr[0] = stage_boundary_values[m]; 2077 qr[1] = xmom_boundary_values[m]; 2078 qr[2] = ymom_boundary_values[m]; 2079 zr = zl; // Extend bed elevation to boundary 2080 } 2081 else { 2082 // Neighbour is a real triangle 2083 m = neighbour_edges[ki]; 2084 nm = n * 3 + m; // Linear index (triangle n, edge m) 2085 2086 qr[0] = stage_edge_values[nm]; 2087 qr[1] = xmom_edge_values[nm]; 2088 qr[2] = ymom_edge_values[nm]; 2089 zr = bed_edge_values[nm]; 2090 } 2091 2092 // Now we have values for this edge  both from left and right side. 2093 2094 if (optimise_dry_cells) { 2095 // Check if flux calculation is necessary across this edge 2096 // This check will exclude dry cells. 2097 // This will also optimise cases where zl != zr as 2098 // long as both are dry 2099 2100 if (fabs(ql[0]  zl) < epsilon && 2101 fabs(qr[0]  zr) < epsilon) { 2102 // Cell boundary is dry 2103 2104 already_computed_flux[ki] = call; // #k Done 2105 if (n >= 0) { 2106 already_computed_flux[nm] = call; // #n Done 2107 } 2108 2109 max_speed = 0.0; 2110 continue; 2111 } 2112 } 2113 2114 2115 if (fabs(zlzr)>1.0e10) { 2116 report_python_error(AT,"Discontinuous Elevation"); 2117 return 0.0; 2118 } 2119 2120 // Outward pointing normal vector (domain.normals[k, 2*i:2*i+2]) 2121 ki2 = 2 * ki; //k*6 + i*2 2122 2123 // Edge flux computation (triangle k, edge i) 2124 _flux_function_central(ql, qr, zl, zr, 2125 normals[ki2], normals[ki2 + 1], 2126 epsilon, h0, limiting_threshold, g, 2127 edgeflux, &max_speed); 2128 2129 2130 // Multiply edgeflux by edgelength 2131 length = edgelengths[ki]; 2132 edgeflux[0] *= length; 2133 edgeflux[1] *= length; 2134 edgeflux[2] *= length; 2135 2136 2137 // Update triangle k with flux from edge i 2138 stage_explicit_update[k] = edgeflux[0]; 2139 xmom_explicit_update[k] = edgeflux[1]; 2140 ymom_explicit_update[k] = edgeflux[2]; 2141 2142 already_computed_flux[ki] = call; // #k Done 2143 2144 2145 // Update neighbour n with same flux but reversed sign 2146 if (n >= 0) { 2147 stage_explicit_update[n] += edgeflux[0]; 2148 xmom_explicit_update[n] += edgeflux[1]; 2149 ymom_explicit_update[n] += edgeflux[2]; 2150 2151 already_computed_flux[nm] = call; // #n Done 2152 } 2153 2154 // Update timestep based on edge i and possibly neighbour n 2155 if (tri_full_flag[k] == 1) { 2156 if (max_speed > epsilon) { 2157 // Apply CFL condition for triangles joining this edge (triangle k and triangle n) 2158 2159 // CFL for triangle k 2160 timestep = min(timestep, radii[k] / max_speed); 2161 2162 if (n >= 0) { 2163 // Apply CFL condition for neigbour n (which is on the ith edge of triangle k) 2164 timestep = min(timestep, radii[n] / max_speed); 2165 } 2166 2167 // Ted Rigby's suggested less conservative version 2168 //if (n>=0) { 2169 // timestep = min(timestep, (radii[k]+radii[n])/max_speed); 2170 //} else { 2171 // timestep = min(timestep, radii[k]/max_speed); 2172 // } 2173 } 2174 } 2175 2176 } // End edge i (and neighbour n) 2177 2178 2179 // Normalise triangle k by area and store for when all conserved 2180 // quantities get updated 2181 inv_area = 1.0 / areas[k]; 2182 stage_explicit_update[k] *= inv_area; 2183 xmom_explicit_update[k] *= inv_area; 2184 ymom_explicit_update[k] *= inv_area; 2185 2186 2187 // Keep track of maximal speeds 2188 max_speed_array[k] = max_speed; 2189 2190 } // End triangle k 2191 2192 return timestep; 2200 2193 } 2201 2194 … … 2264 2257 // Convert Python arguments to C 2265 2258 if (!PyArg_ParseTuple(args, "dOOOO", ×tep, &domain, &stage, &xmom, &ymom, &bed )) { 2266 PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");2267 return NULL;2259 report_python_error(AT, "could not parse input arguments"); 2260 return NULL; 2268 2261 } 2269 2262 … … 2450 2443 &max_speed_array, 2451 2444 &optimise_dry_cells)) { 2452 PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");2445 report_python_error(AT, "could not parse input arguments"); 2453 2446 return NULL; 2454 2447 } … … 2610 2603 &ymom_explicit_update, 2611 2604 &already_computed_flux)) { 2612 PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");2605 report_python_error(AT, "could not parse input arguments"); 2613 2606 return NULL; 2614 2607 } … … 2716 2709 &epsilon, 2717 2710 &wc, &zc, &xmomc, &ymomc)) { 2718 PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: protectcould not parse input arguments");2711 report_python_error(AT, "could not parse input arguments"); 2719 2712 return NULL; 2720 2713 } … … 2775 2768 &wv, &zv, &hvbar, 2776 2769 &xmomc, &ymomc, &xmomv, &ymomv)) { 2777 PyErr_SetString(PyExc_RuntimeError, 2778 "shallow_water_ext.c: balance_deep_and_shallow could not parse input arguments"); 2770 report_python_error(AT, "could not parse input arguments"); 2779 2771 return NULL; 2780 2772 } … … 2788 2780 Tmp = PyObject_GetAttrString(domain, "alpha_balance"); 2789 2781 if (!Tmp) { 2790 PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallowcould not obtain object alpha_balance from domain");2782 report_python_error(AT, "could not obtain object alpha_balance from domain"); 2791 2783 return NULL; 2792 2784 } … … 2797 2789 Tmp = PyObject_GetAttrString(domain, "H0"); 2798 2790 if (!Tmp) { 2799 PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallowcould not obtain object H0 from domain");2791 report_python_error(AT, "could not obtain object H0 from domain"); 2800 2792 return NULL; 2801 2793 } … … 2806 2798 Tmp = PyObject_GetAttrString(domain, "tight_slope_limiters"); 2807 2799 if (!Tmp) { 2808 PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallowcould not obtain object tight_slope_limiters from domain");2800 report_python_error(AT, "could not obtain object tight_slope_limiters from domain"); 2809 2801 return NULL; 2810 2802 } … … 2815 2807 Tmp = PyObject_GetAttrString(domain, "use_centroid_velocities"); 2816 2808 if (!Tmp) { 2817 PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallowcould not obtain object use_centroid_velocities from domain");2809 report_python_error(AT, "could not obtain object use_centroid_velocities from domain"); 2818 2810 return NULL; 2819 2811 } … … 2869 2861 &s_vec, &phi_vec, 2870 2862 &cw)) { 2871 PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: assign_windfield_valuescould not parse input arguments");2863 report_python_error(AT, "could not parse input arguments"); 2872 2864 return NULL; 2873 2865 } 
trunk/anuga_core/source/anuga/utilities/polygon_ext.c
r7276 r8017 20 20 #include "util_ext.h" 21 21 22 #define YES 1 23 #define NO 0 24 22 25 23 26 double dist(double x, … … 89 92 90 93 94 // public domain function by Darel Rex Finley, 2006 95 // http://www.alienryderflex.com/intersect/ 96 // 97 // Determines the intersection point of the line segment defined by points A and B 98 // with the line segment defined by points C and D. 99 // 100 // Returns YES if the intersection point was found, and stores that point in X,Y. 101 // Returns NO if there is no determinable intersection point, in which case X,Y will 102 // be unmodified. 103 104 int __lineSegmentIntersection( 105 double Ax, double Ay, 106 double Bx, double By, 107 double Cx, double Cy, 108 double Dx, double Dy, 109 double *X, double *Y) { 110 111 double distAB, theCos, theSin, newX, ABpos; 112 113 // Fail if either line segment is zerolength. 114 if ( (Ax == Bx && Ay == By)  (Cx == Dx && Cy == Dy) ) return NO ; 115 116 // Fail if the segments share an endpoint. 117 if ( (Ax == Cx && Ay == Cy)  (Bx == Cx && By == Cy) 118  (Ax == Dx && Ay == Dy)  (Bx == Dx && By == Dy) ) { 119 return NO; 120 } 121 122 // (1) Translate the system so that point A is on the origin. 123 Bx = Ax; 124 By = Ay; 125 Cx = Ax; 126 Cy = Ay; 127 Dx = Ax; 128 Dy = Ay; 129 130 // Discover the length of segment AB. 131 distAB = sqrt(Bx * Bx + By * By); 132 133 // (2) Rotate the system so that point B is on the positive X axis. 134 theCos = Bx / distAB; 135 theSin = By / distAB; 136 newX = Cx * theCos + Cy*theSin; 137 Cy = Cy * theCos  Cx*theSin; 138 Cx = newX; 139 newX = Dx * theCos + Dy*theSin; 140 Dy = Dy * theCos  Dx*theSin; 141 Dx = newX; 142 143 // Fail if segment CD doesn't cross line AB. 144 if ( (Cy < 0. && Dy < 0.)  (Cy >= 0. && Dy >= 0.) ) return NO; 145 146 // (3) Discover the position of the intersection point along line AB. 147 ABpos = Dx + (Cx  Dx) * Dy / (Dy  Cy); 148 149 // Fail if segment CD crosses line AB outside of segment AB. 150 if (ABpos < 0.  ABpos > distAB) return NO; 151 152 // (4) Apply the discovered position to line AB in the original coordinate system. 153 *X = Ax + ABpos*theCos; 154 *Y = Ay + ABpos*theSin; 155 156 // Success. 157 return YES; 158 } 91 159 92 160 /* 
trunk/anuga_core/source/anuga/utilities/util_ext.h
r7886 r8017 14 14 #include "math.h" 15 15 16 #include <stdio.h> 17 #define STRINGIFY(x) #x 18 #define TOSTRING(x) STRINGIFY(x) 19 #define AT __FILE__ ":" TOSTRING(__LINE__) 20 #define P_ERROR_BUFFER_SIZE 65 21 void report_python_error(const char *location, const char *msg) 22 { 23 24 char buf[P_ERROR_BUFFER_SIZE]; 25 int n; 26 n = snprintf(buf, P_ERROR_BUFFER_SIZE, "Error at %s: %s\n", location, msg); 27 28 PyErr_SetString(PyExc_RuntimeError, buf); 29 } 30 31 16 32 17 33 double max(double x, double y) {
Note: See TracChangeset
for help on using the changeset viewer.