Index: /anuga_core/documentation/user_manual/anuga_user_manual.tex
===================================================================
--- /anuga_core/documentation/user_manual/anuga_user_manual.tex	(revision 7133)
+++ /anuga_core/documentation/user_manual/anuga_user_manual.tex	(revision 7134)
@@ -216,5 +216,5 @@
 regions using symbolic tags.  These tags may then be used to set the
 actual boundary conditions and attributes for different regions
-(e.g.\ the Manning friction coefficient) for each simulation.
+(e.g. the Manning friction coefficient) for each simulation.
 
 Most \anuga components are written in the object-oriented programming
@@ -237,5 +237,5 @@
   \item The mathematical model is the 2D shallow water wave equation.
   As such it cannot resolve vertical convection and consequently not breaking
-  waves or 3D turbulence (e.g.\ vorticity).
+  waves or 3D turbulence (e.g. vorticity).
   %\item The surface is assumed to be open, e.g. \anuga cannot model
   %flow under ceilings or in pipes
@@ -251,5 +251,5 @@
   \item Frictional resistance is implemented using Manning's formula, but
   \anuga has not yet been fully validated in regard to bottom roughness.
-%\item ANUGA contains no tsunami-genic functionality relating to earthquakes.
+%\item \anuga contains no tsunami-genic functionality relating to earthquakes.
 \end{itemize}
 
@@ -523,5 +523,5 @@
           not steady state as replication of momentum at the boundary
           may cause numerical instabilities propagating into the domain and 
-          eventually causing ANUGA to crash. If this occurs,
+          eventually causing \anuga to crash. If this occurs,
           consider using e.g. a Dirichlet boundary condition with a stage value 
           less than the elevation at the boundary.
@@ -574,5 +574,5 @@
 uses different boundary tags -- in general, the possible tags are entirely selectable by the user when generating the mesh and not
 limited to 'left', 'right', 'top' and 'bottom' as in this example.)
-All segments in bounding polygon must be tagged. If a tag is not supplied, the default tag name 'exterior' will be assigned by ANUGA.
+All segments in bounding polygon must be tagged. If a tag is not supplied, the default tag name 'exterior' will be assigned by \anuga.
 
 Using the boundary objects described above, we assign a boundary
@@ -629,5 +629,5 @@
 The output is a NetCDF file with the extension \code{.sww}. It
 contains stage and momentum information and can be used with the
-ANUGA viewer \code{animate} to generate a visual
+\anuga viewer \code{animate} to generate a visual
 display (see Section \ref{sec:animate}). See Section \ref{sec:file formats}
 (page \pageref{sec:file formats}) for more on NetCDF and other file
@@ -666,5 +666,5 @@
   \centerline{\includegraphics[width=75mm, height=75mm]
     {graphics/bedslopestart.jpg}}
-  \caption{Runup example viewed with the ANUGA viewer}
+  \caption{Runup example viewed with the \anuga viewer}
   \label{fig:runupstart}
 \end{figure}
@@ -727,5 +727,5 @@
 \end{verbatim}
 
-returns an \code{m x n} mesh similar to the one used in the previous example, except that now the
+returns an \code{mxn} mesh similar to the one used in the previous example, except that now the
 extent in the x and y directions are given by the value of \code{length} and \code{width}
 respectively.
@@ -763,5 +763,5 @@
   \centerline{\includegraphics[height=75mm]
     {graphics/channel1.png}}%
-  \caption{Simple channel example viewed with the ANUGA viewer.}
+  \caption{Simple channel example viewed with the \anuga viewer.}
   \label{fig:channel1}
 \end{figure}
@@ -928,5 +928,6 @@
 maximal area of a triangle used for triangulation -- and a triangular
 mesh is created inside the polygon using a mesh generation engine.
-On any given platform, the same mesh will be returned.
+On any given platform, the same mesh will be returned each time the
+script is run.
 
 Boundary tags are not restricted to \code{'left'}, \code{'bottom'},
@@ -1056,12 +1057,13 @@
 
 Since we used \code{create_domain_from_regions} to create the mesh file, we do not need to 
-create the domain explicitly, as the above function also does that.
+create the domain explicitly, as the above function does both mesh and domain creation.
 
 The following statements specify a basename and data directory, and
-sets a minimum storable height, which helps with visualisation.
-
-\begin{verbatim}
-domain.set_name('cairns_' + project.scenario) # Name of sww file
-domain.set_datadir('.')                       # Store sww output here
+sets a minimum storable height, which helps with visualisation and post-processing
+if one wants to remove water less than 1cm deep (for instance).
+
+\begin{verbatim}
+domain.set_name('cairns_' + project.scenario) # Name of SWW file
+domain.set_datadir('.')                       # Store SWW output here
 domain.set_minimum_storable_height(0.01)      # Store only depth > 1cm
 \end{verbatim}
@@ -1073,9 +1075,9 @@
 in this case, many of the values are read from the auxiliary file
 \file{project.py} or, in the case of \code{elevation}, from an
-ancillary points file.
+auxiliary points file.
 
 \subsubsection{Stage}
 
-The stage is initially set to 0.0 by the following statements:
+The stage is initially set to 0.0 (i.e. Mean Sea Level) by the following statements:
 
 \begin{verbatim}
@@ -1083,4 +1085,6 @@
 domain.set_quantity('stage', tide)
 \end{verbatim}
+
+It could also take the value of the highest astronomical tide.
 
 %For the scenario we are modelling in this case, we use a callable
@@ -1109,5 +1113,6 @@
 \subsubsection{Elevation}
 
-The elevation is specified by reading data from a file:
+The elevation is specified by reading data from a file with a name derived from
+\code{project.demname} with the \code{.pts} extension:
 
 \begin{verbatim}
@@ -1118,4 +1123,13 @@
                     alpha=0.1)
 \end{verbatim}
+
+The \code{alpha} parameter controls how smooth the elevation surface
+should be.  See section \ref{class:alpha_shape}, page \pageref{class:alpha_shape}.
+
+Setting \code{cache=True} allows \anuga to save the result in order
+to make subsequent runs faster.
+
+Using \code{verbose=True} tells the function to write diagnostics to
+the screen.
 
 \subsection{Boundary Conditions}\index{boundary conditions}
@@ -1165,15 +1179,14 @@
 \begin{verbatim}
 if project.scenario == 'slide':
+    # Initial run without any event
     for t in domain.evolve(yieldstep=10, finaltime=60):
         print domain.timestepping_statistics()
         print domain.boundary_statistics(tags='ocean_east')
 
-    # Add slide
-    thisstagestep = domain.get_quantity('stage')
+    # Add slide to water surface
     if allclose(t, 60):
-        slide = Quantity(domain)
-        slide.set_values(tsunami_source)
-        domain.set_quantity('stage', slide + thisstagestep)
-
+        domain.add_quantity('stage', tsunami_source)
+
+    # Continue propagating wave
     for t in domain.evolve(yieldstep=10, finaltime=5000,
                            skip_initial_step=True):
@@ -1464,12 +1477,16 @@
 \label{sec:meshgeneration}
 
-\begin{funcdesc} {create\_mesh\_from\_regions}{bounding_polygon,
-                             boundary_tags,
-                             maximum_triangle_area,
-                             filename=None,
-                             interior_regions=None,
-                             poly_geo_reference=None,
-                             mesh_geo_reference=None,
-                             minimum_triangle_angle=28.0}
+\begin{funcdesc}{create_mesh_from_regions}{bounding_polygon,
+                                           boundary_tags,
+                                           maximum_triangle_area=None,
+                                           filename=None,
+                                           interior_regions=None,
+                                           interior_holes=None,
+                                           poly_geo_reference=None,
+                                           mesh_geo_reference=None,
+                                           minimum_triangle_angle=28.0,
+                                           fail_if_polygons_outside=True,
+                                           use_cache=False,
+                                           verbose=True}
 Module: \module{pmesh.mesh\_interface}
 
@@ -1503,4 +1520,78 @@
 \end{funcdesc}
 
+\begin{funcdesc}{create_domain_from_regions}{bounding_polygon,
+                                             boundary_tags,
+                                             maximum_triangle_area=None,
+                                             mesh_filename=None,
+                                             interior_regions=None,
+                                             interior_holes=None,
+                                             poly_geo_reference=None,
+                                             mesh_geo_reference=None,
+                                             minimum_triangle_angle=28.0,
+                                             fail_if_polygons_outside=True,
+                                             use_cache=False,
+                                             verbose=True}
+
+Module: \module{interface.py}
+
+This higher-level function allows a user to create a domain (and associated mesh)
+inside a specified polygon.
+
+\code{bounding_polygon} is a list of points in Eastings and Northings,
+relative to the zone stated in \code{poly_geo_reference} if specified.
+Otherwise points are just x, y coordinates with no particular
+association to any location.
+
+\code{boundary_tags} is a dictionary of symbolic tags. For every tag there
+is a list of indices referring to segments associated with that tag.
+If a segment is omitted it will be assigned the default tag ''.
+
+\code{maximum_triangle_area} is the maximal area per triangle
+for the bounding polygon, excluding the interior regions.
+
+\code{mesh_filename} is the name of the file to contain the generated 
+mesh data.
+
+\code{interior_regions} is a list of tuples consisting of (polygon,
+resolution) for each region to be separately refined. Do not have
+polygon lines cross or be on-top of each other.  Also do not have
+polygons close to each other.
+
+\code{poly_geo_reference} is the geo_reference of the bounding polygon and
+the interior polygons.
+If none, assume absolute.  Please pass one though, since absolute
+references have a zone.
+
+\code{mesh_geo_reference} is the geo_reference of the mesh to be created.
+If none is given one will be automatically generated.  It will use
+the lower left hand corner of  bounding_polygon (absolute)
+as the x and y values for the geo_ref.
+
+\code{minimum_triangle_angle} is the minimum angle allowed for each generated triangle.
+This controls the \emph{slimness} allowed for a triangle.
+
+\code{fail_if_polygons_outside} -- if True (the default) an Exception in thrown
+if interior polygons fall outside the bounding polygon. If False, these
+will be ignored and execution continues.
+
+\textbf{WARNING}. Note that the dictionary structure used for the
+parameter \code{boundary_tags} is different from that used for the
+variable \code{boundary} that occurs in the specification of a mesh.
+In the case of \code{boundary}, the tags are the \emph{values} of
+the dictionary, whereas in the case of \code{boundary_tags}, the
+tags are the \emph{keys} and the \emph{value} corresponding to a
+particular tag is a list of numbers identifying boundary edges
+labelled with that tag. Because of this, it is theoretically
+possible to assign the same edge to more than one tag. However, an
+attempt to do this will cause an error.
+
+\textbf{WARNING}. Do not have polygon lines cross or be on-top of each
+    other. This can result in regions of unspecified resolutions. Do
+    not have polygon close to each other. This can result in the area
+    between the polygons having small triangles.  For more control
+    over the mesh outline use the methods described below.
+
+\end{funcdesc}
+
 \subsection{Advanced mesh generation}
 
@@ -1508,30 +1599,41 @@
 methods of the class \class{Mesh}.
 
-
-\begin{classdesc}  {Mesh}{userSegments=None,
-                 userVertices=None,
-                 holes=None,
-                 regions=None}
+\begin{classdesc}{Mesh}{userSegments=None,
+                        userVertices=None,
+                        holes=None,
+                        regions=None,
+                        geo_reference=None}
 Module: \module{pmesh.mesh}
 
 A class used to build a mesh outline and generate a two-dimensional
 triangular mesh. The mesh outline is used to describe features on the
-mesh, such as the mesh boundary. Many of this classes methods are used
-to build a mesh outline, such as \code{add\_vertices} and
-\code{add\_region\_from\_polygon}.
+mesh, such as the mesh boundary. Many of this class's methods are used
+to build a mesh outline, such as \code{add_vertices()} and
+\code{add_region_from_polygon()}.
+
+\code{userSegments} and \code{userVertices} define the outline enclosing the mesh.
+
+\code{holes} describes any regions inside the mesh that are not to be included in the mesh.
+
+\code{geo_reference} defines the geo_reference to which all point information is relative.
+If \code{geo_reference} is \code{None} then the default geo_reference is used.
 \end{classdesc}
 
 \subsubsection{Key Methods of Class Mesh}
 
-\begin{methoddesc} {add\_hole}{self, x, y}
-Module: \module{pmesh.mesh},  Class: \class{Mesh}
-
-This method is used to build the mesh outline.  It defines a hole,
-when the boundary of the hole has already been defined, by selecting a
-point within the boundary.
-\end{methoddesc}
-
-\begin{methoddesc}  {add\_hole\_from\_polygon}{self, polygon, tags=None}
-Module: \module{pmesh.mesh},  Class: \class{Mesh}
+\begin{methoddesc}{\emph{<mesh>}.add_hole}{x, y, geo_reference=None}
+Module: \module{pmesh.mesh}
+
+This method adds a hole to the mesh outline.
+
+\code{x} and \code{y} define a point on the already defined hole boundary.
+
+If \code{geo_reference} is not supplied the points are assumed to be absolute.
+\end{methoddesc}
+
+\begin{methoddesc}{\emph{<mesh>}.add_hole_from_polygon}{polygon,
+                                                        segment_tags=None,
+                                                        geo_reference=None}
+Module: \module{pmesh.mesh}
 
 This method is used to add a 'hole' within a region -- that is, to
@@ -1540,52 +1642,79 @@
 the polygon passed in.  Additionally, the user specifies a list of
 boundary tags, one for each edge of the bounding polygon.
-\end{methoddesc}
-
-\begin{methoddesc}  {add\_points_and_segments}{self, points, segments,
-    segment\_tags=None}
-Module: \module{pmesh.mesh},  Class: \class{Mesh}
-
-This method is used to build the mesh outline. It adds points and
-segments connecting the points.  Points is a list of points. Segments
-is a list of segments.  Each segment is defined by the start and end
-of the line by it's point index, e.g. use \code{segments =
-[[0,1],[1,2]]} to make a polyline between points 0, 1 and 2. A tag for
-each segment can optionally be added.
-\end{methoddesc}
-
-\begin{methoddesc} {add\_region}{self, x, y}
-Module: \module{pmesh.mesh},  Class: \class{Mesh}
-
-This method is used to build the mesh outline.  It defines a region,
-when the boundary of the region has already been defined, by selecting
-a point within the boundary.  A region instance is returned.  This can
-be used to set the resolution.
-\end{methoddesc}
-
-\begin{methoddesc}  {add\_region\_from\_polygon}{self, polygon,
-                     segment_tags=None, region_tag=None, max_triangle_area=None}
-Module: \module{pmesh.mesh},  Class: \class{Mesh}
-
-This method is used to build the mesh outline.  It adds a region to a
+
+\code{polygon} is the polygon that defines the hole to be added to the \code{<mesh>}.
+
+\code{segment_tags} -- ??
+
+If \code{geo_reference} is \code{None} then the default \code{geo_reference} is used.
+\end{methoddesc}
+
+\begin{methoddesc}{\emph{<mesh>}.add_points_and_segments}{points,
+                                                          segments=None,
+                                                          segment_tags=None}
+Module: \module{pmesh.mesh}
+
+This adds points and segments connecting the points to a mesh.
+
+\code{points} is a list of points.
+
+\code{segments} is a list of segments.  Each segment is defined by the start and end
+of the line by its point index, e.g. use \code{segments = [[0,1],[1,2]]} to make a
+polyline between points 0, 1 and 2.
+
+\code{segment_tags} may be used to optionally define a tag for each segment.
+\end{methoddesc}
+
+\begin{methoddesc}{\emph{<mesh>}.add_region}{x,y, geo_reference=None, tag=None}
+Module: \module{pmesh.mesh}
+
+This method adds a region to a mesh outline.
+
+\code{x} and \code{y} define a point on the already-defined region that is to 
+be added to the mesh.
+
+If \code{geo_reference} is not supplied the points data is assumed to be absolute.
+
+\code{tag} -- ??
+
+A region instance is returned.  This can be used to set the resolution of the added region.
+\end{methoddesc}
+
+\begin{methoddesc}{\emph{<mesh>}.add_region_from_polygon}{polygon,
+                                                          segment_tags=None,
+                                                          max_triangle_area=None,
+                                                          geo_reference=None,
+                                                          region_tag=None}
+Module: \module{pmesh.mesh}
+
+This method adds a region to a
 \class{Mesh} instance.  Regions are commonly used to describe an area
-with an increased density of triangles, by setting
-\code{max_triangle_area}.  The
-region boundary is described by the input \code{polygon}.  Additionally, the
-user specifies a list of segment tags, one for each edge of the
-bounding polygon.  The regional tag is set using  \code{region}.
-\end{methoddesc}
-
-\begin{methoddesc} {add\_vertices}{self, point_data}
-Module: \module{pmesh.mesh},  Class: \class{Mesh}
-
-Add user vertices. The point_data can be a list of (x,y) values, a numeric
-array or a geospatial_data instance.
-\end{methoddesc}
-
-\begin{methoddesc} {auto\_segment}{self, raw_boundary=raw_boundary,
-                    remove_holes=remove_holes,
-                    smooth_indents=smooth_indents,
-                    expand_pinch=expand_pinch}
-Module: \module{pmesh.mesh},  Class: \class{Mesh}
+with an increased density of triangles by setting \code{max_triangle_area}.
+
+\code{polygon} describes the region boundary to add to the \code{<mesh>}.
+
+\code{segment_tags} specifies a list of segment tags, one for each edge of the
+bounding polygon.
+
+If \code{geo_reference} is not supplied the points data is assumed to be absolute.
+
+\code{region_tag} sets the region tag.
+\end{methoddesc}
+
+\begin{methoddesc}{\emph{<mesh>}.add_vertices}{point_data}
+Module: \module{pmesh.mesh}
+
+Add user vertices to a mesh.
+
+\code{point_data} is the list of point data, and can be a list of (x,y) values,
+a numeric array or a geospatial_data instance.
+\end{methoddesc}
+
+\begin{methoddesc}{\emph{<mesh>}.auto_segment}{alpha=None,
+                                               raw_boundary=True,
+                                               remove_holes=False,
+                                               smooth_indents=False,
+                                               expand_pinch=False}
+Module: \module{pmesh.mesh}
 
 Add segments between some of the user vertices to give the vertices an
@@ -1593,41 +1722,63 @@
 useful since a set of user vertices need to be outlined by segments
 before generate_mesh is called.
-\end{methoddesc}
-
-\begin{methoddesc}  {export\_mesh_file}{self, ofile}
-Module: \module{pmesh.mesh},  Class: \class{Mesh}
-
-This method is used to save the mesh to a file. \code{ofile} is the
-name of the mesh file to be written, including the extension.  Use
-the extension \code{.msh} for the file to be in NetCDF format and
+
+\code{alpha} determines the $smoothness$ of the alpha shape.
+
+\code{raw_boundary}, if \code{True} instructs the function to return the raw
+boundary, i.e. the regular edges of the alpha shape.
+
+\code{remove_holes}, if \code{True} enables a filter to remove small holes
+(small is defined by  boundary_points_fraction).
+
+\code{smooth_indents}, if \code{True} removes sharp triangular indents
+in the boundary.
+
+\code{expand_pinch}, if \code{True} tests for pinch-off and corrects
+(i.e. a boundary vertex with more than two edges).
+\end{methoddesc}
+
+\begin{methoddesc}{\emph{<mesh>}.export_mesh_file}{ofile}
+Module: \module{pmesh.mesh}
+
+This method is used to save a mesh to a file.
+
+\code{ofile} is the name of the mesh file to be written, including the extension.
+Use the extension \code{.msh} for the file to be in NetCDF format and
 \code{.tsh} for the file to be ASCII format.
 \end{methoddesc}
 
-\begin{methoddesc}  {generate\_mesh}{self,
-                      maximum_triangle_area=None,
-                      minimum_triangle_angle=28.0,
-                      verbose=False}
-Module: \module{pmesh.mesh},  Class: \class{Mesh}
-
-This method is used to generate the triangular mesh.  The  maximal
-area of any triangle in the mesh can be specified, which is used to
-control the triangle density, along with the
-minimum angle in any triangle.
-\end{methoddesc}
-
-\begin{methoddesc}  {import_ungenerate_file}{self, ofile, tag=None,
-region_tag=None}
-Module: \module{pmesh.mesh},  Class: \class{Mesh}
+\begin{methoddesc}{\emph{<mesh>}.generate_mesh}{maximum_triangle_area="",
+                                                minimum_triangle_angle=28.0,
+                                                verbose=True}
+Module: \module{pmesh.mesh}
+
+This method is used to generate the triangular mesh.
+
+\code{maximum_triangle_area} sets the maximum area of any triangle in the mesh.
+
+\code{minimum_triangle_angle} sets the minimum area of any triangle in the mesh.
+
+These two parameters can be used to control the triangle density.
+\end{methoddesc}
+
+\begin{methoddesc}{\emph{<mesh>}.import_ungenerate_file}{ofile,
+                                                         tag=None,
+                                                         region_tag=None}
+Module: \module{pmesh.mesh}
 
 This method is used to import a polygon file in the ungenerate format,
 which is used by arcGIS. The polygons from the file are converted to
-vertices and segments. \code{ofile} is the name of the polygon file.
+vertices and segments.
+
+\code{ofile} is the name of the polygon file.
+
 \code{tag} is the tag given to all the polygon's segments.
+If \code{tag} is not supplied then the segment will not effect the water
+flow, it will only effect the mesh generation.
+
 \code{region_tag} is the tag given to all the polygon's segments.  If
-it is a string the one value will be assigned to all regions.  If it
+it is a string the tag will be assigned to all regions.  If it
 is a list the first value in the list will be applied to the first
