Index: uga_core/source/anuga_parallel/build_commun.py
===================================================================
--- /anuga_core/source/anuga_parallel/build_commun.py	(revision 7446)
+++ 	(revision )
@@ -1,352 +1,0 @@
-#########################################################
-#
-# Handle the communication between the host machine
-# (processor 0) and the processors. The host machine is
-# responsible for the doing the initial grid partitioning.
-#
-# The routines given below should be moved to the
-# build_submesh.py and build_local.py file to allow
-# overlapping of  communication and computation.
-# This should be done after more debugging.
-#
-#
-#  Author: Linda Stals, June 2005
-#  Modified: Linda Stals, Nov 2005 (optimise python code)
-#            Steve Roberts, Aug 2009 (update to numpy)
-#
-#
-#########################################################
-
-#from Numeric import array, Int, Float, zeros
-
-import numpy as num
-
-import logging, logging.config
-logger = logging.getLogger('parallel')
-logger.setLevel(logging.WARNING)
-
-try:
-    logging.config.fileConfig('log.ini')
-except:
-    pass
-
-import sys
-
-import pypar
-
-
-from build_local import build_local_mesh
-
-#########################################################
-#
-# Send the submesh to processor p.
-#
-# *) The order and form is strongly coupled with
-# rec_submesh.
-#
-# -------------------------------------------------------
-#
-# *) All of the information has been sent to processor p.
-#
-#########################################################
-
-def send_submesh(submesh, triangles_per_proc, p):
-
-    myid = pypar.rank()
-    print 'process %d sending submesh to process %d' %(myid, p)
-    
-    # build and send the tagmap for the boundary conditions
-    
-    tagmap = {}
-    counter = 1
-    for b in submesh["full_boundary"][p]:
-         bkey = submesh["full_boundary"][p][b]
-         if not tagmap.has_key(bkey):
-             tagmap[bkey] = counter
-             counter = counter+1
-    for b in submesh["ghost_boundary"][p]:
-         bkey = submesh["ghost_boundary"][p][b]
-         if not tagmap.has_key(bkey):
-             tagmap[bkey] = counter
-             counter = counter+1
-
-    pypar.send(tagmap, p)
-
-    # send the quantities key information
-
-    pypar.send(submesh["full_quan"].keys(), p)
-    
-    # send the number of triangles per processor
-
-    pypar.send(triangles_per_proc, p)
-
-    # compress full_commun
-
-    flat_full_commun = []
-
-    for c in submesh["full_commun"][p]:
-        for i in range(len(submesh["full_commun"][p][c])):
-            flat_full_commun.append([c,submesh["full_commun"][p][c][i]])
-
-    print 'flat_full_commun', flat_full_commun
-
-    # send the array sizes so memory can be allocated
-
-    setup_array = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
-    setup_array[0] = len(submesh["full_nodes"][p])
-    setup_array[1] = len(submesh["ghost_nodes"][p])
-    setup_array[2] = len(submesh["full_triangles"][p])
-    setup_array[3] = len(submesh["ghost_triangles"][p])
-    setup_array[4] = len(submesh["full_boundary"][p])
-    setup_array[5] = len(submesh["ghost_boundary"][p])
-    setup_array[6] = len(submesh["ghost_commun"][p])
-    setup_array[7] = len(flat_full_commun)
-
-    pypar.send(num.array(setup_array, num.int), p)
-    
-    # send the nodes
-
-    pypar.send(num.array(submesh["full_nodes"][p], num.float), p)
-    pypar.send(num.array(submesh["ghost_nodes"][p], num.float),p)
-
-    # send the triangles
-
-    pypar.send(num.array(submesh["full_triangles"][p],  num.int), p)
-    pypar.send(num.array(submesh["ghost_triangles"][p], num.int), p)
-
-    # send the boundary
-
-    bc = []
-    for b in submesh["full_boundary"][p]:
-        bc.append([b[0], b[1], tagmap[submesh["full_boundary"][p][b]]])
-
-
-    pypar.send(num.array(bc, num.int), p)
-
-    bc = []
-    for b in submesh["ghost_boundary"][p]:
-        bc.append([b[0], b[1], tagmap[submesh["ghost_boundary"][p][b]]])
-
-    pypar.send(num.array(bc, num.int), p)
-
-    # send the communication pattern
-
-    pypar.send(submesh["ghost_commun"][p], p)
-
-    pypar.send(num.array(flat_full_commun, num.int), p)
-
-    # send the quantities
-    
-    for k in submesh["full_quan"]:
-        pypar.send(num.array(submesh["full_quan"][k][p], num.float), p)
-        
-    for k in submesh["ghost_quan"]:
-        pypar.send(num.array(submesh["ghost_quan"][k][p], num.float),p)
-        
-
-#########################################################
-#
-# Receive the submesh from processor p.
-#
-# *) The order and form is strongly coupled with
-# send_submesh.
-#
-# -------------------------------------------------------
-#
-# *) All of the information has been received by the
-# processor p and passed into build_local.
-#
-# *) The information is returned in a form needed by the
-# GA datastructure.
-#
-#########################################################
-
-def rec_submesh_flat(p):
-    
-    numproc = pypar.size()
-    myid = pypar.rank()
-
-    submesh_cell = {}
-    
-    print 'process %d receiving submesh from process %d' %(myid, p)
-
-    # receive the tagmap for the boundary conditions
-    
-    tagmap = pypar.receive(p)
-
-    itagmap = {}
-    for t in tagmap:
-        itagmap[tagmap[t]]=t
-
-    # receive the quantities key information
-
-    qkeys = pypar.receive(p)
-
-    # receive the number of triangles per processor
-
-    triangles_per_proc = pypar.receive(p)
-
-    # recieve information about the array sizes
-
-    setup_array = pypar.receive(p)
-
-    no_full_nodes = setup_array[0]
-    no_ghost_nodes = setup_array[1]
-    no_full_triangles = setup_array[2]
-    no_ghost_triangles = setup_array[3]
-    no_full_boundary = setup_array[4]
-    no_ghost_boundary = setup_array[5]
-    no_ghost_commun = setup_array[6]
-    no_full_commun = setup_array[7]
-    no_quantities = len(qkeys)
-    
-    # receive the full nodes
-
-    submesh_cell["full_nodes"] = pypar.receive(p)
-
-    # receive the ghost nodes
-
-    submesh_cell["ghost_nodes"] = pypar.receive(p)
-    
-    # receive the full triangles
-
-    submesh_cell["full_triangles"] = pypar.receive(p)
-    
-    # receive the ghost triangles
-
-    submesh_cell["ghost_triangles"] = pypar.receive(p)
-
-    # receive the full boundary
-
-    bnd_c = pypar.receive(p)
-
-    submesh_cell["full_boundary"] = {}
-    for b in bnd_c:
-        submesh_cell["full_boundary"][b[0],b[1]]=itagmap[b[2]]
-
-    # receive the ghost boundary
-
-    bnd_c = pypar.receive(p)
-
-    submesh_cell["ghost_boundary"] = {}
-    for b in bnd_c:
-        submesh_cell["ghost_boundary"][b[0],b[1]]=itagmap[b[2]]
-
-    # receive the ghost communication pattern
-
-    submesh_cell["ghost_commun"] = pypar.receive(p)
-    
-    # receive the full communication pattern
-
-    full_commun = pypar.receive(p)
-
-    submesh_cell["full_commun"] = {}
-    for c in full_commun:
-        submesh_cell["full_commun"][c[0]] = []
-    for c in full_commun:
-        submesh_cell["full_commun"][c[0]].append(c[1])
-
-    # receive the quantities
-
-    submesh_cell["full_quan"]={}
-    
-    for i in range(no_quantities):
-        tmp = pypar.receive(p)
-        submesh_cell["full_quan"][qkeys[i]]=num.zeros((no_full_triangles,3), num.float)
-        submesh_cell["full_quan"][qkeys[i]][:] = tmp[:]
-
-    submesh_cell["ghost_quan"]={}
-    for i in range(no_quantities):
-        tmp = pypar.receive(p)
-        submesh_cell["ghost_quan"][qkeys[i]]= num.zeros((no_ghost_triangles,3), num.float)
-        submesh_cell["ghost_quan"][qkeys[i]][:] = tmp[:]
-    
-    return submesh_cell, triangles_per_proc,\
-           no_full_nodes, no_full_triangles
-
-
-
-#########################################################
-#
-# Receive the submesh from processor p.
-#
-# *) The order and form is strongly coupled with
-# send_submesh.
-#
-# -------------------------------------------------------
-#
-# *) All of the information has been received by the
-# processor p and passed into build_local.
-#
-# *) The information is returned in a form needed by the
-# GA datastructure.
-#
-#########################################################
-
-def rec_submesh(p):
-    
-    numproc = pypar.size()
-    myid = pypar.rank()
-
-    [submesh_cell, triangles_per_proc,\
-     number_of_full_nodes, number_of_full_triangles] = rec_submesh_flat(p)
-    
-    # find the full triangles assigned to this processor
-
-    lower_t = 0
-    for i in range(myid):
-        lower_t = lower_t+triangles_per_proc[i]
-    upper_t = lower_t+triangles_per_proc[myid]
-
-    # convert the information into a form needed by the GA
-    # datastructure
-
-    [GAnodes, GAtriangles, boundary, quantities, ghost_rec, full_send] = \
-              build_local_mesh(submesh_cell, lower_t, upper_t, \
-                               numproc)
-    
-    return GAnodes, GAtriangles, boundary, quantities,\
-           ghost_rec, full_send,\
-           number_of_full_nodes, number_of_full_triangles
-
-
-#########################################################
-#
-# Extract the submesh that will belong to the
-# "host processor" (i.e. processor zero)
-#
-#  *) See the documentation for build_submesh
-#
-# -------------------------------------------------------
-#
-#  *) A dictionary containing the full_triangles,
-# full_nodes, full_boundary, ghost_triangles, ghost_nodes,
-# ghost_boundary, ghost_commun and full_commun belonging
-# to processor zero are returned.
-#
-#########################################################
-def extract_hostmesh(submesh, triangles_per_proc):
-
-    submesh_cell = {}
-    submesh_cell["full_nodes"] = submesh["full_nodes"][0]
-    submesh_cell["ghost_nodes"] = submesh["ghost_nodes"][0]
-    submesh_cell["full_triangles"] = submesh["full_triangles"][0]
-    submesh_cell["ghost_triangles"] = submesh["ghost_triangles"][0]
-    submesh_cell["full_boundary"] = submesh["full_boundary"][0]
-    submesh_cell["ghost_boundary"] = submesh["ghost_boundary"][0]
-    submesh_cell["ghost_commun"] = submesh["ghost_commun"][0]
-    submesh_cell["full_commun"] = submesh["full_commun"][0]
-    submesh_cell["full_quan"] ={}
-    submesh_cell["ghost_quan"]={}
-    for k in submesh["full_quan"]:
-        submesh_cell["full_quan"][k] = submesh["full_quan"][k][0]
-        submesh_cell["ghost_quan"][k] = submesh["ghost_quan"][k][0]
-
-    numprocs = pypar.size()
-    points, vertices, boundary, quantities, ghost_recv_dict, full_send_dict = \
-            build_local_mesh(submesh_cell, 0, triangles_per_proc[0], numprocs)
-    return  points, vertices, boundary, quantities, ghost_recv_dict, \
-           full_send_dict
-
-
-
Index: uga_core/source/anuga_parallel/build_local.py
===================================================================
--- /anuga_core/source/anuga_parallel/build_local.py	(revision 7446)
+++ 	(revision )
@@ -1,228 +1,0 @@
-#########################################################
-#
-#  Given the subdivision of the grid assigned to the
-# current processor convert it into a form that is
-# appropriate for the GA datastructure.
-#
-#  The main function of these modules is to change the
-# node numbering. The GA datastructure assumes they
-# are numbered consecutively from 0.
-#
-#  The module also changes the communication pattern
-# datastructure into a form needed by parallel_advection
-#
-#  Authors: Linda Stals and Matthew Hardy, June 2005
-#  Modified: Linda Stals, Nov 2005 (optimise python code)
-#            Steve Roberts, Aug 2009 (updating to numpy)
-#
-#
-#########################################################
-
-#from Numeric import  zeros, Float, Int, concatenate, \
-#     take, arrayrange, put, sort, compress, equal
-
-
-import numpy as num
- 
-#########################################################
-# Convert the format of the data to that used by ANUGA
-#
-#
-# *) Change the nodes global ID's to an integer value,
-#starting from 0.
-#
-# *) The triangles and boundary edges must also be
-# updated accordingly.
-#
-# -------------------------------------------------------
-#
-# *) The nodes, triangles and boundary edges defined by
-# the new numbering scheme are returned
-#
-#########################################################
-
-def build_local_GA(nodes, triangles, boundaries, tri_index):
-
-    Nnodes =len(nodes)
-    Ntriangles = len(triangles)
-    
-    # Extract the nodes (using the local ID)
-    
-    GAnodes = num.take(nodes, (1, 2), 1)
-
-    # Build a global ID to local ID mapping
-
-    NGlobal = 0
-    for i in range(Nnodes):
-        if nodes[i][0] > NGlobal:
-            NGlobal = nodes[i][0]
-    index = num.zeros(int(NGlobal)+1, num.int)
-    num.put(index, num.take(nodes, (0,), 1).astype(num.int), \
-        num.arange(Nnodes))
-        
-    # Change the global IDs in the triangles to the local IDs
-
-    GAtriangles = num.zeros((Ntriangles, 3), num.int)
-    GAtriangles[:,0] = num.take(index, triangles[:,0])
-    GAtriangles[:,1] = num.take(index, triangles[:,1])
-    GAtriangles[:,2] = num.take(index, triangles[:,2])
-
-    # Change the triangle numbering in the boundaries
-
-    GAboundaries = {}
-    for b in boundaries:
-        GAboundaries[tri_index[b[0]], b[1]] = boundaries[b]
-        
-    del (index)
-    
-    return GAnodes, GAtriangles, GAboundaries
-
-
-#########################################################
-# Change the communication format to that needed by the
-# parallel advection file.
-#
-# *) The index contains [global triangle no,
-# local triangle no.]
-#
-# -------------------------------------------------------
-#
-# *) The ghost_recv and full_send dictionaries are
-# returned.
-#
-# *) ghost_recv dictionary is local id, global id, value
-#
-# *) full_recv dictionary is local id, global id, value
-#
-# *) The information is ordered by the global id. This
-# means that the communication order is predetermined and
-# local and global id do not need to be
-# compared when the information is sent/received.
-#
-#########################################################
-
-def build_local_commun(index, ghostc, fullc, nproc):
-
-    # Initialise
-
-    full_send = {}
-    ghost_recv = {}
-
-    # Build the ghost_recv dictionary (sort the
-    # information by the global numbering)
-    
-    ghostc = num.sort(ghostc, 0)
-    
-    for c in range(nproc):
-        s = ghostc[:,0]
-        d = num.compress(num.equal(ghostc[:,1],c), s)
-        if len(d) > 0:
-            ghost_recv[c] = [0, 0]
-            ghost_recv[c][1] = d
-            ghost_recv[c][0] = num.take(index, d)
-            
-    # Build a temporary copy of the full_send dictionary
-    # (this version allows the information to be stored
-    # by the global numbering)
-
-    tmp_send = {}
-    for global_id in fullc:
-        for i in range(len(fullc[global_id])):
-            neigh = fullc[global_id][i]
-            if not tmp_send.has_key(neigh):
-                tmp_send[neigh] = []
-            tmp_send[neigh].append([global_id, \
-                                    index[global_id]])
-
-    # Extract the full send information and put it in the form
-    # required for the full_send dictionary
-
-    for neigh in tmp_send:
-        neigh_commun = num.sort(tmp_send[neigh], 0)
-        full_send[neigh] = [0, 0]
-        full_send[neigh][0] = neigh_commun[:,1]
-        full_send[neigh][1] = neigh_commun[:,0]
-
-    return ghost_recv, full_send
-
-
-#########################################################
-# Convert the format of the data to that used by ANUGA
-#
-#
-# *) Change the nodes global ID's to an integer value,
-# starting from 0. The node numbering in the triangles
-# must also be updated to take this into account.
-#
-# *) The triangle number will also change, which affects
-# the boundary tag information and the communication
-# pattern.
-#
-# -------------------------------------------------------
-#
-# *) The nodes, triangles, boundary edges and communication
-# pattern defined by the new numbering scheme are returned
-#
-#########################################################
-
-def build_local_mesh(submesh, lower_t, upper_t, nproc):
-
-    # Combine the full nodes and ghost nodes
-
-    nodes = num.concatenate((submesh["full_nodes"], \
-                         submesh["ghost_nodes"]))
-    
-    # Combine the full triangles and ghost triangles
-
-    gtri =  num.take(submesh["ghost_triangles"],(1, 2, 3),1)
-    triangles = num.concatenate((submesh["full_triangles"], gtri))
-
-    # Combine the full boundaries and ghost boundaries
-
-    boundaries = submesh["full_boundary"]
-    for b in submesh["ghost_boundary"]:
-        boundaries[b]=submesh["ghost_boundary"][b]
-
-    # Make note of the new triangle numbers, including the ghost
-    # triangles
-
-    NGlobal = upper_t
-    for i in range(len(submesh["ghost_triangles"])):
-        id = submesh["ghost_triangles"][i][0]
-        if id > NGlobal:
-            NGlobal = id
-    index = num.zeros(int(NGlobal)+1, num.int)
-    index[lower_t:upper_t]=num.arange(upper_t-lower_t)
-    for i in range(len(submesh["ghost_triangles"])):
-        index[submesh["ghost_triangles"][i][0]] = i+upper_t-lower_t
-    
-    # Change the node numbering (and update the numbering in the
-    # triangles)
-
-    [GAnodes, GAtriangles, GAboundary] = build_local_GA(nodes, triangles, boundaries, index)
-
-    # Extract the local quantities
-    
-    quantities ={}
-    for k in submesh["full_quan"]:
-        Nf = len(submesh["full_quan"][k])
-        Ng = len(submesh["ghost_quan"][k])
-        quantities[k] = num.zeros((Nf+Ng, 3), num.float)
-        quantities[k][0:Nf] = submesh["full_quan"][k] 
-        quantities[k][Nf:Nf+Ng] = submesh["ghost_quan"][k]
-                             
-    # Change the communication pattern into a form needed by
-    # the parallel_adv
-
-    gcommun = submesh["ghost_commun"]
-    fcommun = submesh["full_commun"]
-    [ghost_rec, full_send] = \
-                build_local_commun(index, gcommun, fcommun, nproc)
-
-    # Clean up before exiting
-
-    del(index)
-
-    return GAnodes, GAtriangles, GAboundary, quantities, ghost_rec, \
-           full_send
-
Index: uga_core/source/anuga_parallel/build_submesh.py
===================================================================
--- /anuga_core/source/anuga_parallel/build_submesh.py	(revision 7446)
+++ 	(revision )
@@ -1,559 +1,0 @@
-#########################################################
-#
-# Subdivide the domain. This module is primarily
-# responsible for building the ghost layer and
-# communication pattern
-#
-#
-#  Author: Linda Stals, June 2005
-#  Modified: Linda Stals, Nov 2005 (optimise python code)
-#            Steve Roberts, Aug 2009 (convert to numpy)
-#
-#
-#########################################################
-
-import sys
-
-#from Numeric import zeros, Float, Int, concatenate, \
-#     reshape, arrayrange, take, nonzero
-
-import numpy as num
-
-from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
-
-
-
-#########################################################
-#
-# Subdivide the triangles into non-overlapping domains.
-#
-#  *)  The subdivision is controlled by triangles_per_proc.
-# The first triangles_per_proc[0] triangles are assigned
-# to the first processor, the second triangles_per_proc[1]
-# are assigned to the second processor etc.
-#
-#  *) nodes, triangles and boundary contains all of the
-# nodes, triangles and boundary tag information for the
-# whole domain. The triangles should be orientated in the
-# correct way and the nodes number consecutively from 0.
-#
-# -------------------------------------------------------
-#
-#  *) A dictionary containing the full_nodes, full_triangles
-# and full_boundary information for each processor is
-# returned. The node information consists of
-# [global_id, x_coord, y_coord].
-#
-#########################################################
-
-def submesh_full(nodes, triangles, boundary, triangles_per_proc):
-
-    # Initialise
-
-    tlower = 0
-    nproc = len(triangles_per_proc)
-    nnodes = len(nodes)
-    node_list = []
-    triangle_list = []
-    boundary_list = []
-    submesh = {}
-    node_range = num.reshape(num.arange(nnodes),(nnodes,1))
-
-    #print node_range
-    tsubnodes = num.concatenate((node_range, nodes), 1)
-
-
-    # Loop over processors
-
-    for p in range(nproc):
-
-        # Find triangles on processor p
-
-        tupper = triangles_per_proc[p]+tlower
-        subtriangles = triangles[tlower:tupper]
-        triangle_list.append(subtriangles)
-
-        # Find the boundary edges on processor p
-
-        subboundary = {}
-        for k in boundary:
-            if (k[0] >=tlower and k[0] < tupper):
-                subboundary[k]=boundary[k]
-        boundary_list.append(subboundary)
-
-        # Find nodes in processor p
-
-        nodemap = num.zeros(nnodes, 'i')
-        for t in subtriangles:
-            nodemap[t[0]]=1
-            nodemap[t[1]]=1
-            nodemap[t[2]]=1
-
-        
-        node_list.append(tsubnodes.take(num.flatnonzero(nodemap),axis=0))
-
-        # Move to the next processor
-
-        tlower = tupper
-
-    # Put the results in a dictionary
-
-    submesh["full_nodes"] = node_list
-    submesh["full_triangles"] = triangle_list
-    submesh["full_boundary"] = boundary_list
-
-    # Clean up before exiting
-
-    del (nodemap)
-
-    return submesh
-
-
-#########################################################
-#
-# Build the ghost layer of triangles
-#
-#  *) Given the triangle subpartion for the processor
-# build a ghost layer of triangles. The ghost layer
-# consists of two layers of neighbouring triangles.
-#
-#  *) The vertices in the ghost triangles must also
-# be added to the node list for the current processor
-#
-#
-# -------------------------------------------------------
-#
-#  *) The extra triangles and nodes are returned.
-#
-#  *)  The node information consists of
-# [global_id, x_coord, y_coord].
-#
-#  *) The triangle information consists of
-# [triangle number, t], where t = [v1, v2, v3].
-#
-#########################################################
-
-def ghost_layer(submesh, mesh, p, tupper, tlower):
-
-    ncoord = mesh.number_of_nodes
-    ntriangles = mesh.number_of_triangles
-
-    # Find the first layer of boundary triangles
-
-    trianglemap = num.zeros(ntriangles, 'i')
-    for t in range(tlower, tupper):
-        
-        n = mesh.neighbours[t, 0]
-
-        if n >= 0:
-            if n < tlower or n >= tupper:
-                trianglemap[n] = 1
-        n = mesh.neighbours[t, 1]
-        if n >= 0:
-            if n < tlower or n >= tupper:
-                trianglemap[n] = 1
-        n = mesh.neighbours[t, 2]
-        if n >= 0:
-            if n < tlower or n >= tupper:
-                trianglemap[n] = 1
-
-    # Find the second layer of boundary triangles
-
-    for t in range(len(trianglemap)):
-        if trianglemap[t]==1:
-            n = mesh.neighbours[t, 0]
-            if n >= 0:
-                if (n < tlower or n >= tupper) and trianglemap[n] == 0:
-                    trianglemap[n] = 2
-            n = mesh.neighbours[t, 1]
-            if n >= 0:
-                if (n < tlower or n >= tupper) and trianglemap[n] == 0:
-                    trianglemap[n] = 2
-            n = mesh.neighbours[t, 2]
-            if n >= 0:
-                if (n < tlower or n >= tupper) and trianglemap[n] == 0:
-                    trianglemap[n] = 2
-
-    # Build the triangle list and make note of the vertices
-
-    nodemap = num.zeros(ncoord, 'i')
-    fullnodes = submesh["full_nodes"][p]
-
-    subtriangles = []
-    for i in range(len(trianglemap)):
-        if trianglemap[i] != 0:
-            t = list(mesh.triangles[i])
-            nodemap[t[0]] = 1
-            nodemap[t[1]] = 1
-            nodemap[t[2]] = 1
-
-    trilist = num.reshape(num.arange(ntriangles),(ntriangles,1))
-    tsubtriangles = num.concatenate((trilist, mesh.triangles), 1)
-    subtriangles = tsubtriangles.take(num.flatnonzero(trianglemap),axis=0)
-
-    
-    # Keep a record of the triangle vertices, if they are not already there
-
-    subnodes = []
-    for n in fullnodes:
-        nodemap[int(n[0])] = 0
-
-    nodelist = num.reshape(num.arange(ncoord),(ncoord,1))
-    tsubnodes = num.concatenate((nodelist, mesh.get_nodes()), 1)
-    subnodes = tsubnodes.take(num.flatnonzero(nodemap),axis=0)
-
-    # Clean up before exiting
-
-    del (nodelist)
-    del (trilist)
-    del (tsubnodes)
-    del (nodemap)
-    del (trianglemap)
-
-    # Return the triangles and vertices sitting on the boundary layer
-
-    return subnodes, subtriangles
-
-#########################################################
-#
-# Find the edges of the ghost trianlges that do not
-# have a neighbour in the current cell. These are
-# treated as a special type of boundary edge.
-#
-#  *) Given the ghost triangles in a particular
-# triangle, use the mesh to find its neigbours. If
-# the neighbour is not in the processor set it to
-# be a boundary edge
-#
-#  *) The vertices in the ghost triangles must also
-# be added to the node list for the current processor
-#
-#  *) The boundary edges for the ghost triangles are
-# ignored.
-#
-# -------------------------------------------------------
-#
-#  *) The type assigned to the ghost boundary edges is 'ghost'
-#
-#  *)  The boundary information is returned as a directorier
-# with the key = (triangle id, edge no) and the values
-# assigned to the key is 'ghost'
-#
-#
-#########################################################
-def is_in_processor(ghost_list, tlower, tupper, n):
-
-    return num.equal(ghost_list,n).any() or (tlower <= n and tupper > n)
-
-
-def ghost_bnd_layer(ghosttri, tlower, tupper, mesh, p):
-
-    ghost_list = []
-    subboundary = {}
-
-
-    for t in ghosttri:
-        ghost_list.append(t[0])
-    
-    for t in ghosttri:
-
-        n = mesh.neighbours[t[0], 0]
-        if not is_in_processor(ghost_list, tlower, tupper, n):
-            subboundary[t[0], 0] = 'ghost'
-
-        n = mesh.neighbours[t[0], 1]
-        if not is_in_processor(ghost_list, tlower, tupper, n):
-            subboundary[t[0], 1] = 'ghost'
-
-        n = mesh.neighbours[t[0], 2]
-        if not is_in_processor(ghost_list, tlower, tupper, n):
-            subboundary[t[0], 2] = 'ghost'
-            
-    return subboundary
-
-#########################################################
-#
-# The ghost triangles on the current processor will need
-# to get updated information from the neighbouring
-# processor containing the corresponding full triangles.
-#
-#  *) The tri_per_proc is used to determine which
-# processor contains the full node copy.
-#
-# -------------------------------------------------------
-#
-#  *) The ghost communication pattern consists of
-# [global node number, neighbour processor number].
-#
-#########################################################
-
-def ghost_commun_pattern(subtri, p, tri_per_proc):
-
-    # Loop over the ghost triangles
-
-    ghost_commun = num.zeros((len(subtri), 2), num.int)
-
-    for i in range(len(subtri)):
-        global_no = subtri[i][0]
-
-        # Find which processor contains the full triangle
-
-        nproc = len(tri_per_proc)
-        neigh = nproc-1
-        sum = 0
-        for q in range(nproc-1):
-            if (global_no < sum+tri_per_proc[q]):
-                neigh = q
-                break
-            sum = sum+tri_per_proc[q]
-
-        # Keep a copy of the neighbour processor number
-
-        ghost_commun[i] = [global_no, neigh]
-
-    return ghost_commun
-
-#########################################################
-#
-# The full triangles in this processor must communicate
-# updated information to neighbouring processor that
-# contain ghost triangles
-#
-#  *) The ghost communication pattern for all of the
-# processor must be built before calling this processor.
-#
-#  *) The full communication pattern is found by looping
-# through the ghost communication pattern for all of the
-# processors. Recall that this information is stored in
-# the form [global node number, neighbour processor number].
-# The full communication for the neighbour processor is
-# then updated.
-#
-# -------------------------------------------------------
-#
-#  *) The full communication pattern consists of
-# [global id, [p1, p2, ...]], where p1, p2 etc contain
-# a ghost node copy of the triangle global id.
-#
-#########################################################
-
-def full_commun_pattern(submesh, tri_per_proc):
-    tlower = 0
-    nproc = len(tri_per_proc)
-    full_commun = []
-
-    # Loop over the processor
-
-    for p in range(nproc):
-
-        # Loop over the full triangles in the current processor
-        # and build an empty dictionary
-
-        fcommun = {}
-        tupper = tri_per_proc[p]+tlower
-        for i in range(tlower, tupper):
-            fcommun[i] = []
-        full_commun.append(fcommun)
-        tlower = tupper
-
-    # Loop over the processor again
-
-    for p in range(nproc):
-
-        # Loop over the ghost triangles in the current processor,
-        # find which processor contains the corresponding full copy
-        # and note that the processor must send updates to this
-        # processor
-
-        for g in submesh["ghost_commun"][p]:
-            neigh = g[1]
-            full_commun[neigh][g[0]].append(p)
-
-    return full_commun
-
-
-#########################################################
-#
-# Given the non-overlapping grid partition, an extra layer
-# of triangles are included to help with the computations.
-# The triangles in this extra layer are not updated by
-# the processor, their updated values must be sent by the
-# processor containing the original, full, copy of the
-# triangle. The communication pattern that controls these
-# updates must also be built.
-#
-#  *) Assumes that full triangles, nodes etc have already
-# been found and stored in submesh
-#
-#  *) See the documentation for ghost_layer,
-# ghost_commun_pattern and full_commun_pattern
-#
-# -------------------------------------------------------
-#
-#  *) The additional information is added to the submesh
-# dictionary. See the documentation for ghost_layer,
-# ghost_commun_pattern and full_commun_pattern
-#
-#  *) The ghost_triangles, ghost_nodes, ghost_boundary,
-# ghost_commun and full_commun is added to submesh
-#########################################################
-
-def submesh_ghost(submesh, mesh, triangles_per_proc):
-
-    nproc = len(triangles_per_proc)
-    tlower = 0
-    ghost_triangles = []
-    ghost_nodes = []
-    ghost_commun = []
-    ghost_bnd = []
-
-    # Loop over the processors
-
-    for p in range(nproc):
-
-        # Find the full triangles in this processor
-
-        tupper = triangles_per_proc[p]+tlower
-
-        # Build the ghost boundary layer
-
-        [subnodes, subtri] = \
-                   ghost_layer(submesh, mesh, p, tupper, tlower)
-        ghost_triangles.append(subtri)
-        ghost_nodes.append(subnodes)
-
-
-        # Find the boundary layer formed by the ghost triangles
-        
-        subbnd = ghost_bnd_layer(subtri, tlower, tupper, mesh, p)
-        ghost_bnd.append(subbnd)
-        
-        # Build the communication pattern for the ghost nodes
-
-        gcommun = \
-                ghost_commun_pattern(subtri, p, triangles_per_proc)
-        ghost_commun.append(gcommun)
-
-        # Move to the next processor
-
-        tlower = tupper
-
-
-    # Record the ghost layer and communication pattern
-
-    submesh["ghost_nodes"] = ghost_nodes
-    submesh["ghost_triangles"] = ghost_triangles
-    submesh["ghost_commun"] = ghost_commun
-    submesh["ghost_boundary"] = ghost_bnd
-    
-    # Build the communication pattern for the full triangles
-
-    full_commun = full_commun_pattern(submesh, triangles_per_proc)
-    submesh["full_commun"] = full_commun
-
-    # Return the submesh
-
-    return submesh
-
-
-#########################################################
-#
-# Certain quantities may be assigned to the triangles,
-# these quantities must be subdivided in the same way
-# as the triangles
-#
-#  *) The quantities are ordered in the same way as the
-# triangles
-#
-# -------------------------------------------------------
-#
-#  *) The quantites attached to the full triangles are
-# stored in full_quan
-#
-#  *) The quantities attached to the ghost triangles are
-# stored in ghost_quan
-#########################################################
-
-def submesh_quantities(submesh, quantities, triangles_per_proc):
-
-    nproc = len(triangles_per_proc)
-
-    lower = 0
-
-    # Build an empty dictionary to hold the quantites
-
-    submesh["full_quan"] = {}
-    submesh["ghost_quan"] = {}
-    for k in quantities:
-        submesh["full_quan"][k] = []
-        submesh["ghost_quan"][k] = []
-
-    # Loop trough the subdomains
-
-    for p in range(nproc):
-        upper =   lower+triangles_per_proc[p]
-
-        # Find the global ID of the ghost triangles
-
-        global_id = []
-        M = len(submesh["ghost_triangles"][p])
-        for j in range(M):
-            global_id.append(submesh["ghost_triangles"][p][j][0])
-
-        # Use the global ID to extract the quantites information from
-        # the full domain
-
-        for k in quantities:
-            submesh["full_quan"][k].append(quantities[k][lower:upper])
-            submesh["ghost_quan"][k].append(num.zeros( (M,3) , num.float))
-            for j in range(M):
-                submesh["ghost_quan"][k][p][j] = \
-                                               quantities[k][global_id[j]]
-
-        lower = upper
-
-    return submesh
-
-#########################################################
-#
-# Build the grid partition on the host.
-#
-#  *) See the documentation for submesh_ghost and
-# submesh_full
-#
-# -------------------------------------------------------
-#
-#  *) A dictionary containing the full_triangles,
-# full_nodes, full_boundary, ghost_triangles, ghost_nodes,
-# ghost_boundary, ghost_commun and full_commun and true boundary polygon is returned.
-#
-#########################################################
-
-def build_submesh(nodes, triangles, edges, quantities,
-                  triangles_per_proc):
-
-    # Temporarily build the mesh to find the neighbouring
-    # triangles and true boundary polygon
-
-    mesh = Mesh(nodes, triangles)
-    boundary_polygon = mesh.get_boundary_polygon()
-    
-
-    # Subdivide into non-overlapping partitions
-
-    submeshf = submesh_full(nodes, triangles, edges, \
-                            triangles_per_proc)
-    
-    # Add any extra ghost boundary layer information
-
-    submeshg = submesh_ghost(submeshf, mesh, triangles_per_proc)
-
-    # Order the quantities information to be the same as the triangle
-    # information
-
-    submesh = submesh_quantities(submeshg, quantities, \
-                                 triangles_per_proc)
-
-    submesh["boundary_polygon"] = boundary_polygon
-    return submesh
-
Index: /anuga_core/source/anuga_parallel/parallel_advection.py
===================================================================
--- /anuga_core/source/anuga_parallel/parallel_advection.py	(revision 7446)
+++ /anuga_core/source/anuga_parallel/parallel_advection.py	(revision 7447)
@@ -90,5 +90,5 @@
         #print self.processor, ltimestep, gtimestep
         
