Index: anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py
===================================================================
--- anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py	(revision 3927)
+++ anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py	(revision 3928)
@@ -47,6 +47,6 @@
                  processor=0,
                  numproc=1,
-                 number_of_full_nodes=0,
-                 number_of_full_triangles=0):
+                 number_of_full_nodes=None,
+                 number_of_full_triangles=None):
 
 
@@ -87,6 +87,11 @@
 
         # Initialise underlying mesh structure
-        Mesh.__init__(self, coordinates, triangles, boundary,
-                      tagged_elements, geo_reference, use_inscribed_circle,
+        Mesh.__init__(self, coordinates, triangles,
+                      boundary=boundary,
+                      tagged_elements=tagged_elements,
+                      geo_reference=geo_reference,
+                      use_inscribed_circle=use_inscribed_circle,
+                      number_of_full_nodes=number_of_full_nodes,
+                      number_of_full_triangles=number_of_full_triangles,
                       verbose=verbose)
 
@@ -136,7 +141,4 @@
         self.numproc = numproc
 
-        self.number_of_full_nodes=number_of_full_nodes
-        self.number_of_full_triangles=number_of_full_triangles
-        
 
         # Setup Communication Buffers
@@ -156,5 +158,5 @@
         # =1 for full
         # =0 for ghost
-        N=self.number_of_elements
+        N = len(self) #number_of_elements
         self.tri_full_flag = ones(N, Int)
         for i in self.ghost_recv_dict.keys():
@@ -206,5 +208,5 @@
         #MH310505 To avoid calculating the flux across each edge twice, keep an integer (boolean) array,
         #to be used during the flux calculation
-        N=self.number_of_elements
+        N = len(self) #number_of_triangles
         self.already_computed_flux = zeros((N, 3), Int)
 
@@ -923,5 +925,5 @@
         from Numeric import ones, sum, equal, Float
 
-        N = self.number_of_elements
+        N = len(self) #number_of_triangles
         d = len(self.conserved_quantities)
 
Index: anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py
===================================================================
--- anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py	(revision 3927)
+++ anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py	(revision 3928)
@@ -56,4 +56,6 @@
     #input
     def __init__(self, coordinates, triangles,
+                 number_of_full_nodes=None,
+                 number_of_full_triangles=None,                 
                  geo_reference=None,
                  verbose=False):
@@ -65,4 +67,11 @@
         origin is a 3-tuple consisting of UTM zone, easting and northing.
         If specified coordinates are assumed to be relative to this origin.
+
+
+        number_of_full_nodes and number_of_full_triangles relate to
+        parallelism when each mesh has an extra layer of ghost points and
+        ghost triangles attached to the end of the two arrays.
+        In this case it is usefull to specify the number of real (called full)
+        nodes and triangles. If omitted they will default to all.
         """
 
@@ -71,5 +80,29 @@
         self.triangles = array(triangles,Int)
         self.coordinates = array(coordinates,Float)
-        
+
+        # Register number of elements and nodes 
+        self.number_of_triangles = N = self.triangles.shape[0]
+        self.number_of_nodes = self.coordinates.shape[0]        
+
+        
+
+        if number_of_full_nodes is None:
+            self.number_of_full_nodes=self.number_of_nodes
+        else:
+            assert int(number_of_full_nodes)
+            self.number_of_full_nodes=number_of_full_nodes            
+
+
+        if number_of_full_triangles is None:
+            self.number_of_full_triangles=self.number_of_triangles           
+        else:
+            assert int(number_of_full_triangles)            
+            self.number_of_full_triangles=number_of_full_triangles
+        
+        
+        #print self.number_of_full_nodes, self.number_of_nodes
+        #print self.number_of_full_triangles, self.number_of_triangles
+        
+            
 
         # FIXME: this stores a geo_reference, but when coords are returned
@@ -94,7 +127,4 @@
         assert max(max(self.triangles)) <= self.coordinates.shape[0], msg
 
-
-        # Register number of elements (N)
-        self.number_of_elements = N = self.triangles.shape[0]
 
         # FIXME: Maybe move to statistics?
@@ -175,5 +205,6 @@
 
     def __len__(self):
-        return self.number_of_elements
+        return self.number_of_triangles
+    
 
     def __repr__(self):
@@ -213,6 +244,8 @@
         """
 
+        N = self.number_of_full_nodes
+
         if unique is True:
-            V = self.coordinates
+            V = self.coordinates[:N,:]
             if absolute is True:
                 if not self.geo_reference.is_absolute():
@@ -227,7 +260,7 @@
             if not self.geo_reference.is_absolute():
             
-                V0 = self.geo_reference.get_absolute(V[:,0:2])
-                V1 = self.geo_reference.get_absolute(V[:,2:4])
-                V2 = self.geo_reference.get_absolute(V[:,4:6])
+                V0 = self.geo_reference.get_absolute(V[:N,0:2])
+                V1 = self.geo_reference.get_absolute(V[:N,2:4])
+                V2 = self.geo_reference.get_absolute(V[:N,4:6])
 
                 # This does double the memory need 
@@ -262,5 +295,5 @@
         #FIXME (Ole) - oh yes they should
 
-        N = self.number_of_elements
+        N = self.number_of_triangles
         vertex_coordinates = zeros((N, 6), Float)
 
@@ -279,8 +312,12 @@
         """
 
-        if (indices ==  None):
-            indices = range(len(self))  #len(self)=number of elements
+        N = self.number_of_full_triangles
+
+        if indices is None:
+            #indices = range(len(self))  #len(self)=number of elements
+            indices = range(N)
 
         return  take(self.triangles, indices)
+    
 
     #FIXME - merge these two (get_vertices and get_triangles)
@@ -314,4 +351,5 @@
         return unique_verts.keys()
 
+
     def build_vertexlist(self):
         """Build vertexlist index by vertex ids and for each entry (point id)
@@ -327,5 +365,5 @@
 
         vertexlist = [None]*len(self.coordinates)
-        for i in range(self.number_of_elements):
+        for i in range(self.number_of_triangles):
 
             a = self.triangles[i, 0]
Index: anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py
===================================================================
--- anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py	(revision 3927)
+++ anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py	(revision 3928)
@@ -67,4 +67,6 @@
                  tagged_elements=None,
                  geo_reference=None,
+                 number_of_full_nodes=None,
+                 number_of_full_triangles=None,
                  use_inscribed_circle=False,
                  verbose=False):