-polygon etc.  If \code{tag} is not given a value it defaults to None,
-which means the segement will not effect the water flow, it will only
-effect the mesh generation.
+polygon etc.
 
 This function can be used to import building footprints.
@@ -1638,26 +1789,24 @@
 \label{sec:initialising the domain}
 
-%Include description of the class Domain and the module domain.
-
-%FIXME (Ole): This is also defined in a later chapter
-%\declaremodule{standard}{...domain}
-
-\begin{classdesc} {Domain} {source=None,
-                 triangles=None,
-                 boundary=None,
-                 conserved_quantities=None,
-                 other_quantities=None,
-                 tagged_elements=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}
+\begin{classdesc}{Domain}{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}
 Module: \refmodule{abstract_2d_finite_volumes.domain}
 
-This class is used to create an instance of a data structure used to
+This class is used to create an instance of a structure used to
 store and manipulate data associated with a mesh. The mesh is
 specified either by assigning the name of a mesh file to
@@ -1668,188 +1817,207 @@
 \subsection{Key Methods of Domain}
 
-\begin{methoddesc} {set\_name}{name}
-    Module: \refmodule{abstract\_2d\_finite\_volumes.domain},
-    page \pageref{mod:domain}
-
-    Assigns the name \code{name} to the domain.
-\end{methoddesc}
-
-\begin{methoddesc} {get\_name}{}
-    Module: \module{abstract\_2d\_finite\_volumes.domain}
-
-    Returns the name assigned to the domain by \code{set\_name}. If no name has been
-    assigned, returns \code{'domain'}.
-\end{methoddesc}
-
-\begin{methoddesc} {set\_datadir}{name}
-    Module: \module{abstract\_2d\_finite\_volumes.domain}
-
-    Specifies the directory used for SWW files, assigning it to the
-    pathname \code{name}. The default value, before
-    \code{set\_datadir} has been run, is the value \code{default\_datadir}
-    specified in \code{config.py}.
-
-    Since different operating systems use different formats for specifying pathnames,
-    it is necessary to specify path separators using the Python code \code{os.sep}, rather than
-    the operating-specific ones such as '$\slash$' or '$\backslash$'.
-    For this to work you will need to include the statement \code{import os}
-    in your code, before the first appearance of \code{set\_datadir}.
-
-    For example, to set the data directory to a subdirectory
-    \code{data} of the directory \code{project}, you could use
-    the statements:
-
-    \begin{verbatim}
+\begin{methoddesc}{\emph{<domain>}.set_name}{name}
+Module: \refmodule{abstract_2d_finite_volumes.domain},
+page \pageref{mod:domain}
+
+\code{name} is used to name the domain.  The \code{name} is also used to identify the output SWW file.
+If no name is assigned to a domain, the assumed name is \code{'domain'}.
+\end{methoddesc}
+
+\begin{methoddesc}{\emph{<domain>}.get_name}{}
+Module: \module{abstract_2d_finite_volumes.domain}
+
+Returns the name assigned to the domain by \code{set_name()}. If no name has been
+assigned, returns \code{'domain'}.
+\end{methoddesc}
+
+\begin{methoddesc}{\emph{<domain>}.set_datadir}{path}
+Module: \module{abstract_2d_finite_volumes.domain}
+
+\code{path} specifies the path to the directory used to store SWW files.
+
+Before this method is used to set the SWW directory path, the assumed directory
+path is \code{default_datadir} specified in \code{config.py}.
+
+Since different operating systems use different formats for specifying pathnames
+it is necessary to specify path separators using the Python code \code{os.sep} rather than
+the operating-specific ones such as '$\slash$' or '$\backslash$'.
+For this to work you will need to include the statement \code{import os}
+in your code, before the first use of \code{set_datadir()}.
+
+For example, to set the data directory to a subdirectory
+\code{data} of the directory \code{project}, you could use
+the statements:
+
+\begin{verbatim}
 import os
 domain.set_datadir{'project' + os.sep + 'data'}
-    \end{verbatim}
-\end{methoddesc}
-
-\begin{methoddesc} {get\_datadir}{}
-    Module: \module{abstract\_2d\_finite\_volumes.domain}
-
-    Returns the data directory set by \code{set\_datadir} or,
-    if \code{set\_datadir} has not
-    been run, returns the value \code{default\_datadir} specified in
-    \code{config.py}.
-\end{methoddesc}
-
-\begin{methoddesc} {set\_minimum_allowed_height}{}
-    Module: \module{shallow\_water.shallow\_water\_domain}
-
-    Set the minimum depth (in meters) that will be recognised in
-    the numerical scheme (including limiters and flux computations)
-
-    Default value is $10^{-3}$ m, but by setting this to a greater value,
-    e.g.\ for large scale simulations, the computation time can be
-    significantly reduced.
-\end{methoddesc}
-
-\begin{methoddesc} {set\_minimum_storable_height}{}
-    Module: \module{shallow\_water.shallow\_water\_domain}
-
-    Sets the minimum depth that will be recognised when writing
-    to an sww file. This is useful for removing thin water layers
-    that seems to be caused by friction creep.
-\end{methoddesc}
-
-\begin{methoddesc} {set\_maximum_allowed_speed}{}
-    Module: \module{shallow\_water.shallow\_water\_domain}
-
-    Set the maximum particle speed that is allowed in water
-    shallower than minimum_allowed_height. This is useful for
-    controlling speeds in very thin layers of water and at the same time
-    allow some movement avoiding pooling of water.
-\end{methoddesc}
-
-\begin{methoddesc} {set\_time}{time=0.0}
-    Module: \module{abstract\_2d\_finite\_volumes.domain}
-
-    Sets the initial time, in seconds, for the simulation. The
-    default is 0.0.
-\end{methoddesc}
-
-\begin{methoddesc} {set\_default\_order}{n}
-    Sets the default (spatial) order to the value specified by
-    \code{n}, which must be either 1 or 2. (Assigning any other value
-    to \code{n} will cause an error.)
-\end{methoddesc}
-
-\begin{methoddesc} {set\_store\_vertices\_uniquely}{flag}
+\end{verbatim}
+\end{methoddesc}
+
+\begin{methoddesc}{\emph{<domain>}.get_datadir}{}
+Module: \module{abstract_2d_finite_volumes.domain}
+
+Returns the path to the directory where SWW files will be stored.
+
+If the path has not previously been set with \code{set_datadir()} this method
+will return the value \code{default_datadir} specified in \code{config.py}.
+\end{methoddesc}
+
+\begin{methoddesc}{\emph{<domain>}.set_minimum_allowed_height}{minimum_allowed_height}
+Module: \module{shallow_water.shallow_water_domain}
+
+Set the minimum depth (in metres) that will be recognised in
+the numerical scheme (including limiters and flux computations)
+
+\code{minimum_allowed_height} is the new minimum allowed height value.
+
+Default value is $10^{-3}$ metre, but by setting this to a greater value,
+e.g. for large scale simulations, the computation time can be
+significantly reduced.
+\end{methoddesc}
+
+\begin{methoddesc}{\emph{<domain>}.set_minimum_storable_height}{minimum_storable_height}
+Module: \module{shallow_water.shallow_water_domain}
+
+Sets the minimum depth that will be recognised when writing
+to an SWW file. This is useful for removing thin water layers
+that seems to be caused by friction creep.
+
+\code{minimum_storable_height} is the new minimum storable height value.
+\end{methoddesc}
+
+\begin{methoddesc}{\emph{<domain>}.set_maximum_allowed_speed}{maximum_allowed_speed}
+Module: \module{shallow_water.shallow_water_domain}
+
+Set the maximum particle speed that is allowed in water
+shallower than \code{minimum_allowed_height}. This is useful for
+controlling speeds in very thin layers of water and at the same time
+allow some movement avoiding pooling of water.
+
+\code{maximum_allowed_speed} sets the maximum allowed speed value.
+\end{methoddesc}
+
+\begin{methoddesc}{\emph{<domain>}.set_time}{time=0.0}
+Module: \module{abstract_2d_finite_volumes.domain}
+
+\code{time} sets the initial time, in seconds, for the simulation. The
+default is 0.0.
+\end{methoddesc}
+
+\begin{methoddesc}{\emph{<domain>}.set_default_order}{n}
+Module: \module{abstract_2d_finite_volumes.domain}
+
+Sets the default (spatial) order to the value specified by
+\code{n}, which must be either 1 or 2. (Assigning any other value
+to \code{n} will cause an error.)
+\end{methoddesc}
+
+\begin{methoddesc}{\emph{<domain>}.set_store_vertices_uniquely}{flag, reduction=None}
+Module: \module{shallow_water.shallow_water_domain}
+
 Decide whether vertex values should be stored uniquely as
 computed in the model or whether they should be reduced to one
 value per vertex using averaging.
 
-Triangles stored in the sww file can be discontinuous reflecting
+\code{flag} may be \code{True} (meaning do not smooth vertex values) or \code{False}.
+
+\code{reduction} defines the smoothing operation if \code{flag} is \code{True}.  If not
+supplied, \code{reduction} is assumed to be \code{mean}.
+
+Triangles stored in the SWW file can be discontinuous reflecting
 the internal representation of the finite-volume scheme
-(this is a feature allowing for arbitrary steepness of the water surface gradient as well as the momentum gradients).
-However, for visual purposes and also for use with \code{Field\_boundary}
-(and \code{File\_boundary}) it is often desirable to store triangles
+(this is a feature allowing for arbitrary steepness of the water surface gradient
+as well as the momentum gradients).
+However, for visual purposes and also for use with \code{Field_boundary}
+(and \code{File_boundary}), it is often desirable to store triangles
 with values at each vertex point as the average of the potentially
 discontinuous numbers found at vertices of different triangles sharing the
 same vertex location.
 
-Storing one way or the other is controlled in ANUGA through the method
-\code{domain.store\_vertices\_uniquely}. Options are
+Storing one way or the other is controlled in \anuga through the method
+\code{<domain>.store_vertices_uniquely()}. Options are:
 \begin{itemize}
-  \item \code{domain.store\_vertices\_uniquely(True)}: Allow discontinuities in the sww file
-  \item \code{domain.store\_vertices\_uniquely(False)}: (Default).
+  \item \code{<domain>.store_vertices_uniquely(True)}: Allow discontinuities in the SWW file
+  \item \code{<domain>.store_vertices_uniquely(False)}: (Default).
   Average values
-  to ensure continuity in sww file. The latter also makes for smaller
-  sww files.
+  to ensure continuity in SWW file. The latter also makes for smaller
+  SWW files.
 \end{itemize}
 
-Note that when model data in the sww file are averaged (i.e. not stored uniquely), 
-then there will most likely be a small discrepancy between values extracted from the sww 
+Note that when model data in the SWW file are averaged (i.e. not stored uniquely), 
+then there will most likely be a small discrepancy between values extracted from the SWW 
 file and the same data stored in the model domain. This must be borne in mind when comparing
-data from the sww files with that of the model internally. 
+data from the SWW files with that of the model internally. 
 \end{methoddesc}
 
 % Structural methods
-\begin{methoddesc}{get\_nodes}{absolute=False}
-    Return x,y coordinates of all nodes in mesh.
-
-    The nodes are ordered in an Nx2 array where N is the number of nodes.
-    This is the same format they were provided in the constructor
-    i.e. without any duplication.
-
-    Boolean keyword argument absolute determines whether coordinates
-    are to be made absolute by taking georeference into account
-    Default is False as many parts of ANUGA expects relative coordinates.
-\end{methoddesc}
-
-\begin{methoddesc}{get\_vertex\_coordinates}{absolute=False}
-    \label{pg:get vertex coordinates}
-    Return vertex coordinates for all triangles.
-
-    Return all vertex coordinates for all triangles as a 3*M x 2 array
-    where the jth vertex of the ith triangle is located in row 3*i+j and
-    M the number of triangles in the mesh.
-
-    Boolean keyword argument absolute determines whether coordinates
-    are to be made absolute by taking georeference into account
-    Default is False as many parts of ANUGA expects relative coordinates.
-\end{methoddesc}
-
-\begin{methoddesc}{get\_centroid\_coordinates}{absolute=False}
-
-    Return centroid coordinates for all triangles.
-
-    Return all centroid coordinates for all triangles as a M x 2 array
+\begin{methoddesc}{\emph{<domain>}.get_nodes}{absolute=False}
+Module: \module{abstract_2d_finite_volumes.domain}
+
+Return x,y coordinates of all nodes in the domain mesh.  The nodes are ordered
+in an \code{Nx2} array where N is the number of nodes.  This is the same format
+they were provided in the constructor i.e. without any duplication.
+
+\code{absolute} is a boolean which determines whether coordinates
+are to be made absolute by taking georeference into account.
+Default is \code{False} as many parts of \anuga expect relative coordinates.
+\end{methoddesc}
+
+\begin{methoddesc}{\emph{<domain>}.get_vertex_coordinates}{absolute=False}
+Module: \module{abstract_2d_finite_volumes.domain}
+
+\label{pg:get vertex coordinates}
+Return vertex coordinates for all triangles as a \code{3*Mx2} array
+where the jth vertex of the ith triangle is located in row 3*i+j and
+M is the number of triangles in the mesh.
+
+\code{absolute} is a boolean which determines whether coordinates
+are to be made absolute by taking georeference into account.
+Default is \code{False} as many parts of \anuga expect relative coordinates.
+\end{methoddesc}
+
+\begin{methoddesc}{\emph{<domain>}.get_centroid_coordinates}{absolute=False}
+Module: \module{abstract_2d_finite_volumes.domain}
+
+Return centroid coordinates for all triangles as an \code{Mx2} array.
     
-    Boolean keyword argument absolute determines whether coordinates
-    are to be made absolute by taking georeference into account
-    Default is False as many parts of ANUGA expects relative coordinates.
-\end{methoddesc}
-
-\begin{methoddesc}{get\_triangles}{indices=None}
-
-    Return Mx3 integer array where M is the number of triangles.
-    Each row corresponds to one triangle and the three entries are
-    indices into the mesh nodes which can be obtained using the method
-    get\_nodes()
-
-    Optional argument, indices is the set of triangle ids of interest.
-\end{methoddesc}
-
-\begin{methoddesc}{get\_disconnected\_triangles}{}
-    Get mesh based on nodes obtained from get_vertex_coordinates.
-
-    Return array Mx3 array of integers where each row corresponds to
-    a triangle. A triangle is a triplet of indices into
-    point coordinates obtained from get_vertex_coordinates and each
-    index appears only once.\\
-
-    This provides a mesh where no triangles share nodes
-    (hence the name disconnected triangles) and different
-    nodes may have the same coordinates.\\
-
-    This version of the mesh is useful for storing meshes with
-    discontinuities at each node and is e.g. used for storing
-    data in sww files.\\
-
-    The triangles created will have the format:
-
-    \begin{verbatim}
+\code{absolute} is a boolean which determines whether coordinates
+are to be made absolute by taking georeference into account.
+Default is \code{False} as many parts of \anuga expect relative coordinates.
+\end{methoddesc}
+
+\begin{methoddesc}{\emph{<domain>}.get_triangles}{indices=None}
+Module: \module{abstract_2d_finite_volumes.domain}
+
+Return an \code{Mx3} integer array where M is the number of triangles.
+Each row corresponds to one triangle and the three entries are
+indices into the mesh nodes which can be obtained using the method
+\code{get_nodes()}.
+
+\code{indices}, if specified, is the set of triangle \code{id}s of interest.
+\end{methoddesc}
+
+\begin{methoddesc}{\emph{<domain>}.get_disconnected_triangles}{}
+Module: \module{abstract_2d_finite_volumes.domain}
+
+Get the domain mesh based on nodes obtained from \code{get_vertex_coordinates()}.
+
+Returns an \code{Mx3} array of integers where each row corresponds to
+a triangle. A triangle is a triplet of indices into
+point coordinates obtained from \code{get_vertex_coordinates()} and each
+index appears only once.
+
+This provides a mesh where no triangles share nodes
+(hence the name disconnected triangles) and different
+nodes may have the same coordinates.
+
+This version of the mesh is useful for storing meshes with
+discontinuities at each node and is e.g. used for storing
+data in SWW files.
+
+The triangles created will have the format:
+
+\begin{verbatim}
 [[0,1,2],
  [3,4,5],
@@ -1857,5 +2025,5 @@
  ...
  [3*M-3 3*M-2 3*M-1]]
-    \end{verbatim}
+\end{verbatim}
 \end{methoddesc}
 
@@ -1871,93 +2039,94 @@
 of initial conditions to refer to initial values for variables which need
 prescription to solve the shallow water wave equation. Further, it must be noted
-that \code{set_quantity} does not necessarily have to be used in the initial
+that \code{set_quantity()} does not necessarily have to be used in the initial
 condition setting; it can be used at any time throughout the simulation.
 
