Changeset 7815


Ignore:
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:
2 edited

Legend:

Unmodified
Added
Removed
  • anuga_work/development/anuga_1d/quantity.py

    r7793 r7815  
    271271            self.centroid_values[i] = (v0 + v1)/2.0
    272272
     273
     274
     275
    273276    def set_values(self, X, location='vertices'):
    274277        """Set values for quantity
     
    294297        """
    295298
    296         if location not in ['vertices', 'centroids']:
    297             msg = 'Invalid location: %s, (possible choices vertices, centroids)' %location
     299        if location not in ['vertices', 'centroids', 'unique_vertices']:
     300            msg = 'Invalid location: %s, (possible choices vertices, centroids, unique_vertices)' %location
    298301            raise msg
    299302
     
    309312         
    310313        elif type(X) in [types.FloatType, types.IntType, types.LongType]:
    311             if location == 'centroids':
    312                 self.centroid_values[:] = X
    313             else:
    314                 self.vertex_values[:] = X
     314           
     315            self.centroid_values[:,] = float(X)
     316            self.vertex_values[:,:] = float(X)
     317
     318        elif isinstance(X, Quantity):
     319            self.set_array_values(X.vertex_values, location)
    315320
    316321        else:
     
    318323            self.set_array_values(X, location)
    319324
    320         if location == 'vertices':
    321             #Intialise centroid and edge values
     325        if location == 'vertices' or location == 'unique_vertices':
     326            #Intialise centroid
    322327            self.interpolate()
     328
     329        if location == 'centroid':
     330            self.extrapolate_first_order()
    323331
    324332
     
    353361        values: Numeric array
    354362        location: Where values are to be stored.
    355                   Permissible options are: vertices, centroid, edges
     363                  Permissible options are: vertices, centroid, unique_vertices
    356364                  Default is "vertices"
    357365
    358366        In case of location == 'centroid' the dimension values must
    359367        be a list of a Numerical array of length N, N being the number
    360         of elements in the mesh. Otherwise it must be of dimension Nx2
     368        of elements in the mesh. If location == 'unique_vertices' the
     369        dimension values must be a list or a Numeric array of length N+1.
     370        Otherwise it must be of dimension Nx2
    361371
    362372        The values will be stored in elements following their
     
    375385        N = self.centroid_values.shape[0]
    376386
    377         msg = 'Number of values must match number of elements'
    378         assert values.shape[0] == N, msg
    379387
    380388        if location == 'centroids':
     389            msg = 'Number of values must match number of elements'
     390            assert values.shape[0] == N, msg
    381391            assert len(values.shape) == 1, 'Values array must be 1d'
    382             self.centroid_values = values
    383         #elif location == 'edges':
    384         #    assert len(values.shape) == 2, 'Values array must be 2d'
    385         #    msg = 'Array must be N x 2'
    386         #    self.edge_values = values
    387         else:
     392
     393            for i in range(N):
     394                self.centroid_values[i] = values[i]
     395
     396            self.vertex_values[:,0] = values
     397            self.vertex_values[:,1] = values
     398 
     399        if location == 'vertices':
     400            msg = 'Number of values must match number of elements'
     401            assert values.shape[0] == N, msg
    388402            assert len(values.shape) == 2, 'Values array must be 2d'
    389403            msg = 'Array must be N x 2'
    390404            assert values.shape[1] == 2, msg
    391405
    392             self.vertex_values[:] = values
     406            self.vertex_values[:,:] = values[:,:]
     407
     408        if location == 'unique_vertices':
     409            msg = 'Number of values must match number of elements +1'
     410            assert values.shape[0] == N+1, msg
     411            assert len(values.shape) == 1, 'Values array must be 1d'
     412
     413            self.vertex_values[:,0] = values[0:-1]
     414            self.vertex_values[:,1] = values[1:N+1] 
     415
     416
     417           
    393418
    394419
  • anuga_work/development/anuga_1d/test_quantity.py

    r7793 r7815  
    108108        quantity = Quantity(self.domain4)
    109109
    110 
    111         quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
    112                             location = 'vertices')
     110        # Test vertices
     111        quantity.set_values([[1,2], [5,5], [0,9], [-6, 3]], location = 'vertices')
    113112        assert allclose(quantity.vertex_values,
    114                         [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
    115         assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
    116         assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
    117                                                [5., 5., 5.],
    118                                                [4.5, 4.5, 0.],
    119                                                [3.0, -1.5, -1.5]])
    120 
    121 
    122         # Test default
    123         quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
     113                        [[1,2], [5,5], [0,9], [-6, 3]])
     114        assert allclose(quantity.centroid_values, [1.5, 5., 4.5, -1.5]) #Centroid
     115
     116
     117
     118        # Test unique_vertices
     119        quantity.set_values([1,2,4,-5,6], location='unique_vertices')
    124120        assert allclose(quantity.vertex_values,
    125                         [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
    126         assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
    127         assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
    128                                                [5., 5., 5.],
    129                                                [4.5, 4.5, 0.],
    130                                                [3.0, -1.5, -1.5]])
     121                        [[1,2], [2,4], [4,-5], [-5,6]])
     122        assert allclose(quantity.centroid_values, [1.5, 3., -.5, .5]) #Centroid
     123
    131124
    132125        # Test centroids
     
    136129        # Test exceptions
    137130        try:
    138             quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
     131            quantity.set_values([[1,3], [5,5], [0,9], [-6, 3]],
    139132                                location = 'bas kamel tuba')
    140133        except:
     
    143136
    144137        try:
    145             quantity.set_values([[1,2,3], [0,0,9]])
     138            quantity.set_values([[1,3], [0,9]])
    146139        except AssertionError:
    147140            pass
     
    155148
    156149        quantity.set_values(1.0, location = 'vertices')
     150
     151        assert allclose(quantity.vertex_values, [[1,1], [1,1], [1,1], [1, 1]])
     152        assert allclose(quantity.centroid_values, [1, 1, 1, 1]) #Centroid
     153
     154
     155        quantity.set_values(2.0, location = 'centroids')
     156
     157        assert allclose(quantity.centroid_values, [2, 2, 2, 2])
     158
     159
     160    def test_set_values_func(self):
     161        quantity = Quantity(self.domain4)
     162
     163        def f(x):
     164            return x*x
     165
     166
     167
     168        quantity.set_values(f, location = 'vertices')
     169
     170
    157171        assert allclose(quantity.vertex_values,
    158                         [[1,1,1], [1,1,1], [1,1,1], [1, 1, 1]])
    159 
    160         assert allclose(quantity.centroid_values, [1, 1, 1, 1]) #Centroid
    161         assert allclose(quantity.edge_values, [[1, 1, 1],
    162                                                [1, 1, 1],
    163                                                [1, 1, 1],
    164                                                [1, 1, 1]])
    165 
    166 
    167         quantity.set_values(2.0, location = 'centroids')
    168         assert allclose(quantity.centroid_values, [2, 2, 2, 2])
    169 
    170 
    171     def test_set_values_func(self):
    172         quantity = Quantity(self.domain4)
    173 
    174         def f(x, y):
    175             return x+y
    176 
    177         quantity.set_values(f, location = 'vertices')
    178         #print "quantity.vertex_values",quantity.vertex_values
    179         assert allclose(quantity.vertex_values,
    180                         [[2,0,2], [2,2,4], [4,2,4], [4,2,4]])
     172                        [[0,1], [1,6.25], [6.25,9], [9,25]])
     173
    181174        assert allclose(quantity.centroid_values,
    182                         [4.0/3, 8.0/3, 10.0/3, 10.0/3])
    183         assert allclose(quantity.edge_values,
    184                         [[1,2,1], [3,3,2], [3,4,3], [3,4,3]])
     175                        [0.5, 3.625, 7.625, 34*0.5])
    185176
    186177
    187178        quantity.set_values(f, location = 'centroids')
     179
     180       
    188181        assert allclose(quantity.centroid_values,
    189                         [4.0/3, 8.0/3, 10.0/3, 10.0/3])
     182                        [0.25, 3.0625, 7.5625, 16.0])
    190183
    191184
     
    222215
    223216
    224     def test_set_vertex_values(self):
    225         quantity = Quantity(self.domain4)
    226         quantity.set_vertex_values([0,1,2,3,4,5])
    227 
    228         assert allclose(quantity.vertex_values,
    229                         [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
    230         assert allclose(quantity.centroid_values,
    231                         [1., 7./3, 11./3, 8./3]) #Centroid
    232         assert allclose(quantity.edge_values, [[1., 1.5, 0.5],
    233                                                [3., 2.5, 1.5],
    234                                                [3.5, 4.5, 3.],
    235                                                [2.5, 3.5, 2]])
    236 
    237 
    238     def test_set_vertex_values_subset(self):
    239         quantity = Quantity(self.domain4)
    240         quantity.set_vertex_values([0,1,2,3,4,5])
    241         quantity.set_vertex_values([0,20,30,50], indices = [0,2,3,5])
    242 
    243         assert allclose(quantity.vertex_values,
    244                         [[1,0,20], [1,20,4], [4,20,50], [30,1,4]])
    245 
    246 
    247     def test_set_vertex_values_using_general_interface(self):
    248         quantity = Quantity(self.domain4)
    249 
    250 
    251         quantity.set_values([0,1,2,3,4,5])
    252 
    253 
    254         assert allclose(quantity.vertex_values,
    255                         [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
    256 
    257         #Centroid
    258         assert allclose(quantity.centroid_values, [1., 7./3, 11./3, 8./3])
    259 
    260         assert allclose(quantity.edge_values, [[1., 1.5, 0.5],
    261                                                [3., 2.5, 1.5],
    262                                                [3.5, 4.5, 3.],
    263                                                [2.5, 3.5, 2]])
    264 
    265 
    266 
    267     def test_set_vertex_values_using_general_interface_with_subset(self):
    268         """test_set_vertex_values_using_general_interface_with_subset(self):
    269        
    270         Test that indices and polygon works (for constants values)
    271         """
    272        
    273         quantity = Quantity(self.domain4)
    274 
    275 
    276         quantity.set_values([0,2,3,5], indices=[0,2,3,5])
    277         assert allclose(quantity.vertex_values,
    278                         [[0,0,2], [0,2,0], [0,2,5], [3,0,0]])
    279 
    280 
    281         # Constant
    282         quantity.set_values(0.0)
    283         quantity.set_values(3.14, indices=[0,2], location='vertices')
    284 
    285         # Indices refer to triangle numbers
    286         assert allclose(quantity.vertex_values,
    287                         [[3.14,3.14,3.14], [0,0,0],
    288                          [3.14,3.14,3.14], [0,0,0]])       
    289        
    290 
    291 
    292         # Now try with polygon (pick points where y>2)
    293         polygon = [[0,2.1], [4,2.1], [4,7], [0,7]]
    294         quantity.set_values(0.0)
    295         quantity.set_values(3.14, polygon=polygon)
    296        
    297         assert allclose(quantity.vertex_values,
    298                         [[0,0,0], [0,0,0], [0,0,0],
    299                          [3.14,3.14,3.14]])               
    300 
    301 
    302         # Another polygon (pick triangle 1 and 2 (rightmost triangles)
    303         # using centroids
    304         polygon = [[2.1, 0.0], [3.5,0.1], [2,2.2], [0.2,2]]
    305         quantity.set_values(0.0)
    306         quantity.set_values(3.14, location='centroids', polygon=polygon)
    307         assert allclose(quantity.vertex_values,
    308                         [[0,0,0],
    309                          [3.14,3.14,3.14],
    310                          [3.14,3.14,3.14],                         
    311                          [0,0,0]])               
    312 
    313 
    314         # Same polygon now use vertices (default)
    315         polygon = [[2.1, 0.0], [3.5,0.1], [2,2.2], [0.2,2]]
    316         quantity.set_values(0.0)
    317         #print 'Here 2'
    318         quantity.set_values(3.14, polygon=polygon)
    319         assert allclose(quantity.vertex_values,
    320                         [[0,0,0],
    321                          [3.14,3.14,3.14],
    322                          [3.14,3.14,3.14],                         
    323                          [0,0,0]])               
    324        
    325 
    326         # Test input checking
    327         try:
    328             quantity.set_values(3.14, polygon=polygon, indices = [0,2])
    329         except:
    330             pass
    331         else:
    332             msg = 'Should have caught this'
    333             raise msg
    334 
    335 
    336 
    337 
    338 
    339     def test_set_vertex_values_using_general_interface_subset_and_geo(self):
    340         """test_set_vertex_values_using_general_interface_with_subset(self):
    341         Test that indices and polygon works using georeferencing
    342         """
    343        
    344         quantity = Quantity(self.domain4)
    345         G = Geo_reference(56, 10, 100)
    346         quantity.domain.geo_reference = G
    347 
    348         #print quantity.domain.get_nodes(absolute=True)
    349 
    350 
    351         # Constant
    352         quantity.set_values(0.0)
    353         quantity.set_values(3.14, indices=[0,2], location='vertices')
    354 
    355         # Indices refer to triangle numbers here - not vertices (why?)
    356         assert allclose(quantity.vertex_values,
    357                         [[3.14,3.14,3.14], [0,0,0],
    358                          [3.14,3.14,3.14], [0,0,0]])       
    359        
    360 
    361 
    362         # Now try with polygon (pick points where y>2)
    363         polygon = array([[0,2.1], [4,2.1], [4,7], [0,7]])
    364         polygon += [G.xllcorner, G.yllcorner]
    365        
    366         quantity.set_values(0.0)
    367         quantity.set_values(3.14, polygon=polygon, location='centroids')
    368        
    369         assert allclose(quantity.vertex_values,
    370                         [[0,0,0], [0,0,0], [0,0,0],
    371                          [3.14,3.14,3.14]])               
    372 
    373 
    374         # Another polygon (pick triangle 1 and 2 (rightmost triangles)
    375         polygon = array([[2.1, 0.0], [3.5,0.1], [2,2.2], [0.2,2]])
    376         polygon += [G.xllcorner, G.yllcorner]
    377        
    378         quantity.set_values(0.0)
    379         quantity.set_values(3.14, polygon=polygon)
    380 
    381         assert allclose(quantity.vertex_values,
    382                         [[0,0,0],
    383                          [3.14,3.14,3.14],
    384                          [3.14,3.14,3.14],                         
    385                          [0,0,0]])               
    386 
    387 
    388217
    389218
     
    391220
    392221        quantity1 = Quantity(self.domain4)
    393         quantity1.set_vertex_values([0,1,2,3,4,5])
     222        quantity1.set_values([0,1,2,3,4], location='unique_vertices')
    394223
    395224        assert allclose(quantity1.vertex_values,
    396                         [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
     225                        [[0,1], [1,2], [2,3], [3,4]])
    397226
    398227
    399228        quantity2 = Quantity(self.domain4)
    400         quantity2.set_values(quantity=quantity1)
     229        quantity2.set_values(quantity1)
    401230        assert allclose(quantity2.vertex_values,
    402                         [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
    403 
    404         quantity2.set_values(quantity = 2*quantity1)
     231                        [[0,1], [1,2], [2,3], [3,4]])
     232
     233        quantity2.set_values(2*quantity1)
     234
     235        print quantity2.vertex_values
    405236        assert allclose(quantity2.vertex_values,
    406                         [[2,0,4], [2,4,8], [8,4,10], [6,2,8]])
    407 
    408         quantity2.set_values(quantity = 2*quantity1 + 3)
    409         assert allclose(quantity2.vertex_values,
    410                         [[5,3,7], [5,7,11], [11,7,13], [9,5,11]])
    411 
    412 
    413         #Check detection of quantity as first orgument
     237                        [[0,2], [2,4], [4,6], [6,8]])
     238
    414239        quantity2.set_values(2*quantity1 + 3)
    415240        assert allclose(quantity2.vertex_values,
    416                         [[5,3,7], [5,7,11], [11,7,13], [9,5,11]])
    417 
    418 
     241                        [[3,5], [5,7], [7,9], [9,11]])
    419242
    420243
     
    573396        quantity.compute_gradients()
    574397
    575         a = quantity.x_gradient
    576         b = quantity.y_gradient
    577        
    578         #print a, b
    579 
    580         assert allclose(a[1], 3.0)
    581         assert allclose(b[1], 1.0)
     398        a = quantity.gradients
     399       
     400
     401        assert allclose(a[1], 2.8)
    582402
    583403        #Work out the others
    584404
    585405        quantity.extrapolate_second_order()
    586 
    587         #print quantity.vertex_values
    588         assert allclose(quantity.vertex_values[1,0], 2.0)
    589         assert allclose(quantity.vertex_values[1,1], 6.0)
    590         assert allclose(quantity.vertex_values[1,2], 8.0)
     406       
     407
     408        assert allclose(quantity.vertex_values[1,0], 3.33333333)
     409        assert allclose(quantity.vertex_values[1,1], 7.33333333)
     410
     411        assert allclose(quantity.centroid_values[1], 0.5*(7.33333333+3.33333333) )
     412
    591413
    592414
     
    611433        quantity.saxpy_centroid_values(2.0, 3.0)
    612434
    613         assert(quantity.centroid_values, 2.0*d_values + 3.0*c_values)
     435        assert allclose(quantity.centroid_values, 2.0*d_values + 3.0*c_values)
    614436
    615437
     
    645467
    646468        quantity.extrapolate_second_order()
    647         quantity.limit()
    648 
    649 
    650         #Assert that central triangle is limited by neighbours
    651         assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0]
    652         assert quantity.vertex_values[1,0] >= quantity.vertex_values[3,1]
    653 
    654         assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1]
    655         assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
    656 
    657         assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
    658         assert quantity.vertex_values[1,2] >= quantity.vertex_values[3,1]
    659469
    660470
    661471        #Assert that quantities are conserved
    662         from Numeric import sum
     472        from numpy import sum
    663473        for k in range(quantity.centroid_values.shape[0]):
    664474            assert allclose (quantity.centroid_values[k],
    665                              sum(quantity.vertex_values[k,:])/3)
     475                             sum(quantity.vertex_values[k,:])/2.0)
    666476
    667477
     
    945755#-------------------------------------------------------------
    946756if __name__ == "__main__":
    947     suite = unittest.makeSuite(Test_Quantity, 'test')
     757    suite = unittest.makeSuite(Test_Quantity, 'test_set')
    948758    runner = unittest.TextTestRunner()
    949759    runner.run(suite)
Note: See TracChangeset for help on using the changeset viewer.