Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py	(revision 6180)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py	(revision 6181)
@@ -38,26 +38,45 @@
 
 
+##
+# @brief Basic Domain class
+# @note Subclasses Mesh.
 class Domain(Mesh):
 
-
-    def __init__(self,
-                 source=None,
-                 triangles=None,
-                 boundary=None,
-                 conserved_quantities=None,
-                 other_quantities=None,
-                 tagged_elements=None,
-                 geo_reference=None,
-                 use_inscribed_circle=False,
-                 mesh_filename=None,
-                 use_cache=False,
-                 verbose=False,
-                 full_send_dict=None,
-                 ghost_recv_dict=None,
-                 processor=0,
-                 numproc=1,
-                 number_of_full_nodes=None,
-                 number_of_full_triangles=None):
-
+    ##
+    # @brief Generic computational Domain constructor.
+    # @param source Name of mesh file or coords of mesh vertices.
+    # @param triangles Mesh connectivity (see mesh.py for more information).
+    # @param boundary (see mesh.py for more information)
+    # @param conserved_quantities List of names of quantities to be conserved.
+    # @param other_quantities List of names of other quantities.
+    # @param tagged_elements
+    # @param geo_reference
+    # @param use_inscribed_circle
+    # @param mesh_filename
+    # @param use_cache
+    # @param verbose
+    # @param full_send_dict
+    # @param ghost_recv_dict
+    # @param processor
+    # @param numproc
+    # @param number_of_full_nodes
+    # @param number_of_full_triangles
+    def __init__(self, source=None,
+                       triangles=None,
+                       boundary=None,
+                       conserved_quantities=None,
+                       other_quantities=None,
+                       tagged_elements=None,
+                       geo_reference=None,
+                       use_inscribed_circle=False,
+                       mesh_filename=None,
+                       use_cache=False,
+                       verbose=False,
+                       full_send_dict=None,
+                       ghost_recv_dict=None,
+                       processor=0,
+                       numproc=1,
+                       number_of_full_nodes=None,
+                       number_of_full_triangles=None):
 
         """Instantiate generic computational Domain.
@@ -76,6 +95,4 @@
           tagged_elements:
           ...
-
-
         """
 
@@ -85,5 +102,4 @@
         else:
             coordinates = source
-
 
         # In case a filename has been specified, extract content
@@ -94,5 +110,4 @@
                                          use_cache=use_cache,
                                          verbose=verbose)
-
 
         # Initialise underlying mesh structure
@@ -108,6 +123,5 @@
         if verbose: print 'Initialising Domain'
 
-        # List of quantity names entering
-        # the conservation equations
+        # List of quantity names entering the conservation equations
         if conserved_quantities is None:
             self.conserved_quantities = []
@@ -121,5 +135,4 @@
             self.other_quantities = other_quantities
 
-
         # Build dictionary of Quantity instances keyed by quantity names
         self.quantities = {}
@@ -138,5 +151,5 @@
             self.full_send_dict = {}
         else:
-            self.full_send_dict  = full_send_dict
+            self.full_send_dict = full_send_dict
 
         # List of other quantity names
@@ -149,5 +162,4 @@
         self.numproc = numproc
 
-
         # Setup Communication Buffers
         if verbose: print 'Domain: Set up communication buffers (parallel)'
@@ -155,9 +167,11 @@
         for key in self.full_send_dict:
             buffer_shape = self.full_send_dict[key][0].shape[0]
-            self.full_send_dict[key].append(num.zeros((buffer_shape, self.nsys), num.Float))
+            self.full_send_dict[key].append(num.zeros((buffer_shape, self.nsys),
+                                                      num.Float))
 
         for key in self.ghost_recv_dict:
             buffer_shape = self.ghost_recv_dict[key][0].shape[0]
-            self.ghost_recv_dict[key].append(num.zeros((buffer_shape, self.nsys), num.Float))
+            self.ghost_recv_dict[key].append(num.zeros((buffer_shape, self.nsys),
+                                             num.Float))
 
         # Setup cell full flag
@@ -173,7 +187,7 @@
         # Test the assumption that all full triangles are store before
         # the ghost triangles.
-        if not num.allclose(self.tri_full_flag[:self.number_of_full_nodes],1):
-            print 'WARNING:  Not all full triangles are store before ghost triangles'
-                        
+        if not num.allclose(self.tri_full_flag[:self.number_of_full_nodes], 1):
+            print ('WARNING:  '
+                   'Not all full triangles are store before ghost triangles')
 
         # Defaults
@@ -182,9 +196,11 @@
         from anuga.config import timestepping_method
         from anuga.config import protect_against_isolated_degenerate_timesteps
-        from anuga.config import default_order        
+        from anuga.config import default_order
+
         self.beta_w = beta_w
         self.epsilon = epsilon
-        self.protect_against_isolated_degenerate_timesteps = protect_against_isolated_degenerate_timesteps
-        
+        self.protect_against_isolated_degenerate_timesteps = \
+                        protect_against_isolated_degenerate_timesteps
+
         self.set_default_order(default_order)
 
@@ -196,6 +212,5 @@
         self.set_timestepping_method(timestepping_method)
         self.set_beta(beta_w)
-        self.boundary_map = None  # Will be populated by set_boundary        
-        
+        self.boundary_map = None  # Will be populated by set_boundary
 
         # Model time
@@ -209,46 +224,47 @@
 
         self.last_walltime = walltime()
-        
+
         # Monitoring
         self.quantities_to_be_monitored = None
         self.monitor_polygon = None
-        self.monitor_time_interval = None                
+        self.monitor_time_interval = None
         self.monitor_indices = None
 
         # Checkpointing and storage
         from anuga.config import default_datadir
+
         self.datadir = default_datadir
         self.simulation_name = 'domain'
         self.checkpoint = False
 
-        # To avoid calculating the flux across each edge twice,
-        # keep an integer (boolean) array, to be used during the flux
-        # calculation
+        # To avoid calculating the flux across each edge twice, keep an integer
+        # (boolean) array, to be used during the flux calculation.
         N = len(self) # Number_of_triangles
         self.already_computed_flux = num.zeros((N, 3), num.Int)
 
         # Storage for maximal speeds computed for each triangle by
-        # compute_fluxes
+        # compute_fluxes.
         # This is used for diagnostics only (reset at every yieldstep)
         self.max_speed = num.zeros(N, num.Float)
 
         if mesh_filename is not None:
-            # If the mesh file passed any quantity values
-            # , initialise with these values.
+            # If the mesh file passed any quantity values,
+            # initialise with these values.
             if verbose: print 'Domain: Initialising quantity values'
             self.set_quantity_vertices_dict(vertex_quantity_dict)
 
-
         if verbose: print 'Domain: Done'
 
