source: branches/inundation-numpy-branch/fit_interpolate/test_interpolate.py @ 7248

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

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

\anuga_core\source\anuga

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

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

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

Cheers
Duncan

File size: 24.6 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 anuga.pyvolution.util import mean
23from anuga.pyvolution.data_manager import get_dataobject
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
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
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,
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, 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
299        #results = interp._build_interpolation_matrix_A(data).todense()
300
301        assert allclose(A, answer)
302
303
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]
308
309        vertices = [v0, v1, v2]
310        triangles = [ [1,0,2] ]   #bac
311
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)
571        assert allclose(interp._A_can_be_reused, True)
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)
579        assert allclose(interp._A_can_be_reused, False)
580
581        #A is recalculated
582        z = interp.interpolate(f)
583        assert allclose(z, answer)
584        assert allclose(interp._A_can_be_reused, True)
585       
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
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)
625        from anuga.pyvolution.util import mean       
626        I = Interpolation_interface(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_interface(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_interface(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 BADtest_interpolate_sww(self):
797        """Not a unit test, rather a system test for interpolate_sww
798        This function is obsolete
799        """
800       
801        self.domain.filename = 'datatest' + str(time.time())
802        self.domain.format = 'sww'
803        self.domain.smooth = True
804        self.domain.reduction = mean
805
806        sww = get_dataobject(self.domain)
807        sww.store_connectivity()
808        sww.store_timestep('stage')
809        self.domain.time = 2.
810        sww.store_timestep('stage')
811
812        #print "self.domain.filename",self.domain.filename
813        interp = interpolate_sww(sww.filename, [0.0, 2.0],
814                                 [[0,1],[0.5,0.5]],
815                                 ['stage'])
816        #assert allclose(interp,[0.0,2.0])
817
818        #Cleanup
819        os.remove(sww.filename)
820       
821#-------------------------------------------------------------
822if __name__ == "__main__":
823    suite = unittest.makeSuite(Test_Interpolate,'test')
824    runner = unittest.TextTestRunner(verbosity=1)
825    runner.run(suite)
826
827
828
829
830
Note: See TracBrowser for help on using the repository browser.