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

Last change on this file since 3900 was 3900, checked in by ole, 17 years ago

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