@@ -80,9 +82,14 @@
 
         General_mesh.__init__(self, coordinates, triangles,
-                              geo_reference, verbose=verbose)
+                              number_of_full_nodes=\
+                              number_of_full_nodes,
+                              number_of_full_triangles=\
+                              number_of_full_triangles,
+                              geo_reference=geo_reference,
+                              verbose=verbose)
 
         if verbose: print 'Initialising mesh'         
 
-        N = self.number_of_elements
+        N = len(self) #Number_of_triangles
 
         self.use_inscribed_circle = use_inscribed_circle
@@ -181,5 +188,5 @@
     def set_to_inscribed_circle(self,safety_factor = 1):
         #FIXME phase out eventually
-        N = self.number_of_elements
+        N = self.number_of_triangles
         V = self.vertex_coordinates
 
@@ -225,5 +232,5 @@
         #to left hand side edge (facing neighbouring triangle)
 
-        N = self.number_of_elements
+        N = len(self) #Number_of_triangles
         neighbourdict = {}
         for i in range(N):
@@ -289,5 +296,5 @@
         """
 
-        N = self.number_of_elements
+        N = len(self) #Number of triangles
         for i in range(N):
             #Find all neighbouring volumes that are not boundaries
@@ -314,5 +321,5 @@
         if boundary is None:
             boundary = {}
-            for vol_id in range(self.number_of_elements):
+            for vol_id in range(len(self)):
                 for edge_id in range(0, 3):
                     if self.neighbours[vol_id, edge_id] < 0:
@@ -330,5 +337,5 @@
 
             #Check that all boundary segments are assigned a tag
-            for vol_id in range(self.number_of_elements):
+            for vol_id in range(len(self)):
                 for edge_id in range(0, 3):
                     if self.neighbours[vol_id, edge_id] < 0:
@@ -374,5 +381,5 @@
 
                 msg = 'Not all elements exist. '
-                assert max(tagged_elements[tag]) < self.number_of_elements, msg
+                assert max(tagged_elements[tag]) < len(self), msg
         #print "tagged_elements", tagged_elements
         self.tagged_elements = tagged_elements
@@ -584,5 +591,5 @@
         from anuga.utilities.numerical_tools import anglediff
 
-        N = self.number_of_elements
+        N = len(self)
         #Get x,y coordinates for all vertices for all triangles
         V = self.get_vertex_coordinates()
@@ -732,5 +739,5 @@
         str =  '------------------------------------------------\n'
         str += 'Mesh statistics:\n'
-        str += '  Number of triangles = %d\n' %self.number_of_elements
+        str += '  Number of triangles = %d\n' %len(self)
         str += '  Extent [m]:\n'
         str += '    x in [%f, %f]\n' %(min(x), max(x))
Index: anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py
===================================================================
--- anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py	(revision 3927)
+++ anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py	(revision 3928)
@@ -29,5 +29,5 @@
 
         if vertex_values is None:
-            N = domain.number_of_elements
+            N = len(domain) # number_of_elements
             self.vertex_values = zeros((N, 3), Float)
         else:
@@ -42,7 +42,7 @@
                   %N
             msg += 'number of elements in specified domain (%d).'\
-                   %domain.number_of_elements
-
-            assert N == domain.number_of_elements, msg
+                   %len(domain)
+
+            assert N == len(domain), msg
 
         self.domain = domain
@@ -1135,5 +1135,5 @@
         """
         integral = 0