-
-
-
-
-
-    #---------------------------    
-    # Public interface to Domain
-    #---------------------------        
-    def get_conserved_quantities(self, vol_id, vertex=None, edge=None):
+    ##
+    # @brief Get conserved quantities for a volume.
+    # @param vol_id ID of the volume we want the conserved quantities for.
+    # @param vertex If specified, use as index for edge values.
+    # @param edge If specified, use as index for edge values.
+    # @return Vector of conserved quantities.
+    # @note If neither 'vertex' or 'edge' specified, use centroid values.
+    # @note If both 'vertex' and 'edge' specified, raise exception.
+    def get_conserved_quantities(self, vol_id,
+                                       vertex=None,
+                                       edge=None):
         """Get conserved quantities at volume vol_id
 
@@ -259,5 +275,4 @@
 
         Return value: Vector of length == number_of_conserved quantities
-
         """
 
@@ -265,5 +280,5 @@
             msg = 'Values for both vertex and edge was specified.'
             msg += 'Only one (or none) is allowed.'
-            raise msg
+            raise Exception, msg
 
         q = num.zeros(len(self.conserved_quantities), num.Float)
@@ -280,4 +295,7 @@
         return q
 
+    ##
+    # @brief Set the relative model time.
+    # @param time The new model time (seconds).
     def set_time(self, time=0.0):
         """Set the model time (seconds)"""
@@ -287,4 +305,7 @@
         self.time = time
 
+    ##
+    # @brief Get the model time.
+    # @return The absolute model time (seconds).
     def get_time(self):
         """Get the absolute model time (seconds)"""
@@ -292,7 +313,9 @@
         return self.time + self.starttime
 
-    def set_beta(self,beta):
-        """Set default beta for limiting
-        """
+    ##
+    # @brief Set the default beta for limiting.
+    # @param beta The new beta value.
+    def set_beta(self, beta):
+        """Set default beta for limiting"""
 
         self.beta = beta
@@ -302,15 +325,20 @@
             Q.set_beta(beta)
 
+    ##
+    # @brief Get the beta value used for limiting.
+    # @return The beta value used for limiting.
     def get_beta(self):
-        """Get default beta for limiting
-        """
+        """Get default beta for limiting"""
 
         return self.beta
 
+    ##
+    # @brief Set default (spatial) order.
+    # @param n The new spatial order value.
+    # @note If 'n' is not 1 or 2, raise exception.
     def set_default_order(self, n):
-        """Set default (spatial) order to either 1 or 2
-        """
-
-        msg = 'Default order must be either 1 or 2. I got %s' %n
+        """Set default (spatial) order to either 1 or 2"""
+
+        msg = 'Default order must be either 1 or 2. I got %s' % n
         assert n in [1,2], msg
 
@@ -327,18 +355,18 @@
             #self.set_timestepping_method('rk2')
             #self.set_all_limiters(beta_rk2)
-        
-
+
+    ##
+    # @brief Set values of named quantities.
+    # @param quantity_dict Dictionary containing name/value pairs.
     def set_quantity_vertices_dict(self, quantity_dict):
         """Set values for named quantities.
-        The index is the quantity
-
-        name: Name of quantity
-        X: Compatible list, Numeric array, const or function (see below)
-
-        The values will be stored in elements following their
-        internal ordering.
-
-        """
-        
+        Supplied dictionary contains name/value pairs:
+
+        name:  Name of quantity
+        value: Compatible list, Numeric array, const or function (see below)
+
+        The values will be stored in elements following their internal ordering.
+        """
+
         # FIXME: Could we name this a bit more intuitively
         # E.g. set_quantities_from_dictionary
@@ -346,8 +374,13 @@
             self.set_quantity(key, quantity_dict[key], location='vertices')
 
-
-    def set_quantity(self, name, *args, **kwargs):
+    ##
+    # @brief Set value(s) for a named quantity.
+    # @param name Name of quantity to be updated.
+    # @param args Positional args.
+    # @param kwargs Keyword args.
+    # @note If 'kwargs' dict has 'expression' key, evaluate expression.
+    def set_quantity(self, name,
+                           *args, **kwargs):
         """Set values for named quantity
-
 
         One keyword argument is documented here:
@@ -370,36 +403,55 @@
         # Assign values
         self.quantities[name].set_values(*args, **kwargs)
-        
-    def add_quantity(self, name, *args, **kwargs):
+
+    ##
+    # @brief Add to a named quantity value.
+    # @param name Name of quantity to be added to.
+    # @param args Positional args.
+    # @param kwargs Keyword args.
+    # @note If 'kwargs' dict has 'expression' key, evaluate expression.
+    def add_quantity(self, name,
+                           *args, **kwargs):
         """Add values to a named quantity
-        
+
         E.g add_quantity('elevation', X)
-        
+
         Option are the same as in set_quantity.
         """
-            
+
         # Do the expression stuff
         if kwargs.has_key('expression'):
             expression = kwargs['expression']
             Q2 = self.create_quantity_from_expression(expression)
-        else:    
+        else:
             # Create new temporary quantity
             Q2 = Quantity(self)
-        
+
             # Assign specified values to temporary quantity
             Q2.set_values(*args, **kwargs)
-            
+
         # Add temporary quantity to named quantity
         Q1 = self.get_quantity(name)
         self.set_quantity(name, Q1 + Q2)
-        
-
+
+    ##
+    # @brief Get list of quantity names for the Domain.
+    # @return List of quantity names.
     def get_quantity_names(self):
         """Get a list of all the quantity names that this domain is aware of.
         Any value in the result should be a valid input to get_quantity.
         """
+
         return self.quantities.keys()
 
-    def get_quantity(self, name, location='vertices', indices = None):
+    ##
+    # @brief Get a quantity object.
+    # @param name Name of the quantity value.
+    # @param location ??
+    # @param indices ??
+    # @return The quantity value object.
+    # @note 'location' and 'indices' are unused.
+    def get_quantity(self, name,
+                           location='vertices',
+                           indices = None):
         """Get pointer to quantity object.
 
@@ -413,8 +465,11 @@
         return self.quantities[name] #.get_values( location, indices = indices)
 
-
-
+    ##
+    # @brief Create a quantity value from an expression.
+    # @param expression The expression (string) to be evaluated.
+    # @return The expression value, evaluated from this Domain's quantities.
+    # @note Valid expression operators are as defined in class Quantity.
     def create_quantity_from_expression(self, expression):
-        """Create new quantity from other quantities using arbitrary expression
+        """Create new quantity from other quantities using arbitrary expression.
 
         Combine existing quantities in domain using expression and return
@@ -426,19 +481,19 @@
 
         Examples creating derived quantities:
-
             Depth = domain.create_quantity_from_expression('stage-elevation')
-
-            exp = '(xmomentum*xmomentum + ymomentum*ymomentum)**0.5'        
+            exp = '(xmomentum*xmomentum + ymomentum*ymomentum)**0.5'
             Absolute_momentum = domain.create_quantity_from_expression(exp)
-
         """
 
         from anuga.abstract_2d_finite_volumes.util import\
              apply_expression_to_dictionary
-        
+
         return apply_expression_to_dictionary(expression, self.quantities)
 
