source: anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py @ 5349

Last change on this file since 5349 was 5349, checked in by duncan, 16 years ago

Minor fixes/ comments

File size: 78.5 KB
Line 
1#!/usr/bin/env python
2
3import unittest
4from math import sqrt, pi
5import tempfile
6
7from quantity import *
8from anuga.config import epsilon
9from Numeric import allclose, array, ones, Float
10
11from anuga.fit_interpolate.fit import fit_to_mesh
12#from anuga.pyvolution.least_squares import fit_to_mesh         
13from domain import Domain
14from anuga.geospatial_data.geospatial_data import Geospatial_data
15from anuga.coordinate_transforms.geo_reference import Geo_reference
16
17#Aux for fit_interpolate.fit example
18def linear_function(point):
19    point = array(point)
20    return point[:,0]+point[:,1]
21
22
23class Test_Quantity(unittest.TestCase):
24    def setUp(self):
25        from domain import Domain
26
27        a = [0.0, 0.0]
28        b = [0.0, 2.0]
29        c = [2.0, 0.0]
30        d = [0.0, 4.0]
31        e = [2.0, 2.0]
32        f = [4.0, 0.0]
33
34        points = [a, b, c, d, e, f]
35
36        #bac, bce, ecf, dbe
37        elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
38
39        self.mesh1 = Domain(points[:3], [elements[0]])
40        self.mesh1.check_integrity()
41
42        self.mesh4 = Domain(points, elements)
43        self.mesh4.check_integrity()
44
45        # UTM round Onslow
46        a = [240000, 7620000]
47        b = [240000, 7680000]
48        c = [300000, 7620000]
49
50        points = [a, b, c]
51        elements = [[0,2,1]]
52       
53        self.mesh_onslow = Domain(points, elements)
54        self.mesh_onslow.check_integrity()
55       
56    def tearDown(self):
57        pass
58        #print "  Tearing down"
59
60
61    def test_creation(self):
62
63        quantity = Quantity(self.mesh1, [[1,2,3]])
64        assert allclose(quantity.vertex_values, [[1.,2.,3.]])
65
66        try:
67            quantity = Quantity()
68        except:
69            pass
70        else:
71            raise 'Should have raised empty quantity exception'
72
73
74        try:
75            quantity = Quantity([1,2,3])
76        except AssertionError:
77            pass
78        except:
79            raise 'Should have raised "mising mesh object" error'
80
81
82    def test_creation_zeros(self):
83
84        quantity = Quantity(self.mesh1)
85        assert allclose(quantity.vertex_values, [[0.,0.,0.]])
86
87
88        quantity = Quantity(self.mesh4)
89        assert allclose(quantity.vertex_values, [[0.,0.,0.], [0.,0.,0.],
90                                                 [0.,0.,0.], [0.,0.,0.]])
91
92
93    def test_interpolation(self):
94        quantity = Quantity(self.mesh1, [[1,2,3]])
95        assert allclose(quantity.centroid_values, [2.0]) #Centroid
96
97        assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5]])
98
99
100    def test_interpolation2(self):
101        quantity = Quantity(self.mesh4,
102                            [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
103        assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
104
105
106        quantity.extrapolate_second_order()
107
108        #print quantity.vertex_values
109        assert allclose(quantity.vertex_values, [[3.5, -1.0, 3.5],
110                                                 [3.+2./3, 6.+2./3, 4.+2./3],
111                                                 [4.6, 3.4, 1.],
112                                                 [-5.0, 1.0, 4.0]])
113
114        #print quantity.edge_values
115        assert allclose(quantity.edge_values, [[1.25, 3.5, 1.25],
116                                               [5. + 2/3.0, 4.0 + 1.0/6, 5.0 + 1.0/6],
117                                               [2.2, 2.8, 4.0],
118                                               [2.5, -0.5, -2.0]])
119
120
121        #assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
122        #                                       [5., 5., 5.],
123        #                                       [4.5, 4.5, 0.],
124        #                                       [3.0, -1.5, -1.5]])
125
126    def test_get_extrema_1(self):
127        quantity = Quantity(self.mesh4,
128                                      [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
129        assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroids
130
131        v = quantity.get_maximum_value()
132        assert v == 5
133
134        v = quantity.get_minimum_value()
135        assert v == 0       
136
137        i = quantity.get_maximum_index()
138        assert i == 1
139
140        i = quantity.get_minimum_index()
141        assert i == 3       
142       
143        x,y = quantity.get_maximum_location()
144        xref, yref = 4.0/3, 4.0/3
145        assert x == xref
146        assert y == yref
147
148        v = quantity.get_values(interpolation_points = [[x,y]])
149        assert allclose(v, 5)
150
151
152        x,y = quantity.get_minimum_location()
153        v = quantity.get_values(interpolation_points = [[x,y]])
154        assert allclose(v, 0)
155
156
157    def test_get_maximum_2(self):
158
159        a = [0.0, 0.0]
160        b = [0.0, 2.0]
161        c = [2.0,0.0]
162        d = [0.0, 4.0]
163        e = [2.0, 2.0]
164        f = [4.0,0.0]
165
166        points = [a, b, c, d, e, f]
167        #bac, bce, ecf, dbe
168        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
169
170        domain = Domain(points, vertices)
171
172        quantity = Quantity(domain)
173        quantity.set_values(lambda x, y: x+2*y) #2 4 4 6
174       
175        v = quantity.get_maximum_value()
176        assert v == 6
177
178        v = quantity.get_minimum_value()
179        assert v == 2       
180
181        i = quantity.get_maximum_index()
182        assert i == 3
183
184        i = quantity.get_minimum_index()
185        assert i == 0       
186       
187        x,y = quantity.get_maximum_location()
188        xref, yref = 2.0/3, 8.0/3
189        assert x == xref
190        assert y == yref
191
192        v = quantity.get_values(interpolation_points = [[x,y]])
193        assert allclose(v, 6)
194
195        x,y = quantity.get_minimum_location()       
196        v = quantity.get_values(interpolation_points = [[x,y]])
197        assert allclose(v, 2)
198
199        #Multiple locations for maximum -
200        #Test that the algorithm picks the first occurrence       
201        v = quantity.get_maximum_value(indices=[0,1,2])
202        assert allclose(v, 4)
203
204        i = quantity.get_maximum_index(indices=[0,1,2])
205        assert i == 1
206       
207        x,y = quantity.get_maximum_location(indices=[0,1,2])
208        xref, yref = 4.0/3, 4.0/3
209        assert x == xref
210        assert y == yref
211
212        v = quantity.get_values(interpolation_points = [[x,y]])
213        assert allclose(v, 4)       
214
215        # More test of indices......
216        v = quantity.get_maximum_value(indices=[2,3])
217        assert allclose(v, 6)
218
219        i = quantity.get_maximum_index(indices=[2,3])
220        assert i == 3
221       
222        x,y = quantity.get_maximum_location(indices=[2,3])
223        xref, yref = 2.0/3, 8.0/3
224        assert x == xref
225        assert y == yref
226
227        v = quantity.get_values(interpolation_points = [[x,y]])
228        assert allclose(v, 6)       
229
230       
231
232    def test_boundary_allocation(self):
233        quantity = Quantity(self.mesh4,
234                            [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
235
236        assert quantity.boundary_values.shape[0] == len(self.mesh4.boundary)
237
238
239    def test_set_values(self):
240        quantity = Quantity(self.mesh4)
241
242
243        quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
244                            location = 'vertices')
245        assert allclose(quantity.vertex_values,
246                        [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
247        assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
248        assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
249                                               [5., 5., 5.],
250                                               [4.5, 4.5, 0.],
251                                               [3.0, -1.5, -1.5]])
252
253
254        #Test default
255        quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
256        assert allclose(quantity.vertex_values,
257                        [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
258        assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
259        assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
260                                               [5., 5., 5.],
261                                               [4.5, 4.5, 0.],
262                                               [3.0, -1.5, -1.5]])
263
264        #Test centroids
265        quantity.set_values([1,2,3,4], location = 'centroids')
266        assert allclose(quantity.centroid_values, [1., 2., 3., 4.]) #Centroid
267
268        #Test edges
269        quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
270                            location = 'edges')
271        assert allclose(quantity.edge_values,
272                        [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
273
274        #Test exceptions
275        try:
276            quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
277                                location = 'bas kamel tuba')
278        except:
279            pass
280
281
282        try:
283            quantity.set_values([[1,2,3], [0,0,9]])
284        except AssertionError:
285            pass
286        except:
287            raise 'should have raised Assertionerror'
288
289
290
291    def test_set_values_const(self):
292        quantity = Quantity(self.mesh4)
293
294        quantity.set_values(1.0, location = 'vertices')
295        assert allclose(quantity.vertex_values,
296                        [[1,1,1], [1,1,1], [1,1,1], [1, 1, 1]])
297
298        assert allclose(quantity.centroid_values, [1, 1, 1, 1]) #Centroid
299        assert allclose(quantity.edge_values, [[1, 1, 1],
300                                               [1, 1, 1],
301                                               [1, 1, 1],
302                                               [1, 1, 1]])
303
304
305        quantity.set_values(2.0, location = 'centroids')
306        assert allclose(quantity.centroid_values, [2, 2, 2, 2])
307
308        quantity.set_values(3.0, location = 'edges')
309        assert allclose(quantity.edge_values, [[3, 3, 3],
310                                               [3, 3, 3],
311                                               [3, 3, 3],
312                                               [3, 3, 3]])
313
314
315    def test_set_values_func(self):
316        quantity = Quantity(self.mesh4)
317
318        def f(x, y):
319            return x+y
320
321        quantity.set_values(f, location = 'vertices')
322        #print "quantity.vertex_values",quantity.vertex_values
323        assert allclose(quantity.vertex_values,
324                        [[2,0,2], [2,2,4], [4,2,4], [4,2,4]])
325        assert allclose(quantity.centroid_values,
326                        [4.0/3, 8.0/3, 10.0/3, 10.0/3])
327        assert allclose(quantity.edge_values,
328                        [[1,2,1], [3,3,2], [3,4,3], [3,4,3]])
329
330
331        quantity.set_values(f, location = 'centroids')
332        assert allclose(quantity.centroid_values,
333                        [4.0/3, 8.0/3, 10.0/3, 10.0/3])
334
335
336    def test_integral(self):
337        quantity = Quantity(self.mesh4)
338
339        #Try constants first
340        const = 5
341        quantity.set_values(const, location = 'vertices')
342        #print 'Q', quantity.get_integral()
343
344        assert allclose(quantity.get_integral(), self.mesh4.get_area() * const)
345
346        #Try with a linear function
347        def f(x, y):
348            return x+y
349
350        quantity.set_values(f, location = 'vertices')
351
352
353        ref_integral = (4.0/3 + 8.0/3 + 10.0/3 + 10.0/3) * 2
354
355        assert allclose (quantity.get_integral(), ref_integral)
356
357
358
359    def test_set_vertex_values(self):
360        quantity = Quantity(self.mesh4)
361        quantity.set_vertex_values([0,1,2,3,4,5])
362
363        assert allclose(quantity.vertex_values,
364                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
365        assert allclose(quantity.centroid_values,
366                        [1., 7./3, 11./3, 8./3]) #Centroid
367        assert allclose(quantity.edge_values, [[1., 1.5, 0.5],
368                                               [3., 2.5, 1.5],
369                                               [3.5, 4.5, 3.],
370                                               [2.5, 3.5, 2]])
371
372
373    def test_set_vertex_values_subset(self):
374        quantity = Quantity(self.mesh4)
375        quantity.set_vertex_values([0,1,2,3,4,5])
376        quantity.set_vertex_values([0,20,30,50], indices = [0,2,3,5])
377
378        assert allclose(quantity.vertex_values,
379                        [[1,0,20], [1,20,4], [4,20,50], [30,1,4]])
380
381
382    def test_set_vertex_values_using_general_interface(self):
383        quantity = Quantity(self.mesh4)
384
385
386        quantity.set_values([0,1,2,3,4,5])
387
388
389        assert allclose(quantity.vertex_values,
390                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
391
392        #Centroid
393        assert allclose(quantity.centroid_values, [1., 7./3, 11./3, 8./3])
394
395        assert allclose(quantity.edge_values, [[1., 1.5, 0.5],
396                                               [3., 2.5, 1.5],
397                                               [3.5, 4.5, 3.],
398                                               [2.5, 3.5, 2]])
399
400
401
402    def test_set_vertex_values_using_general_interface_with_subset(self):
403        """test_set_vertex_values_using_general_interface_with_subset(self):
404        Test that indices and polygon works
405        """
406       
407        quantity = Quantity(self.mesh4)
408
409
410        quantity.set_values([0,2,3,5], indices=[0,2,3,5])
411        assert allclose(quantity.vertex_values,
412                        [[0,0,2], [0,2,0], [0,2,5], [3,0,0]])
413
414
415        # Constant
416        quantity.set_values(0.0)
417        quantity.set_values(3.14, indices=[0,2], location='vertices')
418
419        # Indices refer to triangle numbers here - not vertices (why?)
420        assert allclose(quantity.vertex_values,
421                        [[3.14,3.14,3.14], [0,0,0],
422                         [3.14,3.14,3.14], [0,0,0]])       
423       
424
425
426        # Now try with polygon (pick points where y>2)
427        polygon = [[0,2.1], [4,2.1], [4,7], [0,7]]
428        quantity.set_values(0.0)
429        quantity.set_values(3.14, polygon=polygon)
430       
431        assert allclose(quantity.vertex_values,
432                        [[0,0,0], [0,0,0], [0,0,0],
433                         [3.14,3.14,3.14]])               
434
435
436        # Another polygon (pick triangle 1 and 2 (rightmost triangles)
437        polygon = [[2.1, 0.0], [3.5,0.1], [2,2.2], [0.2,2]]
438        quantity.set_values(0.0)
439        quantity.set_values(3.14, polygon=polygon)
440
441        assert allclose(quantity.vertex_values,
442                        [[0,0,0],
443                         [3.14,3.14,3.14],
444                         [3.14,3.14,3.14],                         
445                         [0,0,0]])               
446
447
448
449        # Test input checking
450        try:
451            quantity.set_values(3.14, polygon=polygon, indices = [0,2])
452        except:
453            pass
454        else:
455            msg = 'Should have caught this'
456            raise msg
457
458
459
460
461
462    def test_set_vertex_values_using_general_interface_subset_and_geo(self):
463        """test_set_vertex_values_using_general_interface_with_subset(self):
464        Test that indices and polygon works using georeferencing
465        """
466       
467        quantity = Quantity(self.mesh4)
468        G = Geo_reference(56, 10, 100)
469        quantity.domain.geo_reference = G
470
471        #print quantity.domain.get_nodes(absolute=True)
472
473
474        # Constant
475        quantity.set_values(0.0)
476        quantity.set_values(3.14, indices=[0,2], location='vertices')
477
478        # Indices refer to triangle numbers here - not vertices (why?)
479        assert allclose(quantity.vertex_values,
480                        [[3.14,3.14,3.14], [0,0,0],
481                         [3.14,3.14,3.14], [0,0,0]])       
482       
483
484
485        # Now try with polygon (pick points where y>2)
486        polygon = array([[0,2.1], [4,2.1], [4,7], [0,7]])
487        polygon += [G.xllcorner, G.yllcorner]
488       
489        quantity.set_values(0.0)
490        quantity.set_values(3.14, polygon=polygon)
491       
492        assert allclose(quantity.vertex_values,
493                        [[0,0,0], [0,0,0], [0,0,0],
494                         [3.14,3.14,3.14]])               
495
496
497        # Another polygon (pick triangle 1 and 2 (rightmost triangles)
498        polygon = array([[2.1, 0.0], [3.5,0.1], [2,2.2], [0.2,2]])
499        polygon += [G.xllcorner, G.yllcorner]
500       
501        quantity.set_values(0.0)
502        quantity.set_values(3.14, polygon=polygon)
503
504        assert allclose(quantity.vertex_values,
505                        [[0,0,0],
506                         [3.14,3.14,3.14],
507                         [3.14,3.14,3.14],                         
508                         [0,0,0]])               
509
510
511
512    def test_set_values_using_fit(self):
513
514
515        quantity = Quantity(self.mesh4)
516
517        #Get (enough) datapoints
518        data_points = [[ 0.66666667, 0.66666667],
519                       [ 1.33333333, 1.33333333],
520                       [ 2.66666667, 0.66666667],
521                       [ 0.66666667, 2.66666667],
522                       [ 0.0, 1.0],
523                       [ 0.0, 3.0],
524                       [ 1.0, 0.0],
525                       [ 1.0, 1.0],
526                       [ 1.0, 2.0],
527                       [ 1.0, 3.0],
528                       [ 2.0, 1.0],
529                       [ 3.0, 0.0],
530                       [ 3.0, 1.0]]
531
532        z = linear_function(data_points)
533
534        #Use built-in fit_interpolate.fit
535        quantity.set_values( Geospatial_data(data_points, z), alpha = 0 )
536        #quantity.set_values(points = data_points, values = z, alpha = 0)
537
538
539        answer = linear_function(quantity.domain.get_vertex_coordinates())
540        #print quantity.vertex_values, answer
541        assert allclose(quantity.vertex_values.flat, answer)
542
543
544        #Now try by setting the same values directly
545        vertex_attributes = fit_to_mesh(data_points,
546                                        quantity.domain.get_nodes(),
547                                        quantity.domain.triangles, #FIXME
548                                        point_attributes=z,
549                                        alpha = 0,
550                                        verbose=False)
551
552        #print vertex_attributes
553        quantity.set_values(vertex_attributes)
554        assert allclose(quantity.vertex_values.flat, answer)
555
556
557
558
559
560    def test_test_set_values_using_fit_w_geo(self):
561
562
563        #Mesh
564        vertex_coordinates = [[0.76, 0.76],
565                              [0.76, 5.76],
566                              [5.76, 0.76]]
567        triangles = [[0,2,1]]
568
569        mesh_georef = Geo_reference(56,-0.76,-0.76)
570        mesh1 = Domain(vertex_coordinates, triangles,
571                       geo_reference = mesh_georef)
572        mesh1.check_integrity()
573
574        #Quantity
575        quantity = Quantity(mesh1)
576
577        #Data
578        data_points = [[ 201.0, 401.0],
579                       [ 201.0, 403.0],
580                       [ 203.0, 401.0]]
581
582        z = [2, 4, 4]
583
584        data_georef = Geo_reference(56,-200,-400)
585
586
587        #Reference
588        ref = fit_to_mesh(data_points, vertex_coordinates, triangles,
589                          point_attributes=z,
590                          data_origin = data_georef.get_origin(),
591                          mesh_origin = mesh_georef.get_origin(),
592                          alpha = 0)
593
594        assert allclose( ref, [0,5,5] )
595
596
597        #Test set_values
598
599        quantity.set_values( Geospatial_data(data_points, z, data_georef), alpha = 0 )
600
601        #quantity.set_values(points = data_points,
602        #                    values = z,
603        #                    data_georef = data_georef,
604        #                    alpha = 0)
605
606
607        #quantity.set_values(points = data_points,
608        #                    values = z,
609        #                    data_georef = data_georef,
610        #                    alpha = 0)
611        assert allclose(quantity.vertex_values.flat, ref)
612
613
614
615        #Test set_values using geospatial data object
616        quantity.vertex_values[:] = 0.0
617
618        geo = Geospatial_data(data_points, z, data_georef)
619
620
621        quantity.set_values(geospatial_data = geo, alpha = 0)
622        assert allclose(quantity.vertex_values.flat, ref)
623
624
625
626    def test_set_values_from_file1(self):
627        quantity = Quantity(self.mesh4)
628
629        #Get (enough) datapoints
630        data_points = [[ 0.66666667, 0.66666667],
631                       [ 1.33333333, 1.33333333],
632                       [ 2.66666667, 0.66666667],
633                       [ 0.66666667, 2.66666667],
634                       [ 0.0, 1.0],
635                       [ 0.0, 3.0],
636                       [ 1.0, 0.0],
637                       [ 1.0, 1.0],
638                       [ 1.0, 2.0],
639                       [ 1.0, 3.0],
640                       [ 2.0, 1.0],
641                       [ 3.0, 0.0],
642                       [ 3.0, 1.0]]
643
644        data_geo_spatial = Geospatial_data(data_points,
645                         geo_reference = Geo_reference(56, 0, 0))
646        data_points_absolute = data_geo_spatial.get_data_points(absolute=True)
647        attributes = linear_function(data_points_absolute)
648        att = 'spam_and_eggs'
649       
650        #Create .txt file
651        ptsfile = tempfile.mktemp(".txt")
652        file = open(ptsfile,"w")
653        file.write(" x,y," + att + " \n")
654        for data_point, attribute in map(None, data_points_absolute
655                                         ,attributes):
656            row = str(data_point[0]) + ',' + str(data_point[1]) \
657                  + ',' + str(attribute)
658            file.write(row + "\n")
659        file.close()
660
661
662        #Check that values can be set from file
663        quantity.set_values(filename = ptsfile,
664                            attribute_name = att, alpha = 0)
665        answer = linear_function(quantity.domain.get_vertex_coordinates())
666
667        #print quantity.vertex_values.flat
668        #print answer
669
670
671        assert allclose(quantity.vertex_values.flat, answer)
672
673
674        #Check that values can be set from file using default attribute
675        quantity.set_values(filename = ptsfile, alpha = 0)
676        assert allclose(quantity.vertex_values.flat, answer)
677
678        #Cleanup
679        import os
680        os.remove(ptsfile)
681
682    def Cache_cache_test_set_values_from_file(self):
683        quantity = Quantity(self.mesh4)
684
685        #Get (enough) datapoints
686        data_points = [[ 0.66666667, 0.66666667],
687                       [ 1.33333333, 1.33333333],
688                       [ 2.66666667, 0.66666667],
689                       [ 0.66666667, 2.66666667],
690                       [ 0.0, 1.0],
691                       [ 0.0, 3.0],
692                       [ 1.0, 0.0],
693                       [ 1.0, 1.0],
694                       [ 1.0, 2.0],
695                       [ 1.0, 3.0],
696                       [ 2.0, 1.0],
697                       [ 3.0, 0.0],
698                       [ 3.0, 1.0]]
699
700        data_geo_spatial = Geospatial_data(data_points,
701                         geo_reference = Geo_reference(56, 0, 0))
702        data_points_absolute = data_geo_spatial.get_data_points(absolute=True)
703        attributes = linear_function(data_points_absolute)
704        att = 'spam_and_eggs'
705       
706        #Create .txt file
707        ptsfile = tempfile.mktemp(".txt")
708        file = open(ptsfile,"w")
709        file.write(" x,y," + att + " \n")
710        for data_point, attribute in map(None, data_points_absolute
711                                         ,attributes):
712            row = str(data_point[0]) + ',' + str(data_point[1]) \
713                  + ',' + str(attribute)
714            file.write(row + "\n")
715        file.close()
716
717
718        #Check that values can be set from file
719        quantity.set_values(filename = ptsfile,
720                            attribute_name = att, alpha = 0, use_cache=True,
721                            verbose=True)
722        answer = linear_function(quantity.domain.get_vertex_coordinates())
723
724        #print quantity.vertex_values.flat
725        #print answer
726
727
728        assert allclose(quantity.vertex_values.flat, answer)
729
730
731        #Check that values can be set from file using default attribute
732        quantity.set_values(filename = ptsfile, alpha = 0)
733        assert allclose(quantity.vertex_values.flat, answer)
734
735        #checking cache
736        #quantity.set_values(filename = ptsfile,
737         #                   attribute_name = att, alpha = 0, use_cache=True,
738          #                  verbose=True)
739        #Cleanup
740        import os
741        os.remove(ptsfile)
742
743    def test_set_values_from_lat_long(self):
744        quantity = Quantity(self.mesh_onslow)
745
746        #Get (enough) datapoints
747        data_points = [[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
748                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6]]
749
750        data_geo_spatial = Geospatial_data(data_points,
751                                           points_are_lats_longs=True)
752        points_UTM = data_geo_spatial.get_data_points(absolute=True)
753        attributes = linear_function(points_UTM)
754        att = 'elevation'
755       
756        #Create .txt file
757        txt_file = tempfile.mktemp(".txt")
758        file = open(txt_file,"w")
759        file.write(" lat,long," + att + " \n")
760        for data_point, attribute in map(None, data_points, attributes):
761            row = str(data_point[0]) + ',' + str(data_point[1]) \
762                  + ',' + str(attribute)
763            #print "row", row
764            file.write(row + "\n")
765        file.close()
766
767
768        #Check that values can be set from file
769        quantity.set_values(filename = txt_file,
770                            attribute_name = att, alpha = 0)
771        answer = linear_function(quantity.domain.get_vertex_coordinates())
772
773        #print "quantity.vertex_values.flat", quantity.vertex_values.flat
774        #print "answer",answer
775
776        assert allclose(quantity.vertex_values.flat, answer)
777
778
779        #Check that values can be set from file using default attribute
780        quantity.set_values(filename = txt_file, alpha = 0)
781        assert allclose(quantity.vertex_values.flat, answer)
782
783        #Cleanup
784        import os
785        os.remove(txt_file)
786         
787    def test_set_values_from_lat_long(self):
788        quantity = Quantity(self.mesh_onslow)
789
790        #Get (enough) datapoints
791        data_points = [[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
792                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6]]
793
794        data_geo_spatial = Geospatial_data(data_points,
795                                           points_are_lats_longs=True)
796        points_UTM = data_geo_spatial.get_data_points(absolute=True)
797        attributes = linear_function(points_UTM)
798        att = 'elevation'
799       
800        #Create .txt file
801        txt_file = tempfile.mktemp(".txt")
802        file = open(txt_file,"w")
803        file.write(" lat,long," + att + " \n")
804        for data_point, attribute in map(None, data_points, attributes):
805            row = str(data_point[0]) + ',' + str(data_point[1]) \
806                  + ',' + str(attribute)
807            #print "row", row
808            file.write(row + "\n")
809        file.close()
810
811
812        #Check that values can be set from file
813        quantity.set_values(filename = txt_file,
814                            attribute_name = att, alpha = 0)
815        answer = linear_function(quantity.domain.get_vertex_coordinates())
816
817        #print "quantity.vertex_values.flat", quantity.vertex_values.flat
818        #print "answer",answer
819
820        assert allclose(quantity.vertex_values.flat, answer)
821
822
823        #Check that values can be set from file using default attribute
824        quantity.set_values(filename = txt_file, alpha = 0)
825        assert allclose(quantity.vertex_values.flat, answer)
826
827        #Cleanup
828        import os
829        os.remove(txt_file)
830       
831    def test_set_values_from_UTM_pts(self):
832        quantity = Quantity(self.mesh_onslow)
833
834        #Get (enough) datapoints
835        data_points = [[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
836                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6]]
837
838        data_geo_spatial = Geospatial_data(data_points,
839                                           points_are_lats_longs=True)
840        points_UTM = data_geo_spatial.get_data_points(absolute=True)
841        attributes = linear_function(points_UTM)
842        att = 'elevation'
843       
844        #Create .txt file
845        txt_file = tempfile.mktemp(".txt")
846        file = open(txt_file,"w")
847        file.write(" x,y," + att + " \n")
848        for data_point, attribute in map(None, points_UTM, attributes):
849            row = str(data_point[0]) + ',' + str(data_point[1]) \
850                  + ',' + str(attribute)
851            #print "row", row
852            file.write(row + "\n")
853        file.close()
854
855
856        pts_file = tempfile.mktemp(".pts")       
857        convert = Geospatial_data(txt_file)
858        convert.export_points_file(pts_file)
859
860        #Check that values can be set from file
861        quantity.set_values_from_file(pts_file, att, 0,
862                                      'vertices', None)
863        answer = linear_function(quantity.domain.get_vertex_coordinates())
864        #print "quantity.vertex_values.flat", quantity.vertex_values.flat
865        #print "answer",answer
866        assert allclose(quantity.vertex_values.flat, answer)
867
868        #Check that values can be set from file
869        quantity.set_values(filename = pts_file,
870                            attribute_name = att, alpha = 0)
871        answer = linear_function(quantity.domain.get_vertex_coordinates())
872        #print "quantity.vertex_values.flat", quantity.vertex_values.flat
873        #print "answer",answer
874        assert allclose(quantity.vertex_values.flat, answer)
875
876
877        #Check that values can be set from file using default attribute
878        quantity.set_values(filename = txt_file, alpha = 0)
879        assert allclose(quantity.vertex_values.flat, answer)
880
881        #Cleanup
882        import os
883        os.remove(txt_file)
884        os.remove(pts_file)
885       
886    def test_set_values_from_UTM_pts(self):
887        quantity = Quantity(self.mesh_onslow)
888
889        VERBOSE = False # Change this to True to see output
890        #Get (enough) datapoints
891        data_points = [[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
892                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
893                       [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
894                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
895                       [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
896                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
897                       [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
898                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
899                       [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
900                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
901                       [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
902                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
903                       [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
904                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
905                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
906                       [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
907                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
908                       [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
909                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
910                       [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
911                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
912                       [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
913                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
914                       [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
915                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
916                       [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
917                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
918                       ]
919
920        data_geo_spatial = Geospatial_data(data_points,
921                                           points_are_lats_longs=True)
922        points_UTM = data_geo_spatial.get_data_points(absolute=True)
923        attributes = linear_function(points_UTM)
924        att = 'elevation'
925       
926        #Create .txt file
927        txt_file = tempfile.mktemp(".txt")
928        file = open(txt_file,"w")
929        file.write(" x,y," + att + " \n")
930        for data_point, attribute in map(None, points_UTM, attributes):
931            row = str(data_point[0]) + ',' + str(data_point[1]) \
932                  + ',' + str(attribute)
933            #print "row", row
934            file.write(row + "\n")
935        file.close()
936
937
938        pts_file = tempfile.mktemp(".pts")       
939        convert = Geospatial_data(txt_file)
940        convert.export_points_file(pts_file)
941
942        #Check that values can be set from file
943        quantity.set_values_from_file(pts_file, att, 0,
944                                      'vertices', None, verbose = VERBOSE,
945                                      max_read_lines=2)
946        answer = linear_function(quantity.domain.get_vertex_coordinates())
947        #print "quantity.vertex_values.flat", quantity.vertex_values.flat
948        #print "answer",answer
949        assert allclose(quantity.vertex_values.flat, answer)
950
951        #Check that values can be set from file
952        quantity.set_values(filename = pts_file, verbose = VERBOSE,
953                            attribute_name = att, alpha = 0)
954        answer = linear_function(quantity.domain.get_vertex_coordinates())
955        #print "quantity.vertex_values.flat", quantity.vertex_values.flat
956        #print "answer",answer
957        assert allclose(quantity.vertex_values.flat, answer)
958
959
960        #Check that values can be set from file using default attribute
961        quantity.set_values(filename = txt_file, alpha = 0)
962        assert allclose(quantity.vertex_values.flat, answer)
963
964        #Check that values can be block read from a text file
965        quantity.set_values_from_file(txt_file, att, 0,
966                                      'vertices', None, verbose = VERBOSE,
967                                      max_read_lines=2)
968       
969        #Cleanup
970        import os
971        os.remove(txt_file)
972        os.remove(pts_file)
973       
974    def test_set_values_from_file_with_georef1(self):
975
976        #Mesh in zone 56 (absolute coords)
977
978        x0 = 314036.58727982
979        y0 = 6224951.2960092
980
981        a = [x0+0.0, y0+0.0]
982        b = [x0+0.0, y0+2.0]
983        c = [x0+2.0, y0+0.0]
984        d = [x0+0.0, y0+4.0]
985        e = [x0+2.0, y0+2.0]
986        f = [x0+4.0, y0+0.0]
987
988        points = [a, b, c, d, e, f]
989
990        #bac, bce, ecf, dbe
991        elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
992
993        #absolute going in ..
994        mesh4 = Domain(points, elements,
995                       geo_reference = Geo_reference(56, 0, 0))
996        mesh4.check_integrity()
997        quantity = Quantity(mesh4)
998
999        #Get (enough) datapoints (relative to georef)
1000        data_points_rel = [[ 0.66666667, 0.66666667],
1001                       [ 1.33333333, 1.33333333],
1002                       [ 2.66666667, 0.66666667],
1003                       [ 0.66666667, 2.66666667],
1004                       [ 0.0, 1.0],
1005                       [ 0.0, 3.0],
1006                       [ 1.0, 0.0],
1007                       [ 1.0, 1.0],
1008                       [ 1.0, 2.0],
1009                       [ 1.0, 3.0],
1010                       [ 2.0, 1.0],
1011                       [ 3.0, 0.0],
1012                       [ 3.0, 1.0]]
1013
1014        data_geo_spatial = Geospatial_data(data_points_rel,
1015                         geo_reference = Geo_reference(56, x0, y0))
1016        data_points_absolute = data_geo_spatial.get_data_points(absolute=True)
1017        attributes = linear_function(data_points_absolute)
1018        att = 'spam_and_eggs'
1019       
1020        #Create .txt file
1021        ptsfile = tempfile.mktemp(".txt")
1022        file = open(ptsfile,"w")
1023        file.write(" x,y," + att + " \n")
1024        for data_point, attribute in map(None, data_points_absolute
1025                                         ,attributes):
1026            row = str(data_point[0]) + ',' + str(data_point[1]) \
1027                  + ',' + str(attribute)
1028            file.write(row + "\n")
1029        file.close()
1030
1031        #file = open(ptsfile, 'r')
1032        #lines = file.readlines()
1033        #file.close()
1034     
1035
1036        #Check that values can be set from file
1037        quantity.set_values(filename = ptsfile,
1038                            attribute_name = att, alpha = 0)
1039        answer = linear_function(quantity.domain.get_vertex_coordinates())
1040
1041        assert allclose(quantity.vertex_values.flat, answer)
1042
1043
1044        #Check that values can be set from file using default attribute
1045        quantity.set_values(filename = ptsfile, alpha = 0)
1046        assert allclose(quantity.vertex_values.flat, answer)
1047
1048        #Cleanup
1049        import os
1050        os.remove(ptsfile)
1051
1052
1053    def test_set_values_from_file_with_georef2(self):
1054
1055        #Mesh in zone 56 (relative coords)
1056
1057        x0 = 314036.58727982
1058        y0 = 6224951.2960092
1059        #x0 = 0.0
1060        #y0 = 0.0
1061
1062        a = [0.0, 0.0]
1063        b = [0.0, 2.0]
1064        c = [2.0, 0.0]
1065        d = [0.0, 4.0]
1066        e = [2.0, 2.0]
1067        f = [4.0, 0.0]
1068
1069        points = [a, b, c, d, e, f]
1070
1071        #bac, bce, ecf, dbe
1072        elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
1073
1074        mesh4 = Domain(points, elements,
1075                       geo_reference = Geo_reference(56, x0, y0))
1076        mesh4.check_integrity()
1077        quantity = Quantity(mesh4)
1078
1079        #Get (enough) datapoints
1080        data_points = [[ x0+0.66666667, y0+0.66666667],
1081                       [ x0+1.33333333, y0+1.33333333],
1082                       [ x0+2.66666667, y0+0.66666667],
1083                       [ x0+0.66666667, y0+2.66666667],
1084                       [ x0+0.0, y0+1.0],
1085                       [ x0+0.0, y0+3.0],
1086                       [ x0+1.0, y0+0.0],
1087                       [ x0+1.0, y0+1.0],
1088                       [ x0+1.0, y0+2.0],
1089                       [ x0+1.0, y0+3.0],
1090                       [ x0+2.0, y0+1.0],
1091                       [ x0+3.0, y0+0.0],
1092                       [ x0+3.0, y0+1.0]]
1093
1094
1095        data_geo_spatial = Geospatial_data(data_points,
1096                         geo_reference = Geo_reference(56, 0, 0))
1097        data_points_absolute = data_geo_spatial.get_data_points(absolute=True)
1098        attributes = linear_function(data_points_absolute)
1099        att = 'spam_and_eggs'
1100       
1101        #Create .txt file
1102        ptsfile = tempfile.mktemp(".txt")
1103        file = open(ptsfile,"w")
1104        file.write(" x,y," + att + " \n")
1105        for data_point, attribute in map(None, data_points_absolute
1106                                         ,attributes):
1107            row = str(data_point[0]) + ',' + str(data_point[1]) \
1108                  + ',' + str(attribute)
1109            file.write(row + "\n")
1110        file.close()
1111
1112
1113        #Check that values can be set from file
1114        quantity.set_values(filename = ptsfile,
1115                            attribute_name = att, alpha = 0)
1116        answer = linear_function(quantity.domain. \
1117                                 get_vertex_coordinates(absolute=True))
1118
1119
1120        assert allclose(quantity.vertex_values.flat, answer)
1121
1122
1123        #Check that values can be set from file using default attribute
1124        quantity.set_values(filename = ptsfile, alpha = 0)
1125        assert allclose(quantity.vertex_values.flat, answer)
1126
1127        #Cleanup
1128        import os
1129        os.remove(ptsfile)
1130
1131
1132
1133
1134    def test_set_values_from_quantity(self):
1135
1136        quantity1 = Quantity(self.mesh4)
1137        quantity1.set_vertex_values([0,1,2,3,4,5])
1138
1139        assert allclose(quantity1.vertex_values,
1140                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
1141
1142
1143        quantity2 = Quantity(self.mesh4)
1144        quantity2.set_values(quantity = quantity1)
1145        assert allclose(quantity2.vertex_values,
1146                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
1147
1148        quantity2.set_values(quantity = 2*quantity1)
1149        assert allclose(quantity2.vertex_values,
1150                        [[2,0,4], [2,4,8], [8,4,10], [6,2,8]])
1151
1152        quantity2.set_values(quantity = 2*quantity1 + 3)
1153        assert allclose(quantity2.vertex_values,
1154                        [[5,3,7], [5,7,11], [11,7,13], [9,5,11]])
1155
1156
1157        #Check detection of quantity as first orgument
1158        quantity2.set_values(2*quantity1 + 3)
1159        assert allclose(quantity2.vertex_values,
1160                        [[5,3,7], [5,7,11], [11,7,13], [9,5,11]])
1161
1162
1163
1164
1165
1166    def test_overloading(self):
1167
1168        quantity1 = Quantity(self.mesh4)
1169        quantity1.set_vertex_values([0,1,2,3,4,5])
1170
1171        assert allclose(quantity1.vertex_values,
1172                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
1173
1174
1175        quantity2 = Quantity(self.mesh4)
1176        quantity2.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
1177                             location = 'vertices')
1178
1179
1180
1181        quantity3 = Quantity(self.mesh4)
1182        quantity3.set_values([[2,2,2], [7,8,9], [7,6,3], [3, 8, -8]],
1183                             location = 'vertices')
1184
1185
1186        #Negation
1187        Q = -quantity1
1188        assert allclose(Q.vertex_values, -quantity1.vertex_values)
1189        assert allclose(Q.centroid_values, -quantity1.centroid_values)
1190        assert allclose(Q.edge_values, -quantity1.edge_values)
1191
1192        #Addition
1193        Q = quantity1 + 7
1194        assert allclose(Q.vertex_values, quantity1.vertex_values + 7)
1195        assert allclose(Q.centroid_values, quantity1.centroid_values + 7)
1196        assert allclose(Q.edge_values, quantity1.edge_values + 7)
1197
1198        Q = 7 + quantity1
1199        assert allclose(Q.vertex_values, quantity1.vertex_values + 7)
1200        assert allclose(Q.centroid_values, quantity1.centroid_values + 7)
1201        assert allclose(Q.edge_values, quantity1.edge_values + 7)
1202
1203        Q = quantity1 + quantity2
1204        assert allclose(Q.vertex_values,
1205                        quantity1.vertex_values + quantity2.vertex_values)
1206        assert allclose(Q.centroid_values,
1207                        quantity1.centroid_values + quantity2.centroid_values)
1208        assert allclose(Q.edge_values,
1209                        quantity1.edge_values + quantity2.edge_values)
1210
1211
1212        Q = quantity1 + quantity2 - 3
1213        assert allclose(Q.vertex_values,
1214                        quantity1.vertex_values + quantity2.vertex_values - 3)
1215
1216        Q = quantity1 - quantity2
1217        assert allclose(Q.vertex_values,
1218                        quantity1.vertex_values - quantity2.vertex_values)
1219
1220        #Scaling
1221        Q = quantity1*3
1222        assert allclose(Q.vertex_values, quantity1.vertex_values*3)
1223        assert allclose(Q.centroid_values, quantity1.centroid_values*3)
1224        assert allclose(Q.edge_values, quantity1.edge_values*3)
1225        Q = 3*quantity1
1226        assert allclose(Q.vertex_values, quantity1.vertex_values*3)
1227
1228        #Multiplication
1229        Q = quantity1 * quantity2
1230        #print Q.vertex_values
1231        #print Q.centroid_values
1232        #print quantity1.centroid_values
1233        #print quantity2.centroid_values
1234
1235        assert allclose(Q.vertex_values,
1236                        quantity1.vertex_values * quantity2.vertex_values)
1237
1238        #Linear combinations
1239        Q = 4*quantity1 + 2
1240        assert allclose(Q.vertex_values,
1241                        4*quantity1.vertex_values + 2)
1242
1243        Q = quantity1*quantity2 + 2
1244        assert allclose(Q.vertex_values,
1245                        quantity1.vertex_values * quantity2.vertex_values + 2)
1246
1247        Q = quantity1*quantity2 + quantity3
1248        assert allclose(Q.vertex_values,
1249                        quantity1.vertex_values * quantity2.vertex_values +
1250                        quantity3.vertex_values)
1251        Q = quantity1*quantity2 + 3*quantity3
1252        assert allclose(Q.vertex_values,
1253                        quantity1.vertex_values * quantity2.vertex_values +
1254                        3*quantity3.vertex_values)
1255        Q = quantity1*quantity2 + 3*quantity3 + 5.0
1256        assert allclose(Q.vertex_values,
1257                        quantity1.vertex_values * quantity2.vertex_values +
1258                        3*quantity3.vertex_values + 5)
1259
1260        Q = quantity1*quantity2 - quantity3
1261        assert allclose(Q.vertex_values,
1262                        quantity1.vertex_values * quantity2.vertex_values -
1263                        quantity3.vertex_values)
1264        Q = 1.5*quantity1*quantity2 - 3*quantity3 + 5.0
1265        assert allclose(Q.vertex_values,
1266                        1.5*quantity1.vertex_values * quantity2.vertex_values -
1267                        3*quantity3.vertex_values + 5)
1268
1269        #Try combining quantities and arrays and scalars
1270        Q = 1.5*quantity1*quantity2.vertex_values -\
1271            3*quantity3.vertex_values + 5.0
1272        assert allclose(Q.vertex_values,
1273                        1.5*quantity1.vertex_values * quantity2.vertex_values -
1274                        3*quantity3.vertex_values + 5)
1275
1276
1277        #Powers
1278        Q = quantity1**2
1279        assert allclose(Q.vertex_values, quantity1.vertex_values**2)
1280
1281        Q = quantity1**2 +quantity2**2
1282        assert allclose(Q.vertex_values,
1283                        quantity1.vertex_values**2 + \
1284                        quantity2.vertex_values**2)
1285
1286        Q = (quantity1**2 +quantity2**2)**0.5
1287        assert allclose(Q.vertex_values,
1288                        (quantity1.vertex_values**2 + \
1289                        quantity2.vertex_values**2)**0.5)
1290
1291
1292
1293
1294
1295
1296
1297    def test_compute_gradient(self):
1298        quantity = Quantity(self.mesh4)
1299
1300        #Set up for a gradient of (2,0) at mid triangle
1301        quantity.set_values([2.0, 4.0, 6.0, 2.0],
1302                            location = 'centroids')
1303
1304
1305        #Gradients
1306        quantity.compute_gradients()
1307
1308        a = quantity.x_gradient
1309        b = quantity.y_gradient
1310        #print self.mesh4.centroid_coordinates
1311        #print a, b
1312
1313        #The central triangle (1)
1314        #(using standard gradient based on neigbours controid values)
1315        assert allclose(a[1], 2.0)
1316        assert allclose(b[1], 0.0)
1317
1318
1319        #Left triangle (0) using two point gradient
1320        #q0 = q1 + a*(x0-x1) + b*(y0-y1)  <=>
1321        #2  = 4  + a*(-2/3)  + b*(-2/3)
1322        assert allclose(a[0] + b[0], 3)
1323        #From orthogonality (a*(y0-y1) + b*(x0-x1) == 0)
1324        assert allclose(a[0] - b[0], 0)
1325
1326
1327        #Right triangle (2) using two point gradient
1328        #q2 = q1 + a*(x2-x1) + b*(y2-y1)  <=>
1329        #6  = 4  + a*(4/3)  + b*(-2/3)
1330        assert allclose(2*a[2] - b[2], 3)
1331        #From orthogonality (a*(y1-y2) + b*(x2-x1) == 0)
1332        assert allclose(a[2] + 2*b[2], 0)
1333
1334
1335        #Top triangle (3) using two point gradient
1336        #q3 = q1 + a*(x3-x1) + b*(y3-y1)  <=>
1337        #2  = 4  + a*(-2/3)  + b*(4/3)
1338        assert allclose(a[3] - 2*b[3], 3)
1339        #From orthogonality (a*(y1-y3) + b*(x3-x1) == 0)
1340        assert allclose(2*a[3] + b[3], 0)
1341
1342
1343
1344        #print a, b
1345        quantity.extrapolate_second_order()
1346
1347        #Apply q(x,y) = qc + a*(x-xc) + b*(y-yc)
1348        assert allclose(quantity.vertex_values[0,:], [3., 0.,  3.])
1349        assert allclose(quantity.vertex_values[1,:], [4./3, 16./3,  16./3])
1350
1351
1352        #a = 1.2, b=-0.6
1353        #q(4,0) = 6 + a*(4 - 8/3) + b*(-2/3)
1354        assert allclose(quantity.vertex_values[2,2], 8)
1355
1356    def test_get_gradients(self):
1357        quantity = Quantity(self.mesh4)
1358
1359        #Set up for a gradient of (2,0) at mid triangle
1360        quantity.set_values([2.0, 4.0, 6.0, 2.0],
1361                            location = 'centroids')
1362
1363
1364        #Gradients
1365        quantity.compute_gradients()
1366
1367        a, b = quantity.get_gradients()
1368        #print self.mesh4.centroid_coordinates
1369        #print a, b
1370
1371        #The central triangle (1)
1372        #(using standard gradient based on neigbours controid values)
1373        assert allclose(a[1], 2.0)
1374        assert allclose(b[1], 0.0)
1375
1376
1377        #Left triangle (0) using two point gradient
1378        #q0 = q1 + a*(x0-x1) + b*(y0-y1)  <=>
1379        #2  = 4  + a*(-2/3)  + b*(-2/3)
1380        assert allclose(a[0] + b[0], 3)
1381        #From orthogonality (a*(y0-y1) + b*(x0-x1) == 0)
1382        assert allclose(a[0] - b[0], 0)
1383
1384
1385        #Right triangle (2) using two point gradient
1386        #q2 = q1 + a*(x2-x1) + b*(y2-y1)  <=>
1387        #6  = 4  + a*(4/3)  + b*(-2/3)
1388        assert allclose(2*a[2] - b[2], 3)
1389        #From orthogonality (a*(y1-y2) + b*(x2-x1) == 0)
1390        assert allclose(a[2] + 2*b[2], 0)
1391
1392
1393        #Top triangle (3) using two point gradient
1394        #q3 = q1 + a*(x3-x1) + b*(y3-y1)  <=>
1395        #2  = 4  + a*(-2/3)  + b*(4/3)
1396        assert allclose(a[3] - 2*b[3], 3)
1397        #From orthogonality (a*(y1-y3) + b*(x3-x1) == 0)
1398        assert allclose(2*a[3] + b[3], 0)
1399
1400
1401    def test_second_order_extrapolation2(self):
1402        quantity = Quantity(self.mesh4)
1403
1404        #Set up for a gradient of (3,1), f(x) = 3x+y
1405        quantity.set_values([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3],
1406                            location = 'centroids')
1407
1408        #Gradients
1409        quantity.compute_gradients()
1410
1411        a = quantity.x_gradient
1412        b = quantity.y_gradient
1413       
1414        #print a, b
1415
1416        assert allclose(a[1], 3.0)
1417        assert allclose(b[1], 1.0)
1418
1419        #Work out the others
1420
1421        quantity.extrapolate_second_order()
1422
1423        #print quantity.vertex_values
1424        assert allclose(quantity.vertex_values[1,0], 2.0)
1425        assert allclose(quantity.vertex_values[1,1], 6.0)
1426        assert allclose(quantity.vertex_values[1,2], 8.0)
1427
1428
1429
1430    def test_backup_saxpy_centroid_values(self):
1431        quantity = Quantity(self.mesh4)
1432
1433        #Set up for a gradient of (3,1), f(x) = 3x+y
1434        c_values = array([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3])
1435        d_values = array([1.0, 2.0, 3.0, 4.0])
1436        quantity.set_values(c_values, location = 'centroids')
1437
1438        #Backup
1439        quantity.backup_centroid_values()
1440
1441        #print quantity.vertex_values
1442        assert allclose(quantity.centroid_values, quantity.centroid_backup_values)
1443
1444
1445        quantity.set_values(d_values, location = 'centroids')
1446
1447        quantity.saxpy_centroid_values(2.0, 3.0)
1448
1449        assert(quantity.centroid_values, 2.0*d_values + 3.0*c_values)
1450
1451
1452
1453    def test_first_order_extrapolator(self):
1454        quantity = Quantity(self.mesh4)
1455
1456        #Test centroids
1457        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
1458        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
1459
1460        #Extrapolate
1461        quantity.extrapolate_first_order()
1462
1463        #Check that gradient is zero
1464        a,b = quantity.get_gradients()
1465        assert allclose(a, [0,0,0,0])
1466        assert allclose(b, [0,0,0,0])
1467
1468        #Check vertices but not edge values
1469        assert allclose(quantity.vertex_values,
1470                        [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
1471
1472
1473    def test_second_order_extrapolator(self):
1474        quantity = Quantity(self.mesh4)
1475
1476        #Set up for a gradient of (3,0) at mid triangle
1477        quantity.set_values([2.0, 4.0, 8.0, 2.0],
1478                            location = 'centroids')
1479
1480
1481
1482        quantity.extrapolate_second_order()
1483        quantity.limit()
1484
1485
1486        #Assert that central triangle is limited by neighbours
1487        assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0]
1488        assert quantity.vertex_values[1,0] >= quantity.vertex_values[3,1]
1489
1490        assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1]
1491        assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
1492
1493        assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
1494        assert quantity.vertex_values[1,2] >= quantity.vertex_values[3,1]
1495
1496
1497        #Assert that quantities are conserved
1498        from Numeric import sum
1499        for k in range(quantity.centroid_values.shape[0]):
1500            assert allclose (quantity.centroid_values[k],
1501                             sum(quantity.vertex_values[k,:])/3)
1502
1503
1504
1505
1506
1507    def test_limit_vertices_by_all_neighbours(self):
1508        quantity = Quantity(self.mesh4)
1509
1510        #Create a deliberate overshoot (e.g. from gradient computation)
1511        quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
1512
1513
1514        #Limit
1515        quantity.limit_vertices_by_all_neighbours()
1516
1517        #Assert that central triangle is limited by neighbours
1518        assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0]
1519        assert quantity.vertex_values[1,0] <= quantity.vertex_values[3,1]
1520
1521        assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1]
1522        assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
1523
1524        assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
1525        assert quantity.vertex_values[1,2] <= quantity.vertex_values[3,1]
1526
1527
1528
1529        #Assert that quantities are conserved
1530        from Numeric import sum
1531        for k in range(quantity.centroid_values.shape[0]):
1532            assert allclose (quantity.centroid_values[k],
1533                             sum(quantity.vertex_values[k,:])/3)
1534
1535
1536
1537    def test_limit_edges_by_all_neighbours(self):
1538        quantity = Quantity(self.mesh4)
1539
1540        #Create a deliberate overshoot (e.g. from gradient computation)
1541        quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
1542
1543
1544        #Limit
1545        quantity.limit_edges_by_all_neighbours()
1546
1547        #Assert that central triangle is limited by neighbours
1548        assert quantity.edge_values[1,0] <= quantity.centroid_values[2]
1549        assert quantity.edge_values[1,0] >= quantity.centroid_values[0]
1550
1551        assert quantity.edge_values[1,1] <= quantity.centroid_values[2]
1552        assert quantity.edge_values[1,1] >= quantity.centroid_values[0]
1553
1554        assert quantity.edge_values[1,2] <= quantity.centroid_values[2]
1555        assert quantity.edge_values[1,2] >= quantity.centroid_values[0]
1556
1557
1558
1559        #Assert that quantities are conserved
1560        from Numeric import sum
1561        for k in range(quantity.centroid_values.shape[0]):
1562            assert allclose (quantity.centroid_values[k],
1563                             sum(quantity.vertex_values[k,:])/3)
1564
1565
1566    def test_limit_edges_by_neighbour(self):
1567        quantity = Quantity(self.mesh4)
1568
1569        #Create a deliberate overshoot (e.g. from gradient computation)
1570        quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
1571
1572
1573        #Limit
1574        quantity.limit_edges_by_neighbour()
1575
1576        #Assert that central triangle is limited by neighbours
1577        assert quantity.edge_values[1,0] <= quantity.centroid_values[3]
1578        assert quantity.edge_values[1,0] >= quantity.centroid_values[1]
1579
1580        assert quantity.edge_values[1,1] <= quantity.centroid_values[2]
1581        assert quantity.edge_values[1,1] >= quantity.centroid_values[1]
1582
1583        assert quantity.edge_values[1,2] <= quantity.centroid_values[1]
1584        assert quantity.edge_values[1,2] >= quantity.centroid_values[0]
1585
1586
1587
1588        #Assert that quantities are conserved
1589        from Numeric import sum
1590        for k in range(quantity.centroid_values.shape[0]):
1591            assert allclose (quantity.centroid_values[k],
1592                             sum(quantity.vertex_values[k,:])/3)
1593
1594    def test_limiter2(self):
1595        """Taken from test_shallow_water
1596        """
1597        quantity = Quantity(self.mesh4)
1598        quantity.domain.beta_w = 0.9
1599       
1600        #Test centroids
1601        quantity.set_values([2.,4.,8.,2.], location = 'centroids')
1602        assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid
1603
1604
1605        #Extrapolate
1606        quantity.extrapolate_second_order()
1607
1608        assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6])
1609
1610        #Limit
1611        quantity.limit()
1612
1613        # limited value for beta_w = 0.9
1614       
1615        assert allclose(quantity.vertex_values[1,:], [2.2, 4.9, 4.9])
1616        # limited values for beta_w = 0.5
1617        #assert allclose(quantity.vertex_values[1,:], [3.0, 4.5, 4.5])
1618
1619
1620        #Assert that quantities are conserved
1621        from Numeric import sum
1622        for k in range(quantity.centroid_values.shape[0]):
1623            assert allclose (quantity.centroid_values[k],
1624                             sum(quantity.vertex_values[k,:])/3)
1625
1626
1627
1628
1629
1630    def test_distribute_first_order(self):
1631        quantity = Quantity(self.mesh4)
1632
1633        #Test centroids
1634        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
1635        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
1636
1637
1638        #Extrapolate from centroid to vertices and edges
1639        quantity.extrapolate_first_order()
1640
1641        #Interpolate
1642        #quantity.interpolate_from_vertices_to_edges()
1643
1644        assert allclose(quantity.vertex_values,
1645                        [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
1646        assert allclose(quantity.edge_values, [[1,1,1], [2,2,2],
1647                                               [3,3,3], [4, 4, 4]])
1648
1649
1650    def test_interpolate_from_vertices_to_edges(self):
1651        quantity = Quantity(self.mesh4)
1652
1653        quantity.vertex_values = array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]],Float)
1654
1655        quantity.interpolate_from_vertices_to_edges()
1656
1657        assert allclose(quantity.edge_values, [[1., 1.5, 0.5],
1658                                               [3., 2.5, 1.5],
1659                                               [3.5, 4.5, 3.],
1660                                               [2.5, 3.5, 2]])
1661
1662
1663    def test_interpolate_from_edges_to_vertices(self):
1664        quantity = Quantity(self.mesh4)
1665
1666        quantity.edge_values = array([[1., 1.5, 0.5],
1667                                [3., 2.5, 1.5],
1668                                [3.5, 4.5, 3.],
1669                                [2.5, 3.5, 2]],Float)
1670
1671        quantity.interpolate_from_edges_to_vertices()
1672
1673        assert allclose(quantity.vertex_values,
1674                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
1675
1676
1677
1678    def test_distribute_second_order(self):
1679        quantity = Quantity(self.mesh4)
1680
1681        #Test centroids
1682        quantity.set_values([2.,4.,8.,2.], location = 'centroids')
1683        assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid
1684
1685
1686        #Extrapolate
1687        quantity.extrapolate_second_order()
1688
1689        assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6])
1690
1691
1692    def test_update_explicit(self):
1693        quantity = Quantity(self.mesh4)
1694
1695        #Test centroids
1696        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
1697        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
1698
1699        #Set explicit_update
1700        quantity.explicit_update = array( [1.,1.,1.,1.] )
1701
1702        #Update with given timestep
1703        quantity.update(0.1)
1704
1705        x = array([1, 2, 3, 4]) + array( [.1,.1,.1,.1] )
1706        assert allclose( quantity.centroid_values, x)
1707
1708    def test_update_semi_implicit(self):
1709        quantity = Quantity(self.mesh4)
1710
1711        #Test centroids
1712        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
1713        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
1714
1715        #Set semi implicit update
1716        quantity.semi_implicit_update = array([1.,1.,1.,1.])
1717
1718        #Update with given timestep
1719        timestep = 0.1
1720        quantity.update(timestep)
1721
1722        sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])
1723        denom = ones(4, Float)-timestep*sem
1724
1725        x = array([1, 2, 3, 4])/denom
1726        assert allclose( quantity.centroid_values, x)
1727
1728
1729    def test_both_updates(self):
1730        quantity = Quantity(self.mesh4)
1731
1732        #Test centroids
1733        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
1734        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
1735
1736        #Set explicit_update
1737        quantity.explicit_update = array( [4.,3.,2.,1.] )
1738
1739        #Set semi implicit update
1740        quantity.semi_implicit_update = array( [1.,1.,1.,1.] )
1741
1742        #Update with given timestep
1743        timestep = 0.1
1744        quantity.update(0.1)
1745
1746        sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])
1747        denom = ones(4, Float)-timestep*sem
1748
1749        x = array([1., 2., 3., 4.])
1750        x /= denom
1751        x += timestep*array( [4.0, 3.0, 2.0, 1.0] )
1752
1753        assert allclose( quantity.centroid_values, x)
1754
1755
1756
1757
1758    #Test smoothing
1759    def test_smoothing(self):
1760
1761        from mesh_factory import rectangular
1762        from shallow_water import Domain, Transmissive_boundary
1763        from Numeric import zeros, Float
1764        from anuga.utilities.numerical_tools import mean
1765
1766        #Create basic mesh
1767        points, vertices, boundary = rectangular(2, 2)
1768
1769        #Create shallow water domain
1770        domain = Domain(points, vertices, boundary)
1771        domain.default_order=2
1772        domain.reduction = mean
1773
1774
1775        #Set some field values
1776        domain.set_quantity('elevation', lambda x,y: x)
1777        domain.set_quantity('friction', 0.03)
1778
1779
1780        ######################
1781        # Boundary conditions
1782        B = Transmissive_boundary(domain)
1783        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
1784
1785
1786        ######################
1787        #Initial condition - with jumps
1788
1789        bed = domain.quantities['elevation'].vertex_values
1790        stage = zeros(bed.shape, Float)
1791
1792        h = 0.03
1793        for i in range(stage.shape[0]):
1794            if i % 2 == 0:
1795                stage[i,:] = bed[i,:] + h
1796            else:
1797                stage[i,:] = bed[i,:]
1798
1799        domain.set_quantity('stage', stage)
1800
1801        stage = domain.quantities['stage']
1802
1803        #Get smoothed stage
1804        A, V = stage.get_vertex_values(xy=False, smooth=True)
1805        Q = stage.vertex_values
1806
1807
1808        assert A.shape[0] == 9
1809        assert V.shape[0] == 8
1810        assert V.shape[1] == 3
1811
1812        #First four points
1813        assert allclose(A[0], (Q[0,2] + Q[1,1])/2)
1814        assert allclose(A[1], (Q[1,0] + Q[3,1] + Q[2,2])/3)
1815        assert allclose(A[2], Q[3,0])
1816        assert allclose(A[3], (Q[0,0] + Q[5,1] + Q[4,2])/3)
1817
1818        #Center point
1819        assert allclose(A[4], (Q[0,1] + Q[1,2] + Q[2,0] +\
1820                               Q[5,0] + Q[6,2] + Q[7,1])/6)
1821
1822
1823        #Check V
1824        assert allclose(V[0,:], [3,4,0])
1825        assert allclose(V[1,:], [1,0,4])
1826        assert allclose(V[2,:], [4,5,1])
1827        assert allclose(V[3,:], [2,1,5])
1828        assert allclose(V[4,:], [6,7,3])
1829        assert allclose(V[5,:], [4,3,7])
1830        assert allclose(V[6,:], [7,8,4])
1831        assert allclose(V[7,:], [5,4,8])
1832
1833        #Get smoothed stage with XY
1834        X, Y, A1, V1 = stage.get_vertex_values(xy=True, smooth=True)
1835
1836        assert allclose(A, A1)
1837        assert allclose(V, V1)
1838
1839        #Check XY
1840        assert allclose(X[4], 0.5)
1841        assert allclose(Y[4], 0.5)
1842
1843        assert allclose(X[7], 1.0)
1844        assert allclose(Y[7], 0.5)
1845
1846
1847
1848
1849    def test_vertex_values_no_smoothing(self):
1850
1851        from mesh_factory import rectangular
1852        from shallow_water import Domain, Transmissive_boundary
1853        from Numeric import zeros, Float
1854        from anuga.utilities.numerical_tools import mean
1855
1856
1857        #Create basic mesh
1858        points, vertices, boundary = rectangular(2, 2)
1859
1860        #Create shallow water domain
1861        domain = Domain(points, vertices, boundary)
1862        domain.default_order=2
1863        domain.reduction = mean
1864
1865
1866        #Set some field values
1867        domain.set_quantity('elevation', lambda x,y: x)
1868        domain.set_quantity('friction', 0.03)
1869
1870
1871        ######################
1872        #Initial condition - with jumps
1873
1874        bed = domain.quantities['elevation'].vertex_values
1875        stage = zeros(bed.shape, Float)
1876
1877        h = 0.03
1878        for i in range(stage.shape[0]):
1879            if i % 2 == 0:
1880                stage[i,:] = bed[i,:] + h
1881            else:
1882                stage[i,:] = bed[i,:]
1883
1884        domain.set_quantity('stage', stage)
1885
1886        #Get stage
1887        stage = domain.quantities['stage']
1888        A, V = stage.get_vertex_values(xy=False, smooth=False)
1889        Q = stage.vertex_values.flat
1890
1891        for k in range(8):
1892            assert allclose(A[k], Q[k])
1893
1894
1895        for k in range(8):
1896            assert V[k, 0] == 3*k
1897            assert V[k, 1] == 3*k+1
1898            assert V[k, 2] == 3*k+2
1899
1900
1901
1902        X, Y, A1, V1 = stage.get_vertex_values(xy=True, smooth=False)
1903
1904
1905        assert allclose(A, A1)
1906        assert allclose(V, V1)
1907
1908        #Check XY
1909        assert allclose(X[1], 0.5)
1910        assert allclose(Y[1], 0.5)
1911        assert allclose(X[4], 0.0)
1912        assert allclose(Y[4], 0.0)
1913        assert allclose(X[12], 1.0)
1914        assert allclose(Y[12], 0.0)
1915
1916
1917
1918    def set_array_values_by_index(self):
1919
1920        from mesh_factory import rectangular
1921        from shallow_water import Domain
1922        from Numeric import zeros, Float
1923
1924        #Create basic mesh
1925        points, vertices, boundary = rectangular(1, 1)
1926
1927        #Create shallow water domain
1928        domain = Domain(points, vertices, boundary)
1929        #print "domain.number_of_elements ",domain.number_of_elements
1930        quantity = Quantity(domain,[[1,1,1],[2,2,2]])
1931        value = [7]
1932        indices = [1]
1933        quantity.set_array_values_by_index(value,
1934                                           location = 'centroids',
1935                                           indices = indices)
1936        #print "quantity.centroid_values",quantity.centroid_values
1937
1938        assert allclose(quantity.centroid_values, [1,7])
1939
1940        quantity.set_array_values([15,20,25], indices = indices)
1941        assert allclose(quantity.centroid_values, [1,20])
1942
1943        quantity.set_array_values([15,20,25], indices = indices)
1944        assert allclose(quantity.centroid_values, [1,20])
1945
1946    def test_setting_some_vertex_values(self):
1947        """
1948        set values based on triangle lists.
1949        """
1950        from mesh_factory import rectangular
1951        from shallow_water import Domain
1952        from Numeric import zeros, Float
1953
1954        #Create basic mesh
1955        points, vertices, boundary = rectangular(1, 3)
1956        #print "vertices",vertices
1957        #Create shallow water domain
1958        domain = Domain(points, vertices, boundary)
1959        #print "domain.number_of_elements ",domain.number_of_elements
1960        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
1961                                    [4,4,4],[5,5,5],[6,6,6]])
1962
1963
1964        # Check that constants work
1965        value = 7
1966        indices = [1]
1967        quantity.set_values(value,
1968                            location = 'centroids',
1969                            indices = indices)
1970        #print "quantity.centroid_values",quantity.centroid_values
1971        assert allclose(quantity.centroid_values, [1,7,3,4,5,6])
1972       
1973        value = [7]
1974        indices = [1]
1975        quantity.set_values(value,
1976                            location = 'centroids',
1977                            indices = indices)
1978        #print "quantity.centroid_values",quantity.centroid_values
1979        assert allclose(quantity.centroid_values, [1,7,3,4,5,6])
1980
1981        value = [[15,20,25]]
1982        quantity.set_values(value, indices = indices)
1983        #print "1 quantity.vertex_values",quantity.vertex_values
1984        assert allclose(quantity.vertex_values[1], value[0])
1985
1986
1987        #print "quantity",quantity.vertex_values
1988        values = [10,100,50]
1989        quantity.set_values(values, indices = [0,1,5], location = 'centroids')
1990        #print "2 quantity.vertex_values",quantity.vertex_values
1991        assert allclose(quantity.vertex_values[0], [10,10,10])
1992        assert allclose(quantity.vertex_values[5], [50,50,50])
1993        #quantity.interpolate()
1994        #print "quantity.centroid_values",quantity.centroid_values
1995        assert allclose(quantity.centroid_values, [10,100,3,4,5,50])
1996
1997
1998        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
1999                                    [4,4,4],[5,5,5],[6,6,6]])
2000        values = [10,100,50]
2001        #this will be per unique vertex, indexing the vertices
2002        #print "quantity.vertex_values",quantity.vertex_values
2003        quantity.set_values(values, indices = [0,1,5])
2004        #print "quantity.vertex_values",quantity.vertex_values
2005        assert allclose(quantity.vertex_values[0], [1,50,10])
2006        assert allclose(quantity.vertex_values[5], [6,6,6])
2007        assert allclose(quantity.vertex_values[1], [100,10,50])
2008
2009        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
2010                                    [4,4,4],[5,5,5],[6,6,6]])
2011        values = [[31,30,29],[400,400,400],[1000,999,998]]
2012        quantity.set_values(values, indices = [3,3,5])
2013        quantity.interpolate()
2014        assert allclose(quantity.centroid_values, [1,2,3,400,5,999])
2015
2016        values = [[1,1,1],[2,2,2],[3,3,3],
2017                                    [4,4,4],[5,5,5],[6,6,6]]
2018        quantity.set_values(values)
2019
2020        # testing the standard set values by vertex
2021        # indexed by vertex_id in general_mesh.coordinates
2022        values = [0,1,2,3,4,5,6,7]
2023
2024        quantity.set_values(values)
2025        #print "1 quantity.vertex_values",quantity.vertex_values
2026        assert allclose(quantity.vertex_values,[[ 4.,  5.,  0.],
2027                                                [ 1.,  0.,  5.],
2028                                                [ 5.,  6.,  1.],
2029                                                [ 2.,  1.,  6.],
2030                                                [ 6.,  7.,  2.],
2031                                                [ 3.,  2.,  7.]])
2032
2033    def test_setting_unique_vertex_values(self):
2034        """
2035        set values based on unique_vertex lists.
2036        """
2037        from mesh_factory import rectangular
2038        from shallow_water import Domain
2039        from Numeric import zeros, Float
2040
2041        #Create basic mesh
2042        points, vertices, boundary = rectangular(1, 3)
2043        #print "vertices",vertices
2044        #Create shallow water domain
2045        domain = Domain(points, vertices, boundary)
2046        #print "domain.number_of_elements ",domain.number_of_elements
2047        quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
2048                                    [4,4,4],[5,5,5]])
2049        value = 7
2050        indices = [1,5]
2051        quantity.set_values(value,
2052                            location = 'unique vertices',
2053                            indices = indices)
2054        #print "quantity.centroid_values",quantity.centroid_values
2055        assert allclose(quantity.vertex_values[0], [0,7,0])
2056        assert allclose(quantity.vertex_values[1], [7,1,7])
2057        assert allclose(quantity.vertex_values[2], [7,2,7])
2058
2059
2060    def test_get_values(self):
2061        """
2062        get values based on triangle lists.
2063        """
2064        from mesh_factory import rectangular
2065        from shallow_water import Domain
2066        from Numeric import zeros, Float
2067
2068        #Create basic mesh
2069        points, vertices, boundary = rectangular(1, 3)
2070
2071        #print "points",points
2072        #print "vertices",vertices
2073        #print "boundary",boundary
2074
2075        #Create shallow water domain
2076        domain = Domain(points, vertices, boundary)
2077        #print "domain.number_of_elements ",domain.number_of_elements
2078        quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
2079                                    [4,4,4],[5,5,5]])
2080
2081        #print "quantity.get_values(location = 'unique vertices')", \
2082        #      quantity.get_values(location = 'unique vertices')
2083
2084        #print "quantity.get_values(location = 'unique vertices')", \
2085        #      quantity.get_values(indices=[0,1,2,3,4,5,6,7], \
2086        #                          location = 'unique vertices')
2087
2088        answer = [0.5,2,4,5,0,1,3,4.5]
2089        assert allclose(answer,
2090                        quantity.get_values(location = 'unique vertices'))
2091
2092        indices = [0,5,3]
2093        answer = [0.5,1,5]
2094        assert allclose(answer,
2095                        quantity.get_values(indices=indices, \
2096                                            location = 'unique vertices'))
2097        #print "quantity.centroid_values",quantity.centroid_values
2098        #print "quantity.get_values(location = 'centroids') ",\
2099        #      quantity.get_values(location = 'centroids')
2100
2101
2102
2103
2104    def test_get_values_2(self):
2105        """Different mesh (working with domain object) - also check centroids.
2106        """
2107
2108       
2109        a = [0.0, 0.0]
2110        b = [0.0, 2.0]
2111        c = [2.0,0.0]
2112        d = [0.0, 4.0]
2113        e = [2.0, 2.0]
2114        f = [4.0,0.0]
2115
2116        points = [a, b, c, d, e, f]
2117        #bac, bce, ecf, dbe
2118        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2119
2120        domain = Domain(points, vertices)
2121
2122        quantity = Quantity(domain)
2123        quantity.set_values(lambda x, y: x+2*y) #2 4 4 6
2124       
2125        assert allclose(quantity.get_values(location='centroids'), [2,4,4,6])
2126        assert allclose(quantity.get_values(location='centroids', indices=[1,3]), [4,6])
2127
2128
2129        assert allclose(quantity.get_values(location='vertices'), [[4,0,2],
2130                                                                   [4,2,6],
2131                                                                   [6,2,4],
2132                                                                   [8,4,6]])
2133       
2134        assert allclose(quantity.get_values(location='vertices', indices=[1,3]), [[4,2,6],
2135                                                                                  [8,4,6]])
2136
2137
2138        assert allclose(quantity.get_values(location='edges'), [[1,3,2],
2139                                                                [4,5,3],
2140                                                                [3,5,4],
2141                                                                [5,7,6]])
2142        assert allclose(quantity.get_values(location='edges', indices=[1,3]),
2143                        [[4,5,3],
2144                         [5,7,6]])       
2145
2146        # Check averaging over vertices
2147        #a: 0
2148        #b: (4+4+4)/3
2149        #c: (2+2+2)/3
2150        #d: 8
2151        #e: (6+6+6)/3       
2152        #f: 4
2153        assert(quantity.get_values(location='unique vertices'), [0, 4, 2, 8, 6, 4])       
2154                                                                                 
2155       
2156
2157
2158
2159
2160    def test_get_interpolated_values(self):
2161
2162        from mesh_factory import rectangular
2163        from shallow_water import Domain
2164        from Numeric import zeros, Float
2165
2166        #Create basic mesh
2167        points, vertices, boundary = rectangular(1, 3)
2168        domain = Domain(points, vertices, boundary)
2169
2170        #Constant values
2171        quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
2172                                    [4,4,4],[5,5,5]])
2173
2174       
2175
2176        # Get interpolated values at centroids
2177        interpolation_points = domain.get_centroid_coordinates()
2178        answer = quantity.get_values(location='centroids')
2179
2180       
2181        #print quantity.get_values(points=interpolation_points)
2182        assert allclose(answer, quantity.get_values(interpolation_points=interpolation_points))
2183
2184
2185        #Arbitrary values
2186        quantity = Quantity(domain,[[0,1,2],[3,1,7],[2,1,2],[3,3,7],
2187                                    [1,4,-9],[2,5,0]])
2188
2189
2190        # Get interpolated values at centroids
2191        interpolation_points = domain.get_centroid_coordinates()
2192        answer = quantity.get_values(location='centroids')
2193        #print answer
2194        #print quantity.get_values(points=interpolation_points)
2195        assert allclose(answer, quantity.get_values(interpolation_points=interpolation_points))       
2196                       
2197
2198        #FIXME TODO
2199        #indices = [0,5,3]
2200        #answer = [0.5,1,5]
2201        #assert allclose(answer,
2202        #                quantity.get_values(indices=indices, \
2203        #                                    location = 'unique vertices'))
2204
2205
2206
2207
2208    def test_get_interpolated_values_2(self):
2209        a = [0.0, 0.0]
2210        b = [0.0, 2.0]
2211        c = [2.0,0.0]
2212        d = [0.0, 4.0]
2213        e = [2.0, 2.0]
2214        f = [4.0,0.0]
2215
2216        points = [a, b, c, d, e, f]
2217        #bac, bce, ecf, dbe
2218        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2219
2220        domain = Domain(points, vertices)
2221
2222        quantity = Quantity(domain)
2223        quantity.set_values(lambda x, y: x+2*y) #2 4 4 6
2224
2225        #First pick one point
2226        x, y = 2.0/3, 8.0/3
2227        v = quantity.get_values(interpolation_points = [[x,y]])
2228        assert allclose(v, 6)       
2229
2230        # Then another to test that algorithm won't blindly
2231        # reuse interpolation matrix
2232        x, y = 4.0/3, 4.0/3
2233        v = quantity.get_values(interpolation_points = [[x,y]])
2234        assert allclose(v, 4)       
2235
2236
2237
2238
2239    def test_getting_some_vertex_values(self):
2240        """
2241        get values based on triangle lists.
2242        """
2243        from mesh_factory import rectangular
2244        from shallow_water import Domain
2245        from Numeric import zeros, Float
2246
2247        #Create basic mesh
2248        points, vertices, boundary = rectangular(1, 3)
2249
2250        #print "points",points
2251        #print "vertices",vertices
2252        #print "boundary",boundary
2253
2254        #Create shallow water domain
2255        domain = Domain(points, vertices, boundary)
2256        #print "domain.number_of_elements ",domain.number_of_elements
2257        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
2258                                    [4,4,4],[5,5,5],[6,6,6]])
2259        value = [7]
2260        indices = [1]
2261        quantity.set_values(value,
2262                            location = 'centroids',
2263                            indices = indices)
2264        #print "quantity.centroid_values",quantity.centroid_values
2265        #print "quantity.get_values(location = 'centroids') ",\
2266        #      quantity.get_values(location = 'centroids')
2267        assert allclose(quantity.centroid_values,
2268                        quantity.get_values(location = 'centroids'))
2269
2270
2271        value = [[15,20,25]]
2272        quantity.set_values(value, indices = indices)
2273        #print "1 quantity.vertex_values",quantity.vertex_values
2274        assert allclose(quantity.vertex_values, quantity.get_values())
2275
2276        assert allclose(quantity.edge_values,
2277                        quantity.get_values(location = 'edges'))
2278
2279        # get a subset of elements
2280        subset = quantity.get_values(location='centroids', indices=[0,5])
2281        answer = [quantity.centroid_values[0],quantity.centroid_values[5]]
2282        assert allclose(subset, answer)
2283
2284
2285        subset = quantity.get_values(location='edges', indices=[0,5])
2286        answer = [quantity.edge_values[0],quantity.edge_values[5]]
2287        #print "subset",subset
2288        #print "answer",answer
2289        assert allclose(subset, answer)
2290
2291        subset = quantity.get_values( indices=[1,5])
2292        answer = [quantity.vertex_values[1],quantity.vertex_values[5]]
2293        #print "subset",subset
2294        #print "answer",answer
2295        assert allclose(subset, answer)
2296
2297    def test_smooth_vertex_values(self):
2298        """
2299        get values based on triangle lists.
2300        """
2301        from mesh_factory import rectangular
2302        from shallow_water import Domain
2303        from Numeric import zeros, Float
2304
2305        #Create basic mesh
2306        points, vertices, boundary = rectangular(2, 2)
2307
2308        #print "points",points
2309        #print "vertices",vertices
2310        #print "boundary",boundary
2311
2312        #Create shallow water domain
2313        domain = Domain(points, vertices, boundary)
2314        #print "domain.number_of_elements ",domain.number_of_elements
2315        quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
2316                                    [4,4,4],[5,5,5],[6,6,6],[7,7,7]])
2317
2318        #print "quantity.get_values(location = 'unique vertices')", \
2319        #      quantity.get_values(location = 'unique vertices')
2320
2321        #print "quantity.get_values(location = 'unique vertices')", \
2322        #      quantity.get_values(indices=[0,1,2,3,4,5,6,7], \
2323        #                          location = 'unique vertices')
2324
2325        #print quantity.get_values(location = 'unique vertices')
2326        #print quantity.domain.number_of_triangles_per_node
2327        #print quantity.vertex_values
2328       
2329        #answer = [0.5, 2, 3, 3, 3.5, 4, 4, 5, 6.5]
2330        #assert allclose(answer,
2331        #                quantity.get_values(location = 'unique vertices'))
2332
2333        quantity.smooth_vertex_values()
2334
2335        #print quantity.vertex_values
2336
2337
2338        answer_vertex_values = [[3,3.5,0.5],[2,0.5,3.5],[3.5,4,2],[3,2,4],
2339                                [4,5,3],[3.5,3,5],[5,6.5,3.5],[4,3.5,6.5]]
2340       
2341        assert allclose(answer_vertex_values,
2342                        quantity.vertex_values)
2343        #print "quantity.centroid_values",quantity.centroid_values
2344        #print "quantity.get_values(location = 'centroids') ",\
2345        #      quantity.get_values(location = 'centroids')
2346
2347
2348
2349#-------------------------------------------------------------
2350if __name__ == "__main__":
2351    suite = unittest.makeSuite(Test_Quantity, 'test')
2352
2353    #suite = unittest.makeSuite(Test_Quantity, 'test_set_values_from_UTM')
2354    #print "restricted test"
2355    #suite = unittest.makeSuite(Test_Quantity,'verbose_test_set_values_from_UTM_pts')
2356    runner = unittest.TextTestRunner() #verbosity=2)
2357    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.