Changeset 7793

Timestamp:

Jun 6, 2010, 10:48:41 PM (14 years ago)

Author:

steve

Message:

Setting up test function for Quantity

Location:

anuga_work/development/anuga_1d

Files:

• domain.py (modified) (2 diffs)
• quantity.py (modified) (5 diffs)
• test_quantity.py (modified) (33 diffs)

anuga_work/development/anuga_1d/domain.py

@@ -112 +112 @@
         self.build_tagged_elements_dictionary(tagged_elements)
         
-        from quantity import Quantity, Conserved_quantity
-        #from quantity_domain import Quantity, Conserved_quantity
+        from quantity import Quantity
+
         
         #List of quantity names entering
@@ -592 +592 @@
             return self.vertices[elem_id,vertex]
 
-    def get_area(self, elem_id):
-        """Return area of element id
-        """
-
-        return self.areas[elem_id]
+    def get_area(self, elem_id=None):
+        """Return total domain area or area of element id
+        """
+
+        if elem_id is None:
+            return sum(self.areas)
+        else:
+            return self.areas[elem_id]
 
     def get_quantity(self, name, location='vertices', indices = None):

anuga_work/development/anuga_1d/quantity.py

@@ -1 @@
-"""Class Quantity - Implements values at each 1d element
+"""
+Class Quantity - Implements values at each 1d element
 
 To create:
@@ -85 +86 @@
         """
         return self.centroid_values.shape[0]
-    
+
+    def __neg__(self):
+        """Negate all values in this quantity giving meaning to the
+        expression -Q where Q is an instance of class Quantity
+        """
+
+        Q = Quantity(self.domain)
+        Q.set_values_from_numeric(-self.vertex_values)
+        return Q
+
+
+
+    def __add__(self, other):
+        """Add to self anything that could populate a quantity
+
+        E.g other can be a constant, an array, a function, another quantity
+        (except for a filename or points, attributes (for now))
+        - see set_values_from_numeric for details
+        """
+
+        Q = Quantity(self.domain)
+        Q.set_values_from_numeric(other)
+
+        result = Quantity(self.domain)
+        result.set_values_from_numeric(self.vertex_values + Q.vertex_values)
+        return result
+
+    def __radd__(self, other):
+        """Handle cases like 7+Q, where Q is an instance of class Quantity
+        """
+
+        return self + other
+
+    def __sub__(self, other):
+        return self + -other            # Invoke self.__neg__()
+
+    def __mul__(self, other):
+        """Multiply self with anything that could populate a quantity
+
+        E.g other can be a constant, an array, a function, another quantity
+        (except for a filename or points, attributes (for now))
+        - see set_values_from_numeric for details
+        """
+
+        if isinstance(other, Quantity):
+            Q = other
+        else:
+            Q = Quantity(self.domain)
+            Q.set_values_from_numeric(other)
+
+        result = Quantity(self.domain)
+
+        # The product of vertex_values, edge_values and centroid_values
+        # are calculated and assigned directly without using
+        # set_values_from_numeric (which calls interpolate). Otherwise
+        # centroid values wouldn't be products from q1 and q2
+        result.vertex_values = self.vertex_values * Q.vertex_values
+        result.centroid_values = self.centroid_values * Q.centroid_values
+
+        return result
+
+    def __rmul__(self, other):
+        """Handle cases like 3*Q, where Q is an instance of class Quantity
+        """
+
+        return self * other
+
+    def __div__(self, other):
+        """Divide self with anything that could populate a quantity
+
+        E.g other can be a constant, an array, a function, another quantity
+        (except for a filename or points, attributes (for now))
+        - see set_values_from_numeric for details
+
+        Zero division is dealt with by adding an epsilon to the divisore
+        FIXME (Ole): Replace this with native INF once we migrate to NumPy
+        """
+
+        if isinstance(other, Quantity):
+            Q = other
+        else:
+            Q = Quantity(self.domain)
+            Q.set_values_from_numeric(other)
+
+        result = Quantity(self.domain)
+
+        # The quotient of vertex_values, edge_values and centroid_values
+        # are calculated and assigned directly without using
+        # set_values_from_numeric (which calls interpolate). Otherwise
+        # centroid values wouldn't be quotient of q1 and q2
+        result.vertex_values = self.vertex_values/(Q.vertex_values + epsilon)
+        result.centroid_values = self.centroid_values/(Q.centroid_values + epsilon)
+
+        return result
+
+    def __rdiv__(self, other):
+        """Handle cases like 3/Q, where Q is an instance of class Quantity
+        """
+
+        return self / other
+
+    def __pow__(self, other):
+        """Raise quantity to (numerical) power
+
+        As with __mul__ vertex values are processed entry by entry
+        while centroid and edge values are re-interpolated.
+
+        Example using __pow__:
+          Q = (Q1**2 + Q2**2)**0.5
+        """
+
+        if isinstance(other, Quantity):
+            Q = other
+        else:
+            Q = Quantity(self.domain)
+            Q.set_values_from_numeric(other)
+
+        result = Quantity(self.domain)
+
+        # The power of vertex_values, edge_values and centroid_values
+        # are calculated and assigned directly without using
+        # set_values_from_numeric (which calls interpolate). Otherwise
+        # centroid values wouldn't be correct
+        result.vertex_values = self.vertex_values ** other
+        result.centroid_values = self.centroid_values ** other
+
+        return result
+
+    def set_values_from_numeric(self, numeric):
+        import Numeric
+        x = Numeric.array([1.0,2.0])
+        y = [1.0,2.0]
+
+        if type(numeric) == type(y):
+            self.set_values_from_array(numeric)
+        elif type(numeric) == type(x):
+            self.set_values_from_array(numeric)
+        elif callable(numeric):
+            self.set_values_from_function(numeric)
+        elif isinstance(numeric, Quantity):
+            self.set_values_from_quantity(numeric)
+        else:   # see if it's coercible to a float (float, int or long, etc)
+            try:
+                numeric = float(numeric)
+            except ValueError:
+                msg = ("Illegal type for variable 'numeric': %s" % type(numeric))
+                raise Exception, msg
+            self.set_values_from_constant(numeric)
+
+    def set_values_from_constant(self,numeric):
+
+        self.vertex_values[:]   = numeric
+        self.centroid_values[:] = numeric
+
+
+    def set_values_from_array(self,numeric):
+
+        self.vertex_values[:]   = numeric
+        self.interpolate()
+
+
+    def set_values_from_quantity(self,quantity):
+
+        self.vertex_values[:]   = quantity.vertex_values
+        self.centroid_values[:] = quantity.centroid_values
+
+    def set_values_from_function(self,function):
+
+        self.vertex_values[:]   = map(function, self.domain.vertices)
+        self.centroid_values[:] = map(function, self.domain.centroids)
+
+
     def interpolate(self):
         """
@@ -218 +390 @@
             assert values.shape[1] == 2, msg
 
