# source:anuga_core/source/anuga/fit_interpolate/search_functions.py@6174

Last change on this file since 6174 was 6174, checked in by rwilson, 15 years ago

Changed .array(A) to .array(A, num.Int) where appropriate. Helps when going to numpy.

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"""
8import time
9
10from anuga.utilities.numerical_tools import get_machine_precision
11from anuga.config import max_float
12
13import Numeric as num
14
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,[(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], num.Int),      #array default#
25                        num.array([0,0], num.Int),      #array default#
26                        num.array([-1.1,-1.1]))]]]
27
28def search_tree_of_vertices(root, mesh, x):
29    """
30    Find the triangle (element) that the point x is in.
31
32    Inputs:
33        root: A quad tree of the vertices
34        mesh: The underlying mesh
35        x:    The point being placed
36
37    Return:
38        element_found, sigma0, sigma1, sigma2, k
39
40        where
41        element_found: True if a triangle containing x was found
42        sigma0, sigma1, sigma2: The interpolated values
43        k: Index of triangle (if found)
44
45    """
46    global search_one_cell_time
47    global search_more_cells_time
48
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
58    # 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
64
65    #print "last_triangle", last_triangle
66    if element_found is True:
67        #print "last_triangle", last_triangle
68        return element_found, sigma0, sigma1, sigma2, k
69
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()
93    # Get triangles in the cell that the point is in.
94    # Triangle is a list, first element triangle_id,
95    # second element the triangle
96    triangles = root.search(x[0], x[1])
97    is_more_elements = True
98
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()
105    while not element_found and is_more_elements:
106        triangles, branch = root.expand_search()
107        if branch == []:
108            # Searching all the verts from the root cell that haven't
109            # been searched.  This is the last try
110            element_found, sigma0, sigma1, sigma2, k = \
111                           _search_triangles_of_vertices(mesh, triangles, x)
112            is_more_elements = False
113        else:
114            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
118
119    return element_found, sigma0, sigma1, sigma2, k
120
121
122def _search_triangles_of_vertices(mesh, triangles, x):
123    """Search for triangle containing x amongs candidate_vertices in mesh
124
125    This is called by search_tree_of_vertices once the appropriate node
126    has been found from the quad tree.
127
128
129    This function is responsible for most of the compute time in
130    fit and interpolate.
131    """
132    global last_triangle
133
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
139    sigma2 = -10.0
140    sigma0 = -10.0
141    sigma1 = -10.0
142    k = -10
143
144    #For all vertices in same cell as point x
145    for k, tri_verts_norms in triangles:
146        tri = tri_verts_norms[0]
147        n0, n1, n2 = tri_verts_norms[1]
148        # k is the triangle index
149        # tri is a list of verts (x, y), representing a tringle
150        # 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:
154            # Don't look for any other triangles in the triangle list
155            last_triangle = [[k,tri_verts_norms]]
156            break
157    return element_found, sigma0, sigma1, sigma2, k
158
159
160
161def 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
164    """
165
166    # Get the three vertex_points of candidate triangle k
167    xi0, xi1, xi2 = triangle
168
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
193    sigma0 = num.dot((x-xi1), n0)/num.dot((xi0-xi1), n0)
194    if sigma0 < -epsilon:
195        return False,0,0,0
196    sigma1 = num.dot((x-xi2), n1)/num.dot((xi1-xi2), n1)
197    if sigma1 < -epsilon:
198        return False,0,0,0
199    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
209
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
219
220def set_last_triangle():
221    global last_triangle
222    last_triangle = LAST_TRIANGLE
223
224def search_times():
225
226    global search_one_cell_time
227    global search_more_cells_time
228
229    return search_one_cell_time, search_more_cells_time
230
231def reset_search_times():
232
233    global search_one_cell_time
234    global search_more_cells_time
235    search_one_cell_time = initial_search_value
236    search_more_cells_time = initial_search_value
Note: See TracBrowser for help on using the repository browser.