source: anuga_core/source/anuga/fit_interpolate/test_search_functions.py @ 7711

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

Refactored geometry classes to live in their own folder.

File size: 8.5 KB
Line 
1#!/usr/bin/env python
2
3
4import unittest
5from search_functions import search_tree_of_vertices, set_last_triangle
6from search_functions import _search_triangles_of_vertices
7from search_functions import _trilist_from_indices
8from search_functions import compute_interpolation_values, MeshQuadtree
9
10from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
11from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular
12from anuga.geometry.quad import Cell
13from anuga.geometry.aabb import AABB
14from anuga.utilities.numerical_tools import ensure_numeric
15from anuga.geometry.polygon import is_inside_polygon, is_inside_triangle   
16
17import numpy as num
18
19class Test_search_functions(unittest.TestCase):
20    def setUp(self):
21        pass
22
23    def tearDown(self):
24        pass
25
26    def NOtest_that_C_extension_compiles(self):
27        FN = 'search_functions_ext.c'
28        try:
29            import search_functions_ext
30        except:
31            from compile import compile
32
33            try:
34                compile(FN)
35            except:
36                raise 'Could not compile %s' %FN
37            else:
38                import search_functions_ext
39
40
41
42    def test_off_and_boundary(self):
43        """test_off: Test a point off the mesh
44        """
45
46        points, vertices, boundary = rectangular(1, 1, 1, 1)
47        mesh = Mesh(points, vertices, boundary)
48
49        #Test that points are arranged in a counter clock wise order
50        mesh.check_integrity()
51
52        root = MeshQuadtree(mesh)
53        set_last_triangle()
54
55        found, s0, s1, s2, k = search_tree_of_vertices(root, mesh, [-0.2, 10.7])
56        assert found is False
57
58        found, s0, s1, s2, k = search_tree_of_vertices(root, mesh, [0, 0])
59        assert found is True
60       
61       
62    def test_small(self):
63        """test_small: Two triangles
64        """
65
66        points, vertices, boundary = rectangular(1, 1, 1, 1)
67        mesh = Mesh(points, vertices, boundary)
68
69        #Test that points are arranged in a counter clock wise order
70        mesh.check_integrity()
71
72        root = MeshQuadtree(mesh)
73        set_last_triangle()
74
75        x = [0.2, 0.7]
76        found, s0, s1, s2, k = search_tree_of_vertices(root, mesh, x)
77        assert k == 1 # Triangle one
78        assert found is True       
79       
80    def test_bigger(self):
81        """test_bigger
82       
83        test larger mesh
84        """
85
86        points, vertices, boundary = rectangular(4, 4, 1, 1)
87        mesh = Mesh(points, vertices, boundary)
88
89        #Test that points are arranged in a counter clock wise order
90        mesh.check_integrity()
91
92        root = MeshQuadtree(mesh)
93        set_last_triangle()
94
95        for x in [[0.6, 0.3], [0.1, 0.2], [0.7,0.7],
96                  [0.1,0.9], [0.4,0.6], [0.9,0.1],
97                  [10, 3]]:
98           
99            found, s0, s1, s2, k = search_tree_of_vertices(root, mesh,
100                                                           ensure_numeric(x))                                   
101                                                           
102            if k >= 0:
103                V = mesh.get_vertex_coordinates(k) # nodes for triangle k
104                assert is_inside_polygon(x, V)
105                assert found is True
106                #print k, x
107            else:
108                assert found is False               
109
110       
111
112    def test_large(self):
113        """test_larger mesh and different quad trees
114        """
115
116        points, vertices, boundary = rectangular(10, 12, 1, 1)
117        mesh = Mesh(points, vertices, boundary)
118
119        #Test that points are arranged in a counter clock wise order
120        mesh.check_integrity()
121
122       
123
124        root = MeshQuadtree(mesh)
125        set_last_triangle()
126        #print m, root.show()
127
128        for x in [[0.6, 0.3], [0.1, 0.2], [0.7,0.7],
129                  [0.1,0.9], [0.4,0.6], [0.9,0.1],
130                  [10, 3]]:
131           
132            found, s0, s1, s2, k = search_tree_of_vertices(root, mesh, x)
133
134            if k >= 0:
135                V = mesh.get_vertex_coordinates(k) # nodes for triangle k
136                assert is_inside_triangle(x, V, closed=True)
137                assert is_inside_polygon(x, V)
138                assert found is True
139            else:
140                assert found is False               
141
142       
143            if k == 0: return   
144
145    def test_underlying_function(self):
146        """test_larger mesh and different quad trees
147        """
148
149        points, vertices, boundary = rectangular(2, 2, 1, 1)
150        mesh = Mesh(points, vertices, boundary)
151
152        root = MeshQuadtree(mesh)
153        set_last_triangle()
154
155        # One point
156        x = ensure_numeric([0.5, 0.5])
157
158        triangles = _trilist_from_indices(mesh, root.search(x[0], x[1]))
159   
160        found, sigma0, sigma1, sigma2, k = \
161               _search_triangles_of_vertices(triangles, x)
162
163        if k >= 0:
164            V = mesh.get_vertex_coordinates(k) # nodes for triangle k
165            assert is_inside_polygon(x, V)
166            assert found is True
167        else:
168            assert found is False               
169
170       
171
172        # More points   
173        for x in [[0.6, 0.3], [0.1, 0.2], [0.7,0.7],
174                  [0.1,0.9], [0.4,0.6], [0.9,0.1],
175                  [10, 3]]:
176               
177            triangles = _trilist_from_indices(mesh, root.search(x[0], x[1]))
178
179            #print x, candidate_vertices
180            found, sigma0, sigma1, sigma2, k = \
181                   _search_triangles_of_vertices(triangles,
182                                                 ensure_numeric(x))
183            if k >= 0:
184                V = mesh.get_vertex_coordinates(k) # nodes for triangle k
185                assert is_inside_polygon(x, V)
186                assert found is True
187            else:
188                assert found is False
189
190               
191
192    def expanding_search(self):
193        """test_larger mesh and different quad trees
194        """
195       
196        p0 = [2,1]
197        p1 = [4,1]
198        p2 = [4.,4]
199        p3 = [2,4]
200        p4 = [5,4]
201
202        p5 = [-1,-1]
203        p6 = [1,-1]
204        p7 = [1,1]
205        p8 = [-1,1]
206
207        points = [p0,p1,p2, p3,p4,p5,p6,p7,p8]
208        #
209        vertices = [[0,1,2],[0,2,3],[1,4,2],[5,6,7], [5,7,8]]
210        mesh = Mesh(points, vertices)
211
212        # Don't do this, want to control the max and mins
213        #root = build_quadtree(mesh, max_points_per_cell=4)
214   
215
216        root = Cell(-3, 9, -3, 9,
217                    max_points_per_cell = 4)
218        #Insert indices of all vertices
219        root.insert( range(mesh.number_of_nodes) )
220
221        #Build quad tree and return
222        root.split()
223       
224        # One point
225        #x = [3.5, 1.5]
226        x = [2.5, 1.5]
227        element_found, sigma0, sigma1, sigma2, k = \
228                       search_tree_of_vertices(root, mesh, x)
229        # One point
230        x = [3.00005, 2.999994]
231        element_found, sigma0, sigma1, sigma2, k = \
232                       search_tree_of_vertices(root, mesh, x)
233        assert element_found is True
234        assert k == 1
235       
236       
237    def test_compute_interpolation_values(self):
238        """test_compute_interpolation_values
239       
240        Test that interpolation values are correc
241       
242        This test used to check element_found as output from
243        find_triangle_compute_interpolation(triangle, n0, n1, n2, x)
244        and that failed before 18th March 2009.
245       
246        Now this function no longer returns this flag, so the test
247        is merely checknig the sigmas.
248       
249        """
250       
251        triangle = num.array([[306951.77151059, 6194462.14596986],
252                              [306952.58403545, 6194459.65001246],
253                              [306953.55109034, 6194462.0041216]])
254
255
256        n0 = [0.92499377, -0.37998227]
257        n1 = [0.07945684,  0.99683831]
258        n2 = [-0.95088404, -0.30954732]
259       
260        x = [306953.344, 6194461.5]
261       
262        # Test that point is indeed inside triangle
263        assert is_inside_polygon(x, triangle, 
264                                 closed=True, verbose=False)
265        assert is_inside_triangle(x, triangle, 
266                                  closed=True, verbose=False)                                 
267       
268        sigma0, sigma1, sigma2 = \
269            compute_interpolation_values(triangle, n0, n1, n2, x)
270           
271        msg = 'Point which is clearly inside triangle was not found'
272        assert abs(sigma0 + sigma1 + sigma2 - 1) < 1.0e-6, msg
273
274#-------------------------------------------------------------
275if __name__ == "__main__":
276    suite = unittest.makeSuite(Test_search_functions, 'test')
277    runner = unittest.TextTestRunner(verbosity=1)
278    runner.run(suite)
279   
Note: See TracBrowser for help on using the repository browser.