Changeset 7815

Timestamp:

Jun 9, 2010, 5:34:19 PM (14 years ago)

Author:

steve

Message:

Almost got the unit test working for 1d quantity

Location:

anuga_work/development/anuga_1d

Files:

• quantity.py (modified) (6 diffs)
• test_quantity.py (modified) (10 diffs)

anuga_work/development/anuga_1d/quantity.py

@@ -271 +271 @@
             self.centroid_values[i] = (v0 + v1)/2.0
 
+
+
+
     def set_values(self, X, location='vertices'):
         """Set values for quantity
@@ -294 +297 @@
         """
 
-        if location not in ['vertices', 'centroids']:
-            msg = 'Invalid location: %s, (possible choices vertices, centroids)' %location
+        if location not in ['vertices', 'centroids', 'unique_vertices']:
+            msg = 'Invalid location: %s, (possible choices vertices, centroids, unique_vertices)' %location
             raise msg
 
@@ -309 +312 @@
           
         elif type(X) in [types.FloatType, types.IntType, types.LongType]:
-            if location == 'centroids':
-                self.centroid_values[:] = X
-            else:
-                self.vertex_values[:] = X
+            
+            self.centroid_values[:,] = float(X)
+            self.vertex_values[:,:] = float(X)
+
+        elif isinstance(X, Quantity):
+            self.set_array_values(X.vertex_values, location)
 
         else:
@@ -318 +323 @@
             self.set_array_values(X, location)
 
-        if location == 'vertices':
-            #Intialise centroid and edge values
+        if location == 'vertices' or location == 'unique_vertices':
+            #Intialise centroid
             self.interpolate()
+
+        if location == 'centroid':
+            self.extrapolate_first_order()
 
 
@@ -353 +361 @@
         values: Numeric array
         location: Where values are to be stored.
-                  Permissible options are: vertices, centroid, edges
+                  Permissible options are: vertices, centroid, unique_vertices
                   Default is "vertices"
 
         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 Nx2
+        of elements in the mesh. If location == 'unique_vertices' the
+        dimension values must be a list or a Numeric array of length N+1.
+        Otherwise it must be of dimension Nx2
 
         The values will be stored in elements following their
@@ -375 +385 @@
         N = self.centroid_values.shape[0]
 
-        msg = 'Number of values must match number of elements'
-        assert values.shape[0] == N, msg
 
         if location == 'centroids':
+            msg = 'Number of values must match number of elements'
+            assert values.shape[0] == N, msg
             assert len(values.shape) == 1, 'Values array must be 1d'
-            self.centroid_values = values
-        #elif location == 'edges':
-        #    assert len(values.shape) == 2, 'Values array must be 2d'
-        #    msg = 'Array must be N x 2'
-        #    self.edge_values = values
-        else:
+
+            for i in range(N):
+                self.centroid_values[i] = values[i]
+
+            self.vertex_values[:,0] = values
+            self.vertex_values[:,1] = values
+ 
+        if location == 'vertices':
+            msg = 'Number of values must match number of elements'
+            assert values.shape[0] == N, msg
             assert len(values.shape) == 2, 'Values array must be 2d'
             msg = 'Array must be N x 2'
             assert values.shape[1] == 2, msg
 
-            self.vertex_values[:] = values
+            self.vertex_values[:,:] = values[:,:]
+
+        if location == 'unique_vertices':
+            msg = 'Number of values must match number of elements +1'
+            assert values.shape[0] == N+1, msg
+            assert len(values.shape) == 1, 'Values array must be 1d'
+
+            self.vertex_values[:,0] = values[0:-1]
+            self.vertex_values[:,1] = values[1:N+1]  
+
+
+            
 
 

anuga_work/development/anuga_1d/test_quantity.py

@@ -108 +108 @@
         quantity = Quantity(self.domain4)
 
