# Changeset 5862

Ignore:
Timestamp:
Oct 24, 2008, 4:18:35 PM (15 years ago)
Message:

Made get_intersecting_segments cacheable and used that
with get_energy_through_cross_section and similar functions.

This helped the culvert routine along heaps, but I still believe
the interpolation could be made much more efficient.

Location:
anuga_core/source/anuga
Files:
8 edited

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

 r5729 from general_mesh import General_mesh from anuga.caching import cache from math import pi, sqrt from Numeric import array, allclose def _get_intersecting_segments(self, line, verbose=False): """Find edges intersected by line Input: line - list of two points forming a segmented line verbose Output: list of instances of class Triangle_intersection This method is used by the public method get_intersecting_segments(self, polyline) which also contains more documentation. """ from anuga.utilities.polygon import intersection from anuga.utilities.polygon import is_inside_polygon msg = 'Line segment must contain exactly two points' assert len(line) == 2, msg # Origin of intersecting line to be used for # establishing direction xi0 = line eta0 = line # Check intersection with edge segments for all triangles # FIXME (Ole): This should be implemented in C V = self.get_vertex_coordinates() N = len(self) triangle_intersections={} # Keep track of segments already done for i in range(N): # Get nodes and edge segments for each triangle x0, y0 = V[3*i, :] x1, y1 = V[3*i+1, :] x2, y2 = V[3*i+2, :] edge_segments = [[[x0,y0], [x1, y1]], [[x1,y1], [x2, y2]], [[x2,y2], [x0, y0]]] # Find segments that are intersected by line intersections = {} # Use dictionary to record points only once for edge in edge_segments: status, value = intersection(line, edge) #if value is not None: print 'Triangle %d, Got' %i, status, value if status == 1: # Normal intersection of one edge or vertex intersections[tuple(value)] = i # Exclude singular intersections with vertices #if not(allclose(value, edge) or\ #       allclose(value, edge)): #    intersections.append(value) if status == 2: # Edge is sharing a segment with line # This is usually covered by the two # vertices that would have been picked up # under status == 1. # However, if coinciding line stops partway # along this edge, it will be recorded here. intersections[tuple(value[0,:])] = i intersections[tuple(value[1,:])] = i if len(intersections) == 1: # Check if either line end point lies fully within this triangle # If this is the case accept that as one end of the intersecting # segment poly = V[3*i:3*i+3] if is_inside_polygon(line, poly, closed=False): intersections[tuple(line)] = i elif is_inside_polygon(line, poly, closed=False): intersections[tuple(line)] = i else: # Ignore situations where one vertex is touch, for instance continue msg = 'There can be only two or no intersections' assert len(intersections) in [0,2], msg if len(intersections) == 2: # Calculate attributes for this segment # End points of intersecting segment points = intersections.keys() x0, y0 = points x1, y1 = points # Determine which end point is closer to the origin of the line # This is necessary for determining the direction of # the line and the normals # Distances from line origin to the two intersections z0 = array([x0 - xi0, y0 - eta0]) z1 = array([x1 - xi0, y1 - eta0]) d0 = sqrt(sum(z0**2)) d1 = sqrt(sum(z1**2)) if d1 < d0: # Swap xi, eta = x0, y0 x0, y0 = x1, y1 x1, y1 = xi, eta # (x0,y0) is now the origin of the intersecting segment # Normal direction: # Right hand side relative to line direction vector = array([x1 - x0, y1 - y0]) # Segment vector length = sqrt(sum(vector**2))      # Segment length normal = array([vector, -vector])/length segment = ((x0,y0), (x1, y1)) T = Triangle_intersection(segment=segment, normal=normal, length=length, triangle_id=i) # Add segment unless it was done earlier if not triangle_intersections.has_key(segment): triangle_intersections[segment] = T # Return segments as a list return triangle_intersections.values() def get_intersecting_segments(self, polyline, use_cache=False, verbose=False): """Find edges intersected by polyline Input: polyline - list of points forming a segmented line verbose Output: list of instances of class Triangle_intersection The polyline may break inside any triangle causing multiple segments per triangle - consequently the same triangle may appear in several entries. If a polyline segment coincides with a triangle edge, the the entire shared segment will be used. Onle one of the triangles thus intersected will be used and that is the first one encoutered. Intersections with single vertices are ignored. Resulting segments are unsorted """ msg = 'Polyline must contain at least two points' assert len(polyline) >= 2, msg # For all segments in polyline 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] triangle_intersections += self._get_intersecting_segments(line) return triangle_intersections """Find edges intersected by polyline Input: polyline - list of points forming a segmented line use_cache verbose Output: list of instances of class Triangle_intersection The polyline may break inside any triangle causing multiple segments per triangle - consequently the same triangle may appear in several entries. If a polyline segment coincides with a triangle edge, the the entire shared segment will be used. Onle one of the triangles thus intersected will be used and that is the first one encountered. Intersections with single vertices are ignored. Resulting segments are unsorted """ V = self.get_vertex_coordinates() N = len(self) # Adjust polyline to mesh spatial origin polyline = self.geo_reference.get_relative(polyline) if use_cache is True: segments = cache(get_intersecting_segments, (V, N, polyline), {'verbose': verbose}, verbose=verbose) else: segments = get_intersecting_segments(V, N, polyline, verbose=verbose) return segments return s def _get_intersecting_segments(V, N, line, verbose=False): """Find edges intersected by line Input: V: Vertex coordinates as obtained by mesh.get_vertex_coordinates() N: Number of triangles in mesh line - list of two points forming a segmented line verbose Output: list of instances of class Triangle_intersection This method is used by the public method get_intersecting_segments(self, polyline) which also contains more documentation. """ from anuga.utilities.polygon import intersection from anuga.utilities.polygon import is_inside_polygon msg = 'Line segment must contain exactly two points' assert len(line) == 2, msg # Origin of intersecting line to be used for # establishing direction xi0 = line eta0 = line # Check intersection with edge segments for all triangles # FIXME (Ole): This should be implemented in C triangle_intersections={} # Keep track of segments already done for i in range(N): # Get nodes and edge segments for each triangle x0, y0 = V[3*i, :] x1, y1 = V[3*i+1, :] x2, y2 = V[3*i+2, :] edge_segments = [[[x0,y0], [x1, y1]], [[x1,y1], [x2, y2]], [[x2,y2], [x0, y0]]] # Find segments that are intersected by line intersections = {} # Use dictionary to record points only once for edge in edge_segments: status, value = intersection(line, edge) #if value is not None: print 'Triangle %d, Got' %i, status, value if status == 1: # Normal intersection of one edge or vertex intersections[tuple(value)] = i # Exclude singular intersections with vertices #if not(allclose(value, edge) or\ #       allclose(value, edge)): #    intersections.append(value) if status == 2: # Edge is sharing a segment with line # This is usually covered by the two # vertices that would have been picked up # under status == 1. # However, if coinciding line stops partway # along this edge, it will be recorded here. intersections[tuple(value[0,:])] = i intersections[tuple(value[1,:])] = i if len(intersections) == 1: # Check if either line end point lies fully within this triangle # If this is the case accept that as one end of the intersecting # segment poly = V[3*i:3*i+3] if is_inside_polygon(line, poly, closed=False): intersections[tuple(line)] = i elif is_inside_polygon(line, poly, closed=False): intersections[tuple(line)] = i else: # Ignore situations where one vertex is touch, for instance continue msg = 'There can be only two or no intersections' assert len(intersections) in [0,2], msg if len(intersections) == 2: # Calculate attributes for this segment # End points of intersecting segment points = intersections.keys() x0, y0 = points x1, y1 = points # Determine which end point is closer to the origin of the line # This is necessary for determining the direction of # the line and the normals # Distances from line origin to the two intersections z0 = array([x0 - xi0, y0 - eta0]) z1 = array([x1 - xi0, y1 - eta0]) d0 = sqrt(sum(z0**2)) d1 = sqrt(sum(z1**2)) if d1 < d0: # Swap xi, eta = x0, y0 x0, y0 = x1, y1 x1, y1 = xi, eta # (x0,y0) is now the origin of the intersecting segment # Normal direction: # Right hand side relative to line direction vector = array([x1 - x0, y1 - y0]) # Segment vector length = sqrt(sum(vector**2))      # Segment length normal = array([vector, -vector])/length segment = ((x0,y0), (x1, y1)) T = Triangle_intersection(segment=segment, normal=normal, length=length, triangle_id=i) # Add segment unless it was done earlier if not triangle_intersections.has_key(segment): triangle_intersections[segment] = T # Return segments as a list return triangle_intersections.values() def get_intersecting_segments(V, N, polyline, verbose=False): """Internal function to find edges intersected by Polyline Input: V: Vertex coordinates as obtained by mesh.get_vertex_coordinates() N: Number of triangles in mesh polyline - list of points forming a segmented line verbose Output: list of instances of class Triangle_intersection This method is used by the public method get_intersecting_segments(self, polyline) which also contains more documentation. """ msg = 'Polyline must contain at least two points' assert len(polyline) >= 2, msg # For all segments in polyline 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] triangle_intersections += _get_intersecting_segments(V, N, line, verbose=verbose) msg = 'No segments found' assert len(triangle_intersections) > 0, msg return triangle_intersections def segment_midpoints(segments): """Calculate midpoints of all segments
• ## anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

 r5855 # This variant will cause Mesh object to be recreated # in fit_to_mesh thus doubling up on the neighbour structure coordinates = self.domain.get_nodes(absolute=True) triangles = self.domain.triangles      #FIXME nodes = self.domain.get_nodes(absolute=True) triangles = self.domain.triangles vertex_attributes = fit_to_mesh(filename, coordinates, triangles, mesh=None, nodes, triangles, mesh=None, alpha=alpha, attribute_name=attribute_name, # Speed is also a consideration here. if isinstance(interpolation_points, Geospatial_data): # Ensure interpolation points are in absolute UTM coordinates use_cache=True, verbose=verbose)
• ## anuga_core/source/anuga/culvert_flows/Test_Culvert_Flat_Water_Lev.py

 r5777 #------------------------------------------------------------------------------ print 'Setting up domain' # Open file for storing some specific results... fid = open('Culvert_Headwall', 'w') length = 40. #------------------------------------------------------------------------------ print 'Setting Boundary Conditions' #Bi = Dirichlet_boundary([0.5, 0.0, 0.0])          # Inflow based on Flow Depth (0.5m) and Approaching Momentum !!! Bi = Dirichlet_boundary([0.0, 0.0, 0.0])          # Inflow based on Flow Depth and Approaching Momentum !!! Br = Reflective_boundary(domain)              # Solid reflective wall domain.set_boundary({'left': Btus, 'right': Btds, 'top': Br, 'bottom': Br}) #------------------------------------------------------------------------------ # Setup Application of  specialised forcing terms #------------------------------------------------------------------------------ print 'Setting up Forcing Terms if Required' # This is the new element implemented by Ole to allow direct input of Inflow in m^3/s #fixed_flow = Inflow(domain, center=(2.1, 2.1), radius=1.261566, flow=1.00)   #   Fixed Flow Value Over Area of 5m2 at 1m/s = 5m^3/s #            flow=file_function('Q/QPMF_Rot_Sub13.tms'))        # Read Time Series in  from File #             flow=lambda t: min(0.01*t, 0.01942))                             # Time Varying Function   Tap turning up #domain.forcing_terms.append(fixed_flow) #------------------------------------------------------------------------------ # Evolve system through time #------------------------------------------------------------------------------ print 'STARTING EVOLVE ======>' #for t in domain.evolve(yieldstep = 0.01, finaltime = 45): #     if int(domain.time*100) % 100 == 0: #            domain.write_time() for t in domain.evolve(yieldstep = 0.01, finaltime = 45): if int(domain.time*100) % 100 == 0: domain.write_time() # Profiling code import time t0 = time.time() #if domain.get_time() >= 4 and tap.flow != 0.0: #    print 'Turning tap off' #    tap.flow = 0.0 #if domain.get_time() >= 3 and sink.flow < 0.0: #    print 'Turning drain on' #    sink.flow = -0.8 # Close s = 'for t in domain.evolve(yieldstep = 0.5, finaltime = 2): domain.write_time()' fid.close() import profile, pstats FN = 'profile.dat' profile.run(s, FN) print 'That took %.2f seconds' %(time.time()-t0) S = pstats.Stats(FN) #S.sort_stats('time').print_stats(20) s = S.sort_stats('cumulative').print_stats(30) print s
