Index: /anuga_work/development/anuga_1d/quantity.py
===================================================================
--- /anuga_work/development/anuga_1d/quantity.py	(revision 5857)
+++ /anuga_work/development/anuga_1d/quantity.py	(revision 5858)
@@ -16,10 +16,11 @@
 
 
+
 class Quantity:
 
+    
     def __init__(self, domain, vertex_values=None):
-
-        #from domain import Domain
-        #from domain_order2 import Domain
+        #Initialise Quantity using optional vertex values.
+        
         from domain import Domain
         from Numeric import array, zeros, Float
@@ -78,10 +79,14 @@
         self.beta = domain.beta        
 
-    #Methods for operator overloading
+
     def __len__(self):
+        """
+        Returns number of intervals.
+        """
         return self.centroid_values.shape[0]
     
     def interpolate(self):
-        """Compute interpolated values at centroid
+        """
+        Compute interpolated values at centroid
         Pre-condition: vertex_values have been set
         """
@@ -104,5 +109,5 @@
         In case of location == 'centroid' the dimension values must
         be a list of a Numerical array of length N, N being the number
-        of elements in the mesh. Otherwise it must be of dimension Nx3
+        of elements in the mesh. Otherwise it must be of dimension Nx2
 
         The values will be stored in elements following their
@@ -117,5 +122,5 @@
         """
 
-        if location not in ['vertices', 'centroids']:#, 'edges']:
+        if location not in ['vertices', 'centroids']:
             msg = 'Invalid location: %s, (possible choices vertices, centroids)' %location
             raise msg
@@ -143,4 +148,5 @@
             #Intialise centroid and edge values
             self.interpolate()
+
 
 
@@ -217,5 +223,5 @@
         return X, Compatible list, Numeric array (see below)
         location: Where values are to be stored.
-                  Permissible options are: vertices, edges, centroid
+                  Permissible options are: vertices, centroid
                   and unique vertices. Default is 'vertices'
 
@@ -233,5 +239,4 @@
         from Numeric import take
 
-        #if location not in ['vertices', 'centroids', 'edges', 'unique vertices']:
         if location not in ['vertices', 'centroids', 'unique vertices']:
             msg = 'Invalid location: %s' %location
@@ -247,8 +252,4 @@
                 indices = range(len(self))
             return take(self.centroid_values,indices)
-        #elif location == 'edges':
-        #    if (indices ==  None):
-        #        indices = range(len(self))
-        #    return take(self.edge_values,indices)
         elif location == 'unique vertices':
             if (indices ==  None):
@@ -257,17 +258,17 @@
             #Go through list of unique vertices
             for unique_vert_id in indices:
-                triangles = self.domain.vertexlist[unique_vert_id]
+                cells = self.domain.vertexlist[unique_vert_id]
 
                 #In case there are unused points
-                if triangles is None:
-                    msg = 'Unique vertex not associated with triangles'
+                if cells is None:
+                    msg = 'Unique vertex not associated with cells'
                     raise msg
 
-                # Go through all triangle, vertex pairs
+                # Go through all cells, vertex pairs
                 # Average the values
                 sum = 0
-                for triangle_id, vertex_id in triangles:
-                    sum += self.vertex_values[triangle_id, vertex_id]
-                vert_values.append(sum/len(triangles))
+                for cell_id, vertex_id in cells:
+                    sum += self.vertex_values[cell_id, vertex_id]
+                vert_values.append(sum/len(cells))
             return Numeric.array(vert_values)
         else:
@@ -276,10 +277,6 @@
             return take(self.vertex_values,indices)
 
