# Changeset 5321

Ignore:
Timestamp:
May 14, 2008, 2:25:10 PM (15 years ago)
Message:

Investigated performance profile of get_flow_through_cross_sections.
Most time appear to be spent in the interpolation routines.

Location:
anuga_core/source/anuga
Files:
6 edited

Unmodified
Removed
• ## anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py

 r5288 def _get_intersecting_segments(self, line): def _get_intersecting_segments(self, line, verbose=False): """Find edges intersected by line Input: line - list of two points forming a segmented line line - list of two points forming a segmented line verbose Output: list of instances of class Triangle_intersection def get_intersecting_segments(self, polyline): def get_intersecting_segments(self, polyline, verbose=False): """Find edges intersected by polyline Input: polyline - list of points forming a segmented line verbose Output: triangle_intersections = [] for i, point0 in enumerate(polyline[:-1]): point1 = polyline[i+1] if verbose: print 'Extracting mesh intersections from line:', print '(%.2f, %.2f) - (%.2f, %.2f)' %(point0, point0, point1, point1) line = [point0, point1]
• ## anuga_core/source/anuga/fit_interpolate/general_fit_interpolate.py

 r4859 if verbose: print 'FitInterpolate: Building mesh' self.mesh = Mesh(vertex_coordinates, triangles) self.mesh.check_integrity() #self.mesh.check_integrity() # Time consuming else: self.mesh = mesh
• ## anuga_core/source/anuga/fit_interpolate/interpolate.py

 r5288 # method is called even if interpolation points are unchanged. #print "point_coordinates interpolate.interpolate", point_coordinates #print "point_coordinates interpolate.interpolate", point_coordinates if verbose: print 'Build intepolation object' if isinstance(point_coordinates, Geospatial_data): point_coordinates = point_coordinates.get_data_points( \ msg = 'ERROR (interpolate.py): No point_coordinates inputted' raise Exception(msg) if point_coordinates is not None: if point_coordinates is not None: self._point_coordinates = point_coordinates if len(point_coordinates) < start_blocking_len or \ f = ensure_numeric(f, Float) self._A = self._build_interpolation_matrix_A(point_coordinates, verbose=verbose) """ #print 'Building interpolation matrix' if verbose: print 'Building interpolation matrix' # Convert point_coordinates to Numeric arrays, in case it was a list. result = interpol.interpolate(Q, point_coordinates=\ self.interpolation_points) self.interpolation_points, verbose=False) # Don't clutter #assert len(result), len(interpolation_points)
• ## anuga_core/source/anuga/shallow_water/data_manager.py

 r5307 # Find all intersections and associated triangles. segments = mesh.get_intersecting_segments(polyline) segments = mesh.get_intersecting_segments(polyline, verbose=verbose) #print #for x in segments: interpolation_points.append(midpoint) if verbose: print 'Interpolating - ', print 'total number of interpolation points = %d'\ %len(interpolation_points) I = Interpolation_function(time, verbose=verbose) if verbose: print 'Computing hydrograph' # Compute hydrograph Q = []
• ## anuga_core/source/anuga/utilities/polygon.py

 r5286 if point_on_line([x3, y3], line0): line1_ends_on_line0 = True #print line0_starts_on_line1 #print line1_starts_on_line0 #print line1_ends_on_line0 #print line0_ends_on_line1 if not(line0_starts_on_line1 or line0_ends_on_line1\ # No intersection return 0, None def NEW_C_intersection(line0, line1): #FIXME(Ole): To write in C """Returns intersecting point between two line segments or None (if parallel or no intersection is found). However, if parallel lines coincide partly (i.e. shara a common segment, the line segment where lines coincide is returned Inputs: line0, line1: Each defined by two end points as in: [[x0, y0], [x1, y1]] A line can also be a 2x2 numeric array with each row corresponding to a point. Output: status, value where status is interpreted as follows status == 0: no intersection with value set to None status == 1: One intersection point found and returned in value as [x,y] status == 2: Coinciding line segment found. Value taks the form [[x0,y0], [x1,y1]] status == 3: Lines would coincide but only if extended. Value set to None status == 4: Lines are parallel with a fixed distance apart. Value set to None. """ line0 = ensure_numeric(line0, Float) line1 = ensure_numeric(line1, Float) status, value = _intersection(line0[0,0], line0[0,1], line0[1,0], line0[1,1], line1[0,0], line1[0,1], line1[1,0], line1[1,1]) return status, value from polygon_ext import _point_on_line from polygon_ext import _separate_points_by_polygon #from polygon_ext import _intersection else:
