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

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

Moved shallow water out from the old pyvolution directory.
All tests pass, most examples and okushiri works again.
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 anuga.shallow_water import Domain, Transmissive_boundary
23from anuga.utilities.numerical_tools import mean, NAN
24from anuga.shallow_water.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.