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

Last change on this file since 4599 was 4599, checked in by ole, 17 years ago

More work on search functions

File size: 5.1 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    #Find vertices near x
40    candidate_vertices = root.search(x[0], x[1])
41    is_more_elements = True
42
43    element_found, sigma0, sigma1, sigma2, k = \
44                   _search_triangles_of_vertices(mesh,
45                                                 candidate_vertices, x)
46    while not element_found and is_more_elements:
47        candidate_vertices, branch = root.expand_search()
48        if branch == []:
49            # Searching all the verts from the root cell that haven't
50            # been searched.  This is the last try
51            element_found, sigma0, sigma1, sigma2, k = \
52                           _search_triangles_of_vertices(mesh,
53                                                         candidate_vertices, x)
54            is_more_elements = False
55        else:
56            element_found, sigma0, sigma1, sigma2, k = \
57                       _search_triangles_of_vertices(mesh,
58                                                     candidate_vertices, x)
59
60    return element_found, sigma0, sigma1, sigma2, k
61
62
63def _search_triangles_of_vertices(mesh, candidate_vertices, 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 v in candidate_vertices:
85       
86        #FIXME (DSG-DSG): this catches verts with no triangle.
87        #Currently pmesh is producing these.
88        #this should be stopped,
89
90        if mesh.number_of_triangles_per_node[v] == 0:
91            continue
92       
93        # Get all triangles which has v as a vertex
94        # The list has elements (triangle, vertex), but only the
95        # first component will be used here
96        triangle_list = mesh.get_triangles_and_vertices_per_node(node=v)
97
98        # Find triangle that contains x (if any) and interpolate
99        element_found, sigma0, sigma1, sigma2, k =\
100                       find_triangle_compute_interpolation(mesh,
101                                                           triangle_list,
102                                                           x)
103
104        if element_found is True:
105            # Don't look for any other triangle           
106            break
107
108       
109    return element_found, sigma0, sigma1, sigma2, k
110
111
112           
113def find_triangle_compute_interpolation(mesh, triangle_list, x):
114    """Compute linear interpolation of point x and triangle k in mesh.
115    It is assumed that x belongs to triangle k.
116    """
117
118    element_found = False
119    for k, _ in triangle_list:
120        # Get the three vertex_points of candidate triangle k
121        xi0, xi1, xi2 = mesh.get_vertex_coordinates(triangle_id=k)           
122       
123        # Get the three normals
124        n0 = mesh.get_normal(k, 0)
125        n1 = mesh.get_normal(k, 1)
126        n2 = mesh.get_normal(k, 2)           
127       
128       
129        # Compute interpolation
130        sigma2 = dot((x-xi0), n2)/dot((xi2-xi0), n2)
131        sigma0 = dot((x-xi1), n0)/dot((xi0-xi1), n0)
132        sigma1 = dot((x-xi2), n1)/dot((xi1-xi2), n1)
133
134        # Integrity check - machine precision is too hard
135        # so we use hardwired single precision
136        epsilon = 1.0e-6
137        delta = abs(sigma0+sigma1+sigma2-1.0) # Should be close to zero
138        msg = 'abs(sigma0+sigma1+sigma2-1) = %.15e, eps = %.15e'\
139              %(delta, epsilon)
140        assert delta < epsilon, msg
141
142
143        # Check that this triangle contains the data point
144        # Sigmas are allowed to get negative within
145        # machine precision on some machines (e.g. nautilus)
146        epsilon = get_machine_precision() * 2
147        if sigma0 >= -epsilon and sigma1 >= -epsilon and sigma2 >= -epsilon:
148            element_found = True
149            break
150           
151    return element_found, sigma0, sigma1, sigma2, k
Note: See TracBrowser for help on using the repository browser.