-        for k in range(self.domain.number_of_elements):
+        for k in range(len(self.domain)):
             area = self.domain.areas[k]
             qc = self.centroid_values[k]
@@ -1165,5 +1165,5 @@
         #flux calculations and forcing functions
 
-        N = domain.number_of_elements
+        N = len(domain) # number_of_triangles
         self.explicit_update = zeros(N, Float )
         self.semi_implicit_update = zeros(N, Float )
@@ -1293,5 +1293,5 @@
 
     #Check each triangle
-    for k in range(quantity.domain.number_of_elements):
+    for k in range(len(quantity.domain)):
         #Centroid coordinates
         x, y = quantity.domain.centroid_coordinates[k]
@@ -1389,5 +1389,5 @@
     from Numeric import zeros, Float
 
-    N = quantity.domain.number_of_elements
+    N = len(quantity.domain)
 
     beta_w = quantity.domain.beta_w
Index: anuga_core/source/anuga/abstract_2d_finite_volumes/test_general_mesh.py
===================================================================
--- anuga_core/source/anuga/abstract_2d_finite_volumes/test_general_mesh.py	(revision 3927)
+++ anuga_core/source/anuga/abstract_2d_finite_volumes/test_general_mesh.py	(revision 3928)
@@ -28,5 +28,4 @@
         domain = General_mesh(points, vertices)
 
-        print domain.__class__
 
         assert allclose(domain.get_vertex_coordinates(unique=True), domain.coordinates)
Index: anuga_core/source/anuga/advection/advection.py
===================================================================
--- anuga_core/source/anuga/advection/advection.py	(revision 3927)
+++ anuga_core/source/anuga/advection/advection.py	(revision 3928)
@@ -169,5 +169,5 @@
         from anuga.config import max_timestep
 
-        N = self.number_of_elements
+        N = len(self)
 
         neighbours = self.neighbours
@@ -260,5 +260,5 @@
         from weave import converters
 
-        N = self.number_of_elements
+        N = len(self)
 
 
Index: anuga_core/source/anuga/shallow_water/data_manager.py
===================================================================
--- anuga_core/source/anuga/shallow_water/data_manager.py	(revision 3927)
+++ anuga_core/source/anuga/shallow_water/data_manager.py	(revision 3928)
@@ -213,10 +213,18 @@
         #Create filename
         self.filename = create_filename(domain.get_datadir(),
-            domain.get_name(), extension)
+                                        domain.get_name(), extension)
 
         #print 'F', self.filename
         self.timestep = 0
-        self.number_of_volumes = len(domain)
         self.domain = domain
+        
+
+
+        # Exclude ghosts in case this is a parallel domain
+        self.number_of_nodes = domain.number_of_full_nodes        
+        self.number_of_volumes = domain.number_of_full_triangles
+        #self.number_of_volumes = len(domain)        
+
+
 
 
@@ -277,5 +285,5 @@
             if hasattr(domain, 'texture'):
                 fid.texture = domain.texture
-        #else:
+            #else:
             #    fid.texture = 'None'
 
