Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py	(revision 5846)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py	(revision 5847)
@@ -199,5 +199,5 @@
         self.CFL = CFL
         self.set_timestepping_method(timestepping_method)
-        
+        self.set_beta(beta_w)
         self.boundary_map = None  # Will be populated by set_boundary        
         
@@ -297,4 +297,20 @@
 
         return self.time + self.starttime
+
+    def set_beta(self,beta):
+        """Set default beta for limiting
+        """
+
+        self.beta = beta
+        for name in self.quantities:
+            #print 'setting beta for quantity ',name
+            Q = self.quantities[name]
+            Q.set_beta(beta)
+
+    def get_beta(self):
+        """Get default beta for limiting
+        """
+
+        return self.beta
 
     def set_default_order(self, n):
@@ -1206,4 +1222,8 @@
         self.update_ghosts()
 
+
+        # Update time
+        self.time += self.timestep
+
         # Update vertex and edge values
         self.distribute_to_vertices_and_edges()
@@ -1212,6 +1232,4 @@
         self.update_boundary()
 
-        # Update time
-        self.time += self.timestep
 
         
@@ -1243,4 +1261,7 @@
         self.update_ghosts()
 
+        # Update time
+        self.time += self.timestep
+
         # Update vertex and edge values
         self.distribute_to_vertices_and_edges()
@@ -1248,7 +1269,4 @@
         # Update boundary values
         self.update_boundary()
-
-        # Update time
-        self.time += self.timestep
 
         #------------------------------------
@@ -1309,4 +1327,7 @@
         self.update_ghosts()
 
+        # Update time
+        self.time += self.timestep
+
         # Update vertex and edge values
         self.distribute_to_vertices_and_edges()
@@ -1315,6 +1336,4 @@
         self.update_boundary()
 
-        # Update time
-        self.time += self.timestep
 
         #------------------------------------
@@ -1339,4 +1358,8 @@
         self.update_ghosts()
 
+
+        # Set substep time
+        self.time = initial_time + self.timestep*0.5
+
         # Update vertex and edge values
         self.distribute_to_vertices_and_edges()
@@ -1345,6 +1368,4 @@
         self.update_boundary()
 
-        # Set substep time
-        self.time = initial_time + self.timestep*0.5
 
         #------------------------------------
@@ -1369,4 +1390,7 @@
         self.update_ghosts()
 
+        # Set new time
+        self.time = initial_time + self.timestep       
+
         # Update vertex and edge values
         self.distribute_to_vertices_and_edges()
@@ -1375,11 +1399,7 @@
         self.update_boundary()
 
-        # Set new time
-        self.time = initial_time + self.timestep       
-        
-
+        
     def evolve_to_end(self, finaltime = 1.0):
-        """Iterate evolve all the way to the end
-        """
+        """Iterate evolve all the way to the end      """
 
         for _ in self.evolve(yieldstep=None, finaltime=finaltime):
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py	(revision 5846)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py	(revision 5847)
@@ -82,5 +82,5 @@
         self.centroid_backup_values = zeros(N, Float)
 
-        self.beta = 1.0
+        self.set_beta(1.0)
 
 
@@ -223,5 +223,21 @@
     #    """
     #    pass
-        
+
+    def set_beta(self,beta):
+        """Set default beta value for limiting
+        """
+
+        if beta < 0.0:
+            print 'WARNING: setting beta < 0.0'
+        if beta > 2.0:
+            print 'WARNING: setting beta > 2.0'
+            
+        self.beta = beta
+
+    def get_beta(self):
+        """Get default beta value for limiting
+        """
+
+        return self.beta
 
     def interpolate(self):
Index: /anuga_core/source/anuga/advection/advection.py
===================================================================
--- /anuga_core/source/anuga/advection/advection.py	(revision 5846)
+++ /anuga_core/source/anuga/advection/advection.py	(revision 5847)
@@ -73,6 +73,7 @@
 
         #Only first is implemented for advection
-        self.default_order = self.order = 1
-
+        self.set_default_order(1)
+        self.set_beta(1.0)
+        
         self.smooth = True
 
@@ -82,4 +83,27 @@
         msg = 'Conserved quantity must be "stage"'
         assert self.conserved_quantities[0] == 'stage', msg
