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

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

Ticket 327 - marginal 5% speed increase by removing is_inside_triangle wrapper to C function and just calling directly.

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