-
-        quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
-                            location = 'vertices')
+        # Test vertices
+        quantity.set_values([[1,2], [5,5], [0,9], [-6, 3]], location = 'vertices')
         assert allclose(quantity.vertex_values,
-                        [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
-        assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
-        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]])
-
-
-        # Test default
-        quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
+                        [[1,2], [5,5], [0,9], [-6, 3]])
+        assert allclose(quantity.centroid_values, [1.5, 5., 4.5, -1.5]) #Centroid
+
+
+
+        # Test unique_vertices
+        quantity.set_values([1,2,4,-5,6], location='unique_vertices')
         assert allclose(quantity.vertex_values,
-                        [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
-        assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
-        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]])
+                        [[1,2], [2,4], [4,-5], [-5,6]])
+        assert allclose(quantity.centroid_values, [1.5, 3., -.5, .5]) #Centroid
+
 
         # Test centroids
@@ -136 +129 @@
         # Test exceptions
         try:
-            quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
+            quantity.set_values([[1,3], [5,5], [0,9], [-6, 3]],
                                 location = 'bas kamel tuba')
         except:
@@ -143 +136 @@
 
         try:
-            quantity.set_values([[1,2,3], [0,0,9]])
+            quantity.set_values([[1,3], [0,9]])
         except AssertionError:
             pass
@@ -155 +148 @@
 
         quantity.set_values(1.0, location = 'vertices')
+
+        assert allclose(quantity.vertex_values, [[1,1], [1,1], [1,1], [1, 1]])
+        assert allclose(quantity.centroid_values, [1, 1, 1, 1]) #Centroid
+
+
+        quantity.set_values(2.0, location = 'centroids')
+
+        assert allclose(quantity.centroid_values, [2, 2, 2, 2])
+
+
+    def test_set_values_func(self):
+        quantity = Quantity(self.domain4)
+
+        def f(x):
+            return x*x
+
+
+
+        quantity.set_values(f, location = 'vertices')
+
+
         assert allclose(quantity.vertex_values,
-                        [[1,1,1], [1,1,1], [1,1,1], [1, 1, 1]])
-
-        assert allclose(quantity.centroid_values, [1, 1, 1, 1]) #Centroid
-        assert allclose(quantity.edge_values, [[1, 1, 1],
-                                               [1, 1, 1],
-                                               [1, 1, 1],
-                                               [1, 1, 1]])
-
-
-        quantity.set_values(2.0, location = 'centroids')
-        assert allclose(quantity.centroid_values, [2, 2, 2, 2])
-
-
-    def test_set_values_func(self):
-        quantity = Quantity(self.domain4)
-
-        def f(x, y):
-            return x+y
-
-        quantity.set_values(f, location = 'vertices')
-        #print "quantity.vertex_values",quantity.vertex_values
-        assert allclose(quantity.vertex_values,
-                        [[2,0,2], [2,2,4], [4,2,4], [4,2,4]])
+                        [[0,1], [1,6.25], [6.25,9], [9,25]])
+
         assert allclose(quantity.centroid_values,
-                        [4.0/3, 8.0/3, 10.0/3, 10.0/3])
-        assert allclose(quantity.edge_values,
-                        [[1,2,1], [3,3,2], [3,4,3], [3,4,3]])
+                        [0.5, 3.625, 7.625, 34*0.5])
 
 
         quantity.set_values(f, location = 'centroids')
