Index: /inundation/pyvolution/domain.py
===================================================================
--- /inundation/pyvolution/domain.py	(revision 1825)
+++ /inundation/pyvolution/domain.py	(revision 1826)
@@ -322,19 +322,47 @@
 
     def write_time(self):
+        print self.timestepping_statistics()
+
+        #Old version
+        #if self.min_timestep == self.max_timestep:
+        #    print 'Time = %.4f, delta t = %.8f, steps=%d (%d)'\
+        #          %(self.time, self.min_timestep, self.number_of_steps,
+        #            self.number_of_first_order_steps)
+        #elif self.min_timestep > self.max_timestep:
+        #    print 'Time = %.4f, steps=%d (%d)'\
+        #          %(self.time, self.number_of_steps,
+        #            self.number_of_first_order_steps)
+        #else:
+        #    print 'Time = %.4f, delta t in [%.8f, %.8f], steps=%d (%d)'\
+        #          %(self.time, self.min_timestep,
+        #            self.max_timestep, self.number_of_steps,
+        #            self.number_of_first_order_steps)        
+        
+    def timestepping_statistics(self):
+        """Return string with time stepping statistics for printing or logging
+        """
+
+        msg = ''
         if self.min_timestep == self.max_timestep:
-            print 'Time = %.4f, delta t = %.8f, steps=%d (%d)'\
-                  %(self.time, self.min_timestep, self.number_of_steps,
-                    self.number_of_first_order_steps)
+            msg += 'Time = %.4f, delta t = %.8f, steps=%d (%d)\n'\
+                   %(self.time, self.min_timestep, self.number_of_steps,
+                     self.number_of_first_order_steps)
         elif self.min_timestep > self.max_timestep:
-            print 'Time = %.4f, steps=%d (%d)'\
-                  %(self.time, self.number_of_steps,
-                    self.number_of_first_order_steps)
+            msg += 'Time = %.4f, steps=%d (%d)\n'\
+                   %(self.time, self.number_of_steps,
+                     self.number_of_first_order_steps)
         else:
-            print 'Time = %.4f, delta t in [%.8f, %.8f], steps=%d (%d)'\
-                  %(self.time, self.min_timestep,
-                    self.max_timestep, self.number_of_steps,
-                    self.number_of_first_order_steps)
-
-    def boundary_stats(self, quantities = None, tags = None):
+            msg += 'Time = %.4f, delta t in [%.8f, %.8f], steps=%d (%d)\n'\
+                   %(self.time, self.min_timestep,
+                     self.max_timestep, self.number_of_steps,
+                     self.number_of_first_order_steps)
+            
+        return msg    
+
+            
+    def write_boundary(self):
+        print self.boundary_statistics()
+        
+    def boundary_statistics(self, quantities = None, tags = None):
         """Output statistics about boundary forcing at each timestep
 
@@ -351,9 +379,8 @@
         If tags are specified only report on those, otherwise take all tags.
 
-        #FIXME: Should return text string
         """
 
         #Input checks
-        import types
+        import types, string
         
         if quantities is None:
@@ -376,63 +403,39 @@
         assert type(tags) == types.ListType, msg            
 
-
-        #Output stats    
-
-        print 'Boundary values at time %.4f:' %self.time
-
+        #Determine width of longest quantity name (for cosmetic purposes)
+        maxwidth = 0
+        for name in quantities:
+            w = len(name) 
+            if w > maxwidth:
+                maxwidth = w
+            
+        #Output stats
+        msg = 'Boundary values at time %.4f:\n' %self.time
         for tag in tags:
-            print '    Tag: %s' %tag
-
+            msg += '    %s:\n' %tag
+          
             for name in quantities:
                 q = self.quantities[name]
-
+        
                 #Find range of boundary values for tag and q 
                 maxval = minval = None
-                for i, ((vol_id, edge_id), B) in enumerate(self.boundary_objects):
+                for i, ((vol_id, edge_id), B) in\
+                        enumerate(self.boundary_objects):
                     if self.boundary[(vol_id, edge_id)] == tag:
                         v = q.boundary_values[i]
                         if minval is None or v < minval: minval = v
                         if maxval is None or v > maxval: maxval = v
-
+         
                 if minval is None or maxval is None:
-                    print 'Sorry no information about tag %s and quantity %s' %(tag, name)
+                    msg += '        Sorry no information available about' +\
+                           ' tag %s and quantity %s\n' %(tag, name)
                 else:
-                    print '        %s\t in [%12.8f, %12.8f]'\
-                          %(name, minval, maxval)
-
-
-            
-        
-        #for name in quantities:
-        #    q = self.quantities[name]
-        #
-        #    if tag is None:
-        #        #Take entire boundary
-        #        print '    Quantity %s: min = %12.8f, max = %12.8f'\
-        #              %(name, min(q.boundary_values), max(q.boundary_values))
-        #    else:
-        #        #Take only boundary associated with tag
-        #        maxval = minval = None
-        #        for i, ((vol_id, edge_id), B) in enumerate(self.boundary_objects):
-        #            if self.boundary[(vol_id, edge_id)] == tag:
-        #                v = q.boundary_values[i]
-        #                if minval is None or v < minval: minval = v
-        #                if maxval is None or v > maxval: maxval = v
-        # 
-        #        if minval is None or maxval is None:
-        #            print 'Sorry no information about tag %s' %tag
-        #        else:
-        #            print '    Quantity %s, tag %s: min = %12.8f, max = %12.8f'\
-        #                  %(name, tag, minval, maxval)
-
-
-
-
-
-
-
-
-
-
+                    msg += '        %s in [%12.8f, %12.8f]\n'\
+                           %(string.ljust(name, maxwidth), minval, maxval)
+
+
+        return msg
+
+    
     def get_name(self):
         return self.filename
Index: /inundation/pyvolution/netherlands.py
===================================================================
--- /inundation/pyvolution/netherlands.py	(revision 1825)
+++ /inundation/pyvolution/netherlands.py	(revision 1826)
@@ -169,4 +169,6 @@
 for t in domain.evolve(yieldstep = 0.02, finaltime = 15.0):
     domain.write_time()
+    domain.write_boundary()    
+    
     print domain.quantities['stage'].get_values(location='centroids',
                                                 indices=[0])