-        pypar.reduce(ltimestep, pypar.MIN, 0, buffer=gtimestep)
+        gtimestep = pypar.reduce(ltimestep, pypar.MIN, 0, buffer=gtimestep)
 
         #print self.processor, ltimestep, gtimestep
@@ -182,32 +182,33 @@
 
 
-    def write_time(self):
-        if self.min_timestep == self.max_timestep:
-            print 'Processor %d, Time = %.4f, delta t = %.8f, steps=%d (%d)'\
-                  %(self.processor, self.time, self.min_timestep, self.number_of_steps,
-                    self.number_of_first_order_steps)
-        elif self.min_timestep > self.max_timestep:
-            print 'Processor %d, Time = %.4f, steps=%d (%d)'\
-                  %(self.processor, self.time, self.number_of_steps,
-                    self.number_of_first_order_steps)
-        else:
-            print 'Processor %d, Time = %.4f, delta t in [%.8f, %.8f], steps=%d (%d)'\
-                  %(self.processor, self.time, self.min_timestep,
-                    self.max_timestep, self.number_of_steps,
-                    self.number_of_first_order_steps)
-
-
-
-    def evolve(self, yieldstep = None, finaltime = None):
-        """Specialisation of basic evolve method from parent class
-        """
-
-        #Initialise real time viz if requested
-        if self.time == 0.0:
-            pass
-
-        #Call basic machinery from parent class
-        for t in Domain.evolve(self, yieldstep, finaltime):
-
-            #Pass control on to outer loop for more specific actions
-            yield(t)
+    ## def write_time(self):
+    ##     if self.min_timestep == self.max_timestep:
+    ##         print 'Processor %d, Time = %.4f, delta t = %.8f, steps=%d (%d)'\
+    ##               %(self.processor, self.time, self.min_timestep, self.number_of_steps,
+    ##                 self.number_of_first_order_steps)
+    ##     elif self.min_timestep > self.max_timestep:
+    ##         print 'Processor %d, Time = %.4f, steps=%d (%d)'\
+    ##               %(self.processor, self.time, self.number_of_steps,
+    ##                 self.number_of_first_order_steps)
+    ##     else:
+    ##         print 'Processor %d, Time = %.4f, delta t in [%.8f, %.8f], steps=%d (%d)'\
+    ##               %(self.processor, self.time, self.min_timestep,
+    ##                 self.max_timestep, self.number_of_steps,
+    ##                 self.number_of_first_order_steps)
+
+
+
+    ## def evolve(self, yieldstep = None, finaltime = None):
+    ##     """Specialisation of basic evolve method from parent class
+    ##     """
+
+    ##     #Initialise real time viz if requested
+    ##     if self.time == 0.0:
+    ##         pass
+
+    ##     #Call basic machinery from parent class
+    ##     for t in Domain.evolve(self, yieldstep, finaltime):
+
+    ##         #Pass control on to outer loop for more specific actions
+    ##         yield(t)
+
Index: /anuga_core/source/anuga_parallel/parallel_meshes.py
===================================================================
--- /anuga_core/source/anuga_parallel/parallel_meshes.py	(revision 7446)
+++ /anuga_core/source/anuga_parallel/parallel_meshes.py	(revision 7447)
@@ -43,6 +43,6 @@
     numproc   = pypar.size()
 
