source: anuga_core/source/anuga/fit_interpolate/search_functions.py @ 5571

Last change on this file since 5571 was 4872, checked in by duncan, 16 years ago

fitting benchmark update. Investigating adding neighbour triangles in fitting.

File size: 7.9 KB
Line 
1"""
2General functions used in fit and interpolate.
3
4   Ole Nielsen, Stephen Roberts, Duncan Gray
5   Geoscience Australia, 2006.
6
7"""
8from Numeric import dot
9import time
10
11from Numeric import array
12
13from anuga.utilities.numerical_tools import get_machine_precision
14from anuga.config import max_float
15
16initial_search_value = 'uncomment search_functions code first'#0
17search_one_cell_time = initial_search_value
18search_more_cells_time = initial_search_value
19
20#FIXME test what happens if a
21LAST_TRIANGLE = [[-10,[(array([max_float,max_float]),
22                        array([max_float,max_float]),
23                        array([max_float,max_float])),
24                       (array([1,1]),array([0,0]),array([-1.1,-1.1]))]]]
25
26def search_tree_of_vertices(root, mesh, x):
27    """
28    Find the triangle (element) that the point x is in.
29
30    Inputs:
31        root: A quad tree of the vertices
32        mesh: The underlying mesh
33        x:    The point being placed
34   
35    Return:
36        element_found, sigma0, sigma1, sigma2, k
37
38        where
39        element_found: True if a triangle containing x was found
40        sigma0, sigma1, sigma2: The interpolated values
41        k: Index of triangle (if found)
42
43    """
44    global search_one_cell_time
45    global search_more_cells_time
46
47    #Find triangle containing x:
48    element_found = False
49
50    # This will be returned if element_found = False
51    sigma2 = -10.0
52    sigma0 = -10.0
53    sigma1 = -10.0
54    k = -10.0
55
56    # Search the last triangle first
57    element_found, sigma0, sigma1, sigma2, k = \
58                   _search_triangles_of_vertices(mesh,last_triangle, x)
59    #print "last_triangle", last_triangle
60    if element_found is True:
61        #print "last_triangle", last_triangle
62        return element_found, sigma0, sigma1, sigma2, k
63
64    # This was only slightly faster than just checking the
65    # last triangle and it significantly slowed down
66    # non-gridded fitting
67    # If the last element was a dud, search its neighbours
68    #print "last_triangle[0][0]", last_triangle[0][0]
69    #neighbours = mesh.get_triangle_neighbours(last_triangle[0][0])
70    #print "neighbours", neighbours
71    #neighbours = []
72  #   for k in neighbours:
73#         if k >= 0:
74#             tri = mesh.get_vertex_coordinates(k,
75#                                                    absolute=True)
76#             n0 = mesh.get_normal(k, 0)
77#             n1 = mesh.get_normal(k, 1)
78#             n2 = mesh.get_normal(k, 2)
79#             triangle =[[k,(tri, (n0, n1, n2))]]
80#             element_found, sigma0, sigma1, sigma2, k = \
81#                            _search_triangles_of_vertices(mesh,
82#                                                          triangle, x)
83#             if element_found is True:
84#                 return element_found, sigma0, sigma1, sigma2, k
85           
86    #t0 = time.time()
87    # Get triangles in the cell that the point is in.
88    # Triangle is a list, first element triangle_id,
89    # second element the triangle
90    triangles = root.search(x[0], x[1])
91    is_more_elements = True
92   
93    element_found, sigma0, sigma1, sigma2, k = \
94                   _search_triangles_of_vertices(mesh,
95                                                 triangles, x)
96    #search_one_cell_time += time.time()-t0
97    #print "search_one_cell_time",search_one_cell_time
98    #t0 = time.time()
99    while not element_found and is_more_elements:
100        triangles, branch = root.expand_search()
101        if branch == []:
102            # Searching all the verts from the root cell that haven't
103            # been searched.  This is the last try
104            element_found, sigma0, sigma1, sigma2, k = \
105                           _search_triangles_of_vertices(mesh,triangles, x)
106            is_more_elements = False
107        else:
108            element_found, sigma0, sigma1, sigma2, k = \
109                       _search_triangles_of_vertices(mesh,triangles, x)
110        #search_more_cells_time += time.time()-t0
111    #print "search_more_cells_time", search_more_cells_time
112       
113    return element_found, sigma0, sigma1, sigma2, k
114
115
116def _search_triangles_of_vertices(mesh, triangles, x):
117    """Search for triangle containing x amongs candidate_vertices in mesh
118
119    This is called by search_tree_of_vertices once the appropriate node
120    has been found from the quad tree.
121   
122
123    This function is responsible for most of the compute time in
124    fit and interpolate.
125    """
126    global last_triangle
127   
128    # these statments are needed if triangles is empty
129    #Find triangle containing x:
130    element_found = False
131
132    # This will be returned if element_found = False
133    sigma2 = -10.0
134    sigma0 = -10.0
135    sigma1 = -10.0
136    k = -10
137   
138    #For all vertices in same cell as point x
139    for k, tri_verts_norms in triangles:
140        tri = tri_verts_norms[0]
141        n0, n1, n2 = tri_verts_norms[1]
142        # k is the triangle index
143        # tri is a list of verts (x, y), representing a tringle
144        # Find triangle that contains x (if any) and interpolate
145        element_found, sigma0, sigma1, sigma2 =\
146                       find_triangle_compute_interpolation(tri, n0, n1, n2, x)
147        if element_found is True:
148            # Don't look for any other triangles in the triangle list
149            last_triangle = [[k,tri_verts_norms]]
150            break
151    return element_found, sigma0, sigma1, sigma2, k
152
153
154           
155def find_triangle_compute_interpolation(triangle, n0, n1, n2, x):
156    """Compute linear interpolation of point x and triangle k in mesh.
157    It is assumed that x belongs to triangle k.max_float
158    """
159
160    # Get the three vertex_points of candidate triangle k
161    xi0, xi1, xi2 = triangle
162
163    # this is where we can call some fast c code.
164     
165    # Integrity check - machine precision is too hard
166    # so we use hardwired single precision
167    epsilon = 1.0e-6
168   
169    if  x[0] > max(xi0[0], xi1[0], xi2[0]) + epsilon:
170        # print "max(xi0[0], xi1[0], xi2[0])", max(xi0[0], xi1[0], xi2[0])
171        return False,0,0,0
172    if  x[0] < min(xi0[0], xi1[0], xi2[0]) - epsilon:
173        return False,0,0,0
174    if  x[1] > max(xi0[1], xi1[1], xi2[1]) + epsilon:
175        return False,0,0,0
176    if  x[1] < min(xi0[1], xi1[1], xi2[1]) - epsilon:
177        return False,0,0,0
178   
179    # machine precision on some machines (e.g. nautilus)
180    epsilon = get_machine_precision() * 2
181   
182    # Compute interpolation - return as soon as possible
183    #  print "(xi0-xi1)", (xi0-xi1)
184    # print "n0", n0
185    # print "dot((xi0-xi1), n0)", dot((xi0-xi1), n0)
186   
187    sigma0 = dot((x-xi1), n0)/dot((xi0-xi1), n0)
188    if sigma0 < -epsilon:
189        return False,0,0,0
190    sigma1 = dot((x-xi2), n1)/dot((xi1-xi2), n1)
191    if sigma1 < -epsilon:
192        return False,0,0,0
193    sigma2 = dot((x-xi0), n2)/dot((xi2-xi0), n2)
194    if sigma2 < -epsilon:
195        return False,0,0,0
196   
197    # epsilon = 1.0e-6
198    # we want to speed this up, so don't do assertions
199    #delta = abs(sigma0+sigma1+sigma2-1.0) # Should be close to zero
200    #msg = 'abs(sigma0+sigma1+sigma2-1) = %.15e, eps = %.15e'\
201    #      %(delta, epsilon)
202    #assert delta < epsilon, msg
203
204
205    # Check that this triangle contains the data point
206    # Sigmas are allowed to get negative within
207    # machine precision on some machines (e.g. nautilus)
208    #if sigma0 >= -epsilon and sigma1 >= -epsilon and sigma2 >= -epsilon:
209    #    element_found = True
210    #else:
211    #    element_found = False
212    return True, sigma0, sigma1, sigma2
213
214def set_last_triangle():
215    global last_triangle
216    last_triangle = LAST_TRIANGLE
217   
218def search_times():
219
220    global search_one_cell_time
221    global search_more_cells_time
222
223    return search_one_cell_time, search_more_cells_time
224
225def reset_search_times():
226
227    global search_one_cell_time
228    global search_more_cells_time
229    search_one_cell_time = initial_search_value
230    search_more_cells_time = initial_search_value
Note: See TracBrowser for help on using the repository browser.