Changeset 7793


Ignore:
Timestamp:
Jun 6, 2010, 10:48:41 PM (14 years ago)
Author:
steve
Message:

Setting up test function for Quantity

Location:
anuga_work/development/anuga_1d
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • anuga_work/development/anuga_1d/domain.py

    r6694 r7793  
    112112        self.build_tagged_elements_dictionary(tagged_elements)
    113113       
    114         from quantity import Quantity, Conserved_quantity
    115         #from quantity_domain import Quantity, Conserved_quantity
     114        from quantity import Quantity
     115
    116116       
    117117        #List of quantity names entering
     
    592592            return self.vertices[elem_id,vertex]
    593593
    594     def get_area(self, elem_id):
    595         """Return area of element id
    596         """
    597 
    598         return self.areas[elem_id]
     594    def get_area(self, elem_id=None):
     595        """Return total domain area or area of element id
     596        """
     597
     598        if elem_id is None:
     599            return sum(self.areas)
     600        else:
     601            return self.areas[elem_id]
    599602
    600603    def get_quantity(self, name, location='vertices', indices = None):
  • anuga_work/development/anuga_1d/quantity.py

    r6453 r7793  
    1 """Class Quantity - Implements values at each 1d element
     1"""
     2Class Quantity - Implements values at each 1d element
    23
    34To create:
     
    8586        """
    8687        return self.centroid_values.shape[0]
    87    
     88
     89    def __neg__(self):
     90        """Negate all values in this quantity giving meaning to the
     91        expression -Q where Q is an instance of class Quantity
     92        """
     93
     94        Q = Quantity(self.domain)
     95        Q.set_values_from_numeric(-self.vertex_values)
     96        return Q
     97
     98
     99
     100    def __add__(self, other):
     101        """Add to self anything that could populate a quantity
     102
     103        E.g other can be a constant, an array, a function, another quantity
     104        (except for a filename or points, attributes (for now))
     105        - see set_values_from_numeric for details
     106        """
     107
     108        Q = Quantity(self.domain)
     109        Q.set_values_from_numeric(other)
     110
     111        result = Quantity(self.domain)
     112        result.set_values_from_numeric(self.vertex_values + Q.vertex_values)
     113        return result
     114
     115    def __radd__(self, other):
     116        """Handle cases like 7+Q, where Q is an instance of class Quantity
     117        """
     118
     119        return self + other
     120
     121    def __sub__(self, other):
     122        return self + -other            # Invoke self.__neg__()
     123
     124    def __mul__(self, other):
     125        """Multiply self with anything that could populate a quantity
     126
     127        E.g other can be a constant, an array, a function, another quantity
     128        (except for a filename or points, attributes (for now))
     129        - see set_values_from_numeric for details
     130        """
     131
     132        if isinstance(other, Quantity):
     133            Q = other
     134        else:
     135            Q = Quantity(self.domain)
     136            Q.set_values_from_numeric(other)
     137
     138        result = Quantity(self.domain)
     139
     140        # The product of vertex_values, edge_values and centroid_values
     141        # are calculated and assigned directly without using
     142        # set_values_from_numeric (which calls interpolate). Otherwise
     143        # centroid values wouldn't be products from q1 and q2
     144        result.vertex_values = self.vertex_values * Q.vertex_values
     145        result.centroid_values = self.centroid_values * Q.centroid_values
     146
     147        return result
     148
     149    def __rmul__(self, other):
     150        """Handle cases like 3*Q, where Q is an instance of class Quantity
     151        """
     152
     153        return self * other
     154
     155    def __div__(self, other):
     156        """Divide self with anything that could populate a quantity
     157
     158        E.g other can be a constant, an array, a function, another quantity
     159        (except for a filename or points, attributes (for now))
     160        - see set_values_from_numeric for details
     161
     162        Zero division is dealt with by adding an epsilon to the divisore
     163        FIXME (Ole): Replace this with native INF once we migrate to NumPy
     164        """
     165
     166        if isinstance(other, Quantity):
     167            Q = other
     168        else:
     169            Q = Quantity(self.domain)
     170            Q.set_values_from_numeric(other)
     171
     172        result = Quantity(self.domain)
     173
     174        # The quotient of vertex_values, edge_values and centroid_values
     175        # are calculated and assigned directly without using
     176        # set_values_from_numeric (which calls interpolate). Otherwise
     177        # centroid values wouldn't be quotient of q1 and q2
     178        result.vertex_values = self.vertex_values/(Q.vertex_values + epsilon)
     179        result.centroid_values = self.centroid_values/(Q.centroid_values + epsilon)
     180
     181        return result
     182
     183    def __rdiv__(self, other):
     184        """Handle cases like 3/Q, where Q is an instance of class Quantity
     185        """
     186
     187        return self / other
     188
     189    def __pow__(self, other):
     190        """Raise quantity to (numerical) power
     191
     192        As with __mul__ vertex values are processed entry by entry
     193        while centroid and edge values are re-interpolated.
     194
     195        Example using __pow__:
     196          Q = (Q1**2 + Q2**2)**0.5
     197        """
     198
     199        if isinstance(other, Quantity):
     200            Q = other
     201        else:
     202            Q = Quantity(self.domain)
     203            Q.set_values_from_numeric(other)
     204
     205        result = Quantity(self.domain)
     206
     207        # The power of vertex_values, edge_values and centroid_values
     208        # are calculated and assigned directly without using
     209        # set_values_from_numeric (which calls interpolate). Otherwise
     210        # centroid values wouldn't be correct
     211        result.vertex_values = self.vertex_values ** other
     212        result.centroid_values = self.centroid_values ** other
     213
     214        return result
     215
     216    def set_values_from_numeric(self, numeric):
     217        import Numeric
     218        x = Numeric.array([1.0,2.0])
     219        y = [1.0,2.0]
     220
     221        if type(numeric) == type(y):
     222            self.set_values_from_array(numeric)
     223        elif type(numeric) == type(x):
     224            self.set_values_from_array(numeric)
     225        elif callable(numeric):
     226            self.set_values_from_function(numeric)
     227        elif isinstance(numeric, Quantity):
     228            self.set_values_from_quantity(numeric)
     229        else:   # see if it's coercible to a float (float, int or long, etc)
     230            try:
     231                numeric = float(numeric)
     232            except ValueError:
     233                msg = ("Illegal type for variable 'numeric': %s" % type(numeric))
     234                raise Exception, msg
     235            self.set_values_from_constant(numeric)
     236
     237    def set_values_from_constant(self,numeric):
     238
     239        self.vertex_values[:]   = numeric
     240        self.centroid_values[:] = numeric
     241
     242
     243    def set_values_from_array(self,numeric):
     244
     245        self.vertex_values[:]   = numeric
     246        self.interpolate()
     247
     248
     249    def set_values_from_quantity(self,quantity):
     250
     251        self.vertex_values[:]   = quantity.vertex_values
     252        self.centroid_values[:] = quantity.centroid_values
     253
     254    def set_values_from_function(self,function):
     255
     256        self.vertex_values[:]   = map(function, self.domain.vertices)
     257        self.centroid_values[:] = map(function, self.domain.centroids)
     258
     259
    88260    def interpolate(self):
    89261        """
     
    218390            assert values.shape[1] == 2, msg
    219391
    220             self.vertex_values = values
     392            self.vertex_values[:] = values
    221393
    222394
     
    428600        """
    429601
    430         #print 'compute_gradient'
    431602
    432603        from Numeric import array, zeros, Float
     
    9081079        self.centroid_values[:] = (a*self.centroid_values + b*self.centroid_backup_values).astype('f')
    9091080       
    910 class Conserved_quantity(Quantity):
    911     """Class conserved quantity adds to Quantity:
    912 
    913     storage and method for updating, and
    914     methods for extrapolation from centropid to vertices inluding
    915     gradients and limiters
    916     """
    917 
    918     def __init__(self, domain, vertex_values=None):
    919         Quantity.__init__(self, domain, vertex_values)
    920 
    921         print "Use Quantity instead of Conserved_quantity"
    922 
    923 
    924 ##
    925 ##def newLinePlot(title='Simple Plot'):
    926 ##    import Gnuplot
    927 ##    g = Gnuplot.Gnuplot()
    928 ##    g.title(title)
    929 ##    g('set data style linespoints')
    930 ##    g.xlabel('x')
    931 ##    g.ylabel('y')
    932 ##    return g
    933 ##
    934 ##def linePlot(g,x,y):
    935 ##    import Gnuplot
    936 ##    g.plot(Gnuplot.PlotItems.Data(x.flat,y.flat))
     1081
    9371082
    9381083def newLinePlot(title='Simple Plot'):
  • anuga_work/development/anuga_1d/test_quantity.py

    r5858 r7793  
     1import quantity
    12#!/usr/bin/env python
    23
    34import unittest
    45from math import sqrt, pi
    5 import tempfile
     6
    67
    78from quantity import *
    8 from config import epsilon
    9 from Numeric import allclose, array, ones, Float
    10        
    11 from domain import Domain
    12 
    13 
    14 #Aux for fit_interpolate.fit example
    15 def linear_function(point):
    16     point = array(point)
    17     return point[:,0]+point[:,1]
     9from domain import *
     10
     11
     12from Numeric import allclose, array, ones, Float, maximum, zeros
    1813
    1914
    2015class Test_Quantity(unittest.TestCase):
    2116    def setUp(self):
    22         from domain import Domain
    23 
    24         a = [0.0, 0.0]
    25         b = [0.0, 2.0]
    26         c = [2.0, 0.0]
    27         d = [0.0, 4.0]
    28         e = [2.0, 2.0]
    29         f = [4.0, 0.0]
    30 
    31         points = [a, b, c, d, e, f]
    32 
    33         #bac, bce, ecf, dbe
    34         elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
    35 
    36         self.mesh1 = Domain(points[:3], [elements[0]])
    37         self.mesh1.check_integrity()
    38 
    39         self.mesh4 = Domain(points, elements)
    40         self.mesh4.check_integrity()
    41 
    42         # UTM round Onslow
    43         a = [240000, 7620000]
    44         b = [240000, 7680000]
    45         c = [300000, 7620000]
    46 
    47         points = [a, b, c]
    48         elements = [[0,2,1]]
    49        
    50         self.mesh_onslow = Domain(points, elements)
    51         self.mesh_onslow.check_integrity()
    52        
     17        self.points3        = [0.0, 1.0, 2.0, 3.0]
     18        self.vertex_values3 = [[1.0,2.0],[4.0,5.0],[-1.0,2.0]]
     19        self.domain3        = Domain(self.points3)
     20
     21
     22       
     23        self.points4          = [0.0, 1.0, 2.5, 3.0, 5.0]
     24        self.vertex_values4   = [[1.0,2.0],[4.0,5.0],[-1.0,2.0],[3.0,6.0]]
     25        self.centroid_values4 = [1.5, 4.5, 0.5, 4.5]
     26        self.boundary4        = {(0, 0): 'left', (3, 1): 'right'}
     27        self.domain4          = Domain(self.points4,self.boundary4)
     28
     29        self.points10 = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]
     30        self.domain10 = Domain(self.points10)
     31
     32        self.points6 = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0]
     33        self.domain6 = Domain(self.points6)
     34
     35
    5336    def tearDown(self):
    5437        pass
     
    5639
    5740
     41    def test_creat_with_boundary(self):
     42
     43        assert self.domain4.boundary  == {(0, 0): 'left', (3, 1): 'right'}
     44       
    5845    def test_creation(self):
    5946
    60         quantity = Quantity(self.mesh1, [[1,2,3]])
    61         assert allclose(quantity.vertex_values, [[1.,2.,3.]])
    62 
     47        quantity = Quantity(self.domain3)
     48        assert allclose(quantity.vertex_values, [[0.0,0.0],[0.0,0.0],[0.0,0.0]])
     49
     50       
    6351        try:
    6452            quantity = Quantity()
     
    7462            pass
    7563        except:
    76             raise 'Should have raised "mising mesh object" error'
     64            raise 'Should have raised "mising domain object" error'
    7765
    7866
    7967    def test_creation_zeros(self):
    8068
    81         quantity = Quantity(self.mesh1)
    82         assert allclose(quantity.vertex_values, [[0.,0.,0.]])
    83 
    84 
    85         quantity = Quantity(self.mesh4)
    86         assert allclose(quantity.vertex_values, [[0.,0.,0.], [0.,0.,0.],
    87                                                  [0.,0.,0.], [0.,0.,0.]])
     69        quantity = Quantity(self.domain3)
     70        assert allclose(quantity.centroid_values, [[0.,0.,0.]])
     71
     72
     73        quantity = Quantity(self.domain4)
     74        assert allclose(quantity.vertex_values, [[0.,0.], [0.,0.],
     75                                                 [0.,0.], [0.,0.]])
    8876
    8977
    9078    def test_interpolation(self):
    91         quantity = Quantity(self.mesh1, [[1,2,3]])
    92         assert allclose(quantity.centroid_values, [2.0]) #Centroid
    93 
    94         assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5]])
     79        quantity = Quantity(self.domain4, self.vertex_values4)
     80        assert allclose(quantity.centroid_values, self.centroid_values4) #Centroid
     81
    9582
    9683
    9784    def test_interpolation2(self):
    98         quantity = Quantity(self.mesh4,
    99                             [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
    100         assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
    101 
     85        quantity = Quantity(self.domain4, self.vertex_values4)
     86        assert allclose(quantity.centroid_values, self.centroid_values4) #Centroid
    10287
    10388        quantity.extrapolate_second_order()
    10489
    10590        #print quantity.vertex_values
    106         assert allclose(quantity.vertex_values, [[3.5, -1.0, 3.5],
    107                                                  [3.+2./3, 6.+2./3, 4.+2./3],
    108                                                  [4.6, 3.4, 1.],
    109                                                  [-5.0, 1.0, 4.0]])
    110 
    111         #print quantity.edge_values
    112         assert allclose(quantity.edge_values, [[1.25, 3.5, 1.25],
    113                                                [5. + 2/3.0, 4.0 + 1.0/6, 5.0 + 1.0/6],
    114                                                [2.2, 2.8, 4.0],
    115                                                [2.5, -0.5, -2.0]])
    116 
    117 
    118         #assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
    119         #                                       [5., 5., 5.],
    120         #                                       [4.5, 4.5, 0.],
    121         #                                       [3.0, -1.5, -1.5]])
    122 
    123     def test_get_extrema_1(self):
    124         quantity = Quantity(self.mesh4,
    125                                       [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
    126         assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroids
    127 
    128         v = quantity.get_maximum_value()
    129         assert v == 5
    130 
    131         v = quantity.get_minimum_value()
    132         assert v == 0       
    133 
    134         i = quantity.get_maximum_index()
    135         assert i == 1
    136 
    137         i = quantity.get_minimum_index()
    138         assert i == 3       
    139        
    140         x,y = quantity.get_maximum_location()
    141         xref, yref = 4.0/3, 4.0/3
    142         assert x == xref
    143         assert y == yref
    144 
    145         v = quantity.get_values(interpolation_points = [[x,y]])
    146         assert allclose(v, 5)
    147 
    148 
    149         x,y = quantity.get_minimum_location()
    150         v = quantity.get_values(interpolation_points = [[x,y]])
    151         assert allclose(v, 0)
    152 
    153 
    154     def test_get_maximum_2(self):
    155 
    156         a = [0.0, 0.0]
    157         b = [0.0, 2.0]
    158         c = [2.0,0.0]
    159         d = [0.0, 4.0]
    160         e = [2.0, 2.0]
    161         f = [4.0,0.0]
    162 
    163         points = [a, b, c, d, e, f]
    164         #bac, bce, ecf, dbe
    165         vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    166 
    167         domain = Domain(points, vertices)
    168 
    169         quantity = Quantity(domain)
    170         quantity.set_values(lambda x, y: x+2*y) #2 4 4 6
    171        
    172         v = quantity.get_maximum_value()
    173         assert v == 6
    174 
    175         v = quantity.get_minimum_value()
    176         assert v == 2       
    177 
    178         i = quantity.get_maximum_index()
    179         assert i == 3
    180 
    181         i = quantity.get_minimum_index()
    182         assert i == 0       
    183        
    184         x,y = quantity.get_maximum_location()
    185         xref, yref = 2.0/3, 8.0/3
    186         assert x == xref
    187         assert y == yref
    188 
    189         v = quantity.get_values(interpolation_points = [[x,y]])
    190         assert allclose(v, 6)
    191 
    192         x,y = quantity.get_minimum_location()       
    193         v = quantity.get_values(interpolation_points = [[x,y]])
    194         assert allclose(v, 2)
    195 
    196         #Multiple locations for maximum -
    197         #Test that the algorithm picks the first occurrence       
    198         v = quantity.get_maximum_value(indices=[0,1,2])
    199         assert allclose(v, 4)
    200 
    201         i = quantity.get_maximum_index(indices=[0,1,2])
    202         assert i == 1
    203        
    204         x,y = quantity.get_maximum_location(indices=[0,1,2])
    205         xref, yref = 4.0/3, 4.0/3
    206         assert x == xref
    207         assert y == yref
    208 
    209         v = quantity.get_values(interpolation_points = [[x,y]])
    210         assert allclose(v, 4)       
    211 
    212         # More test of indices......
    213         v = quantity.get_maximum_value(indices=[2,3])
    214         assert allclose(v, 6)
    215 
    216         i = quantity.get_maximum_index(indices=[2,3])
    217         assert i == 3
    218        
    219         x,y = quantity.get_maximum_location(indices=[2,3])
    220         xref, yref = 2.0/3, 8.0/3
    221         assert x == xref
    222         assert y == yref
    223 
    224         v = quantity.get_values(interpolation_points = [[x,y]])
    225         assert allclose(v, 6)       
    226 
    227        
    228 
     91        assert allclose(quantity.vertex_values,[[ 0.3,  2.7],
     92                                                [ 4.5,  4.5],
     93                                                [ 0.5,  0.5],
     94                                                [ 1.3,  7.7]])
     95
     96 
     97
     98 
    22999    def test_boundary_allocation(self):
    230         quantity = Quantity(self.mesh4,
    231                             [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
    232 
    233         assert quantity.boundary_values.shape[0] == len(self.mesh4.boundary)
     100        quantity = Quantity(self.domain4,
     101                            [[1,2], [5,5], [0,9], [-6, 3]])
     102
     103
     104        assert quantity.boundary_values.shape[0] == len(self.domain4.boundary)
    234105
    235106
    236107    def test_set_values(self):
    237         quantity = Quantity(self.mesh4)
     108        quantity = Quantity(self.domain4)
    238109
    239110
     
    281152
    282153    def test_set_values_const(self):
    283         quantity = Quantity(self.mesh4)
     154        quantity = Quantity(self.domain4)
    284155
    285156        quantity.set_values(1.0, location = 'vertices')
     
    299170
    300171    def test_set_values_func(self):
    301         quantity = Quantity(self.mesh4)
     172        quantity = Quantity(self.domain4)
    302173
    303174        def f(x, y):
     
    320191
    321192    def test_integral(self):
    322         quantity = Quantity(self.mesh4)
     193        quantity = Quantity(self.domain4)
    323194
    324195        #Try constants first
     
    327198        #print 'Q', quantity.get_integral()
    328199
    329         assert allclose(quantity.get_integral(), self.mesh4.get_area() * const)
     200        assert allclose(quantity.get_integral(), self.domain4.get_area() * const)
    330201
    331202        #Try with a linear function
    332         def f(x, y):
    333             return x+y
     203        def f(x):
     204            return x
    334205
    335206        quantity.set_values(f, location = 'vertices')
    336207
    337 
    338         ref_integral = (4.0/3 + 8.0/3 + 10.0/3 + 10.0/3) * 2
     208 
     209        assert allclose (quantity.centroid_values,
     210        [ 0.5,   1.75,  2.75,  4.  ])
     211
     212        assert allclose (quantity.vertex_values, [[ 0.,   1. ],
     213                                                  [ 1.,   2.5],
     214                                                  [ 2.5,  3. ],
     215                                                  [ 3.,   5. ]])
     216
     217
     218        ref_integral = 0.5 + 1.5*1.75 + 0.5*2.75 + 2.0*4.0
    339219
    340220        assert allclose (quantity.get_integral(), ref_integral)
     
    343223
    344224    def test_set_vertex_values(self):
    345         quantity = Quantity(self.mesh4)
     225        quantity = Quantity(self.domain4)
    346226        quantity.set_vertex_values([0,1,2,3,4,5])
    347227
     
    357237
    358238    def test_set_vertex_values_subset(self):
    359         quantity = Quantity(self.mesh4)
     239        quantity = Quantity(self.domain4)
    360240        quantity.set_vertex_values([0,1,2,3,4,5])
    361241        quantity.set_vertex_values([0,20,30,50], indices = [0,2,3,5])
     
    366246
    367247    def test_set_vertex_values_using_general_interface(self):
    368         quantity = Quantity(self.mesh4)
     248        quantity = Quantity(self.domain4)
    369249
    370250
     
    391271        """
    392272       
    393         quantity = Quantity(self.mesh4)
     273        quantity = Quantity(self.domain4)
    394274
    395275
     
    462342        """
    463343       
    464         quantity = Quantity(self.mesh4)
     344        quantity = Quantity(self.domain4)
    465345        G = Geo_reference(56, 10, 100)
    466346        quantity.domain.geo_reference = G
     
    507387
    508388
    509     def test_set_values_using_fit(self):
    510 
    511 
    512         quantity = Quantity(self.mesh4)
    513 
    514         #Get (enough) datapoints
    515         data_points = [[ 0.66666667, 0.66666667],
    516                        [ 1.33333333, 1.33333333],
    517                        [ 2.66666667, 0.66666667],
    518                        [ 0.66666667, 2.66666667],
    519                        [ 0.0, 1.0],
    520                        [ 0.0, 3.0],
    521                        [ 1.0, 0.0],
    522                        [ 1.0, 1.0],
    523                        [ 1.0, 2.0],
    524                        [ 1.0, 3.0],
    525                        [ 2.0, 1.0],
    526                        [ 3.0, 0.0],
    527                        [ 3.0, 1.0]]
    528 
    529         z = linear_function(data_points)
    530 
    531         #Use built-in fit_interpolate.fit
    532         quantity.set_values( Geospatial_data(data_points, z), alpha = 0 )
    533         #quantity.set_values(points = data_points, values = z, alpha = 0)
    534 
    535 
    536         answer = linear_function(quantity.domain.get_vertex_coordinates())
    537         #print quantity.vertex_values, answer
    538         assert allclose(quantity.vertex_values.flat, answer)
    539 
    540 
    541         #Now try by setting the same values directly
    542         vertex_attributes = fit_to_mesh(data_points,
    543                                         quantity.domain.get_nodes(),
    544                                         quantity.domain.triangles, #FIXME
    545                                         point_attributes=z,
    546                                         alpha = 0,
    547                                         verbose=False)
    548 
    549         #print vertex_attributes
    550         quantity.set_values(vertex_attributes)
    551         assert allclose(quantity.vertex_values.flat, answer)
    552 
    553 
    554 
    555 
    556 
    557     def test_test_set_values_using_fit_w_geo(self):
    558 
    559 
    560         #Mesh
    561         vertex_coordinates = [[0.76, 0.76],
    562                               [0.76, 5.76],
    563                               [5.76, 0.76]]
    564         triangles = [[0,2,1]]
    565 
    566         mesh_georef = Geo_reference(56,-0.76,-0.76)
    567         mesh1 = Domain(vertex_coordinates, triangles,
    568                        geo_reference = mesh_georef)
    569         mesh1.check_integrity()
    570 
    571         #Quantity
    572         quantity = Quantity(mesh1)
    573 
    574         #Data
    575         data_points = [[ 201.0, 401.0],
    576                        [ 201.0, 403.0],
    577                        [ 203.0, 401.0]]
    578 
    579         z = [2, 4, 4]
    580 
    581         data_georef = Geo_reference(56,-200,-400)
    582 
    583 
    584         #Reference
    585         ref = fit_to_mesh(data_points, vertex_coordinates, triangles,
    586                           point_attributes=z,
    587                           data_origin = data_georef.get_origin(),
    588                           mesh_origin = mesh_georef.get_origin(),
    589                           alpha = 0)
    590 
    591         assert allclose( ref, [0,5,5] )
    592 
    593 
    594         #Test set_values
    595 
    596         quantity.set_values( Geospatial_data(data_points, z, data_georef), alpha = 0 )
    597 
    598         #quantity.set_values(points = data_points,
    599         #                    values = z,
    600         #                    data_georef = data_georef,
    601         #                    alpha = 0)
    602 
    603 
    604         #quantity.set_values(points = data_points,
    605         #                    values = z,
    606         #                    data_georef = data_georef,
    607         #                    alpha = 0)
    608         assert allclose(quantity.vertex_values.flat, ref)
    609 
    610 
    611 
    612         #Test set_values using geospatial data object
    613         quantity.vertex_values[:] = 0.0
    614 
    615         geo = Geospatial_data(data_points, z, data_georef)
    616 
    617 
    618         quantity.set_values(geospatial_data = geo, alpha = 0)
    619         assert allclose(quantity.vertex_values.flat, ref)
    620 
    621 
    622 
    623     def test_set_values_from_file1(self):
    624         quantity = Quantity(self.mesh4)
    625 
    626         #Get (enough) datapoints
    627         data_points = [[ 0.66666667, 0.66666667],
    628                        [ 1.33333333, 1.33333333],
    629                        [ 2.66666667, 0.66666667],
    630                        [ 0.66666667, 2.66666667],
    631                        [ 0.0, 1.0],
    632                        [ 0.0, 3.0],
    633                        [ 1.0, 0.0],
    634                        [ 1.0, 1.0],
    635                        [ 1.0, 2.0],
    636                        [ 1.0, 3.0],
    637                        [ 2.0, 1.0],
    638                        [ 3.0, 0.0],
    639                        [ 3.0, 1.0]]
    640 
    641         data_geo_spatial = Geospatial_data(data_points,
    642                          geo_reference = Geo_reference(56, 0, 0))
    643         data_points_absolute = data_geo_spatial.get_data_points(absolute=True)
    644         attributes = linear_function(data_points_absolute)
    645         att = 'spam_and_eggs'
    646        
    647         #Create .txt file
    648         ptsfile = tempfile.mktemp(".txt")
    649         file = open(ptsfile,"w")
    650         file.write(" x,y," + att + " \n")
    651         for data_point, attribute in map(None, data_points_absolute
    652                                          ,attributes):
    653             row = str(data_point[0]) + ',' + str(data_point[1]) \
    654                   + ',' + str(attribute)
    655             file.write(row + "\n")
    656         file.close()
    657 
    658 
    659         #Check that values can be set from file
    660         quantity.set_values(filename = ptsfile,
    661                             attribute_name = att, alpha = 0)
    662         answer = linear_function(quantity.domain.get_vertex_coordinates())
    663 
    664         #print quantity.vertex_values.flat
    665         #print answer
    666 
    667 
    668         assert allclose(quantity.vertex_values.flat, answer)
    669 
    670 
    671         #Check that values can be set from file using default attribute
    672         quantity.set_values(filename = ptsfile, alpha = 0)
    673         assert allclose(quantity.vertex_values.flat, answer)
    674 
    675         #Cleanup
    676         import os
    677         os.remove(ptsfile)
    678 
    679 
    680 
    681     def Xtest_set_values_from_file_using_polygon(self):
    682         """test_set_values_from_file_using_polygon(self):
    683        
    684         Test that polygon restriction works for general points data
    685         """
    686        
    687         quantity = Quantity(self.mesh4)
    688 
    689         #Get (enough) datapoints
    690         data_points = [[ 0.66666667, 0.66666667],
    691                        [ 1.33333333, 1.33333333],
    692                        [ 2.66666667, 0.66666667],
    693                        [ 0.66666667, 2.66666667],
    694                        [ 0.0, 1.0],
    695                        [ 0.0, 3.0],
    696                        [ 1.0, 0.0],
    697                        [ 1.0, 1.0],
    698                        [ 1.0, 2.0],
    699                        [ 1.0, 3.0],
    700                        [ 2.0, 1.0],
    701                        [ 3.0, 0.0],
    702                        [ 3.0, 1.0]]
    703 
    704         data_geo_spatial = Geospatial_data(data_points,
    705                          geo_reference = Geo_reference(56, 0, 0))
    706         data_points_absolute = data_geo_spatial.get_data_points(absolute=True)
    707         attributes = linear_function(data_points_absolute)
    708         att = 'spam_and_eggs'
    709        
    710         #Create .txt file
    711         ptsfile = tempfile.mktemp(".txt")
    712         file = open(ptsfile,"w")
    713         file.write(" x,y," + att + " \n")
    714         for data_point, attribute in map(None, data_points_absolute
    715                                          ,attributes):
    716             row = str(data_point[0]) + ',' + str(data_point[1]) \
    717                   + ',' + str(attribute)
    718             file.write(row + "\n")
    719         file.close()
    720 
    721         # Create restricting polygon (containing node #4 (2,2) and
    722         # centroid of triangle #1 (bce)
    723         polygon = [[1.0, 1.0], [4.0, 1.0],
    724                    [4.0, 4.0], [1.0, 4.0]]
    725 
    726         #print self.mesh4.nodes
    727         #print inside_polygon(self.mesh4.nodes, polygon)
    728         assert allclose(inside_polygon(self.mesh4.nodes, polygon), 4)
    729 
    730         #print quantity.domain.get_vertex_coordinates()
    731         #print quantity.domain.get_nodes()       
    732        
    733         # Check that values can be set from file
    734         quantity.set_values(filename=ptsfile,
    735                             polygon=polygon,
    736                             location='unique vertices',
    737                             alpha=0)
    738 
    739         # Get indices for vertex coordinates in polygon
    740         indices = inside_polygon(quantity.domain.get_vertex_coordinates(),
    741                                  polygon)
    742         points = take(quantity.domain.get_vertex_coordinates(), indices)
    743        
    744         answer = linear_function(points)
    745 
    746         #print quantity.vertex_values.flat
    747         #print answer
    748 
    749         # Check vertices in polygon have been set
    750         assert allclose(take(quantity.vertex_values.flat, indices),
    751                              answer)
    752 
    753         # Check vertices outside polygon are zero
    754         indices = outside_polygon(quantity.domain.get_vertex_coordinates(),
    755                                   polygon)       
    756         assert allclose(take(quantity.vertex_values.flat, indices),
    757                              0.0)       
    758 
    759         #Cleanup
    760         import os
    761         os.remove(ptsfile)
    762 
    763 
    764        
    765 
    766     def test_cache_test_set_values_from_file(self):
    767         # FIXME (Ole): What is this about?
    768         # I don't think it checks anything new
    769         quantity = Quantity(self.mesh4)
    770 
    771         #Get (enough) datapoints
    772         data_points = [[ 0.66666667, 0.66666667],
    773                        [ 1.33333333, 1.33333333],
    774                        [ 2.66666667, 0.66666667],
    775                        [ 0.66666667, 2.66666667],
    776                        [ 0.0, 1.0],
    777                        [ 0.0, 3.0],
    778                        [ 1.0, 0.0],
    779                        [ 1.0, 1.0],
    780                        [ 1.0, 2.0],
    781                        [ 1.0, 3.0],
    782                        [ 2.0, 1.0],
    783                        [ 3.0, 0.0],
    784                        [ 3.0, 1.0]]
    785 
    786         georef = Geo_reference(56, 0, 0)
    787         data_geo_spatial = Geospatial_data(data_points,
    788                                            geo_reference=georef)
    789                                            
    790         data_points_absolute = data_geo_spatial.get_data_points(absolute=True)
    791         attributes = linear_function(data_points_absolute)
    792         att = 'spam_and_eggs'
    793        
    794         # Create .txt file
    795         ptsfile = tempfile.mktemp(".txt")
    796         file = open(ptsfile,"w")
    797         file.write(" x,y," + att + " \n")
    798         for data_point, attribute in map(None, data_points_absolute
    799                                          ,attributes):
    800             row = str(data_point[0]) + ',' + str(data_point[1]) \
    801                   + ',' + str(attribute)
    802             file.write(row + "\n")
    803         file.close()
    804 
    805 
    806         # Check that values can be set from file
    807         quantity.set_values(filename=ptsfile,
    808                             attribute_name=att,
    809                             alpha=0,
    810                             use_cache=True,
    811                             verbose=False)
    812         answer = linear_function(quantity.domain.get_vertex_coordinates())
    813         assert allclose(quantity.vertex_values.flat, answer)
    814 
    815 
    816         # Check that values can be set from file using default attribute
    817         quantity.set_values(filename=ptsfile,
    818                             alpha=0)
    819         assert allclose(quantity.vertex_values.flat, answer)
    820 
    821         # Check cache
    822         quantity.set_values(filename=ptsfile,
    823                             attribute_name=att,
    824                             alpha=0,
    825                             use_cache=True,
    826                             verbose=False)
    827        
    828        
    829         #Cleanup
    830         import os
    831         os.remove(ptsfile)
    832 
    833     def test_set_values_from_lat_long(self):
    834         quantity = Quantity(self.mesh_onslow)
    835 
    836         #Get (enough) datapoints
    837         data_points = [[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
    838                        [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6]]
    839 
    840         data_geo_spatial = Geospatial_data(data_points,
    841                                            points_are_lats_longs=True)
    842         points_UTM = data_geo_spatial.get_data_points(absolute=True)
    843         attributes = linear_function(points_UTM)
    844         att = 'elevation'
    845        
    846         #Create .txt file
    847         txt_file = tempfile.mktemp(".txt")
    848         file = open(txt_file,"w")
    849         file.write(" lat,long," + att + " \n")
    850         for data_point, attribute in map(None, data_points, attributes):
    851             row = str(data_point[0]) + ',' + str(data_point[1]) \
    852                   + ',' + str(attribute)
    853             #print "row", row
    854             file.write(row + "\n")
    855         file.close()
    856 
    857 
    858         #Check that values can be set from file
    859         quantity.set_values(filename=txt_file,
    860                             attribute_name=att,
    861                             alpha=0)
    862         answer = linear_function(quantity.domain.get_vertex_coordinates())
    863 
    864         #print "quantity.vertex_values.flat", quantity.vertex_values.flat
    865         #print "answer",answer
    866 
    867         assert allclose(quantity.vertex_values.flat, answer)
    868 
    869 
    870         #Check that values can be set from file using default attribute
    871         quantity.set_values(filename=txt_file, alpha=0)
    872         assert allclose(quantity.vertex_values.flat, answer)
    873 
    874         #Cleanup
    875         import os
    876         os.remove(txt_file)
    877          
    878     def test_set_values_from_lat_long(self):
    879         quantity = Quantity(self.mesh_onslow)
    880 
    881         #Get (enough) datapoints
    882         data_points = [[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
    883                        [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6]]
    884 
    885         data_geo_spatial = Geospatial_data(data_points,
    886                                            points_are_lats_longs=True)
    887         points_UTM = data_geo_spatial.get_data_points(absolute=True)
    888         attributes = linear_function(points_UTM)
    889         att = 'elevation'
    890        
    891         #Create .txt file
    892         txt_file = tempfile.mktemp(".txt")
    893         file = open(txt_file,"w")
    894         file.write(" lat,long," + att + " \n")
    895         for data_point, attribute in map(None, data_points, attributes):
    896             row = str(data_point[0]) + ',' + str(data_point[1]) \
    897                   + ',' + str(attribute)
    898             #print "row", row
    899             file.write(row + "\n")
    900         file.close()
    901 
    902 
    903         #Check that values can be set from file
    904         quantity.set_values(filename=txt_file,
    905                             attribute_name=att, alpha=0)
    906         answer = linear_function(quantity.domain.get_vertex_coordinates())
    907 
    908         #print "quantity.vertex_values.flat", quantity.vertex_values.flat
    909         #print "answer",answer
    910 
    911         assert allclose(quantity.vertex_values.flat, answer)
    912 
    913 
    914         #Check that values can be set from file using default attribute
    915         quantity.set_values(filename=txt_file, alpha=0)
    916         assert allclose(quantity.vertex_values.flat, answer)
    917 
    918         #Cleanup
    919         import os
    920         os.remove(txt_file)
    921        
    922     def test_set_values_from_UTM_pts(self):
    923         quantity = Quantity(self.mesh_onslow)
    924 
    925         #Get (enough) datapoints
    926         data_points = [[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
    927                        [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6]]
    928 
    929         data_geo_spatial = Geospatial_data(data_points,
    930                                            points_are_lats_longs=True)
    931         points_UTM = data_geo_spatial.get_data_points(absolute=True)
    932         attributes = linear_function(points_UTM)
    933         att = 'elevation'
    934        
    935         #Create .txt file
    936         txt_file = tempfile.mktemp(".txt")
    937         file = open(txt_file,"w")
    938         file.write(" x,y," + att + " \n")
    939         for data_point, attribute in map(None, points_UTM, attributes):
    940             row = str(data_point[0]) + ',' + str(data_point[1]) \
    941                   + ',' + str(attribute)
    942             #print "row", row
    943             file.write(row + "\n")
    944         file.close()
    945 
    946 
    947         pts_file = tempfile.mktemp(".pts")       
    948         convert = Geospatial_data(txt_file)
    949         convert.export_points_file(pts_file)
    950 
    951         #Check that values can be set from file
    952         quantity.set_values_from_file(pts_file, att, 0,
    953                                       'vertices', None)
    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         #Check that values can be set from file
    960         quantity.set_values(filename=pts_file,
    961                             attribute_name=att, alpha=0)
    962         answer = linear_function(quantity.domain.get_vertex_coordinates())
    963         #print "quantity.vertex_values.flat", quantity.vertex_values.flat
    964         #print "answer",answer
    965         assert allclose(quantity.vertex_values.flat, answer)
    966 
    967 
    968         #Check that values can be set from file using default attribute
    969         quantity.set_values(filename=txt_file, alpha=0)
    970         assert allclose(quantity.vertex_values.flat, answer)
    971 
    972         #Cleanup
    973         import os
    974         os.remove(txt_file)
    975         os.remove(pts_file)
    976        
    977     def verbose_test_set_values_from_UTM_pts(self):
    978         quantity = Quantity(self.mesh_onslow)
    979 
    980         #Get (enough) datapoints
    981         data_points = [[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
    982                        [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
    983                        [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
    984                        [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
    985                        [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
    986                        [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
    987                        [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
    988                        [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
    989                        [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
    990                        [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
    991                        [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
    992                        [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
    993                        [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
    994                        [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
    995                        [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
    996                        [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
    997                        [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
    998                        [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
    999                        [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
    1000                        [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
    1001                        [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
    1002                        [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
    1003                        [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
    1004                        [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
    1005                        [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
    1006                        [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
    1007                        [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
    1008                        ]
    1009 
    1010         data_geo_spatial = Geospatial_data(data_points,
    1011                                            points_are_lats_longs=True)
    1012         points_UTM = data_geo_spatial.get_data_points(absolute=True)
    1013         attributes = linear_function(points_UTM)
    1014         att = 'elevation'
    1015        
    1016         #Create .txt file
    1017         txt_file = tempfile.mktemp(".txt")
    1018         file = open(txt_file,"w")
    1019         file.write(" x,y," + att + " \n")
    1020         for data_point, attribute in map(None, points_UTM, attributes):
    1021             row = str(data_point[0]) + ',' + str(data_point[1]) \
    1022                   + ',' + str(attribute)
    1023             #print "row", row
    1024             file.write(row + "\n")
    1025         file.close()
    1026 
    1027 
    1028         pts_file = tempfile.mktemp(".pts")       
    1029         convert = Geospatial_data(txt_file)
    1030         convert.export_points_file(pts_file)
    1031 
    1032         #Check that values can be set from file
    1033         quantity.set_values_from_file(pts_file, att, 0,
    1034                                       'vertices', None, verbose = True,
    1035                                       max_read_lines=2)
    1036         answer = linear_function(quantity.domain.get_vertex_coordinates())
    1037         #print "quantity.vertex_values.flat", quantity.vertex_values.flat
    1038         #print "answer",answer
    1039         assert allclose(quantity.vertex_values.flat, answer)
    1040 
    1041         #Check that values can be set from file
    1042         quantity.set_values(filename=pts_file,
    1043                             attribute_name=att, alpha=0)
    1044         answer = linear_function(quantity.domain.get_vertex_coordinates())
    1045         #print "quantity.vertex_values.flat", quantity.vertex_values.flat
    1046         #print "answer",answer
    1047         assert allclose(quantity.vertex_values.flat, answer)
    1048 
    1049 
    1050         #Check that values can be set from file using default attribute
    1051         quantity.set_values(filename=txt_file, alpha=0)
    1052         assert allclose(quantity.vertex_values.flat, answer)
    1053 
    1054         #Cleanup
    1055         import os
    1056         os.remove(txt_file)
    1057         os.remove(pts_file)
    1058        
    1059     def test_set_values_from_file_with_georef1(self):
    1060 
    1061         #Mesh in zone 56 (absolute coords)
    1062 
    1063         x0 = 314036.58727982
    1064         y0 = 6224951.2960092
    1065 
    1066         a = [x0+0.0, y0+0.0]
    1067         b = [x0+0.0, y0+2.0]
    1068         c = [x0+2.0, y0+0.0]
    1069         d = [x0+0.0, y0+4.0]
    1070         e = [x0+2.0, y0+2.0]
    1071         f = [x0+4.0, y0+0.0]
    1072 
    1073         points = [a, b, c, d, e, f]
    1074 
    1075         #bac, bce, ecf, dbe
    1076         elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
    1077 
    1078         #absolute going in ..
    1079         mesh4 = Domain(points, elements,
    1080                        geo_reference = Geo_reference(56, 0, 0))
    1081         mesh4.check_integrity()
    1082         quantity = Quantity(mesh4)
    1083 
    1084         #Get (enough) datapoints (relative to georef)
    1085         data_points_rel = [[ 0.66666667, 0.66666667],
    1086                        [ 1.33333333, 1.33333333],
    1087                        [ 2.66666667, 0.66666667],
    1088                        [ 0.66666667, 2.66666667],
    1089                        [ 0.0, 1.0],
    1090                        [ 0.0, 3.0],
    1091                        [ 1.0, 0.0],
    1092                        [ 1.0, 1.0],
    1093                        [ 1.0, 2.0],
    1094                        [ 1.0, 3.0],
    1095                        [ 2.0, 1.0],
    1096                        [ 3.0, 0.0],
    1097                        [ 3.0, 1.0]]
    1098 
    1099         data_geo_spatial = Geospatial_data(data_points_rel,
    1100                          geo_reference = Geo_reference(56, x0, y0))
    1101         data_points_absolute = data_geo_spatial.get_data_points(absolute=True)
    1102         attributes = linear_function(data_points_absolute)
    1103         att = 'spam_and_eggs'
    1104        
    1105         #Create .txt file
    1106         ptsfile = tempfile.mktemp(".txt")
    1107         file = open(ptsfile,"w")
    1108         file.write(" x,y," + att + " \n")
    1109         for data_point, attribute in map(None, data_points_absolute
    1110                                          ,attributes):
    1111             row = str(data_point[0]) + ',' + str(data_point[1]) \
    1112                   + ',' + str(attribute)
    1113             file.write(row + "\n")
    1114         file.close()
    1115 
    1116         #file = open(ptsfile, 'r')
    1117         #lines = file.readlines()
    1118         #file.close()
    1119      
    1120 
    1121         #Check that values can be set from file
    1122         quantity.set_values(filename=ptsfile,
    1123                             attribute_name=att, alpha=0)
    1124         answer = linear_function(quantity.domain.get_vertex_coordinates())
    1125 
    1126         assert allclose(quantity.vertex_values.flat, answer)
    1127 
    1128 
    1129         #Check that values can be set from file using default attribute
    1130         quantity.set_values(filename=ptsfile, alpha=0)
    1131         assert allclose(quantity.vertex_values.flat, answer)
    1132 
    1133         #Cleanup
    1134         import os
    1135         os.remove(ptsfile)
    1136 
    1137 
    1138     def test_set_values_from_file_with_georef2(self):
    1139 
    1140         #Mesh in zone 56 (relative coords)
    1141 
    1142         x0 = 314036.58727982
    1143         y0 = 6224951.2960092
    1144         #x0 = 0.0
    1145         #y0 = 0.0
    1146 
    1147         a = [0.0, 0.0]
    1148         b = [0.0, 2.0]
    1149         c = [2.0, 0.0]
    1150         d = [0.0, 4.0]
    1151         e = [2.0, 2.0]
    1152         f = [4.0, 0.0]
    1153 
    1154         points = [a, b, c, d, e, f]
    1155 
    1156         #bac, bce, ecf, dbe
    1157         elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
    1158 
    1159         mesh4 = Domain(points, elements,
    1160                        geo_reference = Geo_reference(56, x0, y0))
    1161         mesh4.check_integrity()
    1162         quantity = Quantity(mesh4)
    1163 
    1164         #Get (enough) datapoints
    1165         data_points = [[ x0+0.66666667, y0+0.66666667],
    1166                        [ x0+1.33333333, y0+1.33333333],
    1167                        [ x0+2.66666667, y0+0.66666667],
    1168                        [ x0+0.66666667, y0+2.66666667],
    1169                        [ x0+0.0, y0+1.0],
    1170                        [ x0+0.0, y0+3.0],
    1171                        [ x0+1.0, y0+0.0],
    1172                        [ x0+1.0, y0+1.0],
    1173                        [ x0+1.0, y0+2.0],
    1174                        [ x0+1.0, y0+3.0],
    1175                        [ x0+2.0, y0+1.0],
    1176                        [ x0+3.0, y0+0.0],
    1177                        [ x0+3.0, y0+1.0]]
    1178 
    1179 
    1180         data_geo_spatial = Geospatial_data(data_points,
    1181                          geo_reference = Geo_reference(56, 0, 0))
    1182         data_points_absolute = data_geo_spatial.get_data_points(absolute=True)
    1183         attributes = linear_function(data_points_absolute)
    1184         att = 'spam_and_eggs'
    1185        
    1186         #Create .txt file
    1187         ptsfile = tempfile.mktemp(".txt")
    1188         file = open(ptsfile,"w")
    1189         file.write(" x,y," + att + " \n")
    1190         for data_point, attribute in map(None, data_points_absolute
    1191                                          ,attributes):
    1192             row = str(data_point[0]) + ',' + str(data_point[1]) \
    1193                   + ',' + str(attribute)
    1194             file.write(row + "\n")
    1195         file.close()
    1196 
    1197 
    1198         #Check that values can be set from file
    1199         quantity.set_values(filename=ptsfile,
    1200                             attribute_name=att, alpha=0)
    1201         answer = linear_function(quantity.domain. \
    1202                                  get_vertex_coordinates(absolute=True))
    1203 
    1204 
    1205         assert allclose(quantity.vertex_values.flat, answer)
    1206 
    1207 
    1208         #Check that values can be set from file using default attribute
    1209         quantity.set_values(filename=ptsfile, alpha=0)
    1210         assert allclose(quantity.vertex_values.flat, answer)
    1211 
    1212         #Cleanup
    1213         import os
    1214         os.remove(ptsfile)
    1215 
    1216 
    1217 
    1218389
    1219390    def test_set_values_from_quantity(self):
    1220391
    1221         quantity1 = Quantity(self.mesh4)
     392        quantity1 = Quantity(self.domain4)
    1222393        quantity1.set_vertex_values([0,1,2,3,4,5])
    1223394
     
    1226397
    1227398
    1228         quantity2 = Quantity(self.mesh4)
     399        quantity2 = Quantity(self.domain4)
    1229400        quantity2.set_values(quantity=quantity1)
    1230401        assert allclose(quantity2.vertex_values,
     
    1247418
    1248419
    1249     def Xtest_set_values_from_quantity_using_polygon(self):
    1250         """test_set_values_from_quantity_using_polygon(self):
    1251        
    1252         Check that polygon can be used to restrict set_values when
    1253         using another quantity as argument.
    1254         """
    1255        
    1256         # Create restricting polygon (containing node #4 (2,2) and
    1257         # centroid of triangle #1 (bce)
    1258         polygon = [[1.0, 1.0], [4.0, 1.0],
    1259                    [4.0, 4.0], [1.0, 4.0]]
    1260         assert allclose(inside_polygon(self.mesh4.nodes, polygon), 4)                   
    1261        
    1262         quantity1 = Quantity(self.mesh4)
    1263         quantity1.set_vertex_values([0,1,2,3,4,5])
     420
     421
     422    def test_overloading(self):
     423
     424        quantity1 = Quantity(self.domain4)
     425        quantity1.set_values( [[0,1],[1,2],[2,3],[3,4]],
     426                            location = 'vertices')
    1264427
    1265428        assert allclose(quantity1.vertex_values,
    1266                         [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
    1267 
    1268 
    1269         quantity2 = Quantity(self.mesh4)
    1270         quantity2.set_values(quantity=quantity1,
    1271                              polygon=polygon)
    1272                              
    1273         msg = 'Only node #4(e) at (2,2) should have values applied '
    1274         assert allclose(quantity2.vertex_values,
    1275                         [[0,0,0], [0,0,4], [4,0,0], [0,0,4]]), msg       
    1276                         #bac,     bce,     ecf,     dbe
    1277                        
    1278 
    1279 
    1280     def test_overloading(self):
    1281 
    1282         quantity1 = Quantity(self.mesh4)
    1283         quantity1.set_vertex_values([0,1,2,3,4,5])
    1284 
    1285         assert allclose(quantity1.vertex_values,
    1286                         [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
    1287 
    1288 
    1289         quantity2 = Quantity(self.mesh4)
    1290         quantity2.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
     429                        [[0,1], [1,2], [2,3], [3,4]])
     430
     431
     432        quantity2 = Quantity(self.domain4)
     433        quantity2.set_values([[1,2], [5,5], [0,9], [-6, 3]],
    1291434                             location = 'vertices')
    1292435
    1293436
    1294437
    1295         quantity3 = Quantity(self.mesh4)
    1296         quantity3.set_values([[2,2,2], [7,8,9], [7,6,3], [3, 8, -8]],
     438        quantity3 = Quantity(self.domain4)
     439        quantity3.set_values([[2,2], [7,8], [7,6], [3, -8]],
    1297440                             location = 'vertices')
    1298441
     
    1302445        assert allclose(Q.vertex_values, -quantity1.vertex_values)
    1303446        assert allclose(Q.centroid_values, -quantity1.centroid_values)
    1304         assert allclose(Q.edge_values, -quantity1.edge_values)
     447
    1305448
    1306449        # Addition
     
    1308451        assert allclose(Q.vertex_values, quantity1.vertex_values + 7)
    1309452        assert allclose(Q.centroid_values, quantity1.centroid_values + 7)
    1310         assert allclose(Q.edge_values, quantity1.edge_values + 7)
    1311453
    1312454        Q = 7 + quantity1
    1313455        assert allclose(Q.vertex_values, quantity1.vertex_values + 7)
    1314456        assert allclose(Q.centroid_values, quantity1.centroid_values + 7)
    1315         assert allclose(Q.edge_values, quantity1.edge_values + 7)
    1316457
    1317458        Q = quantity1 + quantity2
     
    1320461        assert allclose(Q.centroid_values,
    1321462                        quantity1.centroid_values + quantity2.centroid_values)
    1322         assert allclose(Q.edge_values,
    1323                         quantity1.edge_values + quantity2.edge_values)
    1324 
    1325463
    1326464        Q = quantity1 + quantity2 - 3
     
    1336474        assert allclose(Q.vertex_values, quantity1.vertex_values*3)
    1337475        assert allclose(Q.centroid_values, quantity1.centroid_values*3)
    1338         assert allclose(Q.edge_values, quantity1.edge_values*3)
     476
    1339477        Q = 3*quantity1
    1340478        assert allclose(Q.vertex_values, quantity1.vertex_values*3)
     
    1342480        #Multiplication
    1343481        Q = quantity1 * quantity2
    1344         #print Q.vertex_values
    1345         #print Q.centroid_values
    1346         #print quantity1.centroid_values
    1347         #print quantity2.centroid_values
    1348 
    1349482        assert allclose(Q.vertex_values,
    1350483                        quantity1.vertex_values * quantity2.vertex_values)
     
    1403536                        quantity2.vertex_values**2)**0.5)
    1404537
    1405 
    1406 
    1407 
    1408 
    1409 
    1410 
    1411538    def test_compute_gradient(self):
    1412         quantity = Quantity(self.mesh4)
     539        quantity = Quantity(self.domain6)
    1413540
    1414541        #Set up for a gradient of (2,0) at mid triangle
    1415         quantity.set_values([2.0, 4.0, 6.0, 2.0],
     542        quantity.set_values([2.0, 4.0, 4.0, 5.0, 10.0, 12.0],
    1416543                            location = 'centroids')
    1417544
     
    1420547        quantity.compute_gradients()
    1421548
    1422         a = quantity.x_gradient
    1423         b = quantity.y_gradient
    1424         #print self.mesh4.centroid_coordinates
    1425         #print a, b
    1426 
    1427         #The central triangle (1)
    1428         #(using standard gradient based on neigbours controid values)
    1429         assert allclose(a[1], 2.0)
    1430         assert allclose(b[1], 0.0)
    1431 
    1432 
    1433         #Left triangle (0) using two point gradient
    1434         #q0 = q1 + a*(x0-x1) + b*(y0-y1)  <=>
    1435         #2  = 4  + a*(-2/3)  + b*(-2/3)
    1436         assert allclose(a[0] + b[0], 3)
    1437         #From orthogonality (a*(y0-y1) + b*(x0-x1) == 0)
    1438         assert allclose(a[0] - b[0], 0)
    1439 
    1440 
    1441         #Right triangle (2) using two point gradient
    1442         #q2 = q1 + a*(x2-x1) + b*(y2-y1)  <=>
    1443         #6  = 4  + a*(4/3)  + b*(-2/3)
    1444         assert allclose(2*a[2] - b[2], 3)
    1445         #From orthogonality (a*(y1-y2) + b*(x2-x1) == 0)
    1446         assert allclose(a[2] + 2*b[2], 0)
    1447 
    1448 
    1449         #Top triangle (3) using two point gradient
    1450         #q3 = q1 + a*(x3-x1) + b*(y3-y1)  <=>
    1451         #2  = 4  + a*(-2/3)  + b*(4/3)
    1452         assert allclose(a[3] - 2*b[3], 3)
    1453         #From orthogonality (a*(y1-y3) + b*(x3-x1) == 0)
    1454         assert allclose(2*a[3] + b[3], 0)
    1455 
    1456 
    1457 
    1458         #print a, b
     549        a = quantity.gradients
     550
     551        assert allclose(a, [ 3., 1., 0.5, 3., 3.5, 0.5])
     552       
    1459553        quantity.extrapolate_second_order()
    1460554
    1461         #Apply q(x,y) = qc + a*(x-xc) + b*(y-yc)
    1462         assert allclose(quantity.vertex_values[0,:], [3., 0.,  3.])
    1463         assert allclose(quantity.vertex_values[1,:], [4./3, 16./3,  16./3])
    1464 
    1465 
    1466         #a = 1.2, b=-0.6
    1467         #q(4,0) = 6 + a*(4 - 8/3) + b*(-2/3)
    1468         assert allclose(quantity.vertex_values[2,2], 8)
    1469 
    1470     def test_get_gradients(self):
    1471         quantity = Quantity(self.mesh4)
    1472 
    1473         #Set up for a gradient of (2,0) at mid triangle
    1474         quantity.set_values([2.0, 4.0, 6.0, 2.0],
    1475                             location = 'centroids')
    1476 
    1477 
    1478         #Gradients
    1479         quantity.compute_gradients()
    1480 
    1481         a, b = quantity.get_gradients()
    1482         #print self.mesh4.centroid_coordinates
    1483         #print a, b
    1484 
    1485         #The central triangle (1)
    1486         #(using standard gradient based on neigbours controid values)
    1487         assert allclose(a[1], 2.0)
    1488         assert allclose(b[1], 0.0)
    1489 
    1490 
    1491         #Left triangle (0) using two point gradient
    1492         #q0 = q1 + a*(x0-x1) + b*(y0-y1)  <=>
    1493         #2  = 4  + a*(-2/3)  + b*(-2/3)
    1494         assert allclose(a[0] + b[0], 3)
    1495         #From orthogonality (a*(y0-y1) + b*(x0-x1) == 0)
    1496         assert allclose(a[0] - b[0], 0)
    1497 
    1498 
    1499         #Right triangle (2) using two point gradient
    1500         #q2 = q1 + a*(x2-x1) + b*(y2-y1)  <=>
    1501         #6  = 4  + a*(4/3)  + b*(-2/3)
    1502         assert allclose(2*a[2] - b[2], 3)
    1503         #From orthogonality (a*(y1-y2) + b*(x2-x1) == 0)
    1504         assert allclose(a[2] + 2*b[2], 0)
    1505 
    1506 
    1507         #Top triangle (3) using two point gradient
    1508         #q3 = q1 + a*(x3-x1) + b*(y3-y1)  <=>
    1509         #2  = 4  + a*(-2/3)  + b*(4/3)
    1510         assert allclose(a[3] - 2*b[3], 3)
    1511         #From orthogonality (a*(y1-y3) + b*(x3-x1) == 0)
    1512         assert allclose(2*a[3] + b[3], 0)
    1513 
     555
     556        assert allclose(quantity.vertex_values, [[ 1., 3. ],
     557                                                 [ 4., 4. ],
     558                                                 [ 4., 4. ],
     559                                                 [ 4., 6.],
     560                                                 [ 8.25, 11.75],
     561                                                 [ 11.,  13. ]])
     562
     563                                                 
    1514564
    1515565    def test_second_order_extrapolation2(self):
    1516         quantity = Quantity(self.mesh4)
     566        quantity = Quantity(self.domain4)
    1517567
    1518568        #Set up for a gradient of (3,1), f(x) = 3x+y
     
    1543593
    1544594    def test_backup_saxpy_centroid_values(self):
    1545         quantity = Quantity(self.mesh4)
     595        quantity = Quantity(self.domain4)
    1546596
    1547597        #Set up for a gradient of (3,1), f(x) = 3x+y
     
    1566616
    1567617    def test_first_order_extrapolator(self):
    1568         quantity = Quantity(self.mesh4)
    1569 
    1570         #Test centroids
    1571         quantity.set_values([1.,2.,3.,4.], location = 'centroids')
    1572         assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
     618        quantity = Quantity(self.domain4)
     619
     620        centroid_values = array([1.,2.,3.,4.])
     621        quantity.set_values(centroid_values, location = 'centroids')
     622        assert allclose(quantity.centroid_values, centroid_values) #Centroid
    1573623
    1574624        #Extrapolate
     
    1576626
    1577627        #Check that gradient is zero
    1578         a,b = quantity.get_gradients()
     628        a = quantity.gradients
    1579629        assert allclose(a, [0,0,0,0])
    1580         assert allclose(b, [0,0,0,0])
     630       
    1581631
    1582632        #Check vertices but not edge values
    1583633        assert allclose(quantity.vertex_values,
    1584                         [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
     634                        [[1,1], [2,2], [3,3], [4,4]])
    1585635
    1586636
    1587637    def test_second_order_extrapolator(self):
    1588         quantity = Quantity(self.mesh4)
     638        quantity = Quantity(self.domain4)
    1589639
    1590640        #Set up for a gradient of (3,0) at mid triangle
     
    1616666
    1617667
    1618 
    1619 
    1620 
    1621     def test_limit_vertices_by_all_neighbours(self):
    1622         quantity = Quantity(self.mesh4)
     668    def test_limit(self):
     669        quantity = Quantity(self.domain4)
    1623670
    1624671        #Create a deliberate overshoot (e.g. from gradient computation)
    1625         quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
    1626 
     672        quantity.set_values([[0,0], [2,20], [-20,3], [8,3]])
    1627673
    1628674        #Limit
    1629         quantity.limit_vertices_by_all_neighbours()
    1630 
    1631         #Assert that central triangle is limited by neighbours
    1632         assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0]
    1633         assert quantity.vertex_values[1,0] <= quantity.vertex_values[3,1]
    1634 
    1635         assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1]
    1636         assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
    1637 
    1638         assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
    1639         assert quantity.vertex_values[1,2] <= quantity.vertex_values[3,1]
    1640 
    1641 
    1642 
    1643         #Assert that quantities are conserved
    1644         from Numeric import sum
    1645         for k in range(quantity.centroid_values.shape[0]):
    1646             assert allclose (quantity.centroid_values[k],
    1647                              sum(quantity.vertex_values[k,:])/3)
    1648 
    1649 
    1650 
    1651     def test_limit_edges_by_all_neighbours(self):
    1652         quantity = Quantity(self.mesh4)
    1653 
    1654         #Create a deliberate overshoot (e.g. from gradient computation)
    1655         quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
    1656 
    1657 
    1658         #Limit
    1659         quantity.limit_edges_by_all_neighbours()
    1660 
    1661         #Assert that central triangle is limited by neighbours
    1662         assert quantity.edge_values[1,0] <= quantity.centroid_values[2]
    1663         assert quantity.edge_values[1,0] >= quantity.centroid_values[0]
    1664 
    1665         assert quantity.edge_values[1,1] <= quantity.centroid_values[2]
    1666         assert quantity.edge_values[1,1] >= quantity.centroid_values[0]
    1667 
    1668         assert quantity.edge_values[1,2] <= quantity.centroid_values[2]
    1669         assert quantity.edge_values[1,2] >= quantity.centroid_values[0]
    1670 
    1671 
    1672 
    1673         #Assert that quantities are conserved
    1674         from Numeric import sum
    1675         for k in range(quantity.centroid_values.shape[0]):
    1676             assert allclose (quantity.centroid_values[k],
    1677                              sum(quantity.vertex_values[k,:])/3)
    1678 
    1679 
    1680     def test_limit_edges_by_neighbour(self):
    1681         quantity = Quantity(self.mesh4)
    1682 
    1683         #Create a deliberate overshoot (e.g. from gradient computation)
    1684         quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
    1685 
    1686 
    1687         #Limit
    1688         quantity.limit_edges_by_neighbour()
    1689 
    1690         #Assert that central triangle is limited by neighbours
    1691         assert quantity.edge_values[1,0] <= quantity.centroid_values[3]
    1692         assert quantity.edge_values[1,0] >= quantity.centroid_values[1]
    1693 
    1694         assert quantity.edge_values[1,1] <= quantity.centroid_values[2]
    1695         assert quantity.edge_values[1,1] >= quantity.centroid_values[1]
    1696 
    1697         assert quantity.edge_values[1,2] <= quantity.centroid_values[1]
    1698         assert quantity.edge_values[1,2] >= quantity.centroid_values[0]
    1699 
    1700 
    1701 
    1702         #Assert that quantities are conserved
    1703         from Numeric import sum
    1704         for k in range(quantity.centroid_values.shape[0]):
    1705             assert allclose (quantity.centroid_values[k],
    1706                              sum(quantity.vertex_values[k,:])/3)
    1707 
    1708     def test_limiter2(self):
    1709         """Taken from test_shallow_water
    1710         """
    1711         quantity = Quantity(self.mesh4)
    1712         quantity.domain.beta_w = 0.9
    1713        
     675        quantity.limit_minmod()
     676
     677       
     678        #cells 1 and 2 are local max and min
     679        assert quantity.vertex_values[1][0] == quantity.centroid_values[1]
     680        assert quantity.vertex_values[1][1] == quantity.centroid_values[1]
     681
     682        assert quantity.vertex_values[2][0] == quantity.centroid_values[2]
     683        assert quantity.vertex_values[2][1] == quantity.centroid_values[2]
     684
     685
     686
     687    def test_distribute_first_order(self):
     688        quantity = Quantity(self.domain4)
     689
    1714690        #Test centroids
    1715         quantity.set_values([2.,4.,8.,2.], location = 'centroids')
    1716         assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid
     691        centroid_values = array([1.,2.,3.,4.])
     692        quantity.set_values(centroid_values, location = 'centroids')
     693        assert allclose(quantity.centroid_values, centroid_values) #Centroid
     694
     695
     696        #Extrapolate from centroid to vertices and edges
     697        quantity.extrapolate_first_order()
     698       
     699        assert allclose(quantity.vertex_values,[[ 1.,  1.],
     700                                                [ 2.,  2.],
     701                                                [ 3.,  3.],
     702                                                [ 4.,  4.]])
     703 
     704
     705
     706    def test_distribute_second_order(self):
     707        quantity = Quantity(self.domain4)
     708
     709        #Test centroids
     710        centroid_values = array([2.,4.,8.,2.])
     711        quantity.set_values(centroid_values, location = 'centroids')
     712        assert allclose(quantity.centroid_values, centroid_values) #Centroid
    1717713
    1718714
     
    1720716        quantity.extrapolate_second_order()
    1721717
    1722         assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6])
    1723 
    1724         #Limit
    1725         quantity.limit()
    1726 
    1727         # limited value for beta_w = 0.9
    1728        
    1729         assert allclose(quantity.vertex_values[1,:], [2.2, 4.9, 4.9])
    1730         # limited values for beta_w = 0.5
    1731         #assert allclose(quantity.vertex_values[1,:], [3.0, 4.5, 4.5])
    1732 
    1733 
    1734         #Assert that quantities are conserved
    1735         from Numeric import sum
    1736         for k in range(quantity.centroid_values.shape[0]):
    1737             assert allclose (quantity.centroid_values[k],
    1738                              sum(quantity.vertex_values[k,:])/3)
    1739 
    1740 
    1741 
    1742 
    1743 
    1744     def test_distribute_first_order(self):
    1745         quantity = Quantity(self.mesh4)
     718        assert allclose(quantity.vertex_values, [[ 1.2,  2.8],
     719                                                 [ 2.,  6. ],
     720                                                 [ 8.,   8. ],
     721                                                 [ 6.8, -2.8]])
     722
     723
     724    def test_update_explicit(self):
     725        quantity = Quantity(self.domain4)
    1746726
    1747727        #Test centroids
     
    1749729        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
    1750730
    1751 
    1752         #Extrapolate from centroid to vertices and edges
    1753         quantity.extrapolate_first_order()
    1754 
    1755         #Interpolate
    1756         #quantity.interpolate_from_vertices_to_edges()
    1757 
    1758         assert allclose(quantity.vertex_values,
    1759                         [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
    1760         assert allclose(quantity.edge_values, [[1,1,1], [2,2,2],
    1761                                                [3,3,3], [4, 4, 4]])
    1762 
    1763 
    1764     def test_interpolate_from_vertices_to_edges(self):
    1765         quantity = Quantity(self.mesh4)
    1766 
    1767         quantity.vertex_values = array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]],Float)
    1768 
    1769         quantity.interpolate_from_vertices_to_edges()
    1770 
    1771         assert allclose(quantity.edge_values, [[1., 1.5, 0.5],
    1772                                                [3., 2.5, 1.5],
    1773                                                [3.5, 4.5, 3.],
    1774                                                [2.5, 3.5, 2]])
    1775 
    1776 
    1777     def test_interpolate_from_edges_to_vertices(self):
    1778         quantity = Quantity(self.mesh4)
    1779 
    1780         quantity.edge_values = array([[1., 1.5, 0.5],
    1781                                 [3., 2.5, 1.5],
    1782                                 [3.5, 4.5, 3.],
    1783                                 [2.5, 3.5, 2]],Float)
    1784 
    1785         quantity.interpolate_from_edges_to_vertices()
    1786 
    1787         assert allclose(quantity.vertex_values,
    1788                         [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
    1789 
    1790 
    1791 
    1792     def test_distribute_second_order(self):
    1793         quantity = Quantity(self.mesh4)
    1794 
    1795         #Test centroids
    1796         quantity.set_values([2.,4.,8.,2.], location = 'centroids')
    1797         assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid
    1798 
    1799 
    1800         #Extrapolate
    1801         quantity.extrapolate_second_order()
    1802 
    1803         assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6])
    1804 
    1805 
    1806     def test_update_explicit(self):
    1807         quantity = Quantity(self.mesh4)
     731        #Set explicit_update
     732        quantity.explicit_update = array( [1.,1.,1.,1.] )
     733
     734        #Update with given timestep
     735        quantity.update(0.1)
     736
     737        x = array([1, 2, 3, 4]) + array( [.1,.1,.1,.1] )
     738        assert allclose( quantity.centroid_values, x)
     739
     740    def test_update_semi_implicit(self):
     741        quantity = Quantity(self.domain4)
    1808742
    1809743        #Test centroids
     
    1811745        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
    1812746
    1813         #Set explicit_update
    1814         quantity.explicit_update = array( [1.,1.,1.,1.] )
    1815 
    1816         #Update with given timestep
    1817         quantity.update(0.1)
    1818 
    1819         x = array([1, 2, 3, 4]) + array( [.1,.1,.1,.1] )
    1820         assert allclose( quantity.centroid_values, x)
    1821 
    1822     def test_update_semi_implicit(self):
    1823         quantity = Quantity(self.mesh4)
    1824 
    1825         #Test centroids
    1826         quantity.set_values([1.,2.,3.,4.], location = 'centroids')
    1827         assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
    1828 
    1829747        #Set semi implicit update
    1830748        quantity.semi_implicit_update = array([1.,1.,1.,1.])
     
    1842760
    1843761    def test_both_updates(self):
    1844         quantity = Quantity(self.mesh4)
     762        quantity = Quantity(self.domain4)
    1845763
    1846764        #Test centroids
    1847         quantity.set_values([1.,2.,3.,4.], location = 'centroids')
    1848         assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
     765        centroid_values = array( [1, 2, 3, 4] )
     766        quantity.set_values(centroid_values, location = 'centroids')
     767        assert allclose(quantity.centroid_values, centroid_values) #Centroid
    1849768
    1850769        #Set explicit_update
    1851         quantity.explicit_update = array( [4.,3.,2.,1.] )
     770        explicit_update = array( [4.,3.,2.,1.] )
     771        quantity.explicit_update[:] = explicit_update
    1852772
    1853773        #Set semi implicit update
    1854         quantity.semi_implicit_update = array( [1.,1.,1.,1.] )
     774        semi_implicit_update = array( [1.,1.,1.,1.] )
     775        quantity.semi_implicit_update[:] = semi_implicit_update
    1855776
    1856777        #Update with given timestep
     
    1858779        quantity.update(0.1)
    1859780
    1860         sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])
    1861         denom = ones(4, Float)-timestep*sem
    1862 
    1863         x = array([1., 2., 3., 4.])
    1864         x /= denom
    1865         x += timestep*array( [4.0, 3.0, 2.0, 1.0] )
    1866 
     781        denom = 1.0-timestep*semi_implicit_update
     782        x = (centroid_values + timestep*explicit_update)/denom
     783 
    1867784        assert allclose( quantity.centroid_values, x)
    1868 
    1869 
    1870 
    1871785
    1872786    #Test smoothing
    1873787    def test_smoothing(self):
    1874788
    1875         from mesh_factory import rectangular
     789       
    1876790        from shallow_water import Domain, Transmissive_boundary
    1877791        from Numeric import zeros, Float
    1878792        from anuga.utilities.numerical_tools import mean
    1879793
    1880         #Create basic mesh
    1881         points, vertices, boundary = rectangular(2, 2)
    1882794
    1883795        #Create shallow water domain
    1884         domain = Domain(points, vertices, boundary)
     796        domain = Domain(points10)
    1885797        domain.default_order=2
    1886798        domain.reduction = mean
     
    1888800
    1889801        #Set some field values
    1890         domain.set_quantity('elevation', lambda x,y: x)
     802        domain.set_quantity('elevation', lambda x: x)
    1891803        domain.set_quantity('friction', 0.03)
    1892804
     
    2030942
    2031943
    2032     def set_array_values_by_index(self):
    2033 
    2034         from mesh_factory import rectangular
    2035         from shallow_water import Domain
    2036         from Numeric import zeros, Float
    2037 
    2038         #Create basic mesh
    2039         points, vertices, boundary = rectangular(1, 1)
    2040 
    2041         #Create shallow water domain
    2042         domain = Domain(points, vertices, boundary)
    2043         #print "domain.number_of_elements ",domain.number_of_elements
    2044         quantity = Quantity(domain,[[1,1,1],[2,2,2]])
    2045         value = [7]
    2046         indices = [1]
    2047         quantity.set_array_values_by_index(value,
    2048                                            location = 'centroids',
    2049                                            indices = indices)
    2050         #print "quantity.centroid_values",quantity.centroid_values
    2051 
    2052         assert allclose(quantity.centroid_values, [1,7])
    2053 
    2054         quantity.set_array_values([15,20,25], indices = indices)
    2055         assert allclose(quantity.centroid_values, [1,20])
    2056 
    2057         quantity.set_array_values([15,20,25], indices = indices)
    2058         assert allclose(quantity.centroid_values, [1,20])
    2059 
    2060     def test_setting_some_vertex_values(self):
    2061         """
    2062         set 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         #print "vertices",vertices
    2071         #Create shallow water domain
    2072         domain = Domain(points, vertices, boundary)
    2073         #print "domain.number_of_elements ",domain.number_of_elements
    2074         quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
    2075                                     [4,4,4],[5,5,5],[6,6,6]])
    2076 
    2077 
    2078         # Check that constants work
    2079         value = 7
    2080         indices = [1]
    2081         quantity.set_values(value,
    2082                             location = 'centroids',
    2083                             indices = indices)
    2084         #print "quantity.centroid_values",quantity.centroid_values
    2085         assert allclose(quantity.centroid_values, [1,7,3,4,5,6])
    2086        
    2087         value = [7]
    2088         indices = [1]
    2089         quantity.set_values(value,
    2090                             location = 'centroids',
    2091                             indices = indices)
    2092         #print "quantity.centroid_values",quantity.centroid_values
    2093         assert allclose(quantity.centroid_values, [1,7,3,4,5,6])
    2094 
    2095         value = [[15,20,25]]
    2096         quantity.set_values(value, indices = indices)
    2097         #print "1 quantity.vertex_values",quantity.vertex_values
    2098         assert allclose(quantity.vertex_values[1], value[0])
    2099 
    2100 
    2101         #print "quantity",quantity.vertex_values
    2102         values = [10,100,50]
    2103         quantity.set_values(values, indices = [0,1,5], location = 'centroids')
    2104         #print "2 quantity.vertex_values",quantity.vertex_values
    2105         assert allclose(quantity.vertex_values[0], [10,10,10])
    2106         assert allclose(quantity.vertex_values[5], [50,50,50])
    2107         #quantity.interpolate()
    2108         #print "quantity.centroid_values",quantity.centroid_values
    2109         assert allclose(quantity.centroid_values, [10,100,3,4,5,50])
    2110 
    2111 
    2112         quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
    2113                                     [4,4,4],[5,5,5],[6,6,6]])
    2114         values = [10,100,50]
    2115         #this will be per unique vertex, indexing the vertices
    2116         #print "quantity.vertex_values",quantity.vertex_values
    2117         quantity.set_values(values, indices = [0,1,5])
    2118         #print "quantity.vertex_values",quantity.vertex_values
    2119         assert allclose(quantity.vertex_values[0], [1,50,10])
    2120         assert allclose(quantity.vertex_values[5], [6,6,6])
    2121         assert allclose(quantity.vertex_values[1], [100,10,50])
    2122 
    2123         quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
    2124                                     [4,4,4],[5,5,5],[6,6,6]])
    2125         values = [[31,30,29],[400,400,400],[1000,999,998]]
    2126         quantity.set_values(values, indices = [3,3,5])
    2127         quantity.interpolate()
    2128         assert allclose(quantity.centroid_values, [1,2,3,400,5,999])
    2129 
    2130         values = [[1,1,1],[2,2,2],[3,3,3],
    2131                                     [4,4,4],[5,5,5],[6,6,6]]
    2132         quantity.set_values(values)
    2133 
    2134         # testing the standard set values by vertex
    2135         # indexed by vertex_id in general_mesh.coordinates
    2136         values = [0,1,2,3,4,5,6,7]
    2137 
    2138         quantity.set_values(values)
    2139         #print "1 quantity.vertex_values",quantity.vertex_values
    2140         assert allclose(quantity.vertex_values,[[ 4.,  5.,  0.],
    2141                                                 [ 1.,  0.,  5.],
    2142                                                 [ 5.,  6.,  1.],
    2143                                                 [ 2.,  1.,  6.],
    2144                                                 [ 6.,  7.,  2.],
    2145                                                 [ 3.,  2.,  7.]])
    2146 
    2147     def test_setting_unique_vertex_values(self):
    2148         """
    2149         set values based on unique_vertex lists.
    2150         """
    2151         from mesh_factory import rectangular
    2152         from shallow_water import Domain
    2153         from Numeric import zeros, Float
    2154 
    2155         #Create basic mesh
    2156         points, vertices, boundary = rectangular(1, 3)
    2157         #print "vertices",vertices
    2158         #Create shallow water domain
    2159         domain = Domain(points, vertices, boundary)
    2160         #print "domain.number_of_elements ",domain.number_of_elements
    2161         quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
    2162                                     [4,4,4],[5,5,5]])
    2163         value = 7
    2164         indices = [1,5]
    2165         quantity.set_values(value,
    2166                             location = 'unique vertices',
    2167                             indices = indices)
    2168         #print "quantity.centroid_values",quantity.centroid_values
    2169         assert allclose(quantity.vertex_values[0], [0,7,0])
    2170         assert allclose(quantity.vertex_values[1], [7,1,7])
    2171         assert allclose(quantity.vertex_values[2], [7,2,7])
    2172 
    2173 
    2174     def test_get_values(self):
    2175         """
    2176         get values based on triangle lists.
    2177         """
    2178         from mesh_factory import rectangular
    2179         from shallow_water import Domain
    2180         from Numeric import zeros, Float
    2181 
    2182         #Create basic mesh
    2183         points, vertices, boundary = rectangular(1, 3)
    2184 
    2185         #print "points",points
    2186         #print "vertices",vertices
    2187         #print "boundary",boundary
    2188 
    2189         #Create shallow water domain
    2190         domain = Domain(points, vertices, boundary)
    2191         #print "domain.number_of_elements ",domain.number_of_elements
    2192         quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
    2193                                     [4,4,4],[5,5,5]])
    2194 
    2195         #print "quantity.get_values(location = 'unique vertices')", \
    2196         #      quantity.get_values(location = 'unique vertices')
    2197 
    2198         #print "quantity.get_values(location = 'unique vertices')", \
    2199         #      quantity.get_values(indices=[0,1,2,3,4,5,6,7], \
    2200         #                          location = 'unique vertices')
    2201 
    2202         answer = [0.5,2,4,5,0,1,3,4.5]
    2203         assert allclose(answer,
    2204                         quantity.get_values(location = 'unique vertices'))
    2205 
    2206         indices = [0,5,3]
    2207         answer = [0.5,1,5]
    2208         assert allclose(answer,
    2209                         quantity.get_values(indices=indices, \
    2210                                             location = 'unique vertices'))
    2211         #print "quantity.centroid_values",quantity.centroid_values
    2212         #print "quantity.get_values(location = 'centroids') ",\
    2213         #      quantity.get_values(location = 'centroids')
    2214 
    2215 
    2216 
    2217 
    2218     def test_get_values_2(self):
    2219         """Different mesh (working with domain object) - also check centroids.
    2220         """
    2221 
    2222        
    2223         a = [0.0, 0.0]
    2224         b = [0.0, 2.0]
    2225         c = [2.0,0.0]
    2226         d = [0.0, 4.0]
    2227         e = [2.0, 2.0]
    2228         f = [4.0,0.0]
    2229 
    2230         points = [a, b, c, d, e, f]
    2231         #bac, bce, ecf, dbe
    2232         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    2233 
    2234         domain = Domain(points, vertices)
    2235 
    2236         quantity = Quantity(domain)
    2237         quantity.set_values(lambda x, y: x+2*y) #2 4 4 6
    2238        
    2239         assert allclose(quantity.get_values(location='centroids'), [2,4,4,6])
    2240         assert allclose(quantity.get_values(location='centroids', indices=[1,3]), [4,6])
    2241 
    2242 
    2243         assert allclose(quantity.get_values(location='vertices'), [[4,0,2],
    2244                                                                    [4,2,6],
    2245                                                                    [6,2,4],
    2246                                                                    [8,4,6]])
    2247        
    2248         assert allclose(quantity.get_values(location='vertices', indices=[1,3]), [[4,2,6],
    2249                                                                                   [8,4,6]])
    2250 
    2251 
    2252         assert allclose(quantity.get_values(location='edges'), [[1,3,2],
    2253                                                                 [4,5,3],
    2254                                                                 [3,5,4],
    2255                                                                 [5,7,6]])
    2256         assert allclose(quantity.get_values(location='edges', indices=[1,3]),
    2257                         [[4,5,3],
    2258                          [5,7,6]])       
    2259 
    2260         # Check averaging over vertices
    2261         #a: 0
    2262         #b: (4+4+4)/3
    2263         #c: (2+2+2)/3
    2264         #d: 8
    2265         #e: (6+6+6)/3       
    2266         #f: 4
    2267         assert(quantity.get_values(location='unique vertices'), [0, 4, 2, 8, 6, 4])       
    2268                                                                                  
    2269        
    2270 
    2271 
    2272 
    2273 
    2274     def test_get_interpolated_values(self):
    2275 
    2276         from mesh_factory import rectangular
    2277         from shallow_water import Domain
    2278         from Numeric import zeros, Float
    2279 
    2280         #Create basic mesh
    2281         points, vertices, boundary = rectangular(1, 3)
    2282         domain = Domain(points, vertices, boundary)
    2283 
    2284         #Constant values
    2285         quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
    2286                                     [4,4,4],[5,5,5]])
    2287 
    2288        
    2289 
    2290         # Get interpolated values at centroids
    2291         interpolation_points = domain.get_centroid_coordinates()
    2292         answer = quantity.get_values(location='centroids')
    2293 
    2294        
    2295         #print quantity.get_values(points=interpolation_points)
    2296         assert allclose(answer, quantity.get_values(interpolation_points=interpolation_points))
    2297 
    2298 
    2299         #Arbitrary values
    2300         quantity = Quantity(domain,[[0,1,2],[3,1,7],[2,1,2],[3,3,7],
    2301                                     [1,4,-9],[2,5,0]])
    2302 
    2303 
    2304         # Get interpolated values at centroids
    2305         interpolation_points = domain.get_centroid_coordinates()
    2306         answer = quantity.get_values(location='centroids')
    2307         #print answer
    2308         #print quantity.get_values(points=interpolation_points)
    2309         assert allclose(answer, quantity.get_values(interpolation_points=interpolation_points))       
    2310                        
    2311 
    2312         #FIXME TODO
    2313         #indices = [0,5,3]
    2314         #answer = [0.5,1,5]
    2315         #assert allclose(answer,
    2316         #                quantity.get_values(indices=indices, \
    2317         #                                    location = 'unique vertices'))
    2318 
    2319 
    2320 
    2321 
    2322     def test_get_interpolated_values_2(self):
    2323         a = [0.0, 0.0]
    2324         b = [0.0, 2.0]
    2325         c = [2.0,0.0]
    2326         d = [0.0, 4.0]
    2327         e = [2.0, 2.0]
    2328         f = [4.0,0.0]
    2329 
    2330         points = [a, b, c, d, e, f]
    2331         #bac, bce, ecf, dbe
    2332         vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    2333 
    2334         domain = Domain(points, vertices)
    2335 
    2336         quantity = Quantity(domain)
    2337         quantity.set_values(lambda x, y: x+2*y) #2 4 4 6
    2338 
    2339         #First pick one point
    2340         x, y = 2.0/3, 8.0/3
    2341         v = quantity.get_values(interpolation_points = [[x,y]])
    2342         assert allclose(v, 6)       
    2343 
    2344         # Then another to test that algorithm won't blindly
    2345         # reuse interpolation matrix
    2346         x, y = 4.0/3, 4.0/3
    2347         v = quantity.get_values(interpolation_points = [[x,y]])
    2348         assert allclose(v, 4)       
    2349 
    2350 
    2351 
    2352 
    2353     def test_getting_some_vertex_values(self):
    2354         """
    2355         get values based on triangle lists.
    2356         """
    2357         from mesh_factory import rectangular
    2358         from shallow_water import Domain
    2359         from Numeric import zeros, Float
    2360 
    2361         #Create basic mesh
    2362         points, vertices, boundary = rectangular(1, 3)
    2363 
    2364         #print "points",points
    2365         #print "vertices",vertices
    2366         #print "boundary",boundary
    2367 
    2368         #Create shallow water domain
    2369         domain = Domain(points, vertices, boundary)
    2370         #print "domain.number_of_elements ",domain.number_of_elements
    2371         quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
    2372                                     [4,4,4],[5,5,5],[6,6,6]])
    2373         value = [7]
    2374         indices = [1]
    2375         quantity.set_values(value,
    2376                             location = 'centroids',
    2377                             indices = indices)
    2378         #print "quantity.centroid_values",quantity.centroid_values
    2379         #print "quantity.get_values(location = 'centroids') ",\
    2380         #      quantity.get_values(location = 'centroids')
    2381         assert allclose(quantity.centroid_values,
    2382                         quantity.get_values(location = 'centroids'))
    2383 
    2384 
    2385         value = [[15,20,25]]
    2386         quantity.set_values(value, indices = indices)
    2387         #print "1 quantity.vertex_values",quantity.vertex_values
    2388         assert allclose(quantity.vertex_values, quantity.get_values())
    2389 
    2390         assert allclose(quantity.edge_values,
    2391                         quantity.get_values(location = 'edges'))
    2392 
    2393         # get a subset of elements
    2394         subset = quantity.get_values(location='centroids', indices=[0,5])
    2395         answer = [quantity.centroid_values[0],quantity.centroid_values[5]]
    2396         assert allclose(subset, answer)
    2397 
    2398 
    2399         subset = quantity.get_values(location='edges', indices=[0,5])
    2400         answer = [quantity.edge_values[0],quantity.edge_values[5]]
    2401         #print "subset",subset
    2402         #print "answer",answer
    2403         assert allclose(subset, answer)
    2404 
    2405         subset = quantity.get_values( indices=[1,5])
    2406         answer = [quantity.vertex_values[1],quantity.vertex_values[5]]
    2407         #print "subset",subset
    2408         #print "answer",answer
    2409         assert allclose(subset, answer)
    2410 
    2411     def test_smooth_vertex_values(self):
    2412         """
    2413         get values based on triangle lists.
    2414         """
    2415         from mesh_factory import rectangular
    2416         from shallow_water import Domain
    2417         from Numeric import zeros, Float
    2418 
    2419         #Create basic mesh
    2420         points, vertices, boundary = rectangular(2, 2)
    2421 
    2422         #print "points",points
    2423         #print "vertices",vertices
    2424         #print "boundary",boundary
    2425 
    2426         #Create shallow water domain
    2427         domain = Domain(points, vertices, boundary)
    2428         #print "domain.number_of_elements ",domain.number_of_elements
    2429         quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
    2430                                     [4,4,4],[5,5,5],[6,6,6],[7,7,7]])
    2431 
    2432         #print "quantity.get_values(location = 'unique vertices')", \
    2433         #      quantity.get_values(location = 'unique vertices')
    2434 
    2435         #print "quantity.get_values(location = 'unique vertices')", \
    2436         #      quantity.get_values(indices=[0,1,2,3,4,5,6,7], \
    2437         #                          location = 'unique vertices')
    2438 
    2439         #print quantity.get_values(location = 'unique vertices')
    2440         #print quantity.domain.number_of_triangles_per_node
    2441         #print quantity.vertex_values
    2442        
    2443         #answer = [0.5, 2, 3, 3, 3.5, 4, 4, 5, 6.5]
    2444         #assert allclose(answer,
    2445         #                quantity.get_values(location = 'unique vertices'))
    2446 
    2447         quantity.smooth_vertex_values()
    2448 
    2449         #print quantity.vertex_values
    2450 
    2451 
    2452         answer_vertex_values = [[3,3.5,0.5],[2,0.5,3.5],[3.5,4,2],[3,2,4],
    2453                                 [4,5,3],[3.5,3,5],[5,6.5,3.5],[4,3.5,6.5]]
    2454        
    2455         assert allclose(answer_vertex_values,
    2456                         quantity.vertex_values)
    2457         #print "quantity.centroid_values",quantity.centroid_values
    2458         #print "quantity.get_values(location = 'centroids') ",\
    2459         #      quantity.get_values(location = 'centroids')
    2460 
    2461 
    2462 
     944 
    2463945#-------------------------------------------------------------
    2464946if __name__ == "__main__":
    2465     suite = unittest.makeSuite(Test_Quantity, 'test')   
    2466     #suite = unittest.makeSuite(Test_Quantity, 'test_set_values_from_file_using_polygon')
    2467 
    2468     #suite = unittest.makeSuite(Test_Quantity, 'test_set_vertex_values_using_general_interface_with_subset')
    2469     #print "restricted test"
    2470     #suite = unittest.makeSuite(Test_Quantity,'verbose_test_set_values_from_UTM_pts')
     947    suite = unittest.makeSuite(Test_Quantity, 'test')
    2471948    runner = unittest.TextTestRunner()
    2472949    runner.run(suite)
Note: See TracChangeset for help on using the changeset viewer.