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

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

Refactored stofrage of quantities to use new concept of static and dynamic quantities.
This is in preparation for flexibly allowing quantities such as elevation or friction
to be time dependent.

All tests and the okushiri validation pass.

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