Index: development/pyvolution-1d/dam.py
===================================================================
--- development/pyvolution-1d/dam.py	(revision 3335)
+++ development/pyvolution-1d/dam.py	(revision 3362)
@@ -1,2 +1,3 @@
+
 import os
 from math import sqrt, pi
@@ -84,5 +85,5 @@
 
 def stage(x):
-    y = zeros(len(x))
+    y = zeros(len(x),Float)
     for i in range(len(x)):
         if x[i]<=1000.0:
@@ -123,5 +124,5 @@
         StageQ = domain.quantities['stage'].vertex_values
         y = analytical_sol(X.flat,domain.time)
-        #linePlot(plot1,X,StageQ,X,y)
+        linePlot(plot1,X,StageQ,X,y)
     #raw_input('press_return')
     #pass
Index: development/pyvolution-1d/dam2.py
===================================================================
--- development/pyvolution-1d/dam2.py	(revision 3362)
+++ development/pyvolution-1d/dam2.py	(revision 3362)
@@ -0,0 +1,146 @@
+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 = 0.1     # 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 = 1.7117807
+    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 II -- 0.1m Deep Downstream"
+
+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),Float)
+    for i in range(len(x)):
+        if x[i]<=1000.0:
+            y[i] = 10.0
+        else:
+            y[i] = 0.1
+    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
+raw_input('press enter')
+
+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/dry_dam.py
===================================================================
--- development/pyvolution-1d/dry_dam.py	(revision 3362)
+++ development/pyvolution-1d/dry_dam.py	(revision 3362)
@@ -0,0 +1,130 @@
+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 = 0.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
+
+    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] <= 2*t*sqrt(g*h1) ):
+            u[i] = u3 
+            h[i] = h3 
+        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)
+
+debug = False
+
+print "TEST 1D-SOLUTION III -- DRY BED"
+
+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),Float)
+    for i in range(len(x)):
+        if x[i]<=1000.0:
+            y[i] = 10.0
+        else:
+            y[i] = 0.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 = 30.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/quantity.py
===================================================================
--- development/pyvolution-1d/quantity.py	(revision 3335)
+++ development/pyvolution-1d/quantity.py	(revision 3362)
@@ -513,5 +513,9 @@
         """
 
+        Z = self.gradients
+        #print "gradients 1",Z
+        
         self.compute_gradients()
+        #print "gradients 2", Z
 
         G = self.gradients
@@ -527,10 +531,11 @@
             #vertex coordinates
             x0, x1 = V[k,:]
+
             #Extrapolate
             Qv[k,0] = Qc[k] + G[k]*(x0-x)
             Qv[k,1] = Qc[k] + G[k]*(x1-x)
 
-    def limit(self):
-        """Limit slopes for each volume to eliminate artificial variance
+    """    def limit(self):
+        Limit slopes for each volume to eliminate artificial variance
         introduced by e.g. second order extrapolator
 
@@ -542,5 +547,5 @@
         postcondition:
         vertex values are updated
-        """
+        
         print "in Q.limit"
         from Numeric import zeros, Float
@@ -591,5 +596,68 @@
             for i in range(2):
                 qv[k,i] = qc[k] + phi*dq[i]
-
+    """
+
+    def limit(self):
+        """Limit slopes for each volume to eliminate artificial variance
+        introduced by e.g. second order extrapolator
+
+        This is an unsophisticated limiter as it does not take into
+        account dependencies among quantities.
+
+        precondition:
+        vertex values are estimated from gradient
+        postcondition:
+        vertex values are updated
+        """
+
+        from Numeric import zeros, Float
+
+        N = self.domain.number_of_elements
+        beta = self.beta
+
+        qc = self.centroid_values
+        qv = self.vertex_values
+
+        #Find min and max of this and neighbour's centroid values
+        qmax = self.qmax
+        qmin = self.qmin
+
+        for k in range(N):
+            qmax[k] = qmin[k] = qc[k]
+            for i in range(2):
+                n = self.domain.neighbours[k,i]
+                if n >= 0:
+                    qn = qc[n] #Neighbour's centroid value
+
+                    qmin[k] = min(qmin[k], qn)
+                    qmax[k] = max(qmax[k], qn)
+
+
+        #Diffences between centroids and maxima/minima
+        dqmax = qmax - qc
+        dqmin = qmin - qc
+
+        #Deltas between vertex and centroid values
+        dq = zeros(qv.shape, Float)
+        for i in range(2):
+            dq[:,i] = qv[:,i] - qc
+
+        #Phi limiter
+        for k in range(N):
+
+            #Find the gradient limiter (phi) across vertices
+            phi = 1.0
+            for i in range(2):
+                r = 1.0
+                if (dq[k,i] > 0): r = dqmax[k]/dq[k,i]
+                if (dq[k,i] < 0): r = dqmin[k]/dq[k,i]
+
+                phi = min( min(r*beta, 1), phi )
+
+            #Then update using phi limiter
+            for i in range(2):
+                qv[k,i] = qc[k] + phi*dq[k,i]
+
+        
 
 def newLinePlot(title='Simple Plot'):
Index: development/pyvolution-1d/shallow_water_1d.py
===================================================================
--- development/pyvolution-1d/shallow_water_1d.py	(revision 3335)
+++ development/pyvolution-1d/shallow_water_1d.py	(revision 3362)
@@ -1,3 +1,3 @@
-"""Class Domain -
+u"""Class Domain -
 1D interval domains for finite-volume computations of
 the shallow water wave equation.
@@ -165,8 +165,8 @@
         #self.check_integrity()
 
-        msg = 'Parameter beta_h must be in the interval [0, 1)'
-        assert 0 <= self.beta_h < 1.0, msg
-        msg = 'Parameter beta_w must be in the interval [0, 1)'
-        assert 0 <= self.beta_w < 1.0, msg
+        #msg = 'Parameter beta_h must be in the interval [0, 1)'
+        #assert 0 <= self.beta_h < 1.0, msg
+        #msg = 'Parameter beta_w must be in the interval [0, 1)'
+        #assert 0 <= self.beta_w < 1.0, msg
 
 
@@ -605,5 +605,5 @@
 
     #Remove very thin layers of water
-    #protect_against_infinitesimal_and_negative_heights(domain)
+    protect_against_infinitesimal_and_negative_heights(domain)
 
     #Extrapolate all conserved quantities
@@ -634,5 +634,5 @@
 
     #Take bed elevation into account when water heights are small
-    #balance_deep_and_shallow(domain)
+    balance_deep_and_shallow(domain)
 
     #Compute edge values by interpolation
@@ -662,5 +662,5 @@
 
             #Control momentum
-#            xmomc[k] = ymomc[k] = 0.0
+            #xmomc[k] = ymomc[k] = 0.0
             xmomc[k] = 0.0
 
