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

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

Added exception and helpful diagnostics if a boundary condition tries
to access file_boundary outside its area defined by the underlying
sww file.

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