Index: /anuga_core/source/anuga_parallel/build_submesh.py
===================================================================
--- /anuga_core/source/anuga_parallel/build_submesh.py	(revision 3817)
+++ /anuga_core/source/anuga_parallel/build_submesh.py	(revision 3818)
@@ -512,5 +512,5 @@
 #  *) A dictionary containing the full_triangles,
 # full_nodes, full_boundary, ghost_triangles, ghost_nodes,
-# ghost_boundary, ghost_commun and full_commun is returned.
+# ghost_boundary, ghost_commun and full_commun and true boundary polygon is returned.
 #
 #########################################################
@@ -520,7 +520,9 @@
 
     # Temporarily build the mesh to find the neighbouring
-    # triangles
+    # triangles and true boundary polygon
 
     mesh = Mesh(nodes, triangles)
+    boundary_polygon = mesh.get_boundary_polygon()
+    
 
     # Subdivide into non-overlapping partitions
@@ -539,4 +541,5 @@
                                  triangles_per_proc)
 
+    submesh["boundary_polygon"] = boundary_polygon
     return submesh
 
Index: /anuga_core/source/anuga_parallel/parallel_api.py
===================================================================
--- /anuga_core/source/anuga_parallel/parallel_api.py	(revision 3817)
+++ /anuga_core/source/anuga_parallel/parallel_api.py	(revision 3818)
@@ -39,4 +39,17 @@
     """ Distribute the domain to all processes
     """
+
+
+    # FIXME: Dummy assignment (until boundaries are refactored to
+    # be independent of domains until they are applied)
+    bdmap = {}
+    for tag in domain.get_boundary_tags():
+        bdmap[tag] = None
+    
+    
+    domain.set_boundary(bdmap)
+
+
+
 
     if not pypar_available: return domain # Bypass
@@ -142,6 +155,7 @@
            pmesh_divide_metis(domain, numprocs)
 
+
     # Build the mesh that should be assigned to each processor,