-
-
+    ##
+    # @brief Associate boundary objects with tagged boundary segments.
+    # @param boundary_map A dict of boundary objects keyed by symbolic tags to
+    #                     matched against tags in the internal dictionary
+    #                     self.boundary.
     def set_boundary(self, boundary_map):
         """Associate boundary objects with tagged boundary segments.
@@ -461,5 +516,4 @@
         self.boundary_objects:  ((vol_id, edge_id), boundary_object)
 
-
         Pre-condition:
           self.boundary has been built.
@@ -471,24 +525,20 @@
         exception is raised.
         However, if a tag is not used to the domain, no error is thrown.
-        FIXME: This would lead to implementation of a
-        default boundary condition
+        FIXME: This would lead to implementation of a default boundary condition
 
         Note: If a segment is listed in the boundary dictionary and if it is
-        not None, it *will* become a boundary -
-        even if there is a neighbouring triangle.
-        This would be the case for internal boundaries
+        not None, it *will* become a boundary - even if there is a neighbouring
+        triangle.  This would be the case for internal boundaries.
 
         Boundary objects that are None will be skipped.
 
-        If a boundary_map has already been set
-        (i.e. set_boundary has been called before), the old boundary map
-        will be updated with new values. The new map need not define all
-        boundary tags, and can thus change only those that are needed.
+        If a boundary_map has already been set (i.e. set_boundary has been
+        called before), the old boundary map will be updated with new values.
+        The new map need not define all boundary tags, and can thus change only
+        those that are needed.
 
         FIXME: If set_boundary is called multiple times and if Boundary
         object is changed into None, the neighbour structure will not be
         restored!!!
-
-
         """
 
@@ -496,13 +546,11 @@
             # This the first call to set_boundary. Store
             # map for later updates and for use with boundary_stats.
-            self.boundary_map = boundary_map 
-        else:    
+            self.boundary_map = boundary_map
+        else:
             # This is a modification of an already existing map
             # Update map an proceed normally
-
             for key in boundary_map.keys():
                 self.boundary_map[key] = boundary_map[key]
-                
-        
+
         # FIXME (Ole): Try to remove the sorting and fix test_mesh.py
         x = self.boundary.keys()
@@ -511,7 +559,7 @@
         # Loop through edges that lie on the boundary and associate them with
         # callable boundary objects depending on their tags
-        self.boundary_objects = []        
+        self.boundary_objects = []
         for k, (vol_id, edge_id) in enumerate(x):
-            tag = self.boundary[ (vol_id, edge_id) ]
+            tag = self.boundary[(vol_id, edge_id)]
 
             if self.boundary_map.has_key(tag):
@@ -519,12 +567,10 @@
 
                 if B is not None:
-                    self.boundary_objects.append( ((vol_id, edge_id), B) )
+                    self.boundary_objects.append(((vol_id, edge_id), B))
                     self.neighbours[vol_id, edge_id] = \
-                                            -len(self.boundary_objects)
+                                        -len(self.boundary_objects)
                 else:
                     pass
                     #FIXME: Check and perhaps fix neighbour structure
-
-
             else:
                 msg = 'ERROR (domain.py): Tag "%s" has not been ' %tag
@@ -533,16 +579,18 @@
                 msg += 'in set_boundary.\n'
                 msg += 'The tags are: %s' %self.get_boundary_tags()
-                raise msg
-
-
+                raise Exception, msg
+
+    ##
+    # @brief Set quantities based on a regional tag.
+    # @param args 
+    # @param kwargs 
     def set_region(self, *args, **kwargs):
-        """
-        This method is used to set quantities based on a regional tag.
-        
+        """Set quantities based on a regional tag.
+
         It is most often called with the following parameters;
         (self, tag, quantity, X, location='vertices')
-        tag: the name of the regional tag used to specify the region 
+        tag:      the name of the regional tag used to specify the region
         quantity: Name of quantity to change
-        X: const or function - how the quantity is changed
+        X:        const or function - how the quantity is changed
         location: Where values are to be stored.
             Permissible options are: vertices, centroid and unique vertices
@@ -558,24 +606,31 @@
             func = region_set_region(*args, **kwargs)
             self._set_region(func)
-            
-        
+
+    ##
+    # @brief ??
+    # @param functions A list or tuple of ??
     def _set_region(self, functions):
+        # coerce to an iterable (list or tuple)
+        if type(functions) not in [types.ListType, types.TupleType]:
+            functions = [functions]
+
         # The order of functions in the list is used.
-        if type(functions) not in [types.ListType,types.TupleType]:
-            functions = [functions]
         for function in functions:
             for tag in self.tagged_elements.keys():
                 function(tag, self.tagged_elements[tag], self)
 
-
-
-
+    ##
+    # @brief Specify the quantities which will be monitored for extrema.
+    # @param q Single or list of quantity names to monitor.
+    # @param polygon If specified, monitor only triangles inside polygon.
+    # @param time_interval If specified, monitor only timesteps inside interval.
+    # @note If 'q' is None, do no monitoring.
     def set_quantities_to_be_monitored(self, q,
-                                       polygon=None,
-                                       time_interval=None):
+                                             polygon=None,
+                                             time_interval=None):
         """Specify which quantities will be monitored for extrema.
 
         q must be either:
-          - the name of a quantity or a derived quantity such as 'stage-elevation'
+          - the name of a quantity or derived quantity such as 'stage-elevation'
           - a list of quantity names
           - None
@@ -583,16 +638,16 @@
         In the two first cases, the named quantities will be monitored at
         each internal timestep
-        
+
         If q is None, monitoring will be switched off altogether.
 
-        polygon (if specified) will restrict monitoring to triangles inside polygon.
+        polygon (if specified) will only monitor triangles inside polygon.
         If omitted all triangles will be included.
 
-        time_interval (if specified) will restrict monitoring to time steps in
+        time_interval, if specified, will restrict monitoring to time steps in
         that interval. If omitted all timesteps will be included.
         """
 
         from anuga.abstract_2d_finite_volumes.util import\
-             apply_expression_to_dictionary        
+             apply_expression_to_dictionary
 
         if q is None:
@@ -600,16 +655,16 @@
             self.monitor_polygon = None
             self.monitor_time_interval = None
-            self.monitor_indices = None            
+            self.monitor_indices = None
             return
 
+        # coerce 'q' to a list if it's a string
         if isinstance(q, basestring):
-            q = [q] # Turn argument into a list
-
-        # Check correcness and initialise
+            q = [q]
+
+        # Check correctness and initialise
         self.quantities_to_be_monitored = {}
         for quantity_name in q:
-            msg = 'Quantity %s is not a valid conserved quantity'\
-                  %quantity_name
-            
+            msg = 'Quantity %s is not a valid conserved quantity' \
+                      % quantity_name
 
             if not quantity_name in self.quantities:
@@ -622,28 +677,25 @@
                           'min_location': None, # Argmin (x, y)
                           'max_location': None, # Argmax (x, y)
