source: inundation/fit_interpolate/test_interpolate.py @ 2577

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

changing interpolate to use geospatial_data objects
File size: 25.0 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 #, mean
22from utilities.numerical_tools import mean
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
107    def test_quad_tree(self):
108        p0 = [-10.0, -10.0]
109        p1 = [20.0, -10.0]
110        p2 = [-10.0, 20.0]
111        p3 = [10.0, 50.0]
112        p4 = [30.0, 30.0]
113        p5 = [50.0, 10.0]
114        p6 = [40.0, 60.0]
115        p7 = [60.0, 40.0]
116        p8 = [-66.0, 20.0]
117        p9 = [10.0, -66.0]
118
119        points = [p0, p1, p2, p3, p4, p5, p6, p7, p8, p9]
120        triangles = [ [0, 1, 2],
121                      [3, 2, 4],
122                      [4, 2, 1],
123                      [4, 1, 5],
124                      [3, 4, 6],
125                      [6, 4, 7],
126                      [7, 4, 5],
127                      [8, 0, 2],
128                      [0, 9, 1]]
129
130        data = [ [4,4] ]
131        interp = Interpolate(points, triangles,
132                               max_vertices_per_cell = 4)
133        #print "PDSG - interp.get_A()", interp.get_A()
134        answer =  [ [ 0.06666667,  0.46666667,  0.46666667,  0.,
135                      0., 0. , 0., 0., 0., 0.]]
136
137       
138        assert allclose(interp._build_interpolation_matrix_A(data).todense(),
139                        answer)
140        #interp.set_point_coordinates([[-30, -30]]) #point outside of mesh
141        #print "PDSG - interp.get_A()", interp.get_A()
142        data = [[-30, -30]]
143        answer =  [ [ 0.0,  0.0,  0.0,  0.,
144                      0., 0. , 0., 0., 0., 0.]]
145       
146        assert allclose(interp._build_interpolation_matrix_A(data).todense(),
147                        answer)
148
149
150        #point outside of quad tree root cell
151        #interp.set_point_coordinates([[-70, -70]])
152        #print "PDSG - interp.get_A()", interp.get_A()
153        data = [[-70, -70]]
154        answer =  [ [ 0.0,  0.0,  0.0,  0.,
155                      0., 0. , 0., 0., 0., 0.]]
156        assert allclose(interp._build_interpolation_matrix_A(data).todense(),
157                        answer)
158
159    def test_datapoints_at_vertices(self):
160        """Test that data points coinciding with vertices yield a diagonal matrix
161        """
162
163        a = [0.0, 0.0]
164        b = [0.0, 2.0]
165        c = [2.0,0.0]
166        points = [a, b, c]
167        vertices = [ [1,0,2] ]   #bac
168
169        data = points #Use data at vertices
170
171        interp = Interpolate(points, vertices)
172        answer = [[1., 0., 0.],
173                   [0., 1., 0.],
174                   [0., 0., 1.]]
175        assert allclose(interp._build_interpolation_matrix_A(data).todense(),
176                        answer)
177
178
179    def test_datapoints_on_edge_midpoints(self):
180        """Try datapoints midway on edges -
181        each point should affect two matrix entries equally
182        """
183
184        a = [0.0, 0.0]
185        b = [0.0, 2.0]
186        c = [2.0,0.0]
187        points = [a, b, c]
188        vertices = [ [1,0,2] ]   #bac
189
190        data = [ [0., 1.], [1., 0.], [1., 1.] ]
191        answer =  [[0.5, 0.5, 0.0],  #Affects vertex 1 and 0
192                    [0.5, 0.0, 0.5],  #Affects vertex 0 and 2
193                    [0.0, 0.5, 0.5]]
194        interp = Interpolate(points, vertices, data)
195
196        assert allclose(interp._build_interpolation_matrix_A(data).todense(),
197                        answer)
198
199    def test_datapoints_on_edges(self):
200        """Try datapoints on edges -
201        each point should affect two matrix entries in proportion
202        """
203
204        a = [0.0, 0.0]
205        b = [0.0, 2.0]
206        c = [2.0,0.0]
207        points = [a, b, c]
208        vertices = [ [1,0,2] ]   #bac
209
210        data = [ [0., 1.5], [1.5, 0.], [1.5, 0.5] ]
211        answer =  [[0.25, 0.75, 0.0],  #Affects vertex 1 and 0
212                   [0.25, 0.0, 0.75],  #Affects vertex 0 and 2
213                   [0.0, 0.25, 0.75]]
214
215        interp = Interpolate(points, vertices, data)
216
217        assert allclose(interp._build_interpolation_matrix_A(data).todense(),
218                        answer)
219
220
221    def test_arbitrary_datapoints(self):
222        """Try arbitrary datapoints
223        """
224
225        from Numeric import sum
226
227        a = [0.0, 0.0]
228        b = [0.0, 2.0]
229        c = [2.0,0.0]
230        points = [a, b, c]
231        vertices = [ [1,0,2] ]   #bac
232
233        data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44] ]
234
235        interp = Interpolate(points, vertices, data)
236        #print "interp.get_A()", interp.get_A()
237        results = interp._build_interpolation_matrix_A(data).todense()
238        assert allclose(sum(results, axis=1), 1.0)
239
240        #FIXME - have to change this test to check default info
241    def NO_test_arbitrary_datapoints_some_outside(self):
242        """Try arbitrary datapoints one outside the triangle.
243        That one should be ignored
244        """
245
246        from Numeric import sum
247
248        a = [0.0, 0.0]
249        b = [0.0, 2.0]
250        c = [2.0,0.0]
251        points = [a, b, c]
252        vertices = [ [1,0,2] ]   #bac
253
254        data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44], [5.0, 7.0]]
255
256
257        interp = Interpolate(points, vertices, data, precrop = True)
258       
259        results = interp._build_interpolation_matrix_A(data).todense()
260        assert allclose(sum(results, axis=1), 1.0)
261
262        interp = Interpolate(points, vertices, data, precrop = False)
263        results = interp._build_interpolation_matrix_A(data).todense()
264        assert allclose(sum(results, axis=1), [1,1,1,0])
265
266
267
268    # this causes a memory error in scipy.sparse
269    def test_more_triangles(self):
270
271        a = [-1.0, 0.0]
272        b = [3.0, 4.0]
273        c = [4.0,1.0]
274        d = [-3.0, 2.0] #3
275        e = [-1.0,-2.0]
276        f = [1.0, -2.0] #5
277
278        points = [a, b, c, d,e,f]
279        triangles = [[0,1,3],[1,0,2],[0,4,5], [0,5,2]] #abd bac aef afc
280
281        #Data points
282        data = [ [-3., 2.0], [-2, 1], [0.0, 1], [0, 3], [2, 3], [-1.0/3,-4./3] ]
283        interp = Interpolate(points, triangles)
284
285        answer = [[0.0, 0.0, 0.0, 1.0, 0.0, 0.0],    #Affects point d
286                  [0.5, 0.0, 0.0, 0.5, 0.0, 0.0],    #Affects points a and d
287                  [0.75, 0.25, 0.0, 0.0, 0.0, 0.0],  #Affects points a and b
288                  [0.0, 0.5, 0.0, 0.5, 0.0, 0.0],    #Affects points a and d
289                  [0.25, 0.75, 0.0, 0.0, 0.0, 0.0],  #Affects points a and b
290                  [1./3, 0.0, 0.0, 0.0, 1./3, 1./3]] #Affects points a, e and f
291
292
293        A = interp._build_interpolation_matrix_A(data).todense()
294        for i in range(A.shape[0]):
295            for j in range(A.shape[1]):
296                if not allclose(A[i,j], answer[i][j]):
297                    print i,j,':',A[i,j], answer[i][j]
298
299
300        #results = interp._build_interpolation_matrix_A(data).todense()
301
302        assert allclose(A, answer)
303
304
305    def test_interpolate_attributes_to_points(self):
306        v0 = [0.0, 0.0]
307        v1 = [0.0, 5.0]
308        v2 = [5.0, 0.0]
309
310        vertices = [v0, v1, v2]
311        triangles = [ [1,0,2] ]   #bac
312
313        d0 = [1.0, 1.0]
314        d1 = [1.0, 2.0]
315        d2 = [3.0, 1.0]
316        point_coords = [ d0, d1, d2]
317
318        interp = Interpolate(vertices, triangles, point_coords)
319        f = linear_function(vertices)
320        z = interp.interpolate(f, point_coords)
321        answer = linear_function(point_coords)
322
323
324        assert allclose(z, answer)
325
326
327
328    def test_interpolate_attributes_to_pointsII(self):
329        a = [-1.0, 0.0]
330        b = [3.0, 4.0]
331        c = [4.0, 1.0]
332        d = [-3.0, 2.0] #3
333        e = [-1.0, -2.0]
334        f = [1.0, -2.0] #5
335
336        vertices = [a, b, c, d,e,f]
337        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
338
339
340        point_coords = [[-2.0, 2.0],
341                        [-1.0, 1.0],
342                        [0.0, 2.0],
343                        [1.0, 1.0],
344                        [2.0, 1.0],
345                        [0.0, 0.0],
346                        [1.0, 0.0],
347                        [0.0, -1.0],
348                        [-0.2, -0.5],
349                        [-0.9, -1.5],
350                        [0.5, -1.9],
351                        [3.0, 1.0]]
352
353        interp = Interpolate(vertices, triangles)
354        f = linear_function(vertices)
355        z = interp.interpolate(f, point_coords)
356        answer = linear_function(point_coords)
357        #print "z",z
358        #print "answer",answer
359        assert allclose(z, answer)
360
361    def test_interpolate_attributes_to_pointsIII(self):
362        """Test linear interpolation of known values at vertices to
363        new points inside a triangle
364        """
365        a = [0.0, 0.0]
366        b = [0.0, 5.0]
367        c = [5.0, 0.0]
368        d = [5.0, 5.0]
369
370        vertices = [a, b, c, d]
371        triangles = [ [1,0,2], [2,3,1] ]   #bac, cdb
372
373        #Points within triangle 1
374        d0 = [1.0, 1.0]
375        d1 = [1.0, 2.0]
376        d2 = [3.0, 1.0]
377
378        #Point within triangle 2
379        d3 = [4.0, 3.0]
380
381        #Points on common edge
382        d4 = [2.5, 2.5]
383        d5 = [4.0, 1.0]
384
385        #Point on common vertex
386        d6 = [0., 5.]
387       
388        point_coords = [d0, d1, d2, d3, d4, d5, d6]
389
390        interp = Interpolate(vertices, triangles)
391
392        #Known values at vertices
393        #Functions are x+y, x+2y, 2x+y, x-y-5
394        f = [ [0., 0., 0., -5.],        # (0,0)
395              [5., 10., 5., -10.],      # (0,5)
396              [5., 5., 10.0, 0.],       # (5,0)
397              [10., 15., 15., -5.]]     # (5,5)
398
399        z = interp.interpolate(f, point_coords)
400        answer = [ [2., 3., 3., -5.],   # (1,1)
401                   [3., 5., 4., -6.],   # (1,2)
402                   [4., 5., 7., -3.],   # (3,1)
403                   [7., 10., 11., -4.], # (4,3)
404                   [5., 7.5, 7.5, -5.], # (2.5, 2.5)
405                   [5., 6., 9., -2.],   # (4,1)
406                   [5., 10., 5., -10.]]  # (0,5)
407
408        #print "***********"
409        #print "z",z
410        #print "answer",answer
411        #print "***********"
412
413        #Should an error message be returned if points are outside
414        # of the mesh? Not currently. 
415
416        assert allclose(z, answer)
417
418
419    def test_interpolate_point_outside_of_mesh(self):
420        """Test linear interpolation of known values at vertices to
421        new points inside a triangle
422        """
423        a = [0.0, 0.0]
424        b = [0.0, 5.0]
425        c = [5.0, 0.0]
426        d = [5.0, 5.0]
427
428        vertices = [a, b, c, d]
429        triangles = [ [1,0,2], [2,3,1] ]   #bac, cdb
430
431        #Far away point
432        d7 = [-1., -1.]
433       
434        point_coords = [ d7]
435        interp = Interpolate(vertices, triangles)
436
437        #Known values at vertices
438        #Functions are x+y, x+2y, 2x+y, x-y-5
439        f = [ [0., 0., 0., -5.],        # (0,0)
440              [5., 10., 5., -10.],      # (0,5)
441              [5., 5., 10.0, 0.],       # (5,0)
442              [10., 15., 15., -5.]]     # (5,5)
443
444        z = interp.interpolate(f, point_coords)
445        answer = [ [0., 0., 0., 0.]] # (-1,-1)
446
447        #print "***********"
448        #print "z",z
449        #print "answer",answer
450        #print "***********"
451
452        #Should an error message be returned if points are outside
453        # of the mesh? Not currently. 
454
455        assert allclose(z, answer)
456
457    def test_interpolate_attributes_to_pointsIV(self):
458        a = [-1.0, 0.0]
459        b = [3.0, 4.0]
460        c = [4.0, 1.0]
461        d = [-3.0, 2.0] #3
462        e = [-1.0, -2.0]
463        f = [1.0, -2.0] #5
464
465        vertices = [a, b, c, d,e,f]
466        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
467
468
469        point_coords = [[-2.0, 2.0],
470                        [-1.0, 1.0],
471                        [0.0, 2.0],
472                        [1.0, 1.0],
473                        [2.0, 1.0],
474                        [0.0, 0.0],
475                        [1.0, 0.0],
476                        [0.0, -1.0],
477                        [-0.2, -0.5],
478                        [-0.9, -1.5],
479                        [0.5, -1.9],
480                        [3.0, 1.0]]
481
482        interp = Interpolate(vertices, triangles)
483        f = array([linear_function(vertices),2*linear_function(vertices) ])
484        f = transpose(f)
485        #print "f",f
486        z = interp.interpolate(f, point_coords)
487        answer = [linear_function(point_coords),
488                  2*linear_function(point_coords) ]
489        answer = transpose(answer)
490        #print "z",z
491        #print "answer",answer
492        assert allclose(z, answer)
493
494
495    def test_interpolate_blocking(self):
496        a = [-1.0, 0.0]
497        b = [3.0, 4.0]
498        c = [4.0, 1.0]
499        d = [-3.0, 2.0] #3
500        e = [-1.0, -2.0]
501        f = [1.0, -2.0] #5
502
503        vertices = [a, b, c, d,e,f]
504        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
505
506
507        point_coords = [[-2.0, 2.0],
508                        [-1.0, 1.0],
509                        [0.0, 2.0],
510                        [1.0, 1.0],
511                        [2.0, 1.0],
512                        [0.0, 0.0],
513                        [1.0, 0.0],
514                        [0.0, -1.0],
515                        [-0.2, -0.5],
516                        [-0.9, -1.5],
517                        [0.5, -1.9],
518                        [3.0, 1.0]]
519
520        interp = Interpolate(vertices, triangles)
521        f = array([linear_function(vertices),2*linear_function(vertices) ])
522        f = transpose(f)
523        #print "f",f
524        for blocking_max in range(len(point_coords)+2):
525        #if True:
526         #   blocking_max = 5
527            z = interp.interpolate(f, point_coords,
528                                   start_blocking_len=blocking_max)
529            answer = [linear_function(point_coords),
530                      2*linear_function(point_coords) ]
531            answer = transpose(answer)
532            #print "z",z
533            #print "answer",answer
534            assert allclose(z, answer)
535
536    def test_interpolate_reuse(self):
537        a = [-1.0, 0.0]
538        b = [3.0, 4.0]
539        c = [4.0, 1.0]
540        d = [-3.0, 2.0] #3
541        e = [-1.0, -2.0]
542        f = [1.0, -2.0] #5
543
544        vertices = [a, b, c, d,e,f]
545        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
546
547
548        point_coords = [[-2.0, 2.0],
549                        [-1.0, 1.0],
550                        [0.0, 2.0],
551                        [1.0, 1.0],
552                        [2.0, 1.0],
553                        [0.0, 0.0],
554                        [1.0, 0.0],
555                        [0.0, -1.0],
556                        [-0.2, -0.5],
557                        [-0.9, -1.5],
558                        [0.5, -1.9],
559                        [3.0, 1.0]]
560
561        interp = Interpolate(vertices, triangles)
562        f = array([linear_function(vertices),2*linear_function(vertices) ])
563        f = transpose(f)
564        z = interp.interpolate(f, point_coords,
565                               start_blocking_len=20)
566        answer = [linear_function(point_coords),
567                  2*linear_function(point_coords) ]
568        answer = transpose(answer)
569        #print "z",z
570        #print "answer",answer
571        assert allclose(z, answer)
572        assert allclose(interp._A_can_be_reused, True)
573
574        z = interp.interpolate(f)
575        assert allclose(z, answer)
576       
577        # This causes blocking to occur.
578        z = interp.interpolate(f, start_blocking_len=10)
579        assert allclose(z, answer)
580        assert allclose(interp._A_can_be_reused, False)
581
582        #A is recalculated
583        z = interp.interpolate(f)
584        assert allclose(z, answer)
585        assert allclose(interp._A_can_be_reused, True)
586       
587        interp = Interpolate(vertices, triangles)
588        #Must raise an exception, no points specified
589        try:
590            z = interp.interpolate(f)
591        except:
592            pass
593       
594
595
596    def test_interpolation_interface_time_only(self):
597        """Test spatio-temporal interpolation
598        Test that spatio temporal function performs the correct
599        interpolations in both time and space
600        """
601
602
603        #Three timesteps
604        time = [1.0, 5.0, 6.0]
605       
606
607        #One quantity
608        Q = zeros( (3,6), Float )
609
610        #Linear in time and space
611        a = [0.0, 0.0]
612        b = [0.0, 2.0]
613        c = [2.0, 0.0]
614        d = [0.0, 4.0]
615        e = [2.0, 2.0]
616        f = [4.0, 0.0]
617
618        points = [a, b, c, d, e, f]
619       
620        for i, t in enumerate(time):
621            Q[i, :] = t*linear_function(points)
622
623           
624        #Check basic interpolation of one quantity using averaging
625        #(no interpolation points or spatial info)
626        I = Interpolation_function(time, [mean(Q[0,:]),
627                                          mean(Q[1,:]),
628                                          mean(Q[2,:])])
629
630
631
632        #Check temporal interpolation
633        for i in [0,1,2]:
634            assert allclose(I(time[i]), mean(Q[i,:]))
635
636        #Midway   
637        assert allclose(I( (time[0] + time[1])/2 ),
638                        (I(time[0]) + I(time[1]))/2 )
639
640        assert allclose(I( (time[1] + time[2])/2 ),
641                        (I(time[1]) + I(time[2]))/2 )
642
643        assert allclose(I( (time[0] + time[2])/2 ),
644                        (I(time[0]) + I(time[2]))/2 )                 
645
646        #1/3
647        assert allclose(I( (time[0] + time[2])/3 ),
648                        (I(time[0]) + I(time[2]))/3 )                         
649
650
651        #Out of bounds checks
652        try:
653            I(time[0]-1) 
654        except:
655            pass
656        else:
657            raise 'Should raise exception'
658
659        try:
660            I(time[-1]+1) 
661        except:
662            pass
663        else:
664            raise 'Should raise exception'       
665
666
667       
668
669    def test_interpolation_interface_spatial_only(self):
670        """Test spatio-temporal interpolation with constant time
671        """
672
673        #Three timesteps
674        time = [1.0, 5.0, 6.0]
675       
676       
677        #Setup mesh used to represent fitted function
678        a = [0.0, 0.0]
679        b = [0.0, 2.0]
680        c = [2.0, 0.0]
681        d = [0.0, 4.0]
682        e = [2.0, 2.0]
683        f = [4.0, 0.0]
684
685        points = [a, b, c, d, e, f]
686        #bac, bce, ecf, dbe
687        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
688
689
690        #New datapoints where interpolated values are sought
691        interpolation_points = [[ 0.0, 0.0],
692                                [ 0.5, 0.5],
693                                [ 0.7, 0.7],
694                                [ 1.0, 0.5],
695                                [ 2.0, 0.4],
696                                [ 2.8, 1.2]]
697
698
699        #One quantity linear in space
700        Q = linear_function(points)
701
702
703        #Check interpolation of one quantity using interpolaton points
704        I = Interpolation_function(time, Q,
705                                   vertex_coordinates = points,
706                                   triangles = triangles, 
707                                   interpolation_points = interpolation_points,
708                                   verbose = False)
709
710
711        answer = linear_function(interpolation_points)
712
713        t = time[0]
714        for j in range(50): #t in [1, 6]
715            for id in range(len(interpolation_points)):
716                assert allclose(I(t, id), answer[id])
717
718            t += 0.1   
719
720
721        try:   
722            I(1)
723        except:
724            pass
725        else:
726            raise 'Should raise exception'
727           
728
729
730    def test_interpolation_interface(self):
731        """Test spatio-temporal interpolation
732        Test that spatio temporal function performs the correct
733        interpolations in both time and space
734        """
735
736
737        #Three timesteps
738        time = [1.0, 5.0, 6.0]
739       
740
741        #Setup mesh used to represent fitted function
742        a = [0.0, 0.0]
743        b = [0.0, 2.0]
744        c = [2.0, 0.0]
745        d = [0.0, 4.0]
746        e = [2.0, 2.0]
747        f = [4.0, 0.0]
748
749        points = [a, b, c, d, e, f]
750        #bac, bce, ecf, dbe
751        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
752
753
754        #New datapoints where interpolated values are sought
755        interpolation_points = [[ 0.0, 0.0],
756                                [ 0.5, 0.5],
757                                [ 0.7, 0.7],
758                                [ 1.0, 0.5],
759                                [ 2.0, 0.4],
760                                [ 2.8, 1.2]]
761
762
763        #One quantity
764        Q = zeros( (3,6), Float )
765
766        #Linear in time and space
767        for i, t in enumerate(time):
768            Q[i, :] = t*linear_function(points)
769
770
771        #Check interpolation of one quantity using interpolaton points)
772        I = Interpolation_function(time, Q,
773                                   vertex_coordinates = points,
774                                   triangles = triangles, 
775                                   interpolation_points = interpolation_points,
776                                   verbose = False)
777
778
779        answer = linear_function(interpolation_points)
780       
781        t = time[0]
782        for j in range(50): #t in [1, 6]
783            for id in range(len(interpolation_points)):
784                assert allclose(I(t, id), t*answer[id])
785
786            t += 0.1   
787
788        try:   
789            I(1)
790        except:
791            pass
792        else:
793            raise 'Should raise exception'
794
795
796    def qtest_interpolation_interface(self):
797        a = [-1.0, 0.0]
798        b = [3.0, 4.0]
799        c = [4.0, 1.0]
800        d = [-3.0, 2.0] #3
801        e = [-1.0, -2.0]
802        f = [1.0, -2.0] #5
803
804        vertices = [a, b, c, d,e,f]
805        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
806
807        point_coords = [[-2.0, 2.0],
808                        [-1.0, 1.0],
809                        [0.0, 2.0],
810                        [1.0, 1.0],
811                        [2.0, 1.0],
812                        [0.0, 0.0],
813                        [1.0, 0.0],
814                        [0.0, -1.0],
815                        [-0.2, -0.5],
816                        [-0.9, -1.5],
817                        [0.5, -1.9],
818                        [999999, 9999999]]
819        geo_data = Geospatial_data(data_points = point_coords)
820
821        interp = Interpolate(vertices, triangles)
822        f = array([linear_function(vertices),2*linear_function(vertices) ])
823        f = transpose(f)
824        #print "f",f
825        z = interp.interpolate(f, geo_data)
826        #z = interp.interpolate(f, point_coords)
827        answer = [linear_function(point_coords),
828                  2*linear_function(point_coords) ]
829        answer = transpose(answer)
830        #print "z",z
831        #print "answer",answer
832        assert allclose(z, answer)
833
834       
835#-------------------------------------------------------------
836if __name__ == "__main__":
837    suite = unittest.makeSuite(Test_Interpolate,'qtest')
838    runner = unittest.TextTestRunner(verbosity=1)
839    runner.run(suite)
840
841
842
843
844
Note: See TracBrowser for help on using the repository browser.