+
+
+    def distribute_to_vertices_and_edges(self):
+        """Extrapolate conserved quantities from centroid to
+        vertices and edge-midpoints for each volume
+
+        Default implementation is straight first order,
+        i.e. constant values throughout each element and
+        no reference to non-conserved quantities.
+        """
+
+        for name in self.conserved_quantities:
+            Q = self.quantities[name]
+            if self._order_ == 1:
+                Q.extrapolate_first_order()
+            elif self._order_ == 2:
+                Q.extrapolate_second_order_and_limit_by_edge()
+                #Q.limit()
+            else:
+                raise 'Unknown order'
+            #Q.interpolate_from_vertices_to_edges()
+
+
 
 
@@ -161,22 +185,22 @@
 
 
-    def evolve(self,
-               yieldstep = None,
-               finaltime = None,
-               duration = None,
-               skip_initial_step = False):
-
-        """Specialisation of basic evolve method from parent class
-        """
-
-        #Call basic machinery from parent class
-        for t in Generic_domain.evolve(self,
-                                       yieldstep=yieldstep,
-                                       finaltime=finaltime,
-                                       duration=duration,
-                                       skip_initial_step=skip_initial_step):
-
-            #Pass control on to outer loop for more specific actions
-            yield(t)
+##     def evolve(self,
+##                yieldstep = None,
+##                finaltime = None,
+##                duration = None,
+##                skip_initial_step = False):
+
+##         """Specialisation of basic evolve method from parent class
+##         """
+
+##         #Call basic machinery from parent class
+##         for t in Generic_domain.evolve(self,
+##                                        yieldstep=yieldstep,
+##                                        finaltime=finaltime,
+##                                        duration=duration,
+##                                        skip_initial_step=skip_initial_step):
+
+##             #Pass control on to outer loop for more specific actions
+##             yield(t)
 
 
Index: /anuga_core/source/anuga/shallow_water/shallow_water_ext.c
===================================================================
--- /anuga_core/source/anuga/shallow_water/shallow_water_ext.c	(revision 5846)
+++ /anuga_core/source/anuga/shallow_water/shallow_water_ext.c	(revision 5847)
@@ -1877,5 +1877,5 @@
     radii             = get_consecutive_array(domain, "radii");    
     areas             = get_consecutive_array(domain, "areas");    
-    tri_full_flag     = get_consecutive_array(domain, "normals");
+    tri_full_flag     = get_consecutive_array(domain, "tri_full_flag");
     already_computed_flux  = get_consecutive_array(domain, "already_computed_flux");
     max_speed_array   = get_consecutive_array(domain, "max_speed");
Index: /anuga_core/source/anuga/shallow_water/test_system.py
===================================================================
--- /anuga_core/source/anuga/shallow_water/test_system.py	(revision 5846)
+++ /anuga_core/source/anuga/shallow_water/test_system.py	(revision 5847)
@@ -150,5 +150,5 @@
         for t in domain.evolve(yieldstep = 5, finaltime = 9.0):
             pass
-            #print domain.write_time()
+            #domain.write_time()
             #print "domain.time", domain.time
 
@@ -157,4 +157,5 @@
         times = fid.variables['time'][:]
         stage = fid.variables['stage'][:]
+        #print stage
         #print "times", times
         #print "fid.starttime", fid.starttime
Index: /anuga_core/source/anuga_parallel/parallel_api.py
===================================================================
--- /anuga_core/source/anuga_parallel/parallel_api.py	(revision 5846)
+++ /anuga_core/source/anuga_parallel/parallel_api.py	(revision 5847)
@@ -43,10 +43,11 @@
     # 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 myid == 0:
+        bdmap = {}
+        for tag in domain.get_boundary_tags():
+            bdmap[tag] = None
+    
+    
+        domain.set_boundary(bdmap)
 
 
Index: /anuga_core/source/anuga_parallel/test_parallel_sw_runup.py
===================================================================
--- /anuga_core/source/anuga_parallel/test_parallel_sw_runup.py	(revision 5846)
+++ /anuga_core/source/anuga_parallel/test_parallel_sw_runup.py	(revision 5847)
@@ -32,13 +32,15 @@
 
 #--------------------------------------------------------------------------
-# Setup computational domain
+# Setup computational domain on processor 0
 #--------------------------------------------------------------------------
-points, vertices, boundary = rectangular_cross(10, 10) # Basic mesh
+if myid == 0 :
+    points, vertices, boundary = rectangular_cross(10, 10) # Basic mesh
 
-#if myid == 0:
-#    print 'points', points
-#    print 'vertices', vertices
+    #print 'points', points
+    #print 'vertices', vertices
     
-domain = Domain(points, vertices, boundary) # Create domain
+    domain = Domain(points, vertices, boundary) # Create domain
+else:
+    domain = None
 
 #--------------------------------------------------------------------------
@@ -49,7 +51,9 @@
     return -x/2                              # linear bed slope
 
-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
+
+if myid == 0:
+    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
 
 
@@ -129,28 +133,30 @@
     stage = domain.get_quantity('stage')
     w = stage.get_values(interpolation_points=interpolation_points)