• ## anuga_core/source/anuga/utilities/polygon_ext.c

 r5225 } /* WORK IN PROGRESS TO OPTIMISE INTERSECTION int __intersection(double x0, double y0, double x1, double y1) { x0 = line0[0,0]; y0 = line0[0,1] x1 = line0[1,0]; y1 = line0[1,1] x2 = line1[0,0]; y2 = line1[0,1] x3 = line1[1,0]; y3 = line1[1,1] denom = (y3-y2)*(x1-x0) - (x3-x2)*(y1-y0) u0 = (x3-x2)*(y0-y2) - (y3-y2)*(x0-x2) u1 = (x2-x0)*(y1-y0) - (y2-y0)*(x1-x0) if allclose(denom, 0.0): # Lines are parallel - check if they coincide on a shared a segment if allclose( [u0, u1], 0.0 ): # We now know that the lines if continued coincide # The remaining check will establish if the finite lines share a segment line0_starts_on_line1 = line0_ends_on_line1 =\ line1_starts_on_line0 = line1_ends_on_line0 = False if point_on_line([x0, y0], line1): line0_starts_on_line1 = True if point_on_line([x1, y1], line1): line0_ends_on_line1 = True if point_on_line([x2, y2], line0): line1_starts_on_line0 = True if point_on_line([x3, y3], line0): line1_ends_on_line0 = True if not(line0_starts_on_line1 or line0_ends_on_line1\ or line1_starts_on_line0 or line1_ends_on_line0): # Lines are parallel and would coincide if extended, but not as they are. return 3, None # One line fully included in the other. Use direction of included line if line0_starts_on_line1 and line0_ends_on_line1: # Shared segment is line0 fully included in line1 segment = array([[x0, y0], [x1, y1]]) if line1_starts_on_line0 and line1_ends_on_line0: # Shared segment is line1 fully included in line0 segment = array([[x2, y2], [x3, y3]]) # Overlap with lines are oriented the same way if line0_starts_on_line1 and line1_ends_on_line0: # Shared segment from line0 start to line 1 end segment = array([[x0, y0], [x3, y3]]) if line1_starts_on_line0 and line0_ends_on_line1: # Shared segment from line1 start to line 0 end segment = array([[x2, y2], [x1, y1]]) # Overlap in opposite directions - use direction of line0 if line0_starts_on_line1 and line1_starts_on_line0: # Shared segment from line0 start to line 1 end segment = array([[x0, y0], [x2, y2]]) if line0_ends_on_line1 and line1_ends_on_line0: # Shared segment from line0 start to line 1 end segment = array([[x3, y3], [x1, y1]]) return 2, segment else: # Lines are parallel but they do not coincide return 4, None #FIXME (Ole): Add distance here instead of None else: # Lines are not parallel or coinciding u0 = u0/denom u1 = u1/denom x = x0 + u0*(x1-x0) y = y0 + u0*(y1-y0) # Sanity check - can be removed to speed up if needed assert allclose(x, x2 + u1*(x3-x2)) assert allclose(y, y2 + u1*(y3-y2)) # Check if point found lies within given line segments if 0.0 <= u0 <= 1.0 and 0.0 <= u1 <= 1.0: # We have intersection return 1, array([x, y]) else: # No intersection return 0, None } */ /* PyObject *_intersection(PyObject *self, PyObject *args) { // // intersection(x0, y0, x1, y1) // double x, y, x0, y0, x1, y1; int res; PyObject *result; // Convert Python arguments to C if (!PyArg_ParseTuple(args, "dddddd", &x, &y, &x0, &y0, &x1, &y1)) { PyErr_SetString(PyExc_RuntimeError, "point_on_line could not parse input"); return NULL; } // Call underlying routine res = __intersection(x, y, x0, y0, x1, y1); // Return values a and b result = Py_BuildValue("i", res); return result; } */ PyObject *_separate_points_by_polygon(PyObject *self, PyObject *args) { {"_point_on_line", _point_on_line, METH_VARARGS, "Print out"}, //{"_intersection", _intersection, METH_VARARGS, "Print out"}, {"_separate_points_by_polygon", _separate_points_by_polygon, METH_VARARGS, "Print out"},
Note: See TracChangeset for help on using the changeset viewer.