-\begin{methoddesc}{set\_quantity}{name,
-    numeric=None,
-    quantity=None,
-    function=None,
-    geospatial_data=None,
-    expression=None,    
-    filename=None,
-    attribute_name=None,
-    alpha=None,
-    location='vertices',
-    indices=None,
-    verbose=False,
-    use_cache=False}
-  Module: \module{abstract\_2d\_finite\_volumes.domain}
-  (see also \module{abstract\_2d\_finite\_volumes.quantity.set\_values})
-
-  This function is used to assign values to individual quantities for a
-  domain. It is very flexible and can be used with many data types: a
-  statement of the form \code{domain.set\_quantity(name, x)} can be used
-  to define a quantity having the name \code{name}, where the other
-  argument \code{x} can be any of the following:
-
-  \begin{itemize}
-    \item a number, in which case all vertices in the mesh gets that for
-          the quantity in question.
-    \item a list of numbers or a numeric array ordered the same way as the mesh vertices.
-    \item a function (e.g.\ see the samples introduced in Chapter 2)
-    \item an expression composed of other quantities and numbers, arrays, lists (for
-          example, a linear combination of quantities, such as
-    \code{domain.set\_quantity('stage','elevation'+x))}
-    \item the name of a file from which the data can be read. In this case, the optional
-          argument attribute\_name will select which attribute to use from the file. If left out,
-          set\_quantity will pick one. This is useful in cases where there is only one attribute.
-    \item a geospatial dataset (See Section \ref{sec:geospatial}).
-          Optional argument attribute\_name applies here as with files.
-  \end{itemize}
-
-  Exactly one of the arguments
-    numeric, quantity, function, points, filename
-  must be present.
-
-  \code{set_quantity} will look at the type of the second argument (\code{numeric}) and
-  determine what action to take.
-
-  Values can also be set using the appropriate keyword arguments.
-  If x is a function, for example, \code{domain.set\_quantity(name, x)}, \code{domain.set\_quantity(name, numeric=x)},
-  and \code{domain.set\_quantity(name, function=x)} are all equivalent.
-
-  Other optional arguments are
-  \begin{itemize}
-    \item \code{indices} which is a list of ids of triangles to which set\_quantity should apply its assignment of values.
-    \item \code{location} determines which part of the triangles to assign
-                          to. Options are 'vertices' (default), 'edges', 'unique vertices', and 'centroids'.
-                          If 'vertices' are use, edge and centroid values are automatically computed as the
-                          appropriate averages. This option ensures continuity of the surface.
-                          If, on the other hand, 'centroids' is used vertex and edge values will be set to the
-                          same value effectively creating a piecewise constant surface with possible discontinuities at the edges.
-  \end{itemize}
-
-  \anuga provides a number of predefined initial conditions to be used
-  with \code{set\_quantity}. See for example callable object \code{slump\_tsunami} below.
-\end{methoddesc}
-
-\begin{methoddesc}{add\_quantity}{name,
-    numeric=None,
-    quantity=None,
-    function=None,
-    geospatial_data=None,
-    expression=None,
-    filename=None,
-    attribute_name=None,
-    alpha=None,
-    location='vertices',
-    indices=None,
-    verbose=False,
-    use_cache=False}
-  Module: \module{abstract\_2d\_finite\_volumes.domain}
-  (see also \module{abstract\_2d\_finite\_volumes.domain.set\_quantity})
-
-  \label{add quantity}
-  This function is used to \emph{add} values to individual quantities for a
-  domain. It has the same syntax as \code{domain.set\_quantity(name, x)}.
-
-  A typical use of this function is to add structures to an existing elevation model:
-
-  \begin{verbatim} 
+\begin{methoddesc}{\emph{<domain>}.set_quantity}{numeric=None,
+                                                 quantity=None,
+                                                 function=None,
+                                                 geospatial_data=None,
+                                                 filename=None,
+                                                 attribute_name=None,
+                                                 alpha=None,
+                                                 location='vertices',
+                                                 polygon=None,
+                                                 indices=None,
+                                                 smooth=False,
+                                                 verbose=False,
+                                                 use_cache=False}
+Module: \module{abstract_2d_finite_volumes.domain} \\
+(This method passes off to \module{abstract_2d_finite_volumes.quantity.set_values()})
+
+This function is used to assign values to individual quantities for a
+domain. It is very flexible and can be used with many data types: a
+statement of the form \code{\emph{<domain>}.set_quantity(name, x)} can be used
+to define a quantity having the name \code{name}, where the other
+argument \code{x} can be any of the following:
+
+\begin{itemize}
+  \item a number, in which case all vertices in the mesh gets that for
+        the quantity in question
+  \item a list of numbers or a numeric array ordered the same way as the mesh vertices
+  \item a function (e.g. see the samples introduced in Chapter 2)
+  \item an expression composed of other quantities and numbers, arrays, lists (for
+        example, a linear combination of quantities, such as
+        \code{\emph{<domain>}.set_quantity('stage','elevation'+x))}
+  \item the name of a file from which the data can be read. In this case, the optional
+        argument \code{attribute_name} will select which attribute to use from the file. If left out,
+        \code{set_quantity()} will pick one. This is useful in cases where there is only one attribute
+  \item a geospatial dataset (See Section \ref{sec:geospatial}).
+        Optional argument \code{attribute_name} applies here as with files
+\end{itemize}
+
+Exactly one of the arguments \code{numeric}, \code{quantity}, \code{function},
+\code{geospatial_data} and \code{filename} must be present.
+
+\code{set_quantity()} will look at the type of the \code{numeric} and
+determine what action to take.
+
+Values can also be set using the appropriate keyword arguments.
+If \code{x} is a function, for example, \code{domain.set_quantity(name, x)}, \code{domain.set_quantity(name, numeric=x)},
+and \code{domain.set_quantity(name, function=x)} are all equivalent.
+
+Other optional arguments are:
+\begin{itemize}
+  \item \code{indices} which is a list of ids of triangles to which \code{set_quantity()}
+        should apply its assignment of values.
+  \item \code{location} determines which part of the triangles to assign to.
+        Options are 'vertices' (the default), 'edges', 'unique vertices', and 'centroids'.
+        If 'vertices' is used, edge and centroid values are automatically computed as the
+        appropriate averages. This option ensures continuity of the surface.
+        If, on the other hand, 'centroids' is used, vertex and edge values will be set to the
+        same value effectively creating a piecewise constant surface with possible
+        discontinuities at the edges.
+\end{itemize}
+
+\anuga provides a number of predefined initial conditions to be used
+with \code{set_quantity()}. See for example callable object \code{slump_tsunami} below.
+\end{methoddesc}
+
+\begin{methoddesc}{\emph{<domain>}.add_quantity}{numeric=None,
+                                                 quantity=None,
+                                                 function=None,
+                                                 geospatial_data=None,
+                                                 filename=None,
+                                                 attribute_name=None,
+                                                 alpha=None,
+                                                 location='vertices',
+                                                 polygon=None,
+                                                 indices=None,
+                                                 smooth=False,
+                                                 verbose=False,
+                                                 use_cache=False}
+Module: \module{abstract_2d_finite_volumes.domain} \\
+(passes off to \module{abstract_2d_finite_volumes.domain.set_quantity()})
+
+\label{add quantity}
+This function is used to \emph{add} values to individual quantities for a
+domain. It has the same syntax as \code{\emph{<domain>}.set_quantity(name, x)}.
+
+A typical use of this function is to add structures to an existing elevation model:
+
+\begin{verbatim} 
 # Create digital elevation model from points file
 domain.set_quantity('elevation', filename='elevation_file.pts, verbose=True)
@@ -1973,148 +2142,143 @@
 
 domain.add_quantity('elevation', Polygon_function(B, default=0.0))
-  \end{verbatim}
-\end{methoddesc}
-
-\begin{funcdesc}{set_region}{tag, quantity, X, location='vertices'}
-  Module: \module{abstract\_2d\_finite\_volumes.domain}
-
-  (see also \module{abstract\_2d\_finite\_volumes.quantity.set\_values})
-
-  This function is used to assign values to individual quantities given
-  a regional tag.   It is similar to \code{set\_quantity}.
-  For example, if in the mesh-generator a regional tag of 'ditch' was
-  used, set\_region can be used to set elevation of this region to
-  -10m. X is the constant or function to be applied to the quantity,
-  over the tagged region.  Location describes how the values will be
-  applied.  Options are 'vertices' (default), 'edges', 'unique
-  vertices', and 'centroids'.
-
-  This method can also be called with a list of region objects.  This is
-  useful for adding quantities in regions, and having one quantity
-  value based on another quantity. See  \module{abstract\_2d\_finite\_volumes.region} for
-  more details.
+\end{verbatim}
+\end{methoddesc}
+
+\begin{methoddesc}{\emph{<domain>}.set_region}{tag, quantity, X, location='vertices'}
+Module: \module{abstract_2d_finite_volumes.domain} \\
+(see also \module{abstract_2d_finite_volumes.quantity.set_values})
+
+This function is used to assign values to individual quantities given
+a regional tag.   It is similar to \code{set_quantity()}.
+
+For example, if in the mesh-generator a regional tag of 'ditch' was
+used, \code{set_region()} can be used to set elevation of this region to
+-10m. \code{X} is the constant or function to be applied to the \code{quantity},
+over the tagged region.  \code{location} describes how the values will be
+applied.  Options are 'vertices' (the default), 'edges', 'unique
+vertices', and 'centroids'.
+
+This method can also be called with a list of region objects.  This is
+useful for adding quantities in regions, and having one quantity
+value based on another quantity. See  \module{abstract_2d_finite_volumes.region} for
+more details.
+\end{methoddesc}
+
+\begin{funcdesc}{slump_tsunami}{length, depth, slope, width=None, thickness=None,
+                                radius=None, dphi=0.48, x0=0.0, y0=0.0, alpha=0.0,
+                                gravity=9.8, gamma=1.85,
+                                massco=1, dragco=1, frictionco=0,
+                                dx=None, kappa=3.0, kappad=1.0, zsmall=0.01, scale=None,
+                                domain=None,
+                                verbose=False}
+Module: \module{shallow\_water.smf}
+
+This function returns a callable object representing an initial water
+displacement generated by a submarine sediment failure. These failures can take the form of
+a submarine slump or slide. In the case of a slide, use \code{slide_tsunami} instead.
+
+\code{length} is the length of the slide or slump.
+
+\code{depth} is the water depth to the centre of the sediment mass.
+
+\code{slope} is the bathymetric slope.
+
+Other slump or slide parameters can be included if they are known.
 \end{funcdesc}
 
-\begin{funcdesc}{slump_tsunami}{length, depth, slope, width=None, thickness=None,
-                x0=0.0, y0=0.0, alpha=0.0,
-                gravity=9.8, gamma=1.85,
-                massco=1, dragco=1, frictionco=0, psi=0,
-                dx=None, kappa=3.0, kappad=0.8, zsmall=0.01,
-                domain=None,
-                verbose=False}
-  Module: \module{shallow\_water.smf}
-
-  This function returns a callable object representing an initial water
-  displacement generated by a submarine sediment failure. These failures can take the form of
-  a submarine slump or slide. In the case of a slide, use \code{slide_tsunami} instead.
-
-  The arguments include as a minimum, the slump or slide length, the water depth to the centre of sediment
-  mass, and the bathymetric slope. Other slump or slide parameters can be included if they are known.
+\begin{funcdesc}{\emph{<callable_object>} = file_function}{filename,
+                                domain=None,
+                                quantities=None,
+                                interpolation_points=None,
+                                time_thinning=1,
+                                time_limit=None,
+                                verbose=False,
+                                use_cache=False,
+                                boundary_polygon=None}
+Module: \module{abstract_2d_finite_volumes.util}
+
+Reads the time history of spatial data for specified interpolation points from
+a NetCDF file and returns a callable object.  Values returned from the \code{\emph{<callable_object>}}
+are interpolated values based on the input file using the underlying \code{interpolation_function}.
+
+\code{filename} is the name of the input file.  Could be an SWW or STS file.
+
+\code{quantities} is either the name of a single quantity to be
+interpolated or a list of such quantity names. In the second case, the resulting
+function will return a tuple of values -- one for each quantity.
+
+\code{interpolation_points} is a list of absolute coordinates or a
+geospatial object for points at which values are sought.
+
+\code{boundary_polygon} is a list of coordinates specifying the vertices of the boundary.
+This must be the same polygon as used when calling \code{create_mesh_from_regions()}.
+This argument can only be used when reading boundary data from an STS format file.
+
+The model time stored within the file function can be accessed using
+the method \code{\emph{<callable_object>}.get_time()}
+
+The underlying algorithm used is as follows:\\
+Given a time series (i.e. a series of values associated with
+different times), whose values are either just numbers, a set of
+numbers defined at the vertices of a triangular mesh (such as those
+stored in SWW files) or a set of
+numbers defined at a number of points on the boundary (such as those
+stored in STS files), \code{Interpolation_function()} is used to
+create a callable object that interpolates a value for an arbitrary
+time \code{t} within the model limits and possibly a point \code{(x, y)}
+within a mesh region.
+
+The actual time series at which data is available is specified by
+means of an array \code{time} of monotonically increasing times. The
+quantities containing the values to be interpolated are specified in
+an array -- or dictionary of arrays (used in conjunction with the
+optional argument \code{quantity_names}) -- called
+\code{quantities}. The optional arguments \code{vertex_coordinates}
+and \code{triangles} represent the spatial mesh associated with the
+quantity arrays. If omitted the function must be created using an STS file
+or a TMS file.
+
+Since, in practice, values need to be computed at specified points,
+the syntax allows the user to specify, once and for all, a list
+\code{interpolation_points} of points at which values are required.
+In this case, the function may be called using the form \code{\emph{<callable_object>}(t, id)},
+where \code{id} is an index for the list \code{interpolation_points}.
 \end{funcdesc}
 
-\begin{funcdesc}{file\_function}{filename,
-    domain=None,
-    quantities=None,
-    interpolation_points=None,
-    verbose=False,
-    use_cache=False}
-  Module: \module{abstract\_2d\_finite\_volumes.util}
-
-  Reads the time history of spatial data for
-  specified interpolation points from a NetCDF file (\code{filename})
-  and returns
-  a callable object. \code{filename} could be a \code{sww} or \code{sts} file.
-  Returns interpolated values based on the input
-  file using the underlying \code{interpolation\_function}.
-
-  \code{quantities} is either the name of a single quantity to be
-  interpolated or a list of such quantity names. In the second case, the resulting
-  function will return a tuple of values -- one for each quantity.
-
-  \code{interpolation\_points} is a list of absolute coordinates or a
-  geospatial object
-  for points at which values are sought.
-
-  \code{boundary_polygon} is a list of coordinates specifying the vertices of the boundary.
-  This must be the same polygon as used when calling \code{create_mesh_from_regions}.
-  This argument can only be used when reading boundary data from the STS format.
-
-  The model time stored within the file function can be accessed using
-  the method \code{f.get\_time()}
-
-  The underlying algorithm used is as follows:\\
-  Given a time series (i.e.\ a series of values associated with
-  different times), whose values are either just numbers, a set of
-  numbers defined at the vertices of a triangular mesh (such as those
-  stored in SWW files) or a set of
-  numbers defined at a number of points on the boundary (such as those
-  stored in STS files), \code{Interpolation\_function} is used to
-  create a callable object that interpolates a value for an arbitrary
-  time \code{t} within the model limits and possibly a point \code{(x,
-  y)} within a mesh region.
-
-  The actual time series at which data is available is specified by
-  means of an array \code{time} of monotonically increasing times. The
-  quantities containing the values to be interpolated are specified in
-  an array -- or dictionary of arrays (used in conjunction with the
-  optional argument \code{quantity\_names}) -- called
-  \code{quantities}. The optional arguments \code{vertex\_coordinates}
-  and \code{triangles} represent the spatial mesh associated with the
-  quantity arrays. If omitted the function must be created using an STS file
-  or a TMS file.
-
-  Since, in practice, values need to be computed at specified points,
-  the syntax allows the user to specify, once and for all, a list
-  \code{interpolation\_points} of points at which values are required.
-  In this case, the function may be called using the form \code{f(t,
-  id)}, where \code{id} is an index for the list
-  \code{interpolation\_points}.
-\end{funcdesc}
-
-% FIXME (OLE): Why has this been commented out?
-%%%
-%% \begin{classdesc}{Interpolation\_function}{self,
-%%     time,
-%%     quantities,
-%%     quantity_names = None,
-%%     vertex_coordinates = None,
-%%     triangles = None,
-%%     interpolation_points = None,
-%%     verbose = False}
-%% Module: \module{abstract\_2d\_finite\_volumes.least\_squares}
-
-%% Given a time series (i.e.\ a series of values associated with
-%% different times), whose values are either just numbers or a set of
-%% numbers defined at the vertices of a triangular mesh (such as those
-%% stored in SWW files), \code{Interpolation\_function} is used to
-%% create a callable object that interpolates a value for an arbitrary
-%% time \code{t} within the model limits and possibly a point \code{(x,
-%% y)} within a mesh region.
-
-%% The actual time series at which data is available is specified by
-%% means of an array \code{time} of monotonically increasing times. The
-%% quantities containing the values to be interpolated are specified in
-%% an array -- or dictionary of arrays (used in conjunction with the
-%% optional argument \code{quantity\_names}) -- called
-%% \code{quantities}. The optional arguments \code{vertex\_coordinates}
-%% and \code{triangles} represent the spatial mesh associated with the
-%% quantity arrays. If omitted the function created by
-%% \code{Interpolation\_function} will be a function of \code{t} only.
-
-%% Since, in practice, values need to be computed at specified points,
-%% the syntax allows the user to specify, once and for all, a list
-%% \code{interpolation\_points} of points at which values are required.
-%% In this case, the function may be called using the form \code{f(t,
-%% id)}, where \code{id} is an index for the list
-%% \code{interpolation\_points}.
-
-%% \end{classdesc}
-
-%%%
-%\begin{funcdesc}{set\_region}{functions}
-%[Low priority. Will be merged into set\_quantity]
-
-%Module:\module{abstract\_2d\_finite\_volumes.domain}
-%\end{funcdesc}
+\begin{classdesc}{\emph{<callable_object>} = Interpolation_function}{time,
+                                          quantities,
+                                          quantity_names=None,
+                                          vertex_coordinates=None,
+                                          triangles=None,
+                                          interpolation_points=None,
+                                          time_thinning=1,
+                                          verbose=False,
+                                          gauge_neighbour_id=None}
+Module: \module{fit_interpolate.interpolate}
+
+Given a time series (i.e. a series of values associated with
+different times) whose values are either just numbers or a set of
+numbers defined at the vertices of a triangular mesh (such as those
+stored in SWW files), \code{Interpolation_function} is used to
+create a callable object that interpolates a value for an arbitrary
+time \code{t} within the model limits and possibly a point \code{(x, y)}
+within a mesh region.
+
+The actual time series at which data is available is specified by
+means of an array \code{time} of monotonically increasing times. The
+quantities containing the values to be interpolated are specified in
+an array -- or dictionary of arrays (used in conjunction with the
+optional argument \code{quantity\_names}) -- called
+\code{quantities}. The optional arguments \code{vertex_coordinates}
+and \code{triangles} represent the spatial mesh associated with the
+quantity arrays. If omitted the function created by
+\code{Interpolation_function} will be a function of \code{t} only.
+
+Since, in practice, values need to be computed at specified points,
+the syntax allows the user to specify, once and for all, a list
+\code{interpolation_points} of points at which values are required.
+In this case, the function may be called using the form \code{f(t, id)},
+where \code{id} is an index for the list \code{interpolation_points}.
+\end{classdesc}
 
 
@@ -2123,155 +2287,173 @@
 
 \anuga provides a large number of predefined boundary conditions,
-represented by objects such as \code{Reflective\_boundary(domain)} and
-\code{Dirichlet\_boundary([0.2, 0.0, 0.0])}, described in the examples
+represented by objects such as \code{Reflective_boundary(domain)} and
+\code{Dirichlet_boundary([0.2, 0.0, 0.0])}, described in the examples
 in Chapter 2. Alternatively, you may prefer to ''roll your own'',
 following the method explained in Section \ref{sec:roll your own}.
 
-These boundary objects may be used with the function \code{set\_boundary} described below
+These boundary objects may be used with the function \code{set_boundary} described below
 to assign boundary conditions according to the tags used to label boundary segments.
 
-\begin{methoddesc}{set\_boundary}{boundary_map}
-  Module: \module{abstract\_2d\_finite\_volumes.domain}
-
-  This function allows you to assign a boundary object (corresponding to a
-  pre-defined or user-specified boundary condition) to every boundary segment that
-  has been assigned a particular tag.
-
-  This is done by specifying a dictionary \code{boundary\_map}, whose values are the boundary objects
-  and whose keys are the symbolic tags.
-\end{methoddesc}
-
-\begin{methoddesc} {get\_boundary\_tags}{}
-  Module: \module{abstract\_2d\_finite\_volumes.domain}
-
-  Returns a list of the available boundary tags.
+\begin{methoddesc}{\emph{<domain>}.set_boundary}{boundary_map}
+Module: \module{abstract_2d_finite_volumes.domain}
+
+This function allows you to assign a boundary object (corresponding to a
+pre-defined or user-specified boundary condition) to every boundary segment that
+has been assigned a particular tag.
+
+\code{boundary_map} is a dictionary of boundary objects keyed by symbolic tags.
+\end{methoddesc}
+
+\begin{methoddesc} {\emph{<domain>}.get_boundary_tags}{}
+Module: \module{abstract\_2d\_finite\_volumes.domain}
+
+Returns a list of the available boundary tags.
 \end{methoddesc}
 
 \subsection{Predefined boundary conditions}
 
-\begin{classdesc}{Reflective\_boundary}{Boundary}
-  Module: \module{shallow\_water}
-
-  Reflective boundary returns same conserved quantities as those present in
-  the neighbouring volume but reflected.
-
-  This class is specific to the shallow water equation as it works with the
-  momentum quantities assumed to be the second and third conserved quantities.
+\begin{classdesc}{Reflective_boundary}{domain=None}
+Module: \module{shallow_water}
+
+Reflective boundary returns same conserved quantities as those present in
+the neighbouring volume but reflected.
+
+This class is specific to the shallow water equation as it works with the
+momentum quantities assumed to be the second and third conserved quantities.
 \end{classdesc}
 
-\begin{classdesc}{Transmissive\_boundary}{domain=None}
+\begin{classdesc}{Transmissive_boundary}{domain=None}
   \label{pg: transmissive boundary}
-  Module: \module{abstract\_2d\_finite\_volumes.generic\_boundary\_conditions}
-
-  A transmissive boundary returns the same conserved quantities as
-  those present in the neighbouring volume.
-
-  The underlying domain must be specified when the boundary is instantiated.
+Module: \module{abstract_2d_finite_volumes.generic_boundary_conditions}
+
+A transmissive boundary returns the same conserved quantities as
+those present in the neighbouring volume.
+
+The underlying domain must be specified when the boundary is instantiated.
 \end{classdesc}
 
-\begin{classdesc}{Dirichlet\_boundary}{conserved_quantities=None}
-  Module: \module{abstract\_2d\_finite\_volumes.generic\_boundary\_conditions}
-
-  A Dirichlet boundary returns constant values for each of conserved
-  quantities. In the example of \code{Dirichlet\_boundary([0.2, 0.0, 0.0])},
-  the \code{stage} value at the boundary is 0.2 and the \code{xmomentum} and
-  \code{ymomentum} at the boundary are set to 0.0. The list must contain
-  a value for each conserved quantity.
+\begin{classdesc}{Dirichlet_boundary}{conserved_quantities=None}
+Module: \module{abstract_2d_finite_volumes.generic_boundary_conditions}
+
+A Dirichlet boundary returns constant values for each of conserved
+quantities. In the example of \code{Dirichlet_boundary([0.2, 0.0, 0.0])},
+the \code{stage} value at the boundary is 0.2 and the \code{xmomentum} and
+\code{ymomentum} at the boundary are set to 0.0. The list must contain
+a value for each conserved quantity.
 \end{classdesc}
 
-\begin{classdesc}{Time\_boundary}{domain = None, f = None}
-  Module: \module{abstract\_2d\_finite\_volumes.generic\_boundary\_conditions}
-
-  A time-dependent boundary returns values for the conserved
-  quantities as a function \code{f(t)} of time. The user must specify
-  the domain to get access to the model time.
-
-  Optional argument \code{default\_boundary} can be used to specify another boundary object
-  to be used in case model time exceeds the time availabel in the file used by \code{File\_boundary}.
-  The \code{default\_boundary} could be a simple Dirichlet condition or
-  even another \code{Time\_boundary} typically using data pertaining to another time interval.
+\begin{classdesc}{Time_boundary}{domain=None,
+                                 function=None,
+                                 default_boundary=None,
+                                 verbose=False}
+Module: \module{abstract_2d_finite_volumes.generic_boundary_conditions}
+
+A time-dependent boundary returns values for the conserved
+quantities as a function of time \code{function(t)}. The user must specify
+the domain to get access to the model time.
+
+Optional argument \code{default_boundary} can be used to specify another boundary object
+to be used in case model time exceeds the time available in the file used by \code{File_boundary}.
+The \code{default_boundary} could be a simple Dirichlet condition or
+even another \code{Time_boundary} typically using data pertaining to another time interval.
 \end{classdesc}
 
-\begin{classdesc}{File\_boundary}{Boundary}
-  Module: \module{abstract\_2d\_finite\_volumes.generic\_boundary\_conditions}
-
-  This method may be used if the user wishes to apply a SWW file, STS file or
-  a time series file (TMS) to a boundary segment or segments.
-  The boundary values are obtained from a file and interpolated to the
-  appropriate segments for each conserved quantity.
-
-  Optional argument \code{default\_boundary} can be used to specify another boundary object
-  to be used in case model time exceeds the time availabel in the file used by \code{File\_boundary}.
-  The \code{default\_boundary} could be a simple Dirichlet condition or
-  even another \code{File\_boundary} typically using data pertaining to another time interval.
+\begin{classdesc}{File_boundary}{filename,
+                                 domain,
+                                 time_thinning=1,
+                                 time_limit=None,
+                                 boundary_polygon=None,
+                                 default_boundary=None,
+                                 use_cache=False,
+                                 verbose=False}
+Module: \module{abstract_2d_finite_volumes.generic_boundary_conditions}
+
+This method may be used if the user wishes to apply a SWW file, STS file or
+a time series file (TMS) to a boundary segment or segments.
+The boundary values are obtained from a file and interpolated to the
+appropriate segments for each conserved quantity.
+
+Optional argument \code{default_boundary} can be used to specify another boundary object
+to be used in case model time exceeds the time available in the file used by \code{File_boundary}.
+The \code{default_boundary} could be a simple Dirichlet condition or
+even another \code{File_boundary} typically using data pertaining to another time interval.
 \end{classdesc}
 
-\begin{classdesc}{Field\_boundary}{Boundary}
-  Module: \module{shallow\_water.shallow\_water\_domain}
-
-  This method works in the same way as \code{File\_boundary} except that it
-  allows for the value of stage to be offset by a constant specified in the
-  keyword argument \code{mean\_stage}.
-
-  This functionality allows for models to be run for a range of tides using
-  the same boundary information (from .sts, .sww or .tms files). The tidal value
-  for each run would then be specified in the keyword argument
-  \code{mean\_stage}.
-  If \code{mean\_stage} = 0.0, \code{Field\_boundary} and \code{File\_boundary}
-  behave identically.
-
-  Note that if the optional argument \code{default\_boundary} is specified
-  it's stage value will be adjusted by \code{mean\_stage} just like the values
-  obtained from the file.
-
-  See \code{File\_boundary} for further details.
+\begin{classdesc}{Field_boundary}{filename,
+                                  domain,
+                                  mean_stage=0.0,
+                                  time_thinning=1,
+                                  time_limit=None,
+                                  boundary_polygon=None,
+                                  default_boundary=None,
+                                  use_cache=False,
+                                  verbose=False}
+Module: \module{shallow_water.shallow_water_domain}
+
+This method works in the same way as \code{File_boundary} except that it
+allows for the value of stage to be offset by a constant specified in the
+keyword argument \code{mean_stage}.
+
+This functionality allows for models to be run for a range of tides using
+the same boundary information (from STS, SWW or TMS files). The tidal value
+for each run would then be specified in the keyword argument \code{mean_stage}.
+If \code{mean_stage} = 0.0, \code{Field_boundary} and \code{File_boundary}
+behave identically.
+
+Note that if the optional argument \code{default_boundary} is specified
+its stage value will be adjusted by \code{mean_stage} just like the values
+obtained from the file.
+
+See \code{File_boundary} for further details.
 \end{classdesc}
 
-\begin{classdesc}{Transmissive\_momentum\_set\_stage\_boundary}{Boundary}
-  Module: \module{shallow\_water}
-  \label{pg: transmissive momentum set stage boundary}
-
-  This boundary returns same momentum conserved quantities as
-  those present in its neighbour volume but sets stage as in a Time\_boundary.
-  The underlying domain must be specified when boundary is instantiated
-
-  This type of boundary is useful when stage is known at the boundary as a
-  function of time, but momenta (or speeds) aren't.
-
-  This class is specific to the shallow water equation as it works with the
-  momentum quantities assumed to be the second and third conserved quantities.
-
-  In some circumstances, this boundary condition may cause numerical instabilities for similar
-  reasons as what has been observed with the simple fully transmissive boundary condition
-  (see Page \pageref{pg: transmissive boundary}). 
-  This could for example be the case if a planar wave is reflected out through this boundary.
+\begin{classdesc}{Transmissive_momentum_set_stage_boundary}{domain=None,
+                                                            function=None}
+Module: \module{shallow_water.shallow_water_domain}
+\label{pg: transmissive momentum set stage boundary}
+
+This boundary returns the same momentum conserved quantities as
+those present in its neighbour volume but sets stage as in a \code{Time_boundary}.
+The underlying domain must be specified when boundary is instantiated.
+
+This type of boundary is useful when stage is known at the boundary as a
+function of time, but momenta (or speeds) aren't.
+
+This class is specific to the shallow water equation as it works with the
+momentum quantities assumed to be the second and third conserved quantities.
+
+In some circumstances, this boundary condition may cause numerical instabilities for similar
+reasons as what has been observed with the simple fully transmissive boundary condition
+(see Page \pageref{pg: transmissive boundary}). 
+This could for example be the case if a planar wave is reflected out through this boundary.
 \end{classdesc}
 
-\begin{classdesc}{Transmissive\_stage\_zero\_momentum\_boundary}{Boundary}
-  Module: \module{shallow\_water}
-  \label{pg: transmissive stage zero momentum boundary}
-
-  This boundary returns same stage conserved quantities as
-  those present in its neighbour volume but sets momentum to zero.
-  The underlying domain must be specified when boundary is instantiated
-
-  This type of boundary is useful when stage is known at the boundary as a
-  function of time, but momentum should be set to zero. This is for example
-  the case where a boundary is needed in the ocean on the two sides perpendicular
-  to the coast to maintain the wave height of the incoming wave.
-
-  This class is specific to the shallow water equation as it works with the
-  momentum quantities assumed to be the second and third conserved quantities.
-
-  This boundary condition should not cause the numerical instabilities seen with the transmissive momentum 
-  boundary conditions (see Page \pageref{pg: transmissive boundary} and
-  Page \pageref{pg: transmissive momentum set stage boundary}).
+\begin{classdesc}{Transmissive_stage_zero_momentum_boundary}{domain=None}
+Module: \module{shallow_water}
+\label{pg: transmissive stage zero momentum boundary}
+
+This boundary returns same stage conserved quantities as
+those present in its neighbour volume but sets momentum to zero.
+The underlying domain must be specified when boundary is instantiated
+
+This type of boundary is useful when stage is known at the boundary as a
+function of time, but momentum should be set to zero. This is for example
+the case where a boundary is needed in the ocean on the two sides perpendicular
+to the coast to maintain the wave height of the incoming wave.
+
+This class is specific to the shallow water equation as it works with the
+momentum quantities assumed to be the second and third conserved quantities.
+
+This boundary condition should not cause the numerical instabilities seen with the transmissive momentum 
+boundary conditions (see Page \pageref{pg: transmissive boundary} and
+Page \pageref{pg: transmissive momentum set stage boundary}).
 \end{classdesc}
 
-\begin{classdesc}{Dirichlet\_discharge\_boundary}{Boundary}
-  Module: \module{shallow\_water}
-
-  Sets stage (stage0)
-  Sets momentum (wh0) in the inward normal direction.
+\begin{classdesc}{Dirichlet_discharge_boundary}{domain=None, stage0=None, wh0=None}
+Module: \module{shallow_water.shallow_water_domain}
+
+\code{stage0} sets the stage.
+
+\code{wh0} sets momentum in the inward normal direction.
 \end{classdesc}
 
@@ -2280,7 +2462,7 @@
 
 All boundary classes must inherit from the generic boundary class
-\code{Boundary} and have a method called \code{evaluate} which must
-take as inputs \code{self, vol\_id, edge\_id} where self refers to the
-object itself and vol\_id and edge\_id are integers referring to
+\code{Boundary} and have a method called \code{evaluate()} which must
+take as inputs \code{self}, \code{vol_id} and \code{edge_id} where \code{self} refers to the
+object itself and \code{vol_id} and \code{edge_id} are integers referring to
 particular edges. The method must return a list of three floating point
 numbers representing values for \code{stage},
@@ -2288,5 +2470,5 @@
 
 The constructor of a particular boundary class may be used to specify
-particular values or flags to be used by the \code{evaluate} method.
+particular values or flags to be used by the \code{evaluate()} method.
 Please refer to the source code for the existing boundary conditions
 for examples of how to implement boundary conditions.
@@ -2299,33 +2481,40 @@
 Gravity and friction are always calculated using the elevation and friction quantities,
 but the user may additionally add forcing terms to the list
-\code{domain.forcing\_terms} and have them affect the model.
-
-Currently, predefined forcing terms are:
-\begin{funcdesc}{General\_forcing}{}
-  Module: \module{shallow\_water.shallow\_water\_domain}
-
-  This is a general class to modify any quantity according to a given rate of change.
-  Other specific forcing terms are based on this class but it can be used by itself as well (e.g.\ to modify momentum).
-
-  The General\_forcing class takes as input:
-  \begin{itemize}
-    \item \code{domain}: a reference to the domain being evolved
-    \item \code{quantity\_name}: The name of the quantity that will be affected by this forcing term
-    \item \code{rate}: The rate at which the quantity should change. The parameter \code{rate} can be eithe a constant or a
-                function of time. Positive values indicate increases,
-                negative values indicate decreases.
-                The parametr \code{rate} can be \code{None} at initialisation but must be specified
-                before forcing term is applied (i.e. simulation has started).
-                The default value is 0.0 -- i.e.\ no forcing.
-    \item \code{center, radius}: Optionally restrict forcing to a circle with given center and radius.
-    \item \code{polygon}: Optionally restrict forcing to an area enclosed by given polygon.
-  \end{itemize}
-  Note specifying both center, radius and polygon will cause an exception to be thrown.
-  Moreover, if the specified polygon or circle does not lie fully within the mesh boundary an Exception will be thrown.
-
-  \bigskip
-  Example:
-
-  \begin{verbatim}
+\code{domain.forcing_terms} and have them affect the model.
+
+Currently, predefined forcing terms are: \\
+\begin{classdesc}{General_forcing}{domain,
+                                  quantity_name,
+                                  rate=0.0,
+                                  center=None,
+                                  radius=None,
+                                  polygon=None,
+                                  default_rate=None,
+                                  verbose=False}
+Module: \module{shallow_water.shallow_water_domain}
+
+This is a general class to modify any quantity according to a given rate of change.
+Other specific forcing terms are based on this class but it can be used by itself as well (e.g. to modify momentum).
+
+\code{domain} is the domain being evolved.
+
+\code{quantity_name} is the name of the quantity that will be affected by this forcing term.
+
+\code{rate} is the rate at which the quantity should change. This can be either a constant or a
+function of time. Positive values indicate increases, negative values indicate decreases.
+The parameter \code{rate} can be \code{None} at initialisation but must be specified
+before a forcing term is applied (i.e. simulation has started).
+The default value is 0.0 -- i.e. no forcing.
+
+\code{center} and \code{ radius} optionally restrict forcing to a circle with given center and radius.
+
+\code{polygon} optionally restricts forcing to an area enclosed by the given polygon.
+
+Note: specifying \code{center}, \code{radius} and \code{polygon} will cause an exception to be thrown.
+Moreover, if the specified polygon or circle does not lie fully within the mesh boundary an Exception will be thrown.
+
+Example:
+
+\begin{verbatim}
 P = [[x0, y0], [x1, y0], [x1, y1], [x0, y1]]    # Square polygon
 
@@ -2338,113 +2527,147 @@
 domain.forcing_terms.append(xmom)
 domain.forcing_terms.append(ymom)
-  \end{verbatim}
-
-  Here, \code{f}, \code{g} are assumed to be defined as functions of time providing
-  a time dependent rate of change for xmomentum and ymomentum respectively.
-  P is assumed to be polygon, specified as a list of points.
-\end{funcdesc}
-
-\begin{funcdesc}{Inflow}{}
-  Module: \module{shallow\_water.shallow\_water\_domain}
-
-  This is a general class for inflow and abstraction of water according to a given rate of change.
-  This class will always modify the \code{stage} quantity.
-
-  Inflow is based on the General_forcing class so the functionality is similar.
-
-  The Inflow class takes as input:
-  \begin{itemize}
-    \item \code{domain}: a reference to the domain being evolved
-    \item \code{rate}: The flow rate in $m^3/s$ at which the stage should change. The parameter \code{rate} can be eithe a constant or a
-                function of time. Positive values indicate inflow,
-                negative values indicate outflow.
-
-                Note: The specified flow will be divided by the area of
-                the inflow region and then applied to update the
-                stage quantity.
-    \item \code{center, radius}: Optionally restrict forcing to a circle with given center and radius.
-    \item \code{polygon}: Optionally restrict forcing to an area enclosed by given polygon.
-  \end{itemize}
-
-  \bigskip
-  Example:
-
-  \begin{verbatim}
-hydrograph = Inflow(center=(320, 300), radius=10, rate=file_function('QPMF_Rot_Sub13.tms'))
+\end{verbatim}
+
+Here, \code{f} and \code{g} are assumed to be defined as functions of time providing
+a time dependent rate of change for xmomentum and ymomentum respectively.
+\code{P} is assumed to be the polygon, specified as a list of points.
+\end{classdesc}
+
+\begin{classdesc}{Inflow}{domain,
+                         rate=0.0,
+                         center=None, radius=None,
+                         polygon=None,
+                         default_rate=None,
+                         verbose=False}
+Module: \module{shallow_water.shallow_water_domain}
+
+This is a general class for inflow and abstraction of water according to a given rate of change.
+This class will always modify the \code{stage} quantity.
+
+Inflow is based on the \code{General_forcing} class so the functionality is similar.
+
+\code{domain} is the domain being evolved.
+
+\code{rate} is the flow rate ($m^3/s$) at which the quantity should change. This can be either a constant or a
+function of time. Positive values indicate inflow, negative values indicate outflow.
+Note: The specified flow will be divided by the area of the inflow region and then applied to update the
+stage quantity.
+
+\code{center} and \code{ radius} optionally restrict forcing to a circle with given center and radius.
+
+\code{polygon} optionally restricts forcing to an area enclosed by the given polygon.
+
+Example:
+
+\begin{verbatim}
+hydrograph = Inflow(center=(320, 300), radius=10,
+                    rate=file_function('QPMF_Rot_Sub13.tms'))
 
 domain.forcing_terms.append(hydrograph)
-  \end{verbatim}
-
-  Here, \code{'QPMF_Rot_Sub13.tms'} is assumed to be a NetCDF file in the format \code{tms} defining a timeseries for a hydrograph.
-\end{funcdesc}
-
-\begin{funcdesc}{Rainfall}{}
-  Module: \module{shallow\_water.shallow\_water\_domain}
-
-  This is a general class for implementing rainfall over the domain, possibly restricted to a given circle or polygon.
-  This class will always modify the \code{stage} quantity.
-
-  Rainfall is based on the General_forcing class so the functionality is similar.
-
-  The Rainfall class takes as input:
-  \begin{itemize}
-    \item \code{domain}: a reference to the domain being evolved
-    \item \code{rate}: Total rain rate over the specified domain.
-                  Note: Raingauge Data needs to reflect the time step.
-                  For example: if rain gauge is mm read every \code{dt} seconds, then the input
-                  here is as \code{mm/dt} so 10 mm in 5 minutes becomes
-                  10/(5x60) = 0.0333mm/s.
-
-                  This parameter can be either a constant or a
-                  function of time. Positive values indicate rain being added (or be used for general infiltration),
-                  negative values indicate outflow at the specified rate (presumably this could model evaporation or abstraction).
-    \item \code{center, radius}: Optionally restrict forcing to a circle with given center and radius.
-    \item \code{polygon}: Optionally restrict forcing to an area enclosed by given polygon.
-  \end{itemize}
-
-  \bigskip
-  Example:
-
-  \begin{verbatim}
+\end{verbatim}
+
+Here, \code{'QPMF_Rot_Sub13.tms'} is assumed to be a NetCDF file in the TMS format defining a timeseries for a hydrograph.
+\end{classdesc}
+
+\begin{classdesc}{Rainfall}{domain,
+                           rate=0.0,
+                           center=None,
+                           radius=None,
+                           polygon=None,
+                           default_rate=None,
+                           verbose=False}
+Module: \module{shallow_water.shallow_water_domain}
+
+This is a general class for implementing rainfall over the domain, possibly restricted to a given circle or polygon.
+This class will always modify the \code{stage} quantity.
+
+Rainfall is based on the \code{General_forcing} class so the functionality is similar.
+
+\code{domain} is the domain being evolved.
+
+\code{rate} is the total rain rate over the specified domain.
+Note: Raingauge Data needs to reflect the time step.
+For example, if rain gauge is \code{mm} read every \code{dt} seconds, then the input
+here is as \code{mm/dt} so 10 mm in 5 minutes becomes
+10/(5x60) = 0.0333mm/s.  This parameter can be either a constant or a
+function of time. Positive values indicate rain being added (or be used for general infiltration),
+negative values indicate outflow at the specified rate (presumably this could model evaporation or abstraction).
+
+\code{center} and \code{ radius} optionally restrict forcing to a circle with given center and radius.
+
+\code{polygon} optionally restricts forcing to an area enclosed by the given polygon.
+
+Example:
+
+\begin{verbatim}
 catchmentrainfall = Rainfall(rate=file_function('Q100_2hr_Rain.tms'))
 domain.forcing_terms.append(catchmentrainfall)
-  \end{verbatim}
-
-  Here, \code{'Q100_2hr_Rain.tms'} is assumed to be a NetCDF file in the format \code{tms} defining a timeseries for the rainfall.
-\end{funcdesc}
-
-\begin{funcdesc}{Culvert\_flow}{}
-  Module: \module{culver\_flows.culvert\_class}
-
-  This is a general class for implementing flow through a culvert.
-  This class modifies the quantities \code{stage, xmomentum, ymomentum} in areas at both ends of the culvert.
-
-  The Culvert\_flow forcing term uses \code{Inflow} and {General\_forcing} to update the quantities.
-  The flow drection is determined on-the-fly so openings are referenced simple as opening0 and opening1
-  with either being able to take the role as Inflow and Outflow.
-
-  The Culvert\_flow class takes as input:
-  \begin{itemize}
-    \item \code{domain}: a reference to the domain being evolved
-    \item \code{label}: Short text naming the culvert
-    \item \code{description}: Text describing it
-    \item \code{end_point0}: Coordinates of one opening
-    \item \code{end_point1}: Coordinates of other opening
-    \item \code{width}:
-    \item \code{height}:
-    \item \code{diameter}:
-    \item \code{manning}: Mannings Roughness for Culvert
-    \item \code{invert_level0}: Invert level if not the same as the Elevation on the Domain
-    \item \code{invert_level1}: Invert level if not the same as the Elevation on the Domain
-    \item \code{culvert_routine}: Function specifying the calculation of flow based on energy difference between the two openings (see below)
-    \item \code{number_of_barrels}: Number of identical pipes in the culvert. Default == 1.
-  \end{itemize}
-
-  The user can specify different culvert routines. Hower ANUGA currently provides only one, namely the
-  \code{boyd\_generalised\_culvert\_model} as used in the example below.
-
-  Example:
-
-  \begin{verbatim}
+\end{verbatim}
+
+Here, \code{'Q100_2hr_Rain.tms'} is assumed to be a NetCDF file in the TMS format defining a timeseries for the rainfall.
+\end{classdesc}
+
+\begin{classdesc}{Culvert_flow}{domain,
+                 culvert_description_filename=None,
+                 culvert_routine=None,
+                 end_point0=None,
+                 end_point1=None,
+                 enquiry_point0=None,
+                 enquiry_point1=None,
+                 type='box',
+                 width=None,
+                 height=None,
+                 length=None,
+                 number_of_barrels=1,
+                 trigger_depth=0.01,
+                 manning=None,
+                 sum_loss=None,
+                 use_velocity_head=False,
+                 use_momentum_jet=False,
+                 label=None,
+                 description=None,
+                 update_interval=None,
+                 log_file=False,
+                 discharge_hydrograph=False,
+                 verbose=False}
+Module: \module{culvert_flows.culvert_class}
+
+This is a general class for implementing flow through a culvert.
+This class modifies the quantities \code{stage}, \code{xmomentum} and \code{ymomentum} in areas at both ends of the culvert.
+
+The \code{Culvert_flow} forcing term uses \code{Inflow} and \code{General_forcing} to update the quantities.
+The flow direction is determined on-the-fly so openings are referenced simple as opening0 and opening1
+with either being able to take the role as Inflow or Outflow.
+
+The Culvert_flow class takes as input:
+\begin{itemize}
+  \item \code{domain}: a reference to the domain being evolved
+  \item \code{culvert_description_filename}:
+  \item \code{culvert_routine}:
+  \item \code{end_point0}: Coordinates of one opening
+  \item \code{end_point1}: Coordinates of other opening
+  \item \code{enquiry_point0}:
+  \item \code{enquiry_point1}:
+  \item \code{type}: (default is 'box')
+  \item \code{width}:
+  \item \code{height}:
+  \item \code{length}:
+  \item \code{number_of_barrels}: Number of identical pipes in the culvert (default is 1)
+  \item \code{trigger_depth}: (default is 0.01)
+  \item \code{manning}: Mannings Roughness for Culvert
+  \item \code{sum_loss}:
+  \item \code{use_velocity_head}:
+  \item \code{use_momentum_jet}:
+  \item \code{label}: Short text naming the culvert
+  \item \code{description}: Text describing the culvert
+  \item \code{update_interval}:
+  \item \code{log_file}:
+  \item \code{discharge_hydrograph}:
+\end{itemize}
+
+The user can specify different culvert routines. Hower \anuga currently provides only one, namely the
+\code{boyd_generalised_culvert_model} as used in the example below:
+
+\begin{verbatim}
 from anuga.culvert_flows.culvert_class import Culvert_flow
 from anuga.culvert_flows.culvert_routines import boyd_generalised_culvert_model
@@ -2455,5 +2678,6 @@
                         end_point0=[9.0, 2.5],
                         end_point1=[13.0, 2.5],
-                        width=1.20,height=0.75,
+                        width=1.20,
+                        height=0.75,
                         culvert_routine=boyd_generalised_culvert_model,
                         number_of_barrels=1,
@@ -2473,6 +2697,6 @@
 domain.forcing_terms.append(culvert1)
 domain.forcing_terms.append(culvert2)
-  \end{verbatim}
-\end{funcdesc}
+\end{verbatim}
+\end{classdesc}
 
 
@@ -2480,31 +2704,44 @@
 \label{sec:evolution}
 
-\begin{methoddesc}{evolve}{yieldstep=None, finaltime=None, duration=None, skip_initial_step=False}
-  Module: \module{abstract\_2d\_finite\_volumes.domain}
-
-  This function (a method of \class{domain}) is invoked once all the
-  preliminaries have been completed, and causes the model to progress
-  through successive steps in its evolution, storing results and
-  outputting statistics whenever a user-specified period
-  \code{yieldstep} is completed (generally during this period the
-  model will evolve through several steps internally
-  as the method forces the water speed to be calculated
-  on successive new cells). 
-  
-  The code specified by the user in the block following the evolve statement is only executed once every \code{yieldstep} even though 
-  ANUGA typically will take many more internal steps behind the scenes.
-  
-  The user
-  specifies the total time period over which the evolution is to take
-  place, by specifying values (in seconds) for either \code{duration}
-  or \code{finaltime}, as well as the interval in seconds after which
-  results are to be stored and statistics output.
-
-  You can include \method{evolve} in a statement of the type:
-
-  \begin{verbatim}
+\begin{methoddesc}{\emph{<domain>}.evolve}{yieldstep=None,
+                                           finaltime=None,
+                                           duration=None,
+                                           skip_initial_step=False}
+Module: \module{abstract_2d_finite_volumes.domain}
+
+This method is invoked once all the
+preliminaries have been completed, and causes the model to progress
+through successive steps in its evolution, storing results and
+outputting statistics whenever a user-specified period
+\code{yieldstep} is completed.  Generally during this period the
+model will evolve through several steps internally
+as the method forces the water speed to be calculated
+on successive new cells. 
+
+\code{yieldstep} is the interval in seconds between yields where results are
+stored, statistics written and the domain is inspected or possibly modified.
+If omitted an internal predefined \code{yieldstep} is used.  Internally, smaller
+timesteps may be taken.
+
+\code{duration} is the duration of the simulation in seconds.
+
+\code{finaltime} is the time in seconds where simulation should end. This is currently
+relative time, so it's the same as \code{duration}.  If both \code{duration} and
+\code{finaltime} are given an exception is thrown.
+
+\code{skip_initial_step} is a 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.
+
+The code specified by the user in the block following the evolve statement is
+only executed once every \code{yieldstep} even though 
+\anuga typically will take many more internal steps behind the scenes.
+
+You can include \method{evolve} in a statement of the type:
+
+\begin{verbatim}
 for t in domain.evolve(yieldstep, finaltime):
     <Do something with domain and t>
-  \end{verbatim}
+\end{verbatim}
 \end{methoddesc}
 
@@ -2512,33 +2749,44 @@
 \label{sec:diagnostics}
 
-\begin{funcdesc}{statistics}{}
-  Module: \module{abstract\_2d\_finite\_volumes.domain}
-
-\end{funcdesc}
-
-\begin{funcdesc}{timestepping\_statistics}{}
-  Module: \module{abstract\_2d\_finite\_volumes.domain}
-
-  Returns a string of the following type for each timestep:
-  \code{Time = 0.9000, delta t in [0.00598964, 0.01177388], steps=12
-  (12)}
-
-  Here the numbers in \code{steps=12 (12)} indicate the number of steps taken and
-  the number of first-order steps, respectively.\\
-
-  The optional keyword argument \code{track_speeds=True} will
-  generate a histogram of speeds generated by each triangle. The
-  speeds relate to the size of the timesteps used by ANUGA and
-  this diagnostics may help pinpoint problem areas where excessive speeds
-  are generated.
-\end{funcdesc}
-
-\begin{funcdesc}{boundary\_statistics}{quantities=None, tags=None}
-  Module: \module{abstract\_2d\_finite\_volumes.domain}
-
-  Returns a string of the following type when \code{quantities = 'stage'}
-  and \code{tags = ['top', 'bottom']}:
-
-  \begin{verbatim}
+\begin{methoddesc}{\emph{<domain>}.statistics}{}
+Module: \module{abstract\_2d\_finite\_volumes.domain}
+
+Outputs statistics about the mesh within the \code{Domain}.
+\end{methoddesc}
+
+\begin{methoddesc}{\emph{<domain>}.timestepping_statistics}{track_speeds=False, triangle_id=None}
+Module: \module{abstract_2d_finite_volumes.domain}
+
+Returns a string of the following type for each timestep:\\
+\code{Time = 0.9000, delta t in [0.00598964, 0.01177388], steps=12 (12)}
+
+Here the numbers in \code{steps=12 (12)} indicate the number of steps taken and
+the number of first-order steps, respectively.
+
+The optional keyword argument \code{track_speeds} will
+generate a histogram of speeds generated by each triangle if set to \code{True}. The
+speeds relate to the size of the timesteps used by \anuga and
+this diagnostics may help pinpoint problem areas where excessive speeds
+are generated.
+
+The optional keyword argument \code{triangle_id} can be used to specify a particular
+triangle rather than the one with the largest speed.
+\end{methoddesc}
+
+\begin{methoddesc}{\emph{<domain>}.boundary_statistics}{quantities=None,
+                                                      tags=None}
+Module: \module{abstract_2d_finite_volumes.domain}
+
+Generates output about boundary forcing at each timestep.
+
+\code{quantities} names the quantities to be reported -- may be \code{None},
+a string or a list of strings.
+
+\code{tags} names the tags to be reported -- may be either None, a string or a list of strings.
+
+When \code{quantities = 'stage'} and \code{tags = ['top', 'bottom']}
+will return a string like:
+
+\begin{verbatim}
 Boundary values at time 0.5000:
     top:
@@ -2546,36 +2794,55 @@
     bottom:
         stage in [ -0.27098821,  -0.02499974]
-  \end{verbatim}
-\end{funcdesc}
-
-\begin{funcdesc}{get\_quantity}{name, location='vertices', indices=None}
-  Module: \module{abstract\_2d\_finite\_volumes.domain}
-
-  This function returns a Quantity object Q.
-  Access to its values should be done through Q.get\__values documented on Page \pageref{pg:get values}.
-\end{funcdesc}
-
-\begin{funcdesc}{set\_quantities\_to\_be\_monitored}{}
-  Module: \module{abstract\_2d\_finite\_volumes.domain}
-
-  Selects quantities and derived quantities for which extrema attained at internal timesteps
-  will be collected.
-
-  Information can be tracked in the evolve loop by printing \code{quantity\_statistics} and
-  collected data will be stored in the sww file.
-
-  Optional parameters \code{polygon} and \code{time\_interval} may be specified to restrict the
-  extremum computation.
-\end{funcdesc}
-
-\begin{funcdesc}{quantity\_statistics}{}
-  Module: \module{abstract\_2d\_finite\_volumes.domain}
-
-  Reports on extrema attained by selected quantities.
-
-  Returns a string of the following type for each
-  timestep:
-
-  \begin{verbatim}
+\end{verbatim}
+\end{methoddesc}
+
+\begin{methoddesc}{Q = \emph{<domain>}.get_quantity}{name,
+                                               location='vertices',
+                                               indices=None}
+Module: \module{abstract_2d_finite_volumes.domain}
+
+This function returns a Quantity object Q.
+Access to its values should be done through \code{Q.get_values()} documented on Page \pageref{pg:get values}.
+
+\code{name} is the name of the quantity to retrieve.
+
+\code{location} is
+
+\code{indices} is
+\end{methoddesc}
+
+\begin{methoddesc}{\emph{<domain>}.set_quantities_to_be_monitored}{quantity,
+                                             polygon=None,
+                                             time_interval=None}
+Module: \module{abstract\_2d\_finite\_volumes.domain}
+
+Selects quantities and derived quantities for which extrema attained at internal timesteps
+will be collected.
+
+\code{quantity} specifies the quantity or quantities to be monitored and must be either:
+\begin{itemize}
+  \item the name of a quantity or derived quantity such as 'stage-elevation',
+  \item a list of quantity names, or
+  \item \code{None}.
+\end{itemize}
+
+\code{polygon} can be used to monitor only triangles inside the polygon. Otherwise
+all triangles will be included.
+
+\code{time_interval} will restrict monitoring to time steps in the interval. Otherwise
+all timesteps will be included.
+
+Information can be tracked in the evolve loop by printing \code{quantity_statistics} and
+collected data will be stored in the SWW file.
+\end{methoddesc}
+
+\begin{methoddesc}{\emph{<domain>}.quantity_statistics}{precision='\%.4f'}
+Module: \module{abstract_2d_finite_volumes.domain}
+
+Reports on extrema attained by selected quantities.
+
+Returns a string of the following type for each timestep:
+
+\begin{verbatim}
 Monitored quantities at time 1.0000:
   stage-elevation:
@@ -2591,246 +2858,290 @@
     minimum attained at time = 0.00000000, location = (0.16666667, 0.33333333)
     maximum attained at time = 0.35103646, location = (0.83333333, 0.16666667)
-  \end{verbatim}
-
-  The quantities (and derived quantities) listed here must be selected at model
-  initialisation using the method \code{domain.set_quantities_to_be_monitored}.\\
-
-  The optional keyword argument \code{precision='\%.4f'} will
-  determine the precision used for floating point values in the output.
-  This diagnostics helps track extrema attained by the selected quantities
-  at every internal timestep.
-
-  These values are also stored in the SWW file for post processing.
-\end{funcdesc}
-
-\begin{funcdesc}{get\_values}{location='vertices', indices=None}
-  \label{pg:get values}
-  Module: \module{abstract\_2d\_finite\_volumes.quantity}
-
-  Extract values for quantity as a numeric array.
-
-  \begin{verbatim}
-Inputs:
-  interpolation_points: List of x, y coordinates where value is
-                        sought (using interpolation). If points 
-                        are given, values of location and indices 
-                        are ignored. Assume either absolute UTM
-                        coordinates or geospatial data object.
+\end{verbatim}
+
+The quantities (and derived quantities) listed here must be selected at model
+initialisation time using the method \code{domain.set_quantities_to_be_monitored()}.
+
+The optional keyword argument \code{precision='\%.4f'} will
+determine the precision used for floating point values in the output.
+This diagnostics helps track extrema attained by the selected quantities
+at every internal timestep.
+
+These values are also stored in the SWW file for post-processing.
+\end{methoddesc}
+
+\begin{methoddesc}{\emph{<quantity>}.get_values}{interpolation_points=None,
+                   location='vertices',
+                   indices=None,
+                   use_cache=False,
+                   verbose=False}
+\label{pg:get values}
+Module: \module{abstract_2d_finite_volumes.quantity}
+
+Extract values for quantity as a numeric array.
+
+\code{interpolation_points} is a list of (x, y) coordinates where the value is
+sought (using interpolation). If \code{interpolation_points} is given, values
+for \code{location} and \code{indices} are ignored.
+Assume either an absolute UTM coordinates or geospatial data object.
     
-  location: Where values are to be stored.
-            Permissible options are: vertices, edges, centroids
-            and unique vertices. Default is 'vertices'
-  \end{verbatim}
-
-  The returned values will have the leading dimension equal to length of the indices list or
-  N (all values) if indices is None.
-
-  In case of location == 'centroids' the dimension of returned
-  values will be a list or a numeric 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
-
-  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" 
-        
-  Indices is the set of element ids that the operation applies to.
-
-  The values will be stored in elements following their
-  internal ordering.
-\end{funcdesc}
-
-\begin{funcdesc}{set\_values}{location='vertices', indices=None}
-  Module: \module{abstract\_2d\_finite\_volumes.quantity}
-
-  Assign values to a quantity object.
-  This method works the same way as \code{set\_quantity} except that it doesn't take
-  a quantity name as the first argument. The reason to use \code{set\_values} is for
-  example to assign values to a new quantity that has been created but which is
-  not part of the domain's predefined quantities.
-
-  The method \code{set\_values} is always called by \code{set\_quantity}
-  behind the scenes.
-\end{funcdesc}
-
-\begin{funcdesc}{get\_integral}{}
-  Module: \module{abstract\_2d\_finite\_volumes.quantity}
-
-  Return computed integral over entire domain for this quantity
-\end{funcdesc}
-
-\begin{funcdesc}{get\_maximum\_value}{indices=None}
-  Module: \module{abstract\_2d\_finite\_volumes.quantity}
-
-  Return maximum value of quantity (on centroids)
-
-  Optional argument indices is the set of element ids that
-  the operation applies to. If omitted all elements are considered.
-
-  We do not seek the maximum at vertices as each vertex can
-  have multiple values -- one for each triangle sharing it.
-\end{funcdesc}
-
-\begin{funcdesc}{get\_maximum\_location}{indices=None}
-  Module: \module{abstract\_2d\_finite\_volumes.quantity}
-
-  Return location of maximum value of quantity (on centroids)
-
-  Optional argument indices is the set of element ids that
-  the operation applies to.
-
-  We do not seek the maximum at vertices as each vertex can
-  have multiple values -- one for each triangle sharing it.
-
-  If there are multiple cells with same maximum value, the
-  first cell encountered in the triangle array is returned.
-\end{funcdesc}
-
-\begin{funcdesc}{get\_wet\_elements}{indices=None}
-  Module: \module{shallow\_water.shallow\_water\_domain}
-
-  Return indices for elements where h $>$ minimum_allowed_height
-  Optional argument indices is the set of element ids that the operation applies to.
-\end{funcdesc}
-
-\begin{funcdesc}{get\_maximum\_inundation\_elevation}{indices=None}
-  Module: \module{shallow\_water.shallow\_water\_domain}
-
-  Return highest elevation where h $>$ 0.\\
-  Optional argument indices is the set of element ids that the operation applies to.\\
-
-  Example to find maximum runup elevation:\\
-     z = domain.get_maximum_inundation_elevation()
-\end{funcdesc}
-
-\begin{funcdesc}{get\_maximum\_inundation\_location}{indices=None}
-  Module: \module{shallow\_water.shallow\_water\_domain}
-
-  Return location (x,y) of highest elevation where h $>$ 0.\\
-  Optional argument indices is the set of element ids that the operation applies to.\\
-
-  Example to find maximum runup location:\\
-     x, y = domain.get_maximum_inundation_location()
-\end{funcdesc}
+\code{location} specifies where values are to be stored.
+Permissible options are 'vertices', 'edges', 'centroids' or 'unique vertices'.
+
+The returned values will have the leading dimension equal to length of the \code{indices} list or
+N (all values) if \code{indices} is \code{None}.
+
+If \code{location} is 'centroids' the dimension of returned
+values will be a list or a numeric array of length N, N being
+the number of elements.
+      
+If \code{location} is 'vertices' or 'edges' the dimension of
+returned values will be of dimension \code{Nx3}.
+
+If \code{location} is '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" 
+      
+\code{indices} is the set of element ids that the operation applies to.
+
+The values will be stored in elements following their internal ordering.
+\end{methoddesc}
+
+\begin{methoddesc}{\emph{<quantity>}.set_values}{numeric=None,
+                               quantity=None,
+                               function=None,
+                               geospatial_data=None,
+                               filename=None,
+                               attribute_name=None,
+                               alpha=None,
+                               location='vertices',
+                               polygon=None,
+                               indices=None,
+                               smooth=False,
+                               verbose=False,
+                               use_cache=False}
+Module: \module{abstract_2d_finite_volumes.quantity}
+
+Assign values to a quantity object.
+
+This method works the same way as \code{set_quantity()} except that it doesn't take
+a quantity name as the first argument since it is applied directly to the quantity.
+Use \code{set_values} is used to assign
+values to a new quantity that has been created but which is
+not part of the domain's predefined quantities.
+
+\code{location} is ??
+
+\code{indices} is ??
+
+The method \code{set_values()} is always called by \code{set_quantity()}
+behind the scenes.
+\end{methoddesc}
+
+\begin{methoddesc}{\emph{<quantity>}.get_integral}{}
+Module: \module{abstract_2d_finite_volumes.quantity}
+
+Return the computed integral over the entire domain for the quantity.
+\end{methoddesc}
+
+\begin{methoddesc}{\emph{<quantity>}.get_maximum_value}{indices=None}
+Module: \module{abstract_2d_finite_volumes.quantity}
+
+Return the maximum value of a quantity (on centroids).
+
+\code{indices} is the optional set of element \code{id}s that
+the operation applies to.
+
+We do not seek the maximum at vertices as each vertex can
+have multiple values -- one for each triangle sharing it.
+\end{methoddesc}
+
+\begin{methoddesc}{\emph{<quantity>}.get_maximum_location}{indices=None}
+Module: \module{abstract_2d_finite_volumes.quantity}
+
+Return the location of the maximum value of a quantity (on centroids).
+
+\code{indices} is the optional set of element \code{id}s that
+the operation applies to.
+
+We do not seek the maximum at vertices as each vertex can
+have multiple values -- one for each triangle sharing it.
+
+If there are multiple cells with the same maximum value, the
+first cell encountered in the triangle array is returned.
+\end{methoddesc}
+
+\begin{methoddesc}{\emph{<domain>}.get_wet_elements}{indices=None}
+Module: \module{shallow_water.shallow_water_domain}
+
+Returns the indices for elements where h $>$ minimum_allowed_height
+
+\code{indices} is the optional set of element \code{id}s that
+the operation applies to.
+\end{methoddesc}
+
+\begin{methoddesc}{\emph{<domain>}.get_maximum_inundation_elevation}{indices=None}
+Module: \module{shallow_water.shallow_water_domain}
+
+Return highest elevation where h $>$ 0.
+
+\code{indices} is the optional set of element \code{id}s that
+the operation applies to.
+
+Example to find maximum runup elevation:
+\begin{verbatim}
+z = domain.get_maximum_inundation_elevation()
+\end{verbatim}
+\end{methoddesc}
+
+\begin{methoddesc}{\emph{<domain>}.get_maximum_inundation_location}{indices=None}
+Module: \module{shallow_water.shallow_water_domain}
+
+Return location (x,y) of highest elevation where h $>$ 0.
+
+\code{indices} is the optional set of element \code{id}s that
+the operation applies to.
+
+Example to find maximum runup location:
+\begin{verbatim}
+x, y = domain.get_maximum_inundation_location()
+\end{verbatim}
+\end{methoddesc}
 
 
 \section{Queries of SWW model output files}
-After a model has been run, it is often useful to extract various information from the sww
+After a model has been run, it is often useful to extract various information from the SWW
 output file (see Section \ref{sec:sww format}). This is typically more convenient than using the
 diagnostics described in Section \ref{sec:diagnostics} which rely on the model running -- something
-that can be very time consuming. The sww files are easy and quick to read and offer much information
+that can be very time consuming. The SWW files are easy and quick to read and offer information
 about the model results such as runup heights, time histories of selected quantities,
 flow through cross sections and much more.
 
-\begin{funcdesc}{get\_maximum\_inundation\_elevation}{filename, polygon=None,
-    time_interval=None, verbose=False}
-  Module: \module{shallow\_water.data\_manager}
-
-  Return highest elevation where depth is positive ($h > 0$)
-
-  Example to find maximum runup elevation:
-
-  \begin{verbatim}
-max_runup = get_maximum_inundation_elevation(filename,
-                                             polygon=None,
-                                             time_interval=None,
-                                             verbose=False)
-  \end{verbatim}
-
-  filename is a NetCDF sww file containing ANUGA model output.
-  Optional arguments polygon and time_interval restricts the maximum runup calculation
-  to a points that lie within the specified polygon and time interval.
-
-  If no inundation is found within polygon and time_interval the return value
-  is None signifying "No Runup" or "Everything is dry".
-
-  See doc string for general function get_maximum_inundation_data for details.
+\begin{funcdesc}{elevation = get_maximum_inundation_elevation}{filename,
+                                     polygon=None,
+                                     time_interval=None,
+                                     verbose=False}
+Module: \module{shallow_water.data_manager}
+
+Return the highest elevation where depth is positive ($h > 0$).
+
+\code{filename} is the path to a NetCDF SWW file containing \anuga model output.
+
+\code{polygon} restricts the query to the points that lie within the polygon.
+
+\code {time_interval} restricts the query to within the time interval.
+
+Example to find maximum runup elevation:
+
+\begin{verbatim}
+max_runup = get_maximum_inundation_elevation(filename)
+\end{verbatim}
+
+If no inundation is found (within the \code{polygon} and \code{time_interval}, if specified)
+the return value is \code{None}. This indicates "No Runup" or "Everything is dry".
 \end{funcdesc}
 
-\begin{funcdesc}{get\_maximum\_inundation\_location}{filename, polygon=None,
-    time_interval=None, verbose=False}
-  Module: \module{shallow\_water.data\_manager}
-
-  Return location (x,y) of highest elevation where depth is positive ($h > 0$)
-
-  Example to find maximum runup location:
-
-  \begin{verbatim}
-max_runup_location = get_maximum_inundation_location(filename,
-                                                     polygon=None,
-                                                     time_interval=None,
-                                                     verbose=False)
-  \end{verbatim}
-
-  \code{filename} is a NetCDF sww file containing ANUGA model output.
-  Optional arguments \code{polygon} and \code{time_interval} restricts the maximum runup calculation
-  to a points that lie within the specified polygon and time interval.
-
-  If no inundation is found within polygon and time_interval the return value
-  is None signifying "No Runup" or "Everything is dry".
-
-  See the 'doc' string for general function \code{get_maximum_inundation_data} for details.
+\begin{funcdesc}{(x, y) = get_maximum_inundation_location}{filename,
+                                    polygon=None,
+                                    time_interval=None,
+                                    verbose=False}
+Module: \module{shallow_water.data_manager}
+
+Return location (x,y) of the highest elevation where depth is positive ($h > 0$).
+
+\code{filename} is the path to a NetCDF SWW file containing \anuga model output.
+
+\code{polygon} restricts the query to the points that lie within the polygon.
+
+\code {time_interval} restricts the query to within the time interval.
+
+Example to find maximum runup location:
+
+\begin{verbatim}
+max_runup_location = get_maximum_inundation_location(filename)
+\end{verbatim}
+
+If no inundation is found (within the \code{polygon} and \code{time_interval}, if specified)
+the return value is \code{None}. This indicates "No Runup" or "Everything is dry".
+is \code{None}. This indicates "No Runup" or "Everything is dry".
 \end{funcdesc}
 
-\begin{funcdesc}{sww2timeseries}{swwfiles, gauge_filename, production_dirs, report=None, reportname=None,
-plot_quantity=None, generate_fig=False, surface=None, time_min=None, time_max=None, time_thinning=1,
-time_unit=None, title_on=None, use_cache=False, verbose=False}
-  Module: \module{anuga.abstract\_2d\_finite\_volumes.util}
-
-  Return csv files for the location in the \code{gauge_filename} and can also return plots of them
-
-  See the 'doc' string for general function \code{sww2timeseries} for details.
+\begin{funcdesc}{sww2timeseries}{swwfiles,
+                                 gauge_filename,
+                                 production_dirs,
+                                 report=None,
+                                 reportname=None,
+                                 plot_quantity=None,
+                                 generate_fig=False,
+                                 surface=None,
+                                 time_min=None,
+                                 time_max=None,
+                                 time_thinning=1,
+                                 time_unit=None,
+                                 title_on=None,
+                                 use_cache=False,
+                                 verbose=False}
+Module: \module{abstract_2d_finite_volumes.util}
+
+Read a set of SWW files and plot the time series for the prescribed quantities
+at defined gauge locations and prescribed time range.
+
+\code{swwfiles} is a dictionary of SWW files with keys of \code{label_id}.
+
+\code{gauge_filename} is the path to a file containing gauge data.
+
+THIS NEEDS MORE WORK.  WORK ON FUNCTION __DOC__ STRING, IF NOTHING ELSE!
 \end{funcdesc}
 
-\begin{funcdesc}{get\_flow\_through\_cross\_section}{filename, cross\_section, verbose=False}
-  Module: \module{shallow\_water.data\_manager}
-
-  Obtain flow $[m^3/s]$ perpendicular to specified cross section.
-
-  Inputs:
-  \begin{itemize}
-        \item filename: Name of sww file containing ANUGA model output.
-        \item polyline: Representation of desired cross section -- it may contain multiple
-                        sections allowing for complex shapes. Assume absolute UTM coordinates.
-  \end{itemize}
-
-  Output:
-  \begin{itemize}
-    \item time: All stored times in sww file
-    \item Q: Hydrograph of total flow across given segments for all stored times.
-  \end{itemize}
-
-  The normal flow is computed for each triangle intersected by the polyline and
-  added up.  Multiple segments at different angles are specified the normal flows
-  may partially cancel each other.
-
-  Example to find flow through cross section:
-
-  \begin{verbatim}
+\begin{funcdesc}{(time, Q) = get_flow_through_cross_section}{filename, polyline, verbose=False}
+Module: \module{shallow_water.data_manager}
+
+Obtain flow ($m^3/s$) perpendicular to specified cross section.
+
+\code{filename} is the path to the SWW file.
+
+\code{polyline} is the representation of the desired cross section -- it may contain
+multiple sections allowing for complex shapes. Assumes absolute UTM coordinates.
+
+Returns a tuple \code{time, Q} where:
+
+\code{time} is all the stored times in the SWW file.
+
+\code{Q} is a hydrograph of total flow across the given segments for all stored times.
+
+The normal flow is computed for each triangle intersected by the \code{polyline} and
+added up.  If multiple segments at different angles are specified the normal flows
+may partially cancel each other.
+
+Example to find flow through cross section:
+
+\begin{verbatim}
 cross_section = [[x, 0], [x, width]]
-time, Q = get_flow_through_cross_section(filename,
-                                         cross_section,
-                                         verbose=False)
-  \end{verbatim}
-
-  See the 'doc' string for general function \code{get_maximum_inundation_data} for details.
+time, Q = get_flow_through_cross_section(filename, cross_section)
+\end{verbatim}
 \end{funcdesc}
 
 
 \section{Other}
-\begin{funcdesc}{domain.create\_quantity\_from\_expression}{string}
-  Handy for creating derived quantities on-the-fly as for example:
-
-  \begin{verbatim}
+\begin{methoddesc}{quantity = \emph{<domain>}.create_quantity_from_expression}{string}
+Module: \module{abstract_2d_finite_volumes.domain}
+
+Create a new quantity from other quantities in the domain using an arbitrary expression.
+
+\code{string} is a string containing an arbitrary quantity expression.
+
+Returns the new \code{Quantity} object.
+
+Handy for creating derived quantities on-the-fly:
+
+\begin{verbatim}
 Depth = domain.create_quantity_from_expression('stage-elevation')
 
 exp = '(xmomentum*xmomentum + ymomentum*ymomentum)**0.5'
 Absolute_momentum = domain.create_quantity_from_expression(exp)
-  \end{verbatim}
-
-  %See also \file{Analytical\_solution\_circular\_hydraulic\_jump.py} for an example of use.
-\end{funcdesc}
+\end{verbatim}
+
+%See also \file{Analytical_solution_circular_hydraulic_jump.py} for an example.
+\end{methoddesc}
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -2845,5 +3156,4 @@
 detail in the paragraphs below.
 
-\bigskip
 \begin{center}
 \begin{tabular}{|ll|}  \hline
@@ -2851,7 +3161,7 @@
   \hline\hline
   \code{.sww} & NetCDF format for storing model output with mesh information \code{f(t,x,y)}\\
-  \code{.sts} & NetCDF format for storing model ouput \code{f(t,x,y)} without any mesh information\\
+  \code{.sts} & NetCDF format for storing model ouput \code{f(t,x,y)} without mesh information\\
   \code{.tms} & NetCDF format for storing time series \code{f(t)}\\
-  \code{.csv/.txt} & ASCII format called points csv for storing arbitrary points and associated attributes\\
+  \code{.csv/.txt} & ASCII format for storing arbitrary points and associated attributes\\
   \code{.pts} & NetCDF format for storing arbitrary points and associated attributes\\
   \code{.asc} & ASCII format of regular DEMs as output from ArcView\\
@@ -2863,5 +3173,4 @@
   \code{.nc} & Native ferret NetCDF format\\
   \code{.geo} & Houdinis ASCII geometry format (?) \\  \par \hline
-  %\caption{File formats used by \anuga}
 \end{tabular}
 \end{center}
@@ -2870,9 +3179,11 @@
 formats of files. However, typically, in referring to a format we
 capitalise the extension and omit the initial full stop -- thus, we
-refer, for example, to 'SWW files' or 'PRJ files'.
+refer to 'SWW files' or 'PRJ files', not 'sww files' or '.prj files'.
 
 \bigskip
 
 A typical dataflow can be described as follows:
+
+SOMETHING MISSING HERE!?
 
 \subsection{Manually Created Files}
@@ -2898,5 +3209,5 @@
 \subsection{SWW, STS and TMS Formats}
 \label{sec:sww format}
-The SWW, STS and TMS formats are all NetCDF formats, and are of key importance for \anuga.
+The SWW, STS and TMS formats are all NetCDF formats and are of key importance for \anuga.
 
 An SWW file is used for storing \anuga output and therefore pertains
@@ -2907,22 +3218,19 @@
 \begin{itemize}
   \item \code{x} and \code{y}: coordinates of the points, represented as numeric arrays
-  \item \code{elevation}, a numeric array storing bed-elevations
-  \item \code{volumes}, a list specifying the points at the vertices of each of the triangles
+  \item \code{elevation}: a numeric array storing bed-elevations
+  \item \code{volumes}: a list specifying the points at the vertices of each of the triangles
     % Refer here to the example to be provided in describing the simple example
-  \item \code{time}, a numeric array containing times for model evaluation
+  \item \code{time}: a numeric array containing times for model evaluation
 \end{itemize}
 
-The contents of an SWW file may be viewed using the anuga viewer
-\code{animate}, which creates an on-screen geometric
-representation. See section \ref{sec:animate} (page
-\pageref{sec:animate}) in Appendix \ref{ch:supportingtools} for more
-on \code{animate}.
+The contents of an SWW file may be viewed using the anuga viewer \code{animate},
+which creates an on-screen visialisation.  See section \ref{sec:animate}
+(page \pageref{sec:animate}) in Appendix \ref{ch:supportingtools} for more on \code{animate}.
 
 Alternatively, there are tools, such as \code{ncdump}, that allow
-you to convert an NetCDF file into a readable format such as the
+you to convert a NetCDF file into a readable format such as the
 Class Definition Language (CDL). The following is an excerpt from a
 CDL representation of the output file \file{runup.sww} generated
-from running the simple example \file{runup.py} of
-Chapter \ref{ch:getstarted}:
+from running the simple example \file{runup.py} of Chapter \ref{ch:getstarted}:
 
 %FIXME (Ole): Should put in example with nonzero xllcorner, yllcorner
@@ -2933,11 +3241,11 @@
 \function{file\_function}, described in Chapter \ref{ch:interface}.
 
-An STS file is used for storing a set of points and and associated set of times.
+An STS file is used for storing a set of points and associated times.
 It contains, in addition to dimension information, the following
 variables:
 \begin{itemize}
   \item \code{x} and \code{y}: coordinates of the points, represented as numeric arrays
-  \item \code{permutation}: Original indices of the points as specified by the optional \code{ordering\_file} 
-                            (see the function \code{urs2sts} in Section \ref{sec:basicfileconversions}).
+  \item \code{permutation}: Original indices of the points as specified by the optional \code{ordering_file} 
+                            (see the function \code{urs2sts()} in Section \ref{sec:basicfileconversions})
   \item \code{elevation}: a numeric array storing bed-elevations
     % Refer here to the example to be provided in describing the simple example
@@ -2959,5 +3267,5 @@
 two formats: it can be either a TSH file, which is an ASCII file, or
 an MSH file, which is a NetCDF file. A mesh file can be generated
-from the function \function{create\_mesh\_from\_regions} (see
+from the function \function{create_mesh_from_regions()} (see
 Section \ref{sec:meshgeneration}) and be used to initialise a domain.
 
@@ -2976,5 +3284,5 @@
 
 A mesh file can also contain a georeference, which describes an
-offset to be applied to $x$ and $y$ values -- e.g.\ to the vertices.
+offset to be applied to $x$ and $y$ values -- e.g. to the vertices.
 
 \subsection{Formats for Storing Arbitrary Points and Attributes}
@@ -2991,12 +3299,12 @@
 
 \begin{verbatim}
-x, y, elevation, friction}
-0.6, 0.7, 4.9, 0.3}
-1.9, 2.8, 5, 0.3}
-2.7, 2.4, 5.2, 0.3}
+x, y, elevation, friction
+0.6, 0.7, 4.9, 0.3
+1.9, 2.8, 5.0, 0.3
+2.7, 2.4, 5.2, 0.3
 \end{verbatim}
 
 The delimiter is a comma. The first two columns are assumed to
