Changeset 5862

Timestamp:

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

Author:

ole

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:

• abstract_2d_finite_volumes/neighbour_mesh.py (modified) (3 diffs)
• abstract_2d_finite_volumes/quantity.py (modified) (3 diffs)
• culvert_flows/Test_Culvert_Flat_Water_Lev.py (modified) (3 diffs)
• culvert_flows/culvert_routines.py (modified) (1 diff)
• fit_interpolate/interpolate.py (modified) (1 diff)
• shallow_water/data_manager.py (modified) (1 diff)
• shallow_water/shallow_water_domain.py (modified) (2 diffs)
• shallow_water/test_shallow_water_domain.py (modified) (1 diff)

anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py

@@ -9 +9 @@
 
 from general_mesh import General_mesh
+from anuga.caching import cache
 from math import pi, sqrt
 from Numeric import array, allclose
@@ -904 +905 @@
 
 
-    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[0][0]
-      eta0 = line[0][1]
-
-      
-      # 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[0]) or\
-                  #       allclose(value, edge[1])):
-                  #    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[1], poly, closed=False):
-                  intersections[tuple(line[1])] = i
-              elif is_inside_polygon(line[0], poly, closed=False):
-                  intersections[tuple(line[0])] = 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[0]
-              x1, y1 = points[1]
-
-
-              # 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[1], -vector[0]])/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[0], point0[1],
-                                                    point1[0], point1[1])
-             
-          
-          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
+        
   
 
@@ -1149 +1004 @@
         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[0][0]
+    eta0 = line[0][1]
+
+  
+    # 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[0]) or\
+                #       allclose(value, edge[1])):
+                #    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[1], poly, closed=False):
+                intersections[tuple(line[1])] = i
+            elif is_inside_polygon(line[0], poly, closed=False):
+                intersections[tuple(line[0])] = 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[0]
+            x1, y1 = points[1]
+
+
+            # 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[1], -vector[0]])/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[0], point0[1],
+                                                  point1[0], point1[1])
+           
+        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

@@ -873 +873 @@
             # 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,
@@ -1054 +1054 @@
         # Speed is also a consideration here.
         
+        
         if isinstance(interpolation_points, Geospatial_data):        
             # Ensure interpolation points are in absolute UTM coordinates
@@ -1077 +1078 @@
                            use_cache=True, 
                            verbose=verbose)
+                           
 
 

anuga_core/source/anuga/culvert_flows/Test_Culvert_Flat_Water_Lev.py

@@ -30 +30 @@
 #------------------------------------------------------------------------------
 print 'Setting up domain'
-
-# Open file for storing some specific results...
-fid = open('Culvert_Headwall', 'w')
 
 length = 40.
@@ -149 +146 @@
 #------------------------------------------------------------------------------
 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
@@ -157 +153 @@
 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

@@ -61 +61 @@
     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

@@ -463 +463 @@
         
         # 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

@@ -6377 +6377 @@
     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

@@ -399 +399 @@
         """       
         
-        # 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)       
@@ -472 +468 @@
                                          
         
-        # 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

@@ -6009 +6009 @@
     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')