Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py	(revision 6225)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py	(revision 6226)
@@ -8,7 +8,9 @@
 """
 
+import types
+from time import time as walltime
+
 from anuga.config import epsilon
 from anuga.config import beta_euler, beta_rk2
-
 from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
 from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
@@ -22,22 +24,16 @@
 from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
      import Transmissive_boundary
-
 from anuga.abstract_2d_finite_volumes.pmesh2domain import pmesh_to_domain
 from anuga.abstract_2d_finite_volumes.region\
      import Set_region as region_set_region
-
 from anuga.utilities.polygon import inside_polygon
 from anuga.abstract_2d_finite_volumes.util import get_textual_float
-
 from quantity import Quantity
 
-import types
-from time import time as walltime
-
 import Numeric as num
 
 
 ##
-# @brief Basic Domain class
+# @brief Generic Domain class
 class Domain:
 
@@ -49,16 +45,16 @@
     # @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
+    # @param tagged_elements ??
+    # @param geo_reference ??
+    # @param use_inscribed_circle ??
+    # @param mesh_filename ??
+    # @param use_cache ??
+    # @param verbose True if this method is to be 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,
@@ -119,27 +115,27 @@
                          number_of_full_triangles=number_of_full_triangles,
                          verbose=verbose)
-                         
+
         # Expose Mesh attributes (FIXME: Maybe turn into methods)
         self.centroid_coordinates = self.mesh.centroid_coordinates
-        self.vertex_coordinates = self.mesh.vertex_coordinates        
+        self.vertex_coordinates = self.mesh.vertex_coordinates
         self.boundary = self.mesh.boundary
         self.neighbours = self.mesh.neighbours
-        self.surrogate_neighbours = self.mesh.surrogate_neighbours        
+        self.surrogate_neighbours = self.mesh.surrogate_neighbours
         self.neighbour_edges = self.mesh.neighbour_edges
         self.normals = self.mesh.normals
-        self.edgelengths = self.mesh.edgelengths        
-        self.radii = self.mesh.radii                
-        self.areas = self.mesh.areas                        
-                
-        self.number_of_boundaries = self.mesh.number_of_boundaries        
+        self.edgelengths = self.mesh.edgelengths
+        self.radii = self.mesh.radii
+        self.areas = self.mesh.areas
+
+        self.number_of_boundaries = self.mesh.number_of_boundaries
         self.number_of_full_nodes = self.mesh.number_of_full_nodes
-        self.number_of_full_triangles = self.mesh.number_of_full_triangles        
+        self.number_of_full_triangles = self.mesh.number_of_full_triangles
         self.number_of_triangles_per_node = self.mesh.number_of_triangles_per_node
 
         self.vertex_value_indices = self.mesh.vertex_value_indices
-        self.number_of_triangles = self.mesh.number_of_triangles        
+        self.number_of_triangles = self.mesh.number_of_triangles
 
         self.geo_reference = self.mesh.geo_reference
-        
+
         if verbose: print 'Initialising Domain'
 
@@ -276,8 +272,9 @@
 
         if verbose: print 'Domain: Done'
-        
-     
-
+
+    ######
     # Expose underlying Mesh functionality
+    ######
+
     def __len__(self):
         return len(self.mesh)
@@ -285,10 +282,10 @@
     def get_centroid_coordinates(self, *args, **kwargs):
         return self.mesh.get_centroid_coordinates(*args, **kwargs)
-        
+
     def get_radii(self, *args, **kwargs):
-        return self.mesh.get_radii(*args, **kwargs)        
-        
+        return self.mesh.get_radii(*args, **kwargs)
+
     def get_areas(self, *args, **kwargs):
-        return self.mesh.get_areas(*args, **kwargs)                
+        return self.mesh.get_areas(*args, **kwargs)
 
     def get_area(self, *args, **kwargs):
@@ -296,24 +293,24 @@
 
     def get_vertex_coordinates(self, *args, **kwargs):
-        return self.mesh.get_vertex_coordinates(*args, **kwargs)                
+        return self.mesh.get_vertex_coordinates(*args, **kwargs)
 
     def get_triangles(self, *args, **kwargs):
-        return self.mesh.get_triangles(*args, **kwargs)                
-        
+        return self.mesh.get_triangles(*args, **kwargs)
+
     def get_nodes(self, *args, **kwargs):
         return self.mesh.get_nodes(*args, **kwargs)
-        
+
     def get_number_of_nodes(self, *args, **kwargs):
         return self.mesh.get_number_of_nodes(*args, **kwargs)
-        
+
     def get_normal(self, *args, **kwargs):
-        return self.mesh.get_normal(*args, **kwargs)        
-        
+        return self.mesh.get_normal(*args, **kwargs)
+
     def get_intersecting_segments(self, *args, **kwargs):
         return self.mesh.get_intersecting_segments(*args, **kwargs)
-        
+
     def get_disconnected_triangles(self, *args, **kwargs):
         return self.mesh.get_disconnected_triangles(*args, **kwargs)
-        
+
     def get_boundary_tags(self, *args, **kwargs):
         return self.mesh.get_boundary_tags(*args, **kwargs)
@@ -321,37 +318,35 @@
     def get_boundary_polygon(self, *args, **kwargs):
         return self.mesh.get_boundary_polygon(*args, **kwargs)
-                
+
     def get_number_of_triangles_per_node(self, *args, **kwargs):
         return self.mesh.get_number_of_triangles_per_node(*args, **kwargs)
-        
+
     def get_triangles_and_vertices_per_node(self, *args, **kwargs):
         return self.mesh.get_triangles_and_vertices_per_node(*args, **kwargs)
-        
+
     def get_interpolation_object(self, *args, **kwargs):
-        return self.mesh.get_interpolation_object(*args, **kwargs)        
-        
+        return self.mesh.get_interpolation_object(*args, **kwargs)
+
     def get_tagged_elements(self, *args, **kwargs):
-        return self.mesh.get_tagged_elements(*args, **kwargs)                
-        
+        return self.mesh.get_tagged_elements(*args, **kwargs)
+
     def get_lone_vertices(self, *args, **kwargs):
-        return self.mesh.get_lone_vertices(*args, **kwargs)        
-        
+        return self.mesh.get_lone_vertices(*args, **kwargs)
+
     def get_unique_vertices(self, *args, **kwargs):
-        return self.mesh.get_unique_vertices(*args, **kwargs)                
+        return self.mesh.get_unique_vertices(*args, **kwargs)
 
     def get_georeference(self, *args, **kwargs):
         return self.mesh.get_georeference(*args, **kwargs)
-        
+
     def set_georeference(self, *args, **kwargs):
-        self.mesh.set_georeference(*args, **kwargs)                    
-        
+        self.mesh.set_georeference(*args, **kwargs)
+
     def build_tagged_elements_dictionary(self, *args, **kwargs):
         self.mesh.build_tagged_elements_dictionary(*args, **kwargs)
-        
+
     def statistics(self, *args, **kwargs):
-        return self.mesh.statistics(*args, **kwargs)        
-                
-        
-        
+        return self.mesh.statistics(*args, **kwargs)
+
     ##
     # @brief Get conserved quantities for a volume.
@@ -365,5 +360,5 @@
                                        vertex=None,
                                        edge=None):
-        """Get conserved quantities at volume vol_id
+        """Get conserved quantities at volume vol_id.
 
         If vertex is specified use it as index for vertex values
@@ -397,5 +392,6 @@
     # @param time The new model time (seconds).
     def set_time(self, time=0.0):
-        """Set the model time (seconds)"""
+        """Set the model time (seconds)."""
+
         # FIXME: this is setting the relative time
         # Note that get_time and set_time are now not symmetric
@@ -407,5 +403,5 @@
     # @return The absolute model time (seconds).
     def get_time(self):
-        """Get the absolute model time (seconds)"""
+        """Get the absolute model time (seconds)."""
 
         return self.time + self.starttime
@@ -415,9 +411,8 @@
     # @param beta The new beta value.
     def set_beta(self, beta):
-        """Set default beta for limiting"""
+        """Set default beta for limiting."""
 
         self.beta = beta
         for name in self.quantities:
-            #print 'setting beta for quantity ',name
             Q = self.quantities[name]
             Q.set_beta(beta)
@@ -427,5 +422,5 @@
     # @return The beta value used for limiting.
     def get_beta(self):
-        """Get default beta for limiting"""
+        """Get default beta for limiting."""
 
         return self.beta
@@ -436,5 +431,5 @@
     # @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"""
