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

Merged trunk into numpy, all tests and validations work.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • 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():
Note: See TracChangeset for help on using the changeset viewer.