+
+    print 'P%d:w(%f) = '%(myid,domain.get_time()),w
     
     for i, _ in enumerate(interpolation_points):
         gauge_values[i].append(w[i])
 
-if myid == 0:
-    for i, (x,y) in enumerate(interpolation_points):
+## if myid == 0:
+##     for i, (x,y) in enumerate(interpolation_points):
 
-        try:
-            from pylab import *
-        except:
-            pass
-        else:
-            ion()
-            hold(False)
-            plot(time, gauge_values[i], 'r.-')
-            #time, predicted_gauge_values[i], 'k-')
+##         try:
+##             from pylab import *
+##         except:
+##             pass
+##         else:
+##             ion()
+##             hold(False)
+##             plot(time, gauge_values[i], 'r.-')
+##             #time, predicted_gauge_values[i], 'k-')
 
-            title('Gauge %d (%f,%f)' %(i,x,y))
-            xlabel('time(s)')
-            ylabel('stage (m)')    
-            #legend(('Observed', 'Modelled'), shadow=True, loc='upper left')
-            #savefig('Gauge_%d.png' %i, dpi = 300)
+##             title('Gauge %d (%f,%f)' %(i,x,y))
+##             xlabel('time(s)')
+##             ylabel('stage (m)')    
+##             #legend(('Observed', 'Modelled'), shadow=True, loc='upper left')
+##             #savefig('Gauge_%d.png' %i, dpi = 300)
 
-        raw_input('Next')
+##         raw_input('Next')
         
 
Index: /anuga_core/source/obsolete_code/netherlands-chris.py
===================================================================
--- /anuga_core/source/obsolete_code/netherlands-chris.py	(revision 5846)
+++ /anuga_core/source/obsolete_code/netherlands-chris.py	(revision 5847)
@@ -97,5 +97,5 @@
 N = 60
 
-N = 4
+N = 40
 #N = 140
 #N = 15
@@ -165,4 +165,18 @@
 domain.set_quantity('stage', Constant_height(Z, 0.))
 
+
+visualize = True
+if visualize:
+    from anuga.visualiser import RealtimeVisualiser
+    vis = RealtimeVisualiser(domain)
+    vis.render_quantity_height("elevation", zScale=100, offset = 5.0, dynamic=False)
+    vis.render_quantity_height("stage", zScale =100, dynamic=True, opacity = 0.3)
+    vis.colour_height_quantity('stage', (lambda q:q['stage'], -0.5, 0.5))
+    vis.start()
+
+time.sleep(2.0)
+
+
+
 #Evolve
 import time
@@ -172,5 +186,7 @@
     domain.write_time()
 
-
+    if visualize: vis.update()
+    
+if visualize: vis.evolveFinished()
 
     
Index: /anuga_core/source/pypar-numeric/mandelbrot_example/compile.py
===================================================================
--- /anuga_core/source/pypar-numeric/mandelbrot_example/compile.py	(revision 5846)
+++ /anuga_core/source/pypar-numeric/mandelbrot_example/compile.py	(revision 5847)
@@ -149,6 +149,6 @@
     if os.name == 'posix' and os.uname()[4] == 'x86_64':
       #Extra flags for 64 bit architectures
-      #s += ' -fPIC -m64' #gcc
-      s += ' -fPIC -tp amd64' #pgcc
+      s += ' -fPIC -m64' #gcc
+      #s += ' -fPIC -tp amd64' #pgcc
       
     
Index: /anuga_core/source/pypar-numeric/mandelbrot_example/mandel_sequential.py
===================================================================
--- /anuga_core/source/pypar-numeric/mandelbrot_example/mandel_sequential.py	(revision 5846)
+++ /anuga_core/source/pypar-numeric/mandelbrot_example/mandel_sequential.py	(revision 5847)
@@ -12,5 +12,5 @@
 # User definable parameters
 kmax = 2**15  # Maximal number of iterations (=number of colors)
-M = N = 700   # width = height = N (200, 400, 600, 700 are good)
+M = N = 200   # width = height = N (200, 400, 600, 700 are good)
 
 # Region in complex plane
