# Changeset 6544

Ignore:
Timestamp:
Mar 18, 2009, 2:34:46 PM (15 years ago)
Message:

Optimised fitting a bit more and cleaned up in search_functions.

Location:
anuga_core/source/anuga
Files:
5 edited

Unmodified
Removed
• ## anuga_core/source/anuga/fit_interpolate/search_functions.py

 r6541 search_more_cells_time = initial_search_value #FIXME test what happens if a LAST_TRIANGLE = [[-10,[(num.array([max_float,max_float]), num.array([max_float,max_float]), num.array([max_float,max_float])), (num.array([1,1], num.Int),      #array default# num.array([0,0], num.Int),      #array default# num.array([-1.1,-1.1]))]]] LAST_TRIANGLE = [[-10, (num.array([[max_float,max_float], [max_float,max_float], [max_float,max_float]]), (num.array([1,1], num.Int),      #array default# num.array([0,0], num.Int),      #array default# num.array([-1.1,-1.1])))]] def search_tree_of_vertices(root, mesh, x): global search_more_cells_time #Find triangle containing x: element_found = False # This will be returned if element_found = False element_found = False sigma2 = -10.0 sigma0 = -10.0 element_found, sigma0, sigma1, sigma2, k = \ _search_triangles_of_vertices(mesh, last_triangle, x) except: print 'This should never happen:', last_triangle element_found = False #print "last_triangle", last_triangle if element_found is True: #print "last_triangle", last_triangle return element_found, sigma0, sigma1, sigma2, k # This was only slightly faster than just checking the # last triangle and it significantly slowed down # non-gridded fitting # If the last element was a dud, search its neighbours #print "last_triangle[0][0]", last_triangle[0][0] #neighbours = mesh.get_triangle_neighbours(last_triangle[0][0]) #print "neighbours", neighbours #neighbours = [] #   for k in neighbours: #         if k >= 0: #             tri = mesh.get_vertex_coordinates(k, #                                                    absolute=True) #             n0 = mesh.get_normal(k, 0) #             n1 = mesh.get_normal(k, 1) #             n2 = mesh.get_normal(k, 2) #             triangle =[[k,(tri, (n0, n1, n2))]] #             element_found, sigma0, sigma1, sigma2, k = \ #                            _search_triangles_of_vertices(mesh, #                                                          triangle, x) #             if element_found is True: #                 return element_found, sigma0, sigma1, sigma2, k #t0 = time.time() # Get triangles in the cell that the point is in. # Triangle is a list, first element triangle_id, triangles = root.search(x[0], x[1]) is_more_elements = True element_found, sigma0, sigma1, sigma2, k = \ _search_triangles_of_vertices(mesh, global last_triangle # these statments are needed if triangles is empty #Find triangle containing x: # These statments are needed if triangles is empty element_found = False # This will be returned if element_found = False sigma2 = -10.0 sigma0 = -10.0 sigma1 = -10.0 k = -10 #For all vertices in same cell as point x # For all vertices in same cell as point x for k, tri_verts_norms in triangles: tri = tri_verts_norms[0] # tri is a list of verts (x, y), representing a tringle # Find triangle that contains x (if any) and interpolate if is_inside_triangle(x, tri, closed=True): # Input check disabled to speed things up. if is_inside_triangle(x, tri, closed=True, check_inputs=False): n0, n1, n2 = tri_verts_norms[1] element_found, sigma0, sigma1, sigma2 =\ find_triangle_compute_interpolation(tri, n0, n1, n2, x) else: element_found = False sigma0, sigma1, sigma2 =\ compute_interpolation_values(tri, n0, n1, n2, x) element_found = True # Don't look for any other triangles in the triangle list last_triangle = [[k, tri_verts_norms]] break if element_found is True: # Don't look for any other triangles in the triangle list last_triangle = [[k,tri_verts_norms]] break return element_found, sigma0, sigma1, sigma2, k def find_triangle_compute_interpolation(triangle, n0, n1, n2, x): """Compute linear interpolation of point x and triangle k in mesh. It is assumed that x belongs to triangle k.max_float def compute_interpolation_values(triangle, n0, n1, n2, x): """Compute linear interpolation of point x and triangle. n0, n1, n2 are normal to the tree edges. """ xi0, xi1, xi2 = triangle # this is where we can call some fast c code. # Integrity check - machine precision is too hard # so we use hardwired single precision #epsilon = 1.0e-6 sigma0 = num.dot((x-xi1), n0)/num.dot((xi0-xi1), n0) sigma1 = num.dot((x-xi2), n1)/num.dot((xi1-xi2), n1) sigma2 = num.dot((x-xi0), n2)/num.dot((xi2-xi0), n2) return sigma0, sigma1, sigma2 #if  x[0] > max(xi0[0], xi1[0], xi2[0]) + epsilon: #    return False,0,0,0 #if  x[0] < min(xi0[0], xi1[0], xi2[0]) - epsilon: #    return False,0,0,0 #if  x[1] > max(xi0[1], xi1[1], xi2[1]) + epsilon: #    return False,0,0,0 #if  x[1] < min(xi0[1], xi1[1], xi2[1]) - epsilon: #    return False,0,0,0 # machine precision on some machines (e.g. nautilus) epsilon = get_machine_precision() * 2 # Compute interpolation - return as soon as possible #  print "(xi0-xi1)", (xi0-xi1) # print "n0", n0 # print "dot((xi0-xi1), n0)", dot((xi0-xi1), n0) sigma0 = num.dot((x-xi1), n0)/num.dot((xi0-xi1), n0) #if sigma0 < -epsilon: #print 'sigma0', sigma0 #sigma0 = 0.0 #return False,0,0,0 sigma1 = num.dot((x-xi2), n1)/num.dot((xi1-xi2), n1) #if sigma1 < -epsilon: #print 'sigma1', sigma1 #sigma1 = 0.0 #return False,0,0,0 sigma2 = num.dot((x-xi0), n2)/num.dot((xi2-xi0), n2) #if sigma2 < -epsilon: #print 'sigma2', sigma2 #sigma2 = 0.0 #return False,0,0,0 # epsilon = 1.0e-6 # we want to speed this up, so don't do assertions #delta = abs(sigma0+sigma1+sigma2-1.0) # Should be close to zero #msg = 'abs(sigma0+sigma1+sigma2-1) = %.15e, eps = %.15e'\ #      %(delta, epsilon) #assert delta < epsilon, msg # Check that this triangle contains the data point # Sigmas are allowed to get negative within # machine precision on some machines (e.g. nautilus) #if sigma0 >= -epsilon and sigma1 >= -epsilon and sigma2 >= -epsilon: #    element_found = True #else: #    element_found = False return True, sigma0, sigma1, sigma2 def set_last_triangle(): global last_triangle
