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

Last change on this file since 7707 was 7707, checked in by hudson, 14 years ago

New quadtree implementation - unoptimised and no tree balancing. A couple of failing unit tests to fix.

File size: 6.0 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.utilities.numerical_tools import ensure_numeric
12from anuga.config import max_float
13
14from anuga.utilities import compile
15if compile.can_use_C_extension('polygon_ext.c'):
16    # Underlying C implementations can be accessed
17    from polygon_ext import _is_inside_triangle       
18else:
19    msg = 'C implementations could not be accessed by %s.\n ' %__file__
20    msg += 'Make sure compile_all.py has been run as described in '
21    msg += 'the ANUGA installation guide.'
22    raise Exception, msg
23
24import numpy as num
25
26
27initial_search_value = 'uncomment search_functions code first'#0
28search_one_cell_time = initial_search_value
29search_more_cells_time = initial_search_value
30
31# FIXME(Ole): Could we come up with a less confusing structure?
32# FIXME(James): remove this global var
33LAST_TRIANGLE = [[-10,
34                   (num.array([[max_float, max_float],
35                               [max_float, max_float],
36                               [max_float, max_float]]),
37                    (num.array([1.,1.]),     
38                     num.array([0.,0.]),     
39                     num.array([-1.1,-1.1])))]]
40
41last_triangle = LAST_TRIANGLE                                   
42                                         
43def search_tree_of_vertices(root, mesh, x):
44    """
45    Find the triangle (element) that the point x is in.
46
47    Inputs:
48        root: A quad tree of the vertices
49        mesh: The mesh which the quad tree indexes into
50        x:    The point being placed
51   
52    Return:
53        element_found, sigma0, sigma1, sigma2, k
54
55        where
56        element_found: True if a triangle containing x was found
57        sigma0, sigma1, sigma2: The interpolated values
58        k: Index of triangle (if found)
59
60    """
61    global search_one_cell_time
62    global search_more_cells_time
63
64    # Search the last triangle first
65    element_found, sigma0, sigma1, sigma2, k = \
66        _search_triangles_of_vertices(last_triangle, x)
67                   
68    if element_found is True:
69        return element_found, sigma0, sigma1, sigma2, k
70
71   
72    # Get triangles in the cell that the point is in.
73    tri_indices = root.search(x[0], x[1])
74    triangles = _trilist_from_indices(mesh, tri_indices)
75   
76    element_found, sigma0, sigma1, sigma2, k = \
77                   _search_triangles_of_vertices(triangles, x)
78
79    is_more_elements = True
80   
81    # while not element_found and is_more_elements:
82        # triangles, branch = root.expand_search()
83        # if branch == []:
84            # # Searching all the verts from the root cell that haven't
85            # # been searched.  This is the last try
86            # element_found, sigma0, sigma1, sigma2, k = \
87                           # _search_triangles_of_vertices(triangles, x)
88            # is_more_elements = False
89        # else:
90            # element_found, sigma0, sigma1, sigma2, k = \
91                       # _search_triangles_of_vertices(triangles, x)
92                       
93       
94    return element_found, sigma0, sigma1, sigma2, k
95
96
97def _search_triangles_of_vertices(triangles, x):
98    """Search for triangle containing x amongs candidate_vertices in triangles
99
100    This is called by search_tree_of_vertices once the appropriate node
101    has been found from the quad tree.
102   
103
104    This function is responsible for most of the compute time in
105    fit and interpolate.
106    """
107    global last_triangle
108
109    x = ensure_numeric(x, num.float)     
110   
111    # These statments are needed if triangles is empty
112    sigma2 = -10.0
113    sigma0 = -10.0
114    sigma1 = -10.0
115    k = -10
116    # For all vertices in same cell as point x
117    element_found = False   
118    for k, tri_verts_norms in triangles:
119        tri = tri_verts_norms[0]
120        tri = ensure_numeric(tri)       
121        # k is the triangle index
122        # tri is a list of verts (x, y), representing a triangle
123        # Find triangle that contains x (if any) and interpolate
124       
125        # Input check disabled to speed things up.   
126        if bool(_is_inside_triangle(x, tri, int(True), 1.0e-12, 1.0e-12)):
127           
128            n0, n1, n2 = tri_verts_norms[1]       
129            sigma0, sigma1, sigma2 =\
130                compute_interpolation_values(tri, n0, n1, n2, x)
131               
132            element_found = True   
133
134            # Don't look for any other triangles in the triangle list
135            last_triangle = [[k, tri_verts_norms]]
136            break           
137           
138    return element_found, sigma0, sigma1, sigma2, k
139
140
141def _trilist_from_indices(mesh, indices):
142    """return a list of lists. For the inner lists,
143    The first element is the triangle index,
144    the second element is a list.for this list
145       the first element is a list of three (x, y) vertices,
146       the following elements are the three triangle normals.
147
148    """
149
150    ret_list = []
151    for i in indices:
152        vertices = mesh.get_vertex_coordinates(triangle_id=i, absolute=True)
153        n0 = mesh.get_normal(i, 0)
154        n1 = mesh.get_normal(i, 1)
155        n2 = mesh.get_normal(i, 2) 
156        ret_list.append([i, [vertices, (n0, n1, n2)]])
157    return ret_list
158               
159           
160def compute_interpolation_values(triangle, n0, n1, n2, x):
161    """Compute linear interpolation of point x and triangle.
162   
163    n0, n1, n2 are normal to the tree edges.
164    """
165
166    # Get the three vertex_points of candidate triangle k
167    xi0, xi1, xi2 = triangle
168
169    sigma0 = num.dot((x-xi1), n0)/num.dot((xi0-xi1), n0)
170    sigma1 = num.dot((x-xi2), n1)/num.dot((xi1-xi2), n1)
171    sigma2 = num.dot((x-xi0), n2)/num.dot((xi2-xi0), n2)
172
173    return sigma0, sigma1, sigma2
174
175def set_last_triangle():
176    global last_triangle
177    last_triangle = LAST_TRIANGLE
178   
179def search_times():
180
181    global search_one_cell_time
182    global search_more_cells_time
183
184    return search_one_cell_time, search_more_cells_time
185
186def reset_search_times():
187
188    global search_one_cell_time
189    global search_more_cells_time
190    search_one_cell_time = initial_search_value
191    search_more_cells_time = initial_search_value
Note: See TracBrowser for help on using the repository browser.