-            self.vertex_values = values
+            self.vertex_values[:] = values
 
 
@@ -428 +600 @@
         """
 
-        #print 'compute_gradient'
 
         from Numeric import array, zeros, Float
@@ -908 +1079 @@
         self.centroid_values[:] = (a*self.centroid_values + b*self.centroid_backup_values).astype('f')
         
-class Conserved_quantity(Quantity):
-    """Class conserved quantity adds to Quantity:
-
-    storage and method for updating, and
-    methods for extrapolation from centropid to vertices inluding
-    gradients and limiters
-    """
-
-    def __init__(self, domain, vertex_values=None):
-        Quantity.__init__(self, domain, vertex_values)
-
-        print "Use Quantity instead of Conserved_quantity"
-
-
-##
-##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'):

anuga_work/development/anuga_1d/test_quantity.py

@@ +1 @@
+import quantity
 #!/usr/bin/env python
 
 import unittest
 from math import sqrt, pi
-import tempfile
+
 
 from quantity import *
-from config import epsilon
-from Numeric import allclose, array, ones, Float
-       
-from domain import Domain
-
-
-#Aux for fit_interpolate.fit example
-def linear_function(point):
-    point = array(point)
-    return point[:,0]+point[:,1]
+from domain import *
+
+
+from Numeric import allclose, array, ones, Float, maximum, zeros
 
 
 class Test_Quantity(unittest.TestCase):
     def setUp(self):
-        from domain import Domain
-
-        a = [0.0, 0.0]
-        b = [0.0, 2.0]
-        c = [2.0, 0.0]
-        d = [0.0, 4.0]
-        e = [2.0, 2.0]
-        f = [4.0, 0.0]
-
-        points = [a, b, c, d, e, f]
-
-        #bac, bce, ecf, dbe
-        elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
-
-        self.mesh1 = Domain(points[:3], [elements[0]])
-        self.mesh1.check_integrity()
-
-        self.mesh4 = Domain(points, elements)
-        self.mesh4.check_integrity()
-
-        # UTM round Onslow
-        a = [240000, 7620000]
-        b = [240000, 7680000]
-        c = [300000, 7620000]
-
-        points = [a, b, c]
-        elements = [[0,2,1]]
-        
-        self.mesh_onslow = Domain(points, elements)
-        self.mesh_onslow.check_integrity()
-        
+        self.points3        = [0.0, 1.0, 2.0, 3.0]
+        self.vertex_values3 = [[1.0,2.0],[4.0,5.0],[-1.0,2.0]]
+        self.domain3        = Domain(self.points3)
+
+
+        
+        self.points4          = [0.0, 1.0, 2.5, 3.0, 5.0]
+        self.vertex_values4   = [[1.0,2.0],[4.0,5.0],[-1.0,2.0],[3.0,6.0]]
+        self.centroid_values4 = [1.5, 4.5, 0.5, 4.5]
+        self.boundary4        = {(0, 0): 'left', (3, 1): 'right'}
+        self.domain4          = Domain(self.points4,self.boundary4)
+
+        self.points10 = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]
+        self.domain10 = Domain(self.points10)
+
+        self.points6 = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0]
+        self.domain6 = Domain(self.points6)
+
+
     def tearDown(self):
         pass
@@ -56 +39 @@
 
 
+    def test_creat_with_boundary(self):
+
+        assert self.domain4.boundary  == {(0, 0): 'left', (3, 1): 'right'}
+        
     def test_creation(self):
 
-        quantity = Quantity(self.mesh1, [[1,2,3]])
-        assert allclose(quantity.vertex_values, [[1.,2.,3.]])
-
+        quantity = Quantity(self.domain3)
+        assert allclose(quantity.vertex_values, [[0.0,0.0],[0.0,0.0],[0.0,0.0]])
+
+        
         try:
             quantity = Quantity()
@@ -74 +62 @@
             pass
         except:
-            raise 'Should have raised "mising mesh object" error'
+            raise 'Should have raised "mising domain object" error'
 
 
     def test_creation_zeros(self):
 
-        quantity = Quantity(self.mesh1)
-        assert allclose(quantity.vertex_values, [[0.,0.,0.]])
-
-
-        quantity = Quantity(self.mesh4)
-        assert allclose(quantity.vertex_values, [[0.,0.,0.], [0.,0.,0.],
-                                                 [0.,0.,0.], [0.,0.,0.]])
+        quantity = Quantity(self.domain3)
+        assert allclose(quantity.centroid_values, [[0.,0.,0.]])
+
+
+        quantity = Quantity(self.domain4)
+        assert allclose(quantity.vertex_values, [[0.,0.], [0.,0.],
+                                                 [0.,0.], [0.,0.]])
 
 
     def test_interpolation(self):
-        quantity = Quantity(self.mesh1, [[1,2,3]])
-        assert allclose(quantity.centroid_values, [2.0]) #Centroid
-
-        assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5]])
+        quantity = Quantity(self.domain4, self.vertex_values4)
+        assert allclose(quantity.centroid_values, self.centroid_values4) #Centroid
+
 
 
     def test_interpolation2(self):
-        quantity = Quantity(self.mesh4,
-                            [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
-        assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
-
+        quantity = Quantity(self.domain4, self.vertex_values4)
+        assert allclose(quantity.centroid_values, self.centroid_values4) #Centroid
 
         quantity.extrapolate_second_order()
 
         #print quantity.vertex_values
-        assert allclose(quantity.vertex_values, [[3.5, -1.0, 3.5],
-                                                 [3.+2./3, 6.+2./3, 4.+2./3],
-                                                 [4.6, 3.4, 1.],
-                                                 [-5.0, 1.0, 4.0]])
-
-        #print quantity.edge_values
-        assert allclose(quantity.edge_values, [[1.25, 3.5, 1.25],
-                                               [5. + 2/3.0, 4.0 + 1.0/6, 5.0 + 1.0/6],
-                                               [2.2, 2.8, 4.0],
-                                               [2.5, -0.5, -2.0]])
-
-
-        #assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
-        #                                       [5., 5., 5.],
-        #                                       [4.5, 4.5, 0.],
-        #                                       [3.0, -1.5, -1.5]])
-
-    def test_get_extrema_1(self):
-        quantity = Quantity(self.mesh4,
-                                      [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
-        assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroids
-
-        v = quantity.get_maximum_value()
-        assert v == 5
-
-        v = quantity.get_minimum_value()
-        assert v == 0        
-
-        i = quantity.get_maximum_index()
-        assert i == 1
-
-        i = quantity.get_minimum_index()
-        assert i == 3        
-        
-        x,y = quantity.get_maximum_location()
-        xref, yref = 4.0/3, 4.0/3
-        assert x == xref
-        assert y == yref
-
-        v = quantity.get_values(interpolation_points = [[x,y]])
-        assert allclose(v, 5)
-
-
-        x,y = quantity.get_minimum_location()
-        v = quantity.get_values(interpolation_points = [[x,y]])
-        assert allclose(v, 0)
-
-
-    def test_get_maximum_2(self):
-
-        a = [0.0, 0.0]
-        b = [0.0, 2.0]
-        c = [2.0,0.0]
-        d = [0.0, 4.0]
-        e = [2.0, 2.0]
-        f = [4.0,0.0]
-
-        points = [a, b, c, d, e, f]
-        #bac, bce, ecf, dbe
-        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
-
-        domain = Domain(points, vertices)
-
-        quantity = Quantity(domain)
-        quantity.set_values(lambda x, y: x+2*y) #2 4 4 6
-        
-        v = quantity.get_maximum_value()
-        assert v == 6
-
-        v = quantity.get_minimum_value()
-        assert v == 2        
-
-        i = quantity.get_maximum_index()
-        assert i == 3
-
-        i = quantity.get_minimum_index()
-        assert i == 0        
-        
-        x,y = quantity.get_maximum_location()
-        xref, yref = 2.0/3, 8.0/3
-        assert x == xref
-        assert y == yref
-
-        v = quantity.get_values(interpolation_points = [[x,y]])
-        assert allclose(v, 6)
-
-        x,y = quantity.get_minimum_location()        
-        v = quantity.get_values(interpolation_points = [[x,y]])
-        assert allclose(v, 2)
-
-        #Multiple locations for maximum -
-        #Test that the algorithm picks the first occurrence        
-        v = quantity.get_maximum_value(indices=[0,1,2])
-        assert allclose(v, 4)
-
-        i = quantity.get_maximum_index(indices=[0,1,2])
-        assert i == 1
-        
-        x,y = quantity.get_maximum_location(indices=[0,1,2])
-        xref, yref = 4.0/3, 4.0/3
-        assert x == xref
-        assert y == yref
-
-        v = quantity.get_values(interpolation_points = [[x,y]])
-        assert allclose(v, 4)        
-
-        # More test of indices......
-        v = quantity.get_maximum_value(indices=[2,3])
-        assert allclose(v, 6)
-
-        i = quantity.get_maximum_index(indices=[2,3])
-        assert i == 3
-        
-        x,y = quantity.get_maximum_location(indices=[2,3])
-        xref, yref = 2.0/3, 8.0/3
-        assert x == xref
-        assert y == yref
-
-        v = quantity.get_values(interpolation_points = [[x,y]])
-        assert allclose(v, 6)        
-
-        
-
+        assert allclose(quantity.vertex_values,[[ 0.3,  2.7],
+                                                [ 4.5,  4.5],
+                                                [ 0.5,  0.5],
+                                                [ 1.3,  7.7]])
+
+ 
+
+ 
     def test_boundary_allocation(self):
-        quantity = Quantity(self.mesh4,
-                            [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
-
-        assert quantity.boundary_values.shape[0] == len(self.mesh4.boundary)
+        quantity = Quantity(self.domain4,
+                            [[1,2], [5,5], [0,9], [-6, 3]])
+
+
+        assert quantity.boundary_values.shape[0] == len(self.domain4.boundary)
 
 
     def test_set_values(self):
-        quantity = Quantity(self.mesh4)
+        quantity = Quantity(self.domain4)
 
 
@@ -281 +152 @@
 
     def test_set_values_const(self):
-        quantity = Quantity(self.mesh4)
+        quantity = Quantity(self.domain4)
 
         quantity.set_values(1.0, location = 'vertices')
@@ -299 +170 @@
 
     def test_set_values_func(self):
-        quantity = Quantity(self.mesh4)
+        quantity = Quantity(self.domain4)
 
         def f(x, y):
@@ -320 +191 @@
 
     def test_integral(self):
-        quantity = Quantity(self.mesh4)
+        quantity = Quantity(self.domain4)
 
         #Try constants first
@@ -327 +198 @@
         #print 'Q', quantity.get_integral()
 
-        assert allclose(quantity.get_integral(), self.mesh4.get_area() * const)
+        assert allclose(quantity.get_integral(), self.domain4.get_area() * const)
 
         #Try with a linear function
-        def f(x, y):
-            return x+y
+        def f(x):
+            return x
 
         quantity.set_values(f, location = 'vertices')
 
-
-        ref_integral = (4.0/3 + 8.0/3 + 10.0/3 + 10.0/3) * 2
+ 
+        assert allclose (quantity.centroid_values,
+        [ 0.5,   1.75,  2.75,  4.  ])
+
+        assert allclose (quantity.vertex_values, [[ 0.,   1. ],
+                                                  [ 1.,   2.5],
+                                                  [ 2.5,  3. ],
+                                                  [ 3.,   5. ]])
+
+
+        ref_integral = 0.5 + 1.5*1.75 + 0.5*2.75 + 2.0*4.0
 
         assert allclose (quantity.get_integral(), ref_integral)
@@ -343 +223 @@
 
     def test_set_vertex_values(self):
-        quantity = Quantity(self.mesh4)
+        quantity = Quantity(self.domain4)
         quantity.set_vertex_values([0,1,2,3,4,5])
 
@@ -357 +237 @@
 
     def test_set_vertex_values_subset(self):
-        quantity = Quantity(self.mesh4)
+        quantity = Quantity(self.domain4)
         quantity.set_vertex_values([0,1,2,3,4,5])
         quantity.set_vertex_values([0,20,30,50], indices = [0,2,3,5])
@@ -366 +246 @@
 
     def test_set_vertex_values_using_general_interface(self):
-        quantity = Quantity(self.mesh4)
+        quantity = Quantity(self.domain4)
 
 
@@ -391 +271 @@
         """
         
