source: trunk/anuga_core/source/anuga/fit_interpolate/test_interpolate.py @ 7876

Last change on this file since 7876 was 7796, checked in by hudson, 14 years ago

ANUGA core modified to fix errors with new API modules not being found.

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