source: anuga_core/source/anuga/fit_interpolate/test_interpolate.py @ 3560

Last change on this file since 3560 was 3560, checked in by ole, 18 years ago

Renamed pyvolution to abstract_2d_finite_volumes. This is
one step towards pulling pyvolution apart. More to follow.
All unit tests pass and most examples fixed up.




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