Index: /inundation/parallel/pmesh_divide.py
===================================================================
--- /inundation/parallel/pmesh_divide.py	(revision 2054)
+++ /inundation/parallel/pmesh_divide.py	(revision 2055)
@@ -20,5 +20,5 @@
 
 from math import floor
-from Numeric import zeros, Float, Int
+from Numeric import zeros, Float, Int, reshape
 
 #########################################################
@@ -151,5 +151,4 @@
         for i in range(N):
             q_reord[k][index[i]]=quantities[k].vertex_values[i]
-
     del index
 
@@ -268,5 +267,5 @@
     tri_index = {}
     N = domain.number_of_elements
-
+ 
     #sort by x coordinate of centroid
     from Numeric import argsort
@@ -286,5 +285,4 @@
 
     #print tri_list
-
     #print tri_index
 
@@ -321,2 +319,78 @@
 
     return nodes, triangles, boundary, triangles_per_proc, quantities
+
+
+from os import sep
+from sys import path
+
+path.append('..' + sep + 'pymetis')
+
+from metis import partMeshNodal
+
+def pmesh_divide_metis(domain, n_procs):
+    # initialise the lists
+    triangles_per_proc = []
+    tri_list = []
+    proc_sum = []
+    for i in range(n_procs):
+        tri_list.append([])
+        triangles_per_proc.append(0)
+        proc_sum.append([])
+
+    # Prepare variables for the metis call
+    n_tri = len(domain.triangles)
+    n_vert = len(domain.coordinates)
+    t_list = domain.triangles.copy()
+    t_list = reshape(t_list, (-1,))
+
+    # The 1 here is for triangular mesh elements.
+    edgecut, epart, npart = partMeshNodal(n_tri, n_vert, t_list, 1, n_procs)
+
+    # print edgecut
+    # print epart
+    # print npart
+
+    del edgecut
+    del npart
+
+    # Assign triangles to processes, according to what metis told us.
+    
+    tri_index = {}
+    triangles = []        
+    for i in range(n_tri):
+        triangles_per_proc[epart[i]] = triangles_per_proc[epart[i]] + 1
+        tri_list[epart[i]].append(domain.triangles[i])
+        tri_index[i] = ([epart[i], len(tri_list[epart[i]]) - 1])
+
+    # print tri_list
+    # print triangles_per_proc
+
+    # order the triangle list so that all of the triangles belonging
+    # to processor i are listed before those belonging to processor
+    # i+1
+
+    for i in range(n_procs):
+        for t in tri_list[i]:
+            triangles.append(t)
+            
+    # the boundary labels have to changed in accoradance with the
+    # new triangle ordering, proc_sum and tri_index help with this
+
+    proc_sum[0] = 0
+    for i in range(n_procs - 1):
+        proc_sum[i+1]=proc_sum[i]+triangles_per_proc[i]
+
+    # relabel the boundary elements to fit in with the new triangle
+    # ordering
+
+    boundary = {}
+    for b in domain.boundary:
+        t =  tri_index[b[0]]
+        boundary[proc_sum[t[0]]+t[1], b[1]] = domain.boundary[b]
+
+    quantities = reorder(domain.quantities, tri_index, proc_sum)
+
+    # extract the node list
+    nodes = domain.coordinates.copy()
+
+    return nodes, triangles, boundary, triangles_per_proc, quantities
Index: /inundation/parallel/run_parallel_sw_merimbula_metis.py
===================================================================
--- /inundation/parallel/run_parallel_sw_merimbula_metis.py	(revision 2055)
+++ /inundation/parallel/run_parallel_sw_merimbula_metis.py	(revision 2055)
@@ -0,0 +1,198 @@
+#!/usr/bin/env python
+#########################################################
+#
+#  Main file for parallel mesh testing.
+#
+#  This is a modification of the run_parallel_advection.py
+# file.
+#
+#  *) The test files currently avaliable are of the form
+# test*.out, eg test_5l_4c.out. The term infront of the l
+# corresponds to the number of levels of refinement
+# required to build the grid, i.e. a higher number
+# corresponds to a finer grid. The term infront of the c
+# corresponds to the number of processors.
+#
+# *) The (new) files that have been added to manage the
+# grid partitioning are
+#    +) mg2ga.py: read in the test files.
+#    +) pmesh_divide.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
+#
+# *) Things still to do:
+#    +) Overlap the communication and computation: The
+# communication routines in build_commun.py should be
+# interdispersed in the build_submesh.py and build_local.py
+# files. This will overlap the communication and
+# computation and will be far more efficient. This should
+# be done after more testing and there more confidence in
+# the subpartioning.
+#    +) Much more testing especially with large numbers of
+# 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')
+
+from Numeric import array
+# pmesh
+
+#from shallow_water import Domain
+
+from shallow_water import Domain
+from parallel_shallow_water import Parallel_Domain
+
+# mesh partition routines
+
+from pmesh_divide import pmesh_divide_metis
+from build_submesh import *
+from build_local import *
+from build_commun import *
+from pmesh2domain import pmesh_to_domain_instance
+
+# read in the processor information
+
+numprocs = pypar.size()
+myid = pypar.rank()
+processor_name = pypar.Get_processor_name()
+
+#-------
+# Domain
+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))
+
+
+if myid == 0:
+
+    # read in the test files
+
+#    filename = 'test-100.tsh'
+    filename = 'merimbula_10785_1.tsh'
+
+    domain_full = pmesh_to_domain_instance(filename, Domain)
+
+#    domain_full.set_quantity('stage', Set_Stage(200.0,300.0,1.0))
+    domain_full.set_quantity('stage', Set_Stage(756000.0,756500.0,2.0))
+
+    nodes, triangles, boundary, triangles_per_proc, quantities = \
+         pmesh_divide_metis(domain_full, numprocs)
+
+    rect = array(domain_full.xy_extent, Float)
+
+    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)
+
+    hostmesh = extract_hostmesh(submesh)
+    points, vertices, boundary, quantities, ghost_recv_dict, full_send_dict = \
+             build_local_mesh(hostmesh, 0, triangles_per_proc[0], numprocs)
+
+# read in the mesh partition that belongs to this
+# processor (note that the information is in the
+# correct form for the GA data structure
+
+else:
+    points, vertices, boundary, quantities, ghost_recv_dict, full_send_dict \
+            = rec_submesh(0)
+
+#if myid == 0:
+#    print 'ghost'
+#    print ghost_recv_dict
+#processor_name
+#if myid == 0:
+#    print 'full'
+#    print full_send_dict
+
+
+pypar.broadcast(rect,0)
+#print rect
+
+domain = Parallel_Domain(points, vertices, boundary,
+                         full_send_dict  = full_send_dict,
+                         ghost_recv_dict = ghost_recv_dict)
+
+
+try:
+    domain.initialise_visualiser(rect=rect)
+    #domain.visualiser.coloring['stage'] = True
+    domain.visualiser.scale_z['stage'] = 0.2
+    domain.visualiser.scale_z['elevation'] = 0.05
+except:
+    print 'No visualiser'
+
+
+domain.default_order = 1
+
+#Boundaries
+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} )
+
+
+domain.set_quantity('stage', quantities['stage'])
+domain.set_quantity('elevation', quantities['elevation'])
+
+#domain.store = True
+#domain.filename = 'merimbula-%d' %domain.processor
+
+#---------
+# Evolution
+t0 = time.time()
+
+print 'Processor %d on %s: No of elements %d'%(domain.processor,processor_name,domain.number_of_elements)
+yieldstep = 0.05
+finaltime = 500.0
+
+#yieldstep = 1
+#finaltime = 1
+#processor_name
+for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
+    if myid == 0:
+        domain.write_time()
+        #print 'Processor %d, Integral of stage %d'%\
+        #       (domain.processor,domain.quantities['stage'].get_integral())
+
+
+#print 'P%d: That took %.2f seconds' %(myid, time.time()-t0)
+#print 'P%d: Communication time %.2f seconds' %(myid, domain.communication_time)
+#print 'P%d: Reduction Communication time %.2f seconds' %(myid, domain.communication_reduce_time)
+#print 'P%d: Broadcast time %.2f seconds' %(myid, domain.communication_broadcast_time)
+
+
+
+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