+        """Set default (spatial) order to either 1 or 2."""
 
         msg = 'Default order must be either 1 or 2. I got %s' % n
@@ -443,14 +438,4 @@
         self.default_order = n
         self._order_ = self.default_order
-
-        if self.default_order == 1:
-            pass
-            #self.set_timestepping_method('euler')
-            #self.set_all_limiters(beta_euler)
-
-        if self.default_order == 2:
-            pass
-            #self.set_timestepping_method('rk2')
-            #self.set_all_limiters(beta_rk2)
 
     ##
@@ -681,6 +666,6 @@
     ##
     # @brief Set quantities based on a regional tag.
-    # @param args 
-    # @param kwargs 
+    # @param args
+    # @param kwargs
     def set_region(self, *args, **kwargs):
         """Set quantities based on a regional tag.
@@ -1228,7 +1213,7 @@
         self.starttime = float(time)
 
-#--------------------------
+################################################################################
 # Main components of evolve
-#--------------------------
+################################################################################
 
     ##
@@ -1275,6 +1260,6 @@
                'evolving system, '
                'e.g. using the method set_boundary.\n'
-               'This system has the boundary tags %s ') \
-                   % self.get_boundary_tags()
+               'This system has the boundary tags %s '
+                   % self.get_boundary_tags())
         assert hasattr(self, 'boundary_objects'), msg
 
@@ -1287,5 +1272,4 @@
 
         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 Exception, msg
@@ -1296,6 +1280,6 @@
                 self.finaltime = self.starttime + float(duration)
 
-        N = len(self) # Number of triangles
-        self.yieldtime = 0.0 # Track time between 'yields'
+        N = len(self)            # Number of triangles
+        self.yieldtime = 0.0     # Track time between 'yields'
 
         # Initialise interval of timestep sizes (for reporting only)
@@ -1322,5 +1306,5 @@
 
         if skip_initial_step is False:
-            yield(self.time)  # Yield initial values
+            yield(self.time)      # Yield initial values
 
         while True:
@@ -1347,6 +1331,6 @@
                 if self.time > finaltime:
                     # FIXME (Ole, 30 April 2006): Do we need this check?
-                    # Probably not (Ole, 18 September 2008). Now changed to
-                    # Exception
+                    # Probably not (Ole, 18 September 2008).
+                    # Now changed to Exception.
                     msg = ('WARNING (domain.py): time overshot finaltime. '
                            'Contact Ole.Nielsen@ga.gov.au')
@@ -1417,7 +1401,7 @@
         self.backup_conserved_quantities()
 
-        #--------------------------------------
+        ######
         # First euler step
-        #--------------------------------------
+        ######
 
         # Compute fluxes across each element edge
@@ -1442,7 +1426,7 @@
         self.update_boundary()
 
-        #------------------------------------
+        ######
         # Second Euler step
-        #------------------------------------
+        ######
 
         # Compute fluxes across each element edge
@@ -1452,8 +1436,8 @@
         self.update_conserved_quantities()
 
-        #------------------------------------
+        ######
         # Combine initial and final values
         # of conserved quantities and cleanup
-        #------------------------------------
+        ######
 
         # Combine steps
@@ -1484,7 +1468,7 @@
         initial_time = self.time
 
-        #--------------------------------------
+        ######
         # First euler step
-        #--------------------------------------
+        ######
 
         # Compute fluxes across each element edge
@@ -1509,7 +1493,7 @@
         self.update_boundary()
 
-        #------------------------------------
+        ######
         # Second Euler step
-        #------------------------------------
+        ######
 
         # Compute fluxes across each element edge
@@ -1519,8 +1503,8 @@
         self.update_conserved_quantities()
 
-        #------------------------------------
-        #Combine steps to obtain intermediate
-        #solution at time t^n + 0.5 h
-        #------------------------------------
+        ######
+        # Combine steps to obtain intermediate
+        # solution at time t^n + 0.5 h
+        ######
 
         # Combine steps
@@ -1539,7 +1523,7 @@
         self.update_boundary()
 
-        #------------------------------------
+        ######
         # Third Euler step
-        #------------------------------------
+        ######
 
         # Compute fluxes across each element edge
@@ -1549,8 +1533,8 @@
         self.update_conserved_quantities()
 
-        #------------------------------------
+        ######
         # Combine final and initial values
         # and cleanup
-        #------------------------------------
+        ######
 
         # Combine steps
@@ -1631,7 +1615,7 @@
 
     ##
-    # @brief 
-    # @param yieldstep 
-    # @param finaltime 
+    # @brief
+    # @param yieldstep
+    # @param finaltime
     def update_timestep(self, yieldstep, finaltime):
         from anuga.config import min_timestep, max_timestep
@@ -1779,17 +1763,15 @@
             elif self._order_ == 2:
                 Q.extrapolate_second_order()
-                #Q.limit()
             else:
                 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 
+    # @param quantity
+    # @param normfunc
     def centroid_norm(self, quantity, normfunc):
-        """Calculate the norm of the centroid values
-        of a specific quantity, using normfunc.
+        """Calculate the norm of the centroid values of a specific quantity,
+        using normfunc.
 
         normfunc should take a list to a float.
@@ -1801,10 +1783,11 @@
 
 
-#------------------
+######
 # Initialise module
-#------------------
+######
 
 # Optimisation with psyco
 from anuga.config import use_psyco
+
 if use_psyco:
     try:
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py	(revision 6225)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py	(revision 6226)
@@ -17,28 +17,31 @@
 from anuga.utilities.numerical_tools import ensure_numeric, is_scalar
 from anuga.utilities.polygon import inside_polygon
-
 from anuga.geospatial_data.geospatial_data import Geospatial_data
 from anuga.fit_interpolate.fit import fit_to_mesh
 from anuga.config import points_file_block_line_size as default_block_line_size
 from anuga.config import epsilon
-from anuga.caching import cache
-
-
 
 import Numeric as num
 
 
+##
+# @brief Implement values at each triangular element.
 class Quantity:
 
-    def __init__(self, domain, vertex_values=None):
-    
+    ##
+    # @brief Construct values art each triangular element.
+    # @param domain ??
+    # @param vertex_values ??
+    def __init__(self, domain,
+                       vertex_values=None):
         from anuga.abstract_2d_finite_volumes.domain import Domain
+
         msg = 'First argument in Quantity.__init__ '
         msg += 'must be of class Domain (or a subclass thereof). '
-        msg += 'I got %s.' %str(domain.__class__)
+        msg += 'I got %s.' % str(domain.__class__)
         assert isinstance(domain, Domain), msg
 
         if vertex_values is None:
-            N = len(domain) # number_of_elements
+            N = len(domain)             # number_of_elements
             self.vertex_values = num.zeros((N, 3), num.Float)
         else:
@@ -46,13 +49,8 @@
 
             N, V = self.vertex_values.shape
-            assert V == 3,\
-                   'Three vertex values per element must be specified'
-
-
-            msg = 'Number of vertex values (%d) must be consistent with'\
-                  %N
-            msg += 'number of elements in specified domain (%d).'\
-                   %len(domain)
-
+            assert V == 3, 'Three vertex values per element must be specified'
+
+            msg = 'Number of vertex values (%d) must be consistent with' % N
+            msg += 'number of elements in specified domain (%d).' % len(domain)
             assert N == len(domain), msg
 
@@ -68,5 +66,5 @@
 
         # Allocate space for Limiter Phi
-        self.phi = num.zeros(N, num.Float)        
+        self.phi = num.zeros(N, num.Float)
 
         # Intialise centroid and edge_values
@@ -87,11 +85,16 @@
         self.set_beta(1.0)
 
-
-
+    #
     # Methods for operator overloading
+    #
+
+    ##
+    # @brief Get length of object.
+    # @return Return axis length of array.
     def __len__(self):
         return self.centroid_values.shape[0]
 
-
+    ##
+    # @brief Negate all values in this quantity.
     def __neg__(self):
         """Negate all values in this quantity giving meaning to the
