Index: /trunk/anuga_core/source/anuga/error_api.py
===================================================================
--- /trunk/anuga_core/source/anuga/error_api.py	(revision 8260)
+++ /trunk/anuga_core/source/anuga/error_api.py	(revision 8260)
@@ -0,0 +1,130 @@
+import matplotlib
+import matplotlib.pyplot as plt
+import matplotlib.tri as tri
+from pylab import *
+
+# The abstract Python-MPI interface
+from anuga_parallel.parallel_abstraction import size, rank, get_processor_name
+from anuga_parallel.parallel_abstraction import finalize, send, receive
+from anuga_parallel.parallel_abstraction import pypar_available, barrier
+
+import random
+
+import numpy as num
+
+'''
+Adds random noise to stage quantities (centroid) in the domain.
+domain - domain where noise is to be added
+samples - number of centroids to add noise to
+dthr - minimum depth threshold, stages corresponding to depths less than this value are not modified
+mag - (relative) magnitude of the noise, uniformly distributed in [-mag, mag]
+'''
+def addRandomNoiseToStage(domain, samples = 50, dthr = 0.1, mag = 1E-15):
+
+    # Calculate depth
+    depth = domain.get_quantity('stage').centroid_values[domain.tri_full_flag == 1] \
+        - domain.get_quantity('elevation').centroid_values[domain.tri_full_flag == 1]
+
+    stage = domain.get_quantity('stage')
+
+    # Sample centroid values with depth greater than dthr
+    n_tri = int(num.sum(domain.tri_full_flag))
+    tri_index = num.array(range(0,n_tri))
+    submerged = tri_index[depth > dthr]
+    submerged = random.sample(submerged, min(samples, len(submerged)))
+
+    # Add random noise to sampled centroid values
+    if len(submerged) > 0:
+        for i in range(len(submerged)):
+            new_val = (1.0 + random.uniform(-1.0, 1.0)*mag)*stage.centroid_values[submerged[i]]
+            domain.get_quantity('stage').centroid_values.put(submerged[i], new_val)
+
+'''
+Plot errors in centroid quantities, triangles with errors are denoted in blue, while
+triangles without errors are denoted in green.
+
+domain - domain for which the error plots will be displayed
+control_data - control data to compare quantity values against (must have serial indexing)
+rthr - relative error threshold
+athr - absolute error threshold
+quantity - quantity e.g stage
+filename - output filename
+'''
+
+def plotCentroidError(domain, control_data, rthr = 1E-7, athr = 1E-12, 
+                      quantity = 'stage', filename = 'centroid_error.png'):
+
+    n_triangles = num.sum(domain.tri_full_flag)
+    
+    if size() > 1:
+        # If parallel, translate control data to parallel indexing
+        local_control_data = num.zeros(n_triangles)
+        inv_tri_map = domain.get_inv_tri_map()
+        
+        for i in range(n_triangles):
+            local_control_data[i] = control_data[inv_tri_map[(rank(), i)]]
+    else:
+        local_control_data = control_data
+
+    # Evaluate absolute and relative difference between control and actual values
+    stage = domain.get_quantity(quantity)
+    actual_data = stage.centroid_values[:n_triangles]
+    adiff = num.fabs((actual_data - local_control_data))
+    rdiff = adiff/num.fabs(local_control_data)
+
+    # Compute masks for error (err_mask) and non-error (acc_mask) vertex indices based on thresholds
+    vertices = domain.get_vertex_coordinates()
+    err_mask = rdiff > rthr
+    err_mask[adiff <= athr] = False    
+    err_mask = num.repeat(err_mask, 3)
+
+    acc_mask = ~err_mask
+    inv_tri_map = domain.get_inv_tri_map()
+
+    # Plot error and non-error triangle
+    if rank() == 0:
+        fx = {}
+        fy = {}
+        gx = {}
+        gy = {}
+
+        fx[0] = vertices[acc_mask,0]
+        fy[0] = vertices[acc_mask,1]
+        gx[0] = vertices[err_mask,0]
+        gy[0] = vertices[err_mask,1]
+
+        # Receive vertex indices of non-error triangles (fx, fy) and error triangles (gx, gy)
+        for i in range(1,size()):
+            fx[i] = receive(i)
+            fy[i] = receive(i)
+            gx[i] = receive(i)
+            gy[i] = receive(i)
+
+        # Plot non-error triangles in green
+        for i in range(0,size()):
+            n = int(len(fx[i])/3)
+            triang = num.array(range(0,3*n))
+            triang.shape = (n, 3)
+
+            if len(fx[i]) > 0:
+                plt.triplot(fx[i], fy[i], triang, 'g-')
+
+        # Plot error triangles in blue
+        for i in range(0,size()):
+            n = int(len(gx[i])/3)
+                            
+            triang = num.array(range(0,3*n))
+            triang.shape = (n, 3)
+
+            if len(gx[i]) > 0: 
+                plt.triplot(gx[i], gy[i], triang, 'b--')
+                
+        # Save plot
+        plt.savefig(filename)
+        
+    else:
+        # Send error and non-error vertex indices to Proc 0
+        send(vertices[acc_mask,0], 0)
+        send(vertices[acc_mask,1], 0)
+        send(vertices[err_mask,0], 0)
+        send(vertices[err_mask,1], 0)      
Index: /trunk/anuga_core/source/anuga_parallel/distribute_mesh.py
===================================================================
--- /trunk/anuga_core/source/anuga_parallel/distribute_mesh.py	(revision 8259)
+++ /trunk/anuga_core/source/anuga_parallel/distribute_mesh.py	(revision 8260)
@@ -101,4 +101,12 @@
 
 def pmesh_divide_metis(domain, n_procs):
