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

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

Various refactorings, all unit tests pass.
Domain renamed to generic domain.

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