source: inundation/fit_interpolate/test_interpolate.py @ 2563

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

get interpolate tests running in test_all

File size: 24.5 KB
RevLine 
[2187]1#!/usr/bin/env python
2
3#TEST
[2394]4
5#import time, os
6
7
[2187]8import sys
[2394]9import os
[2187]10import unittest
11from math import sqrt
[2394]12import tempfile
[2187]13
[2394]14from Scientific.IO.NetCDF import NetCDFFile
15from Numeric import allclose, array, transpose, zeros, Float
[2187]16
[2394]17
18# ANUGA code imports
[2187]19from interpolate import *
20from coordinate_transforms.geo_reference import Geo_reference
[2394]21from shallow_water import Domain, Transmissive_boundary #, mean
[2563]22from utilities.numerical_tools import mean
[2394]23from data_manager import get_dataobject
[2187]24
25def distance(x, y):
26    return sqrt( sum( (array(x)-array(y))**2 ))
27
28def linear_function(point):
29    point = array(point)
30    return point[:,0]+point[:,1]
31
32
33class Test_Interpolate(unittest.TestCase):
34
35    def setUp(self):
36
[2394]37        import time
38        from mesh_factory import rectangular
39
40
41        #Create basic mesh
42        points, vertices, boundary = rectangular(2, 2)
43
44        #Create shallow water domain
45        domain = Domain(points, vertices, boundary)
46        domain.default_order=2
47        domain.beta_h = 0
48
49
50        #Set some field values
51        domain.set_quantity('elevation', lambda x,y: -x)
52        domain.set_quantity('friction', 0.03)
53
54
55        ######################
56        # Boundary conditions
57        B = Transmissive_boundary(domain)
58        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
59
60
61        ######################
62        #Initial condition - with jumps
63
64        bed = domain.quantities['elevation'].vertex_values
65        stage = zeros(bed.shape, Float)
66
67        h = 0.3
68        for i in range(stage.shape[0]):
69            if i % 2 == 0:
70                stage[i,:] = bed[i,:] + h
71            else:
72                stage[i,:] = bed[i,:]
73
74        domain.set_quantity('stage', stage)
75
76        domain.distribute_to_vertices_and_edges()
77
78
79        self.domain = domain
80
81        C = domain.get_vertex_coordinates()
82        self.X = C[:,0:6:2].copy()
83        self.Y = C[:,1:6:2].copy()
84
85        self.F = bed
86
87
88
[2187]89    def tearDown(self):
90        pass
91
92    def test_datapoint_at_centroid(self):
93        a = [0.0, 0.0]
94        b = [0.0, 2.0]
95        c = [2.0,0.0]
96        points = [a, b, c]
97        vertices = [ [1,0,2] ]   #bac
98
99        data = [ [2.0/3, 2.0/3] ] #Use centroid as one data point
100
101        interp = Interpolate(points, vertices)
102        assert allclose(interp._build_interpolation_matrix_A(data).todense(),
103                        [[1./3, 1./3, 1./3]])
104
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,
[2201]131                               max_vertices_per_cell = 4)
[2187]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, data)
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, data)
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, data)
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        #FIXME - have to change this test to check default info
240    def NO_test_arbitrary_datapoints_some_outside(self):
241        """Try arbitrary datapoints one outside the triangle.
242        That one should be ignored
243        """
244
245        from Numeric import sum
246
247        a = [0.0, 0.0]
248        b = [0.0, 2.0]
249        c = [2.0,0.0]
250        points = [a, b, c]
251        vertices = [ [1,0,2] ]   #bac
252
253        data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44], [5.0, 7.0]]
254
255
256        interp = Interpolate(points, vertices, data, precrop = True)
257       
258        results = interp._build_interpolation_matrix_A(data).todense()
259        assert allclose(sum(results, axis=1), 1.0)
260
261        interp = Interpolate(points, vertices, data, precrop = False)
262        results = interp._build_interpolation_matrix_A(data).todense()
263        assert allclose(sum(results, axis=1), [1,1,1,0])
264
265
266
267    # this causes a memory error in scipy.sparse
268    def test_more_triangles(self):
269
270        a = [-1.0, 0.0]
271        b = [3.0, 4.0]
272        c = [4.0,1.0]
273        d = [-3.0, 2.0] #3
274        e = [-1.0,-2.0]
275        f = [1.0, -2.0] #5
276
277        points = [a, b, c, d,e,f]
278        triangles = [[0,1,3],[1,0,2],[0,4,5], [0,5,2]] #abd bac aef afc
279
280        #Data points
281        data = [ [-3., 2.0], [-2, 1], [0.0, 1], [0, 3], [2, 3], [-1.0/3,-4./3] ]
282        interp = Interpolate(points, triangles)
283
284        answer = [[0.0, 0.0, 0.0, 1.0, 0.0, 0.0],    #Affects point d
285                  [0.5, 0.0, 0.0, 0.5, 0.0, 0.0],    #Affects points a and d
286                  [0.75, 0.25, 0.0, 0.0, 0.0, 0.0],  #Affects points a and b
287                  [0.0, 0.5, 0.0, 0.5, 0.0, 0.0],    #Affects points a and d
288                  [0.25, 0.75, 0.0, 0.0, 0.0, 0.0],  #Affects points a and b
289                  [1./3, 0.0, 0.0, 0.0, 1./3, 1./3]] #Affects points a, e and f
290
291
292        A = interp._build_interpolation_matrix_A(data).todense()
293        for i in range(A.shape[0]):
294            for j in range(A.shape[1]):
295                if not allclose(A[i,j], answer[i][j]):
296                    print i,j,':',A[i,j], answer[i][j]
297
298
[2190]299        #results = interp._build_interpolation_matrix_A(data).todense()
[2187]300
[2190]301        assert allclose(A, answer)
[2187]302
[2190]303
[2189]304    def test_interpolate_attributes_to_points(self):
305        v0 = [0.0, 0.0]
306        v1 = [0.0, 5.0]
307        v2 = [5.0, 0.0]
[2187]308
[2189]309        vertices = [v0, v1, v2]
310        triangles = [ [1,0,2] ]   #bac
[2187]311
[2189]312        d0 = [1.0, 1.0]
313        d1 = [1.0, 2.0]
314        d2 = [3.0, 1.0]
315        point_coords = [ d0, d1, d2]
316
317        interp = Interpolate(vertices, triangles, point_coords)
318        f = linear_function(vertices)
319        z = interp.interpolate(f, point_coords)
320        answer = linear_function(point_coords)
321
322
323        assert allclose(z, answer)
324
325
326
327    def test_interpolate_attributes_to_pointsII(self):
328        a = [-1.0, 0.0]
329        b = [3.0, 4.0]
330        c = [4.0, 1.0]
331        d = [-3.0, 2.0] #3
332        e = [-1.0, -2.0]
333        f = [1.0, -2.0] #5
334
335        vertices = [a, b, c, d,e,f]
336        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
337
338
339        point_coords = [[-2.0, 2.0],
340                        [-1.0, 1.0],
341                        [0.0, 2.0],
342                        [1.0, 1.0],
343                        [2.0, 1.0],
344                        [0.0, 0.0],
345                        [1.0, 0.0],
346                        [0.0, -1.0],
347                        [-0.2, -0.5],
348                        [-0.9, -1.5],
349                        [0.5, -1.9],
350                        [3.0, 1.0]]
351
352        interp = Interpolate(vertices, triangles)
353        f = linear_function(vertices)
354        z = interp.interpolate(f, point_coords)
355        answer = linear_function(point_coords)
356        #print "z",z
357        #print "answer",answer
358        assert allclose(z, answer)
359
360    def test_interpolate_attributes_to_pointsIII(self):
361        """Test linear interpolation of known values at vertices to
362        new points inside a triangle
363        """
364        a = [0.0, 0.0]
365        b = [0.0, 5.0]
366        c = [5.0, 0.0]
367        d = [5.0, 5.0]
368
369        vertices = [a, b, c, d]
370        triangles = [ [1,0,2], [2,3,1] ]   #bac, cdb
371
372        #Points within triangle 1
373        d0 = [1.0, 1.0]
374        d1 = [1.0, 2.0]
375        d2 = [3.0, 1.0]
376
377        #Point within triangle 2
378        d3 = [4.0, 3.0]
379
380        #Points on common edge
381        d4 = [2.5, 2.5]
382        d5 = [4.0, 1.0]
383
384        #Point on common vertex
385        d6 = [0., 5.]
386       
387        point_coords = [d0, d1, d2, d3, d4, d5, d6]
388
389        interp = Interpolate(vertices, triangles)
390
391        #Known values at vertices
392        #Functions are x+y, x+2y, 2x+y, x-y-5
393        f = [ [0., 0., 0., -5.],        # (0,0)
394              [5., 10., 5., -10.],      # (0,5)
395              [5., 5., 10.0, 0.],       # (5,0)
396              [10., 15., 15., -5.]]     # (5,5)
397
398        z = interp.interpolate(f, point_coords)
399        answer = [ [2., 3., 3., -5.],   # (1,1)
400                   [3., 5., 4., -6.],   # (1,2)
401                   [4., 5., 7., -3.],   # (3,1)
402                   [7., 10., 11., -4.], # (4,3)
403                   [5., 7.5, 7.5, -5.], # (2.5, 2.5)
404                   [5., 6., 9., -2.],   # (4,1)
405                   [5., 10., 5., -10.]]  # (0,5)
406
407        #print "***********"
408        #print "z",z
409        #print "answer",answer
410        #print "***********"
411
412        #Should an error message be returned if points are outside
413        # of the mesh? Not currently. 
414
415        assert allclose(z, answer)
416
417
418    def test_interpolate_point_outside_of_mesh(self):
419        """Test linear interpolation of known values at vertices to
420        new points inside a triangle
421        """
422        a = [0.0, 0.0]
423        b = [0.0, 5.0]
424        c = [5.0, 0.0]
425        d = [5.0, 5.0]
426
427        vertices = [a, b, c, d]
428        triangles = [ [1,0,2], [2,3,1] ]   #bac, cdb
429
430        #Far away point
431        d7 = [-1., -1.]
432       
433        point_coords = [ d7]
434        interp = Interpolate(vertices, triangles)
435
436        #Known values at vertices
437        #Functions are x+y, x+2y, 2x+y, x-y-5
438        f = [ [0., 0., 0., -5.],        # (0,0)
439              [5., 10., 5., -10.],      # (0,5)
440              [5., 5., 10.0, 0.],       # (5,0)
441              [10., 15., 15., -5.]]     # (5,5)
442
443        z = interp.interpolate(f, point_coords)
444        answer = [ [0., 0., 0., 0.]] # (-1,-1)
445
446        #print "***********"
447        #print "z",z
448        #print "answer",answer
449        #print "***********"
450
451        #Should an error message be returned if points are outside
452        # of the mesh? Not currently. 
453
454        assert allclose(z, answer)
455
456    def test_interpolate_attributes_to_pointsIV(self):
457        a = [-1.0, 0.0]
458        b = [3.0, 4.0]
459        c = [4.0, 1.0]
460        d = [-3.0, 2.0] #3
461        e = [-1.0, -2.0]
462        f = [1.0, -2.0] #5
463
464        vertices = [a, b, c, d,e,f]
465        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
466
467
468        point_coords = [[-2.0, 2.0],
469                        [-1.0, 1.0],
470                        [0.0, 2.0],
471                        [1.0, 1.0],
472                        [2.0, 1.0],
473                        [0.0, 0.0],
474                        [1.0, 0.0],
475                        [0.0, -1.0],
476                        [-0.2, -0.5],
477                        [-0.9, -1.5],
478                        [0.5, -1.9],
479                        [3.0, 1.0]]
480
481        interp = Interpolate(vertices, triangles)
482        f = array([linear_function(vertices),2*linear_function(vertices) ])
483        f = transpose(f)
484        #print "f",f
485        z = interp.interpolate(f, point_coords)
486        answer = [linear_function(point_coords),
487                  2*linear_function(point_coords) ]
488        answer = transpose(answer)
489        #print "z",z
490        #print "answer",answer
491        assert allclose(z, answer)
492
493
494    def test_interpolate_blocking(self):
495        a = [-1.0, 0.0]
496        b = [3.0, 4.0]
497        c = [4.0, 1.0]
498        d = [-3.0, 2.0] #3
499        e = [-1.0, -2.0]
500        f = [1.0, -2.0] #5
501
502        vertices = [a, b, c, d,e,f]
503        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
504
505
506        point_coords = [[-2.0, 2.0],
507                        [-1.0, 1.0],
508                        [0.0, 2.0],
509                        [1.0, 1.0],
510                        [2.0, 1.0],
511                        [0.0, 0.0],
512                        [1.0, 0.0],
513                        [0.0, -1.0],
514                        [-0.2, -0.5],
515                        [-0.9, -1.5],
516                        [0.5, -1.9],
517                        [3.0, 1.0]]
518
519        interp = Interpolate(vertices, triangles)
520        f = array([linear_function(vertices),2*linear_function(vertices) ])
521        f = transpose(f)
522        #print "f",f
523        for blocking_max in range(len(point_coords)+2):
524        #if True:
525         #   blocking_max = 5
526            z = interp.interpolate(f, point_coords,
527                                   start_blocking_len=blocking_max)
528            answer = [linear_function(point_coords),
529                      2*linear_function(point_coords) ]
530            answer = transpose(answer)
531            #print "z",z
532            #print "answer",answer
533            assert allclose(z, answer)
534
535    def test_interpolate_reuse(self):
536        a = [-1.0, 0.0]
537        b = [3.0, 4.0]
538        c = [4.0, 1.0]
539        d = [-3.0, 2.0] #3
540        e = [-1.0, -2.0]
541        f = [1.0, -2.0] #5
542
543        vertices = [a, b, c, d,e,f]
544        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
545
546
547        point_coords = [[-2.0, 2.0],
548                        [-1.0, 1.0],
549                        [0.0, 2.0],
550                        [1.0, 1.0],
551                        [2.0, 1.0],
552                        [0.0, 0.0],
553                        [1.0, 0.0],
554                        [0.0, -1.0],
555                        [-0.2, -0.5],
556                        [-0.9, -1.5],
557                        [0.5, -1.9],
558                        [3.0, 1.0]]
559
560        interp = Interpolate(vertices, triangles)
561        f = array([linear_function(vertices),2*linear_function(vertices) ])
562        f = transpose(f)
563        z = interp.interpolate(f, point_coords,
564                               start_blocking_len=20)
565        answer = [linear_function(point_coords),
566                  2*linear_function(point_coords) ]
567        answer = transpose(answer)
568        #print "z",z
569        #print "answer",answer
570        assert allclose(z, answer)
[2201]571        assert allclose(interp._A_can_be_reused, True)
[2189]572
573        z = interp.interpolate(f)
574        assert allclose(z, answer)
575       
576        # This causes blocking to occur.
577        z = interp.interpolate(f, start_blocking_len=10)
578        assert allclose(z, answer)
[2201]579        assert allclose(interp._A_can_be_reused, False)
[2189]580
581        #A is recalculated
582        z = interp.interpolate(f)
583        assert allclose(z, answer)
[2201]584        assert allclose(interp._A_can_be_reused, True)
585       
[2189]586        interp = Interpolate(vertices, triangles)
587        #Must raise an exception, no points specified
588        try:
589            z = interp.interpolate(f)
590        except:
591            pass
592       
593
[2394]594
595    def test_interpolation_interface_time_only(self):
596        """Test spatio-temporal interpolation
597        Test that spatio temporal function performs the correct
598        interpolations in both time and space
599        """
600
601
602        #Three timesteps
603        time = [1.0, 5.0, 6.0]
604       
605
606        #One quantity
607        Q = zeros( (3,6), Float )
608
609        #Linear in time and space
610        a = [0.0, 0.0]
611        b = [0.0, 2.0]
612        c = [2.0, 0.0]
613        d = [0.0, 4.0]
614        e = [2.0, 2.0]
615        f = [4.0, 0.0]
616
617        points = [a, b, c, d, e, f]
618       
619        for i, t in enumerate(time):
620            Q[i, :] = t*linear_function(points)
621
622           
623        #Check basic interpolation of one quantity using averaging
624        #(no interpolation points or spatial info)
[2563]625        I = Interpolation_function(time, [mean(Q[0,:]),
[2394]626                                          mean(Q[1,:]),
627                                          mean(Q[2,:])])
628
629
630
631        #Check temporal interpolation
632        for i in [0,1,2]:
633            assert allclose(I(time[i]), mean(Q[i,:]))
634
635        #Midway   
636        assert allclose(I( (time[0] + time[1])/2 ),
637                        (I(time[0]) + I(time[1]))/2 )
638
639        assert allclose(I( (time[1] + time[2])/2 ),
640                        (I(time[1]) + I(time[2]))/2 )
641
642        assert allclose(I( (time[0] + time[2])/2 ),
643                        (I(time[0]) + I(time[2]))/2 )                 
644
645        #1/3
646        assert allclose(I( (time[0] + time[2])/3 ),
647                        (I(time[0]) + I(time[2]))/3 )                         
648
649
650        #Out of bounds checks
651        try:
652            I(time[0]-1) 
653        except:
654            pass
655        else:
656            raise 'Should raise exception'
657
658        try:
659            I(time[-1]+1) 
660        except:
661            pass
662        else:
663            raise 'Should raise exception'       
664
665
666       
667
668    def test_interpolation_interface_spatial_only(self):
669        """Test spatio-temporal interpolation with constant time
670        """
671
672        #Three timesteps
673        time = [1.0, 5.0, 6.0]
674       
675       
676        #Setup mesh used to represent fitted function
677        a = [0.0, 0.0]
678        b = [0.0, 2.0]
679        c = [2.0, 0.0]
680        d = [0.0, 4.0]
681        e = [2.0, 2.0]
682        f = [4.0, 0.0]
683
684        points = [a, b, c, d, e, f]
685        #bac, bce, ecf, dbe
686        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
687
688
689        #New datapoints where interpolated values are sought
690        interpolation_points = [[ 0.0, 0.0],
691                                [ 0.5, 0.5],
692                                [ 0.7, 0.7],
693                                [ 1.0, 0.5],
694                                [ 2.0, 0.4],
695                                [ 2.8, 1.2]]
696
697
698        #One quantity linear in space
699        Q = linear_function(points)
700
701
702        #Check interpolation of one quantity using interpolaton points
[2563]703        I = Interpolation_function(time, Q,
[2394]704                                   vertex_coordinates = points,
705                                   triangles = triangles, 
706                                   interpolation_points = interpolation_points,
707                                   verbose = False)
708
709
710        answer = linear_function(interpolation_points)
711
712        t = time[0]
713        for j in range(50): #t in [1, 6]
714            for id in range(len(interpolation_points)):
715                assert allclose(I(t, id), answer[id])
716
717            t += 0.1   
718
719
720        try:   
721            I(1)
722        except:
723            pass
724        else:
725            raise 'Should raise exception'
726           
727
728
729    def test_interpolation_interface(self):
730        """Test spatio-temporal interpolation
731        Test that spatio temporal function performs the correct
732        interpolations in both time and space
733        """
734
735
736        #Three timesteps
737        time = [1.0, 5.0, 6.0]
738       
739
740        #Setup mesh used to represent fitted function
741        a = [0.0, 0.0]
742        b = [0.0, 2.0]
743        c = [2.0, 0.0]
744        d = [0.0, 4.0]
745        e = [2.0, 2.0]
746        f = [4.0, 0.0]
747
748        points = [a, b, c, d, e, f]
749        #bac, bce, ecf, dbe
750        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
751
752
753        #New datapoints where interpolated values are sought
754        interpolation_points = [[ 0.0, 0.0],
755                                [ 0.5, 0.5],
756                                [ 0.7, 0.7],
757                                [ 1.0, 0.5],
758                                [ 2.0, 0.4],
759                                [ 2.8, 1.2]]
760
761
762        #One quantity
763        Q = zeros( (3,6), Float )
764
765        #Linear in time and space
766        for i, t in enumerate(time):
767            Q[i, :] = t*linear_function(points)
768
769
770        #Check interpolation of one quantity using interpolaton points)
[2563]771        I = Interpolation_function(time, Q,
[2394]772                                   vertex_coordinates = points,
773                                   triangles = triangles, 
774                                   interpolation_points = interpolation_points,
775                                   verbose = False)
776
777
778        answer = linear_function(interpolation_points)
779       
780        t = time[0]
781        for j in range(50): #t in [1, 6]
782            for id in range(len(interpolation_points)):
783                assert allclose(I(t, id), t*answer[id])
784
785            t += 0.1   
786
787        try:   
788            I(1)
789        except:
790            pass
791        else:
792            raise 'Should raise exception'
793
794       
795    def BADtest_interpolate_sww(self):
796        """Not a unit test, rather a system test for interpolate_sww
797        This function is obsolete
798        """
799       
800        self.domain.filename = 'datatest' + str(time.time())
801        self.domain.format = 'sww'
802        self.domain.smooth = True
803        self.domain.reduction = mean
804
805        sww = get_dataobject(self.domain)
806        sww.store_connectivity()
807        sww.store_timestep('stage')
808        self.domain.time = 2.
809        sww.store_timestep('stage')
810
811        #print "self.domain.filename",self.domain.filename
812        interp = interpolate_sww(sww.filename, [0.0, 2.0],
813                                 [[0,1],[0.5,0.5]],
814                                 ['stage'])
815        #assert allclose(interp,[0.0,2.0])
816
817        #Cleanup
818        os.remove(sww.filename)
819       
[2187]820#-------------------------------------------------------------
821if __name__ == "__main__":
822    suite = unittest.makeSuite(Test_Interpolate,'test')
823    runner = unittest.TextTestRunner(verbosity=1)
824    runner.run(suite)
825
826
827
828
829
Note: See TracBrowser for help on using the repository browser.