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

Last change on this file since 4653 was 4653, checked in by duncan, 17 years ago

checking in for benchmarking. When fitting cell data - triangle vertices and norms - are calculated the first time a point is looked for in a cell. This is to speed thing up.

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"""
8from Numeric import dot
9
10from anuga.utilities.numerical_tools import get_machine_precision
11
12def search_tree_of_vertices(root, mesh, x):
13    """
14    Find the triangle (element) that the point x is in.
15
16    Inputs:
17        root: A quad tree of the vertices
18        mesh: The underlying mesh
19        x:    The point being placed
20   
21    Return:
22        element_found, sigma0, sigma1, sigma2, k
23
24        where
25        element_found: True if a triangle containing x was found
26        sigma0, sigma1, sigma2: The interpolated values
27        k: Index of triangle (if found)
28
29    """
30    #Find triangle containing x:
31    element_found = False
32
33    # This will be returned if element_found = False
34    sigma2 = -10.0
35    sigma0 = -10.0
36    sigma1 = -10.0
37    k = -10.0
38           
39    #Get triangles in the cell that the point is in.
40    # Triangle is a list, first element triangle_id,
41    # second element the triangle
42    triangles = root.search(x[0], x[1])
43    is_more_elements = True
44
45    element_found, sigma0, sigma1, sigma2, k = \
46                   _search_triangles_of_vertices(mesh,
47                                                 triangles, x)
48    while not element_found and is_more_elements:
49        triangles, branch = root.expand_search()
50        if branch == []:
51            # Searching all the verts from the root cell that haven't
52            # been searched.  This is the last try
53            element_found, sigma0, sigma1, sigma2, k = \
54                           _search_triangles_of_vertices(mesh,triangles, x)
55            is_more_elements = False
56        else:
57            element_found, sigma0, sigma1, sigma2, k = \
58                       _search_triangles_of_vertices(mesh,triangles, x)
59
60    return element_found, sigma0, sigma1, sigma2, k
61
62
63def _search_triangles_of_vertices(mesh, triangles, x):
64    """Search for triangle containing x amongs candidate_vertices in mesh
65
66    This is called by search_tree_of_vertices once the appropriate node
67    has been found from the quad tree.
68   
69
70    This function is responsible for most of the compute time in
71    fit and interpolate.
72    """
73   
74    #Find triangle containing x:
75    element_found = False
76
77    # This will be returned if element_found = False
78    sigma2 = -10.0
79    sigma0 = -10.0
80    sigma1 = -10.0
81    k = -10
82   
83    #For all vertices in same cell as point x
84    for k, tri_verts_norms in triangles:
85        tri = tri_verts_norms[0]
86        n0, n1, n2 = tri_verts_norms[1]
87        # k is the triangle index
88        # tri is a list of verts (x, y), representing a tringle
89        # Find triangle that contains x (if any) and interpolate
90        element_found, sigma0, sigma1, sigma2 =\
91                       find_triangle_compute_interpolation(tri, n0, n1, n2, x)
92        if element_found is True:
93            # Don't look for any other triangles in the triangle list
94            break
95    return element_found, sigma0, sigma1, sigma2, k
96
97
98           
99def find_triangle_compute_interpolation(triangle, n0, n1, n2, x):
100    """Compute linear interpolation of point x and triangle k in mesh.
101    It is assumed that x belongs to triangle k.
102    """
103
104    # Get the three vertex_points of candidate triangle k
105    xi0, xi1, xi2 = triangle
106
107    # this is where we can call some fast c code.
108     
109    # Integrity check - machine precision is too hard
110    # so we use hardwired single precision
111    epsilon = 1.0e-6
112   
113    xmax = max(xi0[0], xi1[0], xi2[0])
114    xmin = min(xi0[0], xi1[0], xi2[0])
115    ymax = max(xi0[1], xi1[1], xi2[1])
116    ymin = min(xi0[1], xi1[1], xi2[1])
117
118   
119    if  x[0] > xmax + epsilon:
120        return False,0,0,0
121    if  x[0] < xmin - epsilon:
122        return False,0,0,0
123    if  x[1] > ymax + epsilon:
124        return False,0,0,0
125    if  x[1] < ymin - epsilon:
126        return False,0,0,0
127   
128    # Get the three normals
129    #n0 = norms[0] 
130    #n1 = norms[1]
131    #n2 = norms[2]
132       
133    # Compute interpolation
134    sigma2 = dot((x-xi0), n2)/dot((xi2-xi0), n2)
135    sigma0 = dot((x-xi1), n0)/dot((xi0-xi1), n0)
136    sigma1 = dot((x-xi2), n1)/dot((xi1-xi2), n1)
137
138    delta = abs(sigma0+sigma1+sigma2-1.0) # Should be close to zero
139    msg = 'abs(sigma0+sigma1+sigma2-1) = %.15e, eps = %.15e'\
140          %(delta, epsilon)
141    assert delta < epsilon, msg
142
143
144    # Check that this triangle contains the data point
145    # Sigmas are allowed to get negative within
146    # machine precision on some machines (e.g. nautilus)
147    epsilon = get_machine_precision() * 2
148    if sigma0 >= -epsilon and sigma1 >= -epsilon and sigma2 >= -epsilon:
149        element_found = True
150    else:
151        element_found = False 
152    return element_found, sigma0, sigma1, sigma2
Note: See TracBrowser for help on using the repository browser.