Index: development/pyvolution-1d/dam.py
===================================================================
--- development/pyvolution-1d/dam.py	(revision 3335)
+++ development/pyvolution-1d/dam.py	(revision 3335)
@@ -0,0 +1,145 @@
+import os
+from math import sqrt, pi
+from shallow_water_1d import *
+from Numeric import allclose, array, zeros, ones, Float, take, sqrt
+from config import g, epsilon
+
+def analytical_sol(C,t):
+    
+    #t  = 0.0     # time (s)
+    g  = 9.81    # gravity (m/s^2)
+    h1 = 10.0    # depth upstream (m)
+    h0 = 5.0     # depth downstream (m)
+    L = 2000.0   # length of stream/domain (m)
+    n = len(C)    # number of cells
+
+    u = zeros(n,Float)
+    h = zeros(n,Float)
+    x = C-L/2
+    #x = zeros(n,Float)
+    #for i in range(n):
+    #    x[i] = C[i]-1000.0
+
+    # Upstream and downstream boundary conditions are set to the intial water
+    # depth for all time.
+
+    # Calculate Shock Speed
+    h2 = 7.2692044
+    S2 = 2*h2/(h2-h0)*(sqrt(g*h1)-sqrt(g*h2))
+    u2 = S2 - g*h0/(4*S2)*(1+sqrt(1+8*S2*S2/(g*h0)))
+    
+    #t=50 
+    #x = (-L/2:L/2) 
+    for i in range(n):
+        # Calculate Analytical Solution at time t > 0
+        u3 = 2/3*(sqrt(g*h1)+x[i]/t) 
+        h3 = 4/(9*g)*(sqrt(g*h1)-x[i]/(2*t))*(sqrt(g*h1)-x[i]/(2*t)) 
+
+        if ( x[i] <= -t*sqrt(g*h1) ):
+            u[i] = 0.0 
+            h[i] = h1 
+        elif ( x[i] <= t*(u2-sqrt(g*h2)) ):
+            u[i] = u3 
+            h[i] = h3 
+        elif ( x[i] < t*S2 ):
+            u[i] = u2 
+            h[i] = h2 
+        else:
+            u[i] = 0.0 
+            h[i] = h0 
+
+    return h
+
+
+def newLinePlot(title='Simple Plot'):
+    import Gnuplot
+    g = Gnuplot.Gnuplot(persist=1)
+    g.title(title)
+    g('set data style linespoints') 
+    g.xlabel('x')
+    g.ylabel('y')
+    return g
+
+def linePlot(g,x1,y1,x2,y2):
+    import Gnuplot
+    plot1 = Gnuplot.PlotItems.Data(x1.flat,y1.flat,with="linespoints")
+    plot2 = Gnuplot.PlotItems.Data(x2.flat,y2.flat, with="lines 3")
+    g.plot(plot1,plot2)
+    #g.plot(Gnuplot.PlotItems.Data(x1.flat,y1.flat),with="linespoints")
+    #g.plot(Gnuplot.PlotItems.Data(x2.flat,y2.flat), with="lines")
+
+debug = False
+
+print "TEST 1D-SOLUTION I"
+
+L = 2000.0     # Length of channel (m)
+N = 100        # Number of compuational cells
+cell_len = L/N # Origin = 0.0
+
+points = zeros(N+1,Float)
+for i in range(N+1):
+    points[i] = i*cell_len
+
+domain = Domain(points)
+
+def stage(x):
+    y = zeros(len(x))
+    for i in range(len(x)):
+        if x[i]<=1000.0:
+            y[i] = 10.0
+        else:
+            y[i] = 5.0
+    return y
+
+domain.set_quantity('stage', stage)
+
+domain.order = 2
+domain.default_order = 2
+print "domain.order", domain.order
+
+if (debug == True):
+    area = domain.areas
+    for i in range(len(area)):
+        if area != 20:
+            print "Cell Areas are Wrong"
+            
+    L = domain.quantities['stage'].vertex_values
+    print "Initial Stage"
+    print L
+
+domain.set_boundary({'exterior': Reflective_boundary(domain)})
+
+X = domain.vertices
+C = domain.centroids
+plot1 = newLinePlot("Stage")
+
+import time
+t0 = time.time()
+yieldstep = 1.0
+finaltime = 50.0
+for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
+    domain.write_time()
+    if t > 0.0:
+        StageQ = domain.quantities['stage'].vertex_values
+        y = analytical_sol(X.flat,domain.time)
+        #linePlot(plot1,X,StageQ,X,y)
+    #raw_input('press_return')
+    #pass
+
+print 'That took %.2f seconds' %(time.time()-t0)
+
+C = domain.quantities['stage'].centroid_values
+
+if (debug == True):
+    L = domain.quantities['stage'].vertex_values
+    print "Final Stage Vertex Values"
+    print L
+    print "Final Stage Centroid Values"
+    print C
+
+#f = file('test_solution_I.out', 'w')
+#for i in range(N):
+#    f.write(str(C[i]))
+#    f.write("\n")
+#f.close
+
Index: development/pyvolution-1d/domain.py
===================================================================
--- development/pyvolution-1d/domain.py	(revision 3322)
+++ development/pyvolution-1d/domain.py	(revision 3335)
@@ -557,14 +557,7 @@
             Q = self.create_quantity_from_expression(expression)
             kwargs['quantity'] = Q
-
-
-        print "self.vertices1"
-        print self.vertices
         
         #Assign values
         self.quantities[name].set_values(*args, **kwargs)
-
-        print "self.vertices2"
-        print self.vertices
 
     def set_boundary(self, boundary_map):
@@ -663,4 +656,28 @@
         ##assert hasattr(self, 'boundary_objects')
 