@@ -103,5 +106,7 @@
         return Q
 
-
+    ##
+    # @brief Add two quantities.
+    # @param other The other quantity (to the right).
     def __add__(self, other):
         """Add to self anything that could populate a quantity
@@ -119,13 +124,22 @@
         return result
 
+    ##
+    # @brief Add two quantities.
+    # @param other The other quantity (to the left).
     def __radd__(self, other):
         """Handle cases like 7+Q, where Q is an instance of class Quantity
         """
+
         return self + other
 
-
+    ##
+    # @brief Subtract two quantities.
+    # @param other The other quantity (to the right).
     def __sub__(self, other):
-        return self + -other  #Invoke __neg__
-
+        return self + -other            # Invoke self.__neg__()
+
+    ##
+    # @brief Multiply two quantities.
+    # @param other The other quantity (to the right).
     def __mul__(self, other):
         """Multiply self with anything that could populate a quantity
@@ -138,5 +152,5 @@
         if isinstance(other, Quantity):
             Q = other
-        else:    
+        else:
             Q = Quantity(self.domain)
             Q.set_values(other)
@@ -151,12 +165,19 @@
         result.edge_values = self.edge_values * Q.edge_values
         result.centroid_values = self.centroid_values * Q.centroid_values
-        
+
         return result
 
+    ##
+    # @brief Multiply two quantities.
+    # @param other The other quantity (to the left).
     def __rmul__(self, other):
         """Handle cases like 3*Q, where Q is an instance of class Quantity
         """
+
         return self * other
 
+    ##
+    # @brief Divide two quantities.
+    # @param other The other quantity (to the right).
     def __div__(self, other):
         """Divide self with anything that could populate a quantity
@@ -172,5 +193,5 @@
         if isinstance(other, Quantity):
             Q = other
-        else:    
+        else:
             Q = Quantity(self.domain)
             Q.set_values(other)
@@ -188,9 +209,16 @@
         return result
 
+    ##
+    # @brief Divide two quantities.
+    # @param other The other quantity (to the left).
     def __rdiv__(self, other):
         """Handle cases like 3/Q, where Q is an instance of class Quantity
         """
+
         return self / other
 
+    ##
+    # @brief Raise a quaintity to some power.
+    # @param other The power to raise quantity.
     def __pow__(self, other):
         """Raise quantity to (numerical) power
@@ -201,10 +229,9 @@
         Example using __pow__:
           Q = (Q1**2 + Q2**2)**0.5
-
         """
 
         if isinstance(other, Quantity):
             Q = other
-        else:    
+        else:
             Q = Quantity(self.domain)
             Q.set_values(other)
@@ -222,12 +249,9 @@
         return result
 
-    #def __sqrt__(self, other):
-    #    """Define in terms of x**0.5
-    #    """
-    #    pass
-
-    def set_beta(self,beta):
-        """Set default beta value for limiting
-        """
+    ##
+    # @brief Set default beta value for limiting.
+    # @param beta ??
+    def set_beta(self, beta):
+        """Set default beta value for limiting """
 
         if beta < 0.0:
@@ -235,23 +259,28 @@
         if beta > 2.0:
             print 'WARNING: setting beta > 2.0'
-            
+
         self.beta = beta
 
+    ##
+    # @brief Get the current beta value.
+    # @return The current beta value.
     def get_beta(self):
-        """Get default beta value for limiting
-        """
+        """Get default beta value for limiting"""
 
         return self.beta
 
+    ##
+    # @brief Compute interpolated values at edges and centroid.
+    # @note vertex_values must be set before calling this.
     def interpolate(self):
         """Compute interpolated values at edges and centroid
         Pre-condition: vertex_values have been set
         """
-        
+
         # FIXME (Ole): Maybe this function
         # should move to the C-interface?
         # However, it isn't called by validate_all.py, so it
         # may not be that important to optimise it?
-        
+
         N = self.vertex_values.shape[0]
         for i in range(N):
@@ -264,35 +293,50 @@
         self.interpolate_from_vertices_to_edges()
 
-
+    ##
+    # @brief ??
     def interpolate_from_vertices_to_edges(self):
-        # Call correct module function
-        # (either from this module or C-extension)
+        # Call correct module function (either from this module or C-extension)
         interpolate_from_vertices_to_edges(self)
 
+    ##
+    # @brief ??
     def interpolate_from_edges_to_vertices(self):
-        # Call correct module function
-        # (either from this module or C-extension)
+        # Call correct module function (either from this module or C-extension)
         interpolate_from_edges_to_vertices(self)
-
-
-
 
     #---------------------------------------------
     # Public interface for setting quantity values
     #---------------------------------------------
