1 | """ |
---|
2 | General functions used in fit and interpolate. |
---|
3 | |
---|
4 | Ole Nielsen, Stephen Roberts, Duncan Gray |
---|
5 | Geoscience Australia, 2006. |
---|
6 | |
---|
7 | """ |
---|
8 | import time |
---|
9 | |
---|
10 | from anuga.utilities.numerical_tools import get_machine_precision |
---|
11 | from anuga.utilities.numerical_tools import ensure_numeric |
---|
12 | from anuga.config import max_float |
---|
13 | |
---|
14 | from anuga.geometry.quad import Cell |
---|
15 | from anuga.geometry.aabb import AABB |
---|
16 | |
---|
17 | from anuga.utilities import compile |
---|
18 | if compile.can_use_C_extension('polygon_ext.c'): |
---|
19 | # Underlying C implementations can be accessed |
---|
20 | from polygon_ext import _is_inside_triangle |
---|
21 | else: |
---|
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 | |
---|
27 | import numpy as num |
---|
28 | |
---|
29 | |
---|
30 | LAST_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 | |
---|
39 | class 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 | |
---|