Changeset 6184


Ignore:
Timestamp:
Jan 17, 2009, 2:48:28 PM (10 years ago)
Author:
ole
Message:

Refactored computationally intensive part of interpolate_polyline out

File:
1 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/fit_interpolate/interpolate.py

    r6173 r6184  
    3838from anuga.abstract_2d_finite_volumes.util import file_function
    3939from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
     40from utilities.polygon import point_on_line
     41
    4042
    4143import Numeric as num
     
    246248          Interpolated values at inputted points (z).
    247249        """
     250       
     251        #print f
     252        #print vertex_coordinates
     253        #print gauge_neighbour_id
     254        #print point_coordinates
     255
    248256
    249257        # FIXME: There is an option of passing a tolerance into this
     
    252260            point_coordinates = point_coordinates.get_data_points(absolute=True)
    253261 
    254         from utilities.polygon import point_on_line
    255 
    256262        z = num.ones(len(point_coordinates), num.Float)
    257263
     
    274280       
    275281        elif n > 1:
    276             for i in range(len(point_coordinates)):
    277                 found = False
    278                 for j in range(n):
    279                     if gauge_neighbour_id[j] >= 0:
    280                         if point_on_line(point_coordinates[i],
    281                                          [vertex_coordinates[j],
    282                                           vertex_coordinates[gauge_neighbour_id[j]]],
    283                                          rtol=1.0e-6):
    284                             found = True
    285                             x0 = vertex_coordinates[j][0]
    286                             y0 = vertex_coordinates[j][1]
    287                             x1 = vertex_coordinates[gauge_neighbour_id[j]][0]
    288                             y1 = vertex_coordinates[gauge_neighbour_id[j]][1]
    289                             x2 = point_coordinates[i][0]
    290                             y2 = point_coordinates[i][1]
    291 
    292                             segment_len = sqrt((x1-x0)**2 + (y1-y0)**2)
    293                             dist = sqrt((x2-x0)**2 + (y2-y0)**2)
    294                             z[i] = (f[gauge_neighbour_id[j]] - f[j]) \
    295                                    / segment_len*dist + f[j]
    296                             break
    297                                  
    298                 if not found:
    299                     z[i] = 0.0
     282            _interpolate_polyline_aux(n, z, f, point_coordinates, vertex_coordinates, gauge_neighbour_id)
     283           
    300284        return z
    301285
     286   
    302287
    303288    ##
     
    361346                point_coordinates = self._point_coordinates
    362347            else:
    363                 #There are no good point_coordinates. import sys; sys.exit()
     348                # There are no good point_coordinates. import sys; sys.exit()
    364349                msg = 'ERROR (interpolate.py): No point_coordinates inputted'
    365350                raise Exception(msg)
     
    373358                                           verbose=verbose)
    374359            else:
    375                 #Handle blocking
     360                # Handle blocking
    376361                self._A_can_be_reused = False
    377362                start = 0
     
    588573
    589574
     575       
     576       
     577       
     578       
     579def _interpolate_polyline_aux(n, z, f, point_coordinates, vertex_coordinates, gauge_neighbour_id):
     580    """Auxiliary function
     581    """
     582    for i in range(len(point_coordinates)):
     583   
     584        x2 = point_coordinates[i][0]
     585        y2 = point_coordinates[i][1]
     586   
     587        found = False
     588        for j in range(n):
     589           
     590            neighbour_id = gauge_neighbour_id[j]
     591            if neighbour_id >= 0:
     592                x0 = vertex_coordinates[j][0]
     593                y0 = vertex_coordinates[j][1]
     594                x1 = vertex_coordinates[neighbour_id][0]
     595                y1 = vertex_coordinates[neighbour_id][1]
     596           
     597                if point_on_line([x2,y2],
     598                                 [[x0, y0], [x1, y1]],
     599                                 rtol=1.0e-6):
     600                                 
     601                    found = True
     602                    segment_len = sqrt((x1-x0)**2 + (y1-y0)**2)
     603                    dist = sqrt((x2-x0)**2 + (y2-y0)**2)
     604                    z[i] = (f[neighbour_id] - f[j])*dist/segment_len + f[j]
     605                    break
     606                                 
     607        if not found:
     608            z[i] = 0.0               
     609       
     610       
    590611##
    591612# @brief ??
Note: See TracChangeset for help on using the changeset viewer.