Changeset 7815 for anuga_work/development/anuga_1d
 Timestamp:
 Jun 9, 2010, 5:34:19 PM (14 years ago)
 Location:
 anuga_work/development/anuga_1d
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

anuga_work/development/anuga_1d/quantity.py
r7793 r7815 271 271 self.centroid_values[i] = (v0 + v1)/2.0 272 272 273 274 275 273 276 def set_values(self, X, location='vertices'): 274 277 """Set values for quantity … … 294 297 """ 295 298 296 if location not in ['vertices', 'centroids' ]:297 msg = 'Invalid location: %s, (possible choices vertices, centroids )' %location299 if location not in ['vertices', 'centroids', 'unique_vertices']: 300 msg = 'Invalid location: %s, (possible choices vertices, centroids, unique_vertices)' %location 298 301 raise msg 299 302 … … 309 312 310 313 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) 315 320 316 321 else: … … 318 323 self.set_array_values(X, location) 319 324 320 if location == 'vertices' :321 #Intialise centroid and edge values325 if location == 'vertices' or location == 'unique_vertices': 326 #Intialise centroid 322 327 self.interpolate() 328 329 if location == 'centroid': 330 self.extrapolate_first_order() 323 331 324 332 … … 353 361 values: Numeric array 354 362 location: Where values are to be stored. 355 Permissible options are: vertices, centroid, edges363 Permissible options are: vertices, centroid, unique_vertices 356 364 Default is "vertices" 357 365 358 366 In case of location == 'centroid' the dimension values must 359 367 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 361 371 362 372 The values will be stored in elements following their … … 375 385 N = self.centroid_values.shape[0] 376 386 377 msg = 'Number of values must match number of elements'378 assert values.shape[0] == N, msg379 387 380 388 if location == 'centroids': 389 msg = 'Number of values must match number of elements' 390 assert values.shape[0] == N, msg 381 391 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 388 402 assert len(values.shape) == 2, 'Values array must be 2d' 389 403 msg = 'Array must be N x 2' 390 404 assert values.shape[1] == 2, msg 391 405 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 393 418 394 419 
anuga_work/development/anuga_1d/test_quantity.py
r7793 r7815 108 108 quantity = Quantity(self.domain4) 109 109 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') 113 112 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') 124 120 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 131 124 132 125 # Test centroids … … 136 129 # Test exceptions 137 130 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]], 139 132 location = 'bas kamel tuba') 140 133 except: … … 143 136 144 137 try: 145 quantity.set_values([[1, 2,3], [0,0,9]])138 quantity.set_values([[1,3], [0,9]]) 146 139 except AssertionError: 147 140 pass … … 155 148 156 149 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 157 171 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 181 174 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]) 185 176 186 177 187 178 quantity.set_values(f, location = 'centroids') 179 180 188 181 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]) 190 183 191 184 … … 222 215 223 216 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]) #Centroid232 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 #Centroid258 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 # Constant282 quantity.set_values(0.0)283 quantity.set_values(3.14, indices=[0,2], location='vertices')284 285 # Indices refer to triangle numbers286 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 centroids304 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 checking327 try:328 quantity.set_values(3.14, polygon=polygon, indices = [0,2])329 except:330 pass331 else:332 msg = 'Should have caught this'333 raise msg334 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 georeferencing342 """343 344 quantity = Quantity(self.domain4)345 G = Geo_reference(56, 10, 100)346 quantity.domain.geo_reference = G347 348 #print quantity.domain.get_nodes(absolute=True)349 350 351 # Constant352 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 388 217 389 218 … … 391 220 392 221 quantity1 = Quantity(self.domain4) 393 quantity1.set_v ertex_values([0,1,2,3,4,5])222 quantity1.set_values([0,1,2,3,4], location='unique_vertices') 394 223 395 224 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]]) 397 226 398 227 399 228 quantity2 = Quantity(self.domain4) 400 quantity2.set_values(quantity =quantity1)229 quantity2.set_values(quantity1) 401 230 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 405 236 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 414 239 quantity2.set_values(2*quantity1 + 3) 415 240 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]]) 419 242 420 243 … … 573 396 quantity.compute_gradients() 574 397 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) 582 402 583 403 #Work out the others 584 404 585 405 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 591 413 592 414 … … 611 433 quantity.saxpy_centroid_values(2.0, 3.0) 612 434 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) 614 436 615 437 … … 645 467 646 468 quantity.extrapolate_second_order() 647 quantity.limit()648 649 650 #Assert that central triangle is limited by neighbours651 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]659 469 660 470 661 471 #Assert that quantities are conserved 662 from Numericimport sum472 from numpy import sum 663 473 for k in range(quantity.centroid_values.shape[0]): 664 474 assert allclose (quantity.centroid_values[k], 665 sum(quantity.vertex_values[k,:])/ 3)475 sum(quantity.vertex_values[k,:])/2.0) 666 476 667 477 … … 945 755 # 946 756 if __name__ == "__main__": 947 suite = unittest.makeSuite(Test_Quantity, 'test ')757 suite = unittest.makeSuite(Test_Quantity, 'test_set') 948 758 runner = unittest.TextTestRunner() 949 759 runner.run(suite)
Note: See TracChangeset
for help on using the changeset viewer.