Changeset 5321

Timestamp:

May 14, 2008, 2:25:10 PM (17 years ago)

Author:

ole

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:

• abstract_2d_finite_volumes/neighbour_mesh.py (modified) (3 diffs)
• fit_interpolate/general_fit_interpolate.py (modified) (1 diff)
• fit_interpolate/interpolate.py (modified) (5 diffs)
• shallow_water/data_manager.py (modified) (3 diffs)
• utilities/polygon.py (modified) (3 diffs)
• utilities/polygon_ext.c (modified) (3 diffs)

anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py

@@ -904 +904 @@
 
 
-    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
@@ -1047 +1049 @@
 
 
-    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:
@@ -1076 +1080 @@
       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]

anuga_core/source/anuga/fit_interpolate/general_fit_interpolate.py

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

@@ -123 +123 @@
         # 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( \
@@ -144 +144 @@
                 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 \
@@ -202 +201 @@
         f = ensure_numeric(f, Float)        
             
-
         self._A = self._build_interpolation_matrix_A(point_coordinates,
                                                      verbose=verbose)
@@ -262 +260 @@
         """
 
-        #print 'Building interpolation matrix'
+        if verbose: print 'Building interpolation matrix'
 
         # Convert point_coordinates to Numeric arrays, in case it was a list.
@@ -670 +668 @@
                     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

@@ -5655 +5655 @@
     
     # 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:
@@ -5668 +5668 @@
         interpolation_points.append(midpoint)
 
+
+    if verbose:
+        print 'Interpolating - ',        
+        print 'total number of interpolation points = %d'\
+              %len(interpolation_points)
 
     I = Interpolation_function(time,
@@ -5677 +5682 @@
                                verbose=verbose)
 
+    if verbose: print 'Computing hydrograph'
     # Compute hydrograph
     Q = []

anuga_core/source/anuga/utilities/polygon.py

@@ -106 +106 @@
             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\
@@ -173 +168 @@
             # 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
+
 
 
@@ -914 +950 @@
     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

@@ -58 +58 @@
 }
 
+
+
+/*
+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
+
+
+} 
+*/
 
 
@@ -176 +281 @@
 
 
+/*
+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) {
@@ -266 +398 @@
 
   {"_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"},