• ## anuga_core/source/anuga/fit_interpolate/test_search_functions.py

 r6541 from search_functions import search_tree_of_vertices, set_last_triangle from search_functions import _search_triangles_of_vertices from search_functions import find_triangle_compute_interpolation from search_functions import compute_interpolation_values from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh from anuga.utilities.polygon import is_inside_polygon from anuga.utilities.quad import build_quadtree, Cell from anuga.utilities.numerical_tools import ensure_numeric from anuga.utilities.polygon import is_inside_polygon, is_inside_triangle x = [0.2, 0.7] found, s0, s1, s2, k = search_tree_of_vertices(root, mesh, x) found, s0, s1, s2, k = search_tree_of_vertices(root, mesh, ensure_numeric(x)) assert k == 1 # Triangle one assert found is True def test_bigger(self): """test_larger mesh """test_bigger test larger mesh """ [0.1,0.9], [0.4,0.6], [0.9,0.1], [10, 3]]: found, s0, s1, s2, k = search_tree_of_vertices(root, mesh, x) found, s0, s1, s2, k = search_tree_of_vertices(root, mesh, ensure_numeric(x)) if k >= 0: mesh.check_integrity() for m in range(8): root = build_quadtree(mesh, max_points_per_cell = m) [10, 3]]: found, s0, s1, s2, k = search_tree_of_vertices(root, mesh, x) found, s0, s1, s2, k = search_tree_of_vertices(root, mesh, x) if k >= 0: V = mesh.get_vertex_coordinates(k) # nodes for triangle k assert is_inside_triangle(x, V, closed=True) assert is_inside_polygon(x, V) assert found is True else: assert found is False if m == k == 0: return def test_underlying_function(self): # One point x = [0.5, 0.5] x = ensure_numeric([0.5, 0.5]) candidate_vertices = root.search(x[0], x[1]) #print x, candidate_vertices # print x, candidate_vertices found, sigma0, sigma1, sigma2, k = \ _search_triangles_of_vertices(mesh, _search_triangles_of_vertices(mesh, candidate_vertices, x) ensure_numeric(x)) if k >= 0: V = mesh.get_vertex_coordinates(k) # nodes for triangle k def test_triangle_compute_interpolation(self): """test_triangle_compute_interpolation Test that triangle can be found if point is inside it def test_compute_interpolation_values(self): """test_compute_interpolation_values Test that interpolation values are correc This test used to check element_found as output from find_triangle_compute_interpolation(triangle, n0, n1, n2, x) and that failed before 18th March 2009. Now this function no longer returns this flag, so the test is merely checknig the sigmas. """ closed=True, verbose=False) assert is_inside_triangle(x, triangle, closed=True, verbose=False) element_found, sigma0, sigma1, sigma2 = \ find_triangle_compute_interpolation(triangle, n0, n1, n2, x) closed=True, verbose=False) sigma0, sigma1, sigma2 = \ compute_interpolation_values(triangle, n0, n1, n2, x) msg = 'Point which is clearly inside triangle was not found' assert element_found is True, msg #print sigma0, sigma1, sigma2 #print sigma0 + sigma1 + sigma2 assert abs(sigma0 + sigma1 + sigma2 - 1) < 1.0e-6 assert abs(sigma0 + sigma1 + sigma2 - 1) < 1.0e-6, msg #-------------------------------------------------------------
