Changeset 2650


Ignore:
Timestamp:
Apr 2, 2006, 9:32:18 PM (18 years ago)
Author:
steve
Message:
 
Location:
inundation/pyvolution
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • inundation/pyvolution/test_data_manager.py

    r2648 r2650  
    25812581
    25822582    def test_sww2domain1(self):
    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
     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
    25922592
    25932593        yiel=0.01
    2594         points, vertices, boundary = rectangular(10,10)
    2595 
    2596         #Create shallow water domain
    2597         domain = Domain(points, vertices, boundary)
     2594        points, vertices, boundary = rectangular(10,10)
     2595
     2596        #Create shallow water domain
     2597        domain = Domain(points, vertices, boundary)
    25982598        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
    26362635        import os
    26372636
    2638         filename = domain.datadir+os.sep+domain.filename+'.sww'
    2639         domain2 = sww2domain(filename,None,fail_if_NaN=False,verbose = False)
    2640         #points, vertices, boundary = rectangular(15,15)
     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)
    26412640        #domain2.boundary = boundary
    26422641        ###################
    2643         ##NOW TEST IT!!!
     2642        ##NOW TEST IT!!!
    26442643        ###################
    26452644
     
    26472646        os.remove(filename)
    26482647
    2649         bits = ['vertex_coordinates']
    2650         for quantity in ['elevation']+domain.quantities_to_be_stored:
    2651             bits.append('get_quantity("%s").get_integral()' %quantity)
    2652             bits.append('get_quantity("%s").get_values()' %quantity)
    2653 
    2654         for bit in bits:
    2655             #print 'testing that domain.'+bit+' has been restored'
     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'
    26562655            #print bit
    2657         #print 'done'
    2658             assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
     2656            #print 'done'
     2657            assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
    26592658
    26602659        ######################################
     
    26632662        visualise = False
    26642663        #visualise = True
    2665         domain.visualise = visualise
     2664        domain.visualise = visualise
    26662665        domain.time = 0.
    26672666        from time import sleep
     
    26692668        final = .1
    26702669        domain.set_quantity('friction', 0.1)
    2671         domain.store = False
    2672         domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
    2673 
    2674 
    2675         for t in domain.evolve(yieldstep = yiel, finaltime = final):
     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):
    26762675            if visualise: sleep(1.)
    2677             #domain.write_time()
    2678             pass
     2676            #domain.write_time()
     2677            pass
    26792678
    26802679        final = final - (domain2.starttime-domain.starttime)
     
    26822681        final = final + (domain2.starttime-domain.starttime)
    26832682
    2684         domain2.smooth = False
    2685         domain2.visualise = visualise
    2686         domain2.store = False
    2687         domain2.default_order=2
    2688         domain2.set_quantity('friction', 0.1)
    2689         #Bed-slope and friction
    2690         # Boundary conditions
     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
    26912690        Bd2=Dirichlet_boundary([0.2,0.,0.])
    26922691        domain2.boundary = domain.boundary
    26932692        #print 'domain2.boundary'
    26942693        #print domain2.boundary
    2695         domain2.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
    2696         #domain2.set_boundary({'exterior': Bd})
    2697 
    2698         domain2.check_integrity()
    2699 
    2700         for t in domain2.evolve(yieldstep = yiel, finaltime = final):
     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):
    27012700            if visualise: sleep(1.)
    2702             #domain2.write_time()
    2703             pass
     2701            #domain2.write_time()
     2702            pass
    27042703
    27052704        ###################
    2706         ##NOW TEST IT!!!
     2705        ##NOW TEST IT!!!
    27072706        ##################
    27082707
    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:
    27162716            #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))
    27182720
    27192721
     
    45824584
    45834585
    4584 
  • inundation/pyvolution/test_domain.py

    r2491 r2650  
    139139    def test_create_quantity_from_expression(self):
    140140        """Quantity created from other quantities using arbitrary expression
    141        
     141
    142142        """
    143143
     
    162162        domain.set_quantity('elevation', -1)
    163163
    164        
     164
    165165        domain.set_quantity('stage', [[1,2,3], [5,5,5],
    166166                                      [0,0,9], [-6, 3, 3]])
     
    170170
    171171        domain.set_quantity('ymomentum', [[3,3,3], [4,2,1],
    172                                           [2,4,-1], [1, 0, 1]])       
     172                                          [2,4,-1], [1, 0, 1]])
    173173
    174174        domain.check_integrity()
     
    186186
    187187        X = domain.quantities['xmomentum'].vertex_values
    188         Y = domain.quantities['ymomentum'].vertex_values       
     188        Y = domain.quantities['ymomentum'].vertex_values
    189189
    190190        assert allclose(Q.vertex_values, (X**2 + Y**2)**0.5)
    191                                      
     191
    192192
    193193
     
    196196    def test_set_quantity_from_expression(self):
    197197        """Quantity set using arbitrary expression
    198        
     198
    199199        """
    200200
     
    219219        domain.set_quantity('elevation', -1)
    220220
    221        
     221
    222222        domain.set_quantity('stage', [[1,2,3], [5,5,5],
    223223                                      [0,0,9], [-6, 3, 3]])
     
    227227
    228228        domain.set_quantity('ymomentum', [[3,3,3], [4,2,1],
    229                                           [2,4,-1], [1, 0, 1]])       
     229                                          [2,4,-1], [1, 0, 1]])
    230230
    231231
     
    244244                                      [1,1,10], [-5, 4, 4]])
    245245
    246        
    247 
    248                                      
     246
     247
     248
    249249
    250250
     
    426426        denom = ones(4, Float)-domain.timestep*sem
    427427
    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.])
    429432        x /= denom
     433        x += domain.timestep*array( [4,3,2,1] )
    430434
    431435        for name in domain.conserved_quantities:
  • inundation/pyvolution/test_quantity.py

    r2526 r2650  
    9898        assert allclose(quantity.vertex_values[1,:],[3.+2./3, 6.+2./3, 4.+2./3])
    9999        #FIXME: Work out the others
    100        
    101                                          
     100
     101
    102102        #print quantity.edge_values
    103103        assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
     
    276276
    277277
    278        
     278
    279279
    280280    def test_set_values_using_least_squares(self):
    281281
    282282        from geospatial_data.geospatial_data import Geospatial_data
    283        
     283
    284284        quantity = Quantity(self.mesh4)
    285285
     
    300300
    301301        z = linear_function(data_points)
    302        
     302
    303303        #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 )
    305305        #quantity.set_values(points = data_points, values = z, alpha = 0)
    306306
    307        
     307
    308308        answer = linear_function(quantity.domain.get_vertex_coordinates(obj = True))
    309309        #print quantity.vertex_values, answer
     
    312312
    313313        #Now try by setting the same values directly
    314         from least_squares import fit_to_mesh       
     314        from least_squares import fit_to_mesh
    315315        vertex_attributes = fit_to_mesh(quantity.domain.coordinates,
    316316                                        quantity.domain.triangles,
     
    331331
    332332        from domain import Domain
    333         from geospatial_data.geospatial_data import Geospatial_data       
     333        from geospatial_data.geospatial_data import Geospatial_data
    334334        from coordinate_transforms.geo_reference import Geo_reference
    335335        from least_squares import fit_to_mesh
     
    345345                       geo_reference = mesh_georef)
    346346        mesh1.check_integrity()
    347        
    348         #Quantity                       
     347
     348        #Quantity
    349349        quantity = Quantity(mesh1)
    350350
    351         #Data                       
     351        #Data
    352352        data_points = [[ 201.0, 401.0],
    353353                       [ 201.0, 403.0],
     
    355355
    356356        z = [2, 4, 4]
    357        
     357
    358358        data_georef = Geo_reference(56,-200,-400)
    359        
     359
    360360
    361361        #Reference
     
    364364                          mesh_origin = mesh_georef.get_origin(),
    365365                          alpha = 0)
    366        
    367         assert allclose( ref, [0,5,5] )       
     366
     367        assert allclose( ref, [0,5,5] )
    368368
    369369
    370370        #Test set_values
    371371
    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 )
    373373
    374374        #quantity.set_values(points = data_points,
     
    377377        #                    alpha = 0)
    378378
    379        
     379
    380380        #quantity.set_values(points = data_points,
    381381        #                    values = z,
     
    388388        #Test set_values using geospatial data object
    389389        quantity.vertex_values[:] = 0.0
    390        
     390
    391391        from geospatial_data.geospatial_data import Geospatial_data
    392392        geo = Geospatial_data(data_points, z, data_georef)
    393393
    394                                          
     394
    395395        quantity.set_values(geospatial_data = geo, alpha = 0)
    396         assert allclose(quantity.vertex_values.flat, ref)       
     396        assert allclose(quantity.vertex_values.flat, ref)
    397397
    398398
     
    417417
    418418        z = linear_function(data_points)
    419        
     419
    420420
    421421        #Create pts file
     
    425425        points_dict = {'pointlist': data_points,
    426426                       'attributelist': {att: z}}
    427                      
    428         export_points_file(ptsfile, points_dict) 
     427
     428        export_points_file(ptsfile, points_dict)
    429429
    430430        #Check that values can be set from file
     
    436436        #print answer
    437437
    438        
     438
    439439        assert allclose(quantity.vertex_values.flat, answer)
    440440
     
    442442        #Check that values can be set from file using default attribute
    443443        quantity.set_values(filename = ptsfile, alpha = 0)
    444         assert allclose(quantity.vertex_values.flat, answer)       
     444        assert allclose(quantity.vertex_values.flat, answer)
    445445
    446446        #Cleanup
    447447        import os
    448448        os.remove(ptsfile)
    449        
     449
    450450
    451451    def test_set_values_from_file_with_georef1(self):
    452        
     452
    453453        #Mesh in zone 56 (absolute coords)
    454454        from domain import Domain
    455455        from coordinate_transforms.geo_reference import Geo_reference
    456        
     456
    457457        x0 = 314036.58727982
    458458        y0 = 6224951.2960092
     
    491491
    492492        z = linear_function(data_points)
    493        
     493
    494494
    495495        #Create pts file
    496         from load_mesh.loadASCII import export_points_file       
     496        from load_mesh.loadASCII import export_points_file
    497497
    498498        ptsfile = 'testptsfile.pts'
     
    504504                                                      xllcorner = x0,
    505505                                                      yllcorner = y0)}
    506                      
     506
    507507        export_points_file(ptsfile, points_dict)
    508        
     508
    509509
    510510        #Check that values can be set from file
     
    512512                            attribute_name = att, alpha = 0)
    513513        answer = linear_function(quantity.domain.get_vertex_coordinates(obj = True) - [x0, y0])
    514        
     514
    515515        assert allclose(quantity.vertex_values.flat, answer)
    516516
     
    518518        #Check that values can be set from file using default attribute
    519519        quantity.set_values(filename = ptsfile, alpha = 0)
    520         assert allclose(quantity.vertex_values.flat, answer)       
     520        assert allclose(quantity.vertex_values.flat, answer)
    521521
    522522        #Cleanup
     
    526526
    527527    def test_set_values_from_file_with_georef2(self):
    528        
     528
    529529        #Mesh in zone 56 (relative coords)
    530530        from domain import Domain
    531531        from coordinate_transforms.geo_reference import Geo_reference
    532        
     532
    533533        x0 = 314036.58727982
    534534        y0 = 6224951.2960092
     
    567567
    568568        z = linear_function(data_points)
    569        
     569
    570570
    571571        #Create pts file
    572         from load_mesh.loadASCII import export_points_file       
     572        from load_mesh.loadASCII import export_points_file
    573573
    574574        ptsfile = 'testptsfile.pts'
     
    580580                                                      xllcorner = 0,
    581581                                                      yllcorner = 0)}
    582                      
     582
    583583        export_points_file(ptsfile, points_dict)
    584        
     584
    585585
    586586        #Check that values can be set from file
     
    589589        answer = linear_function(quantity.domain.get_vertex_coordinates(obj = True) + [x0, y0])
    590590
    591        
     591
    592592        assert allclose(quantity.vertex_values.flat, answer)
    593593
     
    595595        #Check that values can be set from file using default attribute
    596596        quantity.set_values(filename = ptsfile, alpha = 0)
    597         assert allclose(quantity.vertex_values.flat, answer)       
     597        assert allclose(quantity.vertex_values.flat, answer)
    598598
    599599        #Cleanup
    600600        import os
    601601        os.remove(ptsfile)
    602        
     602
    603603
    604604
     
    612612                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
    613613
    614        
     614
    615615        quantity2 = Quantity(self.mesh4)
    616616        quantity2.set_values(quantity = quantity1)
     
    621621        assert allclose(quantity2.vertex_values,
    622622                        [[2,0,4], [2,4,8], [8,4,10], [6,2,8]])
    623        
     623
    624624        quantity2.set_values(quantity = 2*quantity1 + 3)
    625625        assert allclose(quantity2.vertex_values,
     
    630630        quantity2.set_values(2*quantity1 + 3)
    631631        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]])
    633633
    634634
     
    644644                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
    645645
    646        
     646
    647647        quantity2 = Quantity(self.mesh4)
    648648        quantity2.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
     
    653653        quantity3 = Quantity(self.mesh4)
    654654        quantity3.set_values([[2,2,2], [7,8,9], [7,6,3], [3, 8, -8]],
    655                              location = 'vertices')       
     655                             location = 'vertices')
    656656
    657657
     
    660660        assert allclose(Q.vertex_values, -quantity1.vertex_values)
    661661        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
    664664        #Addition
    665665        Q = quantity1 + 7
    666666        assert allclose(Q.vertex_values, quantity1.vertex_values + 7)
    667667        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
    670670        Q = 7 + quantity1
    671         assert allclose(Q.vertex_values, quantity1.vertex_values + 7)       
     671        assert allclose(Q.vertex_values, quantity1.vertex_values + 7)
    672672        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
    675675        Q = quantity1 + quantity2
    676676        assert allclose(Q.vertex_values,
     
    679679                        quantity1.centroid_values + quantity2.centroid_values)
    680680        assert allclose(Q.edge_values,
    681                         quantity1.edge_values + quantity2.edge_values)       
    682        
     681                        quantity1.edge_values + quantity2.edge_values)
     682
    683683
    684684        Q = quantity1 + quantity2 - 3
     
    688688        Q = quantity1 - quantity2
    689689        assert allclose(Q.vertex_values,
    690                         quantity1.vertex_values - quantity2.vertex_values)   
     690                        quantity1.vertex_values - quantity2.vertex_values)
    691691
    692692        #Scaling
     
    694694        assert allclose(Q.vertex_values, quantity1.vertex_values*3)
    695695        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)
    697697        Q = 3*quantity1
    698698        assert allclose(Q.vertex_values, quantity1.vertex_values*3)
     
    704704        #print quantity1.centroid_values
    705705        #print quantity2.centroid_values
    706        
     706
    707707        assert allclose(Q.vertex_values,
    708708                        quantity1.vertex_values * quantity2.vertex_values)
     
    712712        assert allclose(Q.vertex_values,
    713713                        4*quantity1.vertex_values + 2)
    714        
     714
    715715        Q = quantity1*quantity2 + 2
    716716        assert allclose(Q.vertex_values,
    717717                        quantity1.vertex_values * quantity2.vertex_values + 2)
    718        
     718
    719719        Q = quantity1*quantity2 + quantity3
    720720        assert allclose(Q.vertex_values,
    721721                        quantity1.vertex_values * quantity2.vertex_values +
    722                         quantity3.vertex_values)       
     722                        quantity3.vertex_values)
    723723        Q = quantity1*quantity2 + 3*quantity3
    724724        assert allclose(Q.vertex_values,
    725725                        quantity1.vertex_values * quantity2.vertex_values +
    726                         3*quantity3.vertex_values)               
     726                        3*quantity3.vertex_values)
    727727        Q = quantity1*quantity2 + 3*quantity3 + 5.0
    728728        assert allclose(Q.vertex_values,
    729729                        quantity1.vertex_values * quantity2.vertex_values +
    730                         3*quantity3.vertex_values + 5)                       
    731        
     730                        3*quantity3.vertex_values + 5)
     731
    732732        Q = quantity1*quantity2 - quantity3
    733733        assert allclose(Q.vertex_values,
    734734                        quantity1.vertex_values * quantity2.vertex_values -
    735                         quantity3.vertex_values)                               
     735                        quantity3.vertex_values)
    736736        Q = 1.5*quantity1*quantity2 - 3*quantity3 + 5.0
    737737        assert allclose(Q.vertex_values,
    738738                        1.5*quantity1.vertex_values * quantity2.vertex_values -
    739                         3*quantity3.vertex_values + 5)                               
     739                        3*quantity3.vertex_values + 5)
    740740
    741741        #Try combining quantities and arrays and scalars
    742742        Q = 1.5*quantity1*quantity2.vertex_values -\
    743             3*quantity3.vertex_values + 5.0 
     743            3*quantity3.vertex_values + 5.0
    744744        assert allclose(Q.vertex_values,
    745745                        1.5*quantity1.vertex_values * quantity2.vertex_values -
    746                         3*quantity3.vertex_values + 5)                                       
    747        
     746                        3*quantity3.vertex_values + 5)
     747
    748748
    749749        #Powers
     
    751751        assert allclose(Q.vertex_values, quantity1.vertex_values**2)
    752752
    753         Q = quantity1**2 +quantity2**2 
     753        Q = quantity1**2 +quantity2**2
    754754        assert allclose(Q.vertex_values,
    755755                        quantity1.vertex_values**2 + \
    756756                        quantity2.vertex_values**2)
    757757
    758         Q = (quantity1**2 +quantity2**2)**0.5 
     758        Q = (quantity1**2 +quantity2**2)**0.5
    759759        assert allclose(Q.vertex_values,
    760760                        (quantity1.vertex_values**2 + \
    761761                        quantity2.vertex_values**2)**0.5)
    762        
    763 
    764 
    765        
    766        
    767        
     762
     763
     764
     765
     766
     767
    768768
    769769    def test_gradient(self):
     
    791791        #2  = 4  + a*(-2/3)  + b*(-2/3)
    792792        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)
    795795
    796796
    797797        #Right triangle (2) using two point gradient
    798798        #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)
    800800        assert allclose(2*a[2] - b[2], 3)
    801801        #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)
    803803
    804804
    805805        #Top triangle (3) using two point gradient
    806806        #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)
    808808        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)
    811811
    812812
     
    822822        #a = 1.2, b=-0.6
    823823        #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)
    825825
    826826
     
    839839        assert allclose(a[1], 3.0)
    840840        assert allclose(b[1], 1.0)
    841        
     841
    842842        #Work out the others
    843        
     843
    844844        quantity.extrapolate_second_order()
    845845
     
    10521052        denom = ones(4, Float)-timestep*sem
    10531053
    1054         x = array([1, 2, 3, 4]) + array( [.4,.3,.2,.1] )
     1054        x = array([1., 2., 3., 4.])
    10551055        x /= denom
     1056        x += timestep*array( [4.0, 3.0, 2.0, 1.0] )
     1057
    10561058        assert allclose( quantity.centroid_values, x)
    10571059
Note: See TracChangeset for help on using the changeset viewer.