source: inundation/fit_interpolate/test_interpolate.py @ 3330

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

interpolate: adding warnings
File size: 47.8 KB
Line 
1#!/usr/bin/env python
2
3#TEST
4
5#import time, os
6
7
8import sys
9import os
10import unittest
11from math import sqrt
12import tempfile
13
14from Scientific.IO.NetCDF import NetCDFFile
15from Numeric import allclose, array, transpose, zeros, Float
16
17
18# ANUGA code imports
19from interpolate import *
20from coordinate_transforms.geo_reference import Geo_reference
21from shallow_water import Domain, Transmissive_boundary
22from utilities.numerical_tools import mean, INF
23from data_manager import get_dataobject
24from geospatial_data.geospatial_data import Geospatial_data
25
26def distance(x, y):
27    return sqrt( sum( (array(x)-array(y))**2 ))
28
29def linear_function(point):
30    point = array(point)
31    return point[:,0]+point[:,1]
32
33
34class Test_Interpolate(unittest.TestCase):
35
36    def setUp(self):
37
38        import time
39        from mesh_factory import rectangular
40
41
42        #Create basic mesh
43        points, vertices, boundary = rectangular(2, 2)
44
45        #Create shallow water domain
46        domain = Domain(points, vertices, boundary)
47        domain.default_order=2
48        domain.beta_h = 0
49
50
51        #Set some field values
52        domain.set_quantity('elevation', lambda x,y: -x)
53        domain.set_quantity('friction', 0.03)
54
55
56        ######################
57        # Boundary conditions
58        B = Transmissive_boundary(domain)
59        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
60
61
62        ######################
63        #Initial condition - with jumps
64
65        bed = domain.quantities['elevation'].vertex_values
66        stage = zeros(bed.shape, Float)
67
68        h = 0.3
69        for i in range(stage.shape[0]):
70            if i % 2 == 0:
71                stage[i,:] = bed[i,:] + h
72            else:
73                stage[i,:] = bed[i,:]
74
75        domain.set_quantity('stage', stage)
76
77        domain.distribute_to_vertices_and_edges()
78
79
80        self.domain = domain
81
82        C = domain.get_vertex_coordinates()
83        self.X = C[:,0:6:2].copy()
84        self.Y = C[:,1:6:2].copy()
85
86        self.F = bed
87
88
89
90    def tearDown(self):
91        pass
92
93    def test_datapoint_at_centroid(self):
94        a = [0.0, 0.0]
95        b = [0.0, 2.0]
96        c = [2.0,0.0]
97        points = [a, b, c]
98        vertices = [ [1,0,2] ]   #bac
99
100        data = [ [2.0/3, 2.0/3] ] #Use centroid as one data point
101
102        interp = Interpolate(points, vertices)
103        assert allclose(interp._build_interpolation_matrix_A(data).todense(),
104                        [[1./3, 1./3, 1./3]])
105
106    def test_quad_tree(self):
107        p0 = [-10.0, -10.0]
108        p1 = [20.0, -10.0]
109        p2 = [-10.0, 20.0]
110        p3 = [10.0, 50.0]
111        p4 = [30.0, 30.0]
112        p5 = [50.0, 10.0]
113        p6 = [40.0, 60.0]
114        p7 = [60.0, 40.0]
115        p8 = [-66.0, 20.0]
116        p9 = [10.0, -66.0]
117
118        points = [p0, p1, p2, p3, p4, p5, p6, p7, p8, p9]
119        triangles = [ [0, 1, 2],
120                      [3, 2, 4],
121                      [4, 2, 1],
122                      [4, 1, 5],
123                      [3, 4, 6],
124                      [6, 4, 7],
125                      [7, 4, 5],
126                      [8, 0, 2],
127                      [0, 9, 1]]
128
129        data = [ [4,4] ]
130        interp = Interpolate(points, triangles,
131                               max_vertices_per_cell = 4)
132        #print "PDSG - interp.get_A()", interp.get_A()
133        answer =  [ [ 0.06666667,  0.46666667,  0.46666667,  0.,
134                      0., 0. , 0., 0., 0., 0.]]
135
136       
137        assert allclose(interp._build_interpolation_matrix_A(data).todense(),
138                        answer)
139        #interp.set_point_coordinates([[-30, -30]]) #point outside of mesh
140        #print "PDSG - interp.get_A()", interp.get_A()
141        data = [[-30, -30]]
142        answer =  [ [ 0.0,  0.0,  0.0,  0.,
143                      0., 0. , 0., 0., 0., 0.]]
144       
145        assert allclose(interp._build_interpolation_matrix_A(data).todense(),
146                        answer)
147
148
149        #point outside of quad tree root cell
150        #interp.set_point_coordinates([[-70, -70]])
151        #print "PDSG - interp.get_A()", interp.get_A()
152        data = [[-70, -70]]
153        answer =  [ [ 0.0,  0.0,  0.0,  0.,
154                      0., 0. , 0., 0., 0., 0.]]
155        assert allclose(interp._build_interpolation_matrix_A(data).todense(),
156                        answer)
157
158    def test_datapoints_at_vertices(self):
159        #Test that data points coinciding with vertices yield a diagonal matrix
160       
161
162        a = [0.0, 0.0]
163        b = [0.0, 2.0]
164        c = [2.0,0.0]
165        points = [a, b, c]
166        vertices = [ [1,0,2] ]   #bac
167
168        data = points #Use data at vertices
169
170        interp = Interpolate(points, vertices)
171        answer = [[1., 0., 0.],
172                   [0., 1., 0.],
173                   [0., 0., 1.]]
174        assert allclose(interp._build_interpolation_matrix_A(data).todense(),
175                        answer)
176
177
178    def test_datapoints_on_edge_midpoints(self):
179        #Try datapoints midway on edges -
180        #each point should affect two matrix entries equally
181       
182
183        a = [0.0, 0.0]
184        b = [0.0, 2.0]
185        c = [2.0,0.0]
186        points = [a, b, c]
187        vertices = [ [1,0,2] ]   #bac
188
189        data = [ [0., 1.], [1., 0.], [1., 1.] ]
190        answer =  [[0.5, 0.5, 0.0],  #Affects vertex 1 and 0
191                    [0.5, 0.0, 0.5],  #Affects vertex 0 and 2
192                    [0.0, 0.5, 0.5]]
193        interp = Interpolate(points, vertices)
194
195        assert allclose(interp._build_interpolation_matrix_A(data).todense(),
196                        answer)
197
198    def test_datapoints_on_edges(self):
199        #Try datapoints on edges -
200        #each point should affect two matrix entries in proportion
201       
202
203        a = [0.0, 0.0]
204        b = [0.0, 2.0]
205        c = [2.0,0.0]
206        points = [a, b, c]
207        vertices = [ [1,0,2] ]   #bac
208
209        data = [ [0., 1.5], [1.5, 0.], [1.5, 0.5] ]
210        answer =  [[0.25, 0.75, 0.0],  #Affects vertex 1 and 0
211                   [0.25, 0.0, 0.75],  #Affects vertex 0 and 2
212                   [0.0, 0.25, 0.75]]
213
214        interp = Interpolate(points, vertices)
215
216        assert allclose(interp._build_interpolation_matrix_A(data).todense(),
217                        answer)
218
219
220    def test_arbitrary_datapoints(self):
221        #Try arbitrary datapoints
222       
223
224        from Numeric import sum
225
226        a = [0.0, 0.0]
227        b = [0.0, 2.0]
228        c = [2.0,0.0]
229        points = [a, b, c]
230        vertices = [ [1,0,2] ]   #bac
231
232        data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44] ]
233
234        interp = Interpolate(points, vertices)
235        #print "interp.get_A()", interp.get_A()
236        results = interp._build_interpolation_matrix_A(data).todense()
237        assert allclose(sum(results, axis=1), 1.0)
238
239    def test_arbitrary_datapoints_some_outside(self):
240        #Try arbitrary datapoints one outside the triangle.
241        #That one should be ignored
242       
243
244        from Numeric import sum
245
246        a = [0.0, 0.0]
247        b = [0.0, 2.0]
248        c = [2.0,0.0]
249        points = [a, b, c]
250        vertices = [ [1,0,2] ]   #bac
251
252        data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44], [5.0, 7.0]]
253
254        interp = Interpolate(points, vertices)
255        results = interp._build_interpolation_matrix_A(data).todense()
256        assert allclose(sum(results, axis=1), [1,1,1,0])
257
258
259
260    # this causes a memory error in scipy.sparse
261    def test_more_triangles(self):
262
263        a = [-1.0, 0.0]
264        b = [3.0, 4.0]
265        c = [4.0,1.0]
266        d = [-3.0, 2.0] #3
267        e = [-1.0,-2.0]
268        f = [1.0, -2.0] #5
269
270        points = [a, b, c, d,e,f]
271        triangles = [[0,1,3],[1,0,2],[0,4,5], [0,5,2]] #abd bac aef afc
272
273        #Data points
274        data = [ [-3., 2.0], [-2, 1], [0.0, 1], [0, 3], [2, 3], [-1.0/3,-4./3] ]
275        interp = Interpolate(points, triangles)
276
277        answer = [[0.0, 0.0, 0.0, 1.0, 0.0, 0.0],    #Affects point d
278                  [0.5, 0.0, 0.0, 0.5, 0.0, 0.0],    #Affects points a and d
279                  [0.75, 0.25, 0.0, 0.0, 0.0, 0.0],  #Affects points a and b
280                  [0.0, 0.5, 0.0, 0.5, 0.0, 0.0],    #Affects points a and d
281                  [0.25, 0.75, 0.0, 0.0, 0.0, 0.0],  #Affects points a and b
282                  [1./3, 0.0, 0.0, 0.0, 1./3, 1./3]] #Affects points a, e and f
283
284
285        A = interp._build_interpolation_matrix_A(data).todense()
286        for i in range(A.shape[0]):
287            for j in range(A.shape[1]):
288                if not allclose(A[i,j], answer[i][j]):
289                    print i,j,':',A[i,j], answer[i][j]
290
291
292        #results = interp._build_interpolation_matrix_A(data).todense()
293
294        assert allclose(A, answer)
295   
296    def test_geo_ref(self):
297        v0 = [0.0, 0.0]
298        v1 = [0.0, 5.0]
299        v2 = [5.0, 0.0]
300
301        vertices_absolute = [v0, v1, v2]
302        triangles = [ [1,0,2] ]   #bac
303
304        geo = Geo_reference(57,100, 500)
305
306        vertices = geo.change_points_geo_ref(vertices_absolute)
307        #print "vertices",vertices
308       
309        d0 = [1.0, 1.0]
310        d1 = [1.0, 2.0]
311        d2 = [3.0, 1.0]
312        point_coords = [ d0, d1, d2]
313
314        interp = Interpolate(vertices, triangles, mesh_origin=geo)
315        f = linear_function(vertices_absolute)
316        z = interp.interpolate(f, point_coords)
317        answer = linear_function(point_coords)
318
319        #print "z",z
320        #print "answer",answer
321        assert allclose(z, answer)
322
323       
324        z = interp.interpolate(f, point_coords, start_blocking_len = 2)
325        answer = linear_function(point_coords)
326
327        #print "z",z
328        #print "answer",answer
329        assert allclose(z, answer)
330       
331       
332    def test_Geospatial_verts(self):
333        v0 = [0.0, 0.0]
334        v1 = [0.0, 5.0]
335        v2 = [5.0, 0.0]
336
337        vertices_absolute = [v0, v1, v2]
338        triangles = [ [1,0,2] ]   #bac
339
340        geo = Geo_reference(57,100, 500)
341        vertices = geo.change_points_geo_ref(vertices_absolute)
342        geopoints = Geospatial_data(vertices,geo_reference = geo)
343        #print "vertices",vertices
344       
345        d0 = [1.0, 1.0]
346        d1 = [1.0, 2.0]
347        d2 = [3.0, 1.0]
348        point_coords = [ d0, d1, d2]
349
350        interp = Interpolate(geopoints, triangles)
351        f = linear_function(vertices_absolute)
352        z = interp.interpolate(f, point_coords)
353        answer = linear_function(point_coords)
354
355        #print "z",z
356        #print "answer",answer
357        assert allclose(z, answer)
358       
359        z = interp.interpolate(f, point_coords, start_blocking_len = 2)
360        answer = linear_function(point_coords)
361
362        #print "z",z
363        #print "answer",answer
364        assert allclose(z, answer)
365       
366    def test_interpolate_attributes_to_points(self):
367        v0 = [0.0, 0.0]
368        v1 = [0.0, 5.0]
369        v2 = [5.0, 0.0]
370
371        vertices = [v0, v1, v2]
372        triangles = [ [1,0,2] ]   #bac
373
374        d0 = [1.0, 1.0]
375        d1 = [1.0, 2.0]
376        d2 = [3.0, 1.0]
377        point_coords = [ d0, d1, d2]
378
379        interp = Interpolate(vertices, triangles)
380        f = linear_function(vertices)
381        z = interp.interpolate(f, point_coords)
382        answer = linear_function(point_coords)
383
384        #print "z",z
385        #print "answer",answer
386        assert allclose(z, answer)
387
388
389        z = interp.interpolate(f, point_coords, start_blocking_len = 2)
390        answer = linear_function(point_coords)
391
392        #print "z",z
393        #print "answer",answer
394        assert allclose(z, answer)
395
396    def test_interpolate_attributes_to_pointsII(self):
397        a = [-1.0, 0.0]
398        b = [3.0, 4.0]
399        c = [4.0, 1.0]
400        d = [-3.0, 2.0] #3
401        e = [-1.0, -2.0]
402        f = [1.0, -2.0] #5
403
404        vertices = [a, b, c, d,e,f]
405        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
406
407
408        point_coords = [[-2.0, 2.0],
409                        [-1.0, 1.0],
410                        [0.0, 2.0],
411                        [1.0, 1.0],
412                        [2.0, 1.0],
413                        [0.0, 0.0],
414                        [1.0, 0.0],
415                        [0.0, -1.0],
416                        [-0.2, -0.5],
417                        [-0.9, -1.5],
418                        [0.5, -1.9],
419                        [3.0, 1.0]]
420
421        interp = Interpolate(vertices, triangles)
422        f = linear_function(vertices)
423        z = interp.interpolate(f, point_coords)
424        answer = linear_function(point_coords)
425        #print "z",z
426        #print "answer",answer
427        assert allclose(z, answer)
428
429        z = interp.interpolate(f, point_coords, start_blocking_len = 2)
430        answer = linear_function(point_coords)
431
432        #print "z",z
433        #print "answer",answer
434        assert allclose(z, answer)
435       
436    def test_interpolate_attributes_to_pointsIII(self):
437        #Test linear interpolation of known values at vertices to
438        #new points inside a triangle
439       
440        a = [0.0, 0.0]
441        b = [0.0, 5.0]
442        c = [5.0, 0.0]
443        d = [5.0, 5.0]
444
445        vertices = [a, b, c, d]
446        triangles = [ [1,0,2], [2,3,1] ]   #bac, cdb
447
448        #Points within triangle 1
449        d0 = [1.0, 1.0]
450        d1 = [1.0, 2.0]
451        d2 = [3.0, 1.0]
452
453        #Point within triangle 2
454        d3 = [4.0, 3.0]
455
456        #Points on common edge
457        d4 = [2.5, 2.5]
458        d5 = [4.0, 1.0]
459
460        #Point on common vertex
461        d6 = [0., 5.]
462       
463        point_coords = [d0, d1, d2, d3, d4, d5, d6]
464
465        interp = Interpolate(vertices, triangles)
466
467        #Known values at vertices
468        #Functions are x+y, x+2y, 2x+y, x-y-5
469        f = [ [0., 0., 0., -5.],        # (0,0)
470              [5., 10., 5., -10.],      # (0,5)
471              [5., 5., 10.0, 0.],       # (5,0)
472              [10., 15., 15., -5.]]     # (5,5)
473
474        z = interp.interpolate(f, point_coords)
475        answer = [ [2., 3., 3., -5.],   # (1,1)
476                   [3., 5., 4., -6.],   # (1,2)
477                   [4., 5., 7., -3.],   # (3,1)
478                   [7., 10., 11., -4.], # (4,3)
479                   [5., 7.5, 7.5, -5.], # (2.5, 2.5)
480                   [5., 6., 9., -2.],   # (4,1)
481                   [5., 10., 5., -10.]]  # (0,5)
482
483        #print "***********"
484        #print "z",z
485        #print "answer",answer
486        #print "***********"
487
488        assert allclose(z, answer)
489
490
491        z = interp.interpolate(f, point_coords, start_blocking_len = 2)
492
493        #print "z",z
494        #print "answer",answer
495        assert allclose(z, answer)
496       
497    def test_interpolate_point_outside_of_mesh(self):
498        #Test linear interpolation of known values at vertices to
499        #new points inside a triangle
500       
501        a = [0.0, 0.0]
502        b = [0.0, 5.0]
503        c = [5.0, 0.0]
504        d = [5.0, 5.0]
505
506        vertices = [a, b, c, d]
507        triangles = [ [1,0,2], [2,3,1] ]   #bac, cdb
508
509        #Far away point
510        d7 = [-1., -1.]
511       
512        point_coords = [ d7]
513        interp = Interpolate(vertices, triangles)
514
515        #Known values at vertices
516        #Functions are x+y, x+2y, 2x+y, x-y-5
517        f = [ [0., 0., 0., -5.],        # (0,0)
518              [5., 10., 5., -10.],      # (0,5)
519              [5., 5., 10.0, 0.],       # (5,0)
520              [10., 15., 15., -5.]]     # (5,5)
521
522        z = interp.interpolate(f, point_coords) #, verbose=True)
523        answer = array([ [INF, INF, INF, INF]]) # (-1,-1)
524
525        #print "***********"
526        #print "z",z
527        #print "answer",answer
528        #print "***********"
529
530        #Should an error message be returned if points are outside
531        # of the mesh?
532        # A warning message is printed, if verbose is on.
533
534        for i in range(4):
535            self.failUnless( z[0,i] == answer[0,i], 'Fail!')
536       
537        z = interp.interpolate(f, point_coords, start_blocking_len = 2)
538
539        #print "z",z
540        #print "answer",answer
541       
542        for i in range(4):
543            self.failUnless( z[0,i] == answer[0,i], 'Fail!')
544       
545       
546    def test_interpolate_attributes_to_pointsIV(self):
547        a = [-1.0, 0.0]
548        b = [3.0, 4.0]
549        c = [4.0, 1.0]
550        d = [-3.0, 2.0] #3
551        e = [-1.0, -2.0]
552        f = [1.0, -2.0] #5
553
554        vertices = [a, b, c, d,e,f]
555        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
556
557
558        point_coords = [[-2.0, 2.0],
559                        [-1.0, 1.0],
560                        [0.0, 2.0],
561                        [1.0, 1.0],
562                        [2.0, 1.0],
563                        [0.0, 0.0],
564                        [1.0, 0.0],
565                        [0.0, -1.0],
566                        [-0.2, -0.5],
567                        [-0.9, -1.5],
568                        [0.5, -1.9],
569                        [3.0, 1.0]]
570
571        interp = Interpolate(vertices, triangles)
572        f = array([linear_function(vertices),2*linear_function(vertices) ])
573        f = transpose(f)
574        #print "f",f
575        z = interp.interpolate(f, point_coords)
576        answer = [linear_function(point_coords),
577                  2*linear_function(point_coords) ]
578        answer = transpose(answer)
579        #print "z",z
580        #print "answer",answer
581        assert allclose(z, answer)
582
583        z = interp.interpolate(f, point_coords, start_blocking_len = 2)
584
585        #print "z",z
586        #print "answer",answer
587        assert allclose(z, answer)
588
589    def test_interpolate_blocking(self):
590        a = [-1.0, 0.0]
591        b = [3.0, 4.0]
592        c = [4.0, 1.0]
593        d = [-3.0, 2.0] #3
594        e = [-1.0, -2.0]
595        f = [1.0, -2.0] #5
596
597        vertices = [a, b, c, d,e,f]
598        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
599
600
601        point_coords = [[-2.0, 2.0],
602                        [-1.0, 1.0],
603                        [0.0, 2.0],
604                        [1.0, 1.0],
605                        [2.0, 1.0],
606                        [0.0, 0.0],
607                        [1.0, 0.0],
608                        [0.0, -1.0],
609                        [-0.2, -0.5],
610                        [-0.9, -1.5],
611                        [0.5, -1.9],
612                        [3.0, 1.0]]
613
614        interp = Interpolate(vertices, triangles)
615        f = array([linear_function(vertices),2*linear_function(vertices) ])
616        f = transpose(f)
617        #print "f",f
618        for blocking_max in range(len(point_coords)+2):
619        #if True:
620         #   blocking_max = 5
621            z = interp.interpolate(f, point_coords,
622                                   start_blocking_len=blocking_max)
623            answer = [linear_function(point_coords),
624                      2*linear_function(point_coords) ]
625            answer = transpose(answer)
626            #print "z",z
627            #print "answer",answer
628            assert allclose(z, answer)
629           
630        f = array([linear_function(vertices),2*linear_function(vertices),
631                   2*linear_function(vertices) - 100  ])
632        f = transpose(f)
633        #print "f",f
634        for blocking_max in range(len(point_coords)+2):
635        #if True:
636         #   blocking_max = 5
637            z = interp.interpolate(f, point_coords,
638                                   start_blocking_len=blocking_max)
639            answer = array([linear_function(point_coords),
640                      2*linear_function(point_coords) ,
641                      2*linear_function(point_coords)-100 ])
642            z = transpose(z)
643            #print "z",z
644            #print "answer",answer
645            assert allclose(z, answer)
646
647    def test_interpolate_geo_spatial(self):
648        a = [-1.0, 0.0]
649        b = [3.0, 4.0]
650        c = [4.0, 1.0]
651        d = [-3.0, 2.0] #3
652        e = [-1.0, -2.0]
653        f = [1.0, -2.0] #5
654
655        vertices = [a, b, c, d,e,f]
656        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
657
658
659        point_coords_absolute = [[-2.0, 2.0],
660                        [-1.0, 1.0],
661                        [0.0, 2.0],
662                        [1.0, 1.0],
663                        [2.0, 1.0],
664                        [0.0, 0.0],
665                        [1.0, 0.0],
666                        [0.0, -1.0],
667                        [-0.2, -0.5],
668                        [-0.9, -1.5],
669                        [0.5, -1.9],
670                        [3.0, 1.0]]
671
672        geo = Geo_reference(57,100, 500)
673        point_coords = geo.change_points_geo_ref(point_coords_absolute)
674        point_coords = Geospatial_data(point_coords,geo_reference = geo)
675       
676        interp = Interpolate(vertices, triangles)
677        f = array([linear_function(vertices),2*linear_function(vertices) ])
678        f = transpose(f)
679        #print "f",f
680        for blocking_max in range(14):
681        #if True:
682         #   blocking_max = 5
683            z = interp.interpolate(f, point_coords,
684                                   start_blocking_len=blocking_max)
685            answer = [linear_function(point_coords.get_data_points( \
686                      absolute = True)),
687                      2*linear_function(point_coords.get_data_points( \
688                      absolute = True)) ]
689            answer = transpose(answer)
690            #print "z",z
691            #print "answer",answer
692            assert allclose(z, answer)
693           
694        f = array([linear_function(vertices),2*linear_function(vertices),
695                   2*linear_function(vertices) - 100  ])
696        f = transpose(f)
697        #print "f",f
698        for blocking_max in range(14):
699        #if True:
700         #   blocking_max = 5
701            z = interp.interpolate(f, point_coords,
702                                   start_blocking_len=blocking_max)
703            answer = array([linear_function(point_coords.get_data_points( \
704                      absolute = True)),
705                      2*linear_function(point_coords.get_data_points( \
706                      absolute = True)) ,
707                      2*linear_function(point_coords.get_data_points( \
708                      absolute = True))-100 ])
709            z = transpose(z)
710            #print "z",z
711            #print "answer",answer
712            assert allclose(z, answer)
713
714        z = interp.interpolate(f, point_coords, start_blocking_len = 2)
715
716        #print "z",z
717        #print "answer",answer
718        assert allclose(z, answer)
719       
720    def test_interpolate_geo_spatial(self):
721        a = [-1.0, 0.0]
722        b = [3.0, 4.0]
723        c = [4.0, 1.0]
724        d = [-3.0, 2.0] #3
725        e = [-1.0, -2.0]
726        f = [1.0, -2.0] #5
727
728        vertices = [a, b, c, d,e,f]
729        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
730
731
732        point_coords_absolute = [[-2.0, 2.0],
733                        [-1.0, 1.0],
734                        [0.0, 2.0],
735                        [1.0, 1.0],
736                        [2.0, 1.0],
737                        [0.0, 0.0],
738                        [1.0, 0.0],
739                        [0.0, -1.0],
740                        [-0.2, -0.5],
741                        [-0.9, -1.5],
742                        [0.5, -1.9],
743                        [3.0, 1.0]]
744
745        geo = Geo_reference(57,100, 500)
746        point_coords = geo.change_points_geo_ref(point_coords_absolute)
747        point_coords = Geospatial_data(point_coords,geo_reference = geo)
748       
749        interp = Interpolate(vertices, triangles)
750        f = array([linear_function(vertices),2*linear_function(vertices) ])
751        f = transpose(f)
752        #print "f",f
753        z = interp.interpolate_block(f, point_coords)
754        answer = [linear_function(point_coords.get_data_points( \
755                      absolute = True)),
756                  2*linear_function(point_coords.get_data_points( \
757                      absolute = True)) ]
758        answer = transpose(answer)
759        #print "z",z
760        #print "answer",answer
761        assert allclose(z, answer)
762           
763        z = interp.interpolate(f, point_coords, start_blocking_len = 2)
764
765        #print "z",z
766        #print "answer",answer
767        assert allclose(z, answer)
768
769       
770    def test_interpolate_reuse(self):
771        a = [-1.0, 0.0]
772        b = [3.0, 4.0]
773        c = [4.0, 1.0]
774        d = [-3.0, 2.0] #3
775        e = [-1.0, -2.0]
776        f = [1.0, -2.0] #5
777
778        vertices = [a, b, c, d,e,f]
779        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
780
781
782        point_coords = [[-2.0, 2.0],
783                        [-1.0, 1.0],
784                        [0.0, 2.0],
785                        [1.0, 1.0],
786                        [2.0, 1.0],
787                        [0.0, 0.0],
788                        [1.0, 0.0],
789                        [0.0, -1.0],
790                        [-0.2, -0.5],
791                        [-0.9, -1.5],
792                        [0.5, -1.9],
793                        [3.0, 1.0]]
794
795        interp = Interpolate(vertices, triangles)
796        f = array([linear_function(vertices),2*linear_function(vertices) ])
797        f = transpose(f)
798        z = interp.interpolate(f, point_coords,
799                               start_blocking_len=20)
800        answer = [linear_function(point_coords),
801                  2*linear_function(point_coords) ]
802        answer = transpose(answer)
803        #print "z",z
804        #print "answer",answer
805        assert allclose(z, answer)
806        assert allclose(interp._A_can_be_reused, True)
807
808        z = interp.interpolate(f)
809        assert allclose(z, answer)
810       
811        # This causes blocking to occur.
812        z = interp.interpolate(f, start_blocking_len=10)
813        assert allclose(z, answer)
814        assert allclose(interp._A_can_be_reused, False)
815
816        #A is recalculated
817        z = interp.interpolate(f)
818        assert allclose(z, answer)
819        assert allclose(interp._A_can_be_reused, True)
820       
821        interp = Interpolate(vertices, triangles)
822        #Must raise an exception, no points specified
823        try:
824            z = interp.interpolate(f)
825        except:
826            pass
827       
828
829
830    def test_interpolation_interface_time_only(self):
831
832        # Test spatio-temporal interpolation
833        # Test that spatio temporal function performs the correct
834        # interpolations in both time and space
835       
836
837
838        #Three timesteps
839        time = [1.0, 5.0, 6.0]
840       
841
842        #One quantity
843        Q = zeros( (3,6), Float )
844
845        #Linear in time and space
846        a = [0.0, 0.0]
847        b = [0.0, 2.0]
848        c = [2.0, 0.0]
849        d = [0.0, 4.0]
850        e = [2.0, 2.0]
851        f = [4.0, 0.0]
852
853        points = [a, b, c, d, e, f]
854       
855        for i, t in enumerate(time):
856            Q[i, :] = t*linear_function(points)
857
858           
859        #Check basic interpolation of one quantity using averaging
860        #(no interpolation points or spatial info)
861        I = Interpolation_function(time, [mean(Q[0,:]),
862                                          mean(Q[1,:]),
863                                          mean(Q[2,:])])
864
865
866
867        #Check temporal interpolation
868        for i in [0,1,2]:
869            assert allclose(I(time[i]), mean(Q[i,:]))
870
871        #Midway   
872        assert allclose(I( (time[0] + time[1])/2 ),
873                        (I(time[0]) + I(time[1]))/2 )
874
875        assert allclose(I( (time[1] + time[2])/2 ),
876                        (I(time[1]) + I(time[2]))/2 )
877
878        assert allclose(I( (time[0] + time[2])/2 ),
879                        (I(time[0]) + I(time[2]))/2 )                 
880
881        #1/3
882        assert allclose(I( (time[0] + time[2])/3 ),
883                        (I(time[0]) + I(time[2]))/3 )                         
884
885
886        #Out of bounds checks
887        try:
888            I(time[0]-1) 
889        except:
890            pass
891        else:
892            raise 'Should raise exception'
893
894        try:
895            I(time[-1]+1) 
896        except:
897            pass
898        else:
899            raise 'Should raise exception'       
900
901
902       
903
904    def test_interpolation_interface_spatial_only(self):
905        # Test spatio-temporal interpolation with constant time
906       
907        #Three timesteps
908        time = [1.0, 5.0, 6.0]
909       
910       
911        #Setup mesh used to represent fitted function
912        a = [0.0, 0.0]
913        b = [0.0, 2.0]
914        c = [2.0, 0.0]
915        d = [0.0, 4.0]
916        e = [2.0, 2.0]
917        f = [4.0, 0.0]
918
919        points = [a, b, c, d, e, f]
920        #bac, bce, ecf, dbe
921        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
922
923
924        #New datapoints where interpolated values are sought
925        interpolation_points = [[ 0.0, 0.0],
926                                [ 0.5, 0.5],
927                                [ 0.7, 0.7],
928                                [ 1.0, 0.5],
929                                [ 2.0, 0.4],
930                                [ 2.8, 1.2]]
931
932
933        #One quantity linear in space
934        Q = linear_function(points)
935
936
937        #Check interpolation of one quantity using interpolaton points
938        I = Interpolation_function(time, Q,
939                                   vertex_coordinates = points,
940                                   triangles = triangles, 
941                                   interpolation_points = interpolation_points,
942                                   verbose = False)
943
944
945        answer = linear_function(interpolation_points)
946
947        t = time[0]
948        for j in range(50): #t in [1, 6]
949            for id in range(len(interpolation_points)):
950                assert allclose(I(t, id), answer[id])
951
952            t += 0.1   
953
954
955        try:   
956            I(1)
957        except:
958            pass
959        else:
960            raise 'Should raise exception'
961           
962
963
964    def test_interpolation_interface(self):
965        # Test spatio-temporal interpolation
966        # Test that spatio temporal function performs the correct
967        # interpolations in both time and space
968   
969        #Three timesteps
970        time = [1.0, 5.0, 6.0]   
971
972        #Setup mesh used to represent fitted function
973        a = [0.0, 0.0]
974        b = [0.0, 2.0]
975        c = [2.0, 0.0]
976        d = [0.0, 4.0]
977        e = [2.0, 2.0]
978        f = [4.0, 0.0]
979
980        points = [a, b, c, d, e, f]
981        #bac, bce, ecf, dbe
982        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
983
984
985        #New datapoints where interpolated values are sought
986        interpolation_points = [[ 0.0, 0.0],
987                                [ 0.5, 0.5],
988                                [ 0.7, 0.7],
989                                [ 1.0, 0.5],
990                                [ 2.0, 0.4],
991                                [ 2.8, 1.2]]
992
993        #One quantity
994        Q = zeros( (3,6), Float )
995
996        #Linear in time and space
997        for i, t in enumerate(time):
998            Q[i, :] = t*linear_function(points)
999
1000        #Check interpolation of one quantity using interpolaton points)
1001        I = Interpolation_function(time, Q,
1002                                   vertex_coordinates = points,
1003                                   triangles = triangles, 
1004                                   interpolation_points = interpolation_points,
1005                                   verbose = False)
1006
1007        answer = linear_function(interpolation_points)
1008       
1009        t = time[0]
1010        for j in range(50): #t in [1, 6]
1011            for id in range(len(interpolation_points)):
1012                assert allclose(I(t, id), t*answer[id])
1013            t += 0.1   
1014
1015        try:   
1016            I(1)
1017        except:
1018            pass
1019        else:
1020            raise 'Should raise exception'
1021
1022
1023    def test_interpolation_precompute_points(self):
1024        # looking at a discrete mesh
1025        #
1026   
1027        #Three timesteps
1028        time = [0.0, 60.0]   
1029
1030        #Setup mesh used to represent fitted function
1031        points = [[ 15., -20.],
1032                  [ 15.,  10.],
1033                  [  0., -20.],
1034                  [  0.,  10.],
1035                  [  0., -20.],
1036                  [ 15.,  10.]]
1037       
1038        triangles = [[0, 1, 2],
1039                     [3, 4, 5]]
1040
1041        #New datapoints where interpolated values are sought
1042        interpolation_points = [[ 1.,  0.], [0.,1.]]
1043
1044        #One quantity
1045        Q = zeros( (2,6), Float )
1046
1047        #Linear in time and space
1048        for i, t in enumerate(time):
1049            Q[i, :] = t*linear_function(points)
1050        #print "Q", Q
1051
1052
1053       
1054        interp = Interpolate(points, triangles)
1055        f = array([linear_function(points),2*linear_function(points) ])
1056        f = transpose(f)
1057        #print "f",f
1058        z = interp.interpolate(f, interpolation_points)
1059        answer = [linear_function(interpolation_points),
1060                  2*linear_function(interpolation_points) ]
1061        answer = transpose(answer)
1062        #print "z",z
1063        #print "answer",answer
1064        assert allclose(z, answer)
1065
1066
1067        #Check interpolation of one quantity using interpolaton points)
1068        I = Interpolation_function(time, Q,
1069                                   vertex_coordinates = points,
1070                                   triangles = triangles, 
1071                                   interpolation_points = interpolation_points,
1072                                   verbose = False)
1073       
1074        #print "I.precomputed_values", I.precomputed_values
1075
1076        msg = 'Interpolation failed'
1077        assert allclose(I.precomputed_values['Attribute'][1], [60, 60]), msg
1078        #self.failUnless( I.precomputed_values['Attribute'][1] == 60.0,
1079        #                ' failed')
1080       
1081    def test_interpolation_function_outside_point(self):
1082        # Test spatio-temporal interpolation
1083        # Test that spatio temporal function performs the correct
1084        # interpolations in both time and space
1085   
1086        #Three timesteps
1087        time = [1.0, 5.0, 6.0]   
1088
1089        #Setup mesh used to represent fitted function
1090        a = [0.0, 0.0]
1091        b = [0.0, 2.0]
1092        c = [2.0, 0.0]
1093        d = [0.0, 4.0]
1094        e = [2.0, 2.0]
1095        f = [4.0, 0.0]
1096
1097        points = [a, b, c, d, e, f]
1098        #bac, bce, ecf, dbe
1099        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1100
1101
1102        #New datapoints where interpolated values are sought
1103        interpolation_points = [[ 0.0, 0.0],
1104                                [ 0.5, 0.5],
1105                                [ 0.7, 0.7],
1106                                [ 1.0, 0.5],
1107                                [ 2.0, 0.4],
1108                                [ 545354534, 4354354353]] # outside the mesh
1109
1110        #One quantity
1111        Q = zeros( (3,6), Float )
1112
1113        #Linear in time and space
1114        for i, t in enumerate(time):
1115            Q[i, :] = t*linear_function(points)
1116
1117        #Check interpolation of one quantity using interpolaton points)
1118        I = Interpolation_function(time, Q,
1119                                   vertex_coordinates = points,
1120                                   triangles = triangles, 
1121                                   interpolation_points = interpolation_points,
1122                                   verbose = False)
1123
1124        answer = linear_function(interpolation_points)
1125         
1126        t = time[0]
1127        for j in range(50): #t in [1, 6]
1128            for id in range(len(interpolation_points)-1):
1129                assert allclose(I(t, id), t*answer[id])
1130            t += 0.1
1131           
1132        # Now test the point outside the mesh
1133        t = time[0]
1134        for j in range(50): #t in [1, 6]
1135            self.failUnless(I(t, 5) == INF, 'Fail!')
1136            t += 0.1 
1137           
1138        try:   
1139            I(1)
1140        except:
1141            pass
1142        else:
1143            raise 'Should raise exception'
1144
1145
1146    def test_interpolation_function_time(self):
1147        #Test a long time series with an error in it (this did cause an
1148        #error once)
1149       
1150
1151        time = array(\
1152        [0.00000000e+00, 5.00000000e-02, 1.00000000e-01,   1.50000000e-01,
1153        2.00000000e-01,   2.50000000e-01,   3.00000000e-01,   3.50000000e-01,
1154        4.00000000e-01,   4.50000000e-01,   5.00000000e-01,   5.50000000e-01,
1155        6.00000000e-01,   6.50000000e-01,   7.00000000e-01,   7.50000000e-01,
1156        8.00000000e-01,   8.50000000e-01,   9.00000000e-01,   9.50000000e-01,
1157        1.00000000e-00,   1.05000000e+00,   1.10000000e+00,   1.15000000e+00,
1158        1.20000000e+00,   1.25000000e+00,   1.30000000e+00,   1.35000000e+00,
1159        1.40000000e+00,   1.45000000e+00,   1.50000000e+00,   1.55000000e+00,
1160        1.60000000e+00,   1.65000000e+00,   1.70000000e+00,   1.75000000e+00,
1161        1.80000000e+00,   1.85000000e+00,   1.90000000e+00,   1.95000000e+00,
1162        2.00000000e+00,   2.05000000e+00,   2.10000000e+00,   2.15000000e+00,
1163        2.20000000e+00,   2.25000000e+00,   2.30000000e+00,   2.35000000e+00,
1164        2.40000000e+00,   2.45000000e+00,   2.50000000e+00,   2.55000000e+00,
1165        2.60000000e+00,   2.65000000e+00,   2.70000000e+00,   2.75000000e+00,
1166        2.80000000e+00,   2.85000000e+00,   2.90000000e+00,   2.95000000e+00,
1167        3.00000000e+00,   3.05000000e+00,   9.96920997e+36,   3.15000000e+00,
1168        3.20000000e+00,   3.25000000e+00,   3.30000000e+00,   3.35000000e+00,
1169        3.40000000e+00,   3.45000000e+00,   3.50000000e+00,   3.55000000e+00,
1170        3.60000000e+00,   3.65000000e+00,   3.70000000e+00,   3.75000000e+00,
1171        3.80000000e+00,   3.85000000e+00,   3.90000000e+00,   3.95000000e+00,
1172        4.00000000e+00,   4.05000000e+00,   4.10000000e+00,   4.15000000e+00,
1173        4.20000000e+00,   4.25000000e+00,   4.30000000e+00,   4.35000000e+00,
1174        4.40000000e+00,   4.45000000e+00,   4.50000000e+00,   4.55000000e+00,
1175        4.60000000e+00,   4.65000000e+00,   4.70000000e+00,   4.75000000e+00,
1176        4.80000000e+00,   4.85000000e+00,   4.90000000e+00,   4.95000000e+00,
1177        5.00000000e+00,   5.05000000e+00,   5.10000000e+00,   5.15000000e+00,
1178        5.20000000e+00,   5.25000000e+00,   5.30000000e+00,   5.35000000e+00,
1179        5.40000000e+00,   5.45000000e+00,   5.50000000e+00,   5.55000000e+00,
1180        5.60000000e+00,   5.65000000e+00,   5.70000000e+00,   5.75000000e+00,
1181        5.80000000e+00,   5.85000000e+00,   5.90000000e+00,   5.95000000e+00,
1182        6.00000000e+00,   6.05000000e+00,   6.10000000e+00,   6.15000000e+00,
1183        6.20000000e+00,   6.25000000e+00,   6.30000000e+00,   6.35000000e+00,
1184        6.40000000e+00,   6.45000000e+00,   6.50000000e+00,   6.55000000e+00,
1185        6.60000000e+00,   6.65000000e+00,   6.70000000e+00,   6.75000000e+00,
1186        6.80000000e+00,   6.85000000e+00,   6.90000000e+00,   6.95000000e+00,
1187        7.00000000e+00,   7.05000000e+00,   7.10000000e+00,   7.15000000e+00,
1188        7.20000000e+00,   7.25000000e+00,   7.30000000e+00,   7.35000000e+00,
1189        7.40000000e+00,   7.45000000e+00,   7.50000000e+00,   7.55000000e+00,
1190        7.60000000e+00,   7.65000000e+00,   7.70000000e+00,   7.75000000e+00,
1191        7.80000000e+00,   7.85000000e+00,   7.90000000e+00,   7.95000000e+00,
1192        8.00000000e+00,   8.05000000e+00,   8.10000000e+00,   8.15000000e+00,
1193        8.20000000e+00,   8.25000000e+00,   8.30000000e+00,   8.35000000e+00,
1194        8.40000000e+00,   8.45000000e+00,   8.50000000e+00,   8.55000000e+00,
1195        8.60000000e+00,   8.65000000e+00,   8.70000000e+00,   8.75000000e+00,
1196        8.80000000e+00,   8.85000000e+00,   8.90000000e+00,   8.95000000e+00,
1197        9.00000000e+00,   9.05000000e+00,   9.10000000e+00,   9.15000000e+00,
1198        9.20000000e+00,   9.25000000e+00,   9.30000000e+00,   9.35000000e+00,
1199        9.40000000e+00,   9.45000000e+00,   9.50000000e+00,   9.55000000e+00,
1200        9.60000000e+00,   9.65000000e+00,   9.70000000e+00,   9.75000000e+00,
1201        9.80000000e+00,   9.85000000e+00,   9.90000000e+00,   9.95000000e+00,
1202        1.00000000e+01,   1.00500000e+01,   1.01000000e+01,   1.01500000e+01,
1203        1.02000000e+01,   1.02500000e+01,   1.03000000e+01,   1.03500000e+01,
1204        1.04000000e+01,   1.04500000e+01,   1.05000000e+01,   1.05500000e+01,
1205        1.06000000e+01,   1.06500000e+01,   1.07000000e+01,   1.07500000e+01,
1206        1.08000000e+01,   1.08500000e+01,   1.09000000e+01,   1.09500000e+01,
1207        1.10000000e+01,   1.10500000e+01,   1.11000000e+01,   1.11500000e+01,
1208        1.12000000e+01,   1.12500000e+01,   1.13000000e+01,   1.13500000e+01,
1209        1.14000000e+01,   1.14500000e+01,   1.15000000e+01,   1.15500000e+01,
1210        1.16000000e+01,   1.16500000e+01,   1.17000000e+01,   1.17500000e+01,
1211        1.18000000e+01,   1.18500000e+01,   1.19000000e+01,   1.19500000e+01,
1212        1.20000000e+01,   1.20500000e+01,   1.21000000e+01,   1.21500000e+01,
1213        1.22000000e+01,   1.22500000e+01,   1.23000000e+01,   1.23500000e+01,
1214        1.24000000e+01,   1.24500000e+01,   1.25000000e+01,   1.25500000e+01,
1215        1.26000000e+01,   1.26500000e+01,   1.27000000e+01,   1.27500000e+01,
1216        1.28000000e+01,   1.28500000e+01,   1.29000000e+01,   1.29500000e+01,
1217        1.30000000e+01,   1.30500000e+01,   1.31000000e+01,   1.31500000e+01,
1218        1.32000000e+01,   1.32500000e+01,   1.33000000e+01,   1.33500000e+01,
1219        1.34000000e+01,   1.34500000e+01,   1.35000000e+01,   1.35500000e+01,
1220        1.36000000e+01,   1.36500000e+01,   1.37000000e+01,   1.37500000e+01,
1221        1.38000000e+01,   1.38500000e+01,   1.39000000e+01,   1.39500000e+01,
1222        1.40000000e+01,   1.40500000e+01,   1.41000000e+01,   1.41500000e+01,
1223        1.42000000e+01,   1.42500000e+01,   1.43000000e+01,   1.43500000e+01,
1224        1.44000000e+01,   1.44500000e+01,   1.45000000e+01,   1.45500000e+01,
1225        1.46000000e+01,   1.46500000e+01,   1.47000000e+01,   1.47500000e+01,
1226        1.48000000e+01,   1.48500000e+01,   1.49000000e+01,   1.49500000e+01,
1227        1.50000000e+01,   1.50500000e+01,   1.51000000e+01,   1.51500000e+01,
1228        1.52000000e+01,   1.52500000e+01,   1.53000000e+01,   1.53500000e+01,
1229        1.54000000e+01,   1.54500000e+01,   1.55000000e+01,   1.55500000e+01,
1230        1.56000000e+01,   1.56500000e+01,   1.57000000e+01,   1.57500000e+01,
1231        1.58000000e+01,   1.58500000e+01,   1.59000000e+01,   1.59500000e+01,
1232        1.60000000e+01,   1.60500000e+01,   1.61000000e+01,   1.61500000e+01,
1233        1.62000000e+01,   1.62500000e+01,   1.63000000e+01,   1.63500000e+01,
1234        1.64000000e+01,   1.64500000e+01,   1.65000000e+01,   1.65500000e+01,
1235        1.66000000e+01,   1.66500000e+01,   1.67000000e+01,   1.67500000e+01,
1236        1.68000000e+01,   1.68500000e+01,   1.69000000e+01,   1.69500000e+01,
1237        1.70000000e+01,   1.70500000e+01,   1.71000000e+01,   1.71500000e+01,
1238        1.72000000e+01,   1.72500000e+01,   1.73000000e+01,   1.73500000e+01,
1239        1.74000000e+01,   1.74500000e+01,   1.75000000e+01,   1.75500000e+01,
1240        1.76000000e+01,   1.76500000e+01,   1.77000000e+01,   1.77500000e+01,
1241        1.78000000e+01,   1.78500000e+01,   1.79000000e+01,   1.79500000e+01,
1242        1.80000000e+01,   1.80500000e+01,   1.81000000e+01,   1.81500000e+01,
1243        1.82000000e+01,   1.82500000e+01,   1.83000000e+01,   1.83500000e+01,
1244        1.84000000e+01,   1.84500000e+01,   1.85000000e+01,   1.85500000e+01,
1245        1.86000000e+01,   1.86500000e+01,   1.87000000e+01,   1.87500000e+01,
1246        1.88000000e+01,   1.88500000e+01,   1.89000000e+01,   1.89500000e+01,
1247        1.90000000e+01,   1.90500000e+01,   1.91000000e+01,   1.91500000e+01,
1248        1.92000000e+01,   1.92500000e+01,   1.93000000e+01,   1.93500000e+01,
1249        1.94000000e+01,   1.94500000e+01,   1.95000000e+01,   1.95500000e+01,
1250        1.96000000e+01,   1.96500000e+01,   1.97000000e+01,   1.97500000e+01,
1251        1.98000000e+01,   1.98500000e+01,   1.99000000e+01,   1.99500000e+01,
1252        2.00000000e+01,   2.00500000e+01,   2.01000000e+01,   2.01500000e+01,
1253        2.02000000e+01,   2.02500000e+01,   2.03000000e+01,   2.03500000e+01,
1254        2.04000000e+01,   2.04500000e+01,   2.05000000e+01,   2.05500000e+01,
1255        2.06000000e+01,   2.06500000e+01,   2.07000000e+01,   2.07500000e+01,
1256        2.08000000e+01,   2.08500000e+01,   2.09000000e+01,   2.09500000e+01,
1257        2.10000000e+01,   2.10500000e+01,   2.11000000e+01,   2.11500000e+01,
1258        2.12000000e+01,   2.12500000e+01,   2.13000000e+01,   2.13500000e+01,
1259        2.14000000e+01,   2.14500000e+01,   2.15000000e+01,   2.15500000e+01,
1260        2.16000000e+01,   2.16500000e+01,   2.17000000e+01,   2.17500000e+01,
1261        2.18000000e+01,   2.18500000e+01,   2.19000000e+01,   2.19500000e+01,
1262        2.20000000e+01,   2.20500000e+01,   2.21000000e+01,   2.21500000e+01,
1263        2.22000000e+01,   2.22500000e+01,   2.23000000e+01,   2.23500000e+01,
1264        2.24000000e+01,   2.24500000e+01,   2.25000000e+01])
1265
1266        #print 'Diff', time[1:] - time[:-1]
1267
1268        #Setup mesh used to represent fitted function
1269        a = [0.0, 0.0]
1270        b = [0.0, 2.0]
1271        c = [2.0, 0.0]
1272        d = [0.0, 4.0]
1273        e = [2.0, 2.0]
1274        f = [4.0, 0.0]
1275
1276        points = [a, b, c, d, e, f]
1277        #bac, bce, ecf, dbe
1278        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1279
1280
1281        #New datapoints where interpolated values are sought
1282        interpolation_points = [[ 0.0, 0.0],
1283                                [ 0.5, 0.5],
1284                                [ 0.7, 0.7],
1285                                [ 1.0, 0.5],
1286                                [ 2.0, 0.4],
1287                                [ 545354534, 4354354353]] # outside the mesh
1288
1289        #One quantity
1290        Q = zeros( (len(time),6), Float )
1291
1292        #Linear in time and space
1293        for i, t in enumerate(time):
1294            Q[i, :] = t*linear_function(points)
1295
1296        #Check interpolation of one quantity using interpolaton points)
1297        try:
1298            I = Interpolation_function(time, Q,
1299                                       vertex_coordinates = points,
1300                                       triangles = triangles, 
1301                                       interpolation_points = interpolation_points,
1302                                       verbose = False)
1303        except:
1304            pass
1305        else:
1306            raise 'Should raise exception due to time being non-monotoneous'           
1307     
1308
1309    def test_points_outside_the_polygon(self):
1310        a = [-1.0, 0.0]
1311        b = [3.0, 4.0]
1312        c = [4.0, 1.0]
1313        d = [-3.0, 2.0] #3
1314        e = [-1.0, -2.0]
1315        f = [1.0, -2.0] #5
1316
1317        vertices = [a, b, c, d,e,f]
1318        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
1319
1320        point_coords = [[-2.0, 2.0],
1321                        [-1.0, 1.0],
1322                        [9999.0, 9999.0], # point Outside poly
1323                        [-9999.0, 1.0], # point Outside poly
1324                        [2.0, 1.0],
1325                        [0.0, 0.0],
1326                        [1.0, 0.0],
1327                        [0.0, -1.0],
1328                        [-0.2, -0.5],
1329                        [-0.9, -1.5],
1330                        [0.5, -1.9],
1331                        [999999, 9999999]] # point Outside poly
1332        geo_data = Geospatial_data(data_points = point_coords)
1333
1334        interp = Interpolate(vertices, triangles)
1335        f = array([linear_function(vertices),2*linear_function(vertices) ])
1336        f = transpose(f)
1337        #print "f",f
1338        z = interp.interpolate(f, geo_data)
1339        #z = interp.interpolate(f, point_coords)
1340        answer = [linear_function(point_coords),
1341                  2*linear_function(point_coords) ]
1342        answer = transpose(answer)
1343        answer[2,:] = [INF, INF]
1344        answer[3,:] = [INF, INF]
1345        answer[11,:] = [INF, INF]
1346        #print "z",z
1347        #print "answer _ fixed",answer
1348        assert allclose(z[0:1], answer[0:1])
1349        assert allclose(z[4:10], answer[4:10])
1350        for i in [2,3,11]:
1351            self.failUnless( z[i,1] == answer[11,1], 'Fail!')
1352            self.failUnless( z[i,0] == answer[11,0], 'Fail!')
1353       
1354#-------------------------------------------------------------
1355if __name__ == "__main__":
1356
1357    suite = unittest.makeSuite(Test_Interpolate,'test')
1358    #suite = unittest.makeSuite(Test_Interpolate,'test_interpolation_precompute_points')
1359    runner = unittest.TextTestRunner(verbosity=1)
1360    runner.run(suite)
1361
1362
1363
1364
1365
Note: See TracBrowser for help on using the repository browser.