-                          'min_time': None,     # Argmin (t) 
+                          'min_time': None,     # Argmin (t)
                           'max_time': None}     # Argmax (t)
-            
+
             self.quantities_to_be_monitored[quantity_name] = info_block
-
-            
 
         if polygon is not None:
             # Check input
             if isinstance(polygon, basestring):
-
                 # Check if multiple quantities were accidentally
                 # given as separate argument rather than a list.
-                msg = 'Multiple quantities must be specified in a list. '
-                msg += 'Not as multiple arguments. '
-                msg += 'I got "%s" as a second argument' %polygon
-                
+                msg = ('Multiple quantities must be specified in a list. '
+                       'Not as multiple arguments. '
+                       'I got "%s" as a second argument') % polygon
+
                 if polygon in self.quantities:
                     raise Exception, msg
-                
+
                 try:
                     apply_expression_to_dictionary(polygon, self.quantities)
                 except:
-                    # At least polygon wasn't an expression involving quantitites
+                    # At least polygon wasn't expression involving quantitites
                     pass
                 else:
@@ -651,25 +703,21 @@
 
                 # In any case, we don't allow polygon to be a string
-                msg = 'argument "polygon" must not be a string: '
-                msg += 'I got polygon=\'%s\' ' %polygon
+                msg = ('argument "polygon" must not be a string: '
+                       'I got polygon="%s"') % polygon
                 raise Exception, msg
-
 
             # Get indices for centroids that are inside polygon
             points = self.get_centroid_coordinates(absolute=True)
             self.monitor_indices = inside_polygon(points, polygon)
-            
 
         if time_interval is not None:
             assert len(time_interval) == 2
 
-        
         self.monitor_polygon = polygon
-        self.monitor_time_interval = time_interval        
-        
-
-    #--------------------------    
-    # Miscellaneous diagnostics
-    #--------------------------        
+        self.monitor_time_interval = time_interval
+
+    ##
+    # @brief Check Domain integrity.
+    # @note Raises an exception if integrity breached.
     def check_integrity(self):
         Mesh.check_integrity(self)
@@ -680,14 +728,19 @@
 
         ##assert hasattr(self, 'boundary_objects')
-            
-
+
+    ##
+    # @brief Print timestep stats to stdout.
+    # @param track_speeds If True, print smallest track speed.
     def write_time(self, track_speeds=False):
         print self.timestepping_statistics(track_speeds)
 
-
-    def timestepping_statistics(self,
-                                track_speeds=False,
-                                triangle_id=None):
-        """Return string with time stepping statistics for printing or logging
+    ##
+    # @brief Get timestepping stats string.
+    # @param track_speeds If True, report location of smallest timestep.
+    # @param triangle_id If specified, use specific triangle.
+    # @return A string containing timestep stats.
+    def timestepping_statistics(self, track_speeds=False,
+                                      triangle_id=None):
+        """Return string with time stepping statistics
 
         Optional boolean keyword track_speeds decides whether to report
@@ -696,13 +749,11 @@
 
         Optional keyword triangle_id can be used to specify a particular
-        triangle rather than the one with the largest speed. 
+        triangle rather than the one with the largest speed.
         """
 
         from anuga.utilities.numerical_tools import histogram, create_bins
-        
-
-        # qwidth determines the text field used for quantities
+
+        # qwidth determines the the width of the text field used for quantities
         qwidth = self.qwidth = 12
-
 
         msg = ''
@@ -723,20 +774,19 @@
         model_time = self.get_time()
         if self.min_timestep == self.max_timestep:
-            msg += 'Time = %.4f, delta t = %.8f, steps=%d'\
-                   %(model_time, self.min_timestep, self.number_of_steps)
+            msg += 'Time = %.4f, delta t = %.8f, steps=%d' \
+                       % (model_time, self.min_timestep, self.number_of_steps)
         elif self.min_timestep > self.max_timestep:
-            msg += 'Time = %.4f, steps=%d'\
-                   %(model_time, self.number_of_steps)
+            msg += 'Time = %.4f, steps=%d' \
+                       % (model_time, self.number_of_steps)
         else:
-            msg += 'Time = %.4f, delta t in [%.8f, %.8f], steps=%d'\
-                   %(model_time, self.min_timestep,
-                     self.max_timestep, self.number_of_steps)
-                                          
-        msg += ' (%ds)' %(walltime() - self.last_walltime)    
-        self.last_walltime = walltime()            
-        
+            msg += 'Time = %.4f, delta t in [%.8f, %.8f], steps=%d' \
+                       % (model_time, self.min_timestep,
+                          self.max_timestep, self.number_of_steps)
+
+        msg += ' (%ds)' % (walltime() - self.last_walltime)
+        self.last_walltime = walltime()
+
         if track_speeds is True:
             msg += '\n'
-
 
             # Setup 10 bins for speed histogram
@@ -745,6 +795,6 @@
 
             msg += '------------------------------------------------\n'
-            msg += '  Speeds in [%f, %f]\n' %(min(self.max_speed),
-                                              max(self.max_speed))
+            msg += '  Speeds in [%f, %f]\n' % (min(self.max_speed),
+                                               max(self.max_speed))
             msg += '  Histogram:\n'
 
@@ -753,12 +803,11 @@
                 lo = hi
                 if i+1 < len(bins):
-                    # Open upper interval                
+                    # Open upper interval
                     hi = bins[i+1]
-                    msg += '    [%f, %f[: %d\n' %(lo, hi, count)                
+                    msg += '    [%f, %f[: %d\n' % (lo, hi, count)
                 else:
                     # Closed upper interval
                     hi = max(self.max_speed)
-                    msg += '    [%f, %f]: %d\n' %(lo, hi, count)
-
+                    msg += '    [%f, %f]: %d\n' % (lo, hi, count)
 
             N = len(self.max_speed)
@@ -770,45 +819,40 @@
                 k = 0
                 lower = min(speed)
-                for i, a in enumerate(speed):        
+                for i, a in enumerate(speed):
                     if i % (N/10) == 0 and i != 0:
                         # For every 10% of the sorted speeds
-                        msg += '    %d speeds in [%f, %f]\n' %(i-k, lower, a)
+                        msg += '    %d speeds in [%f, %f]\n' % (i-k, lower, a)
                         lower = a
                         k = i
-                        
+
                 msg += '    %d speeds in [%f, %f]\n'\
-                       %(N-k, lower, max(speed))                    
-                
-                      
-
-                
-            
+                           % (N-k, lower, max(speed))
+
             # Find index of largest computed flux speed
             if triangle_id is None:
                 k = self.k = num.argmax(self.max_speed)
             else:
-                errmsg = 'Triangle_id %d does not exist in mesh: %s' %(triangle_id,
-                                                                    str(self))
+                errmsg = 'Triangle_id %d does not exist in mesh: %s' \
+                             % (triangle_id, str(self))
                 assert 0 <= triangle_id < len(self), errmsg
                 k = self.k = triangle_id
-            
 
             x, y = self.get_centroid_coordinates()[k]
             radius = self.get_radii()[k]
-            area = self.get_areas()[k]      
-            max_speed = self.max_speed[k]           
-
-            msg += '  Triangle #%d with centroid (%.4f, %.4f), ' %(k, x, y)
-            msg += 'area = %.4f and radius = %.4f ' %(area, radius)
+            area = self.get_areas()[k]
+            max_speed = self.max_speed[k]
+
+            msg += '  Triangle #%d with centroid (%.4f, %.4f), ' % (k, x, y)
+            msg += 'area = %.4f and radius = %.4f ' % (area, radius)
             if triangle_id is None:
-                msg += 'had the largest computed speed: %.6f m/s ' %(max_speed)
+                msg += 'had the largest computed speed: %.6f m/s ' % (max_speed)
             else:
-                msg += 'had computed speed: %.6f m/s ' %(max_speed)
-                
+                msg += 'had computed speed: %.6f m/s ' % (max_speed)
+
             if max_speed > 0.0:
-                msg += '(timestep=%.6f)\n' %(radius/max_speed)
+                msg += '(timestep=%.6f)\n' % (radius/max_speed)
             else:
-                msg += '(timestep=%.6f)\n' %(0)     
-            
+                msg += '(timestep=%.6f)\n' % (0)
+
             # Report all quantity values at vertices, edges and centroid
             msg += '    Quantity'
@@ -819,14 +863,14 @@
                 V = q.get_values(location='vertices', indices=[k])[0]
                 E = q.get_values(location='edges', indices=[k])[0]
-                C = q.get_values(location='centroids', indices=[k])                 
-                
-                s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
-                     %(name.ljust(qwidth), V[0], V[1], V[2])
-
-                s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
-                     %(name.ljust(qwidth), E[0], E[1], E[2])
-
-                s += '    %s: centroid_value = %.4f\n'\
-                     %(name.ljust(qwidth), C[0])                                
+                C = q.get_values(location='centroids', indices=[k])
+
+                s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n' \
+                         % (name.ljust(qwidth), V[0], V[1], V[2])
+
+                s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n' \
+                         % (name.ljust(qwidth), E[0], E[1], E[2])
+
+                s += '    %s: centroid_value = %.4f\n' \
+                         % (name.ljust(qwidth), C[0])
 
                 msg += s
@@ -834,11 +878,18 @@
         return msg
 
-
-    def write_boundary_statistics(self, quantities = None, tags = None):
+    ##
+    # @brief Print boundary forcing stats at each timestep to stdout.
+    # @param quantities A name or list of names of quantities to report.
+    # @param tags A name or list of names of tags to report.
+    def write_boundary_statistics(self, quantities=None, tags=None):
         print self.boundary_statistics(quantities, tags)
 
-    def boundary_statistics(self, quantities = None, tags = None):
+    # @brief Get a string containing boundary forcing stats at each timestep.
+    # @param quantities A name or list of names of quantities to report.
+    # @param tags A name or list of names of tags to report.
+    # @note If 'quantities' is None, report all.  Same for 'tags'.
+    def boundary_statistics(self, quantities=None,
+                                  tags=None):
         """Output statistics about boundary forcing at each timestep
