Changeset 8017


Ignore:
Timestamp:
Sep 16, 2010, 10:38:49 AM (8 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
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

    r7737 r8017  
    13091309        self.set_vertex_values(A, use_cache=use_cache, verbose=verbose)
    13101310
     1311
     1312
    13111313    ############################################################################
    13121314    # Methods for outputting model results
  • trunk/anuga_core/source/anuga/shallow_water/shallow_water_domain.py

    r7967 r8017  
    243243
    244244
    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.
    247247        0 Corresponds to first order, where as larger values make use of
    248248        the second order scheme.
  • trunk/anuga_core/source/anuga/shallow_water/shallow_water_ext.c

    r8009 r8017  
    2222// Shared code snippets
    2323#include "util_ext.h"
     24
     25
     26
    2427
    2528
     
    10451048            &normal, &ql, &qr, &zl, &zr, &edgeflux,
    10461049            &epsilon, &g, &H0)) {
    1047     PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: flux_function_central could not parse input arguments");
    1048     return NULL;
     1050      report_python_error(AT, "could not parse input arguments");
     1051      return NULL;
    10491052  }
    10501053
     
    10901093            &xmom, &ymom)) {
    10911094    //&epsilon)) {
    1092     PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: gravity could not parse input arguments");
    1093     return NULL;
     1095      report_python_error(AT, "could not parse input arguments");
     1096      return NULL;
    10941097  }
    10951098
     
    11541157  if (!PyArg_ParseTuple(args, "ddOOOOOOOO",
    11551158                        &g, &eps, &x, &w, &uh, &vh, &z,  &eta, &xmom, &ymom)) {
    1156     PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: manning_friction_new could not parse input arguments");
    1157     return NULL;
     1159      report_python_error(AT, "could not parse input arguments");
     1160      return NULL;
    11581161  }
    11591162
     
    11961199  if (!PyArg_ParseTuple(args, "ddOOOOOOOO",
    11971200                        &g, &eps, &x, &w, &uh, &vh, &z,  &chezy, &xmom, &ymom)) {
    1198     PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: chezy_friction could not parse input arguments");
    1199     return NULL;
     1201      report_python_error(AT, "could not parse input arguments");
     1202      return NULL;
    12001203  }
    12011204
     
    12381241  if (!PyArg_ParseTuple(args, "ddOOOOOOO",
    12391242                        &g, &eps, &w, &uh, &vh, &z,  &eta, &xmom, &ymom)) {
    1240     PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: manning_friction_old could not parse input arguments");
    1241     return NULL;
     1243      report_python_error(AT, "could not parse input arguments");
     1244      return NULL;
    12421245  }
    12431246
     
    14381441      if (area2 <= 0)
    14391442      {
    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 anti-clockwise.");
    1442    
    1443         return -1;
     1443          report_python_error(AT, "Negative triangle area");
     1444          return -1;
    14441445      } 
    14451446     
     
    16301631      {
    16311632        // 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");
    16341634        return -1;
    16351635      }
     
    18661866            &elevation_vertex_values,
    18671867            &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");
    18711870    return NULL;
    18721871  }
     
    19571956  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OO|i", argnames,
    19581957                   &Q, &Normal, &direction)) {
    1959     PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: rotate could not parse input arguments");
     1958    report_python_error(AT, "could not parse input arguments");
    19601959    return NULL;
    19611960  } 
     
    19691968
    19701969  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");
    19721971    return NULL;
    19731972  }
     
    20001999
    20012000// Computational function for flux computation
     2001
    20022002double _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(zl-zr)>1.0e-10) {
     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;
    22002193}
    22012194
     
    22642257    // Convert Python arguments to C
    22652258    if (!PyArg_ParseTuple(args, "dOOOO", &timestep, &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;
    22682261    }
    22692262
     
    24502443            &max_speed_array,
    24512444            &optimise_dry_cells)) {
    2452     PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
     2445    report_python_error(AT, "could not parse input arguments");
    24532446    return NULL;
    24542447  }
     
    26102603            &ymom_explicit_update,
    26112604            &already_computed_flux)) {
    2612     PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
     2605    report_python_error(AT, "could not parse input arguments");
    26132606    return NULL;
    26142607  }
     
    27162709            &epsilon,
    27172710            &wc, &zc, &xmomc, &ymomc)) {
    2718     PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: protect could not parse input arguments");
     2711    report_python_error(AT, "could not parse input arguments");
    27192712    return NULL;
    27202713  } 
     
    27752768            &wv, &zv, &hvbar,
    27762769            &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");
    27792771    return NULL;
    27802772  } 
     
    27882780  Tmp = PyObject_GetAttrString(domain, "alpha_balance");
    27892781  if (!Tmp) {
    2790     PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not obtain object alpha_balance from domain");
     2782    report_python_error(AT, "could not obtain object alpha_balance from domain");
    27912783    return NULL;
    27922784  } 
     
    27972789  Tmp = PyObject_GetAttrString(domain, "H0");
    27982790  if (!Tmp) {
    2799     PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not obtain object H0 from domain");
     2791    report_python_error(AT, "could not obtain object H0 from domain");
    28002792    return NULL;
    28012793  } 
     
    28062798  Tmp = PyObject_GetAttrString(domain, "tight_slope_limiters");
    28072799  if (!Tmp) {
    2808     PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not obtain object tight_slope_limiters from domain");
     2800    report_python_error(AT, "could not obtain object tight_slope_limiters from domain");
    28092801    return NULL;
    28102802  } 
     
    28152807  Tmp = PyObject_GetAttrString(domain, "use_centroid_velocities");
    28162808  if (!Tmp) {
    2817     PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not obtain object use_centroid_velocities from domain");
     2809    report_python_error(AT, "could not obtain object use_centroid_velocities from domain");
    28182810    return NULL;
    28192811  } 
     
    28692861            &s_vec, &phi_vec,
    28702862            &cw)) {
    2871     PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: assign_windfield_values could not parse input arguments");
     2863    report_python_error(AT, "could not parse input arguments");
    28722864    return NULL;
    28732865  } 
  • trunk/anuga_core/source/anuga/utilities/polygon_ext.c

    r7276 r8017  
    2020#include "util_ext.h"
    2121
     22#define YES 1
     23#define NO 0
     24
    2225
    2326double dist(double x,
     
    8992
    9093
     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
     104int __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 zero-length.
     114    if ( (Ax == Bx && Ay == By) || (Cx == Dx && Cy == Dy) ) return NO ;
     115
     116    //  Fail if the segments share an end-point.
     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 A-B.
     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 C-D doesn't cross line A-B.
     144    if ( (Cy < 0. && Dy < 0.) || (Cy >= 0. && Dy >= 0.) ) return NO;
     145
     146    //  (3) Discover the position of the intersection point along line A-B.
     147    ABpos = Dx + (Cx - Dx) * Dy / (Dy - Cy);
     148
     149    //  Fail if segment C-D crosses line A-B outside of segment A-B.
     150    if (ABpos < 0. || ABpos > distAB) return NO;
     151
     152    //  (4) Apply the discovered position to line A-B in the original coordinate system.
     153    *X = Ax + ABpos*theCos;
     154    *Y = Ay + ABpos*theSin;
     155
     156    //  Success.
     157    return YES;
     158}
    91159
    92160/*
  • trunk/anuga_core/source/anuga/utilities/util_ext.h

    r7886 r8017  
    1414#include "math.h"
    1515
     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
     21void 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
    1632
    1733double max(double x, double y) { 
Note: See TracChangeset for help on using the changeset viewer.