Index: /inundation/ga/storm_surge/pyvolution/domain.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/domain.py	(revision 850)
+++ /inundation/ga/storm_surge/pyvolution/domain.py	(revision 851)
@@ -270,4 +270,5 @@
                 raise msg
 
+		
     def set_region(self, functions):
         # The order of functions in the list is used.
Index: /inundation/ga/storm_surge/pyvolution/generic_boundary_conditions.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/generic_boundary_conditions.py	(revision 850)
+++ /inundation/ga/storm_surge/pyvolution/generic_boundary_conditions.py	(revision 851)
@@ -187,6 +187,10 @@
 	#any tagged boundary later on.
 	
-        midpoint_coordinates = zeros( (len(domain.boundary), 2), Float)
+        self.midpoint_coordinates = zeros( (len(domain.boundary), 2), Float)
 	boundary_keys = domain.boundary.keys()
+	
+	#Record ordering #FIXME: should this happen in domain.py
+        self.boundary_indices = {} 
+	
         for i, (vol_id, edge_id) in enumerate(boundary_keys):
             
@@ -199,16 +203,19 @@
             if edge_id == 1: m = array([(x0 + x2)/2, (y0 + y2)/2])
             if edge_id == 2: m = array([(x1 + x0)/2, (y1 + y0)/2])            
-            midpoint_coordinates[i,:] = m
-        
-        self.midpoint_coordinates = midpoint_coordinates    
+            self.midpoint_coordinates[i,:] = m
+	    
+	    self.boundary_indices[(vol_id, edge_id)] = i
+        
 	
-        self.F = file_function(filename, domain, interpolation_points=midpoint_coordinates)
+        self.F = file_function(filename, domain, 
+	                       interpolation_points=self.midpoint_coordinates)
         self.domain = domain
 
         #Test
-        q = self.F(0, 0, 0)
+        q = self.F(0, point_id=0)
 
         d = len(domain.conserved_quantities)
-        msg = 'Values specified in file %s must be a list or an array of length %d' %(filename, d)
+        msg = 'Values specified in file %s must be ' %filename 
+	msg += ' a list or an array of length %d' %d
         assert len(q) == d, msg
 
@@ -226,8 +233,8 @@
 
         if vol_id is not None and edge_id is not None:
-            x, y = self.midpoint_coordinates[ vol_id, edge_id ]
-            return self.F(t, x, y)
+            i = self.boundary_indices[ vol_id, edge_id ]
+            return self.F(t, point_id = i)
         else:
-            return self.F(t)
+            return self.F(t) #What should the semantics be?
 
 
Index: /inundation/ga/storm_surge/pyvolution/shallow_water.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/shallow_water.py	(revision 850)
+++ /inundation/ga/storm_surge/pyvolution/shallow_water.py	(revision 851)
@@ -941,4 +941,5 @@
         except Exception, e:
             msg = 'Function %s could not be executed:\n%s' %(f, e)
+	    #FIXME: Reconsider this semantics
             raise msg
                     
Index: /inundation/ga/storm_surge/pyvolution/test_shallow_water.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/test_shallow_water.py	(revision 850)
+++ /inundation/ga/storm_surge/pyvolution/test_shallow_water.py	(revision 851)
@@ -2415,5 +2415,5 @@
         Bd = Dirichlet_boundary([0.3,0,0])	
 	Bt = Transmissive_boundary(domain1)
-        domain1.set_boundary({'left': Bd, 'right': Bt, 'top': Br, 'bottom': Br})
+        domain1.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
 
         #Initial condition
@@ -2452,5 +2452,5 @@
 	domain2.reduction = domain1.reduction
         domain2.smooth = False
-        domain2.visualise = True
+        domain2.visualise = False
         domain2.default_order = 2
 	
@@ -2465,10 +2465,7 @@
 	
 
-	#FIXME: Continue from here!!!!!!!
-	return	
-		
 	Bf = Spatio_temporal_boundary(domain1.filename + '.' + domain1.format,
 				      domain2)
-        domain2.set_boundary({'right':Bt, 'bottom':Br, 'diagonal':Bf})
+        domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
 
         #Initial condition 
@@ -2486,19 +2483,23 @@
         cv2 = domain2.quantities['stage'].centroid_values
 	
-        print take(cv1, (12,14,16))  #Right
-	print take(cv2, (4,6,8))
-	
+        #print take(cv1, (12,14,16))  #Right
+	#print take(cv2, (4,6,8))
         #print take(cv1, (0,6,12))  #Bottom	
-	
-	assert allclose( take(cv1, (12,14,16)), take(cv2, (4,6,8))) #RHS
-	assert allclose( take(cv1, (0,6,12)), take(cv2, (0,1,4))) #Bottom	
+	#print take(cv2, (0,1,4))
+        #print take(cv1, (0,8,16))  #Diag	
+	#print take(cv2, (0,3,8))	
+	#assert allclose( take(cv1, (12,14,16)), take(cv2, (4,6,8))) #RHS
+	#assert allclose( take(cv1, (0,6,12)), take(cv2, (0,1,4))) #Bottom	
 	
         #Cleanup