-
 
         Input:
@@ -847,5 +898,4 @@
           tags:       either None, a string or a list of strings naming the
                       tags to be reported
-
 
         Example output:
@@ -856,14 +906,12 @@
         Tag 'ocean'
 
-
         If quantities are specified only report on those. Otherwise take all
         conserved quantities.
         If tags are specified only report on those, otherwise take all tags.
-
-        """
+        """
+
+        import types, string
 
         # Input checks
-        import types, string
-
         if quantities is None:
             quantities = self.conserved_quantities
@@ -871,8 +919,7 @@
             quantities = [quantities] #Turn it into a list
 
-        msg = 'Keyword argument quantities must be either None, '
-        msg += 'string or list. I got %s' %str(quantities)
+        msg = ('Keyword argument quantities must be either None, '
+               'string or list. I got %s') % str(quantities)
         assert type(quantities) == types.ListType, msg
-
 
         if tags is None:
@@ -881,6 +928,6 @@
             tags = [tags] #Turn it into a list
 
-        msg = 'Keyword argument tags must be either None, '
-        msg += 'string or list. I got %s' %str(tags)
+        msg = ('Keyword argument tags must be either None, '
+               'string or list. I got %s') % str(tags)
         assert type(tags) == types.ListType, msg
 
@@ -893,7 +940,7 @@
 
         # Output statistics
-        msg = 'Boundary values at time %.4f:\n' %self.get_time()
+        msg = 'Boundary values at time %.4f:\n' % self.get_time()
         for tag in tags:
-            msg += '    %s:\n' %tag
+            msg += '    %s:\n' % tag
 
             for name in quantities:
@@ -902,6 +949,5 @@
                 # 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]
@@ -910,14 +956,14 @@
 
                 if minval is None or maxval is None:
-                    msg += '        Sorry no information available about' +\
-                           ' tag %s and quantity %s\n' %(tag, name)
+                    msg += ('        Sorry no information available about'
+                            ' tag %s and quantity %s\n') % (tag, name)
                 else:
-                    msg += '        %s in [%12.8f, %12.8f]\n'\
-                           %(string.ljust(name, maxwidth), minval, maxval)
-
+                    msg += '        %s in [%12.8f, %12.8f]\n' \
+                               % (string.ljust(name, maxwidth), minval, maxval)
 
         return msg
 
-
+    ##
+    # @brief Update extrema if requested by set_quantities_to_be_monitored.
     def update_extrema(self):
         """Update extrema if requested by set_quantities_to_be_monitored.
@@ -929,5 +975,5 @@
         # Define a tolerance for extremum computations
         from anuga.config import single_precision as epsilon
-        
+
         if self.quantities_to_be_monitored is None:
             return
@@ -954,5 +1000,5 @@
             # of the epsilon)
             maxval = Q.get_maximum_value(self.monitor_indices)
-            if info_block['max'] is None or\
+            if info_block['max'] is None or \
                    maxval > info_block['max'] + epsilon:
                 info_block['max'] = maxval
@@ -961,17 +1007,18 @@
                 info_block['max_time'] = self.time
 
-
             # Update minimum
             minval = Q.get_minimum_value(self.monitor_indices)
-            if info_block['min'] is None or\
+            if info_block['min'] is None or \
                    minval < info_block['min'] - epsilon:
-                info_block['min'] = minval                
+                info_block['min'] = minval
                 minloc = Q.get_minimum_location()
                 info_block['min_location'] = minloc
-                info_block['min_time'] = self.time                
-        
-
-
-    def quantity_statistics(self, precision = '%.4f'):
+                info_block['min_time'] = self.time
+
+    ##
+    # @brief Return string with statistics about quantities
+    # @param precision A format string to use for float values.
+    # @return The stats string.
+    def quantity_statistics(self, precision='%.4f'):
         """Return string with statistics about quantities for
         printing or logging
@@ -980,5 +1027,4 @@
 
            set_quantities_to_be_monitored