+
+        
         assert allclose(quantity.centroid_values,
-                        [4.0/3, 8.0/3, 10.0/3, 10.0/3])
+                        [0.25, 3.0625, 7.5625, 16.0])
 
 
@@ -222 +215 @@
 
 
-    def test_set_vertex_values(self):
-        quantity = Quantity(self.domain4)
-        quantity.set_vertex_values([0,1,2,3,4,5])
-
-        assert allclose(quantity.vertex_values,
-                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
-        assert allclose(quantity.centroid_values,
-                        [1., 7./3, 11./3, 8./3]) #Centroid
-        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_set_vertex_values_subset(self):
-        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])
-
-        assert allclose(quantity.vertex_values,
-                        [[1,0,20], [1,20,4], [4,20,50], [30,1,4]])
-
-
-    def test_set_vertex_values_using_general_interface(self):
-        quantity = Quantity(self.domain4)
-
-
-        quantity.set_values([0,1,2,3,4,5])
-
-
-        assert allclose(quantity.vertex_values,
-                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
-
-        #Centroid
-        assert allclose(quantity.centroid_values, [1., 7./3, 11./3, 8./3])
-
-        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_set_vertex_values_using_general_interface_with_subset(self):
-        """test_set_vertex_values_using_general_interface_with_subset(self):
-        
-        Test that indices and polygon works (for constants values)
-        """
-        
-        quantity = Quantity(self.domain4)
-
-
-        quantity.set_values([0,2,3,5], indices=[0,2,3,5])
-        assert allclose(quantity.vertex_values,
-                        [[0,0,2], [0,2,0], [0,2,5], [3,0,0]])
-
-
-        # Constant
-        quantity.set_values(0.0)
-        quantity.set_values(3.14, indices=[0,2], location='vertices')
-
-        # Indices refer to triangle numbers
-        assert allclose(quantity.vertex_values,
-                        [[3.14,3.14,3.14], [0,0,0],
-                         [3.14,3.14,3.14], [0,0,0]])        
-        
-
-
-        # Now try with polygon (pick points where y>2)
-        polygon = [[0,2.1], [4,2.1], [4,7], [0,7]]
-        quantity.set_values(0.0)
-        quantity.set_values(3.14, polygon=polygon)
-        
-        assert allclose(quantity.vertex_values,
-                        [[0,0,0], [0,0,0], [0,0,0],
-                         [3.14,3.14,3.14]])                
-
-
-        # Another polygon (pick triangle 1 and 2 (rightmost triangles) 
-        # using centroids
-        polygon = [[2.1, 0.0], [3.5,0.1], [2,2.2], [0.2,2]]
-        quantity.set_values(0.0)
-        quantity.set_values(3.14, location='centroids', polygon=polygon)
-        assert allclose(quantity.vertex_values,
-                        [[0,0,0],
-                         [3.14,3.14,3.14],
-                         [3.14,3.14,3.14],                         
-                         [0,0,0]])                
-
-
-        # Same polygon now use vertices (default)
-        polygon = [[2.1, 0.0], [3.5,0.1], [2,2.2], [0.2,2]]
-        quantity.set_values(0.0)
-        #print 'Here 2'
-        quantity.set_values(3.14, polygon=polygon)
-        assert allclose(quantity.vertex_values,
-                        [[0,0,0],
-                         [3.14,3.14,3.14],
-                         [3.14,3.14,3.14],                         
-                         [0,0,0]])                
-        
-
-        # Test input checking
-        try:
-            quantity.set_values(3.14, polygon=polygon, indices = [0,2])
-        except:
-            pass
-        else:
-            msg = 'Should have caught this'
-            raise msg
-
-
-
-
-
-    def test_set_vertex_values_using_general_interface_subset_and_geo(self):
-        """test_set_vertex_values_using_general_interface_with_subset(self):
-        Test that indices and polygon works using georeferencing
-        """
-        
-        quantity = Quantity(self.domain4)
-        G = Geo_reference(56, 10, 100)
-        quantity.domain.geo_reference = G
-
-        #print quantity.domain.get_nodes(absolute=True)
-
-
-        # Constant
-        quantity.set_values(0.0)
-        quantity.set_values(3.14, indices=[0,2], location='vertices')
-
-        # Indices refer to triangle numbers here - not vertices (why?)
-        assert allclose(quantity.vertex_values,
-                        [[3.14,3.14,3.14], [0,0,0],
-                         [3.14,3.14,3.14], [0,0,0]])        
-        
-
-
-        # Now try with polygon (pick points where y>2)
-        polygon = array([[0,2.1], [4,2.1], [4,7], [0,7]])
-        polygon += [G.xllcorner, G.yllcorner]
-        
-        quantity.set_values(0.0)
-        quantity.set_values(3.14, polygon=polygon, location='centroids')
-        
-        assert allclose(quantity.vertex_values,
-                        [[0,0,0], [0,0,0], [0,0,0],
-                         [3.14,3.14,3.14]])                
-
-
-        # Another polygon (pick triangle 1 and 2 (rightmost triangles)
-        polygon = array([[2.1, 0.0], [3.5,0.1], [2,2.2], [0.2,2]])
-        polygon += [G.xllcorner, G.yllcorner]
-        
-        quantity.set_values(0.0)
-        quantity.set_values(3.14, polygon=polygon)
-
-        assert allclose(quantity.vertex_values,
-                        [[0,0,0],
-                         [3.14,3.14,3.14],
-                         [3.14,3.14,3.14],                         
-                         [0,0,0]])                
-
-
 
 
@@ -391 +220 @@
 
         quantity1 = Quantity(self.domain4)
-        quantity1.set_vertex_values([0,1,2,3,4,5])
+        quantity1.set_values([0,1,2,3,4], location='unique_vertices')
 
         assert allclose(quantity1.vertex_values,
-                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
+                        [[0,1], [1,2], [2,3], [3,4]])
 
 
         quantity2 = Quantity(self.domain4)
-        quantity2.set_values(quantity=quantity1)
+        quantity2.set_values(quantity1)
         assert allclose(quantity2.vertex_values,
-                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
-
-        quantity2.set_values(quantity = 2*quantity1)
+                        [[0,1], [1,2], [2,3], [3,4]])
+
+        quantity2.set_values(2*quantity1)
+
+        print quantity2.vertex_values
         assert allclose(quantity2.vertex_values,
-                        [[2,0,4], [2,4,8], [8,4,10], [6,2,8]])
-
-        quantity2.set_values(quantity = 2*quantity1 + 3)
-        assert allclose(quantity2.vertex_values,
-                        [[5,3,7], [5,7,11], [11,7,13], [9,5,11]])
-
-
-        #Check detection of quantity as first orgument
+                        [[0,2], [2,4], [4,6], [6,8]])
+
         quantity2.set_values(2*quantity1 + 3)
         assert allclose(quantity2.vertex_values,
-                        [[5,3,7], [5,7,11], [11,7,13], [9,5,11]])
-
-
+                        [[3,5], [5,7], [7,9], [9,11]])
 
 
@@ -573 +396 @@
         quantity.compute_gradients()
 
