Index: /inundation/parallel/documentation/RunParallelAdvection.py
===================================================================
--- /inundation/parallel/documentation/RunParallelAdvection.py	(revision 2767)
+++ /inundation/parallel/documentation/RunParallelAdvection.py	(revision 2767)
@@ -0,0 +1,115 @@
+#!/usr/bin/env python
+###
+#########################################################
+#
+#  Main file for parallel mesh testing. Runs an advection
+# flow simulation using a rectangular mesh
+#
+#
+#  Authors: Steve Roberts June 2005
+#  Modified by Linda Stals April 2006
+#
+#
+#########################################################
+
+
+import sys
+from os import sep
+sys.path.append('..'+sep+'pyvolution')
+
+# Parallel communication routines
+
+import pypar
+
+#  Mesh partition routines
+
+from parallel_meshes import parallel_rectangle
+
+# Parallel Domain
+ 
+from parallel_advection import Parallel_Domain
+from parallel_advection import Transmissive_boundary
+
+############################
+# Set the initial conditions
+############################
+class Set_Stage:
+    """Set an initial condition with constant water height, for x<x0
+    """
+
+    def __init__(self, x0=0.25, x1=0.5, h=1.0):
+        self.x0 = x0
+        self.x1 = x1
+        self.h  = h
+
+    def __call__(self, x, y):
+        return self.h*((x>self.x0)&(x<self.x1))
+
+###############################
+# Read in processor information
+###############################
+
+numprocs = pypar.size()
+myid = pypar.rank()
+processor_name = pypar.Get_processor_name()
+
+N = 20
+M = 20
+
+#######################
+# Partition the domain
+#######################
+
+# Build a unit mesh, subdivide it over numproces processors with each
+# submesh containing M*N nodes
+
+points, vertices, boundary, full_send_dict, ghost_recv_dict =  \
+    parallel_rectangle(N, M, len1_g=1.0)
+
+###########################################
+# Start the computations on each subpartion
+###########################################
+
+# Create advection domain with direction (1,-1)
+# Initial condition is zero by default
+
+domain = Parallel_Domain(points, vertices, boundary,
+                         full_send_dict, ghost_recv_dict, velocity=[1.0, 0.0])
+
+# Turn on the visualisation
+
+rect = [0.0, 0.0, 1.0, 1.0]
+domain.initialise_visualiser(rect=rect)
+
+# Boundaries
+
+T = Transmissive_boundary(domain)
+domain.default_order = 2
+domain.set_boundary( {'left': T, 'right': T, 'bottom': T, 'top': T, 'ghost': None} )
+
+# Ensure that the domain definitions make sense
+
+domain.check_integrity()
+
+# Set the inititial conditions
+
+domain.set_quantity('stage', Set_Stage(0.2,0.4,1.0))
+
+# Let processor 0 output some timing information
+
+if myid == 0:
+    import time
+    t0 = time.time()
+
+# Start the evolve computions
+
+for t in domain.evolve(yieldstep = 0.1, finaltime = 3.0):
+    if myid == 0:
+        domain.write_time()
+
+# Output some computation statistics
+
+if myid == 0:
+    print 'That took %.2f seconds' %(time.time()-t0)
+    print 'Communication time %.2f seconds'%domain.communication_time
+    print 'Reduction Communication time %.2f seconds'%domain.communication_reduce_time
Index: /inundation/parallel/documentation/RunParallelMerimbulaMetis.py
===================================================================
--- /inundation/parallel/documentation/RunParallelMerimbulaMetis.py	(revision 2767)
+++ /inundation/parallel/documentation/RunParallelMerimbulaMetis.py	(revision 2767)
@@ -0,0 +1,178 @@
+#!/usr/bin/env python
+#########################################################
+#
+#  Main file for parallel mesh testing. Runs an advection
+# flow simulation using a rectangular mesh
+#
+#  This is a modification of the run_parallel_advection.py
+# file
+#
+#
+# *) The (new) files that have been added to manage the
+# grid partitioning are
+#    +) pmesh_divide_metis.py: subdivide a pmesh
+#    +) build_submesh.py: build the submeshes on the host
+# processor.
+#    +) build_local.py: build the GA mesh datastructure
+# on each processor.
+#    +) build_commun.py: handle the communication between
+# the host and processors
+#
+#  Authors: Linda Stals, Steve Roberts and Matthew Hardy,
+# June 2005
+#
+#
+#
+#########################################################
+import sys
+import pypar    # The Python-MPI interface
+import time
+
+
+from os import sep
+sys.path.append('..'+sep+'pyvolution')
+
+# Numeric arrays
+
+from Numeric import array, zeros, Float
+
+# pmesh
+
+from pmesh2domain import pmesh_to_domain_instance
+from advection import Domain as Advection_Domain
+from parallel_advection import Parallel_Domain
+
+from generic_boundary_conditions import Transmissive_boundary
+
+# Mesh partition routines
+
+from pmesh_divide  import pmesh_divide_metis
+from build_submesh import build_submesh, extract_hostmesh
+from build_local   import build_local_mesh
+from build_commun  import send_submesh, rec_submesh
+
+
+############################
+# Set the initial conditions
+############################
+
+class Set_Stage:
+    """Set an initial condition with constant water height, for x<x0
+    """
+
+    def __init__(self, x0=0.25, x1=0.5, h=1.0):
+        self.x0 = x0
+        self.x1 = x1
+        self.h  = h
+
+    def __call__(self, x, y):
+        return self.h*((x>self.x0)&(x<self.x1))
+
+rect = zeros( 4, Float) # Buffer for results
+
+###############################
+# Read in processor information
+###############################
+
+numprocs = pypar.size()
+myid = pypar.rank()
+processor_name = pypar.Get_processor_name()
+
+#######################
+# Partition the domain
+#######################
+
+if myid == 0:
+
+    # Read in the test files
+
+    filename = 'merimbula_10785.tsh'
+
+    domain_full = pmesh_to_domain_instance(filename, Advection_Domain)
+    domain_full.set_quantity('stage', Set_Stage(756000.0,756500.0,4.0))
+
+    # Define the domain boundaries for visualisation
+   
+    rect = array(domain_full.xy_extent, Float)
+
+    # Subdivide the mes
+
+    nodes, triangles, boundary, triangles_per_proc, quantities  =\
+            pmesh_divide_metis(domain_full, numprocs)
+
+    # Build the mesh that should be assigned to each processor,
+    # this includes ghost nodes and the communicaiton pattern
+    
+    submesh = build_submesh(nodes, triangles, boundary, quantities, \
+                            triangles_per_proc)
+
+    # Send the mesh partition to the appropriate processor
+
+    for p in range(1, numprocs):
+      send_submesh(submesh, triangles_per_proc, p)
+
+    # Build the local mesh for processor 0
+
+    hostmesh = extract_hostmesh(submesh)
+    points, vertices, boundary, quantities, ghost_recv_dict, full_send_dict = \
+             build_local_mesh(hostmesh, 0, triangles_per_proc[0], numprocs)
+
+else:
+    # Read in the mesh partition that belongs to this
+    # processor (note that the information is in the
+    # correct form for the GA data structure
+
+    points, vertices, boundary, quantities, ghost_recv_dict, full_send_dict = \
+             rec_submesh(0)
+
+###########################################
+# Start the computations on each subpartion
+###########################################
+
+# The visualiser needs to know the size of the whole domain
+
+pypar.broadcast(rect,0)
+
+# Build the domain for this processor
+
+domain = Parallel_Domain(points, vertices, boundary,
+                         full_send_dict  = full_send_dict,
+                         ghost_recv_dict = ghost_recv_dict,
+                         velocity = [0.1,0.0])
+
+# Visualise the domain
+
+try:
+    domain.initialise_visualiser(rect=rect)
+    #domain.visualiser.coloring['stage'] = True
+    #domain.visualiser.scale_z['stage'] = 0.2
+except:
+    print 'No visualiser'
+    
+# Define the boundaries, including the ghost boundary
+
+T = Transmissive_boundary(domain)
+#R = Reflective_boundary(domain)
+domain.set_boundary( {'outflow': T, 'inflow': T, 'inner':T, 'exterior': T, 'open':T, 'ghost':None} )
+
+# Set the initial quantities
+
+
+domain.set_quantity('stage', quantities['stage'])
+
+# Set the number of time steps, as well as the start and end time
+
+t0 = time.time()
+yieldstep = 1000
+finaltime = 50000
+
+# Start the evolve calculations
+
+for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
+    if myid == 0:
+        domain.write_time()
+
+# Print some timing statistics
+
+if myid == 0:
+    print 'That took %.2f seconds' %(time.time()-t0)
Index: /inundation/parallel/documentation/RunParallelSwMerimbulaMetis.py
===================================================================
--- /inundation/parallel/documentation/RunParallelSwMerimbulaMetis.py	(revision 2767)
+++ /inundation/parallel/documentation/RunParallelSwMerimbulaMetis.py	(revision 2767)
@@ -0,0 +1,192 @@
+#!/usr/bin/env python
+###
+#########################################################
+#
+#  Main file for parallel mesh testing. Runs a shallow
+# water simulation using the merimbula mesh
+#
+#
+#
+# *) The (new) files that have been added to manage the
+# grid partitioning are
+#    +) pmesh_divide_metis.py: subdivide a pmesh
+#    +) build_submesh.py: build the submeshes on the host
+# processor.
+#    +) build_local.py: build the GA mesh datastructure
+# on each processor.
+#    +) build_commun.py: handle the communication between
+# the host and processors
+#
+#  Authors: Linda Stals, Steve Roberts and Matthew Hardy,
+# June 2005
+#
+#
+#
+#########################################################
+import sys
+import pypar    # The Python-MPI interface
+import time
+
+from os import sep
+sys.path.append('..'+sep+'pyvolution')
+
+# Numeric arrays
+
+from Numeric import array, zeros, Float
+
+# pmesh
+
+from shallow_water import Domain
+from parallel_shallow_water import Parallel_Domain
+from pmesh2domain import pmesh_to_domain_instance
+
+# Mesh partition routines
+
+from pmesh_divide import pmesh_divide_metis
+from build_submesh import build_submesh, extract_hostmesh
+from build_local   import build_local_mesh
+from build_commun  import send_submesh, rec_submesh
+
+###############################
+# Read in processor information
+###############################
+
+numprocs = pypar.size()
+myid = pypar.rank()
+processor_name = pypar.Get_processor_name()
+
+############################
+# Set the initial conditions
+############################
+
+rect = zeros( 4, Float) # Buffer for results
+
+class Set_Stage:
+    """Set an initial condition with constant water height, for x<x0
+    """
+
+    def __init__(self, x0=0.25, x1=0.5, h=1.0):
+        self.x0 = x0
+        self.x1 = x1
+        self.h  = h
+
+    def __call__(self, x, y):
+        return self.h*((x>self.x0)&(x<self.x1))
+
+#######################
+# Partition the domain
+#######################
+
+if myid == 0:
+
+    # Read in the test files
+
+    filename = 'merimbula_10785_1.tsh'
+
+    # Build the whole domain
+    
+    domain_full = pmesh_to_domain_instance(filename, Domain)
+
+    # Define the domain boundaries for visualisation
+
+    rect = array(domain_full.xy_extent, Float)
+
+    # Initialise the wave
+
+    domain_full.set_quantity('stage', Set_Stage(756000.0,756500.0,2.0))
+
+    # Subdivide the mesh
+    
+    nodes, triangles, boundary, triangles_per_proc, quantities = \
+         pmesh_divide_metis(domain_full, numprocs)
+
+    # Build the mesh that should be assigned to each processor,
+    # this includes ghost nodes and the communicaiton pattern
+    
+    submesh = build_submesh(nodes, triangles, boundary,\
+                            quantities, triangles_per_proc)
+
+    # Send the mesh partition to the appropriate processor
+
+    for p in range(1, numprocs):
+      send_submesh(submesh, triangles_per_proc, p)
+
+    # Build the local mesh for processor 0
+
+    hostmesh = extract_hostmesh(submesh)
+    points, vertices, boundary, quantities, ghost_recv_dict, full_send_dict = \
+             build_local_mesh(hostmesh, 0, triangles_per_proc[0], numprocs)
+
+else:
+    
+    # Read in the mesh partition that belongs to this
+    # processor (note that the information is in the
+    # correct form for the GA data structure)
+
+    points, vertices, boundary, quantities, ghost_recv_dict, full_send_dict \
+            = rec_submesh(0)
+
+
+###########################################
+# Start the computations on each subpartion
+###########################################
+
+# The visualiser needs to know the size of the whole domain
+
+pypar.broadcast(rect,0)
+
+# Build the domain for this processor
+
+domain = Parallel_Domain(points, vertices, boundary,
+                         full_send_dict  = full_send_dict,
+                         ghost_recv_dict = ghost_recv_dict)
+
+# Visualise the domain
+
+try:
+    domain.initialise_visualiser(rect=rect)
+    domain.visualiser.scale_z['stage'] = 0.2
+    domain.visualiser.scale_z['elevation'] = 0.05
+except:
+    print 'No visualiser'
+
+
+domain.default_order = 1
+
+# Define the boundaries, including the ghost boundary
+
+from parallel_shallow_water import Transmissive_boundary, Reflective_boundary
+
+T = Transmissive_boundary(domain)
+R = Reflective_boundary(domain)
+domain.set_boundary( {'outflow': R, 'inflow': R, 'inner':R, 'exterior': R, \
+                      'open':R, 'ghost':None} )
+
+
+# Set the initial quantities
+
+domain.set_quantity('stage', quantities['stage'])
+domain.set_quantity('elevation', quantities['elevation'])
+
+domain.store = False
+
+# Set the number of time steps, as well as the start and end time
+
+t0 = time.time()
+yieldstep = 1
+finaltime = 90
+
+
+# Start the evolve calculations
+
+for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
+    if myid == 0:
+        domain.write_time()
+
+# Print some timing statistics
+
+if myid == 0:
+    print 'That took %.2f seconds' %(time.time()-t0)
+    print 'Communication time %.2f seconds'%domain.communication_time
+    print 'Reduction Communication time %.2f seconds'%domain.communication_reduce_time
+    print 'Broadcast time %.2f seconds'%domain.communication_broadcast_time
Index: /inundation/parallel/documentation/parallel.tex
===================================================================
--- /inundation/parallel/documentation/parallel.tex	(revision 2766)
+++ /inundation/parallel/documentation/parallel.tex	(revision 2767)
@@ -1,17 +1,17 @@
 \chapter{Running the Code in Parallel}
 
