Changeset 5222


Ignore:
Timestamp:
Apr 21, 2008, 5:13:27 PM (15 years ago)
Author:
ole
Message:

Work done during Water Down Under 2008.
Input checks in interpolate_block and comments.

Location:
anuga_core/source/anuga
Files:
2 edited

Legend:

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

    r5196 r5222  
    120120        """
    121121
    122        
    123         # FIXME (Ole): Need an input check that dimensions are compatible
    124 
    125122        # FIXME (Ole): Why is the interpolation matrix rebuilt everytime the
    126123        # method is called even if interpolation points are unchanged.
    127        
     124
     125
    128126        #print "point_coordinates interpolate.interpolate", point_coordinates
    129127        if isinstance(point_coordinates, Geospatial_data):
     
    146144                msg = 'ERROR (interpolate.py): No point_coordinates inputted'
    147145                raise Exception(msg)
     146
    148147       
    149148        if point_coordinates is not None:   
     
    184183                z = concatenate((z,t))
    185184        return z
    186 
    187     def interpolate_block(self, f, point_coordinates=None, verbose=False):
     185   
     186
     187    def interpolate_block(self, f, point_coordinates, verbose=False):
    188188        """
    189189        Call this if you want to control the blocking or make sure blocking
     
    195195        """
    196196        if isinstance(point_coordinates,Geospatial_data):
    197             point_coordinates = point_coordinates.get_data_points( \
    198                 absolute = True)
    199         if point_coordinates is not None:
    200             self._A = self._build_interpolation_matrix_A(point_coordinates,
    201                                                          verbose=verbose)
     197            point_coordinates = point_coordinates.get_data_points(\
     198                absolute=True)
     199
     200        # Convert lists to Numeric arrays if necessary
     201        point_coordinates = ensure_numeric(point_coordinates, Float)
     202        f = ensure_numeric(f, Float)       
     203           
     204
     205        self._A = self._build_interpolation_matrix_A(point_coordinates,
     206                                                     verbose=verbose)
     207
     208
     209        # Check that input dimensions are compatible
     210        msg = 'Two colums must be specified in point coordinates. I got shape=%s'\
     211              %(str(point_coordinates.shape))
     212        assert point_coordinates.shape[1] == 2, msg
     213
     214        msg = 'The number of rows in matrix A must be the same as the number of points supplied.'
     215        msg += ' I got %d points and %d matrix rows.'\
     216               %(point_coordinates.shape[0], self._A.shape[0])
     217        assert point_coordinates.shape[0] == self._A.shape[0], msg       
     218
     219        msg = 'The number of columns in matrix A must be the same as the number of mesh vertices.'
     220        msg += ' I got %d vertices and %d matrix columns.'\
     221               %(f.shape[0], self._A.shape[1])       
     222        assert self._A.shape[1] == f.shape[0], msg
     223
     224        # Compute Matrix vector product and return
    202225        return self._get_point_data_z(f)
     226   
    203227
    204228    def _get_point_data_z(self, f, verbose=False):
     
    209233        The _A matrix has been created
    210234        """
     235
    211236        z = self._A * f
    212237        # Taking into account points outside the mesh.
     
    238263
    239264        #print 'Building interpolation matrix'
    240        
    241         #Convert point_coordinates to Numeric arrays, in case it was a list.
     265
     266        # Convert point_coordinates to Numeric arrays, in case it was a list.
    242267        point_coordinates = ensure_numeric(point_coordinates, Float)
     268       
    243269       
    244270        if verbose: print 'Getting indices inside mesh boundary'
     
    624650                # Interpolate quantities at this timestep
    625651                if verbose and i%((p+10)/10)==0:
    626                     print ' time step %d of %d' %(i, p)
     652                    print '  time step %d of %d' %(i, p)
    627653                   
    628654                for name in quantity_names:
     
    631657                    else:
    632658                        Q = quantities[name][:]   # No time dependency
    633                        
     659
     660                    if verbose and i%((p+10)/10)==0:
     661                        print '    quantity %s, size=%d' %(name, len(Q))
     662                       
    634663                    # Interpolate   
    635664                    result = interpol.interpolate(Q,
  • anuga_core/source/anuga/geospatial_data/geospatial_data.py

    r5146 r5222  
    994994                file_pointer,
    995995                header,
    996                 max_read_lines=1e30) #If the file is bigger that this, block..
     996                max_read_lines=1e30) #If the file is bigger that this, block..  # FIXME (Ole) What's up here?
     997               
    997998    except ANUGAError:
    998999        file_pointer.close()
Note: See TracChangeset for help on using the changeset viewer.