Changeset 7718
- Timestamp:
- May 11, 2010, 8:12:13 PM (15 years ago)
- Location:
- anuga_core/source/anuga
- Files:
-
- 6 edited
-
fit_interpolate/benchmark_least_squares.py (modified) (2 diffs)
-
fit_interpolate/test_search_functions.py (modified) (2 diffs)
-
geometry/mesh_quadtree.py (modified) (9 diffs)
-
geometry/quad.py (modified) (5 diffs)
-
geometry/test_meshquad.py (modified) (1 diff)
-
shallow_water/test_shallow_water_domain.py (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/fit_interpolate/benchmark_least_squares.py
r7317 r7718 23 23 import profile , pstats 24 24 from math import sqrt 25 26 from anuga.fit_interpolate.search_functions import search_times, \27 reset_search_times28 25 29 26 from anuga.fit_interpolate.interpolate import Interpolate … … 262 259 # backing up. 263 260 264 search_one_cell_time, search_more_cells_time = search_times()265 reset_search_times()261 #search_one_cell_time, search_more_cells_time = search_times() 262 #reset_search_times() 266 263 #print "bench - search_one_cell_time",search_one_cell_time 267 264 #print "bench - search_more_cells_time", search_more_cells_time -
anuga_core/source/anuga/fit_interpolate/test_search_functions.py
r7716 r7718 3 3 4 4 import unittest 5 from mesh_quadtree import MeshQuadtree , compute_interpolation_values5 from mesh_quadtree import MeshQuadtree 6 6 7 7 from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh … … 228 228 assert k == 1 229 229 230 231 def test_compute_interpolation_values(self):232 """test_compute_interpolation_values233 234 Test that interpolation values are correc235 236 This test used to check element_found as output from237 find_triangle_compute_interpolation(triangle, n0, n1, n2, x)238 and that failed before 18th March 2009.239 240 Now this function no longer returns this flag, so the test241 is merely checknig the sigmas.242 243 """244 245 triangle = num.array([[306951.77151059, 6194462.14596986],246 [306952.58403545, 6194459.65001246],247 [306953.55109034, 6194462.0041216]])248 249 250 n0 = [0.92499377, -0.37998227]251 n1 = [0.07945684, 0.99683831]252 n2 = [-0.95088404, -0.30954732]253 254 x = [306953.344, 6194461.5]255 256 # Test that point is indeed inside triangle257 assert is_inside_polygon(x, triangle,258 closed=True, verbose=False)259 assert is_inside_triangle(x, triangle,260 closed=True, verbose=False)261 262 sigma0, sigma1, sigma2 = \263 compute_interpolation_values(triangle, n0, n1, n2, x)264 265 msg = 'Point which is clearly inside triangle was not found'266 assert abs(sigma0 + sigma1 + sigma2 - 1) < 1.0e-6, msg267 230 268 231 #------------------------------------------------------------- -
anuga_core/source/anuga/geometry/mesh_quadtree.py
r7716 r7718 27 27 28 28 29 # FIXME(Ole): Could we come up with a less confusing structure?30 # FIXME(James): remove this global var31 29 LAST_TRIANGLE = [[-10, -10, 32 30 (num.array([[max_float, max_float], … … 57 55 N = len(mesh) 58 56 self.mesh = mesh 59 self. last_triangle = LAST_TRIANGLE57 self.set_last_triangle() 60 58 61 59 # Get x,y coordinates for all vertices for all triangles … … 91 89 92 90 """ 91 92 x = ensure_numeric(x, num.float) 93 93 94 if self.last_triangle[0][1] != -10: 94 95 # check the last triangle found first 95 96 element_found, sigma0, sigma1, sigma2, k = \ 96 97 self._search_triangles_of_vertices(self.last_triangle, x) 97 98 98 if element_found: 99 99 return True, sigma0, sigma1, sigma2, k 100 100 101 101 branch = self.last_triangle[0][1] 102 103 if branch == -10:104 branch = self105 102 106 103 # test neighbouring tris … … 112 109 return True, sigma0, sigma1, sigma2, k 113 110 114 # search to bottom of tree from last found leaf115 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)119 if element_found:120 return True, sigma0, sigma1, sigma2, k121 122 111 # search rest of tree 123 112 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: 113 next_search = [branch] 114 while branch: 115 for sibling in next_search: 132 116 tri_data = sibling.search(x) 133 117 triangles = self._trilist_from_data(tri_data) … … 136 120 if element_found: 137 121 return True, sigma0, sigma1, sigma2, k 138 122 123 next_search = branch.get_siblings() 139 124 branch = branch.parent 140 125 if branch: … … 144 129 self._search_triangles_of_vertices(triangles, x) 145 130 if element_found: 146 return True, sigma0, sigma1, sigma2, k 131 return True, sigma0, sigma1, sigma2, k 147 132 148 133 return element_found, sigma0, sigma1, sigma2, k … … 158 143 This function is responsible for most of the compute time in 159 144 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 145 """ 146 171 147 for k, node, tri_verts_norms in triangles: 172 148 tri = tri_verts_norms[0] … … 178 154 # Input check disabled to speed things up. 179 155 if bool(_is_inside_triangle(x, tri, int(True), 1.0e-12, 1.0e-12)): 180 181 156 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 157 xi0, xi1, xi2 = tri 158 159 sigma0 = num.dot((x-xi1), n0)/num.dot((xi0-xi1), n0) 160 sigma1 = num.dot((x-xi2), n1)/num.dot((xi1-xi2), n1) 161 sigma2 = num.dot((x-xi0), n2)/num.dot((xi2-xi0), n2) 186 162 187 163 # Don't look for any other triangles in the triangle list 188 164 self.last_triangle = [[k, node, tri_verts_norms]] 189 break 190 191 return element_found, sigma0, sigma1, sigma2, k 165 return True, sigma0, sigma1, sigma2, k 166 return False, -1, -1, -1, -10 192 167 193 168 194 169 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 list198 the first element is a list of three (x, y) vertices,199 the following elements are the three triangle normals.200 201 """202 203 170 ret_list = [] 204 171 for i, node in indices: … … 212 179 def set_last_triangle(self): 213 180 self.last_triangle = LAST_TRIANGLE 181 self.last_triangle[0][1] = self # point at root by default 214 182 215 183 216 def compute_interpolation_values(triangle, n0, n1, n2, x):217 """Compute linear interpolation of point x and triangle.218 219 n0, n1, n2 are normal to the tree edges.220 """221 184 222 # Get the three vertex_points of candidate triangle k223 xi0, xi1, xi2 = triangle224 225 sigma0 = num.dot((x-xi1), n0)/num.dot((xi0-xi1), n0)226 sigma1 = num.dot((x-xi2), n1)/num.dot((xi1-xi2), n1)227 sigma2 = num.dot((x-xi0), n2)/num.dot((xi2-xi0), n2)228 229 return sigma0, sigma1, sigma2230 185 -
anuga_core/source/anuga/geometry/quad.py
r7714 r7718 19 19 """ 20 20 21 def __init__(self, extents, parent, depth = 0,21 def __init__(self, extents, parent, 22 22 name = 'cell'): 23 23 … … 27 27 self.extents = extents 28 28 self.parent = parent 29 self.searched = [-10, -15]30 self.depth = depth31 29 32 30 # The points in this cell … … 36 34 37 35 def __repr__(self): 38 str = ' (%d)%s: leaves: %d' \39 % (self. depth, self.name , len(self.leaves))36 str = '%s: leaves: %d' \ 37 % (self.name , len(self.leaves)) 40 38 if self.children: 41 39 str += ', children: %d' % (len(self.children)) … … 84 82 # try splitting this cell and see if we get a trivial in 85 83 subregion1, subregion2 = self.extents.split() 86 if subregion1.is_trivial_in(new_region): 87 self.children = [Cell(subregion1, self, depth=self.depth+1), Cell(subregion2, self, depth=self.depth+1)] 88 self.children[0]._insert(new_leaf) 89 return 90 elif subregion2.is_trivial_in(new_region): 91 self.children = [Cell(subregion1, self, depth=self.depth+1), Cell(subregion2, self, depth=self.depth+1)] 92 self.children[1]._insert(new_leaf) 93 return 84 85 # option 1 - try splitting 4 ways 86 subregion11, subregion12 = subregion1.split() 87 subregion21, subregion22 = subregion2.split() 88 regions = [subregion11, subregion12, subregion21, subregion22] 89 for region in regions: 90 if region.is_trivial_in(new_region): 91 self.children = [Cell(x, parent=self) for x in regions] 92 self._insert(new_leaf) 93 return 94 95 # option 2 - try splitting 2 ways 96 #if subregion1.is_trivial_in(new_region): 97 #self.children = [Cell(subregion1, self), Cell(subregion2, self)] 98 #self.children[0]._insert(new_leaf) 99 #return 100 #elif subregion2.is_trivial_in(new_region): 101 #self.children = [Cell(subregion1, self), Cell(subregion2, self)] 102 #self.children[1]._insert(new_leaf) 103 #return 94 104 95 105 # recursion ended without finding a fit, so attach it as a leaf … … 128 138 if depth == 0: 129 139 log.critical() 130 print '%s%s' % (' '*depth, self.name), self.extents,' [', self.leaves, ']'140 print '%s%s' % (' '*depth, self.name), self.extents,' [', self.leaves, ']' 131 141 if self.children: 132 142 log.critical() 133 143 for child in self.children: 134 144 child.show(depth+1) 135 136 145 137 146 def search(self, x): 138 147 """return a list of possible intersections with geometry""" 139 148 140 149 intersecting_regions = self.test_leaves(x) 141 150 -
anuga_core/source/anuga/geometry/test_meshquad.py
r7716 r7718 145 145 self.assertEqual(idx, 2) 146 146 147 def NOtest_ get_parent():148 """ Get a child's parent.147 def NOtest_num_visits(self): 148 """ Test optimisation code. 149 149 """ 150 cell = Cell(AABB(0, 6, 0, 6), 'cell') 150 a = [3, 7] 151 b = [5, 7] 152 c = [5, 5] 153 d = [7, 7] 154 e = [15, 15] 155 f = [15, 30] 156 g = [30, 10] 157 h = [30, 30] 151 158 152 p0 = [2,1] 153 p1 = [4,1] 154 p2 = [4,4] 155 p3 = [2,4] 156 p4 = [5,4] 157 158 points = [p0,p1,p2, p3, p4] 159 vertices = [[0,1,2],[0,2,3],[1,4,2]] 159 points = [a, b, c, d, e, f, g, h] 160 161 #bac, bce, ecf, dbe, daf, dae 162 vertices = [[1,0,2], [1,3,4], [1,2,3], [5,4,7], [4,6,7]] 160 163 161 164 mesh = Mesh(points, vertices) 162 165 163 Q = MeshQuadtree(mesh) 164 results = Q.search([4.5, 3]) 165 assert len(results) == 1 166 self.assertEqual(results[0], 2) 167 results = Q.search([5,4.]) 168 self.assertEqual(len(results),1) 169 self.assertEqual(results[0], 2) 166 Q = MeshQuadtree(mesh) 167 168 169 results = Q.search_fast([5.5, 5.5]) 170 print 'visits: ', Q.count_visits() 171 172 Q.clear_visits() 173 results = Q.search_fast([30, 10]) 174 print 'visits: ', Q.count_visits() 175 176 print 'second time:' 177 178 Q.clear_visits() 179 results = Q.search_fast([5.5, 5.5]) 180 print 'visits: ', Q.count_visits() 170 181 ################################################################################ 171 182 -
anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py
r7712 r7718 1812 1812 gauge_values[i].append(w[i]) 1813 1813 1814 # Reference (nautilus 26/6/2008)1815 1816 G0 = [-0.20000000000000001, -0.20000000000000001, -0. 20000000000000001, -0.19729755264559332, -0.18143006781428656, -0.17264603126739803, -0.16490767093739186, -0.15580226478343825, -0.15341783464373857, -0.14986976789238646, -0.1483294354194323, -0.19385112169384708, -0.19907249891161938, -0.19911053885072424, -0.19913592111915288, -0.1991524040791994, -0.19916445999802629, -0.19917790655716, -0.19919280664859376, -0.19920527965237472, -0.19921381524768964, -0.1992178338453264, -0.19921637002152079, -0.19920783637167225, -0.19921088191840772, -0.19919037392917549, -0.19918812396862254, -0.19918939056535054, -0.19919336800395021, -0.19919832209289323, -0.19920329732763048, -0.19920824490004829, -0.19921301385072837, -0.19921748460975705, -0.1992215655528298, -0.19922516461276893, -0.19922820250665138, -0.19923065047412888, -0.19923254166903823, -0.1992339741831877, -0.19923510970454786, -0.19923612584283332, -0.19923716297385921, -0.199238301846201, -0.19923956605537632, -0.19924092368339749, -0.19924231511786381, -0.19924367609907154, -0.19924495328100339, -0.19924611576332413, -0.19924716014245455]1814 # Captured data from code manually inspected for correctness 11/5/2010 1815 1816 G0 = [-0.20000000000000001, -0.20000000000000001, -0.19999999639817381, -0.19900000000000001, -0.19787180620524364, -0.18775357859677891, -0.17717419377267535, -0.17038758570131016, -0.1632648895200452, -0.16084615293651969, -0.15796677906889614, -0.15610187077034218, -0.19890877268012952, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19767101056475031, -0.19664020763833073, -0.19625551309686509, -0.19622800460530576, -0.19670420796840654, -0.19757556964453959, -0.19853453054821305, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001] 1817 1817 G1 = [-0.29999999999999993, -0.29999999999999993, -0.29312012778165086, -0.26736324399674027, -0.23794021958953263, -0.21878933349240881, -0.2063216250007715, -0.19630760483269205, -0.18778429435995911, -0.18241377240055098, -0.17860448861463396, -0.17352695896444634, -0.16449192186501546, -0.18953645096562932, -0.20557835913635514, -0.21071182844586539, -0.21191011610076083, -0.21068538990761979, -0.20737339293005266, -0.20407505794124087, -0.20190406545920248, -0.19987295968833538, -0.19832486432402194, -0.19774281454937118, -0.19763203737080062, -0.19769585527535841, -0.19795493497088903, -0.19841633252351859, -0.19886909990019369, -0.19940203310127189, -0.19990812587783294, -0.20016182994265733, -0.20026318580583596, -0.20022855664242351, -0.20010463117095048, -0.199930824934102, -0.19973147245180795, -0.19957882377401115, -0.19948188719957555, -0.19944190810764481, -0.19945196574982568, -0.19949690579414533, -0.19955788251561279, -0.19961580710612042, -0.19965570581001513, -0.19966924072268322, -0.19966164565962002, -0.1996429044375452, -0.19962419292036049, -0.19961203640143024, -0.19960927744740292] 1818 1818 G2 = [-0.40000000000000002, -0.39653752812451043, -0.32890250043784469, -0.29967800724628818, -0.27572700190320054, -0.25893877811152488, -0.2422572785756886, -0.22758771043596362, -0.21433402198450011, -0.20266876405969875, -0.19393588052745447, -0.18766846928826031, -0.18182722392707523, -0.16865402349317951, -0.18110962350013771, -0.19547215504627991, -0.20229670797148233, -0.20494160396870145, -0.20614855426298931, -0.20616358034144364, -0.2050749944506553, -0.20329491765401383, -0.20174840479413605, -0.20040155802388543, -0.19945041001985703, -0.19898209259951644, -0.19868944730468457, -0.19849574784901269, -0.19855116942355155, -0.19885365879108613, -0.19925381355621719, -0.19971339601933866, -0.20000203652933318, -0.20011337898513715, -0.2001188815225943, -0.20007116187977508, -0.20000646342088527, -0.19993034594397555, -0.19985512250322437, -0.19978530424723981, -0.19973800656091345, -0.19971854873923284, -0.19972327598147538, -0.19974189986453444, -0.19976101408429672, -0.19977705650936498, -0.19978787136068635, -0.19979206946039513, -0.19978865745584601, -0.19978181317853855, -0.19977875136946421]
Note: See TracChangeset
for help on using the changeset viewer.