-    def set_values(self,
-                   numeric=None,    # List, numeric array or constant
-                   quantity=None,   # Another quantity
-                   function=None,   # Callable object: f(x,y)
-                   geospatial_data=None, # Arbitrary dataset
-                   filename=None, attribute_name=None, # Input from file
-                   alpha=None,
-                   location='vertices',
-                   polygon=None,
-                   indices=None,
-                   smooth=False,
-                   verbose=False,
-                   use_cache=False):
-
+
+    ##
+    # @brief Set values for quantity based on different sources.
+    # @param numeric A num array, list or constant value.
+    # @param quantity Another Quantity.
+    # @param function Any callable object that takes two 1d arrays.
+    # @param geospatial_data Arbitrary instance of class Geospatial_data
+    # @param filename Path to a points file.
+    # @param attribute_name If specified any array using that name will be used.
+    # @param alpha Smoothing parameter to be used with fit_interpolate.fit.
+    # @param location Where to store values (vertices, edges, centroids).
+    # @param polygon Restrict update to locations that fall inside polygon.
+    # @param indices Restrict update to locations specified by this.
+    # @param smooth If True, smooth vertex values.
+    # @param verbose True if this method is to be verbose.
+    # @param use_cache If True cache results for fit_interpolate.fit.
+    # @note Exactly one of 'numeric', 'quantity', 'function', 'filename'
+    #       must be present.
+    def set_values(self, numeric=None,         # List, numeric array or constant
+                         quantity=None,        # Another quantity
+                         function=None,        # Callable object: f(x,y)
+                         geospatial_data=None, # Arbitrary dataset
+                         filename=None,
+                         attribute_name=None,  # Input from file
+                         alpha=None,
+                         location='vertices',
+                         polygon=None,
+                         indices=None,
+                         smooth=False,
+                         verbose=False,
+                         use_cache=False):
         """Set values for quantity based on different sources.
 
@@ -366,8 +410,8 @@
                  constant values so far.
 
-        indices: Restrict update of quantity to locations that are 
+        indices: Restrict update of quantity to locations that are
                  identified by indices (e.g. node ids if location
                  is 'unique vertices' or triangle ids otherwise).
-        
+
         verbose: True means that output to stdout is generated
 
@@ -391,9 +435,9 @@
         # perhaps the notion of location and indices simplified
 
-        # FIXME (Ole): Need to compute indices based on polygon 
+        # FIXME (Ole): Need to compute indices based on polygon
         # (and location) and use existing code after that.
-        
+
         # See ticket:275, ticket:250, ticeket:254 for refactoring plan
-        
+
         if polygon is not None:
             if indices is not None:
@@ -405,104 +449,71 @@
             msg += 'a constant.'
             if numeric is None:
-                raise Exception, msg            
+                raise Exception, msg
             else:
                 # Check that numeric is as constant
                 assert type(numeric) in [FloatType, IntType, LongType], msg
 
-
             location = 'centroids'
-
 
             points = self.domain.get_centroid_coordinates(absolute=True)
             indices = inside_polygon(points, polygon)
-            
-            self.set_values_from_constant(numeric,
-                                          location, indices, verbose)
-
+
+            self.set_values_from_constant(numeric, location, indices, verbose)
 
             self.extrapolate_first_order()
 
             if smooth:
-                self.smooth_vertex_values(use_cache=use_cache,
-                                          verbose=verbose)
-
-                
+                self.smooth_vertex_values()
+
             return
-        
-        
-
-
-
 
         # General input checks
         L = [numeric, quantity, function, geospatial_data, filename]
-        msg = 'Exactly one of the arguments '+\
-              'numeric, quantity, function, geospatial_data, '+\
-              'or filename must be present.'
+        msg = ('Exactly one of the arguments numeric, quantity, function, '
+               'geospatial_data, or filename must be present.')
         assert L.count(None) == len(L)-1, msg
-
 
         if location == 'edges':
             msg = 'edges has been deprecated as valid location'
             raise Exception, msg
-            
+
         if location not in ['vertices', 'centroids', 'unique vertices']:
-            msg = 'Invalid location: %s' %location
+            msg = 'Invalid location: %s' % location
             raise Exception, msg
-
 
         msg = 'Indices must be a list or None'
         assert type(indices) in [ListType, NoneType, num.ArrayType], msg
 
-
-
         # Determine which 'set_values_from_...' to use
-
         if numeric is not None:
             if type(numeric) in [FloatType, IntType, LongType]:
-                self.set_values_from_constant(numeric,
-                                              location, indices, verbose)
+                self.set_values_from_constant(numeric, location,
+                                              indices, verbose)
             elif type(numeric) in [num.ArrayType, ListType]:
-                self.set_values_from_array(numeric,
-                                           location, indices,
-                                           use_cache=use_cache,
-                                           verbose=verbose)
+                self.set_values_from_array(numeric, location, indices, verbose)
             elif callable(numeric):
-                self.set_values_from_function(numeric,
-                                              location, indices,
-                                              use_cache=use_cache,
-                                              verbose=verbose)
+                self.set_values_from_function(numeric, location,
+                                              indices, verbose)
             elif isinstance(numeric, Quantity):
-                self.set_values_from_quantity(numeric,
-                                              location, indices,
-                                              verbose=verbose)
+                self.set_values_from_quantity(numeric, location,
+                                              indices, verbose)
             elif isinstance(numeric, Geospatial_data):
-                self.set_values_from_geospatial_data(numeric,
-                                                     alpha,
+                self.set_values_from_geospatial_data(numeric, alpha, location,
+                                                     indices, verbose=verbose,
+                                                     use_cache=use_cache)
+            else:
+                msg = 'Illegal type for argument numeric: %s' % str(numeric)
+                raise msg
+        elif quantity is not None:
+            self.set_values_from_quantity(quantity, location, indices, verbose)
+        elif function is not None:
+            msg = 'Argument function must be callable'
+            assert callable(function), msg
+            self.set_values_from_function(function, location, indices, verbose)
+        elif geospatial_data is not None:
+                self.set_values_from_geospatial_data(geospatial_data, alpha,
                                                      location, indices,
                                                      verbose=verbose,
                                                      use_cache=use_cache)
-            else:
-                msg = 'Illegal type for argument numeric: %s' %str(numeric)
-                raise msg
-
-        elif quantity is not None:
-            self.set_values_from_quantity(quantity,
-                                          location, indices, verbose)
-        elif function is not None:
-            msg = 'Argument function must be callable'
-            assert callable(function), msg
-            self.set_values_from_function(function,
-                                          location,
-                                          indices,
-                                          use_cache=use_cache,
-                                          verbose=verbose)
-        elif geospatial_data is not None:
-                self.set_values_from_geospatial_data(geospatial_data,
-                                                     alpha,
-                                                     location, indices,
-                                                     verbose=verbose,
-                                                     use_cache=use_cache)
-            
         elif filename is not None:
             if hasattr(self.domain, 'points_file_block_line_size'):
@@ -510,13 +521,10 @@
             else:
                 max_read_lines = default_block_line_size
-            self.set_values_from_file(filename, attribute_name, alpha,
-                                      location, indices,
-                                      verbose=verbose,
+            self.set_values_from_file(filename, attribute_name, alpha, location,
+                                      indices, verbose=verbose,
                                       max_read_lines=max_read_lines,
                                       use_cache=use_cache)
         else:
-            raise Exception, 'This can\'t happen :-)'
-
-
+            raise Exception, "This can't happen :-)"
 
         # Update all locations in triangles
@@ -531,17 +539,20 @@
 
 
-    #-------------------------------------------------------------        
+    #---------------------------------------------------------------------------
     # Specific internal functions for setting values based on type
-    #-------------------------------------------------------------            
-    
-    def set_values_from_constant(self, X,
-                                 location, indices, verbose):
-        """Set quantity values from specified constant X
-        """
+    #---------------------------------------------------------------------------
+
+    ##
+    # @brief Set quantity values from specified constant.
+    # @param X The constant to set quantity values to.
+    # @param location 
+    # @param indices 
+    # @param verbose 
+    def set_values_from_constant(self, X, location, indices, verbose):
+        """Set quantity values from specified constant X"""
 
         # FIXME (Ole): Somehow indices refer to centroids
         # rather than vertices as default. See unit test
         # test_set_vertex_values_using_general_interface_with_subset(self):
-        
 
         if location == 'centroids':
@@ -552,27 +563,17 @@
                 for i in indices:
                     self.centroid_values[i] = X
-
-        #elif location == 'edges':
-        #    if indices is None:
-        #        self.edge_values[:] = X
-        #    else:
-        #        # Brute force
-        #        for i in indices:
-        #            self.edge_values[i] = X
-
         elif location == 'unique vertices':
             if indices is None:
                 self.edge_values[:] = X  #FIXME (Ole): Shouldn't this be vertex_values?
             else:
-
                 # Go through list of unique vertices
                 for unique_vert_id in indices:
-
-                    triangles = self.domain.get_triangles_and_vertices_per_node(node=unique_vert_id)
-                    
+                    triangles = \
+                        self.domain.get_triangles_and_vertices_per_node(node=unique_vert_id)
+
                     # In case there are unused points
                     if len(triangles) == 0:
                         continue
-                    
+
                     # Go through all triangle, vertex pairs
                     # and set corresponding vertex value
@@ -590,12 +591,14 @@
                     self.vertex_values[i_vertex] = X
 
-
-
-
+    ##
+    # @brief Set values for a quantity.
+    # @param values Array of values.
+    # @param location Where values are to be stored.
+    # @param indices Limit update to these indices.
+    # @param verbose True if this method is to be verbose.
     def set_values_from_array(self, values,
-                              location='vertices',
-                              indices=None,
-                              use_cache=False,
-                              verbose=False):
+                                    location='vertices',
+                                    indices=None,
+                                    verbose=False):
         """Set values for quantity
 
@@ -628,7 +631,7 @@
         if indices is not None:
             indices = num.array(indices, num.Int)
-            msg = 'Number of values must match number of indices:'
-            msg += ' You specified %d values and %d indices'\
-                   %(values.shape[0], indices.shape[0])
+            msg = ('Number of values must match number of indices: You '
+                   'specified %d values and %d indices'
+                   % (values.shape[0], indices.shape[0]))
             assert values.shape[0] == indices.shape[0], msg
 
@@ -650,27 +653,15 @@
                 for i in range(len(indices)):
                     self.centroid_values[indices[i]] = values[i]
-
         elif location == 'unique vertices':
-            assert len(values.shape) == 1 or num.allclose(values.shape[1:], 1),\
-                   'Values array must be 1d'
-
-            self.set_vertex_values(values.flat,
-                                   indices=indices,
-                                   use_cache=use_cache,
-                                   verbose=verbose)
-            
+            assert (len(values.shape) == 1 or num.allclose(values.shape[1:], 1),
+                    'Values array must be 1d')
+
+            self.set_vertex_values(values.flat, indices=indices)
         else:
             # Location vertices
             if len(values.shape) == 1:
-                # This is the common case arising from fitted
-                # values (e.g. from pts file).
-                self.set_vertex_values(values,
-                                       indices=indices,
-                                       use_cache=use_cache,
-                                       verbose=verbose)
-
+                self.set_vertex_values(values, indices=indices)
             elif len(values.shape) == 2:
                 # Vertex values are given as a triplet for each triangle
-
                 msg = 'Array must be N x 3'
                 assert values.shape[1] == 3, msg
@@ -683,9 +674,13 @@
             else:
                 msg = 'Values array must be 1d or 2d'