-    print 'numproc',numproc
-    print 'processor ',processor
+    #print 'numproc',numproc
+    #print 'processor ',processor
 
     m_low, m_high = pypar.balance(m_g, numproc, processor)
@@ -52,5 +52,5 @@
     m_high = m_high+1
 
-    print 'm_low, m_high', m_low, m_high
+    #print 'm_low, m_high', m_low, m_high
     m = m_high - m_low
 
@@ -168,18 +168,18 @@
         Idgr.extend(Idgl)
 
-        print Idfl
-        print Idgr
+        #print Idfl
+        #print Idgr
         
         Idfl = num.array(Idfl,num.int)
         Idgr = num.array(Idgr,num.int)
 
-        print Idfl
-        print Idgr
+        #print Idfl
+        #print Idgr
         
         full_send_dict[processor]  = [Idfl, Idfl]
         ghost_recv_dict[processor] = [Idgr, Idgr]
 
-        print  full_send_dict[processor]
-        print ghost_recv_dict[processor]
+        #print  full_send_dict[processor]
+        #print ghost_recv_dict[processor]
     elif numproc == 2:
         Idfl.extend(Idfr)
@@ -203,6 +203,6 @@
 
 
-    print full_send_dict
-    print ghost_recv_dict        
+    #print full_send_dict
+    #print ghost_recv_dict        
     
     return  points, elements, boundary, full_send_dict, ghost_recv_dict
