Changeset 6186


Ignore:
Timestamp:
Jan 17, 2009, 3:24:28 PM (10 years ago)
Author:
ole
Message:

A bit of cleanup and renaming

Location:
anuga_core/source/anuga/fit_interpolate
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/fit_interpolate/interpolate.py

    r6185 r6186  
    516516# @param verbose True if this function is to be verbose.
    517517def interpolate_polyline(f,
    518                          vertex_coordinates,
     518                         polyline_nodes,
    519519                         gauge_neighbour_id,
    520520                         point_coordinates=None,
    521521                         verbose=False):
    522     """Interpolate linearly between values f on nodes (vertex coordinates)
     522    """Interpolate linearly between values f on polyline nodes
    523523    of a polyline to midpoints of triangles of boundary.
    524524
     
    539539      Interpolated values at inputted points (z).
    540540    """
    541    
    542     #print f
    543     #print vertex_coordinates
    544     #print gauge_neighbour_id
    545     #print point_coordinates
    546 
    547 
    548541    # FIXME: There is an option of passing a tolerance into this
    549542   
     
    551544        point_coordinates = point_coordinates.get_data_points(absolute=True)
    552545
    553     z = num.ones(len(point_coordinates), num.Float)
    554 
    555     # input sanity check
     546    z = num.zeros(len(point_coordinates), num.Float)
     547
     548    f = ensure_numeric(f, num.Float)
     549    polyline_nodes = ensure_numeric(polyline_nodes, num.Float)
     550    point_coordinates = ensure_numeric(point_coordinates, num.Float)
     551    gauge_neighbour_id = ensure_numeric(gauge_neighbour_id, num.Int)
     552
     553    n = polyline_nodes.shape[0] # Number of nodes in polyline       
     554    # Input sanity check
    556555    msg = 'point coordinates are not given (interpolate.py)'
    557556    assert point_coordinates is not None, msg
    558557    msg = 'function value must be specified at every interpolation node'
    559     assert f.shape[0]==vertex_coordinates.shape[0], msg
     558    assert f.shape[0]==polyline_nodes.shape[0], msg
    560559    msg = 'Must define function value at one or more nodes'
    561560    assert f.shape[0]>0, msg
    562561
    563     n = f.shape[0]
     562
    564563    if n == 1:
    565         z = f*z
    566         msg = 'Polyline contained only one point. I need more. ', str(f)
     564        msg = 'Polyline contained only one point. I need more. ' + str(f)
    567565        raise Exception, msg
    568        
    569     # FIXME (John): add unit test for only One vertex point.
    570     #               Exception should be thrown.
    571    
    572566    elif n > 1:
    573         _interpolate_polyline_aux(n, z, f, point_coordinates, vertex_coordinates, gauge_neighbour_id)
     567        _interpolate_polyline_aux(point_coordinates, len(point_coordinates),
     568                                  polyline_nodes, len(polyline_nodes),
     569                                  f,
     570                                  gauge_neighbour_id,
     571                                  z)
    574572       
    575573    return z
    576574
    577575       
    578 def _interpolate_polyline_aux(n, z, f, point_coordinates, vertex_coordinates, gauge_neighbour_id):
    579     """Auxiliary function
     576def _interpolate_polyline_aux(point_coordinates, number_of_points,
     577                              polyline_nodes, number_of_nodes,
     578                              f, gauge_neighbour_id, z):
     579    """Auxiliary function used by interpolate_polyline
    580580    """
    581     for i in range(len(point_coordinates)):
    582    
    583         x2 = point_coordinates[i][0]
    584         y2 = point_coordinates[i][1]
    585    
     581   
     582    for i in range(number_of_points):
     583
     584        x2, y2 = point_coordinates[i,:]
    586585        found = False
    587         for j in range(n):
     586   
     587        for j in range(number_of_nodes):           
    588588           
    589589            neighbour_id = gauge_neighbour_id[j]
    590590            if neighbour_id >= 0:
    591                 x0 = vertex_coordinates[j][0]
    592                 y0 = vertex_coordinates[j][1]
    593                 x1 = vertex_coordinates[neighbour_id][0]
    594                 y1 = vertex_coordinates[neighbour_id][1]
     591                x0, y0 = polyline_nodes[j,:]
     592                x1, y1 = polyline_nodes[neighbour_id,:]
    595593           
    596                 if point_on_line([x2,y2],
     594                if point_on_line([x2, y2],
    597595                                 [[x0, y0], [x1, y1]],
    598596                                 rtol=1.0e-6):
  • anuga_core/source/anuga/fit_interpolate/test_interpolate.py

    r6185 r6186  
    18961896        #print z
    18971897        assert num.allclose(z, z_ref)               
     1898       
     1899        # Test exception thrown for one point
     1900        f = num.array([5])               
     1901        vertex_coordinates = num.array([[4., 4.]])             
     1902        try: 
     1903            z = interpolate_polyline(f, vertex_coordinates, gauge_neighbour_id, point_coordinates)
     1904        except Exception:
     1905            pass
     1906        else:
     1907            raise Exception, 'One point should have raised exception'
    18981908                       
    18991909#-------------------------------------------------------------
Note: See TracChangeset for help on using the changeset viewer.