-                raise msg
-            
-
-    def set_values_from_quantity(self, q,
-                                 location, indices, verbose):
+                raise Exception, msg
+
+    ##
+    # @brief Set quantity values from a specified quantity instance.
+    # @param q The quantity instance to take values from.
+    # @param location IGNORED, 'vertices' ALWAYS USED!
+    # @param indices ??
+    # @param verbose True if this method is to be verbose.
+    def set_values_from_quantity(self, q, location, indices, verbose):
         """Set quantity values from specified quantity instance q
 
@@ -702,17 +697,19 @@
 
         self.set_values(A, location='vertices',
-                        indices=indices,
-                        verbose=verbose)
-
-
+                        indices=indices, verbose=verbose)
+
+    ##
+    # @brief Set quantity values from a specified quantity instance.
+    # @param f Callable that takes two 1d array -> 1d array.
+    # @param location Where values are to be stored.
+    # @param indices ??
+    # @param verbose True if this method is to be verbose.
     def set_values_from_function(self, f,
-                                 location='vertices',
-                                 indices=None,
-                                 use_cache=False,
-                                 verbose=False):
+                                       location='vertices',
+                                       indices=None,
+                                       verbose=False):
         """Set values for quantity using specified function
 
         Input
-        
         f: x, y -> z Function where x, y and z are arrays
         location: Where values are to be stored.
@@ -720,14 +717,11 @@
                   unique vertices
                   Default is "vertices"
-        indices:  
-
-                  
+        indices:
         """
 
         # FIXME: Should check that function returns something sensible and
-        # raise a meaningfull exception if it returns None for example
+        # raise a meaningful exception if it returns None for example
 
         # FIXME: Should supply absolute coordinates
-
 
         # Compute the function values and call set_values again
@@ -735,36 +729,15 @@
             if indices is None:
                 indices = range(len(self))
-                
+
             V = num.take(self.domain.get_centroid_coordinates(), indices)
-
-            x = V[:,0]; y = V[:,1]
-            if use_cache is True:
-                res = cache(f, (x, y),
-                            verbose=verbose)
-            else:
-                res = f(x, y)
-
-            self.set_values(res,
-                            location=location,
-                            indices=indices)
-            
+            self.set_values(f(V[:,0], V[:,1]),
+                            location=location, indices=indices)
         elif location == 'vertices':
-            # This is the default branch taken by set_quantity
-            
             M = self.domain.number_of_triangles
             V = self.domain.get_vertex_coordinates()
 
-            x = V[:,0]; y = V[:,1]
-            if use_cache is True:
-                #print 'Caching function'
-                values = cache(f, (x, y),
-                               verbose=verbose)                
-            else:
-                if verbose is True:
-                    print 'Evaluating function in set_values'
-                values = f(x, y)
-
-            #print 'value', min(x), max(x), min(y), max(y), max(values), len(values)
-                
+            x = V[:,0];
+            y = V[:,1];
+            values = f(x, y)
 
             # FIXME (Ole): This code should replace all the
@@ -773,40 +746,41 @@
             # If that could be resolved this one will be
             # more robust and simple.
-            
-            #values = num.reshape(values, (M,3))
-            #self.set_values(values,
-            #                location='vertices',
-            #                indices=indices)
-
 
             # This should be removed
             if is_scalar(values):
                 # Function returned a constant value
-                self.set_values_from_constant(values,
-                                              location, indices, verbose)
+                self.set_values_from_constant(values, location,
+                                              indices, verbose)
                 return
 
-            # This should be removed            
+            # This should be removed
             if indices is None:
                 for j in range(3):
-                    self.vertex_values[:,j] = values[j::3]                 
-            else:    
+                    self.vertex_values[:, j] = values[j::3]
+            else:
                 # Brute force
                 for i in indices:
                     for j in range(3):
-                        self.vertex_values[i,j] = values[3*i+j]
-
-
+                        self.vertex_values[i, j] = values[3*i + j]
         else:
-            raise 'Not implemented: %s' %location
-
-
-
-    def set_values_from_geospatial_data(self, geospatial_data, alpha,
-                                        location, indices,
-                                        verbose=False,
-                                        use_cache=False):
-        """ Set values based on geo referenced geospatial data object.
-        """
+            raise Exception, 'Not implemented: %s' % location
+
+    ##
+    # @brief Set values based on geo referenced geospatial data object.
+    # @param geospatial_data ??
+    # @param alpha ??
+    # @param location ??
+    # @param indices ??
+    # @param verbose ??
+    # @param use_cache ??
+    def set_values_from_geospatial_data(self, geospatial_data,
+                                              alpha,
+                                              location,
+                                              indices,
+                                              verbose=False,
+                                              use_cache=False):
+        """Set values based on geo referenced geospatial data object."""
+
+        from anuga.coordinate_transforms.geo_reference import Geo_reference
 
         points = geospatial_data.get_data_points(absolute=False)
@@ -814,19 +788,11 @@
         data_georef = geospatial_data.get_geo_reference()
 
-
-        from anuga.coordinate_transforms.geo_reference import Geo_reference
-
-
         points = ensure_numeric(points, num.Float)
         values = ensure_numeric(values, num.Float)
 
         if location != 'vertices':
-            msg = 'set_values_from_points is only defined for '+\
-                  'location=\'vertices\''
-            raise ms
-
-        #coordinates = self.domain.get_nodes()
-        #triangles = self.domain.get_triangles()
-
+            msg = ("set_values_from_points is only defined for "
+                   "location='vertices'")
+            raise Exception, msg
 
         # Take care of georeferencing
@@ -834,10 +800,7 @@
             data_georef = Geo_reference()
 
-
         mesh_georef = self.domain.geo_reference
 
-
         # Call fit_interpolate.fit function
-        # args = (coordinates, triangles, points, values)
         args = (points, )
         kwargs = {'vertex_coordinates': None,
@@ -850,58 +813,72 @@
                   'verbose': verbose}
 
-        vertex_attributes = apply(fit_to_mesh,
-                                  args, kwargs)        
+        vertex_attributes = apply(fit_to_mesh, args, kwargs)
 
         # Call underlying method using array values
-        self.set_values_from_array(vertex_attributes,
-                                   location, indices,
-                                   use_cache=use_cache,
-                                   verbose=verbose)
-
-
-
-    def set_values_from_points(self, points, values, alpha,
-                               location, indices,
-                               data_georef=None,
-                               verbose=False,
-                               use_cache=False):
-        """
-        Set quantity values from arbitray data points using
-        fit_interpolate.fit
-        """
+        self.set_values_from_array(vertex_attributes, location,
+                                   indices, verbose)
+
+    ##
+    # @brief Set quantity values from arbitray data points.
+    # @param points ??
+    # @param values ??
+    # @param alpha ??
+    # @param location ??
+    # @param indices ??
+    # @param data_georef ??
+    # @param verbose True if this method is to be verbose.
+    # @param use_cache ??
+    def set_values_from_points(self, points,
+                                     values,
+                                     alpha,
+                                     location,
+                                     indices,
+                                     data_georef=None,
+                                     verbose=False,
+                                     use_cache=False):
+        """Set quantity values from arbitray data points using fit_interpolate.fit"""
 
         raise Exception, 'set_values_from_points is obsolete, use geospatial data object instead'
-        
-
-    def set_values_from_file(self, filename, attribute_name, alpha,
-                             location, indices,
-                             verbose=False,
-                             use_cache=False,
-                             max_read_lines=None):
-        """Set quantity based on arbitrary points in a points file
-        using attribute_name selects name of attribute
-        present in file.
+
+    ##
+    # @brief Set quantity based on arbitrary points in a points file.
+    # @param filename Path to the points file.
+    # @param attribute_name 
+    # @param alpha 
+    # @param location 
+    # @param indices 
+    # @param verbose True if this method is to be verbose.
+    # @param use_cache 
+    # @param max_read_lines 
+    def set_values_from_file(self, filename,
+                                   attribute_name,
+                                   alpha,
+                                   location,
+                                   indices,
+                                   verbose=False,
+                                   use_cache=False,
+                                   max_read_lines=None):
+        """Set quantity based on arbitrary points in a points file using
+        attribute_name selects name of attribute present in file.
         If attribute_name is not specified, use first available attribute
-	as defined in geospatial_data.
+        as defined in geospatial_data.
         """
 
         from types import StringType
+
         msg = 'Filename must be a text string'
         assert type(filename) == StringType, msg
 
-
         if location != 'vertices':