-    #Method for outputting model results
-    #FIXME: Split up into geometric and numeric stuff.
-    #FIXME: Geometric (X,Y,V) should live in mesh.py
-    #FIXME: STill remember to move XY to mesh
+
     def get_vertex_values(self,
-                          #xy=True,
                           x=True,
                           smooth = None,
@@ -289,9 +286,9 @@
 
         The vertex values are returned as one sequence in the 1D float array A.
-        If requested the coordinates will be returned in 1D arrays X and Y.
+        If requested the coordinates will be returned in 1D arrays X.
 
         The connectivity is represented as an integer array, V, of dimension
-        M x 3, where M is the number of volumes. Each row has three indices
-        into the X, Y, A arrays defining the triangle.
+        M x 2, where M is the number of volumes. Each row has two indices
+        into the X, A arrays defining the element.
 
         if smooth is True, vertex values corresponding to one common
@@ -302,10 +299,10 @@
         If no smoothings is required, vertex coordinates and values will
         be aggregated as a concatenation of values at
-        vertices 0, vertices 1 and vertices 2
+        vertices 0, vertices 1
 
 
         Calling convention
-        if xy is True:
-           X,Y,A,V = get_vertex_values
+        if x is True:
+           X,A,V = get_vertex_values
         else:
            A,V = get_vertex_values
@@ -351,6 +348,4 @@
                 A[k] = reduction(contributions)
 
-
-            #if xy is True:
             if x is True:
                  #X = self.domain.coordinates[:,0].astype(precision)
@@ -380,12 +375,11 @@
 
             #Do vertex coordinates
-            #if xy is True:
             if x is True:
-                C = self.domain.get_vertex_coordinates()
-
-                X = C[:,0:6:2].copy()
-                Y = C[:,1:6:2].copy()
-
-                return X.flat, Y.flat, A, V
+                X = self.domain.get_vertex_coordinates()
+
+                #X = C[:,0:6:2].copy()
+                #Y = C[:,1:6:2].copy()
+
+                return X.flat, A, V
             else:
                 return A, V
@@ -412,13 +406,6 @@
         N = self.centroid_values.shape[0]
 
-        #print "pre update",self.centroid_values
-        #print "timestep",timestep
-        #print min(self.domain.areas)
-        #print "explicit_update", self.explicit_update
-
         #Explicit updates
         self.centroid_values += timestep*self.explicit_update
-
-        #print "post update",self.centroid_values
         
         #Semi implicit updates
@@ -438,4 +425,5 @@
         """
 
+        #print 'compute_gradient'
 
         from Numeric import array, zeros, Float
@@ -457,13 +445,19 @@
                 k0 = k
                 k1 = k+1
+                k2 = k+2
 
                 q0 = Q[k0]
                 q1 = Q[k1]
+                q2 = Q[k2]
 
                 x0 = X[k0] #V0 centroid
                 x1 = X[k1] #V1 centroid
+                x2 = X[k2]
 
                 #Gradient
-                G[k] = (q1 - q0)/(x1 - x0)
+                #G[k] = (q1 - q0)/(x1 - x0)
+                
+                G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0)
+                G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1)
 
             elif k == N-1:
@@ -472,13 +466,28 @@
                 k0 = k
                 k1 = k-1
+                k2 = k-2
 
                 q0 = Q[k0]
                 q1 = Q[k1]
+                q2 = Q[k2]
 
                 x0 = X[k0] #V0 centroid
                 x1 = X[k1] #V1 centroid
+                x2 = X[k2]
 
                 #Gradient
-                G[k] = (q1 - q0)/(x1 - x0)
+                #G[k] = (q1 - q0)/(x1 - x0)
+                
+                G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0)
+                G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1)
+
+##                q0 = Q[k0]
+##                q1 = Q[k1]
+##
+##                x0 = X[k0] #V0 centroid
+##                x1 = X[k1] #V1 centroid
+##
+##                #Gradient
+##                G[k] = (q1 - q0)/(x1 - x0)
 
             else:
@@ -507,4 +516,6 @@
         """
 
+        #print 'compute_minmod_gradients'
+        
         from Numeric import array, zeros, Float,sign
         
@@ -533,13 +544,19 @@
                 k0 = k
                 k1 = k+1
+                k2 = k+2
 
                 q0 = Q[k0]
                 q1 = Q[k1]
+                q2 = Q[k2]
 
                 x0 = X[k0] #V0 centroid
                 x1 = X[k1] #V1 centroid
+                x2 = X[k2]
 
                 #Gradient
-                G[k] = (q1 - q0)/(x1 - x0)
+                #G[k] = (q1 - q0)/(x1 - x0)
+                
+                G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0)
+                G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1)
 
             elif k == N-1:
@@ -548,13 +565,32 @@
                 k0 = k
                 k1 = k-1
+                k2 = k-2
 
                 q0 = Q[k0]
                 q1 = Q[k1]
+                q2 = Q[k2]
 
                 x0 = X[k0] #V0 centroid
                 x1 = X[k1] #V1 centroid
+                x2 = X[k2]
 
                 #Gradient
-                G[k] = (q1 - q0)/(x1 - x0)
+                #G[k] = (q1 - q0)/(x1 - x0)
+                
+                G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0)
+                G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1)
+
+##                #Get data
+##                k0 = k
+##                k1 = k-1
+##
+##                q0 = Q[k0]
+##                q1 = Q[k1]
+##
+##                x0 = X[k0] #V0 centroid
+##                x1 = X[k1] #V1 centroid
+##
+##                #Gradient
+##                G[k] = (q1 - q0)/(x1 - x0)
 
             elif (self.domain.wet_nodes[k,0] == 2) & (self.domain.wet_nodes[k,1] == 2):
