Index: uga_core/source/anuga_parallel/test_distribute_domain.py
===================================================================
--- /anuga_core/source/anuga_parallel/test_distribute_domain.py	(revision 7458)
+++ 	(revision )
@@ -1,221 +1,0 @@
-
-"""Test a run of the sequential shallow water domain against
-a run of the parallel shallow water domain.
-
-WARNING: This assumes that the command to run jobs is mpirun.
-Tested with MPICH and LAM (Ole)
-"""
-
-#------------------------------------------------------------------------------
-# Import necessary modules
-#------------------------------------------------------------------------------
-
-import unittest
-import os
-import sys
-import pypar
-
-import numpy as num
-
-
-
-
-from anuga.utilities.numerical_tools import ensure_numeric
-from anuga.utilities.util_ext        import double_precision
-from anuga.utilities.norms           import l1_norm, l2_norm, linf_norm
-
-from anuga.interface import Domain
-from anuga.interface import Reflective_boundary
-from anuga.interface import Dirichlet_boundary
-from anuga.interface import Time_boundary
-from anuga.interface import Transmissive_boundary
-
-from anuga.interface import rectangular_cross
-from anuga.interface import create_domain_from_file
-
-
-from anuga_parallel.interface import distribute, myid, numprocs
-
-
-#--------------------------------------------------------------------------
-# Setup parameters
-#--------------------------------------------------------------------------
-
-mesh_filename = "merimbula_10785_1.tsh"
-#mesh_filename = "test-100.tsh"
-yieldstep = 1
-finaltime = 20
-quantity = 'stage'
-nprocs = 4
-
-#--------------------------------------------------------------------------
-# Setup procedures
-#--------------------------------------------------------------------------
-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))
-
-#--------------------------------------------------------------------------
-# Setup test
-#--------------------------------------------------------------------------
-def evolution_test(parallel=False):
-
-
-    domain = create_domain_from_file(mesh_filename)
-    domain.set_quantity('stage', Set_Stage(756000.0, 756500.0, 2.0))
-
-    #--------------------------------------------------------------------------
-    # Create parallel domain if requested
-    #--------------------------------------------------------------------------
-
-    if parallel:
-        if myid == 0: print 'DISTRIBUTING PARALLEL DOMAIN'
-        domain = distribute(domain)
-
-    #------------------------------------------------------------------------------
-    # Setup boundary conditions
-    # This must currently happen *after* domain has been distributed
-    #------------------------------------------------------------------------------
-    domain.store = False
-    Br = Reflective_boundary(domain)      # Solid reflective wall
-
-    domain.set_boundary({'outflow' :Br, 'inflow' :Br, 'inner' :Br, 'exterior' :Br, 'open' :Br})
-
-    #------------------------------------------------------------------------------
-    # Setup diagnostic arrays
-    #------------------------------------------------------------------------------
-    l1list = []
-    l2list = []
-    linflist = []
-    l1norm = num.zeros(3, num.float)
-    l2norm = num.zeros(3, num.float)
-    linfnorm = num.zeros(3, num.float)
-    recv_norm = num.zeros(3, num.float)
-
-    #------------------------------------------------------------------------------
-    # Evolution
-    #------------------------------------------------------------------------------
-    if parallel:
-        if myid == 0: print 'PARALLEL EVOLVE'
-    else:
-        print 'SEQUENTIAL EVOLVE'
-        
-    for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
-        edges = domain.quantities[quantity].edge_values.take(num.flatnonzero(domain.tri_full_flag),axis=0)
-        l1norm[0] = l1_norm(edges[:,0])
-        l1norm[1] = l1_norm(edges[:,1])
-        l1norm[2] = l1_norm(edges[:,2])
-        l2norm[0] = l2_norm(edges[:,0])
-        l2norm[1] = l2_norm(edges[:,1])
-        l2norm[2] = l2_norm(edges[:,2])
-        linfnorm[0] = linf_norm(edges[:,0])
-        linfnorm[1] = linf_norm(edges[:,1])
-        linfnorm[2] = linf_norm(edges[:,2])
-        if parallel:
-            l2norm[0] = pow(l2norm[0], 2)
-            l2norm[1] = pow(l2norm[1], 2)
-            l2norm[2] = pow(l2norm[2], 2)
-            if myid == 0:
-                #domain.write_time()
-
-                #print edges[:,1]            
-                for p in range(1, numprocs):
-                    recv_norm = pypar.receive(p)
-                    l1norm += recv_norm
-                    recv_norm = pypar.receive(p)
-                    l2norm += recv_norm
-                    recv_norm = pypar.receive(p)
-                    linfnorm[0] = max(linfnorm[0], recv_norm[0])
-                    linfnorm[1] = max(linfnorm[1], recv_norm[1])
-                    linfnorm[2] = max(linfnorm[2], recv_norm[2])
-
-                l2norm[0] = pow(l2norm[0], 0.5)
-                l2norm[1] = pow(l2norm[1], 0.5)
-                l2norm[2] = pow(l2norm[2], 0.5)
-
-                l1list.append(l1norm)                
-                l2list.append(l2norm)
-                linflist.append(linfnorm)                
-            else:
-                pypar.send(l1norm, 0)
-                pypar.send(l2norm, 0)
-                pypar.send(linfnorm, 0)
-        else:
-            #domain.write_time()
-            l1list.append(l1norm)                
-            l2list.append(l2norm)
-            linflist.append(linfnorm)
-            
-
-    return (l1list, l2list, linflist)
-
-# Test an nprocs-way run of the shallow water equations
-# against the sequential code.
-
-class Test_distribute_domain(unittest.TestCase):
-    def test_distribute_domain(self):
-        print "Expect this test to fail if not run from the parallel directory."
-        result = os.system("mpirun -np %d python test_distribute_domain.py" % nprocs)
-        assert_(result == 0)
-
-# Because we are doing assertions outside of the TestCase class
-# the PyUnit defined assert_ function can't be used.
-def assert_(condition, msg="Assertion Failed"):
-    if condition == False:
-        #pypar.finalize()
-        raise AssertionError, msg
-
-if __name__=="__main__":
-    if numprocs == 1: 
-        runner = unittest.TextTestRunner()
-        suite = unittest.makeSuite(Test_distribute_domain, 'test')
-        runner.run(suite)
-    else:
-
-        pypar.barrier()
-        if myid == 0:
-            print 'SEQUENTIAL START'
-            l1norm_seq, l2norm_seq, linfnorm_seq = evolution_test(parallel=False)
-
-        pypar.barrier()
-        if myid ==0:
-            print 'PARALLEL START'
-        
-        l1norm_par, l2norm_par, linfnorm_par = evolution_test(parallel=True)
-        
-        if myid == 0:
-            assert_(len(l1norm_seq) == len(l1norm_par))
-            assert_(len(l2norm_seq) == len(l2norm_par))
-            assert_(len(linfnorm_seq) == len(linfnorm_par))
-            assert_(len(l1norm_seq) == len(l2norm_seq))
-            assert_(len(l2norm_seq) == len(linfnorm_seq))
-            # Anything smaller than tol we consider to be 0.
-            # This usualy comes out to be 10^-14 (DBL_DIG = 10^-15)
-            # * 10 to cater for rounding error in computation.
-            tol = pow(10, -1 * (double_precision() - 1))
-            for x in range(len(l1norm_seq)):
-                for y in range(3):
-                    # Calculate relative difference in the norms
-                    assert_(abs(l1norm_seq[x][y] - l1norm_par[x][y])/l1norm_seq[x][y] < tol)
-                    assert_(abs(l2norm_seq[x][y] - l2norm_par[x][y])/l2norm_seq[x][y] < tol)
-                    assert_(abs(linfnorm_seq[x][y] - linfnorm_par[x][y])/linfnorm_seq[x][y] < tol)
-                    if x > 0:
-                        # Verify that the quantity is being conserved across iterations.
-                        assert_(abs(l1norm_seq[x][y] - l1norm_seq[x-1][y]) < tol)
-                        assert_(abs(l2norm_seq[x][y] - l2norm_seq[x-1][y]) < tol)
-                        assert_(abs(linfnorm_seq[x][y] - linfnorm_seq[x-1][y]) < tol)
-                        assert_(abs(l1norm_par[x][y] - l1norm_par[x-1][y]) < tol)
-                        assert_(abs(l2norm_par[x][y] - l2norm_par[x-1][y]) < tol)
-                        assert_(abs(linfnorm_par[x][y] - linfnorm_par[x-1][y]) < tol)
-                
-            print 'Parallel test OK'
-
-
Index: /anuga_core/source/anuga_parallel/test_distribute_mesh.py
===================================================================
--- /anuga_core/source/anuga_parallel/test_distribute_mesh.py	(revision 7458)
+++ /anuga_core/source/anuga_parallel/test_distribute_mesh.py	(revision 7459)
@@ -40,5 +40,8 @@
 
     def test_pmesh_1(self):