-This chapter focuses on how to run the code in parallel. The first section gives an overview of the main steps used to divide the domain over the processors. The second sections gives some example code and the final section talks about how to run the code on some example architectures.
-
-
-\section {Partitioning the Domain}
+This chapter looks at how to run the code in parallel. The first section gives an overview of the main steps required to divide the mesh over the processors. The second sections gives some example code and the final section talks about how to run the code on a few specific architectures.
+
+
+\section {Partitioning the Mesh}
 There are four main steps required to run the code in parallel. They are;
 \begin{enumerate}
-\item subdivide the domain into a set of non-overlapping subdomains
+\item partition the mesh into a set of non-overlapping submeshes
 (\code{pmesh_divide_metis} from {\tt pmesh_divde.py}),
-\item build a \lq ghost\rq\ or communication layer of boundary triangles
-around each subdomain and define the communication pattern (\code{build_submesh} from {\tt build_submesh.py}),
-\item distribute the subdomains over the processors (\code{send_submesh}
+\item build a \lq ghost\rq\ or communication layer of triangles
+around each submesh and define the communication pattern (\code{build_submesh} from {\tt build_submesh.py}),
+\item distribute the submeshes over the processors (\code{send_submesh}
 and \code{rec_submesh} from {\tt build_commun.py}),
-\item and update the numbering scheme for the local domain assigned to a
+\item and update the numbering scheme for each submesh assigned to a
 processor (\code{build_local_mesh} from {\tt build_local.py}).
 \end{enumerate}
@@ -20,19 +20,19 @@
 \begin{figure}[hbtp]
   \centerline{ \includegraphics[scale = 0.75]{domain.eps}}
-  \caption{The main steps in dividing the domain over the processors.}
+  \caption{The main steps used to divide the mesh over the processors.}
   \label{fig:subpart}
 \end{figure}
 
-\subsection {Subdividing the Global Domain}
-
-The first step in parallelising the code is to subdivide the domain
-into equally sized partitions. On a rectangular domain this may be
+\subsection {Subdividing the Global Mesh}
+
+The first step in parallelising the code is to subdivide the mesh
+into, roughly, equally sized partitions. On a rectangular domain this may be
 done by a simple co-ordinate based dissection, but on a complicated
-domain such as the Merimbula grid shown in Figure \ref{fig:subpart}
+domain such as the Merimbula grid shown in Figure \ref{fig:mergrid}
 a more sophisticated approach must be used.  We use pymetis, a
 python wrapper around the Metis
 (\url{http://glaros.dtc.umn.edu/gkhome/metis/metis/overview})
 partitioning library. The \code{pmesh_divide_metis} function defined
-in {\tt pmesh_divide.py} uses Metis to divide the domain for
+in {\tt pmesh_divide.py} uses Metis to divide the mesh for
 parallel computation. Metis was chosen as the partitioner based on
 the results in the paper \cite{gk:metis}.
@@ -44,30 +44,49 @@
 \end{figure}
 
-Figure \ref{fig:mermesh4} shows the Merimbula grid partitioned over four processor. Table \ref{tbl:mermesh4} gives the node distribution over the four processors while Table \ref{tbl:mermesh8} shows the distribution over eight processors. These results imply that Pymetis gives a reasonably well balanced partition of the domain.
-
-\begin{figure}[hbtp]
-  \centerline{ \includegraphics[scale = 0.75]{mermesh4.eps}}
+Figure \ref{fig:mergrid4} shows the Merimbula grid partitioned over four processor. Table \ref{tbl:mer4} gives the node distribution over the four processors while Table \ref{tbl:mer8} shows the distribution over eight processors. These results imply that Pymetis gives a reasonably well balanced partition of the mesh.
+
+\begin{figure}[hbtp]
+  \centerline{ \includegraphics[scale = 0.75]{mermesh4c.eps}
+  \includegraphics[scale = 0.75]{mermesh4a.eps}}
+
+
+  \centerline{ \includegraphics[scale = 0.75]{mermesh4d.eps}
+  \includegraphics[scale = 0.75]{mermesh4b.eps}} 
   \caption{The Merimbula grid partitioned over 4 processors using Metis.}
  \label{fig:mergrid4}
 \end{figure}
 
-\begin{table} \label{tbl:mermesh8}
+\begin{table} 
+\caption{Running an 4-way test of
+{\tt run_parallel_sw_merimbula_metis.py}}\label{tbl:mer4}
+\begin{center}
+\begin{tabular}{c|c c c c c c c c}
+CPU & 0 & 1 & 2 & 3\\
+Elements & 2757 & 2713 & 2761 & 2554\\
+\% & 25.56\% & 25.16\% & 25.60\% & 23.68\%\
+\end{tabular}
+\end{center}
+\end{table}
+
+\begin{table}
 \caption{Running an 8-way test of
-{\tt run_parallel_sw_merimbula_metis.py}}
+{\tt run_parallel_sw_merimbula_metis.py}} \label{tbl:mer8}
+\begin{center}
 \begin{tabular}{c|c c c c c c c c}
 CPU & 0 & 1 & 2 & 3 & 4 & 5 & 6 & 7\\
-Elements & 1292 & 1647 & 1505 & 1623 & 1381 & 1514 & 1605 & 1388\\
-\% & 10.81\% & 13.78\% & 12.59\% & 13.58\% & 11.55\% & 12.66\% & 13.43\% & 11.61\%\\
+Elements & 1229 & 1293 & 1352 & 1341 & 1349 & 1401 & 1413 & 1407\\
+\% & 11.40\% & 11.99\% & 12.54\% & 12.43\% & 12.51\% & 12.99\% & 13.10\% & 13.05\%\\
 \end{tabular}
+\end{center}
 \end{table}
 
-The number of subdomains found by Pymetis is equal to the number of processors; Subdomain $p$ will be assigned to Processor $p$.
+The number of submeshes found by Pymetis is equal to the number of processors; Submesh $p$ will be assigned to Processor $p$.
 
 \subsection {Building the Ghost Layer}
 The function {\tt build_submesh.py} is the work-horse and is responsible for
-setting up the communication pattern as well deciding the local numbering scheme for the subdomains.
-
-Consider the example subpartitioning given in Figure \ref{fig:subdomain}. During the evolve calculations Triangle 3 in Subdomain 0 will need to access its neighbour Triangle 4 stored in Subdomain 1. The standard approach to this problem is to add an extra layer of triangles, which we call ghost triangles. The ghost triangles
-are read-only, they should not be updated during the calculations, they are only there to hold any extra information that a processor may need to complete its calculations. The ghost triangle values are updated through communication calls. Figure \ref{fig:subdomaing} shows the subdomains with the extra layer of ghost triangles.
+setting up the communication pattern as well as assigning the local numbering scheme for the submeshes.
+
+Consider the example subpartitioning given in Figure \ref{fig:subdomain}. During the \code{evolve} calculations Triangle 3 in Submesh 0 will need to access its neighbour Triangle 4 stored in Submesh 1. The standard approach to this problem is to add an extra layer of triangles, which we call ghost triangles. The ghost triangles
+are read-only, they should not be updated during the calculations, they are only there to hold any extra information that a processor may need to complete its calculations. The ghost triangle values are updated through communication calls. Figure \ref{fig:subdomaing} shows the submeshes with the extra layer of ghost triangles.
 
 \begin{figure}[hbtp]
@@ -84,31 +103,93 @@
 \end{figure}
 
-Looking at Figure \ref{fig:subdomaing} we see that after each evolve step Processor 0  will have to send the updated values for Triangle 3 and Triangle 5 to Processor 1, and similarly Processor 1 will have to send the updated values for triangles 4, 7 and 6 (recall that Subdomain $p$ will be assigned to Processor $p$). The \code{build_submesh} function builds a dictionary that defines the communication pattern.
-
-Finally, the ANUGA code assumes that the triangles (and nodes etc.) are numbered consecutively starting from 1. Consequently, if the subdomain defined in Processor 1 in Figure \ref{fig:subdomaing} is passed into the evolve calculations it would crash. The \code{build_submesh} function determines a local numbering scheme for each subdomain, but it does not actually update the numbering, that is left to \code{build_local}.
-
-\subsection {Sending the Subdomains}
-
-All of functions described so far must be run in serial on Processor 0, the next step is to start the parallel computation by spreading the subdomains over the processors. The communication is carried out by
+When partitioning the mesh we introduce new, dummy, boundary edges. For example, Triangle 3 in Submesh 1, from Figure \ref{fig:subdomaing}, originally shared an edge with Triangle 2, but after partitioning that edge becomes a boundary edge. These new boundary edges are are tagged as \code{ghost} and should, in general, be assigned a type of \code{None}. The following piece of code taken from {\tt run_parallel_advection.py} shows an example.  
+ 
+{\small \begin{verbatim}  
+T = Transmissive_boundary(domain)
+domain.default_order = 2
+domain.set_boundary( {'left': T, 'right': T, 'bottom': T, 'top': T, 'ghost': Non
+e} )
+\end{verbatim}}
+
+
+Looking at Figure \ref{fig:subdomaing} we see that after each \code{evolve} step Processor 0  will have to send the updated values for Triangle 3 and Triangle 5 to Processor 1, and similarly Processor 1 will have to send the updated values for triangles 4, 7 and 6 (recall that Submesh $p$ will be assigned to Processor $p$). The \code{build_submesh} function builds a dictionary that defines the communication pattern.
+
+Finally, the ANUGA code assumes that the triangles (and nodes etc.) are numbered consecutively starting from 1. Consequently, if Submesh 1 in Figure \ref{fig:subdomaing} was passed into the \code{evolve} calculations it would crash. The \code{build_submesh} function determines a local numbering scheme for each submesh, but it does not actually update the numbering, that is left to \code{build_local}.
+
+\subsection {Sending the Submeshes}
+
+All of functions described so far must be run in serial on Processor 0, the next step is to start the parallel computation by spreading the submeshes over the processors. The communication is carried out by
 \code{send_submesh} and \code{rec_submesh} defined in {\tt build_commun.py}.
-The \code{send_submesh} function should be called on Processor 0 and sends the Subdomain $p$ to Processor $p$, while \code{rec_submesh} should be called by Processor $p$ to receive Subdomain $p$ from Processor 0. Note that the order of communication is very important, if any changes are made to the \code{send_submesh} function the corresponding change must be made to the \code{rec_submesh} function.
-
-While it is possible to get Processor 0 to communicate it's subdomain to itself, it is an expensive unnessary communication call. The {\tt build_commun.py} file also includes a function called \code{extract_hostmesh} which simply extracts Subdomain $0$.
-
-
-\subsection {Building the Local Domain}
-After using \code{send_submesh} and \code{rec_submesh}, Processor $p$ should have its own local copy of Subdomain $p$, however as stated previously the triangle numbering may be incorrect. The \code{build_local_mesh} function from {\tt build_local.py} primarily focuses on renumbering the information stored with the subdomain; including the nodes, vertices and quantities. Figure \ref{fig:subdomainf} shows what the domain in each processor may look like.
+The \code{send_submesh} function should be called on Processor 0 and sends the Submesh $p$ to Processor $p$, while \code{rec_submesh} should be called by Processor $p$ to receive Submesh $p$ from Processor 0. Note that the order of communication is very important, if any changes are made to the \code{send_submesh} function the corresponding change must be made to the \code{rec_submesh} function.
+
+While it is possible to get Processor 0 to communicate it's submesh to itself, it is an expensive and unnecessary communication call. The {\tt build_commun.py} file also includes a function called \code{extract_hostmesh} that should be called on Processor 0 to extract Submesh 0.
+
+
+\subsection {Building the Local Mesh}
+After using \code{send_submesh} and \code{rec_submesh}, Processor $p$ should have its own local copy of Submesh $p$, however as stated previously the triangle numbering may be incorrect. The \code{build_local_mesh} function from {\tt build_local.py} primarily focuses on renumbering the information stored with the submesh; including the nodes, vertices and quantities. Figure \ref{fig:subdomainf} shows what the mesh in each processor may look like.
 
 
 \begin{figure}[hbtp]
   \centerline{ \includegraphics[scale = 0.6]{subdomainfinal.eps}}
-  \caption{An example subpartioning after the domain has been renumbered
-and defined on each processor.}
+  \caption{An example subpartioning after the submeshes have been renumbered.}
  \label{fig:subdomainf}
 \end{figure}
 
-Note that ghost triangles are not stored any differently to the other triangles, the only way to differentiate them is to look at the communication pattern. This means that the {\tt innundation} code is essentially the same whether it is being run in serial or parallel, the only difference is a communication call at the end of each evolve step and a check to ensure that the ghost triangles are not used in the time stepping calculations.
+Note that ghost triangles are not stored in the domain any differently to the other triangles, the only way to differentiate them is to look at the communication pattern. This means that the {\tt inundation} code is essentially the same whether it is being run in serial or parallel, the only difference is a communication call at the end of each \code{evolve} step and a check to ensure that the ghost triangles are not used in the time stepping calculations.
 
 \section{Some Example Code}
+
+\begin{figure}
+\begin{verbatim}
+if myid == 0:
+
+    # Read in the test files
+
+    filename = 'merimbula_10785_1.tsh'
+
+    # Build the whole domain
+    
+    domain_full = pmesh_to_domain_instance(filename, Domain)
+
+    # Define the domain boundaries for visualisation
+
+    rect = array(domain_full.xy_extent, Float)
+
+    # Initialise the wave
+
+    domain_full.set_quantity('stage', Set_Stage(756000.0,756500.0,2.0))
+
+    # Subdivide the mesh
+    
+    nodes, triangles, boundary, triangles_per_proc, quantities = \
+         pmesh_divide_metis(domain_full, numprocs)
+
+    # Build the mesh that should be assigned to each processor,
+    # this includes ghost nodes and the communicaiton pattern
+    
+    submesh = build_submesh(nodes, triangles, boundary,\
+                            quantities, triangles_per_proc)
+
+    # Send the mesh partition to the appropriate processor
+
+    for p in range(1, numprocs):
+      send_submesh(submesh, triangles_per_proc, p)
+
+    # Build the local mesh for processor 0
+
+    hostmesh = extract_hostmesh(submesh)
+    points, vertices, boundary, quantities, ghost_recv_dict, full_send_dict = \
+             build_local_mesh(hostmesh, 0, triangles_per_proc[0], numprocs)
+
+else:
+    
+    # Read in the mesh partition that belongs to this
+    # processor (note that the information is in the
+    # correct form for the GA data structure)
+
+    points, vertices, boundary, quantities, ghost_recv_dict, full_send_dict \
+            = rec_submesh(0)
+\end{verbatim}
+\end{figure}
 
 \section{Running the Code}
Index: /inundation/parallel/documentation/report.tex
===================================================================
--- /inundation/parallel/documentation/report.tex	(revision 2766)
+++ /inundation/parallel/documentation/report.tex	(revision 2767)
@@ -64,4 +64,5 @@
 \include{results}
 \include{visualisation}
+\include{code}
 
 \begin{thebibliography}{10}