+    # Wrapper for old pmesh_divide_metis which does not return tri_index or r_tri_index
+    nodes, ttriangles, boundary, triangles_per_proc, quantities, tri_index, r_tri_index = pmesh_divide_metis_helper(domain, n_procs)
+    return nodes, ttriangles, boundary, triangles_per_proc, quantities
+
+def pmesh_divide_metis_with_map(domain, n_procs):
+    return pmesh_divide_metis_helper(domain, n_procs)
+
+def pmesh_divide_metis_helper(domain, n_procs):
     
     # Initialise the lists
@@ -110,4 +118,8 @@
     
     tri_list = []
+
+    # Serial to Parallel and Parallel to Serial Triangle index maps
+    tri_index = {}
+    r_tri_index = {} # reverse tri index, parallel to serial triangle index mapping
     
     # List indexed by processor of cumulative total of triangles allocated.
@@ -148,5 +160,4 @@
         # (local to the processor)
         
-        tri_index = {}
         triangles = []        
         for i in range(n_tri):
@@ -154,4 +165,5 @@
             tri_list[epart[i]].append(domain.triangles[i])
             tri_index[i] = ([epart[i], len(tri_list[epart[i]]) - 1])
+            r_tri_index[epart[i], len(tri_list[epart[i]]) - 1] = i
         
         # Order the triangle list so that all of the triangles belonging
@@ -202,6 +214,8 @@
     for i in range(len(triangles)):
         ttriangles[i] = triangles[i]
-    
-    return nodes, ttriangles, boundary, triangles_per_proc, quantities
+
+    #return nodes, ttriangles, boundary, triangles_per_proc, quantities
+    
+    return nodes, ttriangles, boundary, triangles_per_proc, quantities, tri_index, r_tri_index
 
 
Index: /trunk/anuga_core/source/anuga_parallel/parallel_api.py
===================================================================
--- /trunk/anuga_core/source/anuga_parallel/parallel_api.py	(revision 8259)
+++ /trunk/anuga_core/source/anuga_parallel/parallel_api.py	(revision 8260)
@@ -20,5 +20,5 @@
     # Mesh partitioning using Metis
     from anuga_parallel.distribute_mesh import build_submesh
