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

Last change on this file since 6183 was 6183, checked in by ole, 16 years ago

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