-           
         """
 
@@ -986,8 +1032,8 @@
 
         # Output statistics
-        msg = 'Monitored quantities at time %.4f:\n' %self.get_time()
+        msg = 'Monitored quantities at time %.4f:\n' % self.get_time()
         if self.monitor_polygon is not None:
             p_str = str(self.monitor_polygon)
-            msg += '- Restricted by polygon: %s' %p_str[:maxlen]
+            msg += '- Restricted by polygon: %s' % p_str[:maxlen]
             if len(p_str) >= maxlen:
                 msg += '...\n'
@@ -995,83 +1041,107 @@
                 msg += '\n'
 
-
         if self.monitor_time_interval is not None:
-            msg += '- Restricted by time interval: %s\n'\
-                   %str(self.monitor_time_interval)
+            msg += '- Restricted by time interval: %s\n' \
+                       % str(self.monitor_time_interval)
             time_interval_start = self.monitor_time_interval[0]
         else:
             time_interval_start = 0.0
 
-            
         for quantity_name, info in self.quantities_to_be_monitored.items():
-            msg += '    %s:\n' %quantity_name
-
-            msg += '      values since time = %.2f in [%s, %s]\n'\
-                   %(time_interval_start,
-                     get_textual_float(info['min'], precision),
-                     get_textual_float(info['max'], precision))       
-                     
-            msg += '      minimum attained at time = %s, location = %s\n'\
-                   %(get_textual_float(info['min_time'], precision),
-                     get_textual_float(info['min_location'], precision))
-            
-
-            msg += '      maximum attained at time = %s, location = %s\n'\
-                   %(get_textual_float(info['max_time'], precision),
-                     get_textual_float(info['max_location'], precision))
-
+            msg += '    %s:\n' % quantity_name
+
+            msg += '      values since time = %.2f in [%s, %s]\n' \
+                       % (time_interval_start,
+                          get_textual_float(info['min'], precision),
+                          get_textual_float(info['max'], precision))
+
+            msg += '      minimum attained at time = %s, location = %s\n' \
+                       % (get_textual_float(info['min_time'], precision),
+                          get_textual_float(info['min_location'], precision))
+
+            msg += '      maximum attained at time = %s, location = %s\n' \
+                       % (get_textual_float(info['max_time'], precision),
+                          get_textual_float(info['max_location'], precision))
 
         return msg
 
+    ##
+    # @brief Get the timestep method.
+    # @return The timestep method. One of 'euler', 'rk2' or 'rk3'.
     def get_timestepping_method(self):
         return self.timestepping_method
 
-    def set_timestepping_method(self,timestepping_method):
-        
+    ##
+    # @brief Set the tmestep method to be used.
+    # @param timestepping_method One of 'euler', 'rk2' or 'rk3'.
+    # @note Raises exception of method not known.
+    def set_timestepping_method(self, timestepping_method):
         if timestepping_method in ['euler', 'rk2', 'rk3']:
             self.timestepping_method = timestepping_method
             return
 
-        msg = '%s is an incorrect timestepping type'% timestepping_method
+        msg = '%s is an incorrect timestepping type' % timestepping_method
         raise Exception, msg
 
+    ##
+    # @brief Get the Domain simulation name.
+    # @return The simulation name string.
     def get_name(self):
         return self.simulation_name
 
+    ##
+    # @brief Set the simulation name.
+    # @param name The name of the simulation.
+    # @note The simulation name is also used for the output .sww file.
     def set_name(self, name):
         """Assign a name to this simulation.
         This will be used to identify the output sww file.
-
-        """
+        """
+
+        # remove any '.sww' end
         if name.endswith('.sww'):
             name = name[:-4]
-            
+
         self.simulation_name = name
 
+    ##
+    # @brief Get data directory path.
+    # @return The data directory path string.
     def get_datadir(self):
         return self.datadir
 
+    ##
+    # @brief Set data directory path.
+    # @param name The data directory path string.
     def set_datadir(self, name):
         self.datadir = name
 
+    ##
+    # @brief Get the start time value.
+    # @return The start time value (float).
     def get_starttime(self):
         return self.starttime
 
+    ##
+    # @brief Set the start time value.
+    # @param time The start time value.
     def set_starttime(self, time):
-        self.starttime = float(time)        
-
-
-
-    #--------------------------
-    # Main components of evolve
-    #--------------------------    
-
-    def evolve(self,
-               yieldstep = None,
-               finaltime = None,
-               duration = None,
-               skip_initial_step = False):
+        self.starttime = float(time)
+
+#--------------------------
+# Main components of evolve
+#--------------------------
+
+    ##
+    # @brief Evolve the model through time.
+    # @param yieldstep Interval between yields where results are stored, etc.
+    # @param finaltime Time where simulation should end.
+    # @param duration Duration of simulation.
+    # @param skip_initial_step If True, skip the first yield step.
+    def evolve(self, yieldstep=None,
+                     finaltime=None,
+                     duration=None,
+                     skip_initial_step=False):
         """Evolve model through time starting from self.starttime.
-
 
         yieldstep: Interval between yields where results are stored,
@@ -1088,10 +1158,8 @@
         If both duration and finaltime are given an exception is thrown.
 
-
         skip_initial_step: Boolean flag that decides whether the first
         yield step is skipped or not. This is useful for example to avoid
         duplicate steps when multiple evolve processes are dove tailed.
 
-
         Evolve is implemented as a generator and is to be called as such, e.g.
 
@@ -1099,7 +1167,5 @@
             <Do something with domain and t>
 
-
         All times are given in seconds
-
         """
 
@@ -1107,11 +1173,10 @@
 
         # FIXME: Maybe lump into a larger check prior to evolving
-        msg = 'Boundary tags must be bound to boundary objects before '
-        msg += 'evolving system, '
-        msg += 'e.g. using the method set_boundary.\n'
-        msg += 'This system has the boundary tags %s '\
-               %self.get_boundary_tags()
+        msg = ('Boundary tags must be bound to boundary objects before '
+               'evolving system, '
+               'e.g. using the method set_boundary.\n'
+               'This system has the boundary tags %s ') \
+                   % self.get_boundary_tags()
         assert hasattr(self, 'boundary_objects'), msg
-
 
         if yieldstep is None:
@@ -1122,9 +1187,8 @@
         self._order_ = self.default_order
 
-
         if finaltime is not None and duration is not None:
             # print 'F', finaltime, duration
             msg = 'Only one of finaltime and duration may be specified'
-            raise msg
+            raise Exception, msg
         else:
             if finaltime is not None:
@@ -1133,6 +1197,4 @@
                 self.finaltime = self.starttime + float(duration)
 
-
-
         N = len(self) # Number of triangles
         self.yieldtime = 0.0 # Track time between 'yields'
@@ -1144,5 +1206,4 @@
         self.number_of_first_order_steps = 0
 
-
         # Update ghosts
         self.update_ghosts()
@@ -1153,5 +1214,5 @@
         # Update extrema if necessary (for reporting)
         self.update_extrema()
