Index: /anuga_core/documentation/user_manual/anuga_user_manual.tex
===================================================================
--- /anuga_core/documentation/user_manual/anuga_user_manual.tex	(revision 7085)
+++ /anuga_core/documentation/user_manual/anuga_user_manual.tex	(revision 7086)
@@ -29,4 +29,5 @@
 \usepackage[english]{babel}
 \usepackage{datetime}
+\usepackage[hang,small,bf]{caption}
 
 \input{definitions}
@@ -158,5 +159,5 @@
 
 The latest installation instructions may be found at:
-\url{http://datamining.anu.edu.au/~ole/anuga/user_manual/anuga_installation_guide.pdf}.
+\url{http://datamining.anu.edu.au/\~{}ole/anuga/user_manual/anuga_installation_guide.pdf}.
 
 \section{Audience}
@@ -320,6 +321,4 @@
 \subsection{The Code}
 
-%FIXME: we are using the \code function here.
-%This should be used wherever possible
 For reference we include below the complete code listing for
 \file{runup.py}. Subsequent paragraphs provide a
@@ -327,4 +326,5 @@
 significance.
 
+\label{ref:runup_py_code}
 \verbatiminput{demos/runup.py}
 
@@ -384,5 +384,5 @@
 \end{verbatim}
 
-In addition, the following statement now specifies that the
+In addition, the following statement could be used to state that
 quantities \code{stage}, \code{xmomentum} and \code{ymomentum} are
 to be stored:
@@ -391,4 +391,6 @@
 domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
 \end{verbatim}
+
+However, this is not necessary, as the above quantities are always stored by default.
 
 \subsection{Initial Conditions}
@@ -416,5 +418,5 @@
 
 \begin{verbatim}
-def f(x,y):
+def topography(x, y):
     return -x/2
 \end{verbatim}
@@ -425,13 +427,13 @@
 the \code{y} direction.
 
-Once the function \function{f} is specified, the quantity
+Once the function \function{topography} is specified, the quantity
 \code{elevation} is assigned through the simple statement:
 
 \begin{verbatim}
-domain.set_quantity('elevation', f)
+domain.set_quantity('elevation', topography)
 \end{verbatim}
 
 NOTE: If using function to set \code{elevation} it must be vector
-compatible. For example square root will not work.
+compatible. For example, using square root will not work.
 
 \subsubsection{Friction}
@@ -465,5 +467,5 @@
 
 which specifies that the surface level is set to a height of $-0.4$,
-i.e. 0.4 units (m) below the zero level.
+i.e. 0.4 units (metres) below the zero level.
 
 Although it is not necessary for this example, it may be useful to
@@ -496,10 +498,8 @@
 \begin{verbatim}
 Br = Reflective_boundary(domain)
-
 Bt = Transmissive_boundary(domain)
-
 Bd = Dirichlet_boundary([0.2, 0.0, 0.0])
-
-Bw = Time_boundary(domain=domain, f=lambda t: [(0.1*sin(t*2*pi)-0.3), 0.0, 0.0])
+Bw = Time_boundary(domain=domain,
+                   f=lambda t: [(0.1*sin(t*2*pi)-0.3)*exp(-t), 0.0, 0.0])
 \end{verbatim}
 
@@ -511,7 +511,7 @@
 follows:
 \begin{itemize}
-    \item \textbf{Reflective boundary}\label{def:reflective boundary} Returns same \code{stage} as
-          as present in its neighbour volume but momentum vector
-          reversed 180 degrees (reflected).
+    \item \textbf{Reflective boundary}\label{def:reflective boundary}
+          Returns same \code{stage} as in its neighbour volume but momentum
+          vector reversed 180 degrees (reflected).
           Specific to the shallow water equation as it works with the
           momentum quantities assumed to be the second and third conserved
@@ -615,5 +615,5 @@
 
 \begin{verbatim}
-for t in domain.evolve(yieldstep=0.1, duration=4.0):
+for t in domain.evolve(yieldstep=0.1, duration=10.0):
     print domain.timestepping_statistics()
 \end{verbatim}
@@ -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} visualisation package 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
@@ -661,7 +661,5 @@
 on the previously dry bed.
 
-All figures are screenshots from an interactive animation tool called \code{animate}
-which is part of \anuga and distributed as in the package anuga\_viewer.
-\code{animate} is described in more detail is Section \ref{sec:animate}.
+\code{animate} is described in more detail in Section \ref{sec:animate}.
 
 \begin{figure}[htp]
@@ -713,5 +711,7 @@
 
 Here is the code for the first version of the channel flow \file{channel1.py}:
+
 \verbatiminput{demos/channel1.py}
+
 In discussing the details of this example, we follow the outline
 given above, discussing each major step of the code in turn.
@@ -747,5 +747,5 @@
 as we will later in this example, one merely decreases the values of \code{dx} and \code{dy}.
 
-The rest of this script is as in the previous example.
+The rest of this script is similar to the previous example on page \pageref{ref:runup_py_code}.
 % except for an application of the 'expression' form of \code{set\_quantity} where we use
 % the value of \code{elevation} to define the (dry) initial condition for \code{stage}:
@@ -784,5 +784,5 @@
 
 \begin{verbatim}
-for t in domain.evolve(yieldstep = 0.2, finaltime = 40.0):
+for t in domain.evolve(yieldstep=0.2, finaltime=40.0):
     domain.write_time()
 
@@ -793,5 +793,5 @@
 
 \label{sec:change boundary code}
-The \code{if} statement in the timestepping loop (evolve) gets the quantity
+The \code{if} statement in the timestepping loop (\code{evolve}) gets the quantity
 \code{stage} and obtains the interpolated value at the point (10m,
 2.5m) which is on the right boundary. If the stage exceeds 0m a
@@ -923,18 +923,10 @@
 using a mesh-generator.
 
-In its simplest form, the mesh-generator creates the mesh within a single
+The mesh-generator creates the mesh within a single
 polygon whose vertices are at geographical locations specified by
 the user. The user specifies the \emph{resolution} -- that is, the
-maximal area of a triangle used for triangulationR -- and a triangular
+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.
-%Figure
-%\ref{fig:pentagon} shows a simple example of this, in which the
-%triangulation is carried out within a pentagon.
-
-%\begin{figure}[htp]
-%  \caption{Mesh points are created inside the polygon}
-  %\label{fig:pentagon}
-%\end{figure}
 
 Boundary tags are not restricted to \code{'left'}, \code{'bottom'},
@@ -951,9 +943,4 @@
 polygons in which no triangulation is required.
 
-%\begin{figure}[htp]
-%  \caption{Interior meshes with individual resolution}
-%  \label{fig:interior meshes}
-%\end{figure}
-
 In its general form, the mesh-generator takes for its input a bounding
 polygon and (optionally) a list of interior polygons. The user
@@ -963,7 +950,8 @@
 
 The function used to implement this process is
-\function{create\_mesh\_from\_regions}. Its arguments include the
+\function{create\_domain\_from\_regions} which creates a Domain object as
+well as a mesh file.  Its arguments include the
 bounding polygon and its resolution, a list of boundary tags, and a
-list of pairs \code{[polygon, resolution]}, specifying the interior
+list of pairs \code{[polygon, resolution]} specifying the interior
 polygons and their resolutions.
 
@@ -1033,14 +1021,14 @@
 \begin{verbatim}
 remainder_res = 10000000
-create_mesh_from_regions(project.bounding_polygon,
-                         boundary_tags={'top': [0],
-                                        'ocean_east': [1],
-                                        'bottom': [2],
-                                        'onshore': [3]},
-                         maximum_triangle_area=remainder_res,
-                         filename=meshname,
-                         interior_regions=interior_regions,
-                         use_cache=True,
-                         verbose=True)
+domain = create_domain_from_regions(project.bounding_polygon,
+                                    boundary_tags={'top': [0],
+                                                   'ocean_east': [1],
+                                                   'bottom': [2],
+                                                   'onshore': [3]},
+                                    maximum_triangle_area=project.default_res,
+                                    mesh_filename=project.meshname,
+                                    interior_regions=project.interior_regions,
+                                    use_cache=True,
+                                    verbose=True)
 \end{verbatim}
 
@@ -1056,6 +1044,6 @@
 (Here, the values associated with each boundary tag are one-element lists, but they can have as many indices as there are edges)
 If polygons intersect, or edges coincide (or are even very close) the resolution may be undefined in some regions.
-Use the underlying mesh interface for such cases.
-See Section \ref{sec:mesh interface}.
+Use the underlying mesh interface for such cases
+(see Chapter \ref{sec:mesh interface}).
 If a segment is omitted in the tags definition an Exception is raised.
 
@@ -1067,36 +1055,14 @@
 \subsection{Initialising the Domain}
 
-As with \file{runup.py}, once we have created the mesh, the next
-step is to create the data structure \code{domain}. We did this for
-\file{runup.py} by inputting lists of points and triangles and
-specifying the boundary tags directly. However, in the present case,
-we use a method that works directly with the mesh file
-\code{meshname}, as follows:
+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.
+
+The following statements specify a basename and data directory, and
+sets a minimum storable height, which helps with visualisation.
 
 \begin{verbatim}
-domain = Domain(meshname, use_cache=True, verbose=True)
-\end{verbatim}
-
-Providing a filename instead of the lists used in \file{runup.py}
-above causes \code{Domain} to convert a mesh file \code{meshname}
-into an instance of \code{Domain}, allowing us to use methods like
-\method{set\_quantity} to set quantities and to apply other
-operations.
-
-%(In principle, the
-%second argument of \function{pmesh\_to\_domain\_instance} can be any
-%subclass of \class{Domain}, but for applications involving the
-%shallow-water wave equation, the second argument of
-%\function{pmesh\_to\_domain\_instance} can always be set simply to
-%\class{Domain}.)
-
-The following statements specify a basename and data directory, and
-identify quantities to be stored. For the first two, values are
-taken from \file{project.py}.
-
-\begin{verbatim}
-domain.set_name(project.basename)
-domain.set_datadir(project.outputdir)
-domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
+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}
 