@@ -732,4 +768,5 @@
         #print limiter
         
+        #print 'limit_range'
         N = self.domain.number_of_elements
         qc = self.centroid_values
@@ -866,5 +903,4 @@
         # Call correct module function
         # (either from this module or C-extension)
-        #saxpy_centroid_values(self,a,b)
         self.centroid_values[:] = (a*self.centroid_values + b*self.centroid_backup_values).astype('f')
         
@@ -883,22 +919,36 @@
 
 
+##
+##def newLinePlot(title='Simple Plot'):
+##    import Gnuplot
+##    g = Gnuplot.Gnuplot()
+##    g.title(title)
+##    g('set data style linespoints') 
+##    g.xlabel('x')
+##    g.ylabel('y')
+##    return g
+##
+##def linePlot(g,x,y):
+##    import Gnuplot
+##    g.plot(Gnuplot.PlotItems.Data(x.flat,y.flat))
 
 def newLinePlot(title='Simple Plot'):
-    import Gnuplot
-    g = Gnuplot.Gnuplot()
+    import pylab as g
+    g.ion()
+    g.hold(False)
     g.title(title)
-    g('set data style linespoints') 
     g.xlabel('x')
     g.ylabel('y')
-    return g
-
-def linePlot(g,x,y):
-    import Gnuplot
-    g.plot(Gnuplot.PlotItems.Data(x.flat,y.flat))
-
-
-
-
-
+    
+
+def linePlot(x,y):
+    import pylab as g
+    g.plot(x.flat,y.flat)
+
+
+def closePlots():
+    import pylab as g
+    g.close('all')
+    
 if __name__ == "__main__":
     #from domain import Domain
@@ -950,10 +1000,10 @@
     Qc[1,0] = 3
 
-    Q1.beta = 1.0
     Q1.extrapolate_second_order()
-    Q1.limit_minmod()
-
-    g1 = newLinePlot('plot 1')
-    linePlot(g1,Xc,Qc)
+    #Q1.limit_minmod()
+
+    newLinePlot('plots')
+    linePlot(Xc,Qc)
+    raw_input('press return')
 
     points2 = arange(10)
@@ -964,30 +1014,29 @@
     Xc = Q2.domain.vertices
     Qc = Q2.vertex_values
-
-    g2 = newLinePlot('plot 2')
-    linePlot(g2,Xc,Qc)
-
+    linePlot(Xc,Qc)
+    raw_input('press return')
 
     
     Q2.extrapolate_second_order()
-    Q2.limit_minmod()
+    #Q2.limit_minmod()
     Xc = Q2.domain.vertices
     Qc = Q2.vertex_values
-
     print Q2.centroid_values
     print Qc
-    raw_input('press_return')
-    
-    g3 = newLinePlot('plot 3')
-    linePlot(g3,Xc,Qc)
+    linePlot(Xc,Qc)
     raw_input('press return')
 
 
     for i in range(10):
+        import pylab as g
+        g.hold(True)
         fun = FunClass(i/10.0)
-        Q2.set_values(fun,'vertices')
+        Q2.set_values(fun,'centroids')
+        Q2.extrapolate_second_order()
+        #Q2.limit_minmod()
         Qc = Q2.vertex_values
-        linePlot(g3,Xc,Qc)
-
-    raw_input('press return')
-
+        linePlot(Xc,Qc)
+        raw_input('press return')
+
+    raw_input('press return to quit')
+closePlots()
Index: /anuga_work/development/anuga_1d/test_quantity.py
===================================================================
--- /anuga_work/development/anuga_1d/test_quantity.py	(revision 5857)
+++ /anuga_work/development/anuga_1d/test_quantity.py	(revision 5858)
@@ -6,13 +6,9 @@
 
 from quantity import *
-from anuga.config import epsilon
+from config import epsilon
 from Numeric import allclose, array, ones, Float
-
-from anuga.fit_interpolate.fit import fit_to_mesh
-#from anuga.pyvolution.least_squares import fit_to_mesh         
+       
 from domain import Domain
-from anuga.geospatial_data.geospatial_data import Geospatial_data
-from anuga.coordinate_transforms.geo_reference import Geo_reference
-from anuga.utilities.polygon import *
+
 
 #Aux for fit_interpolate.fit example
