Changeset 5899


Ignore:
Timestamp:
Nov 6, 2008, 12:28:22 PM (15 years ago)
Author:
rwilson
Message:

Initial NumPy? changes (again!).

Location:
anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes
Files:
22 edited

Legend:

Unmodified
Added
Removed

  • anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/domain.py

    r5847 r5899  
     1## Automatically adapted for numpy.oldnumeric Oct 28, 2008 by alter_code1.py
     2
    13"""Class Domain - 2D triangular domains for finite-volume computations of
    24   conservation laws.
    … …  
    810"""
    911
    10 from Numeric import allclose, argmax, zeros, Float
     12##from numpy.oldnumeric import allclose, argmax, zeros, Float
     13from numpy import allclose, argmax, zeros, float
    1114from anuga.config import epsilon
    1215from anuga.config import beta_euler, beta_rk2
    … …  
    105108
    106109        if verbose: print 'Initialising Domain'
    107         from Numeric import zeros, Float, Int, ones
     110##        from numpy.oldnumeric import zeros, Float, Int, ones
     111        from numpy import zeros, float, int, ones
    108112        from quantity import Quantity
    109113
    … …  
    155159        for key in self.full_send_dict:
    156160            buffer_shape = self.full_send_dict[key][0].shape[0]
    157             self.full_send_dict[key].append(zeros( (buffer_shape,self.nsys) ,Float))
     161            self.full_send_dict[key].append(zeros( (buffer_shape,self.nsys) ,float))
    158162
    159163
    160164        for key in self.ghost_recv_dict:
    161165            buffer_shape = self.ghost_recv_dict[key][0].shape[0]
    162             self.ghost_recv_dict[key].append(zeros( (buffer_shape,self.nsys) ,Float))
     166            self.ghost_recv_dict[key].append(zeros( (buffer_shape,self.nsys) ,float))
    163167
    164168
    … …  
    168172        N = len(self) #number_of_elements
    169173        self.number_of_elements = N
    170         self.tri_full_flag = ones(N, Int)
     174        self.tri_full_flag = ones(N, int)
    171175        for i in self.ghost_recv_dict.keys():
    172176            for id in self.ghost_recv_dict[i][0]:
    … …  
    230234        # calculation
    231235        N = len(self) # Number_of_triangles
    232         self.already_computed_flux = zeros((N, 3), Int)
     236        self.already_computed_flux = zeros((N, 3), int)
    233237
    234238        # Storage for maximal speeds computed for each triangle by
    235239        # compute_fluxes
    236240        # This is used for diagnostics only (reset at every yieldstep)
    237         self.max_speed = zeros(N, Float)
     241        self.max_speed = zeros(N, float)
    238242
    239243        if mesh_filename is not None:
    … …  
    266270        """
    267271
    268         from Numeric import zeros, Float
     272##        from numpy.oldnumeric import zeros, Float
     273        from numpy import zeros, float
    269274
    270275        if not (vertex is None or edge is None):
    … …  
    273278            raise msg
    274279
    275         q = zeros( len(self.conserved_quantities), Float)
     280        q = zeros( len(self.conserved_quantities), float)
    276281
    277282        for i, name in enumerate(self.conserved_quantities):
    … …  
    12021207                self.number_of_steps = 0
    12031208                self.number_of_first_order_steps = 0
    1204                 self.max_speed = zeros(N, Float)
     1209                self.max_speed = zeros(N, float)
    12051210
    12061211    def evolve_one_euler_step(self, yieldstep, finaltime):
    … …  
    15681573        """
    15691574
    1570         from Numeric import ones, sum, equal, Float
     1575##        from numpy.oldnumeric import ones, sum, equal, Float
     1576        from numpy import ones, sum, equal, float
    15711577
    15721578        N = len(self) # Number_of_triangles

  • anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/ermapper_grids.py

    r2054 r5899  
     1## Automatically adapted for numpy.oldnumeric Oct 28, 2008 by alter_code1.py
     2
    13#!/usr/bin/env python
    24
    35# from os import open, write, read
    4 import Numeric
    5 
    6 celltype_map = {'IEEE4ByteReal': Numeric.Float32, 'IEEE8ByteReal': Numeric.Float64}
     6##import numpy.oldnumeric as Numeric
     7import numpy
     8
     9celltype_map = {'IEEE4ByteReal': numpy.float32, 'IEEE8ByteReal': numpy.float64}
    710
    811
    … …  
    1114    write_ermapper_grid(ofile, data, header = {}):
    1215   
    13     Function to write a 2D Numeric array to an ERMapper grid.  There are a series of conventions adopted within
     16    Function to write a 2D NumPy array to an ERMapper grid.  There are a series of conventions adopted within
    1417    this code, specifically:
    1518    1)  The registration coordinate for the data is the SW (or lower-left) corner of the data
    1619    2)  The registration coordinates refer to cell centres
    17     3)  The data is a 2D Numeric array with the NW-most data in element (0,0) and the SE-most data in element (N,M)
     20    3)  The data is a 2D NumPy array with the NW-most data in element (0,0) and the SE-most data in element (N,M)
    1821        where N is the last line and M is the last column
    1922    4)  There has been no testng of the use of a rotated grid.  Best to keep data in an NS orientation
    … …  
    2225    ofile:      string - filename for output (note the output will consist of two files
    2326                ofile and ofile.ers.  Either of these can be entered into this function
    24     data:       Numeric.array - 2D array containing the data to be output to the grid
     27    data:       NumPy.array - 2D array containing the data to be output to the grid
    2528    header:     dictionary - contains spatial information about the grid, in particular:
    2629                    header['datum'] datum for the data ('"GDA94"')
    … …  
    5760
    5861    # Check that the data is a 2 dimensional array
    59     data_size = Numeric.shape(data)
     62    data_size = numpy.shape(data)
    6063    assert len(data_size) == 2
    6164   
    … …  
    8386    nrofcellsperlines = int(header['nrofcellsperline'])
    8487    data = read_ermapper_data(data_file)
    85     data = Numeric.reshape(data,(nroflines,nrofcellsperlines))
     88    data = numpy.reshape(data,(nroflines,nrofcellsperlines))
    8689    return data
    8790   
    … …  
    163166    return header                     
    164167
    165 def write_ermapper_data(grid, ofile, data_format = Numeric.Float32):
     168def write_ermapper_data(grid, ofile, data_format = numpy.float32):
    166169
    167170
    … …  
    181184       
    182185   
    183     # Convert the array to data_format (default format is Float32)
     186    # Convert the array to data_format (default format is float32)
    184187    grid_as_float = grid.astype(data_format)
    185188
    … …  
    193196
    194197
    195 def read_ermapper_data(ifile, data_format = Numeric.Float32):
     198def read_ermapper_data(ifile, data_format = numpy.float32):
    196199    # open input file in a binary format and read the input string
    197200    fid = open(ifile,'rb')
    … …  
    199202    fid.close()
    200203
    201     # convert input string to required format (Note default format is Numeric.Float32)
    202     grid_as_float = Numeric.fromstring(input_string,data_format)
     204    # convert input string to required format (Note default format is numpy.float32)
     205    grid_as_float = numpy.fromstring(input_string,data_format)
    203206    return grid_as_float
    204207

  • anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/general_mesh.py

    r5842 r5899  
    1 
    2 from Numeric import concatenate, reshape, take, allclose
    3 from Numeric import array, zeros, Int, Float, sqrt, sum, arange
     1## Automatically adapted for numpy.oldnumeric Oct 28, 2008 by alter_code1.py
     2
     3
     4##from numpy.oldnumeric import concatenate, reshape, take, allclose
     5##from numpy.oldnumeric import array, zeros, Int, Float, sqrt, sum, arange
     6from numpy import concatenate, reshape, take, allclose
     7from numpy import array, zeros, int, float, sqrt, sum, arange
    48
    59from anuga.coordinate_transforms.geo_reference import Geo_reference
    … …  
    9094        if verbose: print 'General_mesh: Building basic mesh structure in ANUGA domain'
    9195
    92         self.triangles = array(triangles, Int)
    93         self.nodes = array(nodes, Float)
     96##        self.triangles = array(triangles, Int)
     97##        self.nodes = array(nodes, Float)
     98        self.triangles = array(triangles, int)
     99        self.nodes = array(nodes, float)
    94100
    95101
    … …  
    138144
    139145        msg = 'Vertex indices reference non-existing coordinate sets'
    140         assert max(self.triangles.flat) < self.nodes.shape[0], msg
     146        assert max(self.triangles.ravel()) < self.nodes.shape[0], msg
    141147
    142148
    … …  
    146152                      max(self.nodes[:,0]), max(self.nodes[:,1]) ]
    147153
    148         self.xy_extent = array(xy_extent, Float)
     154##        self.xy_extent = array(xy_extent, Float)
     155        self.xy_extent = array(xy_extent, float)
    149156
    150157
    151158        # Allocate space for geometric quantities
    152         self.normals = zeros((N, 6), Float)
    153         self.areas = zeros(N, Float)
    154         self.edgelengths = zeros((N, 3), Float)
     159##        self.normals = zeros((N, 6), Float)
     160##        self.areas = zeros(N, Float)
     161##        self.edgelengths = zeros((N, 3), Float)
     162        self.normals = zeros((N, 6), float)
     163        self.areas = zeros(N, float)
     164        self.edgelengths = zeros((N, 3), float)
    155165
    156166        # Get x,y coordinates for all triangles and store
    … …  
    311321            i = triangle_id
    312322            msg = 'triangle_id must be an integer'
     323            print 'type(triangle_id)=%s. triangle_id=%s' % (type(triangle_id), str(triangle_id))
    313324            assert int(i) == i, msg
    314325            assert 0 <= i < self.number_of_triangles
    … …  
    349360
    350361        M = self.number_of_triangles
    351         vertex_coordinates = zeros((3*M, 2), Float)
     362##        vertex_coordinates = zeros((3*M, 2), Float)
     363        vertex_coordinates = zeros((3*M, 2), float)
    352364
    353365        for i in range(M):
    … …  
    376388            indices = range(M)
    377389
    378         return take(self.triangles, indices)
     390        return take(self.triangles, indices, axis=0)
    379391   
    380392
    … …  
    409421       
    410422        #T = reshape(array(range(K)).astype(Int), (M,3))
    411         T = reshape(arange(K).astype(Int), (M,3))  # Faster
     423##        T = reshape(arange(K).astype(Int), (M,3))  # Faster
     424        T = reshape(arange(K).astype(int), (M,3))  # Faster
    412425       
    413426        return T     
    … …  
    442455           
    443456            # Get number of triangles for this node
    444             count = self.number_of_triangles_per_node[node]
     457            count = int(self.number_of_triangles_per_node[node])
    445458
    446459            for i in range(count):
    … …  
    582595        Y = C[:,1:6:2].copy()
    583596
    584         xmin = min(X.flat)
    585         xmax = max(X.flat)
    586         ymin = min(Y.flat)
    587         ymax = max(Y.flat)
     597        xmin = min(X.ravel())
     598        xmax = max(X.ravel())
     599        ymin = min(Y.ravel())
     600        ymax = max(Y.ravel())
    588601        #print "C",C
    589602        return xmin, xmax, ymin, ymax

  • anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py

    r5672 r5899  
     1## Automatically adapted for numpy.oldnumeric Oct 28, 2008 by alter_code1.py
     2
    13"""boundary.py - Classes for implementing boundary conditions
    24"""
    … …  
    6769            raise Exception, msg
    6870
    69         from Numeric import array, Float
    70         self.conserved_quantities=array(conserved_quantities).astype(Float)
     71#        from numpy.oldnumeric import array, Float
     72        from numpy import array, float
     73        self.conserved_quantities=array(conserved_quantities).astype(float)
    7174
    7275    def __repr__(self):
    … …  
    97100
    98101
    99         from Numeric import array, Float
     102#        from numpy.oldnumeric import array, Float
     103        from numpy import array, float
    100104        try:
    101             q = array(q).astype(Float)
     105            q = array(q).astype(float)
    102106        except:
    103107            msg = 'Return value from time boundary function could '
    … …  
    136140    def __init__(self, filename, domain):
    137141        import time
    138         from Numeric import array
     142#        from numpy.oldnumeric import array
     143        from numpy import array
    139144        from anuga.config import time_format
    140145        from anuga.abstract_2d_finite_volumes.util import File_function
    … …  
    206211
    207212        import time
    208         from Numeric import array, zeros, Float
     213#        from numpy.oldnumeric import array, zeros, Float
     214        from numpy import array, zeros, float
    209215        from anuga.config import time_format
    210216        from anuga.abstract_2d_finite_volumes.util import file_function
    … …  
    222228
    223229        if verbose: print 'Find midpoint coordinates of entire boundary'
    224         self.midpoint_coordinates = zeros( (len(domain.boundary), 2), Float)
     230        self.midpoint_coordinates = zeros( (len(domain.boundary), 2), float)
    225231        boundary_keys = domain.boundary.keys()
    226232

  • anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/mesh_factory.py

    r3678 r5899  
     1## Automatically adapted for numpy.oldnumeric Oct 28, 2008 by alter_code1.py
     2
    13"""Library of standard meshes and facilities for reading various
    24mesh file formats
    … …  
    7375
    7476    from anuga.config import epsilon
    75     from Numeric import zeros, Float, Int
     77#    from numpy.oldnumeric import zeros, Float, Int
     78    from numpy import zeros, float, int
    7679
    7780    delta1 = float(len1)/m
    … …  
    9396    index = Index(n,m)
    9497
    95     points = zeros( (Np,2), Float)
     98    points = zeros( (Np,2), float)
    9699
    97100    for i in range(m+1):
    … …  
    105108
    106109
    107     elements = zeros( (Nt,3), Int)
     110    elements = zeros( (Nt,3), int)
    108111    boundary = {}
    109112    nt = -1
    … …  
    149152
    150153    from anuga.config import epsilon
    151     from Numeric import zeros, Float, Int
     154#    from numpy.oldnumeric import zeros, Float, Int
     155    from numpy import zeros, float, int
    152156
    153157    delta1 = float(len1)/m
    … …  
    208212    """
    209213
    210     from Numeric import array
     214#    from numpy.oldnumeric import array
     215    from numpy import array
    211216    import math
    212217
    … …  
    510515    """
    511516
    512     from Numeric import array
     517#    from numpy.oldnumeric import array
     518    from numpy import array
    513519    import math
    514520
    … …  
    592598    """
    593599
    594     from Numeric import array
     600#    from numpy.oldnumeric import array
     601    from numpy import array
    595602    import math
    596603
    … …  
    686693    """
    687694
    688     from Numeric import array
     695#    from numpy.oldnumeric import array
     696    from numpy import array
    689697    import math
    690698

  • anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/neighbour_mesh.py

    r5866 r5899  
     1## Automatically adapted for numpy.oldnumeric Oct 28, 2008 by alter_code1.py
     2
    13"""Classes implementing general 2D triangular mesh with neighbour structure.
    24
    … …  
    1113from anuga.caching import cache
    1214from math import pi, sqrt
    13 from Numeric import array, allclose
     15#from numpy.oldnumeric import array, allclose
     16from numpy import array, allclose
    1417       
    1518
    … …  
    8184
    8285
    83         from Numeric import array, zeros, Int, Float, maximum, sqrt, sum
     86#        from numpy.oldnumeric import array, zeros, Int, Float, maximum, sqrt, sum
     87        from numpy import array, zeros, int, float, maximum, sqrt, sum
    8488
    8589        General_mesh.__init__(self, coordinates, triangles,
    … …  
    98102
    99103        #Allocate space for geometric quantities
    100         self.centroid_coordinates = zeros((N, 2), Float)
    101 
    102         self.radii = zeros(N, Float)
    103 
    104         self.neighbours = zeros((N, 3), Int)
    105         self.neighbour_edges = zeros((N, 3), Int)
    106         self.number_of_boundaries = zeros(N, Int)
    107         self.surrogate_neighbours = zeros((N, 3), Int)
     104        self.centroid_coordinates = zeros((N, 2), float)
     105
     106        self.radii = zeros(N, float)
     107
     108        self.neighbours = zeros((N, 3), int)
     109        self.neighbour_edges = zeros((N, 3), int)
     110        self.number_of_boundaries = zeros(N, int)
     111        self.surrogate_neighbours = zeros((N, 3), int)
    108112
    109113        #Get x,y coordinates for all triangles and store
    … …  
    388392            self.element_tag is defined
    389393        """
    390         from Numeric import array, Int
     394#        from numpy.oldnumeric import array, Int
     395        from numpy import array, int
    391396
    392397        if tagged_elements is None:
    … …  
    395400            #Check that all keys in given boundary exist
    396401            for tag in tagged_elements.keys():
    397                 tagged_elements[tag] = array(tagged_elements[tag]).astype(Int)
     402                tagged_elements[tag] = array(tagged_elements[tag]).astype(int)
    398403
    399404                msg = 'Not all elements exist. '
    … …  
    463468        """
    464469       
    465         from Numeric import allclose, sqrt, array, minimum, maximum
     470#        from numpy.oldnumeric import allclose, sqrt, array, minimum, maximum
     471        from numpy import allclose, sqrt, array, minimum, maximum
    466472        from anuga.utilities.numerical_tools import angle, ensure_numeric     
    467473
    … …  
    629635        from anuga.utilities.numerical_tools import anglediff
    630636
    631         from Numeric import sort, allclose
     637#        from numpy.oldnumeric import sort, allclose
     638        from numpy import sort, allclose
    632639
    633640        N = len(self)
    … …  
    805812        """
    806813
    807         from Numeric import arange
     814#        from numpy.oldnumeric import arange
     815        from numpy import arange
    808816        from anuga.utilities.numerical_tools import histogram, create_bins
    809817

  • anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/pmesh2domain.py

    r4902 r5899  
     1## Automatically adapted for numpy.oldnumeric Oct 28, 2008 by alter_code1.py
     2
    13"""Class pmesh2domain - Converting .tsh files to doamains
    24
    … …  
    161163    """
    162164
    163     from Numeric import transpose
     165#    from numpy.oldnumeric import transpose
     166    from numpy import transpose
    164167    from load_mesh.loadASCII import import_mesh_file
    165168
    … …  
    176179    geo_reference  = mesh_dict['geo_reference']
    177180    if point_atts != None:
    178         for quantity, value_vector in map (None, point_titles, point_atts):
     181        print 'type(point_atts)=%s' % type(point_atts)
     182        for quantity, value_vector in map(None, point_titles, point_atts):
    179183            vertex_quantity_dict[quantity] = value_vector
    180184    tag_dict = pmesh_dict_to_tag_dict(mesh_dict)

  • anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/quantity.py

    r5866 r5899  
     1## Automatically adapted for numpy.oldnumeric Oct 28, 2008 by alter_code1.py
     2
    13"""Class Quantity - Implements values at each triangular element
    24
    … …  
    1517"""
    1618
    17 from Numeric import array, zeros, Float, less, concatenate, NewAxis,\
    18      argmax, argmin, allclose, take, reshape, alltrue
     19#from numpy.oldnumeric import array, zeros, Float, less, concatenate, NewAxis,\
     20#     argmax, argmin, allclose, take, reshape, alltrue, Int
     21from numpy import array, zeros, float, less, concatenate, \
     22     argmax, argmin, allclose, take, reshape, alltrue, int
    1923
    2024from anuga.utilities.numerical_tools import ensure_numeric, is_scalar
    … …  
    3842        if vertex_values is None:
    3943            N = len(domain) # number_of_elements
    40             self.vertex_values = zeros((N, 3), Float)
     44            self.vertex_values = zeros((N, 3), float)
    4145        else:
    42             self.vertex_values = array(vertex_values).astype(Float)
     46            self.vertex_values = array(vertex_values).astype(float)
    4347
    4448            N, V = self.vertex_values.shape
    … …  
    5761
    5862        # Allocate space for other quantities
    59         self.centroid_values = zeros(N, Float)
    60         self.edge_values = zeros((N, 3), Float)
     63        self.centroid_values = zeros(N, float)
     64        self.edge_values = zeros((N, 3), float)
    6165
    6266        # Allocate space for Gradient
    63         self.x_gradient = zeros(N, Float)
    64         self.y_gradient = zeros(N, Float)
     67        self.x_gradient = zeros(N, float)
     68        self.y_gradient = zeros(N, float)
    6569
    6670        # Allocate space for Limiter Phi
    67         self.phi = zeros(N, Float)       
     71        self.phi = zeros(N, float)       
    6872
    6973        # Intialise centroid and edge_values
    … …  
    7276        # Allocate space for boundary values
    7377        L = len(domain.boundary)
    74         self.boundary_values = zeros(L, Float)
     78        self.boundary_values = zeros(L, float)
    7579
    7680        # Allocate space for updates of conserved quantities by
    … …  
    7882
    7983        # Allocate space for update fields
    80         self.explicit_update = zeros(N, Float )
    81         self.semi_implicit_update = zeros(N, Float )
    82         self.centroid_backup_values = zeros(N, Float)
     84        self.explicit_update = zeros(N, float )
     85        self.semi_implicit_update = zeros(N, float )
     86        self.centroid_backup_values = zeros(N, float)
    8387
    8488        self.set_beta(1.0)
    … …  
    382386        from anuga.geospatial_data.geospatial_data import Geospatial_data
    383387        from types import FloatType, IntType, LongType, ListType, NoneType
    384         from Numeric import ArrayType
     388##NumPy        from numpy.oldnumeric import ArrayType
     389        from numpy import ndarray, float
    385390
    386391        # Treat special case: Polygon situation
    … …  
    448453
    449454        msg = 'Indices must be a list or None'
    450         assert type(indices) in [ListType, NoneType, ArrayType], msg
     455##NumPy        assert type(indices) in [ListType, NoneType, ArrayType], msg
     456        assert type(indices) in [ListType, NoneType, ndarray], msg
    451457
    452458
    … …  
    458464                self.set_values_from_constant(numeric,
    459465                                              location, indices, verbose)
    460             elif type(numeric) in [ArrayType, ListType]:
     466##NumPy            elif type(numeric) in [ArrayType, ListType]:
     467            elif type(numeric) in [ndarray, ListType]:
    461468                self.set_values_from_array(numeric,
    462469                                           location, indices, verbose)
    … …  
    474481                                                     use_cache=use_cache)
    475482            else:
    476                 msg = 'Illegal type for argument numeric: %s' %str(numeric)
    477                 raise msg
     483                msg = 'Illegal type for argument numeric: %s' % str(numeric)
     484                raise TypeError, msg
    478485
    479486        elif quantity is not None:
    … …  
    610617        """
    611618
    612         from Numeric import array, Float, Int, allclose
    613 
    614         values = array(values).astype(Float)
     619#        from numpy.oldnumeric import array, Float, Int, allclose
     620        from numpy import array, float, int, allclose
     621
     622        values = array(values).astype(float)
    615623
    616624        if indices is not None:
    617             indices = array(indices).astype(Int)
     625            indices = array(indices).astype(int)
    618626            msg = 'Number of values must match number of indices:'
    619627            msg += ' You specified %d values and %d indices'\
    … …  
    643651                   'Values array must be 1d'
    644652
    645             self.set_vertex_values(values.flat, indices=indices)
     653            self.set_vertex_values(values.ravel(), indices=indices)
    646654           
    647655        else:
    … …  
    676684        A = q.vertex_values
    677685
    678         from Numeric import allclose
     686#        from numpy.oldnumeric import allclose
     687        from numpy import allclose
    679688        msg = 'Quantities are defined on different meshes. '+\
    680689              'This might be a case for implementing interpolation '+\
    … …  
    717726               
    718727            V = take(self.domain.get_centroid_coordinates(), indices)
     728            print 'V=%s' % str(V)
    719729            self.set_values(f(V[:,0], V[:,1]),
    720730                            location=location,
    … …  
    780790
    781791
    782         points = ensure_numeric(points, Float)
    783         values = ensure_numeric(values, Float)
     792        points = ensure_numeric(points, float)
     793        values = ensure_numeric(values, float)
    784794
    785795        if location != 'vertices':
    … …  
    11181128        """
    11191129       
    1120         from Numeric import take
     1130#        from numpy.oldnumeric import take
     1131        from numpy import take
    11211132
    11221133        # FIXME (Ole): I reckon we should have the option of passing a
    … …  
    11431154            raise msg
    11441155
    1145         import types, Numeric
    1146         assert type(indices) in [types.ListType, types.NoneType,
    1147                                  Numeric.ArrayType],\
     1156#        import types, numpy.oldnumeric as Numeric
     1157        import types
     1158        from numpy import ndarray
     1159       
     1160        assert type(indices) in [types.ListType, types.NoneType, ndarray], \
    11481161                                 'Indices must be a list or None'
    11491162
    … …  
    11961209        """
    11971210
    1198         from Numeric import array, Float
     1211#        from numpy.oldnumeric import array, Float
     1212        from numpy import array, float
    11991213
    12001214        # Assert that A can be converted to a Numeric array of appropriate dim
    1201         A = ensure_numeric(A, Float)
     1215        A = ensure_numeric(A, float)
    12021216
    12031217        # print 'SHAPE A', A.shape
    … …  
    12731287        """
    12741288
    1275         from Numeric import concatenate, zeros, Float, Int, array, reshape
     1289#        from numpy.oldnumeric import concatenate, zeros, Float, Int, array, reshape
     1290        from numpy import concatenate, zeros, float, int, array, reshape
    12761291
    12771292
    … …  
    12841299
    12851300        if precision is None:
    1286             precision = Float
     1301            precision = float
    12871302           
    12881303
    … …  
    12931308            V = self.domain.get_triangles()
    12941309            N = self.domain.number_of_full_nodes # Ignore ghost nodes if any
    1295             A = zeros(N, Float)
     1310            A = zeros(N, float)
    12961311            points = self.domain.get_nodes()           
    12971312           
    … …  
    13431358            V = self.domain.get_disconnected_triangles()
    13441359            points = self.domain.get_vertex_coordinates()
    1345             A = self.vertex_values.flat.astype(precision)
     1360            A = self.vertex_values.ravel().astype(precision)
    13461361
    13471362

  • anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/quantity_ext.c

    r5306 r5899  
    1111
    1212#include "Python.h"
    13 #include "Numeric/arrayobject.h"
     13#include "numpy/arrayobject.h"
    1414#include "math.h"
    1515

  • anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/region.py

    r5208 r5899  
     1## Automatically adapted for numpy.oldnumeric Oct 28, 2008 by alter_code1.py
     2
    13"""region.py - Classes for implementing region conditions
    24
    … …  
    810# FIXME (DSG-DSG) add better comments
    911
    10 from Numeric import average
     12from numpy.oldnumeric import average
    1113class Region:
    1214    """Base class for modifying quantities based on a region.

  • anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/show_balanced_limiters.py

    r4204 r5899  
     1## Automatically adapted for numpy.oldnumeric Oct 28, 2008 by alter_code1.py
     2
    13"""Example of shallow water wave equation.
    24
    … …  
    1820
    1921from mesh_factory import rectangular
    20 from Numeric import array
     22from numpy.oldnumeric import array
    2123   
    2224
    … …  
    7880domain.set_quantity('stage', Z)
    7981
    80 from Numeric import allclose
     82from numpy.oldnumeric import allclose
    8183
    8284#Evolve

  • anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/test_domain.py

    r4932 r5899  
     1## Automatically adapted for numpy.oldnumeric Oct 28, 2008 by alter_code1.py
     2
    13#!/usr/bin/env python
    24
    … …  
    68from domain import *
    79from anuga.config import epsilon
    8 from Numeric import allclose, array, ones, Float
     10##from numpy.oldnumeric import allclose, array, ones, Float, alltrue
     11from numpy import allclose, array, ones, float, alltrue
    912
    1013
    … …  
    6467
    6568
    66         assert domain.get_conserved_quantities(0, edge=1) == 0.
     69        assert alltrue(domain.get_conserved_quantities(0, edge=1) == 0.)
    6770
    6871
    … …  
    536539
    537540        sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])
    538         denom = ones(4, Float)-domain.timestep*sem
     541        denom = ones(4, float)-domain.timestep*sem
    539542
    540543#        x = array([1, 2, 3, 4]) + array( [.4,.3,.2,.1] )
    … …  
    610613        from mesh_factory import rectangular
    611614        from shallow_water import Domain
    612         from Numeric import zeros, Float
     615##        from numpy.oldnumeric import zeros, Float
     616        from numpy import zeros, float
    613617
    614618        #Create basic mesh
    … …  
    627631        from mesh_factory import rectangular
    628632        from shallow_water import Domain
    629         from Numeric import zeros, Float
     633##        from numpy.oldnumeric import zeros, Float
     634        from numpy import zeros, float
    630635
    631636        #Create basic mesh
    … …  
    670675        from mesh_factory import rectangular
    671676        from shallow_water import Domain
    672         from Numeric import zeros, Float
     677##        from numpy.oldnumeric import zeros, Float
     678        from numpy import zeros, float
    673679
    674680        #Create basic mesh

  • anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/test_ermapper.py

    r1813 r5899  
     1## Automatically adapted for numpy.oldnumeric Oct 28, 2008 by alter_code1.py
     2
    13#!/usr/bin/env python
    24
    … …  
    57
    68import ermapper_grids
    7 import Numeric
     9##import numpy.oldnumeric as Numeric
     10import numpy
    811from os import remove
    912
    … …  
    1922        data_filename = 'test_write_ermapper_grid'
    2023       
    21         original_grid = Numeric.array([[0.0, 0.1, 1.0], [2.0, 3.0, 4.0]])
     24        original_grid = numpy.array([[0.0, 0.1, 1.0], [2.0, 3.0, 4.0]])
    2225
    2326        # Check that the function works when passing the filename without
    … …  
    2528        ermapper_grids.write_ermapper_grid(data_filename, original_grid)
    2629        new_grid = ermapper_grids.read_ermapper_grid(data_filename)
    27         assert Numeric.allclose(original_grid, new_grid)
     30        assert numpy.allclose(original_grid, new_grid)
    2831
    2932        # Check that the function works when passing the filename with
    … …  
    3134        ermapper_grids.write_ermapper_grid(header_filename, original_grid)
    3235        new_grid = ermapper_grids.read_ermapper_grid(header_filename)
    33         assert Numeric.allclose(original_grid, new_grid)
     36        assert numpy.allclose(original_grid, new_grid)
    3437
    3538        # Clean up created files
    … …  
    4043        # Setup test data
    4144        filename = 'test_write_ermapper_grid'
    42         original_grid = Numeric.array([0.0, 0.1, 1.0, 2.0, 3.0, 4.0])
     45        original_grid = numpy.array([0.0, 0.1, 1.0, 2.0, 3.0, 4.0])
    4346
    4447        # Write test data
    45         ermapper_grids.write_ermapper_data(original_grid, filename, Numeric.Float64)
     48        ermapper_grids.write_ermapper_data(original_grid, filename, numpy.float64)
    4649
    4750        # Read in the test data
    48         new_grid = ermapper_grids.read_ermapper_data(filename, Numeric.Float64)
     51        new_grid = ermapper_grids.read_ermapper_data(filename, numpy.float64)
    4952
    5053        # Check that the test data that has been read in matches the original data
    51         assert Numeric.allclose(original_grid, new_grid)
     54        assert numpy.allclose(original_grid, new_grid)
    5255
    5356        # Clean up created files
    … …  
    5760        # Setup test data
    5861        filename = 'test_write_ermapper_grid'
    59         original_grid = Numeric.array([0.0, 0.1, 1.0, 2.0, 3.0, 4.0])
     62        original_grid = numpy.array([0.0, 0.1, 1.0, 2.0, 3.0, 4.0])
    6063
    6164        # Write test data
    … …  
    6669
    6770        # Check that the test data that has been read in matches the original data
    68         assert Numeric.allclose(original_grid, new_grid)
     71        assert numpy.allclose(original_grid, new_grid)
    6972
    7073        # Clean up created files
    … …  
    7578
    7679        # setup test data
    77         original_grid = Numeric.array([[0, 1, 2, 3], [4, 5, 6, 7], [8, 9, 10, 11]])
     80        original_grid = numpy.array([[0, 1, 2, 3], [4, 5, 6, 7], [8, 9, 10, 11]])
    7881        # Write test data
    7982        ermapper_grids.write_ermapper_data(original_grid, data_filename)

  • anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/test_general_mesh.py

    r5843 r5899  
     1## Automatically adapted for numpy.oldnumeric Oct 28, 2008 by alter_code1.py
     2
    13#!/usr/bin/env python
    24
    … …  
    79
    810from anuga.config import epsilon
    9 from Numeric import allclose, array, ones, Float
     11##from numpy.oldnumeric import allclose, array, ones, Float, alltrue
     12from numpy import allclose, array, ones, float, alltrue
    1013from general_mesh import General_mesh
    1114from anuga.coordinate_transforms.geo_reference import Geo_reference
    … …  
    2326    def test_get_vertex_coordinates(self):
    2427        from mesh_factory import rectangular
    25         from Numeric import zeros, Float
     28##        from numpy.oldnumeric import zeros, Float
     29        from numpy import zeros, float
    2630
    2731        #Create basic mesh
    … …  
    8690        """
    8791        from mesh_factory import rectangular
    88         from Numeric import zeros, Float
     92##        from numpy.oldnumeric import zeros, Float
     93        from numpy import zeros, float
    8994
    9095        #Create basic mesh
    … …  
    114119        """
    115120        from mesh_factory import rectangular
    116         from Numeric import zeros, Float
     121##        from numpy.oldnumeric import zeros, Float
     122        from numpy import zeros, float
    117123
    118124        #Create basic mesh
    … …  
    130136        """
    131137        from mesh_factory import rectangular
    132         from Numeric import zeros, Float, array
     138##        from numpy.oldnumeric import zeros, Float, array
     139        from numpy import zeros, float, array
    133140
    134141        a = [0.0, 0.0]
    … …  
    204211        # One node
    205212        L = domain.get_triangles_and_vertices_per_node(node=2)
     213        print 'L=%s' % str(L)
    206214        assert allclose(L[0], [0, 2])
    207215        assert allclose(L[1], [1, 1])
    … …  
    224232        from mesh_factory import rectangular
    225233        from shallow_water import Domain
    226         from Numeric import zeros, Float
     234##        from numpy.oldnumeric import zeros, Float
     235        from numpy import zeros, float
    227236
    228237        #Create basic mesh
    … …  
    239248        from mesh_factory import rectangular
    240249        from shallow_water import Domain
    241         from Numeric import zeros, Float
     250##        from numpy.oldnumeric import zeros, Float
     251        from numpy import zeros, float
    242252
    243253        #Create basic mesh
    … …  
    283293                       geo_reference = geo)
    284294        node = domain.get_node(2)       
    285         self.assertEqual(c, node)
     295        self.assertTrue(alltrue(c == node))
    286296       
    287297        node = domain.get_node(2, absolute=True)     
    288         self.assertEqual(nodes_absolute[2], node)
     298        self.assertTrue(alltrue(nodes_absolute[2] == node))
    289299       
    290300        node = domain.get_node(2, absolute=True)     
    291         self.assertEqual(nodes_absolute[2], node)
     301        self.assertTrue(alltrue(nodes_absolute[2] == node))
    292302       
    293303
    … …  
    333343if __name__ == "__main__":
    334344    suite = unittest.makeSuite(Test_General_Mesh,'test')
    335     #suite = unittest.makeSuite(Test_General_Mesh,'test_get_node')   
    336345    runner = unittest.TextTestRunner()
    337346    runner.run(suite)

  • anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/test_generic_boundary_conditions.py

    r5666 r5899  
     1## Automatically adapted for numpy.oldnumeric Oct 28, 2008 by alter_code1.py
     2
    13#!/usr/bin/env python
    24
    … …  
    68from generic_boundary_conditions import *
    79from anuga.config import epsilon
    8 from Numeric import allclose, array
     10#from numpy.oldnumeric import allclose, array
     11from numpy import allclose, array
    912
    1013

  • anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/test_ghost.py

    r3514 r5899  
     1## Automatically adapted for numpy.oldnumeric Oct 28, 2008 by alter_code1.py
     2
    13#!/usr/bin/env python
    24
    … …  
    68from domain import *
    79from anuga.config import epsilon
    8 from Numeric import allclose, array, ones, Float
     10from numpy.oldnumeric import allclose, array, ones, Float, alltrue
    911
    1012
    … …  
    4345
    4446
    45         assert domain.get_conserved_quantities(0, edge=1) == 0.
     47        assert alltrue(domain.get_conserved_quantities(0, edge=1) == 0.)
    4648
    4749

  • anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/test_neighbour_mesh.py

    r5288 r5899  
     1## Automatically adapted for numpy.oldnumeric Oct 28, 2008 by alter_code1.py
     2
    13#!/usr/bin/env python
    24
    … …  
    1315from mesh_factory import rectangular
    1416from anuga.config import epsilon
    15 from Numeric import allclose, array, Int
     17#from numpy.oldnumeric import allclose, array, Int
     18from numpy import allclose, array, int
    1619
    1720from anuga.coordinate_transforms.geo_reference import Geo_reference
    … …  
    698701        """
    699702        from mesh_factory import rectangular
    700         from Numeric import zeros, Float
     703#        from numpy.oldnumeric import zeros, Float
     704        from numpy import zeros, float
    701705
    702706        #Create basic mesh
    … …  
    717721        from mesh_factory import rectangular
    718722        #from mesh import Mesh
    719         from Numeric import zeros, Float
     723#        from numpy.oldnumeric import zeros, Float
     724        from numpy import zeros, float
    720725
    721726        #Create basic mesh
    … …  
    736741
    737742    def test_boundary_polygon_II(self):
    738         from Numeric import zeros, Float
     743#        from numpy.oldnumeric import zeros, Float
     744        from numpy import zeros, float
    739745       
    740746
    … …  
    775781        """
    776782
    777         from Numeric import zeros, Float
     783#        from numpy.oldnumeric import zeros, Float
     784        from numpy import zeros, float
    778785
    779786
    … …  
    816823        """
    817824
    818         from Numeric import zeros, Float
     825#        from numpy.oldnumeric import zeros, Float
     826        from numpy import zeros, float
    819827
    820828
    … …  
    864872        """
    865873
    866         from Numeric import zeros, Float
     874#        from numpy.oldnumeric import zeros, Float
     875        from numpy import zeros, float
    867876        from mesh_factory import rectangular       
    868877
    … …  
    907916       
    908917        """
    909         from Numeric import zeros, Float
     918#        from numpy.oldnumeric import zeros, Float
     919        from numpy import zeros, float
    910920       
    911921
    … …  
    9981008                  [  52341.70703125,  38563.39453125]]
    9991009
    1000         ##points = ensure_numeric(points, Int)/1000  # Simplify for ease of interpretation       
     1010        ##points = ensure_numeric(points, int)/1000  # Simplify for ease of interpretation       
    10011011
    10021012        triangles = [[19, 0,15],
    … …  
    11011111                  [  35406.3359375 ,  79332.9140625 ]]
    11021112
    1103         scaled_points = ensure_numeric(points, Int)/1000  # Simplify for ease of interpretation
     1113        scaled_points = ensure_numeric(points, int)/1000  # Simplify for ease of interpretation
    11041114
    11051115        triangles = [[ 0, 1, 2],

  • anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/test_pmesh2domain.py

    r3563 r5899  
     1## Automatically adapted for numpy.oldnumeric Oct 28, 2008 by alter_code1.py
     2
    13#!/usr/bin/env python
    24#
    … …  
    46import unittest
    57
    6 from Numeric import allclose, array
     8#from numpy.oldnumeric import allclose, array
     9from numpy import allclose, array
    710
    811

  • anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/test_quantity.py

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