@@ -291,4 +299,10 @@
             if domain.smooth is True:
                 fid.createDimension('number_of_points', len(domain.vertexlist))
+                #fid.createDimension('number_of_points', self.number_of_nodes)
+
+                # FIXME(Ole): This will cause sww files for paralle domains to
+                # have ghost nodes stored (but not used by triangles).
+                # To clean this up, we have to change get_vertex_values and friends in
+                # quantity.py (but I can't be bothered right now)
             else:
                 fid.createDimension('number_of_points', 3*self.number_of_volumes)
Index: anuga_core/source/anuga/shallow_water/shallow_water_domain.py
===================================================================
--- anuga_core/source/anuga/shallow_water/shallow_water_domain.py	(revision 3927)
+++ anuga_core/source/anuga/shallow_water/shallow_water_domain.py	(revision 3928)
@@ -112,6 +112,6 @@
                  processor=0,
                  numproc=1,
-                 number_of_full_nodes=0,
-                 number_of_full_triangles=0):
+                 number_of_full_nodes=None,
+                 number_of_full_triangles=None):
 
 
@@ -746,5 +746,5 @@
     import sys
 
-    N = domain.number_of_elements
+    N = len(domain) # number_of_triangles
 
     #Shortcuts
@@ -820,5 +820,5 @@
     import sys
 
-    N = domain.number_of_elements
+    N = len(domain) # number_of_triangles
 
     #Shortcuts
@@ -848,5 +848,5 @@
     import sys
 
-    N = domain.number_of_elements
+    N = len(domain) # number_of_triangles
 
     #Shortcuts
@@ -960,6 +960,6 @@
 
     #Update
-    #FIXME: Modify accroditg to c-version - or discard altogether.
-    for k in range(domain.number_of_elements):
+    #FIXME: Modify according to c-version - or discard altogether.
+    for k in range(len(domain)):
 
         if hc[k] < domain.minimum_allowed_height:
@@ -999,5 +999,5 @@
     """
 
-    N = domain.number_of_elements
+    N = len(domain)
     beta_h = domain.beta_h
 
@@ -1069,5 +1069,5 @@
     """
 
-    N = domain.number_of_elements
+    N = len(domain) # number_of_triangles
     beta_h = domain.beta_h
 
@@ -1114,5 +1114,5 @@
     hvbar = h_limiter(domain)
 
-    for k in range(domain.number_of_elements):
+    for k in range(len(domain)):
         #Compute maximal variation in bed elevation
         #  This quantitiy is
@@ -1402,5 +1402,5 @@
     g = domain.g
 
-    for k in range(domain.number_of_elements):
+    for k in range(len(domain)):
         avg_h = sum( h[k,:] )/3
 
@@ -1454,5 +1454,5 @@
     ymom_update = domain.quantities['ymomentum'].semi_implicit_update
 
-    N = domain.number_of_elements
+    N = len(domain)
     eps = domain.minimum_allowed_height
     g = domain.g
@@ -1489,5 +1489,5 @@
     ymom_update = ymom.semi_implicit_update
 
-    N = domain.number_of_elements
+    N = len(domain)
     eps = domain.minimum_allowed_height
     g = domain.g
@@ -1516,5 +1516,5 @@
     ymom_update = ymom.explicit_update
 
-    N = domain.number_of_elements
+    N = len(domain)
     eps = domain.minimum_allowed_height
     g = domain.g
@@ -1543,5 +1543,5 @@
     ymom_update = domain.quantities['ymomentum'].semi_implicit_update
 
-    N = domain.number_of_elements
+    N = len(domain) # number_of_triangles
     eps = domain.minimum_allowed_height
     g = domain.g #Not necessary? Why was this added?
@@ -1681,5 +1681,5 @@
         ymom_update = domain.quantities['ymomentum'].explicit_update
 
-        N = domain.number_of_elements
+        N = len(domain) # number_of_triangles
         t = domain.time
 
@@ -1767,5 +1767,5 @@
     #Computed linear combination between constant stages and and
     #stages parallel to the bed elevation.
-    for k in range(domain.number_of_elements):
+    for k in range(len(domain)):
         #Compute maximal variation in bed elevation
         #  This quantitiy is
@@ -2034,5 +2034,5 @@
     from Numeric import zeros, Float
 
-    N = domain.number_of_elements
+    N = len(domain) # number_of_triangles
 
     #Shortcuts
