Changeset 2650
- Timestamp:
- Apr 2, 2006, 9:32:18 PM (17 years ago)
- Location:
- inundation/pyvolution
- Files:
-
- 3 edited
-
test_data_manager.py (modified) (6 diffs)
-
test_domain.py (modified) (9 diffs)
-
test_quantity.py (modified) (39 diffs)
Legend:
- Unmodified
- Added
- Removed
-
inundation/pyvolution/test_data_manager.py
r2648 r2650 2581 2581 2582 2582 def test_sww2domain1(self): 2583 2584 2585 2586 2587 2588 2589 2590 2591 2583 ################################################ 2584 #Create a test domain, and evolve and save it. 2585 ################################################ 2586 from mesh_factory import rectangular 2587 from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\ 2588 Constant_height, Time_boundary, Transmissive_boundary 2589 from Numeric import array 2590 2591 #Create basic mesh 2592 2592 2593 2593 yiel=0.01 2594 2595 2596 2597 2594 points, vertices, boundary = rectangular(10,10) 2595 2596 #Create shallow water domain 2597 domain = Domain(points, vertices, boundary) 2598 2598 domain.geo_reference = Geo_reference(56,11,11) 2599 domain.smooth = False 2600 domain.visualise = False 2601 domain.store = True 2602 domain.filename = 'bedslope' 2603 domain.default_order=2 2604 #Bed-slope and friction 2605 domain.set_quantity('elevation', lambda x,y: -x/3) 2606 domain.set_quantity('friction', 0.1) 2607 # Boundary conditions 2608 from math import sin, pi 2609 Br = Reflective_boundary(domain) 2610 Bt = Transmissive_boundary(domain) 2611 Bd = Dirichlet_boundary([0.2,0.,0.]) 2612 Bw = Time_boundary(domain=domain, 2613 f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0]) 2614 2615 #domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br}) 2616 domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd}) 2617 2618 domain.quantities_to_be_stored.extend(['xmomentum','ymomentum']) 2619 #Initial condition 2620 h = 0.05 2621 elevation = domain.quantities['elevation'].vertex_values 2622 domain.set_quantity('stage', elevation + h) 2623 2624 domain.check_integrity() 2625 #Evolution 2626 for t in domain.evolve(yieldstep = yiel, finaltime = 0.05): 2627 # domain.write_time() 2628 pass 2629 2630 2631 ########################################## 2632 #Import the example's file as a new domain 2633 ########################################## 2634 from data_manager import sww2domain 2635 from Numeric import allclose 2599 domain.smooth = False 2600 domain.visualise = False 2601 domain.store = True 2602 domain.filename = 'bedslope' 2603 domain.default_order=2 2604 #Bed-slope and friction 2605 domain.set_quantity('elevation', lambda x,y: -x/3) 2606 domain.set_quantity('friction', 0.1) 2607 # Boundary conditions 2608 from math import sin, pi 2609 Br = Reflective_boundary(domain) 2610 Bt = Transmissive_boundary(domain) 2611 Bd = Dirichlet_boundary([0.2,0.,0.]) 2612 Bw = Time_boundary(domain=domain,f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0]) 2613 2614 #domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br}) 2615 domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd}) 2616 2617 domain.quantities_to_be_stored.extend(['xmomentum','ymomentum']) 2618 #Initial condition 2619 h = 0.05 2620 elevation = domain.quantities['elevation'].vertex_values 2621 domain.set_quantity('stage', elevation + h) 2622 2623 domain.check_integrity() 2624 #Evolution 2625 for t in domain.evolve(yieldstep = yiel, finaltime = 0.05): 2626 #domain.write_time() 2627 pass 2628 2629 2630 ########################################## 2631 #Import the example's file as a new domain 2632 ########################################## 2633 from data_manager import sww2domain 2634 from Numeric import allclose 2636 2635 import os 2637 2636 2638 2639 2640 2637 filename = domain.datadir+os.sep+domain.filename+'.sww' 2638 domain2 = sww2domain(filename,None,fail_if_NaN=False,verbose = False) 2639 #points, vertices, boundary = rectangular(15,15) 2641 2640 #domain2.boundary = boundary 2642 2641 ################### 2643 2642 ##NOW TEST IT!!! 2644 2643 ################### 2645 2644 … … 2647 2646 os.remove(filename) 2648 2647 2649 2650 2651 2652 2653 2654 2655 2648 bits = ['vertex_coordinates'] 2649 for quantity in ['elevation']+domain.quantities_to_be_stored: 2650 bits.append('get_quantity("%s").get_integral()' %quantity) 2651 bits.append('get_quantity("%s").get_values()' %quantity) 2652 2653 for bit in bits: 2654 #print 'testing that domain.'+bit+' has been restored' 2656 2655 #print bit 2657 #print 'done'2658 2656 #print 'done' 2657 assert allclose(eval('domain.'+bit),eval('domain2.'+bit)) 2659 2658 2660 2659 ###################################### … … 2663 2662 visualise = False 2664 2663 #visualise = True 2665 2664 domain.visualise = visualise 2666 2665 domain.time = 0. 2667 2666 from time import sleep … … 2669 2668 final = .1 2670 2669 domain.set_quantity('friction', 0.1) 2671 2672 2673 2674 2675 2670 domain.store = False 2671 domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd}) 2672 2673 2674 for t in domain.evolve(yieldstep = yiel, finaltime = final): 2676 2675 if visualise: sleep(1.) 2677 2678 2676 #domain.write_time() 2677 pass 2679 2678 2680 2679 final = final - (domain2.starttime-domain.starttime) … … 2682 2681 final = final + (domain2.starttime-domain.starttime) 2683 2682 2684 2685 2686 2687 2688 2689 2690 2683 domain2.smooth = False 2684 domain2.visualise = visualise 2685 domain2.store = False 2686 domain2.default_order=2 2687 domain2.set_quantity('friction', 0.1) 2688 #Bed-slope and friction 2689 # Boundary conditions 2691 2690 Bd2=Dirichlet_boundary([0.2,0.,0.]) 2692 2691 domain2.boundary = domain.boundary 2693 2692 #print 'domain2.boundary' 2694 2693 #print domain2.boundary 2695 2696 2697 2698 2699 2700 2694 domain2.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd}) 2695 #domain2.set_boundary({'exterior': Bd}) 2696 2697 domain2.check_integrity() 2698 2699 for t in domain2.evolve(yieldstep = yiel, finaltime = final): 2701 2700 if visualise: sleep(1.) 2702 2703 2701 #domain2.write_time() 2702 pass 2704 2703 2705 2704 ################### 2706 2705 ##NOW TEST IT!!! 2707 2706 ################## 2708 2707 2709 bits = [ 'vertex_coordinates'] 2710 2711 for quantity in ['elevation','xmomentum','ymomentum']:#+domain.quantities_to_be_stored: 2712 bits.append('get_quantity("%s").get_integral()' %quantity) 2713 bits.append('get_quantity("%s").get_values()' %quantity) 2714 2715 for bit in bits: 2708 bits = [ 'vertex_coordinates'] 2709 2710 for quantity in ['elevation','xmomentum','ymomentum']:#+domain.quantities_to_be_stored: 2711 bits.append('get_quantity("%s").get_integral()' %quantity) 2712 bits.append('get_quantity("%s").get_values()' %quantity) 2713 2714 #print bits 2715 for bit in bits: 2716 2716 #print bit 2717 assert allclose(eval('domain.'+bit),eval('domain2.'+bit)) 2717 #print eval('domain.'+bit) 2718 #print eval('domain2.'+bit) 2719 assert allclose(eval('domain.'+bit),eval('domain2.'+bit)) 2718 2720 2719 2721 … … 4582 4584 4583 4585 4584 -
inundation/pyvolution/test_domain.py
r2491 r2650 139 139 def test_create_quantity_from_expression(self): 140 140 """Quantity created from other quantities using arbitrary expression 141 141 142 142 """ 143 143 … … 162 162 domain.set_quantity('elevation', -1) 163 163 164 164 165 165 domain.set_quantity('stage', [[1,2,3], [5,5,5], 166 166 [0,0,9], [-6, 3, 3]]) … … 170 170 171 171 domain.set_quantity('ymomentum', [[3,3,3], [4,2,1], 172 [2,4,-1], [1, 0, 1]]) 172 [2,4,-1], [1, 0, 1]]) 173 173 174 174 domain.check_integrity() … … 186 186 187 187 X = domain.quantities['xmomentum'].vertex_values 188 Y = domain.quantities['ymomentum'].vertex_values 188 Y = domain.quantities['ymomentum'].vertex_values 189 189 190 190 assert allclose(Q.vertex_values, (X**2 + Y**2)**0.5) 191 191 192 192 193 193 … … 196 196 def test_set_quantity_from_expression(self): 197 197 """Quantity set using arbitrary expression 198 198 199 199 """ 200 200 … … 219 219 domain.set_quantity('elevation', -1) 220 220 221 221 222 222 domain.set_quantity('stage', [[1,2,3], [5,5,5], 223 223 [0,0,9], [-6, 3, 3]]) … … 227 227 228 228 domain.set_quantity('ymomentum', [[3,3,3], [4,2,1], 229 [2,4,-1], [1, 0, 1]]) 229 [2,4,-1], [1, 0, 1]]) 230 230 231 231 … … 244 244 [1,1,10], [-5, 4, 4]]) 245 245 246 247 248 246 247 248 249 249 250 250 … … 426 426 denom = ones(4, Float)-domain.timestep*sem 427 427 428 x = array([1, 2, 3, 4]) + array( [.4,.3,.2,.1] ) 428 # x = array([1, 2, 3, 4]) + array( [.4,.3,.2,.1] ) 429 # x /= denom 430 431 x = array([1., 2., 3., 4.]) 429 432 x /= denom 433 x += domain.timestep*array( [4,3,2,1] ) 430 434 431 435 for name in domain.conserved_quantities: -
inundation/pyvolution/test_quantity.py
r2526 r2650 98 98 assert allclose(quantity.vertex_values[1,:],[3.+2./3, 6.+2./3, 4.+2./3]) 99 99 #FIXME: Work out the others 100 101 100 101 102 102 #print quantity.edge_values 103 103 assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5], … … 276 276 277 277 278 278 279 279 280 280 def test_set_values_using_least_squares(self): 281 281 282 282 from geospatial_data.geospatial_data import Geospatial_data 283 283 284 284 quantity = Quantity(self.mesh4) 285 285 … … 300 300 301 301 z = linear_function(data_points) 302 302 303 303 #Use built-in least squares fit 304 quantity.set_values( Geospatial_data(data_points, z), alpha = 0 ) 304 quantity.set_values( Geospatial_data(data_points, z), alpha = 0 ) 305 305 #quantity.set_values(points = data_points, values = z, alpha = 0) 306 306 307 307 308 308 answer = linear_function(quantity.domain.get_vertex_coordinates(obj = True)) 309 309 #print quantity.vertex_values, answer … … 312 312 313 313 #Now try by setting the same values directly 314 from least_squares import fit_to_mesh 314 from least_squares import fit_to_mesh 315 315 vertex_attributes = fit_to_mesh(quantity.domain.coordinates, 316 316 quantity.domain.triangles, … … 331 331 332 332 from domain import Domain 333 from geospatial_data.geospatial_data import Geospatial_data 333 from geospatial_data.geospatial_data import Geospatial_data 334 334 from coordinate_transforms.geo_reference import Geo_reference 335 335 from least_squares import fit_to_mesh … … 345 345 geo_reference = mesh_georef) 346 346 mesh1.check_integrity() 347 348 #Quantity 347 348 #Quantity 349 349 quantity = Quantity(mesh1) 350 350 351 #Data 351 #Data 352 352 data_points = [[ 201.0, 401.0], 353 353 [ 201.0, 403.0], … … 355 355 356 356 z = [2, 4, 4] 357 357 358 358 data_georef = Geo_reference(56,-200,-400) 359 359 360 360 361 361 #Reference … … 364 364 mesh_origin = mesh_georef.get_origin(), 365 365 alpha = 0) 366 367 assert allclose( ref, [0,5,5] ) 366 367 assert allclose( ref, [0,5,5] ) 368 368 369 369 370 370 #Test set_values 371 371 372 quantity.set_values( Geospatial_data(data_points, z, data_georef), alpha = 0 ) 372 quantity.set_values( Geospatial_data(data_points, z, data_georef), alpha = 0 ) 373 373 374 374 #quantity.set_values(points = data_points, … … 377 377 # alpha = 0) 378 378 379 379 380 380 #quantity.set_values(points = data_points, 381 381 # values = z, … … 388 388 #Test set_values using geospatial data object 389 389 quantity.vertex_values[:] = 0.0 390 390 391 391 from geospatial_data.geospatial_data import Geospatial_data 392 392 geo = Geospatial_data(data_points, z, data_georef) 393 393 394 394 395 395 quantity.set_values(geospatial_data = geo, alpha = 0) 396 assert allclose(quantity.vertex_values.flat, ref) 396 assert allclose(quantity.vertex_values.flat, ref) 397 397 398 398 … … 417 417 418 418 z = linear_function(data_points) 419 419 420 420 421 421 #Create pts file … … 425 425 points_dict = {'pointlist': data_points, 426 426 'attributelist': {att: z}} 427 428 export_points_file(ptsfile, points_dict) 427 428 export_points_file(ptsfile, points_dict) 429 429 430 430 #Check that values can be set from file … … 436 436 #print answer 437 437 438 438 439 439 assert allclose(quantity.vertex_values.flat, answer) 440 440 … … 442 442 #Check that values can be set from file using default attribute 443 443 quantity.set_values(filename = ptsfile, alpha = 0) 444 assert allclose(quantity.vertex_values.flat, answer) 444 assert allclose(quantity.vertex_values.flat, answer) 445 445 446 446 #Cleanup 447 447 import os 448 448 os.remove(ptsfile) 449 449 450 450 451 451 def test_set_values_from_file_with_georef1(self): 452 452 453 453 #Mesh in zone 56 (absolute coords) 454 454 from domain import Domain 455 455 from coordinate_transforms.geo_reference import Geo_reference 456 456 457 457 x0 = 314036.58727982 458 458 y0 = 6224951.2960092 … … 491 491 492 492 z = linear_function(data_points) 493 493 494 494 495 495 #Create pts file 496 from load_mesh.loadASCII import export_points_file 496 from load_mesh.loadASCII import export_points_file 497 497 498 498 ptsfile = 'testptsfile.pts' … … 504 504 xllcorner = x0, 505 505 yllcorner = y0)} 506 506 507 507 export_points_file(ptsfile, points_dict) 508 508 509 509 510 510 #Check that values can be set from file … … 512 512 attribute_name = att, alpha = 0) 513 513 answer = linear_function(quantity.domain.get_vertex_coordinates(obj = True) - [x0, y0]) 514 514 515 515 assert allclose(quantity.vertex_values.flat, answer) 516 516 … … 518 518 #Check that values can be set from file using default attribute 519 519 quantity.set_values(filename = ptsfile, alpha = 0) 520 assert allclose(quantity.vertex_values.flat, answer) 520 assert allclose(quantity.vertex_values.flat, answer) 521 521 522 522 #Cleanup … … 526 526 527 527 def test_set_values_from_file_with_georef2(self): 528 528 529 529 #Mesh in zone 56 (relative coords) 530 530 from domain import Domain 531 531 from coordinate_transforms.geo_reference import Geo_reference 532 532 533 533 x0 = 314036.58727982 534 534 y0 = 6224951.2960092 … … 567 567 568 568 z = linear_function(data_points) 569 569 570 570 571 571 #Create pts file 572 from load_mesh.loadASCII import export_points_file 572 from load_mesh.loadASCII import export_points_file 573 573 574 574 ptsfile = 'testptsfile.pts' … … 580 580 xllcorner = 0, 581 581 yllcorner = 0)} 582 582 583 583 export_points_file(ptsfile, points_dict) 584 584 585 585 586 586 #Check that values can be set from file … … 589 589 answer = linear_function(quantity.domain.get_vertex_coordinates(obj = True) + [x0, y0]) 590 590 591 591 592 592 assert allclose(quantity.vertex_values.flat, answer) 593 593 … … 595 595 #Check that values can be set from file using default attribute 596 596 quantity.set_values(filename = ptsfile, alpha = 0) 597 assert allclose(quantity.vertex_values.flat, answer) 597 assert allclose(quantity.vertex_values.flat, answer) 598 598 599 599 #Cleanup 600 600 import os 601 601 os.remove(ptsfile) 602 602 603 603 604 604 … … 612 612 [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 613 613 614 614 615 615 quantity2 = Quantity(self.mesh4) 616 616 quantity2.set_values(quantity = quantity1) … … 621 621 assert allclose(quantity2.vertex_values, 622 622 [[2,0,4], [2,4,8], [8,4,10], [6,2,8]]) 623 623 624 624 quantity2.set_values(quantity = 2*quantity1 + 3) 625 625 assert allclose(quantity2.vertex_values, … … 630 630 quantity2.set_values(2*quantity1 + 3) 631 631 assert allclose(quantity2.vertex_values, 632 [[5,3,7], [5,7,11], [11,7,13], [9,5,11]]) 632 [[5,3,7], [5,7,11], [11,7,13], [9,5,11]]) 633 633 634 634 … … 644 644 [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 645 645 646 646 647 647 quantity2 = Quantity(self.mesh4) 648 648 quantity2.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]], … … 653 653 quantity3 = Quantity(self.mesh4) 654 654 quantity3.set_values([[2,2,2], [7,8,9], [7,6,3], [3, 8, -8]], 655 location = 'vertices') 655 location = 'vertices') 656 656 657 657 … … 660 660 assert allclose(Q.vertex_values, -quantity1.vertex_values) 661 661 assert allclose(Q.centroid_values, -quantity1.centroid_values) 662 assert allclose(Q.edge_values, -quantity1.edge_values) 663 662 assert allclose(Q.edge_values, -quantity1.edge_values) 663 664 664 #Addition 665 665 Q = quantity1 + 7 666 666 assert allclose(Q.vertex_values, quantity1.vertex_values + 7) 667 667 assert allclose(Q.centroid_values, quantity1.centroid_values + 7) 668 assert allclose(Q.edge_values, quantity1.edge_values + 7) 669 668 assert allclose(Q.edge_values, quantity1.edge_values + 7) 669 670 670 Q = 7 + quantity1 671 assert allclose(Q.vertex_values, quantity1.vertex_values + 7) 671 assert allclose(Q.vertex_values, quantity1.vertex_values + 7) 672 672 assert allclose(Q.centroid_values, quantity1.centroid_values + 7) 673 assert allclose(Q.edge_values, quantity1.edge_values + 7) 674 673 assert allclose(Q.edge_values, quantity1.edge_values + 7) 674 675 675 Q = quantity1 + quantity2 676 676 assert allclose(Q.vertex_values, … … 679 679 quantity1.centroid_values + quantity2.centroid_values) 680 680 assert allclose(Q.edge_values, 681 quantity1.edge_values + quantity2.edge_values) 682 681 quantity1.edge_values + quantity2.edge_values) 682 683 683 684 684 Q = quantity1 + quantity2 - 3 … … 688 688 Q = quantity1 - quantity2 689 689 assert allclose(Q.vertex_values, 690 quantity1.vertex_values - quantity2.vertex_values) 690 quantity1.vertex_values - quantity2.vertex_values) 691 691 692 692 #Scaling … … 694 694 assert allclose(Q.vertex_values, quantity1.vertex_values*3) 695 695 assert allclose(Q.centroid_values, quantity1.centroid_values*3) 696 assert allclose(Q.edge_values, quantity1.edge_values*3) 696 assert allclose(Q.edge_values, quantity1.edge_values*3) 697 697 Q = 3*quantity1 698 698 assert allclose(Q.vertex_values, quantity1.vertex_values*3) … … 704 704 #print quantity1.centroid_values 705 705 #print quantity2.centroid_values 706 706 707 707 assert allclose(Q.vertex_values, 708 708 quantity1.vertex_values * quantity2.vertex_values) … … 712 712 assert allclose(Q.vertex_values, 713 713 4*quantity1.vertex_values + 2) 714 714 715 715 Q = quantity1*quantity2 + 2 716 716 assert allclose(Q.vertex_values, 717 717 quantity1.vertex_values * quantity2.vertex_values + 2) 718 718 719 719 Q = quantity1*quantity2 + quantity3 720 720 assert allclose(Q.vertex_values, 721 721 quantity1.vertex_values * quantity2.vertex_values + 722 quantity3.vertex_values) 722 quantity3.vertex_values) 723 723 Q = quantity1*quantity2 + 3*quantity3 724 724 assert allclose(Q.vertex_values, 725 725 quantity1.vertex_values * quantity2.vertex_values + 726 3*quantity3.vertex_values) 726 3*quantity3.vertex_values) 727 727 Q = quantity1*quantity2 + 3*quantity3 + 5.0 728 728 assert allclose(Q.vertex_values, 729 729 quantity1.vertex_values * quantity2.vertex_values + 730 3*quantity3.vertex_values + 5) 731 730 3*quantity3.vertex_values + 5) 731 732 732 Q = quantity1*quantity2 - quantity3 733 733 assert allclose(Q.vertex_values, 734 734 quantity1.vertex_values * quantity2.vertex_values - 735 quantity3.vertex_values) 735 quantity3.vertex_values) 736 736 Q = 1.5*quantity1*quantity2 - 3*quantity3 + 5.0 737 737 assert allclose(Q.vertex_values, 738 738 1.5*quantity1.vertex_values * quantity2.vertex_values - 739 3*quantity3.vertex_values + 5) 739 3*quantity3.vertex_values + 5) 740 740 741 741 #Try combining quantities and arrays and scalars 742 742 Q = 1.5*quantity1*quantity2.vertex_values -\ 743 3*quantity3.vertex_values + 5.0 743 3*quantity3.vertex_values + 5.0 744 744 assert allclose(Q.vertex_values, 745 745 1.5*quantity1.vertex_values * quantity2.vertex_values - 746 3*quantity3.vertex_values + 5) 747 746 3*quantity3.vertex_values + 5) 747 748 748 749 749 #Powers … … 751 751 assert allclose(Q.vertex_values, quantity1.vertex_values**2) 752 752 753 Q = quantity1**2 +quantity2**2 753 Q = quantity1**2 +quantity2**2 754 754 assert allclose(Q.vertex_values, 755 755 quantity1.vertex_values**2 + \ 756 756 quantity2.vertex_values**2) 757 757 758 Q = (quantity1**2 +quantity2**2)**0.5 758 Q = (quantity1**2 +quantity2**2)**0.5 759 759 assert allclose(Q.vertex_values, 760 760 (quantity1.vertex_values**2 + \ 761 761 quantity2.vertex_values**2)**0.5) 762 763 764 765 766 767 762 763 764 765 766 767 768 768 769 769 def test_gradient(self): … … 791 791 #2 = 4 + a*(-2/3) + b*(-2/3) 792 792 assert allclose(a[0] + b[0], 3) 793 #From orthogonality (a*(y0-y1) + b*(x0-x1) == 0) 794 assert allclose(a[0] - b[0], 0) 793 #From orthogonality (a*(y0-y1) + b*(x0-x1) == 0) 794 assert allclose(a[0] - b[0], 0) 795 795 796 796 797 797 #Right triangle (2) using two point gradient 798 798 #q2 = q1 + a*(x2-x1) + b*(y2-y1) <=> 799 #6 = 4 + a*(4/3) + b*(-2/3) 799 #6 = 4 + a*(4/3) + b*(-2/3) 800 800 assert allclose(2*a[2] - b[2], 3) 801 801 #From orthogonality (a*(y1-y2) + b*(x2-x1) == 0) 802 assert allclose(a[2] + 2*b[2], 0) 802 assert allclose(a[2] + 2*b[2], 0) 803 803 804 804 805 805 #Top triangle (3) using two point gradient 806 806 #q3 = q1 + a*(x3-x1) + b*(y3-y1) <=> 807 #2 = 4 + a*(-2/3) + b*(4/3) 807 #2 = 4 + a*(-2/3) + b*(4/3) 808 808 assert allclose(a[3] - 2*b[3], 3) 809 #From orthogonality (a*(y1-y3) + b*(x3-x1) == 0) 810 assert allclose(2*a[3] + b[3], 0) 809 #From orthogonality (a*(y1-y3) + b*(x3-x1) == 0) 810 assert allclose(2*a[3] + b[3], 0) 811 811 812 812 … … 822 822 #a = 1.2, b=-0.6 823 823 #q(4,0) = 6 + a*(4 - 8/3) + b*(-2/3) 824 assert allclose(quantity.vertex_values[2,2], 8) 824 assert allclose(quantity.vertex_values[2,2], 8) 825 825 826 826 … … 839 839 assert allclose(a[1], 3.0) 840 840 assert allclose(b[1], 1.0) 841 841 842 842 #Work out the others 843 843 844 844 quantity.extrapolate_second_order() 845 845 … … 1052 1052 denom = ones(4, Float)-timestep*sem 1053 1053 1054 x = array([1 , 2, 3, 4]) + array( [.4,.3,.2,.1])1054 x = array([1., 2., 3., 4.]) 1055 1055 x /= denom 1056 x += timestep*array( [4.0, 3.0, 2.0, 1.0] ) 1057 1056 1058 assert allclose( quantity.centroid_values, x) 1057 1059
Note: See TracChangeset
for help on using the changeset viewer.