-be x, y coordinates.
+be $x$ and $y$ coordinates.
 
 A PTS file is a NetCDF representation of the data held in an points CSV
@@ -3038,5 +3346,6 @@
 \label{sec:basicfileconversions}
 
-\begin{funcdesc}{sww2dem}{basename_in, basename_out=None,
+\begin{funcdesc}{sww2dem}{(basename_in,
+            basename_out=None,
             quantity=None,
             timestep=None,
@@ -3049,18 +3358,19 @@
             northing_min=None,
             northing_max=None,
-            expand_search=False,
             verbose=False,
             origin=None,
             datum='WGS84',
-            format='ers'}
-  Module: \module{shallow\_water.data\_manager}
-
-  Takes data from an SWW file \code{basename_in} and converts it to DEM format (ASC or
-  ERS) of a desired grid size \code{cellsize} in metres. The user can select how
-  many the decimal places the output data can be written to using \code{number_of_decimal_places},
-  with the default being 3.
-  The easting and northing values are used if the user wished to determine the output
-  within a specified rectangular area. The \code{reduction} input refers to a function
-  to reduce the quantities over all time step of the SWW file, example, maximum.
+            format='ers',
+            block_size=None}
+Module: \module{shallow_water.data_manager}
+
+Takes data from an SWW file \code{basename_in} and converts it to DEM format (ASC or
+ERS) of a desired grid size \code{cellsize} in metres. The user can select how
+many decimal places the output data is represented with by using \code{number_of_decimal_places},
+with the default being 3.
+
+The $easting$ and $northing$ values are used if the user wishes to determine the output
+within a specified rectangular area. The \code{reduction} input refers to a function
+to reduce the quantities over all time step of the SWW file, e.g. maximum.
 \end{funcdesc}
 
@@ -3302,5 +3612,5 @@
 $z$ the bed elevation. The latter are assumed to be relative to the
 same height datum.
-The conserved quantities tracked by ANUGA are momentum in the
+The conserved quantities tracked by \anuga are momentum in the
 $x$-direction ($\mu = uh$), momentum in the $y$-direction ($\nu = vh$)
 and depth ($h = w-z$).
@@ -3338,5 +3648,5 @@
 \]
 
