# 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
9
10from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
11from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular
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
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
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
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
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
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
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.