@@ -1111,18 +1077,25 @@
 \subsubsection{Stage}
 
-For the scenario we are modelling in this case, we use a callable
-object \code{tsunami\_source}, assigned by means of a function
-\function{slide\_tsunami}. This is similar to how we set elevation in
-\file{runup.py} using a function -- however, in this case the
-function is both more complex and more interesting.
-
-The function returns the water displacement for all \code{x} and
-\code{y} in the domain. The water displacement is a double Gaussian
-function that depends on the characteristics of the slide (length,
-width, thickness, slope, etc), its location (origin) and the depth at that
-location. For this example, we choose to apply the slide function
-at a specified time into the simulation. {\bf Note, the parameters used
-in this example have been deliberately chosen to generate a suitably
-large amplitude tsunami which would inundate the Cairns region.}
+The stage is initially set to 0.0 by the following statements:
+
+\begin{verbatim}
+tide = 0.0
+domain.set_quantity('stage', tide)
+\end{verbatim}
+
+%For the scenario we are modelling in this case, we use a callable
+%object \code{tsunami_source}, assigned by means of a function
+%\function{slide\_tsunami}. This is similar to how we set elevation in
+%\file{runup.py} using a function -- however, in this case the
+%function is both more complex and more interesting.
+
+%The function returns the water displacement for all \code{x} and
+%\code{y} in the domain. The water displacement is a double Gaussian
+%function that depends on the characteristics of the slide (length,
+%width, thickness, slope, etc), its location (origin) and the depth at that
+%location. For this example, we choose to apply the slide function
+%at a specified time into the simulation. {\bf Note, the parameters used
+%in this example have been deliberately chosen to generate a suitably
+%large amplitude tsunami which would inundate the Cairns region.}
 
 \subsubsection{Friction}
