Changeset 6184

Timestamp:

Jan 17, 2009, 2:48:28 PM (16 years ago)

Author:

ole

Message:

Refactored computationally intensive part of interpolate_polyline out

File:

• anuga_core/source/anuga/fit_interpolate/interpolate.py (modified) (7 diffs)

anuga_core/source/anuga/fit_interpolate/interpolate.py

@@ -38 +38 @@
 from anuga.abstract_2d_finite_volumes.util import file_function
 from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
+from utilities.polygon import point_on_line
+
 
 import Numeric as num
@@ -246 +248 @@
           Interpolated values at inputted points (z).
         """
+        
+        #print f
+        #print vertex_coordinates
+        #print gauge_neighbour_id
+        #print point_coordinates
+
 
         # FIXME: There is an option of passing a tolerance into this
@@ -252 +260 @@
             point_coordinates = point_coordinates.get_data_points(absolute=True)
  
-        from utilities.polygon import point_on_line
-
         z = num.ones(len(point_coordinates), num.Float)
 
@@ -274 +280 @@
         
         elif n > 1:
-            for i in range(len(point_coordinates)):
-                found = False
-                for j in range(n):
-                    if gauge_neighbour_id[j] >= 0:
-                        if point_on_line(point_coordinates[i],
-                                         [vertex_coordinates[j],
-                                          vertex_coordinates[gauge_neighbour_id[j]]],
-                                         rtol=1.0e-6):
-                            found = True
-                            x0 = vertex_coordinates[j][0]
-                            y0 = vertex_coordinates[j][1]
-                            x1 = vertex_coordinates[gauge_neighbour_id[j]][0]
-                            y1 = vertex_coordinates[gauge_neighbour_id[j]][1]
-                            x2 = point_coordinates[i][0]
-                            y2 = point_coordinates[i][1]
-
-                            segment_len = sqrt((x1-x0)**2 + (y1-y0)**2)
-                            dist = sqrt((x2-x0)**2 + (y2-y0)**2)
-                            z[i] = (f[gauge_neighbour_id[j]] - f[j]) \
-                                   / segment_len*dist + f[j]
-                            break
-                                  
-                if not found:
-                    z[i] = 0.0
+            _interpolate_polyline_aux(n, z, f, point_coordinates, vertex_coordinates, gauge_neighbour_id)
+            
         return z
 
+    
 
     ##
@@ -361 +346 @@
                 point_coordinates = self._point_coordinates
             else:
-                #There are no good point_coordinates. import sys; sys.exit()
+                # There are no good point_coordinates. import sys; sys.exit()
                 msg = 'ERROR (interpolate.py): No point_coordinates inputted'
                 raise Exception(msg)
@@ -373 +358 @@
                                            verbose=verbose)
             else:
-                #Handle blocking
+                # Handle blocking
                 self._A_can_be_reused = False
                 start = 0
@@ -588 +573 @@
 
 
+        
+        
+        
+        
+def _interpolate_polyline_aux(n, z, f, point_coordinates, vertex_coordinates, gauge_neighbour_id):
+    """Auxiliary function
+    """
+    for i in range(len(point_coordinates)):
+    
+        x2 = point_coordinates[i][0]
+        y2 = point_coordinates[i][1]
+    
+        found = False
+        for j in range(n):
+            
+            neighbour_id = gauge_neighbour_id[j]
+            if neighbour_id >= 0:
+                x0 = vertex_coordinates[j][0]
+                y0 = vertex_coordinates[j][1]
+                x1 = vertex_coordinates[neighbour_id][0]
+                y1 = vertex_coordinates[neighbour_id][1]
+            
+                if point_on_line([x2,y2], 
+                                 [[x0, y0], [x1, y1]], 
+                                 rtol=1.0e-6):
+                                 
+                    found = True
+                    segment_len = sqrt((x1-x0)**2 + (y1-y0)**2)
+                    dist = sqrt((x2-x0)**2 + (y2-y0)**2)
+                    z[i] = (f[neighbour_id] - f[j])*dist/segment_len + f[j]
+                    break
+                                  
+        if not found:
+            z[i] = 0.0                
+        
+        
 ##
 # @brief ??