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

Last change on this file since 6545 was 6545, checked in by ole, 15 years ago

Removed more cluttter in search functions.

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