Index: /anuga_work/development/anuga_1d/comp_flux_ext.c
===================================================================
--- /anuga_work/development/anuga_1d/comp_flux_ext.c	(revision 5831)
+++ /anuga_work/development/anuga_1d/comp_flux_ext.c	(revision 5832)
@@ -188,10 +188,10 @@
 				// Update timestep based on edge i and possibly neighbour n
 				if (max_speed > epsilon) {
-	                // Original CFL calculation
-					
-					timestep = min(timestep, 0.5*cfl*areas[k]/max_speed); //Here, CFL=1.0 is assumed. ?????????????????????????????????????????????
-					if (n>=0) {
-						timestep = min(timestep, 0.5*cfl*areas[n]/max_speed); //Here, CFL=1.0 is assumed. ?????????????????????????????????????????????
-					}
+				    // Original CFL calculation
+				    
+				    timestep = min(timestep, 0.5*cfl*areas[k]/max_speed); 
+				    if (n>=0) {
+					timestep = min(timestep, 0.5*cfl*areas[n]/max_speed); 
+				    }
 				}
             } // End edge i (and neighbour n)
@@ -263,6 +263,5 @@
   // Call underlying flux computation routine and update 
   // the explicit update arrays 
-  timestep = _compute_fluxes_ext(
-                 cfl,
+  timestep = _compute_fluxes_ext(cfl,
 				 timestep,
 				 epsilon,
@@ -329,5 +328,5 @@
     g                 = get_python_double(domain,"g");
     h0                = get_python_double(domain,"h0");
-	cfl                = get_python_double(domain,"cfl");
+    cfl               = get_python_double(domain,"CFL");
  
     
@@ -358,22 +357,22 @@
     // the explicit update arrays 
     timestep = _compute_fluxes_ext(
-	               cfl,
-				   timestep,
-				   epsilon,
-				   g,
-				   h0,
-				   (long*) neighbours -> data,
-				   (long*) neighbour_vertices -> data,
-				   (double*) normals -> data,
-				   (double*) areas -> data,
-				   (double*) stage_vertex_values -> data,
-				   (double*) xmom_vertex_values -> data,
-				   (double*) bed_vertex_values -> data,
-				   (double*) stage_boundary_values -> data,
-				   (double*) xmom_boundary_values -> data,
-				   (double*) stage_explicit_update -> data,
-				   (double*) xmom_explicit_update -> data,
-				   number_of_elements,
-				   (double*) max_speed_array -> data);
+	cfl,
+	timestep,
+	epsilon,
+	g,
+	h0,
+	(long*) neighbours -> data,
+	(long*) neighbour_vertices -> data,
+	(double*) normals -> data,
+	(double*) areas -> data,
+	(double*) stage_vertex_values -> data,
+	(double*) xmom_vertex_values -> data,
+	(double*) bed_vertex_values -> data,
+	(double*) stage_boundary_values -> data,
+	(double*) xmom_boundary_values -> data,
+	(double*) stage_explicit_update -> data,
+	(double*) xmom_explicit_update -> data,
+	number_of_elements,
+	(double*) max_speed_array -> data);
 
 
Index: /anuga_work/development/anuga_1d/comp_flux_ext_wellbalanced.c
===================================================================
--- /anuga_work/development/anuga_1d/comp_flux_ext_wellbalanced.c	(revision 5831)
+++ /anuga_work/development/anuga_1d/comp_flux_ext_wellbalanced.c	(revision 5832)
@@ -31,5 +31,5 @@
 		int i;
         double z_left, z_right;
-		double ql[2], qr[2], flux_left[2], flux_right[2];
+		double flux_left[2], flux_right[2];
 		double z_star, w_left, h_left, h_left_star, uh_left, soundspeed_left, u_left;
 		double w_right, h_right, h_right_star, uh_right, soundspeed_right, u_right;                          //h_right,
@@ -136,5 +136,5 @@
 		
 		double flux[2], ql[3], qr[3], edgeflux[2];
-		double zl, zr, max_speed, normal;
+		double max_speed, normal;
 		int k, i, ki, n, m, nm=0;
 		
Index: /anuga_work/development/anuga_1d/domain.py
===================================================================
--- /anuga_work/development/anuga_1d/domain.py	(revision 5831)
+++ /anuga_work/development/anuga_1d/domain.py	(revision 5832)
@@ -7,11 +7,16 @@
    Geoscience Australia
 """
+
 from generic_boundary_conditions import *
-#from coordinate_transforms.geo_reference import Geo_reference
+
 
 class Domain:
 
-    def __init__(self, coordinates, boundary = None,
-                 conserved_quantities = None, other_quantities = None,
+    def __init__(self,
+                 coordinates,
+                 boundary = None,
+                 conserved_quantities = None,
+                 evolved_quantities = None,
+                 other_quantities = None,
                  tagged_elements = None):
         """
@@ -23,14 +28,8 @@
         from config import timestepping_method
         from config import CFL
-
-
         
         #Store Points
         self.coordinates = array(coordinates)
 
-# #        if geo_reference is None:
-# #            self. = Geo_reference() #Use defaults 
-# #        else:
-# #            self.geo_reference = geo_reference
 
         #Register number of Elements
@@ -38,6 +37,6 @@
 
         self.beta = 1.0
-        self.limiter = "minmod_kurganov"
-        self.CFL = CFL
+        self.set_limiter("minmod_kurganov")
+        self.set_CFL(CFL)
         self.set_timestepping_method(timestepping_method)
 
@@ -124,4 +123,9 @@
             self.conserved_quantities = conserved_quantities
 
+        if evolved_quantities is None:
+            self.evolved_quantities = self.conserved_quantities
+        else:
+            self.evolved_quantities = evolved_quantities
+
         if other_quantities is None:
             self.other_quantities = []
@@ -133,6 +137,10 @@
         self.quantities = {}
 
+        #print self.conserved_quantities
+        #print self.evolved_quantities
+        
+
         #FIXME: remove later - maybe OK, though....
-        for name in self.conserved_quantities:
+        for name in self.evolved_quantities:
             self.quantities[name] = Quantity(self)
         for name in self.other_quantities:
@@ -160,5 +168,4 @@
         self.number_of_steps = 0
         self.number_of_first_order_steps = 0
-        self.CFL = CFL
 
         #Model time
@@ -169,5 +176,5 @@
         #Checkpointing and storage
         from config import default_datadir
-        self.datadir = default_datadir
+        self.set_datadir(default_datadir)
         self.filename = 'domain'
         self.checkpoint = False
@@ -490,4 +497,37 @@
 
         return q
+
+
+    def get_evolved_quantities(self, vol_id, vertex=None):#, edge=None):
+        """Get evolved quantities at volume vol_id
+
+        If vertex is specified use it as index for vertex values
+        If edge is specified use it as index for edge values
+        If neither are specified use centroid values
+        If both are specified an exeception is raised
+
+        Return value: Vector of length == number_of_evolved quantities
+
+        """
+
+        from Numeric import zeros, Float
+
+        #if not (vertex is None):# or edge is None):
+        #    msg = 'Values for both vertex and edge was specified.'
+        #    msg += 'Only one (or none) is allowed.'
+       #     raise msg
+
+        q = zeros( len(self.evolved_quantities), Float)
+
+        for i, name in enumerate(self.evolved_quantities):
+            Q = self.quantities[name]
+            if vertex is not None:
+                q[i] = Q.vertex_values[vol_id, vertex]
+            #elif edge is not None:
+            #    q[i] = Q.edge_values[vol_id, edge]
+            else:
+                q[i] = Q.centroid_values[vol_id]
+
+        return q
         
 
@@ -677,8 +717,15 @@
     def check_integrity(self):
         #Mesh.check_integrity(self)
+
+        #print self.quantities
+        #print self.conserved_quantities
         
         for quantity in self.conserved_quantities:
             msg = 'Conserved quantities must be a subset of all quantities'
             assert quantity in self.quantities, msg
+
+        for quantity in self.evolved_quantities:
+            msg = 'Evolved quantities must be a subset of all quantities'
+            assert quantity in self.quantities, msg            
 
         # #assert hasattr(self, 'boundary_objects')
@@ -720,5 +767,35 @@
         self.datadir = name
 
-
+    def set_CFL(self, cfl):
+        if cfl > 1.0:
+            print 'WARNING: Setting CFL condition to %f which is greater than 1' % cfl
+        self.CFL = cfl
+
+    def get_CFL(self):
+        return self.CFL
+    
+    def set_filename(self, name):
+        self.filename = name
+
+    def get_filename(self):
+        return self.filename
+
+    def get_limiter(self):
+        return self.limiter
+
+    def set_limiter(self,limiter):
+
+        possible_limiters = \
+        ['pyvolution', 'steve_minmod', 'minmod', 'minmod_kurganov', 'superbee', 'vanleer', 'vanalbada']
+
+        if limiter in possible_limiters:
+            self.limiter = limiter
+            return
+
+        msg = '%s is an incorrect limiter type.\n'% limiter
+        msg += 'Possible types are: '+ ", ".join(["%s" % el for el in possible_limiters])
+        raise Exception, msg
+
+        
     #--------------------------
     # Main components of evolve
@@ -915,5 +992,5 @@
         """
 
-        # Save initial initial conserved quantities values
+        # Save initial conserved quantities values
         self.backup_conserved_quantities()            
 
@@ -1382,5 +1459,5 @@
             q = B.evaluate(vol_id, vertex_id)
 
-            for j, name in enumerate(self.conserved_quantities):
+            for j, name in enumerate(self.evolved_quantities):
                 Q = self.quantities[name]
                 Q.boundary_values[i] = q[j]
Index: /anuga_work/development/anuga_1d/dry_dam_sudi.py
===================================================================
--- /anuga_work/development/anuga_1d/dry_dam_sudi.py	(revision 5831)
+++ /anuga_work/development/anuga_1d/dry_dam_sudi.py	(revision 5832)
@@ -81,8 +81,8 @@
 import time
 
-finaltime = 20.0
-yieldstep = 20.0
+finaltime = 5.0
+yieldstep = finaltime
 L = 2000.0     # Length of channel (m)
-number_of_cells = [810]#,200,500,1000,2000,5000,10000,20000]
+number_of_cells = [1610]#,200,500,1000,2000,5000,10000,20000]
 h_error = zeros(len(number_of_cells),Float)
 uh_error = zeros(len(number_of_cells),Float)
@@ -101,5 +101,5 @@
     domain.set_boundary({'exterior': Reflective_boundary(domain)})
     domain.order = 2
-    domain.set_timestepping_method('rk2')
+    domain.set_timestepping_method('rk3')
     domain.cfl = 1.0
     domain.limiter = "vanleer"
Index: /anuga_work/development/anuga_1d/shallow_water_domain.py
===================================================================
--- /anuga_work/development/anuga_1d/shallow_water_domain.py	(revision 5831)
+++ /anuga_work/development/anuga_1d/shallow_water_domain.py	(revision 5832)
@@ -54,7 +54,10 @@
         conserved_quantities = ['stage', 'xmomentum']
         other_quantities = ['elevation', 'friction']
-        Generic_Domain.__init__(self, coordinates, boundary,
-                                conserved_quantities, other_quantities,
-                                tagged_elements)
+        Generic_Domain.__init__(self,
+                                coordinates = coordinates,
+                                boundary = boundary,
+                                conserved_quantities = conserved_quantities,
+                                other_quantities = other_quantities,
+                                tagged_elements = tagged_elements)
         
         from config import minimum_allowed_height, g, h0
@@ -68,11 +71,4 @@
         #print "\nI have Removed forcing terms line 64 1dsw"
 
-        #Realtime visualisation
-        self.visualiser = None
-        self.visualise  = False
-        self.visualise_color_stage = False
-        self.visualise_stage_range = 1.0
-        self.visualise_timer = True
-        self.visualise_range_z = None
         
         #Stored output
@@ -82,5 +78,5 @@
 
         #Evolve parametrs
-        self.cfl = 1.0
+        self.set_CFL(1.0)
         
         #Reduction operation for get_vertex_values
@@ -89,8 +85,9 @@
         #self.reduction = min  #Looks better near steep slopes
 
-        self.quantities_to_be_stored = ['stage','xmomentum']
+        self.set_quantities_to_be_stored(['stage','xmomentum'])
 
 	self.__doc__ = 'shallow_water_domain'
 
+        self.check_integrity()
 
     def set_quantities_to_be_stored(self, q):
@@ -125,21 +122,6 @@
         
 
-    def initialise_visualiser(self,scale_z=1.0,rect=None):
-        #Realtime visualisation
-        if self.visualiser is None:
-            from realtime_visualisation_new import Visualiser
-            self.visualiser = Visualiser(self,scale_z,rect)
-            self.visualiser.setup['elevation']=True
-            self.visualiser.updating['stage']=True
-        self.visualise = True
-        if self.visualise_color_stage == True:
-            self.visualiser.coloring['stage'] = True
-            self.visualiser.qcolor['stage'] = (0.0, 0.0, 0.8)
-        print 'initialise visualiser'
-        print self.visualiser.setup
-        print self.visualiser.updating
-
     def check_integrity(self):
-        Generic_Domain.check_integrity(self)
+
         #Check that we are solving the shallow water wave equation
 
@@ -149,4 +131,11 @@
         assert self.conserved_quantities[1] == 'xmomentum', msg
 
+        msg = 'First evolved quantity must be "stage"'
+        assert self.evolved_quantities[0] == 'stage', msg
+        msg = 'Second evolved quantity must be "xmomentum"'
+        assert self.evolved_quantities[1] == 'xmomentum', msg        
+
+        Generic_Domain.check_integrity(self)
+        
     def extrapolate_second_order_sw(self):
         #Call correct module function
@@ -584,4 +573,5 @@
 	import sys
 
+        cfl = domain.get_CFL()
 	
 	timestep = float(sys.maxint)
@@ -596,4 +586,6 @@
 	#print 'g=',g
 	#print 'The type of g is',type(g)
+
+        h0 = domain.h0
 	
 	neighbours = domain.neighbours
@@ -652,20 +644,22 @@
         from comp_flux_ext import compute_fluxes_ext
 	
-	domain.flux_timestep = compute_fluxes_ext(timestep,
-                                  epsilon,
-                                  g,
-                                  neighbours,
-                                  neighbour_vertices,
-                                  normals,
-                                  areas, 
-                                  stage_edge_values,
-                                  xmom_edge_values,
-                                  bed_edge_values,
-                                  stage_boundary_values,
-                                  xmom_boundary_values,
-                                  stage_explicit_update,
-                                  xmom_explicit_update,
-                                  number_of_elements,
-                                  max_speed_array)
+	domain.flux_timestep = compute_fluxes_ext(cfl,
+                                                  timestep,
+                                                  epsilon,
+                                                  g,
+                                                  h0,
+                                                  neighbours,
+                                                  neighbour_vertices,
+                                                  normals,
+                                                  areas, 
+                                                  stage_edge_values,
+                                                  xmom_edge_values,
+                                                  bed_edge_values,
+                                                  stage_boundary_values,
+                                                  xmom_boundary_values,
+                                                  stage_explicit_update,
+                                                  xmom_explicit_update,
+                                                  number_of_elements,
+                                                  max_speed_array)
 
 
Index: /anuga_work/development/anuga_1d/shallow_water_domain_suggestion1.py
===================================================================
--- /anuga_work/development/anuga_1d/shallow_water_domain_suggestion1.py	(revision 5831)
+++ /anuga_work/development/anuga_1d/shallow_water_domain_suggestion1.py	(revision 5832)
@@ -53,8 +53,13 @@
 
         conserved_quantities = ['stage', 'xmomentum']
-        other_quantities = ['elevation', 'friction', 'height', 'velocity']
-        Generic_Domain.__init__(self, coordinates, boundary,
-                                conserved_quantities, other_quantities,
-                                tagged_elements)
+        evolved_quantities = ['stage', 'xmomentum', 'elevation', 'height', 'velocity']
+        other_quantities = ['friction']
+        Generic_Domain.__init__(self,
+                                coordinates = coordinates,
+                                boundary = boundary,
+                                conserved_quantities = conserved_quantities,
+                                evolved_quantities = evolved_quantities,
+                                other_quantities = other_quantities,
+                                tagged_elements = tagged_elements)
         
         from config import minimum_allowed_height, g, h0
@@ -68,11 +73,4 @@
         #print "\nI have Removed forcing terms line 64 1dsw"
 
-        #Realtime visualisation
-        self.visualiser = None
-        self.visualise  = False
-        self.visualise_color_stage = False
-        self.visualise_stage_range = 1.0
-        self.visualise_timer = True
-        self.visualise_range_z = None
         
         #Stored output
@@ -81,6 +79,4 @@
         self.smooth = True
 
-        #Evolve parametrs
-        self.cfl = 1.0
         
         #Reduction operation for get_vertex_values
@@ -89,7 +85,9 @@
         #self.reduction = min  #Looks better near steep slopes
 
-        self.quantities_to_be_stored = ['stage','xmomentum']
+        self.set_quantities_to_be_stored(['stage','xmomentum'])
 
 	self.__doc__ = 'shallow_water_domain'
+
+        self.check_integrity()
 
 
@@ -125,21 +123,6 @@
         
 
-    def initialise_visualiser(self,scale_z=1.0,rect=None):
-        #Realtime visualisation
-        if self.visualiser is None:
-            from realtime_visualisation_new import Visualiser
-            self.visualiser = Visualiser(self,scale_z,rect)
-            self.visualiser.setup['elevation']=True
-            self.visualiser.updating['stage']=True
-        self.visualise = True
-        if self.visualise_color_stage == True:
-            self.visualiser.coloring['stage'] = True
-            self.visualiser.qcolor['stage'] = (0.0, 0.0, 0.8)
-        print 'initialise visualiser'
-        print self.visualiser.setup
-        print self.visualiser.updating
-
     def check_integrity(self):
-        Generic_Domain.check_integrity(self)
+
         #Check that we are solving the shallow water wave equation
 
@@ -148,4 +131,17 @@
         msg = 'Second conserved quantity must be "xmomentum"'
         assert self.conserved_quantities[1] == 'xmomentum', msg
+
+        msg = 'First evolved quantity must be "stage"'
+        assert self.evolved_quantities[0] == 'stage', msg
+        msg = 'Second evolved quantity must be "xmomentum"'
+        assert self.evolved_quantities[1] == 'xmomentum', msg
+        msg = 'Third evolved quantity must be "elevation"'
+        assert self.evolved_quantities[2] == 'elevation', msg
+        msg = 'Fourth evolved quantity must be "height"'
+        assert self.evolved_quantities[3] == 'height', msg
+        msg = 'Fifth evolved quantity must be "velocity"'
+        assert self.evolved_quantities[4] == 'velocity', msg
+
+        Generic_Domain.check_integrity(self)
 
     def extrapolate_second_order_sw(self):
@@ -1080,14 +1076,14 @@
 
         #Handy shorthands
-        #self.stage   = domain.quantities['stage'].edge_values
-        #self.xmom    = domain.quantities['xmomentum'].edge_values
-        #self.ymom    = domain.quantities['ymomentum'].edge_values
-        self.normals = domain.normals
-        self.stage   = domain.quantities['stage'].vertex_values
-        self.xmom    = domain.quantities['xmomentum'].vertex_values
+        self.normals  = domain.normals
+        self.stage    = domain.quantities['stage'].vertex_values
+        self.xmom     = domain.quantities['xmomentum'].vertex_values
+        self.bed      = domain.quantities['elevation'].vertex_values
+        self.height   = domain.quantities['height'].vertex_values
+        self.velocity = domain.quantities['velocity'].vertex_values
 
         from Numeric import zeros, Float
         #self.conserved_quantities = zeros(3, Float)
-        self.conserved_quantities = zeros(2, Float)
+        self.evolved_quantities = zeros(5, Float)
 
     def __repr__(self):
@@ -1100,22 +1096,11 @@
         """
 
-        q = self.conserved_quantities
+        q = self.evolved_quantities
         q[0] = self.stage[vol_id, edge_id]
-        q[1] = self.xmom[vol_id, edge_id]
-        #q[2] = self.ymom[vol_id, edge_id]
-        #normal = self.normals[vol_id, 2*edge_id:2*edge_id+2]
-        #normal = self.normals[vol_id, 2*edge_id:2*edge_id+1]
-        normal = self.normals[vol_id,edge_id]
-
-        #r = rotate(q, normal, direction = 1)
-        #r[1] = -r[1]
-        #q = rotate(r, normal, direction = -1)
-        r = q
-        r[1] = normal*r[1]
-        r[1] = -r[1]
-        r[1] = normal*r[1]
-        q = r
-        #For start interval there is no outward momentum so do not need to
-        #reverse direction in this case
+        q[1] = -self.xmom[vol_id, edge_id]
+        q[2] = self.bed[vol_id, edge_id]
+        q[3] = self.height[vol_id, edge_id]
+        q[1] = -self.velocity[vol_id, edge_id]
+
 
         return q
@@ -1127,13 +1112,13 @@
 
 
-    def __init__(self, conserved_quantities=None):
+    def __init__(self, evolved_quantities=None):
         Boundary.__init__(self)
 
-        if conserved_quantities is None:
+        if evolved_quantities is None:
             msg = 'Must specify one value for each conserved quantity'
             raise msg
 
         from Numeric import array, Float
-        self.conserved_quantities=array(conserved_quantities).astype(Float)
+        self.evolved_quantities=array(evolved_quantities).astype(Float)
 
     def __repr__(self):
@@ -1141,5 +1126,5 @@
 
     def evaluate(self, vol_id=None, edge_id=None):
-        return self.conserved_quantities
+        return self.evolved_quantities
 
 
Index: /anuga_work/development/anuga_1d/test_shallow_water.py
===================================================================
--- /anuga_work/development/anuga_1d/test_shallow_water.py	(revision 5831)
+++ /anuga_work/development/anuga_1d/test_shallow_water.py	(revision 5832)
@@ -5,6 +5,9 @@
 
 
-from shallow_water_domain_new import *
-from shallow_water_domain_new import flux_function as domain_flux_function
+#from shallow_water_domain import *
+#from shallow_water_domain import flux_function as domain_flux_function
+
+from shallow_water_domain_suggestion1 import *
+from shallow_water_domain_suggestion1 import flux_function as domain_flux_function
 
 from Numeric import allclose, array, ones, Float, maximum, zeros
@@ -36,10 +39,12 @@
 
         stage_ud, xmom_ud = local_compute_fluxes(domain)
-        
+
+        domain.distribute_to_vertices_and_edges()
         domain.compute_fluxes()
 
-        #print domain.quantities['stage'].explicit_update
-        #print domain.quantities['xmomentum'].explicit_update
-        #print stage_ud
+        print domain.quantities['stage'].explicit_update
+        print domain.quantities['xmomentum'].explicit_update
+        print stage_ud
+        print xmom_ud
 
         assert allclose( domain.quantities['stage'].explicit_update, stage_ud )
@@ -316,5 +321,5 @@
             try:
                 #timestep = min(timestep, 0.5*domain.radii[k]/max_speed)
-                timestep = min(timestep, domain.cfl*0.5*domain.areas[k]/max_speed)
+                timestep = min(timestep, domain.CFL*0.5*domain.areas[k]/max_speed)
             except ZeroDivisionError:
                 pass