-        os.remove(domain1.filename + '.' + domain1.format)
-    
+        try:
+            os.remove(domain1.filename + '.' + domain1.format)
+        except:
+	    pass
+	     
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    #suite = unittest.makeSuite(TestCase,'test')
-    suite = unittest.makeSuite(TestCase,'test_spatio_temporal_boundary')
+    suite = unittest.makeSuite(TestCase,'test')
+    #suite = unittest.makeSuite(TestCase,'test_spatio_temporal_boundary')
     runner = unittest.TextTestRunner()
     runner.run(suite)
Index: /inundation/ga/storm_surge/pyvolution/util.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/util.py	(revision 850)
+++ /inundation/ga/storm_surge/pyvolution/util.py	(revision 851)
@@ -123,5 +123,9 @@
     quantities = ['stage', 'xmomentum', 'ymomentum']
     
-    interpolation_points decides at which points interpolated quantities are to be computed
+    interpolation_points decides at which points interpolated 
+    quantities are to be computed whenever object is called.
+    
+    
+    
     If None, return average value
     """
@@ -133,6 +137,4 @@
         will be checked and possibly modified
 	
-
-
         All times are assumed to be in UTC
         """
@@ -213,5 +215,5 @@
         msg = 'File %s must list time as a monotonuosly ' %self.filename
         msg += 'increasing sequence'
-        assert alltrue( time[1:] - time[:-1] > 0 ), msg	
+        assert alltrue(time[1:] - time[:-1] > 0 ), msg	
 	
 
@@ -227,6 +229,14 @@
 	    
             ####self.Q = zeros( (len(time),
-	     
-	    for i, t in enumerate(time[:]):
+	    self.T = time[:] 
+            self.index = 0 #Initial index	
+
+            self.values = {}	    	    
+	    for name in self.quantities:
+                self.values[name] = zeros( (len(self.T), 
+		                            len(interpolation_points)), 
+					    Float)	    	    
+	    
+	    for i, t in enumerate(self.T):
     
                 #Interpolate quantities at this timestep 
@@ -238,40 +248,70 @@
 					 verbose = False)
 
-                for name in self.quantities:				
-		    print t, name, interpol.interpolate(fid.variables[name][i,:])		     
+                for name in self.quantities:
+		    self.values[name][i, :] =\
+		    interpol.interpolate(fid.variables[name][i,:])   
 
         else:
             pass 
 
-        self.T = time[:]
-	
-	return
-	
-        fid = open(self.filename)        
-        lines = fid.readlines()
-        fid.close()
-        
-        N = len(lines)
-        d = self.number_of_values
-
-        T = zeros(N, Float)       #Time
-        Q = zeros((N, d), Float)  #Values
-        
-        for i, line in enumerate(lines):
-            fields = line.split(',')
-            realtime = calendar.timegm(time.strptime(fields[0], time_format))
-
-            T[i] = realtime - self.starttime
-
-            for j, value in enumerate(fields[1].split()):
-                Q[i, j] = float(value)
-
-
-
-        self.T = T     #Time
-        self.Q = Q     #Values
-        self.index = 0 #Initial index
-
-        
+    def __repr__(self):
+        return 'File function'
+
+    def __call__(self, t, x=None, y=None, point_id = None):
+        """Evaluate f(t, point_id)
+	
+	Inputs:
+	  t: time - Model time (tau = domain.starttime-self.starttime+t)  
+	            must lie within existing timesteps   
+	  point_id: index of one of the preprocessed points. 
+	            If point_id is None all preprocessed points are computed
+
+	  FIXME: point_id could also be a slice	    
+	  FIXME: One could allow arbitrary x, y coordinates 
+	  (requires computation of a new interpolator) 
+	  FIXME: What if x and y are vectors? 
+        """
+
+        from math import pi, cos, sin, sqrt
+        from Numeric import zeros, Float
+
+        #Find time tau such that
+	#
+	# domain.starttime + t == self.starttime + tau
+	
+	if self.domain is not None:
+            tau = self.domain.starttime-self.starttime+t
+	else:
+            tau = t
+	        
+        msg = 'Time interval derived from file %s (%s:%s) does not match model time: %s'\
+              %(self.filename, self.T[0], self.T[1], tau)
+        if tau < self.T[0]: raise msg
+        if tau > self.T[-1]: raise msg        
+
+        oldindex = self.index #Time index
+
+        #Find time slot
+        while tau > self.T[self.index]: self.index += 1
+        while tau < self.T[self.index]: self.index -= 1
+
+        #t is now between index and index+1
+        ratio = (tau - self.T[self.index])/\
+                (self.T[self.index+1] - self.T[self.index])
+
+		
+        #Compute interpolated values
+        q = zeros( len(self.quantities), Float)
+	
+	for i, name in enumerate(self.quantities):
+            Q = self.values[name]    
+
+            q[i] = Q[self.index, point_id] +\
+                 ratio*(Q[self.index+1, point_id] - Q[self.index, point_id]) 
+
+	    
+        #Return vector of interpolated values
+        return q
+	        
     	     
 class File_function_ASCII:
@@ -430,7 +470,5 @@
         return 'File function'
 
-        
-
-    def __call__(self, t, x=None, y=None):
+    def __call__(self, t, x=None, y=None, point_id=None):
         """Evaluate f(t,x,y)
 
@@ -438,4 +476,6 @@
         result is a vector of same length as x and y
         FIXME: Naaaa
+	
+	FIXME: Rethink semantics when x,y are vectors.
         """
 