Index: /anuga_core/source/anuga_parallel/parallel_shallow_water.py
===================================================================
--- /anuga_core/source/anuga_parallel/parallel_shallow_water.py	(revision 7446)
+++ /anuga_core/source/anuga_parallel/parallel_shallow_water.py	(revision 7447)
@@ -22,5 +22,5 @@
 
 
-from anuga.shallow_water.shallow_water_domain import *
+from anuga.interface import Domain
 
 
@@ -52,20 +52,4 @@
         N = len(self) # number_of_triangles
 
-#        self.processor = pypar.rank()
-#        self.numproc   = pypar.size()
-#
-#        # Setup Communication Buffers
-#        self.nsys = 3
-#        for key in full_send_dict:
-#            buffer_shape = full_send_dict[key][0].shape[0]
-#            full_send_dict[key].append(num.zeros( (buffer_shape,self.nsys) ,num.loat))
-#
-#
-#        for key in ghost_recv_dict:
-#            buffer_shape = ghost_recv_dict[key][0].shape[0]
-#            ghost_recv_dict[key].append(num.zeros( (buffer_shape,self.nsys) ,num.float))
-#
-#        self.full_send_dict  = full_send_dict
-        self.ghost_recv_dict = ghost_recv_dict
 
         # Buffers for synchronisation of timesteps
@@ -73,5 +57,5 @@
         self.global_timestep = num.zeros(1, num.float)
 
-        self.local_timesteps = num.zeros(self.numproc, num.loat)
+        self.local_timesteps = num.zeros(self.numproc, num.float)
 
 
@@ -150,4 +134,5 @@
 
         import time
+
         #Compute minimal timestep across all processes
         self.local_timestep[0] = self.flux_timestep
@@ -167,6 +152,5 @@
                 for i in range(1, self.numproc):
                     pypar.receive(i,
-                                  buffer=self.local_timestep,
-                                  bypass=True)
+                                  buffer=self.local_timestep)
 
                     if self.local_timestep[0] < self.global_timestep[0]:
@@ -174,5 +158,5 @@
             else:
                 pypar.send(self.local_timestep, 0,
-                           use_buffer=True, bypass=True)
+                           use_buffer=True)
 
 
@@ -187,5 +171,5 @@
         self.communication_broadcast_time += time.time()-t0
 
-        #old_timestep = self.flux_timestep
+        old_timestep = self.flux_timestep
         self.flux_timestep = self.global_timestep[0]
         #print 'Flux Timestep %15.5e %15.5e P%d_%d' %(self.flux_timestep, old_timestep, self.processor, self.numproc)
@@ -209,6 +193,5 @@
         # the separate processors
 
