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

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

More swb tests passing. Cleaned up some pylint errors.

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