-            msg = 'set_values_from_file is only defined for '+\
-                  'location=\'vertices\''
-            raise msg
-
-        if True:  
+            msg = "set_values_from_file is only defined for location='vertices'"
+            raise Exception, msg
+
+        if True:
             # Use mesh as defined by domain
-            # This used to cause problems for caching 
-            # due to quantities changing, but
-            # it now works using the appropriate Mesh object.
+            # This used to cause problems for caching due to quantities
+            # changing, but it now works using the appropriate Mesh object.
             # This should close ticket:242
             vertex_attributes = fit_to_mesh(filename,
-                                            mesh=self.domain.mesh,  
+                                            mesh=self.domain.mesh,
                                             alpha=alpha,
                                             attribute_name=attribute_name,
@@ -911,10 +888,10 @@
         else:
             # This variant will cause Mesh object to be recreated
-            # in fit_to_mesh thus doubling up on the neighbour structure 
+            # in fit_to_mesh thus doubling up on the neighbour structure
             # FIXME(Ole): This is now obsolete 19 Jan 2009.
             nodes = self.domain.get_nodes(absolute=True)
-            triangles = self.domain.get_triangles()      
+            triangles = self.domain.get_triangles()
             vertex_attributes = fit_to_mesh(filename,
-                                            nodes, triangles, 
+                                            nodes, triangles,
                                             mesh=None,
                                             alpha=alpha,
@@ -923,16 +900,13 @@
                                             verbose=verbose,
                                             max_read_lines=max_read_lines)
-                                            
+
         # Call underlying method using array values
-        if verbose:
-            print 'Applying fitted data to domain'
-        self.set_values_from_array(vertex_attributes,
-                                   location, indices,
-                                   use_cache=use_cache,
-                                   verbose=verbose)
-
-    
-    
-    #-----------------------------------------------------    
+        self.set_values_from_array(vertex_attributes, location,
+                                   indices, verbose)
+
+    ##
+    # @brief Get index for maximum or minimum value of quantity.
+    # @param mode Either 'max' or 'min'.
+    # @param indices Set of IDs of elements to work on.
     def get_extremum_index(self, mode=None, indices=None):
         """Return index for maximum or minimum value of quantity (on centroids)
@@ -958,8 +932,9 @@
         if mode is None or mode == 'max':
             i = num.argmax(V)
-        elif mode == 'min':    
+        elif mode == 'min':
             i = num.argmin(V)
-
-            
+        else:
+            raise ValueError, 'Bad mode value, got: %s' % str(mode)
+
         if indices is None:
             return i
@@ -967,14 +942,15 @@
             return indices[i]
 
-
+    ##
+    # @brief Get index for maximum value of quantity.
+    # @param indices Set of IDs of elements to work on.
     def get_maximum_index(self, indices=None):
-        """See get extreme index for details
-        """
-
-        return self.get_extremum_index(mode='max',
-                                       indices=indices)
-
-
-        
+        """See get extreme index for details"""
+
+        return self.get_extremum_index(mode='max', indices=indices)
+
+    ##
+    # @brief Return maximum value of quantity (on centroids).
+    # @param indices Set of IDs of elements to work on.
     def get_maximum_value(self, indices=None):
         """Return maximum value of quantity (on centroids)
@@ -987,14 +963,15 @@
 
         Note, we do not seek the maximum at vertices as each vertex can
-        have multiple values - one for each triangle sharing it            
-        """
-
+        have multiple values - one for each triangle sharing it
+        """
 
         i = self.get_maximum_index(indices)
-        V = self.get_values(location='centroids') #, indices=indices)
-        
+        V = self.get_values(location='centroids')      #, indices=indices)
+
         return V[i]
-        
-
+
+    ##
+    # @brief Get location of maximum value of quantity (on centroids).
+    # @param indices Set of IDs of elements to work on.
     def get_maximum_location(self, indices=None):
         """Return location of maximum value of quantity (on centroids)
@@ -1005,5 +982,4 @@
         Usage:
             x, y = get_maximum_location()
-
 
         Notes:
@@ -1012,5 +988,5 @@
 
             If there are multiple cells with same maximum value, the
-            first cell encountered in the triangle array is returned.       
+            first cell encountered in the triangle array is returned.
         """
 
@@ -1020,13 +996,15 @@
         return x, y
 
-
+    ##
+    # @brief  Get index for minimum value of quantity.
+    # @param indices Set of IDs of elements to work on.
     def get_minimum_index(self, indices=None):
-        """See get extreme index for details
-        """        
-
-        return self.get_extremum_index(mode='min',
-                                       indices=indices)
-
-
+        """See get extreme index for details"""
+
+        return self.get_extremum_index(mode='min', indices=indices)
+
+    ##
+    # @brief Return minimum value of quantity (on centroids).
+    # @param indices Set of IDs of elements to work on.
     def get_minimum_value(self, indices=None):
         """Return minimum value of quantity (on centroids)
@@ -1038,14 +1016,16 @@
             v = get_minimum_value()
 
-        See get_maximum_value for more details.    
-        """
-
+        See get_maximum_value for more details.
+        """
 
         i = self.get_minimum_index(indices)
         V = self.get_values(location='centroids')
-        
+
         return V[i]
-        
-
+
+
+    ##
+    # @brief Get location of minimum value of quantity (on centroids).
+    # @param indices Set of IDs of elements to work on.
     def get_minimum_location(self, indices=None):
         """Return location of minimum value of quantity (on centroids)
@@ -1056,5 +1036,4 @@
         Usage:
             x, y = get_minimum_location()
-
 
         Notes:
@@ -1063,5 +1042,5 @@
 
             If there are multiple cells with same maximum value, the
-            first cell encountered in the triangle array is returned.       
+            first cell encountered in the triangle array is returned.
         """
 
@@ -1071,70 +1050,75 @@
         return x, y
 
-
-
+    ##
+    # @brief Get values at interpolation points.
+    # @param interpolation_points List of UTM coords or geospatial data object.
+    # @param use_cache ??
+    # @param verbose True if this method is to be verbose.
     def get_interpolated_values(self, interpolation_points,
-                                use_cache=False,
-                                verbose=False):
-        """ Get values at interpolation points
-        
-        The argument interpolation points must be given as either a 
+                                      use_cache=False,
+                                      verbose=False):
+        """Get values at interpolation points
+
+        The argument interpolation points must be given as either a
         list of absolute UTM coordinates or a geospatial data object.
         """
-        
 
         # FIXME (Ole): Points might be converted to coordinates relative to mesh origin
-        # This could all be refactored using the 
-        # 'change_points_geo_ref' method of Class geo_reference. 
+        # This could all be refactored using the
+        # 'change_points_geo_ref' method of Class geo_reference.
         # The purpose is to make interpolation points relative
         # to the mesh origin.
         #
         # Speed is also a consideration here.
-        
-        
-        # Ensure that interpolation points is either a list of 
+
+        # Ensure that interpolation points is either a list of
         # points, Nx2 array, or geospatial and convert to Numeric array
-        if isinstance(interpolation_points, Geospatial_data):        
+        if isinstance(interpolation_points, Geospatial_data):
             # Ensure interpolation points are in absolute UTM coordinates
-            interpolation_points = interpolation_points.get_data_points(absolute=True)
-                
+            interpolation_points = \
+                interpolation_points.get_data_points(absolute=True)
+
         # Reconcile interpolation points with georeference of domain
-        interpolation_points = self.domain.geo_reference.get_relative(interpolation_points) 
+        interpolation_points = \
+            self.domain.geo_reference.get_relative(interpolation_points)
         interpolation_points = ensure_numeric(interpolation_points)
 
-        
+
         # Get internal representation (disconnected) of vertex values
         vertex_values, triangles = self.get_vertex_values(xy=False,
-                                                          smooth=False)                
-    
+                                                          smooth=False)
+
         # Get possibly precomputed interpolation object
         I = self.domain.get_interpolation_object()
 
-        # Call interpolate method with interpolation points                
+        # Call interpolate method with interpolation points
         result = I.interpolate_block(vertex_values, interpolation_points,
-                                     use_cache=use_cache,
-                                     verbose=verbose)
-                               
+                                     use_cache=use_cache, verbose=verbose)
+
         return result