-ANUGA has a global parameter $H_0$ that controls the minimal depth which
+\anuga has a global parameter $H_0$ that controls the minimal depth which
 is considered in the various equations. This parameter is typically set to
 $10^{-3}$. Setting
@@ -3352,15 +3662,9 @@
 \]
 In general, for multiples of the minimal depth $N H_0$ one obtains
-\begin{equation}
+\[
   \left[ \frac{\mu}{h + h_0/h} \right]_{h = N H_0} =
-  \frac{\mu}{N H_0 + H_0/N} =  
-  \frac{\mu}{h (1 + 1/N^2)}
-  \label{eq:flux limit multiple}
-\end{equation}
+  \frac{\mu}{H_0 (1 + 1/N^2)}
+\]
 which converges quadratically to the true value with the multiple N.
-
-Although this equation can be used for any depth, we have restricted its use to depths less than $10 * H_0$ (or 1 cm) to computational resources. 
-According to Equation \ref{eq:flux limit multiple} this cutoff   
-affects the calculated velocity by less than 1 \%.   
 
 %The developed numerical model has been applied to several test cases
@@ -3485,7 +3789,7 @@
 using some relative time (e.g. 0).
 
-The ANUGA domain has an attribute \code{starttime} which is used in cases where the
+The \anuga domain has an attribute \code{starttime} which is used in cases where the
 simulation should be started later than the beginning of some input data such as those
