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

Last change on this file since 8709 was 8709, checked in by wilsonp, 11 years ago

Changes to fitsmooth.c to use PyCObject when Python version is detected as lower than 2.7.3. Also, small fixes to caching of quad tree.

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