-
-        from Numeric import take,put
+        import numpy as num
         import time
         t0 = time.time()
@@ -227,8 +210,7 @@
                             #print 'Send',i,q
                             Q_cv =  self.quantities[q].centroid_values
-                            Xout[:,i] = take(Q_cv, Idf)
-
-                        pypar.send(Xout, send_proc,
-                                   use_buffer=True, bypass = True)
+                            Xout[:,i] = num.take(Q_cv, Idf)
+
+                        pypar.send(Xout, int(send_proc), use_buffer=True)
 
 
@@ -238,12 +220,12 @@
 
                     Idg = self.ghost_recv_dict[iproc][0]
-                    X = self.ghost_recv_dict[iproc][2]
-
-                    X = pypar.receive(iproc, buffer=X, bypass = True)
+                    X   = self.ghost_recv_dict[iproc][2]
+
+                    X = pypar.receive(int(iproc), buffer=X)
 
                     for i, q in enumerate(self.conserved_quantities):
                         #print 'Receive',i,q
                         Q_cv =  self.quantities[q].centroid_values
-                        put(Q_cv, Idg, X[:,i])
+                        num.put(Q_cv, Idg, X[:,i])
 
         #local update of ghost cells
@@ -262,30 +244,30 @@
                 #print 'LOCAL SEND RECEIVE',i,q
                 Q_cv =  self.quantities[q].centroid_values
-                put(Q_cv,     Idg, take(Q_cv,     Idf))
+                num.put(Q_cv, Idg, num.take(Q_cv, Idf))
 
         self.communication_time += time.time()-t0
 
-'''
-This was removed due to not beening required to be redefined in parallel_shallow_water
-the original "write_time" is good... however might need some small edits to work properly 
-with parallel- Nick and Ole April 2007
-    def write_time(self):
-        if self.min_timestep == self.max_timestep:
-            print 'Processor %d/%d, Time = %.4f, delta t = %.8f, steps=%d (%d)'\
-                  %(self.processor, self.numproc,
-                    self.time, self.min_timestep, self.number_of_steps,
-                    self.number_of_first_order_steps)
-        elif self.min_timestep > self.max_timestep:
-            print 'Processor %d/%d, Time = %.4f, steps=%d (%d)'\
-                  %(self.processor, self.numproc,
-                    self.time, self.number_of_steps,
-                    self.number_of_first_order_steps)
-        else:
-            print 'Processor %d/%d, Time = %.4f, delta t in [%.8f, %.8f], steps=%d (%d)'\
-                  %(self.processor, self.numproc,
-                    self.time, self.min_timestep,
-                    self.max_timestep, self.number_of_steps,
-                    self.number_of_first_order_steps)
-'''
+
+## This was removed due to not beening required to be redefined in parallel_shallow_water
+## the original "write_time" is good... however might need some small edits to work properly 
+## with parallel- Nick and Ole April 2007
+##     def write_time(self):
+##         if self.min_timestep == self.max_timestep:
+##             print 'Processor %d/%d, Time = %.4f, delta t = %.8f, steps=%d (%d)'\
+##                   %(self.processor, self.numproc,
+##                     self.time, self.min_timestep, self.number_of_steps,
+##                     self.number_of_first_order_steps)
+##         elif self.min_timestep > self.max_timestep:
+##             print 'Processor %d/%d, Time = %.4f, steps=%d (%d)'\
+##                   %(self.processor, self.numproc,
+##                     self.time, self.number_of_steps,
+##                     self.number_of_first_order_steps)
+##         else:
+##             print 'Processor %d/%d, Time = %.4f, delta t in [%.8f, %.8f], steps=%d (%d)'\
+##                   %(self.processor, self.numproc,
+##                     self.time, self.min_timestep,
+##                     self.max_timestep, self.number_of_steps,
+##                     self.number_of_first_order_steps)
+
 
 # commented out on the 7/11/06
Index: /anuga_core/source/anuga_parallel/run_parallel_advection.py
===================================================================
--- /anuga_core/source/anuga_parallel/run_parallel_advection.py	(revision 7446)
+++ /anuga_core/source/anuga_parallel/run_parallel_advection.py	(revision 7447)
@@ -140,5 +140,5 @@
     t0 = time.time()
 
-for t in domain.evolve(yieldstep = 0.1, finaltime = 30.0):
+for t in domain.evolve(yieldstep = 5.0, finaltime = 30.0):
     if myid == 0:
         domain.write_time()
Index: /anuga_core/source/anuga_parallel/run_parallel_sw_rectangle.py
===================================================================
--- /anuga_core/source/anuga_parallel/run_parallel_sw_rectangle.py	(revision 7446)
+++ /anuga_core/source/anuga_parallel/run_parallel_sw_rectangle.py	(revision 7447)
@@ -1,3 +1,2 @@
-#!/usr/bin/env python
 #########################################################
 #
@@ -18,6 +17,7 @@
 import time
 
-from Numeric import array
+#from Numeric import array
 # pmesh
+import numpy as num
 
 from print_stats import print_test_stats, build_full_flag
@@ -35,5 +35,5 @@
 numprocs = pypar.size()
 myid = pypar.rank()
-processor_name = pypar.Get_processor_name()
+processor_name = pypar.get_processor_name()
 
 M = 50
@@ -126,5 +126,5 @@
 
 
-yieldstep = 0.01
+yieldstep = 0.1
 finaltime = 1.0
 