-    # this includes ghost nodes and the communicaiton pattern
+    # this includes ghost nodes and the communication pattern
     print 'Build submeshes'    
     submesh = build_submesh(nodes, triangles, boundary,\
Index: /anuga_core/source/anuga_parallel/test_parallel_sw_runup.py
===================================================================
--- /anuga_core/source/anuga_parallel/test_parallel_sw_runup.py	(revision 3817)
+++ /anuga_core/source/anuga_parallel/test_parallel_sw_runup.py	(revision 3818)
@@ -15,7 +15,11 @@
 #------------------------------------------------------------------------------
 
+from Numeric import allclose
+
 from anuga.pmesh.mesh_interface import create_mesh_from_regions
 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
 from anuga.utilities.numerical_tools import ensure_numeric
+from anuga.utilities.polygon import is_inside_polygon
+
 from anuga.shallow_water import Domain
 from anuga.shallow_water import Reflective_boundary
@@ -48,6 +52,13 @@
 
 
+#--------------------------------------------------------------------------
+# Create the parallel domain
+#--------------------------------------------------------------------------
+
+domain = distribute(domain, verbose=True)
+
+
 #------------------------------------------------------------------------------
-# Setup boundary conditions
+# Setup boundary conditions (MUST currently happen after domain has been distributed)
 #------------------------------------------------------------------------------
 
@@ -56,14 +67,5 @@
 
 # Associate boundary tags with boundary objects
-domain.set_boundary({'left': Br, 'right': Bd, 'top': Br, 'bottom': Br})
-
-
-#--------------------------------------------------------------------------
-# Create the parallel domain
-#--------------------------------------------------------------------------
-
-# If this is commented out, the code works. Otherwise the interpolated
-# timeseries are wrong - even when run on one processor.
-# domain = distribute(domain, verbose=True)
+domain.modify_boundary({'left': Br, 'right': Bd, 'top': Br, 'bottom': Br})
 
 
@@ -74,11 +76,21 @@
 interpolation_points = [[0.4,0.5], [0.6,0.5], [0.8,0.5], [0.9,0.5]]
 gauge_values = []
-for _ in interpolation_points:
-    gauge_values.append([])
+local_interpolation_points = []
+for i, point in enumerate(interpolation_points):
+    gauge_values.append([]) # Empty list for timeseries
+
+    if is_inside_polygon(point, domain.get_boundary_polygon()):
+
+        # FIXME: One point appears on multiple processes
+        # Need to get true boundary somehow
+        
+        print 'P%d: point=[%f,%f]' %(myid, point[0], point[1])
+        local_interpolation_points.append(i)
+
 
 time = []
 
 for t in domain.evolve(yieldstep = 0.1, finaltime = 5.0):
-    domain.write_time()
+    #domain.write_time()
 
     
@@ -94,8 +106,4 @@
 
 for i, (x,y) in enumerate(interpolation_points):
-    gauge_values[i] = gauge_values[i][:51] # In case an extra step was taken
-#    print i, gauge_values[i]
-#    print 
-
         
     try:
@@ -125,5 +133,4 @@
 G0 = ensure_numeric([-0.20000000000000001, -0.19999681443389281, -0.1986192343695303, -0.19147413648863046, -0.19132688908678019, -0.17642317476621105, -0.167376262630034, -0.16192452887426961, -0.15609171725778803, -0.15127107084302249, -0.14048864340360018, -0.19296484125327093, -0.19997006390580363, -0.19999999999937063, -0.19999999999937063, -0.19999999999938772, -0.19999999999938772, -0.19999999999938772, -0.19999999999938772, -0.19974288463035494, -0.19951636867991712, -0.19966301435195755, -0.19981082259800226, -0.19978575003960128, -0.19992942471933109, -0.19999999931029933, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989])
 
-
 G1 = ensure_numeric([-0.29999999999999993, -0.29988962537199199, -0.29293904425532025, -0.28329367722887888, -0.25999146407696289, -0.22613875068011896, -0.21190705052094994, -0.19900707995208217, -0.18876305176191882, -0.18132447501091936, -0.17395459512711151, -0.15562414200985644, -0.16212999953643359, -0.18964422820514618, -0.20871181844346975, -0.21672207791083464, -0.21774940291862779, -0.21482868050219833, -0.21057786776704043, -0.20649663432591045, -0.20294932949211578, -0.19974459897911329, -0.19733648772704043, -0.19641404599824669, -0.19654095699184146, -0.19709942852191994, -0.19780873983410741, -0.19853259125123518, -0.19916495938961168, -0.19965391267799168, -0.19993539587158982, -0.2001383705551133, -0.20029344332295113, -0.20035349748150011, -0.20029886541561631, -0.20015541958920294, -0.19997273066429103, -0.19979879448668514, -0.19966016997024041, -0.19957558009501869, -0.19955725674938532, -0.19958083002853366, -0.19961752462568647, -0.19965296611330258, -0.19968998132634594, -0.19972532942208607, -0.19975372922008239, -0.19977196116929855, -0.19977951443660594, -0.19977792107284789, -0.19976991595502003])
 
@@ -132,19 +139,13 @@
 G3 = ensure_numeric([-0.45000000000000001, -0.37713945714588398, -0.33029565026933816, -0.30598209033945367, -0.28847101155177313, -0.27211191064563195, -0.25701544058818926, -0.24298945948410997, -0.23010402733784807, -0.21820351802867713, -0.20709938367218383, -0.19719881806182216, -0.18568281604361933, -0.16828653906676322, -0.16977310768235579, -0.1832707289594605, -0.19483524345250974, -0.20233480051649216, -0.20630757214159207, -0.20763927857964531, -0.20724458160595791, -0.20599191745446047, -0.20438329669495012, -0.20256105512496606, -0.20071269486729407, -0.19934403619901719, -0.19866860191898347, -0.19849975056296071, -0.19860870923007437, -0.19885838217851401, -0.19916422433758982, -0.19946861981642039, -0.19972267778871666, -0.19993013816258154, -0.20011063428833351, -0.20024891930311628, -0.20031882555219671, -0.20031326268593497, -0.20024881068472311, -0.20015443214902759, -0.20005669097631221, -0.19997542564643309, -0.19992564006223304, -0.19990746148869892, -0.19990923999172872, -0.19991956416813192, -0.19993484556273733, -0.1999538628054662, -0.19997381636620407, -0.19999130900268777, -0.20000388227457688])
 
-from Numeric import allclose
-msg = 'Computed time series and reference time series are different: %s'\
-      %(gauge_values[0]-G0)
-assert allclose(gauge_values[0], G0), msg
 
-msg = 'Computed time series and reference time series are different: %s'\
-      %(gauge_values[1]-G1)
-assert allclose(gauge_values[1], G1), msg
+# Only compare those that belong to this process id
+G = [G0, G1, G2, G3]
 
-msg = 'Computed time series and reference time series are different: %s'\
-      %(gauge_values[2]-G2)
-assert allclose(gauge_values[2], G2), msg
+for i in local_interpolation_points:
+    msg = 'P%d, point #%d: Computed time series and reference time series are different: %s'\
+          %(myid, i, gauge_values[i]-G[i])
+    assert allclose(gauge_values[i], G[i]), msg
 
-msg = 'Computed time series and reference time series are different: %s'\
-      %(gauge_values[3]-G3)
-assert allclose(gauge_values[3], G3), msg        
+print 'P%d completed succesfully using points = %s' %(myid, local_interpolation_points)
 