-    from anuga_parallel.distribute_mesh import pmesh_divide_metis
+    from anuga_parallel.distribute_mesh import pmesh_divide_metis_with_map
 
     from anuga_parallel.parallel_shallow_water import Parallel_domain
@@ -97,7 +97,12 @@
         points, vertices, boundary, quantities,\
                 ghost_recv_dict, full_send_dict,\
-                number_of_full_nodes, number_of_full_triangles =\
+                number_of_full_nodes, number_of_full_triangles,\
+                s2p_map, p2s_map =\
                 distribute_mesh(domain, verbose=verbose)
 
+        # Send serial to parallel (s2p) and parallel to serial (p2s) triangle mapping to proc 1 .. numprocs
+        for p in range(1, numprocs):
+            send(s2p_map, p)
+            send(p2s_map, p)
 
         if verbose: print 'Communication done'
@@ -112,4 +117,8 @@
                 rec_submesh(0, verbose)
 
+        # Recieve serial to parallel (s2p) and parallel to serial (p2s) triangle mapping
+        s2p_map = receive(0)
+        p2s_map = receive(0)
+
 
     #------------------------------------------------------------------------
@@ -125,5 +134,7 @@
                              number_of_full_nodes=number_of_full_nodes,
                              number_of_full_triangles=number_of_full_triangles,
-                             geo_reference=georef) ## jj added this
+                             geo_reference=georef,
+                             tri_map = s2p_map,
+                             inv_tri_map = p2s_map) ## jj added this
 
     #------------------------------------------------------------------------
@@ -165,8 +176,9 @@
     # Subdivide the mesh
     if verbose: print 'Subdivide mesh'
-    nodes, triangles, boundary, triangles_per_proc, quantities = \
-           pmesh_divide_metis(domain, numprocs)
-
-
+    nodes, triangles, boundary, triangles_per_proc, quantities, s2p_map, p2s_map = \
+           pmesh_divide_metis_with_map(domain, numprocs)
+
+    #PETE: s2p_map (maps serial domain triangles to parallel domain triangles)
+    #p2_map (maps parallel domain triangles to domain triangles)
 
 
@@ -227,7 +239,7 @@
     return points, vertices, boundary, quantities,\
            ghost_recv_dict, full_send_dict,\
-           number_of_full_nodes, number_of_full_triangles
-    
-
-
-
+           number_of_full_nodes, number_of_full_triangles, s2p_map, p2s_map
+    
+
+
+
Index: /trunk/anuga_core/source/anuga_parallel/parallel_shallow_water.py
===================================================================
--- /trunk/anuga_core/source/anuga_parallel/parallel_shallow_water.py	(revision 8259)
+++ /trunk/anuga_core/source/anuga_parallel/parallel_shallow_water.py	(revision 8260)
@@ -20,4 +20,10 @@
 
 
+#Import matplotlib
+
+import matplotlib
+matplotlib.use('Agg')
+import matplotlib.pyplot as plt
+import matplotlib.tri as tri
 
 class Parallel_domain(Domain):
@@ -29,5 +35,7 @@
                  number_of_full_nodes=None,
                  number_of_full_triangles=None,