+    def write_time(self):
+        print self.timestepping_statistics()
+
+    def timestepping_statistics(self):
+        """Return string with time stepping statistics for printing or logging
+        """
+
+        msg = ''
+        if self.min_timestep == self.max_timestep:
+            msg += 'Time = %.4f, delta t = %.8f, steps=%d (%d)'\
+                   %(self.time, self.min_timestep, self.number_of_steps,
+                     self.number_of_first_order_steps)
+        elif self.min_timestep > self.max_timestep:
+            msg += 'Time = %.4f, steps=%d (%d)'\
+                   %(self.time, self.number_of_steps,
+                     self.number_of_first_order_steps)
+        else:
+            msg += 'Time = %.4f, delta t in [%.8f, %.8f], steps=%d (%d)'\
+                   %(self.time, self.min_timestep,
+                     self.max_timestep, self.number_of_steps,
+                     self.number_of_first_order_steps)
+
+        return msg
+
     def get_name(self):
         return self.filename
@@ -676,5 +693,4 @@
 
 #Main components of evolve
-
     def evolve(self, yieldstep = None, finaltime = None,
                skip_initial_step = False):
@@ -721,8 +737,8 @@
         #update ghosts
         #self.update_ghosts()
-
+    
         #Initial update of vertex and edge values
         self.distribute_to_vertices_and_edges()
-
+        
         #Initial update boundary values
         self.update_boundary()
@@ -861,4 +877,5 @@
                 else:
                     #Try to overcome situation by switching to 1 order
+                    print "changing Order 1"
                     self.order = 1
 
Index: development/pyvolution-1d/quantity.py
===================================================================
--- development/pyvolution-1d/quantity.py	(revision 3322)
+++ development/pyvolution-1d/quantity.py	(revision 3335)
@@ -134,11 +134,7 @@
                   Default is "vertices"
         """
-
-        print "In set_function value"
-        print self.domain.vertices
         
         if location == 'centroids':
-
-            print "In set at centroid"            
+            
             P = self.domain.centroids
             self.set_values(f(P), location)
@@ -146,24 +142,7 @@
             #Vertices
             P = self.domain.get_vertices()
-
-            print "P"
-            print P
-            print P[:,0]
-            print f(P[:,0])
-
-            print P[:,1]
-            print f(P[:,1])
+            
             for i in range(2):
-                self.vertex_values[:,i] = f(P[:,i])
-
-
-            
-
-            print "In set at vertices"
-            print self.vertex_values
-
-        print "Out set_function value"
-        print self.domain.vertices        
-
+                self.vertex_values[:,i] = f(P[:,i])        
 
     def set_array_values(self, values, location='vertices'):
@@ -564,5 +543,5 @@
         vertex values are updated
         """
-
+        print "in Q.limit"
         from Numeric import zeros, Float
 
@@ -706,5 +685,6 @@
     print Q2.centroid_values
     print Qc
-
+    raw_input('press_return')
+    
     g3 = newLinePlot('plot 3')
     linePlot(g3,Xc,Qc)
@@ -720,8 +700,2 @@
     raw_input('press return')
 
-
-
-
-
-
-
Index: development/pyvolution-1d/shallow_water_1d.py
===================================================================
--- development/pyvolution-1d/shallow_water_1d.py	(revision 3322)
+++ development/pyvolution-1d/shallow_water_1d.py	(revision 3335)
@@ -63,6 +63,6 @@
 
         #forcing terms not included in 1d domain ?WHy?
-        #self.forcing_terms.append(gravity)
-        #self.forcing_terms.append(manning_friction)
+        self.forcing_terms.append(gravity)
+        self.forcing_terms.append(manning_friction)
         #print "\nI have Removed forcing terms line 64 1dsw"
 
@@ -174,5 +174,5 @@
         #and or visualisation
         self.distribute_to_vertices_and_edges()
-
+        
         #Initialise real time viz if requested
         #if self.visualise is True and self.time == 0.0:
@@ -432,4 +432,5 @@
     #Loop
     timestep = float(sys.maxint)
+    enter = True
     for k in range(N):
 
@@ -473,23 +474,25 @@
             # flux = edgefluxleft - edgefluxright
             flux -= edgeflux #* domain.edgelengths[k,i]
-
             #Update optimal_timestep
             try:
                 #timestep = min(timestep, 0.5*domain.radii[k]/max_speed)
                 #timestep = 0.01
+            
                 timestep = min(timestep, 0.5*domain.areas[k]/max_speed)
-                if timestep < 0.00001:
-                    #print 'max_speed', max_speed
+                if (timestep < 1e-6) & (enter == True):
+                    #print "domain.order", domain.order
+                    #domain.write_time()
+                    print "cell number", k
+                    print "timestep", timestep
+                    print 'max_speed', max_speed
                     s = domain.quantities['stage']
                     s = s.centroid_values
                     xm = domain.quantities['xmomentum']
                     xm = xm.centroid_values
-                    #print 'h', s[k]
-                    #print 'xm', xm[k]
-                    #print 'u', xm[k]/s[k]
-                    #break
-                #timestep = 0.01
-                #print 'areas', domain.areas[k]
-                #print "timestep", timestep
+                    print 'h', s[k]
+                    print 'xm', xm[k]
+                    print 'u', xm[k]/s[k]
+                    enter = False
+            
             except ZeroDivisionError:
                 pass
@@ -497,5 +500,5 @@
         #Normalise by area and store for when all conserved
         #quantities get updated
-        #flux /= domain.areas[k]
+        flux /= domain.areas[k]
         # ADD ABOVE LINE AGAIN
         Stage.explicit_update[k] = flux[0]
@@ -599,5 +602,5 @@
     """
 
-    from config import optimised_gradient_limiter
+    #from config import optimised_gradient_limiter
 
     #Remove very thin layers of water
