Ignore:
Timestamp:
Feb 10, 2009, 11:11:04 AM (16 years ago)
Author:
rwilson
Message:

Initial commit of numpy changes. Still a long way to go.

Location:
branches/numpy
Files:
21 edited
1 copied

Legend:

Unmodified
Added
Removed

  • branches/numpy/anuga/abstract_2d_finite_volumes/domain.py

    r6246 r6304  
    3131from quantity import Quantity
    3232
    33 import Numeric as num
     33import numpy as num
    3434
    3535
    … …  
    185185            buffer_shape = self.full_send_dict[key][0].shape[0]
    186186            self.full_send_dict[key].append(num.zeros((buffer_shape, self.nsys),
    187                                                       num.Float))
     187                                                      num.float))
    188188
    189189        for key in self.ghost_recv_dict:
    190190            buffer_shape = self.ghost_recv_dict[key][0].shape[0]
    191191            self.ghost_recv_dict[key].append(num.zeros((buffer_shape, self.nsys),
    192                                              num.Float))
     192                                             num.float))
    193193
    194194        # Setup cell full flag
    … …  
    197197        N = len(self) #number_of_elements
    198198        self.number_of_elements = N
    199         self.tri_full_flag = num.ones(N, num.Int)
     199        self.tri_full_flag = num.ones(N, num.int)
    200200        for i in self.ghost_recv_dict.keys():
    201201            for id in self.ghost_recv_dict[i][0]:
    … …  
    258258        # (boolean) array, to be used during the flux calculation.
    259259        N = len(self) # Number_of_triangles
    260         self.already_computed_flux = num.zeros((N, 3), num.Int)
     260        self.already_computed_flux = num.zeros((N, 3), num.int)
    261261
    262262        # Storage for maximal speeds computed for each triangle by
    263263        # compute_fluxes.
    264264        # This is used for diagnostics only (reset at every yieldstep)
    265         self.max_speed = num.zeros(N, num.Float)
     265        self.max_speed = num.zeros(N, num.float)
    266266
    267267        if mesh_filename is not None:
    … …  
    375375            raise Exception, msg
    376376
    377         q = num.zeros(len(self.conserved_quantities), num.Float)
     377        q = num.zeros(len(self.conserved_quantities), num.float)
    378378
    379379        for i, name in enumerate(self.conserved_quantities):
    … …  
    447447
    448448        name:  Name of quantity
    449         value: Compatible list, Numeric array, const or function (see below)
     449        value: Compatible list, numeric array, const or function (see below)
    450450
    451451        The values will be stored in elements following their internal ordering.
    … …  
    879879
    880880            msg += '------------------------------------------------\n'
    881             msg += '  Speeds in [%f, %f]\n' % (min(self.max_speed),
    882                                                max(self.max_speed))
     881            msg += '  Speeds in [%f, %f]\n' % (min(self.max_speed.flat),
     882                                               max(self.max_speed.flat))
    883883            msg += '  Histogram:\n'
    884884
    … …  
    892892                else:
    893893                    # Closed upper interval
    894                     hi = max(self.max_speed)
     894                    hi = max(self.max_speed.flat)
    895895                    msg += '    [%f, %f]: %d\n' % (lo, hi, count)
    896896
    897             N = len(self.max_speed)
     897            N = len(self.max_speed.flat)
    898898            if N > 10:
    899899                msg += '  Percentiles (10%):\n'
    … …  
    13571357                self.number_of_steps = 0
    13581358                self.number_of_first_order_steps = 0
    1359                 self.max_speed = num.zeros(N, num.Float)
     1359                self.max_speed = num.zeros(N, num.float)
    13601360
    13611361    ##
    … …  
    16261626        # momentum
    16271627        if self.protect_against_isolated_degenerate_timesteps is True and \
    1628                self.max_speed > 10.0: # FIXME (Ole): Make this configurable
     1628               num.max(self.max_speed) > 10.0: # FIXME (Ole): Make this configurable
    16291629
    16301630            # Setup 10 bins for speed histogram
    … …  
    16411641
    16421642                    # Find triangles in last bin
    1643                     # FIXME - speed up using Numeric
     1643                    # FIXME - speed up using numeric package
    16441644                    d = 0
    16451645                    for i in range(self.number_of_full_triangles):

  • branches/numpy/anuga/abstract_2d_finite_volumes/ermapper_grids.py

    r6161 r6304  
    22
    33# from os import open, write, read
    4 import Numeric as num
    5 
    6 celltype_map = {'IEEE4ByteReal': num.Float32, 'IEEE8ByteReal': num.Float64}
     4import numpy as num
     5
     6celltype_map = {'IEEE4ByteReal': num.float32, 'IEEE8ByteReal': num.float64}
    77
    88
    … …  
    1111    write_ermapper_grid(ofile, data, header = {}):
    1212   
    13     Function to write a 2D Numeric array to an ERMapper grid.  There are a series of conventions adopted within
     13    Function to write a 2D numeric array to an ERMapper grid.  There are a series of conventions adopted within
    1414    this code, specifically:
    1515    1)  The registration coordinate for the data is the SW (or lower-left) corner of the data
    1616    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)
     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)
    1818        where N is the last line and M is the last column
    1919    4)  There has been no testng of the use of a rotated grid.  Best to keep data in an NS orientation
    … …  
    163163    return header                     
    164164
    165 def write_ermapper_data(grid, ofile, data_format = num.Float32):
     165def write_ermapper_data(grid, ofile, data_format=num.float32):
    166166
    167167
    … …  
    193193
    194194
    195 def read_ermapper_data(ifile, data_format = num.Float32):
     195def read_ermapper_data(ifile, data_format = num.float32):
    196196    # open input file in a binary format and read the input string
    197197    fid = open(ifile,'rb')
    … …  
    199199    fid.close()
    200200
    201     # convert input string to required format (Note default format is num.Float32)
     201    # convert input string to required format (Note default format is num.float32)
    202202    grid_as_float = num.fromstring(input_string,data_format)
    203203    return grid_as_float

  • branches/numpy/anuga/abstract_2d_finite_volumes/general_mesh.py

    r6191 r6304  
    1 import Numeric as num
     1import numpy as num
    22
    33from anuga.coordinate_transforms.geo_reference import Geo_reference
    … …  
    1010    which is represented as a coordinate set (x,y).
    1111    Vertices from different triangles can point to the same node.
    12     The nodes are implemented as an Nx2 Numeric array containing the
     12    The nodes are implemented as an Nx2 numeric array containing the
    1313    x and y coordinates.
    1414
    … …  
    1919    where
    2020
    21       nodes is either a list of 2-tuples or an Nx2 Numeric array of
     21      nodes is either a list of 2-tuples or an Nx2 numeric array of
    2222      floats representing all x, y coordinates in the mesh.
    2323
    24       triangles is either a list of 3-tuples or an Mx3 Numeric array of
     24      triangles is either a list of 3-tuples or an Mx3 numeric array of
    2525      integers representing indices of all vertices in the mesh.
    2626      Each vertex is identified by its index i in [0, N-1].
    … …  
    6868       
    6969          nodes: x,y coordinates represented as a sequence of 2-tuples or
    70                  a Nx2 Numeric array of floats.
     70                 a Nx2 numeric array of floats.
    7171                 
    72           triangles: sequence of 3-tuples or Mx3 Numeric array of
     72          triangles: sequence of 3-tuples or Mx3 numeric array of
    7373                     non-negative integers representing indices into
    7474                     the nodes array.
    … …  
    8888        if verbose: print 'General_mesh: Building basic mesh structure in ANUGA domain'
    8989
    90         self.triangles = num.array(triangles, num.Int)
    91         self.nodes = num.array(nodes, num.Float)
     90        self.triangles = num.array(triangles, num.int)
     91        self.nodes = num.array(nodes, num.float)
    9292
    9393
    … …  
    125125
    126126        # Input checks
    127         msg = 'Triangles must an Mx3 Numeric array or a sequence of 3-tuples. '
     127        msg = 'Triangles must an Mx3 numeric array or a sequence of 3-tuples. '
    128128        msg += 'The supplied array has the shape: %s'\
    129129               %str(self.triangles.shape)
    130130        assert len(self.triangles.shape) == 2, msg
    131131
    132         msg = 'Nodes must an Nx2 Numeric array or a sequence of 2-tuples'
     132        msg = 'Nodes must an Nx2 numeric array or a sequence of 2-tuples'
    133133        msg += 'The supplied array has the shape: %s'\
    134134               %str(self.nodes.shape)
    … …  
    144144                      max(self.nodes[:,0]), max(self.nodes[:,1]) ]
    145145
    146         self.xy_extent = num.array(xy_extent, num.Float)
     146        self.xy_extent = num.array(xy_extent, num.float)
    147147
    148148
    149149        # Allocate space for geometric quantities
    150         self.normals = num.zeros((N, 6), num.Float)
    151         self.areas = num.zeros(N, num.Float)
    152         self.edgelengths = num.zeros((N, 3), num.Float)
     150        self.normals = num.zeros((N, 6), num.float)
     151        self.areas = num.zeros(N, num.float)
     152        self.edgelengths = num.zeros((N, 3), num.float)
    153153
    154154        # Get x,y coordinates for all triangles and store
    … …  
    185185            #   - Stored as six floats n0x,n0y,n1x,n1y,n2x,n2y per triangle
    186186
    187             n0 = num.array([x2 - x1, y2 - y1], num.Float)
     187            n0 = num.array([x2 - x1, y2 - y1], num.float)
    188188            l0 = num.sqrt(num.sum(n0**2))
    189189
    190             n1 = num.array([x0 - x2, y0 - y2], num.Float)
     190            n1 = num.array([x0 - x2, y0 - y2], num.float)
    191191            l1 = num.sqrt(num.sum(n1**2))
    192192
    193             n2 = num.array([x1 - x0, y1 - y0], num.Float)
     193            n2 = num.array([x1 - x0, y1 - y0], num.float)
    194194            l2 = num.sqrt(num.sum(n2**2))
    195195
    … …  
    275275            if not self.geo_reference.is_absolute():
    276276                return V + num.array([self.geo_reference.get_xllcorner(),
    277                                       self.geo_reference.get_yllcorner()], num.Float)
     277                                      self.geo_reference.get_yllcorner()], num.float)
    278278            else:
    279279                return V
    … …  
    317317            if absolute is True and not self.geo_reference.is_absolute():
    318318                offset=num.array([self.geo_reference.get_xllcorner(),
    319                                   self.geo_reference.get_yllcorner()], num.Float)
     319                                  self.geo_reference.get_yllcorner()], num.float)
    320320                return num.array([V[i3,:]+offset,
    321321                                  V[i3+1,:]+offset,
    322                                   V[i3+2,:]+offset], num.Float)
     322                                  V[i3+2,:]+offset], num.float)
    323323            else:
    324                 return num.array([V[i3,:], V[i3+1,:], V[i3+2,:]], num.Float)
     324                return num.array([V[i3,:], V[i3+1,:], V[i3+2,:]], num.float)
    325325               
    326326
    … …  
    349349
    350350        M = self.number_of_triangles
    351         vertex_coordinates = num.zeros((3*M, 2), num.Float)
     351        vertex_coordinates = num.zeros((3*M, 2), num.float)
    352352
    353353        for i in range(M):
    … …  
    409409        K = 3*M       # Total number of unique vertices
    410410       
    411         T = num.reshape(num.arange(K, typecode=num.Int), (M,3))
     411        T = num.reshape(num.arange(K, dtype=num.int), (M,3))
    412412       
    413413        return T     
    … …  
    416416
    417417    def get_unique_vertices(self,  indices=None):
    418         """FIXME(Ole): This function needs a docstring
    419         """
     418        """FIXME(Ole): This function needs a docstring"""
     419
    420420        triangles = self.get_triangles(indices=indices)
    421421        unique_verts = {}
    422422        for triangle in triangles:
    423            unique_verts[triangle[0]] = 0
    424            unique_verts[triangle[1]] = 0
    425            unique_verts[triangle[2]] = 0
     423            #print 'triangle(%s)=%s' % (type(triangle), str(triangle))
     424            unique_verts[triangle[0]] = 0
     425            unique_verts[triangle[1]] = 0
     426            unique_verts[triangle[2]] = 0
    426427        return unique_verts.keys()
    427428
    … …  
    452453                triangle_list.append( (volume_id, vertex_id) )
    453454
    454             triangle_list = num.array(triangle_list, num.Int)    #array default#
     455            triangle_list = num.array(triangle_list, num.int)    #array default#
    455456        else:
    456457            # Get info for all nodes recursively.
    … …  
    529530
    530531        # Count number of triangles per node
    531         number_of_triangles_per_node = num.zeros(self.number_of_full_nodes)
     532        number_of_triangles_per_node = num.zeros(self.number_of_full_nodes,
     533                                                 num.int)       #array default#
    532534        for volume_id, triangle in enumerate(self.get_triangles()):
    533535            for vertex_id in triangle:
    … …  
    536538        # Allocate space for inverted structure
    537539        number_of_entries = num.sum(number_of_triangles_per_node)
    538         vertex_value_indices = num.zeros(number_of_entries)
     540        vertex_value_indices = num.zeros(number_of_entries, num.int) #array default#
    539541
    540542        # Register (triangle, vertex) indices for each node

  • branches/numpy/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py

    r6174 r6304  
    88from anuga.fit_interpolate.interpolate import Modeltime_too_early
    99
    10 import Numeric as num
     10import numpy as num
    1111
    1212
    … …  
    6969            raise Exception, msg
    7070
    71         self.conserved_quantities=num.array(conserved_quantities, num.Float)
     71        self.conserved_quantities=num.array(conserved_quantities, num.float)
    7272
    7373    def __repr__(self):
    … …  
    9999
    100100        try:
    101             q = num.array(q, num.Float)
     101            q = num.array(q, num.float)
    102102        except:
    103103            msg = 'Return value from time boundary function could '
    104             msg += 'not be converted into a Numeric array of floats.\n'
     104            msg += 'not be converted into a numeric array of floats.\n'
    105105            msg += 'Specified function should return either list or array.\n'
    106106            msg += 'I got %s' %str(q)
    … …  
    180180
    181181        if verbose: print 'Find midpoint coordinates of entire boundary'
    182         self.midpoint_coordinates = num.zeros((len(domain.boundary), 2), num.Float)
     182        self.midpoint_coordinates = num.zeros((len(domain.boundary), 2), num.float)
    183183        boundary_keys = domain.boundary.keys()
    184184
    … …  
    200200           
    201201            # Compute midpoints
    202             if edge_id == 0: m = num.array([(x1 + x2)/2, (y1 + y2)/2], num.Float)
    203             if edge_id == 1: m = num.array([(x0 + x2)/2, (y0 + y2)/2], num.Float)
    204             if edge_id == 2: m = num.array([(x1 + x0)/2, (y1 + y0)/2], num.Float)
     202            if edge_id == 0: m = num.array([(x1 + x2)/2, (y1 + y2)/2], num.float)
     203            if edge_id == 1: m = num.array([(x0 + x2)/2, (y0 + y2)/2], num.float)
     204            if edge_id == 2: m = num.array([(x1 + x0)/2, (y1 + y0)/2], num.float)
    205205
    206206            # Convert to absolute UTM coordinates
    … …  
    309309                        # See http://docs.python.org/lib/warning-filter.html
    310310                        self.default_boundary_invoked = True
    311                    
    312 
    313             if res == NAN:
     311           
     312            if num.any(res == NAN):
    314313                x,y=self.midpoint_coordinates[i,:]
    315314                msg = 'NAN value found in file_boundary at '

  • branches/numpy/anuga/abstract_2d_finite_volumes/mesh_factory.py

    r6145 r6304  
    33"""
    44
    5 import Numeric as num
     5import numpy as num
    66
    77
    … …  
    9494    index = Index(n,m)
    9595
    96     points = num.zeros((Np, 2), num.Float)
     96    points = num.zeros((Np, 2), num.float)
    9797
    9898    for i in range(m+1):
    … …  
    106106
    107107
    108     elements = num.zeros((Nt, 3), num.Int)
     108    elements = num.zeros((Nt, 3), num.int)
    109109    boundary = {}
    110110    nt = -1

  • branches/numpy/anuga/abstract_2d_finite_volumes/neighbour_mesh.py

    r6191 r6304  
    1212from math import pi, sqrt
    1313
    14 import Numeric as num
     14import numpy as num
    1515       
    1616
    … …  
    2424    coordinate set.
    2525
    26     Coordinate sets are implemented as an N x 2 Numeric array containing
     26    Coordinate sets are implemented as an N x 2 numeric array containing
    2727    x and y coordinates.
    2828
    … …  
    3333    where
    3434
    35       coordinates is either a list of 2-tuples or an Mx2 Numeric array of
     35      coordinates is either a list of 2-tuples or an Mx2 numeric array of
    3636      floats representing all x, y coordinates in the mesh.
    3737
    38       triangles is either a list of 3-tuples or an Nx3 Numeric array of
     38      triangles is either a list of 3-tuples or an Nx3 numeric array of
    3939      integers representing indices of all vertices in the mesh.
    4040      Each vertex is identified by its index i in [0, M-1].
    … …  
    7676        """
    7777        Build triangles from x,y coordinates (sequence of 2-tuples or
    78         Mx2 Numeric array of floats) and triangles (sequence of 3-tuples
    79         or Nx3 Numeric array of non-negative integers).
     78        Mx2 numeric array of floats) and triangles (sequence of 3-tuples
     79        or Nx3 numeric array of non-negative integers).
    8080        """
    8181
    … …  
    9898
    9999        #Allocate space for geometric quantities
    100         self.centroid_coordinates = num.zeros((N, 2), num.Float)
    101 
    102         self.radii = num.zeros(N, num.Float)
    103 
    104         self.neighbours = num.zeros((N, 3), num.Int)
    105         self.neighbour_edges = num.zeros((N, 3), num.Int)
    106         self.number_of_boundaries = num.zeros(N, num.Int)
    107         self.surrogate_neighbours = num.zeros((N, 3), num.Int)
     100        self.centroid_coordinates = num.zeros((N, 2), num.float)
     101
     102        self.radii = num.zeros(N, num.float)
     103
     104        self.neighbours = num.zeros((N, 3), num.int)
     105        self.neighbour_edges = num.zeros((N, 3), num.int)
     106        self.number_of_boundaries = num.zeros(N, num.int)
     107        self.surrogate_neighbours = num.zeros((N, 3), num.int)
    108108
    109109        #Get x,y coordinates for all triangles and store
    … …  
    124124
    125125            #Compute centroid
    126             centroid = num.array([(x0 + x1 + x2)/3, (y0 + y1 + y2)/3], num.Float)
     126            centroid = num.array([(x0 + x1 + x2)/3, (y0 + y1 + y2)/3], num.float)
    127127            self.centroid_coordinates[i] = centroid
    128128
    … …  
    133133
    134134                #Midpoints
    135                 m0 = num.array([(x1 + x2)/2, (y1 + y2)/2], num.Float)
    136                 m1 = num.array([(x0 + x2)/2, (y0 + y2)/2], num.Float)
    137                 m2 = num.array([(x1 + x0)/2, (y1 + y0)/2], num.Float)
     135                m0 = num.array([(x1 + x2)/2, (y1 + y2)/2], num.float)
     136                m1 = num.array([(x0 + x2)/2, (y0 + y2)/2], num.float)
     137                m2 = num.array([(x1 + x0)/2, (y1 + y0)/2], num.float)
    138138
    139139                #The radius is the distance from the centroid of
    … …  
    394394            #Check that all keys in given boundary exist
    395395            for tag in tagged_elements.keys():
    396                 tagged_elements[tag] = num.array(tagged_elements[tag], num.Int)
     396                tagged_elements[tag] = num.array(tagged_elements[tag], num.int)
    397397
    398398                msg = 'Not all elements exist. '
    … …  
    506506            # Sanity check
    507507            if p0 is None:
    508                 raise Exception('Impossible')
    509 
     508                msg = 'Impossible: p0 is None!?'
     509                raise Exception, msg
    510510
    511511            # Register potential paths from A to B
    … …  
    613613                point_registry[tuple(p1)] = len(point_registry)
    614614           
    615             polygon.append(list(p1)) # De-Numeric each point :-)
     615            polygon.append(list(p1)) # De-numeric each point :-)
    616616            p0 = p1
    617617
    … …  
    11391139
    11401140            # Distances from line origin to the two intersections
    1141             z0 = num.array([x0 - xi0, y0 - eta0], num.Float)
    1142             z1 = num.array([x1 - xi0, y1 - eta0], num.Float)
     1141            z0 = num.array([x0 - xi0, y0 - eta0], num.float)
     1142            z1 = num.array([x1 - xi0, y1 - eta0], num.float)
    11431143            d0 = num.sqrt(num.sum(z0**2))
    11441144            d1 = num.sqrt(num.sum(z1**2))
    … …  
    11551155            # Normal direction:
    11561156            # Right hand side relative to line direction
    1157             vector = num.array([x1 - x0, y1 - y0], num.Float) # Segment vector
     1157            vector = num.array([x1 - x0, y1 - y0], num.float) # Segment vector
    11581158            length = num.sqrt(num.sum(vector**2))      # Segment length
    1159             normal = num.array([vector[1], -vector[0]], num.Float)/length
     1159            normal = num.array([vector[1], -vector[0]], num.float)/length
    11601160
    11611161
    … …  
    12371237        assert isinstance(segment, Triangle_intersection), msg
    12381238       
    1239         midpoint = num.sum(num.array(segment.segment, num.Float))/2
     1239        midpoint = num.sum(num.array(segment.segment, num.float))/2
    12401240        midpoints.append(midpoint)
    12411241

  • branches/numpy/anuga/abstract_2d_finite_volumes/pmesh2domain.py

    r6145 r6304  
    1 """Class pmesh2domain - Converting .tsh files to doamains
     1"""Class pmesh2domain - Converting .tsh files to domains
    22
    33
    … …  
    88
    99import sys
    10 
    11 import Numeric as num
    12 
    13 
    14 def pmesh_instance_to_domain_instance(mesh,
    15                                       DomainClass):
    16     """
    17     Convert a pmesh instance/object into a domain instance.
    18 
    19     Use pmesh_to_domain_instance to convert a mesh file to a domain instance.
    20     """
    21 
    22     vertex_coordinates, vertices, tag_dict, vertex_quantity_dict \
    23                         ,tagged_elements_dict, geo_reference = \
    24                         pmesh_to_domain(mesh_instance=mesh)
     10import numpy as num
     11
     12
     13##
     14# @brief Convert a pmesh instance to a domain instance.
     15# @param mesh The pmesh instance to convert.
     16# @param DomainClass The class to instantiate and return.
     17# @return The converted pmesh instance (as a 'DomainClass' instance).
     18def pmesh_instance_to_domain_instance(mesh, DomainClass):
     19    """Convert a pmesh instance/object into a domain instance.
     20
     21    Uses pmesh_to_domain_instance to convert a mesh file to a domain instance.
     22    """
     23
     24    (vertex_coordinates, vertices, tag_dict, vertex_quantity_dict,
     25     tagged_elements_dict, geo_reference) = pmesh_to_domain(mesh_instance=mesh)
     26
     27    # NOTE(Ole): This import cannot be at the module level
     28    #            due to mutual dependency with domain.py
     29    from anuga.abstract_2d_finite_volumes.domain import Domain
     30
     31    # ensure that the required 'DomainClass' actually is an instance of Domain
     32    msg = ('The class %s is not a subclass of the generic domain class %s'
     33           % (DomainClass, Domain))
     34    assert issubclass(DomainClass, Domain), msg
     35
     36    # instantiate the result class
     37    result = DomainClass(coordinates=vertex_coordinates,
     38                         vertices=vertices,
     39                         boundary=tag_dict,
     40                         tagged_elements=tagged_elements_dict,
     41                         geo_reference=geo_reference)
     42
     43    # set the water stage to be the elevation
     44    if (vertex_quantity_dict.has_key('elevation') and
     45        not vertex_quantity_dict.has_key('stage')):
     46        vertex_quantity_dict['stage'] = vertex_quantity_dict['elevation']
     47    result.set_quantity_vertices_dict(vertex_quantity_dict)
     48
     49    return result
     50
     51
     52##
     53# @brief Convert a mesh file to a Domain instance.
     54# @param file_name Name of the file to convert (TSH or MSH).
     55# @param DomainClass Class of return instance.
     56# @param use_cache True if caching is to be used.
     57# @param verbose True if this function is to be verbose.
     58# @return An instance of 'DomainClass' containing the file data.
     59def pmesh_to_domain_instance(file_name, DomainClass, use_cache=False,
     60                             verbose=False):
     61    """Converts a mesh file(.tsh or .msh), to a Domain instance.
     62
     63    file_name is the name of the mesh file to convert, including the extension
     64
     65    DomainClass is the Class that will be returned.
     66    It must be a subclass of Domain, with the same interface as domain.
     67
     68    use_cache: True means that caching is attempted for the computed domain.   
     69    """
     70
     71    if use_cache is True:
     72        from caching import cache
     73        result = cache(_pmesh_to_domain_instance, (file_name, DomainClass),
     74                       dependencies=[file_name], verbose=verbose)
     75    else:
     76        result = apply(_pmesh_to_domain_instance, (file_name, DomainClass))       
     77       
     78    return result
     79
     80
     81##
     82# @brief Convert a mesh file to a Domain instance.
     83# @param file_name Name of the file to convert (TSH or MSH).
     84# @param DomainClass Class of return instance.
     85# @return The DomainClass instance containing the file data.
     86def _pmesh_to_domain_instance(file_name, DomainClass):
     87    """Converts a mesh file(.tsh or .msh), to a Domain instance.
     88
     89    Internal function. See public interface pmesh_to_domain_instance for details
     90    """
     91   
     92    (vertex_coordinates, vertices, tag_dict, vertex_quantity_dict,
     93     tagged_elements_dict, geo_reference) = pmesh_to_domain(file_name=file_name)
    2594
    2695    # NOTE(Ole): This import cannot be at the module level due to mutual
    … …  
    2897    from anuga.abstract_2d_finite_volumes.domain import Domain
    2998
    30  
    31 
    32 
    33     msg = 'The class %s is not a subclass of the generic domain class %s'\
    34           %(DomainClass, Domain)
     99    # ensure the required class is a subclass of Domain
     100    msg = ('The class %s is not a subclass of the generic domain class %s'
     101           % (DomainClass, Domain))
    35102    assert issubclass(DomainClass, Domain), msg
    36 
    37103
    38104    domain = DomainClass(coordinates = vertex_coordinates,
    … …  
    42108                         geo_reference = geo_reference )
    43109
    44     # set the water stage to be the elevation
    45     if vertex_quantity_dict.has_key('elevation') and not vertex_quantity_dict.has_key('stage'):
    46         vertex_quantity_dict['stage'] = vertex_quantity_dict['elevation']
    47 
    48     domain.set_quantity_vertices_dict(vertex_quantity_dict)
    49     #print "vertex_quantity_dict",vertex_quantity_dict
    50     return domain
    51 
    52 
    53 
    54 def pmesh_to_domain_instance(file_name, DomainClass, use_cache = False, verbose = False):
    55     """
    56     Converts a mesh file(.tsh or .msh), to a Domain instance.
    57 
    58     file_name is the name of the mesh file to convert, including the extension
    59 
    60     DomainClass is the Class that will be returned.
    61     It must be a subclass of Domain, with the same interface as domain.
    62 
    63     use_cache: True means that caching is attempted for the computed domain.   
    64     """
    65 
    66     if use_cache is True:
    67         from caching import cache
    68         result = cache(_pmesh_to_domain_instance, (file_name, DomainClass),
    69                        dependencies = [file_name],                     
    70                        verbose = verbose)
    71 
    72     else:
    73         result = apply(_pmesh_to_domain_instance, (file_name, DomainClass))       
    74        
    75     return result
    76 
    77 
    78 
    79 
    80 def _pmesh_to_domain_instance(file_name, DomainClass):
    81     """
    82     Converts a mesh file(.tsh or .msh), to a Domain instance.
    83 
    84     Internal function. See public interface pmesh_to_domain_instance for details
    85     """
    86    
    87     vertex_coordinates, vertices, tag_dict, vertex_quantity_dict, \
    88                         tagged_elements_dict, geo_reference = \
    89                         pmesh_to_domain(file_name=file_name)
    90 
    91 
    92     # NOTE(Ole): This import cannot be at the module level due to mutual
    93     # dependency with domain.py
    94     from anuga.abstract_2d_finite_volumes.domain import Domain
    95 
    96 
    97     msg = 'The class %s is not a subclass of the generic domain class %s'\
    98           %(DomainClass, Domain)
    99     assert issubclass(DomainClass, Domain), msg
    100 
    101 
    102 
    103     domain = DomainClass(coordinates = vertex_coordinates,
    104                          vertices = vertices,
    105                          boundary = tag_dict,
    106                          tagged_elements = tagged_elements_dict,
    107                          geo_reference = geo_reference )
    108 
    109 
    110 
    111     #FIXME (Ole): Is this really the right place to apply the a default
    112     #value specific to the shallow water wave equation?
    113     #The 'assert' above indicates that any subclass of Domain is acceptable.
    114     #Suggestion - module shallow_water.py will eventually take care of this
    115     #(when I get around to it) so it should be removed from here.
     110    # FIXME (Ole): Is this really the right place to apply a default
     111    # value specific to the shallow water wave equation?
     112    # The 'assert' above indicates that any subclass of Domain is acceptable.
     113    # Suggestion - module shallow_water.py will eventually take care of this
     114    # (when I get around to it) so it should be removed from here.
    116115
    117116    # This doesn't work on the domain instance.
    118117    # This is still needed so -ve elevations don't cuase 'lakes'
    119118    # The fixme we discussed was to only create a quantity when its values
    120     #are set.
     119    # are set.
    121120    # I think that's the way to go still
    122121
    123122    # set the water stage to be the elevation
    124     if vertex_quantity_dict.has_key('elevation') and not vertex_quantity_dict.has_key('stage'):
     123    if (vertex_quantity_dict.has_key('elevation') and
     124        not vertex_quantity_dict.has_key('stage')):
    125125        vertex_quantity_dict['stage'] = vertex_quantity_dict['elevation']
    126 
    127126    domain.set_quantity_vertices_dict(vertex_quantity_dict)
    128     #print "vertex_quantity_dict",vertex_quantity_dict
     127
    129128    return domain
    130129
    131130
    132 def pmesh_to_domain(file_name=None,
    133                     mesh_instance=None,
    134                     use_cache=False,
     131##
     132# @brief Convert pmesh file/instance to list(s) that can instantiate a Domain.
     133# @param file_name Path to file to convert.
     134# @param mesh_instance Instance to convert.
     135# @param use_cache True if we are to cache.
     136# @param verbose True if this function is to be verbose.
     137# @return ??
     138def pmesh_to_domain(file_name=None, mesh_instance=None, use_cache=False,
    135139                    verbose=False):
    136     """
    137     Convert a pmesh file or a pmesh mesh instance to a bunch of lists
     140    """Convert a pmesh file or a pmesh mesh instance to a bunch of lists
    138141    that can be used to instanciate a domain object.
    139142
    … …  
    144147        from caching import cache
    145148        result = cache(_pmesh_to_domain, (file_name, mesh_instance),
    146                        dependencies = [file_name],                     
    147                        verbose = verbose)
     149                       dependencies=[file_name], verbose=verbose)
    148150
    149151    else:
    … …  
    153155
    154156
    155 def _pmesh_to_domain(file_name=None,
    156                      mesh_instance=None,
    157                      use_cache=False,
     157##
     158# @brief Convert pmesh file/instance to list(s) that can instantiate a Domain.
     159# @param file_name Path to file to convert.
     160# @param mesh_instance Instance to convert.
     161# @param use_cache True if we are to cache.
     162# @param verbose True if this function is to be verbose.
     163# @return ??
     164def _pmesh_to_domain(file_name=None, mesh_instance=None, use_cache=False,
    158165                     verbose=False):
    159     """
    160     Convert a pmesh file or a pmesh mesh instance to a bunch of lists
     166    """Convert a pmesh file or a pmesh mesh instance to a bunch of lists
    161167    that can be used to instanciate a domain object.
    162168    """
    … …  
    164170    from load_mesh.loadASCII import import_mesh_file
    165171
     172    # get data from mesh instance or file
    166173    if file_name is None:
    167174        mesh_dict = mesh_instance.Mesh2IODict()
    168175    else:
    169176        mesh_dict = import_mesh_file(file_name)
    170     #print "mesh_dict",mesh_dict
     177
     178    # extract required data from the mesh dictionary
    171179    vertex_coordinates = mesh_dict['vertices']
    172180    volumes = mesh_dict['triangles']
    173181    vertex_quantity_dict = {}
    174     point_atts = num.transpose(mesh_dict['vertex_attributes'])
     182    point_atts = mesh_dict['vertex_attributes']
    175183    point_titles  = mesh_dict['vertex_attribute_titles']
    176184    geo_reference  = mesh_dict['geo_reference']
    177     if point_atts != None:
    178         for quantity, value_vector in map (None, point_titles, point_atts):
     185    if point_atts is not None:
     186        point_atts = num.transpose(mesh_dict['vertex_attributes'])
     187        for quantity, value_vector in map(None, point_titles, point_atts):
    179188            vertex_quantity_dict[quantity] = value_vector
    180189    tag_dict = pmesh_dict_to_tag_dict(mesh_dict)
    181190    tagged_elements_dict = build_tagged_elements_dictionary(mesh_dict)
    182     return vertex_coordinates, volumes, tag_dict, vertex_quantity_dict, tagged_elements_dict, geo_reference
    183 
     191
     192    return (vertex_coordinates, volumes, tag_dict, vertex_quantity_dict,
     193            tagged_elements_dict, geo_reference)
    184194
    185195
    186196def build_tagged_elements_dictionary(mesh_dict):
    187197    """Build the dictionary of element tags.
     198
    188199    tagged_elements is a dictionary of element arrays,
    189     keyed by tag:
    190     { (tag): [e1, e2, e3..] }
    191     """
     200    keyed by tag: { (tag): [e1, e2, e3..] }
     201    """
     202
    192203    tri_atts = mesh_dict['triangle_tags']
    193204    tagged_elements = {}
    … …  
    198209            tagged_elements.setdefault(tri_atts[tri_att_index],
    199210                                       []).append(tri_att_index)
     211
    200212    return tagged_elements
     213
    201214
    202215def pmesh_dict_to_tag_dict(mesh_dict):
    … …  
    204217    to a dictionary of tags, indexed with volume id and face number.
    205218    """
     219
    206220    triangles = mesh_dict['triangles']
    207221    sides = calc_sides(triangles)
    208222    tag_dict = {}
    209     for seg, tag in map(None, mesh_dict['segments'],
    210                         mesh_dict['segment_tags']):
     223    for seg, tag in map(None, mesh_dict['segments'], mesh_dict['segment_tags']):
    211224        v1 = int(seg[0])
    212225        v2 = int(seg[1])
    … …  
    222235
    223236def calc_sides(triangles):
    224     #Build dictionary mapping from sides (2-tuple of points)
    225     #to left hand side neighbouring triangle
     237    '''Build dictionary mapping from sides (2-tuple of points)
     238    to left hand side neighbouring triangle
     239    '''
     240
    226241    sides = {}
     242
    227243    for id, triangle in enumerate(triangles):
    228244        a = int(triangle[0])
    229245        b = int(triangle[1])
    230246        c = int(triangle[2])
     247
    231248        sides[a,b] = (id, 2) #(id, face)
    232249        sides[b,c] = (id, 0) #(id, face)
    233250        sides[c,a] = (id, 1) #(id, face)
     251
    234252    return sides
     253

  • branches/numpy/anuga/abstract_2d_finite_volumes/quantity.py

    r6244 r6304  
    2323from anuga.caching import cache
    2424
    25 import Numeric as num
     25import anuga.utilities.numerical_tools as aunt
     26
     27import numpy as num
    2628
    2729
    … …  
    4345        if vertex_values is None:
    4446            N = len(domain)             # number_of_elements
    45             self.vertex_values = num.zeros((N, 3), num.Float)
     47            self.vertex_values = num.zeros((N, 3), num.float)
    4648        else:
    47             self.vertex_values = num.array(vertex_values, num.Float)
     49            self.vertex_values = num.array(vertex_values, num.float)
    4850
    4951            N, V = self.vertex_values.shape
    … …  
    5759
    5860        # Allocate space for other quantities
    59         self.centroid_values = num.zeros(N, num.Float)
    60         self.edge_values = num.zeros((N, 3), num.Float)
     61        self.centroid_values = num.zeros(N, num.float)
     62        self.edge_values = num.zeros((N, 3), num.float)
    6163
    6264        # Allocate space for Gradient
    63         self.x_gradient = num.zeros(N, num.Float)
    64         self.y_gradient = num.zeros(N, num.Float)
     65        self.x_gradient = num.zeros(N, num.float)
     66        self.y_gradient = num.zeros(N, num.float)
    6567
    6668        # Allocate space for Limiter Phi
    67         self.phi = num.zeros(N, num.Float)
     69        self.phi = num.zeros(N, num.float)
    6870
    6971        # Intialise centroid and edge_values
    … …  
    7274        # Allocate space for boundary values
    7375        L = len(domain.boundary)
    74         self.boundary_values = num.zeros(L, num.Float)
     76        self.boundary_values = num.zeros(L, num.float)
    7577
    7678        # Allocate space for updates of conserved quantities by
    … …  
    7880
    7981        # Allocate space for update fields
    80         self.explicit_update = num.zeros(N, num.Float )
    81         self.semi_implicit_update = num.zeros(N, num.Float )
    82         self.centroid_backup_values = num.zeros(N, num.Float)
     82        self.explicit_update = num.zeros(N, num.float )
     83        self.semi_implicit_update = num.zeros(N, num.float )
     84        self.centroid_backup_values = num.zeros(N, num.float)
    8385
    8486        self.set_beta(1.0)
    … …  
    316318
    317319        numeric:
    318           Compatible list, Numeric array (see below) or constant.
     320          Compatible list, numeric array (see below) or constant.
    319321          If callable it will treated as a function (see below)
    320322          If instance of another Quantity it will be treated as such.
    … …  
    357359
    358360                  In case of location == 'centroids' the dimension values must
    359                   be a list of a Numerical array of length N,
     361                  be a list of a numerical array of length N,
    360362                  N being the number of elements.
    361363                  Otherwise it must be of dimension Nx3
    … …  
    458460
    459461        msg = 'Indices must be a list or None'
    460         assert type(indices) in [ListType, NoneType, num.ArrayType], msg
     462        assert (type(indices) in [ListType, NoneType]
     463                or isinstance(indices, num.ndarray)), msg
    461464
    462465        # Determine which 'set_values_from_...' to use
    463466        if numeric is not None:
    464             if type(numeric) in [FloatType, IntType, LongType]:
    465                 self.set_values_from_constant(numeric, location,
    466                                               indices, verbose)
    467             elif type(numeric) in [num.ArrayType, ListType]:
     467            if isinstance(numeric, num.ndarray) or isinstance(numeric, list):
    468468                self.set_values_from_array(numeric, location, indices,
    469469                                           use_cache=use_cache, verbose=verbose)
    … …  
    479479                                                     indices, verbose=verbose,
    480480                                                     use_cache=use_cache)
     481#            if (isinstance(numeric, float) or isinstance(numeric, int)
     482#                or isinstance(numeric, long)):
    481483            else:
    482                 msg = 'Illegal type for argument numeric: %s' % str(numeric)
    483                 raise msg
     484                try:
     485                    numeric = float(numeric)
     486                except ValueError:
     487                    msg = ("Illegal type for variable 'numeric': "
     488                           "%s, shape=%s\n%s"
     489                           % (type(numeric), numeric.shape, str(numeric)))
     490                    msg += ('\ntype(numeric)==FloatType is %s'
     491                            % str(type(numeric)==FloatType))
     492                    msg += ('\nisinstance(numeric, float)=%s'
     493                            % str(isinstance(numeric, float)))
     494                    msg += ('\nisinstance(numeric, num.ndarray)=%s'
     495                            % str(isinstance(numeric, num.ndarray)))
     496                    raise Exception, msg
     497                self.set_values_from_constant(numeric, location,
     498                                              indices, verbose)
     499        elif quantity is not None:
     500            if type(numeric) in [FloatType, IntType, LongType]:
     501                self.set_values_from_constant(numeric, location,
     502                                              indices, verbose)
     503            else:
     504                msg = ("Illegal type for variable 'numeric': %s, shape=%s\n%s"
     505                       % (type(numeric), numeric.shape, str(numeric)))
     506                msg += ('\ntype(numeric)==FloatType is %s'
     507                        % str(type(numeric)==FloatType))
     508                msg += ('\nisinstance(numeric, float)=%s'
     509                        % str(isinstance(numeric, float)))
     510                msg += ('\nisinstance(numeric, num.ndarray)=%s'
     511                        % str(isinstance(numeric, num.ndarray)))
     512                raise Exception, msg
    484513        elif quantity is not None:
    485514            self.set_values_from_quantity(quantity, location, indices, verbose)
    … …  
    583612        """Set values for quantity
    584613
    585         values: Numeric array
     614        values: numeric array
    586615        location: Where values are to be stored.
    587616        Permissible options are: vertices, centroid, unique vertices
    … …  
    593622
    594623        In case of location == 'centroid' the dimension values must
    595         be a list of a Numerical array of length N, N being the number
     624        be a list of a numerical array of length N, N being the number
    596625        of elements.
    597626
    … …  
    607636        """
    608637
    609         values = num.array(values, num.Float)
     638        values = num.array(values, num.float)
    610639
    611640        if indices is not None:
    612             indices = num.array(indices, num.Int)
     641            indices = num.array(indices, num.int)
    613642            msg = ('Number of values must match number of indices: You '
    614643                   'specified %d values and %d indices'
    … …  
    637666                    'Values array must be 1d')
    638667
    639             self.set_vertex_values(values.flat, indices=indices,
     668            self.set_vertex_values(values.flatten(), indices=indices,
    640669                                   use_cache=use_cache, verbose=verbose)
    641670        else:
    … …  
    788817        from anuga.coordinate_transforms.geo_reference import Geo_reference
    789818
    790         points = ensure_numeric(points, num.Float)
    791         values = ensure_numeric(values, num.Float)
     819        points = ensure_numeric(points, num.float)
     820        values = ensure_numeric(values, num.float)
    792821
    793822        if location != 'vertices':
    … …  
    10801109
    10811110        # Ensure that interpolation points is either a list of
    1082         # points, Nx2 array, or geospatial and convert to Numeric array
     1111        # points, Nx2 array, or geospatial and convert to numeric array
    10831112        if isinstance(interpolation_points, Geospatial_data):
    10841113            # Ensure interpolation points are in absolute UTM coordinates
    … …  
    11191148        """Get values for quantity
    11201149
    1121         Extract values for quantity as a Numeric array.
     1150        Extract values for quantity as a numeric array.
    11221151
    11231152        Inputs:
    … …  
    11371166
    11381167        In case of location == 'centroids' the dimension of returned
    1139         values will be a list or a Numerical array of length N, N being
     1168        values will be a list or a numerical array of length N, N being
    11401169        the number of elements.
    11411170
    … …  
    11731202                            'edges', 'unique vertices']:
    11741203            msg = 'Invalid location: %s' % location
    1175             raise msg
     1204            raise Exception, msg
    11761205
    11771206        import types
    11781207
    11791208        assert type(indices) in [types.ListType, types.NoneType,
    1180                                  num.ArrayType], \
    1181                    'Indices must be a list or None'
     1209                                 num.ndarray], \
     1210                   'Indices must be a list or None'     #??#
    11821211
    11831212        if location == 'centroids':
    … …  
    12101239                    sum += self.vertex_values[triangle_id, vertex_id]
    12111240                vert_values.append(sum / len(triangles))
    1212             return num.array(vert_values, num.Float)
     1241            return num.array(vert_values, num.float)
    12131242        else:
    12141243            if (indices is None):
    … …  
    12361265
    12371266        # Check that A can be converted to array and is of appropriate dim
    1238         A = ensure_numeric(A, num.Float)
     1267        A = ensure_numeric(A, num.float)
    12391268        assert len(A.shape) == 1
    12401269
    … …  
    13361365
    13371366        if precision is None:
    1338             precision = num.Float
     1367            precision = num.float
    13391368
    13401369        if smooth is True:
    … …  
    13421371            V = self.domain.get_triangles()
    13431372            N = self.domain.number_of_full_nodes # Ignore ghost nodes if any
    1344             A = num.zeros(N, num.Float)
     1373            A = num.zeros(N, num.float)
    13451374            points = self.domain.get_nodes()
    13461375
    … …  
    13601389                    if current_node == N:
    13611390                        msg = 'Current node exceeding number of nodes (%d) ' % N
    1362                         raise msg
     1391                        raise Exception, msg
    13631392
    13641393                    k += 1
    … …  
    13811410            V = self.domain.get_disconnected_triangles()
    13821411            points = self.domain.get_vertex_coordinates()
    1383             A = self.vertex_values.flat.astype(precision)
     1412            A = self.vertex_values.flatten().astype(precision)
    13841413
    13851414        # Return

  • branches/numpy/anuga/abstract_2d_finite_volumes/quantity_ext.c

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

  • branches/numpy/anuga/abstract_2d_finite_volumes/region.py

    r6145 r6304  
    88# FIXME (DSG-DSG) add better comments
    99
    10 import Numeric as num
     10import numpy as num
    1111
    1212

  • branches/numpy/anuga/abstract_2d_finite_volumes/test_domain.py

    r6195 r6304  
    77from anuga.config import epsilon
    88
    9 import Numeric as num
     9import numpy as num
    1010
    1111
    … …  
    6565
    6666
    67         assert domain.get_conserved_quantities(0, edge=1) == 0.
     67        assert num.alltrue(domain.get_conserved_quantities(0, edge=1) == 0.)
    6868
    6969
    … …  
    345345
    346346
    347         A = num.array([[1,2,3], [5,5,-5], [0,0,9], [-6,3,3]], num.Float)
    348         B = num.array([[2,4,4], [3,2,1], [6,-3,4], [4,5,-1]], num.Float)
     347        A = num.array([[1,2,3], [5,5,-5], [0,0,9], [-6,3,3]], num.float)
     348        B = num.array([[2,4,4], [3,2,1], [6,-3,4], [4,5,-1]], num.float)
    349349       
    350350        #print A
    … …  
    613613
    614614        sem = num.array([1.,1.,1.,1.])/num.array([1, 2, 3, 4])
    615         denom = num.ones(4, num.Float)-domain.timestep*sem
     615        denom = num.ones(4, num.float)-domain.timestep*sem
    616616
    617617#        x = array([1, 2, 3, 4]) + array( [.4,.3,.2,.1] )

  • branches/numpy/anuga/abstract_2d_finite_volumes/test_ermapper.py

    r6145 r6304  
    77from os import remove
    88
    9 import Numeric as num
     9import numpy as num
    1010
    1111
    … …  
    4545
    4646        # Write test data
    47         ermapper_grids.write_ermapper_data(original_grid, filename, num.Float64)
     47        ermapper_grids.write_ermapper_data(original_grid, filename, num.float64)
    4848
    4949        # Read in the test data
    50         new_grid = ermapper_grids.read_ermapper_data(filename, num.Float64)
     50        new_grid = ermapper_grids.read_ermapper_data(filename, num.float64)
    5151
    5252        # Check that the test data that has been read in matches the original data

  • branches/numpy/anuga/abstract_2d_finite_volumes/test_general_mesh.py

    r6174 r6304  
    1010from anuga.coordinate_transforms.geo_reference import Geo_reference
    1111
    12 import Numeric as num
     12import numpy as num
    1313
    1414
    … …  
    5959       
    6060        #bac, bce, ecf, dbe, daf, dae
    61         triangles = num.array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]], num.Int)
     61        triangles = num.array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]], num.int)
    6262
    6363        domain = General_mesh(nodes, triangles,
    … …  
    277277                       geo_reference = geo)
    278278        node = domain.get_node(2)       
    279         self.assertEqual(c, node)
     279        #self.assertEqual(c, node)
     280        self.failUnless(num.alltrue(c == node))
    280281       
    281282        node = domain.get_node(2, absolute=True)     
    282         self.assertEqual(nodes_absolute[2], node)
     283        #self.assertEqual(nodes_absolute[2], node)
     284        self.failUnless(num.alltrue(nodes_absolute[2] == node))
    283285       
    284286        node = domain.get_node(2, absolute=True)     
    285         self.assertEqual(nodes_absolute[2], node)
     287        #self.assertEqual(nodes_absolute[2], node)
     288        self.failUnless(num.alltrue(nodes_absolute[2] == node))
    286289       
    287290

  • branches/numpy/anuga/abstract_2d_finite_volumes/test_generic_boundary_conditions.py

    r6195 r6304  
    88from anuga.config import epsilon
    99
    10 import Numeric as num
     10import numpy as num
    1111
    1212

  • branches/numpy/anuga/abstract_2d_finite_volumes/test_ghost.py

    r6195 r6304  
    66from anuga.abstract_2d_finite_volumes.domain import *
    77from anuga.config import epsilon
     8
     9import numpy as num
    810
    911
    … …  
    4042
    4143
    42         assert domain.get_conserved_quantities(0, edge=1) == 0.
     44        assert num.alltrue(domain.get_conserved_quantities(0, edge=1) == 0.)
    4345
    4446

  • branches/numpy/anuga/abstract_2d_finite_volumes/test_neighbour_mesh.py

    r6174 r6304  
    1818from anuga.utilities.numerical_tools import ensure_numeric
    1919
    20 import Numeric as num
     20import numpy as num
    2121
    2222
    … …  
    988988                  [  75735.4765625 ,  23762.00585938],
    989989                  [  52341.70703125,  38563.39453125]]
    990 
    991         ##points = ensure_numeric(points, Int)/1000  # Simplify for ease of interpretation       
    992990
    993991        triangles = [[19, 0,15],
    … …  
    10921090                  [  35406.3359375 ,  79332.9140625 ]]
    10931091
    1094         scaled_points = ensure_numeric(points, num.Int)/1000  # Simplify for ease of interpretation
     1092        scaled_points = ensure_numeric(points, num.int)/1000  # Simplify for ease of interpretation
    10951093
    10961094        triangles = [[ 0, 1, 2],

  • branches/numpy/anuga/abstract_2d_finite_volumes/test_pmesh2domain.py

    r6191 r6304  
    1414from anuga.pmesh.mesh import importMeshFromFile
    1515
    16 import Numeric as num
     16import numpy as num
    1717
    1818

  • branches/numpy/anuga/abstract_2d_finite_volumes/test_quantity.py

    r6195 r6304  
    1515from anuga.utilities.polygon import *
    1616
    17 import Numeric as num
     17import numpy as num
    1818
    1919
    … …  
    17651765        quantity = Quantity(self.mesh4)
    17661766
    1767         quantity.vertex_values = num.array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]], num.Float)
     1767        quantity.vertex_values = num.array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]], num.float)
    17681768
    17691769        quantity.interpolate_from_vertices_to_edges()
    … …  
    17811781                                          [3., 2.5, 1.5],
    17821782                                          [3.5, 4.5, 3.],
    1783                                           [2.5, 3.5, 2]], num.Float)
     1783                                          [2.5, 3.5, 2]], num.float)
    17841784
    17851785        quantity.interpolate_from_edges_to_vertices()
    … …  
    18351835
    18361836        sem = num.array([1.,1.,1.,1.])/num.array([1, 2, 3, 4])
    1837         denom = num.ones(4, num.Float)-timestep*sem
     1837        denom = num.ones(4, num.float)-timestep*sem
    18381838
    18391839        x = num.array([1, 2, 3, 4])/denom
    … …  
    18591859
    18601860        sem = num.array([1.,1.,1.,1.])/num.array([1, 2, 3, 4])
    1861         denom = num.ones(4, num.Float)-timestep*sem
     1861        denom = num.ones(4, num.float)-timestep*sem
    18621862
    18631863        x = num.array([1., 2., 3., 4.])
    … …  
    19011901
    19021902        bed = domain.quantities['elevation'].vertex_values
    1903         stage = num.zeros(bed.shape, num.Float)
     1903        stage = num.zeros(bed.shape, num.float)
    19041904
    19051905        h = 0.03
    … …  
    19851985
    19861986        bed = domain.quantities['elevation'].vertex_values
    1987         stage = num.zeros(bed.shape, num.Float)
     1987        stage = num.zeros(bed.shape, num.float)
    19881988
    19891989        h = 0.03
    … …  
    19991999        stage = domain.quantities['stage']
    20002000        A, V = stage.get_vertex_values(xy=False, smooth=False)
    2001         Q = stage.vertex_values.flat
     2001        Q = stage.vertex_values.flatten()
    20022002
    20032003        for k in range(8):
    … …  
    25112511    suite = unittest.makeSuite(Test_Quantity, 'test')   
    25122512    #suite = unittest.makeSuite(Test_Quantity, 'test_set_values_from_file_using_polygon')
    2513 
    2514     #suite = unittest.makeSuite(Test_Quantity, 'test_set_vertex_values_using_general_interface_subset_and_geo')
    2515     #print "restricted test"
    2516     #suite = unittest.makeSuite(Test_Quantity,'verbose_test_set_values_from_UTM_pts')
    25172513    runner = unittest.TextTestRunner()
    25182514    runner.run(suite)

  • branches/numpy/anuga/abstract_2d_finite_volumes/test_region.py

    r6145 r6304  
    88#from anuga.config import epsilon
    99
    10 import Numeric as num
     10import numpy as num
    1111
    1212
    … …  
    4242
    4343    def test_region_tags(self):
    44         """
    45         get values based on triangle lists.
    46         """
    47         from mesh_factory import rectangular
    48         from shallow_water import Domain
    49 
    50         #Create basic mesh
    51         points, vertices, boundary = rectangular(1, 3)
    52 
    53         #Create shallow water domain
    54         domain = Domain(points, vertices, boundary)
    55         domain.build_tagged_elements_dictionary({'bottom':[0,1],
    56                                                  'top':[4,5],
    57                                                  'all':[0,1,2,3,4,5]})
    58 
     44        """get values based on triangle lists."""
     45
     46        from mesh_factory import rectangular
     47        from shallow_water import Domain
     48
     49        #Create basic mesh
     50        points, vertices, boundary = rectangular(1, 3)
     51
     52        #Create shallow water domain
     53        domain = Domain(points, vertices, boundary)
     54        domain.build_tagged_elements_dictionary({'bottom': [0,1],
     55                                                 'top': [4,5],
     56                                                 'all': [0,1,2,3,4,5]})
    5957
    6058        #Set friction
    … …  
    6563        b = Set_region('top', 'friction', 1.0)
    6664        domain.set_region([a, b])
    67         #print domain.quantities['friction'].get_values()
    68         assert num.allclose(domain.quantities['friction'].get_values(),\
    69                             [[ 0.09,  0.09,  0.09],
    70                              [ 0.09,  0.09,  0.09],
    71                              [ 0.07,  0.07,  0.07],
    72                              [ 0.07,  0.07,  0.07],
    73                              [ 1.0,  1.0,  1.0],
    74                              [ 1.0,  1.0,  1.0]])
     65
     66        expected = [[ 0.09,  0.09,  0.09],
     67                    [ 0.09,  0.09,  0.09],
     68                    [ 0.07,  0.07,  0.07],
     69                    [ 0.07,  0.07,  0.07],
     70                    [ 1.0,  1.0,  1.0],
     71                    [ 1.0,  1.0,  1.0]]
     72        msg = ("\ndomain.quantities['friction']=%s\nexpected value=%s"
     73               % (str(domain.quantities['friction'].get_values()),
     74                  str(expected)))
     75        assert num.allclose(domain.quantities['friction'].get_values(),
     76                            expected), msg
    7577
    7678        #c = Add_Value_To_region('all', 'friction', 10.0)
    … …  
    126128       
    127129    def test_unique_vertices(self):
    128         """
    129         get values based on triangle lists.
    130         """
     130        """get values based on triangle lists."""
     131
    131132        from mesh_factory import rectangular
    132133        from shallow_water import Domain
    … …  
    147148        a = Set_region('bottom', 'friction', 0.09, location = 'unique vertices')
    148149        domain.set_region(a)
    149         #print domain.quantities['friction'].get_values()
    150         assert num.allclose(domain.quantities['friction'].get_values(),\
     150        assert num.allclose(domain.quantities['friction'].get_values(),
    151151                            [[ 0.09,  0.09,  0.09],
    152152                             [ 0.09,  0.09,  0.09],
    … …  
    217217        #print domain.quantities['friction'].get_values()
    218218        frict_points = domain.quantities['friction'].get_values()
    219         assert num.allclose(frict_points[0],\
    220                             [ calc_frict, calc_frict, calc_frict])
    221         assert num.allclose(frict_points[1],\
    222                             [ calc_frict, calc_frict, calc_frict])
     219        expected = [calc_frict, calc_frict, calc_frict]
     220        msg = ('frict_points[0]=%s\nexpected=%s' % (str(frict_points[0]),
     221                                                    str(expected)))
     222        assert num.allclose(frict_points[0], expected), msg
     223        msg = ('frict_points[1]=%s\nexpected=%s' % (str(frict_points[1]),
     224                                                    str(expected)))
     225        assert num.allclose(frict_points[1], expected), msg
    223226 
    224227    def test_unique_vertices_average_loc_unique_vert(self):
    … …  
    262265#-------------------------------------------------------------
    263266if __name__ == "__main__":
    264     suite = unittest.makeSuite(Test_Region,'test')
     267    suite = unittest.makeSuite(Test_Region, 'test')   
    265268    runner = unittest.TextTestRunner()
    266269    runner.run(suite)

  • branches/numpy/anuga/abstract_2d_finite_volumes/test_util.py

    r6175 r6304  
    2222import string
    2323
    24 import Numeric as num
     24import numpy as num
    2525
    2626
    … …  
    630630                    q1 = F(t+60, point_id=id)
    631631
    632                 if q0 == NAN:
     632                if num.alltrue(q0 == NAN):
    633633                    actual = q0
    634634                else:
    … …  
    641641                #print "actual", actual
    642642                #print
    643                 if q0 == NAN:
    644                      self.failUnless( q == actual, 'Fail!')
     643                if num.alltrue(q0 == NAN):
     644                     self.failUnless(num.alltrue(q == actual), 'Fail!')
    645645                else:
    646646                    assert num.allclose(q, actual)
    … …  
    11811181        #FIXME: Division is not expected to work for integers.
    11821182        #This must be caught.
    1183         foo = num.array([[1,2,3], [4,5,6]], num.Float)
    1184 
    1185         bar = num.array([[-1,0,5], [6,1,1]], num.Float)                 
     1183        foo = num.array([[1,2,3], [4,5,6]], num.float)
     1184
     1185        bar = num.array([[-1,0,5], [6,1,1]], num.float)                 
    11861186
    11871187        D = {'X': foo, 'Y': bar}
    … …  
    12011201
    12021202        # make an error for zero on zero
    1203         # this is really an error in Numeric, SciPy core can handle it
     1203        # this is really an error in numeric, SciPy core can handle it
    12041204        # Z = apply_expression_to_dictionary('0/Y', D)
    12051205

  • branches/numpy/anuga/abstract_2d_finite_volumes/util.py

    r6194 r6304  
    2727from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
    2828
    29 import Numeric as num
     29import numpy as num
    3030
    3131
    … …  
    554554                expression e.g. by overloading
    555555
    556     Due to a limitation with Numeric, this can not evaluate 0/0
     556    Due to a limitation with numeric, this can not evaluate 0/0
    557557    In general, the user can fix by adding 1e-30 to the numerator.
    558558    SciPy core can handle this situation.
    … …  
    789789                        - easting, northing, name , elevation?
    790790                    - OR (this is not yet done)
    791                         - structure which can be converted to a Numeric array,
     791                        - structure which can be converted to a numeric array,
    792792                          such as a geospatial data object
    793793                     
    … …  
    12551255    n0 = int(n0)
    12561256    m = len(locations)
    1257     model_time = num.zeros((n0, m, p), num.Float)
    1258     stages = num.zeros((n0, m, p), num.Float)
    1259     elevations = num.zeros((n0, m, p), num.Float)
    1260     momenta = num.zeros((n0, m, p), num.Float)
    1261     xmom = num.zeros((n0, m, p), num.Float)
    1262     ymom = num.zeros((n0, m, p), num.Float)
    1263     speed = num.zeros((n0, m, p), num.Float)
    1264     bearings = num.zeros((n0, m, p), num.Float)
    1265     due_east = 90.0*num.ones((n0, 1), num.Float)
    1266     due_west = 270.0*num.ones((n0, 1), num.Float)
    1267     depths = num.zeros((n0, m, p), num.Float)
    1268     eastings = num.zeros((n0, m, p), num.Float)
     1257    model_time = num.zeros((n0, m, p), num.float)
     1258    stages = num.zeros((n0, m, p), num.float)
     1259    elevations = num.zeros((n0, m, p), num.float)
     1260    momenta = num.zeros((n0, m, p), num.float)
     1261    xmom = num.zeros((n0, m, p), num.float)
     1262    ymom = num.zeros((n0, m, p), num.float)
     1263    speed = num.zeros((n0, m, p), num.float)
     1264    bearings = num.zeros((n0, m, p), num.float)
     1265    due_east = 90.0*num.ones((n0, 1), num.float)
     1266    due_west = 270.0*num.ones((n0, 1), num.float)
     1267    depths = num.zeros((n0, m, p), num.float)
     1268    eastings = num.zeros((n0, m, p), num.float)
    12691269    min_stages = []
    12701270    max_stages = []
    … …  
    12781278    min_speeds = []   
    12791279    max_depths = []
    1280     model_time_plot3d = num.zeros((n0, m), num.Float)
    1281     stages_plot3d = num.zeros((n0, m), num.Float)
    1282     eastings_plot3d = num.zeros((n0, m),num.Float)
     1280    model_time_plot3d = num.zeros((n0, m), num.float)
     1281    stages_plot3d = num.zeros((n0, m), num.float)
     1282    eastings_plot3d = num.zeros((n0, m),num.float)
    12831283    if time_unit is 'mins': scale = 60.0
    12841284    if time_unit is 'hours': scale = 3600.0
    … …  
    18321832    """
    18331833       
    1834     xc = num.zeros(triangles.shape[0], num.Float) # Space for centroid info
     1834    xc = num.zeros(triangles.shape[0], num.float) # Space for centroid info
    18351835   
    18361836    for k in range(triangles.shape[0]):
    … …  
    21192119                #add tide to stage if provided
    21202120                if quantity == 'stage':
    2121                      quantity_value[quantity] = num.array(quantity_value[quantity], num.Float) \
     2121                     quantity_value[quantity] = num.array(quantity_value[quantity], num.float) \
    21222122                                                          + directory_add_tide
    21232123
    … …  
    24692469       
    24702470    #convert to array for file_function
    2471     points_array = num.array(points,num.Float)
     2471    points_array = num.array(points,num.float)
    24722472       
    24732473    points_array = ensure_absolute(points_array)
Note: See TracChangeset for help on using the changeset viewer.