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

Last change on this file since 7673 was 7673, checked in by hudson, 13 years ago

Ticket 113 - added an output_centroids parameter to allow centroid data to be written to gauges, rather than the data at the sampled point.

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