-        
+
         # Initial update boundary values
         self.update_boundary()
@@ -1165,18 +1226,16 @@
 
         while True:
-
             # Evolve One Step, using appropriate timestepping method
             if self.get_timestepping_method() == 'euler':
-                self.evolve_one_euler_step(yieldstep,finaltime)
-                
+                self.evolve_one_euler_step(yieldstep, finaltime)
+
             elif self.get_timestepping_method() == 'rk2':
-                self.evolve_one_rk2_step(yieldstep,finaltime)
+                self.evolve_one_rk2_step(yieldstep, finaltime)
 
             elif self.get_timestepping_method() == 'rk3':
-                self.evolve_one_rk3_step(yieldstep,finaltime)               
-            
+                self.evolve_one_rk3_step(yieldstep, finaltime)
+
             # Update extrema if necessary (for reporting)
-            self.update_extrema()            
-
+            self.update_extrema()
 
             self.yieldtime += self.timestep
@@ -1187,13 +1246,11 @@
             # Yield results
             if finaltime is not None and self.time >= finaltime-epsilon:
-
                 if self.time > finaltime:
                     # FIXME (Ole, 30 April 2006): Do we need this check?
                     # Probably not (Ole, 18 September 2008). Now changed to
                     # Exception
-                    msg = 'WARNING (domain.py): time overshot finaltime. '
-                    msg += 'Contact Ole.Nielsen@ga.gov.au'
+                    msg = ('WARNING (domain.py): time overshot finaltime. '
+                           'Contact Ole.Nielsen@ga.gov.au')
                     raise Exception, msg
-                    
 
                 # Yield final time and stop
@@ -1202,8 +1259,6 @@
                 break
 
-
             if self.yieldtime >= yieldstep:
                 # Yield (intermediate) time and allow inspection of domain
-
                 if self.checkpoint is True:
                     self.store_checkpoint()
@@ -1221,7 +1276,10 @@
                 self.max_speed = num.zeros(N, num.Float)
 
+    ##
+    # @brief 'Euler' time step method.
+    # @param yieldstep The reporting time step.
+    # @param finaltime The simulation final time.
     def evolve_one_euler_step(self, yieldstep, finaltime):
-        """
-        One Euler Time Step
+        """One Euler Time Step
         Q^{n+1} = E(h) Q^n
         """
@@ -1248,16 +1306,15 @@
         self.update_boundary()
 
-
-        
-
-
+    ##
+    # @brief 'rk2' time step method.
+    # @param yieldstep The reporting time step.
+    # @param finaltime The simulation final time.
     def evolve_one_rk2_step(self, yieldstep, finaltime):
-        """
-        One 2nd order RK timestep
+        """One 2nd order RK timestep
         Q^{n+1} = 0.5 Q^n + 0.5 E(h)^2 Q^n
         """
 
         # Save initial initial conserved quantities values
-        self.backup_conserved_quantities()            
+        self.backup_conserved_quantities()
 
         #--------------------------------------
@@ -1289,5 +1346,5 @@
         # Second Euler step
         #------------------------------------
-            
+
         # Compute fluxes across each element edge
         self.compute_fluxes()
@@ -1300,8 +1357,8 @@
         # of conserved quantities and cleanup
         #------------------------------------
-        
+
         # Combine steps
         self.saxpy_conserved_quantities(0.5, 0.5)
- 
+
         # Update ghosts
         self.update_ghosts()
@@ -1313,9 +1370,10 @@
         self.update_boundary()
 
-
-
+    ##
+    # @brief 'rk3' time step method.
+    # @param yieldstep The reporting time step.
+    # @param finaltime The simulation final time.
     def evolve_one_rk3_step(self, yieldstep, finaltime):
-        """
-        One 3rd order RK timestep
+        """One 3rd order RK timestep
         Q^(1) = 3/4 Q^n + 1/4 E(h)^2 Q^n  (at time t^n + h/2)
         Q^{n+1} = 1/3 Q^n + 2/3 E(h) Q^(1) (at time t^{n+1})
@@ -1323,8 +1381,8 @@
 
         # Save initial initial conserved quantities values
-        self.backup_conserved_quantities()            
+        self.backup_conserved_quantities()
 
         initial_time = self.time
-        
+
         #--------------------------------------
         # First euler step
@@ -1352,9 +1410,8 @@
         self.update_boundary()
 
-
         #------------------------------------
         # Second Euler step
         #------------------------------------
-            
+
         # Compute fluxes across each element edge
         self.compute_fluxes()
@@ -1370,9 +1427,8 @@
         # Combine steps
         self.saxpy_conserved_quantities(0.25, 0.75)
- 
+
         # Update ghosts
         self.update_ghosts()
 
-
         # Set substep time
         self.time = initial_time + self.timestep*0.5
@@ -1383,10 +1439,9 @@
         # Update boundary values
         self.update_boundary()
-
 
         #------------------------------------
         # Third Euler step
         #------------------------------------
-            
+
         # Compute fluxes across each element edge
         self.compute_fluxes()
@@ -1399,13 +1454,13 @@
         # and cleanup
         #------------------------------------
-        
+
         # Combine steps
         self.saxpy_conserved_quantities(2.0/3.0, 1.0/3.0)
- 
+
         # Update ghosts
         self.update_ghosts()
 
         # Set new time
-        self.time = initial_time + self.timestep       
+        self.time = initial_time + self.timestep
 
         # Update vertex and edge values
@@ -1415,12 +1470,15 @@
         self.update_boundary()
 
-        
-    def evolve_to_end(self, finaltime = 1.0):
-        """Iterate evolve all the way to the end      """
+    ##
+    # @brief Evolve simulation to a final time.
+    # @param finaltime Sinulation final time.
+    def evolve_to_end(self, finaltime=1.0):
+        """Iterate evolve all the way to the end."""
 
         for _ in self.evolve(yieldstep=None, finaltime=finaltime):
             pass
 
-
+    ##
+    # @brief Backup conserved quantities.
     def backup_conserved_quantities(self):
         N = len(self) # Number_of_triangles
@@ -1429,7 +1487,11 @@
         for name in self.conserved_quantities:
             Q = self.quantities[name]
-            Q.backup_centroid_values()        
-
-    def saxpy_conserved_quantities(self,a,b):
+            Q.backup_centroid_values()
+
+    ##
+    # @brief ??
+    # @param a ??
+    # @param b ??
+    def saxpy_conserved_quantities(self, a, b):
         N = len(self) #number_of_triangles
 
@@ -1437,12 +1499,13 @@
         for name in self.conserved_quantities:
             Q = self.quantities[name]