-        a = quantity.x_gradient
-        b = quantity.y_gradient
-        
-        #print a, b
-
-        assert allclose(a[1], 3.0)
-        assert allclose(b[1], 1.0)
+        a = quantity.gradients
+        
+
+        assert allclose(a[1], 2.8)
 
         #Work out the others
 
         quantity.extrapolate_second_order()
-
-        #print quantity.vertex_values
-        assert allclose(quantity.vertex_values[1,0], 2.0)
-        assert allclose(quantity.vertex_values[1,1], 6.0)
-        assert allclose(quantity.vertex_values[1,2], 8.0)
+        
+
+        assert allclose(quantity.vertex_values[1,0], 3.33333333)
+        assert allclose(quantity.vertex_values[1,1], 7.33333333)
+
+        assert allclose(quantity.centroid_values[1], 0.5*(7.33333333+3.33333333) )
+
 
 
@@ -611 +433 @@
         quantity.saxpy_centroid_values(2.0, 3.0)
 
-        assert(quantity.centroid_values, 2.0*d_values + 3.0*c_values)
+        assert allclose(quantity.centroid_values, 2.0*d_values + 3.0*c_values)
 
 
@@ -645 +467 @@
 
         quantity.extrapolate_second_order()
-        quantity.limit()
-
-
-        #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
+        from numpy import sum
         for k in range(quantity.centroid_values.shape[0]):
             assert allclose (quantity.centroid_values[k],
-                             sum(quantity.vertex_values[k,:])/3)
+                             sum(quantity.vertex_values[k,:])/2.0)
 
 
@@ -945 +755 @@
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    suite = unittest.makeSuite(Test_Quantity, 'test')
+    suite = unittest.makeSuite(Test_Quantity, 'test_set')
     runner = unittest.TextTestRunner()
     runner.run(suite)