source: inundation/ga/storm_surge/pyvolution/test_least_squares.py @ 626

Last change on this file since 626 was 611, checked in by duncan, 20 years ago

Optimised the least_squares algorithm for building A matrix
File size: 22.5 KB
Line 
1#!/usr/bin/env python
2
3#TEST
4
5import unittest
6from math import sqrt
7
8
9from least_squares import *
10from Numeric import allclose, array, transpose
11
12def distance(x, y):
13    return sqrt( sum( (array(x)-array(y))**2 ))         
14
15def linear_function(point):
16    point = array(point)
17    return point[:,0]+point[:,1]
18   
19       
20class TestCase(unittest.TestCase):
21
22    def setUp(self):
23        pass
24       
25    def tearDown(self):
26        pass
27
28
29    def test_datapoint_at_centroid(self):
30        a = [0.0, 0.0]
31        b = [0.0, 2.0]
32        c = [2.0,0.0]
33        points = [a, b, c]
34        vertices = [ [1,0,2] ]   #bac
35
36        data = [ [2.0/3, 2.0/3] ] #Use centroid as one data point       
37       
38        interp = Interpolation(points, vertices, data)
39        assert allclose(interp.get_A(), [[1./3, 1./3, 1./3]]) 
40       
41
42    def test_quad_tree(self):
43        p0 = [-10.0, -10.0]
44        p1 = [20.0, -10.0]
45        p2 = [-10.0, 20.0]
46        p3 = [10.0, 50.0]
47        p4 = [30.0, 30.0]
48        p5 = [50.0, 10.0]
49        p6 = [40.0, 60.0]
50        p7 = [60.0, 40.0]
51        p8 = [-66.0, 20.0]
52        p9 = [10.0, -66.0]
53       
54        points = [p0, p1, p2, p3, p4, p5, p6, p7, p8, p9]
55        triangles = [ [0, 1, 2],
56                      [3, 2, 4],
57                      [4, 2, 1],
58                      [4, 1, 5],
59                      [3, 4, 6],
60                      [6, 4, 7],
61                      [7, 4, 5],
62                      [8, 0, 2],
63                      [0, 9, 1]]   
64
65        data = [ [4,4] ] 
66       
67        interp = Interpolation(points, triangles, data, alpha = 0.0)
68        #print "PDSG - interp.get_A()", interp.get_A()
69        answer =  [ [ 0.06666667,  0.46666667,  0.46666667,  0.,
70                      0., 0. , 0., 0., 0., 0.]]
71        assert allclose(interp.get_A(), answer)
72        interp.set_point_coordinates([[-30, -30]]) #point outside of mesh
73        #print "PDSG - interp.get_A()", interp.get_A()
74        answer =  [ [ 0.0,  0.0,  0.0,  0.,
75                      0., 0. , 0., 0., 0., 0.]]
76        assert allclose(interp.get_A(), answer)
77
78       
79        #point outside of quad tree root cell
80        interp.set_point_coordinates([[-70, -70]])
81        #print "PDSG - interp.get_A()", interp.get_A()
82        answer =  [ [ 0.0,  0.0,  0.0,  0.,
83                      0., 0. , 0., 0., 0., 0.]]
84        assert allclose(interp.get_A(), answer)
85       
86
87    def test_datapoints_at_vertices(self):
88        """Test that data points coinciding with vertices yield a diagonal matrix
89        """
90       
91        a = [0.0, 0.0]
92        b = [0.0, 2.0]
93        c = [2.0,0.0]
94        points = [a, b, c]
95        vertices = [ [1,0,2] ]   #bac
96
97        data = points #Use data at vertices       
98       
99        interp = Interpolation(points, vertices, data)
100        assert allclose(interp.get_A(), [[1., 0., 0.],
101                                   [0., 1., 0.],
102                                   [0., 0., 1.]]) 
103       
104
105
106    def test_datapoints_on_edge_midpoints(self):
107        """Try datapoints midway on edges -
108        each point should affect two matrix entries equally
109        """
110       
111        a = [0.0, 0.0]
112        b = [0.0, 2.0]
113        c = [2.0,0.0]
114        points = [a, b, c]
115        vertices = [ [1,0,2] ]   #bac
116
117        data = [ [0., 1.], [1., 0.], [1., 1.] ]
118       
119        interp = Interpolation(points, vertices, data)
120       
121        assert allclose(interp.get_A(), [[0.5, 0.5, 0.0],  #Affects vertex 1 and 0   
122                                   [0.5, 0.0, 0.5],  #Affects vertex 0 and 2
123                                   [0.0, 0.5, 0.5]]) #Affects vertex 1 and 2
124
125
126    def test_datapoints_on_edges(self):
127        """Try datapoints on edges -
128        each point should affect two matrix entries in proportion
129        """
130       
131        a = [0.0, 0.0]
132        b = [0.0, 2.0]
133        c = [2.0,0.0]
134        points = [a, b, c]
135        vertices = [ [1,0,2] ]   #bac
136
137        data = [ [0., 1.5], [1.5, 0.], [1.5, 0.5] ]
138       
139        interp = Interpolation(points, vertices, data)
140       
141        assert allclose(interp.get_A(), [[0.25, 0.75, 0.0],  #Affects vertex 1 and 0
142                                   [0.25, 0.0, 0.75],  #Affects vertex 0 and 2
143                                   [0.0, 0.25, 0.75]]) #Affects vertex 1 and 2
144
145    def test_arbitrary_datapoints(self):
146        """Try arbitrary datapoints
147        """
148
149        from Numeric import sum
150       
151        a = [0.0, 0.0]
152        b = [0.0, 2.0]
153        c = [2.0,0.0]
154        points = [a, b, c]
155        vertices = [ [1,0,2] ]   #bac
156
157        data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44] ]
158       
159        interp = Interpolation(points, vertices, data)
160
161        assert allclose(sum(interp.get_A(), axis=1), 1.0)
162       
163
164    # this causes a memory error in scipy.sparce       
165    def test_more_triangles(self):
166       
167        a = [-1.0, 0.0]
168        b = [3.0, 4.0]
169        c = [4.0,1.0]
170        d = [-3.0, 2.0] #3
171        e = [-1.0,-2.0]
172        f = [1.0, -2.0] #5
173
174        points = [a, b, c, d,e,f]
175        triangles = [[0,1,3],[1,0,2],[0,4,5], [0,5,2]] #abd bac aef afc
176
177        #Data points
178        data_points = [ [-3., 2.0], [-2, 1], [0.0, 1], [0, 3], [2, 3], [-1.0/3,-4./3] ]
179        interp = Interpolation(points, triangles, data_points) 
180
181        answer = [[0.0, 0.0, 0.0, 1.0, 0.0, 0.0],    #Affects point d     
182                  [0.5, 0.0, 0.0, 0.5, 0.0, 0.0],    #Affects points a and d
183                  [0.75, 0.25, 0.0, 0.0, 0.0, 0.0],  #Affects points a and b
184                  [0.0, 0.5, 0.0, 0.5, 0.0, 0.0],    #Affects points a and d
185                  [0.25, 0.75, 0.0, 0.0, 0.0, 0.0],  #Affects points a and b
186                  [1./3, 0.0, 0.0, 0.0, 1./3, 1./3]] #Affects points a, e and f
187
188
189        #FIXME: This test failed on our Linux box, but not on Windows.
190        A = interp.get_A()
191        for i in range(A.shape[0]):
192            for j in range(A.shape[1]):
193                if not allclose(A[i,j], answer[i][j]):
194                    print i,j,':',A[i,j], answer[i][j] 
195
196
197        assert allclose(interp.get_A(), answer)
198
199
200
201
202    def test_smooth_attributes_to_mesh(self):
203        a = [0.0, 0.0]
204        b = [0.0, 5.0]
205        c = [5.0, 0.0]
206        points = [a, b, c]
207        triangles = [ [1,0,2] ]   #bac
208         
209        d1 = [1.0, 1.0]
210        d2 = [1.0, 3.0]
211        d3 = [3.0,1.0]
212        z1 = 2
213        z2 = 4
214        z3 = 4
215        data_coords = [d1, d2, d3] 
216       
217        interp = Interpolation(points, triangles, data_coords, alpha=5.0e-20)
218        z = [z1, z2, z3]
219        f = interp.fit(z)
220        answer = [0, 5., 5.]
221
222        #print "f\n",f
223        #print "answer\n",answer
224       
225        assert allclose(f, answer, atol=1e-7)
226       
227       
228    def test_smooth_att_to_meshII(self):
229       
230        a = [0.0, 0.0]
231        b = [0.0, 5.0]
232        c = [5.0, 0.0]
233        points = [a, b, c]
234        triangles = [ [1,0,2] ]   #bac
235         
236        d1 = [1.0, 1.0]
237        d2 = [1.0, 2.0]
238        d3 = [3.0,1.0]
239        data_coords = [d1, d2, d3]       
240        z = linear_function(data_coords)
241        interp = Interpolation(points, triangles, data_coords, alpha=0.0)
242        f = interp.fit(z)
243        answer = linear_function(points)
244       
245        assert allclose(f, answer)
246   
247    def test_smooth_attributes_to_meshIII(self):
248       
249        a = [-1.0, 0.0]
250        b = [3.0, 4.0]
251        c = [4.0,1.0]
252        d = [-3.0, 2.0] #3
253        e = [-1.0,-2.0]
254        f = [1.0, -2.0] #5
255
256        vertices = [a, b, c, d,e,f]
257        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
258
259        point_coords = [[-2.0, 2.0],
260                        [-1.0, 1.0],
261                        [0.0,2.0],
262                        [1.0, 1.0],
263                        [2.0, 1.0],
264                        [0.0,0.0],
265                        [1.0, 0.0],
266                        [0.0, -1.0],
267                        [-0.2,-0.5],
268                        [-0.9, -1.5],
269                        [0.5, -1.9],
270                        [3.0,1.0]]
271
272        z = linear_function(point_coords)
273        interp = Interpolation(vertices, triangles, point_coords, alpha=0.0)
274
275        #print 'z',z
276        f = interp.fit(z)
277        answer = linear_function(vertices)
278        #print "f\n",f
279        #print "answer\n",answer
280        assert allclose(f, answer)
281       
282
283    def test_smooth_attributes_to_meshIV(self):
284        """ Testing 2 attributes smoothed to the mesh
285        """
286       
287        a = [0.0, 0.0]
288        b = [0.0, 5.0]
289        c = [5.0, 0.0]
290        points = [a, b, c]
291        triangles = [ [1,0,2] ]   #bac
292         
293        d1 = [1.0, 1.0]
294        d2 = [1.0, 3.0]
295        d3 = [3.0, 1.0]
296        z1 = [2, 4]
297        z2 = [4, 8]
298        z3 = [4, 8]
299        data_coords = [d1, d2, d3] 
300       
301        interp = Interpolation(points, triangles, data_coords, alpha=0.0)
302        z = [z1, z2, z3]
303        f =  interp.fit_points(z)
304        answer = [[0,0], [5., 10.], [5., 10.]]
305        assert allclose(f, answer)
306       
307    def test_interpolate_attributes_to_points(self):
308        v0 = [0.0, 0.0]
309        v1 = [0.0, 5.0]
310        v2 = [5.0, 0.0]
311       
312        vertices = [v0, v1, v2]
313        triangles = [ [1,0,2] ]   #bac
314         
315        d0 = [1.0, 1.0]
316        d1 = [1.0, 2.0]
317        d2 = [3.0, 1.0]
318        point_coords = [ d0, d1, d2]
319       
320        interp = Interpolation(vertices, triangles, point_coords)
321        f = linear_function(vertices)
322        z = interp.interpolate(f)
323        answer = linear_function(point_coords)
324       
325
326        assert allclose(z, answer)
327
328       
329    def test_interpolate_attributes_to_pointsII(self):
330        a = [-1.0, 0.0]
331        b = [3.0, 4.0]
332        c = [4.0, 1.0]
333        d = [-3.0, 2.0] #3
334        e = [-1.0, -2.0]
335        f = [1.0, -2.0] #5
336
337        vertices = [a, b, c, d,e,f]
338        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
339
340         
341        point_coords = [[-2.0, 2.0],
342                        [-1.0, 1.0],
343                        [0.0, 2.0],
344                        [1.0, 1.0],
345                        [2.0, 1.0],
346                        [0.0, 0.0],
347                        [1.0, 0.0],
348                        [0.0, -1.0],
349                        [-0.2, -0.5],
350                        [-0.9, -1.5],
351                        [0.5, -1.9],
352                        [3.0, 1.0]]
353       
354        interp = Interpolation(vertices, triangles, point_coords)
355        f = linear_function(vertices)
356        z = interp.interpolate(f)
357        answer = linear_function(point_coords)
358        #print "z",z
359        #print "answer",answer
360        assert allclose(z, answer)
361   
362    def test_interpolate_attributes_to_pointsIII(self):
363        v0 = [0.0, 0.0]
364        v1 = [0.0, 5.0]
365        v2 = [5.0, 0.0]
366       
367        vertices = [v0, v1, v2]
368        triangles = [ [1,0,2] ]   #bac
369         
370        d0 = [1.0, 1.0]
371        d1 = [1.0, 2.0]
372        d2 = [3.0,1.0]
373        point_coords = [ d0, d1, d2]
374       
375        interp = Interpolation(vertices, triangles, point_coords)
376        f = [ [0., 0.],  [5., 10.],  [5., 10.]] #linear_function(vertices)
377        #print "f",f
378       
379        z = interp.interpolate(f)
380        answer = [ [2., 4.],  [3., 6.],  [4., 8.]]
381
382        #print "***********"
383        #print "z",z
384        #print "answer",answer
385        #print "***********"
386
387        assert allclose(z, answer)
388   
389    def test_interpolate_attributes_to_pointsIV(self):
390        a = [-1.0, 0.0]
391        b = [3.0, 4.0]
392        c = [4.0, 1.0]
393        d = [-3.0, 2.0] #3
394        e = [-1.0, -2.0]
395        f = [1.0, -2.0] #5
396
397        vertices = [a, b, c, d,e,f]
398        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
399
400         
401        point_coords = [[-2.0, 2.0],
402                        [-1.0, 1.0],
403                        [0.0, 2.0],
404                        [1.0, 1.0],
405                        [2.0, 1.0],
406                        [0.0, 0.0],
407                        [1.0, 0.0],
408                        [0.0, -1.0],
409                        [-0.2, -0.5],
410                        [-0.9, -1.5],
411                        [0.5, -1.9],
412                        [3.0, 1.0]]
413       
414        interp = Interpolation(vertices, triangles, point_coords)
415        f = array([linear_function(vertices),2*linear_function(vertices) ])
416        f = transpose(f)
417        #print "f",f
418        z = interp.interpolate(f)
419        answer = [linear_function(point_coords),
420             
421     2*linear_function(point_coords) ]
422        answer = transpose(answer)
423        #print "z",z
424        #print "answer",answer
425        assert allclose(z, answer)
426       
427    def test_smooth_attributes_to_mesh_function(self):
428        """ Testing 2 attributes smoothed to the mesh
429        """
430        """Test multiple attributes
431        """
432       
433        a = [0.0, 0.0]
434        b = [0.0, 5.0]
435        c = [5.0, 0.0]
436        points = [a, b, c]
437        triangles = [ [1,0,2] ]   #bac
438         
439        d1 = [1.0, 1.0]
440        d2 = [1.0, 3.0]
441        d3 = [3.0, 1.0]
442        z1 = [2, 4]
443        z2 = [4, 8]
444        z3 = [4, 8]
445        data_coords = [d1, d2, d3]
446        z = [z1, z2, z3]
447       
448        f = fit_to_mesh(points, triangles, data_coords, z, alpha=0.0)
449        answer = [[0, 0], [5., 10.], [5., 10.]]
450       
451        assert allclose(f, answer)
452
453
454
455    #Tests of smoothing matrix   
456    def test_smoothing_matrix_one_triangle(self):
457        from Numeric import dot
458        a = [0.0, 0.0]
459        b = [0.0, 2.0]
460        c = [2.0,0.0]
461        points = [a, b, c]
462       
463        vertices = [ [1,0,2] ]   #bac
464
465        interp = Interpolation(points, vertices)
466
467        assert allclose(interp.get_D(), [[1, -0.5, -0.5],
468                                   [-0.5, 0.5, 0],
469                                   [-0.5, 0, 0.5]])
470
471        #Define f(x,y) = x
472        f = array([0,0,2]) #Value at global vertex 2
473
474        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
475        #           int 1 dx dy = area = 2
476        assert dot(dot(f, interp.get_D()), f) == 2
477       
478        #Define f(x,y) = y
479        f = array([0,2,0])  #Value at global vertex 1
480
481        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
482        #           int 1 dx dy = area = 2
483        assert dot(dot(f, interp.get_D()), f) == 2
484
485        #Define f(x,y) = x+y
486        f = array([0,2,2])  #Values at global vertex 1 and 2
487
488        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
489        #           int 2 dx dy = 2*area = 4
490        assert dot(dot(f, interp.get_D()), f) == 4       
491       
492
493
494    def test_smoothing_matrix_more_triangles(self):
495        from Numeric import dot
496       
497        a = [0.0, 0.0]
498        b = [0.0, 2.0]
499        c = [2.0,0.0]
500        d = [0.0, 4.0]
501        e = [2.0, 2.0]
502        f = [4.0,0.0]
503
504        points = [a, b, c, d, e, f]
505        #bac, bce, ecf, dbe, daf, dae
506        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
507
508        interp = Interpolation(points, vertices)
509
510
511        #assert allclose(interp.get_D(), [[1, -0.5, -0.5],
512        #                           [-0.5, 0.5, 0],
513        #                           [-0.5, 0, 0.5]])
514
515        #Define f(x,y) = x
516        f = array([0,0,2,0,2,4]) #f evaluated at points a-f
517
518        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
519        #           int 1 dx dy = total area = 8
520        assert dot(dot(f, interp.get_D()), f) == 8
521       
522        #Define f(x,y) = y
523        f = array([0,2,0,4,2,0]) #f evaluated at points a-f
524
525        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
526        #           int 1 dx dy = area = 8
527        assert dot(dot(f, interp.get_D()), f) == 8
528
529        #Define f(x,y) = x+y
530        f = array([0,2,2,4,4,4])  #f evaluated at points a-f   
531
532        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
533        #           int 2 dx dy = 2*area = 16
534        assert dot(dot(f, interp.get_D()), f) == 16       
535
536
537    def test_fit_and_interpolation(self):
538        from mesh import Mesh
539       
540        a = [0.0, 0.0]
541        b = [0.0, 2.0]
542        c = [2.0, 0.0]
543        d = [0.0, 4.0]
544        e = [2.0, 2.0]
545        f = [4.0, 0.0]
546
547        points = [a, b, c, d, e, f]
548        #bac, bce, ecf, dbe, daf, dae
549        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
550
551        #Get (enough) datapoints
552        data_points = [[ 0.66666667, 0.66666667],
553                       [ 1.33333333, 1.33333333],
554                       [ 2.66666667, 0.66666667],
555                       [ 0.66666667, 2.66666667],
556                       [ 0.0, 1.0],
557                       [ 0.0, 3.0],
558                       [ 1.0, 0.0],
559                       [ 1.0, 1.0],
560                       [ 1.0, 2.0],
561                       [ 1.0, 3.0],                                         
562                       [ 2.0, 1.0],
563                       [ 3.0, 0.0],
564                       [ 3.0, 1.0]]                                         
565       
566        interp = Interpolation(points, triangles, data_points, alpha=0.0)
567
568        z = linear_function(data_points)
569        answer = linear_function(points)
570
571        f = interp.fit(z)
572
573        #print "f",f
574        #print "answer",answer
575        assert allclose(f, answer)
576
577        #Map back
578        z1 = interp.interpolate(f)
579        #print "z1\n", z1
580        #print "z\n",z
581        assert allclose(z, z1)
582
583
584    def test_smoothing_and_interpolation(self):
585       
586        a = [0.0, 0.0]
587        b = [0.0, 2.0]
588        c = [2.0, 0.0]
589        d = [0.0, 4.0]
590        e = [2.0, 2.0]
591        f = [4.0, 0.0]
592
593        points = [a, b, c, d, e, f]
594        #bac, bce, ecf, dbe, daf, dae
595        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
596
597        #Get (too few!) datapoints
598        data_points = [[ 0.66666667, 0.66666667],
599                       [ 1.33333333, 1.33333333],
600                       [ 2.66666667, 0.66666667],
601                       [ 0.66666667, 2.66666667]]
602
603        z = linear_function(data_points)
604        answer = linear_function(points)       
605
606        #Make interpolator with too few data points and no smoothing
607        interp = Interpolation(points, triangles, data_points, alpha=0.0)
608        #Must raise an exception       
609        try:
610            f = interp.fit(z)
611        except:
612            pass
613
614        #Now try with smoothing parameter
615        interp = Interpolation(points, triangles, data_points, alpha=1.0e-13)
616       
617        f = interp.fit(z)
618        #f will be different from answerr due to smoothing
619        assert allclose(f, answer,atol=5)
620
621        #Map back
622        z1 = interp.interpolate(f)
623        assert allclose(z, z1)
624
625    def test_fit_to_mesh_file(self):
626        from load_mesh.loadASCII import mesh_file_to_mesh_dictionary, \
627             export_trianglulation_file   
628        import tempfile
629        import os
630       
631        # create a .tsh file, no user outline
632        mesh_dic = {}
633        mesh_dic['generatedpointlist'] = [[0.0, 0.0],
634                                          [0.0, 5.0],
635                                          [5.0, 0.0]]
636        mesh_dic['generatedtrianglelist'] =  [[0, 2, 1]]
637        mesh_dic['generatedsegmentlist'] = [[0, 1], [2, 0], [1, 2]]
638        mesh_dic['generatedtriangleattributelist'] = [['']]
639        mesh_dic['generatedpointattributelist'] = [[], [], []]
640        mesh_dic['generatedpointattributetitlelist'] = []
641        mesh_dic['generatedtriangleneighborlist'] = [[-1, -1, -1]]
642        mesh_dic['generatedsegmentmarkerlist'] = ['external',
643                                                  'external',
644                                                  'external']       
645        mesh_file = tempfile.mktemp(".tsh")
646        export_trianglulation_file(mesh_file,mesh_dic)
647       
648        # create an .xya file
649        point_file = tempfile.mktemp(".xya")
650        fd = open(point_file,'w')
651        fd.write("elevation, stage \n 1.0, 1.0,2.,4 \n 1.0, 3.0,4,8 \n 3.0,1.0,4.,8 \n")
652        fd.close()
653       
654        mesh_output_file = "new_trianlge.tsh"
655        fit_to_mesh_file(mesh_file,
656                         point_file,
657                         mesh_output_file,
658                         alpha = 0.0) 
659        # load in the .tsh file we just wrote
660        mesh_dic = mesh_file_to_mesh_dictionary(mesh_output_file)
661
662        assert allclose(mesh_dic['generatedpointattributelist'],
663                        [[0.0, 0.0],
664                         [5.0, 10.0],
665                         [5.0,10.0]])
666       
667        self.failUnless(mesh_dic['generatedpointattributetitlelist']  ==
668                        ['elevation','stage'],
669                        'test_fit_to_mesh_file failed')
670       
671        #clean up
672        os.remove(mesh_file)
673        os.remove(point_file)
674       
675    def test_fit_to_mesh_fileII(self):
676        from load_mesh.loadASCII import mesh_file_to_mesh_dictionary, \
677             export_trianglulation_file   
678        import tempfile
679        import os
680       
681        # create a .tsh file, no user outline
682        mesh_dic = {}
683        mesh_dic['generatedpointlist'] = [[0.0, 0.0],
684                                          [0.0, 5.0],
685                                          [5.0, 0.0]]
686        mesh_dic['generatedtrianglelist'] =  [[0, 2, 1]]
687        mesh_dic['generatedsegmentlist'] = [[0, 1], [2, 0], [1, 2]]
688        mesh_dic['generatedtriangleattributelist'] = [['']]
689        mesh_dic['generatedpointattributelist'] = [[1,2], [1,2], [1,2]]
690        mesh_dic['generatedpointattributetitlelist'] = ['density', 'temp']
691        mesh_dic['generatedtriangleneighborlist'] = [[-1, -1, -1]]
692        mesh_dic['generatedsegmentmarkerlist'] = ['external',
693                                                  'external',
694                                                  'external']       
695        mesh_file = tempfile.mktemp(".tsh")
696        export_trianglulation_file(mesh_file,mesh_dic)
697       
698        # create an .xya file
699        point_file = tempfile.mktemp(".xya")
700        fd = open(point_file,'w')
701        fd.write("elevation, stage \n 1.0, 1.0,2.,4 \n 1.0, 3.0,4,8 \n 3.0,1.0,4.,8 \n")
702        fd.close()
703       
704        mesh_output_file = "new_triangle.tsh"
705        fit_to_mesh_file(mesh_file,
706                         point_file,
707                         mesh_output_file,
708                         alpha = 0.0) 
709        # load in the .tsh file we just wrote
710        mesh_dic = mesh_file_to_mesh_dictionary(mesh_output_file)
711
712        assert allclose(mesh_dic['generatedpointattributelist'],
713                        [[1.0, 2.0,0.0, 0.0],
714                         [1.0, 2.0,5.0, 10.0],
715                         [1.0, 2.0,5.0,10.0]])
716       
717        self.failUnless(mesh_dic['generatedpointattributetitlelist']  ==
718                        ['density', 'temp','elevation','stage'],
719                        'test_fit_to_mesh_file failed')
720       
721        #clean up
722        os.remove(mesh_file)
723        os.remove(mesh_output_file)       
724        os.remove(point_file)
725       
726#-------------------------------------------------------------
727if __name__ == "__main__":
728    suite = unittest.makeSuite(TestCase,'test')
729    runner = unittest.TextTestRunner(verbosity=1)
730    runner.run(suite)
731
732   
733   
734
735
736
Note: See TracBrowser for help on using the repository browser.