-            Q.saxpy_centroid_values(a,b)        
-    
-
+            Q.saxpy_centroid_values(a, b)
+
+    ##
+    # @brief Update boundary values for all conserved quantities.
     def update_boundary(self):
         """Go through list of boundary objects and update boundary values
         for all conserved quantities on boundary.
-        It is assumed that the ordering of conserved quantities is 
-        consistent between the domain and the boundary object, i.e. 
+        It is assumed that the ordering of conserved quantities is
+        consistent between the domain and the boundary object, i.e.
         the jth element of vector q must correspond to the jth conserved
         quantity in domain.
@@ -1453,6 +1516,7 @@
         for i, ((vol_id, edge_id), B) in enumerate(self.boundary_objects):
             if B is None:
-                print 'WARNING: Ignored boundary segment %d (None)'
+                print 'WARNING: Ignored boundary segment (None)'
             else:
+                print 'XXXX'
                 q = B.evaluate(vol_id, edge_id)
 
@@ -1461,36 +1525,37 @@
                     Q.boundary_values[i] = q[j]
 
-
+    ##
+    # @brief Compute fluxes.
+    # @note MUST BE OVERRIDEN IN SUBCLASS!
     def compute_fluxes(self):
         msg = 'Method compute_fluxes must be overridden by Domain subclass'
-        raise msg
-
-
+        raise Exception, msg
+
+    ##
+    # @brief 
+    # @param yieldstep 
+    # @param finaltime 
     def update_timestep(self, yieldstep, finaltime):
-
         from anuga.config import min_timestep, max_timestep
 
-        
-        
         # Protect against degenerate timesteps arising from isolated
         # triangles
         # FIXME (Steve): This should be in shallow_water as it assumes x and y
         # momentum
-        if self.protect_against_isolated_degenerate_timesteps is True and\
+        if self.protect_against_isolated_degenerate_timesteps is True and \
                self.max_speed > 10.0: # FIXME (Ole): Make this configurable
 
             # Setup 10 bins for speed histogram
             from anuga.utilities.numerical_tools import histogram, create_bins
-        
+
             bins = create_bins(self.max_speed, 10)
             hist = histogram(self.max_speed, bins)
 
             # Look for characteristic signature
-            if len(hist) > 1 and\
-                hist[-1] > 0 and\
+            if len(hist) > 1 and hist[-1] > 0 and \
                 hist[4] == hist[5] == hist[6] == hist[7] == hist[8] == 0:
                     # Danger of isolated degenerate triangles
-                    # print self.timestepping_statistics(track_speeds=True)  
-                    
+                    # print self.timestepping_statistics(track_speeds=True)
+
                     # Find triangles in last bin
                     # FIXME - speed up using Numeric
@@ -1498,22 +1563,19 @@
                     for i in range(self.number_of_full_triangles):
                         if self.max_speed[i] > bins[-1]:
-                            msg = 'Time=%f: Ignoring isolated high ' %self.time
-                            msg += 'speed triangle ' 
-                            msg += '#%d of %d with max speed=%f'\
-                                  %(i, self.number_of_full_triangles,
-                                    self.max_speed[i])
-                    
+                            msg = 'Time=%f: Ignoring isolated high ' % self.time
+                            msg += 'speed triangle '
+                            msg += '#%d of %d with max speed=%f' \
+                                      % (i, self.number_of_full_triangles,
+                                         self.max_speed[i])
+
                             # print 'Found offending triangle', i,
-                            # self.max_speed[i] 
-                            self.get_quantity('xmomentum').set_values(0.0, indices=[i])
-                            self.get_quantity('ymomentum').set_values(0.0, indices=[i])
+                            # self.max_speed[i]
+                            self.get_quantity('xmomentum').\
+                                            set_values(0.0, indices=[i])
+                            self.get_quantity('ymomentum').\
+                                            set_values(0.0, indices=[i])
                             self.max_speed[i]=0.0
                             d += 1
-                    
-                    #print 'Adjusted %d triangles' %d        
-                    #print self.timestepping_statistics(track_speeds=True)      
-                
-
-                        
+
         # self.timestep is calculated from speed of characteristics
         # Apply CFL condition here
@@ -1524,9 +1586,6 @@
         self.min_timestep = min(timestep, self.min_timestep)
 
-           
- 
         # Protect against degenerate time steps
         if timestep < min_timestep:
-
             # Number of consecutive small steps taken b4 taking action
             self.smallsteps += 1
@@ -1536,8 +1595,8 @@
 
                 if self._order_ == 1:
-                    msg = 'WARNING: Too small timestep %.16f reached '\
-                          %timestep
-                    msg += 'even after %d steps of 1 order scheme'\
-                           %self.max_smallsteps
+                    msg = 'WARNING: Too small timestep %.16f reached ' \
+                              % timestep
+                    msg += 'even after %d steps of 1 order scheme' \
+                               % self.max_smallsteps
                     print msg
                     timestep = min_timestep  # Try enforcing min_step
@@ -1549,5 +1608,4 @@
                     # Try to overcome situation by switching to 1 order
                     self._order_ = 1
-
         else:
             self.smallsteps = 0
@@ -1555,5 +1613,4 @@
                 self._order_ = 2
 
-
         # Ensure that final time is not exceeded
         if finaltime is not None and self.time + timestep > finaltime :
@@ -1566,6 +1623,6 @@
         self.timestep = timestep
 
-
-
+    ##
+    # @brief Compute forcing terms, if any.
     def compute_forcing_terms(self):
         """If there are any forcing functions driving the system
@@ -1577,5 +1634,6 @@
             f(self)
 
-
+    ##
+    # @brief Update vectors of conserved quantities.
     def update_conserved_quantities(self):
         """Update vectors of conserved quantities using previously
@@ -1597,11 +1655,15 @@
 
             # Note that Q.explicit_update is reset by compute_fluxes
-            # Where is Q.semi_implicit_update reset? 
+            # Where is Q.semi_implicit_update reset?
             # It is reset in quantity_ext.c
-            
-
+
+    ##
+    # @brief ??
     def update_ghosts(self):
         pass
 
+    ##
+    # @brief Extrapolate conserved quantities from centroid to vertices
+    #        and edge-midpoints for each volume.
     def distribute_to_vertices_and_edges(self):
         """Extrapolate conserved quantities from centroid to
@@ -1621,8 +1683,12 @@
                 #Q.limit()
             else:
-                raise 'Unknown order'
+                raise Exception, 'Unknown order'
             #Q.interpolate_from_vertices_to_edges()
 
-
+    ##
+    # @brief Calculate the norm of the centroid values of a specific quantity,
+    #        using normfunc.
+    # @param quantity 
+    # @param normfunc 
     def centroid_norm(self, quantity, normfunc):
         """Calculate the norm of the centroid values
@@ -1633,6 +1699,6 @@
         common normfuncs are provided in the module utilities.norms
         """
+
         return normfunc(self.quantities[quantity].centroid_values)
-
 
 
@@ -1652,6 +1718,6 @@
             # Psyco isn't supported on 64 bit systems, but it doesn't matter
         else:
-            msg = 'WARNING: psyco (speedup) could not import'+\
-                  ', you may want to consider installing it'
+            msg = ('WARNING: psyco (speedup) could not be imported, '
+                   'you may want to consider installing it')
             print msg
     else:
