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

Last change on this file since 3689 was 3689, checked in by ole, 17 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()
138
139        I = Interpolate(vertex_coordinates, triangles)
140        result = I.interpolate(vertex_values, interpolation_points)
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()
158
159        I = Interpolate(vertex_coordinates, triangles)
160        result = I.interpolate(vertex_values, interpolation_points)
162
163
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(),
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(),
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(),
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(),
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(),
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(),
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]):
348
349
350        #results = interp._build_interpolation_matrix_A(data).todense()
351
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)
376
377        #print "z",z
380
381
382        z = interp.interpolate(f, point_coords, start_blocking_len = 2)
384
385        #print "z",z
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)
412
413        #print "z",z
416
417        z = interp.interpolate(f, point_coords, start_blocking_len = 2)
419
420        #print "z",z
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)
441
442        #print "z",z
445
446
447        z = interp.interpolate(f, point_coords, start_blocking_len = 2)
449
450        #print "z",z
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)
483        #print "z",z
486
487        z = interp.interpolate(f, point_coords, start_blocking_len = 2)
489
490        #print "z",z
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
544        #print "***********"
545
547
548
549        z = interp.interpolate(f, point_coords, start_blocking_len = 2)
550
551        #print "z",z
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
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
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)
635                  2*linear_function(point_coords) ]
637        #print "z",z
640
641        z = interp.interpolate(f, point_coords, start_blocking_len = 2)
642
643        #print "z",z
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)
682                      2*linear_function(point_coords) ]
684            #print "z",z
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)
698                      2*linear_function(point_coords) ,
699                      2*linear_function(point_coords)-100 ])
700            z = transpose(z)
701            #print "z",z
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)
744                      absolute = True)),
745                      2*linear_function(point_coords.get_data_points( \
746                      absolute = True)) ]
748            #print "z",z
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)
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
771
772        z = interp.interpolate(f, point_coords, start_blocking_len = 2)
773
774        #print "z",z
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)
813                      absolute = True)),
814                  2*linear_function(point_coords.get_data_points( \
815                      absolute = True)) ]
817        #print "z",z
820
821        z = interp.interpolate(f, point_coords, start_blocking_len = 2)
822
823        #print "z",z
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)
859                  2*linear_function(point_coords) ]
861        #print "z",z
864        assert allclose(interp._A_can_be_reused, True)
865
866        z = interp.interpolate(f)
868
869        # This causes blocking to occur.
870        z = interp.interpolate(f, start_blocking_len=10)
872        assert allclose(interp._A_can_be_reused, False)
873
874        #A is recalculated
875        z = interp.interpolate(f)
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)
921                  2*linear_function(point_coords) ]
923
925        assert allclose(interp._A_can_be_reused, True)
926
927
928        z = interp.interpolate(f)    # None
930        z = interp.interpolate(f, point_coords) # Repeated (not really a test)
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
1051
1052        t = time[0]
1053        for j in range(50): #t in [1, 6]
1054            for id in range(len(interpolation_points)):
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
1113
1114        t = time[0]
1115        for j in range(50): #t in [1, 6]
1116            for id in range(len(interpolation_points)):
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)
1165                  2*linear_function(interpolation_points) ]
1167        #print "z",z
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
1230
1231        t = time[0]
1232        for j in range(50): #t in [1, 6]
1233            for id in range(len(interpolation_points)-1):
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)
1446                  2*linear_function(point_coords) ]
1451        #print "z",z
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()
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)
1529        velocity_x_file_handle = file(velocity_x_file)
1537            for i in range(len(depths)):
1538                #print "depths",depths[i]
1540                #print "velocitys",velocitys[i]
1542                msg = 'Interpolation failed'
1545
1546        velocity_y_file_handle = file(velocity_y_file)
1549        for velocitys, velocity_answers in map(None,
1552            for i in range(len(depths)):
1553                #print "depths",depths[i]
1555                #print "velocitys",velocitys[i]
1557                msg = 'Interpolation failed'
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.