Index: uga_core/source/anuga_parallel/test_distribute_mesh.py
===================================================================
--- /anuga_core/source/anuga_parallel/test_distribute_mesh.py	(revision 7446)
+++ 	(revision )
@@ -1,755 +1,0 @@
-#!/usr/bin/env python
-
-import unittest
-import sys
-from math import sqrt
-
-
-from anuga.interface import Domain
-from anuga.interface import rectangular_cross
-
-from parallel_shallow_water import Parallel_Domain
-
-from anuga_parallel.pmesh_divide import pmesh_divide_metis
-from anuga_parallel.build_submesh import build_submesh
-from anuga_parallel.build_submesh import submesh_full, submesh_ghost, submesh_quantities
-from anuga_parallel.build_commun import extract_hostmesh, rec_submesh, send_submesh
-
-#import numpy as num
-import Numeric as num
-
-def topography(x,y): 
-    return -x/2
-
-
-def xcoord(x,y):
-    return x
-
-def ycoord(x,y):
-    return y
-
-
-
-class Test_Distribute_Mesh(unittest.TestCase):
-    def setUp(self):
-        pass
-
-
-    def tearDown(self):
-        pass
-
-
-
-    def test_pmesh_1(self):
-
-        points, vertices, boundary = rectangular_cross(2,2)
-
-
-        true_points = [[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]]
-
-        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]]
-
-
-        assert num.allclose(points,true_points)
-        assert num.allclose(vertices,true_vertices)
-
-        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) #
-
-
-        #print domain.quantities['ymomentum'].centroid_values
- 
-        nodes, triangles, boundary, triangles_per_proc, quantities = pmesh_divide_metis(domain,1)
-
-
-        true_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]]
-
-        true_triangles = [[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]]
-
-
-        assert num.allclose(nodes,true_nodes)
-        assert num.allclose(triangles,true_triangles)
-
-
-        assert num.allclose(triangles_per_proc,[16])
-        
-        
-
-    def test_pmesh_2(self):
-
-        points, vertices, boundary = rectangular_cross(2,2)
-
-
-        true_points = [[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]]
-
-        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]]
-
-
-        assert num.allclose(points,true_points)
-        assert num.allclose(vertices,true_vertices)
-
-        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) #
-
-
-        #print domain.quantities['ymomentum'].centroid_values
- 
-        nodes, triangles, boundary, triangles_per_proc, quantities = pmesh_divide_metis(domain,2)
-
-
-        true_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]]
-
-
-        true_triangles = [[0, 9, 1], [3, 9, 0], [4, 9, 3], [1, 9, 4], [4, 10, 1], [3, 11, 4], [4, 11, 7], [4, 12, 5], [1, 10, 2], [5, 10, 4], [2, 10, 5], [6, 11, 3], [7, 11, 6], [7, 12, 4], [8, 12, 7], [5, 12, 8]]
-
-
-        assert num.allclose(nodes,true_nodes)
-        assert num.allclose(triangles,true_triangles)
-
-
-        assert num.allclose(triangles_per_proc,[8,8])
-
-
-    def test_distibute_3(self):
-
-        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)
-
-            #----------------------------------------------------------------------------------
-            # 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)    
-
-            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]) 
-
-                """
-        par_domain = Parallel_Domain(points, vertices, boundary,
-                                     full_send_dict  = full_send_dict,
-                                     ghost_recv_dict = ghost_recv_dict)
-                                     """
-
-
-
-
-
-
-            
-    def test_build_submesh_3(self):
-
-
-        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]]
-
-
-        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]]
-
-
-        edges = {(13, 1): 'bottom', (7, 1): 'left', (3, 1): 'right', (14, 1): 'right', (11, 1): 'bottom', (10, 1): 'top', (5, 1): 'left', (4, 1): 'top'}
-
-        triangles_per_proc = [5, 6, 5]
-
-
-        quantities = {'stage': num.array([[-0.25 , -0.125, -0.25 ],
-       [-0.25 , -0.375, -0.25 ],
-       [-0.5  , -0.375, -0.25 ],
-       [-0.5  , -0.375, -0.5  ],
-       [-0.25 , -0.375, -0.5  ],
-       [-0.   , -0.125, -0.   ],
-       [-0.   , -0.125, -0.25 ],
-       [-0.   , -0.125, -0.   ],
-       [-0.25 , -0.125, -0.   ],
-       [-0.25 , -0.125, -0.25 ],
-       [-0.   , -0.125, -0.25 ],
-       [-0.25 , -0.125, -0.   ],
-       [-0.25 , -0.375, -0.25 ],
-       [-0.5  , -0.375, -0.25 ],
-       [-0.5  , -0.375, -0.5  ],
-       [-0.25 , -0.375, -0.5  ]]),  'elevation': num.array([[-0.25 , -0.125, -0.25 ],
-       [-0.25 , -0.375, -0.25 ],
-       [-0.5  , -0.375, -0.25 ],
-       [-0.5  , -0.375, -0.5  ],
-       [-0.25 , -0.375, -0.5  ],
-       [-0.   , -0.125, -0.   ],
-       [-0.   , -0.125, -0.25 ],
-       [-0.   , -0.125, -0.   ],
-       [-0.25 , -0.125, -0.   ],
-       [-0.25 , -0.125, -0.25 ],
-       [-0.   , -0.125, -0.25 ],
-       [-0.25 , -0.125, -0.   ],
-       [-0.25 , -0.375, -0.25 ],
-       [-0.5  , -0.375, -0.25 ],
-       [-0.5  , -0.375, -0.5  ],
-       [-0.25 , -0.375, -0.5  ]]),  'ymomentum': num.array([[ 0.5 ,  0.25,  0.  ],
-       [ 0.5 ,  0.75,  1.  ],
-       [ 0.5 ,  0.75,  0.5 ],
-       [ 1.  ,  0.75,  0.5 ],
-       [ 1.  ,  0.75,  1.  ],
-       [ 0.  ,  0.25,  0.5 ],
-       [ 0.5 ,  0.25,  0.5 ],
-       [ 0.5 ,  0.75,  1.  ],
-       [ 0.5 ,  0.75,  0.5 ],
-       [ 1.  ,  0.75,  0.5 ],
-       [ 1.  ,  0.75,  1.  ],
-       [ 0.  ,  0.25,  0.  ],
-       [ 0.  ,  0.25,  0.5 ],
-       [ 0.  ,  0.25,  0.  ],
-       [ 0.5 ,  0.25,  0.  ],
-       [ 0.5 ,  0.25,  0.5 ]]),  'friction': num.array([[ 0.,  0.,  0.],
-       [ 0.,  0.,  0.],
-       [ 0.,  0.,  0.],
-       [ 0.,  0.,  0.],
-       [ 0.,  0.,  0.],
-       [ 0.,  0.,  0.],
-       [ 0.,  0.,  0.],
-       [ 0.,  0.,  0.],
-       [ 0.,  0.,  0.],
-       [ 0.,  0.,  0.],
-       [ 0.,  0.,  0.],
-       [ 0.,  0.,  0.],
-       [ 0.,  0.,  0.],
-       [ 0.,  0.,  0.],
-       [ 0.,  0.,  0.],
-       [ 0.,  0.,  0.]]),  'xmomentum': num.array([[ 2.,  2.,  2.],
-       [ 2.,  2.,  2.],
-       [ 2.,  2.,  2.],
-       [ 2.,  2.,  2.],
-       [ 2.,  2.,  2.],
-       [ 2.,  2.,  2.],
-       [ 2.,  2.,  2.],
-       [ 2.,  2.,  2.],
-       [ 2.,  2.,  2.],
-       [ 2.,  2.,  2.],
-       [ 2.,  2.,  2.],
-       [ 2.,  2.,  2.],
-       [ 2.,  2.,  2.],
-       [ 2.,  2.,  2.],
-       [ 2.,  2.,  2.],
-       [ 2.,  2.,  2.]])}
-
-
-        
-        #----------------------------------------------------------------------------------
-        # Test build_submesh
-        #----------------------------------------------------------------------------------
-        #submesh = build_submesh(nodes, triangles, edges, quantities, triangles_per_proc)
-
-        # Temporarily build the mesh to find the neighbouring
-        # triangles and true boundary polygon
-
-        true_submesh = {'full_boundary': [{(3, 1): 'right', (4, 1): 'top'}, {(5, 1): 'left', (10, 1): 'top', (7, 1): 'left'}, {(13, 1): 'bottom', (14, 1): 'right', (11, 1): 'bottom'}],
-                        'ghost_nodes': [num.array([[  0.  ,   0.  ,   0.  ],
-       [  1.  ,   0.  ,   0.5 ],
-       [  2.  ,   0.  ,   1.  ],
-       [  6.  ,   1.  ,   0.  ],
-       [ 10.  ,   0.25,   0.75],
-       [ 11.  ,   0.75,   0.25]]), num.array([[  3.  ,   0.5 ,   0.  ],
-       [  7.  ,   1.  ,   0.5 ],
-       [  8.  ,   1.  ,   1.  ],
-       [ 11.  ,   0.75,   0.25],
-       [ 12.  ,   0.75,   0.75]]), num.array([[  1.  ,   0.  ,   0.5 ],
-       [  5.  ,   0.5 ,   1.  ],
-       [  8.  ,   1.  ,   1.  ],
-       [ 12.  ,   0.75,   0.75]])],
-                        'full_nodes': [num.array([[  3.  ,   0.5 ,   0.  ],
-       [  4.  ,   0.5 ,   0.5 ],
-       [  5.  ,   0.5 ,   1.  ],
-       [  7.  ,   1.  ,   0.5 ],
-       [  8.  ,   1.  ,   1.  ],
-       [  9.  ,   0.25,   0.25],
-       [ 12.  ,   0.75,   0.75]]), num.array([[  0.  ,   0.  ,   0.  ],
-       [  1.  ,   0.  ,   0.5 ],
-       [  2.  ,   0.  ,   1.  ],
-       [  4.  ,   0.5 ,   0.5 ],
-       [  5.  ,   0.5 ,   1.  ],
-       [  9.  ,   0.25,   0.25],
-       [ 10.  ,   0.25,   0.75]]), num.array([[  0.  ,   0.  ,   0.  ],
-       [  3.  ,   0.5 ,   0.  ],
-       [  4.  ,   0.5 ,   0.5 ],
-       [  6.  ,   1.  ,   0.  ],
-       [  7.  ,   1.  ,   0.5 ],
-       [  9.  ,   0.25,   0.25],
-       [ 11.  ,   0.75,   0.25]])],
-                        '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]])],
-                        'ghost_boundary': [{(13, 1): 'ghost', (8, 0): 'ghost', (14, 1): 'ghost', (11, 1): 'ghost', (10, 1): 'ghost', (5, 1): 'ghost', (10, 2): 'ghost'}, {(12, 2): 'ghost', (12, 0): 'ghost', (2, 1): 'ghost', (11, 1): 'ghost', (2, 2): 'ghost', (4, 1): 'ghost', (4, 0): 'ghost'}, {(3, 2): 'ghost', (6, 1): 'ghost', (3, 1): 'ghost', (5, 1): 'ghost', (1, 0): 'ghost', (1, 1): 'ghost'}],
-                        'full_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]]],
-                        '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]}],
-                        '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]])],  'ghost_quan': {'stage': [num.array([[-0.   , -0.125, -0.   ],
-       [-0.   , -0.125, -0.25 ],
-       [-0.25 , -0.125, -0.   ],
-       [-0.25 , -0.125, -0.25 ],
-       [-0.   , -0.125, -0.25 ],
-       [-0.25 , -0.125, -0.   ],
-       [-0.25 , -0.375, -0.25 ],
-       [-0.5  , -0.375, -0.25 ],
-       [-0.5  , -0.375, -0.5  ],
-       [-0.25 , -0.375, -0.5  ]]), num.array([[-0.25 , -0.125, -0.25 ],
-       [-0.25 , -0.375, -0.25 ],
-       [-0.5  , -0.375, -0.25 ],
-       [-0.25 , -0.375, -0.5  ],
-       [-0.25 , -0.125, -0.   ],
-       [-0.25 , -0.375, -0.25 ]]), num.array([[-0.25 , -0.125, -0.25 ],
-       [-0.25 , -0.375, -0.25 ],
-       [-0.5  , -0.375, -0.25 ],
-       [-0.5  , -0.375, -0.5  ],
-       [-0.   , -0.125, -0.   ],
-       [-0.   , -0.125, -0.25 ]])],  'elevation': [num.array([[-0.   , -0.125, -0.   ],
-       [-0.   , -0.125, -0.25 ],
-       [-0.25 , -0.125, -0.   ],
-       [-0.25 , -0.125, -0.25 ],
-       [-0.   , -0.125, -0.25 ],
-       [-0.25 , -0.125, -0.   ],
-       [-0.25 , -0.375, -0.25 ],
-       [-0.5  , -0.375, -0.25 ],
-       [-0.5  , -0.375, -0.5  ],
-       [-0.25 , -0.375, -0.5  ]]), num.array([[-0.25 , -0.125, -0.25 ],
-       [-0.25 , -0.375, -0.25 ],
-       [-0.5  , -0.375, -0.25 ],
-       [-0.25 , -0.375, -0.5  ],
-       [-0.25 , -0.125, -0.   ],
-       [-0.25 , -0.375, -0.25 ]]), num.array([[-0.25 , -0.125, -0.25 ],
-       [-0.25 , -0.375, -0.25 ],
-       [-0.5  , -0.375, -0.25 ],
-       [-0.5  , -0.375, -0.5  ],
-       [-0.   , -0.125, -0.   ],
-       [-0.   , -0.125, -0.25 ]])],  'ymomentum': [num.array([[ 0.  ,  0.25,  0.5 ],
-       [ 0.5 ,  0.25,  0.5 ],
-       [ 0.5 ,  0.75,  0.5 ],
-       [ 1.  ,  0.75,  0.5 ],
-       [ 1.  ,  0.75,  1.  ],
-       [ 0.  ,  0.25,  0.  ],
-       [ 0.  ,  0.25,  0.5 ],
-       [ 0.  ,  0.25,  0.  ],
-       [ 0.5 ,  0.25,  0.  ],
-       [ 0.5 ,  0.25,  0.5 ]]), num.array([[ 0.5 ,  0.25,  0.  ],
-       [ 0.5 ,  0.75,  1.  ],
-       [ 0.5 ,  0.75,  0.5 ],
-       [ 1.  ,  0.75,  1.  ],
-       [ 0.  ,  0.25,  0.  ],
-       [ 0.  ,  0.25,  0.5 ]]), num.array([[ 0.5 ,  0.25,  0.  ],
-       [ 0.5 ,  0.75,  1.  ],
-       [ 0.5 ,  0.75,  0.5 ],
-       [ 1.  ,  0.75,  0.5 ],
-       [ 0.  ,  0.25,  0.5 ],
-       [ 0.5 ,  0.25,  0.5 ]])],  'friction': [num.array([[ 0.,  0.,  0.],
-       [ 0.,  0.,  0.],
-       [ 0.,  0.,  0.],
-       [ 0.,  0.,  0.],
-       [ 0.,  0.,  0.],
-       [ 0.,  0.,  0.],
-       [ 0.,  0.,  0.],
-       [ 0.,  0.,  0.],
-       [ 0.,  0.,  0.],
-       [ 0.,  0.,  0.]]), num.array([[ 0.,  0.,  0.],
-       [ 0.,  0.,  0.],
-       [ 0.,  0.,  0.],
-       [ 0.,  0.,  0.],
-       [ 0.,  0.,  0.],
-       [ 0.,  0.,  0.]]), num.array([[ 0.,  0.,  0.],
-       [ 0.,  0.,  0.],
-       [ 0.,  0.,  0.],
-       [ 0.,  0.,  0.],
-       [ 0.,  0.,  0.],
-       [ 0.,  0.,  0.]])], 'xmomentum': [num.array([[ 2.,  2.,  2.],
-       [ 2.,  2.,  2.],
-       [ 2.,  2.,  2.],
-       [ 2.,  2.,  2.],
-       [ 2.,  2.,  2.],
-       [ 2.,  2.,  2.],
-       [ 2.,  2.,  2.],
-       [ 2.,  2.,  2.],
-       [ 2.,  2.,  2.],
-       [ 2.,  2.,  2.]]), num.array([[ 2.,  2.,  2.],
-       [ 2.,  2.,  2.],
-       [ 2.,  2.,  2.],
-       [ 2.,  2.,  2.],
-       [ 2.,  2.,  2.],
-       [ 2.,  2.,  2.]]), num.array([[ 2.,  2.,  2.],
-       [ 2.,  2.,  2.],
-       [ 2.,  2.,  2.],
-       [ 2.,  2.,  2.],
-       [ 2.,  2.,  2.],
-       [ 2.,  2.,  2.]])]},  'full_quan': {'stage': [num.array([[-0.25 , -0.125, -0.25 ],
-       [-0.25 , -0.375, -0.25 ],
-       [-0.5  , -0.375, -0.25 ],
-       [-0.5  , -0.375, -0.5  ],
-       [-0.25 , -0.375, -0.5  ]]), num.array([[-0.   , -0.125, -0.   ],
-       [-0.   , -0.125, -0.25 ],
-       [-0.   , -0.125, -0.   ],
-       [-0.25 , -0.125, -0.   ],
-       [-0.25 , -0.125, -0.25 ],
-       [-0.   , -0.125, -0.25 ]]), num.array([[-0.25 , -0.125, -0.   ],
-       [-0.25 , -0.375, -0.25 ],
-       [-0.5  , -0.375, -0.25 ],
-       [-0.5  , -0.375, -0.5  ],
-       [-0.25 , -0.375, -0.5  ]])],  'elevation': [num.array([[-0.25 , -0.125, -0.25 ],
-       [-0.25 , -0.375, -0.25 ],
-       [-0.5  , -0.375, -0.25 ],
-       [-0.5  , -0.375, -0.5  ],
-       [-0.25 , -0.375, -0.5  ]]), num.array([[-0.   , -0.125, -0.   ],
-       [-0.   , -0.125, -0.25 ],
-       [-0.   , -0.125, -0.   ],
-       [-0.25 , -0.125, -0.   ],
-       [-0.25 , -0.125, -0.25 ],
-       [-0.   , -0.125, -0.25 ]]), num.array([[-0.25 , -0.125, -0.   ],
-       [-0.25 , -0.375, -0.25 ],
-       [-0.5  , -0.375, -0.25 ],
-       [-0.5  , -0.375, -0.5  ],
-       [-0.25 , -0.375, -0.5  ]])],  'ymomentum': [num.array([[ 0.5 ,  0.25,  0.  ],
-       [ 0.5 ,  0.75,  1.  ],
-       [ 0.5 ,  0.75,  0.5 ],
-       [ 1.  ,  0.75,  0.5 ],
-       [ 1.  ,  0.75,  1.  ]]), num.array([[ 0.  ,  0.25,  0.5 ],
-       [ 0.5 ,  0.25,  0.5 ],
-       [ 0.5 ,  0.75,  1.  ],
-       [ 0.5 ,  0.75,  0.5 ],
-       [ 1.  ,  0.75,  0.5 ],
-       [ 1.  ,  0.75,  1.  ]]), num.array([[ 0.  ,  0.25,  0.  ],
-       [ 0.  ,  0.25,  0.5 ],
-       [ 0.  ,  0.25,  0.  ],
-       [ 0.5 ,  0.25,  0.  ],
-       [ 0.5 ,  0.25,  0.5 ]])],  'friction': [num.array([[ 0.,  0.,  0.],
-       [ 0.,  0.,  0.],
-       [ 0.,  0.,  0.],
-       [ 0.,  0.,  0.],
-       [ 0.,  0.,  0.]]), num.array([[ 0.,  0.,  0.],
-       [ 0.,  0.,  0.],
-       [ 0.,  0.,  0.],
-       [ 0.,  0.,  0.],
-       [ 0.,  0.,  0.],
-       [ 0.,  0.,  0.]]), num.array([[ 0.,  0.,  0.],
-       [ 0.,  0.,  0.],
-       [ 0.,  0.,  0.],
-       [ 0.,  0.,  0.],
-       [ 0.,  0.,  0.]])],  'xmomentum': [num.array([[ 2.,  2.,  2.],
-       [ 2.,  2.,  2.],
-       [ 2.,  2.,  2.],
-       [ 2.,  2.,  2.],
-       [ 2.,  2.,  2.]]), num.array([[ 2.,  2.,  2.],
-       [ 2.,  2.,  2.],
-       [ 2.,  2.,  2.],
-       [ 2.,  2.,  2.],
-       [ 2.,  2.,  2.],
-       [ 2.,  2.,  2.]]), num.array([[ 2.,  2.,  2.],
-       [ 2.,  2.,  2.],
-       [ 2.,  2.,  2.],
-       [ 2.,  2.,  2.],
-       [ 2.,  2.,  2.]])]}}
-                                  
-
-        from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
-        
-        mesh = Mesh(nodes, triangles)
-        boundary_polygon = mesh.get_boundary_polygon()
-
-
-        # Subdivide into non-overlapping partitions
-
-        submesh = submesh_full(nodes, triangles, edges, \
-                            triangles_per_proc)
-
-        #print submesh
-
-
-        for i in range(3):
-            assert num.allclose(true_submesh['full_triangles'][i],submesh['full_triangles'][i])
-            assert num.allclose(true_submesh['full_nodes'][i],submesh['full_nodes'][i])
-        assert true_submesh['full_boundary'] == submesh['full_boundary']
-
-        # Add any extra ghost boundary layer information
-
-        submesh = submesh_ghost(submesh, mesh, triangles_per_proc)
-
-        for i in range(3):
-            assert num.allclose(true_submesh['ghost_triangles'][i],submesh['ghost_triangles'][i])
-            assert num.allclose(true_submesh['ghost_nodes'][i],submesh['ghost_nodes'][i])
-            assert num.allclose(true_submesh['ghost_commun'][i],submesh['ghost_commun'][i])
-
-        assert true_submesh['full_commun'] == submesh['full_commun']
-
-
-        # Order the quantities information to be the same as the triangle
-        # information
-
-
-        submesh = submesh_quantities(submesh, quantities, \
-                                 triangles_per_proc)
-
-
-
-        for key, value in true_submesh['ghost_quan'].iteritems():
-            for i in range(3):
-                assert num.allclose(true_submesh['ghost_quan'][key][i],submesh['ghost_quan'][key][i])
-                assert num.allclose(true_submesh['full_quan'][key][i],submesh['full_quan'][key][i])
-        
-
-        submesh["boundary_polygon"] = boundary_polygon
-
-
-        #print submesh
-
-#-------------------------------------------------------------
-
-if __name__ == "__main__":
-    suite = unittest.makeSuite(Test_Distribute_Mesh,'test')
-    runner = unittest.TextTestRunner()
-    runner.run(suite)
Index: /anuga_core/source/anuga_parallel/test_parallel_sw.py
===================================================================
--- /anuga_core/source/anuga_parallel/test_parallel_sw.py	(revision 7446)
+++ /anuga_core/source/anuga_parallel/test_parallel_sw.py	(revision 7447)
@@ -1,3 +1,2 @@
-#!/usr/bin/env python
 
 """Test a run of the sequential shallow water domain against
@@ -17,9 +16,11 @@
 import pypar
 
-from Numeric import allclose, array, zeros, Float, take, nonzero
+#from Numeric import allclose, array, zeros, Float, take, nonzero
+
+import numpy as num
 
 from anuga.pmesh.mesh_interface import create_mesh_from_regions
 
-from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
+from anuga.interface import rectangular_cross
 from anuga.abstract_2d_finite_volumes.pmesh2domain import pmesh_to_domain_instance
 
@@ -28,9 +29,9 @@
 from anuga.utilities.norms           import l1_norm, l2_norm, linf_norm
 
-from anuga.shallow_water import Domain
-from anuga.shallow_water import Reflective_boundary
-from anuga.shallow_water import Dirichlet_boundary
-from anuga.shallow_water import Time_boundary
-from anuga.shallow_water import Transmissive_boundary
+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
 
 
@@ -69,4 +70,5 @@
 def evolution_test(parallel=False):
 
+
     domain = pmesh_to_domain_instance(mesh_filename, Domain)
     domain.set_quantity('stage', Set_Stage(756000.0, 756500.0, 2.0))
@@ -95,8 +97,8 @@
     l2list = []
     linflist = []
-    l1norm = zeros(3, Float)
-    l2norm = zeros(3, Float)
-    linfnorm = zeros(3, Float)
-    recv_norm = zeros(3, Float)
+    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)
 
     #------------------------------------------------------------------------------
@@ -109,5 +111,5 @@
         
     for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
-        edges = take(domain.quantities[quantity].edge_values, nonzero(domain.tri_full_flag))
+        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])
@@ -128,9 +130,9 @@
                 #print edges[:,1]            
                 for p in range(1, numprocs):
-                    pypar.receive(p, recv_norm)
+                    recv_norm = pypar.receive(p)
                     l1norm += recv_norm
-                    pypar.receive(p, recv_norm)
+                    recv_norm = pypar.receive(p)
                     l2norm += recv_norm
-                    pypar.receive(p, recv_norm)
+                    recv_norm = pypar.receive(p)
                     linfnorm[0] = max(linfnorm[0], recv_norm[0])
                     linfnorm[1] = max(linfnorm[1], recv_norm[1])
Index: /anuga_core/source/anuga_parallel/test_parallel_sw_runup.py
===================================================================
--- /anuga_core/source/anuga_parallel/test_parallel_sw_runup.py	(revision 7446)
+++ /anuga_core/source/anuga_parallel/test_parallel_sw_runup.py	(revision 7447)
@@ -30,5 +30,5 @@
 from anuga.interface import rectangular_cross
 
-from parallel_api import distribute, myid, numprocs
+from anuga_parallel.parallel_api import distribute, myid, numprocs
 
 
@@ -74,10 +74,4 @@
 domain.set_default_order(1)        
 domain.set_quantities_to_be_stored(None)
-domain.set_maximum_allowed_speed(100) #FIXME (Ole): try to remove this
-
-# FIXME (Ole): Need tests where this is commented out
-domain.tight_slope_limiters = 0 # Backwards compatibility (14/4/7)
-domain.H0 = 0 # Backwards compatibility (6/2/7)
-domain.use_centroid_velocities = 0 # Backwards compatibility (7/5/8)
 
 
@@ -108,5 +102,5 @@
     try:
         k = domain.get_triangle_containing_point(point)
-        #print 'KKK',myid, k, domain.tri_full_flag[k]
+        print 'KKK',myid, k, domain.tri_full_flag[k]
         if domain.tri_full_flag[k] == 1:
             tri_ids.append(k)
@@ -116,6 +110,6 @@
         tri_ids.append(-2)
 
-#print myid, tri_ids
-#print myid, domain.tri_full_flag[tri_ids[0]]
+print myid, tri_ids
+print myid, domain.tri_full_flag[tri_ids[0]]
 
 
@@ -144,16 +138,29 @@
 
 
-G0 = [-0.19166666666666665, -0.19166666666666665, -0.1908726751282839, -0.18320314010752042, -0.18294289434947858, -0.17186692881708804, -0.16703039399924297, -0.16196038919122194, -0.15621633053131387, -0.15130021599977714, -0.13930978857215476, -0.18515941024930252, -0.19141974265470424, -0.19166563809769938, -0.19166666621987774, -0.19166666666616614, -0.19166666666616614, -0.19166666666616614, -0.19163936679820112, -0.19092472615522754, -0.19101180444286325, -0.19133150862916881, -0.1914019526842341, -0.19134927149082989, -0.19146947464147604, -0.19165471548476315, -0.19166666444742567, -0.19166666666632831, -0.19166666666632831, -0.19166666666632831, -0.19166666666632831, -0.19166666666632831, -0.19166666666632831, -0.19166666666632831, -0.19166666666632831, -0.19166666666632831, -0.19166666666632831, -0.19166666666632831, -0.19166666666632831, -0.19166666666632831, -0.19166666666632831, -0.19166666666632831, -0.19166666666632831, -0.19166666666632831, -0.19166666666632831, -0.19166666666632831, -0.19166666666632831, -0.19166666666632831, -0.19166666666632831, -0.19166666666632831, -0.19166666666632831]
+G0 = [-0.19166666666666665, -0.19166666666666665, -0.19166666666666665, -0.18995135073536706, -0.17967620310033938, -0.17592302731546394, -0.17196188698352638, -0.16050146192256606, -0.15408782540491983, -0.15120360396251856, -0.15272836248974528, -0.18515612571281359, -0.19089313596572838, -0.19091930090793735, -0.19092564666264678, -0.19093961544356891, -0.19096113497288297, -0.19100903501132768, -0.19101957821753363, -0.19100676219446935, -0.19098186082883836, -0.19095020375973476, -0.19090160340144657, -0.19060035225639324, -0.19067236768901255, -0.19072092369650068, -0.19084725051273566, -0.19087147781405342, -0.19089622376237744, -0.19090954329991847, -0.19091808778328567, -0.19092562302462121, -0.19093120010640158, -0.19093455662311021, -0.1909359528144898, -0.19093587102253487, -0.19093485355514328, -0.19093340185413665, -0.19093189648362133, -0.19093071766205927, -0.19093005702136678, -0.19092994761677365, -0.19093031876104916, -0.19093104326472229, -0.19093197773630957, -0.19093299111521683, -0.19093398080658583, -0.19093487865230929, -0.19093564970049512, -0.1909362865237911, -0.19093680142438879]
 
-G1 = [-0.29166666666666669, -0.29160604748111429, -0.28349212083663766, -0.26637122760395054, -0.25078182919408626, -0.22796566989512818, -0.21386588163016981, -0.20291676395026542, -0.19288397535979634, -0.18529863721918491, -0.17833464440180594, -0.16041289672597714, -0.15844675190030885, -0.18498595153509523, -0.20472977209518117, -0.21430915468432699, -0.21656349614594508, -0.2143710033505706, -0.21047592509949745, -0.20690277124897924, -0.20383621964071927, -0.2009312504622697, -0.19832906347760781, -0.19715218778669685, -0.19713511857310209, -0.19760258124272423, -0.19821523236790017, -0.19886735995941604, -0.19945089045017395, -0.19990124934465239, -0.20016164831872438, -0.20035891530554095, -0.20053280091253389, -0.20062879519403856, -0.20061453716771055, -0.20050516254700712, -0.20034457406881809, -0.2001804893235041, -0.20004787348381434, -0.19996313076460828, -0.19992645686226715, -0.19992808844809359, -0.19995462843450704, -0.19999336711919727, -0.20003430609120498, -0.20007059456815027, -0.20009823569348062, -0.20011560727571823, -0.20012297891084577, -0.20012200982886325, -0.2001151929510443]
+
+
+
+
+G1 = [-0.29166666666666669, -0.29166666666666669, -0.28685941525964553, -0.27057862587984083, -0.24118168337207896, -0.22173239931821093, -0.20976527589038554, -0.19971857867921899, -0.19173249780540819, -0.18487858041944857, -0.1780033786977564, -0.17215410368358283, -0.17155182647714193, -0.18655941293682041, -0.19977612728273281, -0.2068993409588821, -0.2095974217439926, -0.20961276275180968, -0.20804141273268489, -0.20583271698671743, -0.20358656706459732, -0.20159086143490923, -0.19997988026328414, -0.19903422512577046, -0.19862770982003813, -0.1984841836620021, -0.19852978191537585, -0.19875351136877187, -0.19912658475790856, -0.19955814714527378, -0.20002959678495849, -0.20032623461817603, -0.20045056162352595, -0.20047228698692546, -0.20043159564517379, -0.20035498975577987, -0.2002645931437364, -0.20017748770076391, -0.2001040624974407, -0.20005138031775113, -0.20002019418902675, -0.20000780866902876, -0.200010236268512, -0.20002255766703306, -0.2000398188461642, -0.20005785100286463, -0.2000737026875791, -0.20008570712104165, -0.20009331547005016, -0.20009681408763169, -0.20009700971531758]
+
+
+
+
 
         
-G2 = [-0.39166666666666666, -0.38244067921398905, -0.33350466136630463, -0.29771023004255276, -0.27605439066140891, -0.25986156218997503, -0.24502185018573649, -0.23179262432952102, -0.21981564668803996, -0.20870707082936543, -0.19877739883776596, -0.18980922837977954, -0.1730801167400583, -0.16306400164063853, -0.17798470933316624, -0.19295540736943456, -0.20236705173335867, -0.20695767548229582, -0.20841025868691554, -0.20792102171307641, -0.20655350005113712, -0.2049200252815718, -0.20310627030300929, -0.20105983339290284, -0.19937394568595421, -0.1985391750807548, -0.19836389978207072, -0.19850305023816406, -0.19877764028800521, -0.19910928130800573, -0.19943705712074122, -0.19970344172627957, -0.19991076989513878, -0.20010020127344161, -0.20025937786728282, -0.20035087292654716, -0.20035829921368681, -0.20029606557358018, -0.2001960691549278, -0.20009096093557327, -0.20000371608352077, -0.19994495433034862, -0.1999153566524963, -0.19990981826568235, -0.19992106419902292, -0.19994189853498176, -0.19996624091198109, -0.19998946016949451, -0.20000842303436642, -0.20002144460691332, -0.20002815561319481]
+G2 = [-0.39166666666666666, -0.38165411084088413, -0.33031114311718357, -0.29490845791170966, -0.27401187763459972, -0.25720432221568812, -0.24180370475026186, -0.22803625608700034, -0.21607642473560867, -0.2059173963924967, -0.19732372829977357, -0.18899569210124012, -0.18075177888497074, -0.1758458811698074, -0.18266970744342398, -0.19232195386649067, -0.199153652922028, -0.20292588217669899, -0.2046835952536466, -0.20515202619935408, -0.20476365635086183, -0.20384587744356614, -0.20266293067781077, -0.20143368989447288, -0.20039011693774927, -0.19970226213297601, -0.1992938050267305, -0.19906505599862009, -0.19899924116319651, -0.19908893561605648, -0.1993045458861577, -0.19960344595715143, -0.19989216726729134, -0.20008423737677972, -0.20018439752920908, -0.20022167063422625, -0.20021700109190652, -0.20018483837701212, -0.2001375242307058, -0.20008571798103356, -0.20003819123906194, -0.20000089191108147, -0.1999759332393749, -0.19996300730369299, -0.19996026140162992, -0.19996490481461601, -0.19997388431092483, -0.19998445575566867, -0.19999452242082996, -0.20000274678941998, -0.20000850192468142]
+
+
+
         
-G3 = [-0.44166666666666665, -0.37631169657400343, -0.33000044342859486, -0.30586045469008505, -0.28843572253009925, -0.27215308978603797, -0.25712951540331214, -0.2431608296216613, -0.23032023651386366, -0.21845468734566184, -0.20735123704254327, -0.19740397194806383, -0.18598295640643708, -0.16675980728412546, -0.1695157503295662, -0.18328608714846414, -0.19485758921965809, -0.2023136827680253, -0.20625610365424957, -0.20758116234856749, -0.20721445398906116, -0.20603406830101498, -0.20450262810357958, -0.2026769581530177, -0.20074012124164839, -0.19931160538111795, -0.1986360630235105, -0.19848511941161767, -0.1986009104317408, -0.19885490669335967, -0.19916542732465842, -0.19946678238299068, -0.19971209593781808, -0.19991912886144106, -0.20010584307496243, -0.20024959409137078, -0.20032160254398582, -0.2003158316568579, -0.20025051539345248, -0.20015561158283834, -0.20005952955567172, -0.19998144295750273, -0.19992977821660762, -0.19990457708729134, -0.19990104785520091, -0.19991257153955966, -0.19993258231861782, -0.19995548502852453, -0.19997700760886039, -0.19999429663472679, -0.20000588800224417]
+G3 = [-0.44166666666666665, -0.37566630156374725, -0.32907403981602412, -0.30514246507845882, -0.28736474524309508, -0.27062836126322631, -0.25490557254596713, -0.24015946060857596, -0.22672499664646506, -0.21493194009348415, -0.20488909335665775, -0.19594643601723874, -0.18708799192006037, -0.17881263621051871, -0.17829280002483108, -0.18600506642077144, -0.19389819408710718, -0.19941728321954535, -0.20269142693734038, -0.20431177320907121, -0.20475261565846536, -0.20438149942198236, -0.20350472280402784, -0.20238009732809753, -0.20123125943892334, -0.20030079794157421, -0.19967828320080519, -0.19928934067283594, -0.19907912341626033, -0.199034739569321, -0.19914305100399451, -0.19936563222000359, -0.19965218325198875, -0.19990375153648701, -0.20006953873503872, -0.20015884586614877, -0.20019243432639855, -0.20018670241799616, -0.20015514524359709, -0.2001098495335151, -0.2000610924794913, -0.20001718731119694, -0.19998313833853548, -0.19996070482302816, -0.1999494922032494, -0.1999476393706516, -0.1999524591476477, -0.19996109366774711, -0.19997102151684798, -0.19998034301901382, -0.1999878666221423]
+
 
 
 # Only compare those that belong to this process id
 G = [G0, G1, G2, G3]
+
 
 
@@ -166,4 +173,5 @@
 
 
+
 if success:
     print 'Successful completion on processor ',myid
