source: trunk/anuga_core/source/anuga/pmesh/mesh_quadtree.py @ 7778

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

Cleaned up unit tests to use API.

File size: 6.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"""
8import time
9
10from anuga.utilities.numerical_tools import get_machine_precision
11from anuga.utilities.numerical_tools import ensure_numeric
12from anuga.config import max_float
13
14from anuga.geometry.quad import Cell
15from anuga.geometry.aabb import AABB
16
17from anuga.utilities import compile
18if compile.can_use_C_extension('polygon_ext.c'):
19    # Underlying C implementations can be accessed
20    from polygon_ext import _is_inside_triangle       
21else:
22    msg = 'C implementations could not be accessed by %s.\n ' %__file__
23    msg += 'Make sure compile_all.py has been run as described in '
24    msg += 'the ANUGA installation guide.'
25    raise Exception, msg
26
27import numpy as num
28 
29
30LAST_TRIANGLE = [[[-1, num.array([[max_float, max_float],
31                               [max_float, max_float],
32                               [max_float, max_float]]),
33                      num.array([[max_float, max_float],
34                               [max_float, max_float],
35                               [max_float, max_float]])], -10]]
36
37
38
39class MeshQuadtree(Cell):
40    """ A quadtree constructed from the given mesh.
41        This class is the root node of a quadtree,
42        and derives from a Cell.
43        It contains optimisations and search patterns specific to meshes.
44    """
45    def __init__(self, mesh):
46        """Build quad tree for mesh.
47
48        All vertex indices in the mesh are stored in a quadtree.
49        """
50       
51        extents = AABB(*mesh.get_extent(absolute=True))   
52        extents.grow(1.001) # To avoid round off error
53        Cell.__init__(self, extents, None)  # root has no parent
54       
55        N = len(mesh)
56        self.mesh = mesh
57        self.set_last_triangle() 
58
59        # Get x,y coordinates for all vertices for all triangles
60        V = mesh.get_vertex_coordinates(absolute=True)
61       
62        normals = mesh.get_normals()
63       
64        # Check each triangle
65        for i in range(N):
66            i3 = 3*i
67            x0, y0 = V[i3, :]
68            x1, y1 = V[i3+1, :]
69            x2, y2 = V[i3+2, :]
70
71            node_data = [i, V[i3:i3+3,:], normals[i,:]]
72
73            # insert a tuple with an AABB, and the triangle index as data
74            self.insert_item((AABB(min([x0, x1, x2]), max([x0, x1, x2]), \
75                             min([y0, y1, y2]), max([y0, y1, y2])), \
76                             node_data))
77
78    def search_fast(self, x):
79        """
80        Find the triangle (element) that the point x is in.
81
82        Inputs:
83            x:    The point to test
84       
85        Return:
86            element_found, sigma0, sigma1, sigma2, k
87
88            where
89            element_found: True if a triangle containing x was found
90            sigma0, sigma1, sigma2: The interpolated values
91            k: Index of triangle (if found)
92
93        """
94       
95        x = ensure_numeric(x, num.float)
96                 
97        # check the last triangle found first
98        element_found, sigma0, sigma1, sigma2, k = \
99                   self._search_triangles_of_vertices(self.last_triangle, x)
100        if element_found:
101            return True, sigma0, sigma1, sigma2, k
102
103        branch = self.last_triangle[0][1]
104
105        # test neighbouring tris
106        tri_data = branch.test_leaves(x)         
107        element_found, sigma0, sigma1, sigma2, k = \
108                    self._search_triangles_of_vertices(tri_data, x)
109        if element_found:
110            return True, sigma0, sigma1, sigma2, k       
111
112        # search rest of tree
113        element_found = False
114        next_search = [branch]
115        while branch:               
116            for sibling in next_search:
117                tri_data = sibling.search(x)         
118                element_found, sigma0, sigma1, sigma2, k = \
119                            self._search_triangles_of_vertices(tri_data, x)
120                if element_found:
121                    return True, sigma0, sigma1, sigma2, k
122           
123            next_search = branch.get_siblings()                           
124            branch = branch.parent
125            if branch:
126                tri_data = branch.test_leaves(x)     
127                element_found, sigma0, sigma1, sigma2, k = \
128                            self._search_triangles_of_vertices(tri_data, x)
129                if element_found:
130                    return True, sigma0, sigma1, sigma2, k     
131
132        return element_found, sigma0, sigma1, sigma2, k
133
134
135    def _search_triangles_of_vertices(self, triangles, x):
136        """Search for triangle containing x among triangle list
137
138        This is called by search_tree_of_vertices once the appropriate node
139        has been found from the quad tree.
140       
141        Input check disabled to speed things up.
142       
143        x is the point to test
144        triangles is the triangle list
145        return the found triangle and its interpolation sigma.
146        """ 
147
148        for node_data in triangles:             
149            if bool(_is_inside_triangle(x, node_data[0][1], \
150                        int(True), 1.0e-12, 1.0e-12)):
151                normals = node_data[0][2]     
152                n0 = normals[0:2]
153                n1 = normals[2:4]
154                n2 = normals[4:6]         
155                xi0, xi1, xi2 = node_data[0][1]
156
157                sigma0 = num.dot((x-xi1), n0)/num.dot((xi0-xi1), n0)
158                sigma1 = num.dot((x-xi2), n1)/num.dot((xi1-xi2), n1)
159                sigma2 = num.dot((x-xi0), n2)/num.dot((xi2-xi0), n2)
160
161                # Don't look for any other triangles in the triangle list
162                self.last_triangle = [node_data]
163                return True, sigma0, sigma1, sigma2, node_data[0][0] # tri index
164        return False, -1, -1, -1, -10
165
166       
167       
168    def set_last_triangle(self):
169        """ Reset last triangle.
170            The algorithm is optimised to find nearby triangles to the
171            previously found one. This is called to reset the search to
172            the root of the tree.
173        """
174        self.last_triangle = LAST_TRIANGLE
175        self.last_triangle[0][1] = self # point at root by default         
176
177   
178
179               
Note: See TracBrowser for help on using the repository browser.