Changeset 7716
- Timestamp:
- May 11, 2010, 3:12:21 PM (15 years ago)
- Location:
- anuga_core/source/anuga
- Files:
-
- 5 edited
- 2 moved
-
fit_interpolate/fit.py (modified) (2 diffs)
-
fit_interpolate/general_fit_interpolate.py (modified) (2 diffs)
-
fit_interpolate/interpolate.py (modified) (2 diffs)
-
fit_interpolate/test_search_functions.py (modified) (10 diffs)
-
geometry/mesh_quadtree.py (moved) (moved from anuga_core/source/anuga/fit_interpolate/mesh_quadtree.py) (4 diffs)
-
geometry/test_geometry.py (modified) (1 diff)
-
geometry/test_meshquad.py (moved) (moved from anuga_core/source/anuga/fit_interpolate/test_meshquad.py) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/fit_interpolate/fit.py
r7715 r7716 35 35 from anuga.utilities.sparse import Sparse, Sparse_CSR 36 36 from anuga.geometry.polygon import inside_polygon, is_inside_polygon 37 from anuga.fit_interpolate.mesh_quadtree import search_tree_of_vertices37 from anuga.fit_interpolate.mesh_quadtree import MeshQuadtree 38 38 39 39 from anuga.utilities.cg_solve import conjugate_gradient … … 276 276 277 277 element_found, sigma0, sigma1, sigma2, k = \ 278 se arch_tree_of_vertices(self.root, self.mesh,x)278 self.root.search_fast(x) 279 279 280 280 if element_found is True: -
anuga_core/source/anuga/fit_interpolate/general_fit_interpolate.py
r7715 r7716 31 31 from anuga.geospatial_data.geospatial_data import Geospatial_data, \ 32 32 ensure_absolute 33 from anuga.fit_interpolate.mesh_quadtree import set_last_triangle,MeshQuadtree33 from anuga.fit_interpolate.mesh_quadtree import MeshQuadtree 34 34 import anuga.utilities.log as log 35 35 … … 105 105 106 106 build_quadtree_time = time.time()-t0 107 se t_last_triangle()107 self.root.set_last_triangle() 108 108 109 109 def __repr__(self): -
anuga_core/source/anuga/fit_interpolate/interpolate.py
r7715 r7716 33 33 from anuga.geospatial_data.geospatial_data import Geospatial_data 34 34 from anuga.geospatial_data.geospatial_data import ensure_absolute 35 from anuga.fit_interpolate.mesh_quadtree import search_tree_of_vertices35 from anuga.fit_interpolate.mesh_quadtree import MeshQuadtree 36 36 from anuga.fit_interpolate.general_fit_interpolate import FitInterpolate 37 37 from anuga.abstract_2d_finite_volumes.file_function import file_function … … 492 492 493 493 x = point_coordinates[i] 494 element_found, sigma0, sigma1, sigma2, k = \ 495 search_tree_of_vertices(self.root, self.mesh, x) 494 element_found, sigma0, sigma1, sigma2, k = self.root.search_fast(x) 496 495 497 496 # Update interpolation matrix A if necessary -
anuga_core/source/anuga/fit_interpolate/test_search_functions.py
r7715 r7716 3 3 4 4 import unittest 5 from mesh_quadtree import search_tree_of_vertices, set_last_triangle 6 from mesh_quadtree import _search_triangles_of_vertices 7 from mesh_quadtree import _trilist_from_data 8 from mesh_quadtree import compute_interpolation_values, MeshQuadtree 5 from mesh_quadtree import MeshQuadtree, compute_interpolation_values 9 6 10 7 from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh … … 51 48 52 49 root = MeshQuadtree(mesh) 53 set_last_triangle()54 55 found, s0, s1, s2, k = search_tree_of_vertices(root, mesh,[-0.2, 10.7])50 root.set_last_triangle() 51 52 found, s0, s1, s2, k = root.search_fast([-0.2, 10.7]) 56 53 assert found is False 57 54 58 found, s0, s1, s2, k = search_tree_of_vertices(root, mesh,[0, 0])55 found, s0, s1, s2, k = root.search_fast([0, 0]) 59 56 assert found is True 60 57 … … 71 68 72 69 root = MeshQuadtree(mesh) 73 set_last_triangle()70 root.set_last_triangle() 74 71 75 72 x = [0.2, 0.7] 76 found, s0, s1, s2, k = search_tree_of_vertices(root, mesh,x)73 found, s0, s1, s2, k = root.search_fast(x) 77 74 assert k == 1 # Triangle one 78 75 assert found is True … … 91 88 92 89 root = MeshQuadtree(mesh) 93 set_last_triangle()90 root.set_last_triangle() 94 91 95 92 for x in [[0.6, 0.3], [0.1, 0.2], [0.7,0.7], … … 97 94 [10, 3]]: 98 95 99 found, s0, s1, s2, k = search_tree_of_vertices(root, mesh, 100 ensure_numeric(x)) 96 found, s0, s1, s2, k = root.search_fast(ensure_numeric(x)) 101 97 102 98 if k >= 0: … … 123 119 124 120 root = MeshQuadtree(mesh) 125 set_last_triangle()121 root.set_last_triangle() 126 122 #print m, root.show() 127 123 … … 130 126 [10, 3]]: 131 127 132 found, s0, s1, s2, k = search_tree_of_vertices(root, mesh,x)128 found, s0, s1, s2, k = root.search_fast(x) 133 129 134 130 if k >= 0: … … 151 147 152 148 root = MeshQuadtree(mesh) 153 set_last_triangle()149 root.set_last_triangle() 154 150 155 151 # One point 156 152 x = ensure_numeric([0.5, 0.5]) 157 153 158 triangles = _trilist_from_data(mesh,root.search(x))154 triangles = root._trilist_from_data(root.search(x)) 159 155 160 156 found, sigma0, sigma1, sigma2, k = \ 161 _search_triangles_of_vertices(triangles, x)157 root._search_triangles_of_vertices(triangles, x) 162 158 163 159 if k >= 0: … … 175 171 [10, 3]]: 176 172 177 triangles = _trilist_from_data(mesh,root.search(x))173 triangles = root._trilist_from_data(root.search(x)) 178 174 179 175 #print x, candidate_vertices 180 176 found, sigma0, sigma1, sigma2, k = \ 181 _search_triangles_of_vertices(triangles,177 root._search_triangles_of_vertices(triangles, 182 178 ensure_numeric(x)) 183 179 if k >= 0: … … 225 221 #x = [3.5, 1.5] 226 222 x = [2.5, 1.5] 227 element_found, sigma0, sigma1, sigma2, k = \ 228 search_tree_of_vertices(root, mesh, x) 223 element_found, sigma0, sigma1, sigma2, k = root.search_fast(x) 229 224 # One point 230 225 x = [3.00005, 2.999994] 231 element_found, sigma0, sigma1, sigma2, k = \ 232 search_tree_of_vertices(root, mesh, x) 226 element_found, sigma0, sigma1, sigma2, k = root.search_fast(x) 233 227 assert element_found is True 234 228 assert k == 1 -
anuga_core/source/anuga/geometry/mesh_quadtree.py
r7715 r7716 11 11 from anuga.utilities.numerical_tools import ensure_numeric 12 12 from anuga.config import max_float 13 from anuga.geometry.quad import Cell14 from a nuga.geometry.aabb import AABB13 from quad import Cell 14 from aabb import AABB 15 15 16 16 from anuga.utilities import compile … … 27 27 28 28 29 initial_search_value = 'uncomment search_functions code first'#030 search_one_cell_time = initial_search_value31 search_more_cells_time = initial_search_value32 33 29 # FIXME(Ole): Could we come up with a less confusing structure? 34 30 # FIXME(James): remove this global var … … 41 37 num.array([-1.1,-1.1])))]] 42 38 43 last_triangle = LAST_TRIANGLE 44 45 def search_tree_of_vertices(root, mesh, x): 39 40 41 class MeshQuadtree(Cell): 42 """ A quadtree constructed from the given mesh. 43 This class is the root node of a quadtree, 44 and derives from a Cell. 45 It contains optimisations and search patterns specific to meshes. 46 46 """ 47 Find the triangle (element) that the point x is in. 48 49 Inputs: 50 root: A quad tree of the vertices 51 mesh: The mesh which the quad tree indexes into 52 x: The point being placed 53 54 Return: 55 element_found, sigma0, sigma1, sigma2, k 56 57 where 58 element_found: True if a triangle containing x was found 59 sigma0, sigma1, sigma2: The interpolated values 60 k: Index of triangle (if found) 61 62 """ 63 global search_one_cell_time 64 global search_more_cells_time 65 66 if last_triangle[0][1] != -10: 67 # check the last triangle found first 47 def __init__(self, mesh): 48 """Build quad tree for mesh. 49 50 All vertex indices in the mesh are stored in a quadtree. 51 """ 52 53 extents = AABB(*mesh.get_extent(absolute=True)) 54 extents.grow(1.001) # To avoid round off error 55 Cell.__init__(self, extents, None) # root has no parent 56 57 N = len(mesh) 58 self.mesh = mesh 59 self.last_triangle = LAST_TRIANGLE 60 61 # Get x,y coordinates for all vertices for all triangles 62 V = mesh.get_vertex_coordinates(absolute=True) 63 64 # Check each triangle 65 for i in range(N): 66 x0, y0 = V[3*i, :] 67 x1, y1 = V[3*i+1, :] 68 x2, y2 = V[3*i+2, :] 69 70 # insert a tuple with an AABB, and the triangle index as data 71 self._insert((AABB(min([x0, x1, x2]), max([x0, x1, x2]), \ 72 min([y0, y1, y2]), max([y0, y1, y2])), \ 73 i)) 74 75 def search_fast(self, x): 76 """ 77 Find the triangle (element) that the point x is in. 78 79 Inputs: 80 root: A quad tree of the vertices 81 mesh: The mesh which the quad tree indexes into 82 x: The point being placed 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 if self.last_triangle[0][1] != -10: 94 # check the last triangle found first 95 element_found, sigma0, sigma1, sigma2, k = \ 96 self._search_triangles_of_vertices(self.last_triangle, x) 97 98 if element_found: 99 return True, sigma0, sigma1, sigma2, k 100 101 branch = self.last_triangle[0][1] 102 103 if branch == -10: 104 branch = self 105 106 # test neighbouring tris 107 tri_data = branch.test_leaves(x) 108 triangles = self._trilist_from_data(tri_data) 68 109 element_found, sigma0, sigma1, sigma2, k = \ 69 _search_triangles_of_vertices(last_triangle, x) 70 110 self._search_triangles_of_vertices(triangles, x) 111 if element_found: 112 return True, sigma0, sigma1, sigma2, k 113 114 # search to bottom of tree from last found leaf 115 tri_data = branch.search(x) 116 triangles = self._trilist_from_data(tri_data) 117 element_found, sigma0, sigma1, sigma2, k = \ 118 self._search_triangles_of_vertices(triangles, x) 71 119 if element_found: 72 120 return True, sigma0, sigma1, sigma2, k 73 121 74 branch = last_triangle[0][1] 122 # search rest of tree 123 element_found = False 124 while branch: 125 if not branch.parent: 126 # search from top of tree if we are at root 127 siblings = [self] 128 else: 129 siblings = branch.get_siblings() 130 131 for sibling in siblings: 132 tri_data = sibling.search(x) 133 triangles = self._trilist_from_data(tri_data) 134 element_found, sigma0, sigma1, sigma2, k = \ 135 self._search_triangles_of_vertices(triangles, x) 136 if element_found: 137 return True, sigma0, sigma1, sigma2, k 138 139 branch = branch.parent 140 if branch: 141 tri_data = branch.test_leaves(x) 142 triangles = self._trilist_from_data(tri_data) 143 element_found, sigma0, sigma1, sigma2, k = \ 144 self._search_triangles_of_vertices(triangles, x) 145 if element_found: 146 return True, sigma0, sigma1, sigma2, k 147 148 return element_found, sigma0, sigma1, sigma2, k 149 150 151 def _search_triangles_of_vertices(self, triangles, x): 152 """Search for triangle containing x amongs candidate_vertices in triangles 153 154 This is called by search_tree_of_vertices once the appropriate node 155 has been found from the quad tree. 156 157 158 This function is responsible for most of the compute time in 159 fit and interpolate. 160 """ 161 x = ensure_numeric(x, num.float) 162 163 # These statments are needed if triangles is empty 164 sigma2 = -10.0 165 sigma0 = -10.0 166 sigma1 = -10.0 167 k = -10 168 169 # For all vertices in same cell as point x 170 element_found = False 171 for k, node, tri_verts_norms in triangles: 172 tri = tri_verts_norms[0] 173 tri = ensure_numeric(tri) 174 # k is the triangle index 175 # tri is a list of verts (x, y), representing a triangle 176 # Find triangle that contains x (if any) and interpolate 177 178 # Input check disabled to speed things up. 179 if bool(_is_inside_triangle(x, tri, int(True), 1.0e-12, 1.0e-12)): 180 181 n0, n1, n2 = tri_verts_norms[1] 182 sigma0, sigma1, sigma2 =\ 183 compute_interpolation_values(tri, n0, n1, n2, x) 184 185 element_found = True 186 187 # Don't look for any other triangles in the triangle list 188 self.last_triangle = [[k, node, tri_verts_norms]] 189 break 190 191 return element_found, sigma0, sigma1, sigma2, k 192 193 194 def _trilist_from_data(self, indices): 195 """return a list of lists. For the inner lists, 196 The first element is the triangle index, 197 the second element is a list.for this list 198 the first element is a list of three (x, y) vertices, 199 the following elements are the three triangle normals. 200 201 """ 202 203 ret_list = [] 204 for i, node in indices: 205 vertices = self.mesh.get_vertex_coordinates(triangle_id=i, absolute=True) 206 n0 = self.mesh.get_normal(i, 0) 207 n1 = self.mesh.get_normal(i, 1) 208 n2 = self.mesh.get_normal(i, 2) 209 ret_list.append([i, node, [vertices, (n0, n1, n2)]]) 210 return ret_list 211 212 def set_last_triangle(self): 213 self.last_triangle = LAST_TRIANGLE 214 75 215 76 if branch == -10:77 branch = root78 79 # test neighbouring tris80 tri_data = branch.test_leaves(x)81 triangles = _trilist_from_data(mesh, tri_data)82 element_found, sigma0, sigma1, sigma2, k = \83 _search_triangles_of_vertices(triangles, x)84 if element_found:85 return True, sigma0, sigma1, sigma2, k86 87 # search to bottom of tree from last found leaf88 tri_data = branch.search(x)89 triangles = _trilist_from_data(mesh, tri_data)90 element_found, sigma0, sigma1, sigma2, k = \91 _search_triangles_of_vertices(triangles, x)92 if element_found:93 return True, sigma0, sigma1, sigma2, k94 95 # search rest of tree96 element_found = False97 while branch:98 if not branch.parent:99 # search from top of tree if we are at root100 siblings = [root]101 else:102 siblings = branch.get_siblings()103 104 for sibling in siblings:105 tri_data = sibling.search(x)106 triangles = _trilist_from_data(mesh, tri_data)107 element_found, sigma0, sigma1, sigma2, k = \108 _search_triangles_of_vertices(triangles, x)109 if element_found:110 return True, sigma0, sigma1, sigma2, k111 112 branch = branch.parent113 if branch:114 tri_data = branch.test_leaves(x)115 triangles = _trilist_from_data(mesh, tri_data)116 element_found, sigma0, sigma1, sigma2, k = \117 _search_triangles_of_vertices(triangles, x)118 if element_found:119 return True, sigma0, sigma1, sigma2, k120 121 return element_found, sigma0, sigma1, sigma2, k122 123 124 def _search_triangles_of_vertices(triangles, x):125 """Search for triangle containing x amongs candidate_vertices in triangles126 127 This is called by search_tree_of_vertices once the appropriate node128 has been found from the quad tree.129 130 131 This function is responsible for most of the compute time in132 fit and interpolate.133 """134 global last_triangle135 136 x = ensure_numeric(x, num.float)137 138 # These statments are needed if triangles is empty139 sigma2 = -10.0140 sigma0 = -10.0141 sigma1 = -10.0142 k = -10143 144 # For all vertices in same cell as point x145 element_found = False146 for k, node, tri_verts_norms in triangles:147 tri = tri_verts_norms[0]148 tri = ensure_numeric(tri)149 # k is the triangle index150 # tri is a list of verts (x, y), representing a triangle151 # Find triangle that contains x (if any) and interpolate152 153 # Input check disabled to speed things up.154 if bool(_is_inside_triangle(x, tri, int(True), 1.0e-12, 1.0e-12)):155 156 n0, n1, n2 = tri_verts_norms[1]157 sigma0, sigma1, sigma2 =\158 compute_interpolation_values(tri, n0, n1, n2, x)159 160 element_found = True161 162 # Don't look for any other triangles in the triangle list163 last_triangle = [[k, node, tri_verts_norms]]164 break165 166 return element_found, sigma0, sigma1, sigma2, k167 168 169 def _trilist_from_data(mesh, indices):170 """return a list of lists. For the inner lists,171 The first element is the triangle index,172 the second element is a list.for this list173 the first element is a list of three (x, y) vertices,174 the following elements are the three triangle normals.175 176 """177 178 ret_list = []179 for i, node in indices:180 vertices = mesh.get_vertex_coordinates(triangle_id=i, absolute=True)181 n0 = mesh.get_normal(i, 0)182 n1 = mesh.get_normal(i, 1)183 n2 = mesh.get_normal(i, 2)184 ret_list.append([i, node, [vertices, (n0, n1, n2)]])185 return ret_list186 187 188 216 def compute_interpolation_values(triangle, n0, n1, n2, x): 189 217 """Compute linear interpolation of point x and triangle. … … 200 228 201 229 return sigma0, sigma1, sigma2 202 203 def set_last_triangle(): 204 global last_triangle 205 last_triangle = LAST_TRIANGLE 206 207 def search_times(): 208 209 global search_one_cell_time 210 global search_more_cells_time 211 212 return search_one_cell_time, search_more_cells_time 213 214 def reset_search_times(): 215 216 global search_one_cell_time 217 global search_more_cells_time 218 search_one_cell_time = initial_search_value 219 search_more_cells_time = initial_search_value 220 221 222 class MeshQuadtree(Cell): 223 """ A quadtree constructed from the given mesh. 224 This class is actually the root node of a tree, 225 and derives from a Cell. 226 """ 227 def __init__(self, mesh): 228 """Build quad tree for mesh. 229 230 All vertex indices in the mesh are stored in a quadtree. 231 """ 232 233 extents = AABB(*mesh.get_extent(absolute=True)) 234 extents.grow(1.001) # To avoid round off error 235 Cell.__init__(self, extents, None) # root has no parent 236 237 N = len(mesh) 238 239 # Get x,y coordinates for all vertices for all triangles 240 V = mesh.get_vertex_coordinates(absolute=True) 241 242 # Check each triangle 243 for i in range(N): 244 x0, y0 = V[3*i, :] 245 x1, y1 = V[3*i+1, :] 246 x2, y2 = V[3*i+2, :] 247 248 # insert a tuple with an AABB, and the triangle index as data 249 self._insert((AABB(min([x0, x1, x2]), max([x0, x1, x2]), \ 250 min([y0, y1, y2]), max([y0, y1, y2])), \ 251 i)) 252 230 -
anuga_core/source/anuga/geometry/test_geometry.py
r7714 r7716 75 75 (AABB(7, 8, 3, 4), 333), (AABB(1, 10, 0, 1), 444)]) 76 76 77 result = cell.search( x = 8.5, y = 1.5)77 result = cell.search([8.5, 1.5]) 78 78 assert type(result) in [types.ListType,types.TupleType],\ 79 79 'should be a list' -
anuga_core/source/anuga/geometry/test_meshquad.py
r7715 r7716 2 2 import numpy as num 3 3 4 from a nuga.geometry.aabb import AABB5 from anuga.geometry.quad import Cell6 from anuga.fit_interpolate.mesh_quadtree import MeshQuadtree4 from aabb import AABB 5 from quad import Cell 6 from mesh_quadtree import MeshQuadtree 7 7 from anuga.abstract_2d_finite_volumes.general_mesh import General_mesh as Mesh 8 8
Note: See TracChangeset
for help on using the changeset viewer.
