Changeset 5899 for anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/test_quantity.py
 Timestamp:
 Nov 6, 2008, 12:28:22 PM (15 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/test_quantity.py
r5776 r5899 1 ## Automatically adapted for numpy.oldnumeric Oct 28, 2008 by alter_code1.py 2 1 3 #!/usr/bin/env python 2 4 … … 7 9 from quantity import * 8 10 from anuga.config import epsilon 9 from Numeric import allclose, array, ones, Float 11 #from numpy.oldnumeric import allclose, array, ones, Float 12 from numpy import allclose, array, ones, float 10 13 11 14 from anuga.fit_interpolate.fit import fit_to_mesh … … 60 63 61 64 62 def test_creation(self):63 64 quantity = Quantity(self.mesh1, [[1,2,3]])65 assert allclose(quantity.vertex_values, [[1.,2.,3.]])66 67 try:68 quantity = Quantity()69 except:70 pass71 else:72 raise 'Should have raised empty quantity exception'73 74 75 try:76 quantity = Quantity([1,2,3])77 except AssertionError:78 pass79 except:80 raise 'Should have raised "mising mesh object" error'81 82 83 def test_creation_zeros(self):84 85 quantity = Quantity(self.mesh1)86 assert allclose(quantity.vertex_values, [[0.,0.,0.]])87 88 89 quantity = Quantity(self.mesh4)90 assert allclose(quantity.vertex_values, [[0.,0.,0.], [0.,0.,0.],91 [0.,0.,0.], [0.,0.,0.]])92 93 94 def test_interpolation(self):95 quantity = Quantity(self.mesh1, [[1,2,3]])96 assert allclose(quantity.centroid_values, [2.0]) #Centroid97 98 assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5]])99 100 101 def test_interpolation2(self):102 quantity = Quantity(self.mesh4,103 [[1,2,3], [5,5,5], [0,0,9], [6, 3, 3]])104 assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid105 106 107 quantity.extrapolate_second_order()108 109 #print quantity.vertex_values110 assert allclose(quantity.vertex_values, [[3.5, 1.0, 3.5],111 [3.+2./3, 6.+2./3, 4.+2./3],112 [4.6, 3.4, 1.],113 [5.0, 1.0, 4.0]])114 115 #print quantity.edge_values116 assert allclose(quantity.edge_values, [[1.25, 3.5, 1.25],117 [5. + 2/3.0, 4.0 + 1.0/6, 5.0 + 1.0/6],118 [2.2, 2.8, 4.0],119 [2.5, 0.5, 2.0]])120 121 122 #assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],123 # [5., 5., 5.],124 # [4.5, 4.5, 0.],125 # [3.0, 1.5, 1.5]])126 65 ## def test_creation(self): 66 ## 67 ## quantity = Quantity(self.mesh1, [[1,2,3]]) 68 ## assert allclose(quantity.vertex_values, [[1.,2.,3.]]) 69 ## 70 ## try: 71 ## quantity = Quantity() 72 ## except: 73 ## pass 74 ## else: 75 ## raise 'Should have raised empty quantity exception' 76 ## 77 ## 78 ## try: 79 ## quantity = Quantity([1,2,3]) 80 ## except AssertionError: 81 ## pass 82 ## except: 83 ## raise 'Should have raised "mising mesh object" error' 84 ## 85 ## 86 ## def test_creation_zeros(self): 87 ## 88 ## quantity = Quantity(self.mesh1) 89 ## assert allclose(quantity.vertex_values, [[0.,0.,0.]]) 90 ## 91 ## 92 ## quantity = Quantity(self.mesh4) 93 ## assert allclose(quantity.vertex_values, [[0.,0.,0.], [0.,0.,0.], 94 ## [0.,0.,0.], [0.,0.,0.]]) 95 ## 96 ## 97 ## def test_interpolation(self): 98 ## quantity = Quantity(self.mesh1, [[1,2,3]]) 99 ## assert allclose(quantity.centroid_values, [2.0]) #Centroid 100 ## 101 ## assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5]]) 102 ## 103 ## 104 ## def test_interpolation2(self): 105 ## quantity = Quantity(self.mesh4, 106 ## [[1,2,3], [5,5,5], [0,0,9], [6, 3, 3]]) 107 ## assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid 108 ## 109 ## 110 ## quantity.extrapolate_second_order() 111 ## 112 ## #print quantity.vertex_values 113 ## assert allclose(quantity.vertex_values, [[3.5, 1.0, 3.5], 114 ## [3.+2./3, 6.+2./3, 4.+2./3], 115 ## [4.6, 3.4, 1.], 116 ## [5.0, 1.0, 4.0]]) 117 ## 118 ## #print quantity.edge_values 119 ## assert allclose(quantity.edge_values, [[1.25, 3.5, 1.25], 120 ## [5. + 2/3.0, 4.0 + 1.0/6, 5.0 + 1.0/6], 121 ## [2.2, 2.8, 4.0], 122 ## [2.5, 0.5, 2.0]]) 123 ## 124 ## 125 ## #assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5], 126 ## # [5., 5., 5.], 127 ## # [4.5, 4.5, 0.], 128 ## # [3.0, 1.5, 1.5]]) 129 ## 127 130 def test_get_extrema_1(self): 128 131 quantity = Quantity(self.mesh4, … … 156 159 157 160 158 def test_get_maximum_2(self): 159 160 a = [0.0, 0.0] 161 b = [0.0, 2.0] 162 c = [2.0,0.0] 163 d = [0.0, 4.0] 164 e = [2.0, 2.0] 165 f = [4.0,0.0] 166 167 points = [a, b, c, d, e, f] 168 #bac, bce, ecf, dbe 169 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 170 171 domain = Domain(points, vertices) 172 173 quantity = Quantity(domain) 174 quantity.set_values(lambda x, y: x+2*y) #2 4 4 6 175 176 v = quantity.get_maximum_value() 177 assert v == 6 178 179 v = quantity.get_minimum_value() 180 assert v == 2 181 182 i = quantity.get_maximum_index() 183 assert i == 3 184 185 i = quantity.get_minimum_index() 186 assert i == 0 187 188 x,y = quantity.get_maximum_location() 189 xref, yref = 2.0/3, 8.0/3 190 assert x == xref 191 assert y == yref 192 193 v = quantity.get_values(interpolation_points = [[x,y]]) 194 assert allclose(v, 6) 195 196 x,y = quantity.get_minimum_location() 197 v = quantity.get_values(interpolation_points = [[x,y]]) 198 assert allclose(v, 2) 199 200 #Multiple locations for maximum  201 #Test that the algorithm picks the first occurrence 202 v = quantity.get_maximum_value(indices=[0,1,2]) 203 assert allclose(v, 4) 204 205 i = quantity.get_maximum_index(indices=[0,1,2]) 206 assert i == 1 207 208 x,y = quantity.get_maximum_location(indices=[0,1,2]) 209 xref, yref = 4.0/3, 4.0/3 210 assert x == xref 211 assert y == yref 212 213 v = quantity.get_values(interpolation_points = [[x,y]]) 214 assert allclose(v, 4) 215 216 # More test of indices...... 217 v = quantity.get_maximum_value(indices=[2,3]) 218 assert allclose(v, 6) 219 220 i = quantity.get_maximum_index(indices=[2,3]) 221 assert i == 3 222 223 x,y = quantity.get_maximum_location(indices=[2,3]) 224 xref, yref = 2.0/3, 8.0/3 225 assert x == xref 226 assert y == yref 227 228 v = quantity.get_values(interpolation_points = [[x,y]]) 229 assert allclose(v, 6) 230 231 232 233 def test_boundary_allocation(self): 234 quantity = Quantity(self.mesh4, 235 [[1,2,3], [5,5,5], [0,0,9], [6, 3, 3]]) 236 237 assert quantity.boundary_values.shape[0] == len(self.mesh4.boundary) 238 239 240 def test_set_values(self): 241 quantity = Quantity(self.mesh4) 242 243 244 quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [6, 3, 3]], 245 location = 'vertices') 246 assert allclose(quantity.vertex_values, 247 [[1,2,3], [5,5,5], [0,0,9], [6, 3, 3]]) 248 assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid 249 assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5], 250 [5., 5., 5.], 251 [4.5, 4.5, 0.], 252 [3.0, 1.5, 1.5]]) 253 254 255 # Test default 256 quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [6, 3, 3]]) 257 assert allclose(quantity.vertex_values, 258 [[1,2,3], [5,5,5], [0,0,9], [6, 3, 3]]) 259 assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid 260 assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5], 261 [5., 5., 5.], 262 [4.5, 4.5, 0.], 263 [3.0, 1.5, 1.5]]) 264 265 # Test centroids 266 quantity.set_values([1,2,3,4], location = 'centroids') 267 assert allclose(quantity.centroid_values, [1., 2., 3., 4.]) #Centroid 268 269 # Test exceptions 270 try: 271 quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [6, 3, 3]], 272 location = 'bas kamel tuba') 273 except: 274 pass 275 276 277 try: 278 quantity.set_values([[1,2,3], [0,0,9]]) 279 except AssertionError: 280 pass 281 except: 282 raise 'should have raised Assertionerror' 283 284 285 286 def test_set_values_const(self): 287 quantity = Quantity(self.mesh4) 288 289 quantity.set_values(1.0, location = 'vertices') 290 assert allclose(quantity.vertex_values, 291 [[1,1,1], [1,1,1], [1,1,1], [1, 1, 1]]) 292 293 assert allclose(quantity.centroid_values, [1, 1, 1, 1]) #Centroid 294 assert allclose(quantity.edge_values, [[1, 1, 1], 295 [1, 1, 1], 296 [1, 1, 1], 297 [1, 1, 1]]) 298 299 300 quantity.set_values(2.0, location = 'centroids') 301 assert allclose(quantity.centroid_values, [2, 2, 2, 2]) 302 303 304 def test_set_values_func(self): 305 quantity = Quantity(self.mesh4) 306 307 def f(x, y): 308 return x+y 309 310 quantity.set_values(f, location = 'vertices') 311 #print "quantity.vertex_values",quantity.vertex_values 312 assert allclose(quantity.vertex_values, 313 [[2,0,2], [2,2,4], [4,2,4], [4,2,4]]) 314 assert allclose(quantity.centroid_values, 315 [4.0/3, 8.0/3, 10.0/3, 10.0/3]) 316 assert allclose(quantity.edge_values, 317 [[1,2,1], [3,3,2], [3,4,3], [3,4,3]]) 318 319 320 quantity.set_values(f, location = 'centroids') 321 assert allclose(quantity.centroid_values, 322 [4.0/3, 8.0/3, 10.0/3, 10.0/3]) 323 324 325 def test_integral(self): 326 quantity = Quantity(self.mesh4) 327 328 #Try constants first 329 const = 5 330 quantity.set_values(const, location = 'vertices') 331 #print 'Q', quantity.get_integral() 332 333 assert allclose(quantity.get_integral(), self.mesh4.get_area() * const) 334 335 #Try with a linear function 336 def f(x, y): 337 return x+y 338 339 quantity.set_values(f, location = 'vertices') 340 341 342 ref_integral = (4.0/3 + 8.0/3 + 10.0/3 + 10.0/3) * 2 343 344 assert allclose (quantity.get_integral(), ref_integral) 345 346 347 348 def test_set_vertex_values(self): 349 quantity = Quantity(self.mesh4) 350 quantity.set_vertex_values([0,1,2,3,4,5]) 351 352 assert allclose(quantity.vertex_values, 353 [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 354 assert allclose(quantity.centroid_values, 355 [1., 7./3, 11./3, 8./3]) #Centroid 356 assert allclose(quantity.edge_values, [[1., 1.5, 0.5], 357 [3., 2.5, 1.5], 358 [3.5, 4.5, 3.], 359 [2.5, 3.5, 2]]) 360 361 362 def test_set_vertex_values_subset(self): 363 quantity = Quantity(self.mesh4) 364 quantity.set_vertex_values([0,1,2,3,4,5]) 365 quantity.set_vertex_values([0,20,30,50], indices = [0,2,3,5]) 366 367 assert allclose(quantity.vertex_values, 368 [[1,0,20], [1,20,4], [4,20,50], [30,1,4]]) 369 370 371 def test_set_vertex_values_using_general_interface(self): 372 quantity = Quantity(self.mesh4) 373 374 375 quantity.set_values([0,1,2,3,4,5]) 376 377 378 assert allclose(quantity.vertex_values, 379 [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 380 381 #Centroid 382 assert allclose(quantity.centroid_values, [1., 7./3, 11./3, 8./3]) 383 384 assert allclose(quantity.edge_values, [[1., 1.5, 0.5], 385 [3., 2.5, 1.5], 386 [3.5, 4.5, 3.], 387 [2.5, 3.5, 2]]) 388 389 390 391 def test_set_vertex_values_using_general_interface_with_subset(self): 392 """test_set_vertex_values_using_general_interface_with_subset(self): 393 394 Test that indices and polygon works (for constants values) 395 """ 396 397 quantity = Quantity(self.mesh4) 398 399 400 quantity.set_values([0,2,3,5], indices=[0,2,3,5]) 401 assert allclose(quantity.vertex_values, 402 [[0,0,2], [0,2,0], [0,2,5], [3,0,0]]) 403 404 405 # Constant 406 quantity.set_values(0.0) 407 quantity.set_values(3.14, indices=[0,2], location='vertices') 408 409 # Indices refer to triangle numbers 410 assert allclose(quantity.vertex_values, 411 [[3.14,3.14,3.14], [0,0,0], 412 [3.14,3.14,3.14], [0,0,0]]) 413 414 415 416 # Now try with polygon (pick points where y>2) 417 polygon = [[0,2.1], [4,2.1], [4,7], [0,7]] 418 quantity.set_values(0.0) 419 quantity.set_values(3.14, polygon=polygon) 420 421 assert allclose(quantity.vertex_values, 422 [[0,0,0], [0,0,0], [0,0,0], 423 [3.14,3.14,3.14]]) 424 425 426 # Another polygon (pick triangle 1 and 2 (rightmost triangles) 427 # using centroids 428 polygon = [[2.1, 0.0], [3.5,0.1], [2,2.2], [0.2,2]] 429 quantity.set_values(0.0) 430 quantity.set_values(3.14, location='centroids', polygon=polygon) 431 assert allclose(quantity.vertex_values, 432 [[0,0,0], 433 [3.14,3.14,3.14], 434 [3.14,3.14,3.14], 435 [0,0,0]]) 436 437 438 # Same polygon now use vertices (default) 439 polygon = [[2.1, 0.0], [3.5,0.1], [2,2.2], [0.2,2]] 440 quantity.set_values(0.0) 441 #print 'Here 2' 442 quantity.set_values(3.14, polygon=polygon) 443 assert allclose(quantity.vertex_values, 444 [[0,0,0], 445 [3.14,3.14,3.14], 446 [3.14,3.14,3.14], 447 [0,0,0]]) 448 449 450 # Test input checking 451 try: 452 quantity.set_values(3.14, polygon=polygon, indices = [0,2]) 453 except: 454 pass 455 else: 456 msg = 'Should have caught this' 457 raise msg 458 459 460 461 462 463 def test_set_vertex_values_using_general_interface_subset_and_geo(self): 464 """test_set_vertex_values_using_general_interface_with_subset(self): 465 Test that indices and polygon works using georeferencing 466 """ 467 468 quantity = Quantity(self.mesh4) 469 G = Geo_reference(56, 10, 100) 470 quantity.domain.geo_reference = G 471 472 #print quantity.domain.get_nodes(absolute=True) 473 474 475 # Constant 476 quantity.set_values(0.0) 477 quantity.set_values(3.14, indices=[0,2], location='vertices') 478 479 # Indices refer to triangle numbers here  not vertices (why?) 480 assert allclose(quantity.vertex_values, 481 [[3.14,3.14,3.14], [0,0,0], 482 [3.14,3.14,3.14], [0,0,0]]) 483 484 485 486 # Now try with polygon (pick points where y>2) 487 polygon = array([[0,2.1], [4,2.1], [4,7], [0,7]]) 488 polygon += [G.xllcorner, G.yllcorner] 489 490 quantity.set_values(0.0) 491 quantity.set_values(3.14, polygon=polygon, location='centroids') 492 493 assert allclose(quantity.vertex_values, 494 [[0,0,0], [0,0,0], [0,0,0], 495 [3.14,3.14,3.14]]) 496 497 498 # Another polygon (pick triangle 1 and 2 (rightmost triangles) 499 polygon = array([[2.1, 0.0], [3.5,0.1], [2,2.2], [0.2,2]]) 500 polygon += [G.xllcorner, G.yllcorner] 501 502 quantity.set_values(0.0) 503 quantity.set_values(3.14, polygon=polygon) 504 505 assert allclose(quantity.vertex_values, 506 [[0,0,0], 507 [3.14,3.14,3.14], 508 [3.14,3.14,3.14], 509 [0,0,0]]) 510 511 512 513 def test_set_values_using_fit(self): 514 515 516 quantity = Quantity(self.mesh4) 517 518 #Get (enough) datapoints 519 data_points = [[ 0.66666667, 0.66666667], 520 [ 1.33333333, 1.33333333], 521 [ 2.66666667, 0.66666667], 522 [ 0.66666667, 2.66666667], 523 [ 0.0, 1.0], 524 [ 0.0, 3.0], 525 [ 1.0, 0.0], 526 [ 1.0, 1.0], 527 [ 1.0, 2.0], 528 [ 1.0, 3.0], 529 [ 2.0, 1.0], 530 [ 3.0, 0.0], 531 [ 3.0, 1.0]] 532 533 z = linear_function(data_points) 534 535 #Use builtin fit_interpolate.fit 536 quantity.set_values( Geospatial_data(data_points, z), alpha = 0 ) 537 #quantity.set_values(points = data_points, values = z, alpha = 0) 538 539 540 answer = linear_function(quantity.domain.get_vertex_coordinates()) 541 #print quantity.vertex_values, answer 542 assert allclose(quantity.vertex_values.flat, answer) 543 544 545 #Now try by setting the same values directly 546 vertex_attributes = fit_to_mesh(data_points, 547 quantity.domain.get_nodes(), 548 quantity.domain.triangles, #FIXME 549 point_attributes=z, 550 alpha = 0, 551 verbose=False) 552 553 #print vertex_attributes 554 quantity.set_values(vertex_attributes) 555 assert allclose(quantity.vertex_values.flat, answer) 556 557 558 559 560 561 def test_test_set_values_using_fit_w_geo(self): 562 563 564 #Mesh 565 vertex_coordinates = [[0.76, 0.76], 566 [0.76, 5.76], 567 [5.76, 0.76]] 568 triangles = [[0,2,1]] 569 570 mesh_georef = Geo_reference(56,0.76,0.76) 571 mesh1 = Domain(vertex_coordinates, triangles, 572 geo_reference = mesh_georef) 573 mesh1.check_integrity() 574 575 #Quantity 576 quantity = Quantity(mesh1) 577 578 #Data 579 data_points = [[ 201.0, 401.0], 580 [ 201.0, 403.0], 581 [ 203.0, 401.0]] 582 583 z = [2, 4, 4] 584 585 data_georef = Geo_reference(56,200,400) 586 587 588 #Reference 589 ref = fit_to_mesh(data_points, vertex_coordinates, triangles, 590 point_attributes=z, 591 data_origin = data_georef.get_origin(), 592 mesh_origin = mesh_georef.get_origin(), 593 alpha = 0) 594 595 assert allclose( ref, [0,5,5] ) 596 597 598 #Test set_values 599 600 quantity.set_values( Geospatial_data(data_points, z, data_georef), alpha = 0 ) 601 602 #quantity.set_values(points = data_points, 603 # values = z, 604 # data_georef = data_georef, 605 # alpha = 0) 606 607 608 #quantity.set_values(points = data_points, 609 # values = z, 610 # data_georef = data_georef, 611 # alpha = 0) 612 assert allclose(quantity.vertex_values.flat, ref) 613 614 615 616 #Test set_values using geospatial data object 617 quantity.vertex_values[:] = 0.0 618 619 geo = Geospatial_data(data_points, z, data_georef) 620 621 622 quantity.set_values(geospatial_data = geo, alpha = 0) 623 assert allclose(quantity.vertex_values.flat, ref) 624 625 626 627 def test_set_values_from_file1(self): 628 quantity = Quantity(self.mesh4) 629 630 #Get (enough) datapoints 631 data_points = [[ 0.66666667, 0.66666667], 632 [ 1.33333333, 1.33333333], 633 [ 2.66666667, 0.66666667], 634 [ 0.66666667, 2.66666667], 635 [ 0.0, 1.0], 636 [ 0.0, 3.0], 637 [ 1.0, 0.0], 638 [ 1.0, 1.0], 639 [ 1.0, 2.0], 640 [ 1.0, 3.0], 641 [ 2.0, 1.0], 642 [ 3.0, 0.0], 643 [ 3.0, 1.0]] 644 645 data_geo_spatial = Geospatial_data(data_points, 646 geo_reference = Geo_reference(56, 0, 0)) 647 data_points_absolute = data_geo_spatial.get_data_points(absolute=True) 648 attributes = linear_function(data_points_absolute) 649 att = 'spam_and_eggs' 650 651 #Create .txt file 652 ptsfile = tempfile.mktemp(".txt") 653 file = open(ptsfile,"w") 654 file.write(" x,y," + att + " \n") 655 for data_point, attribute in map(None, data_points_absolute 656 ,attributes): 657 row = str(data_point[0]) + ',' + str(data_point[1]) \ 658 + ',' + str(attribute) 659 file.write(row + "\n") 660 file.close() 661 662 663 #Check that values can be set from file 664 quantity.set_values(filename = ptsfile, 665 attribute_name = att, alpha = 0) 666 answer = linear_function(quantity.domain.get_vertex_coordinates()) 667 668 #print quantity.vertex_values.flat 669 #print answer 670 671 672 assert allclose(quantity.vertex_values.flat, answer) 673 674 675 #Check that values can be set from file using default attribute 676 quantity.set_values(filename = ptsfile, alpha = 0) 677 assert allclose(quantity.vertex_values.flat, answer) 678 679 #Cleanup 680 import os 681 os.remove(ptsfile) 682 683 684 685 def Xtest_set_values_from_file_using_polygon(self): 686 """test_set_values_from_file_using_polygon(self): 687 688 Test that polygon restriction works for general points data 689 """ 690 691 quantity = Quantity(self.mesh4) 692 693 #Get (enough) datapoints 694 data_points = [[ 0.66666667, 0.66666667], 695 [ 1.33333333, 1.33333333], 696 [ 2.66666667, 0.66666667], 697 [ 0.66666667, 2.66666667], 698 [ 0.0, 1.0], 699 [ 0.0, 3.0], 700 [ 1.0, 0.0], 701 [ 1.0, 1.0], 702 [ 1.0, 2.0], 703 [ 1.0, 3.0], 704 [ 2.0, 1.0], 705 [ 3.0, 0.0], 706 [ 3.0, 1.0]] 707 708 data_geo_spatial = Geospatial_data(data_points, 709 geo_reference = Geo_reference(56, 0, 0)) 710 data_points_absolute = data_geo_spatial.get_data_points(absolute=True) 711 attributes = linear_function(data_points_absolute) 712 att = 'spam_and_eggs' 713 714 #Create .txt file 715 ptsfile = tempfile.mktemp(".txt") 716 file = open(ptsfile,"w") 717 file.write(" x,y," + att + " \n") 718 for data_point, attribute in map(None, data_points_absolute 719 ,attributes): 720 row = str(data_point[0]) + ',' + str(data_point[1]) \ 721 + ',' + str(attribute) 722 file.write(row + "\n") 723 file.close() 724 725 # Create restricting polygon (containing node #4 (2,2) and 726 # centroid of triangle #1 (bce) 727 polygon = [[1.0, 1.0], [4.0, 1.0], 728 [4.0, 4.0], [1.0, 4.0]] 729 730 #print self.mesh4.nodes 731 #print inside_polygon(self.mesh4.nodes, polygon) 732 assert allclose(inside_polygon(self.mesh4.nodes, polygon), 4) 733 734 #print quantity.domain.get_vertex_coordinates() 735 #print quantity.domain.get_nodes() 736 737 # Check that values can be set from file 738 quantity.set_values(filename=ptsfile, 739 polygon=polygon, 740 location='unique vertices', 741 alpha=0) 742 743 # Get indices for vertex coordinates in polygon 744 indices = inside_polygon(quantity.domain.get_vertex_coordinates(), 745 polygon) 746 points = take(quantity.domain.get_vertex_coordinates(), indices) 747 748 answer = linear_function(points) 749 750 #print quantity.vertex_values.flat 751 #print answer 752 753 # Check vertices in polygon have been set 754 assert allclose(take(quantity.vertex_values.flat, indices), 755 answer) 756 757 # Check vertices outside polygon are zero 758 indices = outside_polygon(quantity.domain.get_vertex_coordinates(), 759 polygon) 760 assert allclose(take(quantity.vertex_values.flat, indices), 761 0.0) 762 763 #Cleanup 764 import os 765 os.remove(ptsfile) 766 767 768 769 770 def test_cache_test_set_values_from_file(self): 771 # FIXME (Ole): What is this about? 772 # I don't think it checks anything new 773 quantity = Quantity(self.mesh4) 774 775 #Get (enough) datapoints 776 data_points = [[ 0.66666667, 0.66666667], 777 [ 1.33333333, 1.33333333], 778 [ 2.66666667, 0.66666667], 779 [ 0.66666667, 2.66666667], 780 [ 0.0, 1.0], 781 [ 0.0, 3.0], 782 [ 1.0, 0.0], 783 [ 1.0, 1.0], 784 [ 1.0, 2.0], 785 [ 1.0, 3.0], 786 [ 2.0, 1.0], 787 [ 3.0, 0.0], 788 [ 3.0, 1.0]] 789 790 georef = Geo_reference(56, 0, 0) 791 data_geo_spatial = Geospatial_data(data_points, 792 geo_reference=georef) 793 794 data_points_absolute = data_geo_spatial.get_data_points(absolute=True) 795 attributes = linear_function(data_points_absolute) 796 att = 'spam_and_eggs' 797 798 # Create .txt file 799 ptsfile = tempfile.mktemp(".txt") 800 file = open(ptsfile,"w") 801 file.write(" x,y," + att + " \n") 802 for data_point, attribute in map(None, data_points_absolute 803 ,attributes): 804 row = str(data_point[0]) + ',' + str(data_point[1]) \ 805 + ',' + str(attribute) 806 file.write(row + "\n") 807 file.close() 808 809 810 # Check that values can be set from file 811 quantity.set_values(filename=ptsfile, 812 attribute_name=att, 813 alpha=0, 814 use_cache=True, 815 verbose=False) 816 answer = linear_function(quantity.domain.get_vertex_coordinates()) 817 assert allclose(quantity.vertex_values.flat, answer) 818 819 820 # Check that values can be set from file using default attribute 821 quantity.set_values(filename=ptsfile, 822 alpha=0) 823 assert allclose(quantity.vertex_values.flat, answer) 824 825 # Check cache 826 quantity.set_values(filename=ptsfile, 827 attribute_name=att, 828 alpha=0, 829 use_cache=True, 830 verbose=False) 831 832 833 #Cleanup 834 import os 835 os.remove(ptsfile) 836 837 def test_set_values_from_lat_long(self): 838 quantity = Quantity(self.mesh_onslow) 839 840 #Get (enough) datapoints 841 data_points = [[21.5, 114.5],[21.4, 114.6],[21.45,114.65], 842 [21.35, 114.65],[21.45, 114.55],[21.45,114.6]] 843 844 data_geo_spatial = Geospatial_data(data_points, 845 points_are_lats_longs=True) 846 points_UTM = data_geo_spatial.get_data_points(absolute=True) 847 attributes = linear_function(points_UTM) 848 att = 'elevation' 849 850 #Create .txt file 851 txt_file = tempfile.mktemp(".txt") 852 file = open(txt_file,"w") 853 file.write(" lat,long," + att + " \n") 854 for data_point, attribute in map(None, data_points, attributes): 855 row = str(data_point[0]) + ',' + str(data_point[1]) \ 856 + ',' + str(attribute) 857 #print "row", row 858 file.write(row + "\n") 859 file.close() 860 861 862 #Check that values can be set from file 863 quantity.set_values(filename=txt_file, 864 attribute_name=att, 865 alpha=0) 866 answer = linear_function(quantity.domain.get_vertex_coordinates()) 867 868 #print "quantity.vertex_values.flat", quantity.vertex_values.flat 869 #print "answer",answer 870 871 assert allclose(quantity.vertex_values.flat, answer) 872 873 874 #Check that values can be set from file using default attribute 875 quantity.set_values(filename=txt_file, alpha=0) 876 assert allclose(quantity.vertex_values.flat, answer) 877 878 #Cleanup 879 import os 880 os.remove(txt_file) 881 882 def test_set_values_from_lat_long(self): 883 quantity = Quantity(self.mesh_onslow) 884 885 #Get (enough) datapoints 886 data_points = [[21.5, 114.5],[21.4, 114.6],[21.45,114.65], 887 [21.35, 114.65],[21.45, 114.55],[21.45,114.6]] 888 889 data_geo_spatial = Geospatial_data(data_points, 890 points_are_lats_longs=True) 891 points_UTM = data_geo_spatial.get_data_points(absolute=True) 892 attributes = linear_function(points_UTM) 893 att = 'elevation' 894 895 #Create .txt file 896 txt_file = tempfile.mktemp(".txt") 897 file = open(txt_file,"w") 898 file.write(" lat,long," + att + " \n") 899 for data_point, attribute in map(None, data_points, attributes): 900 row = str(data_point[0]) + ',' + str(data_point[1]) \ 901 + ',' + str(attribute) 902 #print "row", row 903 file.write(row + "\n") 904 file.close() 905 906 907 #Check that values can be set from file 908 quantity.set_values(filename=txt_file, 909 attribute_name=att, alpha=0) 910 answer = linear_function(quantity.domain.get_vertex_coordinates()) 911 912 #print "quantity.vertex_values.flat", quantity.vertex_values.flat 913 #print "answer",answer 914 915 assert allclose(quantity.vertex_values.flat, answer) 916 917 918 #Check that values can be set from file using default attribute 919 quantity.set_values(filename=txt_file, alpha=0) 920 assert allclose(quantity.vertex_values.flat, answer) 921 922 #Cleanup 923 import os 924 os.remove(txt_file) 925 926 def test_set_values_from_UTM_pts(self): 927 quantity = Quantity(self.mesh_onslow) 928 929 #Get (enough) datapoints 930 data_points = [[21.5, 114.5],[21.4, 114.6],[21.45,114.65], 931 [21.35, 114.65],[21.45, 114.55],[21.45,114.6]] 932 933 data_geo_spatial = Geospatial_data(data_points, 934 points_are_lats_longs=True) 935 points_UTM = data_geo_spatial.get_data_points(absolute=True) 936 attributes = linear_function(points_UTM) 937 att = 'elevation' 938 939 #Create .txt file 940 txt_file = tempfile.mktemp(".txt") 941 file = open(txt_file,"w") 942 file.write(" x,y," + att + " \n") 943 for data_point, attribute in map(None, points_UTM, attributes): 944 row = str(data_point[0]) + ',' + str(data_point[1]) \ 945 + ',' + str(attribute) 946 #print "row", row 947 file.write(row + "\n") 948 file.close() 949 950 951 pts_file = tempfile.mktemp(".pts") 952 convert = Geospatial_data(txt_file) 953 convert.export_points_file(pts_file) 954 955 #Check that values can be set from file 956 quantity.set_values_from_file(pts_file, att, 0, 957 'vertices', None) 958 answer = linear_function(quantity.domain.get_vertex_coordinates()) 959 #print "quantity.vertex_values.flat", quantity.vertex_values.flat 960 #print "answer",answer 961 assert allclose(quantity.vertex_values.flat, answer) 962 963 #Check that values can be set from file 964 quantity.set_values(filename=pts_file, 965 attribute_name=att, alpha=0) 966 answer = linear_function(quantity.domain.get_vertex_coordinates()) 967 #print "quantity.vertex_values.flat", quantity.vertex_values.flat 968 #print "answer",answer 969 assert allclose(quantity.vertex_values.flat, answer) 970 971 972 #Check that values can be set from file using default attribute 973 quantity.set_values(filename=txt_file, alpha=0) 974 assert allclose(quantity.vertex_values.flat, answer) 975 976 #Cleanup 977 import os 978 os.remove(txt_file) 979 os.remove(pts_file) 980 981 def verbose_test_set_values_from_UTM_pts(self): 982 quantity = Quantity(self.mesh_onslow) 983 984 #Get (enough) datapoints 985 data_points = [[21.5, 114.5],[21.4, 114.6],[21.45,114.65], 986 [21.35, 114.65],[21.45, 114.55],[21.45,114.6], 987 [21.5, 114.5],[21.4, 114.6],[21.45,114.65], 988 [21.35, 114.65],[21.45, 114.55],[21.45,114.6], 989 [21.5, 114.5],[21.4, 114.6],[21.45,114.65], 990 [21.35, 114.65],[21.45, 114.55],[21.45,114.6], 991 [21.5, 114.5],[21.4, 114.6],[21.45,114.65], 992 [21.35, 114.65],[21.45, 114.55],[21.45,114.6], 993 [21.5, 114.5],[21.4, 114.6],[21.45,114.65], 994 [21.35, 114.65],[21.45, 114.55],[21.45,114.6], 995 [21.5, 114.5],[21.4, 114.6],[21.45,114.65], 996 [21.35, 114.65],[21.45, 114.55],[21.45,114.6], 997 [21.5, 114.5],[21.4, 114.6],[21.45,114.65], 998 [21.35, 114.65],[21.45, 114.55],[21.45,114.6], 999 [21.35, 114.65],[21.45, 114.55],[21.45,114.6], 1000 [21.5, 114.5],[21.4, 114.6],[21.45,114.65], 1001 [21.35, 114.65],[21.45, 114.55],[21.45,114.6], 1002 [21.5, 114.5],[21.4, 114.6],[21.45,114.65], 1003 [21.35, 114.65],[21.45, 114.55],[21.45,114.6], 1004 [21.5, 114.5],[21.4, 114.6],[21.45,114.65], 1005 [21.35, 114.65],[21.45, 114.55],[21.45,114.6], 1006 [21.5, 114.5],[21.4, 114.6],[21.45,114.65], 1007 [21.35, 114.65],[21.45, 114.55],[21.45,114.6], 1008 [21.5, 114.5],[21.4, 114.6],[21.45,114.65], 1009 [21.35, 114.65],[21.45, 114.55],[21.45,114.6], 1010 [21.5, 114.5],[21.4, 114.6],[21.45,114.65], 1011 [21.35, 114.65],[21.45, 114.55],[21.45,114.6], 1012 ] 1013 1014 data_geo_spatial = Geospatial_data(data_points, 1015 points_are_lats_longs=True) 1016 points_UTM = data_geo_spatial.get_data_points(absolute=True) 1017 attributes = linear_function(points_UTM) 1018 att = 'elevation' 1019 1020 #Create .txt file 1021 txt_file = tempfile.mktemp(".txt") 1022 file = open(txt_file,"w") 1023 file.write(" x,y," + att + " \n") 1024 for data_point, attribute in map(None, points_UTM, attributes): 1025 row = str(data_point[0]) + ',' + str(data_point[1]) \ 1026 + ',' + str(attribute) 1027 #print "row", row 1028 file.write(row + "\n") 1029 file.close() 1030 1031 1032 pts_file = tempfile.mktemp(".pts") 1033 convert = Geospatial_data(txt_file) 1034 convert.export_points_file(pts_file) 1035 1036 #Check that values can be set from file 1037 quantity.set_values_from_file(pts_file, att, 0, 1038 'vertices', None, verbose = True, 1039 max_read_lines=2) 1040 answer = linear_function(quantity.domain.get_vertex_coordinates()) 1041 #print "quantity.vertex_values.flat", quantity.vertex_values.flat 1042 #print "answer",answer 1043 assert allclose(quantity.vertex_values.flat, answer) 1044 1045 #Check that values can be set from file 1046 quantity.set_values(filename=pts_file, 1047 attribute_name=att, alpha=0) 1048 answer = linear_function(quantity.domain.get_vertex_coordinates()) 1049 #print "quantity.vertex_values.flat", quantity.vertex_values.flat 1050 #print "answer",answer 1051 assert allclose(quantity.vertex_values.flat, answer) 1052 1053 1054 #Check that values can be set from file using default attribute 1055 quantity.set_values(filename=txt_file, alpha=0) 1056 assert allclose(quantity.vertex_values.flat, answer) 1057 1058 #Cleanup 1059 import os 1060 os.remove(txt_file) 1061 os.remove(pts_file) 1062 1063 def test_set_values_from_file_with_georef1(self): 1064 1065 #Mesh in zone 56 (absolute coords) 1066 1067 x0 = 314036.58727982 1068 y0 = 6224951.2960092 1069 1070 a = [x0+0.0, y0+0.0] 1071 b = [x0+0.0, y0+2.0] 1072 c = [x0+2.0, y0+0.0] 1073 d = [x0+0.0, y0+4.0] 1074 e = [x0+2.0, y0+2.0] 1075 f = [x0+4.0, y0+0.0] 1076 1077 points = [a, b, c, d, e, f] 1078 1079 #bac, bce, ecf, dbe 1080 elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ] 1081 1082 #absolute going in .. 1083 mesh4 = Domain(points, elements, 1084 geo_reference = Geo_reference(56, 0, 0)) 1085 mesh4.check_integrity() 1086 quantity = Quantity(mesh4) 1087 1088 #Get (enough) datapoints (relative to georef) 1089 data_points_rel = [[ 0.66666667, 0.66666667], 1090 [ 1.33333333, 1.33333333], 1091 [ 2.66666667, 0.66666667], 1092 [ 0.66666667, 2.66666667], 1093 [ 0.0, 1.0], 1094 [ 0.0, 3.0], 1095 [ 1.0, 0.0], 1096 [ 1.0, 1.0], 1097 [ 1.0, 2.0], 1098 [ 1.0, 3.0], 1099 [ 2.0, 1.0], 1100 [ 3.0, 0.0], 1101 [ 3.0, 1.0]] 1102 1103 data_geo_spatial = Geospatial_data(data_points_rel, 1104 geo_reference = Geo_reference(56, x0, y0)) 1105 data_points_absolute = data_geo_spatial.get_data_points(absolute=True) 1106 attributes = linear_function(data_points_absolute) 1107 att = 'spam_and_eggs' 1108 1109 #Create .txt file 1110 ptsfile = tempfile.mktemp(".txt") 1111 file = open(ptsfile,"w") 1112 file.write(" x,y," + att + " \n") 1113 for data_point, attribute in map(None, data_points_absolute 1114 ,attributes): 1115 row = str(data_point[0]) + ',' + str(data_point[1]) \ 1116 + ',' + str(attribute) 1117 file.write(row + "\n") 1118 file.close() 1119 1120 #file = open(ptsfile, 'r') 1121 #lines = file.readlines() 1122 #file.close() 1123 1124 1125 #Check that values can be set from file 1126 quantity.set_values(filename=ptsfile, 1127 attribute_name=att, alpha=0) 1128 answer = linear_function(quantity.domain.get_vertex_coordinates()) 1129 1130 assert allclose(quantity.vertex_values.flat, answer) 1131 1132 1133 #Check that values can be set from file using default attribute 1134 quantity.set_values(filename=ptsfile, alpha=0) 1135 assert allclose(quantity.vertex_values.flat, answer) 1136 1137 #Cleanup 1138 import os 1139 os.remove(ptsfile) 1140 1141 1142 def test_set_values_from_file_with_georef2(self): 1143 1144 #Mesh in zone 56 (relative coords) 1145 1146 x0 = 314036.58727982 1147 y0 = 6224951.2960092 1148 #x0 = 0.0 1149 #y0 = 0.0 1150 1151 a = [0.0, 0.0] 1152 b = [0.0, 2.0] 1153 c = [2.0, 0.0] 1154 d = [0.0, 4.0] 1155 e = [2.0, 2.0] 1156 f = [4.0, 0.0] 1157 1158 points = [a, b, c, d, e, f] 1159 1160 #bac, bce, ecf, dbe 1161 elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ] 1162 1163 mesh4 = Domain(points, elements, 1164 geo_reference = Geo_reference(56, x0, y0)) 1165 mesh4.check_integrity() 1166 quantity = Quantity(mesh4) 1167 1168 #Get (enough) datapoints 1169 data_points = [[ x0+0.66666667, y0+0.66666667], 1170 [ x0+1.33333333, y0+1.33333333], 1171 [ x0+2.66666667, y0+0.66666667], 1172 [ x0+0.66666667, y0+2.66666667], 1173 [ x0+0.0, y0+1.0], 1174 [ x0+0.0, y0+3.0], 1175 [ x0+1.0, y0+0.0], 1176 [ x0+1.0, y0+1.0], 1177 [ x0+1.0, y0+2.0], 1178 [ x0+1.0, y0+3.0], 1179 [ x0+2.0, y0+1.0], 1180 [ x0+3.0, y0+0.0], 1181 [ x0+3.0, y0+1.0]] 1182 1183 1184 data_geo_spatial = Geospatial_data(data_points, 1185 geo_reference = Geo_reference(56, 0, 0)) 1186 data_points_absolute = data_geo_spatial.get_data_points(absolute=True) 1187 attributes = linear_function(data_points_absolute) 1188 att = 'spam_and_eggs' 1189 1190 #Create .txt file 1191 ptsfile = tempfile.mktemp(".txt") 1192 file = open(ptsfile,"w") 1193 file.write(" x,y," + att + " \n") 1194 for data_point, attribute in map(None, data_points_absolute 1195 ,attributes): 1196 row = str(data_point[0]) + ',' + str(data_point[1]) \ 1197 + ',' + str(attribute) 1198 file.write(row + "\n") 1199 file.close() 1200 1201 1202 #Check that values can be set from file 1203 quantity.set_values(filename=ptsfile, 1204 attribute_name=att, alpha=0) 1205 answer = linear_function(quantity.domain. \ 1206 get_vertex_coordinates(absolute=True)) 1207 1208 1209 assert allclose(quantity.vertex_values.flat, answer) 1210 1211 1212 #Check that values can be set from file using default attribute 1213 quantity.set_values(filename=ptsfile, alpha=0) 1214 assert allclose(quantity.vertex_values.flat, answer) 1215 1216 #Cleanup 1217 import os 1218 os.remove(ptsfile) 1219 1220 1221 1222 1223 def test_set_values_from_quantity(self): 1224 1225 quantity1 = Quantity(self.mesh4) 1226 quantity1.set_vertex_values([0,1,2,3,4,5]) 1227 1228 assert allclose(quantity1.vertex_values, 1229 [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 1230 1231 1232 quantity2 = Quantity(self.mesh4) 1233 quantity2.set_values(quantity=quantity1) 1234 assert allclose(quantity2.vertex_values, 1235 [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 1236 1237 quantity2.set_values(quantity = 2*quantity1) 1238 assert allclose(quantity2.vertex_values, 1239 [[2,0,4], [2,4,8], [8,4,10], [6,2,8]]) 1240 1241 quantity2.set_values(quantity = 2*quantity1 + 3) 1242 assert allclose(quantity2.vertex_values, 1243 [[5,3,7], [5,7,11], [11,7,13], [9,5,11]]) 1244 1245 1246 #Check detection of quantity as first orgument 1247 quantity2.set_values(2*quantity1 + 3) 1248 assert allclose(quantity2.vertex_values, 1249 [[5,3,7], [5,7,11], [11,7,13], [9,5,11]]) 1250 1251 1252 1253 def Xtest_set_values_from_quantity_using_polygon(self): 1254 """test_set_values_from_quantity_using_polygon(self): 1255 1256 Check that polygon can be used to restrict set_values when 1257 using another quantity as argument. 1258 """ 1259 1260 # Create restricting polygon (containing node #4 (2,2) and 1261 # centroid of triangle #1 (bce) 1262 polygon = [[1.0, 1.0], [4.0, 1.0], 1263 [4.0, 4.0], [1.0, 4.0]] 1264 assert allclose(inside_polygon(self.mesh4.nodes, polygon), 4) 1265 1266 quantity1 = Quantity(self.mesh4) 1267 quantity1.set_vertex_values([0,1,2,3,4,5]) 1268 1269 assert allclose(quantity1.vertex_values, 1270 [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 1271 1272 1273 quantity2 = Quantity(self.mesh4) 1274 quantity2.set_values(quantity=quantity1, 1275 polygon=polygon) 1276 1277 msg = 'Only node #4(e) at (2,2) should have values applied ' 1278 assert allclose(quantity2.vertex_values, 1279 [[0,0,0], [0,0,4], [4,0,0], [0,0,4]]), msg 1280 #bac, bce, ecf, dbe 1281 1282 1283 1284 def test_overloading(self): 1285 1286 quantity1 = Quantity(self.mesh4) 1287 quantity1.set_vertex_values([0,1,2,3,4,5]) 1288 1289 assert allclose(quantity1.vertex_values, 1290 [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 1291 1292 1293 quantity2 = Quantity(self.mesh4) 1294 quantity2.set_values([[1,2,3], [5,5,5], [0,0,9], [6, 3, 3]], 1295 location = 'vertices') 1296 1297 1298 1299 quantity3 = Quantity(self.mesh4) 1300 quantity3.set_values([[2,2,2], [7,8,9], [7,6,3], [3, 8, 8]], 1301 location = 'vertices') 1302 1303 1304 # Negation 1305 Q = quantity1 1306 assert allclose(Q.vertex_values, quantity1.vertex_values) 1307 assert allclose(Q.centroid_values, quantity1.centroid_values) 1308 assert allclose(Q.edge_values, quantity1.edge_values) 1309 1310 # Addition 1311 Q = quantity1 + 7 1312 assert allclose(Q.vertex_values, quantity1.vertex_values + 7) 1313 assert allclose(Q.centroid_values, quantity1.centroid_values + 7) 1314 assert allclose(Q.edge_values, quantity1.edge_values + 7) 1315 1316 Q = 7 + quantity1 1317 assert allclose(Q.vertex_values, quantity1.vertex_values + 7) 1318 assert allclose(Q.centroid_values, quantity1.centroid_values + 7) 1319 assert allclose(Q.edge_values, quantity1.edge_values + 7) 1320 1321 Q = quantity1 + quantity2 1322 assert allclose(Q.vertex_values, 1323 quantity1.vertex_values + quantity2.vertex_values) 1324 assert allclose(Q.centroid_values, 1325 quantity1.centroid_values + quantity2.centroid_values) 1326 assert allclose(Q.edge_values, 1327 quantity1.edge_values + quantity2.edge_values) 1328 1329 1330 Q = quantity1 + quantity2  3 1331 assert allclose(Q.vertex_values, 1332 quantity1.vertex_values + quantity2.vertex_values  3) 1333 1334 Q = quantity1  quantity2 1335 assert allclose(Q.vertex_values, 1336 quantity1.vertex_values  quantity2.vertex_values) 1337 1338 #Scaling 1339 Q = quantity1*3 1340 assert allclose(Q.vertex_values, quantity1.vertex_values*3) 1341 assert allclose(Q.centroid_values, quantity1.centroid_values*3) 1342 assert allclose(Q.edge_values, quantity1.edge_values*3) 1343 Q = 3*quantity1 1344 assert allclose(Q.vertex_values, quantity1.vertex_values*3) 1345 1346 #Multiplication 1347 Q = quantity1 * quantity2 1348 #print Q.vertex_values 1349 #print Q.centroid_values 1350 #print quantity1.centroid_values 1351 #print quantity2.centroid_values 1352 1353 assert allclose(Q.vertex_values, 1354 quantity1.vertex_values * quantity2.vertex_values) 1355 1356 #Linear combinations 1357 Q = 4*quantity1 + 2 1358 assert allclose(Q.vertex_values, 1359 4*quantity1.vertex_values + 2) 1360 1361 Q = quantity1*quantity2 + 2 1362 assert allclose(Q.vertex_values, 1363 quantity1.vertex_values * quantity2.vertex_values + 2) 1364 1365 Q = quantity1*quantity2 + quantity3 1366 assert allclose(Q.vertex_values, 1367 quantity1.vertex_values * quantity2.vertex_values + 1368 quantity3.vertex_values) 1369 Q = quantity1*quantity2 + 3*quantity3 1370 assert allclose(Q.vertex_values, 1371 quantity1.vertex_values * quantity2.vertex_values + 1372 3*quantity3.vertex_values) 1373 Q = quantity1*quantity2 + 3*quantity3 + 5.0 1374 assert allclose(Q.vertex_values, 1375 quantity1.vertex_values * quantity2.vertex_values + 1376 3*quantity3.vertex_values + 5) 1377 1378 Q = quantity1*quantity2  quantity3 1379 assert allclose(Q.vertex_values, 1380 quantity1.vertex_values * quantity2.vertex_values  1381 quantity3.vertex_values) 1382 Q = 1.5*quantity1*quantity2  3*quantity3 + 5.0 1383 assert allclose(Q.vertex_values, 1384 1.5*quantity1.vertex_values * quantity2.vertex_values  1385 3*quantity3.vertex_values + 5) 1386 1387 #Try combining quantities and arrays and scalars 1388 Q = 1.5*quantity1*quantity2.vertex_values \ 1389 3*quantity3.vertex_values + 5.0 1390 assert allclose(Q.vertex_values, 1391 1.5*quantity1.vertex_values * quantity2.vertex_values  1392 3*quantity3.vertex_values + 5) 1393 1394 1395 #Powers 1396 Q = quantity1**2 1397 assert allclose(Q.vertex_values, quantity1.vertex_values**2) 1398 1399 Q = quantity1**2 +quantity2**2 1400 assert allclose(Q.vertex_values, 1401 quantity1.vertex_values**2 + \ 1402 quantity2.vertex_values**2) 1403 1404 Q = (quantity1**2 +quantity2**2)**0.5 1405 assert allclose(Q.vertex_values, 1406 (quantity1.vertex_values**2 + \ 1407 quantity2.vertex_values**2)**0.5) 1408 1409 1410 1411 1412 1413 1414 1415 def test_compute_gradient(self): 1416 quantity = Quantity(self.mesh4) 1417 1418 #Set up for a gradient of (2,0) at mid triangle 1419 quantity.set_values([2.0, 4.0, 6.0, 2.0], 1420 location = 'centroids') 1421 1422 1423 #Gradients 1424 quantity.compute_gradients() 1425 1426 a = quantity.x_gradient 1427 b = quantity.y_gradient 1428 #print self.mesh4.centroid_coordinates 1429 #print a, b 1430 1431 #The central triangle (1) 1432 #(using standard gradient based on neigbours controid values) 1433 assert allclose(a[1], 2.0) 1434 assert allclose(b[1], 0.0) 1435 1436 1437 #Left triangle (0) using two point gradient 1438 #q0 = q1 + a*(x0x1) + b*(y0y1) <=> 1439 #2 = 4 + a*(2/3) + b*(2/3) 1440 assert allclose(a[0] + b[0], 3) 1441 #From orthogonality (a*(y0y1) + b*(x0x1) == 0) 1442 assert allclose(a[0]  b[0], 0) 1443 1444 1445 #Right triangle (2) using two point gradient 1446 #q2 = q1 + a*(x2x1) + b*(y2y1) <=> 1447 #6 = 4 + a*(4/3) + b*(2/3) 1448 assert allclose(2*a[2]  b[2], 3) 1449 #From orthogonality (a*(y1y2) + b*(x2x1) == 0) 1450 assert allclose(a[2] + 2*b[2], 0) 1451 1452 1453 #Top triangle (3) using two point gradient 1454 #q3 = q1 + a*(x3x1) + b*(y3y1) <=> 1455 #2 = 4 + a*(2/3) + b*(4/3) 1456 assert allclose(a[3]  2*b[3], 3) 1457 #From orthogonality (a*(y1y3) + b*(x3x1) == 0) 1458 assert allclose(2*a[3] + b[3], 0) 1459 1460 1461 1462 #print a, b 1463 quantity.extrapolate_second_order() 1464 1465 #Apply q(x,y) = qc + a*(xxc) + b*(yyc) 1466 assert allclose(quantity.vertex_values[0,:], [3., 0., 3.]) 1467 assert allclose(quantity.vertex_values[1,:], [4./3, 16./3, 16./3]) 1468 1469 1470 #a = 1.2, b=0.6 1471 #q(4,0) = 6 + a*(4  8/3) + b*(2/3) 1472 assert allclose(quantity.vertex_values[2,2], 8) 1473 1474 def test_get_gradients(self): 1475 quantity = Quantity(self.mesh4) 1476 1477 #Set up for a gradient of (2,0) at mid triangle 1478 quantity.set_values([2.0, 4.0, 6.0, 2.0], 1479 location = 'centroids') 1480 1481 1482 #Gradients 1483 quantity.compute_gradients() 1484 1485 a, b = quantity.get_gradients() 1486 #print self.mesh4.centroid_coordinates 1487 #print a, b 1488 1489 #The central triangle (1) 1490 #(using standard gradient based on neigbours controid values) 1491 assert allclose(a[1], 2.0) 1492 assert allclose(b[1], 0.0) 1493 1494 1495 #Left triangle (0) using two point gradient 1496 #q0 = q1 + a*(x0x1) + b*(y0y1) <=> 1497 #2 = 4 + a*(2/3) + b*(2/3) 1498 assert allclose(a[0] + b[0], 3) 1499 #From orthogonality (a*(y0y1) + b*(x0x1) == 0) 1500 assert allclose(a[0]  b[0], 0) 1501 1502 1503 #Right triangle (2) using two point gradient 1504 #q2 = q1 + a*(x2x1) + b*(y2y1) <=> 1505 #6 = 4 + a*(4/3) + b*(2/3) 1506 assert allclose(2*a[2]  b[2], 3) 1507 #From orthogonality (a*(y1y2) + b*(x2x1) == 0) 1508 assert allclose(a[2] + 2*b[2], 0) 1509 1510 1511 #Top triangle (3) using two point gradient 1512 #q3 = q1 + a*(x3x1) + b*(y3y1) <=> 1513 #2 = 4 + a*(2/3) + b*(4/3) 1514 assert allclose(a[3]  2*b[3], 3) 1515 #From orthogonality (a*(y1y3) + b*(x3x1) == 0) 1516 assert allclose(2*a[3] + b[3], 0) 1517 1518 1519 def test_second_order_extrapolation2(self): 1520 quantity = Quantity(self.mesh4) 1521 1522 #Set up for a gradient of (3,1), f(x) = 3x+y 1523 quantity.set_values([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3], 1524 location = 'centroids') 1525 1526 #Gradients 1527 quantity.compute_gradients() 1528 1529 a = quantity.x_gradient 1530 b = quantity.y_gradient 1531 1532 #print a, b 1533 1534 assert allclose(a[1], 3.0) 1535 assert allclose(b[1], 1.0) 1536 1537 #Work out the others 1538 1539 quantity.extrapolate_second_order() 1540 1541 #print quantity.vertex_values 1542 assert allclose(quantity.vertex_values[1,0], 2.0) 1543 assert allclose(quantity.vertex_values[1,1], 6.0) 1544 assert allclose(quantity.vertex_values[1,2], 8.0) 1545 1546 1547 1548 def test_backup_saxpy_centroid_values(self): 1549 quantity = Quantity(self.mesh4) 1550 1551 #Set up for a gradient of (3,1), f(x) = 3x+y 1552 c_values = array([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3]) 1553 d_values = array([1.0, 2.0, 3.0, 4.0]) 1554 quantity.set_values(c_values, location = 'centroids') 1555 1556 #Backup 1557 quantity.backup_centroid_values() 1558 1559 #print quantity.vertex_values 1560 assert allclose(quantity.centroid_values, quantity.centroid_backup_values) 1561 1562 1563 quantity.set_values(d_values, location = 'centroids') 1564 1565 quantity.saxpy_centroid_values(2.0, 3.0) 1566 1567 assert(quantity.centroid_values, 2.0*d_values + 3.0*c_values) 1568 1569 1570 1571 def test_first_order_extrapolator(self): 1572 quantity = Quantity(self.mesh4) 1573 1574 #Test centroids 1575 quantity.set_values([1.,2.,3.,4.], location = 'centroids') 1576 assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid 1577 1578 #Extrapolate 1579 quantity.extrapolate_first_order() 1580 1581 #Check that gradient is zero 1582 a,b = quantity.get_gradients() 1583 assert allclose(a, [0,0,0,0]) 1584 assert allclose(b, [0,0,0,0]) 1585 1586 #Check vertices but not edge values 1587 assert allclose(quantity.vertex_values, 1588 [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]]) 1589 1590 1591 def test_second_order_extrapolator(self): 1592 quantity = Quantity(self.mesh4) 1593 1594 #Set up for a gradient of (3,0) at mid triangle 1595 quantity.set_values([2.0, 4.0, 8.0, 2.0], 1596 location = 'centroids') 1597 1598 1599 1600 quantity.extrapolate_second_order() 1601 quantity.limit() 1602 1603 1604 #Assert that central triangle is limited by neighbours 1605 assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0] 1606 assert quantity.vertex_values[1,0] >= quantity.vertex_values[3,1] 1607 1608 assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1] 1609 assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2] 1610 1611 assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0] 1612 assert quantity.vertex_values[1,2] >= quantity.vertex_values[3,1] 1613 1614 1615 #Assert that quantities are conserved 1616 from Numeric import sum 1617 for k in range(quantity.centroid_values.shape[0]): 1618 assert allclose (quantity.centroid_values[k], 1619 sum(quantity.vertex_values[k,:])/3) 1620 1621 1622 1623 1624 1625 def test_limit_vertices_by_all_neighbours(self): 1626 quantity = Quantity(self.mesh4) 1627 1628 #Create a deliberate overshoot (e.g. from gradient computation) 1629 quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]]) 1630 1631 1632 #Limit 1633 quantity.limit_vertices_by_all_neighbours() 1634 1635 #Assert that central triangle is limited by neighbours 1636 assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0] 1637 assert quantity.vertex_values[1,0] <= quantity.vertex_values[3,1] 1638 1639 assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1] 1640 assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2] 1641 1642 assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0] 1643 assert quantity.vertex_values[1,2] <= quantity.vertex_values[3,1] 1644 1645 1646 1647 #Assert that quantities are conserved 1648 from Numeric import sum 1649 for k in range(quantity.centroid_values.shape[0]): 1650 assert allclose (quantity.centroid_values[k], 1651 sum(quantity.vertex_values[k,:])/3) 1652 1653 1654 1655 def test_limit_edges_by_all_neighbours(self): 1656 quantity = Quantity(self.mesh4) 1657 1658 #Create a deliberate overshoot (e.g. from gradient computation) 1659 quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]]) 1660 1661 1662 #Limit 1663 quantity.limit_edges_by_all_neighbours() 1664 1665 #Assert that central triangle is limited by neighbours 1666 assert quantity.edge_values[1,0] <= quantity.centroid_values[2] 1667 assert quantity.edge_values[1,0] >= quantity.centroid_values[0] 1668 1669 assert quantity.edge_values[1,1] <= quantity.centroid_values[2] 1670 assert quantity.edge_values[1,1] >= quantity.centroid_values[0] 1671 1672 assert quantity.edge_values[1,2] <= quantity.centroid_values[2] 1673 assert quantity.edge_values[1,2] >= quantity.centroid_values[0] 1674 1675 1676 1677 #Assert that quantities are conserved 1678 from Numeric import sum 1679 for k in range(quantity.centroid_values.shape[0]): 1680 assert allclose (quantity.centroid_values[k], 1681 sum(quantity.vertex_values[k,:])/3) 1682 1683 1684 def test_limit_edges_by_neighbour(self): 1685 quantity = Quantity(self.mesh4) 1686 1687 #Create a deliberate overshoot (e.g. from gradient computation) 1688 quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]]) 1689 1690 1691 #Limit 1692 quantity.limit_edges_by_neighbour() 1693 1694 #Assert that central triangle is limited by neighbours 1695 assert quantity.edge_values[1,0] <= quantity.centroid_values[3] 1696 assert quantity.edge_values[1,0] >= quantity.centroid_values[1] 1697 1698 assert quantity.edge_values[1,1] <= quantity.centroid_values[2] 1699 assert quantity.edge_values[1,1] >= quantity.centroid_values[1] 1700 1701 assert quantity.edge_values[1,2] <= quantity.centroid_values[1] 1702 assert quantity.edge_values[1,2] >= quantity.centroid_values[0] 1703 1704 1705 1706 #Assert that quantities are conserved 1707 from Numeric import sum 1708 for k in range(quantity.centroid_values.shape[0]): 1709 assert allclose (quantity.centroid_values[k], 1710 sum(quantity.vertex_values[k,:])/3) 1711 1712 def test_limiter2(self): 1713 """Taken from test_shallow_water 1714 """ 1715 quantity = Quantity(self.mesh4) 1716 quantity.domain.beta_w = 0.9 1717 1718 #Test centroids 1719 quantity.set_values([2.,4.,8.,2.], location = 'centroids') 1720 assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid 1721 1722 1723 #Extrapolate 1724 quantity.extrapolate_second_order() 1725 1726 assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6]) 1727 1728 #Limit 1729 quantity.limit() 1730 1731 # limited value for beta_w = 0.9 1732 1733 assert allclose(quantity.vertex_values[1,:], [2.2, 4.9, 4.9]) 1734 # limited values for beta_w = 0.5 1735 #assert allclose(quantity.vertex_values[1,:], [3.0, 4.5, 4.5]) 1736 1737 1738 #Assert that quantities are conserved 1739 from Numeric import sum 1740 for k in range(quantity.centroid_values.shape[0]): 1741 assert allclose (quantity.centroid_values[k], 1742 sum(quantity.vertex_values[k,:])/3) 1743 1744 1745 1746 1747 1748 def test_distribute_first_order(self): 1749 quantity = Quantity(self.mesh4) 1750 1751 #Test centroids 1752 quantity.set_values([1.,2.,3.,4.], location = 'centroids') 1753 assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid 1754 1755 1756 #Extrapolate from centroid to vertices and edges 1757 quantity.extrapolate_first_order() 1758 1759 #Interpolate 1760 #quantity.interpolate_from_vertices_to_edges() 1761 1762 assert allclose(quantity.vertex_values, 1763 [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]]) 1764 assert allclose(quantity.edge_values, [[1,1,1], [2,2,2], 1765 [3,3,3], [4, 4, 4]]) 1766 1767 1768 def test_interpolate_from_vertices_to_edges(self): 1769 quantity = Quantity(self.mesh4) 1770 1771 quantity.vertex_values = array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]],Float) 1772 1773 quantity.interpolate_from_vertices_to_edges() 1774 1775 assert allclose(quantity.edge_values, [[1., 1.5, 0.5], 1776 [3., 2.5, 1.5], 1777 [3.5, 4.5, 3.], 1778 [2.5, 3.5, 2]]) 1779 1780 1781 def test_interpolate_from_edges_to_vertices(self): 1782 quantity = Quantity(self.mesh4) 1783 1784 quantity.edge_values = array([[1., 1.5, 0.5], 1785 [3., 2.5, 1.5], 1786 [3.5, 4.5, 3.], 1787 [2.5, 3.5, 2]],Float) 1788 1789 quantity.interpolate_from_edges_to_vertices() 1790 1791 assert allclose(quantity.vertex_values, 1792 [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 1793 1794 1795 1796 def test_distribute_second_order(self): 1797 quantity = Quantity(self.mesh4) 1798 1799 #Test centroids 1800 quantity.set_values([2.,4.,8.,2.], location = 'centroids') 1801 assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid 1802 1803 1804 #Extrapolate 1805 quantity.extrapolate_second_order() 1806 1807 assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6]) 1808 1809 1810 def test_update_explicit(self): 1811 quantity = Quantity(self.mesh4) 1812 1813 #Test centroids 1814 quantity.set_values([1.,2.,3.,4.], location = 'centroids') 1815 assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid 1816 1817 #Set explicit_update 1818 quantity.explicit_update = array( [1.,1.,1.,1.] ) 1819 1820 #Update with given timestep 1821 quantity.update(0.1) 1822 1823 x = array([1, 2, 3, 4]) + array( [.1,.1,.1,.1] ) 1824 assert allclose( quantity.centroid_values, x) 1825 1826 def test_update_semi_implicit(self): 1827 quantity = Quantity(self.mesh4) 1828 1829 #Test centroids 1830 quantity.set_values([1.,2.,3.,4.], location = 'centroids') 1831 assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid 1832 1833 #Set semi implicit update 1834 quantity.semi_implicit_update = array([1.,1.,1.,1.]) 1835 1836 #Update with given timestep 1837 timestep = 0.1 1838 quantity.update(timestep) 1839 1840 sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4]) 1841 denom = ones(4, Float)timestep*sem 1842 1843 x = array([1, 2, 3, 4])/denom 1844 assert allclose( quantity.centroid_values, x) 1845 1846 1847 def test_both_updates(self): 1848 quantity = Quantity(self.mesh4) 1849 1850 #Test centroids 1851 quantity.set_values([1.,2.,3.,4.], location = 'centroids') 1852 assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid 1853 1854 #Set explicit_update 1855 quantity.explicit_update = array( [4.,3.,2.,1.] ) 1856 1857 #Set semi implicit update 1858 quantity.semi_implicit_update = array( [1.,1.,1.,1.] ) 1859 1860 #Update with given timestep 1861 timestep = 0.1 1862 quantity.update(0.1) 1863 1864 sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4]) 1865 denom = ones(4, Float)timestep*sem 1866 1867 x = array([1., 2., 3., 4.]) 1868 x /= denom 1869 x += timestep*array( [4.0, 3.0, 2.0, 1.0] ) 1870 1871 assert allclose( quantity.centroid_values, x) 1872 1873 1874 1875 1876 #Test smoothing 1877 def test_smoothing(self): 1878 1879 from mesh_factory import rectangular 1880 from shallow_water import Domain, Transmissive_boundary 1881 from Numeric import zeros, Float 1882 from anuga.utilities.numerical_tools import mean 1883 1884 #Create basic mesh 1885 points, vertices, boundary = rectangular(2, 2) 1886 1887 #Create shallow water domain 1888 domain = Domain(points, vertices, boundary) 1889 domain.default_order=2 1890 domain.reduction = mean 1891 1892 1893 #Set some field values 1894 domain.set_quantity('elevation', lambda x,y: x) 1895 domain.set_quantity('friction', 0.03) 1896 1897 1898 ###################### 1899 # Boundary conditions 1900 B = Transmissive_boundary(domain) 1901 domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B}) 1902 1903 1904 ###################### 1905 #Initial condition  with jumps 1906 1907 bed = domain.quantities['elevation'].vertex_values 1908 stage = zeros(bed.shape, Float) 1909 1910 h = 0.03 1911 for i in range(stage.shape[0]): 1912 if i % 2 == 0: 1913 stage[i,:] = bed[i,:] + h 1914 else: 1915 stage[i,:] = bed[i,:] 1916 1917 domain.set_quantity('stage', stage) 1918 1919 stage = domain.quantities['stage'] 1920 1921 #Get smoothed stage 1922 A, V = stage.get_vertex_values(xy=False, smooth=True) 1923 Q = stage.vertex_values 1924 1925 1926 assert A.shape[0] == 9 1927 assert V.shape[0] == 8 1928 assert V.shape[1] == 3 1929 1930 #First four points 1931 assert allclose(A[0], (Q[0,2] + Q[1,1])/2) 1932 assert allclose(A[1], (Q[1,0] + Q[3,1] + Q[2,2])/3) 1933 assert allclose(A[2], Q[3,0]) 1934 assert allclose(A[3], (Q[0,0] + Q[5,1] + Q[4,2])/3) 1935 1936 #Center point 1937 assert allclose(A[4], (Q[0,1] + Q[1,2] + Q[2,0] +\ 1938 Q[5,0] + Q[6,2] + Q[7,1])/6) 1939 1940 1941 #Check V 1942 assert allclose(V[0,:], [3,4,0]) 1943 assert allclose(V[1,:], [1,0,4]) 1944 assert allclose(V[2,:], [4,5,1]) 1945 assert allclose(V[3,:], [2,1,5]) 1946 assert allclose(V[4,:], [6,7,3]) 1947 assert allclose(V[5,:], [4,3,7]) 1948 assert allclose(V[6,:], [7,8,4]) 1949 assert allclose(V[7,:], [5,4,8]) 1950 1951 #Get smoothed stage with XY 1952 X, Y, A1, V1 = stage.get_vertex_values(xy=True, smooth=True) 1953 1954 assert allclose(A, A1) 1955 assert allclose(V, V1) 1956 1957 #Check XY 1958 assert allclose(X[4], 0.5) 1959 assert allclose(Y[4], 0.5) 1960 1961 assert allclose(X[7], 1.0) 1962 assert allclose(Y[7], 0.5) 1963 1964 1965 1966 1967 def test_vertex_values_no_smoothing(self): 1968 1969 from mesh_factory import rectangular 1970 from shallow_water import Domain, Transmissive_boundary 1971 from Numeric import zeros, Float 1972 from anuga.utilities.numerical_tools import mean 1973 1974 1975 #Create basic mesh 1976 points, vertices, boundary = rectangular(2, 2) 1977 1978 #Create shallow water domain 1979 domain = Domain(points, vertices, boundary) 1980 domain.default_order=2 1981 domain.reduction = mean 1982 1983 1984 #Set some field values 1985 domain.set_quantity('elevation', lambda x,y: x) 1986 domain.set_quantity('friction', 0.03) 1987 1988 1989 ###################### 1990 #Initial condition  with jumps 1991 1992 bed = domain.quantities['elevation'].vertex_values 1993 stage = zeros(bed.shape, Float) 1994 1995 h = 0.03 1996 for i in range(stage.shape[0]): 1997 if i % 2 == 0: 1998 stage[i,:] = bed[i,:] + h 1999 else: 2000 stage[i,:] = bed[i,:] 2001 2002 domain.set_quantity('stage', stage) 2003 2004 #Get stage 2005 stage = domain.quantities['stage'] 2006 A, V = stage.get_vertex_values(xy=False, smooth=False) 2007 Q = stage.vertex_values.flat 2008 2009 for k in range(8): 2010 assert allclose(A[k], Q[k]) 2011 2012 2013 for k in range(8): 2014 assert V[k, 0] == 3*k 2015 assert V[k, 1] == 3*k+1 2016 assert V[k, 2] == 3*k+2 2017 2018 2019 2020 X, Y, A1, V1 = stage.get_vertex_values(xy=True, smooth=False) 2021 2022 2023 assert allclose(A, A1) 2024 assert allclose(V, V1) 2025 2026 #Check XY 2027 assert allclose(X[1], 0.5) 2028 assert allclose(Y[1], 0.5) 2029 assert allclose(X[4], 0.0) 2030 assert allclose(Y[4], 0.0) 2031 assert allclose(X[12], 1.0) 2032 assert allclose(Y[12], 0.0) 2033 2034 2035 2036 def set_array_values_by_index(self): 2037 2038 from mesh_factory import rectangular 2039 from shallow_water import Domain 2040 from Numeric import zeros, Float 2041 2042 #Create basic mesh 2043 points, vertices, boundary = rectangular(1, 1) 2044 2045 #Create shallow water domain 2046 domain = Domain(points, vertices, boundary) 2047 #print "domain.number_of_elements ",domain.number_of_elements 2048 quantity = Quantity(domain,[[1,1,1],[2,2,2]]) 2049 value = [7] 2050 indices = [1] 2051 quantity.set_array_values_by_index(value, 2052 location = 'centroids', 2053 indices = indices) 2054 #print "quantity.centroid_values",quantity.centroid_values 2055 2056 assert allclose(quantity.centroid_values, [1,7]) 2057 2058 quantity.set_array_values([15,20,25], indices = indices) 2059 assert allclose(quantity.centroid_values, [1,20]) 2060 2061 quantity.set_array_values([15,20,25], indices = indices) 2062 assert allclose(quantity.centroid_values, [1,20]) 2063 2064 def test_setting_some_vertex_values(self): 2065 """ 2066 set values based on triangle lists. 2067 """ 2068 from mesh_factory import rectangular 2069 from shallow_water import Domain 2070 from Numeric import zeros, Float 2071 2072 #Create basic mesh 2073 points, vertices, boundary = rectangular(1, 3) 2074 #print "vertices",vertices 2075 #Create shallow water domain 2076 domain = Domain(points, vertices, boundary) 2077 #print "domain.number_of_elements ",domain.number_of_elements 2078 quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3], 2079 [4,4,4],[5,5,5],[6,6,6]]) 2080 2081 2082 # Check that constants work 2083 value = 7 2084 indices = [1] 2085 quantity.set_values(value, 2086 location = 'centroids', 2087 indices = indices) 2088 #print "quantity.centroid_values",quantity.centroid_values 2089 assert allclose(quantity.centroid_values, [1,7,3,4,5,6]) 2090 2091 value = [7] 2092 indices = [1] 2093 quantity.set_values(value, 2094 location = 'centroids', 2095 indices = indices) 2096 #print "quantity.centroid_values",quantity.centroid_values 2097 assert allclose(quantity.centroid_values, [1,7,3,4,5,6]) 2098 2099 value = [[15,20,25]] 2100 quantity.set_values(value, indices = indices) 2101 #print "1 quantity.vertex_values",quantity.vertex_values 2102 assert allclose(quantity.vertex_values[1], value[0]) 2103 2104 2105 #print "quantity",quantity.vertex_values 2106 values = [10,100,50] 2107 quantity.set_values(values, indices = [0,1,5], location = 'centroids') 2108 #print "2 quantity.vertex_values",quantity.vertex_values 2109 assert allclose(quantity.vertex_values[0], [10,10,10]) 2110 assert allclose(quantity.vertex_values[5], [50,50,50]) 2111 #quantity.interpolate() 2112 #print "quantity.centroid_values",quantity.centroid_values 2113 assert allclose(quantity.centroid_values, [10,100,3,4,5,50]) 2114 2115 2116 quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3], 2117 [4,4,4],[5,5,5],[6,6,6]]) 2118 values = [10,100,50] 2119 #this will be per unique vertex, indexing the vertices 2120 #print "quantity.vertex_values",quantity.vertex_values 2121 quantity.set_values(values, indices = [0,1,5]) 2122 #print "quantity.vertex_values",quantity.vertex_values 2123 assert allclose(quantity.vertex_values[0], [1,50,10]) 2124 assert allclose(quantity.vertex_values[5], [6,6,6]) 2125 assert allclose(quantity.vertex_values[1], [100,10,50]) 2126 2127 quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3], 2128 [4,4,4],[5,5,5],[6,6,6]]) 2129 values = [[31,30,29],[400,400,400],[1000,999,998]] 2130 quantity.set_values(values, indices = [3,3,5]) 2131 quantity.interpolate() 2132 assert allclose(quantity.centroid_values, [1,2,3,400,5,999]) 2133 2134 values = [[1,1,1],[2,2,2],[3,3,3], 2135 [4,4,4],[5,5,5],[6,6,6]] 2136 quantity.set_values(values) 2137 2138 # testing the standard set values by vertex 2139 # indexed by vertex_id in general_mesh.coordinates 2140 values = [0,1,2,3,4,5,6,7] 2141 2142 quantity.set_values(values) 2143 #print "1 quantity.vertex_values",quantity.vertex_values 2144 assert allclose(quantity.vertex_values,[[ 4., 5., 0.], 2145 [ 1., 0., 5.], 2146 [ 5., 6., 1.], 2147 [ 2., 1., 6.], 2148 [ 6., 7., 2.], 2149 [ 3., 2., 7.]]) 2150 2151 def test_setting_unique_vertex_values(self): 2152 """ 2153 set values based on unique_vertex lists. 2154 """ 2155 from mesh_factory import rectangular 2156 from shallow_water import Domain 2157 from Numeric import zeros, Float 2158 2159 #Create basic mesh 2160 points, vertices, boundary = rectangular(1, 3) 2161 #print "vertices",vertices 2162 #Create shallow water domain 2163 domain = Domain(points, vertices, boundary) 2164 #print "domain.number_of_elements ",domain.number_of_elements 2165 quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3], 2166 [4,4,4],[5,5,5]]) 2167 value = 7 2168 indices = [1,5] 2169 quantity.set_values(value, 2170 location = 'unique vertices', 2171 indices = indices) 2172 #print "quantity.centroid_values",quantity.centroid_values 2173 assert allclose(quantity.vertex_values[0], [0,7,0]) 2174 assert allclose(quantity.vertex_values[1], [7,1,7]) 2175 assert allclose(quantity.vertex_values[2], [7,2,7]) 2176 2177 2178 def test_get_values(self): 2179 """ 2180 get values based on triangle lists. 2181 """ 2182 from mesh_factory import rectangular 2183 from shallow_water import Domain 2184 from Numeric import zeros, Float 2185 2186 #Create basic mesh 2187 points, vertices, boundary = rectangular(1, 3) 2188 2189 #print "points",points 2190 #print "vertices",vertices 2191 #print "boundary",boundary 2192 2193 #Create shallow water domain 2194 domain = Domain(points, vertices, boundary) 2195 #print "domain.number_of_elements ",domain.number_of_elements 2196 quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3], 2197 [4,4,4],[5,5,5]]) 2198 2199 #print "quantity.get_values(location = 'unique vertices')", \ 2200 # quantity.get_values(location = 'unique vertices') 2201 2202 #print "quantity.get_values(location = 'unique vertices')", \ 2203 # quantity.get_values(indices=[0,1,2,3,4,5,6,7], \ 2204 # location = 'unique vertices') 2205 2206 answer = [0.5,2,4,5,0,1,3,4.5] 2207 assert allclose(answer, 2208 quantity.get_values(location = 'unique vertices')) 2209 2210 indices = [0,5,3] 2211 answer = [0.5,1,5] 2212 assert allclose(answer, 2213 quantity.get_values(indices=indices, \ 2214 location = 'unique vertices')) 2215 #print "quantity.centroid_values",quantity.centroid_values 2216 #print "quantity.get_values(location = 'centroids') ",\ 2217 # quantity.get_values(location = 'centroids') 2218 2219 2220 2221 2222 def test_get_values_2(self): 2223 """Different mesh (working with domain object)  also check centroids. 2224 """ 2225 2226 2227 a = [0.0, 0.0] 2228 b = [0.0, 2.0] 2229 c = [2.0,0.0] 2230 d = [0.0, 4.0] 2231 e = [2.0, 2.0] 2232 f = [4.0,0.0] 2233 2234 points = [a, b, c, d, e, f] 2235 #bac, bce, ecf, dbe 2236 vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2237 2238 domain = Domain(points, vertices) 2239 2240 quantity = Quantity(domain) 2241 quantity.set_values(lambda x, y: x+2*y) #2 4 4 6 2242 2243 assert allclose(quantity.get_values(location='centroids'), [2,4,4,6]) 2244 assert allclose(quantity.get_values(location='centroids', indices=[1,3]), [4,6]) 2245 2246 2247 assert allclose(quantity.get_values(location='vertices'), [[4,0,2], 2248 [4,2,6], 2249 [6,2,4], 2250 [8,4,6]]) 2251 2252 assert allclose(quantity.get_values(location='vertices', indices=[1,3]), [[4,2,6], 2253 [8,4,6]]) 2254 2255 2256 assert allclose(quantity.get_values(location='edges'), [[1,3,2], 2257 [4,5,3], 2258 [3,5,4], 2259 [5,7,6]]) 2260 assert allclose(quantity.get_values(location='edges', indices=[1,3]), 2261 [[4,5,3], 2262 [5,7,6]]) 2263 2264 # Check averaging over vertices 2265 #a: 0 2266 #b: (4+4+4)/3 2267 #c: (2+2+2)/3 2268 #d: 8 2269 #e: (6+6+6)/3 2270 #f: 4 2271 assert(quantity.get_values(location='unique vertices'), [0, 4, 2, 8, 6, 4]) 2272 2273 2274 2275 2276 2277 2278 def test_get_interpolated_values(self): 2279 2280 from mesh_factory import rectangular 2281 from shallow_water import Domain 2282 from Numeric import zeros, Float 2283 2284 #Create basic mesh 2285 points, vertices, boundary = rectangular(1, 3) 2286 domain = Domain(points, vertices, boundary) 2287 2288 #Constant values 2289 quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3], 2290 [4,4,4],[5,5,5]]) 2291 2292 2293 2294 # Get interpolated values at centroids 2295 interpolation_points = domain.get_centroid_coordinates() 2296 answer = quantity.get_values(location='centroids') 2297 2298 2299 #print quantity.get_values(points=interpolation_points) 2300 assert allclose(answer, quantity.get_values(interpolation_points=interpolation_points)) 2301 2302 2303 #Arbitrary values 2304 quantity = Quantity(domain,[[0,1,2],[3,1,7],[2,1,2],[3,3,7], 2305 [1,4,9],[2,5,0]]) 2306 2307 2308 # Get interpolated values at centroids 2309 interpolation_points = domain.get_centroid_coordinates() 2310 answer = quantity.get_values(location='centroids') 2311 #print answer 2312 #print quantity.get_values(interpolation_points=interpolation_points) 2313 assert allclose(answer, quantity.get_values(interpolation_points=interpolation_points, 2314 verbose=False)) 2315 2316 2317 #FIXME TODO 2318 #indices = [0,5,3] 2319 #answer = [0.5,1,5] 2320 #assert allclose(answer, 2321 # quantity.get_values(indices=indices, \ 2322 # location = 'unique vertices')) 2323 2324 2325 2326 2327 def test_get_interpolated_values_2(self): 2328 a = [0.0, 0.0] 2329 b = [0.0, 2.0] 2330 c = [2.0,0.0] 2331 d = [0.0, 4.0] 2332 e = [2.0, 2.0] 2333 f = [4.0,0.0] 2334 2335 points = [a, b, c, d, e, f] 2336 #bac, bce, ecf, dbe 2337 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2338 2339 domain = Domain(points, vertices) 2340 2341 quantity = Quantity(domain) 2342 quantity.set_values(lambda x, y: x+2*y) #2 4 4 6 2343 2344 #First pick one point 2345 x, y = 2.0/3, 8.0/3 2346 v = quantity.get_values(interpolation_points = [[x,y]]) 2347 assert allclose(v, 6) 2348 2349 # Then another to test that algorithm won't blindly 2350 # reuse interpolation matrix 2351 x, y = 4.0/3, 4.0/3 2352 v = quantity.get_values(interpolation_points = [[x,y]]) 2353 assert allclose(v, 4) 2354 2355 2356 2357 def test_get_interpolated_values_with_georef(self): 2358 2359 zone = 56 2360 xllcorner = 308500 2361 yllcorner = 6189000 2362 a = [0.0, 0.0] 2363 b = [0.0, 2.0] 2364 c = [2.0,0.0] 2365 d = [0.0, 4.0] 2366 e = [2.0, 2.0] 2367 f = [4.0,0.0] 2368 2369 points = [a, b, c, d, e, f] 2370 #bac, bce, ecf, dbe 2371 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2372 2373 domain = Domain(points, vertices, 2374 geo_reference=Geo_reference(zone,xllcorner,yllcorner)) 2375 2376 quantity = Quantity(domain) 2377 quantity.set_values(lambda x, y: x+2*y) #2 4 4 6 2378 2379 #First pick one point (and turn it into absolute coordinates) 2380 x, y = 2.0/3, 8.0/3 2381 v = quantity.get_values(interpolation_points = [[x+xllcorner,y+yllcorner]]) 2382 assert allclose(v, 6) 2383 2384 2385 # Then another to test that algorithm won't blindly 2386 # reuse interpolation matrix 2387 x, y = 4.0/3, 4.0/3 2388 v = quantity.get_values(interpolation_points = [[x+xllcorner,y+yllcorner]]) 2389 assert allclose(v, 4) 2390 2391 # Try two points 2392 pts = [[2.0/3 + xllcorner, 8.0/3 + yllcorner], 2393 [4.0/3 + xllcorner, 4.0/3 + yllcorner]] 2394 v = quantity.get_values(interpolation_points=pts) 2395 assert allclose(v, [6, 4]) 2396 2397 # Test it using the geospatial data format with absolute input points and default georef 2398 pts = Geospatial_data(data_points=pts) 2399 v = quantity.get_values(interpolation_points=pts) 2400 assert allclose(v, [6, 4]) 2401 2402 2403 # Test it using the geospatial data format with relative input points 2404 pts = Geospatial_data(data_points=[[2.0/3, 8.0/3], [4.0/3, 4.0/3]], 2405 geo_reference=Geo_reference(zone,xllcorner,yllcorner)) 2406 v = quantity.get_values(interpolation_points=pts) 2407 assert allclose(v, [6, 4]) 2408 2409 2410 2411 2412 def test_getting_some_vertex_values(self): 2413 """ 2414 get values based on triangle lists. 2415 """ 2416 from mesh_factory import rectangular 2417 from shallow_water import Domain 2418 from Numeric import zeros, Float 2419 2420 #Create basic mesh 2421 points, vertices, boundary = rectangular(1, 3) 2422 2423 #print "points",points 2424 #print "vertices",vertices 2425 #print "boundary",boundary 2426 2427 #Create shallow water domain 2428 domain = Domain(points, vertices, boundary) 2429 #print "domain.number_of_elements ",domain.number_of_elements 2430 quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3], 2431 [4,4,4],[5,5,5],[6,6,6]]) 2432 value = [7] 2433 indices = [1] 2434 quantity.set_values(value, 2435 location = 'centroids', 2436 indices = indices) 2437 #print "quantity.centroid_values",quantity.centroid_values 2438 #print "quantity.get_values(location = 'centroids') ",\ 2439 # quantity.get_values(location = 'centroids') 2440 assert allclose(quantity.centroid_values, 2441 quantity.get_values(location = 'centroids')) 2442 2443 2444 value = [[15,20,25]] 2445 quantity.set_values(value, indices = indices) 2446 #print "1 quantity.vertex_values",quantity.vertex_values 2447 assert allclose(quantity.vertex_values, quantity.get_values()) 2448 2449 assert allclose(quantity.edge_values, 2450 quantity.get_values(location = 'edges')) 2451 2452 # get a subset of elements 2453 subset = quantity.get_values(location='centroids', indices=[0,5]) 2454 answer = [quantity.centroid_values[0],quantity.centroid_values[5]] 2455 assert allclose(subset, answer) 2456 2457 2458 subset = quantity.get_values(location='edges', indices=[0,5]) 2459 answer = [quantity.edge_values[0],quantity.edge_values[5]] 2460 #print "subset",subset 2461 #print "answer",answer 2462 assert allclose(subset, answer) 2463 2464 subset = quantity.get_values( indices=[1,5]) 2465 answer = [quantity.vertex_values[1],quantity.vertex_values[5]] 2466 #print "subset",subset 2467 #print "answer",answer 2468 assert allclose(subset, answer) 2469 2470 def test_smooth_vertex_values(self): 2471 """ 2472 get values based on triangle lists. 2473 """ 2474 from mesh_factory import rectangular 2475 from shallow_water import Domain 2476 from Numeric import zeros, Float 2477 2478 #Create basic mesh 2479 points, vertices, boundary = rectangular(2, 2) 2480 2481 #print "points",points 2482 #print "vertices",vertices 2483 #print "boundary",boundary 2484 2485 #Create shallow water domain 2486 domain = Domain(points, vertices, boundary) 2487 #print "domain.number_of_elements ",domain.number_of_elements 2488 quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3], 2489 [4,4,4],[5,5,5],[6,6,6],[7,7,7]]) 2490 2491 #print "quantity.get_values(location = 'unique vertices')", \ 2492 # quantity.get_values(location = 'unique vertices') 2493 2494 #print "quantity.get_values(location = 'unique vertices')", \ 2495 # quantity.get_values(indices=[0,1,2,3,4,5,6,7], \ 2496 # location = 'unique vertices') 2497 2498 #print quantity.get_values(location = 'unique vertices') 2499 #print quantity.domain.number_of_triangles_per_node 2500 #print quantity.vertex_values 2501 2502 #answer = [0.5, 2, 3, 3, 3.5, 4, 4, 5, 6.5] 2503 #assert allclose(answer, 2504 # quantity.get_values(location = 'unique vertices')) 2505 2506 quantity.smooth_vertex_values() 2507 2508 #print quantity.vertex_values 2509 2510 2511 answer_vertex_values = [[3,3.5,0.5],[2,0.5,3.5],[3.5,4,2],[3,2,4], 2512 [4,5,3],[3.5,3,5],[5,6.5,3.5],[4,3.5,6.5]] 2513 2514 assert allclose(answer_vertex_values, 2515 quantity.vertex_values) 2516 #print "quantity.centroid_values",quantity.centroid_values 2517 #print "quantity.get_values(location = 'centroids') ",\ 2518 # quantity.get_values(location = 'centroids') 161 ## def test_get_maximum_2(self): 162 ## 163 ## a = [0.0, 0.0] 164 ## b = [0.0, 2.0] 165 ## c = [2.0,0.0] 166 ## d = [0.0, 4.0] 167 ## e = [2.0, 2.0] 168 ## f = [4.0,0.0] 169 ## 170 ## points = [a, b, c, d, e, f] 171 ## #bac, bce, ecf, dbe 172 ## vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 173 ## 174 ## domain = Domain(points, vertices) 175 ## 176 ## quantity = Quantity(domain) 177 ## quantity.set_values(lambda x, y: x+2*y) #2 4 4 6 178 ## 179 ## v = quantity.get_maximum_value() 180 ## assert v == 6 181 ## 182 ## v = quantity.get_minimum_value() 183 ## assert v == 2 184 ## 185 ## i = quantity.get_maximum_index() 186 ## assert i == 3 187 ## 188 ## i = quantity.get_minimum_index() 189 ## assert i == 0 190 ## 191 ## x,y = quantity.get_maximum_location() 192 ## xref, yref = 2.0/3, 8.0/3 193 ## assert x == xref 194 ## assert y == yref 195 ## 196 ## v = quantity.get_values(interpolation_points = [[x,y]]) 197 ## assert allclose(v, 6) 198 ## 199 ## x,y = quantity.get_minimum_location() 200 ## v = quantity.get_values(interpolation_points = [[x,y]]) 201 ## assert allclose(v, 2) 202 ## 203 ## #Multiple locations for maximum  204 ## #Test that the algorithm picks the first occurrence 205 ## v = quantity.get_maximum_value(indices=[0,1,2]) 206 ## assert allclose(v, 4) 207 ## 208 ## i = quantity.get_maximum_index(indices=[0,1,2]) 209 ## assert i == 1 210 ## 211 ## x,y = quantity.get_maximum_location(indices=[0,1,2]) 212 ## xref, yref = 4.0/3, 4.0/3 213 ## assert x == xref 214 ## assert y == yref 215 ## 216 ## v = quantity.get_values(interpolation_points = [[x,y]]) 217 ## assert allclose(v, 4) 218 ## 219 ## # More test of indices...... 220 ## v = quantity.get_maximum_value(indices=[2,3]) 221 ## assert allclose(v, 6) 222 ## 223 ## i = quantity.get_maximum_index(indices=[2,3]) 224 ## assert i == 3 225 ## 226 ## x,y = quantity.get_maximum_location(indices=[2,3]) 227 ## xref, yref = 2.0/3, 8.0/3 228 ## assert x == xref 229 ## assert y == yref 230 ## 231 ## v = quantity.get_values(interpolation_points = [[x,y]]) 232 ## assert allclose(v, 6) 233 ## 234 ## 235 ## 236 ## def test_boundary_allocation(self): 237 ## quantity = Quantity(self.mesh4, 238 ## [[1,2,3], [5,5,5], [0,0,9], [6, 3, 3]]) 239 ## 240 ## assert quantity.boundary_values.shape[0] == len(self.mesh4.boundary) 241 ## 242 ## 243 ## def test_set_values(self): 244 ## quantity = Quantity(self.mesh4) 245 ## 246 ## 247 ## quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [6, 3, 3]], 248 ## location = 'vertices') 249 ## assert allclose(quantity.vertex_values, 250 ## [[1,2,3], [5,5,5], [0,0,9], [6, 3, 3]]) 251 ## assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid 252 ## assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5], 253 ## [5., 5., 5.], 254 ## [4.5, 4.5, 0.], 255 ## [3.0, 1.5, 1.5]]) 256 ## 257 ## 258 ## # Test default 259 ## quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [6, 3, 3]]) 260 ## assert allclose(quantity.vertex_values, 261 ## [[1,2,3], [5,5,5], [0,0,9], [6, 3, 3]]) 262 ## assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid 263 ## assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5], 264 ## [5., 5., 5.], 265 ## [4.5, 4.5, 0.], 266 ## [3.0, 1.5, 1.5]]) 267 ## 268 ## # Test centroids 269 ## quantity.set_values([1,2,3,4], location = 'centroids') 270 ## assert allclose(quantity.centroid_values, [1., 2., 3., 4.]) #Centroid 271 ## 272 ## # Test exceptions 273 ## try: 274 ## quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [6, 3, 3]], 275 ## location = 'bas kamel tuba') 276 ## except: 277 ## pass 278 ## 279 ## 280 ## try: 281 ## quantity.set_values([[1,2,3], [0,0,9]]) 282 ## except AssertionError: 283 ## pass 284 ## except: 285 ## raise 'should have raised Assertionerror' 286 ## 287 ## 288 ## 289 ## def test_set_values_const(self): 290 ## quantity = Quantity(self.mesh4) 291 ## 292 ## quantity.set_values(1.0, location = 'vertices') 293 ## assert allclose(quantity.vertex_values, 294 ## [[1,1,1], [1,1,1], [1,1,1], [1, 1, 1]]) 295 ## 296 ## assert allclose(quantity.centroid_values, [1, 1, 1, 1]) #Centroid 297 ## assert allclose(quantity.edge_values, [[1, 1, 1], 298 ## [1, 1, 1], 299 ## [1, 1, 1], 300 ## [1, 1, 1]]) 301 ## 302 ## 303 ## quantity.set_values(2.0, location = 'centroids') 304 ## assert allclose(quantity.centroid_values, [2, 2, 2, 2]) 305 ## 306 ## 307 ## def test_set_values_func(self): 308 ## quantity = Quantity(self.mesh4) 309 ## 310 ## def f(x, y): 311 ## return x+y 312 ## 313 ## quantity.set_values(f, location = 'vertices') 314 ## #print "quantity.vertex_values",quantity.vertex_values 315 ## assert allclose(quantity.vertex_values, 316 ## [[2,0,2], [2,2,4], [4,2,4], [4,2,4]]) 317 ## assert allclose(quantity.centroid_values, 318 ## [4.0/3, 8.0/3, 10.0/3, 10.0/3]) 319 ## assert allclose(quantity.edge_values, 320 ## [[1,2,1], [3,3,2], [3,4,3], [3,4,3]]) 321 ## 322 ## 323 ## quantity.set_values(f, location = 'centroids') 324 ## assert allclose(quantity.centroid_values, 325 ## [4.0/3, 8.0/3, 10.0/3, 10.0/3]) 326 ## 327 ## 328 ## def test_integral(self): 329 ## quantity = Quantity(self.mesh4) 330 ## 331 ## #Try constants first 332 ## const = 5 333 ## quantity.set_values(const, location = 'vertices') 334 ## #print 'Q', quantity.get_integral() 335 ## 336 ## assert allclose(quantity.get_integral(), self.mesh4.get_area() * const) 337 ## 338 ## #Try with a linear function 339 ## def f(x, y): 340 ## return x+y 341 ## 342 ## quantity.set_values(f, location = 'vertices') 343 ## 344 ## 345 ## ref_integral = (4.0/3 + 8.0/3 + 10.0/3 + 10.0/3) * 2 346 ## 347 ## assert allclose (quantity.get_integral(), ref_integral) 348 ## 349 ## 350 ## 351 ## def test_set_vertex_values(self): 352 ## quantity = Quantity(self.mesh4) 353 ## quantity.set_vertex_values([0,1,2,3,4,5]) 354 ## 355 ## assert allclose(quantity.vertex_values, 356 ## [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 357 ## assert allclose(quantity.centroid_values, 358 ## [1., 7./3, 11./3, 8./3]) #Centroid 359 ## assert allclose(quantity.edge_values, [[1., 1.5, 0.5], 360 ## [3., 2.5, 1.5], 361 ## [3.5, 4.5, 3.], 362 ## [2.5, 3.5, 2]]) 363 ## 364 ## 365 ## def test_set_vertex_values_subset(self): 366 ## quantity = Quantity(self.mesh4) 367 ## quantity.set_vertex_values([0,1,2,3,4,5]) 368 ## quantity.set_vertex_values([0,20,30,50], indices = [0,2,3,5]) 369 ## 370 ## assert allclose(quantity.vertex_values, 371 ## [[1,0,20], [1,20,4], [4,20,50], [30,1,4]]) 372 ## 373 ## 374 ## def test_set_vertex_values_using_general_interface(self): 375 ## quantity = Quantity(self.mesh4) 376 ## 377 ## 378 ## quantity.set_values([0,1,2,3,4,5]) 379 ## 380 ## 381 ## assert allclose(quantity.vertex_values, 382 ## [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 383 ## 384 ## #Centroid 385 ## assert allclose(quantity.centroid_values, [1., 7./3, 11./3, 8./3]) 386 ## 387 ## assert allclose(quantity.edge_values, [[1., 1.5, 0.5], 388 ## [3., 2.5, 1.5], 389 ## [3.5, 4.5, 3.], 390 ## [2.5, 3.5, 2]]) 391 ## 392 ## 393 ## 394 ## def test_set_vertex_values_using_general_interface_with_subset(self): 395 ## """test_set_vertex_values_using_general_interface_with_subset(self): 396 ## 397 ## Test that indices and polygon works (for constants values) 398 ## """ 399 ## 400 ## quantity = Quantity(self.mesh4) 401 ## 402 ## 403 ## quantity.set_values([0,2,3,5], indices=[0,2,3,5]) 404 ## assert allclose(quantity.vertex_values, 405 ## [[0,0,2], [0,2,0], [0,2,5], [3,0,0]]) 406 ## 407 ## 408 ## # Constant 409 ## quantity.set_values(0.0) 410 ## quantity.set_values(3.14, indices=[0,2], location='vertices') 411 ## 412 ## # Indices refer to triangle numbers 413 ## assert allclose(quantity.vertex_values, 414 ## [[3.14,3.14,3.14], [0,0,0], 415 ## [3.14,3.14,3.14], [0,0,0]]) 416 ## 417 ## 418 ## 419 ## # Now try with polygon (pick points where y>2) 420 ## polygon = [[0,2.1], [4,2.1], [4,7], [0,7]] 421 ## quantity.set_values(0.0) 422 ## quantity.set_values(3.14, polygon=polygon) 423 ## 424 ## assert allclose(quantity.vertex_values, 425 ## [[0,0,0], [0,0,0], [0,0,0], 426 ## [3.14,3.14,3.14]]) 427 ## 428 ## 429 ## # Another polygon (pick triangle 1 and 2 (rightmost triangles) 430 ## # using centroids 431 ## polygon = [[2.1, 0.0], [3.5,0.1], [2,2.2], [0.2,2]] 432 ## quantity.set_values(0.0) 433 ## quantity.set_values(3.14, location='centroids', polygon=polygon) 434 ## assert allclose(quantity.vertex_values, 435 ## [[0,0,0], 436 ## [3.14,3.14,3.14], 437 ## [3.14,3.14,3.14], 438 ## [0,0,0]]) 439 ## 440 ## 441 ## # Same polygon now use vertices (default) 442 ## polygon = [[2.1, 0.0], [3.5,0.1], [2,2.2], [0.2,2]] 443 ## quantity.set_values(0.0) 444 ## #print 'Here 2' 445 ## quantity.set_values(3.14, polygon=polygon) 446 ## assert allclose(quantity.vertex_values, 447 ## [[0,0,0], 448 ## [3.14,3.14,3.14], 449 ## [3.14,3.14,3.14], 450 ## [0,0,0]]) 451 ## 452 ## 453 ## # Test input checking 454 ## try: 455 ## quantity.set_values(3.14, polygon=polygon, indices = [0,2]) 456 ## except: 457 ## pass 458 ## else: 459 ## msg = 'Should have caught this' 460 ## raise msg 461 ## 462 ## 463 ## 464 ## 465 ## 466 ## def test_set_vertex_values_using_general_interface_subset_and_geo(self): 467 ## """test_set_vertex_values_using_general_interface_with_subset(self): 468 ## Test that indices and polygon works using georeferencing 469 ## """ 470 ## 471 ## quantity = Quantity(self.mesh4) 472 ## G = Geo_reference(56, 10, 100) 473 ## quantity.domain.geo_reference = G 474 ## 475 ## #print quantity.domain.get_nodes(absolute=True) 476 ## 477 ## 478 ## # Constant 479 ## quantity.set_values(0.0) 480 ## quantity.set_values(3.14, indices=[0,2], location='vertices') 481 ## 482 ## # Indices refer to triangle numbers here  not vertices (why?) 483 ## assert allclose(quantity.vertex_values, 484 ## [[3.14,3.14,3.14], [0,0,0], 485 ## [3.14,3.14,3.14], [0,0,0]]) 486 ## 487 ## 488 ## 489 ## # Now try with polygon (pick points where y>2) 490 ## polygon = array([[0,2.1], [4,2.1], [4,7], [0,7]]) 491 ## polygon += [G.xllcorner, G.yllcorner] 492 ## 493 ## quantity.set_values(0.0) 494 ## quantity.set_values(3.14, polygon=polygon, location='centroids') 495 ## 496 ## assert allclose(quantity.vertex_values, 497 ## [[0,0,0], [0,0,0], [0,0,0], 498 ## [3.14,3.14,3.14]]) 499 ## 500 ## 501 ## # Another polygon (pick triangle 1 and 2 (rightmost triangles) 502 ## polygon = array([[2.1, 0.0], [3.5,0.1], [2,2.2], [0.2,2]]) 503 ## polygon += [G.xllcorner, G.yllcorner] 504 ## 505 ## quantity.set_values(0.0) 506 ## quantity.set_values(3.14, polygon=polygon) 507 ## 508 ## assert allclose(quantity.vertex_values, 509 ## [[0,0,0], 510 ## [3.14,3.14,3.14], 511 ## [3.14,3.14,3.14], 512 ## [0,0,0]]) 513 ## 514 ## 515 ## 516 ## def test_set_values_using_fit(self): 517 ## 518 ## 519 ## quantity = Quantity(self.mesh4) 520 ## 521 ## #Get (enough) datapoints 522 ## data_points = [[ 0.66666667, 0.66666667], 523 ## [ 1.33333333, 1.33333333], 524 ## [ 2.66666667, 0.66666667], 525 ## [ 0.66666667, 2.66666667], 526 ## [ 0.0, 1.0], 527 ## [ 0.0, 3.0], 528 ## [ 1.0, 0.0], 529 ## [ 1.0, 1.0], 530 ## [ 1.0, 2.0], 531 ## [ 1.0, 3.0], 532 ## [ 2.0, 1.0], 533 ## [ 3.0, 0.0], 534 ## [ 3.0, 1.0]] 535 ## 536 ## z = linear_function(data_points) 537 ## 538 ## #Use builtin fit_interpolate.fit 539 ## quantity.set_values( Geospatial_data(data_points, z), alpha = 0 ) 540 ## #quantity.set_values(points = data_points, values = z, alpha = 0) 541 ## 542 ## 543 ## answer = linear_function(quantity.domain.get_vertex_coordinates()) 544 ## #print quantity.vertex_values, answer 545 ## assert allclose(quantity.vertex_values.ravel(), answer) 546 ## 547 ## 548 ## #Now try by setting the same values directly 549 ## vertex_attributes = fit_to_mesh(data_points, 550 ## quantity.domain.get_nodes(), 551 ## quantity.domain.triangles, #FIXME 552 ## point_attributes=z, 553 ## alpha = 0, 554 ## verbose=False) 555 ## 556 ## #print vertex_attributes 557 ## quantity.set_values(vertex_attributes) 558 ## assert allclose(quantity.vertex_values.ravel(), answer) 559 ## 560 ## 561 ## 562 ## 563 ## 564 ## def test_test_set_values_using_fit_w_geo(self): 565 ## 566 ## 567 ## #Mesh 568 ## vertex_coordinates = [[0.76, 0.76], 569 ## [0.76, 5.76], 570 ## [5.76, 0.76]] 571 ## triangles = [[0,2,1]] 572 ## 573 ## mesh_georef = Geo_reference(56,0.76,0.76) 574 ## mesh1 = Domain(vertex_coordinates, triangles, 575 ## geo_reference = mesh_georef) 576 ## mesh1.check_integrity() 577 ## 578 ## #Quantity 579 ## quantity = Quantity(mesh1) 580 ## 581 ## #Data 582 ## data_points = [[ 201.0, 401.0], 583 ## [ 201.0, 403.0], 584 ## [ 203.0, 401.0]] 585 ## 586 ## z = [2, 4, 4] 587 ## 588 ## data_georef = Geo_reference(56,200,400) 589 ## 590 ## 591 ## #Reference 592 ## ref = fit_to_mesh(data_points, vertex_coordinates, triangles, 593 ## point_attributes=z, 594 ## data_origin = data_georef.get_origin(), 595 ## mesh_origin = mesh_georef.get_origin(), 596 ## alpha = 0) 597 ## 598 ## assert allclose( ref, [0,5,5] ) 599 ## 600 ## 601 ## #Test set_values 602 ## 603 ## quantity.set_values( Geospatial_data(data_points, z, data_georef), alpha = 0 ) 604 ## 605 ## #quantity.set_values(points = data_points, 606 ## # values = z, 607 ## # data_georef = data_georef, 608 ## # alpha = 0) 609 ## 610 ## 611 ## #quantity.set_values(points = data_points, 612 ## # values = z, 613 ## # data_georef = data_georef, 614 ## # alpha = 0) 615 ## assert allclose(quantity.vertex_values.ravel(), ref) 616 ## 617 ## 618 ## 619 ## #Test set_values using geospatial data object 620 ## quantity.vertex_values[:] = 0.0 621 ## 622 ## geo = Geospatial_data(data_points, z, data_georef) 623 ## 624 ## 625 ## quantity.set_values(geospatial_data = geo, alpha = 0) 626 ## assert allclose(quantity.vertex_values.ravel(), ref) 627 ## 628 ## 629 ## 630 ## def test_set_values_from_file1(self): 631 ## quantity = Quantity(self.mesh4) 632 ## 633 ## #Get (enough) datapoints 634 ## data_points = [[ 0.66666667, 0.66666667], 635 ## [ 1.33333333, 1.33333333], 636 ## [ 2.66666667, 0.66666667], 637 ## [ 0.66666667, 2.66666667], 638 ## [ 0.0, 1.0], 639 ## [ 0.0, 3.0], 640 ## [ 1.0, 0.0], 641 ## [ 1.0, 1.0], 642 ## [ 1.0, 2.0], 643 ## [ 1.0, 3.0], 644 ## [ 2.0, 1.0], 645 ## [ 3.0, 0.0], 646 ## [ 3.0, 1.0]] 647 ## 648 ## data_geo_spatial = Geospatial_data(data_points, 649 ## geo_reference = Geo_reference(56, 0, 0)) 650 ## data_points_absolute = data_geo_spatial.get_data_points(absolute=True) 651 ## attributes = linear_function(data_points_absolute) 652 ## att = 'spam_and_eggs' 653 ## 654 ## #Create .txt file 655 ## ptsfile = tempfile.mktemp(".txt") 656 ## file = open(ptsfile,"w") 657 ## file.write(" x,y," + att + " \n") 658 ## for data_point, attribute in map(None, data_points_absolute 659 ## ,attributes): 660 ## row = str(data_point[0]) + ',' + str(data_point[1]) \ 661 ## + ',' + str(attribute) 662 ## file.write(row + "\n") 663 ## file.close() 664 ## 665 ## 666 ## #Check that values can be set from file 667 ## quantity.set_values(filename = ptsfile, 668 ## attribute_name = att, alpha = 0) 669 ## answer = linear_function(quantity.domain.get_vertex_coordinates()) 670 ## 671 ## #print quantity.vertex_values.ravel() 672 ## #print answer 673 ## 674 ## 675 ## assert allclose(quantity.vertex_values.ravel(), answer) 676 ## 677 ## 678 ## #Check that values can be set from file using default attribute 679 ## quantity.set_values(filename = ptsfile, alpha = 0) 680 ## assert allclose(quantity.vertex_values.ravel(), answer) 681 ## 682 ## #Cleanup 683 ## import os 684 ## os.remove(ptsfile) 685 ## 686 ## 687 ## 688 ## def Xtest_set_values_from_file_using_polygon(self): 689 ## """test_set_values_from_file_using_polygon(self): 690 ## 691 ## Test that polygon restriction works for general points data 692 ## """ 693 ## 694 ## quantity = Quantity(self.mesh4) 695 ## 696 ## #Get (enough) datapoints 697 ## data_points = [[ 0.66666667, 0.66666667], 698 ## [ 1.33333333, 1.33333333], 699 ## [ 2.66666667, 0.66666667], 700 ## [ 0.66666667, 2.66666667], 701 ## [ 0.0, 1.0], 702 ## [ 0.0, 3.0], 703 ## [ 1.0, 0.0], 704 ## [ 1.0, 1.0], 705 ## [ 1.0, 2.0], 706 ## [ 1.0, 3.0], 707 ## [ 2.0, 1.0], 708 ## [ 3.0, 0.0], 709 ## [ 3.0, 1.0]] 710 ## 711 ## data_geo_spatial = Geospatial_data(data_points, 712 ## geo_reference = Geo_reference(56, 0, 0)) 713 ## data_points_absolute = data_geo_spatial.get_data_points(absolute=True) 714 ## attributes = linear_function(data_points_absolute) 715 ## att = 'spam_and_eggs' 716 ## 717 ## #Create .txt file 718 ## ptsfile = tempfile.mktemp(".txt") 719 ## file = open(ptsfile,"w") 720 ## file.write(" x,y," + att + " \n") 721 ## for data_point, attribute in map(None, data_points_absolute 722 ## ,attributes): 723 ## row = str(data_point[0]) + ',' + str(data_point[1]) \ 724 ## + ',' + str(attribute) 725 ## file.write(row + "\n") 726 ## file.close() 727 ## 728 ## # Create restricting polygon (containing node #4 (2,2) and 729 ## # centroid of triangle #1 (bce) 730 ## polygon = [[1.0, 1.0], [4.0, 1.0], 731 ## [4.0, 4.0], [1.0, 4.0]] 732 ## 733 ## #print self.mesh4.nodes 734 ## #print inside_polygon(self.mesh4.nodes, polygon) 735 ## assert allclose(inside_polygon(self.mesh4.nodes, polygon), 4) 736 ## 737 ## #print quantity.domain.get_vertex_coordinates() 738 ## #print quantity.domain.get_nodes() 739 ## 740 ## # Check that values can be set from file 741 ## quantity.set_values(filename=ptsfile, 742 ## polygon=polygon, 743 ## location='unique vertices', 744 ## alpha=0) 745 ## 746 ## # Get indices for vertex coordinates in polygon 747 ## indices = inside_polygon(quantity.domain.get_vertex_coordinates(), 748 ## polygon) 749 ## points = take(quantity.domain.get_vertex_coordinates(), indices) 750 ## 751 ## answer = linear_function(points) 752 ## 753 ## #print quantity.vertex_values.ravel() 754 ## #print answer 755 ## 756 ## # Check vertices in polygon have been set 757 ## assert allclose(take(quantity.vertex_values.ravel(), indices), 758 ## answer) 759 ## 760 ## # Check vertices outside polygon are zero 761 ## indices = outside_polygon(quantity.domain.get_vertex_coordinates(), 762 ## polygon) 763 ## assert allclose(take(quantity.vertex_values.ravel(), indices), 764 ## 0.0) 765 ## 766 ## #Cleanup 767 ## import os 768 ## os.remove(ptsfile) 769 ## 770 ## 771 ## 772 ## 773 ## def test_cache_test_set_values_from_file(self): 774 ## # FIXME (Ole): What is this about? 775 ## # I don't think it checks anything new 776 ## quantity = Quantity(self.mesh4) 777 ## 778 ## #Get (enough) datapoints 779 ## data_points = [[ 0.66666667, 0.66666667], 780 ## [ 1.33333333, 1.33333333], 781 ## [ 2.66666667, 0.66666667], 782 ## [ 0.66666667, 2.66666667], 783 ## [ 0.0, 1.0], 784 ## [ 0.0, 3.0], 785 ## [ 1.0, 0.0], 786 ## [ 1.0, 1.0], 787 ## [ 1.0, 2.0], 788 ## [ 1.0, 3.0], 789 ## [ 2.0, 1.0], 790 ## [ 3.0, 0.0], 791 ## [ 3.0, 1.0]] 792 ## 793 ## georef = Geo_reference(56, 0, 0) 794 ## data_geo_spatial = Geospatial_data(data_points, 795 ## geo_reference=georef) 796 ## 797 ## data_points_absolute = data_geo_spatial.get_data_points(absolute=True) 798 ## attributes = linear_function(data_points_absolute) 799 ## att = 'spam_and_eggs' 800 ## 801 ## # Create .txt file 802 ## ptsfile = tempfile.mktemp(".txt") 803 ## file = open(ptsfile,"w") 804 ## file.write(" x,y," + att + " \n") 805 ## for data_point, attribute in map(None, data_points_absolute 806 ## ,attributes): 807 ## row = str(data_point[0]) + ',' + str(data_point[1]) \ 808 ## + ',' + str(attribute) 809 ## file.write(row + "\n") 810 ## file.close() 811 ## 812 ## 813 ## # Check that values can be set from file 814 ## quantity.set_values(filename=ptsfile, 815 ## attribute_name=att, 816 ## alpha=0, 817 ## use_cache=True, 818 ## verbose=False) 819 ## answer = linear_function(quantity.domain.get_vertex_coordinates()) 820 ## assert allclose(quantity.vertex_values.ravel(), answer) 821 ## 822 ## 823 ## # Check that values can be set from file using default attribute 824 ## quantity.set_values(filename=ptsfile, 825 ## alpha=0) 826 ## assert allclose(quantity.vertex_values.ravel(), answer) 827 ## 828 ## # Check cache 829 ## quantity.set_values(filename=ptsfile, 830 ## attribute_name=att, 831 ## alpha=0, 832 ## use_cache=True, 833 ## verbose=False) 834 ## 835 ## 836 ## #Cleanup 837 ## import os 838 ## os.remove(ptsfile) 839 ## 840 ## def test_set_values_from_lat_long(self): 841 ## quantity = Quantity(self.mesh_onslow) 842 ## 843 ## #Get (enough) datapoints 844 ## data_points = [[21.5, 114.5],[21.4, 114.6],[21.45,114.65], 845 ## [21.35, 114.65],[21.45, 114.55],[21.45,114.6]] 846 ## 847 ## data_geo_spatial = Geospatial_data(data_points, 848 ## points_are_lats_longs=True) 849 ## points_UTM = data_geo_spatial.get_data_points(absolute=True) 850 ## attributes = linear_function(points_UTM) 851 ## att = 'elevation' 852 ## 853 ## #Create .txt file 854 ## txt_file = tempfile.mktemp(".txt") 855 ## file = open(txt_file,"w") 856 ## file.write(" lat,long," + att + " \n") 857 ## for data_point, attribute in map(None, data_points, attributes): 858 ## row = str(data_point[0]) + ',' + str(data_point[1]) \ 859 ## + ',' + str(attribute) 860 ## #print "row", row 861 ## file.write(row + "\n") 862 ## file.close() 863 ## 864 ## 865 ## #Check that values can be set from file 866 ## quantity.set_values(filename=txt_file, 867 ## attribute_name=att, 868 ## alpha=0) 869 ## answer = linear_function(quantity.domain.get_vertex_coordinates()) 870 ## 871 ## #print "quantity.vertex_values.ravel()", quantity.vertex_values.ravel() 872 ## #print "answer",answer 873 ## 874 ## assert allclose(quantity.vertex_values.ravel(), answer) 875 ## 876 ## 877 ## #Check that values can be set from file using default attribute 878 ## quantity.set_values(filename=txt_file, alpha=0) 879 ## assert allclose(quantity.vertex_values.ravel(), answer) 880 ## 881 ## #Cleanup 882 ## import os 883 ## os.remove(txt_file) 884 ## 885 ## def test_set_values_from_lat_long(self): 886 ## quantity = Quantity(self.mesh_onslow) 887 ## 888 ## #Get (enough) datapoints 889 ## data_points = [[21.5, 114.5],[21.4, 114.6],[21.45,114.65], 890 ## [21.35, 114.65],[21.45, 114.55],[21.45,114.6]] 891 ## 892 ## data_geo_spatial = Geospatial_data(data_points, 893 ## points_are_lats_longs=True) 894 ## points_UTM = data_geo_spatial.get_data_points(absolute=True) 895 ## attributes = linear_function(points_UTM) 896 ## att = 'elevation' 897 ## 898 ## #Create .txt file 899 ## txt_file = tempfile.mktemp(".txt") 900 ## file = open(txt_file,"w") 901 ## file.write(" lat,long," + att + " \n") 902 ## for data_point, attribute in map(None, data_points, attributes): 903 ## row = str(data_point[0]) + ',' + str(data_point[1]) \ 904 ## + ',' + str(attribute) 905 ## #print "row", row 906 ## file.write(row + "\n") 907 ## file.close() 908 ## 909 ## 910 ## #Check that values can be set from file 911 ## quantity.set_values(filename=txt_file, 912 ## attribute_name=att, alpha=0) 913 ## answer = linear_function(quantity.domain.get_vertex_coordinates()) 914 ## 915 ## #print "quantity.vertex_values.ravel()", quantity.vertex_values.ravel() 916 ## #print "answer",answer 917 ## 918 ## assert allclose(quantity.vertex_values.ravel(), answer) 919 ## 920 ## 921 ## #Check that values can be set from file using default attribute 922 ## quantity.set_values(filename=txt_file, alpha=0) 923 ## assert allclose(quantity.vertex_values.ravel(), answer) 924 ## 925 ## #Cleanup 926 ## import os 927 ## os.remove(txt_file) 928 ## 929 ## def test_set_values_from_UTM_pts(self): 930 ## quantity = Quantity(self.mesh_onslow) 931 ## 932 ## #Get (enough) datapoints 933 ## data_points = [[21.5, 114.5],[21.4, 114.6],[21.45,114.65], 934 ## [21.35, 114.65],[21.45, 114.55],[21.45,114.6]] 935 ## 936 ## data_geo_spatial = Geospatial_data(data_points, 937 ## points_are_lats_longs=True) 938 ## points_UTM = data_geo_spatial.get_data_points(absolute=True) 939 ## attributes = linear_function(points_UTM) 940 ## att = 'elevation' 941 ## 942 ## #Create .txt file 943 ## txt_file = tempfile.mktemp(".txt") 944 ## file = open(txt_file,"w") 945 ## file.write(" x,y," + att + " \n") 946 ## for data_point, attribute in map(None, points_UTM, attributes): 947 ## row = str(data_point[0]) + ',' + str(data_point[1]) \ 948 ## + ',' + str(attribute) 949 ## #print "row", row 950 ## file.write(row + "\n") 951 ## file.close() 952 ## 953 ## 954 ## pts_file = tempfile.mktemp(".pts") 955 ## convert = Geospatial_data(txt_file) 956 ## convert.export_points_file(pts_file) 957 ## 958 ## #Check that values can be set from file 959 ## quantity.set_values_from_file(pts_file, att, 0, 960 ## 'vertices', None) 961 ## answer = linear_function(quantity.domain.get_vertex_coordinates()) 962 ## #print "quantity.vertex_values.ravel()", quantity.vertex_values.ravel() 963 ## #print "answer",answer 964 ## assert allclose(quantity.vertex_values.ravel(), answer) 965 ## 966 ## #Check that values can be set from file 967 ## quantity.set_values(filename=pts_file, 968 ## attribute_name=att, alpha=0) 969 ## answer = linear_function(quantity.domain.get_vertex_coordinates()) 970 ## #print "quantity.vertex_values.ravel()", quantity.vertex_values.ravel() 971 ## #print "answer",answer 972 ## assert allclose(quantity.vertex_values.ravel(), answer) 973 ## 974 ## 975 ## #Check that values can be set from file using default attribute 976 ## quantity.set_values(filename=txt_file, alpha=0) 977 ## assert allclose(quantity.vertex_values.ravel(), answer) 978 ## 979 ## #Cleanup 980 ## import os 981 ## os.remove(txt_file) 982 ## os.remove(pts_file) 983 ## 984 ## def verbose_test_set_values_from_UTM_pts(self): 985 ## quantity = Quantity(self.mesh_onslow) 986 ## 987 ## #Get (enough) datapoints 988 ## data_points = [[21.5, 114.5],[21.4, 114.6],[21.45,114.65], 989 ## [21.35, 114.65],[21.45, 114.55],[21.45,114.6], 990 ## [21.5, 114.5],[21.4, 114.6],[21.45,114.65], 991 ## [21.35, 114.65],[21.45, 114.55],[21.45,114.6], 992 ## [21.5, 114.5],[21.4, 114.6],[21.45,114.65], 993 ## [21.35, 114.65],[21.45, 114.55],[21.45,114.6], 994 ## [21.5, 114.5],[21.4, 114.6],[21.45,114.65], 995 ## [21.35, 114.65],[21.45, 114.55],[21.45,114.6], 996 ## [21.5, 114.5],[21.4, 114.6],[21.45,114.65], 997 ## [21.35, 114.65],[21.45, 114.55],[21.45,114.6], 998 ## [21.5, 114.5],[21.4, 114.6],[21.45,114.65], 999 ## [21.35, 114.65],[21.45, 114.55],[21.45,114.6], 1000 ## [21.5, 114.5],[21.4, 114.6],[21.45,114.65], 1001 ## [21.35, 114.65],[21.45, 114.55],[21.45,114.6], 1002 ## [21.35, 114.65],[21.45, 114.55],[21.45,114.6], 1003 ## [21.5, 114.5],[21.4, 114.6],[21.45,114.65], 1004 ## [21.35, 114.65],[21.45, 114.55],[21.45,114.6], 1005 ## [21.5, 114.5],[21.4, 114.6],[21.45,114.65], 1006 ## [21.35, 114.65],[21.45, 114.55],[21.45,114.6], 1007 ## [21.5, 114.5],[21.4, 114.6],[21.45,114.65], 1008 ## [21.35, 114.65],[21.45, 114.55],[21.45,114.6], 1009 ## [21.5, 114.5],[21.4, 114.6],[21.45,114.65], 1010 ## [21.35, 114.65],[21.45, 114.55],[21.45,114.6], 1011 ## [21.5, 114.5],[21.4, 114.6],[21.45,114.65], 1012 ## [21.35, 114.65],[21.45, 114.55],[21.45,114.6], 1013 ## [21.5, 114.5],[21.4, 114.6],[21.45,114.65], 1014 ## [21.35, 114.65],[21.45, 114.55],[21.45,114.6], 1015 ## ] 1016 ## 1017 ## data_geo_spatial = Geospatial_data(data_points, 1018 ## points_are_lats_longs=True) 1019 ## points_UTM = data_geo_spatial.get_data_points(absolute=True) 1020 ## attributes = linear_function(points_UTM) 1021 ## att = 'elevation' 1022 ## 1023 ## #Create .txt file 1024 ## txt_file = tempfile.mktemp(".txt") 1025 ## file = open(txt_file,"w") 1026 ## file.write(" x,y," + att + " \n") 1027 ## for data_point, attribute in map(None, points_UTM, attributes): 1028 ## row = str(data_point[0]) + ',' + str(data_point[1]) \ 1029 ## + ',' + str(attribute) 1030 ## #print "row", row 1031 ## file.write(row + "\n") 1032 ## file.close() 1033 ## 1034 ## 1035 ## pts_file = tempfile.mktemp(".pts") 1036 ## convert = Geospatial_data(txt_file) 1037 ## convert.export_points_file(pts_file) 1038 ## 1039 ## #Check that values can be set from file 1040 ## quantity.set_values_from_file(pts_file, att, 0, 1041 ## 'vertices', None, verbose = True, 1042 ## max_read_lines=2) 1043 ## answer = linear_function(quantity.domain.get_vertex_coordinates()) 1044 ## #print "quantity.vertex_values.ravel()", quantity.vertex_values.ravel() 1045 ## #print "answer",answer 1046 ## assert allclose(quantity.vertex_values.ravel(), answer) 1047 ## 1048 ## #Check that values can be set from file 1049 ## quantity.set_values(filename=pts_file, 1050 ## attribute_name=att, alpha=0) 1051 ## answer = linear_function(quantity.domain.get_vertex_coordinates()) 1052 ## #print "quantity.vertex_values.ravel()", quantity.vertex_values.ravel() 1053 ## #print "answer",answer 1054 ## assert allclose(quantity.vertex_values.ravel(), answer) 1055 ## 1056 ## 1057 ## #Check that values can be set from file using default attribute 1058 ## quantity.set_values(filename=txt_file, alpha=0) 1059 ## assert allclose(quantity.vertex_values.ravel(), answer) 1060 ## 1061 ## #Cleanup 1062 ## import os 1063 ## os.remove(txt_file) 1064 ## os.remove(pts_file) 1065 ## 1066 ## def test_set_values_from_file_with_georef1(self): 1067 ## 1068 ## #Mesh in zone 56 (absolute coords) 1069 ## 1070 ## x0 = 314036.58727982 1071 ## y0 = 6224951.2960092 1072 ## 1073 ## a = [x0+0.0, y0+0.0] 1074 ## b = [x0+0.0, y0+2.0] 1075 ## c = [x0+2.0, y0+0.0] 1076 ## d = [x0+0.0, y0+4.0] 1077 ## e = [x0+2.0, y0+2.0] 1078 ## f = [x0+4.0, y0+0.0] 1079 ## 1080 ## points = [a, b, c, d, e, f] 1081 ## 1082 ## #bac, bce, ecf, dbe 1083 ## elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ] 1084 ## 1085 ## #absolute going in .. 1086 ## mesh4 = Domain(points, elements, 1087 ## geo_reference = Geo_reference(56, 0, 0)) 1088 ## mesh4.check_integrity() 1089 ## quantity = Quantity(mesh4) 1090 ## 1091 ## #Get (enough) datapoints (relative to georef) 1092 ## data_points_rel = [[ 0.66666667, 0.66666667], 1093 ## [ 1.33333333, 1.33333333], 1094 ## [ 2.66666667, 0.66666667], 1095 ## [ 0.66666667, 2.66666667], 1096 ## [ 0.0, 1.0], 1097 ## [ 0.0, 3.0], 1098 ## [ 1.0, 0.0], 1099 ## [ 1.0, 1.0], 1100 ## [ 1.0, 2.0], 1101 ## [ 1.0, 3.0], 1102 ## [ 2.0, 1.0], 1103 ## [ 3.0, 0.0], 1104 ## [ 3.0, 1.0]] 1105 ## 1106 ## data_geo_spatial = Geospatial_data(data_points_rel, 1107 ## geo_reference = Geo_reference(56, x0, y0)) 1108 ## data_points_absolute = data_geo_spatial.get_data_points(absolute=True) 1109 ## attributes = linear_function(data_points_absolute) 1110 ## att = 'spam_and_eggs' 1111 ## 1112 ## #Create .txt file 1113 ## ptsfile = tempfile.mktemp(".txt") 1114 ## file = open(ptsfile,"w") 1115 ## file.write(" x,y," + att + " \n") 1116 ## for data_point, attribute in map(None, data_points_absolute 1117 ## ,attributes): 1118 ## row = str(data_point[0]) + ',' + str(data_point[1]) \ 1119 ## + ',' + str(attribute) 1120 ## file.write(row + "\n") 1121 ## file.close() 1122 ## 1123 ## #file = open(ptsfile, 'r') 1124 ## #lines = file.readlines() 1125 ## #file.close() 1126 ## 1127 ## 1128 ## #Check that values can be set from file 1129 ## quantity.set_values(filename=ptsfile, 1130 ## attribute_name=att, alpha=0) 1131 ## answer = linear_function(quantity.domain.get_vertex_coordinates()) 1132 ## 1133 ## assert allclose(quantity.vertex_values.ravel(), answer) 1134 ## 1135 ## 1136 ## #Check that values can be set from file using default attribute 1137 ## quantity.set_values(filename=ptsfile, alpha=0) 1138 ## assert allclose(quantity.vertex_values.ravel(), answer) 1139 ## 1140 ## #Cleanup 1141 ## import os 1142 ## os.remove(ptsfile) 1143 ## 1144 ## 1145 ## def test_set_values_from_file_with_georef2(self): 1146 ## 1147 ## #Mesh in zone 56 (relative coords) 1148 ## 1149 ## x0 = 314036.58727982 1150 ## y0 = 6224951.2960092 1151 ## #x0 = 0.0 1152 ## #y0 = 0.0 1153 ## 1154 ## a = [0.0, 0.0] 1155 ## b = [0.0, 2.0] 1156 ## c = [2.0, 0.0] 1157 ## d = [0.0, 4.0] 1158 ## e = [2.0, 2.0] 1159 ## f = [4.0, 0.0] 1160 ## 1161 ## points = [a, b, c, d, e, f] 1162 ## 1163 ## #bac, bce, ecf, dbe 1164 ## elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ] 1165 ## 1166 ## mesh4 = Domain(points, elements, 1167 ## geo_reference = Geo_reference(56, x0, y0)) 1168 ## mesh4.check_integrity() 1169 ## quantity = Quantity(mesh4) 1170 ## 1171 ## #Get (enough) datapoints 1172 ## data_points = [[ x0+0.66666667, y0+0.66666667], 1173 ## [ x0+1.33333333, y0+1.33333333], 1174 ## [ x0+2.66666667, y0+0.66666667], 1175 ## [ x0+0.66666667, y0+2.66666667], 1176 ## [ x0+0.0, y0+1.0], 1177 ## [ x0+0.0, y0+3.0], 1178 ## [ x0+1.0, y0+0.0], 1179 ## [ x0+1.0, y0+1.0], 1180 ## [ x0+1.0, y0+2.0], 1181 ## [ x0+1.0, y0+3.0], 1182 ## [ x0+2.0, y0+1.0], 1183 ## [ x0+3.0, y0+0.0], 1184 ## [ x0+3.0, y0+1.0]] 1185 ## 1186 ## 1187 ## data_geo_spatial = Geospatial_data(data_points, 1188 ## geo_reference = Geo_reference(56, 0, 0)) 1189 ## data_points_absolute = data_geo_spatial.get_data_points(absolute=True) 1190 ## attributes = linear_function(data_points_absolute) 1191 ## att = 'spam_and_eggs' 1192 ## 1193 ## #Create .txt file 1194 ## ptsfile = tempfile.mktemp(".txt") 1195 ## file = open(ptsfile,"w") 1196 ## file.write(" x,y," + att + " \n") 1197 ## for data_point, attribute in map(None, data_points_absolute 1198 ## ,attributes): 1199 ## row = str(data_point[0]) + ',' + str(data_point[1]) \ 1200 ## + ',' + str(attribute) 1201 ## file.write(row + "\n") 1202 ## file.close() 1203 ## 1204 ## 1205 ## #Check that values can be set from file 1206 ## quantity.set_values(filename=ptsfile, 1207 ## attribute_name=att, alpha=0) 1208 ## answer = linear_function(quantity.domain. \ 1209 ## get_vertex_coordinates(absolute=True)) 1210 ## 1211 ## 1212 ## assert allclose(quantity.vertex_values.ravel(), answer) 1213 ## 1214 ## 1215 ## #Check that values can be set from file using default attribute 1216 ## quantity.set_values(filename=ptsfile, alpha=0) 1217 ## assert allclose(quantity.vertex_values.ravel(), answer) 1218 ## 1219 ## #Cleanup 1220 ## import os 1221 ## os.remove(ptsfile) 1222 ## 1223 ## 1224 ## 1225 ## 1226 ## def test_set_values_from_quantity(self): 1227 ## 1228 ## quantity1 = Quantity(self.mesh4) 1229 ## quantity1.set_vertex_values([0,1,2,3,4,5]) 1230 ## 1231 ## assert allclose(quantity1.vertex_values, 1232 ## [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 1233 ## 1234 ## 1235 ## quantity2 = Quantity(self.mesh4) 1236 ## quantity2.set_values(quantity=quantity1) 1237 ## assert allclose(quantity2.vertex_values, 1238 ## [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 1239 ## 1240 ## quantity2.set_values(quantity = 2*quantity1) 1241 ## assert allclose(quantity2.vertex_values, 1242 ## [[2,0,4], [2,4,8], [8,4,10], [6,2,8]]) 1243 ## 1244 ## quantity2.set_values(quantity = 2*quantity1 + 3) 1245 ## assert allclose(quantity2.vertex_values, 1246 ## [[5,3,7], [5,7,11], [11,7,13], [9,5,11]]) 1247 ## 1248 ## 1249 ## #Check detection of quantity as first orgument 1250 ## quantity2.set_values(2*quantity1 + 3) 1251 ## assert allclose(quantity2.vertex_values, 1252 ## [[5,3,7], [5,7,11], [11,7,13], [9,5,11]]) 1253 ## 1254 ## 1255 ## 1256 ## def Xtest_set_values_from_quantity_using_polygon(self): 1257 ## """test_set_values_from_quantity_using_polygon(self): 1258 ## 1259 ## Check that polygon can be used to restrict set_values when 1260 ## using another quantity as argument. 1261 ## """ 1262 ## 1263 ## # Create restricting polygon (containing node #4 (2,2) and 1264 ## # centroid of triangle #1 (bce) 1265 ## polygon = [[1.0, 1.0], [4.0, 1.0], 1266 ## [4.0, 4.0], [1.0, 4.0]] 1267 ## assert allclose(inside_polygon(self.mesh4.nodes, polygon), 4) 1268 ## 1269 ## quantity1 = Quantity(self.mesh4) 1270 ## quantity1.set_vertex_values([0,1,2,3,4,5]) 1271 ## 1272 ## assert allclose(quantity1.vertex_values, 1273 ## [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 1274 ## 1275 ## 1276 ## quantity2 = Quantity(self.mesh4) 1277 ## quantity2.set_values(quantity=quantity1, 1278 ## polygon=polygon) 1279 ## 1280 ## msg = 'Only node #4(e) at (2,2) should have values applied ' 1281 ## assert allclose(quantity2.vertex_values, 1282 ## [[0,0,0], [0,0,4], [4,0,0], [0,0,4]]), msg 1283 ## #bac, bce, ecf, dbe 1284 ## 1285 ## 1286 ## 1287 ## def test_overloading(self): 1288 ## 1289 ## quantity1 = Quantity(self.mesh4) 1290 ## quantity1.set_vertex_values([0,1,2,3,4,5]) 1291 ## 1292 ## assert allclose(quantity1.vertex_values, 1293 ## [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 1294 ## 1295 ## 1296 ## quantity2 = Quantity(self.mesh4) 1297 ## quantity2.set_values([[1,2,3], [5,5,5], [0,0,9], [6, 3, 3]], 1298 ## location = 'vertices') 1299 ## 1300 ## 1301 ## 1302 ## quantity3 = Quantity(self.mesh4) 1303 ## quantity3.set_values([[2,2,2], [7,8,9], [7,6,3], [3, 8, 8]], 1304 ## location = 'vertices') 1305 ## 1306 ## 1307 ## # Negation 1308 ## Q = quantity1 1309 ## assert allclose(Q.vertex_values, quantity1.vertex_values) 1310 ## assert allclose(Q.centroid_values, quantity1.centroid_values) 1311 ## assert allclose(Q.edge_values, quantity1.edge_values) 1312 ## 1313 ## # Addition 1314 ## Q = quantity1 + 7 1315 ## assert allclose(Q.vertex_values, quantity1.vertex_values + 7) 1316 ## assert allclose(Q.centroid_values, quantity1.centroid_values + 7) 1317 ## assert allclose(Q.edge_values, quantity1.edge_values + 7) 1318 ## 1319 ## Q = 7 + quantity1 1320 ## assert allclose(Q.vertex_values, quantity1.vertex_values + 7) 1321 ## assert allclose(Q.centroid_values, quantity1.centroid_values + 7) 1322 ## assert allclose(Q.edge_values, quantity1.edge_values + 7) 1323 ## 1324 ## Q = quantity1 + quantity2 1325 ## assert allclose(Q.vertex_values, 1326 ## quantity1.vertex_values + quantity2.vertex_values) 1327 ## assert allclose(Q.centroid_values, 1328 ## quantity1.centroid_values + quantity2.centroid_values) 1329 ## assert allclose(Q.edge_values, 1330 ## quantity1.edge_values + quantity2.edge_values) 1331 ## 1332 ## 1333 ## Q = quantity1 + quantity2  3 1334 ## assert allclose(Q.vertex_values, 1335 ## quantity1.vertex_values + quantity2.vertex_values  3) 1336 ## 1337 ## Q = quantity1  quantity2 1338 ## assert allclose(Q.vertex_values, 1339 ## quantity1.vertex_values  quantity2.vertex_values) 1340 ## 1341 ## #Scaling 1342 ## Q = quantity1*3 1343 ## assert allclose(Q.vertex_values, quantity1.vertex_values*3) 1344 ## assert allclose(Q.centroid_values, quantity1.centroid_values*3) 1345 ## assert allclose(Q.edge_values, quantity1.edge_values*3) 1346 ## Q = 3*quantity1 1347 ## assert allclose(Q.vertex_values, quantity1.vertex_values*3) 1348 ## 1349 ## #Multiplication 1350 ## Q = quantity1 * quantity2 1351 ## #print Q.vertex_values 1352 ## #print Q.centroid_values 1353 ## #print quantity1.centroid_values 1354 ## #print quantity2.centroid_values 1355 ## 1356 ## assert allclose(Q.vertex_values, 1357 ## quantity1.vertex_values * quantity2.vertex_values) 1358 ## 1359 ## #Linear combinations 1360 ## Q = 4*quantity1 + 2 1361 ## assert allclose(Q.vertex_values, 1362 ## 4*quantity1.vertex_values + 2) 1363 ## 1364 ## Q = quantity1*quantity2 + 2 1365 ## assert allclose(Q.vertex_values, 1366 ## quantity1.vertex_values * quantity2.vertex_values + 2) 1367 ## 1368 ## Q = quantity1*quantity2 + quantity3 1369 ## assert allclose(Q.vertex_values, 1370 ## quantity1.vertex_values * quantity2.vertex_values + 1371 ## quantity3.vertex_values) 1372 ## Q = quantity1*quantity2 + 3*quantity3 1373 ## assert allclose(Q.vertex_values, 1374 ## quantity1.vertex_values * quantity2.vertex_values + 1375 ## 3*quantity3.vertex_values) 1376 ## Q = quantity1*quantity2 + 3*quantity3 + 5.0 1377 ## assert allclose(Q.vertex_values, 1378 ## quantity1.vertex_values * quantity2.vertex_values + 1379 ## 3*quantity3.vertex_values + 5) 1380 ## 1381 ## Q = quantity1*quantity2  quantity3 1382 ## assert allclose(Q.vertex_values, 1383 ## quantity1.vertex_values * quantity2.vertex_values  1384 ## quantity3.vertex_values) 1385 ## Q = 1.5*quantity1*quantity2  3*quantity3 + 5.0 1386 ## assert allclose(Q.vertex_values, 1387 ## 1.5*quantity1.vertex_values * quantity2.vertex_values  1388 ## 3*quantity3.vertex_values + 5) 1389 ## 1390 ## #Try combining quantities and arrays and scalars 1391 ## Q = 1.5*quantity1*quantity2.vertex_values \ 1392 ## 3*quantity3.vertex_values + 5.0 1393 ## assert allclose(Q.vertex_values, 1394 ## 1.5*quantity1.vertex_values * quantity2.vertex_values  1395 ## 3*quantity3.vertex_values + 5) 1396 ## 1397 ## 1398 ## #Powers 1399 ## Q = quantity1**2 1400 ## assert allclose(Q.vertex_values, quantity1.vertex_values**2) 1401 ## 1402 ## Q = quantity1**2 +quantity2**2 1403 ## assert allclose(Q.vertex_values, 1404 ## quantity1.vertex_values**2 + \ 1405 ## quantity2.vertex_values**2) 1406 ## 1407 ## Q = (quantity1**2 +quantity2**2)**0.5 1408 ## assert allclose(Q.vertex_values, 1409 ## (quantity1.vertex_values**2 + \ 1410 ## quantity2.vertex_values**2)**0.5) 1411 ## 1412 ## 1413 ## 1414 ## 1415 ## 1416 ## 1417 ## 1418 ## def test_compute_gradient(self): 1419 ## quantity = Quantity(self.mesh4) 1420 ## 1421 ## #Set up for a gradient of (2,0) at mid triangle 1422 ## quantity.set_values([2.0, 4.0, 6.0, 2.0], 1423 ## location = 'centroids') 1424 ## 1425 ## 1426 ## #Gradients 1427 ## quantity.compute_gradients() 1428 ## 1429 ## a = quantity.x_gradient 1430 ## b = quantity.y_gradient 1431 ## #print self.mesh4.centroid_coordinates 1432 ## #print a, b 1433 ## 1434 ## #The central triangle (1) 1435 ## #(using standard gradient based on neigbours controid values) 1436 ## assert allclose(a[1], 2.0) 1437 ## assert allclose(b[1], 0.0) 1438 ## 1439 ## 1440 ## #Left triangle (0) using two point gradient 1441 ## #q0 = q1 + a*(x0x1) + b*(y0y1) <=> 1442 ## #2 = 4 + a*(2/3) + b*(2/3) 1443 ## assert allclose(a[0] + b[0], 3) 1444 ## #From orthogonality (a*(y0y1) + b*(x0x1) == 0) 1445 ## assert allclose(a[0]  b[0], 0) 1446 ## 1447 ## 1448 ## #Right triangle (2) using two point gradient 1449 ## #q2 = q1 + a*(x2x1) + b*(y2y1) <=> 1450 ## #6 = 4 + a*(4/3) + b*(2/3) 1451 ## assert allclose(2*a[2]  b[2], 3) 1452 ## #From orthogonality (a*(y1y2) + b*(x2x1) == 0) 1453 ## assert allclose(a[2] + 2*b[2], 0) 1454 ## 1455 ## 1456 ## #Top triangle (3) using two point gradient 1457 ## #q3 = q1 + a*(x3x1) + b*(y3y1) <=> 1458 ## #2 = 4 + a*(2/3) + b*(4/3) 1459 ## assert allclose(a[3]  2*b[3], 3) 1460 ## #From orthogonality (a*(y1y3) + b*(x3x1) == 0) 1461 ## assert allclose(2*a[3] + b[3], 0) 1462 ## 1463 ## 1464 ## 1465 ## #print a, b 1466 ## quantity.extrapolate_second_order() 1467 ## 1468 ## #Apply q(x,y) = qc + a*(xxc) + b*(yyc) 1469 ## assert allclose(quantity.vertex_values[0,:], [3., 0., 3.]) 1470 ## assert allclose(quantity.vertex_values[1,:], [4./3, 16./3, 16./3]) 1471 ## 1472 ## 1473 ## #a = 1.2, b=0.6 1474 ## #q(4,0) = 6 + a*(4  8/3) + b*(2/3) 1475 ## assert allclose(quantity.vertex_values[2,2], 8) 1476 ## 1477 ## def test_get_gradients(self): 1478 ## quantity = Quantity(self.mesh4) 1479 ## 1480 ## #Set up for a gradient of (2,0) at mid triangle 1481 ## quantity.set_values([2.0, 4.0, 6.0, 2.0], 1482 ## location = 'centroids') 1483 ## 1484 ## 1485 ## #Gradients 1486 ## quantity.compute_gradients() 1487 ## 1488 ## a, b = quantity.get_gradients() 1489 ## #print self.mesh4.centroid_coordinates 1490 ## #print a, b 1491 ## 1492 ## #The central triangle (1) 1493 ## #(using standard gradient based on neigbours controid values) 1494 ## assert allclose(a[1], 2.0) 1495 ## assert allclose(b[1], 0.0) 1496 ## 1497 ## 1498 ## #Left triangle (0) using two point gradient 1499 ## #q0 = q1 + a*(x0x1) + b*(y0y1) <=> 1500 ## #2 = 4 + a*(2/3) + b*(2/3) 1501 ## assert allclose(a[0] + b[0], 3) 1502 ## #From orthogonality (a*(y0y1) + b*(x0x1) == 0) 1503 ## assert allclose(a[0]  b[0], 0) 1504 ## 1505 ## 1506 ## #Right triangle (2) using two point gradient 1507 ## #q2 = q1 + a*(x2x1) + b*(y2y1) <=> 1508 ## #6 = 4 + a*(4/3) + b*(2/3) 1509 ## assert allclose(2*a[2]  b[2], 3) 1510 ## #From orthogonality (a*(y1y2) + b*(x2x1) == 0) 1511 ## assert allclose(a[2] + 2*b[2], 0) 1512 ## 1513 ## 1514 ## #Top triangle (3) using two point gradient 1515 ## #q3 = q1 + a*(x3x1) + b*(y3y1) <=> 1516 ## #2 = 4 + a*(2/3) + b*(4/3) 1517 ## assert allclose(a[3]  2*b[3], 3) 1518 ## #From orthogonality (a*(y1y3) + b*(x3x1) == 0) 1519 ## assert allclose(2*a[3] + b[3], 0) 1520 ## 1521 ## 1522 ## def test_second_order_extrapolation2(self): 1523 ## quantity = Quantity(self.mesh4) 1524 ## 1525 ## #Set up for a gradient of (3,1), f(x) = 3x+y 1526 ## quantity.set_values([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3], 1527 ## location = 'centroids') 1528 ## 1529 ## #Gradients 1530 ## quantity.compute_gradients() 1531 ## 1532 ## a = quantity.x_gradient 1533 ## b = quantity.y_gradient 1534 ## 1535 ## #print a, b 1536 ## 1537 ## assert allclose(a[1], 3.0) 1538 ## assert allclose(b[1], 1.0) 1539 ## 1540 ## #Work out the others 1541 ## 1542 ## quantity.extrapolate_second_order() 1543 ## 1544 ## #print quantity.vertex_values 1545 ## assert allclose(quantity.vertex_values[1,0], 2.0) 1546 ## assert allclose(quantity.vertex_values[1,1], 6.0) 1547 ## assert allclose(quantity.vertex_values[1,2], 8.0) 1548 ## 1549 ## 1550 ## 1551 ## def test_backup_saxpy_centroid_values(self): 1552 ## quantity = Quantity(self.mesh4) 1553 ## 1554 ## #Set up for a gradient of (3,1), f(x) = 3x+y 1555 ## c_values = array([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3]) 1556 ## d_values = array([1.0, 2.0, 3.0, 4.0]) 1557 ## quantity.set_values(c_values, location = 'centroids') 1558 ## 1559 ## #Backup 1560 ## quantity.backup_centroid_values() 1561 ## 1562 ## #print quantity.vertex_values 1563 ## assert allclose(quantity.centroid_values, quantity.centroid_backup_values) 1564 ## 1565 ## 1566 ## quantity.set_values(d_values, location = 'centroids') 1567 ## 1568 ## quantity.saxpy_centroid_values(2.0, 3.0) 1569 ## 1570 ## assert(quantity.centroid_values, 2.0*d_values + 3.0*c_values) 1571 ## 1572 ## 1573 ## 1574 ## def test_first_order_extrapolator(self): 1575 ## quantity = Quantity(self.mesh4) 1576 ## 1577 ## #Test centroids 1578 ## quantity.set_values([1.,2.,3.,4.], location = 'centroids') 1579 ## assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid 1580 ## 1581 ## #Extrapolate 1582 ## quantity.extrapolate_first_order() 1583 ## 1584 ## #Check that gradient is zero 1585 ## a,b = quantity.get_gradients() 1586 ## assert allclose(a, [0,0,0,0]) 1587 ## assert allclose(b, [0,0,0,0]) 1588 ## 1589 ## #Check vertices but not edge values 1590 ## assert allclose(quantity.vertex_values, 1591 ## [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]]) 1592 ## 1593 ## 1594 ## def test_second_order_extrapolator(self): 1595 ## quantity = Quantity(self.mesh4) 1596 ## 1597 ## #Set up for a gradient of (3,0) at mid triangle 1598 ## quantity.set_values([2.0, 4.0, 8.0, 2.0], 1599 ## location = 'centroids') 1600 ## 1601 ## 1602 ## 1603 ## quantity.extrapolate_second_order() 1604 ## quantity.limit() 1605 ## 1606 ## 1607 ## #Assert that central triangle is limited by neighbours 1608 ## assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0] 1609 ## assert quantity.vertex_values[1,0] >= quantity.vertex_values[3,1] 1610 ## 1611 ## assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1] 1612 ## assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2] 1613 ## 1614 ## assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0] 1615 ## assert quantity.vertex_values[1,2] >= quantity.vertex_values[3,1] 1616 ## 1617 ## 1618 ## #Assert that quantities are conserved 1619 ### from numpy.oldnumeric import sum 1620 ## from numpy import sum 1621 ## for k in range(quantity.centroid_values.shape[0]): 1622 ## assert allclose (quantity.centroid_values[k], 1623 ## sum(quantity.vertex_values[k,:])/3) 1624 ## 1625 ## 1626 ## 1627 ## 1628 ## 1629 ## def test_limit_vertices_by_all_neighbours(self): 1630 ## quantity = Quantity(self.mesh4) 1631 ## 1632 ## #Create a deliberate overshoot (e.g. from gradient computation) 1633 ## quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]]) 1634 ## 1635 ## 1636 ## #Limit 1637 ## quantity.limit_vertices_by_all_neighbours() 1638 ## 1639 ## #Assert that central triangle is limited by neighbours 1640 ## assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0] 1641 ## assert quantity.vertex_values[1,0] <= quantity.vertex_values[3,1] 1642 ## 1643 ## assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1] 1644 ## assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2] 1645 ## 1646 ## assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0] 1647 ## assert quantity.vertex_values[1,2] <= quantity.vertex_values[3,1] 1648 ## 1649 ## 1650 ## 1651 ## #Assert that quantities are conserved 1652 ### from numpy.oldnumeric import sum 1653 ## from numpy import sum 1654 ## for k in range(quantity.centroid_values.shape[0]): 1655 ## assert allclose (quantity.centroid_values[k], 1656 ## sum(quantity.vertex_values[k,:])/3) 1657 ## 1658 ## 1659 ## 1660 ## def test_limit_edges_by_all_neighbours(self): 1661 ## quantity = Quantity(self.mesh4) 1662 ## 1663 ## #Create a deliberate overshoot (e.g. from gradient computation) 1664 ## quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]]) 1665 ## 1666 ## 1667 ## #Limit 1668 ## quantity.limit_edges_by_all_neighbours() 1669 ## 1670 ## #Assert that central triangle is limited by neighbours 1671 ## assert quantity.edge_values[1,0] <= quantity.centroid_values[2] 1672 ## assert quantity.edge_values[1,0] >= quantity.centroid_values[0] 1673 ## 1674 ## assert quantity.edge_values[1,1] <= quantity.centroid_values[2] 1675 ## assert quantity.edge_values[1,1] >= quantity.centroid_values[0] 1676 ## 1677 ## assert quantity.edge_values[1,2] <= quantity.centroid_values[2] 1678 ## assert quantity.edge_values[1,2] >= quantity.centroid_values[0] 1679 ## 1680 ## 1681 ## 1682 ## #Assert that quantities are conserved 1683 ### from numpy.oldnumeric import sum 1684 ## from numpy import sum 1685 ## for k in range(quantity.centroid_values.shape[0]): 1686 ## assert allclose (quantity.centroid_values[k], 1687 ## sum(quantity.vertex_values[k,:])/3) 1688 ## 1689 ## 1690 ## def test_limit_edges_by_neighbour(self): 1691 ## quantity = Quantity(self.mesh4) 1692 ## 1693 ## #Create a deliberate overshoot (e.g. from gradient computation) 1694 ## quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]]) 1695 ## 1696 ## 1697 ## #Limit 1698 ## quantity.limit_edges_by_neighbour() 1699 ## 1700 ## #Assert that central triangle is limited by neighbours 1701 ## assert quantity.edge_values[1,0] <= quantity.centroid_values[3] 1702 ## assert quantity.edge_values[1,0] >= quantity.centroid_values[1] 1703 ## 1704 ## assert quantity.edge_values[1,1] <= quantity.centroid_values[2] 1705 ## assert quantity.edge_values[1,1] >= quantity.centroid_values[1] 1706 ## 1707 ## assert quantity.edge_values[1,2] <= quantity.centroid_values[1] 1708 ## assert quantity.edge_values[1,2] >= quantity.centroid_values[0] 1709 ## 1710 ## 1711 ## 1712 ## #Assert that quantities are conserved 1713 ### from numpy.oldnumeric import sum 1714 ## from numpy import sum 1715 ## for k in range(quantity.centroid_values.shape[0]): 1716 ## assert allclose (quantity.centroid_values[k], 1717 ## sum(quantity.vertex_values[k,:])/3) 1718 ## 1719 ## def test_limiter2(self): 1720 ## """Taken from test_shallow_water 1721 ## """ 1722 ## quantity = Quantity(self.mesh4) 1723 ## quantity.domain.beta_w = 0.9 1724 ## 1725 ## #Test centroids 1726 ## quantity.set_values([2.,4.,8.,2.], location = 'centroids') 1727 ## assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid 1728 ## 1729 ## 1730 ## #Extrapolate 1731 ## quantity.extrapolate_second_order() 1732 ## 1733 ## assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6]) 1734 ## 1735 ## #Limit 1736 ## quantity.limit() 1737 ## 1738 ## # limited value for beta_w = 0.9 1739 ## 1740 ## assert allclose(quantity.vertex_values[1,:], [2.2, 4.9, 4.9]) 1741 ## # limited values for beta_w = 0.5 1742 ## #assert allclose(quantity.vertex_values[1,:], [3.0, 4.5, 4.5]) 1743 ## 1744 ## 1745 ## #Assert that quantities are conserved 1746 ### from numpy.oldnumeric import sum 1747 ## from numpy import sum 1748 ## for k in range(quantity.centroid_values.shape[0]): 1749 ## assert allclose (quantity.centroid_values[k], 1750 ## sum(quantity.vertex_values[k,:])/3) 1751 ## 1752 ## 1753 ## 1754 ## 1755 ## 1756 ## def test_distribute_first_order(self): 1757 ## quantity = Quantity(self.mesh4) 1758 ## 1759 ## #Test centroids 1760 ## quantity.set_values([1.,2.,3.,4.], location = 'centroids') 1761 ## assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid 1762 ## 1763 ## 1764 ## #Extrapolate from centroid to vertices and edges 1765 ## quantity.extrapolate_first_order() 1766 ## 1767 ## #Interpolate 1768 ## #quantity.interpolate_from_vertices_to_edges() 1769 ## 1770 ## assert allclose(quantity.vertex_values, 1771 ## [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]]) 1772 ## assert allclose(quantity.edge_values, [[1,1,1], [2,2,2], 1773 ## [3,3,3], [4, 4, 4]]) 1774 ## 1775 ## 1776 ## def test_interpolate_from_vertices_to_edges(self): 1777 ## quantity = Quantity(self.mesh4) 1778 ## 1779 ## quantity.vertex_values = array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]],float) 1780 ## 1781 ## quantity.interpolate_from_vertices_to_edges() 1782 ## 1783 ## assert allclose(quantity.edge_values, [[1., 1.5, 0.5], 1784 ## [3., 2.5, 1.5], 1785 ## [3.5, 4.5, 3.], 1786 ## [2.5, 3.5, 2]]) 1787 ## 1788 ## 1789 ## def test_interpolate_from_edges_to_vertices(self): 1790 ## quantity = Quantity(self.mesh4) 1791 ## 1792 ## quantity.edge_values = array([[1., 1.5, 0.5], 1793 ## [3., 2.5, 1.5], 1794 ## [3.5, 4.5, 3.], 1795 ## [2.5, 3.5, 2]],float) 1796 ## 1797 ## quantity.interpolate_from_edges_to_vertices() 1798 ## 1799 ## assert allclose(quantity.vertex_values, 1800 ## [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 1801 ## 1802 ## 1803 ## 1804 ## def test_distribute_second_order(self): 1805 ## quantity = Quantity(self.mesh4) 1806 ## 1807 ## #Test centroids 1808 ## quantity.set_values([2.,4.,8.,2.], location = 'centroids') 1809 ## assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid 1810 ## 1811 ## 1812 ## #Extrapolate 1813 ## quantity.extrapolate_second_order() 1814 ## 1815 ## assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6]) 1816 ## 1817 ## 1818 ## def test_update_explicit(self): 1819 ## quantity = Quantity(self.mesh4) 1820 ## 1821 ## #Test centroids 1822 ## quantity.set_values([1.,2.,3.,4.], location = 'centroids') 1823 ## assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid 1824 ## 1825 ## #Set explicit_update 1826 ## quantity.explicit_update = array( [1.,1.,1.,1.] ) 1827 ## 1828 ## #Update with given timestep 1829 ## quantity.update(0.1) 1830 ## 1831 ## x = array([1, 2, 3, 4]) + array( [.1,.1,.1,.1] ) 1832 ## assert allclose( quantity.centroid_values, x) 1833 ## 1834 ## def test_update_semi_implicit(self): 1835 ## quantity = Quantity(self.mesh4) 1836 ## 1837 ## #Test centroids 1838 ## quantity.set_values([1.,2.,3.,4.], location = 'centroids') 1839 ## assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid 1840 ## 1841 ## #Set semi implicit update 1842 ## quantity.semi_implicit_update = array([1.,1.,1.,1.]) 1843 ## 1844 ## #Update with given timestep 1845 ## timestep = 0.1 1846 ## quantity.update(timestep) 1847 ## 1848 ## sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4]) 1849 ## denom = ones(4, float)timestep*sem 1850 ## 1851 ## x = array([1, 2, 3, 4])/denom 1852 ## assert allclose( quantity.centroid_values, x) 1853 ## 1854 ## 1855 ## def test_both_updates(self): 1856 ## quantity = Quantity(self.mesh4) 1857 ## 1858 ## #Test centroids 1859 ## quantity.set_values([1.,2.,3.,4.], location = 'centroids') 1860 ## assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid 1861 ## 1862 ## #Set explicit_update 1863 ## quantity.explicit_update = array( [4.,3.,2.,1.] ) 1864 ## 1865 ## #Set semi implicit update 1866 ## quantity.semi_implicit_update = array( [1.,1.,1.,1.] ) 1867 ## 1868 ## #Update with given timestep 1869 ## timestep = 0.1 1870 ## quantity.update(0.1) 1871 ## 1872 ## sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4]) 1873 ## denom = ones(4, float)timestep*sem 1874 ## 1875 ## x = array([1., 2., 3., 4.]) 1876 ## x /= denom 1877 ## x += timestep*array( [4.0, 3.0, 2.0, 1.0] ) 1878 ## 1879 ## assert allclose( quantity.centroid_values, x) 1880 ## 1881 ## 1882 ## 1883 ## 1884 ## #Test smoothing 1885 ## def test_smoothing(self): 1886 ## 1887 ## from mesh_factory import rectangular 1888 ## from shallow_water import Domain, Transmissive_boundary 1889 ### from numpy.oldnumeric import zeros, Float 1890 ## from numpy import zeros, float 1891 ## from anuga.utilities.numerical_tools import mean 1892 ## 1893 ## #Create basic mesh 1894 ## points, vertices, boundary = rectangular(2, 2) 1895 ## 1896 ## #Create shallow water domain 1897 ## domain = Domain(points, vertices, boundary) 1898 ## domain.default_order=2 1899 ## domain.reduction = mean 1900 ## 1901 ## 1902 ## #Set some field values 1903 ## domain.set_quantity('elevation', lambda x,y: x) 1904 ## domain.set_quantity('friction', 0.03) 1905 ## 1906 ## 1907 ## ###################### 1908 ## # Boundary conditions 1909 ## B = Transmissive_boundary(domain) 1910 ## domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B}) 1911 ## 1912 ## 1913 ## ###################### 1914 ## #Initial condition  with jumps 1915 ## 1916 ## bed = domain.quantities['elevation'].vertex_values 1917 ## stage = zeros(bed.shape, float) 1918 ## 1919 ## h = 0.03 1920 ## for i in range(stage.shape[0]): 1921 ## if i % 2 == 0: 1922 ## stage[i,:] = bed[i,:] + h 1923 ## else: 1924 ## stage[i,:] = bed[i,:] 1925 ## 1926 ## domain.set_quantity('stage', stage) 1927 ## 1928 ## stage = domain.quantities['stage'] 1929 ## 1930 ## #Get smoothed stage 1931 ## A, V = stage.get_vertex_values(xy=False, smooth=True) 1932 ## Q = stage.vertex_values 1933 ## 1934 ## 1935 ## assert A.shape[0] == 9 1936 ## assert V.shape[0] == 8 1937 ## assert V.shape[1] == 3 1938 ## 1939 ## #First four points 1940 ## assert allclose(A[0], (Q[0,2] + Q[1,1])/2) 1941 ## assert allclose(A[1], (Q[1,0] + Q[3,1] + Q[2,2])/3) 1942 ## assert allclose(A[2], Q[3,0]) 1943 ## assert allclose(A[3], (Q[0,0] + Q[5,1] + Q[4,2])/3) 1944 ## 1945 ## #Center point 1946 ## assert allclose(A[4], (Q[0,1] + Q[1,2] + Q[2,0] +\ 1947 ## Q[5,0] + Q[6,2] + Q[7,1])/6) 1948 ## 1949 ## 1950 ## #Check V 1951 ## assert allclose(V[0,:], [3,4,0]) 1952 ## assert allclose(V[1,:], [1,0,4]) 1953 ## assert allclose(V[2,:], [4,5,1]) 1954 ## assert allclose(V[3,:], [2,1,5]) 1955 ## assert allclose(V[4,:], [6,7,3]) 1956 ## assert allclose(V[5,:], [4,3,7]) 1957 ## assert allclose(V[6,:], [7,8,4]) 1958 ## assert allclose(V[7,:], [5,4,8]) 1959 ## 1960 ## #Get smoothed stage with XY 1961 ## X, Y, A1, V1 = stage.get_vertex_values(xy=True, smooth=True) 1962 ## 1963 ## assert allclose(A, A1) 1964 ## assert allclose(V, V1) 1965 ## 1966 ## #Check XY 1967 ## assert allclose(X[4], 0.5) 1968 ## assert allclose(Y[4], 0.5) 1969 ## 1970 ## assert allclose(X[7], 1.0) 1971 ## assert allclose(Y[7], 0.5) 1972 ## 1973 ## 1974 ## 1975 ## 1976 ## def test_vertex_values_no_smoothing(self): 1977 ## 1978 ## from mesh_factory import rectangular 1979 ## from shallow_water import Domain, Transmissive_boundary 1980 ### from numpy.oldnumeric import zeros, Float 1981 ## from numpy import zeros, float 1982 ## from anuga.utilities.numerical_tools import mean 1983 ## 1984 ## 1985 ## #Create basic mesh 1986 ## points, vertices, boundary = rectangular(2, 2) 1987 ## 1988 ## #Create shallow water domain 1989 ## domain = Domain(points, vertices, boundary) 1990 ## domain.default_order=2 1991 ## domain.reduction = mean 1992 ## 1993 ## 1994 ## #Set some field values 1995 ## domain.set_quantity('elevation', lambda x,y: x) 1996 ## domain.set_quantity('friction', 0.03) 1997 ## 1998 ## 1999 ## ###################### 2000 ## #Initial condition  with jumps 2001 ## 2002 ## bed = domain.quantities['elevation'].vertex_values 2003 ## stage = zeros(bed.shape, float) 2004 ## 2005 ## h = 0.03 2006 ## for i in range(stage.shape[0]): 2007 ## if i % 2 == 0: 2008 ## stage[i,:] = bed[i,:] + h 2009 ## else: 2010 ## stage[i,:] = bed[i,:] 2011 ## 2012 ## domain.set_quantity('stage', stage) 2013 ## 2014 ## #Get stage 2015 ## stage = domain.quantities['stage'] 2016 ## A, V = stage.get_vertex_values(xy=False, smooth=False) 2017 ## Q = stage.vertex_values.ravel() 2018 ## 2019 ## for k in range(8): 2020 ## assert allclose(A[k], Q[k]) 2021 ## 2022 ## 2023 ## for k in range(8): 2024 ## assert V[k, 0] == 3*k 2025 ## assert V[k, 1] == 3*k+1 2026 ## assert V[k, 2] == 3*k+2 2027 ## 2028 ## 2029 ## 2030 ## X, Y, A1, V1 = stage.get_vertex_values(xy=True, smooth=False) 2031 ## 2032 ## 2033 ## assert allclose(A, A1) 2034 ## assert allclose(V, V1) 2035 ## 2036 ## #Check XY 2037 ## assert allclose(X[1], 0.5) 2038 ## assert allclose(Y[1], 0.5) 2039 ## assert allclose(X[4], 0.0) 2040 ## assert allclose(Y[4], 0.0) 2041 ## assert allclose(X[12], 1.0) 2042 ## assert allclose(Y[12], 0.0) 2043 ## 2044 ## 2045 ## 2046 ## def set_array_values_by_index(self): 2047 ## 2048 ## from mesh_factory import rectangular 2049 ## from shallow_water import Domain 2050 ### from numpy.oldnumeric import zeros, Float 2051 ## from numpy import zeros, float 2052 ## 2053 ## #Create basic mesh 2054 ## points, vertices, boundary = rectangular(1, 1) 2055 ## 2056 ## #Create shallow water domain 2057 ## domain = Domain(points, vertices, boundary) 2058 ## #print "domain.number_of_elements ",domain.number_of_elements 2059 ## quantity = Quantity(domain,[[1,1,1],[2,2,2]]) 2060 ## value = [7] 2061 ## indices = [1] 2062 ## quantity.set_array_values_by_index(value, 2063 ## location = 'centroids', 2064 ## indices = indices) 2065 ## #print "quantity.centroid_values",quantity.centroid_values 2066 ## 2067 ## assert allclose(quantity.centroid_values, [1,7]) 2068 ## 2069 ## quantity.set_array_values([15,20,25], indices = indices) 2070 ## assert allclose(quantity.centroid_values, [1,20]) 2071 ## 2072 ## quantity.set_array_values([15,20,25], indices = indices) 2073 ## assert allclose(quantity.centroid_values, [1,20]) 2074 ## 2075 ## def test_setting_some_vertex_values(self): 2076 ## """ 2077 ## set values based on triangle lists. 2078 ## """ 2079 ## from mesh_factory import rectangular 2080 ## from shallow_water import Domain 2081 ### from numpy.oldnumeric import zeros, Float 2082 ## from numpy import zeros, float 2083 ## 2084 ## #Create basic mesh 2085 ## points, vertices, boundary = rectangular(1, 3) 2086 ## #print "vertices",vertices 2087 ## #Create shallow water domain 2088 ## domain = Domain(points, vertices, boundary) 2089 ## #print "domain.number_of_elements ",domain.number_of_elements 2090 ## quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3], 2091 ## [4,4,4],[5,5,5],[6,6,6]]) 2092 ## 2093 ## 2094 ## # Check that constants work 2095 ## value = 7 2096 ## indices = [1] 2097 ## quantity.set_values(value, 2098 ## location = 'centroids', 2099 ## indices = indices) 2100 ## #print "quantity.centroid_values",quantity.centroid_values 2101 ## assert allclose(quantity.centroid_values, [1,7,3,4,5,6]) 2102 ## 2103 ## value = [7] 2104 ## indices = [1] 2105 ## quantity.set_values(value, 2106 ## location = 'centroids', 2107 ## indices = indices) 2108 ## #print "quantity.centroid_values",quantity.centroid_values 2109 ## assert allclose(quantity.centroid_values, [1,7,3,4,5,6]) 2110 ## 2111 ## value = [[15,20,25]] 2112 ## quantity.set_values(value, indices = indices) 2113 ## #print "1 quantity.vertex_values",quantity.vertex_values 2114 ## assert allclose(quantity.vertex_values[1], value[0]) 2115 ## 2116 ## 2117 ## #print "quantity",quantity.vertex_values 2118 ## values = [10,100,50] 2119 ## quantity.set_values(values, indices = [0,1,5], location = 'centroids') 2120 ## #print "2 quantity.vertex_values",quantity.vertex_values 2121 ## assert allclose(quantity.vertex_values[0], [10,10,10]) 2122 ## assert allclose(quantity.vertex_values[5], [50,50,50]) 2123 ## #quantity.interpolate() 2124 ## #print "quantity.centroid_values",quantity.centroid_values 2125 ## assert allclose(quantity.centroid_values, [10,100,3,4,5,50]) 2126 ## 2127 ## 2128 ## quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3], 2129 ## [4,4,4],[5,5,5],[6,6,6]]) 2130 ## values = [10,100,50] 2131 ## #this will be per unique vertex, indexing the vertices 2132 ## #print "quantity.vertex_values",quantity.vertex_values 2133 ## quantity.set_values(values, indices = [0,1,5]) 2134 ## #print "quantity.vertex_values",quantity.vertex_values 2135 ## assert allclose(quantity.vertex_values[0], [1,50,10]) 2136 ## assert allclose(quantity.vertex_values[5], [6,6,6]) 2137 ## assert allclose(quantity.vertex_values[1], [100,10,50]) 2138 ## 2139 ## quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3], 2140 ## [4,4,4],[5,5,5],[6,6,6]]) 2141 ## values = [[31,30,29],[400,400,400],[1000,999,998]] 2142 ## quantity.set_values(values, indices = [3,3,5]) 2143 ## quantity.interpolate() 2144 ## assert allclose(quantity.centroid_values, [1,2,3,400,5,999]) 2145 ## 2146 ## values = [[1,1,1],[2,2,2],[3,3,3], 2147 ## [4,4,4],[5,5,5],[6,6,6]] 2148 ## quantity.set_values(values) 2149 ## 2150 ## # testing the standard set values by vertex 2151 ## # indexed by vertex_id in general_mesh.coordinates 2152 ## values = [0,1,2,3,4,5,6,7] 2153 ## 2154 ## quantity.set_values(values) 2155 ## #print "1 quantity.vertex_values",quantity.vertex_values 2156 ## assert allclose(quantity.vertex_values,[[ 4., 5., 0.], 2157 ## [ 1., 0., 5.], 2158 ## [ 5., 6., 1.], 2159 ## [ 2., 1., 6.], 2160 ## [ 6., 7., 2.], 2161 ## [ 3., 2., 7.]]) 2162 ## 2163 ## def test_setting_unique_vertex_values(self): 2164 ## """ 2165 ## set values based on unique_vertex lists. 2166 ## """ 2167 ## from mesh_factory import rectangular 2168 ## from shallow_water import Domain 2169 ### from numpy.oldnumeric import zeros, Float 2170 ## from numpy import zeros, float 2171 ## 2172 ## #Create basic mesh 2173 ## points, vertices, boundary = rectangular(1, 3) 2174 ## #print "vertices",vertices 2175 ## #Create shallow water domain 2176 ## domain = Domain(points, vertices, boundary) 2177 ## #print "domain.number_of_elements ",domain.number_of_elements 2178 ## quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3], 2179 ## [4,4,4],[5,5,5]]) 2180 ## value = 7 2181 ## indices = [1,5] 2182 ## quantity.set_values(value, 2183 ## location = 'unique vertices', 2184 ## indices = indices) 2185 ## #print "quantity.centroid_values",quantity.centroid_values 2186 ## assert allclose(quantity.vertex_values[0], [0,7,0]) 2187 ## assert allclose(quantity.vertex_values[1], [7,1,7]) 2188 ## assert allclose(quantity.vertex_values[2], [7,2,7]) 2189 ## 2190 ## 2191 ## def test_get_values(self): 2192 ## """ 2193 ## get values based on triangle lists. 2194 ## """ 2195 ## from mesh_factory import rectangular 2196 ## from shallow_water import Domain 2197 ### from numpy.oldnumeric import zeros, Float 2198 ## from numpy import zeros, float 2199 ## 2200 ## #Create basic mesh 2201 ## points, vertices, boundary = rectangular(1, 3) 2202 ## 2203 ## #print "points",points 2204 ## #print "vertices",vertices 2205 ## #print "boundary",boundary 2206 ## 2207 ## #Create shallow water domain 2208 ## domain = Domain(points, vertices, boundary) 2209 ## #print "domain.number_of_elements ",domain.number_of_elements 2210 ## quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3], 2211 ## [4,4,4],[5,5,5]]) 2212 ## 2213 ## #print "quantity.get_values(location = 'unique vertices')", \ 2214 ## # quantity.get_values(location = 'unique vertices') 2215 ## 2216 ## #print "quantity.get_values(location = 'unique vertices')", \ 2217 ## # quantity.get_values(indices=[0,1,2,3,4,5,6,7], \ 2218 ## # location = 'unique vertices') 2219 ## 2220 ## answer = [0.5,2,4,5,0,1,3,4.5] 2221 ## assert allclose(answer, 2222 ## quantity.get_values(location = 'unique vertices')) 2223 ## 2224 ## indices = [0,5,3] 2225 ## answer = [0.5,1,5] 2226 ## assert allclose(answer, 2227 ## quantity.get_values(indices=indices, \ 2228 ## location = 'unique vertices')) 2229 ## #print "quantity.centroid_values",quantity.centroid_values 2230 ## #print "quantity.get_values(location = 'centroids') ",\ 2231 ## # quantity.get_values(location = 'centroids') 2232 ## 2233 ## 2234 ## 2235 ## 2236 ## def test_get_values_2(self): 2237 ## """Different mesh (working with domain object)  also check centroids. 2238 ## """ 2239 ## 2240 ## 2241 ## a = [0.0, 0.0] 2242 ## b = [0.0, 2.0] 2243 ## c = [2.0,0.0] 2244 ## d = [0.0, 4.0] 2245 ## e = [2.0, 2.0] 2246 ## f = [4.0,0.0] 2247 ## 2248 ## points = [a, b, c, d, e, f] 2249 ## #bac, bce, ecf, dbe 2250 ## vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2251 ## 2252 ## domain = Domain(points, vertices) 2253 ## 2254 ## quantity = Quantity(domain) 2255 ## quantity.set_values(lambda x, y: x+2*y) #2 4 4 6 2256 ## 2257 ## assert allclose(quantity.get_values(location='centroids'), [2,4,4,6]) 2258 ## assert allclose(quantity.get_values(location='centroids', indices=[1,3]), [4,6]) 2259 ## 2260 ## 2261 ## assert allclose(quantity.get_values(location='vertices'), [[4,0,2], 2262 ## [4,2,6], 2263 ## [6,2,4], 2264 ## [8,4,6]]) 2265 ## 2266 ## assert allclose(quantity.get_values(location='vertices', indices=[1,3]), [[4,2,6], 2267 ## [8,4,6]]) 2268 ## 2269 ## 2270 ## assert allclose(quantity.get_values(location='edges'), [[1,3,2], 2271 ## [4,5,3], 2272 ## [3,5,4], 2273 ## [5,7,6]]) 2274 ## assert allclose(quantity.get_values(location='edges', indices=[1,3]), 2275 ## [[4,5,3], 2276 ## [5,7,6]]) 2277 ## 2278 ## # Check averaging over vertices 2279 ## #a: 0 2280 ## #b: (4+4+4)/3 2281 ## #c: (2+2+2)/3 2282 ## #d: 8 2283 ## #e: (6+6+6)/3 2284 ## #f: 4 2285 ## assert(quantity.get_values(location='unique vertices'), [0, 4, 2, 8, 6, 4]) 2286 ## 2287 ## 2288 ## 2289 ## 2290 ## 2291 ## 2292 ## def test_get_interpolated_values(self): 2293 ## 2294 ## from mesh_factory import rectangular 2295 ## from shallow_water import Domain 2296 ### from numpy.oldnumeric import zeros, Float 2297 ## from numpy import zeros, float 2298 ## 2299 ## #Create basic mesh 2300 ## points, vertices, boundary = rectangular(1, 3) 2301 ## domain = Domain(points, vertices, boundary) 2302 ## 2303 ## #Constant values 2304 ## quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3], 2305 ## [4,4,4],[5,5,5]]) 2306 ## 2307 ## 2308 ## 2309 ## # Get interpolated values at centroids 2310 ## interpolation_points = domain.get_centroid_coordinates() 2311 ## answer = quantity.get_values(location='centroids') 2312 ## 2313 ## 2314 ## #print quantity.get_values(points=interpolation_points) 2315 ## assert allclose(answer, quantity.get_values(interpolation_points=interpolation_points)) 2316 ## 2317 ## 2318 ## #Arbitrary values 2319 ## quantity = Quantity(domain,[[0,1,2],[3,1,7],[2,1,2],[3,3,7], 2320 ## [1,4,9],[2,5,0]]) 2321 ## 2322 ## 2323 ## # Get interpolated values at centroids 2324 ## interpolation_points = domain.get_centroid_coordinates() 2325 ## answer = quantity.get_values(location='centroids') 2326 ## #print answer 2327 ## #print quantity.get_values(interpolation_points=interpolation_points) 2328 ## assert allclose(answer, quantity.get_values(interpolation_points=interpolation_points, 2329 ## verbose=False)) 2330 ## 2331 ## 2332 ## #FIXME TODO 2333 ## #indices = [0,5,3] 2334 ## #answer = [0.5,1,5] 2335 ## #assert allclose(answer, 2336 ## # quantity.get_values(indices=indices, \ 2337 ## # location = 'unique vertices')) 2338 ## 2339 ## 2340 ## 2341 ## 2342 ## def test_get_interpolated_values_2(self): 2343 ## a = [0.0, 0.0] 2344 ## b = [0.0, 2.0] 2345 ## c = [2.0,0.0] 2346 ## d = [0.0, 4.0] 2347 ## e = [2.0, 2.0] 2348 ## f = [4.0,0.0] 2349 ## 2350 ## points = [a, b, c, d, e, f] 2351 ## #bac, bce, ecf, dbe 2352 ## vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2353 ## 2354 ## domain = Domain(points, vertices) 2355 ## 2356 ## quantity = Quantity(domain) 2357 ## quantity.set_values(lambda x, y: x+2*y) #2 4 4 6 2358 ## 2359 ## #First pick one point 2360 ## x, y = 2.0/3, 8.0/3 2361 ## v = quantity.get_values(interpolation_points = [[x,y]]) 2362 ## assert allclose(v, 6) 2363 ## 2364 ## # Then another to test that algorithm won't blindly 2365 ## # reuse interpolation matrix 2366 ## x, y = 4.0/3, 4.0/3 2367 ## v = quantity.get_values(interpolation_points = [[x,y]]) 2368 ## assert allclose(v, 4) 2369 ## 2370 ## 2371 ## 2372 ## def test_get_interpolated_values_with_georef(self): 2373 ## 2374 ## zone = 56 2375 ## xllcorner = 308500 2376 ## yllcorner = 6189000 2377 ## a = [0.0, 0.0] 2378 ## b = [0.0, 2.0] 2379 ## c = [2.0,0.0] 2380 ## d = [0.0, 4.0] 2381 ## e = [2.0, 2.0] 2382 ## f = [4.0,0.0] 2383 ## 2384 ## points = [a, b, c, d, e, f] 2385 ## #bac, bce, ecf, dbe 2386 ## vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2387 ## 2388 ## domain = Domain(points, vertices, 2389 ## geo_reference=Geo_reference(zone,xllcorner,yllcorner)) 2390 ## 2391 ## quantity = Quantity(domain) 2392 ## quantity.set_values(lambda x, y: x+2*y) #2 4 4 6 2393 ## 2394 ## #First pick one point (and turn it into absolute coordinates) 2395 ## x, y = 2.0/3, 8.0/3 2396 ## v = quantity.get_values(interpolation_points = [[x+xllcorner,y+yllcorner]]) 2397 ## assert allclose(v, 6) 2398 ## 2399 ## 2400 ## # Then another to test that algorithm won't blindly 2401 ## # reuse interpolation matrix 2402 ## x, y = 4.0/3, 4.0/3 2403 ## v = quantity.get_values(interpolation_points = [[x+xllcorner,y+yllcorner]]) 2404 ## assert allclose(v, 4) 2405 ## 2406 ## # Try two points 2407 ## pts = [[2.0/3 + xllcorner, 8.0/3 + yllcorner], 2408 ## [4.0/3 + xllcorner, 4.0/3 + yllcorner]] 2409 ## v = quantity.get_values(interpolation_points=pts) 2410 ## assert allclose(v, [6, 4]) 2411 ## 2412 ## # Test it using the geospatial data format with absolute input points and default georef 2413 ## pts = Geospatial_data(data_points=pts) 2414 ## v = quantity.get_values(interpolation_points=pts) 2415 ## assert allclose(v, [6, 4]) 2416 ## 2417 ## 2418 ## # Test it using the geospatial data format with relative input points 2419 ## pts = Geospatial_data(data_points=[[2.0/3, 8.0/3], [4.0/3, 4.0/3]], 2420 ## geo_reference=Geo_reference(zone,xllcorner,yllcorner)) 2421 ## v = quantity.get_values(interpolation_points=pts) 2422 ## assert allclose(v, [6, 4]) 2423 ## 2424 ## 2425 ## 2426 ## 2427 ## def test_getting_some_vertex_values(self): 2428 ## """ 2429 ## get values based on triangle lists. 2430 ## """ 2431 ## from mesh_factory import rectangular 2432 ## from shallow_water import Domain 2433 ### from numpy.oldnumeric import zeros, Float 2434 ## from numpy import zeros, float 2435 ## 2436 ## #Create basic mesh 2437 ## points, vertices, boundary = rectangular(1, 3) 2438 ## 2439 ## #print "points",points 2440 ## #print "vertices",vertices 2441 ## #print "boundary",boundary 2442 ## 2443 ## #Create shallow water domain 2444 ## domain = Domain(points, vertices, boundary) 2445 ## #print "domain.number_of_elements ",domain.number_of_elements 2446 ## quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3], 2447 ## [4,4,4],[5,5,5],[6,6,6]]) 2448 ## value = [7] 2449 ## indices = [1] 2450 ## quantity.set_values(value, 2451 ## location = 'centroids', 2452 ## indices = indices) 2453 ## #print "quantity.centroid_values",quantity.centroid_values 2454 ## #print "quantity.get_values(location = 'centroids') ",\ 2455 ## # quantity.get_values(location = 'centroids') 2456 ## assert allclose(quantity.centroid_values, 2457 ## quantity.get_values(location = 'centroids')) 2458 ## 2459 ## 2460 ## value = [[15,20,25]] 2461 ## quantity.set_values(value, indices = indices) 2462 ## #print "1 quantity.vertex_values",quantity.vertex_values 2463 ## assert allclose(quantity.vertex_values, quantity.get_values()) 2464 ## 2465 ## assert allclose(quantity.edge_values, 2466 ## quantity.get_values(location = 'edges')) 2467 ## 2468 ## # get a subset of elements 2469 ## subset = quantity.get_values(location='centroids', indices=[0,5]) 2470 ## answer = [quantity.centroid_values[0],quantity.centroid_values[5]] 2471 ## assert allclose(subset, answer) 2472 ## 2473 ## 2474 ## subset = quantity.get_values(location='edges', indices=[0,5]) 2475 ## answer = [quantity.edge_values[0],quantity.edge_values[5]] 2476 ## #print "subset",subset 2477 ## #print "answer",answer 2478 ## assert allclose(subset, answer) 2479 ## 2480 ## subset = quantity.get_values( indices=[1,5]) 2481 ## answer = [quantity.vertex_values[1],quantity.vertex_values[5]] 2482 ## #print "subset",subset 2483 ## #print "answer",answer 2484 ## assert allclose(subset, answer) 2485 ## 2486 ## def test_smooth_vertex_values(self): 2487 ## """ 2488 ## get values based on triangle lists. 2489 ## """ 2490 ## from mesh_factory import rectangular 2491 ## from shallow_water import Domain 2492 ### from numpy.oldnumeric import zeros, Float 2493 ## from numpy import zeros, float 2494 ## 2495 ## #Create basic mesh 2496 ## points, vertices, boundary = rectangular(2, 2) 2497 ## 2498 ## #print "points",points 2499 ## #print "vertices",vertices 2500 ## #print "boundary",boundary 2501 ## 2502 ## #Create shallow water domain 2503 ## domain = Domain(points, vertices, boundary) 2504 ## #print "domain.number_of_elements ",domain.number_of_elements 2505 ## quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3], 2506 ## [4,4,4],[5,5,5],[6,6,6],[7,7,7]]) 2507 ## 2508 ## #print "quantity.get_values(location = 'unique vertices')", \ 2509 ## # quantity.get_values(location = 'unique vertices') 2510 ## 2511 ## #print "quantity.get_values(location = 'unique vertices')", \ 2512 ## # quantity.get_values(indices=[0,1,2,3,4,5,6,7], \ 2513 ## # location = 'unique vertices') 2514 ## 2515 ## #print quantity.get_values(location = 'unique vertices') 2516 ## #print quantity.domain.number_of_triangles_per_node 2517 ## #print quantity.vertex_values 2518 ## 2519 ## #answer = [0.5, 2, 3, 3, 3.5, 4, 4, 5, 6.5] 2520 ## #assert allclose(answer, 2521 ## # quantity.get_values(location = 'unique vertices')) 2522 ## 2523 ## quantity.smooth_vertex_values() 2524 ## 2525 ## #print quantity.vertex_values 2526 ## 2527 ## 2528 ## answer_vertex_values = [[3,3.5,0.5],[2,0.5,3.5],[3.5,4,2],[3,2,4], 2529 ## [4,5,3],[3.5,3,5],[5,6.5,3.5],[4,3.5,6.5]] 2530 ## 2531 ## assert allclose(answer_vertex_values, 2532 ## quantity.vertex_values) 2533 ## #print "quantity.centroid_values",quantity.centroid_values 2534 ## #print "quantity.get_values(location = 'centroids') ",\ 2535 ## # quantity.get_values(location = 'centroids') 2519 2536 2520 2537
Note: See TracChangeset
for help on using the changeset viewer.