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

Last change on this file since 5775 was 5775, checked in by ole, 16 years ago

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