Changeset 5222

Timestamp:

Apr 21, 2008, 5:13:27 PM (16 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:

• fit_interpolate/interpolate.py (modified) (8 diffs)
• geospatial_data/geospatial_data.py (modified) (1 diff)

anuga_core/source/anuga/fit_interpolate/interpolate.py

@@ -120 +120 @@
         """
 
-        
-        # FIXME (Ole): Need an input check that dimensions are compatible
-
         # FIXME (Ole): Why is the interpolation matrix rebuilt everytime the
         # method is called even if interpolation points are unchanged.
-        
+
+
         #print "point_coordinates interpolate.interpolate", point_coordinates 
         if isinstance(point_coordinates, Geospatial_data):
@@ -146 +144 @@
                 msg = 'ERROR (interpolate.py): No point_coordinates inputted'
                 raise Exception(msg)
+
         
         if point_coordinates is not None:   
@@ -184 +183 @@
                 z = concatenate((z,t))
         return z
-
-    def interpolate_block(self, f, point_coordinates=None, verbose=False):
+    
+
+    def interpolate_block(self, f, point_coordinates, verbose=False):
         """
         Call this if you want to control the blocking or make sure blocking
@@ -195 +195 @@
         """
         if isinstance(point_coordinates,Geospatial_data):
-            point_coordinates = point_coordinates.get_data_points( \
-                absolute = True)
-        if point_coordinates is not None:
-            self._A = self._build_interpolation_matrix_A(point_coordinates,
-                                                         verbose=verbose)
+            point_coordinates = point_coordinates.get_data_points(\
+                absolute=True)
+
+        # Convert lists to Numeric arrays if necessary
+        point_coordinates = ensure_numeric(point_coordinates, Float)
+        f = ensure_numeric(f, Float)        
+            
+
+        self._A = self._build_interpolation_matrix_A(point_coordinates,
+                                                     verbose=verbose)
+
+
+        # Check that input dimensions are compatible
+        msg = 'Two colums must be specified in point coordinates. I got shape=%s'\
+              %(str(point_coordinates.shape))
+        assert point_coordinates.shape[1] == 2, msg
+
+        msg = 'The number of rows in matrix A must be the same as the number of points supplied.'
+        msg += ' I got %d points and %d matrix rows.'\
+               %(point_coordinates.shape[0], self._A.shape[0])
+        assert point_coordinates.shape[0] == self._A.shape[0], msg        
+
+        msg = 'The number of columns in matrix A must be the same as the number of mesh vertices.'
+        msg += ' I got %d vertices and %d matrix columns.'\
+               %(f.shape[0], self._A.shape[1])        
+        assert self._A.shape[1] == f.shape[0], msg
+
+        # Compute Matrix vector product and return
         return self._get_point_data_z(f)
+    
 
     def _get_point_data_z(self, f, verbose=False):
@@ -209 +233 @@
         The _A matrix has been created
         """
+
         z = self._A * f
         # Taking into account points outside the mesh.
@@ -238 +263 @@
 
         #print 'Building interpolation matrix'
-        
-        #Convert point_coordinates to Numeric arrays, in case it was a list.
+
+        # Convert point_coordinates to Numeric arrays, in case it was a list.
         point_coordinates = ensure_numeric(point_coordinates, Float)
+        
         
         if verbose: print 'Getting indices inside mesh boundary'
@@ -624 +650 @@
                 # Interpolate quantities at this timestep
                 if verbose and i%((p+10)/10)==0:
-                    print ' time step %d of %d' %(i, p)
+                    print '  time step %d of %d' %(i, p)
                     
                 for name in quantity_names:
@@ -631 +657 @@
                     else:
                         Q = quantities[name][:]   # No time dependency
-                       
+
+                    if verbose and i%((p+10)/10)==0:
+                        print '    quantity %s, size=%d' %(name, len(Q))
+                        
                     # Interpolate    
                     result = interpol.interpolate(Q,

anuga_core/source/anuga/geospatial_data/geospatial_data.py

@@ -994 +994 @@
                 file_pointer,
                 header,
-                max_read_lines=1e30) #If the file is bigger that this, block..
+                max_read_lines=1e30) #If the file is bigger that this, block..  # FIXME (Ole) What's up here?
+                
     except ANUGAError:
         file_pointer.close()