source: anuga_core/source/anuga/fit_interpolate/test_fit.py @ 3514

Last change on this file since 3514 was 3514, checked in by duncan, 18 years ago

Hi all,
I'm doing a change in the anuga structure, moving the code to

\anuga_core\source\anuga

After you have done an svn update, the PYTHONPATH has to be changed to;
PYTHONPATH = anuga_core/source/

This is part of changes required to make installation of anuga quicker and reducing the size of our sandpits.

If any imports are broken, try fixing them. With adding anuga. to them for example. If this seems to have really broken things, email/phone me.

Cheers
Duncan

File size: 26.8 KB
Line 
1#!/usr/bin/env python
2
3#TEST
4import sys
5import unittest
6from math import sqrt
7
8from Numeric import zeros, take, compress, Float, Int, dot, concatenate, \
9     ArrayType, allclose, array
10
11from fit import *
12from anuga.pyvolution.neighbour_mesh import Mesh
13from anuga.utilities.sparse import Sparse, Sparse_CSR
14from anuga.coordinate_transforms.geo_reference import Geo_reference
15from anuga.utilities.numerical_tools import ensure_numeric
16from anuga.geospatial_data.geospatial_data import Geospatial_data
17
18def distance(x, y):
19    return sqrt( sum( (array(x)-array(y))**2 ))
20
21def linear_function(point):
22    point = array(point)
23    return point[:,0]+point[:,1]
24
25
26class Test_Fit(unittest.TestCase):
27
28    def setUp(self):
29        pass
30
31    def tearDown(self):
32        pass
33
34
35    def test_smooth_attributes_to_mesh(self):
36        a = [0.0, 0.0]
37        b = [0.0, 5.0]
38        c = [5.0, 0.0]
39        points = [a, b, c]
40        triangles = [ [1,0,2] ]   #bac
41
42        d1 = [1.0, 1.0]
43        d2 = [1.0, 3.0]
44        d3 = [3.0,1.0]
45        z1 = 2
46        z2 = 4
47        z3 = 4
48        data_coords = [d1, d2, d3]
49
50        z = [z1, z2, z3]
51        fit = Fit(points, triangles, alpha=0)
52        #print "interp.get_A()", interp.get_A()
53        fit._build_matrix_AtA_Atz(ensure_numeric(data_coords),
54                                  ensure_numeric(z))
55        #print "Atz - from fit", fit.Atz
56        #print "AtA - from fit", fit.AtA.todense()
57        #print "z",z
58
59        assert allclose(fit.Atz, [2.8, 3.6, 3.6], atol=1e-7)
60
61        f = fit.fit()
62       
63        answer = [0, 5., 5.]
64
65        #print "f\n",f
66        #print "answer\n",answer
67
68        assert allclose(f, answer, atol=1e-7)
69
70    def test_smooth_att_to_meshII(self):
71
72        a = [0.0, 0.0]
73        b = [0.0, 5.0]
74        c = [5.0, 0.0]
75        points = [a, b, c]
76        triangles = [ [1,0,2] ]   #bac
77
78        d1 = [1.0, 1.0]
79        d2 = [1.0, 2.0]
80        d3 = [3.0,1.0]
81        data_coords = [d1, d2, d3]
82        z = linear_function(data_coords)
83        #print "z",z
84
85        interp = Fit(points, triangles, alpha=0.0)
86        f = interp.fit(data_coords, z)
87        answer = linear_function(points)
88        #print "f\n",f
89        #print "answer\n",answer
90
91        assert allclose(f, answer)
92
93    def test_smooth_attributes_to_meshIII(self):
94
95        a = [-1.0, 0.0]
96        b = [3.0, 4.0]
97        c = [4.0,1.0]
98        d = [-3.0, 2.0] #3
99        e = [-1.0,-2.0]
100        f = [1.0, -2.0] #5
101
102        vertices = [a, b, c, d,e,f]
103        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
104
105        point_coords = [[-2.0, 2.0],
106                        [-1.0, 1.0],
107                        [0.0,2.0],
108                        [1.0, 1.0],
109                        [2.0, 1.0],
110                        [0.0,0.0],
111                        [1.0, 0.0],
112                        [0.0, -1.0],
113                        [-0.2,-0.5],
114                        [-0.9, -1.5],
115                        [0.5, -1.9],
116                        [3.0,1.0]]
117
118        z = linear_function(point_coords)
119        interp = Fit(vertices, triangles,
120                                alpha=0.0)
121
122        #print 'z',z
123        f = interp.fit(point_coords,z)
124        answer = linear_function(vertices)
125        #print "f\n",f
126        #print "answer\n",answer
127        assert allclose(f, answer)
128
129   
130    def test_smooth_attributes_to_meshIV(self):
131        """ Testing 2 attributes smoothed to the mesh
132        """
133
134        a = [0.0, 0.0]
135        b = [0.0, 5.0]
136        c = [5.0, 0.0]
137        points = [a, b, c]
138        triangles = [ [1,0,2] ]   #bac
139
140        d1 = [1.0, 1.0]
141        d2 = [1.0, 3.0]
142        d3 = [3.0, 1.0]
143        z1 = [2, 4]
144        z2 = [4, 8]
145        z3 = [4, 8]
146        data_coords = [d1, d2, d3]
147
148        z = [z1, z2, z3]
149        fit = Fit(points, triangles, alpha=0)
150       
151        f =  fit.fit(data_coords,z)
152        answer = [[0,0], [5., 10.], [5., 10.]]
153        assert allclose(f, answer)
154
155    def test_smooth_attributes_to_mesh_build_fit_subset(self):
156
157        a = [-1.0, 0.0]
158        b = [3.0, 4.0]
159        c = [4.0,1.0]
160        d = [-3.0, 2.0] #3
161        e = [-1.0,-2.0]
162        f = [1.0, -2.0] #5
163
164        vertices = [a, b, c, d,e,f]
165        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
166
167        interp = Fit(vertices, triangles,
168                                alpha=0.0)
169
170        point_coords = [[-2.0, 2.0],
171                        [-1.0, 1.0],
172                        [0.0,2.0],
173                        [1.0, 1.0],
174                       ]
175
176        z = linear_function(point_coords)
177
178        f = interp.build_fit_subset(point_coords,z)
179       
180        point_coords = [
181                        [2.0, 1.0],
182                        [0.0,0.0],
183                        [1.0, 0.0],
184                        [0.0, -1.0],
185                        [-0.2,-0.5],
186                        [-0.9, -1.5],
187                        [0.5, -1.9],
188                        [3.0,1.0]]
189       
190        z = linear_function(point_coords)
191
192        f = interp.build_fit_subset(point_coords,z)
193       
194        #print 'z',z
195        f = interp.fit()
196        answer = linear_function(vertices)
197        #print "f\n",f
198        #print "answer\n",answer
199        assert allclose(f, answer)
200
201    def test_fit_and_interpolation(self):
202
203        a = [0.0, 0.0]
204        b = [0.0, 2.0]
205        c = [2.0, 0.0]
206        d = [0.0, 4.0]
207        e = [2.0, 2.0]
208        f = [4.0, 0.0]
209
210        points = [a, b, c, d, e, f]
211        #bac, bce, ecf, dbe, daf, dae
212        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
213
214        #Get (enough) datapoints
215        data_points = [[ 0.66666667, 0.66666667],
216                       [ 1.33333333, 1.33333333],
217                       [ 2.66666667, 0.66666667],
218                       [ 0.66666667, 2.66666667],
219                       [ 0.0, 1.0],
220                       [ 0.0, 3.0],
221                       [ 1.0, 0.0],
222                       [ 1.0, 1.0],
223                       [ 1.0, 2.0],
224                       [ 1.0, 3.0],
225                       [ 2.0, 1.0],
226                       [ 3.0, 0.0],
227                       [ 3.0, 1.0]]
228
229        z = linear_function(data_points)
230        interp = Fit(points, triangles,
231                                alpha=0.0)
232
233        answer = linear_function(points)
234
235        f = interp.fit(data_points, z)
236
237        #print "f",f
238        #print "answer",answer
239        assert allclose(f, answer)
240
241       
242    def test_smoothing_and_interpolation(self):
243
244        a = [0.0, 0.0]
245        b = [0.0, 2.0]
246        c = [2.0, 0.0]
247        d = [0.0, 4.0]
248        e = [2.0, 2.0]
249        f = [4.0, 0.0]
250
251        points = [a, b, c, d, e, f]
252        #bac, bce, ecf, dbe, daf, dae
253        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
254
255        #Get (too few!) datapoints
256        data_points = [[ 0.66666667, 0.66666667],
257                       [ 1.33333333, 1.33333333],
258                       [ 2.66666667, 0.66666667],
259                       [ 0.66666667, 2.66666667]]
260
261        z = linear_function(data_points)
262        answer = linear_function(points)
263
264        #Make interpolator with too few data points and no smoothing
265        interp = Fit(points, triangles, alpha=0.0)
266        #Must raise an exception
267        try:
268            f = interp.fit(data_points,z)
269        except ToFewPointsError:
270            pass
271
272        #Now try with smoothing parameter
273        interp = Fit(points, triangles, alpha=1.0e-13)
274
275        f = interp.fit(data_points,z)
276        #f will be different from answer due to smoothing
277        assert allclose(f, answer,atol=5)
278
279
280    #Tests of smoothing matrix
281    def test_smoothing_matrix_one_triangle(self):
282        from Numeric import dot
283        a = [0.0, 0.0]
284        b = [0.0, 2.0]
285        c = [2.0,0.0]
286        points = [a, b, c]
287
288        vertices = [ [1,0,2] ]   #bac
289
290        interp = Fit(points, vertices)
291
292        assert allclose(interp.get_D(), [[1, -0.5, -0.5],
293                                   [-0.5, 0.5, 0],
294                                   [-0.5, 0, 0.5]])
295
296        #Define f(x,y) = x
297        f = array([0,0,2]) #Value at global vertex 2
298
299        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
300        #           int 1 dx dy = area = 2
301        assert dot(dot(f, interp.get_D()), f) == 2
302
303        #Define f(x,y) = y
304        f = array([0,2,0])  #Value at global vertex 1
305
306        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
307        #           int 1 dx dy = area = 2
308        assert dot(dot(f, interp.get_D()), f) == 2
309
310        #Define f(x,y) = x+y
311        f = array([0,2,2])  #Values at global vertex 1 and 2
312
313        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
314        #           int 2 dx dy = 2*area = 4
315        assert dot(dot(f, interp.get_D()), f) == 4
316
317
318    def test_smoothing_matrix_more_triangles(self):
319        from Numeric import dot
320
321        a = [0.0, 0.0]
322        b = [0.0, 2.0]
323        c = [2.0,0.0]
324        d = [0.0, 4.0]
325        e = [2.0, 2.0]
326        f = [4.0,0.0]
327
328        points = [a, b, c, d, e, f]
329        #bac, bce, ecf, dbe, daf, dae
330        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
331
332        interp = Fit(points, vertices)
333
334
335        #assert allclose(interp.get_D(), [[1, -0.5, -0.5],
336        #                           [-0.5, 0.5, 0],
337        #                           [-0.5, 0, 0.5]])
338
339        #Define f(x,y) = x
340        f = array([0,0,2,0,2,4]) #f evaluated at points a-f
341
342        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
343        #           int 1 dx dy = total area = 8
344        assert dot(dot(f, interp.get_D()), f) == 8
345
346        #Define f(x,y) = y
347        f = array([0,2,0,4,2,0]) #f evaluated at points a-f
348
349        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
350        #           int 1 dx dy = area = 8
351        assert dot(dot(f, interp.get_D()), f) == 8
352
353        #Define f(x,y) = x+y
354        f = array([0,2,2,4,4,4])  #f evaluated at points a-f
355
356        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
357        #           int 2 dx dy = 2*area = 16
358        assert dot(dot(f, interp.get_D()), f) == 16
359
360
361
362    def test_fit_and_interpolation_with_different_origins(self):
363        """Fit a surface to one set of points. Then interpolate that surface
364        using another set of points.
365        This test tests situtaion where points and mesh belong to a different
366        coordinate system as defined by origin.
367        """
368
369        #Setup mesh used to represent fitted function
370        a = [0.0, 0.0]
371        b = [0.0, 2.0]
372        c = [2.0, 0.0]
373        d = [0.0, 4.0]
374        e = [2.0, 2.0]
375        f = [4.0, 0.0]
376
377        points = [a, b, c, d, e, f]
378        #bac, bce, ecf, dbe, daf, dae
379        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
380
381        #Datapoints to fit from
382        data_points1 = [[ 0.66666667, 0.66666667],
383                        [ 1.33333333, 1.33333333],
384                        [ 2.66666667, 0.66666667],
385                        [ 0.66666667, 2.66666667],
386                        [ 0.0, 1.0],
387                        [ 0.0, 3.0],
388                        [ 1.0, 0.0],
389                        [ 1.0, 1.0],
390                        [ 1.0, 2.0],
391                        [ 1.0, 3.0],
392                        [ 2.0, 1.0],
393                        [ 3.0, 0.0],
394                        [ 3.0, 1.0]]
395
396
397        #First check that things are OK when using same origin
398        mesh_origin = (56, 290000, 618000) #zone, easting, northing
399        data_origin = (56, 290000, 618000) #zone, easting, northing
400
401
402        #Fit surface to mesh
403        interp = Fit(points, triangles, 
404                               alpha=0.0,
405                               mesh_origin = mesh_origin)
406
407        data_geo_spatial = Geospatial_data(data_points1,
408                         geo_reference = Geo_reference(56, 290000, 618000))
409        z = linear_function(data_points1) #Example z-values
410        f = interp.fit(data_geo_spatial, z)   #Fitted values at vertices
411
412        #Shift datapoints according to new origins
413        for k in range(len(data_points1)):
414            data_points1[k][0] += mesh_origin[1] - data_origin[1]
415            data_points1[k][1] += mesh_origin[2] - data_origin[2]
416
417
418
419        #Fit surface to mesh
420        interp = Fit(points, triangles, 
421                               alpha=0.0) #,
422                               # mesh_origin = mesh_origin)
423
424        f1 = interp.fit(data_points1,z) #Fitted values at vertices (using same z as before)
425
426        assert allclose(f,f1), 'Fit should have been unaltered'
427
428
429    def test_smooth_attributes_to_mesh_function(self):
430        """ Testing 2 attributes smoothed to the mesh
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    def test_fit_to_mesh_w_georef(self):
454        """Simple check that georef works at the fit_to_mesh level
455        """
456       
457        from anuga.coordinate_transforms.geo_reference import Geo_reference
458
459        #Mesh
460        vertex_coordinates = [[0.76, 0.76],
461                              [0.76, 5.76],
462                              [5.76, 0.76]]
463        triangles = [[0,2,1]]
464
465        mesh_geo = Geo_reference(56,-0.76,-0.76)
466        #print "mesh_geo.get_absolute(vertex_coordinates)", \
467         #     mesh_geo.get_absolute(vertex_coordinates)
468
469        #Data                       
470        data_points = [[ 201.0, 401.0],
471                       [ 201.0, 403.0],
472                       [ 203.0, 401.0]]
473
474        z = [2, 4, 4]
475       
476        data_geo = Geo_reference(56,-200,-400)
477
478        #print "data_geo.get_absolute(data_points)", \
479        #      data_geo.get_absolute(data_points)
480       
481        #Fit
482        zz = fit_to_mesh(vertex_coordinates, triangles, data_points, z,
483                         data_origin = data_geo.get_origin(),
484                         mesh_origin = mesh_geo.get_origin(),
485                         alpha = 0)
486        assert allclose( zz, [0,5,5] )
487
488
489    def test_fit_to_mesh_file(self):
490        from load_mesh.loadASCII import import_mesh_file, \
491             export_mesh_file
492        import tempfile
493        import os
494
495        # create a .tsh file, no user outline
496        mesh_dic = {}
497        mesh_dic['vertices'] = [[0.0, 0.0],
498                                          [0.0, 5.0],
499                                          [5.0, 0.0]]
500        mesh_dic['triangles'] =  [[0, 2, 1]]
501        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
502        mesh_dic['triangle_tags'] = ['']
503        mesh_dic['vertex_attributes'] = [[], [], []]
504        mesh_dic['vertiex_attribute_titles'] = []
505        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
506        mesh_dic['segment_tags'] = ['external',
507                                                  'external',
508                                                  'external']
509        mesh_file = tempfile.mktemp(".tsh")
510        export_mesh_file(mesh_file,mesh_dic)
511
512        # create an .xya file
513        point_file = tempfile.mktemp(".xya")
514        fd = open(point_file,'w')
515        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")
516        fd.close()
517
518        mesh_output_file = tempfile.mktemp(".tsh") 
519        fit_to_mesh_file(mesh_file,
520                         point_file,
521                         mesh_output_file,
522                         alpha = 0.0)
523        # load in the .tsh file we just wrote
524        mesh_dic = import_mesh_file(mesh_output_file)
525        #print "mesh_dic",mesh_dic
526        ans =[[0.0, 0.0],
527              [5.0, 10.0],
528              [5.0,10.0]]
529        assert allclose(mesh_dic['vertex_attributes'],ans)
530
531        self.failUnless(mesh_dic['vertex_attribute_titles']  ==
532                        ['elevation','stage'],
533                        'test_fit_to_mesh_file failed')
534
535        #clean up
536        os.remove(mesh_file)
537        os.remove(point_file)
538        os.remove(mesh_output_file)
539
540
541    def test_fit_to_mesh_file3(self):
542        from load_mesh.loadASCII import import_mesh_file, \
543             export_mesh_file
544        import tempfile
545        import os
546
547        # create a .tsh file, no user outline
548        mesh_dic = {}
549        mesh_dic['vertices'] = [[0.76, 0.76],
550                                          [0.76, 5.76],
551                                          [5.76, 0.76]]
552        mesh_dic['triangles'] =  [[0, 2, 1]]
553        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
554        mesh_dic['triangle_tags'] = ['']
555        mesh_dic['vertex_attributes'] = [[], [], []]
556        mesh_dic['vertiex_attribute_titles'] = []
557        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
558        mesh_dic['segment_tags'] = ['external',
559                                                  'external',
560                                                  'external']
561        mesh_dic['geo_reference'] = Geo_reference(56,-0.76,-0.76)
562        mesh_file = tempfile.mktemp(".tsh")
563        export_mesh_file(mesh_file,mesh_dic)
564
565        # create an .xya file
566        point_file = tempfile.mktemp(".xya")
567        fd = open(point_file,'w')
568        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")
569        fd.close()
570
571        mesh_output_file = tempfile.mktemp(".tsh")
572        fit_to_mesh_file(mesh_file,
573                         point_file,
574                         mesh_output_file,
575                         alpha = 0.0)
576        # load in the .tsh file we just wrote
577        mesh_dic = import_mesh_file(mesh_output_file)
578        #print "mesh_dic",mesh_dic
579        ans =[[0.0, 0.0],
580              [5.0, 10.0],
581              [5.0,10.0]]
582        assert allclose(mesh_dic['vertex_attributes'],ans)
583
584        self.failUnless(mesh_dic['vertex_attribute_titles']  ==
585                        ['elevation','stage'],
586                        'test_fit_to_mesh_file failed')
587
588        #clean up
589        os.remove(mesh_file)
590        os.remove(point_file)
591        os.remove(mesh_output_file)
592
593    def test_fit_to_mesh_file4(self):
594        from load_mesh.loadASCII import import_mesh_file, \
595             export_mesh_file
596        import tempfile
597        import os
598
599        # create a .tsh file, no user outline
600        mesh_dic = {}
601        mesh_dic['vertices'] = [[0.76, 0.76],
602                                [0.76, 5.76],
603                                [5.76, 0.76]]
604        mesh_dic['triangles'] =  [[0, 2, 1]]
605        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
606        mesh_dic['triangle_tags'] = ['']
607        mesh_dic['vertex_attributes'] = [[], [], []]
608        mesh_dic['vertiex_attribute_titles'] = []
609        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
610        mesh_dic['segment_tags'] = ['external',
611                                    'external',
612                                    'external']
613        mesh_dic['geo_reference'] = Geo_reference(56,-0.76,-0.76)
614        mesh_file = tempfile.mktemp(".tsh")
615        export_mesh_file(mesh_file,mesh_dic)
616
617        geo_ref = Geo_reference(56,-200,-400)
618        # create an .xya file
619        point_file = tempfile.mktemp(".xya")
620        fd = open(point_file,'w')
621        fd.write("elevation, stage \n 201.0, 401.0,2.,4 \n 201.0, 403.0,4,8 \n 203.0, 401.0,4.,8 \n")
622        geo_ref.write_ASCII(fd)
623        fd.close()
624
625        mesh_output_file = tempfile.mktemp(".tsh")
626        fit_to_mesh_file(mesh_file,
627                         point_file,
628                         mesh_output_file,
629                         alpha = 0.0)
630        # load in the .tsh file we just wrote
631        mesh_dic = import_mesh_file(mesh_output_file)
632        #print "mesh_dic",mesh_dic
633        ans =[[0.0, 0.0],
634              [5.0, 10.0],
635              [5.0, 10.0]]
636        assert allclose(mesh_dic['vertex_attributes'],ans)
637
638        self.failUnless(mesh_dic['vertex_attribute_titles']  ==
639                        ['elevation','stage'],
640                        'test_fit_to_mesh_file failed')
641
642        #clean up
643        os.remove(mesh_file)
644        os.remove(point_file)
645        os.remove(mesh_output_file)
646
647    def test_fit_to_mesh_fileII(self):
648        from load_mesh.loadASCII import import_mesh_file, \
649             export_mesh_file
650        import tempfile
651        import os
652
653        # create a .tsh file, no user outline
654        mesh_dic = {}
655        mesh_dic['vertices'] = [[0.0, 0.0],
656                                [0.0, 5.0],
657                                [5.0, 0.0]]
658        mesh_dic['triangles'] =  [[0, 2, 1]]
659        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
660        mesh_dic['triangle_tags'] = ['']
661        mesh_dic['vertex_attributes'] = [[1,2], [1,2], [1,2]]
662        mesh_dic['vertex_attribute_titles'] = ['density', 'temp']
663        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
664        mesh_dic['segment_tags'] = ['external',
665                                                  'external',
666                                                  'external']
667        mesh_file = tempfile.mktemp(".tsh")
668        export_mesh_file(mesh_file,mesh_dic)
669
670        # create an .xya file
671        point_file = tempfile.mktemp(".xya")
672        fd = open(point_file,'w')
673        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")
674        fd.close()
675
676        mesh_output_file = "new_triangle.tsh"
677        fit_to_mesh_file(mesh_file,
678                         point_file,
679                         mesh_output_file,
680                         alpha = 0.0)
681        # load in the .tsh file we just wrote
682        mesh_dic = import_mesh_file(mesh_output_file)
683
684        assert allclose(mesh_dic['vertex_attributes'],
685                        [[1.0, 2.0,0.0, 0.0],
686                         [1.0, 2.0,5.0, 10.0],
687                         [1.0, 2.0,5.0,10.0]])
688
689        self.failUnless(mesh_dic['vertex_attribute_titles']  ==
690                        ['density', 'temp','elevation','stage'],
691                        'test_fit_to_mesh_file failed')
692
693        #clean up
694        os.remove(mesh_file)
695        os.remove(mesh_output_file)
696        os.remove(point_file)
697
698    def test_fit_to_mesh_file_errors(self):
699        from load_mesh.loadASCII import import_mesh_file, export_mesh_file
700        import tempfile
701        import os
702
703        # create a .tsh file, no user outline
704        mesh_dic = {}
705        mesh_dic['vertices'] = [[0.0, 0.0],[0.0, 5.0],[5.0, 0.0]]
706        mesh_dic['triangles'] =  [[0, 2, 1]]
707        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
708        mesh_dic['triangle_tags'] = ['']
709        mesh_dic['vertex_attributes'] = [[1,2], [1,2], [1,2]]
710        mesh_dic['vertex_attribute_titles'] = ['density', 'temp']
711        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
712        mesh_dic['segment_tags'] = ['external', 'external','external']
713        mesh_file = tempfile.mktemp(".tsh")
714        export_mesh_file(mesh_file,mesh_dic)
715
716        # create an .xya file
717        point_file = tempfile.mktemp(".xya")
718        fd = open(point_file,'w')
719        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")
720        fd.close()
721
722        mesh_output_file = "new_triangle.tsh"
723        try:
724            fit_to_mesh_file(mesh_file, point_file,
725                             mesh_output_file, display_errors = False)
726        except IOError:
727            pass
728        else:
729            #self.failUnless(0 ==1,  'Bad file did not raise error!')
730            raise 'Bad file did not raise error!'
731           
732        #clean up
733        os.remove(mesh_file)
734        os.remove(point_file)
735
736    def test_fit_to_mesh_file_errorsII(self):
737        from load_mesh.loadASCII import import_mesh_file, export_mesh_file
738        import tempfile
739        import os
740
741        # create a .tsh file, no user outline
742        mesh_file = tempfile.mktemp(".tsh")
743        fd = open(mesh_file,'w')
744        fd.write("unit testing a bad .tsh file \n")
745        fd.close()
746
747        # create an .xya file
748        point_file = tempfile.mktemp(".xya")
749        fd = open(point_file,'w')
750        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")
751        fd.close()
752
753        mesh_output_file = "new_triangle.tsh"
754        try:
755            fit_to_mesh_file(mesh_file, point_file,
756                             mesh_output_file, display_errors = False)
757        except IOError:
758            pass
759        else:
760            raise 'Bad file did not raise error!'
761           
762        #clean up
763        os.remove(mesh_file)
764        os.remove(point_file)
765
766    def test_fit_to_mesh_file_errorsIII(self):
767        from load_mesh.loadASCII import import_mesh_file, export_mesh_file
768        import tempfile
769        import os
770
771        # create a .tsh file, no user outline
772        mesh_dic = {}
773        mesh_dic['vertices'] = [[0.0, 0.0],[0.0, 5.0],[5.0, 0.0]]
774        mesh_dic['triangles'] =  [[0, 2, 1]]
775        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
776        mesh_dic['triangle_tags'] = ['']
777        mesh_dic['vertex_attributes'] = [[1,2], [1,2], [1,2]]
778        mesh_dic['vertex_attribute_titles'] = ['density', 'temp']
779        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
780        mesh_dic['segment_tags'] = ['external', 'external','external']
781        mesh_file = tempfile.mktemp(".tsh")
782        export_mesh_file(mesh_file,mesh_dic)
783
784        # create an .xya file
785        point_file = tempfile.mktemp(".xya")
786        fd = open(point_file,'w')
787        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")
788        fd.close()
789
790        #This a deliberately illegal filename to invoke the error.
791        mesh_output_file = ".../\z\z:ya.tsh"       
792
793        try:
794            fit_to_mesh_file(mesh_file, point_file,
795                             mesh_output_file, display_errors = False)
796        except IOError:
797            pass
798        else:
799            raise 'Bad file did not raise error!'
800       
801        #clean up
802        os.remove(mesh_file)
803        os.remove(point_file)
804
805
806    def Not_yet_test_smooth_att_to_mesh_with_excess_verts(self):
807
808        a = [0.0, 0.0]
809        b = [0.0, 5.0]
810        c = [5.0, 0.0]
811        d = [1.0, 1.0]
812        e = [18.0, 1000.0]
813       
814        points = [a, b, c, d, e]
815        triangles = [ [1,0,2] ]   #bac
816
817        d1 = [1.0, 1.0]
818        d2 = [1.0, 2.0]
819        d3 = [3.0,1.0]
820        d4 = [2.0,3.0]
821        d5 = [2.0,2.0]
822        d6 = [1.0,3.0]
823        data_coords = [d1, d2, d3, d4, d5, d6]
824        z = linear_function(data_coords)
825        #print "z",z
826
827        interp = Fit(points, triangles, alpha=0.0)
828       
829        try:
830            f = interp.fit(data_coords, z)
831        except VertsWithNoTrianglesError:
832            pass
833        else:
834            raise 'Verts with no triangles did not raise error!'
835       
836        #f = interp.fit(data_coords, z)
837        #answer = linear_function(points)
838
839        #  Removing the bad verts that we don't care about
840        #f = f[0:3]
841        #answer = answer[0:3]
842        #assert allclose(f, answer)
843
844#-------------------------------------------------------------
845if __name__ == "__main__":
846    suite = unittest.makeSuite(Test_Fit,'test')
847    #suite = unittest.makeSuite(Test_Fit,'test_smooth_att_to_mesh_with_excess_verts')
848    #suite = unittest.makeSuite(Test_Fit,'test_smooth_attributes_to_mesh_one_point')
849    runner = unittest.TextTestRunner(verbosity=1)
850    runner.run(suite)
851
852
853
854
855
Note: See TracBrowser for help on using the repository browser.