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

Last change on this file since 6541 was 6541, checked in by ole, 15 years ago

Implemented fix to ticket:314.

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