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

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

Ticket 344 completed and unit tests pass.
Stripped out a heap of tabs and replaced them with 4 spaces.
Removed some redundant imports and deprecated functions.

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