• ## anuga_core/source/anuga/utilities/polygon.py

 r6541 """ triangle = ensure_numeric(triangle) point = ensure_numeric(point, num.Float) if check_inputs is True: triangle = ensure_numeric(triangle) point = ensure_numeric(point, num.Float) # Convert point to array of points msg = 'is_inside_triangle must be invoked with one point only' assert num.allclose(point.shape, [2]), msg # Use C-implementation return bool(_is_inside_triangle(point, triangle, int(closed), rtol, atol)) # FIXME (Ole): The rest of this function has been made # obsolete by the C extension. # Quickly reject points that are clearly outside if point[0] < min(triangle[:,0]): return False if point[1] < min(triangle[:,1]): return False if point[1] > max(triangle[:,1]): return False # Use C-implementation return bool(_is_inside_triangle(point, triangle, int(closed), rtol, atol)) # FIXME (Ole): The rest of this function has been made # obsolete by the C extension. # Start search
• ## anuga_core/source/anuga/utilities/polygon_ext.c

 r6535 double a00, a10, a01, a11, b0, b1; double denom, alpha, beta; double x, y; // Point coordinates int i, j, res; x = point[0]; y = point[1]; // Quickly reject points that are clearly outside if ((x < triangle[0]) && (x < triangle[2]) && (x < triangle[4])) return 0; if ((x > triangle[0]) && (x > triangle[2]) && (x > triangle[4])) return 0; if ((y < triangle[1]) && (y < triangle[3]) && (y < triangle[5])) return 0; if ((y > triangle[1]) && (y > triangle[3]) && (y > triangle[5])) return 0; // v0 = C-A v0x = triangle[4]-triangle[0]; if (fabs(denom) > 0.0) { // v = point-A vx = point[0] - triangle[0]; vy = point[1] - triangle[1]; vx = x - triangle[0]; vy = y - triangle[1]; b0 = v0x*vx + v0y*vy; // innerproduct(v0, v) for (i=0; i<3; i++) { j = (i+1) % 3; // Circular index into triangle array res = __point_on_line(point[0], point[1], res = __point_on_line(x, y, triangle[2*i], triangle[2*i+1], triangle[2*j], triangle[2*j+1],
• ## anuga_core/source/anuga/utilities/test_polygon.py

 r6535 assert res is False res = is_inside_triangle([10, 3], [[0.1, 0.], [0.1, 0.08333333], [0., 0.]]) assert res is False #-------------------------------------------------------------
Note: See TracChangeset for help on using the changeset viewer.