-obtained from boundaries or forcing functions (hydrographs, file\_boundary etc).
+obtained from boundaries or forcing functions (hydrographs, file_boundary etc).
 
 The \code{domain.startime} may be adjusted in \code{file_boundary} in the case the
@@ -3501,19 +3805,19 @@
 All locations must consequently be specified in Cartesian coordinates
 (eastings, northings) or (x,y) where the unit is metres.
-Alternative projections can be assumed, but ANUGA does have the concept of a UTM zone 
+Alternative projections can be used, but \anuga does have the concept of a UTM zone 
 that must be the same for all coordinates in the model.
 
 \subsection{Internal coordinates}
-It is important to realise that ANUGA for numerical precision uses coordinates that are relative 
+It is important to realise that for numerical precision \anuga uses coordinates that are relative 
 to the lower left node of the rectangle containing the mesh ($x_{\mbox{min}}$, $y_{\mbox{min}}$). 
 This origin is referred to internally as \code{xllcorner}, \code{yllcorner} following the ESRI ASCII grid notation.  
-The sww file format also includes \code{xllcorner}, \code{yllcorner} and any coordinate in the file should be adjusted 
+The SWW file format also includes \code{xllcorner}, \code{yllcorner} and any coordinate in the file should be adjusted 
 by adding this origin. See Section \ref{sec:sww format}.
 
-Throughout the ANUGA interface, functions have optional boolean arguments \code{absolute} which controls
+Throughout the \anuga interface, functions have optional boolean arguments \code{absolute} which controls
 whether spatial data received is using the internal representation (\code{absolute=False}) or the 
-user coordinate set (\code{absolute=True}). See e.g. \code{get_vertex_coordinates} on \pageref{pg:get vertex coordinates}.
-
-DEMs, meshes and boundary conditions can have different origins. However, the internal representation in ANUGA 
+user coordinate set (\code{absolute=True}). See e.g. \code{get_vertex_coordinates()} on \pageref{pg:get vertex coordinates}.
+
+DEMs, meshes and boundary conditions can have different origins. However, the internal representation in \anuga 
 will use the origin of the mesh.
 
@@ -3546,5 +3850,5 @@
 %                   native dem and pts formats%
 %
-%  prepare_pts.py: Convert netcdf output from MOST to an sww file suitable
+%  prepare_pts.py: Convert netcdf output from MOST to an SWW file suitable
 %                  as boundary condition%
 %
@@ -3596,10 +3900,9 @@
 
 If the function definition changes after a result has been cached, this will be
-detected by examining the functions \code{bytecode (co\_code, co\_consts,
-func\_defaults, co\_argcount)} and the function will be recomputed.
+detected by examining the functions bytecode and the function will be recomputed.
 However, caching will not detect changes in modules used by \code{func}.
-In this case cache must be cleared manually.
-
-Options are set by means of the function \code{set\_option(key, value)},
+In this case the cache must be cleared manually.
+
+Options are set by means of the function \code{set_option(key, value)},
 where \code{key} is a key associated with a
 Python dictionary \code{options}. This dictionary stores settings such as the name of
@@ -3611,5 +3914,5 @@
 %Other features include support for compression and a capability to \ldots
 
-\textbf{USAGE:} \nopagebreak
+USAGE: \nopagebreak
 
 \begin{verbatim}
@@ -3619,12 +3922,13 @@
 
 
-\section{ANUGA viewer -- animate}
+\pagebreak
+\section{\anuga viewer -- animate}
 \label{sec:animate}
 
 The output generated by \anuga may be viewed by
-means of the visualisation tool \code{animate}, which takes the
-\code{SWW} file generated by \anuga and creates a visual representation
+means of the visualisation tool \code{animate}, which takes an
+SWW file generated by \anuga and creates a visual representation
 of the data. Examples may be seen in Figures \ref{fig:runupstart}
-and \ref{fig:runup2}. To view an \code{SWW} file with
+and \ref{fig:runup2}. To view an SWW file with
 \code{animate} in the Windows environment, you can simply drag the
 icon representing the file over an icon on the desktop for the
@@ -3635,5 +3939,5 @@
 
 Upon starting the viewer, you will see an interactive moving-picture
-display. You can use the keyboard and the mouse to slow down, speed up or
+display. You can use the keyboard and mouse to slow down, speed up or
 stop the display, change the viewing position or carry out a number
 of other simple operations. Help is also displayed when you press
@@ -3657,8 +3961,8 @@
 % \vfill
 
-The following table describes how to operate \code{animate} from the command line:
-
-Usage: \code{animate [options] swwfile \ldots}\\  \nopagebreak
-Options:\\  \nopagebreak
+The following describes how to operate \code{animate} from the command line:
+
+Usage: \code{animate [options] swwfile \ldots}\\ \nopagebreak
+where: \\ \nopagebreak
 \begin{tabular}{ll}
   \code{--display <type>} & \code{MONITOR | POWERWALL | REALITY\_CENTER |}\\
@@ -3673,28 +3977,13 @@
   \code{-alphamax <float 0-1>} & Maximum transparency clamp value\\
   \code{-alphamin <float 0-1>} & Transparency value at \code{hmin}\\
-\end{tabular}
-
-\begin{tabular}{ll}
   \code{-cullangle <float angle 0-90>} & Cull triangles steeper than this value\\
   \code{-help} & Display this information\\
   \code{-hmax <float>} & Height above which transparency is set to
                                      \code{alphamax}\\
-\end{tabular}
-
-\begin{tabular}{ll}
   \code{-hmin <float>} & Height below which transparency is set to
                                      zero\\
-\end{tabular}
-
-\begin{tabular}{ll}
   \code{-lightpos <float>,<float>,<float>} & $x,y,z$ of bedslope directional light ($z$ is
                                      up, default is overhead)\\
-\end{tabular}
-
-\begin{tabular}{ll}
   \code{-loop}  & Repeated (looped) playback of \code{.swm} files\\
-\end{tabular}
-
-\begin{tabular}{ll}
   \code{-movie <dirname>} & Save numbered images to named directory and
                                      quit\\
@@ -3714,124 +4003,239 @@
 \refmodindex{utilities.polygon}
 