Index: /anuga_core/source/pypar-numeric/mandelbrot_example/mandelplot.py
===================================================================
--- /anuga_core/source/pypar-numeric/mandelbrot_example/mandelplot.py	(revision 5846)
+++ /anuga_core/source/pypar-numeric/mandelbrot_example/mandelplot.py	(revision 5847)
@@ -18,7 +18,8 @@
     """
     
-    from Tkinter import Frame, Canvas, TOP, NW, PhotoImage
+    #from Tkinter import Frame, Canvas, TOP, NW, PhotoImage
+    from pylab import imshow, show
     from Numeric import transpose, Float
-    from Image import new             # PIL    
+    #from Image import new             # PIL    
     import time
 
@@ -27,6 +28,6 @@
     # User definable parameters
 
-    imtype = 'ppm'      # Image format (e.g 'ppm', 'bmp', 'tiff')
-    filename ='mandel'  # Base filename
+    #imtype = 'ppm'      # Image format (e.g 'ppm', 'bmp', 'tiff')
+    #filename ='mandel'  # Base filename
     
     exponent = 0.998    # Exponent for morphing colors, good with kmax = 2**15
@@ -47,5 +48,5 @@
     A = transpose(A).astype(Float)          
 
-    im = new("RGB", A.shape)
+    #im = new("RGB", A.shape)
 
     
@@ -78,6 +79,8 @@
     
     # Save image to file        
-    im.putdata(L)
-    im.save(filename + '.' + imtype, imtype)
+    imshow(A)
+    show()
+    
+    #im.save(filename + '.' + imtype, imtype)
     print 'Computed plot in %.2f seconds: ' %(time.time()-t0)
 
@@ -88,6 +91,6 @@
     #	sys.exit()
 	
-    import os
-    os.system('xv %s' %(filename + '.' + imtype))
+    #import os
+    #os.system('xv %s' %(filename + '.' + imtype))
 
 
Index: /anuga_validation/okushiri_2005/project.py
===================================================================
--- /anuga_validation/okushiri_2005/project.py	(revision 5846)
+++ /anuga_validation/okushiri_2005/project.py	(revision 5847)
Index: /anuga_validation/okushiri_2005/run_okushiri.py
===================================================================
--- /anuga_validation/okushiri_2005/run_okushiri.py	(revision 5846)
+++ /anuga_validation/okushiri_2005/run_okushiri.py	(revision 5847)
Index: /anuga_validation/okushiri_2005/run_okushiri_parallel.py
===================================================================
--- /anuga_validation/okushiri_2005/run_okushiri_parallel.py	(revision 5846)
+++ /anuga_validation/okushiri_2005/run_okushiri_parallel.py	(revision 5847)
@@ -39,20 +39,26 @@
 #-------------------------
 # Create Domain from mesh
+# on processor 0
 #-------------------------
-domain = Domain(project.mesh_filename, use_cache=use_cache, verbose=True)
-print domain.statistics()
+if myid == 0 :
+    domain = Domain(project.mesh_filename, use_cache=use_cache, verbose=True)
+    print domain.statistics()
 
-
-#-------------------------
-# Initial Conditions
-#-------------------------
-domain.set_quantity('friction', 0.0)
-domain.set_quantity('stage', 0.0)
-domain.set_quantity('elevation',
-                    filename=project.bathymetry_filename,
-                    alpha=0.02,                    
-                    verbose=True,
-                    use_cache=use_cache)
-
+    #-------------------------
+    # Initial Conditions
+    # At present need to do the
+    # fitting of the bathymetry
+    # on a global domain, ie before
+    # distributing the domain
+    #-------------------------
+    domain.set_quantity('friction', 0.0)
+    domain.set_quantity('stage', 0.0)
+    domain.set_quantity('elevation',
+                        filename=project.bathymetry_filename,
+                        alpha=0.02,                    
+                        verbose=True,
+                        use_cache=use_cache)
+else:
+    domain = None
 
 #-------------------------
@@ -86,4 +92,12 @@
 domain.set_boundary({'wave': Bts, 'wall': Br})
 
+# Select triangle containing ch5 for diagnostic output
+# around known gauge
+# This should get triangle id 32833 with centroid (4.5244, 1.1972) on one processor
+try: 
+    triangle_id = domain.get_triangle_containing_point([4.521, 1.196])
+except:
+    triangle_id = None
+
 
 #-------------------------
@@ -94,5 +108,10 @@
 
 for t in domain.evolve(yieldstep = 0.05, finaltime = 22.5):
-    domain.write_time()
+    if myid == 0 :
+        domain.write_time()
+
+    if triangle_id is not None :
+        print domain.timestepping_statistics(track_speeds=False,
+                                             triangle_id=triangle_id)
 
 print 'That took %.2f seconds' %(time.time()-t0)
Index: /anuga_validation/performance_tests/run_profile.py
===================================================================
--- /anuga_validation/performance_tests/run_profile.py	(revision 5846)
+++ /anuga_validation/performance_tests/run_profile.py	(revision 5847)
@@ -57,4 +57,6 @@
 s = 'for t in domain.evolve(yieldstep = 0.02, finaltime = 0.2): domain.write_time()'
 
+#for t in domain.evolve(yieldstep = 0.02, finaltime = 0.2): domain.write_time()
+
 
 import profile, pstats
