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

Last change on this file since 6545 was 6545, checked in by ole, 14 years ago

Removed more cluttter in search functions.

File size: 4.8 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(Ole): Could we come up with a less confusing structure?
22LAST_TRIANGLE = [[-10,
23                   (num.array([[max_float, max_float],
24                               [max_float, max_float],
25                               [max_float, max_float]]),
26                    (num.array([1.,1.]),
27                     num.array([0.,0.]),
28                     num.array([-1.1,-1.1])))]]
29
30def search_tree_of_vertices(root, x):
31    """
32    Find the triangle (element) that the point x is in.
33
34    Inputs:
35        root: A quad tree of the vertices
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    # Search the last triangle first
51    element_found, sigma0, sigma1, sigma2, k = \
52        _search_triangles_of_vertices(last_triangle, x)
53
54    if element_found is True:
55        return element_found, sigma0, sigma1, sigma2, k
56
57
58    # Get triangles in the cell that the point is in.
59    # Triangle is a list, first element triangle_id,
60    # second element the triangle
61    triangles = root.search(x[0], x[1])
62    element_found, sigma0, sigma1, sigma2, k = \
63                   _search_triangles_of_vertices(triangles, x)
64
65    is_more_elements = True
66    while not element_found and is_more_elements:
67        triangles, branch = root.expand_search()
68        if branch == []:
69            # Searching all the verts from the root cell that haven't
70            # been searched.  This is the last try
71            element_found, sigma0, sigma1, sigma2, k = \
72                           _search_triangles_of_vertices(triangles, x)
73            is_more_elements = False
74        else:
75            element_found, sigma0, sigma1, sigma2, k = \
76                       _search_triangles_of_vertices(triangles, x)
77
78
79    return element_found, sigma0, sigma1, sigma2, k
80
81
82def _search_triangles_of_vertices(triangles, x):
83    """Search for triangle containing x amongs candidate_vertices in triangles
84
85    This is called by search_tree_of_vertices once the appropriate node
86    has been found from the quad tree.
87
88    This function is responsible for most of the compute time in
89    fit and interpolate.
90    """
91    global last_triangle
92
93    # These statments are needed if triangles is empty
94    sigma2 = -10.0
95    sigma0 = -10.0
96    sigma1 = -10.0
97    k = -10
98
99    # For all vertices in same cell as point x
100    element_found = False
101    for k, tri_verts_norms in triangles:
102        tri = tri_verts_norms[0]
103
104        # k is the triangle index
105        # tri is a list of verts (x, y), representing a tringle
106        # Find triangle that contains x (if any) and interpolate
107
108        # Input check disabled to speed things up.
109        if is_inside_triangle(x, tri,
110                              closed=True,
111                              check_inputs=False):
112
113            n0, n1, n2 = tri_verts_norms[1]
114            sigma0, sigma1, sigma2 =\
115                compute_interpolation_values(tri, n0, n1, n2, x)
116
117            element_found = True
118
119            # Don't look for any other triangles in the triangle list
120            last_triangle = [[k, tri_verts_norms]]
121            break
122
123
124    return element_found, sigma0, sigma1, sigma2, k
125
126
127def compute_interpolation_values(triangle, n0, n1, n2, x):
128    """Compute linear interpolation of point x and triangle.
129
130    n0, n1, n2 are normal to the tree edges.
131    """
132
133    # Get the three vertex_points of candidate triangle k
134    xi0, xi1, xi2 = triangle
135
136    sigma0 = num.dot((x-xi1), n0)/num.dot((xi0-xi1), n0)
137    sigma1 = num.dot((x-xi2), n1)/num.dot((xi1-xi2), n1)
138    sigma2 = num.dot((x-xi0), n2)/num.dot((xi2-xi0), n2)
139
140    return sigma0, sigma1, sigma2
141
142
143def set_last_triangle():
144    global last_triangle
145    last_triangle = LAST_TRIANGLE
146
147def search_times():
148    global search_one_cell_time
149    global search_more_cells_time
150
151    return search_one_cell_time, search_more_cells_time
152
153def reset_search_times():
154    global search_one_cell_time
155    global search_more_cells_time
156    search_one_cell_time = initial_search_value
157    search_more_cells_time = initial_search_value
Note: See TracBrowser for help on using the repository browser.