Changeset 1670


Ignore:
Timestamp:
Aug 2, 2005, 2:27:25 PM (19 years ago)
Author:
ole
Message:

Refactored Interpolation_function and tested

Location:
inundation/ga/storm_surge/pyvolution
Files:
7 edited

Legend:

Unmodified
Added
Removed
  • inundation/ga/storm_surge/pyvolution/least_squares.py

    r1653 r1670  
    412412        point_coordinates = ensure_numeric(point_coordinates, Float)
    413413
     414        #Keep track of discarded points (if any).
     415        #This is only registered if precrop is True
     416        self.cropped_points = False
    414417
    415418        #Shift data points to same origin as mesh (if specified)
     
    463466            indices = inside_polygon(point_coordinates, P, verbose = verbose)
    464467
    465             if verbose:
    466                 print 'Done'
    467                 if len(indices) != point_coordinates.shape[0]:
    468                     print '%d points outside mesh have been cropped.'\
     468
     469            if len(indices) != point_coordinates.shape[0]:
     470                self.cropped_points = True
     471                if verbose:
     472                    print 'Done - %d points outside mesh have been cropped.'\
    469473                          %(point_coordinates.shape[0] - len(indices))
     474
    470475            point_coordinates = take(point_coordinates, indices)
    471476            self.point_indices = indices
     
    870875
    871876
     877
     878
     879class Interpolation_function:
     880    """Interpolation_function - creates callable object f(t, id) or f(t,x,y)
     881    which is interpolated from time series defined at vertices of
     882    triangular mesh (such as those stored in sww files)
     883
     884    Let m be the number of vertices, n the number of triangles
     885    and p the number of timesteps.
     886
     887    Mandatory input
     888        time:               px1 array of monotonously increasing times (Float)
     889        quantities:         Dictionary of pxm arrays or 1 pxm array (Float)
     890       
     891    Optional input:
     892        quantity_names:     List of keys into quantities dictionary
     893        vertex_coordinates: mx2 array of coordinates (Float)
     894        triangles:          nx3 array of indices into vertex_coordinates (Int)
     895        interpolation_points: array of coordinates to be interpolated to
     896        verbose:            Level of reporting
     897   
     898   
     899    The quantities returned by the callable object are specified by
     900    the list quantities which must contain the names of the
     901    quantities to be returned and also reflect the order, e.g. for
     902    the shallow water wave equation, on would have
     903    quantities = ['stage', 'xmomentum', 'ymomentum']
     904
     905    The parameter interpolation_points decides at which points interpolated
     906    quantities are to be computed whenever object is called.
     907    If None, return average value
     908    """
     909
     910   
     911   
     912    def __init__(self,
     913                 time,
     914                 quantities,
     915                 quantity_names = None, 
     916                 vertex_coordinates = None,
     917                 triangles = None,
     918                 interpolation_points = None,
     919                 verbose = False):
     920        """Initialise object and build spatial interpolation if required
     921        """
     922
     923        from Numeric import array, zeros, Float, alltrue, concatenate,\
     924             reshape, ArrayType
     925
     926        from util import mean, ensure_numeric
     927        from config import time_format
     928        import types
     929
     930
     931
     932        #Check temporal info
     933        time = ensure_numeric(time)       
     934        msg = 'Time must be a monotonuosly '
     935        msg += 'increasing sequence %s' %time
     936        assert alltrue(time[1:] - time[:-1] > 0 ), msg
     937
     938
     939        #Check if quantities is a single array only
     940        if type(quantities) != types.DictType:
     941            quantities = ensure_numeric(quantities)
     942            quantity_names = ['Attribute']
     943
     944            #Make it a dictionary
     945            quantities = {quantity_names[0]: quantities}
     946
     947
     948        #Use keys if no names are specified
     949        if quantity_names is not None:
     950            self.quantity_names = quantity_names
     951        else:
     952            self.quantity_names = quantities.keys()
     953
     954
     955        #Check spatial info
     956        if vertex_coordinates is None:
     957            self.spatial = False
     958        else:   
     959            vertex_coordinates = ensure_numeric(vertex_coordinates)
     960
     961            assert triangles is not None, 'Triangles array must be specified'
     962            triangles = ensure_numeric(triangles)
     963            self.spatial = True           
     964           
     965 
     966        #     
     967        self.interpolation_points = interpolation_points #FIXWME Needed?
     968        self.T = time[:]  #Time assumed to be relative to starttime
     969        self.index = 0    #Initial time index
     970        self.precomputed_values = {}
     971           
     972
     973
     974        #Precomputed spatial interpolation if requested
     975        if interpolation_points is not None:
     976            if self.spatial is False:
     977                raise 'Triangles and vertex_coordinates must be specified'
     978           
     979
     980            try:
     981                interpolation_points = ensure_numeric(interpolation_points)
     982            except:
     983                msg = 'Interpolation points must be an N x 2 Numeric array '+\
     984                      'or a list of points\n'
     985                msg += 'I got: %s.' %( str(interpolation_points)[:60] + '...')
     986                raise msg
     987
     988
     989            for name in quantity_names:
     990                self.precomputed_values[name] =\
     991                                              zeros((len(self.T),
     992                                                     len(interpolation_points)),
     993                                                    Float)
     994
     995            #Build interpolator
     996            interpol = Interpolation(vertex_coordinates,
     997                                     triangles,
     998                                     point_coordinates = interpolation_points,
     999                                     alpha = 0,
     1000                                     precrop = False,
     1001                                     verbose = verbose)
     1002
     1003            #if interpol.cropped_points is True:
     1004            #    raise 'Some interpolation points were outside mesh'
     1005            #FIXME: This will be raised if triangles are listed as
     1006            #discontinuous even though there is no need to stop
     1007            #(precrop = True above)
     1008
     1009            if verbose: print 'Interpolate'
     1010            for i, t in enumerate(self.T):
     1011                #Interpolate quantities at this timestep
     1012                if verbose: print ' time step %d of %d' %(i, len(self.T))
     1013                for name in quantity_names:
     1014                    self.precomputed_values[name][i, :] =\
     1015                    interpol.interpolate(quantities[name][i,:])
     1016
     1017            #Report
     1018            if verbose:
     1019                x = vertex_coordinates[:,0]
     1020                y = vertex_coordinates[:,1]               
     1021           
     1022                print '------------------------------------------------'
     1023                print 'Interpolation_function statistics:'
     1024                print '  Extent:'
     1025                print '    x in [%f, %f], len(x) == %d'\
     1026                      %(min(x), max(x), len(x))
     1027                print '    y in [%f, %f], len(y) == %d'\
     1028                      %(min(y), max(y), len(y))
     1029                print '    t in [%f, %f], len(t) == %d'\
     1030                      %(min(self.T), max(self.T), len(self.T))
     1031                print '  Quantities:'
     1032                for name in quantity_names:
     1033                    q = quantities[name][:].flat
     1034                    print '    %s in [%f, %f]' %(name, min(q), max(q))
     1035                print '  Interpolation points (xi, eta):'\
     1036                      ' number of points == %d ' %interpolation_points.shape[0]
     1037                print '    xi in [%f, %f]' %(min(interpolation_points[:,0]),
     1038                                             max(interpolation_points[:,0]))
     1039                print '    eta in [%f, %f]' %(min(interpolation_points[:,1]),
     1040                                              max(interpolation_points[:,1]))
     1041                print '  Interpolated quantities (over all timesteps):'
     1042               
     1043                for name in quantity_names:
     1044                    q = self.precomputed_values[name][:].flat
     1045                    print '    %s at interpolation points in [%f, %f]'\
     1046                          %(name, min(q), max(q))
     1047                print '------------------------------------------------'
     1048           
     1049        else:
     1050            #Store quantitites as is
     1051            for name in quantity_names:
     1052                self.precomputed_values[name] = quantities[name]
     1053
     1054
     1055        #else:
     1056        #    #Return an average, making this a time series
     1057        #    for name in quantity_names:
     1058        #        self.values[name] = zeros(len(self.T), Float)
     1059        #
     1060        #    if verbose: print 'Compute mean values'
     1061        #    for i, t in enumerate(self.T):
     1062        #        if verbose: print ' time step %d of %d' %(i, len(self.T))
     1063        #        for name in quantity_names:
     1064        #           self.values[name][i] = mean(quantities[name][i,:])
     1065
     1066
     1067
     1068
     1069    def __repr__(self):
     1070        return 'Interpolation function (spation-temporal)'
     1071
     1072    def __call__(self, t, point_id = None, x = None, y = None):
     1073        """Evaluate f(t), f(t, point_id) or f(t, x, y)
     1074
     1075        Inputs:
     1076          t: time - Model time. Must lie within existing timesteps
     1077          point_id: index of one of the preprocessed points.
     1078          x, y:     Overrides location, point_id ignored
     1079         
     1080          If spatial info is present and all of x,y,point_id
     1081          are None an exception is raised
     1082                   
     1083          If no spatial info is present, point_id and x,y arguments are ignored
     1084          making f a function of time only.
     1085
     1086         
     1087          FIXME: point_id could also be a slice
     1088          FIXME: What if x and y are vectors?
     1089          FIXME: What about f(x,y) without t?
     1090        """
     1091
     1092        from math import pi, cos, sin, sqrt
     1093        from Numeric import zeros, Float
     1094        from util import mean       
     1095
     1096        if self.spatial is True:
     1097            if point_id is None:
     1098                if x is None or y is None:
     1099                    msg = 'Either point_id or x and y must be specified'
     1100                    raise msg
     1101            else:
     1102                if self.interpolation_points is None:
     1103                    msg = 'Interpolation_function must be instantiated ' +\
     1104                          'with a list of interpolation points before parameter ' +\
     1105                          'point_id can be used'
     1106                    raise msg
     1107
     1108
     1109        msg = 'Time interval [%s:%s]' %(self.T[0], self.T[1])
     1110        msg += ' does not match model time: %s\n' %t
     1111        if t < self.T[0]: raise msg
     1112        if t > self.T[-1]: raise msg
     1113
     1114        oldindex = self.index #Time index
     1115
     1116        #Find current time slot
     1117        while t > self.T[self.index]: self.index += 1
     1118        while t < self.T[self.index]: self.index -= 1
     1119
     1120        if t == self.T[self.index]:
     1121            #Protect against case where t == T[-1] (last time)
     1122            # - also works in general when t == T[i]
     1123            ratio = 0
     1124        else:
     1125            #t is now between index and index+1
     1126            ratio = (t - self.T[self.index])/\
     1127                    (self.T[self.index+1] - self.T[self.index])
     1128
     1129        #Compute interpolated values
     1130        q = zeros(len(self.quantity_names), Float)
     1131
     1132        for i, name in enumerate(self.quantity_names):
     1133            Q = self.precomputed_values[name]
     1134
     1135            if self.spatial is False:
     1136                #If there is no spatial info               
     1137                assert len(Q.shape) == 1
     1138
     1139                Q0 = Q[self.index]
     1140                if ratio > 0: Q1 = Q[self.index+1]
     1141
     1142            else:
     1143                if x is not None and y is not None:
     1144                    #Interpolate to x, y
     1145                   
     1146                    raise 'x,y interpolation not yet implemented'
     1147                else:
     1148                    #Use precomputed point
     1149                    Q0 = Q[self.index, point_id]
     1150                    if ratio > 0: Q1 = Q[self.index+1, point_id]
     1151
     1152            #Linear temporal interpolation   
     1153            if ratio > 0:
     1154                q[i] = Q0 + ratio*(Q1 - Q0)
     1155            else:
     1156                q[i] = Q0
     1157
     1158            #if ratio > 0:
     1159            #    q[i] = Q[self.index, point_id] +\
     1160            #           ratio*(Q[self.index+1, point_id] -\
     1161            #                  Q[self.index, point_id])
     1162            #else:
     1163            #    q[i] = Q[self.index, point_id]
     1164
     1165        #Return vector of interpolated values
     1166        if len(q) == 1:
     1167            return q[0]
     1168        else:
     1169            return q
     1170
     1171
     1172
     1173
     1174
    8721175#-------------------------------------------------------------
    8731176if __name__ == "__main__":
  • inundation/ga/storm_surge/pyvolution/mesh.py

    r1575 r1670  
    416416    def get_boundary_polygon(self):
    417417        """Return bounding polygon as a list of points
     418
     419        FIXME: If triangles are listed as discontinuous
     420        (e.g vertex coordinates listed multiple times),
     421        this may not work as expected.
    418422        """
    419423        from Numeric import allclose, sqrt, array, minimum, maximum
     
    421425
    422426
    423         V = self.get_vertex_coordinates()
     427        #V = self.get_vertex_coordinates()
    424428        segments = {}
    425429
  • inundation/ga/storm_surge/pyvolution/test_least_squares.py

    r1657 r1670  
    989989
    990990
     991    def test_interpolation_function_time_only(self):
     992        """Test spatio-temporal interpolation
     993        Test that spatio temporal function performs the correct
     994        interpolations in both time and space
     995        """
     996
     997
     998        #Three timesteps
     999        time = [1.0, 5.0, 6.0]
     1000       
     1001
     1002        #One quantity
     1003        Q = zeros( (3,6), Float )
     1004
     1005        #Linear in time and space
     1006        a = [0.0, 0.0]
     1007        b = [0.0, 2.0]
     1008        c = [2.0, 0.0]
     1009        d = [0.0, 4.0]
     1010        e = [2.0, 2.0]
     1011        f = [4.0, 0.0]
     1012
     1013        points = [a, b, c, d, e, f]
     1014       
     1015        for i, t in enumerate(time):
     1016            Q[i, :] = t*linear_function(points)
     1017
     1018           
     1019        #Check basic interpolation of one quantity using averaging
     1020        #(no interpolation points or spatial info)
     1021        from util import mean       
     1022        I = Interpolation_function(time, [mean(Q[0,:]),
     1023                                          mean(Q[1,:]),
     1024                                          mean(Q[2,:])])
     1025
     1026
     1027
     1028        #Check temporal interpolation
     1029        for i in [0,1,2]:
     1030            assert allclose(I(time[i]), mean(Q[i,:]))
     1031
     1032        #Midway   
     1033        assert allclose(I( (time[0] + time[1])/2 ),
     1034                        (I(time[0]) + I(time[1]))/2 )
     1035
     1036        assert allclose(I( (time[1] + time[2])/2 ),
     1037                        (I(time[1]) + I(time[2]))/2 )
     1038
     1039        assert allclose(I( (time[0] + time[2])/2 ),
     1040                        (I(time[0]) + I(time[2]))/2 )                 
     1041
     1042        #1/3
     1043        assert allclose(I( (time[0] + time[2])/3 ),
     1044                        (I(time[0]) + I(time[2]))/3 )                         
     1045
     1046
     1047        #Out of bounds checks
     1048        try:
     1049            I(time[0]-1)
     1050        except:
     1051            pass
     1052        else:
     1053            raise 'Should raise exception'
     1054
     1055        try:
     1056            I(time[-1]+1)
     1057        except:
     1058            pass
     1059        else:
     1060            raise 'Should raise exception'       
     1061
     1062
     1063       
     1064
     1065    def test_interpolation_function(self):
     1066        """Test spatio-temporal interpolation
     1067        Test that spatio temporal function performs the correct
     1068        interpolations in both time and space
     1069        """
     1070
     1071
     1072        #Three timesteps
     1073        time = [1.0, 5.0, 6.0]
     1074       
     1075
     1076        #Setup mesh used to represent fitted function
     1077        a = [0.0, 0.0]
     1078        b = [0.0, 2.0]
     1079        c = [2.0, 0.0]
     1080        d = [0.0, 4.0]
     1081        e = [2.0, 2.0]
     1082        f = [4.0, 0.0]
     1083
     1084        points = [a, b, c, d, e, f]
     1085        #bac, bce, ecf, dbe
     1086        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
     1087
     1088
     1089        #New datapoints where interpolated values are sought
     1090        interpolation_points = [[ 0.0, 0.0],
     1091                                [ 0.5, 0.5],
     1092                                [ 0.7, 0.7],
     1093                                [ 1.0, 0.5],
     1094                                [ 2.0, 0.4],
     1095                                [ 2.8, 1.2]]
     1096
     1097
     1098        #One quantity
     1099        Q = zeros( (3,6), Float )
     1100
     1101        #Linear in time and space
     1102        for i, t in enumerate(time):
     1103            Q[i, :] = t*linear_function(points)
     1104
     1105
     1106        #Check interpolation of one quantity using interpolaton points)
     1107        I = Interpolation_function(time, Q,
     1108                                   vertex_coordinates = points,
     1109                                   triangles = triangles,
     1110                                   interpolation_points = interpolation_points,
     1111                                   verbose = False)
     1112
     1113
     1114        answer = linear_function(interpolation_points)
     1115
     1116        t = time[0]
     1117        for j in range(50): #t in [1, 6]
     1118            for id in range(len(interpolation_points)):
     1119                assert allclose(I(t, id), t*answer[id])
     1120
     1121            t += 0.1   
     1122
     1123
     1124        try:   
     1125            I(1)
     1126        except:
     1127            pass
     1128        else:
     1129            raise 'Should raise exception'
     1130           
     1131        #
     1132        #interpolation_points = [[ 0.0, 0.0],
     1133        #                        [ 0.5, 0.5],
     1134        #                        [ 0.7, 0.7],
     1135        #                        [-13, 65],  #Outside
     1136        #                        [ 1.0, 0.5],
     1137        #                        [ 2.0, 0.4],
     1138        #                        [ 2.8, 1.2]]
     1139        #
     1140        #try:
     1141        #    I = Interpolation_function(time, Q,
     1142        #                               vertex_coordinates = points,
     1143        #                               triangles = triangles,
     1144        #                               interpolation_points = interpolation_points,
     1145        #                               verbose = False)
     1146        #except:
     1147        #    pass
     1148        #else:
     1149        #    raise 'Should raise exception'
    9911150
    9921151
  • inundation/ga/storm_surge/pyvolution/test_mesh.py

    r1632 r1670  
    800800
    801801
     802    def test_boundary_polygon_IV(self):
     803        """Reproduce test test_spatio_temporal_file_function_time
     804        from test_util.py that looked as if it produced the wrong boundary
     805        """
     806
     807        from mesh import Mesh
     808        from Numeric import zeros, Float
     809        from util import inside_polygon
     810        from mesh_factory import rectangular       
     811
     812        #Create a domain to hold test grid
     813        #(0:15, -20:10)
     814        points, vertices, boundary =\
     815                rectangular(4, 4, 15, 30, origin = (0, -20))       
     816
     817        #####
     818        mesh = Mesh(points, vertices)
     819        mesh.check_integrity()
     820
     821        P = mesh.get_boundary_polygon()
     822
     823        #print P
     824        assert len(P) == 16
     825        for p in points:
     826            assert inside_polygon(p, P)
     827
     828
     829
     830        #####
     831        mesh = Mesh(points, vertices, boundary)
     832        mesh.check_integrity()
     833
     834        P = mesh.get_boundary_polygon()
     835
     836       
     837        #print P, len(P)
     838        assert len(P) == 16
     839
     840        for p in points:
     841            assert inside_polygon(p, P)           
     842
     843
     844
     845
    802846
    803847#-------------------------------------------------------------
  • inundation/ga/storm_surge/pyvolution/test_util.py

    r1668 r1670  
    495495
    496496        #Create a domain to hold test grid
    497 
     497        #(0:15, -20:10)
    498498        points, vertices, boundary =\
    499499                rectangular(4, 4, 15, 30, origin = (0, -20))
     
    566566                          quantities = domain.conserved_quantities,
    567567                          interpolation_points = interpolation_points)               
    568         #for id in range(len(interpolation_points)):
    569         for id in [1]:
     568        for id in range(len(interpolation_points)):
    570569            x = interpolation_points[id][0]
    571570            y = interpolation_points[id][1]
     
    617616
    618617        #Now try interpolation with delta offset
    619         for id in [1]:
     618        for id in range(len(interpolation_points)):           
    620619            x = interpolation_points[id][0]
    621620            y = interpolation_points[id][1]
  • inundation/ga/storm_surge/pyvolution/util.py

    r1668 r1670  
    291291   
    292292
    293     return Interpolation_function(vertex_coordinates,
     293    from least_squares import Interpolation_function
     294    return Interpolation_function(time,
     295                                  quantities,
     296                                  quantity_names,
     297                                  vertex_coordinates,
    294298                                  triangles,
    295                                   time,
    296                                   quantities,
    297                                   quantity_names, 
    298299                                  interpolation_points,
    299                                   alpha = 0,
    300300                                  verbose = verbose)
    301301
     
    303303
    304304
    305 
    306 class Interpolation_function:
    307     """Interpolation_function - creates callable object f(t, id) or f(t,x,y)
    308     which is interpolated from time series defined at vertices of
    309     triangular mesh
    310 
    311     Input
    312         vertex_coordinates: mx2 array of coordinates (Float)
    313         triangles:          nx3 array of indices into vertex_coordinates (Int)
    314         time:               px1 array of times (Float)               
    315         quantities:         Dictionary of pxm values
    316         quantity_names:
    317         interpolation_points:
    318         alpha:
    319         verbose:
    320    
    321 
    322    
    323     The quantities returned by the callable object are specified by
    324     the list quantities which must contain the names of the
    325     quantities to be returned and also reflect the order, e.g. for
    326     the shallow water wave equation, on would have
    327     quantities = ['stage', 'xmomentum', 'ymomentum']
    328 
    329     The parameter interpolation_points decides at which points interpolated
    330     quantities are to be computed whenever object is called.
    331     If None, return average value
    332     """
    333 
    334    
    335 
    336    
    337     def __init__(self,
    338                  vertex_coordinates,
    339                  triangles,
    340                  time,
    341                  quantities,
    342                  quantity_names = None, 
    343                  interpolation_points = None,
    344                  alpha = 0,
    345                  verbose = False):
    346         """
    347         """
    348 
    349         from Numeric import array, zeros, Float, alltrue, concatenate,\
    350              reshape, ArrayType
    351 
    352         from config import time_format
    353         from least_squares import Interpolation
    354 
    355         #Check
    356         msg = 'Time must be a monotonuosly '
    357         msg += 'increasing sequence %s' %time
    358         assert alltrue(time[1:] - time[:-1] > 0 ), msg
    359 
    360 
    361         #Check if quantities is a single array only
    362         if type(quantities) == ArrayType:
    363             quantity_names = 'Attribute'
    364             quantities = {quantity_names: quantities}           
    365 
    366         #Use keys if no names specified
    367         if quantity_names is not None:
    368             self.quantity_names = quantity_names
    369         else:
    370             self.quantity_names = quantities.keys()
    371 
    372 
    373         self.interpolation_points = interpolation_points
    374         self.T = time[:]  #Time assumed to be relative to starttime
    375         self.index = 0    #Initial time index
    376         self.values = {}
    377            
    378 
    379         if interpolation_points is not None:
    380             #Spatial interpolation
    381 
    382             try:
    383                 interpolation_points = ensure_numeric(interpolation_points)
    384             except:
    385                 msg = 'Interpolation points must be an N x 2 Numeric array or a list of points\n'
    386                 msg += 'I got: %s.' %( str(interpolation_points)[:60] + '...')
    387                 raise msg
    388 
    389 
    390             for name in quantity_names:
    391                 self.values[name] = zeros( (len(self.T),
    392                                             len(interpolation_points)),
    393                                             Float)
    394 
    395             #Build interpolator
    396             interpol = Interpolation(vertex_coordinates,
    397                                      triangles,
    398                                      point_coordinates = interpolation_points,
    399                                      alpha = 0,
    400                                      verbose = verbose)
    401 
    402 
    403             if verbose: print 'Interpolate'
    404             for i, t in enumerate(self.T):
    405                 #Interpolate quantities at this timestep
    406                 if verbose: print ' time step %d of %d' %(i, len(self.T))
    407                 for name in quantity_names:
    408                     self.values[name][i, :] =\
    409                     interpol.interpolate(quantities[name][i,:])
    410 
    411             #Report
    412             if verbose:
    413                 x = vertex_coordinates[:,0]
    414                 y = vertex_coordinates[:,1]               
    415                
    416                 print '------------------------------------------------'
    417                 print 'Interpolation_function statistics:'
    418                 print '  Extent:'
    419                 print '    x in [%f, %f], len(x) == %d'\
    420                       %(min(x.flat), max(x.flat), len(x.flat))
    421                 print '    y in [%f, %f], len(y) == %d'\
    422                       %(min(y.flat), max(y.flat), len(y.flat))
    423                 print '    t in [%f, %f], len(t) == %d'\
    424                       %(min(self.T), max(self.T), len(self.T))
    425                 print '  Quantities:'
    426                 for name in quantity_names:
    427                     q = quantities[name][:].flat
    428                     print '    %s in [%f, %f]' %(name, min(q), max(q))
    429 
    430 
    431                 print '  Interpolation points (xi, eta):'\
    432                       'number of points == %d ' %interpolation_points.shape[0]
    433                 print '    xi in [%f, %f]' %(min(interpolation_points[:,0]),
    434                                              max(interpolation_points[:,0]))
    435                 print '    eta in [%f, %f]' %(min(interpolation_points[:,1]),
    436                                               max(interpolation_points[:,1]))
    437                 print '  Interpolated quantities (over all timesteps):'
    438                 for name in quantity_names:
    439                     q = self.values[name][:].flat
    440                     print '    %s at interpolation points in [%f, %f]'\
    441                           %(name, min(q), max(q))
    442                 print '------------------------------------------------'
    443         else:
    444             #Return an average, making this a time series
    445             #FIXME: Needs a test
    446             for name in quantity_names:
    447                 self.values[name] = avg(quantities[name][i,:])
    448 
    449 
    450     def __repr__(self):
    451         return 'Interpolation function'
    452 
    453     def __call__(self, t, point_id = None, x = None, y = None):
    454         """Evaluate f(t), f(t, point_id) or f(t, x, y)
    455 
    456         Inputs:
    457           t: time - Model time. Must lie within existing timesteps
    458           point_id: index of one of the preprocessed points.
    459                     If point_id is None all preprocessed points are computed
    460                    
    461           If no x,y info is present, point_id and x,y arguments are ignored
    462           making f a function of time only.
    463          
    464           FIXME: point_id could also be a slice
    465           FIXME: One could allow arbitrary x, y coordinates
    466           (requires computation of a new interpolator)
    467           FIXME: What if x and y are vectors?
    468           FIXME: Does point_id or x,y take precedence?
    469           FIXME: What about f(x,y) without t
    470          
    471         """
    472 
    473         from math import pi, cos, sin, sqrt
    474         from Numeric import zeros, Float
    475 
    476 
    477         msg = 'Time interval [%s:%s]' %(self.T[0], self.T[1])
    478         msg += ' does not match model time: %s\n' %t
    479         if t < self.T[0]: raise msg
    480         if t > self.T[-1]: raise msg
    481 
    482         oldindex = self.index #Time index
    483 
    484         #Find current time slot
    485         while t > self.T[self.index]: self.index += 1
    486         while t < self.T[self.index]: self.index -= 1
    487 
    488         if t == self.T[self.index]:
    489             #Protect against case where t == T[-1] (last time)
    490             # - also works in general when t == T[i]
    491             ratio = 0
    492         else:
    493             #t is now between index and index+1
    494             ratio = (t - self.T[self.index])/\
    495                     (self.T[self.index+1] - self.T[self.index])
    496 
    497         #Compute interpolated values
    498         q = zeros(len(self.quantity_names), Float)
    499 
    500         for i, name in enumerate(self.quantity_names):
    501             Q = self.values[name]
    502 
    503             if point_id is None:
    504                 Q0 = Q[self.index]
    505 
    506                 try:
    507                     Q1 = Q[self.index+1]
    508                 except:
    509                     pass               
    510             else:   
    511                 Q0 = Q[self.index, point_id]
    512 
    513                 try:
    514                     Q1 = Q[self.index+1, point_id]
    515                 except:
    516                     pass
    517                
    518 
    519             #Linear temporal interpolation   
    520             if ratio > 0:
    521                 q[i] = Q0 + ratio*(Q1 - Q0)
    522             else:
    523                 q[i] = Q0
    524            
    525 
    526             #if ratio > 0:
    527             #    q[i] = Q[self.index, point_id] +\
    528             #           ratio*(Q[self.index+1, point_id] -\
    529             #                  Q[self.index, point_id])
    530             #else:
    531             #    q[i] = Q[self.index, point_id]
    532 
    533         #Return vector of interpolated values
    534         return q
    535305
    536306
  • inundation/ga/storm_surge/pyvolution/wiki/issues.txt

    r1652 r1670  
    88 are outside the mesh.
    99Importance: Mid
     10Comment (Ole): If points and triangles are listed as discontinuous e.g. vertex coordinates replicated for each triangle, then the bounding polygon gets mangled and some points appear to be outside.
    1011Who: DSG
    1112
Note: See TracChangeset for help on using the changeset viewer.