Changeset 6553


Ignore:
Timestamp:
Mar 19, 2009, 1:43:34 PM (15 years ago)
Author:
rwilson
Message:

Merged trunk into numpy, all tests and validations work.

Location:
branches/numpy/anuga
Files:
14 added
31 edited

Legend:

Unmodified
Added
Removed

  • branches/numpy/anuga/abstract_2d_finite_volumes/domain.py

    r6533 r6553  
    319319        return self.mesh.get_boundary_polygon(*args, **kwargs)
    320320
     321    # FIXME(Ole): This doesn't seem to be required
    321322    def get_number_of_triangles_per_node(self, *args, **kwargs):
    322323        return self.mesh.get_number_of_triangles_per_node(*args, **kwargs)
    … …  
    348349    def statistics(self, *args, **kwargs):
    349350        return self.mesh.statistics(*args, **kwargs)
     351       
     352    def get_extent(self, *args, **kwargs):
     353        return self.mesh.get_extent(*args, **kwargs)   
    350354
    351355    ##

  • branches/numpy/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py

    r6533 r6553  
    8888    # FIXME (Ole): We should rename f to function to be consistent with
    8989    # Transmissive_Momentum_Set_Stage_Boundary (cf posting by rrraman)
    90     def __init__(self, domain = None, f = None):
     90    def __init__(self, domain=None,
     91                 f=None,
     92                 default_boundary=None,
     93                 verbose=False):
    9194        Boundary.__init__(self)
     95        self.default_boundary = default_boundary
     96        self.default_boundary_invoked = False    # Flag
     97        self.domain = domain
     98        self.verbose = verbose
    9299
    93100        try:
    … …  
    122129    def evaluate(self, vol_id=None, edge_id=None):
    123130        # FIXME (Ole): I think this should be get_time(), see ticket:306
    124         return self.f(self.domain.time)
     131        try:
     132            res = self.f(self.domain.time)
     133        except Modeltime_too_early, e:
     134            raise Modeltime_too_early, e
     135        except Modeltime_too_late, e:
     136            if self.default_boundary is None:
     137                raise Exception, e # Reraise exception
     138            else:
     139                # Pass control to default boundary
     140                res = self.default_boundary.evaluate(vol_id, edge_id)
     141               
     142                # Ensure that result cannot be manipulated
     143                # This is a real danger in case the
     144                # default_boundary is a Dirichlet type
     145                # for instance.
     146                res = res.copy()
     147               
     148                if self.default_boundary_invoked is False:
     149                    if self.verbose:               
     150                        # Issue warning the first time
     151                        msg = '%s' %str(e)
     152                        msg += 'Instead I will use the default boundary: %s\n'\
     153                            %str(self.default_boundary)
     154                        msg += 'Note: Further warnings will be supressed'
     155                        print msg
     156               
     157                    # FIXME (Ole): Replace this crude flag with
     158                    # Python's ability to print warnings only once.
     159                    # See http://docs.python.org/lib/warning-filter.html
     160                    self.default_boundary_invoked = True
     161
     162        return res
    125163
    126164
    … …  
    299337                    if self.default_boundary_invoked is False:
    300338                        # Issue warning the first time
    301                         msg = '%s' %str(e)
    302                         msg += 'Instead I will use the default boundary: %s\n'\
    303                             %str(self.default_boundary)
    304                         msg += 'Note: Further warnings will be supressed'
    305                         warn(msg)
     339                        if self.verbose:
     340                            msg = '%s' %str(e)
     341                            msg += 'Instead I will use the default boundary: %s\n'\
     342                                %str(self.default_boundary)
     343                            msg += 'Note: Further warnings will be supressed'
     344                            print msg
    306345                   
    307346                        # FIXME (Ole): Replace this crude flag with

  • branches/numpy/anuga/abstract_2d_finite_volumes/quantity.py

    r6533 r6553  
    883883
    884884           
    885         # FIXME(Ole): I noticed a couple of examples where this caused
    886         # a crash in fittng, so disabled it until I can investigate further
    887         # Sorry. 23 Jan 2009. Logged as ticket:314
    888         if False:
     885        if True:
    889886            # Use mesh as defined by domain
    890887            # This used to cause problems for caching due to quantities
    891888            # changing, but it now works using the appropriate Mesh object.
    892             # This should close ticket:242
     889            # This addressed ticket:242 and was made to work when bug
     890            # in ticket:314 was fixed 18 March 2009.
    893891            vertex_attributes = fit_to_mesh(filename,
    894892                                            mesh=self.domain.mesh,
    … …  
    901899            # This variant will cause Mesh object to be recreated
    902900            # in fit_to_mesh thus doubling up on the neighbour structure
    903             # FIXME(Ole): This is now obsolete 19 Jan 2009.
     901            # FIXME(Ole): This is now obsolete 19 Jan 2009 except for bug
     902            # (ticket:314) which was fixed 18 March 2009.
    904903            nodes = self.domain.get_nodes(absolute=True)
    905904            triangles = self.domain.get_triangles()
    … …  
    11871186        import types
    11881187
    1189         assert type(indices) in [types.ListType, types.NoneType,
    1190                                  num.ndarray], \
    1191                    'Indices must be a list or None'     #??#
     1188        msg = 'Indices must be a list or None'
     1189        assert indices is None or isinstance(indices, (list, num.ndarray)), msg
     1190#        assert type(indices) in [types.ListType, types.NoneType,
     1191#                                 num.ndarray], \
     1192#                   'Indices must be a list or None'     #??#
    11921193
    11931194        if location == 'centroids':

  • branches/numpy/anuga/abstract_2d_finite_volumes/test_domain.py

    r6533 r6553  
    542542                        conserved_quantities =\
    543543                        ['stage', 'xmomentum', 'ymomentum'])
     544        domain.set_default_order(1)
    544545        domain.check_integrity()
    545546
    … …  
    650651                        conserved_quantities =\
    651652                        ['stage', 'xmomentum', 'ymomentum'])
     653        domain.set_default_order(1)                       
    652654        domain.check_integrity()
    653655
    … …  
    800802                             [ 11.0,  11.0,  11.0],
    801803                             [ 11.0,  11.0,  11.0]])
     804                             
     805    def test_that_mesh_methods_exist(self):
     806        """test_that_mesh_methods_exist
     807       
     808        Test that relavent mesh methods are made available in
     809        domain through composition
     810        """
     811        from mesh_factory import rectangular
     812        from shallow_water import Domain
     813
     814        # Create basic mesh
     815        points, vertices, boundary = rectangular(1, 3)
     816
     817        # Create shallow water domain
     818        domain = Domain(points, vertices, boundary)                             
     819       
     820       
     821        domain.get_centroid_coordinates()
     822        domain.get_radii()
     823        domain.get_areas()
     824        domain.get_area()
     825        domain.get_vertex_coordinates()
     826        domain.get_triangles()
     827        domain.get_nodes()
     828        domain.get_number_of_nodes()
     829        domain.get_normal(0,0)
     830        domain.get_intersecting_segments([[0.0, 0.0], [0.0, 1.0]])
     831        domain.get_disconnected_triangles()
     832        domain.get_boundary_tags()
     833        domain.get_boundary_polygon()
     834        #domain.get_number_of_triangles_per_node()
     835        domain.get_triangles_and_vertices_per_node()
     836        domain.get_interpolation_object()
     837        domain.get_tagged_elements()
     838        domain.get_lone_vertices()
     839        domain.get_unique_vertices()
     840        g = domain.get_georeference()
     841        domain.set_georeference(g)
     842        domain.build_tagged_elements_dictionary()
     843        domain.statistics()
     844        domain.get_extent()
     845
     846       
     847       
    802848
    803849#-------------------------------------------------------------

  • branches/numpy/anuga/abstract_2d_finite_volumes/test_util.py

    r6458 r6553  
    16451645
    16461646#        point1_answers_array = [[0.0,1.0,-5.0,3.0,4.0], [2.0,10.0,-5.0,3.0,4.0]]
    1647         point1_answers_array = [[0.0,1.0,6.0,-5.0,3.0,4.0], [2.0,10.0,15.0,-5.0,3.0,4.0]]
     1647        point1_answers_array = [[0.0,0.0,1.0,6.0,-5.0,3.0,4.0], [2.0,2.0/3600.,10.0,15.0,-5.0,3.0,4.0]]
    16481648        point1_filename = 'gauge_point1.csv'
    16491649        point1_handle = file(point1_filename)
    … …  
    16531653        line=[]
    16541654        for i,row in enumerate(point1_reader):
    1655             #print 'i',i,'row',row
    1656             line.append([float(row[0]),float(row[1]),float(row[2]),float(row[3]),float(row[4]),float(row[5])])
    1657             #print 'assert line',line[i],'point1',point1_answers_array[i]
     1655#            print 'i',i,'row',row
     1656            line.append([float(row[0]),float(row[1]),float(row[2]),float(row[3]),
     1657                         float(row[4]),float(row[5]),float(row[6])])
     1658#            print 'assert line',line[i],'point1',point1_answers_array[i]
    16581659            assert num.allclose(line[i], point1_answers_array[i])
    16591660
    1660         point2_answers_array = [[0.0,1.0,1.5,-0.5,3.0,4.0], [2.0,10.0,10.5,-0.5,3.0,4.0]]
     1661        point2_answers_array = [[0.0,0.0,1.0,1.5,-0.5,3.0,4.0], [2.0,2.0/3600.,10.0,10.5,-0.5,3.0,4.0]]
    16611662        point2_filename = 'gauge_point2.csv'
    16621663        point2_handle = file(point2_filename)
    … …  
    16661667        line=[]
    16671668        for i,row in enumerate(point2_reader):
    1668             #print 'i',i,'row',row
    1669             line.append([float(row[0]),float(row[1]),float(row[2]),float(row[3]),float(row[4]),float(row[5])])
    1670             #print 'assert line',line[i],'point1',point1_answers_array[i]
     1669#            print 'i',i,'row',row
     1670            line.append([float(row[0]),float(row[1]),float(row[2]),float(row[3]),
     1671                         float(row[4]),float(row[5]),float(row[6])])
     1672#            print 'assert line',line[i],'point1',point1_answers_array[i]
    16711673            assert num.allclose(line[i], point2_answers_array[i])
    16721674                         
    … …  
    17731775        for i,row in enumerate(point1_reader):
    17741776#            print 'i',i,'row',row
    1775             line.append([float(row[0]),float(row[1]),float(row[2])])
     1777            # note the 'hole' (element 1) below - skip the new 'hours' field
     1778            line.append([float(row[0]),float(row[2]),float(row[3])])
    17761779            #print 'line',line[i],'point1',point1_answers_array[i]
    17771780            assert num.allclose(line[i], point1_answers_array[i])
    … …  
    17861789        for i,row in enumerate(point2_reader):
    17871790#            print 'i',i,'row',row
    1788             line.append([float(row[0]),float(row[1]),float(row[2])])
     1791            # note the 'hole' (element 1) below - skip the new 'hours' field
     1792            line.append([float(row[0]),float(row[2]),float(row[3])])
    17891793#            print 'line',line[i],'point1',point1_answers_array[i]
    17901794            assert num.allclose(line[i], point2_answers_array[i])
    … …  
    18831887
    18841888#        point1_answers_array = [[0.0,1.0,-5.0,3.0,4.0], [2.0,10.0,-5.0,3.0,4.0]]
    1885         point1_answers_array = [[5.0,1.0,6.0,-5.0,3.0,4.0], [7.0,10.0,15.0,-5.0,3.0,4.0]]
     1889        point1_answers_array = [[5.0,5.0/3600.,1.0,6.0,-5.0,3.0,4.0], [7.0,7.0/3600.,10.0,15.0,-5.0,3.0,4.0]]
    18861890        point1_filename = 'gauge_point1.csv'
    18871891        point1_handle = file(point1_filename)
    … …  
    18921896        for i,row in enumerate(point1_reader):
    18931897            #print 'i',i,'row',row
    1894             line.append([float(row[0]),float(row[1]),float(row[2]),float(row[3]),float(row[4]),float(row[5])])
     1898            line.append([float(row[0]),float(row[1]),float(row[2]),float(row[3]),
     1899                         float(row[4]), float(row[5]), float(row[6])])
    18951900            #print 'assert line',line[i],'point1',point1_answers_array[i]
    18961901            assert num.allclose(line[i], point1_answers_array[i])
    18971902
    1898         point2_answers_array = [[5.0,1.0,1.5,-0.5,3.0,4.0], [7.0,10.0,10.5,-0.5,3.0,4.0]]
     1903        point2_answers_array = [[5.0,5.0/3600.,1.0,1.5,-0.5,3.0,4.0], [7.0,7.0/3600.,10.0,10.5,-0.5,3.0,4.0]]
    18991904        point2_filename = 'gauge_point2.csv'
    19001905        point2_handle = file(point2_filename)
    … …  
    19051910        for i,row in enumerate(point2_reader):
    19061911            #print 'i',i,'row',row
    1907             line.append([float(row[0]),float(row[1]),float(row[2]),float(row[3]),float(row[4]),float(row[5])])
     1912            line.append([float(row[0]),float(row[1]),float(row[2]),float(row[3]),
     1913                         float(row[4]),float(row[5]), float(row[6])])
    19081914            #print 'assert line',line[i],'point1',point1_answers_array[i]
    19091915            assert num.allclose(line[i], point2_answers_array[i])

  • branches/numpy/anuga/abstract_2d_finite_volumes/util.py

    r6533 r6553  
    925925                    verbose = False):   
    926926       
     927    # FIXME(Ole): Shouldn't print statements here be governed by verbose?
    927928    assert type(gauge_filename) == type(''), 'Gauge filename must be a string'
    928929   
    … …  
    971972            raise msg
    972973
    973         print 'swwfile', swwfile
     974        if verbose:
     975            print 'swwfile', swwfile
    974976
    975977        # Extract parent dir name and use as label
    … …  
    13091311            thisfile = file_loc[j] + sep + 'gauges_time_series' + '_' \
    13101312                       + gaugeloc + '.csv'
    1311             fid_out = open(thisfile, 'w')
    1312             s = 'Time, Stage, Momentum, Speed, Elevation, xmom, ymom, Bearing \n'
    1313             fid_out.write(s)
     1313            if j == 0:
     1314                fid_out = open(thisfile, 'w')
     1315                s = 'Time, Stage, Momentum, Speed, Elevation, xmom, ymom, Bearing \n'
     1316                fid_out.write(s)           
    13141317
    13151318            #### generate quantities #######
    … …  
    23682371        file.close()
    23692372
    2370 
    23712373##
    23722374# @brief ??
    2373 # @param sww_file ??
     2375# @param ??
    23742376# @param gauge_file ??
    23752377# @param out_name ??
    … …  
    23832385                               'xmomentum', 'ymomentum'],
    23842386                   verbose=False,
    2385                    use_cache = True):
     2387                   use_cache=True):
    23862388    """
    23872389   
    23882390    Inputs:
    2389        
    23902391        NOTE: if using csv2timeseries_graphs after creating csv file,
    23912392        it is essential to export quantities 'depth' and 'elevation'.
    … …  
    24122413           myfile_2_point1.csv if <out_name> ='myfile_2_'
    24132414           
    2414            
    24152415        They will all have a header
    24162416   
    … …  
    24362436    import string
    24372437    from anuga.shallow_water.data_manager import get_all_swwfiles
    2438 
    2439 #    quantities =  ['stage', 'elevation', 'xmomentum', 'ymomentum']
    2440     #print "points",points
    24412438
    24422439    assert type(gauge_file) == type(''), 'Gauge filename must be a string'
    … …  
    24552452    point_name = []
    24562453   
    2457     #read point info from file
     2454    # read point info from file
    24582455    for i,row in enumerate(point_reader):
    2459         #read header and determine the column numbers to read correcty.
     2456        # read header and determine the column numbers to read correctly.
    24602457        if i==0:
    24612458            for j,value in enumerate(row):
    … …  
    24892486                                 base_name=base,
    24902487                                 verbose=verbose)
     2488    #print 'sww files just after get_all_swwfiles()', sww_files
     2489    # fudge to get SWW files in 'correct' order, oldest on the left
     2490    sww_files.sort()
     2491
     2492    if verbose:
     2493        print 'sww files', sww_files
    24912494   
    24922495    #to make all the quantities lower case for file_function
    … …  
    24972500
    24982501    core_quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum']
    2499    
     2502
     2503    gauge_file = out_name
     2504
     2505    heading = [quantity for quantity in quantities]
     2506    heading.insert(0,'time')
     2507    heading.insert(1,'hours')
     2508
     2509    #create a list of csv writers for all the points and write header
     2510    points_writer = []
     2511    for point_i,point in enumerate(points):
     2512        points_writer.append(writer(file(dir_name + sep + gauge_file
     2513                                         + point_name[point_i] + '.csv', "wb")))
     2514        points_writer[point_i].writerow(heading)
     2515   
     2516    if verbose: print 'Writing csv files'
     2517
     2518    quake_offset_time = None
     2519
    25002520    for sww_file in sww_files:
    25012521        sww_file = join(dir_name, sww_file+'.sww')
    … …  
    25052525                                     verbose=verbose,
    25062526                                     use_cache=use_cache)
    2507     gauge_file = out_name
    2508 
    2509     heading = [quantity for quantity in quantities]
    2510     heading.insert(0,'time')
    2511 
    2512     #create a list of csv writers for all the points and write header
    2513     points_writer = []
    2514     for i,point in enumerate(points):
    2515         points_writer.append(writer(file(dir_name + sep + gauge_file
    2516                                          + point_name[i] + '.csv', "wb")))
    2517         points_writer[i].writerow(heading)
    2518    
    2519     if verbose: print 'Writing csv files'
    2520 
    2521     for time in callable_sww.get_time():
    2522         for point_i, point in enumerate(points_array):
    2523             #add domain starttime to relative time.
    2524             points_list = [time + callable_sww.starttime]
    2525             point_quantities = callable_sww(time,point_i)
    2526            
    2527             for quantity in quantities:
    2528                 if quantity == NAN:
    2529                     print 'quantity does not exist in' % callable_sww.get_name
    2530                 else:
    2531                     if quantity == 'stage':
    2532                         points_list.append(point_quantities[0])
    2533                        
    2534                     if quantity == 'elevation':
    2535                         points_list.append(point_quantities[1])
    2536                        
    2537                     if quantity == 'xmomentum':
    2538                         points_list.append(point_quantities[2])
    2539                        
    2540                     if quantity == 'ymomentum':
    2541                         points_list.append(point_quantities[3])
    2542                        
    2543                     if quantity == 'depth':
    2544                         points_list.append(point_quantities[0]
    2545                                            - point_quantities[1])
    2546 
    2547                     if quantity == 'momentum':
    2548                         momentum = sqrt(point_quantities[2]**2
    2549                                         + point_quantities[3]**2)
    2550                         points_list.append(momentum)
    2551                        
    2552                     if quantity == 'speed':
    2553                         #if depth is less than 0.001 then speed = 0.0
    2554                         if point_quantities[0] - point_quantities[1] < 0.001:
    2555                             vel = 0.0
    2556                         else:
    2557                             if point_quantities[2] < 1.0e6:
    2558                                 momentum = sqrt(point_quantities[2]**2
    2559                                                 + point_quantities[3]**2)
    2560     #                            vel = momentum/depth             
    2561                                 vel = momentum / (point_quantities[0]
    2562                                                   - point_quantities[1])
    2563     #                            vel = momentum/(depth + 1.e-6/depth)
     2527
     2528        if quake_offset_time is None:
     2529            quake_offset_time = callable_sww.starttime
     2530
     2531        for time in callable_sww.get_time():
     2532            for point_i, point in enumerate(points_array):
     2533               # print 'gauge_file = ', str(point_name[point_i])
     2534                #print 'point_i = ', str(point_i), ' point is = ', str(point)
     2535                #add domain starttime to relative time.
     2536                quake_time = time + quake_offset_time
     2537                points_list = [quake_time, quake_time/3600.]# fudge around SWW time bug
     2538                #print 'point list = ', str(points_list)
     2539                point_quantities = callable_sww(time,point_i)
     2540                #print 'point quantities = ', str(point_quantities)
     2541               
     2542                for quantity in quantities:
     2543                    if quantity == NAN:
     2544                        print 'quantity does not exist in' % callable_sww.get_name
     2545                    else:
     2546                        if quantity == 'stage':
     2547                            points_list.append(point_quantities[0])
     2548                           
     2549                        if quantity == 'elevation':
     2550                            points_list.append(point_quantities[1])
     2551                           
     2552                        if quantity == 'xmomentum':
     2553                            points_list.append(point_quantities[2])
     2554                           
     2555                        if quantity == 'ymomentum':
     2556                            points_list.append(point_quantities[3])
     2557                           
     2558                        if quantity == 'depth':
     2559                            points_list.append(point_quantities[0]
     2560                                               - point_quantities[1])
     2561
     2562                        if quantity == 'momentum':
     2563                            momentum = sqrt(point_quantities[2]**2
     2564                                            + point_quantities[3]**2)
     2565                            points_list.append(momentum)
     2566                           
     2567                        if quantity == 'speed':
     2568                            #if depth is less than 0.001 then speed = 0.0
     2569                            if point_quantities[0] - point_quantities[1] < 0.001:
     2570                                vel = 0.0
    25642571                            else:
    2565                                 momentum = 0
    2566                                 vel = 0
     2572                                if point_quantities[2] < 1.0e6:
     2573                                    momentum = sqrt(point_quantities[2]**2
     2574                                                    + point_quantities[3]**2)
     2575                                    vel = momentum / (point_quantities[0]
     2576                                                      - point_quantities[1])
     2577                                else:
     2578                                    momentum = 0
     2579                                    vel = 0
     2580                               
     2581                            points_list.append(vel)
    25672582                           
    2568                         points_list.append(vel)
    2569                        
    2570                     if quantity == 'bearing':
    2571                         points_list.append(calc_bearing(point_quantities[2],
    2572                                                         point_quantities[3]))
    2573 
    2574             points_writer[point_i].writerow(points_list)
    2575        
     2583                        if quantity == 'bearing':
     2584                            points_list.append(calc_bearing(point_quantities[2],
     2585                                                            point_quantities[3]))
     2586
     2587                points_writer[point_i].writerow(points_list)
    25762588
    25772589##

  • branches/numpy/anuga/caching/test_caching.py

    r6533 r6553  
    386386
    387387        # Check that hash value of callable objects don't change
     388        # FIXME (Ole): The hash values do appear to change when OS
     389        # and/or dependencies are upgraded
    388390        if os.name == 'posix' and os.uname()[4] in ['x86_64', 'ia64']:
    389391          # 64 bit hash values

  • branches/numpy/anuga/compile_all.py

    r6533 r6553  
    2727#execfile('test_all.py')
    2828   
     29if sys.platform == 'win32':
     30    raw_input('Press any key')

  • branches/numpy/anuga/config.py

    r6533 r6553  
    9191       
    9292# FIXME (Ole) Maybe get rid of order altogether and use beta_w
    93 # ... and isn't it about time we make the default 2?
    94 default_order = 1
     93default_order = 2
    9594
    9695################################################################################

  • branches/numpy/anuga/coordinate_transforms/geo_reference.py

    r6442 r6553  
    9494        self.units = units
    9595        self.xllcorner = xllcorner
    96         self.yllcorner = yllcorner
    97         #self.is_absolute = num.allclose([self.xllcorner, self.yllcorner], 0)
    98 
     96        self.yllcorner = yllcorner       
     97           
    9998        if NetCDFObject is not None:
    10099            self.read_NetCDF(NetCDFObject)
    … …  
    103102            self.read_ASCII(ASCIIFile, read_title=read_title)
    104103
    105     ##
    106     # @brief Get the X coordinate of the origin of this georef.
     104           
     105        # Set flag for absolute points (used by get_absolute)   
     106        self.absolute = num.allclose([self.xllcorner, self.yllcorner], 0)
     107           
     108
    107109    def get_xllcorner(self):
    108110        return self.xllcorner
    … …  
    234236            self.yllcorner = self.yllcorner[0]
    235237
     238        assert (type(self.xllcorner) == types.FloatType)
     239        assert (type(self.yllcorner) == types.FloatType)
     240        assert (type(self.zone) == types.IntType)
     241
    236242################################################################################
    237243
    … …  
    243249    # @note If 'points' is a list then a changed list is returned.
    244250    def change_points_geo_ref(self, points, points_geo_ref=None):
    245         """Change the geo reference of a list or numeric array of points to
     251        """Change the geo reference of a list or Numeric array of points to
    246252        be this reference.(The reference used for this object)
    247         If the points do not have a geo ref, assume 'absolute' input values
    248         """
    249 
     253        If the points do not have a geo ref, assume 'absolute' values
     254        """
     255        import copy
     256       
    250257        # remember if we got a list
    251258        is_list = isinstance(points, list)
    252259
    253         points = ensure_numeric(copy.copy(points), num.float)
    254 
    255         # sanity checks
     260        #points = ensure_numeric(copy.copy(points), num.float)
     261        points = ensure_numeric(points, num.float)
     262
     263        # sanity checks
     264        if len(points.shape) == 1:
     265            #One point has been passed
     266            msg = 'Single point must have two elements'
     267            assert len(points) == 2, msg
     268            points = num.reshape(points, (1,2))
     269
     270        msg = 'Points array must be two dimensional.\n'
     271        msg += 'I got %d dimensions' %len(points.shape)
     272        assert len(points.shape) == 2, msg
     273
     274        msg = 'Input must be an N x 2 array or list of (x,y) values. '
     275        msg += 'I got an %d x %d array' %points.shape   
     276        assert points.shape[1] == 2, msg               
     277
     278        # FIXME (Ole): Could also check if zone, xllcorner, yllcorner
     279        # are identical in the two geo refs.   
     280        if points_geo_ref is not self:
     281            # If georeferences are different
     282            points = copy.copy(points) # Don't destroy input                   
     283            if not points_geo_ref is None:
     284                # Convert points to absolute coordinates
     285                points[:,0] += points_geo_ref.xllcorner
     286                points[:,1] += points_geo_ref.yllcorner
     287       
     288            # Make points relative to primary geo reference
     289            points[:,0] -= self.xllcorner
     290            points[:,1] -= self.yllcorner
     291
     292        if is_list:
     293            points = points.tolist()
     294           
     295        return points
     296
     297    def is_absolute(self):
     298        """Return True if xllcorner==yllcorner==0 indicating that points
     299        in question are absolute.
     300        """
     301       
     302        # FIXME(Ole): It is unfortunate that decision about whether points
     303        # are absolute or not lies with the georeference object. Ross pointed this out.
     304        # Moreover, this little function is responsible for a large fraction of the time
     305        # using in data fitting (something in like 40 - 50%.
     306        # This was due to the repeated calls to allclose.
     307        # With the flag method fitting is much faster (18 Mar 2009).
     308
     309        # FIXME(Ole): HACK to be able to reuse data already cached (18 Mar 2009).
     310        # Remove at some point
     311        if not hasattr(self, 'absolute'):
     312            self.absolute = num.allclose([self.xllcorner, self.yllcorner], 0)
     313           
     314        # Return absolute flag   
     315        return self.absolute
     316
     317    def get_absolute(self, points):
     318        """Given a set of points geo referenced to this instance,
     319        return the points as absolute values.
     320        """
     321
     322        is_list = False
     323        if type(points) == types.ListType:
     324            is_list = True
     325
     326        points = ensure_numeric(points, num.float)
    256327        if len(points.shape) == 1:
    257328            # One point has been passed
    258329            msg = 'Single point must have two elements'
    259             assert len(points) == 2, msg
    260             points = num.reshape(points, (1,2))
    261 
    262         msg = ('Points array must be two dimensional.\n'
    263                'I got %d dimensions' % len(points.shape))
    264         assert len(points.shape) == 2, msg
    265 
    266         msg = ('Input must be an N x 2 array or list of (x,y) values. '
    267                'I got an %d x %d array' % points.shape)
    268         assert points.shape[1] == 2, msg
    269 
    270         # FIXME (Ole): Could also check if zone, xllcorner, yllcorner
    271         # are identical in the two geo refs.
    272         if points_geo_ref is not self:
    273             # If georeferences are different
    274             if not points_geo_ref is None:
    275                 # Convert points to absolute coordinates
    276                 points[:,0] += points_geo_ref.xllcorner
    277                 points[:,1] += points_geo_ref.yllcorner
    278 
    279             # Make points relative to primary geo reference
    280             points[:,0] -= self.xllcorner
    281             points[:,1] -= self.yllcorner
    282 
    283         # if we got a list, return a list
     330            if not len(points) == 2:
     331                raise ShapeError, msg   
     332
     333
     334        msg = 'Input must be an N x 2 array or list of (x,y) values. '
     335        msg += 'I got an %d x %d array' %points.shape   
     336        if not points.shape[1] == 2:
     337            raise ShapeError, msg   
     338           
     339       
     340        # Add geo ref to points
     341        if not self.is_absolute():
     342            points = copy.copy(points) # Don't destroy input                   
     343            points[:,0] += self.xllcorner
     344            points[:,1] += self.yllcorner
     345
     346       
    284347        if is_list:
    285348            points = points.tolist()
    286 
    287         return points
    288 
    289     ##
    290     # @brief Is this georef absolute?
    291     # @return True if ???
    292     def is_absolute(self):
    293         """Return True if xllcorner==yllcorner==0"""
    294 
    295         return num.allclose([self.xllcorner, self.yllcorner], 0)
    296 
    297     ##
    298     # @brief Convert points to absolute points in this georef instance.
    299     # @param points
    300     # @return
    301     # @note This method changes the input points!
    302     def get_absolute(self, points):
    303         """Given a set of points geo referenced to this instance
    304 
    305         return the points as absolute values.
    306         """
    307 
    308         # remember if we got a list
    309         is_list = isinstance(points, list)
    310 
    311         # convert to numeric array, force a copy
    312         points = ensure_numeric(copy.copy(points), num.float)
    313 
    314         # sanity checks
    315         if len(points.shape) == 1:
    316             #One point has been passed
    317             if not len(points) == 2:
    318                 msg = 'Single point must have two elements'
    319                 raise ShapeError, msg
    320 
    321         if not points.shape[1] == 2:
    322             msg = ('Input must be an N x 2 array or list of (x,y) values. '
    323                    'I got an %d x %d array' % points.shape)
    324             raise ShapeError, msg
    325 
    326         # Add geo ref to points
    327         if not self.is_absolute():
    328             points[:,0] += self.xllcorner
    329             points[:,1] += self.yllcorner
    330 
    331         # if we got a list, return a list
    332         if is_list:
    333             points = points.tolist()
    334 
     349             
    335350        return points
    336351
    … …  
    350365        is_list = isinstance(points, list)
    351366
    352         # convert to a numeric array, force a copy
    353         points = ensure_numeric(copy.copy(points), num.float)
    354 
    355         # sanity checks
     367        points = ensure_numeric(points, num.float)
    356368        if len(points.shape) == 1:
    357369            #One point has been passed
     370            msg = 'Single point must have two elements'
    358371            if not len(points) == 2:
    359                 msg = 'Single point must have two elements'
    360                 raise ShapeError, msg
     372                raise ShapeError, msg   
    361373
    362374        if not points.shape[1] == 2:
    363375            msg = ('Input must be an N x 2 array or list of (x,y) values. '
    364376                   'I got an %d x %d array' % points.shape)
    365             raise ShapeError, msg
     377            raise ShapeError, msg   
    366378
    367379        # Subtract geo ref from points
    368380        if not self.is_absolute():
    369             points[:,0] -= self.xllcorner
     381            points = copy.copy(points) # Don't destroy input                           
     382            points[:,0] -= self.xllcorner
    370383            points[:,1] -= self.yllcorner
    371384
    372         # if we got a list, return a list
    373385        if is_list:
    374386            points = points.tolist()
    375 
     387             
    376388        return points
    377389

  • branches/numpy/anuga/coordinate_transforms/redfearn.py

    r6533 r6553  
    3737    return sign*dd, mm, ss
    3838
    39 def redfearn(lat, lon, false_easting=None, false_northing=None, zone=None):
     39def redfearn(lat, lon, false_easting=None, false_northing=None,
     40             zone=None, central_meridian=None, scale_factor=None):
    4041    """Compute UTM projection using Redfearn's formula
    4142
    … …  
    4748    If zone is specified reproject lat and long to specified zone instead of
    4849    standard zone
     50    If meridian is specified, reproject lat and lon to that instead of zone. In this case
     51    zone will be set to -1 to indicate non-UTM projection
     52
     53    Note that zone and meridian cannot both be specifed
    4954    """
    5055
    … …  
    5762    a = 6378137.0                       #Semi major axis
    5863    inverse_flattening = 298.257222101  #1/f
    59     K0 = 0.9996                         #Central scale factor   
     64    if scale_factor is None:
     65        K0 = 0.9996                         #Central scale factor   
     66    else:
     67        K0 = scale_factor
     68    #print 'scale', K0
    6069    zone_width = 6                      #Degrees
    6170
    … …  
    138147    m = term1 + term2 + term3 + term4 #OK
    139148
    140     #Zone
     149    if zone is not None and central_meridian is not None:
     150        msg = 'You specified both zone and central_meridian. Provide only one of them'
     151        raise Exception, msg
     152   
     153    # Zone
    141154    if zone is None:
    142155        zone = int((lon - longitude_of_western_edge_zone0)/zone_width)
    143156
    144     central_meridian = zone*zone_width+longitude_of_central_meridian_zone0
     157    # Central meridian
     158    if central_meridian is None:
     159        central_meridian = zone*zone_width+longitude_of_central_meridian_zone0
     160    else:
     161        zone = -1
    145162
    146163    omega = (lon-central_meridian)*pi/180 #Relative longitude (radians)

  • branches/numpy/anuga/coordinate_transforms/test_redfearn.py

    r6533 r6553  
    1111from anuga.utilities.anuga_exceptions import ANUGAError
    1212
     13from anuga.utilities.system_tools import get_pathname_from_package
     14from os.path import join
    1315import numpy as num
    1416
    … …  
    122124
    123125    def test_UTM_6_nonstandard_projection(self):
    124         #Test 6 (Geraldton, WA)
    125 
    126         #First test native projection (zone 50)
     126        """test_UTM_6_nonstandard_projection
     127
     128        Test that projections can be forced to
     129        use other than native zone.
     130
     131        Data is from Geraldton, WA
     132        """
     133       
     134
     135        # First test native projection (zone 50)
    127136        zone, easting, northing = redfearn(-29.233299999,114.05)
    128137
    … …  
    267276    #    #assert allclose(northing, 6181725.1724276)
    268277
     278   
     279    def test_nonstandard_meridian_coinciding_with_native(self):
     280        """test_nonstandard_meridian_coinciding_with_native
     281
     282        This test will verify that redfearn can be used to project
     283        points using an arbitrary central meridian that happens to
     284        coincide with the standard meridian at the center of a UTM zone.
     285        This is a preliminary test before testing this functionality
     286        with a truly arbitrary non-standard meridian.
     287        """
     288
     289        # The file projection_test_points_z53.csv contains 10 points
     290        # which straddle the boundary between UTM zones 53 and 54.
     291        # They have been projected to zone 53 irrespective of where they
     292        # belong.
     293
     294        path = get_pathname_from_package('anuga.coordinate_transforms')
     295       
     296        for forced_zone in [53, 54]:
     297       
     298            datafile = join(path, 'projection_test_points_z%d.csv' % forced_zone)
     299            fid = open(datafile)
     300
     301            for line in fid.readlines()[1:]:
     302                fields = line.strip().split(',')
     303               
     304                lon = float(fields[1])
     305                lat = float(fields[2])
     306                x = float(fields[3])
     307                y = float(fields[4])           
     308
     309                zone, easting, northing = redfearn(lat, lon,
     310                                                   zone=forced_zone)
     311               
     312                # Check calculation
     313                assert zone == forced_zone
     314                assert num.allclose(x, easting)
     315                assert num.allclose(y, northing)
     316
     317   
     318   
     319   
     320    def test_nonstandard_meridian(self):
     321        """test_nonstandard_meridian
     322
     323        This test will verify that redfearn can be used to project
     324        points using an arbitrary central meridian.
     325        """
     326
     327        # The file projection_test_points.csv contains 10 points
     328        # which straddle the boundary between UTM zones 53 and 54.
     329        # They have been projected using a central meridian of 137.5
     330        # degrees (the boundary is 138 so it is pretty much right
     331        # in the middle of zones 53 and 54).
     332
     333        path = get_pathname_from_package('anuga.coordinate_transforms')
     334        datafile = join(path, 'projection_test_points.csv')
     335        fid = open(datafile)
     336
     337        for line in fid.readlines()[1:]:
     338            fields = line.strip().split(',')
     339
     340            lon = float(fields[1])
     341            lat = float(fields[2])
     342            x = float(fields[3])
     343            y = float(fields[4])           
     344
     345            zone, easting, northing = redfearn(lat, lon,
     346                                               central_meridian=137.5,
     347                                               scale_factor=0.9996)
     348
     349            assert zone == -1 # Indicates non UTM projection
     350            assert num.allclose(x, easting)
     351            assert num.allclose(y, northing)
     352
     353        # Test that zone and meridian can't both be specified
     354        try:
     355            zone, easting, northing = redfearn(lat, lon,
     356                                               zone=50,
     357                                               central_meridian=137.5)
     358        except:
     359            pass
     360        else:
     361            msg = 'Should have raised exception'
     362            raise Exception, msg
     363
     364           
    269365    def test_convert_lats_longs(self):
    270366

  • branches/numpy/anuga/culvert_flows/Test_Culvert_Flat_Water_Lev.py

    r6533 r6553  
    123123                       number_of_barrels=1,
    124124                       update_interval=2,
     125                       log_file=True,
     126                       discharge_hydrograph=True,
    125127                       verbose=True)
    126128

  • branches/numpy/anuga/culvert_flows/culvert_class.py

    r6533 r6553  
    204204            log_to_file(self.log_filename, description)
    205205            log_to_file(self.log_filename, self.culvert_type)       
     206        else:
     207            self.log_filename = None
    206208
    207209
    … …  
    296298
    297299        # Print some diagnostics to log if requested
    298         if hasattr(self, 'log_filename'):
     300        if self.log_filename is not None:
    299301            s = 'Culvert Effective Length = %.2f m' %(self.length)
    300302            log_to_file(self.log_filename, s)
    … …  
    324326                update = True
    325327
    326         if hasattr(self, 'log_filename'):           
     328           
     329        if self.log_filename is not None:       
    327330            s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
    328331            log_to_file(self.log_filename, s)
    … …  
    439442            if self.verbose is True:
    440443                print msg
    441             if hasattr(self, 'log_filename'):                   
     444               
     445            if self.log_filename is not None:               
    442446                log_to_file(self.log_filename, msg)
    443447       
    … …  
    524528                print '%.2fs - WARNING: Flow is running uphill.' % time
    525529           
    526         if hasattr(self, 'log_filename'):
     530        if self.log_filename is not None:
    527531            s = 'Time=%.2f, inlet stage = %.2f, outlet stage = %.2f'\
    528532                %(time, self.inlet.stage, self.outlet.stage)
    … …  
    573577                    msg += 'for culvert "%s"' % self.label
    574578                   
    575                     if hasattr(self, 'log_filename'):                   
     579                    if self.log_filename is not None:                   
    576580                        log_to_file(self.log_filename, msg)
    577581            else:
    578582                # User culvert routine
    579583                Q, barrel_velocity, culvert_outlet_depth =\
    580                     self.culvert_routine(self, delta_total_energy, g)
     584                    self.culvert_routine(inlet.depth,
     585                                         outlet.depth,
     586                                         inlet.specific_energy,
     587                                         delta_total_energy,
     588                                         g,
     589                                         culvert_length=self.length,
     590                                         culvert_width=self.width,
     591                                         culvert_height=self.height,
     592                                         culvert_type=self.culvert_type,
     593                                         manning=self.manning,
     594                                         sum_loss=self.sum_loss,
     595                                         log_filename=self.log_filename)
    581596           
    582597           
    … …  
    601616
    602617                           
    603    
     618# OBSOLETE (Except for momentum jet in Culvert_flow_energy)   
    604619class Culvert_flow_rating:
    605620    """Culvert flow - transfer water from one hole to another

  • branches/numpy/anuga/culvert_flows/culvert_routines.py

    r6533 r6553  
    88
    99#NOTE:
    10 # Inlet control:  Delta_Et > Es at the inlet
    11 # Outlet control: Delta_Et < Es at the inlet
    12 # where Et is total energy (w + 0.5*v^2/g) and
    13 # Es is the specific energy (h + 0.5*v^2/g)
    14 
    15 
    16 
    17 #   NEW DEFINED CULVERT FLOW---- Flow from INLET 1 ------> INLET 2 (Outlet)
    18 #
    19 # The First Attempt has a Simple Horizontal Circle as a Hole on the Bed
    20 # Flow Is Removed at a Rate of INFLOW
    21 # Downstream there is a similar Circular Hole on the Bed where INFLOW effectively Surcharges
    22 #
    23 # This SHould be changed to a Vertical Opening Both BOX and Circular
    24 # There will be several Culvert Routines such as:
    25 # CULVERT_Boyd_Channel
    26 # CULVERT_Orifice_and_Weir
    27 # CULVERT_Simple_FLOOR
    28 # CULVERT_Simple_WALL
    29 # CULVERT_Eqn_Floor
    30 # CULVERT_Eqn_Wall
    31 # CULVERT_Tab_Floor
    32 # CULVERT_Tab_Wall
    33 # BRIDGES.....
    34 # NOTE NEED TO DEVELOP CONCEPT 1D Model for Linked Pipe System !!!!
    35 
    36 #  COULD USE EPA SWMM Model !!!!
     10# Inlet control:  Delta_total_energy > inlet_specific_energy
     11# Outlet control: Delta_total_energy < inlet_specific_energy
     12# where total energy is (w + 0.5*v^2/g) and
     13# specific energy is (h + 0.5*v^2/g)
    3714
    3815
    … …  
    4017
    4118
    42 def boyd_generalised_culvert_model(culvert, delta_total_energy, g):
    43 
    44     """Boyd's generalisation of the US department of transportation culvert model
    45         # == The quantity of flow passing through a culvert is controlled by many factors
    46         # == It could be that the culvert is controled by the inlet only, with it being Un submerged this is effectively equivalent to the WEIR Equation
    47         # == Else the culvert could be controlled by the inlet, with it being Submerged, this is effectively the Orifice Equation
    48         # == Else it may be controlled by Down stream conditions where depending on the down stream depth, the momentum in the culvert etc. flow is controlled
     19def boyd_generalised_culvert_model(inlet_depth,
     20                                   outlet_depth,
     21                                   inlet_specific_energy,
     22                                   delta_total_energy,
     23                                   g,
     24                                   culvert_length=0.0,
     25                                   culvert_width=0.0,
     26                                   culvert_height=0.0,
     27                                   culvert_type='box',
     28                                   manning=0.0,
     29                                   sum_loss=0.0,
     30                                   log_filename=None):
     31
     32    """Boyd's generalisation of the US department of transportation culvert
     33    model
     34     
     35    The quantity of flow passing through a culvert is controlled by many factors
     36    It could be that the culvert is controlled by the inlet only, with it
     37    being unsubmerged this is effectively equivalent to the weir Equation
     38    Else the culvert could be controlled by the inlet, with it being
     39    submerged, this is effectively the Orifice Equation
     40    Else it may be controlled by down stream conditions where depending on
     41    the down stream depth, the momentum in the culvert etc. flow is controlled
    4942    """
    5043
    … …  
    5346    from anuga.utilities.numerical_tools import safe_acos as acos
    5447
    55     inlet = culvert.inlet
    56     outlet = culvert.outlet       
    57     Q_outlet_tailwater = 0.0
    58     inlet.rate = 0.0
    59     outlet.rate = 0.0
    60     Q_inlet_unsubmerged = 0.0
    61     Q_inlet_submerged = 0.0
    62     Q_outlet_critical_depth = 0.0
    63 
    64     if hasattr(culvert, 'log_filename'):                               
    65         log_filename = culvert.log_filename
    66 
    67     manning = culvert.manning
    68     sum_loss = culvert.sum_loss
    69     length = culvert.length
    70 
    71     if inlet.depth > 0.01:
    72 
    73         E = inlet.specific_energy
    74 
    75         if hasattr(culvert, 'log_filename'):                           
    76             s = 'Specific energy  = %f m' %E
     48
     49    if inlet_depth > 0.01:
     50        # Water has risen above inlet
     51       
     52        if log_filename is not None:                           
     53            s = 'Specific energy  = %f m' % inlet_specific_energy
    7754            log_to_file(log_filename, s)
    7855       
    79         msg = 'Specific energy is negative'
    80         assert E >= 0.0, msg
     56        msg = 'Specific energy at inlet is negative'
     57        assert inlet_specific_energy >= 0.0, msg
    8158                     
    82        
    83         # Water has risen above inlet
    84         if culvert.culvert_type == 'circle':
    85             # Round culvert
     59
     60        if culvert_type == 'circle':
     61            # Round culvert (use width as diameter)
     62            diameter = culvert_width           
    8663
    8764            # Calculate flows for inlet control           
    88             diameter = culvert.diameter
    89 
    90             Q_inlet_unsubmerged = 0.421*g**0.5*diameter**0.87*E**1.63    # Inlet Ctrl Inlet Unsubmerged
    91             Q_inlet_submerged = 0.530*g**0.5*diameter**1.87*E**0.63      # Inlet Ctrl Inlet Submerged
    92 
    93             if hasattr(culvert, 'log_filename'):
    94                 s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' %(Q_inlet_unsubmerged, Q_inlet_submerged)
     65            Q_inlet_unsubmerged = 0.421*g**0.5*diameter**0.87*inlet_specific_energy**1.63 # Inlet Ctrl Inlet Unsubmerged
     66            Q_inlet_submerged = 0.530*g**0.5*diameter**1.87*inlet_specific_energy**0.63   # Inlet Ctrl Inlet Submerged
     67
     68            if log_filename is not None:                                       
     69                s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' % (Q_inlet_unsubmerged, Q_inlet_submerged)
    9570                log_to_file(log_filename, s)
    9671
    97             case = ''
     72
     73            # FIXME(Ole): Are these functions really for inlet control?                   
    9874            if Q_inlet_unsubmerged < Q_inlet_submerged:
    9975                Q = Q_inlet_unsubmerged
    100 
    101                 alpha = acos(1 - inlet.depth/diameter)
     76                alpha = acos(1 - inlet_depth/diameter)
    10277                flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))
    103                 outlet_culvert_depth = inlet.depth
    104                 width = diameter*sin(alpha)
    105                 #perimeter = alpha*diameter               
     78                outlet_culvert_depth = inlet_depth
    10679                case = 'Inlet unsubmerged'
    10780            else:   
    … …  
    10982                flow_area = (diameter/2)**2 * pi
    11083                outlet_culvert_depth = diameter
    111                 width = diameter
    112                 #perimeter = diameter
    11384                case = 'Inlet submerged'                   
    11485
    11586
    11687
    117             if delta_total_energy < E:
     88            if delta_total_energy < inlet_specific_energy:
    11889                # Calculate flows for outlet control
     90               
    11991                # Determine the depth at the outlet relative to the depth of flow in the Culvert
    120 
    121                 if outlet.depth > diameter:        # The Outlet is Submerged
     92                if outlet_depth > diameter:        # The Outlet is Submerged
    12293                    outlet_culvert_depth=diameter
    12394                    flow_area = (diameter/2)**2 * pi  # Cross sectional area of flow in the culvert
    12495                    perimeter = diameter * pi
    125                     width = diameter
    12696                    case = 'Outlet submerged'
    127                 elif outlet.depth==0.0:
    128                     outlet_culvert_depth=inlet.depth  # For purpose of calculation assume the outlet depth = the inlet depth
    129                     alpha = acos(1 - inlet.depth/diameter)
     97                elif outlet_depth==0.0:
     98                    outlet_culvert_depth=inlet_depth # For purpose of calculation assume the outlet depth = the inlet depth
     99                    alpha = acos(1 - inlet_depth/diameter)
    130100                    flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))
    131101                    perimeter = alpha*diameter
    132                     width = diameter*sin(alpha)                   
    133102
    134103                    case = 'Outlet depth is zero'                       
    135104                else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
    136                     outlet_culvert_depth=outlet.depth   # For purpose of calculation assume the outlet depth = the inlet depth
    137                     alpha = acos(1 - outlet.depth/diameter)
     105                    outlet_culvert_depth=outlet_depth   # For purpose of calculation assume the outlet depth = the inlet depth
     106                    alpha = acos(1 - outlet_depth/diameter)
    138107                    flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))
    139108                    perimeter = alpha*diameter
    140                     width = diameter*sin(alpha)                   
    141109                    case = 'Outlet is open channel flow'
    142110
    143111                hyd_rad = flow_area/perimeter
    144112               
    145                 if hasattr(culvert, 'log_filename'):
     113                if log_filename is not None:
    146114                    s = 'hydraulic radius at outlet = %f' %hyd_rad
    147115                    log_to_file(log_filename, s)
    148116
    149117                # Outlet control velocity using tail water
    150                 culvert_velocity = sqrt(delta_total_energy/((sum_loss/2*g)+(manning**2*length)/hyd_rad**1.33333))
     118                culvert_velocity = sqrt(delta_total_energy/((sum_loss/2*g)+(manning**2*culvert_length)/hyd_rad**1.33333))
    151119                Q_outlet_tailwater = flow_area * culvert_velocity
    152120               
    153                 if hasattr(culvert, 'log_filename'):
     121                if log_filename is not None:               
    154122                    s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
    155123                    log_to_file(log_filename, s)
    … …  
    165133
    166134            # Calculate flows for inlet control
    167             height = culvert.height
    168             width = culvert.width           
     135            height = culvert_height
     136            width = culvert_width           
    169137           
    170             Q_inlet_unsubmerged = 0.540*g**0.5*width*E**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
    171             Q_inlet_submerged = 0.702*g**0.5*width*height**0.89*E**0.61  # Flow based on Inlet Ctrl Inlet Submerged
    172 
    173             if hasattr(culvert, 'log_filename'):
     138            Q_inlet_unsubmerged = 0.540*g**0.5*width*inlet_specific_energy**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
     139            Q_inlet_submerged = 0.702*g**0.5*width*height**0.89*inlet_specific_energy**0.61  # Flow based on Inlet Ctrl Inlet Submerged
     140
     141            if log_filename is not None:                           
    174142                s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' %(Q_inlet_unsubmerged, Q_inlet_submerged)
    175143                log_to_file(log_filename, s)
    176144
    177             case = ''
     145
     146            # FIXME(Ole): Are these functions really for inlet control?   
    178147            if Q_inlet_unsubmerged < Q_inlet_submerged:
    179148                Q = Q_inlet_unsubmerged
    180                 flow_area = width*inlet.depth
    181                 outlet_culvert_depth = inlet.depth
    182                 #perimeter=(width+2.0*inlet.depth)               
     149                flow_area = width*inlet_depth
     150                outlet_culvert_depth = inlet_depth
    183151                case = 'Inlet unsubmerged'
    184152            else:   
    … …  
    186154                flow_area = width*height
    187155                outlet_culvert_depth = height
    188                 #perimeter=2.0*(width+height)               
    189156                case = 'Inlet submerged'                   
    190157
    191             if delta_total_energy < E:
     158            if delta_total_energy < inlet_specific_energy:
    192159                # Calculate flows for outlet control
     160               
    193161                # Determine the depth at the outlet relative to the depth of flow in the Culvert
    194 
    195                 if outlet.depth > height:        # The Outlet is Submerged
     162                if outlet_depth > height:        # The Outlet is Submerged
    196163                    outlet_culvert_depth=height
    197164                    flow_area=width*height       # Cross sectional area of flow in the culvert
    198165                    perimeter=2.0*(width+height)
    199166                    case = 'Outlet submerged'
    200                 elif outlet.depth==0.0:
    201                     outlet_culvert_depth=inlet.depth   # For purpose of calculation assume the outlet depth = the inlet depth
    202                     flow_area=width*inlet.depth
    203                     perimeter=(width+2.0*inlet.depth)
     167                elif outlet_depth==0.0:
     168                    outlet_culvert_depth=inlet_depth   # For purpose of calculation assume the outlet depth = the inlet depth
     169                    flow_area=width*inlet_depth
     170                    perimeter=(width+2.0*inlet_depth)
    204171                    case = 'Outlet depth is zero'                       
    205172                else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
    206                     outlet_culvert_depth=outlet.depth
    207                     flow_area=width*outlet.depth
    208                     perimeter=(width+2.0*outlet.depth)
     173                    outlet_culvert_depth=outlet_depth
     174                    flow_area=width*outlet_depth
     175                    perimeter=(width+2.0*outlet_depth)
    209176                    case = 'Outlet is open channel flow'
    210177
    211178                hyd_rad = flow_area/perimeter
    212179               
    213                 if hasattr(culvert, 'log_filename'):               
    214                     s = 'hydraulic radius at outlet = %f' %hyd_rad
     180                if log_filename is not None:                               
     181                    s = 'hydraulic radius at outlet = %f' % hyd_rad
    215182                    log_to_file(log_filename, s)
    216183
    217184                # Outlet control velocity using tail water
    218                 culvert_velocity = sqrt(delta_total_energy/((sum_loss/2*g)+(manning**2*length)/hyd_rad**1.33333))
     185                culvert_velocity = sqrt(delta_total_energy/((sum_loss/2*g)+(manning**2*culvert_length)/hyd_rad**1.33333))
    219186                Q_outlet_tailwater = flow_area * culvert_velocity
    220187
    221                 if hasattr(culvert, 'log_filename'):                           
    222                     s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
     188                if log_filename is not None:                           
     189                    s = 'Q_outlet_tailwater = %.6f' % Q_outlet_tailwater
    223190                    log_to_file(log_filename, s)
    224191                Q = min(Q, Q_outlet_tailwater)
    … …  
    227194                #FIXME(Ole): What about inlet control?
    228195
     196               
    229197        # Common code for circle and square geometries
    230         if hasattr(culvert, 'log_filename'):                                   
    231             log_to_file(log_filename, 'Case: "%s"' %case)
     198        if log_filename is not None:
     199            log_to_file(log_filename, 'Case: "%s"' % case)
    232200           
    233         flow_rate_control=Q
    234 
    235         if hasattr(culvert, 'log_filename'):                                   
    236             s = 'Flow Rate Control = %f' %flow_rate_control
     201        if log_filename is not None:                       
     202            s = 'Flow Rate Control = %f' % Q
    237203            log_to_file(log_filename, s)
    238204
    239         inlet.rate = -flow_rate_control
    240         outlet.rate = flow_rate_control               
    241205           
    242         culv_froude=sqrt(flow_rate_control**2*width/(g*flow_area**3))
    243        
    244         if hasattr(culvert, 'log_filename'):                           
    245             s = 'Froude in Culvert = %f' %culv_froude
     206        culv_froude=sqrt(Q**2*width/(g*flow_area**3))
     207       
     208        if log_filename is not None:                           
     209            s = 'Froude in Culvert = %f' % culv_froude
    246210            log_to_file(log_filename, s)
    247211
    … …  
    250214
    251215
    252     else: #inlet.depth < 0.01:
     216    else: # inlet_depth < 0.01:
    253217        Q = barrel_velocity = outlet_culvert_depth = 0.0
    254218       

  • branches/numpy/anuga/fit_interpolate/fit.py

    r6533 r6553  
    280280           
    281281            element_found, sigma0, sigma1, sigma2, k = \
    282                            search_tree_of_vertices(self.root, self.mesh, x)
     282                           search_tree_of_vertices(self.root, x)
    283283           
    284284            if element_found is True:
    … …  
    301301                        AtA[j,k] += sigmas[j]*sigmas[k]
    302302            else:
    303                 # FIXME(Ole): This is the message referred to in ticket:314
    304                
    305303                flag = is_inside_polygon(x,
    306304                                         self.mesh_boundary_polygon,
    … …  
    310308                assert flag is True, msg               
    311309               
     310                # FIXME(Ole): This is the message referred to in ticket:314
    312311                minx = min(self.mesh_boundary_polygon[:,0])
    313312                maxx = max(self.mesh_boundary_polygon[:,0])               
    … …  
    317316                msg += 'Mesh boundary extent is (%.f, %.f), (%.f, %.f)'\
    318317                    % (minx, maxx, miny, maxy)
    319                
    320 
    321                 raise Exception(msg)
     318                raise RuntimeError, msg
     319
    322320            self.AtA = AtA
    323321
    … …  
    375373                self.build_fit_subset(points, z, verbose=verbose)
    376374
     375                # FIXME(Ole): I thought this test would make sense here
     376                # See test_fitting_example_that_crashed_2 in test_shallow_water_domain.py
     377                # Committed 11 March 2009
     378                msg = 'Matrix AtA was not built'
     379                assert self.AtA is not None, msg
     380               
     381                #print 'Matrix was built OK'
     382
    377383               
    378384            point_coordinates = None
    … …  
    382388        if point_coordinates is None:
    383389            if verbose: print 'Warning: no data points in fit'
    384             #assert self.AtA != None, 'no interpolation matrix'
    385             #assert self.Atz != None
    386             assert not self.AtA is None, 'no interpolation matrix'
    387             assert not self.Atz is None
     390            msg = 'No interpolation matrix.'
     391            assert self.AtA is not None, msg
     392            assert self.Atz is not None
    388393           
    389394            # FIXME (DSG) - do  a message
    … …  
    471476                alpha=DEFAULT_ALPHA,
    472477                verbose=False,
    473                 acceptable_overshoot=1.01,
    474                 # FIXME: Move to config - this value is assumed in caching test
    475                 # FIXME(Ole): Just realised that this was never implemented (29 Jan 2009). I suggest removing it altogether.
    476478                mesh_origin=None,
    477479                data_origin=None,
    … …  
    490492              'alpha': alpha,
    491493              'verbose': verbose,
    492               'acceptable_overshoot': acceptable_overshoot,
    493494              'mesh_origin': mesh_origin,
    494495              'data_origin': data_origin,
    … …  
    546547                 alpha=DEFAULT_ALPHA,
    547548                 verbose=False,
    548                  acceptable_overshoot=1.01,
    549549                 mesh_origin=None,
    550550                 data_origin=None,
    … …  
    572572          alpha: Smoothing parameter.
    573573
    574           acceptable overshoot: NOT IMPLEMENTED
    575           controls the allowed factor by which
    576           fitted values
    577           may exceed the value of input data. The lower limit is defined
    578           as min(z) - acceptable_overshoot*delta z and upper limit
    579           as max(z) + acceptable_overshoot*delta z
    580          
    581 
    582574          mesh_origin: A geo_reference object or 3-tuples consisting of
    583575              UTM zone, easting and northing.
    584576              If specified vertex coordinates are assumed to be
    585577              relative to their respective origins.
    586          
    587578
    588579          point_attributes: Vector or array of data at the

  • branches/numpy/anuga/fit_interpolate/interpolate.py

    r6428 r6553  
    484484            x = point_coordinates[i]
    485485            element_found, sigma0, sigma1, sigma2, k = \
    486                            search_tree_of_vertices(self.root, self.mesh, x)
    487 
    488             # Update interpolation matrix A if necessary
     486                           search_tree_of_vertices(self.root, x)
     487           
     488            # Update interpolation matrix A if necessary
    489489            if element_found is True:
    490490                # Assign values to matrix A

  • branches/numpy/anuga/fit_interpolate/search_functions.py

    r6304 r6553  
    88import time
    99
     10from anuga.utilities.polygon import is_inside_triangle
    1011from anuga.utilities.numerical_tools import get_machine_precision
    1112from anuga.config import max_float
    … …  
    1819search_more_cells_time = initial_search_value
    1920
    20 #FIXME test what happens if a
    21 LAST_TRIANGLE = [[-10,[(num.array([max_float,max_float]),
    22                         num.array([max_float,max_float]),
    23                         num.array([max_float,max_float])),
    24                        (num.array([1,1]),
    25                         num.array([0,0]),
    26                         num.array([-1.1,-1.1]))]]]
     21# FIXME(Ole): Could we come up with a less confusing structure?
     22LAST_TRIANGLE = [[-10,
     23                   (num.array([[max_float, max_float],
     24                               [max_float, max_float],
     25                               [max_float, max_float]]),
     26                    (num.array([1.,1.]),     
     27                     num.array([0.,0.]),     
     28                     num.array([-1.1,-1.1])))]]
    2729
    28 def search_tree_of_vertices(root, mesh, x):
     30def search_tree_of_vertices(root, x):
    2931    """
    3032    Find the triangle (element) that the point x is in.
    … …  
    3234    Inputs:
    3335        root: A quad tree of the vertices
    34         mesh: The underlying mesh
    3536        x:    The point being placed
    3637   
    … …  
    4748    global search_more_cells_time
    4849
    49     #Find triangle containing x:
    50     element_found = False
    51 
    52     # This will be returned if element_found = False
    53     sigma2 = -10.0
    54     sigma0 = -10.0
    55     sigma1 = -10.0
    56     k = -10.0
    57 
    5850    # Search the last triangle first
    59     try:
    60         element_found, sigma0, sigma1, sigma2, k = \
    61             _search_triangles_of_vertices(mesh, last_triangle, x)
    62     except:
    63         element_found = False
     51    element_found, sigma0, sigma1, sigma2, k = \
     52        _search_triangles_of_vertices(last_triangle, x)
    6453                   
    65     #print "last_triangle", last_triangle
    6654    if element_found is True:
    67         #print "last_triangle", last_triangle
    6855        return element_found, sigma0, sigma1, sigma2, k
    6956
    70     # This was only slightly faster than just checking the
    71     # last triangle and it significantly slowed down
    72     # non-gridded fitting
    73     # If the last element was a dud, search its neighbours
    74     #print "last_triangle[0][0]", last_triangle[0][0]
    75     #neighbours = mesh.get_triangle_neighbours(last_triangle[0][0])
    76     #print "neighbours", neighbours
    77     #neighbours = []
    78   #   for k in neighbours:
    79 #         if k >= 0:
    80 #             tri = mesh.get_vertex_coordinates(k,
    81 #                                                    absolute=True)
    82 #             n0 = mesh.get_normal(k, 0)
    83 #             n1 = mesh.get_normal(k, 1)
    84 #             n2 = mesh.get_normal(k, 2)
    85 #             triangle =[[k,(tri, (n0, n1, n2))]]
    86 #             element_found, sigma0, sigma1, sigma2, k = \
    87 #                            _search_triangles_of_vertices(mesh,
    88 #                                                          triangle, x)
    89 #             if element_found is True:
    90 #                 return element_found, sigma0, sigma1, sigma2, k
    91            
    92     #t0 = time.time()
     57   
    9358    # Get triangles in the cell that the point is in.
    9459    # Triangle is a list, first element triangle_id,
    9560    # second element the triangle
    9661    triangles = root.search(x[0], x[1])
     62    element_found, sigma0, sigma1, sigma2, k = \
     63                   _search_triangles_of_vertices(triangles, x)
     64
    9765    is_more_elements = True
    9866   
    99     element_found, sigma0, sigma1, sigma2, k = \
    100                    _search_triangles_of_vertices(mesh,
    101                                                  triangles, x)
    102     #search_one_cell_time += time.time()-t0
    103     #print "search_one_cell_time",search_one_cell_time
    104     #t0 = time.time()
    10567    while not element_found and is_more_elements:
    10668        triangles, branch = root.expand_search()
    … …  
    10971            # been searched.  This is the last try
    11072            element_found, sigma0, sigma1, sigma2, k = \
    111                            _search_triangles_of_vertices(mesh, triangles, x)
     73                           _search_triangles_of_vertices(triangles, x)
    11274            is_more_elements = False
    11375        else:
    11476            element_found, sigma0, sigma1, sigma2, k = \
    115                        _search_triangles_of_vertices(mesh, triangles, x)
    116         #search_more_cells_time += time.time()-t0
    117     #print "search_more_cells_time", search_more_cells_time
     77                       _search_triangles_of_vertices(triangles, x)
     78                       
    11879       
    11980    return element_found, sigma0, sigma1, sigma2, k
    12081
    12182
    122 def _search_triangles_of_vertices(mesh, triangles, x):
    123     """Search for triangle containing x amongs candidate_vertices in mesh
     83def _search_triangles_of_vertices(triangles, x):
     84    """Search for triangle containing x amongs candidate_vertices in triangles
    12485
    12586    This is called by search_tree_of_vertices once the appropriate node
    … …  
    13293    global last_triangle
    13394   
    134     # these statments are needed if triangles is empty
    135     #Find triangle containing x:
    136     element_found = False
    137 
    138     # This will be returned if element_found = False
     95    # These statments are needed if triangles is empty
    13996    sigma2 = -10.0
    14097    sigma0 = -10.0
    14198    sigma1 = -10.0
    14299    k = -10
    143    
    144     #For all vertices in same cell as point x
     100    # For all vertices in same cell as point x
     101    element_found = False   
    145102    for k, tri_verts_norms in triangles:
    146103        tri = tri_verts_norms[0]
    147         n0, n1, n2 = tri_verts_norms[1]
    148104        # k is the triangle index
    149105        # tri is a list of verts (x, y), representing a tringle
    150106        # Find triangle that contains x (if any) and interpolate
    151         element_found, sigma0, sigma1, sigma2 =\
    152                        find_triangle_compute_interpolation(tri, n0, n1, n2, x)
    153         if element_found is True:
     107       
     108        # Input check disabled to speed things up.
     109        if is_inside_triangle(x, tri,
     110                              closed=True,
     111                              check_inputs=False):
     112           
     113            n0, n1, n2 = tri_verts_norms[1]       
     114            sigma0, sigma1, sigma2 =\
     115                compute_interpolation_values(tri, n0, n1, n2, x)
     116               
     117            element_found = True   
     118
    154119            # Don't look for any other triangles in the triangle list
    155             last_triangle = [[k,tri_verts_norms]]
    156             break
     120            last_triangle = [[k, tri_verts_norms]]
     121            break           
     122           
    157123    return element_found, sigma0, sigma1, sigma2, k
    158124
    159125
    160126           
    161 def find_triangle_compute_interpolation(triangle, n0, n1, n2, x):
    162     """Compute linear interpolation of point x and triangle k in mesh.
    163     It is assumed that x belongs to triangle k.max_float
     127def compute_interpolation_values(triangle, n0, n1, n2, x):
     128    """Compute linear interpolation of point x and triangle.
     129   
     130    n0, n1, n2 are normal to the tree edges.
    164131    """
    165132
    … …  
    167134    xi0, xi1, xi2 = triangle
    168135
    169     # this is where we can call some fast c code.
    170      
    171     # Integrity check - machine precision is too hard
    172     # so we use hardwired single precision
    173     epsilon = 1.0e-6
    174    
    175     if  x[0] > max(xi0[0], xi1[0], xi2[0]) + epsilon:
    176         # print "max(xi0[0], xi1[0], xi2[0])", max(xi0[0], xi1[0], xi2[0])
    177         return False,0,0,0
    178     if  x[0] < min(xi0[0], xi1[0], xi2[0]) - epsilon:
    179         return False,0,0,0
    180     if  x[1] > max(xi0[1], xi1[1], xi2[1]) + epsilon:
    181         return False,0,0,0
    182     if  x[1] < min(xi0[1], xi1[1], xi2[1]) - epsilon:
    183         return False,0,0,0
    184    
    185     # machine precision on some machines (e.g. nautilus)
    186     epsilon = get_machine_precision() * 2
    187    
    188     # Compute interpolation - return as soon as possible
    189     #  print "(xi0-xi1)", (xi0-xi1)
    190     # print "n0", n0
    191     # print "dot((xi0-xi1), n0)", dot((xi0-xi1), n0)
    192    
    193136    sigma0 = num.dot((x-xi1), n0)/num.dot((xi0-xi1), n0)
    194     if sigma0 < -epsilon:
    195         return False,0,0,0
    196137    sigma1 = num.dot((x-xi2), n1)/num.dot((xi1-xi2), n1)
    197     if sigma1 < -epsilon:
    198         return False,0,0,0
    199138    sigma2 = num.dot((x-xi0), n2)/num.dot((xi2-xi0), n2)
    200     if sigma2 < -epsilon:
    201         return False,0,0,0
    202    
    203     # epsilon = 1.0e-6
    204     # we want to speed this up, so don't do assertions
    205     #delta = abs(sigma0+sigma1+sigma2-1.0) # Should be close to zero
    206     #msg = 'abs(sigma0+sigma1+sigma2-1) = %.15e, eps = %.15e'\
    207     #      %(delta, epsilon)
    208     #assert delta < epsilon, msg
    209139
    210 
    211     # Check that this triangle contains the data point
    212     # Sigmas are allowed to get negative within
    213     # machine precision on some machines (e.g. nautilus)
    214     #if sigma0 >= -epsilon and sigma1 >= -epsilon and sigma2 >= -epsilon:
    215     #    element_found = True
    216     #else:
    217     #    element_found = False
    218     return True, sigma0, sigma1, sigma2
     140    return sigma0, sigma1, sigma2
    219141
    220142def set_last_triangle():

  • branches/numpy/anuga/fit_interpolate/test_search_functions.py

    r6360 r6553  
    55from search_functions import search_tree_of_vertices, set_last_triangle
    66from search_functions import _search_triangles_of_vertices
     7from search_functions import compute_interpolation_values
    78
    89from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
    … …  
    1011from anuga.utilities.polygon import is_inside_polygon
    1112from anuga.utilities.quad import build_quadtree, Cell
    12 
     13from anuga.utilities.numerical_tools import ensure_numeric
     14from anuga.utilities.polygon import is_inside_polygon, is_inside_triangle   
     15
     16import numpy as num
    1317
    1418class Test_search_functions(unittest.TestCase):
    … …  
    4953
    5054        x = [0.2, 0.7]
    51         found, s0, s1, s2, k = search_tree_of_vertices(root, mesh, x)
     55        found, s0, s1, s2, k = search_tree_of_vertices(root, x)
    5256        assert k == 1 # Triangle one
    5357        assert found is True
    5458
    5559    def test_bigger(self):
    56         """test_larger mesh
     60        """test_bigger
     61       
     62        test larger mesh
    5763        """
    5864
    … …  
    6975                  [0.1,0.9], [0.4,0.6], [0.9,0.1],
    7076                  [10, 3]]:
    71                
    72             found, s0, s1, s2, k = search_tree_of_vertices(root, mesh, x)
     77           
     78            found, s0, s1, s2, k = search_tree_of_vertices(root,
     79                                                           ensure_numeric(x))
    7380
    7481            if k >= 0:
    … …  
    102109                      [10, 3]]:
    103110               
    104                 found, s0, s1, s2, k = search_tree_of_vertices(root, mesh, x)
     111                found, s0, s1, s2, k = search_tree_of_vertices(root, x)
    105112
    106113                if k >= 0:
    … …  
    111118                    assert found is False               
    112119
    113                
     120           
     121                if m == k == 0: return   
    114122
    115123    def test_underlying_function(self):
    … …  
    124132
    125133        # One point
    126         x = [0.5, 0.5]
     134        x = ensure_numeric([0.5, 0.5])
    127135        candidate_vertices = root.search(x[0], x[1])
    128136
    129137        #print x, candidate_vertices
    130138        found, sigma0, sigma1, sigma2, k = \
    131                _search_triangles_of_vertices(mesh,
    132                                              candidate_vertices,
     139               _search_triangles_of_vertices(candidate_vertices,
    133140                                             x)
    134141
    … …  
    151158            #print x, candidate_vertices
    152159            found, sigma0, sigma1, sigma2, k = \
    153                    _search_triangles_of_vertices(mesh,
    154                                                  candidate_vertices,
    155                                                  x)
     160                   _search_triangles_of_vertices(candidate_vertices,
     161                                                 ensure_numeric(x))
    156162            if k >= 0:
    157163                V = mesh.get_vertex_coordinates(k) # nodes for triangle k
    … …  
    199205        x = [2.5, 1.5]
    200206        element_found, sigma0, sigma1, sigma2, k = \
    201                        search_tree_of_vertices(root, mesh, x)
     207                       search_tree_of_vertices(root, x)
    202208        # One point
    203209        x = [3.00005, 2.999994]
    204210        element_found, sigma0, sigma1, sigma2, k = \
    205                        search_tree_of_vertices(root, mesh, x)
     211                       search_tree_of_vertices(root, x)
    206212        assert element_found is True
    207213        assert k == 1
    208214       
    209 ################################################################################
    210 
     215       
     216    def test_compute_interpolation_values(self):
     217        """test_compute_interpolation_values
     218       
     219        Test that interpolation values are correc
     220       
     221        This test used to check element_found as output from
     222        find_triangle_compute_interpolation(triangle, n0, n1, n2, x)
     223        and that failed before 18th March 2009.
     224       
     225        Now this function no longer returns this flag, so the test
     226        is merely checknig the sigmas.
     227       
     228        """
     229       
     230        triangle = num.array([[306951.77151059, 6194462.14596986],
     231                              [306952.58403545, 6194459.65001246],
     232                              [306953.55109034, 6194462.0041216]])
     233
     234
     235        n0 = [0.92499377, -0.37998227]
     236        n1 = [0.07945684,  0.99683831]
     237        n2 = [-0.95088404, -0.30954732]
     238       
     239        x = [306953.344, 6194461.5]
     240       
     241        # Test that point is indeed inside triangle
     242        assert is_inside_polygon(x, triangle,
     243                                 closed=True, verbose=False)
     244        assert is_inside_triangle(x, triangle,
     245                                  closed=True, verbose=False)                                 
     246       
     247        sigma0, sigma1, sigma2 = \
     248            compute_interpolation_values(triangle, n0, n1, n2, x)
     249           
     250        msg = 'Point which is clearly inside triangle was not found'
     251        assert abs(sigma0 + sigma1 + sigma2 - 1) < 1.0e-6, msg
     252
     253#-------------------------------------------------------------
    211254if __name__ == "__main__":
    212     suite = unittest.makeSuite(Test_search_functions,'test')
    213     runner = unittest.TextTestRunner() #verbosity=1)
     255    suite = unittest.makeSuite(Test_search_functions, 'test')
     256    runner = unittest.TextTestRunner(verbosity=1)
    214257    runner.run(suite)
    215258   

  • branches/numpy/anuga/load_mesh/loadASCII.py

    r6428 r6553  
    821821        titles = fid.variables['vertex_attribute_titles'][:]
    822822        mesh['vertex_attribute_titles'] = [x.tostring().strip() for x in titles]
    823     except:
     823    except KeyError:
    824824        pass
    825825
    … …  
    833833        tags = fid.variables['segment_tags'][:]
    834834        mesh['segment_tags'] = [x.tostring().strip() for x in tags]
    835     except:
    836         pass
     835    except KeyError:
     836        for ob in mesh['segments']:
     837            mesh['segment_tags'].append('')
    837838
    838839    try:

  • branches/numpy/anuga/load_mesh/test_loadASCII.py

    r6428 r6553  
    291291            pass
    292292        else:
    293             self.failUnless(0 == 1, 'bad tsh file did not raise error!')
     293            self.fail('bad tsh file did not raise error!')
    294294        os.remove(fileName)
    295295
    … …  
    307307            pass
    308308        else:
    309             self.failUnless(0 == 1, 'bad tsh file did not raise error!')
     309            self.fail('bad tsh file did not raise error!')
    310310        os.remove(fileName)
    311311
    … …  
    456456            pass
    457457        else:
    458             self.failUnless(0 == 1, 'imaginary file did not raise error!')
     458            self.fail('imaginary file did not raise error!')
    459459
    460460    def throws_error_2_screen_test_import_mesh_bad(self):
    … …  
    472472            pass
    473473        else:
    474             self.failUnless(0 == 1, 'bad msh file did not raise error!')
     474            self.fail('bad msh file did not raise error!')
    475475        os.remove(fileName)
    476476

  • branches/numpy/anuga/shallow_water/__init__.py

    r6533 r6553  
    1212     Transmissive_momentum_set_stage_boundary,\
    1313     Dirichlet_discharge_boundary,\
    14      Field_boundary
     14     Field_boundary,\
     15     Transmissive_stage_zero_momentum_boundary
    1516
    1617

  • branches/numpy/anuga/shallow_water/data_manager.py

    r6533 r6553  
    32363236            volumes.append([v4,v3,v2]) #Lower element
    32373237
    3238     volumes = num.array(volumes)
     3238    volumes = num.array(volumes, num.int)      #array default#
    32393239
    32403240    if origin is None:
    … …  
    41524152            volumes.append([v4,v2,v3]) #Lower element
    41534153
    4154     volumes = num.array(volumes)
     4154    volumes = num.array(volumes, num.int)      #array default#
    41554155
    41564156    geo_ref = Geo_reference(refzone, min(x), min(y))
    … …  
    55715571    for i in range(len(quantities)):
    55725572        for file_in in files_in[i]:
    5573             if (os.access(file_in, os.F_OK) == 0):
     5573            if (os.access(file_in, os.R_OK) == 0):
    55745574                msg = 'File %s does not exist or is not accessible' % file_in
    55755575                raise IOError, msg

  • branches/numpy/anuga/shallow_water/shallow_water_domain.py

    r6451 r6553  
    574574    # @brief Run integrity checks on shallow water domain.
    575575    def check_integrity(self):
    576         Generic_Domain.check_integrity(self)    # super integrity check
     576        Generic_Domain.check_integrity(self)
    577577
    578578        #Check that we are solving the shallow water wave equation
    … …  
    21502150            default_rain = None
    21512151
    2152         # pass to parent class
     2152           
     2153           
    21532154        General_forcing.__init__(self,
    21542155                                 domain,

  • branches/numpy/anuga/shallow_water/test_data_manager.py

    r6472 r6553  
    12661266        time = fid.variables['time'][:]
    12671267        stage = fid.variables['stage'][:]
     1268       
     1269        # Check georeferencig: zone, xllcorner and yllcorner
     1270        assert fid.zone[0] == 56
     1271        assert fid.xllcorner[0] == 308500
     1272        assert fid.yllcorner[0] == 6189000
     1273               
    12681274
    12691275        fid.close()

  • branches/numpy/anuga/shallow_water/test_shallow_water_domain.py

    r6533 r6553  
    16721672        # Setup computational domain
    16731673        #-----------------------------------------------------------------
    1674         points, vertices, boundary = rectangular_cross(10, 10)    # Basic mesh
    1675         domain = Domain(points, vertices, boundary)    # Create domain
     1674        points, vertices, boundary = rectangular_cross(10, 10) # Basic mesh
     1675        domain = Domain(points, vertices, boundary) # Create domain
     1676        domain.set_default_order(1)       
    16761677        domain.set_quantities_to_be_stored(None)
    16771678        domain.set_maximum_allowed_speed(100) #FIXME (Ole): try to remove this
    … …  
    18021803
    18031804        assert num.allclose(gauge_values[0], G0)
    1804         assert num.allclose(gauge_values[1], G1)
     1805        # FIXME(Ole): Disabled when ticket:314 was resolved.
     1806        # The slight change might in fact be for the better.
     1807        #assert num.allclose(gauge_values[1], G1)
    18051808        assert num.allclose(gauge_values[2], G2)
    18061809        assert num.allclose(gauge_values[3], G3)
    … …  
    30233026        a = [0.0, 0.0]
    30243027        b = [0.0, 2.0]
    3025         c = [2.0, 0.0]
     3028        c = [2.0,0.0]
    30263029        d = [0.0, 4.0]
    30273030        e = [2.0, 2.0]
    3028         f = [4.0, 0.0]
     3031        f = [4.0,0.0]
    30293032
    30303033        points = [a, b, c, d, e, f]
    3031         #             bac,     bce,     ecf,    dbe
    3032         vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
     3034        #bac, bce, ecf, dbe
     3035        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    30333036
    30343037        domain = Domain(points, vertices)
    3035         val0 = 2. + 2.0/3
    3036         val1 = 4. + 4.0/3
    3037         val2 = 8. + 2.0/3
    3038         val3 = 2. + 8.0/3
     3038        val0 = 2.+2.0/3
     3039        val1 = 4.+4.0/3
     3040        val2 = 8.+2.0/3
     3041        val3 = 2.+8.0/3
    30393042
    30403043        domain.set_quantity('elevation', [[0,0,0], [6,0,0],
    … …  
    30483051        L = domain.quantities['stage'].vertex_values
    30493052
    3050         # Check that some stages are not above elevation (within eps) -
    3051         # so that the limiter has something to work with
    3052         assert not num.alltrue(num.alltrue(num.greater_equal(L, E-epsilon)))
     3053
     3054        #Check that some stages are not above elevation (within eps)
     3055        #- so that the limiter has something to work with
     3056        assert not num.alltrue(num.alltrue(num.greater_equal(L,E-epsilon)))
    30533057
    30543058        domain._order_ = 1
    … …  
    30563060
    30573061        #Check that all stages are above elevation (within eps)
    3058         assert num.alltrue(num.alltrue(num.greater_equal(L, E-epsilon)))
    3059 
    3060     #####################################################
     3062        assert num.alltrue(num.alltrue(num.greater_equal(L,E-epsilon)))
    30613063
    30623064    def test_distribute_basic(self):

  • branches/numpy/anuga/shallow_water/urs_ext.c

    r6304 r6553  
    99#include "math.h"
    1010#include <stdio.h>
     11#include <string.h>
     12#include <errno.h>
    1113#include <float.h>
    1214#include <time.h>
    … …  
    219221        if((fp = fopen(muxFileName, "r")) == NULL)
    220222        {
    221             fprintf(stderr, "cannot open file %s\n", muxFileName);
     223            char *err_msg = strerror(errno);
     224
     225            fprintf(stderr, "cannot open file '%s': %s\n", muxFileName, err_msg);
    222226            return -1; 
    223227        }

  • branches/numpy/anuga/test_all.py

    r6533 r6553  
    215215        #os.remove(filename)
    216216
     217   
     218    if sys.platform == 'win32':
     219        raw_input('Press any key')

  • branches/numpy/anuga/utilities/polygon.py

    r6410 r6553  
    211211    return status, value
    212212
    213 ##
    214 # @brief Determine if one point is inside a polygon.
    215 # @param point The point of interest.
    216 # @param polygon The polygon to test inclusion in.
    217 # @param closed True if points on boundary are considered 'inside' polygon.
    218 # @param verbose True if this function is to be verbose.
    219 # @return True if point is inside the polygon.
    220 # @note Uses inside_polygon() to do the work.
    221 # @note Raises Exception if more than one point supplied.
     213def is_inside_triangle(point, triangle,
     214                       closed=True,
     215                       rtol=1.0e-12,
     216                       atol=1.0e-12,                     
     217                       check_inputs=True,
     218                       verbose=False):
     219    """Determine if one point is inside a triangle
     220   
     221    This uses the barycentric method:
     222   
     223    Triangle is A, B, C
     224    Point P can then be written as
     225   
     226    P = A + alpha * (C-A) + beta * (B-A)
     227    or if we let
     228    v=P-A, v0=C-A, v1=B-A   
     229   
     230    v = alpha*v0 + beta*v1
     231
     232    Dot this equation by v0 and v1 to get two:
     233   
     234    dot(v0, v) = alpha*dot(v0, v0) + beta*dot(v0, v1)
     235    dot(v1, v) = alpha*dot(v1, v0) + beta*dot(v1, v1)   
     236   
     237    or if a_ij = dot(v_i, v_j) and b_i = dot(v_i, v)
     238    the matrix equation:
     239   
     240    a_00 a_01   alpha     b_0
     241                       =
     242    a_10 a_11   beta      b_1
     243   
     244    Solving for alpha and beta yields:
     245   
     246    alpha = (b_0*a_11 - b_1*a_01)/denom
     247    beta =  (b_1*a_00 - b_0*a_10)/denom
     248   
     249    with denom = a_11*a_00 - a_10*a_01
     250   
     251    The point is in the triangle whenever
     252    alpha and beta and their sums are in the unit interval.
     253   
     254    rtol and atol will determine how close the point has to be to the edge
     255    before it is deemed to be on the edge.
     256   
     257    """
     258
     259    triangle = ensure_numeric(triangle)       
     260    point = ensure_numeric(point, num.float)   
     261   
     262    if check_inputs is True:
     263        msg = 'is_inside_triangle must be invoked with one point only'
     264        assert num.allclose(point.shape, [2]), msg
     265   
     266   
     267    # Use C-implementation
     268    return bool(_is_inside_triangle(point, triangle, int(closed), rtol, atol))
     269   
     270   
     271
     272    # FIXME (Ole): The rest of this function has been made
     273    # obsolete by the C extension.
     274   
     275    # Quickly reject points that are clearly outside
     276    if point[0] < min(triangle[:,0]): return False
     277    if point[0] > max(triangle[:,0]): return False   
     278    if point[1] < min(triangle[:,1]): return False
     279    if point[1] > max(triangle[:,1]): return False       
     280
     281
     282    # Start search   
     283    A = triangle[0, :]
     284    B = triangle[1, :]
     285    C = triangle[2, :]
     286   
     287    # Now check if point lies wholly inside triangle
     288    v0 = C-A
     289    v1 = B-A       
     290       
     291    a00 = num.inner(v0, v0)
     292    a10 = a01 = num.inner(v0, v1)
     293    a11 = num.inner(v1, v1)
     294   
     295    denom = a11*a00 - a01*a10
     296   
     297    if abs(denom) > 0.0:
     298        v = point-A
     299        b0 = num.inner(v0, v)       
     300        b1 = num.inner(v1, v)           
     301   
     302        alpha = (b0*a11 - b1*a01)/denom
     303        beta = (b1*a00 - b0*a10)/denom       
     304
     305        if (alpha > 0.0) and (beta > 0.0) and (alpha+beta < 1.0):
     306            return True
     307
     308
     309    if closed is True:
     310        # Check if point lies on one of the edges
     311       
     312        for X, Y in [[A,B], [B,C], [C,A]]:
     313            res = _point_on_line(point[0], point[1],
     314                                 X[0], X[1],
     315                                 Y[0], Y[1],
     316                                 rtol, atol)
     317           
     318            if res:
     319                return True
     320               
     321    return False
     322
     323def is_inside_polygon_quick(point, polygon, closed=True, verbose=False):
     324    """Determine if one point is inside a polygon
     325    Both point and polygon are assumed to be numeric arrays or lists and
     326    no georeferencing etc or other checks will take place.
     327   
     328    As such it is faster than is_inside_polygon
     329    """
     330
     331    # FIXME(Ole): This function isn't being used
     332    polygon = ensure_numeric(polygon, num.float)
     333    points = ensure_numeric(point, num.float) # Convert point to array of points
     334    points = num.ascontiguousarray(points[num.newaxis, :])
     335    msg = ('is_inside_polygon() must be invoked with one point only.\n'
     336           'I got %s and converted to %s' % (str(point), str(points.shape)))
     337    assert points.shape[0] == 1 and points.shape[1] == 2, msg
     338   
     339    indices = num.zeros(1, num.int)
     340
     341    count = _separate_points_by_polygon(points, polygon, indices,
     342                                        int(closed), int(verbose))
     343
     344    return count > 0
     345
     346
    222347def is_inside_polygon(point, polygon, closed=True, verbose=False):
    223348    """Determine if one point is inside a polygon
    … …  
    234359    else:
    235360        msg = 'is_inside_polygon must be invoked with one point only'
    236         raise Exception, msg
     361        raise msg
    237362
    238363##
    … …  
    263388        # If this fails it is going to be because the points can't be
    264389        # converted to a numeric array.
    265         msg = 'Points could not be converted to numeric array'
     390        msg = 'Points could not be converted to Numeric array'
    266391        raise Exception, msg
    267392
     393    polygon = ensure_absolute(polygon)       
    268394    try:
    269395        polygon = ensure_absolute(polygon)
    … …  
    416542# @return (indices, count) where indices are point indices and count is the
    417543#          delimiter index between point inside (on left) and others.
    418 def separate_points_by_polygon(points, polygon, closed=True, verbose=False):
     544def separate_points_by_polygon(points, polygon,
     545                               closed=True,
     546                               check_input=True,
     547                               verbose=False):
    419548    """Determine whether points are inside or outside a polygon
    420549
    … …  
    425554       regarded as belonging to the polygon (closed = True)
    426555       or not (closed = False)
     556       check_input: Allows faster execution if set to False
    427557
    428558    Outputs:
    … …  
    456586    """
    457587
    458     #Input checks
    459     assert isinstance(closed, bool), 'Keyword arg "closed" must be boolean'
    460     assert isinstance(verbose, bool), 'Keyword arg "verbose" must be boolean'
    461 
    462     try:
    463         points = ensure_numeric(points, num.float)
    464     except NameError, e:
    465         raise NameError, e
    466     except:
    467         msg = 'Points could not be converted to numeric array'
    468         raise Exception, msg
    469 
    470     try:
    471         polygon = ensure_numeric(polygon, num.float)
    472     except NameError, e:
    473         raise NameError, e
    474     except:
    475         msg = 'Polygon could not be converted to numeric array'
    476         raise Exception, msg
    477 
    478     msg = 'Polygon array must be a 2d array of vertices'
    479     assert len(polygon.shape) == 2, msg
    480 
    481     msg = 'Polygon array must have two columns'
    482     assert polygon.shape[1] == 2, msg
    483 
    484     msg = ('Points array must be 1 or 2 dimensional. I got %d dimensions'
    485            % len(points.shape))
    486     assert 0 < len(points.shape) < 3, msg
    487 
    488     if len(points.shape) == 1:
    489         # Only one point was passed in.
    490         # Convert to array of points
    491         points = num.reshape(points, (1,2))
    492 
    493     msg = ('Point array must have two columns (x,y). I got points.shape[1]=%d'
    494            % points.shape[1])
    495     assert points.shape[1] == 2, msg
    496 
    497     msg = 'Points array must be a 2d array. I got %s' % str(points[:30])
    498     assert len(points.shape) == 2, msg
    499 
    500     msg = 'Points array must have two columns'
    501     assert points.shape[1] == 2, msg
    502 
    503     N = polygon.shape[0]    # Number of vertices in polygon
    504     M = points.shape[0]     # Number of points
    505 
    506     indices = num.zeros( M, num.int )
     588    if check_input:
     589        #Input checks
     590        assert isinstance(closed, bool), 'Keyword argument "closed" must be boolean'
     591        assert isinstance(verbose, bool), 'Keyword argument "verbose" must be boolean'
     592
     593        try:
     594            points = ensure_numeric(points, num.float)
     595        except NameError, e:
     596            raise NameError, e
     597        except:
     598            msg = 'Points could not be converted to numeric array'
     599            raise msg
     600
     601        try:
     602            polygon = ensure_numeric(polygon, num.float)
     603        except NameError, e:
     604            raise NameError, e
     605        except:
     606            msg = 'Polygon could not be converted to numeric array'
     607            raise msg
     608
     609        msg = 'Polygon array must be a 2d array of vertices'
     610        assert len(polygon.shape) == 2, msg
     611
     612        msg = 'Polygon array must have two columns'
     613        assert polygon.shape[1]==2, msg
     614
     615        msg = ('Points array must be 1 or 2 dimensional. '
     616               'I got %d dimensions' % len(points.shape))
     617        assert 0 < len(points.shape) < 3, msg
     618
     619        if len(points.shape) == 1:
     620            # Only one point was passed in.  Convert to array of points.
     621            points = num.reshape(points, (1,2))
     622   
     623            msg = ('Point array must have two columns (x,y), '
     624                   'I got points.shape[1]=%d' % points.shape[0])
     625            assert points.shape[1]==2, msg
     626
     627       
     628            msg = ('Points array must be a 2d array. I got %s.'
     629                   % str(points[:30]))
     630            assert len(points.shape)==2, msg
     631
     632            msg = 'Points array must have two columns'
     633            assert points.shape[1]==2, msg
     634
     635    N = polygon.shape[0] # Number of vertices in polygon
     636    M = points.shape[0]  # Number of points
     637
     638    indices = num.zeros(M, num.int)
    507639
    508640    count = _separate_points_by_polygon(points, polygon, indices,
    509641                                        int(closed), int(verbose))
    510642
    511     if verbose: print 'Found %d points (out of %d) inside polygon' % (count, M)
     643    if verbose:
     644        print 'Found %d points (out of %d) inside polygon' % (count, M)
    512645
    513646    return indices, count
    … …  
    11891322    from polygon_ext import _point_on_line
    11901323    from polygon_ext import _separate_points_by_polygon
    1191     from polygon_ext import _interpolate_polyline
     1324    from polygon_ext import _interpolate_polyline   
     1325    from polygon_ext import _is_inside_triangle       
    11921326    #from polygon_ext import _intersection
    11931327

  • branches/numpy/anuga/utilities/polygon_ext.c

    r6410 r6553  
    247247
    248248
     249int __is_inside_triangle(double* point,
     250                         double* triangle,
     251                         int closed,
     252                         double rtol,
     253                         double atol) {
     254                         
     255  double vx, vy, v0x, v0y, v1x, v1y;
     256  double a00, a10, a01, a11, b0, b1;
     257  double denom, alpha, beta;
     258 
     259  double x, y; // Point coordinates
     260  int i, j, res;
     261
     262  x = point[0];
     263  y = point[1];
     264 
     265  // Quickly reject points that are clearly outside
     266  if ((x < triangle[0]) &&
     267      (x < triangle[2]) &&
     268      (x < triangle[4])) return 0;       
     269     
     270  if ((x > triangle[0]) &&
     271      (x > triangle[2]) &&
     272      (x > triangle[4])) return 0;             
     273 
     274  if ((y < triangle[1]) &&
     275      (y < triangle[3]) &&
     276      (y < triangle[5])) return 0;       
     277     
     278  if ((y > triangle[1]) &&
     279      (y > triangle[3]) &&
     280      (y > triangle[5])) return 0;             
     281 
     282 
     283  // v0 = C-A
     284  v0x = triangle[4]-triangle[0];
     285  v0y = triangle[5]-triangle[1];
     286 
     287  // v1 = B-A   
     288  v1x = triangle[2]-triangle[0];
     289  v1y = triangle[3]-triangle[1];
     290
     291  // First check if point lies wholly inside triangle
     292  a00 = v0x*v0x + v0y*v0y; // innerproduct(v0, v0)
     293  a01 = v0x*v1x + v0y*v1y; // innerproduct(v0, v1)
     294  a10 = a01;               // innerproduct(v1, v0)
     295  a11 = v1x*v1x + v1y*v1y; // innerproduct(v1, v1)
     296   
     297  denom = a11*a00 - a01*a10;
     298
     299  if (fabs(denom) > 0.0) {
     300    // v = point-A 
     301    vx = x - triangle[0];
     302    vy = y - triangle[1];     
     303   
     304    b0 = v0x*vx + v0y*vy; // innerproduct(v0, v)       
     305    b1 = v1x*vx + v1y*vy; // innerproduct(v1, v)           
     306   
     307    alpha = (b0*a11 - b1*a01)/denom;
     308    beta = (b1*a00 - b0*a10)/denom;       
     309   
     310    if ((alpha > 0.0) && (beta > 0.0) && (alpha+beta < 1.0)) return 1;
     311  }
     312
     313  if (closed) {
     314    // Check if point lies on one of the edges
     315       
     316    for (i=0; i<3; i++) {
     317      j = (i+1) % 3; // Circular index into triangle array
     318      res = __point_on_line(x, y,
     319                            triangle[2*i], triangle[2*i+1],
     320                            triangle[2*j], triangle[2*j+1],                         
     321                            rtol, atol);
     322      if (res) return 1;
     323    }
     324  }
     325               
     326  // Default return if point is outside triangle                         
     327  return 0;                                             
     328}                                                                             
     329
     330
    249331int __separate_points_by_polygon(int M,     // Number of points
    250332                                int N,     // Number of polygon vertices
    … …  
    427509
    428510
     511     
     512PyObject *_is_inside_triangle(PyObject *self, PyObject *args) {
     513  //
     514  // _is_inside_triangle(point, triangle, int(closed), rtol, atol)
     515  //
     516
     517 
     518  PyArrayObject
     519    *point,
     520    *triangle;
     521
     522  double rtol, atol; 
     523  int closed, res;
     524
     525  PyObject *result;
     526     
     527  // Convert Python arguments to C
     528  if (!PyArg_ParseTuple(args, "OOidd",
     529                        &point,
     530                        &triangle,
     531                        &closed,
     532                        &rtol,
     533                        &atol)) {
     534   
     535    PyErr_SetString(PyExc_RuntimeError,
     536                    "_is_inside_triangle could not parse input");
     537    return NULL;
     538  }
     539
     540  // Call underlying routine
     541  res = __is_inside_triangle((double*) point -> data,
     542                             (double*) triangle -> data,
     543                             closed,
     544                             rtol,
     545                             atol);                                                   
     546
     547
     548  // Return result
     549  result = Py_BuildValue("i", res);
     550  return result; 
     551}
     552     
     553
     554
    429555/*
    430556PyObject *_intersection(PyObject *self, PyObject *args) {
    … …  
    554680  {"_interpolate_polyline", _interpolate_polyline,
    555681                                 METH_VARARGS, "Print out"},                             
     682  {"_is_inside_triangle", _is_inside_triangle,
     683                                 METH_VARARGS, "Print out"},
    556684  {NULL, NULL, 0, NULL}   /* sentinel */
    557685};

  • branches/numpy/anuga/utilities/test_polygon.py

    r6441 r6553  
    181181        polygon = [[0,0], [1,0], [1,1], [0,1]]
    182182
     183        assert is_inside_polygon_quick( (0.5, 0.5), polygon )
     184        assert not is_inside_polygon_quick( (0.5, 1.5), polygon )
     185        assert not is_inside_polygon_quick( (0.5, -0.5), polygon )
     186        assert not is_inside_polygon_quick( (-0.5, 0.5), polygon )
     187        assert not is_inside_polygon_quick( (1.5, 0.5), polygon )
     188
     189        # Try point on borders
     190        assert is_inside_polygon_quick( (1., 0.5), polygon, closed=True)
     191        assert is_inside_polygon_quick( (0.5, 1), polygon, closed=True)
     192        assert is_inside_polygon_quick( (0., 0.5), polygon, closed=True)
     193        assert is_inside_polygon_quick( (0.5, 0.), polygon, closed=True)
     194
     195        assert not is_inside_polygon_quick( (0.5, 1), polygon, closed=False)
     196        assert not is_inside_polygon_quick( (0., 0.5), polygon, closed=False)
     197        assert not is_inside_polygon_quick( (0.5, 0.), polygon, closed=False)
     198        assert not is_inside_polygon_quick( (1., 0.5), polygon, closed=False)
     199
     200
     201    def test_inside_polygon_main(self):
     202        # Simplest case: Polygon is the unit square
     203        polygon = [[0,0], [1,0], [1,1], [0,1]]
     204
    183205        # From real example (that failed)
    184206        polygon = [[20,20], [40,20], [40,40], [20,40]]
    … …  
    195217        # More convoluted and non convex polygon
    196218        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
    197         assert is_inside_polygon((0.5, 0.5), polygon)
    198         assert is_inside_polygon((1, -0.5), polygon)
    199         assert is_inside_polygon((1.5, 0), polygon)
    200 
    201         assert not is_inside_polygon((0.5, 1.5), polygon)
    202         assert not is_inside_polygon((0.5, -0.5), polygon)
     219        assert is_inside_polygon( (0.5, 0.5), polygon )
     220        assert is_inside_polygon( (1, -0.5), polygon )
     221        assert is_inside_polygon( (1.5, 0), polygon )
     222
     223        assert not is_inside_polygon( (0.5, 1.5), polygon )
     224        assert not is_inside_polygon( (0.5, -0.5), polygon )
     225
     226        assert is_inside_polygon_quick( (0.5, 0.5), polygon )
     227        assert is_inside_polygon_quick( (1, -0.5), polygon )
     228        assert is_inside_polygon_quick( (1.5, 0), polygon )
     229
     230        assert not is_inside_polygon_quick( (0.5, 1.5), polygon )
     231        assert not is_inside_polygon_quick( (0.5, -0.5), polygon )
    203232
    204233        # Very convoluted polygon
    … …  
    365394
    366395        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
    367         points = [[0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0],
    368                   [0.5, 1.5], [0.5, -0.5]]
    369         res, count = separate_points_by_polygon(points, polygon)
    370 
    371         assert num.allclose( res, [1, 2, 3, 5, 4, 0] )
     396        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
     397        res, count = separate_points_by_polygon( points, polygon )
     398
     399        assert num.allclose( res, [1,2,3,5,4,0] )
    372400        assert count == 3
     401
     402
     403    def test_is_inside_triangle(self):
     404        # Simplest case:
     405        triangle = [[0, 0], [1, 0], [0.5, 1]]
     406
     407        assert is_inside_triangle((0.5, 0.5), triangle)
     408        assert is_inside_triangle((0.9, 0.1), triangle)
     409        assert not is_inside_triangle((0.5, 1.5), triangle)
     410        assert not is_inside_triangle((0.5, -0.5), triangle)
     411        assert not is_inside_triangle((-0.5, 0.5), triangle)
     412        assert not is_inside_triangle((1.5, 0.5), triangle)
     413
     414        # Try point on borders
     415        assert is_inside_triangle((0.5, 0), triangle, closed=True)
     416        assert is_inside_triangle((1, 0), triangle, closed=True)
     417
     418        assert not is_inside_triangle((0.5, 0), triangle, closed=False)
     419        assert not is_inside_triangle((1, 0), triangle, closed=False)
     420
     421        # Try vertices
     422        for P in triangle:
     423            assert is_inside_triangle(P, triangle, closed=True)
     424            assert not is_inside_triangle(P, triangle, closed=False)
     425
     426        # Slightly different
     427        triangle = [[0, 0.1], [1, -0.2], [0.5, 1]]
     428        assert is_inside_triangle((0.5, 0.5), triangle)
     429        assert is_inside_triangle((0.4, 0.1), triangle)
     430        assert not is_inside_triangle((1, 1), triangle)
     431
     432        # Try vertices
     433        for P in triangle:
     434            assert is_inside_triangle(P, triangle, closed=True)
     435            assert not is_inside_triangle(P, triangle, closed=False)
     436
    373437
    374438    def test_populate_polygon(self):
    … …  
    379443        for point in points:
    380444            assert is_inside_polygon(point, polygon)
     445
     446            assert is_inside_polygon_quick(point, polygon)
     447
    381448
    382449        # Very convoluted polygon
    … …  
    386453        for point in points:
    387454            assert is_inside_polygon(point, polygon)
     455            assert is_inside_polygon_quick(point, polygon)
     456
    388457
    389458    def test_populate_polygon_with_exclude(self):
    … …  
    16861755        assert num.allclose(z, z_ref)
    16871756
    1688     def test_inside_polygon_badtest_1(self):
    1689         '''Test failure found in a larger test in shallow_water.
    1690         (test_get_maximum_inundation in test_data_manager.py (line 10574))
    1691        
    1692         Data below, when fed into:
    1693             indices = inside_polygon(points, polygon)
    1694         returns indices == [], and it shouldn't.
    1695 
    1696         However, it works fine here!?!?
    1697         '''
    1698 
    1699         polygon = [[50, 1], [99, 1], [99, 49], [50, 49]]
    1700 
    1701         points = [[  0.,  0.], [  0., 10.], [  0., 20.], [  0., 30.],
    1702                   [  0., 40.], [  0., 50.], [  5.,  0.], [  5., 10.],
    1703                   [  5., 20.], [  5., 30.], [  5., 40.], [  5., 50.],
    1704                   [ 10.,  0.], [ 10., 10.], [ 10., 20.], [ 10., 30.],
    1705                   [ 10., 40.], [ 10., 50.], [ 15.,  0.], [ 15., 10.],
    1706                   [ 15., 20.], [ 15., 30.], [ 15., 40.], [ 15., 50.],
    1707                   [ 20.,  0.], [ 20., 10.], [ 20., 20.], [ 20., 30.],
    1708                   [ 20., 40.], [ 20., 50.], [ 25.,  0.], [ 25., 10.],
    1709                   [ 25., 20.], [ 25., 30.], [ 25., 40.], [ 25., 50.],
    1710                   [ 30.,  0.], [ 30., 10.], [ 30., 20.], [ 30., 30.],
    1711                   [ 30., 40.], [ 30., 50.], [ 35.,  0.], [ 35., 10.],
    1712                   [ 35., 20.], [ 35., 30.], [ 35., 40.], [ 35., 50.],
    1713                   [ 40.,  0.], [ 40., 10.], [ 40., 20.], [ 40., 30.],
    1714                   [ 40., 40.], [ 40., 50.], [ 45.,  0.], [ 45., 10.],
    1715                   [ 45., 20.], [ 45., 30.], [ 45., 40.], [ 45., 50.],
    1716                   [ 50.,  0.], [ 50., 10.], [ 50., 20.], [ 50., 30.],
    1717                   [ 50., 40.], [ 50., 50.], [ 55.,  0.], [ 55., 10.],
    1718                   [ 55., 20.], [ 55., 30.], [ 55., 40.], [ 55., 50.],
    1719                   [ 60.,  0.], [ 60., 10.], [ 60., 20.], [ 60., 30.],
    1720                   [ 60., 40.], [ 60., 50.], [ 65.,  0.], [ 65., 10.],
    1721                   [ 65., 20.], [ 65., 30.], [ 65., 40.], [ 65., 50.],
    1722                   [ 70.,  0.], [ 70., 10.], [ 70., 20.], [ 70., 30.],
    1723                   [ 70., 40.], [ 70., 50.], [ 75.,  0.], [ 75., 10.],
    1724                   [ 75., 20.], [ 75., 30.], [ 75., 40.], [ 75., 50.],
    1725                   [ 80.,  0.], [ 80., 10.], [ 80., 20.], [ 80., 30.],
    1726                   [ 80., 40.], [ 80., 50.], [ 85.,  0.], [ 85., 10.],
    1727                   [ 85., 20.], [ 85., 30.], [ 85., 40.], [ 85., 50.],
    1728                   [ 90.,  0.], [ 90., 10.], [ 90., 20.], [ 90., 30.],
    1729                   [ 90., 40.], [ 90., 50.], [ 95.,  0.], [ 95., 10.],
    1730                   [ 95., 20.], [ 95., 30.], [ 95., 40.], [ 95., 50.],
    1731                   [100.,  0.], [100., 10.], [100., 20.], [100., 30.],
    1732                   [100., 40.], [100., 50.]]
    1733 
    1734         points = num.array(points)
    1735 
    1736         indices = inside_polygon(points, polygon)
    1737 
     1757
     1758    def test_is_inside_triangle_more(self):
     1759
     1760        res = is_inside_triangle([0.5, 0.5], [[ 0.5,  0. ],
     1761                                              [ 0.5,  0.5],
     1762                                              [ 0.,   0. ]])
     1763        assert res is True
     1764
     1765        res = is_inside_triangle([0.59999999999999998, 0.29999999999999999],
     1766                                 [[ 0.5,  0. ], [ 0.5, 0.5], [0., 0.]])
     1767        assert res is False
     1768
     1769        res = is_inside_triangle([0.59999999999999998, 0.29999999999999999],
     1770                                 [[1., 0.], [1., 0.5], [0.5, 0.]])
     1771        assert res is False
     1772
     1773
     1774        res = is_inside_triangle([0.59999999999999998, 0.29999999999999999],
     1775                                 [[0.5, 0.5], [0.5, 0.], [1., 0.5]])
     1776        assert res is True
     1777
     1778
     1779        res = is_inside_triangle([0.10000000000000001, 0.20000000000000001],
     1780                                 [[0.5, 0.], [0.5, 0.5], [0., 0.]])
     1781        assert res is False
     1782
     1783
     1784        res = is_inside_triangle([0.10000000000000001, 0.20000000000000001],
     1785                                 [[0., 0.5], [0., 0.], [0.5, 0.5]])
     1786        assert res is True
     1787
     1788        res = is_inside_triangle([0.69999999999999996, 0.69999999999999996],
     1789                                 [[0.5, 0.], [0.5, 0.5], [0., 0.]])
     1790        assert res is False
     1791
     1792        res = is_inside_triangle([0.59999999999999998, 0.29999999999999999],
     1793                                 [[0.25, 0.5], [0.25, 0.25], [0.5, 0.5]])
     1794        assert res is False
     1795
     1796
     1797        res = is_inside_triangle([10, 3],
     1798                                 [[0.1, 0.],
     1799                                  [0.1, 0.08333333],
     1800                                  [0., 0.]])
     1801        assert res is False
    17381802
    17391803################################################################################
Note: See TracChangeset for help on using the changeset viewer.