-
+        """
+        test distributing with just one processor
+        """
+        
         points, vertices, boundary = rectangular_cross(2,2)
 
@@ -81,5 +84,7 @@
 
     def test_pmesh_2(self):
-
+        """
+        Test 2 way pmesh
+        """
         points, vertices, boundary = rectangular_cross(2,2)
 
@@ -121,227 +126,9 @@
 
 
-    def test_distibute_3(self):
-        """
-        Do a parallel test of distributing a rectangle onto 3 processors
-        """
-
-        import pypar
-
-        myid = pypar.rank()
-        numprocs = pypar.size()
-
-        assert numprocs == 3, 'Should be run on 3 processors'
-
-        if myid == 0:
-
-            points, vertices, boundary = rectangular_cross(2,2)
-
-            domain = Domain(points, vertices, boundary)
-
-
-            domain.set_quantity('elevation', topography) # Use function for elevation
-            domain.set_quantity('friction', 0.0)         # Constant friction 
-            domain.set_quantity('stage', expression='elevation') # Dry initial stage
-            domain.set_quantity('xmomentum', expression='friction + 2.0') # 
-            domain.set_quantity('ymomentum', ycoord) #
-
-            #----------------------------------------------------------------------------------
-            # Test pmesh_divide_metis
-            #----------------------------------------------------------------------------------
-            nodes, triangles, boundary, triangles_per_proc, quantities = pmesh_divide_metis(domain,numprocs)
-
-            assert num.allclose(nodes,points)
-
-            true_vertices = [[0, 9, 1], [3, 9, 0], [4, 9, 3], [1, 9, 4], [1, 10, 2], [4, 10, 1], [5, 10, 4], [2, 10, 5], [3, 11, 4], [6, 11, 3], [7, 11, 6], [4, 11, 7], [4, 12, 5], [7, 12, 4], [8, 12, 7], [5, 12, 8]]
-
-            true_triangles = [[4, 9, 3], [4, 12, 5], [7, 12, 4], [8, 12, 7], [5, 12, 8], [0, 9, 1], [1, 9, 4], [1, 10, 2], [4, 10, 1], [5, 10, 4], [2, 10, 5], [3, 9, 0], [3, 11, 4], [6, 11, 3], [7, 11, 6], [4, 11, 7]]
-
-            assert num.allclose(vertices,true_vertices)
-            assert num.allclose(triangles,true_triangles)
-
-            assert num.allclose(triangles_per_proc,[5,6,5])
-
-
-            #----------------------------------------------------------------------------------
-            # Test build_submesh
-            #----------------------------------------------------------------------------------
-            submesh = build_submesh(nodes, triangles, boundary, quantities, triangles_per_proc)
-
-
-            assert num.allclose(submesh['full_nodes'][0],[[3.0, 0.5, 0.0], [4.0, 0.5, 0.5], [5.0, 0.5, 1.0], [7.0, 1.0, 0.5], [8.0, 1.0, 1.0], [9.0, 0.25, 0.25], [12.0, 0.75, 0.75]])
-            assert num.allclose(submesh['full_nodes'][1],[[0.0, 0.0, 0.0], [1.0, 0.0, 0.5], [2.0, 0.0, 1.0], [4.0, 0.5, 0.5], [5.0, 0.5, 1.0], [9.0, 0.25, 0.25], [10.0, 0.25, 0.75]])
-            assert num.allclose(submesh['full_nodes'][2],[[0.0, 0.0, 0.0], [3.0, 0.5, 0.0], [4.0, 0.5, 0.5], [6.0, 1.0, 0.0], [7.0, 1.0, 0.5], [9.0, 0.25, 0.25], [11.0, 0.75, 0.25]])
-
-
-            assert num.allclose(submesh['ghost_nodes'][0],[[0.0, 0.0, 0.0], [1.0, 0.0, 0.5], [2.0, 0.0, 1.0], [6.0, 1.0, 0.0], [10.0, 0.25, 0.75], [11.0, 0.75, 0.25]])
-            assert num.allclose(submesh['ghost_nodes'][1],[[3.0, 0.5, 0.0], [7.0, 1.0, 0.5], [8.0, 1.0, 1.0], [11.0, 0.75, 0.25], [12.0, 0.75, 0.75]])
-            assert num.allclose(submesh['ghost_nodes'][2],[[1.0, 0.0, 0.5], [5.0, 0.5, 1.0], [8.0, 1.0, 1.0], [12.0, 0.75, 0.75]])
-
-
-
-            true_full_triangles = [num.array([[ 4,  9,  3],
-                                              [ 4, 12,  5],
-                                              [ 7, 12,  4],
-                                              [ 8, 12,  7],
-                                              [ 5, 12,  8]]),
-                                   num.array([[ 0,  9,  1],
-                                              [ 1,  9,  4],
-                                              [ 1, 10,  2],
-                                              [ 4, 10,  1],
-                                              [ 5, 10,  4],
-                                              [ 2, 10,  5]]),
-                                   num.array([[ 3,  9,  0],
-                                              [ 3, 11,  4],
-                                              [ 6, 11,  3],
-                                              [ 7, 11,  6],
-                                              [ 4, 11,  7]])]
-
-
-            assert num.allclose(submesh['full_triangles'][0],true_full_triangles[0])
-            assert num.allclose(submesh['full_triangles'][1],true_full_triangles[1])
-            assert num.allclose(submesh['full_triangles'][2],true_full_triangles[2])
-
-            true_ghost_triangles = [num.array([[ 5,  0,  9,  1],
-                                               [ 6,  1,  9,  4],
-                                               [ 8,  4, 10,  1],
-                                               [ 9,  5, 10,  4],
-                                               [10,  2, 10,  5],
-                                               [11,  3,  9,  0],
-                                               [12,  3, 11,  4],
-                                               [13,  6, 11,  3],
-                                               [14,  7, 11,  6],
-                                               [15,  4, 11,  7]]),
-                                    num.array([[ 0,  4,  9,  3],
-                                               [ 1,  4, 12,  5],
-                                               [ 2,  7, 12,  4],
-                                               [ 4,  5, 12,  8],
-                                               [11,  3,  9,  0],
-                                               [12,  3, 11,  4]]),
-                                    num.array([[ 0,  4,  9,  3],
-                                               [ 1,  4, 12,  5],
-                                               [ 2,  7, 12,  4],
-                                               [ 3,  8, 12,  7],
-                                               [ 5,  0,  9,  1],
-                                               [ 6,  1,  9,  4]])]
-
-
-
-            assert num.allclose(submesh['ghost_triangles'][0],true_ghost_triangles[0])
-            assert num.allclose(submesh['ghost_triangles'][1],true_ghost_triangles[1])
-            assert num.allclose(submesh['ghost_triangles'][2],true_ghost_triangles[2])
-
-            true_full_commun = [{0: [1, 2], 1: [1, 2], 2: [1, 2], 3: [2], 4: [1]}, {5: [0, 2], 6: [0, 2], 7: [], 8: [0], 9: [0], 10: [0]}, {11: [0, 1], 12: [0, 1], 13: [0], 14: [0], 15: [0]}]
-
-            assert true_full_commun == submesh['full_commun']
-
-
-            true_ghost_commun = [num.array([[ 5,  1],
-                                            [ 6,  1],
-                                            [ 8,  1],
-                                            [ 9,  1],
-                                            [10,  1],
-                                            [11,  2],
-                                            [12,  2],
-                                            [13,  2],
-                                            [14,  2],
-                                            [15,  2]]),
-                                 num.array([[ 0,  0],
-                                            [ 1,  0],
-                                            [ 2,  0],
-                                            [ 4,  0],
-                                            [11,  2],
-                                            [12,  2]]),
-                                 num.array([[0, 0],
-                                            [1, 0],
-                                            [2, 0],
-                                            [3, 0],
-                                            [5, 1],
-                                            [6, 1]])]
-
-            assert num.allclose(submesh['ghost_commun'][0],true_ghost_commun[0])
-            assert num.allclose(submesh['ghost_commun'][1],true_ghost_commun[1])
-            assert num.allclose(submesh['ghost_commun'][2],true_ghost_commun[2])
-
-
-            #----------------------------------------------------------------------------------
-            # Test send_submesh
-            #----------------------------------------------------------------------------------
-            for p in range(1, numprocs):
-                send_submesh(submesh, triangles_per_proc, p, verbose=False)
-
-            #----------------------------------------------------------------------------------
-            # Test extract_hostmesh
-            #----------------------------------------------------------------------------------
-            points, vertices, boundary, quantities, ghost_recv_dict, full_send_dict  =\
-            extract_hostmesh(submesh, triangles_per_proc)
-
-
-            true_points =  [[0.5, 0.0], [0.5, 0.5], [0.5, 1.0], [1.0, 0.5], [1.0, 1.0], [0.25, 0.25], [0.75, 0.75], [0.0, 0.0], [0.0, 0.5], [0.0, 1.0], [1.0, 0.0], [0.25, 0.75], [0.75, 0.25]]
-
-            true_vertices = [[1, 5, 0], [1, 6, 2], [3, 6, 1], [4, 6, 3], [2, 6, 4], [7, 5, 8], [8, 5, 1], [1, 11, 8], [2, 11, 1], [9, 11, 2], [0, 5, 7], [0, 12, 1], [10, 12, 0], [3, 12, 10], [1, 12, 3]]
-
-
-            true_ghost_recv = {1: [num.array([5, 6, 7, 8, 9]), num.array([ 5,  6,  8,  9, 10])], 2: [num.array([10, 11, 12, 13, 14]), num.array([11, 12, 13, 14, 15])]}
-
-
-            true_full_send = {1: [num.array([0, 1, 2, 4]), num.array([0, 1, 2, 4])], 2: [num.array([0, 1, 2, 3]), num.array([0, 1, 2, 3])]}
-
-            assert num.allclose(points,   true_points)
-            assert num.allclose(vertices, true_vertices)
-            assert num.allclose(ghost_recv_dict[1],true_ghost_recv[1])
-            assert num.allclose(ghost_recv_dict[2],true_ghost_recv[2])
-            assert num.allclose(full_send_dict[1],true_full_send[1])
-            assert num.allclose(full_send_dict[2],true_full_send[2])
-
-            #print triangles_per_proc
-
-        else:
-            #----------------------------------------------------------------------------------
-            # Test rec_submesh
-            #----------------------------------------------------------------------------------
-            points, vertices, boundary, quantities, ghost_recv_dict, full_send_dict, no_full_nodes, no_full_trigs = rec_submesh(0, verbose=False)    
-
-            if myid == 1:
-
-
-                true_points =  [[0.0, 0.0], [0.0, 0.5], [0.0, 1.0], [0.5, 0.5], [0.5, 1.0], [0.25, 0.25], [0.25, 0.75], [0.5, 0.0], [1.0, 0.5], [1.0, 1.0], [0.75, 0.25], [0.75, 0.75]] 
-
-                true_vertices =  [[0, 5, 1], [1, 5, 3], [1, 6, 2], [3, 6, 1], [4, 6, 3], [2, 6, 4], [3, 5, 7], [3, 11, 4], [8, 11, 3], [4, 11, 9], [7, 5, 0], [7, 10, 3]]
-
-                true_ghost_recv =  {0: [num.array([6, 7, 8, 9]), num.array([0, 1, 2, 4])], 2: [num.array([10, 11]), num.array([11, 12])]}
-
-                true_full_send =  {0: [num.array([0, 1, 3, 4, 5]), num.array([ 5,  6,  8,  9, 10])], 2: [num.array([0, 1]), num.array([5, 6])]}
-
-                assert num.allclose(points,   true_points)
-                assert num.allclose(vertices, true_vertices)
-                assert num.allclose(ghost_recv_dict[0],true_ghost_recv[0])
-                assert num.allclose(ghost_recv_dict[2],true_ghost_recv[2])
-                assert num.allclose(full_send_dict[0],true_full_send[0])
-                assert num.allclose(full_send_dict[2],true_full_send[2]) 
-
-
-            if myid == 2:
-
-                true_points =   [[0.0, 0.0], [0.5, 0.0], [0.5, 0.5], [1.0, 0.0], [1.0, 0.5], [0.25, 0.25], [0.75, 0.25], [0.0, 0.5], [0.5, 1.0], [1.0, 1.0], [0.75, 0.75]]
-
-                true_vertices =  [[1, 5, 0], [1, 6, 2], [3, 6, 1], [4, 6, 3], [2, 6, 4], [2, 5, 1], [2, 10, 8], [4, 10, 2], [9, 10, 4], [0, 5, 7], [7, 5, 2]]
-
-
-                true_ghost_recv =   {0: [num.array([5, 6, 7, 8]), num.array([0, 1, 2, 3])], 1: [num.array([ 9, 10]), num.array([5, 6])]}
-
-                true_full_send =   {0: [num.array([0, 1, 2, 3, 4]), num.array([11, 12, 13, 14, 15])], 1: [num.array([0, 1]), num.array([11, 12])]}
-
-
-                assert num.allclose(points,   true_points)
-                assert num.allclose(vertices, true_vertices)
-                assert num.allclose(ghost_recv_dict[0],true_ghost_recv[0])
-                assert num.allclose(ghost_recv_dict[1],true_ghost_recv[1])
-                assert num.allclose(full_send_dict[0],true_full_send[0])
-                assert num.allclose(full_send_dict[1],true_full_send[1]) 
-
-
             
     def test_build_submesh_3(self):
