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

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

Added functionality for getting arbitrary interpolated values in Quantity as well as calculating inundation height and location. This work was done at SUT during the last week of September 2006.
File size: 56.2 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    def test_interpolation_precompute_points(self):
1129        # looking at a discrete mesh
1130        #
1131   
1132        #Three timesteps
1133        time = [0.0, 60.0]   
1134
1135        #Setup mesh used to represent fitted function
1136        points = [[ 15., -20.],
1137                  [ 15.,  10.],
1138                  [  0., -20.],
1139                  [  0.,  10.],
1140                  [  0., -20.],
1141                  [ 15.,  10.]]
1142       
1143        triangles = [[0, 1, 2],
1144                     [3, 4, 5]]
1145
1146        #New datapoints where interpolated values are sought
1147        interpolation_points = [[ 1.,  0.], [0.,1.]]
1148
1149        #One quantity
1150        Q = zeros( (2,6), Float )
1151
1152        #Linear in time and space
1153        for i, t in enumerate(time):
1154            Q[i, :] = t*linear_function(points)
1155        #print "Q", Q
1156
1157
1158       
1159        interp = Interpolate(points, triangles)
1160        f = array([linear_function(points),2*linear_function(points) ])
1161        f = transpose(f)
1162        #print "f",f
1163        z = interp.interpolate(f, interpolation_points)
1164        answer = [linear_function(interpolation_points),
1165                  2*linear_function(interpolation_points) ]
1166        answer = transpose(answer)
1167        #print "z",z
1168        #print "answer",answer
1169        assert allclose(z, answer)
1170
1171
1172        #Check interpolation of one quantity using interpolaton points)
1173        I = Interpolation_function(time, Q,
1174                                   vertex_coordinates = points,
1175                                   triangles = triangles, 
1176                                   interpolation_points = interpolation_points,
1177                                   verbose = False)
1178       
1179        #print "I.precomputed_values", I.precomputed_values
1180
1181        msg = 'Interpolation failed'
1182        assert allclose(I.precomputed_values['Attribute'][1], [60, 60]), msg
1183        #self.failUnless( I.precomputed_values['Attribute'][1] == 60.0,
1184        #                ' failed')
1185       
1186    def test_interpolation_function_outside_point(self):
1187        # Test spatio-temporal interpolation
1188        # Test that spatio temporal function performs the correct
1189        # interpolations in both time and space
1190   
1191        #Three timesteps
1192        time = [1.0, 5.0, 6.0]   
1193
1194        #Setup mesh used to represent fitted function
1195        a = [0.0, 0.0]
1196        b = [0.0, 2.0]
1197        c = [2.0, 0.0]
1198        d = [0.0, 4.0]
1199        e = [2.0, 2.0]
1200        f = [4.0, 0.0]
1201
1202        points = [a, b, c, d, e, f]
1203        #bac, bce, ecf, dbe
1204        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1205
1206
1207        #New datapoints where interpolated values are sought
1208        interpolation_points = [[ 0.0, 0.0],
1209                                [ 0.5, 0.5],
1210                                [ 0.7, 0.7],
1211                                [ 1.0, 0.5],
1212                                [ 2.0, 0.4],
1213                                [ 545354534, 4354354353]] # outside the mesh
1214
1215        #One quantity
1216        Q = zeros( (3,6), Float )
1217
1218        #Linear in time and space
1219        for i, t in enumerate(time):
1220            Q[i, :] = t*linear_function(points)
1221
1222        #Check interpolation of one quantity using interpolaton points)
1223        I = Interpolation_function(time, Q,
1224                                   vertex_coordinates = points,
1225                                   triangles = triangles, 
1226                                   interpolation_points = interpolation_points,
1227                                   verbose = False)
1228
1229        answer = linear_function(interpolation_points)
1230         
1231        t = time[0]
1232        for j in range(50): #t in [1, 6]
1233            for id in range(len(interpolation_points)-1):
1234                assert allclose(I(t, id), t*answer[id])
1235            t += 0.1
1236           
1237        # Now test the point outside the mesh
1238        t = time[0]
1239        for j in range(50): #t in [1, 6]
1240            self.failUnless(I(t, 5) == NAN, 'Fail!')
1241            t += 0.1 
1242           
1243        try:   
1244            I(1)
1245        except:
1246            pass
1247        else:
1248            raise 'Should raise exception'
1249
1250
1251    def test_interpolation_function_time(self):
1252        #Test a long time series with an error in it (this did cause an
1253        #error once)
1254       
1255
1256        time = array(\
1257        [0.00000000e+00, 5.00000000e-02, 1.00000000e-01,   1.50000000e-01,
1258        2.00000000e-01,   2.50000000e-01,   3.00000000e-01,   3.50000000e-01,
1259        4.00000000e-01,   4.50000000e-01,   5.00000000e-01,   5.50000000e-01,
1260        6.00000000e-01,   6.50000000e-01,   7.00000000e-01,   7.50000000e-01,
1261        8.00000000e-01,   8.50000000e-01,   9.00000000e-01,   9.50000000e-01,
1262        1.00000000e-00,   1.05000000e+00,   1.10000000e+00,   1.15000000e+00,
1263        1.20000000e+00,   1.25000000e+00,   1.30000000e+00,   1.35000000e+00,
1264        1.40000000e+00,   1.45000000e+00,   1.50000000e+00,   1.55000000e+00,
1265        1.60000000e+00,   1.65000000e+00,   1.70000000e+00,   1.75000000e+00,
1266        1.80000000e+00,   1.85000000e+00,   1.90000000e+00,   1.95000000e+00,
1267        2.00000000e+00,   2.05000000e+00,   2.10000000e+00,   2.15000000e+00,
1268        2.20000000e+00,   2.25000000e+00,   2.30000000e+00,   2.35000000e+00,
1269        2.40000000e+00,   2.45000000e+00,   2.50000000e+00,   2.55000000e+00,
1270        2.60000000e+00,   2.65000000e+00,   2.70000000e+00,   2.75000000e+00,
1271        2.80000000e+00,   2.85000000e+00,   2.90000000e+00,   2.95000000e+00,
1272        3.00000000e+00,   3.05000000e+00,   9.96920997e+36,   3.15000000e+00,
1273        3.20000000e+00,   3.25000000e+00,   3.30000000e+00,   3.35000000e+00,
1274        3.40000000e+00,   3.45000000e+00,   3.50000000e+00,   3.55000000e+00,
1275        3.60000000e+00,   3.65000000e+00,   3.70000000e+00,   3.75000000e+00,
1276        3.80000000e+00,   3.85000000e+00,   3.90000000e+00,   3.95000000e+00,
1277        4.00000000e+00,   4.05000000e+00,   4.10000000e+00,   4.15000000e+00,
1278        4.20000000e+00,   4.25000000e+00,   4.30000000e+00,   4.35000000e+00,
1279        4.40000000e+00,   4.45000000e+00,   4.50000000e+00,   4.55000000e+00,
1280        4.60000000e+00,   4.65000000e+00,   4.70000000e+00,   4.75000000e+00,
1281        4.80000000e+00,   4.85000000e+00,   4.90000000e+00,   4.95000000e+00,
1282        5.00000000e+00,   5.05000000e+00,   5.10000000e+00,   5.15000000e+00,
1283        5.20000000e+00,   5.25000000e+00,   5.30000000e+00,   5.35000000e+00,
1284        5.40000000e+00,   5.45000000e+00,   5.50000000e+00,   5.55000000e+00,
1285        5.60000000e+00,   5.65000000e+00,   5.70000000e+00,   5.75000000e+00,
1286        5.80000000e+00,   5.85000000e+00,   5.90000000e+00,   5.95000000e+00,
1287        6.00000000e+00,   6.05000000e+00,   6.10000000e+00,   6.15000000e+00,
1288        6.20000000e+00,   6.25000000e+00,   6.30000000e+00,   6.35000000e+00,
1289        6.40000000e+00,   6.45000000e+00,   6.50000000e+00,   6.55000000e+00,
1290        6.60000000e+00,   6.65000000e+00,   6.70000000e+00,   6.75000000e+00,
1291        6.80000000e+00,   6.85000000e+00,   6.90000000e+00,   6.95000000e+00,
1292        7.00000000e+00,   7.05000000e+00,   7.10000000e+00,   7.15000000e+00,
1293        7.20000000e+00,   7.25000000e+00,   7.30000000e+00,   7.35000000e+00,
1294        7.40000000e+00,   7.45000000e+00,   7.50000000e+00,   7.55000000e+00,
1295        7.60000000e+00,   7.65000000e+00,   7.70000000e+00,   7.75000000e+00,
1296        7.80000000e+00,   7.85000000e+00,   7.90000000e+00,   7.95000000e+00,
1297        8.00000000e+00,   8.05000000e+00,   8.10000000e+00,   8.15000000e+00,
1298        8.20000000e+00,   8.25000000e+00,   8.30000000e+00,   8.35000000e+00,
1299        8.40000000e+00,   8.45000000e+00,   8.50000000e+00,   8.55000000e+00,
1300        8.60000000e+00,   8.65000000e+00,   8.70000000e+00,   8.75000000e+00,
1301        8.80000000e+00,   8.85000000e+00,   8.90000000e+00,   8.95000000e+00,
1302        9.00000000e+00,   9.05000000e+00,   9.10000000e+00,   9.15000000e+00,
1303        9.20000000e+00,   9.25000000e+00,   9.30000000e+00,   9.35000000e+00,
1304        9.40000000e+00,   9.45000000e+00,   9.50000000e+00,   9.55000000e+00,
1305        9.60000000e+00,   9.65000000e+00,   9.70000000e+00,   9.75000000e+00,
1306        9.80000000e+00,   9.85000000e+00,   9.90000000e+00,   9.95000000e+00,
1307        1.00000000e+01,   1.00500000e+01,   1.01000000e+01,   1.01500000e+01,
1308        1.02000000e+01,   1.02500000e+01,   1.03000000e+01,   1.03500000e+01,
1309        1.04000000e+01,   1.04500000e+01,   1.05000000e+01,   1.05500000e+01,
1310        1.06000000e+01,   1.06500000e+01,   1.07000000e+01,   1.07500000e+01,
1311        1.08000000e+01,   1.08500000e+01,   1.09000000e+01,   1.09500000e+01,
1312        1.10000000e+01,   1.10500000e+01,   1.11000000e+01,   1.11500000e+01,
1313        1.12000000e+01,   1.12500000e+01,   1.13000000e+01,   1.13500000e+01,
1314        1.14000000e+01,   1.14500000e+01,   1.15000000e+01,   1.15500000e+01,
1315        1.16000000e+01,   1.16500000e+01,   1.17000000e+01,   1.17500000e+01,
1316        1.18000000e+01,   1.18500000e+01,   1.19000000e+01,   1.19500000e+01,
1317        1.20000000e+01,   1.20500000e+01,   1.21000000e+01,   1.21500000e+01,
1318        1.22000000e+01,   1.22500000e+01,   1.23000000e+01,   1.23500000e+01,
1319        1.24000000e+01,   1.24500000e+01,   1.25000000e+01,   1.25500000e+01,
1320        1.26000000e+01,   1.26500000e+01,   1.27000000e+01,   1.27500000e+01,
1321        1.28000000e+01,   1.28500000e+01,   1.29000000e+01,   1.29500000e+01,
1322        1.30000000e+01,   1.30500000e+01,   1.31000000e+01,   1.31500000e+01,
1323        1.32000000e+01,   1.32500000e+01,   1.33000000e+01,   1.33500000e+01,
1324        1.34000000e+01,   1.34500000e+01,   1.35000000e+01,   1.35500000e+01,
1325        1.36000000e+01,   1.36500000e+01,   1.37000000e+01,   1.37500000e+01,
1326        1.38000000e+01,   1.38500000e+01,   1.39000000e+01,   1.39500000e+01,
1327        1.40000000e+01,   1.40500000e+01,   1.41000000e+01,   1.41500000e+01,
1328        1.42000000e+01,   1.42500000e+01,   1.43000000e+01,   1.43500000e+01,
1329        1.44000000e+01,   1.44500000e+01,   1.45000000e+01,   1.45500000e+01,
1330        1.46000000e+01,   1.46500000e+01,   1.47000000e+01,   1.47500000e+01,
1331        1.48000000e+01,   1.48500000e+01,   1.49000000e+01,   1.49500000e+01,
1332        1.50000000e+01,   1.50500000e+01,   1.51000000e+01,   1.51500000e+01,
1333        1.52000000e+01,   1.52500000e+01,   1.53000000e+01,   1.53500000e+01,
1334        1.54000000e+01,   1.54500000e+01,   1.55000000e+01,   1.55500000e+01,
1335        1.56000000e+01,   1.56500000e+01,   1.57000000e+01,   1.57500000e+01,
1336        1.58000000e+01,   1.58500000e+01,   1.59000000e+01,   1.59500000e+01,
1337        1.60000000e+01,   1.60500000e+01,   1.61000000e+01,   1.61500000e+01,
1338        1.62000000e+01,   1.62500000e+01,   1.63000000e+01,   1.63500000e+01,
1339        1.64000000e+01,   1.64500000e+01,   1.65000000e+01,   1.65500000e+01,
1340        1.66000000e+01,   1.66500000e+01,   1.67000000e+01,   1.67500000e+01,
1341        1.68000000e+01,   1.68500000e+01,   1.69000000e+01,   1.69500000e+01,
1342        1.70000000e+01,   1.70500000e+01,   1.71000000e+01,   1.71500000e+01,
1343        1.72000000e+01,   1.72500000e+01,   1.73000000e+01,   1.73500000e+01,
1344        1.74000000e+01,   1.74500000e+01,   1.75000000e+01,   1.75500000e+01,
1345        1.76000000e+01,   1.76500000e+01,   1.77000000e+01,   1.77500000e+01,
1346        1.78000000e+01,   1.78500000e+01,   1.79000000e+01,   1.79500000e+01,
1347        1.80000000e+01,   1.80500000e+01,   1.81000000e+01,   1.81500000e+01,
1348        1.82000000e+01,   1.82500000e+01,   1.83000000e+01,   1.83500000e+01,
1349        1.84000000e+01,   1.84500000e+01,   1.85000000e+01,   1.85500000e+01,
1350        1.86000000e+01,   1.86500000e+01,   1.87000000e+01,   1.87500000e+01,
1351        1.88000000e+01,   1.88500000e+01,   1.89000000e+01,   1.89500000e+01,
1352        1.90000000e+01,   1.90500000e+01,   1.91000000e+01,   1.91500000e+01,
1353        1.92000000e+01,   1.92500000e+01,   1.93000000e+01,   1.93500000e+01,
1354        1.94000000e+01,   1.94500000e+01,   1.95000000e+01,   1.95500000e+01,
1355        1.96000000e+01,   1.96500000e+01,   1.97000000e+01,   1.97500000e+01,
1356        1.98000000e+01,   1.98500000e+01,   1.99000000e+01,   1.99500000e+01,
1357        2.00000000e+01,   2.00500000e+01,   2.01000000e+01,   2.01500000e+01,
1358        2.02000000e+01,   2.02500000e+01,   2.03000000e+01,   2.03500000e+01,
1359        2.04000000e+01,   2.04500000e+01,   2.05000000e+01,   2.05500000e+01,
1360        2.06000000e+01,   2.06500000e+01,   2.07000000e+01,   2.07500000e+01,
1361        2.08000000e+01,   2.08500000e+01,   2.09000000e+01,   2.09500000e+01,
1362        2.10000000e+01,   2.10500000e+01,   2.11000000e+01,   2.11500000e+01,
1363        2.12000000e+01,   2.12500000e+01,   2.13000000e+01,   2.13500000e+01,
1364        2.14000000e+01,   2.14500000e+01,   2.15000000e+01,   2.15500000e+01,
1365        2.16000000e+01,   2.16500000e+01,   2.17000000e+01,   2.17500000e+01,
1366        2.18000000e+01,   2.18500000e+01,   2.19000000e+01,   2.19500000e+01,
1367        2.20000000e+01,   2.20500000e+01,   2.21000000e+01,   2.21500000e+01,
1368        2.22000000e+01,   2.22500000e+01,   2.23000000e+01,   2.23500000e+01,
1369        2.24000000e+01,   2.24500000e+01,   2.25000000e+01])
1370
1371        #print 'Diff', time[1:] - time[:-1]
1372
1373        #Setup mesh used to represent fitted function
1374        a = [0.0, 0.0]
1375        b = [0.0, 2.0]
1376        c = [2.0, 0.0]
1377        d = [0.0, 4.0]
1378        e = [2.0, 2.0]
1379        f = [4.0, 0.0]
1380
1381        points = [a, b, c, d, e, f]
1382        #bac, bce, ecf, dbe
1383        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1384
1385
1386        #New datapoints where interpolated values are sought
1387        interpolation_points = [[ 0.0, 0.0],
1388                                [ 0.5, 0.5],
1389                                [ 0.7, 0.7],
1390                                [ 1.0, 0.5],
1391                                [ 2.0, 0.4],
1392                                [ 545354534, 4354354353]] # outside the mesh
1393
1394        #One quantity
1395        Q = zeros( (len(time),6), Float )
1396
1397        #Linear in time and space
1398        for i, t in enumerate(time):
1399            Q[i, :] = t*linear_function(points)
1400
1401        #Check interpolation of one quantity using interpolaton points)
1402        try:
1403            I = Interpolation_function(time, Q,
1404                                       vertex_coordinates = points,
1405                                       triangles = triangles, 
1406                                       interpolation_points = interpolation_points,
1407                                       verbose = False)
1408        except:
1409            pass
1410        else:
1411            raise 'Should raise exception due to time being non-monotoneous'           
1412     
1413
1414    def test_points_outside_the_polygon(self):
1415        a = [-1.0, 0.0]
1416        b = [3.0, 4.0]
1417        c = [4.0, 1.0]
1418        d = [-3.0, 2.0] #3
1419        e = [-1.0, -2.0]
1420        f = [1.0, -2.0] #5
1421
1422        vertices = [a, b, c, d,e,f]
1423        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
1424
1425        point_coords = [[-2.0, 2.0],
1426                        [-1.0, 1.0],
1427                        [9999.0, 9999.0], # point Outside poly
1428                        [-9999.0, 1.0], # point Outside poly
1429                        [2.0, 1.0],
1430                        [0.0, 0.0],
1431                        [1.0, 0.0],
1432                        [0.0, -1.0],
1433                        [-0.2, -0.5],
1434                        [-0.9, -1.5],
1435                        [0.5, -1.9],
1436                        [999999, 9999999]] # point Outside poly
1437        geo_data = Geospatial_data(data_points = point_coords)
1438
1439        interp = Interpolate(vertices, triangles)
1440        f = array([linear_function(vertices),2*linear_function(vertices) ])
1441        f = transpose(f)
1442        #print "f",f
1443        z = interp.interpolate(f, geo_data)
1444        #z = interp.interpolate(f, point_coords)
1445        answer = [linear_function(point_coords),
1446                  2*linear_function(point_coords) ]
1447        answer = transpose(answer)
1448        answer[2,:] = [NAN, NAN]
1449        answer[3,:] = [NAN, NAN]
1450        answer[11,:] = [NAN, NAN]
1451        #print "z",z
1452        #print "answer _ fixed",answer
1453        assert allclose(z[0:1], answer[0:1])
1454        assert allclose(z[4:10], answer[4:10])
1455        for i in [2,3,11]:
1456            self.failUnless( z[i,1] == answer[11,1], 'Fail!')
1457            self.failUnless( z[i,0] == answer[11,0], 'Fail!')
1458
1459    def test_interpolate_sww2csv(self):
1460
1461        def elevation_function(x, y):
1462            return -x
1463       
1464        # create mesh
1465        mesh_file = tempfile.mktemp(".tsh")   
1466        points = [[0.0,0.0],[6.0,0.0],[6.0,6.0],[0.0,6.0]]
1467        m = Mesh()
1468        m.add_vertices(points)
1469        m.auto_segment()
1470        m.generate_mesh(verbose=False)
1471        m.export_mesh_file(mesh_file)
1472       
1473        #Create shallow water domain
1474        domain = Domain(mesh_file)
1475        os.remove(mesh_file)
1476       
1477        domain.default_order=2
1478        domain.beta_h = 0
1479
1480        #Set some field values
1481        domain.set_quantity('elevation', elevation_function)
1482        domain.set_quantity('friction', 0.03)
1483        domain.set_quantity('xmomentum', 3.0)
1484        domain.set_quantity('ymomentum', 4.0)
1485
1486        ######################
1487        # Boundary conditions
1488        B = Transmissive_boundary(domain)
1489        domain.set_boundary( {'exterior': B})
1490
1491        # This call mangles the stage values.
1492        domain.distribute_to_vertices_and_edges()
1493        domain.set_quantity('stage', 1.0)
1494
1495        #sww_file = tempfile.mktemp("")
1496        domain.filename = 'datatest' + str(time.time())
1497        #domain.filename = sww_file
1498        #print "domain.filename",domain.filename
1499        domain.format = 'sww'
1500        domain.smooth = True
1501        domain.reduction = mean
1502
1503        sww = get_dataobject(domain)
1504        sww.store_connectivity()
1505        sww.store_timestep(['stage', 'xmomentum', 'ymomentum'])
1506        domain.set_quantity('stage', 10.0) # This is automatically limmited
1507        # so it will not be less than the elevation
1508        domain.time = 2.
1509        sww.store_timestep(['stage', 'xmomentum', 'ymomentum'])
1510
1511        # test the function
1512        points = [[5.0,1.],[0.5,2.]]
1513        depth_file = tempfile.mktemp(".csv") 
1514        velocity_x_file = tempfile.mktemp(".csv") 
1515        velocity_y_file = tempfile.mktemp(".csv") 
1516        interpolate_sww2csv(sww.filename, points, depth_file,
1517                            velocity_x_file,
1518                            velocity_y_file,
1519                            verbose=False)
1520
1521        depth_answers_array = [[0.0, 6.0, 1.5], [2.0, 15., 10.5]] 
1522        velocity_x_answers_array = [[0.0, 3./6.0, 3./1.5],
1523                                    [2.0, 3./15., 3/10.5]]
1524        velocity_y_answers_array = [[0.0, 4./6.0, 4./1.5],
1525                                    [2.0, 4./15., 4./10.5]]
1526        depth_file_handle = file(depth_file)
1527        depth_reader = csv.reader(depth_file_handle)
1528        depth_reader.next()
1529        velocity_x_file_handle = file(velocity_x_file)
1530        velocity_x_reader = csv.reader(velocity_x_file_handle)
1531        velocity_x_reader.next()
1532        for depths, velocitys, depth_answers, velocity_answers in map(None,
1533                                              depth_reader,
1534                                              velocity_x_reader,
1535                                              depth_answers_array,
1536                                              velocity_x_answers_array):
1537            for i in range(len(depths)):
1538                #print "depths",depths[i]
1539                #print "depth_answers",depth_answers[i]
1540                #print "velocitys",velocitys[i]
1541                #print "velocity_answers_array", velocity_answers[i]
1542                msg = 'Interpolation failed'
1543                assert allclose(float(depths[i]), depth_answers[i]), msg
1544                assert allclose(float(velocitys[i]), velocity_answers[i]), msg
1545
1546        velocity_y_file_handle = file(velocity_y_file)
1547        velocity_y_reader = csv.reader(velocity_y_file_handle)
1548        velocity_y_reader.next()
1549        for velocitys, velocity_answers in map(None,
1550                                              velocity_y_reader,
1551                                              velocity_y_answers_array):
1552            for i in range(len(depths)):
1553                #print "depths",depths[i]
1554                #print "depth_answers",depth_answers[i]
1555                #print "velocitys",velocitys[i]
1556                #print "velocity_answers_array", velocity_answers[i]
1557                msg = 'Interpolation failed'
1558                assert allclose(float(depths[i]), depth_answers[i]), msg
1559                assert allclose(float(velocitys[i]), velocity_answers[i]), msg
1560               
1561        # clean up
1562        depth_file_handle.close()
1563        velocity_y_file_handle.close()
1564        velocity_x_file_handle.close()
1565        #print "sww.filename",sww.filename
1566        os.remove(sww.filename)
1567        os.remove(depth_file)
1568        os.remove(velocity_x_file)
1569        os.remove(velocity_y_file)
1570#-------------------------------------------------------------
1571if __name__ == "__main__":
1572
1573    suite = unittest.makeSuite(Test_Interpolate,'test')
1574    #suite = unittest.makeSuite(Test_Interpolate,'test_interpolation_function_outside_point')
1575    runner = unittest.TextTestRunner(verbosity=1)
1576    runner.run(suite)
1577
1578
1579
1580
1581
Note: See TracBrowser for help on using the repository browser.