source: branches/numpy/anuga/fit_interpolate/search_functions.py @ 6883

Last change on this file since 6883 was 6553, checked in by rwilson, 16 years ago

Merged trunk into numpy, all tests and validations work.

File size: 4.8 KB
RevLine 
[6152]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
[6553]10from anuga.utilities.polygon import is_inside_triangle
[6152]11from anuga.utilities.numerical_tools import get_machine_precision
12from anuga.config import max_float
13
[6304]14import numpy as num
[6152]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
[6553]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])))]]
[6152]29
[6553]30def search_tree_of_vertices(root, x):
[6152]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
[6553]51    element_found, sigma0, sigma1, sigma2, k = \
52        _search_triangles_of_vertices(last_triangle, x)
[6152]53                   
54    if element_found is True:
55        return element_found, sigma0, sigma1, sigma2, k
56
[6553]57   
[6152]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])
[6553]62    element_found, sigma0, sigma1, sigma2, k = \
63                   _search_triangles_of_vertices(triangles, x)
64
[6152]65    is_more_elements = True
66   
67    while not element_found and is_more_elements:
68        triangles, branch = root.expand_search()
69        if branch == []:
70            # Searching all the verts from the root cell that haven't
71            # been searched.  This is the last try
72            element_found, sigma0, sigma1, sigma2, k = \
[6553]73                           _search_triangles_of_vertices(triangles, x)
[6152]74            is_more_elements = False
75        else:
76            element_found, sigma0, sigma1, sigma2, k = \
[6553]77                       _search_triangles_of_vertices(triangles, x)
78                       
[6152]79       
80    return element_found, sigma0, sigma1, sigma2, k
81
82
[6553]83def _search_triangles_of_vertices(triangles, x):
84    """Search for triangle containing x amongs candidate_vertices in triangles
[6152]85
86    This is called by search_tree_of_vertices once the appropriate node
87    has been found from the quad tree.
88   
89
90    This function is responsible for most of the compute time in
91    fit and interpolate.
92    """
93    global last_triangle
94   
[6553]95    # These statments are needed if triangles is empty
[6152]96    sigma2 = -10.0
97    sigma0 = -10.0
98    sigma1 = -10.0
99    k = -10
[6553]100    # For all vertices in same cell as point x
101    element_found = False   
[6152]102    for k, tri_verts_norms in triangles:
103        tri = tri_verts_norms[0]
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
[6553]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
[6152]119            # Don't look for any other triangles in the triangle list
[6553]120            last_triangle = [[k, tri_verts_norms]]
121            break           
122           
[6152]123    return element_found, sigma0, sigma1, sigma2, k
124
125
126           
[6553]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.
[6152]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
[6553]140    return sigma0, sigma1, sigma2
[6152]141
142def set_last_triangle():
143    global last_triangle
144    last_triangle = LAST_TRIANGLE
145   
146def search_times():
147
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
155    global search_one_cell_time
156    global search_more_cells_time
157    search_one_cell_time = initial_search_value
158    search_more_cells_time = initial_search_value
Note: See TracBrowser for help on using the repository browser.