-
+        """
+        Test 3 way build_submesh
+        """
 
         nodes = [[0.0, 0.0], [0.0, 0.5], [0.0, 1.0], [0.5, 0.0], [0.5, 0.5], [0.5, 1.0], [1.0, 0.0], [1.0, 0.5], [1.0, 1.0], [0.25, 0.25], [0.25, 0.75], [0.75, 0.25], [0.75, 0.75]]
Index: /anuga_core/source/anuga_parallel/test_parallel_distribute_domain.py
===================================================================
--- /anuga_core/source/anuga_parallel/test_parallel_distribute_domain.py	(revision 7459)
+++ /anuga_core/source/anuga_parallel/test_parallel_distribute_domain.py	(revision 7459)
@@ -0,0 +1,223 @@
+
+"""Test a run of the sequential shallow water domain against
+a run of the parallel shallow water domain.
+
+WARNING: This assumes that the command to run jobs is mpirun.
+Tested with MPICH and LAM (Ole)
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+
+import unittest
+import os
+import sys
+import pypar
+
+import numpy as num
+
+
+
+
+from anuga.utilities.numerical_tools import ensure_numeric
+from anuga.utilities.util_ext        import double_precision
+from anuga.utilities.norms           import l1_norm, l2_norm, linf_norm
+
+from anuga.interface import Domain
+from anuga.interface import Reflective_boundary
+from anuga.interface import Dirichlet_boundary
+from anuga.interface import Time_boundary
+from anuga.interface import Transmissive_boundary
+
+from anuga.interface import rectangular_cross
+from anuga.interface import create_domain_from_file
+
+
+from anuga_parallel.interface import distribute, myid, numprocs
+
+
+#--------------------------------------------------------------------------
+# Setup parameters
+#--------------------------------------------------------------------------
+
+mesh_filename = "merimbula_10785_1.tsh"
+#mesh_filename = "test-100.tsh"
+yieldstep = 1
+finaltime = 20
+quantity = 'stage'
+nprocs = 4
+verbose = True
+
+#--------------------------------------------------------------------------
+# Setup procedures
+#--------------------------------------------------------------------------
+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))
+
+#--------------------------------------------------------------------------
+# Setup test
+#--------------------------------------------------------------------------
+def evolution_test(parallel=False):
+
+
+    domain = create_domain_from_file(mesh_filename)
+    domain.set_quantity('stage', Set_Stage(756000.0, 756500.0, 2.0))
+
+    #--------------------------------------------------------------------------
+    # Create parallel domain if requested
+    #--------------------------------------------------------------------------
+
+    if parallel:
+        if myid == 0 and verbose: print 'DISTRIBUTING PARALLEL DOMAIN'
+        domain = distribute(domain)
+
+    #------------------------------------------------------------------------------
+    # Setup boundary conditions
+    # This must currently happen *after* domain has been distributed
+    #------------------------------------------------------------------------------
+    domain.store = False
+    Br = Reflective_boundary(domain)      # Solid reflective wall
+
+    domain.set_boundary({'outflow' :Br, 'inflow' :Br, 'inner' :Br, 'exterior' :Br, 'open' :Br})
+
+    #------------------------------------------------------------------------------
+    # Setup diagnostic arrays
+    #------------------------------------------------------------------------------
+    l1list = []
+    l2list = []
+    linflist = []
+    l1norm = num.zeros(3, num.float)
+    l2norm = num.zeros(3, num.float)
+    linfnorm = num.zeros(3, num.float)
+    recv_norm = num.zeros(3, num.float)
+
+    #------------------------------------------------------------------------------
+    # Evolution
+    #------------------------------------------------------------------------------
+    if parallel:
+        if myid == 0 and verbose: print 'PARALLEL EVOLVE'
+    else:
+        if verbose: print 'SEQUENTIAL EVOLVE'
+        
+    for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
+        edges = domain.quantities[quantity].edge_values.take(num.flatnonzero(domain.tri_full_flag),axis=0)
+        l1norm[0] = l1_norm(edges[:,0])
+        l1norm[1] = l1_norm(edges[:,1])
+        l1norm[2] = l1_norm(edges[:,2])
+        l2norm[0] = l2_norm(edges[:,0])
+        l2norm[1] = l2_norm(edges[:,1])
+        l2norm[2] = l2_norm(edges[:,2])
+        linfnorm[0] = linf_norm(edges[:,0])
+        linfnorm[1] = linf_norm(edges[:,1])
+        linfnorm[2] = linf_norm(edges[:,2])
+        if parallel:
+            l2norm[0] = pow(l2norm[0], 2)
+            l2norm[1] = pow(l2norm[1], 2)
+            l2norm[2] = pow(l2norm[2], 2)
+            if myid == 0:
+                #domain.write_time()
+
+                #print edges[:,1]            
+                for p in range(1, numprocs):
+                    recv_norm = pypar.receive(p)
+                    l1norm += recv_norm
+                    recv_norm = pypar.receive(p)
+                    l2norm += recv_norm
+                    recv_norm = pypar.receive(p)
+                    linfnorm[0] = max(linfnorm[0], recv_norm[0])
+                    linfnorm[1] = max(linfnorm[1], recv_norm[1])
+                    linfnorm[2] = max(linfnorm[2], recv_norm[2])
+
+                l2norm[0] = pow(l2norm[0], 0.5)
+                l2norm[1] = pow(l2norm[1], 0.5)
+                l2norm[2] = pow(l2norm[2], 0.5)
+
+                l1list.append(l1norm)                
+                l2list.append(l2norm)
+                linflist.append(linfnorm)                
+            else:
+                pypar.send(l1norm, 0)
+                pypar.send(l2norm, 0)
+                pypar.send(linfnorm, 0)
+        else:
+            #domain.write_time()
+            l1list.append(l1norm)                
+            l2list.append(l2norm)
+            linflist.append(linfnorm)
+            
+
+    return (l1list, l2list, linflist)
+
+# Test an nprocs-way run of the shallow water equations
+# against the sequential code.
+
+class Test_parallel_distribute_domain(unittest.TestCase):
+    def test_parallel_distribute_domain(self):
+        print "Expect this test to fail if not run from the parallel directory."
+        result = os.system("mpirun -np %d python test_parallel_distribute_domain.py" % nprocs)
+        assert_(result == 0)
+
+
+# Because we are doing assertions outside of the TestCase class
+# the PyUnit defined assert_ function can't be used.
+def assert_(condition, msg="Assertion Failed"):
+    if condition == False:
+        #pypar.finalize()
+        raise AssertionError, msg
+
+if __name__=="__main__":
+    if numprocs == 1: 
+        runner = unittest.TextTestRunner()
+        suite = unittest.makeSuite(Test_parallel_distribute_domain, 'test')
+        runner.run(suite)
+    else:
+
+        pypar.barrier()
+        if myid == 0:
+            if verbose: print 'SEQUENTIAL START'
+            l1norm_seq, l2norm_seq, linfnorm_seq = evolution_test(parallel=False)
+
+        pypar.barrier()
+        if myid ==0:
+            if verbose: print 'PARALLEL START'
+        
+        l1norm_par, l2norm_par, linfnorm_par = evolution_test(parallel=True)
+        
+        if myid == 0:
+            assert_(len(l1norm_seq) == len(l1norm_par))
+            assert_(len(l2norm_seq) == len(l2norm_par))
+            assert_(len(linfnorm_seq) == len(linfnorm_par))
+            assert_(len(l1norm_seq) == len(l2norm_seq))
+            assert_(len(l2norm_seq) == len(linfnorm_seq))
+            # Anything smaller than tol we consider to be 0.
+            # This usualy comes out to be 10^-14 (DBL_DIG = 10^-15)
+            # * 10 to cater for rounding error in computation.
+            tol = pow(10, -1 * (double_precision() - 1))
+            for x in range(len(l1norm_seq)):
+                for y in range(3):
+                    # Calculate relative difference in the norms
+                    assert_(abs(l1norm_seq[x][y] - l1norm_par[x][y])/l1norm_seq[x][y] < tol)
+                    assert_(abs(l2norm_seq[x][y] - l2norm_par[x][y])/l2norm_seq[x][y] < tol)
+                    assert_(abs(linfnorm_seq[x][y] - linfnorm_par[x][y])/linfnorm_seq[x][y] < tol)
+                    if x > 0:
+                        # Verify that the quantity is being conserved across iterations.
+                        assert_(abs(l1norm_seq[x][y] - l1norm_seq[x-1][y]) < tol)
+                        assert_(abs(l2norm_seq[x][y] - l2norm_seq[x-1][y]) < tol)
+                        assert_(abs(linfnorm_seq[x][y] - linfnorm_seq[x-1][y]) < tol)
+                        assert_(abs(l1norm_par[x][y] - l1norm_par[x-1][y]) < tol)
+                        assert_(abs(l2norm_par[x][y] - l2norm_par[x-1][y]) < tol)
+                        assert_(abs(linfnorm_par[x][y] - linfnorm_par[x-1][y]) < tol)
+                
+            if verbose: print 'Parallel test OK'
+
+