@@ -1140,30 +1113,9 @@
 \begin{verbatim}
 domain.set_quantity('elevation',
-                    filename=project.dem_name+'.pts',
+                    filename=project.demname + '.pts',
                     use_cache=True,
-                    verbose=True)
+                    verbose=True,
+                    alpha=0.1)
 \end{verbatim}
-
-%However, before this step can be executed, some preliminary steps
-%are needed to prepare the file from which the data is taken. Two
-%source files are used for this data -- their names are specified in
-%the file \file{project.py}, in the variables \code{coarsedemname}
-%and \code{finedemname}. They contain 'coarse' and 'fine' data,
-%respectively -- that is, data sampled at widely spaced points over a
-%large region and data sampled at closely spaced points over a
-%smaller subregion. The data in these files is combined through the
-%statement
-
-%\begin{verbatim}
-%combine_rectangular_points_files(project.finedemname + '.pts',
-%                                 project.coarsedemname + '.pts',
-%                                 project.combineddemname + '.pts')
-%\end{verbatim}
-%
-%The effect of this is simply to combine the datasets by eliminating
-%any coarse data associated with points inside the smaller region
-%common to both datasets. The name to be assigned to the resulting
-%dataset is also derived from the name stored in the variable
-%\code{combinedname} in the file \file{project.py}.
 
 \subsection{Boundary Conditions}\index{boundary conditions}