-        
-        
-
-
-    def get_values(self, 
-                   interpolation_points=None,
-                   location='vertices',
-                   indices=None,
-                   use_cache=False,
-                   verbose=False):
-        """get values for quantity
-
-        Extract values for quantity as a Numeric array.        
-        
+
+    ##
+    # @brief Get values as an array.
+    # @param interpolation_points List of coords to get values at.
+    # @param location Where to store results.
+    # @param indices Set of IDs of elements to work on.
+    # @param use_cache 
+    # @param verbose True if this method is to be verbose.
+    def get_values(self, interpolation_points=None,
+                         location='vertices',
+                         indices=None,
+                         use_cache=False,
+                         verbose=False):
+        """Get values for quantity
+
+        Extract values for quantity as a Numeric array.
+
         Inputs:
            interpolation_points: List of x, y coordinates where value is
-                                 sought (using interpolation). If points 
-                                 are given, values of location and indices 
+                                 sought (using interpolation). If points
+                                 are given, values of location and indices
                                  are ignored. Assume either absolute UTM
                                  coordinates or geospatial data object.
-        
+
            location: Where values are to be stored.
                      Permissible options are: vertices, edges, centroids
@@ -1148,5 +1132,5 @@
         values will be a list or a Numerical array of length N, N being
         the number of elements.
-        
+
         In case of location == 'vertices' or 'edges' the dimension of
         returned values will be of dimension Nx3
@@ -1154,6 +1138,6 @@
         In case of location == 'unique vertices' the average value at
         each vertex will be returned and the dimension of returned values
-        will be a 1d array of length "number of vertices" 
-        
+        will be a 1d array of length "number of vertices"
+
         Indices is the set of element ids that the operation applies to.
 
@@ -1161,12 +1145,11 @@
         internal ordering.
         """
-        
 
         # FIXME (Ole): I reckon we should have the option of passing a
         #              polygon into get_values. The question becomes how
         #              resulting values should be ordered.
-        
+
         if verbose is True:
-            print 'Getting values from %s' %location
+            print 'Getting values from %s' % location
 
         if interpolation_points is not None:
@@ -1174,29 +1157,29 @@
                                                 use_cache=use_cache,
                                                 verbose=verbose)
-        
-        
+
         # FIXME (Ole): Consider deprecating 'edges' - but not if it is used
-        # elsewhere in ANUGA. 
+        # elsewhere in ANUGA.
         # Edges have already been deprecated in set_values, see changeset:5521,
         # but *might* be useful in get_values. Any thoughts anyone?
-        
-        if location not in ['vertices', 'centroids', 'edges',
-                            'unique vertices']:
-            msg = 'Invalid location: %s' %location
+
+        if location not in ['vertices', 'centroids',
+                            'edges', 'unique vertices']:
+            msg = 'Invalid location: %s' % location
             raise msg
 
         import types
+
         assert type(indices) in [types.ListType, types.NoneType,
-                                 num.ArrayType],\
-                                 'Indices must be a list or None'
+                                 num.ArrayType], \
+                   'Indices must be a list or None'
 
         if location == 'centroids':
             if (indices ==  None):
                 indices = range(len(self))
-            return num.take(self.centroid_values,indices)
+            return num.take(self.centroid_values, indices)
         elif location == 'edges':
             if (indices ==  None):
                 indices = range(len(self))
-            return num.take(self.edge_values,indices)
+            return num.take(self.edge_values, indices)
         elif location == 'unique vertices':
             if (indices ==  None):
@@ -1207,18 +1190,17 @@
             for unique_vert_id in indices:
                 triangles = self.domain.get_triangles_and_vertices_per_node(node=unique_vert_id)
-                    
+
                 # In case there are unused points
                 if len(triangles) == 0:
                     msg = 'Unique vertex not associated with triangles'
-                    raise msg
+                    raise Exception, msg
 
                 # Go through all triangle, vertex pairs
                 # Average the values
-                
                 # FIXME (Ole): Should we merge this with get_vertex_values
                 sum = 0
                 for triangle_id, vertex_id in triangles:
                     sum += self.vertex_values[triangle_id, vertex_id]
-                vert_values.append(sum/len(triangles))
+                vert_values.append(sum / len(triangles))
             return num.array(vert_values, num.Float)
         else:
@@ -1227,13 +1209,12 @@
             return num.take(self.vertex_values, indices)
 
-
-
-    def set_vertex_values(self, A,
-                          indices=None,
-                          use_cache=False,
-                          verbose=False):
+    ##
+    # @brief Set vertex values for all unique vertices based on array.
+    # @param A Array to set values with.
+    # @param indices Set of IDs of elements to work on.
+    def set_vertex_values(self, A, indices=None):
         """Set vertex values for all unique vertices based on input array A
-        which has one entry per unique vertex, i.e.
-        one value for each row in array self.domain.nodes.
+        which has one entry per unique vertex, i.e. one value for each row in
+        array self.domain.nodes.
 
         indices is the list of vertex_id's that will be set.
@@ -1242,9 +1223,6 @@
         """
 
-
-        # Assert that A can be converted to a Numeric array of appropriate dim
+        # Check that A can be converted to array and is of appropriate dim
         A = ensure_numeric(A, num.Float)
-
-        # print 'SHAPE A', A.shape
         assert len(A.shape) == 1
 
@@ -1256,28 +1234,11 @@
             vertex_list = indices
 
-
-        #FIXME(Ole): This function ought to be faster.
-        # We need to get the triangles_and_vertices list
-        # from domain in one hit, then cache the computation of the
-        # Nx3 array of vertex values that can then be assigned using
-        # set_values_from_array.
-        #
-        # Alternatively, some C code would be handy
-        #
-        self._set_vertex_values(vertex_list, A)
-            
-
-    def _set_vertex_values(self, vertex_list, A):
-        """Go through list of unique vertices
-        This is the common case e.g. when values
-        are obtained from a pts file through fitting
-        """
-        
+        # Go through list of unique vertices
         for i_index, unique_vert_id in enumerate(vertex_list):
-
             triangles = self.domain.get_triangles_and_vertices_per_node(node=unique_vert_id)
-                    
+
             # In case there are unused points
-            if len(triangles) == 0: continue
+            if len(triangles) == 0:
+                continue
 
             # Go through all triangle, vertex pairs
@@ -1289,23 +1250,24 @@
         self.interpolate()
 
-
-    def smooth_vertex_values(self,
-                             use_cache=False,
-                             verbose=False):
-        """ Smooths vertex values.
-        """
-
-        A,V = self.get_vertex_values(xy=False, smooth=True)
-        self.set_vertex_values(A,
-                               use_cache=use_cache,
-                               verbose=verbose)                               
-                               
-
-
+    ##
+    # @brief Smooth vertex values.
+    def smooth_vertex_values(self):
+        """Smooths vertex values."""
+
+        A, V = self.get_vertex_values(xy=False, smooth=True)
+        self.set_vertex_values(A)
+
+    ############################################################################
     # Methods for outputting model results
-    def get_vertex_values(self,
-                          xy=True,
-                          smooth=None,
-                          precision=None):
+    ############################################################################
+
+    ##
+    # @brief Get vertex values like an OBJ format i.e. one value per node.
+    # @param xy True if we return X and Y as well as A and V.
+    # @param smooth True if vertex values are to be smoothed.
+    # @param precision The type of the result values (default float).
+    def get_vertex_values(self, xy=True,
+                                smooth=None,
+                                precision=None):
         """Return vertex values like an OBJ format i.e. one value per node.
 
@@ -1316,27 +1278,23 @@
         Mx3, where M is the number of triangles. Each row has three indices
         defining the triangle and they correspond to elements in the arrays
-        X, Y and A. 
-
-        if smooth is True, vertex values corresponding to one common
-        coordinate set will be smoothed by taking the average of vertex values for each node.
-        In this case vertex coordinates will be
-        de-duplicated corresponding to the original nodes as obtained from
-        the method general_mesh.get_nodes()
-
-        If no smoothings is required, vertex coordinates and values will
-        be aggregated as a concatenation of values at
-        vertices 0, vertices 1 and vertices 2. This corresponds to
-        the node coordinates obtained from the method
-        general_mesh.get_vertex_coordinates()
-
+        X, Y and A.
+
+        If smooth is True, vertex values corresponding to one common coordinate
+        set will be smoothed by taking the average of vertex values for each
+        node.  In this case vertex coordinates will be de-duplicated
+        corresponding to the original nodes as obtained from the method
+        general_mesh.get_nodes()
+
+        If no smoothings is required, vertex coordinates and values will be
+        aggregated as a concatenation of values at vertices 0, vertices 1 and
+        vertices 2.  This corresponds to the node coordinates obtained from the
+        method general_mesh.get_vertex_coordinates()
 
         Calling convention
         if xy is True:
