Changeset 6553 for branches/numpy/anuga/fit_interpolate/search_functions.py
 Timestamp:
 Mar 19, 2009, 1:43:34 PM (15 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

branches/numpy/anuga/fit_interpolate/search_functions.py
r6304 r6553 8 8 import time 9 9 10 from anuga.utilities.polygon import is_inside_triangle 10 11 from anuga.utilities.numerical_tools import get_machine_precision 11 12 from anuga.config import max_float … … 18 19 search_more_cells_time = initial_search_value 19 20 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? 22 LAST_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])))]] 27 29 28 def search_tree_of_vertices(root, mesh,x):30 def search_tree_of_vertices(root, x): 29 31 """ 30 32 Find the triangle (element) that the point x is in. … … 32 34 Inputs: 33 35 root: A quad tree of the vertices 34 mesh: The underlying mesh35 36 x: The point being placed 36 37 … … 47 48 global search_more_cells_time 48 49 49 #Find triangle containing x:50 element_found = False51 52 # This will be returned if element_found = False53 sigma2 = 10.054 sigma0 = 10.055 sigma1 = 10.056 k = 10.057 58 50 # 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) 64 53 65 #print "last_triangle", last_triangle66 54 if element_found is True: 67 #print "last_triangle", last_triangle68 55 return element_found, sigma0, sigma1, sigma2, k 69 56 70 # This was only slightly faster than just checking the 71 # last triangle and it significantly slowed down 72 # nongridded 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 93 58 # Get triangles in the cell that the point is in. 94 59 # Triangle is a list, first element triangle_id, 95 60 # second element the triangle 96 61 triangles = root.search(x[0], x[1]) 62 element_found, sigma0, sigma1, sigma2, k = \ 63 _search_triangles_of_vertices(triangles, x) 64 97 65 is_more_elements = True 98 66 99 element_found, sigma0, sigma1, sigma2, k = \100 _search_triangles_of_vertices(mesh,101 triangles, x)102 #search_one_cell_time += time.time()t0103 #print "search_one_cell_time",search_one_cell_time104 #t0 = time.time()105 67 while not element_found and is_more_elements: 106 68 triangles, branch = root.expand_search() … … 109 71 # been searched. This is the last try 110 72 element_found, sigma0, sigma1, sigma2, k = \ 111 _search_triangles_of_vertices( mesh,triangles, x)73 _search_triangles_of_vertices(triangles, x) 112 74 is_more_elements = False 113 75 else: 114 76 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 118 79 119 80 return element_found, sigma0, sigma1, sigma2, k 120 81 121 82 122 def _search_triangles_of_vertices( mesh,triangles, x):123 """Search for triangle containing x amongs candidate_vertices in mesh83 def _search_triangles_of_vertices(triangles, x): 84 """Search for triangle containing x amongs candidate_vertices in triangles 124 85 125 86 This is called by search_tree_of_vertices once the appropriate node … … 132 93 global last_triangle 133 94 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 139 96 sigma2 = 10.0 140 97 sigma0 = 10.0 141 98 sigma1 = 10.0 142 99 k = 10 143 144 #For all vertices in same cell as point x100 # For all vertices in same cell as point x 101 element_found = False 145 102 for k, tri_verts_norms in triangles: 146 103 tri = tri_verts_norms[0] 147 n0, n1, n2 = tri_verts_norms[1]148 104 # k is the triangle index 149 105 # tri is a list of verts (x, y), representing a tringle 150 106 # 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 154 119 # 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 157 123 return element_found, sigma0, sigma1, sigma2, k 158 124 159 125 160 126 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 127 def 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. 164 131 """ 165 132 … … 167 134 xi0, xi1, xi2 = triangle 168 135 169 # this is where we can call some fast c code.170 171 # Integrity check  machine precision is too hard172 # so we use hardwired single precision173 epsilon = 1.0e6174 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,0178 if x[0] < min(xi0[0], xi1[0], xi2[0])  epsilon:179 return False,0,0,0180 if x[1] > max(xi0[1], xi1[1], xi2[1]) + epsilon:181 return False,0,0,0182 if x[1] < min(xi0[1], xi1[1], xi2[1])  epsilon:183 return False,0,0,0184 185 # machine precision on some machines (e.g. nautilus)186 epsilon = get_machine_precision() * 2187 188 # Compute interpolation  return as soon as possible189 # print "(xi0xi1)", (xi0xi1)190 # print "n0", n0191 # print "dot((xi0xi1), n0)", dot((xi0xi1), n0)192 193 136 sigma0 = num.dot((xxi1), n0)/num.dot((xi0xi1), n0) 194 if sigma0 < epsilon:195 return False,0,0,0196 137 sigma1 = num.dot((xxi2), n1)/num.dot((xi1xi2), n1) 197 if sigma1 < epsilon:198 return False,0,0,0199 138 sigma2 = num.dot((xxi0), n2)/num.dot((xi2xi0), n2) 200 if sigma2 < epsilon:201 return False,0,0,0202 203 # epsilon = 1.0e6204 # we want to speed this up, so don't do assertions205 #delta = abs(sigma0+sigma1+sigma21.0) # Should be close to zero206 #msg = 'abs(sigma0+sigma1+sigma21) = %.15e, eps = %.15e'\207 # %(delta, epsilon)208 #assert delta < epsilon, msg209 139 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 219 141 220 142 def set_last_triangle():
Note: See TracChangeset
for help on using the changeset viewer.