@@ -1174,28 +1126,48 @@
 boundary tag names introduced when we established the mesh. In place of the four
 boundary types introduced for \file{runup.py}, we use the reflective
-boundary for each of the eight tagged segments defined by \code{create_mesh_from_regions}:
+boundary for each of the tagged segments defined by \code{create_domain_from_regions}:
 
 \begin{verbatim}
-Bd = Dirichlet_boundary([0.0,0.0,0.0])
-domain.set_boundary({'ocean_east': Bd, 'bottom': Bd, 'onshore': Bd, 'top': Bd})
+Bd = Dirichlet_boundary([tide,0,0]) # Mean water level
+Bs = Transmissive_stage_zero_momentum_boundary(domain) # Neutral boundary
+
+if project.scenario == 'fixed_wave':
+    # Huge 50m wave starting after 60 seconds and lasting 1 hour.
+    Bw = Time_boundary(domain=domain,
+                       function=lambda t: [(60<t<3660)*50, 0, 0])
+    domain.set_boundary({'ocean_east': Bw,
+                         'bottom': Bs,
+                         'onshore': Bd,
+                         'top': Bs})
+
+if project.scenario == 'slide':
+    # Boundary conditions for slide scenario
+    domain.set_boundary({'ocean_east': Bd,
+                         'bottom': Bd,
+                         'onshore': Bd,
+                         'top': Bd})
 \end{verbatim}
 
+Note that we use different boundary conditions depending on the \code{scenario}
+defined in \file{project.py}.
+
 \subsection{Evolution}
 
 With the basics established, the running of the 'evolve' step is
-very similar to the corresponding step in \file{runup.py}. For the slide
-scenario, the simulation is run for 5000 seconds with the output stored every ten seconds.
-For this example, we choose to apply the slide at 60 seconds into the simulation:
+very similar to the corresponding step in \file{runup.py}, except we have different \code{evolve}
+loops for the two scenarios.
+
+For the slide scenario, the simulation is run for an intial 60 seconds, at which time
+the slide occurs.  We use the function \function{tsunami_source} to adjust \code{stage}
+values.  We then run the simulation until 5000 seconds with the output stored
+every ten seconds.
 
 \begin{verbatim}
-import time
-
-t0 = time.time()
-
-for t in domain.evolve(yieldstep=10, finaltime=60):
-        domain.write_time()
-        domain.write_boundary_statistics(tags = 'ocean_east')
-
-    # add slide
+if project.scenario == 'slide':
+    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')
     if allclose(t, 60):
@@ -1205,8 +1177,19 @@
 
     for t in domain.evolve(yieldstep=10, finaltime=5000,
-                           skip_initial_ste=True):
-        domain.write_time()
-
-    domain.write_boundary_statistics(tags = 'ocean_east')
+                           skip_initial_step=True):
+        print domain.timestepping_statistics()
+        print domain.boundary_statistics(tags='ocean_east')
+
+if project.scenario == 'fixed_wave':
+    # Save every two mins leading up to wave approaching land
+    for t in domain.evolve(yieldstep=120, finaltime=5000):
+        print domain.timestepping_statistics()
+        print domain.boundary_statistics(tags='ocean_east')
+
+    # Save every 30 secs as wave starts inundating ashore
+    for t in domain.evolve(yieldstep=10, finaltime=10000,
+                           skip_initial_step=True):
+        print domain.timestepping_statistics()
+        print domain.boundary_statistics(tags='ocean_east')
 \end{verbatim}
 