-           X,Y,A,V = get_vertex_values
+           X, Y, A, V = get_vertex_values
         else:
-           A,V = get_vertex_values
-
-        """
-
+           A, V = get_vertex_values
+        """
 
         if smooth is None:
@@ -1349,15 +1307,12 @@
         if precision is None:
             precision = num.Float
-            
 
         if smooth is True:
-            # Ensure continuous vertex values by averaging
-            # values at each node
-            
+            # Ensure continuous vertex values by averaging values at each node
             V = self.domain.get_triangles()
             N = self.domain.number_of_full_nodes # Ignore ghost nodes if any
             A = num.zeros(N, num.Float)
-            points = self.domain.get_nodes()            
-            
+            points = self.domain.get_nodes()
+
             if 1:
                 # Fast C version
@@ -1367,8 +1322,6 @@
                                       A)
                 A = A.astype(precision)
-            else:    
-
+            else:
                 # Slow Python version
-                
                 current_node = 0
                 k = 0 # Track triangles touching on node
@@ -1376,53 +1329,42 @@
                 for index in self.domain.vertex_value_indices:
                     if current_node == N:
-                        msg = 'Current node exceeding number of nodes (%d) ' %(N)
+                        msg = 'Current node exceeding number of nodes (%d) ' % N
                         raise msg
-                    
-
-                    
+
                     k += 1
-                    
+
                     volume_id = index / 3
                     vertex_id = index % 3
-                 
-                    #assert V[volume_id, vertex_id] == current_node
-                
+
                     v = self.vertex_values[volume_id, vertex_id]
                     total += v
 
-                    #print 'current_node=%d, index=%d, k=%d, total=%f' %(current_node, index, k, total)
                     if self.domain.number_of_triangles_per_node[current_node] == k:
                         A[current_node] = total/k
-                
-                    
+
                         # Move on to next node
                         total = 0.0
                         k = 0
                         current_node += 1
-
-
-
         else:
-            # Return disconnected internal vertex values 
+            # Return disconnected internal vertex values
             V = self.domain.get_disconnected_triangles()
             points = self.domain.get_vertex_coordinates()
             A = self.vertex_values.flat.astype(precision)
 
-
-        # Return    
+        # Return
         if xy is True:
             X = points[:,0].astype(precision)
             Y = points[:,1].astype(precision)
-            
+
             return X, Y, A, V
         else:
-            return A, V            
-
-
-
+            return A, V
+
+    ##
+    # @brief Extrapolate conserved quantities from centroid.
     def extrapolate_first_order(self):
-        """Extrapolate conserved quantities from centroid to
-        vertices and edges for each volume using
-        first order scheme.
+        """Extrapolate conserved quantities from centroid to vertices and edges
+        for each volume using first order scheme.
         """
 
@@ -1438,8 +1380,10 @@
         self.y_gradient *= 0.0
 
-
+    ##
+    # @brief Compute the integral of quantity across entire domain.
+    # @return The integral.
     def get_integral(self):
-        """Compute the integral of quantity across entire domain
-        """
+        """Compute the integral of quantity across entire domain."""
+
         areas = self.domain.get_areas()
         integral = 0
@@ -1451,11 +1395,14 @@
         return integral
 
+    ##
+    # @brief get the gradients.
     def get_gradients(self):
-        """Provide gradients. Use compute_gradients first
-        """
+        """Provide gradients. Use compute_gradients first."""
 
         return self.x_gradient, self.y_gradient
 
-
+    ##
+    # @brief ??
+    # @param timestep ??
     def update(self, timestep):
         # Call correct module function
@@ -1463,4 +1410,6 @@
         return update(self, timestep)
 
+    ##
+    # @brief ??
     def compute_gradients(self):
         # Call correct module function
@@ -1468,4 +1417,6 @@
         return compute_gradients(self)
 
+    ##
+    # @brief ??
     def limit(self):
         # Call correct module depending on whether
@@ -1473,4 +1424,6 @@
         limit_old(self)
 
+    ##
+    # @brief ??
     def limit_vertices_by_all_neighbours(self):
         # Call correct module function
@@ -1478,4 +1431,6 @@
         limit_vertices_by_all_neighbours(self)
 
+    ##
+    # @brief ??
     def limit_edges_by_all_neighbours(self):
         # Call correct module function
@@ -1483,9 +1438,13 @@
         limit_edges_by_all_neighbours(self)
 
+    ##
+    # @brief ??
     def limit_edges_by_neighbour(self):
         # Call correct module function
         # (either from this module or C-extension)
-        limit_edges_by_neighbour(self)               
-
+        limit_edges_by_neighbour(self)
+
+    ##
+    # @brief ??
     def extrapolate_second_order(self):
         # Call correct module function
@@ -1493,5 +1452,7 @@
         compute_gradients(self)
         extrapolate_from_gradient(self)
-        
+
+    ##
+    # @brief ??
     def extrapolate_second_order_and_limit_by_edge(self):
         # Call correct module function
@@ -1499,4 +1460,6 @@
         extrapolate_second_order_and_limit_by_edge(self)
 
+    ##
+    # @brief ??
     def extrapolate_second_order_and_limit_by_vertex(self):
         # Call correct module function
@@ -1504,4 +1467,7 @@
         extrapolate_second_order_and_limit_by_vertex(self)
 
+    ##
+    # @brief ??
+    # @param bound ??
     def bound_vertices_below_by_constant(self, bound):
         # Call correct module function
@@ -1509,4 +1475,7 @@
         bound_vertices_below_by_constant(self, bound)
 
+    ##
+    # @brief ??
+    # @param quantity ??
     def bound_vertices_below_by_quantity(self, quantity):
         # Call correct module function
@@ -1515,6 +1484,8 @@
         # check consistency
         assert self.domain == quantity.domain
-        bound_vertices_below_by_quantity(self, quantity)                        
-
+        bound_vertices_below_by_quantity(self, quantity)
+
+    ##
+    # @brief ??
     def backup_centroid_values(self):
         # Call correct module function
@@ -1522,28 +1493,32 @@
         backup_centroid_values(self)
 
-    def saxpy_centroid_values(self,a,b):
+    ##
+    # @brief ??
+    # @param a ??
+    # @param b ??
+    def saxpy_centroid_values(self, a, b):
         # Call correct module function
         # (either from this module or C-extension)
-        saxpy_centroid_values(self,a,b)
-    
-#Conserved_quantity = Quantity
-
+        saxpy_centroid_values(self, a, b)
+
+
+##
+# @brief OBSOLETE!
 class Conserved_quantity(Quantity):
-    """Class conserved quantity being removed, use Quantity
-
-    """
+    """Class conserved quantity being removed, use Quantity."""
 
     def __init__(self, domain, vertex_values=None):
-        #Quantity.__init__(self, domain, vertex_values)
-
         msg = 'ERROR: Use Quantity instead of Conserved_quantity'
-
         raise Exception, msg
 
 
+######
+# Prepare the C extensions.
+######
 
 from anuga.utilities import compile
-if compile.can_use_C_extension('quantity_ext.c'):    
-    # Underlying C implementations can be accessed 
+
+if compile.can_use_C_extension('quantity_ext.c'):
+    # Underlying C implementations can be accessed
 
     from quantity_ext import \
@@ -1564,11 +1539,8 @@
          interpolate_from_vertices_to_edges,\
          interpolate_from_edges_to_vertices,\
-         update    
+         update
 else:
-    msg = 'C implementations could not be accessed by %s.\n ' %__file__
+    msg = 'C implementations could not be accessed by %s.\n ' % __file__
     msg += 'Make sure compile_all.py has been run as described in '
     msg += 'the ANUGA installation guide.'
     raise Exception, msg
-
-
-