-\begin{classdesc}{Polygon\_function}{regions, default=0.0, geo_reference=None}
-  Module: \code{utilities.polygon}
-
-  Creates a callable object that returns one of a specified list of values when
-  evaluated at a point \code{x, y}, depending on which polygon, from a specified list of polygons, the
-  point belongs to. The parameter \code{regions} is a list of pairs
-  \code{(P, v)}, where each \code{P} is a polygon and each \code{v}
-  is either a constant value or a function of coordinates \code{x}
-  and \code{y}, specifying the return value for a point inside \code{P}. The
-  optional parameter \code{default} may be used to specify a value
-  (or a function)
-  for a point not lying inside any of the specified polygons. When a
-  point lies in more than one polygon, the return value is taken to
-  be the value for whichever of these polygon appears later in the
-  list.
-  %FIXME (Howard): CAN x, y BE VECTORS?
-  The optional parameter geo\_reference refers to the status of points
-  that are passed into the function. Typically they will be relative to
-  some origin. In ANUGA, a typical call will take the form:
-
-  \begin{verbatim}
+\begin{classdesc}{<callable> = Polygon_function}{regions,
+                                    default=0.0,
+                                    geo_reference=None}
+Module: \code{utilities.polygon}
+
+Creates a callable object that returns one of a specified list of values when
+evaluated at a point \code{(x, y)}, depending on which polygon, from a specified list of polygons, the
+point belongs to.
+
+\code{regions} is a list of pairs
+\code{(P, v)}, where each \code{P} is a polygon and each \code{v}
+is either a constant value or a function of coordinates \code{x}
+and \code{y}, specifying the return value for a point inside \code{P}.
+
+\code{default} may be used to specify a value (or a function)
+for a point not lying inside any of the specified polygons.
+
+When a point lies in more than one polygon, the return value is taken to
+be the value for whichever of these polygon appears later in the list.
+
+%FIXME (Howard): CAN x, y BE VECTORS?
+
+\code{geo_reference} refers to the status of points
+that are passed into the function. Typically they will be relative to
+some origin.
+
+Typical usage may take the form:
+
+\begin{verbatim}
 set_quantity('elevation',
              Polygon_function([(P1, v1), (P2, v2)],
                               default=v3,
                               geo_reference=domain.geo_reference))
-  \end{verbatim}
+\end{verbatim}
 \end{classdesc}
 
-\begin{funcdesc}{read\_polygon}{filename}
-  Module: \code{utilities.polygon}
-
-  Reads the specified file and returns a polygon. Each
-  line of the file must contain exactly two numbers, separated by a comma, which are interpreted
-  as coordinates of one vertex of the polygon.
+\begin{funcdesc}{<polygon> = read_polygon}{filename, split=','}
+Module: \code{utilities.polygon}
+
+Reads the specified file and returns a polygon.
+Each line of the file must contain exactly two numbers, separated by a delimiter, which are interpreted
+as coordinates of one vertex of the polygon.
+
+\code{filename} is the path to the file to read.
+
+\code{split} sets the delimiter character between the numbers on one line of
+the file.  If not specified, the delimiter is the ',' character.
 \end{funcdesc}
 
-\begin{funcdesc}{populate\_polygon}{polygon, number_of_points, seed = None, exclude = None}
-  Module: \code{utilities.polygon}
-
-  Populates the interior of the specified polygon with the specified number of points,
-  selected by means of a uniform distribution function.
+\label{ref:function_populate_polygon}
+\begin{funcdesc}{populate_polygon}{polygon, number_of_points, seed=None, exclude=None}
+Module: \code{utilities.polygon}
+
+Populates the interior of the specified polygon with the specified number of points,
+selected by means of a uniform distribution function.
+
+\code{polygon} is the polygon to populate.
+
+\code{number_of_points} is the (optional) number of points.
+
+\code{seed}is the optional seed for random number generator.
+
+\code{exclude} is a list of polygons (inside main polygon) from where points should be excluded.
 \end{funcdesc}
 
-\begin{funcdesc}{point\_in\_polygon}{polygon, delta=1e-8}
-  Module: \code{utilities.polygon}
-
-  Returns a point inside the specified polygon and close to the edge. The distance between
-  the returned point and the nearest point of the polygon is less than $\sqrt{2}$ times the
-  second argument \code{delta}, which is taken as $10^{-8}$ by default.
+\label{ref:function_point_in_polygon}
+\begin{funcdesc}{<point> = point_in_polygon}{polygon, delta=1e-8}
+Module: \code{utilities.polygon}
+
+Returns a point inside the specified polygon and close to the edge. The distance between
+the returned point and the nearest point of the polygon is less than $\sqrt{2}$ times the
+second argument \code{delta}, which is taken as $10^{-8}$ by default.
 \end{funcdesc}
 
-\begin{funcdesc}{inside\_polygon}{points, polygon, closed = True, verbose = False}
-  Module: \code{utilities.polygon}
-
-  Used to test whether the members of a list of points
-  are inside the specified polygon. Returns a numeric
-  array comprising the indices of the points in the list that lie inside the polygon.
-  (If none of the points are inside, returns \code{zeros((0,), 'l')}.)
-  Points on the edges of the polygon are regarded as inside if
-  \code{closed} is set to \code{True} or omitted; otherwise they are regarded as outside.
+\label{ref:function_inside_polygon}
+\begin{funcdesc}{<array> = inside_polygon}{points, polygon, closed=True, verbose=False}
+Module: \code{utilities.polygon}
+
+Get a set of points that lie inside a given polygon.
+
+\code{points} is the list of points to test.
+
+\code{polygon} is the polygon to test the points against.
+
+\code{closed} specifies if points on the polygon edge are considered to be inside
+or outside the polygon -- \code{True} means they are inside.
+
+Returns a numeric array comprising the indices of the points in the list
+that lie inside the polygon.  If none of the points are inside, returns
+\code{zeros((0,), 'l') (ie, an empty numeric array)}.
+
+Compare with \code{outside_polygon()}, page \pageref{ref:function_outside_polygon}.
 \end{funcdesc}
 
-\begin{funcdesc}{outside\_polygon}{points, polygon, closed = True, verbose = False}
-  Module: \code{utilities.polygon}
-
-  Exactly like \code{inside\_polygon}, but with the words 'inside' and 'outside' interchanged.
+\label{ref:function_outside_polygon}
+\begin{funcdesc}{<array> = outside_polygon}{points, polygon, closed=True, verbose=False}
+Module: \code{utilities.polygon}
+
+Get a set of points that lie outside a given polygon.
+
+\code{points} is the list of points to test.
+
+\code{polygon} is the polygon to test the points against.
+
+\code{closed} specifies if points on the polygon edge are considered to be outside
+or inside the polygon -- \code{True} means they are outside.
+
+Returns a numeric array comprising the indices of the points in the list
+that lie outside the polygon.  If none of the points are outside, returns
+\code{zeros((0,), 'l')} (ie, an empty numeric array).
+
+Compare with \code{inside_polygon()}, page \pageref{ref:function_inside_polygon}.
 \end{funcdesc}
 
-\begin{funcdesc}{is\_inside\_polygon}{point, polygon, closed=True, verbose=False}
-  Module: \code{utilities.polygon}
-
-  Returns \code{True} if \code{point} is inside \code{polygon} or
-  \code{False} otherwise. Points on the edges of the polygon are regarded as inside if
-  \code{closed} is set to \code{True} or omitted; otherwise they are regarded as outside.
+\label{ref:function_is_inside_polygon}
+\begin{funcdesc}{<boolean> = is_inside_polygon}{point, polygon, closed=True, verbose=False}
+Module: \code{utilities.polygon}
+
+Determines if a single point is inside a polygon.
+
+\code{point} is the point to test.
+
+\code{polygon} is the polygon to test \code{point} against.
+
+\code{closed} is a flag that forces the function to consider a point on the polygon
+edge to be inside or outside -- if \code{True} a point on the edge is considered inside the
+polygon.
+
+Returns \code{True} if \code{point} is inside \code{polygon}.
+
+Compare with \code{inside_polygon()}, page \pageref{ref:function_inside_polygon}.
 \end{funcdesc}
 
-\begin{funcdesc}{is\_outside\_polygon}{point, polygon, closed=True, verbose=False}
-  Module: \code{utilities.polygon}
-
-  Exactly like \code{is\_outside\_polygon}, but with the words 'inside' and 'outside' interchanged.
+\label{ref:function_is_outside_polygon}
+\begin{funcdesc}{<boolean> = is_outside_polygon}{point, polygon, closed=True, verbose=False,
+%                                                 points_geo_ref=None, polygon_geo_ref=None
+                                                }
+Module: \code{utilities.polygon}
+
+Determines if a single point is outside a polygon.
+
+\code{point} is the point to test.
+
+\code{polygon} is the polygon to test \code{point} against.
+
+\code{closed}  is a flag that forces the function to consider a point on the polygon
+edge to be outside or inside -- if \code{True} a point on the edge is considered inside the
+polygon.
+
+%\code{points_geo_ref} is ??
+%
+%\code{polygon_geo_ref} is ??
+
+Compare with \code{outside_polygon()}, page \pageref{ref:function_outside_polygon}.
 \end{funcdesc}
 
-\begin{funcdesc}{point\_on\_line}{x, y, x0, y0, x1, y1}
-  Module: \code{utilities.polygon}
-
-  Returns \code{True} or \code{False}, depending on whether the point with coordinates
-  \code{x, y} is on the line passing through the points with coordinates \code{x0, y0}
-  and \code{x1, y1} (extended if necessary at either end).
+\label{ref:function_point_on_line}
+\begin{funcdesc}{<boolean> = point_on_line}{point, line, rtol=1.0e-5, atol=1.0e-8}
+Module: \code{utilities.polygon}
+
+Determine if a point is on a line to some tolerance.  The line is considered to
+extend past its end-points.
+
+\code{point} is the point to test ([$x$, $y$]).
+
+\code{line} is the line to test \code{point} against ([[$x1$,$y1$], [$x2$,$y2$]]).
+
+\code{rtol} is the relative tolerance to use when testing for coincidence.
+
+\code{atol} is the absolute tolerance to use when testing for coincidence.
+
+Returns \code{True} if the point is on the line, else \code{False}.
 \end{funcdesc}
 
-\begin{funcdesc}{separate\_points\_by\_polygon}{points, polygon, closed = True, verbose = False}
-    \indexedcode{separate\_points\_by\_polygon}
-  Module: \code{utilities.polygon}
-
+\label{ref:function_separate_points_by_polygon}
+\begin{funcdesc}{(indices, count) = separate_points_by_polygon}{points, polygon,
+                                             closed=True,
+                                             check_input=True,
+                                             verbose=False}
+\indexedcode{separate_points_by_polygon}
+Module: \code{utilities.polygon}
+
+Separate a set of points into points that are inside and outside a polygon.
+
+\code{points} is a list of points to separate.
+
+\code{polygon} is the polygon used to separate the points.
+
+\code{closed} determines whether points on the polygon edge should be
+regarded as inside or outside the polygon.  \code{True} means they are inside.
+
+\code{check_input} specifies whether the input parameters are checked -- \code{True}
+means check the input parameters.
+
+The function returns a tuple \code{(indices, count)} where \code{indices} is a list of 
+point $indices$ from the input \code{points} list, with the indices of points inside the
+polygon at the left and indices of points outside the polygon listed at the right.  The 
+\code{count} value is the count (from the left) of the indices of the points $inside$ the
+polygon.
 \end{funcdesc}
 