@@ -1216,17 +1199,4 @@
 This functionality is especially convenient as it allows the detailed
 parts of the simulation to be viewed at higher time resolution.
-
-\begin{verbatim}
-# save every two mins leading up to wave approaching land
-for t in domain.evolve(yieldstep=120, finaltime=5000):
-    domain.write_time()
-    domain.write_boundary_statistics(tags='ocean_east')
-
-# save every 30 secs as wave starts inundating ashore
-for t in domain.evolve(yieldstep=10, finaltime=10000, skip_initial_step=True):
-    domain.write_time()
-    domain.write_boundary_statistics(tags='ocean_east')
-\end{verbatim}
-
 
 \section{Exploring the Model Output}
@@ -3730,5 +3700,5 @@
 
 
-/pagebreak
+\pagebreak
 \section{utilities/polygons}
 
Index: /anuga_core/documentation/user_manual/demos/cairns/runcairns.py
===================================================================
--- /anuga_core/documentation/user_manual/demos/cairns/runcairns.py	(revision 7085)
+++ /anuga_core/documentation/user_manual/demos/cairns/runcairns.py	(revision 7086)
@@ -146,5 +146,5 @@
     # Continue propagating wave
     for t in domain.evolve(yieldstep=10, finaltime=5000, 
-                           skip_initial_step = True):
+                           skip_initial_step=True):
         print domain.timestepping_statistics()
         print domain.boundary_statistics(tags='ocean_east')    
Index: /anuga_core/documentation/user_manual/demos/channel1.py
===================================================================
--- /anuga_core/documentation/user_manual/demos/channel1.py	(revision 7085)
+++ /anuga_core/documentation/user_manual/demos/channel1.py	(revision 7086)
@@ -24,5 +24,5 @@
 # Setup initial conditions
 #------------------------------------------------------------------------------
-def topography(x,y):
+def topography(x, y):
     return -x/10                             # linear bed slope
 
@@ -30,5 +30,5 @@
 domain.set_quantity('friction', 0.01)        # Constant friction 
 domain.set_quantity('stage',                 # Dry bed
-                    expression='elevation + 0.0')  
+                    expression='elevation')  
 
 #------------------------------------------------------------------------------
Index: /anuga_core/documentation/user_manual/demos/runup.py
===================================================================
--- /anuga_core/documentation/user_manual/demos/runup.py	(revision 7085)
+++ /anuga_core/documentation/user_manual/demos/runup.py	(revision 7086)
@@ -29,5 +29,5 @@
 # Setup initial conditions
 #------------------------------------------------------------------------------
-def topography(x,y):
+def topography(x, y):
     return -x/2                              # linear bed slope
     #return x*(-(2.0-x)*.5)                  # curved bed slope
@@ -35,5 +35,5 @@
 domain.set_quantity('elevation', topography) # Use function for elevation
 domain.set_quantity('friction', 0.1)         # Constant friction 
-domain.set_quantity('stage', -.4)            # Constant negative initial stage
+domain.set_quantity('stage', -0.4)           # Constant negative initial stage
 
 #------------------------------------------------------------------------------
@@ -44,5 +44,5 @@
 Bd = Dirichlet_boundary([-0.2,0.,0.]) # Constant boundary values
 Bw = Time_boundary(domain=domain,     # Time dependent boundary  
-                   f=lambda t: [(.1*sin(t*2*pi)-0.3) * exp(-t), 0.0, 0.0])
+                   f=lambda t: [(0.1*sin(t*2*pi)-0.3)*exp(-t), 0.0, 0.0])
 
 # Associate boundary tags with boundary objects
@@ -52,4 +52,4 @@
 # Evolve system through time
 #------------------------------------------------------------------------------
-for t in domain.evolve(yieldstep = 0.1, finaltime = 10.0):
+for t in domain.evolve(yieldstep=0.1, finaltime=10.0):
     print domain.timestepping_statistics()
