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

Last change on this file since 5442 was 5442, checked in by ole, 15 years ago

Retired h-limiter and beta_h as per ticket:194.
All unit tests and validation tests pass.

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