Index: anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py
===================================================================
--- anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py	(revision 3927)
+++ anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py	(revision 3928)
@@ -1311,5 +1311,5 @@
         const = eta_w*rho_a/rho_w
 
-        N = domain.number_of_elements
+        N = len(domain) # number_of_triangles
 
         xc = domain.get_centroid_coordinates()
@@ -1427,5 +1427,5 @@
         const = eta_w*rho_a/rho_w
 
-        N = domain.number_of_elements
+        N = len(domain) # number_of_triangles
 
         t = domain.time
@@ -1922,5 +1922,5 @@
         #Assert that quantities are conserved
         from Numeric import sum
-        for k in range(mesh.number_of_elements):
+        for k in range(len(mesh)):
             assert allclose (ref_centroid_values[k],
                              sum(stage.vertex_values[k,:])/3)
@@ -1944,10 +1944,10 @@
 
         #Assert that all vertex quantities have changed
-        for k in range(mesh.number_of_elements):
+        for k in range(len(mesh)):
             #print ref_vertex_values[k,:], stage.vertex_values[k,:]
             assert not allclose (ref_vertex_values[k,:], stage.vertex_values[k,:])
         #and assert that quantities are still conserved
         from Numeric import sum
-        for k in range(mesh.number_of_elements):
+        for k in range(len(mesh)):
             assert allclose (ref_centroid_values[k],
                              sum(stage.vertex_values[k,:])/3)
@@ -2212,12 +2212,4 @@
 
             #print t, volume, initial_volume
-
-
-            #if not allclose (volume, initial_volume):
-            #    print 't==4.05'
-            #    for k in range(domain.number_of_elements):
-            #       pass
-            #       print domain.quantities['stage'].centroid_values[k] -\
-            #             ref_centroid_values[k]
 
             assert allclose (volume, initial_volume)
Index: anuga_core/source/anuga_parallel/parallel_api.py
===================================================================
--- anuga_core/source/anuga_parallel/parallel_api.py	(revision 3927)
+++ anuga_core/source/anuga_parallel/parallel_api.py	(revision 3928)
@@ -171,8 +171,8 @@
 
     for p in range(numprocs):
+        N = len(submesh['ghost_nodes'][p])                
         M = len(submesh['ghost_triangles'][p])
-        print 'There are %d ghost triangles on proc %d' %(M, p)
-        N = len(submesh['ghost_nodes'][p])
-        print 'There are %d ghost nodes on proc %d' %(N, p)
+        print 'There are %d ghost nodes and %d ghost triangles on proc %d'\
+              %(N, M, p)
 
 
Index: anuga_core/source/anuga_parallel/parallel_shallow_water.py
===================================================================
--- anuga_core/source/anuga_parallel/parallel_shallow_water.py	(revision 3927)
+++ anuga_core/source/anuga_parallel/parallel_shallow_water.py	(revision 3928)
@@ -34,6 +34,6 @@
                  full_send_dict=None,
                  ghost_recv_dict=None,
-                 number_of_full_nodes=0,
-                 number_of_full_triangles=0):
+                 number_of_full_nodes=None,
+                 number_of_full_triangles=None):
 
         Domain.__init__(self,
@@ -48,5 +48,5 @@
                         number_of_full_triangles=number_of_full_triangles)
 
-        N = self.number_of_elements
+        N = len(self) # number_of_triangles
 
 #        self.processor = pypar.rank()
Index: anuga_core/source/anuga_parallel/test_parallel_sw_runup.py
===================================================================
--- anuga_core/source/anuga_parallel/test_parallel_sw_runup.py	(revision 3927)
+++ anuga_core/source/anuga_parallel/test_parallel_sw_runup.py	(revision 3928)
@@ -96,15 +96,12 @@
         local_interpolation_points.append(i)
 
-# Hack
+# Hack before we excluded ghosts.
 if numprocs == 2:
     if myid == 0:
         del local_interpolation_points[0]                
         #local_interpolation_points = [1,2,3]
-
 if numprocs == 3:
     if myid == 1:
         del local_interpolation_points[0]
-
-
 if numprocs == 4:
     if myid == 0:
@@ -122,6 +119,5 @@
 
 for t in domain.evolve(yieldstep = 0.1, finaltime = 5.0):
-    #domain.write_time()
-
+    domain.write_time()
     
     # Record time series at known points