• ## anuga_core/source/anuga/culvert_flows/culvert_routines.py

 r5437 if inlet.depth_trigger >= 0.01 and inlet.depth >= 0.01: # Calculate driving energy E = inlet.specific_energy E = inlet.total_energy s = 'Driving energy  = %f m' %E
• ## anuga_core/source/anuga/fit_interpolate/interpolate.py

 r5855 # Build n x m interpolation matrix if verbose and len(self.outside_poly_indices) > 0: if verbose and len(outside_poly_indices) > 0: print '\n WARNING: Points outside mesh boundary. \n'
• ## anuga_core/source/anuga/shallow_water/data_manager.py

 r5773 mesh, quantities, time = X # Adjust polyline to mesh spatial origin polyline = mesh.geo_reference.get_relative(polyline) # Find all intersections and associated triangles. segments = mesh.get_intersecting_segments(polyline, verbose=verbose)
• ## anuga_core/source/anuga/shallow_water/shallow_water_domain.py

 r5773 """ # Adjust polyline to mesh spatial origin polyline = self.geo_reference.get_relative(polyline) # Find all intersections and associated triangles. segments = self.get_intersecting_segments(polyline, verbose=verbose) msg = 'No segments found' assert len(segments) > 0, msg segments = self.get_intersecting_segments(polyline, use_cache=True, verbose=verbose) # Get midpoints midpoints = segment_midpoints(segments) # Adjust polyline to mesh spatial origin polyline = self.geo_reference.get_relative(polyline) # Find all intersections and associated triangles. segments = self.get_intersecting_segments(polyline, verbose=verbose) msg = 'No segments found' assert len(segments) > 0, msg segments = self.get_intersecting_segments(polyline, use_cache=True, verbose=verbose) # Get midpoints midpoints = segment_midpoints(segments)
• ## anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

 r5736 suite = unittest.makeSuite(Test_Shallow_Water,'test') #suite = unittest.makeSuite(Test_Shallow_Water,'test_get_flow_through_cross_section_with_geo') #suite = unittest.makeSuite(Test_Shallow_Water,'test_get_energy_through_cross_section_with_g') #suite = unittest.makeSuite(Test_Shallow_Water,'test_fitting_using_shallow_water_domain') #suite = unittest.makeSuite(Test_Shallow_Water,'test_tight_slope_limiters')
Note: See TracChangeset for help on using the changeset viewer.