-        quantity = Quantity(self.mesh4)
+        quantity = Quantity(self.domain4)
 
 
@@ -462 +342 @@
         """
         
-        quantity = Quantity(self.mesh4)
+        quantity = Quantity(self.domain4)
         G = Geo_reference(56, 10, 100)
         quantity.domain.geo_reference = G
@@ -507 +387 @@
 
 
-    def test_set_values_using_fit(self):
-
-
-        quantity = Quantity(self.mesh4)
-
-        #Get (enough) datapoints
-        data_points = [[ 0.66666667, 0.66666667],
-                       [ 1.33333333, 1.33333333],
-                       [ 2.66666667, 0.66666667],
-                       [ 0.66666667, 2.66666667],
-                       [ 0.0, 1.0],
-                       [ 0.0, 3.0],
-                       [ 1.0, 0.0],
-                       [ 1.0, 1.0],
-                       [ 1.0, 2.0],
-                       [ 1.0, 3.0],
-                       [ 2.0, 1.0],
-                       [ 3.0, 0.0],
-                       [ 3.0, 1.0]]
-
-        z = linear_function(data_points)
-
-        #Use built-in fit_interpolate.fit
-        quantity.set_values( Geospatial_data(data_points, z), alpha = 0 )
-        #quantity.set_values(points = data_points, values = z, alpha = 0)
-
-
-        answer = linear_function(quantity.domain.get_vertex_coordinates())
-        #print quantity.vertex_values, answer
-        assert allclose(quantity.vertex_values.flat, answer)
-
-
-        #Now try by setting the same values directly
-        vertex_attributes = fit_to_mesh(data_points,
-                                        quantity.domain.get_nodes(),
-                                        quantity.domain.triangles, #FIXME
-                                        point_attributes=z,
-                                        alpha = 0,
-                                        verbose=False)
-
-        #print vertex_attributes
-        quantity.set_values(vertex_attributes)
-        assert allclose(quantity.vertex_values.flat, answer)
-
-
-
-
-
-    def test_test_set_values_using_fit_w_geo(self):
-
-
-        #Mesh
-        vertex_coordinates = [[0.76, 0.76],
-                              [0.76, 5.76],
-                              [5.76, 0.76]]
-        triangles = [[0,2,1]]
-
-        mesh_georef = Geo_reference(56,-0.76,-0.76)
-        mesh1 = Domain(vertex_coordinates, triangles,
-                       geo_reference = mesh_georef)
-        mesh1.check_integrity()
-
-        #Quantity
-        quantity = Quantity(mesh1)
-
-        #Data
-        data_points = [[ 201.0, 401.0],
-                       [ 201.0, 403.0],
-                       [ 203.0, 401.0]]
-
-        z = [2, 4, 4]
-
-        data_georef = Geo_reference(56,-200,-400)
-
-
-        #Reference
-        ref = fit_to_mesh(data_points, vertex_coordinates, triangles,
-                          point_attributes=z,
-                          data_origin = data_georef.get_origin(),
-                          mesh_origin = mesh_georef.get_origin(),
-                          alpha = 0)
-
-        assert allclose( ref, [0,5,5] )
-
-
-        #Test set_values
-
-        quantity.set_values( Geospatial_data(data_points, z, data_georef), alpha = 0 )
-
-        #quantity.set_values(points = data_points,
-        #                    values = z,
-        #                    data_georef = data_georef,
-        #                    alpha = 0)
-
-
-        #quantity.set_values(points = data_points,
-        #                    values = z,
-        #                    data_georef = data_georef,
-        #                    alpha = 0)
-        assert allclose(quantity.vertex_values.flat, ref)
-
-
-
-        #Test set_values using geospatial data object
-        quantity.vertex_values[:] = 0.0
-
-        geo = Geospatial_data(data_points, z, data_georef)
-
-
-        quantity.set_values(geospatial_data = geo, alpha = 0)
-        assert allclose(quantity.vertex_values.flat, ref)
-
-
-
-    def test_set_values_from_file1(self):
-        quantity = Quantity(self.mesh4)
-
-        #Get (enough) datapoints
-        data_points = [[ 0.66666667, 0.66666667],
-                       [ 1.33333333, 1.33333333],
-                       [ 2.66666667, 0.66666667],
-                       [ 0.66666667, 2.66666667],
-                       [ 0.0, 1.0],
-                       [ 0.0, 3.0],
-                       [ 1.0, 0.0],
-                       [ 1.0, 1.0],
-                       [ 1.0, 2.0],
-                       [ 1.0, 3.0],
-                       [ 2.0, 1.0],
-                       [ 3.0, 0.0],
-                       [ 3.0, 1.0]]
-
-        data_geo_spatial = Geospatial_data(data_points,
-                         geo_reference = Geo_reference(56, 0, 0))
-        data_points_absolute = data_geo_spatial.get_data_points(absolute=True)
-        attributes = linear_function(data_points_absolute)
-        att = 'spam_and_eggs'
-        
-        #Create .txt file
-        ptsfile = tempfile.mktemp(".txt")
-        file = open(ptsfile,"w")
-        file.write(" x,y," + att + " \n")
-        for data_point, attribute in map(None, data_points_absolute
-                                         ,attributes):
-            row = str(data_point[0]) + ',' + str(data_point[1]) \
-                  + ',' + str(attribute)
-            file.write(row + "\n")
-        file.close()
-
-
-        #Check that values can be set from file
-        quantity.set_values(filename = ptsfile,
-                            attribute_name = att, alpha = 0)
-        answer = linear_function(quantity.domain.get_vertex_coordinates())
-
-        #print quantity.vertex_values.flat
-        #print answer
-
-
-        assert allclose(quantity.vertex_values.flat, answer)
-
-
-        #Check that values can be set from file using default attribute
-        quantity.set_values(filename = ptsfile, alpha = 0)
-        assert allclose(quantity.vertex_values.flat, answer)
-
-        #Cleanup
-        import os
-        os.remove(ptsfile)
-
-
-
-    def Xtest_set_values_from_file_using_polygon(self):
-        """test_set_values_from_file_using_polygon(self):
-        
-        Test that polygon restriction works for general points data
-        """
-        
-        quantity = Quantity(self.mesh4)
-
-        #Get (enough) datapoints
-        data_points = [[ 0.66666667, 0.66666667],
-                       [ 1.33333333, 1.33333333],
-                       [ 2.66666667, 0.66666667],
-                       [ 0.66666667, 2.66666667],
-                       [ 0.0, 1.0],
-                       [ 0.0, 3.0],
-                       [ 1.0, 0.0],
-                       [ 1.0, 1.0],
-                       [ 1.0, 2.0],
-                       [ 1.0, 3.0],
-                       [ 2.0, 1.0],
-                       [ 3.0, 0.0],
-                       [ 3.0, 1.0]]
-
-        data_geo_spatial = Geospatial_data(data_points,
-                         geo_reference = Geo_reference(56, 0, 0))
-        data_points_absolute = data_geo_spatial.get_data_points(absolute=True)
-        attributes = linear_function(data_points_absolute)
-        att = 'spam_and_eggs'
-        
-        #Create .txt file
-        ptsfile = tempfile.mktemp(".txt")
-        file = open(ptsfile,"w")
-        file.write(" x,y," + att + " \n")
-        for data_point, attribute in map(None, data_points_absolute
-                                         ,attributes):
-            row = str(data_point[0]) + ',' + str(data_point[1]) \
-                  + ',' + str(attribute)
-            file.write(row + "\n")
-        file.close()
-
-        # Create restricting polygon (containing node #4 (2,2) and 
-        # centroid of triangle #1 (bce)
-        polygon = [[1.0, 1.0], [4.0, 1.0],
-                   [4.0, 4.0], [1.0, 4.0]]
-
-        #print self.mesh4.nodes
-        #print inside_polygon(self.mesh4.nodes, polygon)
-        assert allclose(inside_polygon(self.mesh4.nodes, polygon), 4)
-
-        #print quantity.domain.get_vertex_coordinates()
-        #print quantity.domain.get_nodes()        
-        
-        # Check that values can be set from file
-        quantity.set_values(filename=ptsfile,
-                            polygon=polygon,
-                            location='unique vertices',
-                            alpha=0)
-
-        # Get indices for vertex coordinates in polygon
-        indices = inside_polygon(quantity.domain.get_vertex_coordinates(), 
-                                 polygon)
-        points = take(quantity.domain.get_vertex_coordinates(), indices)
-        
-        answer = linear_function(points)
-
-        #print quantity.vertex_values.flat
-        #print answer
-
-        # Check vertices in polygon have been set
-        assert allclose(take(quantity.vertex_values.flat, indices),
-                             answer)
-
-        # Check vertices outside polygon are zero
-        indices = outside_polygon(quantity.domain.get_vertex_coordinates(), 
-                                  polygon)        
-        assert allclose(take(quantity.vertex_values.flat, indices),
-                             0.0)        
-
-        #Cleanup
-        import os
-        os.remove(ptsfile)
-
-
-        
-
-    def test_cache_test_set_values_from_file(self):
-        # FIXME (Ole): What is this about?
-        # I don't think it checks anything new
-        quantity = Quantity(self.mesh4)
-
-        #Get (enough) datapoints
-        data_points = [[ 0.66666667, 0.66666667],
-                       [ 1.33333333, 1.33333333],
-                       [ 2.66666667, 0.66666667],
-                       [ 0.66666667, 2.66666667],
-                       [ 0.0, 1.0],
-                       [ 0.0, 3.0],
-                       [ 1.0, 0.0],
-                       [ 1.0, 1.0],
-                       [ 1.0, 2.0],
-                       [ 1.0, 3.0],
-                       [ 2.0, 1.0],
-                       [ 3.0, 0.0],
-                       [ 3.0, 1.0]]
-
-        georef = Geo_reference(56, 0, 0)
-        data_geo_spatial = Geospatial_data(data_points,
-                                           geo_reference=georef)
-                                           
-        data_points_absolute = data_geo_spatial.get_data_points(absolute=True)
-        attributes = linear_function(data_points_absolute)
-        att = 'spam_and_eggs'
-        
-        # Create .txt file
-        ptsfile = tempfile.mktemp(".txt")
-        file = open(ptsfile,"w")
-        file.write(" x,y," + att + " \n")
-        for data_point, attribute in map(None, data_points_absolute
-                                         ,attributes):
-            row = str(data_point[0]) + ',' + str(data_point[1]) \
-                  + ',' + str(attribute)
-            file.write(row + "\n")
-        file.close()
-
-
-        # Check that values can be set from file
-        quantity.set_values(filename=ptsfile,
-                            attribute_name=att, 
-                            alpha=0, 
-                            use_cache=True,
-                            verbose=False)
-        answer = linear_function(quantity.domain.get_vertex_coordinates())
-        assert allclose(quantity.vertex_values.flat, answer)
-
-
-        # Check that values can be set from file using default attribute
-        quantity.set_values(filename=ptsfile, 
-                            alpha=0)
-        assert allclose(quantity.vertex_values.flat, answer)
-
-        # Check cache
-        quantity.set_values(filename=ptsfile,
-                            attribute_name=att, 
-                            alpha=0, 
-                            use_cache=True,
-                            verbose=False)
-        
-        
-        #Cleanup
-        import os
-        os.remove(ptsfile)
-
-    def test_set_values_from_lat_long(self):
-        quantity = Quantity(self.mesh_onslow)
-
-        #Get (enough) datapoints
-        data_points = [[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
-                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6]]
-
-        data_geo_spatial = Geospatial_data(data_points,
-                                           points_are_lats_longs=True)
-        points_UTM = data_geo_spatial.get_data_points(absolute=True)
-        attributes = linear_function(points_UTM)
-        att = 'elevation'
-        
-        #Create .txt file
-        txt_file = tempfile.mktemp(".txt")
-        file = open(txt_file,"w")
-        file.write(" lat,long," + att + " \n")
-        for data_point, attribute in map(None, data_points, attributes):
-            row = str(data_point[0]) + ',' + str(data_point[1]) \
-                  + ',' + str(attribute)
-            #print "row", row 
-            file.write(row + "\n")
-        file.close()
-
-
-        #Check that values can be set from file
-        quantity.set_values(filename=txt_file,
-                            attribute_name=att, 
-                            alpha=0)
-        answer = linear_function(quantity.domain.get_vertex_coordinates())
-
-        #print "quantity.vertex_values.flat", quantity.vertex_values.flat
-        #print "answer",answer
-
-        assert allclose(quantity.vertex_values.flat, answer)
-
-
-        #Check that values can be set from file using default attribute
-        quantity.set_values(filename=txt_file, alpha=0)
-        assert allclose(quantity.vertex_values.flat, answer)
-
-        #Cleanup
-        import os
-        os.remove(txt_file)
-         
-    def test_set_values_from_lat_long(self):
-        quantity = Quantity(self.mesh_onslow)
-
-        #Get (enough) datapoints
-        data_points = [[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
-                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6]]
-
-        data_geo_spatial = Geospatial_data(data_points,
-                                           points_are_lats_longs=True)
-        points_UTM = data_geo_spatial.get_data_points(absolute=True)
-        attributes = linear_function(points_UTM)
-        att = 'elevation'
-        
-        #Create .txt file
-        txt_file = tempfile.mktemp(".txt")
-        file = open(txt_file,"w")
-        file.write(" lat,long," + att + " \n")
-        for data_point, attribute in map(None, data_points, attributes):
-            row = str(data_point[0]) + ',' + str(data_point[1]) \
-                  + ',' + str(attribute)
-            #print "row", row 
-            file.write(row + "\n")
-        file.close()
-
-
-        #Check that values can be set from file
-        quantity.set_values(filename=txt_file,
-                            attribute_name=att, alpha=0)
-        answer = linear_function(quantity.domain.get_vertex_coordinates())
-
-        #print "quantity.vertex_values.flat", quantity.vertex_values.flat
-        #print "answer",answer
-
-        assert allclose(quantity.vertex_values.flat, answer)
-
-
-        #Check that values can be set from file using default attribute
-        quantity.set_values(filename=txt_file, alpha=0)
-        assert allclose(quantity.vertex_values.flat, answer)
-
-        #Cleanup
-        import os
-        os.remove(txt_file)
-        
-    def test_set_values_from_UTM_pts(self):
-        quantity = Quantity(self.mesh_onslow)
-
-        #Get (enough) datapoints
-        data_points = [[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
-                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6]]
-
-        data_geo_spatial = Geospatial_data(data_points,
-                                           points_are_lats_longs=True)
-        points_UTM = data_geo_spatial.get_data_points(absolute=True)
-        attributes = linear_function(points_UTM)
-        att = 'elevation'
-        
-        #Create .txt file
-        txt_file = tempfile.mktemp(".txt")
-        file = open(txt_file,"w")
-        file.write(" x,y," + att + " \n")
-        for data_point, attribute in map(None, points_UTM, attributes):
-            row = str(data_point[0]) + ',' + str(data_point[1]) \
-                  + ',' + str(attribute)
-            #print "row", row 
-            file.write(row + "\n")
-        file.close()
-
-
-        pts_file = tempfile.mktemp(".pts")        
-        convert = Geospatial_data(txt_file)
-        convert.export_points_file(pts_file)
-
-        #Check that values can be set from file
-        quantity.set_values_from_file(pts_file, att, 0,
-                                      'vertices', None)
-        answer = linear_function(quantity.domain.get_vertex_coordinates())
-        #print "quantity.vertex_values.flat", quantity.vertex_values.flat
-        #print "answer",answer
-        assert allclose(quantity.vertex_values.flat, answer)
-
-        #Check that values can be set from file
-        quantity.set_values(filename=pts_file,
-                            attribute_name=att, alpha=0)
-        answer = linear_function(quantity.domain.get_vertex_coordinates())
-        #print "quantity.vertex_values.flat", quantity.vertex_values.flat
-        #print "answer",answer
-        assert allclose(quantity.vertex_values.flat, answer)
-
-
-        #Check that values can be set from file using default attribute
-        quantity.set_values(filename=txt_file, alpha=0)
-        assert allclose(quantity.vertex_values.flat, answer)
-
-        #Cleanup
-        import os
-        os.remove(txt_file)
-        os.remove(pts_file)
-        
-    def verbose_test_set_values_from_UTM_pts(self):
-        quantity = Quantity(self.mesh_onslow)
-
-        #Get (enough) datapoints
-        data_points = [[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
-                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
-                       [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
-                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
-                       [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
-                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
-                       [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
-                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
-                       [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
-                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
-                       [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
-                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
-                       [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
-                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
-                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
-                       [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
-                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
-                       [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
-                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
-                       [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
-                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
-                       [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
-                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
-                       [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
-                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
-                       [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
-                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
-                       ]
-
-        data_geo_spatial = Geospatial_data(data_points,
-                                           points_are_lats_longs=True)
-        points_UTM = data_geo_spatial.get_data_points(absolute=True)
-        attributes = linear_function(points_UTM)
-        att = 'elevation'
-        
-        #Create .txt file
-        txt_file = tempfile.mktemp(".txt")
-        file = open(txt_file,"w")
-        file.write(" x,y," + att + " \n")
-        for data_point, attribute in map(None, points_UTM, attributes):
-            row = str(data_point[0]) + ',' + str(data_point[1]) \
-                  + ',' + str(attribute)
-            #print "row", row 
-            file.write(row + "\n")
-        file.close()
-
-
-        pts_file = tempfile.mktemp(".pts")        
-        convert = Geospatial_data(txt_file)
-        convert.export_points_file(pts_file)
-
-        #Check that values can be set from file
-        quantity.set_values_from_file(pts_file, att, 0,
-                                      'vertices', None, verbose = True,
-                                      max_read_lines=2)
-        answer = linear_function(quantity.domain.get_vertex_coordinates())
-        #print "quantity.vertex_values.flat", quantity.vertex_values.flat
-        #print "answer",answer
-        assert allclose(quantity.vertex_values.flat, answer)
-
-        #Check that values can be set from file
-        quantity.set_values(filename=pts_file,
-                            attribute_name=att, alpha=0)
-        answer = linear_function(quantity.domain.get_vertex_coordinates())
-        #print "quantity.vertex_values.flat", quantity.vertex_values.flat
-        #print "answer",answer
-        assert allclose(quantity.vertex_values.flat, answer)
-
-
-        #Check that values can be set from file using default attribute
-        quantity.set_values(filename=txt_file, alpha=0)
-        assert allclose(quantity.vertex_values.flat, answer)
-
-        #Cleanup
-        import os
-        os.remove(txt_file)
-        os.remove(pts_file)
-        
-    def test_set_values_from_file_with_georef1(self):
-
-        #Mesh in zone 56 (absolute coords)
-
-        x0 = 314036.58727982
-        y0 = 6224951.2960092
-
-        a = [x0+0.0, y0+0.0]
-        b = [x0+0.0, y0+2.0]
-        c = [x0+2.0, y0+0.0]
-        d = [x0+0.0, y0+4.0]
-        e = [x0+2.0, y0+2.0]
-        f = [x0+4.0, y0+0.0]
-
-        points = [a, b, c, d, e, f]
-
-        #bac, bce, ecf, dbe
-        elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
-
-        #absolute going in ..
-        mesh4 = Domain(points, elements,
-                       geo_reference = Geo_reference(56, 0, 0))
-        mesh4.check_integrity()
-        quantity = Quantity(mesh4)
-
-        #Get (enough) datapoints (relative to georef)
-        data_points_rel = [[ 0.66666667, 0.66666667],
-                       [ 1.33333333, 1.33333333],
-                       [ 2.66666667, 0.66666667],
-                       [ 0.66666667, 2.66666667],
-                       [ 0.0, 1.0],
-                       [ 0.0, 3.0],
-                       [ 1.0, 0.0],
-                       [ 1.0, 1.0],
-                       [ 1.0, 2.0],
-                       [ 1.0, 3.0],
-                       [ 2.0, 1.0],
-                       [ 3.0, 0.0],
-                       [ 3.0, 1.0]]
-
-        data_geo_spatial = Geospatial_data(data_points_rel,
-                         geo_reference = Geo_reference(56, x0, y0))
-        data_points_absolute = data_geo_spatial.get_data_points(absolute=True)
-        attributes = linear_function(data_points_absolute)
-        att = 'spam_and_eggs'
-        
-        #Create .txt file
-        ptsfile = tempfile.mktemp(".txt")
-        file = open(ptsfile,"w")
-        file.write(" x,y," + att + " \n")
-        for data_point, attribute in map(None, data_points_absolute
-                                         ,attributes):
-            row = str(data_point[0]) + ',' + str(data_point[1]) \
-                  + ',' + str(attribute)
-            file.write(row + "\n")
-        file.close()
-
-        #file = open(ptsfile, 'r')
-        #lines = file.readlines()
-        #file.close()
-     
-
-        #Check that values can be set from file
-        quantity.set_values(filename=ptsfile,
-                            attribute_name=att, alpha=0)
-        answer = linear_function(quantity.domain.get_vertex_coordinates())
-
-        assert allclose(quantity.vertex_values.flat, answer)
-
-
-        #Check that values can be set from file using default attribute
-        quantity.set_values(filename=ptsfile, alpha=0)
-        assert allclose(quantity.vertex_values.flat, answer)
-
-        #Cleanup
-        import os
-        os.remove(ptsfile)
-
-
-    def test_set_values_from_file_with_georef2(self):
-
-        #Mesh in zone 56 (relative coords)
-
-        x0 = 314036.58727982
-        y0 = 6224951.2960092
-        #x0 = 0.0
-        #y0 = 0.0
-
-        a = [0.0, 0.0]
-        b = [0.0, 2.0]
-        c = [2.0, 0.0]
-        d = [0.0, 4.0]
-        e = [2.0, 2.0]
-        f = [4.0, 0.0]
-
-        points = [a, b, c, d, e, f]
-
-        #bac, bce, ecf, dbe
-        elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
-
-        mesh4 = Domain(points, elements,
-                       geo_reference = Geo_reference(56, x0, y0))
-        mesh4.check_integrity()
-        quantity = Quantity(mesh4)
-
-        #Get (enough) datapoints
-        data_points = [[ x0+0.66666667, y0+0.66666667],
-                       [ x0+1.33333333, y0+1.33333333],
-                       [ x0+2.66666667, y0+0.66666667],
-                       [ x0+0.66666667, y0+2.66666667],
-                       [ x0+0.0, y0+1.0],
-                       [ x0+0.0, y0+3.0],
-                       [ x0+1.0, y0+0.0],
-                       [ x0+1.0, y0+1.0],
-                       [ x0+1.0, y0+2.0],
-                       [ x0+1.0, y0+3.0],
-                       [ x0+2.0, y0+1.0],
-                       [ x0+3.0, y0+0.0],
-                       [ x0+3.0, y0+1.0]]
-
-
-        data_geo_spatial = Geospatial_data(data_points,
-                         geo_reference = Geo_reference(56, 0, 0))
-        data_points_absolute = data_geo_spatial.get_data_points(absolute=True)
-        attributes = linear_function(data_points_absolute)
-        att = 'spam_and_eggs'
-        
-        #Create .txt file
-        ptsfile = tempfile.mktemp(".txt")
-        file = open(ptsfile,"w")
-        file.write(" x,y," + att + " \n")
-        for data_point, attribute in map(None, data_points_absolute
-                                         ,attributes):
-            row = str(data_point[0]) + ',' + str(data_point[1]) \
-                  + ',' + str(attribute)
-            file.write(row + "\n")
-        file.close()
-
-
-        #Check that values can be set from file
-        quantity.set_values(filename=ptsfile,
-                            attribute_name=att, alpha=0)
-        answer = linear_function(quantity.domain. \
-                                 get_vertex_coordinates(absolute=True))
-
-
-        assert allclose(quantity.vertex_values.flat, answer)
-
-
-        #Check that values can be set from file using default attribute
-        quantity.set_values(filename=ptsfile, alpha=0)
-        assert allclose(quantity.vertex_values.flat, answer)
-
-        #Cleanup
-        import os
-        os.remove(ptsfile)
-
-
-
 
     def test_set_values_from_quantity(self):
 
-        quantity1 = Quantity(self.mesh4)
+        quantity1 = Quantity(self.domain4)
         quantity1.set_vertex_values([0,1,2,3,4,5])
 
@@ -1226 +397 @@
 
 
-        quantity2 = Quantity(self.mesh4)
+        quantity2 = Quantity(self.domain4)
         quantity2.set_values(quantity=quantity1)
         assert allclose(quantity2.vertex_values,
@@ -1247 +418 @@
 
 
-    def Xtest_set_values_from_quantity_using_polygon(self):
-        """test_set_values_from_quantity_using_polygon(self):
-        
-        Check that polygon can be used to restrict set_values when
-        using another quantity as argument.
-        """
-        
-        # Create restricting polygon (containing node #4 (2,2) and 
-        # centroid of triangle #1 (bce)
-        polygon = [[1.0, 1.0], [4.0, 1.0],
-                   [4.0, 4.0], [1.0, 4.0]]
-        assert allclose(inside_polygon(self.mesh4.nodes, polygon), 4)                   
-        
-        quantity1 = Quantity(self.mesh4)
-        quantity1.set_vertex_values([0,1,2,3,4,5])
+
+
+    def test_overloading(self):
+
+        quantity1 = Quantity(self.domain4)
+        quantity1.set_values( [[0,1],[1,2],[2,3],[3,4]],
+                            location = 'vertices')
 
         assert allclose(quantity1.vertex_values,
-                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
-
-
-        quantity2 = Quantity(self.mesh4)
-        quantity2.set_values(quantity=quantity1,
-                             polygon=polygon)
-                             
-        msg = 'Only node #4(e) at (2,2) should have values applied '
-        assert allclose(quantity2.vertex_values,
-                        [[0,0,0], [0,0,4], [4,0,0], [0,0,4]]), msg        
-                        #bac,     bce,     ecf,     dbe
-                        
-
-
-    def test_overloading(self):
-
-        quantity1 = Quantity(self.mesh4)
-        quantity1.set_vertex_values([0,1,2,3,4,5])
-
-        assert allclose(quantity1.vertex_values,
-                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
-
-
-        quantity2 = Quantity(self.mesh4)
-        quantity2.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
+                        [[0,1], [1,2], [2,3], [3,4]])
+
+
+        quantity2 = Quantity(self.domain4)
+        quantity2.set_values([[1,2], [5,5], [0,9], [-6, 3]],
                              location = 'vertices')
 
 
 
-        quantity3 = Quantity(self.mesh4)
-        quantity3.set_values([[2,2,2], [7,8,9], [7,6,3], [3, 8, -8]],
+        quantity3 = Quantity(self.domain4)
+        quantity3.set_values([[2,2], [7,8], [7,6], [3, -8]],
                              location = 'vertices')
 
@@ -1302 +445 @@
         assert allclose(Q.vertex_values, -quantity1.vertex_values)
         assert allclose(Q.centroid_values, -quantity1.centroid_values)
-        assert allclose(Q.edge_values, -quantity1.edge_values)
+
 
         # Addition
@@ -1308 +451 @@
         assert allclose(Q.vertex_values, quantity1.vertex_values + 7)
         assert allclose(Q.centroid_values, quantity1.centroid_values + 7)
-        assert allclose(Q.edge_values, quantity1.edge_values + 7)
 
         Q = 7 + quantity1
         assert allclose(Q.vertex_values, quantity1.vertex_values + 7)
         assert allclose(Q.centroid_values, quantity1.centroid_values + 7)
-        assert allclose(Q.edge_values, quantity1.edge_values + 7)
 
         Q = quantity1 + quantity2
@@ -1320 +461 @@
         assert allclose(Q.centroid_values,
                         quantity1.centroid_values + quantity2.centroid_values)
-        assert allclose(Q.edge_values,
-                        quantity1.edge_values + quantity2.edge_values)
-
 
         Q = quantity1 + quantity2 - 3
@@ -1336 +474 @@
         assert allclose(Q.vertex_values, quantity1.vertex_values*3)
         assert allclose(Q.centroid_values, quantity1.centroid_values*3)
-        assert allclose(Q.edge_values, quantity1.edge_values*3)
+
         Q = 3*quantity1
         assert allclose(Q.vertex_values, quantity1.vertex_values*3)
@@ -1342 +480 @@
         #Multiplication
         Q = quantity1 * quantity2
-        #print Q.vertex_values
-        #print Q.centroid_values
-        #print quantity1.centroid_values
-        #print quantity2.centroid_values
-
         assert allclose(Q.vertex_values,
                         quantity1.vertex_values * quantity2.vertex_values)
@@ -1403 +536 @@
                         quantity2.vertex_values**2)**0.5)
 
-
-
-
-
-
-
     def test_compute_gradient(self):
-        quantity = Quantity(self.mesh4)
+        quantity = Quantity(self.domain6)
 
         #Set up for a gradient of (2,0) at mid triangle
-        quantity.set_values([2.0, 4.0, 6.0, 2.0],
+        quantity.set_values([2.0, 4.0, 4.0, 5.0, 10.0, 12.0],
                             location = 'centroids')
 
@@ -1420 +547 @@
         quantity.compute_gradients()
 
-        a = quantity.x_gradient
-        b = quantity.y_gradient
-        #print self.mesh4.centroid_coordinates
-        #print a, b
-
-        #The central triangle (1)
-        #(using standard gradient based on neigbours controid values)
-        assert allclose(a[1], 2.0)
-        assert allclose(b[1], 0.0)
-
-
-        #Left triangle (0) using two point gradient
-        #q0 = q1 + a*(x0-x1) + b*(y0-y1)  <=>
-        #2  = 4  + a*(-2/3)  + b*(-2/3)
-        assert allclose(a[0] + b[0], 3)
-        #From orthogonality (a*(y0-y1) + b*(x0-x1) == 0)
-        assert allclose(a[0] - b[0], 0)
-
-
-        #Right triangle (2) using two point gradient
-        #q2 = q1 + a*(x2-x1) + b*(y2-y1)  <=>
-        #6  = 4  + a*(4/3)  + b*(-2/3)
-        assert allclose(2*a[2] - b[2], 3)
-        #From orthogonality (a*(y1-y2) + b*(x2-x1) == 0)
-        assert allclose(a[2] + 2*b[2], 0)
-
-
-        #Top triangle (3) using two point gradient
-        #q3 = q1 + a*(x3-x1) + b*(y3-y1)  <=>
-        #2  = 4  + a*(-2/3)  + b*(4/3)
-        assert allclose(a[3] - 2*b[3], 3)
-        #From orthogonality (a*(y1-y3) + b*(x3-x1) == 0)
-        assert allclose(2*a[3] + b[3], 0)
-
-
-
-        #print a, b
+        a = quantity.gradients
+
+        assert allclose(a, [ 3., 1., 0.5, 3., 3.5, 0.5])
+        
         quantity.extrapolate_second_order()
 
-        #Apply q(x,y) = qc + a*(x-xc) + b*(y-yc)
-        assert allclose(quantity.vertex_values[0,:], [3., 0.,  3.])
-        assert allclose(quantity.vertex_values[1,:], [4./3, 16./3,  16./3])
-
-
-        #a = 1.2, b=-0.6
-        #q(4,0) = 6 + a*(4 - 8/3) + b*(-2/3)
-        assert allclose(quantity.vertex_values[2,2], 8)
-
-    def test_get_gradients(self):
-        quantity = Quantity(self.mesh4)
-
-        #Set up for a gradient of (2,0) at mid triangle
-        quantity.set_values([2.0, 4.0, 6.0, 2.0],
-                            location = 'centroids')
-
-
-        #Gradients
-        quantity.compute_gradients()
-
-        a, b = quantity.get_gradients()
-        #print self.mesh4.centroid_coordinates
-        #print a, b
-
-        #The central triangle (1)
-        #(using standard gradient based on neigbours controid values)
-        assert allclose(a[1], 2.0)
-        assert allclose(b[1], 0.0)
-
-
-        #Left triangle (0) using two point gradient
-        #q0 = q1 + a*(x0-x1) + b*(y0-y1)  <=>
-        #2  = 4  + a*(-2/3)  + b*(-2/3)
-        assert allclose(a[0] + b[0], 3)
-        #From orthogonality (a*(y0-y1) + b*(x0-x1) == 0)
-        assert allclose(a[0] - b[0], 0)
-
-
-        #Right triangle (2) using two point gradient
-        #q2 = q1 + a*(x2-x1) + b*(y2-y1)  <=>
-        #6  = 4  + a*(4/3)  + b*(-2/3)
-        assert allclose(2*a[2] - b[2], 3)
-        #From orthogonality (a*(y1-y2) + b*(x2-x1) == 0)
-        assert allclose(a[2] + 2*b[2], 0)
-
-
-        #Top triangle (3) using two point gradient
-        #q3 = q1 + a*(x3-x1) + b*(y3-y1)  <=>
-        #2  = 4  + a*(-2/3)  + b*(4/3)
-        assert allclose(a[3] - 2*b[3], 3)
-        #From orthogonality (a*(y1-y3) + b*(x3-x1) == 0)
-        assert allclose(2*a[3] + b[3], 0)
-
+
+        assert allclose(quantity.vertex_values, [[ 1., 3. ],
+                                                 [ 4., 4. ],
+                                                 [ 4., 4. ],
+                                                 [ 4., 6.],
+                                                 [ 8.25, 11.75],
+                                                 [ 11.,  13. ]])
+
+                                                 
 
     def test_second_order_extrapolation2(self):
-        quantity = Quantity(self.mesh4)
+        quantity = Quantity(self.domain4)
 
         #Set up for a gradient of (3,1), f(x) = 3x+y
@@ -1543 +593 @@
 
     def test_backup_saxpy_centroid_values(self):
-        quantity = Quantity(self.mesh4)
+        quantity = Quantity(self.domain4)
 
         #Set up for a gradient of (3,1), f(x) = 3x+y
@@ -1566 +616 @@
 
     def test_first_order_extrapolator(self):
-        quantity = Quantity(self.mesh4)
-
-        #Test centroids
-        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
-        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
+        quantity = Quantity(self.domain4)
+
+        centroid_values = array([1.,2.,3.,4.])
+        quantity.set_values(centroid_values, location = 'centroids')
+        assert allclose(quantity.centroid_values, centroid_values) #Centroid
 
         #Extrapolate
@@ -1576 +626 @@
 
         #Check that gradient is zero
-        a,b = quantity.get_gradients()
+        a = quantity.gradients
         assert allclose(a, [0,0,0,0])
-        assert allclose(b, [0,0,0,0])
+       
 
         #Check vertices but not edge values
         assert allclose(quantity.vertex_values,
-                        [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
+                        [[1,1], [2,2], [3,3], [4,4]])
 
 
     def test_second_order_extrapolator(self):
-        quantity = Quantity(self.mesh4)
+        quantity = Quantity(self.domain4)
 
         #Set up for a gradient of (3,0) at mid triangle
@@ -1616 +666 @@
 
 
-
-
-
-    def test_limit_vertices_by_all_neighbours(self):
-        quantity = Quantity(self.mesh4)
+    def test_limit(self):
+        quantity = Quantity(self.domain4)
 
         #Create a deliberate overshoot (e.g. from gradient computation)
-        quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
-
+        quantity.set_values([[0,0], [2,20], [-20,3], [8,3]])
 
         #Limit
-        quantity.limit_vertices_by_all_neighbours()
-
-        #Assert that central triangle is limited by neighbours
-        assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0]
-        assert quantity.vertex_values[1,0] <= quantity.vertex_values[3,1]
-
-        assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1]
-        assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
-
-        assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
-        assert quantity.vertex_values[1,2] <= quantity.vertex_values[3,1]
-
-
-
-        #Assert that quantities are conserved
-        from Numeric import sum
-        for k in range(quantity.centroid_values.shape[0]):
-            assert allclose (quantity.centroid_values[k],
-                             sum(quantity.vertex_values[k,:])/3)
-
-
-
-    def test_limit_edges_by_all_neighbours(self):
-        quantity = Quantity(self.mesh4)
-
-        #Create a deliberate overshoot (e.g. from gradient computation)
-        quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
-
-
-        #Limit
-        quantity.limit_edges_by_all_neighbours()
-
-        #Assert that central triangle is limited by neighbours
-        assert quantity.edge_values[1,0] <= quantity.centroid_values[2]
-        assert quantity.edge_values[1,0] >= quantity.centroid_values[0]
-
-        assert quantity.edge_values[1,1] <= quantity.centroid_values[2]
-        assert quantity.edge_values[1,1] >= quantity.centroid_values[0]
-
-        assert quantity.edge_values[1,2] <= quantity.centroid_values[2]
-        assert quantity.edge_values[1,2] >= quantity.centroid_values[0]
-
-
-
-        #Assert that quantities are conserved
-        from Numeric import sum
-        for k in range(quantity.centroid_values.shape[0]):
-            assert allclose (quantity.centroid_values[k],
-                             sum(quantity.vertex_values[k,:])/3)
-
-
-    def test_limit_edges_by_neighbour(self):
-        quantity = Quantity(self.mesh4)
-
-        #Create a deliberate overshoot (e.g. from gradient computation)
-        quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
-
-
-        #Limit
-        quantity.limit_edges_by_neighbour()
-
-        #Assert that central triangle is limited by neighbours
-        assert quantity.edge_values[1,0] <= quantity.centroid_values[3]
-        assert quantity.edge_values[1,0] >= quantity.centroid_values[1]
-
-        assert quantity.edge_values[1,1] <= quantity.centroid_values[2]
-        assert quantity.edge_values[1,1] >= quantity.centroid_values[1]
-
-        assert quantity.edge_values[1,2] <= quantity.centroid_values[1]
-        assert quantity.edge_values[1,2] >= quantity.centroid_values[0]
-
-
-
-        #Assert that quantities are conserved
-        from Numeric import sum
-        for k in range(quantity.centroid_values.shape[0]):
-            assert allclose (quantity.centroid_values[k],
-                             sum(quantity.vertex_values[k,:])/3)
-
-    def test_limiter2(self):
-        """Taken from test_shallow_water
-        """
-        quantity = Quantity(self.mesh4)
-        quantity.domain.beta_w = 0.9
-        
+        quantity.limit_minmod()
+
+       
+        #cells 1 and 2 are local max and min
+        assert quantity.vertex_values[1][0] == quantity.centroid_values[1]
+        assert quantity.vertex_values[1][1] == quantity.centroid_values[1]
+
+        assert quantity.vertex_values[2][0] == quantity.centroid_values[2]
+        assert quantity.vertex_values[2][1] == quantity.centroid_values[2]
+
+
+
+    def test_distribute_first_order(self):
+        quantity = Quantity(self.domain4)
+
         #Test centroids
-        quantity.set_values([2.,4.,8.,2.], location = 'centroids')
-        assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid
+        centroid_values = array([1.,2.,3.,4.])
+        quantity.set_values(centroid_values, location = 'centroids')
+        assert allclose(quantity.centroid_values, centroid_values) #Centroid
+
+
+        #Extrapolate from centroid to vertices and edges
+        quantity.extrapolate_first_order()
+        
+        assert allclose(quantity.vertex_values,[[ 1.,  1.],
+                                                [ 2.,  2.],
+                                                [ 3.,  3.],
+                                                [ 4.,  4.]])
+ 
+
+
+    def test_distribute_second_order(self):
+        quantity = Quantity(self.domain4)
+
+        #Test centroids
+        centroid_values = array([2.,4.,8.,2.])
+        quantity.set_values(centroid_values, location = 'centroids')
+        assert allclose(quantity.centroid_values, centroid_values) #Centroid
 
 
@@ -1720 +716 @@
         quantity.extrapolate_second_order()
 
-        assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6])
-
-        #Limit
-        quantity.limit()
-
-        # limited value for beta_w = 0.9
-        
-        assert allclose(quantity.vertex_values[1,:], [2.2, 4.9, 4.9])
-        # limited values for beta_w = 0.5
-        #assert allclose(quantity.vertex_values[1,:], [3.0, 4.5, 4.5])
-
-
-        #Assert that quantities are conserved
-        from Numeric import sum
-        for k in range(quantity.centroid_values.shape[0]):
-            assert allclose (quantity.centroid_values[k],
-                             sum(quantity.vertex_values[k,:])/3)
-
-
-
-
-
-    def test_distribute_first_order(self):
-        quantity = Quantity(self.mesh4)
+        assert allclose(quantity.vertex_values, [[ 1.2,  2.8],
+                                                 [ 2.,  6. ],
+                                                 [ 8.,   8. ],
+                                                 [ 6.8, -2.8]])
+
+
+    def test_update_explicit(self):
+        quantity = Quantity(self.domain4)
 
         #Test centroids
@@ -1749 +729 @@
         assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
 
-
-        #Extrapolate from centroid to vertices and edges
-        quantity.extrapolate_first_order()
-
-        #Interpolate
-        #quantity.interpolate_from_vertices_to_edges()
-
-        assert allclose(quantity.vertex_values,
-                        [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
-        assert allclose(quantity.edge_values, [[1,1,1], [2,2,2],
-                                               [3,3,3], [4, 4, 4]])
-
-
-    def test_interpolate_from_vertices_to_edges(self):
-        quantity = Quantity(self.mesh4)
-
-        quantity.vertex_values = array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]],Float)
-
-        quantity.interpolate_from_vertices_to_edges()
-
-        assert allclose(quantity.edge_values, [[1., 1.5, 0.5],
-                                               [3., 2.5, 1.5],
-                                               [3.5, 4.5, 3.],
-                                               [2.5, 3.5, 2]])
-
-
-    def test_interpolate_from_edges_to_vertices(self):
-        quantity = Quantity(self.mesh4)
-
-        quantity.edge_values = array([[1., 1.5, 0.5],
-                                [3., 2.5, 1.5],
-                                [3.5, 4.5, 3.],
-                                [2.5, 3.5, 2]],Float)
-
-        quantity.interpolate_from_edges_to_vertices()
-
-        assert allclose(quantity.vertex_values,
-                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
-
-
-
-    def test_distribute_second_order(self):
-        quantity = Quantity(self.mesh4)
-
-        #Test centroids
-        quantity.set_values([2.,4.,8.,2.], location = 'centroids')
-        assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid
-
-
-        #Extrapolate
-        quantity.extrapolate_second_order()
-
-        assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6])
-
-
-    def test_update_explicit(self):
-        quantity = Quantity(self.mesh4)
+        #Set explicit_update
+        quantity.explicit_update = array( [1.,1.,1.,1.] )
+
+        #Update with given timestep
+        quantity.update(0.1)
+
+        x = array([1, 2, 3, 4]) + array( [.1,.1,.1,.1] )
+        assert allclose( quantity.centroid_values, x)
+
+    def test_update_semi_implicit(self):
+        quantity = Quantity(self.domain4)
 
         #Test centroids
@@ -1811 +745 @@
         assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
 
-        #Set explicit_update
-        quantity.explicit_update = array( [1.,1.,1.,1.] )
-
-        #Update with given timestep
-        quantity.update(0.1)
-
-        x = array([1, 2, 3, 4]) + array( [.1,.1,.1,.1] )
-        assert allclose( quantity.centroid_values, x)
-
-    def test_update_semi_implicit(self):
-        quantity = Quantity(self.mesh4)
-
-        #Test centroids
-        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
-        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
-
         #Set semi implicit update
         quantity.semi_implicit_update = array([1.,1.,1.,1.])
@@ -1842 +760 @@
 
     def test_both_updates(self):
-        quantity = Quantity(self.mesh4)
+        quantity = Quantity(self.domain4)
 
         #Test centroids
-        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
-        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
+        centroid_values = array( [1, 2, 3, 4] )
+        quantity.set_values(centroid_values, location = 'centroids')
+        assert allclose(quantity.centroid_values, centroid_values) #Centroid
 
         #Set explicit_update
-        quantity.explicit_update = array( [4.,3.,2.,1.] )
+        explicit_update = array( [4.,3.,2.,1.] )
+        quantity.explicit_update[:] = explicit_update
 
         #Set semi implicit update
-        quantity.semi_implicit_update = array( [1.,1.,1.,1.] )
+        semi_implicit_update = array( [1.,1.,1.,1.] )
+        quantity.semi_implicit_update[:] = semi_implicit_update
 
         #Update with given timestep
@@ -1858 +779 @@
         quantity.update(0.1)
 
-        sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])
-        denom = ones(4, Float)-timestep*sem
-
-        x = array([1., 2., 3., 4.])
-        x /= denom
-        x += timestep*array( [4.0, 3.0, 2.0, 1.0] )
-
+        denom = 1.0-timestep*semi_implicit_update
+        x = (centroid_values + timestep*explicit_update)/denom
+ 
         assert allclose( quantity.centroid_values, x)
-
-
-
 
     #Test smoothing
     def test_smoothing(self):
 
-        from mesh_factory import rectangular
+        
         from shallow_water import Domain, Transmissive_boundary
         from Numeric import zeros, Float
         from anuga.utilities.numerical_tools import mean
 
-        #Create basic mesh
-        points, vertices, boundary = rectangular(2, 2)
 
         #Create shallow water domain
-        domain = Domain(points, vertices, boundary)
+        domain = Domain(points10)
         domain.default_order=2
         domain.reduction = mean
@@ -1888 +800 @@
 
         #Set some field values
-        domain.set_quantity('elevation', lambda x,y: x)
+        domain.set_quantity('elevation', lambda x: x)
         domain.set_quantity('friction', 0.03)
 
@@ -2030 +942 @@
 
 
-    def set_array_values_by_index(self):
-
-        from mesh_factory import rectangular
-        from shallow_water import Domain
-        from Numeric import zeros, Float
-
-        #Create basic mesh
-        points, vertices, boundary = rectangular(1, 1)
-
-        #Create shallow water domain
-        domain = Domain(points, vertices, boundary)
-        #print "domain.number_of_elements ",domain.number_of_elements
-        quantity = Quantity(domain,[[1,1,1],[2,2,2]])
-        value = [7]
-        indices = [1]
-        quantity.set_array_values_by_index(value,
-                                           location = 'centroids',
-                                           indices = indices)
-        #print "quantity.centroid_values",quantity.centroid_values
-
-        assert allclose(quantity.centroid_values, [1,7])
-
-        quantity.set_array_values([15,20,25], indices = indices)
-        assert allclose(quantity.centroid_values, [1,20])
-
-        quantity.set_array_values([15,20,25], indices = indices)
-        assert allclose(quantity.centroid_values, [1,20])
-
-    def test_setting_some_vertex_values(self):
-        """
-        set values based on triangle lists.
-        """
-        from mesh_factory import rectangular
-        from shallow_water import Domain
-        from Numeric import zeros, Float
-
-        #Create basic mesh
-        points, vertices, boundary = rectangular(1, 3)
-        #print "vertices",vertices
-        #Create shallow water domain
-        domain = Domain(points, vertices, boundary)
-        #print "domain.number_of_elements ",domain.number_of_elements
-        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
-                                    [4,4,4],[5,5,5],[6,6,6]])
-
-
-        # Check that constants work
-        value = 7
-        indices = [1]
-        quantity.set_values(value,
-                            location = 'centroids',
-                            indices = indices)
-        #print "quantity.centroid_values",quantity.centroid_values
-        assert allclose(quantity.centroid_values, [1,7,3,4,5,6])
-        
-        value = [7]
-        indices = [1]
-        quantity.set_values(value,
-                            location = 'centroids',
-                            indices = indices)
-        #print "quantity.centroid_values",quantity.centroid_values
-        assert allclose(quantity.centroid_values, [1,7,3,4,5,6])
-
-        value = [[15,20,25]]
-        quantity.set_values(value, indices = indices)
-        #print "1 quantity.vertex_values",quantity.vertex_values
-        assert allclose(quantity.vertex_values[1], value[0])
-
-
-        #print "quantity",quantity.vertex_values
-        values = [10,100,50]
-        quantity.set_values(values, indices = [0,1,5], location = 'centroids')
-        #print "2 quantity.vertex_values",quantity.vertex_values
-        assert allclose(quantity.vertex_values[0], [10,10,10])
-        assert allclose(quantity.vertex_values[5], [50,50,50])
-        #quantity.interpolate()
-        #print "quantity.centroid_values",quantity.centroid_values
-        assert allclose(quantity.centroid_values, [10,100,3,4,5,50])
-
-
-        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
-                                    [4,4,4],[5,5,5],[6,6,6]])
-        values = [10,100,50]
-        #this will be per unique vertex, indexing the vertices
-        #print "quantity.vertex_values",quantity.vertex_values
-        quantity.set_values(values, indices = [0,1,5])
-        #print "quantity.vertex_values",quantity.vertex_values
-        assert allclose(quantity.vertex_values[0], [1,50,10])
-        assert allclose(quantity.vertex_values[5], [6,6,6])
-        assert allclose(quantity.vertex_values[1], [100,10,50])
-
-        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
-                                    [4,4,4],[5,5,5],[6,6,6]])
-        values = [[31,30,29],[400,400,400],[1000,999,998]]
-        quantity.set_values(values, indices = [3,3,5])
-        quantity.interpolate()
-        assert allclose(quantity.centroid_values, [1,2,3,400,5,999])
-
-        values = [[1,1,1],[2,2,2],[3,3,3],
-                                    [4,4,4],[5,5,5],[6,6,6]]
-        quantity.set_values(values)
-
-        # testing the standard set values by vertex
-        # indexed by vertex_id in general_mesh.coordinates
-        values = [0,1,2,3,4,5,6,7]
-
-        quantity.set_values(values)
-        #print "1 quantity.vertex_values",quantity.vertex_values
-        assert allclose(quantity.vertex_values,[[ 4.,  5.,  0.],
-                                                [ 1.,  0.,  5.],
-                                                [ 5.,  6.,  1.],
-                                                [ 2.,  1.,  6.],
-                                                [ 6.,  7.,  2.],
-                                                [ 3.,  2.,  7.]])
-
-    def test_setting_unique_vertex_values(self):
-        """
-        set values based on unique_vertex lists.
-        """
-        from mesh_factory import rectangular
-        from shallow_water import Domain
-        from Numeric import zeros, Float
-
-        #Create basic mesh
-        points, vertices, boundary = rectangular(1, 3)
-        #print "vertices",vertices
-        #Create shallow water domain
-        domain = Domain(points, vertices, boundary)
-        #print "domain.number_of_elements ",domain.number_of_elements
-        quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
-                                    [4,4,4],[5,5,5]])
-        value = 7
-        indices = [1,5]
-        quantity.set_values(value,
-                            location = 'unique vertices',
-                            indices = indices)
-        #print "quantity.centroid_values",quantity.centroid_values
-        assert allclose(quantity.vertex_values[0], [0,7,0])
-        assert allclose(quantity.vertex_values[1], [7,1,7])
-        assert allclose(quantity.vertex_values[2], [7,2,7])
-
-
-    def test_get_values(self):
-        """
-        get values based on triangle lists.
-        """
-        from mesh_factory import rectangular
-        from shallow_water import Domain
-        from Numeric import zeros, Float
-
-        #Create basic mesh
-        points, vertices, boundary = rectangular(1, 3)
-
-        #print "points",points
-        #print "vertices",vertices
-        #print "boundary",boundary
-
-        #Create shallow water domain
-        domain = Domain(points, vertices, boundary)
-        #print "domain.number_of_elements ",domain.number_of_elements
-        quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
-                                    [4,4,4],[5,5,5]])
-
-        #print "quantity.get_values(location = 'unique vertices')", \
-        #      quantity.get_values(location = 'unique vertices')
-
-        #print "quantity.get_values(location = 'unique vertices')", \
-        #      quantity.get_values(indices=[0,1,2,3,4,5,6,7], \
-        #                          location = 'unique vertices')
-
-        answer = [0.5,2,4,5,0,1,3,4.5]
-        assert allclose(answer,
-                        quantity.get_values(location = 'unique vertices'))
-
-        indices = [0,5,3]
-        answer = [0.5,1,5]
-        assert allclose(answer,
-                        quantity.get_values(indices=indices, \
-                                            location = 'unique vertices'))
-        #print "quantity.centroid_values",quantity.centroid_values
-        #print "quantity.get_values(location = 'centroids') ",\
-        #      quantity.get_values(location = 'centroids')
-
-
-
-
-    def test_get_values_2(self):
-        """Different mesh (working with domain object) - also check centroids.
-        """
-
-        
-        a = [0.0, 0.0]
-        b = [0.0, 2.0]
-        c = [2.0,0.0]
-        d = [0.0, 4.0]
-        e = [2.0, 2.0]
-        f = [4.0,0.0]
-
-        points = [a, b, c, d, e, f]
-        #bac, bce, ecf, dbe
-        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
-
-        domain = Domain(points, vertices)
-
-        quantity = Quantity(domain)
-        quantity.set_values(lambda x, y: x+2*y) #2 4 4 6
-        
-        assert allclose(quantity.get_values(location='centroids'), [2,4,4,6])
-        assert allclose(quantity.get_values(location='centroids', indices=[1,3]), [4,6])
-
-
-        assert allclose(quantity.get_values(location='vertices'), [[4,0,2],
-                                                                   [4,2,6],
-                                                                   [6,2,4],
-                                                                   [8,4,6]])
-        
-        assert allclose(quantity.get_values(location='vertices', indices=[1,3]), [[4,2,6],
-                                                                                  [8,4,6]])
-
-
-        assert allclose(quantity.get_values(location='edges'), [[1,3,2],
-                                                                [4,5,3],
-                                                                [3,5,4],
-                                                                [5,7,6]])
-        assert allclose(quantity.get_values(location='edges', indices=[1,3]),
-                        [[4,5,3],
-                         [5,7,6]])        
-
-        # Check averaging over vertices
-        #a: 0
-        #b: (4+4+4)/3
-        #c: (2+2+2)/3
-        #d: 8
-        #e: (6+6+6)/3        
-        #f: 4
-        assert(quantity.get_values(location='unique vertices'), [0, 4, 2, 8, 6, 4])        
-                                                                                  
-        
-
-
-
-
-    def test_get_interpolated_values(self):
-
-        from mesh_factory import rectangular
-        from shallow_water import Domain
-        from Numeric import zeros, Float
-
-        #Create basic mesh
-        points, vertices, boundary = rectangular(1, 3)
-        domain = Domain(points, vertices, boundary)
-
-        #Constant values
-        quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
-                                    [4,4,4],[5,5,5]])
-
-        
-
-        # Get interpolated values at centroids
-        interpolation_points = domain.get_centroid_coordinates()
-        answer = quantity.get_values(location='centroids')
-
-        
-        #print quantity.get_values(points=interpolation_points)
-        assert allclose(answer, quantity.get_values(interpolation_points=interpolation_points))
-
-
-        #Arbitrary values
-        quantity = Quantity(domain,[[0,1,2],[3,1,7],[2,1,2],[3,3,7],
-                                    [1,4,-9],[2,5,0]])
-
-
-        # Get interpolated values at centroids
-        interpolation_points = domain.get_centroid_coordinates()
-        answer = quantity.get_values(location='centroids')
-        #print answer
-        #print quantity.get_values(points=interpolation_points)
-        assert allclose(answer, quantity.get_values(interpolation_points=interpolation_points))        
-                        
-
-        #FIXME TODO
-        #indices = [0,5,3]
-        #answer = [0.5,1,5]
-        #assert allclose(answer,
-        #                quantity.get_values(indices=indices, \
-        #                                    location = 'unique vertices'))
-
-
-
-
-    def test_get_interpolated_values_2(self):
-        a = [0.0, 0.0]
-        b = [0.0, 2.0]
-        c = [2.0,0.0]
-        d = [0.0, 4.0]
-        e = [2.0, 2.0]
-        f = [4.0,0.0]
-
-        points = [a, b, c, d, e, f]
-        #bac, bce, ecf, dbe
-        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
-
-        domain = Domain(points, vertices)
-
-        quantity = Quantity(domain)
-        quantity.set_values(lambda x, y: x+2*y) #2 4 4 6
-
-        #First pick one point
-        x, y = 2.0/3, 8.0/3
-        v = quantity.get_values(interpolation_points = [[x,y]])
-        assert allclose(v, 6)        
-
-        # Then another to test that algorithm won't blindly
-        # reuse interpolation matrix
-        x, y = 4.0/3, 4.0/3
-        v = quantity.get_values(interpolation_points = [[x,y]])
-        assert allclose(v, 4)        
-
-
-
-
-    def test_getting_some_vertex_values(self):
-        """
-        get values based on triangle lists.
-        """
-        from mesh_factory import rectangular
-        from shallow_water import Domain
-        from Numeric import zeros, Float
-
-        #Create basic mesh
-        points, vertices, boundary = rectangular(1, 3)
-
-        #print "points",points
-        #print "vertices",vertices
-        #print "boundary",boundary
-
-        #Create shallow water domain
-        domain = Domain(points, vertices, boundary)
-        #print "domain.number_of_elements ",domain.number_of_elements
-        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
-                                    [4,4,4],[5,5,5],[6,6,6]])
-        value = [7]
-        indices = [1]
-        quantity.set_values(value,
-                            location = 'centroids',
-                            indices = indices)
-        #print "quantity.centroid_values",quantity.centroid_values
-        #print "quantity.get_values(location = 'centroids') ",\
-        #      quantity.get_values(location = 'centroids')
-        assert allclose(quantity.centroid_values,
-                        quantity.get_values(location = 'centroids'))
-
-
-        value = [[15,20,25]]
-        quantity.set_values(value, indices = indices)
-        #print "1 quantity.vertex_values",quantity.vertex_values
-        assert allclose(quantity.vertex_values, quantity.get_values())
-
-        assert allclose(quantity.edge_values,
-                        quantity.get_values(location = 'edges'))
-
-        # get a subset of elements
-        subset = quantity.get_values(location='centroids', indices=[0,5])
-        answer = [quantity.centroid_values[0],quantity.centroid_values[5]]
-        assert allclose(subset, answer)
-
-
-        subset = quantity.get_values(location='edges', indices=[0,5])
-        answer = [quantity.edge_values[0],quantity.edge_values[5]]
-        #print "subset",subset
-        #print "answer",answer
-        assert allclose(subset, answer)
-
-        subset = quantity.get_values( indices=[1,5])
-        answer = [quantity.vertex_values[1],quantity.vertex_values[5]]
-        #print "subset",subset
-        #print "answer",answer
-        assert allclose(subset, answer)
-
-    def test_smooth_vertex_values(self):
-        """
-        get values based on triangle lists.
-        """
-        from mesh_factory import rectangular
-        from shallow_water import Domain
-        from Numeric import zeros, Float
-
-        #Create basic mesh
-        points, vertices, boundary = rectangular(2, 2)
-
-        #print "points",points
-        #print "vertices",vertices
-        #print "boundary",boundary
-
-        #Create shallow water domain
-        domain = Domain(points, vertices, boundary)
-        #print "domain.number_of_elements ",domain.number_of_elements
-        quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
-                                    [4,4,4],[5,5,5],[6,6,6],[7,7,7]])
-
-        #print "quantity.get_values(location = 'unique vertices')", \
-        #      quantity.get_values(location = 'unique vertices')
-
-        #print "quantity.get_values(location = 'unique vertices')", \
-        #      quantity.get_values(indices=[0,1,2,3,4,5,6,7], \
-        #                          location = 'unique vertices')
-
-        #print quantity.get_values(location = 'unique vertices')
-        #print quantity.domain.number_of_triangles_per_node
-        #print quantity.vertex_values
-        
-        #answer = [0.5, 2, 3, 3, 3.5, 4, 4, 5, 6.5]
-        #assert allclose(answer,
-        #                quantity.get_values(location = 'unique vertices'))
-
-        quantity.smooth_vertex_values()
-
-        #print quantity.vertex_values
-
-
-        answer_vertex_values = [[3,3.5,0.5],[2,0.5,3.5],[3.5,4,2],[3,2,4],
-                                [4,5,3],[3.5,3,5],[5,6.5,3.5],[4,3.5,6.5]]
-        
-        assert allclose(answer_vertex_values,
-                        quantity.vertex_values)
-        #print "quantity.centroid_values",quantity.centroid_values
-        #print "quantity.get_values(location = 'centroids') ",\
-        #      quantity.get_values(location = 'centroids')
-
-
-
+ 
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    suite = unittest.makeSuite(Test_Quantity, 'test')    
-    #suite = unittest.makeSuite(Test_Quantity, 'test_set_values_from_file_using_polygon')
-
-    #suite = unittest.makeSuite(Test_Quantity, 'test_set_vertex_values_using_general_interface_with_subset')
-    #print "restricted test"
-    #suite = unittest.makeSuite(Test_Quantity,'verbose_test_set_values_from_UTM_pts')
+    suite = unittest.makeSuite(Test_Quantity, 'test')
     runner = unittest.TextTestRunner()
     runner.run(suite)