-                 geo_reference=None): #jj added this
+                 geo_reference=None,
+                 tri_map=None,
+                 inv_tri_map=None): #jj added this
 
         Domain.__init__(self,
@@ -57,4 +65,7 @@
 
         setup_buffers(self)
+
+        self.tri_map = tri_map
+        self.inv_tri_map = inv_tri_map
 
 
@@ -122,2 +133,65 @@
         return N[:self.number_of_full_nodes_tmp,:]
 
+    def get_tri_map(self):
+        return self.tri_map
+
+    def get_inv_tri_map(self):
+        return self.inv_tri_map
+
+    '''
+    Outputs domain triangulation, full triangles are shown in green while ghost triangles are shown in blue.
+    The default filename is "domain.png"
+    '''
+    def dump_triangulation(self, filename="domain.png"):
+        # Get vertex coordinates, partition full and ghost triangles based on self.tri_full_flag
+        vertices = self.get_vertex_coordinates()
+        full_mask = num.repeat(self.tri_full_flag == 1, 3)
+        ghost_mask = num.repeat(self.tri_full_flag == 0, 3)
+        
+        myid = pypar.rank()
+        numprocs = pypar.size()
+
+        if myid == 0:
+            fx = {}
+            fy = {}
+            gx = {}
+            gy = {}
+
+            # Proc 0 gathers full and ghost nodes from self and other processors
+            fx[0] = vertices[full_mask,0]
+            fy[0] = vertices[full_mask,1]
+            gx[0] = vertices[ghost_mask,0]
+            gy[0] = vertices[ghost_mask,1]
+            
+            for i in range(1,numprocs):
+                fx[i] = pypar.receive(i)
+                fy[i] = pypar.receive(i)
+                gx[i] = pypar.receive(i)
+                gy[i] = pypar.receive(i)
+
+            # Plot full triangles
+            for i in range(0, numprocs):
+                n = int(len(fx[i])/3)
+                            
+                triang = num.array(range(0,3*n))
+                triang.shape = (n, 3)
+                plt.triplot(fx[i], fy[i], triang, 'g-')
+
+            # Plot ghost triangles
+            for i in range(0, numprocs):
+                n = int(len(gx[i])/3)
+                            
+                triang = num.array(range(0,3*n))
+                triang.shape = (n, 3)
+                plt.triplot(gx[i], gy[i], triang, 'b--')
+
+            # Save triangulation to location pointed by filename
+            plt.savefig(filename)
+
+        else:
+            # Proc 1..numprocs send full and ghost triangles to Proc 0
+            pypar.send(vertices[full_mask,0], 0)
+            pypar.send(vertices[full_mask,1], 0)
+            pypar.send(vertices[ghost_mask,0], 0)
+            pypar.send(vertices[ghost_mask,1], 0)
+
Index: /trunk/anuga_core/source/anuga_parallel/test_parallel_frac_op.py
===================================================================
--- /trunk/anuga_core/source/anuga_parallel/test_parallel_frac_op.py	(revision 8259)
+++ /trunk/anuga_core/source/anuga_parallel/test_parallel_frac_op.py	(revision 8260)
@@ -105,5 +105,7 @@
 ##-----------------------------------------------------------------------
 
-    if parallel: domain = distribute(domain)
+    if parallel:
+        domain = distribute(domain)
+        domain.dump_triangulation("frac_op_domain.png")
     
 
@@ -147,4 +149,8 @@
 
     ################ Define Fractional Operators ##########################
+
+    inlet0 = None
+    inlet1 = None
+    boyd_box0 = None
     
     inlet0 = Inlet_operator(domain, line0, Q0, debug = False)
@@ -152,4 +158,5 @@
     
     # Enquiry point [ 19.    2.5] is contained in two domains in 4 proc case
+    
     boyd_box0 = Boyd_box_operator(domain,
                                   end_points=[[9.0, 2.5],[19.0, 2.5]],
@@ -161,5 +168,5 @@
                                   manning=0.013,
                                   verbose=False, debug = False)
-
+        
     if inlet0 is not None: inlet0.print_statistics()
     if inlet1 is not None: inlet1.print_statistics()
@@ -202,5 +209,5 @@
     ##-----------------------------------------------------------------------
 
-    for t in domain.evolve(yieldstep = 1.0, finaltime = 38):
+    for t in domain.evolve(yieldstep = 0.1, finaltime = 38):
         domain.write_time()
 