-\begin{funcdesc}{polygon\_area}{polygon}
-  Module: \code{utilities.polygon}
-
-  Returns area of arbitrary polygon (reference http://mathworld.wolfram.com/PolygonArea.html)
+\begin{funcdesc}{<area> = polygon_area}{polygon}
+Module: \code{utilities.polygon}
+
+Returns area of an arbitrary polygon (reference http://mathworld.wolfram.com/PolygonArea.html).
 \end{funcdesc}
 
-\begin{funcdesc}{plot\_polygons}{polygons, style, figname, verbose = False}
-  Module: \code{utilities.polygon}
-
-  Plots each polygon contained in input polygon list, e.g.
-  \code{polygons = [poly1, poly2, poly3]} where \code{poly1 = [[x11,y11],[x12,y12],[x13,y13]]}
-  etc.  Each polygon can be closed for plotting purposes by assigning the style type to each
-  polygon in the list, e.g. \code{style = ['line','line','line']}. The default will be a line
-  type when \code{style = None}.
-  The subsequent plot will be saved to \code{figname} or defaulted to \code{test_image.png}.
-  The function returns a list containing the minimum and maximum of \code{x} and \code{y},
-  i.e. \code{[x_{min}, x_{max}, y_{min}, y_{max}]}.
+\begin{funcdesc}{[$x_{min}$, $x_{max}$, $y_{min}$, $y_{max}$] = plot_polygons}
+                               {polygons_points, style=None,
+                                figname=None, label=None, verbose=False}
+Module: \code{utilities.polygon}
+
+Plot a list of polygons to a file.
+
+\code{polygons_points} is a list of polygons to plot.
+
+\code{style} is a list of style strings to be applied to the corresponding polygon
+in \code{polygons_points}. A polygon can be closed for plotting purposes by assigning
+the style string 'line' to it in the appropriate place in the \code{style} list.
+The default style is 'line'.
+
+\code{figname} is the path to the file to save the plot in.  If not specified, use
+\file{test_image.png}.
+
+The function returns a list containing the minimum and maximum of the points in all the
+input polygons, i.e. \code{[$x_{min}$, $x_{max}$, $y_{min}$, $y_{max}$]}.
 \end{funcdesc}
 
 
-\section{coordinate\_transforms}
-
-
-\section{geospatial\_data}
+\pagebreak
+\section{coordinate_transforms}
+
+
+\pagebreak
+\section{geospatial_data}
 \label{sec:geospatial}
 
@@ -3843,163 +4247,283 @@
 %\refmodindex{geospatial_data}
 
-\begin{classdesc}{Geospatial\_data}
-  {data_points=None,
-    attributes=None,
-    geo_reference=None,
-    default_attribute_name=None,
-    file_name=None}
-  Module: \code{geospatial\_data}
-
-  This class is used to store a set of data points and associated
-  attributes, allowing these to be manipulated by methods defined for
-  the class.
-
-  The data points are specified either by reading them from a NetCDF
-  or CSV file, identified through the parameter \code{file\_name}, or
-  by providing their \code{x}- and \code{y}-coordinates in metres,
-  either as a sequence of 2-tuples of floats or as an $M \times 2$
-  numeric array of floats, where $M$ is the number of points.
-  Coordinates are interpreted relative to the origin specified by the
-  object \code{geo\_reference}, which contains data indicating the UTM
-  zone, easting and northing. If \code{geo\_reference} is not
-  specified, a default is used.
-
-  Attributes are specified through the parameter \code{attributes},
-  set either to a list or array of length $M$ or to a dictionary whose
-  keys are the attribute names and whose values are lists or arrays of
-  length $M$. One of the attributes may be specified as the default
-  attribute, by assigning its name to \code{default\_attribute\_name}.
-  If no value is specified, the default attribute is taken to be the
-  first one.
-
-  Note that the Geospatial\_data object currently reads entire datasets
-  into memory i.e.\ no memomry blocking takes place.
-  For this we refer to the set\_quantity method which will read .pts and .csv
-  files into \anuga using memory blocking allowing large files to be used.
+\begin{classdesc}{Geospatial_data}
+                {data_points=None,
+                 attributes=None,
+                 geo_reference=None,
+                 default_attribute_name=None,
+                 file_name=None,
+                 latitudes=None,
+                 longitudes=None,
+                 points_are_lats_longs=False,
+                 max_read_lines=None,
+                 load_file_now=True,
+                 verbose=False}
+Module: \code{geospatial_data.geospatial_data}
+
+This class is used to store a set of data points and associated
+attributes, allowing these to be manipulated by methods defined for
+the class.
+
+The data points are specified either by reading them from a NetCDF
+or CSV file, identified through the parameter \code{file_name}, or
+by providing their \code{x}- and \code{y}-coordinates in metres,
+either as a sequence of 2-tuples of floats or as an $M \times 2$
+numeric array of floats, where $M$ is the number of points.
+
+Coordinates are interpreted relative to the origin specified by the
+object \code{geo_reference}, which contains data indicating the UTM
+zone, easting and northing. If \code{geo_reference} is not
+specified, a default is used.
+
+Attributes are specified through the parameter \code{attributes},
+set either to a list or array of length $M$ or to a dictionary whose
+keys are the attribute names and whose values are lists or arrays of
+length $M$. One of the attributes may be specified as the default
+attribute, by assigning its name to \code{default_attribute_name}.
+If no value is specified, the default attribute is taken to be the
+first one.
+
+Note that the \code{Geospatial_data} object currently reads entire datasets
+into memory i.e. no memomry blocking takes place.
+For this we refer to the \code{set_quantity()} method which will read PTS and CSV
+files into \anuga using memory blocking allowing large files to be used.
 \end{classdesc}
 
-\begin{methoddesc}{import\_points\_file}{delimiter=None, verbose=False}
-
-\end{methoddesc}
-
-\begin{methoddesc}{export\_points\_file}{ofile, absolute=False}
-
-\end{methoddesc}
-
-\begin{methoddesc}{get\_data\_points}{absolute=True, as\_lat\_long=False}
-    If \code{as\_lat\_long} is\code{True} the point information
-    returned will be in Latitudes and Longitudes.
-\end{methoddesc}
-
-\begin{methoddesc}{set\_attributes}{attributes}
-
-\end{methoddesc}
-
-\begin{methoddesc}{get\_attributes}{attribute_name=None}
-
-\end{methoddesc}
-
-\begin{methoddesc}{get\_all\_attributes}{}
-
-\end{methoddesc}
-
-\begin{methoddesc}{set\_default\_attribute\_name}{default_attribute_name}
-
-\end{methoddesc}
-
-\begin{methoddesc}{set\_geo\_reference}{geo_reference}
-
-\end{methoddesc}
-
-\begin{methoddesc}{add}{}
-
-\end{methoddesc}
-
-\begin{methoddesc}{clip}{}
-  Clip geospatial data by a polygon
-
-  Inputs are \code{polygon} which is either a list of points, an Nx2 array or
-  a Geospatial data object and \code{closed}(optional) which determines
-  whether points on boundary should be regarded as belonging to the polygon
-  (\code{closed=True}) or not (\code{closed=False}).
-  Default is \code{closed=True}.
-
-  Returns new Geospatial data object representing points inside specified polygon.
-\end{methoddesc}
-
-\begin{methoddesc}{clip_outside}{}
-  Clip geospatial data by a polygon
-
-  Inputs are \code{polygon} which is either a list of points, an Nx2 array or
-  a Geospatial data object and \code{closed}(optional) which determines
-  whether points on boundary should be regarded as belonging to the polygon
-  (\code{closed=True}) or not (\code{closed=False}).
-  Default is \code{closed=True}.
-
-  Returns new Geospatial data object representing points \emph{out}side specified polygon.
-\end{methoddesc}
-
-\begin{methoddesc}{split}{factor=0.5, seed_num=None, verbose=False}
-  Returns two geospatial_data object, first is the size of the 'factor'
-  smaller the original and the second is the remainder. The two
-  new object are disjoin set of each other.
-
-  Points of the two new geospatial_data object are selected RANDOMLY.
-
-  Input -- the (\code{factor}) which to split the object, if 0.1 then 10% of the
-  together object will be returned
-
-  Output -- two geospatial_data objects that are disjoint sets of the original
-\end{methoddesc}
-
-\begin{methoddesc}{find_optimal_smoothing_parameter}{data_file, alpha_list=None,
-                   mesh_file=None, boundary_poly=None, mesh_resolution=100000,
-                   north_boundary=None, south_boundary=None, east_boundary=None,
-                   west_boundary=None, plot_name='all_alphas', split_factor=0.1,
-                   seed_num=None, cache=False, verbose=False}
-  Removes a small random sample of points from 'data_file'. Creates
-  models from resulting points in 'data_file' with different alpha values from 'alpha_list' and cross validates
-  the predicted value to the previously removed point data. Returns the
-  alpha value which has the smallest covariance.
-
-  data_file: must not contain points outside the boundaries defined
-  and it either a pts, txt or csv file.
-
-  alpha_list: the alpha values to test in a single list
-
-  mesh_file: name of the created mesh file or if passed in will read it.
-  NOTE, if there is a mesh file mesh_resolution, north_boundary, south... etc will be ignored.
-
-  mesh_resolution: the maximum area size for a triangle
-
-  north_boundary... west_boundary: the value of the boundary
-
-  plot_name: the name for the plot contain the results
-
-  seed_num: the seed to the random number generator
-
-  USAGE:
-  convariance_value, alpha = find_optimal_smoothing_parameter(data_file=fileName,
-                                               alpha_list=[0.0001, 0.01, 1],
-                                               mesh_file=None,
-                                               mesh_resolution=3,
-                                               north_boundary=5,
-                                               south_boundary=-5,
-                                               east_boundary=5,
-                                               west_boundary=-5,
-                                               plot_name='all_alphas',
-                                               seed_num=100000,
-                                               verbose=False)
-
-  OUTPUT: returns the minumum normalised covalance calculate AND the
-  alpha that created it. PLUS writes a plot of the results
-
-  NOTE: code will not work if the data_file extent is greater than the
-  boundary_polygon or any of the boundaries, eg north_boundary...west_boundary
-\end{methoddesc}
-
-
+\begin{methoddesc}{\emph{<Geospatial_data>}.import_points_file}
+        {file_name, delimiter=None, verbose=False}
+Module: \code{geospatial_data.geospatial_data}
+
+Import a TXT, CSV or PTS points data file into a code{Geospatial_data} object.
+
+\code{file_name} is the path to a TXT, CSV or PTS points data file.
+
+\code{delimiter} is currently unused.
+\end{methoddesc}
+
+\begin{methoddesc}{\emph{<Geospatial_data>}.export_points_file}{file_name, absolute=True,
+                                       as_lat_long=False, isSouthHemisphere=True}
+Module: \code{geospatial_data.geospatial_data}
+
+Export a CSV or PTS points data file from a \code{Geospatial_data} object.
+
+\code{file_name} is the path to the CSV or PTS points file to write.
+
+\code{absolute} determines if the exported data is absolute or relative to the 
+\code{Geospatial_data} object geo_reference.  If \code{True} the exported 
+data is absolute.
+
+\code{as_lat_long} exports the points data as latitudes and longitudes if \code{True}.
+
+\code{isSouthHemisphere} has effect only if \code{as_lat_long} is \code{True} and causes 
+latitude/longitude values to be for the southern (\code{True}) or northern hemispheres
+(\code{False}).
+\end{methoddesc}
+
+\begin{methoddesc}{points = \emph{<Geospatial_data>}.get_data_points}
+        {absolute=True, geo_reference=None,
+         as_lat_long=False, isSouthHemisphere=True}
+Module: \code{geospatial_data.geospatial_data}
+
+Get the coordinates for all the data points as an $N \times 2$ array.
+
+\code{absolute} determines if the exported data is absolute or relative to the 
+\code{Geospatial_data} object geo_reference.  If \code{True} the exported 
+data is absolute.
+
+\code{geo_reference} is the geo_reference the points are relative to, if supplied.
+
+\code{as_lat_long} exports the points data as latitudes and longitudes if \code{True}.
+
+\code{isSouthHemisphere} has effect only if \code{as_lat_long} is \code{True} and causes 
+latitude/longitude values to be for the southern (\code{True}) or northern hemispheres
+(\code{False}).
+\end{methoddesc}
+
+\begin{methoddesc}{\emph{<Geospatial_data>}.set_attributes}{attributes}
+Module: \code{geospatial_data.geospatial_data}
+
+Set the attributes for a \code{Geospatial_data} object.
+
+\code{attributes} is the new value for the object's attributes.  May be a dictionary or \code{None}.
+\end{methoddesc}
+
+\begin{methoddesc}{atributes = \emph{<Geospatial_data>}.get_attributes}{attribute_name=None}
+Module: \code{geospatial_data.geospatial_data}
+
+Get a named attribute from a \code{Geospatial_data} object.
+
+\code{attribute_name} is the name of the desired attribute.  If \code{None}, return 
+the default attribute.
+\end{methoddesc}
+
+\begin{methoddesc}{\emph{<Geospatial_data>}.get_all_attributes}{}
+Module: \code{geospatial_data.geospatial_data}
+
+Get all attributes of a \code{Geospatial_data} object.
+
+Returns \code{None} or the attributes dictionary (which may be empty).
+\end{methoddesc}
+
+\begin{methoddesc}{\emph{<Geospatial_data>}.set_default_attribute_name}{default_attribute_name}
+Module: \code{geospatial_data.geospatial_data}
+
+Set the default attribute name of a \code{Geospatial_data} object.
+
+\code{default_attribute_name} is the new default attribute name.
+\end{methoddesc}
+
+\begin{methoddesc}{\emph{<Geospatial_data>}.set_geo_reference}{geo_reference}
+Module: \code{geospatial_data.geospatial_data}
+
+Set the internal geo_reference of a \code{Geospatial_data} object.
+
+\code{geo_reference} is the new internal geo_reference for the object.
+If \code{None} will use the default geo_reference.
+
+If the \code{Geospatial_data} object already has an internal geo_reference
+then the points data will be changed to use the new geo_reference.
+\end{methoddesc}
+
+\begin{methoddesc}{\emph{<Geospatial_data>}.__add__}{other}
+Module: \code{geospatial_data.geospatial_data}
+
+The \code{__add__()} method is defined so it is possible to add two
+\code{Geospatial_data} objects.
+\end{methoddesc}
+
+\label{ref:function_clip}
+\begin{methoddesc}{geospatial = \emph{<Geospatial_data>}.clip}{polygon, closed=True, verbose=False}
+Module: \code{geospatial_data.geospatial_data}
+
+Clip a \code{Geospatial_data} object with a polygon.
+
+\code{polygon} is the polygon to clip the \code{Geospatial_data} object with.
+This may be a list of points, an $N \times 2$ array or a \code{Geospatial_data}
+object.
+
+\code{closed} determines whether points on the \code{polygon} edge are inside (\code{True})
+or outside (\code{False}) the polygon.
+
+Returns a new \code{Geospatial_data} object representing points inside the
+
+Compare with \code{clip_outside()}, page \pageref{ref:function_clip_outside}.
+specified polygon.
+\end{methoddesc}
+
+\label{ref:function_clip_outside}
+\begin{methoddesc}{geospatial = \emph{<Geospatial_data>}.clip_outside}
+        {polygon, closed=True, verbose=False}
+Module: \code{geospatial_data.geospatial_data}
+
+Clip a \code{Geospatial_data} object with a polygon.
+
+\code{polygon} is the polygon to clip the \code{Geospatial_data} object with.
+This may be a list of points, an $N \times 2$ array or a \code{Geospatial_data}
+object.
+
+\code{closed} determines whether points on the \code{polygon} edge are inside (\code{True})
+or outside (\code{False}) the polygon.
+
+Returns a new \code{Geospatial_data} object representing points outside the
+specified polygon.
+
+Compare with \code{clip()}, page \pageref{ref:function_clip}.
+\end{methoddesc}
+
+\begin{methoddesc}{(g1, g2) = \emph{<Geospatial_data>}.split}
+        {factor=0.5, seed_num=None, verbose=False}
+Module: \code{geospatial_data.geospatial_data}
+
+Split a \code{Geospatial_data} object into two objects of predetermined ratios.
+
+\code{factor} is the ratio of the size of the first returned object to the
+original object.  If '0.5' is supplied, the two resulting objects will be
+of equal size.
+
+\code{seed_num}, if supplied, will be the random number generator seed used for
+the split.
+
+Points of the two new geospatial_data object are selected RANDOMLY.
+
+Returns two geospatial_data objects that are disjoint sets of the original.
+\end{methoddesc}
+
+\subsection{Miscellaneous Functions}
+
+The functions here are not \code{Geospatial_data} object methods, but are used with them.
+
+\begin{methoddesc}{X = find_optimal_smoothing_parameter}
+                  {data_file,
+                   alpha_list=None,
+                   mesh_file=None,
+                   boundary_poly=None,
+                   mesh_resolution=100000,
+                   north_boundary=None,
+                   south_boundary=None,
+                   east_boundary=None,
+                   west_boundary=None,
+                   plot_name='all_alphas',
+                   split_factor=0.1,
+                   seed_num=None,
+                   cache=False,
+                   verbose=False}
+Module: \code{geospatial_data.geospatial_data}
+
+Calculate the minimum covariance from a set of points in a file.  It does this
+by removing a small random sample of points from \code{data_file} and creating
+models with different alpha values from \code{alpha_list} and cross validates
+the predicted value to the previously removed point data. Returns the
+alpha value which has the smallest covariance.
+
+\code{data_file} is the input data file and must not contain points outside
+the boundaries defined and is either a PTS, TXT or CSV file.
+
+\code{alpha_list} is the list of alpha values to use.
+
+\code{mesh_file} is the path to a mesh file to create (if supplied).
+If \code{None} a mesh file will be created (named \file{temp.msh}).
+NOTE: if there is a \code{mesh_resolution} defined or any boundaries are defined,
+any input \code{mesh_file} value is ignored.
+
+\code{mesh_resolution} is the maximum area size for a triangle.
+
+\code{north_boundary}\\
+\code{south_boundary}\\
+\code{east_boundary}\\
+\code{west_boundary} are the boundary values to use.
+
+\code{plot_name} is the path name of the plot file to write.
+
+\code{seed_num} is the random number generator seed to use.
+
+The function returns a tuple \code{(min_covar, alpha)} where \code{min_covar} is
+the minumum normalised covariance and \code{alpha} is the alpha value that 
+created it.  A plot file is also written.
+
+This is an example of function usage: \nopagebreak
+
+\begin{verbatim}
+convariance_value, alpha = $\backslash$
+        find_optimal_smoothing_parameter(data_file=fileName,
+                                         alpha_list=[0.0001, 0.01, 1],
+                                         mesh_file=None,
+                                         mesh_resolution=3,
+                                         north_boundary=5,
+                                         south_boundary=-5,
+                                         east_boundary=5,
+                                         west_boundary=-5,
+                                         plot_name='all_alphas',
+                                         seed_num=100000,
+                                         verbose=False)
+\end{verbatim}
+
+NOTE: The function will not work if the \code{data_file} extent is greater than the
+\code{boundary_poly} polygon or any of the boundaries, e.g. \code{north_boundary}, etc.
+\end{methoddesc}
+
+
+\pagebreak
 \section{Graphical Mesh Generator GUI}
-The program \code{graphical_mesh_generator.py} in the pmesh module
+The program \code{graphical_mesh_generator.py} in the \code{pmesh} module
 allows the user to set up the mesh of the problem interactively.
 It can be used to build the outline of a mesh or to visualise a mesh
@@ -4007,5 +4531,5 @@
 
 Graphical Mesh Generator will let the user select various modes. The
-current allowable modes are vertex, segment, hole or region.  The mode
+current allowable modes are $vertex$, $segment$, $hole$ or $region$.  The mode
 describes what sort of object is added or selected in response to
 mouse clicks.  When changing modes any prior selected objects become
@@ -4017,5 +4541,6 @@
 
 
-\section{alpha\_shape}
+\pagebreak
+\section{class Alpha_Shape}
 \emph{Alpha shapes} are used to generate close-fitting boundaries
 around sets of points. The alpha shape algorithm produces a shape
@@ -4039,61 +4564,72 @@
 with instances of this class.
 
-\begin{classdesc}{Alpha\_Shape}{points, alpha=None}
-  Module: \code{alpha\_shape}
-
-  To instantiate this class the user supplies the points from which
-  the alpha shape is to be created (in the form of a list of 2-tuples
-  \code{[[x1, y1],[x2, y2]}\ldots\code{]}, assigned to the parameter
-  \code{points}) and, optionally, a value for the parameter
-  \code{alpha}. The alpha shape is then computed and the user can then
-  retrieve details of the boundary through the attributes defined for
-  the class.
+\label{class:alpha_shape}
+\begin{classdesc}{Alpha_Shape}{points, alpha=None}
+Module: \code{alpha_shape}
+
+Instantiate an instance of the \code{Alpha_Shape} class.
+
+\code{points} is an $N \times 2$ list of points (\code{[[x1, y1],[x2, y2]}\ldots\code{]}).
+
+\code{alpha} is the 'fitting' parameter.
 \end{classdesc}
 
-\begin{funcdesc}{alpha\_shape\_via\_files}{point_file, boundary_file, alpha=None}
-  Module: \code{alpha\_shape}
-
-  This function reads points from the specified point file
-  \code{point\_file}, computes the associated alpha shape (either
-  using the specified value for \code{alpha} or, if no value is
-  specified, automatically setting it to an optimal value) and outputs
-  the boundary to a file named \code{boundary\_file}. This output file
-  lists the coordinates \code{x, y} of each point in the boundary,
-  using one line per point.
+\begin{funcdesc}{alpha_shape_via_files}{point_file, boundary_file, alpha= None}
+Module: \code{alpha_shape}
+
+This function reads points from the specified point file
+\code{point_file}, computes the associated alpha shape (either
+using the specified value for \code{alpha} or, if no value is
+specified, automatically setting it to an optimal value) and outputs
+the boundary to a file named \code{boundary_file}. This output file
+lists the coordinates \code{(x, y)} of each point in the boundary,
+using one line per point.
 \end{funcdesc}
 
-\begin{methoddesc}{set\_boundary\_type}{self,raw_boundary=True,
-                          remove_holes=False,
-                          smooth_indents=False,
-                          expand_pinch=False,
-                          boundary_points_fraction=0.2}
-  Module: \code{alpha\_shape},  Class: \class{Alpha\_Shape}
-
-  This function sets flags that govern the operation of the algorithm
-  that computes the boundary, as follows:
-
-  \code{raw\_boundary = True} returns raw boundary, i.e. the regular edges of the
-                  alpha shape.\\
-  \code{remove\_holes = True} removes small holes ('small' is defined by
-  \code{boundary\_points\_fraction})\\
-  \code{smooth\_indents = True} removes sharp triangular indents in
-  boundary\\
-  \code{expand\_pinch = True} tests for pinch-off and
-  corrects -- preventing a boundary vertex from having more than two edges.
-\end{methoddesc}
-
-\begin{methoddesc}{get\_boundary}{}
-  Module: \code{alpha\_shape},  Class: \class{Alpha\_Shape}
-
-  Returns a list of tuples representing the boundary of the alpha
-  shape. Each tuple represents a segment in the boundary by providing
-  the indices of its two endpoints.
-\end{methoddesc}
-
-\begin{methoddesc}{write\_boundary}{file_name}
-  Module: \code{alpha\_shape},  Class: \class{Alpha\_Shape}
-
-  Writes the list of 2-tuples returned by \code{get\_boundary} to the
-  file \code{file\_name}, using one line per tuple.
+\label{ref:method_set_boundary_type}
+\begin{methoddesc}{\emph{<Alpha_shape>}.set_boundary_type}
+        {raw_boundary=True,
+         remove_holes=False,
+         smooth_indents=False,
+         expand_pinch=False,
+         boundary_points_fraction=0.2}
+Module: \code{alpha_shape}
+
+This function sets internal state that controls how the \code{Alpha_shape}
+boundary is presented or exported.
+
+\code{raw_boundary} sets the type to $raw$ if \code{True},
+i.e. the regular edges of the alpha shape.
+
+\code{remove_holes}, if \code{True} removes small holes ('small' is defined by
+\code{boundary_points_fraction}).
+
+\code{smooth_indents}, if \code{True} removes sharp triangular indents in
+the boundary.
+
+\code{expand_pinch}, if \code{True} tests for pinch-off and
+corrects -- preventing a boundary vertex from having more than two edges.
+\end{methoddesc}
+
+\label{ref:method_get_boundary}
+\begin{methoddesc}{boundary = \emph{<Alpha_shape>}.get_boundary}{}
+Module: \code{alpha_shape}
+
+Returns a list of tuples representing the boundary of the alpha
+shape. Each tuple represents a segment in the boundary by providing
+the indices of its two endpoints.
+
+See \code{set_boundary_type()}, page \pageref{ref:method_set_boundary_type}.
+\end{methoddesc}
+
+\label{ref:method_write_boundary}
+\begin{methoddesc}{\emph{<Alpha_shape>}.write_boundary}{file_name}
+Module: \code{alpha_shape}
+
+Writes the list of 2-tuples returned by \code{get_boundary()} to the
+file \code{file_name}, using one line per tuple.
+
+See \code{set_boundary_type()}, page \pageref{ref:method_set_boundary_type}. \\
+See \code{get_boundary()}, page \pageref{ref:method_get_boundary}.
 \end{methoddesc}
 
@@ -4153,60 +4689,115 @@
 more to come very soon
 
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-\chapter{Modules available in \anuga}
-
-
-\section{\module{abstract\_2d\_finite\_volumes.general\_mesh} }
-\declaremodule[generalmesh]{}{general\_mesh}
-\label{mod:generalmesh}
-
-
-\section{\module{abstract\_2d\_finite\_volumes.neighbour\_mesh} }
-\declaremodule[neighbourmesh]{}{neighbour\_mesh}
-\label{mod:neighbourmesh}
-
-
-\section{\module{abstract\_2d\_finite\_volumes.domain}}
-Generic module for 2D triangular domains for finite-volume computations of conservation laws
-\declaremodule{}{domain}
-\label{mod:domain}
-
-
-\section{\module{abstract\_2d\_finite\_volumes.quantity}}
-\declaremodule{}{quantity}
-\label{mod:quantity}
-
-\begin{verbatim}
-Class Quantity - Implements values at each triangular element
-
-To create:
-
-   Quantity(domain, vertex_values)
-
-   domain: Associated domain structure. Required.
-
-   vertex_values: N x 3 array of values at each vertex for each element.
-                  Default None
-
-   If vertex_values are None Create array of zeros compatible with domain.
-   Otherwise check that it is compatible with dimenions of domain.
-   Otherwise raise an exception
-\end{verbatim}
-
-
-\section{\module{shallow\_water}}
-
-2D triangular domains for finite-volume
-computations of the shallow water wave equation.
-This module contains a specialisation of class Domain from module domain.py consisting of methods specific to the Shallow Water
-Wave Equation
-
-\declaremodule[shallowwater]{}{shallow\_water}
-\label{mod:shallowwater}
+%?% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%?% 
+%?% \chapter{Modules available in \anuga}
+%?% 
+%?% 
+%?% %abstract_2d_finite_volumes
+%?% 
+%?% \section{\module{abstract_2d_finite_volumes.domain}}
+%?% Generic module for 2D triangular domains for finite-volume computations of conservation laws
+%?% \declaremodule[domain]{}{domain}
+%?% \label{mod:domain}
+%?% 
+%?% 
+%?% \section{\module{abstract_2d_finite_volumes.ermapper_grids}}
+%?% \declaremodule[ermappergrids]{}{ermapper_grids}
+%?% \label{mod:ermapper_grids}
+%?% 
+%?% 
+%?% \section{\module{abstract_2d_finite_volumes.general_mesh} }
+%?% \declaremodule[generalmesh]{}{general_mesh}
+%?% \label{mod:general_mesh}
+%?% 
+%?% 
+%?% \section{\module{abstract_2d_finite_volumes.generic_boundary_conditions} }
+%?% \declaremodule[genericboundaryconditions]{}{generic_boundary_conditions}
+%?% \label{mod:generic_boundary_conditions}
+%?% 
+%?% 
+%?% \section{\module{abstract_2d_finite_volumes.mesh_factory} }
+%?% \declaremodule[meshfactory]{}{mesh_factory}
+%?% \label{mod:mesh_factory}
+%?% 
+%?% 
+%?% \section{\module{abstract_2d_finite_volumes.mesh_factory} }
+%?% \declaremodule[meshfactory]{}{mesh_factory}
+%?% \label{mod:mesh_factory}
+%?% 
+%?% 
+%?% \section{\module{abstract_2d_finite_volumes.neighbour_mesh} }
+%?% \declaremodule[neighbourmesh]{}{neighbour_mesh}
+%?% \label{mod:neighbour_mesh}
+%?% 
+%?% 
+%?% \section{\module{abstract_2d_finite_volumes.pmesh2domain} }
+%?% \declaremodule[pmesh2domain]{}{pmesh2domain}
+%?% \label{mod:pmesh2domain}
+%?% 
+%?% 
+%?% \section{\module{abstract_2d_finite_volumes.quantity}}
+%?% \declaremodule{}{quantity}
+%?% \label{mod:quantity}
+%?% 
+%?% \begin{verbatim}
+%?% Class Quantity - Implements values at each triangular element
+%?% 
+%?% To create:
+%?% 
+%?%    Quantity(domain, vertex_values)
+%?% 
+%?%    domain: Associated domain structure. Required.
+%?% 
+%?%    vertex_values: Nx3 array of values at each vertex for each element.
+%?%                   Default None
+%?% 
+%?%    If vertex_values are None Create array of zeros compatible with domain.
+%?%    Otherwise check that it is compatible with dimenions of domain.
+%?%    Otherwise raise an exception
+%?% \end{verbatim}
+%?% 
+%?% 
+%?% \section{\module{abstract_2d_finite_volumes.region} }
+%?% \declaremodule[region]{}{region}
+%?% \label{mod:region}
+%?% 
+%?% 
+%?% \section{\module{abstract_2d_finite_volumes.util} }
+%?% \declaremodule[util]{}{util}
+%?% \label{mod:util}
+%?% 
+%?% 
+%?% advection
+%?% alpha_shape
+%?% caching
+%?% coordinate_transforms
+%?% culvert_flows
+%?% damage_modelling
+%?% euler
+%?% fit_interpolate
+%?% geospatial_data
+%?% lib
+%?% load_mesh
+%?% mesh_engine
+%?% pmesh
+%?% SConstruct
+%?% shallow_water
+%?% utilities
+%?% 
+%?% 
+%?% \section{\module{shallow\_water}}
+%?% 
+%?% 2D triangular domains for finite-volume
+%?% computations of the shallow water wave equation.
+%?% This module contains a specialisation of class Domain from module domain.py consisting of methods specific to the Shallow Water
+%?% Wave Equation
+%?% 
+%?% \declaremodule[shallowwater]{}{shallow\_water}
+%?% \label{mod:shallowwater}
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-\chapter{ANUGA Full-scale Validations}
+\chapter{\anuga Full-scale Validations}
 
 
@@ -4232,5 +4823,5 @@
 
 You may download the validation data files by hand and then run the validation.
-Just go to the ANUGA Sourceforge project download page at
+Just go to the \anuga Sourceforge project download page at
 \module{http://sourceforge.net/project/showfiles.php?group_id=172848} and select
 the \module{validation_data} package, \module{patong-1.0} release.  You need the
@@ -4335,9 +4926,9 @@
 %
 
-%begin{latexonly}
-%\renewcommand{\indexname}{Module Index}
-%end{latexonly}
-\input{mod\jobname.ind}        % Module Index
-%
+%?% %begin{latexonly}
+%?% %\renewcommand{\indexname}{Module Index}
+%?% %end{latexonly}
+%?% \input{mod\jobname.ind}        % Module Index
+
 %begin{latexonly}
 %\renewcommand{\indexname}{Index}
Index: /anuga_core/documentation/user_manual/python.sty
===================================================================
--- /anuga_core/documentation/user_manual/python.sty	(revision 7133)
+++ /anuga_core/documentation/user_manual/python.sty	(revision 7134)
@@ -611,4 +611,6 @@
 }
 \newenvironment{cfuncdesc}[4][\py@badkey]{
+  \linebreak
+  \linebreak
   \begin{fulllineitems}
     \cfuncline{#2}{#3}{#4}
@@ -616,4 +618,6 @@
       \emph{Return value: \textbf{#1}.}\\
     \fi
+  \linebreak
+  \linebreak
 }{\end{fulllineitems}}
 
